[MMTK] Reading of Amber Parameters

Berit Hinnemann Berit.Hinnemann@fysik.dtu.dk
Tue, 06 Nov 2001 15:43:23 +0100


This is a multi-part message in MIME format.
--------------7229F9E3D257C0C5BE30B1E6
Content-Type: text/plain; charset=us-ascii
Content-Transfer-Encoding: 7bit

Hi Konrad,

I have a question concerning the reading of the Amber Parameters, which
is done, when the object Amber94ForceField() is created. If I've
understood the code correctly, the reading of the parameters occurs by
creating an object AmberData.AmberParameters(), where the parameters are
read from the main file and possible modification files. The parameters
are stored in dictionaries, which are attributes of AmberParameters().

Where does the caching of the parameters from the main file
(amber_parm94) take place exactly? As far as I can see, the storage of
the parameters in the dictionaries AmberParameters.atom_types, .bonds,
.bond_angles, etc. does not distinguish between parameters coming from
the main file or from the modification files, and after that, only those
dictionaries are processed. So I don't understand, where the copies of
the parameters of the main file are taken.

The reason for my question is that I've modified AmberData.py and
AmberForceField.py such that when the Amber force field is created, it
also takes a dictionary of Amber parameters as input (I'd like to change
the parameters frequently while running MMTK and with file I/O it takes
too long time). That works fine, only, that when I create a new universe
and force field, the parameters of the list are not reread, and I'd like
them to be reread each time I create a new universe. I've attached the
modified files to this mail.

Thanks a lot,
Berit


-- 
Berit Hinnemann
Center for Atomic-Scale Materials Physics (CAMP)
Dept. of Physics, Technical University of Denmark
Building 307, DK-2800 Lyngby, Denmark 
Phone: +45 4525 3209, Fax: +45 4593 2399
--------------7229F9E3D257C0C5BE30B1E6
Content-Type: text/plain; charset=us-ascii;
 name="AmberListData.py"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="AmberListData.py"

# This module handles input and application of Amber parameter files.
# This is a module modified from AmberData.py and besides reading the
# parameter files are the parameters delivered by the list also set.
#
# Written by Konrad Hinsen
# last revision: 1999-9-13
#
# Modified by Berit Hinnemann
# last revision: 2001-11-06

_undocumented = 1

from MMTK import Units
from Scientific.IO.FortranFormat import FortranFormat, FortranLine
from Scientific.IO.TextFile import TextFile
from Scientific.DictWithDefault import DictWithDefault
import os, string

amber_energy_unit = Units.kcal/Units.mol
amber_length_unit = Units.Ang
amber_angle_unit = Units.deg
amber_fc_angle_unit = Units.rad

