Density problem with Amoeba
Posted: Tue Jan 28, 2014 8:11 am
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()
.
.
.
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()
.
.
.