I'm looking to use NonbondedSoftcoreForce() to "ghost" in an atom (or collection of atoms) during a single MD simulation by changing the value of softcoreLJLambda from 0.0 to 1.0 over a specific amount of simulation time.
I'm using the Python interface and I think I have a clumsy, but correct, first draft. Would it be possible for someone to have a look over it to make sure that I've not missed something fundamental here. I'm worried that I'm misunderstood the concept of context here.
Thanks,
Mark
Code: Select all
# Example illustrating use of Free Energy OpenMM plugin to lambda map in an Atom
# Based upon argon-chemical-potential.py
from simtk.openmm.app import *
from simtk.openmm import *
from simtk.unit import *
from sys import stdout
import numpy
# =============================================================================
# Specify simulation parameters
# =============================================================================
nparticles = 216 # number of particles
mass = 39.9 * amu # mass
sigma = 3.4 * angstrom # Lennard-Jones sigma
epsilon = 0.238 * kilocalories_per_mole # Lennard-Jones well-depth
charge = 0.0 * elementary_charge # argon model has no charge
nequil_steps = 5000 # number of dynamics steps for equilibration
nprod_steps = 50000 # number of dynamics steps per production iteration
reduced_density = 0.85 # reduced density rho*sigma^3
temperature = 300 * kelvin # temperature
pressure = 1 * atmosphere # pressure
collision_rate = 5.0 / picosecond # collision rate for Langevin thermostat
timestep = 1 * femtosecond # integrator timestep
# =============================================================================
# Compute box size.
# =============================================================================
volume = nparticles*(sigma**3)/reduced_density
box_edge = volume**(1.0/3.0)
cutoff = min(box_edge*0.49, 2.5*sigma) # Compute cutoff
print "sigma = %s" % sigma
print "box_edge = %s" % box_edge
print "cutoff = %s" % cutoff
# =============================================================================
# Build system
# =============================================================================
# Create argon system where first particle is alchemically modified by lambda_value.
# It is currently zero
lambda_value = 0.0
system = System()
system.setDefaultPeriodicBoxVectors(Vec3(box_edge, 0, 0), Vec3(0, box_edge, 0), Vec3(0, 0, box_edge))
topology = Topology()
newChain = topology.addChain()
#nonbondedforce = NonbondedForce()
nonbondedforce = NonbondedSoftcoreForce()
nonbondedforce.setNonbondedMethod(NonbondedForce.CutoffPeriodic)
nonbondedforce.setCutoffDistance(cutoff)
for particle_index in range(nparticles):
system.addParticle(mass)
newResidue = topology.addResidue("UNK", newChain)
if (particle_index == 0 ):
# Add alchemically-modified particle.
# I think I'm using addAtom in the wrong way here...
topology.addAtom("MODD", "Ar", newResidue)
nonbondedforce.addParticle(charge, sigma, epsilon, lambda_value)
else:
# Add normal particle.
topology.addAtom("Argo", "ar", newResidue)
nonbondedforce.addParticle(charge, sigma, epsilon)
system.addForce(nonbondedforce)
print "System.getNumParticles() is %i" % system.getNumParticles()
print "system.getNumForces() is %i " % system.getNumForces()
# Create random initial positions.
import numpy.random
positions = Quantity(numpy.random.uniform(high=box_edge/angstroms, size=[nparticles,3]), angstrom)
# Create Integrator
integrator = LangevinIntegrator(temperature, collision_rate, timestep)
simulation = Simulation(topology, system, integrator)
# Initiate from Random set of positions
simulation.context.setPositions(positions)
# Minimize energy from coordinates.
simulation.minimizeEnergy()
print "simulation.system.getNumForces() is %i " % simulation.system.getNumForces()
# Equilibrate.
simulation.step(nequil_steps)
# Run dynamics.
simulation.reporters.append( PDBReporter('output.pdb', 100) )
simulation.reporters.append( StateDataReporter(stdout, 100, step=True, potentialEnergy=True, temperature=True) )
# =============================================================================
# Lambda map in the Argon atom
# =============================================================================
for i in range(10):
print i
simulation.context.getSystem().getForce(0).setParticleParameters(0, charge, sigma, epsilon, i*0.1)
print "Force Object status"
print simulation.context.getSystem().getForce(0).getParticleParameters(0)
simulation.step(1000)