A third-generation dispersion and third-generation hydrogen bonding corrected PM6 method: PM6-D3H+

PM6 implementation in GAMESS

E(PM6) is the molecular PM6 (Stewart, 2007) energy, which is taken as the gas phase energy unless otherwise noted. As part of this work we implemented the PM6 method in the GAMESS program for elements up to neon. The PM6 method also involves d-orbitals for elements past neon but the associated integral code has not yet been implemented. Physical constants in the semi-empirical part of the GAMESS source code were updated, to match those in the current version of MOPAC. All semi-empirical methods in GAMESS uses a finite difference scheme for gradient evaluation.

Dispersion correction E(D3)

E(D3) is the third generation dispersion correction developed by Grimme et al. (2010) DFT-D3 and implemented in GAMESS by R Peverati. Unless otherwise noted E(D3) refers to the pair-wise additive dispersion correction as proposed in Grimme et al. (2010). Only the zero-damping version was used, with dispersion order 6 and 8. The fitting parameters are those obtained by Grimme (2012) for PM6. As described by Grimme, the parameter s₆ is set to unity, α was set to its default value. s₈ and the scaling parameter s_r,6 of the atomic cut-off radii used in the dispersion damping function are fitted parameters as in standard DFT-D3 (see Table 1 for parameters). Thus only s₈ and s_r,6 are optimized by Grimme for PM6-D3H, which is also used for PM6-D3H+. The gradient of the dispersion correction is evaluated numerically, by using a centered finite difference scheme, for three-body calculations, and analytically for two-body calculations.

Hydrogen bond correction E(H+)

E(H+) is a slightly modified version of the third-generation hydrogen bonding correction, H+, by Korth (2010), which is given by: (2)${\displaystyle E\left(\mathrm{H+}\right)=\sum _{AB}\frac{C_{\mathrm{A}}+C_{\mathrm{B}}}{2r_{AB}^2}\cdot f_{\mathrm{geom}}\cdot f_{\mathrm{bond}}\cdot f_{\mathrm{damp}}}$ where the sum runs over all hydrogen bonds involving N and O atoms. r_AB is the donor–acceptor distance for the given hydrogen bond geometry, with A and B being the two possible acceptor/donor electronegative atoms, either oxygen or nitrogen. C_A and C_B are adjustable parameters and refer to either C_N and C_O. C_N and C_O are re-parametrized as part of this work as described below.

The geometrical correction f_geom is defined as (3)${\displaystyle f_{\mathrm{geom}}={cos}^2\theta \cdot {cos}^2{\varphi }_{\mathrm{A}}\cdot {cos}^2{\psi }_{\mathrm{A}}\cdot {cos}^2{\varphi }_{\mathrm{B}}\cdot {cos}^2{\psi }_{\mathrm{B}}}$ where θ is the angle defined by atom A, atom B and the hydrogen (see Figs. 1 and 2). The angle ϕ, and torsion angle ψ are both defined by the hydrogen bonding geometry. The angles ϕ are calculated from the difference between the target angle ϕ_target and the present bond angle in the complex Φ_X. The target angle ϕ_target is the optimum angle for hydrogen bonds. Target angles are defined in a complicated heuristic fashion, please see the source code posted on GitHub for more details (JensenGroup, 2014b). The torsion angles ψ are defined similarly and calculated as the difference between target dihedral angle and the structural angle Ψ. Where Ψ_X is the dihedral angle between R₁R₂X…H, which is used for both the donor and acceptor as seen in Figs. 1 and 2. Here R₁ is defined as the R_x closest to the hydrogen.

The bond damping function f_bond is defined as: (4)${\displaystyle f_{\mathrm{bond}}=1-\frac{1}{1+exp\left[-60\cdot \left(r_{\mathrm{XH}}/1.2-1\right)\right]}}$ where r_XH is the distance between the hydrogen atom and the donor atom, which is defined as the shorter one of the distances r_AH and r_BH. The damping function f_damp is defined as: (5)${\displaystyle f_{\mathrm{damp}}=\left(\frac{1}{1+exp\left[-100\cdot \left(r_{\mathrm{AB}}/2.4-1\right)\right]}\right)\left(1-\frac{1}{1+exp\left[-10\cdot \left(r_{\mathrm{AB}}/7.0-1\right)\right]}\right)}$ where r_AB is the distance between the two electronegative atoms A and B.

