qmdffgen

Introduction

qmdffgen uses the output of an QM calculation for a molecule to generate a QMDFF for it. (see Grimme lab homepage for original QMDFF release) This QMDFF is able to simulate the full configurational landscape of this molecule and further to describe bond breakings correctly. A QMDFF needs a predefine set of information provided by the QM program:

• Optimized geometry
• Wiberg-Mayer bond orders
• Hirshfeld charges
• Hessian matrix

Provided the correct settings, all those can be obtained from a single QM calculation (see below). After doing the QM calculation, qmdffgen can be invoked:

qmdffgen.x

The program is interative and asks you for the Software used for the reference calculation (currently available: orca, Gaussian, Turbomole and CP2K), the number of QMDFFs to generate (one to three), and the prefix of the QM input files (name without file ending). Then, the QMDFF is generated in a black-box procedure and directly usable for other calculations!

List of examples:

• 1: Octane with Orca (QMDFF for an octane molecule from Orca QM reference) qmdffgen/octane_orca
• 2: Guaicol with Gaussian (QMDFF for a guaicol molecule from Gaussian QM reference) qmdffgen/guaicol_gaussian
• 3: Copper with CP2K (QMDFF for a hydrogen atom on a copper surface from CP2K QM reference) qmdffgen/copper_cp2k
• 4: Cyclopropane coordinate analysis (Use qmdffgen as a tool to generate and investigate internal coordinates of cyclopropane) qmdffgen/cyclopropane_coords

Example 1: Octane with Orca

Folder: qmdff/octane_orca

The input as well as the output of the QM reference calculation with orca is placed here. Basically, the keywords for the reference calculation with orca are always the same. Open the file `octane.inp*:

! PBE0 D3 def2-SVP opt freq
* xyz 0 1
C         -8.36972        0.54078       -0.36595
C         -6.84957        0.56234       -0.35222
...
H         -0.30022        5.11409        2.54651
H         -0.26113        4.29744        0.97442
*
%output
        Print[ P_Hirshfeld ] 1
end

A geometry optimization of the octane molecule was done with the PBE0 DFT functional (with Grimme D3 dispersion correction) and the cheap def2-SVP basis set. For useful output, a geometry optimization opt followed by a frequency calculation freq need to be done, further, the section

%output
        Print[ P_Hirshfeld ] 1
end

must be added for the calculation of Hirshfeld charges. The qmdffgen program needs the octane.hess and octane.out files as reference.

Since these are present, simply type

qmdffgen.x

and give the respective commands during the interactive input. Then, the QMDFF octane.qmdff will be generated.

During the calculation, the central quality parameters are printed into the terminal:

Zero point vibr. energy comparison (FF/true, kcal)  :     150.432     154.650
Mean absolute deviation (MAD) of frequencies (cm-1) :      54.560

The MAD value in particular is a good indicator of the QMDFF-quality. It is calculated by performing frequency calculations for the QM reference and the QMDFF hessian and comparing them. Good QMDFFs achieve values between 20 and 50, if the value is above 100, something might went wrong during the QM reference calculation or the optimization. If the error cannot be located, please let me know. The QMDFF was written into octane.qmdff. This file is the entry-point for all the other programs included into EVB-QMDFF.

A number of additional output files was produced during the QMDFF optimization. these are:

• octane_qmdff.log: Detailed information about the QMDFF generation (QMDFF internal coordinates, fragments, Levenberg-Marquardt fit of the force constants and QMDFF quality). See the QMDFF paper for more details about the whole procedure.
• octane_molden.out: This file can be opened with the molden molecule viewer in order to visualize the normal modes of the different frequencies of the QMDFF hessian.
• coord_def.inp: List of internal coordinates used to construct the QMDFF (bonds, angles, dihedrals) in the formate also used for several other programs inside QMDFF (irc, rpmd, ...).
• coord_analysis.dat: Values of all internal coordinates, of the Wilson matrix for transformation of cartesian gradients into intenal coordinates and the Wilson matrix derivative for the transformation of cartesian Hessians into internal coordinates. These matrices are used for setup of the DG-EVB coupling term in evbopt.

Example 2: Guaicol with Gaussian

Folder: qmdff/guaicol_gaussian

For Gaussian, geometry optimization and frequency calculation must be done separately, since a combined geopt+freq calculation writes out the Hessian in an incomplete formate.Therefore, the file guaicol_opt.gjf performs the geometry optimization and the file guaicol_freq.gjf performs the calculation of the Hessian.

The QMDFF can again be calculated by generated by typing:

qmdffgen.x

Note, that the prefix of the QM reference must be that of the frequency calculation, i.e., guaicol_freq.

The quality of the generated QMDFF can be seen from the lines:

Zero point vibr. energy comparison (FF/true, kcal)  :      85.684      86.016
Mean absolute deviation (MAD) of frequencies (cm-1) :      26.815

The MAD value of 26 suggests a good description of the low-energy part of the potential energy surface of Guaicol with the generated QMDFF (which was written to guaicol_freq.qmdff. The other output files are the same as explained in the first example.

Example 3: Copper with CP2K

Folder: qmdff/copper_cp2k

QMDFFs for periodic systems can be generated as well, if the reference calculations are done with the free open source program CP2K. This example covers a Cu(100) surface slab (3x5 atoms per layer, 4 layers deep) with a single hydrogen atom adsorbed on a hollow position. The same system is treated in one of the dynamic examples, where a second QMDFF is coupled via a dE-EVB coupling term (example 2).

The CP2K calculation can be started with the given input files copper.inp (main input with keywords etc.) and copper.xyz (cartesian geometry of the optimized slab). Additionally, the files BASIS_MOLOPT and GTH_POTENTIALS are needed for definition of the atomic basis sets and the pseudopotentials describing the regions near the atomic cores. More details concerning running the calculation and interpretations of the given input can be found in the [CP2K Wiki])https://www.cp2k.org/howto).

For our pupose, the output file was already calculated. Since it is quite large, it is compressed in the archive outfile.tar.gz. Decompress it first with the command:

tar -xzvf outfile.tar.gz

which generates the file copper.out. Similar to Gaussian and orca, all needed information (Structure, Hessian matrix and Hirshfeld charges) is contained within that one file. Finally, the QMDFF is generated by

qmdffgen.x

Choose C for CP2K and the prefix copper, the QMDFF should be generated within roughly one minute (file copper.qmdff). After the generattion process, the bunch of output files described in the first example are generated, with exception of the coord_analysis.dat file. Since this file becomes increasingly huge for large systems, it is only generated for structures with up to 20 atoms.

Example 4: Cyclopropane coordinate analysis

Folder: qmdff/cyclopropane_coordinates

The qmdffgen program can also be used to generate internal coordinate sets of molecules and analyze them. During this process, no QMDFF is generated. Therefore, no Hessian matrix is needed as input, significantly reducing the cost of the reference calculation. Only Wiberg-Mayer bond orders must be calculated.

In the example, the geometry of the cyclopropane molecule has been optimized with Gaussian 16. The resulting output can now be used to generate the internal coordinate set of cyclopropane as it would be used for a QMDFF of the molecule. Since, however, no Hessian matrix has been calculated, the force constant fit cannot be done and will be skipped.

Besides the Gaussian output files (cyclopropane.log and cyclopropane.chk) an additional keyfile coordinates.key is needed, containing one single line:

check_coord

This command leads qmdffgen to skip the Hessian fit process but nevertheless to initialize the set of internal coordinates for the molecule. The program needs to be called with the keyfile as command line argument:

qmdffgen.x coordinates.key

Then, the user is asked for the used software, the number and QMDFFs and the prefix of the input, as for an usual QMDFF generation. (Note: the coordinate analysis mode can only be done for one QMDFF input at a time!)

After executing, three output files are written:

1. cyclopropane_qmdff.log: The usual logfile of the QMDFF generation, now with limited information, mainly concerning the number and types of internal coordinates.
2. coord_def.inp: The internal coordinate input file in the format as being used for several other applications within Caracal.
3. coord_analysis.dat: Main output of the analysis, containing a list of all generated internal coordinates with their current values, as well as the elements of the Wilson matrix (derivatives of the internal coordinates with respect to Cartesian coordinates, needed for a coordinate transformation) and the derivatives of the Wilson matrix.