#!/bin/env python #============================================================================================= # GLOBAL IMPORTS #============================================================================================= import os import os.path import sys import copy import re import math import time import simtk.openmm from ValidationUtilities import ValidationUtilities #============================================================================================= # Results class for force-energy gradient test #============================================================================================= class ForceEnergyGradientResult(object): """ ForceEnergyGradientResult stores results for force-energy gradient test """ def __init__( self, test, platform, forceNorm, potentialEnergy, perturbedPotentialEnergy, verbose=False ): """ Generate results ARGUMENTS test (ForceEnergyGradientTest) - ForceEnergyGradientTest reference platform (string) - platform name forceNorm (float) - force norm of unperturbed system potentialEnergy (float) - PE of unperturbed system perturbedPotentialEnergy (float) - PE of perturbed system verbose (Boolean) - if true, print logging info """ self._platform = platform self._openMM = simtk.openmm self._verbose = verbose self._validationUtilities = ValidationUtilities(verbose) self._forceStepDelta = test.getForceStepDelta() self._systemName = test.getSystemName() self._testName = test.getTestName() self._tolerance = test.getTolerance() self._forceNorm = forceNorm self._passed = -1 # compute finite difference energy derivative self._deltaPotentialEnergy = abs( potentialEnergy - perturbedPotentialEnergy )/test.getForceStepDelta() if self._verbose: print "PEs %14.7e %14.7e deltaE=%14.7e forceNorm=%14.7e diff=%14.7e" % ( potentialEnergy, perturbedPotentialEnergy, self._deltaPotentialEnergy, forceNorm, abs(self._deltaPotentialEnergy-forceNorm)/forceNorm ) sys.stdout.flush() return #============================================================================================= # Return tolerance #============================================================================================= def getTolerance(self): return self._tolerance #============================================================================================= # Get delta in potential energies #============================================================================================= def getDeltaPotentialEnergy(self): return self._deltaPotentialEnergy #============================================================================================= # Get relative delta between energy gradient and force norm #============================================================================================= def getRelativeDelta(self): return self._relDelta #============================================================================================= # Get system name #============================================================================================= def getSystemName(self): return self._systemName #============================================================================================= # Get platform name #============================================================================================= def getPlatformName(self): return self._platform #============================================================================================= # Get test name #============================================================================================= def getTestName(self): return self._testName #============================================================================================= # Return force step delta used in perturbing coordinates #============================================================================================= def getForceStepDelta(self): return self._forceStepDelta #============================================================================================= # Return 1 if test passsed, else 0 #============================================================================================= def testPassed(self): sys.stdout.flush() if( self._passed > -1 ): return self._passed self._passed = 1 self._delta = abs( self._deltaPotentialEnergy - self._forceNorm ) self._relDelta = self._delta/self._forceNorm if( self._relDelta > self._tolerance ): if( self._verbose ): print "Test failed: relative delta %14.7e is greater than tolerance=%10.3e" % (self._relDelta, self._tolerance ) self._passed = 0 if( self._verbose ): sys.stdout.flush() return self._passed #============================================================================================= # Base class for force-energy gradient tests #============================================================================================= class ForceEnergyGradientTest(object): """Base class for force-energy gradient tests Centralizes common functionality across force-energy gradient tests """ def __init__(self, defaultOverrides, verbose=False): """ ARGUMENTS defaultOverrides (dictionary of overrides) - verbose (Boolean) - if true, be verbose """ # tolerance if( 'Tolerance' in defaultOverrides ): self._tolerance = defaultOverrides['Tolerance'] else: self._tolerance = 4.0e-02 # delta used in perturbing positions in direction of force if( 'ForceStepDelta' in defaultOverrides ): self._forceStepDelta = defaultOverrides['ForceStepDelta'] else: self._forceStepDelta = 1.0e-03 self._isActive = 1 self._openMM = simtk.openmm self._verbose = verbose self._validationUtilities = ValidationUtilities(verbose) # try a range of deltas in direction of force step # in case default value is not successful self._forceStepDeltaHash = {} self._forceStepDeltaHash[1.0e-02] = 0 self._forceStepDeltaHash[5.0e-03] = 0 self._forceStepDeltaHash[1.0e-03] = 0 self._forceStepDeltaHash[5.0e-04] = 0 self._forceStepDeltaHash[1.0e-04] = 0 self._forceStepDeltaHash[5.0e-05] = 0 self._forceStepDeltaHash[1.0e-05] = 0 # device id if( 'DeviceId' in defaultOverrides ): self._deviceId = defaultOverrides['DeviceId'] else: self._deviceId = -1 #============================================================================================= # Return 1 if test is active #============================================================================================= def isActive(self): return self._isActive #============================================================================================= # Return tolerance #============================================================================================= def getTolerance(self): return self._tolerance #============================================================================================= # Return force step delta used in perturbing coordinates #============================================================================================= def getForceStepDelta(self): return self._forceStepDelta #============================================================================================= # Get test name #============================================================================================= def getFullTestName(self): return self._systemName + '_' + self._testName #============================================================================================= # Get system name #============================================================================================= def getSystemName(self): return self._systemName #============================================================================================= # Get system name #============================================================================================= def getTestName(self): return self._testName #============================================================================================= # Get result #============================================================================================= def getResult(self): return self._result 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 def setNonbondedBox( self ): self.getBoxDimensions() minBoxDimension = 1.0e+30 for jj in range( 3 ): if( self._box[jj] < minBoxDimension ): minBoxDimension = self._box[jj]; jj += 1 offset = 0.5 box0 = [ float(self._box[0] + offset), 0.0, 0.0 ] box1 = [ 0.0, float(self._box[1] + offset), 0.0 ] box2 = [ 0.0, 0.0, float(self._box[2] + offset) ] if( minBoxDimension > 0.0 ): if( minBoxDimension > 2.0 ): cutoff = 1.0 else: cutoff = 0.4*minBoxDimension self._cutoff = cutoff self._box0 = box0 self._box1 = box1 self._box2 = box2 else: self._cutoff = None self._box0 = None #============================================================================================= # 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 #============================================================================================= # Get energy for perturbed positions #============================================================================================= def getEnergyForPerturbedPositions(self, forces, context, step ): energyUnit = self._validationUtilities.getEnergyUnit( ) forceUnit = self._validationUtilities.getForceUnit( ) # perturb coordinates ii = 0 self._perturbedPositions = [] while( ii < len( self._positions ) ): coordinates = self._positions[ii] force = forces[ii] / forceUnit xCoord = coordinates[0] - step*force[0] yCoord = coordinates[1] - step*force[1] zCoord = coordinates[2] - step*force[2] coords = [] coords.append( xCoord ) coords.append( yCoord ) coords.append( zCoord ) self._perturbedPositions.append( coords ) ii += 1 context.setPositions( self._perturbedPositions ) state = context.getState(getEnergy=True, getForces=True) # check if test passes w/ current forceStepDelta value perturbedPotentialEnergy = state.getPotentialEnergy() / energyUnit return perturbedPotentialEnergy #============================================================================================= # Run test #============================================================================================= def runTest(self, platformName): if self._verbose: print "runTest: %s Particles=%d NumForces=%d Platform=%s delta=%10.3e" % (self.getFullTestName(), self._system.getNumParticles(), self._system.getNumForces(), platformName, self._forceStepDelta ) sys.stdout.flush() timestep = 0.01 integrator = self._openMM.VerletIntegrator(timestep) platform = self._openMM.Platform.getPlatformByName( platformName ) if( platformName is not 'Reference' and self._deviceId > -1 ): self._validationUtilities.setDeviceId( platform, self._deviceId ) context = self._openMM.Context(self._system, integrator, platform ) context.setPositions( self._positions ) state = context.getState(getEnergy=True, getForces=True) if( self._verbose ): self._validationUtilities.printPlatformProperties( context ) # get units from ValidationUtilities energyUnit = self._validationUtilities.getEnergyUnit( ) forceUnit = self._validationUtilities.getForceUnit( ) # get force norm and perturb coordinates potentialEnergy = state.getPotentialEnergy() / energyUnit forces = state.getForces() forceNorm = 0.0 ii = 0 while( ii < len( forces ) ): #print "runTest: forces=%s" % (str( forces[ii] ) ) force = forces[ii] / forceUnit forceNorm += force[0]*force[0] + force[1]*force[1] + force[2]*force[2] ii += 1 forceNorm = math.sqrt( forceNorm ) if( forceNorm <= 0.0 ): print "Force norm is zero for %s!" % self.getFullTestName() self._result = ForceEnergyGradientResult( self, platformName, state1, state2, self._verbose ) return elif( self._verbose ): print "Force norm is %14.7e" % forceNorm # loop over range of values of forceStepDelta, starting w/ the initial value # if the initial value is not successful, try another value from # self._forceStepDeltaHash. Continue until the test passes or all values in # self._forceStepDeltaHash have been tried minRelDelta = 1.0e+30 minForceStepDelta = 0 passed = 0 self._forceStepDeltaHash[self._forceStepDelta] = 0 for forceStepDelta, value in self._forceStepDeltaHash.iteritems(): self._forceStepDelta = forceStepDelta step = self._forceStepDelta/forceNorm perturbedPotentialEnergy = self.getEnergyForPerturbedPositions(forces, context, step ) self._result = ForceEnergyGradientResult( self, platformName, forceNorm, potentialEnergy, perturbedPotentialEnergy, self._verbose ) self._forceStepDeltaHash[self._forceStepDelta] = 1 if( self._result.testPassed() ): passed = 1 if( self._verbose ): print "%s passed with forceStepDelta=%10.3e\n" % (self.getFullTestName(), self._forceStepDelta) else: if( self._verbose ): print "%s did not pass with forceStepDelta=%10.3e\n" % (self.getFullTestName(), self._forceStepDelta) if( self._result.getRelativeDelta() < minRelDelta ): minRelDelta = self._result.getRelativeDelta() minForceStepDelta = self._forceStepDelta if( passed ): if( self._verbose ): print "%s min forceStepDelta=%10.3e min relatvive delta=%10.3e\n" % (self.getFullTestName(), self._forceStepDelta, minRelDelta) self._forceStepDelta = minForceStepDelta step = self._forceStepDelta/forceNorm perturbedPotentialEnergy = self.getEnergyForPerturbedPositions( forces, context, step ) self._result = ForceEnergyGradientResult( self, platformName, forceNorm, potentialEnergy, perturbedPotentialEnergy, self._verbose ) else: if( self._verbose ): print "%s did not pass for any deltas\n" % (self.getFullTestName() ) # notDone = 1 # while( notDone ): # # step = self._forceStepDelta/forceNorm # # # perturb coordinates # # while( ii < len( self._positions ) ): # coordinates = self._positions[ii] # force = forces[ii] / forceUnit # xCoord = coordinates[0] - step*force[0] # yCoord = coordinates[1] - step*force[1] # zCoord = coordinates[2] - step*force[2] # self._perturbedPositions.append( Vec3( (xCoord, yCoord, zCoord) ) ) # ii += 1 # # context.setPositions( self._perturbedPositions ) # state = context.getState(getEnergy=True, getForces=True) # # # check if test passes w/ current forceStepDelta value # # perturbedPotentialEnergy = state.getPotentialEnergy() / energyUnit # self._result = ForceEnergyGradientResult( self, platformName, forceNorm, potentialEnergy, # perturbedPotentialEnergy, self._verbose ) # self._forceStepDeltaHash[self._forceStepDelta] = 1 # if( self._result.testPassed() ): # notDone = 0 # if( self._verbose ): # print "%s passed with forceStepDelta=%10.3e\n" % (self.getFullTestName(), self._forceStepDelta) # else: # if( self._verbose ): # print "%s did not pass with forceStepDelta=%10.3e\n" % (self.getFullTestName(), self._forceStepDelta) # # # get next forceStepDelta; if all have been tried, then # # set notDone == 0, and exit loop/method # # hit = 0 # for forceStepDelta, value in self._forceStepDeltaHash.iteritems(): # if( value == 0 and hit == 0 ): # hit = 1 # self._forceStepDelta = forceStepDelta # if( hit == 0 ): # notDone = 0 # # if( self._verbose ): # sys.stdout.flush() return #============================================================================================= # Test class for force-energy gradient test for HarmonicBondForce #============================================================================================= class HarmonicBondForceEnergyGradientTest(ForceEnergyGradientTest): """ HarmonicBondForceEnergyGradientTest create system for testing HarmonicBondForce inherits from ForceEnergyGradientTest """ 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 """ ForceEnergyGradientTest.__init__(self, defaultOverrides, verbose ) self._systemName = parameterInput.getName().replace( '.txt', '' ) self._testName = 'HarmonicBondForceEnergyGradientTest' if( parameterInput.getNumberOfHarmonicBonds() == 0 ): if verbose: print "No harmonic bonds for %s" % parameterInput.getName() return harmonicBondForce = parameterInput.getHarmonicBondForce() # build system self._system = self._openMM.System() self._validationUtilities.loadSystemMasses( parameterInput.getMasses(), self._system ) testHarmonicBondForce = self._validationUtilities.copyHarmonicBondForce( harmonicBondForce ) self._system.addForce(testHarmonicBondForce) if self._verbose: print "\n\n" + self.getFullTestName() print "Added %d harmonic bonds" % testHarmonicBondForce.getNumBonds() ii = 0 lengthUnit = self._validationUtilities.getLengthUnit( ) forceConstantUnit = self._validationUtilities.getEnergyUnit( ) / lengthUnit**2 maxPrint = 10 while( ii < testHarmonicBondForce.getNumBonds() ): args = testHarmonicBondForce.getBondParameters(ii) print "HarmonicBond: %5d [%5d %5d] [%14.7e %14.7e]" % (ii, args[0], args[1], args[2]/lengthUnit, args[3]/forceConstantUnit ) ii += 1 if( ii == maxPrint and testHarmonicBondForce.getNumBonds() > maxPrint ):ii = testHarmonicBondForce.getNumBonds() - maxPrint sys.stdout.flush() self._positions = parameterInput.getPositions() #============================================================================================= # Test class for force-energy gradient test for HarmonicAngleForce #============================================================================================= class HarmonicAngleForceEnergyGradientTest(ForceEnergyGradientTest): """ HarmonicAngleForceEnergyGradientTest create system for testing HarmonicAngleForce inherits from ForceEnergyGradientTest """ 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 """ ForceEnergyGradientTest.__init__(self, defaultOverrides, verbose ) self._systemName = parameterInput.getName().replace( '.txt', '' ) self._testName = 'HarmonicAngleForceEnergyGradientTest' print "No. of harmonic angles for %s=%d QQQQ" % (parameterInput.getName(), parameterInput.getNumberOfHarmonicAngles() ) if( parameterInput.getNumberOfHarmonicAngles() == 0 ): if verbose: print "No harmonic angles for %s" % parameterInput.getName() return harmonicAngleForce = parameterInput.getHarmonicAngleForce() # build system self._system = self._openMM.System() self._validationUtilities.loadSystemMasses( parameterInput.getMasses(), self._system ) testHarmonicAngleForce = self._validationUtilities.copyHarmonicAngleForce( harmonicAngleForce ) self._system.addForce(testHarmonicAngleForce) if self._verbose: print "\n\n" + self.getFullTestName() print "Added %d harmonic angles" % testHarmonicAngleForce.getNumAngles() ii = 0 lengthUnit = self._validationUtilities.getRadianUnit( ) forceConstantUnit = self._validationUtilities.getEnergyUnit( ) / lengthUnit**2 maxPrint = 10 while( ii < testHarmonicAngleForce.getNumAngles() ): args = testHarmonicAngleForce.getAngleParameters(ii) print "HarmonicAngle: %5d [%5d %5d %5d] [%14.7e %14.7e]" % (ii, args[0], args[1], args[2], args[3]/lengthUnit, args[4]/forceConstantUnit ) ii += 1 if( ii == maxPrint and testHarmonicAngleForce.getNumAngles() > maxPrint ):ii = testHarmonicAngleForce.getNumAngles() - maxPrint sys.stdout.flush() self._positions = parameterInput.getPositions() #============================================================================================= # Test class for force-energy gradient test for PeriodicTorsionForce #============================================================================================= class PeriodicTorsionForceEnergyGradientTest(ForceEnergyGradientTest): """ PeriodicTorsionForceEnergyGradientTest create system for testing PeriodicTorsionForce inherits from ForceEnergyGradientTest """ 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 """ ForceEnergyGradientTest.__init__(self, defaultOverrides, verbose ) self._systemName = parameterInput.getName().replace( '.txt', '' ) self._testName = 'PeriodicTorsionForceEnergyGradientTest' if( parameterInput.getNumberOfPeriodicTorsions() == 0 ): if verbose: print "No periodic torsions for %s" % parameterInput.getName() return periodicTorsionForce = parameterInput.getPeriodicTorsionForce() # build system self._system = self._openMM.System() self._validationUtilities.loadSystemMasses( parameterInput.getMasses(), self._system ) testPeriodicTorsionForce = self._validationUtilities.copyPeriodicTorsionForce( periodicTorsionForce ) self._system.addForce(testPeriodicTorsionForce) if self._verbose: print "\n\n" + self.getFullTestName() print "Added %d periodic torsions " % testPeriodicTorsionForce.getNumTorsions() ii = 0 lengthUnit = self._validationUtilities.getRadianUnit( ) energyUnit = self._validationUtilities.getEnergyUnit( ) maxPrint = 10 while( ii < testPeriodicTorsionForce.getNumTorsions() ): args = testPeriodicTorsionForce.getTorsionParameters(ii) print "PeriodicTorsion: %5d [%5d %5d %5d %5d] %3d [%14.7e %14.7e]" % (ii, args[0], args[1], args[2], args[3], args[4], args[5]/lengthUnit, args[6]/energyUnit) ii += 1 if( ii == maxPrint and testPeriodicTorsionForce.getNumTorsions() > maxPrint ):ii = testPeriodicTorsionForce.getNumTorsions() - maxPrint sys.stdout.flush() self._positions = parameterInput.getPositions() #============================================================================================= # Test class for force-energy gradient test for RBTorsionForce #============================================================================================= class RBTorsionForceEnergyGradientTest(ForceEnergyGradientTest): """ RBTorsionForceEnergyGradientTest create system for testing RBTorsionForce inherits from ForceEnergyGradientTest """ 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 """ ForceEnergyGradientTest.__init__(self, defaultOverrides, verbose ) self._systemName = parameterInput.getName().replace( '.txt', '' ) self._testName = 'RBTorsionForceEnergyGradientTest' if( parameterInput.getNumberOfRBTorsions() == 0 ): if verbose: print "No rb torsions for %s" % parameterInput.getName() return rbTorsionForce = parameterInput.getRBTorsionForce() # build system self._system = self._openMM.System() self._validationUtilities.loadSystemMasses( parameterInput.getMasses(), self._system ) testRBTorsionForce = self._validationUtilities.copyRBTorsionForce( rbTorsionForce ) self._system.addForce(testRBTorsionForce) if self._verbose: print "\n\n" + self.getFullTestName() print "Added %d rb torsions " % testRBTorsionForce.getNumTorsions() ii = 0 maxPrint = 10 while( ii < testRBTorsionForce.getNumTorsions() ): args = testRBTorsionForce.getTorsionParameters(ii) print "RBTorsion: %5d [%5d %5d %5d %5d] [%14.7e %14.7e %14.7e %14.7e %14.7e %14.7e]" % (ii, args[0], args[1], args[2], args[3], args[4], args[5], args[6], args[7], args[8], args[9]) ii += 1 if( ii == maxPrint and testRBTorsionForce.getNumTorsions() > maxPrint ):ii = testRBTorsionForce.getNumTorsions() - maxPrint sys.stdout.flush() self._positions = parameterInput.getPositions() #============================================================================================= # Base class for force-energy gradient tests for NonbondedForce #============================================================================================= class NonbondedForceEnergyGradientBaseTest(ForceEnergyGradientTest): """ NonbondedForceTest create system for testing NonbondedForce inherits from ForceEnergyGradientTest """ 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 """ ForceEnergyGradientTest.__init__(self, defaultOverrides, verbose ) self._systemName = parameterInput.getName().replace( '.txt', '' ) if( parameterInput.getNumberOfNonbondeds() == 0 ): if verbose: print "No nonboned particles %s" % parameterInput.getName() return nonbondedForce = parameterInput.getNonbondedForce() # build system self._system = self._openMM.System() self._validationUtilities.loadSystemMasses( parameterInput.getMasses(), self._system ) self._testNonbondedForce = self._validationUtilities.copyNonbondedForce( nonbondedForce ) self._positions = parameterInput.getPositions() self.setNonbondedBox( ) if( self._cutoff is not None ): self._testNonbondedForce.setCutoffDistance( self._cutoff ) if( self._box0 is not None ): self._system.setDefaultPeriodicBoxVectors( self._box0, self._box1, self._box2) self._system.addForce(self._testNonbondedForce) #============================================================================================= # Test class for force-energy gradient test for NonbondedForce w/ no cutoff #============================================================================================= class NonbondedForceEnergyGradientNoCutoffTest(NonbondedForceEnergyGradientBaseTest): """ NonbondedForceTest create system for testing NonbondedForce w/ no cutoff inherits from NonbondedForceEnergyGradientBaseTest """ 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 """ NonbondedForceEnergyGradientBaseTest.__init__(self, parameterInput, defaultOverrides, verbose ) self._testName = 'NonbondedForceEnergyGradientNoCutoffTest' self._testNonbondedForce.setNonbondedMethod( self._openMM.NonbondedForce.NoCutoff ) if self._verbose: self._validationUtilities.printNonbonded( self._testName, self._testNonbondedForce ) #============================================================================================= # Test class for force-energy gradient test for NonbondedForce w/ cutoff #============================================================================================= class NonbondedForceEnergyGradientCutoffNonPeriodicTest(NonbondedForceEnergyGradientBaseTest): """ NonbondedForceTest create system for testing NonbondedForce w/ cutoff and non-periodic boundary conditions inherits from NonbondedForceEnergyGradientBaseTest """ 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 """ NonbondedForceEnergyGradientBaseTest.__init__(self, parameterInput, defaultOverrides, verbose ) self._testName = 'NonbondedForceEnergyGradientCutoffNonPeriodicTest' self._testNonbondedForce.setNonbondedMethod( self._openMM.NonbondedForce.CutoffNonPeriodic ) if self._verbose: self._validationUtilities.printNonbonded( self._testName, self._testNonbondedForce ) #============================================================================================= # Test class for force-energy gradient test for NonbondedForce w/ cutoff # and periodic boundary conditions #============================================================================================= class NonbondedForceEnergyGradientCutoffPeriodicTest(NonbondedForceEnergyGradientBaseTest): """ NonbondedForceTest create system for testing NonbondedForce w/ cutoff and periodic boundary conditions inherits from NonbondedForceEnergyGradientBaseTest """ 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 """ NonbondedForceEnergyGradientBaseTest.__init__(self, parameterInput, defaultOverrides, verbose ) self._testName = 'NonbondedForceEnergyGradientCutoffPeriodicTest' self._testNonbondedForce.setNonbondedMethod( self._openMM.NonbondedForce.CutoffPeriodic ) if self._verbose: self._validationUtilities.printNonbonded( self._testName, self._testNonbondedForce ) #============================================================================================= # Test class for force-energy gradient test for NonbondedForce using Ewald #============================================================================================= class NonbondedForceEnergyGradientEwaldTest(NonbondedForceEnergyGradientBaseTest): """ NonbondedForceTest create system for testing NonbondedForce using Ewald inherits from NonbondedForceEnergyGradientBaseTest """ 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 """ NonbondedForceEnergyGradientBaseTest.__init__(self, parameterInput, defaultOverrides, verbose ) self._testName = 'NonbondedForceEnergyGradientEwaldTest' self._testNonbondedForce.setNonbondedMethod( self._openMM.NonbondedForce.Ewald ) if self._verbose: self._validationUtilities.printNonbonded( self._testName, self._testNonbondedForce ) #============================================================================================= # Test class for force-energy gradient test for NonbondedForce w/ PME #============================================================================================= class NonbondedForceEnergyGradientPMETest(NonbondedForceEnergyGradientBaseTest): """ NonbondedForceTest create system for testing NonbondedForce using PME inherits from NonbondedForceEnergyGradientBaseTest """ 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 """ NonbondedForceEnergyGradientBaseTest.__init__(self, parameterInput, defaultOverrides, verbose ) self._testName = 'NonbondedForceEnergyGradientPMETest' self._testNonbondedForce.setNonbondedMethod( self._openMM.NonbondedForce.PME ) if self._verbose: self._validationUtilities.printNonbonded( self._testName, self._testNonbondedForce ) #============================================================================================= # Base class for force-energy gradient tests for GbsaObcForce #============================================================================================= class GbsaObcForceEnergyGradientBaseTest(ForceEnergyGradientTest): """ GbsaObcForceTest create system for testing GbsaObcForce inherits from ForceEnergyGradientTest """ 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 """ ForceEnergyGradientTest.__init__(self, defaultOverrides, verbose ) self._systemName = parameterInput.getName().replace( '.txt', '' ) # include nonbonded and OBC forces if( parameterInput.getNumberOfGbsaObcs() == 0 ): if verbose: print "No Gbsa Obc particles %s" % parameterInput.getName() self._isActive = 0 return if( parameterInput.getNumberOfNonbondeds() == 0 ): if verbose: print "No nonboned particles %s" % parameterInput.getName() self._isActive = 0 return # build system self._system = self._openMM.System() self._validationUtilities.loadSystemMasses( parameterInput.getMasses(), self._system ) nonbondedForce = parameterInput.getNonbondedForce() gbsaObcForce = parameterInput.getGbsaObcForce() self._positions = parameterInput.getPositions() self.setNonbondedBox( ) print "_cutoff XX %12.3f" % self._cutoff sys.stdout.flush() if( self._box0 is not None ): self._system.setDefaultPeriodicBoxVectors( self._box0, self._box1, self._box2) self._testNonbondedForce = self._validationUtilities.copyNonbondedForce( nonbondedForce ) if( self._cutoff is not None ): self._testNonbondedForce.setCutoffDistance( self._cutoff ) self._system.addForce(self._testNonbondedForce) self._testGbsaObcForce = self._validationUtilities.copyGbsaObcForce( gbsaObcForce ) self._system.addForce(self._testGbsaObcForce) self._testGbsaObcForce.setCutoffDistance( self._testNonbondedForce.getCutoffDistance() ) #============================================================================================= # Print GBSA Obc info #============================================================================================= def printGbsaObc(self): dimensionlessUnit = self._validationUtilities.getDimensionlessUnit( ) chargeUnit = self._validationUtilities.getChargeUnit( ) lengthUnit = self._validationUtilities.getLengthUnit( ) print "\n\n" + self.getFullTestName() print "GbsaObc particles %14d" % self._testGbsaObcForce.getNumParticles() print "Gbsa Obc method %14s" % self._testGbsaObcForce.getNonbondedMethod() print "Solute dielectric %14.7f" % (self._testGbsaObcForce.getSoluteDielectric()/dimensionlessUnit) print "Solvent dielectric %14.7f" % (self._testGbsaObcForce.getSolventDielectric()/dimensionlessUnit) print "Cutoff distance %14.7f" % (self._testNonbondedForce.getCutoffDistance()/lengthUnit) ii = 0 maxPrint = 10 while( ii < self._testGbsaObcForce.getNumParticles() ): args = self._testGbsaObcForce.getParticleParameters(ii) print "GbsaObc: %5d [%14.7f %14.7f %14.7f]" % (ii, args[0]/chargeUnit, args[1]/lengthUnit, args[2]/dimensionlessUnit) ii += 1 if( ii == maxPrint and self._testGbsaObcForce.getNumParticles() > maxPrint ):ii = self._testGbsaObcForce.getNumParticles() - maxPrint #============================================================================================= # Test class for force-energy gradient test for GbsaObcForce w/ no cutoff #============================================================================================= class GbsaObcForceEnergyGradientNoCutoffTest(GbsaObcForceEnergyGradientBaseTest): """ GbsaObcForceEnergyGradientNoCutoffTest create system for testing GbsaObcForce w/ no cutoff inherits from GbsaObcForceEnergyGradientBaseTest """ 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 """ GbsaObcForceEnergyGradientBaseTest.__init__(self, parameterInput, defaultOverrides, verbose ) if( self.isActive() == 0 ): return self._testName = 'GbsaObcForceEnergyGradientNoCutoffTest' self._testGbsaObcForce.setNonbondedMethod( self._openMM.NonbondedForce.NoCutoff ) self._testNonbondedForce.setNonbondedMethod( self._openMM.NonbondedForce.NoCutoff ) if self._verbose: self.printGbsaObc() #============================================================================================= # Test class for force-energy gradient test for GbsaObcForce w/ cutoff # and non-periodic boundary conditions #============================================================================================= class GbsaObcForceEnergyGradientCutoffNonPeriodicTest(GbsaObcForceEnergyGradientBaseTest): """ GbsaObcForceTest create system for testing GbsaObcForce w/ cutoff and non-periodic boundary conditions inherits from GbsaObcForceEnergyGradientBaseTest """ 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 """ GbsaObcForceEnergyGradientBaseTest.__init__(self, parameterInput, defaultOverrides, verbose ) if( self.isActive() == 0 ): return self._testName = 'GbsaObcForceEnergyGradientCutoffNonPeriodicTest' self._testNonbondedForce.setNonbondedMethod( self._openMM.NonbondedForce.CutoffNonPeriodic ) self._testGbsaObcForce.setNonbondedMethod( self._openMM.NonbondedForce.CutoffNonPeriodic ) if self._verbose: self.printGbsaObc() #============================================================================================= # Test class for force-energy gradient test for GbsaObcForce w/ cutoff # and periodic boundary conditions #============================================================================================= class GbsaObcForceEnergyGradientCutoffPeriodicTest(GbsaObcForceEnergyGradientBaseTest): """ GbsaObcForceTest create system for testing GbsaObcForce w/ cutoff and periodic boundary conditions inherits from GbsaObcForceEnergyGradientBaseTest """ 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 """ GbsaObcForceEnergyGradientBaseTest.__init__(self, parameterInput, defaultOverrides, verbose ) if( self.isActive() == 0 ): return self._testName = 'GbsaObcForceEnergyGradientCutoffPeriodicTest' self._testNonbondedForce.setNonbondedMethod( self._openMM.NonbondedForce.CutoffPeriodic ) self._testGbsaObcForce.setNonbondedMethod( self._openMM.NonbondedForce.CutoffPeriodic ) if self._verbose: self.printGbsaObc() #============================================================================================= # Base class for force-energy gradient tests for GbviForce #============================================================================================= class GbviForceEnergyGradientBaseTest(ForceEnergyGradientTest): """ GbviForceTest create system for testing GbviForce inherits from ForceEnergyGradientTest """ 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 """ ForceEnergyGradientTest.__init__(self, defaultOverrides, verbose ) self._systemName = parameterInput.getName().replace( '.txt', '' ) # include nonbonded and GB/VI forces if( parameterInput.getNumberOfGbvis() == 0 ): if verbose: print "No Gbvi particles %s" % parameterInput.getName() self._isActive = 0 return if( parameterInput.getNumberOfNonbondeds() == 0 ): if verbose: print "No nonboned particles %s" % parameterInput.getName() self._isActive = 0 return # build system self._system = self._openMM.System() self._validationUtilities.loadSystemMasses( parameterInput.getMasses(), self._system ) nonbondedForce = parameterInput.getNonbondedForce() gbviForce = parameterInput.getGbviForce() self._positions = parameterInput.getPositions() self.setNonbondedBox( ) #print "_cutoff %12.3f" % self._cutoff #sys.stdout.flush() if( self._box0 is not None ): self._system.setDefaultPeriodicBoxVectors( self._box0, self._box1, self._box2) self._testNonbondedForce = self._validationUtilities.copyNonbondedForce( nonbondedForce ) if( self._cutoff is not None ): self._testNonbondedForce.setCutoffDistance( self._cutoff ) self._system.addForce(self._testNonbondedForce) self._testGbviForce = self._validationUtilities.copyGbviForce( gbviForce ) self._system.addForce(self._testGbviForce) self._testGbviForce.setCutoffDistance( self._testNonbondedForce.getCutoffDistance() ) #============================================================================================= # Print GBSA Obc info #============================================================================================= def printGbvi(self): dimensionlessUnit = self._validationUtilities.getDimensionlessUnit( ) energyUnit = self._validationUtilities.getEnergyUnit( ) chargeUnit = self._validationUtilities.getChargeUnit( ) lengthUnit = self._validationUtilities.getLengthUnit( ) print "\n\n" + self.getFullTestName() print "Gbvi particles %14d" % self._testGbviForce.getNumParticles() print "Gbsa method %14s" % self._testGbviForce.getNonbondedMethod() print "Solute dielectric %14.7f" % (self._testGbviForce.getSoluteDielectric()/dimensionlessUnit) print "Solvent dielectric %14.7f" % (self._testGbviForce.getSolventDielectric()/dimensionlessUnit) print "Cutoff distance %14.7f" % (self._testNonbondedForce.getCutoffDistance()/lengthUnit) ii = 0 maxPrint = 10 while( ii < self._testGbviForce.getNumParticles() ): args = self._testGbviForce.getParticleParameters(ii) print "Gbvi: %5d [%14.7f %14.7f %14.7f]" % (ii, args[0]/chargeUnit, args[1]/lengthUnit, args[2]/energyUnit) ii += 1 if( ii == maxPrint and self._testGbviForce.getNumParticles() > maxPrint ):ii = self._testGbviForce.getNumParticles() - maxPrint ii = 0 while( ii < self._testGbviForce.getNumBonds() ): args = self._testGbviForce.getBondParameters(ii) print "Gbvi bonds: %5d [%6d %6d %14.7f]" % (ii, args[0], args[1], args[2]/lengthUnit) ii += 1 if( ii == maxPrint and self._testGbviForce.getNumBonds() > maxPrint ):ii = self._testGbviForce.getNumBonds() - maxPrint #============================================================================================= # Test class for force-energy gradient test for GbviForce w/ no cutoff #============================================================================================= class GbviForceEnergyGradientNoCutoffTest(GbviForceEnergyGradientBaseTest): """ GbviForceEnergyGradientNoCutoffTest create system for testing GbviForce w/ no cutoff inherits from GbviForceEnergyGradientBaseTest """ 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 """ GbviForceEnergyGradientBaseTest.__init__(self, parameterInput, defaultOverrides, verbose ) if( self.isActive() == 0 ): return self._testName = 'GbviForceEnergyGradientNoCutoffTest' self._testGbviForce.setNonbondedMethod( self._openMM.NonbondedForce.NoCutoff ) self._testNonbondedForce.setNonbondedMethod( self._openMM.NonbondedForce.NoCutoff ) if self._verbose: self.printGbvi()