Metadynamics implementation

[junlang-liu]

Hi,

I am trying to use OpenMM to run a metadynamics simulation, in which i would like to get the free energy during the dissemble of alpha-helix to random coil. Here is the code and error. Could you please help me with this? Thank you!

fname = 'protein'
pdb = PDBxFile(f'{fname}_amber14_tip3pfb_minimized.mmcif')
forcefield = ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')
platform = Platform.getPlatformByName('OpenCL')

# Prepare the system
system = forcefield.createSystem(pdb.topology, nonbondedMethod=PME, nonbondedCutoff=1*nanometer)

# Select the backbone atoms
backbone_atoms = [atom.index for atom in pdb.topology.atoms() if atom.name in ['N', 'CA', 'C']]
ref_positions = pdb.getPositions(asNumpy=False)

# Set up the integrator and simulation
integrator = LangevinMiddleIntegrator(363*kelvin, 1/picosecond, 0.002*picoseconds)
simulation = Simulation(pdb.topology, system, integrator, platform)
simulation.context.setPositions(pdb.positions)

# Create the metadynamics force
rmsd_force = RMSDForce(ref_positions, backbone_atoms)

# Run the metadynamics simulation
rmsd = BiasVariable(force=rmsd_force, minValue=0 * unit.angstrom, maxValue=10 * unit.angstrom, biasWidth=1 * unit.angstrom, periodic=False)
meta = Metadynamics(system=system, variables=[rmsd],
                    temperature=363.0, biasFactor=5.0,
                    height=1 * unit.kilocalories_per_mole,
                    frequency=50, saveFrequency=50, biasDir='./')

# set simulation reporters
simulation.reporters.append(DCDReporter('output.dcd', 10))
simulation.reporters.append(StateDataReporter(stdout, 10, step=True,
        potentialEnergy=True, temperature=True, progress=True, remainingTime=True, totalSteps=5000000, separator='\t'))

# Run small-scale simulation and plot the free energy landscape
meta.step(simulation, 500)
plot.imshow(meta.getFreeEnergy())
plot.show()

Here is the error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Input In [51], in <cell line: 27>()
     25 # Run the metadynamics simulation
     26 rmsd = BiasVariable(force=rmsd_force, minValue=0 * unit.angstrom, maxValue=10 * unit.angstrom, biasWidth=1 * unit.angstrom, periodic=False)
---> 27 meta = Metadynamics(system=system, variables=[rmsd],
     28                     temperature=363.0, biasFactor=5.0,
     29                     height=1 * unit.kilocalories_per_mole,
     30                     frequency=50, saveFrequency=50, biasDir='./')
     32 # set simulation reporters
     33 simulation.reporters.append(DCDReporter('output.dcd', 10))

File ~/opt/anaconda3/envs/research/lib/python3.9/site-packages/openmm/app/metadynamics.py:126, in Metadynamics.__init__(self, system, variables, temperature, biasFactor, height, frequency, saveFrequency, biasDir)
    124 self._totalBias = np.zeros(tuple(v.gridWidth for v in reversed(variables)))
    125 self._loadedBiases = {}
--> 126 self._syncWithDisk()
    127 self._deltaT = temperature*(biasFactor-1)
    128 varNames = ['cv%d' % i for i in range(len(variables))]

File ~/opt/anaconda3/envs/research/lib/python3.9/site-packages/openmm/app/metadynamics.py:269, in Metadynamics._syncWithDisk(self)
    267 self._totalBias = np.copy(self._selfBias)
    268 for bias in self._loadedBiases.values():
--> 269     self._totalBias += bias.bias

ValueError: operands could not be broadcast together with shapes (50,) (250,) (50,) 

[peastman]

As you run a metadynamics simulation, it saves biases to disk in biasDir. When you start a simulation, it loads any biases already present in that directory so that it can continue where it left off. In this case, it looks like you have files in that directory from an incompatible simulation that used different settings for the CVs.

[junlang-liu]

Thank you so much! It works!

Then I run into another dilemma:

