Mods

platforms/brook/src/BrookBonded.cpp

@@ -62,10 +62,6 @@ BrookBonded::BrookBonded( ){
    _setupCompleted            = 0;
    _numberOfParticles         = 0;
 
-   //_ljScale                   = (BrookOpenMMFloat) 0.83333333;
-   _ljScale                   = 1.0;
-   //_coulombFactor             = 332.0;
-
    _coulombFactor             = (BrookOpenMMFloat) 138.935485;
 
    _particleIndicesStream     = NULL;
@@ -195,17 +191,6 @@ int BrookBonded::getInverseMapStreamWidth( void ) const {
    return _inverseMapStreamWidth;
 }
 
-/** 
- * Get LJ 14 scaling parameter
- * 
- * @return LJ 14 scaling parameter
- *
- */
-
-BrookOpenMMFloat BrookBonded::getLJ_14Scale( void ) const {
-   return _ljScale;
-}
-
 /** 
  * Get Coulomb factor
  * 
@@ -981,7 +966,6 @@ int saveIbond = ibonded;
  * @param bonded14Indices           each element contains the indices of two particles whose nonbonded interactions should be reduced since
  *                                  they form a bonded 1-4 pair
  * @param nonbondedParameters       the nonbonded force parameters (charge, sigma, epsilon) for each particle
- * @param lj14Scale                 the factor by which van der Waals interactions should be reduced for bonded 1-4 pairs
  * @param log                       log reference
  *
  * @return nonzero value if error
@@ -991,8 +975,7 @@ int saveIbond = ibonded;
 int BrookBonded::addPairs( int *nbondeds, int *particles, BrookOpenMMFloat* params[],
                            BrookOpenMMFloat* charges,
                            const std::vector<std::vector<int> >& bonded14Indices,
-                           const std::vector<std::vector<double> >& nonbondedParameters,
-                           double lj14Scale, double coulombScale ){
+                           const std::vector<std::vector<double> >& nonbondedParameters ){
 
 // ---------------------------------------------------------------------------------------
 
@@ -1192,8 +1175,6 @@ int BrookBonded::loadInvMaps( int nbondeds, int nparticles, int *particles, int
  * @param rbTorsionParameters          vector of vector of rb torsion                 bond parameters -- one entry each bond (5 parameters)
  * @param bonded14Indices              vector of vector of Lennard-Jones 14           particle indices    -- one entry each bond (2 particles     )
  * @param nonbondedParameters          vector of vector of Lennard-Jones 14           parameters      -- one entry each bond (3 parameters)
- * @param lj14Scale                    scaling factor for 1-4 ixns
- * @param coulombScale                 Coulomb scaling factor for 1-4 ixns
  * @param platform                     Brook platform reference
  *
  * @return always 1
@@ -1214,7 +1195,7 @@ int BrookBonded::setup( int numberOfParticles,
                         BrookBondParameters* periodicTorsionBrookBondParameters,
                         BrookBondParameters* rbTorsionBrookBondParameters,
                         BrookBondParameters* nonBonded14ForceParameters,  
-                        double lj14Scale, double coulombScale, int particleStreamWidth, int particleStreamSize ){
+                        int particleStreamWidth, int particleStreamSize ){
 
 // ---------------------------------------------------------------------------------------
 
@@ -1226,8 +1207,9 @@ int BrookBonded::setup( int numberOfParticles,
 
    if( PrintOn && getLog() ){
       (void) fprintf( getLog(), "%s particles=%d\n   [%p %p %p %p %p] (bond, angle, pd, rb, 14)\n"
-                      "14Scale=%f %f StreamW=%d StreamSz=%d\n", methodName.c_str(), numberOfParticles, harmonicBondBrookBondParameters, harmonicAngleBrookBondParameters,
-                      periodicTorsionBrookBondParameters, rbTorsionBrookBondParameters, nonBonded14ForceParameters, lj14Scale,coulombScale,
+                      "StreamW=%d StreamSz=%d\n", methodName.c_str(), numberOfParticles, harmonicBondBrookBondParameters,
+                      harmonicAngleBrookBondParameters,
+                      periodicTorsionBrookBondParameters, rbTorsionBrookBondParameters, nonBonded14ForceParameters,
                       particleStreamWidth, particleStreamSize ); fflush( getLog() );
    }
 
@@ -1293,7 +1275,7 @@ int BrookBonded::setup( int numberOfParticles,
       addBonds(      &nbondeds, particles, params, harmonicBondBrookBondParameters->getParticleIndices(),      harmonicBondBrookBondParameters->getBondParameters() );
    }
    if( nonBonded14ForceParameters ){
-      addPairs(      &nbondeds, particles, params, charges, nonBonded14ForceParameters->getParticleIndices(), nonBonded14ForceParameters->getBondParameters(), lj14Scale, coulombScale );
+      addPairs(      &nbondeds, particles, params, charges, nonBonded14ForceParameters->getParticleIndices(), nonBonded14ForceParameters->getBondParameters() );
    }
 
 // ---------------------------------------------------------------------------------------
@@ -1419,11 +1401,6 @@ int BrookBonded::setup( int numberOfParticles,
    delete[] particles;
    delete[] charges;
 
-   // set the fudge factors
-
-   //_ljScale        = (BrookOpenMMFloat) lj14Scale;
-   //_coulombFactor  = (BrookOpenMMFloat) coulombScale;
-
    // initialize output streams
 
    for( int ii = 0; ii < getNumberOfForceStreams(); ii++ ){
@@ -1473,12 +1450,6 @@ std::string BrookBonded::getContentsString( int level ) const {
    (void) LOCAL_SPRINTF( value, "%d", getNumberOfParticles() );
    message << _getLine( tab, "Number of particles:", value ); 
 
-   (void) LOCAL_SPRINTF( value, "%.5f", getLJ_14Scale() );
-   message << _getLine( tab, "LJ 14 scaling:", value ); 
-
-   (void) LOCAL_SPRINTF( value, "%.5f", getCoulombFactor() );
-   message << _getLine( tab, "Coulomb factor:", value ); 
-
    (void) LOCAL_SPRINTF( value, "%d", getInverseMapStreamWidth() );
    message << _getLine( tab, "Inverse map stream width:", value ); 
 
@@ -1835,7 +1806,7 @@ void BrookBonded::computeForces( BrookStreamImpl& positionStream, BrookStreamImp
 
    // bonded
 
-   float epsfac                                        = (float) (getLJ_14Scale()*getCoulombFactor());
+   float epsfac                                        = (float) (getCoulombFactor());
    float width                                         = (float) (getInverseMapStreamWidth());
 
    // bonded forces
@@ -1872,7 +1843,7 @@ void BrookBonded::computeForces( BrookStreamImpl& positionStream, BrookStreamImp
 
       int countPrintInvMap[4] = { 3, 5, 2, 4 }; 
 
-      (void) fprintf( getLog(), "\nPost kbonded_CDLJ: epsFac=%.6f %.6f %.6f", epsfac, getLJ_14Scale(), getCoulombFactor());
+      (void) fprintf( getLog(), "\nPost kbonded_CDLJ: epsFac=%.6f %.6f", epsfac, getCoulombFactor());
       (void) fprintf( getLog(), "\nParticle indices stream\n" );
       getParticleIndicesStream()->printToFile( getLog() );
 
@@ -1922,6 +1893,20 @@ void BrookBonded::computeForces( BrookStreamImpl& positionStream, BrookStreamImp
 
                          forceStream.getBrookStream(), forceStream.getBrookStream() );
 
+   } else if( getInverseMapStreamCount( I_Stream ) == 2 && getInverseMapStreamCount( K_Stream ) == 2 ){
+
+      kinvmap_gather2_2( width,
+
+                         inverseStreamMaps[I_Stream][0]->getBrookStream(),
+                         inverseStreamMaps[I_Stream][1]->getBrookStream(),
+                         bondedForceStreams[I_Stream]->getBrookStream(),
+
+                         inverseStreamMaps[K_Stream][0]->getBrookStream(),
+                         inverseStreamMaps[K_Stream][1]->getBrookStream(),
+                         bondedForceStreams[K_Stream]->getBrookStream(),
+
+                         forceStream.getBrookStream(), forceStream.getBrookStream() );
+
    } else if( getInverseMapStreamCount( I_Stream ) == 3 && getInverseMapStreamCount( K_Stream ) == 4 ){
 
       kinvmap_gather3_4( width,

platforms/brook/src/BrookBonded.h

@@ -79,8 +79,6 @@ class BrookBonded : public BrookCommon {
        * @param bonded14Indices           each element contains the indices of two particles whose nonbonded interactions should be reduced since
        *                                  they form a bonded 1-4 pair
        * @param nonbondedParameters       the nonbonded force parameters (charge, sigma, epsilon) for each particle
-       * @param lj14Scale                 the factor by which van der Waals interactions should be reduced for bonded 1-4 pairs
-       * @param coulomb14Scale            the factor by which Coulomb interactions should be reduced for bonded 1-4 pairs
        * @param log                       log reference
        *
        * @return nonzero value if error
@@ -94,7 +92,7 @@ class BrookBonded : public BrookCommon {
                  const std::vector<std::vector<int> >& periodicTorsionIndices, const std::vector<std::vector<double> >& periodicTorsionParameters,
                  const std::vector<std::vector<int> >& rbTorsionIndices,       const std::vector<std::vector<double> >& rbTorsionParameters,
                  const std::vector<std::vector<int> >& bonded14Indices,        const std::vector<std::vector<double> >& nonbondedParameters,
-                 double lj14Scale, double coulomb14Scale, const Platform& platform );
+                 const Platform& platform );
 */
 
       int setup( int numberOfParticles,
@@ -102,7 +100,7 @@ class BrookBonded : public BrookCommon {
                  BrookBondParameters* harmonicAngleBrookBondParameters,
                  BrookBondParameters* periodicTorsionBrookBondParameters,
                  BrookBondParameters* rbTorsionBrookBondParameters,
-                 BrookBondParameters* nonBonded14ForceParameters,  double lj14Scale, double coulombScale, int particleStreamWidth, int particleStreamSize );
+                 BrookBondParameters* nonBonded14ForceParameters,  int particleStreamWidth, int particleStreamSize );
 
       /**
        * Get inverse map stream width
@@ -168,15 +166,6 @@ class BrookBonded : public BrookCommon {
 
       BrookFloatStreamInternal* getBrookParticleIndices( void ) const; 
 
-      /**
-       * Get LJ 14 scale factor
-       *
-       * @return  LJ 14 scaling (fudge) factor
-       * 
-       */
-
-       BrookOpenMMFloat getLJ_14Scale( void ) const;
-
       /** 
        * Get Coulomb factor
        * 
@@ -312,7 +301,6 @@ class BrookBonded : public BrookCommon {
 
       // scale factors for 1-4 ixn
 
-      BrookOpenMMFloat _ljScale;
       BrookOpenMMFloat _coulombFactor;
 
       // streams
@@ -408,15 +396,13 @@ class BrookBonded : public BrookCommon {
        * @param bonded14Indices           each element contains the indices of two particles whose nonbonded interactions should be reduced since
        *                                  they form a bonded 1-4 pair
        * @param nonbondedParameters       the nonbonded force parameters (charge, sigma, epsilon) for each particle
-       * @param lj14Scale                 the factor by which van der Waals interactions should be reduced for bonded 1-4 pairs
        *
        * @return nonzero value if error
        *
        */
 
       int addPairs( int *nbondeds, int *particles, BrookOpenMMFloat* params[], BrookOpenMMFloat* charges,
-                    const std::vector<std::vector<int> >& bonded14Indices, const std::vector<std::vector<double> >& nonbondedParameters,
-                    double lj14Scale, double coulombScale );
+                    const std::vector<std::vector<int> >& bonded14Indices, const std::vector<std::vector<double> >& nonbondedParameters );
       
       /**
        * Create and load inverse maps for bonded ixns

platforms/brook/src/BrookBrownianDynamics.cpp

@@ -524,10 +524,9 @@ int BrookBrownianDynamics::update( Stream& positions, Stream& velocities,
 
    if( brookShakeAlgorithm.getNumberOfConstraints() > 0 ){
       kshakeh_fix1( 
-                    10.0f,
+                    25.0f,
                     (float) getBrownianDynamicsParticleStreamWidth(),
                     brookShakeAlgorithm.getInverseHydrogenMass(),
-                    omega, 
                     brookShakeAlgorithm.getShakeParticleIndicesStream()->getBrookStream(),
                     positionStream.getBrookStream(),
                     getXPrimeStream()->getBrookStream(),

platforms/brook/src/BrookCalcHarmonicAngleForceKernel.h

@@ -152,7 +152,6 @@ class BrookCalcHarmonicAngleForceKernel : public CalcHarmonicAngleForceKernel {
 
       FILE* _log;
 
-
       // Brook bond parameters
 
       BrookBondParameters* _brookBondParameters;

platforms/brook/src/BrookCalcNonbondedForceKernel.cpp

@@ -240,7 +240,7 @@ void BrookCalcNonbondedForceKernel::initialize14Interactions( const System& syst
    
          _brookBondParameters->setBond( ii, particles, parameters );
       }   
-      _openMMBrookInterface.setNonBonded14ForceParameters( _brookBondParameters, 1.0, 1.0 );
+      _openMMBrookInterface.setNonBonded14ForceParameters( _brookBondParameters );
    
       if( log ){
          std::string contents = _brookBondParameters->getContentsString( ); 

platforms/brook/src/BrookIntegrateLangevinStepKernel.cpp

@@ -58,6 +58,12 @@ BrookIntegrateLangevinStepKernel::BrookIntegrateLangevinStepKernel( std::string
    _brookLangevinDynamics          = NULL;
    _brookShakeAlgorithm            = NULL;
    _brookRandomNumberGenerator     = NULL;
+   _log                            = NULL;
+
+   const BrookPlatform brookPlatform        = dynamic_cast<const BrookPlatform&> (platform);
+   if( brookPlatform.getLog() != NULL ){
+      setLog( brookPlatform.getLog() );
+   }
 
 }
 
@@ -80,6 +86,31 @@ BrookIntegrateLangevinStepKernel::~BrookIntegrateLangevinStepKernel( ){
 
 }
 
+/** 
+ * Get log file reference
+ * 
+ * @return  log file reference
+ *
+ */
+
+FILE* BrookIntegrateLangevinStepKernel::getLog( void ) const {
+   return _log;
+}
+
+/** 
+ * Set log file reference
+ * 
+ * @param  log file reference
+ *
+ * @return  DefaultReturnValue
+ *
+ */
+
+int BrookIntegrateLangevinStepKernel::setLog( FILE* log ){
+   _log = log;
+   return DefaultReturnValue;
+}
+
 /** 
  * Initialize the kernel, setting up all parameters related to integrator.
  * 
@@ -92,10 +123,17 @@ void BrookIntegrateLangevinStepKernel::initialize( const System& system, const L
 
 // ---------------------------------------------------------------------------------------
 
-   // static const std::string methodName      = "BrookIntegrateLangevinStepKernel::initialize";
+   static const int printOn                  = 1;
+   static const std::string methodName       = "BrookIntegrateLangevinStepKernel::initialize";
 
 // ---------------------------------------------------------------------------------------
    
+   FILE* log             = getLog();
+   if( printOn && log ){
+      (void) fprintf( log, "%s\n", methodName.c_str() );
+      (void) fflush( log );
+   }
+
    int numberOfParticles = system.getNumParticles();
 
    // masses
@@ -103,14 +141,24 @@ void BrookIntegrateLangevinStepKernel::initialize( const System& system, const L
    std::vector<double> masses;
    masses.resize( numberOfParticles );
 
+   if( printOn && log ){
+      (void) fprintf( log, "%s %d\n", methodName.c_str(), numberOfParticles );
+      (void) fflush( log );
+   }
+
    for( int ii = 0; ii < numberOfParticles; ii++ ){
-      masses[ii] = static_cast<RealOpenMM>(system.getParticleMass(ii));
+      masses[ii] = static_cast<double>(system.getParticleMass(ii));
    }
 
    // constraints
 
    int numberOfConstraints = system.getNumConstraints();
 
+   if( printOn && log ){
+      (void) fprintf( log, "%s const=%d\n", methodName.c_str(), numberOfConstraints );
+      (void) fflush( log );
+   }
+
    std::vector<std::vector<int> > constraintIndicesVector;
    constraintIndicesVector.resize( numberOfConstraints );
 
@@ -122,13 +170,15 @@ void BrookIntegrateLangevinStepKernel::initialize( const System& system, const L
       int particle1, particle2;
       double distance;
 
+(void) fprintf( log, "%s shake setup const=%d ", methodName.c_str(), ii ); fflush( log );
       system.getConstraintParameters( ii, particle1, particle2, distance );
-      std::vector<int> constraintIndices;
       constraintIndicesVector[ii]  = constraintIndices;
 
-      constraintIndices[0]         = particle1;
-      constraintIndices[1]         = particle2;
-      constraintLengths[ii]        = static_cast<RealOpenMM>(distance);
+(void) fprintf( log, "[ %d %d %f]\n", methodName.c_str(), particle1, particle2, distance ); fflush( log );
+      constraintIndicesVector[ii].push_back( particle1 );
+      constraintIndicesVector[ii].push_back([particle2 );
+      constraintIndicesVector[ii].push_back([particle2 );
+      constraintLengths.push_back( distance );
    }
 
    _brookLangevinDynamics        = new BrookLangevinDynamics( );
@@ -139,9 +189,13 @@ void BrookIntegrateLangevinStepKernel::initialize( const System& system, const L
 
    // assert( (_brookShakeAlgorithm->getNumberOfConstraints() > 0) );
 
+   if( printOn && log ){
+      (void) fprintf( log, "%s done shake setup const=%d\n", methodName.c_str(), numberOfConstraints );
+      (void) fflush( log );
+   }
+
    _brookRandomNumberGenerator   = new BrookRandomNumberGenerator( );
    _brookRandomNumberGenerator->setup( (int) masses.size(), getPlatform() );
-//   _brookRandomNumberGenerator->setVerbosity( 1 );
 
 }
 
@@ -174,10 +228,10 @@ void BrookIntegrateLangevinStepKernel::execute( OpenMMContextImpl& context, cons
    differences[1] = integrator.getFriction()    - (double) _brookLangevinDynamics->getFriction();
    differences[2] = integrator.getStepSize()    - (double) _brookLangevinDynamics->getStepSize();
    if( fabs( differences[0] ) > epsilon || fabs( differences[1] ) > epsilon || fabs( differences[2] ) > epsilon ){
-//printf( "%s calling updateParameters\n", methodName.c_str() );
+printf( "%s calling updateParameters\n", methodName.c_str() );
       _brookLangevinDynamics->updateParameters( integrator.getTemperature(), integrator.getFriction(), integrator.getStepSize() );
    } else {
-//printf( "%s NOT calling updateParameters\n", methodName.c_str() );
+printf( "%s NOT calling updateParameters\n", methodName.c_str() );
 }
 
    _brookLangevinDynamics->update( *(_openMMBrookInterface.getParticlePositions()), *(_openMMBrookInterface.getParticleVelocities()), 

platforms/brook/src/BrookIntegrateLangevinStepKernel.h

@@ -89,8 +89,41 @@ class BrookIntegrateLangevinStepKernel : public IntegrateLangevinStepKernel {
 
       void execute( OpenMMContextImpl& context, const LangevinIntegrator& integrator );
 
+      /** 
+       * Set log file reference
+       * 
+       * @param  log file reference
+       *
+       * @return DefaultReturnValue
+       *
+       */
+          
+      int setLog( FILE* log );
+
+      /*  
+       * Get contents of object
+       *
+       * @param level of dump
+       *
+       * @return string containing contents
+       *
+       * */
+          
+      std::string getContents( int level ) const;
+
+      /** 
+       * Get log file reference
+       * 
+       * @return  log file reference
+       *
+       */
+          
+      FILE* getLog( void ) const;
+
    protected:
 
