Skip to content

GitLab

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

Code to map torques onto forces

parent 50a15fb0
Hide whitespace changes
Inline Side-by-side
Showing with 177 additions and 9 deletions
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceForce.cpp
View file @ a33839b0
......@@ -133,7 +133,7 @@ RealOpenMM AmoebaReferenceForce::getNorm3( const RealOpenMM* inputVector ){
return SQRT( inputVector[0]*inputVector[0] + inputVector[1]*inputVector[1] + inputVector[2]*inputVector[2] );
}
void AmoebaReferenceForce::normalizeVector3( RealOpenMM* inputVector ){
RealOpenMM AmoebaReferenceForce::normalizeVector3( RealOpenMM* inputVector ){
// ---------------------------------------------------------------------------------------
......@@ -143,11 +143,13 @@ void AmoebaReferenceForce::normalizeVector3( RealOpenMM* inputVector ){
RealOpenMM norm = SQRT( inputVector[0]*inputVector[0] + inputVector[1]*inputVector[1] + inputVector[2]*inputVector[2] );
if( norm > 0.0 ){
norm = 1.0/norm;
inputVector[0] *= norm;
inputVector[1] *= norm;
inputVector[2] *= norm;
RealOpenMM normI = 1.0/norm;
inputVector[0] *= normI;
inputVector[1] *= normI;
inputVector[2] *= normI;
}
return norm;
}
/**---------------------------------------------------------------------------------------
......@@ -198,6 +200,19 @@ RealOpenMM AmoebaReferenceForce::getDotProduct3( const RealOpenMM* xVector, cons
return xVector[0]*yVector[0] + xVector[1]*yVector[1] + xVector[2]*yVector[2];
}
RealOpenMM AmoebaReferenceForce::getDotProduct3( const RealOpenMM* xVector, const OpenMM::Vec3& yVector ){
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "AmoebaReferenceForce::getDotProduct3";
// ---------------------------------------------------------------------------------------
// get dot product
return xVector[0]*yVector[0] + xVector[1]*yVector[1] + xVector[2]*yVector[2];
}
/**---------------------------------------------------------------------------------------
Calculate dot product of 3d vectors
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceForce.h
View file @ a33839b0
......@@ -26,6 +26,7 @@
#define __AmoebaReferenceForce_H__
#include "SimTKUtilities/SimTKOpenMMRealType.h"
#include "openmm/Vec3.h"
#include <vector>
// ---------------------------------------------------------------------------------------
......@@ -95,10 +96,12 @@ public:
Normalize 3d vector
@param inputVector vector to normalize
@return norm
--------------------------------------------------------------------------------------- */
static void normalizeVector3( RealOpenMM* inputVector );
static RealOpenMM normalizeVector3( RealOpenMM* inputVector );
/**---------------------------------------------------------------------------------------
......@@ -113,6 +116,7 @@ public:
static RealOpenMM getDotProduct3( const std::vector<RealOpenMM>& xVector, const std::vector<RealOpenMM>& yVector );
static RealOpenMM getDotProduct3( const RealOpenMM* xVector, const RealOpenMM* yVector );
static RealOpenMM getDotProduct3( const RealOpenMM* xVector, const OpenMM::Vec3& yVector );
static RealOpenMM getDotProduct3( unsigned int vectorOffset, const std::vector<RealOpenMM>& xVector, const RealOpenMM* yVector );
/**---------------------------------------------------------------------------------------
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceMultipoleForce.cpp
View file @ a33839b0
......@@ -1342,8 +1342,110 @@ RealOpenMM AmoebaReferenceMultipoleForce::calculateNoCutoffElectrostaticPairIxn(
return energy;
}
RealOpenMM AmoebaReferenceMultipoleForce::calculateNoCutoffElectrostatic( std::vector<MultipoleParticleData>& particleData, RealOpenMM** forces ) const {
/**---------------------------------------------------------------------------------------
Map torques to forces
--------------------------------------------------------------------------------------- */
void AmoebaReferenceMultipoleForce::mapTorqueToForce( const MultipoleParticleData& particleI,
const MultipoleParticleData& particleU,
const MultipoleParticleData& particleV,
int axisType, const Vec3& torque, RealOpenMM** forces ) const {
// ---------------------------------------------------------------------------------------
static const std::string methodName = "mapTorqueToForce";
static const int U = 0;
static const int V = 1;
static const int W = 2;
static const RealOpenMM one = 1.0;
static const RealOpenMM two = 2.0;
RealOpenMM* vector[3];
// ---------------------------------------------------------------------------------------
// get coordinates of this atom and the z & x axis atoms
// compute the vector between the atoms and 1/sqrt(d2), d2 is distance between
// this atom and the axis atom
RealOpenMM workspace[9];
vector[U] = workspace;
vector[V] = workspace + 3;
vector[W] = workspace + 6;
RealOpenMM norms[3];
getDelta( particleI, particleU, vector[U] );
norms[U] = AmoebaReferenceForce::normalizeVector3( vector[U] );
getDelta( particleI, particleV, vector[V] );
norms[V] = AmoebaReferenceForce::normalizeVector3( vector[V] );
AmoebaReferenceForce::getCrossProduct( vector[U], vector[V], vector[W] );
AmoebaReferenceForce::normalizeVector3( vector[W] );
RealOpenMM diff[3];
diff[0] = vector[V][0] - vector[U][0];
diff[1] = vector[V][1] - vector[U][1];
diff[2] = vector[V][2] - vector[U][2];
RealOpenMM dotDu = AmoebaReferenceForce::getDotProduct3( vector[U], diff );
RealOpenMM dotDv = AmoebaReferenceForce::getDotProduct3( vector[V], diff );
RealOpenMM up[3], vp[3];
for( int jj = 0; jj < 3; jj++ ){
up[jj] = diff[jj] - dotDu*vector[U][jj];
vp[jj] = diff[jj] - dotDv*vector[V][jj];
}
AmoebaReferenceForce::normalizeVector3( up );
AmoebaReferenceForce::normalizeVector3( vp );
RealOpenMM dphi[3];
for( int ii = 0; ii < 3; ii++ ){
dphi[ii] = AmoebaReferenceForce::getDotProduct3( vector[ii], torque );
}
RealOpenMM c = AmoebaReferenceForce::getDotProduct3( vector[U], vector[V] );
RealOpenMM s = SQRT( one - (c*c) );
RealOpenMM uvdis = norms[U]*s;
RealOpenMM vudis = norms[V]*s;
// branch based on axis type
if( axisType == AmoebaMultipoleForce::Bisector ){
// bisector
for( int ii = 0; ii < 3; ii++ ){
RealOpenMM du = -vector[W][ii]*dphi[V]/uvdis + up[ii]*dphi[W]/(two*norms[U]);
RealOpenMM dv = vector[W][ii]*dphi[U]/vudis + vp[ii]*dphi[W]/(two*norms[V]);
forces[particleU.particleIndex][ii] += du;
forces[particleV.particleIndex][ii] += dv;
forces[particleI.particleIndex][ii] -= (dv + du);
}
} else if( axisType == AmoebaMultipoleForce::ZThenX ){
for( int ii = 0; ii < 3; ii++ ){
RealOpenMM du = -vector[W][ii]*dphi[V]/uvdis + up[ii]*dphi[W]/norms[U];
RealOpenMM dv = vector[W][ii]*dphi[U]/vudis;
forces[particleU.particleIndex][ii] += du;
forces[particleV.particleIndex][ii] += dv;
forces[particleI.particleIndex][ii] -= (dv + du);
}
}
return;
}
RealOpenMM AmoebaReferenceMultipoleForce::calculateNoCutoffElectrostatic( std::vector<MultipoleParticleData>& particleData,
const std::vector<int>& axisTypes,
const std::vector<int>& multipoleAtomId1s,
const std::vector<int>& multipoleAtomId2s,
RealOpenMM** forces ) const {
// ---------------------------------------------------------------------------------------
......@@ -1400,6 +1502,8 @@ RealOpenMM AmoebaReferenceMultipoleForce::calculateNoCutoffElectrostatic( std::v
}
}
// diagnostics
if( log ){
RealOpenMM conversion = 1.0/4.184;
RealOpenMM torqueConversion = conversion;
......@@ -1416,6 +1520,29 @@ RealOpenMM AmoebaReferenceMultipoleForce::calculateNoCutoffElectrostatic( std::v
}
}
// map torques to forces
for( unsigned int ii = 0; ii < particleData.size(); ii++ ){
mapTorqueToForce( particleData[ii], particleData[multipoleAtomId1s[ii]], particleData[multipoleAtomId2s[ii]], axisTypes[ii], torques[ii], forces );
}
// diagnostics
if( log ){
RealOpenMM conversion = 1.0/4.184;
RealOpenMM torqueConversion = conversion;
(void) fprintf( log, "Force post torque mapping energy=%15.7e\n", 10.0*energy*conversion );
conversion *= -0.1;
unsigned int maxPrint = 10000;
for( unsigned int ii = 0; ii < particleData.size(); ii++ ){
(void) fprintf( log, "%6u [%15.7e %15.7e %15.7e]\n", ii,
conversion*forces[ii][0], conversion*forces[ii][1], conversion*forces[ii][2] );
if( ii == maxPrint && 2*maxPrint < particleData.size() ){
ii = particleData.size() - maxPrint;
}
}
}
return energy;
}
......@@ -1510,7 +1637,8 @@ RealOpenMM AmoebaReferenceMultipoleForce::calculateNoCutoffForceAndEnergy( unsig
throw OpenMMException(message.str());
}
RealOpenMM energy = calculateNoCutoffElectrostatic( particleData, forces );
RealOpenMM energy = calculateNoCutoffElectrostatic( particleData, axisTypes, multipoleAtomId1s, multipoleAtomId2s, forces );
return energy;
}
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceMultipoleForce.h
View file @ a33839b0
......@@ -442,7 +442,28 @@ private:
const MultipoleParticleData& particleK,
RealOpenMM* scalingFactors, RealOpenMM** forces, std::vector<Vec3>& torque ) const;
RealOpenMM calculateNoCutoffElectrostatic( std::vector<MultipoleParticleData>& particleData, RealOpenMM** forces ) const;
/**---------------------------------------------------------------------------------------
Map torque to forces
@param particleI particle whose torque is to be mapped
@param particleU particle1 of lab frame for particleI
@param particleV particle2 of lab frame for particleI
@param axisType axis type (Bisector/Z-then-X, ...)
@param torque torque on particle I
--------------------------------------------------------------------------------------- */
void mapTorqueToForce( const MultipoleParticleData& particleI,
const MultipoleParticleData& particleU,
const MultipoleParticleData& particleV,
int axisType, const Vec3& torque, RealOpenMM** forces ) const;
RealOpenMM calculateNoCutoffElectrostatic( std::vector<MultipoleParticleData>& particleData,
const std::vector<int>& axisTypes,
const std::vector<int>& multipoleAtomId1s,
const std::vector<int>& multipoleAtomId2s,
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