Initial unit tests for StandardMMForceField pass

platforms/brook/src/BrookBonded.cpp

@@ -59,8 +59,9 @@ BrookBonded::BrookBonded( ){
 // ---------------------------------------------------------------------------------------
 
    _numberOfAtoms             = 0;
-   _ljScale                   = 1.0;
-   _coulombFactor             = 332.0;
+   _ljScale                   = (BrookOpenMMFloat) 0.83333333;
+   //_coulombFactor             = 332.0;
+   _coulombFactor             = (BrookOpenMMFloat) 138.935485;
 
    _atomIndicesStream         = NULL;
    _chargeStream              = NULL;
@@ -354,6 +355,24 @@ BrookFloatStreamInternal* BrookBonded::getAtomIndicesStream( void ) const {
    return _atomIndicesStream;
 }
 
+/** 
+ * Get bonded charge stream
+ * 
+ * @return  charge stream
+ *
+ */
+
+BrookFloatStreamInternal* BrookBonded::getChargeStream( void ) const {
+
+// ---------------------------------------------------------------------------------------
+
+   // static const std::string methodName = "BrookBonded::getChargeStream";
+
+// ---------------------------------------------------------------------------------------
+
+   return _chargeStream;
+}
+
 /** 
  * Get array of bonded parameter streams
  * 
@@ -593,8 +612,7 @@ int BrookBonded::matchBond( int i, int j, int nbondeds, int *atoms, int *flag ){
             *flag = 1;
             return n;
          }
-      }
-      else if(  ( i == ATOMS(n, 1) && j == ATOMS(n, 2) ) ||( i == ATOMS(n, 2) && j == ATOMS(n, 1) ) ){
+      } else if(  ( i == ATOMS(n, 1) && j == ATOMS(n, 2) ) ||( i == ATOMS(n, 2) && j == ATOMS(n, 1) ) ){
          *flag = 1;
          return n;
       }
@@ -880,7 +898,7 @@ int BrookBonded::addBonds( int *nbondeds, int *atoms, float *params[], const vec
 // ---------------------------------------------------------------------------------------
 
    static const std::string methodName = "BrookBonded::addBonds";
-   static const int debug              = 0;
+   static const int debug              = 1;
 
 // ---------------------------------------------------------------------------------------
 
@@ -899,8 +917,17 @@ int BrookBonded::addBonds( int *nbondeds, int *atoms, float *params[], const vec
       int i     = atomsIndices[index++];
       int j     = atomsIndices[index++];
 
+      // insure i < j
+
+      if( i > j ){
+         int k = i;
+         i     = j;
+         j     = k;
+      }
+
       int flag;
       int ibonded = matchBond( i, j, *nbondeds, atoms, &flag );
+int saveIbond = ibonded;
       if( ibonded < 0 ){
          ibonded = *nbondeds;
          ATOMS( ibonded, 0 ) = i;
@@ -918,7 +945,7 @@ int BrookBonded::addBonds( int *nbondeds, int *atoms, float *params[], const vec
       }
 
       if( debug && getLog() ){
-         (void) fprintf( getLog(), "   %d [%d %d ] flag=%d %.3e %.3e\n", ibonded, i, j, flag,
+         (void) fprintf( getLog(), "   %d (%5d) [%6d %6d ] flag=%2d %.3e %.3e\n", ibonded, saveIbond, i, j, flag,
                          bndParameters[0], bndParameters[1] );
       }
    }
@@ -951,9 +978,9 @@ int BrookBonded::addPairs( int *nbondeds, int *atoms, BrookOpenMMFloat* params[]
 
 // ---------------------------------------------------------------------------------------
 
-   static const std::string methodName = "BrookBonded::BrookBonded";
+   static const std::string methodName = "BrookBonded::addPairs";
    static const double oneSixth        = 1.0/6.0;
-   static const int debug              = 0;
+   static const int debug              = 1;
 
 // ---------------------------------------------------------------------------------------
 
@@ -980,13 +1007,16 @@ int BrookBonded::addPairs( int *nbondeds, int *atoms, BrookOpenMMFloat* params[]
       vector<double> iParameters  = nonbondedParameters[i];
       vector<double> jParameters  = nonbondedParameters[j];
 
-      double c6                   = iParameters[0] + jParameters[0];
-      double c12                  = lj14Scale*(iParameters[1] * jParameters[1]);
+      double c6                   = 0.5*(iParameters[1] + jParameters[1]);
+      //double c12                  = lj14Scale*(iParameters[2] * jParameters[2]);
+      double c12                  = 2.0*(iParameters[2] * jParameters[2]);
 
       double sig, eps;
       if( c12 != 0.0 ){
-         eps = c6*c6/c12;
-         sig = pow( c12/c6, oneSixth );
+         //eps = c6*c6/c12;
+         //sig = pow( c12/c6, oneSixth );
+         eps = c12;
+         sig = c6;
       } else {
          eps = 0.0;
          sig = 1.0;
@@ -997,11 +1027,11 @@ int BrookBonded::addPairs( int *nbondeds, int *atoms, BrookOpenMMFloat* params[]
 
       // a little wasteful, but ...
 
-      charges[i] = (BrookOpenMMFloat) iParameters[2];
-      charges[j] = (BrookOpenMMFloat) jParameters[2];
+      charges[i] = (BrookOpenMMFloat) iParameters[0];
+      charges[j] = (BrookOpenMMFloat) jParameters[0];
 
       if( debug ){
-         (void) fprintf( getLog(), "\n %d [%d %d ] %.3e %.3e", ibonded, i, j, sig, eps );
+         (void) fprintf( getLog(), "   %d [%d %d ] %.3e %.3e\n", ibonded, i, j, sig, eps );
       }
    }
 
@@ -1037,22 +1067,32 @@ int BrookBonded::loadInvMaps( int nbondeds, int natoms, int *atoms, const BrookP
    int atomStreamSize                           = brookPlatform.getStreamSize( getNumberOfAtoms(), atomStreamWidth, NULL );
    _inverseMapStreamWidth                       = atomStreamWidth;
    
+// ---------------------------------------------------------------------------------------
+
    // allocate temp memory
 
    float4** invmaps = new float4*[getMaxInverseMapStreamCount()];
    float* block     = new float[4*getMaxInverseMapStreamCount()*atomStreamSize];
+
+   //memset( block, 0, 4*getMaxInverseMapStreamCount()*atomStreamSize*sizeof( float ) );
+
+   float* blockPtr = block;
    for( int ii = 0; ii < getMaxInverseMapStreamCount(); ii++ ){
-      invmaps[ii]  = (float4*) block;
-      block       += 4*atomStreamSize;
+      invmaps[ii]  = (float4*) blockPtr;
+      blockPtr    += 4*atomStreamSize;
    }
    int* counts = new int[atomStreamSize];
 
+// ---------------------------------------------------------------------------------------
+
    // get inverse maps and load into streams
 
    // create streams
    // done independently from loading since for test cases some stream counts == 0, but kernels expect stream to
    // have been created even though unused
 
+   // load streams -- initialize all streams even if unused since gather methods will still pick up
+
    for( int ii = 0; ii < getNumberOfForceStreams(); ii++ ){
       for( int jj = 0; jj < getMaxInverseMapStreamCount(ii); jj++ ){
          _inverseStreamMaps[ii][jj] = new BrookFloatStreamInternal( BrookCommon::BondedInverseMapStreams, atomStreamSize,
@@ -1060,28 +1100,60 @@ int BrookBonded::loadInvMaps( int nbondeds, int natoms, int *atoms, const BrookP
       }
    }
 
+   if( getLog() ){
+      (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),
+                      atomStreamSize, atomStreamWidth );
+      (void) fflush( getLog() );
+   }
+
    // load data
 
    for( int ii = 0; ii < getNumberOfForceStreams(); ii++ ){
       gpuCalcInvMap( ii, 4, nbondeds, natoms, atoms, getInverseMapStreamCount( ii ), counts, invmaps, &(_inverseMapStreamCount[ii]) );
+//gpuPrintInvMaps( _inverseMapStreamCount[ii], natoms, counts, invmaps, getLog() );
       validateInverseMapStreamCount( ii, _inverseMapStreamCount[ii] ); 
       for( int jj = 0; jj < _inverseMapStreamCount[ii]; jj++ ){
          _inverseStreamMaps[ii][jj]->loadFromArray( invmaps[jj] );
+
+         (void) fprintf( getLog(), "%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 < atomStreamSize; 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  );
+         }
+*/
+
+      }    
+
+      // for small systems, must all initialize inverse maps to negative values in order to 
+      // keep invalid entries from being included in forces
+
+      if( _inverseMapStreamCount[ii] < getMaxInverseMapStreamCount( ii ) ){
+         for( int jj = 0; jj < 4*atomStreamSize; jj++ ){
+            block[jj] = -1.0f;
+         }
+         for( int jj = _inverseMapStreamCount[ii]; jj < getMaxInverseMapStreamCount( ii ); jj++ ){
+             _inverseStreamMaps[ii][jj]->loadFromArray( invmaps[0] );
+         }
       }    
 
    }
 
+   // diagnostics 
+
    if( getLog() ){
       (void) fprintf( getLog(), "%s done\n", methodName.c_str() );
       (void) fflush( getLog() );
-      //gpuPrintInvMaps( bp->nimaps, natoms, counts, invmaps, gpu->log );
    }
 
    // free memory
 
-   delete counts;
-   delete invmaps[0];
-   delete invmaps;
+   delete[] counts;
+   delete[] invmaps[0];
+   delete[] invmaps;
 
    return DefaultReturnValue;
 }
