JOYCE Tutorial. Part (III / III)
Author: Anna Piras, Pablo M. Martinez
Contact: anna.piras@pi.iccom.cnr.it
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 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 iii) 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 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. This is the most fundamental test we can perform to validate our force field and internal coordinate selection.
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. In Part II/III, the QM Hessian was used to fit the rigid degrees of freedom. Here, we quantify how good of a representation the QMD-FF offers as compared to the DFT reference. We will compare the predicted eigenvectors and eigenvalues of the Hessian at QM-DFT and QMD-FF level. JOYCE output conveniently does it for us at the end of the .out file. in Step 1.
=================================================
Compare Norm Modes from QM and FF
=================================================
Standard deviation (cm-1) 81.91
iFF Freq/FF ig09 freq/g09 overl err
27 3199.2 27 3202.5 0.725 -3.3
26 3198.4 26 3193.6 0.972 4.9
25 3183.3 25 3177.6 0.748 5.7
24 3149.2 24 3150.5 0.946 -1.3
23 3147.6 23 3148.1 0.971 -0.5
22 1649.0 22 1635.1 0.972 13.9
21 1645.9 21 1627.5 0.817 18.4
19 1515.5 20 1506.8 0.989 8.7
18 1424.7 19 1465.2 0.623 -40.5
17 1261.8 18 1369.2 0.532 -107.4
20 1561.1 17 1305.9 0.531 255.2
15 1120.3 16 1235.4 0.853 -115.1
16 1140.9 15 1158.9 0.964 -18.1
9 934.3 14 1089.1 0.693 -154.8
11 952.0 13 1075.9 0.984 -123.9
7 797.1 12 1047.8 0.593 -250.7
13 1019.7 11 1017.5 0.995 2.2
14 1096.5 10 1008.8 0.592 87.6
12 1009.5 9 1001.4 0.996 8.1
10 948.1 8 956.9 0.990 -8.8
8 888.8 7 896.5 0.997 -7.7
6 784.6 6 762.8 0.806 21.8
5 694.5 5 718.3 0.816 -23.8
4 656.0 4 665.2 0.964 -9.2
3 603.9 3 610.9 0.974 -7.0
2 438.2 2 419.4 0.999 18.8
1 409.5 1 385.3 0.998 24.2
Standard deviation (cm-1) 87.08
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: 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.