[MMTK] Strategies for GAFF support into MMTK-2.5.18

[Mark Williamson]

Dear All,

I've been trying to "incorporate" GAFF 
(http://dx.doi.org/10.1002/jcc.20035) support into the MMTK-2.5.18

I'm using the following additions to 
MMTK/ForceFields/Amber/AmberForceField.py :

def readGAFF(main_file = None, mod_files = None):
      global GAFF
      if main_file is None and mod_files is None:
          if GAFF is None:
              paramfile = os.path.join(this_directory, "gaff.dat")
              GAFF = AmberData.AmberParameters(paramfile, [(modfile, 
'MOD4')])
              GAFF.lennard_jones_1_4 = 0.5
              GAFF.electrostatic_1_4 = 1./1.2
              GAFF.default_ljpar_set = GAFF.ljpar_sets['MOD4']
              GAFF.atom_type_property = 'GAFF_atom_type'
              GAFF.charge_property = 'GAFF_charge'
          return GAFF
      else:
          if main_file is None:
              main_file = os.path.join(this_directory, "gaff.dat")
          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 = 'GAFF_atom_type'
          params.charge_property = 'GAFF_charge'
          return params


class GAFFForceField(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 = readGAFF(main_file, mod_files)
          MMForceField.MMForceField.__init__(self, 'GAFF', parameters,
                                             lj_options, es_options,
                                             bonded_scale_factor)


I've also deposited the gaff.dat forcefile file in that directory. Now, 
I've hit an interesting problem due to GAFF's dependence of lower case 
for atom types. As a matter of internal consistency, _normalizeName 
coverts all the atoms types to upper case when AmberData.py parses in 
the relevant data file for that particular AMBER force field. One simple 
hack is to disable this case conversion:

in MMTK/ForceFields/Amber/AmberData.py

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


But, this is hacky and would break other things. One side effect of this 
would be that GAFF atom types would not be able to co-exist with AMBER 
FF94 types in the same simulation.

So, I'm after some advice on how to do this "the right way". One idea I 
had was to prepend every GAFF atom type with a unique string, hence the 
GAFF atom type "ca" would become, internally to MMTK, "GAFF_CA"

"ca-ca" would become "GAFF_CA-GAFF_CA"

I know this is a bit sketchy at the moment, but is this a good way forward?

regards,

Mark