adding functions to CustomCompoundBondForce

[Aron Broom]

I've been setting up angle and dihedral center of mass or geometry restraints using the CustomCompoundBondForce method.

In general it works well, but there are a few issues which I'm hoping are solvable.

The big issue seems to be instability near 0 or 180 degree angles (or -180.0 for dihedrals), which I imagine at least in the case of the angle results from having to use arccosine to calculate the angle. The result of this instability seems to be that some atoms which I suppose are particulary close to 0 or 180 receive massive forces.

To solve the above I tried som min/max statement combos, but it doesn't really help, and I think the simple solution would be to use arctan2, which can accomplish the same goal as arccos but without the ill-definition near 0 and 180. I checked that CUDA and OpenCL both support the atan2 function, and thought to add it into the code. Following the pattern for the CudaExpressionUtilities.cpp was simple enough, and similarly for the lepton/Operations.h, but when I got to adding it into lepton/Operations.cpp, I ran up against a wall in that I don't even partially grasp what is happening with the expression trees.

I noted also that there was a custom expression option (which I take is seperate from the tabulated values option).

Anyway, I'm wondering if anyone has come across similar problems and found a solution, or can recommend how to use the custom expression or how to add things with these expression trees. Also, if there is a particular reason atan2 couldn't be added, that would be good to know also.

Thanks

~Aron

[Peter Eastman]

Hi Aron,

Computing angles in cases like this is a bit tricky. For example, here's the code that CustomCompoundBondForce uses to evaluate a dihedral term:

Code: Select all

/**
 * Compute the angle between two vectors.  The w component of each vector should contain the squared magnitude.
 */
real ccb_computeAngle(real4 vec1, real4 vec2) {
    real dotProduct = vec1.x*vec2.x + vec1.y*vec2.y + vec1.z*vec2.z;
    real cosine = dotProduct*RSQRT(vec1.w*vec2.w);
    real angle;
    if (cosine > 0.99f || cosine < -0.99f) {
        // We're close to the singularity in acos(), so take the cross product and use asin() instead.

        real4 crossProduct = cross(vec1, vec2);
        real scale = vec1.w*vec2.w;
        angle = asin(SQRT(dot(crossProduct, crossProduct)/scale));
        if (cosine < 0)
            angle = M_PI-angle;
    }
    else
       angle = acos(cosine);
    return angle;
}

I'm not sure how you would reformulate something like this in terms of atan(), and in any case, I don't think atan2() has any significant numerical advantages over atan(). The main reason for using it is to preserve sign information so you can tell which quadrant you're in.

What is the custom expression you're currently using?

Peter

[Aron Broom]

Hi Peter,

I saw that code or something very similar as part of the torsionalForces code I think. I actually modelled my angle and dihedral expressions off of what you guys had in that platform code, but maybe I'm missing some of the functionality of the CustomCompoundBondForce.

In particular, my version of what you have is one giant string, and instead of if statements, I tried to use a series of min/max nests to get the right result (which you can imagine leads to an insanely large and slow string). Is it possible to use normal code, and in particular, if statements within the CustomCompoundBondForce?

Also, I was manually calculating the cross-product, which is probably a minor slowdown aswell.

So for me atan2 seemed useful because I could avoid the min/max nests, but if there is a way to just use something like what you've given, that would be outrageously ideal by comparison.

~Aron

[Aron Broom]

I've been digging around and can't exactly find where or how I would add something similar to what you've posted. I see that in "ReferenceCustomCompoundBondIxn.cpp" there appear to be references to the built in distance, angle and dihedral options of the CustomCompoundBondForce, but it's not clear to me how they get from the expression string to there, as I can't find mention of them in places like Parer.cpp or CustomExpressionUtilities.cpp.

Anyway, what I'm currently using for a Center of Geometry dihedral restraint for instance, looks like this:

Code: Select all

