Issues with Gromacs Input Files

Hello,

I'm running BSS version 2022.1.0, 8.g3f6b1c18 on linux, installed using conda. I'm attempting to run a short simulation with the gromacs input files supplied here. All relevant inputs/ notebooks are here.

I've run into a few issues:

Failure of getSystem(block=True). During equilibration steps (see bdr4_equilibration.ipynb) this command sometimes fails the first time it is run. However, it runs correctly the second time (without repeating the simulation). The system was also unexpectedly unstable (LINCS warnings) and exploded during "PMEMD NPT equilibration for 2ns without restraints" (Fatal error: 3 particles communicated to PME rank 1 are more than 2/3 times the cut-off out of the domain decomposition cell of their charge group in dimension x. This usually means that your system is not well equilibrated.) . I had similar issues when using sander. The input was the ligand 1 complex for brd4.
Single point energy calculation works with Gromacs but fails with Amber (see gh_issue.ipynb). Error reads: "EXTRA POINTS: nnb too small!". Input was the ligand 1 complex for jnk1.
BSS fails to read Amber inputs generated with ParmEd (see bottom of gh_issue.ipynb). I converted the jnk1 input (ligand 1 complex) from gro/top to parm7/rst7 according to the discussion here. However, upon attempting to reload the parm7/rst files, I get [OSError: Failed to read molecules from ['jnk1_system.parm7', 'jnk1_system.rst7']. It looks like you failed to include a topology file.]()

Thanks very much.

Hi there,

For part of question 1): I'm not sure why this happens but I've recently fixed this (last) week by simply re-trying the trjconv command if it happens to fail the first time. (Weirdly the failures only ever happen within getSystem() when calling a subprocesses, i.e. using the exact same commands on the command-line works every time!).

I should be able to check the other issues at some point tomorrow.

Cheers.

Hi,

Thanks for the quick response. Great - I'll update to the latest version.

Cheers

Out of interest, is there are reason why you need to use ParmEd to convert the GROMACS files, rather than using BioSImSpace directly, i.e.:

import BioSImSpace as BSS

# Load GROMACS system.
system = BSS.IO.readMolecules(["jnk1_complex.gro", "jnk1_complex.top"])

# Write to AMBER format.
BSS.IO.saveMolecules("test", system, ["prm7", "rst7"])

(This is what would be done behind the scenes if you loaded GROMACS files and tried to run a simulation with AMBER.)

This works, but the resulting topology file looks quite different to the one generated by ParmEd. A single-point calculation would test whether the information is self-consistent, despite the differences.

@jmichel80

The original issue I was having was instability during equilibration using sander through BSS and @jmichel80 mentioned that there had been issues previously with file conversion. To rule this out, I wanted to repeat the file conversion using ParmEd and re-run the equilibration to see if I got the same issues. Unfortunately I've not been able to run a single-point calculation using the BSS-converted files (see first part of see gh_issue.ipynb), or to load in the files produced by ParmEd to run a single-point calculation on those.

Hi again,

I'm really not sure what's going on with those GROMACS files provided in the repository you linked to. Out of interest, I looked at their example notebook which shows how the files are generated. If I reproduce their code, but instead write to AMBER format files rather than GROMACS, then everything works, i.e. I can load back the files with BioSimSpace and I get the same results when performing single-point energy calculations with AMBER and GROMACS.

For reference, they are just parameterising the protein using FF14SB and "openff_unconstrained-2.0.0" for the ligand. It isn't possible to do this directly in BioSimSpace since the protein PDB contains missing atom types and the SDF file is required in order to provide the required stereochemistry information for the ligand.

Here's my script. (You'll just need to clone the abfe-benchmark repo to your working directory.)

import BioSimSpace as BSS

import parmed

from openff.toolkit.topology import Molecule
from openff.toolkit.typing.engines.smirnoff import ForceField as OFF_ForceField

try:
    from openmm import XmlSerializer, app, unit
    from openmm.app import HBonds, NoCutoff, PDBFile, ForceField
except ImportError:
    from simtk import unit
    from simtk.openmm import XmlSerializer, app
    from simtk.openmm.app import HBonds, NoCutoff, PDBFile

# Load in the ligand.
ligand_molecule = Molecule("abfe-benchmark/structures/jnk1/ligand1.sdf")

# Specify the "Sage" forcefield.
force_field = OFF_ForceField("openff_unconstrained-2.0.0.offxml")

# Parametrize the ligand molecule by creating a Topology object from it.
ligand_system = force_field.create_openmm_system(ligand_molecule.to_topology())

# Read in the coordinates of the ligand from the PDB file.
ligand_pdbfile = PDBFile("abfe-benchmark/structures/jnk1/ligand1.pdb")

