Commit 98652498 authored by peastman's avatar peastman
Browse files

Merge pull request #667 from peastman/gb

Optimizations to the custom GB models
parents 55bae85f 5bed3590
......@@ -43,7 +43,7 @@ import simtk.unit as u
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.internal.customgbforces import (GBSAHCTForce,
GBSAOBC1Force, GBSAOBC2Force, GBSAGBnForce, GBSAGBn2Force)
GBSAOBC1Force, GBSAOBC2Force, GBSAGBnForce, GBSAGBn2Force, convertParameters)
# CHARMM imports
from simtk.openmm.app.internal.charmm.topologyobjects import (
ResidueList, AtomList, TrackedList, Bond, Angle, Dihedral,
......@@ -1374,6 +1374,7 @@ class CharmmPsfFile(object):
if implicitSolvent is not None:
if verbose: print('Adding GB parameters...')
gb_parms = self._get_gb_params(implicitSolvent)
gb_parms = convertParameters(gb_parms, str(implicitSolvent))
# If implicitSolventKappa is None, compute it from salt
# concentration
......
......@@ -1032,6 +1032,8 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
if units.is_quantity(cutoff):
cutoff = cutoff.value_in_unit(units.nanometers)
gb_parms = prmtop.getGBParms(gbmodel, elements)
if gbmodel != 'OBC2' or implicitSolventKappa != 0:
gb_parms = customgb.convertParameters(gb_parms, gbmodel)
if gbmodel == 'HCT':
gb = customgb.GBSAHCTForce(solventDielectric, soluteDielectric, 'ACE', cutoff, implicitSolventKappa)
elif gbmodel == 'OBC1':
......
......@@ -189,36 +189,39 @@ for i in range (len(d0)):
m0[i]=m0[i]*10
def _createEnergyTerms(force, SA, cutoff, kappa):
def _createEnergyTerms(force, solventDielectric, soluteDielectric, SA, cutoff, kappa, offset):
# Add the energy terms to the CustomGBForce. These are identical for all the GB models.
params = "; solventDielectric=%.16g; soluteDielectric=%.16g; kappa=%.16g; offset=%.16g" % (solventDielectric, soluteDielectric, kappa, offset)
if cutoff is not None:
params += "; cutoff=%.16g" % cutoff
if kappa > 0:
force.addEnergyTerm("-0.5*138.935485*(1/soluteDielectric-exp(-kappa*B)/solventDielectric)*q^2/B",
force.addEnergyTerm("-0.5*138.935485*(1/soluteDielectric-exp(-kappa*B)/solventDielectric)*q^2/B"+params,
CustomGBForce.SingleParticle)
elif kappa < 0:
# Do kappa check here to avoid repeating code everywhere
raise ValueError('kappa/ionic strength must be >= 0')
else:
force.addEnergyTerm("-0.5*138.935485*(1/soluteDielectric-1/solventDielectric)*q^2/B",
force.addEnergyTerm("-0.5*138.935485*(1/soluteDielectric-1/solventDielectric)*q^2/B"+params,
CustomGBForce.SingleParticle)
if SA=='ACE':
force.addEnergyTerm("28.3919551*(radius+0.14)^2*(radius/B)^6", CustomGBForce.SingleParticle)
force.addEnergyTerm("28.3919551*(radius+0.14)^2*(radius/B)^6; radius=or+offset"+params, CustomGBForce.SingleParticle)
elif SA is not None:
raise ValueError('Unknown surface area method: '+SA)
if cutoff is None:
if kappa > 0:
force.addEnergyTerm("-138.935485*(1/soluteDielectric-exp(-kappa*f)/solventDielectric)*q1*q2/f;"
"f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))", CustomGBForce.ParticlePairNoExclusions)
"f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))"+params, CustomGBForce.ParticlePairNoExclusions)
else:
force.addEnergyTerm("-138.935485*(1/soluteDielectric-1/solventDielectric)*q1*q2/f;"
"f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))", CustomGBForce.ParticlePairNoExclusions)
"f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))"+params, CustomGBForce.ParticlePairNoExclusions)
else:
if kappa > 0:
force.addEnergyTerm("-138.935485*(1/soluteDielectric-kappa/solventDielectric)*q1*q2*(1/f-"+str(1/cutoff)+");"
"f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))", CustomGBForce.ParticlePairNoExclusions)
"f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))"+params, CustomGBForce.ParticlePairNoExclusions)
else:
force.addEnergyTerm("-138.935485*(1/soluteDielectric-1/solventDielectric)*q1*q2*(1/f-"+str(1/cutoff)+");"
"f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))", CustomGBForce.ParticlePairNoExclusions)
"f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))"+params, CustomGBForce.ParticlePairNoExclusions)
"""
Amber Equivalent: igb = 1
......@@ -231,22 +234,15 @@ def GBSAHCTForce(solventDielectric=78.5, soluteDielectric=1, SA=None,
custom = CustomGBForce()
custom.addPerParticleParameter("q")
custom.addPerParticleParameter("radius")
custom.addPerParticleParameter("scale")
if kappa > 0: custom.addGlobalParameter('kappa', kappa)
custom.addGlobalParameter("solventDielectric", solventDielectric)
custom.addGlobalParameter("soluteDielectric", soluteDielectric)
custom.addGlobalParameter("offset", 0.009)
custom.addPerParticleParameter("or") # Offset radius
custom.addPerParticleParameter("sr") # Scaled offset radius
custom.addComputedValue("I", "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r);"
"U=r+sr2;"
"L=max(or1, D);"
"D=abs(r-sr2);"
"sr2 = scale2*or2;"
"or1 = radius1-offset; or2 = radius2-offset", CustomGBForce.ParticlePairNoExclusions)
"D=abs(r-sr2)", CustomGBForce.ParticlePairNoExclusions)
custom.addComputedValue("B", "1/(1/or-I);"
"or=radius-offset", CustomGBForce.SingleParticle)
_createEnergyTerms(custom, SA, cutoff, kappa)
custom.addComputedValue("B", "1/(1/or-I)", CustomGBForce.SingleParticle)
_createEnergyTerms(custom, solventDielectric, soluteDielectric, SA, cutoff, kappa, 0.009)
return custom
"""
......@@ -258,22 +254,16 @@ def GBSAOBC1Force(solventDielectric=78.5, soluteDielectric=1, SA=None,
custom = CustomGBForce()
custom.addPerParticleParameter("q")
custom.addPerParticleParameter("radius")
custom.addPerParticleParameter("scale")
if kappa > 0: custom.addGlobalParameter('kappa', kappa)
custom.addGlobalParameter("solventDielectric", solventDielectric)
custom.addGlobalParameter("soluteDielectric", soluteDielectric)
custom.addGlobalParameter("offset", 0.009)
custom.addPerParticleParameter("or") # Offset radius
custom.addPerParticleParameter("sr") # Scaled offset radius
custom.addComputedValue("I", "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r);"
"U=r+sr2;"
"L=max(or1, D);"
"D=abs(r-sr2);"
"sr2 = scale2*or2;"
"or1 = radius1-offset; or2 = radius2-offset", CustomGBForce.ParticlePairNoExclusions)
"D=abs(r-sr2)", CustomGBForce.ParticlePairNoExclusions)
custom.addComputedValue("B", "1/(1/or-tanh(0.8*psi+2.909125*psi^3)/radius);"
"psi=I*or; or=radius-offset", CustomGBForce.SingleParticle)
_createEnergyTerms(custom, SA, cutoff, kappa)
"psi=I*or; radius=or+offset; offset=0.009", CustomGBForce.SingleParticle)
_createEnergyTerms(custom, solventDielectric, soluteDielectric, SA, cutoff, kappa, 0.009)
return custom
"""
......@@ -285,22 +275,16 @@ def GBSAOBC2Force(solventDielectric=78.5, soluteDielectric=1, SA=None,
custom = CustomGBForce()
custom.addPerParticleParameter("q")
custom.addPerParticleParameter("radius")
custom.addPerParticleParameter("scale")
if kappa > 0: custom.addGlobalParameter('kappa', kappa)
custom.addGlobalParameter("solventDielectric", solventDielectric)
custom.addGlobalParameter("soluteDielectric", soluteDielectric)
custom.addGlobalParameter("offset", 0.009)
custom.addPerParticleParameter("or") # Offset radius
custom.addPerParticleParameter("sr") # Scaled offset radius
custom.addComputedValue("I", "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r);"
"U=r+sr2;"
"L=max(or1, D);"
"D=abs(r-sr2);"
"sr2 = scale2*or2;"
"or1 = radius1-offset; or2 = radius2-offset", CustomGBForce.ParticlePairNoExclusions)
"D=abs(r-sr2)", CustomGBForce.ParticlePairNoExclusions)
custom.addComputedValue("B", "1/(1/or-tanh(psi-0.8*psi^2+4.85*psi^3)/radius);"
"psi=I*or; or=radius-offset", CustomGBForce.SingleParticle)
_createEnergyTerms(custom, SA, cutoff, kappa)
"psi=I*or; radius=or+offset; offset=0.009", CustomGBForce.SingleParticle)
_createEnergyTerms(custom, solventDielectric, soluteDielectric, SA, cutoff, kappa, 0.009)
return custom
"""
......@@ -321,15 +305,8 @@ def GBSAGBnForce(solventDielectric=78.5, soluteDielectric=1, SA=None,
custom = CustomGBForce()
custom.addPerParticleParameter("q")
custom.addPerParticleParameter("radius")
custom.addPerParticleParameter("scale")
if kappa > 0: custom.addGlobalParameter('kappa', kappa)
custom.addGlobalParameter("solventDielectric", solventDielectric)
custom.addGlobalParameter("soluteDielectric", soluteDielectric)
custom.addGlobalParameter("offset", 0.009)
custom.addGlobalParameter("neckScale", 0.361825)
custom.addGlobalParameter("neckCut", 0.68)
custom.addPerParticleParameter("or") # Offset radius
custom.addPerParticleParameter("sr") # Scaled offset radius
custom.addTabulatedFunction("getd0", Discrete1DFunction(d0))
custom.addTabulatedFunction("getm0", Discrete1DFunction(m0))
......@@ -341,12 +318,12 @@ def GBSAGBnForce(solventDielectric=78.5, soluteDielectric=1, SA=None,
"U=r+sr2;"
"L=max(or1, D);"
"D=abs(r-sr2);"
"sr2 = scale2*or2;"
"or1 = radius1-offset; or2 = radius2-offset", CustomGBForce.ParticlePairNoExclusions)
"radius1=or1+offset; radius2=or2+offset;"
"neckScale=0.361825; neckCut=0.68; offset=0.009", CustomGBForce.ParticlePairNoExclusions)
custom.addComputedValue("B", "1/(1/or-tanh(1.09511284*psi-1.907992938*psi^2+2.50798245*psi^3)/radius);"
"psi=I*or; or=radius-offset", CustomGBForce.SingleParticle)
_createEnergyTerms(custom, SA, cutoff, kappa)
"psi=I*or; radius=or+offset; offset=0.009", CustomGBForce.SingleParticle)
_createEnergyTerms(custom, solventDielectric, soluteDielectric, SA, cutoff, kappa, 0.009)
return custom
"""
......@@ -367,19 +344,12 @@ def GBSAGBn2Force(solventDielectric=78.5, soluteDielectric=1, SA=None,
custom = CustomGBForce()
custom.addPerParticleParameter("q")
custom.addPerParticleParameter("radius")
custom.addPerParticleParameter("scale")
custom.addPerParticleParameter("or") # Offset radius
custom.addPerParticleParameter("sr") # Scaled offset radius
custom.addPerParticleParameter("alpha")
custom.addPerParticleParameter("beta")
custom.addPerParticleParameter("gamma")
if kappa > 0: custom.addGlobalParameter('kappa', kappa)
custom.addGlobalParameter("solventDielectric", solventDielectric)
custom.addGlobalParameter("soluteDielectric", soluteDielectric)
custom.addGlobalParameter("offset", 0.0195141)
custom.addGlobalParameter("neckScale", 0.826836)
custom.addGlobalParameter("neckCut", 0.68)
custom.addTabulatedFunction("getd0", Discrete1DFunction(d0))
custom.addTabulatedFunction("getm0", Discrete1DFunction(m0))
......@@ -390,10 +360,24 @@ def GBSAGBn2Force(solventDielectric=78.5, soluteDielectric=1, SA=None,
"U=r+sr2;"
"L=max(or1, D);"
"D=abs(r-sr2);"
"sr2 = scale2*or2;"
"or1 = radius1-offset; or2 = radius2-offset", CustomGBForce.ParticlePairNoExclusions)
"radius1=or1+offset; radius2=or2+offset;"
"neckScale=0.826836; neckCut=0.68; offset=0.0195141", CustomGBForce.ParticlePairNoExclusions)
custom.addComputedValue("B", "1/(1/or-tanh(alpha*psi-beta*psi^2+gamma*psi^3)/radius);"
"psi=I*or; or=radius-offset", CustomGBForce.SingleParticle)
_createEnergyTerms(custom, SA, cutoff, kappa)
"psi=I*or; radius=or+offset; offset=0.0195141", CustomGBForce.SingleParticle)
_createEnergyTerms(custom, solventDielectric, soluteDielectric, SA, cutoff, kappa, 0.0195141)
return custom
def convertParameters(params, gbmodel):
"""Convert the GB parameters from the file into the values expected by the appropriate CustomGBForce."""
newparams = [None]*len(params)
if gbmodel == 'GBn2':
offset = 0.0195141
else:
offset = 0.009
for i in range(len(params)):
newparams[i] = list(params[i])
newparams[i][0] -= offset
newparams[i][1] *= newparams[i][0]
return newparams
\ No newline at end of file
......@@ -100,24 +100,13 @@ class TestAmberPrmtopFile(unittest.TestCase):
for method in methodMap:
system = prmtop2.createSystem(implicitSolvent=implicitSolvent_value,
solventDielectric=50.0, soluteDielectric=0.9, nonbondedMethod=method)
found_matching_solvent_dielectric=False
found_matching_solute_dielectric=False
if implicitSolvent_value in set([HCT, OBC1, GBn]):
for force in system.getForces():
if isinstance(force, CustomGBForce):
self.assertEqual(force.getNonbondedMethod(), methodMap[method])
for j in range(force.getNumGlobalParameters()):
if (force.getGlobalParameterName(j) == 'solventDielectric' and
force.getGlobalParameterDefaultValue(j) == 50.0):
found_matching_solvent_dielectric = True
if (force.getGlobalParameterName(j) == 'soluteDielectric' and
force.getGlobalParameterDefaultValue(j) == 0.9):
found_matching_solute_dielectric = True
if isinstance(force, NonbondedForce):
self.assertEqual(force.getReactionFieldDielectric(), 1.0)
self.assertEqual(force.getNonbondedMethod(), methodMap[method])
self.assertTrue(found_matching_solvent_dielectric and
found_matching_solute_dielectric)
else:
for force in system.getForces():
if isinstance(force, GBSAOBCForce):
......@@ -241,5 +230,24 @@ class TestAmberPrmtopFile(unittest.TestCase):
# Amber using this force field.
self.assertAlmostEqual(-7307.2735621/ene, 1, places=3)
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']
pdb = PDBFile('systems/alanine-dipeptide-implicit.pdb')
for i in range(5):
system = prmtop2.createSystem(implicitSolvent=solventType[i], nonbondedMethod=nonbondedMethod[i], implicitSolventSaltConc=salt[i])
integrator = VerletIntegrator(0.001)
context = Context(system, integrator, Platform.getPlatformByName("CPU"))
context.setPositions(pdb.positions)
state1 = context.getState(getForces=True)
state2 = XmlSerializer.deserialize(open('systems/alanine-dipeptide-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__':
unittest.main()
......@@ -69,21 +69,9 @@ class TestCharmmFiles(unittest.TestCase):
system = self.psf_x.createSystem(self.params, implicitSolvent=GBn,
solventDielectric=50.0,
soluteDielectric = 0.9)
found_matching_solvent_dielectric=False
found_matching_solute_dielectric=False
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):
self.assertEqual(force.getReactionFieldDielectric(), 1.0)
self.assertTrue(found_matching_solvent_dielectric and
found_matching_solute_dielectric)
def test_HydrogenMass(self):
"""Test that altering the mass of hydrogens works correctly."""
......@@ -127,6 +115,26 @@ class TestCharmmFiles(unittest.TestCase):
ene = state.getPotentialEnergy().value_in_unit(kilocalories_per_mole)
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__':
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>
<?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="125.03448115289211" y="43.32996869087219" z="-.32721251249313354"/>
<Force x="341.6234832406044" y="347.617433026433" z="-.8306254320777953"/>
<Force x="1.465131115168333" y="-51.00541362911463" z="-38.11954540386796"/>
<Force x="2.7838090173900127" y="-50.54201374202967" z="37.45753550902009"/>
<Force x="-694.3128572702408" y="-40.11703059077263" z="141.99350219080225"/>
<Force x="-507.0770571231842" y="-589.0909288525581" z="389.5440996333491"/>
<Force x="30.07020881283097" y="-865.4282982449513" z="-253.73720738082193"/>
<Force x="-46.64711940288544" y="-11.553307116031647" z="-1.14782332489267"/>
<Force x="412.9440284720622" y="393.1155642759986" z="63.61426969873719"/>
<Force x="-27.9816143065691" y="24.851514041423798" z="52.195745466742665"/>
<Force x="7.146615496603772" y="-242.70665084780194" z="-305.01612206571735"/>
<Force x="25.575431060045958" y="-15.982945490628481" z="-64.17020486621186"/>
<Force x="93.36584681086242" y="-178.10863384790719" z="-121.07110506272875"/>
<Force x="418.07980923354626" y="303.29553681379184" z="-342.0748033244163"/>
<Force x="214.27255564928055" y="-102.44735606014729" z="291.2024629737716"/>
<Force x="-420.52409839630127" y="128.17991256713867" z="2.5522894316818565"/>
<Force x="-241.46922235935926" y="416.7907632468268" z="152.65439392975532"/>
<Force x="-52.96275818347931" y="-151.0294795036316" z="-1.817022401606664"/>
<Force x="246.45249317586422" y="187.45216877479106" z="-1.6977270720526576"/>
<Force x="76.24567371606827" y="320.694318652153" z="-.5214127898216248"/>
<Force x="-1.9259837809950113" y="67.05363237857819" z="8.533471319824457"/>
<Force x="-2.1587012950330973" y="65.63135588169098" z="-9.216962072998285"/>
</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="116.92173205316067" y="35.74797606468201" z="-.31320083141326904"/>
<Force x="346.4335723519325" y="347.91397757828236" z="-.7274334602989256"/>
<Force x="-3.8193766362965107" y="-50.163410760462284" z="-34.13527924194932"/>
<Force x="-2.5832430608570576" y="-49.72161541134119" z="33.492162469774485"/>
<Force x="-678.0981813669205" y="-43.59890314936638" z="144.39658980676904"/>
<Force x="-506.10328698158264" y="-589.8690050244331" z="388.1552735015284"/>
<Force x="22.236811901209876" y="-862.8642609219532" z="-255.22323748259805"/>
<Force x="-45.493120074272156" y="-14.552918016910553" z="-1.145501604769379"/>
<Force x="411.3508411399089" y="390.2969187987037" z="66.42137118685059"/>
<Force x="-25.265792086720467" y="27.495839297771454" z="54.14952242234722"/>
<Force x="3.733357897726819" y="-243.36420928896405" z="-297.93731590383686"/>
<Force x="25.430794905871153" y="-14.971891451627016" z="-70.3868496301584"/>
<Force x="97.456183610484" y="-176.72447475977242" z="-122.43685823655687"/>
<Force x="422.9076767414808" y="300.37319069867954" z="-346.16921898908913"/>
<Force x="203.8950589299202" y="-90.0563190728426" z="292.4585710193496"/>
<Force x="-417.8394365310669" y="132.75754928588867" z="2.3947719985153526"/>
<Force x="-232.2461846396327" y="410.49045501928777" z="152.14686514786445"/>
<Force x="-57.324524998664856" y="-151.90891790390015" z="-2.435667592799291"/>
<Force x="241.87253235280514" y="184.33851583022624" z="-1.5684559224173427"/>
<Force x="75.61852604150772" y="324.36484706401825" z="-.5052977800369263"/>
<Force x=".5770832356065512" y="67.66817677021027" z="11.481010410934687"/>
<Force x=".3389774616807699" y="66.34845554828644" z="-12.111703846603632"/>
</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="119.87015162408352" y="35.32908082008362" z="-.2700755298137665"/>
<Force x="317.75814777612686" y="347.0047370940447" z="-.36405922705307603"/>
<Force x="-1.5870520360767841" y="-49.64923916012049" z="-32.16182732954621"/>
<Force x="-.6230551488697529" y="-49.28765354305506" z="31.62008560076356"/>
<Force x="-662.1668611764908" y="-41.318858593702316" z="142.827594331"/>
<Force x="-499.81304955482483" y="-628.789857685566" z="396.9687448188197"/>
<Force x="9.883120800135657" y="-860.5162613491993" z="-256.9060352721717"/>
<Force x="-47.38457667827606" y="-20.238594591617584" z="-1.099657526705414"/>
<Force x="425.8578887931071" y="385.30582945467904" z="62.49412080156617"/>
<Force x="-27.185586169362068" y="24.652142703533173" z="56.237679122481495"/>
<Force x="6.795258990256116" y="-237.71378814638592" z="-303.22915097349323"/>
<Force x="23.833146285265684" y="-15.633165407925844" z="-65.40645025996491"/>
<Force x="94.45341509394348" y="-180.38257107324898" z="-123.42125898576342"/>
<Force x="425.9342289417982" y="308.1854618168436" z="-349.6926793772727"/>
<Force x="189.02624255418777" y="-87.48103220760822" z="291.035372557817"/>
<Force x="-438.08852100372314" y="149.00048446655273" z="1.9025959426071495"/>
<Force x="-214.21626215428114" y="412.8787807105109" z="150.88754918030463"/>
<Force x="-49.139330983161926" y="-146.01251077651978" z=".8237868959549814"/>
<Force x="250.3118935674429" y="196.49128634948283" z="-1.326529967598617"/>
<Force x="79.39004474878311" y="327.3214110136032" z="-.41108113527297974"/>
<Force x="-1.35963093675673" y="65.87895834445953" z="11.501466009765863"/>
<Force x="-1.5495685283094645" y="64.97552835941315" z="-12.010117743164301"/>
</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="118.75197611749172" y="37.31967115402222" z="-.27059119939804077"/>
<Force x="319.6333022713661" y="341.64345021545887" z="-.5288973026908934"/>
<Force x=".2122395746409893" y="-48.4974714294076" z="-33.63206410780549"/>
<Force x="1.172732662409544" y="-48.143215753138065" z="33.092639449983835"/>
<Force x="-684.7492891550064" y="-22.71616694331169" z="142.82319979975"/>
<Force x="-480.4688265323639" y="-646.0529344677925" z="396.2935791655909"/>
<Force x="19.13356330501847" y="-860.4653923611622" z="-259.51502784877084"/>
<Force x="-55.92518413066864" y="-27.145088732242584" z="-.7603530795313418"/>
<Force x="429.2403093329631" y="382.1770090353675" z="60.34707065927796"/>
<Force x="-26.955144122242928" y="24.935925662517548" z="56.52590095857158"/>
<Force x="3.497140398947522" y="-239.83717452944256" z="-301.7230455505196"/>
<Force x="24.14258260652423" y="-15.982009936124086" z="-62.18898533610627"/>
<Force x="92.46831530146301" y="-180.48632702417672" z="-122.92353540635668"/>
<Force x="424.0495262593031" y="312.01640206342563" z="-348.84565827436745"/>
<Force x="204.3266766667366" y="-102.10523875057697" z="290.7208612586837"/>
<Force x="-454.11176013946533" y="168.45798873901367" z="1.5088972502853721"/>
<Force x="-223.03344660252333" y="420.29524646978825" z="149.7152159397956"/>
<Force x="-39.53340256214142" y="-150.4378457069397" z="1.7263585741166025"/>
<Force x="248.08668254315853" y="195.30956656951457" z="-1.4342415807768703"/>
<Force x="83.9709010720253" y="328.7274760007858" z="-.4131510257720947"/>
<Force x="-1.8617023173719645" y="65.95496952533722" z="9.68197262659669"/>
<Force x="-2.0471442881971598" y="65.03140318393707" z="-10.200121019035578"/>
</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="118.3535195440054" y="35.969449162483215" z="-.281414270401001"/>
<Force x="314.4263282418251" y="344.44179029762745" z="-.434397476259619"/>
<Force x="-1.4066988714039326" y="-49.76428470760584" z="-32.37674808874726"/>
<Force x="-.40957896783947945" y="-49.397793389856815" z="31.823710206896067"/>
<Force x="-657.3909059762955" y="-38.007091015577316" z="141.8415627987124"/>
<Force x="-493.7057650089264" y="-643.0512407422066" z="398.3075271293055"/>
<Force x="7.296756054041907" y="-855.5627281765919" z="-257.9519661345985"/>
<Force x="-53.50385844707489" y="-24.334915697574615" z="-1.1155095011927187"/>
<Force x="429.9477740279399" y="379.88140909792855" z="60.92471059667878"/>
<Force x="-26.743043139576912" y="21.195145785808563" z="59.63117956975475"/>
<Force x="5.737775317160413" y="-236.0452548211906" z="-302.6426021682564"/>
<Force x="23.517739962786436" y="-16.610735941678286" z="-63.854551303666085"/>
<Force x="94.2509424071759" y="-182.99365315027535" z="-125.87754732347094"/>
<Force x="424.38538555800915" y="311.03262597089633" z="-349.64496991224587"/>
<Force x="177.18182879686356" y="-78.61917002499104" z="287.95103532331996"/>
<Force x="-442.6554002761841" y="158.27599716186523" z="2.4879631453659385"/>
<Force x="-211.27339106053114" y="403.45306001883" z="150.30480124405585"/>
<Force x="-37.10933029651642" y="-140.19635248184204" z="3.2473717445973307"/>
<Force x="249.92661701142788" y="200.21128232497722" z="-1.3900556145235896"/>
<Force x="83.43381077051163" y="328.46902453899384" z="-.4314945936203003"/>
<Force x="-2.0295280162245035" y="66.29071056842804" z="10.737368497997522"/>
<Force x="-2.2310355845838785" y="65.36275541782379" z="-11.255987975746393"/>
</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