External force applied but not effective

[Ross Gunderson]

Hey everyone,

Sorry to keep posting, but it seems when I respond the forum seems to die. Does anyone have an idea as to why an applied force on a set of atoms would be ineffective to changing the velocity? In my simulation, I am applying a constant force to all water in X, which I can see from the output force text file, but this does not impact the velocity in either a) the output velocity text file (from a velocity reporter), or b) when taking the derivative of the positions.

I've tried different ensembles, different integrators, but the same problem persists. I will attach to files, one showing the difference in applied force (X vs Y direction) and the other showing the velocities in X and Y direction.

Code: Select all

Here is the added force: 
shearForce = CustomExternalForce('fc_x*x')
        shearForce.addPerParticleParameter('fc_x')
        for residue in psf.residue_list:
            if residue.resname == 'TIP3':
                for atom in residue.atoms:
                    shearForce.addParticle(atom.idx, [inputs.fc_x])
        system.addForce(shearForce)

Then I am also applying a positional restraint on the center water atoms in a 5 angstrom width. Their indices are read from a text file. The reason we are doing this is because we DID notice with the exact same shearforce code, a simulation with a DOPC bilayer did show solvent moving with an increased velocity when applied.

Code: Select all

holdForce = CustomExternalForce('fc_divide*(px^2+py^2+pz^2); \
                                    px=min(dx, boxlx-dx); \
                                    py=min(dy, boxly-dy); \
                                    pz=min(dz, boxlz-dz); \
                                    dx=abs(x-x0); \
                                    dy=abs(y-y0); \
                                    dz=abs(z-z0);')
        holdForce.addGlobalParameter('boxlx', boxlx)
        holdForce.addGlobalParameter('boxly', boxly)
        holdForce.addGlobalParameter('boxlz', boxlz)
        holdForce.addPerParticleParameter('fc_divide')
        holdForce.addPerParticleParameter('x0')
        holdForce.addPerParticleParameter('y0')
        holdForce.addPerParticleParameter('z0')
        # now we check every atom and if the atom is in a certain z range, we add it.
        for line in open('restraints/water_pos.txt', 'r'):
            segments = line.strip().split()
            atom1 = int(segments[0])
            xpos  = crd.positions[atom1].value_in_unit(angstroms)[0]/10
            ypos  = crd.positions[atom1].value_in_unit(angstroms)[1]/10
            zpos  = crd.positions[atom1].value_in_unit(angstroms)[2]/10
            holdForce.addParticle(atom1, [100, xpos, ypos, zpos])
        system.addForce(holdForce)

Any thoughts would be greatly appreciated!

[Peter Eastman]

if residue.resname == 'TIP3':

I think you meant that to be 'HOH', not 'TIP3'. Topology objects always use standard PDB names for atoms and residues. That means water is called HOH.

For debugging problems like this, it's good to test each part of your code to make sure it's really doing what you intend. For example, add the line

Code: Select all

print(shearForce.getNumParticles())

to verify that you're actually applying it to the number of particles you intended.

[Ross Gunderson]

Well, I am sure from the graph, the force is applied. When running in equilibrium the two forces are centered on zero.

[Peter Eastman]

Can you post a complete working example? It's very hard to tell exactly what you're doing. For example, what does height have to do with it? Your energy doesn't depend on y or z. It's just a constant force applied along the x axis, independent of position. What exactly are you doing, what do you expect to happen, and what is happening instead?

[Ross Gunderson]

Sure! Basically we are trying to apply a solvent flow by applying a constant force to every atom in the solvent. We expect the profile to be similar to a Couette flow model, where the velocity of the flow varies with the height above the stationary plane [which explains why I plot them vs z position]. When we were working with a bilayer (our plane, but not stationary), we applied a force to all of the water molecules using the "shearForce" code. I then calculated the velocities with the derivative of the positions from the DCD file. (Attached are two plots of the bilayer before and after flow applied.) This gave us results we were happy with, but we wanted to take a step back to configure how much the flow rate depends on the applied force. Force this we wanted to use a waterbox. With a positional restraint we held the middle of the box still, forming a stationary plane, and applied the constant force to the solvent. The results were shown in the figures before. Even when visualizing with VMD, we do not notice solvent movement. (Although, perhaps the result is more turbulent than we anticipate.)

[Peter Eastman]

Can you provide a complete working example? Without that, I can't tell exactly what you're doing.

[Ross Gunderson]

Sure. Below is the file I use to run the simulation. Let me know if theres a detail you think is missing.

Code: Select all

"""
Generated by CHARMM-GUI (http://www.charmm-gui.org)

openmm_run.py

This program is OpenMM running scripts written in python.

Correspondance: jul316@lehigh.edu or wonpil@lehigh.edu
Last update: March 29, 2017
"""

from __future__ import print_function
import argparse
import sys
import os

from omm_readinputs import *
from omm_readparams import *
from omm_vfswitch import *
from omm_barostat import *
from omm_restraints import *
from omm_rewrap import *
from omm_reporters import *

from simtk.unit import *
from simtk.openmm import *
from simtk.openmm.app import *

parser = argparse.ArgumentParser()
parser.add_argument('-i', dest='inpfile', help='Input parameter file', required=True)
parser.add_argument('-p', dest='psffile', help='Input CHARMM PSF file', required=True)
parser.add_argument('-c', dest='crdfile', help='Input CHARMM CRD file', required=True)
parser.add_argument('-t', dest='toppar', help='Input CHARMM-GUI toppar stream file', required=True)
parser.add_argument('-b', dest='sysinfo', help='Input CHARMM-GUI sysinfo stream file (optional)', default=None)
parser.add_argument('-icrst', metavar='RSTFILE', dest='icrst', help='Input CHARMM RST file (optional)', default=None)
parser.add_argument('-irst', metavar='RSTFILE', dest='irst', help='Input restart file (optional)', default=None)
parser.add_argument('-ichk', metavar='CHKFILE', dest='ichk', help='Input checkpoint file (optional)', default=None)
parser.add_argument('-opdb', metavar='PDBFILE', dest='opdb', help='Output PDB file (optional)', default=None)
parser.add_argument('-orst', metavar='RSTFILE', dest='orst', help='Output restart file (optional)', default=None)
parser.add_argument('-ochk', metavar='CHKFILE', dest='ochk', help='Output checkpoint file (optional)', default=None)
parser.add_argument('-odcd', metavar='DCDFILE', dest='odcd', help='Output trajectory file (optional)', default=None)
parser.add_argument('-rewrap', dest='rewrap', help='Re-wrap the coordinates in a molecular basis (optional)', action='store_true', default=False)
args = parser.parse_args()

# Load parameters
print("Loading parameters")
inputs = read_inputs(args.inpfile)
params = read_params(args.toppar)
psf = read_psf(args.psffile)
crd = read_crd(args.crdfile)
if args.sysinfo:
    psf = read_box(psf, args.sysinfo)
else:
    psf = gen_box(psf, crd)
vel_file = args.odcd[:-3]+"txt"
# Build system
system = psf.createSystem(params, nonbondedMethod=inputs.coulomb,
                          nonbondedCutoff=inputs.r_off*nanometers,
                          constraints=inputs.cons,
                          ewaldErrorTolerance=inputs.ewald_Tol)

if inputs.vdw == 'Force-switch': system = vfswitch(system, psf, inputs)
if inputs.pcouple == 'yes':      system = barostat(system, inputs)
if inputs.rest == 'yes':         system = restraints(system, psf, crd, inputs)
integrator = LangevinIntegrator(inputs.temp*kelvin, inputs.fric_coeff/picosecond, inputs.dt*picoseconds)

exit

# Set platform
platform = Platform.getPlatformByName('CUDA')
prop = dict(CudaPrecision='single', DisablePmeStream='true')

# Build simulation context
simulation = Simulation(psf.topology, system, integrator, platform, prop)
simulation.context.setPositions(crd.positions)
if args.icrst:
    charmm_rst = read_charmm_rst(args.icrst)
    simulation.context.setPositions(charmm_rst.positions)
    simulation.context.setVelocities(charmm_rst.velocities)
    simulation.context.setPeriodicBoxVectors(charmm_rst.box[0], charmm_rst.box[1], charmm_rst.box[2])
if args.irst:
    with open(args.irst, 'r') as f:
        simulation.context.setState(XmlSerializer.deserialize(f.read()))
if args.ichk:
    with open(args.ichk, 'rb') as f:
        simulation.context.loadCheckpoint(f.read())

# Re-wrap
if args.rewrap:
    simulation = rewrap(simulation)

# Calculate initial system energy
print("\nInitial system energy")
print(simulation.context.getState(getEnergy=True).getPotentialEnergy())

# Energy minimization
if inputs.mini_nstep > 0:
    print("\nEnergy minimization: %s steps" % inputs.mini_nstep)
    simulation.minimizeEnergy(tolerance=inputs.mini_Tol*kilojoule/mole, maxIterations=inputs.mini_nstep)
    print(simulation.context.getState(getEnergy=True).getPotentialEnergy())

# Generate initial velocities
if inputs.gen_vel == 'yes':
    print("\nGenerate initial velocities")
    if inputs.gen_seed:
        simulation.context.setVelocitiesToTemperature(inputs.gen_temp, inputs.gen_seed)
    else:
        simulation.context.setVelocitiesToTemperature(inputs.gen_temp)
# set up velocity reporter
velocity_dis = simulation.context.getState( getVelocities=True )
velocity_reporter = VelReporter(vel_file, inputs.nstdcd)

# Production
if inputs.nstep > 0:
    print("\nMD run: %s steps" % inputs.nstep)
    if inputs.nstdcd > 0:
        if not args.odcd: args.odcd = 'output.dcd'
        simulation.reporters.append(DCDReporter(args.odcd, inputs.nstdcd))
        simulation.reporters.append(velocity_reporter)
    simulation.reporters.append(
        StateDataReporter(sys.stdout, inputs.nstout, step=True, time=True, potentialEnergy=True, temperature=True, progress=True,
                          remainingTime=True, speed=True, totalSteps=inputs.nstep, separator='\t')
    )
    simulation.step(inputs.nstep)

# Write restart file
if not (args.orst or args.ochk): args.orst = 'output.rst'
if args.orst:
    state = simulation.context.getState( getPositions=True )
    with open(args.orst, 'w') as f:
        f.write(XmlSerializer.serialize(state))
if args.ochk:
    with open(args.ochk, 'wb') as f:
        f.write(simulation.context.createCheckpoint())
if args.opdb:
    crd = simulation.context.getState(getPositions=True).getPositions()
    PDBFile.writeFile(psf.topology, crd, open(args.opdb, 'w'))

Here is my input file:

Code: Select all

nstep       = 500000                            # Number of steps to run
dt          = 0.002                             # Time-step (ps)

nstout      = 1000                              # Writing output frequency (steps)
nstdcd      = 10000                             # Writing coordinates trajectory frequency (steps)

coulomb     = PME                               # Electrostatic cut-off method
ewaldTol    = 0.0005                            # Ewald error tolerance
vdw         = Force-switch                      # vdW cut-off method
r_on        = 1.0                               # Switch-on distance (nm)
r_off       = 1.2                               # Switch-off distance (nm)

temp        = 310.15                            # Temperature (K)
fric_coeff  = 1                                 # Friction coefficient for Langevin dynamics