#
# The force field parameter file
#
class AmberListParameters:

    def __init__(self, amber_params, filename, modifications=[]):
	file = TextFile(filename)
	title = file.readline()[:-1]

	self.atom_types = DictWithDefault(None)
        self._readAtomTypes(file)

	format = FortranFormat('20(A2,2X)')
	done = 0
	while not done:
	    l = FortranLine(file.readline()[:-1], format)
	    for entry in l:
		name = _normalizeName(entry)
		if len(name) == 0:
		    done = 1
		    break
		try: # ignore errors for now
		    self.atom_types[name].hydrophylic = 1
		except: pass

	self.bonds = {}
        self._readBondParameters(file)

	self.bond_angles = {}
        self._readAngleParameters(file)

	self.dihedrals = {}
	self.dihedrals_2 = {}
        self._readDihedralParameters(file)

	self.impropers = {}
	self.impropers_1 = {}
	self.impropers_2 = {}
        self._readImproperParameters(file)

	self.hbonds = {}
        self._readHbondParameters(file)

	self.lj_equivalent = {}
	format = FortranFormat('20(A2,2X)')
	while 1:
	    l = FortranLine(file.readline()[:-1], format)
	    if l.isBlank(): break
	    name1 = _normalizeName(l[0])
	    for s in l[1:]:
		name2 = _normalizeName(s)
		self.lj_equivalent[name2] = name1

	self.ljpar_sets = {}
	while 1:
	    l = FortranLine(file.readline()[:-1], 'A4,6X,A2')
	    if l[0] == 'END ': break
	    set_name = _normalizeName(l[0])
	    ljpar_set = AmberLJParameterSet(set_name, l[1])
	    self.ljpar_sets[set_name] = ljpar_set
            self._readLJParameters(file, ljpar_set)

	file.close()

        for mod, ljname in modifications:
            file = TextFile(mod)
            title = file.readline()[:-1]
            blank = file.readline()[:-1]
            while 1:
                keyword = file.readline()
                if not keyword: break
                keyword = string.strip(keyword)[:4]
                if keyword == 'MASS':
                    self._readAtomTypes(file)
                elif keyword == 'BOND':
                    self._readBondParameters(file)
                elif keyword == 'ANGL':
                    self._readAngleParameters(file)
                elif keyword == 'DIHE':
                    self._readDihedralParameters(file)
                elif keyword == 'IMPR':
                    self._readImproperParameters(file)
                elif keyword == 'HBON':
                    self._readHBondParameters(file)
                elif keyword == 'NONB':
                    self._readLJParameters(file, self.ljpar_sets[ljname])
        #Added: set the parameters given to AmberListParameters
        self._setAtomTypes(amber_params['atom_types'])
        self._setBondParameters(amber_params['bonds'])
        self._setAngleParameters(amber_params['angles'])
        self._setDihedralParameters(amber_params['dihedrals'])
        self._setImproperParameters(amber_params['impropers'])
        self._setLJParameters(amber_params['lennard_jones'], self.ljpar_sets[ljname])

    def _readAtomTypes(self, file):
	format = FortranFormat('A2,1X,F10.2')
	while 1:
	    l = FortranLine(file.readline()[:-1], format)
	    if l.isBlank(): break
	    name = _normalizeName(l[0])
	    self.atom_types[name] = AmberAtomType(name, l[1])

    def _readBondParameters(self, file):
	format = FortranFormat('A2,1X,A2,2F10.2')
	while 1:
	    l = FortranLine(file.readline()[:-1], format)
	    if l.isBlank(): break
	    name1 = _normalizeName(l[0])
	    name2 = _normalizeName(l[1])
	    name1, name2 = _sort(name1, name2)
	    self.bonds[(name1, name2)] = \
			       AmberBondParameters(self.atom_types[name1],
						   self.atom_types[name2],
						   l[2], l[3])

    def _readAngleParameters(self, file):
	format = FortranFormat('A2,1X,A2,1X,A2,2F10.2')
	while 1:
	    l = FortranLine(file.readline()[:-1], format)
	    if l.isBlank(): break
	    name1 = _normalizeName(l[0])
	    name2 = _normalizeName(l[1])
	    name3 = _normalizeName(l[2])
	    name1, name3 = _sort(name1, name3)
	    self.bond_angles[(name1, name2, name3)] = \
		   AmberBondAngleParameters(self.atom_types[name1],
					    self.atom_types[name2],
					    self.atom_types[name3],
					    l[3], l[4])

    def _readDihedralParameters(self, file):
	format = FortranFormat('A2,1X,A2,1X,A2,1X,A2,I4,3F15.2')
	append = None
	while 1:
	    l = FortranLine(file.readline()[:-1], format)
	    if l.isBlank(): break
	    name1 = _normalizeName(l[0])
	    name2 = _normalizeName(l[1])
	    name3 = _normalizeName(l[2])
	    name4 = _normalizeName(l[3])
	    name1, name2, name3, name4 = _sort4(name1, name2, name3, name4)
	    if append is not None:
		append.addTerm(l[4], l[5], l[6], l[7])
		if l[7] >= 0: append = None
	    else:
		p = AmberDihedralParameters(self.atom_types[name1],
					    self.atom_types[name2],
					    self.atom_types[name3],
					    self.atom_types[name4],
					    l[4], l[5], l[6], l[7])
		if name1 == 'X' and name4 == 'X':
		    self.dihedrals_2[(name2, name3)] = p
		else:
		    self.dihedrals[(name1, name2, name3, name4)] = p
		if l[7] < 0: append = p

    def _readImproperParameters(self, file):
	format = FortranFormat('A2,1X,A2,1X,A2,1X,A2,I4,3F15.2')
	while 1:
	    l = FortranLine(file.readline()[:-1], format)
	    if l.isBlank(): break
	    name1 = _normalizeName(l[0])
	    name2 = _normalizeName(l[1])
	    name3 = _normalizeName(l[2])
	    name4 = _normalizeName(l[3])
	    if name1 == 'X':
		if name2 == 'X':
		    name1 = name3
		    name2 = name4
		    name3 = 'X'
		    name4 = 'X'
		else:
		    name1 = name3
		    name2, name3 = _sort(name2, name4)
		    name4 = 'X'
	    else:
		name1, name2, name3, name4 = \
		       (name3, ) + _sort3(name1, name2, name4)
	    p = AmberImproperParameters(self.atom_types[name1],
					self.atom_types[name2],
					self.atom_types[name3],
					self.atom_types[name4],
					l[4], l[5], l[6], l[7])
	    if name4 == 'X':
		if name3 == 'X':
		    self.impropers_2[(name1, name2)] = p
		else:
		    self.impropers_1[(name1, name2, name3)] = p
	    else:
		self.impropers[(name1, name2, name3, name4)] = p

    def _readHbondParameters(self, file):
	format = FortranFormat('2X,A2,2X,A2,2X,2F10.2')
	while 1:
	    l = FortranLine(file.readline()[:-1], format)
	    if l.isBlank(): break
	    name1 = _normalizeName(l[0])
	    name2 = _normalizeName(l[1])
	    name1, name2 = _sort(name1, name2)
	    self.hbonds[(name1, name2)] = \
				AmberHbondParameters(self.atom_types[name1],
						     self.atom_types[name2],
						     l[2], l[3])

    def _readLJParameters(self, file, ljpar_set):
	format = FortranFormat('2X,A2,6X,3F10.6')
        while 1:
            l = FortranLine(file.readline()[:-1], format)
            if l.isBlank(): break
            name = _normalizeName(l[0])
            ljpar_set.addEntry(name, l[1], l[2], l[3])

    #Those are the functions added by BH
    def _setAtomTypes(self, atom_dict):
        for i in atom_dict.keys():
            name = _normalizeName(i)
            mass = atom_dict[i]
            self.atom_types[name] = AmberAtomType(name, mass)

    def _setBondParameters(self, bond_dict):
        for i in bond_dict.keys():
            name1 = _normalizeName(string.split(i, '-')[0])
            name2 = _normalizeName(string.split(i, '-')[1])
            name1, name2 = _sort(name1, name2)
            force = bond_dict[i][0]
            equi = bond_dict[i][1]
	    self.bonds[(name1, name2)] = \
			       AmberBondParameters(self.atom_types[name1],
						   self.atom_types[name2],
						   force, equi)

    def _setAngleParameters(self, angle_dict):
        for i in angle_dict.keys():
            name1 = _normalizeName(string.split(i, '-')[0])
	    name2 = _normalizeName(string.split(i, '-')[1])
	    name3 = _normalizeName(string.split(i, '-')[2])
	    name1, name3 = _sort(name1, name3)
            force = angle_dict[i][0]
            equi = angle_dict[i][1]
	    self.bond_angles[(name1, name2, name3)] = \
		   AmberBondAngleParameters(self.atom_types[name1],
					    self.atom_types[name2],
					    self.atom_types[name3],
					    force, equi)

    def _setDihedralParameters(self, dihedral_dict):
        for i in dihedral_dict.keys():
            name1 = _normalizeName(string.split(i, '-')[0])
	    name2 = _normalizeName(string.split(i, '-')[1])
	    name3 = _normalizeName(string.split(i, '-')[2])
            name4 = _normalizeName(string.split(i, '-')[3])
            name1, name2, name3, name4 = _sort4(name1, name2, name3, name4)
            par1 = dihedral_dict[i][0]
            par2 = dihedral_dict[i][1]
            par3 = dihedral_dict[i][2]
            par4 = dihedral_dict[i][3]
            p = AmberDihedralParameters(self.atom_types[name1],
                                        self.atom_types[name2],
                                        self.atom_types[name3],
                                        self.atom_types[name4],
                                        par1, par2, par3, par4)
            if name1 == 'X' and name4 == 'X':
                self.dihedrals_2[(name2, name3)] = p
            else:
                self.dihedrals[(name1, name2, name3, name4)] = p
            ct= 4
            while 1:
                if par4 >= 0: break
                else:
                    par1 = dihedral_dict[i][ct+0]
                    par2 = dihedral_dict[i][ct+1]
                    par3 = dihedral_dict[i][ct+2]
                    par4 = dihedral_dict[i][ct+3]
                    p.addTerm(par1, par2, par3, par4)
                    ct = ct + 4

    def _setImproperParameters(self, improper_dict):
        for i in improper_dict.keys():
            name1 = _normalizeName(string.split(i, '-')[0])
	    name2 = _normalizeName(string.split(i, '-')[1])
	    name3 = _normalizeName(string.split(i, '-')[2])
	    name4 = _normalizeName(string.split(i, '-')[3])
            par1 = improper_dict[i][0]
            par2 = improper_dict[i][1]
            par3 = improper_dict[i][2]
            par4 = improper_dict[i][3]
	    if name1 == 'X':
		if name2 == 'X':
		    name1 = name3
		    name2 = name4
		    name3 = 'X'
		    name4 = 'X'
		else:
		    name1 = name3
		    name2, name3 = _sort(name2, name4)
		    name4 = 'X'
	    else:
		name1, name2, name3, name4 = \
		       (name3, ) + _sort3(name1, name2, name4)
	    p = AmberImproperParameters(self.atom_types[name1],
					self.atom_types[name2],
					self.atom_types[name3],
					self.atom_types[name4],
					par1, par2, par3, par4)
	    if name4 == 'X':
		if name3 == 'X':
		    self.impropers_2[(name1, name2)] = p
		else:
		    self.impropers_1[(name1, name2, name3)] = p
	    else:
		self.impropers[(name1, name2, name3, name4)] = p

    def _setLJParameters(self, lj_dict, ljpar_set):
        for i in lj_dict.keys():
            name = _normalizeName(i)
            par1 = lj_dict[i][0]
            par2 = lj_dict[i][1]
            par3 = 0.0
            ljpar_set.addEntry(name, par1, par2, par3)
            

    def bondParameters(self, name1, name2):
	name1 = _normalizeName(name1)
	name2 = _normalizeName(name2)
	name1, name2 = _sort(name1, name2)
	p = self.bonds[(name1, name2)]
	return (p.l*amber_length_unit,
		p.k*amber_energy_unit/amber_length_unit**2)

    def bondAngleParameters(self, name1, name2, name3):
	name1 = _normalizeName(name1)
	name2 = _normalizeName(name2)
	name3 = _normalizeName(name3)
	name1, name3 = _sort(name1, name3)
	p = self.bond_angles[(name1, name2, name3)]
	return (p.angle*amber_angle_unit,
		p.k*amber_energy_unit/amber_fc_angle_unit)

    def dihedralParameters(self, name1, name2, name3, name4):
	name1 = _normalizeName(name1)
	name2 = _normalizeName(name2)
	name3 = _normalizeName(name3)
	name4 = _normalizeName(name4)
	name1, name2, name3, name4 = _sort4(name1, name2, name3, name4)
	try:
	    p = self.dihedrals[(name1, name2, name3, name4)]
	except KeyError:
	    try:
		p = self.dihedrals_2[(name2, name3)]
	    except KeyError: p = None
	if p is None:
	    return p
	else:
	    return map(lambda p: (p[2], p[1]*amber_angle_unit,
				  p[0]*amber_energy_unit),
		       p.terms)

    def improperParameters(self, name1, name2, name3, name4):
	name1 = _normalizeName(name1)
	name2 = _normalizeName(name2)
	name3 = _normalizeName(name3)
	name4 = _normalizeName(name4)
	name2, name3, name4 = _sort3(name2, name3, name4)
	p = None
	try:
	    p = self.impropers[(name1, name2, name3, name4)]
	except KeyError:
	    for names in [(name2, name3), (name2, name4), (name3, name4)]:
		try:
		    p = self.impropers_1[(name1,) + names]
		except KeyError: pass
	    if p is None:
		for name in [name2, name3, name4]:
		    try:
			p = self.impropers_2[(name1, name)]
		    except KeyError: pass
	if p is None:
	    return p
	else:
	    return map(lambda p: (p[2], p[1]*amber_angle_unit,
				  p[0]*amber_energy_unit),
		       p.terms)

    def ljParameters(self, name, parameter_set = None):
	if parameter_set is None:
	    ljp = self.default_ljpar_set
	else:
	    ljp = self.ljpar_sets[parameter_set]
	try:
	    name = self.lj_equivalent[name]
	except KeyError: pass
	if ljp.type == 'RE':
	    p = ljp.getEntry(name)
	    return (p[1]*amber_energy_unit,
		    1.78179743628*p[0]*amber_length_unit, 0)
        elif ljp.type == 'AC':
	    p = ljp.getEntry(name)
            a = p[0]*amber_energy_unit*amber_length_unit**12
            c = p[1]*amber_energy_unit*amber_length_unit**6
            if c == 0.:
                eps = 0.
            else:
                eps = 0.25*c*c/a
            if a == 0.:
                sigma = 0.
            else:
                sigma = (a/c)**(1./6.)
	    return (eps, sigma, 1)
	else:
	    raise ValueError, 'Type ' + ljp.type + ' not supported.'

#
# Atom type class
#
class AmberAtomType:

    def __init__(self, name, mass):
	self.name = name
	self.mass = mass
	self.hydrophylic = 0

#
# Bond parameter class
#
class AmberBondParameters:

    def __init__(self, a1, a2, k, l):
	self.a1 = a1
	self.a2 = a2
	self.k = k
	self.l = l

#
# Bond angle parameter class
#
class AmberBondAngleParameters:

    def __init__(self, a1, a2, a3, k, angle):
	self.a1 = a1
	self.a2 = a2
	self.a3 = a3
	self.k = k
	self.angle = angle

#
# Dihedral angle parameter class
#
class AmberDihedralParameters:

    def __init__(self, a1, a2, a3, a4, divf, k, delta, n):
	self.a1 = a1
	self.a2 = a2
	self.a3 = a3
	self.a4 = a4
	self.terms = [(k/divf, delta, int(abs(n)))]

    def addTerm(self, divf, k, delta, n):
	self.terms.append((k/divf, delta, int(abs(n))))

    def __repr__(self):
	if self.a1 is None and self.a4 is None:
	    return 'X-' + self.a2.name + '-' + self.a3.name + '-X'
	else:
	    return self.a1.name + '-' + self.a2.name + '-' + \
		   self.a3.name + '-' + self.a4.name
    __str__ = __repr__

