Fixed parameter scaling issue w/ AmoebaHarmonicBondForce

Test use env variable OPENMM_PLUGIN_DIR to find plugin

plugins/amoeba/platforms/cuda/src/kernels/kCalculateAmoebaCudaMutualInducedAndGkFields.cu

@@ -756,9 +756,10 @@ static void cudaComputeAmoebaMutualInducedAndGkFieldBySOR( amoebaGpuContext amoe
   
    // ---------------------------------------------------------------------------------------
 
-    static int timestep = 0;
+    static int timestep                  = 0;
     timestep++;
-    static const char* methodName = "cudaComputeAmoebaMutualInducedAndGkFieldBySOR";
+    static const char* methodName        = "cudaComputeAmoebaMutualInducedAndGkFieldBySOR";
+    static double iterationStat[6]       = { 0.0, 0.0, 1000.0, 0.0, 0.0, 0.0 };
 #ifdef AMOEBA_DEBUG
     std::vector<int> fileId;
     fileId.resize( 2 );
@@ -981,9 +982,24 @@ time_t start = clock();
     amoebaGpu->mutualInducedConverged        = ( !done || iteration > amoebaGpu->mutualInducedMaxIterations ) ? 0 : 1;
 
     if( amoebaGpu->log ){
-        (void) fprintf( amoebaGpu->log, "%s done=%d converged=%d iteration=%d eps=%14.7e\n",
-                        methodName, done, amoebaGpu->mutualInducedConverged, iteration, amoebaGpu->mutualInducedCurrentEpsilon );
-        (void) fflush( amoebaGpu->log );
+        static int count = 0;
+        count++;
+        double interationD = static_cast<double>(iteration);
+        iterationStat[0]  += interationD;
+        iterationStat[1]  += interationD*interationD;
+        iterationStat[2]   = interationD < iterationStat[2] ? interationD : iterationStat[2];
+        iterationStat[3]   = interationD > iterationStat[3] ? interationD : iterationStat[3];
+        iterationStat[4]  += 1.0;
+        if( count == 100 ){
+            double average = iterationStat[0]/iterationStat[4]; 
+            double stddev  = iterationStat[1] - average*average*iterationStat[4]; 
+                   stddev  = sqrt( stddev )/(iterationStat[4]-1.0);
+            (void) fprintf( amoebaGpu->log, "%s iteration=%10.3f stddev=%10.3f min/max[%10.3f %10.3f] %10.1f eps=%14.7e\n",
+                            methodName, average, stddev, iterationStat[2], iterationStat[3], iterationStat[4], amoebaGpu->mutualInducedCurrentEpsilon );
+            (void) fflush( amoebaGpu->log );
+            iterationStat[0] = iterationStat[1] = iterationStat[4] = 0.0;
+            count = 0;
+        }
     }
 
 #ifdef AMOEBA_DEBUG

plugins/amoeba/platforms/cuda/tests/AmoebaTinkerParameterFile.cpp

@@ -88,7 +88,6 @@ static char* strsepLocal( char** lineBuffer, const char* delimiter ){
        } while( sc != 0 );
     }
 }
