Skip to content

Instantly share code, notes, and snippets.

@synapticarbors

synapticarbors/gen_bstate.py

Created November 4, 2014 02:06
Show Gist options
  • Save synapticarbors/2cfceba28d4674d81989 to your computer and use it in GitHub Desktop.
Save synapticarbors/2cfceba28d4674d81989 to your computer and use it in GitHub Desktop.
Download ZIP
NetCDFReporter problem
Raw
gen_bstate.py
import os
import numpy as np
import argparse
import simtk
import simtk.unit as units
import simtk.openmm as openmm
import wcadimer
def minimize(platform, system, positions):
# Create a Context.
timestep = 1.0 * units.femtoseconds
integrator = openmm.VerletIntegrator(timestep)
context = openmm.Context(system, integrator, platform)
# Set coordinates.
context.setPositions(positions)
# Compute initial energy.
state = context.getState(getEnergy=True)
initial_potential = state.getPotentialEnergy()
print "initial potential: %12.3f kcal/mol" % (initial_potential / units.kilocalories_per_mole)
# Minimize.
openmm.LocalEnergyMinimizer.minimize(context)
# Compute final energy.
state = context.getState(getEnergy=True, getPositions=True)
final_potential = state.getPotentialEnergy()
positions = state.getPositions(asNumpy=True)
# Report
print "final potential : %12.3f kcal/mol" % (final_potential / units.kilocalories_per_mole)
return positions
def norm(n01):
return n01.unit * np.sqrt(np.dot(n01/n01.unit, n01/n01.unit))
def run(platform_name, deviceid, two):
not_obs = [True, True]
system, coordinates = wcadimer.WCADimer()
print "Time step: ", (wcadimer.stable_timestep * 2.0).in_units_of(units.femtoseconds)
# Minimization
platform = openmm.Platform.getPlatformByName('Reference')
print 'Minimizing energy...'
coordinates = minimize(platform, system, coordinates)
print 'Separation distance: {}'.format(norm(coordinates[1,:] - coordinates[0,:]) / units.angstroms)
print 'Equilibrating...'
platform = openmm.Platform.getPlatformByName(platform_name)
if platform_name == 'CUDA':
platform.setPropertyDefaultValue('CudaDeviceIndex', '%d' % deviceid)
platform.setPropertyDefaultValue('CudaPrecision', 'mixed')
elif platform_name == 'OpenCL':
platform.setPropertyDefaultValue('OpenCLDeviceIndex', '%d' % deviceid)
platform.setPropertyDefaultValue('OpenCLPrecision', 'mixed')
integrator = openmm.LangevinIntegrator(wcadimer.temperature,
wcadimer.collision_rate,
2.0 * wcadimer.stable_timestep)
context = openmm.Context(system, integrator, platform)
context.setPositions(coordinates)
context.setVelocitiesToTemperature(wcadimer.temperature)
if two:
while not_obs[0] or not_obs[1]:
integrator.step(5000)
state = context.getState(getPositions=True)
coordinates = state.getPositions(asNumpy=True)
sep_dist = norm(coordinates[1,:] - coordinates[0,:]) / units.angstroms
print 'Separation distance: {}'.format(sep_dist)
if sep_dist < 5.7:
not_obs[0] = False
tag = '_a'
sep_dist_a = sep_dist
else:
not_obs[1] = False
tag = '_b'
sep_dist_b = sep_dist
if not os.path.isdir('bstates'):
os.makedirs('bstates')
np.save('bstates/init_coords{}.npy'.format(tag), coordinates / units.nanometers)
print sep_dist_a, sep_dist_b
else:
integrator.step(5000)
state = context.getState(getPositions=True)
coordinates = state.getPositions(asNumpy=True)
print 'Separation distance: {}'.format(norm(coordinates[1,:] - coordinates[0,:]) / units.angstroms)
if not os.path.isdir('bstates'):
os.makedirs('bstates')
np.save('bstates/init_coords.npy', coordinates / units.nanometers)
if __name__ == '__main__':
available_platforms = [openmm.Platform.getPlatform(i).getName() for i in range(openmm.Platform.getNumPlatforms())]
print 'Available platforms:'
for platform_index, platform in enumerate(available_platforms):
print "%5d %s" % (platform_index, platform)
print
parser = argparse.ArgumentParser()
parser.add_argument('-p', '--platform', choices=available_platforms,
default='CUDA',
help='Platform to use for equilibration')
parser.add_argument('-d', '--deviceid', type=int, default=0, help='Device id for equilibration if using a GPU')
parser.add_argument('--two', action="store_true", default=False, help='Create initial states in both basins')
args = parser.parse_args()
run(args.platform, args.deviceid, args.two)
Raw
integrator.xml
<?xml version="1.0" ?>
<Integrator constraintTolerance="1e-05" friction=".46509600382126837" randomSeed="-1990315844" stepSize=".004300187452843778" temperature="98.88" type="LangevinIntegrator" version="1"/>
Raw
openmm_app.py
import numpy as np
from simtk.openmm import app
from simtk.openmm import Vec3
import simtk.openmm.openmm as openmm
import simtk.unit as units
import wcadimer
from mdtraj.reporters import NetCDFReporter, DCDReporter, HDF5Reporter
import random
import argparse
from sys import stdout
def run(opts):
system_xml_file = opts.system
integrator_xml_file = opts.integrator
coords_f = opts.coords
velocs_f = opts.velocs
restart = opts.restart
platform_name = opts.platform
deviceid = opts.device
write_freq = opts.write_freq
output = opts.output
nsteps = opts.nsteps
platform_properties = {'OpenCLPrecision': 'mixed',
'OpenCLPlatformIndex': '0',
'OpenCLDeviceIndex': '0',
'CudaPrecision': 'mixed',
'CudaDeviceIndex': '0'}
platform_properties['CudaDeviceIndex'] = deviceid
platform_properties['OpenCLDeviceIndex'] = deviceid
with open(system_xml_file, 'r') as f:
system = openmm.XmlSerializer.deserialize(f.read())
with open(integrator_xml_file, 'r') as f:
integrator = openmm.XmlSerializer.deserialize(f.read())
integrator.setRandomNumberSeed(random.randint(0, 2**16))
platform = openmm.Platform.getPlatformByName(platform_name)
properties = {key: platform_properties[key] for key in platform_properties if key.lower().startswith(platform_name.lower())}
print properties
# Create dummy topology to satisfy Simulation object
topology = app.Topology()
volume = wcadimer.natoms / wcadimer.density
length = volume**(1.0/3.0)
L = length.value_in_unit(units.nanometer)
topology.setUnitCellDimensions(Vec3(L, L, L)*units.nanometer)
simulation = app.Simulation(topology, system,
integrator, platform, properties)
coords = units.Quantity(np.load(coords_f), units.nanometer)
simulation.context.setPositions(coords)
if restart:
velocs = units.Quantity(np.load(velocs_f), units.nanometer / units.picosecond)
simulation.context.setVelocities(velocs)
else:
simulation.context.setVelocitiesToTemperature(wcadimer.temperature)
# Attach reporters
simulation.reporters.append(NetCDFReporter(output + '.nc', write_freq, atomSubset=[0,1]))
#simulation.reporters.append(DCDReporter(output + '.dcd', write_freq, atomSubset=[0,1]))
#simulation.reporters.append(HDF5Reporter(output + '.h5', write_freq, atomSubset=[0,1]))
simulation.reporters.append(app.StateDataReporter(stdout, 20*write_freq, step=True,
potentialEnergy=True, temperature=True, progress=True, remainingTime=True,
speed=True, totalSteps=nsteps, separator='\t'))
# Run segment
simulation.step(nsteps)
# Write restart data
state = simulation.context.getState(getPositions=True, getVelocities=True)
coords = state.getPositions(asNumpy=True)
velocs = state.getVelocities(asNumpy=True)
np.save(output + '_coor.npy', coords)
np.save(output + '_veloc.npy', velocs)
if __name__ == '__main__':
parser = argparse.ArgumentParser()
parser.add_argument('-c', '--coords', required=True, help='Coordinates to initialize simulation from')
parser.add_argument('-v', '--velocs', help='Velocities to initialize simulation from')
parser.add_argument('-r', '--restart', action='store_true')
parser.add_argument('-s', '--system', required=True, help='OpenMM system xml file')
parser.add_argument('-i', '--integrator', required=True, help='OpenMM integrator xml file')
parser.add_argument('-p', '--platform', default='CUDA',
choices=['CUDA', 'OpenCL', 'Reference', 'CPU'], help='OpenMM platform name')
parser.add_argument('-d', '--device', default='0', help='Device ID')
parser.add_argument('-w', '--write_freq', required=True, type=int, help='write frequency for output files')
parser.add_argument('-o', '--output', required=True, help='base name for output files')
parser.add_argument('-n', '--nsteps', type=int, required=True, help='Number of steps to run')
opts = parser.parse_args()
if opts.restart and opts.velocs is None:
parser.error('Most supply velocities if restart flag present')
run(opts)
Raw
system.xml
<?xml version="1.0" ?>
<System type="System" version="1">
<PeriodicBoxVectors>
<A x="2.067948818206787" y="0" z="0"/>
<B x="0" y="2.067948818206787" z="0"/>
<C x="0" y="0" z="2.067948818206787"/>
</PeriodicBoxVectors>
<Particles>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
<Particle mass="39.9"/>
</Particles>
<Constraints/>
<Forces>
<Force cutoff=".38163709642518684" energy="4.0*epsilon*((sigma/r)^12 - (sigma/r)^6) + epsilon" forceGroup="0" method="2" switchingDistance="-1" type="CustomNonbondedForce" useLongRangeCorrection="0" useSwitchingFunction="0" version="1">
<PerParticleParameters/>
<GlobalParameters>
<Parameter default=".9977366965464258" name="epsilon"/>
<Parameter default=".34" name="sigma"/>
</GlobalParameters>
<Particles>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
<Particle/>
</Particles>
<Exclusions>
<Exclusion p1="0" p2="1"/>
</Exclusions>
<Functions/>
<InteractionGroups/>
</Force>
<Force energy="h*(1-((r-r0-w)/w)^2)^2;" forceGroup="0" type="CustomBondForce" version="1">
<PerBondParameters/>
<GlobalParameters>
<Parameter default="4.110675189771274" name="h"/>
<Parameter default=".38163709642518684" name="r0"/>
<Parameter default=".19081854821259342" name="w"/>
</GlobalParameters>
<Bonds>
<Bond p1="0" p2="1"/>
</Bonds>
</Force>
</Forces>
</System>
Raw
wcadimer.py
#!/usr/local/bin/env python
#=============================================================================================
# MODULE DOCSTRING
#=============================================================================================
"""
WCA fluid and WCA dimer systems.
DESCRIPTION
COPYRIGHT
@author John D. Chodera <jchodera@gmail.com>
All code in this repository is released under the GNU General Public License.
This program is free software: you can redistribute it and/or modify it under
the terms of the GNU General Public License as published by the Free Software
Foundation, either version 3 of the License, or (at your option) any later
version.
This program is distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
PARTICULAR PURPOSE. See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with
this program. If not, see <http://www.gnu.org/licenses/>.
TODO
"""
#=============================================================================================
# GLOBAL IMPORTS
#=============================================================================================
import numpy
import simtk
import simtk.unit as units
import simtk.openmm as openmm
#=============================================================================================
# CONSTANTS
#=============================================================================================
kB = units.BOLTZMANN_CONSTANT_kB * units.AVOGADRO_CONSTANT_NA
#=============================================================================================
# DEFAULT PARAMETERS
#=============================================================================================
natoms = 216 # number of particles
# WCA fluid parameters (argon).
mass = 39.9 * units.amu # reference mass
sigma = 3.4 * units.angstrom # reference lengthscale
epsilon = 120.0 * units.kelvin * kB # reference energy
r_WCA = 2.**(1./6.) * sigma # interaction truncation range
tau = units.sqrt(sigma**2 * mass / epsilon) # characteristic timescale
# Simulation conditions.
temperature = 0.824 / (kB / epsilon) # temperature
kT = kB * temperature # thermal energy
beta = 1.0 / kT # inverse temperature
density = 0.96 / sigma**3 # default density
stable_timestep = 0.001 * tau # stable timestep
collision_rate = 1 / tau # collision rate for Langevin interator
# Dimer potential parameters.
h = 5.0 * kT # barrier height
r0 = r_WCA # compact state distance
w = 0.5 * r_WCA # extended state distance is r0 + 2*w
#=============================================================================================
# WCA Fluid
#=============================================================================================
def WCAFluid(N=natoms, density=density, mm=None, mass=mass, epsilon=epsilon, sigma=sigma):
"""
Create a Weeks-Chandler-Andersen system.
OPTIONAL ARGUMENTS
N (int) - total number of atoms (default: 150)
density (float) - N sigma^3 / V (default: 0.96)
sigma
epsilon
"""
# Choose OpenMM package.
if mm is None:
mm = openmm
# Create system
system = mm.System()
# Compute total system volume.
volume = N / density
# Make system cubic in dimension.
length = volume**(1.0/3.0)
# TODO: Can we change this to use tuples or 3x3 array?
a = units.Quantity(numpy.array([1.0, 0.0, 0.0], numpy.float32), units.nanometer) * length/units.nanometer
b = units.Quantity(numpy.array([0.0, 1.0, 0.0], numpy.float32), units.nanometer) * length/units.nanometer
c = units.Quantity(numpy.array([0.0, 0.0, 1.0], numpy.float32), units.nanometer) * length/units.nanometer
print "box edge length = %s" % str(length)
system.setDefaultPeriodicBoxVectors(a, b, c)
# Add particles to system.
for n in range(N):
system.addParticle(mass)
# Create nonbonded force term implementing Kob-Andersen two-component Lennard-Jones interaction.
energy_expression = '4.0*epsilon*((sigma/r)^12 - (sigma/r)^6) + epsilon'
# Create force.
force = mm.CustomNonbondedForce(energy_expression)
# Set epsilon and sigma global parameters.
force.addGlobalParameter('epsilon', epsilon)
force.addGlobalParameter('sigma', sigma)
# Add particles
for n in range(N):
force.addParticle([])
# Set periodic boundary conditions with cutoff.
force.setNonbondedMethod(mm.CustomNonbondedForce.CutoffNonPeriodic)
print "setting cutoff distance to %s" % str(r_WCA)
force.setCutoffDistance(r_WCA)
# Add nonbonded force term to the system.
system.addForce(force)
# Create initial coordinates using random positions.
coordinates = units.Quantity(numpy.random.rand(N,3), units.nanometer) * (length / units.nanometer)
# Return system and coordinates.
return [system, coordinates]
#=============================================================================================
# WCA dimer plus fluid
#=============================================================================================
def WCADimer(N=natoms, density=density, mm=None, mass=mass, epsilon=epsilon, sigma=sigma, h=h, r0=r0, w=w):
"""
Create a bistable bonded pair of particles (indices 0 and 1) optionally surrounded by a Weeks-Chandler-Andersen fluid.
The bistable potential has form
U(r) = h*(1-((r-r0-w)/w)^2)^2
where r0 is the compact state separation, r0+2w is the extended state separation, and h is the barrier height.
The WCA potential has form
U(r) = 4 epsilon [ (sigma/r)^12 - (sigma/r)^6 ] + epsilon (r < r*)
= 0 (r >= r*)
where r* = 2^(1/6) sigma.
OPTIONAL ARGUMENTS
N (int) - total number of atoms (default: 2)
density (float) - number density of particles (default: 0.96 / sigma**3)
mass (simtk.unit.Quantity of mass) - particle mass (default: 39.948 amu)
sigma (simtk.unit.Quantity of length) - Lennard-Jones sigma parameter (default: 0.3405 nm)
epsilon (simtk.unit.Quantity of energy) - Lennard-Jones well depth (default: (119.8 Kelvin)*kB)
h (simtk.unit.Quantity of energy) - bistable potential barrier height (default: ???)
r0 (simtk.unit.Quantity of length) - bistable potential compact state separation (default: ???)
w (simtk.unit.Quantity of length) - bistable potential extended state separation is r0+2*w (default: ???)
"""
# Choose OpenMM package.
if mm is None:
mm = openmm
# Compute cutoff for WCA fluid.
r_WCA = 2.**(1./6.) * sigma # cutoff at minimum of potential
# Create system
system = mm.System()
# Compute total system volume.
volume = N / density
# Make system cubic in dimension.
length = volume**(1.0/3.0)
a = units.Quantity(numpy.array([1.0, 0.0, 0.0], numpy.float32), units.nanometer) * length/units.nanometer
b = units.Quantity(numpy.array([0.0, 1.0, 0.0], numpy.float32), units.nanometer) * length/units.nanometer
c = units.Quantity(numpy.array([0.0, 0.0, 1.0], numpy.float32), units.nanometer) * length/units.nanometer
print "box edge length = %s" % str(length)
system.setDefaultPeriodicBoxVectors(a, b, c)
# Add particles to system.
for n in range(N):
system.addParticle(mass)
# WCA: Lennard-Jones truncated at minim and shifted so potential is zero at cutoff.
energy_expression = '4.0*epsilon*((sigma/r)^12 - (sigma/r)^6) + epsilon'
# Create force.
force = mm.CustomNonbondedForce(energy_expression)
# Set epsilon and sigma global parameters.
force.addGlobalParameter('epsilon', epsilon)
force.addGlobalParameter('sigma', sigma)
# Add particles
for n in range(N):
force.addParticle([])
# Add exclusion between bonded particles.
force.addExclusion(0,1)
# Set periodic boundary conditions with cutoff.
if (N > 2):
force.setNonbondedMethod(mm.CustomNonbondedForce.CutoffPeriodic)
else:
force.setNonbondedMethod(mm.CustomNonbondedForce.CutoffNonPeriodic)
print "setting cutoff distance to %s" % str(r_WCA)
force.setCutoffDistance(r_WCA)
# Add nonbonded force term to the system.
system.addForce(force)
# Add dimer potential to first two particles.
dimer_force = openmm.CustomBondForce('h*(1-((r-r0-w)/w)^2)^2;')
dimer_force.addGlobalParameter('h', h) # barrier height
dimer_force.addGlobalParameter('r0', r0) # compact state separation
dimer_force.addGlobalParameter('w', w) # second minimum is at r0 + 2*w
dimer_force.addBond(0, 1, [])
system.addForce(dimer_force)
# Create initial coordinates using random positions.
coordinates = units.Quantity(numpy.random.rand(N,3), units.nanometer) * (length / units.nanometer)
# Reposition dimer particles at compact minimum.
coordinates[0,:] *= 0.0
coordinates[1,:] *= 0.0
coordinates[1,0] = r0
# Return system and coordinates.
return [system, coordinates]
#=============================================================================================
# WCA dimer in vacuum
#=============================================================================================
def WCADimerVacuum(mm=None, mass=mass, epsilon=epsilon, sigma=sigma, h=h, r0=r0, w=w):
"""
Create a bistable dimer.
OPTIONAL ARGUMENTS
"""
# Choose OpenMM package.
if mm is None:
mm = openmm
# Create system
system = mm.System()
# Add particles to system.
for n in range(2):
system.addParticle(mass)
# Add dimer potential to first two particles.
dimer_force = openmm.CustomBondForce('h*(1-((r-r0-w)/w)^2)^2;')
dimer_force.addGlobalParameter('h', h) # barrier height
dimer_force.addGlobalParameter('r0', r0) # compact state separation
dimer_force.addGlobalParameter('w', w) # second minimum is at r0 + 2*w
dimer_force.addBond(0, 1, [])
system.addForce(dimer_force)
# Create initial coordinates using random positions.
coordinates = units.Quantity(numpy.zeros([2,3], numpy.float64), units.nanometer)
# Reposition dimer particles at compact minimum.
coordinates[0,:] *= 0.0
coordinates[1,:] *= 0.0
coordinates[1,0] = r0
# Return system and coordinates.
return [system, coordinates]
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment