PM6 in GAMESS, Part 2

Okay, so I'm still working on implementing PM6 integrals in GAMESS.

I got the source code from MOPAC 7.1 which includes d-integrals for the MNDO-D method (which is what Jimmy Stewart is using for PM6 in the newest MOPAC (hopefully), which originates from a program written by Walter Thiel).

So the strategy is simply to 'export' the subroutines / modules from MOPAC 7.1 needed to replicate the d-integrals in GAMESS (written in Fortran 90),  and 'import' them into GAMESS-US.
Now, the semi-emperical part of GAMESS-US is actually based on a older version of MOPAC (written in Fortran 77) so the subroutines should be very similar to the code I'll be trying to import.

First part of this mission is too map the relevant subroutines in both GAMESS and MOPAC. And hopefully I'll be able too se a pattern and merge the 'trees'.

The map for MOPAC is:

And the map for GAMESS is:

Now I just need an idea for merging the two trees. Since Stewart based his d-integrals on code he got from Thiel it seems like most of the subroutines is collected in a single file called mndod.F90 (fitting name, lol).

This means I 'just' (I beginning to hate that work) need to copy-and-paste the file into GAMESS and make sure the file is hooked to GAMESS common blocks instead of the fortran 90 modules from MOPAC. So step 1: Include the file and make it compile (which is a lot of rewriting interfaces and modules into the actual file so it is a standalone solution.)

The highlighted area is only the first part of the problem though. After the fock matrix has been put together with the new and cool d-items the matrix needs to be solved and we need the fockd1 and fockd2 for that. They are conveniently also put in the same file with the rest of the subroutines.

Furthermore I have been told by +Jan Jensen that I need to watch out for 'guessmo' subroutine when implementing the new integrals. As described in his figure;

So to recap, implementation in 5 easy steps (said in a very television kitchen accent):

Step 1: Get mndod.f90 compiled with gamess (using gamess common blocks instead of mopac modules and interfaces)
Step 2: Integration: Make IF(PM6) and run the mndod code instead of gamess with pm6 parameters and more.
Step 3: Integration: Make IF(PM6) and run the fock-d 1 and 2 instead in mpcg()
Step 4: Find out why it does not work and solve the problem.
Step 5: Celebration.

To be continued!

PM6 in GAMESS, Part 1

Okay, so I'm working on implementing the semi-empirical method PM6 (by Jimmy "Mopac" Stewart) in GAMESS-US.

The status is; (before I started working) GAMESS has up to and including PM3 already implemented. So the idea is just to update the SE parameters and substitute the subroutines necessary to get PM6 working. Without prior knowledge to GAMESS this really did not sound like a big deal, as the differences between PM6 and PM3 only lies in the way the parameters are used (roughly). The parameterization of PM3 (and AM1) is utilised in the core-core repulsion term (nuclear repulsion) of the Heat of Formation to compensate for the aproximations made in SE methods. Heat of Formation is calcuated acordingly:

$$\Delta H_f = E_{\rm Elect} + E_{\rm Core} - \sum_{A}^{} E_{el}^{A} + \sum_{A}^{} \Delta H_{f}^{A}$$

The emperical parameters from PM3 is fitted via a scaleing factor on the core-core term to fit the experimental heat of formation for the molecule. Fitting of the data and derivation of the parameters was done my Jimmy Stewart, for his program MOPAC where the methods were original implemented. The PM3 core-core repulsion term looks like this;

$$E_n(A,B) = Z_A Z_B \langle s_A s_A | s_B s_B \rangle \left ( 1 + e^{-\alpha_A R_{AB}} + e^{-\alpha_B R_{AB}} \right ) $$

which is then summed over all nuclear repulsions/interactions between any atom A and B. This core-core term needs to be substituted with the new term from PM6:

$$E_n(A,B) = Z_A Z_B \langle s_A s_A | s_B s_B \rangle \left ( 1 + x_{AB} e^{-\alpha_{AB} (R_{AB} + 3 \cdot 10^{-4} R_{AB}^6)} \right )$$

Note that the $\alpha$ parameter is now a di-atomic parameter unlike the mono-atomic parameter is PM3. Another parameter $x$ is also introduced, but that is 'pritty much it'. (there are also a Lennard Jones Term and a van der Waals term, but that is for another blog post). The parameters are all located in the PM6 article, but Jimmy Stewart was kind of enough to send his files including his implementation of PM6 core equation and the list of all parameters. This saved me alot of pointless typing time, so thanks!

Okay, so after implementing the new PM6 specific code and the corresponding parameters, I discovered that the result did not match its MOPAC equivalent. In fact nether electronic, core or the total energy fit the MOPAC value. This was happening for single point energy calculation even for very small molecules (even water). For reference I did a similar test for PM3 and AM1, and found that the already implemented methods results did not fit it's rightfull energies (MOPAC energy) with the same order of magnitude as PM6, which was quickly discovered to be size dependt. This is clearly shown in the below figure, which shows single point energy calculations on a carbon chain from 1 carbon to 20 carbons
.

Energy difference $\Delta E$ is calculated from mopac energy minus the corresponding gamess energy.

Arrgghh! How am I going to implement a new method, when the already implemented methods varies this much from the original program?

Okay, so the problem was that the SE part of GAMESS was based on a very old version of MOPAC, and so we figured alot of the energy deviation must be originating from the lack of update on physical constants. The MOPAC integrals use two physical constants to calculate the integrals in atomic units, namely bohr radius and electron volts, so by using grep I found all the places where the constants/variables were defined (which was alot!), and then updated acording to the constant defined on MOPAC's website using a common block, instead of alot of local instances.

This resulted in

Okay, but is this better? Hell yeah! The total energy is clearly more stable compared to MOPAC energy, which is the energy that matters most. The deviation in the nuclear and electronic energy looks very much linear which hints to more constants needed to be updated. Note I have only updated the constants located in MOPAC part of GAMESS and therefor only effects the semi-empirical part of GAMESS.

However the effect is there, and even though the energy is working now, it will prove a problem for people who wants to reproduce data already calculated with GAMESS. So be warned GAMESS users, keep a copy of your GAMESS when the PM6 update is integrated in GAMESS-US.

PM6 Gradient:

The integration of gradient was actually really easy, because GAMESS only uses numerical gradients for semi-emperical calculations.

Am I done? Unfortunately no. To get PM6 fully working I need to implement the d-integrals from MOPAC. As it is now only s- and p-orbitals are used for calculating the integrals. Is that easy? No.

To be continued...

tldr; PM6, PM3 and AM1 did not work as expected, which was partially fixed by updating physical constants in the semi-empirical part of GAMESS. PM6 energy and gradient now works up to including Ne, but will need d-orbitals before it is fully operational.