pcouple     = yes                               # Turn on/off pressure coupling
p_ref       = 1.0                               # Pressure (Pref or Pxx, Pyy, Pzz; bar)
p_type      = membrane                          # MonteCarloBarostat type
p_XYMode    = XYIsotropic                       # For MonteCarloMembraneBarostat
p_ZMode     = ZFree                             # For MonteCarloMembraneBarostat
p_tens      = 0.0                               # Sulface tension for MonteCarloMembraneBarostat (dyne/cm)
p_freq      = 100                               # Pressure coupling frequency (steps)

cons        = HBonds                            # Constraints mehtod

rest        = yes                               # Turn on/off restraints

fc_mpos	    = 0
fc_x        = 1
divide      = 1

[Ross Gunderson]

I also added the print statement to verify the custom forces are applied, and the numbers are correct.

[Peter Eastman]

I'm confused. There's nothing in this script to add a CustomExternalForce. What does it have to do with the problem you were asking about?

Once again, I need a complete working example. If I don't know what your code is doing, I can't tell you why it doesn't produce the result you expect.

[Ross Gunderson]

Okay. My apologies. The custom forces are loaded from an additional file along with some of the other settings. I will include them below. Almost all of these were obtained through CHARMM-GUI.

inputs are given in a separate file:

Code: Select all

nstep       = 500000                            # Number of steps to run
dt          = 0.002                             # Time-step (ps)
nstout      = 1000                              # Writing output frequency (steps)
nstdcd      = 10000                             # Writing coordinates trajectory frequency (steps)
coulomb     = PME                               # Electrostatic cut-off method
ewaldTol    = 0.0005                            # Ewald error tolerance
vdw         = Force-switch                      # vdW cut-off method
r_on        = 1.0                               # Switch-on distance (nm)
r_off       = 1.2                               # Switch-off distance (nm)
temp        = 310.15                            # Temperature (K)
fric_coeff  = 1                                 # Friction coefficient for Langevin dynamics
pcouple     = yes                               # Turn on/off pressure coupling
p_ref       = 1.0                               # Pressure (Pref or Pxx, Pyy, Pzz; bar)
p_type      = membrane                          # MonteCarloBarostat type
p_XYMode    = XYIsotropic                       # For MonteCarloMembraneBarostat
p_ZMode     = ZFree                             # For MonteCarloMembraneBarostat
p_tens      = 0.0                               # Sulface tension for MonteCarloMembraneBarostat (dyne/cm)
p_freq      = 100                               # Pressure coupling frequency (steps)
cons        = HBonds                            # Constraints mehtod
rest        = yes                               # Turn on/off restraints
fc_mpos	    = 0
fc_x        = 1
divide      = 1

Then read by another file:
omm_readinputs

Code: Select all

from simtk.unit import *
from simtk.openmm import *
from simtk.openmm.app import *