#
# Improper dihedral angle parameter class
#
class AmberImproperParameters:

    def __init__(self, a1, a2, a3, a4, divf, k, delta, n):
	self.a1 = a1
	self.a2 = a2
	self.a3 = a3
	self.a4 = a4
	self.terms = [(0.5*k, delta, int(abs(n)))]

    def addTerm(self, divf, k, delta, n):
	self.terms.append((0.5*k, delta, int(abs(n))))

    def __repr__(self):
	if self.a4 is None:
	    if self.a3 is None:
		return self.a1.name + '-' + self.a2.name + '-X-X'
	    else:
		return self.a1.name + '-' + self.a2.name + '-' + \
		       self.a3.name + '-X'
	else:
	    return self.a1.name + '-' + self.a2.name + '-' + \
		   self.a3.name + '-' + self.a4.name
    __str__ = __repr__

#
# H-bond parameter class
#
class AmberHbondParameters:

    def __init__(self, a1, a2, a, b):
	self.a1 = a1
	self.a2 = a2
	self.a = a
	self.b = b

#
# Lennard-Jones parameter sets
#
class AmberLJParameterSet:

    def __init__(self, name, type):
	self.name = name
	self.type = type
	self.entries = {}

    def addEntry(self, name, p1, p2, p3):
	if self.type == 'SK':
	    self.entries[name] = (p1, p2, p3)
	else:
	    self.entries[name] = (p1, p2)

    def getEntry(self, name):
	return self.entries[name]

