Commit 5bed3590 authored by Peter Eastman's avatar Peter Eastman
Browse files

Fixed CHARMM parser to work with changes to custom GB models

parent 71b67421
...@@ -43,7 +43,7 @@ import simtk.unit as u ...@@ -43,7 +43,7 @@ import simtk.unit as u
from simtk.openmm.app import (forcefield as ff, Topology, element) from simtk.openmm.app import (forcefield as ff, Topology, element)
from simtk.openmm.app.amberprmtopfile import HCT, OBC1, OBC2, GBn, GBn2 from simtk.openmm.app.amberprmtopfile import HCT, OBC1, OBC2, GBn, GBn2
from simtk.openmm.app.internal.customgbforces import (GBSAHCTForce, from simtk.openmm.app.internal.customgbforces import (GBSAHCTForce,
GBSAOBC1Force, GBSAOBC2Force, GBSAGBnForce, GBSAGBn2Force) GBSAOBC1Force, GBSAOBC2Force, GBSAGBnForce, GBSAGBn2Force, convertParameters)
# CHARMM imports # CHARMM imports
from simtk.openmm.app.internal.charmm.topologyobjects import ( from simtk.openmm.app.internal.charmm.topologyobjects import (
ResidueList, AtomList, TrackedList, Bond, Angle, Dihedral, ResidueList, AtomList, TrackedList, Bond, Angle, Dihedral,
...@@ -1374,6 +1374,7 @@ class CharmmPsfFile(object): ...@@ -1374,6 +1374,7 @@ class CharmmPsfFile(object):
if implicitSolvent is not None: if implicitSolvent is not None:
if verbose: print('Adding GB parameters...') if verbose: print('Adding GB parameters...')
gb_parms = self._get_gb_params(implicitSolvent) gb_parms = self._get_gb_params(implicitSolvent)
gb_parms = convertParameters(gb_parms, str(implicitSolvent))
# If implicitSolventKappa is None, compute it from salt # If implicitSolventKappa is None, compute it from salt
# concentration # concentration
......
...@@ -69,21 +69,9 @@ class TestCharmmFiles(unittest.TestCase): ...@@ -69,21 +69,9 @@ class TestCharmmFiles(unittest.TestCase):
system = self.psf_x.createSystem(self.params, implicitSolvent=GBn, system = self.psf_x.createSystem(self.params, implicitSolvent=GBn,
solventDielectric=50.0, solventDielectric=50.0,
soluteDielectric = 0.9) soluteDielectric = 0.9)
found_matching_solvent_dielectric=False
found_matching_solute_dielectric=False
for force in system.getForces(): for force in system.getForces():
if isinstance(force, CustomGBForce):
for i in range(force.getNumGlobalParameters()):
if force.getGlobalParameterName(i) == 'solventDielectric':
if force.getGlobalParameterDefaultValue(i) == 50.0:
found_matching_solvent_dielectric = True
elif force.getGlobalParameterName(i) == 'soluteDielectric':
if force.getGlobalParameterDefaultValue(i) == 0.9:
found_matching_solute_dielectric = True
if isinstance(force, NonbondedForce): if isinstance(force, NonbondedForce):
self.assertEqual(force.getReactionFieldDielectric(), 1.0) self.assertEqual(force.getReactionFieldDielectric(), 1.0)
self.assertTrue(found_matching_solvent_dielectric and
found_matching_solute_dielectric)
def test_HydrogenMass(self): def test_HydrogenMass(self):
"""Test that altering the mass of hydrogens works correctly.""" """Test that altering the mass of hydrogens works correctly."""
...@@ -127,6 +115,26 @@ class TestCharmmFiles(unittest.TestCase): ...@@ -127,6 +115,26 @@ class TestCharmmFiles(unittest.TestCase):
ene = state.getPotentialEnergy().value_in_unit(kilocalories_per_mole) ene = state.getPotentialEnergy().value_in_unit(kilocalories_per_mole)
self.assertAlmostEqual(ene, 15490.0033559, delta=0.05) self.assertAlmostEqual(ene, 15490.0033559, delta=0.05)
def test_ImplicitSolventForces(self):
"""Compute forces for different implicit solvent types, and compare them to ones generated with a previous version of OpenMM to ensure they haven't changed."""
solventType = [HCT, OBC1, OBC2, GBn, GBn2]
nonbondedMethod = [NoCutoff, CutoffNonPeriodic, CutoffNonPeriodic, NoCutoff, NoCutoff]
salt = [0.0, 0.0, 0.5, 0.5, 0.0]*(moles/liter)
file = ['HCT_NoCutoff', 'OBC1_NonPeriodic', 'OBC2_NonPeriodic_Salt', 'GBn_NoCutoff_Salt', 'GBn2_NoCutoff']
for i in range(5):
system = self.psf_c.createSystem(self.params, implicitSolvent=solventType[i], nonbondedMethod=nonbondedMethod[i], implicitSolventSaltConc=salt[i])
integrator = VerletIntegrator(0.001)
context = Context(system, integrator, Platform.getPlatformByName("CPU"))
context.setPositions(self.pdb.positions)
state1 = context.getState(getForces=True)
#out = open('systems/ala-ala-ala-implicit-forces/'+file[i]+'.xml', 'w')
#out.write(XmlSerializer.serialize(state1))
#out.close()
state2 = XmlSerializer.deserialize(open('systems/ala-ala-ala-implicit-forces/'+file[i]+'.xml').read())
for f1, f2, in zip(state1.getForces().value_in_unit(kilojoules_per_mole/nanometer), state2.getForces().value_in_unit(kilojoules_per_mole/nanometer)):
diff = norm(f1-f2)
self.assertTrue(diff < 0.1 or diff/norm(f1) < 1e-4)
if __name__ == '__main__': if __name__ == '__main__':
unittest.main() unittest.main()
<?xml version="1.0" ?>
<State time="0" type="State" version="1">
<PeriodicBoxVectors>
<A x="2" y="0" z="0"/>
<B x="0" y="2" z="0"/>
<C x="0" y="0" z="2"/>
</PeriodicBoxVectors>
<Forces>
<Force x="183.12361104909448" y="155.26784211796442" z="34.691159962151325"/>
<Force x="58.02808994091839" y="161.1479771333285" z="-6.905719184803445"/>
<Force x="53.116947200286404" y="-197.39924721297302" z="-191.13127107770634"/>
<Force x="46.218727097648014" y="-212.66428112142592" z="171.32495750170045"/>
<Force x="-143.51220792224802" y="-77.6749878542526" z="106.54286260678916"/>
<Force x="-122.2580981087889" y="-56.69313081273846" z="155.02475224419112"/>
<Force x="256.861210924291" y="4.991885101995734" z="56.89967589615151"/>
<Force x="-50.134423267799164" y="26.69582029770774" z="-148.0309316306637"/>
<Force x="93.82349928954996" y="-98.75344995978912" z="-30.84769608742689"/>
<Force x="-41.91077517758289" y="59.289304779748626" z="-110.90263298668452"/>
<Force x="-402.5528947784169" y="-502.5128785408756" z="146.8525667854756"/>
<Force x="140.35927800593402" y="706.2615945572692" z="-116.71709633774552"/>
<Force x="336.6168876754595" y="357.247339062576" z="-34.749492380483105"/>
<Force x="-415.7108187560625" y="-212.13687180023163" z="100.590812009496"/>
<Force x="-2064.878698363039" y="-1169.1702011127431" z="-141.28708209576615"/>
<Force x="-194.9400953435316" y="-24.670932939610708" z="11.635635952592624"/>
<Force x="-733.9608254321217" y="-237.08948410189265" z="-552.6426811401486"/>
<Force x="-1060.8507100087338" y="-17.22971806538561" z="-834.9535646446606"/>
<Force x="-11.966996404264762" y="-67.77194077117352" z="-66.3175636555273"/>
<Force x="-5.938422468557391" y="-.8846999866901797" z="-79.25483041064173"/>
<Force x="-535049.450575994" y="-759418.9456168906" z="-62586.54371173267"/>
<Force x="-811.1877273151022" y="-752.1488295390167" z="-559.2811902937815"/>
<Force x="-59122.721636109374" y="-11224.366289486006" z="4161.654215282473"/>
<Force x="8491.852921433589" y="-4987.486057146457" z="198.97782166846662"/>
<Force x="-92073.12391443842" y="-81863.60475973453" z="-3121.980801468735"/>
<Force x="629672.3177092219" y="832975.8356195987" z="70299.92572789358"/>
<Force x="137697.00593367853" y="-29842.060582353926" z="13288.173983467228"/>
<Force x="-16855.75815442956" y="-20321.644238398876" z="-2810.1468347196355"/>
<Force x="-35343.17474223994" y="19898.531265412435" z="9385.986414694597"/>
<Force x="-6566690.1522480175" y="1625423.158741199" z="-3339704.996853593"/>
<Force x="8959691.989050966" y="469395.856032878" z="202260.29458456024"/>
<Force x="-1344238.7618171684" y="1727648.777500562" z="3128189.6150764455"/>
<Force x="-1081288.5392807235" y="-3765528.083662639" z="-17472.013005055705"/>
</Forces>
</State>
<?xml version="1.0" ?>
<State time="0" type="State" version="1">
<PeriodicBoxVectors>
<A x="2" y="0" z="0"/>
<B x="0" y="2" z="0"/>
<C x="0" y="0" z="2"/>
</PeriodicBoxVectors>
<Forces>
<Force x="181.30306966725854" y="159.29243928593317" z="35.637069462273395"/>
<Force x="66.09966892041058" y="128.05425765090664" z="-9.903758430408914"/>
<Force x="10.526981379973904" y="-196.22108559187927" z="-167.87643404157353"/>
<Force x="44.156288132804264" y="-215.63223766439467" z="132.1042391179114"/>
<Force x="-165.5816354124824" y="-73.37206274438932" z="116.79308385923056"/>
<Force x="-93.0862116646971" y="-66.09492906102948" z="148.86181889458175"/>
<Force x="244.00171385397852" y="8.38662997504261" z="48.56955962417885"/>
<Force x="-47.67093182797983" y="20.03437181900901" z="-143.35024849467737"/>
<Force x="97.01566847900308" y="-101.16990427497467" z="-22.730537197656382"/>
<Force x="-28.89086245785633" y="50.526134537439056" z="-99.99265984215327"/>
<Force x="-425.5357133819325" y="-510.30285656821934" z="146.09318263020216"/>
<Force x="168.28850774226214" y="748.6720865006286" z="-113.89239884323868"/>
<Force x="329.44051133268607" y="372.46074467170683" z="-40.64059803233857"/>
<Force x="-421.84491299434376" y="-220.10242126922577" z="99.85900811545304"/>
<Force x="-2079.694860472414" y="-1170.8042343158681" z="-140.53129901471146"/>
<Force x="-197.0285963200941" y="-18.209816148839224" z="10.138367275834812"/>
<Force x="-737.922312248528" y="-243.0193547073614" z="-549.1476860229611"/>
<Force x="-1058.3116475087338" y="-16.151417589311393" z="-834.6196413048168"/>
<Force x="-9.61508676681359" y="-66.0211377773991" z="-67.43522417798823"/>
<Force x="-5.07464530791286" y="5.05774386096607" z="-82.80545144335657"/>
<Force x="-535008.3027793631" y="-759464.7106467978" z="-62553.71453387623"/>
<Force x="-862.5694411822897" y="-750.0902968241729" z="-623.4501966414377"/>
<Force x="-59105.97902517798" y="-11220.1010917321" z="4132.6760353508325"/>
<Force x="8500.659516282221" y="-4979.370342119601" z="193.63300355323224"/>
<Force x="-92088.56655855768" y="-81859.31729559269" z="-3117.676529770737"/>
<Force x="629680.3119413996" y="832994.1522089542" z="70286.58493046927"/>
<Force x="137692.58820906916" y="-29845.6824085258" z="13291.340731514103"/>
<Force x="-16865.853125132686" y="-20307.415448115673" z="-2800.674971926667"/>
<Force x="-35340.09261943721" y="19914.31762710677" z="9365.387474798967"/>
<Force x="-6566805.4022480175" y="1625380.033741199" z="-3339664.496853593"/>
<Force x="8959673.99140082" y="469389.32107117755" z="202279.90824727996"/>
<Force x="-1344166.3868171684" y="1727692.902500562" z="3128204.6150764455"/>
<Force x="-1081175.4142807235" y="-3765539.333662639" z="-17459.43121794633"/>
</Forces>
</State>
<?xml version="1.0" ?>
<State time="0" type="State" version="1">
<PeriodicBoxVectors>
<A x="2" y="0" z="0"/>
<B x="0" y="2" z="0"/>
<C x="0" y="0" z="2"/>
</PeriodicBoxVectors>
<Forces>
<Force x="222.31410177175073" y="158.23418122929255" z="62.23766207549361"/>
<Force x="31.58283084668011" y="56.528867025906635" z="-49.24259510033079"/>
<Force x="-68.39497754336594" y="-183.3738489585785" z="-144.30000276715947"/>
<Force x="16.829643235343326" y="-185.04890941732435" z="53.43717063646609"/>
<Force x="-78.42200940540233" y="-22.341671051274083" z="164.26952428891806"/>
<Force x="-92.5467943978514" y="-56.300773940179866" z="163.112368210988"/>
<Force x="242.50880919089258" y="18.500357544622688" z="43.11554579971596"/>
<Force x="-38.76807214017465" y="14.818074669350807" z="-144.67003838115198"/>
<Force x="96.32516059020914" y="-104.26505808356842" z="-14.603355190820444"/>
<Force x="-24.334008209809454" y="38.8018417823121" z="-96.48091820396968"/>
<Force x="-411.4825899078114" y="-455.03906262290684" z="140.87504298176466"/>
<Force x="154.93731145319964" y="685.3672464127379" z="-110.25670815415177"/>
<Force x="308.8163200485064" y="345.19005756843535" z="-30.099754068715527"/>
<Force x="-403.72101162715626" y="-212.12824676971405" z="92.38245324484757"/>
<Force x="-2072.498571409914" y="-1146.2246444721181" z="-137.43201007428178"/>
<Force x="-192.43813458913706" y="-22.078983476720083" z="15.487461666703952"/>
<Force x="-734.999216545403" y="-239.14217270052546" z="-553.3522148315549"/>
<Force x="-1061.3953877431088" y="-17.592316484575065" z="-834.3818483360668"/>
<Force x="-13.522195073698356" y="-70.5011468563786" z="-65.48352968946284"/>
<Force x="-6.612387922658954" y=".5695389049113828" z="-80.99231820361048"/>
<Force x="-535051.2465965018" y="-759438.2090675131" z="-62597.107310365485"/>
<Force x="-884.2599929401022" y="-797.3043471171417" z="-632.1942762312815"/>
<Force x="-59132.63463049414" y="-11137.957750667647" z="4081.7088112297383"/>
<Force x="8528.217522568842" y="-5016.619579180148" z="204.09790452369123"/>
<Force x="-92081.11046381586" y="-81866.61880163517" z="-3128.77401435936"/>
<Force x="629647.9728300715" y="833005.4992853702" z="70272.57769780325"/>
<Force x="137667.60212508478" y="-29817.4577991508" z="13267.939608467228"/>
<Force x="-16849.90088514245" y="-20264.64237682661" z="-2859.725020754792"/>
<Force x="-35329.194647330274" y="19910.60506073104" z="9347.950028204728"/>
<Force x="-6566760.6522480175" y="1625379.283741199" z="-3339633.246853593"/>
<Force x="8959741.435523134" y="469410.801597148" z="202249.88889913543"/>
<Force x="-1344212.3868171684" y="1727628.777500562" z="3128339.3650764455"/>
<Force x="-1081157.9142807235" y="-3765599.583662639" z="-17386.979313649455"/>
</Forces>
</State>
<?xml version="1.0" ?>
<State time="0" type="State" version="1">
<PeriodicBoxVectors>
<A x="2" y="0" z="0"/>
<B x="0" y="2" z="0"/>
<C x="0" y="0" z="2"/>
</PeriodicBoxVectors>
<Forces>
<Force x="223.10238607350854" y="161.07696229618708" z="53.57150912139205"/>
<Force x="23.494802892578548" y="41.735654135281635" z="-52.5400346755261"/>
<Force x="-92.55845257998703" y="-192.42526344832459" z="-143.55398531110478"/>
<Force x="-.5067253254476896" y="-197.50359081380873" z="36.47189964037234"/>
<Force x="-71.12311166217724" y="-22.923248352299474" z="176.16349859311728"/>
<Force x="-96.18238399746078" y="-49.78676637182049" z="182.7166696636247"/>
<Force x="259.34113035788477" y="3.8626164557555" z="52.02276091812416"/>
<Force x="-36.2982378124769" y="14.8250612873928" z="-145.59947411113245"/>
<Force x="96.56693991760172" y="-105.06556849006256" z="-10.534827969140757"/>
<Force x="-28.052346222504767" y="34.63957409737558" z="-98.39121559776851"/>
<Force x="-407.87485522763563" y="-474.4734498299381" z="148.8828402229756"/>
<Force x="147.76729802546527" y="736.635572218402" z="-121.47636223740372"/>
<Force x="316.56870041960013" y="358.55342655647246" z="-36.919048044545605"/>
<Force x="-409.68670986934376" y="-228.9934544418332" z="90.34555749289444"/>
<Force x="-2068.203893675539" y="-1136.9321640033681" z="-138.67511310162553"/>
<Force x="-190.68837110036753" y="-21.274162462315786" z="20.88912067854477"/>
<Force x="-735.502878654778" y="-239.51714218294734" z="-553.5976371948361"/>
<Force x="-1062.0531025868588" y="-15.81420979512194" z="-835.9012575157543"/>
<Force x="-12.003510686491325" y="-71.06613022429852" z="-68.69160998121089"/>
<Force x="-3.519080427541766" y="3.9758651988566953" z="-83.97223763720423"/>
<Force x="-535082.8585349784" y="-759414.3580314413" z="-62615.07069308643"/>
<Force x="-836.7894118854147" y="-806.9301405741729" z="-616.891663926594"/>
<Force x="-59128.21766119482" y="-11171.795824398116" z="4101.884592479739"/>
<Force x="8512.672028489253" y="-4991.504176958469" z="198.77263673194318"/>
<Force x="-92080.33764429681" y="-81868.46521429386" z="-3130.532299881821"/>
<Force x="629659.817190423" y="832994.4047113956" z="70281.15389256644"/>
<Force x="137692.96052352228" y="-29826.63016243205" z="13279.40994538129"/>
<Force x="-16853.26160291589" y="-20269.398999141064" z="-2820.911116946198"/>
<Force x="-35330.140265006055" y="19900.15206161605" z="9353.563283705522"/>
<Force x="-6566768.9022480175" y="1625391.033741199" z="-3339651.496853593"/>
<Force x="8959779.915076356" y="469416.5334315596" z="202232.03704671844"/>
<Force x="-1344235.2618171684" y="1727582.902500562" z="3128285.3650764455"/>
<Force x="-1081183.2892807235" y="-3765535.833662639" z="-17369.176579274455"/>
</Forces>
</State>
<?xml version="1.0" ?>
<State time="0" type="State" version="1">
<PeriodicBoxVectors>
<A x="2" y="0" z="0"/>
<B x="0" y="2" z="0"/>
<C x="0" y="0" z="2"/>
</PeriodicBoxVectors>
<Forces>
<Force x="230.53816183034448" y="161.56808931034723" z="57.52342867705611"/>
<Force x="12.564077794922298" y="-.47807572311680246" z="-67.73089160911985"/>
<Force x="-138.11146161319016" y="-187.28762917098084" z="-128.55946321637822"/>
<Force x="-14.216667189461361" y="-184.42361950033217" z="-7.755578386971408"/>
<Force x="-51.715805944708976" y="-8.607643188725255" z="189.75546331479697"/>
<Force x="-92.82982967617171" y="-37.94701798925213" z="201.6270395366716"/>
<Force x="254.05585203757227" y="10.254599487982063" z="45.8494172596769"/>
<Force x="-27.659651768165375" y="8.348966087807838" z="-142.57089821391565"/>
<Force x="101.01881369689859" y="-109.75517229560455" z="-.21244790322278817"/>
<Force x="-21.38037204037586" y="22.86551289437265" z="-92.96600655235835"/>
<Force x="-397.09125537412" y="-424.17272961509434" z="135.01295191731154"/>
<Force x="156.04946355280902" y="705.0998513932067" z="-112.9041481108168"/>
<Force x="306.5652214156939" y="346.46968632209746" z="-27.505218241078808"/>
<Force x="-391.75561856075" y="-226.95233582000702" z="83.27667931906632"/>
<Force x="-2067.734655394289" y="-1129.4901474018056" z="-140.06254901471146"/>
<Force x="-191.13938275808238" y="-18.663742235265005" z="22.03099880659653"/>
<Force x="-729.4196877368092" y="-242.06033981478328" z="-550.7472953979611"/>
<Force x="-1061.6405049306088" y="-16.047501421098502" z="-835.6902589805981"/>
<Force x="-12.191029759977653" y="-71.33537155730633" z="-67.96942675367183"/>
<Force x="-3.4285805496120787" y="3.4221542613566953" z="-83.54534249560267"/>
<Force x="-535059.354506658" y="-759427.0708915488" z="-62607.30990893726"/>
<Force x="-865.6928542682272" y="-791.6607924296417" z="-636.7355970320627"/>
<Force x="-59134.540453248046" y="-11164.06852947624" z="4084.1721748771993"/>
<Force x="8523.705742783686" y="-5007.283077104953" z="201.13411363013654"/>
<Force x="-92085.95982714288" y="-81874.13444189396" z="-3131.0444535073093"/>
<Force x="629671.6843168879" y="832997.032213837" z="70286.08367161917"/>
<Force x="137691.52277938166" y="-29833.1218616508" z="13271.79861725629"/>
<Force x="-16859.409597911006" y="-20265.883709834423" z="-2793.2637780790105"/>
<Force x="-35331.66387798213" y="19901.063774262533" z="9356.077337288667"/>
<Force x="-6566795.1522480175" y="1625389.908741199" z="-3339649.746853593"/>
<Force x="8959758.524878602" y="469416.7061763106" z="202234.18939046844"/>
<Force x="-1344222.3868171684" y="1727597.777500562" z="3128292.1150764455"/>
<Force x="-1081153.1642807235" y="-3765540.083662639" z="-17383.96051482133"/>
</Forces>
</State>
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