Skip to content

GitLab

Commit d91ae250 authored by Mark Friedrichs's avatar Mark Friedrichs
Browse files

AmoebaReferencePiTorsionForce and ReferenceAmoebaTorsionForce

parent 967968ad
Hide whitespace changes
Inline Side-by-side
Showing with 1453 additions and 188 deletions
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernelFactory.cpp
View file @ d91ae250
...@@ -42,9 +42,9 @@ fprintf( stderr,"In registerKernelFactories AmoebaReferenceKernelFactory\n" ); f ...@@ -42,9 +42,9 @@ fprintf( stderr,"In registerKernelFactories AmoebaReferenceKernelFactory\n" ); f
platform.registerKernelFactory(CalcAmoebaHarmonicBondForceKernel::Name(), factory); platform.registerKernelFactory(CalcAmoebaHarmonicBondForceKernel::Name(), factory);
platform.registerKernelFactory(CalcAmoebaHarmonicAngleForceKernel::Name(), factory); platform.registerKernelFactory(CalcAmoebaHarmonicAngleForceKernel::Name(), factory);
platform.registerKernelFactory(CalcAmoebaHarmonicInPlaneAngleForceKernel::Name(), factory); platform.registerKernelFactory(CalcAmoebaHarmonicInPlaneAngleForceKernel::Name(), factory);
/*
platform.registerKernelFactory(CalcAmoebaTorsionForceKernel::Name(), factory); platform.registerKernelFactory(CalcAmoebaTorsionForceKernel::Name(), factory);
platform.registerKernelFactory(CalcAmoebaPiTorsionForceKernel::Name(), factory); platform.registerKernelFactory(CalcAmoebaPiTorsionForceKernel::Name(), factory);
/*
platform.registerKernelFactory(CalcAmoebaStretchBendForceKernel::Name(), factory); platform.registerKernelFactory(CalcAmoebaStretchBendForceKernel::Name(), factory);
platform.registerKernelFactory(CalcAmoebaOutOfPlaneBendForceKernel::Name(), factory); platform.registerKernelFactory(CalcAmoebaOutOfPlaneBendForceKernel::Name(), factory);
platform.registerKernelFactory(CalcAmoebaTorsionTorsionForceKernel::Name(), factory); platform.registerKernelFactory(CalcAmoebaTorsionTorsionForceKernel::Name(), factory);
...@@ -71,13 +71,13 @@ KernelImpl* AmoebaReferenceKernelFactory::createKernelImpl(std::string name, con ...@@ -71,13 +71,13 @@ KernelImpl* AmoebaReferenceKernelFactory::createKernelImpl(std::string name, con
if (name == CalcAmoebaHarmonicInPlaneAngleForceKernel::Name()) if (name == CalcAmoebaHarmonicInPlaneAngleForceKernel::Name())
return new ReferenceCalcAmoebaHarmonicInPlaneAngleForceKernel(name, platform, context.getSystem()); return new ReferenceCalcAmoebaHarmonicInPlaneAngleForceKernel(name, platform, context.getSystem());
/*
if (name == CalcAmoebaTorsionForceKernel::Name()) if (name == CalcAmoebaTorsionForceKernel::Name())
return new ReferenceCalcAmoebaTorsionForceKernel(name, platform, context.getSystem()); return new ReferenceCalcAmoebaTorsionForceKernel(name, platform, context.getSystem());
if (name == CalcAmoebaPiTorsionForceKernel::Name()) if (name == CalcAmoebaPiTorsionForceKernel::Name())
return new ReferenceCalcAmoebaPiTorsionForceKernel(name, platform, context.getSystem()); return new ReferenceCalcAmoebaPiTorsionForceKernel(name, platform, context.getSystem());
/*
if (name == CalcAmoebaStretchBendForceKernel::Name()) if (name == CalcAmoebaStretchBendForceKernel::Name())
return new ReferenceCalcAmoebaStretchBendForceKernel(name, platform, context.getSystem()); return new ReferenceCalcAmoebaStretchBendForceKernel(name, platform, context.getSystem());
......
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.cpp
View file @ d91ae250
...@@ -28,6 +28,8 @@ ...@@ -28,6 +28,8 @@
#include "AmoebaReferenceHarmonicBondForce.h" #include "AmoebaReferenceHarmonicBondForce.h"
#include "AmoebaReferenceHarmonicAngleForce.h" #include "AmoebaReferenceHarmonicAngleForce.h"
#include "AmoebaReferenceHarmonicInPlaneAngleForce.h" #include "AmoebaReferenceHarmonicInPlaneAngleForce.h"
#include "AmoebaReferenceTorsionForce.h"
#include "AmoebaReferencePiTorsionForce.h"
#include "ReferencePlatform.h" #include "ReferencePlatform.h"
#include "openmm/internal/ContextImpl.h" #include "openmm/internal/ContextImpl.h"
//#include "internal/AmoebaMultipoleForceImpl.h" //#include "internal/AmoebaMultipoleForceImpl.h"
...@@ -207,98 +209,123 @@ double ReferenceCalcAmoebaHarmonicInPlaneAngleForceKernel::execute(ContextImpl& ...@@ -207,98 +209,123 @@ double ReferenceCalcAmoebaHarmonicInPlaneAngleForceKernel::execute(ContextImpl&
return static_cast<double>(energy); return static_cast<double>(energy);
} }
//ReferenceCalcAmoebaTorsionForceKernel::ReferenceCalcAmoebaTorsionForceKernel(std::string name, const Platform& platform, System& system) : ReferenceCalcAmoebaTorsionForceKernel::ReferenceCalcAmoebaTorsionForceKernel(std::string name, const Platform& platform, System& system) :
// CalcAmoebaTorsionForceKernel(name, platform), system(system) { CalcAmoebaTorsionForceKernel(name, platform), system(system) {
// data.incrementKernelCount(); }
//}
// ReferenceCalcAmoebaTorsionForceKernel::~ReferenceCalcAmoebaTorsionForceKernel() {
//ReferenceCalcAmoebaTorsionForceKernel::~ReferenceCalcAmoebaTorsionForceKernel() { }
// data.decrementKernelCount();
//} void ReferenceCalcAmoebaTorsionForceKernel::initialize(const System& system, const AmoebaTorsionForce& force) {
//
//void ReferenceCalcAmoebaTorsionForceKernel::initialize(const System& system, const AmoebaTorsionForce& force) { numTorsions = force.getNumTorsions();
// torsionParameters1.resize( numTorsions );
// data.setAmoebaLocalForcesKernel( this ); torsionParameters2.resize( numTorsions );
// numTorsions = force.getNumTorsions(); torsionParameters3.resize( numTorsions );
// std::vector<int> particle1(numTorsions); for (int ii = 0; ii < numTorsions; ii++) {
// std::vector<int> particle2(numTorsions);
// std::vector<int> particle3(numTorsions); int particle1Index, particle2Index, particle3Index, particle4Index;
// std::vector<int> particle4(numTorsions); std::vector<double> torsionParameter1;
// std::vector<double> torsionParameter2;
// std::vector< std::vector<RealOpenMM> > torsionParameters1(numTorsions); std::vector<double> torsionParameter3;
// std::vector< std::vector<RealOpenMM> > torsionParameters2(numTorsions);
// std::vector< std::vector<RealOpenMM> > torsionParameters3(numTorsions); std::vector<RealOpenMM> torsionParameters1F(3);
// std::vector<RealOpenMM> torsionParameters2F(3);
// for (int i = 0; i < numTorsions; i++) { std::vector<RealOpenMM> torsionParameters3F(3);
//
// std::vector<double> torsionParameter1; force.getTorsionParameters(ii, particle1Index, particle2Index, particle3Index, particle4Index, torsionParameter1, torsionParameter2, torsionParameter3 );
// std::vector<double> torsionParameter2; particle1.push_back( particle1Index );
// std::vector<double> torsionParameter3; particle2.push_back( particle2Index );
// particle3.push_back( particle3Index );
// std::vector<RealOpenMM> torsionParameters1F(3); particle4.push_back( particle4Index );
// std::vector<RealOpenMM> torsionParameters2F(3); torsionParameters1[ii].resize( torsionParameter1.size() );
// std::vector<RealOpenMM> torsionParameters3F(3); torsionParameters2[ii].resize( torsionParameter2.size() );
// torsionParameters3[ii].resize( torsionParameter3.size() );
// force.getTorsionParameters(i, particle1[i], particle2[i], particle3[i], particle4[i], torsionParameter1, torsionParameter2, torsionParameter3 ); for ( unsigned int jj = 0; jj < torsionParameter1.size(); jj++) {
// for ( unsigned int jj = 0; jj < torsionParameter1.size(); jj++) { torsionParameters1[ii][jj] = static_cast<RealOpenMM>(torsionParameter1[jj]);
// torsionParameters1F[jj] = static_cast<RealOpenMM>(torsionParameter1[jj]); torsionParameters2[ii][jj] = static_cast<RealOpenMM>(torsionParameter2[jj]);
// torsionParameters2F[jj] = static_cast<RealOpenMM>(torsionParameter2[jj]); torsionParameters3[ii][jj] = static_cast<RealOpenMM>(torsionParameter3[jj]);
// torsionParameters3F[jj] = static_cast<RealOpenMM>(torsionParameter3[jj]); }
// } }
// torsionParameters1[i] = torsionParameters1F; }
// torsionParameters2[i] = torsionParameters2F;
// torsionParameters3[i] = torsionParameters3F; double ReferenceCalcAmoebaTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
// } RealOpenMM** posData = extractPositions(context);
// gpuSetAmoebaTorsionParameters(data.getAmoebaGpu(), particle1, particle2, particle3, particle4, torsionParameters1, torsionParameters2, torsionParameters3 ); RealOpenMM** forceData = extractForces(context);
// RealOpenMM energy = 0.0;
//} for (unsigned int ii = 0; ii < numTorsions; ii++) {
//
//double ReferenceCalcAmoebaTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) { int particle1Index = particle1[ii];
// if( data.getAmoebaLocalForcesKernel() == this ){ int particle2Index = particle2[ii];
// computeAmoebaLocalForces( data ); int particle3Index = particle3[ii];
// } int particle4Index = particle4[ii];
// return 0.0;
//} RealOpenMM* forces[4];
// forces[0] = forceData[particle1Index];
//ReferenceCalcAmoebaPiTorsionForceKernel::ReferenceCalcAmoebaPiTorsionForceKernel(std::string name, const Platform& platform, System& system) : forces[1] = forceData[particle2Index];
// CalcAmoebaPiTorsionForceKernel(name, platform), system(system) { forces[2] = forceData[particle3Index];
// data.incrementKernelCount(); forces[3] = forceData[particle4Index];
//}
// energy += AmoebaReferenceTorsionForce::calculateForceAndEnergy(
//ReferenceCalcAmoebaPiTorsionForceKernel::~ReferenceCalcAmoebaPiTorsionForceKernel() { posData[particle1Index], posData[particle2Index], posData[particle3Index], posData[particle4Index],
// data.decrementKernelCount(); torsionParameters1[ii], torsionParameters2[ii], torsionParameters3[ii], forces );
//} }
// return static_cast<double>(energy);
//void ReferenceCalcAmoebaPiTorsionForceKernel::initialize(const System& system, const AmoebaPiTorsionForce& force) { }
//
// data.setAmoebaLocalForcesKernel( this ); ReferenceCalcAmoebaPiTorsionForceKernel::ReferenceCalcAmoebaPiTorsionForceKernel(std::string name, const Platform& platform, System& system) :
// numPiTorsions = force.getNumPiTorsions(); CalcAmoebaPiTorsionForceKernel(name, platform), system(system) {
// }
// std::vector<int> particle1(numPiTorsions);
// std::vector<int> particle2(numPiTorsions); ReferenceCalcAmoebaPiTorsionForceKernel::~ReferenceCalcAmoebaPiTorsionForceKernel() {
// std::vector<int> particle3(numPiTorsions); }
// std::vector<int> particle4(numPiTorsions);
// std::vector<int> particle5(numPiTorsions); void ReferenceCalcAmoebaPiTorsionForceKernel::initialize(const System& system, const AmoebaPiTorsionForce& force) {
// std::vector<int> particle6(numPiTorsions);
// numPiTorsions = force.getNumPiTorsions();
// std::vector<RealOpenMM> torsionKParameters(numPiTorsions); for (int ii = 0; ii < numPiTorsions; ii++) {
//
// for (int i = 0; i < numPiTorsions; i++) { int particle1Index, particle2Index, particle3Index, particle4Index, particle5Index, particle6Index;
// double kTorsionParameter;
// double torsionKParameter; force.getPiTorsionParameters(ii, particle1Index, particle2Index, particle3Index, particle4Index, particle5Index, particle6Index, kTorsionParameter );
// particle1.push_back( particle1Index );
// force.getPiTorsionParameters(i, particle1[i], particle2[i], particle3[i], particle4[i], particle5[i], particle6[i], torsionKParameter); particle2.push_back( particle2Index );
// torsionKParameters[i] = static_cast<RealOpenMM>(torsionKParameter); particle3.push_back( particle3Index );
// } particle4.push_back( particle4Index );
// gpuSetAmoebaPiTorsionParameters(data.getAmoebaGpu(), particle1, particle2, particle3, particle4, particle5, particle6, torsionKParameters); particle5.push_back( particle5Index );
//} particle6.push_back( particle6Index );
// kTorsion.push_back( static_cast<RealOpenMM>(kTorsionParameter) );
//double ReferenceCalcAmoebaPiTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) { }
// if( data.getAmoebaLocalForcesKernel() == this ){ }
// computeAmoebaLocalForces( data );
// } double ReferenceCalcAmoebaPiTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
// return 0.0; RealOpenMM** posData = extractPositions(context);
//} RealOpenMM** forceData = extractForces(context);
RealOpenMM energy = 0.0;
for( unsigned int ii = 0; ii < numPiTorsions; ii++ ){
int particle1Index = particle1[ii];
int particle2Index = particle2[ii];
int particle3Index = particle3[ii];
int particle4Index = particle4[ii];
int particle5Index = particle5[ii];
int particle6Index = particle6[ii];
RealOpenMM* forces[6];
forces[0] = forceData[particle1Index];
forces[1] = forceData[particle2Index];
forces[2] = forceData[particle3Index];
forces[3] = forceData[particle4Index];
forces[4] = forceData[particle5Index];
forces[5] = forceData[particle6Index];
energy += AmoebaReferencePiTorsionForce::calculateForceAndEnergy(
posData[particle1Index], posData[particle2Index], posData[particle3Index], posData[particle4Index],
posData[particle5Index], posData[particle6Index], kTorsion[ii], forces );
}
return static_cast<double>(energy);
}
// //
//ReferenceCalcAmoebaStretchBendForceKernel::ReferenceCalcAmoebaStretchBendForceKernel(std::string name, const Platform& platform, System& system) : //ReferenceCalcAmoebaStretchBendForceKernel::ReferenceCalcAmoebaStretchBendForceKernel(std::string name, const Platform& platform, System& system) :
// CalcAmoebaStretchBendForceKernel(name, platform), system(system) { // CalcAmoebaStretchBendForceKernel(name, platform), system(system) {
......
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.h
View file @ d91ae250
...@@ -144,62 +144,76 @@ private: ...@@ -144,62 +144,76 @@ private:
System& system; System& system;
}; };
// /** /**
// * This kernel is invoked by AmoebaTorsionForce to calculate the forces acting on the system and the energy of the system. * This kernel is invoked by AmoebaTorsionForce to calculate the forces acting on the system and the energy of the system.
// */ */
// class ReferenceCalcAmoebaTorsionForceKernel : public CalcAmoebaTorsionForceKernel { class ReferenceCalcAmoebaTorsionForceKernel : public CalcAmoebaTorsionForceKernel {
// public: public:
// ReferenceCalcAmoebaTorsionForceKernel(std::string name, const Platform& platform, System& system); ReferenceCalcAmoebaTorsionForceKernel(std::string name, const Platform& platform, System& system);
// ~ReferenceCalcAmoebaTorsionForceKernel(); ~ReferenceCalcAmoebaTorsionForceKernel();
// /** /**
// * Initialize the kernel. * Initialize the kernel.
// * *
// * @param system the System this kernel will be applied to * @param system the System this kernel will be applied to
// * @param force the AmoebaTorsionForce this kernel will be used for * @param force the AmoebaTorsionForce this kernel will be used for
// */ */
// void initialize(const System& system, const AmoebaTorsionForce& force); void initialize(const System& system, const AmoebaTorsionForce& force);
// /** /**
// * Execute the kernel to calculate the forces and/or energy. * Execute the kernel to calculate the forces and/or energy.
// * *
// * @param context the context in which to execute this kernel * @param context the context in which to execute this kernel
// * @param includeForces true if forces should be calculated * @param includeForces true if forces should be calculated
// * @param includeEnergy true if the energy should be calculated * @param includeEnergy true if the energy should be calculated
// * @return the potential energy due to the force * @return the potential energy due to the force
// */ */
// double execute(ContextImpl& context, bool includeForces, bool includeEnergy); double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
// private: private:
// int numTorsions; int numTorsions;
// System& system; std::vector<int> particle1;
// }; std::vector<int> particle2;
// std::vector<int> particle3;
// /** std::vector<int> particle4;
// * This kernel is invoked by AmoebaPiTorsionForce to calculate the forces acting on the system and the energy of the system. std::vector< std::vector<RealOpenMM> > torsionParameters1;
// */ std::vector< std::vector<RealOpenMM> > torsionParameters2;
// class ReferenceCalcAmoebaPiTorsionForceKernel : public CalcAmoebaPiTorsionForceKernel { std::vector< std::vector<RealOpenMM> > torsionParameters3;
// public: System& system;
// ReferenceCalcAmoebaPiTorsionForceKernel(std::string name, const Platform& platform, System& system); };
// ~ReferenceCalcAmoebaPiTorsionForceKernel();
// /** /**
// * Initialize the kernel. * This kernel is invoked by AmoebaPiTorsionForce to calculate the forces acting on the system and the energy of the system.
// * */
// * @param system the System this kernel will be applied to class ReferenceCalcAmoebaPiTorsionForceKernel : public CalcAmoebaPiTorsionForceKernel {
// * @param force the AmoebaPiTorsionForce this kernel will be used for public:
// */ ReferenceCalcAmoebaPiTorsionForceKernel(std::string name, const Platform& platform, System& system);
// void initialize(const System& system, const AmoebaPiTorsionForce& force); ~ReferenceCalcAmoebaPiTorsionForceKernel();
// /** /**
// * Execute the kernel to calculate the forces and/or energy. * Initialize the kernel.
// * *
// * @param context the context in which to execute this kernel * @param system the System this kernel will be applied to
// * @param includeForces true if forces should be calculated * @param force the AmoebaPiTorsionForce this kernel will be used for
// * @param includeEnergy true if the energy should be calculated */
// * @return the potential energy due to the force void initialize(const System& system, const AmoebaPiTorsionForce& force);
// */ /**
// double execute(ContextImpl& context, bool includeForces, bool includeEnergy); * Execute the kernel to calculate the forces and/or energy.
// private: *
// int numPiTorsions; * @param context the context in which to execute this kernel
// System& system; * @param includeForces true if forces should be calculated
// }; * @param includeEnergy true if the energy should be calculated
// * @return the potential energy due to the force
*/
double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
private:
int numPiTorsions;
std::vector<int> particle1;
std::vector<int> particle2;
std::vector<int> particle3;
std::vector<int> particle4;
std::vector<int> particle5;
std::vector<int> particle6;
std::vector<RealOpenMM> kTorsion;
System& system;
};
// /** // /**
// * This kernel is invoked by AmoebaStretchBendForce to calculate the forces acting on the system and the energy of the system. // * This kernel is invoked by AmoebaStretchBendForce to calculate the forces acting on the system and the energy of the system.
// */ // */
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceForce.h
View file @ d91ae250
...@@ -22,8 +22,8 @@ ...@@ -22,8 +22,8 @@
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/ */
#ifndef __AmoebaReferenceHarmonicBondForce_H__ #ifndef __AmoebaReferenceForce_H__
#define __AmoebaReferenceHarmonicBondForce_H__ #define __AmoebaReferenceForce_H__
#include "SimTKUtilities/SimTKOpenMMRealType.h" #include "SimTKUtilities/SimTKOpenMMRealType.h"
#include <vector> #include <vector>
...@@ -118,4 +118,4 @@ public: ...@@ -118,4 +118,4 @@ public:
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
#endif // _AmoebaReferenceHarmonicBondForce___ #endif // _AmoebaReferenceForce___
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceHarmonicAngleForce.cpp
View file @ d91ae250
...@@ -111,8 +111,8 @@ RealOpenMM AmoebaReferenceHarmonicAngleForce::getPrefactorsGivenAngleCosine( Rea ...@@ -111,8 +111,8 @@ RealOpenMM AmoebaReferenceHarmonicAngleForce::getPrefactorsGivenAngleCosine( Rea
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
RealOpenMM AmoebaReferenceHarmonicAngleForce::calculateForceAndEnergy( RealOpenMM* positionAtomA, RealOpenMM* positionAtomB, RealOpenMM AmoebaReferenceHarmonicAngleForce::calculateForceAndEnergy( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
RealOpenMM* positionAtomC, const RealOpenMM* positionAtomC,
RealOpenMM angle, RealOpenMM angleK, RealOpenMM angle, RealOpenMM angleK,
RealOpenMM angleCubic, RealOpenMM angleQuartic, RealOpenMM angleCubic, RealOpenMM angleQuartic,
RealOpenMM anglePentic, RealOpenMM angleSextic, RealOpenMM anglePentic, RealOpenMM angleSextic,
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceHarmonicAngleForce.h
View file @ d91ae250
...@@ -91,8 +91,8 @@ public: ...@@ -91,8 +91,8 @@ public:
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
static RealOpenMM calculateForceAndEnergy( RealOpenMM* positionAtomA, RealOpenMM* positionAtomB, static RealOpenMM calculateForceAndEnergy( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
RealOpenMM* positionAtomC, const RealOpenMM* positionAtomC,
RealOpenMM angle, RealOpenMM angleK, RealOpenMM angle, RealOpenMM angleK,
RealOpenMM angleCubic, RealOpenMM angleQuartic, RealOpenMM angleCubic, RealOpenMM angleQuartic,
RealOpenMM anglePentic, RealOpenMM angleSextic, RealOpenMM anglePentic, RealOpenMM angleSextic,
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceHarmonicBondForce.cpp
View file @ d91ae250
...@@ -23,6 +23,7 @@ ...@@ -23,6 +23,7 @@
*/ */
#include "AmoebaReferenceHarmonicBondForce.h" #include "AmoebaReferenceHarmonicBondForce.h"
#include "AmoebaReferenceForce.h"
#include <vector> #include <vector>
/**--------------------------------------------------------------------------------------- /**---------------------------------------------------------------------------------------
...@@ -41,7 +42,7 @@ ...@@ -41,7 +42,7 @@
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
RealOpenMM AmoebaReferenceHarmonicBondForce::calculateForceAndEnergy( RealOpenMM* positionAtomA, RealOpenMM* positionAtomB, RealOpenMM AmoebaReferenceHarmonicBondForce::calculateForceAndEnergy( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
RealOpenMM bondLength, RealOpenMM bondK, RealOpenMM bondLength, RealOpenMM bondK,
RealOpenMM bondCubic, RealOpenMM bondQuartic, RealOpenMM bondCubic, RealOpenMM bondQuartic,
RealOpenMM** forces ){ RealOpenMM** forces ){
...@@ -60,11 +61,8 @@ RealOpenMM AmoebaReferenceHarmonicBondForce::calculateForceAndEnergy( RealOpenMM ...@@ -60,11 +61,8 @@ RealOpenMM AmoebaReferenceHarmonicBondForce::calculateForceAndEnergy( RealOpenMM
// get deltaR, R2, and R between 2 atoms // get deltaR, R2, and R between 2 atoms
std::vector<RealOpenMM> deltaR; std::vector<RealOpenMM> deltaR;
deltaR.push_back( positionAtomB[0] - positionAtomA[0] ); AmoebaReferenceForce::loadDeltaR( positionAtomA, positionAtomB, deltaR );
deltaR.push_back( positionAtomB[1] - positionAtomA[1] ); RealOpenMM r = AmoebaReferenceForce::getNorm3( deltaR );
deltaR.push_back( positionAtomB[2] - positionAtomA[2] );
RealOpenMM r = SQRT( deltaR[0]*deltaR[0] + deltaR[1]*deltaR[1] + deltaR[2]*deltaR[2] );
// deltaIdeal = r - r_0 // deltaIdeal = r - r_0
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceHarmonicBondForce.h
View file @ d91ae250
...@@ -31,25 +31,23 @@ ...@@ -31,25 +31,23 @@
class AmoebaReferenceHarmonicBondForce { class AmoebaReferenceHarmonicBondForce {
private: public:
public: /**---------------------------------------------------------------------------------------
/**---------------------------------------------------------------------------------------
Constructor Constructor
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
AmoebaReferenceHarmonicBondForce( ){}; AmoebaReferenceHarmonicBondForce( ){};
/**--------------------------------------------------------------------------------------- /**---------------------------------------------------------------------------------------
Destructor Destructor
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
~AmoebaReferenceHarmonicBondForce( ); ~AmoebaReferenceHarmonicBondForce( );
/**--------------------------------------------------------------------------------------- /**---------------------------------------------------------------------------------------
...@@ -68,7 +66,7 @@ class AmoebaReferenceHarmonicBondForce { ...@@ -68,7 +66,7 @@ class AmoebaReferenceHarmonicBondForce {
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
static RealOpenMM calculateForceAndEnergy( RealOpenMM* positionAtomA, RealOpenMM* positionAtomB, static RealOpenMM calculateForceAndEnergy( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
RealOpenMM bondLength, RealOpenMM bondK, RealOpenMM bondLength, RealOpenMM bondK,
RealOpenMM bondCubic, RealOpenMM bondQuartic, RealOpenMM bondCubic, RealOpenMM bondQuartic,
RealOpenMM** forces ); RealOpenMM** forces );
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceHarmonicInPlaneAngleForce.cpp
View file @ d91ae250
...@@ -112,11 +112,11 @@ RealOpenMM AmoebaReferenceHarmonicInPlaneAngleForce::getPrefactorsGivenAngleCosi ...@@ -112,11 +112,11 @@ RealOpenMM AmoebaReferenceHarmonicInPlaneAngleForce::getPrefactorsGivenAngleCosi
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
RealOpenMM AmoebaReferenceHarmonicInPlaneAngleForce::calculateForceAndEnergy( RealOpenMM* positionAtomA, RealOpenMM* positionAtomB, RealOpenMM AmoebaReferenceHarmonicInPlaneAngleForce::calculateForceAndEnergy( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
RealOpenMM* positionAtomC, RealOpenMM* positionAtomD, const RealOpenMM* positionAtomC, const RealOpenMM* positionAtomD,
RealOpenMM angle, RealOpenMM angleK, RealOpenMM angle, RealOpenMM angleK,
RealOpenMM angleCubic, RealOpenMM angleQuartic, RealOpenMM angleCubic, RealOpenMM angleQuartic,
RealOpenMM anglePentic, RealOpenMM angleSextic, RealOpenMM anglePentic, RealOpenMM angleSextic,
RealOpenMM** forces ){ RealOpenMM** forces ){
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -134,6 +134,7 @@ RealOpenMM AmoebaReferenceHarmonicInPlaneAngleForce::calculateForceAndEnergy( Re ...@@ -134,6 +134,7 @@ RealOpenMM AmoebaReferenceHarmonicInPlaneAngleForce::calculateForceAndEnergy( Re
// AP = A - P // AP = A - P
// CP = A - P // CP = A - P
// M = CP x AP // M = CP x AP
enum { AD, BD, CD, T, AP, P, CP, M, APxM, CPxM, ADxBD, BDxCD, TxCD, ADxT, dBxAD, CDxdB, LastDeltaAtomIndex }; enum { AD, BD, CD, T, AP, P, CP, M, APxM, CPxM, ADxBD, BDxCD, TxCD, ADxT, dBxAD, CDxdB, LastDeltaAtomIndex };
std::vector<RealOpenMM> deltaR[LastDeltaAtomIndex]; std::vector<RealOpenMM> deltaR[LastDeltaAtomIndex];
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceHarmonicInPlaneAngleForce.h
View file @ d91ae250
...@@ -80,20 +80,20 @@ public: ...@@ -80,20 +80,20 @@ public:
@param positionAtomB Cartesian coordinates of atom B @param positionAtomB Cartesian coordinates of atom B
@param positionAtomC Cartesian coordinates of atom C @param positionAtomC Cartesian coordinates of atom C
@param positionAtomD Cartesian coordinates of atom D @param positionAtomD Cartesian coordinates of atom D
@param angleLength angle @param angle angle
@param angleK quadratic angle force @param angleK quadratic angle force parameter
@param angleCubic cubic angle force parameter @param angleCubic cubic angle force parameter
@param angleQuartic quartic angle force parameter @param angleQuartic quartic angle force parameter
@param anglePentic pentic angle force parameter @param anglePentic pentic angle force parameter
@param angleSextic sextic angle force parameter @param angleSextic sextic angle force parameter
@param forces force vector @param forces force vector
@return energy @return energy
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
static RealOpenMM calculateForceAndEnergy( RealOpenMM* positionAtomA, RealOpenMM* positionAtomB, static RealOpenMM calculateForceAndEnergy( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
RealOpenMM* positionAtomC, RealOpenMM* positionAtomD, const RealOpenMM* positionAtomC, const RealOpenMM* positionAtomD,
RealOpenMM angle, RealOpenMM angleK, RealOpenMM angle, RealOpenMM angleK,
RealOpenMM angleCubic, RealOpenMM angleQuartic, RealOpenMM angleCubic, RealOpenMM angleQuartic,
RealOpenMM anglePentic, RealOpenMM angleSextic, RealOpenMM anglePentic, RealOpenMM angleSextic,
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferencePiTorsionForce.cpp 0 → 100644
View file @ d91ae250
/* Portions copyright (c) 2006 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject
* to the following conditions:
*
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
* LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
* OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#include "AmoebaReferenceForce.h"
#include "AmoebaReferencePiTorsionForce.h"
#include <vector>
/**---------------------------------------------------------------------------------------
Calculate Amoeba pi-torsion ixn (force and energy)
@param positionAtomA Cartesian coordinates of atom A
@param positionAtomB Cartesian coordinates of atom B
@param positionAtomC Cartesian coordinates of atom C
@param positionAtomD Cartesian coordinates of atom D
@param positionAtomE Cartesian coordinates of atom E
@param positionAtomF Cartesian coordinates of atom F
@param piTorsionK force constant
@param forces force vector
@return energy
--------------------------------------------------------------------------------------- */
RealOpenMM AmoebaReferencePiTorsionForce::calculateForceAndEnergy( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
const RealOpenMM* positionAtomC, const RealOpenMM* positionAtomD,
const RealOpenMM* positionAtomE, const RealOpenMM* positionAtomF,
RealOpenMM piTorsionK, RealOpenMM** forces ){
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "AmoebaReferencePiTorsionForce::calculateHarmonicForce";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM two = 2.0;
// ---------------------------------------------------------------------------------------
enum { AD, BD, EC, FC, P, Q, CP, DC, QD, T, U, TU, DP, QC, dT, dU, dP, dQ, dC1, dC2, dD1, dD2, LastDeltaIndex };
std::vector<RealOpenMM> deltaR[LastDeltaIndex];
for( unsigned int ii = 0; ii < LastDeltaIndex; ii++ ){
deltaR[ii].resize(3);
}
AmoebaReferenceForce::loadDeltaR( positionAtomD, positionAtomA, deltaR[AD] );
AmoebaReferenceForce::loadDeltaR( positionAtomD, positionAtomB, deltaR[BD] );
AmoebaReferenceForce::loadDeltaR( positionAtomC, positionAtomE, deltaR[EC] );
AmoebaReferenceForce::loadDeltaR( positionAtomC, positionAtomF, deltaR[FC] );
enum { A, B, C, D, E, F, LastAtomIndex };
std::vector<RealOpenMM> d[LastAtomIndex];
for( unsigned int ii = 0; ii < LastAtomIndex; ii++ ){
d[ii].resize(3);
}
AmoebaReferenceForce::getCrossProduct( deltaR[AD], deltaR[BD], deltaR[P] );
AmoebaReferenceForce::getCrossProduct( deltaR[EC], deltaR[FC], deltaR[Q] );
for( int ii = 0; ii < 3; ii++ ){
deltaR[CP][ii] = -deltaR[P][ii];
deltaR[DC][ii] = positionAtomD[ii] - positionAtomC[ii];
deltaR[QD][ii] = deltaR[Q][ii];
deltaR[P][ii] += positionAtomC[ii];
deltaR[Q][ii] += positionAtomD[ii];
}
AmoebaReferenceForce::getCrossProduct( deltaR[CP], deltaR[DC], deltaR[T] );
AmoebaReferenceForce::getCrossProduct( deltaR[DC], deltaR[QD], deltaR[U] );
AmoebaReferenceForce::getCrossProduct( deltaR[T], deltaR[U], deltaR[TU] );
RealOpenMM rT2 = AmoebaReferenceForce::getNormSquared3( deltaR[T] );
RealOpenMM rU2 = AmoebaReferenceForce::getNormSquared3( deltaR[U] );
RealOpenMM rTrU = SQRT( rT2*rU2 );
if( rTrU <= zero ){
return zero;
}
RealOpenMM rDC = AmoebaReferenceForce::getNorm3( deltaR[DC] );
RealOpenMM cosine = AmoebaReferenceForce::getDotProduct3( deltaR[T], deltaR[U] );
cosine /= rTrU;
RealOpenMM sine = AmoebaReferenceForce::getDotProduct3( deltaR[DC], deltaR[TU] );
sine /= ( rDC*rTrU );
RealOpenMM cosine2 = cosine*cosine - sine*sine;
RealOpenMM sine2 = two*cosine*sine;
RealOpenMM phi2 = one - cosine2;
RealOpenMM dphi2 = two*sine2;
RealOpenMM dedphi = piTorsionK*dphi2;
for( unsigned int ii = 0; ii < 3; ii++ ){
deltaR[DP][ii] = positionAtomD[ii] - deltaR[P][ii];
deltaR[QC][ii] = deltaR[Q][ii] - positionAtomC[ii];
}
RealOpenMM factorT = dedphi/( rDC*rT2 );
RealOpenMM factorU = -dedphi/( rDC*rU2 );
AmoebaReferenceForce::getCrossProduct( deltaR[T], deltaR[DC], deltaR[dT] );
AmoebaReferenceForce::getCrossProduct( deltaR[U], deltaR[DC], deltaR[dU] );
for( int ii = 0; ii < 3; ii++ ){
deltaR[dT][ii] *= factorT;
deltaR[dU][ii] *= factorU;
}
AmoebaReferenceForce::getCrossProduct( deltaR[dT], deltaR[DC], deltaR[dP] );
AmoebaReferenceForce::getCrossProduct( deltaR[dU], deltaR[DC], deltaR[dQ] );
AmoebaReferenceForce::getCrossProduct( deltaR[DP], deltaR[dT], deltaR[dC1] );
AmoebaReferenceForce::getCrossProduct( deltaR[dU], deltaR[QD], deltaR[dC2] );
AmoebaReferenceForce::getCrossProduct( deltaR[dT], deltaR[CP], deltaR[dD1] );
AmoebaReferenceForce::getCrossProduct( deltaR[QC], deltaR[dU], deltaR[dD2] );
AmoebaReferenceForce::getCrossProduct( deltaR[BD], deltaR[dP], d[A] );
AmoebaReferenceForce::getCrossProduct( deltaR[dP], deltaR[AD], d[B] );
AmoebaReferenceForce::getCrossProduct( deltaR[FC], deltaR[dQ], d[E] );
AmoebaReferenceForce::getCrossProduct( deltaR[dQ], deltaR[EC], d[F] );
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[D][ii] = deltaR[dD1][ii] + deltaR[dD2][ii] + deltaR[dQ][ii] - d[A][ii] - d[B][ii];
}
// ---------------------------------------------------------------------------------------
// add in forces
forces[0][0] -= d[0][0];
forces[0][1] -= d[0][1];
forces[0][2] -= d[0][2];
forces[1][0] -= d[1][0];
forces[1][1] -= d[1][1];
forces[1][2] -= d[1][2];
forces[2][0] -= d[2][0];
forces[2][1] -= d[2][1];
forces[2][2] -= d[2][2];
forces[3][0] -= d[3][0];
forces[3][1] -= d[3][1];
forces[3][2] -= d[3][2];
forces[4][0] -= d[4][0];
forces[4][1] -= d[4][1];
forces[4][2] -= d[4][2];
forces[5][0] -= d[5][0];
forces[5][1] -= d[5][1];
forces[5][2] -= d[5][2];
// ---------------------------------------------------------------------------------------
// calculate energy if 'energy' is set
return piTorsionK*phi2;
}
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferencePiTorsionForce.h 0 → 100644
View file @ d91ae250
/* Portions copyright (c) 2006 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject
* to the following conditions:
*
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
* LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
* OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#ifndef __AmoebaReferencePiTorsionForce_H__
#define __AmoebaReferencePiTorsionForce_H__
#include "SimTKUtilities/SimTKOpenMMRealType.h"
// ---------------------------------------------------------------------------------------
class AmoebaReferencePiTorsionForce {
public:
/**---------------------------------------------------------------------------------------
Constructor
--------------------------------------------------------------------------------------- */
AmoebaReferencePiTorsionForce( );
/**---------------------------------------------------------------------------------------
Destructor
--------------------------------------------------------------------------------------- */
~AmoebaReferencePiTorsionForce( );
/**---------------------------------------------------------------------------------------
Calculate Amoeba harmonic angle ixn (force and energy)
@param positionAtomA Cartesian coordinates of atom A
@param positionAtomB Cartesian coordinates of atom B
@param positionAtomC Cartesian coordinates of atom C
@param positionAtomD Cartesian coordinates of atom D
@param positionAtomE Cartesian coordinates of atom E
@param positionAtomF Cartesian coordinates of atom F
@param kTorsion k-torsion parameter
@param forces force vector
@return energy
--------------------------------------------------------------------------------------- */
static RealOpenMM calculateForceAndEnergy( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
const RealOpenMM* positionAtomC, const RealOpenMM* positionAtomD,
const RealOpenMM* positionAtomE, const RealOpenMM* positionAtomF,
RealOpenMM kTorsion, RealOpenMM** forces );
};
// ---------------------------------------------------------------------------------------
#endif // _AmoebaReferencePiTorsionForce___
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceTorsionForce.cpp 0 → 100644
View file @ d91ae250
/* Portions copyright (c) 2006 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject
* to the following conditions:
*
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
* LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
* OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#include "AmoebaReferenceForce.h"
#include "AmoebaReferenceTorsionForce.h"
#include <vector>
/**---------------------------------------------------------------------------------------
Calculate Amoeba torsion ixn (force and energy)
@param positionAtomA Cartesian coordinates of atom A
@param positionAtomB Cartesian coordinates of atom B
@param positionAtomC Cartesian coordinates of atom C
@param positionAtomD Cartesian coordinates of atom D
@param torsion1 vector of torsion params for first index (amplitude, phase, fold)
@param torsion2 vector of torsion params for second index (amplitude, phase, fold)
@param torsion3 vector of torsion params for third index (amplitude, phase, fold)
@param forces force vector
@return energy
--------------------------------------------------------------------------------------- */
RealOpenMM AmoebaReferenceTorsionForce::calculateForceAndEnergy( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
const RealOpenMM* positionAtomC, const RealOpenMM* positionAtomD,
const std::vector<RealOpenMM>& torsionParameters1,
const std::vector<RealOpenMM>& torsionParameters2,
const std::vector<RealOpenMM>& torsionParameters3,
RealOpenMM** forces ){
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "AmoebaReferenceTorsionForce::calculateHarmonicForce";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM two = 2.0;
static const RealOpenMM three = 3.0;
// ---------------------------------------------------------------------------------------
enum { BA, CB, DC, CA, DB, LastDeltaIndex };
std::vector<RealOpenMM> deltaR[LastDeltaIndex];
AmoebaReferenceForce::loadDeltaR( positionAtomA, positionAtomB, deltaR[BA] );
AmoebaReferenceForce::loadDeltaR( positionAtomB, positionAtomC, deltaR[CB] );
AmoebaReferenceForce::loadDeltaR( positionAtomC, positionAtomD, deltaR[DC] );
AmoebaReferenceForce::loadDeltaR( positionAtomA, positionAtomC, deltaR[CA] );
AmoebaReferenceForce::loadDeltaR( positionAtomB, positionAtomD, deltaR[DB] );
enum { Xt, Xu, Xtu, LastXtIndex };
std::vector<RealOpenMM> crossProducts[LastXtIndex];
for( unsigned int ii = 0; ii < LastXtIndex; ii++ ){
crossProducts[ii].resize(3);
}
AmoebaReferenceForce::getCrossProduct( deltaR[BA], deltaR[CB], crossProducts[Xt] );
AmoebaReferenceForce::getCrossProduct( deltaR[CB], deltaR[DC], crossProducts[Xu] );
AmoebaReferenceForce::getCrossProduct( crossProducts[Xt], crossProducts[Xu], crossProducts[Xtu] );
RealOpenMM rT2 = AmoebaReferenceForce::getNormSquared3( crossProducts[Xt] );
RealOpenMM rU2 = AmoebaReferenceForce::getNormSquared3( crossProducts[Xu] );
RealOpenMM rTrU = SQRT( rT2*rU2 );
if( rTrU <= zero ){
return zero;
}
RealOpenMM rCB = AmoebaReferenceForce::getNorm3( deltaR[CB] );
// ---------------------------------------------------------------------------------------
// cos(w), cos(2w), cos(3w), ...
// sin(w), sin(2w), sin(3w), ...
RealOpenMM cosine = AmoebaReferenceForce::getDotProduct3( crossProducts[Xt], crossProducts[Xu] );
cosine /= rTrU;
RealOpenMM sine = AmoebaReferenceForce::getDotProduct3( deltaR[CB], crossProducts[Xtu] );
sine /= (rCB*rTrU);
RealOpenMM cosine2 = cosine*cosine - sine*sine;
RealOpenMM sine2 = two*cosine*sine;
RealOpenMM cosine3 = cosine*cosine2 - sine*sine2;
RealOpenMM sine3 = cosine*sine2 + sine*cosine2;
/*
RealOpenMM cosine4 = cosine*cosine3 - sine*sine3;
RealOpenMM sine4 = cosine*sine3 + sine*cosine3;
RealOpenMM cosine5 = cosine*cosine4 - sine*sine4;
RealOpenMM sine5 = cosine*sine4 + sine*cosine4;
RealOpenMM cosine6 = cosine*cosine5 - sine*sine5;
RealOpenMM sine6 = cosine*sine5 + sine*cosine5;
*/
// ---------------------------------------------------------------------------------------
// dEdPhi prefactor
RealOpenMM dEdPhi = torsionParameters1[0]* (cosine* sin( torsionParameters1[1] ) - sine* cos( torsionParameters1[1] ) );
dEdPhi += torsionParameters2[0]*two* (cosine2*sin( torsionParameters2[1] ) - sine2*cos( torsionParameters2[1] ) );
dEdPhi += torsionParameters3[0]*three*(cosine3*sin( torsionParameters3[1] ) - sine3*cos( torsionParameters3[1] ) );
// ---------------------------------------------------------------------------------------
// dEdtu[0] = dEdT
// dEdtu[1] = dEdU
// tempVector[0] == dEdA: dEdT x CB
// tempVector[1] == dEdB: (CA x dEdT) + (dEdU x DC)
// tempVector[2] == dEdC: (dEdT x BA) + (DB x dEdU)
// tempVector[3] == dEdD: (dEdU x CB)
std::vector<RealOpenMM> dEdtu[2];
std::vector<RealOpenMM> tempVector[6];
for( unsigned int ii = 0; ii < 2; ii++ ){
dEdtu[ii].resize(3);
}
for( unsigned int ii = 0; ii < 6; ii++ ){
tempVector[ii].resize(3);
}
// dEdT & dEdU
AmoebaReferenceForce::getCrossProduct( crossProducts[Xt], deltaR[CB], tempVector[0] );
AmoebaReferenceForce::getCrossProduct( crossProducts[Xu], deltaR[CB], tempVector[1] );
RealOpenMM norm[2] = { dEdPhi/(rT2*rCB ), -dEdPhi/(rU2*rCB ) };
dEdtu[0][0] = norm[0]*tempVector[0][0];
dEdtu[0][1] = norm[0]*tempVector[0][1];
dEdtu[0][2] = norm[0]*tempVector[0][2];
dEdtu[1][0] = norm[1]*tempVector[1][0];
dEdtu[1][1] = norm[1]*tempVector[1][1];
dEdtu[1][2] = norm[1]*tempVector[1][2];
// dEdA
AmoebaReferenceForce::getCrossProduct( dEdtu[0], deltaR[CB], tempVector[0] );
// dEdB
AmoebaReferenceForce::getCrossProduct( deltaR[CA], dEdtu[0], tempVector[4] );
AmoebaReferenceForce::getCrossProduct( dEdtu[1], deltaR[DC], tempVector[1] );
// dEdC
AmoebaReferenceForce::getCrossProduct( dEdtu[0], deltaR[BA], tempVector[5] );
AmoebaReferenceForce::getCrossProduct( deltaR[DB], dEdtu[1], tempVector[2] );
tempVector[1][0] += tempVector[4][0];
tempVector[2][0] += tempVector[5][0];
tempVector[1][1] += tempVector[4][1];
tempVector[2][1] += tempVector[5][1];
tempVector[1][2] += tempVector[4][2];
tempVector[2][2] += tempVector[5][2];
// dEdD
AmoebaReferenceForce::getCrossProduct( dEdtu[1], deltaR[CB], tempVector[3] );
// ---------------------------------------------------------------------------------------
// forces
forces[0][0] -= tempVector[0][0];
forces[0][1] -= tempVector[0][1];
forces[0][2] -= tempVector[0][2];
forces[1][0] -= tempVector[1][0];
forces[1][1] -= tempVector[1][1];
forces[1][2] -= tempVector[1][2];
forces[2][0] -= tempVector[2][0];
forces[2][1] -= tempVector[2][1];
forces[2][2] -= tempVector[2][2];
forces[3][0] -= tempVector[3][0];
forces[3][1] -= tempVector[3][1];
forces[3][2] -= tempVector[3][2];
// ---------------------------------------------------------------------------------------
// calculate energy
RealOpenMM energy = torsionParameters1[0]*( one + cosine* cos( torsionParameters1[1] ) + sine *sin( torsionParameters1[1] ) );
energy += torsionParameters2[0]*( one + cosine2*cos( torsionParameters2[1] ) + sine2*sin( torsionParameters2[1] ) );
energy += torsionParameters3[0]*( one + cosine3*cos( torsionParameters3[1] ) + sine3*sin( torsionParameters3[1] ) );
return energy;
}
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceTorsionForce.h 0 → 100644
View file @ d91ae250
/* Portions copyright (c) 2006 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject
* to the following conditions:
*
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
* LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
* OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#ifndef __AmoebaReferenceTorsionForce_H__
#define __AmoebaReferenceTorsionForce_H__
#include "SimTKUtilities/SimTKOpenMMRealType.h"
// ---------------------------------------------------------------------------------------
class AmoebaReferenceTorsionForce {
public:
/**---------------------------------------------------------------------------------------
Constructor
--------------------------------------------------------------------------------------- */
AmoebaReferenceTorsionForce( );
/**---------------------------------------------------------------------------------------
Destructor
--------------------------------------------------------------------------------------- */
~AmoebaReferenceTorsionForce( );
/**---------------------------------------------------------------------------------------
Calculate Amoeba harmonic angle ixn (force and energy)
@param positionAtomA Cartesian coordinates of atom A
@param positionAtomB Cartesian coordinates of atom B
@param positionAtomC Cartesian coordinates of atom C
@param positionAtomD Cartesian coordinates of atom D
@param torsion1 vector of torsion params for first index (amplitude, phase, fold)
@param torsion2 vector of torsion params for second index (amplitude, phase, fold)
@param torsion3 vector of torsion params for third index (amplitude, phase, fold)
@param forces force vector
@return energy
--------------------------------------------------------------------------------------- */
static RealOpenMM calculateForceAndEnergy( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
const RealOpenMM* positionAtomC, const RealOpenMM* positionAtomD,
const std::vector<RealOpenMM>& torsionParameters1,
const std::vector<RealOpenMM>& torsionParameters2,
const std::vector<RealOpenMM>& torsionParameters3,
RealOpenMM** forces );
};
// ---------------------------------------------------------------------------------------
#endif // _AmoebaReferenceTorsionForce___
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaPiTorsionForce.cpp 0 → 100644
View file @ d91ae250
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the Reference implementation of ReferenceAmoebaPiTorsionForce.
*/
#include "../../../tests/AssertionUtilities.h"
//#include "AmoebaTinkerParameterFile.h"
#include "openmm/Context.h"
#include "AmoebaOpenMM.h"
#include "openmm/System.h"
#include "openmm/LangevinIntegrator.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
const double TOL = 1e-5;
#define PI_M 3.141592653589
#define RADIAN 57.29577951308
/* ---------------------------------------------------------------------------------------
Compute cross product of two 3-vectors and place in 3rd vector
vectorZ = vectorX x vectorY
@param vectorX x-vector
@param vectorY y-vector
@param vectorZ z-vector
@return vector is vectorZ
--------------------------------------------------------------------------------------- */
static void crossProductVector3( double* vectorX, double* vectorY, double* vectorZ ){
vectorZ[0] = vectorX[1]*vectorY[2] - vectorX[2]*vectorY[1];
vectorZ[1] = vectorX[2]*vectorY[0] - vectorX[0]*vectorY[2];
vectorZ[2] = vectorX[0]*vectorY[1] - vectorX[1]*vectorY[0];
return;
}
static double dotVector3( double* vectorX, double* vectorY ){
return vectorX[0]*vectorY[0] + vectorX[1]*vectorY[1] + vectorX[2]*vectorY[2];
}
static void computeAmoebaPiTorsionForce(int bondIndex, std::vector<Vec3>& positions, AmoebaPiTorsionForce& amoebaPiTorsionForce,
std::vector<Vec3>& forces, double* energy, FILE* log ) {
int particle1, particle2, particle3, particle4, particle5, particle6;
double kTorsion;
amoebaPiTorsionForce.getPiTorsionParameters(bondIndex, particle1, particle2, particle3, particle4, particle5, particle6, kTorsion);
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 );
}
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];
enum { A, B, C, D, E, F, LastAtomIndex };
double d[LastAtomIndex][3];
for( int ii = 0; ii < 3; ii++ ){
deltaR[AD][ii] = positions[particle1][ii] - positions[particle4][ii];
deltaR[BD][ii] = positions[particle2][ii] - positions[particle4][ii];
deltaR[EC][ii] = positions[particle5][ii] - positions[particle3][ii];
deltaR[FC][ii] = positions[particle6][ii] - positions[particle3][ii];
}
crossProductVector3( deltaR[AD], deltaR[BD], deltaR[P] );
crossProductVector3( deltaR[EC], deltaR[FC], deltaR[Q] );
for( int ii = 0; ii < 3; ii++ ){
deltaR[CP][ii] = -deltaR[P][ii];
deltaR[DC][ii] = positions[particle4][ii] - positions[particle3][ii];
deltaR[QD][ii] = deltaR[Q][ii];
deltaR[P][ii] += positions[particle3][ii];
deltaR[Q][ii] += positions[particle4][ii];
}
crossProductVector3( deltaR[CP], deltaR[DC], deltaR[T] );
crossProductVector3( deltaR[DC], deltaR[QD], deltaR[U] );
crossProductVector3( deltaR[T], deltaR[U], deltaR[TU] );
double rT2 = dotVector3( deltaR[T], deltaR[T] );
double rU2 = dotVector3( deltaR[U], deltaR[U] );
double rTrU = sqrt( rT2*rU2 );
if( rTrU <= 0.0 ){
return;
}
double rDC = dotVector3( deltaR[DC], deltaR[DC] );
rDC = sqrt( rDC );
double cosine = dotVector3( deltaR[T], deltaR[U] );
cosine /= rTrU;
double sine = dotVector3( deltaR[DC], deltaR[TU] );
sine /= ( rDC*rTrU );
double cosine2 = cosine*cosine - sine*sine;
double sine2 = 2.0*cosine*sine;
double phi2 = 1.0 - cosine2;
double dphi2 = 2.0*sine2;
double dedphi = kTorsion*dphi2;
for( int ii = 0; ii < 3; ii++ ){
deltaR[DP][ii] = positions[particle4][ii] - deltaR[P][ii];
deltaR[QC][ii] = deltaR[Q][ii] - positions[particle3][ii];
}
double factorT = dedphi/( rDC*rT2 );
double factorU = -dedphi/( rDC*rU2 );
crossProductVector3( deltaR[T], deltaR[DC], deltaR[dT] );
crossProductVector3( deltaR[U], deltaR[DC], deltaR[dU] );
for( int ii = 0; ii < 3; ii++ ){
deltaR[dT][ii] *= factorT;
deltaR[dU][ii] *= factorU;
}
crossProductVector3( deltaR[dT], deltaR[DC], deltaR[dP] );
crossProductVector3( deltaR[dU], deltaR[DC], deltaR[dQ] );
crossProductVector3( deltaR[DP], deltaR[dT], deltaR[dC1] );
crossProductVector3( deltaR[dU], deltaR[QD], deltaR[dC2] );
crossProductVector3( deltaR[dT], deltaR[CP], deltaR[dD1] );
crossProductVector3( deltaR[QC], deltaR[dU], deltaR[dD2] );
crossProductVector3( deltaR[BD], deltaR[dP], d[A] );
crossProductVector3( deltaR[dP], deltaR[AD], d[B] );
crossProductVector3( deltaR[FC], deltaR[dQ], d[E] );
crossProductVector3( deltaR[dQ], deltaR[EC], d[F] );
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[D][ii] = deltaR[dD1][ii] + deltaR[dD2][ii] + deltaR[dQ][ii] - d[A][ii] - d[B][ii];
}
// ---------------------------------------------------------------------------------------
// accumulate forces and energy
forces[particle1][0] -= d[0][0];
forces[particle1][1] -= d[0][1];
forces[particle1][2] -= d[0][2];
forces[particle2][0] -= d[1][0];
forces[particle2][1] -= d[1][1];
forces[particle2][2] -= d[1][2];
forces[particle3][0] -= d[2][0];
forces[particle3][1] -= d[2][1];
forces[particle3][2] -= d[2][2];
forces[particle4][0] -= d[3][0];
forces[particle4][1] -= d[3][1];
forces[particle4][2] -= d[3][2];
forces[particle5][0] -= d[4][0];
forces[particle5][1] -= d[4][1];
forces[particle5][2] -= d[4][2];
forces[particle6][0] -= d[5][0];
forces[particle6][1] -= d[5][1];
forces[particle6][2] -= d[5][2];
*energy += kTorsion*phi2;
return;
}
static void computeAmoebaPiTorsionForces( Context& context, AmoebaPiTorsionForce& amoebaPiTorsionForce,
std::vector<Vec3>& expectedForces, double* expectedEnergy, FILE* log ) {
// get positions and zero forces
State state = context.getState(State::Positions);
std::vector<Vec3> positions = state.getPositions();
expectedForces.resize( positions.size() );
for( unsigned int ii = 0; ii < expectedForces.size(); ii++ ){
expectedForces[ii][0] = expectedForces[ii][1] = expectedForces[ii][2] = 0.0;
}
// calculates forces/energy
*expectedEnergy = 0.0;
for( int ii = 0; ii < amoebaPiTorsionForce.getNumPiTorsions(); ii++ ){
computeAmoebaPiTorsionForce(ii, positions, amoebaPiTorsionForce, expectedForces, expectedEnergy, log );
}
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 );
}
return;
}
void compareWithExpectedForceAndEnergy( Context& context, AmoebaPiTorsionForce& amoebaPiTorsionForce,
double tolerance, const std::string& idString, FILE* log) {
std::vector<Vec3> expectedForces;
double expectedEnergy;
computeAmoebaPiTorsionForces( context, amoebaPiTorsionForce, expectedForces, &expectedEnergy, log );
State state = context.getState(State::Forces | State::Energy);
const std::vector<Vec3> forces = state.getForces();
if( log ){
(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 );
}
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 testOnePiTorsion( FILE* log ) {
System system;
int numberOfParticles = 6;
for( int ii = 0; ii < numberOfParticles; ii++ ){
system.addParticle(1.0);
}
LangevinIntegrator integrator(0.0, 0.1, 0.01);
AmoebaPiTorsionForce* amoebaPiTorsionForce = new AmoebaPiTorsionForce();
double kTorsion = 6.85;
amoebaPiTorsionForce->addPiTorsion(0, 1, 2, 3, 4, 5, kTorsion );
system.addForce(amoebaPiTorsionForce);
Context context(system, integrator, Platform::getPlatformByName( "Reference"));
std::vector<Vec3> positions(numberOfParticles);
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[2] = Vec3( 0.269130000E+02, 0.266390000E+02, 0.353100000E+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[5] = Vec3( 0.254124630E+02, 0.234691880E+02, 0.773335400E+01 );
context.setPositions(positions);
compareWithExpectedForceAndEnergy( context, *amoebaPiTorsionForce, TOL, "testOnePiTorsion", log );
}
int main( int numberOfArguments, char* argv[] ) {
try {
std::cout << "TestReferenceAmoebaPiTorsionForce running test..." << std::endl;
Platform::loadPluginsFromDirectory( Platform::getDefaultPluginsDirectory() );
//FILE* log = NULL;
FILE* log = stderr;
//FILE* log = fopen( "AmoebaPiTorsionForce1.log", "w" );;
testOnePiTorsion( log );
if( log && log != stderr )
(void) fclose( log );
}
catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl;
std::cout << "FAIL - ERROR. Test failed." << std::endl;
return 1;
}
std::cout << "PASS - Test succeeded." << std::endl;
return 0;
}
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaTorsionForce.cpp 0 → 100644
View file @ d91ae250
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the Reference implementation of ReferenceAmoebaTorsionForce.
*/
#include "../../../tests/AssertionUtilities.h"
//#include "AmoebaTinkerParameterFile.h"
#include "openmm/Context.h"
#include "AmoebaOpenMM.h"
#include "openmm/System.h"
#include "openmm/LangevinIntegrator.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
const double TOL = 1e-3;
#define PI_M 3.141592653589
#define RADIAN 57.29577951308
/* ---------------------------------------------------------------------------------------
Compute cross product of two 3-vectors and place in 3rd vector
vectorZ = vectorX x vectorY
@param vectorX x-vector
@param vectorY y-vector
@param vectorZ z-vector
@return vector is vectorZ
--------------------------------------------------------------------------------------- */
static void crossProductVector3( double* vectorX, double* vectorY, double* vectorZ ){
vectorZ[0] = vectorX[1]*vectorY[2] - vectorX[2]*vectorY[1];
vectorZ[1] = vectorX[2]*vectorY[0] - vectorX[0]*vectorY[2];
vectorZ[2] = vectorX[0]*vectorY[1] - vectorX[1]*vectorY[0];
return;
}
static double dotVector3( double* vectorX, double* vectorY ){
return vectorX[0]*vectorY[0] + vectorX[1]*vectorY[1] + vectorX[2]*vectorY[2];
}
static void computeAmoebaTorsionForce(int bondIndex, std::vector<Vec3>& positions, AmoebaTorsionForce& amoebaTorsionForce,
std::vector<Vec3>& forces, double* energy, FILE* log ) {
int particle1, particle2, particle3, particle4;
std::vector<double> torsion1;
std::vector<double> torsion2;
std::vector<double> torsion3;
torsion1.resize(3);
torsion2.resize(3);
torsion3.resize(3);
amoebaTorsionForce.getTorsionParameters(bondIndex, particle1, particle2, particle3, particle4, torsion1, torsion2, torsion3);
std::vector< std::vector<double> > torsions;
torsions.push_back( torsion1 );
torsions.push_back( torsion2 );
torsions.push_back( torsion3 );
if( log ){
(void) fprintf( log, "computeAmoebaTorsionForce: bond %d [%d %d %d %d]\n",
bondIndex, particle1, particle2, particle3, particle4 );
for( unsigned int ii = 0; ii < 3; ii++ ){
(void) fprintf( log, " [%10.3e %10.3e %10.3e]\n", torsions[ii][0], torsions[ii][1], torsions[ii][2] );
}
(void) fflush( log );
}
enum { BA, CB, DC, CA, DB, LastDeltaIndex };
double deltaR[LastDeltaIndex][3];
for( int ii = 0; ii < 3; ii++ ){
deltaR[BA][ii] = positions[particle2][ii] - positions[particle1][ii];
deltaR[CB][ii] = positions[particle3][ii] - positions[particle2][ii];
deltaR[DC][ii] = positions[particle4][ii] - positions[particle3][ii];
deltaR[CA][ii] = positions[particle3][ii] - positions[particle1][ii];
deltaR[DB][ii] = positions[particle4][ii] - positions[particle2][ii];
}
enum { Xt, Xu, Xtu, LastXtIndex };
double crossProducts[LastXtIndex][3];
crossProductVector3( deltaR[BA], deltaR[CB], crossProducts[Xt] );
crossProductVector3( deltaR[CB], deltaR[DC], crossProducts[Xu] );
crossProductVector3( crossProducts[Xt], crossProducts[Xu], crossProducts[Xtu] );
double rT2 = dotVector3( crossProducts[Xt], crossProducts[Xt] );
double rU2 = dotVector3( crossProducts[Xu], crossProducts[Xu] );
double rTrU = sqrt( rT2*rU2 );
if( rTrU <= 0.0 ){
return;
}
double rCB = dotVector3( deltaR[CB], deltaR[CB] );
rCB = sqrt( rCB );
// ---------------------------------------------------------------------------------------
// cos(w), cos(2w), cos(3w), ...
// sin(w), sin(2w), sin(3w), ...
double cosine[6], sine[6];
cosine[0] = dotVector3( crossProducts[Xt], crossProducts[Xu] );
cosine[0] /= rTrU;
sine[0] = dotVector3( deltaR[CB], crossProducts[Xtu] );
sine[0] /= (rCB*rTrU);
for( int ii = 1; ii < 3; ii++ ){
cosine[ii] = cosine[0]*cosine[ii-1] - sine[0]* sine[ii-1];
sine[ii] = cosine[0]* sine[ii-1] + sine[0]*cosine[ii-1];
}
// ---------------------------------------------------------------------------------------
// dEdPhi prefactor
double dEdPhi = 0.0;
for( int ii = 0; ii < 3; ii++ ){
dEdPhi += torsions[ii][0]*((double) (ii+1))*( cosine[ii]*sin( torsions[ii][1] ) - sine[ii]*cos( torsions[ii][1] ) );
}
fprintf( stderr,"Tst: dEdPhi=%15.7e\n", dEdPhi ); fflush( stderr );
// ---------------------------------------------------------------------------------------
// dEdtu[0] = dEdT
// dEdtu[1] = dEdU
// tempVector[0] == dEdA: dEdT x CB
// tempVector[1] == dEdB: (CA x dEdT) + (dEdU x DC)
// tempVector[2] == dEdC: (dEdT x BA) + (DB x dEdU)
// tempVector[3] == dEdD: (dEdU x CB)
double dEdtu[2][3];
double tempVector[6][3];
// dEdT & dEdU
crossProductVector3( crossProducts[Xt], deltaR[CB], tempVector[0] );
crossProductVector3( crossProducts[Xu], deltaR[CB], tempVector[1] );
double norm[2] = { dEdPhi/(rT2*rCB ), -dEdPhi/(rU2*rCB ) };
for( int jj = 0; jj < 2; jj++ ){
for( int ii = 0; ii < 3; ii++ ){
dEdtu[jj][ii] = norm[jj]*tempVector[jj][ii];
}
}
// dEdA
crossProductVector3( dEdtu[0], deltaR[CB], tempVector[0] );
// dEdB
crossProductVector3( deltaR[CA], dEdtu[0], tempVector[4] );
crossProductVector3( dEdtu[1], deltaR[DC], tempVector[1] );
// dEdC
crossProductVector3( dEdtu[0], deltaR[BA], tempVector[5] );
crossProductVector3( deltaR[DB], dEdtu[1], tempVector[2] );
for( int jj = 0; jj < 3; jj++ ){
tempVector[1][jj] += tempVector[4][jj];
tempVector[2][jj] += tempVector[5][jj];
}
// dEdD
crossProductVector3( dEdtu[1], deltaR[CB], tempVector[3] );
// ---------------------------------------------------------------------------------------
// accumulate forces and energy
forces[particle1][0] -= tempVector[0][0];
forces[particle1][1] -= tempVector[0][1];
forces[particle1][2] -= tempVector[0][2];
forces[particle2][0] -= tempVector[1][0];
forces[particle2][1] -= tempVector[1][1];
forces[particle2][2] -= tempVector[1][2];
forces[particle3][0] -= tempVector[2][0];
forces[particle3][1] -= tempVector[2][1];
forces[particle3][2] -= tempVector[2][2];
forces[particle4][0] -= tempVector[3][0];
forces[particle4][1] -= tempVector[3][1];
forces[particle4][2] -= tempVector[3][2];
double energyTerm = 0.0;
for( int ii = 0; ii < 3; ii++ ){
energyTerm += torsions[ii][0]*( 1.0 + cosine[ii]*cos( torsions[ii][1] ) + sine[ii]*sin( torsions[ii][1] ) );
}
*energy += energyTerm;
}
static void computeAmoebaTorsionForces( Context& context, AmoebaTorsionForce& amoebaTorsionForce,
std::vector<Vec3>& expectedForces, double* expectedEnergy, FILE* log ) {
// get positions and zero forces
State state = context.getState(State::Positions);
std::vector<Vec3> positions = state.getPositions();
expectedForces.resize( positions.size() );
for( unsigned int ii = 0; ii < expectedForces.size(); ii++ ){
expectedForces[ii][0] = expectedForces[ii][1] = expectedForces[ii][2] = 0.0;
}
// calculates forces/energy
*expectedEnergy = 0.0;
for( int ii = 0; ii < amoebaTorsionForce.getNumTorsions(); ii++ ){
computeAmoebaTorsionForce(ii, positions, amoebaTorsionForce, expectedForces, expectedEnergy, log );
}
if( log ){
(void) fprintf( log, "computeAmoebaTorsionForces: 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 );
}
return;
}
void compareWithExpectedForceAndEnergy( Context& context, AmoebaTorsionForce& amoebaTorsionForce,
double tolerance, const std::string& idString, FILE* log) {
std::vector<Vec3> expectedForces;
double expectedEnergy;
computeAmoebaTorsionForces( context, amoebaTorsionForce, expectedForces, &expectedEnergy, log );
State state = context.getState(State::Forces | State::Energy);
const std::vector<Vec3> forces = state.getForces();
if( log ){
(void) fprintf( log, "computeAmoebaTorsionForces: 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 );
}
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 testOneTorsion( FILE* log ) {
System system;
int numberOfParticles = 4;
for( int ii = 0; ii < numberOfParticles; ii++ ){
system.addParticle(1.0);
}
LangevinIntegrator integrator(0.0, 0.1, 0.01);
AmoebaTorsionForce* amoebaTorsionForce = new AmoebaTorsionForce();
std::vector<double> torsion1;
torsion1.push_back( 0.619500000E+00 );
torsion1.push_back( 0.000000000E+00 );
std::vector<double> torsion2;
torsion2.push_back( -0.202500000E+00 );
torsion2.push_back( 0.180000000E+03 );
std::vector<double> torsion3;
torsion3.push_back( 0.175000000E-01 );
torsion3.push_back( 0.000000000E+00 );
amoebaTorsionForce->addTorsion(0, 1, 2, 3, torsion1, torsion2, torsion3 );
system.addForce(amoebaTorsionForce);
Context context(system, integrator, Platform::getPlatformByName( "Reference"));
std::vector<Vec3> positions(numberOfParticles);
positions[0] = Vec3( 0.278860000E+01, 0.264630000E+01, 0.426300000E+00 );
positions[1] = Vec3( 0.273400000E+01, 0.244300000E+01, 0.261400000E+00 );
positions[2] = Vec3( 0.262660000E+01, 0.254130000E+01, 0.284200000E+00 );
positions[3] = Vec3( 0.269130000E+01, 0.266390000E+01, 0.353100000E+00 );
context.setPositions(positions);
compareWithExpectedForceAndEnergy( context, *amoebaTorsionForce, TOL, "testOneTorsion", log );
}
int main( int numberOfArguments, char* argv[] ) {
try {
std::cout << "TestReferenceAmoebaTorsionForce running test..." << std::endl;
Platform::loadPluginsFromDirectory( Platform::getDefaultPluginsDirectory() );
FILE* log = stderr;
//FILE* log = NULL;
//FILE* log = fopen( "AmoebaTorsionForce.log", "w" );;
testOneTorsion( log );
if( log && log != stderr )
(void) fclose( log );
}
catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl;
std::cout << "FAIL - ERROR. Test failed." << std::endl;
return 1;
}
std::cout << "PASS - Test succeeded." << std::endl;
return 0;
}
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