resolved conflict

platforms/reference/src/ReferenceKernels.cpp

@@ -280,7 +280,7 @@ void ReferenceUpdateStateDataKernel::getPeriodicBoxVectors(ContextImpl& context,
     c = vectors[2];
 }
 
-void ReferenceUpdateStateDataKernel::setPeriodicBoxVectors(ContextImpl& context, const Vec3& a, const Vec3& b, const Vec3& c) const {
+void ReferenceUpdateStateDataKernel::setPeriodicBoxVectors(ContextImpl& context, const Vec3& a, const Vec3& b, const Vec3& c) {
     RealVec& box = extractBoxSize(context);
     box[0] = (RealOpenMM) a[0];
     box[1] = (RealOpenMM) b[1];
@@ -370,6 +370,7 @@ void ReferenceCalcHarmonicBondForceKernel::initialize(const System& system, cons
         bondParamArray[i][0] = (RealOpenMM) length;
         bondParamArray[i][1] = (RealOpenMM) k;
     }
+    usePeriodic = force.usesPeriodicBoundaryConditions();
 }
 
 double ReferenceCalcHarmonicBondForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
@@ -378,6 +379,8 @@ double ReferenceCalcHarmonicBondForceKernel::execute(ContextImpl& context, bool
     RealOpenMM energy = 0;
     ReferenceBondForce refBondForce;
     ReferenceHarmonicBondIxn harmonicBond;
+    if (usePeriodic)
+        harmonicBond.setPeriodic(extractBoxVectors(context));
     refBondForce.calculateForce(numBonds, bondIndexArray, posData, bondParamArray, forceData, includeEnergy ? &energy : NULL, harmonicBond);
     return energy;
 }
@@ -409,6 +412,7 @@ ReferenceCalcCustomBondForceKernel::~ReferenceCalcCustomBondForceKernel() {
 void ReferenceCalcCustomBondForceKernel::initialize(const System& system, const CustomBondForce& force) {
     numBonds = force.getNumBonds();
     int numParameters = force.getNumPerBondParameters();
+    usePeriodic = force.usesPeriodicBoundaryConditions();
 
     // Build the arrays.
 
@@ -448,8 +452,10 @@ double ReferenceCalcCustomBondForceKernel::execute(ContextImpl& context, bool in
     for (int i = 0; i < (int) globalParameterNames.size(); i++)
         globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
     ReferenceBondForce refBondForce;
-    ReferenceCustomBondIxn harmonicBond(energyExpression, forceExpression, parameterNames, globalParameters);
-    refBondForce.calculateForce(numBonds, bondIndexArray, posData, bondParamArray, forceData, includeEnergy ? &energy : NULL, harmonicBond);
+    ReferenceCustomBondIxn bond(energyExpression, forceExpression, parameterNames, globalParameters);
+    if (usePeriodic)
+        bond.setPeriodic(extractBoxVectors(context));
+    refBondForce.calculateForce(numBonds, bondIndexArray, posData, bondParamArray, forceData, includeEnergy ? &energy : NULL, bond);
     return energy;
 }
 
@@ -490,6 +496,7 @@ void ReferenceCalcHarmonicAngleForceKernel::initialize(const System& system, con
         angleParamArray[i][0] = (RealOpenMM) angle;
         angleParamArray[i][1] = (RealOpenMM) k;
     }
+    usePeriodic = force.usesPeriodicBoundaryConditions();
 }
 
 double ReferenceCalcHarmonicAngleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
@@ -498,6 +505,8 @@ double ReferenceCalcHarmonicAngleForceKernel::execute(ContextImpl& context, bool
     RealOpenMM energy = 0;
     ReferenceBondForce refBondForce;
     ReferenceAngleBondIxn angleBond;
+    if (usePeriodic)
+        angleBond.setPeriodic(extractBoxVectors(context));
     refBondForce.calculateForce(numAngles, angleIndexArray, posData, angleParamArray, forceData, includeEnergy ? &energy : NULL, angleBond);
     return energy;
 }
@@ -527,6 +536,7 @@ ReferenceCalcCustomAngleForceKernel::~ReferenceCalcCustomAngleForceKernel() {
 void ReferenceCalcCustomAngleForceKernel::initialize(const System& system, const CustomAngleForce& force) {
     numAngles = force.getNumAngles();
     int numParameters = force.getNumPerAngleParameters();
+    usePeriodic = force.usesPeriodicBoundaryConditions();
 
     // Build the arrays.
 
@@ -568,6 +578,8 @@ double ReferenceCalcCustomAngleForceKernel::execute(ContextImpl& context, bool i
         globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
     ReferenceBondForce refBondForce;
     ReferenceCustomAngleIxn customAngle(energyExpression, forceExpression, parameterNames, globalParameters);
+    if (usePeriodic)
+        customAngle.setPeriodic(extractBoxVectors(context));
     refBondForce.calculateForce(numAngles, angleIndexArray, posData, angleParamArray, forceData, includeEnergy ? &energy : NULL, customAngle);
     return energy;
 }
@@ -611,6 +623,7 @@ void ReferenceCalcPeriodicTorsionForceKernel::initialize(const System& system, c
         torsionParamArray[i][1] = (RealOpenMM) phase;
         torsionParamArray[i][2] = (RealOpenMM) periodicity;
     }
+    usePeriodic = force.usesPeriodicBoundaryConditions();
 }
 
 double ReferenceCalcPeriodicTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
@@ -619,6 +632,8 @@ double ReferenceCalcPeriodicTorsionForceKernel::execute(ContextImpl& context, bo
     RealOpenMM energy = 0;
     ReferenceBondForce refBondForce;
     ReferenceProperDihedralBond periodicTorsionBond;
+    if (usePeriodic)
+        periodicTorsionBond.setPeriodic(extractBoxVectors(context));
     refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, includeEnergy ? &energy : NULL, periodicTorsionBond);
     return energy;
 }
@@ -665,6 +680,7 @@ void ReferenceCalcRBTorsionForceKernel::initialize(const System& system, const R
         torsionParamArray[i][4] = (RealOpenMM) c4;
         torsionParamArray[i][5] = (RealOpenMM) c5;
     }
+    usePeriodic = force.usesPeriodicBoundaryConditions();
 }
 
 double ReferenceCalcRBTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
@@ -673,6 +689,8 @@ double ReferenceCalcRBTorsionForceKernel::execute(ContextImpl& context, bool inc
     RealOpenMM energy = 0;
     ReferenceBondForce refBondForce;
     ReferenceRbDihedralBond rbTorsionBond;
+    if (usePeriodic)
+        rbTorsionBond.setPeriodic(extractBoxVectors(context));
     refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, includeEnergy ? &energy : NULL, rbTorsionBond);
     return energy;
 }
@@ -722,6 +740,7 @@ void ReferenceCalcCMAPTorsionForceKernel::initialize(const System& system, const
         force.getTorsionParameters(i, torsionMaps[i], torsionIndices[i][0], torsionIndices[i][1], torsionIndices[i][2],
             torsionIndices[i][3], torsionIndices[i][4], torsionIndices[i][5], torsionIndices[i][6], torsionIndices[i][7]);
     }
+    usePeriodic = force.usesPeriodicBoundaryConditions();
 }
 
 double ReferenceCalcCMAPTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
@@ -729,6 +748,8 @@ double ReferenceCalcCMAPTorsionForceKernel::execute(ContextImpl& context, bool i
     vector<RealVec>& forceData = extractForces(context);
     RealOpenMM totalEnergy = 0;
     ReferenceCMAPTorsionIxn torsion(coeff, torsionMaps, torsionIndices);
+    if (usePeriodic)
+        torsion.setPeriodic(extractBoxVectors(context));
     torsion.calculateIxn(posData, forceData, &totalEnergy);
     return totalEnergy;
 }
@@ -775,6 +796,7 @@ ReferenceCalcCustomTorsionForceKernel::~ReferenceCalcCustomTorsionForceKernel()
 void ReferenceCalcCustomTorsionForceKernel::initialize(const System& system, const CustomTorsionForce& force) {
     numTorsions = force.getNumTorsions();
     int numParameters = force.getNumPerTorsionParameters();
+    usePeriodic = force.usesPeriodicBoundaryConditions();
 
     // Build the arrays.
 
@@ -817,6 +839,8 @@ double ReferenceCalcCustomTorsionForceKernel::execute(ContextImpl& context, bool
         globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
     ReferenceBondForce refBondForce;
     ReferenceCustomTorsionIxn customTorsion(energyExpression, forceExpression, parameterNames, globalParameters);
+    if (usePeriodic)
+        customTorsion.setPeriodic(extractBoxVectors(context));
     refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, includeEnergy ? &energy : NULL, customTorsion);
     return energy;
 }