class _OpenMMReadInputs():
    def __init__(self):
        self.mini_nstep       = 0                                         # Number of steps for minimization
        self.mini_Tol         = 1.0                                       # Minimization energy tolerance
        self.gen_vel          = 'no'                                      # Generate initial velocities
        self.gen_temp         = 300.0                                     # Temperature for generating initial velocities (K)
        self.gen_seed         = None                                      # Seed for generating initial velocities
        self.nstep            = 0                                         # Number of steps to run
        self.dt               = 0.002                                     # Time-step (ps)
        self.nstout           = 100                                       # Writing output frequency (steps)
        self.nstdcd           = 0                                         # Wrtiing coordinates trajectory frequency (steps)
        self.coulomb          = PME                                       # Electrostatic cut-off method
        self.ewald_Tol        = 0.0005                                    # Ewald error tolerance
        self.vdw              = 'Force-switch'                            # vdW cut-off method
        self.r_on             = 1.0                                       # Switch-on distance (nm)
        self.r_off            = 1.2                                       # Switch-off distance (nm)
        self.temp             = 300.0                                     # Temperature (K)
        self.fric_coeff       = 1                                         # Friction coefficient for Langevin dynamics
        self.pcouple          = 'no'                                      # Turn on/off pressure coupling
        self.p_ref            = 1.0                                       # Pressure (Pref or Pxx, Pyy, Pzz; bar)
        self.p_type           = 'membrane'                                # MonteCarloBarotat type
        self.p_scale          = True, True, True                          # For MonteCarloAnisotropicBarostat
        self.p_XYMode         = MonteCarloMembraneBarostat.XYIsotropic    # For MonteCarloMembraneBarostat
        self.p_ZMode          = MonteCarloMembraneBarostat.ZFree          # For MonteCarloMembraneBarostat
        self.p_tens           = 0.0                                       # Sulface tension for MonteCarloMembraneBarostat (dyne/cm)
        self.p_freq           = 15                                        # Pressure coupling frequency (steps)
        self.cons             = HBonds                                    # Constraints method
        self.rest             = 'no'                                      # Turn on/off restraints
        self.fc_x             = 0.0                                       # shear force in x direction (kJ/mol/nm^2)
        self.fc_bb            = 0.0                                       # Positional restraint force constant for protein backbone (kJ/mol/nm^2)
        self.fc_sc            = 0.0                                       # Positional restraint force constant for protein side-chain (kJ/mol/nm^2)
        self.fc_mpos          = 0.0                                       # Positional restraint force constant for micelle lipids (kJ/mol/nm^2)
        self.fc_lpos          = 0.0                                       # Positional restraint force constant for lipids (kJ/mol/nm^2)
        self.fc_ldih          = 0.0                                       # Dihedral restraint force constant for lipids (kJ/mol/rad^2)
        self.fc_cdih          = 0.0                                       # Dihedral restraint force constant for carbohydrates (kJ/mol/rad^2)
        self.divide           = 0.0                                       # restraint force for the middle of the box.

    def read(self, inputFile):
        for line in open(inputFile, 'r'):
            if line.find('#') >= 0: line = line.split('#')[0]
            line = line.strip()
            if len(line) > 0:
                segments = line.split('=')
                input_param = segments[0].upper().strip()
                try:    input_value = segments[1].strip()
                except: input_value = None
                if input_value:
                    if input_param == 'MINI_NSTEP':                     self.mini_nstep       = int(input_value)
                    if input_param == 'MINI_TOL':                       self.mini_Tol         = float(input_value)
                    if input_param == 'GEN_VEL':
                        if input_value.upper() == 'YES':                self.gen_vel          = 'yes'
                        if input_value.upper() == 'NO':                 self.gen_vel          = 'no'
                    if input_param == 'GEN_TEMP':                       self.gen_temp         = float(input_value)
                    if input_param == 'GEN_SEED':                       self.gen_seed         = int(input_value)
                    if input_param == 'NSTEP':                          self.nstep            = int(input_value)
                    if input_param == 'DT':                             self.dt               = float(input_value)
                    if input_param == 'NSTOUT':                         self.nstout           = int(input_value)
                    if input_param == 'NSTDCD':                         self.nstdcd           = int(input_value)
                    if input_param == 'COULOMB':
                        if input_value.upper() == 'NOCUTOFF':           self.coulomb          = NoCutoff
                        if input_value.upper() == 'CUTOFFNONPERIODIC':  self.coulomb          = CutoffNonPeriodic
                        if input_value.upper() == 'CUTOFFPERIODIC':     self.coulomb          = CutoffPeriodic
                        if input_value.upper() == 'EWALD':              self.coulomb          = Ewald
                        if input_value.upper() == 'PME':                self.coulomb          = PME
                    if input_param == 'EWALD_TOL':                      self.ewald_Tol        = float(input_value)
                    if input_param == 'VDW':
                        if input_value.upper() == 'FORCE-SWITCH':       self.vdw              = 'Force-switch'
                    if input_param == 'R_ON':                           self.r_on             = float(input_value)
                    if input_param == 'R_OFF':                          self.r_off            = float(input_value)
                    if input_param == 'TEMP':                           self.temp             = float(input_value)
                    if input_param == 'FRIC_COEFF':                     self.fric_coeff       = float(input_value)
                    if input_param == 'PCOUPLE':
                        if input_value.upper() == 'YES':                self.pcouple          = 'yes'
                        if input_value.upper() == 'NO':                 self.pcouple          = 'no'
                    if input_param == 'P_REF':
                        if input_value.find(',') < 0:
                            self.p_ref = float(input_value)
                        else:
                            Pxx = float(input_value.split(',')[0])
                            Pyy = float(input_value.split(',')[1])
                            Pzz = float(input_value.split(',')[2])
                            self.p_ref = Pxx, Pyy, Pzz
                    if input_param == 'P_TYPE':
                        if input_value.upper() == 'ISOTROPIC':          self.p_type           = 'isotropic'
                        if input_value.upper() == 'MEMBRANE':           self.p_type           = 'membrane'
                        if input_value.upper() == 'ANISOTROPIC':        self.p_type           = 'anisotropic'
                    if input_param == 'P_SCALE':
                        scaleX = True
                        scaleY = True
                        scaleZ = True
                        if input_value.upper().find('X') < 0: scaleX = False
                        if input_value.upper().find('Y') < 0: scaleY = False
                        if input_value.upper().find('Z') < 0: scaleZ = False
                        self.p_scale = scaleX, scaleY, scaleZ
                    if input_param == 'P_XYMODE':
                        if input_value.upper() == 'XYISOTROPIC':        self.p_XYMode         = MonteCarloMembraneBarostat.XYIsotropic
                        if input_value.upper() == 'XYANISOTROPIC':      self.p_XYMode         = MonteCarloMembraneBarostat.XYAnisotropic
                    if input_param == 'P_ZMODE':
                        if input_value.upper() == 'ZFREE':              self.p_ZMode          = MonteCarloMembraneBarostat.ZFree
                        if input_value.upper() == 'ZFIXED':             self.p_ZMode          = MonteCarloMembraneBarostat.ZFixed
                        if input_value.upper() == 'CONSTANTVOLUME':     self.p_ZMode          = MonteCarloMembraneBarostat.ConstantVolume
                    if input_param == 'P_TENS':                         self.p_tens           = float(input_value)
                    if input_param == 'P_FREQ':                         self.p_freq           = int(input_value)
                    if input_param == 'CONS':
                        if input_value.upper() == 'NONE':               self.cons             = None
                        if input_value.upper() == 'HBONDS':             self.cons             = HBonds
                        if input_value.upper() == 'ALLBONDS':           self.cons             = AllBonds
                        if input_value.upper() == 'HANGLES':            self.cons             = HAngles
                    if input_param == 'REST':
                        if input_value.upper() == 'YES':                self.rest             = 'yes'
                        if input_value.upper() == 'NO':                 self.rest             = 'no'
                    if input_param == 'FC_X':                           self.fc_x             = float(input_value)
                    if input_param == 'FC_BB':                          self.fc_bb            = float(input_value)
                    if input_param == 'FC_BB':                          self.fc_bb            = float(input_value)
                    if input_param == 'FC_SC':                          self.fc_sc            = float(input_value)
                    if input_param == 'FC_MPOS':                        self.fc_mpos          = float(input_value)
                    if input_param == 'FC_LPOS':                        self.fc_lpos          = float(input_value)
                    if input_param == 'FC_LDIH':                        self.fc_ldih          = float(input_value)
                    if input_param == 'FC_CDIH':                        self.fc_cdih          = float(input_value)
                    if input_param == 'DIVIDE':                         self.divide           = float(input_value)
        return self
def read_inputs(inputFile):
    return _OpenMMReadInputs().read(inputFile)

We read the parameters of the system:
omm_readparams

Code: Select all

import os
from simtk.unit import *
from simtk.openmm import *
from simtk.openmm.app import *

def read_psf(filename):
    psf = CharmmPsfFile(filename)
    return psf

def read_crd(filename):
    crd = CharmmCrdFile(filename)
    return crd

