Skip to content

GitLab

Commit 80d8311e authored by Mark Friedrichs's avatar Mark Friedrichs
Browse files

Improved accuracy of OutOfPlaneBend force

parent 2af97193
Hide whitespace changes
Inline Side-by-side
Showing with 35 additions and 46 deletions
plugins/amoeba/platforms/cuda/src/kernels/AmoebaGpu.cpp
View file @ 80d8311e
......@@ -1033,7 +1033,14 @@ exit(0);
extern "C"
void gpuSetAmoebaBondOffsets(amoebaGpuContext amoebaGpu )
{
// make sure only flip once
static int flipped = 0;
if( amoebaGpu && flipped ){
return;
}
flipped = 1;
_gpuContext* gpu = amoebaGpu->gpuContext;
......
plugins/amoeba/platforms/cuda/src/kernels/amoebaGpuTypes.h
View file @ 80d8311e
......@@ -343,6 +343,9 @@ void gpuSetAmoebaSASAParameters( amoebaGpuContext amoebaGpu , float probeRadius,
extern "C"
void amoebaGpuSetConstants(amoebaGpuContext gpu);
extern "C"
void gpuSetAmoebaBondOffsets(amoebaGpuContext gpu);
/*
extern "C"
void gpuSetDihedralParameters(gpuContext gpu, const std::vector<int>& atom1, const std::vector<int>& atom2, const std::vector<int>& atom3, const std::vector<int>& atom4,
......
plugins/amoeba/platforms/cuda/src/kernels/kCalculateAmoebaLocalForces.cu
View file @ 80d8311e
......@@ -1172,51 +1172,21 @@ void kCalculateAmoebaLocalForces_kernel()
float bkk2 = (cc != 0.0f ) ? (ee*ee)/(cc) : 0.0f;
float bkk3 = (rdb2 != 0.0f ) ? bkk2/rdb2 : 0.0f;
bkk2 = rdb2 - bkk2;
float cosine;
if( fabs( bkk3 ) < 0.98f ){
bkk3 = 1.0f - bkk3;
cosine = bkk3 > 0.0f ? sqrtf(bkk3) : 0.0f;
cosine = (cosine > 1.0f) ? 0.0f : acos(cosine);
float angle;
if( bkk3 < 1.0e-6f ){
angle = sqrtf( bkk3 );
} else if( bkk3 < 0.98f ){
float cosine = sqrtf(1.0f - bkk3);
angle = acos(cosine);
} else {
cosine = sqrtf(bkk3);
cosine = asin(cosine);
float sin = sqrtf(bkk3);
angle = asin(sin);
}
/*
c
c W-D-C angle between A-B-C plane and B-D vector for D-B<AC
c
if (opbtyp .eq. 'W-D-C') then
rab2 = xab*xab + yab*yab + zab*zab
rcb2 = xcb*xcb + ycb*ycb + zcb*zcb
dot = xab*xcb+yab*ycb+zab*zcb
cc = rab2*rcb2 - dot*dot
c
c Allinger angle between A-C-D plane and D-B vector for D-B<AC
c
else if (opbtyp .eq. 'ALLINGER') then
rad2 = xad*xad + yad*yad + zad*zad
rcd2 = xcd*xcd + ycd*ycd + zcd*zcd
dot = xad*xcd + yad*ycd + zad*zcd
cc = rad2*rcd2 - dot*dot
end if
c
c find the out-of-plane angle bending energy
c
ee = xdb*(yab*zcb-zab*ycb) + ydb*(zab*xcb-xab*zcb)
& + zdb*(xab*ycb-yab*xcb)
rdb2 = xdb*xdb + ydb*ydb + zdb*zdb
if (rdb2.ne.0.0d0 .and. cc.ne.0.0d0) then
bkk2 = rdb2 - ee*ee/cc
cosine = sqrt(bkk2/rdb2)
cosine = min(1.0d0,max(-1.0d0,cosine))
angle = radian * acos(cosine)
*/
// find the out-of-plane energy and master chain rule terms
float dt = LOCAL_HACK_RADIAN_D*cosine;
float dt = LOCAL_HACK_RADIAN_D*angle;
float dt2 = dt * dt;
float dt3 = dt2 * dt;
float dt4 = dt2 * dt2;
......@@ -1284,7 +1254,7 @@ c
float dedyib = -dedyia - dedyic - dedyid;
float dedzib = -dedzia - dedzic - dedzid;
// increment the out-of-plane bending energy and gradient;
// increment the out-of-plane bending gradient
int4 atom2 = cAmoebaSim.pAmoebaOutOfPlaneBendID2[pos1];
......
plugins/amoeba/platforms/cuda/tests/AmoebaTinkerParameterFile.cpp
View file @ 80d8311e
......@@ -4808,12 +4808,16 @@ void testUsingAmoebaTinkerParameterFile( const std::string& amoebaTinkerParamete
for( unsigned int ii = 0; ii < fileList.size(); ii++ ){
FILE* filePtr = fileList[ii];
(void) fprintf( filePtr, "\n" );
(void) fprintf( filePtr, "%s: conversion factors %15.7e %15.7e tolerance=%15.7e %s\n", methodName.c_str(), energyConversion, forceConversion, tolerance, amoebaTinkerParameterFileName.c_str() );
(void) fprintf( filePtr, "%s: conversion factors %15.7e %15.7e tolerance=%15.7e %s\n",
methodName.c_str(), energyConversion, forceConversion, tolerance, amoebaTinkerParameterFileName.c_str() );
double deltaE = fabs( expectedEnergy - energyConversion*state.getPotentialEnergy());
double denom = fabs( expectedEnergy ) + fabs( energyConversion*state.getPotentialEnergy());
if( denom > 0.0 )deltaE *= 2.0/denom;
(void) fprintf( filePtr, "expectedE %10.3e %15.7e %15.7e %20s %30s\n", deltaE, expectedEnergy, energyConversion*state.getPotentialEnergy(), amoebaTinkerParameterFileName.c_str(), activeForceNames.c_str() );
(void) fprintf( filePtr, "%s: %u %u Active forces: %s\n", methodName.c_str(), expectedForces.size(), forces.size(), activeForceNames.c_str() );
(void) fprintf( filePtr, "expectedE %10.3e %15.7e %15.7e %20s %30s\n",
deltaE, expectedEnergy, energyConversion*state.getPotentialEnergy(),
amoebaTinkerParameterFileName.c_str(), activeForceNames.c_str() );
(void) fprintf( filePtr, "%s: %u %u Active forces: %s\n",
methodName.c_str(), expectedForces.size(), forces.size(), activeForceNames.c_str() );
double maxRelativeDelta = -1.0e+30;
unsigned int maxRelativeDeltaIndex = -1;
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
......@@ -4824,6 +4828,7 @@ void testUsingAmoebaTinkerParameterFile( const std::string& amoebaTinkerParamete
double sumNorms = 0.5*(normF1 + normF2);
double relativeDelta = sumNorms > 0.0 ? fabs( normF1 - normF2 )/sumNorms : 0.0;
bool badMatch = (cutoffDelta < relativeDelta) && (sumNorms > 0.1) ? true : false;
badMatch = badMatch || (normF1 == 0.0 && normF2 > 0.0) || (normF2 == 0.0 && normF1 > 0.0);
if( badMatch || showAll ){
(void) fprintf( filePtr, "%6u %10.3e %10.3e [%15.7e %15.7e %15.7e] [%15.7e %15.7e %15.7e] %s\n", ii, relativeDelta, delta,
expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2], forceConversion*forces[ii][0], forceConversion*forces[ii][1], forceConversion*forces[ii][2], ( (showAll && badMatch) ? " XXX" : "") );
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaHarmonicBondForce.cpp
View file @ 80d8311e
......@@ -52,8 +52,8 @@ static void computeAmoebaHarmonicBondForce(int bondIndex, std::vector<Vec3>& po
int particle1, particle2;
double bondLength;
double quadraticK;
double cubicK;
double quarticK;
double cubicK = amoebaHarmonicBondForce.getAmoebaGlobalHarmonicBondCubic();
double quarticK = amoebaHarmonicBondForce.getAmoebaGlobalHarmonicBondQuartic();
amoebaHarmonicBondForce.getBondParameters(bondIndex, particle1, particle2, bondLength, quadraticK );
double deltaR[3];
......@@ -152,6 +152,8 @@ void testOneBond( FILE* log ) {
double quadraticK = 1.0;
double cubicK = 2.0;
double quarticicK = 3.0;
amoebaHarmonicBondForce->setAmoebaGlobalHarmonicBondCubic( cubicK );
amoebaHarmonicBondForce->setAmoebaGlobalHarmonicBondQuartic( quarticicK );
amoebaHarmonicBondForce->addBond(0, 1, bondLength, quadraticK);
system.addForce(amoebaHarmonicBondForce);
......@@ -181,6 +183,8 @@ void testTwoBond( FILE* log ) {
double quadraticK = 1.0;
double cubicK = 2.0;
double quarticicK = 3.0;
amoebaHarmonicBondForce->setAmoebaGlobalHarmonicBondCubic( cubicK );
amoebaHarmonicBondForce->setAmoebaGlobalHarmonicBondQuartic( quarticicK );
amoebaHarmonicBondForce->addBond(0, 1, bondLength, quadraticK);
amoebaHarmonicBondForce->addBond(1, 2, bondLength, quadraticK);
......
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