A basic implementation of virtual sites in SOMD/OpenMM

[bieniekmateusz]

Hi

I created a branch vsites2 and uploaded the changes contained in one commit here:
https://github.com/bieniekmateusz/Sire/tree/vsites2

There are some useful code snippets I left in the code that I'll have to remove before making a PR.

It doesn't depart from Sofia's code, and is only implemented for the "discharge" stage. The virtual sites are simply not assembled in the vanish stage.

Many thanks, Mat

[bieniekmateusz]

More context:

We initially started by ensuring that the single point energies are the same for lambda 0 and 1 for three molecules. This was done in the vacuum discharge leg. For example at l=0 dimethylsulfide in openmm gave 3.110 kcal/mol, and SOMD gave 3.109 kcal/mol, and for l=1 it was -0.093 in both. The energies were also equivalent for methylbenzoate and chloromethane. Whereas for lambda=0.5 we found different values, according to our checks, it causes no issues for the final results - the scaling of electrostatics is "accelerated".

We computed the full SOMD protocol for several molecules. Our reference is absolv and openmmtools (which we also partly verified against gromacs). With the following results (kcal/mol):

• DMSO absolv vs SOMD, agreement to two decimal places: -10.15.
• dimethylsulfide: SOMD: -2.429477 absolv: -2.54
• chloromethane: SOMD: 1.150995 absolv: ~ 1.083
• ethyl propionate: SOMD: -3.203012 absolv: -3.293
  Which are the good cases. But we also have some bad cases:
• ethanol: SOMD: -8.155983, absolv: -8.547
• methylbenzoate: SOMD -1.791851, absolv: -4.449

After that we focused on methylbenzoate. We suspect that the issue goes back to the size of the molecule and the number of intramolecular interactions, where methylbenzoate dominates. Specifically, we are now focusing on the reaction field correction script FUNC.py as the potential culprit.

As an approximation, we removed the virtual sites in methylbenzoate and moved the charges to the oxygen manually, and took the new value from the FUNC.py. That shows a very good agreement: SOMD: -4.390081, absolv: -4.449.

So the question now is how to best adjust FUNC.py, which currently doesn't load any virtual sites, or if you think that the issue might be somewhere else.

[jmichel80]

hi @bieniekmateusz as a sanity check you can repeat the vacuum calculations using a non periodic cutoff and a reaction field with the same dielectric constant as in the solvated leg. This will give a consistent description of intramolecular non-bonded electrostatic interactions between the two phases and make FUNC.py redundant. It probably doesn't make a big difference to use either approach for such small molecules in the gas phase. FUNC is not needed for binding free eneryg cycles.

[bieniekmateusz]

Thanks @jmichel80 , switching to the reaction field to avoid the correction has resolved the problem for methylbenzoate. Therefore it appears that indeed the issue was in the FUNC.py.