OPLS-AA adding new rules

EMC
1

Hi all-

I am interested in using OPLS-AA to run a comparison MD simulation with PCFF for the PIM-1 polymer of intrinsic microporosity. I have been carefully studying the atom types in the most recently updated OPLS/2020 .par and .sb force field files, which I have attached here for reference, and I think I’ve correctly identified all the atom types that I need to build the molecule.

PIM-1 is a ladder polymer, which cannot easily be built as a polymer in EMC due to its connectivity; my workaround for this is to create a very long SMILES string, such that it represents an entire 16mer or the like. For now I am interested in building a dimer to test that my system is working.

I have updated a .define file copied over from OPLS-AA bundled with EMC, and added the necessary atom types and rules for typing the atoms in PIM-1. When I execute setup.sh in my directory, a pair of .prm and .top files are generated which include the atom types and, crucially, have imported the correct charges and Lennard-Jones parameters from the source. However, the only bond, angle, dihedral, or impropers that are included are coming from either Torsion_Auto or are for the benzene ring bonds and angles between C and H. Additionally, if I try to run a setup.esh file using these parameters, I get an error, core/field.c:521 FieldIndexInit: No rule for type ‘?’. I haven’t set the dummy atom in the Masses section because that wasn’t done in the original .define, but perhaps it needs to be there?

I think my confusion lies in that I’m not sure how to appropriately use the Extra, Bond, Angle, Torsion, and Improper sections in the .define file to ensure that EMC identifies the right pairs/triplets/quadruplets in the source. Once I figure this part out I think it will be fairly straightforward for me to continue expanding this OPLS/2020 implementation for EMC (and as I do it, I’d be happy to continue sharing it here).

Would anyone (@veld , perhaps) with more expertise on this know what I might be missing?

Thank you all so much,

Sam

opls-aa.define (5.2 KB)
opls-aa.prm (4.6 KB)
opls-aa.top (1.5 KB)
oplsaa.par (238.3 KB)
oplsaa.sb (163.9 KB)
polymer.esh (1.6 KB)

2

Hi Sam,

I noticed, that you are trying to use the ITEM REPAIR section to rename existing types. Unfortunately, this section was only intended for repairing types of ions, as you can see in the original opls-aa.define file in $EMC_ROOT/field/opls/2012. Ions only have nonbonded entries. The section was not intended to rename types and therefore does not alter the types in the bonded entries. This then results in your renamed types not being reflected by the bonded entries in the resulting opls-aa.prm.

I think what you are trying to do is covered by the ITEM REDEFINE section. The original opls-aa.define redefines original OPLS types with * in it to avoid conflicts in EMC, since an * is interpreted as a wildcard (e.g. N* is redefined to NG). This, however, redefines all occurring entries.

Adding new types based on copies of old ones is a whole different ball of wax. Here you would need a duplication operator, which is currently not supported by the emc_opls.pl script. That is not to say, that one could not add the functionality.

3

Hi Pieter,

After playing around with this for a bit today, I think I have figured out a solution to my issue. Essentially, what I want to do is create a new subtype for the different versions of each of the “basis” types in OPLS (i.e. CT, CA, HC, HA, OS, CZ, NZ, etc.) to account for the fact that some of the different versions listed in oplsaa.par have not only differing partial charges, but also differing Lennard-Jones parameters. These new versions should still have the same bond, angle, torsion, and improper equivalences to the base forms for bringing in that information from the second half of oplsaa.par and oplsaa.sb, though there may need to be some additional rules defined in those four sections to ensure that the right parameters are used for the cases described in the comments of those files.

What I have done to accomplish this is such:
As you suggested, I have used the Extra section to define my new “subtypes” and include their partial charges, with the correct basis type defined.

ITEM EXTRA
# index basis type charge torsions … (index[:occurence])
54 ct ct3 -0.180
57 ct ct2 -0.120
60 hc hc 0.060
64 ct ct4 0.115
64.1 ct ct4 0.230 # for the central sp3 C(c)(c)(C)(C) carbon connecting two benzene rings
145 ca ca -0.115
146 ha ha 0.115
166 ca cah 0.100
167 oh oha -0.530
168 ho hoa 0.430
177 os os -0.170
199 ca cao 0.085
260 ca caz 0.035
261 cz cz 0.395
262 nz nz -0.430
ITEM END

Then, for those subtypes which have different L-J parameters than the basis type, I have marked in the Equivalence that their pair and incr info should come from their own subtype, not the basis atom. For all of these new subtypes, I have noted that the bonded parameters should all come from the basis atom.

ITEM EQUIVALENCE
# type pair incr bond angle torsion improper
ct3 ct3 ct3 ct ct ct ct # new LJ
ct2 ct2 ct2 ct ct ct ct # new LJ
hc hc hc hc hc hc hc
ct4 ct4 ct4 ct ct ct ct # new LJ
cah ca ca ca ca ca ca
oha oha oha oh oh oh oh # new LJ
hoa ho ho ho ho ho ho
os os os os os os os
cao ca ca ca ca ca ca
caz ca ca ca ca ca ca
cz cz cz cz cz cz cz # new LJ
nz nz nz nz nz nz nz # new LJ
ITEM END

Finally, I have manually included the new LJ paramters for the indicated types in the NONBOND section. I have also added a LJ-sigma for the alcohol hydrogen so as to not trip an error in EMC when building.

ITEM NONBOND
# type1 type2 sigma epsilon
ct3 ct3 3.550 0.066 # 54
ct2 ct2 3.510 0.066 # 57
ct4 ct4 3.500 0.066 # 64
oha oha 3.070 0.170 # 167
cz cz 3.650 0.150 # 261
nz nz 3.200 0.170 # 262
ho ho 0.945 0.0 # From original opls-aa.define
ITEM END

This opls-pim-aa.define file, which I have attached, can be converted into a force field with an unmodified convert.sh, and then correctly builds the PIM-1 dimer with the appropriate atom types and nonbonded parameters. There is only one bonded (torsional) parameter not accounted for, [ca-ca-ct-ca], which does not appear in the OPLS-AA force field; tentatively, I have set this to zero.

For future expansion of OPLS/2020 implementation in EMC, I think that the appropriate strategy may be:

  1. Identify the atom index that EMC chooses for each basis type.
  2. Create subtypes in Masses for each type that has different Lennard-Jones parameters.
  3. Use Extra to set the partial charges of all atom types with different charges from the basis atom.
  4. Use Nonbond to set the Lennard-Jones parameters of all atom types with different parameters from the basis atom.

This seems like something that could be pretty readily automated, and would not require modification to the EMC source code. For the time being I don’t have the capacity to do this myself, however in the coming months I think that I will probably try to build a master .define file that works for the OPLS/2020 field and share it here.

Thank you for the tip on my misuse of the Repair section that helped me get on the right track to do this!

Sam

opls-pim-aa.define (6.0 KB)
opls-pim-aa.prm (7.9 KB)
opls-pim-aa.top (1.6 KB)

4

Hi Sam,

Thank you for the extended report. I guess it wasn’t clear to me, that you intended to create new types. Your description of the route to solution is indeed the proper one. It makes use of the concept of equivalences, which saves a great number of entries in subsequent bonded parameter tables. I would like to note, that – if you so feel inclined – you could look into the concept of bond increments, which is used for assigning partial charges. It negates the need to have a partial charge per type, which generalizes the eventual force field. OPLS is the only force field within the EMC family of force fields, which does not use bond increments, because it was one of the first force fields I added and did not end up using much. Let me know, if you need support.