Examples

Assign atom types using customized typing rule

Atom type is the most fundamental element in a force field. It determines which parameters should be used for describing the potential energy of each topological element (bond, angle, etc…) in a simulation system. The assignment of atom types is usually the most cumbersome task in a simulation workflow. For bio simulations, template matching is commonly used, because almost all bio-polymers are made of several kinds of repeating units. However, the template matching won’t work well for other molecules because there is no common repeating units.

mstk provide a typing engine SmartsTyper for assigning atom types. It is based on local chemical environment defined by SMARTS pattern and a hierarchical rule described in a type definition file.

In order to write the type definition file correctly, knowledge about SMARTS pattern is required. The idea of the hierarchical rule is described in this article.

Herein, an example of using SmartsTyper to assign atom types is provided. A type definition file should be fed to this typing engine. The text below shows the type definition for linear alkane in TEAM force field. Let’s copy its content and save it to file alkane.smt.

TypeDefinition

h_1    [H][CX4]
c_4    [CX4]
c_4h2  [CX4;H2]
c_4h3  [CX4;H3]

HierarchicalTree

h_1
c_4
    c_4h2
    c_4h3

The script below will assign atom types for butane using this type definition.

from mstk.forcefield.typer import SmartsTyper
from mstk.topology import Molecule

# Construct a typing engine from type definition file
typer = SmartsTyper('alkane.smt')

# Create a molecule from SMILES string
butane = Molecule.from_smiles('CCCC')

# Assign atom types based on SMARTS pattern
typer.type(butane)
for atom in butane.atoms:
    print(atom.name, atom.type)

A predefined typer typer_primitive is shipped with mstk for demonstration purpose.

Generate input files for simulation engines

Note: this example requires that Packmol is installed.

In order to generate input files for simulation engines, a force filed file is required. Currently, four formats are supported for storing force field:

1. .zff - format used by mstk. Refer to Zff for details.

2. .zfp - deprecated XML format used by mstk. Refer to Zfp for details.

3. .ppf - format used by DFF program. Refer to Ppf for details.

4. .ff - format used by Agilio Padua’s fftool. Refer to Padua for details.

Herein, an example force field file for linear alkanes is given in ZFF format. The parameters are taken from TEAM_MGI force field. Let’s copy its content and save it to file alkane.zff.

Setting          vdw_cutoff       1.0
Setting          vdw_long_range   correct
Setting          lj_mixing_rule   lorentz-berthelot
Setting          scale_14_vdw     0.5
Setting          scale_14_coulomb 0.8333333333333334

#AtomType        name           mass    charge eqt_bci   eqt_vdw   eqt_bond  eqt_ang_c eqt_ang_s eqt_dih_c eqt_dih_s eqt_imp_c eqt_imp_s eqt_polar
AtomType         c_4h2       12.0110    0.0000 c_4       c_4h2     c_4       c_4       c_4       c_4       c_4       c_4       c_4       c_4h2
AtomType         c_4h3       12.0110    0.0000 c_4       c_4h3     c_4       c_4       c_4       c_4       c_4       c_4       c_4       c_4h3
AtomType         h_1          1.0079    0.0000 h_1       h_1       h_1       h_1       h_1       h_1       h_1       h_1       h_1       h_1
#ChargeIncrement type1     type2         value
ChargeIncrement  c_4       h_1         -0.0600
#LJ126           type1     type2       epsilon     sigma
LJ126            c_4h2     c_4h2        0.3138    0.3492
LJ126            c_4h3     c_4h3        0.3284    0.3645
LJ126            h_1       h_1          0.0891    0.2434
#HarmonicBond    type1     type2        length         k  fixed
HarmonicBond     c_4       c_4          0.1533   89309.6  False
HarmonicBond     c_4       h_1          0.1097  150741.0   True
#HarmonicAngle   type1     type2     type3         theta         k  fixed
HarmonicAngle    c_4       c_4       c_4        114.1800  260.5306  False
HarmonicAngle    c_4       c_4       h_1        110.1300  190.5402  False
HarmonicAngle    h_1       c_4       h_1        106.6800  149.1466  False
#OplsDihedral    type1     type2     type3     type4            k1        k2        k3        k4
OplsDihedral     c_4       c_4       c_4       c_4         -1.6891   -0.0075    0.2657    0.0000
OplsDihedral     c_4       c_4       c_4       h_1         -1.0544    0.7448    0.7870    0.0000
OplsDihedral     h_1       c_4       c_4       h_1          0.6485    1.0678    0.6226    0.0000

The workflow is as follows:

from mstk.topology import Molecule, Topology
from mstk.forcefield import ForceField, SmartsTyper
from mstk.simsys import System
from mstk.wrapper import Packmol

typer = SmartsTyper('alkane.smt')

butane = Molecule.from_smiles('CCCC    butane')
typer.type(butane)

# Load force field parameters from ZFF file
ff = ForceField.open('alkane.zff')

# Initialize a topology with periodic boundary condition
# For now, it contains only one butane molecule
top = Topology([butane])
top.cell.set_box([3, 3, 3])

# Assign atomic charges based on the charge parameters in the force field
ff.assign_charge(top)

# Call Packmol to build a configuration containing 100 butane molecules
packmol = Packmol(r'/path/of/packmol')
top.scale_with_packmol([100], packmol=packmol)

# Associate force field parameters to the topology
# And then export input files for simulation engines
system = System(top, ff)
system.export_gromacs()
system.export_lammps()
system.export_namd()