Skip to content

GitLab

Commit 54e47d2d authored by peastman's avatar peastman
Browse files

Converted lots of ExpressionPrograms to CompiledExpressions

parent 8b3bfc92
Hide whitespace changes
Inline Side-by-side
Showing with 285 additions and 170 deletions
openmmapi/include/openmm/internal/CustomNonbondedForceImpl.h
View file @ 54e47d2d
...@@ -9,7 +9,7 @@ ...@@ -9,7 +9,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for * * Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. * * Medical Research, grant U54 GM072970. See https://simtk.org. *
* * * *
* Portions copyright (c) 2008-2012 Stanford University and the Authors. * * Portions copyright (c) 2008-2013 Stanford University and the Authors. *
* Authors: Peter Eastman * * Authors: Peter Eastman *
* Contributors: * * Contributors: *
* * * *
...@@ -35,7 +35,7 @@ ...@@ -35,7 +35,7 @@
#include "ForceImpl.h" #include "ForceImpl.h"
#include "openmm/CustomNonbondedForce.h" #include "openmm/CustomNonbondedForce.h"
#include "openmm/Kernel.h" #include "openmm/Kernel.h"
#include "lepton/ExpressionProgram.h" #include "lepton/CompiledExpression.h"
#include <utility> #include <utility>
#include <map> #include <map>
#include <string> #include <string>
...@@ -68,7 +68,7 @@ public: ...@@ -68,7 +68,7 @@ public:
static double calcLongRangeCorrection(const CustomNonbondedForce& force, const Context& context); static double calcLongRangeCorrection(const CustomNonbondedForce& force, const Context& context);
private: private:
class TabulatedFunction; class TabulatedFunction;
static double integrateInteraction(const Lepton::ExpressionProgram& expression, const std::vector<double>& params1, const std::vector<double>& params2, static double integrateInteraction(Lepton::CompiledExpression& expression, const std::vector<double>& params1, const std::vector<double>& params2,
const CustomNonbondedForce& force, const Context& context); const CustomNonbondedForce& force, const Context& context);
const CustomNonbondedForce& owner; const CustomNonbondedForce& owner;
Kernel kernel; Kernel kernel;
......
openmmapi/src/CustomNonbondedForceImpl.cpp
View file @ 54e47d2d
...@@ -182,7 +182,7 @@ double CustomNonbondedForceImpl::calcLongRangeCorrection(const CustomNonbondedFo ...@@ -182,7 +182,7 @@ double CustomNonbondedForceImpl::calcLongRangeCorrection(const CustomNonbondedFo
force.getFunctionParameters(i, name, values, min, max); force.getFunctionParameters(i, name, values, min, max);
functions[name] = new TabulatedFunction(min, max, values); functions[name] = new TabulatedFunction(min, max, values);
} }
Lepton::ExpressionProgram expression = Lepton::Parser::parse(force.getEnergyFunction(), functions).optimize().createProgram(); Lepton::CompiledExpression expression = Lepton::Parser::parse(force.getEnergyFunction(), functions).createCompiledExpression();
// Identify all particle classes (defined by parameters), and record the class of each particle. // Identify all particle classes (defined by parameters), and record the class of each particle.
...@@ -251,25 +251,35 @@ double CustomNonbondedForceImpl::calcLongRangeCorrection(const CustomNonbondedFo ...@@ -251,25 +251,35 @@ double CustomNonbondedForceImpl::calcLongRangeCorrection(const CustomNonbondedFo
return 2*M_PI*numParticles*numParticles*sum; return 2*M_PI*numParticles*numParticles*sum;
} }
double CustomNonbondedForceImpl::integrateInteraction(const Lepton::ExpressionProgram& expression, const vector<double>& params1, const vector<double>& params2, double CustomNonbondedForceImpl::integrateInteraction(Lepton::CompiledExpression& expression, const vector<double>& params1, const vector<double>& params2,
const CustomNonbondedForce& force, const Context& context) { const CustomNonbondedForce& force, const Context& context) {
map<std::string, double> variables; const set<string>& variables = expression.getVariables();
for (int i = 0; i < force.getNumPerParticleParameters(); i++) { for (int i = 0; i < force.getNumPerParticleParameters(); i++) {
stringstream name1, name2; stringstream name1, name2;
name1 << force.getPerParticleParameterName(i) << 1; name1 << force.getPerParticleParameterName(i) << 1;
name2 << force.getPerParticleParameterName(i) << 2; name2 << force.getPerParticleParameterName(i) << 2;
variables[name1.str()] = params1[i]; if (variables.find(name1.str()) != variables.end())
variables[name2.str()] = params2[i]; expression.getVariableReference(name1.str()) = params1[i];
if (variables.find(name2.str()) != variables.end())
expression.getVariableReference(name2.str()) = params2[i];
} }
for (int i = 0; i < force.getNumGlobalParameters(); i++) { for (int i = 0; i < force.getNumGlobalParameters(); i++) {
const string& name = force.getGlobalParameterName(i); const string& name = force.getGlobalParameterName(i);
variables[name] = context.getParameter(name); if (variables.find(name) != variables.end())
expression.getVariableReference(name) = context.getParameter(name);
} }
// To integrate from r_cutoff to infinity, make the change of variables x=r_cutoff/r and integrate from 0 to 1. // To integrate from r_cutoff to infinity, make the change of variables x=r_cutoff/r and integrate from 0 to 1.
// This introduces another r^2 into the integral, which along with the r^2 in the formula for the correction // This introduces another r^2 into the integral, which along with the r^2 in the formula for the correction
// means we multiply the function by r^4. Use the midpoint method. // means we multiply the function by r^4. Use the midpoint method.
double* rPointer;
try {
rPointer = &expression.getVariableReference("r");
}
catch (exception& ex) {
throw OpenMMException("CustomNonbondedForce: Cannot use long range correction with a force that does not depend on r.");
}
double cutoff = force.getCutoffDistance(); double cutoff = force.getCutoffDistance();
double sum = 0; double sum = 0;
int numPoints = 1; int numPoints = 1;
...@@ -281,9 +291,9 @@ double CustomNonbondedForceImpl::integrateInteraction(const Lepton::ExpressionPr ...@@ -281,9 +291,9 @@ double CustomNonbondedForceImpl::integrateInteraction(const Lepton::ExpressionPr
continue; continue;
double x = (i+0.5)/numPoints; double x = (i+0.5)/numPoints;
double r = cutoff/x; double r = cutoff/x;
variables["r"] = r; *rPointer = r;
double r2 = r*r; double r2 = r*r;
newSum += expression.evaluate(variables)*r2*r2; newSum += expression.evaluate()*r2*r2;
} }
sum = newSum/numPoints + oldSum/3; sum = newSum/numPoints + oldSum/3;
if (iteration > 2 && (fabs((sum-oldSum)/sum) < 1e-5 || sum == 0)) if (iteration > 2 && (fabs((sum-oldSum)/sum) < 1e-5 || sum == 0))
...@@ -309,8 +319,8 @@ double CustomNonbondedForceImpl::integrateInteraction(const Lepton::ExpressionPr ...@@ -309,8 +319,8 @@ double CustomNonbondedForceImpl::integrateInteraction(const Lepton::ExpressionPr
double x = (i+0.5)/numPoints; double x = (i+0.5)/numPoints;
double r = rswitch+x*(cutoff-rswitch); double r = rswitch+x*(cutoff-rswitch);
double switchValue = x*x*x*(10+x*(-15+x*6)); double switchValue = x*x*x*(10+x*(-15+x*6));
variables["r"] = r; *rPointer = r;
newSum += switchValue*expression.evaluate(variables)*r*r; newSum += switchValue*expression.evaluate()*r*r;
} }
sum2 = newSum/numPoints + oldSum/3; sum2 = newSum/numPoints + oldSum/3;
if (iteration > 2 && (fabs((sum2-oldSum)/sum2) < 1e-5 || sum2 == 0)) if (iteration > 2 && (fabs((sum2-oldSum)/sum2) < 1e-5 || sum2 == 0))
......
platforms/reference/include/ReferenceCustomAngleIxn.h
View file @ 54e47d2d
/* Portions copyright (c) 2010 Stanford University and Simbios. /* Portions copyright (c) 2010-2013 Stanford University and Simbios.
* Contributors: Peter Eastman * Contributors: Peter Eastman
* *
* Permission is hereby granted, free of charge, to any person obtaining * Permission is hereby granted, free of charge, to any person obtaining
...@@ -25,17 +25,20 @@ ...@@ -25,17 +25,20 @@
#define __ReferenceCustomAngleIxn_H__ #define __ReferenceCustomAngleIxn_H__
#include "ReferenceBondIxn.h" #include "ReferenceBondIxn.h"
#include "lepton/ExpressionProgram.h" #include "lepton/CompiledExpression.h"
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
class ReferenceCustomAngleIxn : public ReferenceBondIxn { class ReferenceCustomAngleIxn : public ReferenceBondIxn {
private: private:
Lepton::ExpressionProgram energyExpression; Lepton::CompiledExpression energyExpression;
Lepton::ExpressionProgram forceExpression; Lepton::CompiledExpression forceExpression;
std::vector<std::string> paramNames; std::vector<double*> energyParams;
std::map<std::string, double> globalParameters; std::vector<double*> forceParams;
double* energyTheta;
double* forceTheta;
int numParameters;
public: public:
...@@ -45,7 +48,7 @@ class ReferenceCustomAngleIxn : public ReferenceBondIxn { ...@@ -45,7 +48,7 @@ class ReferenceCustomAngleIxn : public ReferenceBondIxn {
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
ReferenceCustomAngleIxn(const Lepton::ExpressionProgram& energyExpression, const Lepton::ExpressionProgram& forceExpression, ReferenceCustomAngleIxn(const Lepton::CompiledExpression& energyExpression, const Lepton::CompiledExpression& forceExpression,
const std::vector<std::string>& parameterNames, std::map<std::string, double> globalParameters); const std::vector<std::string>& parameterNames, std::map<std::string, double> globalParameters);
/**--------------------------------------------------------------------------------------- /**---------------------------------------------------------------------------------------
......
platforms/reference/include/ReferenceCustomBondIxn.h
View file @ 54e47d2d
/* Portions copyright (c) 2009 Stanford University and Simbios. /* Portions copyright (c) 2009-2013 Stanford University and Simbios.
* Contributors: Peter Eastman * Contributors: Peter Eastman
* *
* Permission is hereby granted, free of charge, to any person obtaining * Permission is hereby granted, free of charge, to any person obtaining
...@@ -26,17 +26,20 @@ ...@@ -26,17 +26,20 @@
#define __ReferenceCustomBondIxn_H__ #define __ReferenceCustomBondIxn_H__
#include "ReferenceBondIxn.h" #include "ReferenceBondIxn.h"
#include "lepton/ExpressionProgram.h" #include "lepton/CompiledExpression.h"
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
class ReferenceCustomBondIxn : public ReferenceBondIxn { class ReferenceCustomBondIxn : public ReferenceBondIxn {
private: private:
Lepton::ExpressionProgram energyExpression; Lepton::CompiledExpression energyExpression;
Lepton::ExpressionProgram forceExpression; Lepton::CompiledExpression forceExpression;
std::vector<std::string> paramNames; std::vector<double*> energyParams;
std::map<std::string, double> globalParameters; std::vector<double*> forceParams;
double* energyR;
double* forceR;
int numParameters;
public: public:
...@@ -46,7 +49,7 @@ class ReferenceCustomBondIxn : public ReferenceBondIxn { ...@@ -46,7 +49,7 @@ class ReferenceCustomBondIxn : public ReferenceBondIxn {
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
ReferenceCustomBondIxn(const Lepton::ExpressionProgram& energyExpression, const Lepton::ExpressionProgram& forceExpression, ReferenceCustomBondIxn(const Lepton::CompiledExpression& energyExpression, const Lepton::CompiledExpression& forceExpression,
const std::vector<std::string>& parameterNames, std::map<std::string, double> globalParameters); const std::vector<std::string>& parameterNames, std::map<std::string, double> globalParameters);
/**--------------------------------------------------------------------------------------- /**---------------------------------------------------------------------------------------
......
platforms/reference/include/ReferenceCustomExternalIxn.h
View file @ 54e47d2d
/* Portions copyright (c) 2009 Stanford University and Simbios. /* Portions copyright (c) 2009-2013 Stanford University and Simbios.
* Contributors: Peter Eastman * Contributors: Peter Eastman
* *
* Permission is hereby granted, free of charge, to any person obtaining * Permission is hereby granted, free of charge, to any person obtaining
...@@ -26,19 +26,26 @@ ...@@ -26,19 +26,26 @@
#define __ReferenceCustomExternalIxn_H__ #define __ReferenceCustomExternalIxn_H__
#include "ReferenceCustomExternalIxn.h" #include "ReferenceCustomExternalIxn.h"
#include "lepton/ExpressionProgram.h" #include "lepton/CompiledExpression.h"
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
class ReferenceCustomExternalIxn { class ReferenceCustomExternalIxn {
private: private:
Lepton::ExpressionProgram energyExpression; Lepton::CompiledExpression energyExpression;
Lepton::ExpressionProgram forceExpressionX; Lepton::CompiledExpression forceExpressionX;
Lepton::ExpressionProgram forceExpressionY; Lepton::CompiledExpression forceExpressionY;
Lepton::ExpressionProgram forceExpressionZ; Lepton::CompiledExpression forceExpressionZ;
std::vector<std::string> paramNames; std::vector<double*> energyParams;
std::map<std::string, double> globalParameters; std::vector<double*> forceXParams;
std::vector<double*> forceYParams;
std::vector<double*> forceZParams;
double *energyX, *energyY, *energyZ;
double *forceXX, *forceXY, *forceXZ;
double *forceYX, *forceYY, *forceYZ;
double *forceZX, *forceZY, *forceZZ;
int numParameters;
public: public:
...@@ -48,8 +55,8 @@ class ReferenceCustomExternalIxn { ...@@ -48,8 +55,8 @@ class ReferenceCustomExternalIxn {
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
ReferenceCustomExternalIxn(const Lepton::ExpressionProgram& energyExpression, const Lepton::ExpressionProgram& forceExpressionX, ReferenceCustomExternalIxn(const Lepton::CompiledExpression& energyExpression, const Lepton::CompiledExpression& forceExpressionX,
const Lepton::ExpressionProgram& forceExpressionY, const Lepton::ExpressionProgram& forceExpressionZ, const Lepton::CompiledExpression& forceExpressionY, const Lepton::CompiledExpression& forceExpressionZ,
const std::vector<std::string>& parameterNames, std::map<std::string, double> globalParameters); const std::vector<std::string>& parameterNames, std::map<std::string, double> globalParameters);
/**--------------------------------------------------------------------------------------- /**---------------------------------------------------------------------------------------
......
platforms/reference/include/ReferenceCustomNonbondedIxn.h
View file @ 54e47d2d
/* Portions copyright (c) 2009 Stanford University and Simbios. /* Portions copyright (c) 2009-2013 Stanford University and Simbios.
* Contributors: Peter Eastman * Contributors: Peter Eastman
* *
* Permission is hereby granted, free of charge, to any person obtaining * Permission is hereby granted, free of charge, to any person obtaining
...@@ -27,7 +27,7 @@ ...@@ -27,7 +27,7 @@
#include "ReferencePairIxn.h" #include "ReferencePairIxn.h"
#include "ReferenceNeighborList.h" #include "ReferenceNeighborList.h"
#include "lepton/ExpressionProgram.h" #include "lepton/CompiledExpression.h"
#include <map> #include <map>
#include <set> #include <set>
#include <utility> #include <utility>
...@@ -45,10 +45,13 @@ class ReferenceCustomNonbondedIxn { ...@@ -45,10 +45,13 @@ class ReferenceCustomNonbondedIxn {
const OpenMM::NeighborList* neighborList; const OpenMM::NeighborList* neighborList;
RealOpenMM periodicBoxSize[3]; RealOpenMM periodicBoxSize[3];
RealOpenMM cutoffDistance, switchingDistance; RealOpenMM cutoffDistance, switchingDistance;
Lepton::ExpressionProgram energyExpression; Lepton::CompiledExpression energyExpression;
Lepton::ExpressionProgram forceExpression; Lepton::CompiledExpression forceExpression;
std::vector<std::string> paramNames; std::vector<std::string> paramNames;
std::vector<std::string> particleParamNames; std::vector<double*> energyParticleParams;
std::vector<double*> forceParticleParams;
double* energyR;
double* forceR;
std::vector<std::pair<std::set<int>, std::set<int> > > interactionGroups; std::vector<std::pair<std::set<int>, std::set<int> > > interactionGroups;
/**--------------------------------------------------------------------------------------- /**---------------------------------------------------------------------------------------
...@@ -65,9 +68,8 @@ class ReferenceCustomNonbondedIxn { ...@@ -65,9 +68,8 @@ class ReferenceCustomNonbondedIxn {
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
void calculateOneIxn( int atom1, int atom2, std::vector<OpenMM::RealVec>& atomCoordinates, void calculateOneIxn( int atom1, int atom2, std::vector<OpenMM::RealVec>& atomCoordinates, std::vector<OpenMM::RealVec>& forces,
std::map<std::string, double>& variables, std::vector<OpenMM::RealVec>& forces, RealOpenMM* energyByAtom, RealOpenMM* totalEnergy );
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const;
public: public:
...@@ -78,7 +80,7 @@ class ReferenceCustomNonbondedIxn { ...@@ -78,7 +80,7 @@ class ReferenceCustomNonbondedIxn {
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
ReferenceCustomNonbondedIxn(const Lepton::ExpressionProgram& energyExpression, const Lepton::ExpressionProgram& forceExpression, ReferenceCustomNonbondedIxn(const Lepton::CompiledExpression& energyExpression, const Lepton::CompiledExpression& forceExpression,
const std::vector<std::string>& parameterNames); const std::vector<std::string>& parameterNames);
/**--------------------------------------------------------------------------------------- /**---------------------------------------------------------------------------------------
...@@ -153,7 +155,7 @@ class ReferenceCustomNonbondedIxn { ...@@ -153,7 +155,7 @@ class ReferenceCustomNonbondedIxn {
void calculatePairIxn( int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates, void calculatePairIxn( int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
RealOpenMM** atomParameters, std::vector<std::set<int> >& exclusions, RealOpenMM** atomParameters, std::vector<std::set<int> >& exclusions,
RealOpenMM* fixedParameters, const std::map<std::string, double>& globalParameters, RealOpenMM* fixedParameters, const std::map<std::string, double>& globalParameters,
std::vector<OpenMM::RealVec>& forces, RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const; std::vector<OpenMM::RealVec>& forces, RealOpenMM* energyByAtom, RealOpenMM* totalEnergy );
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
......
platforms/reference/include/ReferenceCustomTorsionIxn.h
View file @ 54e47d2d
/* Portions copyright (c) 2010 Stanford University and Simbios. /* Portions copyright (c) 2010-2013 Stanford University and Simbios.
* Contributors: Peter Eastman * Contributors: Peter Eastman
* *
* Permission is hereby granted, free of charge, to any person obtaining * Permission is hereby granted, free of charge, to any person obtaining
...@@ -25,17 +25,20 @@ ...@@ -25,17 +25,20 @@
#define __ReferenceCustomTorsionIxn_H__ #define __ReferenceCustomTorsionIxn_H__
#include "ReferenceBondIxn.h" #include "ReferenceBondIxn.h"
#include "lepton/ExpressionProgram.h" #include "lepton/CompiledExpression.h"
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
class ReferenceCustomTorsionIxn : public ReferenceBondIxn { class ReferenceCustomTorsionIxn : public ReferenceBondIxn {
private: private:
Lepton::ExpressionProgram energyExpression; Lepton::CompiledExpression energyExpression;
Lepton::ExpressionProgram forceExpression; Lepton::CompiledExpression forceExpression;
std::vector<std::string> paramNames; std::vector<double*> energyParams;
std::map<std::string, double> globalParameters; std::vector<double*> forceParams;
double* energyTheta;
double* forceTheta;
int numParameters;
public: public:
...@@ -45,7 +48,7 @@ class ReferenceCustomTorsionIxn : public ReferenceBondIxn { ...@@ -45,7 +48,7 @@ class ReferenceCustomTorsionIxn : public ReferenceBondIxn {
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
ReferenceCustomTorsionIxn(const Lepton::ExpressionProgram& energyExpression, const Lepton::ExpressionProgram& forceExpression, ReferenceCustomTorsionIxn(const Lepton::CompiledExpression& energyExpression, const Lepton::CompiledExpression& forceExpression,
const std::vector<std::string>& parameterNames, std::map<std::string, double> globalParameters); const std::vector<std::string>& parameterNames, std::map<std::string, double> globalParameters);
/**--------------------------------------------------------------------------------------- /**---------------------------------------------------------------------------------------
......
platforms/reference/include/ReferenceForce.h
View file @ 54e47d2d
/* Portions copyright (c) 2006-2013 Stanford University and Simbios.
/* Portions copyright (c) 2006 Stanford University and Simbios.
* Contributors: Pande Group * Contributors: Pande Group
* *
* Permission is hereby granted, free of charge, to any person obtaining * Permission is hereby granted, free of charge, to any person obtaining
...@@ -25,6 +24,7 @@ ...@@ -25,6 +24,7 @@
#ifndef __ReferenceForce_H__ #ifndef __ReferenceForce_H__
#define __ReferenceForce_H__ #define __ReferenceForce_H__
#include "lepton/CompiledExpression.h"
#include "openmm/internal/windowsExport.h" #include "openmm/internal/windowsExport.h"
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -109,6 +109,18 @@ class OPENMM_EXPORT ReferenceForce { ...@@ -109,6 +109,18 @@ class OPENMM_EXPORT ReferenceForce {
static void getDeltaROnly( const RealOpenMM* atomCoordinatesI, const RealOpenMM* atomCoordinatesJ, static void getDeltaROnly( const RealOpenMM* atomCoordinatesI, const RealOpenMM* atomCoordinatesJ,
RealOpenMM* deltaR ); RealOpenMM* deltaR );
/**
* Get a pointer to the memory for setting a variable in a CompiledExpression. If the expression
* does not use the specified variable, return NULL.
*/
static double* getVariablePointer(Lepton::CompiledExpression& expression, const std::string& name);
/**
* Set the value of a variable in a CompiledExpression, using the pointer that was returned by getVariablePointer().
*/
static void setVariable(double* pointer, double value);
}; };
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
......
platforms/reference/include/ReferenceKernels.h
View file @ 54e47d2d
...@@ -36,6 +36,7 @@ ...@@ -36,6 +36,7 @@
#include "openmm/kernels.h" #include "openmm/kernels.h"
#include "SimTKOpenMMRealType.h" #include "SimTKOpenMMRealType.h"
#include "ReferenceNeighborList.h" #include "ReferenceNeighborList.h"
#include "lepton/CompiledExpression.h"
#include "lepton/ExpressionProgram.h" #include "lepton/ExpressionProgram.h"
class CpuObc; class CpuObc;
...@@ -316,7 +317,7 @@ private: ...@@ -316,7 +317,7 @@ private:
int numBonds; int numBonds;
int **bondIndexArray; int **bondIndexArray;
RealOpenMM **bondParamArray; RealOpenMM **bondParamArray;
Lepton::ExpressionProgram energyExpression, forceExpression; Lepton::CompiledExpression energyExpression, forceExpression;
std::vector<std::string> parameterNames, globalParameterNames; std::vector<std::string> parameterNames, globalParameterNames;
}; };
...@@ -392,7 +393,7 @@ private: ...@@ -392,7 +393,7 @@ private:
int numAngles; int numAngles;
int **angleIndexArray; int **angleIndexArray;
RealOpenMM **angleParamArray; RealOpenMM **angleParamArray;
Lepton::ExpressionProgram energyExpression, forceExpression; Lepton::CompiledExpression energyExpression, forceExpression;
std::vector<std::string> parameterNames, globalParameterNames; std::vector<std::string> parameterNames, globalParameterNames;
}; };
...@@ -534,7 +535,7 @@ private: ...@@ -534,7 +535,7 @@ private:
int numTorsions; int numTorsions;
int **torsionIndexArray; int **torsionIndexArray;
RealOpenMM **torsionParamArray; RealOpenMM **torsionParamArray;
Lepton::ExpressionProgram energyExpression, forceExpression; Lepton::CompiledExpression energyExpression, forceExpression;
std::vector<std::string> parameterNames, globalParameterNames; std::vector<std::string> parameterNames, globalParameterNames;
}; };
...@@ -621,7 +622,7 @@ private: ...@@ -621,7 +622,7 @@ private:
CustomNonbondedForce* forceCopy; CustomNonbondedForce* forceCopy;
std::map<std::string, double> globalParamValues; std::map<std::string, double> globalParamValues;
std::vector<std::set<int> > exclusions; std::vector<std::set<int> > exclusions;
Lepton::ExpressionProgram energyExpression, forceExpression; Lepton::CompiledExpression energyExpression, forceExpression;
std::vector<std::string> parameterNames, globalParameterNames; std::vector<std::string> parameterNames, globalParameterNames;
std::vector<std::pair<std::set<int>, std::set<int> > > interactionGroups; std::vector<std::pair<std::set<int>, std::set<int> > > interactionGroups;
NonbondedMethod nonbondedMethod; NonbondedMethod nonbondedMethod;
...@@ -781,7 +782,7 @@ private: ...@@ -781,7 +782,7 @@ private:
int numParticles; int numParticles;
std::vector<int> particles; std::vector<int> particles;
RealOpenMM **particleParamArray; RealOpenMM **particleParamArray;
Lepton::ExpressionProgram energyExpression, forceExpressionX, forceExpressionY, forceExpressionZ; Lepton::CompiledExpression energyExpression, forceExpressionX, forceExpressionY, forceExpressionZ;
std::vector<std::string> parameterNames, globalParameterNames; std::vector<std::string> parameterNames, globalParameterNames;
}; };
......
platforms/reference/src/ReferenceKernels.cpp
View file @ 54e47d2d
...@@ -420,8 +420,8 @@ void ReferenceCalcCustomBondForceKernel::initialize(const System& system, const ...@@ -420,8 +420,8 @@ void ReferenceCalcCustomBondForceKernel::initialize(const System& system, const
// Parse the expression used to calculate the force. // Parse the expression used to calculate the force.
Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction()).optimize(); Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction()).optimize();
energyExpression = expression.createProgram(); energyExpression = expression.createCompiledExpression();
forceExpression = expression.differentiate("r").optimize().createProgram(); forceExpression = expression.differentiate("r").createCompiledExpression();
for (int i = 0; i < numParameters; i++) for (int i = 0; i < numParameters; i++)
parameterNames.push_back(force.getPerBondParameterName(i)); parameterNames.push_back(force.getPerBondParameterName(i));
for (int i = 0; i < force.getNumGlobalParameters(); i++) for (int i = 0; i < force.getNumGlobalParameters(); i++)
...@@ -534,8 +534,8 @@ void ReferenceCalcCustomAngleForceKernel::initialize(const System& system, const ...@@ -534,8 +534,8 @@ void ReferenceCalcCustomAngleForceKernel::initialize(const System& system, const
// Parse the expression used to calculate the force. // Parse the expression used to calculate the force.
Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction()).optimize(); Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction()).optimize();
energyExpression = expression.createProgram(); energyExpression = expression.createCompiledExpression();
forceExpression = expression.differentiate("theta").optimize().createProgram(); forceExpression = expression.differentiate("theta").createCompiledExpression();
for (int i = 0; i < numParameters; i++) for (int i = 0; i < numParameters; i++)
parameterNames.push_back(force.getPerAngleParameterName(i)); parameterNames.push_back(force.getPerAngleParameterName(i));
for (int i = 0; i < force.getNumGlobalParameters(); i++) for (int i = 0; i < force.getNumGlobalParameters(); i++)
...@@ -744,8 +744,8 @@ void ReferenceCalcCustomTorsionForceKernel::initialize(const System& system, con ...@@ -744,8 +744,8 @@ void ReferenceCalcCustomTorsionForceKernel::initialize(const System& system, con
// Parse the expression used to calculate the force. // Parse the expression used to calculate the force.
Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction()).optimize(); Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction()).optimize();
energyExpression = expression.createProgram(); energyExpression = expression.createCompiledExpression();
forceExpression = expression.differentiate("theta").optimize().createProgram(); forceExpression = expression.differentiate("theta").createCompiledExpression();
for (int i = 0; i < numParameters; i++) for (int i = 0; i < numParameters; i++)
parameterNames.push_back(force.getPerTorsionParameterName(i)); parameterNames.push_back(force.getPerTorsionParameterName(i));
for (int i = 0; i < force.getNumGlobalParameters(); i++) for (int i = 0; i < force.getNumGlobalParameters(); i++)
...@@ -1028,8 +1028,8 @@ void ReferenceCalcCustomNonbondedForceKernel::initialize(const System& system, c ...@@ -1028,8 +1028,8 @@ void ReferenceCalcCustomNonbondedForceKernel::initialize(const System& system, c
// Parse the various expressions used to calculate the force. // Parse the various expressions used to calculate the force.
Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction(), functions).optimize(); Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction(), functions).optimize();
energyExpression = expression.createProgram(); energyExpression = expression.createCompiledExpression();
forceExpression = expression.differentiate("r").optimize().createProgram(); forceExpression = expression.differentiate("r").createCompiledExpression();
for (int i = 0; i < numParameters; i++) for (int i = 0; i < numParameters; i++)
parameterNames.push_back(force.getPerParticleParameterName(i)); parameterNames.push_back(force.getPerParticleParameterName(i));
for (int i = 0; i < force.getNumGlobalParameters(); i++) { for (int i = 0; i < force.getNumGlobalParameters(); i++) {
...@@ -1426,10 +1426,10 @@ void ReferenceCalcCustomExternalForceKernel::initialize(const System& system, co ...@@ -1426,10 +1426,10 @@ void ReferenceCalcCustomExternalForceKernel::initialize(const System& system, co
// Parse the expression used to calculate the force. // Parse the expression used to calculate the force.
Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction()).optimize(); Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction()).optimize();
energyExpression = expression.createProgram(); energyExpression = expression.createCompiledExpression();
forceExpressionX = expression.differentiate("x").optimize().createProgram(); forceExpressionX = expression.differentiate("x").createCompiledExpression();
forceExpressionY = expression.differentiate("y").optimize().createProgram(); forceExpressionY = expression.differentiate("y").createCompiledExpression();
forceExpressionZ = expression.differentiate("z").optimize().createProgram(); forceExpressionZ = expression.differentiate("z").createCompiledExpression();
for (int i = 0; i < numParameters; i++) for (int i = 0; i < numParameters; i++)
parameterNames.push_back(force.getPerParticleParameterName(i)); parameterNames.push_back(force.getPerParticleParameterName(i));
for (int i = 0; i < force.getNumGlobalParameters(); i++) for (int i = 0; i < force.getNumGlobalParameters(); i++)
......
platforms/reference/src/SimTKReference/ReferenceCustomAngleIxn.cpp
View file @ 54e47d2d
/* Portions copyright (c) 2010 Stanford University and Simbios. /* Portions copyright (c) 2010-2013 Stanford University and Simbios.
* Contributors: Peter Eastman * Contributors: Peter Eastman
* *
* Permission is hereby granted, free of charge, to any person obtaining * Permission is hereby granted, free of charge, to any person obtaining
...@@ -39,16 +39,21 @@ using namespace std; ...@@ -39,16 +39,21 @@ using namespace std;
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
ReferenceCustomAngleIxn::ReferenceCustomAngleIxn(const Lepton::ExpressionProgram& energyExpression, ReferenceCustomAngleIxn::ReferenceCustomAngleIxn(const Lepton::CompiledExpression& energyExpression,
const Lepton::ExpressionProgram& forceExpression, const vector<string>& parameterNames, map<string, double> globalParameters) : const Lepton::CompiledExpression& forceExpression, const vector<string>& parameterNames, map<string, double> globalParameters) :
energyExpression(energyExpression), forceExpression(forceExpression), paramNames(parameterNames), globalParameters(globalParameters) { energyExpression(energyExpression), forceExpression(forceExpression) {
// --------------------------------------------------------------------------------------- energyTheta = ReferenceForce::getVariablePointer(this->energyExpression, "theta");
forceTheta = ReferenceForce::getVariablePointer(this->forceExpression, "theta");
// static const char* methodName = "\nReferenceCustomAngleIxn::ReferenceCustomAngleIxn"; numParameters = parameterNames.size();
for (int i = 0; i < (int) numParameters; i++) {
// --------------------------------------------------------------------------------------- energyParams.push_back(ReferenceForce::getVariablePointer(this->energyExpression, parameterNames[i]));
forceParams.push_back(ReferenceForce::getVariablePointer(this->forceExpression, parameterNames[i]));
}
for (map<string, double>::const_iterator iter = globalParameters.begin(); iter != globalParameters.end(); ++iter) {
ReferenceForce::setVariable(ReferenceForce::getVariablePointer(this->energyExpression, iter->first), iter->second);
ReferenceForce::setVariable(ReferenceForce::getVariablePointer(this->forceExpression, iter->first), iter->second);
}
} }
/**--------------------------------------------------------------------------------------- /**---------------------------------------------------------------------------------------
...@@ -91,9 +96,10 @@ void ReferenceCustomAngleIxn::calculateBondIxn( int* atomIndices, ...@@ -91,9 +96,10 @@ void ReferenceCustomAngleIxn::calculateBondIxn( int* atomIndices,
static const RealOpenMM one = 1.0; static const RealOpenMM one = 1.0;
RealOpenMM deltaR[2][ReferenceForce::LastDeltaRIndex]; RealOpenMM deltaR[2][ReferenceForce::LastDeltaRIndex];
map<string, double> variables = globalParameters; for (int i = 0; i < numParameters; i++) {
for (int i = 0; i < (int) paramNames.size(); ++i) ReferenceForce::setVariable(energyParams[i], parameters[i]);
variables[paramNames[i]] = parameters[i]; ReferenceForce::setVariable(forceParams[i], parameters[i]);
}
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -118,12 +124,13 @@ void ReferenceCustomAngleIxn::calculateBondIxn( int* atomIndices, ...@@ -118,12 +124,13 @@ void ReferenceCustomAngleIxn::calculateBondIxn( int* atomIndices,
angle = PI_M; angle = PI_M;
else else
angle = ACOS(cosine); angle = ACOS(cosine);
variables["theta"] = angle; ReferenceForce::setVariable(energyTheta, angle);
ReferenceForce::setVariable(forceTheta, angle);
// Compute the force and energy, and apply them to the atoms. // Compute the force and energy, and apply them to the atoms.
RealOpenMM energy = (RealOpenMM) energyExpression.evaluate(variables); RealOpenMM energy = (RealOpenMM) energyExpression.evaluate();
RealOpenMM dEdR = (RealOpenMM) forceExpression.evaluate(variables); RealOpenMM dEdR = (RealOpenMM) forceExpression.evaluate();
RealOpenMM termA = dEdR/(deltaR[0][ReferenceForce::R2Index]*rp); RealOpenMM termA = dEdR/(deltaR[0][ReferenceForce::R2Index]*rp);
RealOpenMM termC = -dEdR/(deltaR[1][ReferenceForce::R2Index]*rp); RealOpenMM termC = -dEdR/(deltaR[1][ReferenceForce::R2Index]*rp);
......
platforms/reference/src/SimTKReference/ReferenceCustomBondIxn.cpp
View file @ 54e47d2d
/* Portions copyright (c) 2009 Stanford University and Simbios. /* Portions copyright (c) 2009-2013 Stanford University and Simbios.
* Contributors: Peter Eastman * Contributors: Peter Eastman
* *
* Permission is hereby granted, free of charge, to any person obtaining * Permission is hereby granted, free of charge, to any person obtaining
...@@ -40,16 +40,20 @@ using namespace OpenMM; ...@@ -40,16 +40,20 @@ using namespace OpenMM;
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
ReferenceCustomBondIxn::ReferenceCustomBondIxn(const Lepton::ExpressionProgram& energyExpression, ReferenceCustomBondIxn::ReferenceCustomBondIxn(const Lepton::CompiledExpression& energyExpression,
const Lepton::ExpressionProgram& forceExpression, const vector<string>& parameterNames, map<string, double> globalParameters) : const Lepton::CompiledExpression& forceExpression, const vector<string>& parameterNames, map<string, double> globalParameters) :
energyExpression(energyExpression), forceExpression(forceExpression), paramNames(parameterNames), globalParameters(globalParameters) { energyExpression(energyExpression), forceExpression(forceExpression) {
energyR = ReferenceForce::getVariablePointer(this->energyExpression, "r");
// --------------------------------------------------------------------------------------- forceR = ReferenceForce::getVariablePointer(this->forceExpression, "r");
numParameters = parameterNames.size();
// static const char* methodName = "\nReferenceCustomBondIxn::ReferenceCustomBondIxn"; for (int i = 0; i < (int) numParameters; i++) {
energyParams.push_back(ReferenceForce::getVariablePointer(this->energyExpression, parameterNames[i]));
// --------------------------------------------------------------------------------------- forceParams.push_back(ReferenceForce::getVariablePointer(this->forceExpression, parameterNames[i]));
}
for (map<string, double>::const_iterator iter = globalParameters.begin(); iter != globalParameters.end(); ++iter) {
ReferenceForce::setVariable(ReferenceForce::getVariablePointer(this->energyExpression, iter->first), iter->second);
ReferenceForce::setVariable(ReferenceForce::getVariablePointer(this->forceExpression, iter->first), iter->second);
}
} }
/**--------------------------------------------------------------------------------------- /**---------------------------------------------------------------------------------------
...@@ -95,9 +99,10 @@ void ReferenceCustomBondIxn::calculateBondIxn( int* atomIndices, ...@@ -95,9 +99,10 @@ void ReferenceCustomBondIxn::calculateBondIxn( int* atomIndices,
static const RealOpenMM half = 0.5; static const RealOpenMM half = 0.5;
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex]; RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
map<string, double> variables = globalParameters; for (int i = 0; i < numParameters; i++) {
for (int i = 0; i < (int) paramNames.size(); ++i) ReferenceForce::setVariable(energyParams[i], parameters[i]);
variables[paramNames[i]] = parameters[i]; ReferenceForce::setVariable(forceParams[i], parameters[i]);
}
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -106,8 +111,10 @@ void ReferenceCustomBondIxn::calculateBondIxn( int* atomIndices, ...@@ -106,8 +111,10 @@ void ReferenceCustomBondIxn::calculateBondIxn( int* atomIndices,
int atomAIndex = atomIndices[0]; int atomAIndex = atomIndices[0];
int atomBIndex = atomIndices[1]; int atomBIndex = atomIndices[1];
ReferenceForce::getDeltaR( atomCoordinates[atomAIndex], atomCoordinates[atomBIndex], deltaR ); ReferenceForce::getDeltaR( atomCoordinates[atomAIndex], atomCoordinates[atomBIndex], deltaR );
variables["r"] = deltaR[ReferenceForce::RIndex];
RealOpenMM dEdR = (RealOpenMM) forceExpression.evaluate(variables); ReferenceForce::setVariable(energyR, deltaR[ReferenceForce::RIndex]);
ReferenceForce::setVariable(forceR, deltaR[ReferenceForce::RIndex]);
RealOpenMM dEdR = (RealOpenMM) forceExpression.evaluate();
dEdR = deltaR[ReferenceForce::RIndex] > zero ? (dEdR/deltaR[ReferenceForce::RIndex]) : zero; dEdR = deltaR[ReferenceForce::RIndex] > zero ? (dEdR/deltaR[ReferenceForce::RIndex]) : zero;
forces[atomAIndex][0] += dEdR*deltaR[ReferenceForce::XIndex]; forces[atomAIndex][0] += dEdR*deltaR[ReferenceForce::XIndex];
...@@ -119,5 +126,5 @@ void ReferenceCustomBondIxn::calculateBondIxn( int* atomIndices, ...@@ -119,5 +126,5 @@ void ReferenceCustomBondIxn::calculateBondIxn( int* atomIndices,
forces[atomBIndex][2] -= dEdR*deltaR[ReferenceForce::ZIndex]; forces[atomBIndex][2] -= dEdR*deltaR[ReferenceForce::ZIndex];
if (totalEnergy != NULL) if (totalEnergy != NULL)
*totalEnergy += (RealOpenMM) energyExpression.evaluate(variables); *totalEnergy += (RealOpenMM) energyExpression.evaluate();
} }
platforms/reference/src/SimTKReference/ReferenceCustomExternalIxn.cpp
View file @ 54e47d2d
/* Portions copyright (c) 2009 Stanford University and Simbios. /* Portions copyright (c) 2009-2013 Stanford University and Simbios.
* Contributors: Peter Eastman * Contributors: Peter Eastman
* *
* Permission is hereby granted, free of charge, to any person obtaining * Permission is hereby granted, free of charge, to any person obtaining
...@@ -40,18 +40,37 @@ using namespace OpenMM; ...@@ -40,18 +40,37 @@ using namespace OpenMM;
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
ReferenceCustomExternalIxn::ReferenceCustomExternalIxn(const Lepton::ExpressionProgram& energyExpression, ReferenceCustomExternalIxn::ReferenceCustomExternalIxn(const Lepton::CompiledExpression& energyExpression,
const Lepton::ExpressionProgram& forceExpressionX, const Lepton::ExpressionProgram& forceExpressionY, const Lepton::CompiledExpression& forceExpressionX, const Lepton::CompiledExpression& forceExpressionY,
const Lepton::ExpressionProgram& forceExpressionZ, const vector<string>& parameterNames, map<string, double> globalParameters) : const Lepton::CompiledExpression& forceExpressionZ, const vector<string>& parameterNames, map<string, double> globalParameters) :
energyExpression(energyExpression), forceExpressionX(forceExpressionX), forceExpressionY(forceExpressionY), energyExpression(energyExpression), forceExpressionX(forceExpressionX), forceExpressionY(forceExpressionY),
forceExpressionZ(forceExpressionZ), paramNames(parameterNames), globalParameters(globalParameters) { forceExpressionZ(forceExpressionZ) {
// --------------------------------------------------------------------------------------- energyX = ReferenceForce::getVariablePointer(this->energyExpression, "x");
energyY = ReferenceForce::getVariablePointer(this->energyExpression, "y");
// static const char* methodName = "\nReferenceCustomExternalIxn::ReferenceCustomExternalIxn"; energyZ = ReferenceForce::getVariablePointer(this->energyExpression, "z");
forceXX = ReferenceForce::getVariablePointer(this->forceExpressionX, "x");
// --------------------------------------------------------------------------------------- forceXY = ReferenceForce::getVariablePointer(this->forceExpressionX, "y");
forceXZ = ReferenceForce::getVariablePointer(this->forceExpressionX, "z");
forceYX = ReferenceForce::getVariablePointer(this->forceExpressionY, "x");
forceYY = ReferenceForce::getVariablePointer(this->forceExpressionY, "y");
forceYZ = ReferenceForce::getVariablePointer(this->forceExpressionY, "z");
forceZX = ReferenceForce::getVariablePointer(this->forceExpressionZ, "x");
forceZY = ReferenceForce::getVariablePointer(this->forceExpressionZ, "y");
forceZZ = ReferenceForce::getVariablePointer(this->forceExpressionZ, "z");
numParameters = parameterNames.size();
for (int i = 0; i < (int) numParameters; i++) {
energyParams.push_back(ReferenceForce::getVariablePointer(this->energyExpression, parameterNames[i]));
forceXParams.push_back(ReferenceForce::getVariablePointer(this->forceExpressionX, parameterNames[i]));
forceYParams.push_back(ReferenceForce::getVariablePointer(this->forceExpressionY, parameterNames[i]));
forceZParams.push_back(ReferenceForce::getVariablePointer(this->forceExpressionZ, parameterNames[i]));
}
for (map<string, double>::const_iterator iter = globalParameters.begin(); iter != globalParameters.end(); ++iter) {
ReferenceForce::setVariable(ReferenceForce::getVariablePointer(this->energyExpression, iter->first), iter->second);
ReferenceForce::setVariable(ReferenceForce::getVariablePointer(this->forceExpressionX, iter->first), iter->second);
ReferenceForce::setVariable(ReferenceForce::getVariablePointer(this->forceExpressionY, iter->first), iter->second);
ReferenceForce::setVariable(ReferenceForce::getVariablePointer(this->forceExpressionZ, iter->first), iter->second);
}
} }
/**--------------------------------------------------------------------------------------- /**---------------------------------------------------------------------------------------
...@@ -90,18 +109,30 @@ void ReferenceCustomExternalIxn::calculateForce( int atomIndex, ...@@ -90,18 +109,30 @@ void ReferenceCustomExternalIxn::calculateForce( int atomIndex,
static const std::string methodName = "\nReferenceCustomExternalIxn::calculateBondIxn"; static const std::string methodName = "\nReferenceCustomExternalIxn::calculateBondIxn";
map<string, double> variables = globalParameters; for (int i = 0; i < numParameters; i++) {
for (int i = 0; i < (int) paramNames.size(); ++i) ReferenceForce::setVariable(energyParams[i], parameters[i]);
variables[paramNames[i]] = parameters[i]; ReferenceForce::setVariable(forceXParams[i], parameters[i]);
variables["x"] = atomCoordinates[atomIndex][0]; ReferenceForce::setVariable(forceYParams[i], parameters[i]);
variables["y"] = atomCoordinates[atomIndex][1]; ReferenceForce::setVariable(forceZParams[i], parameters[i]);
variables["z"] = atomCoordinates[atomIndex][2]; }
ReferenceForce::setVariable(energyX, atomCoordinates[atomIndex][0]);
ReferenceForce::setVariable(energyY, atomCoordinates[atomIndex][1]);
ReferenceForce::setVariable(energyZ, atomCoordinates[atomIndex][2]);
ReferenceForce::setVariable(forceXX, atomCoordinates[atomIndex][0]);
ReferenceForce::setVariable(forceXY, atomCoordinates[atomIndex][1]);
ReferenceForce::setVariable(forceXZ, atomCoordinates[atomIndex][2]);
ReferenceForce::setVariable(forceYX, atomCoordinates[atomIndex][0]);
ReferenceForce::setVariable(forceYY, atomCoordinates[atomIndex][1]);
ReferenceForce::setVariable(forceYZ, atomCoordinates[atomIndex][2]);
ReferenceForce::setVariable(forceZX, atomCoordinates[atomIndex][0]);
ReferenceForce::setVariable(forceZY, atomCoordinates[atomIndex][1]);
ReferenceForce::setVariable(forceZZ, atomCoordinates[atomIndex][2]);
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
forces[atomIndex][0] -= (RealOpenMM) forceExpressionX.evaluate(variables); forces[atomIndex][0] -= (RealOpenMM) forceExpressionX.evaluate();
forces[atomIndex][1] -= (RealOpenMM) forceExpressionY.evaluate(variables); forces[atomIndex][1] -= (RealOpenMM) forceExpressionY.evaluate();
forces[atomIndex][2] -= (RealOpenMM) forceExpressionZ.evaluate(variables); forces[atomIndex][2] -= (RealOpenMM) forceExpressionZ.evaluate();
if (energy != NULL) if (energy != NULL)
*energy += (RealOpenMM) energyExpression.evaluate(variables); *energy += (RealOpenMM) energyExpression.evaluate();
} }
platforms/reference/src/SimTKReference/ReferenceCustomNonbondedIxn.cpp
View file @ 54e47d2d
/* Portions copyright (c) 2009 Stanford University and Simbios. /* Portions copyright (c) 2009-2013 Stanford University and Simbios.
* Contributors: Peter Eastman * Contributors: Peter Eastman
* *
* Permission is hereby granted, free of charge, to any person obtaining * Permission is hereby granted, free of charge, to any person obtaining
...@@ -45,8 +45,8 @@ using OpenMM::RealVec; ...@@ -45,8 +45,8 @@ using OpenMM::RealVec;
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
ReferenceCustomNonbondedIxn::ReferenceCustomNonbondedIxn(const Lepton::ExpressionProgram& energyExpression, ReferenceCustomNonbondedIxn::ReferenceCustomNonbondedIxn(const Lepton::CompiledExpression& energyExpression,
const Lepton::ExpressionProgram& forceExpression, const vector<string>& parameterNames) : const Lepton::CompiledExpression& forceExpression, const vector<string>& parameterNames) :
cutoff(false), useSwitch(false), periodic(false), energyExpression(energyExpression), forceExpression(forceExpression), paramNames(parameterNames) { cutoff(false), useSwitch(false), periodic(false), energyExpression(energyExpression), forceExpression(forceExpression), paramNames(parameterNames) {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -55,11 +55,14 @@ ReferenceCustomNonbondedIxn::ReferenceCustomNonbondedIxn(const Lepton::Expressio ...@@ -55,11 +55,14 @@ ReferenceCustomNonbondedIxn::ReferenceCustomNonbondedIxn(const Lepton::Expressio
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
energyR = ReferenceForce::getVariablePointer(this->energyExpression, "r");
forceR = ReferenceForce::getVariablePointer(this->forceExpression, "r");
for (int i = 0; i < (int) paramNames.size(); i++) { for (int i = 0; i < (int) paramNames.size(); i++) {
for (int j = 1; j < 3; j++) { for (int j = 1; j < 3; j++) {
stringstream name; stringstream name;
name << paramNames[i] << j; name << paramNames[i] << j;
particleParamNames.push_back(name.str()); energyParticleParams.push_back(ReferenceForce::getVariablePointer(this->energyExpression, name.str()));
forceParticleParams.push_back(ReferenceForce::getVariablePointer(this->forceExpression, name.str()));
} }
} }
} }
...@@ -166,9 +169,12 @@ void ReferenceCustomNonbondedIxn::setUseSwitchingFunction( RealOpenMM distance ) ...@@ -166,9 +169,12 @@ void ReferenceCustomNonbondedIxn::setUseSwitchingFunction( RealOpenMM distance )
void ReferenceCustomNonbondedIxn::calculatePairIxn( int numberOfAtoms, vector<RealVec>& atomCoordinates, void ReferenceCustomNonbondedIxn::calculatePairIxn( int numberOfAtoms, vector<RealVec>& atomCoordinates,
RealOpenMM** atomParameters, vector<set<int> >& exclusions, RealOpenMM** atomParameters, vector<set<int> >& exclusions,
RealOpenMM* fixedParameters, const map<string, double>& globalParameters, vector<RealVec>& forces, RealOpenMM* fixedParameters, const map<string, double>& globalParameters, vector<RealVec>& forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const { RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) {
map<string, double> variables = globalParameters; for (map<string, double>::const_iterator iter = globalParameters.begin(); iter != globalParameters.end(); ++iter) {
ReferenceForce::setVariable(ReferenceForce::getVariablePointer(energyExpression, iter->first), iter->second);
ReferenceForce::setVariable(ReferenceForce::getVariablePointer(forceExpression, iter->first), iter->second);
}
if (interactionGroups.size() > 0) { if (interactionGroups.size() > 0) {
// The user has specified interaction groups, so compute only the requested interactions. // The user has specified interaction groups, so compute only the requested interactions.
...@@ -182,10 +188,12 @@ void ReferenceCustomNonbondedIxn::calculatePairIxn( int numberOfAtoms, vector<Re ...@@ -182,10 +188,12 @@ void ReferenceCustomNonbondedIxn::calculatePairIxn( int numberOfAtoms, vector<Re
if (*atom1 > *atom2 && set1.find(*atom2) != set1.end() && set2.find(*atom1) != set2.end()) if (*atom1 > *atom2 && set1.find(*atom2) != set1.end() && set2.find(*atom1) != set2.end())
continue; // Both atoms are in both sets, so skip duplicate interactions. continue; // Both atoms are in both sets, so skip duplicate interactions.
for (int j = 0; j < (int) paramNames.size(); j++) { for (int j = 0; j < (int) paramNames.size(); j++) {
variables[particleParamNames[j*2]] = atomParameters[*atom1][j]; ReferenceForce::setVariable(energyParticleParams[j*2], atomParameters[*atom1][j]);
variables[particleParamNames[j*2+1]] = atomParameters[*atom2][j]; ReferenceForce::setVariable(energyParticleParams[j*2+1], atomParameters[*atom2][j]);
ReferenceForce::setVariable(forceParticleParams[j*2], atomParameters[*atom1][j]);
ReferenceForce::setVariable(forceParticleParams[j*2+1], atomParameters[*atom2][j]);
} }
calculateOneIxn(*atom1, *atom2, atomCoordinates, variables, forces, energyByAtom, totalEnergy); calculateOneIxn(*atom1, *atom2, atomCoordinates, forces, energyByAtom, totalEnergy);
} }
} }
} }
...@@ -196,10 +204,12 @@ void ReferenceCustomNonbondedIxn::calculatePairIxn( int numberOfAtoms, vector<Re ...@@ -196,10 +204,12 @@ void ReferenceCustomNonbondedIxn::calculatePairIxn( int numberOfAtoms, vector<Re
for (int i = 0; i < (int) neighborList->size(); i++) { for (int i = 0; i < (int) neighborList->size(); i++) {
OpenMM::AtomPair pair = (*neighborList)[i]; OpenMM::AtomPair pair = (*neighborList)[i];
for (int j = 0; j < (int) paramNames.size(); j++) { for (int j = 0; j < (int) paramNames.size(); j++) {
variables[particleParamNames[j*2]] = atomParameters[pair.first][j]; ReferenceForce::setVariable(energyParticleParams[j*2], atomParameters[pair.first][j]);
variables[particleParamNames[j*2+1]] = atomParameters[pair.second][j]; ReferenceForce::setVariable(energyParticleParams[j*2+1], atomParameters[pair.second][j]);
ReferenceForce::setVariable(forceParticleParams[j*2], atomParameters[pair.first][j]);
ReferenceForce::setVariable(forceParticleParams[j*2+1], atomParameters[pair.second][j]);
} }
calculateOneIxn(pair.first, pair.second, atomCoordinates, variables, forces, energyByAtom, totalEnergy); calculateOneIxn(pair.first, pair.second, atomCoordinates, forces, energyByAtom, totalEnergy);
} }
} }
else { else {
...@@ -209,10 +219,12 @@ void ReferenceCustomNonbondedIxn::calculatePairIxn( int numberOfAtoms, vector<Re ...@@ -209,10 +219,12 @@ void ReferenceCustomNonbondedIxn::calculatePairIxn( int numberOfAtoms, vector<Re
for (int jj = ii+1; jj < numberOfAtoms; jj++) { for (int jj = ii+1; jj < numberOfAtoms; jj++) {
if (exclusions[jj].find(ii) == exclusions[jj].end()) { if (exclusions[jj].find(ii) == exclusions[jj].end()) {
for (int j = 0; j < (int) paramNames.size(); j++) { for (int j = 0; j < (int) paramNames.size(); j++) {
variables[particleParamNames[j*2]] = atomParameters[ii][j]; ReferenceForce::setVariable(energyParticleParams[j*2], atomParameters[ii][j]);
variables[particleParamNames[j*2+1]] = atomParameters[jj][j]; ReferenceForce::setVariable(energyParticleParams[j*2+1], atomParameters[jj][j]);
ReferenceForce::setVariable(forceParticleParams[j*2], atomParameters[ii][j]);
ReferenceForce::setVariable(forceParticleParams[j*2+1], atomParameters[jj][j]);
} }
calculateOneIxn(ii, jj, atomCoordinates, variables, forces, energyByAtom, totalEnergy); calculateOneIxn(ii, jj, atomCoordinates, forces, energyByAtom, totalEnergy);
} }
} }
} }
...@@ -233,9 +245,8 @@ void ReferenceCustomNonbondedIxn::calculatePairIxn( int numberOfAtoms, vector<Re ...@@ -233,9 +245,8 @@ void ReferenceCustomNonbondedIxn::calculatePairIxn( int numberOfAtoms, vector<Re
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
void ReferenceCustomNonbondedIxn::calculateOneIxn( int ii, int jj, vector<RealVec>& atomCoordinates, void ReferenceCustomNonbondedIxn::calculateOneIxn( int ii, int jj, vector<RealVec>& atomCoordinates, vector<RealVec>& forces,
map<string, double>& variables, vector<RealVec>& forces, RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) {
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -266,9 +277,10 @@ void ReferenceCustomNonbondedIxn::calculateOneIxn( int ii, int jj, vector<RealVe ...@@ -266,9 +277,10 @@ void ReferenceCustomNonbondedIxn::calculateOneIxn( int ii, int jj, vector<RealVe
// accumulate forces // accumulate forces
variables["r"] = r; ReferenceForce::setVariable(energyR, r);
RealOpenMM dEdR = (RealOpenMM) (forceExpression.evaluate(variables)/(deltaR[ReferenceForce::RIndex])); ReferenceForce::setVariable(forceR, r);
RealOpenMM energy = (RealOpenMM) energyExpression.evaluate(variables); RealOpenMM dEdR = (RealOpenMM) (forceExpression.evaluate()/(deltaR[ReferenceForce::RIndex]));
RealOpenMM energy = (RealOpenMM) energyExpression.evaluate();
if (useSwitch) { if (useSwitch) {
if (r > switchingDistance) { if (r > switchingDistance) {
RealOpenMM t = (r-switchingDistance)/(cutoffDistance-switchingDistance); RealOpenMM t = (r-switchingDistance)/(cutoffDistance-switchingDistance);
......
platforms/reference/src/SimTKReference/ReferenceCustomTorsionIxn.cpp
View file @ 54e47d2d
/* Portions copyright (c) 2010 Stanford University and Simbios. /* Portions copyright (c) 2010-2013 Stanford University and Simbios.
* Contributors: Peter Eastman * Contributors: Peter Eastman
* *
* Permission is hereby granted, free of charge, to any person obtaining * Permission is hereby granted, free of charge, to any person obtaining
...@@ -39,16 +39,21 @@ using namespace OpenMM; ...@@ -39,16 +39,21 @@ using namespace OpenMM;
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
ReferenceCustomTorsionIxn::ReferenceCustomTorsionIxn(const Lepton::ExpressionProgram& energyExpression, ReferenceCustomTorsionIxn::ReferenceCustomTorsionIxn(const Lepton::CompiledExpression& energyExpression,
const Lepton::ExpressionProgram& forceExpression, const vector<string>& parameterNames, map<string, double> globalParameters) : const Lepton::CompiledExpression& forceExpression, const vector<string>& parameterNames, map<string, double> globalParameters) :
energyExpression(energyExpression), forceExpression(forceExpression), paramNames(parameterNames), globalParameters(globalParameters) { energyExpression(energyExpression), forceExpression(forceExpression) {
// --------------------------------------------------------------------------------------- energyTheta = ReferenceForce::getVariablePointer(this->energyExpression, "theta");
forceTheta = ReferenceForce::getVariablePointer(this->forceExpression, "theta");
// static const char* methodName = "\nReferenceCustomTorsionIxn::ReferenceCustomTorsionIxn"; numParameters = parameterNames.size();
for (int i = 0; i < (int) numParameters; i++) {
// --------------------------------------------------------------------------------------- energyParams.push_back(ReferenceForce::getVariablePointer(this->energyExpression, parameterNames[i]));
forceParams.push_back(ReferenceForce::getVariablePointer(this->forceExpression, parameterNames[i]));
}
for (map<string, double>::const_iterator iter = globalParameters.begin(); iter != globalParameters.end(); ++iter) {
ReferenceForce::setVariable(ReferenceForce::getVariablePointer(this->energyExpression, iter->first), iter->second);
ReferenceForce::setVariable(ReferenceForce::getVariablePointer(this->forceExpression, iter->first), iter->second);
}
} }
/**--------------------------------------------------------------------------------------- /**---------------------------------------------------------------------------------------
...@@ -91,9 +96,10 @@ void ReferenceCustomTorsionIxn::calculateBondIxn( int* atomIndices, ...@@ -91,9 +96,10 @@ void ReferenceCustomTorsionIxn::calculateBondIxn( int* atomIndices,
static const RealOpenMM one = 1.0; static const RealOpenMM one = 1.0;
RealOpenMM deltaR[3][ReferenceForce::LastDeltaRIndex]; RealOpenMM deltaR[3][ReferenceForce::LastDeltaRIndex];
map<string, double> variables = globalParameters; for (int i = 0; i < numParameters; i++) {
for (int i = 0; i < (int) paramNames.size(); ++i) ReferenceForce::setVariable(energyParams[i], parameters[i]);
variables[paramNames[i]] = parameters[i]; ReferenceForce::setVariable(forceParams[i], parameters[i]);
}
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -118,13 +124,13 @@ void ReferenceCustomTorsionIxn::calculateBondIxn( int* atomIndices, ...@@ -118,13 +124,13 @@ void ReferenceCustomTorsionIxn::calculateBondIxn( int* atomIndices,
RealOpenMM dotDihedral; RealOpenMM dotDihedral;
RealOpenMM signOfAngle; RealOpenMM signOfAngle;
variables["theta"] = getDihedralAngleBetweenThreeVectors(deltaR[0], deltaR[1], deltaR[2], RealOpenMM angle = getDihedralAngleBetweenThreeVectors(deltaR[0], deltaR[1], deltaR[2], crossProduct, &dotDihedral, deltaR[0], &signOfAngle, 1);
crossProduct, &dotDihedral, deltaR[0], ReferenceForce::setVariable(energyTheta, angle);
&signOfAngle, 1); ReferenceForce::setVariable(forceTheta, angle);
// evaluate delta angle, dE/d(angle) // evaluate delta angle, dE/d(angle)
RealOpenMM dEdAngle = (RealOpenMM) forceExpression.evaluate(variables); RealOpenMM dEdAngle = (RealOpenMM) forceExpression.evaluate();
// compute force // compute force
...@@ -166,6 +172,6 @@ void ReferenceCustomTorsionIxn::calculateBondIxn( int* atomIndices, ...@@ -166,6 +172,6 @@ void ReferenceCustomTorsionIxn::calculateBondIxn( int* atomIndices,
// accumulate energies // accumulate energies
if (totalEnergy != NULL) if (totalEnergy != NULL)
*totalEnergy += (RealOpenMM) energyExpression.evaluate(variables); *totalEnergy += (RealOpenMM) energyExpression.evaluate();
} }
platforms/reference/src/SimTKReference/ReferenceForce.cpp
View file @ 54e47d2d
/* Portions copyright (c) 2006 Stanford University and Simbios. /* Portions copyright (c) 2006-2013 Stanford University and Simbios.
* Contributors: Pande Group * Contributors: Pande Group
* *
* Permission is hereby granted, free of charge, to any person obtaining * Permission is hereby granted, free of charge, to any person obtaining
...@@ -161,3 +161,14 @@ void ReferenceForce::getDeltaROnly( const RealOpenMM* atomCoordinatesI, ...@@ -161,3 +161,14 @@ void ReferenceForce::getDeltaROnly( const RealOpenMM* atomCoordinatesI,
deltaR[YIndex] = atomCoordinatesJ[1] - atomCoordinatesI[1]; deltaR[YIndex] = atomCoordinatesJ[1] - atomCoordinatesI[1];
deltaR[ZIndex] = atomCoordinatesJ[2] - atomCoordinatesI[2]; deltaR[ZIndex] = atomCoordinatesJ[2] - atomCoordinatesI[2];
} }
double* ReferenceForce::getVariablePointer(Lepton::CompiledExpression& expression, const std::string& name) {
if (expression.getVariables().find(name) == expression.getVariables().end())
return NULL;
return &expression.getVariableReference(name);
}
void ReferenceForce::setVariable(double* pointer, double value) {
if (pointer != NULL)
*pointer = value;
}
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