+ * CustomCentroidBondForce* force = new CustomCentroidBondForce(2, "0.5*k*distance(g1,g2)^2");
+ * force->addPerBondParameter("k");
+ * force->addGroup(particles1);
+ * force->addGroup(particles2);
+ * vector bondGroups;
+ * bondGroups.push_back(0);
+ * bondGroups.push_back(1);
+ * vector bondParameters;
+ * bondParameters.push_back(k);
+ * force->addBond(bondGroups, bondParameters);
+ *
+ *
+ * Expressions may involve the operators + (add), - (subtract), * (multiply), / (divide), and ^ (power), and the following
+ * functions: sqrt, exp, log, sin, cos, sec, csc, tan, cot, asin, acos, atan, sinh, cosh, tanh, erf, erfc, min, max, abs, floor, ceil, step, delta, select. All trigonometric functions
+ * are defined in radians, and log is the natural logarithm. step(x) = 0 if x is less than 0, 1 otherwise. delta(x) = 1 if x is 0, 0 otherwise.
+ * select(x,y,z) = z if x = 0, y otherwise.
+ *
+ * In addition, you can call addTabulatedFunction() to define a new function based on tabulated values. You specify the function by
+ * creating a TabulatedFunction object. That function can then appear in the expression.
+ */
+
+class OPENMM_EXPORT CustomCentroidBondForce : public Force {
+public:
+ /**
+ * Create a CustomCentroidBondForce.
+ *
+ * @param numGroups the number of groups used to define each bond
+ * @param energy an algebraic expression giving the interaction energy of each bond as a function
+ * of particle positions, inter-particle distances, angles, and dihedrals, and any global
+ * and per-bond parameters
+ */
+ explicit CustomCentroidBondForce(int numGroups, const std::string& energy);
+ ~CustomCentroidBondForce();
+ /**
+ * Get the number of groups used to define each bond.
+ */
+ int getNumGroupsPerBond() const {
+ return groupsPerBond;
+ }
+ /**
+ * Get the number of particle groups that have been defined.
+ */
+ int getNumGroups() const {
+ return groups.size();
+ }
+ /**
+ * Get the number of bonds for which force field parameters have been defined.
+ */
+ int getNumBonds() const {
+ return bonds.size();
+ }
+ /**
+ * Get the number of per-bond parameters that the interaction depends on.
+ */
+ int getNumPerBondParameters() const {
+ return bondParameters.size();
+ }
+ /**
+ * Get the number of global parameters that the interaction depends on.
+ */
+ int getNumGlobalParameters() const {
+ return globalParameters.size();
+ }
+ /**
+ * Get the number of tabulated functions that have been defined.
+ */
+ int getNumTabulatedFunctions() const {
+ return functions.size();
+ }
+ /**
+ * Get the number of tabulated functions that have been defined.
+ *
+ * @deprecated This method exists only for backward compatibility. Use getNumTabulatedFunctions() instead.
+ */
+ int getNumFunctions() const {
+ return functions.size();
+ }
+ /**
+ * Get the algebraic expression that gives the interaction energy of each bond
+ */
+ const std::string& getEnergyFunction() const;
+ /**
+ * Set the algebraic expression that gives the interaction energy of each bond
+ */
+ void setEnergyFunction(const std::string& energy);
+ /**
+ * Add a new per-bond parameter that the interaction may depend on.
+ *
+ * @param name the name of the parameter
+ * @return the index of the parameter that was added
+ */
+ int addPerBondParameter(const std::string& name);
+ /**
+ * Get the name of a per-bond parameter.
+ *
+ * @param index the index of the parameter for which to get the name
+ * @return the parameter name
+ */
+ const std::string& getPerBondParameterName(int index) const;
+ /**
+ * Set the name of a per-bond parameter.
+ *
+ * @param index the index of the parameter for which to set the name
+ * @param name the name of the parameter
+ */
+ void setPerBondParameterName(int index, const std::string& name);
+ /**
+ * Add a new global parameter that the interaction may depend on.
+ *
+ * @param name the name of the parameter
+ * @param defaultValue the default value of the parameter
+ * @return the index of the parameter that was added
+ */
+ int addGlobalParameter(const std::string& name, double defaultValue);
+ /**
+ * Get the name of a global parameter.
+ *
+ * @param index the index of the parameter for which to get the name
+ * @return the parameter name
+ */
+ const std::string& getGlobalParameterName(int index) const;
+ /**
+ * Set the name of a global parameter.
+ *
+ * @param index the index of the parameter for which to set the name
+ * @param name the name of the parameter
+ */
+ void setGlobalParameterName(int index, const std::string& name);
+ /**
+ * Get the default value of a global parameter.
+ *
+ * @param index the index of the parameter for which to get the default value
+ * @return the parameter default value
+ */
+ double getGlobalParameterDefaultValue(int index) const;
+ /**
+ * Set the default value of a global parameter.
+ *
+ * @param index the index of the parameter for which to set the default value
+ * @param name the default value of the parameter
+ */
+ void setGlobalParameterDefaultValue(int index, double defaultValue);
+ /**
+ * Add a particle group.
+ *
+ * @param particles the indices of the particles to include in the group
+ * @param weights the weight to use for each particle when computing the center position.
+ * If this is omitted, then particle masses will be used as weights.
+ * @return the index of the group that was added
+ */
+ int addGroup(const std::vector