Missing OPLS-AA parameters for poly styrene-butadiene (PSB) rubber

Hello,

I am building cis-polybutadiene and was missing a torsion parameter in the opls-aa.prm for cm-ct-ct-cm (=CH-CH2-CH2-CH=)
Following the paper of Fröhlich 2013 (https://onlinelibrary.wiley.com/action/downloadSupplement?doi=10.1002%2Fmats.201300103&file=mats201300103-sm-0001-SuppTables.pdf), I suggest adding this line to file field/opls/2012/opls-aa.prm:
cm ct ct cm 1.3 1 0 -0.05 2 0 0.2 3 0

Also, using the OPLS-UA is not preserving the stereochemistry of cis-butadiene even though I put SMILES code “*C/C=C\C*”, only OPLS-AA is preserving. Is there any torsion parameters missing in OPLS-UA too ?

Best,

Eddy BARRAUD

I am also modeling styrene butadiene rubber, and OPLS-aa 2012 is missing the definition of styrene cycles. If I try the OPLS-aa 2024 then it is missing torsion parameters:
Warning: no torsion coefficients found for [c4:c4:c3=:h1=].
Warning: no torsion coefficients found for [c3a:c4:c4:c3=].
Warning: no torsion coefficients found for [c3=:c4:c4:c3=].

Dear user,

I would suggest using the 2024 instead of the 2012 version of OPLS. The 2012 version uses direct charge assignment, while the 2024 version incorporates bond increments, which makes it a more general force field. I checked out Fröhlich and Sewell, Macromol. Theor. Sim. 2013, 22, 344–353 and noticed, that their suggested dihedral contribution is based on the standard CT-CT-CT-CT or – in 2024 vernacular – c4-c4-c4-c4 (see $EMC_ROOT/field/opls/2024/src/opls-aa.xref for naming convention). I added a generalized rule to the 2024 version (c3*-c4-c4-c*) using the parameters for c4-c4-c4-c4. The resulting 2024 version is given by:

opls-aa.prm (66.7 KB)
opls-aa.top (20.1 KB)

To answer your observation with respect to stereochemistry: EMC forces initial placement of atoms in the correct stereochemistry during the growth process. However, based on the force field, this initial placement could be destroyed during the relaxation step when the accompanying intramolecular structure of the immediate surrounding displays force constants which are unable of maintaining the desired stereochemistry. It seems that this is the case for OPLS-UA.

On a scientific note: styrene-butadiene rubbers seem to consist of a random mix of monomers on the backbone (see Styrene-butadiene - Wikipedia). To simulate this with EMC, one could do the following:

#!/usr/bin/env emc.pl
#
#  script:	setup.esh
#  author:	Pieter J. in 't Veld
#  date:	July 16, 2026
#  purpose:	Example for polymers using EMC setup chemistry mode
#

# Options section

ITEM	OPTIONS

replace		true
field		2024/opls-aa
number		true
ntotal		1000
density		0.1
focus		true
pdb_licorice	true
emc_execute	true

ITEM	END	# OPTIONS

# Groups section

ITEM	GROUPS

# monomers

bd12		*C(C=C)C*, &
		1,bd12:2, 2,bd12:1, &
		1,bd14:2, 2,bd14:1 &
		1,styrene:2, 2,styrene:1, &
		1,term:1, 2,term:1

bd14		*C/C=C\C*, &
		1,bd12:2, 2,bd12:1, &
		1,bd14:2, 2,bd14:1 &
		1,styrene:2, 2,styrene:1, &
		1,term:1, 2,term:1

styrene		*C(c1ccccc1)C*, &
		1,bd12:2, 2,bd12:1, &
		1,bd14:2, 2,bd14:1 &
		1,styrene:2, 2,styrene:1, &
		1,term:1, 2,term:1

term		*C

# polymer

poly		random

ITEM	END	# GROUPS

# Clusters section

ITEM	CLUSTERS

polymer		poly,10

ITEM	END	# CLUSTERS

# Polymers section

ITEM	POLYMERS

poly
log:avg=20,pdi=1.01	bd12:bd14:styrene=1:1:2,$@, term,2

ITEM	END	# POLYMERS

In the above example I chose a low density of 0.1 g/cc for no other reason than a quick test building only 10 molecules. A standard simulation would be built at 0.94 g/cc (as stated below) and should include over a 100 molecules. The latter serves as a rule of thumb to avoid system size effects. The choice of defining a group as a polymer instead of a cluster – i.e. the word random appears in a group and not in a cluster – allows for randomly picking from a distribution, rather than predefining the number of clusters a priori. Check out the difference in the resulting build.emc.

I have added a log-normal distribution to EMC Setup. I also included the above mentioned force field additions. The result I uploaded to https://montecarlo.sourceforge.net/ for the latest version of EMC.

1 Like

Thank you very much for all that @veld. Sorry for the delayed answer. I will investigate this model later on. First, benchmarking against cis-polybutadiene because I have really good models for it. The input script is really welcome and similar to what I had. Are you voluntarily setting a density of 0.1 ? I tried a density of 0.94, close to the target density, and it seemed to build without many difficulties. I have the molar percentage of each monomer, but I am not sure of the log-normal distribution.

This was my code before:

#!/usr/bin/env emc_setup.pl

ITEM OPTIONS

seed            34220
field           opls/2012/opls-aa
density         0.94
temperature     300
ntotal          20000
replace         true
build_dir       .
radius          5.0
nrelax          100
pressure        1.0
precision       0.00001
polymer         bias=accumulative,niterations=5000
emc_run         nblocks=100,ncycles=2000,nequil=1000,clusters=all,groups=all,sites=all
emc_traject     append=true,frequency=1

ITEM END

ITEM GROUPS

# Atactic styrene repeat unit

S   *C(c1ccccc1)C*,1,S:2,1,BV:2,1,BT:2,1,BC:2,2,BV:1,2,BT:1,2,BC:1

# 1,2-vinyl butadiene repeat unit

BV  *C(C=C)C*,1,BV:2,1,BT:2,1,BC:2,2,BT:1,2,BC:1

# trans 1,4 butadiene repeat unit

BT  *C/C=C/C*,1,BT:2,1,BC:2,2,BC:1

# cis 1,4 butadiene repeat unit

BC  *C/C=C\C*,1,BC:2,1,BV:2,1,BT:2,2,BT:1,2,BV:1

# terminator CH3

T   *C,1,S:1,1,S:2,1,BV:1,1,BV:2,1,BT:1,1,BT:2,1,BC:1,1,BC:2

ITEM END

ITEM CLUSTERS

sbr     random,1

ITEM END

ITEM POLYMERS

sbr
1       S:BV:BT:BC=0.158:0.362:0.260:0.220,90,T,2

ITEM END

Best,

Eddy BARRAUD

Dear user,

Thank you for sharing your input script. I took the liberty of answering your question by altering my response in the message above. There were no other reasons than illustrative ones for choosing 0.1 g/cc. You can swap this out with the standard condensed matter condition of 0.94 g/cc without issue. The energetic building approach of EMC will intuitively avoid possible ring concatenations when incorporating styrene monomers under condensed matter conditions.

Subsequently, I would like to emphasize the difference between defining a polymer under ITEM CLUSTERS as opposed to defining a polymer under ITEM GROUPS. The former will result in hard-coding the number of clusters per distribution entry before building, meaning that any probability resulting in a number of chains smaller than 1 will not be build. Polymers defined under ITEM GROUPS on the other hand will form a distribution from which – during building – molecules are drawn based on the provided distribution, i.e. in a Monte Carlo fashion. This means that long chains with a low probability could be built when defining a polymer in ITEM GROUPS, while they will be excluded when defining it in ITEM CLUSTERS.

I am using the log-normal approach as a quick way to introduce a distribution. You can also provide your own distribution by explicitly defining the probability for each chain entry under ITEM POLYMERS.

As an aside, you set polymer bias=accumulative in ITEM OPTIONS in your above response. This setting only has effect when defining your polymers through ITEM GROUPS.

I hope this helps understanding the nuances of how to build polymers and incorporate their distributions.

Great, thank you for the explanations. Indeed, ITEM GROUPS is clearly the best way to build realistic polymer chains with varying chain lengths. I am still discovering the EMC software, and your messages helped me better understand it. It would be awesome to have the code open source so that we could see the Monte Carlo algorithm under the hood. Still, this free software is amazing!

Thank you for your kind words. Good to hear, that you find EMC helpful for your research. The ITEM CLUSTERS solution is actually a hold-over from the original solution. This solution proved to have its limitations, esp. when for polymer distributions consisting of prepolymers with distributions, hence the newer ITEM GROUPS approach. Good luck with your research, and please let us know, if you end up with further questions!