Minor edits; force duplicationFactor to 1 if no. particles < 65; for system w/...

Minor edits; force duplicationFactor to 1 if no. particles < 65; for system w/ 32 atoms, NAN were appearing in nonbonded force term; plan to fix later

platforms/brook/src/BrookCalcPeriodicTorsionForceKernel.cpp

@@ -167,7 +167,7 @@ void BrookCalcPeriodicTorsionForceKernel::initialize( const System& system, cons
 
    if( log ){
       std::string contents = _brookBondParameters->getContentsString( ); 
-      (void) fprintf( log, "%s brookGbsa::contents\n%s", methodName.c_str(), contents.c_str() );
+      (void) fprintf( log, "%s contents\n%s", methodName.c_str(), contents.c_str() );
       (void) fflush( log );
    }
 

platforms/brook/src/BrookNonBonded.cpp

@@ -1346,22 +1346,34 @@ nonbondedForceStreams[3]->fillWithValue( &zerof );
    if( printOn ){
    //static int step = 0;
    //if( step++ < 1 ){
-      //FILE* log = getLog();
-      FILE* log = stderr;
+      FILE* log = getLog();
+      //FILE* log = stderr;
       (void) fprintf( log, "%s\n", methodName.c_str() ); (void) fflush( log );
-      (void) fprintf( log, "\nPost knbforce_CDLJ4: particles=%6d ceiling=%3d dupFac=%3d", getNumberOfParticles(),  
+      (void) fprintf( log, "Post knbforce_CDLJ4: particles=%6d ceiling=%3d dupFac=%3d\n", getNumberOfParticles(),  
                                                                                           getParticleSizeCeiling(),
                                                                                           getDuplicationFactor()  );
 
-      (void) fprintf( log, "\n                      hght=%6d   width=%3d   jWid=%3d", getParticleStreamHeight( ),
-                                                                                      getParticleStreamWidth( ),
-                                                                                      getJStreamWidth( ) );
-      (void) fprintf( log, "\n                      pFrc=%6d     eps=%12.5e\n",       getPartialForceStreamWidth( ), epsfac );
+      (void) fprintf( log, "                        hght=%6d   width=%3d   jWid=%3d\n", getParticleStreamHeight( ),
+                                                                                        getParticleStreamWidth( ),
+                                                                                        getJStreamWidth( ) );
+      (void) fprintf( log, "                        pFrc=%6d     eps=%12.5e\n",         getPartialForceStreamWidth( ), epsfac );
 
-      (void) fprintf( log, "\nFinal NB & bonded forces" );
+      (void) fprintf( log, "%s Final NB & bonded forces\n", methodName.c_str() );
+      (void) fprintf( log, "Forces here should be zero (pregather)\n" );
       BrookStreamInternal* brookStreamInternalF   = forceStream.getBrookStreamInternal();
       brookStreamInternalF->printToFile( log );
  
+      (void) fprintf( log, "Coords\n" );
+      brookStreamInternalF   = positionStream.getBrookStreamInternal();
+      brookStreamInternalF->printToFile( log );
+/*
+void* dataArrayV = brookStreamInternalF->getData( 1 );
+float* ddd = (float*) dataArrayV;
+int idx1 = 0;
+int idx2 = 3;
+      (void) fprintf( log, "%d %d %.5e\n", idx1/3, idx2/3, computeDistance( idx1, idx2, ddd ) );
+*/
+ 
       (void) fprintf( log, "\nOuterVdwStreamd\n" );
       getOuterVdwStream()->printToFile( log );
 
@@ -1413,3 +1425,32 @@ nonbondedForceStreams[3]->fillWithValue( &zerof );
 
    return;
 }
+
+/*  
+ * Utility to compute distances between two points
+ *
+ * @param idx1       index of point1 into array 
+ * @param idx2       index of point2 into array 
+ * @param coords     array of points p1={coord[0], coord[1], coord[2] },
+ *                                   p2={coord[3], coord[4], coord[5] },
+ *                                   p3={coord[6], coord[7], coord[8] }
+ *
+ * @return distance
+ *
+ **/
+    
+
+float BrookNonBonded::_computeDistance( int idx1, int idx2, float* coords ){
+
+// ---------------------------------------------------------------------------------------
+
+   //static const std::string methodName      = "BrookNonBonded::_computeDistance";
+
+   float d1 = sqrtf( 
+                      (coords[idx1]-coords[idx2])*(coords[idx1]-coords[idx2]) +
+                      (coords[idx1+1]-coords[idx2+1])*(coords[idx1+1]-coords[idx2+1]) +
+                      (coords[idx1+2]-coords[idx2+2])*(coords[idx1+2]-coords[idx2+2]) );
+   return d1;
+
+// ---------------------------------------------------------------------------------------
+}

platforms/brook/src/BrookNonBonded.h

@@ -548,6 +548,20 @@ class BrookNonBonded : public BrookCommon {
       
       int _initializePartialForceStreamSize( int particleStreamSize, int particleStreamWidth );
       
+      /* 
+       * Utility to compute distances between two points
+       *
+       * @param idx1       index of point1 into array 
+       * @param idx2       index of point2 into array 
+       * @param coord      array of points p1={coord[0], coord[1], coord[2] },
+       *                                   p2={coord[3], coord[4], coord[5] },
+       *                                   p3={coord[6], coord[7], coord[8] }
+       *
+       * @return distance
+       *
+       * */
+      
+      float _computeDistance( int idx1, int idx2, float* coord );
 };
 
 } // namespace OpenMM

platforms/brook/src/BrookPlatform.cpp

@@ -208,7 +208,6 @@ BrookPlatform::BrookPlatform( ){
    _particleStreamWidth  = DefaultParticleStreamWidth;
    _minSuggestedThreads  = -1;
    _log                  = NULL;
-//_log                  = stderr;
 
    // get Brook runtime
 
@@ -386,13 +385,18 @@ int BrookPlatform::getDuplicationFactor( int numberOfParticles ) const {
 
 // ---------------------------------------------------------------------------------------
 
-   // static const std::string methodName      = "BrookPlatform::getDuplicationFactor";
+   static const std::string methodName      = "BrookPlatform::getDuplicationFactor";
 
 // ---------------------------------------------------------------------------------------
 
    // default value
 
    int duplicationFactor = 4;
+   if( numberOfParticles < 65 ){
+      duplicationFactor = 1;
+      (void) fprintf( stderr, "%s forcing duplicationFactor=%d since number of particles=%d < 65\n",
+                      methodName.c_str(), duplicationFactor, numberOfParticles );
+   }
 
    // set only if _minSuggestedThreads is available from board
 /*

platforms/brook/src/OpenMMBrookInterface.cpp

@@ -68,6 +68,7 @@ OpenMMBrookInterface::OpenMMBrookInterface( int streamWidth, int duplicationFact
    _forces                                  = NULL;
 
    _log                                     = NULL;
+//_log = stderr;
 
    _particleStreamSize                      = -1;