Mods

platforms/brook/src/BrookBonded.cpp

@@ -787,18 +787,23 @@ int BrookBonded::_addPTorsions( int *nbondeds, int *particles, BrookOpenMMFloat*
 
 // ---------------------------------------------------------------------------------------
 
-   static const std::string methodName = "BrookBonded::_addPTorsion";
-   static const int debug              = 0;
+   static const BrookOpenMMFloat DEG2RAD =  (BrookOpenMMFloat) (M_PI/180.0);
+
+   static const std::string methodName   = "BrookBonded::_addPTorsion";
+   static const int debug                = 0;
 
 // ---------------------------------------------------------------------------------------
 
-   if( debug && getLog() ){
-      (void) fprintf( getLog(), "%s npdih=%d\n", methodName.c_str(), periodicTorsionIndices.size() );
+   FILE* log = getLog();
+   if( debug ){
+      log = log ? log : stderr;
+      (void) fprintf( log, "%s npdih=%d\n", methodName.c_str(), periodicTorsionIndices.size() );
+      (void) fflush( log );
    }
 
    for( unsigned int ii = 0; ii < periodicTorsionIndices.size(); ii++ ){
 
-      vector<int> particlesIndices      = periodicTorsionIndices[ii];
+      vector<int> particlesIndices  = periodicTorsionIndices[ii];
       vector<double> pTParameters   = periodicTorsionParameters[ii];
 
       int index = 0;
@@ -821,12 +826,13 @@ int BrookBonded::_addPTorsions( int *nbondeds, int *particles, BrookOpenMMFloat*
       // note: parameters 0 & 2 switched
 
       PARAMS( ibonded, 1, 1 ) = (BrookOpenMMFloat) pTParameters[2];
-      PARAMS( ibonded, 1, 2 ) = (BrookOpenMMFloat) pTParameters[1];
+      PARAMS( ibonded, 1, 2 ) = (BrookOpenMMFloat) pTParameters[1]*DEG2RAD;
       PARAMS( ibonded, 1, 3 ) = (BrookOpenMMFloat) pTParameters[0];
 
-      if( debug && getLog() ){
-         (void) fprintf( getLog(), "   %d [%d %d %d %d] %.3e %.3e %.3e\n", ibonded, i, j, k, l,
+      if( debug ){
+         (void) fprintf( log, "   %d [%d %d %d %d] %.3e %.3e %.3e\n", ibonded, i, j, k, l,
                          pTParameters[0], pTParameters[1], pTParameters[2] );
+        (void) fflush( log );
       }
    }
 
@@ -1065,12 +1071,12 @@ int BrookBonded::_loadInvMaps( int nbondeds, int nparticles, int *particles, int
 // ---------------------------------------------------------------------------------------
 
    static const std::string methodName = "BrookBonded::_loadInvMaps";
-   static       int PrintOn            = 0;
+   static       int printOn            = 0;
    double dangleValue                  = 0.0;
 
 // ---------------------------------------------------------------------------------------
 
-   PrintOn = (PrintOn && getLog()) ? 1 : 0;
+   printOn = (printOn && getLog()) ? 1 : 0;
 
    // get particle stream size
 
@@ -1114,7 +1120,7 @@ int BrookBonded::_loadInvMaps( int nbondeds, int nparticles, int *particles, int
       }
    }
 
-   if( PrintOn ){
+   if( printOn ){
       (void) fprintf( getLog(), "%s force stream strms=%d nbondeds=%d max counts=[%d %d %d %d] strSz&Wd=%d %d\n", methodName.c_str(), getNumberOfForceStreams(),
                       nbondeds, getMaxInverseMapStreamCount(0), getMaxInverseMapStreamCount(1), getMaxInverseMapStreamCount(2), getMaxInverseMapStreamCount(3),
                       particleStreamSize, particleStreamWidth );
@@ -1127,21 +1133,26 @@ int BrookBonded::_loadInvMaps( int nbondeds, int nparticles, int *particles, int
       for( int jj = 0; jj < 4*getMaxInverseMapStreamCount()*particleStreamSize; jj++ ){
          block[jj] = -1.0f;
       }
-//(void) fprintf( getLog(), "%s _gpuCalcInvMap %d getInverseMapStreamCount=%d\n", methodName.c_str(), ii, getInverseMapStreamCount( ii ) ); (void) fflush( getLog() );
+
+(void) fprintf( stderr, "%s _gpuCalcInvMap %d getInverseMapStreamCount=%d\n", methodName.c_str(), ii, getInverseMapStreamCount( ii ) ); (void) fflush( getLog() );
 
       _gpuCalcInvMap( ii, 4, nbondeds, nparticles, particles, getMaxInverseMapStreamCount( ii ), counts, invmaps, &(_inverseMapStreamCount[ii]) );
-//gpuPrintInvMaps( _inverseMapStreamCount[ii], nparticles, counts, invmaps, getLog() );
+
+_gpuPrintInvMaps( _inverseMapStreamCount[ii], nparticles, counts, invmaps, stderr );
+
       _validateInverseMapStreamCount( ii, _inverseMapStreamCount[ii] ); 
       for( int jj = 0; jj < _inverseMapStreamCount[ii]; jj++ ){
          _inverseStreamMaps[ii][jj]->loadFromArray( invmaps[jj] );
 
-         if( PrintOn ){
-            (void) fprintf( getLog(), "%s inverseMap stream strms=%d count=%d index=%d %d InverseMapStreamCount[ii]=%d max=%d\n",
+         //if( printOn ){
+         if( 1 ){
+            FILE* log = stderr;
+            (void) fprintf( log, "%s inverseMap stream strms=%d count=%d index=%d %d InverseMapStreamCount[ii]=%d max=%d\n",
                             methodName.c_str(), getNumberOfForceStreams(), _inverseMapStreamCount[ii], ii, jj,
                             getInverseMapStreamCount( ii ), getMaxInverseMapStreamCount( ii ) );
 
             for( int kk = 0; kk < particleStreamSize; kk++ ){
-               (void) fprintf( getLog(), "%8d [ %.1f %.1f %.1f %.1f]\n", kk, invmaps[jj][kk].x, invmaps[jj][kk].y, invmaps[jj][kk].z, invmaps[jj][kk].w  );
+               (void) fprintf( log, "%8d [ %.1f %.1f %.1f %.1f]\n", kk, invmaps[jj][kk].x, invmaps[jj][kk].y, invmaps[jj][kk].z, invmaps[jj][kk].w  );
             }
          }
 
@@ -1163,7 +1174,7 @@ int BrookBonded::_loadInvMaps( int nbondeds, int nparticles, int *particles, int
 
    // diagnostics 
 
-   if( PrintOn ){
+   if( printOn ){
       (void) fprintf( getLog(), "%s done\n", methodName.c_str() );
       (void) fflush( getLog() );
    }
@@ -1216,19 +1227,22 @@ int BrookBonded::setup( int numberOfParticles,
 // ---------------------------------------------------------------------------------------
 
    static const std::string methodName = "BrookBonded::setup";
-   static       int PrintOn            = 0;
+   static       int printOn            = 0;
    double dangleValue                  = 0.0;
+   FILE* log;
 
 // ---------------------------------------------------------------------------------------
 
-   PrintOn = (PrintOn && getLog()) ? 1 : 0;
+//setLog( stderr );
+   printOn = (printOn && getLog()) ? printOn : 0;
 
-   if( PrintOn ){
-      (void) fprintf( getLog(), "%s particles=%d\n   [%p %p %p %p %p] (bond, angle, pd, rb, 14)\n"
+   if( printOn ){
+      log = getLog();
+      (void) fprintf( log, "%s particles=%d\n   [%p %p %p %p %p] (bond, angle, pd, rb, 14)\n"
                       "StreamW=%d StreamSz=%d\n", methodName.c_str(), numberOfParticles, harmonicBondBrookBondParameters,
                       harmonicAngleBrookBondParameters,
                       periodicTorsionBrookBondParameters, rbTorsionBrookBondParameters, nonBonded14ForceParameters,
-                      particleStreamWidth, particleStreamSize ); fflush( getLog() );
+                      particleStreamWidth, particleStreamSize ); fflush( log );
    }
 
    _numberOfParticles = numberOfParticles;
@@ -1309,16 +1323,16 @@ int BrookBonded::setup( int numberOfParticles,
 
       // return if no bonds
 
-      if( PrintOn || getLog() ){
+      if( printOn || getLog() ){
          (void) fprintf( getLog(), "%s WARNING: particles=%d number of bonds=%d maxBonds=%d\n", methodName.c_str(), numberOfParticles, nbondeds, maxBonds );
       }
 
       _setupCompleted = 1;
       return DefaultReturnValue;
 
-   } else if( PrintOn ){
-      (void) fprintf( getLog(), "%s particles=%d number of bonds=%d maxBonds=%d\n", methodName.c_str(), numberOfParticles, nbondeds, maxBonds );
-      (void) fflush( getLog() );
+   } else if( printOn ){
+      (void) fprintf( log, "%s particles=%d number of bonds=%d maxBonds=%d\n", methodName.c_str(), numberOfParticles, nbondeds, maxBonds );
+      (void) fflush( log );
    }
 
 // ---------------------------------------------------------------------------------------
@@ -1366,12 +1380,13 @@ int BrookBonded::setup( int numberOfParticles,
 
    // debug stuff
 
-   if( PrintOn ){
+   if( printOn ){
 
-      (void) fprintf( getLog(), "%s nbondeds=%d strDim [%d %d ] sz=%d\n", methodName.c_str(), nbondeds,
+      (void) fprintf( log, "%s nbondeds=%d strDim [%d %d ] sz=%d\n", methodName.c_str(), nbondeds,
                       _particleIndicesStream->getStreamWidth(),
                       _particleIndicesStream->getStreamHeight(), 
                       _particleIndicesStream->getStreamSize() );
+      (void) fflush( log );
 
       int kIndex = 0;
       int jIndex = 1;
@@ -1387,11 +1402,11 @@ int BrookBonded::setup( int numberOfParticles,
        * float4 parm4   (x0, k, sig14, eps14) bond k-l, i-l
        */
       
-      (void) fprintf( getLog(), "\nParams\n" );
+      (void) fprintf( log, "\nParams\n" );
       int index = 0; 
       for( int ii = 0; ii < 4*nbondeds; ii += 4, index++ ){
          //  #define PARAMS(X,Y,Z) (params[(Y)][4*(X) + Z])
-         (void) fprintf( getLog(), "\n%4d   [%4d %4d %4d %4d]\n"
+         (void) fprintf( log, "\n%4d   [%4d %4d %4d %4d]\n"
                               "     rb[%10.6f %10.6f %10.6f %10.6f %10.6f]\n"
                               "    phi[%10.6f %6.1f %5.1f]\n"
                               "    ang[%6.3f %8.2f] [%6.3f %8.2f]\n"
@@ -1408,8 +1423,8 @@ int BrookBonded::setup( int numberOfParticles,
 
    // load inverse maps to streams
 
-   if( PrintOn ){
-      (void) fprintf( getLog(), "\n_loadInvMaps %d %d %d %d\n", nbondeds, getNumberOfParticles(), particleStreamWidth, particleStreamSize ); (void) fflush( getLog() );
+   if( printOn ){
+      (void) fprintf( log, "\n_loadInvMaps %d %d %d %d\n", nbondeds, getNumberOfParticles(), particleStreamWidth, particleStreamSize ); (void) fflush( getLog() );
    }
 
    _loadInvMaps( nbondeds, getNumberOfParticles(), particles, particleStreamWidth, particleStreamSize );
@@ -1572,7 +1587,7 @@ int BrookBonded::_gpuCalcInvMap( int posflag, int niparticles, int nints, int np
    //char value[MAX_LINE_CHARS];
    static const char* Set                   = "Set";
    static const char* NotSet                = "Not set";
-   static const int PrintOn                 = 0;
+   static const int printOn                 = 0;
 
    int particleRange[2]                     = { 90000000, -90000000 };
    int mapnumRange[2]                       = { 90000000, -90000000 };
@@ -1593,7 +1608,7 @@ int BrookBonded::_gpuCalcInvMap( int posflag, int niparticles, int nints, int np
 
    //Now note down the positions where each particle occurs
 
-   if( PrintOn && getLog() ){
+   if( printOn && getLog() ){
       (void) fprintf( getLog(), "%s: pos=%d ni=%d nints=%d nparticles=%d nmaps=<%d>\n", methodName.c_str(), posflag, niparticles, nints, nparticles, nmaps ); 
       (void) fflush( getLog() );
    }
@@ -1615,7 +1630,7 @@ if( particle > particleRange[1] ){
 }
       //Check to make sure we're inside the limits
       if ( counts[particle] > nmaps * 4 ){
-         if( PrintOn && getLog() ){
+         if( printOn && getLog() ){
             (void) fprintf( getLog(), "%s Particle %d has too many proper torsions (%d, max %d)\n",
                             methodName.c_str(), particle, counts[particle], nmaps*4 );
             (void) fflush( getLog() );
@@ -1644,7 +1659,7 @@ if( particle > particleRange[1] ){
          case 2: invmaps[mapnum][particle].z = (float) i; break;
          case 3: invmaps[mapnum][particle].w = (float) i; break;
          default:
-            if( PrintOn && getLog() ){
+            if( printOn && getLog() ){
                (void) fprintf( getLog(), "mapcomp %d invalid -- impossible!\n", mapcomp );
                (void) fflush( getLog() );
             }
@@ -1669,7 +1684,7 @@ if( mapnum > mapnumRange[1] ){
 
    (*nimaps)++;
 
-if( PrintOn && getLog() ){
+if( printOn && getLog() ){
    (void) fprintf( getLog(), "%s mnmaps=%d Ranges: particle [%d %d] mapnum [%d %d]\n",
                    methodName.c_str(), *nimaps, particleRange[0], particleRange[1], mapnumRange[0], mapnumRange[1] );
    (void) fflush( getLog() );
@@ -1814,7 +1829,7 @@ void BrookBonded::computeForces( BrookStreamImpl& positionStream, BrookStreamImp
 
    static const std::string methodName      = "BrookBonded::computeForces";
 
-   static int PrintOn                       = 1;
+   static int printOn                       = 0;
 
    static const int I_Stream                = 0;
    static const int J_Stream                = 1;
@@ -1824,13 +1839,9 @@ void BrookBonded::computeForces( BrookStreamImpl& positionStream, BrookStreamImp
    static const int MaxErrorMessages        = 2;
    static       int ErrorMessages           = 0;
 
-   static const float4 dummyParameters( 0.0, 0.0, 0.0, 0.0 );
-
-   // static const int debug                   = 1;
-
 // ---------------------------------------------------------------------------------------
 
-   PrintOn = (PrintOn && getLog()) ? 1 : 0;
+   printOn = (printOn && getLog()) ? printOn : 0;
 
    // bonded
 
@@ -1850,7 +1861,6 @@ void BrookBonded::computeForces( BrookStreamImpl& positionStream, BrookStreamImp
 
    kbonded_CDLJ( epsfac, 
                  (float) bondedForceStreams[0]->getStreamWidth(),
-                 dummyParameters,
                  positionStream.getBrookStream(),
                  getChargeStream()->getBrookStream(),
                  getParticleIndicesStream()->getBrookStream(),
@@ -1867,10 +1877,9 @@ void BrookBonded::computeForces( BrookStreamImpl& positionStream, BrookStreamImp
 
    // diagnostics
 
-   //if( 1 && PrintOn ){
-   if( 1 ){
-      //FILE* log = getLog();
-      FILE* log = stderr;
+   if( printOn ){
+      FILE* log = getLog();
+      //FILE* log = stderr;
 
       int countPrintInvMap[4] = { 3, 5, 2, 4 }; 
 
@@ -2015,8 +2024,7 @@ void BrookBonded::computeForces( BrookStreamImpl& positionStream, BrookStreamImp
 
    // diagnostics
 
-   //if( 1 && PrintOn ){
-   if( 1 ){
+   if( printOn ){
 
       FILE* log =  getLog() ?  getLog() : stderr;
       (void) fprintf( log, "\nPost 3_4/3_5 && NB forces" );
@@ -2086,8 +2094,7 @@ void BrookBonded::computeForces( BrookStreamImpl& positionStream, BrookStreamImp
 
    // diagnostics
 
-   //if( 1 && PrintOn ){
-   if( 1 ){
+   if( printOn ){
 
       FILE* log =  getLog() ?  getLog() : stderr;
       (void) fprintf( log, "\nFinal NB & bonded forces" );

platforms/brook/src/BrookFloatStreamInternal.cpp

@@ -59,6 +59,9 @@ BrookFloatStreamInternal::BrookFloatStreamInternal( const std::string& name, int
 
 // ---------------------------------------------------------------------------------------
 
+//fprintf( stderr,"%s %s\n", methodName.c_str(), getName().c_str() );
+//fflush( stderr );
+
    // set base type (currently only FLOAT supported)
 
    switch( type ){

platforms/brook/src/BrookNonBonded.cpp

@@ -1329,7 +1329,8 @@ nonbondedForceStreams[3]->fillWithValue( &zerof );
    // diagnostics
 
    //if( 1 && PrintOn ){
-   if( 1 ){
+static int step = 0;
+   if( step++ < 1 ){
       //FILE* log = getLog();
       FILE* log = stderr;
       (void) fprintf( log, "%s\n", methodName.c_str() ); (void) fflush( log );
@@ -1386,7 +1387,8 @@ nonbondedForceStreams[3]->fillWithValue( &zerof );
 
    // diagnostics
 
-   if( PrintOn ){
+   //if( PrintOn ){
+   if( 1 ){
 
       //FILE* log = getLog();
       FILE* log = stderr;

platforms/brook/src/OpenMMBrookInterface.cpp

@@ -34,8 +34,6 @@
 #include "OpenMMException.h"
 #include <sstream>
 
-// used for energy calculartion
-
 #include "LangevinIntegrator.h"
 #include "ReferencePlatform.h"
 #include "internal/OpenMMContextImpl.h"
@@ -533,23 +531,28 @@ void OpenMMBrookInterface::computeForces( OpenMMContextImpl& context ){
 
    static const std::string methodName      = "OpenMMBrookInterface::executeForces";
 
-   static       int PrintOn                 = 1;
+   static       int printOn                 = 0;
    static const int MaxErrorMessages        = 2;
    static       int ErrorMessages           = 0;
-
-   static const float4 dummyParameters( 0.0, 0.0, 0.0, 0.0 );
-
-   // static const int debug                   = 1;
+   FILE* log;
 
 // ---------------------------------------------------------------------------------------
 
-   PrintOn = (PrintOn && getLog()) ? 1 : 0;
-
-   // nonbonded forces
+//setLog( stderr );
+   printOn = (printOn && getLog()) ? printOn : 0;
 
    BrookStreamImpl* positions = getParticlePositions();
    BrookStreamImpl* forces    = getParticleForces();
 
+   // info
+
+   //if( printOn > 1 ){
+   if( 1 ){
+      printForcesToFile( context );
+   }
+
+   // nonbonded forces
+
    if( _brookNonBonded.isActive() ){
       _brookNonBonded.computeForces( *positions, *forces );
    }
@@ -558,9 +561,8 @@ void OpenMMBrookInterface::computeForces( OpenMMContextImpl& context ){
 
    // bonded forces
 
-   //if( PrintOn ){
-   if( 1 ){
-      FILE* log = stderr;
+   if( printOn ){
+      log = getLog();
       (void) fprintf( log, "%s done nonbonded: bonded=%d completed=%d\n", methodName.c_str(),
                       _brookBonded.isActive(), _brookBonded.isSetupCompleted() ); 
       (void) fflush( log );
@@ -587,13 +589,11 @@ void OpenMMBrookInterface::computeForces( OpenMMContextImpl& context ){
    
       // diagnostics
    
-      //if( 1 && PrintOn ){
-      if( 1 ){
+      if( printOn ){
          FILE* log = stderr;
          static int step           = 0;
          static const int stopStep = 10;
-         (void) fprintf( log, "%s done bonded computeForces\n", methodName.c_str(),
-                         _brookBonded.isActive(), _brookBonded.isSetupCompleted() ); 
+         (void) fprintf( log, "%s done bonded computeForces\n", methodName.c_str() );
          (void) fflush( log );
    
    /*
@@ -625,6 +625,120 @@ void OpenMMBrookInterface::computeForces( OpenMMContextImpl& context ){
    // ---------------------------------------------------------------------------------------
 }
 
+/** 
+ * Print forces to file 
+ * 
+ * @param context     context
+ *
+ */
+
+void OpenMMBrookInterface::printForcesToFile( OpenMMContextImpl& context ){
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName      = "OpenMMBrookInterface::printForcesToFile";
+   float zero                               = 0.0f;
+   static int step                          = 0;
+
+// ---------------------------------------------------------------------------------------
+
+   // first step only?
+
+   if( step > 0 ){
+      return;
+   }
+
+   std::stringstream fileNameBaseS;
+   //fileNameBase << "Brook_" << context.getTime() << "_"; 
+   fileNameBaseS << "Brook_" << step++ << "_"; 
+   std::string fileNameBase = fileNameBaseS.str();
+
+   // create vector of streams for output
+
+   BrookStreamImpl* positions = getParticlePositions();
+   BrookStreamImpl* forces    = getParticleForces();
+
+   std::vector<BrookStreamImpl*> streams;
+   streams.push_back( positions );
+   streams.push_back( forces );
+
+// ---------------------------------------------------------------------------------------
+
+   // nonbonded forces
+
+   if( _brookNonBonded.isActive() ){
+      forces->fillWithValue( &zero );
+      _brookNonBonded.computeForces( *positions, *forces );
+      std::string fileName = fileNameBase + "NonBonded.txt";
+      printStreamsToFile( fileName, streams );
+   }
+
+// ---------------------------------------------------------------------------------------
+
+   // bonded forces
+
+   if( _brookBonded.isActive() ){
+
+      // perform setup first time through
+
+      if( _brookBonded.isSetupCompleted() == 0 ){
+
+         _brookBonded.setup( getNumberOfParticles(), getHarmonicBondForceParameters(), getHarmonicAngleForceParameters(),
+                             getPeriodicTorsionForceParameters(), getRBTorsionForceParameters(),
+                             getNonBonded14ForceParameters(),
+                             getParticleStreamWidth(), getParticleStreamSize() );
+      }
+
+      forces->fillWithValue( &zero );
+      _brookBonded.computeForces( *positions, *forces );
+      std::string fileName = fileNameBase + "Bonded.txt";
+      printStreamsToFile( fileName, streams );
+   
+   }
+
+// ---------------------------------------------------------------------------------------
+
+   // GBSA OBC forces
+
+   if( _brookGbsa.isActive() ){
+
+      forces->fillWithValue( &zero );
+
+      _brookGbsa.computeForces( *positions, *forces );
+
+      std::string fileName = fileNameBase + "Obc.txt";
+      printStreamsToFile( fileName, streams );
+
+   }
+   forces->fillWithValue( &zero );
+
+// ---------------------------------------------------------------------------------------
+
+   // all forces
+
+   if( 1 ){
+
+      if( _brookNonBonded.isActive() ){
+         _brookNonBonded.computeForces( *positions, *forces );
+      }
+
+      if( _brookBonded.isActive() ){
+         _brookBonded.computeForces( *positions, *forces );
+      }
+
+      if( _brookGbsa.isActive() ){
+         _brookGbsa.computeForces( *positions, *forces );
+      }
+
+      std::string fileName = fileNameBase + "AllF.txt";
+      printStreamsToFile( fileName, streams );
+   }
+
+   forces->fillWithValue( &zero );
+
+   // ---------------------------------------------------------------------------------------
+}
+
 /** 
  * Compute energy
  * 
@@ -653,7 +767,7 @@ double OpenMMBrookInterface::computeEnergy( OpenMMContextImpl& context, System&
    positions.saveToArray(posData);
    vector<Vec3> pos(positions.getSize());
 
-   for( int ii = 0; ii < pos.size(); ii++ ){
+   for( unsigned int ii = 0; ii < pos.size(); ii++ ){
       pos[ii] = Vec3(posData[3*ii], posData[3*ii+1], posData[3*ii+2]);
    }
    delete[] posData;
@@ -662,3 +776,113 @@ double OpenMMBrookInterface::computeEnergy( OpenMMContextImpl& context, System&
    return refContext.getState(State::Energy).getPotentialEnergy();
 
 }
+
+/* 
+ * Print contents of object to file
+ *
+ * @param fileName     file name
+ * @param streams      streams to print
+ *
+ * @return DefaultReturnValue
+ *
+ * */
+
+int OpenMMBrookInterface::printStreamsToFile( std::string fileName, std::vector<BrookStreamImpl*>& streams ){
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName      = "OpenMMBrookInterface::printStreamsToFile";
+
+// ---------------------------------------------------------------------------------------
+
+   FILE* filePtr = fopen( fileName.c_str(), "w" );
+   if( !filePtr ){
+      (void) fprintf( stderr, "%s coud not open file=<%s>\n", methodName.c_str(), fileName.c_str() );
+      (void) fflush( stderr );
+      return ErrorReturnValue;
+   }
+
+   // gather arrays, widths for eah stream, and set index for each stream
+   // also set minimum of stream sizes
+
+   int minIndex    = 10000000;
+   float** arrays  = new float*[streams.size()];
+   float** sums    = new float*[streams.size()];
+   int* widths     = new int[streams.size()];
+   int* indices    = new int[streams.size()];
+   for( unsigned int ii = 0; ii < streams.size(); ii++ ){
+      BrookStreamImpl* stream  = streams[ii];
+      void* dataArrayV         = stream->getData( 1 );
+      arrays[ii]               = (float*) dataArrayV;
+      widths[ii]               =  stream->getWidth();
+      indices[ii]              = 0;
+      sums[ii]                 = new float[4];
+      sums[ii][0] = sums[ii][1] = sums[ii][2] = sums[ii][3] = 0.0f;
+      if( minIndex > stream->getStreamSize() ){
+         minIndex = stream->getStreamSize();
+      }
+   }
+
+   // sum columns 
+
+   for( int ii = 0; ii < minIndex; ii++ ){
+      for( unsigned int kk = 0; kk < streams.size(); kk++ ){
+         for( int jj = 0; jj < widths[kk]; jj++ ){
+            sums[kk][jj] += arrays[kk][indices[kk]++];
+         }
+      }
+   }
+
+   // reinitialize indices
+
+   for( unsigned int kk = 0; kk < streams.size(); kk++ ){
+      indices[kk] = 0;
+   }
+
+   // show column sums
+
+   for( unsigned int kk = 0; kk < streams.size(); kk++ ){
+      (void) fprintf( filePtr, "Sms " );
+      for( int jj = 0; jj < widths[kk]; jj++ ){
+         (void) fprintf( filePtr, "%12.5e ", sums[kk][jj] );
+      }
+   }
+   (void) fprintf( filePtr, "\n" );
+
+   for( int ii = 0; ii < minIndex; ii++ ){
+      (void) fprintf( filePtr, "%d ", ii );
+
+      // streams
+
+      for( unsigned int kk = 0; kk < streams.size(); kk++ ){
+
+//         (void) fprintf( filePtr, "[ " );
+ 
+         // ii elements of stream kk
+
+         for( int jj = 0; jj < widths[kk]; jj++ ){
+            (void) fprintf( filePtr, "%12.5e ", arrays[kk][indices[kk]++] );
+         }
+
+ //        (void) fprintf( filePtr, " ]", ii );
+      }
+      (void) fprintf( filePtr, "\n", ii );
+   }
+
+
+   // cleanup
+
+   (void) fclose( filePtr );
+
+   delete[] arrays;
+   delete[] widths;
+   delete[] indices;
+   for( int ii = 0; ii < 4; ii++ ){
+      delete[] sums[ii];
+   }
+   delete[] sums;
+
+   return DefaultReturnValue;
+
+}
+

platforms/brook/src/OpenMMBrookInterface.h

@@ -381,6 +381,27 @@ class OpenMMBrookInterface {
       
       int setParticleForces( BrookStreamImpl* forces );
 
+      /** 
+       * Print forces to file 
+       * 
+       * @param context     context
+       *
+       */
+      
+      void printForcesToFile( OpenMMContextImpl& context );
+      
+      /* 
+       * Print contents of object to file
+       *
+       * @param fileName     file name
+       * @param streams      streams to print
+       *
+       * @return DefaultReturnValue
+       *
+       * */
+      
+      int printStreamsToFile( std::string fileName, std::vector<BrookStreamImpl*>& streams );
+      
    private:
    
       static const int DefaultReturnValue = 0;

platforms/brook/src/gpu/kbonds.br

-
-/* -------------------------------------------------------------------------- *
- *                                   OpenMM                                   *
- * -------------------------------------------------------------------------- *
- * This is part of the OpenMM molecular simulation toolkit originating from   *
- * Simbios, the NIH National Center for Physics-Based Simulation of           *
- * Biological Structures at Stanford, funded under the NIH Roadmap for        *
- * Medical Research, grant U54 GM072970. See https://simtk.org.               *
- *                                                                            *
- * Portions copyright (c) 2008 Stanford University and the Authors.           *
- * Authors: Peter Eastman, Mark Friedrichs, Chris Bruns                       *
- * Contributors:                                                              *
- *                                                                            *
- * Permission is hereby granted, free of charge, to any person obtaining a    *
- * copy of this software and associated documentation files (the "Software"), *
- * to deal in the Software without restriction, including without limitation  *
- * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
- * and/or sell copies of the Software, and to permit persons to whom the      *
- * Software is furnished to do so, subject to the following conditions:       *
- *                                                                            *
- * The above copyright notice and this permission notice shall be included in *
- * all copies or substantial portions of the Software.                        *
- *                                                                            *
- * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
- * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
- * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
- * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
- * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
- * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
- * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
- * -------------------------------------------------------------------------- */
-
-//Harmonic bonds kernel
-//Input is a stream of i, j pairs
-//parms is float2( b0, kA )
-//Output is two streams of forces fi, fj
-//Can be optimized as necessary
-kernel void kbonds_harmonic(
-		float xstrwidth, //atom stream width
-		float2 atoms<>,
-		float2 parms<>,
-		float4 posq[][],
-		out float3 fi<>,
-		out float3 fj<>
-		) {
-	float2 ai, aj;
-	float3 rij;
-	float rinv;
-	
-	ai.y = floor( atoms.x / xstrwidth );
-	ai.x = atoms.x - ai.y * xstrwidth;
-
-	aj.y = floor( atoms.y / xstrwidth );
-	aj.x = atoms.y - aj.y * xstrwidth;
-
-	rij = posq[ai].xyz - posq[aj].xyz; //3
-	rinv = rsqrt( dot(rij, rij) ); //6
-	fi = -parms.y * ( 1.0f - parms.x * rinv ) * rij; //6
-	fj = -fi; 
-
-	//Total: 15 flops
-}

platforms/brook/src/gpu/kforce.h

@@ -317,7 +317,6 @@ typedef void (*gpuBondedFunction)(
 
 void  kbonded_CDLJ (const float  epsfac,
 		const float  xstrwidth,
-		const float4 params,
 		::brook::stream posq,
 		::brook::stream charge,
 		::brook::stream atoms,

platforms/brook/src/gpu/kforce_CDLJ.br

@@ -43,7 +43,7 @@ kernel float4 scalar_force_CDLJ( float4 qq, float epsfac, float4 sig, float4 eps
  return f;
 }
 
-kernel float scalar_force_single_CDLJ( float qq, float epsfac, float sig, float eps, float r2, float4 params ) {
+kernel float scalar_force_single_CDLJ( float qq, float epsfac, float sig, float eps, float r2 ) {
  float invr, invrsig2, invrsig6;
  float f;
 
@@ -572,7 +572,7 @@ kernel void kforce14_CDLJ(
 
  r2 = dot( d1, d1 );
 
- fs = scalar_force_single_CDLJ( qq, epsfac, sigeps.x, sigeps.y, r2, params );
+ fs = scalar_force_single_CDLJ( qq, epsfac, sigeps.x, sigeps.y, r2 );
 
  fi = fs.x * d1;
  fj = -fi;
@@ -732,7 +732,6 @@ kernel void kAcos( float xIn, out float acosVal<> ){
 kernel void kbonded_CDLJ (
       float epsfac,
       float xstrwidth,
-      float4 nbparams,
       float3 posq[][],
       float  charge<>,
       float4 atoms<>,
@@ -923,7 +922,7 @@ fl = float3( ddphi, 0.0f, 0.0f );
   
    if( parm4.z > -0.5f ){
       r2        = dot( ril, ril );
-      fs        = scalar_force_single_CDLJ( charge, epsfac, parm4.z, parm4.w, r2, nbparams );
+      fs        = scalar_force_single_CDLJ( charge, epsfac, parm4.z, parm4.w, r2 );
       fi_pair   = fs * ril;
       fi       += fi_pair;
       fl       -= fi_pair;