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

Konrad Hinsen research at khinsen.fastmail.net
Mon May 17 09:19:14 UTC 2010


On 14.05.2010, at 00:45, bfaucher at ualberta.ca wrote:

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

That should work if the precision of finite-difference gradients is good enough for your application. However, don't forget to reset the position of each atom back to its original value:

> 	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)
				a.setPosition(x)
> 			gradients = ParticleVector(self.universe)
> 			gradients[a]=Vector(grad[0], grad[1], grad[2])
> 		results['gradients'] = gradients

Konrad.
--
---------------------------------------------------------------------
Konrad Hinsen
Centre de Biophysique Moléculaire, CNRS Orléans
Synchrotron Soleil - Division Expériences
Saint Aubin - BP 48
91192 Gif sur Yvette Cedex, France
Tel. +33-1 69 35 97 15
E-Mail: research at khinsen dot fastmail dot net
---------------------------------------------------------------------










More information about the mmtk mailing list