The E(H+) implementation differs slightly from the one originally proposed by Korth (2010). Changes were made to avoid problems with optimization of hydrogen bond complexes involving particular configurations, including especially ketone (C=O) groups interacting with amide-like (NR3) groups. In the original implementation, optimization problems can originate from target angle calculation based on the torsion angle of the NR3 group. Target angles are the optimal (text-book) angles for a given H-bond arrangement. H-bond energies are computed based on the deviation of all angular coordinates from their respective target (optimal) angles, see Korth (2011) for a detailed explanation. The target angle would switch during optimization steps as the definition of the torsion angle would switch, and never find a minimum, as the torsion angle is defined as seen in Fig. 1. The model was updated with new target angles for tetragonal NR3 configuration case, and the estimation of target angles for NR3 groups now based on the hydrogen bonding configuration (with a double bond indicating a planar structure).

The analytical gradient is done using internal coordinates from the energy model (angles and distances in Eq. (2)), and an algorithm for converting the gradient to the Cartesian atomic coordinates.

Source code for the H+ module, including gradient code, is available on GitHub (JensenGroup, 2014b). The PM6-D3H+ will be made available in the official version of GAMESS as soon as possible.

Parameterization of correction terms

Because we use a different dispersion energy function than in the previous DH+ model and make modification to the original hydrogen bonding correction model, it is necessary to determine new optimum values for the C_N and C_O. The parameters for H+ are parameterized to minimize the root-mean-square deviation (RMSD) between the interaction energies for PM6 with dispersion correction only (PM6-D3) for a subset of structures from the S22 and S66 data sets (1–7 and 1–23, respectively), plus the H+ term and the estimated CCSD(T)/CBS reference interaction energy. The C_N and C_O parameters are then scanned in ranges from −0.2 to 0.0, around the original optimum. A global optimum was found at C_N = −0.11 and C_O = −0.12, with a RMSD of 1.11 kcal/mol, as seen in Fig. 3 and Table 2.

This was done using both two and three-body dispersion, but including three-body dispersion did not make any substantial difference in the resulting optimum, and the default was set to two-body for PM6-D3H+, because of the extra computational time associated with three-body. The computational cost becomes a time consuming issue for protein-sized molecules. The final set of parameters for both dispersion and hydrogen bond correction terms can be seen in Table 1.

Interaction energies

Table 2 shows results of PM6, PM6-DH+ and PM6-D3H+ for the full, dispersion and hydrogen bond dominant complexes sets of the S22 and S66 from BEGDB. Root-mean-square deviation (RMSD), mean absolute deviation (MAD) and maximum error span (Max) with respect to the benchmark estimated CCSD(T)/CBS interaction energies are given in kcal/mol. The PM6-D3H+ method was tested using both two and three-body dispersion.

Overall, the accuracy of PM6-D3H+ is very similar to PM6-DH2 and PM6-DH+, with RMSD and MAD values within 0.02 kcal/mol of one another. The main difference is that the maximum error for PM6-D3H+ is 1.42 and 0.36 kcal/mol smaller than for PM6-DH2 and PM6-DH+, respectively. The maximum error for PM6-DH2, -DH+, and -D3H+ were observed for the S66-19, S66-60, and S66-65 dimer, respectively and, in general, we did not notice any particular dimer that resulted in unusually large errors for all three corrections. All interaction energies can be found in Supplemental Information. The differences in RMSD and MAD between methods are slightly larger (up to 0.13 kcal/mol) for subsets where dispersion and hydrogen-bonding dominate. Including three-body dispersion correction had no substantial effect on accuracy, but might play a role for large systems.

Next, we test PM6-DH3+ on two sets of molecules not in the training set. Table 3 lists computed interaction energies for formamide dimer, pentamer–monomer, and trimer–trimer (Fig. 4) computed with various methods. Compared to MP2/TZVP PM6-DH2 performs best for this particular system, while PM6-DH+ and PM6-D3H+ appear to perform roughly similarly, with mean absolute deviations (MAD) of 0.8 and 1.3 kcal/mol, respectively. However, it is interesting to note that the decrease in interaction energy on going from the dimer to the pentamer–monomer predicted by PM6-DH+ (3.6 kcal/mol) is somewhat lower than that predicted by other methods corrections and MP2/TZV (4.1–4.6 kcal/mol). This decrease comes primarily from cooperative polarization effects that are accounted for by the underlying PM6 method, and PM6, PM6-DH2, and PM6-D3H+ all predict similar decreases. It is not clear why the DH+ terms leads to an underestimation of the cooperative effect.

Table 4 contains RMSD, MAD, mean-deviation (MD) and maximum deviation relative to estimated CCSD(T)/CBS//MP2/pVTZ interaction energies computed for 12 hydrogen bonded base pair complexes (Table S1) from the JSCH-2005 (Jurečka et al., 2006) set from BEGDB. The 12 complexes represent all the complexes in the JSCH-2005 set with hydrogen bonds involving N and O atoms and for which interaction energies have been computed at a level similar to that used in the parameterization of PM6-D3H+ [i.e., CCSD(T)/CBS// MP2/pVTZ]. For this set all three corrections offer very significant increases in accuracy (e.g., a ca 8 kcal/mol decrease in the MAD) compared to PM6. As for the training set (Table 2) the accuracy of PM6-DH2, PM6-DH+, and PM6-D3H+ are very similar, with MADs between 0.7 and 1.1 kcal/mol.