def CenterOfGeometryDihedralRestraintString(g1, g2, g3, g4):	
	# where g1 contains the atom indices for the atoms in the first center of geometry

	# setup our common characters for string insertion
	x = 'x'
	y = 'y'
	z = 'z'

	# get the center of the first group
	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 += ')'

	# geometrically average the centers
	M1 = str(1.0 / float(len(g1)))
	x1 = '(' + x1 + '*' + M1 + ')'
	y1 = '(' + y1 + '*' + M1 + ')'
	z1 = '(' + z1 + '*' + M1 + ')'


	# get the center of the second group
	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 += ')'

	# geometrically average the centers
	M2 = str(1.0 / float(len(g2)))
	x2 = '(' + x2 + '*' + M2 + ')'
	y2 = '(' + y2 + '*' + M2 + ')'
	z2 = '(' + z2 + '*' + M2 + ')'


	# get the center of the third group
	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 += ')'

	# geometrically average the centers
	M3 = str(1.0 / float(len(g3)))
	x3 = '(' + x3 + '*' + M3 + ')'
	y3 = '(' + y3 + '*' + M3 + ')'
	z3 = '(' + z3 + '*' + M3 + ')'


	# get the center of the fourth group
	x4 = '('
	y4 = '('
	z4 = '('
	
	for i in range(len(g1)+len(g2)+len(g3),len(g1)+len(g2)+len(g3)+len(g4)):
		if i != len(g1)+len(g2)+len(g3):
			x4 += ' + '
			y4 += ' + '
			z4 += ' + '
		k = str(i+1)
		x4 += x + k
		y4 += y + k
		z4 += z + k
	x4 += ')'
	y4 += ')'
	z4 += ')'

	# geometrically average the centers
	M4 = str(1.0 / float(len(g4)))
	x4 = '(' + x4 + '*' + M4 + ')'
	y4 = '(' + y4 + '*' + M4 + ')'
	z4 = '(' + z4 + '*' + M4 + ')'


	# setup the calculations involving the centers we've calculated
	# get the vector unique to the first plane
	v0x = '((' + x1 + ' * ' + M1 + ') - (' + x2 + ' * ' + M2 + '))' 
	v0y = '((' + y1 + ' * ' + M1 + ') - (' + y2 + ' * ' + M2 + '))' 
	v0z = '((' + z1 + ' * ' + M1 + ') - (' + z2 + ' * ' + M2 + '))' 
	
	# get the common vector
	v1x = '((' + x3 + ' * ' + M3 + ') - (' + x2 + ' * ' + M2 + '))' 
	v1y = '((' + y3 + ' * ' + M3 + ') - (' + y2 + ' * ' + M2 + '))' 
	v1z = '((' + z3 + ' * ' + M3 + ') - (' + z2 + ' * ' + M2 + '))' 

	# get the vector unique to the second plane
	v2x = '((' + x3 + ' * ' + M3 + ') - (' + x4 + ' * ' + M4 + '))' 
	v2y = '((' + y3 + ' * ' + M3 + ') - (' + y4 + ' * ' + M4 + '))' 
	v2z = '((' + z3 + ' * ' + M3 + ') - (' + z4 + ' * ' + M4 + '))' 


	# calculate the normals to each of the two plains using the vector pairs above (would be the cross product)
	# get the normal to the first plane
	c0x = '((' + v0z + '*' + v1y + ')' + ' - ' + '(' + v0y + '*' + v1z + '))'
	c0y = '((' + v0x + '*' + v1z + ')' + ' - ' + '(' + v0z + '*' + v1x + '))'
	c0z = '((' + v0y + '*' + v1x + ')' + ' - ' + '(' + v0x + '*' + v1y + '))'
	
	# get the normal to the second plane
	c1x = '((' + v1z + '*' + v2y + ')' + ' - ' + '(' + v1y + '*' + v2z + '))'
	c1y = '((' + v1x + '*' + v2z + ')' + ' - ' + '(' + v1z + '*' + v2x + '))'
	c1z = '((' + v1y + '*' + v2x + ')' + ' - ' + '(' + v1x + '*' + v2y + '))'

	# calculate the absolute magnitude of the vectors and the inverse of their product
	r0 = 'sqrt(' + c0x + '*' + c0x + '+' + c0y + '*' + c0y + '+' + c0z + '*' + c0z + ')'
	r1 = 'sqrt(' + c1x + '*' + c1x + '+' + c1y + '*' + c1y + '+' + c1z + '*' + c1z + ')'
	r0r1 = '(1.0 / (' + r0 + '*' + r1 + '))'

	# calculate the dot product of the vectors
	dot = '(' + c0x + '*' + c1x + '+' + c0y + '*' + c1y + '+' + c0z + '*' + c1z + ')'

	# calculate the cosine of the angle
	cosine = '(' + dot + '*' + r0r1 + ')'

	# calculate the angle between the two planes using arccosine
	dihedral = 'acos(' + cosine + ')'

	# since arccosine cannot return values in the range 0:-180, need to do that part ourselves with some ugly math
	# first calculate the dot product of the first vector and the second cross product
	signmod = '(' + v0x + '*' + c1x + '+' + v0y + '*' + c1y + '+' + v0z + '*' + c1z + ')'

	# next calculate the absolute value of the above
	abssignmod = 'abs(' + signmod + ')'

	# finally the real dihedral is the the arccosine dihedral multiplied by -1 if the above is -ve, otherwise left as is
	sign = '(' + signmod + '/' + abssignmod + ')'
	realdihedral = '(' + dihedral + '*' + sign + '*' + '-1.0' + ')'

	# because torsion is periodic and continuous, we need to ensure that the smallest value is taken as the difference from the minimum
	minimumDifference = 'min(' + '(' + '(' + 'psi0' + '-' + realdihedral + ')^2' + ')' + ',' + '(' + '(' + '360.0' + '-' + 'abs(' + realdihedral + '-' + 'psi0' + ')' + ')^2' + ')' + ')'

	# finally return the string
	return_string = '0.5 * ' + 'dihedralK' + ' * ' + '(' + minimumDifference + ')^2'