def read_charmm_rst(filename):
    charmm_rst = CharmmRstFile(filename)
    for i, line in enumerate(charmm_rst.header):
        line = line.strip()
        words = line.split()
        if len(line) != 0:
            if words[0] == 'CRYSTAL' or words[0] == '!CRYSTAL':
                line1 = charmm_rst.header[i+1]
                line2 = charmm_rst.header[i+2]
                boxlx = Vec3(float(line1.split()[0].replace("D", "E")), 0.0, 0.0)
                boxly = Vec3(0.0, float(line1.split()[2].replace("D", "E")), 0.0)
                boxlz = Vec3(0.0, 0.0, float(line2.split()[2].replace("D", "E")))
                box = (boxlx, boxly, boxlz)
                break

    positions = charmm_rst.positionsold
    new_positions = []
    for position in positions:
        oldx = position[0]/angstrom
        oldy = position[1]/angstrom
        oldz = position[2]/angstrom
        newx = oldx + boxlx[0]/2.0
        newy = oldy + boxly[1]/2.0
        newz = oldz + boxlz[2]/2.0
        new_positions.append(Vec3(newx, newy, newz))
    charmm_rst.box = Quantity(box, angstroms)
    charmm_rst.positions = Quantity(new_positions, angstroms)
    return charmm_rst

def read_params(filename):
    extlist = ['rtf', 'prm', 'str']
    parFiles = ()
    for line in open(filename, 'r'):
        if '!' in line: line = line.split('!')[0]
        parfile = line.strip()
        if len(parfile) != 0:
            ext = parfile.lower().split('.')[-1]
            if not ext in extlist: continue
            parFiles += ( parfile, )
    params = CharmmParameterSet( *parFiles )
    return params

def read_box(psf, filename):
    for line in open(filename, 'r'):
        segments = line.split('=')
        if segments[0].strip() == "BOXLX": boxlx = float(segments[1])
        if segments[0].strip() == "BOXLY": boxly = float(segments[1])
        if segments[0].strip() == "BOXLZ": boxlz = float(segments[1])
    psf.setBox(boxlx*angstroms, boxly*angstroms, boxlz*angstroms)
    return psf

def gen_box(psf, crd):
    coords = crd.positions
    min_crds = [coords[0][0], coords[0][1], coords[0][2]]
    max_crds = [coords[0][0], coords[0][1], coords[0][2]]
    for coord in coords:
        min_crds[0] = min(min_crds[0], coord[0])
        min_crds[1] = min(min_crds[1], coord[1])
        min_crds[2] = min(min_crds[2], coord[2])
        max_crds[0] = max(max_crds[0], coord[0])
        max_crds[1] = max(max_crds[1], coord[1])
        max_crds[2] = max(max_crds[2], coord[2])
    boxlx = max_crds[0]-min_crds[0]
    boxly = max_crds[1]-min_crds[1]
    boxlz = max_crds[2]-min_crds[2]
    psf.setBox(boxlx, boxly, boxlz)
    return psf

The file that handles PBC:
omm_rewrap

Code: Select all

from __future__ import print_function
import os
from math import *
from simtk.unit import *
from simtk.openmm import *
from simtk.openmm.app import *

def rewrap(simulation):
    bonds = simulation.topology.bonds()
    positions = simulation.context.getState(getPositions=True).getPositions()
    box = simulation.context.getState().getPeriodicBoxVectors()
    boxlx = box[0][0]/angstrom
    boxly = box[1][1]/angstrom
    boxlz = box[2][2]/angstrom
    min_crds = [positions[0][0]/angstrom, positions[0][1]/angstrom, positions[0][2]/angstrom]
    max_crds = [positions[0][0]/angstrom, positions[0][1]/angstrom, positions[0][2]/angstrom]
    for position in positions:
        min_crds[0] = min(min_crds[0], position[0]/angstrom)
        min_crds[1] = min(min_crds[1], position[1]/angstrom)
        min_crds[2] = min(min_crds[2], position[2]/angstrom)
        max_crds[0] = max(max_crds[0], position[0]/angstrom)
        max_crds[1] = max(max_crds[1], position[1]/angstrom)
        max_crds[2] = max(max_crds[2], position[2]/angstrom)
    xcen = (max_crds[0] + min_crds[0]) / 2.0
    ycen = (max_crds[1] + min_crds[1]) / 2.0
    zcen = (max_crds[2] + min_crds[2]) / 2.0
    for bond in bonds:
        atom1 = bond[0]
        atom2 = bond[1]
        atom1id = atom1.index
        atom2id = atom2.index
        res1 = atom1.residue
        res2 = atom2.residue
        x1, y1, z1 = positions[atom1id]
        x2, y2, z2 = positions[atom2id]
        dx = fabs(x1/angstrom - x2/angstrom)
        dy = fabs(y1/angstrom - y2/angstrom)
        dz = fabs(z1/angstrom - z2/angstrom)
        if dx > boxlx/2 or dy > boxly/2 or dz > boxlz/2:
            for atom in res2.atoms():
                oldx = positions[atom.index][0]/angstrom
                oldy = positions[atom.index][1]/angstrom
                oldz = positions[atom.index][2]/angstrom
                if dx > boxlx/2.0:
                    if oldx < xcen: newx = oldx + boxlx
                    else: newx = oldx - boxlx
                else:
                    newx = oldx
                if dy > boxly/2.0:
                    if oldy < ycen: newy = oldy + boxly
                    else: newy = oldy - boxly
                else:
                    newy = oldy
                if dz > boxlz/2.0:
                    if oldz < zcen: newz = oldz + boxlz
                    else: newz = oldz - boxlz
                else:
                    newz = oldz
                new_position = Vec3(newx, newy, newz)
                positions[atom.index] = Quantity(new_position, angstroms)
    simulation.context.setPositions(positions)
    return simulation

omm_restraints:

Code: Select all

