Hi,
I got a problem by simulating peptides and proteins with the Amoeba FF.
If i start with a minimizied structure (density close to 1) and then heat my system up with a ramp the density decrease to a value of 0.7 and after a while reach a value of 0.95 to 0.98. I am quite sure that this density is not correct.
Can anyone tell me what is happening or what kind of protocol you use for Amoeba to equilibrate the system?
Thanks,
Michael
Probably you can have a look at my input file:
from __future__ import print_function
from simtk.openmm.app import *
from simtk.openmm import *
from simtk.unit import *
from sys import stdout
import os
import time
import datetime
import sys
from progressreporter import ProgressReporter
import pickle
Platform.loadPluginsFromDirectory(Platform.getDefaultPluginsDirectory())
pdb=PDBFile('peptide_wat.pdb')
forcefield = ForceField('amoeba2009.xml')
platform = Platform.getPlatformByName('CUDA')
properties = {'CudaPrecision': 'mixed'}
system = forcefield.createSystem(pdb.topology, nonbondedMethod=PME, constraints=None, rigidWater=False, polarization='mutual',mutualInducedTargetEpsilon=0.005, vdwCutoff=0.9*nanometer,nonbondedCutoff=0.7*nanometer)
system.addForce(MonteCarloBarostat(1*bar,300*kelvin))
integrator = LangevinIntegrator(300*kelvin, 20/picosecond, 1*femtoseconds)
simulation = Simulation(pdb.topology, system, integrator,platform, properties)
simulation.context.setPositions(pdb.positions)
simulation.reporters.append(StateDataReporter('e0.dat', 2000, step=True, time=True, potentialEnergy=True, kineticEnergy=True, totalEnergy=True, temperature=True, volume=True, density=True))
simulation.reporters.append(DCDReporter('e0.dcd', 2000))
multipole = [f for f in system.getForces() if isinstance(f, AmoebaMultipoleForce)][0]
f1=open('./multipoles.dat', 'w')
#Equilibration
print('Minimizing')
simulation.minimizeEnergy()
print('Heating Up')
for i in range(300):
integrator.setTemperature(i*kelvin)
simulation.step(20000)
print ('Equilibrating')
for i in range(100):
integrator.setTemperature(300*kelvin)
simulation.step(100000)
print(multipole.getInducedDipoles(simulation.context),file=f1)
f1.close()
.
.
.
Density problem with Amoeba
- Lee-Ping Wang
- Posts: 102
- Joined: Sun Jun 19, 2011 5:14 pm
Re: Density problem with Amoeba
Hi Michael,
Thanks for your question. I can see a few things in your script that you can modify to troubleshoot the problem:
1) Try a time step of 0.5 fs. Due to the flexible bonds, a 1.0 fs simulation corresponds to artificial heating of about 10 degrees.
2) Try setting mutualInducedTargetEpsilon to 1e-5. The current criterion is a bit too loose and leads to eneregy conservation issues.
I know this might slow down your simulation, but let me know if this helps. We can speed things back up later with a multiple timestep integrator.
Thanks,
- Lee-Ping
Thanks for your question. I can see a few things in your script that you can modify to troubleshoot the problem:
1) Try a time step of 0.5 fs. Due to the flexible bonds, a 1.0 fs simulation corresponds to artificial heating of about 10 degrees.
2) Try setting mutualInducedTargetEpsilon to 1e-5. The current criterion is a bit too loose and leads to eneregy conservation issues.
I know this might slow down your simulation, but let me know if this helps. We can speed things back up later with a multiple timestep integrator.
Thanks,
- Lee-Ping
- Michael Schauperl
- Posts: 22
- Joined: Wed Sep 04, 2013 6:10 am
Re: Density problem with Amoeba
Hi Lee-Ping,
Thanks for your quick response. I started simulations with the new parameters and at first glance it seems to work but it is nearly a magnitude slower.. I will report you in a few hours if the system is still stabile.
Thanks in advance.
Michael
Thanks for your quick response. I started simulations with the new parameters and at first glance it seems to work but it is nearly a magnitude slower.. I will report you in a few hours if the system is still stabile.
Thanks in advance.
Michael
- Lee-Ping Wang
- Posts: 102
- Joined: Sun Jun 19, 2011 5:14 pm
Re: Density problem with Amoeba
Hi Michael,
It's good to hear that things might be working now. I understand it's a lot slower now but at least it's not running in double precision. Here's a quote from Jay Ponder discussing convergence criteria:
>> A remaining decision is the level of convergence we should adopt as a
>> standard default value. I already have evidence 0.000001 RMS is good
>> enough (ie, negligible bias) for everything I will want to do.
>> The 0.00001 level may also be fine-- especially since the RMS change
>> metric used by AMOEBA is overly conservative at tight convergence.
>> Iterative solvers tend to oscillate near the solution. So at an RMS
>> of 0.00001 between iterations, the error from the exact solution is
>> often much closer to 0.000001 RMS.
We should look at implementing improved SCF convergence algorithms for the dipoles, that might give you a factor of 2 or so. Jay Ponder already has this in TINKER so I'll see if we can put this in. (OpenMM is still a fair bit faster than TINKER on multiple cores, even though the latter has the improved convergence algorithm.)
Also, here is a good multiple timestep Langevin integrator (code from John Chodera) which should allow you to go to 2.0 fs time steps. Let me know if it helps.
Thanks,
- Lee-Ping
It's good to hear that things might be working now. I understand it's a lot slower now but at least it's not running in double precision. Here's a quote from Jay Ponder discussing convergence criteria:
>> A remaining decision is the level of convergence we should adopt as a
>> standard default value. I already have evidence 0.000001 RMS is good
>> enough (ie, negligible bias) for everything I will want to do.
>> The 0.00001 level may also be fine-- especially since the RMS change
>> metric used by AMOEBA is overly conservative at tight convergence.
>> Iterative solvers tend to oscillate near the solution. So at an RMS
>> of 0.00001 between iterations, the error from the exact solution is
>> often much closer to 0.000001 RMS.
We should look at implementing improved SCF convergence algorithms for the dipoles, that might give you a factor of 2 or so. Jay Ponder already has this in TINKER so I'll see if we can put this in. (OpenMM is still a fair bit faster than TINKER on multiple cores, even though the latter has the improved convergence algorithm.)
Also, here is a good multiple timestep Langevin integrator (code from John Chodera) which should allow you to go to 2.0 fs time steps. Let me know if it helps.
Code: Select all
def MTSVVVRIntegrator(temperature, collision_rate, timestep, system, ninnersteps=4):
"""
Create a multiple timestep velocity verlet with velocity randomization (VVVR) integrator.
ARGUMENTS
temperature (Quantity compatible with kelvin) - the temperature
collision_rate (Quantity compatible with 1/picoseconds) - the collision rate
timestep (Quantity compatible with femtoseconds) - the integration timestep
system (simtk.openmm.System) - system whose forces will be partitioned
ninnersteps (int) - number of inner timesteps (default: 4)
RETURNS
integrator (openmm.CustomIntegrator) - a VVVR integrator
NOTES
This integrator is equivalent to a Langevin integrator in the velocity Verlet discretization with a
timestep correction to ensure that the field-free diffusion constant is timestep invariant. The inner
velocity Verlet discretization is transformed into a multiple timestep algorithm.
REFERENCES
VVVR Langevin integrator:
* http://arxiv.org/abs/1301.3800
* http://arxiv.org/abs/1107.2967 (to appear in PRX 2013)
TODO
Move initialization of 'sigma' to setting the per-particle variables.
"""
# Multiple timestep Langevin integrator.
for i in system.getForces():
if i.__class__.__name__ in ["NonbondedForce", "CustomNonbondedForce", "AmoebaVdwForce", "AmoebaMultipoleForce"]:
# Slow force.
logger.info(i.__class__.__name__ + " is a Slow Force\n")
i.setForceGroup(1)
else:
logger.info(i.__class__.__name__ + " is a Fast Force\n")
# Fast force.
i.setForceGroup(0)
kB = BOLTZMANN_CONSTANT_kB * AVOGADRO_CONSTANT_NA
kT = kB * temperature
integrator = CustomIntegrator(timestep)
integrator.addGlobalVariable("dt_fast", timestep/float(ninnersteps)) # fast inner timestep
integrator.addGlobalVariable("kT", kT) # thermal energy
integrator.addGlobalVariable("a", np.exp(-collision_rate*timestep)) # velocity mixing parameter
integrator.addGlobalVariable("b", np.sqrt((2/(collision_rate*timestep)) * np.tanh(collision_rate*timestep/2))) # timestep correction parameter
integrator.addPerDofVariable("sigma", 0)
integrator.addPerDofVariable("x1", 0) # position before application of constraints
#
# Pre-computation.
# This only needs to be done once, but it needs to be done for each degree of freedom.
# Could move this to initialization?
#
integrator.addComputePerDof("sigma", "sqrt(kT/m)")
#
# Velocity perturbation.
#
integrator.addComputePerDof("v", "sqrt(a)*v + sqrt(1-a)*sigma*gaussian")
integrator.addConstrainVelocities();
#
# Symplectic inner multiple timestep.
#
integrator.addUpdateContextState();
integrator.addComputePerDof("v", "v + 0.5*b*dt*f1/m")
for innerstep in range(ninnersteps):
# Fast inner symplectic timestep.
integrator.addComputePerDof("v", "v + 0.5*b*dt_fast*f0/m")
integrator.addComputePerDof("x", "x + v*b*dt_fast")
integrator.addComputePerDof("x1", "x")
integrator.addConstrainPositions();
integrator.addComputePerDof("v", "v + 0.5*b*dt_fast*f0/m + (x-x1)/dt_fast")
integrator.addComputePerDof("v", "v + 0.5*b*dt*f1/m") # TODO: Additional velocity constraint correction?
integrator.addConstrainVelocities();
#
# Velocity randomization
#
integrator.addComputePerDof("v", "sqrt(a)*v + sqrt(1-a)*sigma*gaussian")
integrator.addConstrainVelocities();
return integrator
- Lee-Ping
- Peter Eastman
- Posts: 2573
- Joined: Thu Aug 09, 2007 1:25 pm
Re: Density problem with Amoeba
There's also a multiple time step integrator in the code repository: http://wiki.simtk.org/openmm/VirtualRepository
Peter
Peter