Skip to content

GitLab

Commit 8a57342c authored by Mark Friedrichs's avatar Mark Friedrichs
Browse files

Reference TorsionTorsion and Reference/Cuda tests for TorsionTorsion

parent 7f6c8bbc
Show whitespace changes
Inline Side-by-side
Showing with 3566 additions and 350 deletions
plugins/amoeba/openmmapi/include/AmoebaTorsionTorsionForce.h
View file @ 8a57342c
...@@ -34,6 +34,7 @@ ...@@ -34,6 +34,7 @@
#include "openmm/Force.h" #include "openmm/Force.h"
#include <vector> #include <vector>
#include <cmath>
#include "openmm/internal/windowsExport.h" #include "openmm/internal/windowsExport.h"
namespace OpenMM { namespace OpenMM {
...@@ -41,29 +42,30 @@ namespace OpenMM { ...@@ -41,29 +42,30 @@ namespace OpenMM {
typedef std::vector< std::vector< std::vector<double> > > TorsionTorsionGrid; typedef std::vector< std::vector< std::vector<double> > > TorsionTorsionGrid;
/** /**
* This class implements the Amoeba Out-of-plane bend interaction * This class implements the Amoeba torsion-torsion interaction
* To use it, create a TorsionTorsionForce object then call addTorsionTorsion() once for each torsionTorsion. After * To use it, create a TorsionTorsionForce object then call addTorsionTorsion() once for each torsionTorsion. After
* a torsionTorsion has been added, you can modify its force field parameters by calling setTorsionTorsionParameters(). * a torsionTorsion has been added, you can modify its force field parameters by calling setTorsionTorsionParameters().
*/ */
class OPENMM_EXPORT AmoebaTorsionTorsionForce : public Force { class OPENMM_EXPORT AmoebaTorsionTorsionForce : public Force {
public: public:
/** /**
* Create a Amoeba TorsionTorsionForce. * Create a Amoeba TorsionTorsionForce.
*/ */
AmoebaTorsionTorsionForce(); AmoebaTorsionTorsionForce( void );
/** /**
* Get the number of torsionTorsion terms in the potential function * Get the number of torsionTorsion terms in the potential function
*/ */
int getNumTorsionTorsions() const { int getNumTorsionTorsions( void ) const {
return torsionTorsions.size(); return torsionTorsions.size();
} }
/** /**
* Get the number of torsionTorsion grids * Get the number of torsionTorsion grids
*/ */
int getNumTorsionTorsionGrids() const { int getNumTorsionTorsionGrids( void ) const {
return torsionTorsionGrids.size(); return torsionTorsionGrids.size();
} }
...@@ -76,7 +78,7 @@ public: ...@@ -76,7 +78,7 @@ public:
* @param particle4 the index of the fourth particle connected by the torsionTorsion * @param particle4 the index of the fourth particle connected by the torsionTorsion
* @param particle5 the index of the fifth particle connected by the torsionTorsion * @param particle5 the index of the fifth particle connected by the torsionTorsion
* @param chiralCheckAtomIndex the index of the particle connected to particle3, but not particle2 or particle4 to be used in chirality check * @param chiralCheckAtomIndex the index of the particle connected to particle3, but not particle2 or particle4 to be used in chirality check
* @param gridIndex the grid index * @param gridIndex the index to the grid to be used
* @return the index of the torsionTorsion that was added * @return the index of the torsionTorsion that was added
*/ */
int addTorsionTorsion(int particle1, int particle2, int particle3, int particle4, int particle5, int chiralCheckAtomIndex, int gridIndex ); int addTorsionTorsion(int particle1, int particle2, int particle3, int particle4, int particle5, int chiralCheckAtomIndex, int gridIndex );
...@@ -93,7 +95,7 @@ public: ...@@ -93,7 +95,7 @@ public:
* @param chiralCheckAtomIndex the index of the particle connected to particle3, but not particle2 or particle4 to be used in chirality check * @param chiralCheckAtomIndex the index of the particle connected to particle3, but not particle2 or particle4 to be used in chirality check
* @param gridIndex the grid index * @param gridIndex the grid index
*/ */
void getTorsionTorsionParameters(int index, int& particle1, int& particle2, int& particle3, int& particle4, int& particle5, int& chiralCheckAtomIndex, int& gridIndex ) const; void getTorsionTorsionParameters( int index, int& particle1, int& particle2, int& particle3, int& particle4, int& particle5, int& chiralCheckAtomIndex, int& gridIndex ) const;
/** /**
* Set the force field parameters for a torsionTorsion term. * Set the force field parameters for a torsionTorsion term.
...@@ -107,15 +109,15 @@ public: ...@@ -107,15 +109,15 @@ public:
* @param chiralCheckAtomIndex the index of the particle connected to particle3, but not particle2 or particle4 to be used in chirality check * @param chiralCheckAtomIndex the index of the particle connected to particle3, but not particle2 or particle4 to be used in chirality check
* @param gridIndex the grid index * @param gridIndex the grid index
*/ */
void setTorsionTorsionParameters(int index, int particle1, int particle2, int particle3, int particle4, int particle5, int chiralCheckAtomIndex, int gridIndex ); void setTorsionTorsionParameters( int index, int particle1, int particle2, int particle3, int particle4, int particle5, int chiralCheckAtomIndex, int gridIndex );
/** /**
* Get the torsion-torsion grid at the specified index * Get the torsion-torsion grid at the specified index
* *
* @param gridIndex the grid index * @param gridIndex the grid index
* @param grid the grid * @return grid return grid reference
*/ */
void getTorsionTorsionGrid(int index, TorsionTorsionGrid& grid ) const; const TorsionTorsionGrid& getTorsionTorsionGrid( int index ) const;
/** /**
* Set the torsion-torsion grid at the specified index * Set the torsion-torsion grid at the specified index
...@@ -136,6 +138,7 @@ protected: ...@@ -136,6 +138,7 @@ protected:
private: private:
class TorsionTorsionInfo; class TorsionTorsionInfo;
class TorsionTorsionGridInfo;
// Retarded visual studio compiler complains about being unable to // Retarded visual studio compiler complains about being unable to
// export private stl class members. // export private stl class members.
...@@ -145,7 +148,7 @@ private: ...@@ -145,7 +148,7 @@ private:
#pragma warning(disable:4251) #pragma warning(disable:4251)
#endif #endif
std::vector<TorsionTorsionInfo> torsionTorsions; std::vector<TorsionTorsionInfo> torsionTorsions;
std::vector<TorsionTorsionGrid> torsionTorsionGrids; std::vector<TorsionTorsionGridInfo> torsionTorsionGrids;
#if defined(_MSC_VER) #if defined(_MSC_VER)
#pragma warning(pop) #pragma warning(pop)
#endif #endif
...@@ -168,6 +171,63 @@ public: ...@@ -168,6 +171,63 @@ public:
} }
}; };
class AmoebaTorsionTorsionForce::TorsionTorsionGridInfo {
public:
TorsionTorsionGridInfo( ) {
_size[0] = _size[1] = 0;
_startValues[0] = _startValues[1] = 0.0;
_spacing[0] = _spacing[1] = 1.0;
}
TorsionTorsionGridInfo( const TorsionTorsionGrid& grid ) {
_grid.resize( grid.size() );
for( unsigned int kk = 0; kk < grid.size(); kk++ ){
_grid[kk].resize( grid[kk].size() );
for( unsigned int jj = 0; jj < grid[kk].size(); jj++ ){
_grid[kk][jj].resize( grid[kk][jj].size() );
for( unsigned int ii = 0; ii < grid[kk][jj].size(); ii++ ){
_grid[kk][jj][ii] = grid[kk][jj][ii];
}
}
}
_startValues[0] = _grid[0][0][0];
_startValues[1] = _grid[0][0][1];
//_spacing[0] = fabs( _grid[1][0][0] - _grid[0][0][0] );
//_spacing[1] = fabs( _grid[0][1][1] - _grid[0][0][1] );
_spacing[0] = static_cast<double>(_grid.size()-1)/360.0;
_spacing[1] = static_cast<double>(grid.size()-1)/360.0;
_size[0] = static_cast<int>(grid.size());
_size[1] = static_cast<int>(grid[0].size());
}
const TorsionTorsionGrid& getTorsionTorsionGrid( void ) const {
return _grid;
}
int getDimensionSize( int index ) const {
return _size[index];
}
double getStartValue( int index ) const {
return _startValues[index];
}
double getSpacing( int index ) const {
return _spacing[index];
}
private:
TorsionTorsionGrid _grid;
int _size[2];
double _startValues[2];
double _spacing[2];
};
} // namespace OpenMM } // namespace OpenMM
#endif /*OPENMM_AMOEBA_TORSION_TORSION_FORCE_H_*/ #endif /*OPENMM_AMOEBA_TORSION_TORSION_FORCE_H_*/
plugins/amoeba/openmmapi/src/AmoebaTorsionTorsionForce.cpp
View file @ 8a57342c
...@@ -68,8 +68,8 @@ void AmoebaTorsionTorsionForce::setTorsionTorsionParameters(int index, int parti ...@@ -68,8 +68,8 @@ void AmoebaTorsionTorsionForce::setTorsionTorsionParameters(int index, int parti
torsionTorsions[index].gridIndex = gridIndex; torsionTorsions[index].gridIndex = gridIndex;
} }
void AmoebaTorsionTorsionForce::getTorsionTorsionGrid(int index, TorsionTorsionGrid& grid ) const { const TorsionTorsionGrid& AmoebaTorsionTorsionForce::getTorsionTorsionGrid(int index ) const {
grid = torsionTorsionGrids[index]; return torsionTorsionGrids[index].getTorsionTorsionGrid();
} }
void AmoebaTorsionTorsionForce::setTorsionTorsionGrid(int index, TorsionTorsionGrid& grid ) { void AmoebaTorsionTorsionForce::setTorsionTorsionGrid(int index, TorsionTorsionGrid& grid ) {
......
plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.cpp
View file @ 8a57342c
...@@ -398,20 +398,19 @@ void CudaCalcAmoebaTorsionTorsionForceKernel::initialize(const System& system, c ...@@ -398,20 +398,19 @@ void CudaCalcAmoebaTorsionTorsionForceKernel::initialize(const System& system, c
std::vector< std::vector< std::vector< std::vector<float> > > > floatGrids; std::vector< std::vector< std::vector< std::vector<float> > > > floatGrids;
floatGrids.resize(numTorsionTorsionGrids); floatGrids.resize(numTorsionTorsionGrids);
for (int i = 0; i < numTorsionTorsionGrids; i++) { for (int gridIndex = 0; gridIndex < numTorsionTorsionGrids; gridIndex++) {
TorsionTorsionGrid grid; const TorsionTorsionGrid& grid = force.getTorsionTorsionGrid( gridIndex );
force.getTorsionTorsionGrid(i, grid );
floatGrids[i].resize( grid.size() ); floatGrids[gridIndex].resize( grid.size() );
for (unsigned int ii = 0; ii < grid.size(); ii++) { for (unsigned int ii = 0; ii < grid.size(); ii++) {
floatGrids[i][ii].resize( grid[ii].size() ); floatGrids[gridIndex][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++) {
floatGrids[i][ii][jj].resize( grid[ii][jj].size() ); floatGrids[gridIndex][ii][jj].resize( grid[ii][jj].size() );
for (unsigned int kk = 0; kk < grid[ii][kk].size(); kk++) { for (unsigned int kk = 0; kk < grid[ii][kk].size(); kk++) {
floatGrids[i][ii][jj][kk] = static_cast<float>(grid[ii][jj][kk]); floatGrids[gridIndex][ii][jj][kk] = static_cast<float>(grid[ii][jj][kk]);
} }
} }
} }
......
plugins/amoeba/platforms/cuda/src/kernels/AmoebaGpu.cpp
View file @ 8a57342c
...@@ -950,9 +950,9 @@ void gpuSetAmoebaTorsionTorsionGrids(amoebaGpuContext amoebaGpu, const std::vect ...@@ -950,9 +950,9 @@ void gpuSetAmoebaTorsionTorsionGrids(amoebaGpuContext amoebaGpu, const std::vect
// 4 (grids) * (25 *25 grid)*(2 +4 a1, a2, f, f1,f2, f12) = 15000 // 4 (grids) * (25 *25 grid)*(2 +4 a1, a2, f, f1,f2, f12) = 15000
for (unsigned int ii = 0; ii < floatGrids.size(); ii++) { for (unsigned int ii = 0; ii < floatGrids.size(); ii++) {
amoebaGpu->amoebaSim.amoebaTorTorGridOffset[ii] = (totalGridEntries/4); amoebaGpu->amoebaSim.amoebaTorTorGridOffset[ii] = (totalGridEntries/4);
amoebaGpu->amoebaSim.amoebaTorTorGridBegin[ii] = -180.0f; amoebaGpu->amoebaSim.amoebaTorTorGridBegin[ii] = floatGrids[ii][0][0][0];
amoebaGpu->amoebaSim.amoebaTorTorGridDelta[ii] = 15.0f; amoebaGpu->amoebaSim.amoebaTorTorGridDelta[ii] = 360.0f/static_cast<float>(floatGrids[ii].size()-1);
amoebaGpu->amoebaSim.amoebaTorTorGridNy[ii] = 25; amoebaGpu->amoebaSim.amoebaTorTorGridNy[ii] = floatGrids[ii].size();
for (unsigned int jj = 0; jj < floatGrids[ii].size(); jj++) { for (unsigned int jj = 0; jj < floatGrids[ii].size(); jj++) {
for (unsigned int kk = 0; kk < floatGrids[ii][jj].size(); kk++) { for (unsigned int kk = 0; kk < floatGrids[ii][jj].size(); kk++) {
totalGridEntries += (floatGrids[ii][jj][kk].size() - 2); totalGridEntries += (floatGrids[ii][jj][kk].size() - 2);
...@@ -966,7 +966,11 @@ void gpuSetAmoebaTorsionTorsionGrids(amoebaGpuContext amoebaGpu, const std::vect ...@@ -966,7 +966,11 @@ void gpuSetAmoebaTorsionTorsionGrids(amoebaGpuContext amoebaGpu, const std::vect
amoebaGpu->amoebaSim.pAmoebaTorsionTorsionGrids = psTorsionTorsionGrids->_pDevStream[0]; amoebaGpu->amoebaSim.pAmoebaTorsionTorsionGrids = psTorsionTorsionGrids->_pDevStream[0];
if( amoebaGpu->log ){ if( amoebaGpu->log ){
(void) fprintf( amoebaGpu->log, "totalGridEntries=%u totalFloat4 entries=%u\n", totalGridEntries, totalEntries ); (void) fprintf( amoebaGpu->log, "Grids %u totalGridEntries=%u totalFloat4 entries=%u\n", torsionTorsionGrids, totalGridEntries, totalEntries );
for (unsigned int ii = 0; ii < floatGrids.size(); ii++) {
(void) fprintf( amoebaGpu->log, "Grid %u offset=%d begin=%.3f delta=%.3f Ny=%d\n", ii, amoebaGpu->amoebaSim.amoebaTorTorGridOffset[ii],
amoebaGpu->amoebaSim.amoebaTorTorGridBegin[ii], amoebaGpu->amoebaSim.amoebaTorTorGridDelta[ii], amoebaGpu->amoebaSim.amoebaTorTorGridNy[ii]);
}
} }
unsigned int index = 0; unsigned int index = 0;
...@@ -982,7 +986,7 @@ void gpuSetAmoebaTorsionTorsionGrids(amoebaGpuContext amoebaGpu, const std::vect ...@@ -982,7 +986,7 @@ void gpuSetAmoebaTorsionTorsionGrids(amoebaGpuContext amoebaGpu, const std::vect
#define DUMP_PARAMETERS 5 #define DUMP_PARAMETERS 5
#if (DUMP_PARAMETERS > 0 ) #if (DUMP_PARAMETERS > 0 )
if( (index < DUMP_PARAMETERS || index > totalEntries - (DUMP_PARAMETERS + 1)) && amoebaGpu->log ) if( (index < DUMP_PARAMETERS || index > totalEntries - (DUMP_PARAMETERS + 1)) && amoebaGpu->log )
fprintf( amoebaGpu->log, "TorsionTorsionGrid: %d %5d [%5d %5d ] [%10.3f %10.3f] [%15.7e %15.7e %15.7e %15.7e]\n", index, ii, jj, kk, (void) fprintf( amoebaGpu->log, "TorsionTorsionGrid: %5d %5d [%5d %5d ] [%10.3f %10.3f] [%15.7e %15.7e %15.7e %15.7e]\n", index, ii, jj, kk,
floatGrids[ii][jj][kk][0], floatGrids[ii][jj][kk][1], floatGrids[ii][jj][kk][0], floatGrids[ii][jj][kk][1],
(*psTorsionTorsionGrids)[index].x, (*psTorsionTorsionGrids)[index].y, (*psTorsionTorsionGrids)[index].z, (*psTorsionTorsionGrids)[index].w ); (*psTorsionTorsionGrids)[index].x, (*psTorsionTorsionGrids)[index].y, (*psTorsionTorsionGrids)[index].z, (*psTorsionTorsionGrids)[index].w );
#endif #endif
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaTorsionTorsionForce.cpp
View file @ 8a57342c
This source diff could not be displayed because it is too large. You can view the blob instead.
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernelFactory.cpp
View file @ 8a57342c
...@@ -46,8 +46,8 @@ fprintf( stderr,"In registerKernelFactories AmoebaReferenceKernelFactory\n" ); f ...@@ -46,8 +46,8 @@ fprintf( stderr,"In registerKernelFactories AmoebaReferenceKernelFactory\n" ); f
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);
/*
platform.registerKernelFactory(CalcAmoebaMultipoleForceKernel::Name(), factory); platform.registerKernelFactory(CalcAmoebaMultipoleForceKernel::Name(), factory);
platform.registerKernelFactory(CalcAmoebaGeneralizedKirkwoodForceKernel::Name(), factory); platform.registerKernelFactory(CalcAmoebaGeneralizedKirkwoodForceKernel::Name(), factory);
platform.registerKernelFactory(CalcAmoebaVdwForceKernel::Name(), factory); platform.registerKernelFactory(CalcAmoebaVdwForceKernel::Name(), factory);
...@@ -83,9 +83,9 @@ KernelImpl* AmoebaReferenceKernelFactory::createKernelImpl(std::string name, con ...@@ -83,9 +83,9 @@ KernelImpl* AmoebaReferenceKernelFactory::createKernelImpl(std::string name, con
if (name == CalcAmoebaOutOfPlaneBendForceKernel::Name()) if (name == CalcAmoebaOutOfPlaneBendForceKernel::Name())
return new ReferenceCalcAmoebaOutOfPlaneBendForceKernel(name, platform, context.getSystem()); return new ReferenceCalcAmoebaOutOfPlaneBendForceKernel(name, platform, context.getSystem());
/*
if (name == CalcAmoebaTorsionTorsionForceKernel::Name()) if (name == CalcAmoebaTorsionTorsionForceKernel::Name())
return new ReferenceCalcAmoebaTorsionTorsionForceKernel(name, platform, context.getSystem()); return new ReferenceCalcAmoebaTorsionTorsionForceKernel(name, platform, context.getSystem());
/*
if (name == CalcAmoebaMultipoleForceKernel::Name()) if (name == CalcAmoebaMultipoleForceKernel::Name())
return new ReferenceCalcAmoebaMultipoleForceKernel(name, platform, context.getSystem()); return new ReferenceCalcAmoebaMultipoleForceKernel(name, platform, context.getSystem());
......
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.cpp
View file @ 8a57342c
...@@ -32,6 +32,7 @@ ...@@ -32,6 +32,7 @@
#include "AmoebaReferencePiTorsionForce.h" #include "AmoebaReferencePiTorsionForce.h"
#include "AmoebaReferenceStretchBendForce.h" #include "AmoebaReferenceStretchBendForce.h"
#include "AmoebaReferenceOutOfPlaneBendForce.h" #include "AmoebaReferenceOutOfPlaneBendForce.h"
#include "AmoebaReferenceTorsionTorsionForce.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"
...@@ -427,72 +428,95 @@ double ReferenceCalcAmoebaOutOfPlaneBendForceKernel::execute(ContextImpl& contex ...@@ -427,72 +428,95 @@ double ReferenceCalcAmoebaOutOfPlaneBendForceKernel::execute(ContextImpl& contex
return static_cast<double>(energy); return static_cast<double>(energy);
} }
//ReferenceCalcAmoebaTorsionTorsionForceKernel::ReferenceCalcAmoebaTorsionTorsionForceKernel(std::string name, const Platform& platform, System& system) : ReferenceCalcAmoebaTorsionTorsionForceKernel::ReferenceCalcAmoebaTorsionTorsionForceKernel(std::string name, const Platform& platform, System& system) :
// CalcAmoebaTorsionTorsionForceKernel(name, platform), system(system) { CalcAmoebaTorsionTorsionForceKernel(name, platform), system(system) {
// data.incrementKernelCount(); }
//}
// ReferenceCalcAmoebaTorsionTorsionForceKernel::~ReferenceCalcAmoebaTorsionTorsionForceKernel() {
//ReferenceCalcAmoebaTorsionTorsionForceKernel::~ReferenceCalcAmoebaTorsionTorsionForceKernel() { }
// data.decrementKernelCount();
//} void ReferenceCalcAmoebaTorsionTorsionForceKernel::initialize(const System& system, const AmoebaTorsionTorsionForce& force) {
//
//void ReferenceCalcAmoebaTorsionTorsionForceKernel::initialize(const System& system, const AmoebaTorsionTorsionForce& force) { numTorsionTorsions = force.getNumTorsionTorsions();
//
// data.setAmoebaLocalForcesKernel( this ); // torsion-torsion parameters
// numTorsionTorsions = force.getNumTorsionTorsions();
// for (int ii = 0; ii < numTorsionTorsions; ii++) {
// // torsion-torsion parameters int particle1Index, particle2Index, particle3Index, particle4Index, particle5Index, chiralCheckAtomIndex, gridIndex;
// force.getTorsionTorsionParameters(ii, particle1Index, particle2Index, particle3Index,
// std::vector<int> particle1(numTorsionTorsions); particle4Index, particle5Index, chiralCheckAtomIndex, gridIndex);
// std::vector<int> particle2(numTorsionTorsions); particle1.push_back( particle1Index );
// std::vector<int> particle3(numTorsionTorsions); particle2.push_back( particle2Index );
// std::vector<int> particle4(numTorsionTorsions); particle3.push_back( particle3Index );
// std::vector<int> particle5(numTorsionTorsions); particle4.push_back( particle4Index );
// std::vector<int> chiralCheckAtomIndex(numTorsionTorsions); particle5.push_back( particle5Index );
// std::vector<int> gridIndices(numTorsionTorsions); chiralCheckAtom.push_back( chiralCheckAtomIndex );
// gridIndices.push_back( gridIndex );
// for (int i = 0; i < numTorsionTorsions; i++) { }
// force.getTorsionTorsionParameters(i, particle1[i], particle2[i], particle3[i],
// particle4[i], particle5[i], // torsion-torsion grids
// chiralCheckAtomIndex[i], gridIndices[i]);
// } numTorsionTorsionGrids = force.getNumTorsionTorsionGrids();
// gpuSetAmoebaTorsionTorsionParameters(data.getAmoebaGpu(), particle1, particle2, particle3, particle4, particle5, chiralCheckAtomIndex, gridIndices ); torsionTorsionGrids.resize(numTorsionTorsionGrids);
// for (int ii = 0; ii < numTorsionTorsionGrids; ii++) {
// // torsion-torsion grids
// const TorsionTorsionGrid grid = force.getTorsionTorsionGrid( ii );
// numTorsionTorsionGrids = force.getNumTorsionTorsionGrids();
// std::vector< std::vector< std::vector< std::vector<RealOpenMM> > > > RealOpenMMGrids; torsionTorsionGrids[ii].resize( grid.size() );
// for (unsigned int kk = 0; kk < grid.size(); kk++) {
// RealOpenMMGrids.resize(numTorsionTorsionGrids);
// for (int i = 0; i < numTorsionTorsionGrids; i++) { torsionTorsionGrids[ii][kk].resize( grid[kk].size() );
// for (unsigned int jj = 0; jj < grid[kk].size(); jj++) {
// TorsionTorsionGrid grid;
// force.getTorsionTorsionGrid(i, grid ); torsionTorsionGrids[ii][kk][jj].resize( grid[kk][jj].size() );
// for (unsigned int ll = 0; ll < grid[ll][jj].size(); ll++) {
// RealOpenMMGrids[i].resize( grid.size() ); torsionTorsionGrids[ii][kk][jj][ll] = static_cast<RealOpenMM>(grid[kk][jj][ll]);
// for (unsigned int ii = 0; ii < grid.size(); ii++) { }
// }
// RealOpenMMGrids[i][ii].resize( grid[ii].size() ); }
// for (unsigned int jj = 0; jj < grid[ii].size(); jj++) { }
// }
// RealOpenMMGrids[i][ii][jj].resize( grid[ii][jj].size() );
// for (unsigned int kk = 0; kk < grid[ii][kk].size(); kk++) { double ReferenceCalcAmoebaTorsionTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
// RealOpenMMGrids[i][ii][jj][kk] = static_cast<RealOpenMM>(grid[ii][jj][kk]);
// } RealOpenMM** posData = extractPositions(context);
// } RealOpenMM** forceData = extractForces(context);
// } RealOpenMM energy = 0.0;
// }
// gpuSetAmoebaTorsionTorsionGrids(data.getAmoebaGpu(), RealOpenMMGrids ); for( unsigned int ii = 0; ii < numTorsionTorsions; ii++ ){
//
//} int particle1Index = particle1[ii];
// int particle2Index = particle2[ii];
//double CudaCalcAmoebaTorsionTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) { int particle3Index = particle3[ii];
// if( data.getAmoebaLocalForcesKernel() == this ){ int particle4Index = particle4[ii];
// computeAmoebaLocalForces( data ); int particle5Index = particle5[ii];
// }
// return 0.0; int chiralCheckAtomIndex = chiralCheckAtom[ii];
//}
// int gridIndex = gridIndices[ii];
RealOpenMM* forces[5];
forces[0] = forceData[particle1Index];
forces[1] = forceData[particle2Index];
forces[2] = forceData[particle3Index];
forces[3] = forceData[particle4Index];
forces[4] = forceData[particle5Index];
RealOpenMM* chiralCheckAtom;
if( chiralCheckAtomIndex >= 0 ){
chiralCheckAtom = posData[chiralCheckAtomIndex];
} else {
chiralCheckAtom = NULL;
}
energy += AmoebaReferenceTorsionTorsionForce::calculateForceAndEnergy(
posData[particle1Index], posData[particle2Index],
posData[particle3Index], posData[particle4Index],
posData[particle5Index], chiralCheckAtom, torsionTorsionGrids[gridIndex],
forces );
}
return static_cast<double>(energy);
}
///* -------------------------------------------------------------------------- * ///* -------------------------------------------------------------------------- *
// * AmoebaMultipole * // * AmoebaMultipole *
// * -------------------------------------------------------------------------- */ // * -------------------------------------------------------------------------- */
......
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.h
View file @ 8a57342c
...@@ -286,35 +286,45 @@ private: ...@@ -286,35 +286,45 @@ private:
System& system; System& system;
}; };
// /** /**
// * This kernel is invoked by AmoebaTorsionTorsionForce to calculate the forces acting on the system and the energy of the system. * This kernel is invoked by AmoebaTorsionTorsionForce to calculate the forces acting on the system and the energy of the system.
// */ */
// class ReferenceCalcAmoebaTorsionTorsionForceKernel : public CalcAmoebaTorsionTorsionForceKernel { class ReferenceCalcAmoebaTorsionTorsionForceKernel : public CalcAmoebaTorsionTorsionForceKernel {
// public: public:
// ReferenceCalcAmoebaTorsionTorsionForceKernel(std::string name, const Platform& platform, System& system); ReferenceCalcAmoebaTorsionTorsionForceKernel(std::string name, const Platform& platform, System& system);
// ~ReferenceCalcAmoebaTorsionTorsionForceKernel(); ~ReferenceCalcAmoebaTorsionTorsionForceKernel();
// /** /**
// * 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 AmoebaTorsionTorsionForce this kernel will be used for * @param force the AmoebaTorsionTorsionForce this kernel will be used for
// */ */
// void initialize(const System& system, const AmoebaTorsionTorsionForce& force); void initialize(const System& system, const AmoebaTorsionTorsionForce& 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 numTorsionTorsions; int numTorsionTorsions;
// int numTorsionTorsionGrids; std::vector<int> particle1;
// System& system; std::vector<int> particle2;
// }; std::vector<int> particle3;
// std::vector<int> particle4;
std::vector<int> particle5;
std::vector<int> chiralCheckAtom;
std::vector<int> gridIndices;
int numTorsionTorsionGrids;
std::vector< std::vector< std::vector< std::vector<RealOpenMM> > > > torsionTorsionGrids;
System& system;
};
// /** // /**
// * This kernel is invoked by AmoebaMultipoleForce to calculate the forces acting on the system and the energy of the system. // * This kernel is invoked by AmoebaMultipoleForce to calculate the forces acting on the system and the energy of the system.
// */ // */
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceTorsionTorsionForce.cpp 0 → 100644
View file @ 8a57342c
/* 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 "AmoebaReferenceTorsionTorsionForce.h"
#include <vector>
/**---------------------------------------------------------------------------------------
Load grid values from rectangle enclosing angles
@param grid grid values [ angle1, angle2, f, df1, df2, df12 ]
@param angle1 angle in first dimension
@param angle2 angle in second dimension
@param corners on return contains the coordinates of the rectangle
@param fValues on return contains the values of the function at the vertices of the rectangle
@param fValues1 on return contains the values of the derivative of the function wrt first dimension
@param fValues2 on return contains the values of the derivative of the function wrt second dimension
@param fValues12 on return contains the values of the derivative of the function wrt first & second dimension
@return AmoebaCommon::DefaultReturn
On first call a check is performed to see if the grid is valid
--------------------------------------------------------------------------------------- */
void AmoebaReferenceTorsionTorsionForce::loadGridValuesFromEnclosingRectangle(
const std::vector< std::vector< std::vector<RealOpenMM> > >& grid,
RealOpenMM angle1, RealOpenMM angle2, RealOpenMM corners[2][2],
RealOpenMM* fValues, RealOpenMM* fValues1, RealOpenMM* fValues2, RealOpenMM* fValues12 ){
// ---------------------------------------------------------------------------------------
// static const std::string methodName = "loadGridValuesFromEnclosingRectangle";
// ---------------------------------------------------------------------------------------
// get 2 opposing grid indices for rectangle
unsigned int gridSize = grid.size();
double gridSpacingI = static_cast<double>(gridSize-1)/360.0;
const int X_gridIndex = (int) ( (angle1 - grid[0][0][0] )*gridSpacingI + 1.0e-06);
const int Y_gridIndex = (int) ( (angle2 - grid[0][0][1] )*gridSpacingI + 1.0e-06);
// get coordinates of corner indices
corners[0][0] = grid[X_gridIndex][Y_gridIndex ][0];
corners[0][1] = grid[X_gridIndex+1][Y_gridIndex][0];
corners[1][0] = grid[X_gridIndex][Y_gridIndex ][1];
corners[1][1] = grid[X_gridIndex+1][Y_gridIndex+1][1];
for( int ii = 0; ii < 4; ii++ ){
int gridX = X_gridIndex;
int gridY = Y_gridIndex;
if( ii == 1 ){
gridX++;
} else if( ii == 2 ){
gridX++;
gridY++;
} else if( ii == 3 ){
gridY++;
}
fValues[ii] = grid[gridX][gridY][2];
fValues1[ii] = grid[gridX][gridY][3];
fValues2[ii] = grid[gridX][gridY][4];
fValues12[ii] = grid[gridX][gridY][5];
}
return;
}
/**---------------------------------------------------------------------------------------
Determines the coefficient matrix needed for bicubic
interpolation of a function, gradients and cross derivatives
Reference:
W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P.
Flannery, Numerical Recipes (Fortran), 2nd Ed., Cambridge
University Press, 1992, Section 3.6
@param y y values
@param y1 dy/dx1 values
@param y2 dy/dx2 values
@param y12 d2y/dx1dx2 values
@param d1 d1Upper - d1Lower
@param d2 d2Upper - d2Lower
@param c 4x4 return coefficient matrix
--------------------------------------------------------------------------------------- */
void AmoebaReferenceTorsionTorsionForce::getBicubicCoefficientMatrix( const RealOpenMM* y,
const RealOpenMM* y1, const RealOpenMM* y2, const RealOpenMM* y12, const RealOpenMM d1, const RealOpenMM d2,
RealOpenMM c[4][4] ){
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "AmoebaReferenceTorsionTorsionForce::getBicubicCoefficientMatrix";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM two = 2.0;
static const RealOpenMM three = 3.0;
static const RealOpenMM four = 4.0;
static const RealOpenMM five = 5.0;
static const RealOpenMM six = 6.0;
static const RealOpenMM nine = 9.0;
// transpose of matrix in Tinker due to difference in C/Fotran row/column major
// change indices when multiplying by weightMatrix
static const RealOpenMM weightMatrix[16][16] = {
{ one, zero, -three, two, zero, zero, zero, zero, -three, zero, nine, -six, two, zero, -six, four },
{ zero, zero, zero, zero, zero, zero, zero, zero, three, zero, -nine, six, -two, zero, six, -four },
{ zero, zero, zero, zero, zero, zero, zero, zero, zero, zero, nine, -six, zero, zero, -six, four },
{ zero, zero, three, -two, zero, zero, zero, zero, zero, zero, -nine, six, zero, zero, six, -four },
{ zero, zero, zero, zero, one, zero, -three, two, -two, zero, six, -four, one, zero, -three, two },
{ zero, zero, zero, zero, zero, zero, zero, zero, -one, zero, three, -two, one, zero, -three, two },
{ zero, zero, zero, zero, zero, zero, zero, zero, zero, zero, -three, two, zero, zero, three, -two },
{ zero, zero, zero, zero, zero, zero, three, -two, zero, zero, -six, four, zero, zero, three, -two },
{ zero, one, -two, one, zero, zero, zero, zero, zero, -three, six, -three, zero, two, -four, two },
{ zero, zero, zero, zero, zero, zero, zero, zero, zero, three, -six, three, zero, -two, four, -two },
{ zero, zero, zero, zero, zero, zero, zero, zero, zero, zero, -three, three, zero, zero, two, -two },
{ zero, zero, -one, one, zero, zero, zero, zero, zero, zero, three, -three, zero, zero, -two, two },
{ zero, zero, zero, zero, zero, one, -two, one, zero, -two, four, -two, zero, one, -two, one },
{ zero, zero, zero, zero, zero, zero, zero, zero, zero, -one, two, -one, zero, one, -two, one },
{ zero, zero, zero, zero, zero, zero, zero, zero, zero, zero, one, -one, zero, zero, -one, one },
{ zero, zero, zero, zero, zero, zero, -one, one, zero, zero, two, -two, zero, zero, -one, one } };
// ---------------------------------------------------------------------------------------
// pack y, y1, y2, y12 into single vector of dimension 16
std::vector<RealOpenMM> x( 16 );
RealOpenMM d1d2 = d1*d2;
for( int ii = 0; ii < 4; ii++ ){
x[ii] = y[ii];
x[ii+4] = y1[ii]*d1;
x[ii+8] = y2[ii]*d2;
x[ii+12] = y12[ii]*d1d2;
}
// matrix multiply by the transpose of the stored weight table
int rowIndex = 0;
int colIndex = 0;
for( int ii = 0; ii < 16; ii++ ){
RealOpenMM sum = weightMatrix[0][ii]*x[0];
for( int jj = 1; jj < 16; jj++ ){
sum += weightMatrix[jj][ii]*x[jj];
}
c[rowIndex][colIndex++] = sum;
if( !(colIndex % 4) ){
colIndex = 0;
rowIndex++;
}
}
}
/**---------------------------------------------------------------------------------------
Determines the coefficient matrix needed for bicubic
interpolation of a function, gradients and cross derivatives
Reference:
W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P.
Flannery, Numerical Recipes (Fortran), 2nd Ed., Cambridge
University Press, 1992, Section 3.6
@param y y values
@param y1 dy/dx1 values
@param y2 dy/dx2 values
@param y12 d2y/dx1dx2 values
@param x1Upper upper x1
@param x1Lower lower x1
@param x2Upper upper x2
@param x2Lower lower x2
@param gridValue1 grid value 1
@param gridValue2 grid value 2
@param functionValue function value (energy)
@param functionValue1 d(energy)/dx1
@param functionValue2 d(energy)/dx2
--------------------------------------------------------------------------------------- */
void AmoebaReferenceTorsionTorsionForce::getBicubicValues(
const RealOpenMM* y, const RealOpenMM* y1, const RealOpenMM* y2, const RealOpenMM* y12,
const RealOpenMM x1Lower, const RealOpenMM x1Upper,
const RealOpenMM x2Lower, const RealOpenMM x2Upper,
const RealOpenMM gridValue1, const RealOpenMM gridValue2,
RealOpenMM* functionValue, RealOpenMM* functionValue1, RealOpenMM* functionValue2 ){
// ---------------------------------------------------------------------------------------
// static const std::string methodName = "getBicubicValues";
static const RealOpenMM zero = 0.0;
static const RealOpenMM two = 2.0;
static const RealOpenMM three = 3.0;
// ---------------------------------------------------------------------------------------
// get coefficent matrix
RealOpenMM coefficientMatrix[4][4];
getBicubicCoefficientMatrix( y, y1, y2, y12, x1Upper-x1Lower, x2Upper-x2Lower, coefficientMatrix );
// apply coefficent matrix
RealOpenMM t = (gridValue1 - x1Lower)/(x1Upper - x1Lower);
RealOpenMM u = (gridValue2 - x2Lower)/(x2Upper - x2Lower);
*functionValue = zero;
*functionValue1 = zero;
*functionValue2 = zero;
for( int ii = 3; ii >= 0; ii-- ){
*functionValue = t*(*functionValue) + ( ( coefficientMatrix[ii][3]*u + coefficientMatrix[ii][2] )*u + coefficientMatrix[ii][1])*u + coefficientMatrix[ii][0];
*functionValue1 = u*(*functionValue1) + ( three*coefficientMatrix[3][ii]*t + two*coefficientMatrix[2][ii] )*t + coefficientMatrix[1][ii];
*functionValue2 = t*(*functionValue2) + ( three*coefficientMatrix[ii][3]*u + two*coefficientMatrix[ii][2] )*u + coefficientMatrix[ii][1];
}
*functionValue1 /= (x1Upper - x1Lower);
*functionValue2 /= (x2Upper - x2Lower);
return;
}
/**---------------------------------------------------------------------------------------
Tests the attached atoms at a torsion-torsion central
site and returns -1 if the site is chiral; else return 1
@param atomA a-atom
@param atomB b-atom
@param atomC c-atom (central atom)
@param atomD d-atom
@return 1.0 or -1.0 depending on whether chirality has an inverted sign
--------------------------------------------------------------------------------------- */
int AmoebaReferenceTorsionTorsionForce::checkTorsionSign(
const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
const RealOpenMM* positionAtomC, const RealOpenMM* positionAtomD ){
// ---------------------------------------------------------------------------------------
// static const std::string methodName = "AmoebaReferenceTorsionTorsionForce::checkTorsionSign";
static const RealOpenMM zero = 0.0;
static const int one = 1;
// ---------------------------------------------------------------------------------------
// compute parallelpiped volume at atomC and return sign based on sign of volume
enum { CA, CB, CD, LastDeltaIndex };
std::vector<RealOpenMM> deltaR[LastDeltaIndex];
for( unsigned int ii = 0; ii < LastDeltaIndex; ii++ ){
deltaR[ii].resize(3);
}
AmoebaReferenceForce::loadDeltaR( positionAtomC, positionAtomA, deltaR[CA] );
AmoebaReferenceForce::loadDeltaR( positionAtomC, positionAtomB, deltaR[CB] );
AmoebaReferenceForce::loadDeltaR( positionAtomC, positionAtomD, deltaR[CD] );
RealOpenMM volume = deltaR[CA][0]*( deltaR[CB][1]*deltaR[CD][2] - deltaR[CB][2]*deltaR[CD][1] ) +
deltaR[CB][0]*( deltaR[CD][1]*deltaR[CA][2] - deltaR[CD][2]*deltaR[CA][1] ) +
deltaR[CD][0]*( deltaR[CA][1]*deltaR[CB][2] - deltaR[CA][2]*deltaR[CB][1] );
return (volume >= zero ? one : -one);
}
/**---------------------------------------------------------------------------------------
Calculate Amoeba torsion-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 positionChiralCheckAtom Cartesian coordinates of atom to be used in chiral check;
if NULL, then no check is performed
@param grid grid vector
@param forces force vector
@return energy
--------------------------------------------------------------------------------------- */
RealOpenMM AmoebaReferenceTorsionTorsionForce::calculateForceAndEnergy( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
const RealOpenMM* positionAtomC, const RealOpenMM* positionAtomD,
const RealOpenMM* positionAtomE, const RealOpenMM* positionChiralCheckAtom,
const std::vector< std::vector< std::vector<RealOpenMM> > >& grid,
RealOpenMM** forces ){
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "AmoebaReferenceTorsionTorsionForce::calculateForceAndEnergy";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
// ---------------------------------------------------------------------------------------
enum { A, B, C, D, E, LastAtomIndex };
// get deltaR between various combinations of the 4 atoms
// and various intermediate terms
enum { BA, CB, DC, ED, T, U, V, UxV, CA, DB, EC, dT, dU, dU2, dV2, LastDeltaIndex };
std::vector<RealOpenMM> deltaR[LastDeltaIndex];
for( unsigned int ii = 0; ii < LastDeltaIndex; ii++ ){
deltaR[ii].resize(3);
}
AmoebaReferenceForce::loadDeltaR( positionAtomA, positionAtomB, deltaR[BA] );
AmoebaReferenceForce::loadDeltaR( positionAtomB, positionAtomC, deltaR[CB] );
AmoebaReferenceForce::loadDeltaR( positionAtomC, positionAtomD, deltaR[DC] );
AmoebaReferenceForce::loadDeltaR( positionAtomD, positionAtomE, deltaR[ED] );
AmoebaReferenceForce::loadDeltaR( positionAtomA, positionAtomC, deltaR[CA] );
AmoebaReferenceForce::loadDeltaR( positionAtomB, positionAtomD, deltaR[DB] );
AmoebaReferenceForce::loadDeltaR( positionAtomC, positionAtomE, deltaR[EC] );
std::vector<RealOpenMM> d[LastAtomIndex];
for( unsigned int ii = 0; ii < LastAtomIndex; ii++ ){
d[ii].resize(3);
}
AmoebaReferenceForce::getCrossProduct( deltaR[BA], deltaR[CB], deltaR[T] );
AmoebaReferenceForce::getCrossProduct( deltaR[CB], deltaR[DC], deltaR[U] );
AmoebaReferenceForce::getCrossProduct( deltaR[DC], deltaR[ED], deltaR[V] );
AmoebaReferenceForce::getCrossProduct( deltaR[U], deltaR[V], deltaR[UxV] );
RealOpenMM rT2 = AmoebaReferenceForce::getNormSquared3( deltaR[T] );
RealOpenMM rU2 = AmoebaReferenceForce::getNormSquared3( deltaR[U] );
RealOpenMM rV2 = AmoebaReferenceForce::getNormSquared3( deltaR[V] );
RealOpenMM rUrV = SQRT( rU2*rV2 );
RealOpenMM rTrU = SQRT( rT2*rU2 );
if( rTrU <= zero || rUrV <= zero ){
return zero;
}
RealOpenMM rCB = AmoebaReferenceForce::getNorm3( deltaR[CB] );
RealOpenMM cosine1 = AmoebaReferenceForce::getDotProduct3( deltaR[T], deltaR[U] );
cosine1 /= rTrU;
RealOpenMM angle1;
if( cosine1 <= -one ){
angle1 = PI_M*RADIAN;
} else if( cosine1 >= one ){
angle1 = zero;
} else {
angle1 = RADIAN*ACOS(cosine1);
}
RealOpenMM sign = AmoebaReferenceForce::getDotProduct3( deltaR[BA], deltaR[U] );
if( sign < zero ){
angle1 = -angle1;
}
// value1 = angle1;
RealOpenMM rDC = AmoebaReferenceForce::getNorm3( deltaR[DC] );
RealOpenMM cosine2 = AmoebaReferenceForce::getDotProduct3( deltaR[U], deltaR[V] );
cosine2 /= rUrV;
RealOpenMM angle2;
if( cosine2 <= -one ){
angle2 = PI_M*RADIAN;
} else if( cosine1 >= one ){
angle2 = zero;
} else {
angle2 = RADIAN*ACOS(cosine2);
}
sign = AmoebaReferenceForce::getDotProduct3( deltaR[CB], deltaR[V] );
if( sign < zero ){
angle2 = -angle2;
}
// swap signs of angles if chirality at central atom
// is 'negative'
if( positionChiralCheckAtom ){
sign = checkTorsionSign( positionChiralCheckAtom, positionAtomB, positionAtomC, positionAtomD );
if( sign < zero ){
angle1 = -angle1;
angle2 = -angle2;
}
} else {
sign = one;
}
// bicubic interpolation
RealOpenMM corners[2][2];
RealOpenMM eValues[4][4];
enum { E0, E1, E2, E12, LastEIndex };
loadGridValuesFromEnclosingRectangle( grid, angle1, angle2, corners, eValues[E0], eValues[E1], eValues[E2], eValues[E12] );
// get coordinates of point in grid closest to angle1 & angle2
// get corners of grid encompassing point
// get width/height of encompassing rectangle
RealOpenMM gridEnergy;
RealOpenMM dEdAngle1;
RealOpenMM dEdAngle2;
AmoebaReferenceTorsionTorsionForce::getBicubicValues(
eValues[E0], eValues[E1], eValues[E2], eValues[E12],
corners[0][0], corners[0][1], corners[1][0], corners[1][1],
angle1, angle2, &gridEnergy, &dEdAngle1, &dEdAngle2 );
dEdAngle1 = sign*RADIAN*dEdAngle1;
dEdAngle2 = sign*RADIAN*dEdAngle2;
AmoebaReferenceForce::getCrossProduct( deltaR[T], deltaR[CB], deltaR[dT] );
AmoebaReferenceForce::getCrossProduct( deltaR[U], deltaR[CB], deltaR[dU] );
RealOpenMM factorT = dEdAngle1/(rCB*rT2);
RealOpenMM factorU = -dEdAngle1/(rCB*rU2);
deltaR[dT][0] *= factorT;
deltaR[dT][1] *= factorT;
deltaR[dT][2] *= factorT;
deltaR[dU][0] *= factorU;
deltaR[dU][1] *= factorU;
deltaR[dU][2] *= factorU;
AmoebaReferenceForce::getCrossProduct( deltaR[dT], deltaR[CB], d[A] );
AmoebaReferenceForce::getCrossProduct( deltaR[dU], deltaR[CB], d[D] );
std::vector<RealOpenMM> tmp[3];
for( unsigned int ii = 0; ii < 3; ii++ ){
tmp[ii].resize(3);
}
AmoebaReferenceForce::getCrossProduct( deltaR[CA], deltaR[dT], d[B] );
AmoebaReferenceForce::getCrossProduct( deltaR[dU], deltaR[DC], tmp[0] );
AmoebaReferenceForce::getCrossProduct( deltaR[dT], deltaR[BA], d[C] );
AmoebaReferenceForce::getCrossProduct( deltaR[DB], deltaR[dU], tmp[1]);
d[B][0] += tmp[0][0];
d[B][1] += tmp[0][1];
d[B][2] += tmp[0][2];
d[C][0] += tmp[1][0];
d[C][1] += tmp[1][1];
d[C][2] += tmp[1][2];
// angle2
AmoebaReferenceForce::getCrossProduct( deltaR[U], deltaR[DC], deltaR[dU2] );
AmoebaReferenceForce::getCrossProduct( deltaR[V], deltaR[DC], deltaR[dV2] );
RealOpenMM factorU2 = dEdAngle2/(rDC*rU2);
RealOpenMM factorV2 = -dEdAngle2/(rDC*rV2);
deltaR[dU2][0] *= factorU2;
deltaR[dU2][1] *= factorU2;
deltaR[dU2][2] *= factorU2;
deltaR[dV2][0] *= factorV2;
deltaR[dV2][1] *= factorV2;
deltaR[dV2][2] *= factorV2;
// dB2
AmoebaReferenceForce::getCrossProduct( deltaR[dU2], deltaR[DC], tmp[0] );
// dC2
AmoebaReferenceForce::getCrossProduct( deltaR[DB], deltaR[dU2], tmp[1] );
AmoebaReferenceForce::getCrossProduct( deltaR[dV2], deltaR[ED], tmp[2] );
d[B][0] += tmp[0][0];
d[B][1] += tmp[0][1];
d[B][2] += tmp[0][2];
d[C][0] += tmp[1][0] + tmp[2][0];
d[C][1] += tmp[1][1] + tmp[2][1];
d[C][2] += tmp[1][2] + tmp[2][2];
// dD2
AmoebaReferenceForce::getCrossProduct( deltaR[dU2], deltaR[CB], tmp[0] );
AmoebaReferenceForce::getCrossProduct( deltaR[EC], deltaR[dV2], tmp[1] );
d[D][0] += tmp[0][0] + tmp[1][0];
d[D][1] += tmp[0][1] + tmp[1][1];
d[D][2] += tmp[0][2] + tmp[1][2];
// dE
AmoebaReferenceForce::getCrossProduct( deltaR[dV2], deltaR[DC], d[E] );
// ---------------------------------------------------------------------------------------
// add in forces
for( int jj = 0; jj < LastAtomIndex; jj++ ){
forces[jj][0] -= d[jj][0];
forces[jj][1] -= d[jj][1];
forces[jj][2] -= d[jj][2];
}
// ---------------------------------------------------------------------------------------
return gridEnergy;
}
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceTorsionTorsionForce.h 0 → 100644
View file @ 8a57342c
/* 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 __AmoebaReferenceTorsionTorsionForce_H__
#define __AmoebaReferenceTorsionTorsionForce_H__
#include "SimTKUtilities/SimTKOpenMMRealType.h"
// ---------------------------------------------------------------------------------------
class AmoebaReferenceTorsionTorsionForce {
public:
/**---------------------------------------------------------------------------------------
Constructor
--------------------------------------------------------------------------------------- */
AmoebaReferenceTorsionTorsionForce( );
/**---------------------------------------------------------------------------------------
Destructor
--------------------------------------------------------------------------------------- */
~AmoebaReferenceTorsionTorsionForce( );
/**---------------------------------------------------------------------------------------
Load grid values from rectangle enclosing angles
@param angle1 angle in first dimension
@param angle2 angle in second dimension
@param corners on return contains the coordinates of the rectangle
@param fValues on return contains the values of the function at the vertices of the rectangle
@param fValues1 on return contains the values of the derivative of the function wrt first dimension
@param fValues2 on return contains the values of the derivative of the function wrt second dimension
@param fValues12 on return contains the values of the derivative of the function wrt first & second dimension
On first call a check is performed to see if the grid is valid
--------------------------------------------------------------------------------------- */
static void loadGridValuesFromEnclosingRectangle(
const std::vector< std::vector< std::vector<RealOpenMM> > >& grid,
RealOpenMM angle1, RealOpenMM angle2, RealOpenMM corners[2][2],
RealOpenMM* fValues, RealOpenMM* fValues1, RealOpenMM* fValues2, RealOpenMM* fValues12 );
/**---------------------------------------------------------------------------------------
Determines the coefficient matrix needed for bicubic
interpolation of a function, gradients and cross derivatives
Reference:
W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P.
Flannery, Numerical Recipes (Fortran), 2nd Ed., Cambridge
University Press, 1992, Section 3.6
@param y y values
@param y1 dy/dx1 values
@param y2 dy/dx2 values
@param y12 d2y/dx1dx2 values
@param d1 d1Upper - d1Lower
@param d2 d2Upper - d2Lower
@param c 4x4 return coefficient matrix
--------------------------------------------------------------------------------------- */
static void getBicubicCoefficientMatrix( const RealOpenMM* y,
const RealOpenMM* y1, const RealOpenMM* y2, const RealOpenMM* y12, const RealOpenMM d1, const RealOpenMM d2,
RealOpenMM c[4][4] );
/**---------------------------------------------------------------------------------------
Determines the coefficient matrix needed for bicubic
interpolation of a function, gradients and cross derivatives
Reference:
W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P.
Flannery, Numerical Recipes (Fortran), 2nd Ed., Cambridge
University Press, 1992, Section 3.6
@param y y values
@param y1 dy/dx1 values
@param y2 dy/dx2 values
@param y12 d2y/dx1dx2 values
@param x1Upper upper x1
@param x1Lower lower x1
@param x2Upper upper x2
@param x2Lower lower x2
@param gridValue1 grid value 1
@param gridValue2 grid value 2
@param functionValue function value (energy)
@param functionValue1 d(energy)/dx1
@param functionValue2 d(energy)/dx2
@return AmoebaCommon::DefaultReturn
--------------------------------------------------------------------------------------- */
static void getBicubicValues(
const RealOpenMM* y, const RealOpenMM* y1, const RealOpenMM* y2, const RealOpenMM* y12,
const RealOpenMM x1Lower, const RealOpenMM x1Upper,
const RealOpenMM x2Lower, const RealOpenMM x2Upper,
const RealOpenMM gridValue1, const RealOpenMM gridValue2,
RealOpenMM* functionValue, RealOpenMM* functionValue1, RealOpenMM* functionValue2 );
/**---------------------------------------------------------------------------------------
Tests the attached atoms at a torsion-torsion central
site and returns -1 if the site is chiral; else return 1
@param atomA a-atom
@param atomB b-atom
@param atomC c-atom (central atom)
@param atomD d-atom
@return 1.0 or -1.0 depending on whether chirality has an inverted sign
--------------------------------------------------------------------------------------- */
static int checkTorsionSign(
const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
const RealOpenMM* positionAtomC, const RealOpenMM* positionAtomD );
/**---------------------------------------------------------------------------------------
Calculate Amoeba torsion-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 positionChiralCheckAtom Cartesian coordinates of atom to be used in chiral check;
if NULL, then no check is performed
@param grid grid vector
@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* chiralCheckAtom,
const std::vector< std::vector< std::vector<RealOpenMM> > >& grid,
RealOpenMM** forces );
};
// ---------------------------------------------------------------------------------------
#endif // _AmoebaReferenceTorsionTorsionForce___
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