from simtk.unit import *
from simtk.openmm import *
from simtk.openmm.app import *

def restraints(system, psf, crd, inputs):

    boxlx = system.getDefaultPeriodicBoxVectors()[0][0].value_in_unit(nanometers)
    boxly = system.getDefaultPeriodicBoxVectors()[1][1].value_in_unit(nanometers)
    boxlz = system.getDefaultPeriodicBoxVectors()[2][2].value_in_unit(nanometers)

    if inputs.fc_x > 0:
        # position restraints for protein
        shearForce = CustomExternalForce('fc_x*x')
        shearForce.addPerParticleParameter('fc_x')
        for residue in psf.residue_list:
            if residue.resname == 'TIP3':
                for atom in residue.atoms:
                    shearForce.addParticle(atom.idx, [inputs.fc_x])
        system.addForce(shearForce)
        
    if inputs.divide > 0:
    # we divide the box in two by applying a strong force on the atoms in the middle of the box.
        holdForce = CustomExternalForce('fc_divide*(px^2+py^2+pz^2); \
                                    px=min(dx, boxlx-dx); \
                                    py=min(dy, boxly-dy); \
                                    pz=min(dz, boxlz-dz); \
                                    dx=abs(x-x0); \
                                    dy=abs(y-y0); \
                                    dz=abs(z-z0);')
        holdForce.addGlobalParameter('boxlx', boxlx)
        holdForce.addGlobalParameter('boxly', boxly)
        holdForce.addGlobalParameter('boxlz', boxlz)
        holdForce.addPerParticleParameter('fc_divide')
        holdForce.addPerParticleParameter('x0')
        holdForce.addPerParticleParameter('y0')
        holdForce.addPerParticleParameter('z0')
        # now we check every atom and if the atom is in a certain z range, we add it.
        for line in open('restraints/water_pos.txt', 'r'):
            segments = line.strip().split()
            atom1 = int(segments[0])
            xpos  = crd.positions[atom1].value_in_unit(angstroms)[0]/10
            ypos  = crd.positions[atom1].value_in_unit(angstroms)[1]/10
            zpos  = crd.positions[atom1].value_in_unit(angstroms)[2]/10
            holdForce.addParticle(atom1, [100, xpos, ypos, zpos])
        system.addForce(holdForce)
        
    return system

omm_barostat:

Code: Select all

from simtk.unit import *
from simtk.openmm import *
from simtk.openmm.app import *

def barostat(system, inputs):
    if inputs.p_type == 'isotropic':
        barostat = MonteCarloBarostat( inputs.p_ref*bar, inputs.temp*kelvin )
    if inputs.p_type == 'membrane':
        inputs.p_tens = inputs.p_tens*10.0
        barostat = MonteCarloMembraneBarostat( inputs.p_ref*bar, inputs.p_tens*bar*nanometers, inputs.temp*kelvin, inputs.p_XYMode, inputs.p_ZMode, inputs.p_freq )
    if inputs.p_type == 'anisotropic':
        barostat = MonteCarloAnisotropicBarostat( inputs.p_ref*bar, inputs.temp*kelvin, inputs.p_scale[0], inputs.p_scale[1], inputs.p_scale[2], inputs.p_freq )
    system.addForce(barostat)
    return system

The file to write force and velocity:
omm_reporters

Code: Select all

from simtk.unit import *
from simtk.openmm import *
from simtk.openmm.app import *

class ForceReporter(object):
    def __init__(self, file, reportInterval):
        self._out = open(file, 'w')
        self._reportInterval = reportInterval
    def __del__(self):
        self._out.close()
    def describeNextReport(self, simulation):
        steps = self._reportInterval - simulation.currentStep%self._reportInterval
        return (steps, False, False, True, False, None)
    def report(self, simulation, state):
        forces = state.getForces().value_in_unit(kilojoules/mole/nanometer)
        for f in forces:
            self._out.write('%g %g %g\n' % (f[0], f[1], f[2]))

class VelReporter(object):
    def __init__(self, file, reportInterval):
        self._out = open(file, 'w')
        self._reportInterval = reportInterval
    def __del__(self):
        self._out.close()
    def describeNextReport(self, simulation):
        steps = self._reportInterval - simulation.currentStep%self._reportInterval
        return (steps, False, True, False, False, None)
    def report(self, simulation, state):
        velocities = state.getVelocities().value_in_unit(nanometer/picosecond)
        for v in velocities:
            self._out.write('%g %g %g\n' % (v[0], v[1], v[2]))

vfswitching file:
omm_vfswitch

Code: Select all

from simtk.unit import *
from simtk.openmm import *
from simtk.openmm.app import *

