JOYCE Tutorial. Part (II / III)

Author: Anna Piras, Pablo M. Martinez

Keywords: Tutorial, FF, JOYCE

Once we have obtained the QM Hessian, we are ready to begin with the JOYCE parametrization.

Step 0. Atom types and internal coordinates.

Create a folder and call it Step0. You must include the following files:

Step0

pyr.step0.top

joyce.pyr.inp

The first one is an empty GROMACS topology file. Download it here. It will define the atom types discussed in part I/III and help JOYCE to propose the necessary internal coordinates for the parametrization. The second file (joyce.pyr.inp) contains the instructions for the fitting. In Step0, it looks like this:

Pyridine - step0                            ; a comment

$title Pyridine                             ; the name
$equil ../../QMdata/opt+freq.fcc            ; the location of the QM Hessian
$forcefield gromacs pyr.step0.top           ; the type and location of topology
$zero 1.d-16                                ; precision
$whess 5000. 2500.0                         ; weight of diagonal and off-diagonal Hessian elements
$generate                                   ; instruction

You can also download it from here.

Run JOYCE as described in the manual.

JOYCE will then generate the file generated.IC.txt with the suggested internal coordinates. It will also generate a .xyz file containing the optimized geometry for visualization , although it will not be used for the fitting.

Step 1. Fitting the stiff (harmonic) ICs.

A new folder named Step1 is created, containing a starting topology file (pyr.step1.top) and a Joyce input. The internal coordinates generated in Step0 are inserted in the topology.

Before proceding with the fitting, let us review what type of ICs have been generated.

Stretchings. These are guessed from the molecular topology. No changes are needed.
Bendings. These are also guessed from the molecular topology. No changes are needed.
Torsions. You will find three types of torsion
- Proper stiff torsions. For instance, 4 consecutive carbon atoms within a ring. You should remove all of them involving a terminal H, i.e. (C-C-C-H) because the next type better captures the H motion.
- Improper (star) stiff torsions. They control the out of plane flexibility of a given atom with respect to the plane defined by the other three. Thus, the dihedral is not linear. No changes are needed here.
Pairs. These should be removed by default for a rigid system such as this one.
Exclusions. No changes needed.

All the molecular motion will be captured by a combination of these terms. For a simple molecule, i.e. pyridine, not a lot of additional work is required on the automatically generated internal coordinates.

Highly encouraged activity

Even for such a simple system it is advisable to look at every IC and understand its meaning. It becomes necessary for more complex molecules

You can download the ordered topology file from here.

Then, we begin the fitting. Change the JOYCE input to:

Pyridine - step1            ; a comment

$title pyridine                             ; the name
$equil ../../opt+freq.fcc                   ; the location of the QM Hessian
$forcefield gromacs pyr.step1.top           ; the type and location of topology
$zero 1.d-16                                ; precision
$whess 5000. 3000.0                         ; weight of diagonal and off-diagonal Hessian elements

Removing the $generate keyword indicates to JOYCE that it should fit the specified ICs. It disregards any value in the starting topology, just takes the ICs defined and considers them for the fitting procedure.

You can run JOYCE and check the resulting joyce.new.top file. Equilibrium values are extracted from the Gaussian optimized structure, and the force constants are the output from JOYCE. You will see that some of them are negative. A negative force constant is sign of overfitting. You can check the JOYCE output at the .out file. Search for "RICOs", there you will see the comparison between the expected and the actual number of ICs defined. Remember that the ideal scenario is having only $3N-6$ ICs. The fact that we obtain negative force constants stems from this overfitting issue. We should reduce the number of independent ICs as much as possible and think of a) dropping unnecessary ICs and/or b) introduce constraints between ICs.

JOYCE automatically generates the suggdeps.txt file which can be copied right after the $whess keyword in the joyce.pyr.inp file. Thus, the force constants of some ICs are forced to be the same by symmetry arguments.

The final version of the JOYCE

$equil $forcefield $zero $whess $gracefreq $dependence                        $end

input file will be: name="__codelineno-2-1" href="#__codelineno-2-1">$title Pyridine - Step 1 - ../../QMdata/opt+freq.fcc gromacs pyr.step1.top 1.d-15 5000. 2500.0 joyce.pyr_freqchk.agr 1.2 5 = 1*1.d0 ; n=c1 = n=c1 6 = 2*1.d0 ; c1=c2 = c1=c2 4 = 3*1.d0 ; c2=c3 = c2=c3 11 = 7*1.d0 ; c1-h1 = c1-h1 10 = 8*1.d0 ; c2-h2 = c2-h2 24 = 12*1.d0 ; n=c1=c2 = n=c1=c2 25 = 14*1.d0 ; n=c1-h1 = n=c1-h1 21 = 15*1.d0 ; c1=c2=c3 = c1=c2=c3 27 = 16*1.d0 ; c2=c1-h1 = c2=c1-h1 26 = 17*1.d0 ; c1=c2-h2 = c1=c2-h2 23 = 19*1.d0 ; c3=c2-h2 = c3=c2-h2 22 = 20*1.d0 ; c2=c3-h3 = c2=c3-h3 42 = 28*1.d0 ; n=c1=c2=c3 = n=c1=c2=c3 31 = 29*1.d0 ; c1=n=c1=c2 = c1=n=c1=c2 46 = 30*1.d0 ; n=c1=c2-h2 = n=c1=c2-h2 33 = 32*1.d0 ; c1=n=c1-h1 = c1=n=c1-h1 38 = 34*1.d0 ; c1=c2=c3=c2 = c1=c2=c3=c2 44 = 35*1.d0 ; h1-c1=c2=c3 = h1-c1=c2=c3 43 = 36*1.d0 ; c1=c2=c3-h3 = c1=c2=c3-h3 47 = 37*1.d0 ; h1-c1=c2-h2 = h1-c1=c2-h2 40 = 39*1.d0 ; h2-c2=c3=c2 = h2-c2=c3=c2 45 = 41*1.d0 ; h2-c2=c3-h3 = h2-c2=c3-h3 52 = 48*1.d0 ; h1--n--c2=c1 = h1--n--c2=c1

The gracefreq keyword is used for the force field validation (see here for details).

You may download the final result here.

Continue with the tutorial here to validate the final QMD-FF.