Skip to content

GitLab

Commit d5ab0e76 authored by Mark Friedrichs's avatar Mark Friedrichs
Browse files

Mods

parent 7a51ccf6
Hide whitespace changes
Inline Side-by-side
Showing with 60 additions and 28 deletions
platforms/brook/src/BrookKernelFactory.cpp
View file @ d5ab0e76
...@@ -38,7 +38,7 @@ ...@@ -38,7 +38,7 @@
#include "BrookCalcNonbondedForceKernel.h" #include "BrookCalcNonbondedForceKernel.h"
#include "BrookIntegrateLangevinStepKernel.h" #include "BrookIntegrateLangevinStepKernel.h"
#include "BrookIntegrateVerletStepKernel.h" #include "BrookIntegrateVerletStepKernel.h"
#include "BrookIntegrateBrownianStepKernel.h" //#include "BrookIntegrateBrownianStepKernel.h"
#include "BrookCalcKineticEnergyKernel.h" #include "BrookCalcKineticEnergyKernel.h"
#include "BrookCalcGBSAOBCForceKernel.h" #include "BrookCalcGBSAOBCForceKernel.h"
#include "BrookRemoveCMMotionKernel.h" #include "BrookRemoveCMMotionKernel.h"
......
platforms/brook/src/BrookShakeAlgorithm.cpp
View file @ d5ab0e76
...@@ -48,8 +48,8 @@ BrookShakeAlgorithm::BrookShakeAlgorithm( ){ ...@@ -48,8 +48,8 @@ BrookShakeAlgorithm::BrookShakeAlgorithm( ){
//static const std::string methodName = "BrookShakeAlgorithm::BrookShakeAlgorithm"; //static const std::string methodName = "BrookShakeAlgorithm::BrookShakeAlgorithm";
BrookOpenMMFloat zero = (BrookOpenMMFloat) 0.0; BrookOpenMMFloat zero = static_cast<BrookOpenMMFloat>( 0.0 );
BrookOpenMMFloat one = (BrookOpenMMFloat) 1.0; BrookOpenMMFloat one = static_cast<BrookOpenMMFloat>( 1.0 );
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -348,7 +348,7 @@ int BrookShakeAlgorithm::_initializeStreams( const Platform& platform ){ ...@@ -348,7 +348,7 @@ int BrookShakeAlgorithm::_initializeStreams( const Platform& platform ){
//static const std::string methodName = "BrookShakeAlgorithm::_initializeStreams"; //static const std::string methodName = "BrookShakeAlgorithm::_initializeStreams";
BrookOpenMMFloat dangleValue = (BrookOpenMMFloat) 0.0; BrookOpenMMFloat dangleValue = static_cast<BrookOpenMMFloat>( 0.0 );
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -405,19 +405,22 @@ int BrookShakeAlgorithm::_initializeStreams( const Platform& platform ){ ...@@ -405,19 +405,22 @@ int BrookShakeAlgorithm::_initializeStreams( const Platform& platform ){
* *
*/ */
int BrookShakeAlgorithm::_setShakeStreams( const std::vector<double>& masses, const std::vector< std::vector<int> >& constraintIndices, int BrookShakeAlgorithm::_setShakeStreams( const std::vector<double>& masses,
const std::vector<double>& constraintLengths, const Platform& platform ){ const std::vector< std::vector<int> >& constraintIndices,
const std::vector<double>& constraintLengths,
const Platform& platform ){
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
BrookOpenMMFloat one = (BrookOpenMMFloat) 1.0; BrookOpenMMFloat one = static_cast<BrookOpenMMFloat>( 1.0 );
BrookOpenMMFloat half = (BrookOpenMMFloat) 0.5; BrookOpenMMFloat half = static_cast<BrookOpenMMFloat>( 0.5 );
static const std::string methodName = "BrookShakeAlgorithm::_updateSdStreams"; static const std::string methodName = "BrookShakeAlgorithm::_updateSdStreams";
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
FILE* log = getLog(); FILE* log = getLog();
// check that number of constraints for two input vectors is consistent // check that number of constraints for two input vectors is consistent
if( constraintIndices.size() != constraintLengths.size() ){ if( constraintIndices.size() != constraintLengths.size() ){
...@@ -542,7 +545,7 @@ int BrookShakeAlgorithm::_setShakeStreams( const std::vector<double>& masses, co ...@@ -542,7 +545,7 @@ int BrookShakeAlgorithm::_setShakeStreams( const std::vector<double>& masses, co
} }
shakeParameters[constraintIndex] = static_cast<BrookOpenMMFloat>( cluster._centralInvMass ); shakeParameters[constraintIndex] = static_cast<BrookOpenMMFloat>( cluster._centralInvMass );
shakeParameters[constraintIndex+1] = half/( static_cast<BrookOpenMMFloat>( cluster._centralInvMass + cluster._peripheralInvMass) ); shakeParameters[constraintIndex+1] = half/( static_cast<BrookOpenMMFloat>( cluster._centralInvMass + cluster._peripheralInvMass ) );
shakeParameters[constraintIndex+2] = static_cast<BrookOpenMMFloat>( cluster._distance*cluster._distance ); shakeParameters[constraintIndex+2] = static_cast<BrookOpenMMFloat>( cluster._distance*cluster._distance );
shakeParameters[constraintIndex+3] = static_cast<BrookOpenMMFloat>( cluster._peripheralInvMass ); shakeParameters[constraintIndex+3] = static_cast<BrookOpenMMFloat>( cluster._peripheralInvMass );
......
platforms/brook/src/BrookVerletDynamics.cpp
View file @ d5ab0e76
...@@ -474,12 +474,12 @@ int BrookVerletDynamics::update( BrookStreamImpl& positionStream, BrookStreamImp ...@@ -474,12 +474,12 @@ int BrookVerletDynamics::update( BrookStreamImpl& positionStream, BrookStreamImp
_internalStepCount++; _internalStepCount++;
//setLog( stderr ); setLog( stderr );
printOn = (printOn && getLog()) ? printOn : 0; printOn = (printOn && getLog()) ? printOn : 0;
BrookStreamImpl& forceStream = const_cast<BrookStreamImpl&> (forceStreamC); BrookStreamImpl& forceStream = const_cast<BrookStreamImpl&> (forceStreamC);
if( printOn ){ if( 1 || printOn ){
static int showAux = 1; static int showAux = 1;
log = getLog(); log = getLog();
...@@ -490,7 +490,8 @@ int BrookVerletDynamics::update( BrookStreamImpl& positionStream, BrookStreamImp ...@@ -490,7 +490,8 @@ int BrookVerletDynamics::update( BrookStreamImpl& positionStream, BrookStreamImp
std::string contents = _brookVelocityCenterOfMassRemoval->getContentsString( ); std::string contents = _brookVelocityCenterOfMassRemoval->getContentsString( );
(void) fprintf( log, "%s VelocityCenterOfMassRemoval contents\n%s", methodName, contents.c_str() ); (void) fprintf( log, "%s VelocityCenterOfMassRemoval contents\n%s", methodName, contents.c_str() );
*/ */
(void) fprintf( log, "%s step=%d Shake contents\n%s", methodName.c_str(), _internalStepCount, brookShakeAlgorithm.getContentsString().c_str() ); (void) fprintf( log, "%s step=%d \n%s\n\nShake contents\n%s", methodName.c_str(), _internalStepCount, getContentsString().c_str(),
brookShakeAlgorithm.getContentsString().c_str() );
(void) fflush( log ); (void) fflush( log );
} }
...@@ -679,7 +680,7 @@ int BrookVerletDynamics::update( BrookStreamImpl& positionStream, BrookStreamImp ...@@ -679,7 +680,7 @@ int BrookVerletDynamics::update( BrookStreamImpl& positionStream, BrookStreamImp
// diagnostics // diagnostics
if( (_internalStepCount % 1) == 0 ){ if( (_internalStepCount % 10000) == 0 ){
FILE* log1 = stderr; FILE* log1 = stderr;
float epsilon = 1.0e-01f; float epsilon = 1.0e-01f;
......
platforms/brook/src/gpu/kshakeh.br
View file @ d5ab0e76
...@@ -135,8 +135,7 @@ kshakeh_fix1( ...@@ -135,8 +135,7 @@ kshakeh_fix1(
i = 0.0f; i = 0.0f;
converged = 1.0f; converged = 1.0f;
//tolerance = 2.0f*inputTolerance; tolerance = 2.0f*inputTolerance;
tolerance = inputTolerance;
while( i < maxIterations && converged > 0.0f ){ while( i < maxIterations && converged > 0.0f ){
//First hydrogen //First hydrogen
...@@ -144,8 +143,8 @@ kshakeh_fix1( ...@@ -144,8 +143,8 @@ kshakeh_fix1(
rpij = xpi - xpj1; //This is really rpij - rij rpij = xpi - xpj1; //This is really rpij - rij
rpsqij = dot( rpij, rpij ); //This is really deltar ^ 2 rpsqij = dot( rpij, rpij ); //This is really deltar ^ 2
rrpr = dot( rij1, rpij ); //This is r.deltar rrpr = dot( rij1, rpij ); //This is r.deltar
acor = ( ld1 - 2 * rrpr - rpsqij ) * params.y / ( rrpr + rij1sq ) ; acor = ( ld1 - 2.0f * rrpr - rpsqij ) * params.y / ( rrpr + rij1sq ) ;
diff = abs( ld1 - 2 * rrpr - rpsqij ) / (params.z * tolerance ); diff = abs( ld1 - 2.0f * rrpr - rpsqij ) / (params.z * tolerance );
acor = (diff < 1.0f) ? 0.0f : acor; acor = (diff < 1.0f) ? 0.0f : acor;
converged = abs( acor ); converged = abs( acor );
......
platforms/brook/src/gpu/kupdatesd.br
View file @ d5ab0e76
...@@ -109,17 +109,18 @@ kernel void kupdate_sd1_fix1( ...@@ -109,17 +109,18 @@ kernel void kupdate_sd1_fix1(
out float3 vnew<>, out float3 vnew<>,
out float3 posqp<> out float3 posqp<>
){ ){
float c1;
float3 Vmh; float3 Vmh;
float3 fg1, fg2; float3 fg1, fg2;
float linind_gauss; //linear index into the random numbers float linind_gauss; //linear index into the random numbers
float2 igauss; //2d index float2 igauss; //2d index
igauss = indexof( posq ); igauss = indexof( posq ).xy;
linind_gauss = igauss.y * xstrwidth + igauss.x; linind_gauss = igauss.y * xstrwidth + igauss.x;
//2 random vectors used in each fragment //2 random vectors used in each fragment
linind_gauss = 2 * linind_gauss + goffset; linind_gauss = 2.0f * linind_gauss + goffset;
igauss.y = round( (linind_gauss - fmod( linind_gauss, gstrwidth))/gstrwidth); igauss.y = round( (linind_gauss - fmod( linind_gauss, gstrwidth))/gstrwidth);
igauss.x = linind_gauss - igauss.y * gstrwidth; igauss.x = linind_gauss - igauss.y * gstrwidth;
fg1 = fgauss[ igauss ]; fg1 = fgauss[ igauss ];
...@@ -129,10 +130,27 @@ kernel void kupdate_sd1_fix1( ...@@ -129,10 +130,27 @@ kernel void kupdate_sd1_fix1(
igauss.x = linind_gauss - igauss.y * gstrwidth; igauss.x = linind_gauss - igauss.y * gstrwidth;
fg2 = fgauss[ igauss ]; fg2 = fgauss[ igauss ];
Vmh = (sd2X*pc3) + (sdpc.x*fg1); //Vmh = (sd2X*pc3) + (sdpc.x*fg1);
sd1V = sdpc.y*fg2;
vnew = (v*cem) + (invmass*f*pc1) + sd1V - (Vmh*cem); Vmh.x = pc3*sd2X.x + sdpc.x*fg1.x;
posqp = vnew * pc2; Vmh.y = pc3*sd2X.y + sdpc.x*fg1.y;
Vmh.z = pc3*sd2X.z + sdpc.x*fg1.z;
//sd1V = sdpc.y*fg2;
sd1V.x = sdpc.y*fg2.x;
sd1V.y = sdpc.y*fg2.y;
sd1V.z = sdpc.y*fg2.z;
//vnew = (v*cem) + (invmass*f*pc1) + sd1V - (Vmh*cem);
c1 = invmass*pc1;
vnew.x = cem*(v.x - Vmh.x) + c1*f.x + sd1V.x;
vnew.y = cem*(v.y - Vmh.y) + c1*f.y + sd1V.y;
vnew.z = cem*(v.z - Vmh.z) + c1*f.z + sd1V.z;
posqp.x = pc2*vnew.x;
posqp.y = pc2*vnew.y;
posqp.z = pc2*vnew.z;
} }
kernel void kupdate_sd2_fix1( kernel void kupdate_sd2_fix1(
...@@ -167,9 +185,9 @@ kernel void kupdate_sd2_fix1( ...@@ -167,9 +185,9 @@ kernel void kupdate_sd2_fix1(
float3 fg1, fg2; float3 fg1, fg2;
float3 Xmh; float3 Xmh;
igauss = indexof( posq ); igauss = indexof( posq ).xy;
linind_gauss = igauss.y * xstrwidth + igauss.x; linind_gauss = igauss.y * xstrwidth + igauss.x;
linind_gauss = 2 * linind_gauss + goffset; linind_gauss = 2.0f * linind_gauss + goffset;
igauss.y = round( (linind_gauss - fmod( linind_gauss, gstrwidth))/gstrwidth); igauss.y = round( (linind_gauss - fmod( linind_gauss, gstrwidth))/gstrwidth);
igauss.x = linind_gauss - igauss.y * gstrwidth; igauss.x = linind_gauss - igauss.y * gstrwidth;
...@@ -180,8 +198,19 @@ kernel void kupdate_sd2_fix1( ...@@ -180,8 +198,19 @@ kernel void kupdate_sd2_fix1(
igauss.x = linind_gauss - igauss.y * gstrwidth; igauss.x = linind_gauss - igauss.y * gstrwidth;
fg2 = fgauss[ igauss ]; fg2 = fgauss[ igauss ];
v = pc1 * posqp; v.x = pc1*posqp.x;
Xmh = sd1V * pc2 + sdpc.x * fg1; v.y = pc1*posqp.y;
sd2X = sdpc.y * fg2; v.z = pc1*posqp.z;
//Xmh = sd1V * pc2 + sdpc.x * fg1;
Xmh.x = pc2*sd1V.x + sdpc.x*fg1.x;
Xmh.y = pc2*sd1V.y + sdpc.x*fg1.y;
Xmh.z = pc2*sd1V.z + sdpc.x*fg1.z;
//sd2X = sdpc.y * fg2;
sd2X.x = sdpc.y*fg2.x;
sd2X.y = sdpc.y*fg2.y;
sd2X.z = sdpc.y*fg2.z;
posqp2 = posqp + sd2X - Xmh; posqp2 = posqp + sd2X - Xmh;
} }
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment