#!/bin/env python #============================================================================================= # GLOBAL IMPORTS #============================================================================================= import os import os.path import sys import copy import re import math import traceback import doctest import time import simtk.openmm from ValidationUtilities import ValidationUtilities #============================================================================================= # Results class for force/energy platform comparison tests #============================================================================================= class FilePlatformComparisonResult(object): """ FilePlatformComparisonResult generates/computes results for ForceEnergy comparisons """ def __init__(self, test, platform, state1, state2, verbose=False): """ Generate results summary ARGUMENTS test (FilePlatformTest) - FilePlatformTest reference platforms (string list) - list of two string platform names state1 (OpenMM::State) - state object to retrieve forces/energy for first platform state2 (OpenMM::State) - state object to retrieve forces/energy for second platform verbose (Boolean) - if true, print logging info """ self._platform = platform self._states = [] self._states.append( state1 ) self._states.append( state2 ) self._openMM = simtk.openmm self._verbose = verbose self._validationUtilities = ValidationUtilities(verbose) self._minNormCutoff = 1.0e-02 self._systemName = test.getSystemName() self._testName = test.getTestName() self._tolerance = test.getTolerance() self._minForceNormCutoff = test._minForceNormCutoff self._minEnergyNormCutoff = test._minEnergyNormCutoff self._numberOfParticles = test.getNumberOfParticles() self._passed = -1 # get units from ValidationUtilities lengthUnit = self._validationUtilities.getLengthUnit( ) energyUnit = self._validationUtilities.getEnergyUnit( ) # check for differences in potential energies try: potential1 = state1.getPotentialEnergy() / energyUnit if( 'Energy' in state2 and self._validationUtilities.isNumber( str(state2['Energy']) ) ): potential2 = state2['Energy'] self._deltaPotentialEnergy = abs( potential2 - potential1 ) self._averagePotentials = 0.5*( abs( potential2 ) + abs( potential1 ) ) self._relativeDeltaPotentialEnergy = abs( potential2 - potential1 ) if( self._averagePotentials > 0.0 ): self._relativeDeltaPotentialEnergy = self._relativeDeltaPotentialEnergy/self._averagePotentials if self._verbose: print "PE %14.7e %14.7e %14.7e %14.7e " % ( potential1, potential2, self._deltaPotentialEnergy, self._relativeDeltaPotentialEnergy ) else: self._deltaPotentialEnergy = 0.0 self._relativeDeltaPotentialEnergy = 0.0 self._averagePotentials = 0.0 if self._verbose: print "PE %14.7e (not available from file)." % ( potential1 ) # check for differences in forces # compute delta and relative deltas # compute maximum deleta and maximum relative delta and track the indies where max occur forces1 = state1.getForces() forces2 = state2['Force'] forceUnit = energyUnit / lengthUnit self._maxDeltaInForce = -1.0 self._maxDeltaInForceIndex = -1 self._cntRelativeDeltaInForce = 0.0 self._avgRelativeDeltaInForce = 0.0 self._stdRelativeDeltaInForce = 0.0 self._maxRelativeDeltaInForce = -1.0 self._maxRelativeDeltaInForceIndex = -1 self._maxRelativeDeltaInForceNormSum = -1.0 self._oppositeSigns = 0 self._threshold = 0.1 self._count = 0 for ii in range( len( forces1 ) ): force1 = forces1[ii] / forceUnit force2 = forces2[ii] force1Norm = math.sqrt( force1[0]*force1[0] + force1[1]*force1[1] + force1[2]*force1[2] ) #print "%s\n%s fff" % (str(force2[0]), repr( force2 ) ) force2Norm = math.sqrt( force2[0]*force2[0] + force2[1]*force2[1] + force2[2]*force2[2] ) if force1Norm > 0.0 or force2Norm > 0.0: deltaNorm = abs( force2Norm - force1Norm ) relativeDeltaNorm = 2.0*deltaNorm/(force2Norm + force1Norm ) logOfRelativeDeltaNorm = math.log( relativeDeltaNorm ) self._avgRelativeDeltaInForce += logOfRelativeDeltaNorm self._stdRelativeDeltaInForce += logOfRelativeDeltaNorm*logOfRelativeDeltaNorm self._cntRelativeDeltaInForce += 1.0 self._count += 1.0 change = 0 sign = '' for jj in range( len(force1) ): if( ( (force1[jj]*force2[jj]) < 0.0) and ( (abs( force1[jj] ) > self._threshold ) or (abs( force2[jj] ) > self._threshold ) ) ): self._oppositeSigns = 1 if sign == '' : sign = "Signs opposite: " sign += "%d " % (jj) change = 1 if deltaNorm > self._maxDeltaInForce: self._maxDeltaInForce = deltaNorm self._maxDeltaInForceIndex = ii change = 1 if relativeDeltaNorm > self._maxRelativeDeltaInForce: self._maxRelativeDeltaInForce = relativeDeltaNorm self._maxRelativeDeltaInForceIndex = ii self._maxRelativeDeltaInForceNormSum = 0.5*(force2Norm + force1Norm ) change = 1 change = 1 if self._verbose and change: print "Force: %5d %14.7e %14.7e [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e] Nrms %14.7e %14.7e %s" % (ii, deltaNorm, relativeDeltaNorm, force1[0], force1[1], force1[2], force2[0], force2[1], force2[2], force1Norm, force2Norm, sign) if( self._cntRelativeDeltaInForce > 1 ): self._avgRelativeDeltaInForce /= self._cntRelativeDeltaInForce self._stdRelativeDeltaInForce = self._stdRelativeDeltaInForce - self._cntRelativeDeltaInForce*self._avgRelativeDeltaInForce*self._avgRelativeDeltaInForce print "FilePlatformComparison Result 0 %d %15.7e" % (self._cntRelativeDeltaInForce, self._stdRelativeDeltaInForce) sys.stdout.flush() self._stdRelativeDeltaInForce = math.sqrt( self._stdRelativeDeltaInForce/(self._cntRelativeDeltaInForce-1) ) self._avgRelativeDeltaInForce = math.exp( self._avgRelativeDeltaInForce ) self._stdRelativeDeltaInForce = math.exp( -self._stdRelativeDeltaInForce ) if self._verbose: self._avgRelativeDeltaInForce += relativeDeltaNorm print "MaxDeltaForce: %5d %14.7e %20s %20s %10s" % (self._maxDeltaInForceIndex, self._maxDeltaInForce, self.getTestName(), self.getSystemName(), self._platform) print "MaxRelativeDeltaForce: %5d %14.7e %20s %20s %10s" % (self._maxRelativeDeltaInForceIndex, self._maxRelativeDeltaInForce, self.getTestName(), self.getSystemName(), self._platform ) print "LogAvgRelativeDeltaForce: %14.7e %20s %20s %10s" % (self._avgRelativeDeltaInForce, self.getTestName(), self.getSystemName(), self._platform) print "LogStdRelativeDeltaForce: %14.7e %20s %20s %10s" % (self._stdRelativeDeltaInForce, self.getTestName(), self.getSystemName(), self._platform) sys.stdout.flush() except: print "FilePlatformComparisonResult problem: %s %s" % (self._systemName, self._testName ) print repr(traceback.print_exception( sys.exc_info()[0], sys.exc_info()[1], sys.exc_info()[2])) sys.stdout.flush() #raise #============================================================================================= # Get number of particles #============================================================================================= def getNumberOfParticles(self): return self._numberOfParticles #============================================================================================= # Return tolerance #============================================================================================= def getTolerance(self): return self._tolerance #============================================================================================= # Get platform #============================================================================================= def getPlatform(self): return self._platform #============================================================================================= # Get delta in potential energies #============================================================================================= def getDeltaPotentialEnergy(self): return self._deltaPotentialEnergy #============================================================================================= # Get relative delta potential energy #============================================================================================= def getRelativeDeltaPotentialEnergy(self): return self._relativeDeltaPotentialEnergy #============================================================================================= # Get average relative delta in force #============================================================================================= def getAvgRelativeDeltaInForce(self): return self._avgRelativeDeltaInForce #============================================================================================= # Get std dev relative delta in force #============================================================================================= def getStdRelativeDeltaInForce(self): return self._stdRelativeDeltaInForce #============================================================================================= # Get count for relative delta in force #============================================================================================= def getCntRelativeDeltaInForce(self): return self._cntRelativeDeltaInForce #============================================================================================= # Get MaxDelta #============================================================================================= def getMaxDeltaInForce(self): return self._maxDeltaInForce #============================================================================================= # Get MaxDeltaIndex #============================================================================================= def getMaxDeltaInForceIndex(self): return self._maxDeltaInForceIndex #============================================================================================= # Get MaxRelativeDelta #============================================================================================= def getMaxRelativeDeltaInForce(self): return self._maxRelativeDeltaInForce #============================================================================================= # Get MaxRelativeDeltaIndex #============================================================================================= def getMaxRelativeDeltaInForceIndex(self): return self._maxDeltaInForceIndex #============================================================================================= # Get MaxRelativeDeltaNormSum #============================================================================================= def getMaxRelativeDeltaInForceNormSum(self): return self._maxRelativeDeltaInForceNormSum #============================================================================================= # Get OppositeForceSigns #============================================================================================= def getOppositeForceSigns(self): return self._oppositeSigns; #============================================================================================= # Get system name #============================================================================================= def getSystemName(self): return self._systemName #============================================================================================= # Get system name #============================================================================================= def getTestName(self): return self._testName #============================================================================================= # Return 1 if test passsed, else 0 #============================================================================================= def testPassed(self): if( self._passed > -1 ): return self._passed sys.stdout.flush() forceTolerance = self._tolerance energyTolerance = forceTolerance self._passed = 1 if( self.getOppositeForceSigns() ): if( self._verbose ): print "Test failed: forces have opposite signs." self._passed = 0 if( self._maxRelativeDeltaInForce > forceTolerance and self._maxRelativeDeltaInForceNormSum > self._minForceNormCutoff ): if( self._verbose ): print "Test failed on force comparison: maxRelativeDelta %14.7e is greater than tolerance=%10.3e and NormSum=%14.7e is greater than threshold=%10.3e" % (self._maxRelativeDeltaInForce, forceTolerance, self._maxRelativeDeltaInForceNormSum, self._minForceNormCutoff) self._passed = 0 if( self._relativeDeltaPotentialEnergy > energyTolerance and self._averagePotentials > self._minEnergyNormCutoff ): if( self._verbose ): print "Test failed on energy comparison: maxRelativeDelta %14.7e is greater than tolerance=%10.3e and NormSum=%14.7e is greater than threshold=%10.3e" % (self._relativeDeltaPotentialEnergy, energyTolerance, self._averagePotentials, self._minNormCutoff ) self._passed = 0 if( self._verbose ): sys.stdout.flush() return self._passed #============================================================================================= # Base class for ForceEnergy comparison tests #============================================================================================= class FilePlatformTest(object): """Base class for ForceEnergy comparison tests Centralizes common functionality across ForceEnergy comparison 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 = 1.0e-02 # min cutoffs for norms if( 'MinForceNormCutoff' in defaultOverrides ): self._minForceNormCutoff = defaultOverrides['MinForceNormCutoff'] else: self._minForceNormCutoff = 1.0e-01 if( 'MinEnergyNormCutoff' in defaultOverrides ): self._minEnergyNormCutoff = defaultOverrides['MinEnergyNormCutoff'] else: self._minEnergyNormCutoff = 1.0e-01 self._isActive = 1 self._system = None self._systems = None self._result = None self._openMM = simtk.openmm self._verbose = verbose self._validationUtilities = ValidationUtilities(verbose) self._deviceId = -1 #============================================================================================= # Return 1 if test is active #============================================================================================= def isActive(self): return self._isActive #============================================================================================= # Return tolerance #============================================================================================= def getTolerance(self): return self._tolerance #============================================================================================= # Return minimum force norm cutoff #============================================================================================= def getMinForceNormCutoff(self): return self._minForceNormCutoff #============================================================================================= # Return minimum energy norm cutoff #============================================================================================= def getMinEnergyNormCutoff(self): return self._minEnergyNormCutoff #============================================================================================= # 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 #============================================================================================= # 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 #============================================================================================= # Include file forces in list #============================================================================================= def includeFileForces(self, forceList, parameterInput ): summedForces = None for force in forceList: forceVector = parameterInput.getVec3Array(force) if( forceVector is not None ): if( summedForces is None ): summedForces = forceVector else: self._validationUtilities.addForce( summedForces, forceVector ) return summedForces #============================================================================================= # Get number of particles #============================================================================================= def getNumberOfParticles(self): return self._numberOfParticles #============================================================================================= # Include file energies in list #============================================================================================= def includeFileEnergies(self, energyList, parameterInput ): summedEnergies = 0.0 for energy in energyList: energyValue = parameterInput.getScalar(energy) if( energyValue is not None ): #print "E %s %15.7e" % ( energy, energyValue) summedEnergies += energyValue return summedEnergies #============================================================================================= # Run test #============================================================================================= def runTest(self, platformName): if self._system is not None: self._numberOfParticles = self._system.getNumParticles() if self._verbose: print "runTest: %s %s" % (self.getFullTestName(), platformName ) if self._system is not None: print " Particles=%d NumForces=%d " % (self._system.getNumParticles(), self._system.getNumForces() ) if self._system is None and self._systems is None: print "runTest: %s no system available" % (self.getFullTestName() ) return timestep = 0.01 integrator = self._openMM.VerletIntegrator(timestep) system = self._system try: platform = self._openMM.Platform.getPlatformByName( platformName ) if( self._deviceId > -1 ): self._validationUtilities.setDeviceId( platform, self._deviceId ) context1 = self._openMM.Context(system, integrator, platform ) except: print "runTest: %s context not created for platform %s" % (self.getFullTestName(), platformName ) sys.stdout.flush() raise context1.setPositions( self._positions ) if self._verbose: print "Positions %d" % len( self._positions ) maxPrint = 10 ii = 0 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( range[0][jj] > coordinates[jj] ): range[0][jj] = coordinates[jj] if( range[1][jj] < coordinates[jj] ): 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 jj = 0 while( jj < 3): print "Enclosing Box %2d [%14.7f %14.7f] %14.7f" % ( jj, range[0][jj], range[1][jj], (range[1][jj] - range[0][jj]) ) jj += 1 self._validationUtilities.printPlatformProperties( context1 ) state1 = context1.getState(getEnergy=True, getForces=True) self._result = FilePlatformComparisonResult( self, platformName, state1, self._compareState, self._verbose ) if( self._verbose ): if( self._result.testPassed() ): print "%s passed." % self.getFullTestName() else: print "%s did not pass." % self.getFullTestName() sys.stdout.flush() return #============================================================================================= # Test class for force/energy platform comparison tests for HarmonicBondForce test #============================================================================================= class CMAPTorsionForceTest(FilePlatformTest): """ CMAPTorsionTest create system for testing CMAPTorsionFirce inherits from FilePlatformTest """ 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 """ FilePlatformTest.__init__(self, defaultOverrides, verbose ) self._systemName = parameterInput.getName().replace( '.txt', '' ) self._testName = 'CMAPTorsionForceTest' if( parameterInput.getNumberOfCMAPTorsions() == 0 ): if verbose: print "No CMAPTorsionForceTest entries for %s" % parameterInput.getName() return CMAPTorsionForce = parameterInput.getCMAPTorsionForce() # build system self._system = self._openMM.System() self._validationUtilities.loadSystemMasses( parameterInput.getMasses(), self._system ) testCMAPTorsionForce = self._validationUtilities.copyCMAPTorsionForce( CMAPTorsionForce ) self._system.addForce(testCMAPTorsionForce) if self._verbose: print "\n\n" + self.getFullTestName() print "Added %d CMAPTorsionForce entries " % testCMAPTorsionForce.getNumTorsions() ii = 0 lengthUnit = self._validationUtilities.getLengthUnit( ) forceConstantUnit = self._validationUtilities.getEnergyUnit( ) / lengthUnit**2 maxPrint = 10 while( ii < testCMAPTorsionForce.getNumTorsions() ): args = testCMAPTorsionForce.getTorsionParameters(ii) print "CMAP Torsion: %5d gridIndex=%2d [%5d %5d %5d %5d] [%5d %5d %5d %5d] " % (ii, args[0], args[1], args[2], args[3], args[4], args[5], args[6], args[7], args[8] ) ii += 1 if( ii == maxPrint and testCMAPTorsionForce.getNumTorsions() > 2*maxPrint ):ii = testCMAPTorsionForce.getNumTorsions() - maxPrint sys.stdout.flush() self._positions = parameterInput.getVec3Array('Positions') self._compareState = {} self._compareState['Force'] = parameterInput.getVec3Array('TorsionTorsionForce') self._compareState['Energy'] = parameterInput._scalars['TorsionTorsionEnergy'] #============================================================================================= # Test class for force/energy comparison for HarmonicBondForce w/ results from Gromacs stored in a file #============================================================================================= class HarmonicBondForceGromacsTest(FilePlatformTest): """ create system for testing HarmonicBondForce and comparing w/ Gromacs results inherits from FilePlatformTest """ 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 """ FilePlatformTest.__init__(self, defaultOverrides, verbose ) self._systemName = parameterInput.getName().replace( '.txt', '' ) self._testName = 'HarmonicBondForceGromacsTest' 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 harmonicBondForce entries " % 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() > 2*maxPrint ):ii = testHarmonicBondForce.getNumBonds() - maxPrint sys.stdout.flush() self._positions = parameterInput.getVec3Array('Positions') self._compareState = {} self._compareState['Force'] = parameterInput.getVec3Array('GromacsHarmonicBondForce') self._compareState['Energy'] = parameterInput.getScalar('GromacsHarmonicBondEnergy') #============================================================================================= # Test class for force/energy comparison for HarmonicAngleForce w/ results from Gromacs stored in a file #============================================================================================= class HarmonicAngleForceGromacsTest(FilePlatformTest): """ create system for testing HarmonicAngleForce and comparing w/ Gromacs results inherits from FilePlatformTest """ 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 """ FilePlatformTest.__init__(self, defaultOverrides, verbose ) self._systemName = parameterInput.getName().replace( '.txt', '' ) self._testName = 'HarmonicAngleForceGromacsTest' 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: maxPrint = 10 self._validationUtilities.printHarmonicAngleForce( self._testName, testHarmonicAngleForce, maxPrint ) self._positions = parameterInput.getVec3Array('Positions') self._compareState = {} self._compareState['Force'] = parameterInput.getVec3Array('GromacsHarmonicAngleForce') self._compareState['Energy'] = parameterInput.getScalar('GromacsHarmonicAngleEnergy') #============================================================================================= # Test class for force/energy comparison for PeriodicTorsionForce w/ results from Gromacs stored in a file #============================================================================================= class PeriodicTorsionForceGromacsTest(FilePlatformTest): """ create system for testing PeriodicTorsionForce and comparing w/ Gromacs results inherits from FilePlatformTest """ 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 """ FilePlatformTest.__init__(self, defaultOverrides, verbose ) self._systemName = parameterInput.getName().replace( '.txt', '' ) self._testName = 'PeriodicTorsionForceGromacsTest' 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: maxPrint = 10 self._validationUtilities.printPeriodicTorsionForce( self.getFullTestName(), testPeriodicTorsionForce, maxPrint ) self._positions = parameterInput.getPositions() self._compareState = {} self._compareState['Force'] = parameterInput.GromacsHarmonicBondForce('GromacsPeriodicTorsionForce') self._compareState['Energy'] = parameterInput.getScalar('GromacsPeriodicTorsionEnergy') #============================================================================================= # Test class for force/energy comparison for PeriodicTorsionForce w/ results from Gromacs stored in a file # The comparison is based on the sum of the GromacsProperDihedralForce and GromacsProperIDihedralForce #============================================================================================= class PeriodicTorsionForceGromacsDihedralTest(FilePlatformTest): """ create system for testing PeriodicTorsionForce and comparing w/ Gromacs results inherits from FilePlatformTest """ 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 """ FilePlatformTest.__init__(self, defaultOverrides, verbose ) self._systemName = parameterInput.getName().replace( '.txt', '' ) self._testName = 'PeriodicTorsionForceGromacsDihedralTest' 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: maxPrint = 10 self._validationUtilities.printPeriodicTorsionForce( self.getFullTestName(), testPeriodicTorsionForce, maxPrint ) self._positions = parameterInput.getPositions() self._compareState = {} forceList = [ 'GromacsProperDihedralForce', 'GromacsProperIDihedralForce' ] self._compareState['Force'] = self.includeFileForces( forceList, parameterInput ) energyList = [ 'GromacsProperDihedralEnergy', 'GromacsProperIDihedralEnergy' ] self._compareState['Energy'] = self.includeFileEnergies( energyList, parameterInput ) #============================================================================================= # Test class for force/energy comparison for RBTorsionForce w/ results from Gromacs stored in a file #============================================================================================= class RBTorsionForceGromacsTest(FilePlatformTest): """ create system for testing RBTorsionForce and comparing w/ Gromacs results inherits from FilePlatformTest """ 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 """ FilePlatformTest.__init__(self, defaultOverrides, verbose ) self._systemName = parameterInput.getName().replace( '.txt', '' ) self._testName = 'RBTorsionForceGromacsTest' 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() > 2*maxPrint ):ii = testRBTorsionForce.getNumTorsions() - maxPrint sys.stdout.flush() self._positions = parameterInput.getPositions() self._compareState = {} self._compareState['Force'] = parameterInput.getVec3Array('GromacsRBTorsionForce') self._compareState['Energy'] = parameterInput.getScalar('GromacsRBTorsionEnergy') #============================================================================================= # Test class for force/energy comparison for NonbondedForce w/ results from Gromacs stored in a file #============================================================================================= class NonbondedForceGromacsTest(FilePlatformTest): """ create system for testing NonbondedForce and comparing w/ Gromacs results inherits from FilePlatformTest """ 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 """ FilePlatformTest.__init__(self, defaultOverrides, verbose ) self._systemName = parameterInput.getName().replace( '.txt', '' ) self._testName = 'NonbondedForceGromacsTest' if( parameterInput.getNumberOfNonbondeds() == 0 ): if verbose: print "No nonbonded ixns for %s" % parameterInput.getName() return nonbondedForce = parameterInput.getNonbondedForce() # build system self._system = self._openMM.System() self._validationUtilities.loadSystemMasses( parameterInput.getMasses(), self._system ) self._validationUtilities.copyDefaultBoxVectors( parameterInput.getDefaultBoxVectors(), self._system ) testNonbondedForce = self._validationUtilities.copyNonbondedForce( nonbondedForce, 1 ) self._system.addForce(testNonbondedForce) if self._verbose: self._validationUtilities.printNonbonded( self._testName, testNonbondedForce ) self._validationUtilities.printDefaultBoxVectors( self._testName, self._system ) self._positions = parameterInput.getPositions() self._compareState = {} forceList = [ 'GromacsNonbondedForce', 'GromacsNonbondedForceExceptions', 'GromacsCOUL_LRForce', 'Gromacs14Force' ] #forceList = [ 'GromacsNonbondedForce', 'GromacsNonbondedForceExceptions', 'Gromacs14Force' ] self._compareState['Force'] = self.includeFileForces( forceList, parameterInput ) energyList = [ 'GromacsNonbondedEnergy', 'GromacsNonbondedEnergyExceptions', 'GromacsCOUL_LREnergy', 'Gromacs14Energy' ] self._compareState['Energy'] = self.includeFileEnergies( energyList, parameterInput ) #============================================================================================= # Test class for force/energy comparison for NonbondedForce w/ results from Gromacs stored in a file #============================================================================================= class NonbondedLROnlyForceGromacsTest(FilePlatformTest): """ create system for testing NonbondedForce and comparing w/ Gromacs results inherits from FilePlatformTest """ 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 """ FilePlatformTest.__init__(self, defaultOverrides, verbose ) self._systemName = parameterInput.getName().replace( '.txt', '' ) self._testName = 'NonbondedLROnlyForceGromacsTest' if( parameterInput.getNumberOfNonbondeds() == 0 ): if verbose: print "No nonbonded ixns for %s" % parameterInput.getName() return nonbondedForce = parameterInput.getNonbondedForce() # build system self._system = self._openMM.System() self._validationUtilities.loadSystemMasses( parameterInput.getMasses(), self._system ) self._validationUtilities.copyDefaultBoxVectors( parameterInput.getDefaultBoxVectors(), self._system ) # zero out 14 ixns testNonbondedForce = self._validationUtilities.copyNonbondedForce( nonbondedForce, 2 ) self._system.addForce(testNonbondedForce) if self._verbose: self._validationUtilities.printNonbonded( self._testName, testNonbondedForce ) self._validationUtilities.printDefaultBoxVectors( self._testName, self._system ) self._positions = parameterInput.getPositions() self._compareState = {} forceList = [ 'GromacsNonbondedForce', 'GromacsCOUL_LRForce', ] #forceList = [ 'GromacsNonbondedForce', 'GromacsNonbondedForceExceptions', 'Gromacs14Force' ] self._compareState['Force'] = self.includeFileForces( forceList, parameterInput ) energyList = [ 'GromacsNonbondedEnergy', 'GromacsCOUL_LREnergy' ] self._compareState['Energy'] = self.includeFileEnergies( energyList, parameterInput ) #============================================================================================= # Test class for force/energy comparison for 14 ixns of NonbondedForce w/ results from Gromacs stored in a file #============================================================================================= class Nonbonded14OnlyForceGromacsTest(FilePlatformTest): """ create system for testing NonbondedForce and comparing w/ Gromacs results inherits from FilePlatformTest """ 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 """ FilePlatformTest.__init__(self, defaultOverrides, verbose ) self._systemName = parameterInput.getName().replace( '.txt', '' ) self._testName = 'Nonbonded14OnlyForceGromacsTest' if( parameterInput.getNumberOfNonbondeds() == 0 ): if verbose: print "No nonbonded ixns for %s" % parameterInput.getName() return nonbondedForce = parameterInput.getNonbondedForce() # build system self._system = self._openMM.System() self._validationUtilities.loadSystemMasses( parameterInput.getMasses(), self._system ) self._validationUtilities.copyDefaultBoxVectors( parameterInput.getDefaultBoxVectors(), self._system ) testNonbondedForce = self._validationUtilities.copyNonbondedForce( nonbondedForce, -2 ) self._system.addForce(testNonbondedForce) if self._verbose: self._validationUtilities.printNonbonded( self._testName, testNonbondedForce ) self._validationUtilities.printDefaultBoxVectors( self._testName, self._system ) self._positions = parameterInput.getPositions() self._compareState = {} forceList = [ 'GromacsNonbondedForceExceptions', 'Gromacs14Force' ] self._compareState['Force'] = self.includeFileForces( forceList, parameterInput ) energyList = [ 'GromacsNonbondedEnergyExceptions', 'Gromacs14Energy' ] self._compareState['Energy'] = self.includeFileEnergies( energyList, parameterInput ) #============================================================================================= # Test class for force/energy comparison for ObcForce w/ results from Gromacs stored in a file #============================================================================================= class ObcForceGromacsTest(FilePlatformTest): """ create system for testing ObcForce and comparing w/ Gromacs results inherits from FilePlatformTest """ 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 """ FilePlatformTest.__init__(self, defaultOverrides, verbose ) self._systemName = parameterInput.getName().replace( '.txt', '' ) self._testName = 'ObcForceGromacsTest' if( parameterInput.getNumberOfGbsaObcs() == 0 ): if verbose: print "No obc ixns for %s" % parameterInput.getName() return obcForce = parameterInput.getGbsaObcForce() # build system self._system = self._openMM.System() self._validationUtilities.loadSystemMasses( parameterInput.getMasses(), self._system ) self._validationUtilities.copyDefaultBoxVectors( parameterInput.getDefaultBoxVectors(), self._system ) testObcForce = self._validationUtilities.copyGbsaObcForce( obcForce ) self._system.addForce(testObcForce) if self._verbose: self._validationUtilities.printGbsaObc( self._testName, testObcForce ) self._positions = parameterInput.getPositions() self._compareState = {} self._compareState['Force'] = parameterInput.getVec3Array('GromacsGbsaObcForce') self._validationUtilities.subtractForce( self._compareState['Force'], parameterInput.getVec3Array('GromacsNonbondedForce') ) # factor 0.9875 was due bug in Gromacs 4.1 from which this file was generated # Obc forces obtained from ~/bptiComparisonGmx4.1Obc self._validationUtilities.multiplyForce( self._compareState['Force'], 0.9875 ) #============================================================================================= # Test class for force/energy comparison for ObcForce w/ results from Gromacs stored in a file #============================================================================================= class ObcAndNonbondedForceGromacsTest(FilePlatformTest): """ create system for testing Nonbonded + ObcForce and comparing w/ Gromacs results inherits from FilePlatformTest """ 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 """ FilePlatformTest.__init__(self, defaultOverrides, verbose ) self._systemName = parameterInput.getName().replace( '.txt', '' ) self._testName = 'ObcAndNonbondedForceGromacsTest' if( parameterInput.getNumberOfGbsaObcs() == 0 ): if verbose: print "No obc ixns for %s" % parameterInput.getName() return # build system: nonbonded + obcGBSA self._system = self._openMM.System() self._validationUtilities.loadSystemMasses( parameterInput.getMasses(), self._system ) obcForce = parameterInput.getGbsaObcForce() testObcForce = self._validationUtilities.copyGbsaObcForce( obcForce ) testObcForce.setNonbondedMethod( self._openMM.NonbondedForce.NoCutoff ) self._system.addForce(testObcForce) # includeExceptions == 2 -> eps and q exception values are zeroed nonbondedForce = parameterInput.getNonbondedForce() testNonbondedForce = self._validationUtilities.copyNonbondedForce( nonbondedForce, 2 ) testNonbondedForce.setNonbondedMethod( self._openMM.NonbondedForce.NoCutoff ) self._system.addForce(testNonbondedForce) if self._verbose: self._validationUtilities.printNonbonded( self._testName, testNonbondedForce ) self._validationUtilities.printGbsaObc( self._testName, testObcForce ) self._positions = parameterInput.getPositions() self._compareState = {} #forceList = [ 'GromacsNonbondedForce', 'GromacsNonbondedForceExceptions', 'GromacsCOUL_LRForce', 'Gromacs14Force', 'GromacsGBPolForce' ] forceList = [ 'GromacsNonbondedForce', 'GromacsNonbondedForceExceptions', 'GromacsGBPolForce', 'GromacsCOUL_LRForce' ] self._compareState['Force'] = self.includeFileForces( forceList, parameterInput ) energyList = [ 'GromacsNonbondedForce', 'GromacsNonbondedEnergyExceptions', 'GromacsCOUL_LREnergy', 'GromacsGBPolEnergy' ] self._compareState['Energy'] = self.includeFileEnergies( energyList, parameterInput )