Skip to content

GitLab

Commit 2c6cff12 authored by Mark Friedrichs's avatar Mark Friedrichs
Browse files

Update

parent f7f79b04
Show whitespace changes
Inline Side-by-side
Showing with 449 additions and 237 deletions
plugins/amoeba/platforms/cuda/src/AmoebaCudaData.cpp
View file @ 2c6cff12
......@@ -95,6 +95,7 @@ KernelImpl* AmoebaCudaData::getAmoebaLocalForcesKernel( void ) const {
void AmoebaCudaData::setLog( FILE* inputLog ) {
log = inputLog;
amoebaGpu->log = inputLog;
}
FILE* AmoebaCudaData::getLog( void ) const {
......
plugins/amoeba/platforms/cuda/src/AmoebaCudaKernelFactory.cpp
View file @ 2c6cff12
......@@ -103,7 +103,7 @@ KernelImpl* AmoebaCudaKernelFactory::createKernelImpl(std::string name, const Pl
if( mapIterator == contextToAmoebaDataMap.end() ){
amoebaCudaData = new AmoebaCudaData( cudaPlatformData );
contextToAmoebaDataMap[&context] = amoebaCudaData;
//amoebaCudaData->setLog( stderr );
amoebaCudaData->setLog( stderr );
amoebaCudaData->setContextImpl( static_cast<void*>(&context) );
//(void) fprintf( stderr, "AmoebaCudaKernelFactory::createKernelImpl amoebaCudaDataV=%p\n", static_cast<void*>(amoebaCudaData) );
} else {
......
plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.cpp
View file @ 2c6cff12
......@@ -46,8 +46,8 @@ using namespace std;
static void computeAmoebaLocalForces( AmoebaCudaData& data ) {
amoebaGpuContext gpu = data.getAmoebaGpu();
if( gpu->log ){
(void) fprintf( gpu->log, "computeAmoebaLocalForces\n" ); (void) fflush( gpu->log );
if( 0 && data.getLog() ){
(void) fprintf( data.getLog(), "computeAmoebaLocalForces\n" ); (void) fflush( data.getLog() );
}
data.initializeGpu();
......@@ -485,7 +485,7 @@ static void computeAmoebaMultipoleForce( AmoebaCudaData& data ) {
amoebaGpuContext gpu = data.getAmoebaGpu();
data.initializeGpu();
if( data.getLog() ){
if( 0 && data.getLog() ){
(void) fprintf( data.getLog(), "computeAmoebaMultipoleForce\n" );
(void) fflush( data.getLog());
}
......@@ -510,7 +510,7 @@ static void computeAmoebaMultipoleForce( AmoebaCudaData& data ) {
kCalculateAmoebaKirkwood(gpu);
}
if( data.getLog() ){
if( 0 && data.getLog() ){
(void) fprintf( data.getLog(), "completed computeAmoebaMultipoleForce\n" );
(void) fflush( data.getLog());
}
......@@ -821,13 +821,13 @@ double CudaCalcAmoebaVdwForceKernel::executeEnergy(ContextImpl& context) {
static void computeAmoebaWcaDispersionForce( AmoebaCudaData& data ) {
data.initializeGpu();
if( data.getLog() ){
if( 0 && data.getLog() ){
(void) fprintf( data.getLog(), "Calling computeAmoebaWcaDispersionForce " ); (void) fflush( data.getLog() );
}
kCalculateAmoebaWcaDispersionForces( data.getAmoebaGpu() );
if( data.getLog() ){
if( 0 && data.getLog() ){
(void) fprintf( data.getLog(), " -- completed\n" ); (void) fflush( data.getLog() );
}
}
......
plugins/amoeba/platforms/cuda/src/kernels/AmoebaGpu.cpp
View file @ 2c6cff12
......@@ -2115,6 +2115,7 @@ void amoebaGpuBuildVdwExclusionList( amoebaGpuContext amoebaGpu, const std::vec
// ---------------------------------------------------------------------------------------
static const std::string methodName = "amoebaGpuBuildVdwExclusionList";
static const int debugOn = 0;
// ---------------------------------------------------------------------------------------
......@@ -2155,7 +2156,7 @@ void amoebaGpuBuildVdwExclusionList( amoebaGpuContext amoebaGpu, const std::vec
// diagnostics
if( amoebaGpu->log ){
if( debugOn && amoebaGpu->log ){
(void) fprintf( amoebaGpu->log, "%s min/max cell indices:\n", methodName.c_str() );
for (int ii = 0; ii < dim; ii++)
{
......@@ -2202,7 +2203,7 @@ void amoebaGpuBuildVdwExclusionList( amoebaGpuContext amoebaGpu, const std::vec
// diagnostics
if( amoebaGpu->log ){
if( debugOn && amoebaGpu->log ){
(void) fprintf( amoebaGpu->log, "%s %d cells w/ exclusions\n", methodName.c_str(), numWithExclusionIndices );
for (int ii = 0; ii < cells; ii++)
{
......@@ -2270,7 +2271,7 @@ void amoebaGpuBuildVdwExclusionList( amoebaGpuContext amoebaGpu, const std::vec
// diagnostics
if( amoebaGpu->log ){
if( debugOn && amoebaGpu->log ){
(void) fprintf( amoebaGpu->log, "%s Echo exclusions\n", methodName.c_str() );
(void) fflush( amoebaGpu->log );
......@@ -2890,6 +2891,7 @@ void amoebaGpuBuildScalingList( amoebaGpuContext amoebaGpu )
// ---------------------------------------------------------------------------------------
static const std::string methodName = "amoebaGpuBuildScalingList";
static const int debugOn = 0;
// ---------------------------------------------------------------------------------------
......@@ -2944,7 +2946,7 @@ void amoebaGpuBuildScalingList( amoebaGpuContext amoebaGpu )
// diagnostics
if( amoebaGpu->log ){
if( debugOn && amoebaGpu->log ){
(void) fprintf( amoebaGpu->log, "%s min/max cell indices:\n", methodName.c_str() );
for (int ii = 0; ii < dim; ii++)
{
......@@ -2980,7 +2982,7 @@ void amoebaGpuBuildScalingList( amoebaGpuContext amoebaGpu )
// diagnostics
#if 1
#if 0
if( 0 && amoebaGpu->log ){
(void) fprintf( amoebaGpu->log, "%s %d cells\n",
methodName.c_str(), numWithScalingIndices );
......@@ -3000,7 +3002,7 @@ void amoebaGpuBuildScalingList( amoebaGpuContext amoebaGpu )
(void) fflush( amoebaGpu->log );
}
#else
if( amoebaGpu->log ){
if( debugOn && amoebaGpu->log ){
(void) fprintf( amoebaGpu->log, "%s %d cells w/ exclusions\n",
methodName.c_str(), numWithScalingIndices );
for (int ii = 0; ii < cells; ii++)
......@@ -3158,7 +3160,7 @@ static unsigned int targetAtoms[2] = { 0, 1};
// diagnostics
if( amoebaGpu->log && 0 ){
if( debugOn && amoebaGpu->log ){
float* pScaleCheckSum = (float*) malloc( sizeof( float )*paddedAtoms );
float* dScaleCheckSum = (float*) malloc( sizeof( float )*paddedAtoms );
......
plugins/amoeba/platforms/cuda/src/kernels/kCalculateAmoebaCudaElectrostatic.h
View file @ 2c6cff12
......@@ -27,6 +27,7 @@
#include "amoebaScaleFactors.h"
__global__
/*
#if (__CUDA_ARCH__ >= 200)
__launch_bounds__(GF1XX_NONBOND_THREADS_PER_BLOCK, 1)
#elif (__CUDA_ARCH__ >= 130)
......@@ -34,6 +35,7 @@ __launch_bounds__(GT2XX_NONBOND_THREADS_PER_BLOCK, 1)
#else
__launch_bounds__(G8X_NONBOND_THREADS_PER_BLOCK, 1)
#endif
*/
void METHOD_NAME(kCalculateAmoebaCudaElectrostatic, Forces_kernel)(
unsigned int* workUnit,
float4* atomCoord,
......
plugins/amoeba/platforms/cuda/src/kernels/kCalculateAmoebaCudaKirkwood.cu
View file @ 2c6cff12
......@@ -2342,7 +2342,6 @@ void kCalculateAmoebaKirkwood( amoebaGpuContext amoebaGpu )
for( int ii = 0; ii < amoebaGpu->gpuContext->sim.paddedNumberOfAtoms; ii++ ){
(void) fprintf( amoebaGpu->log,"Born %6d %16.9e\n", ii,
gpu->psBornRadii->_pSysStream[0][ii] );
}
#endif
......@@ -2350,8 +2349,26 @@ void kCalculateAmoebaKirkwood( amoebaGpuContext amoebaGpu )
if( threadsPerBlock == 0 ){
threadsPerBlock = getThreadsPerBlock( amoebaGpu, sizeof(KirkwoodParticle));
threadsPerBlock = 32;
//unsigned int eDiffhreadsPerBlock = getThreadsPerBlock( amoebaGpu, sizeof(KirkwoodEDiffParticle));
//unsigned int maxThreadsPerBlock = threadsPerBlock> eDiffhreadsPerBlock ? threadsPerBlock : eDiffhreadsPerBlock;
if( amoebaGpu->log ){
#if (__CUDA_ARCH__ >= 200)
unsigned int maxThreads = GF1XX_NONBOND_THREADS_PER_BLOCK;
#elif (__CUDA_ARCH__ >= 130)
unsigned int maxThreads = GT2XX_NONBOND_THREADS_PER_BLOCK;
#else
unsigned int maxThreads = G8X_NONBOND_THREADS_PER_BLOCK;
#endif
(void) fprintf( amoebaGpu->log, "kCalculateAmoebaCudaKirkwood: blcks=%u tds=%u %u bPrWrp=%u atm=%u shrd=%u Ebuf=%u ixnCt=%u workUnits=%u\n",
amoebaGpu->nonbondBlocks, threadsPerBlock, maxThreads, amoebaGpu->bOutputBufferPerWarp,
sizeof(KirkwoodParticle), sizeof(KirkwoodParticle)*threadsPerBlock,
amoebaGpu->energyOutputBuffers, (*gpu->psInteractionCount)[0], gpu->sim.workUnits );
(void) fflush( amoebaGpu->log );
}
}
kClearFields_1( amoebaGpu );
......@@ -2531,7 +2548,6 @@ void kCalculateAmoebaKirkwood( amoebaGpuContext amoebaGpu )
}
delete debugArray;
#endif
// map torques to forces
......
plugins/amoeba/platforms/cuda/src/kernels/kCalculateAmoebaCudaKirkwood.h
View file @ 2c6cff12
......@@ -27,6 +27,7 @@
#include "amoebaScaleFactors.h"
__global__
/*
#if (__CUDA_ARCH__ >= 200)
__launch_bounds__(GF1XX_NONBOND_THREADS_PER_BLOCK, 1)
#elif (__CUDA_ARCH__ >= 130)
......@@ -34,6 +35,7 @@ __launch_bounds__(GT2XX_NONBOND_THREADS_PER_BLOCK, 1)
#else
__launch_bounds__(G8X_NONBOND_THREADS_PER_BLOCK, 1)
#endif
*/
void METHOD_NAME(kCalculateAmoebaCudaKirkwood, Forces_kernel)(
unsigned int* workUnit,
float4* atomCoord,
......
plugins/amoeba/platforms/cuda/src/kernels/kCalculateAmoebaCudaKirkwoodEDiff.cu
View file @ 2c6cff12
......@@ -1156,11 +1156,10 @@ void kCalculateAmoebaKirkwoodEDiff( amoebaGpuContext amoebaGpu )
static const char* methodName = "kCalculateAmoebaKirkwoodEDiff";
static unsigned int threadsPerBlock = 0;
#ifdef AMOEBA_DEBUG
static int timestep = 0;
std::vector<int> fileId;
timestep++;
#ifdef AMOEBA_DEBUG
std::vector<int> fileId;
fileId.resize( 2 );
fileId[0] = timestep;
fileId[1] = 1;
......@@ -1188,21 +1187,21 @@ void kCalculateAmoebaKirkwoodEDiff( amoebaGpuContext amoebaGpu )
unsigned int targetAtom = 0;
#endif
if( amoebaGpu->log ){
(void) fprintf( amoebaGpu->log, "%s %d maxCovalentDegreeSz=%d"
" gamma=%.3e scalingDistanceCutoff=%.3f ZZZ\n",
methodName, gpu->natoms,
amoebaGpu->maxCovalentDegreeSz, amoebaGpu->pGamma,
amoebaGpu->scalingDistanceCutoff );
gpuPrintCudaAmoebaGmxSimulation(amoebaGpu, amoebaGpu->log );
(void) fflush( amoebaGpu->log );
}
kClearFields_3( amoebaGpu, 6 );
if( threadsPerBlock == 0 ){
threadsPerBlock = getThreadsPerBlock( amoebaGpu, sizeof(KirkwoodEDiffParticle));
}
if( amoebaGpu->log && timestep == 1 ){
(void) fprintf( amoebaGpu->log, "kCalculateAmoebaCudaKirkwoodEDiffN2Forces: numBlocks=%u numThreads=%u bufferPerWarp=%u atm=%u shrd=%u Ebuf=%u ixnCt=%u workUnits=%u\n",
amoebaGpu->nonbondBlocks, threadsPerBlock, amoebaGpu->bOutputBufferPerWarp,
sizeof(KirkwoodEDiffParticle), sizeof(KirkwoodEDiffParticle)*threadsPerBlock,
amoebaGpu->energyOutputBuffers, (*gpu->psInteractionCount)[0], gpu->sim.workUnits );
//gpuPrintCudaAmoebaGmxSimulation(amoebaGpu, amoebaGpu->log );
(void) fflush( amoebaGpu->log );
}
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",
......@@ -1414,7 +1413,6 @@ void kCalculateAmoebaKirkwoodEDiff( amoebaGpuContext amoebaGpu )
}
}
#endif
// ---------------------------------------------------------------------------------------
......
plugins/amoeba/platforms/cuda/src/kernels/kCalculateAmoebaCudaMutualInducedAndGkFields.cu
View file @ 2c6cff12
......@@ -756,12 +756,11 @@ static void cudaComputeAmoebaMutualInducedAndGkFieldBySOR( amoebaGpuContext amoe
// ---------------------------------------------------------------------------------------
#ifdef AMOEBA_DEBUG
static const char* methodName = "cudaComputeAmoebaMutualInducedAndGkFieldBySOR";
static int timestep = 0;
std::vector<int> fileId;
timestep++;
static const char* methodName = "cudaComputeAmoebaMutualInducedAndGkFieldBySOR";
#ifdef AMOEBA_DEBUG
std::vector<int> fileId;
fileId.resize( 2 );
fileId[0] = timestep;
fileId[1] = 1;
......@@ -780,7 +779,7 @@ static void cudaComputeAmoebaMutualInducedAndGkFieldBySOR( amoebaGpuContext amoe
if( (numOfElems % numThreads) != 0 )numBlocks++;
#ifdef AMOEBA_DEBUG
if( amoebaGpu->log ){
if( amoebaGpu->log && timestep == 1 ){
(void) fprintf( amoebaGpu->log, "%s %d numOfElems=%d numThreads=%d numBlocks=%d "
"maxIterations=%d targetEpsilon=%.3e\n",
methodName, gpu->natoms, numOfElems, numThreads, numBlocks,
......@@ -981,6 +980,12 @@ time_t start = clock();
amoebaGpu->mutualInducedDone = done;
amoebaGpu->mutualInducedConverged = ( !done || iteration > amoebaGpu->mutualInducedMaxIterations ) ? 0 : 1;
if( amoebaGpu->log ){
(void) fprintf( amoebaGpu->log, "%s done=%d converged=%d iteration=%d eps=%14.7e\n",
methodName, done, amoebaGpu->mutualInducedConverged, iteration, amoebaGpu->mutualInducedCurrentEpsilon );
(void) fflush( amoebaGpu->log );
}
#ifdef AMOEBA_DEBUG
if( 1 ){
std::vector<int> fileId;
......
plugins/amoeba/platforms/cuda/src/kernels/kCalculateAmoebaLocalForces.cu
View file @ 2c6cff12
......@@ -1613,8 +1613,12 @@ void kCalculateAmoebaLocalForces(amoebaGpuContext gpu)
{
if( gpu->log ){
static int call = 0;
if( call == 0 ){
(void) fprintf( gpu->log,"kCalculateAmoebaLocalForces: blks=%u thrds/blk=%u\n",
gpu->gpuContext->sim.blocks, gpu->gpuContext->sim.localForces_threads_per_block); fflush( gpu->log );
call++;
}
}
kCalculateAmoebaLocalForces_kernel<<<gpu->gpuContext->sim.blocks, gpu->gpuContext->sim.localForces_threads_per_block>>>();
......
plugins/amoeba/platforms/cuda/tests/AmoebaTinkerParameterFile.cpp
View file @ 2c6cff12
......@@ -2181,6 +2181,18 @@ static int readAmoebaMultipoleParameters( FILE* filePtr, MapStringInt& forceMap,
}
}
// set parameters if available
MapStringStringI isPresent = inputArgumentMap.find( MUTUAL_INDUCED_MAX_ITERATIONS );
if( isPresent != inputArgumentMap.end() ){
multipoleForce->setMutualInducedMaxIterations( atoi( isPresent->second.c_str() ) );
}
isPresent = inputArgumentMap.find( MUTUAL_INDUCED_TARGET_EPSILON );
if( isPresent != inputArgumentMap.end() ){
multipoleForce->setMutualInducedTargetEpsilon( atof( isPresent->second.c_str() ) );
}
// diagnostics
if( log ){
......@@ -2204,8 +2216,11 @@ static int readAmoebaMultipoleParameters( FILE* filePtr, MapStringInt& forceMap,
//static const unsigned int maxPrint = MAX_PRINT;
static const unsigned int maxPrint = 15;
unsigned int arraySize = static_cast<unsigned int>(multipoleForce->getNumMultipoles());
(void) fprintf( log, "%s: %u sample of AmoebaMultipoleForce parameters\n",
methodName.c_str(), arraySize );
(void) fprintf( log, "%s: %u maxIter=%d targetEps=%14.7e\n",
methodName.c_str(), arraySize,
multipoleForce->getMutualInducedMaxIterations(),
multipoleForce->getMutualInducedTargetEpsilon() );
(void) fprintf( log, "Sample of AmoebaMultipoleForce parameters\n" );
(void) fflush( log );
for( unsigned int ii = 0; ii < arraySize; ii++ ){
......@@ -3461,8 +3476,7 @@ static int addForces( std::vector<Vec3>& forceToAdd, std::vector<Vec3>& forceAcc
return static_cast<int>(forceAccumulator.size());
}
static void getForceStrings( System& system, StringVector& forceStringArray, FILE* log ){
static void getStringForceMap( System& system, MapStringForce& forceMap, FILE* log ){
// print active forces and relevant parameters
......@@ -3477,7 +3491,7 @@ static void getForceStrings( System& system, StringVector& forceStringArray, FIL
try {
AmoebaHarmonicBondForce& harmonicBondForce = dynamic_cast<AmoebaHarmonicBondForce&>(force);
forceStringArray.push_back( AMOEBA_HARMONIC_BOND_FORCE );
forceMap[AMOEBA_HARMONIC_BOND_FORCE] = &force;
hit++;
} catch( std::bad_cast ){
}
......@@ -3489,7 +3503,7 @@ static void getForceStrings( System& system, StringVector& forceStringArray, FIL
try {
AmoebaMultipoleForce& multipoleForce = dynamic_cast<AmoebaMultipoleForce &>(force);
forceStringArray.push_back( AMOEBA_MULTIPOLE_FORCE );
forceMap[AMOEBA_MULTIPOLE_FORCE] = &force;
hit++;
} catch( std::bad_cast ){
}
......@@ -3501,7 +3515,7 @@ static void getForceStrings( System& system, StringVector& forceStringArray, FIL
try {
AmoebaOutOfPlaneBendForce& outOfPlaneBend = dynamic_cast<AmoebaOutOfPlaneBendForce&>(force);
forceStringArray.push_back( AMOEBA_OUT_OF_PLANE_BEND_FORCE );
forceMap[AMOEBA_OUT_OF_PLANE_BEND_FORCE] = &force;
hit++;
} catch( std::bad_cast ){
}
......@@ -3513,7 +3527,7 @@ static void getForceStrings( System& system, StringVector& forceStringArray, FIL
try {
AmoebaPiTorsionForce& piTorsion = dynamic_cast<AmoebaPiTorsionForce&>(force);
forceStringArray.push_back( AMOEBA_PI_TORSION_FORCE );
forceMap[AMOEBA_PI_TORSION_FORCE] = &force;
hit++;
} catch( std::bad_cast ){
}
......@@ -3525,7 +3539,7 @@ static void getForceStrings( System& system, StringVector& forceStringArray, FIL
try {
AmoebaTorsionTorsionForce& torsionTorsion = dynamic_cast<AmoebaTorsionTorsionForce&>(force);
forceStringArray.push_back( AMOEBA_TORSION_TORSION_FORCE );
forceMap[AMOEBA_TORSION_TORSION_FORCE] = &force;
hit++;
} catch( std::bad_cast ){
}
......@@ -3537,7 +3551,7 @@ static void getForceStrings( System& system, StringVector& forceStringArray, FIL
try {
AmoebaTorsionForce& torsion = dynamic_cast<AmoebaTorsionForce&>(force);
forceStringArray.push_back( AMOEBA_TORSION_FORCE );
forceMap[AMOEBA_TORSION_FORCE] = &force;
hit++;
} catch( std::bad_cast ){
}
......@@ -3549,7 +3563,7 @@ static void getForceStrings( System& system, StringVector& forceStringArray, FIL
try {
AmoebaSASAForce& sasaForce = dynamic_cast<AmoebaSASAForce&>(force);
forceStringArray.push_back( AMOEBA_SASA_FORCE );
forceMap[AMOEBA_SASA_FORCE] = &force;
hit++;
} catch( std::bad_cast ){
}
......@@ -3561,7 +3575,7 @@ static void getForceStrings( System& system, StringVector& forceStringArray, FIL
try {
AmoebaStretchBendForce& stretchBend = dynamic_cast<AmoebaStretchBendForce&>(force);
forceStringArray.push_back( AMOEBA_STRETCH_BEND_FORCE );
forceMap[AMOEBA_STRETCH_BEND_FORCE] = &force;
hit++;
} catch( std::bad_cast ){
}
......@@ -3573,7 +3587,7 @@ static void getForceStrings( System& system, StringVector& forceStringArray, FIL
try {
AmoebaVdwForce& vdw = dynamic_cast<AmoebaVdwForce&>(force);
forceStringArray.push_back( AMOEBA_VDW_FORCE );
forceMap[AMOEBA_VDW_FORCE] = &force;
hit++;
} catch( std::bad_cast ){
}
......@@ -3585,7 +3599,7 @@ static void getForceStrings( System& system, StringVector& forceStringArray, FIL
try {
AmoebaWcaDispersionForce& wcaDispersionForce = dynamic_cast<AmoebaWcaDispersionForce&>(force);
forceStringArray.push_back( AMOEBA_WCA_DISPERSION_FORCE );
forceMap[AMOEBA_WCA_DISPERSION_FORCE] = &force;
hit++;
} catch( std::bad_cast ){
}
......@@ -3597,7 +3611,7 @@ static void getForceStrings( System& system, StringVector& forceStringArray, FIL
try {
AmoebaHarmonicAngleForce & harmonicAngleForce = dynamic_cast<AmoebaHarmonicAngleForce&>(force);
forceStringArray.push_back( AMOEBA_HARMONIC_ANGLE_FORCE );
forceMap[AMOEBA_HARMONIC_ANGLE_FORCE] = &force;
hit++;
} catch( std::bad_cast ){
}
......@@ -3609,7 +3623,7 @@ static void getForceStrings( System& system, StringVector& forceStringArray, FIL
try {
AmoebaHarmonicInPlaneAngleForce & harmonicAngleForce = dynamic_cast<AmoebaHarmonicInPlaneAngleForce&>(force);
forceStringArray.push_back( AMOEBA_HARMONIC_IN_PLANE_ANGLE_FORCE );
forceMap[AMOEBA_HARMONIC_IN_PLANE_ANGLE_FORCE] = &force;
hit++;
} catch( std::bad_cast ){
}
......@@ -3620,7 +3634,7 @@ static void getForceStrings( System& system, StringVector& forceStringArray, FIL
if( !hit ){
try {
AmoebaGeneralizedKirkwoodForce& kirkwoodForce = dynamic_cast<AmoebaGeneralizedKirkwoodForce&>(force);
forceStringArray.push_back( AMOEBA_GK_FORCE );
forceMap[AMOEBA_GK_FORCE] = &force;
hit++;
} catch( std::bad_cast ){
}
......@@ -3646,6 +3660,15 @@ static void getForceStrings( System& system, StringVector& forceStringArray, FIL
return;
}
static void getForceStrings( System& system, StringVector& forceStringArray, FILE* log ){
MapStringForce forceMap;
getStringForceMap( system, forceMap, log );
for( MapStringForceI ii = forceMap.begin(); ii != forceMap.end(); ii++ ) {
forceStringArray.push_back( ii->first );
}
}
/**---------------------------------------------------------------------------------------
Read parameter file
......@@ -3922,11 +3945,17 @@ Integrator* readAmoebaParameterFile( const std::string& inputParameterFile, MapS
returnIntegrator = new VerletIntegrator(0.001);
}
if( log ){
(void) fprintf( log, "Potential energies\n" );
double totalPotentialEnergy = 0.0;
if( log )(void) fprintf( log, "Potential energies\n" );
for( MapStringDoubleI ii = potentialEnergy.begin(); ii != potentialEnergy.end(); ii++ ){
(void) fprintf( log, "%30s %14.7e\n", ii->first.c_str(), ii->second );
totalPotentialEnergy += ii->second;
if( log )(void) fprintf( log, "%30s %14.7e\n", ii->first.c_str(), ii->second );
}
potentialEnergy["AllForces"] = totalPotentialEnergy;
if( log ){
(void) fprintf( log, "Total PE %14.7e\n", totalPotentialEnergy );
(void) fprintf( log, "Read %d lines from file=<%s>\n", lineCount, inputParameterFile.c_str() );
(void) fflush( log );
}
......@@ -4322,44 +4351,6 @@ static std::string getIntegratorName( Integrator* integrator ){
return "NotFound";
}
/**---------------------------------------------------------------------------------------
* Set velocities based on temperature
*
* @param system System reference -- retrieve particle masses
* @param velocities array of Vec3 for velocities (size must be set)
* @param temperature temperature
* @param log optional log reference
*
* @return DefaultReturnValue
*
--------------------------------------------------------------------------------------- */
static void setVelocitiesBasedOnTemperature( const System& system, std::vector<Vec3>& velocities, double temperature, FILE* log ) {
// ---------------------------------------------------------------------------------------
static const std::string methodName = "setVelocitiesBasedOnTemperature";
double randomValues[3];
// ---------------------------------------------------------------------------------------
// set velocities based on temperature
temperature *= BOLTZ;
double randMax = static_cast<double>(RAND_MAX);
randMax = 1.0/randMax;
for( unsigned int ii = 0; ii < velocities.size(); ii++ ){
double velocityScale = std::sqrt( temperature/system.getParticleMass(ii) );
randomValues[0] = randMax*( (double) rand() );
randomValues[1] = randMax*( (double) rand() );
randomValues[2] = randMax*( (double) rand() );
velocities[ii] = Vec3( randomValues[0]*velocityScale, randomValues[1]*velocityScale, randomValues[2]*velocityScale );
}
return;
}
/**---------------------------------------------------------------------------------------
* Print Integrator info to log
*
......@@ -4905,6 +4896,164 @@ void testEnergyForcesConsistent( std::string parameterFileName, MapStringInt& fo
}
/**
* Check that energy and force are consistent
*
* @return DefaultReturnValue or ErrorReturnValue
*
*/
System* getCopyOfSystem( System& system, FILE* log ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "getCopyOfSystem";
// ---------------------------------------------------------------------------------------
System* newSystem = new System();
for( int ii = 0; ii < system.getNumParticles(); ii++ ){
newSystem->addParticle( system.getParticleMass( ii ) );
}
for( int ii = 0; ii < system.getNumConstraints(); ii++ ){
int particle1, particle2;
double distance;
system.getConstraintParameters( ii, particle1, particle2, distance );
newSystem->addConstraint( particle1, particle2, distance );
}
return newSystem;
}
/**
* Check that energy and force are consistent
*
* @return DefaultReturnValue or ErrorReturnValue
*
*/
double getEnergyForceBreakdown( Context& context, MapStringDouble& mapEnergies, FILE* log ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "getEnergyForceBreakdown";
// ---------------------------------------------------------------------------------------
int allTypes = State::Positions | State::Velocities | State::Forces | State::Energy;
State state = context.getState( allTypes );
std::vector<Vec3> coordinates = state. getPositions();
System& system = context.getSystem();
MapStringForce forceMap;
getStringForceMap( system, forceMap, log );
MapStringForceI gkIsPresent = forceMap.find( AMOEBA_GK_FORCE );
bool gkIsActive = gkIsPresent == forceMap.end() ? false : true;
double totalEnergy = 0.0;
for( MapStringForceI ii = forceMap.begin(); ii != forceMap.end(); ii++ ){
Force* force = ii->second;
int addForce = 1;
if( gkIsActive ){
if( ii->first == AMOEBA_MULTIPOLE_FORCE ){
addForce = 0;
} else if( ii->first == AMOEBA_GK_FORCE ){
addForce = 2;
}
}
if( addForce ){
System* newSystem = getCopyOfSystem( system, log );
newSystem->addForce( force );
if( addForce == 2 ){
newSystem->addForce( forceMap[AMOEBA_MULTIPOLE_FORCE] );
}
Platform& platform = Platform::getPlatformByName( "Cuda");
platform.setPropertyDefaultValue( "CudaDevice", "3");
//Context newContext = Context( *newSystem, context.getIntegrator(), Platform::getPlatformByName( "Cuda"));
Context newContext = Context( *newSystem, context.getIntegrator(), platform );
newContext.setPositions(coordinates);
State newState = newContext.getState( allTypes );
mapEnergies[ii->first] = newState.getPotentialEnergy();
totalEnergy += newState.getPotentialEnergy();
}
}
return totalEnergy;
}
/**---------------------------------------------------------------------------------------
* Set velocities based on temperature
*
* @param system System reference -- retrieve particle masses
* @param velocities array of Vec3 for velocities (size must be set)
* @param temperature temperature
* @param log optional log reference
*
* @return DefaultReturnValue
*
--------------------------------------------------------------------------------------- */
static void setVelocitiesBasedOnTemperature( const System& system, std::vector<Vec3>& velocities, double temperature, FILE* log ) {
// ---------------------------------------------------------------------------------------
static const std::string methodName = "setVelocitiesBasedOnTemperature";
double randomValues[3];
// ---------------------------------------------------------------------------------------
// set velocities based on temperature
temperature *= BOLTZ;
double randMax = static_cast<double>(RAND_MAX);
randMax = 1.0/randMax;
for( unsigned int ii = 0; ii < velocities.size(); ii++ ){
double velocityScale = std::sqrt( temperature/system.getParticleMass(ii) );
randomValues[0] = randMax*( (double) rand() );
randomValues[1] = randMax*( (double) rand() );
randomValues[2] = randMax*( (double) rand() );
velocities[ii] = Vec3( randomValues[0]*velocityScale, randomValues[1]*velocityScale, randomValues[2]*velocityScale );
}
return;
}
/**
* Check that energy and force are consistent
*
* @return DefaultReturnValue or ErrorReturnValue
*
*/
void writeIntermediateStateFile( Context& context, int currentStep, FILE* intermediateStateFile, FILE* log ){
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "writeIntermediateStateFile";
// ---------------------------------------------------------------------------------------
int allTypes = State::Positions | State::Velocities | State::Forces | State::Energy;
State state = context.getState( allTypes );
const std::vector<Vec3> positions = state.getPositions();
const std::vector<Vec3> velocities = state.getVelocities();
const std::vector<Vec3> forces = state.getForces();
(void) fprintf( intermediateStateFile, "%7u %7d %14.7e %14.7e %14.7e State\n",
positions.size(), currentStep, state.getKineticEnergy(), state.getPotentialEnergy(),
state.getKineticEnergy() + state.getPotentialEnergy() );
for( unsigned int ii = 0; ii < positions.size(); ii++ ){
(void) fprintf( intermediateStateFile, "%7u %14.7e %14.7e %14.7e %14.7e %14.7e %14.7e %14.7e %14.7e %14.7e\n", ii,
positions[ii][0], positions[ii][1], positions[ii][2],
velocities[ii][0], velocities[ii][1], velocities[ii][2],
forces[ii][0], forces[ii][1], forces[ii][2] );
}
return;
}
/**
* Check that energy and force are consistent
*
......@@ -4945,6 +5094,17 @@ void testEnergyConservation( std::string parameterFileName, MapStringInt& forceM
setIntFromMap( inputArgumentMap, "equilibrationTotalSteps", equilibrationTotalSteps );
setIntFromMap( inputArgumentMap, "simulationTotalSteps", simulationTotalSteps );
std::string intermediateStateFileName = "NA";
setStringFromMap( inputArgumentMap, "intermediateStateFileName", intermediateStateFileName );
FILE* intermediateStateFile = NULL;
if( intermediateStateFileName != "NA" ){
#ifdef _MSC_VER
fopen_s( &intermediateStateFile, intermediateStateFileName.c_str(), "w"));
#else
intermediateStateFile = fopen( intermediateStateFileName.c_str(), "w" );
#endif
}
// ---------------------------------------------------------------------------------------
......@@ -4960,8 +5120,10 @@ void testEnergyConservation( std::string parameterFileName, MapStringInt& forceM
System& system = context->getSystem();
int numberOfAtoms = system.getNumParticles();
std::vector<Vec3> velocities;
//velocities.resize( numberOfAtoms );
//_setVelocitiesBasedOnTemperature( system, velocities, initialTemperature, log );
velocities.resize( numberOfAtoms );
double initialT = 300.0;
setVelocitiesBasedOnTemperature( system, velocities, initialT, log );
context->setVelocities(velocities);
// get integrator for equilibration and context
......@@ -4976,6 +5138,7 @@ void testEnergyConservation( std::string parameterFileName, MapStringInt& forceM
int currentStep = 0;
int equilibrationStepsBetweenReports = static_cast<int>(static_cast<double>(equilibrationTotalSteps)*equilibrationStepsBetweenReportsRatio);
if( equilibrationStepsBetweenReports < 1 )equilibrationStepsBetweenReports = 1;
setIntFromMap( inputArgumentMap, "equilibrationStepsBetweenReports", equilibrationStepsBetweenReports );
if( log ){
(void) fprintf( log, "equilibrationTotalSteps=%d equilibrationStepsBetweenReports=%d ratio=%.4f\n",
......@@ -4994,6 +5157,7 @@ void testEnergyConservation( std::string parameterFileName, MapStringInt& forceM
cpuTime = clock();
integrator.step(equilibrationStepsBetweenReports);
//integrator.step(1);
totalEquilibrationTime += clock() - cpuTime;
currentStep += equilibrationStepsBetweenReports;
......@@ -5004,6 +5168,12 @@ void testEnergyConservation( std::string parameterFileName, MapStringInt& forceM
double kineticEnergy = state.getKineticEnergy();
double potentialEnergy = state.getPotentialEnergy();
double totalEnergy = kineticEnergy + potentialEnergy;
if( intermediateStateFile ){
writeIntermediateStateFile( *context, currentStep, intermediateStateFile, log );
}
//MapStringDouble mapEnergies;
//double calculatedPE = getEnergyForceBreakdown( *context, mapEnergies, log );
if( log ){
(void) fprintf( log, "%6d KE=%14.7e PE=%14.7e E=%14.7e\n", currentStep, kineticEnergy, potentialEnergy, totalEnergy );
(void) fflush( log );
......@@ -5075,6 +5245,8 @@ void testEnergyConservation( std::string parameterFileName, MapStringInt& forceM
int simulationStepsBetweenReports = static_cast<int>(static_cast<double>(simulationTotalSteps)*simulationStepsBetweenReportsRatio);
if( simulationStepsBetweenReports < 1 )simulationStepsBetweenReports = 1;
setIntFromMap( inputArgumentMap, "simulationStepsBetweenReports", simulationStepsBetweenReports );
int equilibrationSteps = currentStep;
currentStep = 0;
if( log ){
......@@ -5100,6 +5272,10 @@ void testEnergyConservation( std::string parameterFileName, MapStringInt& forceM
currentStep += simulationStepsBetweenReports;
if( intermediateStateFile ){
writeIntermediateStateFile( *context, (equilibrationSteps+currentStep), intermediateStateFile, log );
}
// record energies
State state = context->getState( State::Energy );
......
plugins/amoeba/platforms/cuda/tests/AmoebaTinkerParameterFile.h
View file @ 2c6cff12
......@@ -84,7 +84,9 @@ static std::string AMOEBA_FIXED_E_GK = "AmoebaF
static std::string AMOEBA_INDUCDED_DIPOLES = "AmoebaInducedDipoles";
static std::string AMOEBA_INDUCDED_DIPOLES_GK = "AmoebaInducedDipoles_GK";
static std::string INCLUDE_OBC_CAVITY_TERM = "INCLUDE_OBC_CAVITY_TERM";
static std::string INCLUDE_OBC_CAVITY_TERM = "includeObcCavityTerm";
static std::string MUTUAL_INDUCED_MAX_ITERATIONS = "mutualInducedMaxIterations";
static std::string MUTUAL_INDUCED_TARGET_EPSILON = "mutualInducedTargetEpsilon";
#define AmoebaHarmonicBondIndex 0
#define AmoebaHarmonicAngleIndex 1
......@@ -147,6 +149,10 @@ typedef std::map< std::string, double > MapStringDouble;
typedef MapStringDouble::iterator MapStringDoubleI;
typedef MapStringDouble::const_iterator MapStringDoubleCI;
typedef std::map< std::string, Force*> MapStringForce;
typedef MapStringForce::iterator MapStringForceI;
typedef MapStringForce::const_iterator MapStringForceCI;
// default return value from methods
static const int DefaultReturnValue = 0;
......
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