def vfswitch(system, psf, inputs):
    r_on = inputs.r_on
    r_off = inputs.r_off

    # custom nonbonded force for force-switch
    chknbfix = False
    for force in system.getForces():
        if isinstance(force, NonbondedForce):
            nonbonded = force
        if isinstance(force, CustomNonbondedForce):
            nbfix     = force
            chknbfix  = True

    # vfswitch
    vfswitch = CustomNonbondedForce('step(Ron-r)*(ccnba*tr6*tr6-ccnbb*tr6+ccnbb*onoff3-ccnba*onoff6) \
                                     +step(r-Ron)*step(Roff-r)*(cr12*rjunk6 - cr6*rjunk3) \
                                     -step(r-Ron)*step(Ron-r)*(cr12*rjunk6 - cr6*rjunk3); \
                                     cr6  = ccnbb*ofdif3*rjunk3; \
                                     cr12 = ccnba*ofdif6*rjunk6; \
                                     rjunk3 = r3-recof3; \
                                     rjunk6 = tr6-recof6; \
                                     r3 = r1*tr2; \
                                     r1 = sqrt(tr2); \
                                     tr6 = tr2 * tr2 * tr2; \
                                     tr2 = 1.0/s2; \
                                     s2 = r*r; \
                                     ccnbb = 4.0*epsilon*sigma^6; \
                                     ccnba = 4.0*epsilon*sigma^12; \
                                     sigma = sigma1+sigma2; \
                                     epsilon = epsilon1*epsilon2; \
                                     onoff3 = recof3/on3; \
                                     onoff6 = recof6/on6; \
                                     ofdif3 = off3/(off3 - on3); \
                                     ofdif6 = off6/(off6 - on6); \
                                     recof3 = 1.0/off3; \
                                     on6 = on3*on3; \
                                     on3 = c2onnb*Ron; \
                                     recof6 = 1.0/off6; \
                                     off6 = off3*off3; \
                                     off3 = c2ofnb*Roff; \
                                     c2ofnb = Roff*Roff; \
                                     c2onnb = Ron*Ron; \
                                     Ron  = %f; \
                                     Roff = %f;' % (r_on, r_off) )
    vfswitch.addPerParticleParameter('sigma')
    vfswitch.addPerParticleParameter('epsilon')
    vfswitch.setNonbondedMethod(vfswitch.CutoffPeriodic)
    vfswitch.setCutoffDistance(nonbonded.getCutoffDistance())
    for i in range(nonbonded.getNumParticles()):
        chg, sig, eps = nonbonded.getParticleParameters(i)
        nonbonded.setParticleParameters(i, chg, 0.0, 0.0) # zero-out LJ
        sig = sig*0.5
        eps = eps**0.5
        vfswitch.addParticle([sig, eps])
    for i in range(nonbonded.getNumExceptions()):
        atom1, atom2 = nonbonded.getExceptionParameters(i)[:2]
        vfswitch.addExclusion(atom1, atom2)
    vfswitch.setForceGroup(psf.NONBONDED_FORCE_GROUP)
    system.addForce(vfswitch)

    # vfswitch14
    vfswitch14 = CustomBondForce('step(Ron-r)*(ccnba*tr6*tr6-ccnbb*tr6+ccnbb*onoff3-ccnba*onoff6) \
                                  +step(r-Ron)*step(Roff-r)*(cr12*rjunk6 - cr6*rjunk3) \
                                  -step(r-Ron)*step(Ron-r)*(cr12*rjunk6 - cr6*rjunk3); \
                                  cr6  = ccnbb*ofdif3*rjunk3; \
                                  cr12 = ccnba*ofdif6*rjunk6; \
                                  rjunk3 = r3-recof3; \
                                  rjunk6 = tr6-recof6; \
                                  r3 = r1*tr2; \
                                  r1 = sqrt(tr2); \
                                  tr6 = tr2 * tr2 * tr2; \
                                  tr2 = 1.0/s2; \
                                  s2 = r*r; \
                                  ccnbb = 4.0*epsilon*sigma^6; \
                                  ccnba = 4.0*epsilon*sigma^12; \
                                  onoff3 = recof3/on3; \
                                  onoff6 = recof6/on6; \
                                  ofdif3 = off3/(off3 - on3); \
                                  ofdif6 = off6/(off6 - on6); \
                                  recof3 = 1.0/off3; \
                                  on6 = on3*on3; \
                                  on3 = c2onnb*Ron; \
                                  recof6 = 1.0/off6; \
                                  off6 = off3*off3; \
                                  off3 = c2ofnb*Roff; \
                                  c2ofnb = Roff*Roff; \
                                  c2onnb = Ron*Ron; \
                                  Ron  = %f; \
                                  Roff = %f;' % (r_on, r_off) )
    vfswitch14.addPerBondParameter('sigma')
    vfswitch14.addPerBondParameter('epsilon')
    for i in range(nonbonded.getNumExceptions()):
        atom1, atom2, chg, sig, eps = nonbonded.getExceptionParameters(i)
        nonbonded.setExceptionParameters(i, atom1, atom2, chg, 0.0, 0.0) # zero-out LJ14
        vfswitch14.addBond(atom1, atom2, [sig, eps])
    system.addForce(vfswitch14)

    # vfswitch_NBFIX
    if chknbfix:
        nbfix.setEnergyFunction('step(Ron-r)*(ccnba*tr6*tr6-ccnbb*tr6+ccnbb*onoff3-ccnba*onoff6) \
                                 +step(r-Ron)*step(Roff-r)*(cr12*rjunk6 - cr6*rjunk3) \
                                 -step(r-Ron)*step(Ron-r)*(cr12*rjunk6 - cr6*rjunk3); \
                                 cr6  = ccnbb*ofdif3*rjunk3; \
                                 cr12 = ccnba*ofdif6*rjunk6; \
                                 rjunk3 = r3-recof3; \
                                 rjunk6 = tr6-recof6; \
                                 r3 = r1*tr2; \
                                 r1 = sqrt(tr2); \
                                 tr6 = tr2 * tr2 * tr2; \
                                 tr2 = 1.0/s2; \
                                 s2 = r*r; \
                                 ccnbb = bcoef(type1, type2); \
                                 ccnba = acoef(type1, type2)^2; \
                                 onoff3 = recof3/on3; \
                                 onoff6 = recof6/on6; \
                                 ofdif3 = off3/(off3 - on3); \
                                 ofdif6 = off6/(off6 - on6); \
                                 recof3 = 1.0/off3; \
                                 on6 = on3*on3; \
                                 on3 = c2onnb*Ron; \
                                 recof6 = 1.0/off6; \
                                 off6 = off3*off3; \
                                 off3 = c2ofnb*Roff; \
                                 c2ofnb = Roff*Roff; \
                                 c2onnb = Ron*Ron; \
                                 Ron  = %f; \
                                 Roff = %f;' % (r_on, r_off) )

    return system

Then, I run the simulation with this file:
openmm_run

Code: Select all

from __future__ import print_function
import argparse
import sys
import os

from omm_readinputs import *
from omm_readparams import *
from omm_vfswitch import *
from omm_barostat import *
from omm_restraints import *
from omm_rewrap import *
from omm_reporters import *

