adding exotic atomsto molecules in moltemplate (not defined in the force field)

I want to use Moltemplate to generate SDS data file using OPLS force field, and I have re-defined the parameters of the head group which is not included in OPLS FF. Can you give me some suggestion about how to add these parameters?

Hi Jun

Yes it is possible.

  1. Augmenting a force-field file like “oplsaa.lt” to include new atom types is possible. But if the new atom type(s) only appear in once place, then this might be more work than necessary. Editing force-field files is never fun. For completeness, I will explain how to do this in some detail below. But in your case, I think doing this is probably overkill.

  2. Alternatively, you can mix explicit bonded interactions with force-field-defined bonds in the same molecule. This is probably the cleanest and easiest way to solve this kind of problem (in moltemplate).

In other words, you could derive a new molecule type (eg. “SDS”) containing all of the atoms in your system. You have to explicitly define bonds, angles, dihedrals, and impropers (if applicable) for only the new portion of your molecule. However the ordinary (OPLSAA) force-field rules will automatically generate angles, dihedrals, and impropers for the rest of your molecule.

Here is a sketch of what such a molecule definition might look like:

-------- “sds.lt” file -------

SDS inherits OPLSAA {

write_once(“Data Masses”) {

Define the masses for the new atom type(s) here:

@atom:Na 22.9898 # Na description “sodium ion bound to O-”
}

write(“Data Atoms”) {

List all of the 43 atoms in the entire SDS molecule here.

Use the standard 7-column format (for atom_style “full”)

$atom:c12 $mol @atom:C4 0.15 -11.674 0.8164904 0.0
$atom:c11 $mol @atom:C4 0.15 -10.451 -0.8164904 0.0
:
}

write(“Data Bond List”) {
#List only the bonds which don’t include the new atom type(s) here.

#Omit the bond type (@bond) for each entry.

#This section of the file has 3-columns.
#See details below, as well as this example.
}

write(“Data Bonds”) {

List only the bonds containing the new atom type(s) here

Include the bond type (@bond) for each entry.

This section of the file has 4-columns.

In your case, I suspect it only has one line:

$bond:O_Na @bond:Om_Na $atom:Om $atom:Na

}

write(“Data Angles”) {

List only the angles containing the new atom type(s) here

$angle:S_Om_Na @angle:S_Om_Na $atom:S $atom:Om $atom:Na
}

write(“Data Dihedrals”) {

List only the dihedrals containing the new atom type(s) here

$dihedral:O1_S_Om_Na @dihedral:Oc_S_Om_Na $atom:O1 $atom:S $atom:Om $atom:Na
$dihedral:O2_S_Om_Na @dihedral:O_S_Om_Na $atom:O2 $atom:S $atom:Om $atom:Na
$dihedral:O3_S_Om_Na @dihedral:O_S_Om_Na $atom:O3 $atom:S $atom:Om $atom:Na
}

Now define the new force-field parameters referred to above:

write_once(“In Settings”) {
pair_coeff @atom:Na @atom:Na # new pair parameters go here
bond_coeff @bond:Om_Na # new bond parameters go here
angle_coeff @angle:S_Om_Na # new angle parameters go here
dihedral_coeff @dihedral:Oc_S_Om_Na # new parameters go here
dihedral_coeff @dihedral:O_S_Om_Na # new parameters go here
}

} # SDS molecule

Note: I was too lazy to lookup appropriate names for the atom types in your molecule (using the naming conventions in the “oplsaa.lt” file). Instead I refer to them here with names like “@atom:C4”, “@atom:O”, “@atom:S”, “@atom:Om”, etc… Currently the “oplsaa.lt” file uses numeric names for the atom types (unfortunately), and you will have to lookup the correct (numeric) name for each @…8471… type and substitute it into the text above.

Note that the “Data Bond List” section includes the list the bonds between atoms which are defined by the force-field you are using (eg OPLSAA). It uses exactly the same 3-column format that you would normally use when defining a new molecule that uses OPLSAA (or any other type of force-field for that matter). There is an example of that format here (the context explaining where this file is used is located here):

In contrast, the “Data Bonds” section includes only the new bonds which you want to explictly add to the SDS molecule connecting to the new atom type you created (eg “@atom:Na”). Unlike the “Data Bond List” section, the “Data Bonds” section contains 4 columns and explicitly refers to the type of bond
ou would list any additional bonds, using the “Data Bonds” section.

Warning: I wrote this example in a hurry without testing it, and there may be obvious syntax errors and other problems. But hopefully this gives you an idea how to proceed. If you spot a mistake, or if this doesn’t work, please post a followup message and I will reply.

Using explicit bonds, angles, dihedrals, and impropers is usually only done when you have a small molecule (such as water) that is simple, so that the list of these interactions is short. In your case, you have a much more complicated molecule, however you only need to define explicit angles, dihedrals, and impropers for the small portion of the molecule containing the new exotic atom type (eg “@atom:Na”).

