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

Added missing N2ByWarp calls

parent ce0cc4e0
...@@ -1907,6 +1907,13 @@ void kCalculateAmoebaKirkwood( amoebaGpuContext amoebaGpu ) ...@@ -1907,6 +1907,13 @@ void kCalculateAmoebaKirkwood( amoebaGpuContext amoebaGpu )
kClearFields_3( amoebaGpu, 6 ); kClearFields_3( amoebaGpu, 6 );
if (gpu->bOutputBufferPerWarp){ if (gpu->bOutputBufferPerWarp){
kCalculateAmoebaCudaKirkwoodN2Forces_kernel<<<amoebaGpu->nonbondBlocks, threadsPerBlock, sizeof(KirkwoodParticle)*threadsPerBlock>>>(
amoebaGpu->psWorkUnit->_pDevStream[0]
#ifdef AMOEBA_DEBUG
, debugArray->_pDevStream[0], targetAtom );
#else
);
#endif
} else { } else {
#ifdef AMOEBA_DEBUG #ifdef AMOEBA_DEBUG
...@@ -1917,7 +1924,7 @@ void kCalculateAmoebaKirkwood( amoebaGpuContext amoebaGpu ) ...@@ -1917,7 +1924,7 @@ void kCalculateAmoebaKirkwood( amoebaGpuContext amoebaGpu )
(void) fflush( amoebaGpu->log ); (void) fflush( amoebaGpu->log );
#endif #endif
kCalculateAmoebaCudaKirkwoodN2Forces_kernel<<<amoebaGpu->nonbondBlocks, threadsPerBlock, sizeof(KirkwoodParticle)*threadsPerBlock>>>( kCalculateAmoebaCudaKirkwoodN2ByWarpForces_kernel<<<amoebaGpu->nonbondBlocks, threadsPerBlock, sizeof(KirkwoodParticle)*threadsPerBlock>>>(
amoebaGpu->psWorkUnit->_pDevStream[0] amoebaGpu->psWorkUnit->_pDevStream[0]
#ifdef AMOEBA_DEBUG #ifdef AMOEBA_DEBUG
, debugArray->_pDevStream[0], targetAtom ); , debugArray->_pDevStream[0], targetAtom );
......
...@@ -1070,28 +1070,22 @@ void kCalculateAmoebaKirkwoodEDiff( amoebaGpuContext amoebaGpu ) ...@@ -1070,28 +1070,22 @@ void kCalculateAmoebaKirkwoodEDiff( amoebaGpuContext amoebaGpu )
} }
if (gpu->bOutputBufferPerWarp){ if (gpu->bOutputBufferPerWarp){
#if 0
(void) fprintf( amoebaGpu->log, "kCalculateAmoebaCudaKirkwoodEDiffN2Forces warp: numBlocks=%u numThreads=%u bufferPerWarp=%u atm=%u shrd=%u Ebuf=%u ixnCt=%u workUnits=%u\n", kCalculateAmoebaCudaKirkwoodEDiffN2ByWarpForces_kernel<<<amoebaGpu->nonbondBlocks, threadsPerBlock, sizeof(KirkwoodEDiffParticle)*threadsPerBlock>>>(
amoebaGpu->nonbondBlocks, amoebaGpu->nonbondElectrostaticThreadsPerBlock, amoebaGpu->bOutputBufferPerWarp, amoebaGpu->psWorkUnit->_pDevStream[0],
sizeof(KirkwoodEDiffParticle), sizeof(KirkwoodEDiffParticle)*amoebaGpu->nonbondElectrostaticThreadsPerBlock, amoebaGpu->energyOutputBuffers, (*gpu->psInteractionCount)[0], gpu->sim.workUnits ); gpu->psPosq4->_pDevStream[0],
(void) fflush( amoebaGpu->log ); amoebaGpu->psLabFrameDipole->_pDevStream[0],
kCalculateAmoebaCudaKirkwoodEDiffN2ByWarpForces_kernel<<<amoebaGpu->nonbondBlocks, amoebaGpu->nonbondElectrostaticThreadsPerBlock, amoebaGpu->psLabFrameQuadrupole->_pDevStream[0],
sizeof(KirkwoodEDiffParticle)*amoebaGpu->nonbondElectrostaticThreadsPerBlock>>>( amoebaGpu->psInducedDipole->_pDevStream[0],
amoebaGpu->psInducedDipolePolar->_pDevStream[0],
amoebaGpu->psWorkUnit->_pDevStream[0], amoebaGpu->psInducedDipoleS->_pDevStream[0],
gpu->psPosq4->_pDevStream[0], amoebaGpu->psInducedDipolePolarS->_pDevStream[0],
amoebaGpu->psLabFrameDipole->_pDevStream[0], amoebaGpu->psWorkArray_3_1->_pDevStream[0],
amoebaGpu->psLabFrameQuadrupole->_pDevStream[0],
amoebaGpu->psInducedDipole->_pDevStream[0],
amoebaGpu->psInducedDipolePolar->_pDevStream[0],
amoebaGpu->psWorkArray_3_1->_pDevStream[0],
amoebaGpu->psWorkArray_3_2->_pDevStream[0],
#ifdef AMOEBA_DEBUG #ifdef AMOEBA_DEBUG
amoebaGpu->psWorkArray_3_2->_pDevStream[0], amoebaGpu->psWorkArray_3_2->_pDevStream[0],
debugArray->_pDevStream[0], targetAtom ); debugArray->_pDevStream[0], targetAtom );
#else #else
amoebaGpu->psWorkArray_3_2->_pDevStream[0] ); amoebaGpu->psWorkArray_3_2->_pDevStream[0] );
#endif
#endif #endif
} else { } else {
......
...@@ -259,9 +259,7 @@ static void cudaComputeAmoebaMutualInducedFieldMatrixMultiply( amoebaGpuContext ...@@ -259,9 +259,7 @@ static void cudaComputeAmoebaMutualInducedFieldMatrixMultiply( amoebaGpuContext
kClearFields_3( amoebaGpu, 2 ); kClearFields_3( amoebaGpu, 2 );
if (gpu->bOutputBufferPerWarp){ if (gpu->bOutputBufferPerWarp){
kCalculateAmoebaMutualInducedFieldN2ByWarp_kernel<<<amoebaGpu->nonbondBlocks, kCalculateAmoebaMutualInducedFieldN2ByWarp_kernel<<<amoebaGpu->nonbondBlocks, amoebaGpu->nonbondThreadsPerBlock, sizeof(MutualInducedParticle)*amoebaGpu->nonbondThreadsPerBlock>>>(
amoebaGpu->nonbondThreadsPerBlock,
sizeof(MutualInducedParticle)*amoebaGpu->nonbondThreadsPerBlock>>>(
amoebaGpu->psWorkUnit->_pDevStream[0], amoebaGpu->psWorkUnit->_pDevStream[0],
amoebaGpu->psWorkArray_3_1->_pDevStream[0], amoebaGpu->psWorkArray_3_1->_pDevStream[0],
#ifdef AMOEBA_DEBUG #ifdef AMOEBA_DEBUG
...@@ -281,9 +279,7 @@ static void cudaComputeAmoebaMutualInducedFieldMatrixMultiply( amoebaGpuContext ...@@ -281,9 +279,7 @@ static void cudaComputeAmoebaMutualInducedFieldMatrixMultiply( amoebaGpuContext
amoebaGpu->energyOutputBuffers, (*gpu->psInteractionCount)[0], gpu->sim.workUnits ); amoebaGpu->energyOutputBuffers, (*gpu->psInteractionCount)[0], gpu->sim.workUnits );
(void) fflush( amoebaGpu->log ); (void) fflush( amoebaGpu->log );
#endif #endif
kCalculateAmoebaMutualInducedFieldN2_kernel<<<amoebaGpu->nonbondBlocks, kCalculateAmoebaMutualInducedFieldN2_kernel<<<amoebaGpu->nonbondBlocks, amoebaGpu->nonbondThreadsPerBlock, sizeof(MutualInducedParticle)*amoebaGpu->nonbondThreadsPerBlock>>>(
amoebaGpu->nonbondThreadsPerBlock,
sizeof(MutualInducedParticle)*amoebaGpu->nonbondThreadsPerBlock>>>(
amoebaGpu->psWorkUnit->_pDevStream[0], amoebaGpu->psWorkUnit->_pDevStream[0],
amoebaGpu->psWorkArray_3_1->_pDevStream[0], amoebaGpu->psWorkArray_3_1->_pDevStream[0],
#ifdef AMOEBA_DEBUG #ifdef AMOEBA_DEBUG
......
...@@ -654,7 +654,6 @@ void kCalculateAmoebaVdw14_7Forces( amoebaGpuContext amoebaGpu, int applyCutoff ...@@ -654,7 +654,6 @@ void kCalculateAmoebaVdw14_7Forces( amoebaGpuContext amoebaGpu, int applyCutoff
LAUNCHERROR("kCalculateAmoebaVdw14_7N2"); LAUNCHERROR("kCalculateAmoebaVdw14_7N2");
} }
#ifdef AMOEBA_DEBUG_PRINT #ifdef AMOEBA_DEBUG_PRINT
if( amoebaGpu->log ){ if( amoebaGpu->log ){
......
...@@ -419,6 +419,17 @@ void kCalculateAmoebaWcaDispersionForces( amoebaGpuContext amoebaGpu ) ...@@ -419,6 +419,17 @@ void kCalculateAmoebaWcaDispersionForces( amoebaGpuContext amoebaGpu )
} }
if (gpu->bOutputBufferPerWarp){ if (gpu->bOutputBufferPerWarp){
kCalculateAmoebaWcaDispersionN2ByWarp_kernel<<<amoebaGpu->nonbondBlocks, threadsPerBlock, sizeof(WcaDispersionParticle)*threadsPerBlock>>>(
amoebaGpu->psWorkUnit->_pDevStream[0],
gpu->psPosq4->_pDevStream[0],
amoebaGpu->psWcaDispersionRadiusEpsilon->_pDevStream[0],
#ifdef AMOEBA_DEBUG
amoebaGpu->psWorkArray_3_1->_pDevStream[0],
debugArray->_pDevStream[0], targetAtom );
#else
amoebaGpu->psWorkArray_3_1->_pDevStream[0] );
#endif
} else { } else {
#ifdef AMOEBA_DEBUG #ifdef AMOEBA_DEBUG
......
...@@ -803,7 +803,36 @@ static int readAmoebaHarmonicBondParameters( FILE* filePtr, MapStringInt& forceM ...@@ -803,7 +803,36 @@ static int readAmoebaHarmonicBondParameters( FILE* filePtr, MapStringInt& forceM
} }
} }
// convert units to kJ-nm from kCal-Angstrom? // check if ZERO_HARMONIC_BOND_IXN is set; used to allow ixn in for
MapStringStringI zeroIxnI = inputArgumentMap.find( ZERO_HARMONIC_BOND_IXN );
int zeroIxn;
if( zeroIxnI != inputArgumentMap.end() ){
zeroIxn = atoi( zeroIxnI->second.c_str() );
} else {
zeroIxn = 0;
}
if( log ){
(void) fprintf( log, "zero harmonic bond ixn=%d.\n", zeroIxn ); (void) fflush( log );
}
// zero ixn
if( zeroIxn ){
double cubic = 0.0;
double quartic = 0.0;
bondForce->setAmoebaGlobalHarmonicBondCubic( cubic );
bondForce->setAmoebaGlobalHarmonicBondQuartic( quartic );
for( int ii = 0; ii < bondForce->getNumBonds(); ii++ ){
int particle1, particle2;
double length, k;
bondForce->getBondParameters( ii, particle1, particle2, length, k );
k = 0.0;
bondForce->setBondParameters( ii, particle1, particle2, length, k );
}
}
if( useOpenMMUnits ){ if( useOpenMMUnits ){
...@@ -2152,6 +2181,7 @@ static int readAmoebaMultipoleParameters( FILE* filePtr, int version, MapStringI ...@@ -2152,6 +2181,7 @@ static int readAmoebaMultipoleParameters( FILE* filePtr, int version, MapStringI
// usePme, aewald, cutoffDistance added w/ Version 1 // usePme, aewald, cutoffDistance added w/ Version 1
if( version > 0 ){ if( version > 0 ){
usePme = atoi( tokens[tokenIndex++].c_str() ); usePme = atoi( tokens[tokenIndex++].c_str() );
aewald = atof( tokens[tokenIndex++].c_str() ); aewald = atof( tokens[tokenIndex++].c_str() );
cutoffDistance = atof( tokens[tokenIndex++].c_str() ); cutoffDistance = atof( tokens[tokenIndex++].c_str() );
...@@ -2169,21 +2199,34 @@ static int readAmoebaMultipoleParameters( FILE* filePtr, int version, MapStringI ...@@ -2169,21 +2199,34 @@ static int readAmoebaMultipoleParameters( FILE* filePtr, int version, MapStringI
grid[2] = atoi( tokens[tokenIndex++].c_str() ); grid[2] = atoi( tokens[tokenIndex++].c_str() );
} }
if( usePme ){ MapStringStringI applyN2I = inputArgumentMap.find( APPLY_N2 );
int applyN2;
if( applyN2I != inputArgumentMap.end() ){
applyN2 = atoi( applyN2I->second.c_str() );
} else {
applyN2 = 0;
}
if( log ){
(void) fprintf( log, "applyN2=%d.\n", applyN2 ); (void) fflush( log );
}
if( usePme && !applyN2 ){
multipoleForce->setNonbondedMethod( AmoebaMultipoleForce::PME ); multipoleForce->setNonbondedMethod( AmoebaMultipoleForce::PME );
multipoleForce->setCutoffDistance( cutoffDistance );
multipoleForce->setAEwald( aewald );
multipoleForce->setPmeBSplineOrder( bsOrder );
multipoleForce->setPmeGridDimensions( grid );
system.setDefaultPeriodicBoxVectors( Vec3(box[0], 0.0, 0.0), Vec3(0.0, box[1], 0.0), Vec3(0.0, 0.0, box[2]) );
} else { } else {
multipoleForce->setNonbondedMethod( AmoebaMultipoleForce::NoCutoff ); multipoleForce->setNonbondedMethod( AmoebaMultipoleForce::NoCutoff );
} }
multipoleForce->setCutoffDistance( cutoffDistance );
multipoleForce->setAEwald( aewald );
multipoleForce->setPmeBSplineOrder( bsOrder );
multipoleForce->setPmeGridDimensions( grid );
system.setDefaultPeriodicBoxVectors( Vec3(box[0], 0.0, 0.0), Vec3(0.0, box[1], 0.0), Vec3(0.0, 0.0, box[2]) );
if( log ){ if( log ){
(void) fprintf( log, "%s number of MultipoleParameter terms=%d usePme=%d aewald=%15.7e cutoffDistance=%12.4f\n", (void) fprintf( log, "%s number of MultipoleParameter terms=%d usePme=%d aewald=%15.7e cutoffDistance=%12.4f\n",
methodName.c_str(), numberOfMultipoles, usePme, aewald, cutoffDistance ); methodName.c_str(), numberOfMultipoles, usePme, aewald, cutoffDistance );
(void) fflush( log ); (void) fflush( log );
} }
for( int ii = 0; ii < numberOfMultipoles; ii++ ){ for( int ii = 0; ii < numberOfMultipoles; ii++ ){
StringVector lineTokens; StringVector lineTokens;
int isNotEof = readLine( filePtr, lineTokens, lineCount, log ); int isNotEof = readLine( filePtr, lineTokens, lineCount, log );
...@@ -2342,8 +2385,8 @@ static int readAmoebaMultipoleParameters( FILE* filePtr, int version, MapStringI ...@@ -2342,8 +2385,8 @@ static int readAmoebaMultipoleParameters( FILE* filePtr, int version, MapStringI
(useOpenMMUnits ? "OpenMM" : "Amoeba") ); (useOpenMMUnits ? "OpenMM" : "Amoeba") );
std::string nonbondedMethod = multipoleForce->getNonbondedMethod( ) == AmoebaMultipoleForce::PME ? "PME" : "NoCutoff"; std::string nonbondedMethod = multipoleForce->getNonbondedMethod( ) == AmoebaMultipoleForce::PME ? "PME" : "NoCutoff";
(void) fprintf( log, "NonbondedMethod=%s aEwald=%15.7e cutoff=%15.7e.\n", nonbondedMethod.c_str(), (void) fprintf( log, "NonbondedMethod=%s aEwald=%15.7e cutoff=%15.7e tol=%15.7e.\n", nonbondedMethod.c_str(),
multipoleForce->getAEwald(), multipoleForce->getCutoffDistance() ); multipoleForce->getAEwald(), multipoleForce->getCutoffDistance(), multipoleForce->getEwaldErrorTolerance() );
Vec3 a,b,c; Vec3 a,b,c;
system.getDefaultPeriodicBoxVectors( a, b, c ); system.getDefaultPeriodicBoxVectors( a, b, c );
...@@ -4743,6 +4786,18 @@ void testUsingAmoebaTinkerParameterFile( const std::string& amoebaTinkerParamete ...@@ -4743,6 +4786,18 @@ void testUsingAmoebaTinkerParameterFile( const std::string& amoebaTinkerParamete
std::vector<Vec3> coordinates = state.getPositions(); std::vector<Vec3> coordinates = state.getPositions();
std::vector<Vec3> forces = state.getForces(); std::vector<Vec3> forces = state.getForces();
// check if ZERO_HARMONIC_BOND_IXN is set; used to allow ixn in for
MapStringStringI zeroIxnI = inputArgumentMap.find( ZERO_HARMONIC_BOND_IXN );
int zeroIxn;
if( zeroIxnI != inputArgumentMap.end() ){
zeroIxn = atoi( zeroIxnI->second.c_str() );
} else {
zeroIxn = 0;
}
if( log ){
(void) fprintf( log, "zero harmonic bond ixn=%d.\n", zeroIxn ); (void) fflush( log );
}
// get list of forces and then accumulate expected energies/forces // get list of forces and then accumulate expected energies/forces
StringVector forceList; StringVector forceList;
...@@ -4753,8 +4808,10 @@ void testUsingAmoebaTinkerParameterFile( const std::string& amoebaTinkerParamete ...@@ -4753,8 +4808,10 @@ void testUsingAmoebaTinkerParameterFile( const std::string& amoebaTinkerParamete
forceList.push_back( AMOEBA_GK_CAVITY_FORCE ); forceList.push_back( AMOEBA_GK_CAVITY_FORCE );
activeForceNames += AMOEBA_GK_CAVITY_FORCE + ":"; activeForceNames += AMOEBA_GK_CAVITY_FORCE + ":";
} else { } else {
forceList.push_back( ii->first ); if( !zeroIxn || ii->first != AMOEBA_HARMONIC_BOND_FORCE ){
activeForceNames += ii->first + ":"; forceList.push_back( ii->first );
activeForceNames += ii->first + ":";
}
} }
} }
} }
...@@ -4826,7 +4883,7 @@ void testUsingAmoebaTinkerParameterFile( const std::string& amoebaTinkerParamete ...@@ -4826,7 +4883,7 @@ void testUsingAmoebaTinkerParameterFile( const std::string& amoebaTinkerParamete
bool badMatch = (cutoffDelta < relativeDelta) && (sumNorms > 0.1) ? true : false; bool badMatch = (cutoffDelta < relativeDelta) && (sumNorms > 0.1) ? true : false;
badMatch = badMatch || (normF1 == 0.0 && normF2 > 0.0) || (normF2 == 0.0 && normF1 > 0.0); badMatch = badMatch || (normF1 == 0.0 && normF2 > 0.0) || (normF2 == 0.0 && normF1 > 0.0);
if( badMatch || showAll ){ if( badMatch || showAll ){
(void) fprintf( filePtr, "%6u %10.3e %10.3e [%15.7e %15.7e %15.7e] [%15.7e %15.7e %15.7e] %s\n", ii, relativeDelta, delta, (void) fprintf( filePtr, "%6u %10.3e %10.3e [%15.7e %15.7e %15.7e] [%15.7e %15.7e %15.7e] %s\n", ii, relativeDelta, delta,
expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2], expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2],
forceConversion*forces[ii][0], forceConversion*forces[ii][1], forceConversion*forces[ii][2], forceConversion*forces[ii][0], forceConversion*forces[ii][1], forceConversion*forces[ii][2],
( (showAll && badMatch) ? " XXX" : "") ); ( (showAll && badMatch) ? " XXX" : "") );
......
...@@ -91,6 +91,8 @@ static std::string AMOEBA_INDUCDED_DIPOLES_GK = "AmoebaI ...@@ -91,6 +91,8 @@ static std::string AMOEBA_INDUCDED_DIPOLES_GK = "AmoebaI
static std::string INCLUDE_OBC_CAVITY_TERM = "includeObcCavityTerm"; static std::string INCLUDE_OBC_CAVITY_TERM = "includeObcCavityTerm";
static std::string MUTUAL_INDUCED_MAX_ITERATIONS = "mutualInducedMaxIterations"; static std::string MUTUAL_INDUCED_MAX_ITERATIONS = "mutualInducedMaxIterations";
static std::string MUTUAL_INDUCED_TARGET_EPSILON = "mutualInducedTargetEpsilon"; static std::string MUTUAL_INDUCED_TARGET_EPSILON = "mutualInducedTargetEpsilon";
static std::string APPLY_N2 = "applyN2";
static std::string ZERO_HARMONIC_BOND_IXN = "zeroHarmonicBondIxn";
#define AmoebaHarmonicBondIndex 0 #define AmoebaHarmonicBondIndex 0
#define AmoebaHarmonicAngleIndex 1 #define AmoebaHarmonicAngleIndex 1
......
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