Reproducibility of restarting metadynamics

[Miyazaki-Yu]

I am writing a program to resume metadynamics (MTD) in openmm-plumed after it has been interrupted. I am facing a problem where the potential immediately after resumption does not reproduce the potential at the time of interruption. The following code is the minimal code to reproduce that condition.

import openmmtools
from openmm import  Platform, VerletIntegrator, LangevinMiddleIntegrator
from openmm.app import CheckpointReporter, HBonds, Simulation
from openmm.app.pdbfile import PDBFile
from openmm.vec3 import Vec3
from sys import stdout
from openmmplumed import PlumedForce
from simtk.unit import (amu, angstroms, femtosecond, kelvin,
                        kilojoule_per_mole, picosecond)
import numpy as np
import shutil

ala2 = openmmtools.testsystems.AlanineDipeptideVacuum(constraints=HBonds)
temperature = 300 * kelvin
frictionCoeff = 1 / picosecond
timeStep = 2 * femtosecond


platform = Platform.getPlatformByName("Reference")   

script = """
RESTART NO 
UNITS LENGTH=A ENERGY=eV
psi: TORSION ATOMS=7,9,15,17
METAD ...
    LABEL=metad
    ARG=psi 
    PACE=100
    HEIGHT=0.1
    SIGMA=0.35
    FILE=HILLS_restart
    GRID_MIN=-pi
    GRID_MAX=pi
    GRID_SPACING=0.1
    BIASFACTOR=6.0
    TEMP=300.0
    CALC_RCT
... METAD
PRINT STRIDE=100 ARG=* FILE=COLVARS_restart
FLUSH STRIDE=100
"""

ala2.system.addForce(PlumedForce(script))
integrator = VerletIntegrator(timeStep)
simulation = Simulation(ala2.topology, ala2.system, integrator, platform)
simulation.context.setPositions(ala2.positions)
reporter_chk = CheckpointReporter(
    file="mtd.chk",
    reportInterval= 100,
    writeState=False
)
simulation.reporters.append(reporter_chk)
simulation.step(200)

state = simulation.context.getState(getPositions=True, getVelocities=True, getForces=True, getEnergy=True)
positions = state.getPositions(asNumpy=True)._value
velocities = state.getVelocities(asNumpy=True)._value
simulation.context.setPositions(state.getPositions())
forces = state.getForces(asNumpy=True)._value
kin = state.getKineticEnergy()._value
pot = state.getPotentialEnergy()._value

ala2_re = openmmtools.testsystems.AlanineDipeptideVacuum(constraints=HBonds)

platform_re = Platform.getPlatformByName("Reference")

script_re = """
RESTART
UNITS LENGTH=A ENERGY=eV
psi: TORSION ATOMS=7,9,15,17
METAD ...
    LABEL=metad
    ARG=psi 
    PACE=100
    HEIGHT=0.1
    SIGMA=0.35
    FILE=HILLS_restart
    GRID_MIN=-pi
    GRID_MAX=pi
    GRID_SPACING=0.1
    BIASFACTOR=6.0
    TEMP=300.0
    CALC_RCT
... METAD
PRINT STRIDE=100 ARG=* FILE=COLVARS_restart
FLUSH STRIDE=100
"""

ala2_re.system.addForce(PlumedForce(script_re))
integrator_re = VerletIntegrator(timeStep)
simulation_re = Simulation(ala2_re.topology, ala2_re.system, integrator_re, platform_re)
simulation_re.loadCheckpoint("mtd.chk")

state_re = simulation_re.context.getState(getPositions=True, getVelocities=True, getForces=True, getEnergy=True)
positions_re = state_re.getPositions(asNumpy=True)._value
velocities_re = state_re.getVelocities(asNumpy=True)._value
forces_re = state_re.getForces(asNumpy=True)._value
kin_re = state_re.getKineticEnergy()._value
pot_re = state_re.getPotentialEnergy()._value

print("---- error ----")
print("pos: ", np.abs(positions_re-positions).max())
print("vel: ", np.abs(velocities_re-velocities).max())
print("for: ", np.abs(forces_re-forces).max())
print("kin: ", kin_re-kin)
print("pot: ", pot_re-pot)

And the following outputs are obtained.

---- error ----
pos:  0.0
vel:  0.0
for:  0.0593004667935233
kin:  -1.65732980619282e-05
pot:  9.64877287559841

In the above code, I compare the values (positions, velocities, forces, kinetic energy, potential energy) immediately after advancing 200 steps of MTD-MD, with the values when that checkpoint is loaded. As a result, despite the positions and velocities matching exactly, the potentials and forces are significantly off.
issue.zip

Thanks!

[peastman]

This is going beyond my knowledge of Plumed, but I just skimmed the documentation at https://www.plumed.org/doc-v2.8/user-doc/html/_m_e_t_a_d.html, and it looks to me like you may need to add the GRID_RFILE option? Otherwise it doesn't know where to load the initial biases from.

[Miyazaki-Yu]

Thank you. As far as I've tried, it seems that HILLS and COLVARS will be loaded if they are set with the same name. I also tried the GRID_RFILE option, but there is still a discrepancy.

[Miyazaki-Yu]

