How to define user force fields with EMC

EMC
1

I’ve been using Enhanced Monte Carlo (EMC) to create polymer membranes system as an input for LAMMPS , and up until now, the process has been rather simple. I have successfully created Poly norbornene system with butyl quaternary ammonium as it’s side chain. the system also has equal number of Hydroxide ions as the quaternary ammoniums to neutralize the system .It has certain amount of water molecules as well. However the resultant polymer.data file does not have different atom types for the carbons in the linking site with other monomers, in norbornene ring, the carbons in butyl side chain and the carbons in quaternary ammonium. Whereas in literature the difference in these type carbons has been reported. Even the partial charges for the atom types in the EMC output file is different from what has been reported.
the following is the EMC .esh script that I have been using, where it generates polymer chain with 2 monomers, 2 hydroxide ions and 6 water molecules

#!/usr/bin/env emc.pl
#

# OPTIONS
# -------------------------
ITEM OPTIONS
field           opls/2024/opls-aa
number          true
density         0.95
ewald           true
replace         true
pdb_licorice    true
emc_execute     true

ITEM END

ITEM GROUPS
# monomer with built-in [N+] is good
monomer     C1(C(C2C(CC1C2)CCCC[N+](C)(C)C)(*))(*), 1, monomer:2
terminator  *C, 1, monomer:1, 1, monomer:2
OH_ion      [O-]
# Prefer an OPLS water template if accessible; if not, keep this and set charges in LAMMPS
Water       O
ITEM END

ITEM CLUSTERS
polymer     random, 1
OH_ion      OH_ion, 2
Water       Water, 6
ITEM END

ITEM POLYMERS
polymer
2           monomer, 2, terminator, 2
ITEM END

the above script generated polymer.data which can be used as input for LAMMPS

LAMMPS output created by EMC v9.4.4, build Jul 31 2025 20:33:49

         112  atoms
         107  bonds
         198  angles
         315  dihedrals

           7  atom types
           5  bond types
           7  angle types
           7  dihedral types

             0    10.09696439  xlo xhi
             0    10.09696439  ylo yhi
             0    10.09696439  zlo zhi

Masses

       1   12.01100  # c4
       2    1.00790  # h1
       3    1.00790  # h1o
       4    1.00790  # h1w
       5   14.00670  # n4
       6   15.99940  # o1h-
       7   15.99940  # o2w

 Pair Coeffs

	  1    0.06600    3.50000  # c4,c4
	  2    0.03000    2.50000  # h1,h1
	  3    0.00000    2.50000  # h1o,h1o
	  4    0.00000    2.50000  # h1w,h1w
	  5    0.17000    3.25000  # n4,n4
	  6    0.25000    3.20000  # o1h-,o1h-
	  7    0.15540    3.16557  # o2w,o2w

 Bond Coeffs

	1  268.00000    1.52900  # c4,c4
	2  340.00000    1.09000  # c4,h1
	3  367.00000    1.47100  # c4,n4
	4  553.00000    0.94500  # h1o,o1h-
	5  600.00000    0.95720  # h1w,o2w

Angle Coeffs

	1   58.35000  112.70000  # c4,c4,c4
	2   37.50000  110.70000  # c4,c4,h1
	3   80.00000  111.20000  # c4,c4,n4
	4   50.00000  113.00000  # c4,n4,c4
	5   33.00000  107.80000  # h1,c4,h1
	6   35.00000  109.50000  # h1,c4,n4
	7   75.00000  109.50000  # h1w,o2w,h1w

Dihedral Coeffs

	1    0.70000    0.35000    0.05000    0.40000    0.00000  # c4,c4,c4,c4
	2    0.15000   -0.45000    0.00000    0.60000    0.00000  # c4,c4,c4,h1
	3    1.37950    0.63850    0.22900    0.97000    0.00000  # c4,c4,c4,n4
	4    0.72710    0.32310    0.12380    0.52780    0.00000  # c4,c4,n4,c4
	5    0.15085   -0.45255    0.00000    0.60340    0.00000  # c4,n4,c4,h1
	6    0.15000   -0.45000    0.00000    0.60000    0.00000  # h1,c4,c4,h1
	7    0.19200   -0.57600    0.00000    0.76800    0.00000  # h1,c4,c4,n4

