Use the PM6-DH2 procedure of: M. Korth, M. Pitonák, J. Rezác, and P. Hobza, "A Transferable H-bonding Correction For Semiempirical Quantum-Chemical Methods", J. Chem. Theory Comp. 2010, 6, 344–352 and J. Rezác, J. Fanfrlik, D. Salahub and P. Hobza, " Semiempirical Quantum Chemical PM6 Method Augmented by Dispersion and H-Bonding Correction Terms Reliably Describes Various Types of Noncovalent Complexes " J. Chem. Theory Comp. 2009, 5, 1749-1760. See: Abstract.
The PM6-DH2 method was parameterized to reproduce interaction energies for geometries obtained from high-level quantum mechanical calculations. While it is possible to use it for geometry optimization, please be aware of the following limitation and that the method should be used with extra care. Users are recommended to optimize the geometry with the PM6 or PM6-D methods and calculate the final energy or interaction energy using PM6-DH2.
KNOWN LIMITATION OF PM6-DH2: The gradients obtained by the current implementation of the -H correction do not include the term containing the derivative of atomic charges with change of the coordinates. The structure obtained by minimization using this gradient is not the exact minimum of PM6-DH2 energy. This error is negligible in the case of weaker H-bonds, but in case of a strongly bound formic acid dimer, the error is 0.30 kcal/mol. This error can be avoided by using NOANCI, but then the calculations will then take much longer.
The PM6-DH2 procedure corrects binding errors in the PM6 method. It can be used with geometry optimization or with a single point (1SCF) calculation. Normally, two or three calculations would be needed to get the binding energy.
(A) Optimize the geometry of the two species (call these R1
and R2) separately using PM6, i.e., not using PM6-DH2.
(B) Optimize the geometry of the adduct R1 non-covalently bound to R2,
using PM6.
(C) Calculate the heats of formation (DHf)
of each of the three systems, R1, R2,
and R1 - R2, using 1SCF and PM6-DH2.
The binding or intermolecular interaction energy is then given by DHf(R1 - R2) - DHf(R1) - DHf(R2) Ignore the quantities: "PM6-DH2 Dispersion contribution" and "PM6-DH2 H-bond contribution," these are only of use if you want to know the component contributions.
The rational for this strategy is as follows: Geometry optimization using PM6-DH2 has small errors. These errors cannot easily be corrected (they can be corrected by using NOANCI, but the jobs would then take a long time), however by using PM6 a well-defined stationary point can be easily obtained. By optimizing the various structures separately, the energy penalty that would result from the geometry change when the species form the noncovalent complex is avoided. The binding energy is thus the energy released when the two species form the complex in the gas phase; it is not the energy of interaction of R1 with R2.
Frozen geometries were used during the development of PM6-DH2. By contrast, in the proposed strategy, fully optimized geometries are used. No inconsistency is involved - by sketching a simple Born cycle, it becomes apparent that any errors arising from optimizing the PM6-DH2 parameters using frozen geometries and using those same parameters when calculating the binding energy using the above strategy would be very small; they would contribute only second order perturbation effects, and would be completely negligible.
Without the optional arguments, PM6-DH2 will apply to all atoms. With optional arguments, PM6-DH2 will apply only to the interaction of one part of the system with a second part. The optional argument defines atoms used in one of the parts. Atom numbers are specified by numbers separated by commas, and ranges of the type 645-670 or 645:670. To specify atoms 600, 610, 611, 612, and 630 the keyword would be PM6-DH2(600,610:612,630) or PM6-DH2(600,630,610-612)
If only the dispersion contribution is wanted, replace PM6-DH2 with PM6-D2. If only the hydrogen-bonding term is wanted, replace PM6-DH2 with PM6-H2.
Consider the dimer of formic acid: The binding energy would be the DHf of the formic acid dimer minus the DHf of two well-separated formic acid molecules. Here, two calculations are needed. For this work, the accurate geometries from the S22 database are used.
First calculation:
PM6-DH2(1:5) 1scf C -1.888896000 1 -0.179692000 1 0.000000000 1 O -1.493280000 1 1.073689000 1 0.000000000 1 O -1.170435000 1 -1.166590000 1 0.000000000 1 H -2.979488000 1 -0.258829000 1 0.000000000 1 H -0.498833000 1 1.107195000 1 0.000000000 1 C 1.888896000 1 0.179692000 1 0.000000000 1 O 1.493280000 1 -1.073689000 1 0.000000000 1 O 1.170435000 1 1.166590000 1 0.000000000 1 H 2.979488000 1 0.258829000 1 0.000000000 1 H 0.498833000 1 -1.107195000 1 0.000000000 1
Second calculation:
PM6-DH2(1:5) 1scf C -1.888896000 1 -0.179692000 1 99.000000000 1 O -1.493280000 1 1.073689000 1 99.000000000 1 O -1.170435000 1 -1.166590000 1 99.000000000 1 H -2.979488000 1 -0.258829000 1 99.000000000 1 H -0.498833000 1 1.107195000 1 99.000000000 1 C 1.888896000 1 0.179692000 1 0.000000000 1 O 1.493280000 1 -1.073689000 1 0.000000000 1 O 1.170435000 1 1.166590000 1 0.000000000 1 H 2.979488000 1 0.258829000 1 0.000000000 1 H 0.498833000 1 -1.107195000 1 0.000000000 1
Consider a substrate docked in an enzyme. The binding energy would be the DHf of the docked assembly minus the DHf of the isolated substrate minus the DHf of the isolated enzyme. As in vivo enzymes exist in aqueous media, solvation will be necessary.