RustOpenMM/openmm
1
0
Fork 0
mirror of https://github.com/openmm/openmm synced 2026-06-03 06:39:48 +09:00
Code Issues Packages Projects Releases Wiki Activity
Files
master
openmm/devtools/forcefield-scripts/processTinkerForceField.py
Christian Clauss 82703dff73 Fix typos discovered by codespell in Python files (#5173)
2025-12-30 12:21:00 -08:00

1253 lines
55 KiB
Python
Raw Permalink Blame History

#=============================================================================================
# 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()
Reference in New Issue View Git Blame Copy Permalink