The attached is the schematic diagram of the poly norbornene system with different carbon atom types
And also to test the EMC .esh script I have generated the Poly styrene system as well, due to the presence of aromatic ring there are different carbon atom types. This shows the EMC ability to differentiate the carbon atom types.
Ultimately I need to generate a polymer system with different carbon atoms for the linkers,for the carbons in the ring, in the side chain and finally in the quaternary ammonium. Please guide me for building such system

image
image720×460 70.6 KB

Error while running EMC setup.esh script
2

Dear user,

Thank you for the report. Can you share with us the cited literature and the expected typing and partial charges? Thanks in advance.

3

Dear user,

I assume you added some types to your provided image. Original OPLS-AA does not have any CB1, CB2, or CB3 types. I assume these are derivatives of the CB type, which typically is used for aromatic SP2 carbons. Please note that EMC does not use the original OPLS-AA types, but has redefined these types to be more consistent and transferable. A translation table can be found at $EMC_ROOT/field/opls/2024/src/opls-aa.xref. You state the following:

Could you provide the mentioned literature? This would not only include the necessary types, but also the partial charges.

4

https://www.nature.com/articles/s41467-024-55621-z

Thank you so much for your reply, it really means alot.

In this paper if you check the supplementary information for simulation section, you’ll see the PBBNB structure with different carbon atoms

image
image1948×1516 226 KB

the screenshot of the partial charges is

image
image1948×1516 176 KB

thanks again

5

Dear user,

I checked the literature you cited and see that the force field you refer to, is a force field specifically developed for these type of systems, reflected in the cited paper by Pádua:

  • Canongia Lopes, J. N., Deschamps, J. & Pádua, A. A. H. Modeling Ionic Liquids Using a Systematic All-Atom Force Field. Journal of Physical Chemistry B 108, 2038–2047 (2004).

It uses the OPLS form, but not necessarily the OPLS parameters. This means you would have to create the force field files for this application using the tabulated parameters. I also noticed, that under ’ Non-bonded Interaction Parameters’ there is no mention of inclusion of 1-4 nonbonds with a factor of 1/2, which is standard for OPLS. 1-4 nonbonds are the atoms that are separated by three bonds, i.e. that partake in torsion calculations. However, judging from the torsional parameters, I would assume the 1-4 nonbonds should be included, since the torsional parameters resemble the standard OPLS-AA parameters. It would be a question to the authors of the Nature Communication paper.

Having said all that, you can create and include this specific force field by using ITEM FIELD. Alternatively, you can create the separate force field file in the same manner as the provided $EMC_ROOT/field/opls/2024/opls-aa.[prm|top].

I have chosen to do the latter. After inspection of the provided force field parameters, it becomes clear that the only difference between the CB types is their charge. I chose to collapse these onto one type as not to make the force field needlessly complicated. The same applies to HCs and HN when compared to HC. Charges will be addressed through bond increments. As a basis, I used norbornene.define for my definitions, which I compiled into norbornene.top and norbornene.prm.

norbornene.define (3.9 KB)
norbornene.prm (2.8 KB)
norbornene.top (1007 Bytes)

I altered my setup.esh such as to include this user force field

setup.esh (632 Bytes)

The resulting setup.esh looks like

#!/usr/bin/env emc.pl

# Options section

ITEM	OPTIONS

field_location	.
field		norbornene
system_charge	false

temperature	300
density		0.95
ntotal		2000
mass		true
focus		polymer

pdb_licorice	true
emc_execute	true
replace		true

ITEM	END	# OPTIONS

# Groups section

ITEM	GROUPS

monomer		*C1C(C2C(CC1C2)CCCC[N+](C)(C)C)*.[Br-], 1, monomer:2
terminator	*C, 1, monomer:1, 1, monomer:2

ion		[Br-]

water		O

ITEM	END	# GROUPS

# Clusters section

ITEM	CLUSTERS

polymer		random, 0.1
water		water, 0.9

ion		ion, 0

ITEM	END	# CLUSTERS