# Convert the ligand system to a ParmEd object.
ligand_parmed_structure = parmed.openmm.load_topology(ligand_pdbfile.topology,
											          ligand_system,
											          ligand_pdbfile.positions)

# Write the ligand structure to AMBER format.
ligand_parmed_structure.save("ligand1.rst7", overwrite=True)
ligand_parmed_structure.save("ligand1.parm7", overwrite=True)

# Parse the protein PDB file.
protein_pdbfile = PDBFile("abfe-benchmark/structures/jnk1/protein.pdb")

# Load the AMBER protein force field through OpenMM.
omm_forcefield = app.ForceField("amber14/protein.ff14SB.xml", "amber14/tip3p.xml")

# Parameterize the protein.
protein_system = omm_forcefield.createSystem(protein_pdbfile.topology,
                                             nonbondedCutoff=1*unit.nanometer,
                                             nonbondedMethod=app.NoCutoff,
                                             constraints=None,
                                             rigidWater=False)

# Convert the protein System into a ParmEd Structure.
protein_parmed_structure = parmed.openmm.load_topology(
                                        protein_pdbfile.topology,
                                        protein_system,
                                        xyz=protein_pdbfile.positions)

# Write the protein structure to AMBER format.
protein_parmed_structure.save("protein.rst7", overwrite=True)
protein_parmed_structure.save("protein.parm7", overwrite=True)

# Load the ligand and protein with BioSimSpace.
ligand = BSS.IO.readMolecules("ligand1.*7")[0]
protein = BSS.IO.readMolecules("protein.*7")[0]

# Combine to form a complex.
complx = (ligand + protein).toSystem()

# Create a single-step minimisation protocol.
protocol = BSS.Protocol.Minimisation(steps=1)

# Create processes for AMBER and GROMACS simulations.
process_amb = BSS.Process.Amber(complx, protocol)
process_gmx = BSS.Process.Gromacs(complx, protocol)

# Modify the GROMACS config to perform zero steps.
config = process_gmx.getConfig()
config[1] = "nsteps = 0"
process_gmx.setConfig(config)

# Run the processes.
process_amb.start()
process_amb.wait()
process_gmx.start()
process_gmx.wait()

# Compare energies. (Where direct comparisons can be made.)
print(f"Bond energies... AMBER = {process_amb.getBondEnergy()}, GROMACS = {process_gmx.getBondEnergy()}")
print(f"Angle energies... AMBER = {process_amb.getAngleEnergy()}, GROMACS = {process_gmx.getAngleEnergy()}")
print(f"Dihedral energies... AMBER = {process_amb.getDihedralEnergy()}, GROMACS = {process_gmx.getDihedralEnergy()}")

This outputs:

Bond energies... AMBER = 202.0904 kcal/mol, GROMACS = 202.0662 kcal/mol
Angle energies... AMBER = 1111.9015 kcal/mol, GROMACS = 1111.8045 kcal/mol
Dihedral energies... AMBER = 4336.2379 kcal/mol, GROMACS = 4315.6788 kcal/mol

I'm not yet sure if the issue with the GROMACS files is our ability to read them, or whether they were invalid in the first place. I'll try to reproduce the above writing to GROMACS format, i.e. using the output of the script, rather what was uploaded to the benchmark repository.

Cheers.

I've just re-run the above script, instead using ParmEd to write to GROMACS format. As you've found, on conversion to AMBER format within BioSimSpace, any AMBER simulation fails with the EXTRA POINTS: nnb too small!" message. While the GROMACS simulations work, the single point energies are different (but close) to those that I obtained above using the the AMBER files written by ParmEd. The output was as follows:

Bond energies... AMBER = None, GROMACS = 241.7304 kcal/mol
Angle energies... AMBER = None, GROMACS = 1124.9594 kcal/mol
Dihedral energies... AMBER = None, GROMACS = 4337.7390 kcal/mol

Given the inconsistencies, I'm not sure whether the issue is with ParmEd or us.

@lohedges

That's very helpful @lohedges ! @fjclark could you check whether the AMBER prm files produced by @lohedges 's script work with SOMD ? We can feedback to Bayer next week the above issues.

@lohedges

Thanks very much @lohedges. @jmichel80 will do.

@jmichel80

@jmichel80 I can confirm that the AMBER files produced by @lohedges's script work with SOMD. However, test alchemical simulations crash when I use the same config file I'm using for my current work, giving "NaN or Inf has been generated along the simulation" during minimisation. This turned out to be due to the use of "constraint = allbonds". Do you have any suggestions as to why this is happening, given that the system appeared well equilibrated (I have used an equilibration scheme almost identical to the one which produced input which works with this config file).

All relevant input is here. I have:

-Produced AMBER format files for brd4, lig1 complex using @lohedges 's script (this failed at first due to this issue - despite having compatible versions of ParmEd and OpenMM in my environment, an incompatible version of ParmEd was bundled with AmberTools, which worked after manually modifying decorators.py).

