Density problem with Amoeba

[Michael Schauperl]

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()

.
.
.

[Lee-Ping Wang]

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

[Michael Schauperl]

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

[Lee-Ping Wang]

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.

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

Thanks,

- Lee-Ping

[Peter Eastman]

There's also a multiple time step integrator in the code repository: http://wiki.simtk.org/openmm/VirtualRepository

Peter