Mods

platforms/brook/src/BrookCalcKineticEnergyKernel.cpp

@@ -127,29 +127,33 @@ double BrookCalcKineticEnergyKernel::execute( OpenMMContextImpl& context ){
 
 // ---------------------------------------------------------------------------------------
 
-   void* dataV                            = _openMMBrookInterface.getParticleVelocities()->getData( );
-   float* velocity                        = (float*) dataV;
-
-   double energy                          = 0.0;
-   int index                              = 0;
-
    if( _masses == NULL ){
       std::stringstream message;
       message << methodName << " masses not set.";
       throw OpenMMException( message.str() );
    }    
 
-/*
-printf( "   BrookCalcKineticEnergyKernel Masses=%12.5e %12.5e", _masses[0], _masses[1] );
-printf( " [%12.5e %12.5e %12.5e]", velocity[index], velocity[index+1], velocity[index+2] );
-index += 3;
-printf( " [%12.5e %12.5e %12.5e]\n", velocity[index], velocity[index+1], velocity[index+2] );
+   void* dataV                            = _openMMBrookInterface.getParticleVelocities()->getData( 1 );
+   float* velocity                        = (float*) dataV;
+   double energy                          = 0.0;
+   int index                              = 0;
+
+printf( "BrookCalcKineticEnergyKernel:\n" );
+double com[3] = { 0.0, 0.0, 0.0 };
+for ( int ii = 0; ii < _numberOfParticles; ii++, index += 3 ){
+   com[0] += velocity[index];
+   com[1] += velocity[index+1];
+   com[2] += velocity[index+2];
+   printf( "  %d %.3f [%12.5e %12.5e %12.5e]\n", ii, _masses[ii], velocity[index], velocity[index+1], velocity[index+2] );
+}
+printf( "Com [%12.5e %12.5e %12.5e]\n", com[0], com[1], com[2] );
+
 index = 0;
-*/
 
    for ( int ii = 0; ii < _numberOfParticles; ii++, index += 3 ){
       energy += _masses[ii]*(velocity[index]*velocity[index] + velocity[index + 1]*velocity[index + 1] + velocity[index + 2]*velocity[index + 2]);
    }
 
+printf( " Ke=%12.5e\n", 0.5*energy );
    return 0.5*energy;
 }

platforms/brook/src/BrookStreamImpl.cpp

@@ -271,6 +271,17 @@ void* BrookStreamImpl::getData( void ){
    return _brookStreamInternal->getData( );
 }
 
+/** 
+ * Set every element of this stream to the same value.
+ *
+ * @param readFromBoard if set, read data from board 
+ *
+ * @return data array
+ */
+void* BrookStreamImpl::getData( int readFromBoard ){
+   return _brookStreamInternal->getData( readFromBoard );
+}
+
 /** 
  * Get Brook stream
  * 

platforms/brook/src/BrookStreamImpl.h

@@ -114,6 +114,7 @@ class BrookStreamImpl : public StreamImpl {
        * @return data array
        */
       void* getData( void );
+      void* getData( int readFromBoard );
 
       /**
        * Get Brook stream

platforms/brook/src/BrookVelocityCenterOfMassRemoval.cpp

@@ -174,15 +174,16 @@ int BrookVelocityCenterOfMassRemoval::removeVelocityCenterOfMass( BrookStreamImp
    // subtract it (/totalMass) from velocities
 
    kCalculateLinearMomentum( getMassStream()->getBrookStream(), velocityStream.getBrookStream(), getWorkStream()->getBrookStream() );
-   kSumLinearMomentum( (float) getComParticleStreamWidth(), (float) getNumberOfParticles(), getWorkStream()->getBrookStream(), getLinearMomentumStream()->getBrookStream() );
-   kScale( (float) getTotalInverseMass(), getLinearMomentumStream()->getBrookStream(), getLinearMomentumStream()->getBrookStream() );
+   kSumLinearMomentum( (float) getComParticleStreamWidth(), (float) getNumberOfParticles(), (float) getTotalInverseMass(),
+                        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( 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] );
+      (void) fprintf( getLog(), "%s strW=%d iatm=%d invMass=%.4e Post removal com: [%12.5e %12.5e %12.5e]", methodName,
+                      getComParticleStreamWidth(), getNumberOfParticles(), getTotalInverseMass(),  com[0], com[1], com[2] );
 
       void* linMoV = getLinearMomentumStream()->getData( 1 );
       float* linMo = (float*) linMoV;
@@ -227,7 +228,6 @@ int BrookVelocityCenterOfMassRemoval::getVelocityCenterOfMass( BrookStreamImpl&
 
    // static const char* methodName  = "\nBrookVelocityCenterOfMassRemoval::getVelocityCenterOfMass";
 
-   static int debug         = 0;
    BrookOpenMMFloat zero    = (BrookOpenMMFloat) 0.0;
 
    // ---------------------------------------------------------------------------------------
@@ -240,7 +240,7 @@ int BrookVelocityCenterOfMassRemoval::getVelocityCenterOfMass( BrookStreamImpl&
    void* velV                           = velocityStream->getData( 1 );
    const float* vArray                  = (float*) velV;
 
-   void* massV                          = getMassStream()->getData( 1);
+   void* massV                          = getMassStream()->getData( 1 );
    const float* mArray                  = (float*) massV;
 
    int numberOfParticles                = getNumberOfParticles();
@@ -440,6 +440,12 @@ int BrookVelocityCenterOfMassRemoval::setup( const std::vector<double>& masses,
 
    _setMasses( masses );
 
+   if( 1 && getLog() ){
+      std::string contents = getContentsString( 0 );
+      (void) fprintf( getLog(), "%s contents:\n%s\n", methodName.c_str(), contents.c_str() );
+      (void) fflush( getLog() );
+   }
+   
    return DefaultReturnValue;
 }
 

platforms/brook/src/gpu/kcom.br

@@ -85,7 +85,7 @@ reduce void kReduceLinearMomentum( float3 momentum<>, reduce float3 linearMoment
 
    --------------------------------------------------------------------------------------- */
 
-kernel void kSumLinearMomentum( float atomStrWidth, float numberOfAtoms, float3 momentum[][],
+kernel void kSumLinearMomentum( float atomStrWidth, float numberOfAtoms, float scale, float3 momentum[][],
                                 out float3 linearMomentum<> ){
 
    float atomCount;
@@ -108,6 +108,7 @@ kernel void kSumLinearMomentum( float atomStrWidth, float numberOfAtoms, float3
 
       atomCount += 1.0f;
    } 
+   linearMomentum *= scale;
 }
 
 /**---------------------------------------------------------------------------------------

platforms/brook/src/gpu/kcom.h

@@ -42,11 +42,12 @@ void kSum( ::brook::stream array, ::brook::stream sum );
 
    @param atomStrWidth   atom stream width
    @param numberOfAtoms  number of atoms
+   @param scale          sum of inverse masses
    @param momentum       momentum
    @param linearMomentum total momentum
 
    --------------------------------------------------------------------------------------- */
 
-void kSumLinearMomentum( float atomStrWidth, float numberOfAtoms, ::brook::stream momentum, ::brook::stream linearMomentum );
+void kSumLinearMomentum( float atomStrWidth, float numberOfAtoms, float scale, ::brook::stream momentum, ::brook::stream linearMomentum );