MB-Pol spurious energy minima?

water14_no_local_energy_minimization.py

from simtk.openmm import app
import simtk.openmm as mm
from simtk import unit
import mbpol
import numpy as np
import mdtraj as md
import matplotlib.pyplot as plt

# create simulation
pdb = app.PDBFile("water14_cluster.pdb")
forcefield = app.ForceField(mbpol.__file__.replace('mbpol.py', 'mbpol.xml'))
system = forcefield.createSystem(pdb.topology, nonbondedMethod=app.CutoffNonPeriodic, nonbondedCutoff=1e3*unit.nanometer)
dt = 0.1*unit.femtoseconds
integrator = mm.VerletIntegrator(dt)
platform = mm.Platform.getPlatformByName('Reference')
simulation = app.Simulation(pdb.topology, system, integrator, platform)
simulation.context.setPositions(pdb.positions)
simulation.context.computeVirtualSites()

# convenient functions
def get_potential_energy():
    state = simulation.context.getState(getEnergy=True)
    potential_energy = state.getPotentialEnergy()
    return potential_energy / (unit.kilocalorie_per_mole)

def get_positions():
    return simulation.context.getState(getPositions=True).getPositions(asNumpy=True) / unit.nanometer

## minimize energy
#from simtk.openmm import LocalEnergyMinimizer
#LocalEnergyMinimizer.minimize(simulation.context)

# run until potential energy exceeds 0
simulation.context.setVelocitiesToTemperature(300*unit.kelvin)
traj = [get_positions()]
potential_energy_trace = [get_potential_energy()]
max_steps = 1000
while (potential_energy_trace[-1] < 0) and (len(traj) < max_steps):
    simulation.step(1)
    U = get_potential_energy()
    print('potential energy: {:.4} kcal/mol'.format(U))
    potential_energy_trace.append(get_potential_energy())
    traj.append(get_positions())

# save to pdb
top = md.load_pdb("water14_cluster.pdb").topology
xyz = np.array(traj)
traj = md.Trajectory(xyz, top)
traj.save_pdb("water14_trajectory_no_minimization.pdb")

# plot potential energy trace
y = np.array(potential_energy_trace[:-1])
x = (np.arange(len(y)) * dt) / unit.femtosecond
plt.plot(x, y)
plt.xlabel('time (fs)')
plt.ylabel('potential energy (kcal/mol)')
plt.title(r'velocity verlet ($\Delta t$={}fs)'.format(dt / unit.femtosecond))
plt.savefig('potential_energy_trace_no_minimization.png', dpi=300, bbox_inches='tight')

water14_with_local_energy_minimization.py

from simtk.openmm import app
import simtk.openmm as mm
from simtk import unit
import mbpol
import numpy as np
import mdtraj as md
import matplotlib.pyplot as plt

# create simulation
pdb = app.PDBFile("water14_cluster.pdb")
forcefield = app.ForceField(mbpol.__file__.replace('mbpol.py', 'mbpol.xml'))
system = forcefield.createSystem(pdb.topology, nonbondedMethod=app.CutoffNonPeriodic, nonbondedCutoff=1e3*unit.nanometer)
dt = 0.00001*unit.femtoseconds
integrator = mm.VerletIntegrator(dt)
platform = mm.Platform.getPlatformByName('Reference')
simulation = app.Simulation(pdb.topology, system, integrator, platform)
simulation.context.setPositions(pdb.positions)
simulation.context.computeVirtualSites()

# convenient functions
def get_potential_energy():
    state = simulation.context.getState(getEnergy=True)
    potential_energy = state.getPotentialEnergy()
    return potential_energy / (unit.kilocalorie_per_mole)

def get_positions():
    return simulation.context.getState(getPositions=True).getPositions(asNumpy=True) / unit.nanometer

# minimize energy
from simtk.openmm import LocalEnergyMinimizer
LocalEnergyMinimizer.minimize(simulation.context)

# run until potential energy exceeds 0
simulation.context.setVelocitiesToTemperature(300*unit.kelvin)
traj = [get_positions()]
potential_energy_trace = [get_potential_energy()]
max_steps = 1000
while (potential_energy_trace[-1] < 0) and (len(traj) < max_steps):
    simulation.step(1)
    U = get_potential_energy()
    print('potential energy: {:.4} kcal/mol'.format(U))
    potential_energy_trace.append(get_potential_energy())
    traj.append(get_positions())

# save to pdb
top = md.load_pdb("water14_cluster.pdb").topology
xyz = np.array(traj)
traj = md.Trajectory(xyz, top)
traj.save_pdb("water14_trajectory_with_minimization.pdb")

# plot potential energy trace
y = np.array(potential_energy_trace[:-1])
x = (np.arange(len(y)) * dt) / unit.femtosecond
plt.plot(x, y)
plt.xlabel('time (fs)')
plt.ylabel('potential energy (kcal/mol)')
plt.title(r'velocity verlet ($\Delta t$={}fs)'.format(dt / unit.femtosecond))
plt.savefig('potential_energy_trace_with_minimization.png', dpi=300, bbox_inches='tight')