Converting from Amber to LAMMPS with flexible water

[nickhine]

I am trying to use Intermol's convert.py to convert sets of GAFF parameters for organic molecules in solvent into inputs for LAMMPS. This appears to be mostly working, in terms of outputs from Amber and LAMMPS agreeing, apart from the small discrepancies associated with Coulomb parameters noted in your 2016 paper.

However, I run into a sudden problem when the solvent is water, and the numbers do not remotely agree between Amber and LAMMPS. I believe I can identify the cause: the water is being treated as rigid in LAMMPS, and thus not including any bonded terms, whereas it is flexible (ntf=1, ntc=1) in Amber.

I can see from the code that there is a #ifdef FLEXIBLE inside the intermediate GROMACS file which is meant to switch between flexible and rigid water, but I can't figure out how it is meant to be used/activated and it does not seem to be mentioned in the documentation. Can anyone advise on how to ensure the end result in LAMMPS comes out flexible?

[mrshirts]

That's a good question. That doesn't quite make sense. Can you upload a minimal file that causes the problem, and I can try to debug.

[nickhine]

Thanks - the most minimal version of this I can come up with is to create a single water molecule from a pdb then use antechamber to parameterise it (ie treat it as an unknown organic molecule, bypassing any cleverness in Amber related to water) then Intermol to convert to LAMMPS via GROMACS.

This script does all of that (paths will need tweaking on the python command at the end) and I attach the output I get from it from Intermol comparing energies. No agreement at all, for fairly obvious reasons when you look at the lammps output: it has applied SHAKE constraints (which I did not intentionally request) to the water.

I am running AmberTools 16 and the most recent versions of LAMMPS and Intermol.

lammps_stdout.txt
amber.out.txt
output.txt
script.txt

[nickhine]

I discovered a way to hack it to produce a LAMMPS input that looks mostly correct, namely by using the preprocessor to do "cpp -DFLEXIBLE" on the intermediate gromacs outputs and then reconverting from gromacs to lammps. Unfortunately then I can't do a side-by-side comparison, which is one of the nicest features of Intermol.

I imagine there must be a canonical way of requesting the intermol gromacs module to apply the "flexible" #define, but I can't see from the code how it should work.

Even then, with water I see discrepancies in the nonbonded terms that don't occur for systems containing other solvents. These are very small for a single water molecule but become enormous for a large box full of water (~9000 atoms):

Energy group summary
type input(amber) output (lammps) diff (lammps)

Not comparable energies: are likely not to be the same

   coulomb total     26026.76697440   978748.55968000   952721.79270560
      coulomb-14       420.94261680        0.00000000     -420.94261680
          proper        32.05571600       32.05567416       -0.00004184
        vdw (LR)         0.00000000     -615.18464944     -615.18464944
       vdw total    -12999.35035120   -14451.32537744    -1451.97502624
          vdw-14        51.70294320        0.00000000      -51.70294320

---

Comparable energy terms: these should be very close

           angle      9289.08751680     9289.08709840       -0.00041840
            bond       336.29318400      336.29310869       -0.00007531
          bonded      9657.43641680     9657.43588125       -0.00053555
        dihedral        32.05571600       32.05567416       -0.00004184
       nonbonded     13027.41662320   -24421.83478240   -37449.25140560
       potential     22684.85304000   -14764.39920240   -37449.25224240

---------------- Total Potential Energy Comparison --------------------
Input amber potential energy: 22684.85304000
Difference in potential energy from amber=>lammps conversion: -37449.25224240

[mrshirts]

Can you post the actual amber input files? That will help ensure replicatability. Somehow, it's using a TIP3P file that includes an #ifdef statement, which is likely one from the GROMACS installation, or that ParmEd is generating.

This may actually be an issue with ParmEd, since InterMol uses ParmEd to convert from AMBER to GROMACS - you might look at the instructions there as well.

[nickhine]

The file "script.txt" in the post above generates an amber input file for water (the script writes the pdb then the tleap file).

The primary issue is that the GROMACS output hedges its bets as to whether you want rigid water or flexible water, but when doing a Amber -> Gromacs -> LAMMPS conversion you have no option to choose and the final LAMMPS output doesn't have the option.

However, the issue with the nonbonded terms is different, and more alarming as it still occurs even when the rigid/flexible issue is fixed by the method I mentioned above. Here are some amber inputs which show that up. I converted them just now by hard-wiring defines={'FLEXIBLE': True} into the Intermol gromacs module (but that is obviously a hack!).

AQS_watr_solv.prmtop.txt

AQS_watr_solv.crd.txt

[mrshirts]

The file "script.txt" in the post above generates an amber input file for water (the script writes the pdb then the tleap file).

Yeah, I wanted to avoid any potential issues with differences there.

The primary issue is that the GROMACS output hedges its bets as to whether you want rigid water or flexible water, but when doing a Amber -> Gromacs -> LAMMPS conversion you have no option to choose and the final LAMMPS output doesn't have the option.

Right, the question is whether this is something ParmEd hardcodes or not. It shouldn't be making assumptions about what the GROMACS .mdp may or may not have . . . but one might still have to work around that.

However, the issue with the nonbonded terms is different, and more alarming as it still occurs even when the rigid/flexible issue is fixed by the method I mentioned above.

The most likely issue is that the exclusions are not being set correctly in either AMBER or LAMMPS. PROBABLY the AMBER energy isn't being calculated correctly, since a TIP3P or SPC/E water box should not have positive potential energy, and including H-H interactions would create big positive energies.

I'll see what I can look at in terms of actual debugging, but in the middle of COVID-19 class preparation, so debugging time limited.