Pressure issue?

[Dan-Burns]

Really great work here, thank you! Very nice to have such an easy way to run hremd with openmm.

I generated my openmm system outside of femto and ran a short hremd simulation following your tutorial. I noticed that the water diffuses far outside of the box boundaries. The system has a MonteCarloBarostat set to 1 bar and the box vectors are correct.

I went back and ran a single simulation with the following

structure = pmd.load_file('structure.pdb', structure=True)

system = pmd.load_file('system.xml')

# easier for me to get the atom indices with MDAnalysis
u = mda.Universe('structure.pdb')

import femto.md.simulate
stages = [femto.md.config.Simulation(
        integrator=femto.md.config.LangevinIntegrator(
            timestep=4.0 * openmm.unit.femtosecond,
            
        ),
        temperature=293.15 * openmm.unit.kelvin,
        pressure=1.0 * openmm.unit.bar,
        n_steps=50000,
    )]
rest_config = femto.md.config.REST(scale_torsions=True, scale_nonbonded=True)

indices = u.select_atoms('protein').atoms.ix
solute_idxs = set(indices)

femto.md.rest.apply_rest(system, solute_idxs, rest_config)
state = {femto.md.rest.REST_CTX_PARAM: 1.0}

final_coords = femto.md.simulate.simulate_state(
    system, structure, state, stages, femto.md.constants.OpenMMPlatform.CUDA
) # pass empty dictionary instead of state
with open(f'output/test.pdb', 'w') as f:
    openmm.app.PDBFile.writeFile(structure.topology, final_coords.getPositions(), f)

The system still drifts apart with this. I ran the same system using openmm directly and it behaves as expected. Do you have any idea what's going on?

Thank you,
Dan

Starting Structure

100ps directly with openmm

200ps with femto code above

[SimonBoothroyd]

Hi @Dan-Burns, I'm glad the framework looks like it will be useful to you!

I think there could be two things going on:

1 - the system wasn't simulated with PBC switched on

I don't think this will be the case here, but for reference you can easily check with something like:

method_lookup = {
    openmm.NonbondedForce.NoCutoff: "NoCutoff",
    openmm.NonbondedForce.CutoffNonPeriodic: "CutoffNonPeriodic",
    openmm.NonbondedForce.CutoffPeriodic: "CutoffPeriodic",
    openmm.NonbondedForce.Ewald: "Ewald",
    openmm.NonbondedForce.PME: "PME",
    openmm.NonbondedForce.LJPME: "LJPME",

}
nonbonded_methods = [
    method_lookup[f.getNonbondedMethod()]
    for f in system.getForces()
    if isinstance(f, openmm.NonbondedForce)
]

print(nonbonded_methods)

and the output should be something like ['PME']

2 - by default OpenMM does not enforce PBC when getting the current coordinates from a context, and so even though PBC are being correctly applied during the simulation when computing energies, the returned coordinates will effectively be unwrapped.

In this case, you will currently likely need to apply the PBC to the coordinates manually. A quick an dirty way to do this purely using OpenMM would be:

context_tmp = openmm.Context(system, openmm.VerletIntegrator(0.001))
context_tmp.setPeriodicBoxVectors(*final_coords.getPeriodicBoxVectors())
context_tmp.setPositions(final_coords.getPositions())

wrapped_coords = context_tmp.getState(getPositions=True, enforcePeriodicBox=True)

with open('wrapped.pdb', 'w') as f:
    openmm.app.PDBFile.writeFile(topology, wrapped_coords.getPositions(), f)

I can also looking into exposing the enforcePeriodicBox flag on future versions of the simulate_state method if that would be helpful?

Let me know if this answers your question!

[Dan-Burns]

Thank you @SimonBoothroyd - the latter solution took care of it. I guess when using a DCDReporter, the coordinates are adjusted by default but femto is handling the writing to DCDFile. That would be great to have the option.