from simtk.unit import *
from simtk.openmm import *
from simtk.openmm.app import *

parser = argparse.ArgumentParser()
parser.add_argument('-i', dest='inpfile', help='Input parameter file', required=True)
parser.add_argument('-p', dest='psffile', help='Input CHARMM PSF file', required=True)
parser.add_argument('-c', dest='crdfile', help='Input CHARMM CRD file', required=True)
parser.add_argument('-t', dest='toppar', help='Input CHARMM-GUI toppar stream file', required=True)
parser.add_argument('-b', dest='sysinfo', help='Input CHARMM-GUI sysinfo stream file (optional)', default=None)
parser.add_argument('-icrst', metavar='RSTFILE', dest='icrst', help='Input CHARMM RST file (optional)', default=None)
parser.add_argument('-irst', metavar='RSTFILE', dest='irst', help='Input restart file (optional)', default=None)
parser.add_argument('-ichk', metavar='CHKFILE', dest='ichk', help='Input checkpoint file (optional)', default=None)
parser.add_argument('-opdb', metavar='PDBFILE', dest='opdb', help='Output PDB file (optional)', default=None)
parser.add_argument('-orst', metavar='RSTFILE', dest='orst', help='Output restart file (optional)', default=None)
parser.add_argument('-ochk', metavar='CHKFILE', dest='ochk', help='Output checkpoint file (optional)', default=None)
parser.add_argument('-odcd', metavar='DCDFILE', dest='odcd', help='Output trajectory file (optional)', default=None)
parser.add_argument('-rewrap', dest='rewrap', help='Re-wrap the coordinates in a molecular basis (optional)', action='store_true', default=False)
args = parser.parse_args()

# Load parameters
print("Loading parameters")
inputs = read_inputs(args.inpfile)
params = read_params(args.toppar)
psf = read_psf(args.psffile)
crd = read_crd(args.crdfile)
if args.sysinfo:
    psf = read_box(psf, args.sysinfo)
else:
    psf = gen_box(psf, crd)
vel_file = args.odcd[:-3]+"txt"
# Build system
system = psf.createSystem(params, nonbondedMethod=inputs.coulomb,
                          nonbondedCutoff=inputs.r_off*nanometers,
                          constraints=inputs.cons,
                          ewaldErrorTolerance=inputs.ewald_Tol)

if inputs.vdw == 'Force-switch': system = vfswitch(system, psf, inputs)
if inputs.pcouple == 'yes':      system = barostat(system, inputs)
if inputs.rest == 'yes':         system = restraints(system, psf, crd, inputs)
integrator = LangevinIntegrator(inputs.temp*kelvin, inputs.fric_coeff/picosecond, inputs.dt*picoseconds)

exit

# Set platform
platform = Platform.getPlatformByName('CUDA')
prop = dict(CudaPrecision='single', DisablePmeStream='true')

# Build simulation context
simulation = Simulation(psf.topology, system, integrator, platform, prop)
simulation.context.setPositions(crd.positions)
if args.icrst:
    charmm_rst = read_charmm_rst(args.icrst)
    simulation.context.setPositions(charmm_rst.positions)
    simulation.context.setVelocities(charmm_rst.velocities)
    simulation.context.setPeriodicBoxVectors(charmm_rst.box[0], charmm_rst.box[1], charmm_rst.box[2])
if args.irst:
    with open(args.irst, 'r') as f:
        simulation.context.setState(XmlSerializer.deserialize(f.read()))
if args.ichk:
    with open(args.ichk, 'rb') as f:
        simulation.context.loadCheckpoint(f.read())

# Re-wrap
if args.rewrap:
    simulation = rewrap(simulation)

# Calculate initial system energy
print("\nInitial system energy")
print(simulation.context.getState(getEnergy=True).getPotentialEnergy())

# Energy minimization
if inputs.mini_nstep > 0:
    print("\nEnergy minimization: %s steps" % inputs.mini_nstep)
    simulation.minimizeEnergy(tolerance=inputs.mini_Tol*kilojoule/mole, maxIterations=inputs.mini_nstep)
    print(simulation.context.getState(getEnergy=True).getPotentialEnergy())

# Generate initial velocities
if inputs.gen_vel == 'yes':
    print("\nGenerate initial velocities")
    if inputs.gen_seed:
        simulation.context.setVelocitiesToTemperature(inputs.gen_temp, inputs.gen_seed)
    else:
        simulation.context.setVelocitiesToTemperature(inputs.gen_temp)
# set up velocity reporter
velocity_dis = simulation.context.getState( getVelocities=True )
velocity_reporter = VelReporter(vel_file, inputs.nstdcd)

# Production
if inputs.nstep > 0:
    print("\nMD run: %s steps" % inputs.nstep)
    if inputs.nstdcd > 0:
        if not args.odcd: args.odcd = 'output.dcd'
        simulation.reporters.append(DCDReporter(args.odcd, inputs.nstdcd))
        simulation.reporters.append(velocity_reporter)
    simulation.reporters.append(
        StateDataReporter(sys.stdout, inputs.nstout, step=True, time=True, potentialEnergy=True, temperature=True, progress=True,
                          remainingTime=True, speed=True, totalSteps=inputs.nstep, separator='\t')
    )
    simulation.step(inputs.nstep)

# Write restart file
if not (args.orst or args.ochk): args.orst = 'output.rst'
if args.orst:
    state = simulation.context.getState( getPositions=True )
    with open(args.orst, 'w') as f:
        f.write(XmlSerializer.serialize(state))
if args.ochk:
    with open(args.ochk, 'wb') as f:
        f.write(simulation.context.createCheckpoint())
if args.opdb:
    crd = simulation.context.getState(getPositions=True).getPositions()
    PDBFile.writeFile(psf.topology, crd, open(args.opdb, 'w'))