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

bfaucher at ualberta.ca bfaucher at ualberta.ca
Thu May 13 15:23:45 UTC 2010

Dr. Hinsen;

Thank you. My problem stems from MMTK's indexing its gradients on an  
atom-by-atom basis while I compute energies molecule-by-molecule*. For  
now, I will designate the "global" gradient to a single atom and set  
0s for the gradients of other atoms; hopefully this will suffice.

*I am calculating energies by extending Python with iOpenShell's Black  
Box-ish C++ library, which runs as "atomic (cartesian) coordinates in  
-> molecular PE (at the origin) out".

Thank you again!
Bryan Faucher

Quoting "Konrad Hinsen" <research at khinsen.fastmail.net>:

> On 13 May 2010, at 00:56, bfaucher at ualberta.ca wrote:
>> I am currently attempting to develop an MMTK force field which  
>> implements iOpenShell's library of Potential Energy Surfaces  
>> (http://iopenshell.usc.edu/downloads/ezpes/). By editing the  
>> included Harmonic Oscillator example I've been able to develop an  
>> ff that returns a proper EnergyTerm, but I'm having some trouble  
>> getting the derivatives into the right places. Would someone mind  
>> walking me through my "if do_gradients" block?
> I am assuming you are looking at the Python version of  
> HarmonicOscillator, right? Then this is the piece of code of interest:
>         if do_gradients:
>             gradients = ParticleVector(self.universe)
>             gradients[self.atom_index] = self.force_constant*d
>             results['gradients'] = gradients
> In fact, the only line you need to change (and perhaps replace by  
> several lines) is
>             gradients[self.atom_index] = self.force_constant*d
> This is where the calculated gradients get put in the right place.  
> "gradients" holds a 3-element-vector for each atom. For each atom,  
> you have to assign the corresponding gradient vector (an object of  
> type Scientific.Geometry.Vector) to the corresponding entry in  
> "gradients". You can use the atom object itself for indexing, e.g.:
> 	for atom in self.universe.atomList():
> 		gradients[atom] = Vector(0., 0., 0.)
> to set all gradients to zero. You can also use the atom indices, as  
> shown in the example, if that is more convenient for you.
> I hope this helps. If you need more details, it would help if you  
> would provide the interface code you have, so that we can see how  
> you do the energy calculations as a function of atom positions.
> 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