Implement custom force field via the FIELD keyword

Dear all,

I am trying to create an initial configuration of a mixture of water and tetrafluoroethylene using a custom parameterization of the OPLS force field. I am able to create an initial configuration using the “standard” OPLS FF which is shipped together with EMC. The custom parameters are defined in the attached standard.esh file in the “ITEM FIELD” section. When running the command emc_setup.pl standard.esh I get a message of successful completion and the build.emc file is created. When I run the emc_linux build.emc command including also the standard.top file in the same directory I get the following error message:

Error: core/chemistry/ring.c:651 ChemistryRings:
Found nrings != nlinks (0 != 1); increase scanning depth?
Program aborted.

Could you give me any advice on how to properly format the “ITEM FIELD” section and/or modify the standard.top file? Is there a utility to create automatically the .top file? Lastly, I would appreciate any comment on the formatting of the dihedral parameters in the FIELD section (for instance what is the expected input functional form).

Thanks, Evangelos

standard.esh (2.5 KB)
standard.top (725 Bytes)

Hi Evangelos,

For illustrative purposes, your original script read as follows:

#!/usr/bin/env emc_setup.pl

# Options

ITEM	OPTIONS

ntotal		1000
temperature	353.15
density		0.1
field_name      standard
field_location  .
mass		true
replace		true
build_dir	.

ITEM	END

# Shorthand


# Groups

ITEM    GROUPS

water         O
tfe           *C2(F)(F)C2(F)(F)*,1,tfe:2
htfe          *C2(F)(F)C3(F)(F)F,1,tfe:1,1,tfe:2

ITEM    END

# Clusters

ITEM    CLUSTERS

water           water,0.054523562
poly            alternate,1

ITEM    END

# Polymers

ITEM    POLYMERS

poly
100             tfe,98,htfe,2

ITEM    END

# Force Field specification according to the B3LYP/6-13G* data

ITEM    FIELD

# Definitions

ITEM    DEFINE

FFNAME  STANDARD
FFTYPE  ATOMISTIC
LENGTH  ANGSTROM
ENERGY  KCAL/MOL
DENSITY G/CC
MIX     BERTHELOT
CUTOFF  9.5
NBONDED 2
ANGLE   COMPLETE
TORSION COMPLETE

ITEM    END

# Masses

ITEM    MASS

# type  mass    element ncons   charge  comment

C2      12.011  C       4       0.378   C atom for -CF2-
C3      12.011  C       4       0.567   C atom for -CFE
F       18.998  F       1       -0.189  F atom
O       15.9994 O       2       -0.8476 O atom for water       
H       1.00790 H       1       0.4238  H atom for water

ITEM    END

# Nonbond parameters

ITEM    NONBOND

# type1 type2   sigma   eps

F       F       2.94    0.061  
C2      F2      3.50    0.066  
C3      F3      3.50    0.066  
O       O       3.16556 0.15539
H       H       1.0     0.0

ITEM    END

# Bond parameters

ITEM    BOND

# type1 type2   k       l0

C2      C2      259.93  1.54
C2      C3      259.93  1.54
F       C2      304.82  1.38
F       C3      304.82  1.38
O       H       600.00  1.00

ITEM    END

# Angle parameters

ITEM    ANGLE

# type1 type2   type3   k       theta0

F       C2      F       107.64  108.72
F       C3      F       107.64  108.72
F       C2      C2      56.39   107.90
F       C2      C3      56.39   107.90
F       C3      C2      56.39   107.90
F       C3      C3      56.39   107.90
C2      C2      C2      56.39   107.90
C2      C2      C3      56.39   107.90
C2      C3      C2      56.39   107.90
C2      C3      C3      56.39   107.90
C3      C2      C3      56.39   107.90
C3      C3      C3      56.39   107.90
H       O       H       75.00   109.47

ITEM    END

# Torsion parameters

ITEM    TORSION

# type1 type2   type3   type3   k       n       delta   [...]
F       C2      C2      F       1       1       1 
F       C2      C3      F       1       1       1 
*       *       *       *       1       1       1 

ITEM    END

ITEM    END     # FIELD

As a general remark: this script mainly contains spaces as separating characters. EMC Setup and derived modules function better when you use TABs.