@@ -1191,7 +1263,13 @@ int BrookBonded::setup( int numberOfAtoms,
    for( int ii = 0; ii < getNumberOfParameterStreams(); ii++ ){
       params[ii] = new BrookOpenMMFloat[4*maxBonds];
    }
-   BrookOpenMMFloat* charges = new BrookOpenMMFloat[numberOfAtoms];
+
+// ---------------------------------------------------------------------------------------
+
+   // build streams
+
+   const BrookStreamFactory& brookStreamFactory = dynamic_cast<const BrookStreamFactory&> (brookPlatform.getDefaultStreamFactory() );
+   int atomStreamWidth                          = brookStreamFactory.getDefaultAtomStreamWidth();
 
    // Initialize all atom indices to -1 to indicate empty slots
    // All parameters must be initialized to values that will 
@@ -1222,6 +1300,19 @@ int BrookBonded::setup( int numberOfAtoms,
    addAngles(      &nbondeds, atoms, params, angleIndices,           angleParameters           );
    addBonds(       &nbondeds, atoms, params, bondIndices,            bondParameters            );
 
+// ---------------------------------------------------------------------------------------
+
+   // charge stream
+
+   _chargeStream         = new BrookFloatStreamInternal( BrookCommon::BondedChargeStream, numberOfAtoms, atomStreamWidth,
+                                                         BrookStreamInternal::Float, dangleValue );
+   BrookOpenMMFloat* charges = new BrookOpenMMFloat[_chargeStream->getStreamSize()];
+   memset( charges, 0, _chargeStream->getStreamSize()*sizeof( BrookOpenMMFloat ) );
+
+// ---------------------------------------------------------------------------------------
+
+//(void) fprintf( getLog(), "%s Post addBonds atoms=%d number of bonds=%d maxBonds=%d\n", methodName.c_str(), numberOfAtoms, nbondeds, maxBonds );
+
    addPairs( &nbondeds, atoms, params, charges, bonded14Indices, nonbondedParameters, lj14Scale, coulombScale );
 
    // check that number of bonds not too large for memory allocated
@@ -1235,25 +1326,42 @@ int BrookBonded::setup( int numberOfAtoms,
       (void) fflush( getLog() );
    }
 
-   const BrookStreamFactory& brookStreamFactory = dynamic_cast<const BrookStreamFactory&> (brookPlatform.getDefaultStreamFactory() );
-   int atomStreamWidth                          = brookStreamFactory.getDefaultAtomStreamWidth();
+// ---------------------------------------------------------------------------------------
 
-   // build streams
+   // atom indices stream
 
    _atomIndicesStream    = new BrookFloatStreamInternal( BrookCommon::BondedAtomIndicesStream, nbondeds, atomStreamWidth,
                                                          BrookStreamInternal::Float4, dangleValue );
-   _atomIndicesStream->loadFromArray( atoms, BrookStreamInternal::Integer ); 
 
-   _chargeStream         = new BrookFloatStreamInternal( BrookCommon::BondedChargeStream, numberOfAtoms, atomStreamWidth,
-                                                         BrookStreamInternal::Float, dangleValue );
+   int* buffer           = new int[4*_atomIndicesStream->getStreamSize()];
+   memset( buffer, 0, sizeof( int )*4*_atomIndicesStream->getStreamSize() );
+
+   int index             = 0;
+   for( int ii = 0; ii < nbondeds; ii++ ){
+      for( int jj = 0; jj < 4; jj++ ){
+         buffer[index++] = ATOMS( ii, jj );
+//(void) fprintf( getLog(), "%s atomIndices %d %d  %d buffer=%d atoms=%d\n", methodName.c_str(), ii, jj, index, buffer[index-1], ATOMS( ii, jj ) );
+      }
+   }
+   _atomIndicesStream->loadFromArray( buffer, BrookStreamInternal::Integer ); 
+   delete[] buffer;
+
+// ---------------------------------------------------------------------------------------
+
    _chargeStream->loadFromArray( charges ); 
 
+// ---------------------------------------------------------------------------------------
+
+   // bonded parameters
+
    for( int ii = 0; ii < getNumberOfParameterStreams(); ii++ ){
       _bondedParameters[ii]  = new BrookFloatStreamInternal( BrookCommon::BondedParametersStream, nbondeds, atomStreamWidth,
                                                              BrookStreamInternal::Float4, dangleValue );
       _bondedParameters[ii]->loadFromArray( params[ii] );
    }
 
+// ---------------------------------------------------------------------------------------
+
    // debug stuff
 
    if( 1 && getLog() ){
@@ -1300,19 +1408,22 @@ int BrookBonded::setup( int numberOfAtoms,
 
    loadInvMaps( nbondeds, getNumberOfAtoms(), atoms, brookPlatform );
    
+// ---------------------------------------------------------------------------------------
+
    // free memory
 
    for( int ii = 0; ii < getNumberOfParameterStreams(); ii++ ){
       delete[] params[ii];
    }
+
    delete[] params;
    delete[] atoms;
    delete[] charges;
 
    // set the fudge factors
 
-   _ljScale        = (BrookOpenMMFloat) lj14Scale;
-   _coulombFactor  = (BrookOpenMMFloat) coulombScale;
+   //_ljScale        = (BrookOpenMMFloat) lj14Scale;
+   //_coulombFactor  = (BrookOpenMMFloat) coulombScale;
 
    // initialize output streams
 

platforms/brook/src/BrookBonded.h

@@ -186,6 +186,15 @@ class BrookBonded : public BrookCommon {
       
       BrookFloatStreamInternal* getAtomIndicesStream( void ) const;
       
+      /** 
+       * Get bonded charge stream
+       * 
+       * @return  charge stream
+       *
+       */
+      
+      BrookFloatStreamInternal* getChargeStream( void ) const;
+      
       /** 
        * Get array of bonded parameter streams
        * 

platforms/brook/src/BrookCalcStandardMMForceFieldKernel.cpp

@@ -223,19 +223,109 @@ void BrookCalcStandardMMForceFieldKernel::executeForces( const Stream& positions
    static const int K_Stream                = 2;
    static const int L_Stream                = 3;
 
+   static const int PrintOn                 = 0;
+
    static const float4 dummyParameters( 0.0, 0.0, 0.0, 0.0 );
 
    // static const int debug                   = 1;
 
 // ---------------------------------------------------------------------------------------
 
-   // bonded
+   // ---------------------------------------------------------------------------------------
 
    const BrookStreamImpl& positionStreamC              = dynamic_cast<const BrookStreamImpl&> (positions.getImpl());
    BrookStreamImpl& positionStream                     = const_cast<BrookStreamImpl&>         (positionStreamC);
    BrookStreamImpl& forceStream                        = dynamic_cast<BrookStreamImpl&>       (forces.getImpl());
    
-   float epsfac                                        = (float) (_brookBonded->getLJ_14Scale()*_brookBonded->getCoulombFactor());
+   // nonbonded forces
+
+   // added charge stream to knbforce_CDLJ4
+   // libs generated from ~/src/gmxgpu-nsqOpenMM
+
+   BrookFloatStreamInternal**  nonbondedForceStreams = _brookNonBonded->getForceStreams();
+
+   float epsfac                                      = 138.935485f;
+
+   knbforce_CDLJ4(
+             (float) _brookNonBonded->getNumberOfAtoms(),
+             (float) _brookNonBonded->getAtomSizeCeiling(),
+             (float) _brookNonBonded->getDuplicationFactor(),
+             (float) _brookNonBonded->getAtomStreamHeight( ),
+             (float) _brookNonBonded->getAtomStreamWidth( ),
+             (float) _brookNonBonded->getJStreamWidth( ),
+             (float) _brookNonBonded->getPartialForceStreamWidth( ),
+             epsfac,
+             dummyParameters,
+             positionStream.getBrookStream(),
+             _brookNonBonded->getChargeStream()->getBrookStream(),
+             _brookNonBonded->getOuterVdwStream()->getBrookStream(),
+             _brookNonBonded->getInnerSigmaStream()->getBrookStream(),
+             _brookNonBonded->getInnerEpsilonStream()->getBrookStream(),
+             _brookNonBonded->getExclusionStream()->getBrookStream(),
+             nonbondedForceStreams[0]->getBrookStream(),
+             nonbondedForceStreams[1]->getBrookStream(),
+             nonbondedForceStreams[2]->getBrookStream(),
+             nonbondedForceStreams[3]->getBrookStream()
+           );
+
+/*
+float zerof = 0.0f;
+nonbondedForceStreams[0]->fillWithValue( &zerof );
+nonbondedForceStreams[1]->fillWithValue( &zerof );
+nonbondedForceStreams[2]->fillWithValue( &zerof );
+nonbondedForceStreams[3]->fillWithValue( &zerof );
+*/
+
+   // diagnostics
+
+   if( PrintOn ){
+      (void) fprintf( getLog(), "\nPost knbforce_CDLJ4: atoms=%6d ceiling=%3d dupFac=%3d", _brookNonBonded->getNumberOfAtoms(),  
+                                                                                           _brookNonBonded->getAtomSizeCeiling(),
+                                                                                           _brookNonBonded->getDuplicationFactor()  );
+
+      (void) fprintf( getLog(), "\n                      hght=%6d   width=%3d   jWid=%3d", _brookNonBonded->getAtomStreamHeight( ),
+                                                                                           _brookNonBonded->getAtomStreamWidth( ),
+                                                                                           _brookNonBonded->getJStreamWidth( ) );
+      (void) fprintf( getLog(), "\n                      pFrc=%6d     eps=%12.5e\n",       _brookNonBonded->getPartialForceStreamWidth( ), epsfac );
+
+      (void) fprintf( getLog(), "\nOuterVdwStreamd\n" );
+      _brookNonBonded->getOuterVdwStream()->printToFile( getLog() );
+
+      (void) fprintf( getLog(), "\nInnerSigmaStream\n" );
+      _brookNonBonded->getInnerSigmaStream()->printToFile( getLog() );
+
+      (void) fprintf( getLog(), "\nInnerEpsilonStream\n" );
+      _brookNonBonded->getInnerEpsilonStream()->printToFile( getLog() );
+
+      (void) fprintf( getLog(), "\nExclusionStream\n" );
+      _brookNonBonded->getExclusionStream()->printToFile( getLog() );
+
+      (void) fprintf( getLog(), "\nChargeStream\n" );
+      _brookNonBonded->getChargeStream()->printToFile( getLog() );
+
+      for( int ii = 0; ii < 4; ii++ ){
+         (void) fprintf( getLog(), "\nForce stream %d\n", ii );
+         nonbondedForceStreams[ii]->printToFile( getLog() );
+      }
+   }
+
+   // gather forces
+
+   kMergeFloat3_4_nobranch( (float) _brookNonBonded->getDuplicationFactor(),
+                            (float) _brookNonBonded->getAtomStreamWidth(),
+                            (float) _brookNonBonded->getPartialForceStreamWidth(),
+                            (float) _brookNonBonded->getNumberOfAtoms(),
+                            (float) _brookNonBonded->getAtomSizeCeiling(),
+                            (float) _brookNonBonded->getOuterLoopUnroll(),
+                            nonbondedForceStreams[0]->getBrookStream(),
+                            nonbondedForceStreams[1]->getBrookStream(),
+                            nonbondedForceStreams[2]->getBrookStream(),
+                            nonbondedForceStreams[3]->getBrookStream(),
+                            forceStream.getBrookStream() );
+
+   // bonded
+
+         epsfac                                        = (float) (_brookBonded->getLJ_14Scale()*_brookBonded->getCoulombFactor());
    float width                                         = (float) (_brookBonded->getInverseMapStreamWidth());
 
    // bonded forces
@@ -253,6 +343,7 @@ void BrookCalcStandardMMForceFieldKernel::executeForces( const Stream& positions
                  (float) bondedForceStreams[0]->getStreamWidth(),
                  dummyParameters,
                  positionStream.getBrookStream(),
+                 _brookBonded->getChargeStream()->getBrookStream(),
                  _brookBonded->getAtomIndicesStream()->getBrookStream(),
                  bondedParameters[0]->getBrookStream(),
                  bondedParameters[1]->getBrookStream(),
@@ -264,6 +355,37 @@ void BrookCalcStandardMMForceFieldKernel::executeForces( const Stream& positions
                  bondedForceStreams[2]->getBrookStream(),
                  bondedForceStreams[3]->getBrookStream() );
 
+
+   // diagnostics
+
+   if( 1 && PrintOn ){
+
+      int countPrintInvMap[4] = { 3, 5, 2, 4 }; 
+
+      (void) fprintf( getLog(), "\nPost kbonded_CDLJ: epsFac=%.6f %.6f %.6f", epsfac, _brookBonded->getLJ_14Scale(), _brookBonded->getCoulombFactor());
+      (void) fprintf( getLog(), "\nAtom indices stream\n" );
+      _brookBonded->getAtomIndicesStream()->printToFile( getLog() );
+
+      (void) fprintf( getLog(), "\nCharge stream\n" );
+      _brookBonded->getChargeStream()->printToFile( getLog() );
+
+      for( int ii = 0; ii < 5; ii++ ){
+         (void) fprintf( getLog(), "\nParam stream %d\n", ii );
+         bondedParameters[ii]->printToFile( getLog() );
+      }
+      for( int ii = 0; ii < 4; ii++ ){
+         (void) fprintf( getLog(), "\nForce stream %d\n", ii );
+         bondedForceStreams[ii]->printToFile( getLog() );
+      }
+      (void) fprintf( getLog(), "\nInverse map streams\n" );
+      for( int ii = 0; ii < 4; ii++ ){
+         for( int jj = 0; jj < countPrintInvMap[ii]; jj++ ){
+            (void) fprintf( getLog(), "\n   Inverse map streams index=%d %d\n", ii, jj );
+            inverseStreamMaps[ii][jj]->printToFile( getLog() );
+         }
+      }
+   }
+
    // gather forces
 
    if( _brookBonded->getInverseMapStreamCount( K_Stream ) <= 4 ){
@@ -329,50 +451,6 @@ void BrookCalcStandardMMForceFieldKernel::executeForces( const Stream& positions
                       forceStream.getBrookStream(), forceStream.getBrookStream() );
 
    // ---------------------------------------------------------------------------------------
-
-   // nonbonded forces
-
-   BrookFloatStreamInternal**  nonbondedForceStreams = _brookNonBonded->getForceStreams();
-
-   knbforce_CDLJ4(
-             (float) _brookNonBonded->getNumberOfAtoms(),
-             (float) _brookNonBonded->getAtomSizeCeiling(),
-             (float) _brookNonBonded->getDuplicationFactor(),
-             (float) _brookNonBonded->getAtomStreamHeight( ),
-             (float) _brookNonBonded->getAtomStreamWidth( ),
-             (float) _brookNonBonded->getJStreamWidth( ),
-             (float) _brookNonBonded->getPartialForceStreamWidth( ),
-             epsfac,
-             dummyParameters,
-             positionStream.getBrookStream(),
-             _brookNonBonded->getOuterVdwStream()->getBrookStream(),
-             _brookNonBonded->getInnerSigmaStream()->getBrookStream(),
-             _brookNonBonded->getInnerEpsilonStream()->getBrookStream(),
-             _brookNonBonded->getExclusionStream()->getBrookStream(),
-             nonbondedForceStreams[0]->getBrookStream(),
-             nonbondedForceStreams[1]->getBrookStream(),
-             nonbondedForceStreams[2]->getBrookStream(),
-             nonbondedForceStreams[3]->getBrookStream()
-           );
-
-   // strF is assummed to be zero here
-
-   // gather forces
-
-   kMergeFloat3_4_nobranch( (float) _brookNonBonded->getDuplicationFactor(),
-                            (float) _brookNonBonded->getAtomStreamWidth(),
-                            (float) _brookNonBonded->getPartialForceStreamWidth(),
-                            (float) _brookNonBonded->getNumberOfAtoms(),
-                            (float) _brookNonBonded->getAtomSizeCeiling(),
-                            (float) _brookNonBonded->getOuterLoopUnroll(),
-                            nonbondedForceStreams[0]->getBrookStream(),
-                            nonbondedForceStreams[1]->getBrookStream(),
-                            nonbondedForceStreams[2]->getBrookStream(),
-                            nonbondedForceStreams[3]->getBrookStream(),
-                            forceStream.getBrookStream() );
-
-
-   // ---------------------------------------------------------------------------------------
 }
 
 /**

platforms/brook/src/BrookFloatStreamInternal.cpp

@@ -54,7 +54,6 @@ BrookFloatStreamInternal::BrookFloatStreamInternal( const std::string& name, int
 // ---------------------------------------------------------------------------------------
 
    static const std::string methodName      = "BrookFloatStreamInternal::BrookFloatStreamInternal";
-   // static const int debug                   = 1;
 
 // ---------------------------------------------------------------------------------------
 
@@ -119,6 +118,7 @@ BrookFloatStreamInternal::BrookFloatStreamInternal( const std::string& name, int
    // set stream height based on specified stream _width
 
    if( streamWidth < 1 ){
+
       std::stringstream message;
       message << methodName << " stream=" << name << " input stream width=" << streamWidth << " is less than 1.";
       throw OpenMMException( message.str() );
@@ -155,6 +155,9 @@ BrookFloatStreamInternal::BrookFloatStreamInternal( const std::string& name, int
 
    _data = new float[streamSize*_width];
 
+//printf( "%s %s data=%u stream=%d [%d %d] width=%d\n", methodName.c_str(), getName().c_str(), (unsigned int) _data, streamSize, _streamHeight, _streamWidth, _width );
+//fflush( stdout );
+
 }
 
 /** 
@@ -166,12 +169,15 @@ BrookFloatStreamInternal::~BrookFloatStreamInternal( ){
 
 // ---------------------------------------------------------------------------------------
 
-   //static const std::string methodName      = "BrookFloatStreamInternal::~BrookFloatStreamInternal";
-   // static const int debug                   = 1;
+   static const std::string methodName      = "BrookFloatStreamInternal::~BrookFloatStreamInternal";
 
 // ---------------------------------------------------------------------------------------
 
    //delete _aStream;
+
+//printf( "%s %s data=%u stream=%d [%d %d] width=%d\n", methodName.c_str(), getName().c_str(), (unsigned int) _data, getStreamSize(), _streamHeight, _streamWidth, _width );
+//fflush( stdout );
+
    delete[] _data;
 
 }
@@ -218,7 +224,7 @@ void BrookFloatStreamInternal::loadFromArray( const void* array ){
  * the values should be packed into a single array: all the values for the first element, followed by all the values
  * for the next element, etc.
  *
- * @throw exception if baseType not float or double
+ * @throw exception if baseType not float, double, or integer
  *
  */
 
@@ -227,11 +233,11 @@ void BrookFloatStreamInternal::loadFromArray( const void* array, BrookStreamInte
 // ---------------------------------------------------------------------------------------
 
    static const std::string methodName      = "BrookFloatStreamInternal::loadFromArray";
-   // static const int debug                   = 1;
 
 // ---------------------------------------------------------------------------------------
 
    int totalSize                           = getSize()*getWidth();
+   //int totalSize                           = getSize();
 
    if( baseType == BrookStreamInternal::Float ){
 
@@ -293,6 +299,8 @@ void BrookFloatStreamInternal::saveToArray( void* array ){
 
    if( baseType == BrookStreamInternal::Float ){
 
+//printf( "%s Basetype is float\n", methodName.c_str() );
+//fflush( stdout );
       memcpy( array, _data, sizeof( float )*totalSize );
 /*
       float* arrayData = (float*) array;
@@ -303,6 +311,9 @@ void BrookFloatStreamInternal::saveToArray( void* array ){
 
    } else if( baseType == BrookStreamInternal::Double ){
 
+//printf( "%s Basetype is double\n", methodName.c_str() );
+//fflush( stdout );
+
       double* arrayData = (double*) array;
       for( int ii = 0; ii < totalSize; ii++ ){
          arrayData[ii] = (double) _data[ii];
@@ -310,6 +321,9 @@ void BrookFloatStreamInternal::saveToArray( void* array ){
 
    } else if( baseType == BrookStreamInternal::Integer ){
 
+//printf( "%s Basetype is int\n", methodName.c_str() );
+//fflush( stdout );
+
       int* arrayData = (int*) array;
       for( int ii = 0; ii < totalSize; ii++ ){
          arrayData[ii] = (int) _data[ii];
@@ -410,8 +424,7 @@ void BrookFloatStreamInternal::_loadDanglingValues( float danglingValue ){
 
 // ---------------------------------------------------------------------------------------
 
-   // static const std::string methodName      = "BrookFloatStreamInternal::_loadDanglingValues";
-   // static const int debug                   = 1;
+    static const std::string methodName      = "BrookFloatStreamInternal::_loadDanglingValues";
 
 // ---------------------------------------------------------------------------------------
 
@@ -420,6 +433,9 @@ void BrookFloatStreamInternal::_loadDanglingValues( float danglingValue ){
    int arraySize  = getSize()*width;
    int streamSize = getStreamSize()*width;
 
+//printf( "%s array=%d stream=%d width=%d %s\n", methodName.c_str(), arraySize, streamSize, width, getName().c_str() );
+//fflush( stdout );
+
    if( arraySize < streamSize ){
       for( int ii = arraySize; ii < streamSize; ii++ ){
           _data[ii] = danglingValue;
@@ -489,3 +505,68 @@ const std::string BrookFloatStreamInternal::getContentsString( int level ) const
    return message.str();
 }
 
+/* 
+ * Print array contents of object to file
+ *
+ * @param log         file to print to
+ *
+ * @return DefaultReturnValue
+ *
+ * */
+
+int BrookFloatStreamInternal::_bodyPrintToFile( FILE* log ){
+
+// ---------------------------------------------------------------------------------------
+
+   //static const std::string methodName      = "BrookStreamInternal::_bodyPrintToFile";
+
+// ---------------------------------------------------------------------------------------
+
+   void* dataArrayV = getDataArray( );
+   saveToArray( dataArrayV );
+
+   int streamSize   = getStreamSize();
+   int width        = getWidth();
+   int index        = 0;
+    float* dataArray = (float*) dataArrayV;
+   for( int ii = 0; ii < streamSize; ii++ ){
+      std::stringstream message;
+      message.width( 15 );
+      message.precision( 5 );
+      message << ii << " [ ";
+      for( int jj = 0; jj < width; jj++ ){
+         message << dataArray[index++] << " ";
+      }
+      message << "]\n";
+      (void) fprintf( log, "%s", message.str().c_str() );    
+   }
+
+   delete[] dataArrayV;
+
+   return DefaultReturnValue;
+
+}
+
+/** 
+ * Get array of appropritate size for loading data
+ *
+ * @return data array -- user's responsibility to free
+ */
+
+void* BrookFloatStreamInternal::getDataArray( void ){
+
+// ---------------------------------------------------------------------------------------
+
+   //static const std::string methodName      = "BrookStreamInternal::getDataArray";
+
+// ---------------------------------------------------------------------------------------
+
+   int totalSize                          = getStreamSize()*getWidth();
+   BrookStreamInternal::DataType baseType = getBaseDataType();
+   
+   if( baseType == Double || baseType == Double2 ||  baseType == Double3 ||  baseType == Double4 ){
+      totalSize *= 2;
+   }
+   return new float[totalSize];
+}
+

platforms/brook/src/BrookFloatStreamInternal.h

@@ -132,6 +132,13 @@ class BrookFloatStreamInternal : public BrookStreamInternal {
       
       void* getData( int readFromBoard );
 
+      /** 
+       * Get array of appropritate size for loading data
+       *
+       * @return data array -- user's responsibility to free
+       */
+      void* getDataArray( void );
+  
       /** 
        * Get dangle value
        * 
@@ -163,6 +170,16 @@ class BrookFloatStreamInternal : public BrookStreamInternal {
       void _loadDanglingValues( void  );
       void _loadDanglingValues( float );
 
+      /*  
+       * Print array to file
+       *
+       * @param log  log file
+       *
+       * @return  DefaultReturnValue
+       *
+       * */
+
+      int _bodyPrintToFile( FILE* log );
 };
 
 } // namespace OpenMM

platforms/brook/src/BrookIntStreamInternal.cpp

@@ -194,6 +194,24 @@ void BrookIntStreamInternal::fillWithValue( void* value ){
    }
 }
 
+/** 
+ * Get array of appropritate size for loading data
+ *
+ * @return data array -- user's responsibility to free
+ */
+
+void* BrookIntStreamInternal::getDataArray( void ){
+
+// ---------------------------------------------------------------------------------------
+
+   //static const std::string methodName      = "BrookIntStreamInternal::getDataArray";
+
+// ---------------------------------------------------------------------------------------
+
+   int totalSize                          = getStreamSize()*getWidth();
+   return new int[totalSize];
+}  
+
 /** 
  * Get data
  * 
@@ -205,3 +223,119 @@ void* BrookIntStreamInternal::getData( void ){
    return _data;
 }
 
+/** 
+ * Get data
+ * 
+ * @param readFromBoard if set, read values on board 
+ *
+ * @return data array
+ *
+ */
+
+void* BrookIntStreamInternal::getData( int readFromBoard ){
+
+// ---------------------------------------------------------------------------------------
+
+   // static const std::string methodName      = "BrookIntStreamInternal::getData";
+
+// ---------------------------------------------------------------------------------------
+
+   if( readFromBoard ){
+      _aStream.write( _data );
+   }
+
+   return (void*) _data;
+}
+
+/* 
+ * Print array contents of object to file
+ *
+ * @param log         file to print to
+ *
+ * @return DefaultReturnValue
+ *
+ * */
+
+int BrookIntStreamInternal::_bodyPrintToFile( FILE* log ){
+
+// ---------------------------------------------------------------------------------------
+
+   //static const std::string methodName      = "BrookIntStreamInternal::_bodyPrintToFile";
+
+// ---------------------------------------------------------------------------------------
+
+   void* dataArrayV = getDataArray( );
+   saveToArray( dataArrayV );
+
+   int streamSize   = getStreamSize();
+   int width        = getWidth();
+   int index        = 0;
+   int* dataArray   = (int*) dataArrayV;
+   for( int ii = 0; ii < streamSize; ii++ ){
+      std::stringstream message;
+      message.width( 10 );
+      message << ii << " [ ";
+      for( int jj = 0; jj < width; jj++ ){
+         message << dataArray[index++] << " ";
+      }   
+      message << "]\n";
+      (void) fprintf( log, "%s", message.str().c_str() );    
+   }   
+
+   delete[] dataArrayV;
+
+   return DefaultReturnValue;
+
+}
+
+/* 
+ * Get contents of object
+ *
+ * @param level   level of dump
+ *
+ * @return string containing contents
+ *
+ * */
+
+const std::string BrookIntStreamInternal::getContentsString( int level ) const {
+
+// ---------------------------------------------------------------------------------------
+
+   // static const std::string methodName      = "BrookIntStreamInternal::getContentsString";
+
+   static const unsigned int MAX_LINE_CHARS = 256;
+   char value[MAX_LINE_CHARS];
+   //static const char* Set                   = "Set";
+   //static const char* NotSet                = "Not set";
+
+
+// ---------------------------------------------------------------------------------------
+
+   std::stringstream message;
+   std::string tab   = "   ";
+
+#ifdef WIN32
+#define LOCAL_SPRINTF(a,b,c) sprintf_s( (a), MAX_LINE_CHARS, (b), (c) );   
+#else
+#define LOCAL_SPRINTF(a,b,c) sprintf( (a), (b), (c) );   
+#endif
+
+   (void) LOCAL_SPRINTF( value, "%s", getName().c_str() );
+   message << _getLine( tab, "Name:", value );
+
+   (void) LOCAL_SPRINTF( value, "%d", getWidth() );
+   message << _getLine( tab, "Width:", value );
+
+   (void) LOCAL_SPRINTF( value, "%d", getStreamSize() );
+   message << _getLine( tab, "Stream size:", value );
+
+   (void) LOCAL_SPRINTF( value, "%d", getStreamWidth() );
+   message << _getLine( tab, "Stream width:", value );
+
+   (void) LOCAL_SPRINTF( value, "%d", getStreamHeight() );
+   message << _getLine( tab, "Stream height:", value );
+
+   return message.str();
+}
+
+

platforms/brook/src/BrookIntStreamInternal.h

@@ -123,11 +123,53 @@ public:
 
       void* getData( void );
 
+      /** 
+       * Get data
+       * 
+       * @param readFromBoard if set, read values on board 
+       *
+       * @return data array
+       *
+       */
+      
+      void* getData( int readFromBoard );
+      
+      /** 
+       * Get array of appropritate size for loading data
+       *
+       * @return data array -- user's responsibility to free
+       */
+      
+      void* getDataArray( void );
+      
+      /*  
+       * Get contents of object
+       *
+       *
+       * @param level   level of dump
+       *
+       * @return string containing contents
+       *
+       * */
+
+      const std::string getContentsString( int level = 0 ) const;
+
 private:
 
     int _dangleValue;
 
     int* _data;
+
+      /*  
+       * Print array to file
+       *
+       * @param log  log file
+       *
+       * @return  DefaultReturnValue
+       *
+       * */
+
+      int _bodyPrintToFile( FILE* log );
 };
 
 } // namespace OpenMM

platforms/brook/src/BrookIntegrateLangevinStepKernel.cpp

@@ -103,7 +103,7 @@ void BrookIntegrateLangevinStepKernel::execute( Stream& positions, Stream& veloc
 // ---------------------------------------------------------------------------------------
 
    double epsilon = 1.0e-04;
-   //static const std::string methodName      = "BrookIntegrateLangevinStepKernel::execute";
+   static const std::string methodName      = "BrookIntegrateLangevinStepKernel::execute";
 
 // ---------------------------------------------------------------------------------------
    
@@ -121,7 +121,17 @@ void BrookIntegrateLangevinStepKernel::execute( Stream& positions, Stream& veloc
    if( fabs( differences[0] ) < epsilon || fabs( differences[1] ) < epsilon || fabs( differences[2] ) < epsilon ){
       _brookStochasticDynamics->updateParameters( temperature, friction, stepSize );
    } 
-   _brookStochasticDynamics->update( positions, velocities, forces, (_brookShakeAlgorithm ? *_brookShakeAlgorithm : NULL),
-                                     *_brookRandomNumberGenerator );
+
+   assert( _brookShakeAlgorithm );
+/*
+   if( _brookShakeAlgorithm == NULL ){
+      std::stringstream message;
+      message << methodName << " _brookShakeAlgorithm is not set -- case not handled.";
+      throw OpenMMException( message.str() );
+      return NULL;
+   }
+*/
+
+   _brookStochasticDynamics->update( positions, velocities, forces, *_brookShakeAlgorithm, *_brookRandomNumberGenerator );
 
 }

platforms/brook/src/BrookIntegrateLangevinStepKernel.h

@@ -86,7 +86,7 @@ class BrookIntegrateLangevinStepKernel : public IntegrateLangevinStepKernel {
        *
        */
 
-      void execute( Stream& positions, Stream& velocities, const Stream& forces, double temperature, double friction, double stepSize);
+      void execute( Stream& positions, Stream& velocities, const Stream& forces, double temperature, double friction, double stepSize );
 
    protected:
 

platforms/brook/src/BrookNonBonded.cpp

@@ -50,7 +50,6 @@ BrookNonBonded::BrookNonBonded(  ){
 // ---------------------------------------------------------------------------------------
 
    // static const std::string methodName      = "BrookNonBonded::BrookNonBonded";
-   // static const int debug                   = 1;
 
 // ---------------------------------------------------------------------------------------
 
@@ -177,6 +176,7 @@ int BrookNonBonded::getAtomSizeCeiling( void ) const {
       if( localThis->_atomSizeCeiling ){
          localThis->_atomSizeCeiling = localThis->getOuterLoopUnroll() - localThis->_atomSizeCeiling;
       }   
+      localThis->_atomSizeCeiling += localThis->getNumberOfAtoms();
    }
 
    return _atomSizeCeiling;
@@ -277,7 +277,6 @@ int BrookNonBonded::getPartialForceStreamWidth( void ) const {
 // ---------------------------------------------------------------------------------------
 
    // static const std::string methodName      = "BrookNonBonded::getPartialForceStreamWidth";
-   // static const int debug                   = 1;
 
 // ---------------------------------------------------------------------------------------
 
@@ -473,11 +472,11 @@ int BrookNonBonded::isForceStreamSet( int index ) const {
  *
  */
 
-int BrookNonBonded::initializeStreamSizes( int numberOfAtoms, const Platform& platform ){
+int BrookNonBonded::_initializeStreamSizes( int numberOfAtoms, const Platform& platform ){
 
 // ---------------------------------------------------------------------------------------
 
-   static const std::string methodName      = "BrookNonBonded::initializeStreamSizes";
+   static const std::string methodName      = "BrookNonBonded::_initializeStreamSizes";
    // static const int debug                   = 1;
 
 // ---------------------------------------------------------------------------------------
@@ -485,10 +484,10 @@ int BrookNonBonded::initializeStreamSizes( int numberOfAtoms, const Platform& pl
    int atomStreamSize                  = getAtomStreamSize( platform );
    int atomStreamWidth                 = getAtomStreamWidth( platform );
 
-   initializeExclusionStreamSize( atomStreamSize, atomStreamWidth );
-   initializeJStreamSize( atomStreamSize, atomStreamWidth );
-   initializeOuterVdwStreamSize( atomStreamSize, atomStreamWidth );
-   initializePartialForceStreamSize( atomStreamSize, atomStreamWidth );
+   _initializeExclusionStreamSize( atomStreamSize, atomStreamWidth );
+   _initializeJStreamSize( atomStreamSize, atomStreamWidth );
+   _initializeOuterVdwStreamSize( atomStreamSize, atomStreamWidth );
+   _initializePartialForceStreamSize( atomStreamSize, atomStreamWidth );
 
    return DefaultReturnValue;
 }
@@ -502,11 +501,11 @@ int BrookNonBonded::initializeStreamSizes( int numberOfAtoms, const Platform& pl
  *
  */
 
-int BrookNonBonded::initializeStreams( const Platform& platform ){
+int BrookNonBonded::_initializeStreams( const Platform& platform ){
 
 // ---------------------------------------------------------------------------------------
 
-   static const std::string methodName      = "BrookNonBonded::initializeStreams";
+   static const std::string methodName      = "BrookNonBonded::_initializeStreams";
    static const double dangleValue          = 0.0;
    // static const int debug                   = 1;
 
@@ -539,13 +538,14 @@ int BrookNonBonded::initializeStreams( const Platform& platform ){
                                                                                 getAtomStreamWidth(), BrookStreamInternal::Float, dangleValue );
 
    
-   // partial force stream
+   // partial force streams
+
    std::string partialForceStream = BrookCommon::PartialForceStream;
    for( int ii = 0; ii < getNumberOfForceStreams(); ii++ ){
       std::stringstream name;
       name << partialForceStream << ii;
       _nonbondedForceStreams[ii] = new BrookFloatStreamInternal( name.str(), getPartialForceStreamSize(),
-                                                                 getPartialForceStreamWidth(), BrookStreamInternal::Float4, dangleValue );
+                                                                 getPartialForceStreamWidth(), BrookStreamInternal::Float3, dangleValue );
    }
 
    return DefaultReturnValue;
@@ -563,12 +563,11 @@ int BrookNonBonded::initializeStreams( const Platform& platform ){
  *
  */
 
-int BrookNonBonded::setExclusion( int i, int j, int exclusionStreamWidth, BrookOpenMMFloat* exclusion ){
+int BrookNonBonded::_setExclusion( int i, int j, int exclusionStreamWidth, BrookOpenMMFloat* exclusion ){
 
 // ---------------------------------------------------------------------------------------
 
-   static const std::string methodName      = "BrookNonBonded::setExclusion";
-   // static const int debug                   = 1;
+   static const std::string methodName      = "BrookNonBonded::_setExclusion";
 
 // ---------------------------------------------------------------------------------------
 
@@ -612,23 +611,23 @@ int BrookNonBonded::setExclusion( int i, int j, int exclusionStreamWidth, BrookO
  *
  */
 
-int BrookNonBonded::initializeExclusions( const std::vector<std::set<int> >& exclusionsVector, const Platform& platform ){
+int BrookNonBonded::_initializeExclusions( const std::vector<std::set<int> >& exclusionsVector, const Platform& platform ){
 
 // ---------------------------------------------------------------------------------------
 
-   static const std::string methodName      = "BrookNonBonded::initializeExclusions";
+   static const std::string methodName      = "BrookNonBonded::_initializeExclusions";
    static const int debug                   = 1;
 
 // ---------------------------------------------------------------------------------------
 
    if( debug && getLog() ){
-      (void) fprintf( getLog(), "In initializeExclusions exclusions vector size=%d\n", exclusionsVector.size() );
+      (void) fprintf( getLog(), "%s exclusions vector size=%d\n", methodName.c_str(), exclusionsVector.size() );
    }
 
    int exclusionStreamSize = getExclusionStreamSize();
    if( exclusionStreamSize < 1 ){
       std::stringstream message;
-      message << "BrookNonBonded::initializeExclusions exclusionStreamSize=" << exclusionStreamSize << " is not set.";
+      message << methodName << " exclusionStreamSize=" << exclusionStreamSize << " is not set.";
       throw OpenMMException( message.str() );
       return ErrorReturnValue;
    }
@@ -637,7 +636,7 @@ int BrookNonBonded::initializeExclusions( const std::vector<std::set<int> >& exc
 
    BrookOpenMMFloat* exclusions = new BrookOpenMMFloat[exclusionStreamSize];
    for( unsigned int ii = 0; ii < (unsigned int) exclusionStreamSize; ii++ ){
-      exclusions[ii] = (BrookOpenMMFloat) 210.0;
+      exclusions[ii] = (BrookOpenMMFloat) 210.0f;
    }
 
    // pack in values
@@ -648,14 +647,22 @@ int BrookNonBonded::initializeExclusions( const std::vector<std::set<int> >& exc
    for( unsigned int ii = 0; ii < exclusionsVector.size(); ii++ ){
 
       set<int> exclusionIndices  = exclusionsVector[ii];
+
+//(void) fprintf( getLog(), "%s atoms=%d excludeSz=%d\n", methodName.c_str(), ii, exclusionIndices.size() );
+
+      int hitII = 0;
       for( set<int>::const_iterator jj = exclusionIndices.begin(); jj != exclusionIndices.end(); jj++ ){
-         setExclusion( ii, *jj, exclusionStreamWidth, exclusions );
+//(void) fprintf( getLog(), "%s atoms=%d exclude=%d\n", methodName.c_str(), ii, *jj );
+         _setExclusion( ii, *jj, exclusionStreamWidth, exclusions );
+         if( *jj == ii )hitII = 1;
+      }
+      if( !hitII ){
+         _setExclusion( ii, ii, exclusionStreamWidth, exclusions );
       }
 
       // explicitly exclude junk atoms from interacting w/ atom ii
-
       for( int jj = numberOfAtoms; jj < atomStreamSize; jj++ ){
-         setExclusion( ii, jj, exclusionStreamWidth, exclusions );
+         _setExclusion( ii, jj, exclusionStreamWidth, exclusions );
       }
    }
 
@@ -663,8 +670,40 @@ int BrookNonBonded::initializeExclusions( const std::vector<std::set<int> >& exc
 
    for( int ii = numberOfAtoms; ii < atomStreamSize; ii++ ){
       for( int jj = 0; jj < atomStreamSize; jj++ ){
-         setExclusion( ii, jj, exclusionStreamWidth, exclusions );
+         _setExclusion( ii, jj, exclusionStreamWidth, exclusions );
+      }
+   }
+
+   // diagnostics
+
+   if( 1 && debug && getLog() ){
+      FILE* log = getLog();
+      (void) fprintf( log, "%s atoms=%d excl=%d exStrSz=%d w=%d atmStrSz=%d\n", methodName.c_str(), numberOfAtoms, exclusionsVector.size(), 
+                      exclusionStreamSize, exclusionStreamWidth, atomStreamSize );
+/*
+      for( int ii = 0; ii < exclusionStreamSize; ii++ ){
+int index  = ii/exclusionStreamWidth;
+int offset = ii - index*exclusionStreamWidth;
+         (void) fprintf( log, "   %6d %6d [%6d %6d] %10.0f\n", ii, offset, index*4, index*4+3, exclusions[ii] );
       }
+*/
+
+      for( int ii = 0; ii < numberOfAtoms; ii++ ){
+         (void) fprintf( log, "%6d ", ii );
+         int count = 0;
+         for( int jj = 0; jj <= numberOfAtoms/4; jj++ ){
+            int index = jj*exclusionStreamWidth + ii;
+            (void) fprintf( log, " [%4d %4d] %5.1f", jj*4, jj*4+3, exclusions[index] );
+            int excludeValue = (int) (exclusions[index] + 0.01);
+            if( (excludeValue % 2) )count++;
+            if( (jj*4+1) < numberOfAtoms && (excludeValue % 3) )count++;
+            if( (jj*4+2) < numberOfAtoms && (excludeValue % 5) )count++;
+            if( (jj*4+3) < numberOfAtoms && (excludeValue % 7) )count++;
+         }
+         (void) fprintf( log, " TtlExcl=%d\n", count );
+      }
+         
+      (void) fflush( log );
    }
 
    // load stream
@@ -687,12 +726,11 @@ int BrookNonBonded::initializeExclusions( const std::vector<std::set<int> >& exc
  *
  */
 
-int BrookNonBonded::setSigmaEpsilon( double c6, double c12, double* sigma , double* epsilon ){
+int BrookNonBonded::_setSigmaEpsilon( double c6, double c12, double* sigma , double* epsilon ){
 
 // ---------------------------------------------------------------------------------------
 
-   // static const std::string methodName      = "BrookNonBonded::setSigmaEpsilon";
-   // static const int debug                   = 1;
+   // static const std::string methodName      = "BrookNonBonded::_setSigmaEpsilon";
 
 // ---------------------------------------------------------------------------------------
 
@@ -700,8 +738,10 @@ int BrookNonBonded::setSigmaEpsilon( double c6, double c12, double* sigma , doub
       *sigma    = 1.0;
       *epsilon  = 0.0;
    } else {
-      *sigma    = 0.5*pow( c12/c6, 1.0/6.0 );
-      *epsilon  = sqrt( c6*c6/c12 );
+      //*sigma    = 0.5*pow( c12/c6, 1.0/6.0 );
+      //*epsilon  = sqrt( c6*c6/c12 );
+      *sigma    = 0.5*c6;
+      *epsilon  = 2.0*sqrt( c12 );
    }
 
    return DefaultReturnValue;
@@ -717,17 +757,17 @@ int BrookNonBonded::setSigmaEpsilon( double c6, double c12, double* sigma , doub
  *
  */
 
-int BrookNonBonded::initializeVdwAndCharge( const vector<vector<double> >& nonbondedParameters, const Platform& platform ){
+int BrookNonBonded::_initializeVdwAndCharge( const vector<vector<double> >& nonbondedParameters, const Platform& platform ){
 
 // ---------------------------------------------------------------------------------------
 
-   static const std::string methodName      = "BrookNonBonded::initializeVdwAndCharge";
+   static const std::string methodName      = "BrookNonBonded::_initializeVdwAndCharge";
    static const int debug                   = 1;
 
 // ---------------------------------------------------------------------------------------
 
    if( debug && getLog() ){
-      (void) fprintf( getLog(), "In initializeVdwAndCharge nonbonded vector size=%d\n", nonbondedParameters.size() );
+      (void) fprintf( getLog(), "%s nonbonded vector size=%d\n", methodName.c_str(), nonbondedParameters.size() );
    }
 
    int atomStreamSize        = getAtomStreamSize();
@@ -741,23 +781,31 @@ int BrookNonBonded::initializeVdwAndCharge( const vector<vector<double> >& nonbo
    BrookOpenMMFloat* charges       = new BrookOpenMMFloat[atomStreamSize];
    memset( charges, 0, sizeof( BrookOpenMMFloat )*atomStreamSize );
    
+// ---------------------------------------------------------------------------------------
+
    // pack in values
 
    int  numberOfAtoms        = getNumberOfAtoms();
-   double sigma, epsilon;
    int trackingIndex         = 0;
+   double sigma, epsilon;
    for( unsigned int ii = 0; ii < nonbondedParameters.size(); ii++ ){
 
       vector<double> nonbondedParameterVector 
                                      = nonbondedParameters[ii];
 
-      double c6                      = nonbondedParameterVector[0];
-      double c12                     = nonbondedParameterVector[1];
-      double charge                  = nonbondedParameterVector[2];
+      double c6                      = nonbondedParameterVector[1];
+      double c12                     = nonbondedParameterVector[2];
+      double charge                  = nonbondedParameterVector[0];
+
+/*
+      double charge                  = nonbondedParameterVector[0];
+      double sigma                   = nonbondedParameterVector[1];
+      double epsilon                 = nonbondedParameterVector[2];
+*/
 
       // geometric combination rules
 
-      setSigmaEpsilon( c6, c12, &sigma, &epsilon );
+      _setSigmaEpsilon( c6, c12, &sigma, &epsilon );
 
       vdwParameters[trackingIndex++] = (BrookOpenMMFloat) sigma;
       vdwParameters[trackingIndex++] = (BrookOpenMMFloat) epsilon;
@@ -772,17 +820,126 @@ int BrookNonBonded::initializeVdwAndCharge( const vector<vector<double> >& nonbo
       vdwParameters[trackingIndex++] = (BrookOpenMMFloat) 0.0;
    }
 
+// ---------------------------------------------------------------------------------------
+
    // load stream
 
    _nonbondedStreams[OuterVdwStream]->loadFromArray( vdwParameters ); 
    _nonbondedStreams[ChargeStream]->loadFromArray( charges ); 
 
+   // diagnostics
+
+   if( 1 && debug && getLog() ){
+      FILE* log = getLog();
+      int trackingIndex         = 0;
+      (void) fprintf( log, "%s atoms=%d strSz=%d\n   vdw[Sig Eps] Q\n", methodName.c_str(), numberOfAtoms, atomStreamSize );
+      for( int ii = 0; ii < atomStreamSize; ii++, trackingIndex += 2 ){
+         (void) fprintf( log, "   %d %10.3f %10.3f %10.3f\n", ii, vdwParameters[trackingIndex], vdwParameters[trackingIndex+1], charges[ii] );
+      }
+      (void) fflush( log );
+   }
+
    delete[] vdwParameters;
    delete[] charges;
 
    return DefaultReturnValue;
 }
 
+
+/** 
+ * Initialize vdw & charge
+ * 
+ * @param exclusions                vector of sets containing exclusions (1 set entry for every atom)
+ * @param platform                  platform
+ *
+ * @return ErrorReturnValue if error, else DefaultReturnValue
+ *
+ */
+
+int BrookNonBonded::_initializeJStreamVdw( const vector<vector<double> >& nonbondedParameters, const Platform& platform ){
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName      = "BrookNonBonded::_initializeJStreamVdw";
+   static const int debug                   = 1;
+
+// ---------------------------------------------------------------------------------------
+
+   if( debug && getLog() ){
+      (void) fprintf( getLog(), "%s nonbonded vector size=%d\n", methodName.c_str(), nonbondedParameters.size() );
+   }
+
+   int jStreamSize        = getJStreamSize();
+
+   // allocate and initialize vdw & charge array
+
+   unsigned int jParametersSize        = 4*jStreamSize;
+
+   BrookOpenMMFloat* sigmaParameters   = new BrookOpenMMFloat[jParametersSize];
+   //memset( sigmaParameters, 0, sizeof( BrookOpenMMFloat )*jParametersSize );
+   
+   BrookOpenMMFloat* epsilonParameters = new BrookOpenMMFloat[jParametersSize];
+   //memset( epsilonParameters, 0, sizeof( BrookOpenMMFloat )*jParametersSize );
+   
+// ---------------------------------------------------------------------------------------
+
+   // pack in values
+
+   int  numberOfAtoms        = getNumberOfAtoms();
+   int trackingIndex         = 0;
+   double sigma, epsilon;
+   for( unsigned int ii = 0; ii < nonbondedParameters.size(); ii++ ){
+
+      vector<double> nonbondedParameterVector 
+                                     = nonbondedParameters[ii];
+
+      double c6                      = nonbondedParameterVector[1];
+      double c12                     = nonbondedParameterVector[2];
+
+      //double sigma                   = nonbondedParameterVector[1];
+      //double epsilon                 = nonbondedParameterVector[2];
+
+      // geometric combination rules
+
+      _setSigmaEpsilon( c6, c12, &sigma, &epsilon );
+
+      sigmaParameters[ii]   = (BrookOpenMMFloat) sigma;
+      epsilonParameters[ii] = (BrookOpenMMFloat) epsilon;
+         
+   }
+
+   // set outlier atoms
+
+   for( unsigned int ii = nonbondedParameters.size(); ii < (unsigned int) jParametersSize; ii++ ){
+      sigmaParameters[ii]   = (BrookOpenMMFloat) 1.0;
+      epsilonParameters[ii] = (BrookOpenMMFloat) 0.0;
+   }
+
+// ---------------------------------------------------------------------------------------
+
+   // load streams
+
+   _nonbondedStreams[InnerSigmaStream]->loadFromArray( sigmaParameters ); 
+   _nonbondedStreams[InnerEpsilonStream]->loadFromArray( epsilonParameters ); 
+
+   // diagnostics
+
+   if( 1 && debug && getLog() ){
+      FILE* log = getLog();
+      int trackingIndex         = 0;
+      (void) fprintf( log, "%s atoms=%d strSz=%d\n   Innervdw[Sig Eps]\n", methodName.c_str(), numberOfAtoms, jStreamSize );
+      for( unsigned int ii = 0; ii < jParametersSize; ii++ ){
+         (void) fprintf( log, "   %d %10.3f %10.3f\n", ii, sigmaParameters[ii], epsilonParameters[ii] );
+      }
+      (void) fflush( log );
+   }
+
+   delete[] sigmaParameters;
+   delete[] epsilonParameters;
+
+   return DefaultReturnValue;
+}
+
 /* 
  * Setup of nonbonded ixns
  *
@@ -808,12 +965,12 @@ int BrookNonBonded::setup( int numberOfAtoms, const std::vector<std::vector<doub
 
    setNumberOfAtoms( numberOfAtoms );
 
-   initializeStreamSizes( numberOfAtoms, platform );
-   initializeStreams( platform );
+   _initializeStreamSizes( numberOfAtoms, platform );
+   _initializeStreams( platform );
 
-   initializeExclusions( exclusions, platform );
-   initializeVdwAndCharge( nonbondedParameters, platform );
-   //initializeJStreamVdw( nonbondedParameters, platform );
+   _initializeExclusions( exclusions, platform );
+   _initializeVdwAndCharge( nonbondedParameters, platform );
+   _initializeJStreamVdw( nonbondedParameters, platform );
 
    return DefaultReturnValue;
 }
@@ -828,11 +985,11 @@ int BrookNonBonded::setup( int numberOfAtoms, const std::vector<std::vector<doub
  *
  * */
 
-int BrookNonBonded::initializeExclusionStreamSize( int atomStreamSize, int atomStreamWidth ){
+int BrookNonBonded::_initializeExclusionStreamSize( int atomStreamSize, int atomStreamWidth ){
 
 // ---------------------------------------------------------------------------------------
 
-   static const std::string methodName      = "BrookNonBonded::initializeExclusionStreamSize";
+   static const std::string methodName      = "BrookNonBonded::_initializeExclusionStreamSize";
    //static const int debug                   = 1;
 
 // ---------------------------------------------------------------------------------------
@@ -864,11 +1021,11 @@ int BrookNonBonded::initializeExclusionStreamSize( int atomStreamSize, int atomS
  *
  * */
 
-int BrookNonBonded::initializeJStreamSize( int atomStreamSize, int atomStreamWidth ){
+int BrookNonBonded::_initializeJStreamSize( int atomStreamSize, int atomStreamWidth ){
 
 // ---------------------------------------------------------------------------------------
 
-   static const std::string methodName      = "BrookNonBonded::initializeJStreamSize";
+   static const std::string methodName      = "BrookNonBonded::_initializeJStreamSize";
 
 // ---------------------------------------------------------------------------------------
 
@@ -912,7 +1069,7 @@ int BrookNonBonded::initializeJStreamSize( int atomStreamSize, int atomStreamWid
  *
  * */
 
-int BrookNonBonded::initializeOuterVdwStreamSize( int atomStreamSize, int atomStreamWidth ){
+int BrookNonBonded::_initializeOuterVdwStreamSize( int atomStreamSize, int atomStreamWidth ){
 
 // ---------------------------------------------------------------------------------------
 
@@ -959,12 +1116,11 @@ int BrookNonBonded::initializeOuterVdwStreamSize( int atomStreamSize, int atomSt
  *
  * */
 
-int BrookNonBonded::initializePartialForceStreamSize( int atomStreamSize, int atomStreamWidth ){
+int BrookNonBonded::_initializePartialForceStreamSize( int atomStreamSize, int atomStreamWidth ){
 
 // ---------------------------------------------------------------------------------------
 
-   static const std::string methodName      = "BrookNonBonded::initializePartialForceStreamSize";
-   //static const int debug                   = 1;
+   static const std::string methodName      = "BrookNonBonded::_initializePartialForceStreamSize";
 
 // ---------------------------------------------------------------------------------------
 
@@ -991,16 +1147,13 @@ int BrookNonBonded::initializePartialForceStreamSize( int atomStreamSize, int at
 }
 
 /* 
- * Setup of j-stream dimensions
- * 
  * Get contents of object
  *
- *
  * @param level   level of dump
  *
  * @return string containing contents
  *
- * */
+ **/
 
 std::string BrookNonBonded::getContentsString( int level ) const {
 

platforms/brook/src/BrookNonBonded.h

@@ -391,7 +391,7 @@ class BrookNonBonded : public BrookCommon {
        *
        * */
       
-      int initializeExclusionStreamSize( int atomStreamSize, int atomStreamWidth );
+      int _initializeExclusionStreamSize( int atomStreamSize, int atomStreamWidth );
       
       /** 
        * Initialize exclusion stream dimensions and stream
@@ -402,7 +402,7 @@ class BrookNonBonded : public BrookCommon {
        *
        */
       
-      int initializeExclusionStream( const Platform& platform );
+      int _initializeExclusionStream( const Platform& platform );
             
       /** 
        * Set exclusion (4x4)
@@ -416,7 +416,7 @@ class BrookNonBonded : public BrookCommon {
        *
        */
       
-      int setExclusion( int i, int j, int exclusionStreamWidth, BrookOpenMMFloat* exclusion );
+      int _setExclusion( int i, int j, int exclusionStreamWidth, BrookOpenMMFloat* exclusion );
       
       /** 
        * Initialize exclusions
@@ -428,7 +428,7 @@ class BrookNonBonded : public BrookCommon {
        *
        */
       
-      int initializeExclusions( const std::vector<std::set<int> >& exclusionsVector, const Platform& platform );
+      int _initializeExclusions( const std::vector<std::set<int> >& exclusionsVector, const Platform& platform );
 
       /** 
        * Initialize stream dimensions
@@ -440,7 +440,7 @@ class BrookNonBonded : public BrookCommon {
        *
        */
       
-      int initializeStreamSizes( int numberOfAtoms, const Platform& platform );
+      int _initializeStreamSizes( int numberOfAtoms, const Platform& platform );
       
       /** 
        * Initialize stream dimensions and streams
@@ -451,7 +451,7 @@ class BrookNonBonded : public BrookCommon {
        *
        */
       
-      int initializeStreams( const Platform& platform );
+      int _initializeStreams( const Platform& platform );
       
       /* 
        * Setup of j-stream dimensions
@@ -465,7 +465,7 @@ class BrookNonBonded : public BrookCommon {
        *
        * */
       
-      int initializeJStreamSize( int atomStreamSize, int atomStreamWidth );
+      int _initializeJStreamSize( int atomStreamSize, int atomStreamWidth );
 
       /* 
        * Setup of outer vdw stream size
@@ -479,7 +479,7 @@ class BrookNonBonded : public BrookCommon {
        *
        * */
       
-      int initializeOuterVdwStreamSize( int atomStreamSize, int atomStreamWidth );
+      int _initializeOuterVdwStreamSize( int atomStreamSize, int atomStreamWidth );
 
       /** 
        * Set sigma & epsilon given c6 & c12 (geometric rule)
@@ -493,7 +493,7 @@ class BrookNonBonded : public BrookCommon {
        *
        */
       
-      int setSigmaEpsilon( double c6, double c12, double* sigma , double* epsilon );
+      int _setSigmaEpsilon( double c6, double c12, double* sigma , double* epsilon );
 
       /** 
        * Initialize vdw & charge
@@ -505,7 +505,19 @@ class BrookNonBonded : public BrookCommon {
        *
        */
       
-      int initializeVdwAndCharge( const std::vector<std::vector<double> >& nonbondedParameters, const Platform& platform );
+      int _initializeVdwAndCharge( const std::vector<std::vector<double> >& nonbondedParameters, const Platform& platform );
+      
+      /** 
+       * Initialize vdw & charge
+       * 
+       * @param exclusions                vector of sets containing exclusions (1 set entry for every atom)
+       * @param platform                  platform
+       *
+       * @return nonzero value if error
+       *
+       */
+      
+      int _initializeJStreamVdw( const std::vector<std::vector<double> >& nonbondedParameters, const Platform& platform );
       
       /* 
        * Setup of stream dimensions for partial force streams
@@ -517,7 +529,7 @@ class BrookNonBonded : public BrookCommon {
        *
        * */
       
-      int initializePartialForceStreamSize( int atomStreamSize, int atomStreamWidth );
+      int _initializePartialForceStreamSize( int atomStreamSize, int atomStreamWidth );
       
 };
 

platforms/brook/src/BrookStreamFactory.cpp

@@ -40,7 +40,7 @@ const std::string BrookStreamFactory::AtomPositions              = "atomPosition
 const std::string BrookStreamFactory::AtomVelocities             = "atomVelocities";
 const std::string BrookStreamFactory::AtomForces                 = "atomForces";
 
-const double DefaultDangleValue                                  = 1.0e+38;
+const double DefaultDangleValue                                  = 1.0e+08;
 /** 
  * BrookStreamFactory constructor
  * 
@@ -55,7 +55,7 @@ BrookStreamFactory::BrookStreamFactory( void ){
 
 // ---------------------------------------------------------------------------------------
 
-	_defaultDangleValue                      = 1.0e+38;
+	_defaultDangleValue                      = 1.0e+08;
 	_defaultAtomStreamWidth                  = DefaultStreamAtomWidth;
    _defaultStreamRandomNumberWidth          = DefaultStreamRandomNumberWidth;
    _defaultStreamRandomNumberSize           = DefaultStreamRandomNumberSize;

platforms/brook/src/BrookStreamInternal.cpp

@@ -281,3 +281,82 @@ std::string BrookStreamInternal::_getLine( const std::string& tab,
 
 }
 
+/* 
+ * Print contents of object to file
+ *
+ * @param log         file to print to
+ *
+ * @return DefaultReturnValue
+ *
+ * */
+
+int BrookStreamInternal::printToFile( FILE* log ){
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName      = "BrookStreamInternal::printToFile";
+
+// ---------------------------------------------------------------------------------------
+
+   if( log == NULL ){
+      log = stderr;
+   }
+   std::string contents = getContentsString();
+   (void) fprintf( log, "%s\n", contents.c_str() );
+
+   _bodyPrintToFile( log );
+
+   return DefaultReturnValue;
+
+}
+
+/* 
+ * Get contents of object
+ *
+ * @param level   level of dump
+ *
+ * @return string containing contents
+ *
+ * */
+
+const std::string BrookStreamInternal::getContentsString( int level ) const {
+
+// ---------------------------------------------------------------------------------------
+
+   // static const std::string methodName      = "BrookIntStreamInternal::getContentsString";
+
+   static const unsigned int MAX_LINE_CHARS = 256;
+   char value[MAX_LINE_CHARS];
+   //static const char* Set                   = "Set";
+   //static const char* NotSet                = "Not set";
+
+
+// ---------------------------------------------------------------------------------------
+
+   std::stringstream message;
+   std::string tab   = "   ";
+
+#ifdef WIN32
+#define LOCAL_SPRINTF(a,b,c) sprintf_s( (a), MAX_LINE_CHARS, (b), (c) );   
+#else
+#define LOCAL_SPRINTF(a,b,c) sprintf( (a), (b), (c) );   
+#endif
+
+   (void) LOCAL_SPRINTF( value, "%s", getName().c_str() );
+   message << _getLine( tab, "Name:", value );
+
+   (void) LOCAL_SPRINTF( value, "%d", getWidth() );
+   message << _getLine( tab, "Width:", value );
+
+   (void) LOCAL_SPRINTF( value, "%d", getStreamSize() );
+   message << _getLine( tab, "Stream size:", value );
+
+   (void) LOCAL_SPRINTF( value, "%d", getStreamWidth() );
+   message << _getLine( tab, "Stream width:", value );
+
+   (void) LOCAL_SPRINTF( value, "%d", getStreamHeight() );
+   message << _getLine( tab, "Stream height:", value );
+
+   return message.str();
+}
+

platforms/brook/src/BrookStreamInternal.h

@@ -50,6 +50,13 @@ class BrookStreamInternal {
         */
       enum DataType { Float, Float2, Float3, Float4, Double, Double2, Double3, Double4, Integer, Integer2, Integer3, Integer4, Unknown };
 
+      // return values
+
+      static const int DefaultReturnValue = 0;
+      static const int ErrorReturnValue   = -1; 
+
+      // ---------------------------------------------------------------------------------------
+
       /**
        * BrookStreamInternal constructor
        * 
@@ -101,7 +108,7 @@ class BrookStreamInternal {
        * the values should be packed into a single array: all the values for the first element, followed by all the values
        * for the next element, etc.
        */
-      virtual void loadFromArray(const void* array) = 0;
+      virtual void loadFromArray( const void* array ) = 0;
   
       /**
        * Copy the contents of an array into this stream.
@@ -121,14 +128,14 @@ class BrookStreamInternal {
        * the values should be packed into a single array: all the values for the first element, followed by all the values
        * for the next element, etc.
        */
-      virtual void saveToArray(void* array) = 0;
+      virtual void saveToArray( void* array ) = 0;
   
       /**
        * Set every element of this stream to the same value.
        * 
        * @param a pointer to the value.  It is assumed to be of the correct data type for this stream.
        */
-      virtual void fillWithValue(void* value) = 0;
+      virtual void fillWithValue( void* value ) = 0;
   
       /**
        * Get data
@@ -137,6 +144,22 @@ class BrookStreamInternal {
        */
       virtual void* getData( void ) = 0;
   
+      /**
+       * Get data
+       *
+       * @param readFromBoard  read data from board
+       * 
+       * @return data array
+       */
+      virtual void* getData( int readFromBoard ) = 0;
+  
+      /**
+       * Get array of appropritate size for loading data
+       *
+       * @return data array -- user's responsibility to free
+       */
+      virtual void* getDataArray( void ) = 0;
+  
       /** 
        * Get type string
        *
@@ -200,6 +223,17 @@ class BrookStreamInternal {
 
       const std::string getContentsString( int level = 0 ) const;
 
+      /* 
+       * Print to file
+       *
+       * @param log  log file
+       *
+       * @return  DefaultReturnValue
+       *
+       * */
+
+      int printToFile( FILE* log );
+
    protected:
 
       std::string _name;
@@ -229,6 +263,16 @@ class BrookStreamInternal {
       std::string _getLine( const std::string& tab, const std::string& description, 
                             const std::string& value ) const;
       
+      /* 
+       * Print array to file
+       *
+       * @param log  log file
+       *
+       * @return  DefaultReturnValue
+       *
+       * */
+
+      virtual int _bodyPrintToFile( FILE* log ) = 0;
 };
 
 } // namespace OpenMM

platforms/brook/src/gpu/invmap.cpp

@@ -57,6 +57,9 @@ gpuCalcInvMap(
 	//This will hold the number of imaps actually used
 	*nimaps = -1;
 
+   printf( "gpuCalcInvMap: posflag=%d niatoms=%d nints=%d natoms=%d nmaps=%d\n",
+           posflag, niatoms, nints, natoms, nmaps );
+
 	//Now note down the positions where each atom occurs
 	for ( i = 0; i < nints; i++ ) {
 		//This is our atom
@@ -94,6 +97,9 @@ gpuCalcInvMap(
 		}
 		
 		counts[atom]++;
+
+printf( "Atom %d count=%d max %d mapcomp=%d val=%d mapnum=%d\n", atom,  counts[atom],
+        nmaps * 4, mapcomp, i, mapnum );
 	}
 
 	(*nimaps)++;

platforms/brook/src/gpu/kforce.h

@@ -65,6 +65,7 @@ void  knbforce_CDLJ4(
 		const float  epsfac,
 		const float4 params,
 		::brook::stream posq,
+		::brook::stream charge,
 		::brook::stream isigeps,
 		::brook::stream sigma,
 		::brook::stream epsilon,
@@ -115,6 +116,7 @@ void  knbforce_CDLJ4NoEx(
 		const float  epsfac,
 		const float4 params,
 		::brook::stream posq,
+		::brook::stream charge,
 		::brook::stream isigeps,
 		::brook::stream sigma,
 		::brook::stream epsilon,
@@ -286,6 +288,7 @@ void  kbonded_CDLJ (const float  epsfac,
 		const float  xstrwidth,
 		const float4 params,
 		::brook::stream posq,
+		::brook::stream charge,
 		::brook::stream atoms,
 		::brook::stream parm0,
 		::brook::stream parm1,

platforms/brook/src/gpu/kforce_CDLJ.br

@@ -12,7 +12,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 ){
+kernel float scalar_force_single_CDLJ( float qq, float epsfac, float sig, float eps, float r2, float4 params ) {
  float invr, invrsig2, invrsig6;
  float f;
 
@@ -409,6 +409,263 @@ kernel void knbforce_CDLJ2(
 
 }
 
+kernel void knbforce_CDLJ4Ex(
+  float natoms,
+  float nAtomsCeiling,
+  float dupfac,
+  float strheight,
+  float strwidth,
+  float jstrwidth,
+  float fstrwidth,
+  float epsfac,
+  float4 params,
+  float3 posq[][],
+  float  charge[][],
+  float4 isigeps[][],
+  float4 sigma[][],
+  float4 epsilon[][],
+  float excl[][],
+  out float3 outforce1<>,
+  out float3 outforce2<>,
+  out float3 outforce3<>,
+  out float3 outforce4<> ) {
+
+ float4 jsig, jeps;
+ float3 ipos1, ipos2, ipos3, ipos4;
+ float4 jpos;
+ float3 jpos1, jpos2, jpos3, jpos4;
+ float qi1, qi2, qi3, qi4;
+
+ float isig1, isig2, isig3, isig4;
+ float ieps1, ieps2, ieps3, ieps4;
+
+ float4 qj, qq;
+ float3 d1, d2, d3, d4;
+ float2 exclind;
+ float4 sig, eps;
+ float  exclusions;
+ float4 r2;
+ float3 pad;
+ float4 exclmask;
+ float4 exclconst;
+ float2 iatom;
+ float2 jstrindex;
+ float a_iatom;
+ float linind;
+
+ float breakflag;
+ float2 jatom;
+ float4 fs;
+ float jend, jstart;
+ float which_rep;
+
+ exclconst  = float4( 2.0f, 3.0f, 5.0f, 7.0f );
+
+ outforce1  = float4( 0.0f, 0.0f, 0.0f, 0.0f );
+ outforce2  = float4( 0.0f, 0.0f, 0.0f, 0.0f );
+ outforce3  = float4( 0.0f, 0.0f, 0.0f, 0.0f );
+ outforce4  = float4( 0.0f, 0.0f, 0.0f, 0.0f );
+
+ a_iatom    = 4.0f*((indexof outforce1).x + (indexof outforce1).y * fstrwidth);
+
+ // linind     = fmod( a_iatom, natoms );
+ linind     = fmod( a_iatom, nAtomsCeiling );
+
+ iatom.x    = fmod( linind, strwidth );
+ iatom.y    = round( (linind - fmod(linind, strwidth))/strwidth );
+
+/* --------------------------------------------------------------------- */
+
+ ipos1      = posq[ iatom ];
+ qi1        = charge[ iatom ];
+
+ jstrindex  = isigeps[ iatom ];
+ isig1      = jstrindex.x;
+ ieps1      = jstrindex.y;
+
+/* --------------------------------------------------------------------- */
+
+ iatom.x   += 1;
+ ipos2      = posq[ iatom ];
+ qi2        = charge[ iatom ];
+
+ jstrindex  = isigeps[ iatom ];
+ isig2      = jstrindex.x;
+ ieps2      = jstrindex.y;
+
+/* --------------------------------------------------------------------- */
+
+ iatom.x   += 1;
+ ipos3      = posq[ iatom ];
+ qi3        = charge[ iatom ];
+
+ jstrindex  = isigeps[ iatom ];
+ isig3      = jstrindex.x;
+ ieps3      = jstrindex.y;
+
+/* --------------------------------------------------------------------- */
+
+ iatom.x   += 1;
+ ipos4      = posq[ iatom ];
+ qi4        = charge[ iatom ];
+
+ jstrindex  = isigeps[ iatom ];
+ isig4      = jstrindex.x;
+ ieps4      = jstrindex.y;
+
+/* --------------------------------------------------------------------- */
+
+ breakflag  = 1.0f;
+ jatom.y    = 0.0f;
+
+ // which_rep  = round( (a_iatom - fmod(a_iatom, natoms))/natoms );
+ which_rep  = round( (a_iatom - fmod(a_iatom, nAtomsCeiling ))/nAtomsCeiling );
+
+ jend       = 1 + round( (natoms - fmod(natoms, (dupfac * strwidth )))/(dupfac * strwidth ) );
+ jstart     = which_rep * jend;
+ if ( which_rep > dupfac - 1 - 0.5f ) {
+  jend = 1e6f;
+ }
+ jend      += jstart;
+
+ // exclind.x  = linind;
+ // exclind.y  = jstart * strwidth/4.0f;
+
+ jatom.y    = jstart;
+ while ( jatom.y < jend && breakflag > 0.0f ) {
+  jatom.x = 0;
+  while ( jatom.x < strwidth && breakflag > 0.0f ) {
+
+   //exclusions  = excl[ exclind ];
+   //exclmask    = fmod( exclusions, exclconst );
+
+   linind      = (jatom.x + jatom.y*strwidth)/4.0f;
+
+   jstrindex.y = round( (linind - fmod(linind, jstrwidth))/jstrwidth );
+   jstrindex.x = linind - jstrindex.y * jstrwidth;
+
+   /* --------------------------------------------------------------------- */
+
+   jsig        = sigma[ jstrindex ];
+   jeps        = epsilon[ jstrindex ];
+
+   jpos1       = posq[ jatom ];
+   qj.x        = charge[ jatom ];
+   jatom.x    += 1.0f;
+
+   jpos2       = posq[ jatom ];
+   qj.y        = charge[ jatom ];
+   jatom.x    += 1.0f;
+
+   jpos3       = posq[ jatom ];
+   qj.z        = charge[ jatom ];
+   jatom.x    += 1.0f;
+
+   jpos4       = posq[ jatom ];
+   qj.w        = charge[ jatom ];
+   jatom.x    += 1.0f;
+
+   /* --------------------------------------------------------------------- */
+
+   sig         = isig1 + jsig;
+   eps         = ieps1 * jeps;
+   qq          = qi1*qj;
+
+   d1          = ipos1 - jpos1;
+   d2          = ipos1 - jpos2;
+   d3          = ipos1 - jpos3;
+   d4          = ipos1 - jpos4;
+
+   r2          = get_r2_CDLJ( d1, d2, d3, d4 );
+   fs          = scalar_force_CDLJ( qq, epsfac, sig, eps, r2, params );
+
+   outforce1  += fs.x * d1;
+   outforce1  += fs.y * d2;
+   outforce1  += fs.z * d3;
+   outforce1  += fs.w * d4;
+
+   /* --------------------------------------------------------------------- */
+
+   //exclind.x  += 1.0f;
+   //exclusions  = excl[ exclind ];
+   //exclmask    = fmod( exclusions, exclconst );
+
+   sig         = isig2 + jsig;
+   eps         = ieps2 * jeps;
+   qq          = qi2 * qj;
+
+   d1          = ipos2 - jpos1;
+   d2          = ipos2 - jpos2;
+   d3          = ipos2 - jpos3;
+   d4          = ipos2 - jpos4;
+
+   r2          = get_r2_CDLJ( d1, d2, d3, d4 );
+   fs          = scalar_force_CDLJ( qq, epsfac, sig, eps, r2, params );
+
+   outforce2  += fs.x * d1;
+   outforce2  += fs.y * d2;
+   outforce2  += fs.z * d3;
+   outforce2  += fs.w * d4;
+
+   /* --------------------------------------------------------------------- */
+
+   //exclind.x  += 1.0f;
+   //exclusions  = excl[ exclind ];
+   //exclmask    = fmod( exclusions, exclconst );
+
+   sig         = isig3 + jsig;
+   eps         = ieps3 * jeps;
+   qq          = qi3   * qj;
+
+   d1          = ipos3 - jpos1;
+   d2          = ipos3 - jpos2;
+   d3          = ipos3 - jpos3;
+   d4          = ipos3 - jpos4;
+
+   r2          = get_r2_CDLJ( d1, d2, d3, d4 );
+   fs          = scalar_force_CDLJ( qq, epsfac, sig, eps, r2, params );
+
+   outforce3  += fs.x * d1;
+   outforce3  += fs.y * d2;
+   outforce3  += fs.z * d3;
+   outforce3  += fs.w * d4;
+
+   /* --------------------------------------------------------------------- */
+
+   //exclind.x  += 1.0f;
+   //exclusions  = excl[ exclind ];
+   //exclmask    = fmod( exclusions, exclconst );
+
+   sig         = isig4 + jsig;
+   eps         = ieps4 * jeps;
+   qq          = qi4   * qj;
+
+   d1          = ipos4 - jpos1;
+   d2          = ipos4 - jpos2;
+   d3          = ipos4 - jpos3;
+   d4          = ipos4 - jpos4;
+
+   r2          = get_r2_CDLJ( d1, d2, d3, d4 );
+   fs          = scalar_force_CDLJ( qq, epsfac, sig, eps, r2, params );
+
+   outforce4  += fs.x * d1;
+   outforce4  += fs.y * d2;
+   outforce4  += fs.z * d3;
+   outforce4  += fs.w * d4;
+
+   /* --------------------------------------------------------------------- */
+
+   //exclind.x  -= 3.0f;
+   //exclind.y  += 1.0f;
+
+   if ( natoms - (jatom.y*strwidth + jatom.x) < 0.9f )
+    breakflag = -1.0f;
+  }
+
+  jatom.y += 1.0f;
+ }
+
+}
 
 kernel void knbforce_CDLJ4(
   float natoms,
@@ -420,7 +677,8 @@ kernel void knbforce_CDLJ4(
   float fstrwidth,
   float epsfac,
   float4 params,
-  float4 posq[][],
+  float3 posq[][],
+  float  charge[][],
   float4 isigeps[][],
   float4 sigma[][],
   float4 epsilon[][],
@@ -458,6 +716,9 @@ kernel void knbforce_CDLJ4(
  float4 fs;
  float jend, jstart;
  float which_rep;
+ float maskMultiplier;
+
+ maskMultiplier = 1.0e+06f;
 
  exclconst  = float4( 2.0f, 3.0f, 5.0f, 7.0f );
 
@@ -476,9 +737,8 @@ kernel void knbforce_CDLJ4(
 
 /* --------------------------------------------------------------------- */
 
- jpos       = posq[ iatom ];
- ipos1      = jpos.xyz;
- qi1        = jpos.w;
+ ipos1      = posq[ iatom ];
+ qi1        = charge[ iatom ];
 
  jstrindex  = isigeps[ iatom ];
  isig1      = jstrindex.x;
@@ -486,10 +746,9 @@ kernel void knbforce_CDLJ4(
 
 /* --------------------------------------------------------------------- */
 
- iatom.x   += 1.0f;
- jpos       = posq[ iatom ];
- ipos2      = jpos.xyz;
- qi2        = jpos.w;
+ iatom.x   += 1;
+ ipos2      = posq[ iatom ];
+ qi2        = charge[ iatom ];
 
  jstrindex  = isigeps[ iatom ];
  isig2      = jstrindex.x;
@@ -497,10 +756,9 @@ kernel void knbforce_CDLJ4(
 
 /* --------------------------------------------------------------------- */
 
- iatom.x   += 1.0f;
- jpos       = posq[ iatom ];
- ipos3      = jpos.xyz;
- qi3        = jpos.w;
+ iatom.x   += 1;
+ ipos3      = posq[ iatom ];
+ qi3        = charge[ iatom ];
 
  jstrindex  = isigeps[ iatom ];
  isig3      = jstrindex.x;
@@ -508,10 +766,10 @@ kernel void knbforce_CDLJ4(
 
 /* --------------------------------------------------------------------- */
 
- iatom.x   += 1.0f;
- jpos       = posq[ iatom ];
- ipos4      = jpos.xyz;
- qi4        = jpos.w;
+ iatom.x   += 1;
+ ipos4      = posq[ iatom ];
+ qi4        = charge[ iatom ];
+
 
  jstrindex  = isigeps[ iatom ];
  isig4      = jstrindex.x;
@@ -556,29 +814,21 @@ kernel void knbforce_CDLJ4(
    jsig        = sigma[ jstrindex ];
    jeps        = epsilon[ jstrindex ];
 
-   jpos        = posq[ jatom ];
-   jpos1       = jpos.xyz;
-   qj.x        = jpos.w;
+   jpos1       = posq[ jatom ];
+   qj.x        = charge[ jatom ];
    jatom.x     = floor( jatom.x + 1.2f );
-   //jatom.x    += 1.0f;
 
-   jpos        = posq[ jatom ];
-   jpos2       = jpos.xyz;
-   qj.y        = jpos.w;
+   jpos2       = posq[ jatom ];
+   qj.y        = charge[ jatom ];
    jatom.x     = floor( jatom.x + 1.2f );
-   //jatom.x    += 1.0f;
 
-   jpos        = posq[ jatom ];
-   jpos3       = jpos.xyz;
-   qj.z        = jpos.w;
+   jpos3       = posq[ jatom ];
+   qj.z        = charge[ jatom ];
    jatom.x     = floor( jatom.x + 1.2f );
-   //jatom.x    += 1.0f;
 
-   jpos        = posq[ jatom ];
-   jpos4       = jpos.xyz;
-   qj.w        = jpos.w;
+   jpos4       = posq[ jatom ];
+   qj.w        = charge[ jatom ];
    jatom.x     = floor( jatom.x + 1.2f );
-   //jatom.x    += 1.0f;
 
    /* --------------------------------------------------------------------- */
 
@@ -591,7 +841,7 @@ kernel void knbforce_CDLJ4(
    d3          = ipos1 - jpos3;
    d4          = ipos1 - jpos4;
 
-   r2          = get_r2_CDLJ( d1, d2, d3, d4 ) + exclmask*10000.0f;
+   r2          = get_r2_CDLJ( d1, d2, d3, d4 ) + (exclmask*maskMultiplier);
    fs          = scalar_force_CDLJ( qq, epsfac, sig, eps, r2, params );
 
    outforce1  += fs.x * d1;
@@ -615,7 +865,7 @@ kernel void knbforce_CDLJ4(
    d3          = ipos2 - jpos3;
    d4          = ipos2 - jpos4;
 
-   r2          = get_r2_CDLJ( d1, d2, d3, d4 ) + exclmask * 10000.0f;
+   r2          = get_r2_CDLJ( d1, d2, d3, d4 ) + (exclmask*maskMultiplier);
    fs          = scalar_force_CDLJ( qq, epsfac, sig, eps, r2, params );
 
    outforce2  += fs.x * d1;
@@ -639,7 +889,7 @@ kernel void knbforce_CDLJ4(
    d3          = ipos3 - jpos3;
    d4          = ipos3 - jpos4;
 
-   r2          = get_r2_CDLJ( d1, d2, d3, d4 ) + exclmask * 10000.0f;
+   r2          = get_r2_CDLJ( d1, d2, d3, d4 ) + (exclmask*maskMultiplier);
    fs          = scalar_force_CDLJ( qq, epsfac, sig, eps, r2, params );
 
    outforce3  += fs.x * d1;
@@ -663,7 +913,7 @@ kernel void knbforce_CDLJ4(
    d3          = ipos4 - jpos3;
    d4          = ipos4 - jpos4;
 
-   r2          = get_r2_CDLJ( d1, d2, d3, d4 ) + exclmask*10000.0f;
+   r2          = get_r2_CDLJ( d1, d2, d3, d4 ) + (exclmask*maskMultiplier);
    fs          = scalar_force_CDLJ( qq, epsfac, sig, eps, r2, params );
 
    outforce4  += fs.x * d1;
@@ -688,6 +938,210 @@ kernel void knbforce_CDLJ4(
 
 }
 
+kernel void knbforce_CDLJ_1(
+  float natoms,
+  float nAtomsCeiling,
+  float strwidth,
+  float epsfac,
+  float4 posq[][],
+  float2 isigeps[][],
+  float excl[][],
+  out float3 outforce<> ){
+
+ float  jsig, jeps;
+ float3 ipos;
+ float3 jpos;
+
+ float qi, qj;
+ float3 d1;
+ float sig, eps;
+ float isig, ieps;
+
+ float2 exclind;
+ float  exclusion;
+ float r2, invr, invrsig2, invrsig6;
+ float4 exclmask;
+ float2 iatom;
+ float a_iatom;
+ float linind;
+
+ float2 jatom;
+ float  fs;
+
+ outforce   = float3( 0.0f, 0.0f, 0.0f );
+
+ iatom.x    = (indexof outforce).x;
+ iatom.y    = (indexof outforce).y;
+ a_iatom    = iatom.x + iatom.y*strwidth;
+
+/* --------------------------------------------------------------------- */
+
+ ipos       = posq[ iatom ].xyz;
+ qi         = epsfac*posq[ iatom ].w; 
+
+ isig       = isigeps[ iatom ].x; 
+ ieps       = isigeps[ iatom ].y; 
+
+/* --------------------------------------------------------------------- */
+
+ jatom      = float2( 0.0f, 0.0f );
+ linind     = 0.0f;
+
+ exclind.x  = 0.0f;
+ exclind.y  = a_iatom;
+
+ while ( linind < natoms ){
+   
+    exclusion   = excl[ exclind ];
+ 
+    /* --------------------------------------------------------------------- */
+ 
+    jsig       = isigeps[ jatom ].x;
+    jeps       = isigeps[ jatom ].y;
+ 
+    jpos       = posq[ jatom ].xyz;
+    qj         = posq[ jatom ].w;
+ 
+    /* --------------------------------------------------------------------- */
+ 
+    sig         = isig + jsig;
+    eps         = ieps * jeps;
+ 
+    d1          = ipos - jpos;
+ 
+    r2          = dot( d1, d1 ) + exclusion*10000.0f;
+    invr        = rsqrt( r2 );
+    
+    invrsig2    = invr*sig;
+    invrsig2    = invrsig2*invrsig2;
+    invrsig6    = invrsig2*invrsig2*invrsig2;
+    fs          = eps* ( 12.0f*invrsig6 - 6.0f ) * invrsig6;
+ 
+    fs         += qi*qj*invr;
+ 
+    fs         *= invr*invr;
+ 
+    outforce   += fs*d1;
+ 
+    /* --------------------------------------------------------------------- */
+ 
+    exclind.x  = floor(  exclind.x  + 1.05f );
+ 
+    jatom.x   += 1.0f;
+    if( abs( jatom.x - strwidth ) < 0.1f ){
+       jatom.x  = 0.0f;
+       jatom.y += 1.0f;
+    }
+ 
+    linind += 1.0f;
+ }
+
+}
+
+
+kernel void knbforce_CDLJ_1z(
+  float natoms,
+  float nAtomsCeiling,
+  float strwidth,
+  float epsfac,
+  float4 posq[][],
+  float2 isigeps[][],
+  float excl[][],
+  out float3 outforce<> ){
+
+ float  jsig, jeps;
+ float3 ipos;
+ float3 jpos;
+
+ float qi, qj;
+ float3 d1;
+ float sig, eps;
+ float isig, ieps;
+
+ float2 exclind;
+ float  exclusion;
+ float r2, invr, invrsig2, invrsig6;
+ float4 exclmask;
+ float2 iatom;
+ float a_iatom;
+ float linind;
+
+ float2 jatom;
+ float  fs;
+
+ outforce   = float3( 0.0f, 0.0f, 0.0f );
+
+ iatom.x    = (indexof outforce).x;
+ iatom.y    = (indexof outforce).y;
+ a_iatom    = iatom.x + iatom.y*strwidth;
+
+/* --------------------------------------------------------------------- */
+
+ ipos       = posq[ iatom ].xyz;
+ qi         = epsfac*posq[ iatom ].w; 
+
+ isig       = isigeps[ iatom ].x; 
+ ieps       = isigeps[ iatom ].y; 
+
+/* --------------------------------------------------------------------- */
+
+ jatom      = float2( 0.0f, 0.0f );
+ linind     = 0.0f;
+
+ exclind.x  = 0.0f;
+ exclind.y  = a_iatom;
+
+//outforce.y = exclind.y;
+//outforce.z = a_iatom;
+
+   while ( linind < natoms ){
+   
+      exclusion   = excl[ exclind ];
+   
+      /* --------------------------------------------------------------------- */
+ 
+      jsig       = isigeps[ jatom ].x;
+      jeps       = isigeps[ jatom ].y;
+   
+      jpos       = posq[ jatom ].xyz;
+      qj         = posq[ jatom ].w;
+ 
+      /* --------------------------------------------------------------------- */
+   
+      eps         = ieps * jeps;
+   
+      d1          = ipos - jpos;
+   
+      r2          = dot( d1, d1 ) + exclusion*10000.0f;
+      invr        = rsqrt( r2 );
+      
+      invrsig2    = invr*sig;
+      invrsig2    = invrsig2*invrsig2;
+      invrsig6    = invrsig2*invrsig2*invrsig2;
+      fs          = eps* ( 12.0f*invrsig6 - 6.0f ) * invrsig6;
+   
+      fs         += qi*qj*invr;
+   
+      fs         *= invr*invr;
+   
+      outforce   += fs*d1;
+      //outforce.x += exclusion;
+   
+      /* --------------------------------------------------------------------- */
+   
+      exclind.x  = floor(  exclind.x  + 1.05f );
+   
+      jatom.x   += 1.0f;
+      if( abs( jatom.x - strwidth ) < 0.1f ){
+         jatom.x  = 0.0f;
+         jatom.y += 1.0f;
+      }
+   
+      linind += 1.0f;
+   }
+
+}
+
 kernel void knbforce_CDLJ4Debug(
   float natoms,
   float nAtomsCeiling,
@@ -741,8 +1195,10 @@ kernel void knbforce_CDLJ4Debug(
  float4 zero4;
  float4 frac4;
 float sum;
+float bigValue;
 
  exclconst  = float4( 2.0f, 3.0f, 5.0f, 7.0f );
+ bigValue   = 1000000.0f;
 
  outforce1  = float3( 0.0f, 0.0f, 0.0f );
  outforce2  = float3( 0.0f, 0.0f, 0.0f );
@@ -828,7 +1284,7 @@ float sum;
 /*
 outforce1   = float3( linind, exclind.y, jstart );
 outforce2   = float3( linind, exclind.y, jend );
-outforce3   = float3( linind, exclind.y, strwidth );
+outforce3   = float3( linind, exclind.y, which_rep );
 outforce4   = float3( linind, exclind.y, exclind.y );
 */
 
@@ -883,9 +1339,13 @@ outforce4   = float3( linind, exclind.y, exclind.y );
    exclusions  = excl[ exclind ];
    exclmask    = fmod( exclusions, exclconst );
 
-   r2          = get_r2_CDLJ( d1, d2, d3, d4 ) + exclmask*10000.0f;
+   r2          = get_r2_CDLJ( d1, d2, d3, d4 );
+   r2         += exclmask*10000.0f;
    fs          = scalar_force_CDLJ( qq, epsfac, sig, eps, r2, params );
 
+//r2 = 1.0f/r2;
+//outforce1.x += r2.x + r2.y + r2.z + r2.w;
+//outforce1.y += 4.0f;
    //    outforce1  += fs.x * d1;
    //    outforce1  += fs.y * d2;
    //    outforce1  += fs.z * d3;
@@ -895,19 +1355,14 @@ outforce4   = float3( linind, exclind.y, exclind.y );
 
    exclmask    = exclmask > frac4 ? one4 : zero4;
    sum         = floor( exclmask.x + exclmask.y + exclmask.z + exclmask.w + 0.25f );
-/*
+
+   //if( exclind.y < 5.9f && exclind.y > 4.1f ){
    if( sum > 0.25f ){
        outforce1.x  += 1.0f;
        outforce1.y  += exclusions;
+       //outforce1.y  += linind;
        outforce1.z  += sum;
-        // outforce1.y  += exclind.x + strwidth*exclind.y;
-        //outforce1.y  += 4.0;
-    }
-*/
-   if( exclind.y < 5.9f && exclind.y > 4.1f ){
-       outforce1.y  += exclind.x + 576.0f*exclind.y;
-       outforce1.z  += exclusions;
-       outforce1.x  += exclind.y;
+       //outforce1.z  += exclind.x + exclind.y*576;
     }    
 
    // --------------------------------------------------------------------- 
@@ -925,22 +1380,39 @@ outforce4   = float3( linind, exclind.y, exclind.y );
    d3          = ipos2 - jpos3;
    d4          = ipos2 - jpos4;
 
-   r2          = get_r2_CDLJ( d1, d2, d3, d4 ) + exclmask*10000.0f;
+   r2          = get_r2_CDLJ( d1, d2, d3, d4 );
+   r2         += exclmask*10000.0f;
    fs          = scalar_force_CDLJ( qq, epsfac, sig, eps, r2, params );
 
-   // outforce2  += fs.x * d1;
-   // outforce2  += fs.y * d2;
-   // outforce2  += fs.z * d3;
-   // outforce2  += fs.w * d4;
+//r2 =1.0f/r2;
+//outforce2.x += r2.x + r2.y + r2.z + r2.w;
+//outforce2.y += 4.0f;
+
+//
+//   outforce2  += fs.x * d1;
+//   outforce2  += fs.y * d2;
+//   outforce2  += fs.z * d3;
+//   outforce2  += fs.w * d4;
+//
    // outforce2  += float3( 4.0f, 4.0f, 4.0f);
     exclmask    = exclmask > frac4 ? one4 : zero4;
     sum         = floor( exclmask.x + exclmask.y + exclmask.z + exclmask.w + 0.25f );
     if( sum > 0.25f ){
         outforce2.x  += 1.0f;
-        outforce2.y  += exclusions;
-        outforce2.z  += sum;
+        outforce2.y  += linind;
+       // outforce2.z  += sum;
+       outforce2.z  += exclind.x + exclind.y*576;
     }
 
+//
+//    sum         = abs( fs.x ) + abs( fs.y ) + abs( fs.z ) + abs( fs.w );
+//    if( sum > bigValue ){
+//        outforce2.x   = linind;
+//        outforce2.y   = linind1;
+//        outforce2.z   = sum;
+//    }
+//
+//
    // --------------------------------------------------------------------- 
 
    exclind.x   = floor( exclind.x + 1.2f );
@@ -956,24 +1428,38 @@ outforce4   = float3( linind, exclind.y, exclind.y );
    d3          = ipos3 - jpos3;
    d4          = ipos3 - jpos4;
 
-   r2          = get_r2_CDLJ( d1, d2, d3, d4 ) + exclmask * 10000.0f;
+   r2          = get_r2_CDLJ( d1, d2, d3, d4 );
+   r2         += exclmask * 10000.0f;
    fs          = scalar_force_CDLJ( qq, epsfac, sig, eps, r2, params );
 
+//r2 =1.0f/r2;
+//outforce3.x += r2.x + r2.y + r2.z + r2.w;
+//outforce3.y += 4.0f;
+
    // outforce3  += fs.x * d1;
    // outforce3  += fs.y * d2;
    // outforce3  += fs.z * d3;
    // outforce3  += fs.w * d4;
    // outforce3  += float3( 4.0f, 4.0f, 4.0f);
    // outforce3  += 1.0f/r2;
+
     exclmask    = exclmask > frac4 ? one4 : zero4;
     sum         = floor( exclmask.x + exclmask.y + exclmask.z + exclmask.w + 0.25f );
     if( sum > 0.25f ){
         outforce3.x  += 1.0f;
-        outforce3.y  += exclusions;
-        outforce3.z  += sum;
-        // outforce3.y  += exclind.x + strwidth*exclind.y;
+        outforce3.y  += linind;
+        //outforce3.z  += sum;
+        outforce3.z  += exclind.x + strwidth*exclind.y;
         //outforce3.y  += 4.0;
     }
+//
+//    sum         = abs( fs.x ) + abs( fs.y ) + abs( fs.z ) + abs( fs.w );
+//    if( sum > bigValue ){
+//        outforce3.x   = linind;
+//        outforce3.y   = linind1;
+//        outforce3.z   = sum;
+//    }
+//
 
    // --------------------------------------------------------------------- 
 
@@ -990,25 +1476,38 @@ outforce4   = float3( linind, exclind.y, exclind.y );
    d3          = ipos4 - jpos3;
    d4          = ipos4 - jpos4;
 
-   r2          = get_r2_CDLJ( d1, d2, d3, d4 ) + exclmask * 10000.0f;
+   r2          = get_r2_CDLJ( d1, d2, d3, d4 );
+   r2          = exclmask * 10000.0f;
    fs          = scalar_force_CDLJ( qq, epsfac, sig, eps, r2, params );
 
+//r2 =1.0f/r2;
+//outforce4.x += r2.x + r2.y + r2.z + r2.w;
+//outforce4.y += 4.0f;
+
    // outforce4  += fs.x * d1;
    // outforce4  += fs.y * d2;
    // outforce4  += fs.z * d3;
    // outforce4  += fs.w * d4;
    // outforce4  += float3( 4.0f, 4.0f, 4.0f);
    // outforce4  += 1.0f/r2;
-
     exclmask    = exclmask > frac4 ? one4 : zero4;
     sum         = floor( exclmask.x + exclmask.y + exclmask.z + exclmask.w + 0.25f );
     if( sum > 0.25f ){
         outforce4.x  += 1.0f;
-        outforce4.y  += exclusions;
-        outforce4.z  += sum;
-        // outforce4.y  += exclind.x + strwidth*exclind.y;
+        //outforce4.y  += exclusions;
+        outforce4.y  += linind;
+        //outforce4.z  += sum;
+        outforce4.z  += exclind.x + strwidth*exclind.y;
         //outforce4.y  += 4.0;
     }
+//
+//    sum         = abs( fs.x ) + abs( fs.y ) + abs( fs.z ) + abs( fs.w );
+//    if( sum > bigValue ){
+//        outforce4.x   = linind;
+//        outforce4.y   = linind1;
+//        outforce4.z   = sum;
+//    }
+//
 
    // --------------------------------------------------------------------- 
 
@@ -1024,9 +1523,48 @@ outforce4   = float3( linind, exclind.y, exclind.y );
 
 }
 
+kernel void kforce14_CDLJ(
+  float natoms,
+  float strwidth,
+  float epsfac,
+  float4 params,
+  float2 atoms<>,
+  float4 posq[][],
+  float2 sigeps<>,
+  out float3 fi<>,
+  out float3 fj<>
+  )
+{
+ float2 ai, aj;
+ float3 d1;
+ float r2;
+ float qq, fs;
+
+ //ai.y = floor( atoms.x / strwidth );
+ ai.y = round( (atoms.x - fmod( atoms.x, strwidth ))/strwidth );
+ ai.x = atoms.x - ai.y * strwidth;
+
+ //aj.y = floor( atoms.y / strwidth );
+ aj.y = round( (atoms.y - fmod( atoms.y, strwidth ))/strwidth );
+ aj.x = atoms.y - aj.y * strwidth;
+
+ d1 = posq[ai].xyz - posq[aj].xyz;
+ qq = posq[ai].w * posq[aj].w;
+
+ r2 = dot( d1, d1 );
+
+ fs = scalar_force_single_CDLJ( qq, epsfac, sigeps.x, sigeps.y, r2, params );
+
+ fi = fs.x * d1;
+ fj = -fi;
+}
+
 kernel void kbonded_CDLJ (
+  float epsfac,
   float xstrwidth,
-  float4 posq[][],
+  float4 nbparams,
+  float3 posq[][],
+  float  charge[][],
   float4 atoms<>,
   float4 parm0<>,
   float4 parm1<>,
@@ -1061,12 +1599,24 @@ kernel void kbonded_CDLJ (
  float torsion_numerator;
  float theta1, costheta1, invsintheta1;
  float theta2, costheta2, invsintheta2;
- float4 posqi, posqj, posqk, posql;
+ float3 posqi, posqj, posqk, posql;
+ float  chargei, chargel;
  float sig, eps;
 
+
+
+
+
+
    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 ) {
 
@@ -1091,11 +1641,13 @@ kernel void kbonded_CDLJ (
       al.y = round( (atoms.w - fmod( atoms.w, xstrwidth ))/xstrwidth );
   al.x = atoms.w - al.y * xstrwidth;
   posql = posq[ al ];
- }
- else
+  chargel = charge[ al ];
+ } else {
   posql = float4( 1.0f, 1.0f, 1.0f, 0.0f );
+  chargel = 0.0f;
+}
 
- qq = posqi.w * posql.w;
+ qq = chargei*chargel;
  rij = posqi.xyz - posqj.xyz;
  rkj = posqk.xyz - posqj.xyz;
  rkl = posqk.xyz - posql.xyz;
@@ -1218,7 +1770,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 );
+  fs = scalar_force_single_CDLJ( qq, epsfac, parm4.z, parm4.w, r2, nbparams );
   fi_pair = fs * ril;
   fi += fi_pair;
   fl -= fi_pair;
@@ -1440,7 +1992,7 @@ kernel void kbonded_CDLJDebug (
 /*
  if ( parm4.z > -0.5f ) {
   r2 = dot( ril, ril );
-  fs = scalar_force_single_CDLJ( qq, epsfac, parm4.z, parm4.w, r2 );
+  fs = scalar_force_single_CDLJ( qq, epsfac, parm4.z, parm4.w, r2, nbparams );
   fi_pair = fs * ril;
 //fi_pair = float3( posqi.w + posql.w, posqi.w, posql.w );
 //fi_pair = float3( posqi.w + posql.w, posqi.w, posql.w );