#
# Utility functions
#
def _normalizeName(name):
    return string.upper(string.strip(name))

def _sort(s1, s2):
    if s2 < s1:
	return s2, s1
    else:
	return s1, s2

def _sort3(s1, s2, s3):
    if s2 < s1:
	s1, s2 = s2, s1
    if s3 < s2:
	s2, s3 = s3, s2
    if s2 < s1:
	s1, s2 = s2, s1
    return s1, s2, s3

def _sort4(s1, s2, s3, s4):
    if s3 < s2:
	return s4, s3, s2, s1
    else:
	return s1, s2, s3, s4

--------------7229F9E3D257C0C5BE30B1E6
Content-Type: text/plain; charset=us-ascii;
 name="AmberListForceField.py"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="AmberListForceField.py"

# This file provides the Amber force field, using Amber parameter files.
# This module is modified from AmberForceField.py, so that the Amber94-
# ForceField also takes a dictionary of amber parameters as input
#
# Written by Konrad Hinsen
# last revision: 2000-9-25
#
# Modified by Berit Hinnemann
# last revision: 2001-11-06

_undocumented = 1

from MMTK.ForceFields import MMForceField
import AmberListData
import os

#
# Read parameter files
#
Amber94 = None
Amber91 = None
OPLS = None

this_directory = os.path.split(__file__)[0]