-Equilibrated using pmemd and pmemd.cuda - see input.

-Run a short production run using SOMD through BSS - successful

-Run single alchemical windows using my standard config file with (failure) and without (success) constraint set to "allbonds".

Thanks

This is odd, constraints are disabled before energy minimisation. See https://github.com/michellab/Sire/blob/devel/wrapper/Tools/OpenMMMD.py#L1507-L1520

Also your single point energy before minimisation is identical before minimisation between the two runs, so I expect the same input was passed to integrator.minimiseEnergy()

Is the error reproducible ?

You could try disabling minimisation, that step may not be necessary for an ABFE run started from a structure equilibrated at discharge = 0.0

I've repeated the calculations, only adding or removing "set_constraints = allbonds" to the config. Adding this consistently causes the simulations to crash during minimisation. I'm running my modified version of SOMD with Boresch restraints, but I've repeated this with the latest version of sire as contained within BSS-dev - I observe the same thing.

When disabling minimisation, I get:

  2 NaN or Inf has been generated along the simulation
  3 Starting /home/finlayclark/sire_boresch_2022.app/bin/somd-freenrg: number of threads equals 20
  4 srun: error: ishikawa: task 0: Exited with exit code 255

Very similar to with minimisation, only no output related to setting up the system is displayed.

Out of interest, do the original GROMACS files work with SOMD? (For example, if you just perform minimisation and equilibration with GROMACS rather than PMEMD.) It would be interesting to know if you experience the same constraint related crashes with those, too.

Could you setup from scratch the input files with BioSimSpace (starting from pdb files extracted from the pmemd equilibration and solvating in BioSimSpace) to verify if SOMD runs without errors in that case ? If so comparing the prm7 files should indicate a discrepancy in the handling of bonded terms.

@lohedges

@lohedges I'm not able to equilibrate the original GROMACS files in BSS using GROMACS - I get failures during my final NPT equilibration. Based on this, it seems likely that the issue is with the original files.

@lohedges