Regarding your error: in general, the error occurs when a ring closure cannot be found, as is the case for the cyclodextrin question on this board. Specifically, your error has to do with how you defined your SMILES under ITEM GROUPS. You defined your tfe group with tfe *C2(F)(F)C(F)(F)*. You provided an opening link in your SMILES by providing a 2, but you never provided the accompanying closing link, hence the error. If you would have wanted not to use rules, you could have defined this group by using your types: the *[c2](f)(f)c(f)(f). Note that the field module converts used types to lower case.

Regarding using ITEM FIELD: EMC comes with an included example in $EMC_ROOT/examples/setup/chemistry/field/user, where a user defined polymer force field is used. EMC Setup uses a force field module to obtain your standard.prm and standard.top. Be aware, that this force field requires a TAB to be the separator between fields. Using spaces will cause the module to function incorrectly. Both parameters and rules can be included in the ITEM FIELD section. I included partial wild cards for your parameters. With these changes, your script looks as follows:

#!/usr/bin/env emc_setup.pl

# Options section

ITEM	OPTIONS

ntotal		1000
temperature	353.15
density		1.0
field_name	standard
field_location	.
field_format	%g
mass		true
replace		true
build_dir	.
emc_execute	true

ITEM	END	# OPTIONS

# Groups section

ITEM	GROUPS

water		O
tfe		*C(F)(F)C(F)(F)*,1,tfe:2
htfe		*C(F)(F)C(F)(F)F,1,tfe:1,1,tfe:2

ITEM	END	# GROUPS

# Clusters section

ITEM	CLUSTERS

water		water,0.054523562
poly		alternate,1

ITEM	END	# CLUSTERS

# Polymers section

ITEM	POLYMERS

poly
100		tfe,98,htfe,2

ITEM	END	# POLYMERS

# Force Field specification according to the B3LYP/6-13G* data

ITEM	FIELD

# Definitions

ITEM	DEFINE

FFNAME		STANDARD
FFTYPE		ATOMISTIC
FFDEPTH		2
LENGTH		ANGSTROM
ENERGY		KCAL/MOL
DENSITY 	G/CC
MIX		BERTHELOT
CUTOFF		9.5
NBONDED 	2
ANGLE		ERROR
TORSION 	ERROR

ITEM	END

# Masses

ITEM	MASS

# type	mass	element ncons	charge	comment

c2	12.011	C	4	0	C atom for -CF2-
c3	12.011	C	4	0	C atom for -CFE
f	18.998	F	1	0	F atom
o	15.9994	O	2	0	O atom for water	
h	1.00790	H	1	0	H atom for water

ITEM	END

# Rule precedences

ITEM	PRECEDENCE

(
  (c2)
  (c3)
  (f)
  (o)
  (h)
)

ITEM	END

# Rules

ITEM	RULES

# type	charge	rule

c2	0.378	C(F)(F)(C)(C)
c3	0.567	C(F)(F)(F)(C)
f	-0.189	F(C)
o	-0.8476	[O](H)(H)
h	0.4238	H(O(H))

ITEM	END

# Nonbond parameters

ITEM	NONBOND

# type1 type2	sigma	eps

f	f	2.94	0.061	
c2	c2	3.50	0.066	
c3	c3	3.50	0.066	
o	o	3.16556	0.15539
h	h	1.0	0.0

ITEM	END

# Bond parameters

ITEM	BOND

# type1 type2	k	l0

c*	c*	259.93	1.54
f	c*	304.82	1.38
o	h	600.00	1.00

ITEM	END

# Angle parameters

ITEM	ANGLE

# type1 type2	type3	k	theta0

f	c*	f	107.64	108.72
*	c*	c*	56.39	107.90
h	o	h	75.00	109.47

ITEM	END

# Torsion parameters

ITEM	TORSION

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

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

ITEM	END

ITEM	END	# FIELD

I added ffdepth 2 since your rules have a maximum recursion depth of 2 (see your rules for water oxygen and hydrogen). I added field_format %g to change the numerical format of the resulting force field files. For completeness sake, I also attached the resulting script.

standard.esh (1.8 KB)

Thanks a lot for your prompt and detailed response Pieter! If i can contribute to the manual please let me know.
Best regards Evangelos