#!/bin/env python """ Tools for constructing systems from flat files. @author Mark Friedrichs """ #============================================================================================= # GLOBAL IMPORTS #============================================================================================= import os import os.path import sys import copy import re import math import traceback import simtk.openmm import time from ValidationUtilities import ValidationUtilities #============================================================================================= # Results class for energy conservation tests #============================================================================================= class EnergyConservationResult(object): """ EnergyConservationResult generates/computes results for energy conservation tests """ def __init__(self, integratorTest, platform , verbose=False): """ Generate results summary ARGUMENTS integratorTest(EnergyConservationTest) - EnergyConservationTest platform (string) - platform name verbose (Boolean) - if true, print logging info """ self._platform = platform self._integratorTest = integratorTest self._openMM = simtk.openmm self._verbose = verbose self._permanentVerbose = verbose self._validationUtilities = ValidationUtilities(verbose) self._minNormCutoff = 1.0e-02 self._systemName = integratorTest.getSystemName() self._testName = integratorTest.getTestName() self._tolerance = integratorTest.getTolerance() self._totalTime = integratorTest.getTotalTime() self._totalSteps = integratorTest.getTotalSteps() self._stepSize = integratorTest.getStepSize() self._passed = -1 self._numberOfConstraintChecks = integratorTest._numberOfConstraintChecks self._constraintViolations = integratorTest._constraintViolations self._constraintTolerance = integratorTest._constraintTolerance self._maximumConstraintViolation = integratorTest._maximumConstraintViolation self._indexOfMaximumConstraintViolation = integratorTest._indexOfMaximumConstraintViolation self._stepIndexOfMaximumConstraintViolation = integratorTest._stepIndexOfMaximumConstraintViolation self._averageVolume = integratorTest.getAverageVolume() self._stddevVolume = integratorTest.getStddevVolume() # get units from ValidationUtilities lengthUnit = self._validationUtilities.getLengthUnit( ) energyUnit = self._validationUtilities.getEnergyUnit( ) dof = integratorTest.getDegreesOfFreedom() Boltzmann = (1.380658e-23)*(6.0221367e23)*0.001 conversionFactor = 2.0/(dof*Boltzmann) # compute drift in energy for Verlet integrators # compute temperature stability for Langevin integrators integratorName = self._integratorTest.getIntegratorName() if( integratorName == 'VerletIntegrator' or integratorName == 'VariableVerletIntegrator' ): self._isVerlet = 1 else: self._isVerlet = 0 if( self._isVerlet ): # compute average/std dev of total energy mean = 0.0 stddev = 0.0 meanKE = 0.0 count = len( self._integratorTest._kineticEnergies ) for ii in range( count ): totalEnergy = self._integratorTest._kineticEnergies[ii] + self._integratorTest._potentialEnergies[ii] mean += totalEnergy stddev += totalEnergy*totalEnergy meanKE += self._integratorTest._kineticEnergies[ii] if( count > 0 ): mean = mean/count meanKE = meanKE/count if( count > 1 ): stddev = stddev - mean*mean*count stddev = math.sqrt( stddev/(count - 1) ) # compute drift (std dev energy)/(kT/dof/ns of simulation) self._temperature = meanKE*conversionFactor; kT = self._temperature*Boltzmann; self._drift = stddev*1000.0/(kT*dof*self._integratorTest._times[count-1]) if self._verbose: print "Verlet: drift=%14.7e stddev=%14.7e meanE=%14.7e meanKE=%14.7e " % ( self._drift, stddev, mean, meanKE) sys.stdout.flush() else: # compute average temperature meanKE = 0.0 stddevKE = 0.0 count = len( self._integratorTest._kineticEnergies ) for ii in range( count ): meanKE += self._integratorTest._kineticEnergies[ii] stddevKE += (self._integratorTest._kineticEnergies[ii]*self._integratorTest._kineticEnergies[ii]) if( count > 0 ): meanKE = meanKE/count if( count > 1 ): stddevKE = stddevKE - meanKE*meanKE*count stddevKE = math.sqrt( stddevKE/(count - 1) ) self._temperature = meanKE*conversionFactor; self._temperatureStddev = stddevKE*conversionFactor self._deltaT = abs( self._temperature - self._integratorTest._simulationTemperature ) self._relativeDeltaT = self._deltaT/( self._temperature + self._integratorTest._simulationTemperature ) if self._verbose: print "Stochastic: Calc.T=%14.7e specifiedT=%14.7e delta=%14.7e %14.7e" % ( self._temperature, self._integratorTest._simulationTemperature, self._deltaT, self._relativeDeltaT) sys.stdout.flush() # report average volume if non-zero if( self._verbose and self._averageVolume > 0.0 ): print "Average Volume=%14.7e stddev=%14.7e" % ( self._averageVolume, self._stddevVolume) sys.stdout.flush() #============================================================================================= # Turn off verbosity #============================================================================================= def verboseOff(self): self._verbose = False return #============================================================================================= # Return number of constraint violations #============================================================================================= def getNumberOfConstraintViolations(self): return self._constraintViolations #============================================================================================= # Return number of times constraints were checked #============================================================================================= def getNumberOfConstraintChecks(self): return self._numberOfConstraintChecks #============================================================================================= # Return constraint tolerance #============================================================================================= def getConstraintTolerance(self): return self._constraintTolerance #============================================================================================= # Return maximum constraint violation #============================================================================================= def getMaximumConstraintViolation(self): return self._maximumConstraintViolation #============================================================================================= # Return maximum constraint violation #============================================================================================= def getConstraintIndexOfMaximumConstraintViolation(self): return self._indexOfMaximumConstraintViolation #============================================================================================= # Return step on which maximum constraint violation was observed #============================================================================================= def getConstraintIndexOfMaximumConstraintViolation(self): return self._stepIndexOfMaximumConstraintViolation #============================================================================================= # Get tolerance #============================================================================================= def getTolerance(self): return self._tolerance #============================================================================================= # Get platform #============================================================================================= def getPlatform(self): return self._platform #============================================================================================= # Restore verbosity #============================================================================================= def verboseRestore(self): self._verbose = self._permanentVerbose return #============================================================================================= # Get energy drift #============================================================================================= def getDrift(self): return self._drift #============================================================================================= # Get total time #============================================================================================= def getTotalTime(self): return self._totalTime #============================================================================================= # Get total number of steps #============================================================================================= def getTotalSteps(self): return self._totalSteps #============================================================================================= # Get time/step #============================================================================================= def getTimePerStep(self): if( self._totalSteps > 0 ): return self._totalTime/self._totalSteps else: return 0.0 #============================================================================================= # Get ns/day #============================================================================================= def getNsPerDay(self): if( self._totalTime> 0 ): simulationTime = self._totalSteps*self._stepSize nsPerDay = 86.4*simulationTime/self._totalTime return nsPerDay else: return 0.0 #============================================================================================= # Get delta temperature #============================================================================================= def getDeltaT(self): return self._deltaT #============================================================================================= # Get relative delta T #============================================================================================= def getRelativeDeltaT(self): return self._relativeDeltaT #============================================================================================= # Get flag signalling whether test used Verlet integrator #============================================================================================= def isVerlet(self): return self._isVerlet #============================================================================================= # Get average temeperature #============================================================================================= def getMeanTemperature(self): return self._temperature #============================================================================================= # Get system name #============================================================================================= def getSystemName(self): return self._systemName #============================================================================================= # Get test type #============================================================================================= def getTestName(self): return self._testName #============================================================================================= # Return 1 if test passsed, else 0 #============================================================================================= def testPassed(self): if( self._passed > -1 ): return self._passed self._passed = 1 if( self.isVerlet() and self.getDrift() > self._tolerance ): if( self._verbose ): print "%30s %30s failed: drift=%14.7e is greater than tolerance=%10.3e" % (self.getTestName(), self.getSystemName(), self.getDrift(), self._tolerance ) self._passed = 0 if( self.isVerlet() == 0 and self.getRelativeDeltaT() > self._tolerance ): if( self._verbose ): print "%30s %30s failed: relative delta T=%14.7e is greater than tolerance=%10.3e" % (self.getTestName(), self.getSystemName(), self.getRelativeDeltaT(), self._tolerance ) self._passed = 0 if( self._verbose and self._passed ): print "%30s %30s self._passed" % (self.getTestName(), self.getSystemName() ) return self._passed #============================================================================================= # Base class for ForceEnergy comparison tests #============================================================================================= class EnergyConservationTest(object): """Base class for energy conservation tests Centralizes common functionality across energy conservation tests """ def __init__(self, parameterInput, defaultOverrides, verbose=False): """ ARGUMENTS defaultOverrides (dictionary of overrides) - verbose (Boolean) - if true, be verbose """ self._isActive = 1 self._openMM = simtk.openmm self._verbose = verbose self._validationUtilities = ValidationUtilities(verbose) self._dof = 0 self._integratorName = 'NA' self._totalTime = 0 self._totalSteps = 0 self._stepSize = 0 # constraint-related fields self._numberOfConstraintChecks = 0 self._constraintViolations = 0 self._constraintTolerance = 1.0e-04 self._maximumConstraintViolation = -1.0e+10 self._indexOfMaximumConstraintViolation = -1 self._stepIndexOfMaximumConstraintViolation = 0 self._calculateVolume = 0 self._volume = 0.0 self._volume2 = 0.0 self._volumeCount = 0 # set parameters if( 'DeviceId' in defaultOverrides ): self._deviceId = defaultOverrides['DeviceId'] else: self._deviceId = -1 if( 'EquilibrationTotalSteps' in defaultOverrides ): self._equilibrationTotalSteps = defaultOverrides['EquilibrationTotalSteps'] else: self._equilibrationTotalSteps = 2000 if( 'SimulationStepSize' in defaultOverrides ): self._simulationStepSize = defaultOverrides['SimulationStepSize'] else: self._simulationStepSize = 0.001 if( 'SimulationTotalSteps' in defaultOverrides ): self._simulationTotalSteps = defaultOverrides['SimulationTotalSteps'] else: self._simulationTotalSteps = 100000 if( 'SimulationStepsBetweenReportsRatio' in defaultOverrides ): self._simulationReportRatio = defaultOverrides['SimulationStepsBetweenReportsRatio'] else: self._simulationReportRatio = 0.01 if( 'SimulationErrorTolerance' in defaultOverrides ): self._simulationErrorTolerance = defaultOverrides['SimulationErrorTolerance'] else: self._simulationErrorTolerance = 1.0e-06 if( 'SimulationConstraintTolerance' in defaultOverrides ): self._equilibrationConstraintTolerance = defaultOverrides['SimulationConstraintTolerance'] else: self._equilibrationConstraintTolerance = 1.0e-06 self._systemName = parameterInput.getName().replace( '.txt', '' ) self._testName = defaultOverrides['TestName'] self._positions = parameterInput.getPositions() # build system self._system = self._validationUtilities.copySystem( parameterInput ) self.editSystemBasedOnTestName() self._dof = self._validationUtilities.getDegreesOfFreedom( parameterInput ) #============================================================================================= # Return nonzero if test is active #============================================================================================= def isActive(self): return self._isActive #============================================================================================= # Get full test name (system tested and test) #============================================================================================= def getFullTestName(self): return self._systemName + '_' + self._testName #============================================================================================= # Get system name #============================================================================================= def getSystemName(self): return self._systemName #============================================================================================= # Get test type #============================================================================================= def getTestName(self): return self._testName #============================================================================================= # Get result #============================================================================================= def getResult(self): return self._result #============================================================================================= # Get degrees-of-freedom #============================================================================================= def getDegreesOfFreedom(self): return self._dof #============================================================================================= # Get integrator name #============================================================================================= def getIntegratorName(self): return self._integratorName #============================================================================================= # Get tolerance #============================================================================================= def getTolerance(self): return self._tolerance #============================================================================================= # Get SimulationStepSize #============================================================================================= def getSimulationStepSize(self): return self._simulationStepSize #============================================================================================= # Get total time #============================================================================================= def getTotalTime(self): return self._totalTime #============================================================================================= # Get step size #============================================================================================= def getStepSize(self): return self._stepSize #============================================================================================= # Get average volume #============================================================================================= def getAverageVolume(self): if( self._volumeCount > 0 ): return (self._volume/self._volumeCount) else: return 0.0 #============================================================================================= # Get volume std dev #============================================================================================= def getStddevVolume(self): if( self._volumeCount > 1 ): average = (self._volume/self._volumeCount) stddev = self._volume2 - average*average*self._volumeCount stddev = math.sqrt( stddev/(self._volumeCount-1 ) ) return stddev else: return 0.0 #============================================================================================= # Get total number of steps #============================================================================================= def getTotalSteps(self): return self._totalSteps #============================================================================================= # Get time/step #============================================================================================= def getTimePerStep(self): if( self._totalSteps > 0 ): return self._totalTime/self._totalSteps else: return 0.0 #============================================================================================= # Get dimensions of box enclosing system #============================================================================================= def getEnclosingBox(self): ii = 0 self._range = [ [ 1.0e+30, 1.0e+30, 1.0e+30 ], [ -1.0e+30 , -1.0e+30, -1.0e+30 ] ] while( ii < len( self._positions ) ): coordinates = self._positions[ii] jj = 0 while( jj < 3): if( self._range[0][jj] > coordinates[jj] ): self._range[0][jj] = coordinates[jj] if( self._range[1][jj] < coordinates[jj] ): self._range[1][jj] = coordinates[jj] jj += 1 #print "%5d [%14.7f %14.7f %14.7f]" % ( ii, coordinates[0], coordinates[1], coordinates[2] ) ii += 1 #if( ii == maxPrint and len( self._positions ) > 2*maxPrint ):ii = len( self._positions ) - maxPrint if self._verbose: jj = 0 while( jj < 3): print "Enclosing Box %2d [%14.7f %14.7f] %14.7f" % ( jj, self._range[0][jj], self._range[1][jj], (self._range[1][jj] - self._range[0][jj]) ) jj += 1 def getBoxDimensions(self): self.getEnclosingBox() self._box = [ 0, 0, 0 ] jj = 0 while( jj < 3): self._box[jj] = self._range[1][jj] - self._range[0][jj] jj += 1 #============================================================================================= # Check constraints; return number of violations #============================================================================================= def checkConstraints(self, context, step ): # any active constraints if( self._system.getNumConstraints() < 1 ): return 0 state = context.getState(getPositions=True) positions = state.getPositions() if( self._verbose and (self._numberOfConstraintChecks == 0) ): print "checkConstraints: %s particles=%d constraints=%d" % (self.getFullTestName(), self._system.getNumParticles(), self._system.getNumForces() ) self._numberOfConstraintChecks += 1 numberOfViolations = 0 lengthUnit = self._validationUtilities.getLengthUnit( ) lengthUnit2 = self._validationUtilities.getLengthUnit( )*self._validationUtilities.getLengthUnit( ) for ii in range( self._system.getNumConstraints() ): args = self._system.getConstraintParameters( ii ) particle1 = args[0] particle2 = args[1] expectedDistance = args[2]/lengthUnit currentDistance = (positions[particle1][0] - positions[particle2][0])*(positions[particle1][0] - positions[particle2][0]) + (positions[particle1][1] - positions[particle2][1])*(positions[particle1][1] - positions[particle2][1]) + (positions[particle1][2] - positions[particle2][2])*(positions[particle1][2] - positions[particle2][2]) currentDistance /= lengthUnit2 currentDistance = math.sqrt( currentDistance ) delta = abs( currentDistance - expectedDistance ) if( delta > self._constraintTolerance ): self._constraintViolations += 1 numberOfViolations += 1 if( delta > self._maximumConstraintViolation ): self._maximumConstraintViolation = delta self._indexOfMaximumConstraintViolation = ii self._stepIndexOfMaximumConstraintViolation = step return numberOfViolations #============================================================================================= # Accumulate volume of enclosing box #============================================================================================= def accumulateVolume(self, context ): lengthUnit = self._validationUtilities.getLengthUnit( ) box = context.getState(0).getPeriodicBoxVectors() volume = (box[0][0]/lengthUnit)*(box[1][1]/lengthUnit)*(box[2][2]/lengthUnit) self._volume += volume self._volume2 += volume*volume self._volumeCount += 1 #if self._verbose: print "accumulateVolume: volume=%12.3e %12.3e %d [%12.3e %12.3e %12.3e]" % ( volume, self._volume, self._volumeCount, box[0][0]/lengthUnit,box[1][1]/lengthUnit,box[2][2]/lengthUnit ) return #============================================================================================= # Serialize #============================================================================================= def serialize(self, serializeDirectory ): if self._system is None: print "serialize: %s no system available" % (self.getFullTestName() ) return # serialize system self._validationUtilities.serialize( self.getFullTestName(), self._system, self._positions, serializeDirectory, self.getFullTestName() ) return #============================================================================================= # Run test #============================================================================================= def runTest(self, platformName ): if self._verbose: print "runTest: %s Particles=%d NumForces=%d" % (self.getFullTestName(), self._system.getNumParticles(), self._system.getNumForces() ) sys.stdout.flush() energyUnit = self._validationUtilities.getEnergyUnit() timeUnit = self._validationUtilities.getTimeUnit() self._platformName = platformName; self._integrator = self.getIntegrator( self._integratorName ) platform = self._openMM.Platform.getPlatformByName( platformName ) deviceId = self._deviceId if( deviceId > -1 ): self._validationUtilities.setDeviceId( platform, deviceId ) context = self._openMM.Context(self._system, self._integrator, platform ) context.setPositions( self._positions ) self._validationUtilities.printPlatformProperties( context ) self._stepSize = self._integrator.getStepSize()/timeUnit # steps between recording of energies self._stepsPerChunk = int( self._simulationTotalSteps*self._simulationReportRatio ) if( self._stepsPerChunk < 1 ): self._stepsPerChunk= 1 numberOfChunks = int( self._simulationTotalSteps/self._stepsPerChunk ) if self._verbose: print "runTest: %s Total steps=%d stepsPerChunk=%d numberOfChunks=%d" % (self.getFullTestName(), self._simulationTotalSteps, self._stepsPerChunk, numberOfChunks ) sys.stdout.flush() totalSteps = 0 self._kineticEnergies = [] self._potentialEnergies = [] self._times = [] self._steps = [] # equilibration -- collect initial energies self._integrator.step( self._equilibrationTotalSteps ) state = context.getState(getEnergy=True) baseTime = state.getTime()/timeUnit self._times.append( 0 ) kineticEnergy = state.getKineticEnergy()/energyUnit self._kineticEnergies.append( kineticEnergy ) potentialEnergy = state.getPotentialEnergy()/energyUnit self._potentialEnergies.append( potentialEnergy ) self._steps.append( 0 ) if( self._verbose ): totalEnergy = kineticEnergy + potentialEnergy print "Equilibrate: %30s %8d %14.7e %14.7e %14.7e %14.7e" % (self.getFullTestName(), self._equilibrationTotalSteps, baseTime, kineticEnergy, potentialEnergy, totalEnergy) sys.stdout.flush() # simulation chunk = 0 printChunk = int(0.1*numberOfChunks) totalTime = 0 timeErrors = 0 localtime = time.strftime("%b %d %Y %H:%M:%S", time.localtime(time.time())) print "Local current time :", localtime loopStartTime = time.time() while( totalSteps < self._simulationTotalSteps ): if( (totalSteps + self._stepsPerChunk) > self._simulationTotalSteps ): self._stepsPerChunk = self._simulationTotalSteps - totalSteps try: startTime = time.time() except: startTime = 0.0 # get state before timing for OpenCL runs self._integrator.step( self._stepsPerChunk ) state = context.getState(getEnergy=True) try: now = time.time() totalTime += now - startTime #print "Time now=%s start=%s delta=%s ttl=%12.3e" % (str(now), str(startTime), str(now-startTime), totalTime) except: if( timeErrors < 10 ): print "Time.time() not working" print repr(traceback.print_exception( sys.exc_info()[0], sys.exc_info()[1], sys.exc_info()[2])) timeErrors = timeErrors + 1 simulationTime = state.getTime()/timeUnit - baseTime self._times.append( simulationTime) kineticEnergy = state.getKineticEnergy()/energyUnit self._kineticEnergies.append( kineticEnergy ) potentialEnergy = state.getPotentialEnergy()/energyUnit self._potentialEnergies.append( potentialEnergy ) totalSteps += self._stepsPerChunk self._steps.append( totalSteps ) # check for constraint violations constraintViolations = self.checkConstraints( context, totalSteps ) # accumulate volume if( self._calculateVolume ): self.accumulateVolume( context ) chunk += 1 if( self._verbose and ( (numberOfChunks <= 10) or (chunk == printChunk) ) ): totalEnergy = kineticEnergy + potentialEnergy if( totalSteps <= 0 ): totalSteps = 1 if( totalTime > 0 ): if( self.getFullTestName().startswith('Variable') ): nsPerDay = (86.4*simulationTime)/totalTime else: nsPerDay = (86.4*self._stepSize*totalSteps)/totalTime else: nsPerDay = 0.0 print "Simulation : %30s %8d %14.7e %14.7e %14.7e %14.7e Cnst=%6d Time=%10.3e ns/day=%10.3e" % (self.getFullTestName(), totalSteps, simulationTime, kineticEnergy, potentialEnergy, totalEnergy, constraintViolations, totalTime, nsPerDay) # print "Simulation : %30s %8d %14.7e %14.7e %14.7e %14.7e Cnst=%6d" % (self.getFullTestName(), totalSteps, simulationTime, kineticEnergy, potentialEnergy, totalEnergy, constraintViolations ) chunk = 0 sys.stdout.flush() loopStopTime = time.time() localtime = time.strftime("%b %d %Y %H:%M:%S", time.localtime(time.time())) print "Local current time :", localtime # get result self._totalTime = totalTime self._totalLoopTime = loopStopTime - loopStartTime self._totalSteps = totalSteps if( self._totalTime > 0 ): nsPerDay = (86.4*self._stepSize*totalSteps)/self._totalTime else: nsPerDay = 0.0 if( self._totalLoopTime > 0 ): nsPerDayLoop = (86.4*self._stepSize*totalSteps)/self._totalLoopTime else: nsPerDayLoop = 0.0 print "Loop time delta=%s ns/day=%12.3e ns/day(Loop)=%12.3e" % (str(loopStopTime-loopStartTime), nsPerDay, nsPerDayLoop) sys.stdout.flush() self._result = EnergyConservationResult( self, platformName, self._verbose ) return #============================================================================================= # Run test on variable integrator #============================================================================================= def runVariableTest(self, platformName ): if self._verbose: print "runVariableTest: %s Particles=%d NumForces=%d" % (self.getFullTestName(), self._system.getNumParticles(), self._system.getNumForces() ) sys.stdout.flush() energyUnit = self._validationUtilities.getEnergyUnit() timeUnit = self._validationUtilities.getTimeUnit() self._platformName = platformName; self._integrator = self.getIntegrator( self._integratorName ) platform = self._openMM.Platform.getPlatformByName( platformName ) deviceId = self._deviceId if( deviceId > -1 ): self._validationUtilities.setDeviceId( platform, deviceId ) context = self._openMM.Context(self._system, self._integrator, platform ) context.setPositions( self._positions ) self._validationUtilities.printPlatformProperties( context ) self._stepSize = self._integrator.getStepSize()/timeUnit # steps between recording of energies self._totalSimulationTime = self._simulationStepSize*self._simulationTotalSteps self._timePerChunk = int( self._totalSimulationTime*self._simulationReportRatio ) numberOfChunks = int( self._totalSimulationTime/self._timePerChunk ) if self._verbose: print "runTest: %s Total time=%d timePerChunk=%d numberOfChunks=%d" % (self.getFullTestName(), self._totalSimulationTime, self._timePerChunk, numberOfChunks ) sys.stdout.flush() self._kineticEnergies = [] self._potentialEnergies = [] self._times = [] self._steps = [] totalSteps = 0 # equilibration -- collect initial energies self._integrator.step( self._equilibrationTotalSteps ) state = context.getState(getEnergy=True) baseTime = state.getTime()/timeUnit self._times.append( 0 ) kineticEnergy = state.getKineticEnergy()/energyUnit self._kineticEnergies.append( kineticEnergy ) potentialEnergy = state.getPotentialEnergy()/energyUnit self._potentialEnergies.append( potentialEnergy ) self._steps.append( 0 ) if( self._verbose ): totalEnergy = kineticEnergy + potentialEnergy print "Equilibrate: %30s %8d %14.7e %14.7e %14.7e %14.7e" % (self.getFullTestName(), self._equilibrationTotalSteps, baseTime, kineticEnergy, potentialEnergy, totalEnergy) sys.stdout.flush() totalSimulationTime = baseTime; # simulation chunk = 0 printChunk = int(0.1*numberOfChunks) totalTime = 0 timeErrors = 0 localtime = time.strftime("%b %d %Y %H:%M:%S", time.localtime(time.time())) print "Local current time :", localtime loopStartTime = time.time() while( totalSimulationTime < (self._totalSimulationTime + baseTime) ): totalSimulationTime += self._timePerChunk; if( (totalSimulationTime - baseTime) > self._totalSimulationTime ): totalSimulationTime = self._totalSimulationTime + baseTime; try: startTime = time.time() except: startTime = 0.0 # get state before timing for OpenCL runs self._integrator.stepTo( totalSimulationTime ) state = context.getState(getEnergy=True) try: now = time.time() totalTime += now - startTime #print "Time now=%s start=%s delta=%s ttl=%12.3e" % (str(now), str(startTime), str(now-startTime), totalTime) except: if( timeErrors < 10 ): print "Time.time() not working" print repr(traceback.print_exception( sys.exc_info()[0], sys.exc_info()[1], sys.exc_info()[2])) timeErrors = timeErrors + 1 simulationTime = state.getTime()/timeUnit - baseTime self._times.append( simulationTime) kineticEnergy = state.getKineticEnergy()/energyUnit self._kineticEnergies.append( kineticEnergy ) potentialEnergy = state.getPotentialEnergy()/energyUnit self._potentialEnergies.append( potentialEnergy ) totalSteps += self._timePerChunk self._steps.append( totalSteps ) # check for constraint violations constraintViolations = self.checkConstraints( context, totalSteps ) # accumulate volume if( self._calculateVolume ): self.accumulateVolume( context ) chunk += 1 if( self._verbose and ( (numberOfChunks <= 10) or (chunk == printChunk) ) ): totalEnergy = kineticEnergy + potentialEnergy if( totalTime > 0 ): nsPerDay = (86.4*simulationTime)/totalTime else: nsPerDay = 0.0 print "Simulation : %30s %14.7e %14.7e %14.7e %14.7e Cnst=%6d Time=%10.3e ns/day=%10.3e" % (self.getFullTestName(), simulationTime, kineticEnergy, potentialEnergy, totalEnergy, constraintViolations, totalTime, nsPerDay) chunk = 0 sys.stdout.flush() loopStopTime = time.time() localtime = time.strftime("%b %d %Y %H:%M:%S", time.localtime(time.time())) print "Local current time :", localtime # get result self._totalTime = totalTime self._totalLoopTime = loopStopTime - loopStartTime self._totalSteps = totalSteps if( self._totalTime > 0 ): nsPerDay = (86.4*simulationTime)/self._totalTime else: nsPerDay = 0.0 if( self._totalLoopTime > 0 ): nsPerDayLoop = (86.4*simulationTime)/self._totalLoopTime else: nsPerDayLoop = 0.0 print "Loop time delta=%s ns/day=%12.3e ns/day(Loop)=%12.3e" % (str(loopStopTime-loopStartTime), nsPerDay, nsPerDayLoop) sys.stdout.flush() self._result = EnergyConservationResult( self, platformName, self._verbose ) return #============================================================================================= # Create integrator #============================================================================================= def getIntegrator(self, integratorName ): if( integratorName.startswith( 'Variable' ) ): errorTolerance = self._simulationErrorTolerance else: timestep = self._simulationStepSize # VerletIntegrator if( integratorName == 'VerletIntegrator' ): integrator = self._openMM.VerletIntegrator(timestep) if( self._verbose ): print "%30s %30s Timestep=%10.3e" % ( self.getFullTestName(), integratorName, timestep ) return integrator # VariableVerletIntegrator if( integratorName == 'VariableVerletIntegrator' ): integrator = self._openMM.VariableVerletIntegrator(errorTolerance) if( self._verbose ): print "%30s %30s ErrorTolerance=%10.3e" % ( self.getFullTestName(), integratorName, errorTolerance ) return integrator # Stochastic integrators friction = self._simulationFriction temperature = self._simulationTemperature seed = self._randomNumberSeed # Langevin integrator if( integratorName == 'LangevinIntegrator' ): integrator = self._openMM.LangevinIntegrator(temperature,friction,timestep) if( self._verbose ): print "%30s %30s Timestep=%10.3e T=%10.f Friction=%10.3f Seed=%d" % ( self.getFullTestName(), integratorName, timestep, temperature, friction, seed ) # Brownian integrator elif( integratorName == 'BrownianIntegrator' ): integrator = self._openMM.BrownianIntegrator(temperature,friction,timestep) if( self._verbose ): print "%30s %30s Timestep=%10.3e T=%10.f Friction=%10.3f Seed=%d" % ( self.getFullTestName(), integratorName, timestep, temperature, friction, seed ) # VariableLangevin integrator elif( integratorName == 'VariableLangevinIntegrator' ): integrator = self._openMM.VariableLangevinIntegrator(temperature, friction, errorTolerance) if( self._verbose ): print "%30s %30s ErrorTolerance=%10.3e T=%10.f Friction=%10.3f Seed=%d" % ( self.getFullTestName(), integratorName, errorTolerance, temperature, friction, seed ) # woops else: print "\nIntegrator %s not recognized.\n" % integratorName; os.exit(-1) if( seed ): integrator.setRandomNumberSeed( seed ) return integrator #============================================================================================= # Edit system based on test name; e.g., change nonbonded method #============================================================================================= def editSystemBasedOnTestName( self ): self.getBoxDimensions() sys.stdout.flush() minBoxDimension = 1.0e+30 maxBoxDimension = -1.0e+30 for jj in range( len( self._box) ): if( self._box[jj] < minBoxDimension ): minBoxDimension = self._box[jj]; if( self._box[jj] > maxBoxDimension ): maxBoxDimension = self._box[jj]; jj += 1 sys.stdout.flush() offset = 1.5 boxSize = maxBoxDimension + offset box0 = [ boxSize, 0.0, 0.0 ] box1 = [ 0.0, boxSize, 0.0 ] box2 = [ 0.0, 0.0, boxSize ] if( boxSize > 2.0 ): cutoff = 1.0 else: cutoff = 0.4*boxSize if( "NoCutoff" in self._testName ): self.editNonbondedMethod( self._openMM.NonbondedForce.NoCutoff, cutoff, box0, box1, box2 ) if( self._verbose ): print "Setting nonbonded method to NoCutoff\n" elif( "CutoffNonPeriodic" in self._testName ): self.editNonbondedMethod( self._openMM.NonbondedForce.CutoffNonPeriodic, cutoff, box0, box1, box2 ) if( self._verbose ): print "Setting nonbonded method to CutoffNonPeriodic\n" elif( "CutoffPeriodic" in self._testName ): print "In editSystemBasedOnTestName %s CutoffPeriodic\n" % self._testName self.editNonbondedMethod( self._openMM.NonbondedForce.CutoffPeriodic, cutoff, box0, box1, box2 ) if( self._verbose ): print "Setting nonbonded method to CutoffPeriodic\n" elif( "Ewald" in self._testName ): self.editNonbondedMethod( self._openMM.NonbondedForce.Ewald, cutoff, box0, box1, box2 ) if( self._verbose ): print "Setting nonbonded method to Ewald\n" elif( "PME" in self._testName ): self.editNonbondedMethod( self._openMM.NonbondedForce.PME, cutoff, box0, box1, box2 ) if( self._verbose ): print "Setting nonbonded method to PME\n" #self._system.setDefaultPeriodicBoxVectors( box0, box1, box2 ) box = self._system.getDefaultPeriodicBoxVectors( ) if( self._verbose ): print "Box: %s " % (str(box)) nonbondedForce = self._validationUtilities.getForceByName( self._system, "NonbondedForce" ) if( nonbondedForce is not None ): self._validationUtilities.printNonbonded( self.getFullTestName(), nonbondedForce ) sys.stdout.flush() #============================================================================================= # Find nonbonded force and edit nonbonded method #============================================================================================= def editNonbondedMethod( self, nonbondedMethod, cutoff, box0, box1, box2 ): nonbondedForce = self._validationUtilities.getForceByName( self._system, "NonbondedForce" ) if( nonbondedForce is not None ): nonbondedForce.setNonbondedMethod( nonbondedMethod ) nonbondedForce.setCutoffDistance( cutoff ) self._system.setDefaultPeriodicBoxVectors( box0, box1, box2 ) #============================================================================================= # Test class for Verlet integrator test #============================================================================================= class VerletIntegratorTest(EnergyConservationTest): """ Verlet Integrator test inherits from EnergyConservationTest """ def __init__(self, parameterInput, defaultOverrides, verbose=False): """ ARGUMENTS parameterInput (ParameterFileParser class) - ParameterFileParser class containing parameters and defaultOverrides (dictionary of overrides) - verbose (Boolean) - if true, be verbose """ EnergyConservationTest.__init__(self, parameterInput, defaultOverrides, verbose ) self._integratorName = 'VerletIntegrator' self._tolerance = defaultOverrides['VerletTolerance'] if self._verbose: print "\n\n" + self.getFullTestName() #============================================================================================= # Test class for Verlet integrator test w/ constraints converted into bonds #============================================================================================= class VerletIntegratorNoConstraintsTest(EnergyConservationTest): """ Verlet Integrator test inherits from EnergyConservationTest """ def __init__(self, parameterInput, defaultOverrides, verbose=False): """ ARGUMENTS parameterInput (ParameterFileParser class) - ParameterFileParser class containing parameters and defaultOverrides (dictionary of overrides) - verbose (Boolean) - if true, be verbose """ EnergyConservationTest.__init__(self, parameterInput, defaultOverrides, verbose ) # build system self._system = self._validationUtilities.copySystem( parameterInput, 1 ) self.editSystemBasedOnTestName() self._dof = self._validationUtilities.getDegreesOfFreedom( parameterInput ) self._integratorName = 'VerletIntegrator' self._tolerance = defaultOverrides['VerletTolerance'] if self._verbose: print "\n\n" + self.getFullTestName() #============================================================================================= # Test class for VariableVerlet integrator test #============================================================================================= class VariableVerletIntegratorTest(EnergyConservationTest): """ VariableVerlet Integrator test inherits from EnergyConservationTest """ def __init__(self, parameterInput, defaultOverrides, verbose=False): """ ARGUMENTS parameterInput (ParameterFileParser class) - ParameterFileParser class containing parameters and defaultOverrides (dictionary of overrides) - verbose (Boolean) - if true, be verbose """ EnergyConservationTest.__init__(self, parameterInput, defaultOverrides, verbose ) # self._systemName = parameterInput.getName().replace( '.txt', '' ) # self._testName = defaultOverrides['TestName'] # # # build system # # self._system = self._validationUtilities.copySystem( parameterInput ) # self.editSystemBasedOnTestName() # self._dof = self._validationUtilities.getDegreesOfFreedom( parameterInput ) self._integratorName = 'VariableVerletIntegrator' self._tolerance = defaultOverrides['VerletTolerance'] if self._verbose: print "\n\n" + self.getFullTestName() # self._positions = parameterInput.getPositions() #============================================================================================= # Base class for Langevin energy conservation tests #============================================================================================= class LangevinEnergyConservationTest(EnergyConservationTest): """Base class for Langevin energy conservation tests Centralizes common functionality across Langevin energy conservation tests """ def __init__(self, parameterInput, defaultOverrides, verbose=False): """ ARGUMENTS defaultOverrides (dictionary of overrides) - verbose (Boolean) - if true, be verbose """ EnergyConservationTest.__init__(self, parameterInput, defaultOverrides, verbose ) # set parameters if( 'SimulationTemperature' in defaultOverrides ): self._simulationTemperature = defaultOverrides['SimulationTemperature'] else: self._simulationTemperature = 300 if( 'SimulationFriction' in defaultOverrides ): self._simulationFriction = defaultOverrides['SimulationFriction'] else: self._simulationFriction = 91 if( 'SimulationRandomNumberSeed' in defaultOverrides ): self._randomNumberSeed = defaultOverrides['SimulationRandomNumberSeed'] else: self._randomNumberSeed = 1934 self._tolerance = defaultOverrides['LangevinTolerance'] #============================================================================================= # Test class for Langevin integrator test #============================================================================================= class LangevinIntegratorTest(LangevinEnergyConservationTest): """ Langevin Integrator test inherits from LangevinEnergyConservationTest """ def __init__(self, parameterInput, defaultOverrides, verbose=False): """ ARGUMENTS parameterInput (ParameterFileParser class) - ParameterFileParser class containing parameters and coordinates defaultOverrides (dictionary of overrides) - verbose (Boolean) - if true, be verbose """ LangevinEnergyConservationTest.__init__(self, parameterInput, defaultOverrides, verbose ) self._integratorName = 'LangevinIntegrator' if self._verbose: print "\n\n" + self.getFullTestName() # self._positions = parameterInput.getPositions() #============================================================================================= # Test class for VariableLangevin integrator test #============================================================================================= class VariableLangevinIntegratorTest(LangevinEnergyConservationTest): """ VariableLangevin Integrator test inherits from LangevinEnergyConservationTest """ def __init__(self, parameterInput, defaultOverrides, verbose=False): """ ARGUMENTS parameterInput (ParameterFileParser class) - ParameterFileParser class containing parameters and defaultOverrides (dictionary of overrides) - verbose (Boolean) - if true, be verbose """ LangevinEnergyConservationTest.__init__(self, parameterInput, defaultOverrides, verbose ) self._integratorName = 'VariableLangevinIntegrator' if self._verbose: print "\n\n" + self.getFullTestName() # self._positions = parameterInput.getPositions() #============================================================================================= # Test class for VariableLangevin integrator test #============================================================================================= class MonteCarloBarostatLangevinIntegratorTest(LangevinEnergyConservationTest): """ MonteCarloBarostatLangevinIntegratorTest Integrator test inherits from VariableLangevinEnergyConservationTest """ def __init__(self, parameterInput, defaultOverrides, verbose=False): """ ARGUMENTS parameterInput (ParameterFileParser class) - ParameterFileParser class containing parameters and defaultOverrides (dictionary of overrides) - verbose (Boolean) - if true, be verbose """ LangevinEnergyConservationTest.__init__(self, parameterInput, defaultOverrides, verbose ) self._integratorName = 'LangevinIntegrator' pressure = 3.0 if( 'Pressure' in defaultOverrides ): pressure = float( defaultOverrides['Pressure'] ) if( 'RandomNumberSeed' in defaultOverrides ): randomNumberSeed = int( defaultOverrides['RandomNumberSeed'] ) else: randomNumberSeed = -1 monteCarloBarostat = self._openMM.MonteCarloBarostat(pressure, self._simulationTemperature ) if( randomNumberSeed > 0 ): monteCarloBarostat.setRandomNumberSeed( randomNumberSeed ) self._system.addForce(monteCarloBarostat) self._calculateVolume = 1 if self._verbose: print "\n\n" + self.getFullTestName() print "Pressure %s" % (str(monteCarloBarostat.getDefaultPressure())) print "Temperature %s" % (str(monteCarloBarostat.getTemperature())) print "Frequency %d" % (monteCarloBarostat.getFrequency()) print "Random number seed %12d" % (monteCarloBarostat.getRandomNumberSeed())