Mods

platforms/brook/src/BrookBrownianDynamics.cpp

@@ -523,10 +523,12 @@ int BrookBrownianDynamics::update( Stream& positions, Stream& velocities,
    // Shake
 
    if( brookShakeAlgorithm.getNumberOfConstraints() > 0 ){
+/*
       kshakeh_fix1( 
-                    25.0f,
+                    (float) brookShakeAlgorithm.getMaxIterations(),
                     (float) getBrownianDynamicsParticleStreamWidth(),
                     brookShakeAlgorithm.getInverseHydrogenMass(),
+                    brookShakeAlgorithm.getShakeTolerance(),
                     brookShakeAlgorithm.getShakeParticleIndicesStream()->getBrookStream(),
                     positionStream.getBrookStream(),
                     getXPrimeStream()->getBrookStream(),
@@ -549,7 +551,8 @@ int BrookBrownianDynamics::update( Stream& positions, Stream& velocities,
                     brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream(),
                     positionStream.getBrookStream() );
    
-
+*/
+fprintf( stderr, "\nBrownian shake off!!\n" );
       // diagnostics
    
       if( PrintOn ){

platforms/brook/src/BrookIntegrateLangevinStepKernel.cpp

@@ -161,24 +161,21 @@ void BrookIntegrateLangevinStepKernel::initialize( const System& system, const L
 
    std::vector<std::vector<int> > constraintIndicesVector;
    constraintIndicesVector.resize( numberOfConstraints );
-
    std::vector<double> constraintLengths;
-   constraintLengths.resize( numberOfConstraints );
 
    for( int ii = 0; ii < numberOfConstraints; ii++ ){
 
       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 );
-      constraintIndicesVector[ii]  = constraintIndices;
 
-(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 );
+      constraintIndicesVector[ii].push_back( particle2 );
       constraintLengths.push_back( distance );
+
+//(void) fprintf( log, "%s shake setup const=%d ", methodName.c_str(), ii ); fflush( log );
+//(void) fprintf( log, "[ %d %d %f]\n", particle1, particle2, distance ); fflush( log );
    }
 
    _brookLangevinDynamics        = new BrookLangevinDynamics( );
@@ -187,6 +184,13 @@ void BrookIntegrateLangevinStepKernel::initialize( const System& system, const L
    _brookShakeAlgorithm          = new BrookShakeAlgorithm( );
    _brookShakeAlgorithm->setup( masses, constraintIndicesVector, constraintLengths, getPlatform() );
 
+   // tolerance
+
+   BrookOpenMMFloat tolerance = static_cast<BrookOpenMMFloat>( integrator.getConstraintTolerance() );
+   _brookShakeAlgorithm->setShakeTolerance( tolerance );
+   _brookShakeAlgorithm->setMaxIterations( 30 );
+   _brookShakeAlgorithm->setLog( log );
+
    // assert( (_brookShakeAlgorithm->getNumberOfConstraints() > 0) );
 
    if( printOn && log ){
@@ -211,7 +215,7 @@ void BrookIntegrateLangevinStepKernel::execute( OpenMMContextImpl& context, cons
 
 // ---------------------------------------------------------------------------------------
 
-   double epsilon                           = 1.0e-04;
+   double epsilon                           = 1.0e-06;
    static const std::string methodName      = "BrookIntegrateLangevinStepKernel::execute";
 
 // ---------------------------------------------------------------------------------------
@@ -228,11 +232,11 @@ 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()), 
                                    *(_openMMBrookInterface.getParticleForces()), *_brookShakeAlgorithm, *_brookRandomNumberGenerator );

platforms/brook/src/BrookLangevinDynamics.cpp

@@ -369,13 +369,9 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
 
 // ---------------------------------------------------------------------------------------
 
-   // unused Shake parameter
-
-   float numberOfIterations      = 25.0f;
-
    static const char* methodName = "\nBrookLangevinDynamics::update";
 
-   static const int PrintOn      = 1;
+   static const int PrintOn      = 0;
 
 // ---------------------------------------------------------------------------------------
 
@@ -482,9 +478,10 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
 
    if( brookShakeAlgorithm.getNumberOfConstraints() > 0 ){
       kshakeh_fix1( 
-                    numberOfIterations,
+                    (float) brookShakeAlgorithm.getMaxIterations(),
                     (float) getLangevinDynamicsParticleStreamWidth(),
                     brookShakeAlgorithm.getInverseHydrogenMass(),
+                    brookShakeAlgorithm.getShakeTolerance(),
                     brookShakeAlgorithm.getShakeParticleIndicesStream()->getBrookStream(),
                     positionStream.getBrookStream(),
                     getXPrimeStream()->getBrookStream(),
@@ -498,11 +495,9 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
    
       kshakeh_update1_fix1(
                     (float) getLangevinDynamicsParticleStreamWidth(),
-                    derivedParameters[Sd2pc1],
                     brookShakeAlgorithm.getShakeInverseMapStream()->getBrookStream(),
                     positionStream.getBrookStream(),
                     getXPrimeStream()->getBrookStream(),
-                    getVPrimeStream()->getBrookStream(),
                     brookShakeAlgorithm.getShakeXCons0Stream()->getBrookStream(),
                     brookShakeAlgorithm.getShakeXCons1Stream()->getBrookStream(),
                     brookShakeAlgorithm.getShakeXCons2Stream()->getBrookStream(),
@@ -574,9 +569,10 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
 
    if( brookShakeAlgorithm.getNumberOfConstraints() > 0 ){
       kshakeh_fix1( 
-                    numberOfIterations,
+                    (float) brookShakeAlgorithm.getMaxIterations(),
                     (float) getLangevinDynamicsParticleStreamWidth(),
                     brookShakeAlgorithm.getInverseHydrogenMass(),
+                    brookShakeAlgorithm.getShakeTolerance(),
                     brookShakeAlgorithm.getShakeParticleIndicesStream()->getBrookStream(),
                     positionStream.getBrookStream(),
                     getXPrimeStream()->getBrookStream(),

platforms/brook/src/BrookShakeAlgorithm.cpp

@@ -57,6 +57,8 @@ BrookShakeAlgorithm::BrookShakeAlgorithm( ){
 
    _numberOfParticles            = -1;
    _numberOfConstraints          = -1;
+   _maxIterations                = 25;
+   _shakeTolerance               = static_cast<BrookOpenMMFloat>(1.0e-04);
 
    // mark stream dimension variables as unset
 
@@ -112,6 +114,56 @@ int BrookShakeAlgorithm::getNumberOfConstraints( void ) const {
    return _numberOfConstraints;
 }
 
+/** 
+ * Get max iterations
+ * 
+ * @return  max iterations
+ *
+ */
+         
+int BrookShakeAlgorithm::getMaxIterations( void ) const {
+   return _maxIterations;
+}
+         
+/** 
+ * Set max iterations
+ * 
+ * @param  max iterations
+ *
+ * @return DefaultReturnValue
+ *
+ */
+         
+int BrookShakeAlgorithm::setMaxIterations( int maxIterations ){
+   _maxIterations = maxIterations;
+   return DefaultReturnValue;
+}
+    
+/** 
+ * Get SHAKE tolerance
+ * 
+ * @return  SHAKE tolerance  
+ *
+ */
+         
+BrookOpenMMFloat BrookShakeAlgorithm::getShakeTolerance( void ) const {
+   return _shakeTolerance;
+}
+         
+/** 
+ * Set SHAKE tolerance
+ * 
+ * @param  SHAKE tolerance  
+ *
+ * @return DefaultReturnValue
+ *
+ */
+         
+int BrookShakeAlgorithm::setShakeTolerance( BrookOpenMMFloat tolerance ){
+   _shakeTolerance = tolerance;
+   return DefaultReturnValue;
+}
+    
 /** 
  * Get inverse of hydrogen mass
  * 
@@ -585,6 +637,12 @@ std::string BrookShakeAlgorithm::getContentsString( int level ) const {
    (void) LOCAL_SPRINTF( value, "%d", getNumberOfParticles() );
    message << _getLine( tab, "Number of particles:", value ); 
 
+   (void) LOCAL_SPRINTF( value, "%d", getMaxIterations() );
+   message << _getLine( tab, "Max iterations:", value ); 
+
+   (void) LOCAL_SPRINTF( value, "%2e", getShakeTolerance() );
+   message << _getLine( tab, "Tolerance:", value ); 
+
    (void) LOCAL_SPRINTF( value, "%d", getParticleStreamWidth() );
    message << _getLine( tab, "Particle stream width:", value ); 
 

platforms/brook/src/BrookShakeAlgorithm.h

@@ -74,6 +74,46 @@ class BrookShakeAlgorithm : public BrookCommon {
       
       int getNumberOfConstraints( void ) const;
       
+      /** 
+       * Get max iterations
+       * 
+       * @return   max iterations
+       *
+       */
+      
+      int getMaxIterations( void ) const;
+      
+      /** 
+       * Set  max iterations
+       * 
+       * @param   max iterations
+       *
+       * @return DefaultReturnValue
+       *
+       */
+      
+      int setMaxIterations( int maxIterations );
+      
+      /** 
+       * Get SHAKE tolerance
+       * 
+       * @return  SHAKE tolerance  
+       *
+       */
+      
+      BrookOpenMMFloat getShakeTolerance( void ) const;
+      
+      /** 
+       * Set SHAKE tolerance
+       * 
+       * @param  SHAKE tolerance  
+       *
+       * @return DefaultReturnValue
+       *
+       */
+      
+      int setShakeTolerance( BrookOpenMMFloat tolerance );
+      
       /**
        * Get inverse of hydrogen mass
        *
@@ -247,6 +287,10 @@ class BrookShakeAlgorithm : public BrookCommon {
 
       int _numberOfConstraints;
 
+      // max iterations
+
+      int _maxIterations;
+
       // inverse of H mass
 
       BrookOpenMMFloat _inverseHydrogenMass;
@@ -263,6 +307,10 @@ class BrookShakeAlgorithm : public BrookCommon {
       int _shakeConstraintStreamWidth;
       int _shakeConstraintStreamHeight;
 
+      // SHAKE tolerance
+
+      BrookOpenMMFloat _shakeTolerance;
+
       // inverse sqrt masses
 
       BrookOpenMMFloat* _inverseSqrtMasses;

platforms/brook/src/BrookVerletDynamics.cpp

@@ -260,10 +260,12 @@ int BrookVerletDynamics::update( Stream& positions, Stream& velocities,
 
    if( brookShakeAlgorithm.getNumberOfConstraints() > 0 ){
 
+/*
       kshakeh_fix1( 
-                    25.0f,
+                    (float) brookShakeAlgorithm.getMaxIterations(),
                     (float) getVerletDynamicsParticleStreamWidth(),
                     brookShakeAlgorithm.getInverseHydrogenMass(),
+                    brookShakeAlgorithm.getShakeTolerance(),
                     brookShakeAlgorithm.getShakeParticleIndicesStream()->getBrookStream(),
                     positionStream.getBrookStream(),
                     getXPrimeStream()->getBrookStream(),
@@ -272,6 +274,8 @@ int BrookVerletDynamics::update( Stream& positions, Stream& velocities,
                     brookShakeAlgorithm.getShakeXCons1Stream()->getBrookStream(),
                     brookShakeAlgorithm.getShakeXCons2Stream()->getBrookStream(),
                     brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream() );
+*/
+fprintf( stderr, "\nVerlet shake off!!\n" );
    
       if( (1|| PrintOn) && getLog() ){
 

platforms/brook/src/OpenMMBrookInterface.cpp

@@ -551,13 +551,17 @@ void OpenMMBrookInterface::computeForces( OpenMMContextImpl& context ){
 
    // bonded forces
 
-(void) fprintf( stderr, "%s done nonbonded: bonded=%d completed=%d\n", methodName.c_str(), _brookBonded.isActive(), _brookBonded.isSetupCompleted() ); (void) fflush( stderr );
+   if( PrintOn && getLog() ){
+      (void) fprintf( getLog(), "%s done nonbonded: bonded=%d completed=%d\n", methodName.c_str(),
+                      _brookBonded.isActive(), _brookBonded.isSetupCompleted() ); 
+      (void) fflush( getLog() );
+   }
+
    if( _brookBonded.isActive() ){
 
       // perform setup first time through
 
-      //if( _brookBonded.isSetupCompleted() == 0 ){
-      if( _brookBonded.isSetupCompleted() >= -1 ){
+      if( _brookBonded.isSetupCompleted() == 0 ){
          _brookBonded.setup( getNumberOfParticles(), getHarmonicBondForceParameters(), getHarmonicAngleForceParameters(),
                              getPeriodicTorsionForceParameters(), getRBTorsionForceParameters(),
                              getNonBonded14ForceParameters(),