JOYCE Tutorial. Part (IV / IV)

Author: Anna Piras, Pablo M. Martinez

Keywords: Tutorial, FF, JOYCE

The final step of the parametrization is the validation. This is a project-dependent step so any relevant and/or convenient test can be performed. It depends also on the system, systems presenting flexible domains would need more validation tests. We propose different tests that will ensure the correct mapping of the QM potential energy surface into the selected internal coordinates. These tests concern:

i) The absolute minima of the PES.

ii) The local harmonic approximation (LHA).

iii) The validity of the FIRA.

iv) The energy fluctuations through a sample MD run.

These are presented in ascending order of complexity. It is relatively easy to pass i), since the equilibrium values are withdrawn from a converged first principles calculation. A good agreement between QM and MM descriptions for iv) is less straighforward, as every degree of freedom is tested simulatenously.

We will make use of the validation toolset available (download from here).

Preparations

Follow the next steps:

Untar the toolset folder: tar xvf joycetools-v0.1.tar.gz
Get into the folder and place the binaries into Bin/ with: make precompiled
Update you PATH to include Bin/ folder (assuming bash): export PATH=$PWD/Bin:$PATH
Take the Step0, Step1, Step2 folders and move their contents to the folder Tests/JOYCE.
Similarly, move the QM results to the Tests/Gaussian and Tests/QMdata to replace the blueprint molecule with your own data.

Test 1: Ground state geometry

First, we may compare the QM minimized structure with the QMD-FF. We use a modified version of GROMACS interfaced with Gaussian to make use of their minimization algorithms. Go to the folder Tests/Gaussian/. Then, remove all the Scan folders and change the Makefile accordingly, i.e. remove all references to scans. Then, run make fchk

You will find the QM and MM minimized geometries in FCHK/joyce.g0_fix.xyz and FCHK/joyce.g0_min.xyz respectively.

We observe that the absolute minimum of the molecule is well captured by our force field.

Test 2: Rigid degrees of freedom

Secondly, we may extend the validation beyond the static picture to assess the goodness of fit within the local harmonic approximation (LHA), i.e. small oscillation limit. This is conveniently given does it for us at the end of the .out file. in Step 1. Remember the $gracefreq keyword produces the plot for us in the standard input file.

         ================================================= 
              Compare Norm Modes from QM and FF 
         ================================================= 

        iFF Freq/FF   ig09 freq/g09      overl    err 
        35   3107.8    36   3096.1       0.989    11.8 
        36   3108.0    35   3095.3       0.988    12.7 
        34   3107.0    34   3090.5       0.916    16.5 
        33   3106.8    33   3087.0       0.974    19.7 
        31   3047.5    32   3046.8       0.917     0.7 
        32   3053.9    31   3024.2       0.974    29.7 
        30   3000.2    30   3021.4       0.969   -21.2 
        29   2999.6    29   3019.9       0.950   -20.3 
        27   2985.2    28   3008.7       0.947   -23.6 
        28   2992.3    27   3002.1       0.974    -9.8 
        25   1491.7    26   1484.8       0.653     6.8 
        18   1414.8    25   1478.3       0.715   -63.5 
        21   1417.8    24   1470.6       0.979   -52.7 
        22   1418.3    23   1468.7       0.992   -50.4 
        19   1414.8    22   1463.0       0.864   -48.2 
        24   1480.7    21   1460.0       0.842    20.7 
        23   1449.1    20   1399.0       0.750    50.1 
        26   1561.8    19   1398.0       0.773   163.8 
        20   1415.4    18   1387.3       0.539    28.1 
        15   1224.6    17   1321.0       0.978   -96.4 
        17   1355.0    16   1310.6       0.938    44.5 
        16   1232.2    15   1277.6       0.883   -45.4 
        13   1029.6    14   1197.0       0.988  -167.4 
        14   1065.1    13   1163.2       0.973   -98.1 
        10    907.1    12   1074.4       0.940  -167.3 
         9    874.3    11   1024.5       0.799  -150.2 
        11    951.1    10    976.5       0.793   -25.3 
        12    952.0     9    955.3       0.902    -3.3 
         8    812.2     8    844.9       0.937   -32.7 
         7    804.4     7    809.4       0.979    -5.0 
         6    689.5     6    741.9       0.989   -52.4 
         5    405.1     5    424.8       0.982   -19.6 
         3    258.7     4    267.9       0.998    -9.2 
         4    318.2     3    257.5       1.000    60.7 
         2    230.3     2    228.1       0.991     2.1 
         1    125.7     1    123.0       0.995     2.7 
                  Standard deviation (cm-1)       65.80

Below, we plot the results taken from the JOYCE output for both eigenvectors and eigenvalues, obtained with the $gracefreq keyword.

Note

Pairing the eigenmodes found in both descriptions is not straightforward and does not have a unique answer. JOYCE starts with the highest frequency eigenmode at the MM level and finds the largest projection with respect to all the QM modes. The best match pairs those two modes. The process repeats itself until only two modes are left, which complete the pairing. Once all the modes are paired, we may look at the eigenfrequencies and check how well they correlate. Note that this algorithm does not ensure maximizing the overall overlap. For example, consider two modes match up to 85%. Perhaps, later on, one mode would have reached a 95% with the previous one, but as they where removed from the candidate pool, the better pairing is not performed. In this manner, this test yields "pessimistic results".

Test 3: Validity of FIRA

The felxible degrees of freedom were obtained through the FIRA. We may relax this assumption and perform a relaxed scan to verify its accuracy. The JOYCEtoolset is the best way for this. The easiest way is to perform a:

make test

in the home directory. The result will be inside each of the test/MMscan/Scan? folders. This tool makes use of GROMACS performing a full relaxation around each of the flexible torsions. The dihed.ndx should be changed with the atomic indexes of each scan so that the procedure is performed sequentially. Below, find the relaxed scans in blue lines compared with the reference QM data in the form of red points. Despite the assumption of FIRA during the procedure, the fit is still very accurate for a relexed scan and therefore valid in the simulations.

Test 4: Virial test

Finally, we test the accuracy of the force field within the statistical physics realm. Naturally, this is the most exhaustive test of all as every comformation is sampled under the influence of thermal fluctuations. We will perform an in-vacuo MD simulation where we will extract several snapshots whose energy will be evaluated at the QM-DFT level. The purpose of this test is twofold:

Verify that the total average energy complies with statistical mechanics through application of the virial theorem.
Assess how closely fluctuations at the MM level follow the QM level.

For that end, we take an in-vacuo MD simulation with GROMACS and extracting various snapshots through the gmx command:

gmx trjconv -f output.trr -o snap.g96 -s topol.tpr -dt 100 -sep -nzero 2 -center -pbc mol

Then, each snapshot is translated into a Gaussian/ORCA file single point evaluation.
One can compare the fluctuations of the QMD-FF intramolecular energy around the mean value predicted by the virial theorem. This is the most strict test of the four, since all degrees of freedom are excited simulaneously.