[MMTK] Implementing Gradients and Force Constants for a Molecular Force Field

bfaucher at ualberta.ca bfaucher at ualberta.ca
Thu May 13 22:45:24 UTC 2010


Update:

I had clearly been over complicating my problem; my apologies. I've  
settled on the very simple finite difference code:

	if do_gradients:
		atoms = self.universe.atomList()
		for a in atoms:
			grad=[]
			for v in [ex, ey, ez]:
				delta = 0.001
				x = a.position()
				a.setPosition(x+delta*v)
				eplus1 = self.universe.energy()
				a.setPosition(x+2.*delta*v)
				eplus2 = self.universe.energy()
				a.setPosition(x-delta*v)
				emin1 = self.universe.energy()
				a.setPosition(x-2.*delta*v)
				emin2 = self.universe.energy()
				grad.append(((1.0/12.0)*emin2-(2.0/3.0)*emin1+(2.0/3.0)*eplus1-(1.0/12.0)*eplus2)/delta)
			gradients = ParticleVector(self.universe)
			gradients[a]=Vector(grad[0], grad[1], grad[2])
		results['gradients'] = gradients

Which appears to be working for my test geometries. Thanks again!

Bryan Faucher



More information about the mmtk mailing list