def readAmber94List(amber_params = None, main_file = None, mod_files = None):
    global Amber94
    if main_file is None and mod_files is None:
        if Amber94 is None:
            paramfile = os.path.join(this_directory, "amber_parm94")
            modfile = os.path.join(this_directory, "amber_parm94.heme")
            Amber94 = AmberListData.AmberListParameters(amber_params, paramfile, [(modfile, 'MOD4')])
            Amber94.lennard_jones_1_4 = 0.5
            Amber94.electrostatic_1_4 = 1./1.2
            Amber94.default_ljpar_set = Amber94.ljpar_sets['MOD4']
            Amber94.atom_type_property = 'amber_atom_type'
            Amber94.charge_property = 'amber_charge'
        return Amber94
    else:
        if main_file is None:
            main_file = os.path.join(this_directory, "amber_parm94")
        mod_files = map(lambda mf: (mf, 'MOD4'), mod_files)
        params = AmberListData.AmberListParameters(amber_params, main_file, mod_files)
        params.lennard_jones_1_4 = 0.5
        params.electrostatic_1_4 = 1./1.2
        params.default_ljpar_set = params.ljpar_sets['MOD4']
        params.atom_type_property = 'amber_atom_type'
        params.charge_property = 'amber_charge'
        return params

def readAmber91():
    global Amber91
    if Amber91 is None:
        Amber91 = AmberData.AmberParameters(os.path.join(this_directory,
                                                         "amber_parm91"))
        Amber91.lennard_jones_1_4 = 0.5
        Amber91.electrostatic_1_4 = 0.5
        Amber91.default_ljpar_set = Amber91.ljpar_sets['STDA']
        Amber91.atom_type_property = 'amber91_atom_type'
        Amber91.charge_property = 'amber91_charge'

