[MMTK] Adding/modifying forcefields

[konrad.hinsen at laposte.net]

On Dec 17, 2004, at 16:49, david.evans at ulsop.ac.uk wrote:

> Thanks for the reply. I was hoping that I could use the existing
> Amber94 classes, and just change the input parameter files, i.e. use
> something instead of 'amber_parm94'. The functional form of Amber99
> and gaff is identical to Amber94, I believe, so I thought it  would
> just be a question of ensuring the correct data was read in.

It is, but it is a bit more work than you suspect, which is why I 
didn't do it already.

The Python code you need is minor, and just a clone of the Amber 94 
code. I have appended it below (save as Amber99ForceField.py). The 
parameter file needs to be called amber_parm99 and must be located in 
the same directory as the module, but both details are easy to change. 
Something to verify is the 1-4 factors, if they are different from 
their 94 value, you must change the respective lines.

However, this is just the first step of using Amber 99. The second one 
is simple in principle but quite time-consuming: you must put the atom 
type and charge parameters for Amber 99 into all the group and molecule 
description files that you plan to use. If you work with proteins, 
that's a lot! You need to define amber99_atom_type and amber99_charge 
in the same way as amber_atom_type and amber_charge are defined for the 
94 version. You can either add the data to the existing database files, 
or make a copy of the files you plan to modify and put them into 
~/.mmtk/Database. You can find all the information in the Amber 
topology files.

The original MMTK database files were machine-generated from the Amber 
94 topology files to reduce the effort (and mistakes). In principle, 
that could be done again, but if I remember correctly Amber has changed 
the format of the topology files, so a new converter would have to be 
written. Moreover, a lot of manual fine-tuning was applied to the MMTK 
files in order to be able to read PDB files of various origins, all 
that would be lost with a simple regeneration.

Still, I think an automatic approach is the only viable solution unless 
you plan to use just a few molecules.

Konrad.
--
---------------------------------------------------------------------
Konrad Hinsen
Laboratoire Léon Brillouin, CEA Saclay,
91191 Gif-sur-Yvette Cedex, France
Tel.: +33-1 69 08 79 25
Fax: +33-1 69 08 82 61
E-Mail: hinsen at llb.saclay.cea.fr
---------------------------------------------------------------------

# Amber 99 force field
#
# Written by Konrad Hinsen
# last revision: 2004-12-17
#

from MMTK.ForceFields import MMForceField
from MMTK.ForceFields.Amber import AmberData
import os

#
# Read parameter files
#
Amber99 = None

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

def readAmber99(main_file = None, mod_files = None):
     global Amber99
     if main_file is None and mod_files is None:
         if Amber99 is None:
             paramfile = os.path.join(this_directory, "amber_parm99")
             Amber99 = AmberData.AmberParameters(paramfile)
             Amber99.lennard_jones_1_4 = 0.5
             Amber99.electrostatic_1_4 = 1./1.2
             Amber99.default_ljpar_set = Amber99.ljpar_sets['MOD4']
             Amber99.atom_type_property = 'amber99_atom_type'
             Amber99.charge_property = 'amber99_charge'
         return Amber99
     else:
         if main_file is None:
             main_file = os.path.join(this_directory, "amber_parm99")
         mod_files = map(lambda mf: (mf, 'MOD4'), mod_files)
         params = AmberData.AmberParameters(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 = 'amber99_atom_type'
         params.charge_property = 'amber99_charge'
         return params

#
# The force field class
#
class Amber99ForceField(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)
         main_file = kwargs.get('parameter_file', None)
         mod_files = kwargs.get('mod_files', None)
         parameters = readAmber99(main_file, mod_files)
         MMForceField.MMForceField.__init__(self, 'Amber99', parameters,
                                            lj_options, es_options,
                                            bonded_scale_factor)