# Polymers section

ITEM	POLYMERS

polymer
1		monomer, 5, terminator, 2

ITEM	END	# POLYMERS

Comparing the results with standard OPLS-AA, one sees that the differences arise from the ammonium graft of the molecule, which has a nominal charge of +0.8. To compensate, the bromium ion has a counter charge of -0.8. Partial charges of the chemical environment also changed according to Pádua. All these alterations are reflected by the provided user force field.

To reiterate: your original statement, that you expect types CB1, CB2, and CB3 to be part of OPLS-AA is incorrect since the force field as presented by the Nature Communications paper is a user derivation of OPLS-AA. This creates the need for defining your own force field. Furthermore, using three different types to define partial charges only makes sense, when their Lennard-Jones parameters are different. This is not the case here. Typically, different partial charges are controlled by using the correct bond increments instead. These can be found under ITEM INCREMENT in norbornene.prm.

The original OPLS-AA typing gives an integer charge for the ammonium group. As a result thereof, the partial charges of the chemical environment around this group are different. Pádua reports an improved set of parameters, thus deviating from standard OPLS-AA.

6

Thank you so much for your reply.

As someone who is new to this, I’m just curious how you generated the forcefield files by yourself. Until now I thought that through EMC only we could generate FF parameters.
It would be a great help for me to learn that. Thank you.

7

Dear user,

Force fields typically consist of two parts:

  1. A set of parameters
  2. A set of rules describing how to assign parameters

EMC separates these two components into two different files – a .prm and a .top – although they could be combined into one file as well.

The set of parameters describes force field settings by sections identified by DEFINE, MASS, EQUIVALENCE, INCREMENT, NONBOND, BOND, ANGLE, TORSION, and IMPROPER. You can check out how parameters are expressed by comparing page 34 f.f. of the supplement of the Nature Communications paper with norbornene.prm.

Chemical environment description in the form of typing rules are described by sections identified by PRECEDENCE and RULES. Rules describe the minimal chemical environment needed to destinguish each type, while precedence controls the order in which to interpret the given rules. You can study the norbornene.top file to understand the concept of rules better.

The EMC manual describes the above in somewhat more detail under `Simulation Setup => File Formats => Field File’.

As an aside, please note that force fields are not equal force fields, meaning that one should not blindly mix and match parameters. Parameters should only be combined if all parameters have been obtained by the same in itself consistent parameterization scheme.

8

Thank you so much for the reply.

after reading your reply I got some clarity about the concept of force fields. However, I still did not get how to generate the files like

norbornene.define (3.9 KB)
norbornene.prm (2.8 KB)
norbornene.top (1007 Bytes)

I downloaded the files from your reply and run the .esh script. It ran successfully. But I just want to learn how to create the above files from scratch.

I apologize if I’m asking too many questions.Thank you so much for your time.

9

Hi Veld,

hope you are doing well

I’m trying to build the same polymer, But instead of Br ion I want to use Hydroxide ions (OH-) .I edited the file norbornene.define specifying OH- ions and deleting Br ions and was able to generate .prm and .top file too and ran the setup.esh but getting errors.

norbornene_OH.top (414 Bytes)

norbornene_OH.prm (2.1 KB)

norbornene_OH.define (4.2 KB)

the setup.esh file That I have used looks like

Summary 
#!/usr/bin/env emc.pl

# Options section

ITEM OPTIONS

field_location .
field norbornene
system_charge false

temperature 300
density 0.95
ntotal 2000
mass true
focus polymer

pdb_licorice true
emc_execute true
replace true

ITEM END # OPTIONS

# Groups section

ITEM GROUPS

monomer *C1C(C2C(CC1C2)CCCC[N+](C)(C)C)*, 1, monomer:2
terminator *C, 1, monomer:1, 1, monomer:2

ion [O-]

water O

ITEM END # GROUPS

# Clusters section

ITEM CLUSTERS

polymer random, 1
water water, 150

ion ion,50

ITEM END # CLUSTERS

# Polymers section

ITEM POLYMERS

polymer
1 monomer, 50, terminator, 2

ITEM END # POLYMERS

Please guide me

Thank you
Amala