Skip to content

GitLab

Commit 032e18de authored by Mark Friedrichs's avatar Mark Friedrichs
Browse files

Removed debug code 

In Born radii calculation for particle i, removed factor of lambda associated w/ particle i
parent 076eb95e
Show whitespace changes
Inline Side-by-side
Showing with 929 additions and 1855 deletions
plugins/freeEnergy/platforms/cuda/src/kernels/kCalculateCDLJObcGbsaSoftcoreForces1.cu
View file @ 032e18de
......@@ -38,7 +38,6 @@
#include <sstream>
#define USE_SOFTCORE_LJ
//#define DEBUG
struct Atom {
float x;
......@@ -134,52 +133,11 @@ void kCalculateCDLJObcGbsaSoftcoreForces1( freeEnergyGpuContext freeEnergyGpu )
}
threadsPerBlock = threadsPerBlockPerMethod[methodIndex];
#ifdef DEBUG
fprintf( stderr, "kCalculateCDLJObcGbsaSoftcoreForces1 blks=%u thread/block=%u %u shMem=%u nbMethod==%d warp=%u\n",
gpu->sim.nonbond_blocks, threadsPerBlock, gpu->sim.nonbond_threads_per_block, sizeof(Atom)*threadsPerBlock,
freeEnergyGpu->freeEnergySim.nonbondedMethod, gpu->bOutputBufferPerWarp);
int psize = gpu->sim.paddedNumberOfAtoms;
CUDAStream<float4>* pdE1 = new CUDAStream<float4>( psize, 1, "pdE");
CUDAStream<float4>* pdE2 = new CUDAStream<float4>( psize, 1, "pdE");
float bF,bR;
float bF1,b2;
float ratio;
float atomicRadii;
showWorkUnitsFreeEnergy( freeEnergyGpu, 1 );
for( int ii = 0; ii < psize; ii++ ){
pdE1->_pSysData[ii].x = 0.0f;
pdE1->_pSysData[ii].y = 0.001f;
pdE1->_pSysData[ii].z = 0.001f;
pdE1->_pSysData[ii].w = 0.001f;
pdE2->_pSysData[ii].x = 0.0f;
pdE2->_pSysData[ii].y = 0.001f;
pdE2->_pSysData[ii].z = 0.001f;
pdE2->_pSysData[ii].w = 0.001f;
}
pdE1->Upload();
pdE2->Upload();
#endif
switch( freeEnergyGpu->freeEnergySim.nonbondedMethod )
{
case FREE_ENERGY_NO_CUTOFF:
// use softcore LJ potential
#ifdef DEBUG
(void) fprintf( stderr, "kCalculateCDLJObcGbsaSoftcoreForces1 ver=%u blks=%u threadsPerBlock=%u shMem=%u %u wrp=%u\n", gpu->sm_version,
gpu->sim.nonbond_blocks, threadsPerBlock, sizeof(Atom)*threadsPerBlock, gpu->sharedMemoryPerBlock, gpu->bOutputBufferPerWarp ); fflush( stderr );
if (gpu->bOutputBufferPerWarp)
kCalculateCDLJObcGbsaSoftcoreN2ByWarpForces1_kernel<<<gpu->sim.nonbond_blocks, threadsPerBlock,
sizeof(Atom)*threadsPerBlock>>>(gpu->sim.pWorkUnit,pdE1->_pDevData, pdE2->_pDevData);
else
kCalculateCDLJObcGbsaSoftcoreN2Forces1_kernel<<<gpu->sim.nonbond_blocks, threadsPerBlock,
sizeof(Atom)*threadsPerBlock>>>(gpu->sim.pWorkUnit, pdE1->_pDevData, pdE2->_pDevData);
#else
if (gpu->bOutputBufferPerWarp)
kCalculateCDLJObcGbsaSoftcoreN2ByWarpForces1_kernel<<<gpu->sim.nonbond_blocks, threadsPerBlock,
sizeof(Atom)*threadsPerBlock>>>(gpu->sim.pWorkUnit );
......@@ -187,33 +145,17 @@ pdE2->Upload();
kCalculateCDLJObcGbsaSoftcoreN2Forces1_kernel<<<gpu->sim.nonbond_blocks, threadsPerBlock,
sizeof(Atom)*threadsPerBlock>>>(gpu->sim.pWorkUnit );
#endif
LAUNCHERROR("kCalculateCDLJObcGbsaSoftcoreForces1");
break;
case FREE_ENERGY_CUTOFF:
#ifdef DEBUG
(void) fprintf( stderr, "kCalculateCDLJObcGbsaSoftcoreCutoffForces1 %6d blks=%u nonbond_threads_per_block=%5u shMem=%5u\n",
gpu->natoms, gpu->sim.nonbond_blocks, threadsPerBlock, (sizeof(Atom)+sizeof(float))*threadsPerBlock);
(void) fflush( stderr );
if (gpu->bOutputBufferPerWarp)
kCalculateCDLJObcGbsaSoftcoreCutoffByWarpForces1_kernel<<<gpu->sim.nonbond_blocks, threadsPerBlock,
(sizeof(Atom)+sizeof(float))*threadsPerBlock>>>(gpu->sim.pInteractingWorkUnit, pdE1->_pDevData, pdE2->_pDevData);
else
kCalculateCDLJObcGbsaSoftcoreCutoffForces1_kernel<<<gpu->sim.nonbond_blocks, threadsPerBlock,
(sizeof(Atom)+sizeof(float))*threadsPerBlock>>>(gpu->sim.pInteractingWorkUnit, pdE1->_pDevData, pdE2->_pDevData);
#else
if (gpu->bOutputBufferPerWarp)
kCalculateCDLJObcGbsaSoftcoreCutoffByWarpForces1_kernel<<<gpu->sim.nonbond_blocks, threadsPerBlock,
(sizeof(Atom)+sizeof(float))*threadsPerBlock>>>(gpu->sim.pInteractingWorkUnit );
else
kCalculateCDLJObcGbsaSoftcoreCutoffForces1_kernel<<<gpu->sim.nonbond_blocks, threadsPerBlock,
(sizeof(Atom)+sizeof(float))*threadsPerBlock>>>(gpu->sim.pInteractingWorkUnit );
#endif
LAUNCHERROR("kCalculateCDLJObcGbsaSoftcoreCutoffForces1");
......@@ -221,85 +163,15 @@ pdE2->Upload();
case FREE_ENERGY_PERIODIC:
#ifdef DEBUG
if (gpu->bOutputBufferPerWarp)
kCalculateCDLJObcGbsaSoftcorePeriodicByWarpForces1_kernel<<<gpu->sim.nonbond_blocks, threadsPerBlock,
(sizeof(Atom)+sizeof(float))*threadsPerBlock>>>(gpu->sim.pInteractingWorkUnit, pdE1->_pDevData, pdE2->_pDevData);
else
kCalculateCDLJObcGbsaSoftcorePeriodicForces1_kernel<<<gpu->sim.nonbond_blocks, threadsPerBlock,
(sizeof(Atom)+sizeof(float))*threadsPerBlock>>>(gpu->sim.pInteractingWorkUnit, pdE1->_pDevData, pdE2->_pDevData);
#else
if (gpu->bOutputBufferPerWarp)
kCalculateCDLJObcGbsaSoftcorePeriodicByWarpForces1_kernel<<<gpu->sim.nonbond_blocks, threadsPerBlock,
(sizeof(Atom)+sizeof(float))*threadsPerBlock>>>(gpu->sim.pInteractingWorkUnit);
else
kCalculateCDLJObcGbsaSoftcorePeriodicForces1_kernel<<<gpu->sim.nonbond_blocks, threadsPerBlock,
(sizeof(Atom)+sizeof(float))*threadsPerBlock>>>(gpu->sim.pInteractingWorkUnit);
#endif
LAUNCHERROR("kCalculateCDLJObcGbsaSoftcorePeriodicForces1");
break;
}
#ifdef DEBUG
/*
gpu->psBornForce->Download();
freeEnergyGpu->psSwitchDerivative->Download();
gpu->psBornRadii->Download();
fprintf( stderr, "XX bR=%15.7e swd=%15.7e\n", gpu->psBornRadii->_pSysData[0], freeEnergyGpu->psSwitchDerivative->_pSysData[0] );
for( int ii = 0; ii < gpu->sim.nonbondOutputBuffers; ii++ ){
fprintf( stderr, "strx %4d %15.7e %15.7e %15.7e %15.7e\n", ii,
gpu->psBornForce->_pSysStream[ii][0],
gpu->psBornForce->_pSysStream[ii][1],
gpu->psBornForce->_pSysStream[ii][2],
gpu->psBornForce->_pSysStream[ii][3] );
if( gpu->natoms > 1984 ){
int idx = 1983;
fprintf( stderr, "stry %4d %15.7e %15.7e %15.7e %15.7e %5d\n", ii,
gpu->psBornForce->_pSysStream[ii][idx+0],
gpu->psBornForce->_pSysStream[ii][idx+1],
gpu->psBornForce->_pSysStream[ii][idx+2],
gpu->psBornForce->_pSysStream[ii][idx+3], idx );
}
}
int bufferI = 62;
if( bufferI < gpu->sim.nonbondOutputBuffers ){
fprintf( stderr, "BufferI %4d \n", bufferI );
for( int ii = 0; ii < gpu->sim.paddedNumberOfAtoms; ii++ ){
fprintf( stderr, "strz %4d %15.7e\n", ii, gpu->psBornForce->_pSysStream[bufferI][ii] );
}
}
*/
pdE1->Download();
pdE2->Download();
gpu->psPosq4->Download();
gpu->psGBVIData->Download();
gpu->psBornRadii->Download();
gpu->psBornForce->Download();
freeEnergyGpu->psSwitchDerivative->Download();
fprintf( stderr, "PdeCud %d\n", TARGET );
bF = 0.0;
int count =0;
for( int ii = 0; ii < psize; ii++ ){
bF += pdE1->_pSysData[ii].x;
if( fabs( pdE1->_pSysData[ii].w ) > 1.0e-03 && fabs( pdE1->_pSysData[ii].x ) > 0.0 ){
count++;
fprintf( stderr, "%4d %4d %15.7e %15.7e %15.7e %15.7e %15.7e %15.7e %15.7e %15.7e bF=%15.7e swd=%15.7e\n", count, ii,
pdE1->_pSysData[ii].x, pdE1->_pSysData[ii].y, pdE1->_pSysData[ii].z, pdE1->_pSysData[ii].w,
pdE2->_pSysData[ii].x, pdE2->_pSysData[ii].y, pdE2->_pSysData[ii].z, pdE2->_pSysData[ii].w, gpu->psBornForce->_pSysData[ii],
freeEnergyGpu->psSwitchDerivative->_pSysData[ii] );
}
}
bR = gpu->psBornRadii->_pSysData[TARGET];
atomicRadii = gpu->psGBVIData->_pSysData[TARGET].x;
ratio = (atomicRadii/bR);
bF1 = bF + (3.0f*gpu->psGBVIData->_pSysData[TARGET].z*ratio*ratio*ratio)/bR;
b2 = bR*bR;
bF1 *= (1.0f/3.0f)*b2*b2;
fprintf( stderr, "sumbF Cud %6d count=%d %15.7e %15.7e %15.7e\n", TARGET, count, bF, bF1, bR);
#endif
}
plugins/freeEnergy/platforms/cuda/src/kernels/kCalculateCDLJObcGbsaSoftcoreForces1.h
View file @ 032e18de
......@@ -36,9 +36,6 @@
#include "kSoftcoreLJ.h"
#endif
#undef TARGET
#define TARGET 1926
__global__
#if (__CUDA_ARCH__ >= 200)
__launch_bounds__(GF1XX_NONBOND_THREADS_PER_BLOCK, 1)
......@@ -57,15 +54,12 @@ void METHOD_NAME(kCalculateCDLJObcGbsaSoftcore, Forces1_kernel)(unsigned int* wo
unsigned int totalWarps = gridDim.x*blockDim.x/GRID;
unsigned int warp = (blockIdx.x*blockDim.x+threadIdx.x)/GRID;
//unsigned int totalWarps = cSim.nonbond_blocks*cSim.nonbond_threads_per_block/GRID;
//unsigned int warp = (blockIdx.x*blockDim.x+threadIdx.x)/GRID;
unsigned int numWorkUnits = cSim.pInteractionCount[0];
unsigned int pos = warp*numWorkUnits/totalWarps;
unsigned int end = (warp+1)*numWorkUnits/totalWarps;
float CDLJObcGbsa_energy;
float energy = 0.0f;
#ifdef USE_CUTOFF
//float* tempBuffer = (float*) &sA[cSim.nonbond_threads_per_block];
float* tempBuffer = (float*) &sA[blockDim.x];
#endif
......@@ -171,24 +165,6 @@ void METHOD_NAME(kCalculateCDLJObcGbsaSoftcore, Forces1_kernel)(unsigned int* wo
dGpol_dalpha2_ij = 0.0f;
}
af.w += dGpol_dalpha2_ij * psA[j].br;
#ifdef DEBUG
int jIdx = j;
if( i == TARGET ){
int tjj = y+jIdx;
pdE1[tjj].x = dGpol_dalpha2_ij * psA[jIdx].br;
pdE1[tjj].y = sqrt(r2);
pdE1[tjj].w = 1.0f;
}
if( (y+jIdx) == TARGET ){
int tjj = i;
pdE1[tjj].x = dGpol_dalpha2_ij * psA[jIdx].br;
pdE1[tjj].y = sqrt(r2);
pdE1[tjj].w = -1.0f;
}
#endif
energy += 0.5f*CDLJObcGbsa_energy;
// Add Forces
......@@ -278,22 +254,6 @@ pdE1[tjj].w = -1.0f;
dGpol_dalpha2_ij = 0.0f;
}
#ifdef DEBUG
int jIdx = j;
if( i == TARGET ){
int tjj = (y+jIdx);
pdE1[tjj].x = dGpol_dalpha2_ij * psA[jIdx].br;
pdE1[tjj].y = sqrt(r2);
pdE1[tjj].w = 2.0f;
}
if( (y+jIdx) == TARGET ){
int tjj = i;
pdE1[tjj].x = dGpol_dalpha2_ij * psA[jIdx].br;
pdE1[tjj].y = sqrt(r2);
pdE1[tjj].w = -2.0f;
}
#endif
af.w += dGpol_dalpha2_ij * psA[j].br;
energy += 0.5f*CDLJObcGbsa_energy;
......@@ -428,26 +388,6 @@ pdE1[tjj].w = -2.0f;
af.w += dGpol_dalpha2_ij * psA[tj].br;
energy += CDLJObcGbsa_energy;
#ifdef DEBUG
int jIdx = tj;
if( i == TARGET ){
int tjj = y+jIdx;
pdE1[tjj].x = dGpol_dalpha2_ij * psA[jIdx].br;
pdE1[tjj].y = sqrt(r2);
pdE1[tjj].w = 3.0f;
}
if( (y+jIdx) == TARGET ){
int tjj = i;
pdE1[tjj].x = dGpol_dalpha2_ij * br;
pdE1[tjj].y = sqrt(r2);
pdE1[tjj].w = -3.0f;
}
#endif
// Add forces
dx *= dEdR;
......@@ -545,24 +485,6 @@ pdE1[tjj].w = -3.0f;
if (tgx == 0)
psA[j].fb += tempBuffer[threadIdx.x] + tempBuffer[threadIdx.x+16];
#ifdef DEBUG
int jIdx = j;
if( i == TARGET ){
int tjj = y+jIdx;
pdE1[tjj].x = dGpol_dalpha2_ij * psA[j].br;
pdE1[tjj].y = sqrt(r2);
pdE1[tjj].w = 4.0f;
}
if( (y+jIdx) == TARGET ){
int tjj = i;
pdE1[tjj].x = tempBuffer[threadIdx.x] + tempBuffer[threadIdx.x+16];
pdE1[tjj].y = sqrt(r2);
pdE1[tjj].w = -4.0f;
}
#endif
energy += CDLJObcGbsa_energy;
// Add forces
......@@ -678,26 +600,6 @@ pdE1[tjj].w = -4.0f;
dGpol_dalpha2_ij = 0.0f;
}
#ifdef DEBUG
int jIdx = tj;
if( i == TARGET ){
int tjj = y+jIdx;
pdE1[tjj].x = dGpol_dalpha2_ij * psA[tj].br;
pdE1[tjj].y = sqrt(r2);
pdE1[tjj].z = dGpol_dalpha2_ij;
pdE1[tjj].w = 6.0f;
}
if( (y+jIdx) == TARGET ){
int tjj = i;
pdE1[tjj].x = dGpol_dalpha2_ij * br;
pdE1[tjj].y = sqrt(r2);
pdE1[tjj].z = dGpol_dalpha2_ij;
pdE1[tjj].w = -6.0f;
}
#endif
af.w += dGpol_dalpha2_ij * psA[tj].br;
psA[tj].fb += dGpol_dalpha2_ij * br;
energy += CDLJObcGbsa_energy;
......
plugins/freeEnergy/platforms/cuda/src/kernels/kCalculateGBVISoftcoreBornSum.cu
View file @ 032e18de
......@@ -39,7 +39,6 @@
#define PARAMETER_PRINT 0
#define MAX_PARAMETER_PRINT 10
//#define DEBUG
static __constant__ cudaGmxSimulation cSim;
static __constant__ cudaFreeEnergyGmxSimulation gbviSimDev;
......@@ -525,35 +524,6 @@ void kCalculateGBVISoftcoreBornSum( freeEnergyGpuContext freeEnergyGpu )
}
threadsPerBlock = threadsPerBlockPerMethod[methodIndex];
#ifdef DEBUG
fprintf( stderr, "kCalculateGBVISoftcoreBornSum blks=%u threadsPerBlock=%u %u shMem=%u\n",
gpu->sim.nonbond_blocks, threadsPerBlock, gpu->sim.nonbond_threads_per_block, (sizeof(Atom)+sizeof(float))*threadsPerBlock ); fflush( stderr );
int psize = gpu->sim.paddedNumberOfAtoms;
CUDAStream<float4>* pdE1 = new CUDAStream<float4>( psize, 1, "pdE");
CUDAStream<float4>* pdE2 = new CUDAStream<float4>( psize, 1, "pdE");
float bF;
float bF1;
showWorkUnitsFreeEnergy( freeEnergyGpu, 1 );
for( int ii = 0; ii < psize; ii++ ){
pdE1->_pSysData[ii].x = 0.0f;
pdE1->_pSysData[ii].y = 0.001f;
pdE1->_pSysData[ii].z = 0.001f;
pdE1->_pSysData[ii].w = 0.001f;
pdE2->_pSysData[ii].x = 0.001f;
pdE2->_pSysData[ii].y = 0.001f;
pdE2->_pSysData[ii].z = 0.001f;
pdE2->_pSysData[ii].w = 0.001f;
}
pdE1->Upload();
pdE2->Upload();
#endif
kClearGBVISoftcoreBornSum( gpu );
LAUNCHERROR("kClearGBVIBornSum from kCalculateGBVISoftcoreBornSum");
......@@ -561,15 +531,6 @@ pdE2->Upload();
{
case FREE_ENERGY_NO_CUTOFF:
#ifdef DEBUG
if (gpu->bOutputBufferPerWarp){
kCalculateGBVISoftcoreN2ByWarpBornSum_kernel<<<gpu->sim.nonbond_blocks, threadsPerBlock,
sizeof(Atom)*threadsPerBlock>>>(gpu->sim.pWorkUnit, pdE1->_pDevData, pdE2->_pDevData);
} else {
kCalculateGBVISoftcoreN2BornSum_kernel<<<gpu->sim.nonbond_blocks, threadsPerBlock,
sizeof(Atom)*threadsPerBlock>>>(gpu->sim.pWorkUnit, pdE1->_pDevData, pdE2->_pDevData);
}
#else
if (gpu->bOutputBufferPerWarp){
kCalculateGBVISoftcoreN2ByWarpBornSum_kernel<<<gpu->sim.nonbond_blocks, threadsPerBlock,
sizeof(Atom)*threadsPerBlock>>>(gpu->sim.pWorkUnit);
......@@ -577,8 +538,6 @@ pdE2->Upload();
kCalculateGBVISoftcoreN2BornSum_kernel<<<gpu->sim.nonbond_blocks, threadsPerBlock,
sizeof(Atom)*threadsPerBlock>>>(gpu->sim.pWorkUnit);
}
#endif
break;
case FREE_ENERGY_CUTOFF:
......@@ -591,27 +550,13 @@ pdE2->Upload();
kFindInteractionsWithinBlocksCutoff_kernel<<<gpu->sim.nonbond_blocks, threadsPerBlock,
sizeof(unsigned int)*threadsPerBlock>>>(gpu->sim.pInteractingWorkUnit);
#ifdef DEBUG
(void) fprintf( stderr, "kCalculateGBVISoftcoreBornSum cutoff=%15.7e warp=%u GridBoundingBox.length=%u interaction_blocks=%u interaction_threads_per_block=%u nonbond_blocks=%u nonbond_threads_per_block=%u\n",
gpu->sim.nonbondedCutoffSqr, gpu->bOutputBufferPerWarp, gpu->psGridBoundingBox->_length, gpu->sim.interaction_blocks,
gpu->sim.interaction_threads_per_block, gpu->sim.nonbond_blocks, threadsPerBlock ); fflush( stderr );
if (gpu->bOutputBufferPerWarp)
kCalculateGBVISoftcoreCutoffByWarpBornSum_kernel<<<gpu->sim.nonbond_blocks, threadsPerBlock,
(sizeof(Atom)+sizeof(float))*threadsPerBlock>>>(gpu->sim.pInteractingWorkUnit, pdE1->_pDevData, pdE2->_pDevData);
else
kCalculateGBVISoftcoreCutoffBornSum_kernel<<<gpu->sim.nonbond_blocks, threadsPerBlock,
(sizeof(Atom)+sizeof(float))*threadsPerBlock>>>(gpu->sim.pInteractingWorkUnit, pdE1->_pDevData, pdE2->_pDevData );
#else
if (gpu->bOutputBufferPerWarp)
kCalculateGBVISoftcoreCutoffByWarpBornSum_kernel<<<gpu->sim.nonbond_blocks, threadsPerBlock,
(sizeof(Atom)+sizeof(float))*threadsPerBlock>>>(gpu->sim.pInteractingWorkUnit);
else
kCalculateGBVISoftcoreCutoffBornSum_kernel<<<gpu->sim.nonbond_blocks, threadsPerBlock,
(sizeof(Atom)+sizeof(float))*threadsPerBlock>>>(gpu->sim.pInteractingWorkUnit );
#endif
break;
case FREE_ENERGY_PERIODIC:
......@@ -624,21 +569,13 @@ pdE2->Upload();
kFindInteractionsWithinBlocksPeriodic_kernel<<<gpu->sim.nonbond_blocks, threadsPerBlock,
sizeof(unsigned int)*threadsPerBlock>>>(gpu->sim.pInteractingWorkUnit);
#ifdef DEBUG
if (gpu->bOutputBufferPerWarp)
kCalculateGBVISoftcorePeriodicByWarpBornSum_kernel<<<gpu->sim.nonbond_blocks, threadsPerBlock,
(sizeof(Atom)+sizeof(float))*threadsPerBlock>>>(gpu->sim.pInteractingWorkUnit, pdE1->_pDevData, pdE2->_pDevData );
else
kCalculateGBVISoftcorePeriodicBornSum_kernel<<<gpu->sim.nonbond_blocks, threadsPerBlock,
(sizeof(Atom)+sizeof(float))*threadsPerBlock>>>(gpu->sim.pInteractingWorkUnit, pdE1->_pDevData, pdE2->_pDevData );
#else
if (gpu->bOutputBufferPerWarp)
kCalculateGBVISoftcorePeriodicByWarpBornSum_kernel<<<gpu->sim.nonbond_blocks, threadsPerBlock,
(sizeof(Atom)+sizeof(float))*threadsPerBlock>>>(gpu->sim.pInteractingWorkUnit );
else
kCalculateGBVISoftcorePeriodicBornSum_kernel<<<gpu->sim.nonbond_blocks, threadsPerBlock,
(sizeof(Atom)+sizeof(float))*threadsPerBlock>>>(gpu->sim.pInteractingWorkUnit );
#endif
break;
default:
......@@ -647,23 +584,4 @@ pdE2->Upload();
}
LAUNCHERROR("kCalculateGBVISoftcoreBornSum");
#ifdef DEBUG
pdE1->Download();
pdE2->Download();
fprintf( stderr, "bSum Cud method=%u warp=%u\n", freeEnergyGpu->freeEnergySim.nonbondedMethod, gpu->bOutputBufferPerWarp );
bF = 0.0;
bF1 = 0.0;
for( int ii = 0; ii < gpu->natoms; ii++ ){
if( fabsf( pdE1->_pSysData[ii].w ) > 0.002 ){
bF1 += pdE1->_pSysData[ii].x;
if( fabsf( pdE1->_pSysData[ii].x ) > 0.001 ){
fprintf( stderr, "%4d %15.7e %15.7e %15.7e %15.7e %15.7e %15.7e %15.7e %15.7e\n", ii,
pdE1->_pSysData[ii].x, pdE1->_pSysData[ii].y, pdE1->_pSysData[ii].z, pdE1->_pSysData[ii].w,
pdE2->_pSysData[ii].x, pdE2->_pSysData[ii].y, pdE2->_pSysData[ii].z, pdE2->_pSysData[ii].w );
}
}
bF += pdE1->_pSysData[ii].x;
}
fprintf( stderr, "bSum Cud %6d %15.7e %15.7e\n", TARGET, bF, bF1 );
#endif
}
plugins/freeEnergy/platforms/cuda/src/kernels/kCalculateGBVISoftcoreBornSum.h
View file @ 032e18de
......@@ -37,9 +37,6 @@
#include "kCalculateGBVIAux.h"
#undef TARGET
//#define TARGET 39
__global__
#if (__CUDA_ARCH__ >= 200)
__launch_bounds__(GF1XX_NONBOND_THREADS_PER_BLOCK, 1)
......@@ -48,11 +45,7 @@ __launch_bounds__(GT2XX_NONBOND_THREADS_PER_BLOCK, 1)
#else
__launch_bounds__(G8X_NONBOND_THREADS_PER_BLOCK, 1)
#endif
#ifdef DEBUG
void METHOD_NAME(kCalculateGBVISoftcore, BornSum_kernel)(unsigned int* workUnit, float4* pdE1, float4* pdE2 )
#else
void METHOD_NAME(kCalculateGBVISoftcore, BornSum_kernel)(unsigned int* workUnit)
#endif
{
extern __shared__ Atom sA[];
......@@ -119,28 +112,6 @@ void METHOD_NAME(kCalculateGBVISoftcore, BornSum_kernel)(unsigned int* workUnit)
{
bSum += psA[j].bornRadiusScaleFactor*getGBVI_Volume( sqrt(r2), ar.x, psA[j].sr );
#ifdef DEBUG
int jIdx = j;
if( i == TARGET ){
int tjj = y+jIdx;
pdE1[tjj].x = psA[jIdx].bornRadiusScaleFactor*getGBVI_Volume( sqrt(r2), ar.x, psA[jIdx].sr );
pdE1[tjj].y = psA[jIdx].bornRadiusScaleFactor;
pdE1[tjj].z = ar.x;
pdE1[tjj].w = 1.0f;
pdE2[tjj].x = sqrt(r2);
pdE2[tjj].y = psA[jIdx].sr;
pdE2[tjj].z = ar.x;
pdE2[tjj].w = 1.0f;
}
if( (y+jIdx) == TARGET ){
int tjj = i;
pdE1[tjj].x = psA[jIdx].bornRadiusScaleFactor*getGBVI_Volume( sqrt(r2), ar.x, psA[jIdx].sr );
pdE1[tjj].y = psA[jIdx].bornRadiusScaleFactor;
pdE1[tjj].z = ar.x;
pdE1[tjj].w = -1.0f;
}
#endif
}
}
......@@ -214,32 +185,6 @@ pdE1[tjj].w = -1.0f;
apos.w += psA[tj].bornRadiusScaleFactor*getGBVI_Volume( r, ar.x, psA[tj].sr );
psA[tj].sum += ar.w*getGBVI_Volume( r, psA[tj].r, ar.y );
#ifdef DEBUG
int jIdx = tj;
if( i == TARGET ){
int tjj = y+jIdx;
pdE1[tjj].x = psA[jIdx].bornRadiusScaleFactor*getGBVI_Volume( r, ar.x, psA[jIdx].sr );
pdE1[tjj].y = psA[jIdx].bornRadiusScaleFactor;
pdE1[tjj].z = ar.x;
pdE1[tjj].w = 2.0f;
float R = ar.x;
float S = psA[tj].sr;
pdE2[tjj].x = getGBVI_L( r, (r + S), S );
pdE2[tjj].y = -getGBVI_L( r, (r - S), S );
pdE2[tjj].z = -getGBVI_L( r, R, S );
pdE2[tjj].w = (1.0f/(R*R*R));
}
if( (y+jIdx) == TARGET ){
int tjj = i;
pdE1[tjj].x = ar.w*getGBVI_Volume( r, psA[jIdx].r, ar.y );
pdE1[tjj].y = ar.w;
pdE1[tjj].z = psA[jIdx].r;
pdE1[tjj].w = -2.0f;
}
#endif
}
tj = (tj - 1) & (GRID - 1);
}
......
plugins/freeEnergy/platforms/cuda/src/kernels/kCalculateGBVISoftcoreForces2.cu
View file @ 032e18de
......@@ -68,66 +68,6 @@ void SetCalculateGBVISoftcoreForces2Sim( freeEnergyGpuContext gpu)
#include "kCalculateGBVIAux.h"
/**
* This file contains the kernel for evalauating the second stage of GBSA. It is included
* several times in kCalculateGBVIForces2.cu with different #defines to generate
* different versions of the kernels.
*/
__global__ void kCalculateGBVISoftcoreForces2a_kernel()
{
unsigned int pos = (blockIdx.x * blockDim.x + threadIdx.x);
if( pos >= cSim.atoms )return;
float4 apos = cSim.pPosq[pos];
float4 ar = cSim.pGBVIData[pos];
float fb = cSim.pBornForce[pos];
unsigned int posJ = 0;
float4 force;
force.x = force.y = force.z = force.w = 0.0f;
while ( posJ < cSim.atoms )
{
float4 aposJ = cSim.pPosq[posJ];
float4 arJ = cSim.pGBVIData[posJ];
float fbJ = cSim.pBornForce[posJ];
float dx = aposJ.x - apos.x;
float dy = aposJ.y - apos.y;
float dz = aposJ.z - apos.z;
float r2 = dx * dx + dy * dy + dz * dz;
float r = sqrt(r2);
float dE = getGBVI_dE2( r, ar.x, arJ.y, fb );
dE = r > 1.0e-08f ? dE : 0.0f;
//dx = dy = dz = 1.0f;
float d = dx*dE;
force.x -= d;
d = dy*dE;
force.y -= d;
d = dz*dE;
force.z -= d;
#if 1
dE = getGBVI_dE2( r, arJ.x, ar.y, fbJ );
dE = r > 1.0e-08f ? dE : 0.0f;
d = dx*dE;
force.x -= d;
d = dy*dE;
force.y -= d;
d = dz*dE;
force.z -= d;
#endif
posJ += 1;
}
// Write results
cSim.pForce4[pos] = force;
}
// Include versions of the kernels for N^2 calculations.
#define METHOD_NAME(a, b) a##N2##b
......@@ -180,46 +120,6 @@ void kCalculateGBVISoftcoreForces2( freeEnergyGpuContext freeEnergyGpu )
}
threadsPerBlock = threadsPerBlockPerMethod[methodIndex];
#ifdef DEBUG
fprintf( stderr,"kCalculateGBVISoftcoreForces2 nonbondedMethod=%d bornForce2_blocks=%u threadsPerBlock=%u shMem=%u\n",
freeEnergyGpu->freeEnergySim.nonbondedMethod,
gpu->sim.bornForce2_blocks, threadsPerBlock, (sizeof(Atom)+sizeof(float3))*threadsPerBlock ); fflush( stderr );
int psize = 64;
CUDAStream<float4>* pdE1 = new CUDAStream<float4>( psize, 1, "pdE");
CUDAStream<float4>* pdE2 = new CUDAStream<float4>( psize, 1, "pdE");
for( int ii = 0; ii < 32; ii++ ){
pdE1->_pSysData[ii].x = 0.0f;
pdE1->_pSysData[ii].y = 0.0f;
pdE1->_pSysData[ii].z = 0.0f;
pdE1->_pSysData[ii].w = 0.0f;
pdE2->_pSysData[ii].x = 0.0f;
pdE2->_pSysData[ii].y = 0.0f;
pdE2->_pSysData[ii].z = 0.0f;
pdE2->_pSysData[ii].w = 0.0f;
}
pdE1->Upload();
pdE2->Upload();
if (gpu->bOutputBufferPerWarp)
kCalculateGBVISoftcoreN2ByWarpForces2_kernel<<<gpu->sim.bornForce2_blocks, threadsPerBlock,
sizeof(Atom)*threadsPerBlock>>>(gpu->sim.pWorkUnit, gpu->sim.workUnits, pdE1->_pDevData, pdE2->_pDevData);
else
kCalculateGBVISoftcoreN2Forces2_kernel<<<gpu->sim.bornForce2_blocks, threadsPerBlock,
sizeof(Atom)*threadsPerBlock>>>(gpu->sim.pWorkUnit, gpu->sim.workUnits, pdE1->_pDevData, pdE2->_pDevData);
pdE1->Download();
pdE2->Download();
fprintf( stderr, "Pde\n" );
for( int ii = 0; ii < 32; ii++ ){
fprintf( stderr, "%4d %15.7e %15.7e %15.7e %15.7e %15.7e %15.7e %15.7e %15.7e\n", ii,
pdE1->_pSysData[ii].x, pdE1->_pSysData[ii].y, pdE1->_pSysData[ii].z, pdE1->_pSysData[ii].w,
pdE2->_pSysData[ii].x, pdE2->_pSysData[ii].y, pdE2->_pSysData[ii].z, pdE2->_pSysData[ii].w );
}
break;
#endif
switch (freeEnergyGpu->freeEnergySim.nonbondedMethod)
{
case FREE_ENERGY_NO_CUTOFF:
......
plugins/freeEnergy/platforms/cuda/src/kernels/kCalculateGBVISoftcoreForces2.h
View file @ 032e18de
......@@ -47,18 +47,14 @@ __launch_bounds__(G8X_BORNFORCE2_THREADS_PER_BLOCK, 1)
#endif
void METHOD_NAME(kCalculateGBVISoftcore, Forces2_kernel)(unsigned int* workUnit )
{
//METHOD_NAME(kCalculateGBVISoftcore, Forces2_kernel)(unsigned int* workUnit, float4* pdE1, float4* pdE2 )
extern __shared__ Atom sA[];
unsigned int totalWarps = gridDim.x*blockDim.x/GRID;
unsigned int warp = (blockIdx.x*blockDim.x+threadIdx.x)/GRID;
//unsigned int totalWarps = cSim.bornForce2_blocks*cSim.bornForce2_threads_per_block/GRID;
//unsigned int warp = (blockIdx.x*blockDim.x+threadIdx.x)/GRID;
unsigned int numWorkUnits = cSim.pInteractionCount[0];
unsigned int pos = warp*numWorkUnits/totalWarps;
unsigned int end = (warp+1)*numWorkUnits/totalWarps;
#ifdef USE_CUTOFF
//float3* tempBuffer = (float3*) &sA[cSim.bornForce2_threads_per_block];
float3* tempBuffer = (float3*) &sA[blockDim.x];
#endif
......@@ -126,19 +122,6 @@ void METHOD_NAME(kCalculateGBVISoftcore, Forces2_kernel)(unsigned int* workUnit
dE = 0.0f;
}
/*
if( i == TARGET ){
pdE1[x+j].x = dE;
pdE1[x+j].y = psA[j].bornRadiusScaleFactor;
pdE1[x+j].z = r;
pdE1[x+j].w = dE1;
}
if( (x+j) == TARGET ){
pdE2[i].x = dE;
pdE2[i].y = psA[j].bornRadiusScaleFactor;
pdE2[i].z = r;
pdE2[i].w = psA[j].sr-ar.x;
}*/
float d = dx * dE;
af.x -= d;
psA[j].fx += d;
......
plugins/freeEnergy/platforms/cuda/src/kernels/kCalculateLocalSoftcoreForces.cu
View file @ 032e18de
......@@ -29,7 +29,7 @@
#include <cudatypes.h>
#include "kSoftcoreLJ.h"
#define PARAMETER_PRINT 1
#define PARAMETER_PRINT 0
#define MAX_PARAMETER_PRINT 10
static __constant__ cudaGmxSimulation cSim;
......
plugins/freeEnergy/platforms/cuda/src/kernels/kCalculateNonbondedSoftcore.cu
View file @ 032e18de
......@@ -117,9 +117,6 @@ void kCalculateCDLJSoftcoreForces( freeEnergyGpuContext freeEnergyGpu )
{
gpuContext gpu = freeEnergyGpu->gpuContext;
// (void) fprintf( stderr,"kCalculateCDLJCutoffForces %d warp=%u nonbond_blocks=%u nonbond_threads_per_block=%u rfK=%15.7e rfC=%15.7e\n", freeEnergyGpu->freeEnergySim.nonbondedMethod,
// gpu->bOutputBufferPerWarp, gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block, gpu->sim.reactionFieldK, gpu->sim.reactionFieldC); fflush( stderr );
switch (freeEnergyGpu->freeEnergySim.nonbondedMethod)
{
case FREE_ENERGY_NO_CUTOFF:
......
plugins/freeEnergy/platforms/cuda/src/kernels/kCalculateObcGbsaSoftcoreBornSum.cu
View file @ 032e18de
......@@ -34,7 +34,6 @@
#define PARAMETER_PRINT 0
#define MAX_PARAMETER_PRINT 10
//#define DEBUG
struct Atom {
float x;
......@@ -211,7 +210,6 @@ __global__ void kReduceObcGbsaSoftcoreBornSum_kernel()
// Now calculate Born radius and OBC term.
sum *= 0.5f * atom.x;
sum *= gbsaSimDev.pNonPolarScalingFactors[pos];
float sum2 = sum * sum;
float sum3 = sum * sum2;
float tanhSum = tanh(cSim.alphaOBC * sum - cSim.betaOBC * sum2 + cSim.gammaOBC * sum3);
......@@ -377,31 +375,6 @@ void kCalculateObcGbsaSoftcoreBornSum( freeEnergyGpuContext freeEnergyGpu )
// printf("kCalculateObcGbsaSoftcoreBornSum\n");
gpuContext gpu = freeEnergyGpu->gpuContext;
#ifdef DEBUG
fprintf( stderr, "kCalculateObcGbsaSoftcoreBornSum cutoff=%15.7e\n", gpu->sim.nonbondedCutoffSqr );
int psize = gpu->sim.paddedNumberOfAtoms;
CUDAStream<float4>* pdE1 = new CUDAStream<float4>( psize, 1, "pdE");
CUDAStream<float4>* pdE2 = new CUDAStream<float4>( psize, 1, "pdE");
float bF;
float bF1;
showWorkUnitsFreeEnergy( freeEnergyGpu, 1 );
for( int ii = 0; ii < psize; ii++ ){
pdE1->_pSysData[ii].x = 0.0f;
pdE1->_pSysData[ii].y = 0.001f;
pdE1->_pSysData[ii].z = 0.001f;
pdE1->_pSysData[ii].w = 0.001f;
pdE2->_pSysData[ii].x = 0.001f;
pdE2->_pSysData[ii].y = 0.001f;
pdE2->_pSysData[ii].z = 0.001f;
pdE2->_pSysData[ii].w = 0.001f;
}
pdE1->Upload();
pdE2->Upload();
#endif
kClearObcGbsaSoftcoreBornSum(gpu);
LAUNCHERROR("kClearBornSum from kCalculateObcGbsaSoftcoreBornSum");
......@@ -409,22 +382,13 @@ pdE2->Upload();
{
case FREE_ENERGY_NO_CUTOFF:
#ifdef DEBUG
if (gpu->bOutputBufferPerWarp)
kCalculateObcGbsaSoftcoreN2ByWarpBornSum_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block,
sizeof(Atom)*gpu->sim.nonbond_threads_per_block>>>(gpu->sim.pWorkUnit, pdE1->_pDevData, pdE2->_pDevData);
else
kCalculateObcGbsaSoftcoreN2BornSum_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block,
sizeof(Atom)*gpu->sim.nonbond_threads_per_block>>>(gpu->sim.pWorkUnit, pdE1->_pDevData, pdE2->_pDevData);
#else
if (gpu->bOutputBufferPerWarp)
kCalculateObcGbsaSoftcoreN2ByWarpBornSum_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block,
sizeof(Atom)*gpu->sim.nonbond_threads_per_block>>>(gpu->sim.pWorkUnit);
else
kCalculateObcGbsaSoftcoreN2BornSum_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block,
sizeof(Atom)*gpu->sim.nonbond_threads_per_block>>>(gpu->sim.pWorkUnit);
#endif
break;
case FREE_ENERGY_CUTOFF:
......@@ -437,14 +401,6 @@ pdE2->Upload();
kFindInteractionsWithinBlocksCutoff_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block,
sizeof(unsigned int)*gpu->sim.nonbond_threads_per_block>>>(gpu->sim.pInteractingWorkUnit);
#ifdef DEBUG
if (gpu->bOutputBufferPerWarp)
kCalculateObcGbsaSoftcoreCutoffByWarpBornSum_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block,
(sizeof(Atom)+sizeof(float))*gpu->sim.nonbond_threads_per_block>>>(gpu->sim.pInteractingWorkUnit, pdE1->_pDevData, pdE2->_pDevData);
else
kCalculateObcGbsaSoftcoreCutoffBornSum_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block,
(sizeof(Atom)+sizeof(float))*gpu->sim.nonbond_threads_per_block>>>(gpu->sim.pInteractingWorkUnit, pdE1->_pDevData, pdE2->_pDevData);
#else
if (gpu->bOutputBufferPerWarp)
kCalculateObcGbsaSoftcoreCutoffByWarpBornSum_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block,
......@@ -452,7 +408,6 @@ pdE2->Upload();
else
kCalculateObcGbsaSoftcoreCutoffBornSum_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block,
(sizeof(Atom)+sizeof(float))*gpu->sim.nonbond_threads_per_block>>>(gpu->sim.pInteractingWorkUnit);
#endif
break;
......@@ -466,16 +421,6 @@ pdE2->Upload();
kFindInteractionsWithinBlocksPeriodic_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block,
sizeof(unsigned int)*gpu->sim.nonbond_threads_per_block>>>(gpu->sim.pInteractingWorkUnit);
#ifdef DEBUG
if (gpu->bOutputBufferPerWarp)
kCalculateObcGbsaSoftcorePeriodicByWarpBornSum_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block,
(sizeof(Atom)+sizeof(float))*gpu->sim.nonbond_threads_per_block>>>(gpu->sim.pInteractingWorkUnit, pdE1->_pDevData, pdE2->_pDevData);
else
kCalculateObcGbsaSoftcorePeriodicBornSum_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block,
(sizeof(Atom)+sizeof(float))*gpu->sim.nonbond_threads_per_block>>>(gpu->sim.pInteractingWorkUnit, pdE1->_pDevData, pdE2->_pDevData);
#else
if (gpu->bOutputBufferPerWarp)
kCalculateObcGbsaSoftcorePeriodicByWarpBornSum_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block,
(sizeof(Atom)+sizeof(float))*gpu->sim.nonbond_threads_per_block>>>(gpu->sim.pInteractingWorkUnit);
......@@ -483,7 +428,7 @@ pdE2->Upload();
kCalculateObcGbsaSoftcorePeriodicBornSum_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block,
(sizeof(Atom)+sizeof(float))*gpu->sim.nonbond_threads_per_block>>>(gpu->sim.pInteractingWorkUnit);
#endif
break;
default:
......@@ -492,23 +437,4 @@ pdE2->Upload();
}
LAUNCHERROR("kCalculateObcGbsaSoftcoreBornSum");
#ifdef DEBUG
pdE1->Download();
pdE2->Download();
//gpu->psBornRadii->Download();
//gpu->psObcData->Download();
fprintf( stderr, "bL Obc Cud\n" );
bF = 0.0;
bF1 = 0.0;
for( int ii = 0; ii < gpu->natoms; ii++ ){
bF1 += pdE1->_pSysData[ii].x;
fprintf( stderr, "%4d %15.7e %15.7e %15.7e %15.7e %15.7e %15.7e %15.7e %15.7e\n", ii,
pdE1->_pSysData[ii].x, pdE1->_pSysData[ii].y, pdE1->_pSysData[ii].z, pdE1->_pSysData[ii].w,
pdE2->_pSysData[ii].x, pdE2->_pSysData[ii].y, pdE2->_pSysData[ii].z, pdE2->_pSysData[ii].w );
bF += pdE1->_pSysData[ii].x;
}
fprintf( stderr, "bS Obc Cud %6d %15.7e %15.7e\n", TARGET, bF, bF1 );
#endif
}
plugins/freeEnergy/platforms/cuda/src/kernels/kCalculateObcGbsaSoftcoreBornSum.h
View file @ 032e18de
......@@ -41,11 +41,7 @@ __launch_bounds__(GT2XX_NONBOND_THREADS_PER_BLOCK, 1)
#else
__launch_bounds__(G8X_NONBOND_THREADS_PER_BLOCK, 1)
#endif
#ifdef DEBUG
void METHOD_NAME(kCalculateObcGbsaSoftcore, BornSum_kernel)(unsigned int* workUnit, float4* pdE1, float4* pdE2)
#else
void METHOD_NAME(kCalculateObcGbsaSoftcore, BornSum_kernel)(unsigned int* workUnit)
#endif
{
extern __shared__ Atom sA[];
unsigned int totalWarps = cSim.nonbond_blocks*cSim.nonbond_threads_per_block/GRID;
......@@ -135,32 +131,6 @@ void METHOD_NAME(kCalculateObcGbsaSoftcore, BornSum_kernel)(unsigned int* workUn
term += 2.0f * ((1.0f / ar.x) - l_ij);
}
apos.w += psA[j].polarScaleData*term;
#ifdef DEBUG
int jIdx = j;
if( i == TARGET ){
int tjj = y+jIdx;
pdE1[tjj].x = term;
pdE1[tjj].y = r;
pdE1[tjj].z = ar.x;
pdE1[tjj].w = 1.0f;
pdE2[tjj].x = r;
pdE2[tjj].y = l_ij;
pdE2[tjj].z = rj;
pdE2[tjj].w = 1.0f;
}
/*
if( (y+jIdx) == TARGET ){
int tjj = i;
pdE2[tjj].x = sum;
pdE2[tjj].y = psA[jIdx].polarScaleData;
pdE2[tjj].z = ar.x;
pdE2[tjj].w = -1.0f;
} */
#endif
}
}
}
......@@ -254,37 +224,6 @@ pdE2[tjj].w = -1.0f;
}
apos.w += (scale*term);
#ifdef DEBUG
int jIdx = tj;
if( i == TARGET ){
int tjj = y+jIdx;
pdE1[tjj].x = term;
pdE1[tjj].y = r;
pdE1[tjj].z = ar.x;
pdE1[tjj].w = 2.0f;
/*
pdE2[tjj].x = r;
pdE2[tjj].y = l_ij;
pdE2[tjj].z = rj;
pdE2[tjj].w = 2.0f;
*/
}
if( (y+jIdx) == TARGET ){
int tjj = i;
/*
pdE1[tjj].x = term;
pdE1[tjj].y = r;
pdE1[tjj].z = ar.x;
pdE1[tjj].w = -2.0f;
*/
pdE2[tjj].x = term;
pdE2[tjj].y = r;
pdE2[tjj].z = ar.x;
pdE2[tjj].w = -2.0f;
}
#endif
}
float rScaledRadiusI = r + ar.y;
if (psA[tj].r < rScaledRadiusI)
......@@ -307,25 +246,6 @@ pdE2[tjj].w = -2.0f;
}
psA[tj].sum += polarScaleDataI*term;
#ifdef DEBUG
int jIdx = tj;
if( i == TARGET ){
int tjj = y+jIdx;
pdE1[tjj].x = term;
pdE1[tjj].y = r;
pdE1[tjj].z = ar.x;
pdE1[tjj].w = 3.0f;
}
if( (y+jIdx) == TARGET ){
int tjj = i;
pdE2[tjj].x = term;
pdE2[tjj].y = r;
pdE2[tjj].z = ar.x;
pdE2[tjj].w = -3.0f;
}
#endif
}
}
tj = (tj - 1) & (GRID - 1);
......
plugins/freeEnergy/platforms/cuda/src/kernels/kCalculateObcGbsaSoftcoreForces2.h
View file @ 032e18de
......@@ -44,14 +44,10 @@ void METHOD_NAME(kCalculateObcGbsaSoftcore, Forces2_kernel)(unsigned int* workUn
unsigned int totalWarps = gridDim.x*blockDim.x/GRID;
unsigned int warp = (blockIdx.x*blockDim.x+threadIdx.x)/GRID;
//unsigned int totalWarps = cSim.bornForce2_blocks*cSim.bornForce2_threads_per_block/GRID;
//unsigned int warp = (blockIdx.x*blockDim.x+threadIdx.x)/GRID;
unsigned int numWorkUnits = cSim.pInteractionCount[0];
unsigned int pos = warp*numWorkUnits/totalWarps;
unsigned int end = (warp+1)*numWorkUnits/totalWarps;
#ifdef USE_CUTOFF
//float3* tempBuffer = (float3*) &sA[cSim.bornForce2_threads_per_block];
float3* tempBuffer = (float3*) &sA[blockDim.x];
#endif
......@@ -124,7 +120,7 @@ void METHOD_NAME(kCalculateObcGbsaSoftcore, Forces2_kernel)(unsigned int* workUn
// Born Forces term
float term = 0.125f * (1.000f + psA[j].sr * psA[j].sr * r2Inverse) * t3 +
0.250f * t1 * r2Inverse;
term *= psA[j].npScale*nonPolarScaleDataI;
term *= psA[j].npScale;
float dE = fb * term;
#if defined USE_CUTOFF
......@@ -231,7 +227,7 @@ void METHOD_NAME(kCalculateObcGbsaSoftcore, Forces2_kernel)(unsigned int* workUn
float term = 0.125f *
(1.000f + psA[tj].sr * psA[tj].sr * r2Inverse) * t3J +
0.250f * t1J * r2Inverse;
term *= psA[tj].npScale*nonPolarScaleDataI;
term *= psA[tj].npScale;
float dE = fb * term;
#if defined USE_CUTOFF
......@@ -257,7 +253,7 @@ void METHOD_NAME(kCalculateObcGbsaSoftcore, Forces2_kernel)(unsigned int* workUn
term = 0.125f *
(1.000f + sr2 * r2Inverse) * t3I +
0.250f * t1I * r2Inverse;
term *= psA[tj].npScale*nonPolarScaleDataI;
term *= nonPolarScaleDataI;
dE = psA[tj].fb * term;
float rj = psA[tj].r;
......@@ -293,11 +289,11 @@ void METHOD_NAME(kCalculateObcGbsaSoftcore, Forces2_kernel)(unsigned int* workUn
float dx = psA[j].x - apos.x;
float dy = psA[j].y - apos.y;
float dz = psA[j].z - apos.z;
#ifdef USE_PERIODIC
#ifdef USE_PERIODIC
dx -= floor(dx/cSim.periodicBoxSizeX+0.5f)*cSim.periodicBoxSizeX;
dy -= floor(dy/cSim.periodicBoxSizeY+0.5f)*cSim.periodicBoxSizeY;
dz -= floor(dz/cSim.periodicBoxSizeZ+0.5f)*cSim.periodicBoxSizeZ;
#endif
#endif
float r2 = dx * dx + dy * dy + dz * dz;
float r = sqrt(r2);
......@@ -327,16 +323,16 @@ void METHOD_NAME(kCalculateObcGbsaSoftcore, Forces2_kernel)(unsigned int* workUn
float term = 0.125f *
(1.000f + psA[j].sr * psA[j].sr * r2Inverse) * t3J +
0.250f * t1J * r2Inverse;
term *= psA[j].npScale*nonPolarScaleDataI;
term *= psA[j].npScale;
float dE = fb * term;
#if defined USE_PERIODIC
#if defined USE_PERIODIC
if (a.x >= rScaledRadiusJ || i >= cSim.atoms || y+j >= cSim.atoms || r2 > cSim.nonbondedCutoffSqr)
#elif defined USE_CUTOFF
#elif defined USE_CUTOFF
if (a.x >= rScaledRadiusJ || r2 > cSim.nonbondedCutoffSqr)
#else
#else
if (a.x >= rScaledRadiusJ)
#endif
#endif
{
dE = 0.0f;
}
......@@ -355,17 +351,17 @@ void METHOD_NAME(kCalculateObcGbsaSoftcore, Forces2_kernel)(unsigned int* workUn
term = 0.125f *
(1.000f + sr2 * r2Inverse) * t3I +
0.250f * t1I * r2Inverse;
term *= psA[j].npScale*nonPolarScaleDataI;
term *= nonPolarScaleDataI;
dE = psA[j].fb * term;
float rj = psA[j].r;
#ifdef USE_PERIODIC
#ifdef USE_PERIODIC
if (rj >= rScaledRadiusI || i >= cSim.atoms || y+j >= cSim.atoms || r2 > cSim.nonbondedCutoffSqr)
#elif defined USE_CUTOFF
#elif defined USE_CUTOFF
if (rj >= rScaledRadiusI || r2 > cSim.nonbondedCutoffSqr)
#else
#else
if (rj >= rScaledRadiusI)
#endif
#endif
{
dE = 0.0f;
}
......
plugins/freeEnergy/platforms/cuda/tests/CMakeLists.txt
View file @ 032e18de
......@@ -15,14 +15,12 @@ SET( INCLUDE_SERIALIZATION FALSE )
#SET( INCLUDE_SERIALIZATION TRUE )
SET( SHARED_OPENMM_TARGET OpenMMFreeEnergy)
SET( STATIC_OPENMM_TARGET OpenMMFreeEnergy_static)
SET( SHARED_CUDA_TARGET OpenMMCuda)
SET( STATIC_CUDA_TARGET OpenMMCuda_static)
IF( INCLUDE_SERIALIZATION )
INCLUDE_DIRECTORIES(${OPENMM_DIR}/serialization/include)
SET( SHARED_OPENMM_SERIALIZATION OpenMMSerialization )
SET( SHARED_FREE_ENERGY_SERIALIZATION FreeEnergySerialization )
SET( SHARED_OPENMM_SERIALIZATION "OpenMMSerialization" )
SET( SHARED_FREE_ENERGY_SERIALIZATION "FreeEnergySerialization" )
ENDIF( INCLUDE_SERIALIZATION )
IF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
......@@ -32,8 +30,6 @@ IF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
SET(SHARED_OPENMM_SERIALIZATION ${SHARED_OPENMM_SERIALIZATION}_d)
SET(SHARED_FREE_ENERGY_SERIALIZATION ${SHARED_FREE_ENERGY_SERIALIZATION}_d)
ENDIF( INCLUDE_SERIALIZATION )
SET(STATIC_CUDA_TARGET ${STATIC_CUDA_TARGET}_d)
SET(STATIC_OPENMM_TARGET ${STATIC_OPENMM_TARGET}_d)
ENDIF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
# Automatically create tests using files named "Test*.cpp"
......@@ -54,13 +50,52 @@ FOREACH(TEST_PROG ${TEST_PROGS})
SET(DEFINE_STRING "${DEFINE_STRING} -DOPENMM_SERIALIZE")
ENDIF( INCLUDE_SERIALIZATION )
IF( ${TEST_ROOT} STREQUAL "TestCudaOBCSoftcoreForce" )
SET(DEFINE_STRING "${DEFINE_STRING} -DIMPLICIT_SOLVENT=1")
ENDIF( ${TEST_ROOT} STREQUAL "TestCudaOBCSoftcoreForce" )
IF( ${TEST_ROOT} STREQUAL "TestCudaGBVISoftcoreForce2" )
# serialize
SET(DEFINE_STRING "-DTEST_PLATFORM=0 ")
IF( INCLUDE_SERIALIZATION )
SET(DEFINE_STRING "${DEFINE_STRING} -DOPENMM_SERIALIZE ")
SET(TARGET_LINK_LIBRARIES_STRING "${SHARED_TARGET} ${SHARED_OPENMM_SERIALIZATION}")
ELSE( INCLUDE_SERIALIZATION )
SET(TARGET_LINK_LIBRARIES_STRING "${SHARED_TARGET}")
ENDIF( INCLUDE_SERIALIZATION )
# obc
SET(OBC_DEFINE_STRING "${DEFINE_STRING} -DIMPLICIT_SOLVENT=1")
SET(OBC_TEST "TestCudaGBSAOBCSoftcoreForce2")
CUDA_ADD_EXECUTABLE(${OBC_TEST} ${TEST_PROG})
IF( INCLUDE_SERIALIZATION )
TARGET_LINK_LIBRARIES(${OBC_TEST} ${SHARED_TARGET} ${SHARED_OPENMM_SERIALIZATION} )
ELSE( INCLUDE_SERIALIZATION )
TARGET_LINK_LIBRARIES(${OBC_TEST} ${SHARED_TARGET})
ENDIF( INCLUDE_SERIALIZATION )
SET_TARGET_PROPERTIES(${OBC_TEST} PROPERTIES COMPILE_FLAGS ${OBC_DEFINE_STRING} )
ADD_TEST(${OBC_TEST} ${EXECUTABLE_OUTPUT_PATH}/${OBC_TEST})
# nonbond
SET(NONBOND_DEFINE_STRING "${DEFINE_STRING} -DIMPLICIT_SOLVENT=0")
SET(NONBOND_TEST "TestCudaNonbondSoftcoreForce2")
CUDA_ADD_EXECUTABLE(${NONBOND_TEST} ${TEST_PROG})
IF( INCLUDE_SERIALIZATION )
TARGET_LINK_LIBRARIES(${NONBOND_TEST} ${SHARED_TARGET} ${SHARED_OPENMM_SERIALIZATION} )
ELSE( INCLUDE_SERIALIZATION )
TARGET_LINK_LIBRARIES(${NONBOND_TEST} ${SHARED_TARGET})
ENDIF( INCLUDE_SERIALIZATION )
SET_TARGET_PROPERTIES(${NONBOND_TEST} PROPERTIES COMPILE_FLAGS ${NONBOND_DEFINE_STRING} )
ADD_TEST(${NONBOND_TEST} ${EXECUTABLE_OUTPUT_PATH}/${NONBOND_TEST})
# gbvi
IF( ${TEST_ROOT} STREQUAL "TestCudaGBVISoftcoreForce" )
SET(DEFINE_STRING "${DEFINE_STRING} -DIMPLICIT_SOLVENT=2")
ENDIF( ${TEST_ROOT} STREQUAL "TestCudaGBVISoftcoreForce" )
SET_TARGET_PROPERTIES(${TEST_ROOT} PROPERTIES COMPILE_FLAGS ${DEFINE_STRING} )
ENDIF( ${TEST_ROOT} STREQUAL "TestCudaGBVISoftcoreForce2" )
#MESSAGE( "${TEST_ROOT} ${DEFINE_STRING}" )
SET_TARGET_PROPERTIES(${TEST_ROOT} PROPERTIES COMPILE_FLAGS ${DEFINE_STRING} )
......
plugins/freeEnergy/platforms/cuda/tests/TestCudaGBVISoftcoreForce.cpp
View file @ 032e18de
......@@ -33,27 +33,44 @@
* This tests the reference implementation of GBVIForce.
*/
#include "TestCudaSoftcoreForce.h"
//#define USE_SOFTCORE
//#define IMPLICIT_SOLVENT GBVI
//#define IMPLICIT_SOLVENT OBC
#define OBC_FLAG 1
#define GBVI_FLAG 2
#include "../../../tests/AssertionUtilities.h"
#include "openmm/System.h"
#include "openmm/Context.h"
#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
#include "OpenMM.h"
#include "openmm/GBVIForce.h"
#include "openmm/GBSAOBCForce.h"
#include "openmm/NonbondedForce.h"
#ifdef USE_SOFTCORE
#include "openmm/GBVISoftcoreForce.h"
#include "openmm/GBSAOBCSoftcoreForce.h"
#include "openmm/NonbondedSoftcoreForce.h"
#endif
#include "OpenMMFreeEnergy.h"
#include "openmm/freeEnergyKernels.h"
#include "sfmt/SFMT.h"
#include "openmm/VerletIntegrator.h"
#include <iostream>
#include <vector>
#include <algorithm>
#include <cstdio>
#include <typeinfo>
#include <iomanip>
typedef std::vector<double> DoubleVector;
typedef DoubleVector::iterator DoubleVectorI;
typedef DoubleVector::const_iterator DoubleVectorCI;
typedef std::vector<DoubleVector> DoubleVectorVector;
typedef std::pair<int, double> IntDoublePair;
typedef std::vector<IntDoublePair> IntDoublePairVector;
typedef IntDoublePairVector::iterator IntDoublePairVectorI;
typedef IntDoublePairVector::const_iterator IntDoublePairVectorCI;
extern "C" void registerFreeEnergyCudaKernelFactories();
using namespace OpenMM;
using namespace std;
void testSingleParticle( FILE* log ) {
System system;
......@@ -96,6 +113,141 @@ void testSingleParticle( FILE* log ) {
ASSERT_EQUAL_TOL((bornEnergy+nonpolarEnergy), state.getPotentialEnergy(), 0.01);
}
/**
* Predicate for sorting <int,double> pair
*
* @param d1 first IntDoublePair to compare
* @param d2 second IntDoublePair to compare
*
*/
bool TestIntDoublePair( const IntDoublePair& d1, const IntDoublePair& d2 ){
return d1.second < d2.second;
}
/**---------------------------------------------------------------------------------------
*
* Get relative difference between two forces
*
* @param f1 force1
* @param f2 force2
* @param forceNorm1 output norm of force1
* @param forceNorm2 output norm of force2
* @param relativeDiff output relative difference between force norms
* @param log if set, output forces
*
*
--------------------------------------------------------------------------------------- */
static void getForceRelativeDifference( const Vec3& f1, const Vec3& f2, double& forceNorm1, double& forceNorm2,
double& relativeDiff, FILE* log ) {
double diff = (f1[0] - f2[0])*(f1[0] - f2[0]) +
(f1[1] - f2[1])*(f1[1] - f2[1]) +
(f1[2] - f2[2])*(f1[2] - f2[2]);
forceNorm1 = sqrt( f1[0]*f1[0] + f1[1]*f1[1] + f1[2]*f1[2] );
forceNorm2 = sqrt( f2[0]*f2[0] + f2[1]*f2[1] + f2[2]*f2[2] );
if( forceNorm1 > 0.0 || forceNorm2 > 0.0 ){
relativeDiff = 2.0*sqrt( diff )/(forceNorm1+forceNorm2);
} else {
relativeDiff = 0.0;
}
return;
}
/**---------------------------------------------------------------------------------------
*
* Compare forces from two states
*
* @param state1 state1
* @param state2 state2
* @param relativeTolerance relative tolerance
* @param log if set, output forces
*
* @return number of entries with relative difference > tolerance
*
--------------------------------------------------------------------------------------- */
int compareForcesOfTwoStates( State& state1, State& state2, double relativeTolerance,
DoubleVector& stats, FILE* log ) {
int error = 0;
vector<Vec3> force1 = state1.getForces();
vector<Vec3> force2 = state2.getForces();
double averageRelativeDifference = 0.0;
double count = 0.0;
DoubleVector medians1( force1.size() );
DoubleVector medians2( force1.size() );
IntDoublePairVector relativeDifferences;
for( unsigned int ii = 0; ii < force1.size(); ii++ ){
double forceNorm1;
double forceNorm2;
double relativeDiff;
getForceRelativeDifference( force1[ii], force2[ii], forceNorm1, forceNorm2, relativeDiff, log );
medians1[ii] = forceNorm1;
medians2[ii] = forceNorm2;
relativeDifferences.push_back( IntDoublePair(ii, relativeDiff ) );
averageRelativeDifference += relativeDiff;
count += 1.0;
if( relativeDiff > relativeTolerance ){
error++;
}
if( log ){
(void) fprintf( log, "F %6u %15.7e [%15.7e %15.7e %15.7e] [%15.7e %15.7e %15.7e] %15.7e %15.7e %s\n", static_cast<unsigned int>(ii),
relativeDiff, force1[ii][0], force1[ii][1], force1[ii][2], force2[ii][0], force2[ii][1], force2[ii][2],
forceNorm1, forceNorm2, (relativeDiff < relativeTolerance ? "":"XXXXXX") );
}
}
// sort relative differences
std::sort( relativeDifferences.begin(), relativeDifferences.end(), TestIntDoublePair );
if( log ){
(void) fprintf( log, "\nEntries w/ largest relative differences.\n" );
for( unsigned int ii = relativeDifferences.size()-1; ii >= relativeDifferences.size()-10 && ii >= 0; ii-- ){
double forceNorm1;
double forceNorm2;
double relativeDiff;
int index = relativeDifferences[ii].first;
getForceRelativeDifference( force1[index], force2[index], forceNorm1, forceNorm2, relativeDiff, log );
(void) fprintf( log, "Fs %6u %15.7e [%15.7e %15.7e %15.7e] [%15.7e %15.7e %15.7e] %15.7e %15.7e %s\n",
static_cast<unsigned int>(index), relativeDiff,
force1[index][0], force1[index][1], force1[index][2],
force2[index][0], force2[index][1], force2[index][2],
forceNorm1, forceNorm2, (relativeDiff < relativeTolerance ? "":"XXXXXX") );
}
}
if( count > 0.0 ){
averageRelativeDifference /= count;
}
std::sort( medians1.begin(), medians1.end() );
std::sort( medians2.begin(), medians2.end() );
double median1 = medians1[medians1.size()/2];
double median2 = medians2[medians2.size()/2];
stats.resize( 4 );
stats[0] = averageRelativeDifference;
IntDoublePair pair = relativeDifferences[relativeDifferences.size()-1];
stats[1] = pair.second;
stats[2] = static_cast<double>(pair.first);
stats[3] = median1 < median2 ? median1 : median2;
return error;
}
void testEnergyEthaneSwitchingFunction( int useSwitchingFunction, FILE* log ) {
std::string methodName = "testEnergyEthaneSwitchingFunction";
......@@ -250,7 +402,7 @@ void testEnergyEthaneSwitchingFunction( int useSwitchingFunction, FILE* log ) {
// Take a small step in the direction of the energy gradient.
DoubleVector stats;
std::vector<double> stats;
if( compareForcesOfTwoStates( state, referenceState, 0.001, stats, log ) ){
ASSERT_EQUAL_TOL(0.0, 1.0, 0.01)
}
......@@ -333,555 +485,20 @@ void testEnergyEthaneSwitchingFunction( int useSwitchingFunction, FILE* log ) {
}
}
static GBVISoftcoreForce* copyGbviSoftcoreForce( const GBVISoftcoreForce& gbviSoftcoreForce ){
GBVISoftcoreForce* copyGbviSoftcoreForce = new GBVISoftcoreForce(gbviSoftcoreForce);
/*
GBVISoftcoreForce* copyGbviSoftcoreForce = new GBVISoftcoreForce();
copyGbviSoftcoreForce->setNonbondedMethod( gbviSoftcoreForce.getNonbondedMethod() );
copyGbviSoftcoreForce->setCutoffDistance( gbviSoftcoreForce.getCutoffDistance() );
copyGbviSoftcoreForce->setSolventDielectric( gbviSoftcoreForce.getSolventDielectric() );
copyGbviSoftcoreForce->setSoluteDielectric( gbviSoftcoreForce.getSoluteDielectric() );
copyGbviSoftcoreForce->setBornRadiusScalingMethod( gbviSoftcoreForce.getBornRadiusScalingMethod() );
copyGbviSoftcoreForce->setQuinticLowerLimitFactor( gbviSoftcoreForce.getQuinticLowerLimitFactor() );
copyGbviSoftcoreForce->setQuinticUpperBornRadiusLimit( gbviSoftcoreForce.getQuinticUpperBornRadiusLimit() );
// particle parameters
for( unsigned int ii = 0; ii < gbviSoftcoreForce.getNumParticles(); ii++ ){
double charge;
double sigma;
double gamma;
double softcoreLJLambda;
gbviSoftcoreForce.getParticleParameters(ii, charge, sigma, gamma, softcoreLJLambda);
copyGbviSoftcoreForce->addParticle( charge, sigma, gamma, softcoreLJLambda);
}
// bonds
for( unsigned int ii = 0; ii < gbviSoftcoreForce.getNumBonds(); ii++ ){
int particle1, particle2;
double distance;
gbviSoftcoreForce.getBondParameters( ii, particle1, particle2, distance);
copyGbviSoftcoreForce->addBond( particle1, particle2, distance );
}
*/
return copyGbviSoftcoreForce;
}
static GBVIForce* copyGbviForce( const GBVIForce& gbviForce ){
return new GBVIForce(gbviForce);
}
static GBSAOBCSoftcoreForce* copyGBSAOBCSoftcoreForce( const GBSAOBCSoftcoreForce& gbviSoftcoreForce ){
return new GBSAOBCSoftcoreForce(gbviSoftcoreForce);
}
static GBSAOBCForce* copyGbsaObcForce( const GBSAOBCForce& gbviForce ){
return new GBSAOBCForce(gbviForce);
}
void testGbviSoftcore( MapStringToDouble& inputArgumentMap, FILE* log ){
double lambda1 = 1.0;
double lambda2 = 1.0;
int nonbondedMethod = 0;
int numMolecules = 1;
int numParticlesPerMolecule = 2;
int useQuinticSpline = 1;
int applyAssert = 1;
int positionPlacementMethod = 0;
int serialize = 0;
double boxSize = 10.0;
double relativeTolerance = 1.0e-04;
setDoubleFromMapStringToDouble( inputArgumentMap, "lambda1", lambda1 );
setDoubleFromMapStringToDouble( inputArgumentMap, "lambda2", lambda2 );
setDoubleFromMapStringToDouble( inputArgumentMap, "boxSize", boxSize );
double cutoffDistance = boxSize*0.4;;
setDoubleFromMapStringToDouble( inputArgumentMap, "cutoffDistance", cutoffDistance);
setDoubleFromMapStringToDouble( inputArgumentMap, "relativeTolerance", relativeTolerance );
setIntFromMapStringToDouble( inputArgumentMap, "positionPlacementMethod", positionPlacementMethod ) ;
setIntFromMapStringToDouble( inputArgumentMap, "nonbondedMethod", nonbondedMethod );
setIntFromMapStringToDouble( inputArgumentMap, "numMolecules", numMolecules );
setIntFromMapStringToDouble( inputArgumentMap, "numParticlesPerMolecule", numParticlesPerMolecule );
setIntFromMapStringToDouble( inputArgumentMap, "serialize", serialize );
if( nonbondedMethod == 2 && cutoffDistance > boxSize*0.5 ){
cutoffDistance = boxSize*0.5;
}
int numParticles = numMolecules*numParticlesPerMolecule;
int includeGbvi = 1;
double reactionFieldDielectric = 80.0;
if( log ){
double particleDensity = static_cast<double>(numParticles)/(boxSize*boxSize*boxSize);
double particleCube = pow( particleDensity, (-1.0/3.0) );
(void) fprintf( log, "\n--------------------------------------------------------------------------------------\n" );
(void) fprintf( log, "Input arguments\n" );
(void) fflush( log );
(void) fprintf( log, " includeGbvi %d\n", includeGbvi );
(void) fprintf( log, " nonbondedMethod %d\n", nonbondedMethod );
(void) fprintf( log, " numParticles %d\n", numParticles );
(void) fprintf( log, " numMolecules %d\n", numMolecules );
(void) fprintf( log, " numParticlesPerMolecule %d\n", numParticlesPerMolecule );
(void) fprintf( log, " useQuinticSpline %d\n", useQuinticSpline );
(void) fprintf( log, " positionPlacementMethod %d\n", positionPlacementMethod);
#ifdef USE_SOFTCORE
(void) fprintf( log, " lambda1 %8.3f\n", lambda1 );
(void) fprintf( log, " lambda2 %8.3f\n", lambda2 );
#endif
(void) fprintf( log, " boxSize %8.3f\n", boxSize );
(void) fprintf( log, " cutoffDistance %8.3f\n", cutoffDistance );
(void) fprintf( log, " reactionFieldDielectric %8.3f\n", reactionFieldDielectric );
(void) fprintf( log, " relativeTolerance %8.1e\n", relativeTolerance );
(void) fprintf( log, " particleDensity %8.2e\n", particleDensity );
(void) fprintf( log, " particleCube %8.2e\n", particleCube );
}
// Create two systems: one with GbviSoftcoreForce NonbondedSoftcoreForce forces, and one using a CustomNonbondedForce, CustomGBVI force to implement the same interaction.
System standardSystem;
for (int i = 0; i < numParticles; i++) {
standardSystem.addParticle(1.0);
}
standardSystem.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
#ifdef USE_SOFTCORE
NonbondedSoftcoreForce* nonbondedSoftcoreForce = new NonbondedSoftcoreForce();
if( nonbondedMethod == NoCutoff ){
nonbondedSoftcoreForce->setNonbondedMethod( NonbondedSoftcoreForce::NoCutoff );
} else {
if( nonbondedMethod == CutoffNonPeriodic ){
nonbondedSoftcoreForce->setNonbondedMethod( NonbondedSoftcoreForce::CutoffNonPeriodic );
} else {
nonbondedSoftcoreForce->setNonbondedMethod( NonbondedSoftcoreForce::CutoffPeriodic );
}
}
#else
NonbondedForce* nonbondedSoftcoreForce = new NonbondedForce();
if( nonbondedMethod == NoCutoff ){
nonbondedSoftcoreForce->setNonbondedMethod( NonbondedForce::NoCutoff );
} else {
if( nonbondedMethod == CutoffNonPeriodic ){
nonbondedSoftcoreForce->setNonbondedMethod( NonbondedForce::CutoffNonPeriodic );
} else {
nonbondedSoftcoreForce->setNonbondedMethod( NonbondedForce::CutoffPeriodic );
}
}
#endif
nonbondedSoftcoreForce->setCutoffDistance( cutoffDistance );
nonbondedSoftcoreForce->setReactionFieldDielectric( reactionFieldDielectric );
#ifdef USE_SOFTCORE
#if IMPLICIT_SOLVENT == GBVI_FLAG
GBVISoftcoreForce* gbviSoftcoreForce = new GBVISoftcoreForce();
if( nonbondedMethod == NoCutoff ){
gbviSoftcoreForce->setNonbondedMethod( GBVISoftcoreForce::NoCutoff );
} else {
if( nonbondedMethod == CutoffNonPeriodic ){
gbviSoftcoreForce->setNonbondedMethod( GBVISoftcoreForce::CutoffNonPeriodic );
} else {
gbviSoftcoreForce->setNonbondedMethod( GBVISoftcoreForce::CutoffPeriodic );
}
}
#else
GBSAOBCSoftcoreForce* gbviSoftcoreForce = new GBSAOBCSoftcoreForce();
if( nonbondedMethod == NoCutoff ){
gbviSoftcoreForce->setNonbondedMethod( GBSAOBCSoftcoreForce::NoCutoff );
} else {
if( nonbondedMethod == CutoffNonPeriodic ){
gbviSoftcoreForce->setNonbondedMethod( GBSAOBCSoftcoreForce::CutoffNonPeriodic );
} else {
gbviSoftcoreForce->setNonbondedMethod( GBSAOBCSoftcoreForce::CutoffPeriodic );
}
}
#endif
#else
#if IMPLICIT_SOLVENT == GBVI_FLAG
GBVIForce* gbviSoftcoreForce = new GBVIForce();
if( nonbondedMethod == NoCutoff ){
gbviSoftcoreForce->setNonbondedMethod( GBVIForce::NoCutoff );
} else {
if( nonbondedMethod == CutoffNonPeriodic ){
gbviSoftcoreForce->setNonbondedMethod( GBVIForce::CutoffNonPeriodic );
} else {
gbviSoftcoreForce->setNonbondedMethod( GBVIForce::CutoffPeriodic );
}
}
#else
GBSAOBCForce* gbviSoftcoreForce = new GBSAOBCForce();
if( nonbondedMethod == NoCutoff ){
gbviSoftcoreForce->setNonbondedMethod( GBSAOBCForce::NoCutoff );
} else {
if( nonbondedMethod == CutoffNonPeriodic ){
gbviSoftcoreForce->setNonbondedMethod( GBSAOBCForce::CutoffNonPeriodic );
} else {
gbviSoftcoreForce->setNonbondedMethod( GBSAOBCForce::CutoffPeriodic );
}
}
#endif
#endif
#if IMPLICIT_SOLVENT == GBVI_FLAG
#ifdef USE_SOFTCORE
if( useQuinticSpline ){
gbviSoftcoreForce->setBornRadiusScalingMethod( GBVISoftcoreForce::QuinticSpline );
} else {
gbviSoftcoreForce->setBornRadiusScalingMethod( GBVISoftcoreForce::NoScaling );
}
#else
if( useQuinticSpline ){
gbviSoftcoreForce->setBornRadiusScalingMethod( GBVIForce::QuinticSpline );
} else {
gbviSoftcoreForce->setBornRadiusScalingMethod( GBVIForce::NoScaling );
}
#endif
#endif
gbviSoftcoreForce->setSolventDielectric( 78.3 );
//gbviSoftcoreForce->setSolventDielectric( 1.0e+10 );
//gbviSoftcoreForce->setSolventDielectric( 1.0 );
gbviSoftcoreForce->setSoluteDielectric( 1.0 );
gbviSoftcoreForce->setCutoffDistance( nonbondedSoftcoreForce->getCutoffDistance( ) );
std::vector<Vec3> positions(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
PositionGenerator positionGenerator( numMolecules, numParticlesPerMolecule, boxSize );
if( log ){
positionGenerator.setLog( log );
}
if( positionPlacementMethod == 1 ){
positionGenerator.setPositions( PositionGenerator::SimpleGrid, sfmt, positions );
} else {
positionGenerator.setBondDistance( 0.3 );
positionGenerator.setPositions( PositionGenerator::Random, sfmt, positions );
}
// show info on particle positions
if( log ){
Vec3 box[2];
positionGenerator.getEnclosingBox( positions, box );
(void) fprintf( log, "Enclosing Box (in A): [%15.7e %15.7e] [%15.7e %15.7e] [%15.7e %15.7e] [%15.7e %15.7e %15.7e]\n",
box[0][0], box[1][0], box[0][1], box[1][1], box[0][2], box[1][2],
(box[1][0] - box[0][0]), (box[1][1] - box[0][1]), (box[1][2] - box[0][2]) );
int showIndex = 5;
int periodicBoundaryConditions = (nonbondedMethod == 2) ? 1 : 0;
IntVector positionIndexVector;
positionIndexVector.push_back( 0 );
positionIndexVector.push_back( static_cast<int>(positions.size())-1 );
//positionIndexVector.push_back( 542 );
for( unsigned int ii = 0; ii < positionIndexVector.size(); ii++ ){
if( positionIndexVector[ii] < positions.size() ){
int positionIndex = positionIndexVector[ii];
IntDoublePairVector sortVector;
positionGenerator.getSortedDistances( periodicBoundaryConditions, positionIndex, positions, sortVector );
(void) fprintf( log, "Min/max distance from %6d:\n ", positionIndex );
for( unsigned int jj = 0; jj < sortVector.size() && jj < showIndex; jj++ ){
IntDoublePair pair = sortVector[jj];
(void) fprintf( log, "[%6d %15.7e] ", pair.first, pair.second);
}
(void) fprintf( log, "\n " );
for( unsigned int jj = (sortVector.size() - showIndex); jj < sortVector.size() && jj >= 0; jj++ ){
IntDoublePair pair = sortVector[jj];
(void) fprintf( log, "[%6d %15.7e] ", pair.first, pair.second);
}
(void) fprintf( log, "\n" );
}
}
IntIntPairVector pairs;
pairs.push_back( IntIntPair( 732, 0 ) );
pairs.push_back( IntIntPair( 732, 1 ) );
pairs.push_back( IntIntPair( 732, 2 ) );
pairs.push_back( IntIntPair( 732, 3 ) );
pairs.push_back( IntIntPair( 732, 4 ) );
for( IntIntPairVectorCI ii = pairs.begin(); ii != pairs.end(); ii++ ){
if( ii->first < positions.size() && ii->second < positions.size() ){
double d = positionGenerator.getDistance( ii->first, ii->second, positions );
(void) fprintf( log, "Distance %6d %6d %15.7e d2=%15.7e\n", ii->first, ii->second, d, d*d );
}
}
}
const int numberOfParameters = 5;
const int ChargeIndex = 0;
const int SigmaIndex = 1;
const int EpsIndex = 2;
const int GammaIndex = 3;
const int LambdaIndex = 4;
std::vector<double> parameterLowerBound( numberOfParameters, 0.0 );
double fixedCharge = 1.0;
parameterLowerBound[ChargeIndex] = fixedCharge; // charge
parameterLowerBound[SigmaIndex] = 0.1; // sigma
parameterLowerBound[EpsIndex] = 0.5; // eps
parameterLowerBound[GammaIndex] = 0.1; // gamma
parameterLowerBound[LambdaIndex] = lambda1; // lambda
std::vector<double> parameterUpperBound( parameterLowerBound );
parameterUpperBound[ChargeIndex] = fixedCharge; // charge
parameterUpperBound[SigmaIndex] = 0.3; // sigma
parameterUpperBound[EpsIndex] = 40.0; // eps
parameterUpperBound[GammaIndex] = 40.0; // gamma
#if IMPLICIT_SOLVENT == OBC_FLAG
parameterLowerBound[GammaIndex] = 0.1; // overlap factor
parameterUpperBound[GammaIndex] = 1.5;
#endif
std::vector<double> parameters( numberOfParameters );
double charge = fixedCharge;
for( int ii = 0; ii < numMolecules; ii++) {
charge *= -1.0;
double lambda = ii < (numMolecules/2) ? lambda1 : lambda2;
randomizeParameters( parameterLowerBound, parameterUpperBound, sfmt, parameters );
#ifdef USE_SOFTCORE
nonbondedSoftcoreForce->addParticle( charge, parameters[SigmaIndex], parameters[EpsIndex], lambda );
gbviSoftcoreForce->addParticle( charge, parameters[SigmaIndex], parameters[GammaIndex], lambda );
#else
nonbondedSoftcoreForce->addParticle( charge, parameters[SigmaIndex], parameters[EpsIndex] );
gbviSoftcoreForce->addParticle( charge, parameters[SigmaIndex], parameters[GammaIndex] );
#endif
int baseParticleIndex = ii*numParticlesPerMolecule;
for( int jj = 1; jj < numParticlesPerMolecule; jj++) {
// alternate charges
charge *= -1.0;
randomizeParameters( parameterLowerBound, parameterUpperBound, sfmt, parameters );
#ifdef USE_SOFTCORE
nonbondedSoftcoreForce->addParticle( charge, parameters[SigmaIndex], parameters[EpsIndex], lambda );
gbviSoftcoreForce->addParticle( charge, parameters[SigmaIndex], parameters[GammaIndex], lambda );
#else
nonbondedSoftcoreForce->addParticle( charge, parameters[SigmaIndex], parameters[EpsIndex] );
gbviSoftcoreForce->addParticle( charge, parameters[SigmaIndex], parameters[GammaIndex] );
#endif
nonbondedSoftcoreForce->addException( baseParticleIndex, baseParticleIndex+jj, 0.0f, 1.0, 0.0f );
#if IMPLICIT_SOLVENT == GBVI_FLAG
double bondDistance = positionGenerator.getDistance( baseParticleIndex, baseParticleIndex+jj, positions );
gbviSoftcoreForce->addBond( baseParticleIndex, baseParticleIndex+jj, bondDistance );
#endif
}
// alternate charge if numParticlesPerMolecule is odd
if( (numParticlesPerMolecule % 2) ){
charge *= -1.0;
}
}
standardSystem.addForce(nonbondedSoftcoreForce);
if( includeGbvi ){
standardSystem.addForce(gbviSoftcoreForce);
}
// copy system and forces
System* systemCopy = copySystem( standardSystem );
#ifdef USE_SOFTCORE
NonbondedSoftcoreForce* nonbondedSoftcoreForceCopy;
nonbondedSoftcoreForceCopy = copyNonbondedSoftcoreForce( *nonbondedSoftcoreForce );
#else
NonbondedForce* nonbondedSoftcoreForceCopy;
nonbondedSoftcoreForceCopy = copyNonbondedForce( *nonbondedSoftcoreForce );
#endif
systemCopy->addForce( nonbondedSoftcoreForceCopy );
std::stringstream baseFileName;
if( includeGbvi ){
#ifdef USE_SOFTCORE
#if IMPLICIT_SOLVENT == GBVI_FLAG
GBVISoftcoreForce* gBVISoftcoreForceCopy = copyGbviSoftcoreForce( *gbviSoftcoreForce );
baseFileName << "GBVISoftcore";
#endif
#if IMPLICIT_SOLVENT == OBC_FLAG
baseFileName << "GBSAObcSoftcore";
GBSAOBCSoftcoreForce* gBVISoftcoreForceCopy = copyGBSAOBCSoftcoreForce( *gbviSoftcoreForce );
#endif
baseFileName << "_lbda" << std::fixed << setprecision(2) << lambda2;
#else
#if IMPLICIT_SOLVENT == GBVI_FLAG
GBVIForce* gBVISoftcoreForceCopy = copyGbviForce( *gbviSoftcoreForce );
baseFileName << "Gbvi";
#endif
#if IMPLICIT_SOLVENT == OBC_FLAG
GBSAOBCForce* gBVISoftcoreForceCopy = copyGbsaObcForce( *gbviSoftcoreForce );
baseFileName << "GBSAOBC";
#endif
#endif
systemCopy->addForce( gBVISoftcoreForceCopy );
}
// perform comparison
std::stringstream idString;
idString << "Nb " << nonbondedMethod << " l2 " << std::fixed << setprecision(2) << lambda2;
runSystemComparisonTest( standardSystem, *systemCopy, "Cuda", "Reference", positions, inputArgumentMap, idString.str(), log );
// serialize
baseFileName << "_N" << positions.size();
baseFileName << "_Nb" << nonbondedMethod;
serializeSystemAndPositions( standardSystem, positions, baseFileName.str(), log);
delete systemCopy;
}
int main() {
try {
registerFreeEnergyCudaKernelFactories( );
VectorOfMapStringToDouble vectorOfMapStringToDouble;
MapStringToDouble inputArgumentMap;
MapStringToDoubleVector generativeArgumentMaps;
//FILE* log = stderr;
FILE* log = NULL;
/*
testSingleParticle( log );
testEnergyEthaneSwitchingFunction( 0, log );
testEnergyEthaneSwitchingFunction( 1, log );
*/
inputArgumentMap["lambda2"] = 1.0;
inputArgumentMap["nonbondedMethod"] = 0;
inputArgumentMap["numMolecules"] = 10;
inputArgumentMap["boxSize"] = 5.0;
inputArgumentMap["positionPlacementMethod"] = 0;
inputArgumentMap["cutoffDistance"] = 0.3*inputArgumentMap["boxSize"];
//inputArgumentMap["cutoffDistance"] = 1.0;
inputArgumentMap["relativeTolerance"] = 5.0e-04;
inputArgumentMap["serialize"] = 1;
//inputArgumentMap["numParticlesPerMolecule"] = 2;
#ifdef USE_SOFTCORE
DoubleVector lamda2;
lamda2.push_back( 1.0 );
//lamda2.push_back( 0.5 );
//lamda2.push_back( 0.0 );
if( lamda2.size() > 0 ){
generativeArgumentMaps["lambda2"] = lamda2;
inputArgumentMap["lambda2"] = lamda2[0];
}
#endif
DoubleVector numberOfMolecules;
numberOfMolecules.push_back( 10 );
numberOfMolecules.push_back( 100 );
numberOfMolecules.push_back( 1000 );
//numberOfMolecules.push_back( 2000 );
//numberOfMolecules.push_back( 4000 );
//numberOfMolecules.push_back( 8000 );
if( numberOfMolecules.size() > 0 ){
generativeArgumentMaps["numMolecules"] = numberOfMolecules;
inputArgumentMap["numMolecules"] = numberOfMolecules[0];
}
DoubleVector nonbondedMethod;
nonbondedMethod.push_back( 0 );
nonbondedMethod.push_back( 1 );
nonbondedMethod.push_back( 2 );
if( nonbondedMethod.size() > 0 ){
generativeArgumentMaps["nonbondedMethod"] = nonbondedMethod;
inputArgumentMap["nonbondedMethod"] = nonbondedMethod[0];
}
vectorOfMapStringToDouble.push_back( inputArgumentMap );
generateInputArgumentMapsFromStringVectors( generativeArgumentMaps, vectorOfMapStringToDouble );
// big box/many particle tests
//bool bigBox = true;
bool bigBox = false;
if( bigBox ){
MapStringToDouble inputArgumentMapBig;
VectorOfMapStringToDouble vectorOfMapStringToDoubleBig;
inputArgumentMapBig["lambda2"] = 1.0;
inputArgumentMapBig["nonbondedMethod"] = 1;
inputArgumentMapBig["numMolecules"] = 10;
inputArgumentMapBig["boxSize"] = 20.0;
inputArgumentMapBig["relativeTolerance"] = 6.0e-04;
vectorOfMapStringToDoubleBig.push_back( inputArgumentMapBig );
//MapStringToDoubleVector generativeArgumentMapsBig;
numberOfMolecules.resize( 0 );
numberOfMolecules.push_back( 4000 );
generativeArgumentMaps["numMolecules"] = numberOfMolecules;
nonbondedMethod.resize( 0 );
nonbondedMethod.push_back( 1 );
nonbondedMethod.push_back( 2 );
generativeArgumentMaps["nonbondedMethod"] = nonbondedMethod;
generateInputArgumentMapsFromStringVectors( generativeArgumentMaps, vectorOfMapStringToDoubleBig );
vectorOfMapStringToDouble.resize( 0 );
vectorOfMapStringToDouble.insert( vectorOfMapStringToDouble.end(), vectorOfMapStringToDoubleBig.begin(), vectorOfMapStringToDoubleBig.end() );
}
if( log ){
MapStringToInt exclude;
exclude["lambda1"] = 1;
exclude["numParticlesPerMolecule"] = 1;
std::stringstream outputStream;
std::sort( vectorOfMapStringToDouble.begin(), vectorOfMapStringToDouble.end(), TestMapSortPredicate);
StringVector printOrder;
printOrder.push_back( "numMolecules" );
printOrder.push_back( "nonbondedMethod" );
printOrder.push_back( "lambda2" );
printOrder.push_back( "boxSize" );
for( unsigned int kk = 0; kk < vectorOfMapStringToDouble.size(); kk++ ){
streamArgumentMapOneLine( vectorOfMapStringToDouble[kk], exclude, printOrder, kk, outputStream );
}
(void) fprintf( log, "Initial argument maps: %u\n%s", static_cast<unsigned int>(vectorOfMapStringToDouble.size()), outputStream.str().c_str() );
}
// run tests
for( unsigned int kk = 0; kk < vectorOfMapStringToDouble.size(); kk++ ){
testGbviSoftcore( vectorOfMapStringToDouble[kk], log );
sleep(2);
}
} catch(const exception& e) {
cout << "exception: " << e.what() << endl;
......
plugins/freeEnergy/platforms/cuda/tests/TestCudaOBCSoftcoreForce.cpp
View file @ 032e18de
......@@ -30,42 +30,52 @@
* -------------------------------------------------------------------------- */
/**
* This tests the reference implementation of GBVIForce.
* This tests the Cuda implementations of GBSAOBCSoftcoreForce
*/
#include "TestCudaSoftcoreForce.h"
#include "openmm/System.h"
#include "../../../tests/AssertionUtilities.h"
//#define USE_SOFTCORE
//#define IMPLICIT_SOLVENT GBVI
//#define IMPLICIT_SOLVENT OBC
#include "sfmt/SFMT.h"
#include "openmm/Context.h"
#define OBC_FLAG 1
#define GBVI_FLAG 2
#include "openmm/GBVIForce.h"
#include "openmm/GBSAOBCForce.h"
#include "openmm/NonbondedForce.h"
#ifdef USE_SOFTCORE
#include "openmm/GBVISoftcoreForce.h"
#include "openmm/GBSAOBCSoftcoreForce.h"
#include "openmm/NonbondedSoftcoreForce.h"
#endif
#include "openmm/GBSAOBCSoftcoreForce.h"
#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
#include "openmm/VerletIntegrator.h"
#include <iostream>
#include <vector>
#include <algorithm>
#include <iomanip>
typedef std::vector<double> DoubleVector;
typedef DoubleVector::iterator DoubleVectorI;
typedef DoubleVector::const_iterator DoubleVectorCI;
typedef std::vector<DoubleVector> DoubleVectorVector;
typedef std::pair<int, double> IntDoublePair;
typedef std::vector<IntDoublePair> IntDoublePairVector;
typedef IntDoublePairVector::iterator IntDoublePairVectorI;
typedef IntDoublePairVector::const_iterator IntDoublePairVectorCI;
extern "C" void registerFreeEnergyCudaKernelFactories();
using namespace OpenMM;
using namespace std;
void testSingleParticle( FILE* log ) {
System system;
system.addParticle(2.0);
VerletIntegrator integrator(0.01);
GBVISoftcoreForce* forceField = new GBVISoftcoreForce;
GBSAOBCSoftcoreForce* forceField = new GBSAOBCSoftcoreForce();
double charge = 1.0;
double charge = 0.5;
double radius = 0.15;
double gamma = 1.0;
forceField->addParticle(charge, radius, gamma);
double scaleFactor = 1.0;
forceField->addParticle(charge, radius, scaleFactor);
system.addForce(forceField);
NonbondedSoftcoreForce* nonbonded = new NonbondedSoftcoreForce();
......@@ -79,12 +89,13 @@ void testSingleParticle( FILE* log ) {
context.setPositions(positions);
State state = context.getState(State::Energy);
double bornRadius = radius;
double bornRadius = radius-0.009;
double eps0 = EPSILON0;
double tau = (1.0/forceField->getSoluteDielectric()-1.0/forceField->getSolventDielectric());
double bornEnergy = (-charge*charge/(8*PI_M*eps0))*tau/bornRadius;
double nonpolarEnergy = -gamma*tau*std::pow( radius/bornRadius, 3.0);
double bornEnergy = (-0.5*0.5/(8*PI_M*eps0))*(1.0/forceField->getSoluteDielectric()-1.0/forceField->getSolventDielectric())/bornRadius;
double extendedRadius = bornRadius+0.14; // probe radius
double nonpolarEnergy = CAL2JOULE*PI_M*0.0216*(10*extendedRadius)*(10*extendedRadius)*std::pow(0.15/bornRadius, 6.0); // Where did this formula come from? Just copied it from CpuImplicitSolvent.cpp
double expectedE = (bornEnergy+nonpolarEnergy);
double obtainedE = state.getPotentialEnergy();
......@@ -96,6 +107,142 @@ void testSingleParticle( FILE* log ) {
ASSERT_EQUAL_TOL((bornEnergy+nonpolarEnergy), state.getPotentialEnergy(), 0.01);
}
/**
* Predicate for sorting <int,double> pair
*
* @param d1 first IntDoublePair to compare
* @param d2 second IntDoublePair to compare
*
*/
bool TestIntDoublePair( const IntDoublePair& d1, const IntDoublePair& d2 ){
return d1.second < d2.second;
}
/**---------------------------------------------------------------------------------------
*
* Get relative difference between two forces
*
* @param f1 force1
* @param f2 force2
* @param forceNorm1 output norm of force1
* @param forceNorm2 output norm of force2
* @param relativeDiff output relative difference between force norms
* @param log if set, output forces
*
*
--------------------------------------------------------------------------------------- */
static void getForceRelativeDifference( const Vec3& f1, const Vec3& f2, double& forceNorm1, double& forceNorm2,
double& relativeDiff, FILE* log ) {
double diff = (f1[0] - f2[0])*(f1[0] - f2[0]) +
(f1[1] - f2[1])*(f1[1] - f2[1]) +
(f1[2] - f2[2])*(f1[2] - f2[2]);
forceNorm1 = sqrt( f1[0]*f1[0] + f1[1]*f1[1] + f1[2]*f1[2] );
forceNorm2 = sqrt( f2[0]*f2[0] + f2[1]*f2[1] + f2[2]*f2[2] );
if( forceNorm1 > 0.0 || forceNorm2 > 0.0 ){
relativeDiff = 2.0*sqrt( diff )/(forceNorm1+forceNorm2);
} else {
relativeDiff = 0.0;
}
return;
}
/**---------------------------------------------------------------------------------------
*
* Compare forces from two states
*
* @param state1 state1
* @param state2 state2
* @param relativeTolerance relative tolerance
* @param log if set, output forces
*
* @return number of entries with relative difference > tolerance
*
--------------------------------------------------------------------------------------- */
int compareForcesOfTwoStates( State& state1, State& state2, double relativeTolerance,
DoubleVector& stats, FILE* log ) {
int error = 0;
vector<Vec3> force1 = state1.getForces();
vector<Vec3> force2 = state2.getForces();
double averageRelativeDifference = 0.0;
double count = 0.0;
DoubleVector medians1( force1.size() );
DoubleVector medians2( force1.size() );
IntDoublePairVector relativeDifferences;
for( unsigned int ii = 0; ii < force1.size(); ii++ ){
double forceNorm1;
double forceNorm2;
double relativeDiff;
getForceRelativeDifference( force1[ii], force2[ii], forceNorm1, forceNorm2, relativeDiff, log );
medians1[ii] = forceNorm1;
medians2[ii] = forceNorm2;
relativeDifferences.push_back( IntDoublePair(ii, relativeDiff ) );
averageRelativeDifference += relativeDiff;
count += 1.0;
if( relativeDiff > relativeTolerance ){
error++;
}
if( log ){
(void) fprintf( log, "F %6u %15.7e [%15.7e %15.7e %15.7e] [%15.7e %15.7e %15.7e] %15.7e %15.7e %s\n", static_cast<unsigned int>(ii),
relativeDiff, force1[ii][0], force1[ii][1], force1[ii][2], force2[ii][0], force2[ii][1], force2[ii][2],
forceNorm1, forceNorm2, (relativeDiff < relativeTolerance ? "":"XXXXXX") );
}
}
// sort relative differences
std::sort( relativeDifferences.begin(), relativeDifferences.end(), TestIntDoublePair );
if( log ){
(void) fprintf( log, "\nEntries w/ largest relative differences.\n" );
for( unsigned int ii = relativeDifferences.size()-1; ii >= relativeDifferences.size()-10 && ii >= 0; ii-- ){
double forceNorm1;
double forceNorm2;
double relativeDiff;
int index = relativeDifferences[ii].first;
getForceRelativeDifference( force1[index], force2[index], forceNorm1, forceNorm2, relativeDiff, log );
(void) fprintf( log, "Fs %6u %15.7e [%15.7e %15.7e %15.7e] [%15.7e %15.7e %15.7e] %15.7e %15.7e %s\n",
static_cast<unsigned int>(index), relativeDiff,
force1[index][0], force1[index][1], force1[index][2],
force2[index][0], force2[index][1], force2[index][2],
forceNorm1, forceNorm2, (relativeDiff < relativeTolerance ? "":"XXXXXX") );
}
}
if( count > 0.0 ){
averageRelativeDifference /= count;
}
std::sort( medians1.begin(), medians1.end() );
std::sort( medians2.begin(), medians2.end() );
double median1 = medians1[medians1.size()/2];
double median2 = medians2[medians2.size()/2];
stats.resize( 4 );
stats[0] = averageRelativeDifference;
IntDoublePair pair = relativeDifferences[relativeDifferences.size()-1];
stats[1] = pair.second;
stats[2] = static_cast<double>(pair.first);
stats[3] = median1 < median2 ? median1 : median2;
return error;
}
void testEnergyEthaneSwitchingFunction( int useSwitchingFunction, FILE* log ) {
std::string methodName = "testEnergyEthaneSwitchingFunction";
......@@ -150,7 +297,7 @@ void testEnergyEthaneSwitchingFunction( int useSwitchingFunction, FILE* log ) {
bornRadiusScaleFactorsEven, bornRadiusScaleFactorsOdd);
}
GBVISoftcoreForce* forceField = new GBVISoftcoreForce();
GBSAOBCSoftcoreForce* forceField = new GBSAOBCSoftcoreForce();
for( int i = 0; i < numParticles; i++ ){
forceField->addParticle( H_charge, H_radius, H_gamma, (i%2) ? bornRadiusScaleFactorsOdd : bornRadiusScaleFactorsEven);
nonbonded->addParticle( H_charge, H_radius, 0.0);
......@@ -165,19 +312,6 @@ void testEnergyEthaneSwitchingFunction( int useSwitchingFunction, FILE* log ) {
// forceField->setParticleParameters( 8, C_charge, (C_radius+0.5), C_gamma, bornRadiusScaleFactorsEven);
// nonbonded->setParticleParameters( 8, C_charge, C_radius, 0.0);
if( useSwitchingFunction ){
forceField->setBornRadiusScalingMethod( GBVISoftcoreForce::QuinticSpline );
} else {
forceField->setBornRadiusScalingMethod( GBVISoftcoreForce::NoScaling );
}
forceField->addBond( 0, 1, C_HBondDistance );
forceField->addBond( 2, 1, C_HBondDistance );
forceField->addBond( 3, 1, C_HBondDistance );
forceField->addBond( 1, 4, C_CBondDistance );
forceField->addBond( 5, 4, C_HBondDistance );
forceField->addBond( 6, 4, C_HBondDistance );
forceField->addBond( 7, 4, C_HBondDistance );
std::vector<pair<int, int> > bonds;
std::vector<double> bondDistances;
......@@ -306,465 +440,7 @@ void testEnergyEthaneSwitchingFunction( int useSwitchingFunction, FILE* log ) {
(void) fprintf( log, "r48=%14.6e r28=%14.6e r24=%14.6e\n", positions[8][2]-positions[4][2], positions[8][2], positions[4][2] );
}
}
#if 0
int carbonIndex = 1;
int hydrogenIndex = 0;
while( hydrogenIndex < 8 ){
Vec3 carbonDelta;
for( int kk = 0; kk < 3; kk++ ){
positions[hydrogenIndex][kk] += positionIncrement*(positions[carbonIndex][kk] - positions[hydrogenIndex][kk] );
}
double dist = 0.0;
for( int kk = 0; kk < 3; kk++ ){
dist += (positions[carbonIndex][kk] - positions[hydrogenIndex][kk] )*(positions[carbonIndex][kk] - positions[hydrogenIndex][kk]);
}
(void) fprintf( log, "H=%d C=%d r=%14.6e\n", hydrogenIndex, carbonIndex, dist );
hydrogenIndex++;
if( hydrogenIndex == carbonIndex ){
hydrogenIndex++;
}
if( carbonIndex == 1 && hydrogenIndex == 4 ){
carbonIndex = 4;
hydrogenIndex = 5;
}
}
#endif
}
}
static GBVISoftcoreForce* copyGbviSoftcoreForce( const GBVISoftcoreForce& gbviSoftcoreForce ){
GBVISoftcoreForce* copyGbviSoftcoreForce = new GBVISoftcoreForce(gbviSoftcoreForce);
/*
GBVISoftcoreForce* copyGbviSoftcoreForce = new GBVISoftcoreForce();
copyGbviSoftcoreForce->setNonbondedMethod( gbviSoftcoreForce.getNonbondedMethod() );
copyGbviSoftcoreForce->setCutoffDistance( gbviSoftcoreForce.getCutoffDistance() );
copyGbviSoftcoreForce->setSolventDielectric( gbviSoftcoreForce.getSolventDielectric() );
copyGbviSoftcoreForce->setSoluteDielectric( gbviSoftcoreForce.getSoluteDielectric() );
copyGbviSoftcoreForce->setBornRadiusScalingMethod( gbviSoftcoreForce.getBornRadiusScalingMethod() );
copyGbviSoftcoreForce->setQuinticLowerLimitFactor( gbviSoftcoreForce.getQuinticLowerLimitFactor() );
copyGbviSoftcoreForce->setQuinticUpperBornRadiusLimit( gbviSoftcoreForce.getQuinticUpperBornRadiusLimit() );
// particle parameters
for( unsigned int ii = 0; ii < gbviSoftcoreForce.getNumParticles(); ii++ ){
double charge;
double sigma;
double gamma;
double softcoreLJLambda;
gbviSoftcoreForce.getParticleParameters(ii, charge, sigma, gamma, softcoreLJLambda);
copyGbviSoftcoreForce->addParticle( charge, sigma, gamma, softcoreLJLambda);
}
// bonds
for( unsigned int ii = 0; ii < gbviSoftcoreForce.getNumBonds(); ii++ ){
int particle1, particle2;
double distance;
gbviSoftcoreForce.getBondParameters( ii, particle1, particle2, distance);
copyGbviSoftcoreForce->addBond( particle1, particle2, distance );
}
*/
return copyGbviSoftcoreForce;
}
static GBVIForce* copyGbviForce( const GBVIForce& gbviForce ){
return new GBVIForce(gbviForce);
}
static GBSAOBCSoftcoreForce* copyGBSAOBCSoftcoreForce( const GBSAOBCSoftcoreForce& gbviSoftcoreForce ){
return new GBSAOBCSoftcoreForce(gbviSoftcoreForce);
}
static GBSAOBCForce* copyGbsaObcForce( const GBSAOBCForce& gbviForce ){
return new GBSAOBCForce(gbviForce);
}
void testGbviSoftcore( MapStringToDouble& inputArgumentMap, FILE* log ){
double lambda1 = 1.0;
double lambda2 = 1.0;
int nonbondedMethod = 0;
int numMolecules = 1;
int numParticlesPerMolecule = 2;
int useQuinticSpline = 1;
int applyAssert = 1;
int positionPlacementMethod = 0;
int serialize = 0;
double boxSize = 10.0;
double relativeTolerance = 1.0e-04;
setDoubleFromMapStringToDouble( inputArgumentMap, "lambda1", lambda1 );
setDoubleFromMapStringToDouble( inputArgumentMap, "lambda2", lambda2 );
setDoubleFromMapStringToDouble( inputArgumentMap, "boxSize", boxSize );
double cutoffDistance = boxSize*0.4;;
setDoubleFromMapStringToDouble( inputArgumentMap, "cutoffDistance", cutoffDistance);
setDoubleFromMapStringToDouble( inputArgumentMap, "relativeTolerance", relativeTolerance );
setIntFromMapStringToDouble( inputArgumentMap, "positionPlacementMethod", positionPlacementMethod ) ;
setIntFromMapStringToDouble( inputArgumentMap, "nonbondedMethod", nonbondedMethod );
setIntFromMapStringToDouble( inputArgumentMap, "numMolecules", numMolecules );
setIntFromMapStringToDouble( inputArgumentMap, "numParticlesPerMolecule", numParticlesPerMolecule );
setIntFromMapStringToDouble( inputArgumentMap, "serialize", serialize );
if( nonbondedMethod == 2 && cutoffDistance > boxSize*0.5 ){
cutoffDistance = boxSize*0.5;
}
int numParticles = numMolecules*numParticlesPerMolecule;
int includeGbvi = 1;
double reactionFieldDielectric = 80.0;
if( log ){
double particleDensity = static_cast<double>(numParticles)/(boxSize*boxSize*boxSize);
double particleCube = pow( particleDensity, (-1.0/3.0) );
(void) fprintf( log, "\n--------------------------------------------------------------------------------------\n" );
(void) fprintf( log, "Input arguments\n" );
(void) fflush( log );
(void) fprintf( log, " includeGbvi %d\n", includeGbvi );
(void) fprintf( log, " nonbondedMethod %d\n", nonbondedMethod );
(void) fprintf( log, " numParticles %d\n", numParticles );
(void) fprintf( log, " numMolecules %d\n", numMolecules );
(void) fprintf( log, " numParticlesPerMolecule %d\n", numParticlesPerMolecule );
(void) fprintf( log, " useQuinticSpline %d\n", useQuinticSpline );
(void) fprintf( log, " positionPlacementMethod %d\n", positionPlacementMethod);
#ifdef USE_SOFTCORE
(void) fprintf( log, " lambda1 %8.3f\n", lambda1 );
(void) fprintf( log, " lambda2 %8.3f\n", lambda2 );
#endif
(void) fprintf( log, " boxSize %8.3f\n", boxSize );
(void) fprintf( log, " cutoffDistance %8.3f\n", cutoffDistance );
(void) fprintf( log, " reactionFieldDielectric %8.3f\n", reactionFieldDielectric );
(void) fprintf( log, " relativeTolerance %8.1e\n", relativeTolerance );
(void) fprintf( log, " particleDensity %8.2e\n", particleDensity );
(void) fprintf( log, " particleCube %8.2e\n", particleCube );
}
// Create two systems: one with GbviSoftcoreForce NonbondedSoftcoreForce forces, and one using a CustomNonbondedForce, CustomGBVI force to implement the same interaction.
System standardSystem;
for (int i = 0; i < numParticles; i++) {
standardSystem.addParticle(1.0);
}
standardSystem.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
#ifdef USE_SOFTCORE
NonbondedSoftcoreForce* nonbondedSoftcoreForce = new NonbondedSoftcoreForce();
if( nonbondedMethod == NoCutoff ){
nonbondedSoftcoreForce->setNonbondedMethod( NonbondedSoftcoreForce::NoCutoff );
} else {
if( nonbondedMethod == CutoffNonPeriodic ){
nonbondedSoftcoreForce->setNonbondedMethod( NonbondedSoftcoreForce::CutoffNonPeriodic );
} else {
nonbondedSoftcoreForce->setNonbondedMethod( NonbondedSoftcoreForce::CutoffPeriodic );
}
}
#else
NonbondedForce* nonbondedSoftcoreForce = new NonbondedForce();
if( nonbondedMethod == NoCutoff ){
nonbondedSoftcoreForce->setNonbondedMethod( NonbondedForce::NoCutoff );
} else {
if( nonbondedMethod == CutoffNonPeriodic ){
nonbondedSoftcoreForce->setNonbondedMethod( NonbondedForce::CutoffNonPeriodic );
} else {
nonbondedSoftcoreForce->setNonbondedMethod( NonbondedForce::CutoffPeriodic );
}
}
#endif
nonbondedSoftcoreForce->setCutoffDistance( cutoffDistance );
nonbondedSoftcoreForce->setReactionFieldDielectric( reactionFieldDielectric );
#ifdef USE_SOFTCORE
#if IMPLICIT_SOLVENT == GBVI_FLAG
GBVISoftcoreForce* gbviSoftcoreForce = new GBVISoftcoreForce();
if( nonbondedMethod == NoCutoff ){
gbviSoftcoreForce->setNonbondedMethod( GBVISoftcoreForce::NoCutoff );
} else {
if( nonbondedMethod == CutoffNonPeriodic ){
gbviSoftcoreForce->setNonbondedMethod( GBVISoftcoreForce::CutoffNonPeriodic );
} else {
gbviSoftcoreForce->setNonbondedMethod( GBVISoftcoreForce::CutoffPeriodic );
}
}
#else
GBSAOBCSoftcoreForce* gbviSoftcoreForce = new GBSAOBCSoftcoreForce();
if( nonbondedMethod == NoCutoff ){
gbviSoftcoreForce->setNonbondedMethod( GBSAOBCSoftcoreForce::NoCutoff );
} else {
if( nonbondedMethod == CutoffNonPeriodic ){
gbviSoftcoreForce->setNonbondedMethod( GBSAOBCSoftcoreForce::CutoffNonPeriodic );
} else {
gbviSoftcoreForce->setNonbondedMethod( GBSAOBCSoftcoreForce::CutoffPeriodic );
}
}
#endif
#else
#if IMPLICIT_SOLVENT == GBVI_FLAG
GBVIForce* gbviSoftcoreForce = new GBVIForce();
if( nonbondedMethod == NoCutoff ){
gbviSoftcoreForce->setNonbondedMethod( GBVIForce::NoCutoff );
} else {
if( nonbondedMethod == CutoffNonPeriodic ){
gbviSoftcoreForce->setNonbondedMethod( GBVIForce::CutoffNonPeriodic );
} else {
gbviSoftcoreForce->setNonbondedMethod( GBVIForce::CutoffPeriodic );
}
}
#else
GBSAOBCForce* gbviSoftcoreForce = new GBSAOBCForce();
if( nonbondedMethod == NoCutoff ){
gbviSoftcoreForce->setNonbondedMethod( GBSAOBCForce::NoCutoff );
} else {
if( nonbondedMethod == CutoffNonPeriodic ){
gbviSoftcoreForce->setNonbondedMethod( GBSAOBCForce::CutoffNonPeriodic );
} else {
gbviSoftcoreForce->setNonbondedMethod( GBSAOBCForce::CutoffPeriodic );
}
}
#endif
#endif
#if IMPLICIT_SOLVENT == GBVI_FLAG
#ifdef USE_SOFTCORE
if( useQuinticSpline ){
gbviSoftcoreForce->setBornRadiusScalingMethod( GBVISoftcoreForce::QuinticSpline );
} else {
gbviSoftcoreForce->setBornRadiusScalingMethod( GBVISoftcoreForce::NoScaling );
}
#else
if( useQuinticSpline ){
gbviSoftcoreForce->setBornRadiusScalingMethod( GBVIForce::QuinticSpline );
} else {
gbviSoftcoreForce->setBornRadiusScalingMethod( GBVIForce::NoScaling );
}
#endif
#endif
gbviSoftcoreForce->setSolventDielectric( 78.3 );
//gbviSoftcoreForce->setSolventDielectric( 1.0e+10 );
//gbviSoftcoreForce->setSolventDielectric( 1.0 );
gbviSoftcoreForce->setSoluteDielectric( 1.0 );
gbviSoftcoreForce->setCutoffDistance( nonbondedSoftcoreForce->getCutoffDistance( ) );
std::vector<Vec3> positions(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
PositionGenerator positionGenerator( numMolecules, numParticlesPerMolecule, boxSize );
if( log ){
positionGenerator.setLog( log );
}
if( positionPlacementMethod == 1 ){
positionGenerator.setPositions( PositionGenerator::SimpleGrid, sfmt, positions );
} else {
positionGenerator.setBondDistance( 0.3 );
positionGenerator.setPositions( PositionGenerator::Random, sfmt, positions );
}
// show info on particle positions
if( log ){
Vec3 box[2];
positionGenerator.getEnclosingBox( positions, box );
(void) fprintf( log, "Enclosing Box (in A): [%15.7e %15.7e] [%15.7e %15.7e] [%15.7e %15.7e] [%15.7e %15.7e %15.7e]\n",
box[0][0], box[1][0], box[0][1], box[1][1], box[0][2], box[1][2],
(box[1][0] - box[0][0]), (box[1][1] - box[0][1]), (box[1][2] - box[0][2]) );
int showIndex = 5;
int periodicBoundaryConditions = (nonbondedMethod == 2) ? 1 : 0;
IntVector positionIndexVector;
positionIndexVector.push_back( 0 );
positionIndexVector.push_back( static_cast<int>(positions.size())-1 );
//positionIndexVector.push_back( 542 );
for( unsigned int ii = 0; ii < positionIndexVector.size(); ii++ ){
if( positionIndexVector[ii] < positions.size() ){
int positionIndex = positionIndexVector[ii];
IntDoublePairVector sortVector;
positionGenerator.getSortedDistances( periodicBoundaryConditions, positionIndex, positions, sortVector );
(void) fprintf( log, "Min/max distance from %6d:\n ", positionIndex );
for( unsigned int jj = 0; jj < sortVector.size() && jj < showIndex; jj++ ){
IntDoublePair pair = sortVector[jj];
(void) fprintf( log, "[%6d %15.7e] ", pair.first, pair.second);
}
(void) fprintf( log, "\n " );
for( unsigned int jj = (sortVector.size() - showIndex); jj < sortVector.size() && jj >= 0; jj++ ){
IntDoublePair pair = sortVector[jj];
(void) fprintf( log, "[%6d %15.7e] ", pair.first, pair.second);
}
(void) fprintf( log, "\n" );
}
}
IntIntPairVector pairs;
pairs.push_back( IntIntPair( 732, 0 ) );
pairs.push_back( IntIntPair( 732, 1 ) );
pairs.push_back( IntIntPair( 732, 2 ) );
pairs.push_back( IntIntPair( 732, 3 ) );
pairs.push_back( IntIntPair( 732, 4 ) );
for( IntIntPairVectorCI ii = pairs.begin(); ii != pairs.end(); ii++ ){
if( ii->first < positions.size() && ii->second < positions.size() ){
double d = positionGenerator.getDistance( ii->first, ii->second, positions );
(void) fprintf( log, "Distance %6d %6d %15.7e d2=%15.7e\n", ii->first, ii->second, d, d*d );
}
}
}
const int numberOfParameters = 5;
const int ChargeIndex = 0;
const int SigmaIndex = 1;
const int EpsIndex = 2;
const int GammaIndex = 3;
const int LambdaIndex = 4;
std::vector<double> parameterLowerBound( numberOfParameters, 0.0 );
double fixedCharge = 1.0;
parameterLowerBound[ChargeIndex] = fixedCharge; // charge
parameterLowerBound[SigmaIndex] = 0.1; // sigma
parameterLowerBound[EpsIndex] = 0.5; // eps
parameterLowerBound[GammaIndex] = 0.1; // gamma
parameterLowerBound[LambdaIndex] = lambda1; // lambda
std::vector<double> parameterUpperBound( parameterLowerBound );
parameterUpperBound[ChargeIndex] = fixedCharge; // charge
parameterUpperBound[SigmaIndex] = 0.3; // sigma
parameterUpperBound[EpsIndex] = 40.0; // eps
parameterUpperBound[GammaIndex] = 40.0; // gamma
#if IMPLICIT_SOLVENT == OBC_FLAG
parameterLowerBound[GammaIndex] = 0.1; // overlap factor
parameterUpperBound[GammaIndex] = 1.5;
#endif
std::vector<double> parameters( numberOfParameters );
double charge = fixedCharge;
for( int ii = 0; ii < numMolecules; ii++) {
charge *= -1.0;
double lambda = ii < (numMolecules/2) ? lambda1 : lambda2;
randomizeParameters( parameterLowerBound, parameterUpperBound, sfmt, parameters );
#ifdef USE_SOFTCORE
nonbondedSoftcoreForce->addParticle( charge, parameters[SigmaIndex], parameters[EpsIndex], lambda );
gbviSoftcoreForce->addParticle( charge, parameters[SigmaIndex], parameters[GammaIndex], lambda );
#else
nonbondedSoftcoreForce->addParticle( charge, parameters[SigmaIndex], parameters[EpsIndex] );
gbviSoftcoreForce->addParticle( charge, parameters[SigmaIndex], parameters[GammaIndex] );
#endif
int baseParticleIndex = ii*numParticlesPerMolecule;
for( int jj = 1; jj < numParticlesPerMolecule; jj++) {
// alternate charges
charge *= -1.0;
randomizeParameters( parameterLowerBound, parameterUpperBound, sfmt, parameters );
#ifdef USE_SOFTCORE
nonbondedSoftcoreForce->addParticle( charge, parameters[SigmaIndex], parameters[EpsIndex], lambda );
gbviSoftcoreForce->addParticle( charge, parameters[SigmaIndex], parameters[GammaIndex], lambda );
#else
nonbondedSoftcoreForce->addParticle( charge, parameters[SigmaIndex], parameters[EpsIndex] );
gbviSoftcoreForce->addParticle( charge, parameters[SigmaIndex], parameters[GammaIndex] );
#endif
nonbondedSoftcoreForce->addException( baseParticleIndex, baseParticleIndex+jj, 0.0f, 1.0, 0.0f );
#if IMPLICIT_SOLVENT == GBVI_FLAG
double bondDistance = positionGenerator.getDistance( baseParticleIndex, baseParticleIndex+jj, positions );
gbviSoftcoreForce->addBond( baseParticleIndex, baseParticleIndex+jj, bondDistance );
#endif
}
// alternate charge if numParticlesPerMolecule is odd
if( (numParticlesPerMolecule % 2) ){
charge *= -1.0;
}
}
standardSystem.addForce(nonbondedSoftcoreForce);
if( includeGbvi ){
standardSystem.addForce(gbviSoftcoreForce);
}
// copy system and forces
System* systemCopy = copySystem( standardSystem );
#ifdef USE_SOFTCORE
NonbondedSoftcoreForce* nonbondedSoftcoreForceCopy;
nonbondedSoftcoreForceCopy = copyNonbondedSoftcoreForce( *nonbondedSoftcoreForce );
#else
NonbondedForce* nonbondedSoftcoreForceCopy;
nonbondedSoftcoreForceCopy = copyNonbondedForce( *nonbondedSoftcoreForce );
#endif
systemCopy->addForce( nonbondedSoftcoreForceCopy );
std::stringstream baseFileName;
if( includeGbvi ){
#ifdef USE_SOFTCORE
#if IMPLICIT_SOLVENT == GBVI_FLAG
GBVISoftcoreForce* gBVISoftcoreForceCopy = copyGbviSoftcoreForce( *gbviSoftcoreForce );
baseFileName << "GBVISoftcore";
#endif
#if IMPLICIT_SOLVENT == OBC_FLAG
baseFileName << "GBSAObcSoftcore";
GBSAOBCSoftcoreForce* gBVISoftcoreForceCopy = copyGBSAOBCSoftcoreForce( *gbviSoftcoreForce );
#endif
baseFileName << "_lbda" << std::fixed << setprecision(2) << lambda2;
#else
#if IMPLICIT_SOLVENT == GBVI_FLAG
GBVIForce* gBVISoftcoreForceCopy = copyGbviForce( *gbviSoftcoreForce );
baseFileName << "Gbvi";
#endif
#if IMPLICIT_SOLVENT == OBC_FLAG
GBSAOBCForce* gBVISoftcoreForceCopy = copyGbsaObcForce( *gbviSoftcoreForce );
baseFileName << "GBSAOBC";
#endif
#endif
systemCopy->addForce( gBVISoftcoreForceCopy );
}
// perform comparison
std::stringstream idString;
idString << "Nb " << nonbondedMethod << " l2 " << std::fixed << setprecision(2) << lambda2;
runSystemComparisonTest( standardSystem, *systemCopy, "Cuda", "Reference", positions, inputArgumentMap, idString.str(), log );
// serialize
baseFileName << "_N" << positions.size();
baseFileName << "_Nb" << nonbondedMethod;
serializeSystemAndPositions( standardSystem, positions, baseFileName.str(), log);
delete systemCopy;
}
int main() {
......@@ -773,115 +449,13 @@ int main() {
registerFreeEnergyCudaKernelFactories( );
VectorOfMapStringToDouble vectorOfMapStringToDouble;
MapStringToDouble inputArgumentMap;
MapStringToDoubleVector generativeArgumentMaps;
//FILE* log = stderr;
FILE* log = NULL;
/*
testSingleParticle( log );
testEnergyEthaneSwitchingFunction( 0, log );
testEnergyEthaneSwitchingFunction( 1, log );
*/
inputArgumentMap["lambda2"] = 1.0;
inputArgumentMap["nonbondedMethod"] = 0;
inputArgumentMap["numMolecules"] = 10;
inputArgumentMap["boxSize"] = 5.0;
inputArgumentMap["positionPlacementMethod"] = 0;
inputArgumentMap["cutoffDistance"] = 0.3*inputArgumentMap["boxSize"];
//inputArgumentMap["cutoffDistance"] = 1.0;
inputArgumentMap["relativeTolerance"] = 5.0e-04;
inputArgumentMap["serialize"] = 1;
//inputArgumentMap["numParticlesPerMolecule"] = 2;
#ifdef USE_SOFTCORE
DoubleVector lamda2;
lamda2.push_back( 1.0 );
lamda2.push_back( 0.5 );
lamda2.push_back( 0.0 );
if( lamda2.size() > 0 ){
generativeArgumentMaps["lambda2"] = lamda2;
inputArgumentMap["lambda2"] = lamda2[0];
}
#endif
DoubleVector numberOfMolecules;
numberOfMolecules.push_back( 10 );
numberOfMolecules.push_back( 100 );
numberOfMolecules.push_back( 1000 );
//numberOfMolecules.push_back( 2000 );
//numberOfMolecules.push_back( 4000 );
//numberOfMolecules.push_back( 8000 );
if( numberOfMolecules.size() > 0 ){
generativeArgumentMaps["numMolecules"] = numberOfMolecules;
inputArgumentMap["numMolecules"] = numberOfMolecules[0];
}
DoubleVector nonbondedMethod;
nonbondedMethod.push_back( 0 );
nonbondedMethod.push_back( 1 );
nonbondedMethod.push_back( 2 );
if( nonbondedMethod.size() > 0 ){
generativeArgumentMaps["nonbondedMethod"] = nonbondedMethod;
inputArgumentMap["nonbondedMethod"] = nonbondedMethod[0];
}
vectorOfMapStringToDouble.push_back( inputArgumentMap );
generateInputArgumentMapsFromStringVectors( generativeArgumentMaps, vectorOfMapStringToDouble );
// big box/many particle tests
//bool bigBox = true;
bool bigBox = false;
if( bigBox ){
MapStringToDouble inputArgumentMapBig;
VectorOfMapStringToDouble vectorOfMapStringToDoubleBig;
inputArgumentMapBig["lambda2"] = 1.0;
inputArgumentMapBig["nonbondedMethod"] = 1;
inputArgumentMapBig["numMolecules"] = 10;
inputArgumentMapBig["boxSize"] = 20.0;
inputArgumentMapBig["relativeTolerance"] = 6.0e-04;
vectorOfMapStringToDoubleBig.push_back( inputArgumentMapBig );
//MapStringToDoubleVector generativeArgumentMapsBig;
numberOfMolecules.resize( 0 );
numberOfMolecules.push_back( 4000 );
generativeArgumentMaps["numMolecules"] = numberOfMolecules;
nonbondedMethod.resize( 0 );
nonbondedMethod.push_back( 1 );
nonbondedMethod.push_back( 2 );
generativeArgumentMaps["nonbondedMethod"] = nonbondedMethod;
generateInputArgumentMapsFromStringVectors( generativeArgumentMaps, vectorOfMapStringToDoubleBig );
vectorOfMapStringToDouble.resize( 0 );
vectorOfMapStringToDouble.insert( vectorOfMapStringToDouble.end(), vectorOfMapStringToDoubleBig.begin(), vectorOfMapStringToDoubleBig.end() );
}
if( log ){
MapStringToInt exclude;
exclude["lambda1"] = 1;
exclude["numParticlesPerMolecule"] = 1;
std::stringstream outputStream;
std::sort( vectorOfMapStringToDouble.begin(), vectorOfMapStringToDouble.end(), TestMapSortPredicate);
StringVector printOrder;
printOrder.push_back( "numMolecules" );
printOrder.push_back( "nonbondedMethod" );
printOrder.push_back( "lambda2" );
printOrder.push_back( "boxSize" );
for( unsigned int kk = 0; kk < vectorOfMapStringToDouble.size(); kk++ ){
streamArgumentMapOneLine( vectorOfMapStringToDouble[kk], exclude, printOrder, kk, outputStream );
}
(void) fprintf( log, "Initial argument maps: %u\n%s", static_cast<unsigned int>(vectorOfMapStringToDouble.size()), outputStream.str().c_str() );
}
// run tests
for( unsigned int kk = 0; kk < vectorOfMapStringToDouble.size(); kk++ ){
testGbviSoftcore( vectorOfMapStringToDouble[kk], log );
sleep(2);
}
} catch(const exception& e) {
cout << "exception: " << e.what() << endl;
......
plugins/freeEnergy/platforms/reference/src/gbsa/CpuObcSoftcore.cpp
View file @ 032e18de
......@@ -221,7 +221,7 @@ void CpuObcSoftcore::computeBornRadii( const vector<RealVec>& atomCoordinates,
// OBC-specific code (Eqs. 6-8 in paper)
sum *= nonPolarScaleFactors[atomI]*half*offsetRadiusI;
sum *= half*offsetRadiusI;
RealOpenMM sum2 = sum*sum;
RealOpenMM sum3 = sum*sum2;
RealOpenMM tanhSum = TANH( alphaObc*sum - betaObc*sum2 + gammaObc*sum3 );
......@@ -486,7 +486,7 @@ RealOpenMM CpuObcSoftcore::computeBornEnergyForces( vector<RealVec>& atomCoordin
RealOpenMM r2Inverse = rInverse*rInverse;
RealOpenMM t3 = eighth*(one + scaledRadiusJ2*r2Inverse)*(l_ij2 - u_ij2) + fourth*LN( u_ij/l_ij )*r2Inverse;
t3 *= nonPolarScaleFactors[atomI]*nonPolarScaleFactors[atomJ];
t3 *= nonPolarScaleFactors[atomJ];
RealOpenMM de = bornForces[atomI]*t3*rInverse;
......
plugins/freeEnergy/platforms/reference/tests/TestReferenceGBVISoftcoreForce.cpp
View file @ 032e18de
......@@ -38,6 +38,7 @@
#include "ReferencePlatform.h"
#include "openmm/GBVISoftcoreForce.h"
#include "openmm/GBSAOBCForce.h"
#include "openmm/CustomGBForce.h"
#include "openmm/System.h"
#include "openmm/LangevinIntegrator.h"
#include "openmm/NonbondedForce.h"
......@@ -48,14 +49,19 @@
#include "openmm/freeEnergyKernels.h"
#include "ReferenceFreeEnergyKernelFactory.h"
#include <iostream>
#include <vector>
#define USE_SOFTCORE
#include "TestReferenceSoftcoreForce.h"
#define OBC_FLAG 1
#define GBVI_FLAG 2
#define IMPLICIT_SOLVENT GBVI_FLAG
#include <iomanip>
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
#define PRINT_ON 0
int compareForcesOfTwoStates( int numParticles, State& state1, State& state2, double relativeTolerance, double absoluteTolerance ) {
......@@ -399,13 +405,500 @@ void testEnergyEthaneSwitchingFunction( int useSwitchingFunction ) {
}
}
// computes the scaled radii based on covalent info and atomic radii
#ifdef USE_SOFTCORE
static void findScaledRadii( const GBVISoftcoreForce& gbviForce, std::vector<double> & scaledRadii, FILE* log ) {
#else
static void findScaledRadii( const GBVIForce& gbviForce, std::vector<double> & scaledRadii, FILE* log) {
#endif
int numberOfParticles = gbviForce.getNumParticles();
int numberOfBonds = gbviForce.getNumBonds();
FILE* errorLog = log ? log : stderr;
// load 1-2 atom pairs along w/ bond distance using HarmonicBondForce & constraints
// numberOfBonds < 1, indicating they were not set by the user
if( numberOfBonds < 1 && numberOfParticles > 1 ){
(void) fprintf( errorLog, "Warning: no covalent bonds set for GB/VI force!\n" );
}
std::vector< std::vector<int> > bondIndices( numberOfBonds );
std::vector<double> bondLengths( numberOfBonds );
std::vector<double> radii( numberOfParticles);
scaledRadii.resize(numberOfParticles);
for (int i = 0; i < numberOfParticles; i++) {
double charge, radius, gamma;
#ifdef USE_SOFTCORE
double lambda;
gbviForce.getParticleParameters(i, charge, radius, gamma, lambda);
#else
gbviForce.getParticleParameters(i, charge, radius, gamma);
#endif
radii[i] = radius;
scaledRadii[i] = radius;
}
for (int i = 0; i < numberOfBonds; i++) {
int particle1, particle2;
double bondLength;
gbviForce.getBondParameters(i, particle1, particle2, bondLength);
if (particle1 < 0 || particle1 >= gbviForce.getNumParticles()) {
std::stringstream msg;
msg << "GBVISoftcoreForce: Illegal particle index: ";
msg << particle1;
throw OpenMMException(msg.str());
}
if (particle2 < 0 || particle2 >= gbviForce.getNumParticles()) {
std::stringstream msg;
msg << "GBVISoftcoreForce: Illegal particle index: ";
msg << particle2;
throw OpenMMException(msg.str());
}
if (bondLength < 0 ) {
std::stringstream msg;
msg << "GBVISoftcoreForce: negative bondlength: ";
msg << bondLength;
throw OpenMMException(msg.str());
}
bondIndices[i].push_back( particle1 );
bondIndices[i].push_back( particle2 );
bondLengths[i] = bondLength;
}
// load 1-2 indicies for each atom
std::vector<std::vector<int> > bonded12(numberOfParticles);
for (int i = 0; i < (int) bondIndices.size(); ++i) {
bonded12[bondIndices[i][0]].push_back(i);
bonded12[bondIndices[i][1]].push_back(i);
}
int errors = 0;
// compute scaled radii (Eq. 5 of Labute paper [JCC 29 p. 1693-1698 2008])
for (int j = 0; j < (int) bonded12.size(); ++j){
double radiusJ = radii[j];
double scaledRadiusJ;
if( bonded12[j].size() == 0 ){
if( numberOfParticles > 1 ){
(void) fprintf( errorLog, "Warning GBVIForceImpl::findScaledRadii atom %d has no covalent bonds; using atomic radius=%.3f.\n", j, radiusJ );
}
scaledRadiusJ = radiusJ;
// errors++;
} else {
double rJ2 = radiusJ*radiusJ;
// loop over bonded neighbors of atom j, applying Eq. 5 in Labute
scaledRadiusJ = 0.0;
for (int i = 0; i < (int) bonded12[j].size(); ++i){
int index = bonded12[j][i];
int bondedAtomIndex = (j == bondIndices[index][0]) ? bondIndices[index][1] : bondIndices[index][0];
double radiusI = radii[bondedAtomIndex];
double rI2 = radiusI*radiusI;
double a_ij = (radiusI - bondLengths[index]);
a_ij *= a_ij;
a_ij = (rJ2 - a_ij)/(2.0*bondLengths[index]);
double a_ji = radiusJ - bondLengths[index];
a_ji *= a_ji;
a_ji = (rI2 - a_ji)/(2.0*bondLengths[index]);
scaledRadiusJ += a_ij*a_ij*(3.0*radiusI - a_ij) + a_ji*a_ji*( 3.0*radiusJ - a_ji );
}
scaledRadiusJ = (radiusJ*radiusJ*radiusJ) - 0.125*scaledRadiusJ;
if( scaledRadiusJ > 0.0 ){
scaledRadiusJ = 0.95*pow( scaledRadiusJ, (1.0/3.0) );
} else {
scaledRadiusJ = 0.0;
}
}
if( log ){
//(void) fprintf( stderr, "scaledRadii %d %12.4f\n", j, scaledRadiusJ );
}
scaledRadii[j] = scaledRadiusJ;
}
// abort if errors
if( errors ){
throw OpenMMException("findScaledRadii errors -- aborting");
}
if( log ){
(void) fprintf( log, " R q gamma scaled radii no. bnds\n" );
double totalQ = 0.0;
for( int i = 0; i < (int) scaledRadii.size(); i++ ){
double charge;
double gamma;
double radiusI;
gbviForce.getParticleParameters(i, charge, radiusI, gamma);
totalQ += charge;
(void) fprintf( log, "%4d %14.5e %14.5e %14.5e %14.5e %d\n", i, radiusI, charge, gamma, scaledRadii[i], (int) bonded12[i].size() );
}
(void) fprintf( log, "Total charge=%e\n", totalQ );
(void) fflush( log );
}
return;
}
// load parameters from gbviForce to customGbviForce
// findScaledRadii() is called to calculate the scaled radii (S)
// S is derived quantity in GBVIForce, not a parameter is the case here
#ifdef USE_SOFTCORE
static void loadGbviParameters( const GBVISoftcoreForce& gbviSoftcoreForce, CustomGBForce& customGbviForce, FILE* log ) {
vector<double> params(5);
double lambda;
#else
static void loadGbviParameters( const GBVIForce& gbviSoftcoreForce, CustomGBForce& customGbviForce, FILE* log ) {
vector<double> params(4);
#endif
int numParticles = gbviSoftcoreForce.getNumParticles();
// get scaled radii
std::vector<double> scaledRadii( numParticles );
findScaledRadii( gbviSoftcoreForce, scaledRadii, log);
for( int ii = 0; ii < numParticles; ii++) {
double charge, radius, gamma;
#ifdef USE_SOFTCORE
gbviSoftcoreForce.getParticleParameters(ii, charge, radius, gamma, lambda);
params[4] = lambda;
#else
gbviSoftcoreForce.getParticleParameters(ii, charge, radius, gamma);
#endif
params[0] = charge;
params[1] = radius;
params[2] = scaledRadii[ii];
params[3] = gamma;
customGbviForce.addParticle(params);
}
}
// create custom GB/VI force
#ifdef USE_SOFTCORE
static void createCustomGBVI( CustomGBForce& customGbviForce, const GBVISoftcoreForce& gbviSoftcoreForce, FILE* log ){
#else
static void createCustomGBVI( CustomGBForce& customGbviForce, const GBVIForce& gbviSoftcoreForce, FILE* log ){
#endif
customGbviForce.setCutoffDistance( gbviSoftcoreForce.getCutoffDistance() );
customGbviForce.addPerParticleParameter("q");
customGbviForce.addPerParticleParameter("radius");
customGbviForce.addPerParticleParameter("scaleFactor"); // derived in GBVIForceImpl implmentation, but parameter here
customGbviForce.addPerParticleParameter("gamma");
customGbviForce.addGlobalParameter("solventDielectric", gbviSoftcoreForce.getSolventDielectric() );
customGbviForce.addGlobalParameter("soluteDielectric", gbviSoftcoreForce.getSoluteDielectric() );
#ifdef USE_SOFTCORE
customGbviForce.addPerParticleParameter("lambda");
customGbviForce.addComputedValue("V", " lambda2*(uL - lL + factor3/(radius1*radius1*radius1));"
"uL = 1.5*x2uI*(0.25*rI-0.33333*xuI+0.125*(r2-S2)*rI*x2uI);"
"lL = 1.5*x2lI*(0.25*rI-0.33333*xlI+0.125*(r2-S2)*rI*x2lI);"
"x2lI = 1.0/(xl*xl);"
"xlI = 1.0/(xl);"
"xuI = 1.0/(xu);"
"x2uI = 1.0/(xu*xu);"
"xu = (r+scaleFactor2);"
"rI = 1.0/(r);"
"r2 = (r*r);"
"xl = factor1*lMax + factor2*xuu + factor3*(r-scaleFactor2);"
"xuu = (r+scaleFactor2);"
"S2 = (scaleFactor2*scaleFactor2);"
"factor1 = step(r-absRadiusScaleDiff);"
"absRadiusScaleDiff = abs(radiusScaleDiff);"
"radiusScaleDiff = (radius1-scaleFactor2);"
"factor2 = step(radius1-scaleFactor2-r);"
"factor3 = step(scaleFactor2-radius1-r);"
"lMax = max(radius1,r-scaleFactor2);"
, CustomGBForce::ParticlePairNoExclusions);
customGbviForce.addComputedValue("B", "(1.0/(radius*radius*radius)-V)^(-0.33333333)", CustomGBForce::SingleParticle);
// nonpolar term + polar self energy
customGbviForce.addEnergyTerm("(-138.935485*0.5*((1.0/soluteDielectric)-(1.0/solventDielectric))*q^2/B)-((1.0/soluteDielectric)-(1.0/solventDielectric))*((gamma*(radius/B)^3))", CustomGBForce::SingleParticle);
// polar pair energy
customGbviForce.addEnergyTerm("-138.935485*(1/soluteDielectric-1/solventDielectric)*q1*q2/f;"
"f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))", CustomGBForce::ParticlePairNoExclusions);
#else
customGbviForce.addComputedValue("V", " uL - lL + factor3/(radius1*radius1*radius1);"
"uL = 1.5*x2uI*(0.25*rI-0.33333*xuI+0.125*(r2-S2)*rI*x2uI);"
"lL = 1.5*x2lI*(0.25*rI-0.33333*xlI+0.125*(r2-S2)*rI*x2lI);"
"x2lI = 1.0/(xl*xl);"
"xlI = 1.0/(xl);"
"xuI = 1.0/(xu);"
"x2uI = 1.0/(xu*xu);"
"xu = (r+scaleFactor2);"
"rI = 1.0/(r);"
"r2 = (r*r);"
"xl = factor1*lMax + factor2*xuu + factor3*(r-scaleFactor2);"
"xuu = (r+scaleFactor2);"
"S2 = (scaleFactor2*scaleFactor2);"
"factor1 = step(r-absRadiusScaleDiff);"
"absRadiusScaleDiff = abs(radiusScaleDiff);"
"radiusScaleDiff = (radius1-scaleFactor2);"
"factor2 = step(radius1-scaleFactor2-r);"
"factor3 = step(scaleFactor2-radius1-r);"
"lMax = max(radius1,r-scaleFactor2);"
, CustomGBForce::ParticlePairNoExclusions);
customGbviForce.addComputedValue("B", "(1.0/(radius*radius*radius)-V)^(-0.33333333)", CustomGBForce::SingleParticle);
// nonpolar term + polar self energy
customGbviForce.addEnergyTerm("(-138.935485*0.5*((1.0/soluteDielectric)-(1.0/solventDielectric))*q^2/B)-((1.0/soluteDielectric)-(1.0/solventDielectric))*((gamma*(radius/B)^3))", CustomGBForce::SingleParticle);
// polar pair energy
customGbviForce.addEnergyTerm("-138.935485*(1/soluteDielectric-1/solventDielectric)*q1*q2/f;"
"f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))", CustomGBForce::ParticlePairNoExclusions);
#endif
// load energies
loadGbviParameters( gbviSoftcoreForce, customGbviForce, log );
return;
}
void testGBVISoftcore( MapStringToDouble inputArgumentMap, FILE* log ){
double lambda1 = 1.0;
double lambda2 = 1.0;
int nonbondedMethod = 0;
int numMolecules = 1;
int numParticlesPerMolecule = 2;
int useQuinticSpline = 1;
int applyAssert = 1;
int positionPlacementMethod = 0;
int serialize = 0;
double boxSize = 10.0;
double relativeTolerance = 1.0e-04;
setDoubleFromMapStringToDouble( inputArgumentMap, "lambda1", lambda1 );
setDoubleFromMapStringToDouble( inputArgumentMap, "lambda2", lambda2 );
setDoubleFromMapStringToDouble( inputArgumentMap, "boxSize", boxSize );
double cutoffDistance = boxSize*0.4;;
setDoubleFromMapStringToDouble( inputArgumentMap, "cutoffDistance", cutoffDistance);
setDoubleFromMapStringToDouble( inputArgumentMap, "relativeTolerance", relativeTolerance );
setIntFromMapStringToDouble( inputArgumentMap, "positionPlacementMethod", positionPlacementMethod ) ;
setIntFromMapStringToDouble( inputArgumentMap, "nonbondedMethod", nonbondedMethod );
setIntFromMapStringToDouble( inputArgumentMap, "numMolecules", numMolecules );
setIntFromMapStringToDouble( inputArgumentMap, "numParticlesPerMolecule", numParticlesPerMolecule );
setIntFromMapStringToDouble( inputArgumentMap, "serialize", serialize );
int numParticles = numMolecules*numParticlesPerMolecule;
int includeGbvi = 1;
double reactionFieldDielectric = 80.0;
if( log ){
double particleDensity = static_cast<double>(numParticles)/(boxSize*boxSize*boxSize);
double particleCube = pow( particleDensity, (-1.0/3.0) );
(void) fprintf( log, "\n--------------------------------------------------------------------------------------\n" );
(void) fprintf( log, "Input arguments\n" );
(void) fflush( log );
(void) fprintf( log, " includeGbvi %d\n", includeGbvi );
(void) fprintf( log, " nonbondedMethod %d\n", nonbondedMethod );
(void) fprintf( log, " numParticles %d\n", numParticles );
(void) fprintf( log, " numMolecules %d\n", numMolecules );
(void) fprintf( log, " numParticlesPerMolecule %d\n", numParticlesPerMolecule );
(void) fprintf( log, " useQuinticSpline %d\n", useQuinticSpline );
(void) fprintf( log, " positionPlacementMethod %d\n", positionPlacementMethod);
#ifdef USE_SOFTCORE
(void) fprintf( log, " lambda1 %8.3f\n", lambda1 );
(void) fprintf( log, " lambda2 %8.3f\n", lambda2 );
#endif
(void) fprintf( log, " boxSize %8.3f\n", boxSize );
(void) fprintf( log, " cutoffDistance %8.3f\n", cutoffDistance );
(void) fprintf( log, " reactionFieldDielectric %8.3f\n", reactionFieldDielectric );
(void) fprintf( log, " relativeTolerance %8.1e\n", relativeTolerance );
(void) fprintf( log, " particleDensity %8.2e\n", particleDensity );
(void) fprintf( log, " particleCube %8.2e\n", particleCube );
}
// set positions
vector<Vec3> positions(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
PositionGenerator positionGenerator( numMolecules, numParticlesPerMolecule, boxSize );
if( log ){
positionGenerator.setLog( log );
}
if( positionPlacementMethod == 1 ){
positionGenerator.setPositions( PositionGenerator::SimpleGrid, sfmt, positions );
} else {
positionGenerator.setBondDistance( 0.3 );
positionGenerator.setPositions( PositionGenerator::Random, sfmt, positions );
}
// create GBSAGBVISoftcoreForce and populate w/ parameters
#ifdef USE_SOFTCORE
GBVISoftcoreForce* gbviSoftcoreForce = new GBVISoftcoreForce();
#else
GBVIForce* gbviSoftcoreForce = new GBVIForce();
#endif
gbviSoftcoreForce->setSolventDielectric( 78.3 );
//gbviSoftcoreForce.setSolventDielectric( 1.0e+10 );
//gbviSoftcoreForce.setSolventDielectric( 1.0 );
gbviSoftcoreForce->setSoluteDielectric( 1.0 );
gbviSoftcoreForce->setCutoffDistance( cutoffDistance );
const int numberOfParameters = 4;
const int ChargeIndex = 0;
const int SigmaIndex = 1;
const int GammaIndex = 2;
const int LambdaIndex = 3;
std::vector<double> parameterLowerBound( numberOfParameters, 0.0 );
double fixedCharge = 1.0;
parameterLowerBound[ChargeIndex] = fixedCharge; // charge
parameterLowerBound[SigmaIndex] = 0.1; // sigma
parameterLowerBound[GammaIndex] = 0.1; // gamma
parameterLowerBound[LambdaIndex] = lambda1; // lambda
std::vector<double> parameterUpperBound( parameterLowerBound );
parameterUpperBound[ChargeIndex] = fixedCharge; // charge
parameterUpperBound[SigmaIndex] = 0.3; // sigma
parameterUpperBound[GammaIndex] = 40.0; // gamma
std::vector<double> parameters( numberOfParameters );
double charge = fixedCharge;
for( int ii = 0; ii < numMolecules; ii++) {
charge *= -1.0;
double lambda = ii < (numMolecules/2) ? lambda1 : lambda2;
randomizeParameters( parameterLowerBound, parameterUpperBound, sfmt, parameters );
#ifdef USE_SOFTCORE
gbviSoftcoreForce->addParticle( charge, parameters[SigmaIndex], parameters[GammaIndex], lambda );
#else
gbviSoftcoreForce->addParticle( charge, parameters[SigmaIndex], parameters[GammaIndex] );
#endif
int baseParticleIndex = ii*numParticlesPerMolecule;
for( int jj = 1; jj < numParticlesPerMolecule; jj++) {
// alternate charges
charge *= -1.0;
randomizeParameters( parameterLowerBound, parameterUpperBound, sfmt, parameters );
#ifdef USE_SOFTCORE
gbviSoftcoreForce->addParticle( charge, parameters[SigmaIndex], parameters[GammaIndex], lambda );
#else
gbviSoftcoreForce->addParticle( charge, parameters[SigmaIndex], parameters[GammaIndex] );
#endif
#if IMPLICIT_SOLVENT == GBVI_FLAG
double bondDistance = positionGenerator.getDistance( baseParticleIndex, baseParticleIndex+jj, positions );
gbviSoftcoreForce->addBond( baseParticleIndex, baseParticleIndex+jj, bondDistance );
#endif
}
// alternate charge if numParticlesPerMolecule is odd
if( (numParticlesPerMolecule % 2) ){
charge *= -1.0;
}
}
// Create system
System standardSystem;
for (int i = 0; i < numParticles; i++) {
standardSystem.addParticle(1.0);
}
standardSystem.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
standardSystem.addForce(gbviSoftcoreForce);
// copy system and forces
System systemCopy;
copySystem( standardSystem, systemCopy );
CustomGBForce* customGbviForce = new CustomGBForce();
createCustomGBVI( *customGbviForce, *gbviSoftcoreForce, log );
systemCopy.addForce(customGbviForce);
// perform comparison
std::stringstream idString;
idString << "Nb " << nonbondedMethod << " l2 " << std::fixed << setprecision(2) << lambda2;
runSystemComparisonTest( standardSystem, systemCopy, "Reference", "Reference", positions, inputArgumentMap, idString.str(), log );
// serialize
std::stringstream baseFileName;
if( serialize ){
baseFileName << "_N" << positions.size();
baseFileName << "_Nb" << nonbondedMethod;
serializeSystemAndPositions( standardSystem, positions, baseFileName.str(), log);
}
}
int main() {
try {
testSingleParticle();
testEnergyEthaneSwitchingFunction( 0 );
testEnergyEthaneSwitchingFunction( 1 );
}
catch(const exception& e) {
//FILE* log = stderr;
FILE* log = NULL;
//testSingleParticle();
//testEnergyEthaneSwitchingFunction( 0 );
//testEnergyEthaneSwitchingFunction( 1 );
MapStringToDouble inputArgumentMap;
inputArgumentMap["lambda2"] = 0.5;
inputArgumentMap["nonbondedMethod"] = 0;
inputArgumentMap["numMolecules"] = 10;
inputArgumentMap["boxSize"] = 5.0;
inputArgumentMap["positionPlacementMethod"] = 0;
inputArgumentMap["cutoffDistance"] = 0.3*inputArgumentMap["boxSize"];
//inputArgumentMap["cutoffDistance"] = 1.0;
inputArgumentMap["relativeTolerance"] = 5.0e-04;
inputArgumentMap["serialize"] = 1;
//inputArgumentMap["numParticlesPerMolecule"] = 2;
testGBVISoftcore( inputArgumentMap, log );
} catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
......
plugins/freeEnergy/platforms/reference/tests/TestReferenceOBCSoftcoreForce.cpp
View file @ 032e18de
......@@ -77,7 +77,7 @@ void testOBCSoftcore( double lambda1, double lambda2 ){
custom->addGlobalParameter("solventDielectric", obc->getSolventDielectric());
custom->addGlobalParameter("soluteDielectric", obc->getSoluteDielectric());
custom->addComputedValue("I", "lambda1*lambda2*step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(1/U^2-1/L^2)*(r-sr2*sr2/r)+0.5*log(L/U)/r+C);"
custom->addComputedValue("I", "lambda2*step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(1/U^2-1/L^2)*(r-sr2*sr2/r)+0.5*log(L/U)/r+C);"
"U=r+sr2;"
"C=2*(1/or1-1/L)*step(sr2-r-or1);"
"L=max(or1, D);"
......@@ -86,8 +86,8 @@ void testOBCSoftcore( double lambda1, double lambda2 ){
"or1 = radius1-0.009; or2 = radius2-0.009", CustomGBForce::ParticlePairNoExclusions);
custom->addComputedValue("B", "1/(1/or-tanh(1*psi-0.8*psi^2+4.85*psi^3)/radius);"
"psi=I*or; or=radius-0.009", CustomGBForce::SingleParticle);
custom->addEnergyTerm("lambda*28.3919551*(radius+0.14)^2*(radius/B)^6-lambda*lambda*0.5*138.935485*(1/soluteDielectric-1/solventDielectric)*q^2/B", CustomGBForce::SingleParticle);
custom->addEnergyTerm("-138.935485*lambda1*lambda2*(1/soluteDielectric-1/solventDielectric)*q1*q2/f;"
custom->addEnergyTerm("lambda*28.3919551*(radius+0.14)^2*(radius/B)^6-0.5*138.935485*(1/soluteDielectric-1/solventDielectric)*q^2/B", CustomGBForce::SingleParticle);
custom->addEnergyTerm("-138.935485*(1/soluteDielectric-1/solventDielectric)*q1*q2/f;"
"f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))", CustomGBForce::ParticlePairNoExclusions);
vector<Vec3> positions(numParticles);
......@@ -103,14 +103,14 @@ void testOBCSoftcore( double lambda1, double lambda2 ){
obc->addParticle( charge*lambda1, 0.2, 0.5, lambda1);
obc->addParticle(-charge*lambda1, 0.1, 0.5, lambda1);
params[0] = charge;
params[0] = charge*lambda1;
params[1] = 0.2;
params[2] = 0.5;
params[3] = lambda1;
custom->addParticle(params);
params[0] = -charge;
params[0] = -charge*lambda1;
params[1] = 0.1;
custom->addParticle(params);
......@@ -119,13 +119,13 @@ void testOBCSoftcore( double lambda1, double lambda2 ){
obc->addParticle( charge*lambda2, 0.2, 0.8, lambda2);
obc->addParticle(-charge*lambda2, 0.1, 0.8, lambda2);
params[0] = charge;
params[0] = charge*lambda2;
params[1] = 0.2;
params[2] = 0.8;
params[3] = lambda2;
custom->addParticle(params);
params[0] = -charge;
params[0] = -charge*lambda2;
params[1] = 0.1;
custom->addParticle(params);
}
......
plugins/freeEnergy/platforms/cuda/tests/TestCudaSoftcoreForce.h → plugins/freeEnergy/platforms/reference/tests/TestReferenceSoftcoreForce.h
View file @ 032e18de
......@@ -59,10 +59,7 @@
#include <iostream>
#include <vector>
#include <algorithm>
#include <iostream>
#include <cstdio>
#include <vector>
#include <typeinfo>
extern "C" void registerFreeEnergyCudaKernelFactories();
......@@ -1322,29 +1319,27 @@ static NonbondedForce* copyNonbondedForce( const NonbondedForce& nonbondedForce
*
*/
static System* copySystem( const System& inputSystem ){
System* systemCopy = new System();
static void copySystem( const System& inputSystem, System& systemCopy ){
for( unsigned int ii = 0; ii < inputSystem.getNumParticles(); ii++ ){
systemCopy->addParticle( inputSystem.getParticleMass( static_cast<int>(ii) ) );
systemCopy.addParticle( inputSystem.getParticleMass( static_cast<int>(ii) ) );
}
Vec3 a;
Vec3 b;
Vec3 c;
inputSystem.getDefaultPeriodicBoxVectors( a, b, c );
systemCopy->setDefaultPeriodicBoxVectors( a, b, c );
systemCopy.setDefaultPeriodicBoxVectors( a, b, c );
for( unsigned int ii = 0; ii < inputSystem.getNumConstraints(); ii++ ){
int index;
int particle1, particle2;
double distance;
inputSystem.getConstraintParameters( ii, particle1, particle2, distance);
systemCopy->addConstraint( particle1, particle2, distance);
systemCopy.addConstraint( particle1, particle2, distance);
}
return systemCopy;
return;
}
/**
......@@ -1809,8 +1804,9 @@ void runSystemComparisonTest( System& system1, System& system2,
(void) fprintf( log, "System1: particles=%d forces=%d System2: particles=%d forces=%d\n",
system1.getNumParticles(), system1.getNumForces(),
system2.getNumParticles(), system2.getNumForces() );
(void) fprintf( log, "Positions=%u\n",
static_cast<unsigned int>(positions.size()) );
(void) fprintf( log, "Positions=%u\n", static_cast<unsigned int>(positions.size()) );
(void) fprintf( log, "Platform1=%s Platform2=%s\n", platform1.c_str(), platform2.c_str() );
(void) fprintf( log, "relativeTolerance=%8.2e applyAssert=%d\n", relativeTolerance, applyAssert );
MapStringInt stringForceVector1;
MapStringInt stringForceVector2;
......@@ -1836,7 +1832,7 @@ void runSystemComparisonTest( System& system1, System& system2,
if( system1.getNumParticles() != static_cast<int>(positions.size()) ){
std::stringstream msg;
msg << "Number of partciles for system des not equal size of position array: " << system1.getNumParticles() << " != " << positions.size();
msg << "Number of particles for system does not equal size of position array: " << system1.getNumParticles() << " != " << positions.size();
throw OpenMMException( msg.str() );
}
......@@ -1844,7 +1840,7 @@ void runSystemComparisonTest( System& system1, System& system2,
context1.setPositions(positions);
State state1 = context1.getState(State::Forces | State::Energy);
Context context2( system2, integrator2, Platform::getPlatformByName( "Reference"));
Context context2( system2, integrator2, Platform::getPlatformByName( platform2 ));
context2.setPositions(positions);
State state2 = context2.getState(State::Forces | State::Energy);
......
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