Fixed memory bug in torque/force mapping;

An exception is now thrown if the multipole axis type is not Bisector or Z-then-X
Redid velocity initialization in test code

plugins/amoeba/openmmapi/include/AmoebaMultipoleForce.h

@@ -33,18 +33,21 @@
  * -------------------------------------------------------------------------- */
 
 #include "openmm/Force.h"
-#include <vector>
+#include "openmm/OpenMMException.h"
 #include "openmm/internal/windowsExport.h"
 
+#include <sstream>
+#include <vector>
+
 namespace OpenMM {
 
 /**
- * This class implements the Amoeba pi-torsion interaction
- * To use it, create a MultipoleForce object then call addMultipole() once for each torsion.  After
- * a torsion has been added, you can modify its force field parameters by calling setMultipoleParameters().
+ * This class implements the Amoeba multipole interaction
+ * To use it, create a MultipoleForce object then call addMultipole() once for each atom.  After
+ * a entry has been added, you can modify its force field parameters by calling setMultipoleParameters().
  */
 
- class OPENMM_EXPORT AmoebaMultipoleForce : public Force {
+class OPENMM_EXPORT AmoebaMultipoleForce : public Force {
 public:
  
     enum MultipoleAxisTypes { ZThenX, Bisector };
@@ -281,6 +284,15 @@ public:
         charge(charge), axisType(axisType), multipoleAtomId1(multipoleAtomId1), multipoleAtomId2(multipoleAtomId2),
         thole(thole), dampingFactor(dampingFactor), polarity(polarity) {
 
+       // only 'Z-then-X' or 'Bisector' currently handled
+
+       if( axisType != ZThenX && axisType != Bisector ){
+           std::stringstream buffer;
+           buffer << "AmoebaMultipoleForce: axis type=" << axisType;
+           buffer << " not currently handled - only axisTypes[ " << ZThenX << ", " << Bisector << "] (ZThenX, Bisector) currently handled .";
+           throw OpenMMException(buffer.str());
+       }
+
        covalentInfo.resize( CovalentEnd );
 
        molecularDipole.resize( 3 );

plugins/amoeba/openmmapi/src/AmoebaMultipoleForce.cpp

@@ -40,6 +40,7 @@ const int AmoebaMultipoleForce::CovalentDegrees[8] = { 1, 2, 3, 4, 0, 1, 2, 3 };
 
 
 AmoebaMultipoleForce::AmoebaMultipoleForce() {
+
     mutualInducedIterationMethod  = SOR;
     mutualInducedMaxIterations    = 60;
     mutualInducedTargetEpsilon    = 1.0e-06;

plugins/amoeba/platforms/cuda/src/kernels/AmoebaGpu.cpp

@@ -971,7 +971,8 @@ if( (index < DUMP_PARAMETERS || index > totalEntries - (DUMP_PARAMETERS + 1)) &&
     }
     
 #if 0
-    float* grids = (float*) malloc( totalGridEntries*sizeof( float ) );
+    std::vector<float> grids;
+    grids.resize( totalGridEntries );
     int index    = 0;
     for (unsigned int ii = 0; ii < floatGrids.size(); ii++) {
         for (unsigned int jj = 0; jj < floatGrids[ii].size(); jj++) {
@@ -1022,7 +1023,6 @@ if( (index < DUMP_PARAMETERS || index > totalEntries - (DUMP_PARAMETERS + 1)) &&
     if( !errors ){
         (void) fprintf( amoebaGpu->log, "No errors in grid readback\n" );
     }
-    free( grids );
 exit(0);
 #endif
 
@@ -1458,13 +1458,19 @@ void gpuSetAmoebaMultipoleParameters(amoebaGpuContext amoebaGpu, const std::vect
 #ifdef AMOEBA_DEBUG
     static unsigned int targetAtoms[2] = { 0, 700000 };
 
-    float* pScaleCheckSum = (float*) malloc( sizeof( float )*charges.size() );
-    float* dScaleCheckSum = (float*) malloc( sizeof( float )*charges.size() );
-    float* mScaleCheckSum = (float*) malloc( sizeof( float )*charges.size() );
+    std::vector<float> pScaleCheckSum;
+    pScaleCheckSum.resize( charges.size() );
 
-    memset( pScaleCheckSum, 0, charges.size()*sizeof( float ) );
-    memset( dScaleCheckSum, 0, charges.size()*sizeof( float ) );
-    memset( mScaleCheckSum, 0, charges.size()*sizeof( float ) );
+    std::vector<float> dScaleCheckSum;
+    dScaleCheckSum.resize( charges.size() );
+
+    std::vector<float> mScaleCheckSum;
+    mScaleCheckSum.resize( charges.size() );
+    for( unsigned int ii = 0; ii < charges.size(); ii++ ){
+        pScaleCheckSum[ii] = 0.0;
+        dScaleCheckSum[ii] = 0.0;
+        mScaleCheckSum[ii] = 0.0;
+    }
 
     amoebaGpu->pMapArray = (MapIntFloat**) malloc( sizeof( MapIntFloat* )*amoebaGpu->paddedNumberOfAtoms );
     amoebaGpu->dMapArray = (MapIntFloat**) malloc( sizeof( MapIntFloat* )*amoebaGpu->paddedNumberOfAtoms );
@@ -1502,9 +1508,9 @@ void gpuSetAmoebaMultipoleParameters(amoebaGpuContext amoebaGpu, const std::vect
     unsigned int quadrupoleIndex                                              = 0;
     unsigned int maxPrint                                                     = 5;
     
-    int* maxIndices = (int*) malloc( charges.size()*sizeof(int) );
+    std::vector<int> maxIndices;
     for( unsigned int ii = 0; ii < charges.size(); ii++ ){
-        maxIndices[ii]                                                        = ii;
+        maxIndices.push_back(ii);
         amoebaGpu->psMultipoleParticlesIdsAndAxisType->_pSysStream[0][ii].z   = ii;
     }
 
@@ -1788,9 +1794,6 @@ void gpuSetAmoebaMultipoleParameters(amoebaGpuContext amoebaGpu, const std::vect
             (void) fprintf( filePtr, "%6u %14.6e %14.6e\n", kk, pScaleCheckSum[kk], dScaleCheckSum[kk] );
         }
         (void) fclose( filePtr );
-        free( pScaleCheckSum );
-        free( dScaleCheckSum );
-        free( mScaleCheckSum );
         filePtr = fopen( "oldScaleMap.txt", "w" );
         for( unsigned int kk = 0; kk < charges.size(); kk++ ){
             MapIntFloat* pMap = amoebaGpu->pMapArray[kk];
@@ -1813,24 +1816,24 @@ void gpuSetAmoebaMultipoleParameters(amoebaGpuContext amoebaGpu, const std::vect
     for( unsigned int ii = 0; ii < charges.size(); ii++ ){
 
         // axis type & multipole particles ids
+
         int diff = maxIndices[ii] - amoebaGpu->psMultipoleParticlesIdsAndAxisType->_pSysStream[0][ii].z;
         if( diff > amoebaGpu->maxMapTorqueDifference ){
             amoebaGpu->maxMapTorqueDifference = diff;
         }
     }
+    amoebaGpu->maxMapTorqueDifference += 1;
 
     // set maxMapTorqueDifferencePow2 to smallest power of 2 greater than maxMapTorqueDifference
 
     float logDiff                                = logf( static_cast<float>(amoebaGpu->maxMapTorqueDifference) )/logf(2.0f);    
     int logDiffI                                 = static_cast<int>(logDiff) + 1;
     amoebaGpu->maxMapTorqueDifferencePow2        = static_cast<int>(powf( 2.0f, static_cast<float>(logDiffI) ));
-
+ 
     amoebaGpu->torqueMapForce                    = new CUDAStream<float>(amoebaGpu->paddedNumberOfAtoms*3*amoebaGpu->maxMapTorqueDifference, 1, "torqueMapForce");
     memset( amoebaGpu->torqueMapForce->_pSysStream[0], 0, amoebaGpu->paddedNumberOfAtoms*3*amoebaGpu->maxMapTorqueDifference*sizeof( float ) );
     amoebaGpu->torqueMapForce->Upload();
 
-    free( maxIndices );
-
     amoebaGpuBuildOutputBuffers( amoebaGpu );
     amoebaGpuBuildThreadBlockWorkList( amoebaGpu );
     amoebaGpuBuildScalingList( amoebaGpu );
@@ -3196,14 +3199,20 @@ static unsigned int targetAtoms[2] = { 0, 1};
 
     if( debugOn && amoebaGpu->log ){
 
-        float* pScaleCheckSum = (float*) malloc( sizeof( float )*paddedAtoms );
-        float* dScaleCheckSum = (float*) malloc( sizeof( float )*paddedAtoms );
-        float* mScaleCheckSum = (float*) malloc( sizeof( float )*paddedAtoms );
+        std::vector<float> pScaleCheckSum;
+        pScaleCheckSum.resize( paddedAtoms );
+    
+        std::vector<float> dScaleCheckSum;
+        dScaleCheckSum.resize( paddedAtoms );
+    
+        std::vector<float> mScaleCheckSum;
+        mScaleCheckSum.resize( paddedAtoms );
+        for( unsigned int ii = 0; ii < paddedAtoms; ii++ ){
+            pScaleCheckSum[ii] = 0.0;
+            dScaleCheckSum[ii] = 0.0;
+            mScaleCheckSum[ii] = 0.0;
+        }
     
-        memset( pScaleCheckSum, 0, paddedAtoms*sizeof( float ) );
-        memset( dScaleCheckSum, 0, paddedAtoms*sizeof( float ) );
-        memset( mScaleCheckSum, 0, paddedAtoms*sizeof( float ) );
-
         MapIntFloat** pMapArray = (MapIntFloat**) malloc( sizeof( MapIntFloat* )*paddedAtoms );
         MapIntFloat** dMapArray = (MapIntFloat**) malloc( sizeof( MapIntFloat* )*paddedAtoms );
         for( unsigned int ii = 0; ii < paddedAtoms; ii++ ){
@@ -3404,9 +3413,6 @@ tgx     = 0;
             (void) fprintf( filePtr, "%6u %14.6e %14.6e\n", kk, pScaleCheckSum[kk], dScaleCheckSum[kk] );
         }
         (void) fclose( filePtr );
-        free( pScaleCheckSum );
-        free( dScaleCheckSum );
-        free( mScaleCheckSum );
         filePtr = fopen( "newScaleMap.txt", "w" );
         //char buffer[1024];
 

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

@@ -963,4 +963,3 @@ void cudaComputeAmoebaElectrostatic( amoebaGpuContext amoebaGpu )
 
    // ---------------------------------------------------------------------------------------
 }
-

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

@@ -125,9 +125,9 @@ void amoebaMapTorqueToForce_kernel(
     norm                          = normVector3( vp );
 
     float dphi[3];
-    dphi[0]                       = vector[0][0]*torque[threadId*3] + vector[0][1]*torque[threadId*3+1] + vector[0][2]*torque[threadId*3+2];
-    dphi[1]                       = vector[1][0]*torque[threadId*3] + vector[1][1]*torque[threadId*3+1] + vector[1][2]*torque[threadId*3+2];
-    dphi[2]                       = vector[2][0]*torque[threadId*3] + vector[2][1]*torque[threadId*3+1] + vector[2][2]*torque[threadId*3+2];
+    dphi[0]                       = vector[U][0]*torque[threadId*3] + vector[U][1]*torque[threadId*3+1] + vector[U][2]*torque[threadId*3+2];
+    dphi[1]                       = vector[V][0]*torque[threadId*3] + vector[V][1]*torque[threadId*3+1] + vector[V][2]*torque[threadId*3+2];
+    dphi[2]                       = vector[W][0]*torque[threadId*3] + vector[W][1]*torque[threadId*3+1] + vector[W][2]*torque[threadId*3+2];
 
     // clamp c to interval [-1,1]
 

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

@@ -374,9 +374,7 @@ void kCalculateAmoebaMultipoleForces(amoebaGpuContext amoebaGpu, bool hasAmoebaG
 
     // check if induce dipole calculation converged -- abort if it did not
 
-    if( amoebaGpu->mutualInducedDone ){
-       //cudaComputeAmoebaElectrostatic( amoebaGpuContextGlobal );
-    } else {
+    if( amoebaGpu->mutualInducedDone == 0 ){
        (void) fprintf( amoebaGpu->log, "%s induced dipole calculation did not converge -- aborting!\n", methodName.c_str() );
        (void) fflush( amoebaGpu->log );
        exit(-1);

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

@@ -30,6 +30,7 @@
 #include "kernels/amoebaGpuTypes.h"
 #include "AmoebaCudaData.h"
 #include "openmm/LocalEnergyMinimizer.h"
+#include "../../../../../platforms//reference/src//SimTKUtilities/SimTKOpenMMUtilities.h"
 
 #include <exception>
 
@@ -2769,7 +2770,7 @@ static int readAmoebaWcaDispersionParameters( FILE* filePtr, MapStringInt& force
             (void) fprintf( log, "%8d %10.4f %10.4f", ii, radius, epsilon );
             if( ii < maxDispersionEnergyVector.size() ){
                 wcaDispersionForce->getMaximumDispersionEnergy( ii, maxDispersionEnergy );
-                maxDispersionEnergy /= CalToJoule;
+                if( useOpenMMUnits )maxDispersionEnergy /= CalToJoule;
                 double delta = fabs( maxDispersionEnergy - maxDispersionEnergyVector[ii] );
                 const char* error  = (delta > 1.0e-05) ? "XXX" : "";
                 (void) fprintf( log, " maxDispEDiff=%12.5e %14.7f %14.7f  %s",
@@ -2791,7 +2792,7 @@ static int readAmoebaWcaDispersionParameters( FILE* filePtr, MapStringInt& force
         for( unsigned int ii = 0; ii < arraySize && ii < maxDispersionEnergyVector.size(); ii++ ){
             double maxDispersionEnergy;
             wcaDispersionForce->getMaximumDispersionEnergy( ii, maxDispersionEnergy );
-            maxDispersionEnergy /= CalToJoule;
+            if( useOpenMMUnits )maxDispersionEnergy /= CalToJoule;
             double delta = fabs( maxDispersionEnergy - maxDispersionEnergyVector[ii] );
             if( delta > 1.0e-05 ){
                 (void) fprintf( log, " maxDispEDiff=%12.5e %14.7f %14.7f  XXX\n", delta, maxDispersionEnergy, maxDispersionEnergyVector[ii] );
@@ -3825,7 +3826,7 @@ void checkIntermediateMultipoleQuantities( Context* context, MapStringVectorOfVe
                                            int useOpenMMUnits, FILE* log ) {
 
 // ---------------------------------------------------------------------------------------
-/*
+
     static const std::string methodName      = "checkIntermediateMultipoleQuantities";
  
 // ---------------------------------------------------------------------------------------
@@ -3951,7 +3952,7 @@ void checkIntermediateMultipoleQuantities( Context* context, MapStringVectorOfVe
         unsigned int misses           = 0;
         double tolerance              = 1.0e-03;
         double dipoleConversion       = useOpenMMUnits ? 1.0/AngstromToNm : 1.0;
-        numberOfEntries      = expectedInducedDipoles.size() < numberOfEntries ? expectedInducedDipoles.size() : numberOfEntries;
+        numberOfEntries               = expectedInducedDipoles.size() < numberOfEntries ? expectedInducedDipoles.size() : numberOfEntries;
         for( unsigned int ii = 0; ii < numberOfEntries; ii++ ){
 
             std::vector<double> expectedInducedDipole = expectedInducedDipoles[ii];
@@ -3964,11 +3965,16 @@ void checkIntermediateMultipoleQuantities( Context* context, MapStringVectorOfVe
                 if( diff > 1.0e-04 ){
                     diff = 2.0*diff/(fabs( inducedDipoleValue ) + fabs( expectedInducedDipole[jj] ) );
                 }
+
+                int printDipole = 0;
                 if( diff > tolerance ){
                     misses++;
-                    if( misses == 1 ){
-                        (void) fprintf( log, "%s: induced dipole\n", methodName.c_str() );
-                    }
+                    printDipole = 1;
+                }
+                if( misses == 1 && printDipole ){
+                    (void) fprintf( log, "%s: induced dipoles\n", methodName.c_str() );
+                }
+                if( printDipole ){
                     if( !rowHit ){
                         (void) fprintf( log, "     Row %5u\n", ii );
                         rowHit = 1;
@@ -3987,7 +3993,7 @@ void checkIntermediateMultipoleQuantities( Context* context, MapStringVectorOfVe
         (void) fprintf( log, "Induced dipoles not available %s\n", e.what() ); 
         (void) fflush( log );
     }
-*/
+
 }
 
 void calculateBorn1( System& amoebaSystem, std::vector<Vec3>& tinkerCoordinates, FILE* log ) {
@@ -5235,20 +5241,45 @@ static void setVelocitiesBasedOnTemperature( const System& system, std::vector<V
 
 // ---------------------------------------------------------------------------------------
 
-   // set velocities based on temperature
-
-   temperature   *= BOLTZ;
-   double randMax = static_cast<double>(RAND_MAX);
-   randMax        = 1.0/randMax;
-   for( unsigned int ii = 0; ii < velocities.size(); ii++ ){
-      double velocityScale      = std::sqrt( temperature/system.getParticleMass(ii) );
-      randomValues[0]           = randMax*( (double) rand() );
-      randomValues[1]           = randMax*( (double) rand() );
-      randomValues[2]           = randMax*( (double) rand() );
-      velocities[ii]            = Vec3( randomValues[0]*velocityScale, randomValues[1]*velocityScale, randomValues[2]*velocityScale );
-   }
+    // set velocities based on temperature
 
-   return;
+    temperature          *= BOLTZ;
+    double kineticEnergy  = 0.0;
+    for( unsigned int ii = 0; ii < velocities.size(); ii++ ){
+        double mass               = system.getParticleMass(ii);
+        double velocityScale      = std::sqrt( temperature/mass );
+        randomValues[0]           = SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
+        randomValues[1]           = SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
+        randomValues[2]           = SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
+        velocities[ii]            = Vec3( randomValues[0]*velocityScale, randomValues[1]*velocityScale, randomValues[2]*velocityScale );
+        kineticEnergy            += mass*(velocities[ii][0]*velocities[ii][0] + velocities[ii][1]*velocities[ii][1] + velocities[ii][2]*velocities[ii][2]);
+    }
+ 
+    double degreesOfFreedom  = static_cast<double>(3*velocities.size() - system.getNumConstraints() - 3 );
+    double approximateT      = (kineticEnergy)/(degreesOfFreedom*BOLTZ);
+     if( approximateT > 0.0 ){
+        double scale = sqrt(temperature/approximateT);
+        for( unsigned int ii = 0; ii < velocities.size(); ii++ ){
+            velocities[ii][0] *= scale;
+            velocities[ii][1] *= scale;
+            velocities[ii][2] *= scale;
+        }
+     }
+ 
+    if( log ){
+        double finalKineticEnergy  = 0.0;
+        for( unsigned int ii = 0; ii < velocities.size(); ii++ ){
+            double mass             = system.getParticleMass(ii);
+            finalKineticEnergy     += mass*(velocities[ii][0]*velocities[ii][0] + velocities[ii][1]*velocities[ii][1] + velocities[ii][2]*velocities[ii][2]);
+        }
+    
+        (void) fprintf( log, "%s KE=%15.7e approximateT=%15.7e desiredT=%15.7e dof=%12.3f final KE=%12.3e\n",
+                        methodName.c_str(), kineticEnergy, approximateT, temperature/BOLTZ,
+                        degreesOfFreedom, 0.5*finalKineticEnergy );
+    }
+ 
+ 
+    return;
 }
 
 /** 
@@ -5427,7 +5458,7 @@ double getEnergyDrift( std::vector<double>& totalEnergyArray, std::vector<double
  
     std::vector<double> kineticEnergyStatistics;
     getStatistics( kineticEnergyArray, kineticEnergyStatistics );
-    double temperature  = kineticEnergyStatistics[0]/(degreesOfFreedom*BOLTZ);
+    double temperature  = 2.0*kineticEnergyStatistics[0]/(degreesOfFreedom*BOLTZ);
     double kT           = temperature*BOLTZ;
  
     // compute stddev in units of kT/dof/ns
@@ -5503,6 +5534,15 @@ void testEnergyConservation( std::string parameterFileName, MapStringInt& forceM
         writeIntermediateStateFile( *context, intermediateStateFile, log );
     }
  
+    System& system                     = context->getSystem();
+    int numberOfAtoms                  = system.getNumParticles();
+ 
+    std::vector<Vec3> velocities; 
+    velocities.resize( numberOfAtoms );
+    double initialT = 300.0;
+    setVelocitiesBasedOnTemperature( system, velocities, initialT, log );
+    context->setVelocities(velocities);
+
     // energy minimize
 
     setIntFromMap( inputArgumentMap, "energyMinimize",          energyMinimize );
@@ -5535,14 +5575,6 @@ void testEnergyConservation( std::string parameterFileName, MapStringInt& forceM
  
     // set velocities based on temperature
  
-    System& system                     = context->getSystem();
-    int numberOfAtoms                  = system.getNumParticles();
-    std::vector<Vec3> velocities; 
-    velocities.resize( numberOfAtoms );
-    double initialT = 300.0;
-    setVelocitiesBasedOnTemperature( system, velocities, initialT, log );
-    context->setVelocities(velocities);
- 
     // get integrator 
  
     Integrator& integrator                                  = context->getIntegrator();

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

@@ -106,11 +106,6 @@ static std::string MUTUAL_INDUCED_TARGET_EPSILON                      = "mutualI
 #define SumIndex                                          12
 #define AmoebaLastIndex                                   13
 
-#define BOLTZMANN  (1.380658e-23)                         /* (J/K) */
-#define AVOGADRO   (6.0221367e23)                         /* ()    */
-#define RGAS       (BOLTZMANN*AVOGADRO)                   /* (J/(mol K))  */
-#define BOLTZ      (RGAS/1.0e+03)                         /* (kJ/(mol K)) */
-
 #define AngstromToNm 0.1
 #define CalToJoule   4.184