If I run above code, it gives me the following error when adding the bias.

---------------------------------------------------------------------------
OpenMMException                           Traceback (most recent call last)
Input In [37], in <cell line: 53>()
     49 simulation.reporters.append(StateDataReporter(stdout, 10, step=True,
     50         potentialEnergy=True, temperature=True, progress=True, remainingTime=True, totalSteps=500, separator='\t'))
     52 # Run small-scale simulation and plot the free energy landscape
---> 53 meta.step(simulation, 500)
     54 plot.imshow(meta.getFreeEnergy())
     55 plot.show()

File ~/opt/anaconda3/envs/research/lib/python3.9/site-packages/openmm/app/metadynamics.py:174, in Metadynamics.step(self, simulation, steps)
    172 simulation.step(nextSteps)
    173 if simulation.currentStep % self.frequency == 0:
--> 174     position = self._force.getCollectiveVariableValues(simulation.context)
    175     energy = simulation.context.getState(getEnergy=True, groups={forceGroup}).getPotentialEnergy()
    176     height = self.height*np.exp(-energy/(unit.MOLAR_GAS_CONSTANT_R*self._deltaT))

File ~/opt/anaconda3/envs/research/lib/python3.9/site-packages/openmm/openmm.py:8969, in CustomCVForce.getCollectiveVariableValues(self, context)
   8954 def getCollectiveVariableValues(self, context):
   8955     r"""
   8956     getCollectiveVariableValues(self, context)
   8957     Get the current values of the collective variables in a Context.
   (...)
   8967         the values of the collective variables are computed and stored into this
   8968     """
-> 8969     return _openmm.CustomCVForce_getCollectiveVariableValues(self, context)

OpenMMException: getImplInContext: This Force is not present in the Context

Then I add system.addForce(rmsd_force)

Now it gives me the following error immediately

---------------------------------------------------------------------------
Exception                                 Traceback (most recent call last)
Input In [38], in <cell line: 42>()
     40 # Run the metadynamics simulation
     41 RMSD = BiasVariable(force=custom_cv_force, minValue=0 * unit.angstrom, maxValue=10 * unit.angstrom, biasWidth=1 * unit.angstrom, periodic=False)
---> 42 meta = Metadynamics(system=system, variables=[RMSD],
     43                     temperature=363.0, biasFactor=10.0,
     44                     height=1 * unit.kilocalories_per_mole,
     45                     frequency=50, saveFrequency=100, biasDir=f'{biasDir}')
     47 # set simulation reporters
     48 simulation.reporters.append(DCDReporter('output.dcd', 10))

File ~/opt/anaconda3/envs/research/lib/python3.9/site-packages/openmm/app/metadynamics.py:131, in Metadynamics.__init__(self, system, variables, temperature, biasFactor, height, frequency, saveFrequency, biasDir)
    129 self._force = mm.CustomCVForce('table(%s)' % ', '.join(varNames))
    130 for name, var in zip(varNames, variables):
--> 131     self._force.addCollectiveVariable(name, var.force)
    132 self._widths = [v.gridWidth for v in variables]
    133 self._limits = sum(([v.minValue, v.maxValue] for v in variables), [])

File ~/opt/anaconda3/envs/research/lib/python3.9/site-packages/openmm/openmm.py:8731, in CustomCVForce.addCollectiveVariable(self, name, variable)
   8728 if not variable.thisown:
   8729     s = ("the %s object does not own its corresponding OpenMM object"
   8730          % self.__class__.__name__)
-> 8731     raise Exception(s)
   8734 val = _openmm.CustomCVForce_addCollectiveVariable(self, name, variable)
   8736 variable.thisown=0

Exception: the CustomCVForce object does not own its corresponding OpenMM object

Don't know how to fix it as there seems no similar issues.

Thank you very much!

[peastman]

You need to create the Metadynamics object before the Simulation. In its constructor, it adds a new Force to the System. If the Simulation has already been created before that happens, it won't know about the new Force.

[pearlDingzhen]

have you found out how to fix the bug and run metadynamic?