Note that you don’t need to worry about defining short and long names for the new atom type you are creating. (eg. “@atom:Na”. If you don’t know what this means, don’t worry about it.)

------ Details for Method 1) Modifying the “oplsaa.lt” file: ---------

You will need to become familiar with the format of the “oplsaa.lt” file.
You will need to define new atom types by adding new lines to the following sections of the “oplsaa.lt” file:

“Data Masses”
“In Charges”
“replace{}”
“pair_coeff”

If the new atom types are covalently bonded to other atoms, then unfortunately you will have to define new types of bonds, angles, dihedrals, and (perhaps) impropers which contain the new atom type(s). Generally speaking, you would do this by adding lines to the “Bonds By Type”, “Angles By Type”, “Dihedrals By Type”, “Impropers By Type”, as well as defining new “bond_coeff”, “angle_coeff”, “dihedral_coeff”, and “improper_coeff” commands.

Unfortunately, the original “oplsaa.lt” file is too complicated to understand.

Fortunately, there is a simplified example of the “oplsaa.lt” file located as part of the simple “alkane-chain” example:
http://moltemplate.org/examples/force_fields/index.html#alkane_chain_single

The force-field file (“oplsaa_simple.lt”) for this example is located here:
http://moltemplate.org/examples/alkane_chain_single/oplsaa_simple.lt

Look at the “Bonds By Type” section (line 117)

write_once(“Data Bonds By Type”) {
@bond:13-13 @atom:_b13_a_d*_i* @atom:_b13_a_d*_i*
@bond:13-46 @atom:_b13_a_d*_i* @atom:_b46_a_d*_i*
@bond:13-47 @atom:_b13_a_d*_i* @atom:_b47_a_d*_i*
@bond:46-47 @atom:_b46_a_d*_i* @atom:_b47_a_d*_i*
@bond:47-47 @atom:_b47_a_d*_i* @atom:_b47_a_d*_i*
} #(end of bonds by type)

The first line of this section means:

…Whenever an atom whose type-name contains the string “b13” is bonded to another atom containing “b13”, then assign the type of that bond to be “@bond:13-13”. Elsewhere in in the file (on line 105), we specify that the following force-field parameters should be used for a bond of this type:

bond_coeff @bond:13-13 harmonic 268.0 1.529

The same procedure is used to define “Angles By Type”, “Dihedrals By Type”, and “Impropers By Type”.
(Note: The original “oplsaa.lt” file no longer includes the force-field style names, such as “harmonic” in each “_coeff” command, so don’t include them in your modified “oplsaa.lt” file. I suppose I should update the “oplsaa_simple.lt” file for consistency.)

Notice that there are 3 different types of atoms containing the string “b13

@atom:CT3_b13_a13_d13_i13
@atom:CT2_b13_a13_d13_i13
@atom:CTH_b13_a13_d13_i13

It means that these 3 types of atoms all use the same kind of bonds.
(This makes the force-field file much shorter.)

Unfortunately “@atom:CT2_b13_a13_d13_i13” is long and difficult to type.
For this reason, each of these atom types has a SHORT name:

@atom:CT3
@atom:CT2
@atom:CTH

Incidentally, when a user wants to refer to one of these atom types, they can type the SHORT name instead of the long name. For example, here is the definition of the “CH2” monomer:

http://moltemplate.org/examples/force_fields/alkane_chain_single/ch2group.lt

Notice that the first atom is called “@atom:CT2” instead of “@atom:CT2_b13_a13_d13_i13”.
The same atom type has two different names. That is the purpose of the “replace” command. The “replace” command creates a short-name for every atom type:

replace{ @atom:CT3 @atom:CT3_b13_a13_d13_i13 }
replace{ @atom:CT2 @atom:CT2_b13_a13_d13_i13 }
replace{ @atom:CTH @atom:CTH_b13_a13_d13_i13}

---- defining independent force-field modules ----

You can either change the original “oplsaa.lt” file, OR you can create a new file which augments the original “oplsaa.lt” file without modifying it. This way, you can distribute your modifications to the OPLSAA force field as a separate module which can be loaded optionally. (This is how we implemented the LOPLSAA force field. The “oplsaa.lt” file contains only the original OPLSAA atom types and force field parameters. But the “loplsaa.lt” file defines new atom types, and contains the definitions of bonded interactions between these new atom types and all of the original atom types in the original OPLSAA force field:
https://github.com/jewettaij/moltemplate/blob/master/moltemplate/force_fields/loplsaa.lt

Doing it this way is completely optional.
Don’t worry about it now.
If you get this far, I can help you with this step.

----- It does not have to be this complicated -----

In your case, you don’t have to worry about defining bonds between your new atoms and every possible atom in the OPLSAA force field. You only need to define rules for how to bond your new atom type(s) with one or two atom types that you know it will be bonded to within the SDS molecule.

Whew!
(long email)

Cheers

Andrew

Would you mind if I distribute your SDS molecule with the moltemplate examples once you get it working?