def readOPLS():
    global OPLS
    if OPLS is None:
        OPLS = AmberData.AmberParameters(os.path.join(this_directory,
                                                      "opls_parm"))
        OPLS.lennard_jones_1_4 = 0.125
        OPLS.electrostatic_1_4 = 0.5
        OPLS.default_ljpar_set = OPLS.ljpar_sets['OPLS']
        OPLS.atom_type_property = 'opls_atom_type'
        OPLS.charge_property = 'opls_charge'

#
# The total force field
#
class Amber94ListForceField(MMForceField.MMForceField):

    """Amber 94 force field
    
    Change Berit: This force field is changed, so that it also
    takes a list of AMBER parameters as input, which gives
    an interface to AMBER directly in the program and not
    only through a modification file

    Constructor: Amber94ForceField(|lennard_jones_options|,
                                   |electrostatic_options|,
                                   |amber_params|)

    Arguments:

    |lennard_jones_options| -- parameters for Lennard-Jones interactions;
                               one of:

      * a number, specifying the cutoff
      * 'None', meaning the default method (no cutoff; inclusion of all
        pairs, using the minimum-image conventions for periodic universes)
      * a dictionary with an entry "method" which specifies the
        calculation method as either "direct" (all pair terms)
        or "cutoff", with the cutoff specified by the dictionary
        entry "cutoff".

    |electrostatic_options| -- parameters for electrostatic interactions;
                               one of:

      * a number, specifying the cutoff
      * 'None', meaning the default method (all pairs without cutoff for
        non-periodic system, Ewald summation for periodic systems)
      * a dictionary with an entry "method" which specifies the
        calculation method as either "direct" (all pair terms),
        "cutoff" (with the cutoff specified by the dictionary
        entry "cutoff"), "ewald" (Ewald summation, only for periodic
        universes), "screened" (see below),
        or "multipole" (fast-multipole method).

    Pair interactions in periodic systems are calculated using the
    minimum-image convention; the cutoff should therefore never be
    larger than half the smallest edge length of the elementary
    cell.

    For Lennard-Jones interactions, all terms for pairs whose distance
    exceeds the cutoff are set to zero, without any form of correction.
    For electrostatic interactions, a charge-neutralizing surface charge
    density is added around the cutoff sphere in order to reduce
    cutoff effects [Article:Wolf1999].

    For Ewald summation, there are some additional parameters that can
    be specified by dictionary entries:

    - "beta" specifies the Ewald screening parameter
    - "real_cutoff" specifies the cutoff for the real-space sum.
      It should be significantly larger than 1/beta to ensure that
      the neglected terms are small.
    - "reciprocal_cutoff" specifies the cutoff for the reciprocal-space sum.
      Note that, like the real-space cutoff, this is a distance; it describes
      the smallest wavelength of plane waves to take into account.
      Consequently, a smaller value means a more precise (and more expensive)
      calculation.

    MMTK provides default values for these parameter which are calculated
    as a function of the system size. However, these parameters are
    exaggerated in most cases of practical interest and can lead to excessive
    calculation times for large systems. It is preferable to determine
    suitable values empirically for the specific system to be simulated.

    The method "screened" uses the real-space part of the Ewald sum
    with a charge-neutralizing surface charge density around the
    cutoff sphere, and no reciprocal sum [Article:Wolf1999]. It requires
    the specification of the dictionary entries "cutoff" and "beta".

    The fast-multipole method uses the DPMTA library [Article:DPMTA].
    Note that this method provides only energy and forces, but no
    second-derivative matrix. There are several optional dictionary
    entries for this method, all of which are set to reasonable default
    values. The entries are "spatial_decomposition_levels",
    "multipole_expansion_terms", "use_fft", "fft_blocking_factor",
    "macroscopic_expansion_terms", and "multipole_acceptance".
    For an explanation of these options, refer to the
    DPMTA manual [Article:DPMTA].

    |amber_params| -- dictionary of amber parameters, which are also
                      set in addition to a modification file

    Format of the dictionary:
    
    * under the entry "atom_types" should be a dictionary containing
      the AMBER atom types as keys and the masses as values
    * under the entry "bonds" should be a dictionary with th keys
      being a string which specifies the bond and as values a list
      of the bond parameters, i.e. force constant and equilibrium
      length
    * under the entry "angles" should be a dictionary containing
      the angle specification as key and as values a list containg
      the parameters (force constant and equilibrium angle)
    * under the entry "dihedrals" should be a dictionary with keys
      being dihedral specification and values a list of the
      parameters
    * under the entry "impropers" should be a dictionary with keys
      being the specification of the impropers and as values a list
      of the parameters
    * under the entry "lennard_jones" should be a dictionary with
      the atom type as key and a list of the parameters as values
    
    """

    def __init__(self, lj_options = None, es_options = None,
                 amber_params = None,
                 bonded_scale_factor = 1., **kwargs):
	self.arguments = (lj_options, es_options, bonded_scale_factor)
        main_file = kwargs.get('parameter_file', None)
        mod_files = kwargs.get('mod_files', None)
        parameters = readAmber94List(amber_params, main_file, mod_files)
        MMForceField.MMForceField.__init__(self, 'Amber94', parameters,
                                           lj_options, es_options,
                                           bonded_scale_factor)

