Replace C++ code examples in Python API docs (#5043)

* Replace C++ code examples in Python API docs

* Use newer Python for building docs

* Remove sphinx version pin

* Add Python versions of code examples

* Minor edit

* Add Python version of one more code example

.github/workflows/CI.yml

@@ -538,7 +538,7 @@ jobs:
       - uses: conda-incubator/setup-miniconda@v2
         name: "Prepare base dependencies"
         with:
-          python-version: "3.7"
+          python-version: "3.13"
           activate-environment: build
           environment-file: devtools/ci/gh-actions/conda-envs/docs.yml
           auto-activate-base: false

devtools/ci/gh-actions/conda-envs/docs.yml

@@ -13,6 +13,6 @@ dependencies:
 - swig
 - doxygen
 - setuptools
-- sphinx==4.0.2
+- sphinx
 - sphinxcontrib-bibtex
 - breathe>=4.30,<5.0

docs-source/api-c++/CMakeLists.txt

@@ -19,7 +19,9 @@ configure_file(
 add_custom_command(
     OUTPUT "${WRAPPER_DOXYGEN_DIR}/xml/index.xml"
     COMMAND "${DOXYGEN_EXECUTABLE}"
+    COMMAND "${PYTHON_EXECUTABLE}" "${CMAKE_CURRENT_BINARY_DIR}/process-docstring.py" "${WRAPPER_DOXYGEN_DIR}/xml"
     DEPENDS "${WRAPPER_DOXYGEN_DIR}/Doxyfile"
+            "${CMAKE_CURRENT_BINARY_DIR}/process-docstring.py"
     WORKING_DIRECTORY "${WRAPPER_DOXYGEN_DIR}"
     COMMENT "Parsing OpenMM header files with Doxygen..."
 )

docs-source/api-c++/forces.rst

@@ -76,8 +76,8 @@ be combined in arbitrary ways.
 
 .. _custom-forces:
 
-Customizing ``Force``
-=====================
+Custom Forces
+=============
 
 OpenMM provides a number of classes that make it easier to implement custom
 forces for common scenarios. These classes implement constructors that take an

docs-source/api-c++/process-docstring.py

+import os
+import re
+import sys
+
+def replace_code(match):
+    lines = match.group(1).split('\n')
+    for i, line in enumerate(list(lines)):
+        if line.startswith('*'):
+            lines[i] = line[1:]
+    return f"\n* .. code-block:: c++\n*    {'\n*    '.join(lines)}"
+
+
+def process_file(file):
+    """Edit docstrings in an XML file to remove the Python versions of code examples,
+    and put the correct Sphinx directives around the C++ versions."""
+    with open(file) as input:
+        content = input.read()
+    processed = re.sub(r'<c\+\+>((.|\n)*?)</c\+\+>', replace_code, content)
+    processed = re.sub(r'<python>((.|\n)*?)</python>', '', processed)
+    print(file, processed != content)
+    if processed != content:
+        with open(file, 'w') as output:
+            output.write(processed)
+
+
+dir = sys.argv[1]
+for file in os.listdir(dir):
+    process_file(os.path.join(dir, file))

docs-source/api-python/process-docstring.py

@@ -51,6 +51,12 @@ def process_docstring(app, what, name, obj, options, lines):
                 s = linesep + s
             newline = '.. verbatim::' + linesep
             return newline + '    ' + s.replace(linesep, linesep + '    ')
+    def replace_code(m):
+        s = m.group(1)
+        if not s.startswith(linesep):
+            s = linesep + s
+        newline = linesep + '.. code-block:: python' + linesep
+        return newline + '    ' + s.replace(linesep, linesep + '    ')
     def replace_subscript(m):
         """ Replace subscript tags. """
         return r'\ :sub:`{}`\ '.format(m.group(1))
@@ -68,6 +74,8 @@ def process_docstring(app, what, name, obj, options, lines):
     joined = re.sub(r'<verbatim>(.*?)</verbatim>', repl4, joined)
     joined = re.sub(r'<sub>(.*?)</sub>', replace_subscript, joined)
     joined = re.sub(r'<sup>(.*?)</sup>', replace_superscript, joined)
+    joined = re.sub(r'<c\+\+>(.*?)</c\+\+>', '', joined)
+    joined = re.sub(r'<python>(.*?)</python>', replace_code, joined)
 
     lines[:] = [(l if not l.isspace() else '') for l in joined.split(linesep)]
 
@@ -90,10 +98,10 @@ def convert_type(type):
     if type in substitutions:
         type = substitutions[type]
     if '<' in type:
-        match = re.match('vector<(.*?),.*>', type)
+        match = re.match(r'vector<(.*?),.*>', type)
         if match is not None:
             type = f'tuple[{convert_type(match[1])}]'
-        match = re.match('map<(.*?),(.*?),.*>', type)
+        match = re.match(r'map<(.*?),(.*?),.*>', type)
         if match is not None:
             type = f'Mapping[{convert_type(match[1])}, {convert_type(match[2])}]'
     return type

olla/include/openmm/Platform.h

@@ -53,10 +53,12 @@ class KernelFactory;
  * To get a Platform object, call
  *
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    Platform& platform = Platform::findPlatform(kernelNames);
- *
+ * <c++>
+ * Platform& platform = Platform::findPlatform(kernelNames);
+ * </c++>
+ * <python>
+ * platform = Platform.findPlatform(kernelNames)
+ * </python>
  * \endverbatim
  *
  * passing in the names of all kernels that will be required for the calculation you plan to perform.  It

openmmapi/include/openmm/ATMForce.h

@@ -80,14 +80,22 @@ namespace OpenMM {
  * The energy change is dialed using an alchemical parameter Lambda, which in this case is set to 1/2:
  *
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    ATMForce *atmforce = new ATMForce("u0 + Lambda*(u1 - u0)");
- *    atm->addGlobalParameter("Lambda", 0.5);
- *    atm->addParticle();
- *    atm->addParticle(new ATMForce::FixedDisplacement(Vec3(1, 0, 0)));
- *    CustomBondForce* force = new CustomBondForce("0.5*r^2");
- *    atm->addForce(force);
+ * <c++>
+ * ATMForce *atmforce = new ATMForce("u0 + Lambda*(u1 - u0)");
+ * atm->addGlobalParameter("Lambda", 0.5);
+ * atm->addParticle();
+ * atm->addParticle(new ATMForce::FixedDisplacement(Vec3(1, 0, 0)));
+ * CustomBondForce* force = new CustomBondForce("0.5*r^2");
+ * atm->addForce(force);
+ * </c++>
+ * <python>
+ * atmforce = ATMForce("u0 + Lambda*(u1 - u0)")
+ * atm.addGlobalParameter("Lambda", 0.5)
+ * atm.addParticle()
+ * atm.addParticle(FixedDisplacement(Vec3(1, 0, 0)))
+ * force = CustomBondForce("0.5*r^2")
+ * atm.addForce(force)
+ * </python>
  * \endverbatim
  *
  * Note that calling addParticle() without arguments is equivalent to a zero fixed displacement.
@@ -98,15 +106,24 @@ namespace OpenMM {
  * pos[2]-pos[1] going from the second particle to the third particle, 
  *
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    ATMForce *atmforce = new ATMForce("u0 + Lambda*(u1 - u0)");
- *    atm->addGlobalParameter("Lambda", 0.5);
- *    atm->addParticle(new ATMForce::ParticleOffsetDisplacement(2, 1));
- *    atm->addParticle();
- *    atm->addParticle();
- *    CustomBondForce* force = new CustomBondForce("0.5*r^2");
- *    atm->addForce(force);
+ * <c++>
+ * ATMForce *atmforce = new ATMForce("u0 + Lambda*(u1 - u0)");
+ * atm->addGlobalParameter("Lambda", 0.5);
+ * atm->addParticle(new ATMForce::ParticleOffsetDisplacement(2, 1));
+ * atm->addParticle();
+ * atm->addParticle();
+ * CustomBondForce* force = new CustomBondForce("0.5*r^2");
+ * atm->addForce(force);
+ * </c++>
+ * <python>
+ * atmforce = ATMForce("u0 + Lambda*(u1 - u0)")
+ * atm.addGlobalParameter("Lambda", 0.5)
+ * atm.addParticle(ParticleOffsetDisplacement(2, 1))
+ * atm.addParticle()
+ * atm.addParticle()
+ * force = CustomBondForce("0.5*r^2")
+ * atm.addForce(force)
+ * </python>
  * \endverbatim
  *
  * where ParticleOffsetDisplacement is the class that describes this particular type of coordinate transformation.
@@ -136,9 +153,12 @@ namespace OpenMM {
  * The ATMForce is then added to the System as any other Force
  *
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *  system.addForce(atmforce);
+ * <c++>
+ * system.addForce(atmforce);
+ * </c++>
+ * <python>
+ * system.addForce(atmforce)
+ * </python>
  * \endverbatim
  *
  * after which it will be used for energy/force evaluations for molecular dynamics and energy optimization.
@@ -151,18 +171,24 @@ namespace OpenMM {
  * vectors to the fixed displacement transformation given to addParticle(). For example, with:
  *
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    atm->addParticle(new ATMForce::FixedDisplacement(Vec3(1, 0, 0), Vec3(-1, 0, 0)));
+ * <c++>
+ * atm->addParticle(new ATMForce::FixedDisplacement(Vec3(1, 0, 0), Vec3(-1, 0, 0)));
+ * </c++>
+ * <python>
+ * atm.addParticle(FixedDisplacement(Vec3(1, 0, 0), Vec3(-1, 0, 0)))
+ * </python>
  * \endverbatim
  * 
  * the energy u1 will be computed after displacing the particle in the positive x direction, and
  * u0 will be computed after displacing it in the negative x direction. Similarly,
  *
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    atm->addParticle(new ATMForce::ParticleOffsetDisplacement(4, 3, 2, 1));
+ * <c++>
+ * atm->addParticle(new ATMForce::ParticleOffsetDisplacement(4, 3, 2, 1));
+ * </c++>
+ * <python>
+ * atm.addParticle(ParticleOffsetDisplacement(4, 3, 2, 1))
+ * </python>
  * \endverbatim
  *
  * adds a particle whose position is displaced by pos[4]-pos[3] before calculating u1 and by pos[2]-pos[1] before calculating u0.

openmmapi/include/openmm/CompoundIntegrator.h

@@ -44,35 +44,46 @@ namespace OpenMM {
  * you need, then add all of them to a CustomIntegrator:
  * 
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    CompoundIntegrator compoundIntegrator;
- *    compoundIntegrator.addIntegrator(new VerletIntegrator(0.001));
- *    compoundIntegrator.addIntegrator(new LangevinIntegrator(300.0, 1.0, 0.001));
- *
+ * <c++>
+ * CompoundIntegrator compoundIntegrator;
+ * compoundIntegrator.addIntegrator(new VerletIntegrator(0.001));
+ * compoundIntegrator.addIntegrator(new LangevinIntegrator(300.0, 1.0, 0.001));
+ * </c++>
+ * <python>
+ * compoundIntegrator = CompoundIntegrator()
+ * compoundIntegrator.addIntegrator(VerletIntegrator(0.001))
+ * compoundIntegrator.addIntegrator(LangevinIntegrator(300.0, 1.0, 0.001))
+ * </python>
  * \endverbatim
  * 
  * Next create a Context, specifying the CompoundIntegrator as the Integrator to use for
  * the Context:
  * 
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *     Context context(system, compoundIntegrator);
- *
+ * <c++>
+ * Context context(system, compoundIntegrator);
+ * </c++>
+ * <python>
+ * context = Context(system, compoundIntegrator)
+ * </python>
  * \endverbatim
  * 
  * Finally, call setCurrentIntegrator() to set which Integrator is active.  That one will
  * be used for all calls to step() until the next time you change it.
  * 
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *     compoundIntegrator.setCurrentIntegrator(0);
- *     compoundIntegrator.step(1000); // Take 1000 steps of Verlet dynamics
- *     compoundIntegrator.setCurrentIntegrator(1);
- *     compoundIntegrator.step(1000); // Take 1000 steps of Langevin dynamics
- *
+ * <c++>
+ * compoundIntegrator.setCurrentIntegrator(0);
+ * compoundIntegrator.step(1000); // Take 1000 steps of Verlet dynamics
+ * compoundIntegrator.setCurrentIntegrator(1);
+ * compoundIntegrator.step(1000); // Take 1000 steps of Langevin dynamics
+ * </c++>
+ * <python>
+ * compoundIntegrator.setCurrentIntegrator(0)
+ * compoundIntegrator.step(1000) # Take 1000 steps of Verlet dynamics
+ * compoundIntegrator.setCurrentIntegrator(1)
+ * compoundIntegrator.step(1000) # Take 1000 steps of Langevin dynamics
+ * </python>
  * \endverbatim
  * 
  * When switching between integrators, it is important to make sure they are compatible with

openmmapi/include/openmm/CustomAngleForce.h

@@ -56,20 +56,25 @@ namespace OpenMM {
  * As an example, the following code creates a CustomAngleForce that implements a harmonic potential:
  *
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    CustomAngleForce* force = new CustomAngleForce("0.5*k*(theta-theta0)^2");
- *
+ * <c++>
+ * CustomAngleForce* force = new CustomAngleForce("0.5*k*(theta-theta0)^2");
+ * </c++>
+ * <python>
+ * force = CustomAngleForce("0.5*k*(theta-theta0)^2")
+ * </python>
  * \endverbatim
  *
  * This force depends on two parameters: the spring constant k and equilibrium angle theta0.  The following code defines these parameters:
  *
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    force->addPerAngleParameter("k");
- *    force->addPerAngleParameter("theta0");
- *
+ * <c++>
+ * force->addPerAngleParameter("k");
+ * force->addPerAngleParameter("theta0");
+ * </c++>
+ * <python>
+ * force.addPerAngleParameter("k")
+ * force.addPerAngleParameter("theta0")
+ * </python>
  * \endverbatim
  * 
  * This class also has the ability to compute derivatives of the potential energy with respect to global parameters.

openmmapi/include/openmm/CustomBondForce.h

@@ -56,20 +56,25 @@ namespace OpenMM {
  * As an example, the following code creates a CustomBondForce that implements a harmonic potential:
  *
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    CustomBondForce* force = new CustomBondForce("0.5*k*(r-r0)^2");
- *
+ * <c++>
+ * CustomBondForce* force = new CustomBondForce("0.5*k*(r-r0)^2");
+ * </c++>
+ * <python>
+ * force = CustomBondForce("0.5*k*(r-r0)^2")
+ * </python>
  * \endverbatim
  *
  * This force depends on two parameters: the spring constant k and equilibrium distance r0.  The following code defines these parameters:
  *
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    force->addPerBondParameter("k");
- *    force->addPerBondParameter("r0");
- *
+ * <c++>
+ * force->addPerBondParameter("k");
+ * force->addPerBondParameter("r0");
+ * </c++>
+ * <python>
+ * force.addPerBondParameter("k")
+ * force.addPerBondParameter("r0")
+ * </python>
  * \endverbatim
  * 
  * This class also has the ability to compute derivatives of the potential energy with respect to global parameters.

openmmapi/include/openmm/CustomCentroidBondForce.h

@@ -87,19 +87,27 @@ namespace OpenMM {
  * centers of mass of two groups of particles.
  *
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    CustomCentroidBondForce* force = new CustomCentroidBondForce(2, "0.5*k*distance(g1,g2)^2");
- *    force->addPerBondParameter("k");
- *    force->addGroup(particles1);
- *    force->addGroup(particles2);
- *    vector<int> bondGroups;
- *    bondGroups.push_back(0);
- *    bondGroups.push_back(1);
- *    vector<double> bondParameters;
- *    bondParameters.push_back(k);
- *    force->addBond(bondGroups, bondParameters);
- *
+ * <c++>
+ * CustomCentroidBondForce* force = new CustomCentroidBondForce(2, "0.5*k*distance(g1,g2)^2");
+ * force->addPerBondParameter("k");
+ * force->addGroup(particles1);
+ * force->addGroup(particles2);
+ * vector<int> bondGroups;
+ * bondGroups.push_back(0);
+ * bondGroups.push_back(1);
+ * vector<double> bondParameters;
+ * bondParameters.push_back(k);
+ * force->addBond(bondGroups, bondParameters);
+ * </c++>
+ * <python>
+ * force = CustomCentroidBondForce(2, "0.5*k*distance(g1,g2)^2")
+ * force.addPerBondParameter("k")
+ * force.addGroup(particles1)
+ * force.addGroup(particles2)
+ * bondGroups = [0, 1]
+ * bondParameters = [k]
+ * force.addBond(bondGroups, bondParameters)
+ * </python>
  * \endverbatim
  * 
  * This class also has the ability to compute derivatives of the potential energy with respect to global parameters.
@@ -119,10 +127,12 @@ namespace OpenMM {
  * from group g1 to the midpoint between groups g2 and g3.
  * 
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    CustomCentroidBondForce* force = new CustomCentroidBondForce(3, "pointdistance(x1, y1, z1, (x2+x3)/2, (y2+y3)/2, (z2+z3)/2)");
- *
+ * <c++>
+ * CustomCentroidBondForce* force = new CustomCentroidBondForce(3, "pointdistance(x1, y1, z1, (x2+x3)/2, (y2+y3)/2, (z2+z3)/2)");
+ * </c++>
+ * <python>
+ * force = CustomCentroidBondForce(3, "pointdistance(x1, y1, z1, (x2+x3)/2, (y2+y3)/2, (z2+z3)/2)")
+ * </python>
  * \endverbatim
  *
  * In addition, you can call addTabulatedFunction() to define a new function based on tabulated values.  You specify the function by

openmmapi/include/openmm/CustomCompoundBondForce.h

@@ -78,22 +78,29 @@ namespace OpenMM {
  * p1 and p3.
  *
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    CustomCompoundBondForce* force = new CustomCompoundBondForce(3, "0.5*(kangle*(angle(p1,p2,p3)-theta0)^2+kbond*(distance(p1,p3)-r0)^2)");
- *
+ * <c++>
+ * CustomCompoundBondForce* force = new CustomCompoundBondForce(3, "0.5*(kangle*(angle(p1,p2,p3)-theta0)^2+kbond*(distance(p1,p3)-r0)^2)");
+ * </c++>
+ * <python>
+ * force = CustomCompoundBondForce(3, "0.5*(kangle*(angle(p1,p2,p3)-theta0)^2+kbond*(distance(p1,p3)-r0)^2)")
+ * </python>
  * \endverbatim
  *
  * This force depends on four parameters: kangle, kbond, theta0, and r0.  The following code defines these as per-bond parameters:
  *
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    force->addPerBondParameter("kangle");
- *    force->addPerBondParameter("kbond");
- *    force->addPerBondParameter("theta0");
- *    force->addPerBondParameter("r0");
- *
+ * <c++>
+ * force->addPerBondParameter("kangle");
+ * force->addPerBondParameter("kbond");
+ * force->addPerBondParameter("theta0");
+ * force->addPerBondParameter("r0");
+ * </c++>
+ * <python>
+ * force.addPerBondParameter("kangle")
+ * force.addPerBondParameter("kbond")
+ * force.addPerBondParameter("theta0")
+ * force.addPerBondParameter("r0")
+ * </python>
  * \endverbatim
  *
  * This class also has the ability to compute derivatives of the potential energy with respect to global parameters.
@@ -113,10 +120,12 @@ namespace OpenMM {
  * from particle p1 to the midpoint between particles p2 and p3.
  * 
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    CustomCompoundBondForce* force = new CustomCompoundBondForce(3, "pointdistance(x1, y1, z1, (x2+x3)/2, (y2+y3)/2, (z2+z3)/2)");
- *
+ * <c++>
+ * CustomCompoundBondForce* force = new CustomCompoundBondForce(3, "pointdistance(x1, y1, z1, (x2+x3)/2, (y2+y3)/2, (z2+z3)/2)");
+ * </c++>
+ * <python>
+ * force = CustomCompoundBondForce(3, "pointdistance(x1, y1, z1, (x2+x3)/2, (y2+y3)/2, (z2+z3)/2)")
+ * </python>
  * \endverbatim
  *
  * In addition, you can call addTabulatedFunction() to define a new function based on tabulated values.  You specify the function by

openmmapi/include/openmm/CustomExternalForce.h

@@ -57,22 +57,29 @@ namespace OpenMM {
  * via a harmonic potential:
  *
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    CustomExternalForce* force = new CustomExternalForce("k*((x-x0)^2+(y-y0)^2+(z-z0)^2)");
- *
+ * <c++>
+ * CustomExternalForce* force = new CustomExternalForce("k*((x-x0)^2+(y-y0)^2+(z-z0)^2)");
+ * </c++>
+ * <python>
+ * force = CustomExternalForce("k*((x-x0)^2+(y-y0)^2+(z-z0)^2)")
+ * </python>
  * \endverbatim
  *
  * This force depends on four parameters: the spring constant k and equilibrium coordinates x0, y0, and z0.  The following code defines these parameters:
  *
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    force->addGlobalParameter("k", 100.0);
- *    force->addPerParticleParameter("x0");
- *    force->addPerParticleParameter("y0");
- *    force->addPerParticleParameter("z0");
- *
+ * <c++>
+ * force->addGlobalParameter("k", 100.0);
+ * force->addPerParticleParameter("x0");
+ * force->addPerParticleParameter("y0");
+ * force->addPerParticleParameter("z0");
+ * </c++>
+ * <python>
+ * force.addGlobalParameter("k", 100.0)
+ * force.addPerParticleParameter("x0")
+ * force.addPerParticleParameter("y0")
+ * force.addPerParticleParameter("z0")
+ * </python>
  * \endverbatim
  *
  * Special care is needed in systems that use periodic boundary conditions.  In that case, each particle really represents
@@ -82,10 +89,12 @@ namespace OpenMM {
  * periodic copies of the points (x1, y1, z1) and (x2, y2, z2).  For example, the force given above would be rewritten as
  *
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    CustomExternalForce* force = new CustomExternalForce("k*periodicdistance(x, y, z, x0, y0, z0)^2");
- *
+ * <c++>
+ * CustomExternalForce* force = new CustomExternalForce("k*periodicdistance(x, y, z, x0, y0, z0)^2");
+ * </c++>
+ * <python>
+ * force = CustomExternalForce("k*periodicdistance(x, y, z, x0, y0, z0)^2")
+ * </python>
  * \endverbatim
  *
  * Expressions may involve the operators + (add), - (subtract), * (multiply), / (divide), and ^ (power), and the following

openmmapi/include/openmm/CustomGBForce.h

@@ -83,28 +83,48 @@ namespace OpenMM {
  * of the GB/SA solvation model, using the ACE approximation to estimate surface area:
  *
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    CustomGBForce* custom = new CustomGBForce();
- *    custom->addPerParticleParameter("q");
- *    custom->addPerParticleParameter("radius");
- *    custom->addPerParticleParameter("scale");
- *    custom->addGlobalParameter("solventDielectric", obc->getSolventDielectric());
- *    custom->addGlobalParameter("soluteDielectric", obc->getSoluteDielectric());
- *    custom->addComputedValue("I", "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(1/U^2-1/L^2)*(r-sr2*sr2/r)+0.5*log(L/U)/r+C);"
- *                                  "U=r+sr2;"
- *                                  "C=2*(1/or1-1/L)*step(sr2-r-or1);"
- *                                  "L=max(or1, D);"
- *                                  "D=abs(r-sr2);"
- *                                  "sr2 = scale2*or2;"
- *                                  "or1 = radius1-0.009; or2 = radius2-0.009", CustomGBForce::ParticlePairNoExclusions);
- *    custom->addComputedValue("B", "1/(1/or-tanh(1*psi-0.8*psi^2+4.85*psi^3)/radius);"
- *                                  "psi=I*or; or=radius-0.009", CustomGBForce::SingleParticle);
- *    custom->addEnergyTerm("28.3919551*(radius+0.14)^2*(radius/B)^6-0.5*138.935456*(1/soluteDielectric-1/solventDielectric)*q^2/B",
- *                          CustomGBForce::SingleParticle);
- *    custom->addEnergyTerm("-138.935456*(1/soluteDielectric-1/solventDielectric)*q1*q2/f;"
- *                          "f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))", CustomGBForce::ParticlePair);
- *
+ * <c++>
+ * CustomGBForce* custom = new CustomGBForce();
+ * custom->addPerParticleParameter("q");
+ * custom->addPerParticleParameter("radius");
+ * custom->addPerParticleParameter("scale");
+ * custom->addGlobalParameter("solventDielectric", obc->getSolventDielectric());
+ * custom->addGlobalParameter("soluteDielectric", obc->getSoluteDielectric());
+ * custom->addComputedValue("I", "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(1/U^2-1/L^2)*(r-sr2*sr2/r)+0.5*log(L/U)/r+C);"
+ *                               "U=r+sr2;"
+ *                               "C=2*(1/or1-1/L)*step(sr2-r-or1);"
+ *                               "L=max(or1, D);"
+ *                               "D=abs(r-sr2);"
+ *                               "sr2 = scale2*or2;"
+ *                               "or1 = radius1-0.009; or2 = radius2-0.009", CustomGBForce::ParticlePairNoExclusions);
+ * custom->addComputedValue("B", "1/(1/or-tanh(1*psi-0.8*psi^2+4.85*psi^3)/radius);"
+ *                               "psi=I*or; or=radius-0.009", CustomGBForce::SingleParticle);
+ * custom->addEnergyTerm("28.3919551*(radius+0.14)^2*(radius/B)^6-0.5*138.935456*(1/soluteDielectric-1/solventDielectric)*q^2/B",
+ *                       CustomGBForce::SingleParticle);
+ * custom->addEnergyTerm("-138.935456*(1/soluteDielectric-1/solventDielectric)*q1*q2/f;"
+ *                       "f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))", CustomGBForce::ParticlePair);
+ * </c++>
+ * <python>
+ * custom = CustomGBForce()
+ * custom.addPerParticleParameter("q")
+ * custom.addPerParticleParameter("radius")
+ * custom.addPerParticleParameter("scale")
+ * custom.addGlobalParameter("solventDielectric", obc.getSolventDielectric())
+ * custom.addGlobalParameter("soluteDielectric", obc.getSoluteDielectric())
+ * custom.addComputedValue("I", "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(1/U^2-1/L^2)*(r-sr2*sr2/r)+0.5*log(L/U)/r+C);"
+ *                              "U=r+sr2;"
+ *                              "C=2*(1/or1-1/L)*step(sr2-r-or1);"
+ *                              "L=max(or1, D);"
+ *                              "D=abs(r-sr2);"
+ *                              "sr2 = scale2*or2;"
+ *                              "or1 = radius1-0.009; or2 = radius2-0.009", CustomGBForce.ParticlePairNoExclusions)
+ * custom.addComputedValue("B", "1/(1/or-tanh(1*psi-0.8*psi^2+4.85*psi^3)/radius);"
+ *                              "psi=I*or; or=radius-0.009", CustomGBForce.SingleParticle)
+ * custom.addEnergyTerm("28.3919551*(radius+0.14)^2*(radius/B)^6-0.5*138.935456*(1/soluteDielectric-1/solventDielectric)*q^2/B",
+ *                      CustomGBForce.SingleParticle)
+ * custom.addEnergyTerm("-138.935456*(1/soluteDielectric-1/solventDielectric)*q1*q2/f;"
+ *                      "f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))", CustomGBForce.ParticlePair)
+ * </python>
  * \endverbatim
  *
  * It begins by defining three per-particle parameters (charge, atomic radius, and scale factor) and two global parameters

openmmapi/include/openmm/CustomHbondForce.h

@@ -79,21 +79,27 @@ namespace OpenMM {
  * to keep the distance between a1 and d1, and the angle formed by a1-d1-d2, near ideal values:
  *
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    CustomHbondForce* force = new CustomHbondForce("k*(distance(a1,d1)-r0)^2*(angle(a1,d1,d2)-theta0)^2");
- *
+ * <c++>
+ * CustomHbondForce* force = new CustomHbondForce("k*(distance(a1,d1)-r0)^2*(angle(a1,d1,d2)-theta0)^2");
+ * </c++>
+ * <python>
+ * force = CustomHbondForce("k*(distance(a1,d1)-r0)^2*(angle(a1,d1,d2)-theta0)^2")
+ * </python>
  * \endverbatim
  *
  * This force depends on three parameters: k, r0, and theta0.  The following code defines these as per-donor parameters:
  *
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    force->addPerDonorParameter("k");
- *    force->addPerDonorParameter("r0");
- *    force->addPerDonorParameter("theta0");
- *
+ * <c++>
+ * force->addPerDonorParameter("k");
+ * force->addPerDonorParameter("r0");
+ * force->addPerDonorParameter("theta0");
+ * </c++>
+ * <python>
+ * force.addPerDonorParameter("k")
+ * force.addPerDonorParameter("r0")
+ * force.addPerDonorParameter("theta0")
+ * </python>
  * \endverbatim
  *
  * Expressions may involve the operators + (add), - (subtract), * (multiply), / (divide), and ^ (power), and the following

openmmapi/include/openmm/CustomIntegrator.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) 2011-2020 Stanford University and the Authors.      *
+ * Portions copyright (c) 2011-2025 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -129,13 +129,18 @@ namespace OpenMM {
  * integrator:
  *
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    CustomIntegrator integrator(0.001);
- *    integrator.addComputePerDof("v", "v+0.5*dt*f/m");
- *    integrator.addComputePerDof("x", "x+dt*v");
- *    integrator.addComputePerDof("v", "v+0.5*dt*f/m");
- *
+ * <c++>
+ * CustomIntegrator integrator(0.001);
+ * integrator.addComputePerDof("v", "v+0.5*dt*f/m");
+ * integrator.addComputePerDof("x", "x+dt*v");
+ * integrator.addComputePerDof("v", "v+0.5*dt*f/m");
+ * </c++>
+ * <python>
+ * integrator = CustomIntegrator(0.001)
+ * integrator.addComputePerDof("v", "v+0.5*dt*f/m")
+ * integrator.addComputePerDof("x", "x+dt*v")
+ * integrator.addComputePerDof("v", "v+0.5*dt*f/m")
+ * </python>
  * \endverbatim
  *
  * The first step updates the velocities based on the current forces.
@@ -155,18 +160,28 @@ namespace OpenMM {
  * both these problems, using the RATTLE algorithm to apply constraints:
  *
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    CustomIntegrator integrator(0.001);
- *    integrator.addPerDofVariable("x1", 0);
- *    integrator.addUpdateContextState();
- *    integrator.addComputePerDof("v", "v+0.5*dt*f/m");
- *    integrator.addComputePerDof("x", "x+dt*v");
- *    integrator.addComputePerDof("x1", "x");
- *    integrator.addConstrainPositions();
- *    integrator.addComputePerDof("v", "v+0.5*dt*f/m+(x-x1)/dt");
- *    integrator.addConstrainVelocities();
- *
+ * <c++>
+ * CustomIntegrator integrator(0.001);
+ * integrator.addPerDofVariable("x1", 0);
+ * integrator.addUpdateContextState();
+ * integrator.addComputePerDof("v", "v+0.5*dt*f/m");
+ * integrator.addComputePerDof("x", "x+dt*v");
+ * integrator.addComputePerDof("x1", "x");
+ * integrator.addConstrainPositions();
+ * integrator.addComputePerDof("v", "v+0.5*dt*f/m+(x-x1)/dt");
+ * integrator.addConstrainVelocities();
+ * </c++>
+ * <python>
+ * integrator = CustomIntegrator(0.001)
+ * integrator.addPerDofVariable("x1", 0)
+ * integrator.addUpdateContextState()
+ * integrator.addComputePerDof("v", "v+0.5*dt*f/m")
+ * integrator.addComputePerDof("x", "x+dt*v")
+ * integrator.addComputePerDof("x1", "x")
+ * integrator.addConstrainPositions()
+ * integrator.addComputePerDof("v", "v+0.5*dt*f/m+(x-x1)/dt")
+ * integrator.addConstrainVelocities()
+ * </python>
  * \endverbatim
  *
  * CustomIntegrator can be used to implement multiple time step integrators.  The
@@ -175,17 +190,25 @@ namespace OpenMM {
  * It evaluates the "fast" forces four times as often as the "slow" forces.
  *
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    CustomIntegrator integrator(0.004);
- *    integrator.addComputePerDof("v", "v+0.5*dt*f1/m");
- *    for (int i = 0; i &lt; 4; i++) {
- *        integrator.addComputePerDof("v", "v+0.5*(dt/4)*f0/m");
- *        integrator.addComputePerDof("x", "x+(dt/4)*v");
- *        integrator.addComputePerDof("v", "v+0.5*(dt/4)*f0/m");
- *    }
- *    integrator.addComputePerDof("v", "v+0.5*dt*f1/m");
- *
+ * <c++>
+ * CustomIntegrator integrator(0.004);
+ * integrator.addComputePerDof("v", "v+0.5*dt*f1/m");
+ * for (int i = 0; i < 4; i++) {
+ *     integrator.addComputePerDof("v", "v+0.5*(dt/4)*f0/m");
+ *     integrator.addComputePerDof("x", "x+(dt/4)*v");
+ *     integrator.addComputePerDof("v", "v+0.5*(dt/4)*f0/m");
+ * }
+ * integrator.addComputePerDof("v", "v+0.5*dt*f1/m");
+ * </c++>
+ * <python>
+ * integrator = CustomIntegrator(0.004)
+ * integrator.addComputePerDof("v", "v+0.5*dt*f1/m")
+ * for i in range(4):
+ *     integrator.addComputePerDof("v", "v+0.5*(dt/4)*f0/m")
+ *     integrator.addComputePerDof("x", "x+(dt/4)*v")
+ *     integrator.addComputePerDof("v", "v+0.5*(dt/4)*f0/m")
+ * integrator.addComputePerDof("v", "v+0.5*dt*f1/m")
+ * </python>
  * \endverbatim
  *
  * The sequence of computations in a CustomIntegrator can include flow control in
@@ -203,12 +226,16 @@ namespace OpenMM {
  *
  *
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    integrator.beginIfBlock("uniform < acceptanceProbability");
- *    integrator.addComputePerDof("x", "xnew");
- *    integrator.endBlock();
- *
+ * <c++>
+ * integrator.beginIfBlock("uniform < acceptanceProbability");
+ * integrator.addComputePerDof("x", "xnew");
+ * integrator.endBlock();
+ * </c++>
+ * <python>
+ * integrator.beginIfBlock("uniform < acceptanceProbability")
+ * integrator.addComputePerDof("x", "xnew")
+ * integrator.endBlock()
+ * </python>
  * \endverbatim
  *
  * The condition in an "if" or "while" block is evaluated globally, so it may
@@ -228,10 +255,12 @@ namespace OpenMM {
  * particle and stores it into a per-DOF variable.
  * 
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    integrator.addComputePerDof("angularMomentum", "m*cross(x, v)");
- *
+ * <c++>
+ * integrator.addComputePerDof("angularMomentum", "m*cross(x, v)");
+ * </c++>
+ * <python>
+ * integrator.addComputePerDof("angularMomentum", "m*cross(x, v)")
+ * </python>
  * \endverbatim
  * 
  * Here are two more examples that may be useful as starting points for writing
@@ -241,40 +270,64 @@ namespace OpenMM {
  * constraints once in each time step.
  * 
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    CustomIntegrator integrator(dt);
- *    integrator.addPerDofVariable("x0", 0);
- *    integrator.addUpdateContextState();
- *    integrator.addComputePerDof("x0", "x");
- *    integrator.addComputePerDof("v", "v+dt*f/m");
- *    integrator.addComputePerDof("x", "x+dt*v");
- *    integrator.addConstrainPositions();
- *    integrator.addComputePerDof("v", "(x-x0)/dt");
- *
+ * <c++>
+ * CustomIntegrator integrator(dt);
+ * integrator.addPerDofVariable("x0", 0);
+ * integrator.addUpdateContextState();
+ * integrator.addComputePerDof("x0", "x");
+ * integrator.addComputePerDof("v", "v+dt*f/m");
+ * integrator.addComputePerDof("x", "x+dt*v");
+ * integrator.addConstrainPositions();
+ * integrator.addComputePerDof("v", "(x-x0)/dt");
+ * </c++>
+ * <python>
+ * integrator = CustomIntegrator(dt)
+ * integrator.addPerDofVariable("x0", 0)
+ * integrator.addUpdateContextState()
+ * integrator.addComputePerDof("x0", "x")
+ * integrator.addComputePerDof("v", "v+dt*f/m")
+ * integrator.addComputePerDof("x", "x+dt*v")
+ * integrator.addConstrainPositions()
+ * integrator.addComputePerDof("v", "(x-x0)/dt")
+ * </python>
  * \endverbatim
  * 
  * The second one implements the algorithm used by the standard
  * LangevinMiddleIntegrator class.  kB is Boltzmann's constant.
  * 
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    CustomIntegrator integrator(dt);
- *    integrator.addGlobalVariable("a", exp(-friction*dt));
- *    integrator.addGlobalVariable("b", sqrt(1-exp(-2*friction*dt)));
- *    integrator.addGlobalVariable("kT", kB*temperature);
- *    integrator.addPerDofVariable("x1", 0);
- *    integrator.addUpdateContextState();
- *    integrator.addComputePerDof("v", "v + dt*f/m");
- *    integrator.addConstrainVelocities();
- *    integrator.addComputePerDof("x", "x + 0.5*dt*v");
- *    integrator.addComputePerDof("v", "a*v + b*sqrt(kT/m)*gaussian");
- *    integrator.addComputePerDof("x", "x + 0.5*dt*v");
- *    integrator.addComputePerDof("x1", "x");
- *    integrator.addConstrainPositions();
- *    integrator.addComputePerDof("v", "v + (x-x1)/dt");
- *
+ * <c++>
+ * CustomIntegrator integrator(dt);
+ * integrator.addGlobalVariable("a", exp(-friction*dt));
+ * integrator.addGlobalVariable("b", sqrt(1-exp(-2*friction*dt)));
+ * integrator.addGlobalVariable("kT", kB*temperature);
+ * integrator.addPerDofVariable("x1", 0);
+ * integrator.addUpdateContextState();
+ * integrator.addComputePerDof("v", "v + dt*f/m");
+ * integrator.addConstrainVelocities();
+ * integrator.addComputePerDof("x", "x + 0.5*dt*v");
+ * integrator.addComputePerDof("v", "a*v + b*sqrt(kT/m)*gaussian");
+ * integrator.addComputePerDof("x", "x + 0.5*dt*v");
+ * integrator.addComputePerDof("x1", "x");
+ * integrator.addConstrainPositions();
+ * integrator.addComputePerDof("v", "v + (x-x1)/dt");
+ * </c++>
+ * <python>
+ * integrator = CustomIntegrator(dt)
+ * integrator.addGlobalVariable("a", exp(-friction*dt))
+ * integrator.addGlobalVariable("b", sqrt(1-exp(-2*friction*dt)))
+ * integrator.addGlobalVariable("kT", kB*temperature)
+ * integrator.addPerDofVariable("x1", 0)
+ * integrator.addUpdateContextState()
+ * integrator.addComputePerDof("v", "v + dt*f/m")
+ * integrator.addConstrainVelocities()
+ * integrator.addComputePerDof("x", "x + 0.5*dt*v")
+ * integrator.addComputePerDof("v", "a*v + b*sqrt(kT/m)*gaussian")
+ * integrator.addComputePerDof("x", "x + 0.5*dt*v")
+ * integrator.addComputePerDof("x1", "x")
+ * integrator.addConstrainPositions()
+ * integrator.addComputePerDof("v", "v + (x-x1)/dt")
+ * </python>
  * \endverbatim
  * 
  * Another feature of CustomIntegrator is that it can use derivatives of the
@@ -303,10 +356,12 @@ namespace OpenMM {
  * when using a leapfrog algorithm:
  *
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    integrator.setKineticEnergyExpression("m*v1*v1/2; v1=v+0.5*dt*f/m");
- *
+ * <c++>
+ * integrator.setKineticEnergyExpression("m*v1*v1/2; v1=v+0.5*dt*f/m");
+ * </c++>
+ * <python>
+ * integrator.setKineticEnergyExpression("m*v1*v1/2; v1=v+0.5*dt*f/m")
+ * </python>
  * \endverbatim
  *
  * The kinetic energy expression may depend on the following pre-defined variables:

openmmapi/include/openmm/CustomManyParticleForce.h

@@ -86,23 +86,31 @@ namespace OpenMM {
  * is an interaction between three particles that depends on all three distances and angles formed by the particles.
  *
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    CustomManyParticleForce* force = new CustomManyParticleForce(3,
- *        "C*(1+3*cos(theta1)*cos(theta2)*cos(theta3))/(r12*r13*r23)^3;"
- *        "theta1=angle(p1,p2,p3); theta2=angle(p2,p3,p1); theta3=angle(p3,p1,p2);"
- *        "r12=distance(p1,p2); r13=distance(p1,p3); r23=distance(p2,p3)");
- *    force->setPermutationMode(CustomManyParticleForce::SinglePermutation);
- *
+ * <c++>
+ * CustomManyParticleForce* force = new CustomManyParticleForce(3,
+ *     "C*(1+3*cos(theta1)*cos(theta2)*cos(theta3))/(r12*r13*r23)^3;"
+ *     "theta1=angle(p1,p2,p3); theta2=angle(p2,p3,p1); theta3=angle(p3,p1,p2);"
+ *     "r12=distance(p1,p2); r13=distance(p1,p3); r23=distance(p2,p3)");
+ * force->setPermutationMode(CustomManyParticleForce::SinglePermutation);
+ * </c++>
+ * <python>
+ * force = CustomManyParticleForce(3,
+ *     "C*(1+3*cos(theta1)*cos(theta2)*cos(theta3))/(r12*r13*r23)^3;"
+ *     "theta1=angle(p1,p2,p3); theta2=angle(p2,p3,p1); theta3=angle(p3,p1,p2);"
+ *     "r12=distance(p1,p2); r13=distance(p1,p3); r23=distance(p2,p3)")
+ * force.setPermutationMode(CustomManyParticleForce.SinglePermutation)
+ * </python>
  * \endverbatim
  *
  * This force depends on one parameter, C.  The following code defines it as a global parameter:
  *
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    force->addGlobalParameter("C", 1.0);
- *
+ * <c++>
+ * force->addGlobalParameter("C", 1.0);
+ * </c++>
+ * <python>
+ * force.addGlobalParameter("C", 1.0)
+ * </python>
  * \endverbatim
  *
  * Notice that the expression is symmetric with respect to the particles.  It only depends on the products
@@ -118,13 +126,18 @@ namespace OpenMM {
  * potential:
  *
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    CustomManyParticleForce* force = new CustomManyParticleForce(3,
- *        "L*eps*(cos(theta1)+1/3)^2*exp(sigma*gamma/(r12-a*sigma))*exp(sigma*gamma/(r13-a*sigma));"
- *        "r12 = distance(p1,p2); r13 = distance(p1,p3); theta1 = angle(p3,p1,p2)");
- *    force->setPermutationMode(CustomManyParticleForce::UniqueCentralParticle);
- *
+ * <c++>
+ * CustomManyParticleForce* force = new CustomManyParticleForce(3,
+ *     "L*eps*(cos(theta1)+1/3)^2*exp(sigma*gamma/(r12-a*sigma))*exp(sigma*gamma/(r13-a*sigma));"
+ *     "r12 = distance(p1,p2); r13 = distance(p1,p3); theta1 = angle(p3,p1,p2)");
+ * force->setPermutationMode(CustomManyParticleForce::UniqueCentralParticle);
+ * </c++>
+ * <python>
+ * force = CustomManyParticleForce(3,
+ *     "L*eps*(cos(theta1)+1/3)^2*exp(sigma*gamma/(r12-a*sigma))*exp(sigma*gamma/(r13-a*sigma));"
+ *     "r12 = distance(p1,p2); r13 = distance(p1,p3); theta1 = angle(p3,p1,p2)")
+ * force.setPermutationMode(CustomManyParticleForce.UniqueCentralParticle)
+ * </python>
  * \endverbatim
  *
  * When the permutation mode is set to UniqueCentralParticle, particle p1 is treated as the central particle.  For a set of
@@ -144,18 +157,24 @@ namespace OpenMM {
  *
  * A particle type is an integer that you specify when you call addParticle().  (If you omit the argument, it defaults
  * to 0.)  For the water model, you could specify 0 for all oxygen atoms and 1 for all hydrogen atoms.  You can then
- * call setTypeFilter() to specify the list of allowed types for each of the N particles involved in an interaction:
+ * call setTypeFilter() to specify the set of allowed types for each of the N particles involved in an interaction:
  *
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    set&lt;int&gt; oxygenTypes, hydrogenTypes;
- *    oxygenTypes.insert(0);
- *    hydrogenTypes.insert(1);
- *    force->setTypeFilter(0, oxygenTypes);
- *    force->setTypeFilter(1, hydrogenTypes);
- *    force->setTypeFilter(2, hydrogenTypes);
- *
+ * <c++>
+ * set<int> oxygenTypes, hydrogenTypes;
+ * oxygenTypes.insert(0);
+ * hydrogenTypes.insert(1);
+ * force.setTypeFilter(0, oxygenTypes);
+ * force.setTypeFilter(1, hydrogenTypes);
+ * force.setTypeFilter(2, hydrogenTypes);
+ * </c++>
+ * <python>
+ * oxygenTypes = {0}
+ * hydrogenTypes = {1}
+ * force.setTypeFilter(0, oxygenTypes)
+ * force.setTypeFilter(1, hydrogenTypes)
+ * force.setTypeFilter(2, hydrogenTypes)
+ * </python>
  * \endverbatim
  *
  * This specifies that of the three particles in an interaction, p1 must be oxygen while p2 and p3 must be hydrogen.
@@ -177,10 +196,12 @@ namespace OpenMM {
  * from particle p1 to the midpoint between particles p2 and p3.
  * 
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    CustomManyParticleForce* force = new CustomManyParticleForce(3, "pointdistance(x1, y1, z1, (x2+x3)/2, (y2+y3)/2, (z2+z3)/2)");
- *
+ * <c++>
+ * CustomManyParticleForce* force = new CustomManyParticleForce(3, "pointdistance(x1, y1, z1, (x2+x3)/2, (y2+y3)/2, (z2+z3)/2)");
+ * </c++>
+ * <python>
+ * force = CustomManyParticleForce(3, "pointdistance(x1, y1, z1, (x2+x3)/2, (y2+y3)/2, (z2+z3)/2)")
+ * </python>
  * \endverbatim
  *
  * In addition, you can call addTabulatedFunction() to define a new function based on tabulated values.  You specify the function by

openmmapi/include/openmm/CustomNonbondedForce.h

@@ -67,20 +67,25 @@ namespace OpenMM {
  * As an example, the following code creates a CustomNonbondedForce that implements a 12-6 Lennard-Jones potential:
  *
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    CustomNonbondedForce* force = new CustomNonbondedForce("4*epsilon*((sigma/r)^12-(sigma/r)^6); sigma=0.5*(sigma1+sigma2); epsilon=sqrt(epsilon1*epsilon2)");
- *
+ * <c++>
+ * CustomNonbondedForce* force = new CustomNonbondedForce("4*epsilon*((sigma/r)^12-(sigma/r)^6); sigma=0.5*(sigma1+sigma2); epsilon=sqrt(epsilon1*epsilon2)");
+ * </c++>
+ * <python>
+ * force = CustomNonbondedForce("4*epsilon*((sigma/r)^12-(sigma/r)^6); sigma=0.5*(sigma1+sigma2); epsilon=sqrt(epsilon1*epsilon2)")
+ * </python>
  * \endverbatim
  *
  * This force depends on two parameters: sigma and epsilon.  The following code defines these as per-particle parameters:
  *
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    force->addPerParticleParameter("sigma");
- *    force->addPerParticleParameter("epsilon");
- *
+ * <c++>
+ * force->addPerParticleParameter("sigma");
+ * force->addPerParticleParameter("epsilon");
+ * </c++>
+ * <python>
+ * force.addPerParticleParameter("sigma")
+ * force.addPerParticleParameter("epsilon")
+ * </python>
  * \endverbatim
  *
  * The expression <i>must</i> be symmetric with respect to the two particles.  It typically will only be evaluated once
@@ -94,15 +99,22 @@ namespace OpenMM {
  * code uses a global parameter (lambda) to interpolate between two different sigma values for each particle (sigmaA and sigmaB).
  * 
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    CustomNonbondedForce* force = new CustomNonbondedForce("4*epsilon*((sigma/r)^12-(sigma/r)^6); sigma=0.5*(sigma1+sigma2); epsilon=sqrt(epsilon1*epsilon2)");
- *    force->addComputedValue("sigma", "(1-lambda)*sigmaA + lambda*sigmaB");
- *    force->addGlobalParameter("lambda", 0);
- *    force->addPerParticleParameter("sigmaA");
- *    force->addPerParticleParameter("sigmaB");
- *    force->addPerParticleParameter("epsilon");
- *
+ * <c++>
+ * CustomNonbondedForce* force = new CustomNonbondedForce("4*epsilon*((sigma/r)^12-(sigma/r)^6); sigma=0.5*(sigma1+sigma2); epsilon=sqrt(epsilon1*epsilon2)");
+ * force->addComputedValue("sigma", "(1-lambda)*sigmaA + lambda*sigmaB");
+ * force->addGlobalParameter("lambda", 0);
+ * force->addPerParticleParameter("sigmaA");
+ * force->addPerParticleParameter("sigmaB");
+ * force->addPerParticleParameter("epsilon");
+ * </c++>
+ * <python>
+ * force = CustomNonbondedForce("4*epsilon*((sigma/r)^12-(sigma/r)^6); sigma=0.5*(sigma1+sigma2); epsilon=sqrt(epsilon1*epsilon2)")
+ * force.addComputedValue("sigma", "(1-lambda)*sigmaA + lambda*sigmaB")
+ * force.addGlobalParameter("lambda", 0)
+ * force.addPerParticleParameter("sigmaA")
+ * force.addPerParticleParameter("sigmaB")
+ * force.addPerParticleParameter("epsilon")
+ * </python>
  * \endverbatim
  *
  * You could, of course, embed the computation of sigma directly into the energy expression, but then it would need to be

openmmapi/include/openmm/CustomTorsionForce.h

@@ -57,29 +57,36 @@ namespace OpenMM {
  * As an example, the following code creates a CustomTorsionForce that implements a periodic potential:
  *
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    CustomTorsionForce* force = new CustomTorsionForce("0.5*k*(1-cos(theta-theta0))");
- *
+ * <c++>
+ * CustomTorsionForce* force = new CustomTorsionForce("0.5*k*(1-cos(theta-theta0))");
+ * </c++>
+ * <python>
+ * force = CustomTorsionForce("0.5*k*(1-cos(theta-theta0))")
+ * </python>
  * \endverbatim
  *
  * This force depends on two parameters: the spring constant k and equilibrium angle theta0.  The following code defines these parameters:
  *
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    force->addPerTorsionParameter("k");
- *    force->addPerTorsionParameter("theta0");
- *
+ * <c++>
+ * force->addPerTorsionParameter("k");
+ * force->addPerTorsionParameter("theta0");
+ * </c++>
+ * <python>
+ * force.addPerTorsionParameter("k")
+ * force.addPerTorsionParameter("theta0")
+ * </python>
  * \endverbatim
  *
  * If a harmonic restraint is desired, it is important to be careful of the domain for theta, using an idiom like this:
  *
  * \verbatim embed:rst:leading-asterisk
- * .. code-block:: cpp
- *
- *    CustomTorsionForce* force = new CustomTorsionForce("0.5*k*min(dtheta, 2*pi-dtheta)^2; dtheta = abs(theta-theta0); pi = 3.1415926535");
- *
+ * <c++>
+ * CustomTorsionForce* force = new CustomTorsionForce("0.5*k*min(dtheta, 2*pi-dtheta)^2; dtheta = abs(theta-theta0); pi = 3.1415926535");
+ * </c++>
+ * <python>
+ * force = CustomTorsionForce("0.5*k*min(dtheta, 2*pi-dtheta)^2; dtheta = abs(theta-theta0); pi = 3.1415926535")
+ * </python>
  * \endverbatim
  *
  * This class also has the ability to compute derivatives of the potential energy with respect to global parameters.