[MMTK] GAFF parameters in MMTK

[robjorissen at yahoo.com.au]

G'day

A few people have made contributions to this list discussing the use of the GAFF forcefield. So, I wrote some code which might help. This Python script uses as input your small molecule of interest in MDL/SDF format and then uses the Antechamber suite of programs (from UCSF) to generate appropriate forcefield parameters and then store the resulting forcefield file (and a couple of other files) in a Molecules directory which is generated if it does not already exist.

A few more notes:
Obviously you need to get the Antechamber suite of programs (http://ambermd.org/antechamber/antechamber.html). From memory, if you compile in the PC/Cygwin environment, you should get the package bison which seems to used instead of YACC. It has been a while since I have compiled Antechamber.

And obviously you probably need to change the value of the ACHOME environmental variable (set in this script) to the location you have compiled Antechamber.

I hope people find this useful. I can't guarantee it is bug-free, but it seems to work for me.

Rob Jorissen
Ludwig Insitute for Cancer Research



#! /usr/bin/env python


import os, shutil, string, subprocess, sys

# add ACHOME environmental variable
ACHOME = "C:/Work_programs/antechamber-1.27"
os.environ["ACHOME"] = ACHOME

inFormat = "mdl"
outSuffix = ".prepin"
# atom_type = "amber"        # gaff, amber, bcc and sybyl
atom_type = "gaff"        # gaff, amber, bcc and sybyl

remove_files = "y"        # removes intermediate files
charge_method = 2        # corresponds to AM1-BCC

remove_prepin = "n"


def writeWithBreaks(out, list, startText, endText, perLine):
  out.write(startText)
  s = " " * (len(startText) - 1)
  n1 = len(list) - 1
  for i, l in enumerate(list):
    out.write(l)  # provisional
    if i == n1:
      pass
    elif perLine == 1:
      out.write(",\n%s" % s)
    elif (i+1) % perLine == 0 and i != 0:
      out..write(",\n%s" % s)
    else:
      out.write(", ")
  out..write(endText + "\n")



if len(sys.argv) > 2:
  netCharge = int(sys.argv[2])
else:
  netCharge = 0

# consider adding code to automatically determine the net formal
# charge from the input file


# 1. Run antechamber
exe = ACHOME + "/exe/antechamber"
prefix = (sys.argv[1].split("."))[0]
outName = prefix + "." + atom_type + ".prepin"
cmdList = [exe, "-i", sys.argv[1], "-fi", inFormat, "-o", outName, "-fo ac -c",
           str(charge_method), "-s 0 -pf", remove_files, "-at", atom_type,
           "-nc", str(netCharge)]
cmd = string.join(cmdList, " ")

subprocess.call(cmd, shell=True)
os.remove("mopac.in")
os.remove("mopac.out")


# 2. Run parmchk to find forcefield parameters not in the gaff.dat file

exe = ACHOME + "/exe/parmchk"
inParmchk = outName
outParmchk = prefix + "." + atom_type + ".parmout"

cmdList = [exe, "-i", inParmchk, "-o", outParmchk, "-f ac"]
cmd = string.join(cmdList, " ")

subprocess..call(cmd, shell=True)

# may need to add gaff.dat to this command...


# 3. Read in and process the prepin file
prepinList = open(outName, "r").readlines()
## os.remove(outName)

name = prefix
PDBname = "lig"

i0 = 2
i1 = -1
atomList = []
bondList = []
for i, line in enumerate(prepinList[2: ]):
  if line[0:4] == "ATOM":
    i1 = i
    atomname = line[12: 16].strip()
    ligname = line[17: 20].strip()
    xcoord= float(line[30: 38]..strip())
    ycoord= float(line[38: 46].strip())
    zcoord= float(line[46: 54].strip())
    pcharge = float(line[54: 64].strip())
    type = line[69: 74].strip()

    i = 0
    while atomname[i].isalpha():
      i += 1
    element = atomname[: i]

    atomList.append([element, atomname, pcharge, type,
                     xcoord, ycoord, zcoord])

  elif line[0:4] == "BOND":
    atom1 = line[8: 14].strip()
    atom2 = line[14: 19].strip()
    bType = line[19: 24].strip()
    aName1 = line[27: 32].strip()
    aName2 = line[32: 37].strip()
# still not working for all cases...   
    # print "_%s_%s...%s" % (aName1, aName2, line)
    bondList.append([aName1, aName2])


# 4. write out a parameter file to the Molecules subdirectory

thisDir = os.getcwd()
MoleculesDir = os.path.join(thisDir, "Molecules")
if not os.path.isdir(MoleculesDir):
  os.mkdir(MoleculesDir)


outFile = open(os.path.join(MoleculesDir, prefix), "w")

# title
outFile.write("name = \"%s\"\n\n" % prefix)

# atom names and elements
for atom in atomList:
  outFile.write("%s = Atom(\"%s\")\n" % (atom[1], atom[0]))
  # print atom




bondListWrite = []
for bond in bondList:
  bondListWrite.append("Bond(%s, %s)" %(bond[0], bond[1]))

writeWithBreaks(outFile, bondListWrite, "\nbonds = [", "]", 4)

typeListWrite = []
pdbListWrite = []
coordListWrite = []
for atom in atomList:
  typeListWrite.append("%s: \"%s\"" %(atom[1], atom[3]))
  pdbListWrite.append("\"%s\": %s" % (atom[1], atom[1]))
  coordListWrite.append("%s: (%.3f,  %.3f, %.3f)" %
                         (atom[1], atom[4], atom[5], atom[6]))

# write out the atom types
if atom_type == "gaff":
  atom_type = "GAFF"
thisStr = "\n" + atom_type + "_atom_type = {"
# thisStr = "\n" + "amber" + "_atom_type = {"
writeWithBreaks(outFile, typeListWrite, thisStr, "}", 4)


# write out the pdb names
startTxt = "\npdbmap = [(\"%s\", {" % ligname
writeWithBreaks(outFile, pdbListWrite, startTxt, "})]", 4)


# write out the partial charges
outFile.write("\n%s_charge = {\n" % atom_type)
for atom in atomList:
  outFile.write("    %-5s : %7.4f,\n" % (atom[1], atom[2]))
outFile.write("    }\n\n")


# write out a default set of coordinates to the parameter file
outFile.write("configurations = {\n")
writeWithBreaks(outFile, coordListWrite,
                "    \"default\": Cartesian({","})", 1)
outFile.write("    }\n")

outFile.close()



# move the parmout file into the Molecules subdirectory
parmName = "%s.%s.parmout" % (prefix, atom_type)
shutil.move(parmName, MoleculesDir)


# rename the prepin file, potentially overwriting previous version
newPDBfilename = os.path.join(MoleculesDir, (prefix + ".pdb"))
try:
  os.remove(newPDBfilename)
except:
  pass

os.rename(outName, newPDBfilename)







      Make Yahoo!7 your homepage and win a trip to the Quiksilver Pro. Find out more