[MMTK] Building Force Fields under MMTK

[Guido Wagner]

Hi !

Basicly terms for more than one atom may be programed the same way as the 
example for the harmonic oscilator force field, you just have to pass more 
than one atom object to your term-routine. MMTK also has methods that list 
bonds bond angles and dihedrals. Using these you may build a force filed - 
class, that creates an evaluator list (this is a python list, each element 
consists of the therm-routine and the necesary data for the therm).

For example if you want to do a simple energy-calculation with your new force 
field, the Python script may look like this:

=============================

from MOMOxmlMoleculeFactory import *
from eval import *
from Numeric import *
import sys

# Universe aus MOMO-XML Datei erzeugen

universe = MOMOXMLMoleculeFactory("BENZOL.XML").universe
universe.setForceField(MOMOForceField())
print "Energy:", universe.energyGradientsAndForceConstants() 

==================

The eval-library contains my the force field class (that creates the evaluator 
list):

================================


class MOMOForceField:
    def __init__(self,parameterfile="MOMO.MPL"):
        #print "MOMOFF init..."
        self.param = params.Parameter()
        file = open(parameterfile,"r")
        for line in file.readlines():
            self.param.parse_ParamLine(line)
        file.close()
        print 



    def evaluatorTerms(self, universe, subset1, subset2,globaldata):

	Atomzahl = len(universe.atomList())
        Atomtypen = zeros((Atomzahl))
        AtomFlags = zeros((Atomzahl))
        AtomKoordinaten = zeros((Atomzahl,3),Float)
        MatrixSize = Atomzahl * ((Atomzahl+1.)/2)
        ICONMatrix = zeros(MatrixSize,Int)

        #print "Atome::",Atomzahl,"Matrix:",MatrixSize
              
        i=-1
        ObjektZahl=0
        ErstAtomInObj = []
        AtomeInObjekt = []
        ObjektdesAtoms = []

        for object in universe:

            ObjektZahl = ObjektZahl + 1
            j = i

            atomcount = 0
            erster = 1

            for atom in object.atomList():
                i=i+1
                atomcount = atomcount + 1
                ObjektdesAtoms.append(j)
                
                if erster:
                    erster = 0
                    ErstAtomInObj.append(i)

                Atomtypen[i] = object.getAtomProperty(atom,'MOMOAtomType')
                AtomKoordinaten[i] = atom.position()
                AtomFlags[i] = self.param.C.ITPFLG[Atomtypen[i]]
                #AtomFlags[i]= object.getAtomProperty(atom,'MOMO_std_atflag')

            AtomeInObjekt.append(atomcount)
            
            bonds = []
            angles = []
            dihedrals = []
            impropers = []
	   
            for bu in object.bondedUnits():

                for bond in bu.bonds:
                    # Bindungslängenterm dazu
                    bonds.append((bond.a1.index,bond.a2.index))
                    Substituentenzaehler[bond.a1.index] = 
Substituentenzaehler[bond.a1.index] + 1
                    Substituentenzaehler[bond.a2.index] = 
Substituentenzaehler[bond.a2.index] + 1
                    
                for angle in bu.bonds.bondAngles():
                    # Bindungswinkelterm dazu
                    angles.append((angle.ca.index,angle.a1.index, 
angle.a2.index))


                for angle in bu.bonds.dihedralAngles():
                    if angle.improper:
                        impropers.append((angle.a1.index, angle.a2.index, 
angle.a3.index, angle.a4.index))

                    elif not angle.improper:
                        dihedrals.append((angle.a1.index, angle.a2.index, 
angle.a3.index, angle.a4.index))
                        

        # Bindungslisten aufbauen
        ListBonds = array(bonds)
        ListAngles = array(angles)
        ListDihedrals = array(dihedrals)
        ListImpropers = array(impropers)




..........

lots of stuff to prepare and pack data skiped

........
      # a term using all atoms and bonds for atomic point charge 
calculation ....

	from MOMO_GMcharge import GMchargeTerm
        
        	ChargeTerm = 				
GMchargeTerm(universe._spec,GMparams,ListBonds,QPI,QSIG,SWI,deltaq)
        eval_list.append(ChargeTerm)



 .....

============


I hope this example was helpfull.

CU Guido


Am Sonntag, 20. November 2005 14:16 schrieb Gustavo Mercier:
> Hi!
>
>   I am interested in building a simple fully harmonic force field under
> MMTK. The purpose is two fold: 1) Learn how to build force fields under
> MMTK ; and 2) Assit me with the development of other software.
>
>   I've seen the example under the Forcefield directory (MMTK-2.5.9), but
> remain a bit confused. I am starting simple by looking at the python
> example of the harmonic potential.
>
>   It seems to me that the potential terms are defined per atom. For the two
> examples in the Forcefield directory this is adequate. But how do you do it
> when you need to specify a term that needs the position of multiple
> particles (atoms)? I've browsed through the code of MMTK and I'm getting
> lost in the hierarchy of many objects and their relationships.
>
>   Any one can shed some light on how to build a potential term that depends
> in the position of multiple atoms and done within the MMTK scheme?
>
>   Thanks for your help!
>   Gustavo Mercier, Jr MD,PhD
>   gamercier at yahoo.com
>   gustavom at baylorhealth.edu
>
>
>
> --
> Gustavo A. Mercier, Jr. MD,PhD
> Baylor University Medical Center
> Radiology
> American Radiology Associates
> 712 N. Washington, Suite 101
> Dallas, TX 75246
> 214-826-8822
> 214-826-9792 fax
> gamercier at yahoo.com