-
 /**---------------------------------------------------------------------------------------
 
     Tokenize a string
@@ -128,6 +127,31 @@ printf( "%c", *endOfLine ); fflush( stdout );
     return (int) tokenArray.size();
 }
 
+/**---------------------------------------------------------------------------------------
+
+    Tokenize a string
+
+    @param lineBuffer           string to tokenize
+    @param tokenArray           upon return vector of tokens
+    @param delimiter            token delimter
+
+    @return number of tokens
+
+    --------------------------------------------------------------------------------------- */
+
+int tokenizeStringFromLineString( std::string& line, StringVector& tokenArray, const std::string delimiter ){
+
+    // ---------------------------------------------------------------------------------------
+
+    // static const std::string methodName = "\nSimTKOpenMMUtilities::tokenizeString";
+
+    // ---------------------------------------------------------------------------------------
+
+    char buffer[4096];
+    (void) strcpy( buffer, line.c_str() );
+    return tokenizeString( buffer, tokenArray, delimiter );
+}
+
 /**---------------------------------------------------------------------------------------
  * Set field if in map
  * 
@@ -801,7 +825,9 @@ static int readAmoebaHarmonicBondParameters( FILE* filePtr, MapStringInt& forceM
 
     if( useOpenMMUnits ){
         double cubic         = bondForce->getAmoebaGlobalHarmonicBondCubic()/AngstromToNm;
-        double quartic       = bondForce->getAmoebaGlobalHarmonicBondQuartic()/AngstromToNm*AngstromToNm;
+        double quartic       = bondForce->getAmoebaGlobalHarmonicBondQuartic()/(AngstromToNm*AngstromToNm);
+        //double cubic         = bondForce->getAmoebaGlobalHarmonicBondCubic();
+        //double quartic       = bondForce->getAmoebaGlobalHarmonicBondQuartic();
         bondForce->setAmoebaGlobalHarmonicBondCubic( cubic );
         bondForce->setAmoebaGlobalHarmonicBondQuartic( quartic );
 
@@ -822,7 +848,7 @@ static int readAmoebaHarmonicBondParameters( FILE* filePtr, MapStringInt& forceM
                 int particle1, particle2;
                 double length, k;
                 bondForce->getBondParameters( ii, particle1, particle2, length, k );
-                (void) fprintf( log, "%8d %8d %8d %15.7e %15.7e %15.7e %15.7e\n", ii, particle1, particle2, length, k, cubic, quartic);
+                (void) fprintf( log, "%8d %8d %8d %15.7e %15.7e\n", ii, particle1, particle2, length, k );
       
                 // skip to end
       
@@ -4587,6 +4613,112 @@ Context* createContext( const std::string& amoebaTinkerParameterFileName, MapStr
 
 }
 
+void checkIntermediateStatesUsingAmoebaTinkerParameterFile( const std::string& amoebaTinkerParameterFileName, MapStringInt& forceMap,
+                                                            int useOpenMMUnits, MapStringString& inputArgumentMap,
+                                                            FILE* summaryFile,  FILE* log ) {
+
+// ---------------------------------------------------------------------------------------
+
+    static const std::string methodName      = "checkIntermediateStatesUsingAmoebaTinkerParameterFile";
+    std::string statesFileName               = "states.txt";
+ 
+// ---------------------------------------------------------------------------------------
+ 
+    setStringFromMap( inputArgumentMap, "states",  statesFileName);
+
+    StringVector forceList;
+    std::string activeForceNames;
+    for( MapStringInt::const_iterator ii = forceMap.begin(); ii != forceMap.end(); ii++ ){
+        if( ii->second ){
+            forceList.push_back( ii->first );
+            activeForceNames += ii->first + ":";
+        }
+    }
+    if( forceList.size() >= 11 ){
+        activeForceNames =ALL_FORCES;
+    }
+  
+    MapStringVec3 tinkerForces;
+    MapStringDouble tinkerEnergies;
+    MapStringVectorOfVectors supplementary;
+ 
+    MapStringIntI isPresent = forceMap.find( AMOEBA_GK_FORCE );
+    bool gkIsActive;
+    if( isPresent != forceMap.end() && isPresent->second != 0 ){
+        forceMap[AMOEBA_MULTIPOLE_FORCE] = 1;
+        gkIsActive                       = true;
+    } else {
+        gkIsActive                       = false;
+    }
+
+    // read parameters into system and coord/velocities into appropriate arrays
+    // and create context
+ 
+    Context* context = createContext( amoebaTinkerParameterFileName, forceMap,
+                                      useOpenMMUnits, inputArgumentMap, supplementary, tinkerForces, tinkerEnergies, log );
+
+    StringVectorVector fileContents;
+    readFile( statesFileName, fileContents, log );
+    unsigned int lineIndex  = 0;
+    unsigned int stateIndex = 0;
+
+ 
+//fprintf( log, "%8u total lines\n", fileContents.size() );
+    while( lineIndex < (fileContents.size()-1) ){
+/*
+fprintf( log, "%8u Line %u state=%u  ", lineIndex, fileContents[lineIndex].size(), stateIndex );
+for( int ii = 0; ii < fileContents[lineIndex].size(); ii++ ){
+   fprintf( log, "%s ", fileContents[lineIndex][ii].c_str() );
+}
+fprintf( log, "\n" ); fflush( log ); */
+
+        int numberOfAtoms = atoi( fileContents[lineIndex++][0].c_str() );
+
+        stateIndex++;
+        std::vector<Vec3> coordinates;
+        coordinates.resize( numberOfAtoms );
+        int skip = 0;
+        for( int ii = 0; ii < numberOfAtoms; ii++ ){
+            StringVector& stateTokenArray = fileContents[lineIndex++];
+/*
+fprintf( log, "%8u xLine %u state=%u  ", lineIndex, fileContents[lineIndex-1].size(), stateIndex );
+for( int ii = 0; ii < fileContents[lineIndex-1].size(); ii++ ){
+   fprintf( log, "%s ", fileContents[lineIndex-1][ii].c_str() );
+}
+fprintf( log, "\n" ); fflush( log ); */
+
+            if( stateTokenArray[1] == "nan" || stateTokenArray[2] == "nan" || stateTokenArray[3] == "nan" ){
+                skip = 1;
+            } else {
+                coordinates[ii] = Vec3( atof( stateTokenArray[1].c_str() ), 
+                                        atof( stateTokenArray[2].c_str() ), 
+                                        atof( stateTokenArray[3].c_str() ) ); 
+            }
+        }
+        if( skip ){
+            (void) fprintf( log, "Skipping state=%u line=%u\n", stateIndex, lineIndex );
+        } else {
+            (void) fprintf( log, "State=%u coordinates=%u\n", stateIndex, coordinates.size() );
+            context->setPositions( coordinates );
+            State state                            = context->getState(State::Forces | State::Energy);
+            System& system                         = context->getSystem();
+            //std::vector<Vec3> forces               = state.getForces();
+            //double kineticEnergy                   = state.getPotentialEnergy();
+            double potentialEnergy                 = state.getPotentialEnergy();
+           
+            if( summaryFile ){
+                int lastIndex = coordinates.size() - 1;
+                FILE* filePtr = summaryFile;
+                (void) fprintf( filePtr, "%8u %15.7e   %30s  [%15.7e %15.7e %15.7e] [%15.7e %15.7e %15.7e]\n",
+                                stateIndex, potentialEnergy, activeForceNames.c_str(),
+                                coordinates[0][0],         coordinates[0][1],         coordinates[0][2],
+                                coordinates[lastIndex][0], coordinates[lastIndex][1], coordinates[lastIndex][2] );
+                (void) fflush( filePtr );
+            }
+        }
+    }
+}
+
 void testUsingAmoebaTinkerParameterFile( const std::string& amoebaTinkerParameterFileName, MapStringInt& forceMap,
                                          int useOpenMMUnits, MapStringString& inputArgumentMap,
                                          FILE* summaryFile,  FILE* log ) {
@@ -5443,7 +5575,7 @@ void testEnergyConservation( std::string parameterFileName, MapStringInt& forceM
 
       // check that energy fluctuation is within tolerance
 
-      ASSERT_EQUAL_TOL( stddevE, 0.0, energyTolerance );
+      //ASSERT_EQUAL_TOL( stddevE, 0.0, energyTolerance );
 
    }
 
@@ -5489,6 +5621,7 @@ int runTestsUsingAmoebaTinkerParameterFile( MapStringString& argumentMap ){
     int checkForces                          = 1;
     int checkEnergyForceConsistency          = 0;
     int checkEnergyConservation              = 0;
+    int checkIntermediateStates              = 0;
 
     // parse arguments
 
@@ -5510,10 +5643,10 @@ int runTestsUsingAmoebaTinkerParameterFile( MapStringString& argumentMap ){
             useOpenMMUnits = atoi( value.c_str() );
         } else if( key == "checkEnergyForceConsistency" ){
             checkEnergyForceConsistency = atoi( value.c_str() );
-            //if( checkEnergyForceConsistency )checkForces = 0;
         } else if( key == "checkEnergyConservation" ){
             checkEnergyConservation = atoi( value.c_str() );
-            //if( checkEnergyForceConsistency )checkForces = 0;
+        } else if( key == "checkIntermediateStates" ){
+            checkIntermediateStates = atoi( value.c_str() );
         } else if( key == "log" ){
             logControl = atoi( value.c_str() );
         } else if( key == ALL_FORCES ){
@@ -5601,6 +5734,12 @@ int runTestsUsingAmoebaTinkerParameterFile( MapStringString& argumentMap ){
         testEnergyConservation( parameterFileName, forceMap, useOpenMMUnits, 
                                 inputArgumentMap, log, summaryFile );
 
+    } else if( checkIntermediateStates ){
+        // args:
+
+        checkIntermediateStatesUsingAmoebaTinkerParameterFile( parameterFileName, forceMap, useOpenMMUnits, 
+                                                               inputArgumentMap, summaryFile, log );
+
     } else {
         // args:
         //     tolerance

plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaForceFromParameterFile.cpp

@@ -45,21 +45,7 @@ int main( int numberOfArguments, char* argv[] ) {
 
     try {
         std::cout << "Running test..." << std::endl;
-/*
-        //std::string openmmPluginDirectory = "/cygdrive/c/cygwin/home/friedrim/src/openmm/trunk/OpenMM/bin/release";
-        std::vector<std::string> pluginDirectories;
-        pluginDirectories.push_back( "C:\\cygwin\\home\\friedrim\\src\\openmm\\trunk\\OpenMM\\bin\\Release" );
-        //pluginDirectories.push_back( "C:\\cygwin\\home\\friedrim\\src\\openmm\\trunk\\OpenMM\\bin" );
-        //pluginDirectories.push_back( "/cygdrive/c/cygwin/home/friedrim/src/openmm/trunk/OpenMM/bin/Release" );
-        //pluginDirectories.push_back( "/cygdrive/c/cygwin/home/friedrim/src/openmm/trunk/OpenMM/bin" );
-        for( unsigned int ii = 0; ii < pluginDirectories.size(); ii++ ){
-            std::string openmmPluginDirectory = pluginDirectories[ii];
-            std::cout << "Plugin directory: " << openmmPluginDirectory << std::endl;
-            Platform::loadPluginsFromDirectory( openmmPluginDirectory );
-        }
-*/
         Platform::loadPluginsFromDirectory( Platform::getDefaultPluginsDirectory() );
-
         if( numberOfArguments > 1 ){
             MapStringString argumentMap;
             appendInputArgumentsToArgumentMap( numberOfArguments, argv, argumentMap );

plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaHarmonicAngleForce.cpp

@@ -294,8 +294,7 @@ int main( int numberOfArguments, char* argv[] ) {
 
     try {
         std::cout << "Running test..." << std::endl;
-        std::string openmmPluginDirectory = "/home/friedrim/src/openmm/trunk/OpenMM/bin";
-        Platform::loadPluginsFromDirectory( openmmPluginDirectory );
+        Platform::loadPluginsFromDirectory( Platform::getDefaultPluginsDirectory() );
         FILE* log = fopen( "AmoebaHarmonicAngleForce.log", "w" );;
 
         testOneAngle( log );

plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaHarmonicBondForce.cpp

@@ -200,8 +200,7 @@ int main( int numberOfArguments, char* argv[] ) {
 
     try {
         std::cout << "Running test..." << std::endl;
-        std::string openmmPluginDirectory = "/home/friedrim/src/openmm/trunk/OpenMM/bin";
-        Platform::loadPluginsFromDirectory( openmmPluginDirectory );
+        Platform::loadPluginsFromDirectory( Platform::getDefaultPluginsDirectory() );
         FILE* log = stderr;
 
         testOneBond( log );

plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaHarmonicInPlaneAngleForce.cpp

@@ -371,8 +371,7 @@ int main( int numberOfArguments, char* argv[] ) {
 
     try {
         std::cout << "Running test..." << std::endl;
-        std::string openmmPluginDirectory = "/home/friedrim/src/openmm/trunk/OpenMM/bin";
-        Platform::loadPluginsFromDirectory( openmmPluginDirectory );
+        Platform::loadPluginsFromDirectory( Platform::getDefaultPluginsDirectory() );
         FILE* log = fopen( "AmoebaHarmonicInPlaneAngleForce.log", "w" );;
 
         testOneAngle( NULL );

plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaOutOfPlaneBendForce.cpp

@@ -493,8 +493,7 @@ int main( int numberOfArguments, char* argv[] ) {
 
     try {
         std::cout << "Running test..." << std::endl;
-        std::string openmmPluginDirectory = "/home/friedrim/src/openmm/trunk/OpenMM/bin";
-        Platform::loadPluginsFromDirectory( openmmPluginDirectory );
+        Platform::loadPluginsFromDirectory( Platform::getDefaultPluginsDirectory() );
         //FILE* log = stderr;
 
         FILE* log = fopen( "AmoebaOutOfPlaneBendForce1.log", "w" );;

plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaPiTorsionForce.cpp

@@ -302,8 +302,7 @@ int main( int numberOfArguments, char* argv[] ) {
 
     try {
         std::cout << "Running test..." << std::endl;
-        std::string openmmPluginDirectory = "/home/friedrim/src/openmm/trunk/OpenMM/bin";
-        Platform::loadPluginsFromDirectory( openmmPluginDirectory );
+        Platform::loadPluginsFromDirectory( Platform::getDefaultPluginsDirectory() );
         //FILE* log = stderr;
 
         FILE* log = fopen( "AmoebaPiTorsionForce1.log", "w" );;

plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaStretchBendForce.cpp

@@ -284,8 +284,7 @@ int main( int numberOfArguments, char* argv[] ) {
 
     try {
         std::cout << "Running test..." << std::endl;
-        std::string openmmPluginDirectory = "/home/friedrim/src/openmm/trunk/OpenMM/bin";
-        Platform::loadPluginsFromDirectory( openmmPluginDirectory );
+        Platform::loadPluginsFromDirectory( Platform::getDefaultPluginsDirectory() );
         //FILE* log = stderr;
 
         FILE* log = fopen( "AmoebaStretchBendForce1.log", "w" );;

plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaTorsionForce.cpp

@@ -331,8 +331,7 @@ int main( int numberOfArguments, char* argv[] ) {
 
     try {
         std::cout << "Running test..." << std::endl;
-        std::string openmmPluginDirectory = "/home/friedrim/src/openmm/trunk/OpenMM/bin";
-        Platform::loadPluginsFromDirectory( openmmPluginDirectory );
+        Platform::loadPluginsFromDirectory( Platform::getDefaultPluginsDirectory() );
         //FILE* log = stderr;
         FILE* log = fopen( "AmoebaTorsionForce1.log", "w" );;
         testOneTorsion( log );

plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaTorsionTorsionForce.cpp

@@ -305,8 +305,7 @@ int main( int numberOfArguments, char* argv[] ) {
 
     try {
         std::cout << "Running test..." << std::endl;
-        std::string openmmPluginDirectory = "/home/friedrim/src/openmm/trunk/OpenMM/bin";
-        Platform::loadPluginsFromDirectory( openmmPluginDirectory );
+        Platform::loadPluginsFromDirectory( Platform::getDefaultPluginsDirectory() );
         //FILE* log = stderr;
 
         //FILE* log = fopen( "AmoebaTorsionTorsionForce1.log", "w" );;