Combining NonbondedForce and CustomNonbondedForce

[Joe Napoli]

This question relates to https://simtk.org/forums/viewtopic.php?f=161&t=5299.

I have tried to implement a CustomNonbondedForce in an xml file. I'd like the NonbondedForce of my system to handle just electrostatic interactions, so I've set epsilon = 0 for each particle in NonbondedForce in restrict NonbondedForce to computing electrostatic interactions.

I'd like interactions of the CustomNonbondedForce to be restricted to Oxygen atoms of my system, so I have specified it as:

Code: Select all

<CustomNonbondedForce energy="2*epsilon/(1-3/(gamma+3))*(sigma^6/(sigma^6+r^6))*(3/(gamma+3)*exp(gamma*(1-r/sigma))-1); sigma=0.5*(sigma1+sigma2); epsilon=sqr    t(epsilon1*epsilon2); gamma=gamma1+gamma2">
  <PerParticleParameter name="epsilon"/>
  <PerParticleParameter name="gamma"/>
  <PerParticleParameter name="sigma"/>
  <Atom class="OW" epsilon="1.0334" gamma="6.628" sigma="0.36174"/>
  <Atom class="HW" epsilon="0" gamma="1" sigma="1"/>
  <Atom class="MW" epsilon="0" gamma="1" sigma="1"/>
</CustomNonbondedForce>

Hence I've set epsilon = 0 for "HW" and "MW" atoms in order to eliminate interactions involving those atoms. My system is generated using:

Code: Select all

system = forcefield.createSystem(pdb.topology, nonbondedMethod=PME, nonbondedCutoff=0.85*nanometer, constraints=None, rigidWater=False)

My issue is that an exception is thrown that says that the energy is infinite, but it is unclear to me why this would be the case as the CustomNonbondedForce as written and with the above parameters is similar in form to the standard LJ potential. Therefore something must be wrong with my specification of the force, but I am at a loss as to what the issue is. Any input would be greatly appreciated.

Thank you,
Joe

[Peter Eastman]

Do you have files you could post that I could use to reproduce this? If so, I'll take a look at it. Aside from that, here are a few thoughts.

Are you certain the CustomNonbondedForce is the one producing the infinite energy?

If so, try removing pieces of the custom expression one at a time. Can you narrow down to what piece of the expression is causing the infinity?

Peter

[Joe Napoli]

Dear Peter,

Thanks for the input. My mistake - it's only after I set epsilon = 0 for Oxygen in the normal NonbondedForce that the energy becomes infinite. NonbondedForce with both electrostatic and LJ contributions works fine, as does my CustomNonbondedForce alone. Therefore, this suggests it is the electrostatic contribution that is diverging. I have attached the related xml.




Best,
Joe

[Joe Napoli]

Apologies, the website wouldn't allow me to attach my script and the structure pdb as well. Here is the python script. I can email the pdb file if you'd like to run it.

Thank you,
Joe

run.py

Code: Select all

from __future__ import print_function
from simtk.openmm import *
from simtk.openmm.app import *  
from simtk.unit import *
from sys import stdout
import numpy as np
import sys
import time

# Define the system temperature
sys_temp = 300.0

# Index of cuda device to use
device = sys.argv[1]
copies = int(sys.argv[2])
contraction = int(sys.argv[3])

print('------------------------------------------')
print(' Classical simulation of q-TIP4P/F water model ')
print('------------------------------------------')

# This file will be used to extract the topology. Bead cooridnates taken from equilibrated files (see below)
pdb = PDBFile('qtip4pf_216.pdb')

# Here we use the qTIP4P/F water model
forcefield = ForceField('qtip4pf_custom2.xml')

# Initialize simulation parameters 
system = forcefield.createSystem(pdb.topology, nonbondedMethod=PME, nonbondedCutoff=0.85*nanometer, constraints=None, rigidWater=False)

integrator = LangevinIntegrator(sys_temp*kelvin, 1.0/picosecond, 0.5*femtosecond)

# Set up the simulation, given the input and directions above
platform = Platform.getPlatformByName('CUDA')
properties = {'CudaPrecision': 'mixed', 'CudaDeviceIndex': device}
simulation = Simulation(pdb.topology, system, integrator, platform, properties)

# This line will initialize all copies to the same coordinates
simulation.context.setPositions(pdb.positions)

# Short equilibration (with reset velocities)
simulation.context.setVelocitiesToTemperature(sys_temp*kelvin)

print('# Running Production...')
# This script will print out a separate energy and position file for each bead

# Number of steps between printing
nprint = 2000
# Total number of print cycles
ncycles = 1000
# Total number of steps is nprint*ncycles

ereport = StateDataReporter('energies.dat', nprint, step=True, time=True, 
                              kineticEnergy=True,potentialEnergy=True,totalEnergy=True,
                              temperature=True,density=True,volume=False,speed=True,separator=' ')

for icycle in range(1,ncycles):
    simulation.step(nprint) 
    temp_state = simulation.context.getState(getPositions=True,getVelocities=False,getForces=True,getEnergy=True,getParameters=False,enforcePeriodicBox=True,groups=-1) 
    ereport.report(simulation,temp_state)

[Peter Eastman]

> I can email the pdb file if you'd like to run it.

That would be great. Thanks. It looks like there's also an XML file it needs.

Peter