@@ -1644,6 +1668,7 @@ ReferenceCalcCustomCentroidBondForceKernel::~ReferenceCalcCustomCentroidBondForc
 }
 
 void ReferenceCalcCustomCentroidBondForceKernel::initialize(const System& system, const CustomCentroidBondForce& force) {
+    usePeriodic = force.usesPeriodicBoundaryConditions();
 
     // Build the arrays.
 
@@ -1697,6 +1722,8 @@ double ReferenceCalcCustomCentroidBondForceKernel::execute(ContextImpl& context,
     map<string, double> globalParameters;
     for (int i = 0; i < (int) globalParameterNames.size(); i++)
         globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
+    if (usePeriodic)
+        ixn->setPeriodic(extractBoxVectors(context));
     ixn->calculatePairIxn(posData, bondParamArray, globalParameters, forceData, includeEnergy ? &energy : NULL);
     return energy;
 }
@@ -1728,6 +1755,7 @@ ReferenceCalcCustomCompoundBondForceKernel::~ReferenceCalcCustomCompoundBondForc
 }
 
 void ReferenceCalcCustomCompoundBondForceKernel::initialize(const System& system, const CustomCompoundBondForce& force) {
+    usePeriodic = force.usesPeriodicBoundaryConditions();
 
     // Build the arrays.
 
@@ -1774,6 +1802,8 @@ double ReferenceCalcCustomCompoundBondForceKernel::execute(ContextImpl& context,
     map<string, double> globalParameters;
     for (int i = 0; i < (int) globalParameterNames.size(); i++)
         globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
+    if (usePeriodic)
+        ixn->setPeriodic(extractBoxVectors(context));
     ixn->calculatePairIxn(posData, bondParamArray, globalParameters, forceData, includeEnergy ? &energy : NULL);
     return energy;
 }

platforms/reference/src/ReferenceTabulatedFunction.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2014 Stanford University and the Authors.           *
+ * Portions copyright (c) 2014-2016 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -80,6 +80,13 @@ ReferenceContinuous1DFunction::ReferenceContinuous1DFunction(const Continuous1DF
     SplineFitter::createNaturalSpline(x, values, derivs);
 }
 
+ReferenceContinuous1DFunction::ReferenceContinuous1DFunction(const ReferenceContinuous1DFunction& other) : function(other.function) {
+    function.getFunctionParameters(values, min, max);
+    x = other.x;
+    values = other.values;
+    derivs = other.derivs;
+}
+
 int ReferenceContinuous1DFunction::getNumArguments() const {
     return 1;
 }
@@ -99,7 +106,7 @@ double ReferenceContinuous1DFunction::evaluateDerivative(const double* arguments
 }
 
 CustomFunction* ReferenceContinuous1DFunction::clone() const {
-    return new ReferenceContinuous1DFunction(function);
+    return new ReferenceContinuous1DFunction(*this);
 }
 
 ReferenceContinuous2DFunction::ReferenceContinuous2DFunction(const Continuous2DFunction& function) : function(function) {
@@ -113,6 +120,14 @@ ReferenceContinuous2DFunction::ReferenceContinuous2DFunction(const Continuous2DF
     SplineFitter::create2DNaturalSpline(x, y, values, c);
 }
 
+ReferenceContinuous2DFunction::ReferenceContinuous2DFunction(const ReferenceContinuous2DFunction& other) : function(other.function) {
+    function.getFunctionParameters(xsize, ysize, values, xmin, xmax, ymin, ymax);
+    x = other.x;
+    y = other.y;
+    values = other.values;
+    c = other.c;
+}
+
 int ReferenceContinuous2DFunction::getNumArguments() const {
     return 2;
 }
@@ -144,7 +159,7 @@ double ReferenceContinuous2DFunction::evaluateDerivative(const double* arguments
 }
 
 CustomFunction* ReferenceContinuous2DFunction::clone() const {
-    return new ReferenceContinuous2DFunction(function);
+    return new ReferenceContinuous2DFunction(*this);
 }
 
 ReferenceContinuous3DFunction::ReferenceContinuous3DFunction(const Continuous3DFunction& function) : function(function) {
@@ -161,6 +176,15 @@ ReferenceContinuous3DFunction::ReferenceContinuous3DFunction(const Continuous3DF
     SplineFitter::create3DNaturalSpline(x, y, z, values, c);
 }
 
+ReferenceContinuous3DFunction::ReferenceContinuous3DFunction(const ReferenceContinuous3DFunction& other) : function(other.function) {
+    function.getFunctionParameters(xsize, ysize, zsize, values, xmin, xmax, ymin, ymax, zmin, zmax);
+    x = other.x;
+    y = other.y;
+    z = other.z;
+    values = other.values;
+    c = other.c;
+}
+
 int ReferenceContinuous3DFunction::getNumArguments() const {
     return 3;
 }
@@ -200,7 +224,7 @@ double ReferenceContinuous3DFunction::evaluateDerivative(const double* arguments
 }
 
 CustomFunction* ReferenceContinuous3DFunction::clone() const {
-    return new ReferenceContinuous3DFunction(function);
+    return new ReferenceContinuous3DFunction(*this);
 }
 
 ReferenceDiscrete1DFunction::ReferenceDiscrete1DFunction(const Discrete1DFunction& function) : function(function) {

platforms/reference/src/SimTKReference/ReferenceAngleBondIxn.cpp

@@ -38,7 +38,7 @@ using namespace OpenMM;
 
    --------------------------------------------------------------------------------------- */
 
-ReferenceAngleBondIxn::ReferenceAngleBondIxn() {
+ReferenceAngleBondIxn::ReferenceAngleBondIxn() : usePeriodic(false) {
 
    // ---------------------------------------------------------------------------------------
 
@@ -64,6 +64,13 @@ ReferenceAngleBondIxn::~ReferenceAngleBondIxn() {
 
 }
 
+void ReferenceAngleBondIxn::setPeriodic(OpenMM::RealVec* vectors) {
+    usePeriodic = true;
+    boxVectors[0] = vectors[0];
+    boxVectors[1] = vectors[1];
+    boxVectors[2] = vectors[2];
+}
+
 /**---------------------------------------------------------------------------------------
 
    Get dEdR and energy term for angle bond
@@ -145,8 +152,14 @@ void ReferenceAngleBondIxn::calculateBondIxn(int* atomIndices,
    int atomAIndex = atomIndices[0];
    int atomBIndex = atomIndices[1];
    int atomCIndex = atomIndices[2];
-   ReferenceForce::getDeltaR(atomCoordinates[atomAIndex], atomCoordinates[atomBIndex], deltaR[0]);  
-   ReferenceForce::getDeltaR(atomCoordinates[atomCIndex], atomCoordinates[atomBIndex], deltaR[1]);  
+   if (usePeriodic) {
+      ReferenceForce::getDeltaRPeriodic(atomCoordinates[atomAIndex], atomCoordinates[atomBIndex], boxVectors, deltaR[0]);  
+      ReferenceForce::getDeltaRPeriodic(atomCoordinates[atomCIndex], atomCoordinates[atomBIndex], boxVectors, deltaR[1]);  
+   }
+   else {
+      ReferenceForce::getDeltaR(atomCoordinates[atomAIndex], atomCoordinates[atomBIndex], deltaR[0]);  
+      ReferenceForce::getDeltaR(atomCoordinates[atomCIndex], atomCoordinates[atomBIndex], deltaR[1]);  
+   }
 
    RealOpenMM pVector[threeI];
    SimTKOpenMMUtilities::crossProductVector3(deltaR[0], deltaR[1], pVector);

platforms/reference/src/SimTKReference/ReferenceCMAPTorsionIxn.cpp

 
-/* Portions copyright (c) 2010 Stanford University and Simbios.
+/* Portions copyright (c) 2010-2016 Stanford University and Simbios.
  * Contributors: Peter Eastman
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -36,7 +36,14 @@ using namespace OpenMM;
 
 ReferenceCMAPTorsionIxn::ReferenceCMAPTorsionIxn(const vector<vector<vector<RealOpenMM> > >& coeff,
         const vector<int>& torsionMaps, const vector<vector<int> >& torsionIndices) :
-        coeff(coeff), torsionMaps(torsionMaps), torsionIndices(torsionIndices) {
+        coeff(coeff), torsionMaps(torsionMaps), torsionIndices(torsionIndices), usePeriodic(false) {
+}
+
+void ReferenceCMAPTorsionIxn::setPeriodic(OpenMM::RealVec* vectors) {
+    usePeriodic = true;
+    boxVectors[0] = vectors[0];
+    boxVectors[1] = vectors[1];
+    boxVectors[2] = vectors[2];
 }
 
 /**---------------------------------------------------------------------------------------
@@ -83,13 +90,23 @@ void ReferenceCMAPTorsionIxn::calculateOneIxn(int index, vector<RealVec>& atomCo
     // Compute deltas between the various atoms involved.
 
     RealOpenMM deltaA[3][ReferenceForce::LastDeltaRIndex];
-    ReferenceForce::getDeltaR(atomCoordinates[a2], atomCoordinates[a1], deltaA[0]);
-    ReferenceForce::getDeltaR(atomCoordinates[a2], atomCoordinates[a3], deltaA[1]);
-    ReferenceForce::getDeltaR(atomCoordinates[a4], atomCoordinates[a3], deltaA[2]);
     RealOpenMM deltaB[3][ReferenceForce::LastDeltaRIndex];
-    ReferenceForce::getDeltaR(atomCoordinates[b2], atomCoordinates[b1], deltaB[0]);
-    ReferenceForce::getDeltaR(atomCoordinates[b2], atomCoordinates[b3], deltaB[1]);
-    ReferenceForce::getDeltaR(atomCoordinates[b4], atomCoordinates[b3], deltaB[2]);
+    if (usePeriodic) {
+        ReferenceForce::getDeltaRPeriodic(atomCoordinates[a2], atomCoordinates[a1], boxVectors, deltaA[0]);
+        ReferenceForce::getDeltaRPeriodic(atomCoordinates[a2], atomCoordinates[a3], boxVectors, deltaA[1]);
+        ReferenceForce::getDeltaRPeriodic(atomCoordinates[a4], atomCoordinates[a3], boxVectors, deltaA[2]);
+        ReferenceForce::getDeltaRPeriodic(atomCoordinates[b2], atomCoordinates[b1], boxVectors, deltaB[0]);
+        ReferenceForce::getDeltaRPeriodic(atomCoordinates[b2], atomCoordinates[b3], boxVectors, deltaB[1]);
+        ReferenceForce::getDeltaRPeriodic(atomCoordinates[b4], atomCoordinates[b3], boxVectors, deltaB[2]);
+    }
+    else {
+        ReferenceForce::getDeltaR(atomCoordinates[a2], atomCoordinates[a1], deltaA[0]);
+        ReferenceForce::getDeltaR(atomCoordinates[a2], atomCoordinates[a3], deltaA[1]);
+        ReferenceForce::getDeltaR(atomCoordinates[a4], atomCoordinates[a3], deltaA[2]);
+        ReferenceForce::getDeltaR(atomCoordinates[b2], atomCoordinates[b1], deltaB[0]);
+        ReferenceForce::getDeltaR(atomCoordinates[b2], atomCoordinates[b3], deltaB[1]);
+        ReferenceForce::getDeltaR(atomCoordinates[b4], atomCoordinates[b3], deltaB[2]);
+    }
 
     // Visual Studio complains if crossProduct declared as 'crossProduct[2][3]'
 

platforms/reference/src/SimTKReference/ReferenceCustomAngleIxn.cpp

-/* Portions copyright (c) 2010-2013 Stanford University and Simbios.
+/* Portions copyright (c) 2010-2016 Stanford University and Simbios.
  * Contributors: Peter Eastman
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -39,7 +39,7 @@ using namespace std;
 
 ReferenceCustomAngleIxn::ReferenceCustomAngleIxn(const Lepton::CompiledExpression& energyExpression,
         const Lepton::CompiledExpression& forceExpression, const vector<string>& parameterNames, map<string, double> globalParameters) :
-        energyExpression(energyExpression), forceExpression(forceExpression) {
+        energyExpression(energyExpression), forceExpression(forceExpression), usePeriodic(false) {
     
     energyTheta = ReferenceForce::getVariablePointer(this->energyExpression, "theta");
     forceTheta = ReferenceForce::getVariablePointer(this->forceExpression, "theta");
@@ -61,13 +61,13 @@ ReferenceCustomAngleIxn::ReferenceCustomAngleIxn(const Lepton::CompiledExpressio
    --------------------------------------------------------------------------------------- */
 
 ReferenceCustomAngleIxn::~ReferenceCustomAngleIxn() {
+}
 
-   // ---------------------------------------------------------------------------------------
-
-   // static const char* methodName = "\nReferenceCustomAngleIxn::~ReferenceCustomAngleIxn";
-
-   // ---------------------------------------------------------------------------------------
-
+void ReferenceCustomAngleIxn::setPeriodic(OpenMM::RealVec* vectors) {
+    usePeriodic = true;
+    boxVectors[0] = vectors[0];
+    boxVectors[1] = vectors[1];
+    boxVectors[2] = vectors[2];
 }
 
 /**---------------------------------------------------------------------------------------
@@ -106,8 +106,14 @@ void ReferenceCustomAngleIxn::calculateBondIxn(int* atomIndices,
    int atomAIndex = atomIndices[0];
    int atomBIndex = atomIndices[1];
    int atomCIndex = atomIndices[2];
-   ReferenceForce::getDeltaR(atomCoordinates[atomAIndex], atomCoordinates[atomBIndex], deltaR[0]);
-   ReferenceForce::getDeltaR(atomCoordinates[atomCIndex], atomCoordinates[atomBIndex], deltaR[1]);
+   if (usePeriodic) {
+      ReferenceForce::getDeltaRPeriodic(atomCoordinates[atomAIndex], atomCoordinates[atomBIndex], boxVectors, deltaR[0]);
+      ReferenceForce::getDeltaRPeriodic(atomCoordinates[atomCIndex], atomCoordinates[atomBIndex], boxVectors, deltaR[1]);
+   }
+   else {
+      ReferenceForce::getDeltaR(atomCoordinates[atomAIndex], atomCoordinates[atomBIndex], deltaR[0]);
+      ReferenceForce::getDeltaR(atomCoordinates[atomCIndex], atomCoordinates[atomBIndex], deltaR[1]);
+   }
    RealOpenMM pVector[3];
    SimTKOpenMMUtilities::crossProductVector3(deltaR[0], deltaR[1], pVector);
    RealOpenMM rp = SQRT(DOT3(pVector, pVector));

platforms/reference/src/SimTKReference/ReferenceCustomBondIxn.cpp

 
-/* Portions copyright (c) 2009-2013 Stanford University and Simbios.
+/* Portions copyright (c) 2009-2016 Stanford University and Simbios.
  * Contributors: Peter Eastman
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -40,7 +40,7 @@ using namespace OpenMM;
 
 ReferenceCustomBondIxn::ReferenceCustomBondIxn(const Lepton::CompiledExpression& energyExpression,
         const Lepton::CompiledExpression& forceExpression, const vector<string>& parameterNames, map<string, double> globalParameters) :
-        energyExpression(energyExpression), forceExpression(forceExpression) {
+        energyExpression(energyExpression), forceExpression(forceExpression), usePeriodic(false) {
     energyR = ReferenceForce::getVariablePointer(this->energyExpression, "r");
     forceR = ReferenceForce::getVariablePointer(this->forceExpression, "r");
     numParameters = parameterNames.size();
@@ -61,13 +61,13 @@ ReferenceCustomBondIxn::ReferenceCustomBondIxn(const Lepton::CompiledExpression&
    --------------------------------------------------------------------------------------- */
 
 ReferenceCustomBondIxn::~ReferenceCustomBondIxn() {
+}
 
-   // ---------------------------------------------------------------------------------------
-
-   // static const char* methodName = "\nReferenceCustomBondIxn::~ReferenceCustomBondIxn";
-
-   // ---------------------------------------------------------------------------------------
-
+void ReferenceCustomBondIxn::setPeriodic(OpenMM::RealVec* vectors) {
+    usePeriodic = true;
+    boxVectors[0] = vectors[0];
+    boxVectors[1] = vectors[1];
+    boxVectors[2] = vectors[2];
 }
 
 /**---------------------------------------------------------------------------------------
@@ -108,7 +108,10 @@ void ReferenceCustomBondIxn::calculateBondIxn(int* atomIndices,
 
    int atomAIndex = atomIndices[0];
    int atomBIndex = atomIndices[1];
-   ReferenceForce::getDeltaR(atomCoordinates[atomAIndex], atomCoordinates[atomBIndex], deltaR);
+   if (usePeriodic)
+       ReferenceForce::getDeltaRPeriodic(atomCoordinates[atomAIndex], atomCoordinates[atomBIndex], boxVectors, deltaR);
+   else
+       ReferenceForce::getDeltaR(atomCoordinates[atomAIndex], atomCoordinates[atomBIndex], deltaR);
    
    ReferenceForce::setVariable(energyR, deltaR[ReferenceForce::RIndex]);
    ReferenceForce::setVariable(forceR, deltaR[ReferenceForce::RIndex]);

platforms/reference/src/SimTKReference/ReferenceCustomCentroidBondIxn.cpp

 
-/* Portions copyright (c) 2009-2015 Stanford University and Simbios.
+/* Portions copyright (c) 2009-2016 Stanford University and Simbios.
  * Contributors: Peter Eastman
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -41,7 +41,7 @@ ReferenceCustomCentroidBondIxn::ReferenceCustomCentroidBondIxn(int numGroupsPerB
             const vector<vector<double> >& normalizedWeights, const vector<vector<int> >& bondGroups,
             const Lepton::ParsedExpression& energyExpression, const vector<string>& bondParameterNames,
             const map<string, vector<int> >& distances, const map<string, vector<int> >& angles, const map<string, vector<int> >& dihedrals) :
-            groupAtoms(groupAtoms), normalizedWeights(normalizedWeights), bondGroups(bondGroups), energyExpression(energyExpression.createProgram()), bondParamNames(bondParameterNames) {
+            groupAtoms(groupAtoms), normalizedWeights(normalizedWeights), bondGroups(bondGroups), energyExpression(energyExpression.createProgram()), bondParamNames(bondParameterNames), usePeriodic(false) {
     for (int i = 0; i < numGroupsPerBond; i++) {
         stringstream xname, yname, zname;
         xname << 'x' << (i+1);
@@ -62,6 +62,13 @@ ReferenceCustomCentroidBondIxn::ReferenceCustomCentroidBondIxn(int numGroupsPerB
 ReferenceCustomCentroidBondIxn::~ReferenceCustomCentroidBondIxn() {
 }
 
+void ReferenceCustomCentroidBondIxn::setPeriodic(OpenMM::RealVec* vectors) {
+    usePeriodic = true;
+    boxVectors[0] = vectors[0];
+    boxVectors[1] = vectors[1];
+    boxVectors[2] = vectors[2];
+}
+
 void ReferenceCustomCentroidBondIxn::calculatePairIxn(vector<RealVec>& atomCoordinates, RealOpenMM** bondParameters,
                                              const map<string, double>& globalParameters, vector<RealVec>& forces,
                                              RealOpenMM* totalEnergy) const {
@@ -215,7 +222,10 @@ void ReferenceCustomCentroidBondIxn::calculateOneIxn(int bond, vector<RealVec>&
 }
 
 void ReferenceCustomCentroidBondIxn::computeDelta(int group1, int group2, RealOpenMM* delta, vector<RealVec>& groupCenters) const {
-    ReferenceForce::getDeltaR(groupCenters[group1], groupCenters[group2], delta);
+    if (usePeriodic)
+        ReferenceForce::getDeltaRPeriodic(groupCenters[group1], groupCenters[group2], boxVectors, delta);
+    else
+        ReferenceForce::getDeltaR(groupCenters[group1], groupCenters[group2], delta);
 }
 
 RealOpenMM ReferenceCustomCentroidBondIxn::computeAngle(RealOpenMM* vec1, RealOpenMM* vec2) {

platforms/reference/src/SimTKReference/ReferenceCustomCompoundBondIxn.cpp

 
-/* Portions copyright (c) 2009-2015 Stanford University and Simbios.
+/* Portions copyright (c) 2009-2016 Stanford University and Simbios.
  * Contributors: Peter Eastman
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -46,7 +46,7 @@ using namespace OpenMM;
 ReferenceCustomCompoundBondIxn::ReferenceCustomCompoundBondIxn(int numParticlesPerBond, const vector<vector<int> >& bondAtoms,
             const Lepton::ParsedExpression& energyExpression, const vector<string>& bondParameterNames,
             const map<string, vector<int> >& distances, const map<string, vector<int> >& angles, const map<string, vector<int> >& dihedrals) :
-            bondAtoms(bondAtoms), energyExpression(energyExpression.createProgram()), bondParamNames(bondParameterNames) {
+            bondAtoms(bondAtoms), energyExpression(energyExpression.createProgram()), bondParamNames(bondParameterNames), usePeriodic(false) {
     for (int i = 0; i < numParticlesPerBond; i++) {
         stringstream xname, yname, zname;
         xname << 'x' << (i+1);
@@ -73,6 +73,13 @@ ReferenceCustomCompoundBondIxn::ReferenceCustomCompoundBondIxn(int numParticlesP
 ReferenceCustomCompoundBondIxn::~ReferenceCustomCompoundBondIxn() {
 }
 
+void ReferenceCustomCompoundBondIxn::setPeriodic(OpenMM::RealVec* vectors) {
+    usePeriodic = true;
+    boxVectors[0] = vectors[0];
+    boxVectors[1] = vectors[1];
+    boxVectors[2] = vectors[2];
+}
+
 /**---------------------------------------------------------------------------------------
 
    Calculate custom hbond interaction
@@ -232,7 +239,10 @@ void ReferenceCustomCompoundBondIxn::calculateOneIxn(int bond, vector<RealVec>&
 }
 
 void ReferenceCustomCompoundBondIxn::computeDelta(int atom1, int atom2, RealOpenMM* delta, vector<RealVec>& atomCoordinates) const {
-    ReferenceForce::getDeltaR(atomCoordinates[atom1], atomCoordinates[atom2], delta);
+    if (usePeriodic)
+        ReferenceForce::getDeltaRPeriodic(atomCoordinates[atom1], atomCoordinates[atom2], boxVectors, delta);
+    else
+        ReferenceForce::getDeltaR(atomCoordinates[atom1], atomCoordinates[atom2], delta);
 }
 
 RealOpenMM ReferenceCustomCompoundBondIxn::computeAngle(RealOpenMM* vec1, RealOpenMM* vec2) {

platforms/reference/src/SimTKReference/ReferenceCustomDynamics.cpp

 
-/* Portions copyright (c) 2011-2015 Stanford University and Simbios.
+/* Portions copyright (c) 2011-2016 Stanford University and Simbios.
  * Contributors: Peter Eastman
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -56,13 +56,11 @@ ReferenceCustomDynamics::ReferenceCustomDynamics(int numberOfAtoms, const Custom
         string expression;
         integrator.getComputationStep(i, stepType[i], stepVariable[i], expression);
     }
-    kineticEnergyExpression = Lepton::Parser::parse(integrator.getKineticEnergyExpression()).optimize().createProgram();
+    kineticEnergyExpression = Lepton::Parser::parse(integrator.getKineticEnergyExpression()).optimize().createCompiledExpression();
+    expressionSet.registerExpression(kineticEnergyExpression);
     kineticEnergyNeedsForce = false;
-    for (int i = 0; i < kineticEnergyExpression.getNumOperations(); i++) {
-        const Lepton::Operation& op = kineticEnergyExpression.getOperation(i);
-        if (op.getId() == Lepton::Operation::VARIABLE && op.getName() == "f")
-            kineticEnergyNeedsForce = true;
-    }
+    if (kineticEnergyExpression.getVariables().find("f") != kineticEnergyExpression.getVariables().end())
+        kineticEnergyNeedsForce = true;
 }
 
 /**---------------------------------------------------------------------------------------
@@ -74,6 +72,75 @@ ReferenceCustomDynamics::ReferenceCustomDynamics(int numberOfAtoms, const Custom
 ReferenceCustomDynamics::~ReferenceCustomDynamics() {
 }
 
+void ReferenceCustomDynamics::initialize(ContextImpl& context, vector<RealOpenMM>& masses, map<string, RealOpenMM>& globals) {
+    // Some initialization can't be done in the constructor, since we need a ContextImpl from which to get the list of
+    // Context parameters.  Instead, we do it the first time update() or computeKineticEnergy() is called.
+
+    int numSteps = stepType.size();
+    vector<int> forceGroup;
+    vector<vector<Lepton::ParsedExpression> > expressions;
+    CustomIntegratorUtilities::analyzeComputations(context, integrator, expressions, comparisons, blockEnd, invalidatesForces, needsForces, needsEnergy, computeBothForceAndEnergy, forceGroup);
+    stepExpressions.resize(expressions.size());
+    for (int i = 0; i < numSteps; i++) {
+        stepExpressions[i].resize(expressions[i].size());
+        for (int j = 0; j < (int) expressions[i].size(); j++) {
+            stepExpressions[i][j] = expressions[i][j].createCompiledExpression();
+            expressionSet.registerExpression(stepExpressions[i][j]);
+        }
+        if (stepType[i] == CustomIntegrator::WhileBlockStart)
+            blockEnd[blockEnd[i]] = i; // Record where to branch back to.
+    }
+
+    // Record the variable names and flags for the force and energy in each step.
+
+    forceGroupFlags.resize(numSteps, -1);
+    fIndex = expressionSet.getVariableIndex("f");
+    energyIndex = expressionSet.getVariableIndex("energy");
+    forceVariableIndex.resize(numSteps, fIndex);
+    energyVariableIndex.resize(numSteps, energyIndex);
+    vector<string> forceGroupName;
+    vector<string> energyGroupName;
+    for (int i = 0; i < 32; i++) {
+        stringstream fname;
+        fname << "f" << i;
+        forceGroupName.push_back(fname.str());
+        stringstream ename;
+        ename << "energy" << i;
+        energyGroupName.push_back(ename.str());
+    }
+    for (int i = 0; i < numSteps; i++) {
+        if (needsForces[i] && forceGroup[i] > -1)
+            forceVariableIndex[i] = expressionSet.getVariableIndex(forceGroupName[forceGroup[i]]);
+        if (needsEnergy[i] && forceGroup[i] > -1)
+            energyVariableIndex[i] = expressionSet.getVariableIndex(energyGroupName[forceGroup[i]]);
+        if (forceGroup[i] > -1)
+            forceGroupFlags[i] = 1<<forceGroup[i];
+    }
+
+    // Build the list of inverse masses.
+
+    int numberOfAtoms = masses.size();
+    inverseMasses.resize(numberOfAtoms);
+    for (int i = 0; i < numberOfAtoms; i++) {
+        if (masses[i] == 0.0)
+            inverseMasses[i] = 0.0;
+        else
+            inverseMasses[i] = 1.0/masses[i];
+    }
+
+    // Record indices of other variables.
+
+    xIndex = expressionSet.getVariableIndex("x");
+    vIndex = expressionSet.getVariableIndex("v");
+    mIndex = expressionSet.getVariableIndex("m");
+    gaussianIndex = expressionSet.getVariableIndex("gaussian");
+    uniformIndex = expressionSet.getVariableIndex("uniform");
+    for (int i = 0; i < integrator.getNumPerDofVariables(); i++)
+        perDofVariableIndex.push_back(expressionSet.getVariableIndex(integrator.getPerDofVariableName(i)));
+    for (int i = 0; i < stepVariable.size(); i++)
+        stepVariableIndex.push_back(expressionSet.getVariableIndex(stepVariable[i]));
+}
+
 /**---------------------------------------------------------------------------------------
 
    Update -- driver routine for performing Custom dynamics update of coordinates
@@ -93,58 +160,13 @@ ReferenceCustomDynamics::~ReferenceCustomDynamics() {
 void ReferenceCustomDynamics::update(ContextImpl& context, int numberOfAtoms, vector<RealVec>& atomCoordinates,
                                      vector<RealVec>& velocities, vector<RealVec>& forces, vector<RealOpenMM>& masses,
                                      map<string, RealOpenMM>& globals, vector<vector<RealVec> >& perDof, bool& forcesAreValid, RealOpenMM tolerance) {
+    if (invalidatesForces.size() == 0)
+        initialize(context, masses, globals);
     int numSteps = stepType.size();
     globals.insert(context.getParameters().begin(), context.getParameters().end());
+    for (map<string, RealOpenMM>::const_iterator iter = globals.begin(); iter != globals.end(); ++iter)
+        expressionSet.setVariable(expressionSet.getVariableIndex(iter->first), iter->second);
     oldPos = atomCoordinates;
-    if (invalidatesForces.size() == 0) {
-        // Some initialization can't be done in the constructor, since we need a ContextImpl from which to get the list of
-        // Context parameters.  Instead, we do it the first time this method is called.
-
-        vector<int> forceGroup;
-        vector<vector<Lepton::ParsedExpression> > expressions;
-        CustomIntegratorUtilities::analyzeComputations(context, integrator, expressions, comparisons, blockEnd, invalidatesForces, needsForces, needsEnergy, computeBothForceAndEnergy, forceGroup);
-        stepExpressions.resize(expressions.size());
-        for (int i = 0; i < numSteps; i++) {
-            for (int j = 0; j < (int) expressions[i].size(); j++)
-                stepExpressions[i].push_back(expressions[i][j].createProgram());
-            if (stepType[i] == CustomIntegrator::BeginWhileBlock)
-                blockEnd[blockEnd[i]] = i; // Record where to branch back to.
-        }
-
-        // Record the variable names and flags for the force and energy in each step.
-
-        forceGroupFlags.resize(numSteps, -1);
-        forceName.resize(numSteps, "f");
-        energyName.resize(numSteps, "energy");
-        vector<string> forceGroupName;
-        vector<string> energyGroupName;
-        for (int i = 0; i < 32; i++) {
-            stringstream fname;
-            fname << "f" << i;
-            forceGroupName.push_back(fname.str());
-            stringstream ename;
-            ename << "energy" << i;
-            energyGroupName.push_back(ename.str());
-        }
-        for (int i = 0; i < numSteps; i++) {
-            if (needsForces[i] && forceGroup[i] > -1)
-                forceName[i] = forceGroupName[forceGroup[i]];
-            if (needsEnergy[i] && forceGroup[i] > -1)
-                energyName[i] = energyGroupName[forceGroup[i]];
-            if (forceGroup[i] > -1)
-                forceGroupFlags[i] = 1<<forceGroup[i];
-        }
-        
-        // Build the list of inverse masses.
-        
-        inverseMasses.resize(numberOfAtoms);
-        for (int i = 0; i < numberOfAtoms; i++) {
-            if (masses[i] == 0.0)
-                inverseMasses[i] = 0.0;
-            else
-                inverseMasses[i] = 1.0/masses[i];
-        }
-    }
     
     // Loop over steps and execute them.
     
@@ -160,42 +182,44 @@ void ReferenceCustomDynamics::update(ContextImpl& context, int numberOfAtoms, ve
                 energy = e;
             forcesAreValid = true;
         }
-        globals[energyName[step]] = energy;
+        expressionSet.setVariable(energyVariableIndex[step], energy);
         
         // Execute the step.
 
         int nextStep = step+1;
         switch (stepType[step]) {
             case CustomIntegrator::ComputeGlobal: {
-                map<string, RealOpenMM> variables = globals;
-                variables["uniform"] = SimTKOpenMMUtilities::getUniformlyDistributedRandomNumber();
-                variables["gaussian"] = SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
-                globals[stepVariable[step]] = stepExpressions[step][0].evaluate(variables);
+                expressionSet.setVariable(uniformIndex, SimTKOpenMMUtilities::getUniformlyDistributedRandomNumber());
+                expressionSet.setVariable(gaussianIndex, SimTKOpenMMUtilities::getNormallyDistributedRandomNumber());
+                RealOpenMM result = stepExpressions[step][0].evaluate();
+                globals[stepVariable[step]] = result;
+                expressionSet.setVariable(stepVariableIndex[step], result);
                 break;
             }
             case CustomIntegrator::ComputePerDof: {
                 vector<RealVec>* results = NULL;
-                if (stepVariable[step] == "x")
+                if (stepVariableIndex[step] == xIndex)
                     results = &atomCoordinates;
-                else if (stepVariable[step] == "v")
+                else if (stepVariableIndex[step] == vIndex)
                     results = &velocities;
                 else {
                     for (int j = 0; j < integrator.getNumPerDofVariables(); j++)
-                        if (stepVariable[step] == integrator.getPerDofVariableName(j))
+                        if (stepVariableIndex[step] == perDofVariableIndex[j])
                             results = &perDof[j];
                 }
                 if (results == NULL)
                     throw OpenMMException("Illegal per-DOF output variable: "+stepVariable[step]);
-                computePerDof(numberOfAtoms, *results, atomCoordinates, velocities, forces, masses, globals, perDof, stepExpressions[step][0], forceName[step]);
+                computePerDof(numberOfAtoms, *results, atomCoordinates, velocities, forces, masses, perDof, stepExpressions[step][0], forceVariableIndex[step]);
                 break;
             }
             case CustomIntegrator::ComputeSum: {
-                computePerDof(numberOfAtoms, sumBuffer, atomCoordinates, velocities, forces, masses, globals, perDof, stepExpressions[step][0], forceName[step]);
+                computePerDof(numberOfAtoms, sumBuffer, atomCoordinates, velocities, forces, masses, perDof, stepExpressions[step][0], forceVariableIndex[step]);
                 RealOpenMM sum = 0.0;
                 for (int j = 0; j < numberOfAtoms; j++)
                     if (masses[j] != 0.0)
                         sum += sumBuffer[j][0]+sumBuffer[j][1]+sumBuffer[j][2];
                 globals[stepVariable[step]] = sum;
+                expressionSet.setVariable(stepVariableIndex[step], sum);
                 break;
             }
             case CustomIntegrator::ConstrainPositions: {
@@ -211,19 +235,21 @@ void ReferenceCustomDynamics::update(ContextImpl& context, int numberOfAtoms, ve
                 recordChangedParameters(context, globals);
                 context.updateContextState();
                 globals.insert(context.getParameters().begin(), context.getParameters().end());
+                for (map<string, RealOpenMM>::const_iterator iter = globals.begin(); iter != globals.end(); ++iter)
+                    expressionSet.setVariable(expressionSet.getVariableIndex(iter->first), iter->second);
                 break;
             }
-            case CustomIntegrator::BeginIfBlock: {
-                if (!evaluateCondition(step, globals))
+            case CustomIntegrator::IfBlockStart: {
+                if (!evaluateCondition(step))
                     nextStep = blockEnd[step]+1;
                 break;
             }
-            case CustomIntegrator::BeginWhileBlock: {
-                if (!evaluateCondition(step, globals))
+            case CustomIntegrator::WhileBlockStart: {
+                if (!evaluateCondition(step))
                     nextStep = blockEnd[step]+1;
                 break;
             }
-            case CustomIntegrator::EndBlock: {
+            case CustomIntegrator::BlockEnd: {
                 if (blockEnd[step] != -1)
                     nextStep = blockEnd[step]; // Return to the start of a while block.
                 break;
@@ -240,36 +266,33 @@ void ReferenceCustomDynamics::update(ContextImpl& context, int numberOfAtoms, ve
 
 void ReferenceCustomDynamics::computePerDof(int numberOfAtoms, vector<RealVec>& results, const vector<RealVec>& atomCoordinates,
               const vector<RealVec>& velocities, const vector<RealVec>& forces, const vector<RealOpenMM>& masses,
-              const map<string, RealOpenMM>& globals, const vector<vector<RealVec> >& perDof,
-              const Lepton::ExpressionProgram& expression, const std::string& forceName) {
+              const vector<vector<RealVec> >& perDof, const Lepton::CompiledExpression& expression, int forceIndex) {
     // Loop over all degrees of freedom.
 
-    map<string, RealOpenMM> variables = globals;
     for (int i = 0; i < numberOfAtoms; i++) {
         if (masses[i] != 0.0) {
-            variables["m"] = masses[i];
+            expressionSet.setVariable(mIndex, masses[i]);
             for (int j = 0; j < 3; j++) {
                 // Compute the expression.
 
-                variables["x"] = atomCoordinates[i][j];
-                variables["v"] = velocities[i][j];
-                variables[forceName] = forces[i][j];
-                variables["uniform"] = SimTKOpenMMUtilities::getUniformlyDistributedRandomNumber();
-                variables["gaussian"] = SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
+                expressionSet.setVariable(xIndex, atomCoordinates[i][j]);
+                expressionSet.setVariable(vIndex, velocities[i][j]);
+                expressionSet.setVariable(forceIndex, forces[i][j]);
+                expressionSet.setVariable(uniformIndex, SimTKOpenMMUtilities::getUniformlyDistributedRandomNumber());
+                expressionSet.setVariable(gaussianIndex, SimTKOpenMMUtilities::getNormallyDistributedRandomNumber());
                 for (int k = 0; k < (int) perDof.size(); k++)
-                    variables[integrator.getPerDofVariableName(k)] = perDof[k][i][j];
-                results[i][j] = expression.evaluate(variables);
+                    expressionSet.setVariable(perDofVariableIndex[k], perDof[k][i][j]);
+                results[i][j] = expression.evaluate();
             }
         }
     }
 }
 
-bool ReferenceCustomDynamics::evaluateCondition(int step, map<string, RealOpenMM>& globals) {
-    map<string, RealOpenMM> variables = globals;
-    variables["uniform"] = SimTKOpenMMUtilities::getUniformlyDistributedRandomNumber();
-    variables["gaussian"] = SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
-    double lhs = stepExpressions[step][0].evaluate(variables);
-    double rhs = stepExpressions[step][1].evaluate(variables);
+bool ReferenceCustomDynamics::evaluateCondition(int step) {
+    expressionSet.setVariable(uniformIndex, SimTKOpenMMUtilities::getUniformlyDistributedRandomNumber());
+    expressionSet.setVariable(gaussianIndex, SimTKOpenMMUtilities::getNormallyDistributedRandomNumber());
+    double lhs = stepExpressions[step][0].evaluate();
+    double rhs = stepExpressions[step][1].evaluate();
     switch (comparisons[step]) {
         case CustomIntegratorUtilities::EQUAL:
             return (lhs == rhs);
@@ -318,12 +341,16 @@ void ReferenceCustomDynamics::recordChangedParameters(OpenMM::ContextImpl& conte
 double ReferenceCustomDynamics::computeKineticEnergy(OpenMM::ContextImpl& context, int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
         std::vector<OpenMM::RealVec>& velocities, std::vector<OpenMM::RealVec>& forces, std::vector<RealOpenMM>& masses,
         std::map<std::string, RealOpenMM>& globals, std::vector<std::vector<OpenMM::RealVec> >& perDof, bool& forcesAreValid) {
+    if (invalidatesForces.size() == 0)
+        initialize(context, masses, globals);
     globals.insert(context.getParameters().begin(), context.getParameters().end());
+    for (map<string, RealOpenMM>::const_iterator iter = globals.begin(); iter != globals.end(); ++iter)
+        expressionSet.setVariable(expressionSet.getVariableIndex(iter->first), iter->second);
     if (kineticEnergyNeedsForce) {
         energy = context.calcForcesAndEnergy(true, true, -1);
         forcesAreValid = true;
     }
-    computePerDof(numberOfAtoms, sumBuffer, atomCoordinates, velocities, forces, masses, globals, perDof, kineticEnergyExpression, "f");
+    computePerDof(numberOfAtoms, sumBuffer, atomCoordinates, velocities, forces, masses, perDof, kineticEnergyExpression, fIndex);
     RealOpenMM sum = 0.0;
     for (int j = 0; j < numberOfAtoms; j++)
         if (masses[j] != 0.0)

platforms/reference/src/SimTKReference/ReferenceCustomTorsionIxn.cpp

-/* Portions copyright (c) 2010-2013 Stanford University and Simbios.
+/* Portions copyright (c) 2010-2016 Stanford University and Simbios.
  * Contributors: Peter Eastman
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -39,7 +39,7 @@ using namespace OpenMM;
 
 ReferenceCustomTorsionIxn::ReferenceCustomTorsionIxn(const Lepton::CompiledExpression& energyExpression,
         const Lepton::CompiledExpression& forceExpression, const vector<string>& parameterNames, map<string, double> globalParameters) :
-        energyExpression(energyExpression), forceExpression(forceExpression) {
+        energyExpression(energyExpression), forceExpression(forceExpression), usePeriodic(false) {
 
     energyTheta = ReferenceForce::getVariablePointer(this->energyExpression, "theta");
     forceTheta = ReferenceForce::getVariablePointer(this->forceExpression, "theta");
@@ -61,13 +61,13 @@ ReferenceCustomTorsionIxn::ReferenceCustomTorsionIxn(const Lepton::CompiledExpre
    --------------------------------------------------------------------------------------- */
 
 ReferenceCustomTorsionIxn::~ReferenceCustomTorsionIxn() {
+}
 
-   // ---------------------------------------------------------------------------------------
-
-   // static const char* methodName = "\nReferenceCustomTorsionIxn::~ReferenceCustomTorsionIxn";
-
-   // ---------------------------------------------------------------------------------------
-
+void ReferenceCustomTorsionIxn::setPeriodic(OpenMM::RealVec* vectors) {
+    usePeriodic = true;
+    boxVectors[0] = vectors[0];
+    boxVectors[1] = vectors[1];
+    boxVectors[2] = vectors[2];
 }
 
 /**---------------------------------------------------------------------------------------
@@ -107,9 +107,16 @@ void ReferenceCustomTorsionIxn::calculateBondIxn(int* atomIndices,
    int atomBIndex = atomIndices[1];
    int atomCIndex = atomIndices[2];
    int atomDIndex = atomIndices[3];
-   ReferenceForce::getDeltaR(atomCoordinates[atomBIndex], atomCoordinates[atomAIndex], deltaR[0]);
-   ReferenceForce::getDeltaR(atomCoordinates[atomBIndex], atomCoordinates[atomCIndex], deltaR[1]);
-   ReferenceForce::getDeltaR(atomCoordinates[atomDIndex], atomCoordinates[atomCIndex], deltaR[2]);
+   if (usePeriodic) {
+      ReferenceForce::getDeltaRPeriodic(atomCoordinates[atomBIndex], atomCoordinates[atomAIndex], boxVectors, deltaR[0]);  
+      ReferenceForce::getDeltaRPeriodic(atomCoordinates[atomBIndex], atomCoordinates[atomCIndex], boxVectors, deltaR[1]);  
+      ReferenceForce::getDeltaRPeriodic(atomCoordinates[atomDIndex], atomCoordinates[atomCIndex], boxVectors, deltaR[2]);  
+   }
+   else {
+      ReferenceForce::getDeltaR(atomCoordinates[atomBIndex], atomCoordinates[atomAIndex], deltaR[0]);  
+      ReferenceForce::getDeltaR(atomCoordinates[atomBIndex], atomCoordinates[atomCIndex], deltaR[1]);  
+      ReferenceForce::getDeltaR(atomCoordinates[atomDIndex], atomCoordinates[atomCIndex], deltaR[2]);  
+   }
 
    // Visual Studio complains if crossProduct declared as 'crossProduct[2][3]'
 

platforms/reference/src/SimTKReference/ReferenceForce.cpp

@@ -73,6 +73,10 @@ void ReferenceForce::getDeltaR(const RealVec& atomCoordinatesI, const RealVec& a
    deltaR[RIndex]    = (RealOpenMM) SQRT(deltaR[R2Index]);
 }
 
+RealVec ReferenceForce::getDeltaR(const RealVec& atomCoordinatesI, const RealVec& atomCoordinatesJ) {
+    return atomCoordinatesJ-atomCoordinatesI;
+}
+
 void ReferenceForce::getDeltaRPeriodic(const RealVec& atomCoordinatesI, const RealVec& atomCoordinatesJ,
                                const RealOpenMM* boxSize, RealOpenMM* deltaR) {
    deltaR[XIndex]    = periodicDifference(atomCoordinatesJ[0], atomCoordinatesI[0], boxSize[0]);
@@ -96,6 +100,15 @@ void ReferenceForce::getDeltaRPeriodic(const RealVec& atomCoordinatesI, const Re
     deltaR[RIndex] = SQRT(deltaR[R2Index]);
 }
 
+RealVec ReferenceForce::getDeltaRPeriodic(const RealVec& atomCoordinatesI, const RealVec& atomCoordinatesJ,
+                               const RealVec* boxVectors) {
+    RealVec diff = atomCoordinatesJ-atomCoordinatesI;
+    diff -= boxVectors[2]*floor(diff[2]/boxVectors[2][2]+0.5);
+    diff -= boxVectors[1]*floor(diff[1]/boxVectors[1][1]+0.5);
+    diff -= boxVectors[0]*floor(diff[0]/boxVectors[0][0]+0.5);
+    return diff;
+}
+
 double* ReferenceForce::getVariablePointer(Lepton::CompiledExpression& expression, const std::string& name) {
     if (expression.getVariables().find(name) == expression.getVariables().end())
         return NULL;

platforms/reference/src/SimTKReference/ReferenceHarmonicBondIxn.cpp

 
-/* Portions copyright (c) 2006 Stanford University and Simbios.
+/* Portions copyright (c) 2006-2016 Stanford University and Simbios.
  * Contributors: Pande Group
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -38,14 +38,7 @@ using namespace OpenMM;
 
    --------------------------------------------------------------------------------------- */
 
-ReferenceHarmonicBondIxn::ReferenceHarmonicBondIxn() {
-
-   // ---------------------------------------------------------------------------------------
-
-   // static const char* methodName = "\nReferenceHarmonicBondIxn::ReferenceHarmonicBondIxn";
-
-   // ---------------------------------------------------------------------------------------
-
+ReferenceHarmonicBondIxn::ReferenceHarmonicBondIxn() : usePeriodic(false) {
 }
 
 /**---------------------------------------------------------------------------------------
@@ -55,13 +48,13 @@ ReferenceHarmonicBondIxn::ReferenceHarmonicBondIxn() {
    --------------------------------------------------------------------------------------- */
 
 ReferenceHarmonicBondIxn::~ReferenceHarmonicBondIxn() {
+}
 
-   // ---------------------------------------------------------------------------------------
-
-   // static const char* methodName = "\nReferenceHarmonicBondIxn::~ReferenceHarmonicBondIxn";
-
-   // ---------------------------------------------------------------------------------------
-
+void ReferenceHarmonicBondIxn::setPeriodic(OpenMM::RealVec* vectors) {
+    usePeriodic = true;
+    boxVectors[0] = vectors[0];
+    boxVectors[1] = vectors[1];
+    boxVectors[2] = vectors[2];
 }
 
 /**---------------------------------------------------------------------------------------
@@ -99,7 +92,10 @@ void ReferenceHarmonicBondIxn::calculateBondIxn(int* atomIndices,
 
    int atomAIndex = atomIndices[0];
    int atomBIndex = atomIndices[1];
-   ReferenceForce::getDeltaR(atomCoordinates[atomAIndex], atomCoordinates[atomBIndex], deltaR);  
+   if (usePeriodic)
+       ReferenceForce::getDeltaRPeriodic(atomCoordinates[atomAIndex], atomCoordinates[atomBIndex], boxVectors, deltaR);  
+   else
+       ReferenceForce::getDeltaR(atomCoordinates[atomAIndex], atomCoordinates[atomBIndex], deltaR);  
 
    // deltaIdeal = r - r_0
 

platforms/reference/src/SimTKReference/ReferenceProperDihedralBond.cpp

 
-/* Portions copyright (c) 2006 Stanford University and Simbios.
+/* Portions copyright (c) 2006-2016 Stanford University and Simbios.
  * Contributors: Pande Group
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -38,14 +38,7 @@ using namespace OpenMM;
 
    --------------------------------------------------------------------------------------- */
 
-ReferenceProperDihedralBond::ReferenceProperDihedralBond() {
-
-   // ---------------------------------------------------------------------------------------
-
-   // static const char* methodName = "\nReferenceProperDihedralBond::ReferenceProperDihedralBond";
-
-   // ---------------------------------------------------------------------------------------
-
+ReferenceProperDihedralBond::ReferenceProperDihedralBond() : usePeriodic(false) {
 }
 
 /**---------------------------------------------------------------------------------------
@@ -55,13 +48,13 @@ ReferenceProperDihedralBond::ReferenceProperDihedralBond() {
    --------------------------------------------------------------------------------------- */
 
 ReferenceProperDihedralBond::~ReferenceProperDihedralBond() {
+}
 
-   // ---------------------------------------------------------------------------------------
-
-   // static const char* methodName = "\nReferenceProperDihedralBond::~ReferenceProperDihedralBond";
-
-   // ---------------------------------------------------------------------------------------
-
+void ReferenceProperDihedralBond::setPeriodic(OpenMM::RealVec* vectors) {
+    usePeriodic = true;
+    boxVectors[0] = vectors[0];
+    boxVectors[1] = vectors[1];
+    boxVectors[2] = vectors[2];
 }
 
 /**---------------------------------------------------------------------------------------
@@ -110,9 +103,16 @@ void ReferenceProperDihedralBond::calculateBondIxn(int* atomIndices,
    int atomBIndex = atomIndices[1];
    int atomCIndex = atomIndices[2];
    int atomDIndex = atomIndices[3];
-   ReferenceForce::getDeltaR(atomCoordinates[atomBIndex], atomCoordinates[atomAIndex], deltaR[0]);  
-   ReferenceForce::getDeltaR(atomCoordinates[atomBIndex], atomCoordinates[atomCIndex], deltaR[1]);  
-   ReferenceForce::getDeltaR(atomCoordinates[atomDIndex], atomCoordinates[atomCIndex], deltaR[2]);  
+   if (usePeriodic) {
+      ReferenceForce::getDeltaRPeriodic(atomCoordinates[atomBIndex], atomCoordinates[atomAIndex], boxVectors, deltaR[0]);  
+      ReferenceForce::getDeltaRPeriodic(atomCoordinates[atomBIndex], atomCoordinates[atomCIndex], boxVectors, deltaR[1]);  
+      ReferenceForce::getDeltaRPeriodic(atomCoordinates[atomDIndex], atomCoordinates[atomCIndex], boxVectors, deltaR[2]);  
+   }
+   else {
+      ReferenceForce::getDeltaR(atomCoordinates[atomBIndex], atomCoordinates[atomAIndex], deltaR[0]);  
+      ReferenceForce::getDeltaR(atomCoordinates[atomBIndex], atomCoordinates[atomCIndex], deltaR[1]);  
+      ReferenceForce::getDeltaR(atomCoordinates[atomDIndex], atomCoordinates[atomCIndex], deltaR[2]);  
+   }
 
    RealOpenMM dotDihedral;
    RealOpenMM signOfAngle;

platforms/reference/src/SimTKReference/ReferenceRbDihedralBond.cpp

 
-/* Portions copyright (c) 2006 Stanford University and Simbios.
+/* Portions copyright (c) 2006-2016 Stanford University and Simbios.
  * Contributors: Pande Group
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -38,14 +38,7 @@ using namespace OpenMM;
 
    --------------------------------------------------------------------------------------- */
 
-ReferenceRbDihedralBond::ReferenceRbDihedralBond() {
-
-   // ---------------------------------------------------------------------------------------
-
-   // static const char* methodName = "\nReferenceRbDihedralBond::ReferenceRbDihedralBond";
-
-   // ---------------------------------------------------------------------------------------
-
+ReferenceRbDihedralBond::ReferenceRbDihedralBond() : usePeriodic(false) {
 }
 
 /**---------------------------------------------------------------------------------------
@@ -55,13 +48,13 @@ ReferenceRbDihedralBond::ReferenceRbDihedralBond() {
    --------------------------------------------------------------------------------------- */
 
 ReferenceRbDihedralBond::~ReferenceRbDihedralBond() {
+}
 
-   // ---------------------------------------------------------------------------------------
-
-   // static const char* methodName = "\nReferenceRbDihedralBond::~ReferenceRbDihedralBond";
-
-   // ---------------------------------------------------------------------------------------
-
+void ReferenceRbDihedralBond::setPeriodic(OpenMM::RealVec* vectors) {
+    usePeriodic = true;
+    boxVectors[0] = vectors[0];
+    boxVectors[1] = vectors[1];
+    boxVectors[2] = vectors[2];
 }
 
 /**---------------------------------------------------------------------------------------
@@ -112,9 +105,16 @@ void ReferenceRbDihedralBond::calculateBondIxn(int* atomIndices,
    int atomBIndex = atomIndices[1];
    int atomCIndex = atomIndices[2];
    int atomDIndex = atomIndices[3];
-   ReferenceForce::getDeltaR(atomCoordinates[atomBIndex], atomCoordinates[atomAIndex], deltaR[0]);  
-   ReferenceForce::getDeltaR(atomCoordinates[atomBIndex], atomCoordinates[atomCIndex], deltaR[1]);  
-   ReferenceForce::getDeltaR(atomCoordinates[atomDIndex], atomCoordinates[atomCIndex], deltaR[2]);  
+   if (usePeriodic) {
+      ReferenceForce::getDeltaRPeriodic(atomCoordinates[atomBIndex], atomCoordinates[atomAIndex], boxVectors, deltaR[0]);  
+      ReferenceForce::getDeltaRPeriodic(atomCoordinates[atomBIndex], atomCoordinates[atomCIndex], boxVectors, deltaR[1]);  
+      ReferenceForce::getDeltaRPeriodic(atomCoordinates[atomDIndex], atomCoordinates[atomCIndex], boxVectors, deltaR[2]);  
+   }
+   else {
+      ReferenceForce::getDeltaR(atomCoordinates[atomBIndex], atomCoordinates[atomAIndex], deltaR[0]);  
+      ReferenceForce::getDeltaR(atomCoordinates[atomBIndex], atomCoordinates[atomCIndex], deltaR[1]);  
+      ReferenceForce::getDeltaR(atomCoordinates[atomDIndex], atomCoordinates[atomCIndex], deltaR[2]);  
+   }
 
    RealOpenMM cosPhi;
    RealOpenMM signOfAngle;

platforms/reference/tests/TestReferenceMonteCarloMembraneBarostat.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2014 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2016 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -82,7 +82,7 @@ void testIdealGas(MonteCarloMembraneBarostat::XYMode xymode, MonteCarloMembraneB
     // Test it for three different temperatures.
 
     for (int i = 0; i < 3; i++) {
-        barostat->setTemperature(temp[i]);
+        barostat->setDefaultTemperature(temp[i]);
         LangevinIntegrator integrator(temp[i], 0.1, 0.01);
         Context context(system, integrator, platform);
         context.setPositions(positions);

plugins/amoeba/openmmapi/include/openmm/AmoebaAngleForce.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2012 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2016 Stanford University and the Authors.      *
  * Authors: Mark Friedrichs, Peter Eastman                                    *
  * Contributors:                                                              *
  *                                                                            *
@@ -166,21 +166,25 @@ public:
      * in an angle cannot be changed, nor can new angles be added.
      */
     void updateParametersInContext(Context& context);
+    /**
+     * Set whether this force should apply periodic boundary conditions when calculating displacements.
+     * Usually this is not appropriate for bonded forces, but there are situations when it can be useful.
+     */
+    void setUsesPeriodicBoundaryConditions(bool periodic);
     /**
      * Returns whether or not this force makes use of periodic boundary
      * conditions.
      *
-     * @returns true if nonbondedMethod uses PBC and false otherwise
+     * @returns true if force uses PBC and false otherwise
      */
-    bool usesPeriodicBoundaryConditions() const {
-        return false;
-    }
+    bool usesPeriodicBoundaryConditions() const;
 protected:
     ForceImpl* createImpl() const;
     double _globalCubicK, _globalQuarticK, _globalPenticK, _globalSexticK;
 private:
     class AngleInfo;
     std::vector<AngleInfo> angles;
+    bool usePeriodic;
 };
 
 /**

plugins/amoeba/openmmapi/include/openmm/AmoebaBondForce.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2012 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2016 Stanford University and the Authors.      *
  * Authors: Mark Friedrichs, Peter Eastman                                    *
  * Contributors:                                                              *
  *                                                                            *
@@ -138,21 +138,25 @@ public:
      * in a bond cannot be changed, nor can new bonds be added.
      */
     void updateParametersInContext(Context& context);
+    /**
+     * Set whether this force should apply periodic boundary conditions when calculating displacements.
+     * Usually this is not appropriate for bonded forces, but there are situations when it can be useful.
+     */
+    void setUsesPeriodicBoundaryConditions(bool periodic);
     /**
      * Returns whether or not this force makes use of periodic boundary
      * conditions.
      *
-     * @returns true if nonbondedMethod uses PBC and false otherwise
+     * @returns true if force uses PBC and false otherwise
      */
-    bool usesPeriodicBoundaryConditions() const {
-        return false;
-    }
+    bool usesPeriodicBoundaryConditions() const;
 protected:
     double _globalQuarticK, _globalCubicK;
     ForceImpl* createImpl() const;
 private:
     class BondInfo;
     std::vector<BondInfo> bonds;
+    bool usePeriodic;
 };
 
 /**

plugins/amoeba/openmmapi/include/openmm/AmoebaInPlaneAngleForce.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2012 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2016 Stanford University and the Authors.      *
  * Authors: Mark Friedrichs, Peter Eastman                                    *
  * Contributors:                                                              *
  *                                                                            *
@@ -171,21 +171,25 @@ public:
      * in an angle cannot be changed, nor can new angles be added.
      */
     void updateParametersInContext(Context& context);
+    /**
+     * Set whether this force should apply periodic boundary conditions when calculating displacements.
+     * Usually this is not appropriate for bonded forces, but there are situations when it can be useful.
+     */
+    void setUsesPeriodicBoundaryConditions(bool periodic);
     /**
      * Returns whether or not this force makes use of periodic boundary
      * conditions.
      *
-     * @returns true if nonbondedMethod uses PBC and false otherwise
+     * @returns true if force uses PBC and false otherwise
      */
-    bool usesPeriodicBoundaryConditions() const {
-        return false;
-    }
+    bool usesPeriodicBoundaryConditions() const;
 protected:
     ForceImpl* createImpl() const;
     double _globalCubicK, _globalQuarticK, _globalPenticK, _globalSexticK;
 private:
     class AngleInfo;
     std::vector<AngleInfo> angles;
+    bool usePeriodic;
 };
 
 /**

plugins/amoeba/openmmapi/include/openmm/AmoebaOutOfPlaneBendForce.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2012 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2016 Stanford University and the Authors.      *
  * Authors: Mark Friedrichs, Peter Eastman                                    *
  * Contributors:                                                              *
  *                                                                            *
@@ -163,21 +163,25 @@ public:
      * in a term cannot be changed, nor can new terms be added.
      */
     void updateParametersInContext(Context& context);
+    /**
+     * Set whether this force should apply periodic boundary conditions when calculating displacements.
+     * Usually this is not appropriate for bonded forces, but there are situations when it can be useful.
+     */
+    void setUsesPeriodicBoundaryConditions(bool periodic);
     /**
      * Returns whether or not this force makes use of periodic boundary
      * conditions.
      *
-     * @returns true if nonbondedMethod uses PBC and false otherwise
+     * @returns true if force uses PBC and false otherwise
      */
-    bool usesPeriodicBoundaryConditions() const {
-        return false;
-    }
+    bool usesPeriodicBoundaryConditions() const;
 protected:
     ForceImpl* createImpl() const;
     double _globalCubicK, _globalQuarticK, _globalPenticK, _globalSexticK;
 private:
     class OutOfPlaneBendInfo;
     std::vector<OutOfPlaneBendInfo> outOfPlaneBends;
+    bool usePeriodic;
 };
 
 /**

plugins/amoeba/openmmapi/include/openmm/AmoebaPiTorsionForce.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2012 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2016 Stanford University and the Authors.      *
  * Authors: Mark Friedrichs, Peter Eastman                                    *
  * Contributors:                                                              *
  *                                                                            *
@@ -113,20 +113,24 @@ public:
      * in a torsion cannot be changed, nor can new torsions be added.
      */
     void updateParametersInContext(Context& context);
+    /**
+     * Set whether this force should apply periodic boundary conditions when calculating displacements.
+     * Usually this is not appropriate for bonded forces, but there are situations when it can be useful.
+     */
+    void setUsesPeriodicBoundaryConditions(bool periodic);
     /**
      * Returns whether or not this force makes use of periodic boundary
      * conditions.
      *
-     * @returns true if nonbondedMethod uses PBC and false otherwise
+     * @returns true if force uses PBC and false otherwise
      */
-    bool usesPeriodicBoundaryConditions() const {
-        return false;
-    }
+    bool usesPeriodicBoundaryConditions() const;
 protected:
     ForceImpl* createImpl() const;
 private:
     class PiTorsionInfo;
     std::vector<PiTorsionInfo> piTorsions;
+    bool usePeriodic;
 };
 
 /**