Skip to content

GitLab

Commit 0843c5f3 authored by Peter Eastman's avatar Peter Eastman
Browse files

First stage of a refactoring to clean up the reference platform

parent bacc1eff
Hide whitespace changes
Inline Side-by-side
Showing with 38 additions and 813 deletions
plugins/freeEnergy/platforms/reference/src/SimTKReference/ReferenceFreeEnergyLJCoulombSoftcoreIxn.cpp
View file @ 0843c5f3
......@@ -39,6 +39,7 @@
#include "openmm/internal/MSVC_erfc.h"
using std::vector;
using OpenMM::RealVec;
/**---------------------------------------------------------------------------------------
......@@ -238,9 +239,9 @@ int ReferenceFreeEnergyLJCoulombSoftcoreIxn::getDerivedParameters( RealOpenMM c6
--------------------------------------------------------------------------------------- */
int ReferenceFreeEnergyLJCoulombSoftcoreIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** atomCoordinates,
int ReferenceFreeEnergyLJCoulombSoftcoreIxn::calculateEwaldIxn( int numberOfAtoms, vector<RealVec>& atomCoordinates,
RealOpenMM** atomParameters, int** exclusions,
RealOpenMM* fixedParameters, RealOpenMM** forces,
RealOpenMM* fixedParameters, vector<RealVec>& forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy) const {
#if 0
typedef std::complex<RealOpenMM> d_complex;
......@@ -546,9 +547,9 @@ int ReferenceFreeEnergyLJCoulombSoftcoreIxn::calculateEwaldIxn( int numberOfAtom
--------------------------------------------------------------------------------------- */
int ReferenceFreeEnergyLJCoulombSoftcoreIxn::calculatePMEIxn( int numberOfAtoms, RealOpenMM** atomCoordinates,
int ReferenceFreeEnergyLJCoulombSoftcoreIxn::calculatePMEIxn( int numberOfAtoms, vector<RealVec>& atomCoordinates,
RealOpenMM** atomParameters, int** exclusions,
RealOpenMM* fixedParameters, RealOpenMM** forces,
RealOpenMM* fixedParameters, vector<RealVec>& forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy) const {
......@@ -678,9 +679,9 @@ int ReferenceFreeEnergyLJCoulombSoftcoreIxn::calculatePMEIxn( int numberOfAtoms,
--------------------------------------------------------------------------------------- */
int ReferenceFreeEnergyLJCoulombSoftcoreIxn::calculatePairIxn( int numberOfAtoms, RealOpenMM** atomCoordinates,
int ReferenceFreeEnergyLJCoulombSoftcoreIxn::calculatePairIxn( int numberOfAtoms, vector<RealVec>& atomCoordinates,
RealOpenMM** atomParameters, int** exclusions,
RealOpenMM* fixedParameters, RealOpenMM** forces,
RealOpenMM* fixedParameters, vector<RealVec>& forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const {
if (ewald || pme)
......@@ -739,8 +740,8 @@ int ReferenceFreeEnergyLJCoulombSoftcoreIxn::calculatePairIxn( int numberOfAtoms
--------------------------------------------------------------------------------------- */
int ReferenceFreeEnergyLJCoulombSoftcoreIxn::calculateOneIxn( int ii, int jj, RealOpenMM** atomCoordinates,
RealOpenMM** atomParameters, RealOpenMM** forces,
int ReferenceFreeEnergyLJCoulombSoftcoreIxn::calculateOneIxn( int ii, int jj, vector<RealVec>& atomCoordinates,
RealOpenMM** atomParameters, vector<RealVec>& forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const {
// ---------------------------------------------------------------------------------------
......@@ -826,59 +827,6 @@ int ReferenceFreeEnergyLJCoulombSoftcoreIxn::calculateOneIxn( int ii, int jj, Re
energyByAtom[jj] += energy;
}
}
// debug
if( debug == ii ){
static bool printHeader = false;
std::stringstream message;
message << methodName;
message << std::endl;
int pairArray[2] = { ii, jj };
if( !printHeader ){
printHeader = true;
message << std::endl;
message << methodName.c_str() << " a0 k [c q p s] r1 r2 angle dt rp p[] dot cosine angle dEdR*r F[]" << std::endl;
}
message << std::endl;
for( int kk = 0; kk < 2; kk++ ){
message << " Atm " << pairArray[kk] << " [" << atomCoordinates[pairArray[kk]][0] << " " << atomCoordinates[pairArray[kk]][1] << " " << atomCoordinates[pairArray[kk]][2] << "] ";
}
message << std::endl << " Delta:";
for( int kk = 0; kk < (LastAtomIndex - 1); kk++ ){
message << " [";
for( int jj = 0; jj < ReferenceForce::LastDeltaRIndex; jj++ ){
message << deltaR[kk][jj] << " ";
}
message << "]";
}
message << std::endl;
for( int kk = 0; kk < 2; kk++ ){
message << " p" << pairArray[kk] << " [";
message << atomParameters[pairArray[kk]][0] << " " << atomParameters[pairArray[kk]][1] << " " << atomParameters[pairArray[kk]][2];
message << "]";
}
message << std::endl;
message << " dEdR=" << dEdR;
message << " E=" << energy << " force factors: ";
message << "F=compute force; f=cumulative force";
message << std::endl << " ";
message << " f" << ii << "[";
SimTKOpenMMUtilities::formatRealStringStream( message, deltaR[0], threeI, dEdR );
message << "]";
for( int kk = 0; kk < 2; kk++ ){
message << " F" << pairArray[kk] << " [";
SimTKOpenMMUtilities::formatRealStringStream( message, forces[pairArray[kk]], threeI );
message << "]";
}
//SimTKOpenMMLog::printMessage( message );
}
return ReferenceForce::DefaultReturn;
}
......
plugins/freeEnergy/platforms/reference/src/SimTKReference/ReferenceFreeEnergyLJCoulombSoftcoreIxn.h
View file @ 0843c5f3
......@@ -26,6 +26,7 @@
#include "SimTKReference/ReferencePairIxn.h"
#include "SimTKReference/ReferenceNeighborList.h"
#include <vector>
// ---------------------------------------------------------------------------------------
......@@ -68,8 +69,8 @@ class ReferenceFreeEnergyLJCoulombSoftcoreIxn : public ReferencePairIxn {
--------------------------------------------------------------------------------------- */
int calculateOneIxn( int atom1, int atom2, RealOpenMM** atomCoordinates,
RealOpenMM** atomParameters, RealOpenMM** forces,
int calculateOneIxn( int atom1, int atom2, std::vector<OpenMM::RealVec>& atomCoordinates,
RealOpenMM** atomParameters, std::vector<OpenMM::RealVec>& forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const;
......@@ -194,9 +195,9 @@ class ReferenceFreeEnergyLJCoulombSoftcoreIxn : public ReferencePairIxn {
--------------------------------------------------------------------------------------- */
int calculatePairIxn( int numberOfAtoms, RealOpenMM** atomCoordinates,
int calculatePairIxn( int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
RealOpenMM** atomParameters, int** exclusions,
RealOpenMM* fixedParameters, RealOpenMM** forces,
RealOpenMM* fixedParameters, std::vector<OpenMM::RealVec>& forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const;
private:
......@@ -220,9 +221,9 @@ private:
--------------------------------------------------------------------------------------- */
int calculateEwaldIxn( int numberOfAtoms, RealOpenMM** atomCoordinates,
int calculateEwaldIxn( int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
RealOpenMM** atomParameters, int** exclusions,
RealOpenMM* fixedParameters, RealOpenMM** forces,
RealOpenMM* fixedParameters, std::vector<OpenMM::RealVec>& forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const;
/**---------------------------------------------------------------------------------------
......@@ -245,9 +246,9 @@ private:
--------------------------------------------------------------------------------------- */
int calculatePMEIxn( int numberOfAtoms, RealOpenMM** atomCoordinates,
int calculatePMEIxn( int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
RealOpenMM** atomParameters, int** exclusions,
RealOpenMM* fixedParameters, RealOpenMM** forces,
RealOpenMM* fixedParameters, std::vector<OpenMM::RealVec>& forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const;
/**---------------------------------------------------------------------------------------
......
plugins/freeEnergy/platforms/reference/src/gbsa/CpuGBVISoftcore.cpp
View file @ 0843c5f3
......@@ -32,6 +32,9 @@
#include "../SimTKReference/ReferenceForce.h"
#include <math.h>
using std::vector;
using OpenMM::RealVec;
/**---------------------------------------------------------------------------------------
CpuGBVISoftcore constructor
......@@ -319,7 +322,7 @@ int CpuGBVISoftcore::computeBornRadiiUsingQuinticSpline( RealOpenMM atomicRadius
#define GBVISoftcoreDebug 0
int CpuGBVISoftcore::computeBornRadii( RealOpenMM** atomCoordinates, RealOpenMM* bornRadii, RealOpenMM* switchDeriviative ){
int CpuGBVISoftcore::computeBornRadii( vector<RealVec>& atomCoordinates, RealOpenMM* bornRadii, RealOpenMM* switchDeriviative ){
// ---------------------------------------------------------------------------------------
......@@ -604,7 +607,7 @@ RealOpenMM CpuGBVISoftcore::Sgb( RealOpenMM t ){
--------------------------------------------------------------------------------------- */
RealOpenMM CpuGBVISoftcore::computeBornEnergy( const RealOpenMM* bornRadii, RealOpenMM** atomCoordinates,
RealOpenMM CpuGBVISoftcore::computeBornEnergy( const RealOpenMM* bornRadii, vector<RealVec>& atomCoordinates,
const RealOpenMM* partialCharges ){
// ---------------------------------------------------------------------------------------
......@@ -720,8 +723,8 @@ RealOpenMM e3 = -partialChargeI2*partialCharges[atomJ]*Sgb( t )/deltaR[Reference
--------------------------------------------------------------------------------------- */
int CpuGBVISoftcore::computeBornForces( const RealOpenMM* bornRadii, RealOpenMM** atomCoordinates,
const RealOpenMM* partialCharges, RealOpenMM** inputForces ){
int CpuGBVISoftcore::computeBornForces( const RealOpenMM* bornRadii, vector<RealVec>& atomCoordinates,
const RealOpenMM* partialCharges, vector<RealVec>& inputForces ){
// ---------------------------------------------------------------------------------------
......
plugins/freeEnergy/platforms/reference/src/gbsa/CpuGBVISoftcore.h
View file @ 0843c5f3
......@@ -115,7 +115,7 @@ class CpuGBVISoftcore : public CpuImplicitSolvent {
--------------------------------------------------------------------------------------- */
int computeBornRadii( RealOpenMM** atomCoordinates, RealOpenMM* bornRadii,
int computeBornRadii( std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM* bornRadii,
RealOpenMM* switchDeriviative = NULL );
/**---------------------------------------------------------------------------------------
......@@ -210,7 +210,7 @@ class CpuGBVISoftcore : public CpuImplicitSolvent {
--------------------------------------------------------------------------------------- */
RealOpenMM computeBornEnergy( const RealOpenMM* bornRadii, RealOpenMM** atomCoordinates,
RealOpenMM computeBornEnergy( const RealOpenMM* bornRadii, std::vector<OpenMM::RealVec>& atomCoordinates,
const RealOpenMM* partialCharges );
/**---------------------------------------------------------------------------------------
......@@ -226,8 +226,8 @@ class CpuGBVISoftcore : public CpuImplicitSolvent {
--------------------------------------------------------------------------------------- */
int computeBornForces( const RealOpenMM* bornRadii, RealOpenMM** atomCoordinates,
const RealOpenMM* partialCharges, RealOpenMM** inputForces );
int computeBornForces( const RealOpenMM* bornRadii, std::vector<OpenMM::RealVec>& atomCoordinates,
const RealOpenMM* partialCharges, std::vector<OpenMM::RealVec>& inputForces );
/**---------------------------------------------------------------------------------------
......
plugins/freeEnergy/platforms/reference/src/gbsa/CpuObcSoftcore.cpp
View file @ 0843c5f3
......@@ -33,6 +33,9 @@
#include <cmath>
#include <cstdio>
using std::vector;
using OpenMM::RealVec;
/**---------------------------------------------------------------------------------------
CpuObcSoftcore constructor
......@@ -216,7 +219,7 @@ RealOpenMM* CpuObcSoftcore::getObcChainTemp( void ){
--------------------------------------------------------------------------------------- */
int CpuObcSoftcore::computeBornRadii( RealOpenMM** atomCoordinates, RealOpenMM* bornRadii, RealOpenMM* obcChain ){
int CpuObcSoftcore::computeBornRadii( vector<RealVec>& atomCoordinates, RealOpenMM* bornRadii, RealOpenMM* obcChain ){
// ---------------------------------------------------------------------------------------
......@@ -426,8 +429,8 @@ int CpuObcSoftcore::computeAceNonPolarForce( const ObcSoftcoreParameters* obcSof
--------------------------------------------------------------------------------------- */
int CpuObcSoftcore::computeBornEnergyForces( RealOpenMM* bornRadii, RealOpenMM** atomCoordinates,
const RealOpenMM* partialCharges, RealOpenMM** inputForces ){
int CpuObcSoftcore::computeBornEnergyForces( RealOpenMM* bornRadii, vector<RealVec>& atomCoordinates,
const RealOpenMM* partialCharges, vector<RealVec>& inputForces ){
// ---------------------------------------------------------------------------------------
......@@ -738,662 +741,3 @@ int CpuObcSoftcore::computeBornEnergyForces( RealOpenMM* bornRadii, RealOpenMM**
return 0;
}
/**---------------------------------------------------------------------------------------
Get string w/ state
@param title title (optional)
@return string containing state
--------------------------------------------------------------------------------------- */
std::string CpuObcSoftcore::getStateString( const char* title ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuImplicitSolvent::getStateString";
// ---------------------------------------------------------------------------------------
std::stringstream message;
message << CpuImplicitSolvent::getStateString( title );
return message.str();
}
/**---------------------------------------------------------------------------------------
Write Born energy and forces (Simbios)
@param atomCoordinates atomic coordinates
@param partialCharges partial charges
@param forces forces
@param resultsFileName output file name
@return 0 unless
file cannot be opened
in which case return -1
--------------------------------------------------------------------------------------- */
int CpuObcSoftcore::writeBornEnergyForces( RealOpenMM** atomCoordinates,
const RealOpenMM* partialCharges, RealOpenMM** forces,
const std::string& resultsFileName ) const {
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nCpuObcSoftcore::writeBornEnergyForces";
// ---------------------------------------------------------------------------------------
ImplicitSolventParameters* implicitSolventParameters = getImplicitSolventParameters();
const ObcSoftcoreParameters* obcSoftcoreParameters = static_cast<const ObcSoftcoreParameters*>(implicitSolventParameters);
int numberOfAtoms = obcSoftcoreParameters->getNumberOfAtoms();
const RealOpenMM* atomicRadii = obcSoftcoreParameters->getAtomicRadii();
const RealOpenMM* bornRadii = getBornRadiiConst();
const RealOpenMM* scaledRadii = obcSoftcoreParameters->getScaledRadiusFactors();
const RealOpenMM* obcChain = getObcChainConst();
const RealOpenMM energy = getEnergy();
// ---------------------------------------------------------------------------------------
// open file -- return if unsuccessful
FILE* implicitSolventResultsFile = NULL;
#ifdef _MSC_VER
fopen_s( &implicitSolventResultsFile, resultsFileName.c_str(), "w" );
#else
implicitSolventResultsFile = fopen( resultsFileName.c_str(), "w" );
#endif
// diganostics
std::stringstream message;
message << methodName;
if( implicitSolventResultsFile != NULL ){
std::stringstream message;
message << methodName;
message << " Opened file=<" << resultsFileName << ">.";
//SimTKOpenMMLog::printMessage( message );
(void) fprintf( stderr, "%s", message.str().c_str() );
} else {
std::stringstream message;
message << methodName;
message << " could not open file=<" << resultsFileName << "> -- abort output.";
//SimTKOpenMMLog::printMessage( message );
(void) fprintf( stderr, "%s", message.str().c_str() );
return -1;
}
// header
(void) fprintf( implicitSolventResultsFile, "# %d atoms E=%.7e format: coords(3) bornRadii(input) q atomicRadii scaleFactors forces obcChain\n",
numberOfAtoms, energy );
RealOpenMM forceConversion = (RealOpenMM) 1.0;
RealOpenMM lengthConversion = (RealOpenMM) 1.0;
// output
if( forces != NULL && atomCoordinates != NULL && partialCharges != NULL && atomicRadii != NULL ){
for( int ii = 0; ii < numberOfAtoms; ii++ ){
(void) fprintf( implicitSolventResultsFile, "%.7e %.7e %.7e %.7e %.5f %.5f %.5f %.7e %.7e %.7e %.7e\n",
lengthConversion*atomCoordinates[ii][0],
lengthConversion*atomCoordinates[ii][1],
lengthConversion*atomCoordinates[ii][2],
(bornRadii != NULL ? lengthConversion*bornRadii[ii] : 0.0),
partialCharges[ii], lengthConversion*atomicRadii[ii], scaledRadii[ii],
forceConversion*forces[ii][0],
forceConversion*forces[ii][1],
forceConversion*forces[ii][2],
forceConversion*obcChain[ii]
);
}
}
(void) fclose( implicitSolventResultsFile );
return 0;
}
/**---------------------------------------------------------------------------------------
Write results from first loop
@param numberOfAtoms number of atoms
@param forces forces
@param bornForce Born force prefactor
@param outputFileName output file name
@return 0 unless
file cannot be opened
in which case return -1
--------------------------------------------------------------------------------------- */
int CpuObcSoftcore::writeForceLoop1( int numberOfAtoms, RealOpenMM** forces, const RealOpenMM* bornForce,
const std::string& outputFileName ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObcSoftcore::writeForceLoop1";
// ---------------------------------------------------------------------------------------
int chunkSize;
if( bornForce ){
chunkSize = 3;
} else {
chunkSize = 4;
}
StringVector lineVector;
std::stringstream header;
lineVector.push_back( "# bornF F" );
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
std::stringstream line;
line << (atomI+1) << " ";
SimTKOpenMMUtilities::formatRealStringStream( line, forces[atomI], chunkSize );
if( bornForce ){
line << " " << bornForce[atomI];
}
lineVector.push_back( line.str() );
}
return SimTKOpenMMUtilities::writeFile( lineVector, outputFileName );
}
/**---------------------------------------------------------------------------------------
Write results
@param numberOfAtoms number of atoms
@param chunkSizes vector of chunk sizes for realRealOpenMMVector
@param realRealOpenMMVector vector of RealOpenMM**
@param realVector vector of RealOpenMM*
@param outputFileName output file name
@return SimTKOpenMMCommon::DefaultReturn unless
file cannot be opened
in which case return SimTKOpenMMCommon::ErrorReturn
--------------------------------------------------------------------------------------- */
int CpuObcSoftcore::writeForceLoop( int numberOfAtoms, const IntVector& chunkSizes,
const RealOpenMMPtrPtrVector& realRealOpenMMVector,
const RealOpenMMPtrVector& realVector,
const std::string& outputFileName ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObcSoftcore::writeForceLoop";
static const int maxChunks = 10;
int chunks[maxChunks];
// ---------------------------------------------------------------------------------------
for( int ii = 0; ii < (int) chunkSizes.size(); ii++ ){
chunks[ii] = chunkSizes[ii];
}
for( int ii = (int) chunkSizes.size(); ii < maxChunks; ii++ ){
chunks[ii] = 3;
}
StringVector lineVector;
std::stringstream header;
// lineVector.push_back( "# " );
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
std::stringstream line;
char buffer[128];
(void) sprintf( buffer, "%4d ", atomI );
line << buffer;
int index = 0;
for( RealOpenMMPtrPtrVectorCI ii = realRealOpenMMVector.begin(); ii != realRealOpenMMVector.end(); ii++ ){
RealOpenMM** forces = *ii;
(void) sprintf( buffer, "%11.5f %11.5f %11.5f ", forces[atomI][0], forces[atomI][1], forces[atomI][2] );
line << buffer;
// SimTKOpenMMUtilities::formatRealStringStream( line, forces[atomI], chunks[index++] );
// line << " ";
}
for( RealOpenMMPtrVectorCI ii = realVector.begin(); ii != realVector.end(); ii++ ){
RealOpenMM* array = *ii;
(void) sprintf( buffer, "%11.5f ", array[atomI] );
line << buffer;
}
lineVector.push_back( line.str() );
}
return SimTKOpenMMUtilities::writeFile( lineVector, outputFileName );
}
/**---------------------------------------------------------------------------------------
Get Obc Born energy and forces -- used debugging
@param bornRadii Born radii -- optional; if NULL, then ObcSoftcoreParameters
entry is used
@param atomCoordinates atomic coordinates
@param partialCharges partial charges
@param forces forces
@return SimTKOpenMMCommon::DefaultReturn;
The array bornRadii is also updated and the obcEnergy
--------------------------------------------------------------------------------------- */
int CpuObcSoftcore::computeBornEnergyForcesPrint( RealOpenMM* bornRadii, RealOpenMM** atomCoordinates,
const RealOpenMM* partialCharges, RealOpenMM** forces ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObcSoftcore::computeBornEnergyForcesPrint";
static const RealOpenMM zero = (RealOpenMM) 0.0;
static const RealOpenMM one = (RealOpenMM) 1.0;
static const RealOpenMM two = (RealOpenMM) 2.0;
static const RealOpenMM three = (RealOpenMM) 3.0;
static const RealOpenMM four = (RealOpenMM) 4.0;
static const RealOpenMM half = (RealOpenMM) 0.5;
static const RealOpenMM fourth = (RealOpenMM) 0.25;
static const RealOpenMM eighth = (RealOpenMM) 0.125;
// ---------------------------------------------------------------------------------------
const ObcSoftcoreParameters* obcSoftcoreParameters = getObcSoftcoreParameters();
const int numberOfAtoms = obcSoftcoreParameters->getNumberOfAtoms();
if( bornRadii == NULL ){
bornRadii = getBornRadii();
}
// suppress warning about fopen in Visual Studio
#if defined(_MSC_VER)
#pragma warning(push)
#pragma warning(disable:4996)
#endif
FILE* logFile = NULL;
//FILE* logFile = SimTKOpenMMLog::getSimTKOpenMMLogFile( );
//FILE* logFile = fopen( "bF", "w" );
#if defined(_MSC_VER)
#pragma warning(pop)
#endif
// ---------------------------------------------------------------------------------------
// constants
const RealOpenMM preFactor = obcSoftcoreParameters->getPreFactor();
const RealOpenMM dielectricOffset = obcSoftcoreParameters->getDielectricOffset();
// ---------------------------------------------------------------------------------------
// set energy/forces to zero
RealOpenMM obcEnergy = zero;
const unsigned int arraySzInBytes = sizeof( RealOpenMM )*numberOfAtoms;
for( int ii = 0; ii < numberOfAtoms; ii++ ){
memset( forces[ii], 0, 3*sizeof( RealOpenMM ) );
}
RealOpenMM* bornForces = getBornForce();
memset( bornForces, 0, arraySzInBytes );
// ---------------------------------------------------------------------------------------
// N*( 8 + pow) ACE
// compute the nonpolar solvation via ACE approximation
if( includeAceApproximation() ){
computeAceNonPolarForce( obcSoftcoreParameters, bornRadii, &obcEnergy, bornForces );
if( logFile ){
(void) fprintf( logFile, "\nACE E=%.5e\n", obcEnergy );
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
(void) fprintf( logFile, " %d bR=%.6e bF=%.6e\n", atomI, bornRadii[atomI], bornForces[atomI] );
}
}
}
// ---------------------------------------------------------------------------------------
// first main loop
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
RealOpenMM partialChargeI = preFactor*partialCharges[atomI];
for( int atomJ = atomI; atomJ < numberOfAtoms; atomJ++ ){
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (_obcSoftcoreParameters->getPeriodic())
ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomI], atomCoordinates[atomJ], _obcSoftcoreParameters->getPeriodicBox(), deltaR );
else
ReferenceForce::getDeltaR( atomCoordinates[atomI], atomCoordinates[atomJ], deltaR );
if (_obcSoftcoreParameters->getUseCutoff() && deltaR[ReferenceForce::RIndex] > _obcSoftcoreParameters->getCutoffDistance())
continue;
RealOpenMM r2 = deltaR[ReferenceForce::R2Index];
RealOpenMM deltaX = deltaR[ReferenceForce::XIndex];
RealOpenMM deltaY = deltaR[ReferenceForce::YIndex];
RealOpenMM deltaZ = deltaR[ReferenceForce::ZIndex];
// 3 FLOP
RealOpenMM alpha2_ij = bornRadii[atomI]*bornRadii[atomJ];
RealOpenMM D_ij = r2/(four*alpha2_ij);
// exp + 2 + sqrt FLOP
RealOpenMM expTerm = EXP( -D_ij );
RealOpenMM denominator2 = r2 + alpha2_ij*expTerm;
RealOpenMM denominator = SQRT( denominator2 );
// 6 FLOP
RealOpenMM Gpol = (partialChargeI*partialCharges[atomJ])/denominator;
// dGpol/dr = -1/2*(Gpol/denominator2)*(2r - r/2*exp() )
RealOpenMM dGpol_dr = -Gpol*( one - fourth*expTerm )/denominator2;
// 5 FLOP
RealOpenMM dGpol_dalpha2_ij = -half*Gpol*expTerm*( one + D_ij )/denominator2;
// 11 FLOP
if( atomI != atomJ ){
bornForces[atomJ] += dGpol_dalpha2_ij*bornRadii[atomI];
deltaX *= dGpol_dr;
deltaY *= dGpol_dr;
deltaZ *= dGpol_dr;
forces[atomI][0] += deltaX;
forces[atomI][1] += deltaY;
forces[atomI][2] += deltaZ;
forces[atomJ][0] -= deltaX;
forces[atomJ][1] -= deltaY;
forces[atomJ][2] -= deltaZ;
} else {
Gpol *= half;
}
// 3 FLOP
obcEnergy += Gpol;
bornForces[atomI] += dGpol_dalpha2_ij*bornRadii[atomJ];
//if( logFile && (atomI == -1 || atomJ == -1) ){
// (void) fprintf( logFile, "\nWWX %d %d F[%.6e %.6e %.6e] bF=[%.6e %.6e] Gpl[%.6e %.6e %.6e] rb[%6.4f %7.4f] rs[%6.4f %7.4f] ",
// atomI, atomJ,
// forces[atomI][0], forces[atomI][1], forces[atomI][2],
// bornForces[atomI], bornForces[atomJ],
// Gpol,dGpol_dr,dGpol_dalpha2_ij,
// bornRadii[atomI],bornRadii[atomJ],atomicRadii[atomI],atomicRadii[atomJ] );
//
// (void) fprintf( logFile, "\nWWX %d %d %.1f r2=%.4f q=%.2f bF=[%.6e %.6e] Gpl[%.6e %.6e %.6e] rb[%.5f %.5f] add[%.6e %.6e] ",
// atomI, atomJ, preFactor, r2, partialCharges[atomJ],
// bornForces[atomI], bornForces[atomJ],
// Gpol,dGpol_dr,dGpol_dalpha2_ij,
// bornRadii[atomI], bornRadii[atomJ],
// dGpol_dalpha2_ij*bornRadii[atomJ], dGpol_dalpha2_ij*bornRadii[atomI] );
//}
}
}
if( logFile ){
(void) fprintf( logFile, "\nWXX bF & F E=%.8e preFactor=%.5f", obcEnergy, preFactor );
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
(void) fprintf( logFile, "\nWXX %d q=%.4f bR=%.5e bF=%.3f F[%.6e %.6e %.6e] ",
atomI, partialCharges[atomI], bornRadii[atomI], bornForces[atomI], forces[atomI][0], forces[atomI][1], forces[atomI][2] );
}
}
if( 1 ){
std::string outputFileName = "Loop1Cpu.txt";
CpuObcSoftcore::writeForceLoop1( numberOfAtoms, forces, bornForces, outputFileName );
/*
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
forces[atomI][0] = forces[atomI][1] = forces[atomI][2] = (RealOpenMM) 0.0;
}
*/
}
// ---------------------------------------------------------------------------------------
// second main loop
// initialize Born radii & ObcChain temp arrays -- contain values
// used in next iteration
RealOpenMM* bornRadiiTemp = getBornRadiiTemp();
memset( bornRadiiTemp, 0, arraySzInBytes );
RealOpenMM* obcChainTemp = getObcChainTemp();
memset( obcChainTemp, 0, arraySzInBytes );
RealOpenMM* obcChain = getObcChain();
const RealOpenMM* atomicRadii = obcSoftcoreParameters->getAtomicRadii();
const RealOpenMM alphaObc = obcSoftcoreParameters->getAlphaObc();
const RealOpenMM betaObc = obcSoftcoreParameters->getBetaObc();
const RealOpenMM gammaObc = obcSoftcoreParameters->getGammaObc();
const RealOpenMM* scaledRadiusFactor = obcSoftcoreParameters->getScaledRadiusFactors();
// compute factor that depends only on the outer loop index
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
bornForces[atomI] *= bornRadii[atomI]*bornRadii[atomI]*obcChain[atomI];
}
if( 1 ){
std::string outputFileName = "PostLoop1Cpu.txt";
IntVector chunkVector;
chunkVector.push_back( 3 );
RealOpenMMPtrPtrVector realPtrPtrVector;
realPtrPtrVector.push_back( forces );
RealOpenMMPtrVector realPtrVector;
realPtrVector.push_back( bornRadii );
realPtrVector.push_back( bornForces );
realPtrVector.push_back( obcChain );
CpuObcSoftcore::writeForceLoop( numberOfAtoms, chunkVector, realPtrPtrVector, realPtrVector, outputFileName );
}
RealOpenMM* bornSumArray = (RealOpenMM*) malloc( sizeof( RealOpenMM )*numberOfAtoms );
memset( bornSumArray, 0, sizeof( RealOpenMM )*numberOfAtoms );
/*
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
forces[atomI][0] = 0.0;
forces[atomI][1] = 0.0;
forces[atomI][2] = 0.0;
} */
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
// radius w/ dielectric offset applied
RealOpenMM radiusI = atomicRadii[atomI];
RealOpenMM offsetRadiusI = radiusI - dielectricOffset;
// used to compute Born radius for next iteration
RealOpenMM bornSum = zero;
for( int atomJ = 0; atomJ < numberOfAtoms; atomJ++ ){
if( atomJ != atomI ){
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (_obcSoftcoreParameters->getPeriodic())
ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomI], atomCoordinates[atomJ], _obcSoftcoreParameters->getPeriodicBox(), deltaR );
else
ReferenceForce::getDeltaR( atomCoordinates[atomI], atomCoordinates[atomJ], deltaR );
if (_obcSoftcoreParameters->getUseCutoff() && deltaR[ReferenceForce::RIndex] > _obcSoftcoreParameters->getCutoffDistance())
continue;
RealOpenMM deltaX = deltaR[ReferenceForce::XIndex];
RealOpenMM deltaY = deltaR[ReferenceForce::YIndex];
RealOpenMM deltaZ = deltaR[ReferenceForce::ZIndex];
RealOpenMM r = deltaR[ReferenceForce::RIndex];
// radius w/ dielectric offset applied
RealOpenMM radiusJ = atomicRadii[atomJ] - dielectricOffset;
RealOpenMM scaledRadiusJ = radiusJ*scaledRadiusFactor[atomJ];
RealOpenMM scaledRadiusJ2 = scaledRadiusJ*scaledRadiusJ;
RealOpenMM rScaledRadiusJ = r + scaledRadiusJ;
// L_ij != 1 && U_ij != 1
// dL/dr & dU/dr are zero (this can be shown analytically)
// removed from calculation
if( offsetRadiusI < rScaledRadiusJ ){
RealOpenMM l_ij = offsetRadiusI > FABS( r - scaledRadiusJ ) ? offsetRadiusI : FABS( r - scaledRadiusJ );
l_ij = one/l_ij;
RealOpenMM l_ij2 = l_ij*l_ij;
RealOpenMM u_ij = one/rScaledRadiusJ;
RealOpenMM u_ij2 = u_ij*u_ij;
RealOpenMM rInverse = one/r;
RealOpenMM r2Inverse = rInverse*rInverse;
RealOpenMM logRatio = LN( u_ij/l_ij );
RealOpenMM t3 = eighth*(one + scaledRadiusJ2*r2Inverse)*(l_ij2 - u_ij2) + fourth*logRatio*r2Inverse;
RealOpenMM de = bornForces[atomI]*t3*rInverse;
deltaX *= de;
deltaY *= de;
deltaZ *= de;
forces[atomI][0] -= deltaX;
forces[atomI][1] -= deltaY;
forces[atomI][2] -= deltaZ;
forces[atomJ][0] += deltaX;
forces[atomJ][1] += deltaY;
forces[atomJ][2] += deltaZ;
// Born radius term
RealOpenMM term = l_ij - u_ij + fourth*r*(u_ij2 - l_ij2) + (half*rInverse)*logRatio + (fourth*scaledRadiusJ*scaledRadiusJ*rInverse)*(l_ij2-u_ij2);
if( offsetRadiusI < (scaledRadiusJ - r) ){
term += two*( (one/offsetRadiusI) - l_ij);
}
bornSum += term;
if( atomI == -1 || atomJ == -1 ){
(void) fprintf( logFile, "\nXXY %d %d de=%.6e bF[%.6e %6e] t3=%.6e r=%.6e trm=%.6e bSm=%.6e f[%.6e %.6e %.6e]",
atomI, atomJ, de,
bornForces[atomI], obcChain[atomI],
t3, r, term, bornSum, forces[atomI][0], forces[atomI][1], forces[atomI][2] );
}
}
}
}
bornSumArray[atomI] = bornSum;
// OBC-specific code (Eqs. 6-8 in paper)
bornSum *= half*offsetRadiusI;
RealOpenMM sum2 = bornSum*bornSum;
RealOpenMM sum3 = bornSum*sum2;
RealOpenMM tanhSum = TANH( alphaObc*bornSum - betaObc*sum2 + gammaObc*sum3 );
bornRadiiTemp[atomI] = one/( one/offsetRadiusI - tanhSum/radiusI);
obcChainTemp[atomI] = offsetRadiusI*( alphaObc - two*betaObc*bornSum + three*gammaObc*sum2 );
obcChainTemp[atomI] = (one - tanhSum*tanhSum)*obcChainTemp[atomI]/radiusI;
if( logFile && atomI >= 0 ){
(void) fprintf( logFile, "\nXXX %d bSum[%.6e %.6e %.6e] bRt=[%.6e %6e] obc=%.6e rI=[%.5f %.5f]",
atomI, bornSumArray[atomI], bornSum, tanhSum, bornRadii[atomI], bornRadiiTemp[atomI], obcChainTemp[atomI], radiusI, offsetRadiusI );
}
}
setEnergy( obcEnergy );
if( 1 ){
std::string outputFileName = "Loop2Cpu.txt";
IntVector chunkVector;
chunkVector.push_back( 3 );
RealOpenMMPtrPtrVector realPtrPtrVector;
realPtrPtrVector.push_back( forces );
RealOpenMMPtrVector realPtrVector;
realPtrVector.push_back( bornSumArray );
// realPtrVector.push_back( bornRadiiTemp );
// realPtrVector.push_back( obcChainTemp );
CpuObcSoftcore::writeForceLoop( numberOfAtoms, chunkVector, realPtrPtrVector, realPtrVector, outputFileName );
}
if( bornSumArray ){
free( (char*) bornSumArray );
}
// 6 FLOP
/*
RealOpenMM forceFactor = getForceConversionFactor();
RealOpenMM constantFactor = 1.0f/electricConstant;
if( fabs(forceFactor - 1.0f) > 1.0e-04 ){
constantFactor *= forceFactor;
for( int ii = 0; ii < numberOfAtoms; ii++ ){
forces[ii][0] *= forceFactor;
forces[ii][1] *= forceFactor;
forces[ii][2] *= forceFactor;
}
} */
// copy new Born radii and obcChain values into permanent array
//(void) fprintf( logFile, "\nBorn radii not being updated!!!!" );
memcpy( bornRadii, bornRadiiTemp, arraySzInBytes );
memcpy( obcChain, obcChainTemp, arraySzInBytes );
if( logFile ){
(void) fclose( logFile );
}
return 0;
}
plugins/freeEnergy/platforms/reference/src/gbsa/CpuObcSoftcore.h
View file @ 0843c5f3
......@@ -128,7 +128,7 @@ class CpuObcSoftcore : public CpuImplicitSolvent {
--------------------------------------------------------------------------------------- */
int computeBornRadii( RealOpenMM** atomCoordinates, RealOpenMM* bornRadii,
int computeBornRadii( std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM* bornRadii,
RealOpenMM* obcChain = NULL );
/**---------------------------------------------------------------------------------------
......@@ -162,79 +162,8 @@ class CpuObcSoftcore : public CpuImplicitSolvent {
--------------------------------------------------------------------------------------- */
int computeBornEnergyForces( RealOpenMM* bornRadii, RealOpenMM** atomCoordinates,
const RealOpenMM* partialCharges, RealOpenMM** forces );
int computeBornEnergyForcesPrint( RealOpenMM* bornRadii, RealOpenMM** atomCoordinates,
const RealOpenMM* partialCharges, RealOpenMM** forces );
/**---------------------------------------------------------------------------------------
Get state
title title (optional)
@return state string
--------------------------------------------------------------------------------------- */
std::string getStateString( const char* title ) const;
/**---------------------------------------------------------------------------------------
Write Born energy and forces (Simbios)
@param atomCoordinates atomic coordinates
@param partialCharges partial atom charges
@param forces force array
@param resultsFileName output file name
@return SimTKOpenMMCommon::DefaultReturn if file opened; else return SimTKOpenMMCommon::ErrorReturn
--------------------------------------------------------------------------------------- */
int writeBornEnergyForces( RealOpenMM** atomCoordinates,
const RealOpenMM* partialCharges, RealOpenMM** forces,
const std::string& resultsFileName ) const;
/**---------------------------------------------------------------------------------------
Write results from first loop
@param atomCoordinates atomic coordinates
@param RealOpenMM forces forces
@param outputFileName output file name
@return SimTKOpenMMCommon::DefaultReturn unless
file cannot be opened
in which case return SimTKOpenMMCommon::ErrorReturn
--------------------------------------------------------------------------------------- */
static int writeForceLoop1( int numberOfAtoms, RealOpenMM** forces, const RealOpenMM* bornForce,
const std::string& outputFileName );
/**---------------------------------------------------------------------------------------
Write results
@param numberOfAtoms number of atoms
@param chunkSizes vector of chunk sizes for realRealOpenMMVector
@param realRealOpenMMVector vector of RealOpenMM**
@param realVector vector of RealOpenMM*
@param outputFileName output file name
@return SimTKOpenMMCommon::DefaultReturn unless
file cannot be opened
in which case return SimTKOpenMMCommon::ErrorReturn
--------------------------------------------------------------------------------------- */
static int writeForceLoop( int numberOfAtoms, const IntVector& chunkSizes,
const RealOpenMMPtrPtrVector& realRealOpenMMVector,
const RealOpenMMPtrVector& realVector,
const std::string& outputFileName );
int computeBornEnergyForces( RealOpenMM* bornRadii, std::vector<OpenMM::RealVec>& atomCoordinates,
const RealOpenMM* partialCharges, std::vector<OpenMM::RealVec>& forces );
};
// ---------------------------------------------------------------------------------------
......
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