JOYCE Tutorial. Part (I / IV)

Author: Anna Piras, Pablo M. Martinez

Contact: pablo.martinezm@uam.es

Keywords: Tutorial, FF, JOYCE

---

Once we have obtained the QM Hessian, we are ready to begin with the JOYCE parametrization. This part is analogous to previous tutorials.

Step 0. Atom types and internal coordinates.

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

• Step0
  • but.step0.top
  • joyce.but.inp

Download the topology from here. It will define the atom types discussed in part I/IV and help JOYCE to propose the necessary internal coordinates for the parametrization. The second file (joyce.but.inp) contains the instructions for the fitting. You can also download it from here.

Run JOYCE as described in the manual.

Step 1. Fitting the stiff (harmonic) ICs.

Create a second folder named Step1, containing a starting topology file (but.step1.top) and a Joyce input. The internal coordinates generated in Step0 are inserted in the topology.

As in previous examples, we comment on the ICs that 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 one type of torsion

  • Proper flexible torsions. For instance, 4 consecutive carbon atoms. You should remove all of them involving a terminal H, i.e. (C-C-C-H) since their motion is well-governed by the combination of backbone torsions, bendings and stretchings.

• Pairs. These should be removed by default for a small system such as this one.

• Exclusions. No changes needed.

The key difference with previous cases lies on the type of dihedrals. In this case, all atom quadruplets are not rotationally constrained. These internal coodinates will be fitted though a Fourier series -- see the Manual. Owing to symmetry, we may expect these rotations to display 3 minima. For this reason, we can ask JOYCE to use a third-order cosine function in this step. This may be seen as a coarse-grained preparation for the actual fitting. In any case, direct comparison between the force constants of stretchings/bendings quickly shows how the energy of these "coarse" terms are much finer than the other ICs. For this reason, it is not mandatory to include these flexible terms at this point, but useful to remember how many flexible torsions there are. After removing the unnecessary torsions, we obtain the starting topology, which can be downloaded from here.

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

$title Butane - Step 1 -
$equil ../../QMdata/opt+freq.fcc
$forcefield gromacs but.step1.top
$zero   1.d-15
$whess  5000. 2500.0

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

The final version of the JOYCE input file will be:

$title Butane - Step 1 -
$equil ../../QMdata/opt+freq.fcc
$forcefield gromacs but.step1.top
$zero   1.d-15
$whess  5000. 2500.0
$dependence 1.2
      3 =    1*1.d0   ; C1-C2                = C1-C2               
      5 =    4*1.d0   ; C1-H1                = C1-H1               
      6 =    4*1.d0   ; C1-H1                = C1-H1               
     11 =    4*1.d0   ; C1-H1                = C1-H1               
     12 =    4*1.d0   ; C1-H1                = C1-H1               
     13 =    4*1.d0   ; C1-H1                = C1-H1               
      8 =    7*1.d0   ; C2-H2                = C2-H2               
      9 =    7*1.d0   ; C2-H2                = C2-H2               
     10 =    7*1.d0   ; C2-H2                = C2-H2               
     20 =   14*1.d0   ; C1-C2-C2             = C1-C2-C2            
     16 =   15*1.d0   ; C2-C1-H1             = C2-C1-H1            
     17 =   15*1.d0   ; C2-C1-H1             = C2-C1-H1            
     27 =   15*1.d0   ; C2-C1-H1             = C2-C1-H1            
     28 =   15*1.d0   ; C2-C1-H1             = C2-C1-H1            
     29 =   15*1.d0   ; C2-C1-H1             = C2-C1-H1            
     19 =   18*1.d0   ; C1-C2-H2             = C1-C2-H2            
     25 =   18*1.d0   ; C1-C2-H2             = C1-C2-H2            
     26 =   18*1.d0   ; C1-C2-H2             = C1-C2-H2            
     22 =   21*1.d0   ; C2-C2-H2             = C2-C2-H2            
     23 =   21*1.d0   ; C2-C2-H2             = C2-C2-H2            
     24 =   21*1.d0   ; C2-C2-H2             = C2-C2-H2            
     31 =   30*1.d0   ; H1-C1-H1             = H1-C1-H1            
     32 =   30*1.d0   ; H1-C1-H1             = H1-C1-H1            
     35 =   30*1.d0   ; H1-C1-H1             = H1-C1-H1            
     36 =   30*1.d0   ; H1-C1-H1             = H1-C1-H1            
     37 =   30*1.d0   ; H1-C1-H1             = H1-C1-H1            
     34 =   33*1.d0   ; H2-C2-H2             = H2-C2-H2            
     40 =   39*1.d0   ; H1-C1-C2-C2_n=3      = H1-C1-C2-C2_n=3     
     41 =   39*1.d0   ; H1-C1-C2-C2_n=3      = H1-C1-C2-C2_n=3     
     42 =   39*1.d0   ; H1-C1-C2-C2_n=3      = H1-C1-C2-C2_n=3     
     43 =   39*1.d0   ; H1-C1-C2-C2_n=3      = H1-C1-C2-C2_n=3     
     44 =   39*1.d0   ; H1-C1-C2-C2_n=3      = H1-C1-C2-C2_n=3     
$end 

The selection of atom types was done through symmetry arguments. Then, all ICs with the same labels must be equal. This is done to avoid overfitting, as this may produce unphysical negative force constants and, what is more, no reason can be argued to differentiate them. This will avoid unnecesary complexity in the final force field and improves its clarity.

You may download the final result here.

Continue with the tutorial here to proceed with the flexible ICs.