Added functionality to BrookRemoveCMMotionKernel' fixed problem w/ BrookVelocityCenterOfMassRemoval

platforms/brook/src/BrookIntegrateLangevinStepKernel.cpp

@@ -144,7 +144,7 @@ void BrookIntegrateLangevinStepKernel::execute( Stream& positions, Stream& veloc
    differences[1] = friction    - (double) _brookStochasticDynamics->getFriction();
    differences[2] = stepSize    - (double) _brookStochasticDynamics->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() );
       _brookStochasticDynamics->updateParameters( temperature, friction, stepSize );
    } else {
 //printf( "%s NOT calling updateParameters\n", methodName.c_str() );

platforms/brook/src/BrookKernelFactory.cpp

@@ -35,6 +35,7 @@
 #include "BrookIntegrateVerletStepKernel.h"
 #include "BrookCalcKineticEnergyKernel.h"
 #include "BrookCalcGBSAOBCForceFieldKernel.h"
+#include "BrookRemoveCMMotionKernel.h"
 
 using namespace OpenMM;
 
@@ -82,6 +83,12 @@ KernelImpl* BrookKernelFactory::createKernelImpl( std::string name, const Platfo
 
       return new BrookIntegrateLangevinStepKernel( name, platform );
 
+   // Remove com
+
+	} else if( name == RemoveCMMotionKernel::Name() ){
+
+      return new BrookRemoveCMMotionKernel( name, platform );
+
    // KE calculator
 
    } else if( name == CalcKineticEnergyKernel::Name() ){

platforms/brook/src/BrookPlatform.cpp

@@ -160,6 +160,7 @@ void BrookPlatform::_initializeKernelFactory( void ){
    //registerKernelFactory( IntegrateBrownianStepKernel::Name(),    factory);
    //registerKernelFactory( ApplyAndersenThermostatKernel::Name(),  factory);
    registerKernelFactory( CalcKineticEnergyKernel::Name(),        factory);
+   registerKernelFactory( RemoveCMMotionKernel::Name(),           factory);
    
 }
 

platforms/brook/src/BrookRemoveCMMotionKernel.cpp

@@ -52,6 +52,8 @@ BrookRemoveCMMotionKernel::BrookRemoveCMMotionKernel( std::string name, const Pl
 
 // ---------------------------------------------------------------------------------------
 
+   _brookVelocityCenterOfMassRemoval = NULL;
+
 }
 
 /** 
@@ -67,6 +69,7 @@ BrookRemoveCMMotionKernel::~BrookRemoveCMMotionKernel( ){
 
 // ---------------------------------------------------------------------------------------
 
+   delete _brookVelocityCenterOfMassRemoval;
 }
 
 /** 
@@ -84,11 +87,8 @@ void BrookRemoveCMMotionKernel::initialize( const vector<double>& masses ){
 
 // ---------------------------------------------------------------------------------------
 
-/*
-    this->masses.resize(masses.size());
-    for (size_t i = 0; i < masses.size(); ++i)
-        this->masses[i] = masses[i];
-*/
+   _brookVelocityCenterOfMassRemoval = new BrookVelocityCenterOfMassRemoval();
+   _brookVelocityCenterOfMassRemoval->setup( masses, getPlatform() );
 
    return;
 
@@ -109,28 +109,7 @@ void BrookRemoveCMMotionKernel::execute( Stream& velocities ){
 
 // ---------------------------------------------------------------------------------------
 
-/*
-    RealOpenMM** velData = ((BrookFloatStreamImpl&) velocities.getImpl()).getData();
-    
-    // Calculate the center of mass momentum.
-    
-    RealOpenMM momentum[] = {0.0, 0.0, 0.0};
-    for (size_t i = 0; i < masses.size(); ++i) {
-        momentum[0] += static_cast<RealOpenMM>( masses[i]*velData[i][0] );
-        momentum[1] += static_cast<RealOpenMM>( masses[i]*velData[i][1] );
-        momentum[2] += static_cast<RealOpenMM>( masses[i]*velData[i][2] );
-    }
-    
-    // Adjust the atom velocities.
-    
-    momentum[0] /= static_cast<RealOpenMM>( masses.size() );
-    momentum[1] /= static_cast<RealOpenMM>( masses.size() );
-    momentum[2] /= static_cast<RealOpenMM>( masses.size() );
-    for (size_t i = 0; i < masses.size(); ++i) {
-        velData[i][0] -= static_cast<RealOpenMM>( momentum[0]/masses[i] );
-        velData[i][1] -= static_cast<RealOpenMM>( momentum[1]/masses[i] );
-        velData[i][2] -= static_cast<RealOpenMM>( momentum[2]/masses[i] );
-    }
-*/
+   _brookVelocityCenterOfMassRemoval->removeVelocityCenterOfMass( velocities );
+
    return;
 }

platforms/brook/src/BrookRemoveCMMotionKernel.h

@@ -33,6 +33,7 @@
  * -------------------------------------------------------------------------- */
 
 #include "kernels.h"
+#include "BrookVelocityCenterOfMassRemoval.h"
 
 namespace OpenMM {
 
@@ -80,6 +81,7 @@ class BrookRemoveCMMotionKernel : public RemoveCMMotionKernel {
 
    private:
 
+      BrookVelocityCenterOfMassRemoval* _brookVelocityCenterOfMassRemoval;
 };
 
 } // namespace OpenMM

platforms/brook/src/BrookShakeAlgorithm.cpp

@@ -331,19 +331,19 @@ int BrookShakeAlgorithm::_initializeStreams( const Platform& platform ){
 
     _shakeStreams[ShakeXCons0Stream]         = new BrookFloatStreamInternal( BrookCommon::ShakeXCons0Stream,
                                                                              shakeConstraintStreamSize, shakeConstraintStreamWidth,
-                                                                             BrookStreamInternal::Float4, dangleValue );
+                                                                             BrookStreamInternal::Float3, dangleValue );
 
     _shakeStreams[ShakeXCons1Stream]         = new BrookFloatStreamInternal( BrookCommon::ShakeXCons1Stream,
                                                                              shakeConstraintStreamSize, shakeConstraintStreamWidth,
-                                                                             BrookStreamInternal::Float4, dangleValue );
+                                                                             BrookStreamInternal::Float3, dangleValue );
 
     _shakeStreams[ShakeXCons2Stream]         = new BrookFloatStreamInternal( BrookCommon::ShakeXCons2Stream,
                                                                              shakeConstraintStreamSize, shakeConstraintStreamWidth,
-                                                                             BrookStreamInternal::Float4, dangleValue );
+                                                                             BrookStreamInternal::Float3, dangleValue );
 
     _shakeStreams[ShakeXCons3Stream]         = new BrookFloatStreamInternal( BrookCommon::ShakeXCons3Stream,
                                                                              shakeConstraintStreamSize, shakeConstraintStreamWidth,
-                                                                             BrookStreamInternal::Float4, dangleValue );
+                                                                             BrookStreamInternal::Float3, dangleValue );
 
     _shakeStreams[ShakeInverseMapStream]     = new BrookFloatStreamInternal( BrookCommon::ShakeInverseMapStream,
                                                                              shakeAtomStreamSize, shakeAtomStreamWidth,
@@ -402,6 +402,7 @@ int BrookShakeAlgorithm::_setShakeStreams( const std::vector<double>& masses, co
    std::vector<double>::const_iterator distanceIterator         = constraintLengths.begin();
 
    int constraintIndex                = 0;
+   _numberOfConstraints               = 0;
    while( atomIterator != constraintIndices.end() ){
 
       std::vector<int> atomVector = *atomIterator;
@@ -434,6 +435,21 @@ int BrookShakeAlgorithm::_setShakeStreams( const std::vector<double>& masses, co
          }
       }
 
+      // skip negative indices
+
+      if( atomIndex1 < 0 || atomIndex2 < 0 ){
+         atomIterator++;
+         distanceIterator++;
+         continue;
+      }
+
+      if( atomIndex1 >= (int) masses.size() || atomIndex2 >= (int) masses.size() ){
+         std::stringstream message;
+         message << methodName << " constraint indices=[" << atomIndex1 << ", " << atomIndex2 << "] greater than mass array size=" << masses.size();
+         throw OpenMMException( message.str() );
+         return ErrorReturnValue;
+      }   
+
       // insure heavy atom is first
 
       if( masses[atomIndex1] < masses[atomIndex2] ){
@@ -460,6 +476,7 @@ int BrookShakeAlgorithm::_setShakeStreams( const std::vector<double>& masses, co
       atomIterator++;
       distanceIterator++;
       constraintIndex += 4;
+      _numberOfConstraints++;
    }
 
    // write entries to board

platforms/brook/src/BrookStochasticDynamics.cpp

@@ -91,8 +91,6 @@ BrookStochasticDynamics::BrookStochasticDynamics( ){
 
    _randomNumberSeed  = 1393;
 
-   _brookVelocityCenterOfMassRemoval = NULL;
-
    //_randomNumberSeed = randomNumberSeed ? randomNumberSeed : 1393;
    //SimTKOpenMMUtilities::setRandomNumberSeed( randomNumberSeed );
 }   
@@ -116,8 +114,6 @@ BrookStochasticDynamics::~BrookStochasticDynamics( ){
 
    delete[] _inverseSqrtMasses;
 
-   delete _brookVelocityCenterOfMassRemoval;
-
 }
 
 /** 
@@ -344,7 +340,7 @@ int BrookStochasticDynamics::updateParameters( double temperature, double fricti
 
    if( showUpdate && getLog() && (showUpdate++ < maxShowUpdate) ){
       std::string contents = getContentsString( );
-      (void) fprintf( getLog(), "%s contents %s", methodName, contents.c_str() );
+      (void) fprintf( getLog(), "%s contents\n%s", methodName, contents.c_str() );
       (void) fflush( getLog() );
 
    }
@@ -403,11 +399,14 @@ int BrookStochasticDynamics::update( Stream& positions, Stream& velocities,
    
       if( showAux ){
          showAux = 0;
-         std::string contents = _brookVelocityCenterOfMassRemoval->getContentsString( );
-         (void) fprintf( getLog(), "%s VelocityCenterOfMassRemoval contents\n%s", methodName, contents.c_str() );
-         contents             = brookRandomNumberGenerator.getContentsString( );
+
+         std::string contents             = brookRandomNumberGenerator.getContentsString( );
          (void) fprintf( getLog(), "%s RNG contents\n%s", methodName, contents.c_str() );
-         brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->printToFile( getLog() ); 
+         // brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->printToFile( getLog() ); 
+
+         contents             = brookShakeAlgorithm.getContentsString( );
+         (void) fprintf( getLog(), "%s Shake contents\n%s", methodName, contents.c_str() );
+
          (void) fflush( getLog() );
       }
    }
@@ -612,8 +611,6 @@ int BrookStochasticDynamics::update( Stream& positions, Stream& velocities,
       ksetStr3( getXPrimeStream()->getBrookStream(), positionStream.getBrookStream() );
    }
 
-   _brookVelocityCenterOfMassRemoval->removeVelocityCenterOfMass( velocities );
-
    return DefaultReturnValue;
 
 };
@@ -1059,9 +1056,6 @@ int BrookStochasticDynamics::setup( const std::vector<double>& masses, const Pla
 
    _setInverseSqrtMasses( masses );
 
-   _brookVelocityCenterOfMassRemoval = new BrookVelocityCenterOfMassRemoval();
-   _brookVelocityCenterOfMassRemoval->setup( masses, platform );
-
    return DefaultReturnValue;
 }
 

platforms/brook/src/BrookStochasticDynamics.h

@@ -38,7 +38,6 @@
 #include "BrookFloatStreamInternal.h"
 #include "BrookShakeAlgorithm.h"
 #include "BrookRandomNumberGenerator.h"
-#include "BrookVelocityCenterOfMassRemoval.h"
 #include "BrookPlatform.h"
 #include "BrookCommon.h"
 
@@ -324,10 +323,6 @@ class BrookStochasticDynamics : public BrookCommon {
 
       BrookOpenMMFloat* _inverseSqrtMasses;
 
-      // remove com
-
-      BrookVelocityCenterOfMassRemoval* _brookVelocityCenterOfMassRemoval;
-
       // internal streams
 
       BrookFloatStreamInternal* _sdStreams[LastStreamIndex];

platforms/brook/src/BrookVelocityCenterOfMassRemoval.cpp

@@ -147,29 +147,102 @@ int BrookVelocityCenterOfMassRemoval::removeVelocityCenterOfMass( Stream& veloci
 
    // ---------------------------------------------------------------------------------------
 
-   // static const char* methodName  = "\nBrookVelocityCenterOfMassRemoval::removeVelocityCenterOfMass";
+   static const char* methodName  = "BrookVelocityCenterOfMassRemoval::removeVelocityCenterOfMass";
 
    static int debug               = 0;
-
-   // static char* testName[2] = { "kCalculateLinearMomentum", "kRemoveLinearMomentum" };
-   // char fileName[128];
-
    float zero                     = 0.0f;
-   float one                = 1.0f;
 
    // ---------------------------------------------------------------------------------------
 
+   if( debug && getLog() ){
+      BrookOpenMMFloat com[3];
+      getVelocityCenterOfMass( velocities, com );
+      (void) fprintf( getLog(), "\n%s Pre removal com: [%12.5e %12.5e %12.5e]\n", methodName, com[0], com[1], com[2] );
+   }
+
    // calculate linear momentum via reduction
    // subtract it (/totalMass) from velocities
 
    BrookStreamImpl& velocityStream  = dynamic_cast<BrookStreamImpl&> (velocities.getImpl());
 
    kCalculateLinearMomentum( getMassStream()->getBrookStream(), velocityStream.getBrookStream(), getWorkStream()->getBrookStream() );
-   kScale( zero, getLinearMomentumStream()->getBrookStream(), getLinearMomentumStream()->getBrookStream() );
    kSumLinearMomentum( (float) getComAtomStreamWidth(), (float) getNumberOfAtoms(), getWorkStream()->getBrookStream(), getLinearMomentumStream()->getBrookStream() );
    kScale( (float) getTotalInverseMass(), getLinearMomentumStream()->getBrookStream(), getLinearMomentumStream()->getBrookStream() );
    kRemoveLinearMomentum( getLinearMomentumStream()->getBrookStream(), velocityStream.getBrookStream(), velocityStream.getBrookStream() );
 
+   if( (0 || debug) && getLog() ){
+      BrookOpenMMFloat com[3];
+      getVelocityCenterOfMass( velocities, com );
+      (void) fprintf( getLog(), "\n%s strW=%d iatm=%d Post removal com: [%12.5e %12.5e %12.5e]\n", methodName,
+                      getComAtomStreamWidth(), getNumberOfAtoms(),  com[0], com[1], com[2] );
+/*
+      (void) fprintf( getLog(), "LM [%12.5e %12.5e %12.5e]\n", com[0], com[1], com[2] );
+      BrookStreamImpl& v1                  = dynamic_cast<BrookStreamImpl&> (velocities.getImpl());
+      BrookStreamInternal* velocityStream  = dynamic_cast<BrookStreamInternal*> (v1.getBrookStreamImpl());
+
+      void* velV                       = velocityStream->getData( 1 );
+      const float* vArray              = (float*) velV;
+
+      void* w1                        = getWorkStream()->getData( 1 );
+      const float* w2                 = (float*) w1;
+
+      int index                        = 0;
+      for( int ii = 0; ii < getNumberOfAtoms(); ii++, index += 3 ){
+         (void) fprintf( getLog(), "V %d [%12.5e %12.5e %12.5e] [%12.5e %12.5e %12.5e]\n", ii,
+                         vArray[index], vArray[index+1], vArray[index+2], w2[index], w2[index+1], w2[index+2] );
+      }
+      (void) fflush( getLog() );
+*/
+//exit(0);
+   }
+
+   return DefaultReturnValue;
+}
+
+/**
+ * Get velocity-COM
+ *
+ * @param velocities velocities
+ * @param  velocityCom                 output velocity com
+ *
+ * @return DefaultReturnValue
+ *
+ */
+   
+int BrookVelocityCenterOfMassRemoval::getVelocityCenterOfMass( Stream& velocities, BrookOpenMMFloat velocityCom[3] ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName  = "\nBrookVelocityCenterOfMassRemoval::getVelocityCenterOfMass";
+
+   static int debug         = 0;
+   BrookOpenMMFloat zero    = (BrookOpenMMFloat) 0.0;
+
+   // ---------------------------------------------------------------------------------------
+
+   // calculate linear momentum via reduction
+   // subtract it (/totalMass) from velocities
+
+   BrookStreamImpl& v1                  = dynamic_cast<BrookStreamImpl&> (velocities.getImpl());
+   BrookStreamInternal* velocityStream  = dynamic_cast<BrookStreamInternal*> (v1.getBrookStreamImpl());
+
+   void* velV                           = velocityStream->getData( 1 );
+   const float* vArray                  = (float*) velV;
+
+   void* massV                          = getMassStream()->getData( 1);
+   const float* mArray                  = (float*) massV;
+
+   int numberOfAtoms                    = getNumberOfAtoms();
+   int index                            = 0;
+
+   velocityCom[0] = velocityCom[1] = velocityCom[2] = zero;
+
+   for( int ii = 0; ii < numberOfAtoms; ii++, index += 3 ){
+      velocityCom[0] += mArray[ii]*vArray[index];
+      velocityCom[1] += mArray[ii]*vArray[index+1];
+      velocityCom[2] += mArray[ii]*vArray[index+2];
+   }
+
    return DefaultReturnValue;
 }
 
@@ -254,7 +327,7 @@ int BrookVelocityCenterOfMassRemoval::_initializeStreams( const Platform& platfo
 
     _streams[WorkStream]            = new BrookFloatStreamInternal( BrookCommon::BrookVelocityCenterOfMassRemovalWorkStream,
                                                                     atomStreamSize, atomStreamWidth,
-                                                                    BrookStreamInternal::Float, dangleValue );
+                                                                    BrookStreamInternal::Float3, dangleValue );
 
 
     _streams[MassStream]            = new BrookFloatStreamInternal( BrookCommon::BrookVelocityCenterOfMassRemovalMassStream,
@@ -263,14 +336,14 @@ int BrookVelocityCenterOfMassRemoval::_initializeStreams( const Platform& platfo
 
 
     _streams[LinearMomentumStream]  = new BrookFloatStreamInternal( "LinearMomentumStream",
-                                                                    1, 3, BrookStreamInternal::Float, dangleValue );
+                                                                    1, 1, BrookStreamInternal::Float3, dangleValue );
 
 
    return DefaultReturnValue;
 }
 
 /** 
- * Set inverse masses 
+ * Set masses & _totalInverseMass 
  * 
  * @param masses             atomic masses
  *
@@ -287,36 +360,40 @@ int BrookVelocityCenterOfMassRemoval::_setMasses( const std::vector<double>& mas
 
 // ---------------------------------------------------------------------------------------
 
+   // check that masses vector is not larger than expected
+
    if( (int) masses.size() > getComAtomStreamSize() ){
       std::stringstream message;
       message << methodName << " mass array size=" << masses.size() << " is larger than stream size=" << getComAtomStreamSize();
       throw OpenMMException( message.str() );
    }
 
+   // setup masses
 
-   // setup inverse sqrt masses
-
-   BrookOpenMMFloat* inverseMasses = new BrookOpenMMFloat[getComAtomStreamSize()];
-   memset( inverseMasses, 0, sizeof( BrookOpenMMFloat )*getComAtomStreamSize() );
+   BrookOpenMMFloat* localMasses = new BrookOpenMMFloat[getComAtomStreamSize()];
+   memset( localMasses, 0, sizeof( BrookOpenMMFloat )*getComAtomStreamSize() );
 
    int index          = 0;
    _totalInverseMass  = (BrookOpenMMFloat) 0.0;
    for( std::vector<double>::const_iterator ii = masses.begin(); ii != masses.end(); ii++, index++ ){
       if( *ii != 0.0 ){
          BrookOpenMMFloat value    = static_cast<BrookOpenMMFloat>(*ii);
-         inverseMasses[index]      = value;
+         localMasses[index]        = value;
          _totalInverseMass        += value;
-      } else {
-         inverseMasses[index] = zero;
       }
    }
+
+   // 1/Sum[masses]
+
    if( _totalInverseMass > zero ){
       _totalInverseMass = one/_totalInverseMass;
    }
 
-   _streams[MassStream]->loadFromArray( inverseMasses );
+   // write masses to board
+
+   _streams[MassStream]->loadFromArray( localMasses );
 
-   delete[] inverseMasses;
+   delete[] localMasses;
 
    return DefaultReturnValue;
 }   

platforms/brook/src/BrookVelocityCenterOfMassRemoval.h

@@ -90,13 +90,9 @@ class BrookVelocityCenterOfMassRemoval : public BrookCommon {
       int getComAtomStreamSize( void ) const; 
 
       /** 
-       * Update
+       * Remove velocity center-of-mass
        * 
-       * @param  positions                   atom positions
        * @param  velocities                  atom velocities
-       * @param  forces                      atom forces
-       * @param  brookShakeAlgorithm         BrookShakeAlgorithm reference
-       * @param  brookRandomNumberGenerator  BrookRandomNumberGenerator reference
        *
        * @return  DefaultReturnValue
        *
@@ -104,6 +100,18 @@ class BrookVelocityCenterOfMassRemoval : public BrookCommon {
       
       int removeVelocityCenterOfMass( Stream& velocities );
 
+      /** 
+       * Get velocity center-of-mass
+       * 
+       * @param  velocities                  atom velocities
+       * @param  velocityCom                 output velocity com
+       *
+       * @return  DefaultReturnValue
+       *
+       */
+      
+      int getVelocityCenterOfMass( Stream& velocities, BrookOpenMMFloat velocityCom[3] );
+
       /* 
        * Setup of  parameters
        *

platforms/brook/src/BrookVerletDynamics.cpp

@@ -155,7 +155,7 @@ int BrookVerletDynamics::updateParameters( double stepSize ){
 
    if( showUpdate && getLog() && (showUpdate++ < maxShowUpdate) ){
       std::string contents = getContentsString( );
-      (void) fprintf( getLog(), "%s contents %s", methodName, contents.c_str() );
+      (void) fprintf( getLog(), "%s contents\n%s", methodName, contents.c_str() );
       (void) fflush( getLog() );
 
    }
@@ -189,7 +189,7 @@ int BrookVerletDynamics::update( Stream& positions, Stream& velocities,
 
    static const char* methodName = "\nBrookVerletDynamics::update";
 
-   static const int PrintOn      = 1;
+   static const int PrintOn      = 0;
 
 // ---------------------------------------------------------------------------------------
 
@@ -198,10 +198,26 @@ int BrookVerletDynamics::update( Stream& positions, Stream& velocities,
    const BrookStreamImpl& forceStreamC                 = dynamic_cast<const BrookStreamImpl&> (forces.getImpl());
    BrookStreamImpl& forceStream                        = const_cast<BrookStreamImpl&> (forceStreamC);
 
-   if( PrintOn && getLog() ){
+   if( (1 || PrintOn) && getLog() ){
+
+      static int showAux = 1;
+
       (void) fprintf( getLog(), "%s shake=%d\n", methodName, brookShakeAlgorithm.getNumberOfConstraints() );
       (void) fflush( getLog() );
 
+      if( showAux ){
+         showAux = 0; 
+
+/*
+         std::string contents = _brookVelocityCenterOfMassRemoval->getContentsString( );
+         (void) fprintf( getLog(), "%s VelocityCenterOfMassRemoval contents\n%s", methodName, contents.c_str() );
+*/
+         std::string contents   = brookShakeAlgorithm.getContentsString( );
+         (void) fprintf( getLog(), "%s Shake contents\n%s", methodName, contents.c_str() );
+
+         (void) fflush( getLog() );
+      }    
+
    }
 
    // integration step