As you can see, it will generate quite a massive string if one has each center of geometry composed of say ~5 atoms. The performance actually isn't that bad in general which shocks me (perhaps this is because of the expression tree stuff?). But, currently it doesn't do well near 0 or +/- 180 degrees, often blowing apart the system (works great otherwise though). I was going to try something like your arcsin trick in the code you showed, but even though a large string seems to run well once the simulation starts, the setup time can start becoming multiple minutes even approaching an hour, and so I'd love to avoid that if possible.

Also, I'm not sure if there's a better way to deal with the fact that an angle of -170 is quite close to 170. In the above I have an ugly method that just calculates the two options: either 170 - (-170) = 340, or 360 - (abs( 170 - (-170)) = 20. I couldn't really see how that was handled in the torsionForces.cu for instance.

I guess my thought was atan2 was that by preserving the quadrant information I wouldn't need the extra tricks with the sign thing, and in general by using an atan style of approach, I might not need to try and introduce the whole asin condition.

Any suggestions?

~Aron

[Peter Eastman]

Hi Aron,

The code I quoted was taken from customCompoundBond.cl. The easiest way to do an if/else in an expression is by using a step() function. Compute the angle using both methods and store them in angle1 and angle2. Then you can specify something like "angle = select*angle1 + (1-select)*angle2; select=step(abs(cosine)-0.99)".

Startup time can definitely be a problem with very complicated expressions. If that becomes unacceptable, your best bet is probably a plugin. Ultimately I'd like to create a new custom force class that's specially designed for cases like this, but that's further in the future.

It's also possible, of course, that the expression parser could be made faster. I'll look into that.

Peter

[Aron Broom]

Peter,

Thanks for the suggestion, I'll give that a try and see if it can solve the instability problem.

To be fair to the parser, the final string is about a page long in the case I mention, so I can well imagine it is a task.

~Aron

[Peter Eastman]

I did some profiling to see why it's taking so long. It turns out parsing the expression isn't the problem - that takes a fraction of a second. The expensive part is differentiating it with respect to all the different variables, and optimizing all those derivative expressions.

Peter

[Aron Broom]

Well, it does seem like despite the massive expression, the performance of the actual simulation is hardly effected, so I would say your whole method of allowing these custom expressions has worked really well. I'll just have to try and run longer simulations, or at least, see if I can keep using an existing system setup.

Thanks for looking at it

[Aron Broom]

I've found something rather odd, or perhaps I've totally lost my mind.

I tested out the step idea you suggested, but it seems to not have the intended effect.

I calculate 'angle1' using the acos method, which becomes unstable near 0 and 180, great elsewhere.
I calculate 'angle2' using the asin method I grabbed from your code snipped. It works great at 0 and 180, but near 90 becomes unstable.

Ok, so all good at this point. If I then use:

Code: Select all

select = 'step(abs(' + cosine + ')-0.99)'

as you suggested, which should return 0 when we're in the acos range and 1 when we're in the asin range, and then make the final statement for the angle:

Code: Select all

angle = '(' + select + '*' + angle2 + ')' + '+' + '((' + '1' + '-' + select + ')' + '*' + angle1 + ')'

Rather than getting the intended result of it now working everywhere, it infact now fails not only near 0 and 180, but also 90. I checked quickly if I was confused about what step returns by swapping the angle1 and angle2 in the above expression, but it has no effect, still working well anywhere but 0, 90, and 180.

I guess this potentially means that as anything in the expression becomes undefined, multiplying it by 0 doesn't necessarily result in 0? Or maybe I'm missing something more fundamental.

[Peter Eastman]

If a function is actually returning nan then multiplying it by 0 may not work. The result could still come out nan. But as long as it has a well defined finite value, this should work correctly.

Perhaps roundoff error is causing the arguments to go outside the allowed range? Try replacing "acos(x)" with "acos(min(max(x, -1), 1))". Does that help?

Peter