@lohedges , I was also wondering about comparing bond energies. In your comment above (#289 (comment)) you use these when comparing single-point energies. However, when I convert one of these systems (ligand 1 complex) loaded from the same (AMBER) input to GROMACS, I find that the bond energy changes by over 100 kcal / mol (although the angle and dihedral energies are unaffected).

Bond energies... AMBER = 324.5188 kcal/mol, GROMACS = 84.6532 kcal/mol
Angle energies... AMBER = 421.5360 kcal/mol, GROMACS = 421.5368 kcal/mol
Dihedral energies... AMBER = 1470.6892 kcal/mol, GROMACS = 1464.7347 kcal/mol

Is this unexpected or does this just depend on the details of the interconversion?

Thanks.

Bond energies only compare in vacuum, or for comparable water models. AMBER and GROMACS handle rigid waters differently, one with explicit H-H bonds, the other with constraints.

OK, thank you.

I'm not sure why the AMBER process is returning None for the bond energy, will try to debug that tomorrow. Other than water molecules, the number of bonded terms are the same, i.e.:

import BioSimSpace as BSS

s_gro = BSS.IO.readMolecules("fixed_ligand1/*)

# Check bonds in the ligand.
s_gro[0]._sire_object.property("bond")
TwoAtomFunctions( nFunctions() == 59 )

# Check bonds in the first water molecule. (No explicit H-H bond for GROMACS, rigid water handled by SETTLE.)
s_gro[1]._sire_object.property("bond")
TwoAtomFunctions( nFunctions() == 0 )

# Convert to AMBER format.
BSS.IO.saveMolecules("tmp", s_gro, ["prm7", "rst7"])

# Reload.
s_amb = BSS.IO.readMolecules("tmp.*")

# Check bonds in the ligand. (Same as before.)
s_amb[0]._sire_object.property("bond")
TwoAtomFunctions( nFunctions() == 59 )

# Check bonds in the first water molecule. (AMBER rigid water has three bonds.)
s_amb[1]._sire_object.property("bond")
TwoAtomFunctions( nFunctions() == 3

Just to say that I don't think the file is valid. If I look at the output files (both GROMACS and AMBER) it looks like the energy has blown up. Some terms appear sensible (and similar), e.g. angle and dihedral, but the total energy is huge or NaN. For AMBER, you are seeing None for the bond energy because the energy record is just asterisks, i.e.:

anber.nrg_

   NSTEP       ENERGY          RMS            GMAX         NAME    NUMBER
      1              NaN            NaN     2.8472E+05     O        2736

 BOND    = *************  ANGLE   =      107.9843  DIHED      =       19.3600
 VDWAALS =           NaN  EEL     =           NaN  HBOND      =        0.0000
 1-4 VDW =        0.0000  1-4 EEL =        0.0000  RESTRAINT  =        0.0000

For GROMACS, I needed to do some monkeying to get things to run due to an invalid decomposition on my laptop. The log file gives:

gromacs.log

   Energies (kJ/mol)
           Bond          Angle    Proper Dih.        LJ (SR)  Disper. corr.
    5.43645e+01    4.51814e+02    8.10024e+01    3.55627e+17   -2.30155e+02
   Coulomb (SR)   Coul. recip.      Potential Pres. DC (bar) Pressure (bar)
   -2.82337e+04    5.48558e+03    3.55627e+17   -1.01028e+02   -3.15806e+18

Note the huge potential energy. I wonder one/some of the atoms should actually be dummies, i.e. not part of the energy calculation.

I'd forgotten to update the names for hydrogens in the water molecules, so I think it was matching against the wrong atom type. Will do that now to see if it fixes things.

Bah, doing that ends up back with the original error. I think that atom names were correct in the original file, I had been misreading names as types. However, I think it's the ligand name that's causing issues, so will try with unchanged atom names, and updated molecule name, i.e. Ligand.

Nope, still the original error. Sorry, I think the files that I sent were wrong, but somehow magically loaded. I'll add some debugging statements to the Sire.IO.GroTop parser to work out where the crash occurs. That should (hopefully) tell us which record is causing the issue. The problem occurs with the topology file alone (if you do things just in Sire), so I assume that it's a formatting issue there (or something that Sire can't parser.)

Apologies for going round in circles. The joys of debugging at the end of the day!

Not at all - thanks very much for all the work on this!

Bah, it was a Sire file caching issue. Changing the ligand names in the topology from Ligand 2 and Ligand 8 to Ligan fixes things, i.e. it loads. Will check the energies after lunch 🤦‍♂️

Yes, that fixes things. (Thank goodness.)

process_gmx.getAngleEnergy().kj_per_mol()
451.8140 kJ/mol

process_amb.getAngleEnergy().kj_per_mol()
451.8063 kJ/mol

process_gmx.getDihedralEnergy().kj_per_mol()
81.0024 kJ/mol

process_amb.getDihedralEnergy().kj_per_mol()
81.0022 kJ/mol

And the log files:

GROMACS:

   Energies (kJ/mol)
           Bond          Angle    Proper Dih.        LJ (SR)  Disper. corr.
    5.43645e+01    4.51814e+02    8.10024e+01    3.84003e+04   -2.30155e+02
   Coulomb (SR)   Coul. recip.      Potential Pres. DC (bar) Pressure (bar)
   -4.76565e+04    1.50650e+03   -7.62286e+03   -1.01028e+02    7.40276e+04

AMBER:

   NSTEP       ENERGY          RMS            GMAX         NAME    NUMBER
      1      -1.6796E+03     1.5635E+02     3.5811E+03     O        3294

 BOND    =       46.1225  ANGLE   =      107.9843  DIHED      =       19.3600
 VDWAALS =     9169.2413  EEL     =   -11022.3213  HBOND      =        0.0000
 1-4 VDW =        0.0000  1-4 EEL =        0.0000  RESTRAINT  =        0.0000

So, the answer is to not have multi-word residue names in the topology file, and to match the residue name between the gro and top files. Should be something that is easy to enforce at their end, or to patch afterwards.

Brilliant, thanks very much for sorting this. Will pass that on.

Merge pull request michellab#289 from OpenBioSim/fix_288

Backport fix from PR michellab#289. [ci skip]

Merge pull request michellab#290 from OpenBioSim/backport_287_289

Navigation Menu

Repository navigation

Issues with Gromacs Input Files #289

15 remaining items

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issues with Gromacs Input Files #289

Description

Activity

lohedges commented on Mar 29, 2022

fjclark commented on Mar 29, 2022

lohedges commented on Mar 29, 2022

fjclark commented on Mar 29, 2022

lohedges commented on Mar 30, 2022

lohedges commented on Mar 30, 2022

jmichel80 commented on Mar 30, 2022

fjclark commented on Mar 30, 2022

fjclark commented on Mar 31, 2022

jmichel80 commented on Mar 31, 2022

fjclark commented on Mar 31, 2022

lohedges commented on Mar 31, 2022

jmichel80 commented on Mar 31, 2022

fjclark commented on Mar 31, 2022

15 remaining items

fjclark commented on Oct 19, 2022

lohedges commented on Oct 19, 2022

fjclark commented on Oct 19, 2022

lohedges commented on Oct 19, 2022

lohedges commented on Oct 20, 2022

lohedges commented on Oct 20, 2022

lohedges commented on Oct 20, 2022

lohedges commented on Oct 20, 2022

fjclark commented on Oct 20, 2022

lohedges commented on Oct 20, 2022

lohedges commented on Oct 20, 2022

fjclark commented on Oct 20, 2022

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions