mirror of
https://github.com/openmm/openmm
synced 2026-06-03 06:39:48 +09:00
1253 lines
55 KiB
Python
1253 lines
55 KiB
Python
#=============================================================================================
|
|
# MODULE DOCSTRING
|
|
#=============================================================================================
|
|
|
|
"""
|
|
processTinkerForceField.py
|
|
|
|
Convert TINKER force field files into xml files for use by pyopenmm
|
|
|
|
(1) read residue template file
|
|
(2) read TINKER parameter file
|
|
(3) assign biotypes to each atom in residue template file
|
|
(4) output force-field parameter file
|
|
|
|
"""
|
|
#=============================================================================================
|
|
# GLOBAL IMPORTS
|
|
#=============================================================================================
|
|
|
|
import os
|
|
import xml.etree.ElementTree as etree
|
|
import sys
|
|
import shlex
|
|
import math
|
|
import datetime
|
|
import os.path
|
|
|
|
#=============================================================================================
|
|
# Ion list
|
|
#=============================================================================================
|
|
|
|
# biotype 2003 NA "Sodium Ion" 250
|
|
# biotype 2004 K "Potassium Ion" 251
|
|
# biotype 2005 MG "Magnesium Ion" 255
|
|
# biotype 2006 CA "Calcium Ion" 256
|
|
# biotype 2007 CL "Chloride Ion" 258
|
|
|
|
# Ions for amoeabio18.prm
|
|
ions2018 = { 'Li+' : ['LI', 351],
|
|
'Na+' : ['NA', 352],
|
|
'K+' : ['K', 353],
|
|
'Rb+' : ['RB', 354],
|
|
'Cs+' : ['CS', 355],
|
|
'Be2' : ['BE', 356],
|
|
'Mg2' : ['MG', 357],
|
|
'Ca2' : ['CA', 358],
|
|
'Sr2' : ['SR', 359],
|
|
'Ba2' : ['BA', 360],
|
|
'Zn2' : ['ZN', 361],
|
|
'F-' : ['F', 362],
|
|
'Cl-' : ['Cl', 363],
|
|
'Br-' : ['Br', 364],
|
|
'I-' : ['I', 365]
|
|
}
|
|
|
|
# Ions for amoebapro13.prm
|
|
ions2013 = { 'Li+' : ['LI', 249],
|
|
'Na+' : ['NA', 250],
|
|
'K+' : ['K', 251],
|
|
'Rb+' : ['RB', 252],
|
|
'Cs+' : ['CS', 253],
|
|
'Be2' : ['BE', 254],
|
|
'Mg2' : ['MG', 255],
|
|
'Ca2' : ['CA', 256],
|
|
'Sr2' : ['SR', 257],
|
|
'Ba2' : ['BA', 258],
|
|
'Zn2' : ['ZN', 259],
|
|
'F-' : ['F', 260],
|
|
'Cl-' : ['Cl', 261],
|
|
'Br-' : ['Br', 262],
|
|
'I-' : ['I', 263]
|
|
}
|
|
|
|
# Ions for amoebabio09.prm
|
|
ions2009 = { 'Li+' : ['LI', 404],
|
|
'Na+' : ['NA', 405],
|
|
'K+' : ['K', 406],
|
|
'Rb+' : ['RB', 407],
|
|
'Cs+' : ['CS', 408],
|
|
'Be2' : ['BE', 409],
|
|
'Mg2' : ['MG', 410],
|
|
'Ca2' : ['CA', 411],
|
|
'Zn2' : ['ZN', 412],
|
|
'F-' : ['F', 413],
|
|
'Cl-' : ['Cl', 414],
|
|
'Br-' : ['Br', 415],
|
|
'I-' : ['I', 416]
|
|
}
|
|
|
|
atomTypes = {}
|
|
bioTypes = {}
|
|
|
|
#=============================================================================================
|
|
# Helper function for XML formatting
|
|
#=============================================================================================
|
|
|
|
def xml_float(x, min_decimals=1):
|
|
"""Format a float for XML, cleaning tiny floating-point errors and keeping at least one decimal."""
|
|
f = float(x)
|
|
s = format(f, ".15g")
|
|
# Only add decimal if it is not already scientific and has no decimal
|
|
if "e" not in s and "." not in s:
|
|
s += "." + "0" * min_decimals
|
|
return s
|
|
|
|
#=============================================================================================
|
|
# Default 'constructor' for atoms
|
|
#=============================================================================================
|
|
|
|
def getDefaultAtom( ):
|
|
|
|
atom = dict();
|
|
atom['tinkerLookupName'] = 'XXX'
|
|
atom['type'] = -1
|
|
atom['bonds'] = dict()
|
|
|
|
return atom
|
|
|
|
#=============================================================================================
|
|
# Add bond to atomDict[]; atoms are added to atomDict[] if missing
|
|
#=============================================================================================
|
|
|
|
def addBond( atomDict, atom1, atom2 ):
|
|
if( atom1 not in atomDict ):
|
|
atomDict[atom1] = getDefaultAtom()
|
|
|
|
if( atom2 not in atomDict ):
|
|
atomDict[atom2] = getDefaultAtom()
|
|
|
|
atomDict[atom2]['bonds'][atom1] = 1
|
|
atomDict[atom1]['bonds'][atom2] = 1
|
|
|
|
#=============================================================================================
|
|
# Get atom dictionary from xml atom list
|
|
#=============================================================================================
|
|
|
|
def getXmlAtoms( atoms ):
|
|
|
|
atomInfo = dict();
|
|
for atom in atoms:
|
|
name = atom.attrib['name']
|
|
atomInfo[name] = getDefaultAtom()
|
|
atomInfo[name]['tinkerLookupName'] = atom.attrib['tinkerLookupName']
|
|
|
|
return atomInfo
|
|
|
|
#=============================================================================================
|
|
# Get bond dictionary from xml bond list
|
|
#=============================================================================================
|
|
|
|
def getXmlBonds( bonds ):
|
|
|
|
bondInfo = dict();
|
|
for bond in bonds:
|
|
atom1 = bond.attrib['from']
|
|
atom2 = bond.attrib['to']
|
|
if( atom1 not in bondInfo ):
|
|
bondInfo[atom1] = dict()
|
|
if( atom2 not in bondInfo ):
|
|
bondInfo[atom2] = dict()
|
|
|
|
bondInfo[atom1][atom2] = 1
|
|
bondInfo[atom2][atom1] = 1
|
|
|
|
return bondInfo
|
|
|
|
#=============================================================================================
|
|
# Build entry for protein residue
|
|
#=============================================================================================
|
|
|
|
def buildProteinResidue( residueDict, atoms, bondInfo, abbr, loc, tinkerLookupName, include, residueName, type ):
|
|
|
|
# residueDict[abbr] abbr=ALA, CALA, NALA, ...
|
|
# residueDict[abbr]['atoms'] list if atom dict()
|
|
# residueDict[abbr]['type'] molecule type ('protein', 'nucleic acid', ...)
|
|
# residueDict[abbr]['tinkerLookupName'] Tinker lookup name
|
|
# residueDict[abbr]['residueName'] residueName
|
|
# residueDict[abbr]['include'] include in output
|
|
|
|
residueDict[abbr] = dict()
|
|
residueDict[abbr]['atoms'] = atoms
|
|
residueDict[abbr]['type'] = type
|
|
residueDict[abbr]['loc'] = loc
|
|
residueDict[abbr]['tinkerLookupName'] = tinkerLookupName
|
|
residueDict[abbr]['residueName'] = residueName
|
|
residueDict[abbr]['include'] = include
|
|
|
|
# for each bond, add entry to
|
|
# residueDict[abbr]['atoms'][atom]['bonds']
|
|
# residueDict[abbr]['atoms'][bondedAtom]['bonds']
|
|
|
|
for atom in bondInfo:
|
|
if( atom in residueDict[abbr]['atoms'] ):
|
|
|
|
if( 'bonds' not in residueDict[abbr]['atoms'][atom] ):
|
|
residueDict[abbr]['atoms'][atom]['bonds'] = dict()
|
|
|
|
for bondedAtom in bondInfo[atom]:
|
|
if( bondedAtom in residueDict[abbr]['atoms'] ):
|
|
if( 'bonds' not in residueDict[abbr]['atoms'][bondedAtom] ):
|
|
residueDict[abbr]['atoms'][bondedAtom]['bonds'] = dict()
|
|
residueDict[abbr]['atoms'][bondedAtom]['bonds'][atom] = 1
|
|
residueDict[abbr]['atoms'][atom]['bonds'][bondedAtom] = 1
|
|
else:
|
|
print("Error: bonded atom=%s not in residue=%s" % ( atom, abbr ))
|
|
else:
|
|
print("Error: bonded atom=%s nt in residue=%s" % ( atom, abbr ))
|
|
|
|
return
|
|
|
|
#=============================================================================================
|
|
# Copy a bond (dict() copy)
|
|
#=============================================================================================
|
|
|
|
def copyBonds( bonds ):
|
|
bondCopy = dict()
|
|
for key in bonds.keys():
|
|
bondCopy[key] = bonds[key]
|
|
return bondCopy
|
|
|
|
#=============================================================================================
|
|
# Copy a atom (dict() copy, including the 'bonds' list)
|
|
#=============================================================================================
|
|
|
|
def copyAtom( atom ):
|
|
atomCopy = dict()
|
|
for key in atom.keys():
|
|
if( key != 'bonds' ):
|
|
atomCopy[key] = atom[key]
|
|
else:
|
|
atomCopy['bonds'] = copyBonds( atom[key] )
|
|
return atomCopy
|
|
|
|
#=============================================================================================
|
|
# Copy a residue, including atom list
|
|
#=============================================================================================
|
|
|
|
def copyProteinResidue( residue ):
|
|
|
|
residueCopy = dict()
|
|
residueCopy['atoms'] = dict()
|
|
residueCopy['type'] = residue['type']
|
|
residueCopy['loc'] = residue['loc']
|
|
residueCopy['tinkerLookupName'] = residue['tinkerLookupName']
|
|
residueCopy['residueName'] = residue['residueName']
|
|
residueCopy['include'] = residue['include']
|
|
|
|
for atom in residue['atoms']:
|
|
residueCopy['atoms'][atom] = copyAtom( residue['atoms'][atom] )
|
|
|
|
return residueCopy
|
|
|
|
#=============================================================================================
|
|
# Build residue hash based on xml file
|
|
#=============================================================================================
|
|
|
|
def buildResidueDict( residueXmlFileName ):
|
|
|
|
residueTree = etree.parse(residueXmlFileName)
|
|
print("Read %s" % (residueXmlFileName))
|
|
root = residueTree.getroot()
|
|
residueDict = dict()
|
|
|
|
# residueDict[residueName] = dict()
|
|
# ['loc']
|
|
# ['type']
|
|
# ['atoms'] = dict()
|
|
# [atomName] = dict()
|
|
# ['bonds'] = dict{}
|
|
for residue in root.findall('Residue'):
|
|
|
|
# <Residue abbreviation="MET" loc="middle" type="protein" tinkerLookupName="Methionine" fullName="Methionine">
|
|
abbr = residue.attrib['abbreviation']
|
|
loc = residue.attrib['loc']
|
|
type = residue.attrib['type']
|
|
tinkerName = residue.attrib['tinkerLookupName']
|
|
residueName = residue.attrib['fullName']
|
|
isProtein = False
|
|
isWater = False
|
|
isDNA = False
|
|
isRNA = False
|
|
if type == 'protein':
|
|
isProtein = True
|
|
elif type == 'AmoebaWater':
|
|
isWater = True
|
|
elif type == 'dna':
|
|
isDNA = True
|
|
elif type == 'rna':
|
|
isRNA = True
|
|
else:
|
|
continue
|
|
|
|
atoms = getXmlAtoms( residue.findall('Atom') )
|
|
bondInfo = getXmlBonds( residue.findall('Bond') )
|
|
|
|
# if residue is an amino acid, then create CALA and NALA residues, in addition to non-termianal residue, and include appropriate atoms
|
|
# HXT is excluded from all residues
|
|
|
|
if( isWater ):
|
|
|
|
buildProteinResidue( residueDict, atoms, bondInfo, abbr, 'x', tinkerName, 1, 'HOH', 'water' )
|
|
|
|
elif( isProtein ):
|
|
|
|
buildProteinResidue( residueDict, atoms, bondInfo, abbr, 'm', tinkerName, 1, residueName, 'protein' )
|
|
|
|
cResidueName = 'C' + abbr
|
|
residueDict[cResidueName] = copyProteinResidue( residueDict[abbr] )
|
|
residueDict[cResidueName]['loc'] = 'c'
|
|
if( residueDict[abbr]['tinkerLookupName'].find('(') > -1 ):
|
|
begin = residueDict[abbr]['tinkerLookupName'].find('(')
|
|
end = residueDict[abbr]['tinkerLookupName'].find(')') + 1
|
|
sub = residueDict[abbr]['tinkerLookupName'][begin:end]
|
|
if( sub == '(HD)' or sub == '(HE)' ):
|
|
residueDict[cResidueName]['tinkerLookupName'] = 'C-Terminal ' + 'HIS ' + sub
|
|
else:
|
|
residueDict[cResidueName]['tinkerLookupName'] = 'C-Terminal ' + abbr + ' ' + sub
|
|
print("tinkerLookupName %s %s" % ( abbr, residueDict[cResidueName]['tinkerLookupName']))
|
|
else:
|
|
residueDict[cResidueName]['tinkerLookupName'] = 'C-Terminal ' + abbr
|
|
residueDict[cResidueName]['atoms']['OXT'] = copyAtom( residueDict[abbr]['atoms']['O'] )
|
|
residueDict[cResidueName]['atoms']['OXT']['tinkerLookupName'] = 'OXT'
|
|
residueDict[cResidueName]['atoms']['O']['tinkerLookupName'] = 'OXT'
|
|
residueDict[cResidueName]['parent'] = residueDict[abbr]
|
|
|
|
nResidueName = 'N' + abbr
|
|
residueDict[nResidueName] = copyProteinResidue( residueDict[abbr] )
|
|
residueDict[nResidueName]['loc'] = 'n'
|
|
residueDict[nResidueName]['tinkerLookupName'] = 'N-Terminal ' + abbr
|
|
residueDict[nResidueName]['parent'] = residueDict[abbr]
|
|
|
|
if( abbr == 'PRO' ):
|
|
#<Atom name="H" tinkerLookupName="HN" bonds="1" />
|
|
|
|
residueDict[nResidueName]['atoms']['H2'] = getDefaultAtom()
|
|
residueDict[nResidueName]['atoms']['H3'] = getDefaultAtom()
|
|
|
|
residueDict[nResidueName]['atoms']['H2']['tinkerLookupName'] = 'HN'
|
|
residueDict[nResidueName]['atoms']['H3']['tinkerLookupName'] = 'HN'
|
|
|
|
addBond( residueDict[nResidueName]['atoms'], 'H2', 'N' )
|
|
addBond( residueDict[nResidueName]['atoms'], 'H3', 'N' )
|
|
|
|
else:
|
|
residueDict[nResidueName]['atoms']['H2'] = copyAtom( residueDict[abbr]['atoms']['H'] )
|
|
residueDict[nResidueName]['atoms']['H3'] = copyAtom( residueDict[abbr]['atoms']['H'] )
|
|
|
|
elif isDNA or isRNA:
|
|
|
|
buildProteinResidue( residueDict, atoms, bondInfo, abbr, loc, tinkerName, 1, residueName, type )
|
|
|
|
print("Start Lookup XML FFFFinal\n\n")
|
|
printXml = 1
|
|
if( printXml ):
|
|
print("<Residues>")
|
|
for resName in sorted( residueDict.keys() ):
|
|
if( 'include' in residueDict[resName] and residueDict[resName]['include'] ):
|
|
type = residueDict[resName]['type']
|
|
loc = residueDict[resName]['loc']
|
|
tinkerLookupName = residueDict[resName]['tinkerLookupName']
|
|
fullName = residueDict[resName]['residueName']
|
|
outputString = """ <Residue abbreviation="%s" loc="%s" type="%s" tinkerLookupName="%s" fullName="%s">""" % (resName, loc, type, tinkerLookupName, fullName )
|
|
print("%s" % outputString)
|
|
|
|
atomsInfo = residueDict[resName]['atoms']
|
|
for atomName in sorted( atomsInfo.keys() ):
|
|
tinkerLookupName = atomsInfo[atomName]['tinkerLookupName']
|
|
outputString = """ <Atom name="%s" tinkerLookupName="%s" />""" % (atomName, tinkerLookupName)
|
|
print("%s" % outputString)
|
|
|
|
includedBonds = dict()
|
|
for atomName in sorted( atomsInfo.keys() ):
|
|
bondsInfo = atomsInfo[atomName]['bonds']
|
|
for bondedAtom in bondsInfo:
|
|
if( bondedAtom not in includedBonds or atomName not in includedBonds[bondedAtom] ):
|
|
outputString = """ <Bond from="%s" to="%s" />""" % (atomName, bondedAtom)
|
|
if( atomName not in includedBonds ):
|
|
includedBonds[atomName] = dict()
|
|
if( bondedAtom not in includedBonds ):
|
|
includedBonds[bondedAtom] = dict()
|
|
includedBonds[atomName][bondedAtom] = 1
|
|
includedBonds[bondedAtom][atomName] = 1
|
|
print("%s" % outputString)
|
|
print("</Residue>")
|
|
print("</Residues>")
|
|
|
|
return residueDict
|
|
|
|
#=============================================================================================
|
|
# Set biotype for each atom in residueDict
|
|
#=============================================================================================
|
|
|
|
def setBioTypes( bioTypes, residueDict ):
|
|
|
|
for resname, res in residueDict.items():
|
|
for atom in res['atoms']:
|
|
atomLookup = res['atoms'][atom]['tinkerLookupName']
|
|
resLookup = []
|
|
if res['type'] == 'dna':
|
|
if res['loc'] in ('5', 'N'):
|
|
resLookup.append("5'-Hydroxyl DNA")
|
|
if res['loc'] in ('3', 'N'):
|
|
resLookup.append("3'-Hydroxyl DNA")
|
|
resLookup.append("Phosphodiester DNA")
|
|
if res['type'] == 'rna':
|
|
if res['loc'] in ('5', 'N'):
|
|
resLookup.append("5'-Hydroxyl RNA")
|
|
if res['loc'] in ('3', 'N'):
|
|
resLookup.append("3'-Hydroxyl RNA")
|
|
resLookup.append("Phosphodiester RNA")
|
|
resLookup.append(res['tinkerLookupName'])
|
|
for suffix in resLookup:
|
|
lookupName = atomLookup+'_'+suffix
|
|
if lookupName in bioTypes:
|
|
break
|
|
if lookupName in bioTypes:
|
|
res['atoms'][atom]['type'] = bioTypes[lookupName][3]
|
|
else:
|
|
print("For %s lookupName=%s not in biotype" % (atom,lookupName))
|
|
if( 'parent' in res ):
|
|
lookupName = res['atoms'][atom]['tinkerLookupName'] + '_' + res['parent']['tinkerLookupName']
|
|
if( lookupName in bioTypes ):
|
|
res['atoms'][atom]['type'] = bioTypes[lookupName][3]
|
|
else:
|
|
print("Missing lookupName=%s from biotype" % (lookupName))
|
|
return 0
|
|
|
|
#=============================================================================================
|
|
# Add multipole for forces[]; added entry is a list of axis info [kz, kx, ky] and another
|
|
# list of multipoles [charge, dipole, quadrupole]
|
|
#=============================================================================================
|
|
|
|
def addMultipole( lineIndex, allLines, forces ):
|
|
|
|
if( 'multipole' not in forces ):
|
|
forces['multipole'] = []
|
|
|
|
# axis indices and charge
|
|
|
|
fields = allLines[lineIndex]
|
|
multipoles = [ fields[-1] ]
|
|
axisInfo = fields[1:-1]
|
|
|
|
# dipole
|
|
|
|
lineIndex += 1
|
|
fields = allLines[lineIndex]
|
|
multipoles.append( fields[0] )
|
|
multipoles.append( fields[1] )
|
|
multipoles.append( fields[2] )
|
|
|
|
# quadrupole
|
|
|
|
lineIndex += 1
|
|
fields = allLines[lineIndex]
|
|
multipoles.append( fields[0] )
|
|
|
|
lineIndex += 1
|
|
fields = allLines[lineIndex]
|
|
multipoles.append( fields[0] )
|
|
multipoles.append( fields[1] )
|
|
|
|
lineIndex += 1
|
|
fields = allLines[lineIndex]
|
|
multipoles.append( fields[0] )
|
|
multipoles.append( fields[1] )
|
|
multipoles.append( fields[2] )
|
|
|
|
lineIndex += 1
|
|
|
|
# save info
|
|
|
|
multipoleInfo = [ axisInfo, multipoles ]
|
|
forces['multipole'].append( multipoleInfo )
|
|
|
|
return (lineIndex)
|
|
|
|
#=============================================================================================
|
|
# Add tortor parameters/grid to forces[]; format of each entry is [ first tortor line, grid ]
|
|
#=============================================================================================
|
|
|
|
def addTorTor( lineIndex, allLines, forces ):
|
|
|
|
if( 'tortors' not in forces ):
|
|
forces['tortors'] = []
|
|
|
|
fields = allLines[lineIndex]
|
|
tortorInfo = fields[1:]
|
|
|
|
# read grid lines
|
|
# New format has multiple data points per line: angle1 angle2 f angle1 angle2 f ...
|
|
# Need to calculate number of lines based on total grid points
|
|
|
|
nx = int(fields[6])
|
|
ny = int(fields[7])
|
|
totalGridPoints = nx * ny
|
|
|
|
grid = []
|
|
currentLineIndex = lineIndex + 1
|
|
gridPointsProcessed = 0
|
|
|
|
while gridPointsProcessed < totalGridPoints and currentLineIndex < len(allLines):
|
|
lineFields = allLines[currentLineIndex]
|
|
|
|
# Process triplets of (angle1, angle2, f) from current line
|
|
i = 0
|
|
while i + 3 <= len(lineFields) and gridPointsProcessed < totalGridPoints:
|
|
angle1 = lineFields[i]
|
|
angle2 = lineFields[i + 1]
|
|
f = lineFields[i + 2]
|
|
grid.append([angle1, angle2, f])
|
|
gridPointsProcessed += 1
|
|
i += 3
|
|
|
|
currentLineIndex += 1
|
|
|
|
forces['tortors'].append( [ tortorInfo, grid ] )
|
|
|
|
return (currentLineIndex - 1)
|
|
|
|
|
|
#=============================================================================================
|
|
|
|
# recognizedForces[] contain raw list entries from TINKER parameter file
|
|
|
|
resAtomTypes = {}
|
|
forces = {}
|
|
recognizedForces = {}
|
|
recognizedForces['bond'] = 1
|
|
recognizedForces['angle'] = 1
|
|
recognizedForces['anglep'] = 1
|
|
recognizedForces['strbnd'] = 1
|
|
recognizedForces['ureybrad'] = 1
|
|
recognizedForces['opbend'] = 1
|
|
recognizedForces['torsion'] = 1
|
|
recognizedForces['pitors'] = 1
|
|
recognizedForces['strtors'] = 1
|
|
recognizedForces['angtors'] = 1
|
|
recognizedForces['vdw'] = 1
|
|
recognizedForces['vdwpair'] = 1
|
|
recognizedForces['polarize'] = 1
|
|
recognizedForces['tortors'] = addTorTor
|
|
recognizedForces['multipole'] = addMultipole
|
|
|
|
#=============================================================================================
|
|
|
|
# recognizedScalars[] contain raw scalar entries from TINKER parameter file
|
|
|
|
scalars = {}
|
|
recognizedScalars = {}
|
|
recognizedScalars['forcefield'] = '-2.55'
|
|
recognizedScalars['bond-cubic'] = '-2.55'
|
|
recognizedScalars['bond-quartic'] = '3.793125'
|
|
recognizedScalars['angle-cubic'] = '-0.014'
|
|
recognizedScalars['angle-quartic'] = '0.000056'
|
|
recognizedScalars['angle-pentic'] = '-0.0000007'
|
|
recognizedScalars['angle-sextic'] = '0.000000022'
|
|
recognizedScalars['opbendtype'] = 'ALLINGER'
|
|
recognizedScalars['opbend-cubic'] = '-0.014'
|
|
recognizedScalars['opbend-quartic'] = '0.000056'
|
|
recognizedScalars['opbend-pentic'] = '-0.0000007'
|
|
recognizedScalars['opbend-sextic'] = '0.000000022'
|
|
recognizedScalars['torsionunit'] = '0.5'
|
|
recognizedScalars['vdwtype'] = 'BUFFERED-14-7'
|
|
recognizedScalars['radiusrule'] = 'CUBIC-MEAN'
|
|
recognizedScalars['radiustype'] = 'R-MIN'
|
|
recognizedScalars['radiussize'] = 'DIAMETER'
|
|
recognizedScalars['epsilonrule'] = 'HHG'
|
|
recognizedScalars['dielectric'] = '1.0'
|
|
recognizedScalars['polarization'] = 'MUTUAL'
|
|
recognizedScalars['vdw-13-scale'] = '0.0'
|
|
recognizedScalars['vdw-14-scale'] = '1.0'
|
|
recognizedScalars['vdw-15-scale'] = '1.0'
|
|
recognizedScalars['mpole-12-scale'] = '0.0'
|
|
recognizedScalars['mpole-13-scale'] = '0.0'
|
|
recognizedScalars['mpole-14-scale'] = '0.4'
|
|
recognizedScalars['mpole-15-scale'] = '0.8'
|
|
recognizedScalars['polar-12-scale'] = '0.0'
|
|
recognizedScalars['polar-13-scale'] = '0.0'
|
|
recognizedScalars['polar-14-scale'] = '1.0'
|
|
recognizedScalars['polar-15-scale'] = '1.0'
|
|
recognizedScalars['polar-14-intra'] = '0.5'
|
|
recognizedScalars['direct-11-scale'] = '0.0'
|
|
recognizedScalars['direct-12-scale'] = '1.0'
|
|
recognizedScalars['direct-13-scale'] = '1.0'
|
|
recognizedScalars['direct-14-scale'] = '1.0'
|
|
recognizedScalars['mutual-11-scale'] = '1.0'
|
|
recognizedScalars['mutual-12-scale'] = '1.0'
|
|
recognizedScalars['mutual-13-scale'] = '1.0'
|
|
recognizedScalars['mutual-14-scale'] = '1.0'
|
|
|
|
#=============================================================================================
|
|
# get all 'interesting' lines in file
|
|
|
|
fileName = sys.argv[1]
|
|
if fileName == 'amoebabio18.prm':
|
|
ions = ions2018
|
|
residueXmlFileName = 'residuesFinalAmoeba2018.xml'
|
|
elif fileName == 'amoebapro13.prm':
|
|
ions = ions2013
|
|
residueXmlFileName = 'residuesFinalAmoeba2013.xml'
|
|
elif fileName == 'amoebabio09.prm':
|
|
ions = ions2009
|
|
residueXmlFileName = 'residuesFinalAmoeba2009.xml'
|
|
else:
|
|
print("Error: unrecognized force field file name %s" % (fileName))
|
|
sys.exit()
|
|
|
|
allLines = []
|
|
for line in open(fileName):
|
|
try:
|
|
fields = shlex.split(line)
|
|
except:
|
|
continue
|
|
if len(fields) == 0:
|
|
continue
|
|
if fields[0][0] == '#':
|
|
continue
|
|
allLines.append( fields )
|
|
|
|
#=============================================================================================
|
|
|
|
residueDict = buildResidueDict( residueXmlFileName )
|
|
|
|
#=============================================================================================
|
|
|
|
# load lines in lists/scalar values
|
|
|
|
lineIndex = 0
|
|
while lineIndex < len( allLines ):
|
|
|
|
fields = allLines[lineIndex]
|
|
|
|
if fields[0] == 'atom':
|
|
if( fields[3] in ions ):
|
|
ionInfo = ions[fields[3]]
|
|
element = ionInfo[0]
|
|
ionInfo[1] = int(fields[1])
|
|
else:
|
|
element = fields[3][0]
|
|
atomTypes[int(fields[1])] = (fields[2], element, fields[6])
|
|
lineIndex += 1
|
|
|
|
elif fields[0] == 'biotype':
|
|
|
|
lookUp = fields[2] + '_' + fields[3]
|
|
if lookUp in bioTypes:
|
|
# Workaround for Tinker using the same name but different types for H2', H2'', and for H5', H5''
|
|
lookUp = fields[2]+'*_'+fields[3]
|
|
bioTypes[lookUp] = fields[1:]
|
|
lineIndex += 1
|
|
|
|
elif fields[0] in recognizedForces:
|
|
if( recognizedForces[fields[0]] == 1 ):
|
|
if( fields[0] not in forces ):
|
|
forces[fields[0]] = []
|
|
forces[fields[0]].append( fields[1:] )
|
|
lineIndex += 1
|
|
else:
|
|
lineIndex = recognizedForces[fields[0]]( lineIndex, allLines, forces )
|
|
|
|
elif fields[0] in recognizedScalars:
|
|
scalars[fields[0]] = fields[1]
|
|
lineIndex += 1
|
|
else:
|
|
print("Field %s not recognized: line=<%s>" % ( fields[0], allLines[lineIndex] ))
|
|
lineIndex += 1
|
|
|
|
#=============================================================================================
|
|
|
|
# set biotypes for all atoms
|
|
|
|
setBioTypes( bioTypes, residueDict )
|
|
|
|
#=============================================================================================
|
|
|
|
# open force field xml file for output
|
|
|
|
tinkerXmlFileName = scalars['forcefield']
|
|
tinkerXmlFileName += '.xml'
|
|
tinkerXmlFile = open( tinkerXmlFileName, 'w' )
|
|
print("Opened %s." % (tinkerXmlFileName))
|
|
|
|
gkXmlFileName = scalars['forcefield']
|
|
gkXmlFileName += '_gk.xml'
|
|
gkXmlFile = open( gkXmlFileName, 'w' )
|
|
print("Opened %s." % (gkXmlFileName))
|
|
|
|
today = datetime.date.today().isoformat()
|
|
sourceFile = os.path.basename(sys.argv[1])
|
|
header = """ <Info>
|
|
<Source>%s</Source>
|
|
<DateGenerated>%s</DateGenerated>
|
|
<Reference></Reference>
|
|
</Info>
|
|
""" % (sourceFile, today)
|
|
|
|
gkXmlFile.write( "<ForceField>\n" )
|
|
gkXmlFile.write(header)
|
|
tinkerXmlFile.write( "<ForceField>\n" )
|
|
tinkerXmlFile.write(header)
|
|
tinkerXmlFile.write( " <AtomTypes>\n")
|
|
|
|
if( scalars['forcefield'].find( 'AMOEBA' ) > -1 ):
|
|
isAmoeba = 1
|
|
else:
|
|
isAmoeba = 0
|
|
|
|
#=============================================================================================
|
|
|
|
# atom type/class
|
|
|
|
# atmType class name name atomicNo. mass valence
|
|
# atom 1 1 N "Glycine N" 7 14.003 3
|
|
# atom 2 2 CA "Glycine CA" 6 12.000 4
|
|
# atom 3 3 C "Glycine C" 6 12.000 3
|
|
# atom 4 4 HN "Glycine HN" 1 1.008 1
|
|
# atom 5 5 O "Glycine O" 8 15.995 1
|
|
|
|
# atom 380 73 O "AMOEBA Water O" 8 15.999 2
|
|
# atom 381 74 H "AMOEBA Water H" 1 1.008 1
|
|
# atom 383 76 Na+ "Sodium Ion Na+" 11 22.990 0
|
|
# atom 384 77 K+ "Potassium Ion K+" 19 39.098 0
|
|
# atom 385 78 Rb+ "Rubidium Ion Rb+" 37 85.468 0
|
|
# atom 386 79 Cs+ "Cesium Ion Cs+" 55 132.905 0
|
|
# atom 387 80 Be+ "Beryllium Ion Be+2" 4 9.012 0
|
|
# atom 388 81 Mg+ "Magnesium Ion Mg+2" 12 24.305 0
|
|
# atom 389 82 Ca+ "Calcium Ion Ca+2" 20 40.078 0
|
|
# atom 390 83 Cl- "Chloride Ion Cl-" 17 35.453 0
|
|
|
|
|
|
# biotype atmType
|
|
# biotype 1 N "Glycine" 1
|
|
# biotype 2 CA "Glycine" 2
|
|
# biotype 3 C "Glycine" 3
|
|
# biotype 4 HN "Glycine" 4
|
|
# biotype 5 O "Glycine" 5
|
|
|
|
# biotype 2001 O "Water" 380
|
|
# biotype 2002 H "Water" 381
|
|
# biotype 2003 NA "Sodium Ion" 383
|
|
# biotype 2004 K "Potassium Ion" 384
|
|
# biotype 2005 MG "Magnesium Ion" 388
|
|
# biotype 2006 CA "Calcium Ion" 389
|
|
# biotype 2007 CL "Chloride Ion" 390
|
|
|
|
if( isAmoeba ):
|
|
for type in sorted(atomTypes):
|
|
outputString = """ <Type name="%s" class="%s" element="%s" mass="%s"/>""" % (type, atomTypes[type][0], atomTypes[type][1], atomTypes[type][2])
|
|
tinkerXmlFile.write( "%s\n" % (outputString) )
|
|
else:
|
|
for type in sorted(atomTypes):
|
|
outputString = """ <Type name="%s" class="%s" mass="%s"/>""" % (type, atomTypes[type][0], atomTypes[type][1])
|
|
tinkerXmlFile.write( "%s\n" % (outputString) )
|
|
|
|
tinkerXmlFile.write( " </AtomTypes>\n")
|
|
|
|
#=============================================================================================
|
|
|
|
# residues
|
|
|
|
tinkerXmlFile.write( " <Residues>\n" )
|
|
for resname, res in sorted(residueDict.items()):
|
|
if res['include']:
|
|
if resname == 'HOH':
|
|
tinkerXmlFile.write(f' <Residue name="{resname}" rigidWater="false">\n')
|
|
else:
|
|
tinkerXmlFile.write(f' <Residue name="{resname}">\n')
|
|
atomIndex = dict()
|
|
atomCount = 0
|
|
for atom in sorted( res['atoms'].keys() ):
|
|
type = res['atoms'][atom]['type']
|
|
typeI = int( type )
|
|
if( typeI < 0 ):
|
|
print("Error: type=%s for atom=%s of residue=%s" % (type, atom, resname))
|
|
tag = " <Atom name=\"%s\" type=\"%s\" />" % (atom, type)
|
|
atomIndex[atom] = atomCount
|
|
atomCount += 1
|
|
tinkerXmlFile.write( "%s\n" % (tag) )
|
|
|
|
includedBonds = dict()
|
|
for atomName in sorted( res['atoms'].keys() ):
|
|
bondsInfo = res['atoms'][atomName]['bonds']
|
|
for bondedAtom in bondsInfo:
|
|
if( bondedAtom not in includedBonds or atomName not in includedBonds[bondedAtom] ):
|
|
outputString = """ <Bond from="%s" to="%s" />""" % (str(atomIndex[atomName]), str(atomIndex[bondedAtom]))
|
|
if( atomName not in includedBonds ):
|
|
includedBonds[atomName] = dict()
|
|
if( bondedAtom not in includedBonds ):
|
|
includedBonds[bondedAtom] = dict()
|
|
includedBonds[atomName][bondedAtom] = 1
|
|
includedBonds[bondedAtom][atomName] = 1
|
|
tinkerXmlFile.write( "%s\n" % (outputString) )
|
|
|
|
outputStrings = []
|
|
if res['type'] in ('rna', 'dna'):
|
|
if res['loc'] in ('middle', '3'):
|
|
outputStrings.append( """ <ExternalBond from="%s" />""" % (str(atomIndex['P']) ))
|
|
if res['loc'] in ('middle', '5'):
|
|
outputStrings.append( """ <ExternalBond from="%s" />""" % (str(atomIndex["O3'"]) ))
|
|
else:
|
|
if res['loc'] == 'm':
|
|
outputStrings.append( """ <ExternalBond from="%s" />""" % (str(atomIndex['N']) ))
|
|
outputStrings.append( """ <ExternalBond from="%s" />""" % (str(atomIndex['C']) ) )
|
|
if res['loc'] == 'n':
|
|
outputStrings.append( """ <ExternalBond from="%s" />""" % (str(atomIndex['C']) ) )
|
|
if res['loc'] == 'c':
|
|
outputStrings.append( """ <ExternalBond from="%s" />""" % (str(atomIndex['N']) ) )
|
|
if resname.find('CYX') > -1:
|
|
outputStrings.append( """ <ExternalBond from="%s" />""" % (str(atomIndex['SG']) ) )
|
|
for outputString in outputStrings:
|
|
tinkerXmlFile.write( "%s\n" % (outputString) )
|
|
|
|
tinkerXmlFile.write( " </Residue>\n" )
|
|
|
|
# End caps
|
|
|
|
tinkerXmlFile.write(' <Residue name="ACE">\n')
|
|
tinkerXmlFile.write(' <Atom name="HH31" type="%s"/>\n' % bioTypes['H_Acetyl N-Terminus'][3])
|
|
tinkerXmlFile.write(' <Atom name="CH3" type="%s"/>\n' % bioTypes['CH3_Acetyl N-Terminus'][3])
|
|
tinkerXmlFile.write(' <Atom name="HH32" type="%s"/>\n' % bioTypes['H_Acetyl N-Terminus'][3])
|
|
tinkerXmlFile.write(' <Atom name="HH33" type="%s"/>\n' % bioTypes['H_Acetyl N-Terminus'][3])
|
|
tinkerXmlFile.write(' <Atom name="C" type="%s"/>\n' % bioTypes['C_Acetyl N-Terminus'][3])
|
|
tinkerXmlFile.write(' <Atom name="O" type="%s"/>\n' % bioTypes['O_Acetyl N-Terminus'][3])
|
|
tinkerXmlFile.write(' <Bond from="0" to="1"/>\n')
|
|
tinkerXmlFile.write(' <Bond from="1" to="2"/>\n')
|
|
tinkerXmlFile.write(' <Bond from="1" to="3"/>\n')
|
|
tinkerXmlFile.write(' <Bond from="1" to="4"/>\n')
|
|
tinkerXmlFile.write(' <Bond from="4" to="5"/>\n')
|
|
tinkerXmlFile.write(' <ExternalBond from="4"/>\n')
|
|
tinkerXmlFile.write(' </Residue>\n')
|
|
tinkerXmlFile.write(' <Residue name="NME">\n')
|
|
tinkerXmlFile.write(' <Atom name="N" type="%s"/>\n' % bioTypes['N_N-MeAmide C-Terminus'][3])
|
|
tinkerXmlFile.write(' <Atom name="H" type="%s"/>\n' % bioTypes['HN_N-MeAmide C-Terminus'][3])
|
|
tinkerXmlFile.write(' <Atom name="CH3" type="%s"/>\n' % bioTypes['C_N-MeAmide C-Terminus'][3])
|
|
tinkerXmlFile.write(' <Atom name="HH31" type="%s"/>\n' % bioTypes['H_N-MeAmide C-Terminus'][3])
|
|
tinkerXmlFile.write(' <Atom name="HH32" type="%s"/>\n' % bioTypes['H_N-MeAmide C-Terminus'][3])
|
|
tinkerXmlFile.write(' <Atom name="HH33" type="%s"/>\n' % bioTypes['H_N-MeAmide C-Terminus'][3])
|
|
tinkerXmlFile.write(' <Bond from="0" to="1"/>\n')
|
|
tinkerXmlFile.write(' <Bond from="0" to="2"/>\n')
|
|
tinkerXmlFile.write(' <Bond from="2" to="3"/>\n')
|
|
tinkerXmlFile.write(' <Bond from="2" to="4"/>\n')
|
|
tinkerXmlFile.write(' <Bond from="2" to="5"/>\n')
|
|
tinkerXmlFile.write(' <ExternalBond from="0"/>\n')
|
|
tinkerXmlFile.write(' </Residue>\n')
|
|
|
|
|
|
# ions
|
|
|
|
for ion,ionInfo in ions.items():
|
|
outputString = """ <Residue name="%s">\n""" % (ionInfo[0])
|
|
outputString += """ <Atom name="%s" type="%s"/>\n""" % (ionInfo[0], str(ionInfo[1]))
|
|
outputString += """ </Residue>\n"""
|
|
tinkerXmlFile.write( "%s" % (outputString) )
|
|
|
|
tinkerXmlFile.write( " </Residues>\n" )
|
|
|
|
radian = 57.2957795130
|
|
if( isAmoeba ):
|
|
|
|
#=============================================================================================
|
|
|
|
# AmoebaBondForce
|
|
|
|
cubic = 10.*float(scalars['bond-cubic'])
|
|
quartic = 100.*float(scalars['bond-quartic'])
|
|
outputString = """ <AmoebaBondForce bond-cubic="%s" bond-quartic="%s">""" %( str(cubic), str(quartic) )
|
|
tinkerXmlFile.write( "%s\n" % (outputString ) )
|
|
bonds = forces['bond']
|
|
for bond in bonds:
|
|
length = float(bond[3])*0.1
|
|
k = float(bond[2])*100.0*4.184
|
|
outputString = """ <Bond class1="%s" class2="%s" length="%s" k="%s"/>""" % (bond[0], bond[1], xml_float(length), xml_float(k) )
|
|
tinkerXmlFile.write( "%s\n" % (outputString ) )
|
|
tinkerXmlFile.write( " </AmoebaBondForce>\n" )
|
|
|
|
#=============================================================================================
|
|
|
|
# AmoebaAngleForce
|
|
|
|
cubic = float(scalars['angle-cubic'])
|
|
quartic = float(scalars['angle-quartic'])
|
|
pentic = float(scalars['angle-pentic'])
|
|
sextic = float(scalars['angle-sextic'])
|
|
outputString = """ <AmoebaAngleForce angle-cubic="%s" angle-quartic="%s" angle-pentic="%s" angle-sextic="%s">""" %( str(cubic), str(quartic), str(pentic), str(sextic) )
|
|
tinkerXmlFile.write( "%s\n" % (outputString ) )
|
|
radian = 57.2957795130
|
|
radian2 = 4.184/(radian*radian)
|
|
for set in ['angle', 'anglep']:
|
|
for angle in forces[set]:
|
|
k = float(angle[3])*radian2
|
|
outputString = ' <Angle class1="%s" class2="%s" class3="%s" k="%s" angle1="%s"' % (angle[0], angle[1], angle[2], xml_float(k), xml_float(angle[4]) )
|
|
if( len(angle) > 5 ):
|
|
outputString += ' angle2="%s"' % (angle[5])
|
|
|
|
if( len(angle) > 6 ):
|
|
outputString += ' angle3="%s"' % (angle[6])
|
|
outputString += ' inPlane="%s"/>' % (set == 'anglep')
|
|
|
|
tinkerXmlFile.write( "%s\n" % (outputString ) )
|
|
tinkerXmlFile.write( " </AmoebaAngleForce>\n" )
|
|
|
|
#=============================================================================================
|
|
|
|
# AmoebaOutOfPlaneBendForce
|
|
|
|
cubic = float(scalars['opbend-cubic'])
|
|
quartic = float(scalars['opbend-quartic'])
|
|
pentic = float(scalars['opbend-pentic'])
|
|
sextic = float(scalars['opbend-sextic'])
|
|
type = scalars['opbendtype']
|
|
outputString = """ <AmoebaOutOfPlaneBendForce type="%s" opbend-cubic="%s" opbend-quartic="%s" opbend-pentic="%s" opbend-sextic="%s">""" % (
|
|
type, str(cubic), str(quartic), str(pentic), str(sextic) )
|
|
tinkerXmlFile.write( "%s\n" % (outputString ) )
|
|
opbends = forces['opbend']
|
|
radian2 = 4.184/(radian*radian)
|
|
for opbend in opbends:
|
|
k = float(opbend[4])*radian2
|
|
for i in range(4):
|
|
if opbend[i] == '0':
|
|
opbend[i] = ''
|
|
outputString = """ <Angle class1="%s" class2="%s" class3="%s" class4="%s" k="%s"/>""" % (opbend[0], opbend[1], opbend[2], opbend[3], xml_float(k))
|
|
tinkerXmlFile.write( "%s\n" % (outputString ) )
|
|
tinkerXmlFile.write( " </AmoebaOutOfPlaneBendForce>\n" )
|
|
|
|
#=============================================================================================
|
|
|
|
# AmoebaTorsionForce
|
|
|
|
torsionUnit = float(scalars['torsionunit'])
|
|
outputString = """ <PeriodicTorsionForce>"""
|
|
tinkerXmlFile.write( "%s\n" % (outputString ) )
|
|
torsions = forces['torsion']
|
|
conversion = 4.184*torsionUnit
|
|
for torsion in torsions:
|
|
outputString = """ <Proper class1="%s" class2="%s" class3="%s" class4="%s" """ % (torsion[0], torsion[1], torsion[2], torsion[3])
|
|
startIndex = 4
|
|
for ii in range(0,3):
|
|
torsionSuffix = str(ii+1)
|
|
amplitudeAttributeName = 'k'+torsionSuffix
|
|
angleAttributeName = 'phase'+torsionSuffix
|
|
periodicityAttributeName = 'periodicity'+torsionSuffix
|
|
amplitude = float(torsion[startIndex])*conversion
|
|
angle = float(torsion[startIndex+1])/radian
|
|
periodicity = int(torsion[startIndex+2])
|
|
outputString += """ %s="%s" %s="%s" %s="%d" """ % (amplitudeAttributeName, xml_float(amplitude), angleAttributeName, xml_float(angle), periodicityAttributeName, periodicity)
|
|
startIndex += 3
|
|
outputString += "/>"
|
|
tinkerXmlFile.write( "%s\n" % (outputString ) )
|
|
tinkerXmlFile.write( " </PeriodicTorsionForce>\n" )
|
|
|
|
#=============================================================================================
|
|
|
|
# AmoebaPiTorsionForce
|
|
|
|
piTorsionUnit = 1.0
|
|
outputString = """ <AmoebaPiTorsionForce piTorsionUnit="%s">""" % ( piTorsionUnit )
|
|
tinkerXmlFile.write( "%s\n" % (outputString ) )
|
|
piTorsions = forces['pitors']
|
|
conversion = 4.184*piTorsionUnit
|
|
for piTorsion in piTorsions:
|
|
k = float(piTorsion[2])*conversion
|
|
outputString = """ <PiTorsion class1="%s" class2="%s" k="%s" />""" % (piTorsion[0], piTorsion[1], xml_float(k))
|
|
tinkerXmlFile.write( "%s\n" % (outputString ) )
|
|
tinkerXmlFile.write( " </AmoebaPiTorsionForce>\n" )
|
|
|
|
#=============================================================================================
|
|
|
|
# Stretch torsion
|
|
if 'strtors' in forces:
|
|
tinkerXmlFile.write(' <AmoebaStretchTorsionForce>\n')
|
|
for torsion in forces['strtors']:
|
|
v = [float(x)*10*4.184 for x in torsion[4:]]
|
|
tinkerXmlFile.write(f' <Torsion class1="{torsion[0]}" class2="{torsion[1]}" class3="{torsion[2]}" class4="{torsion[3]}" v11="{xml_float(v[0])}" v12="{xml_float(v[1])}" v13="{xml_float(v[2])}" v21="{xml_float(v[3])}" v22="{xml_float(v[4])}" v23="{xml_float(v[5])}" v31="{xml_float(v[6])}" v32="{xml_float(v[7])}" v33="{xml_float(v[8])}"/>\n')
|
|
tinkerXmlFile.write(' </AmoebaStretchTorsionForce>\n')
|
|
|
|
#=============================================================================================
|
|
|
|
# Angle torsion
|
|
if 'angtors' in forces:
|
|
tinkerXmlFile.write(' <AmoebaAngleTorsionForce>\n')
|
|
for torsion in forces['angtors']:
|
|
v = [float(x)*4.184 for x in torsion[4:]]
|
|
tinkerXmlFile.write(f' <Torsion class1="{torsion[0]}" class2="{torsion[1]}" class3="{torsion[2]}" class4="{torsion[3]}" v11="{xml_float(v[0])}" v12="{xml_float(v[1])}" v13="{xml_float(v[2])}" v21="{xml_float(v[3])}" v22="{xml_float(v[4])}" v23="{xml_float(v[5])}"/>\n')
|
|
tinkerXmlFile.write(' </AmoebaAngleTorsionForce>\n')
|
|
|
|
#=============================================================================================
|
|
|
|
# AmoebaStretchBendForce
|
|
|
|
stretchBendUnit = 1.0
|
|
outputString = """ <AmoebaStretchBendForce stretchBendUnit="%s">""" % ( stretchBendUnit )
|
|
tinkerXmlFile.write( "%s\n" % (outputString ) )
|
|
conversion = 41.84/radian
|
|
stretchBends = forces['strbnd']
|
|
for stretchBend in stretchBends:
|
|
k1 = float(stretchBend[3])*conversion
|
|
k2 = float(stretchBend[4])*conversion
|
|
outputString = """ <StretchBend class1="%s" class2="%s" class3="%s" k1="%s" k2="%s" />""" % (stretchBend[0], stretchBend[1], stretchBend[2], xml_float(k1), xml_float(k2))
|
|
tinkerXmlFile.write( "%s\n" % (outputString ) )
|
|
tinkerXmlFile.write( "</AmoebaStretchBendForce>\n" )
|
|
|
|
#=============================================================================================
|
|
|
|
# AmoebaTorsionTorsionForce
|
|
|
|
torsionTorsionUnit = 1.0
|
|
outputString = """ <AmoebaTorsionTorsionForce >"""
|
|
tinkerXmlFile.write( "%s\n" % (outputString ) )
|
|
torsionTorsions = forces['tortors']
|
|
for (index, torsionTorsion) in enumerate(torsionTorsions):
|
|
torInfo = torsionTorsion[0]
|
|
grid = torsionTorsion[1]
|
|
outputString = """ <TorsionTorsion class1="%s" class2="%s" class3="%s" class4="%s" class5="%s" grid="%s" nx="%s" ny="%s" />""" % (torInfo[0], torInfo[1], torInfo[2], torInfo[3], torInfo[4], str(index),
|
|
torInfo[5], torInfo[6] )
|
|
tinkerXmlFile.write( "%s\n" % (outputString ) )
|
|
|
|
for (index, torsionTorsion) in enumerate(torsionTorsions):
|
|
torInfo = torsionTorsion[0]
|
|
grid = torsionTorsion[1]
|
|
outputString = """ <TorsionTorsionGrid grid="%s" nx="%s" ny="%s" >""" % (str(index), torInfo[5], torInfo[6] )
|
|
tinkerXmlFile.write( "%s\n" % (outputString ) )
|
|
for (gridIndex, gridEntry) in enumerate(grid):
|
|
print("Gxx %d %s" % ( gridIndex, str(gridEntry) ))
|
|
if( len( gridEntry ) > 5 ):
|
|
f = float( gridEntry[2] )*4.184
|
|
fx = float( gridEntry[3] )*4.184
|
|
fy = float( gridEntry[4] )*4.184
|
|
fxy = float( gridEntry[5] )*4.184
|
|
outputString = """ <Grid angle1="%s" angle2="%s" f="%s" fx="%s" fy="%s" fxy="%s" />""" % ( gridEntry[0], gridEntry[1], xml_float(f), xml_float(fx), xml_float(fy), xml_float(fxy) )
|
|
tinkerXmlFile.write( " %s\n" % (outputString ) )
|
|
elif( len( gridEntry ) > 2 ):
|
|
f = float( gridEntry[2] )*4.184
|
|
outputString = """ <Grid angle1="%s" angle2="%s" f="%s" />""" % ( gridEntry[0], gridEntry[1], xml_float(f) )
|
|
tinkerXmlFile.write( " %s\n" % (outputString ) )
|
|
outputString = '</TorsionTorsionGrid >'
|
|
tinkerXmlFile.write( "%s\n" % (outputString ) )
|
|
|
|
tinkerXmlFile.write( "</AmoebaTorsionTorsionForce>\n" )
|
|
|
|
#=============================================================================================
|
|
|
|
# AmoebaVdwForce
|
|
|
|
outputString = """ <AmoebaVdwForce type="%s" radiusrule="%s" radiustype="%s" radiussize="%s" epsilonrule="%s" vdw-13-scale="%s" vdw-14-scale="%s" vdw-15-scale="%s" >""" % (
|
|
scalars['vdwtype'], scalars['radiusrule'], scalars['radiustype'], scalars['radiussize'], scalars['epsilonrule'], scalars['vdw-13-scale'], scalars['vdw-14-scale'], scalars['vdw-15-scale'] )
|
|
tinkerXmlFile.write( "%s\n" % (outputString ) )
|
|
for vdw in forces['vdw']:
|
|
sigma = float(vdw[1])*0.1
|
|
epsilon = float(vdw[2])*4.184
|
|
if( len(vdw) > 3 ):
|
|
reduction = vdw[3]
|
|
else:
|
|
reduction = 1.0
|
|
outputString = """ <Vdw class="%s" sigma="%s" epsilon="%s" reduction="%s" />""" % (vdw[0], xml_float(sigma), xml_float(epsilon), xml_float(reduction))
|
|
tinkerXmlFile.write( "%s\n" % (outputString ) )
|
|
for pair in forces['vdwpair']:
|
|
sigma = float(pair[2])*0.1
|
|
epsilon = float(pair[3])*4.184
|
|
outputString = """ <Pair class1="%s" class2="%s" sigma="%s" epsilon="%s"/>""" % (pair[0], pair[1], xml_float(sigma), xml_float(epsilon))
|
|
tinkerXmlFile.write( "%s\n" % (outputString ) )
|
|
tinkerXmlFile.write( " </AmoebaVdwForce>\n" )
|
|
|
|
#=============================================================================================
|
|
|
|
# AmoebaMultipoleForce
|
|
|
|
scalarList = dict()
|
|
scalarList['mpole12Scale'] = recognizedScalars['mpole-12-scale']
|
|
scalarList['mpole13Scale'] = recognizedScalars['mpole-13-scale']
|
|
scalarList['mpole14Scale'] = recognizedScalars['mpole-14-scale']
|
|
scalarList['mpole15Scale'] = recognizedScalars['mpole-15-scale']
|
|
|
|
scalarList['polar12Scale'] = recognizedScalars['polar-12-scale']
|
|
scalarList['polar13Scale'] = recognizedScalars['polar-13-scale']
|
|
scalarList['polar14Scale'] = recognizedScalars['polar-14-scale']
|
|
scalarList['polar15Scale'] = recognizedScalars['polar-15-scale']
|
|
scalarList['polar14Intra'] = recognizedScalars['polar-14-intra']
|
|
|
|
scalarList['direct11Scale'] = recognizedScalars['direct-11-scale']
|
|
scalarList['direct12Scale'] = recognizedScalars['direct-12-scale']
|
|
scalarList['direct13Scale'] = recognizedScalars['direct-13-scale']
|
|
scalarList['direct14Scale'] = recognizedScalars['direct-14-scale']
|
|
|
|
scalarList['mutual11Scale'] = recognizedScalars['mutual-11-scale']
|
|
scalarList['mutual12Scale'] = recognizedScalars['mutual-12-scale']
|
|
scalarList['mutual13Scale'] = recognizedScalars['mutual-13-scale']
|
|
scalarList['mutual14Scale'] = recognizedScalars['mutual-14-scale']
|
|
|
|
outputString = """ <AmoebaMultipoleForce """
|
|
for key in sorted( scalarList.keys() ):
|
|
outputString += """ %s="%s" """ % ( key, scalarList[key] )
|
|
outputString += """ > """
|
|
tinkerXmlFile.write( "%s\n" % (outputString ) )
|
|
|
|
multipoleArray = forces['multipole']
|
|
bohr = 0.52917720859
|
|
dipoleConversion = 0.1*bohr
|
|
quadrupoleConversion = 0.01*bohr*bohr/3.0
|
|
for multipoleInfo in multipoleArray:
|
|
axisInfo = multipoleInfo[0]
|
|
multipoles = multipoleInfo[1]
|
|
outputString = """ <Multipole type="%s" """ % (axisInfo[0] )
|
|
axisInfoLen = len(axisInfo)
|
|
|
|
if( axisInfoLen > 1 ):
|
|
outputString += """kz="%s" """ % ( axisInfo[1] )
|
|
|
|
if( axisInfoLen > 2 ):
|
|
outputString += """kx="%s" """ % ( axisInfo[2] )
|
|
|
|
if( axisInfoLen > 3 ):
|
|
outputString += """ky="%s" """ % ( axisInfo[3] )
|
|
|
|
outputString += """c0="%s" d1="%s" d2="%s" d3="%s" q11="%s" q21="%s" q22="%s" q31="%s" q32="%s" q33="%s" """ % ( xml_float(multipoles[0]),
|
|
xml_float(dipoleConversion*float(multipoles[1])), xml_float(dipoleConversion*float(multipoles[2])), xml_float(dipoleConversion*float(multipoles[3])),
|
|
xml_float(quadrupoleConversion*float(multipoles[4])), xml_float(quadrupoleConversion*float(multipoles[5])), xml_float(quadrupoleConversion*float(multipoles[6])),
|
|
xml_float(quadrupoleConversion*float(multipoles[7])), xml_float(quadrupoleConversion*float(multipoles[8])), xml_float(quadrupoleConversion*float(multipoles[9])) )
|
|
outputString += "/>"
|
|
tinkerXmlFile.write( "%s\n" % (outputString ) )
|
|
|
|
polarizeArray = forces['polarize']
|
|
|
|
polarityConversion = 0.001
|
|
m = {}
|
|
for polarize in polarizeArray:
|
|
m[polarize[0]] = []
|
|
outputString = """ <Polarize type="%s" polarizability="%s" thole="%s" """ % (polarize[0], xml_float(polarityConversion*float(polarize[1])), xml_float(float(polarize[2])) )
|
|
for ii in range( 3, len(polarize) ):
|
|
outputString += """pgrp%d="%s" """ % (ii-2, polarize[ii])
|
|
m[polarize[0]].append(polarize[ii])
|
|
|
|
outputString += "/>"
|
|
tinkerXmlFile.write( "%s\n" % (outputString ) )
|
|
print(m[polarize[0]])
|
|
for t in sorted(m):
|
|
for k in m[t]:
|
|
if t not in m[k]:
|
|
print(t, k)
|
|
|
|
tinkerXmlFile.write( " </AmoebaMultipoleForce>\n" )
|
|
|
|
#=============================================================================================
|
|
|
|
# AmoebaGeneralizedKirkwoodForce
|
|
|
|
solventDielectric = 78.3
|
|
soluteDielectric = 1.0
|
|
includeCavityTerm = 1
|
|
probeRadius = 0.14
|
|
surfaceAreaFactor = -6.0*3.1415926535*0.0216*1000.0*0.4184
|
|
outputString = """ <AmoebaGeneralizedKirkwoodForce solventDielectric="%s" soluteDielectric="%s" includeCavityTerm="%s" probeRadius="%s" surfaceAreaFactor="%s">""" % (
|
|
solventDielectric, soluteDielectric, includeCavityTerm, probeRadius, surfaceAreaFactor )
|
|
gkXmlFile.write( "%s\n" % (outputString ) )
|
|
|
|
# radii are set in forcefield.py
|
|
|
|
for type in sorted( atomTypes ):
|
|
print("atom type=%s %s" % ( str(type), str(atomTypes[type]) ))
|
|
|
|
for type in sorted( bioTypes ):
|
|
print("bio type=%s %s" % ( str(type), str(bioTypes[type]) ))
|
|
|
|
multipoleArray = forces['multipole']
|
|
for multipoleInfo in multipoleArray:
|
|
axisInfo = multipoleInfo[0]
|
|
multipoles = multipoleInfo[1]
|
|
type = int(axisInfo[0])
|
|
shct = 0.8
|
|
if( type in atomTypes ):
|
|
element = atomTypes[type][1]
|
|
if( element == 'H' ):
|
|
shct = 0.85
|
|
elif( element == 'C' ):
|
|
shct = 0.72
|
|
elif( element == 'N' ):
|
|
shct = 0.79
|
|
elif( element == 'O' ):
|
|
shct = 0.85
|
|
elif( element == 'F' ):
|
|
shct = 0.88
|
|
elif( element == 'P' ):
|
|
shct = 0.86
|
|
elif( element == 'S' ):
|
|
shct = 0.96
|
|
elif( element == 'Fe' ):
|
|
shct = 0.88
|
|
else:
|
|
print("Warning no overlap scale factor for type=%d element=%s" % (type, element))
|
|
else:
|
|
print("Warning no overlap scale factor for type=%d " % (type))
|
|
|
|
outputString = """ <GeneralizedKirkwood type="%s" charge="%s" shct="%s" /> """ % ( axisInfo[0], xml_float(multipoles[0]), xml_float(shct) )
|
|
gkXmlFile.write( "%s\n" % (outputString ) )
|
|
gkXmlFile.write( " </AmoebaGeneralizedKirkwoodForce>\n" )
|
|
|
|
#=============================================================================================
|
|
|
|
# AmoebaWcaDispersionForce
|
|
|
|
epso = 0.1100
|
|
epsh = 0.0135
|
|
rmino = 1.7025
|
|
rminh = 1.3275
|
|
awater = 0.033428
|
|
slevy = 1.0
|
|
dispoff = 0.26
|
|
shctd = 0.81
|
|
|
|
outputString = """ <AmoebaWcaDispersionForce epso="%s" epsh="%s" rmino="%s" rminh="%s" awater="%s" slevy="%s" dispoff="%s" shctd="%s" >""" % (
|
|
epso*4.184, epsh*4.184, rmino*0.1, rminh*0.1, 1000.0*awater, slevy, 0.1*dispoff, shctd )
|
|
gkXmlFile.write( "%s\n" % (outputString ) )
|
|
vdws = forces['vdw']
|
|
convert = 0.1
|
|
if( scalars['radiustype'] == 'SIGMA' ):
|
|
convert *= 1.122462048309372
|
|
|
|
if( scalars['radiussize'] == 'DIAMETER' ):
|
|
convert *= 0.5
|
|
|
|
for vdw in vdws:
|
|
sigma = float(vdw[1])
|
|
sigma *= convert
|
|
epsilon = float(vdw[2])*4.184
|
|
outputString = """ <WcaDispersion class="%s" radius="%s" epsilon="%s" /> """ % ( vdw[0], xml_float(sigma), xml_float(epsilon) )
|
|
gkXmlFile.write( "%s\n" % (outputString ) )
|
|
gkXmlFile.write( " </AmoebaWcaDispersionForce>\n" )
|
|
|
|
#=============================================================================================
|
|
|
|
# AmoebaUreyBradleyForce
|
|
|
|
cubic = 0.0
|
|
quartic = 0.0
|
|
|
|
outputString = """ <AmoebaUreyBradleyForce cubic="%s" quartic="%s" >""" % ( cubic, quartic )
|
|
tinkerXmlFile.write( "%s\n" % (outputString ) )
|
|
ubs = forces['ureybrad']
|
|
for ub in ubs:
|
|
k = float(ub[3])*4.184*100.0
|
|
d = float(ub[4])*0.1
|
|
outputString = """ <UreyBradley class1="%s" class2="%s" class3="%s" k="%s" d="%s" /> """ % ( ub[0], ub[1], ub[2], xml_float(k), xml_float(d) )
|
|
tinkerXmlFile.write( "%s\n" % (outputString ) )
|
|
tinkerXmlFile.write( " </AmoebaUreyBradleyForce>\n" )
|
|
|
|
#=============================================================================================
|
|
|
|
tinkerXmlFile.write("</ForceField>\n")
|
|
gkXmlFile.write("</ForceField>\n")
|
|
tinkerXmlFile.close()
|
|
gkXmlFile.close()
|
|
|