Fixed bug in PME with triclinic boxes

[peastman]

Fixes #1486.

[jchodera]

Can you elaborate on the bug? What cases does it apply to, what is the effect, and do we need to do a bugfix update for 7.0?

[jchodera]

Huh, it looks like the OpenCL is failing on travis:

Test project /home/travis/build/pandegroup/openmm
    Start 86: TestOpenCLCustomExternalForceMixed
*** Error in `/home/travis/build/pandegroup/openmm/TestOpenCLCustomExternalForce': free(): corrupted unsorted chunks: 0x0000000001439900 ***
1/4 Test #86: TestOpenCLCustomExternalForceMixed ....***Exception: Other  2.54 sec

[peastman]

I'm not really sure. I found the problem by tracing through the calculation with the CPU and CUDA platforms and identifying where they diverged. I then spent a while trying to create a simplified test case that would reproduce the problem, but I couldn't get it to fail. So there seems to still be something about it I don't fully understand. I'll continue working on it.

[peastman]

I figured out what I was missing: in my test case, every atom was a separate molecule, so they were getting correctly wrapped to the triclinic box at an earlier stage in the calculation. I had to add constraints between all of them so they would be treated as a single molecule, thus ensuring some positions would be outside the box when they arrived at the charge spreading kernel.

I do think we need to issue a patch for this. I'll open a separate issue on that.

[jchodera]

Can you elaborate on the impact of the bug at this stage?

[peastman]

I'd say that any simulation with a triclinic box, PME, and CUDA or OpenCL could potentially have incorrect forces. So it's a serious bug.

[jchodera]

Do you know which range of version numbers were affected?

[cespos]

Is the bug fixed? I was running MD simulations for 18 small organic molecules in different solvents, namely water, octanol, and POPC lipid bilayer (the input files are pasted below) and while the temperature equilibrated at the set value for simulations in water, it did not for simulations in octanol or lipid bilayer. In particular, I set the temperature to 298.15 K and I obtained an average temperature of 310 K for simulations in octanol and 305 K for simulations in the lipid bilayer (see Figure below). Do you know which might be the problem? I am using the 7.3.0 version of openMM.

`
#------------------------

Water and Octanol

#------------------------
platform = Platform.getPlatformByName('CPU')
system = prmtop.createSystem(nonbondedMethod=app.PME, nonbondedCutoff=1*unit.nanometer, constraints=app.AllBonds) #check contraints Hbonds in the tutorial

thermostat = AndersenThermostat(298.15 * unit.kelvin, 1/unit.picosecond)
system.addForce(thermostat)
barostat = MonteCarloBarostat(1.013 * unit.bar, 298.15 * unit.kelvin)
system.addForce(barostat)
integrator = VerletIntegrator(0.002 * unit.picoseconds)
simulation = Simulation(prmtop.topology, system, integrator, platform)
if inpcrd.boxVectors is not None:
simulation.context.setPeriodicBoxVectors(*inpcrd.boxVectors)
simulation.context.setPositions(inpcrd.positions)

print("Minimizing")
simulation.minimizeEnergy(tolerance=10*unit.kilojoule/unit.mole)

print("NPT Equilibration")
simulation.reporters.append(StateDataReporter('{}_npt.cvs'.format(output_name), 500, step=True, potentialEnergy=True, temperature=True))
simulation.step(25000) #50 ps nvt equilibration
`

`
#------------------------

POPC bilayer

#------------------------
system = prmtop.createSystem(nonbondedMethod=app.PME, nonbondedCutoff=1*unit.nanometer, constraints=app.AllBonds) #check contraints Hbonds in the tutorial

thermostat = AndersenThermostat(298.15 * unit.kelvin, 1/unit.picosecond)
system.addForce(thermostat)
barostat = MonteCarloMembraneBarostat(1.013unit.bar, 0unit.bar*unit.nanometer, 298.15 * unit.kelvin, MonteCarloMembraneBarostat.XYIsotropic, MonteCarloMembraneBarostat.ZFree)
system.addForce(barostat)
integrator = VerletIntegrator(0.002 * unit.picoseconds)
simulation = Simulation(prmtop.topology, system, integrator, platform)
if prmtop.box_vectors is not None:
simulation.context.setPeriodicBoxVectors(*prmtop.box_vectors)
simulation.context.setPositions(prmtop.positions)

print("Minimizing")
simulation.minimizeEnergy(tolerance=10*unit.kilojoule/unit.mole)

print("NPT Equilibration")
simulation.reporters.append(StateDataReporter('{}_npt.cvs'.format(output_name), 500, step=True, potentialEnergy=True, temperature=True, totalEnergy=True, volume=True, density=True)) # !!!
simulation.reporters.append(NetCDFReporter('{}_npt.nc'.format(output_name), 500))
simulation.step(500000) #1 ns npt equilibration
`

[peastman]

Yes, this PR was merged fixing the bug. If you think you might have a problem, could you open a new issue for it rather than tacking it onto an unrelated one? Thanks!

[cespos]

Thanks. I created a new issue #2223