from __future__ import print_function
from simtk.openmm import app
import simtk.openmm as mm
from simtk import unit
import sys
import mbpol
pdb = app.PDBFile("water14_cluster.pdb")
forcefield = app.ForceField("mbpol.xml")
# use tip4p
#forcefield = app.ForceField("tip4pfb.xml")
system = forcefield.createSystem(pdb.topology, nonbondedMethod=app.CutoffNonPeriodic, nonBondedCutoff=1e3*unit.nanometer)
integrator = mm.VerletIntegrator(0.00001*unit.femtoseconds)
platform = mm.Platform.getPlatformByName('Reference')
simulation = app.Simulation(pdb.topology, system, integrator, platform)
simulation.context.setPositions(pdb.positions)
simulation.context.computeVirtualSites()
state = simulation.context.getState(getForces=True, getEnergy=True)
potential_energy = state.getPotentialEnergy()
potential_energy.in_units_of(unit.kilocalorie_per_mole)
Quantity(value=-59.08694620286226, unit=kilocalorie/mole)
kilocalorie_per_mole_per_angstrom = unit.kilocalorie_per_mole/unit.angstrom
for f in state.getForces():
print(f.in_units_of(kilocalorie_per_mole_per_angstrom))
(9.361155484266149, 24.6798696179452, 52.85333122003323) kcal/(A mol) (12.166966073937992, -24.355536441527022, -25.340594147661065) kcal/(A mol) (-16.97652467171154, 3.556694920697774, -26.503543851222314) kcal/(A mol) (-40.10091402032537, 21.61077279055185, -9.649226045382184) kcal/(A mol) (36.48423956423134, 0.8625984852644293, -6.918559745188426) kcal/(A mol) (-49.88129048702903, -1.7979500802603612, 2.096316985531012) kcal/(A mol) (11.208648607179654, 2.0062968057561514, 2.944334483995184) kcal/(A mol) (-28.97149876118144, -5.65740371993149, 35.957040708878466) kcal/(A mol) (-23.378226317254697, 19.827807887879363, 26.3007621639795) kcal/(A mol) (33.1974506708894, -25.060760649092913, -25.03692529296061) kcal/(A mol) (-6.1011578636490125, 4.197580658974067, -6.6099988209806675) kcal/(A mol) (15.617213990925427, -13.783295711575882, -21.70331351341204) kcal/(A mol) (-25.52400183375194, 33.84910786390996, -25.363089216435615) kcal/(A mol) (13.641684734409363, -12.580904291747116, 25.61972380198023) kcal/(A mol) (13.034745650156811, -17.154178238712454, 1.8507497091725407) kcal/(A mol) (-1.3570052883969652, -7.99447684491713, -16.099990726342906) kcal/(A mol) (-32.663102708767184, 64.78082728023968, -34.098537666111355) kcal/(A mol) (40.24412080370862, -49.08121503432236, 32.37658075755621) kcal/(A mol) (-8.480283905570502, -21.046295803842625, 2.0721861116701703) kcal/(A mol) (18.48158680408055, 0.455445163057897, 17.56514529096203) kcal/(A mol) (13.700741104389452, -7.7334176880584575, 21.159216910152768) kcal/(A mol) (-11.804968670715253, 4.7011442493489675, -11.481564130240853) kcal/(A mol) (-2.3886415575171354, -2.3873744728366675, -3.0947339750170646) kcal/(A mol) (-25.168794765268352, 12.877231761339253, -6.021691654504919) kcal/(A mol) (-21.95717429857838, -3.2175407957151068, 36.6977967471649) kcal/(A mol) (14.93046186980213, -12.778116481606775, -31.72494740697525) kcal/(A mol) (0.7078827402935988, 8.042910957388216, -1.7323877739688065) kcal/(A mol) (-9.556117795542123, -42.123960606108305, -1.0272121204079496) kcal/(A mol) (-30.317770384062257, 3.6492059461278035, -6.390030268320828) kcal/(A mol) (28.151283018218024, 6.325360626778724, 2.7587154521137403) kcal/(A mol) (-2.415184791493537, 2.437233181027425, 4.111694492038786) kcal/(A mol) (-7.200843021410636, 20.979492400506164, 0.324349359605217) kcal/(A mol) (4.575353981323404, -10.66571222346671, 41.484245869356585) kcal/(A mol) (-4.3124146735344056, 15.310190060842016, -2.763246643147457) kcal/(A mol) (-0.47954099009264683, -2.9714194205113516, -37.665911463707324) kcal/(A mol) (21.061183786042267, -22.617558992289883, 9.306767817875642) kcal/(A mol) (32.848263285091456, -1.8470061270642046, 0.4135650951173963) kcal/(A mol) (4.123430880386882, -15.84444176207767, -2.6947125951569175) kcal/(A mol) (-25.170667253562563, 25.564398772882008, -4.846602549846862) kcal/(A mol) (-33.64438733218469, -12.526521149742683, -33.80410524804523) kcal/(A mol) (-21.324624984240877, -7.519253849484188, 23.80380564246662) kcal/(A mol) (11.3195974785383, 9.74991209678986, -12.377619581954676) kcal/(A mol) (3.233219265138623, -6.549185122371548, -4.9835102205419135) kcal/(A mol) (7.243503789098807, -11.280415442751403, 4.995197893582755) kcal/(A mol) (-4.133456446725416, -25.7152897794951, 57.42764058983423) kcal/(A mol) (0.18429751687391063, 5.369782471760653, -4.90721508006656) kcal/(A mol) (1.3493442896477301, 17.698938042427734, -52.36591102973547) kcal/(A mol) (5.792925501442907, -11.920684260641053, 23.231884598262514) kcal/(A mol) (-4.612678085586851, 13.531516474348988, -2.88822957679308) kcal/(A mol) (9.540039171847129, -10.659253641080046, -4.786150214298018) kcal/(A mol) (-1.8941845125710348, 2.214247103662966, 1.0648110021114723) kcal/(A mol) (-2.285173162595949, -31.55128554302126, 17.705403880780942) kcal/(A mol) (5.256264033325889, -22.612189218932134, 5.316484195739347) kcal/(A mol) (-35.06924503922657, -4.215684608952684, -8.185698061060739) kcal/(A mol) (29.62594925198498, 17.437102227105534, 2.407758081377877) kcal/(A mol) (-16.373001457359017, 21.569840267429303, -12.244328310134897) kcal/(A mol)
from simtk.openmm import LocalEnergyMinimizer
LocalEnergyMinimizer.minimize(simulation.context, 1e-1)
state = simulation.context.getState(getForces=True, getEnergy=True, getPositions=True)
potential_energy = state.getPotentialEnergy()
potential_energy.in_units_of(unit.kilocalorie_per_mole)
Quantity(value=-125.33267648277747, unit=kilocalorie/mole)
kilocalorie_per_mole_per_angstrom = unit.kilocalorie_per_mole/unit.angstrom
for f in state.getForces():
print(f.in_units_of(kilocalorie_per_mole_per_angstrom))
(0.00029612184499987605, -0.0001276382829680238, -7.588497073174688e-05) kcal/(A mol) (-0.0006066711232409618, 0.00041482779529332153, 0.00029284168195442317) kcal/(A mol) (0.0005401988397150456, 0.0006844449312107604, 0.0002715594458234667) kcal/(A mol) (-35.555351640663375, 59.50384661709677, -35.69139857847326) kcal/(A mol) (-0.0005196879062149462, 7.881823649510822e-05, -0.0001762693611465173) kcal/(A mol) (-0.00013258587698162234, -0.0004048991496414358, -0.00019705154522994556) kcal/(A mol) (-0.00011797576614643462, -0.0007752841098020341, 0.0005987259978678855) kcal/(A mol) (-24.225508598267584, 33.12682949478675, 58.46224443437141) kcal/(A mol) (-0.0005122872309905698, 0.0003934222815791098, -6.956581824633793e-05) kcal/(A mol) (0.000989886294828942, 0.0006468202501747544, 3.091765072980642e-05) kcal/(A mol) (-0.00030155130598996297, -0.0005694219263209035, -0.000529555447740272) kcal/(A mol) (32.60446893332641, -28.691834534100725, -37.652556474157585) kcal/(A mol) (-0.0005773072929318601, -0.00021688228963443013, -6.520947206015218e-05) kcal/(A mol) (8.315372471287164e-05, -0.00015472111634221852, 0.0002651799391316153) kcal/(A mol) (-6.379049098525495e-05, 0.00012264366231299576, 0.00019503136650396395) kcal/(A mol) (33.710254647152844, 5.185527409632133, -27.005893815415547) kcal/(A mol) (-0.00044900431519237955, 0.000304329707400328, -0.00010514654835193214) kcal/(A mol) (-0.0001543650949552137, 0.0001790071721944535, -0.00045011027738431526) kcal/(A mol) (-0.0007934144975490045, 0.00043314431671685386, 0.00012003676456828696) kcal/(A mol) (38.79079584024731, 20.454255046159634, 15.072097175191718) kcal/(A mol) (-0.0001286982193724868, -0.0006675173803852515, 0.0005628388081778345) kcal/(A mol) (0.00033000464061240086, -0.00039331606275272603, -1.896510829985661e-05) kcal/(A mol) (-0.00018648483317002218, 0.00010092614059006212, 0.0002550381296488896) kcal/(A mol) (-27.690539265927946, -6.062437621790124, -23.705225520801932) kcal/(A mol) (-0.00016234103626547254, -0.000245527961418829, 0.00045150258975958655) kcal/(A mol) (7.961114958359942e-05, -0.00026754477442618885, -0.0007327984804255176) kcal/(A mol) (0.0002577446167839126, 0.0011821904378852795, 0.0006046245304185761) kcal/(A mol) (4.93034671665006, -51.34641290599498, -3.6529188101311787) kcal/(A mol) (1.573302014047706e-05, -0.00020785072904015315, -0.0003782979879706562) kcal/(A mol) (-0.00019777823454126303, 0.00030909817164474003, 0.00027241906532358664) kcal/(A mol) (-0.00017908878599850493, -0.0004446079906513416, -0.00011361454723818024) kcal/(A mol) (-29.036612533795104, 25.95813707642048, 25.256533589768107) kcal/(A mol) (0.00042839832815443285, 0.00019108014409356216, -0.00011603227122698524) kcal/(A mol) (-0.0009123787933950912, -0.0002579789412986651, 0.00015213953366295622) kcal/(A mol) (0.00026877521839741735, 0.000614102464517063, -0.00047355193193780536) kcal/(A mol) (30.95900112314338, -20.217378534813285, 12.68474277223199) kcal/(A mol) (-0.0005420678144251588, -0.00033989880160781354, -9.775126966566909e-05) kcal/(A mol) (-0.0006278355813283012, 7.859439794775073e-05, -0.0008953689574874138) kcal/(A mol) (3.841710528603946e-05, 0.00010157332418524101, -0.00024607008679804643) kcal/(A mol) (-36.91352581225864, -22.385930443208807, -28.85943671466232) kcal/(A mol) (4.909607107672652e-05, -0.00018003351374709572, 0.00023922664171383915) kcal/(A mol) (0.0002732630854927443, -0.001640904876445954, -4.5977975835956864e-05) kcal/(A mol) (0.0007633035172479743, 0.0010273560369076814, -0.0004743493214465508) kcal/(A mol) (-9.834369071491595, 24.79432560805018, -34.934555076440674) kcal/(A mol) (0.0007132856138798463, -0.00023930818327756764, 0.00023357506600469807) kcal/(A mol) (0.00046012799631624663, 3.359300444903496e-05, -0.00021864668250570952) kcal/(A mol) (2.5275824524949782e-05, 0.00019026789402062474, 0.0009062645607798349) kcal/(A mol) (14.99869222785024, -31.607986749528003, 22.753477598693554) kcal/(A mol) (-0.00030799334692983003, -0.0004587849986137371, -0.00047319323169277283) kcal/(A mol) (5.666070186669188e-05, -0.00013660770033525896, -0.00019379354707875938) kcal/(A mol) (0.00035609199414866677, -0.00026191878355675456, 0.0009292305520827062) kcal/(A mol) (11.164958424853038, -38.597694311281266, 18.00761344164186) kcal/(A mol) (0.0002586987979597095, 0.0009259170982386248, 0.0003530783336402696) kcal/(A mol) (0.0003581665469156733, -6.808656653868162e-05, -0.0004627695859073767) kcal/(A mol) (0.0008312926140726903, 4.657667095145012e-05, -0.00012425623155255458) kcal/(A mol) (-35.08017697136835, 29.565577400244905, -15.329743617099835) kcal/(A mol)
state.getPositions()
Quantity(value=[(0.003560411538791464, -0.0003382258050843299, -0.005766652216817242), (-0.009917813375329295, -0.09518462361639908, -0.022004942227067264), (0.03205102826239977, 0.007771279116298518, 0.08692282891829316), (0.00676335698898585, -0.018843840364799736, 0.01054447247822038), (-0.24890932435322186, 0.08749172040129859, 0.007675088727051557), (-0.15610996556353152, 0.056143001306187455, 0.0029521876978927677), (-0.2752483602359659, 0.07794280554129548, -0.08564747563060827), (-0.23472978577105377, 0.07876606946821653, -0.013243248687672592), (0.044545320842960257, 0.026547174463420056, 0.26478674114362216), (0.030481233289616796, 0.12211929021916271, 0.27635392109179685), (-0.041106568008704436, -0.014727833363660223, 0.28148227640035567), (0.023270574154276955, 0.03813164927634856, 0.27081668868459235), (0.051249112372636026, 0.20633694238362976, -0.17628490339801928), (0.02566658114793184, 0.1283081190931048, -0.12419405939243018), (0.01086269331019525, 0.2791842918347074, -0.12788610388967658), (0.037174409829715704, 0.2052314571105322, -0.15484509216945863), (-0.13836479590059775, 0.3374058833655197, 0.02016963385576341), (-0.1972342558108169, 0.4104054062869363, 0.0009474997648938342), (-0.19343760381207648, 0.25652564781686893, 0.015836226006853187), (-0.1626747708157006, 0.33572450620896716, 0.015143977928291577), (0.18370780562100777, -0.18241898377775012, -0.21668970172162363), (0.2008675433100884, -0.23929178521648545, -0.14071715209874094), (0.22293826609510684, -0.09750185560139285, -0.1929713654904698), (0.19573884843779407, -0.17643563016775018, -0.19542030818832673), (0.28611312144013984, -0.07199763851935273, 0.1344153151257371), (0.22052161153558292, -0.03441464888699876, 0.1943629225738926), (0.28875920601121763, -0.012899784737546466, 0.057755878307976474), (0.2726834983500335, -0.05137044789815946, 0.130849788824327), (-0.21930758728616118, 0.05931199922950993, -0.2668733280621772), (-0.14648410834664544, 0.11913522299826179, -0.2856981200225727), (-0.17980769920148168, -0.02938007022237302, -0.27530279890311293), (-0.1953430104268432, 0.05315273166531592, -0.2726881188127595), (-0.21883390189002913, -0.10087357118505721, 0.2148318813381277), (-0.255841414292043, -0.040564883064328275, 0.14802882683266164), (-0.2922457598930669, -0.1248716322282096, 0.27139005603262273), (-0.24239225451900048, -0.09312657395318848, 0.2126461009778256), (0.2025366613894322, -0.31244457476190074, 0.04594374470808561), (0.23957103865262855, -0.2304281451111624, 0.08509808734575212), (0.25983525160105164, -0.38387297513768254, 0.0735795206996376), (0.22266292459992523, -0.31018558228268006, 0.060193646342427576), (-0.07238058137976466, -0.26699047967274747, 0.018817518625549225), (0.012610078168482836, -0.3037399979491179, 0.04697247535332153), (-0.12120120329579429, -0.2426605606705672, 0.09828536725585695), (-0.06466357938983149, -0.2696402439137371, 0.04177921447876884), (0.01096995472579537, 0.300257415454047, 0.25453123036489905), (-0.006937264225554583, 0.3653205840421358, 0.32231084645257635), (-0.04154120752216275, 0.3263239748219668, 0.17721329202800984), (-0.004054049212726713, 0.3197002564908554, 0.2524961964916629), (0.29214178220878156, 0.07444034733747983, -0.11894106041828818), (0.22311578992308412, 0.13927446473237273, -0.14123838178214748), (0.3734166785459871, 0.10905246658759625, -0.15581799447313172), (0.29475512803664966, 0.09565754412019287, -0.1315660914304813), (-0.09647527054469407, -0.1923970760588945, -0.26285901249760374), (0.0004845512207284013, -0.1883468743428491, -0.2634296551326173), (-0.11542103627579053, -0.23935181710225986, -0.1809406592939449), (-0.0798307032434038, -0.20155090715889806, -0.24550319853031122)], unit=nanometer)
simulation.context.setVelocitiesToTemperature(300*unit.kelvin)
# Equilibrate
simulation.step(10)
Add a reporter
that prints out the simulation status every 10 steps
simulation.reporters.append(app.StateDataReporter(sys.stdout, 10, step=True,
potentialEnergy=True, temperature=True, progress=True, remainingTime=True,
speed=True, totalSteps=110, separator='\t'))
Add a PDBReporter
that writes molecules positions every 20 steps in a pdb file.
simulation.reporters.append(app.PDBReporter('trajectory.pdb', 20))
Run 100 steps
simulation.step(100)
#"Progress (%)" "Step" "Potential Energy (kJ/mole)" "Temperature (K)" "Speed (ns/day)" "Time Remaining" 18.2% 20 -524.391917919 236.419400636 0 -- 27.3% 30 -524.391917285 236.419399841 3.83e-06 0:18 36.4% 40 -524.391916389 236.419398258 3.84e-06 0:15 45.5% 50 -524.391915232 236.419396147 3.84e-06 0:13 54.5% 60 -524.391913813 236.419393904 3.85e-06 0:11 63.6% 70 -524.391912133 236.419391369 3.85e-06 0:08 72.7% 80 -524.391910192 236.41938822 3.85e-06 0:06 81.8% 90 -524.391907989 236.419384268 3.85e-06 0:04 90.9% 100 -524.391905526 236.419379902 3.85e-06 0:02 100.0% 110 -524.3919028 236.419375026 3.85e-06 0:00
!head trajectory.pdb
REMARK 1 CREATED WITH OPENMM 6.0, 2014-07-15 MODEL 1 ATOM 1 O HOH A 1 0.036 -0.003 -0.058 1.00 0.00 O ATOM 2 H1 HOH A 1 -0.099 -0.952 -0.220 1.00 0.00 H ATOM 3 H2 HOH A 1 0.321 0.078 0.869 1.00 0.00 H ATOM 4 M HOH A 1 0.068 -0.188 0.105 1.00 0.00 ATOM 5 O HOH A 2 -2.489 0.875 0.077 1.00 0.00 O ATOM 6 H1 HOH A 2 -1.561 0.561 0.030 1.00 0.00 H ATOM 7 H2 HOH A 2 -2.752 0.779 -0.856 1.00 0.00 H ATOM 8 M HOH A 2 -2.347 0.788 -0.132 1.00 0.00
!echo Number of lines: `wc -l trajectory.pdb`
Number of lines: 263 trajectory.pdb