+      FILE* _log;
+
       BrookLangevinDynamics*        _brookLangevinDynamics;
       BrookShakeAlgorithm*          _brookShakeAlgorithm;
       BrookRandomNumberGenerator*   _brookRandomNumberGenerator;

platforms/brook/src/BrookKernelFactory.cpp

@@ -126,7 +126,7 @@ KernelImpl* BrookKernelFactory::createKernelImpl( std::string name, const Platfo
 
 	} else if( name == RemoveCMMotionKernel::Name() ){
 
-//      return new BrookRemoveCMMotionKernel( name, platform, openMMBrookInterface );
+      return new BrookRemoveCMMotionKernel( name, platform, openMMBrookInterface, context.getSystem() );
 
    // KE calculator
 

platforms/brook/src/BrookLangevinDynamics.cpp

@@ -371,11 +371,11 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
 
    // unused Shake parameter
 
-   float omega                   = 1.0f;
+   float numberOfIterations      = 25.0f;
 
    static const char* methodName = "\nBrookLangevinDynamics::update";
 
-   static const int PrintOn      = 0;
+   static const int PrintOn      = 1;
 
 // ---------------------------------------------------------------------------------------
 
@@ -482,10 +482,9 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
 
    if( brookShakeAlgorithm.getNumberOfConstraints() > 0 ){
       kshakeh_fix1( 
-                    10.0f,
+                    numberOfIterations,
                     (float) getLangevinDynamicsParticleStreamWidth(),
                     brookShakeAlgorithm.getInverseHydrogenMass(),
-                    omega, 
                     brookShakeAlgorithm.getShakeParticleIndicesStream()->getBrookStream(),
                     positionStream.getBrookStream(),
                     getXPrimeStream()->getBrookStream(),
@@ -575,10 +574,9 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
 
    if( brookShakeAlgorithm.getNumberOfConstraints() > 0 ){
       kshakeh_fix1( 
-                    10.0f,
+                    numberOfIterations,
                     (float) getLangevinDynamicsParticleStreamWidth(),
                     brookShakeAlgorithm.getInverseHydrogenMass(),
-                    omega, 
                     brookShakeAlgorithm.getShakeParticleIndicesStream()->getBrookStream(),
                     positionStream.getBrookStream(),
                     getXPrimeStream()->getBrookStream(),

platforms/brook/src/BrookRemoveCMMotionKernel.cpp

@@ -29,6 +29,8 @@
  * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
  * -------------------------------------------------------------------------- */
 
+#include "internal/OpenMMContextImpl.h"
+#include "System.h"
 #include "BrookRemoveCMMotionKernel.h"
 #include "BrookStreamInternal.h"
 
@@ -40,11 +42,13 @@ using namespace std;
  * 
  * @param name        name of the stream to create
  * @param platform    platform
+ * @param openMMBrookInterface      OpenMM-Brook interface
+ * @param System                    System reference
  *
  */
   
-BrookRemoveCMMotionKernel::BrookRemoveCMMotionKernel( std::string name, const Platform& platform ) :
-                              RemoveCMMotionKernel( name, platform ){
+BrookRemoveCMMotionKernel::BrookRemoveCMMotionKernel( std::string name, const Platform& platform, OpenMMBrookInterface& openMMBrookInterface, System& system ) :
+                              RemoveCMMotionKernel( name, platform ), _openMMBrookInterface( openMMBrookInterface ), _system( system ){
 
 // ---------------------------------------------------------------------------------------
 
@@ -52,8 +56,16 @@ BrookRemoveCMMotionKernel::BrookRemoveCMMotionKernel( std::string name, const Pl
 
 // ---------------------------------------------------------------------------------------
 
+   _frequency                        = 0;
+   _log                              = NULL;
    _brookVelocityCenterOfMassRemoval = NULL;
 
+   const BrookPlatform brookPlatform = dynamic_cast<const BrookPlatform&> (platform);
+   if( brookPlatform.getLog() != NULL ){
+      setLog( brookPlatform.getLog() );
+   }   
+   _openMMBrookInterface.setNumberOfParticles( system.getNumParticles() );
+
 }
 
 /** 
@@ -75,11 +87,12 @@ BrookRemoveCMMotionKernel::~BrookRemoveCMMotionKernel( ){
 /** 
  * Initialize the kernel
  * 
- * @param masses   array of particle masses
+ * @param system   System reference
+ * @param force    CMMotionRemover reference
  *
  */
 
-void BrookRemoveCMMotionKernel::initialize( const vector<double>& masses ){
+void BrookRemoveCMMotionKernel::initialize( const System& system, const CMMotionRemover& force ){
 
 // ---------------------------------------------------------------------------------------
 
@@ -87,6 +100,13 @@ void BrookRemoveCMMotionKernel::initialize( const vector<double>& masses ){
 
 // ---------------------------------------------------------------------------------------
 
+   _frequency = force.getFrequency();
+   std::vector<double> masses;
+   masses.resize( system.getNumParticles() );
+   for( size_t ii = 0; ii < masses.size(); ii++ ){
+      masses[ii] = system.getParticleMass(ii);
+   }
+
    _brookVelocityCenterOfMassRemoval = new BrookVelocityCenterOfMassRemoval();
    _brookVelocityCenterOfMassRemoval->setup( masses, getPlatform() );
 
@@ -94,14 +114,51 @@ void BrookRemoveCMMotionKernel::initialize( const vector<double>& masses ){
 
 }
 
+/** 
+ * Get log file reference
+ * 
+ * @return  log file reference
+ *
+ */
+
+FILE* BrookRemoveCMMotionKernel::getLog( void ) const {
+   return _log;
+}
+
+/** 
+ * Set log file reference
+ * 
+ * @param  log file reference
+ *
+ * @return  DefaultReturnValue
+ *
+ */
+
+int BrookRemoveCMMotionKernel::setLog( FILE* log ){
+   _log = log;
+   return BrookCommon::DefaultReturnValue;
+}
+
+
+/** 
+ * Get COM removal frequency
+ * 
+ * @return frequency 
+ *
+ */
+
+int BrookRemoveCMMotionKernel::getFrequency( void ) const {
+   return _frequency;
+}
+
 /** 
  * Execute kernel
  * 
- * @param velocities  array of particle velocities
+ * @param context    OpenMMContextImpl reference
  *
  */
 
-void BrookRemoveCMMotionKernel::execute( Stream& velocities ){
+void BrookRemoveCMMotionKernel::execute( OpenMMContextImpl& context ){
 
 // ---------------------------------------------------------------------------------------
 
@@ -109,7 +166,14 @@ void BrookRemoveCMMotionKernel::execute( Stream& velocities ){
 
 // ---------------------------------------------------------------------------------------
 
-   _brookVelocityCenterOfMassRemoval->removeVelocityCenterOfMass( velocities );
+   int step      = (int) floor( context.getTime()/context.getIntegrator().getStepSize() );
+   int frequency = getFrequency();
+   if( frequency <= 0 || (step % frequency) != 0 ){
+        return;
+   }
+
+   BrookStreamImpl* velocities = _openMMBrookInterface.getParticleVelocities();
+   _brookVelocityCenterOfMassRemoval->removeVelocityCenterOfMass( *velocities );
 
    return;
 }

platforms/brook/src/BrookRemoveCMMotionKernel.h

@@ -33,6 +33,7 @@
  * -------------------------------------------------------------------------- */
 
 #include "kernels.h"
+#include "OpenMMBrookInterface.h"
 #include "BrookVelocityCenterOfMassRemoval.h"
 
 namespace OpenMM {
@@ -50,10 +51,12 @@ class BrookRemoveCMMotionKernel : public RemoveCMMotionKernel {
        * 
        * @param name                      name of the stream to create
        * @param platform                  platform
+       * @param openMMBrookInterface      OpenMM-Brook interface
+       * @param System                    System reference
        *
        */
   
-      BrookRemoveCMMotionKernel( std::string name, const Platform& platform );
+      BrookRemoveCMMotionKernel( std::string name, const Platform& platform, OpenMMBrookInterface& openMMBrookInterface, System& system );
 
       /**
        * BrookRemoveCMMotionKernel destructor
@@ -62,26 +65,82 @@ class BrookRemoveCMMotionKernel : public RemoveCMMotionKernel {
   
       ~BrookRemoveCMMotionKernel();
 
+      /** 
+       * Get COM removal frequency
+       * 
+       * @return frequency 
+       *
+       */
+      
+      int getFrequency( void ) const;
+       
+      /** 
+       * Set log file reference
+       * 
+       * @param  log file reference
+       *
+       * @return DefaultReturnValue
+       *
+       */
+ 
+      int setLog( FILE* log );
+
+      /*  
+       * Get contents of object
+       *
+       * @param level of dump
+       *
+       * @return string containing contents
+       *
+       * */
+          
+      std::string getContents( int level ) const;
+
+      /** 
+       * Get log file reference
+       * 
+       * @return  log file reference
+       *
+       */
+    
+      FILE* getLog( void ) const;
+
       /** 
        * Initialize the kernel
        * 
-       * @param masses   mass of each particle
+       * @param system   System reference
+       * @param force    CMMotionRemover reference
        *
        */
-      void initialize( const std::vector<double>& masses );
+      void initialize( const System& system, const CMMotionRemover& force );
 
       /** 
        * Execute the kernel.
        * 
-       * @param velocities  stream of particle velocities
+       * @param context  OpenMMContextImpl reference
        *
        */
 
-      void execute( Stream& velocities );
+      void execute( OpenMMContextImpl& context );
 
    private:
 
+      int _frequency;
+
+      // log file reference
+
+      FILE* _log;
+
       BrookVelocityCenterOfMassRemoval* _brookVelocityCenterOfMassRemoval;
+
+      // interface
+
+      OpenMMBrookInterface& _openMMBrookInterface;
+
+      // System reference
+
+      System& _system;
+
 };
 
 } // namespace OpenMM

platforms/brook/src/BrookVelocityCenterOfMassRemoval.cpp

@@ -143,7 +143,7 @@ BrookFloatStreamInternal* BrookVelocityCenterOfMassRemoval::getLinearMomentumStr
  *
  */
    
-int BrookVelocityCenterOfMassRemoval::removeVelocityCenterOfMass( Stream& velocities ){
+int BrookVelocityCenterOfMassRemoval::removeVelocityCenterOfMass( BrookStreamImpl& velocityStream ){
 
    // ---------------------------------------------------------------------------------------
 
@@ -154,12 +154,11 @@ int BrookVelocityCenterOfMassRemoval::removeVelocityCenterOfMass( Stream& veloci
 
    if( debug && getLog() ){
       BrookOpenMMFloat com[3];
-      getVelocityCenterOfMass( velocities, com );
+      getVelocityCenterOfMass( velocityStream, com );
       (void) fprintf( getLog(), "\n%s Pre removal com: [%12.5e %12.5e %12.5e]\n", methodName, com[0], com[1], com[2] );
-      BrookStreamImpl& v1                  = dynamic_cast<BrookStreamImpl&> (velocities.getImpl());
-      BrookStreamInternal* velocityStream  = dynamic_cast<BrookStreamInternal*> (v1.getBrookStreamImpl());
+      BrookStreamInternal* outputVelocityStream = dynamic_cast<BrookStreamInternal*> (velocityStream.getBrookStreamImpl());
 
-      void* velV                       = velocityStream->getData( 1 );
+      void* velV                       = outputVelocityStream->getData( 1 );
       const float* vArray              = (float*) velV;
 
       int index                        = 0;
@@ -174,8 +173,6 @@ int BrookVelocityCenterOfMassRemoval::removeVelocityCenterOfMass( Stream& veloci
    // calculate linear momentum via reduction
    // subtract it (/totalMass) from velocities
 
-   BrookStreamImpl& velocityStream  = dynamic_cast<BrookStreamImpl&> (velocities.getImpl());
-
    kCalculateLinearMomentum( getMassStream()->getBrookStream(), velocityStream.getBrookStream(), getWorkStream()->getBrookStream() );
    kSumLinearMomentum( (float) getComParticleStreamWidth(), (float) getNumberOfParticles(), getWorkStream()->getBrookStream(), getLinearMomentumStream()->getBrookStream() );
    kScale( (float) getTotalInverseMass(), getLinearMomentumStream()->getBrookStream(), getLinearMomentumStream()->getBrookStream() );
@@ -183,7 +180,7 @@ int BrookVelocityCenterOfMassRemoval::removeVelocityCenterOfMass( Stream& veloci
 
    if( (0 || debug) && getLog() ){
       BrookOpenMMFloat com[3];
-      getVelocityCenterOfMass( velocities, com );
+      getVelocityCenterOfMass( velocityStream, com );
       (void) fprintf( getLog(), "%s strW=%d iatm=%d Post removal com: [%12.5e %12.5e %12.5e]", methodName,
                       getComParticleStreamWidth(), getNumberOfParticles(),  com[0], com[1], com[2] );
 
@@ -191,10 +188,9 @@ int BrookVelocityCenterOfMassRemoval::removeVelocityCenterOfMass( Stream& veloci
       float* linMo = (float*) linMoV;
       (void) fprintf( getLog(), "LM [%12.5e %12.5e %12.5e]\n", linMo[0], linMo[1], linMo[2] );
 
-      BrookStreamImpl& v1                  = dynamic_cast<BrookStreamImpl&> (velocities.getImpl());
-      BrookStreamInternal* velocityStream  = dynamic_cast<BrookStreamInternal*> (v1.getBrookStreamImpl());
+      BrookStreamInternal* outputVelocityStream = dynamic_cast<BrookStreamInternal*> (velocityStream.getBrookStreamImpl());
 
-      void* velV                       = velocityStream->getData( 1 );
+      void* velV                       = outputVelocityStream->getData( 1 );
       const float* vArray              = (float*) velV;
 
       void* w1                        = getWorkStream()->getData( 1 );
@@ -225,7 +221,7 @@ int BrookVelocityCenterOfMassRemoval::removeVelocityCenterOfMass( Stream& veloci
  *
  */
    
-int BrookVelocityCenterOfMassRemoval::getVelocityCenterOfMass( Stream& velocities, BrookOpenMMFloat velocityCom[3] ){
+int BrookVelocityCenterOfMassRemoval::getVelocityCenterOfMass( BrookStreamImpl& vStream, BrookOpenMMFloat velocityCom[3] ){
 
    // ---------------------------------------------------------------------------------------
 
@@ -239,8 +235,7 @@ int BrookVelocityCenterOfMassRemoval::getVelocityCenterOfMass( Stream& velocitie
    // calculate linear momentum via reduction
    // subtract it (/totalMass) from velocities
 
-   BrookStreamImpl& v1                  = dynamic_cast<BrookStreamImpl&> (velocities.getImpl());
-   BrookStreamInternal* velocityStream  = dynamic_cast<BrookStreamInternal*> (v1.getBrookStreamImpl());
+   BrookStreamInternal* velocityStream  = dynamic_cast<BrookStreamInternal*> (vStream.getBrookStreamImpl());
 
    void* velV                           = velocityStream->getData( 1 );
    const float* vArray                  = (float*) velV;

platforms/brook/src/BrookVelocityCenterOfMassRemoval.h

@@ -35,7 +35,8 @@
 #include <vector>
 #include <set>
 
-#include "BrookFloatStreamInternal.h"
+#include "BrookStreamImpl.h"
+//#include "BrookFloatStreamInternal.h"
 #include "BrookPlatform.h"
 #include "BrookCommon.h"
 
@@ -98,7 +99,7 @@ class BrookVelocityCenterOfMassRemoval : public BrookCommon {
        *
        */
       
-      int removeVelocityCenterOfMass( Stream& velocities );
+      int removeVelocityCenterOfMass( BrookStreamImpl& velocityStream );
 
       /** 
        * Get velocity center-of-mass
@@ -110,7 +111,7 @@ class BrookVelocityCenterOfMassRemoval : public BrookCommon {
        *
        */
       
-      int getVelocityCenterOfMass( Stream& velocities, BrookOpenMMFloat velocityCom[3] );
+      int getVelocityCenterOfMass( BrookStreamImpl& vStream, BrookOpenMMFloat velocityCom[3] );
 
       /* 
        * Setup of  parameters

platforms/brook/src/BrookVerletDynamics.cpp

@@ -185,7 +185,7 @@ int BrookVerletDynamics::update( Stream& positions, Stream& velocities,
 
    // unused Shake parameter
 
-   float omega                   = 1.0f;
+   float numberOfIterations       = 25.0f;
 
    static const char* methodName  = "\nBrookVerletDynamics::update";
 
@@ -261,10 +261,9 @@ int BrookVerletDynamics::update( Stream& positions, Stream& velocities,
    if( brookShakeAlgorithm.getNumberOfConstraints() > 0 ){
 
       kshakeh_fix1( 
-                    10.0f,
+                    25.0f,
                     (float) getVerletDynamicsParticleStreamWidth(),
                     brookShakeAlgorithm.getInverseHydrogenMass(),
-                    omega, 
                     brookShakeAlgorithm.getShakeParticleIndicesStream()->getBrookStream(),
                     positionStream.getBrookStream(),
                     getXPrimeStream()->getBrookStream(),

platforms/brook/src/OpenMMBrookInterface.cpp

@@ -78,9 +78,6 @@ OpenMMBrookInterface::OpenMMBrookInterface( int streamWidth ) : _particleStreamW
    _referencePlatform                       = NULL;
    _refVerletIntegrator                     = NULL;
 
-   _lj14Scale                               = 1.0;
-   _coulomb14Scale                          = 1.0;
-
    _particleStreamSize                      = -1;
 
    for( int ii = 0; ii < LastBondForce; ii++ ){
@@ -165,28 +162,6 @@ int OpenMMBrookInterface::getParticleStreamSize( void ) const {
    return _particleStreamSize;
 }
 
-/** 
- * Get LJ-14 scale factor
- * 
- * @return    LJ-14 scale factor
- *
- */
-
-double OpenMMBrookInterface::getLj14Scale( void ) const {
-   return _lj14Scale;
-}
-
-/** 
- * Get Coulomb scale factor
- * 
- * @return    Coulomb scale factor
- *
- */
-
-double OpenMMBrookInterface::getCoulomb14Scale( void ) const {
-   return _coulomb14Scale;
-}
-
 /** 
  * Get log file reference
  * 
@@ -367,17 +342,12 @@ BrookBondParameters* OpenMMBrookInterface::getNonBonded14ForceParameters( void )
  * Set BrookBondParameters for LJ 14 force
  * 
  * @param brookBondParameters brookBondParameters for LJ 14 force
- * @param lj14Scale           LJ 14 scale factor
- * @param coulomb14Scale        Coulomb 14 scale factor
  *
  * @return  DefaultReturnValue
  *
  */
 
-int OpenMMBrookInterface::setNonBonded14ForceParameters( BrookBondParameters* brookBondParameters, 
-                                                         double lj14Scale, double coulomb14Scale ){
-   _lj14Scale        = lj14Scale;
-   _coulomb14Scale   = coulomb14Scale;
+int OpenMMBrookInterface::setNonBonded14ForceParameters( BrookBondParameters* brookBondParameters ){
    return _setBondParameters( LJ14Index, brookBondParameters );
 }
 
@@ -584,19 +554,17 @@ void OpenMMBrookInterface::computeForces( OpenMMContextImpl& context ){
 (void) fprintf( stderr, "%s done nonbonded: bonded=%d completed=%d\n", methodName.c_str(), _brookBonded.isActive(), _brookBonded.isSetupCompleted() ); (void) fflush( stderr );
    if( _brookBonded.isActive() ){
 
-//(void) fprintf( stderr, "%s Bonded\n", methodName.c_str() ); (void) fflush( stderr );
       // perform setup first time through
 
+      //if( _brookBonded.isSetupCompleted() == 0 ){
       if( _brookBonded.isSetupCompleted() >= -1 ){
          _brookBonded.setup( getNumberOfParticles(), getHarmonicBondForceParameters(), getHarmonicAngleForceParameters(),
                              getPeriodicTorsionForceParameters(), getRBTorsionForceParameters(),
                              getNonBonded14ForceParameters(),
-                             getLj14Scale(), getCoulomb14Scale(), getParticleStreamWidth(), getParticleStreamSize() );
+                             getParticleStreamWidth(), getParticleStreamSize() );
       }
 
-(void) fprintf( stderr, "%s done setup bonded=%d completed=%d\n", methodName.c_str(), _brookBonded.isActive(), _brookBonded.isSetupCompleted() ); (void) fflush( stderr );
       _brookBonded.computeForces( *positions, *forces );
-(void) fprintf( stderr, "%s done forces bonded=%d completed=%d\n", methodName.c_str(), _brookBonded.isActive(), _brookBonded.isSetupCompleted() ); (void) fflush( stderr );
    
       // diagnostics
    

platforms/brook/src/OpenMMBrookInterface.h

@@ -128,24 +128,6 @@ class OpenMMBrookInterface {
       
       int getParticleStreamSize( void ) const;
       
-      /** 
-       * Get LJ-14 scale factor
-       * 
-       * @return    LJ-14 scale factor
-       *
-       */
-      
-      double getLj14Scale( void ) const;
-      
-      /** 
-       * Get Coulomb scale factor
-       * 
-       * @return    Coulomb scale factor
-       *
-       */
-      
-      double getCoulomb14Scale( void ) const;
-
       /** 
        * Get log file reference
        * 
@@ -332,14 +314,12 @@ class OpenMMBrookInterface {
        * Set BrookBondParameters for LJ 14 force
        * 
        * @param brookBondParameters brookBondParameters for LJ 14 force
-       * @param lj14Scale           LJ 14 scale factor
-       * @param coulomb14Scale        Coulomb 14 scale factor
        *
        * @return  DefaultReturnValue
        *
        */
       
-      int setNonBonded14ForceParameters( BrookBondParameters* brookBondParameters, double lj14Scale, double coulomb14Scale );
+      int setNonBonded14ForceParameters( BrookBondParameters* brookBondParameters );
       
       /** 
        * Get positions stream
@@ -419,9 +399,6 @@ class OpenMMBrookInterface {
 
        int _numberOfParticles;
    
-       double _lj14Scale;
-       double _coulomb14Scale;
-
        // Brook bonded, nonbonded, Gbsa
 
        BrookBonded     _brookBonded;

platforms/brook/src/gpu/kinvmap_gather.br

@@ -249,6 +249,29 @@ kernel void kinvmap_gather6(
 }
 
 
+//Takes three + four inverse maps
+kernel void kinvmap_gather2_2(
+		float strwidth,    //stream width of the dihedral forces
+		float4 invmap3_1<>,   //indices into the dihedral forces
+		float4 invmap3_2<>,   //indices into the dihedral forces
+		float3 forces3[][], //dihedral forces
+		float4 invmap4_1<>,   //indices into the dihedral forces
+		float4 invmap4_2<>,   //indices into the dihedral forces
+		float3 forces4[][], //dihedral forces
+		float3 inforce<>,  //particle forces before
+		out float3 outforce<> //particle forces after
+		)
+{
+	outforce = inforce;
+
+	outforce += do_gather_nobranch( strwidth, invmap3_1, forces3 );
+	outforce += do_gather_nobranch( strwidth, invmap3_2, forces3 );
+    
+	outforce += do_gather_nobranch( strwidth, invmap4_1, forces4 );
+	outforce += do_gather_nobranch( strwidth, invmap4_2, forces4 );
+}
+
+
 //Takes three + four inverse maps
 kernel void kinvmap_gather3_3(
 		float strwidth,    //stream width of the dihedral forces

platforms/brook/src/gpu/kinvmap_gather.h

@@ -81,6 +81,16 @@ void  kinvmap_gather6 (const float  strwidth,
 		::brook::stream outforce); 
 
 
+void  kinvmap_gather2_2 (const float  strwidth,
+		::brook::stream invmap3_1,
+		::brook::stream invmap3_2,
+		::brook::stream forces3,
+		::brook::stream invmap4_1,
+		::brook::stream invmap4_2,
+		::brook::stream forces4,
+		::brook::stream inforce,
+		::brook::stream outforce);
+
 void  kinvmap_gather3_3 (const float  strwidth,
 		::brook::stream invmap3_1,
 		::brook::stream invmap3_2,

platforms/brook/src/gpu/kshakeh.br

@@ -43,169 +43,6 @@
  * to avoid stalling
  * *************************************************************/
 
-/*
-kernel void
-kshakeh(
-		float nit,      //number of iterations
-		float strwidth, //stream width of posq
-		float invmH,    //inverse mass of hydrogen
-		float omega,    //parameter from gromacs, possibly unused
-		float4 atoms<>, //heavy0, h1, h2, h3
-		float4 posq[][], //positions before update
-		float4 posqp[][], //positions after update
-		float4 params<>, // (1/m0, mu/2, blensq, 0.0f )
-		out float4 cposq0<>, //constrained position for heavy atom
-		out float4 cposq1<>, //ditto for h1
-		out float4 cposq2<>, //ditto for h2
-		out float4 cposq3<> //ditto for h3
-		) {
-	float2 ai, aj1, aj2, aj3; //2d indices, can be precalc.
-	float i; //iteration count
-	float3 xi, xj1, xj2, xj3; //coordinates
-	float3 xpi, xpj1, xpj2, xpj3; //coordinates
-	float3 rij1, rij2, rij3, rpij, dr;
-	float rpsqij, rrpr, acor;
-	float mask2, mask3;
-	float diff;
-	
-	ai.y  = floor( atoms.x / strwidth );
-	ai.x  = atoms.x - ai.y * strwidth;
-
-	aj1.y = floor( atoms.y / strwidth );
-	aj1.x = atoms.y - aj1.y * strwidth;
-
-	//If further hydrogens are absent,
-	//just set to the coordinates of the first
-	//so we don't have any junk memory accesses
-	//or nans/infs in the calcs.	
-	if ( atoms.z > -0.5f ) {
-		aj2.y = floor( atoms.z / strwidth );
-		aj2.x = atoms.z - aj2.y * strwidth;
-		mask2 = 1.0f;
-	}
-	else {
-		aj2 = aj1; 
-		mask2 = 0.0f;
-	}
-
-	if ( atoms.w > -0.5f ) {	
-		aj3.y = floor( atoms.w / strwidth );
-		aj3.x = atoms.w - aj3.y * strwidth;
-		mask3 = 1.0f;
-	}
-	else {
-		aj3 = aj1;
-		mask3 = 0.0f;
-	}
-
-
-	cposq0 = posq[ai];
-	cposq1 = posq[aj1];
-	cposq2 = posq[aj2];
-	cposq3 = posq[aj3];
-	
-	
-	xi  = cposq0.xyz;
-	xj1 = cposq1.xyz;
-	xj2 = cposq2.xyz;
-	xj3 = cposq3.xyz;
-	
-	rij1 = xi - xj1;
-	rij2 = xi - xj2;
-	rij3 = xi - xj3;
-
-	xpi  = posqp[ai].xyz;
-	xpj1 = posqp[aj1].xyz;
-	xpj2 = posqp[aj2].xyz;
-	xpj3 = posqp[aj3].xyz;
-
-	i = 0.0f;
-	while ( i < 15.0f ) {
-		//First hydrogen
-		rpij   = xpi - xpj1;
-		rpsqij = dot( rpij, rpij );
-		rrpr   = dot( rij1, rpij );
-		//for debugging only
-		diff = abs( params.z - rpsqij ) / (params.z * 2 * 1e-4 );
-		if ( diff < 1.0f )
-			acor = 0.0f;
-		else
-			//params.y = mu/2, params.z = blen*blen
-			acor = omega * ( params.z - rpsqij ) * params.y / rrpr ;
-		dr   = rij1 * acor;
-		xpi  += dr * params.x;
-		xpj1 -= dr * invmH;
-
-		//Second hydrogen
-		rpij   = xpi - xpj2;
-		rpsqij = dot( rpij, rpij );
-		rrpr   = dot( rij2, rpij );
-		//for debugging only
-		diff = abs( params.z - rpsqij ) / (params.z * 2 * 1e-4 );
-		if ( diff < 1.0f )
-			acor = 0.0f;
-		else
-			//params.y = mu/2, params.z = blen*blen
-			acor = mask2 * omega * ( params.z - rpsqij ) * params.y / rrpr ;
-		dr   = rij2 * acor;
-		xpi  += dr * params.x;
-		xpj2 -= dr * invmH;
-
-		//Third hydrogen
-		rpij   = xpi - xpj3;
-		rpsqij = dot( rpij, rpij );
-		rrpr   = dot( rij3, rpij );
-		//for debugging only
-		diff = abs( params.z - rpsqij ) / (params.z * 2 * 1e-4 );
-		if ( diff < 1.0f )
-			acor = 0.0f;
-		else
-			//params.y = mu/2, params.z = blen*blen
-			acor = mask3 * omega * ( params.z - rpsqij ) * params.y / rrpr ;
-		dr   = rij3 * acor;
-		xpi  += dr * params.x;
-		xpj3 -= dr * invmH;
-		
-		i += 1.0f;
-	}
-
-	cposq0.xyz = xpi;
-	cposq1.xyz = xpj1;
-	cposq2.xyz = xpj2;
-	cposq3.xyz = xpj3;
-} */
-
-/*Applies the updated positions
- *Each atom occurs at one and only one position
- * */
-/*
-kernel void kshakeh_update( 
-		float strwidth, //width of cposq streams
-		float2 invmap<>,    //shakeh inverse map
-		float4 posq<>,     //untouched positions
-		float4 cposq0[][], //constrained position for heavy atom
-		float4 cposq1[][], //ditto for h1
-		float4 cposq2[][], //ditto for h2
-		float4 cposq3[][], //ditto for h3
-		out float4 oposq<> //updated positions
-		){
-	float2 atom;
-	
-	atom.y = floor( invmap.x / strwidth );
-	atom.x = invmap.x - atom.y * strwidth;
-	
-	if ( invmap.y < 0 )
-		oposq = posq;
-	else if ( invmap.y < 0.5f ) 
-		oposq = cposq0[ atom ];
-	else if ( invmap.y < 1.5f )
-		oposq = cposq1[ atom ];
-	else if ( invmap.y < 2.5f )
-		oposq = cposq2[ atom ];
-	else if ( invmap.y < 3.5f )
-		oposq = cposq3[ atom ];
-} */
-
 /* Fix for precision of terms first order in dt
  * To be used with corresponding update kernels
  * The input posqp is now changes to posq rather than
@@ -216,7 +53,6 @@ kshakeh_fix1(
 		float nit,      //number of iterations
 		float strwidth, //stream width of posq
 		float invmH,    //inverse mass of hydrogen
-		float omega,    //parameter from gromacs, possibly unused
 		float4 atoms<>, //heavy0, h1, h2, h3
 		float3 posq[][], //positions before update
 		float3 posqp[][], //changes to positions
@@ -329,7 +165,7 @@ kshakeh_fix1(
 		if ( diff < 1.0f )
 		  acor = 0.0f;
 		 else
-			acor = omega * ( ld1 - 2 * rrpr - rpsqij ) * params.y / ( rrpr + rij1sq ) ;
+			acor = ( ld1 - 2 * rrpr - rpsqij ) * params.y / ( rrpr + rij1sq ) ;
 		dr   = rij1 * acor;
 		xpi  += dr * params.x;
 		xpj1 -= dr * invmH;
@@ -343,7 +179,7 @@ kshakeh_fix1(
 		if ( diff < 1.0f )
 			acor = 0.0f;
 		else
-			acor = mask2 * omega * ( ld2 - 2 * rrpr - rpsqij ) * params.y / ( rrpr + rij2sq ) ;
+			acor = mask2 * ( ld2 - 2 * rrpr - rpsqij ) * params.y / ( rrpr + rij2sq ) ;
 		dr   = rij2 * acor;
 		xpi  += dr * params.x;
 		xpj2 -= dr * invmH;
@@ -357,7 +193,7 @@ kshakeh_fix1(
 		if ( diff < 1.0f )
 			 acor = 0.0f;
 		else
-			acor = mask3 * omega * ( ld3 - 2 * rrpr - rpsqij ) * params.y / ( rrpr + rij3sq ) ;
+			acor = mask3 * ( ld3 - 2 * rrpr - rpsqij ) * params.y / ( rrpr + rij3sq ) ;
 		dr   = rij3 * acor;
 		xpi  += dr * params.x;
 		xpj3 -= dr * invmH;
@@ -378,7 +214,6 @@ kshakeh_fix2(
 		float nit,      //number of iterations
 		float strwidth, //stream width of posq
 		float invmH,    //inverse mass of hydrogen
-		float omega,    //parameter from gromacs, possibly unused
 		float4 atoms<>, //heavy0, h1, h2, h3
 		float3 posq[][], //positions before update
 		float3 posqp[][], //changes to positions
@@ -492,7 +327,7 @@ kshakeh_fix2(
 		if ( diff < 1.0f )
 		  acor = 0.0f;
 		 else
-			acor = omega * ( ld1 - 2 * rrpr - rpsqij ) * params.y / ( rrpr + rij1sq ) ;
+			acor = ( ld1 - 2 * rrpr - rpsqij ) * params.y / ( rrpr + rij1sq ) ;
 		dr   = rij1 * acor;
 		xpi  += dr * params.x;
 		xpj1 -= dr * invmH;
@@ -506,7 +341,7 @@ kshakeh_fix2(
 		if ( diff < 1.0f )
 			acor = 0.0f;
 		else
-			acor = mask2 * omega * ( ld2 - 2 * rrpr - rpsqij ) * params.y / ( rrpr + rij2sq ) ;
+			acor = mask2 * ( ld2 - 2 * rrpr - rpsqij ) * params.y / ( rrpr + rij2sq ) ;
 		dr   = rij2 * acor;
 		xpi  += dr * params.x;
 		xpj2 -= dr * invmH;
@@ -520,7 +355,7 @@ kshakeh_fix2(
 		if ( diff < 1.0f )
 			 acor = 0.0f;
 		else
-			acor = mask3 * omega * ( ld3 - 2 * rrpr - rpsqij ) * params.y / ( rrpr + rij3sq ) ;
+			acor = mask3 * ( ld3 - 2 * rrpr - rpsqij ) * params.y / ( rrpr + rij3sq ) ;
 		dr   = rij3 * acor;
 		xpi  += dr * params.x;
 		xpj3 -= dr * invmH;
@@ -536,36 +371,6 @@ kshakeh_fix2(
 
 }
 
-/*Applies the updated delta's
- *Each atom occurs at one and only one position
- * */
-kernel void kshakeh_update1_fix1Old( 
-		float strwidth, //width of cposq streams
-		float2 invmap<>,    //shakeh inverse map
-		float3 posq<>,     //constrained deltas
-		float3 cposq0[][], //constrained delta for heavy atom
-		float3 cposq1[][], //ditto for h1
-		float3 cposq2[][], //ditto for h2
-		float3 cposq3[][], //ditto for h3
-		out float3 oposq<> //updated deltas
-		){
-	float2 atom;
-	
-	//atom.y = floor( invmap.x / strwidth );
-   atom.y  = round( (invmap.x - fmod( invmap.x, strwidth ))/strwidth );
-	atom.x = invmap.x - atom.y * strwidth;
-	
-	if ( invmap.y < 0 )
-		oposq = posq;
-	else if ( invmap.y < 0.5f ) 
-		oposq = cposq0[ atom ];
-	else if ( invmap.y < 1.5f )
-		oposq = cposq1[ atom ];
-	else if ( invmap.y < 2.5f )
-		oposq = cposq2[ atom ];
-	else if ( invmap.y < 3.5f )
-		oposq = cposq3[ atom ];
-}
 
 kernel void kshakeh_update2_fix1( 
 		float strwidth, //width of cposq streams

platforms/brook/src/gpu/kshakeh.h

@@ -32,46 +32,6 @@
 void  kshakeh_fix1 (const float  nit,
 		const float  strwidth,
 		const float  invmH,
-		const float  omega,
-		::brook::stream atoms,
-		::brook::stream posq,
-		::brook::stream posqp,
-		::brook::stream params,
-		::brook::stream cposq0,
-		::brook::stream cposq1,
-		::brook::stream cposq2,
-		::brook::stream cposq3);
-
-
-void  kshakeh_fix2 (const float  nit,
-		const float  strwidth,
-		const float  invmH,
-		const float  omega,
-		::brook::stream atoms,
-		::brook::stream posq,
-		::brook::stream posqp,
-		::brook::stream params,
-		::brook::stream cposq0,
-		::brook::stream cposq1,
-		::brook::stream cposq2,
-		::brook::stream cposq3);
-
-
-
-void  kshakeh_update (const float  strwidth,
-		::brook::stream invmap,
-		::brook::stream posq,
-		::brook::stream cposq0,
-		::brook::stream cposq1,
-		::brook::stream cposq2,
-		::brook::stream cposq3,
-		::brook::stream oposq) ;
-
-
-void  kshakeh (const float  nit,
-		const float  strwidth,
-		const float  invmH,
-		const float  omega,
 		::brook::stream atoms,
 		::brook::stream posq,
 		::brook::stream posqp,
@@ -94,14 +54,6 @@ void  kshakeh_update1_fix1 (
 		::brook::stream cposq3,
 		::brook::stream oposq); 
 
-void  kshakeh_update1_fix1Old (const float  strwidth,
-		::brook::stream invmap,
-		::brook::stream posq,
-		::brook::stream cposq0,
-		::brook::stream cposq1,
-		::brook::stream cposq2,
-		::brook::stream cposq3,
-		::brook::stream oposq); 
 
 void  kshakeh_update2_fix1 (const float  strwidth,
 		::brook::stream invmap,