Thursday, November 17, 2011

"Missing" MOPAC parameters

A long, long time ago in a land far, far away I implemented MNDO, AM1, and PM3 in GAMESS.  This was done by taking chunks of code from MOPAC that contained contained the parameters, integral code, and Fock matrix builder.  In doing so, I never noticed that the parameter file contained more parameters than where published in the papers describing the method.  This conundrum came back to haunt us recently as we're trying to implement PM6.

Of course, it's not a conundrum at all: the "missing parameters" are "simply" functions of the other parameters, so they are not missing but finding these functions is not so "simple" either, and I couldn't have done it without some very helpful emails from Jimmy Stewart.

To find expressions for two of the parameters, called $DD$ and $QQ$ in GAMESS, you have to dig out the MNDO paper from 1976 where they are given in equations (11) and (12)

$DD = \frac{5}{\sqrt 3 }\frac{{\left( {4 \cdot ZS \cdot ZP} \right)^{5/2} }}{{\left( {ZS + ZP} \right)^6 }}$   (1)

$QQ = \sqrt {\frac{3}{2}} \cdot \frac{1}{ZP} $   (2)

The $AM$ parameter turns out to be related to $\rho_0$ in the MNDO paper, and is just the $GSS$ parameter in Hartree units (1 Hartree = 27.21 eV).

$AM=\frac{GSS}{27.21}$   (3)

Similarly, AD and AQ are then related to $\rho_1$ and $\rho_2$ which "are calculated by numerical methods$^{20}$", where reference 20 is a paper on NDDO integrals.  Here, a few sentences on page 95 provided sufficient clues:

$AD=\frac{1}{2 \rho_1}$   (4)

where $\rho_1$ is the solution to
$\frac{1}{2} \left( 4\rho_1^2 \right)^{-1/2}-\frac{1}{2} \left( 4DD^2+4\rho_1^2 \right)^{-1/2}=\frac{HSP}{27.21}$   (5)

which is equation (56) where $R = 0$ and $D_1^A=D_1^B$.  Back then I assume they wrote some Fortran code to solve the equation iteratively, but now this can be done in a few seconds using a web browser.

Similarly, for $AQ$:

$AQ=\frac{1}{2 \rho_2}$   (6)

$\frac{1}{4} \left( 8QQ^2+4\rho_2^2 \right)^{-1/2}-\frac{1}{2} \left( 4QQ^2+4\rho_2^2 \right)^{-1/2}+\frac{1}{4} \left( 4\rho_2^2 \right)^{-1/2}= \frac{HPP}{27.21}$   (7)

where the latter equation is derived from equation (62) in the NDDO integral paper (more on $HPP$ below).

This leaves us with the last parameter, $EISOL$, which corresponds to $E_{el}^A$ in the MNDO paper: which "are calculated from restricted single-determinantal wave functions using the same approximations and parameters as in molecular NDDO calculations."  Not much to go on if you ask me, inspection of the parameters themselves held some clues.  For example for hydrogen $EISOL=USS$

So, after a bit of fiddling around I was able to reproduce $EISOL$ for beryllium $\left( [He]2s^2 \right)$ (Remember that semiempirical methods ignore the core electrons):

$\begin{aligned}EISOL&=2h_{11}+J_{11}+4J_{12}-K_{12}+2h_{22}+J_{22}\\&=2h_{22}+J_{22}\\&=2USS+GSS\end{aligned}$   (8)
Carbon $\left( [He]2s^22p_x^12p_y^1 \right)$ was a bit trickier to verify, mainly because I still didn't know what $HPP$ was:


Luckily, Jimmy Stewart knew this from the top of his head: $HPP = (GPP-GP2)/2$

Update: I just found some useful pages here and here

Update: For some elements $GPP < GP2$.  In that case use $HPP = min(0.1,HPP)$  (Thanks again, Jimmy!)

Update: Turns out that you can get these parameters printed out in MOPAC using the HCORE keyword:

$DD = DD2$

$QQ = \sqrt{2}DD3$

$AM = \frac{1}{2 PO1}$

$AD = \frac{1}{2 PO2}$

$AQ = \frac{1}{2 PO3}$

(This post was written using MathJax, which is very easy to install on Blogger)


Jan Jensen said...

And yes, you can also use MathJax in the comments:

Geoff Hutchison said...

Do you think you'll implement Jimmy's "expanded" AM1 and PM3? I seem to remember in an ACS meeting (and possibly in the PM6 paper) he extended the previous methods to a few more elements. This served as a comparison between the older methods and PM6 -- which was obviously better.

Jan Jensen said...

No, no plans there. We're implementing PM6 because the resulting TS structures agree much better with QM calculations than AM1 and PM3.

The first milestone will be elements up to Ne, as the rest require d-orbitals, which are not in GAMESS yet.