Skip to content

P
pdbfixer
Commit 6da5bb67 by Peter Eastman Committed by GitHub
Options

Soft core force field includes torsions with wildcards (#308) 

* Soft core force field includes torsions with wildcards

* Bug fixes
parent 2744212e
Showing with 23 additions and 29 deletions
devtools/createSoftForcefield.py
...@@ -51,18 +51,21 @@ print(' <AtomTypes>') ...@@ -51,18 +51,21 @@ print(' <AtomTypes>')
omitTypes = set() omitTypes = set()
omitClasses = set() omitClasses = set()
for atomType in forcefield._atomTypes: for atomType in forcefield._atomTypes:
(atomClass, mass, element) = forcefield._atomTypes[atomType] data = forcefield._atomTypes[atomType]
if element is None or element == elem.hydrogen: if data.element is None or data.element == elem.hydrogen:
omitTypes.add(atomType) omitTypes.add(atomType)
omitClasses.add(atomClass) omitClasses.add(data.atomClass)
else: else:
print(' <Type name="%s" class="%s" element="%s" mass="%g"/>' % (atomType, atomClass, element.symbol, mass)) print(' <Type name="%s" class="%s" element="%s" mass="%g"/>' % (atomType, data.atomClass, data.element.symbol, data.mass))
print(' </AtomTypes>') print(' </AtomTypes>')
# Print the residue templates. # Print the residue templates.
print(' <Residues>') print(' <Residues>')
skipTemplates = ('ASH', 'CYM', 'GLH', 'HID', 'HIE', 'LYN')
for template in forcefield._templates.values(): for template in forcefield._templates.values():
if template.name in skipTemplates or template.name[1:] in skipTemplates:
continue
print(' <Residue name="%s">' % template.name) print(' <Residue name="%s">' % template.name)
atomIndex = {} atomIndex = {}
for i, atom in enumerate(template.atoms): for i, atom in enumerate(template.atoms):
...@@ -86,8 +89,8 @@ for i in range(len(bonds.types1)): ...@@ -86,8 +89,8 @@ for i in range(len(bonds.types1)):
type1 = next(iter(bonds.types1[i])) type1 = next(iter(bonds.types1[i]))
type2 = next(iter(bonds.types2[i])) type2 = next(iter(bonds.types2[i]))
if type1 not in omitTypes and type2 not in omitTypes: if type1 not in omitTypes and type2 not in omitTypes:
class1 = forcefield._atomTypes[type1][0] class1 = forcefield._atomTypes[type1].atomClass
class2 = forcefield._atomTypes[type2][0] class2 = forcefield._atomTypes[type2].atomClass
print(' <Bond class1="%s" class2="%s" length="%g" k="%g"/>' % (class1, class2, bonds.length[i], bondK)) print(' <Bond class1="%s" class2="%s" length="%g" k="%g"/>' % (class1, class2, bonds.length[i], bondK))
print(' </HarmonicBondForce>') print(' </HarmonicBondForce>')
...@@ -100,9 +103,9 @@ for i in range(len(angles.types1)): ...@@ -100,9 +103,9 @@ for i in range(len(angles.types1)):
type2 = next(iter(angles.types2[i])) type2 = next(iter(angles.types2[i]))
type3 = next(iter(angles.types3[i])) type3 = next(iter(angles.types3[i]))
if type1 not in omitTypes and type2 not in omitTypes and type3 not in omitTypes: if type1 not in omitTypes and type2 not in omitTypes and type3 not in omitTypes:
class1 = forcefield._atomTypes[type1][0] class1 = forcefield._atomTypes[type1].atomClass
class2 = forcefield._atomTypes[type2][0] class2 = forcefield._atomTypes[type2].atomClass
class3 = forcefield._atomTypes[type3][0] class3 = forcefield._atomTypes[type3].atomClass
print(' <Angle class1="%s" class2="%s" class3="%s" angle="%g" k="%g"/>' % (class1, class2, class3, angles.angle[i], angleK)) print(' <Angle class1="%s" class2="%s" class3="%s" angle="%g" k="%g"/>' % (class1, class2, class3, angles.angle[i], angleK))
print(' </HarmonicAngleForce>') print(' </HarmonicAngleForce>')
...@@ -110,31 +113,22 @@ print(' </HarmonicAngleForce>') ...@@ -110,31 +113,22 @@ print(' </HarmonicAngleForce>')
print(' <PeriodicTorsionForce>') print(' <PeriodicTorsionForce>')
torsions = [f for f in forcefield._forces if isinstance(f, ff.PeriodicTorsionGenerator)][0] torsions = [f for f in forcefield._forces if isinstance(f, ff.PeriodicTorsionGenerator)][0]
wildcardType = forcefield._atomClasses['']
for torsion in torsions.proper: for torsion in torsions.proper:
type1 = next(iter(torsion.types1)) type = [next(iter(torsion.types1)), next(iter(torsion.types2)), next(iter(torsion.types3)), next(iter(torsion.types4))]
type2 = next(iter(torsion.types2)) wildcard = [torsion.types1 == wildcardType, torsion.types2 == wildcardType, torsion.types3 == wildcardType, torsion.types4 == wildcardType]
type3 = next(iter(torsion.types3)) if all(type[i] not in omitTypes or wildcard[i] for i in range(4)):
type4= next(iter(torsion.types4)) classes = tuple('' if wildcard[i] else forcefield._atomTypes[type[i]].atomClass for i in range(4))
if type1 not in omitTypes and type2 not in omitTypes and type3 not in omitTypes and type4 not in omitTypes: print(' <Proper class1="%s" class2="%s" class3="%s" class4="%s"' % classes, end=' ')
class1 = forcefield._atomTypes[type1][0]
class2 = forcefield._atomTypes[type2][0]
class3 = forcefield._atomTypes[type3][0]
class4 = forcefield._atomTypes[type4][0]
print(' <Proper class1="%s" class2="%s" class3="%s" class4="%s"' % (class1, class2, class3, class4), end=' ')
for i in range(len(torsion.k)): for i in range(len(torsion.k)):
print(' periodicity%d="%d" phase%d="%g" k%d="%g"' % (i+1, torsion.periodicity[i], i+1, torsion.phase[i], i+1, torsion.k[i]), end=' ') print(' periodicity%d="%d" phase%d="%g" k%d="%g"' % (i+1, torsion.periodicity[i], i+1, torsion.phase[i], i+1, torsion.k[i]), end=' ')
print('/>') print('/>')
for torsion in torsions.improper: for torsion in torsions.improper:
type1 = next(iter(torsion.types1)) type = [next(iter(torsion.types1)), next(iter(torsion.types2)), next(iter(torsion.types3)), next(iter(torsion.types4))]
type2 = next(iter(torsion.types2)) wildcard = [torsion.types1 == wildcardType, torsion.types2 == wildcardType, torsion.types3 == wildcardType, torsion.types4 == wildcardType]
type3 = next(iter(torsion.types3)) if all(type[i] not in omitTypes or wildcard[i] for i in range(4)):
type4= next(iter(torsion.types4)) classes = tuple('' if wildcard[i] else forcefield._atomTypes[type[i]].atomClass for i in range(4))
if type1 not in omitTypes and type2 not in omitTypes and type3 not in omitTypes and type4 not in omitTypes: print(' <Improper class1="%s" class2="%s" class3="%s" class4="%s"' % classes, end=' ')
class1 = forcefield._atomTypes[type1][0]
class2 = forcefield._atomTypes[type2][0]
class3 = forcefield._atomTypes[type3][0]
class4 = forcefield._atomTypes[type4][0]
print(' <Improper class1="%s" class2="%s" class3="%s" class4="%s"' % (class1, class2, class3, class4), end=' ')
for i in range(len(torsion.k)): for i in range(len(torsion.k)):
print(' periodicity%d="%d" phase%d="%g" k%d="%g"' % (i+1, torsion.periodicity[i], i+1, torsion.phase[i], i+1, torsion.k[i]), end=' ') print(' periodicity%d="%d" phase%d="%g" k%d="%g"' % (i+1, torsion.periodicity[i], i+1, torsion.phase[i], i+1, torsion.k[i]), end=' ')
print('/>') print('/>')
......
pdbfixer/soft.xml
This source diff could not be displayed because it is too large. You can view the blob instead.
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment