Problem with Gromacs to Amber topology conversion

[bernieleecy]

I tried to convert a solvated and minimised Gromacs topology to Amber format using BioSimSpace, with the water conversion fix described in #52 applied.

import BioSimSpace as BSS
from Sire.IO import setAmberWater
from Sire.Mol import MGName

system = BSS.IO.readMolecules(['complex_ions.gro', 'topol.top'])
waters = setAmberWater(system._sire_object.search("water"), "tip3p")
system.removeWaterMolecules()
for wat in waters:
    system._sire_object.add(wat, MGName("all"))

BSS.IO.saveMolecules("bromo", system, ["rst7", "prm7"])

These files were then used for simulations in Amber 20, i.e. BSS was used to convert the topology only. I equilibrated the system (NVT then NPT, total of 600 ps of restrained equilibration), then ran a short unrestrained MD simulation (1 ns). In the unrestrained MD simulation, the protein was unstable (started to unfold), and the 1-4 EEL value looks odd. The energy values at the end of the simulation were:

 NSTEP =   500000   TIME(PS) =    1600.000  TEMP(K) =   300.95  PRESS =    16.6
 Etot   =    -96104.5192  EKtot   =     20623.6641  EPtot      =   -116728.1833
 BOND   =       364.0951  ANGLE   =      1108.9572  DIHED      =      1233.2389
 1-4 NB =       392.5921  1-4 EEL =       -95.3411  VDWAALS    =     14511.8637
 EELEC  =   -134243.5891  EHBOND  =         0.0000  RESTRAINT  =         0.0000
 EKCMT  =      9654.4261  VIRIAL  =      9531.5242  VOLUME     =    341930.5524
                                                    Density    =         1.0085

I went back and checked the energies of the starting topologies. For the unconverted topology, the energies were as follows (from gmx mdrun -rerun):

Energy                      Average   Err.Est.       RMSD  Tot-Drift
-------------------------------------------------------------------------------
Bond                        267.273         --          0          0  (kJ/mol)
Angle                       904.499         --          0          0  (kJ/mol)
Proper Dih.                  4637.5         --          0          0  (kJ/mol)
Improper Dih.               70.5391         --          0          0  (kJ/mol)
LJ-14                       1726.73         --          0          0  (kJ/mol)
Coulomb-14                  18689.7         --          0          0  (kJ/mol)
LJ (SR)                     69915.4         --          0          0  (kJ/mol)
Disper. corr.              -4258.81         --          0          0  (kJ/mol)
Coulomb (SR)                -638272         --          0          0  (kJ/mol)
Coul. recip.             -2.46847e-14         --          0          0  (kJ/mol)
Potential                   -546319         --          0          0  (kJ/mol)

For the converted topology, the energies were as follows (from sander). The main difference I could see was in the bonded terms, rather than the electrostatic terms.

   NSTEP       ENERGY          RMS            GMAX         NAME    NUMBER
      1      -1.2936E+05     1.4984E+01     6.0450E+01     O        3768

 BOND    =      371.2039  ANGLE   =      216.1800  DIHED      =     1125.2485
 VDWAALS =    15693.9659  EEL     =  -151649.8146  HBOND      =        0.0000
 1-4 VDW =      412.6990  1-4 EEL =     4466.9648  RESTRAINT  =        0.0000

I also tried to convert the Amber topology back to Gromacs format, and checked the energies again. The energies were very similar to the original Gromacs topology.

Energy                      Average   Err.Est.       RMSD  Tot-Drift
-------------------------------------------------------------------------------
Bond                        267.274         --          0          0  (kJ/mol)
Angle                       904.499         --          0          0  (kJ/mol)
Proper Dih.                 4708.04         --          0          0  (kJ/mol)
LJ-14                       1726.73         --          0          0  (kJ/mol)
Coulomb-14                  18690.4         --          0          0  (kJ/mol)
LJ (SR)                     69946.3         --          0          0  (kJ/mol)
Disper. corr.              -4257.32         --          0          0  (kJ/mol)
Coulomb (SR)                -638272         --          0          0  (kJ/mol)
Coul. recip.             -2.46847e-14         --          0          0  (kJ/mol)
Potential                   -546286         --          0          0  (kJ/mol)

From this, I'm not sure what's causing the instability in MD simulations. I originally thought it was the dihedrals, but the single point calculations on the minimised structure suggest that's not the case?

key_files_bss.zip

[lohedges]

Hi there. Thanks for posting.

You shouldn't need to manually convert the water topology, this will be done for you on write by BioSimSpace.IO.saveMolecules.

For your AMBER simulations that are unfolding, are you using sander or pmemd? We recently discovered and fixed an issue with zero periodicity torsions that need to be treated differently in pmemd. In our own simulations we also experienced protein unfolding. See here for details. Could you report what version of Sire you are using, so we know that you have the fixed applied.

import Sire
print(Sire.__version__)

For the difference in the bonded terms: It's not possible to directly compare AMBER and GROMACS bond energies due to the way in which they treat three-point water models. To do so, you'd need to make sure that they are both implementing SHAKE in the same way, which might mean tweaking the config file options. See this thread for a discussion.

I'll try to reproduce your numbers here when I get a chance.

[bernieleecy]

I've been using pmemd and Sire 2020.1.0. I've also been using the stable release of BioSimSpace rather than the latest development version, so that might be why I needed to fix the waters manually.

I'll download the development version and test it and let you know if there are any further issues.

[lohedges]

Thanks for the update. You'll need the latest version of Sire, 2021.1.0. Let me know if this doesn't work for you and I'll look into it further.

Another think that I thought of is that AMBER can have issues if molecules in the topology are wrapped across the periodic image. I've not checked your topology, but it might be worth having a look if it still doesn't work with the most recent version of Sire. MDAnalysis has functionality to unwrap the coordinates if needed.

Cheers

[lohedges]

This also reminds me that I need to push out a new release of BioSimSpace sometime soon!

[bernieleecy]

It's fine with Sire 2021.1.0, I ran a 10 ns simulation and the protein was stable.

Thank you for your help!

[lohedges]

No problem, glad to hear that it's now working. This was a really peculiar issue that caused us a lot of headaches to track down and solve.

Thanks again for reporting!