AmberForceField = Amber94ListForceField

class Amber91ForceField(MMForceField.MMForceField):

    def __init__(self, lj_options = None, es_options = None,
                 bonded_scale_factor=1., **kwargs):
	self.arguments = (lj_options, es_options, bonded_scale_factor)
        readAmber91()
        MMForceField.MMForceField.__init__(self, 'Amber91', Amber91,
                                           lj_options, es_options,
                                           bonded_scale_factor)

class OPLSForceField(MMForceField.MMForceField):

    def __init__(self, lj_options = None, es_options = None,
                 bonded_scale_factor=1., **kwargs):
	self.arguments = (lj_options, es_options, bonded_scale_factor)
        readOPLS()
        MMForceField.MMForceField.__init__(self, 'OPLS', OPLS,
                                           lj_options, es_options,
                                           bonded_scale_factor)


#
# The following classes provides access to individual terms of the
# Amber force field. They are not used normally, but can be useful
# for testing.
#

#
# Bonded interactions
#
class AmberBondedForceField(MMForceField.MMBondedForceField):

    def __init__(self):
        readAmber94()
	MMForceField.MMBondedForceField.__init__(self, 'Amber_bonded', Amber94)

#
# Nonbonded interactions
#
class AmberLJForceField(MMForceField.MMLJForceField):

    def __init__(self, cutoff = None):
        readAmber94()
	MMForceField.MMLJForceField.__init__(self, 'Amber_LJ', Amber94, cutoff)

class AmberESForceField(MMForceField.MMESForceField):

    def __init__(self, cutoff = None):
        readAmber94()
	MMForceField.MMESForceField.__init__(self, 'Amber_ES', Amber94, cutoff)

class AmberEwaldESForceField(MMForceField.MMEwaldESForceField):

    def __init__(self, options = {}):
        readAmber94()
	MMForceField.MMEwaldForceESField.__init__(self, 'Amber_ES',
                                                  Amber94, options)

class AmberMPESForceField(MMForceField.MMMPESForceField):

    def __init__(self, options = {}):
        readAmber94()
	MMForceField.MMMPESForceField.__init__(self, 'Amber_ES',
                                               Amber94, options)

class AmberNonbondedForceField(MMForceField.MMNonbondedForceField):

    def __init__(self, lj_options = None, es_options = None):
        readAmber94()
        MMForceField.MMNonbondedForceField.__init__(self, Amber94,
                                                    lj_options, es_options)

--------------7229F9E3D257C0C5BE30B1E6--