Yeah I'd tested that out originally before I saw how you had done the calculation with asin, and I had the same issue.
I put it in again, for both the acos and asin methods, and tested them alone and together with the step function you had suggested, and the same problem occurs in all cases.
adding functions to CustomCompoundBondForce
-
Aron Broom
- Posts: 54
- Joined: Tue Mar 13, 2012 11:33 am
-
Peter Eastman
- Posts: 2588
- Joined: Thu Aug 09, 2007 1:25 pm
Re: adding functions to CustomCompoundBondForce
It sounds like that isn't the problem then.
To debug this, I suggest creating a minimal test case with just a few atoms and an energy function for which you know the correct result in advance. For example, arrange four atoms to form a particular angle, then set your energy expression equal to angle1 or angle2. Ask OpenMM to evaluate the energy and see whether the result is correct. Once you verify that each method of computing the angle is working correctly in its target domain, set the energy expression equal to "select". Is it correctly picking which of the two angles is more accurate? Once that is working correctly, form the combination of the two angles. Is that being computed correctly?
By doing this, you can verify your expression one piece at a time and track down where the problem is.
Peter
To debug this, I suggest creating a minimal test case with just a few atoms and an energy function for which you know the correct result in advance. For example, arrange four atoms to form a particular angle, then set your energy expression equal to angle1 or angle2. Ask OpenMM to evaluate the energy and see whether the result is correct. Once you verify that each method of computing the angle is working correctly in its target domain, set the energy expression equal to "select". Is it correctly picking which of the two angles is more accurate? Once that is working correctly, form the combination of the two angles. Is that being computed correctly?
By doing this, you can verify your expression one piece at a time and track down where the problem is.
Peter
-
Aron Broom
- Posts: 54
- Joined: Tue Mar 13, 2012 11:33 am
Re: adding functions to CustomCompoundBondForce
Sounds like a plan. I'll give it a shot and report back on the findings.
~Aron
~Aron
-
Aron Broom
- Posts: 54
- Joined: Tue Mar 13, 2012 11:33 am
Re: adding functions to CustomCompoundBondForce
I had a go at trying to isolate the problem, there are some odd things:
First, I tried my expression with just a 3 atom system with no bonds, and it works well under all conditions. In fact, just using the acos method even seems to be stable at 180 degress
Second, I moved to a 3 water system, using all 3 atoms of each as the geometric center for the restraint, and that is where problems come in. Using the simple acos() method, gives nan for the potential energy of the system at 180 degrees. Even if that is placed within a min(max( nest to try and constrain it, it similarly becomes nan, and the same is true for step it seems.
I noticed a really odd thing when I actually tried to run a simulation with an extremely small time step rather than just asking for the energy. The waters stayed in place, but rotated rapidly around their centers of geometry and the entire setup had a very high, but still defined potential energy.
I'm not entirely sure I grasp what is likely happening, and particularly why it works fine with single atoms but not a center.
Anyway, the end conclusion is that with the two methods of calculating the angle available, I can just choose the one that is appropriate. I don't think this is an issue with my expression at this point, so probably beyond the scope of what I want to accomplish.
Thanks for all the help and suggestions Peter.
First, I tried my expression with just a 3 atom system with no bonds, and it works well under all conditions. In fact, just using the acos method even seems to be stable at 180 degress
Second, I moved to a 3 water system, using all 3 atoms of each as the geometric center for the restraint, and that is where problems come in. Using the simple acos() method, gives nan for the potential energy of the system at 180 degrees. Even if that is placed within a min(max( nest to try and constrain it, it similarly becomes nan, and the same is true for step it seems.
I noticed a really odd thing when I actually tried to run a simulation with an extremely small time step rather than just asking for the energy. The waters stayed in place, but rotated rapidly around their centers of geometry and the entire setup had a very high, but still defined potential energy.
I'm not entirely sure I grasp what is likely happening, and particularly why it works fine with single atoms but not a center.
Anyway, the end conclusion is that with the two methods of calculating the angle available, I can just choose the one that is appropriate. I don't think this is an issue with my expression at this point, so probably beyond the scope of what I want to accomplish.
Thanks for all the help and suggestions Peter.
-
Peter Eastman
- Posts: 2588
- Joined: Thu Aug 09, 2007 1:25 pm
Re: adding functions to CustomCompoundBondForce
Very strange. Can you send me your test case? I'd like to try to figure out what's going on.
Peter
Peter
-
Aron Broom
- Posts: 54
- Joined: Tue Mar 13, 2012 11:33 am
Re: adding functions to CustomCompoundBondForce
Peter,
So I was re-writing the code to have something nice and clean to send you, when I realized a potential logic flaw in that the min/max nest I'd used for the angles looked like this:
angle1 = 'acos(min(max(' + cosine + ', -1), 1))'
and a similar thing for the asin portion. But, of course we don't want it to hit +/- 1. When I replaced that with:
angle1 = 'acos(min(max(' + cosine + ', -0.99), 0.99))'
and then use the step function you suggested it works. So I guess the trick was that I hadn't headed off the problem early enough. I'm going to continue to test it out, as I feel a little suspicious that it works, as I'm sure I tried a good many different permutations of things, but I'm hopeful.
If you want the code anyway let me know, and I'll post here again if it turns out to have not actually worked.
~Aron
So I was re-writing the code to have something nice and clean to send you, when I realized a potential logic flaw in that the min/max nest I'd used for the angles looked like this:
angle1 = 'acos(min(max(' + cosine + ', -1), 1))'
and a similar thing for the asin portion. But, of course we don't want it to hit +/- 1. When I replaced that with:
angle1 = 'acos(min(max(' + cosine + ', -0.99), 0.99))'
and then use the step function you suggested it works. So I guess the trick was that I hadn't headed off the problem early enough. I'm going to continue to test it out, as I feel a little suspicious that it works, as I'm sure I tried a good many different permutations of things, but I'm hopeful.
If you want the code anyway let me know, and I'll post here again if it turns out to have not actually worked.
~Aron
-
Peter Eastman
- Posts: 2588
- Joined: Thu Aug 09, 2007 1:25 pm
Re: adding functions to CustomCompoundBondForce
Go ahead and post it. acos of 1 or -1 should still be ok, so I'd like to take a look at it.
Peter
Peter
-
Aron Broom
- Posts: 54
- Joined: Tue Mar 13, 2012 11:33 am
Re: adding functions to CustomCompoundBondForce
Ok, I'm interested to know what is actually happening with the forces.
I've attached a zip of the relevant files, so you can just run AngleTestCase.py
In the restraint string generating function I've left it as -1 and 1, which should fail with NaN errors. I've put the -0.99 and 0.99 lines commented out below, which fix things.
Finally, not sure if this will help, but it's setup to do 0.5 fs timesteps at 100 kelvin, which will fail. If you change to timestep down to 0.001 fs, you can get the simulation to run, but the angle between the waters moves through a large range, while each water molecule itself rotates at an incredible speed around its own center. Weird stuff.
Also, for the record, I've written out below the code.
And here is a module for generating the restraint string and doing a quick angle calculation for some on-screen output. Currently it is expected to be saved as AngleRestraint.py
I've attached a zip of the relevant files, so you can just run AngleTestCase.py
In the restraint string generating function I've left it as -1 and 1, which should fail with NaN errors. I've put the -0.99 and 0.99 lines commented out below, which fix things.
Finally, not sure if this will help, but it's setup to do 0.5 fs timesteps at 100 kelvin, which will fail. If you change to timestep down to 0.001 fs, you can get the simulation to run, but the angle between the waters moves through a large range, while each water molecule itself rotates at an incredible speed around its own center. Weird stuff.
Also, for the record, I've written out below the code.
Code: Select all
# script to test the oddity surrounding COG or COM angular restraints using CustomCompoundBondForce in OpenMM
# import the OpenMM stuff
from simtk.openmm.app import *
from simtk.openmm import *
from simtk.unit import *
# import stdout so that we can toss data to the terminal
from sys import stdout
# import pi
from math import pi
# import our custom angle restraint module
import AngleRestraint
# define a little function for printing out our values of interest
def printStats():
# get the actual angle
state = simulation.context.getState(getPositions=True, getEnergy=True)
currentPositions = state.getPositions()
#print currentPositions
actualAngle = AngleRestraint.CalculateAngle(currentPositions, angleG1, angleG2, angleG3, angleG1M, angleG2M, angleG3M, 'COG')
# calculate the potential energy
energy = state.getPotentialEnergy() / (kilocalorie / mole)
# print stuff out
print "Step", i, ", Restraint angle", restraintCenterDeg, ", Actual angle", actualAngle, ", Potential energy", energy
########### do the basic setup #####################
# set the default values for the simulation
testCase = 'water'
# load the files
if testCase == 'water':
prmtop = AmberPrmtopFile('Water_Trimer_TIP3P.prmtop')
pdb = PDBFile('Water_Trimer_TIP3P.pdb')
elif testCase == 'ions':
prmtop = AmberPrmtopFile('Triple_Sodium.prmtop')
pdb = PDBFile('Triple_Sodium.pdb')
else:
print 'could not understand test case, using the water trimer'
testCase = 'water'
prmtop = AmberPrmtopFile('Water_Trimer_TIP3P.prmtop')
pdb = PDBFile('Water_Trimer_TIP3P.pdb')
# create the system
system = prmtop.createSystem(constraints=HBonds, nonbondedMethod=NoCutoff)
############ setup the angular restraint ############
# set the atom indexes involved
if testCase == 'water':
angleG1 = [0, 1, 2]
angleG2 = [3, 4, 5]
angleG3 = [6, 7, 8]
elif testCase == 'ions':
angleG1 = [0]
angleG2 = [1]
angleG3 = [2]
# set the masses in case we do center of mass restraints
angleG1M = []
angleG2M = []
angleG3M = []
for i in range(system.getNumParticles()):
if i in angleG1:
angleG1M.append(system.getParticleMass(i))
if i in angleG2:
angleG2M.append(system.getParticleMass(i))
if i in angleG3:
angleG3M.append(system.getParticleMass(i))
# setup the restraint force
angleParticles = angleG1 + angleG2 + angleG3
angleNumParticles = len(angleParticles)
restraintForceKcalMolDeg = 1.0
restraintForce = restraintForceKcalMolDeg * 4.184 * 57.2958**2
restraintCenterDeg = 180.0
restraintCenter = (restraintCenterDeg / 360.0) * (2.0 * pi)
restraintCenterType = 'COG'
restraintCalculationType = 'full'
angleRestraintString = AngleRestraint.CreateAngleRestraintString(angleG1, angleG2, angleG3, angleG1M, angleG2M, angleG3M, restraintCenterType, restraintCalculationType)
angleRestraintForce = CustomCompoundBondForce(angleNumParticles, angleRestraintString)
params = []
angleRestraintForce.addGlobalParameter('angleK', restraintForce)
angleRestraintForce.addGlobalParameter('theta0', restraintCenter)
CustomCompoundBondForce.addBond(angleRestraintForce, angleParticles, params)
system.addForce(angleRestraintForce)
print 'Successfully created an angle restraint at', restraintCenterDeg, 'degrees. With forceConstant', restraintForce, 'kJ/mol*rad^2'
#print angleRestraintString
############# setup the system #################
# define some values
timeStep = 0.5 * femtosecond # set to 0.001 in order to avoid NaN, but get odd behavior with rapid internal rotations of the waters
outputFrequency = 10
totalTimeSteps = 1000
temperature = 100.0 * kelvin
# add the integrator
integrator = LangevinIntegrator(temperature, 1/picosecond, timeStep)
# finish setting up the simulation
simulation = Simulation(prmtop.topology, system, integrator)
simulation.context.setPositions(pdb.positions)
# energy minimize
simulation.minimizeEnergy()
# add a DCD reporter to we can see any odd stuff
simulation.reporters.append(DCDReporter('trajectory.dcd', outputFrequency))
# run the simulation in blocks, stopping to output angle and energy information, and changing the angle
print '\n'
for i in range(0, int(totalTimeSteps / outputFrequency)*10):
printStats()
simulation.step(outputFrequency)
Code: Select all
# module to hold the function to make our fairly complicated angular restraint string expression for CustomCompoundBondForce in OpenMM
# import the OpenMM stuff
from simtk.openmm.app import *
from simtk.openmm import *
from simtk.unit import *
# import the math library
import math
# function to create the actual restraint string
def CreateAngleRestraintString(g1,g2,g3,g1m,g2m,g3m,option1,option2):
if option1 == 'COG':
print 'Creating Centre of Geometry Angle Restraint Energy Expression...'
# setup our common characters for string insertion
x = 'x'
y = 'y'
z = 'z'
# get the center of the first point
x1 = '('
y1 = '('
z1 = '('
for i in range(0,len(g1)):
if i != 0:
x1 += ' + '
y1 += ' + '
z1 += ' + '
k = str(i+1)
x1 += x + k
y1 += y + k
z1 += z + k
x1 += ')'
y1 += ')'
z1 += ')'
M1 = str(1.0 / float(len(g1)))
# get the center of the second point
x2 = '('
y2 = '('
z2 = '('
for i in range(len(g1),len(g1)+len(g2)):
if i != len(g1):
x2 += ' + '
y2 += ' + '
z2 += ' + '
k = str(i+1)
x2 += x + k
y2 += y + k
z2 += z + k
x2 += ')'
y2 += ')'
z2 += ')'
M2 = str(1.0 / float(len(g2)))
# get the center of the third point
x3 = '('
y3 = '('
z3 = '('
for i in range(len(g1)+len(g2),len(g1)+len(g2)+len(g3)):
if i != len(g1)+len(g2):
x3 += ' + '
y3 += ' + '
z3 += ' + '
k = str(i+1)
x3 += x + k
y3 += y + k
z3 += z + k
x3 += ')'
y3 += ')'
z3 += ')'
M3 = str(1.0 / float(len(g3)))
elif option1 == 'COM':
print 'Creating Centre of Mass Angle Restraint Energy Expression...'
# setup our common characters for string insertion
x = 'x'
y = 'y'
z = 'z'
# get the center of the first point
x1 = '('
y1 = '('
z1 = '('
for i in range(0,len(g1)):
if i != 0:
x1 += ' + '
y1 += ' + '
z1 += ' + '
M1 += ' + '
k = str(i+1)
x1 += str(g1m[i].value_in_unit(dalton)) + '*' + x + k
y1 += str(g1m[i].value_in_unit(dalton)) + '*' + y + k
z1 += str(g1m[i].value_in_unit(dalton)) + '*' + z + k
x1 += ')'
y1 += ')'
z1 += ')'
# get the center of the second point
x2 = '('
y2 = '('
z2 = '('
for i in range(len(g1),len(g1)+len(g2)):
if i != len(g1):
x2 += ' + '
y2 += ' + '
z2 += ' + '
k = str(i+1)
x2 += str(g2m[i - len(g1)].value_in_unit(dalton)) + '*' + x + k
y2 += str(g2m[i - len(g1)].value_in_unit(dalton)) + '*' + y + k
z2 += str(g2m[i - len(g1)].value_in_unit(dalton)) + '*' + z + k
x2 += ')'
y2 += ')'
z2 += ')'
# get the center of the third point
x3 = '('
y3 = '('
z3 = '('
for i in range(len(g1)+len(g2),len(g1)+len(g2)+len(g3)):
if i != len(g1)+len(g2):
x3 += ' + '
y3 += ' + '
z3 += ' + '
k = str(i+1)
x3 += str(g3m[i - (len(g1) + len(g2))].value_in_unit(dalton)) + '*' + x + k
y3 += str(g3m[i - (len(g1) + len(g2))].value_in_unit(dalton)) + '*' + y + k
z3 += str(g3m[i - (len(g1) + len(g2))].value_in_unit(dalton)) + '*' + z + k
x3 += ')'
y3 += ')'
z3 += ')'
# sum up the masses so we can be more efficient with an invert and multiply
M1Sum = 0.0
M2Sum = 0.0
M3Sum = 0.0
for mass in g1m:
M1Sum += mass.value_in_unit(dalton)
for mass in g2m:
M2Sum += mass.value_in_unit(dalton)
for mass in g3m:
M3Sum += mass.value_in_unit(dalton)
M1 = '(' + str(1.0 / M1Sum) + ')'
M2 = '(' + str(1.0 / M2Sum) + ')'
M3 = '(' + str(1.0 / M3Sum) + ')'
# the remainder of the setup now uses the COM or COG generated centers
# calculate the weighted the centers
x1 = '(' + x1 + '*' + M1 + ')'
y1 = '(' + y1 + '*' + M1 + ')'
z1 = '(' + z1 + '*' + M1 + ')'
x2 = '(' + x2 + '*' + M2 + ')'
y2 = '(' + y2 + '*' + M2 + ')'
z2 = '(' + z2 + '*' + M2 + ')'
x3 = '(' + x3 + '*' + M3 + ')'
y3 = '(' + y3 + '*' + M3 + ')'
z3 = '(' + z3 + '*' + M3 + ')'
# now do the calculations using these new centers as the points
# calculate the first vector
v0x = '(' + x2 + '-' + x1 + ')'
v0y = '(' + y2 + '-' + y1 + ')'
v0z = '(' + z2 + '-' + z1 + ')'
# calculate the second vector
v1x = '(' + x2 + '-' + x3 + ')'
v1y = '(' + y2 + '-' + y3 + ')'
v1z = '(' + z2 + '-' + z3 + ')'
# calculate the normalized magnitude of the vectors and get the inverse of their product
r0 = 'sqrt(' + v0x + '*' + v0x + '+' + v0y + '*' + v0y + '+' + v0z + '*' + v0z + ')'
r1 = 'sqrt(' + v1x + '*' + v1x + '+' + v1y + '*' + v1y + '+' + v1z + '*' + v1z + ')'
r0r1 = '(1.0 / (' + r0 + '*' + r1 + '))'
# calculate the dot product of the vectors
dot = '(' + v0x + '*' + v1x + '+' + v0y + '*' + v1y + '+' + v0z + '*' + v1z + ')'
# calculate the cosine of the angle between the vectors, but adjust for the vector magnitude
cosine = '(' + dot + '*' + r0r1 + ')'
# calculate theta using arccosine, which will be good for all values not near 0 or 180
#angle1 = '(acos(' + cosine + '))'
angle1 = '(acos(min(max(' + cosine + ', -1), 1)))'
#angle1 = '(acos(min(max(' + cosine + ', -0.99), 0.99)))' # uncomment this line and the similar one below to make it work properly
# do another calculation of theta using the asin method
# calculate the cross product of the vectors
c0x = '((' + v0z + '*' + v1y + ')' + ' - ' + '(' + v0y + '*' + v1z + '))'
c0y = '((' + v0x + '*' + v1z + ')' + ' - ' + '(' + v0z + '*' + v1x + '))'
c0z = '((' + v0y + '*' + v1x + ')' + ' - ' + '(' + v0x + '*' + v1y + '))'
# square each vector and multliply the products
scale = '((' + v0x + '*' + v0x + '+' + v0y + '*' + v0y + '+' + v0z + '*' + v0z + ')' + '*' + '(' + v1x + '*' + v1x + '+' + v1y + '*' + v1y + '+' + v1z + '*' + v1z + '))'
# square the cross product
squareCross = '(' + c0x + '*' + c0x + '+' + c0y + '*' + c0y + '+' + c0z + '*' + c0z + ')'
# calculate theta using arcsine (from P.Eastman's code)
#angle2raw = 'asin(' + 'sqrt(' + squareCross + '/' + scale + '))'
angle2raw = 'asin(min(max(' + '(sqrt(' + squareCross + '/' + scale + ')), -1), 1))'
#angle2raw = 'asin(min(max(' + '(sqrt(' + squareCross + '/' + scale + ')), -0.99), 0.99))' # uncomment this line and the similar one above to make it work properly
# fix the angle to be in the right range (from P.Eastman's code)
angle2 = '((step(' + cosine + ')' + '*' + angle2raw + ')' + '+' + '((' + '1' + '-' + 'step(' + cosine + '))' + '*' + '(' + '3.14159265' + '-' + angle2raw + ')' + '))'
# use a step function to select the appropriate method for calculating the angle. step(x) returns 0 if x is less than 0, otherwise 1
select = '(step(abs(' + cosine + ') - 0.99))' # this will return 1 when cosine needs to be calculated by the asin method, 0 otherwise
angle = '((' + select + '*' + angle2 + ')' + '+' + '((' + '1' + '-' + select + ')' + '*' + angle1 + '))'
# allow the user to select the method being used
#if option2 == 'full':
# return_string = '0.5 * ' + 'angleK' + ' * ' + '(' + 'theta0' + ' - ' + angle + ')^2'
#elif option2 == 'acos':
# return_string = '0.5 * ' + 'angleK' + ' * ' + '(' + 'theta0' + ' - ' + angle1 + ')^2'
#elif option2 == 'asin':
# return_string = '0.5 * ' + 'angleK' + ' * ' + '(' + 'theta0' + ' - ' + angle2 + ')^2'
return_string = '0.5 * ' + 'angleK' + ' * ' + '(' + 'theta0' + ' - ' + angle + ')^2'
return return_string
# function to calculate the current angle using the same inputs as the restraint string
def CalculateAngle(Positions, g1, g2, g3, g1m, g2m, g3m, option):
x1 = 0.0
y1 = 0.0
z1 = 0.0
x2 = 0.0
y2 = 0.0
z2 = 0.0
x3 = 0.0
y3 = 0.0
z3 = 0.0
M1 = 0.0
M2 = 0.0
M3 = 0.0
try:
if option == 'COM':
for count, index in enumerate(g1):
x1 += Positions[index][0].value_in_unit(nanometer) * g1m[count].value_in_unit(dalton)
y1 += Positions[index][1].value_in_unit(nanometer) * g1m[count].value_in_unit(dalton)
z1 += Positions[index][2].value_in_unit(nanometer) * g1m[count].value_in_unit(dalton)
M1 += g1m[count].value_in_unit(dalton)
for count, index in enumerate(g2):
x2 += Positions[index][0].value_in_unit(nanometer) * g2m[count].value_in_unit(dalton)
y2 += Positions[index][1].value_in_unit(nanometer) * g2m[count].value_in_unit(dalton)
z2 += Positions[index][2].value_in_unit(nanometer) * g2m[count].value_in_unit(dalton)
M2 += g2m[count].value_in_unit(dalton)
for count, index in enumerate(g3):
x3 += Positions[index][0].value_in_unit(nanometer) * g3m[count].value_in_unit(dalton)
y3 += Positions[index][1].value_in_unit(nanometer) * g3m[count].value_in_unit(dalton)
z3 += Positions[index][2].value_in_unit(nanometer) * g3m[count].value_in_unit(dalton)
M3 += g3m[count].value_in_unit(dalton)
elif option == 'COG':
for count, index in enumerate(g1):
x1 += Positions[index][0].value_in_unit(nanometer)
y1 += Positions[index][1].value_in_unit(nanometer)
z1 += Positions[index][2].value_in_unit(nanometer)
M1 += 1.0
for count, index in enumerate(g2):
x2 += Positions[index][0].value_in_unit(nanometer)
y2 += Positions[index][1].value_in_unit(nanometer)
z2 += Positions[index][2].value_in_unit(nanometer)
M2 += 1.0
for count, index in enumerate(g3):
x3 += Positions[index][0].value_in_unit(nanometer)
y3 += Positions[index][1].value_in_unit(nanometer)
z3 += Positions[index][2].value_in_unit(nanometer)
M3 += 1.0
except AttributeError:
if option == 'COM':
for count, index in enumerate(g1):
x1 += Positions[index][0] * g1m[count].value_in_unit(dalton)
y1 += Positions[index][1] * g1m[count].value_in_unit(dalton)
z1 += Positions[index][2] * g1m[count].value_in_unit(dalton)
M1 += g1m[count].value_in_unit(dalton)
for count, index in enumerate(g2):
x2 += Positions[index][0] * g2m[count].value_in_unit(dalton)
y2 += Positions[index][1] * g2m[count].value_in_unit(dalton)
z2 += Positions[index][2] * g2m[count].value_in_unit(dalton)
M2 += g2m[count].value_in_unit(dalton)
for count, index in enumerate(g3):
x3 += Positions[index][0] * g3m[count].value_in_unit(dalton)
y3 += Positions[index][1] * g3m[count].value_in_unit(dalton)
z3 += Positions[index][2] * g3m[count].value_in_unit(dalton)
M3 += g3m[count].value_in_unit(dalton)
elif option == 'COG':
for count, index in enumerate(g1):
x1 += Positions[index][0]
y1 += Positions[index][1]
z1 += Positions[index][2]
M1 += 1.0
for count, index in enumerate(g2):
x2 += Positions[index][0]
y2 += Positions[index][1]
z2 += Positions[index][2]
M2 += 1.0
for count, index in enumerate(g3):
x3 += Positions[index][0]
y3 += Positions[index][1]
z3 += Positions[index][2]
M3 += 1.0
# calculate the geometry weighted centers
x1 = x1 / M1
y1 = y1 / M1
z1 = z1 / M1
x2 = x2 / M2
y2 = y2 / M2
z2 = z2 / M2
x3 = x3 / M3
y3 = y3 / M3
z3 = z3 / M3
# calculate the two vectors that make up the angle
# calculate the first vector
v0x = x2 - x1
v0y = y2 - y1
v0z = z2 - z1
# calculate the second vector
v1x = x2 - x3
v1y = y2 - y3
v1z = z2 - z3
# calculate the normalized magnitude of the vectors and get the inverse of their product
r0 = math.sqrt( v0x * v0x + v0y * v0y + v0z * v0z )
r1 = math.sqrt( v1x * v1x + v1y * v1y + v1z * v1z )
r0r1 = (1.0 / ( r0 * r1 ))
# calculate the dot product of the vectors
dot = v0x * v1x + v0y * v1y + v0z * v1z
# calculate the cosine of the angle between the vectors
cosine = dot * r0r1
# calculate theta
angle = math.acos( cosine )
angle = angle * (360.0 / (2.0 * math.pi))
return angle
- Attachments
-
- AngleTestCase.tar.bz2
- (7.44 KiB) Downloaded 48 times
-
Peter Eastman
- Posts: 2588
- Joined: Thu Aug 09, 2007 1:25 pm
Re: adding functions to CustomCompoundBondForce
It looks like the problem is actually in the derivatives. d(asin(x))/dx = 1/sqrt(1-x^2), which is infinite when x is 1 or -1. So computing the energy works fine, but computing the force ends up multiplying 0*infinity, which is nan.
Peter
Peter
-
Aron Broom
- Posts: 54
- Joined: Tue Mar 13, 2012 11:33 am
Re: adding functions to CustomCompoundBondForce
I see, good to get to the bottom of it, and understand why it required some dynamics to see the problem.