Skip to content

GitLab

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

Revamp of reference local force classes

parent 3b67aadf
Hide whitespace changes
Inline Side-by-side
Showing with 686 additions and 285 deletions
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.cpp
View file @ a37e41ea
......@@ -97,19 +97,10 @@ void ReferenceCalcAmoebaHarmonicBondForceKernel::initialize(const System& system
double ReferenceCalcAmoebaHarmonicBondForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
RealOpenMM** posData = extractPositions(context);
RealOpenMM** forceData = extractForces(context);
RealOpenMM energy = 0.0;
for (int ii = 0; ii < numBonds; ii++) {
int particle1Index = particle1[ii];
int particle2Index = particle2[ii];
RealOpenMM bondLength = length[ii];
RealOpenMM bondK = kQuadratic[ii];
RealOpenMM* forces[2];
forces[0] = forceData[particle1Index];
forces[1] = forceData[particle2Index];
energy += AmoebaReferenceHarmonicBondForce::calculateForceAndEnergy( posData[particle1Index], posData[particle2Index],
bondLength, bondK, globalHarmonicBondCubic, globalHarmonicBondQuartic,
forces );
}
AmoebaReferenceHarmonicBondForce amoebaReferenceHarmonicBondForce;
RealOpenMM energy = amoebaReferenceHarmonicBondForce.calculateForceAndEnergy( numBonds, posData, particle1, particle2, length, kQuadratic,
globalHarmonicBondCubic, globalHarmonicBondQuartic,
forceData );
return static_cast<double>(energy);
}
......@@ -143,22 +134,9 @@ void ReferenceCalcAmoebaHarmonicAngleForceKernel::initialize(const System& syste
double ReferenceCalcAmoebaHarmonicAngleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
RealOpenMM** posData = extractPositions(context);
RealOpenMM** forceData = extractForces(context);
RealOpenMM energy = 0.0;
for (unsigned int ii = 0; ii < numAngles; ii++) {
int particle1Index = particle1[ii];
int particle2Index = particle2[ii];
int particle3Index = particle3[ii];
RealOpenMM idealAngle = angle[ii];
RealOpenMM angleK = kQuadratic[ii];
RealOpenMM* forces[3];
forces[0] = forceData[particle1Index];
forces[1] = forceData[particle2Index];
forces[2] = forceData[particle3Index];
energy += AmoebaReferenceHarmonicAngleForce::calculateForceAndEnergy(
posData[particle1Index], posData[particle2Index], posData[particle3Index],
idealAngle, angleK, globalHarmonicAngleCubic, globalHarmonicAngleQuartic,
globalHarmonicAnglePentic, globalHarmonicAngleSextic, forces );
}
AmoebaReferenceHarmonicAngleForce amoebaReferenceHarmonicAngleForce;
RealOpenMM energy = amoebaReferenceHarmonicAngleForce.calculateForceAndEnergy( numAngles,
posData, particle1, particle2, particle3, angle, kQuadratic, globalHarmonicAngleCubic, globalHarmonicAngleQuartic, globalHarmonicAnglePentic, globalHarmonicAngleSextic, forceData );
return static_cast<double>(energy);
}
......@@ -190,26 +168,13 @@ void ReferenceCalcAmoebaHarmonicInPlaneAngleForceKernel::initialize(const System
}
double ReferenceCalcAmoebaHarmonicInPlaneAngleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
RealOpenMM** posData = extractPositions(context);
RealOpenMM** forceData = extractForces(context);
RealOpenMM energy = 0.0;
for (unsigned int ii = 0; ii < numAngles; ii++) {
int particle1Index = particle1[ii];
int particle2Index = particle2[ii];
int particle3Index = particle3[ii];
int particle4Index = particle4[ii];
RealOpenMM idealAngle = angle[ii];
RealOpenMM angleK = kQuadratic[ii];
RealOpenMM* forces[4];
forces[0] = forceData[particle1Index];
forces[1] = forceData[particle2Index];
forces[2] = forceData[particle3Index];
forces[3] = forceData[particle4Index];
energy += AmoebaReferenceHarmonicInPlaneAngleForce::calculateForceAndEnergy(
posData[particle1Index], posData[particle2Index], posData[particle3Index], posData[particle4Index],
idealAngle, angleK, globalHarmonicInPlaneAngleCubic, globalHarmonicInPlaneAngleQuartic,
globalHarmonicInPlaneAnglePentic, globalHarmonicInPlaneAngleSextic, forces );
}
AmoebaReferenceHarmonicInPlaneAngleForce amoebaReferenceHarmonicInPlaneAngleForce;
RealOpenMM energy = amoebaReferenceHarmonicInPlaneAngleForce.calculateForceAndEnergy( numAngles, posData, particle1, particle2, particle3, particle4,
angle, kQuadratic, globalHarmonicInPlaneAngleCubic, globalHarmonicInPlaneAngleQuartic,
globalHarmonicInPlaneAnglePentic, globalHarmonicInPlaneAngleSextic, forceData );
return static_cast<double>(energy);
}
......@@ -256,24 +221,9 @@ void ReferenceCalcAmoebaTorsionForceKernel::initialize(const System& system, con
double ReferenceCalcAmoebaTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
RealOpenMM** posData = extractPositions(context);
RealOpenMM** forceData = extractForces(context);
RealOpenMM energy = 0.0;
for (unsigned int ii = 0; ii < numTorsions; ii++) {
int particle1Index = particle1[ii];
int particle2Index = particle2[ii];
int particle3Index = particle3[ii];
int particle4Index = particle4[ii];
RealOpenMM* forces[4];
forces[0] = forceData[particle1Index];
forces[1] = forceData[particle2Index];
forces[2] = forceData[particle3Index];
forces[3] = forceData[particle4Index];
energy += AmoebaReferenceTorsionForce::calculateForceAndEnergy(
posData[particle1Index], posData[particle2Index], posData[particle3Index], posData[particle4Index],
torsionParameters1[ii], torsionParameters2[ii], torsionParameters3[ii], forces );
}
AmoebaReferenceTorsionForce amoebaReferenceTorsionForce;
RealOpenMM energy = amoebaReferenceTorsionForce.calculateForceAndEnergy( numTorsions, posData, particle1, particle2, particle3, particle4,
torsionParameters1, torsionParameters2, torsionParameters3, forceData );
return static_cast<double>(energy);
}
......@@ -305,28 +255,10 @@ void ReferenceCalcAmoebaPiTorsionForceKernel::initialize(const System& system, c
double ReferenceCalcAmoebaPiTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
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 );
}
AmoebaReferencePiTorsionForce amoebaReferencePiTorsionForce;
RealOpenMM energy = amoebaReferencePiTorsionForce.calculateForceAndEnergy( numPiTorsions, posData, particle1, particle2,
particle3, particle4, particle5, particle6,
kTorsion, forceData );
return static_cast<double>(energy);
}
......@@ -357,23 +289,9 @@ void ReferenceCalcAmoebaStretchBendForceKernel::initialize(const System& system,
double ReferenceCalcAmoebaStretchBendForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
RealOpenMM** posData = extractPositions(context);
RealOpenMM** forceData = extractForces(context);
RealOpenMM energy = 0.0;
for( unsigned int ii = 0; ii < numStretchBends; ii++ ){
int particle1Index = particle1[ii];
int particle2Index = particle2[ii];
int particle3Index = particle3[ii];
RealOpenMM* forces[3];
forces[0] = forceData[particle1Index];
forces[1] = forceData[particle2Index];
forces[2] = forceData[particle3Index];
energy += AmoebaReferenceStretchBendForce::calculateForceAndEnergy(
posData[particle1Index], posData[particle2Index], posData[particle3Index],
lengthABParameters[ii], lengthCBParameters[ii],
angleParameters[ii], kParameters[ii], forces );
}
AmoebaReferenceStretchBendForce amoebaReferenceStretchBendForce;
RealOpenMM energy = amoebaReferenceStretchBendForce.calculateForceAndEnergy( numStretchBends, posData, particle1, particle2, particle3,
lengthABParameters, lengthCBParameters, angleParameters, kParameters, forceData );
return static_cast<double>(energy);
}
......@@ -409,23 +327,14 @@ void ReferenceCalcAmoebaOutOfPlaneBendForceKernel::initialize(const System& syst
double ReferenceCalcAmoebaOutOfPlaneBendForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
RealOpenMM** posData = extractPositions(context);
RealOpenMM** forceData = extractForces(context);
RealOpenMM energy = 0.0;
for (unsigned int ii = 0; ii < numOutOfPlaneBends; ii++) {
int particle1Index = particle1[ii];
int particle2Index = particle2[ii];
int particle3Index = particle3[ii];
int particle4Index = particle4[ii];
RealOpenMM angleK = kParameters[ii];
RealOpenMM* forces[4];
forces[0] = forceData[particle1Index];
forces[1] = forceData[particle2Index];
forces[2] = forceData[particle3Index];
forces[3] = forceData[particle4Index];
energy += AmoebaReferenceOutOfPlaneBendForce::calculateForceAndEnergy(
posData[particle1Index], posData[particle2Index], posData[particle3Index], posData[particle4Index],
angleK, globalOutOfPlaneBendAngleCubic, globalOutOfPlaneBendAngleQuartic,
globalOutOfPlaneBendAnglePentic, globalOutOfPlaneBendAngleSextic, forces );
}
AmoebaReferenceOutOfPlaneBendForce amoebaReferenceOutOfPlaneBendForce;
RealOpenMM energy = amoebaReferenceOutOfPlaneBendForce.calculateForceAndEnergy( numOutOfPlaneBends, posData,
particle1, particle2, particle3, particle4,
kParameters,
globalOutOfPlaneBendAngleCubic,
globalOutOfPlaneBendAngleQuartic,
globalOutOfPlaneBendAnglePentic,
globalOutOfPlaneBendAngleSextic, forceData );
return static_cast<double>(energy);
}
......@@ -482,39 +391,10 @@ double ReferenceCalcAmoebaTorsionTorsionForceKernel::execute(ContextImpl& contex
RealOpenMM** posData = extractPositions(context);
RealOpenMM** forceData = extractForces(context);
RealOpenMM energy = 0.0;
for( unsigned int ii = 0; ii < numTorsionTorsions; ii++ ){
int particle1Index = particle1[ii];
int particle2Index = particle2[ii];
int particle3Index = particle3[ii];
int particle4Index = particle4[ii];
int particle5Index = particle5[ii];
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 );
}
AmoebaReferenceTorsionTorsionForce amoebaReferenceTorsionTorsionForce;
RealOpenMM energy = amoebaReferenceTorsionTorsionForce.calculateForceAndEnergy( numTorsionTorsions, posData,
particle1, particle2, particle3, particle4, particle5,
chiralCheckAtom, gridIndices, torsionTorsionGrids, forceData );
return static_cast<double>(energy);
}
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceHarmonicAngleForce.cpp
View file @ a37e41ea
......@@ -24,7 +24,6 @@
#include "AmoebaReferenceForce.h"
#include "AmoebaReferenceHarmonicAngleForce.h"
#include <vector>
/**---------------------------------------------------------------------------------------
......@@ -48,7 +47,7 @@ RealOpenMM AmoebaReferenceHarmonicAngleForce::getPrefactorsGivenAngleCosine( Rea
RealOpenMM idealAngle, RealOpenMM angleK,
RealOpenMM angleCubic, RealOpenMM angleQuartic,
RealOpenMM anglePentic, RealOpenMM angleSextic,
RealOpenMM* dEdR ){
RealOpenMM* dEdR ) const {
// ---------------------------------------------------------------------------------------
......@@ -111,12 +110,12 @@ RealOpenMM AmoebaReferenceHarmonicAngleForce::getPrefactorsGivenAngleCosine( Rea
--------------------------------------------------------------------------------------- */
RealOpenMM AmoebaReferenceHarmonicAngleForce::calculateForceAndEnergy( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
const RealOpenMM* positionAtomC,
RealOpenMM angle, RealOpenMM angleK,
RealOpenMM angleCubic, RealOpenMM angleQuartic,
RealOpenMM anglePentic, RealOpenMM angleSextic,
RealOpenMM** forces ){
RealOpenMM AmoebaReferenceHarmonicAngleForce::calculateAngleIxn( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
const RealOpenMM* positionAtomC,
RealOpenMM angle, RealOpenMM angleK,
RealOpenMM angleCubic, RealOpenMM angleQuartic,
RealOpenMM anglePentic, RealOpenMM angleSextic,
RealOpenMM** forces ) const {
// ---------------------------------------------------------------------------------------
......@@ -181,3 +180,30 @@ RealOpenMM AmoebaReferenceHarmonicAngleForce::calculateForceAndEnergy( const Rea
return energy;
}
RealOpenMM AmoebaReferenceHarmonicAngleForce::calculateForceAndEnergy( int numAngles, RealOpenMM** posData,
const std::vector<int>& particle1,
const std::vector<int>& particle2,
const std::vector<int>& particle3,
const std::vector<RealOpenMM>& angle,
const std::vector<RealOpenMM>& kQuadratic,
RealOpenMM angleCubic,
RealOpenMM angleQuartic,
RealOpenMM anglePentic,
RealOpenMM angleSextic,
RealOpenMM** forceData) const {
RealOpenMM energy = 0.0;
for (unsigned int ii = 0; ii < numAngles; ii++) {
int particle1Index = particle1[ii];
int particle2Index = particle2[ii];
int particle3Index = particle3[ii];
RealOpenMM idealAngle = angle[ii];
RealOpenMM angleK = kQuadratic[ii];
RealOpenMM* forces[3];
forces[0] = forceData[particle1Index];
forces[1] = forceData[particle2Index];
forces[2] = forceData[particle3Index];
energy += calculateAngleIxn( posData[particle1Index], posData[particle2Index], posData[particle3Index],
idealAngle, angleK, angleCubic, angleQuartic, anglePentic, angleSextic, forces );
}
return energy;
}
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceHarmonicAngleForce.h
View file @ a37e41ea
......@@ -26,6 +26,7 @@
#define __AmoebaReferenceHarmonicAngleForce_H__
#include "SimTKUtilities/SimTKOpenMMRealType.h"
#include <vector>
// ---------------------------------------------------------------------------------------
......@@ -39,7 +40,7 @@ public:
--------------------------------------------------------------------------------------- */
AmoebaReferenceHarmonicAngleForce( );
AmoebaReferenceHarmonicAngleForce( ){};
/**---------------------------------------------------------------------------------------
......@@ -47,8 +48,44 @@ public:
--------------------------------------------------------------------------------------- */
~AmoebaReferenceHarmonicAngleForce( );
~AmoebaReferenceHarmonicAngleForce( ){};
/**---------------------------------------------------------------------------------------
Calculate Amoeba harmonic angle ixns (force and energy)
@param numBonds number of angles
@param posData particle positions
@param particle1 particle 1 indices
@param particle2 particle 2 indices
@param particle3 particle 3 indices
@param angle ideal angle
@param angleK angle force constant
@param angleCubic cubic force parameter
@param angleQuartic quartic force parameter
@param anglePentic pentic force parameter
@param angleSexic sextic force parameter
@param forces output force vector
@return total energy
--------------------------------------------------------------------------------------- */
RealOpenMM calculateForceAndEnergy( int numAngles, RealOpenMM** posData,
const std::vector<int>& particle1,
const std::vector<int>& particle2,
const std::vector<int>& particle3,
const std::vector<RealOpenMM>& angle,
const std::vector<RealOpenMM>& kQuadratic,
RealOpenMM globalHarmonicAngleCubic,
RealOpenMM globalHarmonicAngleQuartic,
RealOpenMM globalHarmonicAnglePentic,
RealOpenMM globalHarmonicAngleSextic,
RealOpenMM** forceData ) const;
private:
/**---------------------------------------------------------------------------------------
Get dEdT and energy prefactor given cosine of angle :: the calculation for different
......@@ -67,10 +104,10 @@ public:
--------------------------------------------------------------------------------------- */
static RealOpenMM getPrefactorsGivenAngleCosine( RealOpenMM cosine, RealOpenMM idealAngle, RealOpenMM angleK,
RealOpenMM angleCubic, RealOpenMM angleQuartic,
RealOpenMM anglePentic, RealOpenMM angleSextic,
RealOpenMM* dEdR );
RealOpenMM getPrefactorsGivenAngleCosine( RealOpenMM cosine, RealOpenMM idealAngle, RealOpenMM angleK,
RealOpenMM angleCubic, RealOpenMM angleQuartic,
RealOpenMM anglePentic, RealOpenMM angleSextic,
RealOpenMM* dEdR ) const;
/**---------------------------------------------------------------------------------------
......@@ -91,12 +128,12 @@ public:
--------------------------------------------------------------------------------------- */
static RealOpenMM calculateForceAndEnergy( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
const RealOpenMM* positionAtomC,
RealOpenMM angle, RealOpenMM angleK,
RealOpenMM angleCubic, RealOpenMM angleQuartic,
RealOpenMM anglePentic, RealOpenMM angleSextic,
RealOpenMM** forces );
RealOpenMM calculateAngleIxn( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
const RealOpenMM* positionAtomC,
RealOpenMM angle, RealOpenMM angleK,
RealOpenMM angleCubic, RealOpenMM angleQuartic,
RealOpenMM anglePentic, RealOpenMM angleSextic,
RealOpenMM** forces ) const;
};
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceHarmonicBondForce.cpp
View file @ a37e41ea
......@@ -24,7 +24,6 @@
#include "AmoebaReferenceHarmonicBondForce.h"
#include "AmoebaReferenceForce.h"
#include <vector>
/**---------------------------------------------------------------------------------------
......@@ -42,10 +41,10 @@
--------------------------------------------------------------------------------------- */
RealOpenMM AmoebaReferenceHarmonicBondForce::calculateForceAndEnergy( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
RealOpenMM bondLength, RealOpenMM bondK,
RealOpenMM bondCubic, RealOpenMM bondQuartic,
RealOpenMM** forces ){
RealOpenMM AmoebaReferenceHarmonicBondForce::calculateBondIxn( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
RealOpenMM bondLength, RealOpenMM bondK,
RealOpenMM bondCubic, RealOpenMM bondQuartic,
RealOpenMM** forces ) const {
// ---------------------------------------------------------------------------------------
......@@ -85,3 +84,27 @@ RealOpenMM AmoebaReferenceHarmonicBondForce::calculateForceAndEnergy( const Real
return energy;
}
RealOpenMM AmoebaReferenceHarmonicBondForce::calculateForceAndEnergy( int numBonds,
RealOpenMM** particlePositions,
const std::vector<int>& particle1,
const std::vector<int>& particle2,
const std::vector<RealOpenMM>& length,
const std::vector<RealOpenMM>& kQuadratic,
RealOpenMM globalHarmonicBondCubic,
RealOpenMM globalHarmonicBondQuartic,
RealOpenMM** forceData ) const {
RealOpenMM energy = 0.0;
for( int ii = 0; ii < numBonds; ii++ ){
int particle1Index = particle1[ii];
int particle2Index = particle2[ii];
RealOpenMM bondLength = length[ii];
RealOpenMM bondK = kQuadratic[ii];
RealOpenMM* forces[2];
forces[0] = forceData[particle1Index];
forces[1] = forceData[particle2Index];
energy += calculateBondIxn( particlePositions[particle1Index], particlePositions[particle2Index],
bondLength, bondK, globalHarmonicBondCubic, globalHarmonicBondQuartic,
forces );
}
return energy;
}
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceHarmonicBondForce.h
View file @ a37e41ea
......@@ -26,6 +26,7 @@
#define __AmoebaReferenceHarmonicBondForce_H__
#include "SimTKUtilities/SimTKOpenMMRealType.h"
#include <vector>
// ---------------------------------------------------------------------------------------
......@@ -47,12 +48,40 @@ public:
--------------------------------------------------------------------------------------- */
~AmoebaReferenceHarmonicBondForce( );
~AmoebaReferenceHarmonicBondForce( ){};
/**---------------------------------------------------------------------------------------
Calculate Amoeba harmonic bond ixn (force and energy)
Calculate Amoeba harmonic bond ixns (force and energy)
@param numBonds number of bonds
@param posData particle positions
@param particle1 particle 1 indices
@param particle2 particle 2 indices
@param bondLength bond length
@param bondK bond force
@param bondCubic cubic bond force parameter
@param bondQuartic quartic bond force parameter
@param forces output force vector
@return total energy
--------------------------------------------------------------------------------------- */
RealOpenMM calculateForceAndEnergy( int numBonds, RealOpenMM** posData,
const std::vector<int>& particle1,
const std::vector<int>& particle2,
const std::vector<RealOpenMM>& bondLength,
const std::vector<RealOpenMM>& bondK,
RealOpenMM bondCubic, RealOpenMM bondQuartic,
RealOpenMM** forceData ) const;
private:
/**---------------------------------------------------------------------------------------
Calculate Amoeba harmonic bond ixns (force and energy)
@param positionAtomA Cartesian coordinates of atom A
@param positionAtomB Cartesian coordinates of atom B
......@@ -66,10 +95,10 @@ public:
--------------------------------------------------------------------------------------- */
static RealOpenMM calculateForceAndEnergy( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
RealOpenMM bondLength, RealOpenMM bondK,
RealOpenMM bondCubic, RealOpenMM bondQuartic,
RealOpenMM** forces );
RealOpenMM calculateBondIxn( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
RealOpenMM bondLength, RealOpenMM bondK,
RealOpenMM bondCubic, RealOpenMM bondQuartic,
RealOpenMM** forces ) const;
};
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceHarmonicInPlaneAngleForce.cpp
View file @ a37e41ea
......@@ -24,7 +24,6 @@
#include "AmoebaReferenceForce.h"
#include "AmoebaReferenceHarmonicInPlaneAngleForce.h"
#include <vector>
/**---------------------------------------------------------------------------------------
......@@ -48,7 +47,7 @@ RealOpenMM AmoebaReferenceHarmonicInPlaneAngleForce::getPrefactorsGivenAngleCosi
RealOpenMM idealAngle, RealOpenMM angleK,
RealOpenMM angleCubic, RealOpenMM angleQuartic,
RealOpenMM anglePentic, RealOpenMM angleSextic,
RealOpenMM* dEdR ){
RealOpenMM* dEdR ) const {
// ---------------------------------------------------------------------------------------
......@@ -112,12 +111,12 @@ RealOpenMM AmoebaReferenceHarmonicInPlaneAngleForce::getPrefactorsGivenAngleCosi
--------------------------------------------------------------------------------------- */
RealOpenMM AmoebaReferenceHarmonicInPlaneAngleForce::calculateForceAndEnergy( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
const RealOpenMM* positionAtomC, const RealOpenMM* positionAtomD,
RealOpenMM angle, RealOpenMM angleK,
RealOpenMM angleCubic, RealOpenMM angleQuartic,
RealOpenMM anglePentic, RealOpenMM angleSextic,
RealOpenMM** forces ){
RealOpenMM AmoebaReferenceHarmonicInPlaneAngleForce::calculateAngleIxn( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
const RealOpenMM* positionAtomC, const RealOpenMM* positionAtomD,
RealOpenMM angle, RealOpenMM angleK,
RealOpenMM angleCubic, RealOpenMM angleQuartic,
RealOpenMM anglePentic, RealOpenMM angleSextic,
RealOpenMM** forces ) const {
// ---------------------------------------------------------------------------------------
......@@ -242,3 +241,35 @@ RealOpenMM AmoebaReferenceHarmonicInPlaneAngleForce::calculateForceAndEnergy( co
return energy;
}
RealOpenMM AmoebaReferenceHarmonicInPlaneAngleForce::calculateForceAndEnergy( int numAngles, RealOpenMM** posData,
const std::vector<int>& particle1,
const std::vector<int>& particle2,
const std::vector<int>& particle3,
const std::vector<int>& particle4,
const std::vector<RealOpenMM>& angle,
const std::vector<RealOpenMM>& kQuadratic,
RealOpenMM angleCubic,
RealOpenMM angleQuartic,
RealOpenMM anglePentic,
RealOpenMM angleSextic,
RealOpenMM** forceData) const {
RealOpenMM energy = 0.0;
for (unsigned int ii = 0; ii < numAngles; ii++) {
int particle1Index = particle1[ii];
int particle2Index = particle2[ii];
int particle3Index = particle3[ii];
int particle4Index = particle4[ii];
RealOpenMM idealAngle = angle[ii];
RealOpenMM angleK = kQuadratic[ii];
RealOpenMM* forces[4];
forces[0] = forceData[particle1Index];
forces[1] = forceData[particle2Index];
forces[2] = forceData[particle3Index];
forces[3] = forceData[particle4Index];
energy += calculateAngleIxn( posData[particle1Index], posData[particle2Index], posData[particle3Index], posData[particle4Index],
idealAngle, angleK, angleCubic, angleQuartic, anglePentic, angleSextic, forces );
}
return energy;
}
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceHarmonicInPlaneAngleForce.h
View file @ a37e41ea
......@@ -26,6 +26,7 @@
#define __AmoebaReferenceHarmonicInPlaneAngleForce_H__
#include "SimTKUtilities/SimTKOpenMMRealType.h"
#include <vector>
// ---------------------------------------------------------------------------------------
......@@ -39,7 +40,7 @@ public:
--------------------------------------------------------------------------------------- */
AmoebaReferenceHarmonicInPlaneAngleForce( );
AmoebaReferenceHarmonicInPlaneAngleForce( ){};
/**---------------------------------------------------------------------------------------
......@@ -47,8 +48,46 @@ public:
--------------------------------------------------------------------------------------- */
~AmoebaReferenceHarmonicInPlaneAngleForce( );
~AmoebaReferenceHarmonicInPlaneAngleForce( ){};
/**---------------------------------------------------------------------------------------
Calculate Amoeba harmonic in-plane angle ixns (force and energy)
@param numBonds number of angles
@param posData particle positions
@param particle1 particle 1 indices
@param particle2 particle 2 indices
@param particle3 particle 3 indices
@param particle4 particle 4 indices
@param angle ideal angle
@param angleK angle force constant
@param angleCubic cubic force parameter
@param angleQuartic quartic force parameter
@param anglePentic pentic force parameter
@param angleSexic sextic force parameter
@param forces output force vector
@return total energy
--------------------------------------------------------------------------------------- */
RealOpenMM calculateForceAndEnergy( int numAngles, RealOpenMM** posData,
const std::vector<int>& particle1,
const std::vector<int>& particle2,
const std::vector<int>& particle3,
const std::vector<int>& particle4,
const std::vector<RealOpenMM>& angle,
const std::vector<RealOpenMM>& kQuadratic,
RealOpenMM globalHarmonicAngleCubic,
RealOpenMM globalHarmonicAngleQuartic,
RealOpenMM globalHarmonicAnglePentic,
RealOpenMM globalHarmonicAngleSextic,
RealOpenMM** forceData ) const;
private:
/**---------------------------------------------------------------------------------------
Get dEdT and energy prefactor given cosine of angle :: the calculation for different
......@@ -67,10 +106,10 @@ public:
--------------------------------------------------------------------------------------- */
static RealOpenMM getPrefactorsGivenAngleCosine( RealOpenMM cosine, RealOpenMM idealAngle, RealOpenMM angleK,
RealOpenMM angleCubic, RealOpenMM angleQuartic,
RealOpenMM anglePentic, RealOpenMM angleSextic,
RealOpenMM* dEdR );
RealOpenMM getPrefactorsGivenAngleCosine( RealOpenMM cosine, RealOpenMM idealAngle, RealOpenMM angleK,
RealOpenMM angleCubic, RealOpenMM angleQuartic,
RealOpenMM anglePentic, RealOpenMM angleSextic,
RealOpenMM* dEdR ) const;
/**---------------------------------------------------------------------------------------
......@@ -92,12 +131,12 @@ public:
--------------------------------------------------------------------------------------- */
static RealOpenMM calculateForceAndEnergy( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
const RealOpenMM* positionAtomC, const RealOpenMM* positionAtomD,
RealOpenMM angle, RealOpenMM angleK,
RealOpenMM angleCubic, RealOpenMM angleQuartic,
RealOpenMM anglePentic, RealOpenMM angleSextic,
RealOpenMM** forces );
RealOpenMM calculateAngleIxn( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
const RealOpenMM* positionAtomC, const RealOpenMM* positionAtomD,
RealOpenMM angle, RealOpenMM angleK,
RealOpenMM angleCubic, RealOpenMM angleQuartic,
RealOpenMM anglePentic, RealOpenMM angleSextic,
RealOpenMM** forces ) const;
};
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceOutOfPlaneBendForce.cpp
View file @ a37e41ea
......@@ -24,11 +24,10 @@
#include "AmoebaReferenceForce.h"
#include "AmoebaReferenceOutOfPlaneBendForce.h"
#include <vector>
/**---------------------------------------------------------------------------------------
Calculate Amoeba harmonic angle ixn (force and energy)
Calculate Amoeba Out-Of-Plane-Bend ixn (force and energy)
@param positionAtomA Cartesian coordinates of atom A
@param positionAtomB Cartesian coordinates of atom B
......@@ -46,12 +45,12 @@
--------------------------------------------------------------------------------------- */
RealOpenMM AmoebaReferenceOutOfPlaneBendForce::calculateForceAndEnergy( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
const RealOpenMM* positionAtomC, const RealOpenMM* positionAtomD,
RealOpenMM angleK,
RealOpenMM angleCubic, RealOpenMM angleQuartic,
RealOpenMM anglePentic, RealOpenMM angleSextic,
RealOpenMM** forces ){
RealOpenMM AmoebaReferenceOutOfPlaneBendForce::calculateOutOfPlaneBendIxn( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
const RealOpenMM* positionAtomC, const RealOpenMM* positionAtomD,
RealOpenMM angleK,
RealOpenMM angleCubic, RealOpenMM angleQuartic,
RealOpenMM anglePentic, RealOpenMM angleSextic,
RealOpenMM** forces ) const {
// ---------------------------------------------------------------------------------------
......@@ -193,3 +192,33 @@ RealOpenMM AmoebaReferenceOutOfPlaneBendForce::calculateForceAndEnergy( const Re
return energy;
}
RealOpenMM AmoebaReferenceOutOfPlaneBendForce::calculateForceAndEnergy( int numOutOfPlaneBends, RealOpenMM** posData,
const std::vector<int>& particle1,
const std::vector<int>& particle2,
const std::vector<int>& particle3,
const std::vector<int>& particle4,
const std::vector<RealOpenMM>& kQuadratic,
RealOpenMM angleCubic,
RealOpenMM angleQuartic,
RealOpenMM anglePentic,
RealOpenMM angleSextic,
RealOpenMM** forceData) const {
RealOpenMM energy = 0.0;
for (unsigned int ii = 0; ii < numOutOfPlaneBends; ii++) {
int particle1Index = particle1[ii];
int particle2Index = particle2[ii];
int particle3Index = particle3[ii];
int particle4Index = particle4[ii];
RealOpenMM kAngle = kQuadratic[ii];
RealOpenMM* forces[4];
forces[0] = forceData[particle1Index];
forces[1] = forceData[particle2Index];
forces[2] = forceData[particle3Index];
forces[3] = forceData[particle4Index];
energy += calculateOutOfPlaneBendIxn( posData[particle1Index], posData[particle2Index], posData[particle3Index], posData[particle4Index],
kAngle, angleCubic, angleQuartic, anglePentic, angleSextic, forces );
}
return energy;
}
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceOutOfPlaneBendForce.h
View file @ a37e41ea
......@@ -26,6 +26,7 @@
#define __AmoebaReferenceOutOfPlaneBendForce_H__
#include "SimTKUtilities/SimTKOpenMMRealType.h"
#include <vector>
// ---------------------------------------------------------------------------------------
......@@ -39,7 +40,7 @@ public:
--------------------------------------------------------------------------------------- */
AmoebaReferenceOutOfPlaneBendForce( );
AmoebaReferenceOutOfPlaneBendForce( ){};
/**---------------------------------------------------------------------------------------
......@@ -47,11 +48,47 @@ public:
--------------------------------------------------------------------------------------- */
~AmoebaReferenceOutOfPlaneBendForce( );
~AmoebaReferenceOutOfPlaneBendForce( ){};
/**---------------------------------------------------------------------------------------
Calculate Amoeba out-of-plane-bend angle (force and energy)
@param numOutOfPlaneBends number of angles
@param posData particle positions
@param particle1 particle 1 indices
@param particle2 particle 2 indices
@param particle3 particle 3 indices
@param particle4 particle 4 indices
@param kAngle angle force constant
@param angleCubic cubic force parameter
@param angleQuartic quartic force parameter
@param anglePentic pentic force parameter
@param angleSexic sextic force parameter
@param forces output force vector
@return total energy
--------------------------------------------------------------------------------------- */
RealOpenMM calculateForceAndEnergy( int numOutOfPlaneBends, RealOpenMM** posData,
const std::vector<int>& particle1,
const std::vector<int>& particle2,
const std::vector<int>& particle3,
const std::vector<int>& particle4,
const std::vector<RealOpenMM>& kAngle,
RealOpenMM angleCubic,
RealOpenMM angleQuartic,
RealOpenMM anglePentic,
RealOpenMM angleSextic,
RealOpenMM** forceData) const;
private:
/**---------------------------------------------------------------------------------------
Calculate Amoeba harmonic angle ixn (force and energy)
Calculate Amoeba Out-Of-Plane-Bend ixn (force and energy)
@param positionAtomA Cartesian coordinates of atom A
@param positionAtomB Cartesian coordinates of atom B
......@@ -68,12 +105,12 @@ public:
--------------------------------------------------------------------------------------- */
static RealOpenMM calculateForceAndEnergy( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
const RealOpenMM* positionAtomC, const RealOpenMM* positionAtomD,
RealOpenMM angleK,
RealOpenMM angleCubic, RealOpenMM angleQuartic,
RealOpenMM anglePentic, RealOpenMM angleSextic,
RealOpenMM** forces );
RealOpenMM calculateOutOfPlaneBendIxn( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
const RealOpenMM* positionAtomC, const RealOpenMM* positionAtomD,
RealOpenMM angleK,
RealOpenMM angleCubic, RealOpenMM angleQuartic,
RealOpenMM anglePentic, RealOpenMM angleSextic,
RealOpenMM** forces ) const;
};
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferencePiTorsionForce.cpp
View file @ a37e41ea
......@@ -43,10 +43,10 @@
--------------------------------------------------------------------------------------- */
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 ){
RealOpenMM AmoebaReferencePiTorsionForce::calculatePiTorsionIxn( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
const RealOpenMM* positionAtomC, const RealOpenMM* positionAtomD,
const RealOpenMM* positionAtomE, const RealOpenMM* positionAtomF,
RealOpenMM piTorsionK, RealOpenMM** forces ) const {
// ---------------------------------------------------------------------------------------
......@@ -182,3 +182,37 @@ RealOpenMM AmoebaReferencePiTorsionForce::calculateForceAndEnergy( const RealOpe
return piTorsionK*phi2;
}
RealOpenMM AmoebaReferencePiTorsionForce::calculateForceAndEnergy( int numPiTorsions, RealOpenMM** posData,
const std::vector<int>& particle1,
const std::vector<int>& particle2,
const std::vector<int>& particle3,
const std::vector<int>& particle4,
const std::vector<int>& particle5,
const std::vector<int>& particle6,
const std::vector<RealOpenMM>& kTorsion,
RealOpenMM** forceData ) const {
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 += calculatePiTorsionIxn( posData[particle1Index], posData[particle2Index],
posData[particle3Index], posData[particle4Index],
posData[particle5Index], posData[particle6Index],
kTorsion[ii], forces );
}
return energy;
}
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferencePiTorsionForce.h
View file @ a37e41ea
......@@ -26,6 +26,7 @@
#define __AmoebaReferencePiTorsionForce_H__
#include "SimTKUtilities/SimTKOpenMMRealType.h"
#include <vector>
// ---------------------------------------------------------------------------------------
......@@ -39,7 +40,7 @@ public:
--------------------------------------------------------------------------------------- */
AmoebaReferencePiTorsionForce( );
AmoebaReferencePiTorsionForce( ){};
/**---------------------------------------------------------------------------------------
......@@ -47,11 +48,46 @@ public:
--------------------------------------------------------------------------------------- */
~AmoebaReferencePiTorsionForce( );
~AmoebaReferencePiTorsionForce( ){};
/**---------------------------------------------------------------------------------------
Calculate Amoeba torsion ixns (force and energy)
@param numTorsions number of torsions
@param posData particle positions
@param particle1 particle 1 indices
@param particle2 particle 2 indices
@param particle3 particle 3 indices
@param particle4 particle 4 indices
@param particle5 particle 5 indices
@param particle6 particle 6 indices
@param torsionParameters1 first index torsion parameters (amplitude, phase, fold)
@param torsionParameters2 second index torsion parameters (amplitude, phase, fold)
@param torsionParameters3 third index torsion parameters (amplitude, phase, fold)
@param forces output force vector
@return total energy
--------------------------------------------------------------------------------------- */
RealOpenMM calculateForceAndEnergy( int numPiTorsions, RealOpenMM** posData,
const std::vector<int>& particle1,
const std::vector<int>& particle2,
const std::vector<int>& particle3,
const std::vector<int>& particle4,
const std::vector<int>& particle5,
const std::vector<int>& particle6,
const std::vector<RealOpenMM>& kTorsion,
RealOpenMM** forceData ) const;
private:
/**---------------------------------------------------------------------------------------
Calculate Amoeba harmonic angle ixn (force and energy)
Calculate Amoeba pi-torsion ixn (force and energy)
@param positionAtomA Cartesian coordinates of atom A
@param positionAtomB Cartesian coordinates of atom B
......@@ -66,10 +102,10 @@ public:
--------------------------------------------------------------------------------------- */
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 );
RealOpenMM calculatePiTorsionIxn( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
const RealOpenMM* positionAtomC, const RealOpenMM* positionAtomD,
const RealOpenMM* positionAtomE, const RealOpenMM* positionAtomF,
RealOpenMM kTorsion, RealOpenMM** forces ) const;
};
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceStretchBendForce.cpp
View file @ a37e41ea
......@@ -28,7 +28,7 @@
/**---------------------------------------------------------------------------------------
Calculate Amoeba harmonic angle ixn (force and energy)
Calculate Amoeba stretch bend angle ixn (force and energy)
@param positionAtomA Cartesian coordinates of atom A
@param positionAtomB Cartesian coordinates of atom B
......@@ -46,11 +46,11 @@
--------------------------------------------------------------------------------------- */
RealOpenMM AmoebaReferenceStretchBendForce::calculateForceAndEnergy( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
RealOpenMM AmoebaReferenceStretchBendForce::calculateStretchBendIxn( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
const RealOpenMM* positionAtomC,
RealOpenMM lengthAB, RealOpenMM lengthCB,
RealOpenMM idealAngle, RealOpenMM kParameter,
RealOpenMM** forces ){
RealOpenMM** forces ) const {
// ---------------------------------------------------------------------------------------
......@@ -141,6 +141,34 @@ RealOpenMM AmoebaReferenceStretchBendForce::calculateForceAndEnergy( const RealO
// ---------------------------------------------------------------------------------------
// return energy
return (kParameter*dt*dr);
}
RealOpenMM AmoebaReferenceStretchBendForce::calculateForceAndEnergy( int numStretchBends, RealOpenMM** posData,
const std::vector<int>& particle1,
const std::vector<int>& particle2,
const std::vector<int>& particle3,
const std::vector<RealOpenMM>& lengthABParameters,
const std::vector<RealOpenMM>& lengthCBParameters,
const std::vector<RealOpenMM>& angle,
const std::vector<RealOpenMM>& kQuadratic,
RealOpenMM** forceData) const {
RealOpenMM energy = 0.0;
for (unsigned int ii = 0; ii < numStretchBends; ii++) {
int particle1Index = particle1[ii];
int particle2Index = particle2[ii];
int particle3Index = particle3[ii];
RealOpenMM abLength = lengthABParameters[ii];
RealOpenMM cbLength = lengthCBParameters[ii];
RealOpenMM idealAngle = angle[ii];
RealOpenMM angleK = kQuadratic[ii];
RealOpenMM* forces[3];
forces[0] = forceData[particle1Index];
forces[1] = forceData[particle2Index];
forces[2] = forceData[particle3Index];
energy += calculateStretchBendIxn( posData[particle1Index], posData[particle2Index], posData[particle3Index],
abLength, cbLength, idealAngle, angleK, forces );
}
return energy;
}
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceStretchBendForce.h
View file @ a37e41ea
......@@ -26,6 +26,7 @@
#define __AmoebaReferenceStretchBendForce_H__
#include "SimTKUtilities/SimTKOpenMMRealType.h"
#include <vector>
// ---------------------------------------------------------------------------------------
......@@ -39,7 +40,7 @@ public:
--------------------------------------------------------------------------------------- */
AmoebaReferenceStretchBendForce( );
AmoebaReferenceStretchBendForce( ){};
/**---------------------------------------------------------------------------------------
......@@ -47,11 +48,44 @@ public:
--------------------------------------------------------------------------------------- */
~AmoebaReferenceStretchBendForce( );
~AmoebaReferenceStretchBendForce( ){};
/**---------------------------------------------------------------------------------------
Calculate Amoeba stretch bend ixns (force and energy)
@param numBonds number of angles
@param posData particle positions
@param particle1 particle 1 indices
@param particle2 particle 2 indices
@param particle3 particle 3 indices
@param lengthABParameters ideal AB bond length
@param lengthCBParameters ideal CB bond length
@param angle ideal angle
@param kQuadratic force constant
@param forces output force vector
@return total energy
--------------------------------------------------------------------------------------- */
RealOpenMM calculateForceAndEnergy( int numAngles, RealOpenMM** posData,
const std::vector<int>& particle1,
const std::vector<int>& particle2,
const std::vector<int>& particle3,
const std::vector<RealOpenMM>& lengthABParameters,
const std::vector<RealOpenMM>& lengthCBParameters,
const std::vector<RealOpenMM>& angle,
const std::vector<RealOpenMM>& kQuadratic,
RealOpenMM** forceData ) const;
private:
/**---------------------------------------------------------------------------------------
Calculate Amoeba harmonic angle ixn (force and energy)
Calculate Amoeba stretch bend angle ixn (force and energy)
@param positionAtomA Cartesian coordinates of atom A
@param positionAtomB Cartesian coordinates of atom B
......@@ -65,11 +99,12 @@ public:
@return energy
--------------------------------------------------------------------------------------- */
static RealOpenMM calculateForceAndEnergy( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
const RealOpenMM* positionAtomC,
RealOpenMM lengthAB, RealOpenMM lengthCB,
RealOpenMM idealAngle, RealOpenMM kParameter,
RealOpenMM** forces );
RealOpenMM calculateStretchBendIxn( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
const RealOpenMM* positionAtomC,
RealOpenMM lengthAB, RealOpenMM lengthCB,
RealOpenMM idealAngle, RealOpenMM kParameter,
RealOpenMM** forces ) const;
};
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceTorsionForce.cpp
View file @ a37e41ea
/* Portions copyright (c) 2006 Stanford University and Simbios.
* Contributors: Pande Group
*
......@@ -24,7 +23,6 @@
#include "AmoebaReferenceForce.h"
#include "AmoebaReferenceTorsionForce.h"
#include <vector>
/**---------------------------------------------------------------------------------------
......@@ -43,12 +41,12 @@
--------------------------------------------------------------------------------------- */
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 ){
RealOpenMM AmoebaReferenceTorsionForce::calculateTorsionIxn( 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 ) const {
// ---------------------------------------------------------------------------------------
......@@ -210,3 +208,31 @@ RealOpenMM AmoebaReferenceTorsionForce::calculateForceAndEnergy( const RealOpenM
return energy;
}
RealOpenMM AmoebaReferenceTorsionForce::calculateForceAndEnergy( int numTorsions, RealOpenMM** posData,
const std::vector<int>& particle1,
const std::vector<int>& particle2,
const std::vector<int>& particle3,
const std::vector<int>& particle4,
const std::vector< std::vector<RealOpenMM> >& torsionParameters1,
const std::vector< std::vector<RealOpenMM> >& torsionParameters2,
const std::vector< std::vector<RealOpenMM> >& torsionParameters3,
RealOpenMM** forceData ) const {
RealOpenMM energy = 0.0;
for (unsigned int ii = 0; ii < numTorsions; ii++) {
int particle1Index = particle1[ii];
int particle2Index = particle2[ii];
int particle3Index = particle3[ii];
int particle4Index = particle4[ii];
RealOpenMM* forces[4];
forces[0] = forceData[particle1Index];
forces[1] = forceData[particle2Index];
forces[2] = forceData[particle3Index];
forces[3] = forceData[particle4Index];
energy += calculateTorsionIxn( posData[particle1Index], posData[particle2Index], posData[particle3Index], posData[particle4Index],
torsionParameters1[ii], torsionParameters2[ii], torsionParameters3[ii], forces );
}
return energy;
}
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceTorsionForce.h
View file @ a37e41ea
......@@ -26,6 +26,7 @@
#define __AmoebaReferenceTorsionForce_H__
#include "SimTKUtilities/SimTKOpenMMRealType.h"
#include <vector>
// ---------------------------------------------------------------------------------------
......@@ -39,7 +40,7 @@ public:
--------------------------------------------------------------------------------------- */
AmoebaReferenceTorsionForce( );
AmoebaReferenceTorsionForce( ){};
/**---------------------------------------------------------------------------------------
......@@ -47,8 +48,40 @@ public:
--------------------------------------------------------------------------------------- */
~AmoebaReferenceTorsionForce( );
~AmoebaReferenceTorsionForce( ){};
/**---------------------------------------------------------------------------------------
Calculate Amoeba torsion ixns (force and energy)
@param numTorsions number of torsions
@param posData particle positions
@param particle1 particle 1 indices
@param particle2 particle 2 indices
@param particle3 particle 3 indices
@param particle4 particle 4 indices
@param torsionParameters1 first index torsion parameters (amplitude, phase, fold)
@param torsionParameters2 second index torsion parameters (amplitude, phase, fold)
@param torsionParameters3 third index torsion parameters (amplitude, phase, fold)
@param forces output force vector
@return total energy
--------------------------------------------------------------------------------------- */
RealOpenMM calculateForceAndEnergy( int numTorsions, RealOpenMM** posData,
const std::vector<int>& particle1,
const std::vector<int>& particle2,
const std::vector<int>& particle3,
const std::vector<int>& particle4,
const std::vector< std::vector<RealOpenMM> >& torsionParameters1,
const std::vector< std::vector<RealOpenMM> >& torsionParameters2,
const std::vector< std::vector<RealOpenMM> >& torsionParameters3,
RealOpenMM** forceData ) const;
private:
/**---------------------------------------------------------------------------------------
Calculate Amoeba harmonic angle ixn (force and energy)
......@@ -66,15 +99,16 @@ public:
--------------------------------------------------------------------------------------- */
static RealOpenMM calculateForceAndEnergy( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
RealOpenMM calculateTorsionIxn( 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 );
RealOpenMM** forces ) const;
};
// ---------------------------------------------------------------------------------------
#endif // _AmoebaReferenceTorsionForce___
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceTorsionTorsionForce.cpp
View file @ a37e41ea
......@@ -24,11 +24,10 @@
#include "AmoebaReferenceForce.h"
#include "AmoebaReferenceTorsionTorsionForce.h"
#include <vector>
/**---------------------------------------------------------------------------------------
Load grid values from rectangle enclosing angles
Load grid values from rectenclosing angles
@param grid grid values [ angle1, angle2, f, df1, df2, df12 ]
@param angle1 angle in first dimension
......@@ -50,7 +49,7 @@
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 ){
RealOpenMM* fValues, RealOpenMM* fValues1, RealOpenMM* fValues2, RealOpenMM* fValues12 ) const {
// ---------------------------------------------------------------------------------------
......@@ -119,7 +118,7 @@ void AmoebaReferenceTorsionTorsionForce::loadGridValuesFromEnclosingRectangle(
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] ){
RealOpenMM c[4][4] ) const {
// ---------------------------------------------------------------------------------------
......@@ -221,7 +220,7 @@ void AmoebaReferenceTorsionTorsionForce::getBicubicValues(
const RealOpenMM x1Lower, const RealOpenMM x1Upper,
const RealOpenMM x2Lower, const RealOpenMM x2Upper,
const RealOpenMM gridValue1, const RealOpenMM gridValue2,
RealOpenMM* functionValue, RealOpenMM* functionValue1, RealOpenMM* functionValue2 ){
RealOpenMM* functionValue, RealOpenMM* functionValue1, RealOpenMM* functionValue2 ) const {
// ---------------------------------------------------------------------------------------
......@@ -275,7 +274,7 @@ void AmoebaReferenceTorsionTorsionForce::getBicubicValues(
int AmoebaReferenceTorsionTorsionForce::checkTorsionSign(
const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
const RealOpenMM* positionAtomC, const RealOpenMM* positionAtomD ){
const RealOpenMM* positionAtomC, const RealOpenMM* positionAtomD ) const {
// ---------------------------------------------------------------------------------------
......@@ -324,11 +323,11 @@ int AmoebaReferenceTorsionTorsionForce::checkTorsionSign(
--------------------------------------------------------------------------------------- */
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 ){
RealOpenMM AmoebaReferenceTorsionTorsionForce::calculateTorsionTorsionIxn( 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 ) const {
// ---------------------------------------------------------------------------------------
......@@ -549,3 +548,48 @@ RealOpenMM AmoebaReferenceTorsionTorsionForce::calculateForceAndEnergy( const Re
return gridEnergy;
}
RealOpenMM AmoebaReferenceTorsionTorsionForce::calculateForceAndEnergy( int numTorsionTorsions, RealOpenMM** posData,
const std::vector<int>& particle1,
const std::vector<int>& particle2,
const std::vector<int>& particle3,
const std::vector<int>& particle4,
const std::vector<int>& particle5,
const std::vector<int>& chiralCheckAtom,
const std::vector<int>& gridIndices,
const std::vector< std::vector< std::vector< std::vector<RealOpenMM> > > >& torsionTorsionGrids,
RealOpenMM** forceData ) const {
RealOpenMM energy = 0.0;
for (unsigned int ii = 0; ii < numTorsionTorsions; ii++) {
int particle1Index = particle1[ii];
int particle2Index = particle2[ii];
int particle3Index = particle3[ii];
int particle4Index = particle4[ii];
int particle5Index = particle5[ii];
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 += calculateTorsionTorsionIxn( posData[particle1Index], posData[particle2Index],
posData[particle3Index], posData[particle4Index],
posData[particle5Index], chiralCheckAtom, torsionTorsionGrids[gridIndex],
forces );
}
return energy;
}
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceTorsionTorsionForce.h
View file @ a37e41ea
......@@ -26,6 +26,7 @@
#define __AmoebaReferenceTorsionTorsionForce_H__
#include "SimTKUtilities/SimTKOpenMMRealType.h"
#include <vector>
// ---------------------------------------------------------------------------------------
......@@ -39,7 +40,7 @@ public:
--------------------------------------------------------------------------------------- */
AmoebaReferenceTorsionTorsionForce( );
AmoebaReferenceTorsionTorsionForce( ){};
/**---------------------------------------------------------------------------------------
......@@ -47,7 +48,41 @@ public:
--------------------------------------------------------------------------------------- */
~AmoebaReferenceTorsionTorsionForce( );
~AmoebaReferenceTorsionTorsionForce( ){};
/**---------------------------------------------------------------------------------------
Calculate Amoeba torsion-torsion ixns (force and energy)
@param numTorsions number of torsions
@param posData particle positions
@param particle1 particle 1 indices
@param particle2 particle 2 indices
@param particle3 particle 3 indices
@param particle4 particle 4 indices
@param particle5 particle 5 indices
@param chiralCheckAtom chiralCheckAtom (-1 if no chiral check is to be performed)
@param gridIndices grid indices
@param torsionTorsionGrids torsion-torsion grids
@param forces output force vector
@return total energy
--------------------------------------------------------------------------------------- */
RealOpenMM calculateForceAndEnergy( int numTorsionTorsions, RealOpenMM** posData,
const std::vector<int>& particle1,
const std::vector<int>& particle2,
const std::vector<int>& particle3,
const std::vector<int>& particle4,
const std::vector<int>& particle5,
const std::vector<int>& chiralCheckAtom,
const std::vector<int>& gridIndices,
const std::vector< std::vector< std::vector< std::vector<RealOpenMM> > > >& torsionTorsionGrids,
RealOpenMM** forceData ) const;
private:
/**---------------------------------------------------------------------------------------
......@@ -67,10 +102,10 @@ public:
--------------------------------------------------------------------------------------- */
static void loadGridValuesFromEnclosingRectangle(
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 );
RealOpenMM* fValues, RealOpenMM* fValues1, RealOpenMM* fValues2, RealOpenMM* fValues12 ) const;
/**---------------------------------------------------------------------------------------
......@@ -93,9 +128,8 @@ public:
--------------------------------------------------------------------------------------- */
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] );
void getBicubicCoefficientMatrix( const RealOpenMM* y, const RealOpenMM* y1, const RealOpenMM* y2, const RealOpenMM* y12,
const RealOpenMM d1, const RealOpenMM d2, RealOpenMM c[4][4] ) const;
/**---------------------------------------------------------------------------------------
......@@ -130,12 +164,12 @@ public:
--------------------------------------------------------------------------------------- */
static void getBicubicValues(
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 );
RealOpenMM* functionValue, RealOpenMM* functionValue1, RealOpenMM* functionValue2 ) const;
/**---------------------------------------------------------------------------------------
......@@ -151,9 +185,8 @@ public:
--------------------------------------------------------------------------------------- */
static int checkTorsionSign(
const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
const RealOpenMM* positionAtomC, const RealOpenMM* positionAtomD );
int checkTorsionSign( const RealOpenMM* positionAtomA, const RealOpenMM* positionAtomB,
const RealOpenMM* positionAtomC, const RealOpenMM* positionAtomD ) const;
/**---------------------------------------------------------------------------------------
......@@ -173,11 +206,11 @@ public:
--------------------------------------------------------------------------------------- */
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 );
RealOpenMM calculateTorsionTorsionIxn( 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 ) const;
};
......
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