To verify the conditions under which errors occur, I further conducted the verification by changing the conditions.

In the below program, I clarified the conditions under which discrepancies occur by changing the values of stride, reportInterval, and steps.

• There is no discrepancy from the time a checkpoint report occurs in openmm until a bias addition occurs in plumed next.

• Conversely, there is a discrepancy from the moment bias addition occurs in plumed until a checkpoint report occurs in openmm.

• If the timing of the bias addition in plumed and the checkpoint report in openmm overlaps, there will be a discrepancy until the next checkpoint report occurs.

I think this issue might be related to the previous one.

import openmmtools
from openmm import  Platform, VerletIntegrator, LangevinMiddleIntegrator
from openmm.app import CheckpointReporter, HBonds, Simulation
from openmm.app.pdbfile import PDBFile
from openmm.vec3 import Vec3
from sys import stdout
from openmmplumed import PlumedForce
from simtk.unit import (amu, angstroms, femtosecond, kelvin,
                        kilojoule_per_mole, picosecond)
import numpy as np
import shutil
import os
try:
    os.remove("COLVARS_restart")    
except OSError as e:
    pass
try:
    os.remove("HILLS_restart")
except OSError as e:
    pass
try:
    os.remove("mtd.chk")
except OSError as e:
    pass


ala2 = openmmtools.testsystems.AlanineDipeptideVacuum(constraints=HBonds)
temperature = 300 * kelvin
frictionCoeff = 1 / picosecond
timeStep = 2 * femtosecond


platform = Platform.getPlatformByName("Reference")   
stride = 3; reportInterval = 2; steps = 8

script = f"""
RESTART NO 
UNITS LENGTH=A ENERGY=eV
psi: TORSION ATOMS=7,9,15,17
METAD ...
    LABEL=metad
    ARG=psi 
    PACE={stride}
    HEIGHT=0.1
    SIGMA=0.35
    FILE=HILLS_restart
    GRID_MIN=-pi
    GRID_MAX=pi
    GRID_SPACING=0.1
    BIASFACTOR=6.0
    TEMP=300.0
    CALC_RCT
... METAD
PRINT STRIDE={stride} ARG=* FILE=COLVARS_restart
FLUSH STRIDE={stride}
"""

ala2.system.addForce(PlumedForce(script))
integrator = VerletIntegrator(timeStep)
simulation = Simulation(ala2.topology, ala2.system, integrator, platform)
simulation.context.setPositions(ala2.positions)
reporter_chk = CheckpointReporter(
    file="mtd.chk",
    reportInterval= reportInterval,
    writeState=False
)
simulation.reporters.append(reporter_chk)
simulation.step(steps)

state = simulation.context.getState(getPositions=True, getVelocities=True, getForces=True, getEnergy=True)
positions = state.getPositions(asNumpy=True)._value
velocities = state.getVelocities(asNumpy=True)._value
simulation.context.setPositions(state.getPositions())
forces = state.getForces(asNumpy=True)._value
kin = state.getKineticEnergy()._value
pot = state.getPotentialEnergy()._value

ala2_re = openmmtools.testsystems.AlanineDipeptideVacuum(constraints=HBonds)

platform_re = Platform.getPlatformByName("Reference")

script_re = f"""
RESTART
UNITS LENGTH=A ENERGY=eV
psi: TORSION ATOMS=7,9,15,17
METAD ...
    LABEL=metad
    ARG=psi 
    PACE={stride}
    HEIGHT=0.1
    SIGMA=0.35
    FILE=HILLS_restart
    GRID_MIN=-pi
    GRID_MAX=pi
    GRID_SPACING=0.1
    BIASFACTOR=6.0
    TEMP=300.0
    CALC_RCT
... METAD
PRINT STRIDE={stride} ARG=* FILE=COLVARS_restart
FLUSH STRIDE={stride}
"""

ala2_re.system.addForce(PlumedForce(script_re))
integrator_re = VerletIntegrator(timeStep)
simulation_re = Simulation(ala2_re.topology, ala2_re.system, integrator_re, platform_re)
simulation_re.loadCheckpoint("mtd.chk")
steps_remain = steps % reportInterval
if steps_remain > 0:
    print(f"advance {steps_remain} MD step")
    simulation_re.step(steps_remain)
else:
    print("no MD step")
state_re = simulation_re.context.getState(getPositions=True, getVelocities=True, getForces=True, getEnergy=True)
positions_re = state_re.getPositions(asNumpy=True)._value
velocities_re = state_re.getVelocities(asNumpy=True)._value
forces_re = state_re.getForces(asNumpy=True)._value
kin_re = state_re.getKineticEnergy()._value
pot_re = state_re.getPotentialEnergy()._value
print("---- allclose ----")
print("pos: ", np.allclose(positions_re,positions))
print("vel: ", np.allclose(velocities_re,velocities))
print("for: ", np.allclose(forces_re,forces))
print("---- error ----")
print("pos: ", np.abs(positions_re-positions).max())
print("vel: ", np.abs(velocities_re-velocities).max())
print("for: ", np.abs(forces_re-forces).max())
print("kin: ", kin_re-kin)
print("pot: ", pot_re-pot)
print("steps_remain",steps_remain)