Latest

platforms/brook/src/BrookBondParameters.cpp

@@ -219,6 +219,7 @@ int BrookBondParameters::setBond( int bondIndex, int* particleIndices, double* b
 
    int numberOfParametersInBond = getNumberOfParametersInBond();
    for( int ii = 0; ii < numberOfParametersInBond; ii++ ){
+//fprintf( stderr, "Param %s %d %e\n", getBondName().c_str(), ii, bondParameters[ii] );
       _bondParameters[bondIndex].push_back( bondParameters[ii] );
    }
 

platforms/brook/src/BrookBonded.cpp

@@ -61,8 +61,11 @@ BrookBonded::BrookBonded( ){
 
    _setupCompleted            = 0;
    _numberOfParticles         = 0;
-   _ljScale                   = (BrookOpenMMFloat) 0.83333333;
+
+   //_ljScale                   = (BrookOpenMMFloat) 0.83333333;
+   _ljScale                   = 1.0;
    //_coulombFactor             = 332.0;
+
    _coulombFactor             = (BrookOpenMMFloat) 138.935485;
 
    _particleIndicesStream     = NULL;
@@ -1000,16 +1003,13 @@ int BrookBonded::addPairs( int *nbondeds, int *particles, BrookOpenMMFloat* para
 // ---------------------------------------------------------------------------------------
 
    if( debug && getLog() ){
-      (void) fprintf( getLog(), "%s npairs=%d\n", methodName.c_str(), bonded14Indices.size() );
+      (void) fprintf( getLog(), "%s npairs=%d sz=%u %u\n", methodName.c_str(), bonded14Indices.size(), bonded14Indices.size(), nonbondedParameters.size() ); fflush( getLog() );
    }
 
    for( unsigned int ii = 0; ii < bonded14Indices.size(); ii++ ){
 
-      std::vector<int> particlesIndices        = bonded14Indices[ii];
-
-      int index = 0;
-      int i     = particlesIndices[index++];
-      int j     = particlesIndices[index++];
+      int i     = bonded14Indices[ii][0];
+      int j     = bonded14Indices[ii][1];
 
       int ibonded = matchPair( i, j, *nbondeds, particles );
       if( ibonded < 0 ){
@@ -1019,13 +1019,12 @@ int BrookBonded::addPairs( int *nbondeds, int *particles, BrookOpenMMFloat* para
          (*nbondeds)++;
       }
 
-      vector<double> iParameters  = nonbondedParameters[i];
-      vector<double> jParameters  = nonbondedParameters[j];
-
-      double c6                   = 0.5*(iParameters[1] + jParameters[1]);
-      //double c12                  = lj14Scale*(iParameters[2] * jParameters[2]);
-      double c12                  = 2.0*(iParameters[2] * jParameters[2]);
+//(void) fprintf( getLog(), "%s   %d %ibonded=%d done\n", methodName.c_str(), ii, ibonded ); fflush( getLog() );
 
+      vector<double> iParameters  = nonbondedParameters[ii];
+/*
+      double c6                   = iParameters[1];
+      double c12                  = 4.0*iParameters[2];
       double sig, eps;
       if( c12 != 0.0 ){
          //eps = c6*c6/c12;
@@ -1036,17 +1035,17 @@ int BrookBonded::addPairs( int *nbondeds, int *particles, BrookOpenMMFloat* para
          eps = 0.0;
          sig = 1.0;
       }
+*/
 
-      PARAMS( ibonded, 4, 2 ) = (BrookOpenMMFloat) sig;
-      PARAMS( ibonded, 4, 3 ) = (BrookOpenMMFloat) eps;
+      PARAMS( ibonded, 4, 2 ) = (BrookOpenMMFloat) iParameters[1];
+      PARAMS( ibonded, 4, 3 ) = (BrookOpenMMFloat) (4.0*iParameters[2]);
 
       // a little wasteful, but ...
 
-      charges[i] = (BrookOpenMMFloat) iParameters[0];
-      charges[j] = (BrookOpenMMFloat) jParameters[0];
+      charges[ibonded] = (BrookOpenMMFloat) iParameters[0];
 
       if( debug ){
-         (void) fprintf( getLog(), "   %d [%d %d ] %.3e %.3e\n", ibonded, i, j, sig, eps );
+         (void) fprintf( getLog(), "   %d [%d %d ] %.3e %.3e q=%.4f\n", ibonded, i, j, iParameters[1], iParameters[2], charges[ibonded] );
       }
    }
 
@@ -1220,11 +1219,18 @@ int BrookBonded::setup( int numberOfParticles,
 // ---------------------------------------------------------------------------------------
 
    static const std::string methodName = "BrookBonded::setup";
-   static const int PrintOn            = 0;
+   static const int PrintOn            = 1;
    double dangleValue                  = 0.0;
 
 // ---------------------------------------------------------------------------------------
 
+   if( PrintOn && getLog() ){
+      (void) fprintf( getLog(), "%s particles=%d\n   [%p %p %p %p %p] (bond, angle, pd, rb, 14)\n"
+                      "14Scale=%f %f StreamW=%d StreamSz=%d\n", methodName.c_str(), numberOfParticles, harmonicBondBrookBondParameters, harmonicAngleBrookBondParameters,
+                      periodicTorsionBrookBondParameters, rbTorsionBrookBondParameters, nonBonded14ForceParameters, lj14Scale,coulombScale,
+                      particleStreamWidth, particleStreamSize ); fflush( getLog() );
+   }
+
    _numberOfParticles = numberOfParticles;
 
    // check that particle indices & parameters agree
@@ -1233,6 +1239,7 @@ int BrookBonded::setup( int numberOfParticles,
 
    int maxBonds     = 10*numberOfParticles;
    int* particles   = new int[5*maxBonds];
+   float* charges   = new BrookOpenMMFloat[maxBonds];
 
    BrookOpenMMFloat** params = new BrookOpenMMFloat*[getNumberOfParameterStreams()];
    for( int ii = 0; ii < getNumberOfParameterStreams(); ii++ ){
@@ -1250,6 +1257,7 @@ int BrookBonded::setup( int numberOfParticles,
    // All parameters must be initialized to values that will 
    // produce zero for the corresponding force. 
 
+   memset( charges, 0, maxBonds*sizeof( BrookOpenMMFloat ) );
    for( int ii = 0; ii < maxBonds; ii++ ){
 
       ATOMS( ii, 0 ) = -1;
@@ -1273,46 +1281,50 @@ int BrookBonded::setup( int numberOfParticles,
    int nbondeds = 0;
 
    if( rbTorsionBrookBondParameters ){
-      addRBTorsions( &nbondeds, particles, params, rbTorsionBrookBondParameters->getParticleIndices(),          rbTorsionBrookBondParameters->getBondParameters() );
+      addRBTorsions( &nbondeds, particles, params, rbTorsionBrookBondParameters->getParticleIndices(),         rbTorsionBrookBondParameters->getBondParameters() );
    }
    if( periodicTorsionBrookBondParameters ){
-      addPTorsions( &nbondeds, particles, params, periodicTorsionBrookBondParameters->getParticleIndices(),     periodicTorsionBrookBondParameters->getBondParameters() );
+      addPTorsions(  &nbondeds, particles, params, periodicTorsionBrookBondParameters->getParticleIndices(),   periodicTorsionBrookBondParameters->getBondParameters() );
    }
    if( harmonicAngleBrookBondParameters ){
-      addAngles(      &nbondeds, particles, params, harmonicAngleBrookBondParameters->getParticleIndices(),    harmonicAngleBrookBondParameters->getBondParameters() );
+      addAngles(     &nbondeds, particles, params, harmonicAngleBrookBondParameters->getParticleIndices(),     harmonicAngleBrookBondParameters->getBondParameters() );
    }
    if( harmonicBondBrookBondParameters ){
-      addBonds(       &nbondeds, particles, params, harmonicBondBrookBondParameters->getParticleIndices(),     harmonicBondBrookBondParameters->getBondParameters() );
+      addBonds(      &nbondeds, particles, params, harmonicBondBrookBondParameters->getParticleIndices(),      harmonicBondBrookBondParameters->getBondParameters() );
    }
-
-// ---------------------------------------------------------------------------------------
-
-   // charge stream
-
-   _chargeStream             = new BrookFloatStreamInternal( BrookCommon::BondedChargeStream, numberOfParticles, particleStreamWidth,
-                                                             BrookStreamInternal::Float, dangleValue );
-   BrookOpenMMFloat* charges = new BrookOpenMMFloat[_chargeStream->getStreamSize()];
-   memset( charges, 0, _chargeStream->getStreamSize()*sizeof( BrookOpenMMFloat ) );
-
-// ---------------------------------------------------------------------------------------
-
-//(void) fprintf( getLog(), "%s Post addBonds particles=%d number of bonds=%d maxBonds=%d\n", methodName.c_str(), numberOfParticles, nbondeds, maxBonds );
-
    if( nonBonded14ForceParameters ){
-      addPairs( &nbondeds, particles, params, charges, nonBonded14ForceParameters->getParticleIndices(), nonBonded14ForceParameters->getBondParameters(), lj14Scale, coulombScale );
+      addPairs(      &nbondeds, particles, params, charges, nonBonded14ForceParameters->getParticleIndices(), nonBonded14ForceParameters->getBondParameters(), lj14Scale, coulombScale );
    }
 
+// ---------------------------------------------------------------------------------------
+
    // check that number of bonds not too large for memory allocated
 
    if( nbondeds >= maxBonds ){
       std::stringstream message;
       message << methodName << " number of bonds=" << nbondeds << " is greater than maxBonds=" << maxBonds << " numberOfParticles=" << numberOfParticles;
       throw OpenMMException( message.str() );
+
+   } else if( nbondeds < 1 ){
+
+      // return if no bonds
+
+      (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 && getLog() ){
       (void) fprintf( getLog(), "%s particles=%d number of bonds=%d maxBonds=%d\n", methodName.c_str(), numberOfParticles, nbondeds, maxBonds );
       (void) fflush( getLog() );
    }
 
+// ---------------------------------------------------------------------------------------
+
+   // charge stream
+
+   _chargeStream             = new BrookFloatStreamInternal( BrookCommon::BondedChargeStream, nbondeds, particleStreamWidth,
+                                                             BrookStreamInternal::Float, dangleValue );
+
 // ---------------------------------------------------------------------------------------
 
    // particle indices stream
@@ -1381,13 +1393,13 @@ int BrookBonded::setup( int numberOfParticles,
                               "    phi[%10.6f %6.1f %5.1f]\n"
                               "    ang[%6.3f %8.2f] [%6.3f %8.2f]\n"
                               "      h[%8.5f %14.6e] [%8.5f %14.6e] [%8.5f %14.6e]\n"
-                              "     14[%15.6e %15.6e]\n", index,
+                              "     14[%15.6e %15.6e] q14=%15.6e\n", index,
                          ATOMS(index, 0), ATOMS(index, 1), ATOMS(index, 2), ATOMS(index, 3),
                          params[kIndex][ii], params[kIndex][ii+1], params[kIndex][ii+2], params[kIndex][ii+3], params[jIndex][ii],
                          params[jIndex][ii+1], params[jIndex][ii+2], params[jIndex][ii+3],
                          params[iIndex][ii], params[iIndex][ii+1], params[iIndex][ii+2], params[iIndex][ii+3], 
                          params[lIndex][ii], params[lIndex][ii+1], params[lIndex][ii+2], params[lIndex][ii+3],
-                         params[mIndex][ii], params[mIndex][ii+1], params[mIndex][ii+2], params[mIndex][ii+3] );
+                         params[mIndex][ii], params[mIndex][ii+1], params[mIndex][ii+2], params[mIndex][ii+3], charges[index] );
       }
    }
 

platforms/brook/src/BrookCalcNonbondedForceKernel.cpp

@@ -48,6 +48,8 @@ const std::string BrookCalcNonbondedForceKernel::BondName = "LJ14";
  * 
  * @param name                      kernel name
  * @param platform                  platform
+ * @param openMMBrookInterface      OpenMMBrookInterface reference
+ * @param system                    System reference 
  *
  */
 
@@ -68,12 +70,6 @@ BrookCalcNonbondedForceKernel::BrookCalcNonbondedForceKernel( std::string name,
 
    _log                                     = NULL;
 
-   _refForceField                           = NULL;
-   _refSystem                               = NULL;
-   _refOpenMMContext                        = NULL;
-   _referencePlatform                       = NULL;
-   _refVerletIntegrator                     = NULL;
-
    const BrookPlatform brookPlatform        = dynamic_cast<const BrookPlatform&> (platform);
    if( brookPlatform.getLog() != NULL ){
       setLog( brookPlatform.getLog() );
@@ -94,17 +90,8 @@ BrookCalcNonbondedForceKernel::~BrookCalcNonbondedForceKernel( ){
 
 // ---------------------------------------------------------------------------------------
 
-   //delete _brookBondParameters;
-
-   // deleted w/ kernel delete? If activated, program crashes
+   delete _brookBondParameters;
 
-   //delete _refForceField;
-/*
-   delete _refSystem;
-   delete _refOpenMMContext;
-   delete _referencePlatform;
-   delete _refVerletIntegrator;
-*/
 }
 
 /** 
@@ -133,9 +120,11 @@ int BrookCalcNonbondedForceKernel::setLog( FILE* log ){
 }
 
 /** 
- * Set log file reference
+ * Initialize object 
  * 
- * @param  log file reference
+ * @param  system     System reference (currently not used)
+ * @param  force      NonbondedForce reference -- extract charge, and vdw parameters from this object
+ * @param  exclusions list of execlusions
  *
  * @return  DefaultReturnValue
  *
@@ -150,6 +139,8 @@ void BrookCalcNonbondedForceKernel::initialize( const System& system, const Nonb
 // ---------------------------------------------------------------------------------------
 
     FILE* log                 = getLog();
+(void) fprintf( log, "%s begin\n", methodName.c_str() ); fflush( log );
+
     _numberOfParticles        = force.getNumParticles();
 /*
     nonbondedMethod = CalcNonbondedForceKernel::NonbondedMethod(force.getNonbondedMethod());
@@ -173,14 +164,13 @@ void BrookCalcNonbondedForceKernel::initialize( const System& system, const Nonb
    // charge & LJ parameters
 
    std::vector<std::vector<double> > nonbondedParameters;
+   nonbondedParameters.resize( _numberOfParticles );
    for( int ii = 0; ii < _numberOfParticles; ii++ ){
       double charge, radius, depth;
       force.getParticleParameters( ii, charge, radius, depth );
-      std::vector<double> particleParamArray;
-      nonbondedParameters[ii] = particleParamArray;
-      particleParamArray[0]   = radius;
-      particleParamArray[1]   = depth;
-      particleParamArray[2]   = charge;
+      nonbondedParameters[ii].push_back( charge  );
+      nonbondedParameters[ii].push_back( radius  );
+      nonbondedParameters[ii].push_back( depth   );
    }   
 
    brookNonBonded.setup( _numberOfParticles, nonbondedParameters, exclusions, getPlatform() );
@@ -219,44 +209,53 @@ void BrookCalcNonbondedForceKernel::initialize14Interactions( const System& syst
 
    FILE* log                 = getLog();
 
+//(void) fprintf( log, "%s begin\n", methodName.c_str() ); fflush( log );
+
    // ---------------------------------------------------------------------------------------
 
    // create _brookBondParameters object containing particle indices/parameters
 
    int numberOf14Forces         = force.getNumNonbonded14();
-
-   //delete _brookBondParameters;
-   _brookBondParameters = new BrookBondParameters( BondName, NumberOfParticlesInBond, NumberOfParametersInBond, numberOf14Forces, getLog() );
-
-   for( int ii = 0; ii < numberOf14Forces; ii++ ){
-
-      int particle1, particle2;
-      double  charge, radius, depth;
-
-      int particles[NumberOfParticlesInBond];
-      double parameters[NumberOfParametersInBond];
-
-      force.getNonbonded14Parameters( ii, particle1, particle2,  charge, radius, depth ); 
-      particles[0]    = particle1;
-      particles[1]    = particle2;
- 
-      parameters[0]   = charge;
-      parameters[1]   = radius;
-      parameters[2]   = depth;
-
-      _brookBondParameters->setBond( ii, particles, parameters );
-   }   
-   _openMMBrookInterface.setNonBonded14ForceParameters( _brookBondParameters, 1.0, 1.0 );
-
-   if( log ){
-      std::string contents = _brookBondParameters->getContentsString( ); 
-      (void) fprintf( log, "%s contents:\n%s", methodName.c_str(), contents.c_str() );
+   if( numberOf14Forces > 0 ){
+
+      _brookBondParameters         = new BrookBondParameters( BondName, NumberOfParticlesInBond, NumberOfParametersInBond, numberOf14Forces, getLog() );
+   
+      for( int ii = 0; ii < numberOf14Forces; ii++ ){
+   
+         int particle1, particle2;
+         double  charge, radius, depth;
+   
+         int particles[NumberOfParticlesInBond];
+         double parameters[NumberOfParametersInBond];
+   
+         force.getNonbonded14Parameters( ii, particle1, particle2, charge, radius, depth ); 
+   
+(void) fprintf( log, "%s idx=%d [%d %d] [%f %f %f]\n", methodName.c_str(), ii, particle1, particle2, charge, radius, depth );
+         particles[0]    = particle1;
+         particles[1]    = particle2;
+    
+         parameters[0]   = charge;
+         parameters[1]   = radius;
+         parameters[2]   = depth;
+   
+         _brookBondParameters->setBond( ii, particles, parameters );
+      }   
+      _openMMBrookInterface.setNonBonded14ForceParameters( _brookBondParameters, 1.0, 1.0 );
+   
+      if( log ){
+         std::string contents = _brookBondParameters->getContentsString( ); 
+         (void) fprintf( log, "%s contents:\n%s", methodName.c_str(), contents.c_str() );
+         (void) fflush( log );
+      }
+   } else if( log ){
+      (void) fprintf( log, "%s no 14 ixns\n", methodName.c_str() );
       (void) fflush( log );
    }
 
    // ---------------------------------------------------------------------------------------
     
 }
+
 /** 
  * Execute the kernel to calculate the nonbonded forces
  * 
@@ -304,35 +303,3 @@ double BrookCalcNonbondedForceKernel::executeEnergy( OpenMMContextImpl& context
    }   
 
 }
-
-/**
- * Get reference Context
- * 
- * @param numberOfParticles  number of particles
- *
- * @return  OpenMMContext
- *
- */
-
-OpenMMContext* BrookCalcNonbondedForceKernel::getReferenceOpenMMContext( int numberOfParticles ){
-
-// ---------------------------------------------------------------------------------------
-
-   //static const std::string methodName      = "BrookCalcNonbondedForceKernel::getReferenceOpenMMContext";
-
-// ---------------------------------------------------------------------------------------
-
-   if( _refOpenMMContext == NULL ){
-      
-      _referencePlatform       = new ReferencePlatform();
-      _refSystem               = new System( numberOfParticles, 0 ); 
-		_refVerletIntegrator     = new VerletIntegrator( 0.01 );
-
-      _refSystem->addForce( _refForceField );
-
-      _refOpenMMContext        = new OpenMMContext( *_refSystem, *_refVerletIntegrator, *_referencePlatform );
-   }
-
-   return _refOpenMMContext;
-
-}

platforms/brook/src/BrookCalcNonbondedForceKernel.h

@@ -33,7 +33,6 @@
  * -------------------------------------------------------------------------- */
 
 #include "kernels.h"
-//#include "../../reference/src/SimTKUtilities/SimTKOpenMMRealType.h"
 #include "OpenMMBrookInterface.h"
 #include "NonbondedForce.h"
 
@@ -110,17 +109,6 @@ class BrookCalcNonbondedForceKernel : public CalcNonbondedForceKernel {
   
       double executeEnergy( OpenMMContextImpl& context );
 
-      /**
-       * Get reference Context
-       * 
-       * @param numberOfParticles  number of particles
-       *
-       * @return  OpenMMContext
-       *
-       */
-      
-      OpenMMContext* getReferenceOpenMMContext( int numberOfParticles );
-      
       /** 
        * Set log file reference
        * 
@@ -158,6 +146,8 @@ class BrookCalcNonbondedForceKernel : public CalcNonbondedForceKernel {
 
       static const std::string BondName;
 
+      // number of LJ14 particles/parameters in 'bond'
+
       static const int NumberOfParticlesInBond  = 2;
       static const int NumberOfParametersInBond = 3;
 
@@ -174,14 +164,6 @@ class BrookCalcNonbondedForceKernel : public CalcNonbondedForceKernel {
 
        BrookBondParameters* _brookBondParameters;
 
-       // used to calculate energy
-
-       NonbondedForce* _refForceField;
-       System*               _refSystem;
-       OpenMMContext*        _refOpenMMContext;
-       ReferencePlatform*    _referencePlatform;
-       VerletIntegrator*     _refVerletIntegrator;
-
 };
 
 } // namespace OpenMM

platforms/brook/src/BrookNonBonded.cpp

@@ -45,7 +45,7 @@ using namespace std;
  * 
  */
 
-BrookNonBonded::BrookNonBonded(  ){
+BrookNonBonded::BrookNonBonded( ){
 
 // ---------------------------------------------------------------------------------------
 
@@ -958,8 +958,7 @@ int BrookNonBonded::setup( int numberOfParticles, const std::vector<std::vector<
 
 // ---------------------------------------------------------------------------------------
 
-   static const std::string methodName      = "BrookNonBonded::setup";
-   static const int debug                   = 1;
+   //static const std::string methodName      = "BrookNonBonded::setup";
 
 // ---------------------------------------------------------------------------------------
 
@@ -972,6 +971,8 @@ int BrookNonBonded::setup( int numberOfParticles, const std::vector<std::vector<
    _initializeVdwAndCharge( nonbondedParameters, platform );
    _initializeJStreamVdw( nonbondedParameters, platform );
 
+   setIsActive( 1 );
+
    return DefaultReturnValue;
 }
 
@@ -1325,35 +1326,38 @@ nonbondedForceStreams[3]->fillWithValue( &zerof );
 
    // diagnostics
 
-   if( 1 && PrintOn ){
-      (void) fprintf( getLog(), "\nPost knbforce_CDLJ4: particles=%6d ceiling=%3d dupFac=%3d", getNumberOfParticles(),  
+   if( 1 && PrintOn && getLog() ){
+      FILE* log = getLog();
+      (void) fprintf( log, "%s\n", methodName.c_str() ); (void) fflush( log );
+      (void) fprintf( log, "\nPost knbforce_CDLJ4: particles=%6d ceiling=%3d dupFac=%3d", getNumberOfParticles(),  
                                                                                                getParticleSizeCeiling(),
                                                                                                getDuplicationFactor()  );
 
-      (void) fprintf( getLog(), "\n                      hght=%6d   width=%3d   jWid=%3d", getParticleStreamHeight( ),
+      (void) fprintf( log, "\n                      hght=%6d   width=%3d   jWid=%3d", getParticleStreamHeight( ),
                                                                                            getParticleStreamWidth( ),
                                                                                            getJStreamWidth( ) );
-      (void) fprintf( getLog(), "\n                      pFrc=%6d     eps=%12.5e\n",       getPartialForceStreamWidth( ), epsfac );
+      (void) fprintf( log, "\n                      pFrc=%6d     eps=%12.5e\n",       getPartialForceStreamWidth( ), epsfac );
 
-      (void) fprintf( getLog(), "\nOuterVdwStreamd\n" );
-      getOuterVdwStream()->printToFile( getLog() );
+      (void) fprintf( log, "\nOuterVdwStreamd\n" );
+      getOuterVdwStream()->printToFile( log );
 
-      (void) fprintf( getLog(), "\nInnerSigmaStream\n" );
-      getInnerSigmaStream()->printToFile( getLog() );
+      (void) fprintf( log, "\nInnerSigmaStream\n" );
+      getInnerSigmaStream()->printToFile( log );
 
-      (void) fprintf( getLog(), "\nInnerEpsilonStream\n" );
-      getInnerEpsilonStream()->printToFile( getLog() );
+      (void) fprintf( log, "\nInnerEpsilonStream\n" );
+      getInnerEpsilonStream()->printToFile( log );
 
-      (void) fprintf( getLog(), "\nExclusionStream\n" );
-      getExclusionStream()->printToFile( getLog() );
+      (void) fprintf( log, "\nExclusionStream\n" );
+      getExclusionStream()->printToFile( log );
 
-      (void) fprintf( getLog(), "\nChargeStream\n" );
-      getChargeStream()->printToFile( getLog() );
+      (void) fprintf( log, "\nChargeStream\n" );
+      getChargeStream()->printToFile( log );
 
       for( int ii = 0; ii < 4; ii++ ){
-         (void) fprintf( getLog(), "\nForce stream %d\n", ii );
-         nonbondedForceStreams[ii]->printToFile( getLog() );
+         (void) fprintf( log, "\nForce stream %d\n", ii );
+         nonbondedForceStreams[ii]->printToFile( log );
       }
+      (void) fflush( log );
    }
 
 // ---------------------------------------------------------------------------------------
@@ -1374,11 +1378,12 @@ nonbondedForceStreams[3]->fillWithValue( &zerof );
 
    // diagnostics
 
-   if( 0 && PrintOn ){
+   if( PrintOn && getLog() ){
 
-      (void) fprintf( getLog(), "\nNB forces" );
-      BrookStreamInternal* brookStreamInternalF   = forceStream.getBrookStreamImpl();
-      brookStreamInternalF->printToFile( getLog() );
+      FILE* log = getLog();
+      (void) fprintf( log, "\n%s NB forces", methodName.c_str() ); (void) fflush( log );
+//      forceStream.printToFile( log );
+      (void) fprintf( log, "\n%s Done", methodName.c_str() ); (void) fflush( log );
 
    }
 

platforms/brook/src/BrookPlatform.cpp

@@ -197,7 +197,6 @@ BrookPlatform::BrookPlatform( ){
 
    _particleStreamWidth  = DefaultParticleStreamWidth;
    _log                  = NULL;
-   _openMMBrookInterface = NULL;
 
    // get Brook runtime
 
@@ -234,8 +233,6 @@ BrookPlatform::BrookPlatform( int particleStreamWidth, const std::string& runtim
 
 // ---------------------------------------------------------------------------------------
 
-   _openMMBrookInterface = NULL;
-
    _initializeKernelFactory( );
    _setBrookRuntime( runtime );
 
@@ -476,6 +473,4 @@ void BrookPlatform::contextDestroyed( OpenMMContextImpl& context ) const {
 
 // ---------------------------------------------------------------------------------------
 
-   delete _openMMBrookInterface;
-
 }

platforms/brook/src/OpenMMBrookInterface.cpp

@@ -581,19 +581,22 @@ void OpenMMBrookInterface::computeForces( OpenMMContextImpl& context ){
 
    // bonded forces
 
+(void) fprintf( stderr, "%s done nonbonded: bonded=%d completed=%d\n", methodName.c_str(), _brookBonded.isActive(), _brookBonded.isSetupCompleted() ); (void) fflush( stderr );
    if( _brookBonded.isActive() ){
 
 //(void) fprintf( stderr, "%s Bonded\n", methodName.c_str() ); (void) fflush( stderr );
       // perform setup first time through
 
-      if( _brookBonded.isSetupCompleted() == 0 ){
+      if( _brookBonded.isSetupCompleted() >= -1 ){
          _brookBonded.setup( getNumberOfParticles(), getHarmonicBondForceParameters(), getHarmonicAngleForceParameters(),
-                              getPeriodicTorsionForceParameters(), getRBTorsionForceParameters(),
-                              getNonBonded14ForceParameters(),
-                              getLj14Scale(), getCoulomb14Scale(), getParticleStreamWidth(), getParticleStreamSize() );
+                             getPeriodicTorsionForceParameters(), getRBTorsionForceParameters(),
+                             getNonBonded14ForceParameters(),
+                             getLj14Scale(), getCoulomb14Scale(), getParticleStreamWidth(), getParticleStreamSize() );
       }
 
+(void) fprintf( stderr, "%s done setup bonded=%d completed=%d\n", methodName.c_str(), _brookBonded.isActive(), _brookBonded.isSetupCompleted() ); (void) fflush( stderr );
       _brookBonded.computeForces( *positions, *forces );
+(void) fprintf( stderr, "%s done forces bonded=%d completed=%d\n", methodName.c_str(), _brookBonded.isActive(), _brookBonded.isSetupCompleted() ); (void) fflush( stderr );
    
       // diagnostics
    

platforms/brook/src/gpu/kforce_CDLJ.br

@@ -583,7 +583,7 @@ kernel void kbonded_CDLJ (
   float xstrwidth,
   float4 nbparams,
   float3 posq[][],
-  float  charge[][],
+  float  charge<>,
   float4 atoms<>,
   float4 parm0<>,
   float4 parm1<>,
@@ -596,6 +596,7 @@ kernel void kbonded_CDLJ (
   out float3 fl<>
   )
 {
+
  float3 rij, rkj, rkl, ril;
  float2 ai, aj, ak, al;
  float3 m, n;
@@ -619,23 +620,11 @@ kernel void kbonded_CDLJ (
  float theta1, costheta1, invsintheta1;
  float theta2, costheta2, invsintheta2;
  float3 posqi, posqj, posqk, posql;
- float  chargei, chargel;
  float sig, eps;
 
-
-
-
-
-
-   ai.y = round( (atoms.x - fmod( atoms.x, xstrwidth ))/xstrwidth );
+ ai.y = round( (atoms.x - fmod( atoms.x, xstrwidth ))/xstrwidth );
  ai.x = atoms.x - ai.y * xstrwidth;
  posqi = posq[ ai ];
- chargei = charge[ ai ];
-
-
-
-
-
 
  if ( atoms.y > -0.5f ) {
 
@@ -656,17 +645,13 @@ kernel void kbonded_CDLJ (
   posqk = float4( 0.0f, 1.0f, 1.0f, 0.0f );
 
  if ( atoms.w > -0.5f ) {
-
-      al.y = round( (atoms.w - fmod( atoms.w, xstrwidth ))/xstrwidth );
+  al.y = round( (atoms.w - fmod( atoms.w, xstrwidth ))/xstrwidth );
   al.x = atoms.w - al.y * xstrwidth;
   posql = posq[ al ];
-  chargel = charge[ al ];
  } else {
   posql = float4( 1.0f, 1.0f, 1.0f, 0.0f );
-  chargel = 0.0f;
-}
+ }
 
- qq = chargei*chargel;
  rij = posqi.xyz - posqj.xyz;
  rkj = posqk.xyz - posqj.xyz;
  rkl = posqk.xyz - posql.xyz;
@@ -789,7 +774,7 @@ kernel void kbonded_CDLJ (
 
  if ( parm4.z > -0.5f ) {
   r2 = dot( ril, ril );
-  fs = scalar_force_single_CDLJ( qq, epsfac, parm4.z, parm4.w, r2, nbparams );
+  fs = scalar_force_single_CDLJ( charge, epsfac, parm4.z, parm4.w, r2, nbparams );
   fi_pair = fs * ril;
   fi += fi_pair;
   fl -= fi_pair;

platforms/brook/src/gpu/knlforce.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.                                     *
- * -------------------------------------------------------------------------- */
-
-/* 
- * 3rd attempt at writing a force kernel that uses neighborlists.
- * Keeping things very simple this time - no unrolling at all.
- *
- * Starting with simple Coul-LJ
- * */
-
-//Constant dielectric, LJ (normal forces)
-kernel float4 nl_scalar_force_CDLJ( float4 qq, float epsfac, float4 sig, float4 eps, float4 r2 ) {
-	float4 invr, invrsig2, invrsig6;
-	float4 f;
-
-	invr = rsqrt( r2 ); // 1 or 2 flops?
-	invrsig2 = invr * sig; //1 flop
-	invrsig2 = invrsig2 * invrsig2; //1 flop
-	invrsig6 = invrsig2 * invrsig2 * invrsig2; //2flops
-	f = eps * ( 12.0f*invrsig6 - 6.0f ) * invrsig6; //4flops
-	f += epsfac*qq*invr; //2 flops
-	f *= invr*invr; //2 flops 
-	return f; //Total ? flops
-}
-
-kernel float4 nl_get_r2( float3 d1, float3 d2, float3 d3, float3 d4 ) {
-	float4 r2;
-	r2 = float4( dot(d1, d1), dot( d2, d2 ), dot( d3, d3 ), dot( d4, d4 ) ); //5*4 flops
-	return r2;
-}
-
-kernel void knbforce_nl(
-		float AtomStrWidth,  //Width of atom position stream
-		float epsfac,
-		iter float2 wpos<>,  //pixel position of output (i-atom)
-		float4 posq[][],     //coordinates and charges
-		float3 inforce<>,
-		float4 nlist0<>,     //nbor indices
-		float4 sig0<>,        //LJ sigmas
-		float4 eps0<>,        //LJ epsilons
-		//We'll add more nlists here to reduce number of kernel calls
-		out float3 force<>   //output forces
-		)
-{
-	float2 jind;
-	float4 ind_tmp1, ind_tmp2;
-	float4 qj;
-	float3 jpos1, jpos2, jpos3, jpos4;
-	float4 r2, fs, qq;
-	float3 d1, d2, d3, d4;
-	float3 ipos;
-	float qi;
-	
-	//Since neighbors will not be adjacent in memory
-	//We may want to stagger compute and etch, but 
-	//with fxc, we don't have control over scheduling anyway.
-	//Will try hand scheduling the ps3 code later
-
-	//etch i parameters
-	ipos = posq[ wpos ].xyz;
-	qi   = posq[ wpos ].w;
-
-	//etch j parameters
-	
-	ind_tmp1 = floor( nlist0 / AtomStrWidth );
-	ind_tmp2 = nlist0 - ind_tmp1 * AtomStrWidth;
-
-	force = inforce;
-	
-	//1st set of 4 neighbors
-	jind.y = ind_tmp1.x;
-	jind.x = ind_tmp2.x;
-	jpos1  = posq[ jind ].xyz;
-	qj.x   = posq[ jind ].w;
-	
-	jind.y = ind_tmp1.y;
-	jind.x = ind_tmp2.y;
-	jpos2  = posq[ jind ].xyz;
-	qj.y   = posq[ jind ].w;
-	
-	jind.y = ind_tmp1.z;
-	jind.x = ind_tmp2.z;
-	jpos3  = posq[ jind ].xyz;
-	qj.z   = posq[ jind ].w;
-
-	
-	jind.y = ind_tmp1.w;
-	jind.x = ind_tmp2.w;
-	jpos4  = posq[ jind ].xyz;
-	qj.w   = posq[ jind ].w;
-
-
-	d1 = ipos - jpos1;
-	d2 = ipos - jpos2;
-	d3 = ipos - jpos3;
-	d4 = ipos - jpos4;
-	
-	r2 = nl_get_r2( d1, d2, d3, d4 );
-
-	qq = qi * qj;
-
-	fs = nl_scalar_force_CDLJ( qq, epsfac, sig0, eps0, r2 );
-	
-	force += fs.x * d1;
-	force += fs.y * d2;
-	force += fs.z * d3;
-	force += fs.w * d4;	
-}
-
-

platforms/brook/src/gpu/pre_test.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.                                     *
- * -------------------------------------------------------------------------- */
-
-/* pre-tests, known to be right */
-
-kernel void pre_test0( float input<>, out float output<> )
-{
-  output = input;
-}
-
-kernel void pre_test1( float input<>, out float output<> )
-{
-  output = ( input * input ) + input;
-}