Minor optimizations to custom forces

libraries/lepton/src/ParsedExpression.cpp

@@ -189,6 +189,12 @@ ExpressionTreeNode ParsedExpression::substituteSimplerExpression(const Expressio
                 return ExpressionTreeNode(new Operation::Divide(), children[0], children[1].getChildren()[0]);
             if (children[0].getOperation().getId() == Operation::RECIPROCAL) // (1/a)*b = b/a
                 return ExpressionTreeNode(new Operation::Divide(), children[1], children[0].getChildren()[0]);
+            if (children[0] == children[1])
+                return ExpressionTreeNode(new Operation::Square(), children[0]); // x*x = square(x)
+            if (children[0].getOperation().getId() == Operation::SQUARE && children[0].getChildren()[0] == children[1])
+                return ExpressionTreeNode(new Operation::Cube(), children[1]); // x*x*x = cube(x)
+            if (children[1].getOperation().getId() == Operation::SQUARE && children[1].getChildren()[0] == children[0])
+                return ExpressionTreeNode(new Operation::Cube(), children[0]); // x*x*x = cube(x)
             break;
         }
         case Operation::DIVIDE:

platforms/opencl/src/OpenCLExpressionUtilities.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) 2009 Stanford University and the Authors.           *
+ * Portions copyright (c) 2009-2011 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -48,26 +48,33 @@ string OpenCLExpressionUtilities::intToString(int value) {
 
 string OpenCLExpressionUtilities::createExpressions(const map<string, ParsedExpression>& expressions, const map<string, string>& variables,
         const vector<pair<string, string> >& functions, const string& prefix, const string& functionParams) {
+    vector<pair<ExpressionTreeNode, string> > variableNodes;
+    for (map<string, string>::const_iterator iter = variables.begin(); iter != variables.end(); ++iter)
+        variableNodes.push_back(make_pair(ExpressionTreeNode(new Operation::Variable(iter->first)), iter->second));
+    return createExpressions(expressions, variableNodes, functions, prefix, functionParams);
+}
+
+string OpenCLExpressionUtilities::createExpressions(const map<string, ParsedExpression>& expressions, const vector<pair<ExpressionTreeNode, string> >& variables,
+        const vector<pair<string, string> >& functions, const string& prefix, const string& functionParams) {
     stringstream out;
     vector<ParsedExpression> allExpressions;
     for (map<string, ParsedExpression>::const_iterator iter = expressions.begin(); iter != expressions.end(); ++iter)
         allExpressions.push_back(iter->second);
-    vector<pair<ExpressionTreeNode, string> > temps;
+    vector<pair<ExpressionTreeNode, string> > temps = variables;
     for (map<string, ParsedExpression>::const_iterator iter = expressions.begin(); iter != expressions.end(); ++iter) {
-        processExpression(out, iter->second.getRootNode(), temps, variables, functions, prefix, functionParams, allExpressions);
+        processExpression(out, iter->second.getRootNode(), temps, functions, prefix, functionParams, allExpressions);
         out << iter->first << getTempName(iter->second.getRootNode(), temps) << ";\n";
     }
     return out.str();
 }
 
 void OpenCLExpressionUtilities::processExpression(stringstream& out, const ExpressionTreeNode& node, vector<pair<ExpressionTreeNode, string> >& temps,
-        const map<string, string>& variables, const vector<pair<string, string> >& functions, const string& prefix, const string& functionParams,
-        const vector<ParsedExpression>& allExpressions) {
+        const vector<pair<string, string> >& functions, const string& prefix, const string& functionParams, const vector<ParsedExpression>& allExpressions) {
     for (int i = 0; i < (int) temps.size(); i++)
         if (temps[i].first == node)
             return;
     for (int i = 0; i < (int) node.getChildren().size(); i++)
-        processExpression(out, node.getChildren()[i], temps, variables, functions, prefix, functionParams, allExpressions);
+        processExpression(out, node.getChildren()[i], temps, functions, prefix, functionParams, allExpressions);
     string name = prefix+intToString(temps.size());
     bool hasRecordedNode = false;
     
@@ -77,13 +84,7 @@ void OpenCLExpressionUtilities::processExpression(stringstream& out, const Expre
             out << doubleToString(dynamic_cast<const Operation::Constant*>(&node.getOperation())->getValue());
             break;
         case Operation::VARIABLE:
-        {
-            map<string, string>::const_iterator iter = variables.find(node.getOperation().getName());
-            if (iter == variables.end())
             throw OpenMMException("Unknown variable in expression: "+node.getOperation().getName());
-            out << iter->second;
-            break;
-        }
         case Operation::CUSTOM:
         {
             int i;
@@ -145,8 +146,17 @@ void OpenCLExpressionUtilities::processExpression(stringstream& out, const Expre
             out << getTempName(node.getChildren()[0], temps) << "*" << getTempName(node.getChildren()[1], temps);
             break;
         case Operation::DIVIDE:
+        {
+            bool haveReciprocal = false;
+            for (int i = 0; i < (int) temps.size(); i++)
+                if (temps[i].first.getOperation().getId() == Operation::RECIPROCAL && temps[i].first.getChildren()[0] == node.getChildren()[1]) {
+                    haveReciprocal = true;
+                    out << getTempName(node.getChildren()[0], temps) << "*" << temps[i].second;
+                }
+            if (!haveReciprocal)
                 out << getTempName(node.getChildren()[0], temps) << "/" << getTempName(node.getChildren()[1], temps);
             break;
+        }
         case Operation::POWER:
             out << "pow(" << getTempName(node.getChildren()[0], temps) << ", " << getTempName(node.getChildren()[1], temps) << ")";
             break;

platforms/opencl/src/OpenCLExpressionUtilities.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) 2009 Stanford University and the Authors.           *
+ * Portions copyright (c) 2009-2011 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -49,13 +49,26 @@ public:
      * Generate the source code for calculating a set of expressions.
      *
      * @param expressions    the expressions to generate code for (keys are the variables to store the output values in)
-     * @param variables      defines the source code to generate for each variable that may appear in the expressions
+     * @param variables      defines the source code to generate for each variable that may appear in the expressions.  Keys are
+     *                       variable names, and the values are the code to generate for them.
      * @param functions      defines the variable name for each tabulated function that may appear in the expressions
      * @param prefix         a prefix to put in front of temporary variables
      * @param functionParams the variable name containing the parameters for each tabulated function
      */
     static std::string createExpressions(const std::map<std::string, Lepton::ParsedExpression>& expressions, const std::map<std::string, std::string>& variables,
             const std::vector<std::pair<std::string, std::string> >& functions, const std::string& prefix, const std::string& functionParams);
+    /**
+     * Generate the source code for calculating a set of expressions.
+     *
+     * @param expressions    the expressions to generate code for (keys are the variables to store the output values in)
+     * @param variables      defines the source code to generate for each variable or precomputed sub-expression that may appear in the expressions.
+     *                       Each entry is an ExpressionTreeNode, and the code to generate wherever an identical node appears.
+     * @param functions      defines the variable name for each tabulated function that may appear in the expressions
+     * @param prefix         a prefix to put in front of temporary variables
+     * @param functionParams the variable name containing the parameters for each tabulated function
+     */
+    static std::string createExpressions(const std::map<std::string, Lepton::ParsedExpression>& expressions, const std::vector<std::pair<Lepton::ExpressionTreeNode, std::string> >& variables,
+            const std::vector<std::pair<std::string, std::string> >& functions, const std::string& prefix, const std::string& functionParams);
     /**
      * Calculate the spline coefficients for a tabulated function that appears in expressions.
      *
@@ -76,7 +89,7 @@ public:
     class FunctionPlaceholder;
 private:
     static void processExpression(std::stringstream& out, const Lepton::ExpressionTreeNode& node,
-            std::vector<std::pair<Lepton::ExpressionTreeNode, std::string> >& temps, const std::map<std::string, std::string>& variables,
+            std::vector<std::pair<Lepton::ExpressionTreeNode, std::string> >& temps,
             const std::vector<std::pair<std::string, std::string> >& functions, const std::string& prefix, const std::string& functionParams,
             const std::vector<Lepton::ParsedExpression>& allExpressions);
     static std::string getTempName(const Lepton::ExpressionTreeNode& node, const std::vector<std::pair<Lepton::ExpressionTreeNode, std::string> >& temps);

platforms/opencl/src/OpenCLKernels.cpp

@@ -38,6 +38,7 @@
 #include "OpenCLIntegrationUtilities.h"
 #include "OpenCLNonbondedUtilities.h"
 #include "OpenCLKernelSources.h"
+#include "lepton/ExpressionTreeNode.h"
 #include "lepton/Operation.h"
 #include "lepton/Parser.h"
 #include "lepton/ParsedExpression.h"
@@ -47,6 +48,8 @@
 
 using namespace OpenMM;
 using namespace std;
+using Lepton::ExpressionTreeNode;
+using Lepton::Operation;
 
 static string doubleToString(double value) {
     stringstream s;
@@ -68,6 +71,10 @@ static bool isZeroExpression(const Lepton::ParsedExpression& expression) {
     return (dynamic_cast<const Lepton::Operation::Constant&>(op).getValue() == 0.0);
 }
 
+static pair<ExpressionTreeNode, string> makeVariable(const string& name, const string& value) {
+    return make_pair(ExpressionTreeNode(new Operation::Variable(name)), value);
+}
+
 void OpenCLCalcForcesAndEnergyKernel::initialize(const System& system) {
 }
 
@@ -1368,17 +1375,20 @@ void OpenCLCalcCustomNonbondedForceKernel::initialize(const System& system, cons
 
     // Create the kernels.
 
-    map<string, string> variables;
-    variables["r"] = "r";
+    vector<pair<ExpressionTreeNode, string> > variables;
+    ExpressionTreeNode rnode(new Operation::Variable("r"));
+    variables.push_back(make_pair(rnode, "r"));
+    variables.push_back(make_pair(ExpressionTreeNode(new Operation::Square(), rnode), "r2"));
+    variables.push_back(make_pair(ExpressionTreeNode(new Operation::Reciprocal(), rnode), "invR"));
     for (int i = 0; i < force.getNumPerParticleParameters(); i++) {
         const string& name = force.getPerParticleParameterName(i);
-        variables[name+"1"] = prefix+"params"+params->getParameterSuffix(i, "1");
-        variables[name+"2"] = prefix+"params"+params->getParameterSuffix(i, "2");
+        variables.push_back(makeVariable(name+"1", prefix+"params"+params->getParameterSuffix(i, "1")));
+        variables.push_back(makeVariable(name+"2", prefix+"params"+params->getParameterSuffix(i, "2")));
     }
     for (int i = 0; i < force.getNumGlobalParameters(); i++) {
         const string& name = force.getGlobalParameterName(i);
         string value = "globals["+intToString(i)+"]";
-        variables[name] = prefix+value;
+        variables.push_back(makeVariable(name, prefix+value));
     }
     stringstream compute;
     compute << OpenCLExpressionUtilities::createExpressions(forceExpressions, variables, functionDefinitions, prefix+"temp", prefix+"functionParams");
@@ -1787,20 +1797,23 @@ void OpenCLCalcCustomGBForceKernel::initialize(const System& system, const Custo
     {
         // Create the N2 value kernel.
 
-        map<string, string> variables;
+        vector<pair<ExpressionTreeNode, string> > variables;
         map<string, string> rename;
-        variables["r"] = "r";
+        ExpressionTreeNode rnode(new Operation::Variable("r"));
+        variables.push_back(make_pair(rnode, "r"));
+        variables.push_back(make_pair(ExpressionTreeNode(new Operation::Square(), rnode), "r2"));
+        variables.push_back(make_pair(ExpressionTreeNode(new Operation::Reciprocal(), rnode), "invR"));
         for (int i = 0; i < force.getNumPerParticleParameters(); i++) {
             const string& name = force.getPerParticleParameterName(i);
-            variables[name+"1"] = "params"+params->getParameterSuffix(i, "1");
-            variables[name+"2"] = "params"+params->getParameterSuffix(i, "2");
+            variables.push_back(makeVariable(name+"1", "params"+params->getParameterSuffix(i, "1")));
+            variables.push_back(makeVariable(name+"2", "params"+params->getParameterSuffix(i, "2")));
             rename[name+"1"] = name+"2";
             rename[name+"2"] = name+"1";
         }
         for (int i = 0; i < force.getNumGlobalParameters(); i++) {
             const string& name = force.getGlobalParameterName(i);
             string value = "globals["+intToString(i)+"]";
-            variables[name] = value;
+            variables.push_back(makeVariable(name, value));
         }
         map<string, Lepton::ParsedExpression> n2ValueExpressions;
         stringstream n2ValueSource;
@@ -1901,19 +1914,22 @@ void OpenCLCalcCustomGBForceKernel::initialize(const System& system, const Custo
     {
         // Create the N2 energy kernel.
 
-        map<string, string> variables;
-        variables["r"] = "r";
+        vector<pair<ExpressionTreeNode, string> > variables;
+        ExpressionTreeNode rnode(new Operation::Variable("r"));
+        variables.push_back(make_pair(rnode, "r"));
+        variables.push_back(make_pair(ExpressionTreeNode(new Operation::Square(), rnode), "r2"));
+        variables.push_back(make_pair(ExpressionTreeNode(new Operation::Reciprocal(), rnode), "invR"));
         for (int i = 0; i < force.getNumPerParticleParameters(); i++) {
             const string& name = force.getPerParticleParameterName(i);
-            variables[name+"1"] = "params"+params->getParameterSuffix(i, "1");
-            variables[name+"2"] = "params"+params->getParameterSuffix(i, "2");
+            variables.push_back(makeVariable(name+"1", "params"+params->getParameterSuffix(i, "1")));
+            variables.push_back(makeVariable(name+"2", "params"+params->getParameterSuffix(i, "2")));
         }
         for (int i = 0; i < force.getNumComputedValues(); i++) {
-            variables[computedValueNames[i]+"1"] = "values"+computedValues->getParameterSuffix(i, "1");
-            variables[computedValueNames[i]+"2"] = "values"+computedValues->getParameterSuffix(i, "2");
+            variables.push_back(makeVariable(computedValueNames[i]+"1", "values"+computedValues->getParameterSuffix(i, "1")));
+            variables.push_back(makeVariable(computedValueNames[i]+"2", "values"+computedValues->getParameterSuffix(i, "2")));
         }
         for (int i = 0; i < force.getNumGlobalParameters(); i++)
-            variables[force.getGlobalParameterName(i)] = "globals["+intToString(i)+"]";
+            variables.push_back(makeVariable(force.getGlobalParameterName(i), "globals["+intToString(i)+"]"));
         stringstream n2EnergySource;
         bool anyExclusions = (force.getNumExclusions() > 0);
         for (int i = 0; i < force.getNumEnergyTerms(); i++) {
@@ -2182,19 +2198,22 @@ void OpenCLCalcCustomGBForceKernel::initialize(const System& system, const Custo
     {
         // Create the code to calculate chain rules terms as part of the default nonbonded kernel.
 
-        map<string, string> globalVariables;
+        vector<pair<ExpressionTreeNode, string> > globalVariables;
         for (int i = 0; i < force.getNumGlobalParameters(); i++) {
             const string& name = force.getGlobalParameterName(i);
             string value = "globals["+intToString(i)+"]";
-            globalVariables[name] = prefix+value;
+            globalVariables.push_back(makeVariable(name, prefix+value));
         }
-        map<string, string> variables = globalVariables;
+        vector<pair<ExpressionTreeNode, string> > variables = globalVariables;
         map<string, string> rename;
-        variables["r"] = "r";
+        ExpressionTreeNode rnode(new Operation::Variable("r"));
+        variables.push_back(make_pair(rnode, "r"));
+        variables.push_back(make_pair(ExpressionTreeNode(new Operation::Square(), rnode), "r2"));
+        variables.push_back(make_pair(ExpressionTreeNode(new Operation::Reciprocal(), rnode), "invR"));
         for (int i = 0; i < force.getNumPerParticleParameters(); i++) {
             const string& name = force.getPerParticleParameterName(i);
-            variables[name+"1"] = prefix+"params"+params->getParameterSuffix(i, "1");
-            variables[name+"2"] = prefix+"params"+params->getParameterSuffix(i, "2");
+            variables.push_back(makeVariable(name+"1", prefix+"params"+params->getParameterSuffix(i, "1")));
+            variables.push_back(makeVariable(name+"2", prefix+"params"+params->getParameterSuffix(i, "2")));
             rename[name+"1"] =  name+"2";
             rename[name+"2"] =  name+"1";
         }
@@ -2213,12 +2232,12 @@ void OpenCLCalcCustomGBForceKernel::initialize(const System& system, const Custo
         variables = globalVariables;
         map<string, string> rename1;
         map<string, string> rename2;
-        variables["x1"] = "posq1.x";
-        variables["y1"] = "posq1.y";
-        variables["z1"] = "posq1.z";
-        variables["x2"] = "posq2.x";
-        variables["y2"] = "posq2.y";
-        variables["z2"] = "posq2.z";
+        variables.push_back(makeVariable("x1", "posq1.x"));
+        variables.push_back(makeVariable("y1", "posq1.y"));
+        variables.push_back(makeVariable("z1", "posq1.z"));
+        variables.push_back(makeVariable("x2", "posq2.x"));
+        variables.push_back(makeVariable("y2", "posq2.y"));
+        variables.push_back(makeVariable("z2", "posq2.z"));
         rename1["x"] = "x1";
         rename1["y"] = "y1";
         rename1["z"] = "z1";
@@ -2227,15 +2246,15 @@ void OpenCLCalcCustomGBForceKernel::initialize(const System& system, const Custo
         rename2["z"] = "z2";
         for (int i = 0; i < force.getNumPerParticleParameters(); i++) {
             const string& name = force.getPerParticleParameterName(i);
-            variables[name+"1"] = prefix+"params"+params->getParameterSuffix(i, "1");
-            variables[name+"2"] = prefix+"params"+params->getParameterSuffix(i, "2");
+            variables.push_back(makeVariable(name+"1", prefix+"params"+params->getParameterSuffix(i, "1")));
+            variables.push_back(makeVariable(name+"2", prefix+"params"+params->getParameterSuffix(i, "2")));
             rename1[name] = name+"1";
             rename2[name] = name+"2";
         }
         for (int i = 0; i < force.getNumComputedValues(); i++) {
             const string& name = computedValueNames[i];
-            variables[name+"1"] = prefix+"values"+computedValues->getParameterSuffix(i, "1");
-            variables[name+"2"] = prefix+"values"+computedValues->getParameterSuffix(i, "2");
+            variables.push_back(makeVariable(name+"1", prefix+"values"+computedValues->getParameterSuffix(i, "1")));
+            variables.push_back(makeVariable(name+"2", prefix+"values"+computedValues->getParameterSuffix(i, "2")));
             rename1[name] = name+"1";
             rename2[name] = name+"2";
             if (i == 0)

platforms/opencl/src/kernels/customGBEnergyN2_cpu.cl

@@ -95,7 +95,8 @@ __kernel void computeN2Energy(__global float4* forceBuffers, __global float* ene
 #ifdef USE_CUTOFF
                     if (r2 < CUTOFF_SQUARED) {
 #endif
-                    float r = SQRT(r2);
+                    float invR = RSQRT(r2);
+                    float r = RECIP(invR);
                     unsigned int atom2 = j;
                     LOAD_ATOM2_PARAMETERS
                     atom2 = y*TILE_SIZE+j;
@@ -151,7 +152,8 @@ __kernel void computeN2Energy(__global float4* forceBuffers, __global float* ene
 #endif
                                 float r2 = dot(delta.xyz, delta.xyz);
                                 if (r2 < CUTOFF_SQUARED) {
-                                    float r = SQRT(r2);
+                                    float invR = RSQRT(r2);
+                                    float r = RECIP(invR);
                                     unsigned int atom2 = j;
                                     LOAD_ATOM2_PARAMETERS
                                     atom2 = y*TILE_SIZE+j;
@@ -204,7 +206,8 @@ __kernel void computeN2Energy(__global float4* forceBuffers, __global float* ene
 #ifdef USE_CUTOFF
                         if (r2 < CUTOFF_SQUARED) {
 #endif
-                        float r = SQRT(r2);
+                        float invR = RSQRT(r2);
+                        float r = RECIP(invR);
                         unsigned int atom2 = j;
                         LOAD_ATOM2_PARAMETERS
                         atom2 = y*TILE_SIZE+j;

platforms/opencl/src/kernels/customGBEnergyN2_default.cl

@@ -98,7 +98,8 @@ void computeN2Energy(__global float4* forceBuffers, __global float* energyBuffer
 #ifdef USE_CUTOFF
                 if (r2 < CUTOFF_SQUARED) {
 #endif
-                float r = SQRT(r2);
+                float invR = RSQRT(r2);
+                float r = RECIP(invR);
                 LOAD_ATOM2_PARAMETERS
                 atom2 = y*TILE_SIZE+baseLocalAtom+j;
                 float dEdR = 0.0f;
@@ -173,7 +174,8 @@ void computeN2Energy(__global float4* forceBuffers, __global float* energyBuffer
 #ifdef USE_CUTOFF
                 if (r2 < CUTOFF_SQUARED) {
 #endif
-                float r = SQRT(r2);
+                float invR = RSQRT(r2);
+                float r = RECIP(invR);
                 LOAD_ATOM2_PARAMETERS
                 atom2 = y*TILE_SIZE+baseLocalAtom+tj;
                 float dEdR = 0.0f;

platforms/opencl/src/kernels/customGBEnergyN2_nvidia.cl

@@ -113,7 +113,8 @@ __kernel void computeN2Energy(
 #ifdef USE_CUTOFF
                     if (r2 < CUTOFF_SQUARED) {
 #endif
-                    float r = SQRT(r2);
+                    float invR = RSQRT(r2);
+                    float r = RECIP(invR);
                     LOAD_ATOM2_PARAMETERS
                     atom2 = y*TILE_SIZE+j;
                     float dEdR = 0.0f;
@@ -175,7 +176,8 @@ __kernel void computeN2Energy(
 #ifdef USE_CUTOFF
                         if (r2 < CUTOFF_SQUARED) {
 #endif
-                        float r = SQRT(r2);
+                        float invR = RSQRT(r2);
+                        float r = RECIP(invR);
                         LOAD_ATOM2_PARAMETERS
                         atom2 = y*TILE_SIZE+tj;
                         float dEdR = 0.0f;

platforms/opencl/src/kernels/customGBValueN2_cpu.cl

@@ -92,7 +92,8 @@ __kernel void computeN2Value(__global float4* posq, __local float4* local_posq,
 #ifdef USE_CUTOFF
                     if (r2 < CUTOFF_SQUARED) {
 #endif
-                    float r = SQRT(r2);
+                    float invR = RSQRT(r2);
+                    float r = RECIP(invR);
                     unsigned int atom2 = j;
                     LOAD_ATOM2_PARAMETERS
                     atom2 = y*TILE_SIZE+j;
@@ -148,7 +149,8 @@ __kernel void computeN2Value(__global float4* posq, __local float4* local_posq,
                                 float tempValue1 = 0.0f;
                                 float tempValue2 = 0.0f;
                                 if (r2 < CUTOFF_SQUARED) {
-                                    float r = SQRT(r2);
+                                    float invR = RSQRT(r2);
+                                    float r = RECIP(invR);
                                     unsigned int atom2 = j;
                                     LOAD_ATOM2_PARAMETERS
                                     atom2 = y*TILE_SIZE+j;
@@ -194,7 +196,8 @@ __kernel void computeN2Value(__global float4* posq, __local float4* local_posq,
 #ifdef USE_CUTOFF
                         if (r2 < CUTOFF_SQUARED) {
 #endif
-                        float r = SQRT(r2);
+                        float invR = RSQRT(r2);
+                        float r = RECIP(invR);
                         unsigned int atom2 = j;
                         LOAD_ATOM2_PARAMETERS
                         atom2 = y*TILE_SIZE+j;

platforms/opencl/src/kernels/customGBValueN2_default.cl

@@ -94,7 +94,8 @@ void computeN2Value(__global float4* posq, __local float4* local_posq, __global
 #ifdef USE_CUTOFF
                 if (r2 < CUTOFF_SQUARED) {
 #endif
-                float r = SQRT(r2);
+                float invR = RSQRT(r2);
+                float r = RECIP(invR);
                 LOAD_ATOM2_PARAMETERS
                 atom2 = y*TILE_SIZE+baseLocalAtom+j;
                 float tempValue1 = 0.0f;
@@ -166,7 +167,8 @@ void computeN2Value(__global float4* posq, __local float4* local_posq, __global
 #ifdef USE_CUTOFF
                 if (r2 < CUTOFF_SQUARED) {
 #endif
-                float r = SQRT(r2);
+                float invR = RSQRT(r2);
+                float r = RECIP(invR);
                 LOAD_ATOM2_PARAMETERS
                 atom2 = y*TILE_SIZE+baseLocalAtom+tj;
                 float tempValue1 = 0.0f;

platforms/opencl/src/kernels/customGBValueN2_nvidia.cl

@@ -106,7 +106,8 @@ __kernel void computeN2Value(__global float4* posq, __local float4* local_posq,
 #ifdef USE_CUTOFF
                     if (r2 < CUTOFF_SQUARED) {
 #endif
-                    float r = SQRT(r2);
+                    float invR = RSQRT(r2);
+                    float r = RECIP(invR);
                     LOAD_ATOM2_PARAMETERS
                     atom2 = y*TILE_SIZE+j;
                     float tempValue1 = 0.0f;
@@ -160,7 +161,8 @@ __kernel void computeN2Value(__global float4* posq, __local float4* local_posq,
                                 float tempValue1 = 0.0f;
                                 float tempValue2 = 0.0f;
                                 if (r2 < CUTOFF_SQUARED) {
-                                    float r = SQRT(r2);
+                                    float invR = RSQRT(r2);
+                                    float r = RECIP(invR);
                                     LOAD_ATOM2_PARAMETERS
                                     atom2 = y*TILE_SIZE+j;
                                     if (atom1 < NUM_ATOMS && atom2 < NUM_ATOMS) {
@@ -206,7 +208,8 @@ __kernel void computeN2Value(__global float4* posq, __local float4* local_posq,
 #ifdef USE_CUTOFF
                         if (r2 < CUTOFF_SQUARED) {
 #endif
-                        float r = SQRT(r2);
+                        float invR = RSQRT(r2);
+                        float r = RECIP(invR);
                         LOAD_ATOM2_PARAMETERS
                         atom2 = y*TILE_SIZE+tj;
                         float tempValue1 = 0.0f;

tests/TestParser.cpp

@@ -239,6 +239,9 @@ int main() {
         verifyDerivative("abs(3*x)", "step(3*x)*3+(1-step(3*x))*-3");
         testCustomFunction("custom(x, y)/2", "x*y");
         testCustomFunction("custom(x^2, 1)+custom(2, y-1)", "2*x^2+4*(y-1)");
+        cout << Parser::parse("x*x").optimize() << endl;
+        cout << Parser::parse("x*(x*x)").optimize() << endl;
+        cout << Parser::parse("(x*x)*x").optimize() << endl;
         cout << Parser::parse("2*3*x").optimize() << endl;
         cout << Parser::parse("1/(1+x)").optimize() << endl;
         cout << Parser::parse("x^(1/2)").optimize() << endl;