openmm / pdbfixer

PDBFixer fixes problems in PDB files
Other
453 stars 114 forks source link

Potential Energy mantisa is same for all iteration. #242

Closed yhgon closed 2 years ago

yhgon commented 2 years ago

I found odd report what Potential Energy mantisa is same during whole iteration. X..806077383 as below :

#"Progress (%)","Step","Potential Energy (kJ/mole)","Kinetic Energy (kJ/mole)","Total Energy (kJ/mole)","Temperature (K)","Speed (ns/day)","Elapsed Time (s)","Time Remaining"
0.1%,100,-5566615.806077383,758098.6629617857,-4808517.143115598,222.6030279483443,0,0.0003311634063720703,--
0.1%,200,-5523309.806077383,810809.0738813787,-4712500.732196005,238.0805609507925,15.2,1.1333811283111572,37:44
0.1%,300,-5488829.806077383,851931.4012187873,-4636898.404858597,250.15544648854942,15,2.3056721687316895,38:22
0.2%,400,-5452543.806077383,880249.3672351063,-4572294.438842277,258.4705448900469,15.1,3.4337408542633057,38:04
0.2%,500,-5425341.806077383,902689.0839320947,-4522652.722145289,265.0595934230404,15.2,4.547600746154785,37:48
0.3%,600,-5397853.806077383,918830.4644251866,-4479023.341652197,269.79924058056304,15.3,5.6640002727508545,37:38
0.3%,700,-5373821.806077383,934447.826814159,-4439373.979263225,274.38501856198013,15.4,6.71558952331543,37:10
0.4%,800,-5350193.806077383,943267.2487715428,-4406926.557305841,276.9746947194313,15.5,7.811719179153442,37:02
0.5%,900,-5331357.806077383,953770.0557347691,-4377587.750342614,280.0586688064413,15.6,8.845742464065552,36:41

test code is same as openMM tutorial except I use CUDA single precision.


conda install -y -q -c conda-forge openmm pdbfixer python=3.7 > /dev/null

from openmm.app import *
from openmm import *
from openmm.unit import *
from sys import stdout

platform = Platform.getPlatformByName('CUDA')
properties = {'DeviceIndex': '0', 'Precision': 'single'}
forcefield = ForceField('amber99sb.xml', 'tip3p.xml')

pdb_water='6hqv_A_water.pdb'
pdb_output = '6hqv_A_traj3_300k_001p_total_200K_1k_iter.pdb'

pdb = PDBFile(pdb_water)

system = forcefield.createSystem(pdb.topology, nonbondedMethod=PME, 
                                 nonbondedCutoff=1.0*nanometer, constraints=HBonds)
integrator = LangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds) # temperature, friction, timestep 
simulation1 = Simulation(pdb.topology, system, integrator, platform, properties)
simulation1.context.setPositions(pdb.positions)
simulation1.minimizeEnergy(maxIterations=25)

steps = 200000 # 100k 20min, 10k 2min
saveInterval = 1000
printInterval = 100

simulation1.reporters.append(PDBReporter( file = pdb_output, reportInterval = saveInterval, enforcePeriodicBox =False))
simulation1.reporters.append(StateDataReporter(stdout, printInterval, 
                                              step=True,  
                                              totalSteps = steps, 
                                              potentialEnergy=True, 
                                              kineticEnergy =True, 
                                              totalEnergy =True, 
                                              temperature=True,
                                              elapsedTime =True,
                                              remainingTime =True,
                                              progress =True,
                                              speed=True ))

simulation1.step(steps)
yhgon commented 2 years ago

another simulation result show below.

#"Progress (%)","Step","Potential Energy (kJ/mole)","Kinetic Energy (kJ/mole)","Total Energy (kJ/mole)","Temperature (K)","Speed (ns/day)","Elapsed Time (s)","Time Remaining"
0.1%,100,-121370.49897545716,49133.63804504863,-72236.86093040853,193.83086050003715,0,0.00028634071350097656,--
0.1%,200,-118733.46772545716,55663.14637669065,-63070.321348766505,219.58959258097806,144,0.12048935890197754,4:00
0.1%,300,-116580.84272545716,60875.83453027613,-55705.00819518103,240.1535050150919,145,0.23886895179748535,3:58
0.2%,400,-113773.78022545716,63534.93906139524,-50238.841164061916,250.6435997838489,146,0.3542344570159912,3:55
0.2%,500,-112449.37397545716,66616.5959524457,-45832.77802301146,262.80065207479714,147,0.4690840244293213,3:53
0.3%,600,-111380.03022545716,69195.02884695033,-42185.00137850683,272.9724994398387,148,0.5849928855895996,3:53
0.3%,700,-109970.49897545716,71040.78516078438,-38929.71381467278,280.2539576997748,148,0.6993200778961182,3:52
0.4%,800,-109522.93647545716,72284.3991101383,-37238.53736531886,285.1599807732568,149,0.8135404586791992,3:51
0.5%,900,-108419.09272545716,72012.97321176575,-36406.1195136914,284.08921301542637,149,0.9291989803314209,3:51
0.5%,1000,-108747.96772545716,73945.91190383676,-34802.05582162039,291.71460337700364,119,1.3119151592254639,4:50
0.6%,1100,-108387.65522545716,74383.8260657452,-34003.82915971195,293.4421627885425,121,1.4302194118499756,4:44
0.6%,1200,-108362.46772545716,75495.17901755986,-32867.2887078973,297.8264198381016,123,1.5439362525939941,4:39
0.7%,1300,-107781.31147545716,75015.68742250348,-32765.624052953674,295.93483858806064,125,1.6578710079193115,4:34
0.7%,1400,-108697.03022545716,76548.92161350226,-32148.1086119549,301.98340560678804,127,1.772719383239746,4:30

it seems that inital value is not memset(0) and only X and Y part is updated in the value XXXXX.YYYYYZZZZZZZ 5 digists of Y is related with single precision.

-121370.49897545716,
-118733.46772545716,
peastman commented 2 years ago

That's because single precision only has about eight significant digits. The energy gets converted to double precision when it's returned, but any additional digits beyond that are meaningless.