Regarding dihedrals in EMC input file

EMC
1

Dear Pieter,

  • is it possible to specify a weight value for the 1-4 pairwise Lennard-Jones and Coulombic interactions in the EMC setup script (for instance in the ITEM FIELD section) similar to the weight value specified by the “special_bonds” command in LAMMPS? Several OPLS-derived FFs specify a value of 0.5 for the 1-4 interactions which is neither completely on nor completely off.

  • I have converted the parameters of a Ryckaert-Bellemans dihedral (6 in total since the series of cos(phi) goes up to the fifth power) to the corresponding five Fourier dihedral terms plus the constant term. The generated output LAMMPS files write the data for a multi/harmonic dihedral style which includes only the first five terms. Is EMC reading the entire line with all 6 terms when creating the initial configurations and can I modify the output style to nharmonic so that all 6 terms are written to the LAMMPS parameter file?

Thanks, Evangelos

2

Dear Evangelos,

I will answer your question in the order you provided.

  • EMC deals with 1-4 pair interactions through the chosen force field. Force fields supporting a pairwise 1-4 interaction are: CHARMM, Mike, SDK, OPLS, and TraPPE, for all of which the 1-4 contribution is set to 0.5. EMC does not support the 3/8 contribution as needed by AMBER. EMC translates the 1-4 condition correctly to LAMMPS for all the aforementioned force fields. The 1-4 contribution is controlled by the ITEM DEFINE paragraph in the .prm parameter file by the keywords NONBONDED and PAIR14. NONBONDED sets the exclusion of pair, e.g. the pairs 1-2, 1-3, and 1-4 are excluded when set to 3. PAIR14 controls a separate 1-4 pair interaction function by either EXCLUDE or INCLUDE.

  • Input of dihedral terms are according to the definition of the chosen force field. The EMC manual explains all of the choices per force field. Alternatively, you can check the entries for dihedrals in the force field parameter files as can be found in the $EMC_ROOT/field directory and its derivatives. For instance – when considering CHARMM – the dihedral form is an expansion of k_j (1 ± cos(j phi + delta_j)). Internally, EMC converts this form in a cosine series, such that the expansion function as is available in LAMMPS can be used. Having said that, you are correct, that the multi/harmonic dihedral style stops at order 5. You could use CHARMM as an EMC force field function set to also include order 6 harmonics. If I am correct, CHARMM and OPLS have in essence identical functional forms.

3

Dear Pieter,

Thanks a lot for the detailed response.

By adding the keyword PAIR14 with INCLUDE and setting NONBONDED equal to 3 in the attached script I have indeed created systems with special_bonds of 0 0 0.5 similar to OPLS settings.

standard.esh (1.8 KB)

Still it is not clear to me how to modify the torsion keyword in a custom FF to take care of a Ryckaert-Bellemans dihedral. I have the impression that the allowed input format is a series of Fourier something along the following though shorter:

c* c* c* c* -3.764288196 0 0 13.22827194 1 0 5.349242785 2 0 2.69991817 3 0 -1.070633166 4 0 -1.141931683 5 0

Is EMC reading all parameters & using them for the initial configuration but simply not writing the info to the lammps params file? Can I use a different functional form for the ITEM TORSION section?

Thanks, Evangelos

4

Dear Evangelos,

I noticed, that the manual is missing the number of contributions for torsions taken into account. This number is 4. Input commences through the lines provided in the ITEM TORSION section, e.g. using your paragraph

ITEM	TORSION

# type1 type2	type3	type3	k	n	delta	[...]

f	c*	c*	f	1	1	1 
*	*	*	*	1	1	1 

ITEM	END

means that for all type combinations – since you provided a full wildcard line with * * * * – you set k = 1, n = 1, and delta = 1 for a formula, where

e_torsion[i] = k[i]*(1 - (-1)^n[i]*cos(n_[i]*phi + delta[i]*pi / 180)),

with the above becomes

e_torsion = 1 + cos(phi + pi/180)

for each torsion contribution, where pi = 3.141592.... I assume, that that is not what you want. I would set delta = 0 or perhaps delta = 180.

Each line allows you to supply up to four contributions, where n = 6 is the maximum. Unfortunately, due to the limitation of multiharmonic Tailor expansion used when porting to LAMMPS, only values up to n = 5 are communicated. A valid (fictitious) line could be

*	*	*	*	1.2	1	0	2.0	2	180	1.0	4	0

etc. This results in a sum of

e_torsion = 1.2*(1+cos(phi)) + 2.0*(1-cos(2*phi+pi)) + 1.0*(1-cos(4*phi))

for the torsion contribution. Having said that, this corresponds to the entry line you suggest in your post.

Torsion functional forms are non-changeable, since they are set by the choice of force field. You can go to a different form by changing force field representation. However, this might also change other nonbonded and bonded force field contributions, e.g. an OPLS force field representation significantly differs from a class II CFF representation. Details of force field definitions form the reason for choosing static force field representations. This minimizes operation errors.

1 Like
5

Thanks for the detailed and prompt answer. Everything is crystal clear !