adding functions to CustomCompoundBondForce

The functionality of OpenMM will (eventually) include everything that one would need to run modern molecular simulation.
User avatar
Aron Broom
Posts: 54
Joined: Tue Mar 13, 2012 11:33 am

Re: adding functions to CustomCompoundBondForce

Post by Aron Broom » Wed Jul 24, 2013 1:18 pm

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.

User avatar
Peter Eastman
Posts: 2588
Joined: Thu Aug 09, 2007 1:25 pm

Re: adding functions to CustomCompoundBondForce

Post by Peter Eastman » Wed Jul 24, 2013 2:13 pm

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

User avatar
Aron Broom
Posts: 54
Joined: Tue Mar 13, 2012 11:33 am

Re: adding functions to CustomCompoundBondForce

Post by Aron Broom » Wed Jul 24, 2013 2:44 pm

Sounds like a plan. I'll give it a shot and report back on the findings.

~Aron

User avatar
Aron Broom
Posts: 54
Joined: Tue Mar 13, 2012 11:33 am

Re: adding functions to CustomCompoundBondForce

Post by Aron Broom » Thu Jul 25, 2013 2:46 pm

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.

User avatar
Peter Eastman
Posts: 2588
Joined: Thu Aug 09, 2007 1:25 pm

Re: adding functions to CustomCompoundBondForce

Post by Peter Eastman » Thu Jul 25, 2013 3:46 pm

Very strange. Can you send me your test case? I'd like to try to figure out what's going on.

Peter

User avatar
Aron Broom
Posts: 54
Joined: Tue Mar 13, 2012 11:33 am

Re: adding functions to CustomCompoundBondForce

Post by Aron Broom » Fri Jul 26, 2013 2:34 pm

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

User avatar
Peter Eastman
Posts: 2588
Joined: Thu Aug 09, 2007 1:25 pm

Re: adding functions to CustomCompoundBondForce

Post by Peter Eastman » Fri Jul 26, 2013 2:55 pm

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

User avatar
Aron Broom
Posts: 54
Joined: Tue Mar 13, 2012 11:33 am

Re: adding functions to CustomCompoundBondForce

Post by Aron Broom » Tue Jul 30, 2013 9:04 am

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.

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)
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

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

User avatar
Peter Eastman
Posts: 2588
Joined: Thu Aug 09, 2007 1:25 pm

Re: adding functions to CustomCompoundBondForce

Post by Peter Eastman » Thu Aug 01, 2013 3:41 pm

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

User avatar
Aron Broom
Posts: 54
Joined: Tue Mar 13, 2012 11:33 am

Re: adding functions to CustomCompoundBondForce

Post by Aron Broom » Thu Aug 22, 2013 8:46 am

I see, good to get to the bottom of it, and understand why it required some dynamics to see the problem.

POST REPLY