Skip to content

GitLab

Commit 767ea1bd authored by Mark Friedrichs's avatar Mark Friedrichs
Browse files

Removed sigma4/epsilon4 from Vdw to make consistent w/ latest version of Tinker

parent 5f0e2f92
Hide whitespace changes
Inline Side-by-side
Showing with 63 additions and 306 deletions
plugins/amoeba/openmmapi/include/AmoebaVdwForce.h
View file @ 767ea1bd
......@@ -58,46 +58,6 @@ public:
return parameters.size();
}
/**
* Set table size
*
* @param tableSize table size
*/
void setSigEpsTableSize( int tableSize );
/**
* Get the SigEps table size
*
* @return the table size
*/
int getSigEpsTableSize( void ) const;
/**
* Set a entry to sigma/epsilon table.
*
* @param indexI row index of table
* @param indexJ column index of table
* @param sigma sigma entry
* @param epsilon epsilon entry
* @param sigma4 sigma4 entry
* @param epsilon4 epsilon4 entry
*/
void setSigEpsTableEntry(int indexI, int indexJ, double sigma, double epsilon,
double sigma4, double epsilon4 );
/**
* Get a entry to sigma/epsilon table.
*
* @param indexI row index of table
* @param indexJ column index of table
* @param sigma sigma entry
* @param epsilon epsilon entry
* @param sigma4 sigma4 entry
* @param epsilon4 epsilon4 entry
*/
void getSigEpsTableEntry(int indexI, int indexJ, double& sigma, double& epsilon,
double& sigma4, double& epsilon4 ) const;
/**
* Set the force field parameters for a vdw particle.
*
......@@ -105,13 +65,11 @@ public:
* @param ivIndex the iv index
* @param classIndex the class index into the sig-eps table
* @param sigma vdw sigma
* @param sigma4 1-4 vdw sigma
* @param epsilon vdw epsilon
* @param epsilon4 1-4 vdw epsilon
* @param reductionFactor the reduction factor for 1-4 ixns
* @param reductionFactor the reduction factor
*/
void setParticleParameters(int particleIndex, int ivIndex, int classIndex,
double sigma, double sigma4, double epsilon, double epsilon4, double reductionFactor );
double sigma, double epsilon, double reductionFactor );
/**
* Get the force field parameters for a vdw particle.
......@@ -120,13 +78,11 @@ public:
* @param ivIndex the iv index
* @param classIndex the class index into the sig-eps table
* @param sigma vdw sigma
* @param sigma4 1-4 vdw sigma
* @param epsilon vdw epsilon
* @param epsilon4 1-4 vdw epsilon
* @param reductionFactor the reduction factor for 1-4 ixns
* @param reductionFactor the reduction factor
*/
void getParticleParameters(int particleIndex, int& ivIndex, int& classIndex,
double& sigma, double& sigma4, double& epsilon, double& epsilon4, double& reductionFactor ) const;
double& sigma, double& epsilon, double& reductionFactor ) const;
/**
......@@ -136,14 +92,12 @@ public:
* @param ivIndex the iv index
* @param classIndex the class index into the sig-eps table
* @param sigma vdw sigma
* @param sigma4 1-4 vdw sigma
* @param epsilon vdw epsilon
* @param epsilon4 1-4 vdw epsilon
* @param reductionFactor the reduction factor for 1-4 ixns
* @param reductionFactor the reduction factor
* @return index of added particle
*/
int addParticle(int ivIndex, int classIndex,
double sigma, double sigma4, double epsilon, double epsilon4, double reductionFactor );
double sigma, double epsilon, double reductionFactor );
/**
* Set sigma combining rule
......@@ -216,17 +170,15 @@ private:
class AmoebaVdwForce::VdwInfo {
public:
int ivIndex, classIndex;
double reductionFactor, sigma, sigma4, epsilon, epsilon4;
double reductionFactor, sigma, epsilon;
VdwInfo() {
ivIndex = classIndex = -1;
reductionFactor = 0.0;
sigma = 1.0;
sigma4 = 1.0;
epsilon = 0.0;
epsilon4 = 0.0;
}
VdwInfo(int ivIndex, int classIndex, double sigma, double sigma4, double epsilon, double epsilon4, double reductionFactor ) :
ivIndex(ivIndex), classIndex(classIndex), sigma(sigma), sigma4(sigma4), epsilon(epsilon), epsilon4(epsilon4), reductionFactor(reductionFactor) {
VdwInfo(int ivIndex, int classIndex, double sigma, double epsilon, double reductionFactor ) :
ivIndex(ivIndex), classIndex(classIndex), sigma(sigma), epsilon(epsilon), reductionFactor(reductionFactor) {
}
};
......
plugins/amoeba/openmmapi/src/AmoebaVdwForce.cpp
View file @ 767ea1bd
......@@ -39,59 +39,29 @@ using namespace OpenMM;
AmoebaVdwForce::AmoebaVdwForce() {
}
int AmoebaVdwForce::addParticle(int ivIndex, int classIndex, double sigma, double sigma4, double epsilon, double epsilon4, double reductionFactor ) {
parameters.push_back(VdwInfo(ivIndex, classIndex, sigma, sigma4, epsilon, epsilon4, reductionFactor));
int AmoebaVdwForce::addParticle(int ivIndex, int classIndex, double sigma, double epsilon, double reductionFactor ) {
parameters.push_back(VdwInfo(ivIndex, classIndex, sigma, epsilon, reductionFactor));
return parameters.size()-1;
}
void AmoebaVdwForce::getParticleParameters(int particleIndex, int& ivIndex, int& classIndex,
double& sigma, double& sigma4, double& epsilon, double& epsilon4, double& reductionFactor ) const {
double& sigma, double& epsilon, double& reductionFactor ) const {
ivIndex = parameters[particleIndex].ivIndex;
classIndex = parameters[particleIndex].classIndex;
sigma = parameters[particleIndex].sigma;
sigma4 = parameters[particleIndex].sigma4;
epsilon = parameters[particleIndex].epsilon;
epsilon4 = parameters[particleIndex].epsilon4;
reductionFactor = parameters[particleIndex].reductionFactor;
}
void AmoebaVdwForce::setParticleParameters(int particleIndex, int ivIndex, int classIndex,
double sigma, double sigma4, double epsilon, double epsilon4, double reductionFactor ) {
double sigma, double epsilon, double reductionFactor ) {
parameters[particleIndex].ivIndex = ivIndex;
parameters[particleIndex].classIndex = classIndex;
parameters[particleIndex].sigma = sigma;
parameters[particleIndex].sigma4 = sigma4;
parameters[particleIndex].epsilon = epsilon;
parameters[particleIndex].epsilon4 = epsilon4;
parameters[particleIndex].reductionFactor = reductionFactor;
}
void AmoebaVdwForce::setSigEpsTableSize(int tableSize ) {
sigEpsTable.resize( tableSize );
for( unsigned int ii = 0; ii < tableSize; ii++ ){
sigEpsTable[ii].resize( tableSize );
}
}
int AmoebaVdwForce::getSigEpsTableSize(void ) const {
return static_cast<int>(sigEpsTable.size( ));
}
void AmoebaVdwForce::setSigEpsTableEntry(int indexI, int indexJ, double combinedSigma, double combinedEpsilon, double combinedSigma4, double combinedEpsilon4 ) {
sigEpsTable[indexI][indexJ].resize( 4 );
sigEpsTable[indexI][indexJ][0] = combinedSigma;
sigEpsTable[indexI][indexJ][1] = combinedEpsilon;
sigEpsTable[indexI][indexJ][2] = combinedSigma4;
sigEpsTable[indexI][indexJ][3] = combinedEpsilon4;
}
void AmoebaVdwForce::getSigEpsTableEntry(int indexI, int indexJ, double& combinedSigma, double& combinedEpsilon, double& combinedSigma4, double& combinedEpsilon4 ) const {
combinedSigma = sigEpsTable[indexI][indexJ][0];
combinedEpsilon = sigEpsTable[indexI][indexJ][1];
combinedSigma4 = sigEpsTable[indexI][indexJ][2];
combinedEpsilon4 = sigEpsTable[indexI][indexJ][3];
}
void AmoebaVdwForce::setSigmaCombiningRule(std::string& inputSigmaCombiningRule ) {
sigmaCombiningRule = inputSigmaCombiningRule;
}
......@@ -130,7 +100,6 @@ void AmoebaVdwForce::getParticleExclusions( int particleIndex, std::vector< int
}
ForceImpl* AmoebaVdwForce::createImpl() {
return new AmoebaVdwForceImpl(*this);
}
plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.cpp
View file @ 767ea1bd
......@@ -752,21 +752,20 @@ void CudaCalcAmoebaVdwForceKernel::initialize(const System& system, const Amoeba
// per-particle parameters
int numParticles = system.getNumParticles();
std::vector<int> indexIVs(numParticles);
std::vector<int> indexClasses(numParticles);
std::vector< std::vector<int> > allExclusions(numParticles);
std::vector<float> sigmas(numParticles);
std::vector<float> epsilons(numParticles);
std::vector<float> sigma4s(numParticles);
std::vector<float> epsilon4s(numParticles);
std::vector<float> reductions(numParticles);
for( int ii = 0; ii < numParticles; ii++ ){
int indexIV, indexClass;
double sigma, sigma4, epsilon, epsilon4, reduction;
double sigma, epsilon, reduction;
std::vector<int> exclusions;
force.getParticleParameters( ii, indexIV, indexClass, sigma, sigma4, epsilon, epsilon4, reduction );
force.getParticleParameters( ii, indexIV, indexClass, sigma, epsilon, reduction );
force.getParticleExclusions( ii, exclusions );
for( unsigned int jj = 0; jj < exclusions.size(); jj++ ){
allExclusions[ii].push_back( exclusions[jj] );
......@@ -776,33 +775,11 @@ void CudaCalcAmoebaVdwForceKernel::initialize(const System& system, const Amoeba
indexClasses[ii] = indexClass;
sigmas[ii] = static_cast<float>( sigma );
epsilons[ii] = static_cast<float>( epsilon );
sigma4s[ii] = static_cast<float>( sigma4 );
epsilon4s[ii] = static_cast<float>( epsilon4 );
reductions[ii] = static_cast<float>( reduction );
}
// table
unsigned int tableSize = static_cast<unsigned int>(force.getSigEpsTableSize());
std::vector< std::vector< std::vector<float> > > sigEpsTable;
sigEpsTable.resize( tableSize );
for( unsigned int ii = 0; ii < tableSize; ii++ ){
sigEpsTable[ii].resize( tableSize );
for( unsigned int jj = 0; jj < tableSize; jj++ ){
double combinedSigma, combinedEpsilon, combinedSigma4, combinedEpsilon4;
force.getSigEpsTableEntry( ii, jj, combinedSigma, combinedEpsilon, combinedSigma4, combinedEpsilon4 );
sigEpsTable[ii][jj].resize( 4 );
sigEpsTable[ii][jj][0] = static_cast<float>( combinedSigma );
sigEpsTable[ii][jj][1] = static_cast<float>( combinedEpsilon );
sigEpsTable[ii][jj][2] = static_cast<float>( combinedSigma4 );
sigEpsTable[ii][jj][3] = static_cast<float>( combinedEpsilon4 );
}
}
gpuSetAmoebaVdwParameters( data.getAmoebaGpu(), indexIVs, indexClasses, sigmas, epsilons, sigma4s, epsilon4s, reductions,
force.getSigmaCombiningRule(), force.getEpsilonCombiningRule(), sigEpsTable,
gpuSetAmoebaVdwParameters( data.getAmoebaGpu(), indexIVs, indexClasses, sigmas, epsilons, reductions,
force.getSigmaCombiningRule(), force.getEpsilonCombiningRule(),
allExclusions );
}
......
plugins/amoeba/platforms/cuda/src/kernels/AmoebaGpu.cpp
View file @ 767ea1bd
......@@ -563,7 +563,7 @@ void gpuSetAmoebaInPlaneAngleParameters(amoebaGpuContext amoebaGpu, const std::v
psAngleID2->_pSysData[i].w = gpu->pOutputBufferCounter[psAngleID1->_pSysData[i].w]++;
#undef DUMP_PARAMETERS
#define DUMP_PARAMETERS 5
#define DUMP_PARAMETERS 50000
#if (DUMP_PARAMETERS > 0 )
if( (i < DUMP_PARAMETERS || i > bond_angles - (DUMP_PARAMETERS + 1)) && amoebaGpu->log )
fprintf( amoebaGpu->log, "InPlaneAngles: %5d [%5d %5d %5d %5d] [%5d %5d %5d %5d] A=%15.7e k=%15.7e [%5d %5d %5d %5d]\n", i,
......@@ -806,7 +806,7 @@ void gpuSetAmoebaOutOfPlaneBendParameters(amoebaGpuContext amoebaGpu, const std:
amoebaGpu->amoebaSim.amoebaOutOfPlaneBendSexticK = sexticK;
#undef DUMP_PARAMETERS
#define DUMP_PARAMETERS 5
#define DUMP_PARAMETERS 50000
#if (DUMP_PARAMETERS > 0 )
if( amoebaGpu->log )
fprintf( amoebaGpu->log, "OutOfPlaneBends: global ks[%15.7e %15.7e %15.7e %15.7e]\n", cubicK, quarticK, penticK, sexticK );
......@@ -1942,12 +1942,9 @@ void gpuSetAmoebaVdwParameters( amoebaGpuContext amoebaGpu,
const std::vector<int>& indexClasses,
const std::vector<float>& sigmas,
const std::vector<float>& epsilons,
const std::vector<float>& sigma4s,
const std::vector<float>& epsilon4s,
const std::vector<float>& reductions,
const std::string& vdwSigmaCombiningRule,
const std::string& vdwEpsilonCombiningRule,
const std::vector< std::vector< std::vector<float> > >& sigEpsTable,
const std::vector< std::vector<int> >& allExclusions )
{
// ---------------------------------------------------------------------------------------
......@@ -1959,7 +1956,6 @@ void gpuSetAmoebaVdwParameters( amoebaGpuContext amoebaGpu,
gpuContext gpu = amoebaGpu->gpuContext;
amoebaGpu->paddedNumberOfAtoms = amoebaGpu->gpuContext->sim.paddedNumberOfAtoms;
unsigned int particles = sigmas.size();
amoebaGpu->useVdwTable = 0;
// set sigma combining rule flag
......@@ -1997,40 +1993,26 @@ void gpuSetAmoebaVdwParameters( amoebaGpuContext amoebaGpu,
}
}
if( amoebaGpu->useVdwTable ){
amoebaGpu->vdwTableSize = sigEpsTable.size();
amoebaGpu->psVdwTable = new CUDAStream<float2>( amoebaGpu->vdwTableSize*amoebaGpu->vdwTableSize, 1, "VdwTable");
for (unsigned int ii = 0; ii < amoebaGpu->vdwTableSize; ii++)
{
for (unsigned int jj = 0; jj < amoebaGpu->vdwTableSize; jj++)
{
amoebaGpu->psVdwTable->_pSysStream[0][ii].x = sigEpsTable[ii][jj][0];
amoebaGpu->psVdwTable->_pSysStream[0][ii].y = sigEpsTable[ii][jj][1];
}
}
amoebaGpu->psVdwTable->Upload();
} else {
if( particles < 1 ){
(void) fprintf( stderr, "%s number of particles\n", methodName );
return;
}
amoebaGpu->psVdwSigmaEpsilon = new CUDAStream<float2>(amoebaGpu->paddedNumberOfAtoms, 1, "VdwSigmaEpsilon");
for (unsigned int ii = 0; ii < particles; ii++)
{
amoebaGpu->psVdwSigmaEpsilon->_pSysStream[0][ii].x = sigmas[ii];
amoebaGpu->psVdwSigmaEpsilon->_pSysStream[0][ii].y = epsilons[ii];
}
// Dummy out extra particles data
if( particles < 1 ){
(void) fprintf( stderr, "%s number of particles\n", methodName );
return;
}
for (unsigned int ii = particles; ii < amoebaGpu->paddedNumberOfAtoms; ii++)
{
amoebaGpu->psVdwSigmaEpsilon->_pSysStream[0][ii].x = 1.0f;
amoebaGpu->psVdwSigmaEpsilon->_pSysStream[0][ii].y = 0.0f;
}
amoebaGpu->psVdwSigmaEpsilon->Upload();
}
amoebaGpu->psVdwSigmaEpsilon = new CUDAStream<float2>(amoebaGpu->paddedNumberOfAtoms, 1, "VdwSigmaEpsilon");
for (unsigned int ii = 0; ii < particles; ii++)
{
amoebaGpu->psVdwSigmaEpsilon->_pSysStream[0][ii].x = sigmas[ii];
amoebaGpu->psVdwSigmaEpsilon->_pSysStream[0][ii].y = epsilons[ii];
}
// Dummy out extra particles data
for (unsigned int ii = particles; ii < amoebaGpu->paddedNumberOfAtoms; ii++)
{
amoebaGpu->psVdwSigmaEpsilon->_pSysStream[0][ii].x = 1.0f;
amoebaGpu->psVdwSigmaEpsilon->_pSysStream[0][ii].y = 0.0f;
}
amoebaGpu->psVdwSigmaEpsilon->Upload();
std::vector< std::vector<unsigned int> > ivMapping;
std::vector< unsigned int > ivNonMapping;
......
plugins/amoeba/platforms/cuda/src/kernels/amoebaGpuTypes.h
View file @ 767ea1bd
......@@ -317,12 +317,9 @@ void gpuSetAmoebaVdwParameters( amoebaGpuContext amoebaGpu,
const std::vector<int>& indexClasses,
const std::vector<float>& sigmas,
const std::vector<float>& epsilons,
const std::vector<float>& sigma4s,
const std::vector<float>& epsilon4s,
const std::vector<float>& reductions,
const std::string& sigmaCombiningRule,
const std::string& epsilonCombiningRule,
const std::vector< std::vector< std::vector<float> > >& sigEpsTable,
const std::vector< std::vector<int> >& allExclusions );
extern "C"
void amoebaGpuBuildVdwExclusionList( amoebaGpuContext amoebaGpu, const std::vector< std::vector<int> >& exclusions );
......
plugins/amoeba/platforms/cuda/src/kernels/kCalculateAmoebaLocalForces.cu
View file @ 767ea1bd
......@@ -306,9 +306,7 @@ __launch_bounds__(G8X_LOCALFORCES_THREADS_PER_BLOCK, 1)
void kCalculateAmoebaLocalForces_kernel()
{
unsigned int pos = blockIdx.x * blockDim.x + threadIdx.x;
//Vectors* A = &sV[threadIdx.x];
float energy = 0.0f;
float energy = 0.0f;
while (pos < cAmoebaSim.amoebaBond_offset)
{
......@@ -805,15 +803,6 @@ void kCalculateAmoebaLocalForces_kernel()
int4 atom1 = cAmoebaSim.pAmoebaPiTorsionID1[pos1];
int4 atom2 = cAmoebaSim.pAmoebaPiTorsionID2[pos1];
/*
torsionParam1.x amplitude(1)
torsionParam1.y phase(1)
torsionParam1.z amplitude(2)
torsionParam1.w phase(2)
torsionParam2.x amplitude(3)
torsionParam2.y phase(3)
*/
float4 a1 = cSim.pPosq[atom1.x];
float4 a2 = cSim.pPosq[atom1.y];
float4 a3 = cSim.pPosq[atom1.z];
......
plugins/amoeba/platforms/cuda/tests/AmoebaTinkerParameterFile.cpp
View file @ 767ea1bd
......@@ -2415,42 +2415,7 @@ static int readAmoebaVdwParameters( FILE* filePtr, MapStringInt& forceMap, const
(void) fprintf( log, "Amoeba Vdw force is not being included.\n" );
}
// read sig/eps table
bool tableLoaded = 0;
int numberOfParticles;
if( tokens[0] == "AmoebaVdw14_7SigEpsTable" ){
tableLoaded = 1;
int tableSize = atoi( tokens[1].c_str() );
if( log ){
(void) fprintf( log, "%s vdwForce table size=%d\n", methodName.c_str(), tableSize);
}
vdwForce->setSigEpsTableSize( tableSize );
for( int ii = 0; ii < tableSize; ii++ ){
StringVector lineTokens;
int isNotEof = readLine( filePtr, lineTokens, lineCount, log );
if( lineTokens.size() > 2 ){
int tokenIndex = 0;
int index = atoi( lineTokens[tokenIndex++].c_str() );
for( int jj = 0; jj < tableSize; jj++ ){
double radius = atof( lineTokens[tokenIndex++].c_str() );
double epsilon = atof( lineTokens[tokenIndex++].c_str() );
double radius4 = atof( lineTokens[tokenIndex++].c_str() );
double epsilon4 = atof( lineTokens[tokenIndex++].c_str() );
vdwForce->setSigEpsTableEntry( ii, jj, radius, epsilon, radius4, epsilon4 );
}
} else {
(void) fprintf( log, "%s AmoebaVdw table tokens incomplete at line=%d\n", methodName.c_str(), *lineCount );
exit(-1);
}
}
StringVector lineTokensT;
int isNotEof = readLine( filePtr, lineTokensT, lineCount, log );
numberOfParticles = atoi( lineTokensT[1].c_str() );
} else {
numberOfParticles = atoi( tokens[1].c_str() );
}
int numberOfParticles = atoi( tokens[1].c_str() );
// read in parameters
......@@ -2467,11 +2432,9 @@ static int readAmoebaVdwParameters( FILE* filePtr, MapStringInt& forceMap, const
int indexIV = atoi( lineTokens[tokenIndex++].c_str() );
int indexClass = atoi( lineTokens[tokenIndex++].c_str() );
double sigma = atof( lineTokens[tokenIndex++].c_str() );
double sigma4 = atof( lineTokens[tokenIndex++].c_str() );
double epsilon = atof( lineTokens[tokenIndex++].c_str() );
double epsilon4 = atof( lineTokens[tokenIndex++].c_str() );
double reduction = atof( lineTokens[tokenIndex++].c_str() );
vdwForce->addParticle( indexIV, indexClass, sigma, sigma4, epsilon, epsilon4, reduction );
vdwForce->addParticle( indexIV, indexClass, sigma, epsilon, reduction );
} else {
(void) fprintf( log, "%s AmoebaVdwForce tokens incomplete at line=%d\n", methodName.c_str(), *lineCount );
exit(-1);
......@@ -2559,40 +2522,13 @@ static int readAmoebaVdwParameters( FILE* filePtr, MapStringInt& forceMap, const
// convert units to kJ-nm from kCal-Angstrom?
if( useOpenMMUnits ){
unsigned int tableSize = static_cast<unsigned int>(vdwForce->getSigEpsTableSize());
if( tableLoaded ){
/*
for( unsigned int ii = 0; ii < tableSize; ii++ ){
for( unsigned int jj = 0; jj < tableSize; jj++ ){
double sig, eps, sig4, eps4;
vdwForce->getSigEpsTableEntry( ii, jj, sig, eps, sig4, eps4);
(void) fprintf( log, "%8d %8d %10.4f %10.4f %10.4f %10.4f\n", ii, jj, sig, eps, sig4, eps4 );
// skip to end
if( jj == maxPrint && (tableSize - maxPrint) > jj ){
jj = tableSize - maxPrint - 1;
}
}
// skip to end
if( ii == maxPrint && (tableSize - maxPrint) > ii ){
ii = tableSize - maxPrint - 1;
}
}
*/
}
for( int ii = 0; ii < vdwForce->getNumParticles(); ii++ ){
int indexIV, indexClass;
double sigma, sigma4, epsilon, epsilon4, reduction;
vdwForce->getParticleParameters( ii, indexIV, indexClass, sigma, sigma4, epsilon, epsilon4, reduction );
double sigma, epsilon, reduction;
vdwForce->getParticleParameters( ii, indexIV, indexClass, sigma, epsilon, reduction );
sigma *= AngstromToNm;
sigma4 *= AngstromToNm;
epsilon *= CalToJoule;
vdwForce->setParticleParameters( ii, indexIV, indexClass, sigma, sigma4, epsilon, epsilon4, reduction );
vdwForce->setParticleParameters( ii, indexIV, indexClass, sigma, epsilon, reduction );
}
}
......@@ -2603,44 +2539,18 @@ static int readAmoebaVdwParameters( FILE* filePtr, MapStringInt& forceMap, const
//static const unsigned int maxPrint = MAX_PRINT;
static const int maxPrint = 15;
unsigned int arraySize = static_cast<unsigned int>(vdwForce->getNumParticles());
unsigned int tableSize = static_cast<unsigned int>(vdwForce->getSigEpsTableSize());
(void) fprintf( log, "%s: %u sample of AmoebaVdwForce parameters using %s units; SigEpsTableSize=%d combining rules=[sig=%s eps=%s]\n",
(void) fprintf( log, "%s: %u sample of AmoebaVdwForce parameters using %s units; combining rules=[sig=%s eps=%s]\n",
methodName.c_str(), arraySize, (useOpenMMUnits ? "OpenMM" : "Amoeba"),
vdwForce->getSigEpsTableSize(),
vdwForce->getSigmaCombiningRule().c_str(), vdwForce->getEpsilonCombiningRule().c_str() );
if( tableLoaded ){
for( unsigned int ii = 0; ii < tableSize; ii++ ){
for( unsigned int jj = 0; jj < tableSize; jj++ ){
double sig, eps, sig4, eps4;
vdwForce->getSigEpsTableEntry( ii, jj, sig, eps, sig4, eps4);
(void) fprintf( log, "%8d %8d %10.4f %10.4f %10.4f %10.4f\n", ii, jj, sig, eps, sig4, eps4 );
// skip to end
if( jj == maxPrint && (tableSize - maxPrint) > jj ){
jj = tableSize - maxPrint - 1;
}
}
// skip to end
if( ii == maxPrint && (tableSize - maxPrint) > ii ){
ii = tableSize - maxPrint - 1;
}
}
} else {
(void) fprintf( log, "SigWEps table not loaded\n" );
}
for( int ii = 0; ii < vdwForce->getNumParticles(); ii++ ){
int indexIV, indexClass;
double sigma, sigma4, epsilon, epsilon4, reduction;
double sigma, epsilon, reduction;
std::vector< int > exclusions;
vdwForce->getParticleParameters( ii, indexIV, indexClass, sigma, sigma4, epsilon, epsilon4, reduction );
vdwForce->getParticleParameters( ii, indexIV, indexClass, sigma, epsilon, reduction );
vdwForce->getParticleExclusions( ii, exclusions );
(void) fprintf( log, "%8d %8d %8d sig=[%10.4f %10.4f] eps=[ %10.4f %10.4f] redct=%10.4f ",
ii, indexIV, indexClass, sigma, sigma4, epsilon, epsilon4, reduction );
(void) fprintf( log, "%8d %8d %8d sig=%10.4f eps=%10.4f redct=%10.4f ",
ii, indexIV, indexClass, sigma, epsilon, reduction );
(void) fprintf( log, "Excl=%3u [", static_cast<unsigned int>(exclusions.size()) );
for( unsigned int jj = 0; jj < exclusions.size(); jj++ ){
......@@ -2654,30 +2564,6 @@ static int readAmoebaVdwParameters( FILE* filePtr, MapStringInt& forceMap, const
ii = arraySize - maxPrint - 1;
}
}
// check if sigmas/sgma4s and epslon/epsilon4s same
bool anyDifferentSigma = false;
bool anyDifferentEpsilon = false;
for( unsigned int ii = 0; ii < arraySize && !anyDifferentSigma && !anyDifferentEpsilon; ii++ ){
int indexIV, indexClass;
double sigma, sigma4, epsilon, epsilon4, reduction;
vdwForce->getParticleParameters( ii, indexIV, indexClass, sigma, sigma4, epsilon, epsilon4, reduction );
if( fabs( sigma - sigma4 ) > 1.0e-05 ){
anyDifferentSigma = true;
(void) fprintf( log, "sigmas and sigma4 different at entry %d [%14.7f [%14.7f]\n", ii, sigma, sigma4 );
}
if( fabs( epsilon - epsilon4 ) > 1.0e-05 ){
anyDifferentEpsilon = true;
(void) fprintf( log, "epsilons and epsilon4 different at entry %d [%14.7f [%14.7f]\n", ii, epsilon, epsilon4 );
}
}
if( anyDifferentSigma || anyDifferentEpsilon ){
(void) fprintf( log, "sigmas and sigma4s and epsilons and epsilon4 are NOT identical.\n" );
} else {
(void) fprintf( log, "sigmas and sigma4s and epsilons and epsilon4s are identical.\n" );
}
(void) fflush( log );
}
......@@ -5157,14 +5043,19 @@ void testEnergyForceByFiniteDifference( std::string parameterFileName, MapString
}
std::vector<double> energyForceDeltas;
energyForceDeltas.push_back( 1.0e-02 );
energyForceDeltas.push_back( 5.0e-03 );
energyForceDeltas.push_back( 1.0e-03 );
energyForceDeltas.push_back( 5.0e-04 );
energyForceDeltas.push_back( 1.0e-04 );
energyForceDeltas.push_back( 5.0e-05 );
energyForceDeltas.push_back( 1.0e-05 );
energyForceDeltas.push_back( 5.0e-06 );
int scanEnergyForceDeltas = 0;
if( scanEnergyForceDeltas ){
energyForceDeltas.push_back( 1.0e-02 );
energyForceDeltas.push_back( 5.0e-03 );
energyForceDeltas.push_back( 1.0e-03 );
energyForceDeltas.push_back( 5.0e-04 );
energyForceDeltas.push_back( 1.0e-04 );
energyForceDeltas.push_back( 5.0e-05 );
energyForceDeltas.push_back( 1.0e-05 );
energyForceDeltas.push_back( 5.0e-06 );
} else {
energyForceDeltas.push_back( energyForceDelta );
}
for( unsigned int kk = 0; kk < energyForceDeltas.size(); kk++ ){
energyForceDelta = energyForceDeltas[kk];
std::vector<double> relativeDifferenceStatistics;
......
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