@@ -243,6 +259,7 @@ int BrookVerletDynamics::update( Stream& positions, Stream& velocities,
    // second Shake step
 
    if( brookShakeAlgorithm.getNumberOfConstraints() > 0 ){
+
       kshakeh_fix1( 
                     10.0f,
                     (float) getVerletDynamicsAtomStreamWidth(),
@@ -257,6 +274,38 @@ int BrookVerletDynamics::update( Stream& positions, Stream& velocities,
                     brookShakeAlgorithm.getShakeXCons2Stream()->getBrookStream(),
                     brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream() );
    
+      if( (1|| PrintOn) && getLog() ){
+
+         (void) fprintf( getLog(), "\nPost kshakeh_fix1: atomStrW=%3d invMass_H=%.5f",
+                                   getVerletDynamicsAtomStreamWidth(), brookShakeAlgorithm.getInverseHydrogenMass() );
+   
+         (void) fprintf( getLog(), "\nShakeAtomIndicesStream\n" );
+         brookShakeAlgorithm.getShakeAtomIndicesStream()->printToFile( getLog() );
+   
+         BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamImpl();
+         (void) fprintf( getLog(), "\nPositionStream\n" );
+         brookStreamInternalPos->printToFile( getLog() );
+   
+         (void) fprintf( getLog(), "\nXPrimeStream\n" );
+         getXPrimeStream()->printToFile( getLog() ); 
+
+         (void) fprintf( getLog(), "\nShakeAtomParameterStream\n" );
+         brookShakeAlgorithm.getShakeAtomParameterStream()->printToFile( getLog() ); 
+
+         (void) fprintf( getLog(), "\nShakeXCons0\n" );
+         brookShakeAlgorithm.getShakeXCons0Stream()->printToFile( getLog() ); 
+
+         (void) fprintf( getLog(), "\nShakeXCons1\n" );
+         brookShakeAlgorithm.getShakeXCons1Stream()->printToFile( getLog() ); 
+
+         (void) fprintf( getLog(), "\nShakeXCons2\n" );
+         brookShakeAlgorithm.getShakeXCons2Stream()->printToFile( getLog() ); 
+
+         (void) fprintf( getLog(), "\nShakeXCons3\n" );
+         brookShakeAlgorithm.getShakeXCons3Stream()->printToFile( getLog() ); 
+
+      }   
+
       // second Shake gather
    
       kshakeh_update2_fix1( 
@@ -270,11 +319,42 @@ int BrookVerletDynamics::update( Stream& positions, Stream& velocities,
                     brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream(),
                     positionStream.getBrookStream() );
    
+      if( ( 1 || PrintOn) && getLog() ){
+
+         (void) fprintf( getLog(), "\nPost kshakeh_update2_fix1: atomStrW=%3d",
+                                   getVerletDynamicsAtomStreamWidth() );
+   
+         (void) fprintf( getLog(), "\nShakeInverseMapStream\n" );
+         brookShakeAlgorithm.getShakeInverseMapStream()->printToFile( getLog() );
+   
+         BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamImpl();
+         (void) fprintf( getLog(), "\nPositionStream\n" );
+         brookStreamInternalPos->printToFile( getLog() );
+   
+         (void) fprintf( getLog(), "\nXPrimeStream\n" );
+         getXPrimeStream()->printToFile( getLog() ); 
+
+         (void) fprintf( getLog(), "\nShakeXCons0\n" );
+         brookShakeAlgorithm.getShakeXCons0Stream()->printToFile( getLog() ); 
+
+         (void) fprintf( getLog(), "\nShakeXCons1\n" );
+         brookShakeAlgorithm.getShakeXCons1Stream()->printToFile( getLog() ); 
+
+         (void) fprintf( getLog(), "\nShakeXCons2\n" );
+         brookShakeAlgorithm.getShakeXCons2Stream()->printToFile( getLog() ); 
+
+         (void) fprintf( getLog(), "\nShakeXCons3\n" );
+         brookShakeAlgorithm.getShakeXCons3Stream()->printToFile( getLog() ); 
+
+      }   
+
    } else {
       //kadd3( getXPrimeStream()->getBrookStream(), positionStream.getBrookStream() );
       ksetStr3( getXPrimeStream()->getBrookStream(), positionStream.getBrookStream() );
    }
 
+//_brookVelocityCenterOfMassRemoval->removeVelocityCenterOfMass( velocities );
+
    return DefaultReturnValue;
 
 };

platforms/brook/src/gpu/kshakeh.br

@@ -193,13 +193,13 @@ kshakeh_fix1(
 		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[][], //changes to positions
+		float3 posq[][], //positions before update
+		float3 posqp[][], //changes to positions
 		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
+		out float3 cposq0<>, //constrained position for heavy atom
+		out float3 cposq1<>, //ditto for h1
+		out float3 cposq2<>, //ditto for h2
+		out float3 cposq3<> //ditto for h3
 		) {
 	float2 ai, aj1, aj2, aj3; //2d indices, can be precalc.
 	float i; //iteration count
@@ -246,16 +246,15 @@ kshakeh_fix1(
 		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;
+	xi  = cposq0;
+	xj1 = cposq1;
+	xj2 = cposq2;
+	xj3 = cposq3;
 	
 	rij1 = xi - xj1;
 	rij2 = xi - xj2;
@@ -270,17 +269,17 @@ kshakeh_fix1(
 	ld3 = params.z - rij3sq;
 
 /*
-	xpi  = posqp[ai].xyz;
-	xpj1 = posqp[aj1].xyz;
-	xpj2 = posqp[aj2].xyz;
-	xpj3 = posqp[aj3].xyz;
+	xpi  = posqp[ai];
+	xpj1 = posqp[aj1];
+	xpj2 = posqp[aj2];
+	xpj3 = posqp[aj3];
 */
 
 
-	xpi  = posqp[ai].xyz   - xi;
-	xpj1 = posqp[aj1].xyz  - xj1;
-	xpj2 = posqp[aj2].xyz  - xj2;
-	xpj3 = posqp[aj3].xyz  - xj3;
+	xpi  = posqp[ai]   - xi;
+	xpj1 = posqp[aj1]  - xj1;
+	xpj2 = posqp[aj2]  - xj2;
+	xpj3 = posqp[aj3]  - xj3;
 
 
 // XXX
@@ -342,10 +341,10 @@ kshakeh_fix1(
 	}
 
 	//Output modified delta's
-	cposq0.xyz = xpi;
-	cposq1.xyz = xpj1;
-	cposq2.xyz = xpj2;
-	cposq3.xyz = xpj3;
+	cposq0 = xpi;
+	cposq1 = xpj1;
+	cposq2 = xpj2;
+	cposq3 = xpj3;
 
 }
 
@@ -356,13 +355,13 @@ kshakeh_fix2(
 		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[][], //changes to positions
+		float3 posq[][], //positions before update
+		float3 posqp[][], //changes to positions
 		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
+		out float3 cposq0<>, //constrained position for heavy atom
+		out float3 cposq1<>, //ditto for h1
+		out float3 cposq2<>, //ditto for h2
+		out float3 cposq3<> //ditto for h3
 		) {
 	float2 ai, aj1, aj2, aj3; //2d indices, can be precalc.
 	float i; //iteration count
@@ -415,10 +414,10 @@ kshakeh_fix2(
 	cposq2 = posq[aj2];
 	cposq3 = posq[aj3];
 	
-	xi  = cposq0.xyz;
-	xj1 = cposq1.xyz;
-	xj2 = cposq2.xyz;
-	xj3 = cposq3.xyz;
+	xi  = cposq0;
+	xj1 = cposq1;
+	xj2 = cposq2;
+	xj3 = cposq3;
 	
 	rij1 = xi - xj1;
 	rij2 = xi - xj2;
@@ -432,17 +431,17 @@ kshakeh_fix2(
 	ld2 = params.z - rij2sq;
 	ld3 = params.z - rij3sq;
 
-	xpi  = posqp[ai].xyz;
-	xpj1 = posqp[aj1].xyz;
-	xpj2 = posqp[aj2].xyz;
-	xpj3 = posqp[aj3].xyz;
+	xpi  = posqp[ai];
+	xpj1 = posqp[aj1];
+	xpj2 = posqp[aj2];
+	xpj3 = posqp[aj3];
 
 
 /*
-	xpi  = posqp[ai].xyz   - xi;
-	xpj1 = posqp[aj1].xyz  - xj1;
-	xpj2 = posqp[aj2].xyz  - xj2;
-	xpj3 = posqp[aj3].xyz  - xj3;
+	xpi  = posqp[ai]   - xi;
+	xpj1 = posqp[aj1]  - xj1;
+	xpj2 = posqp[aj2]  - xj2;
+	xpj3 = posqp[aj3]  - xj3;
 
 */
 
@@ -505,10 +504,10 @@ kshakeh_fix2(
 	}
 
 	//Output modified delta's
-	cposq0.xyz = xpi;
-	cposq1.xyz = xpj1;
-	cposq2.xyz = xpj2;
-	cposq3.xyz = xpj3;
+	cposq0 = xpi;
+	cposq1 = xpj1;
+	cposq2 = xpj2;
+	cposq3 = xpj3;
 
 }
 
@@ -518,12 +517,12 @@ kshakeh_fix2(
 kernel void kshakeh_update1_fix1Old( 
 		float strwidth, //width of cposq streams
 		float2 invmap<>,    //shakeh inverse map
-		float4 posq<>,     //constrained deltas
-		float4 cposq0[][], //constrained delta for heavy atom
-		float4 cposq1[][], //ditto for h1
-		float4 cposq2[][], //ditto for h2
-		float4 cposq3[][], //ditto for h3
-		out float4 oposq<> //updated deltas
+		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;
 	
@@ -546,13 +545,13 @@ kernel void kshakeh_update1_fix1Old(
 kernel void kshakeh_update2_fix1( 
 		float strwidth, //width of cposq streams
 		float2 invmap<>,    //shakeh inverse map
-		float4 posq<>,     //old positions
-		float4 posqp<>,    //deltas from sd2
-		float4 cposq0[][], //constrained delta for heavy atom
-		float4 cposq1[][], //ditto for h1
-		float4 cposq2[][], //ditto for h2
-		float4 cposq3[][], //ditto for h3
-		out float4 oposq<> //updated deltas
+		float3 posq<>,     //old positions
+		float3 posqp<>,    //deltas from sd2
+		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;
 	
@@ -561,8 +560,7 @@ kernel void kshakeh_update2_fix1(
 	atom.x = invmap.x - atom.y * strwidth;
 	
 	if ( invmap.y < 0 ){
-		oposq.w = posq.w;
-		oposq.xyz = posqp.xyz - posq.xyz;
+		oposq = posqp - posq;
 	} else if ( invmap.y < 0.5f ) 
 		oposq = cposq0[ atom ];
 	else if ( invmap.y < 1.5f)
@@ -572,21 +570,21 @@ kernel void kshakeh_update2_fix1(
 	else if ( invmap.y < 3.5f )
 		oposq = cposq3[ atom ];
 
-	oposq.xyz += posq.xyz;
+	oposq += posq;
 }
 
 kernel void kshakeh_update1_fix1( 
 		float strwidth, //width of cposq streams
       float sdFactor,
 		float2 invmap<>,    //shakeh inverse map
-		float4 posq<>,     //old positions
-		float4 posqp<>,    //deltas from sd2
+		float3 posq<>,     //old positions
+		float3 posqp<>,    //deltas from sd2
 		float3 vp<>,    //deltas from sd2
-		float4 cposq0[][], //constrained delta for heavy atom
-		float4 cposq1[][], //ditto for h1
-		float4 cposq2[][], //ditto for h2
-		float4 cposq3[][], //ditto for h3
-		out float4 oposq<> //updated 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;
 	
@@ -598,12 +596,12 @@ kernel void kshakeh_update1_fix1(
 	if ( invmap.y < 0 ){
 		oposq = posqp;
 	} else if ( invmap.y < 0.5f ) 
-		oposq.xyz += cposq0[ atom ].xyz;
+		oposq += cposq0[ atom ];
 	else if ( invmap.y < 1.5f )
-		oposq.xyz += cposq1[ atom ].xyz;
+		oposq += cposq1[ atom ];
 	else if ( invmap.y < 2.5f )
-		oposq.xyz += cposq2[ atom ].xyz;
+		oposq += cposq2[ atom ];
 	else if ( invmap.y < 3.5f )
-		oposq.xyz += cposq3[ atom ].xyz;
+		oposq += cposq3[ atom ];
 
 }