Skip to content

GitLab

Commit f11bafd9 authored by peastman's avatar peastman
Browse files

Merge pull request #828 from peastman/cleanup 

Cleanup to AMOEBA code
parents 162d69a2 d95b90b9
Hide whitespace changes
Inline Side-by-side
Showing with 4030 additions and 4530 deletions
plugins/amoeba/openmmapi/src/AmoebaPiTorsionForce.cpp
View file @ f11bafd9
...@@ -44,7 +44,7 @@ int AmoebaPiTorsionForce::addPiTorsion(int particle1, int particle2, int particl ...@@ -44,7 +44,7 @@ int AmoebaPiTorsionForce::addPiTorsion(int particle1, int particle2, int particl
return piTorsions.size()-1; return piTorsions.size()-1;
} }
void AmoebaPiTorsionForce::getPiTorsionParameters(int index, int& particle1, int& particle2, int& particle3, int& particle4, int& particle5, int& particle6, double& k ) const { void AmoebaPiTorsionForce::getPiTorsionParameters(int index, int& particle1, int& particle2, int& particle3, int& particle4, int& particle5, int& particle6, double& k) const {
particle1 = piTorsions[index].particle1; particle1 = piTorsions[index].particle1;
particle2 = piTorsions[index].particle2; particle2 = piTorsions[index].particle2;
particle3 = piTorsions[index].particle3; particle3 = piTorsions[index].particle3;
...@@ -54,7 +54,7 @@ void AmoebaPiTorsionForce::getPiTorsionParameters(int index, int& particle1, int ...@@ -54,7 +54,7 @@ void AmoebaPiTorsionForce::getPiTorsionParameters(int index, int& particle1, int
k = piTorsions[index].k; k = piTorsions[index].k;
} }
void AmoebaPiTorsionForce::setPiTorsionParameters(int index, int particle1, int particle2, int particle3, int particle4, int particle5, int particle6, double k ) { void AmoebaPiTorsionForce::setPiTorsionParameters(int index, int particle1, int particle2, int particle3, int particle4, int particle5, int particle6, double k) {
piTorsions[index].particle1 = particle1; piTorsions[index].particle1 = particle1;
piTorsions[index].particle2 = particle2; piTorsions[index].particle2 = particle2;
piTorsions[index].particle3 = particle3; piTorsions[index].particle3 = particle3;
......
plugins/amoeba/openmmapi/src/AmoebaStretchBendForce.cpp
View file @ f11bafd9
...@@ -46,7 +46,7 @@ int AmoebaStretchBendForce::addStretchBend(int particle1, int particle2, int par ...@@ -46,7 +46,7 @@ int AmoebaStretchBendForce::addStretchBend(int particle1, int particle2, int par
} }
void AmoebaStretchBendForce::getStretchBendParameters(int index, int& particle1, int& particle2, int& particle3, void AmoebaStretchBendForce::getStretchBendParameters(int index, int& particle1, int& particle2, int& particle3,
double& lengthAB, double& lengthCB, double& angle, double& k1, double& k2 ) const { double& lengthAB, double& lengthCB, double& angle, double& k1, double& k2) const {
particle1 = stretchBends[index].particle1; particle1 = stretchBends[index].particle1;
particle2 = stretchBends[index].particle2; particle2 = stretchBends[index].particle2;
particle3 = stretchBends[index].particle3; particle3 = stretchBends[index].particle3;
......
plugins/amoeba/openmmapi/src/AmoebaTorsionTorsionForce.cpp
View file @ f11bafd9
...@@ -70,13 +70,13 @@ void AmoebaTorsionTorsionForce::setTorsionTorsionParameters(int index, int parti ...@@ -70,13 +70,13 @@ void AmoebaTorsionTorsionForce::setTorsionTorsionParameters(int index, int parti
torsionTorsions[index].gridIndex = gridIndex; torsionTorsions[index].gridIndex = gridIndex;
} }
const TorsionTorsionGrid& AmoebaTorsionTorsionForce::getTorsionTorsionGrid(int index ) const { const TorsionTorsionGrid& AmoebaTorsionTorsionForce::getTorsionTorsionGrid(int index) const {
return torsionTorsionGrids[index].getTorsionTorsionGrid(); return torsionTorsionGrids[index].getTorsionTorsionGrid();
} }
void AmoebaTorsionTorsionForce::setTorsionTorsionGrid(int index, const TorsionTorsionGrid& grid ) { void AmoebaTorsionTorsionForce::setTorsionTorsionGrid(int index, const TorsionTorsionGrid& grid) {
if( index >= static_cast<int>(torsionTorsionGrids.size()) ){ if (index >= static_cast<int>(torsionTorsionGrids.size())) {
torsionTorsionGrids.resize( index + 1); torsionTorsionGrids.resize(index + 1);
} }
torsionTorsionGrids[index] = grid; torsionTorsionGrids[index] = grid;
} }
......
plugins/amoeba/openmmapi/src/AmoebaTorsionTorsionForceImpl.cpp
View file @ f11bafd9
...@@ -70,12 +70,11 @@ typedef std::map< double, Map_Double_IntPair > Map_Double_MapDoubleIntPair; ...@@ -70,12 +70,11 @@ typedef std::map< double, Map_Double_IntPair > Map_Double_MapDoubleIntPair;
typedef Map_Double_MapDoubleIntPair::iterator Map_Double_MapDoubleIntPairI; typedef Map_Double_MapDoubleIntPair::iterator Map_Double_MapDoubleIntPairI;
typedef Map_Double_MapDoubleIntPair::const_iterator Map_Double_MapDoubleIntPairCI; typedef Map_Double_MapDoubleIntPair::const_iterator Map_Double_MapDoubleIntPairCI;
void AmoebaTorsionTorsionForceImpl::reorderGrid( const TorsionTorsionGrid& grid, TorsionTorsionGrid& reorderedGrid ){ void AmoebaTorsionTorsionForceImpl::reorderGrid(const TorsionTorsionGrid& grid, TorsionTorsionGrid& reorderedGrid) {
reorderedGrid.resize( grid.size() ); reorderedGrid.resize(grid.size());
std::vector<Map_Double_IntPair> map_Double_IntPair_Vector( grid.size() ); std::vector<Map_Double_IntPair> map_Double_IntPair_Vector(grid.size());
Map_Double_MapDoubleIntPair mapAngles; Map_Double_MapDoubleIntPair mapAngles;
//(void) fprintf( stderr, "AmoebaTorsionTorsionForceImpl::reorder grid\n" );
// (1) set dimensions for reorderd grid // (1) set dimensions for reorderd grid
// (2) build map: // (2) build map:
...@@ -84,27 +83,27 @@ void AmoebaTorsionTorsionForceImpl::reorderGrid( const TorsionTorsionGrid& grid, ...@@ -84,27 +83,27 @@ void AmoebaTorsionTorsionForceImpl::reorderGrid( const TorsionTorsionGrid& grid,
for (unsigned int ii = 0; ii < grid.size(); ii++) { for (unsigned int ii = 0; ii < grid.size(); ii++) {
reorderedGrid[ii].resize( grid[ii].size() ); reorderedGrid[ii].resize(grid[ii].size());
for (unsigned int jj = 0; jj < grid[ii].size(); jj++) { for (unsigned int jj = 0; jj < grid[ii].size(); jj++) {
reorderedGrid[ii][jj].resize( grid[ii][jj].size() ); reorderedGrid[ii][jj].resize(grid[ii][jj].size());
double angleX = grid[ii][jj][0]; double angleX = grid[ii][jj][0];
double angleY = grid[ii][jj][1]; double angleY = grid[ii][jj][1];
if( mapAngles.find( angleX ) == mapAngles.end() ){ if (mapAngles.find(angleX) == mapAngles.end()) {
if( map_Double_IntPair_Vector[ii].size() > 0 ){ if (map_Double_IntPair_Vector[ii].size() > 0) {
char buffer[1024]; char buffer[1024];
(void) sprintf( buffer, "TorsionTorsion grid reorder: x-angle not set correctly: x=%15.7e y=%15.7e size=%u should be zero; ii/jj indies=%u %u.\n", (void) sprintf(buffer, "TorsionTorsion grid reorder: x-angle not set correctly: x=%15.7e y=%15.7e size=%u should be zero; ii/jj indies=%u %u.\n",
angleX, angleY, static_cast<unsigned int>(map_Double_IntPair_Vector[ii].size()), ii, jj ); angleX, angleY, static_cast<unsigned int>(map_Double_IntPair_Vector[ii].size()), ii, jj);
throw OpenMMException(buffer); throw OpenMMException(buffer);
} }
mapAngles[angleX] = map_Double_IntPair_Vector[ii]; mapAngles[angleX] = map_Double_IntPair_Vector[ii];
} }
Map_Double_IntPair& map_Double_IntPair = mapAngles[angleX]; Map_Double_IntPair& map_Double_IntPair = mapAngles[angleX];
if( map_Double_IntPair.find( angleY ) != map_Double_IntPair.end() ){ if (map_Double_IntPair.find(angleY) != map_Double_IntPair.end()) {
char buffer[1024]; char buffer[1024];
(void) sprintf( buffer, "TorsionTorsion grid reorder: angle pair found twice: %15.7e %15.7e %u\n", angleX, angleY, static_cast<unsigned int>(map_Double_IntPair.size()) ); (void) sprintf(buffer, "TorsionTorsion grid reorder: angle pair found twice: %15.7e %15.7e %u\n", angleX, angleY, static_cast<unsigned int>(map_Double_IntPair.size()));
throw OpenMMException(buffer); throw OpenMMException(buffer);
} }
struct IntPair pair; struct IntPair pair;
...@@ -114,24 +113,6 @@ void AmoebaTorsionTorsionForceImpl::reorderGrid( const TorsionTorsionGrid& grid, ...@@ -114,24 +113,6 @@ void AmoebaTorsionTorsionForceImpl::reorderGrid( const TorsionTorsionGrid& grid,
} }
} }
#if 0
(void) fprintf( stderr, "TorsionTorsion grid reorder map\n" );
for( Map_Double_MapDoubleIntPairCI ii = mapAngles.begin(); ii != mapAngles.end(); ii++ ){
double angleX = ii->first;
Map_Double_IntPair map_Double_IntPair = ii->second;
(void) fprintf( stderr, " %15.7e %u \n", angleX, static_cast<unsigned int>(map_Double_IntPair.size()) );
}
for( Map_Double_MapDoubleIntPairCI ii = mapAngles.begin(); ii != mapAngles.end(); ii++ ){
double angleX = ii->first;
Map_Double_IntPair map_Double_IntPair = ii->second;
for( Map_Double_IntPairCI jj = map_Double_IntPair.begin(); jj != map_Double_IntPair.end(); jj++ ){
double angle = jj->first;
struct IntPair pair = jj->second;
(void) fprintf( stderr, " %15.7e %15.7e %d %d\n", angleX, angle, pair.index1, pair.index2 );
}
}
#endif
// load reordered entries // load reordered entries
Map_Double_MapDoubleIntPairCI mapII = mapAngles.begin(); Map_Double_MapDoubleIntPairCI mapII = mapAngles.begin();
...@@ -144,7 +125,6 @@ void AmoebaTorsionTorsionForceImpl::reorderGrid( const TorsionTorsionGrid& grid, ...@@ -144,7 +125,6 @@ void AmoebaTorsionTorsionForceImpl::reorderGrid( const TorsionTorsionGrid& grid,
struct IntPair pair = mapJJ->second; struct IntPair pair = mapJJ->second;
int index1 = pair.index1; int index1 = pair.index1;
int index2 = pair.index2; int index2 = pair.index2;
//(void) fprintf( stderr, " %3d %3d %15.7e %15.7e %3d %3d zzz\n", ii, jj, mapII->first, mapJJ->first, index1, index2 );
for (unsigned int kk = 0; kk < grid[ii][jj].size(); kk++) { for (unsigned int kk = 0; kk < grid[ii][jj].size(); kk++) {
reorderedGrid[ii][jj][kk] = static_cast<float>(grid[index1][index2][kk]); reorderedGrid[ii][jj][kk] = static_cast<float>(grid[index1][index2][kk]);
...@@ -153,12 +133,12 @@ void AmoebaTorsionTorsionForceImpl::reorderGrid( const TorsionTorsionGrid& grid, ...@@ -153,12 +133,12 @@ void AmoebaTorsionTorsionForceImpl::reorderGrid( const TorsionTorsionGrid& grid,
// increment map iterators // increment map iterators
mapJJ++; mapJJ++;
if( mapJJ == map_Double_IntPair.end() ){ if (mapJJ == map_Double_IntPair.end()) {
mapII++; mapII++;
if( mapII == mapAngles.end() ){ if (mapII == mapAngles.end()) {
if( (jj != (grid[ii].size()-1)) && (ii != (grid.size()-1)) ){ if ((jj != (grid[ii].size()-1)) && (ii != (grid.size()-1))) {
char buffer[1024]; char buffer[1024];
(void) sprintf( buffer, "AmoebaTorsionTorsionForceImpl::reorderGrid: error detected with map iterators.\n" ); (void) sprintf(buffer, "AmoebaTorsionTorsionForceImpl::reorderGrid: error detected with map iterators.\n");
throw OpenMMException(buffer); throw OpenMMException(buffer);
} }
} else { } else {
......
plugins/amoeba/openmmapi/src/AmoebaVdwForce.cpp
View file @ f11bafd9
...@@ -41,13 +41,13 @@ using std::vector; ...@@ -41,13 +41,13 @@ using std::vector;
AmoebaVdwForce::AmoebaVdwForce() : nonbondedMethod(NoCutoff), sigmaCombiningRule("CUBIC-MEAN"), epsilonCombiningRule("HHG"), cutoff(1.0e+10), useDispersionCorrection(true) { AmoebaVdwForce::AmoebaVdwForce() : nonbondedMethod(NoCutoff), sigmaCombiningRule("CUBIC-MEAN"), epsilonCombiningRule("HHG"), cutoff(1.0e+10), useDispersionCorrection(true) {
} }
int AmoebaVdwForce::addParticle(int parentIndex, double sigma, double epsilon, double reductionFactor ) { int AmoebaVdwForce::addParticle(int parentIndex, double sigma, double epsilon, double reductionFactor) {
parameters.push_back(VdwInfo(parentIndex, sigma, epsilon, reductionFactor)); parameters.push_back(VdwInfo(parentIndex, sigma, epsilon, reductionFactor));
return parameters.size()-1; return parameters.size()-1;
} }
void AmoebaVdwForce::getParticleParameters(int particleIndex, int& parentIndex, void AmoebaVdwForce::getParticleParameters(int particleIndex, int& parentIndex,
double& sigma, double& epsilon, double& reductionFactor ) const { double& sigma, double& epsilon, double& reductionFactor) const {
parentIndex = parameters[particleIndex].parentIndex; parentIndex = parameters[particleIndex].parentIndex;
sigma = parameters[particleIndex].sigma; sigma = parameters[particleIndex].sigma;
epsilon = parameters[particleIndex].epsilon; epsilon = parameters[particleIndex].epsilon;
...@@ -55,14 +55,14 @@ void AmoebaVdwForce::getParticleParameters(int particleIndex, int& parentIndex, ...@@ -55,14 +55,14 @@ void AmoebaVdwForce::getParticleParameters(int particleIndex, int& parentIndex,
} }
void AmoebaVdwForce::setParticleParameters(int particleIndex, int parentIndex, void AmoebaVdwForce::setParticleParameters(int particleIndex, int parentIndex,
double sigma, double epsilon, double reductionFactor ) { double sigma, double epsilon, double reductionFactor) {
parameters[particleIndex].parentIndex = parentIndex; parameters[particleIndex].parentIndex = parentIndex;
parameters[particleIndex].sigma = sigma; parameters[particleIndex].sigma = sigma;
parameters[particleIndex].epsilon = epsilon; parameters[particleIndex].epsilon = epsilon;
parameters[particleIndex].reductionFactor = reductionFactor; parameters[particleIndex].reductionFactor = reductionFactor;
} }
void AmoebaVdwForce::setSigmaCombiningRule( const std::string& inputSigmaCombiningRule ) { void AmoebaVdwForce::setSigmaCombiningRule(const std::string& inputSigmaCombiningRule) {
sigmaCombiningRule = inputSigmaCombiningRule; sigmaCombiningRule = inputSigmaCombiningRule;
} }
...@@ -70,7 +70,7 @@ const std::string& AmoebaVdwForce::getSigmaCombiningRule() const { ...@@ -70,7 +70,7 @@ const std::string& AmoebaVdwForce::getSigmaCombiningRule() const {
return sigmaCombiningRule; return sigmaCombiningRule;
} }
void AmoebaVdwForce::setEpsilonCombiningRule( const std::string& inputEpsilonCombiningRule ) { void AmoebaVdwForce::setEpsilonCombiningRule(const std::string& inputEpsilonCombiningRule) {
epsilonCombiningRule = inputEpsilonCombiningRule; epsilonCombiningRule = inputEpsilonCombiningRule;
} }
...@@ -78,31 +78,31 @@ const std::string& AmoebaVdwForce::getEpsilonCombiningRule() const { ...@@ -78,31 +78,31 @@ const std::string& AmoebaVdwForce::getEpsilonCombiningRule() const {
return epsilonCombiningRule; return epsilonCombiningRule;
} }
void AmoebaVdwForce::setParticleExclusions( int particleIndex, const std::vector< int >& inputExclusions ) { void AmoebaVdwForce::setParticleExclusions(int particleIndex, const std::vector< int >& inputExclusions) {
if( exclusions.size() < parameters.size() ){ if (exclusions.size() < parameters.size()) {
exclusions.resize( parameters.size() ); exclusions.resize(parameters.size());
} }
if( static_cast<int>(exclusions.size()) < particleIndex ){ if (static_cast<int>(exclusions.size()) < particleIndex) {
exclusions.resize( particleIndex + 10 ); exclusions.resize(particleIndex + 10);
} }
for( unsigned int ii = 0; ii < inputExclusions.size(); ii++ ){ for (unsigned int ii = 0; ii < inputExclusions.size(); ii++) {
exclusions[particleIndex].push_back( inputExclusions[ii] ); exclusions[particleIndex].push_back(inputExclusions[ii]);
} }
} }
void AmoebaVdwForce::getParticleExclusions( int particleIndex, std::vector< int >& outputExclusions ) const { void AmoebaVdwForce::getParticleExclusions(int particleIndex, std::vector< int >& outputExclusions) const {
if( particleIndex < static_cast<int>(exclusions.size()) ){ if (particleIndex < static_cast<int>(exclusions.size())) {
outputExclusions.resize( exclusions[particleIndex].size() ); outputExclusions.resize(exclusions[particleIndex].size());
for( unsigned int ii = 0; ii < exclusions[particleIndex].size(); ii++ ){ for (unsigned int ii = 0; ii < exclusions[particleIndex].size(); ii++) {
outputExclusions[ii] = exclusions[particleIndex][ii]; outputExclusions[ii] = exclusions[particleIndex][ii];
} }
} }
} }
void AmoebaVdwForce::setCutoff( double inputCutoff ){ void AmoebaVdwForce::setCutoff(double inputCutoff) {
cutoff = inputCutoff; cutoff = inputCutoff;
} }
......
plugins/amoeba/openmmapi/src/AmoebaWcaDispersionForce.cpp
View file @ f11bafd9
...@@ -48,17 +48,17 @@ AmoebaWcaDispersionForce::AmoebaWcaDispersionForce() { ...@@ -48,17 +48,17 @@ AmoebaWcaDispersionForce::AmoebaWcaDispersionForce() {
dispoff = 0.26; dispoff = 0.26;
} }
int AmoebaWcaDispersionForce::addParticle( double radius, double epsilon ) { int AmoebaWcaDispersionForce::addParticle(double radius, double epsilon) {
parameters.push_back(WcaDispersionInfo( radius, epsilon)); parameters.push_back(WcaDispersionInfo(radius, epsilon));
return parameters.size()-1; return parameters.size()-1;
} }
void AmoebaWcaDispersionForce::getParticleParameters(int particleIndex, double& radius, double& epsilon ) const { void AmoebaWcaDispersionForce::getParticleParameters(int particleIndex, double& radius, double& epsilon) const {
radius = parameters[particleIndex].radius; radius = parameters[particleIndex].radius;
epsilon = parameters[particleIndex].epsilon; epsilon = parameters[particleIndex].epsilon;
} }
void AmoebaWcaDispersionForce::setParticleParameters(int particleIndex, double radius, double epsilon ) { void AmoebaWcaDispersionForce::setParticleParameters(int particleIndex, double radius, double epsilon) {
parameters[particleIndex].radius = radius; parameters[particleIndex].radius = radius;
parameters[particleIndex].epsilon = epsilon; parameters[particleIndex].epsilon = epsilon;
} }
...@@ -95,35 +95,35 @@ double AmoebaWcaDispersionForce::getSlevy() const { ...@@ -95,35 +95,35 @@ double AmoebaWcaDispersionForce::getSlevy() const {
return slevy; return slevy;
} }
void AmoebaWcaDispersionForce::setEpso( double inputEpso ){ void AmoebaWcaDispersionForce::setEpso(double inputEpso) {
epso = inputEpso; epso = inputEpso;
} }
void AmoebaWcaDispersionForce::setEpsh( double inputEpsh ){ void AmoebaWcaDispersionForce::setEpsh(double inputEpsh) {
epsh = inputEpsh; epsh = inputEpsh;
} }
void AmoebaWcaDispersionForce::setRmino( double inputRmino ){ void AmoebaWcaDispersionForce::setRmino(double inputRmino) {
rmino = inputRmino; rmino = inputRmino;
} }
void AmoebaWcaDispersionForce::setRminh( double inputRminh ){ void AmoebaWcaDispersionForce::setRminh(double inputRminh) {
rminh = inputRminh; rminh = inputRminh;
} }
void AmoebaWcaDispersionForce::setAwater( double inputAwater ){ void AmoebaWcaDispersionForce::setAwater(double inputAwater) {
awater = inputAwater; awater = inputAwater;
} }
void AmoebaWcaDispersionForce::setShctd( double inputShctd ){ void AmoebaWcaDispersionForce::setShctd(double inputShctd) {
shctd = inputShctd; shctd = inputShctd;
} }
void AmoebaWcaDispersionForce::setDispoff( double inputDispoff ){ void AmoebaWcaDispersionForce::setDispoff(double inputDispoff) {
dispoff = inputDispoff; dispoff = inputDispoff;
} }
void AmoebaWcaDispersionForce::setSlevy( double inputSlevy ){ void AmoebaWcaDispersionForce::setSlevy(double inputSlevy) {
slevy = inputSlevy; slevy = inputSlevy;
} }
......
plugins/amoeba/openmmapi/src/AmoebaWcaDispersionForceImpl.cpp
View file @ f11bafd9
...@@ -61,15 +61,15 @@ double AmoebaWcaDispersionForceImpl::calcForcesAndEnergy(ContextImpl& context, b ...@@ -61,15 +61,15 @@ double AmoebaWcaDispersionForceImpl::calcForcesAndEnergy(ContextImpl& context, b
return kernel.getAs<CalcAmoebaWcaDispersionForceKernel>().execute(context, includeForces, includeEnergy); return kernel.getAs<CalcAmoebaWcaDispersionForceKernel>().execute(context, includeForces, includeEnergy);
return 0.0; return 0.0;
} }
void AmoebaWcaDispersionForceImpl::getMaximumDispersionEnergy( const AmoebaWcaDispersionForce& force, int particleIndex, double& maxDispersionEnergy ) { void AmoebaWcaDispersionForceImpl::getMaximumDispersionEnergy(const AmoebaWcaDispersionForce& force, int particleIndex, double& maxDispersionEnergy) {
const double pi = 3.1415926535897; const double pi = 3.1415926535897;
// from last loop in subroutine knp in ksolv.f // from last loop in subroutine knp in ksolv.f
double rdisp, epsi; double rdisp, epsi;
force.getParticleParameters( particleIndex, rdisp, epsi ); force.getParticleParameters(particleIndex, rdisp, epsi);
if( epsi <= 0.0 || rdisp <= 0.0 ){ if (epsi <= 0.0 || rdisp <= 0.0) {
maxDispersionEnergy = 0.0; maxDispersionEnergy = 0.0;
return; return;
} }
...@@ -104,32 +104,27 @@ void AmoebaWcaDispersionForceImpl::getMaximumDispersionEnergy( const AmoebaWcaDi ...@@ -104,32 +104,27 @@ void AmoebaWcaDispersionForceImpl::getMaximumDispersionEnergy( const AmoebaWcaDi
double rdisp11 = rdisp7*rdisp3*rdisp; double rdisp11 = rdisp7*rdisp3*rdisp;
double cdisp; double cdisp;
if( rdisp < rmixh) { if (rdisp < rmixh) {
cdisp = -4.0*pi*emixh*(rmixh3-rdisp3)/3.0 - emixh*18.0/11.0*rmixh3*pi; cdisp = -4.0*pi*emixh*(rmixh3-rdisp3)/3.0 - emixh*18.0/11.0*rmixh3*pi;
} else { } else {
cdisp = 2.0*pi*(2.0*rmixh7-11.0*rdisp7)*ah/ (11.0*rdisp11); cdisp = 2.0*pi*(2.0*rmixh7-11.0*rdisp7)*ah/ (11.0*rdisp11);
} }
cdisp *= 2.0; cdisp *= 2.0;
if (rdisp < rmixo ) { if (rdisp < rmixo) {
cdisp -= 4.0*pi*emixo*(rmixo3-rdisp3)/3.0; cdisp -= 4.0*pi*emixo*(rmixo3-rdisp3)/3.0;
cdisp -= emixo*18.0/11.0*rmixo3*pi; cdisp -= emixo*18.0/11.0*rmixo3*pi;
} else { } else {
cdisp += 2.0*pi*(2.0*rmixo7-11.0*rdisp7) * ao/(11.0*rdisp11); cdisp += 2.0*pi*(2.0*rmixo7-11.0*rdisp7) * ao/(11.0*rdisp11);
} }
maxDispersionEnergy = force.getSlevy()*force.getAwater()*cdisp; maxDispersionEnergy = force.getSlevy()*force.getAwater()*cdisp;
// (void) fprintf( stderr,"Wca %5d %14.7e %14.7e %14.7e %14.7e %14.7e %14.7e %14.7e %14.7e\n",
// particleIndex, rdisp,rmini,epsi, emixh,rmixh,emixo,rmixo,cdisp );
return;
} }
double AmoebaWcaDispersionForceImpl::getTotalMaximumDispersionEnergy( const AmoebaWcaDispersionForce& force ){ double AmoebaWcaDispersionForceImpl::getTotalMaximumDispersionEnergy(const AmoebaWcaDispersionForce& force) {
double totalMaximumDispersionEnergy = 0.0; double totalMaximumDispersionEnergy = 0.0;
for( int ii = 0; ii < force.getNumParticles(); ii++ ){ for (int ii = 0; ii < force.getNumParticles(); ii++) {
double maximumDispersionEnergy; double maximumDispersionEnergy;
getMaximumDispersionEnergy( force, ii, maximumDispersionEnergy ); getMaximumDispersionEnergy(force, ii, maximumDispersionEnergy);
totalMaximumDispersionEnergy += maximumDispersionEnergy; totalMaximumDispersionEnergy += maximumDispersionEnergy;
} }
return totalMaximumDispersionEnergy; return totalMaximumDispersionEnergy;
......
plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.cpp
View file @ f11bafd9
...@@ -774,16 +774,16 @@ public: ...@@ -774,16 +774,16 @@ public:
vector<double> dipole1, dipole2, quadrupole1, quadrupole2; vector<double> dipole1, dipole2, quadrupole1, quadrupole2;
force.getMultipoleParameters(particle1, charge1, dipole1, quadrupole1, axis1, multipole11, multipole21, multipole31, thole1, damping1, polarity1); force.getMultipoleParameters(particle1, charge1, dipole1, quadrupole1, axis1, multipole11, multipole21, multipole31, thole1, damping1, polarity1);
force.getMultipoleParameters(particle2, charge2, dipole2, quadrupole2, axis2, multipole12, multipole22, multipole32, thole2, damping2, polarity2); force.getMultipoleParameters(particle2, charge2, dipole2, quadrupole2, axis2, multipole12, multipole22, multipole32, thole2, damping2, polarity2);
if (charge1 != charge2 || thole1 != thole2 || damping1 != damping2 || polarity1 != polarity2 || axis1 != axis2){ if (charge1 != charge2 || thole1 != thole2 || damping1 != damping2 || polarity1 != polarity2 || axis1 != axis2) {
return false; return false;
} }
for (int i = 0; i < (int) dipole1.size(); ++i){ for (int i = 0; i < (int) dipole1.size(); ++i) {
if (dipole1[i] != dipole2[i]){ if (dipole1[i] != dipole2[i]) {
return false; return false;
} }
} }
for (int i = 0; i < (int) quadrupole1.size(); ++i){ for (int i = 0; i < (int) quadrupole1.size(); ++i) {
if (quadrupole1[i] != quadrupole2[i]){ if (quadrupole1[i] != quadrupole2[i]) {
return false; return false;
} }
} }
......
plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.h
View file @ f11bafd9
...@@ -343,7 +343,7 @@ public: ...@@ -343,7 +343,7 @@ public:
* @param outputElectrostaticPotential output potential * @param outputElectrostaticPotential output potential
*/ */
void getElectrostaticPotential(ContextImpl& context, const std::vector< Vec3 >& inputGrid, void getElectrostaticPotential(ContextImpl& context, const std::vector< Vec3 >& inputGrid,
std::vector< double >& outputElectrostaticPotential ); std::vector< double >& outputElectrostaticPotential);
/** /**
* Get the system multipole moments * Get the system multipole moments
...@@ -353,7 +353,7 @@ public: ...@@ -353,7 +353,7 @@ public:
* dipole_x, dipole_y, dipole_z, * dipole_x, dipole_y, dipole_z,
* quadrupole_xx, quadrupole_xy, quadrupole_xz, * quadrupole_xx, quadrupole_xy, quadrupole_xz,
* quadrupole_yx, quadrupole_yy, quadrupole_yz, * quadrupole_yx, quadrupole_yy, quadrupole_yz,
* quadrupole_zx, quadrupole_zy, quadrupole_zz ) * quadrupole_zx, quadrupole_zy, quadrupole_zz)
*/ */
void getSystemMultipoleMoments(ContextImpl& context, std::vector<double>& outputMultipoleMoments); void getSystemMultipoleMoments(ContextImpl& context, std::vector<double>& outputMultipoleMoments);
/** /**
......
plugins/amoeba/platforms/cuda/src/kernels/bicubic.cu
View file @ f11bafd9
...@@ -26,13 +26,13 @@ __device__ void bicubic(real4 y, real4 y1i, real4 y2i, real4 y12i, real x1, real ...@@ -26,13 +26,13 @@ __device__ void bicubic(real4 y, real4 y1i, real4 y2i, real4 y12i, real x1, real
-2.0f*y12.x - y12.y - y12.z - 2.0f*y12.w; -2.0f*y12.x - y12.y - y12.z - 2.0f*y12.w;
c[3][0] = 2.0f*(y.x - y.y) + y1.x + y1.y; c[3][0] = 2.0f*(y.x - y.y) + y1.x + y1.y;
c[3][1] = 2.0f*(y2.x - y2.y) + y12.x + y12.y; c[3][1] = 2.0f*(y2.x - y2.y) + y12.x + y12.y;
c[3][2] = 6.0f*( y.y - y.x + y.w - y.z) + c[3][2] = 6.0f*(y.y - y.x + y.w - y.z) +
3.0f*(y1.z + y1.w - y1.x - y1.y) + 3.0f*(y1.z + y1.w - y1.x - y1.y) +
2.0f*( 2.0f*(y2.y - y2.x) + y2.z - y2.w) + 2.0f*(2.0f*(y2.y - y2.x) + y2.z - y2.w) +
-2.0f*(y12.x + y12.y) - y12.z - y12.w; -2.0f*(y12.x + y12.y) - y12.z - y12.w;
c[3][3] = 4.0f*( y.x - y.y + y.z - y.w) + c[3][3] = 4.0f*( y.x - y.y + y.z - y.w) +
2.0f*( y1.x + y1.y - y1.z - y1.w) + 2.0f*(y1.x + y1.y - y1.z - y1.w) +
2.0f*( y2.x - y2.y - y2.z + y2.w) + 2.0f*(y2.x - y2.y - y2.z + y2.w) +
y12.x + y12.y + y12.z + y12.w; y12.x + y12.y + y12.z + y12.w;
real t = (x1-x1l) / (x1u-x1l); real t = (x1-x1l) / (x1u-x1l);
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaAngleForce.cpp
View file @ f11bafd9
...@@ -66,7 +66,7 @@ const double TOL = 1e-5; ...@@ -66,7 +66,7 @@ const double TOL = 1e-5;
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
static void crossProductVector3( double* vectorX, double* vectorY, double* vectorZ ){ static void crossProductVector3(double* vectorX, double* vectorY, double* vectorZ) {
vectorZ[0] = vectorX[1]*vectorY[2] - vectorX[2]*vectorY[1]; vectorZ[0] = vectorX[1]*vectorY[2] - vectorX[2]*vectorY[1];
vectorZ[1] = vectorX[2]*vectorY[0] - vectorX[0]*vectorY[2]; vectorZ[1] = vectorX[2]*vectorY[0] - vectorX[0]*vectorY[2];
...@@ -75,26 +75,19 @@ static void crossProductVector3( double* vectorX, double* vectorY, double* vecto ...@@ -75,26 +75,19 @@ static void crossProductVector3( double* vectorX, double* vectorY, double* vecto
return; return;
} }
static void getPrefactorsGivenAngleCosine( double cosine, double idealAngle, double quadraticK, double cubicK, static void getPrefactorsGivenAngleCosine(double cosine, double idealAngle, double quadraticK, double cubicK,
double quarticK, double penticK, double sexticK, double quarticK, double penticK, double sexticK,
double* dEdR, double* energyTerm, FILE* log ) { double* dEdR, double* energyTerm) {
double angle; double angle;
if( cosine >= 1.0 ){ if (cosine >= 1.0) {
angle = 0.0f; angle = 0.0f;
} else if( cosine <= -1.0 ){ } else if (cosine <= -1.0) {
angle = RADIAN*PI_M; angle = RADIAN*PI_M;
} else { } else {
angle = RADIAN*acos(cosine); angle = RADIAN*acos(cosine);
} }
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "getPrefactorsGivenAngleCosine: cosine=%10.3e angle=%10.3e ideal=%10.3e\n", cosine, angle, idealAngle );
(void) fflush( log );
}
#endif
double deltaIdeal = angle - idealAngle; double deltaIdeal = angle - idealAngle;
double deltaIdeal2 = deltaIdeal*deltaIdeal; double deltaIdeal2 = deltaIdeal*deltaIdeal;
double deltaIdeal3 = deltaIdeal*deltaIdeal2; double deltaIdeal3 = deltaIdeal*deltaIdeal2;
...@@ -102,11 +95,11 @@ static void getPrefactorsGivenAngleCosine( double cosine, double idealAngle, dou ...@@ -102,11 +95,11 @@ static void getPrefactorsGivenAngleCosine( double cosine, double idealAngle, dou
// deltaIdeal = r - r_0 // deltaIdeal = r - r_0
*dEdR = ( 2.0 + *dEdR = (2.0 +
3.0*cubicK* deltaIdeal + 3.0*cubicK* deltaIdeal +
4.0*quarticK*deltaIdeal2 + 4.0*quarticK*deltaIdeal2 +
5.0*penticK* deltaIdeal3 + 5.0*penticK* deltaIdeal3 +
6.0*sexticK* deltaIdeal4 ); 6.0*sexticK* deltaIdeal4 );
*dEdR *= RADIAN*quadraticK*deltaIdeal; *dEdR *= RADIAN*quadraticK*deltaIdeal;
...@@ -121,30 +114,22 @@ static void getPrefactorsGivenAngleCosine( double cosine, double idealAngle, dou ...@@ -121,30 +114,22 @@ static void getPrefactorsGivenAngleCosine( double cosine, double idealAngle, dou
} }
static void computeAmoebaAngleForce(int bondIndex, std::vector<Vec3>& positions, AmoebaAngleForce& amoebaAngleForce, static void computeAmoebaAngleForce(int bondIndex, std::vector<Vec3>& positions, AmoebaAngleForce& amoebaAngleForce,
std::vector<Vec3>& forces, double* energy, FILE* log ) { std::vector<Vec3>& forces, double* energy) {
int particle1, particle2, particle3; int particle1, particle2, particle3;
double idealAngle; double idealAngle;
double quadraticK; double quadraticK;
amoebaAngleForce.getAngleParameters(bondIndex, particle1, particle2, particle3, idealAngle, quadraticK ); amoebaAngleForce.getAngleParameters(bondIndex, particle1, particle2, particle3, idealAngle, quadraticK);
double cubicK = amoebaAngleForce.getAmoebaGlobalAngleCubic(); double cubicK = amoebaAngleForce.getAmoebaGlobalAngleCubic();
double quarticK = amoebaAngleForce.getAmoebaGlobalAngleQuartic(); double quarticK = amoebaAngleForce.getAmoebaGlobalAngleQuartic();
double penticK = amoebaAngleForce.getAmoebaGlobalAnglePentic(); double penticK = amoebaAngleForce.getAmoebaGlobalAnglePentic();
double sexticK = amoebaAngleForce.getAmoebaGlobalAngleSextic(); double sexticK = amoebaAngleForce.getAmoebaGlobalAngleSextic();
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaAngleForce: bond %d [%d %d %d] ang=%10.3f k=%10.3f [%10.3e %10.3e %10.3e %10.3e]\n",
bondIndex, particle1, particle2, particle3, idealAngle, quadraticK, cubicK, quarticK, penticK, sexticK );
(void) fflush( log );
}
#endif
double deltaR[2][3]; double deltaR[2][3];
double r2_0 = 0.0; double r2_0 = 0.0;
double r2_1 = 0.0; double r2_1 = 0.0;
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
deltaR[0][ii] = positions[particle1][ii] - positions[particle2][ii]; deltaR[0][ii] = positions[particle1][ii] - positions[particle2][ii];
r2_0 += deltaR[0][ii]*deltaR[0][ii]; r2_0 += deltaR[0][ii]*deltaR[0][ii];
...@@ -155,33 +140,26 @@ static void computeAmoebaAngleForce(int bondIndex, std::vector<Vec3>& positions ...@@ -155,33 +140,26 @@ static void computeAmoebaAngleForce(int bondIndex, std::vector<Vec3>& positions
} }
double pVector[3]; double pVector[3];
crossProductVector3( deltaR[0], deltaR[1], pVector ); crossProductVector3(deltaR[0], deltaR[1], pVector);
double rp = sqrt( pVector[0]*pVector[0] + pVector[1]*pVector[1] + pVector[2]*pVector[2] ); double rp = sqrt(pVector[0]*pVector[0] + pVector[1]*pVector[1] + pVector[2]*pVector[2]);
if( rp < 1.0e-06 ){ if (rp < 1.0e-06) {
rp = 1.0e-06; rp = 1.0e-06;
} }
double dot = deltaR[0][0]*deltaR[1][0] + deltaR[0][1]*deltaR[1][1] + deltaR[0][2]*deltaR[1][2]; double dot = deltaR[0][0]*deltaR[1][0] + deltaR[0][1]*deltaR[1][1] + deltaR[0][2]*deltaR[1][2];
double cosine = dot/sqrt(r2_0*r2_1); double cosine = dot/sqrt(r2_0*r2_1);
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "dot=%10.3e r2_0=%10.3e r2_1=%10.3e\n", dot, r2_0, r2_1 );
(void) fflush( log );
}
#endif
double dEdR; double dEdR;
double energyTerm; double energyTerm;
getPrefactorsGivenAngleCosine( cosine, idealAngle, quadraticK, cubicK, getPrefactorsGivenAngleCosine(cosine, idealAngle, quadraticK, cubicK,
quarticK, penticK, sexticK, &dEdR, &energyTerm, log ); quarticK, penticK, sexticK, &dEdR, &energyTerm);
double termA = -dEdR/(r2_0*rp); double termA = -dEdR/(r2_0*rp);
double termC = dEdR/(r2_1*rp); double termC = dEdR/(r2_1*rp);
double deltaCrossP[3][3]; double deltaCrossP[3][3];
crossProductVector3( deltaR[0], pVector, deltaCrossP[0] ); crossProductVector3(deltaR[0], pVector, deltaCrossP[0]);
crossProductVector3( deltaR[1], pVector, deltaCrossP[2] ); crossProductVector3(deltaR[1], pVector, deltaCrossP[2]);
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
deltaCrossP[0][ii] *= termA; deltaCrossP[0][ii] *= termA;
deltaCrossP[2][ii] *= termC; deltaCrossP[2][ii] *= termC;
deltaCrossP[1][ii] = -1.0*(deltaCrossP[0][ii] + deltaCrossP[2][ii]); deltaCrossP[1][ii] = -1.0*(deltaCrossP[0][ii] + deltaCrossP[2][ii]);
...@@ -202,72 +180,51 @@ static void computeAmoebaAngleForce(int bondIndex, std::vector<Vec3>& positions ...@@ -202,72 +180,51 @@ static void computeAmoebaAngleForce(int bondIndex, std::vector<Vec3>& positions
*energy += energyTerm; *energy += energyTerm;
} }
static void computeAmoebaAngleForces( Context& context, AmoebaAngleForce& amoebaAngleForce, static void computeAmoebaAngleForces(Context& context, AmoebaAngleForce& amoebaAngleForce,
std::vector<Vec3>& expectedForces, double* expectedEnergy, FILE* log ) { std::vector<Vec3>& expectedForces, double* expectedEnergy) {
// get positions and zero forces // get positions and zero forces
State state = context.getState(State::Positions); State state = context.getState(State::Positions);
std::vector<Vec3> positions = state.getPositions(); std::vector<Vec3> positions = state.getPositions();
expectedForces.resize( positions.size() ); expectedForces.resize(positions.size());
for( unsigned int ii = 0; ii < expectedForces.size(); ii++ ){ for (unsigned int ii = 0; ii < expectedForces.size(); ii++) {
expectedForces[ii][0] = expectedForces[ii][1] = expectedForces[ii][2] = 0.0; expectedForces[ii][0] = expectedForces[ii][1] = expectedForces[ii][2] = 0.0;
} }
// calculates forces/energy // calculates forces/energy
*expectedEnergy = 0.0; *expectedEnergy = 0.0;
for( int ii = 0; ii < amoebaAngleForce.getNumAngles(); ii++ ){ for (int ii = 0; ii < amoebaAngleForce.getNumAngles(); ii++) {
computeAmoebaAngleForce(ii, positions, amoebaAngleForce, expectedForces, expectedEnergy, log ); computeAmoebaAngleForce(ii, positions, amoebaAngleForce, expectedForces, expectedEnergy);
}
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaAngleForces: expected energy=%14.7e\n", *expectedEnergy );
for( unsigned int ii = 0; ii < positions.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e]\n", ii, expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2] );
}
(void) fflush( log );
} }
#endif
return; return;
} }
void compareWithExpectedForceAndEnergy( Context& context, AmoebaAngleForce& amoebaAngleForce, void compareWithExpectedForceAndEnergy(Context& context, AmoebaAngleForce& amoebaAngleForce,
double tolerance, const std::string& idString, FILE* log) { double tolerance, const std::string& idString) {
std::vector<Vec3> expectedForces; std::vector<Vec3> expectedForces;
double expectedEnergy; double expectedEnergy;
computeAmoebaAngleForces( context, amoebaAngleForce, expectedForces, &expectedEnergy, log ); computeAmoebaAngleForces(context, amoebaAngleForce, expectedForces, &expectedEnergy);
State state = context.getState(State::Forces | State::Energy); State state = context.getState(State::Forces | State::Energy);
const std::vector<Vec3> forces = state.getForces(); const std::vector<Vec3> forces = state.getForces();
#ifdef AMOEBA_DEBUG for (unsigned int ii = 0; ii < forces.size(); ii++) {
if( log ){ ASSERT_EQUAL_VEC(expectedForces[ii], forces[ii], tolerance);
(void) fprintf( log, "computeAmoebaAngleForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy() );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2], forces[ii][0], forces[ii][1], forces[ii][2] );
}
(void) fflush( log );
}
#endif
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
ASSERT_EQUAL_VEC( expectedForces[ii], forces[ii], tolerance );
} }
ASSERT_EQUAL_TOL( expectedEnergy, state.getPotentialEnergy(), tolerance ); ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), tolerance);
} }
void testOneAngle( FILE* log ) { void testOneAngle() {
System system; System system;
int numberOfParticles = 3; int numberOfParticles = 3;
for( int ii = 0; ii < numberOfParticles; ii++ ){ for (int ii = 0; ii < numberOfParticles; ii++) {
system.addParticle(1.0); system.addParticle(1.0);
} }
...@@ -289,7 +246,7 @@ void testOneAngle( FILE* log ) { ...@@ -289,7 +246,7 @@ void testOneAngle( FILE* log ) {
amoebaAngleForce->setAmoebaGlobalAngleSextic(sexticK); amoebaAngleForce->setAmoebaGlobalAngleSextic(sexticK);
system.addForce(amoebaAngleForce); system.addForce(amoebaAngleForce);
Context context(system, integrator, Platform::getPlatformByName( "CUDA")); Context context(system, integrator, Platform::getPlatformByName("CUDA"));
std::vector<Vec3> positions(numberOfParticles); std::vector<Vec3> positions(numberOfParticles);
...@@ -298,7 +255,7 @@ void testOneAngle( FILE* log ) { ...@@ -298,7 +255,7 @@ void testOneAngle( FILE* log ) {
positions[2] = Vec3(0, 0, 1); positions[2] = Vec3(0, 0, 1);
context.setPositions(positions); context.setPositions(positions);
compareWithExpectedForceAndEnergy( context, *amoebaAngleForce, TOL, "testOneAngle", log ); compareWithExpectedForceAndEnergy(context, *amoebaAngleForce, TOL, "testOneAngle");
// Try changing the angle parameters and make sure it's still correct. // Try changing the angle parameters and make sure it's still correct.
...@@ -306,14 +263,14 @@ void testOneAngle( FILE* log ) { ...@@ -306,14 +263,14 @@ void testOneAngle( FILE* log ) {
bool exceptionThrown = false; bool exceptionThrown = false;
try { try {
// This should throw an exception. // This should throw an exception.
compareWithExpectedForceAndEnergy( context, *amoebaAngleForce, TOL, "testOneAngle", log ); compareWithExpectedForceAndEnergy(context, *amoebaAngleForce, TOL, "testOneAngle");
} }
catch (std::exception ex) { catch (std::exception ex) {
exceptionThrown = true; exceptionThrown = true;
} }
ASSERT(exceptionThrown); ASSERT(exceptionThrown);
amoebaAngleForce->updateParametersInContext(context); amoebaAngleForce->updateParametersInContext(context);
compareWithExpectedForceAndEnergy( context, *amoebaAngleForce, TOL, "testOneAngle", log ); compareWithExpectedForceAndEnergy(context, *amoebaAngleForce, TOL, "testOneAngle");
} }
int main(int argc, char* argv[]) { int main(int argc, char* argv[]) {
...@@ -322,8 +279,7 @@ int main(int argc, char* argv[]) { ...@@ -322,8 +279,7 @@ int main(int argc, char* argv[]) {
registerAmoebaCudaKernelFactories(); registerAmoebaCudaKernelFactories();
if (argc > 1) if (argc > 1)
Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1])); Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1]));
FILE* log = NULL; testOneAngle();
testOneAngle( log );
} catch(const std::exception& e) { } catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl; std::cout << "exception: " << e.what() << std::endl;
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaBondForce.cpp
View file @ f11bafd9
...@@ -48,22 +48,22 @@ extern "C" void registerAmoebaCudaKernelFactories(); ...@@ -48,22 +48,22 @@ extern "C" void registerAmoebaCudaKernelFactories();
const double TOL = 1e-5; const double TOL = 1e-5;
static void computeAmoebaBondForce(int bondIndex, std::vector<Vec3>& positions, AmoebaBondForce& amoebaBondForce, static void computeAmoebaBondForce(int bondIndex, std::vector<Vec3>& positions, AmoebaBondForce& amoebaBondForce,
std::vector<Vec3>& forces, double* energy ) { std::vector<Vec3>& forces, double* energy) {
int particle1, particle2; int particle1, particle2;
double bondLength; double bondLength;
double quadraticK; double quadraticK;
double cubicK = amoebaBondForce.getAmoebaGlobalBondCubic(); double cubicK = amoebaBondForce.getAmoebaGlobalBondCubic();
double quarticK = amoebaBondForce.getAmoebaGlobalBondQuartic(); double quarticK = amoebaBondForce.getAmoebaGlobalBondQuartic();
amoebaBondForce.getBondParameters(bondIndex, particle1, particle2, bondLength, quadraticK ); amoebaBondForce.getBondParameters(bondIndex, particle1, particle2, bondLength, quadraticK);
double deltaR[3]; double deltaR[3];
double r2 = 0.0; double r2 = 0.0;
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
deltaR[ii] = positions[particle2][ii] - positions[particle1][ii]; deltaR[ii] = positions[particle2][ii] - positions[particle1][ii];
r2 += deltaR[ii]*deltaR[ii]; r2 += deltaR[ii]*deltaR[ii];
} }
double r = sqrt( r2 ); double r = sqrt(r2);
double bondDelta = (r - bondLength); double bondDelta = (r - bondLength);
double bondDelta2 = bondDelta*bondDelta; double bondDelta2 = bondDelta*bondDelta;
...@@ -83,66 +83,43 @@ static void computeAmoebaBondForce(int bondIndex, std::vector<Vec3>& positions, ...@@ -83,66 +83,43 @@ static void computeAmoebaBondForce(int bondIndex, std::vector<Vec3>& positions,
} }
static void computeAmoebaBondForces( Context& context, AmoebaBondForce& amoebaBondForce, static void computeAmoebaBondForces(Context& context, AmoebaBondForce& amoebaBondForce,
std::vector<Vec3>& expectedForces, double* expectedEnergy, FILE* log ) { std::vector<Vec3>& expectedForces, double* expectedEnergy) {
// get positions and zero forces // get positions and zero forces
State state = context.getState(State::Positions); State state = context.getState(State::Positions);
std::vector<Vec3> positions = state.getPositions(); std::vector<Vec3> positions = state.getPositions();
expectedForces.resize( positions.size() ); expectedForces.resize(positions.size());
for( unsigned int ii = 0; ii < expectedForces.size(); ii++ ){ for (unsigned int ii = 0; ii < expectedForces.size(); ii++) {
expectedForces[ii][0] = expectedForces[ii][1] = expectedForces[ii][2] = 0.0; expectedForces[ii][0] = expectedForces[ii][1] = expectedForces[ii][2] = 0.0;
} }
// calculates forces/energy // calculates forces/energy
*expectedEnergy = 0.0; *expectedEnergy = 0.0;
for( int ii = 0; ii < amoebaBondForce.getNumBonds(); ii++ ){ for (int ii = 0; ii < amoebaBondForce.getNumBonds(); ii++) {
computeAmoebaBondForce(ii, positions, amoebaBondForce, expectedForces, expectedEnergy ); computeAmoebaBondForce(ii, positions, amoebaBondForce, expectedForces, expectedEnergy);
} }
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaBondForces: expected energy=%15.7e\n", *expectedEnergy );
for( unsigned int ii = 0; ii < positions.size(); ii++ ){
(void) fprintf( log, "%6u [%15.7e %15.7e %15.7e]\n", ii, expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2] );
}
(void) fflush( log );
}
#endif
return;
} }
void compareWithExpectedForceAndEnergy( Context& context, AmoebaBondForce& amoebaBondForce, double tolerance, const std::string& idString, FILE* log) { void compareWithExpectedForceAndEnergy(Context& context, AmoebaBondForce& amoebaBondForce, double tolerance, const std::string& idString) {
std::vector<Vec3> expectedForces; std::vector<Vec3> expectedForces;
double expectedEnergy; double expectedEnergy;
computeAmoebaBondForces( context, amoebaBondForce, expectedForces, &expectedEnergy, NULL ); computeAmoebaBondForces(context, amoebaBondForce, expectedForces, &expectedEnergy);
State state = context.getState(State::Forces | State::Energy); State state = context.getState(State::Forces | State::Energy);
const std::vector<Vec3> forces = state.getForces(); const std::vector<Vec3> forces = state.getForces();
#ifdef AMOEBA_DEBUG for (unsigned int ii = 0; ii < forces.size(); ii++) {
if( log ){ ASSERT_EQUAL_VEC(expectedForces[ii], forces[ii], tolerance);
(void) fprintf( log, "computeAmoebaBondForces: expected energy=%15.7e %15.7e\n", expectedEnergy, state.getPotentialEnergy() );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
(void) fprintf( log, "%6u [%15.7e %15.7e %15.7e] [%15.7e %15.7e %15.7e]\n", ii,
expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2], forces[ii][0], forces[ii][1], forces[ii][2] );
}
(void) fflush( log );
}
#endif
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
ASSERT_EQUAL_VEC( expectedForces[ii], forces[ii], tolerance );
} }
ASSERT_EQUAL_TOL( expectedEnergy, state.getPotentialEnergy(), tolerance ); ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), tolerance);
} }
void testOneBond( FILE* log ) { void testOneBond() {
System system; System system;
...@@ -157,22 +134,22 @@ void testOneBond( FILE* log ) { ...@@ -157,22 +134,22 @@ void testOneBond( FILE* log ) {
double quadraticK = 1.0; double quadraticK = 1.0;
double cubicK = 2.0; double cubicK = 2.0;
double quarticicK = 3.0; double quarticicK = 3.0;
amoebaBondForce->setAmoebaGlobalBondCubic( cubicK ); amoebaBondForce->setAmoebaGlobalBondCubic(cubicK);
amoebaBondForce->setAmoebaGlobalBondQuartic( quarticicK ); amoebaBondForce->setAmoebaGlobalBondQuartic(quarticicK);
amoebaBondForce->addBond(0, 1, bondLength, quadraticK); amoebaBondForce->addBond(0, 1, bondLength, quadraticK);
system.addForce(amoebaBondForce); system.addForce(amoebaBondForce);
Context context(system, integrator, Platform::getPlatformByName( "CUDA")); Context context(system, integrator, Platform::getPlatformByName("CUDA"));
std::vector<Vec3> positions(2); std::vector<Vec3> positions(2);
positions[0] = Vec3(0, 1, 0); positions[0] = Vec3(0, 1, 0);
positions[1] = Vec3(0, 0, 0); positions[1] = Vec3(0, 0, 0);
context.setPositions(positions); context.setPositions(positions);
compareWithExpectedForceAndEnergy( context, *amoebaBondForce, TOL, "testOneBond", log ); compareWithExpectedForceAndEnergy(context, *amoebaBondForce, TOL, "testOneBond");
} }
void testTwoBond( FILE* log ) { void testTwoBond() {
System system; System system;
...@@ -188,14 +165,14 @@ void testTwoBond( FILE* log ) { ...@@ -188,14 +165,14 @@ void testTwoBond( FILE* log ) {
double quadraticK = 1.0; double quadraticK = 1.0;
double cubicK = 2.0; double cubicK = 2.0;
double quarticicK = 3.0; double quarticicK = 3.0;
amoebaBondForce->setAmoebaGlobalBondCubic( cubicK ); amoebaBondForce->setAmoebaGlobalBondCubic(cubicK);
amoebaBondForce->setAmoebaGlobalBondQuartic( quarticicK ); amoebaBondForce->setAmoebaGlobalBondQuartic(quarticicK);
amoebaBondForce->addBond(0, 1, bondLength, quadraticK); amoebaBondForce->addBond(0, 1, bondLength, quadraticK);
amoebaBondForce->addBond(1, 2, bondLength, quadraticK); amoebaBondForce->addBond(1, 2, bondLength, quadraticK);
system.addForce(amoebaBondForce); system.addForce(amoebaBondForce);
Context context(system, integrator, Platform::getPlatformByName( "CUDA")); Context context(system, integrator, Platform::getPlatformByName("CUDA"));
//Context context(system, integrator, platform ); //Context context(system, integrator, platform);
std::vector<Vec3> positions(3); std::vector<Vec3> positions(3);
positions[0] = Vec3(0, 1, 0); positions[0] = Vec3(0, 1, 0);
...@@ -203,7 +180,7 @@ void testTwoBond( FILE* log ) { ...@@ -203,7 +180,7 @@ void testTwoBond( FILE* log ) {
positions[2] = Vec3(1, 0, 1); positions[2] = Vec3(1, 0, 1);
context.setPositions(positions); context.setPositions(positions);
compareWithExpectedForceAndEnergy( context, *amoebaBondForce, TOL, "testTwoBond", log ); compareWithExpectedForceAndEnergy(context, *amoebaBondForce, TOL, "testTwoBond");
// Try changing the bond parameters and make sure it's still correct. // Try changing the bond parameters and make sure it's still correct.
...@@ -212,14 +189,14 @@ void testTwoBond( FILE* log ) { ...@@ -212,14 +189,14 @@ void testTwoBond( FILE* log ) {
bool exceptionThrown = false; bool exceptionThrown = false;
try { try {
// This should throw an exception. // This should throw an exception.
compareWithExpectedForceAndEnergy( context, *amoebaBondForce, TOL, "testTwoBond", log ); compareWithExpectedForceAndEnergy(context, *amoebaBondForce, TOL, "testTwoBond");
} }
catch (std::exception ex) { catch (std::exception ex) {
exceptionThrown = true; exceptionThrown = true;
} }
ASSERT(exceptionThrown); ASSERT(exceptionThrown);
amoebaBondForce->updateParametersInContext(context); amoebaBondForce->updateParametersInContext(context);
compareWithExpectedForceAndEnergy( context, *amoebaBondForce, TOL, "testTwoBond", log ); compareWithExpectedForceAndEnergy(context, *amoebaBondForce, TOL, "testTwoBond");
} }
int main(int argc, char* argv[]) { int main(int argc, char* argv[]) {
...@@ -228,8 +205,7 @@ int main(int argc, char* argv[]) { ...@@ -228,8 +205,7 @@ int main(int argc, char* argv[]) {
registerAmoebaCudaKernelFactories(); registerAmoebaCudaKernelFactories();
if (argc > 1) if (argc > 1)
Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1])); Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1]));
FILE* log = NULL; testTwoBond();
testTwoBond( log );
} catch(const std::exception& e) { } catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl; std::cout << "exception: " << e.what() << std::endl;
std::cout << "FAIL - ERROR. Test failed." << std::endl; std::cout << "FAIL - ERROR. Test failed." << std::endl;
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaGeneralizedKirkwoodForce.cpp
View file @ f11bafd9
...@@ -64,10 +64,10 @@ static void setupMultipoleAmmonia(System& system, AmoebaGeneralizedKirkwoodForce ...@@ -64,10 +64,10 @@ static void setupMultipoleAmmonia(System& system, AmoebaGeneralizedKirkwoodForce
AmoebaMultipoleForce* amoebaMultipoleForce = new AmoebaMultipoleForce();; AmoebaMultipoleForce* amoebaMultipoleForce = new AmoebaMultipoleForce();;
int numberOfParticles = 8; int numberOfParticles = 8;
   
amoebaMultipoleForce->setNonbondedMethod( AmoebaMultipoleForce::NoCutoff ); amoebaMultipoleForce->setNonbondedMethod(AmoebaMultipoleForce::NoCutoff);
amoebaMultipoleForce->setPolarizationType( polarizationType ); amoebaMultipoleForce->setPolarizationType(polarizationType);
amoebaMultipoleForce->setMutualInducedTargetEpsilon( 1.0e-06 ); amoebaMultipoleForce->setMutualInducedTargetEpsilon(1.0e-06);
amoebaMultipoleForce->setMutualInducedMaxIterations( 500 ); amoebaMultipoleForce->setMutualInducedMaxIterations(500);
   
std::vector<double> nitrogenMolecularDipole(3); std::vector<double> nitrogenMolecularDipole(3);
std::vector<double> nitrogenMolecularQuadrupole(9); std::vector<double> nitrogenMolecularQuadrupole(9);
...@@ -88,8 +88,8 @@ static void setupMultipoleAmmonia(System& system, AmoebaGeneralizedKirkwoodForce ...@@ -88,8 +88,8 @@ static void setupMultipoleAmmonia(System& system, AmoebaGeneralizedKirkwoodForce
   
// first N // first N
   
system.addParticle( 1.4007000e+01 ); system.addParticle(1.4007000e+01);
amoebaMultipoleForce->addMultipole( -5.7960000e-01, nitrogenMolecularDipole, nitrogenMolecularQuadrupole, 2, 1, 2, 3, 3.9000000e-01, 3.1996314e-01, 1.0730000e-03 ); amoebaMultipoleForce->addMultipole(-5.7960000e-01, nitrogenMolecularDipole, nitrogenMolecularQuadrupole, 2, 1, 2, 3, 3.9000000e-01, 3.1996314e-01, 1.0730000e-03);
   
// 3 H attached to first N // 3 H attached to first N
   
...@@ -109,181 +109,181 @@ static void setupMultipoleAmmonia(System& system, AmoebaGeneralizedKirkwoodForce ...@@ -109,181 +109,181 @@ static void setupMultipoleAmmonia(System& system, AmoebaGeneralizedKirkwoodForce
hydrogenMolecularQuadrupole[7] = 0.0000000e+00; hydrogenMolecularQuadrupole[7] = 0.0000000e+00;
hydrogenMolecularQuadrupole[8] = 2.4549167e-06; hydrogenMolecularQuadrupole[8] = 2.4549167e-06;
   
system.addParticle( 1.0080000e+00 ); system.addParticle(1.0080000e+00);
system.addParticle( 1.0080000e+00 ); system.addParticle(1.0080000e+00);
system.addParticle( 1.0080000e+00 ); system.addParticle(1.0080000e+00);
amoebaMultipoleForce->addMultipole( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 0, 2, 3, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04 ); amoebaMultipoleForce->addMultipole(1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 0, 2, 3, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04);
amoebaMultipoleForce->addMultipole( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 0, 1, 3, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04 ); amoebaMultipoleForce->addMultipole(1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 0, 1, 3, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04);
amoebaMultipoleForce->addMultipole( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 0, 1, 2, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04 ); amoebaMultipoleForce->addMultipole(1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 0, 1, 2, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04);
   
// second N // second N
   
system.addParticle( 1.4007000e+01 ); system.addParticle( 1.4007000e+01);
amoebaMultipoleForce->addMultipole( -5.7960000e-01, nitrogenMolecularDipole, nitrogenMolecularQuadrupole, 2, 5, 6, 7, 3.9000000e-01, 3.1996314e-01, 1.0730000e-03 ); amoebaMultipoleForce->addMultipole( -5.7960000e-01, nitrogenMolecularDipole, nitrogenMolecularQuadrupole, 2, 5, 6, 7, 3.9000000e-01, 3.1996314e-01, 1.0730000e-03);
   
// 3 H attached to second N // 3 H attached to second N
   
system.addParticle( 1.0080000e+00 ); system.addParticle( 1.0080000e+00);
system.addParticle( 1.0080000e+00 ); system.addParticle( 1.0080000e+00);
system.addParticle( 1.0080000e+00 ); system.addParticle( 1.0080000e+00);
amoebaMultipoleForce->addMultipole( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 4, 6, 7, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04 ); amoebaMultipoleForce->addMultipole( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 4, 6, 7, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04);
amoebaMultipoleForce->addMultipole( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 4, 5, 7, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04 ); amoebaMultipoleForce->addMultipole( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 4, 5, 7, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04);
amoebaMultipoleForce->addMultipole( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 4, 5, 6, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04 ); amoebaMultipoleForce->addMultipole( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 4, 5, 6, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04);
   
// covalent maps // covalent maps
   
std::vector< int > covalentMap; std::vector< int > covalentMap;
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 1 ); covalentMap.push_back(1);
covalentMap.push_back( 2 ); covalentMap.push_back(2);
covalentMap.push_back( 3 ); covalentMap.push_back(3);
amoebaMultipoleForce->setCovalentMap( 0, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap ); amoebaMultipoleForce->setCovalentMap(0, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 0 ); covalentMap.push_back(0);
covalentMap.push_back( 1 ); covalentMap.push_back(1);
covalentMap.push_back( 2 ); covalentMap.push_back(2);
covalentMap.push_back( 3 ); covalentMap.push_back(3);
amoebaMultipoleForce->setCovalentMap( 0, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap ); amoebaMultipoleForce->setCovalentMap(0, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 0 ); covalentMap.push_back(0);
amoebaMultipoleForce->setCovalentMap( 1, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap ); amoebaMultipoleForce->setCovalentMap(1, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 2 ); covalentMap.push_back(2);
covalentMap.push_back( 3 ); covalentMap.push_back(3);
amoebaMultipoleForce->setCovalentMap( 1, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap ); amoebaMultipoleForce->setCovalentMap(1, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 0 ); covalentMap.push_back(0);
covalentMap.push_back( 1 ); covalentMap.push_back(1);
covalentMap.push_back( 2 ); covalentMap.push_back(2);
covalentMap.push_back( 3 ); covalentMap.push_back(3);
amoebaMultipoleForce->setCovalentMap( 1, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap ); amoebaMultipoleForce->setCovalentMap(1, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 0 ); covalentMap.push_back(0);
amoebaMultipoleForce->setCovalentMap( 2, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap ); amoebaMultipoleForce->setCovalentMap(2, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 1 ); covalentMap.push_back(1);
covalentMap.push_back( 3 ); covalentMap.push_back(3);
amoebaMultipoleForce->setCovalentMap( 2, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap ); amoebaMultipoleForce->setCovalentMap(2, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 0 ); covalentMap.push_back(0);
covalentMap.push_back( 1 ); covalentMap.push_back(1);
covalentMap.push_back( 2 ); covalentMap.push_back(2);
covalentMap.push_back( 3 ); covalentMap.push_back(3);
amoebaMultipoleForce->setCovalentMap( 2, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap ); amoebaMultipoleForce->setCovalentMap(2, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 0 ); covalentMap.push_back(0);
amoebaMultipoleForce->setCovalentMap( 3, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap ); amoebaMultipoleForce->setCovalentMap(3, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 1 ); covalentMap.push_back(1);
covalentMap.push_back( 2 ); covalentMap.push_back(2);
amoebaMultipoleForce->setCovalentMap( 3, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap ); amoebaMultipoleForce->setCovalentMap(3, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 0 ); covalentMap.push_back(0);
covalentMap.push_back( 1 ); covalentMap.push_back(1);
covalentMap.push_back( 2 ); covalentMap.push_back(2);
covalentMap.push_back( 3 ); covalentMap.push_back(3);
amoebaMultipoleForce->setCovalentMap( 3, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap ); amoebaMultipoleForce->setCovalentMap(3, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 5 ); covalentMap.push_back(5);
covalentMap.push_back( 6 ); covalentMap.push_back(6);
covalentMap.push_back( 7 ); covalentMap.push_back(7);
amoebaMultipoleForce->setCovalentMap( 4, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap ); amoebaMultipoleForce->setCovalentMap(4, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 4 ); covalentMap.push_back(4);
covalentMap.push_back( 5 ); covalentMap.push_back(5);
covalentMap.push_back( 6 ); covalentMap.push_back(6);
covalentMap.push_back( 7 ); covalentMap.push_back(7);
amoebaMultipoleForce->setCovalentMap( 4, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap ); amoebaMultipoleForce->setCovalentMap(4, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 4 ); covalentMap.push_back(4);
amoebaMultipoleForce->setCovalentMap( 5, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap ); amoebaMultipoleForce->setCovalentMap(5, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 6 ); covalentMap.push_back(6);
covalentMap.push_back( 7 ); covalentMap.push_back(7);
amoebaMultipoleForce->setCovalentMap( 5, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap ); amoebaMultipoleForce->setCovalentMap(5, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 4 ); covalentMap.push_back(4);
covalentMap.push_back( 5 ); covalentMap.push_back(5);
covalentMap.push_back( 6 ); covalentMap.push_back(6);
covalentMap.push_back( 7 ); covalentMap.push_back(7);
amoebaMultipoleForce->setCovalentMap( 5, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap ); amoebaMultipoleForce->setCovalentMap(5, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 4 ); covalentMap.push_back(4);
amoebaMultipoleForce->setCovalentMap( 6, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap ); amoebaMultipoleForce->setCovalentMap(6, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 5 ); covalentMap.push_back(5);
covalentMap.push_back( 7 ); covalentMap.push_back(7);
amoebaMultipoleForce->setCovalentMap( 6, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap ); amoebaMultipoleForce->setCovalentMap(6, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 4 ); covalentMap.push_back(4);
covalentMap.push_back( 5 ); covalentMap.push_back(5);
covalentMap.push_back( 6 ); covalentMap.push_back(6);
covalentMap.push_back( 7 ); covalentMap.push_back(7);
amoebaMultipoleForce->setCovalentMap( 6, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap ); amoebaMultipoleForce->setCovalentMap(6, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 4 ); covalentMap.push_back(4);
amoebaMultipoleForce->setCovalentMap( 7, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap ); amoebaMultipoleForce->setCovalentMap(7, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 5 ); covalentMap.push_back(5);
covalentMap.push_back( 6 ); covalentMap.push_back(6);
amoebaMultipoleForce->setCovalentMap( 7, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap ); amoebaMultipoleForce->setCovalentMap(7, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 4 ); covalentMap.push_back(4);
covalentMap.push_back( 5 ); covalentMap.push_back(5);
covalentMap.push_back( 6 ); covalentMap.push_back(6);
covalentMap.push_back( 7 ); covalentMap.push_back(7);
amoebaMultipoleForce->setCovalentMap( 7, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap ); amoebaMultipoleForce->setCovalentMap(7, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap);
system.addForce(amoebaMultipoleForce); system.addForce(amoebaMultipoleForce);
   
// GK force // GK force
   
amoebaGeneralizedKirkwoodForce->setSolventDielectric( 7.8300000e+01 ); amoebaGeneralizedKirkwoodForce->setSolventDielectric( 7.8300000e+01);
amoebaGeneralizedKirkwoodForce->setSoluteDielectric( 1.0000000e+00 ); amoebaGeneralizedKirkwoodForce->setSoluteDielectric( 1.0000000e+00);
amoebaGeneralizedKirkwoodForce->setIncludeCavityTerm( includeCavityTerm ); amoebaGeneralizedKirkwoodForce->setIncludeCavityTerm(includeCavityTerm);
   
// addParticle: charge, radius, scalingFactor // addParticle: charge, radius, scalingFactor
   
for( unsigned int ii = 0; ii < 2; ii++ ){ for (unsigned int ii = 0; ii < 2; ii++) {
amoebaGeneralizedKirkwoodForce->addParticle( -5.7960000e-01, 1.5965000e-01, 6.9000000e-01 ); amoebaGeneralizedKirkwoodForce->addParticle( -5.7960000e-01, 1.5965000e-01, 6.9000000e-01);
amoebaGeneralizedKirkwoodForce->addParticle( 1.9320000e-01, 1.2360000e-01, 6.9000000e-01 ); amoebaGeneralizedKirkwoodForce->addParticle( 1.9320000e-01, 1.2360000e-01, 6.9000000e-01);
amoebaGeneralizedKirkwoodForce->addParticle( 1.9320000e-01, 1.2360000e-01, 6.9000000e-01 ); amoebaGeneralizedKirkwoodForce->addParticle( 1.9320000e-01, 1.2360000e-01, 6.9000000e-01);
amoebaGeneralizedKirkwoodForce->addParticle( 1.9320000e-01, 1.2360000e-01, 6.9000000e-01 ); amoebaGeneralizedKirkwoodForce->addParticle( 1.9320000e-01, 1.2360000e-01, 6.9000000e-01);
} }
system.addForce(amoebaGeneralizedKirkwoodForce); system.addForce(amoebaGeneralizedKirkwoodForce);
} }
   
static void getForcesEnergyMultipoleAmmonia(Context& context, std::vector<Vec3>& forces, double& energy, FILE* log) { static void getForcesEnergyMultipoleAmmonia(Context& context, std::vector<Vec3>& forces, double& energy) {
std::vector<Vec3> positions(context.getSystem().getNumParticles()); std::vector<Vec3> positions(context.getSystem().getNumParticles());
   
positions[0] = Vec3( 1.5927280e-01, 1.7000000e-06, 1.6491000e-03 ); positions[0] = Vec3( 1.5927280e-01, 1.7000000e-06, 1.6491000e-03);
positions[1] = Vec3( 2.0805540e-01, -8.1258800e-02, 3.7282500e-02 ); positions[1] = Vec3( 2.0805540e-01, -8.1258800e-02, 3.7282500e-02);
positions[2] = Vec3( 2.0843610e-01, 8.0953200e-02, 3.7462200e-02 ); positions[2] = Vec3( 2.0843610e-01, 8.0953200e-02, 3.7462200e-02);
positions[3] = Vec3( 1.7280780e-01, 2.0730000e-04, -9.8741700e-02 ); positions[3] = Vec3( 1.7280780e-01, 2.0730000e-04, -9.8741700e-02);
positions[4] = Vec3( -1.6743680e-01, 1.5900000e-05, -6.6149000e-03 ); positions[4] = Vec3( -1.6743680e-01, 1.5900000e-05, -6.6149000e-03);
positions[5] = Vec3( -2.0428260e-01, 8.1071500e-02, 4.1343900e-02 ); positions[5] = Vec3( -2.0428260e-01, 8.1071500e-02, 4.1343900e-02);
positions[6] = Vec3( -6.7308300e-02, 1.2800000e-05, 1.0623300e-02 ); positions[6] = Vec3( -6.7308300e-02, 1.2800000e-05, 1.0623300e-02);
positions[7] = Vec3( -2.0426290e-01, -8.1231400e-02, 4.1033500e-02 ); positions[7] = Vec3( -2.0426290e-01, -8.1231400e-02, 4.1033500e-02);
   
context.setPositions(positions); context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy); State state = context.getState(State::Forces | State::Energy);
...@@ -293,8 +293,8 @@ static void getForcesEnergyMultipoleAmmonia(Context& context, std::vector<Vec3>& ...@@ -293,8 +293,8 @@ static void getForcesEnergyMultipoleAmmonia(Context& context, std::vector<Vec3>&
   
// setup for villin // setup for villin
   
static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::PolarizationType polarizationType, static void setupAndGetForcesEnergyMultipoleVillin(AmoebaMultipoleForce::PolarizationType polarizationType,
int includeCavityTerm, std::vector<Vec3>& forces, double& energy, FILE* log ){ int includeCavityTerm, std::vector<Vec3>& forces, double& energy) {
   
// beginning of Multipole setup // beginning of Multipole setup
   
...@@ -303,13 +303,13 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari ...@@ -303,13 +303,13 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari
AmoebaMultipoleForce* amoebaMultipoleForce = new AmoebaMultipoleForce();; AmoebaMultipoleForce* amoebaMultipoleForce = new AmoebaMultipoleForce();;
int numberOfParticles = 596; int numberOfParticles = 596;
   
amoebaMultipoleForce->setNonbondedMethod( AmoebaMultipoleForce::NoCutoff ); amoebaMultipoleForce->setNonbondedMethod(AmoebaMultipoleForce::NoCutoff);
amoebaMultipoleForce->setPolarizationType( polarizationType ); amoebaMultipoleForce->setPolarizationType(polarizationType);
amoebaMultipoleForce->setMutualInducedTargetEpsilon( 1.0e-06 ); amoebaMultipoleForce->setMutualInducedTargetEpsilon(1.0e-06);
amoebaMultipoleForce->setMutualInducedMaxIterations( 500 ); amoebaMultipoleForce->setMutualInducedMaxIterations(500);
   
for( unsigned int ii = 0; ii < numberOfParticles; ii++ ){ for (unsigned int ii = 0; ii < numberOfParticles; ii++) {
system.addParticle( 1.0 ); system.addParticle(1.0);
} }
   
static const double multipoleData[] = { static const double multipoleData[] = {
...@@ -927,7 +927,7 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari ...@@ -927,7 +927,7 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari
std::vector<double> quadrupole(9); std::vector<double> quadrupole(9);
unsigned int entriesPerParticle = 21; unsigned int entriesPerParticle = 21;
const double* data = multipoleData; const double* data = multipoleData;
for( unsigned int ii = 0; ii < numberOfParticles; ii++ ){ for (unsigned int ii = 0; ii < numberOfParticles; ii++) {
   
dipole[0] = data[dipoleIndex + 0]; dipole[0] = data[dipoleIndex + 0];
dipole[1] = data[dipoleIndex + 1]; dipole[1] = data[dipoleIndex + 1];
...@@ -943,9 +943,9 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari ...@@ -943,9 +943,9 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari
quadrupole[7] = data[quadrupoleIndex + 7]; quadrupole[7] = data[quadrupoleIndex + 7];
quadrupole[8] = data[quadrupoleIndex + 8]; quadrupole[8] = data[quadrupoleIndex + 8];
   
amoebaMultipoleForce->addMultipole( data[chargeIndex], dipole, quadrupole, static_cast<int>(data[axisTypeIndex]), amoebaMultipoleForce->addMultipole(data[chargeIndex], dipole, quadrupole, static_cast<int>(data[axisTypeIndex]),
static_cast<int>(data[multipoleAtomZIndex]), static_cast<int>(data[multipoleAtomXIndex]), static_cast<int>(data[multipoleAtomYIndex]), static_cast<int>(data[multipoleAtomZIndex]), static_cast<int>(data[multipoleAtomXIndex]), static_cast<int>(data[multipoleAtomYIndex]),
data[tholeIndex], data[dampingFactorIndex], data[polarityIndex] ); data[tholeIndex], data[dampingFactorIndex], data[polarityIndex]);
data += entriesPerParticle; data += entriesPerParticle;
} }
   
...@@ -5725,15 +5725,15 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari ...@@ -5725,15 +5725,15 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari
unsigned int covalentMapDataSize = sizeof(covalentMapData)/sizeof(int); unsigned int covalentMapDataSize = sizeof(covalentMapData)/sizeof(int);
   
unsigned int covalentIndex = 0; unsigned int covalentIndex = 0;
while( covalentIndex < covalentMapDataSize ){ while (covalentIndex < covalentMapDataSize) {
int particleIndex = covalentMapData[covalentIndex++]; int particleIndex = covalentMapData[covalentIndex++];
int typeIndex = covalentMapData[covalentIndex++]; int typeIndex = covalentMapData[covalentIndex++];
int entries = covalentMapData[covalentIndex++]; int entries = covalentMapData[covalentIndex++];
std::vector< int > covalentMap(entries); std::vector< int > covalentMap(entries);
for( unsigned int ii = 0; ii < entries; ii++ ){ for (unsigned int ii = 0; ii < entries; ii++) {
covalentMap[ii] = covalentMapData[covalentIndex++]; covalentMap[ii] = covalentMapData[covalentIndex++];
} }
amoebaMultipoleForce->setCovalentMap( particleIndex, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(typeIndex), covalentMap ); amoebaMultipoleForce->setCovalentMap(particleIndex, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(typeIndex), covalentMap);
} }
system.addForce(amoebaMultipoleForce); system.addForce(amoebaMultipoleForce);
   
...@@ -5742,9 +5742,9 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari ...@@ -5742,9 +5742,9 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari
// GK force // GK force
   
AmoebaGeneralizedKirkwoodForce* amoebaGeneralizedKirkwoodForce = new AmoebaGeneralizedKirkwoodForce(); AmoebaGeneralizedKirkwoodForce* amoebaGeneralizedKirkwoodForce = new AmoebaGeneralizedKirkwoodForce();
amoebaGeneralizedKirkwoodForce->setSolventDielectric( 7.8300000e+01 ); amoebaGeneralizedKirkwoodForce->setSolventDielectric( 7.8300000e+01);
amoebaGeneralizedKirkwoodForce->setSoluteDielectric( 1.0000000e+00 ); amoebaGeneralizedKirkwoodForce->setSoluteDielectric( 1.0000000e+00);
amoebaGeneralizedKirkwoodForce->setIncludeCavityTerm( includeCavityTerm ); amoebaGeneralizedKirkwoodForce->setIncludeCavityTerm(includeCavityTerm);
   
// addParticle: charge, radius, scalingFactor // addParticle: charge, radius, scalingFactor
   
...@@ -6348,8 +6348,8 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari ...@@ -6348,8 +6348,8 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari
}; };
   
const double* gkData = generalizedKirkwoodData; const double* gkData = generalizedKirkwoodData;
for( unsigned int ii = 0; ii < numberOfParticles; ii++ ){ for (unsigned int ii = 0; ii < numberOfParticles; ii++) {
amoebaGeneralizedKirkwoodForce->addParticle( gkData[0], gkData[1], gkData[2] ); amoebaGeneralizedKirkwoodForce->addParticle(gkData[0], gkData[1], gkData[2]);
gkData += 3; gkData += 3;
} }
system.addForce(amoebaGeneralizedKirkwoodForce); system.addForce(amoebaGeneralizedKirkwoodForce);
...@@ -6958,15 +6958,15 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari ...@@ -6958,15 +6958,15 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari
}; };
   
const double* positionsDataPtr = positionsData; const double* positionsDataPtr = positionsData;
for( unsigned int ii = 0; ii < numberOfParticles; ii++ ){ for (unsigned int ii = 0; ii < numberOfParticles; ii++) {
positions[ii] = Vec3( positionsDataPtr[0], positionsDataPtr[1], positionsDataPtr[2] ); positions[ii] = Vec3(positionsDataPtr[0], positionsDataPtr[1], positionsDataPtr[2]);
positionsDataPtr += 3; positionsDataPtr += 3;
} }
   
std::string platformName; std::string platformName;
platformName = "CUDA"; platformName = "CUDA";
LangevinIntegrator integrator(0.0, 0.1, 0.01); LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName( platformName ) ); Context context(system, integrator, Platform::getPlatformByName(platformName));
   
context.setPositions(positions); context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy); State state = context.getState(State::Forces | State::Energy);
...@@ -6976,117 +6976,39 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari ...@@ -6976,117 +6976,39 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari
   
// compare forces and energies // compare forces and energies
   
static void compareForcesEnergy( std::string& testName, double expectedEnergy, double energy, static void compareForcesEnergy(std::string& testName, double expectedEnergy, double energy,
const std::vector<Vec3>& expectedForces, const std::vector<Vec3>& expectedForces,
const std::vector<Vec3>& forces, double tolerance, FILE* log ) { const std::vector<Vec3>& forces, double tolerance) {
for (unsigned int ii = 0; ii < forces.size(); ii++) {
//#define AMOEBA_DEBUG ASSERT_EQUAL_VEC_MOD(expectedForces[ii], forces[ii], tolerance, testName);
#ifdef AMOEBA_DEBUG
if( log ){
double conversion = 1.0/4.184;
double energyAbsDiff = fabs( expectedEnergy - energy );
double energyRelDiff = 2.0*energyAbsDiff/( fabs( expectedEnergy ) + fabs( energy ) + 1.0e-08 );
(void) fprintf( log, "%s: expected energy=%14.7e %14.7e absDiff=%15.7e relDiff=%15.7e\n", testName.c_str(), conversion*expectedEnergy, conversion*energy,
conversion*energyAbsDiff, conversion*energyRelDiff );
if( conversion != 1.0 )conversion *= -0.1;
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
double expectedNorm = sqrt( expectedForces[ii][0]*expectedForces[ii][0] +
expectedForces[ii][1]*expectedForces[ii][1] +
expectedForces[ii][2]*expectedForces[ii][2] );
double norm = sqrt( forces[ii][0]*forces[ii][0] + forces[ii][1]*forces[ii][1] + forces[ii][2]*forces[ii][2] );
double absDiff = fabs( norm - expectedNorm );
double relDiff = 2.0*absDiff/(fabs( norm ) + fabs( expectedNorm ) + 1.0e-08);
(void) fprintf( log, "%6u %15.7e %15.7e [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
conversion*absDiff, conversion*relDiff,
conversion*expectedForces[ii][0], conversion*expectedForces[ii][1], conversion*expectedForces[ii][2],
conversion*forces[ii][0], conversion*forces[ii][1], conversion*forces[ii][2], conversion*expectedNorm, conversion*norm );
}
(void) fflush( log );
conversion = 1.0;
(void) fprintf( log, "\n%s: expected energy=%14.7e %14.7e no conversion\n", testName.c_str(), conversion*expectedEnergy, conversion*energy );
if( conversion != 1.0 )conversion = -1.0;
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
conversion*expectedForces[ii][0], conversion*expectedForces[ii][1], conversion*expectedForces[ii][2],
conversion*forces[ii][0], conversion*forces[ii][1], conversion*forces[ii][2] );
}
(void) fflush( log );
} }
#endif ASSERT_EQUAL_TOL_MOD(expectedEnergy, energy, tolerance, testName);
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
ASSERT_EQUAL_VEC_MOD( expectedForces[ii], forces[ii], tolerance, testName );
}
ASSERT_EQUAL_TOL_MOD( expectedEnergy, energy, tolerance, testName );
} }
   
// compare relative differences in force norms and energies // compare relative differences in force norms and energies
   
static void compareForceNormsEnergy( std::string& testName, double expectedEnergy, double energy, static void compareForceNormsEnergy(std::string& testName, double expectedEnergy, double energy,
std::vector<Vec3>& expectedForces, std::vector<Vec3>& expectedForces,
const std::vector<Vec3>& forces, double tolerance, FILE* log ) { const std::vector<Vec3>& forces, double tolerance) {
for (unsigned int ii = 0; ii < forces.size(); ii++) {
double expectedNorm = sqrt(expectedForces[ii][0]*expectedForces[ii][0] +
//#define AMOEBA_DEBUG expectedForces[ii][1]*expectedForces[ii][1] +
#ifdef AMOEBA_DEBUG expectedForces[ii][2]*expectedForces[ii][2]);
if( log ){
double conversion = 1.0/4.184;
double energyAbsDiff = fabs( expectedEnergy - energy );
double energyRelDiff = 2.0*energyAbsDiff/( fabs( expectedEnergy ) + fabs( energy ) + 1.0e-08 );
(void) fprintf( log, "%s: expected energy=%14.7e %14.7e absDiff=%15.7e relDiff=%15.7e\n", testName.c_str(), conversion*expectedEnergy, conversion*energy,
conversion*energyAbsDiff, conversion*energyRelDiff );
if( conversion != 1.0 )conversion *= -0.1;
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
double expectedNorm = sqrt( expectedForces[ii][0]*expectedForces[ii][0] +
expectedForces[ii][1]*expectedForces[ii][1] +
expectedForces[ii][2]*expectedForces[ii][2] );
double norm = sqrt( forces[ii][0]*forces[ii][0] + forces[ii][1]*forces[ii][1] + forces[ii][2]*forces[ii][2] );
double absDiff = fabs( (norm - expectedNorm) );
double relDiff = 2.0*absDiff/(fabs( norm ) + fabs( expectedNorm ) + 1.0e-08);
(void) fprintf( log, "%6u %15.7e %15.7e [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e] %15.7e %15.7e\n", ii,
fabs(conversion)*absDiff, relDiff,
conversion*expectedForces[ii][0], conversion*expectedForces[ii][1], conversion*expectedForces[ii][2],
conversion*forces[ii][0], conversion*forces[ii][1], conversion*forces[ii][2],
fabs(conversion)*expectedNorm, fabs(conversion)*norm );
}
(void) fflush( log );
conversion = 1.0;
(void) fprintf( log, "\n%s: expected energy=%14.7e %14.7e no conversion\n", testName.c_str(), conversion*expectedEnergy, conversion*energy );
if( conversion != 1.0 )conversion = -1.0;
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
conversion*expectedForces[ii][0], conversion*expectedForces[ii][1], conversion*expectedForces[ii][2],
conversion*forces[ii][0], conversion*forces[ii][1], conversion*forces[ii][2] );
}
(void) fflush( log );
}
#endif
   
for( unsigned int ii = 0; ii < forces.size(); ii++ ){ double norm = sqrt(forces[ii][0]*forces[ii][0] + forces[ii][1]*forces[ii][1] + forces[ii][2]*forces[ii][2]);
double expectedNorm = sqrt( expectedForces[ii][0]*expectedForces[ii][0] + double absDiff = fabs(norm - expectedNorm);
expectedForces[ii][1]*expectedForces[ii][1] + double relDiff = 2.0*absDiff/(fabs(norm) + fabs(expectedNorm) + 1.0e-08);
expectedForces[ii][2]*expectedForces[ii][2] );
   
double norm = sqrt( forces[ii][0]*forces[ii][0] + forces[ii][1]*forces[ii][1] + forces[ii][2]*forces[ii][2] ); if (relDiff > tolerance && absDiff > 0.001) {
double absDiff = fabs( norm - expectedNorm );
double relDiff = 2.0*absDiff/(fabs( norm ) + fabs( expectedNorm ) + 1.0e-08);
if( relDiff > tolerance && absDiff > 0.001 ){
std::stringstream details; std::stringstream details;
details << testName << "Relative difference in norms " << relDiff << " larger than allowed tolerance at particle=" << ii; details << testName << "Relative difference in norms " << relDiff << " larger than allowed tolerance at particle=" << ii;
details << ": norms=" << norm << " expected norm=" << expectedNorm; details << ": norms=" << norm << " expected norm=" << expectedNorm;
throwException(__FILE__, __LINE__, details.str()); throwException(__FILE__, __LINE__, details.str());
} }
} }
double energyAbsDiff = fabs( expectedEnergy - energy ); double energyAbsDiff = fabs(expectedEnergy - energy);
double energyRelDiff = 2.0*energyAbsDiff/( fabs( expectedEnergy ) + fabs( energy ) + 1.0e-08 ); double energyRelDiff = 2.0*energyAbsDiff/(fabs(expectedEnergy) + fabs(energy) + 1.0e-08);
if( energyRelDiff > tolerance ){ if (energyRelDiff > tolerance) {
std::stringstream details; std::stringstream details;
details << testName << "Relative difference in energies " << energyRelDiff << " larger than allowed tolerance."; details << testName << "Relative difference in energies " << energyRelDiff << " larger than allowed tolerance.";
details << "Energies=" << energy << " expected energy=" << expectedEnergy; details << "Energies=" << energy << " expected energy=" << expectedEnergy;
...@@ -7096,7 +7018,7 @@ static void compareForceNormsEnergy( std::string& testName, double expectedEnerg ...@@ -7096,7 +7018,7 @@ static void compareForceNormsEnergy( std::string& testName, double expectedEnerg
   
// test GK direct polarization for system comprised of two ammonia molecules // test GK direct polarization for system comprised of two ammonia molecules
   
static void testGeneralizedKirkwoodAmmoniaDirectPolarization( FILE* log ) { static void testGeneralizedKirkwoodAmmoniaDirectPolarization() {
   
std::string testName = "testGeneralizedKirkwoodAmmoniaDirectPolarization"; std::string testName = "testGeneralizedKirkwoodAmmoniaDirectPolarization";
   
...@@ -7109,27 +7031,27 @@ static void testGeneralizedKirkwoodAmmoniaDirectPolarization( FILE* log ) { ...@@ -7109,27 +7031,27 @@ static void testGeneralizedKirkwoodAmmoniaDirectPolarization( FILE* log ) {
setupMultipoleAmmonia(system, amoebaGeneralizedKirkwoodForce, AmoebaMultipoleForce::Direct, 0); setupMultipoleAmmonia(system, amoebaGeneralizedKirkwoodForce, AmoebaMultipoleForce::Direct, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01); LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName("CUDA")); Context context(system, integrator, Platform::getPlatformByName("CUDA"));
getForcesEnergyMultipoleAmmonia(context, forces, energy, log ); getForcesEnergyMultipoleAmmonia(context, forces, energy);
std::vector<Vec3> expectedForces(numberOfParticles); std::vector<Vec3> expectedForces(numberOfParticles);
   
double expectedEnergy = -7.6636680e+01; double expectedEnergy = -7.6636680e+01;
   
expectedForces[0] = Vec3( -6.9252994e+02, -8.9085133e+00, 9.6489739e+01 ); expectedForces[0] = Vec3( -6.9252994e+02, -8.9085133e+00, 9.6489739e+01);
expectedForces[1] = Vec3( 1.5593797e+02, -6.0331931e+01, 1.5104507e+01 ); expectedForces[1] = Vec3( 1.5593797e+02, -6.0331931e+01, 1.5104507e+01);
expectedForces[2] = Vec3( 1.5870088e+02, 6.1702809e+01, 6.7708985e+00 ); expectedForces[2] = Vec3( 1.5870088e+02, 6.1702809e+01, 6.7708985e+00);
expectedForces[3] = Vec3( 1.4089885e+02, 7.5870617e+00, -1.1362294e+02 ); expectedForces[3] = Vec3( 1.4089885e+02, 7.5870617e+00, -1.1362294e+02);
expectedForces[4] = Vec3( -1.8916205e+02, 2.1465549e-01, -4.3433152e+02 ); expectedForces[4] = Vec3( -1.8916205e+02, 2.1465549e-01, -4.3433152e+02);
expectedForces[5] = Vec3( 1.0208290e+01, 6.2676753e+01, 1.4987953e+02 ); expectedForces[5] = Vec3( 1.0208290e+01, 6.2676753e+01, 1.4987953e+02);
expectedForces[6] = Vec3( 4.0621859e+02, 1.8962203e-01, 1.3021956e+02 ); expectedForces[6] = Vec3( 4.0621859e+02, 1.8962203e-01, 1.3021956e+02);
expectedForces[7] = Vec3( 9.7274235e+00, -6.3130458e+01, 1.4949024e+02 ); expectedForces[7] = Vec3( 9.7274235e+00, -6.3130458e+01, 1.4949024e+02);
   
double tolerance = 1.0e-04; double tolerance = 1.0e-04;
compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log ); compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance);
} }
   
// test GK mutual polarization for system comprised of two ammonia molecules // test GK mutual polarization for system comprised of two ammonia molecules
   
static void testGeneralizedKirkwoodAmmoniaMutualPolarization( FILE* log ) { static void testGeneralizedKirkwoodAmmoniaMutualPolarization() {
   
std::string testName = "testGeneralizedKirkwoodAmmoniaMutualPolarization"; std::string testName = "testGeneralizedKirkwoodAmmoniaMutualPolarization";
   
...@@ -7142,28 +7064,28 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarization( FILE* log ) { ...@@ -7142,28 +7064,28 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarization( FILE* log ) {
setupMultipoleAmmonia(system, amoebaGeneralizedKirkwoodForce, AmoebaMultipoleForce::Mutual, 0); setupMultipoleAmmonia(system, amoebaGeneralizedKirkwoodForce, AmoebaMultipoleForce::Mutual, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01); LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName("CUDA")); Context context(system, integrator, Platform::getPlatformByName("CUDA"));
getForcesEnergyMultipoleAmmonia(context, forces, energy, log ); getForcesEnergyMultipoleAmmonia(context, forces, energy);
std::vector<Vec3> expectedForces(numberOfParticles); std::vector<Vec3> expectedForces(numberOfParticles);
   
double expectedEnergy = -7.8018875e+01; double expectedEnergy = -7.8018875e+01;
   
expectedForces[0] = Vec3( -7.6820301e+02, -1.0102760e+01, 1.0094389e+02 ); expectedForces[0] = Vec3( -7.6820301e+02, -1.0102760e+01, 1.0094389e+02);
expectedForces[1] = Vec3( 1.7037307e+02, -7.5621857e+01, 2.3320365e+01 ); expectedForces[1] = Vec3( 1.7037307e+02, -7.5621857e+01, 2.3320365e+01);
expectedForces[2] = Vec3( 1.7353828e+02, 7.7199741e+01, 1.3965379e+01 ); expectedForces[2] = Vec3( 1.7353828e+02, 7.7199741e+01, 1.3965379e+01);
expectedForces[3] = Vec3( 1.5045244e+02, 8.5784569e+00, -1.3377619e+02 ); expectedForces[3] = Vec3( 1.5045244e+02, 8.5784569e+00, -1.3377619e+02);
expectedForces[4] = Vec3( -2.1811615e+02, -1.6818022e-01, -4.6103163e+02 ); expectedForces[4] = Vec3( -2.1811615e+02, -1.6818022e-01, -4.6103163e+02);
expectedForces[5] = Vec3( 6.2091942e+00, 7.6748687e+01, 1.5883463e+02 ); expectedForces[5] = Vec3( 6.2091942e+00, 7.6748687e+01, 1.5883463e+02);
expectedForces[6] = Vec3( 4.8035662e+02, 4.9704902e-01, 1.3948083e+02 ); expectedForces[6] = Vec3( 4.8035662e+02, 4.9704902e-01, 1.3948083e+02);
expectedForces[7] = Vec3( 5.3895456e+00, -7.7131137e+01, 1.5826273e+02 ); expectedForces[7] = Vec3( 5.3895456e+00, -7.7131137e+01, 1.5826273e+02);
   
double tolerance = 2.0e-04; double tolerance = 2.0e-04;
compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log ); compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance);
} }
   
// test GK mutual polarization for system comprised of two ammonia molecules // test GK mutual polarization for system comprised of two ammonia molecules
// including cavity term // including cavity term
   
static void testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm( FILE* log ) { static void testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm() {
   
std::string testName = "testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm"; std::string testName = "testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm";
   
...@@ -7176,22 +7098,22 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm( FILE ...@@ -7176,22 +7098,22 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm( FILE
setupMultipoleAmmonia(system, amoebaGeneralizedKirkwoodForce, AmoebaMultipoleForce::Mutual, 1); setupMultipoleAmmonia(system, amoebaGeneralizedKirkwoodForce, AmoebaMultipoleForce::Mutual, 1);
LangevinIntegrator integrator(0.0, 0.1, 0.01); LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName("CUDA")); Context context(system, integrator, Platform::getPlatformByName("CUDA"));
getForcesEnergyMultipoleAmmonia(context, forces, energy, log ); getForcesEnergyMultipoleAmmonia(context, forces, energy);
std::vector<Vec3> expectedForces(numberOfParticles); std::vector<Vec3> expectedForces(numberOfParticles);
   
double expectedEnergy = -6.0434582e+01; double expectedEnergy = -6.0434582e+01;
   
expectedForces[0] = Vec3( -7.8323218e+02, -1.0097644e+01, 1.0256890e+02 ); expectedForces[0] = Vec3( -7.8323218e+02, -1.0097644e+01, 1.0256890e+02);
expectedForces[1] = Vec3( 1.7078480e+02, -7.1896701e+01, 2.0840172e+01 ); expectedForces[1] = Vec3( 1.7078480e+02, -7.1896701e+01, 2.0840172e+01);
expectedForces[2] = Vec3( 1.7394089e+02, 7.3488594e+01, 1.1484648e+01 ); expectedForces[2] = Vec3( 1.7394089e+02, 7.3488594e+01, 1.1484648e+01);
expectedForces[3] = Vec3( 1.5169364e+02, 8.5611299e+00, -1.2968050e+02 ); expectedForces[3] = Vec3( 1.5169364e+02, 8.5611299e+00, -1.2968050e+02);
expectedForces[4] = Vec3( -2.1669693e+02, -1.5926823e-01, -4.6636274e+02 ); expectedForces[4] = Vec3( -2.1669693e+02, -1.5926823e-01, -4.6636274e+02);
expectedForces[5] = Vec3( 8.7397444e+00, 7.3330990e+01, 1.6016898e+02 ); expectedForces[5] = Vec3( 8.7397444e+00, 7.3330990e+01, 1.6016898e+02);
expectedForces[6] = Vec3( 4.8684950e+02, 4.8937161e-01, 1.4137061e+02 ); expectedForces[6] = Vec3( 4.8684950e+02, 4.8937161e-01, 1.4137061e+02);
expectedForces[7] = Vec3( 7.9205382e+00, -7.3716473e+01, 1.5960993e+02 ); expectedForces[7] = Vec3( 7.9205382e+00, -7.3716473e+01, 1.5960993e+02);
   
double tolerance = 1.0e-04; double tolerance = 1.0e-04;
compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log ); compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance);
// Try changing the particle parameters and make sure it's still correct. // Try changing the particle parameters and make sure it's still correct.
...@@ -7208,7 +7130,7 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm( FILE ...@@ -7208,7 +7130,7 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm( FILE
bool exceptionThrown = false; bool exceptionThrown = false;
try { try {
// This should throw an exception. // This should throw an exception.
compareForcesEnergy(testName, state2.getPotentialEnergy(), state1.getPotentialEnergy(), state2.getForces(), state1.getForces(), tolerance, log); compareForcesEnergy(testName, state2.getPotentialEnergy(), state1.getPotentialEnergy(), state2.getForces(), state1.getForces(), tolerance);
} }
catch (std::exception ex) { catch (std::exception ex) {
exceptionThrown = true; exceptionThrown = true;
...@@ -7216,12 +7138,12 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm( FILE ...@@ -7216,12 +7138,12 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm( FILE
ASSERT(exceptionThrown); ASSERT(exceptionThrown);
amoebaGeneralizedKirkwoodForce->updateParametersInContext(context); amoebaGeneralizedKirkwoodForce->updateParametersInContext(context);
state1 = context.getState(State::Forces | State::Energy); state1 = context.getState(State::Forces | State::Energy);
compareForcesEnergy(testName, state2.getPotentialEnergy(), state1.getPotentialEnergy(), state2.getForces(), state1.getForces(), tolerance, log); compareForcesEnergy(testName, state2.getPotentialEnergy(), state1.getPotentialEnergy(), state2.getForces(), state1.getForces(), tolerance);
} }
   
// test GK direct polarization for villin system // test GK direct polarization for villin system
   
static void testGeneralizedKirkwoodVillinDirectPolarization( FILE* log ) { static void testGeneralizedKirkwoodVillinDirectPolarization() {
   
std::string testName = "testGeneralizedKirkwoodVillinDirectPolarization"; std::string testName = "testGeneralizedKirkwoodVillinDirectPolarization";
   
...@@ -7229,7 +7151,7 @@ static void testGeneralizedKirkwoodVillinDirectPolarization( FILE* log ) { ...@@ -7229,7 +7151,7 @@ static void testGeneralizedKirkwoodVillinDirectPolarization( FILE* log ) {
std::vector<Vec3> forces; std::vector<Vec3> forces;
double energy; double energy;
   
setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Direct, 0, forces, energy, log ); setupAndGetForcesEnergyMultipoleVillin(AmoebaMultipoleForce::Direct, 0, forces, energy);
std::vector<Vec3> expectedForces(numberOfParticles); std::vector<Vec3> expectedForces(numberOfParticles);
   
double expectedEnergy = -8.4281157e+03; double expectedEnergy = -8.4281157e+03;
...@@ -7834,18 +7756,18 @@ static void testGeneralizedKirkwoodVillinDirectPolarization( FILE* log ) { ...@@ -7834,18 +7756,18 @@ static void testGeneralizedKirkwoodVillinDirectPolarization( FILE* log ) {
}; };
   
const double* forceDataPtr = forceData; const double* forceDataPtr = forceData;
for( unsigned int ii = 0; ii < numberOfParticles; ii++ ){ for (unsigned int ii = 0; ii < numberOfParticles; ii++) {
expectedForces[ii] = Vec3( forceDataPtr[0], forceDataPtr[1], forceDataPtr[2] ); expectedForces[ii] = Vec3(forceDataPtr[0], forceDataPtr[1], forceDataPtr[2]);
forceDataPtr += 3; forceDataPtr += 3;
} }
   
double tolerance = 1.0e-05; double tolerance = 1.0e-05;
compareForceNormsEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log ); compareForceNormsEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance);
} }
   
// test GK mutual polarization for villin system // test GK mutual polarization for villin system
   
static void testGeneralizedKirkwoodVillinMutualPolarization( FILE* log ) { static void testGeneralizedKirkwoodVillinMutualPolarization() {
   
std::string testName = "testGeneralizedKirkwoodVillinMutualPolarization"; std::string testName = "testGeneralizedKirkwoodVillinMutualPolarization";
   
...@@ -7853,7 +7775,7 @@ static void testGeneralizedKirkwoodVillinMutualPolarization( FILE* log ) { ...@@ -7853,7 +7775,7 @@ static void testGeneralizedKirkwoodVillinMutualPolarization( FILE* log ) {
std::vector<Vec3> forces; std::vector<Vec3> forces;
double energy; double energy;
   
setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Mutual, 0, forces, energy, log ); setupAndGetForcesEnergyMultipoleVillin(AmoebaMultipoleForce::Mutual, 0, forces, energy);
std::vector<Vec3> expectedForces(numberOfParticles); std::vector<Vec3> expectedForces(numberOfParticles);
   
double expectedEnergy = -8.6477811e+03; double expectedEnergy = -8.6477811e+03;
...@@ -8458,13 +8380,13 @@ static void testGeneralizedKirkwoodVillinMutualPolarization( FILE* log ) { ...@@ -8458,13 +8380,13 @@ static void testGeneralizedKirkwoodVillinMutualPolarization( FILE* log ) {
}; };
   
const double* forceDataPtr = forceData; const double* forceDataPtr = forceData;
for( unsigned int ii = 0; ii < numberOfParticles; ii++ ){ for (unsigned int ii = 0; ii < numberOfParticles; ii++) {
expectedForces[ii] = Vec3( forceDataPtr[0], forceDataPtr[1], forceDataPtr[2] ); expectedForces[ii] = Vec3(forceDataPtr[0], forceDataPtr[1], forceDataPtr[2]);
forceDataPtr += 3; forceDataPtr += 3;
} }
   
double tolerance = 1.0e-05; double tolerance = 1.0e-05;
compareForceNormsEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log ); compareForceNormsEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance);
} }
   
int main(int argc, char* argv[]) { int main(int argc, char* argv[]) {
...@@ -8474,17 +8396,15 @@ int main(int argc, char* argv[]) { ...@@ -8474,17 +8396,15 @@ int main(int argc, char* argv[]) {
if (argc > 1) if (argc > 1)
Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1])); Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1]));
   
FILE* log = NULL;
// test direct and mutual polarization cases and // test direct and mutual polarization cases and
// mutual polarization w/ the cavity term // mutual polarization w/ the cavity term
   
testGeneralizedKirkwoodAmmoniaDirectPolarization( log ); testGeneralizedKirkwoodAmmoniaDirectPolarization();
testGeneralizedKirkwoodAmmoniaMutualPolarization( log ); testGeneralizedKirkwoodAmmoniaMutualPolarization();
testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm( log ); testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm();
   
testGeneralizedKirkwoodVillinDirectPolarization( log ); testGeneralizedKirkwoodVillinDirectPolarization();
testGeneralizedKirkwoodVillinMutualPolarization( log ); testGeneralizedKirkwoodVillinMutualPolarization();
   
} catch(const std::exception& e) { } catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl; std::cout << "exception: " << e.what() << std::endl;
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaInPlaneAngleForce.cpp
View file @ f11bafd9
...@@ -63,7 +63,7 @@ const double TOL = 1e-5; ...@@ -63,7 +63,7 @@ const double TOL = 1e-5;
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
static void crossProductVector3( double* vectorX, double* vectorY, double* vectorZ ){ static void crossProductVector3(double* vectorX, double* vectorY, double* vectorZ) {
vectorZ[0] = vectorX[1]*vectorY[2] - vectorX[2]*vectorY[1]; vectorZ[0] = vectorX[1]*vectorY[2] - vectorX[2]*vectorY[1];
vectorZ[1] = vectorX[2]*vectorY[0] - vectorX[0]*vectorY[2]; vectorZ[1] = vectorX[2]*vectorY[0] - vectorX[0]*vectorY[2];
...@@ -72,30 +72,23 @@ static void crossProductVector3( double* vectorX, double* vectorY, double* vecto ...@@ -72,30 +72,23 @@ static void crossProductVector3( double* vectorX, double* vectorY, double* vecto
return; return;
} }
static double dotVector3( double* vectorX, double* vectorY ){ static double dotVector3(double* vectorX, double* vectorY) {
return vectorX[0]*vectorY[0] + vectorX[1]*vectorY[1] + vectorX[2]*vectorY[2]; return vectorX[0]*vectorY[0] + vectorX[1]*vectorY[1] + vectorX[2]*vectorY[2];
} }
static void getPrefactorsGivenInPlaneAngleCosine( double cosine, double idealInPlaneAngle, double quadraticK, double cubicK, static void getPrefactorsGivenInPlaneAngleCosine(double cosine, double idealInPlaneAngle, double quadraticK, double cubicK,
double quarticK, double penticK, double sexticK, double quarticK, double penticK, double sexticK,
double* dEdR, double* energyTerm, FILE* log ) { double* dEdR, double* energyTerm) {
double angle; double angle;
if( cosine >= 1.0 ){ if (cosine >= 1.0) {
angle = 0.0f; angle = 0.0f;
} else if( cosine <= -1.0 ){ } else if (cosine <= -1.0) {
angle = RADIAN*PI_M; angle = RADIAN*PI_M;
} else { } else {
angle = RADIAN*acos(cosine); angle = RADIAN*acos(cosine);
} }
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "getPrefactorsGivenInPlaneAngleCosine: cosine=%10.3e angle=%10.3e ideal=%10.3e\n", cosine, angle, idealInPlaneAngle );
(void) fflush( log );
}
#endif
double deltaIdeal = angle - idealInPlaneAngle; double deltaIdeal = angle - idealInPlaneAngle;
double deltaIdeal2 = deltaIdeal*deltaIdeal; double deltaIdeal2 = deltaIdeal*deltaIdeal;
double deltaIdeal3 = deltaIdeal*deltaIdeal2; double deltaIdeal3 = deltaIdeal*deltaIdeal2;
...@@ -103,11 +96,11 @@ static void getPrefactorsGivenInPlaneAngleCosine( double cosine, double idealInP ...@@ -103,11 +96,11 @@ static void getPrefactorsGivenInPlaneAngleCosine( double cosine, double idealInP
// deltaIdeal = r - r_0 // deltaIdeal = r - r_0
*dEdR = ( 2.0 + *dEdR = (2.0 +
3.0*cubicK* deltaIdeal + 3.0*cubicK* deltaIdeal +
4.0*quarticK*deltaIdeal2 + 4.0*quarticK*deltaIdeal2 +
5.0*penticK* deltaIdeal3 + 5.0*penticK* deltaIdeal3 +
6.0*sexticK* deltaIdeal4 ); 6.0*sexticK* deltaIdeal4 );
*dEdR *= RADIAN*quadraticK*deltaIdeal; *dEdR *= RADIAN*quadraticK*deltaIdeal;
...@@ -122,26 +115,18 @@ static void getPrefactorsGivenInPlaneAngleCosine( double cosine, double idealInP ...@@ -122,26 +115,18 @@ static void getPrefactorsGivenInPlaneAngleCosine( double cosine, double idealInP
} }
static void computeAmoebaInPlaneAngleForce(int bondIndex, std::vector<Vec3>& positions, AmoebaInPlaneAngleForce& amoebaInPlaneAngleForce, static void computeAmoebaInPlaneAngleForce(int bondIndex, std::vector<Vec3>& positions, AmoebaInPlaneAngleForce& amoebaInPlaneAngleForce,
std::vector<Vec3>& forces, double* energy, FILE* log ) { std::vector<Vec3>& forces, double* energy) {
int particle1, particle2, particle3, particle4; int particle1, particle2, particle3, particle4;
double idealInPlaneAngle; double idealInPlaneAngle;
double quadraticK; double quadraticK;
amoebaInPlaneAngleForce.getAngleParameters(bondIndex, particle1, particle2, particle3, particle4, idealInPlaneAngle, quadraticK ); amoebaInPlaneAngleForce.getAngleParameters(bondIndex, particle1, particle2, particle3, particle4, idealInPlaneAngle, quadraticK);
double cubicK = amoebaInPlaneAngleForce.getAmoebaGlobalInPlaneAngleCubic(); double cubicK = amoebaInPlaneAngleForce.getAmoebaGlobalInPlaneAngleCubic();
double quarticK = amoebaInPlaneAngleForce.getAmoebaGlobalInPlaneAngleQuartic(); double quarticK = amoebaInPlaneAngleForce.getAmoebaGlobalInPlaneAngleQuartic();
double penticK = amoebaInPlaneAngleForce.getAmoebaGlobalInPlaneAnglePentic(); double penticK = amoebaInPlaneAngleForce.getAmoebaGlobalInPlaneAnglePentic();
double sexticK = amoebaInPlaneAngleForce.getAmoebaGlobalInPlaneAngleSextic(); double sexticK = amoebaInPlaneAngleForce.getAmoebaGlobalInPlaneAngleSextic();
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaInPlaneAngleForce: bond %d [%d %d %d %d] ang=%10.3f k=%10.3f [%10.3e %10.3e %10.3e %10.3e]\n",
bondIndex, particle1, particle2, particle3, particle4, idealInPlaneAngle, quadraticK, cubicK, quarticK, penticK, sexticK );
(void) fflush( log );
}
#endif
// T = AD x CD // T = AD x CD
// P = B + T*delta // P = B + T*delta
// AP = A - P // AP = A - P
...@@ -161,83 +146,75 @@ static void computeAmoebaInPlaneAngleForce(int bondIndex, std::vector<Vec3>& po ...@@ -161,83 +146,75 @@ static void computeAmoebaInPlaneAngleForce(int bondIndex, std::vector<Vec3>& po
// APxM, CPxM, ADxBD, BDxCD, TxCD, ADxT, dBxAD, CDxdB, LastDeltaAtomIndex // APxM, CPxM, ADxBD, BDxCD, TxCD, ADxT, dBxAD, CDxdB, LastDeltaAtomIndex
double deltaR[LastDeltaAtomIndex][3]; double deltaR[LastDeltaAtomIndex][3];
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
deltaR[AD][ii] = positions[particle1][ii] - positions[particle4][ii]; deltaR[AD][ii] = positions[particle1][ii] - positions[particle4][ii];
deltaR[BD][ii] = positions[particle2][ii] - positions[particle4][ii]; deltaR[BD][ii] = positions[particle2][ii] - positions[particle4][ii];
deltaR[CD][ii] = positions[particle3][ii] - positions[particle4][ii]; deltaR[CD][ii] = positions[particle3][ii] - positions[particle4][ii];
} }
crossProductVector3( deltaR[AD], deltaR[CD], deltaR[T] ); crossProductVector3(deltaR[AD], deltaR[CD], deltaR[T]);
double rT2 = dotVector3( deltaR[T], deltaR[T] ); double rT2 = dotVector3(deltaR[T], deltaR[T]);
double delta = dotVector3( deltaR[T], deltaR[BD] ); double delta = dotVector3(deltaR[T], deltaR[BD]);
delta *= -1.0/rT2; delta *= -1.0/rT2;
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
deltaR[P][ii] = positions[particle2][ii] + deltaR[T][ii]*delta; deltaR[P][ii] = positions[particle2][ii] + deltaR[T][ii]*delta;
deltaR[AP][ii] = positions[particle1][ii] - deltaR[P][ii]; deltaR[AP][ii] = positions[particle1][ii] - deltaR[P][ii];
deltaR[CP][ii] = positions[particle3][ii] - deltaR[P][ii]; deltaR[CP][ii] = positions[particle3][ii] - deltaR[P][ii];
} }
double rAp2 = dotVector3( deltaR[AP], deltaR[AP] ); double rAp2 = dotVector3(deltaR[AP], deltaR[AP]);
double rCp2 = dotVector3( deltaR[CP], deltaR[CP] ); double rCp2 = dotVector3(deltaR[CP], deltaR[CP]);
if( rAp2 <= 0.0 && rCp2 <= 0.0 ){ if (rAp2 <= 0.0 && rCp2 <= 0.0) {
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaInPlaneAngleForce: rAp2 or rCp2 <= 0.0\n" );
(void) fflush( log );
}
#endif
return;
} }
crossProductVector3( deltaR[CP], deltaR[AP], deltaR[M] ); crossProductVector3(deltaR[CP], deltaR[AP], deltaR[M]);
double rm = dotVector3( deltaR[M], deltaR[M] ); double rm = dotVector3(deltaR[M], deltaR[M]);
rm = sqrt( rm ); rm = sqrt(rm);
if( rm < 0.000001 ){ if (rm < 0.000001) {
rm = 0.000001; rm = 0.000001;
} }
double dot = dotVector3( deltaR[AP], deltaR[CP] ); double dot = dotVector3(deltaR[AP], deltaR[CP]);
double cosine = dot/sqrt( rAp2*rCp2 ); double cosine = dot/sqrt(rAp2*rCp2);
double dEdR; double dEdR;
double energyTerm; double energyTerm;
getPrefactorsGivenInPlaneAngleCosine( cosine, idealInPlaneAngle, quadraticK, cubicK, getPrefactorsGivenInPlaneAngleCosine(cosine, idealInPlaneAngle, quadraticK, cubicK,
quarticK, penticK, sexticK, &dEdR, &energyTerm, log ); quarticK, penticK, sexticK, &dEdR, &energyTerm);
double termA = -dEdR/(rAp2*rm); double termA = -dEdR/(rAp2*rm);
double termC = dEdR/(rCp2*rm); double termC = dEdR/(rCp2*rm);
crossProductVector3( deltaR[AP], deltaR[M], deltaR[APxM] ); crossProductVector3(deltaR[AP], deltaR[M], deltaR[APxM]);
crossProductVector3( deltaR[CP], deltaR[M], deltaR[CPxM] ); crossProductVector3(deltaR[CP], deltaR[M], deltaR[CPxM]);
// forces will be gathered here // forces will be gathered here
enum { dA, dB, dC, dD, LastDIndex }; enum { dA, dB, dC, dD, LastDIndex };
double forceTerm[LastDIndex][3]; double forceTerm[LastDIndex][3];
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
forceTerm[dA][ii] = deltaR[APxM][ii]*termA; forceTerm[dA][ii] = deltaR[APxM][ii]*termA;
forceTerm[dC][ii] = deltaR[CPxM][ii]*termC; forceTerm[dC][ii] = deltaR[CPxM][ii]*termC;
forceTerm[dB][ii] = -1.0*( forceTerm[dA][ii] + forceTerm[dC][ii] ); forceTerm[dB][ii] = -1.0*(forceTerm[dA][ii] + forceTerm[dC][ii]);
} }
double pTrT2 = dotVector3( forceTerm[dB], deltaR[T] ); double pTrT2 = dotVector3(forceTerm[dB], deltaR[T]);
pTrT2 /= rT2; pTrT2 /= rT2;
crossProductVector3( deltaR[CD], forceTerm[dB], deltaR[CDxdB] ); crossProductVector3(deltaR[CD], forceTerm[dB], deltaR[CDxdB]);
crossProductVector3( forceTerm[dB], deltaR[AD], deltaR[dBxAD] ); crossProductVector3(forceTerm[dB], deltaR[AD], deltaR[dBxAD]);
if( fabs( pTrT2 ) > 1.0e-08 ){ if (fabs(pTrT2) > 1.0e-08) {
double delta2 = delta*2.0; double delta2 = delta*2.0;
crossProductVector3( deltaR[BD], deltaR[CD], deltaR[BDxCD] ); crossProductVector3(deltaR[BD], deltaR[CD], deltaR[BDxCD]);
crossProductVector3( deltaR[T], deltaR[CD], deltaR[TxCD] ); crossProductVector3(deltaR[T], deltaR[CD], deltaR[TxCD] );
crossProductVector3( deltaR[AD], deltaR[BD], deltaR[ADxBD] ); crossProductVector3(deltaR[AD], deltaR[BD], deltaR[ADxBD]);
crossProductVector3( deltaR[AD], deltaR[T], deltaR[ADxT] ); crossProductVector3(deltaR[AD], deltaR[T], deltaR[ADxT] );
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
double term = deltaR[BDxCD][ii] + delta2*deltaR[TxCD][ii]; double term = deltaR[BDxCD][ii] + delta2*deltaR[TxCD][ii];
forceTerm[dA][ii] += delta*deltaR[CDxdB][ii] + term*pTrT2; forceTerm[dA][ii] += delta*deltaR[CDxdB][ii] + term*pTrT2;
...@@ -245,15 +222,15 @@ static void computeAmoebaInPlaneAngleForce(int bondIndex, std::vector<Vec3>& po ...@@ -245,15 +222,15 @@ static void computeAmoebaInPlaneAngleForce(int bondIndex, std::vector<Vec3>& po
term = deltaR[ADxBD][ii] + delta2*deltaR[ADxT][ii]; term = deltaR[ADxBD][ii] + delta2*deltaR[ADxT][ii];
forceTerm[dC][ii] += delta*deltaR[dBxAD][ii] + term*pTrT2; forceTerm[dC][ii] += delta*deltaR[dBxAD][ii] + term*pTrT2;
forceTerm[dD][ii] = -( forceTerm[dA][ii] + forceTerm[dB][ii] + forceTerm[dC][ii] ); forceTerm[dD][ii] = -(forceTerm[dA][ii] + forceTerm[dB][ii] + forceTerm[dC][ii]);
} }
} else { } else {
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
forceTerm[dA][ii] += delta*deltaR[CDxdB][ii]; forceTerm[dA][ii] += delta*deltaR[CDxdB][ii];
forceTerm[dC][ii] += delta*deltaR[dBxAD][ii]; forceTerm[dC][ii] += delta*deltaR[dBxAD][ii];
forceTerm[dD][ii] = -( forceTerm[dA][ii] + forceTerm[dB][ii] + forceTerm[dC][ii] ); forceTerm[dD][ii] = -(forceTerm[dA][ii] + forceTerm[dB][ii] + forceTerm[dC][ii]);
} }
} }
...@@ -279,71 +256,47 @@ static void computeAmoebaInPlaneAngleForce(int bondIndex, std::vector<Vec3>& po ...@@ -279,71 +256,47 @@ static void computeAmoebaInPlaneAngleForce(int bondIndex, std::vector<Vec3>& po
} }
static void computeAmoebaInPlaneAngleForces( Context& context, AmoebaInPlaneAngleForce& amoebaInPlaneAngleForce, static void computeAmoebaInPlaneAngleForces(Context& context, AmoebaInPlaneAngleForce& amoebaInPlaneAngleForce,
std::vector<Vec3>& expectedForces, double* expectedEnergy, FILE* log ) { std::vector<Vec3>& expectedForces, double* expectedEnergy) {
// get positions and zero forces // get positions and zero forces
State state = context.getState(State::Positions); State state = context.getState(State::Positions);
std::vector<Vec3> positions = state.getPositions(); std::vector<Vec3> positions = state.getPositions();
expectedForces.resize( positions.size() ); expectedForces.resize(positions.size());
for( unsigned int ii = 0; ii < expectedForces.size(); ii++ ){ for (unsigned int ii = 0; ii < expectedForces.size(); ii++) {
expectedForces[ii][0] = expectedForces[ii][1] = expectedForces[ii][2] = 0.0; expectedForces[ii][0] = expectedForces[ii][1] = expectedForces[ii][2] = 0.0;
} }
// calculates forces/energy // calculates forces/energy
*expectedEnergy = 0.0; *expectedEnergy = 0.0;
for( int ii = 0; ii < amoebaInPlaneAngleForce.getNumAngles(); ii++ ){ for (int ii = 0; ii < amoebaInPlaneAngleForce.getNumAngles(); ii++) {
computeAmoebaInPlaneAngleForce(ii, positions, amoebaInPlaneAngleForce, expectedForces, expectedEnergy, log ); computeAmoebaInPlaneAngleForce(ii, positions, amoebaInPlaneAngleForce, expectedForces, expectedEnergy);
}
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaInPlaneAngleForces: expected energy=%14.7e\n", *expectedEnergy );
for( unsigned int ii = 0; ii < positions.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e]\n", ii, expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2] );
}
(void) fflush( log );
} }
#endif
return;
} }
void compareWithExpectedForceAndEnergy( Context& context, AmoebaInPlaneAngleForce& amoebaInPlaneAngleForce, void compareWithExpectedForceAndEnergy(Context& context, AmoebaInPlaneAngleForce& amoebaInPlaneAngleForce,
double tolerance, const std::string& idString, FILE* log) { double tolerance, const std::string& idString) {
std::vector<Vec3> expectedForces; std::vector<Vec3> expectedForces;
double expectedEnergy; double expectedEnergy;
computeAmoebaInPlaneAngleForces( context, amoebaInPlaneAngleForce, expectedForces, &expectedEnergy, log ); computeAmoebaInPlaneAngleForces(context, amoebaInPlaneAngleForce, expectedForces, &expectedEnergy);
State state = context.getState(State::Forces | State::Energy); State state = context.getState(State::Forces | State::Energy);
const std::vector<Vec3> forces = state.getForces(); const std::vector<Vec3> forces = state.getForces();
for (unsigned int ii = 0; ii < forces.size(); ii++) {
#ifdef AMOEBA_DEBUG ASSERT_EQUAL_VEC(expectedForces[ii], forces[ii], tolerance);
if( log ){
(void) fprintf( log, "computeAmoebaInPlaneAngleForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy() );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2], forces[ii][0], forces[ii][1], forces[ii][2] );
}
(void) fflush( log );
}
#endif
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
ASSERT_EQUAL_VEC( expectedForces[ii], forces[ii], tolerance );
} }
ASSERT_EQUAL_TOL( expectedEnergy, state.getPotentialEnergy(), tolerance ); ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), tolerance);
} }
void testOneAngle( FILE* log ) { void testOneAngle() {
System system; System system;
int numberOfParticles = 4; int numberOfParticles = 4;
for( int ii = 0; ii < numberOfParticles; ii++ ){ for (int ii = 0; ii < numberOfParticles; ii++) {
system.addParticle(1.0); system.addParticle(1.0);
} }
...@@ -365,7 +318,7 @@ void testOneAngle( FILE* log ) { ...@@ -365,7 +318,7 @@ void testOneAngle( FILE* log ) {
amoebaInPlaneAngleForce->setAmoebaGlobalInPlaneAngleSextic(sexticK); amoebaInPlaneAngleForce->setAmoebaGlobalInPlaneAngleSextic(sexticK);
system.addForce(amoebaInPlaneAngleForce); system.addForce(amoebaInPlaneAngleForce);
Context context(system, integrator, Platform::getPlatformByName( "CUDA")); Context context(system, integrator, Platform::getPlatformByName("CUDA"));
std::vector<Vec3> positions(numberOfParticles); std::vector<Vec3> positions(numberOfParticles);
...@@ -375,7 +328,7 @@ void testOneAngle( FILE* log ) { ...@@ -375,7 +328,7 @@ void testOneAngle( FILE* log ) {
positions[3] = Vec3(1, 1, 1); positions[3] = Vec3(1, 1, 1);
context.setPositions(positions); context.setPositions(positions);
compareWithExpectedForceAndEnergy( context, *amoebaInPlaneAngleForce, TOL, "testOneInPlaneAngle", log ); compareWithExpectedForceAndEnergy(context, *amoebaInPlaneAngleForce, TOL, "testOneInPlaneAngle");
// Try changing the angle parameters and make sure it's still correct. // Try changing the angle parameters and make sure it's still correct.
...@@ -383,14 +336,14 @@ void testOneAngle( FILE* log ) { ...@@ -383,14 +336,14 @@ void testOneAngle( FILE* log ) {
bool exceptionThrown = false; bool exceptionThrown = false;
try { try {
// This should throw an exception. // This should throw an exception.
compareWithExpectedForceAndEnergy( context, *amoebaInPlaneAngleForce, TOL, "testOneInPlaneAngle", log ); compareWithExpectedForceAndEnergy(context, *amoebaInPlaneAngleForce, TOL, "testOneInPlaneAngle");
} }
catch (std::exception ex) { catch (std::exception ex) {
exceptionThrown = true; exceptionThrown = true;
} }
ASSERT(exceptionThrown); ASSERT(exceptionThrown);
amoebaInPlaneAngleForce->updateParametersInContext(context); amoebaInPlaneAngleForce->updateParametersInContext(context);
compareWithExpectedForceAndEnergy( context, *amoebaInPlaneAngleForce, TOL, "testOneInPlaneAngle", log ); compareWithExpectedForceAndEnergy(context, *amoebaInPlaneAngleForce, TOL, "testOneInPlaneAngle");
} }
int main(int argc, char* argv[]) { int main(int argc, char* argv[]) {
...@@ -399,8 +352,7 @@ int main(int argc, char* argv[]) { ...@@ -399,8 +352,7 @@ int main(int argc, char* argv[]) {
registerAmoebaCudaKernelFactories(); registerAmoebaCudaKernelFactories();
if (argc > 1) if (argc > 1)
Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1])); Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1]));
FILE* log = NULL; testOneAngle();
testOneAngle( NULL );
} catch(const std::exception& e) { } catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl; std::cout << "exception: " << e.what() << std::endl;
std::cout << "FAIL - ERROR. Test failed." << std::endl; std::cout << "FAIL - ERROR. Test failed." << std::endl;
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaMultipoleForce.cpp
View file @ f11bafd9
This source diff could not be displayed because it is too large. You can view the blob instead.
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaOutOfPlaneBendForce.cpp
View file @ f11bafd9
...@@ -63,7 +63,7 @@ const double TOL = 1e-3; ...@@ -63,7 +63,7 @@ const double TOL = 1e-3;
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
static void crossProductVector3( double* vectorX, double* vectorY, double* vectorZ ){ static void crossProductVector3(double* vectorX, double* vectorY, double* vectorZ) {
vectorZ[0] = vectorX[1]*vectorY[2] - vectorX[2]*vectorY[1]; vectorZ[0] = vectorX[1]*vectorY[2] - vectorX[2]*vectorY[1];
vectorZ[1] = vectorX[2]*vectorY[0] - vectorX[0]*vectorY[2]; vectorZ[1] = vectorX[2]*vectorY[0] - vectorX[0]*vectorY[2];
...@@ -72,13 +72,13 @@ static void crossProductVector3( double* vectorX, double* vectorY, double* vecto ...@@ -72,13 +72,13 @@ static void crossProductVector3( double* vectorX, double* vectorY, double* vecto
return; return;
} }
static double dotVector3( double* vectorX, double* vectorY ){ static double dotVector3(double* vectorX, double* vectorY) {
return vectorX[0]*vectorY[0] + vectorX[1]*vectorY[1] + vectorX[2]*vectorY[2]; return vectorX[0]*vectorY[0] + vectorX[1]*vectorY[1] + vectorX[2]*vectorY[2];
} }
static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>& positions, AmoebaOutOfPlaneBendForce& amoebaOutOfPlaneBendForce, static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>& positions, AmoebaOutOfPlaneBendForce& amoebaOutOfPlaneBendForce,
std::vector<Vec3>& forces, double* energy, FILE* log ) { std::vector<Vec3>& forces, double* energy) {
double kAngleCubic = amoebaOutOfPlaneBendForce.getAmoebaGlobalOutOfPlaneBendCubic(); double kAngleCubic = amoebaOutOfPlaneBendForce.getAmoebaGlobalOutOfPlaneBendCubic();
...@@ -88,15 +88,7 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>& ...@@ -88,15 +88,7 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>&
int particle1, particle2, particle3, particle4; int particle1, particle2, particle3, particle4;
double kAngleQuadratic; double kAngleQuadratic;
amoebaOutOfPlaneBendForce.getOutOfPlaneBendParameters(bondIndex, particle1, particle2, particle3, particle4, kAngleQuadratic ); amoebaOutOfPlaneBendForce.getOutOfPlaneBendParameters(bondIndex, particle1, particle2, particle3, particle4, kAngleQuadratic);
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaOutOfPlaneBendForce: bond %d [%d %d %d %d] k=[%10.3e %10.3e %10.3e %10.3e %10.3e]\n",
bondIndex, particle1, particle2, particle3, particle4, kAngleQuadratic, kAngleCubic, kAngleQuartic, kAnglePentic, kAngleSextic );
(void) fflush( log );
}
#endif
enum { A, B, C, D, LastAtomIndex }; enum { A, B, C, D, LastAtomIndex };
enum { AB, CB, DB, AD, CD, LastDeltaIndex }; enum { AB, CB, DB, AD, CD, LastDeltaIndex };
...@@ -107,7 +99,7 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>& ...@@ -107,7 +99,7 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>&
// and various intermediate terms // and various intermediate terms
double deltaR[LastDeltaIndex][6]; double deltaR[LastDeltaIndex][6];
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
deltaR[AB][ii] = positions[particle1][ii] - positions[particle2][ii]; deltaR[AB][ii] = positions[particle1][ii] - positions[particle2][ii];
deltaR[CB][ii] = positions[particle3][ii] - positions[particle2][ii]; deltaR[CB][ii] = positions[particle3][ii] - positions[particle2][ii];
deltaR[DB][ii] = positions[particle4][ii] - positions[particle2][ii]; deltaR[DB][ii] = positions[particle4][ii] - positions[particle2][ii];
...@@ -115,37 +107,29 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>& ...@@ -115,37 +107,29 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>&
deltaR[CD][ii] = positions[particle3][ii] - positions[particle4][ii]; deltaR[CD][ii] = positions[particle3][ii] - positions[particle4][ii];
} }
double rDB2 = dotVector3( deltaR[DB], deltaR[DB] ); double rDB2 = dotVector3(deltaR[DB], deltaR[DB]);
double rAD2 = dotVector3( deltaR[AD], deltaR[AD] ); double rAD2 = dotVector3(deltaR[AD], deltaR[AD]);
double rCD2 = dotVector3( deltaR[CD], deltaR[CD] ); double rCD2 = dotVector3(deltaR[CD], deltaR[CD]);
double tempVector[3]; double tempVector[3];
crossProductVector3( deltaR[CB], deltaR[DB], tempVector ); crossProductVector3(deltaR[CB], deltaR[DB], tempVector);
double eE = dotVector3( deltaR[AB], tempVector ); double eE = dotVector3(deltaR[AB], tempVector );
double dot = dotVector3( deltaR[AD], deltaR[CD] ); double dot = dotVector3(deltaR[AD], deltaR[CD]);
double cc = rAD2*rCD2 - dot*dot; double cc = rAD2*rCD2 - dot*dot;
if( rDB2 <= 0.0 || cc == 0.0 ){ if (rDB2 <= 0.0 || cc == 0.0) {
return; return;
} }
double bkk2 = rDB2 - eE*eE/cc; double bkk2 = rDB2 - eE*eE/cc;
double cosine = sqrt(bkk2/rDB2); double cosine = sqrt(bkk2/rDB2);
double angle; double angle;
if( cosine >= 1.0 ){ if (cosine >= 1.0) {
angle = 0.0; angle = 0.0;
} else if( cosine <= -1.0 ){ } else if (cosine <= -1.0) {
angle = PI_M; angle = PI_M;
} else { } else {
angle = RADIAN*acos(cosine); angle = RADIAN*acos(cosine);
} }
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaOutOfPlaneBendForce: bkk2=%14.7e rDB2=%14.7e cos=%14.7e dt=%14.7e]\n",
bkk2, rDB2, cosine, angle );
(void) fflush( log );
}
#endif
// chain rule // chain rule
...@@ -159,23 +143,23 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>& ...@@ -159,23 +143,23 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>&
double dEdCos; double dEdCos;
dEdCos = dEdDt/sqrt(cc*bkk2); dEdCos = dEdDt/sqrt(cc*bkk2);
if( eE > 0.0 ){ if (eE > 0.0) {
dEdCos *= -1.0; dEdCos *= -1.0;
} }
double term = eE/cc; double term = eE/cc;
double dccd[LastAtomIndex][3]; double dccd[LastAtomIndex][3];
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
dccd[A][ii] = (deltaR[AD][ii]*rCD2 - deltaR[CD][ii]*dot)*term; dccd[A][ii] = (deltaR[AD][ii]*rCD2 - deltaR[CD][ii]*dot)*term;
dccd[C][ii] = (deltaR[CD][ii]*rAD2 - deltaR[AD][ii]*dot)*term; dccd[C][ii] = (deltaR[CD][ii]*rAD2 - deltaR[AD][ii]*dot)*term;
dccd[D][ii] = -1.0*(dccd[A][ii] + dccd[C][ii]); dccd[D][ii] = -1.0*(dccd[A][ii] + dccd[C][ii]);
} }
double deed[LastAtomIndex][3]; double deed[LastAtomIndex][3];
crossProductVector3( deltaR[DB], deltaR[CB], deed[A] ); crossProductVector3(deltaR[DB], deltaR[CB], deed[A]);
crossProductVector3( deltaR[AB], deltaR[DB], deed[C] ); crossProductVector3(deltaR[AB], deltaR[DB], deed[C]);
crossProductVector3( deltaR[CB], deltaR[AB], deed[D] ); crossProductVector3(deltaR[CB], deltaR[AB], deed[D]);
term = eE/rDB2; term = eE/rDB2;
deed[D][0] += deltaR[DB][0]*term; deed[D][0] += deltaR[DB][0]*term;
...@@ -187,24 +171,24 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>& ...@@ -187,24 +171,24 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>&
// forces // forces
// calculate forces for atoms a, c, d // calculate forces for atoms a, c, d
// the force for b is then -( a+ c + d) // the force for b is then -(a+ c + d)
double subForce[LastAtomIndex][3]; double subForce[LastAtomIndex][3];
for( int jj = 0; jj < LastAtomIndex; jj++ ){ for (int jj = 0; jj < LastAtomIndex; jj++) {
// A, C, D // A, C, D
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
subForce[jj][ii] = dEdCos*( dccd[jj][ii] + deed[jj][ii] ); subForce[jj][ii] = dEdCos*(dccd[jj][ii] + deed[jj][ii]);
} }
if( jj == 0 )jj++; // skip B if (jj == 0)jj++; // skip B
// now compute B // now compute B
if( jj == 3 ){ if (jj == 3) {
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
subForce[1][ii] = -1.0*(subForce[0][ii] + subForce[2][ii] + subForce[3][ii]); subForce[1][ii] = -1.0*(subForce[0][ii] + subForce[2][ii] + subForce[3][ii]);
} }
} }
...@@ -241,72 +225,47 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>& ...@@ -241,72 +225,47 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>&
return; return;
} }
static void computeAmoebaOutOfPlaneBendForces( Context& context, AmoebaOutOfPlaneBendForce& amoebaOutOfPlaneBendForce, static void computeAmoebaOutOfPlaneBendForces(Context& context, AmoebaOutOfPlaneBendForce& amoebaOutOfPlaneBendForce,
std::vector<Vec3>& expectedForces, double* expectedEnergy, FILE* log ) { std::vector<Vec3>& expectedForces, double* expectedEnergy) {
// get positions and zero forces // get positions and zero forces
State state = context.getState(State::Positions); State state = context.getState(State::Positions);
std::vector<Vec3> positions = state.getPositions(); std::vector<Vec3> positions = state.getPositions();
expectedForces.resize( positions.size() ); expectedForces.resize(positions.size());
for( unsigned int ii = 0; ii < expectedForces.size(); ii++ ){ for (unsigned int ii = 0; ii < expectedForces.size(); ii++) {
expectedForces[ii][0] = expectedForces[ii][1] = expectedForces[ii][2] = 0.0; expectedForces[ii][0] = expectedForces[ii][1] = expectedForces[ii][2] = 0.0;
} }
// calculates forces/energy // calculates forces/energy
*expectedEnergy = 0.0; *expectedEnergy = 0.0;
for( int ii = 0; ii < amoebaOutOfPlaneBendForce.getNumOutOfPlaneBends(); ii++ ){ for (int ii = 0; ii < amoebaOutOfPlaneBendForce.getNumOutOfPlaneBends(); ii++) {
computeAmoebaOutOfPlaneBendForce(ii, positions, amoebaOutOfPlaneBendForce, expectedForces, expectedEnergy, log ); computeAmoebaOutOfPlaneBendForce(ii, positions, amoebaOutOfPlaneBendForce, expectedForces, expectedEnergy);
}
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaOutOfPlaneBendForces: expected energy=%14.7e\n", *expectedEnergy );
for( unsigned int ii = 0; ii < positions.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e]\n", ii, expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2] );
}
(void) fflush( log );
} }
#endif
return;
} }
void compareWithExpectedForceAndEnergy( Context& context, AmoebaOutOfPlaneBendForce& amoebaOutOfPlaneBendForce, void compareWithExpectedForceAndEnergy(Context& context, AmoebaOutOfPlaneBendForce& amoebaOutOfPlaneBendForce,
double tolerance, const std::string& idString, FILE* log) { double tolerance, const std::string& idString) {
std::vector<Vec3> expectedForces; std::vector<Vec3> expectedForces;
double expectedEnergy; double expectedEnergy;
computeAmoebaOutOfPlaneBendForces( context, amoebaOutOfPlaneBendForce, expectedForces, &expectedEnergy, log ); computeAmoebaOutOfPlaneBendForces(context, amoebaOutOfPlaneBendForce, expectedForces, &expectedEnergy);
State state = context.getState(State::Forces | State::Energy); State state = context.getState(State::Forces | State::Energy);
const std::vector<Vec3> forces = state.getForces(); const std::vector<Vec3> forces = state.getForces();
for (unsigned int ii = 0; ii < forces.size(); ii++) {
#ifdef AMOEBA_DEBUG ASSERT_EQUAL_VEC(expectedForces[ii], forces[ii], tolerance);
if( log ){
(void) fprintf( log, "computeAmoebaOutOfPlaneBendForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy() );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2], forces[ii][0], forces[ii][1], forces[ii][2] );
}
(void) fflush( log );
} }
#endif ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), tolerance);
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
ASSERT_EQUAL_VEC( expectedForces[ii], forces[ii], tolerance );
}
ASSERT_EQUAL_TOL( expectedEnergy, state.getPotentialEnergy(), tolerance );
} }
void testOneOutOfPlaneBend( FILE* log ) { void testOneOutOfPlaneBend() {
System system; System system;
int numberOfParticles = 4; int numberOfParticles = 4;
for( int ii = 0; ii < numberOfParticles; ii++ ){ for (int ii = 0; ii < numberOfParticles; ii++) {
system.addParticle(1.0); system.addParticle(1.0);
} }
...@@ -314,27 +273,27 @@ void testOneOutOfPlaneBend( FILE* log ) { ...@@ -314,27 +273,27 @@ void testOneOutOfPlaneBend( FILE* log ) {
AmoebaOutOfPlaneBendForce* amoebaOutOfPlaneBendForce = new AmoebaOutOfPlaneBendForce(); AmoebaOutOfPlaneBendForce* amoebaOutOfPlaneBendForce = new AmoebaOutOfPlaneBendForce();
amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendCubic( -0.1400000E-01 ); amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendCubic( -0.1400000E-01);
amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendQuartic( 0.5600000E-04 ); amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendQuartic( 0.5600000E-04);
amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendPentic( -0.7000000E-06 ); amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendPentic( -0.7000000E-06);
amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendSextic( 0.2200000E-07 ); amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendSextic( 0.2200000E-07);
double kOutOfPlaneBend = 0.328682196E-01; double kOutOfPlaneBend = 0.328682196E-01;
amoebaOutOfPlaneBendForce->addOutOfPlaneBend(0, 1, 2, 3, kOutOfPlaneBend ); amoebaOutOfPlaneBendForce->addOutOfPlaneBend(0, 1, 2, 3, kOutOfPlaneBend);
system.addForce(amoebaOutOfPlaneBendForce); system.addForce(amoebaOutOfPlaneBendForce);
Context context(system, integrator, Platform::getPlatformByName( "CUDA")); Context context(system, integrator, Platform::getPlatformByName("CUDA"));
std::vector<Vec3> positions(numberOfParticles); std::vector<Vec3> positions(numberOfParticles);
positions[0] = Vec3( 0.262660000E+02, 0.254130000E+02, 0.284200000E+01 ); positions[0] = Vec3(0.262660000E+02, 0.254130000E+02, 0.284200000E+01);
positions[1] = Vec3( 0.269130000E+02, 0.266390000E+02, 0.353100000E+01 ); positions[1] = Vec3(0.269130000E+02, 0.266390000E+02, 0.353100000E+01);
positions[2] = Vec3( 0.278860000E+02, 0.264630000E+02, 0.426300000E+01 ); positions[2] = Vec3(0.278860000E+02, 0.264630000E+02, 0.426300000E+01);
positions[3] = Vec3( 0.245568230E+02, 0.250215290E+02, 0.796852800E+01 ); positions[3] = Vec3(0.245568230E+02, 0.250215290E+02, 0.796852800E+01);
context.setPositions(positions); context.setPositions(positions);
compareWithExpectedForceAndEnergy( context, *amoebaOutOfPlaneBendForce, TOL, "testOneOutOfPlaneBend", log ); compareWithExpectedForceAndEnergy(context, *amoebaOutOfPlaneBendForce, TOL, "testOneOutOfPlaneBend");
// Try changing the bend parameters and make sure it's still correct. // Try changing the bend parameters and make sure it's still correct.
...@@ -342,21 +301,21 @@ void testOneOutOfPlaneBend( FILE* log ) { ...@@ -342,21 +301,21 @@ void testOneOutOfPlaneBend( FILE* log ) {
bool exceptionThrown = false; bool exceptionThrown = false;
try { try {
// This should throw an exception. // This should throw an exception.
compareWithExpectedForceAndEnergy( context, *amoebaOutOfPlaneBendForce, TOL, "testOneOutOfPlaneBend", log ); compareWithExpectedForceAndEnergy(context, *amoebaOutOfPlaneBendForce, TOL, "testOneOutOfPlaneBend");
} }
catch (std::exception ex) { catch (std::exception ex) {
exceptionThrown = true; exceptionThrown = true;
} }
ASSERT(exceptionThrown); ASSERT(exceptionThrown);
amoebaOutOfPlaneBendForce->updateParametersInContext(context); amoebaOutOfPlaneBendForce->updateParametersInContext(context);
compareWithExpectedForceAndEnergy( context, *amoebaOutOfPlaneBendForce, TOL, "testOneOutOfPlaneBend", log ); compareWithExpectedForceAndEnergy(context, *amoebaOutOfPlaneBendForce, TOL, "testOneOutOfPlaneBend");
} }
void testOneOutOfPlaneBend2( FILE* log, int setId ) { void testOneOutOfPlaneBend2(int setId) {
System system; System system;
int numberOfParticles = 4; int numberOfParticles = 4;
for( int ii = 0; ii < numberOfParticles; ii++ ){ for (int ii = 0; ii < numberOfParticles; ii++) {
system.addParticle(1.0); system.addParticle(1.0);
} }
...@@ -364,10 +323,10 @@ void testOneOutOfPlaneBend2( FILE* log, int setId ) { ...@@ -364,10 +323,10 @@ void testOneOutOfPlaneBend2( FILE* log, int setId ) {
AmoebaOutOfPlaneBendForce* amoebaOutOfPlaneBendForce = new AmoebaOutOfPlaneBendForce(); AmoebaOutOfPlaneBendForce* amoebaOutOfPlaneBendForce = new AmoebaOutOfPlaneBendForce();
amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendCubic( -0.1400000E-01 ); amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendCubic( -0.1400000E-01);
amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendQuartic( 0.5600000E-04 ); amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendQuartic( 0.5600000E-04);
amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendPentic( -0.7000000E-06 ); amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendPentic( -0.7000000E-06);
amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendSextic( 0.2200000E-07 ); amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendSextic( 0.2200000E-07);
/* /*
285 441 442 443 444 0.328682196E-01 285 441 442 443 444 0.328682196E-01
286 441 442 444 443 0.164493407E-01 286 441 442 444 443 0.164493407E-01
...@@ -387,139 +346,106 @@ void testOneOutOfPlaneBend2( FILE* log, int setId ) { ...@@ -387,139 +346,106 @@ void testOneOutOfPlaneBend2( FILE* log, int setId ) {
*/ */
std::map<int,Vec3> coordinates; std::map<int,Vec3> coordinates;
coordinates[440] = Vec3( 0.893800000E+01, 0.439800000E+01, 0.343100000E+01 ); coordinates[440] = Vec3( 0.893800000E+01, 0.439800000E+01, 0.343100000E+01);
coordinates[441] = Vec3( 0.779100000E+01, 0.614600000E+01, 0.390100000E+01 ); coordinates[441] = Vec3( 0.779100000E+01, 0.614600000E+01, 0.390100000E+01);
coordinates[442] = Vec3( 0.915400000E+01, 0.683900000E+01, 0.389400000E+01 ); coordinates[442] = Vec3( 0.915400000E+01, 0.683900000E+01, 0.389400000E+01);
coordinates[443] = Vec3( 0.101770000E+02, 0.619000000E+01, 0.379900000E+01 ); coordinates[443] = Vec3( 0.101770000E+02, 0.619000000E+01, 0.379900000E+01);
coordinates[444] = Vec3( 0.921000000E+01, 0.813800000E+01, 0.398600000E+01 ); coordinates[444] = Vec3( 0.921000000E+01, 0.813800000E+01, 0.398600000E+01);
coordinates[445] = Vec3( 0.708500000E+01, 0.672900000E+01, 0.332700000E+01 ); coordinates[445] = Vec3( 0.708500000E+01, 0.672900000E+01, 0.332700000E+01);
coordinates[446] = Vec3( 0.744300000E+01, 0.605200000E+01, 0.491900000E+01 ); coordinates[446] = Vec3( 0.744300000E+01, 0.605200000E+01, 0.491900000E+01);
coordinates[447] = Vec3( 0.100820000E+02, 0.859300000E+01, 0.398200000E+01 ); coordinates[447] = Vec3( 0.100820000E+02, 0.859300000E+01, 0.398200000E+01);
coordinates[448] = Vec3( 0.838000000E+01, 0.866100000E+01, 0.406000000E+01 ); coordinates[448] = Vec3( 0.838000000E+01, 0.866100000E+01, 0.406000000E+01);
double kOutOfPlaneBend = 0.328682196E-01; double kOutOfPlaneBend = 0.328682196E-01;
std::vector<int> particleIndices; std::vector<int> particleIndices;
if( setId == 1 ){ if (setId == 1) {
particleIndices.push_back( 441 ); particleIndices.push_back(441);
particleIndices.push_back( 442 ); particleIndices.push_back(442);
particleIndices.push_back( 443 ); particleIndices.push_back(443);
particleIndices.push_back( 444 ); particleIndices.push_back(444);
kOutOfPlaneBend = 0.328682196E-01; kOutOfPlaneBend = 0.328682196E-01;
} else if( setId == 2 ){ } else if (setId == 2) {
particleIndices.push_back( 441 ); particleIndices.push_back(441);
particleIndices.push_back( 442 ); particleIndices.push_back(442);
particleIndices.push_back( 444 ); particleIndices.push_back(444);
particleIndices.push_back( 443 ); particleIndices.push_back(443);
kOutOfPlaneBend = 0.164493407E-01; kOutOfPlaneBend = 0.164493407E-01;
} else if( setId == 3 ){ } else if (setId == 3) {
particleIndices.push_back( 443 ); particleIndices.push_back(443);
particleIndices.push_back( 442 ); particleIndices.push_back(442);
particleIndices.push_back( 444 ); particleIndices.push_back(444);
particleIndices.push_back( 441 ); particleIndices.push_back(441);
kOutOfPlaneBend = 0.636650407E-02; kOutOfPlaneBend = 0.636650407E-02;
} else if( setId == 4 ){ } else if (setId == 4) {
particleIndices.push_back( 442 ); particleIndices.push_back(442);
particleIndices.push_back( 444 ); particleIndices.push_back(444);
particleIndices.push_back( 447 ); particleIndices.push_back(447);
particleIndices.push_back( 448 ); particleIndices.push_back(448);
kOutOfPlaneBend = 0.392956472E-02; kOutOfPlaneBend = 0.392956472E-02;
} else if( setId == 5 ){ } else if (setId == 5) {
particleIndices.push_back( 442 ); particleIndices.push_back(442);
particleIndices.push_back( 444 ); particleIndices.push_back(444);
particleIndices.push_back( 448 ); particleIndices.push_back(448);
particleIndices.push_back( 447 ); particleIndices.push_back(447);
kOutOfPlaneBend = 0.392956472E-02; kOutOfPlaneBend = 0.392956472E-02;
} else if( setId == 6 ){ } else if (setId == 6) {
particleIndices.push_back( 447 ); particleIndices.push_back(447);
particleIndices.push_back( 444 ); particleIndices.push_back(444);
particleIndices.push_back( 448 ); particleIndices.push_back(448);
particleIndices.push_back( 442 ); particleIndices.push_back(442);
kOutOfPlaneBend = 0.214755281E-01; kOutOfPlaneBend = 0.214755281E-01;
} else { } else {
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "Set id %d not recognized.\n", setId );
}
#endif
std::stringstream buffer; std::stringstream buffer;
buffer << "Set id " << setId << " not recognized."; buffer << "Set id " << setId << " not recognized.";
throw OpenMMException( buffer.str() ); throw OpenMMException(buffer.str());
} }
#ifdef AMOEBA_DEBUG amoebaOutOfPlaneBendForce->addOutOfPlaneBend(0, 1, 2, 3, kOutOfPlaneBend);
if( log ){
(void) fprintf( log, "Set id %d.\n", setId );
}
#endif
amoebaOutOfPlaneBendForce->addOutOfPlaneBend(0, 1, 2, 3, kOutOfPlaneBend );
system.addForce(amoebaOutOfPlaneBendForce); system.addForce(amoebaOutOfPlaneBendForce);
Context context(system, integrator, Platform::getPlatformByName( "CUDA")); Context context(system, integrator, Platform::getPlatformByName("CUDA"));
std::vector<Vec3> positions(numberOfParticles); std::vector<Vec3> positions(numberOfParticles);
for( unsigned int ii = 0; ii < numberOfParticles; ii++ ){ for (unsigned int ii = 0; ii < numberOfParticles; ii++) {
if( coordinates.find( particleIndices[ii] ) == coordinates.end() ){ if (coordinates.find(particleIndices[ii]) == coordinates.end()) {
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "Coordinates %d not loaded.", particleIndices[ii] );
}
#endif
std::stringstream buffer; std::stringstream buffer;
buffer << "Coordinates " << particleIndices[ii] << " not loaded."; buffer << "Coordinates " << particleIndices[ii] << " not loaded.";
throw OpenMMException( buffer.str() ); throw OpenMMException(buffer.str());
} }
positions[ii] = coordinates[particleIndices[ii]]; positions[ii] = coordinates[particleIndices[ii]];
} }
context.setPositions(positions); context.setPositions(positions);
compareWithExpectedForceAndEnergy( context, *amoebaOutOfPlaneBendForce, TOL, "testOneOutOfPlaneBend", log ); compareWithExpectedForceAndEnergy(context, *amoebaOutOfPlaneBendForce, TOL, "testOneOutOfPlaneBend");
static int iter = 0; static int iter = 0;
static std::map<int,Vec3> totalForces; static std::map<int,Vec3> totalForces;
static double totalEnergy; static double totalEnergy;
if( iter == 0 ){ if (iter == 0) {
totalForces[441] = Vec3( 0.0, 0.0, 0.0 ); totalForces[441] = Vec3( 0.0, 0.0, 0.0);
totalForces[442] = Vec3( 0.0, 0.0, 0.0 ); totalForces[442] = Vec3( 0.0, 0.0, 0.0);
totalForces[443] = Vec3( 0.0, 0.0, 0.0 ); totalForces[443] = Vec3( 0.0, 0.0, 0.0);
totalForces[444] = Vec3( 0.0, 0.0, 0.0 ); totalForces[444] = Vec3( 0.0, 0.0, 0.0);
totalForces[445] = Vec3( 0.0, 0.0, 0.0 ); totalForces[445] = Vec3( 0.0, 0.0, 0.0);
totalForces[446] = Vec3( 0.0, 0.0, 0.0 ); totalForces[446] = Vec3( 0.0, 0.0, 0.0);
totalForces[447] = Vec3( 0.0, 0.0, 0.0 ); totalForces[447] = Vec3( 0.0, 0.0, 0.0);
totalForces[448] = Vec3( 0.0, 0.0, 0.0 ); totalForces[448] = Vec3( 0.0, 0.0, 0.0);
totalForces[449] = Vec3( 0.0, 0.0, 0.0 ); totalForces[449] = Vec3( 0.0, 0.0, 0.0);
totalEnergy = 0.0; totalEnergy = 0.0;
} }
iter++; iter++;
std::vector<Vec3> forces; std::vector<Vec3> forces;
forces.resize( numberOfParticles ); forces.resize(numberOfParticles);
double energy; double energy;
computeAmoebaOutOfPlaneBendForce( 0, positions, *amoebaOutOfPlaneBendForce, forces, &energy, log ); computeAmoebaOutOfPlaneBendForce(0, positions, *amoebaOutOfPlaneBendForce, forces, &energy);
totalEnergy += energy; totalEnergy += energy;
for( unsigned int ii = 0; ii < numberOfParticles; ii++ ){ for (unsigned int ii = 0; ii < numberOfParticles; ii++) {
for( unsigned int jj = 0; jj < 3; jj++ ){ for (unsigned int jj = 0; jj < 3; jj++) {
totalForces[particleIndices[ii]][jj] += forces[ii][jj]; totalForces[particleIndices[ii]][jj] += forces[ii][jj];
} }
} }
if( iter == 6 ){
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaOutOfPlaneBendForces: energy=%14.7e\n", totalEnergy);
for( std::map<int,Vec3>::iterator ii = totalForces.begin(); ii != totalForces.end(); ii++ ){
int particleIndex = ii->first;
Vec3 forces = ii->second;
(void) fprintf( log, "%6d [%14.7e %14.7e %14.7e] \n", particleIndex,
forces[0], forces[1], forces[2] );
}
(void) fflush( log );
}
#endif
}
} }
int main(int argc, char* argv[]) { int main(int argc, char* argv[]) {
...@@ -528,9 +454,8 @@ int main(int argc, char* argv[]) { ...@@ -528,9 +454,8 @@ int main(int argc, char* argv[]) {
registerAmoebaCudaKernelFactories(); registerAmoebaCudaKernelFactories();
if (argc > 1) if (argc > 1)
Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1])); Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1]));
FILE* log = NULL;
testOneOutOfPlaneBend( log ); testOneOutOfPlaneBend();
} catch(const std::exception& e) { } catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl; std::cout << "exception: " << e.what() << std::endl;
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaPiTorsionForce.cpp
View file @ f11bafd9
...@@ -63,7 +63,7 @@ const double TOL = 1e-5; ...@@ -63,7 +63,7 @@ const double TOL = 1e-5;
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
static void crossProductVector3( double* vectorX, double* vectorY, double* vectorZ ){ static void crossProductVector3(double* vectorX, double* vectorY, double* vectorZ) {
vectorZ[0] = vectorX[1]*vectorY[2] - vectorX[2]*vectorY[1]; vectorZ[0] = vectorX[1]*vectorY[2] - vectorX[2]*vectorY[1];
vectorZ[1] = vectorX[2]*vectorY[0] - vectorX[0]*vectorY[2]; vectorZ[1] = vectorX[2]*vectorY[0] - vectorX[0]*vectorY[2];
...@@ -72,43 +72,35 @@ static void crossProductVector3( double* vectorX, double* vectorY, double* vecto ...@@ -72,43 +72,35 @@ static void crossProductVector3( double* vectorX, double* vectorY, double* vecto
return; return;
} }
static double dotVector3( double* vectorX, double* vectorY ){ static double dotVector3(double* vectorX, double* vectorY) {
return vectorX[0]*vectorY[0] + vectorX[1]*vectorY[1] + vectorX[2]*vectorY[2]; return vectorX[0]*vectorY[0] + vectorX[1]*vectorY[1] + vectorX[2]*vectorY[2];
} }
static void computeAmoebaPiTorsionForce(int bondIndex, std::vector<Vec3>& positions, AmoebaPiTorsionForce& amoebaPiTorsionForce, static void computeAmoebaPiTorsionForce(int bondIndex, std::vector<Vec3>& positions, AmoebaPiTorsionForce& amoebaPiTorsionForce,
std::vector<Vec3>& forces, double* energy, FILE* log ) { std::vector<Vec3>& forces, double* energy) {
int particle1, particle2, particle3, particle4, particle5, particle6; int particle1, particle2, particle3, particle4, particle5, particle6;
double kTorsion; double kTorsion;
amoebaPiTorsionForce.getPiTorsionParameters(bondIndex, particle1, particle2, particle3, particle4, particle5, particle6, kTorsion); amoebaPiTorsionForce.getPiTorsionParameters(bondIndex, particle1, particle2, particle3, particle4, particle5, particle6, kTorsion);
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaPiTorsionForce: bond %d [%d %d %d %d %d %d] k=%10.3e\n",
bondIndex, particle1, particle2, particle3, particle4, particle5, particle6, kTorsion );
(void) fflush( log );
}
#endif
enum { AD, BD, EC, FC, P, Q, CP, DC, QD, T, U, TU, DP, QC, dT, dU, dP, dQ, dC1, dC2, dD1, dD2, LastDeltaIndex }; enum { AD, BD, EC, FC, P, Q, CP, DC, QD, T, U, TU, DP, QC, dT, dU, dP, dQ, dC1, dC2, dD1, dD2, LastDeltaIndex };
double deltaR[LastDeltaIndex][3]; double deltaR[LastDeltaIndex][3];
enum { A, B, C, D, E, F, LastAtomIndex }; enum { A, B, C, D, E, F, LastAtomIndex };
double d[LastAtomIndex][3]; double d[LastAtomIndex][3];
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
deltaR[AD][ii] = positions[particle1][ii] - positions[particle4][ii]; deltaR[AD][ii] = positions[particle1][ii] - positions[particle4][ii];
deltaR[BD][ii] = positions[particle2][ii] - positions[particle4][ii]; deltaR[BD][ii] = positions[particle2][ii] - positions[particle4][ii];
deltaR[EC][ii] = positions[particle5][ii] - positions[particle3][ii]; deltaR[EC][ii] = positions[particle5][ii] - positions[particle3][ii];
deltaR[FC][ii] = positions[particle6][ii] - positions[particle3][ii]; deltaR[FC][ii] = positions[particle6][ii] - positions[particle3][ii];
} }
crossProductVector3( deltaR[AD], deltaR[BD], deltaR[P] ); crossProductVector3(deltaR[AD], deltaR[BD], deltaR[P]);
crossProductVector3( deltaR[EC], deltaR[FC], deltaR[Q] ); crossProductVector3(deltaR[EC], deltaR[FC], deltaR[Q]);
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
deltaR[CP][ii] = -deltaR[P][ii]; deltaR[CP][ii] = -deltaR[P][ii];
deltaR[DC][ii] = positions[particle4][ii] - positions[particle3][ii]; deltaR[DC][ii] = positions[particle4][ii] - positions[particle3][ii];
deltaR[QD][ii] = deltaR[Q][ii]; deltaR[QD][ii] = deltaR[Q][ii];
...@@ -116,25 +108,25 @@ static void computeAmoebaPiTorsionForce(int bondIndex, std::vector<Vec3>& posit ...@@ -116,25 +108,25 @@ static void computeAmoebaPiTorsionForce(int bondIndex, std::vector<Vec3>& posit
deltaR[P][ii] += positions[particle3][ii]; deltaR[P][ii] += positions[particle3][ii];
deltaR[Q][ii] += positions[particle4][ii]; deltaR[Q][ii] += positions[particle4][ii];
} }
crossProductVector3( deltaR[CP], deltaR[DC], deltaR[T] ); crossProductVector3(deltaR[CP], deltaR[DC], deltaR[T] );
crossProductVector3( deltaR[DC], deltaR[QD], deltaR[U] ); crossProductVector3(deltaR[DC], deltaR[QD], deltaR[U] );
crossProductVector3( deltaR[T], deltaR[U], deltaR[TU] ); crossProductVector3(deltaR[T], deltaR[U], deltaR[TU]);
double rT2 = dotVector3( deltaR[T], deltaR[T] ); double rT2 = dotVector3(deltaR[T], deltaR[T]);
double rU2 = dotVector3( deltaR[U], deltaR[U] ); double rU2 = dotVector3(deltaR[U], deltaR[U]);
double rTrU = sqrt( rT2*rU2 ); double rTrU = sqrt(rT2*rU2);
if( rTrU <= 0.0 ){ if (rTrU <= 0.0) {
return; return;
} }
double rDC = dotVector3( deltaR[DC], deltaR[DC] ); double rDC = dotVector3(deltaR[DC], deltaR[DC]);
rDC = sqrt( rDC ); rDC = sqrt(rDC);
double cosine = dotVector3( deltaR[T], deltaR[U] ); double cosine = dotVector3(deltaR[T], deltaR[U]);
cosine /= rTrU; cosine /= rTrU;
double sine = dotVector3( deltaR[DC], deltaR[TU] ); double sine = dotVector3(deltaR[DC], deltaR[TU]);
sine /= ( rDC*rTrU ); sine /= (rDC*rTrU);
double cosine2 = cosine*cosine - sine*sine; double cosine2 = cosine*cosine - sine*sine;
double sine2 = 2.0*cosine*sine; double sine2 = 2.0*cosine*sine;
...@@ -144,37 +136,37 @@ static void computeAmoebaPiTorsionForce(int bondIndex, std::vector<Vec3>& posit ...@@ -144,37 +136,37 @@ static void computeAmoebaPiTorsionForce(int bondIndex, std::vector<Vec3>& posit
double dedphi = kTorsion*dphi2; double dedphi = kTorsion*dphi2;
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
deltaR[DP][ii] = positions[particle4][ii] - deltaR[P][ii]; deltaR[DP][ii] = positions[particle4][ii] - deltaR[P][ii];
deltaR[QC][ii] = deltaR[Q][ii] - positions[particle3][ii]; deltaR[QC][ii] = deltaR[Q][ii] - positions[particle3][ii];
} }
double factorT = dedphi/( rDC*rT2 ); double factorT = dedphi/(rDC*rT2);
double factorU = -dedphi/( rDC*rU2 ); double factorU = -dedphi/(rDC*rU2);
crossProductVector3( deltaR[T], deltaR[DC], deltaR[dT] ); crossProductVector3(deltaR[T], deltaR[DC], deltaR[dT] );
crossProductVector3( deltaR[U], deltaR[DC], deltaR[dU] ); crossProductVector3(deltaR[U], deltaR[DC], deltaR[dU] );
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
deltaR[dT][ii] *= factorT; deltaR[dT][ii] *= factorT;
deltaR[dU][ii] *= factorU; deltaR[dU][ii] *= factorU;
} }
crossProductVector3( deltaR[dT], deltaR[DC], deltaR[dP] ); crossProductVector3(deltaR[dT], deltaR[DC], deltaR[dP] );
crossProductVector3( deltaR[dU], deltaR[DC], deltaR[dQ] ); crossProductVector3(deltaR[dU], deltaR[DC], deltaR[dQ] );
crossProductVector3( deltaR[DP], deltaR[dT], deltaR[dC1] ); crossProductVector3(deltaR[DP], deltaR[dT], deltaR[dC1] );
crossProductVector3( deltaR[dU], deltaR[QD], deltaR[dC2] ); crossProductVector3(deltaR[dU], deltaR[QD], deltaR[dC2] );
crossProductVector3( deltaR[dT], deltaR[CP], deltaR[dD1] ); crossProductVector3(deltaR[dT], deltaR[CP], deltaR[dD1] );
crossProductVector3( deltaR[QC], deltaR[dU], deltaR[dD2] ); crossProductVector3(deltaR[QC], deltaR[dU], deltaR[dD2] );
crossProductVector3( deltaR[BD], deltaR[dP], d[A] ); crossProductVector3(deltaR[BD], deltaR[dP], d[A] );
crossProductVector3( deltaR[dP], deltaR[AD], d[B] ); crossProductVector3(deltaR[dP], deltaR[AD], d[B] );
crossProductVector3( deltaR[FC], deltaR[dQ], d[E] ); crossProductVector3(deltaR[FC], deltaR[dQ], d[E] );
crossProductVector3( deltaR[dQ], deltaR[EC], d[F] ); crossProductVector3(deltaR[dQ], deltaR[EC], d[F] );
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
d[C][ii] = deltaR[dC1][ii] + deltaR[dC2][ii] + deltaR[dP][ii] - d[E][ii] - d[F][ii]; d[C][ii] = deltaR[dC1][ii] + deltaR[dC2][ii] + deltaR[dP][ii] - d[E][ii] - d[F][ii];
d[D][ii] = deltaR[dD1][ii] + deltaR[dD2][ii] + deltaR[dQ][ii] - d[A][ii] - d[B][ii]; d[D][ii] = deltaR[dD1][ii] + deltaR[dD2][ii] + deltaR[dQ][ii] - d[A][ii] - d[B][ii];
} }
...@@ -212,71 +204,48 @@ static void computeAmoebaPiTorsionForce(int bondIndex, std::vector<Vec3>& posit ...@@ -212,71 +204,48 @@ static void computeAmoebaPiTorsionForce(int bondIndex, std::vector<Vec3>& posit
return; return;
} }
static void computeAmoebaPiTorsionForces( Context& context, AmoebaPiTorsionForce& amoebaPiTorsionForce, static void computeAmoebaPiTorsionForces(Context& context, AmoebaPiTorsionForce& amoebaPiTorsionForce,
std::vector<Vec3>& expectedForces, double* expectedEnergy, FILE* log ) { std::vector<Vec3>& expectedForces, double* expectedEnergy) {
// get positions and zero forces // get positions and zero forces
State state = context.getState(State::Positions); State state = context.getState(State::Positions);
std::vector<Vec3> positions = state.getPositions(); std::vector<Vec3> positions = state.getPositions();
expectedForces.resize( positions.size() ); expectedForces.resize(positions.size());
for( unsigned int ii = 0; ii < expectedForces.size(); ii++ ){ for (unsigned int ii = 0; ii < expectedForces.size(); ii++) {
expectedForces[ii][0] = expectedForces[ii][1] = expectedForces[ii][2] = 0.0; expectedForces[ii][0] = expectedForces[ii][1] = expectedForces[ii][2] = 0.0;
} }
// calculates forces/energy // calculates forces/energy
*expectedEnergy = 0.0; *expectedEnergy = 0.0;
for( int ii = 0; ii < amoebaPiTorsionForce.getNumPiTorsions(); ii++ ){ for (int ii = 0; ii < amoebaPiTorsionForce.getNumPiTorsions(); ii++) {
computeAmoebaPiTorsionForce(ii, positions, amoebaPiTorsionForce, expectedForces, expectedEnergy, log ); computeAmoebaPiTorsionForce(ii, positions, amoebaPiTorsionForce, expectedForces, expectedEnergy);
}
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaPiTorsionForces: expected energy=%14.7e\n", *expectedEnergy );
for( unsigned int ii = 0; ii < positions.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e]\n", ii, expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2] );
}
(void) fflush( log );
} }
#endif
return;
} }
void compareWithExpectedForceAndEnergy( Context& context, AmoebaPiTorsionForce& amoebaPiTorsionForce, void compareWithExpectedForceAndEnergy(Context& context, AmoebaPiTorsionForce& amoebaPiTorsionForce,
double tolerance, const std::string& idString, FILE* log) { double tolerance, const std::string& idString) {
std::vector<Vec3> expectedForces; std::vector<Vec3> expectedForces;
double expectedEnergy; double expectedEnergy;
computeAmoebaPiTorsionForces( context, amoebaPiTorsionForce, expectedForces, &expectedEnergy, log ); computeAmoebaPiTorsionForces(context, amoebaPiTorsionForce, expectedForces, &expectedEnergy);
State state = context.getState(State::Forces | State::Energy); State state = context.getState(State::Forces | State::Energy);
const std::vector<Vec3> forces = state.getForces(); const std::vector<Vec3> forces = state.getForces();
#ifdef AMOEBA_DEBUG for (unsigned int ii = 0; ii < forces.size(); ii++) {
if( log ){ ASSERT_EQUAL_VEC(expectedForces[ii], forces[ii], tolerance);
(void) fprintf( log, "computeAmoebaPiTorsionForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy() );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2], forces[ii][0], forces[ii][1], forces[ii][2] );
}
(void) fflush( log );
}
#endif
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
ASSERT_EQUAL_VEC( expectedForces[ii], forces[ii], tolerance );
} }
ASSERT_EQUAL_TOL( expectedEnergy, state.getPotentialEnergy(), tolerance ); ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), tolerance);
} }
void testOnePiTorsion( FILE* log ) { void testOnePiTorsion() {
System system; System system;
int numberOfParticles = 6; int numberOfParticles = 6;
for( int ii = 0; ii < numberOfParticles; ii++ ){ for (int ii = 0; ii < numberOfParticles; ii++) {
system.addParticle(1.0); system.addParticle(1.0);
} }
...@@ -285,23 +254,23 @@ void testOnePiTorsion( FILE* log ) { ...@@ -285,23 +254,23 @@ void testOnePiTorsion( FILE* log ) {
AmoebaPiTorsionForce* amoebaPiTorsionForce = new AmoebaPiTorsionForce(); AmoebaPiTorsionForce* amoebaPiTorsionForce = new AmoebaPiTorsionForce();
double kTorsion = 6.85; double kTorsion = 6.85;
amoebaPiTorsionForce->addPiTorsion(0, 1, 2, 3, 4, 5, kTorsion ); amoebaPiTorsionForce->addPiTorsion(0, 1, 2, 3, 4, 5, kTorsion);
system.addForce(amoebaPiTorsionForce); system.addForce(amoebaPiTorsionForce);
Context context(system, integrator, Platform::getPlatformByName( "CUDA")); Context context(system, integrator, Platform::getPlatformByName("CUDA"));
std::vector<Vec3> positions(numberOfParticles); std::vector<Vec3> positions(numberOfParticles);
positions[0] = Vec3( 0.262660000E+02, 0.254130000E+02, 0.284200000E+01 ); positions[0] = Vec3(0.262660000E+02, 0.254130000E+02, 0.284200000E+01);
positions[1] = Vec3( 0.278860000E+02, 0.264630000E+02, 0.426300000E+01 ); positions[1] = Vec3(0.278860000E+02, 0.264630000E+02, 0.426300000E+01);
positions[2] = Vec3( 0.269130000E+02, 0.266390000E+02, 0.353100000E+01 ); positions[2] = Vec3(0.269130000E+02, 0.266390000E+02, 0.353100000E+01);
positions[3] = Vec3( 0.245568230E+02, 0.250215290E+02, 0.796852800E+01 ); positions[3] = Vec3(0.245568230E+02, 0.250215290E+02, 0.796852800E+01);
positions[4] = Vec3( 0.261000000E+02, 0.292530000E+02, 0.520200000E+01 ); positions[4] = Vec3(0.261000000E+02, 0.292530000E+02, 0.520200000E+01);
positions[5] = Vec3( 0.254124630E+02, 0.234691880E+02, 0.773335400E+01 ); positions[5] = Vec3(0.254124630E+02, 0.234691880E+02, 0.773335400E+01);
context.setPositions(positions); context.setPositions(positions);
compareWithExpectedForceAndEnergy( context, *amoebaPiTorsionForce, TOL, "testOnePiTorsion", log ); compareWithExpectedForceAndEnergy(context, *amoebaPiTorsionForce, TOL, "testOnePiTorsion");
// Try changing the torsion parameters and make sure it's still correct. // Try changing the torsion parameters and make sure it's still correct.
...@@ -309,14 +278,14 @@ void testOnePiTorsion( FILE* log ) { ...@@ -309,14 +278,14 @@ void testOnePiTorsion( FILE* log ) {
bool exceptionThrown = false; bool exceptionThrown = false;
try { try {
// This should throw an exception. // This should throw an exception.
compareWithExpectedForceAndEnergy( context, *amoebaPiTorsionForce, TOL, "testOnePiTorsion", log ); compareWithExpectedForceAndEnergy(context, *amoebaPiTorsionForce, TOL, "testOnePiTorsion");
} }
catch (std::exception ex) { catch (std::exception ex) {
exceptionThrown = true; exceptionThrown = true;
} }
ASSERT(exceptionThrown); ASSERT(exceptionThrown);
amoebaPiTorsionForce->updateParametersInContext(context); amoebaPiTorsionForce->updateParametersInContext(context);
compareWithExpectedForceAndEnergy( context, *amoebaPiTorsionForce, TOL, "testOnePiTorsion", log ); compareWithExpectedForceAndEnergy(context, *amoebaPiTorsionForce, TOL, "testOnePiTorsion");
} }
int main(int argc, char* argv[]) { int main(int argc, char* argv[]) {
...@@ -325,8 +294,7 @@ int main(int argc, char* argv[]) { ...@@ -325,8 +294,7 @@ int main(int argc, char* argv[]) {
registerAmoebaCudaKernelFactories(); registerAmoebaCudaKernelFactories();
if (argc > 1) if (argc > 1)
Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1])); Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1]));
FILE* log = NULL; testOnePiTorsion();
testOnePiTorsion( log );
} catch(const std::exception& e) { } catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl; std::cout << "exception: " << e.what() << std::endl;
std::cout << "FAIL - ERROR. Test failed." << std::endl; std::cout << "FAIL - ERROR. Test failed." << std::endl;
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaStretchBendForce.cpp
View file @ f11bafd9
...@@ -64,7 +64,7 @@ const double DegreesToRadians = PI_M/180.0; ...@@ -64,7 +64,7 @@ const double DegreesToRadians = PI_M/180.0;
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
static void crossProductVector3( double* vectorX, double* vectorY, double* vectorZ ){ static void crossProductVector3(double* vectorX, double* vectorY, double* vectorZ) {
vectorZ[0] = vectorX[1]*vectorY[2] - vectorX[2]*vectorY[1]; vectorZ[0] = vectorX[1]*vectorY[2] - vectorX[2]*vectorY[1];
vectorZ[1] = vectorX[2]*vectorY[0] - vectorX[0]*vectorY[2]; vectorZ[1] = vectorX[2]*vectorY[0] - vectorX[0]*vectorY[2];
...@@ -73,27 +73,19 @@ static void crossProductVector3( double* vectorX, double* vectorY, double* vecto ...@@ -73,27 +73,19 @@ static void crossProductVector3( double* vectorX, double* vectorY, double* vecto
return; return;
} }
static double dotVector3( double* vectorX, double* vectorY ){ static double dotVector3(double* vectorX, double* vectorY) {
return vectorX[0]*vectorY[0] + vectorX[1]*vectorY[1] + vectorX[2]*vectorY[2]; return vectorX[0]*vectorY[0] + vectorX[1]*vectorY[1] + vectorX[2]*vectorY[2];
} }
static void computeAmoebaStretchBendForce(int bondIndex, std::vector<Vec3>& positions, AmoebaStretchBendForce& amoebaStretchBendForce, static void computeAmoebaStretchBendForce(int bondIndex, std::vector<Vec3>& positions, AmoebaStretchBendForce& amoebaStretchBendForce,
std::vector<Vec3>& forces, double* energy, FILE* log ) { std::vector<Vec3>& forces, double* energy) {
int particle1, particle2, particle3; int particle1, particle2, particle3;
double abBondLength, cbBondLength, angleStretchBend, kStretchBend, k2StretchBend; double abBondLength, cbBondLength, angleStretchBend, kStretchBend, k2StretchBend;
amoebaStretchBendForce.getStretchBendParameters(bondIndex, particle1, particle2, particle3, abBondLength, cbBondLength, angleStretchBend, kStretchBend, k2StretchBend); amoebaStretchBendForce.getStretchBendParameters(bondIndex, particle1, particle2, particle3, abBondLength, cbBondLength, angleStretchBend, kStretchBend, k2StretchBend);
angleStretchBend *= RADIAN; angleStretchBend *= RADIAN;
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaStretchBendForce: bond %d [%d %d %d] ab=%10.3e cb=%10.3e angle=%10.3e k1=%10.3e k2=%10.3e\n",
bondIndex, particle1, particle2, particle3, abBondLength, cbBondLength, angleStretchBend, kStretchBend, k2StretchBend );
(void) fflush( log );
}
#endif
enum { A, B, C, LastAtomIndex }; enum { A, B, C, LastAtomIndex };
enum { AB, CB, CBxAB, ABxP, CBxP, LastDeltaIndex }; enum { AB, CB, CBxAB, ABxP, CBxP, LastDeltaIndex };
...@@ -105,31 +97,31 @@ static void computeAmoebaStretchBendForce(int bondIndex, std::vector<Vec3>& pos ...@@ -105,31 +97,31 @@ static void computeAmoebaStretchBendForce(int bondIndex, std::vector<Vec3>& pos
double deltaR[LastDeltaIndex][3]; double deltaR[LastDeltaIndex][3];
double rAB2 = 0.0; double rAB2 = 0.0;
double rCB2 = 0.0; double rCB2 = 0.0;
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
deltaR[AB][ii] = positions[particle1][ii] - positions[particle2][ii]; deltaR[AB][ii] = positions[particle1][ii] - positions[particle2][ii];
rAB2 += deltaR[AB][ii]*deltaR[AB][ii]; rAB2 += deltaR[AB][ii]*deltaR[AB][ii];
deltaR[CB][ii] = positions[particle3][ii] - positions[particle2][ii]; deltaR[CB][ii] = positions[particle3][ii] - positions[particle2][ii];
rCB2 += deltaR[CB][ii]*deltaR[CB][ii]; rCB2 += deltaR[CB][ii]*deltaR[CB][ii];
} }
double rAB = sqrt( rAB2 ); double rAB = sqrt(rAB2);
double rCB = sqrt( rCB2 ); double rCB = sqrt(rCB2);
crossProductVector3( deltaR[CB], deltaR[AB], deltaR[CBxAB] ); crossProductVector3(deltaR[CB], deltaR[AB], deltaR[CBxAB]);
double rP = dotVector3( deltaR[CBxAB], deltaR[CBxAB] ); double rP = dotVector3(deltaR[CBxAB], deltaR[CBxAB]);
rP = sqrt( rP ); rP = sqrt(rP);
if( rP <= 0.0 ){ if (rP <= 0.0) {
return; return;
} }
double dot = dotVector3( deltaR[CB], deltaR[AB] ); double dot = dotVector3(deltaR[CB], deltaR[AB]);
double cosine = dot/(rAB*rCB); double cosine = dot/(rAB*rCB);
double angle; double angle;
if( cosine >= 1.0 ){ if (cosine >= 1.0) {
angle = 0.0; angle = 0.0;
} }
else if( cosine <= -1.0 ){ else if (cosine <= -1.0) {
angle = PI_M; angle = PI_M;
} }
else { else {
...@@ -141,9 +133,9 @@ static void computeAmoebaStretchBendForce(int bondIndex, std::vector<Vec3>& pos ...@@ -141,9 +133,9 @@ static void computeAmoebaStretchBendForce(int bondIndex, std::vector<Vec3>& pos
// P = CBxAB // P = CBxAB
crossProductVector3( deltaR[AB], deltaR[CBxAB], deltaR[ABxP] ); crossProductVector3(deltaR[AB], deltaR[CBxAB], deltaR[ABxP]);
crossProductVector3( deltaR[CB], deltaR[CBxAB], deltaR[CBxP] ); crossProductVector3(deltaR[CB], deltaR[CBxAB], deltaR[CBxP]);
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
deltaR[ABxP][ii] *= termA; deltaR[ABxP][ii] *= termA;
deltaR[CBxP][ii] *= termC; deltaR[CBxP][ii] *= termC;
} }
...@@ -161,14 +153,14 @@ static void computeAmoebaStretchBendForce(int bondIndex, std::vector<Vec3>& pos ...@@ -161,14 +153,14 @@ static void computeAmoebaStretchBendForce(int bondIndex, std::vector<Vec3>& pos
// forces // forces
// calculate forces for atoms a, b, c // calculate forces for atoms a, b, c
// the force for b is then -( a + c) // the force for b is then -(a + c)
double subForce[LastAtomIndex][3]; double subForce[LastAtomIndex][3];
double dt = angle - angleStretchBend; double dt = angle - angleStretchBend;
for( int jj = 0; jj < 3; jj++ ){ for (int jj = 0; jj < 3; jj++) {
subForce[A][jj] = kStretchBend*dt*termA*deltaR[AB][jj] + drkk*deltaR[ABxP][jj]; subForce[A][jj] = kStretchBend*dt*termA*deltaR[AB][jj] + drkk*deltaR[ABxP][jj];
subForce[C][jj] = k2StretchBend*dt*termC*deltaR[CB][jj] + drkk*deltaR[CBxP][jj]; subForce[C][jj] = k2StretchBend*dt*termC*deltaR[CB][jj] + drkk*deltaR[CBxP][jj];
subForce[B][jj] = -( subForce[A][jj] + subForce[C][jj] ); subForce[B][jj] = -(subForce[A][jj] + subForce[C][jj]);
} }
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -188,81 +180,49 @@ static void computeAmoebaStretchBendForce(int bondIndex, std::vector<Vec3>& pos ...@@ -188,81 +180,49 @@ static void computeAmoebaStretchBendForce(int bondIndex, std::vector<Vec3>& pos
forces[particle3][2] -= subForce[2][2]; forces[particle3][2] -= subForce[2][2];
*energy += dt*drkk; *energy += dt*drkk;
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaStretchBendForce: angle=%10.3e dt=%10.3e dr=%10.3e\n", angle, dt, dr );
(void) fflush( log );
}
#endif
return;
} }
static void computeAmoebaStretchBendForces( Context& context, AmoebaStretchBendForce& amoebaStretchBendForce, static void computeAmoebaStretchBendForces(Context& context, AmoebaStretchBendForce& amoebaStretchBendForce,
std::vector<Vec3>& expectedForces, double* expectedEnergy, FILE* log ) { std::vector<Vec3>& expectedForces, double* expectedEnergy) {
// get positions and zero forces // get positions and zero forces
State state = context.getState(State::Positions); State state = context.getState(State::Positions);
std::vector<Vec3> positions = state.getPositions(); std::vector<Vec3> positions = state.getPositions();
expectedForces.resize( positions.size() ); expectedForces.resize(positions.size());
for( unsigned int ii = 0; ii < expectedForces.size(); ii++ ){ for (unsigned int ii = 0; ii < expectedForces.size(); ii++) {
expectedForces[ii][0] = expectedForces[ii][1] = expectedForces[ii][2] = 0.0; expectedForces[ii][0] = expectedForces[ii][1] = expectedForces[ii][2] = 0.0;
} }
// calculates forces/energy // calculates forces/energy
*expectedEnergy = 0.0; *expectedEnergy = 0.0;
for( int ii = 0; ii < amoebaStretchBendForce.getNumStretchBends(); ii++ ){ for (int ii = 0; ii < amoebaStretchBendForce.getNumStretchBends(); ii++) {
computeAmoebaStretchBendForce(ii, positions, amoebaStretchBendForce, expectedForces, expectedEnergy, log ); computeAmoebaStretchBendForce(ii, positions, amoebaStretchBendForce, expectedForces, expectedEnergy);
}
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaStretchBendForces: expected energy=%14.7e\n", *expectedEnergy );
for( unsigned int ii = 0; ii < positions.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e]\n", ii, expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2] );
}
(void) fflush( log );
} }
#endif
return;
} }
void compareWithExpectedForceAndEnergy( Context& context, AmoebaStretchBendForce& amoebaStretchBendForce, void compareWithExpectedForceAndEnergy(Context& context, AmoebaStretchBendForce& amoebaStretchBendForce,
double tolerance, const std::string& idString, FILE* log) { double tolerance, const std::string& idString) {
std::vector<Vec3> expectedForces; std::vector<Vec3> expectedForces;
double expectedEnergy; double expectedEnergy;
computeAmoebaStretchBendForces( context, amoebaStretchBendForce, expectedForces, &expectedEnergy, log ); computeAmoebaStretchBendForces(context, amoebaStretchBendForce, expectedForces, &expectedEnergy);
State state = context.getState(State::Forces | State::Energy); State state = context.getState(State::Forces | State::Energy);
const std::vector<Vec3> forces = state.getForces(); const std::vector<Vec3> forces = state.getForces();
for (unsigned int ii = 0; ii < forces.size(); ii++) {
#ifdef AMOEBA_DEBUG ASSERT_EQUAL_VEC(expectedForces[ii], forces[ii], tolerance);
if( log ){
(void) fprintf( log, "computeAmoebaStretchBendForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy() );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2], forces[ii][0], forces[ii][1], forces[ii][2] );
}
(void) fflush( log );
}
#endif
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
ASSERT_EQUAL_VEC( expectedForces[ii], forces[ii], tolerance );
} }
ASSERT_EQUAL_TOL( expectedEnergy, state.getPotentialEnergy(), tolerance ); ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), tolerance);
} }
void testOneStretchBend( FILE* log ) { void testOneStretchBend() {
System system; System system;
int numberOfParticles = 3; int numberOfParticles = 3;
for( int ii = 0; ii < numberOfParticles; ii++ ){ for (int ii = 0; ii < numberOfParticles; ii++) {
system.addParticle(1.0); system.addParticle(1.0);
} }
...@@ -276,19 +236,19 @@ void testOneStretchBend( FILE* log ) { ...@@ -276,19 +236,19 @@ void testOneStretchBend( FILE* log ) {
//double kStretchBend = 0.750491578E-01; //double kStretchBend = 0.750491578E-01;
double kStretchBend = 1.0; double kStretchBend = 1.0;
amoebaStretchBendForce->addStretchBend(0, 1, 2, abLength, cbLength, angleStretchBend, kStretchBend, kStretchBend ); amoebaStretchBendForce->addStretchBend(0, 1, 2, abLength, cbLength, angleStretchBend, kStretchBend, kStretchBend);
system.addForce(amoebaStretchBendForce); system.addForce(amoebaStretchBendForce);
Context context(system, integrator, Platform::getPlatformByName( "CUDA")); Context context(system, integrator, Platform::getPlatformByName("CUDA"));
std::vector<Vec3> positions(numberOfParticles); std::vector<Vec3> positions(numberOfParticles);
positions[0] = Vec3( 0.262660000E+02, 0.254130000E+02, 0.284200000E+01 ); positions[0] = Vec3(0.262660000E+02, 0.254130000E+02, 0.284200000E+01);
positions[1] = Vec3( 0.273400000E+02, 0.244300000E+02, 0.261400000E+01 ); positions[1] = Vec3(0.273400000E+02, 0.244300000E+02, 0.261400000E+01);
positions[2] = Vec3( 0.269573220E+02, 0.236108860E+02, 0.216376800E+01 ); positions[2] = Vec3(0.269573220E+02, 0.236108860E+02, 0.216376800E+01);
context.setPositions(positions); context.setPositions(positions);
compareWithExpectedForceAndEnergy( context, *amoebaStretchBendForce, TOL, "testOneStretchBend", log ); compareWithExpectedForceAndEnergy(context, *amoebaStretchBendForce, TOL, "testOneStretchBend");
// Try changing the stretch-bend parameters and make sure it's still correct. // Try changing the stretch-bend parameters and make sure it's still correct.
...@@ -296,14 +256,14 @@ void testOneStretchBend( FILE* log ) { ...@@ -296,14 +256,14 @@ void testOneStretchBend( FILE* log ) {
bool exceptionThrown = false; bool exceptionThrown = false;
try { try {
// This should throw an exception. // This should throw an exception.
compareWithExpectedForceAndEnergy( context, *amoebaStretchBendForce, TOL, "testOneStretchBend", log ); compareWithExpectedForceAndEnergy(context, *amoebaStretchBendForce, TOL, "testOneStretchBend");
} }
catch (std::exception ex) { catch (std::exception ex) {
exceptionThrown = true; exceptionThrown = true;
} }
ASSERT(exceptionThrown); ASSERT(exceptionThrown);
amoebaStretchBendForce->updateParametersInContext(context); amoebaStretchBendForce->updateParametersInContext(context);
compareWithExpectedForceAndEnergy( context, *amoebaStretchBendForce, TOL, "testOneStretchBend", log ); compareWithExpectedForceAndEnergy(context, *amoebaStretchBendForce, TOL, "testOneStretchBend");
} }
int main(int argc, char* argv[]) { int main(int argc, char* argv[]) {
...@@ -312,8 +272,7 @@ int main(int argc, char* argv[]) { ...@@ -312,8 +272,7 @@ int main(int argc, char* argv[]) {
registerAmoebaCudaKernelFactories(); registerAmoebaCudaKernelFactories();
if (argc > 1) if (argc > 1)
Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1])); Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1]));
FILE* log = NULL; testOneStretchBend();
testOneStretchBend( log );
} catch(const std::exception& e) { } catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl; std::cout << "exception: " << e.what() << std::endl;
std::cout << "FAIL - ERROR. Test failed." << std::endl; std::cout << "FAIL - ERROR. Test failed." << std::endl;
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaTorsionTorsionForce.cpp
View file @ f11bafd9
...@@ -48,7 +48,7 @@ extern "C" void registerAmoebaCudaKernelFactories(); ...@@ -48,7 +48,7 @@ extern "C" void registerAmoebaCudaKernelFactories();
const double TOL = 1e-4; const double TOL = 1e-4;
TorsionTorsionGrid& getTorsionGrid( int gridIndex ) { TorsionTorsionGrid& getTorsionGrid(int gridIndex) {
static double grid[4][625][6] = { static double grid[4][625][6] = {
{ {
...@@ -2558,22 +2558,22 @@ TorsionTorsionGrid& getTorsionGrid( int gridIndex ) { ...@@ -2558,22 +2558,22 @@ TorsionTorsionGrid& getTorsionGrid( int gridIndex ) {
// static std::vector< std::vector< std::vector<double> > > TorsionTorsionGrid // static std::vector< std::vector< std::vector<double> > > TorsionTorsionGrid
static std::vector<TorsionTorsionGrid> grids; static std::vector<TorsionTorsionGrid> grids;
if( grids.size() == 0 ){ if (grids.size() == 0) {
grids.resize(4); grids.resize(4);
for( int ii = 0; ii < 4; ii++ ){ for (int ii = 0; ii < 4; ii++) {
grids[ii].resize( 25 ); grids[ii].resize(25);
for( int jj = 0; jj < 25; jj++ ){ for (int jj = 0; jj < 25; jj++) {
grids[ii][jj].resize(25); grids[ii][jj].resize(25);
for( int kk = 0; kk < 25; kk++ ){ for (int kk = 0; kk < 25; kk++) {
grids[ii][jj][kk].resize(6); grids[ii][jj][kk].resize(6);
} }
} }
int index = 0; int index = 0;
for( int jj = 0; jj < 25; jj++ ){ for (int jj = 0; jj < 25; jj++) {
for( int kk = 0; kk < 25; kk++ ){ for (int kk = 0; kk < 25; kk++) {
int jjIndex = static_cast<int>(((grid[ii][index][0] + 180.0)/15.0)+1.0e-05); int jjIndex = static_cast<int>(((grid[ii][index][0] + 180.0)/15.0)+1.0e-05);
int kkIndex = static_cast<int>(((grid[ii][index][1] + 180.0)/15.0)+1.0e-05); int kkIndex = static_cast<int>(((grid[ii][index][1] + 180.0)/15.0)+1.0e-05);
for( int ll = 0; ll < 6; ll++ ){ for (int ll = 0; ll < 6; ll++) {
grids[ii][kk][jj][ll] = grid[ii][index][ll]; grids[ii][kk][jj][ll] = grid[ii][index][ll];
} }
index++; index++;
...@@ -2584,11 +2584,11 @@ TorsionTorsionGrid& getTorsionGrid( int gridIndex ) { ...@@ -2584,11 +2584,11 @@ TorsionTorsionGrid& getTorsionGrid( int gridIndex ) {
return grids[gridIndex]; return grids[gridIndex];
} }
void testTorsionTorsion( FILE* log, int systemId ) { void testTorsionTorsion(int systemId) {
System system; System system;
int numberOfParticles = 6; int numberOfParticles = 6;
for( int ii = 0; ii < numberOfParticles; ii++ ){ for (int ii = 0; ii < numberOfParticles; ii++) {
system.addParticle(1.0); system.addParticle(1.0);
} }
LangevinIntegrator integrator(0.0, 0.1, 0.01); LangevinIntegrator integrator(0.0, 0.1, 0.01);
...@@ -2600,84 +2600,71 @@ void testTorsionTorsion( FILE* log, int systemId ) { ...@@ -2600,84 +2600,71 @@ void testTorsionTorsion( FILE* log, int systemId ) {
std::vector<Vec3> positions(numberOfParticles); std::vector<Vec3> positions(numberOfParticles);
std::vector<Vec3> expectedForces(numberOfParticles); std::vector<Vec3> expectedForces(numberOfParticles);
double expectedEnergy; double expectedEnergy;
if( systemId == 0 ){ if (systemId == 0) {
// villin: 2 19 20 21 38 25 grid=2 // villin: 2 19 20 21 38 25 grid=2
chiralCheckAtomIndex = 5; chiralCheckAtomIndex = 5;
gridIndex = 2; gridIndex = 2;
positions[0] = Vec3( -0.422792800E+01, -0.110605910E+02, -0.508156700E+01 ); positions[0] = Vec3(-0.422792800E+01, -0.110605910E+02, -0.508156700E+01);
positions[1] = Vec3( -0.447153100E+01, -0.978627900E+01, -0.466405800E+01 ); positions[1] = Vec3(-0.447153100E+01, -0.978627900E+01, -0.466405800E+01);
positions[2] = Vec3( -0.531878400E+01, -0.940508600E+01, -0.352283100E+01 ); positions[2] = Vec3(-0.531878400E+01, -0.940508600E+01, -0.352283100E+01);
positions[3] = Vec3( -0.679606000E+01, -0.974353100E+01, -0.382975700E+01 ); positions[3] = Vec3(-0.679606000E+01, -0.974353100E+01, -0.382975700E+01);
positions[4] = Vec3( -0.760612300E+01, -0.992590200E+01, -0.275088400E+01 ); positions[4] = Vec3(-0.760612300E+01, -0.992590200E+01, -0.275088400E+01);
positions[5] = Vec3( -0.516893900E+01, -0.788347000E+01, -0.316943000E+01 ); positions[5] = Vec3(-0.516893900E+01, -0.788347000E+01, -0.316943000E+01);
expectedForces[0] = Vec3( 0.906091624E+00, -0.529814945E-01, 0.690384140E+00 ); expectedForces[0] = Vec3( 0.906091624E+00, -0.529814945E-01, 0.690384140E+00);
expectedForces[1] = Vec3( -0.124550232E+01, -0.999341692E+00, -0.590867130E+00 ); expectedForces[1] = Vec3(-0.124550232E+01, -0.999341692E+00, -0.590867130E+00);
expectedForces[2] = Vec3( 0.534419689E+00, 0.612404926E-01, 0.547380310E-01 ); expectedForces[2] = Vec3( 0.534419689E+00, 0.612404926E-01, 0.547380310E-01);
expectedForces[3] = Vec3( -5.732010432E-01, 2.645718463E+00, -1.585204274E-01 ); expectedForces[3] = Vec3(-5.732010432E-01, 2.645718463E+00, -1.585204274E-01);
expectedForces[4] = Vec3( 3.781920539E-01, -1.654635768E+00, 4.265386268E-03 ); expectedForces[4] = Vec3( 3.781920539E-01, -1.654635768E+00, 4.265386268E-03);
expectedForces[5] = Vec3( 0.0, 0.0, 0.0 ); expectedForces[5] = Vec3( 0.0, 0.0, 0.0);
expectedEnergy = -2.699654759E+00; expectedEnergy = -2.699654759E+00;
} else if( systemId == 1 ){ } else if (systemId == 1) {
// villin: 158 176 177 178 183 -1 0 // villin: 158 176 177 178 183 -1 0
chiralCheckAtomIndex = -1; chiralCheckAtomIndex = -1;
gridIndex = 0; gridIndex = 0;
positions[0] = Vec3( -0.105946640E+02, -0.917797000E+00, 0.105486310E+02 ); positions[0] = Vec3(-0.105946640E+02, -0.917797000E+00, 0.105486310E+02);
positions[1] = Vec3( -0.115059090E+02, -0.141876700E+01, 0.966933200E+01 ); positions[1] = Vec3(-0.115059090E+02, -0.141876700E+01, 0.966933200E+01);
positions[2] = Vec3( -0.128314660E+02, -0.876338000E+00, 0.942959800E+01 ); positions[2] = Vec3(-0.128314660E+02, -0.876338000E+00, 0.942959800E+01);
positions[3] = Vec3( -0.130879850E+02, -0.760280000E-01, 0.814732200E+01 ); positions[3] = Vec3(-0.130879850E+02, -0.760280000E-01, 0.814732200E+01);
positions[4] = Vec3( -0.120888080E+02, 0.112050000E-01, 0.722704500E+01 ); positions[4] = Vec3(-0.120888080E+02, 0.112050000E-01, 0.722704500E+01);
positions[5] = Vec3( 0.0, 0.0, 0.0 ); positions[5] = Vec3( 0.0, 0.0, 0.0);
expectedForces[0] = Vec3( 4.165851130E-01, 6.608242922E-01, -8.082168261E-01 ); expectedForces[0] = Vec3( 4.165851130E-01, 6.608242922E-01, -8.082168261E-01);
expectedForces[1] = Vec3( -6.024659721E-01, -8.878744406E-01, 1.322274444E+00 ); expectedForces[1] = Vec3(-6.024659721E-01, -8.878744406E-01, 1.322274444E+00);
expectedForces[2] = Vec3( 3.196925118E-02, -3.137497848E-01, -8.207984001E-01 ); expectedForces[2] = Vec3( 3.196925118E-02, -3.137497848E-01, -8.207984001E-01);
expectedForces[3] = Vec3( 3.842205941E-02, 2.602732089E-01, 1.547586195E-01 ); expectedForces[3] = Vec3( 3.842205941E-02, 2.602732089E-01, 1.547586195E-01);
expectedForces[4] = Vec3( 1.154895485E-01, 2.805267242E-01, 1.519821623E-01 ); expectedForces[4] = Vec3( 1.154895485E-01, 2.805267242E-01, 1.519821623E-01);
expectedForces[5] = Vec3( 0.0, 0.0, 0.0 ); expectedForces[5] = Vec3( 0.0, 0.0, 0.0);
expectedEnergy = -3.372536909E+00; expectedEnergy = -3.372536909E+00;
} }
amoebaTorsionTorsionForce->addTorsionTorsion( 0, 1, 2, 3, 4, chiralCheckAtomIndex, 0 ); amoebaTorsionTorsionForce->addTorsionTorsion(0, 1, 2, 3, 4, chiralCheckAtomIndex, 0);
amoebaTorsionTorsionForce->setTorsionTorsionGrid( 0, getTorsionGrid( gridIndex ) ); amoebaTorsionTorsionForce->setTorsionTorsionGrid(0, getTorsionGrid(gridIndex));
system.addForce(amoebaTorsionTorsionForce); system.addForce(amoebaTorsionTorsionForce);
Context context(system, integrator, Platform::getPlatformByName( "CUDA")); Context context(system, integrator, Platform::getPlatformByName("CUDA"));
context.setPositions(positions); context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy); State state = context.getState(State::Forces | State::Energy);
std::vector<Vec3> forces = state.getForces(); std::vector<Vec3> forces = state.getForces();
const double conversion = -1.0; const double conversion = -1.0;
for( unsigned int ii = 0; ii < forces.size(); ii++ ){ for (unsigned int ii = 0; ii < forces.size(); ii++) {
forces[ii][0] *= conversion; forces[ii][0] *= conversion;
forces[ii][1] *= conversion; forces[ii][1] *= conversion;
forces[ii][2] *= conversion; forces[ii][2] *= conversion;
} }
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaTorsionTorsionForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy() );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2],
forces[ii][0], forces[ii][1], forces[ii][2] );
}
(void) fflush( log );
}
#endif
double tolerance = 1.0e-03; double tolerance = 1.0e-03;
for( unsigned int ii = 0; ii < forces.size(); ii++ ){ for (unsigned int ii = 0; ii < forces.size(); ii++) {
ASSERT_EQUAL_VEC( expectedForces[ii], forces[ii], tolerance ); ASSERT_EQUAL_VEC(expectedForces[ii], forces[ii], tolerance);
} }
ASSERT_EQUAL_TOL( expectedEnergy, state.getPotentialEnergy(), tolerance ); ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), tolerance);
} }
...@@ -2687,8 +2674,7 @@ int main(int argc, char* argv[]) { ...@@ -2687,8 +2674,7 @@ int main(int argc, char* argv[]) {
registerAmoebaCudaKernelFactories(); registerAmoebaCudaKernelFactories();
if (argc > 1) if (argc > 1)
Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1])); Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1]));
FILE* log = NULL; testTorsionTorsion(1);
testTorsionTorsion( log, 1 );
} catch(const std::exception& e) { } catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl; std::cout << "exception: " << e.what() << std::endl;
std::cout << "FAIL - ERROR. Test failed." << std::endl; std::cout << "FAIL - ERROR. Test failed." << std::endl;
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaVdwForce.cpp
View file @ f11bafd9
This source diff could not be displayed because it is too large. You can view the blob instead.
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