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 $k=30$ and $R_0=1.5$ (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.

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.SetPotential(potential)

# define harmonic bond and its parameters for bonds of type 0 sim.SetBondHarmonic(bond_type=0, lbond=0.75, ks=200.0)

sim.Run(20000)Note that when you specify the bond parameters after calling SetPotential contributions from those bonds will automatically be removed from the calculation of the pair potential, otherwise this will not be the case. For this reason it is important that you call SetPotential before calling sim.SetBondHarmonic. If you wish to keep the pair force interactions between the bonded particles you can specify this using

sim.SetBondHarmonic(bond_type=0, lbond=0.75, ks=200.0, exclude=False)
$\Box$

In the case you wish to use the FENE potential you simply use

sim.SetBondFENE(bond_type=0)
to specify that bond type 0 is a FENE bond type. In the FENE potential, the pair interactions between bonded particles are not 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 the SetBondConstraint method

sim.SetBondConstraint(bond_type=0, lbond=0.75)

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

Heine Larsen 2017-07-21