MD error with FESetup output

[ppxasjsm]

Running MD simulations from minimised and equilibrated FESetup output crashes using sander from amber2018.

The error zip file is attached. The code that lead to the error is:

import BioSimSpace as BSS
system = BSS.IO.readMolecules(['FESetup2/_complexes/cats:CatS_155/solvated.rst7', 'FESetup2/_complexes/cats:CatS_155/solvated.parm7'])
protocol = BSS.Protocol.Production(runtime=20*BSS.Units.Time.nanosecond, frames=1000)
process = BSS.MD.run(system, protocol)
print(process.isError())
process.getOutput()

md.zip

[lohedges]

Do you have the unmodified FESetup output files (rather than the ones written by BioSImSpace). If those don't run, then it's an issue with the configurations, rather than BioSImSpace (although we might be able to modify the protocol to help). I usually see crashes like this when the system hasn't been minimised correctly.

[lohedges]

For reference, sander can't even perform a miminisation using these input files (at least with SHAKE on, which is what I was told we want to use).

[lohedges]

It looks like the original configuration is a bunch of junk (I don't know what FESetup uses for minimisation/equilibration). If I disable SHAKE and re-minimise the system using BioSimSpace, then production MD runs just fine.

[skfegan]

If FESetup is doing minimization/equilibration its default is to use amber16/sander although you can use other programs if you specify the mdengine in the input file. However, for each of the three categories of molecules (ligand, protein, and complex) the minimization/equilibration must be set separately, i.e. if you minimize/equilibrate the ligand or protein it does not automatically minimize/equilibrate the complex. fesetup/_complexes/*/solvated.rst is the input for minimization not the result.

[jmichel80]

Hi Lester,

FYI here is the more or less default minimisation/equilibration protocol used by FESetup. Within fesetup this would led to an equilibrated rst7 file called md0006.rst7 or similar. As Sarah pointed out solvated.rst7 is the leap output before any minimisation and equilibration.

Just doing energy minimisation with BSS will still leave you with incorrect densities.

BW

Julien

• For solvated ligands

1) We carry out a restrained energy minimisation first. The idea is to seek to preserve the input binding mode by allowing initially only the solvent to relax around the ligand

min.nsteps = 200 # up to 200 steps of energy minimisation

min.ncyc = 100 # the first 100 steps are steepest descent

min.restr_force = 10.0 # atoms subject to the restraint are restrained at 10 kcal/mol/A**2

min.restraint = notsolvent # all atoms that are not solvent particles are restrained

2) heat the system to the final temperature running NVT

md.heat.nsteps = 1000 # This is 2ps with a 2fs timestep

md.heat.T = 300.0 # Will get to 300 K by the end of the 2 ps

md.heat.restraint = notsolvent # Still using restraints on everything but solvent at this stage

md.heat.restr_force = 5.0

3) fix the density of the system running NpT

md.press.nsteps = 5000 # 10 ps is usually sufficient to get stable box density with pmemd. With openMM (Monte Carlo Barostat) this may not be the case ! May need some tests.

md.press.T = 300.0

md.press.p = 1.0

md.press.restraint = notsolvent

md.press.restr_force = 4.0

4) restraints release in 4 steps, this is a NpT protocol

md.relax.nrestr = 4

md.relax.nsteps = 500

md.relax.T = 300.0

md.relax.p = 1.0

md.relax.restraint = notsolvent

• For protein-ligand complexes

Minimisation

min.nsteps = 2000

min.ncyc = 1000

min.restr_force = 10.0

min.restraint = :LIG # The ligand is restrained initially, we try to preserve the input binding mode. May need adjustments with the OpenMM energy minimiser (check)

heat the system to the final temperature running NVT

md.heat.nsteps = 1000 # Should be enough with an Andersen thermostat in OpenMM to get to 300 K (check)

md.heat.T = 300.0

md.heat.restraint = bb_lig

md.heat.restr_force = 5.0

fix the density of the system running NpT

md.press.nsteps = 5000 # In theory is system dependent, but experience suggests 10 ps is ok to get box density about right with sander or pmemd. Again have to check efficiency of OpenMM Barostat.

md.press.T = 300.0

md.press.p = 1.0

md.press.restraint = :LIG

md.press.restr_force = 4.0

restraints release in 4 steps, this is a NpT protocol

md.relax.nrestr = 4

md.relax.nsteps = 500

md.relax.T = 300.0

md.relax.p = 1.0

md.relax.restraint = :LIG

---

Dr. Julien Michel, Senior Lecturer Room 263, School of Chemistry University of Edinburgh David Brewster road Edinburgh, EH9 3FJ United Kingdom phone: +44 (0)131 650 4797 http://www.julienmichel.net/

On Mon, Oct 22, 2018 at 5:08 PM Sarah Fegan notifications@github.com<mailto:notifications@github.com> wrote:

If FESetup is doing minimization/equilibration its default is to use amber16/sander although you can use other programs if you specify the mdengine in the input file. However, for each of the three categories of molecules (ligand, protein, and complex) the minimization/equilibration must be set separately, i.e. if you minimize/equilibrate the ligand or protein it does not automatically minimize/equilibrate the complex. fesetup/_complexes/*/solvated.rst is the input for minimization not the result.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHubhttps://github.com/michellab/BioSimSpace/issues/36#issuecomment-431881140, or mute the threadhttps://github.com/notifications/unsubscribe-auth/ALLd5Ol9uk_VH0_xM8DFuNmCgNnpUsQBks5une1zgaJpZM4XzfZE.

The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336.

[lohedges]

Thanks for the clear explanation of the protocol, that's interesting.

Sorry for any confusion, I wasn't trying to imply that you can just minimise with BioSimSpace and be good to go, only that minimisation fixes whatever steric clashes are left after the FESetup step. I just wanted to do the simplest thing to see if I could then run any protocol using the re-minimised system (since the FESetup configurations crashed within the first few iterations).

[lohedges]

If I am reading the above correctly, it sounds like we can't rely on SOMD, i.e. OpenMM, to correctly minimise or equilibrate the system. At the moment I'm making using of the minimisation and equilibration options in SOMD during my free energy simulations, since it is currently the only MD engine that can run dynamics with a system containing perturbed molecules.

For the purposes of the D3R challenge I could hack the other engines to write atom properties at lambda = 0, then copy the updated coordinates back into the original system following minimisation/equilibration, but this would take a bit of work. (I discussed something like this with Christopher last week.)

[jmichel80]

Hi Lester,

Until now we have used FESetup to equilibrate unperturbed molecules with the above protocol. Hannes was able to map the above described protocol to implementations in various MD engines (most will support some form of positional restraints and can do energy minimisation, NVT, NPT).

We don't do any minimisation/equilibration of the perturbed molecules within FESetup. FESetup merely uses the equilibrated coordinates to assemble a system around a perturbed molecule. This is why we tend to start all our SOMD simulations from FESetup inputs with an energy minimisation (this time supporting perturbed molecules at any lambda value). We tend not to do 'equilibration' but we discard the first few % of each lambda trajectory, which effectively amounts to an equilibration.

I suggest you use the existing generic MD capabilities of BioSimSpace to implement equilibration pre free energy input generation. Anything that involves equilibration for a perturbed molecule should be folded in the free energy protocol because they will be free energy package specific.

Best wishes,

Julien

---

Dr. Julien Michel, Senior Lecturer Room 263, School of Chemistry University of Edinburgh David Brewster road Edinburgh, EH9 3FJ United Kingdom phone: +44 (0)131 650 4797 http://www.julienmichel.net/

On Tue, Oct 23, 2018 at 9:15 AM Lester Hedges notifications@github.com<mailto:notifications@github.com> wrote:

If I am reading the above correctly, it sounds like we can't rely on SOMD, i.e. OpenMM, to correctly minimise or equilibrate the system. At the moment I'm making using of the minimisation and equilibration options in SOMD during my free energy simulations, since it is currently the only MD engine that can run dynamics with a system containing perturbed molecules.

For the purposes of the D3R challenge I could hack the other engines to write atom properties at lambda = 0, then copy the updated coordinates back into the original system following minimisation/equilibration, but this would take a bit of work. (I discussed something like this with Christopher last week.)

— You are receiving this because you commented. Reply to this email directly, view it on GitHubhttps://github.com/michellab/BioSimSpace/issues/36#issuecomment-432143769, or mute the threadhttps://github.com/notifications/unsubscribe-auth/ALLd5BWJtQBxZRwTIG-EF-rs4Ozl6rzOks5untApgaJpZM4XzfZE.

The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336.

[lohedges]

Great, this is pretty much what I'm doing anyway. It sounds like I'll need to add another keyword option to my protocols to restrain all non-solvent atoms. I'll also need to deal with writing input files for the lambda=0 state (using the property map) then just re-adding the updated solvent coordinates back into the original system.