Geometry optimization

All structures from the S22 and S66 data sets were optimized with PM6, and PM6-DH+ using MOPAC or PM6 and PM6-D3H+ using GAMESS to test how well the methods reproduce the reference MP2/cc-pVTZ geometries and to compare the optimization algorithms in GAMESS and MOPAC.

For the GAMESS optimizations we used the default (quasi Newton–Raphson) geometry optimizer and defined convergence as having a maximum gradient component less than 5 × 10⁻⁴ Hartree/Bohr and an RMS gradient less than 5/3 × 10⁻⁴ Hartree/Bohr. These convergence criteria are five times higher than the default and are chosen because we have found that for large systems these criteria can lead to significantly faster convergence without affecting the structure or final energy significantly. See supporting information for GAMESS examples of input files. For complex 58 in the S66 set it was necessary to re-compute the Hessian every 20 steps to obtain convergence and in the case of complex 22, 51, and 58 it was necessary to skip the projection of translational and rotational degrees of freedom from the gradient to obtain convergence, which was done by settings the keyword PROJCT=.F. in $Force. For 11 of the complexes (see Table S1) it was necessary to decrease convergence criterion to 10⁻⁴ Hartree/Bohr in order to remove imaginary frequencies. In the case of complex 4 and 5 from S22 PM6-D3H+ predicted that the minimum has C₁ symmetry rather than C_s as predicted by MP2, and a deviation in the planarity structure of 0.1 Å was needed (added to the first atom). This is not the case for PM6 and thus a result of the D3H+ energy correction.

For the MOPAC optimization we used the LBFGS geometry optimizer because we found that this is the only optimization algorithm that can be practically applied to optimization of large systems. Using eigenvector following leads to termination of the geometry optimization and the following error message: “TRUST RADIUS NOW LESS THAN 0.00010 OPTIMIZATION TERMINATING”. Based on the output the convergence criterion for the LFBGS optimizer appears to be a change heat of formation of less than ca 0.1 kcal/mol during several consecutive optimization steps. For PM6, this convergence test was not passed after ca 200 geometry optimization steps for complex 10 and 17 from the S22 set and 29, 53, and 54 from the S66 set. For PM6-DH+, this convergence failed after ca 140 geometry optimization steps for complex 11 from the S22 set and 53, 54 and 60 from the S66 set. In all these cases MOPAC terminates the geometry optimization after the mentioned number of steps with the message: “A FAILURE HAS OCCURRED”.

The results are summarized in Table 5. The average RMSD between the MP2/cc-pVTZ and semi-empirical structures are below 0.28 Å for all methods and a factor of two lower for the GAMESS optimizations. The RMSD was calculated using the Kabsch algorithm (Kabsch, 1976), for all the atoms, including hydrogens. For the hydrogen bonding subset RMSD was calculated for the hydrogen bond lengths, which are much lower with GAMESS, and with PM6-D3H+ being the lowest with a RMSD of 0.08 Å. The GAMESS optimizations converge, on average, in 30 steps, while the MOPAC optimization takes 10 times more steps.

Furthermore, MOPAC optimized geometries tend to have a significantly larger RMS gradient, compared to GAMESS. This leads to significantly more imaginary frequencies in a subsequent vibrational analyses compared to those obtained with GAMESS. In the case of MOPAC 54, 17, 15, and 2 geometries result in 0, 1, 2, and ≥3 imaginary frequencies, while the corresponding numbers for GAMESS are 82, 6, 0, and 0 (Tables S2–S4). Using the (default) eigenvector following algorithm in MOPAC for comparison results in 60, 19, 5, and 4 geometries with 0, 1, 2, and ≥3 imaginary frequencies, respectively, with complexes 1 and 3 from S22 and 1 and 20 from S66 failing the optimization.

For four of the six cases where a GAMESS optimization leads to a structure with a single imaginary frequency a convergence criterion of 10⁻⁴ Hartree/Bohr is used, but lowering the convergence criterion further does not remove the imaginary frequencies. In the sixth case, complex 16 in the S66 set (water hydrogen bonded to an amide group— Fig. 5), the optimization stalls, when setting convergence criterion to 10⁻⁴ Hartree/Bohr, with the maximum gradient oscillating between 3 × 10⁻⁴ and 2 × 10⁻⁴ Hartree/Bohr. This is due to the dihedral angle ψ (Eq. (3)) which is defined as R₁R₂X…H (cf. Figs. 1 and 2), where R₁ is defined as the atom closest to the H atom. In the case of the amide–water hydrogen bond, R₁ and R₂ are the two water H atoms, which are approximately equidistant from the amide proton. The oscillation in the maximum gradient is caused by the oscillation between two different definitions of ψ, which has an effect on the gradient direction. The normal mode associated with the imaginary frequency for the structure converged with a convergence criterion of 5 × 10⁻⁴ corresponds to a motion between these two structures, so this is likely the explanation for the imaginary frequency. Similarly, in the case of the complex 1 in the S22 set (ammonia dimer), we believe the imaginary frequency is due to highly symmetric hydrogen configuration, with switching torsion angles (atomic definition of ψ). Since this only affects structures with highly symmetric hydrogen bonds it is unlikely to cause problems in most applications. We note that the PM6-DH+ method has the same problem.

In the remaining four cases where a GAMESS optimization leads to a structure with an imaginary frequency the cause is most likely an extremely flat potential energy surface for the corresponding degrees of freedom: all imaginary frequencies are <31i cm⁻¹. Similarly, the lowest real frequencies for these five cases are all <40 cm⁻¹.

In summary, the PM6-D3H+ method as implemented in GAMESS offers an attractive alternative to PM6-DH+ in MOPAC in cases where the default geometry optimizer fails to find a converged structure and the LBFGS optimizer must be used and a vibrational analysis is needed e.g., when computing vibrational free energies.

Application to protein structure refinement

In this section we test the applicability of the PM6-D3H+ method, combined with the PCM for bulk solvation as implemented in GAMESS, to geometry optimization of large systems such as proteins and compare to corresponding calculations performed using MOPAC.

We optimize the proteins Chignolin (1UAO) and Trp-Cage (1L2Y), which are two small proteins with 138 and 304 atoms, respectively. We optimize the structures using PM6-DH+ in MOPAC (Stewart, 2012), and PM6-D3H+ in GAMESS (Schmidt et al., 1993), with and without implicit solvent models. The optimized semi-empirical structures are compared to the reference structure optimized at the RHF/6-31G(d) level of theory using dispersion correction (DFTD3) and two-body Fragment Molecular Method (FMO2). Previous calculations by Nagata et al. (2012) have shown that this level of theory yields protein structures in good agreement with corresponding MP2 calculations. Optimized reference structures are available on GitHub (JensenGroup, 2014a).

The results are summarized in Table 6. The RMSD values are about 1 Å in the gas phase for both methods, with PM6-DH+ being slightly smaller. The RMSD values for the structures in solution are slightly larger compared to the corresponding gas phase values for PM6-DH+, and slightly smaller for PM6-D3H+. The structural overlap between PM6-D3H+/PCM optimization and the reference structure can been seen in Figs. 6 and 7. For Trp-cage both methods converge in about half the number of steps compared to gas phase. MOPAC requires significantly more optimization steps than GAMESS to converge, but the overall time for optimization of the structures is by far faster than GAMESS. The difference in CPU time per geometry optimization step is significantly larger for optimization in bulk solvent, which indicates that it is the difference in the COSMO and PCM interfaces that differ most in terms of CPU requirements. Despite being significantly slower than PM6-DH+/COSMO, the PM6-D3H+/PCM implementation in GAMESS is sufficiently fast to make geometry optimizations of small proteins feasible.

The number of imaginary frequencies computed for the optimized protein geometries (N_i) are listed in Table 6. Again, the GAMESS optimization leads to significantly fewer imaginary frequencies: 3 and 3 using PM6-D3H+/PCM implemented in GAMESS, compared to 5 and 12 for Chignolin and Trp-cage using PM6-DH+/COSMO implemented in MOPAC. In the case of GAMESS the number of imaginary frequencies can be reduced to 0 for both proteins by decreasing the geometry optimization criterion (OPTTOL) to 1 × 10⁻⁴ aus. This required 205 and 298 additional optimization steps for Chignolin and Trp-cage, respectively.

The relative speedup from running in parallel in solvent is shown on Fig. 8, where no improvement is observed beyond 8 cores for all methods. The timings were done on 24 core AMD Opteron(tm) Processor 6172 @ 2.1 GHz machine for GAMESS and 8 core Intel(R) Xeon(R) CPU X5560 @ 2.80 GHz for MOPAC, because we were unable to get MOPAC running on the AMD ones. Using the dispersion correction and hydrogen bond correction on the PM6 method in GAMESS reduces the relative speedup from 4 to about 2. The correction terms to the PM6 energy only runs in serial, and a modest speedup could be gained by parallelising them. Here we note that the poor scaling of run times with regards to the number of CPUs used is an inherent problem for semi-empirical since the matrix diagonalization in the SCF procedure cannot be efficiently parallelized (Maia et al., 2012).