Simulating molecules

With RUMD you can simulate molecular systems. The intra-molecular force field includes bond, angle and torsion interactions. The total potential energy due to intra-molecular interactions excluding possible pair potentials is given by

$\displaystyle U(\vec{r}) = \frac{1}{2}\sum_{\mbox{\scriptsize {bonds}}}k_s^{(i)... ...)})]^2 + \sum_{\mbox{\scriptsize {dihed}}} \sum_{n=0}^5 c_n^{(i)} \cos^n(\phi),$

(1)

where $k_s^{(i)}$ is the spring force constant for bond type $i$

, $k_\theta^{(i)}$ the angle force constant for angle force type $i$

, and $c_n^{(i)}$ the torsion coefficients for torsional force type $i$

. $l_b^{(i)}$ and $\theta_0^{(i)}$ are the zero force bond length and angle, respectively.

Beside the standard harmonic bond potential RUMD also supports simulation of rigid bonds using constraint method as well as the Finite Extensible Nonlinear Elastic (FENE) potential

$\displaystyle U(\vec{r}) = -\frac{1}{2}kR_0^2\sum_{\mbox{\scriptsize {bonds}}} \ln\left[ 1 -\left(\frac{r_{ij}}{R_0}\right)^2\right],$

(2)

where

and

(in reduced units). At the moment the constraint method is applicable for molecules with few constraints.

Example
In all one starts by creating (or copying and modifying) a topology file (with extension .top). Examples can be found in the subdirectory Conf/Mol, in particular one called mol.top, which is associated with a configuration file ExampleMol.xyz.gz. Copy both of these files to your test directory. They specify a system containing 100 polymeric molecules, each consisting of 10 monomer units. The appropriate lines to include in the python script include one for reading the topology file and one for setting the parameters of the (in this case) single bond-type. Note that this has changed for Version 3.5, so that adding intra-molecular interactions, such as bonds, works similarly to adding ordinary pair interactions.

sim = Simulation("ExampleMol.xyz.gz") # read topology file sim.ReadMoleculeData("mol.top")

# create integrator object itg = IntegratorNVE(timeStep=0.0025) sim.SetIntegrator(itg)

# create pair potential object potential = Pot_LJ_12_6(cutoff_method=ShiftedPotential) potential.SetParams(i=0, j=0, Epsilon=1.0, Sigma=1.0, Rcut=2.5) sim.AddPotential(potential)

# define harmonic bond and its parameters for bonds of type 0 pot_harm = rumd.BondHarmonic() pot_harm.SetParams(bond_type=0, bond_length=0.75, stiffness=200.0, exclude=True) sim.AddPotential(pot_harm)

sim.Run(20000)Specifying exclude=True ensure that contributions from the bonds will be removed from the calculation of the pair potential. If you wish to keep the pair force interactions between the bonded particles you can specify this using

pot_harm.SetParams(bond_type=0, bond_length=0.75, stiffness=200.0, exclude=False)

$\Box$

If there are other bond-types which have different lengths or stiffnesses, you just call SetParams again; it is not necessary to create another BondHarmonic object (this will in fact give an error). In the case you wish to use the FENE potential you simply use

pot_fene = rumd.BondFENE()
pot_fene.SetParams(bond_type=0, bond_length=0.75, stiffness=30., exclude=False)
sim.AddPotential(pot_fene)

to specify that bond type 0 is a FENE bond type. In the FENE potential, the pair interactions between bonded particles should not be excluded.

As noted above, you can also simulate molecules with rigid bonds. To specify that bond type 0 is a rigid bond you add the bond constraint using an object of type ConstraintPotential:

pot_cons = rumd.ConstraintPotential()
pot_cons.SetParams(bond_type=0, bond_length=0.75)
sim.AddPotential(pot_cons)

Details of the format used for the .top files and on tools available for creating them can be found in the user manual.