JOYCE Tutorial. Part (III / IV)
Author: Anna Piras, Pablo M. Martinez
Contact: pablo.martinezm@uam.es
Keywords: Tutorial, FF, JOYCE
After parametrizing the stiff ICs, we move on to the flexible ones. We recommend storing the data corresponding to the flexible scans in a separate folder than the Hessian. The procedure is based on the Frozen Internal Rotation Approximation or FIRA. The results obtained in the previous section will be fixed constant during the fitting. This is not fully accurate since the electronic cloud may suffer changes as the relative posiiton of the nuclei vary through the IC rotation. The validity of this assumption is checked in the last step of the Tutorial.
Step 2. Defining rigid scans in Gaussian
Create a folder and call it Step2. You must include the following files:
- Step2
- but.step2.top
- joyce.but.inp
Download the topology from here. Some modifications have been made from the output of Step1. Observing the PES of the terminal H scans (right), we can see they are 3-fold degenerate around a full rotation.
For this reason, only the order n=3 is specified. For the backbone C-C-C-C torsion (left) the degeneracy is broken and we include additional terms to improve the fitting quality -- see Manual.
The second file (joyce.pyr.inp) contains the instructions for the fitting. There are new sections corresponding dedicated to the fitting of flexible ICs.
-
$geom specifies the location of the reference QM data to map the PES into the selected Fourier series. The quadruplet after the semicolon must refer to the atom indexes of the topology file associated with said torsion. The -f flag avoids repeting the same quadruplet for all points withing a scan.
-
$assign contains the force constants generated for the stiff ICs in the previous step. They are written in the assign.dat file at the end of Step1. This way, we constrain those ICs to encode the same information as we fit the flexible ICs.
-
$keepff preserves the equilibrium value of the flexible ICs. The indexes run through all the flexible ICs. Notice how in the topology, we have manually set them to be 0. This is consistent with the output for Step1, as those fluctuations can be ascribed to numerical errors and have no chemical justification.
-
$scan generates an ASCII file with the resulting fitted Fourier series.
$title Butane - Step 2
$equil ../../QMdata/opt+freq.fcc
$forcefield gromacs but.step2.top
$zero 1.d-12
$whess 5000. 2500.0
$gracefreq joyce.but_freqchk.agr
$gracetors but_torschk
$geom
../../QMdata/Scan1/butane.delta_000.fcc 10.0 0.0 0.0 0.0 ; 1 2 3 4
../../QMdata/Scan1/butane.delta_030.fcc 10.0 0.0 0.0 0.0 -f
../../QMdata/Scan1/butane.delta_060.fcc 10.0 0.0 0.0 0.0 -f
../../QMdata/Scan1/butane.delta_090.fcc 10.0 0.0 0.0 0.0 -f
../../QMdata/Scan1/butane.delta_120.fcc 10.0 0.0 0.0 0.0 -f
../../QMdata/Scan1/butane.delta_150.fcc 10.0 0.0 0.0 0.0 -f
../../QMdata/Scan1/butane.delta_180.fcc 10.0 0.0 0.0 0.0 -f
../../QMdata/Scan2/butane.dch3_000.fcc 10.0 0.0 0.0 0.0 ; 12 4 3 2
../../QMdata/Scan2/butane.dch3_030.fcc 10.0 0.0 0.0 0.0 -f
../../QMdata/Scan2/butane.dch3_060.fcc 10.0 0.0 0.0 0.0 -f
../../QMdata/Scan2/butane.dch3_090.fcc 10.0 0.0 0.0 0.0 -f
../../QMdata/Scan2/butane.dch3_120.fcc 10.0 0.0 0.0 0.0 -f
../../QMdata/Scan2/butane.dch3_150.fcc 10.0 0.0 0.0 0.0 -f
../../QMdata/Scan2/butane.dch3_180.fcc 10.0 0.0 0.0 0.0 -f
$end
$assign
1 = 2254.5152345 C1-C2
2 = 2130.5117241 C2-C2
3 = 2254.5152345 C1-C2
4 = 3100.0797022 C1-H1
5 = 3100.0797022 C1-H1
6 = 3100.0797022 C1-H1
7 = 2994.0724390 C2-H2
8 = 2994.0724390 C2-H2
9 = 2994.0724390 C2-H2
10 = 2994.0724390 C2-H2
11 = 3100.0797022 C1-H1
12 = 3100.0797022 C1-H1
13 = 3100.0797022 C1-H1
14 = 710.6287265 C1-C2-C2
15 = 346.6702831 C2-C1-H1
16 = 346.6702831 C2-C1-H1
17 = 346.6702831 C2-C1-H1
18 = 401.0089354 C1-C2-H2
19 = 401.0089354 C1-C2-H2
20 = 710.6287265 C1-C2-C2
21 = 404.4980676 C2-C2-H2
22 = 404.4980676 C2-C2-H2
23 = 404.4980676 C2-C2-H2
24 = 404.4980676 C2-C2-H2
25 = 401.0089354 C1-C2-H2
26 = 401.0089354 C1-C2-H2
27 = 346.6702831 C2-C1-H1
28 = 346.6702831 C2-C1-H1
29 = 346.6702831 C2-C1-H1
30 = 313.3362781 H1-C1-H1
31 = 313.3362781 H1-C1-H1
32 = 313.3362781 H1-C1-H1
33 = 346.3056172 H2-C2-H2
34 = 346.3056172 H2-C2-H2
35 = 313.3362781 H1-C1-H1
36 = 313.3362781 H1-C1-H1
37 = 313.3362781 H1-C1-H1
$end
$dependence 1.2
44 = 43*1.d0 ; H1-C1-C2-C2_n=2 = H1-C1-C2-C2_n=2
45 = 43*1.d0 ; H1-C1-C2-C2_n=2 = H1-C1-C2-C2_n=2
46 = 43*1.d0 ; H1-C1-C2-C2_n=2 = H1-C1-C2-C2_n=2
47 = 43*1.d0 ; H1-C1-C2-C2_n=2 = H1-C1-C2-C2_n=2
48 = 43*1.d0 ; H1-C1-C2-C2_n=2 = H1-C1-C2-C2_n=2
$end
$keepff 38 - 48
$scan butane.joyce.delta.dat
1 2 3 4 ; -180. 180. 1.
$end
$scan butane.joyce.dch3.dat
12 4 3 2 ; -180. 180. 1.
$end
Notice how it is possible to introduce dependencies for flexible ICs in the same way as for the stiff ones. In the case of butane, all terminal hydrogen torsions are chemically equivalent. This way, the result of a single QM scan can be used to fit all 6 dihedral functions.
Finally, including both keywords $gracefreq & $gracetors will generate data useful to validate the accuracy of the force field, corresponding to the resulting Hessian matrix and rigid torsions accuracy. You can observe the result of the fitting below.
You may download the final result here.
Continue with the tutorial here to validate the force field.