SireAdvanced Molecular Simulation Framework

I've made a decision regarding my QM forcefields - I won't add any MM terms to them! I'm back working with MolproFF and I was thinking about how I would use it during a QM/MM simulation. Molpro can handle QM/MM by allowing point charges to be placed around the QM atoms, thus allowing electrostatic polarisation of the wavefunctions when they are calculated. This polarisation thus incorporates the electrostatic interaction energy between the QM and MM atoms. However, it does not include the van der waals interactions between the QM and MM atoms. These are typically represented by adding on the Lennard Jones energy between the QM and MM atoms. I have been debating whether or not this LJ calculation should be part of MolproFF. The obvious argument in favour of including the calculation is that then QM/MM requires only a single forcefield, and that the LJ calculation could occur at the same time as the calculation of which MM atoms are within the cutoff distance of the QM atoms. This is a good and compelling argument. However, the main argument against is that including LJ in MolproFF would make it more difficult to use MolproFF with other methods of calculating the van der waals energy (e.g. more accurate potentials). Also, and this is the most compelling argument in my opinion, including the LJ energy in MolproFF would mean that the LJ energy is calculated in series with the Molpro QM energy. Using a separate LJ forcefield would allow the user to calculate the LJ energy in parallel with the Molpro QM energy. Indeed, the QM/MM van der waals calculation could then be lumped with the MM forcefield calculations and all be performed together while the Molpro QM calculation is performed on another processor.

So, I was very much leaning towards not including the vdw calculation in MolproFF. Then, the clincher was that I thought about how MolproFF would be used during a free energy calculation. In particular, what if I wanted to calculate the free energy of going from a QM/MM representation of the system to a pure MM representation? The only difference between the two representations would be in the forcefields used to represents the internal energy of the QM region and the interaction between the QM and MM regions;

For the QM/MM forcefield;

while for the pure MM forcefield

The forcefield used for the free energy simulation would therefore be;

Note how the LJFF is independent of lambda. If the vdw calculation was part of MolproFF then there would be no way of separating it out, and we would have had to calculate the LJ energy as part of MolproFF and then again in a separate LJFF that was tied to CoulFF (e.g. perhaps in a combined CLJFF). By separating the vdw calculation from the QM calculation, we not only increase flexibility by allowing easy combination with any vdw forcefield, but we also save time by allowing the vdw calculation to be independent of lambda.

Finally, the really nice thing about the above free energy equation is that the energy from MolproFF appears only once, and both it, and the IntraFF+CoulFF are scaled by lambda. This means that once the energies of MolproFF, IntraFF and CoulFF are evaluated, it then becomes possible to calculate the total energy for that configuration for any value of lambda! This means that the code could, with trivial extra cost, record the average of several perturbed states simultaneously, thereby calculating estimates of the free energy change for many lambda values simultaneously. While the errors would increase as delta lambda increased, it would allow initial estimates of the free energy change to be made, and would allow perhaps a more robust error analysis, as multiple estimates of the change in free energy between two lambda values can be made from multiple simulations at different reference states. Also, as none of MolproFF, IntraFF or CoulFF would depend on lambda, we could also use this forcefield in thermodynamic integration, as the gradient of the energy with respect to lambda is just MolproFF - (IntraFF+CoulFF). So, once again, (assuming the same vdw representation in the QM/MM and pure MM system), the vdw energy is not involved in the free energy change, and therefore keeping it out of MolproFF is the best thing to do.