diff --git a/docs/examples/coreshellnp.py b/docs/examples/coreshellnp.py
index 430726b1..23871a46 100644
--- a/docs/examples/coreshellnp.py
+++ b/docs/examples/coreshellnp.py
@@ -91,14 +91,14 @@ def makeRecipe(stru1, stru2, datname):
     # diameter to twice the shell radius.
     recipe.add_variable(contribution.radius, 15)
     recipe.add_variable(contribution.thickness, 11)
-    recipe.constrain(contribution.psize, "2 * radius")
+    recipe.add_constraint(contribution.psize, "2 * radius")
 
     # Configure the fit variables
     # Start by configuring the scale factor and resolution factors.
     # We want the sum of the phase scale factors to be 1.
     recipe.create_new_variable("scale_CdS", 0.7)
-    recipe.constrain(generator_cds.scale, "scale_CdS")
-    recipe.constrain(generator_zns.scale, "1 - scale_CdS")
+    recipe.add_constraint(generator_cds.scale, "scale_CdS")
+    recipe.add_constraint(generator_zns.scale, "1 - scale_CdS")
     # We also want the resolution factor to be the same on each.
 
     # Vary the global scale as well.
@@ -117,7 +117,7 @@ def makeRecipe(stru1, stru2, datname):
     )
     # Since we know these have stacking disorder, constrain the B33 adps for
     # each atom type.
-    recipe.constrain("B33_1_cds", "B33_0_cds")
+    recipe.add_constraint("B33_1_cds", "B33_0_cds")
     recipe.add_variable(generator_cds.delta2, name="delta2_cds", value=5)
 
     phase_zns = generator_zns.phase
@@ -128,7 +128,7 @@ def makeRecipe(stru1, stru2, datname):
     recipe.add_variable(
         phase_zns.sgpars.xyzpars.z_1, name="z_1_zns", tag="xyz"
     )
-    recipe.constrain("B33_1_zns", "B33_0_zns")
+    recipe.add_constraint("B33_1_zns", "B33_0_zns")
     recipe.add_variable(generator_zns.delta2, name="delta2_zns", value=2.5)
 
     # Give the recipe away so it can be used!
diff --git a/docs/examples/crystalpdfall.py b/docs/examples/crystalpdfall.py
index db5fefa3..1ee8903c 100644
--- a/docs/examples/crystalpdfall.py
+++ b/docs/examples/crystalpdfall.py
@@ -113,15 +113,15 @@ def makeRecipe(
     for par in phase_ni.sgpars:
         recipe.add_variable(par, name=par.name + "_ni")
     delta2_ni = recipe.create_new_variable("delta2_ni", 2.5)
-    recipe.constrain(xgenerator_ni.delta2, delta2_ni)
-    recipe.constrain(ngenerator_ni.delta2, delta2_ni)
-    recipe.constrain(xgenerator_sini_ni.delta2, delta2_ni)
+    recipe.add_constraint(xgenerator_ni.delta2, delta2_ni)
+    recipe.add_constraint(ngenerator_ni.delta2, delta2_ni)
+    recipe.add_constraint(xgenerator_sini_ni.delta2, delta2_ni)
 
     for par in phase_si.sgpars:
         recipe.add_variable(par, name=par.name + "_si")
     delta2_si = recipe.create_new_variable("delta2_si", 2.5)
-    recipe.constrain(xgenerator_si.delta2, delta2_si)
-    recipe.constrain(xgenerator_sini_si.delta2, delta2_si)
+    recipe.add_constraint(xgenerator_si.delta2, delta2_si)
+    recipe.add_constraint(xgenerator_sini_si.delta2, delta2_si)
 
     # Now the experimental parameters
     recipe.add_variable(xgenerator_ni.scale, name="xscale_ni")
@@ -129,8 +129,8 @@ def makeRecipe(
     recipe.add_variable(ngenerator_ni.scale, name="nscale_ni")
     recipe.add_variable(xcontribution_sini.scale, 1.0, "xscale_sini")
     recipe.create_new_variable("pscale_sini_ni", 0.8)
-    recipe.constrain(xgenerator_sini_ni.scale, "pscale_sini_ni")
-    recipe.constrain(xgenerator_sini_si.scale, "1 - pscale_sini_ni")
+    recipe.add_constraint(xgenerator_sini_ni.scale, "pscale_sini_ni")
+    recipe.add_constraint(xgenerator_sini_si.scale, "1 - pscale_sini_ni")
 
     # The qdamp parameters are too correlated to vary so we fix them based on
     # previous measurements.
diff --git a/docs/examples/crystalpdftwodata.py b/docs/examples/crystalpdftwodata.py
index 601cb01b..0b0d6093 100644
--- a/docs/examples/crystalpdftwodata.py
+++ b/docs/examples/crystalpdftwodata.py
@@ -121,8 +121,8 @@ def makeRecipe(ciffile, xdatname, ndatname):
     # delta2 is a non-structual material property. Thus, we constrain together
     # delta2 Parameter from each PDFGenerator.
     delta2 = recipe.create_new_variable("delta2", 2)
-    recipe.constrain(xgenerator.delta2, delta2)
-    recipe.constrain(ngenerator.delta2, delta2)
+    recipe.add_constraint(xgenerator.delta2, delta2)
+    recipe.add_constraint(ngenerator.delta2, delta2)
 
     # We only need to constrain phase properties once since there is a single
     # ObjCrystCrystalParSet for the Crystal.
diff --git a/docs/examples/crystalpdftwophase.py b/docs/examples/crystalpdftwophase.py
index de9cd965..1d6878b2 100644
--- a/docs/examples/crystalpdftwophase.py
+++ b/docs/examples/crystalpdftwophase.py
@@ -90,12 +90,12 @@ def makeRecipe(niciffile, siciffile, datname):
     # Start by configuring the scale factor and resolution factors.
     # We want the sum of the phase scale factors to be 1.
     recipe.create_new_variable("scale_ni", 0.1)
-    recipe.constrain(generator_ni.scale, "scale_ni")
-    recipe.constrain(generator_si.scale, "1 - scale_ni")
+    recipe.add_constraint(generator_ni.scale, "scale_ni")
+    recipe.add_constraint(generator_si.scale, "1 - scale_ni")
     # We also want the resolution factor to be the same on each.
     recipe.create_new_variable("qdamp", 0.03)
-    recipe.constrain(generator_ni.qdamp, "qdamp")
-    recipe.constrain(generator_si.qdamp, "qdamp")
+    recipe.add_constraint(generator_ni.qdamp, "qdamp")
+    recipe.add_constraint(generator_si.qdamp, "qdamp")
 
     # Vary the global scale as well.
     recipe.add_variable(contribution.scale, 1)
@@ -123,18 +123,30 @@ def makeRecipe(niciffile, siciffile, datname):
     # derived has no uncertainty. Thus, we will tell the recipe to scale the
     # residual, which means that it will be weighted as much as the average
     # data point during the fit.
-    recipe.restrain("a_ni", lb=3.527, ub=3.527, scaled=True)
+    recipe.add_soft_bounds(
+        "a_ni", lower_bound=3.527, upper_bound=3.527, scaled=True
+    )
     # Now we do the same with the delta2 and Biso parameters (remember that
     # Biso = 8*pi**2*Uiso)
-    recipe.restrain("delta2_ni", lb=2.22, ub=2.22, scaled=True)
-    recipe.restrain("Biso_0_ni", lb=0.454, ub=0.454, scaled=True)
+    recipe.add_soft_bounds(
+        "delta2_ni", lower_bound=2.22, upper_bound=2.22, scaled=True
+    )
+    recipe.add_soft_bounds(
+        "Biso_0_ni", lower_bound=0.454, upper_bound=0.454, scaled=True
+    )
     #
     # We can do the same with the silicon values. We haven't done a thorough
     # job of measuring the uncertainties in the results, so we'll scale these
     # as well.
-    recipe.restrain("a_si", lb=5.430, ub=5.430, scaled=True)
-    recipe.restrain("delta2_si", lb=3.54, ub=3.54, scaled=True)
-    recipe.restrain("Biso_0_si", lb=0.645, ub=0.645, scaled=True)
+    recipe.add_soft_bounds(
+        "a_si", lower_bound=5.430, upper_bound=5.430, scaled=True
+    )
+    recipe.add_soft_bounds(
+        "delta2_si", lower_bound=3.54, upper_bound=3.54, scaled=True
+    )
+    recipe.add_soft_bounds(
+        "Biso_0_si", lower_bound=0.645, upper_bound=0.645, scaled=True
+    )
 
     # Give the recipe away so it can be used!
     return recipe
diff --git a/docs/examples/debyemodel.py b/docs/examples/debyemodel.py
index cd2dcbbd..7813be4e 100644
--- a/docs/examples/debyemodel.py
+++ b/docs/examples/debyemodel.py
@@ -152,7 +152,7 @@ def makeRecipe():
     # breaking the restraint by the point-average chi^2 value so that the
     # restraint is roughly as significant as any other data point throughout
     # the fit.
-    recipe.restrain(recipe.offset, lb=0, scaled=True)
+    recipe.add_soft_bounds(recipe.offset, lower_bound=0, scaled=True)
 
     # We're done setting up the recipe. We can now do other things with it.
     return recipe
diff --git a/docs/examples/debyemodelII.py b/docs/examples/debyemodelII.py
index 4817d11d..a79df5f6 100644
--- a/docs/examples/debyemodelII.py
+++ b/docs/examples/debyemodelII.py
@@ -89,8 +89,8 @@ def makeRecipeII():
     # We create a new Variable and use the recipe's "constrain" method to
     # associate the Debye temperature parameters with that variable.
     recipe.create_new_variable("thetaD", 100)
-    recipe.constrain(recipe.lowT.thetaD, "thetaD")
-    recipe.constrain(recipe.highT.thetaD, "thetaD")
+    recipe.add_constraint(recipe.lowT.thetaD, "thetaD")
+    recipe.add_constraint(recipe.highT.thetaD, "thetaD")
     return recipe
 
 
diff --git a/docs/examples/npintensity.py b/docs/examples/npintensity.py
index b9d33075..eaf17012 100644
--- a/docs/examples/npintensity.py
+++ b/docs/examples/npintensity.py
@@ -290,15 +290,15 @@ def gaussian(q, q0, width):
     lattice = phase.getLattice()
     a = lattice.a
     recipe.add_variable(a)
-    recipe.constrain(lattice.b, a)
-    recipe.constrain(lattice.c, a)
+    recipe.add_constraint(lattice.b, a)
+    recipe.add_constraint(lattice.c, a)
     # We want to refine the thermal parameters as well. We will add a new
     # Variable that we call "Uiso" and constrain the atomic Uiso values to
     # this. Note that we don't give Uiso an initial value. The initial value
     # will be inferred from the following constraints.
     Uiso = recipe.create_new_variable("Uiso")
     for atom in phase.getScatterers():
-        recipe.constrain(atom.Uiso, Uiso)
+        recipe.add_constraint(atom.Uiso, Uiso)
 
     # Give the recipe away so it can be used!
     return recipe
diff --git a/docs/examples/npintensityII.py b/docs/examples/npintensityII.py
index 92e5e899..5a6415a6 100644
--- a/docs/examples/npintensityII.py
+++ b/docs/examples/npintensityII.py
@@ -173,15 +173,15 @@ def gaussian(q, q0, width):
     a = recipe.add_variable(lattice.a)
     # We want to allow for isotropic expansion, so we'll make constraints for
     # that.
-    recipe.constrain(lattice.b, a)
-    recipe.constrain(lattice.c, a)
+    recipe.add_constraint(lattice.b, a)
+    recipe.add_constraint(lattice.c, a)
     # We want to refine the thermal parameters as well. We will add a new
     # variable that we call "Uiso" and constrain the atomic Uiso values to
     # this. Note that we don't give Uiso an initial value. The initial value
     # will be inferred from the subsequent constraints.
     Uiso = recipe.create_new_variable("Uiso")
     for atom in phase.getScatterers():
-        recipe.constrain(atom.Uiso, Uiso)
+        recipe.add_constraint(atom.Uiso, Uiso)
 
     # Give the recipe away so it can be used!
     return recipe
diff --git a/docs/examples/nppdfobjcryst.py b/docs/examples/nppdfobjcryst.py
index cce39962..03511d19 100644
--- a/docs/examples/nppdfobjcryst.py
+++ b/docs/examples/nppdfobjcryst.py
@@ -73,7 +73,7 @@ def makeRecipe(molecule, datname):
         # has no scattering power. It is only used as a reference point for
         # our bond length. We don't want to constrain it.
         if not atom.isDummy():
-            recipe.constrain(atom.Biso, Biso)
+            recipe.add_constraint(atom.Biso, Biso)
 
     # We need to let the molecule expand. If we were modeling it as a crystal,
     # we could let the unit cell expand. For instruction purposes, we use a
@@ -97,7 +97,7 @@ def makeRecipe(molecule, datname):
         # This creates a Parameter that moves the second atom according to the
         # bond length. Note that each Parameter needs a unique name.
         par = c60.addBondLengthParameter("rad%i" % i, center, atom)
-        recipe.constrain(par, radius)
+        recipe.add_constraint(par, radius)
 
     # Add the correlation term, scale. The scale is too short to effectively
     # determine qdamp.
diff --git a/docs/examples/nppdfsas.py b/docs/examples/nppdfsas.py
index 69386a8b..7bea8069 100644
--- a/docs/examples/nppdfsas.py
+++ b/docs/examples/nppdfsas.py
@@ -113,8 +113,8 @@ def makeRecipe(ciffile, grdata, iqdata):
     # Even though the cfcalculator and sasgenerator depend on the same sas
     # model, we must still constrain the cfcalculator Parameters so that it is
     # informed of changes in the refined parameters.
-    recipe.constrain(cfcalculator.radius_a, "radius_a")
-    recipe.constrain(cfcalculator.radius_b, "radius_b")
+    recipe.add_constraint(cfcalculator.radius_a, "radius_a")
+    recipe.add_constraint(cfcalculator.radius_b, "radius_b")
 
     return recipe
 
@@ -126,9 +126,9 @@ def fitRecipe(recipe):
     recipe.set_weight(recipe.pdf, 0)
     recipe.fix("all")
     recipe.free("radius_a", "radius_b", iqscale=1e8)
-    recipe.constrain("radius_b", "radius_a")
+    recipe.add_constraint("radius_b", "radius_a")
     scipyOptimize(recipe)
-    recipe.unconstrain("radius_b")
+    recipe.remove_constraint("radius_b")
 
     # Tune PDF
     recipe.set_weight(recipe.pdf, 1)
diff --git a/docs/examples/simplepdftwophase.py b/docs/examples/simplepdftwophase.py
index 7b252baf..0b7ee046 100644
--- a/docs/examples/simplepdftwophase.py
+++ b/docs/examples/simplepdftwophase.py
@@ -47,8 +47,8 @@ def makeRecipe(niciffile, siciffile, datname):
     # Start by configuring the scale factor and resolution factors.
     # We want the sum of the phase scale factors to be 1.
     recipe.create_new_variable("scale_ni", 0.1)
-    recipe.constrain(contribution.ni.scale, "scale_ni")
-    recipe.constrain(contribution.si.scale, "1 - scale_ni")
+    recipe.add_constraint(contribution.ni.scale, "scale_ni")
+    recipe.add_constraint(contribution.si.scale, "1 - scale_ni")
     # We also want the resolution factor to be the same on each. This is done
     # for free by the PDFContribution. We simply need to add it to the recipe.
     recipe.add_variable(contribution.qdamp, 0.03)
@@ -82,18 +82,30 @@ def makeRecipe(niciffile, siciffile, datname):
     # derived has no uncertainty. Thus, we will tell the recipe to scale the
     # residual, which means that it will be weighted as much as the average
     # data point during the fit.
-    recipe.restrain("a_ni", lb=3.527, ub=3.527, scaled=True)
+    recipe.add_soft_bounds(
+        "a_ni", lower_bound=3.527, upper_bound=3.527, scaled=True
+    )
     # Now we do the same with the delta2 and Biso parameters (remember that
     # Biso = 8*pi**2*Uiso)
-    recipe.restrain("delta2_ni", lb=2.22, ub=2.22, scaled=True)
-    recipe.restrain("Biso_0_ni", lb=0.454, ub=0.454, scaled=True)
+    recipe.add_soft_bounds(
+        "delta2_ni", lower_bound=2.22, upper_bound=2.22, scaled=True
+    )
+    recipe.add_soft_bounds(
+        "Biso_0_ni", lower_bound=0.454, upper_bound=0.454, scaled=True
+    )
     #
     # We can do the same with the silicon values. We haven't done a thorough
     # job of measuring the uncertainties in the results, so we'll scale these
     # as well.
-    recipe.restrain("a_si", lb=5.430, ub=5.430, scaled=True)
-    recipe.restrain("delta2_si", lb=3.54, ub=3.54, scaled=True)
-    recipe.restrain("Biso_0_si", lb=0.645, ub=0.645, scaled=True)
+    recipe.add_soft_bounds(
+        "a_si", lower_bound=5.430, upper_bound=5.430, scaled=True
+    )
+    recipe.add_soft_bounds(
+        "delta2_si", lower_bound=3.54, upper_bound=3.54, scaled=True
+    )
+    recipe.add_soft_bounds(
+        "Biso_0_si", lower_bound=0.645, upper_bound=0.645, scaled=True
+    )
 
     # Give the recipe away so it can be used!
     return recipe
diff --git a/docs/examples/threedoublepeaks.py b/docs/examples/threedoublepeaks.py
index 4a3d7eb6..ecf01bf4 100644
--- a/docs/examples/threedoublepeaks.py
+++ b/docs/examples/threedoublepeaks.py
@@ -121,9 +121,9 @@ def peakloc(mu):
         return 180 / pi * arcsin(pi / 180 * l2 * sin(mu) / l1)
 
     recipe.register_function(peakloc)
-    recipe.constrain(contribution.mu12, "peakloc(mu11)")
-    recipe.constrain(contribution.mu22, "peakloc(mu21)")
-    recipe.constrain(contribution.mu32, "peakloc(mu31)")
+    recipe.add_constraint(contribution.mu12, "peakloc(mu11)")
+    recipe.add_constraint(contribution.mu22, "peakloc(mu21)")
+    recipe.add_constraint(contribution.mu32, "peakloc(mu31)")
 
     # Vary the width of the peaks. We know the functional form of the peak
     # broadening.
@@ -139,20 +139,20 @@ def sig(sig0, dsig, mu):
     # Now constrain the peak widths to this
     recipe.sig0.value = 0.001
     recipe.dsig.value = 4.0
-    recipe.constrain(contribution.sig11, "sig(sig0, dsig, mu11)")
-    recipe.constrain(
+    recipe.add_constraint(contribution.sig11, "sig(sig0, dsig, mu11)")
+    recipe.add_constraint(
         contribution.sig12,
         "sig(sig0, dsig, mu12)",
         ns={"mu12": contribution.mu12},
     )
-    recipe.constrain(contribution.sig21, "sig(sig0, dsig, mu21)")
-    recipe.constrain(
+    recipe.add_constraint(contribution.sig21, "sig(sig0, dsig, mu21)")
+    recipe.add_constraint(
         contribution.sig22,
         "sig(sig0, dsig, mu22)",
         ns={"mu22": contribution.mu22},
     )
-    recipe.constrain(contribution.sig31, "sig(sig0, dsig, mu31)")
-    recipe.constrain(
+    recipe.add_constraint(contribution.sig31, "sig(sig0, dsig, mu31)")
+    recipe.add_constraint(
         contribution.sig32,
         "sig(sig0, dsig, mu32)",
         ns={"mu32": contribution.mu32},
diff --git a/news/recipeorg-dep2.rst b/news/recipeorg-dep2.rst
new file mode 100644
index 00000000..bcd8a3b0
--- /dev/null
+++ b/news/recipeorg-dep2.rst
@@ -0,0 +1,35 @@
+**Added:**
+
+* Added ``add_constraint`` method to ``RecipeOrganizer``.
+* Added ``remove_constraint`` method to ``RecipeOrganizer``.
+* Added ``add_soft_bounds`` method to ``RecipeOrganizer``.
+* Added ``remove_soft_bounds`` method to ``RecipeOrganizer``.
+* Added ``register_soft_bounds`` method to ``RecipeOrganizer``.
+* Added ``clear_all_soft_bounds`` method to ``RecipeOrganizer``.
+* Added ``get_equation_from_string`` method to ``RecipeOrganizer``.
+
+**Changed:**
+
+* <news item>
+
+**Deprecated:**
+
+* Deprecated ``constrain`` method of ``RecipeOrganizer``. Use ``add_constraint`` instead.
+* Deprecated ``unconstrain`` method of ``RecipeOrganizer``. Use ``remove_constraint`` instead.
+* Deprecated ``restrain`` method of ``RecipeOrganizer``. Use ``add_soft_bounds`` instead.
+* Deprecated ``unrestrain`` methods of ``RecipeOrganizer``. Use ``remove_soft_bounds`` instead.
+* Deprecated ``addRestraint`` method of ``RecipeOrganizer``. Use ``register_soft_bounds`` instead.
+* Deprecate ``clearRestraints`` method of ``RecipeOrganizer``. Use ``clear_all_soft_bounds`` instead.
+* Deprecated ``equationFromString`` method of ``RecipeOrganizer``. Use ``get_equation_from_string`` instead.
+
+**Removed:**
+
+* <news item>
+
+**Fixed:**
+
+* <news item>
+
+**Security:**
+
+* <news item>
diff --git a/src/diffpy/srfit/fitbase/constraint.py b/src/diffpy/srfit/fitbase/constraint.py
index b0a89d19..dab6002b 100644
--- a/src/diffpy/srfit/fitbase/constraint.py
+++ b/src/diffpy/srfit/fitbase/constraint.py
@@ -24,6 +24,24 @@
 
 from diffpy.srfit.exceptions import SrFitError
 from diffpy.srfit.fitbase.validatable import Validatable
+from diffpy.utils._deprecator import build_deprecation_message, deprecated
+
+base = "diffpy.srfit.fitbase.constraint.Constraint"
+removal_version = "4.0.0"
+
+constrain_deprecation_msg = build_deprecation_message(
+    base,
+    "constrain",
+    "add_constraint",
+    removal_version,
+)
+
+unconstrain_deprecation_msg = build_deprecation_message(
+    base,
+    "unconstrain",
+    "remove_constraint",
+    removal_version,
+)
 
 
 class Constraint(Validatable):
@@ -47,13 +65,23 @@ def __init__(self):
         self.eq = None
         return
 
-    def constrain(self, par, eq):
+    def add_constraint(self, par, eq):
         """Constrain a Parameter according to an Equation.
 
         The parameter will be set constant once it is constrained. This
         will keep it from being constrained multiple times.
 
-        Raises a ValueError if par is const.
+        Parameters
+        ----------
+        par : Parameter
+            The Parameter to constrain.
+        eq : Equation
+            The Equation to use to constrain the Parameter.
+
+        Raises
+        ------
+        ValueError
+            If par is constant or already constrained.
         """
 
         if par.const:
@@ -69,13 +97,37 @@ def constrain(self, par, eq):
         self.update()
         return
 
-    def unconstrain(self):
-        """Clear the constraint."""
+    @deprecated(constrain_deprecation_msg)
+    def constrain(self, par, eq):
+        """This function has been deprecated and will be removed in
+        version 4.0.0.
+
+        Please use
+        diffpy.srfit.fitbase.constraint.Constraint.add_constraint
+        instead.
+        """
+        self.add_constraint(par, eq)
+        return
+
+    def remove_constraint(self):
+        """Clear the constraint from a Parameter."""
         self.par.constrained = False
         self.par = None
         self.eq = None
         return
 
+    @deprecated(unconstrain_deprecation_msg)
+    def unconstrain(self):
+        """This function has been deprecated and will be removed in
+        version 4.0.0.
+
+        Please use
+        diffpy.srfit.fitbase.constraint.Constraint.remove_constraint
+        instead.
+        """
+        self.remove_constraint()
+        return
+
     def update(self):
         """Update the parameter according to the equation."""
         # This will be evaluated quickly thanks to the Equation class.
diff --git a/src/diffpy/srfit/fitbase/fitcontribution.py b/src/diffpy/srfit/fitbase/fitcontribution.py
index d5af4709..feda75ba 100644
--- a/src/diffpy/srfit/fitbase/fitcontribution.py
+++ b/src/diffpy/srfit/fitbase/fitcontribution.py
@@ -28,7 +28,7 @@
 from diffpy.srfit.fitbase.parameter import ParameterProxy
 from diffpy.srfit.fitbase.parameterset import ParameterSet
 from diffpy.srfit.fitbase.profile import Profile
-from diffpy.srfit.fitbase.recipeorganizer import equationFromString
+from diffpy.srfit.fitbase.recipeorganizer import get_equation_from_string
 from diffpy.utils._deprecator import build_deprecation_message, deprecated
 
 base = "diffpy.srfit.fitbase.FitContribution"
@@ -296,7 +296,9 @@ def set_equation(self, eqstr, ns={}):
         variable.
         """
         # Build the equation instance.
-        eq = equationFromString(eqstr, self._eqfactory, buildargs=True, ns=ns)
+        eq = get_equation_from_string(
+            eqstr, self._eqfactory, buildargs=True, ns=ns
+        )
         eq.name = "eq"
 
         # Register any new Parameters.
@@ -386,7 +388,7 @@ def set_residual_equation(self, eqstr):
         elif eqstr == "resv":
             eqstr = resvstr
 
-        reseq = equationFromString(eqstr, self._eqfactory)
+        reseq = get_equation_from_string(eqstr, self._eqfactory)
         self._eqfactory.wipeout(self._reseq)
         self._reseq = reseq
 
diff --git a/src/diffpy/srfit/fitbase/fitrecipe.py b/src/diffpy/srfit/fitbase/fitrecipe.py
index 95021384..0961b504 100644
--- a/src/diffpy/srfit/fitbase/fitrecipe.py
+++ b/src/diffpy/srfit/fitbase/fitrecipe.py
@@ -125,6 +125,14 @@
     base, "boundsToRestraints", "convert_bounds_to_restraints", removal_version
 )
 
+constrain_dep_msg = build_deprecation_message(
+    base, "constrain", "add_constraint", removal_version
+)
+
+unconstrain_dep_msg = build_deprecation_message(
+    base, "unconstrain", "remove_constraint", removal_version
+)
+
 
 class FitRecipe(_fitrecipe_interface, RecipeOrganizer):
     """FitRecipe class.
@@ -472,6 +480,11 @@ def removeParameterSet(self, parset):
     def residual(self, p=[]):
         """Calculate the vector residual to be optimized.
 
+        The residual is by default the weighted concatenation of each
+        FitContribution's residual, plus the value of each restraint. The array
+        returned, denoted chiv, is such that
+        dot(chiv, chiv) = chi^2 + restraints.
+
         Parameters
         ----------
         p : list or numpy.ndarray
@@ -481,10 +494,11 @@ def residual(self, p=[]):
             updated in some other way, and the explicit update within this
             function is skipped.
 
-        The residual is by default the weighted concatenation of each
-        FitContribution's residual, plus the value of each restraint. The array
-        returned, denoted chiv, is such that
-        dot(chiv, chiv) = chi^2 + restraints.
+        Return
+        ------
+        chiv : numpy.ndarray
+            The array of residuals to be optimized. The array is such that
+            dot(chiv, chiv) = chi^2 + restraints.
         """
 
         # Prepare, if necessary
@@ -1093,7 +1107,7 @@ def isFree(self, var):
         """
         return self.is_free(var)
 
-    def unconstrain(self, *pars):
+    def remove_constraint(self, *pars):
         """Unconstrain a Parameter.
 
         This removes any constraints on a Parameter. If the Parameter is also a
@@ -1119,7 +1133,7 @@ def unconstrain(self, *pars):
                 raise ValueError("The parameter cannot be found")
 
             if par in self._constraints:
-                self._constraints[par].unconstrain()
+                self._constraints[par].remove_constraint()
                 del self._constraints[par]
                 update = True
 
@@ -1132,7 +1146,18 @@ def unconstrain(self, *pars):
 
         return
 
-    def constrain(self, par, con, ns={}):
+    @deprecated(unconstrain_dep_msg)
+    def unconstrain(self, *pars):
+        """This function has been deprecated and will be removed in
+        version 4.0.0.
+
+        Please use diffpy.srfit.fitbase.FitRecipe.remove_constraint
+        instead.
+        """
+        self.remove_constraint(*pars)
+        return
+
+    def add_constraint(self, par, con, ns={}):
         """Constrain a parameter to an equation.
 
         Note that only one constraint can exist on a Parameter at a time.
@@ -1190,7 +1215,18 @@ def constrain(self, par, con, ns={}):
         if par in self._parameters.values():
             self.fix(par)
 
-        RecipeOrganizer.constrain(self, par, con, ns)
+        RecipeOrganizer.add_constraint(self, par, con, ns)
+        return
+
+    @deprecated(constrain_dep_msg)
+    def constrain(self, par, con, ns={}):
+        """This function has been deprecated and will be removed in
+        version 4.0.0.
+
+        Please use diffpy.srfit.fitbase.FitRecipe.add_constraint
+        instead.
+        """
+        self.add_constraint(par, con, ns)
         return
 
     def get_values(self):
@@ -1720,7 +1756,7 @@ def convert_bounds_to_restraints(self, sig=1, scaled=False):
         if not hasattr(sig, "__iter__"):
             sig = [sig] * len(pars)
         for par, x in zip(pars, sig):
-            self.restrain(
+            self.add_soft_bounds(
                 par, par.bounds[0], par.bounds[1], sig=x, scaled=scaled
             )
         return
diff --git a/src/diffpy/srfit/fitbase/parameter.py b/src/diffpy/srfit/fitbase/parameter.py
index 1c7450d7..6bddf43b 100644
--- a/src/diffpy/srfit/fitbase/parameter.py
+++ b/src/diffpy/srfit/fitbase/parameter.py
@@ -95,10 +95,10 @@ def set_value(self, val):
         ----------
         val
             The value to assign.
-        lb
+        lower_bound : float
             The lower bounds for the bounds list. If this is None
             (default), then the lower bound will not be alterered.
-        ub
+        upper_bound : float
             The upper bounds for the bounds list. If this is None
             (default), then the upper bound will not be alterered.
 
@@ -142,14 +142,14 @@ def setConst(self, const=True, value=None):
             self.set_value(value)
         return self
 
-    def boundRange(self, lb=None, ub=None):
+    def boundRange(self, lower_bound=None, upper_bound=None):
         """Set lower and upper bound of the Parameter.
 
         Attributes
         ----------
-        lb
+        lower_bound : float
             The lower bound for the bounds list.
-        ub
+        upper_bound : float
             The upper bound for the bounds list.
 
         Returns
@@ -157,10 +157,10 @@ def boundRange(self, lb=None, ub=None):
         self
             Returns self so that mutators can be chained.
         """
-        if lb is not None:
-            self.bounds[0] = lb
-        if ub is not None:
-            self.bounds[1] = ub
+        if lower_bound is not None:
+            self.bounds[0] = lower_bound
+        if upper_bound is not None:
+            self.bounds[1] = upper_bound
         return self
 
     def boundWindow(self, lr=0, ur=None):
@@ -182,11 +182,11 @@ def boundWindow(self, lr=0, ur=None):
             Returns self so that mutators can be chained.
         """
         val = self.getValue()
-        lb = val - lr
+        lower_bound = val - lr
         if ur is None:
             ur = lr
-        ub = val + ur
-        self.bounds = [lb, ub]
+        upper_bound = val + ur
+        self.bounds = [lower_bound, upper_bound]
         return self
 
     def _validate(self):
@@ -283,8 +283,8 @@ def setConst(self, const=True, value=None):
         return self.par.setConst(const, value)
 
     @wraps(Parameter.boundRange)
-    def boundRange(self, lb=None, ub=None):
-        return self.par.boundRange(lb, ub)
+    def boundRange(self, lower_bound=None, upper_bound=None):
+        return self.par.boundRange(lower_bound, upper_bound)
 
     @wraps(Parameter.boundWindow)
     def boundWindow(self, lr=0, ur=None):
diff --git a/src/diffpy/srfit/fitbase/recipeorganizer.py b/src/diffpy/srfit/fitbase/recipeorganizer.py
index 495b3706..a3d442f7 100644
--- a/src/diffpy/srfit/fitbase/recipeorganizer.py
+++ b/src/diffpy/srfit/fitbase/recipeorganizer.py
@@ -17,10 +17,10 @@
 RecipeContainer is the base class for organizing Parameters, and other
 RecipeContainers.  RecipeOrganizer is an extended RecipeContainer that
 incorporates equation building, constraints and Restraints.
-equationFromString creates an Equation instance from a string.
+get_equation_from_string creates an Equation instance from a string.
 """
 
-__all__ = ["RecipeContainer", "RecipeOrganizer", "equationFromString"]
+__all__ = ["RecipeContainer", "RecipeOrganizer", "get_equation_from_string"]
 
 import re
 from collections import OrderedDict
@@ -118,6 +118,55 @@
     removal_version,
 )
 
+addRestraint_deprecation_msg = build_deprecation_message(
+    recipeorganizer_base,
+    "addRestraint",
+    "register_soft_bounds",
+    removal_version,
+)
+
+constrain_deprecation_msg = build_deprecation_message(
+    recipeorganizer_base,
+    "constrain",
+    "add_constraint",
+    removal_version,
+)
+
+unconstrain_deprecation_msg = build_deprecation_message(
+    recipeorganizer_base,
+    "unconstrain",
+    "remove_constraint",
+    removal_version,
+)
+
+restrain_deprecation_msg = build_deprecation_message(
+    recipeorganizer_base,
+    "restrain",
+    "add_soft_bounds",
+    removal_version,
+)
+
+unrestrain_deprecation_msg = build_deprecation_message(
+    recipeorganizer_base,
+    "unrestrain",
+    "remove_soft_bounds",
+    removal_version,
+)
+
+clearRestraints_deprecation_msg = build_deprecation_message(
+    recipeorganizer_base,
+    "clearRestraints",
+    "clear_all_soft_bounds",
+    removal_version,
+)
+
+equationFromString_deprecation_msg = build_deprecation_message(
+    recipeorganizer_base,
+    "equationFromString",
+    "get_equation_from_string",
+    removal_version,
+)
+
 
 class RecipeContainer(Observable, Configurable, Validatable):
     """Base class for organizing pieces of a FitRecipe.
@@ -809,7 +858,7 @@ def register_string_function(self, function_str, name, func_params={}):
         """
 
         # Build the equation instance.
-        eq = equationFromString(
+        eq = get_equation_from_string(
             function_str, self._eqfactory, ns=func_params, buildargs=True
         )
         eq.name = name
@@ -865,7 +914,9 @@ def evaluate_equation(self, equation_str, func_params={}):
             If `func_params` uses a name that is already used for a
             variable.
         """
-        eq = equationFromString(equation_str, self._eqfactory, func_params)
+        eq = get_equation_from_string(
+            equation_str, self._eqfactory, func_params
+        )
         try:
             returned_value = eq()
         finally:
@@ -883,7 +934,7 @@ def evaluateEquation(self, eqstr, ns={}):
         """
         return self.evaluate_equation(eqstr, func_params=ns)
 
-    def constrain(self, parameter, constraint_eq, params={}):
+    def add_constraint(self, parameter, constraint_eq, params={}):
         """Constrain a parameter to an equation.
 
         Note that only one constraint can exist on a Parameter at a time.
@@ -927,7 +978,9 @@ def constrain(self, parameter, constraint_eq, params={}):
 
         if isinstance(constraint_eq, str):
             eqstr = constraint_eq
-            eq = equationFromString(constraint_eq, self._eqfactory, params)
+            eq = get_equation_from_string(
+                constraint_eq, self._eqfactory, params
+            )
         else:
             eq = Equation(root=constraint_eq)
             eqstr = constraint_eq.name
@@ -936,7 +989,7 @@ def constrain(self, parameter, constraint_eq, params={}):
 
         # Make and store the constraint
         constraint_eq = Constraint()
-        constraint_eq.constrain(parameter, eq)
+        constraint_eq.add_constraint(parameter, eq)
         # Store the equation string so it can be shown later.
         constraint_eq.eqstr = eqstr
         self._constraints[parameter] = constraint_eq
@@ -946,6 +999,18 @@ def constrain(self, parameter, constraint_eq, params={}):
 
         return
 
+    @deprecated(constrain_deprecation_msg)
+    def constrain(self, parameter, constraint_eq, params={}):
+        """This function has been deprecated and will be removed in
+        version 4.0.0.
+
+        Please use
+        diffpy.srfit.fitbase.recipeorganizer.RecipeOrganizer.add_constraint
+        instead.
+        """
+        self.add_constraint(parameter, constraint_eq, params=params)
+        return
+
     def is_constrained(self, parameter):
         """Determine if a Parameter is constrained in this object.
 
@@ -977,14 +1042,14 @@ def isConstrained(self, parameter):
         """
         return self.is_constrained(parameter)
 
-    def unconstrain(self, *pars):
+    def remove_constraint(self, *pars):
         """Unconstrain a Parameter.
 
         This removes any constraints on a Parameter.
 
         Attributes
         ----------
-        *pars
+        *pars : str or Parameter
             The names of Parameters or Parameters to unconstrain.
 
 
@@ -1000,7 +1065,7 @@ def unconstrain(self, *pars):
                 raise ValueError("The parameter cannot be found")
 
             if parameter in self._constraints:
-                self._constraints[parameter].unconstrain()
+                self._constraints[parameter].remove_constraint()
                 del self._constraints[parameter]
                 update = True
 
@@ -1014,6 +1079,18 @@ def unconstrain(self, *pars):
 
         return
 
+    @deprecated(unconstrain_deprecation_msg)
+    def unconstrain(self, *pars):
+        """This function has been deprecated and will be removed in
+        version 4.0.0.
+
+        Please use
+        diffpy.srfit.fitbase.recipeorganizer.RecipeOrganizer.remove_constraint
+        instead.
+        """
+        self.remove_constraint(*pars)
+        return
+
     def get_constrained_parmeters(self, recurse=False):
         """Get a list of constrained managed Parameters in this object.
 
@@ -1059,7 +1136,7 @@ def clear_all_constraints(self, recurse=False):
             sub-objects are also cleared.
         """
         if self._constraints:
-            self.unconstrain(*self._constraints)
+            self.remove_constraint(*self._constraints)
 
         if recurse:
             for m in filter(_has_clear_constraints, self._iter_managed()):
@@ -1077,28 +1154,37 @@ def clearConstraints(self, recurse=False):
         """
         return self.clear_all_constraints(recurse=recurse)
 
-    def restrain(
-        self, param_or_eq, lb=-inf, ub=inf, sig=1, scaled=False, params={}
+    def add_soft_bounds(
+        self,
+        param_or_eq,
+        lower_bound=-inf,
+        upper_bound=inf,
+        sig=1,
+        scaled=False,
+        params={},
     ):
         """Restrain an expression to specified bounds.
 
+        See Notes for how the penalty is calculated.
+
         Parameters
         ----------
-        param_or_eq : str or Parameter
-            The equation string or a Parameter object to restrain.
-        lb : float, optional
+        param_or_eq : str
+            The equation or parameter to restrain.
+        lower_bound : float, optional
             The lower bound for the restraint evaluation (default is -inf).
-        ub : float, optional
+        upper_bound : float, optional
             The upper bound for the restraint evaluation (default is inf).
         sig : float, optional
             The uncertainty associated with the bounds (default is 1).
+            Please see Notes for how this is used in the penalty calculation.
         scaled : bool, optional
             If True, the restraint penalty is scaled by the unrestrained
             point-average chi^2 (chi^2/numpoints) (default is False).
         params : dict, optional
-            The dictionary of Parameters, indexed by name, that are used in the
-            equation string but are not part of the RecipeOrganizer
-            (default is {}).
+            The dictionary of Parameters, indexed by name, that are used in
+            `param_or_eq` (if an equation string is used) but are not part
+            of the RecipeOrganizer (default is {}).
 
         Returns
         -------
@@ -1111,40 +1197,81 @@ def restrain(
         The penalty is calculated as:
 
         ..
-            (max(0, lb - val, val - ub) / sig) ** 2
+            (max(0, lower_bound - val, val - upper_bound) / sig) ** 2
 
-        where `val` is the value of the evaluated equation.
+        where `val` is the value of the evaluated `param_or_eq`.
         If `scaled` is True, this penalty is multiplied by
         the average chi^2.
 
+        Examples
+        --------
+        Restraining the lattice parameters of an Ni lattice to be
+        approximately 7.4Å (2x the original lattice param)
+        can be done with the following code:
+        ..
+            recipe.add_soft_bounds(
+                "a_ni + b_ni",
+                lower_bound=7.0,
+                upper_bound=7.5,
+                sig=0.1,
+                scaled=True,
+                params={"b_ni": Parameter("b_ni", 3.473)}
+            )
+
         Raises
         ------
         ValueError
-            If `func_params` contains a name that is already used
+            If `params` contains a name that is already used
             for a Parameter.
         ValueError
             If `param_or_eq` depends on a Parameter that is not part of the
-            RecipeOrganizer and is not defined in `func_params`.
+            RecipeOrganizer and is not defined in `params`.
         """
         if isinstance(param_or_eq, str):
             eqstr = param_or_eq
-            eq = equationFromString(param_or_eq, self._eqfactory, params)
+            eq = get_equation_from_string(param_or_eq, self._eqfactory, params)
         else:
             eq = Equation(root=param_or_eq)
             eqstr = param_or_eq.name
 
         # Make and store the restraint
-        param_or_eq = Restraint(eq, lb, ub, sig, scaled)
+        param_or_eq = Restraint(eq, lower_bound, upper_bound, sig, scaled)
         param_or_eq.eqstr = eqstr
-        self.addRestraint(param_or_eq)
+        self.register_soft_bounds(param_or_eq)
         return param_or_eq
 
-    def addRestraint(self, res):
+    @deprecated(restrain_deprecation_msg)
+    def restrain(
+        self,
+        param_or_eq,
+        lower_bound=-inf,
+        upper_bound=inf,
+        sig=1,
+        scaled=False,
+        params={},
+    ):
+        """This function has been deprecated and will be removed in
+        version 4.0.0.
+
+        Please use
+        diffpy.srfit.fitbase.recipeorganizer.RecipeOrganizer.add_soft_bounds
+        instead.
+        """
+        return self.add_soft_bounds(
+            param_or_eq,
+            lower_bound=lower_bound,
+            upper_bound=upper_bound,
+            sig=sig,
+            scaled=scaled,
+            params=params,
+        )
+
+    def register_soft_bounds(self, res):
         """Add a Restraint instance to the RecipeOrganizer.
 
-        Attributes
+        Parameters
         ----------
-        res
+        res : Restraint
             A Restraint instance.
         """
         self._restraints.add(res)
@@ -1152,14 +1279,26 @@ def addRestraint(self, res):
         self._update_configuration()
         return
 
-    def unrestrain(self, *ress):
+    @deprecated(addRestraint_deprecation_msg)
+    def addRestraint(self, res):
+        """This function has been deprecated and will be removed in
+        version 4.0.0.
+
+        Please use
+        diffpy.srfit.fitbase.recipeorganizer.RecipeOrganizer.register_soft_bounds
+        instead.
+        """
+        self.register_soft_bounds(res)
+        return
+
+    def remove_soft_bounds(self, *ress):
         """Remove a Restraint from the RecipeOrganizer.
 
         Attributes
         ----------
-        *ress
-            Restraints returned from the 'restrain' method or added
-            with the 'addRestraint' method.
+        *ress : Restraint
+            The Restraints returned from the 'add_soft_bounds' method or added
+            with the 'register_soft_bounds' method.
         """
         update = False
         restuple = tuple(self._restraints)
@@ -1174,7 +1313,19 @@ def unrestrain(self, *ress):
 
         return
 
-    def clearRestraints(self, recurse=False):
+    @deprecated(unrestrain_deprecation_msg)
+    def unrestrain(self, *ress):
+        """This function has been deprecated and will be removed in
+        version 4.0.0.
+
+        Please use
+        diffpy.srfit.fitbase.recipeorganizer.RecipeOrganizer.remove_soft_bounds
+        instead.
+        """
+        self.remove_soft_bounds(*ress)
+        return
+
+    def clear_all_soft_bounds(self, recurse=False):
         """Clear all restraints.
 
         Attributes
@@ -1183,10 +1334,22 @@ def clearRestraints(self, recurse=False):
             Recurse into managed objects and clear all restraints
             found there as well.
         """
-        self.unrestrain(*self._restraints)
+        self.remove_soft_bounds(*self._restraints)
         if recurse:
             for msg in filter(_has_clear_restraints, self._iter_managed()):
-                msg.clearRestraints(recurse)
+                msg.clear_all_soft_bounds(recurse)
+        return
+
+    @deprecated(clearRestraints_deprecation_msg)
+    def clearRestraints(self, recurse=False):
+        """This function has been deprecated and will be removed in
+        version 4.0.0.
+
+        Please use
+        diffpy.srfit.fitbase.recipeorganizer.RecipeOrganizer.clear_all_soft_bounds
+        instead.
+        """
+        self.clear_all_soft_bounds(recurse=recurse)
         return
 
     def _get_constraints(self, recurse=True):
@@ -1301,12 +1464,15 @@ def _format_restraints(self):
         rset = self._get_restraints()
         rlines = []
         for res in rset:
-            line = "%s: lb = %f, ub = %f, sig = %f, scaled = %s" % (
-                res.eqstr,
-                res.lb,
-                res.ub,
-                res.sig,
-                res.scaled,
+            line = (
+                "%s: lower_bound = %f, upper_bound = %f, sig = %f, scaled = %s"
+                % (
+                    res.eqstr,
+                    res.lower_bound,
+                    res.upper_bound,
+                    res.sig,
+                    res.scaled,
+                )
             )
             rlines.append(line)
         rlines.sort(key=numstr)
@@ -1372,38 +1538,45 @@ def show(self, pattern="", textwidth=78):
 # End RecipeOrganizer
 
 
-def equationFromString(
+def get_equation_from_string(
     eqstr, factory, ns={}, buildargs=False, argclass=Parameter, argkw={}
 ):
-    """Make an equation from a string.
+    """Make an Equation object from a string.
 
     Attributes
     ----------
-    eqstr
+    eqstr : str
         A string representation of the equation. The equation must
         consist of numpy operators and "known" Parameters. Parameters
         are known if they are in ns, or already defined in the factory.
-    factory
+    factory : EquationFactory
         An EquationFactory instance.
-    ns
-        A dictionary of Parameters indexed by name that are used
+    ns : dict, optional
+        The dictionary of Parameters indexed by name that are used
         in the eqstr but not already defined in the factory
         (default {}).
-    buildargs
+    buildargs : bool, optional
         A flag indicating whether missing Parameters can be created
         by the Factory (default False). If False, then the a ValueError
         will be raised if there are undefined arguments in the eqstr.
-    argclass
+    argclass : Parameter class, optional
         Class to use when creating new Arguments (default
         Parameter). The class constructor must accept the 'name' key
         word.
-    argkw
+    argkw : dict, optional
         Key word dictionary to pass to the argclass constructor
         (default {}).
 
+    Returns
+    -------
+    eq : Equation
+        An Equation instance representing the equation in eqstr.
 
-    Raises ValueError if ns uses a name that is already defined in the factory.
-    Raises ValueError if the equation has undefined parameters.
+    Raises
+    ------
+    ValueError
+        If buildargs is False and there are undefined parameters in eqstr
+        or if ns uses a name that is already defined in the factory.
     """
 
     defined = set(factory.builders.keys())
@@ -1425,12 +1598,33 @@ def equationFromString(
     return eq
 
 
+@deprecated(equationFromString_deprecation_msg)
+def equationFromString(
+    eqstr, factory, ns={}, buildargs=False, argclass=Parameter, argkw={}
+):
+    """This function has been deprecated and will be removed in version
+    4.0.0.
+
+    Please use
+    diffpy.srfit.fitbase.recipeorganizer.get_equation_from_string
+    instead.
+    """
+    return get_equation_from_string(
+        eqstr,
+        factory,
+        ns=ns,
+        buildargs=buildargs,
+        argclass=argclass,
+        argkw=argkw,
+    )
+
+
 def _has_clear_constraints(msg):
     return hasattr(msg, "clear_all_constraints")
 
 
 def _has_clear_restraints(msg):
-    return hasattr(msg, "clearRestraints")
+    return hasattr(msg, "clear_all_soft_bounds")
 
 
 def _has_get_restraints(msg):
diff --git a/src/diffpy/srfit/fitbase/restraint.py b/src/diffpy/srfit/fitbase/restraint.py
index 5630d808..7b6ec747 100644
--- a/src/diffpy/srfit/fitbase/restraint.py
+++ b/src/diffpy/srfit/fitbase/restraint.py
@@ -33,28 +33,30 @@ class Restraint(Validatable):
 
     Attributes
     ----------
-    eq
+    eq : Equation
         An equation whose evaluation is compared against the restraint
         bounds.
-    lb
+    lower_bound : float
         The lower bound on the restraint evaluation (default -inf).
-    ub
+    upper_bound : float
         The lower bound on the restraint evaluation (default inf).
-    sig
+    sig : float
         The uncertainty on the bounds (default 1).
-    scaled
+    scaled : bool
         A flag indicating if the restraint is scaled (multiplied) by
         the unrestrained point-average chi^2 (chi^2/numpoints)
         (default False).
 
 
     The penalty is calculated as
-    (max(0, lb - val, val - ub)/sig)**2
+    (max(0, lower_bound - val, val - upper_bound)/sig)**2
     and val is the value of the calculated equation.  This is multiplied by the
     average chi^2 if scaled is True.
     """
 
-    def __init__(self, eq, lb=-inf, ub=inf, sig=1, scaled=False):
+    def __init__(
+        self, eq, lower_bound=-inf, upper_bound=inf, sig=1, scaled=False
+    ):
         """Restrain an equation to specified bounds.
 
         Attributes
@@ -62,10 +64,10 @@ def __init__(self, eq, lb=-inf, ub=inf, sig=1, scaled=False):
         eq
             An equation whose evaluation is compared against the
             restraint bounds.
-        lb
+        lower_bound
             The lower bound on the restraint evaluation (float, default
             -inf).
-        ub
+        upper_bound
             The lower bound on the restraint evaluation (float, default
             inf).
         sig
@@ -76,8 +78,8 @@ def __init__(self, eq, lb=-inf, ub=inf, sig=1, scaled=False):
             (bool, default False).
         """
         self.eq = eq
-        self.lb = float(lb)
-        self.ub = float(ub)
+        self.lower_bound = float(lower_bound)
+        self.upper_bound = float(upper_bound)
         self.sig = float(sig)
         self.scaled = bool(scaled)
         return
@@ -97,7 +99,9 @@ def penalty(self, w=1.0):
             Returns the penalty as a float
         """
         val = self.eq()
-        penalty = (max(0, self.lb - val, val - self.ub) / self.sig) ** 2
+        penalty = (
+            max(0, self.lower_bound - val, val - self.upper_bound) / self.sig
+        ) ** 2
 
         if self.scaled:
             penalty *= w
diff --git a/src/diffpy/srfit/pdf/pdfcontribution.py b/src/diffpy/srfit/pdf/pdfcontribution.py
index 9b2bdd93..07d33c5d 100644
--- a/src/diffpy/srfit/pdf/pdfcontribution.py
+++ b/src/diffpy/srfit/pdf/pdfcontribution.py
@@ -287,8 +287,8 @@ def _setup_generator(self, gen):
         gen.processMetaData()
 
         # Constrain the shared parameters
-        self.constrain(gen.qdamp, self.qdamp)
-        self.constrain(gen.qbroad, self.qbroad)
+        self.add_constraint(gen.qdamp, self.qdamp)
+        self.add_constraint(gen.qbroad, self.qbroad)
         return
 
     # Calculation setup methods
diff --git a/src/diffpy/srfit/structure/sgconstraints.py b/src/diffpy/srfit/structure/sgconstraints.py
index 8e1e198d..2b48ec1d 100644
--- a/src/diffpy/srfit/structure/sgconstraints.py
+++ b/src/diffpy/srfit/structure/sgconstraints.py
@@ -391,7 +391,7 @@ def _clear_constraints(self):
 
             for par in [scatterer.x, scatterer.y, scatterer.z]:
                 if scatterer.is_constrained(par):
-                    scatterer.unconstrain(par)
+                    scatterer.remove_constraint(par)
                 par.setConst(False)
 
         # Clear the lattice
@@ -408,7 +408,7 @@ def _clear_constraints(self):
             ]
             for par in latpars:
                 if lattice.is_constrained(par):
-                    lattice.unconstrain(par)
+                    lattice.remove_constraint(par)
                 par.setConst(False)
 
         # Clear ADPs
@@ -418,14 +418,14 @@ def _clear_constraints(self):
                     par = scatterer.get(isosymbol)
                     if par is not None:
                         if scatterer.is_constrained(par):
-                            scatterer.unconstrain(par)
+                            scatterer.remove_constraint(par)
                         par.setConst(False)
 
                 for pname in adpsymbols:
                     par = scatterer.get(pname)
                     if par is not None:
                         if scatterer.is_constrained(par):
-                            scatterer.unconstrain(par)
+                            scatterer.remove_constraint(par)
                         par.setConst(False)
 
         return
@@ -596,7 +596,7 @@ def _constrain_adps(self, positions):
                     continue
                 isoidx.append(j)
                 scatterer = scatterers[j]
-                scatterer.constrain(
+                scatterer.add_constraint(
                     isosymbol, isoname, params=self._parameters
                 )
 
@@ -691,7 +691,7 @@ def _constrain_tetragonal(lattice):
     lattice.alpha.setConst(True, ang90)
     lattice.beta.setConst(True, ang90)
     lattice.gamma.setConst(True, ang90)
-    lattice.constrain(lattice.b, lattice.a)
+    lattice.add_constraint(lattice.b, lattice.a)
     return
 
 
@@ -708,15 +708,15 @@ def _constrain_trigonal(lattice):
     ang90 = 90.0 * afactor
     ang120 = 120.0 * afactor
     if lattice.gamma.getValue() == ang120:
-        lattice.constrain(lattice.b, lattice.a)
+        lattice.add_constraint(lattice.b, lattice.a)
         lattice.alpha.setConst(True, ang90)
         lattice.beta.setConst(True, ang90)
         lattice.gamma.setConst(True, ang120)
     else:
-        lattice.constrain(lattice.b, lattice.a)
-        lattice.constrain(lattice.c, lattice.a)
-        lattice.constrain(lattice.beta, lattice.alpha)
-        lattice.constrain(lattice.gamma, lattice.alpha)
+        lattice.add_constraint(lattice.b, lattice.a)
+        lattice.add_constraint(lattice.c, lattice.a)
+        lattice.add_constraint(lattice.beta, lattice.alpha)
+        lattice.add_constraint(lattice.gamma, lattice.alpha)
     return
 
 
@@ -731,7 +731,7 @@ def _constrain_hexagonal(lattice):
         afactor = deg2rad
     ang90 = 90.0 * afactor
     ang120 = 120.0 * afactor
-    lattice.constrain(lattice.b, lattice.a)
+    lattice.add_constraint(lattice.b, lattice.a)
     lattice.alpha.setConst(True, ang90)
     lattice.beta.setConst(True, ang90)
     lattice.gamma.setConst(True, ang120)
@@ -748,8 +748,8 @@ def _constrain_cubic(lattice):
     if lattice.angunits == "rad":
         afactor = deg2rad
     ang90 = 90.0 * afactor
-    lattice.constrain(lattice.b, lattice.a)
-    lattice.constrain(lattice.c, lattice.a)
+    lattice.add_constraint(lattice.b, lattice.a)
+    lattice.add_constraint(lattice.c, lattice.a)
     lattice.alpha.setConst(True, ang90)
     lattice.beta.setConst(True, ang90)
     lattice.gamma.setConst(True, ang90)
@@ -811,7 +811,7 @@ def _makeconstraint(parname, formula, scatterer, idx, ns={}):
     # If we got here, then we have a constraint equation
     # Fix any division issues
     formula = formula.replace("/", "*1.0/")
-    scatterer.constrain(par, formula, params=ns)
+    scatterer.add_constraint(par, formula, params=ns)
     return
 
 
diff --git a/tests/conftest.py b/tests/conftest.py
index f6e58e7e..6eb797c3 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -215,11 +215,11 @@ def build_recipe_two_contributions():
     recipe.add_variable(contribution2.d, 0.1)
 
     # ---- Meaningful constraints ----
-    recipe.constrain(contribution2.m, "2*k")
-    recipe.constrain(contribution2.d, contribution1.c)
-    recipe.constrain(contribution2.B, "0.5*A")
-    recipe.restrain(contribution1.A, 0.5, 1.5)
-    recipe.restrain(contribution1.k, 0.8, 1.2)
+    recipe.add_constraint(contribution2.m, "2*k")
+    recipe.add_constraint(contribution2.d, contribution1.c)
+    recipe.add_constraint(contribution2.B, "0.5*A")
+    recipe.add_soft_bounds(contribution1.A, 0.5, 1.5)
+    recipe.add_soft_bounds(contribution1.k, 0.8, 1.2)
 
     return recipe
 
diff --git a/tests/test_constraint.py b/tests/test_constraint.py
index 00da5b14..d045feb5 100644
--- a/tests/test_constraint.py
+++ b/tests/test_constraint.py
@@ -19,11 +19,50 @@
 from diffpy.srfit.equation.builder import EquationFactory
 from diffpy.srfit.fitbase.constraint import Constraint
 from diffpy.srfit.fitbase.parameter import Parameter
-from diffpy.srfit.fitbase.recipeorganizer import equationFromString
+from diffpy.srfit.fitbase.recipeorganizer import get_equation_from_string
 
 
 class TestConstraint(unittest.TestCase):
 
+    def test_constrain_parameter(self):
+        """Test the Constraint class."""
+
+        p1 = Parameter("p1", 1)
+        p2 = Parameter("p2", 2)
+
+        factory = EquationFactory()
+
+        factory.registerArgument("p1", p1)
+        factory.registerArgument("p2", p2)
+
+        c = Constraint()
+        # Constrain p1 = 2*p2
+        eq = get_equation_from_string("2*p2", factory)
+        c.add_constraint(p1, eq)
+
+        self.assertTrue(p1.constrained)
+        self.assertFalse(p2.constrained)
+
+        eq2 = get_equation_from_string("2*p2+1", factory)
+        c2 = Constraint()
+        self.assertRaises(ValueError, c2.constrain, p1, eq2)
+        p2.setConst()
+        eq3 = get_equation_from_string("p1", factory)
+        self.assertRaises(ValueError, c2.constrain, p2, eq3)
+
+        p2.set_value(2.5)
+        c.update()
+        self.assertEqual(5.0, p1.getValue())
+
+        p2.set_value(8.1)
+        self.assertEqual(5.0, p1.getValue())
+        c.update()
+        self.assertEqual(16.2, p1.getValue())
+        return
+
+
+class TestConstraint_deprecated(unittest.TestCase):
+
     def testConstraint(self):
         """Test the Constraint class."""
 
@@ -37,17 +76,17 @@ def testConstraint(self):
 
         c = Constraint()
         # Constrain p1 = 2*p2
-        eq = equationFromString("2*p2", factory)
+        eq = get_equation_from_string("2*p2", factory)
         c.constrain(p1, eq)
 
         self.assertTrue(p1.constrained)
         self.assertFalse(p2.constrained)
 
-        eq2 = equationFromString("2*p2+1", factory)
+        eq2 = get_equation_from_string("2*p2+1", factory)
         c2 = Constraint()
         self.assertRaises(ValueError, c2.constrain, p1, eq2)
         p2.setConst()
-        eq3 = equationFromString("p1", factory)
+        eq3 = get_equation_from_string("p1", factory)
         self.assertRaises(ValueError, c2.constrain, p2, eq3)
 
         p2.set_value(2.5)
diff --git a/tests/test_fitrecipe.py b/tests/test_fitrecipe.py
index 4c928431..4da38ab8 100644
--- a/tests/test_fitrecipe.py
+++ b/tests/test_fitrecipe.py
@@ -119,7 +119,7 @@ def test_variables(self):
 
         # Constrain a parameter to the B-variable to give it a value
         p = Parameter("Bpar", -1)
-        recipe.constrain(recipe.B, p)
+        recipe.add_constraint(recipe.B, p)
         values = recipe.get_values()
         self.assertTrue((values == [2, 1, 0]).all())
         recipe.delete_variable(recipe.B)
@@ -215,11 +215,11 @@ def testResidual(self):
         # Try some constraints
         # Make c = 2*A, A = Avar
         var = self.recipe.create_new_variable("Avar")
-        self.recipe.constrain(
+        self.recipe.add_constraint(
             self.fitcontribution.c, "2*A", {"A": self.fitcontribution.A}
         )
         self.assertEqual(2, self.fitcontribution.c.value)
-        self.recipe.constrain(self.fitcontribution.A, var)
+        self.recipe.add_constraint(self.fitcontribution.A, var)
         self.assertEqual(1, var.getValue())
         self.assertEqual(self.recipe.cont.A.getValue(), var.getValue())
         # c is constrained to a constrained parameter.
@@ -232,7 +232,7 @@ def testResidual(self):
 
         # Now try some restraints. We want c to be exactly zero. It should give
         # a penalty of (c-0)**2, which is 4 in this case
-        r1 = self.recipe.restrain(self.fitcontribution.c, 0, 0, 1)
+        r1 = self.recipe.add_soft_bounds(self.fitcontribution.c, 0, 0, 1)
         self.recipe._ready = False
         res = self.recipe.residual()
         chi2 = 4 + dot(y - self.profile.y, y - self.profile.y)
@@ -240,14 +240,14 @@ def testResidual(self):
 
         # Clear the constraint and restore the value of c to 0. This should
         # give us chi2 = 0 again.
-        self.recipe.unconstrain(self.fitcontribution.c)
+        self.recipe.remove_constraint(self.fitcontribution.c)
         self.fitcontribution.c.set_value(0)
         res = self.recipe.residual([self.recipe.cont.A.getValue()])
         chi2 = 0
         self.assertAlmostEqual(chi2, dot(res, res))
 
         # Remove the restraint and variable
-        self.recipe.unrestrain(r1)
+        self.recipe.remove_soft_bounds(r1)
         self.recipe.delete_variable(self.recipe.Avar)
         self.recipe._ready = False
         res = self.recipe.residual()
@@ -255,7 +255,7 @@ def testResidual(self):
         self.assertAlmostEqual(chi2, dot(res, res))
 
         # Add constraints at the fitcontribution level.
-        self.fitcontribution.constrain(self.fitcontribution.c, "2*A")
+        self.fitcontribution.add_constraint(self.fitcontribution.c, "2*A")
         # This should evaluate to sin(x+2)
         x = self.profile.x
         y = sin(x + 2)
@@ -263,7 +263,9 @@ def testResidual(self):
         self.assertTrue(array_equal(y - self.profile.y, res))
 
         # Add a restraint at the fitcontribution level.
-        r1 = self.fitcontribution.restrain(self.fitcontribution.c, 0, 0, 1)
+        r1 = self.fitcontribution.add_soft_bounds(
+            self.fitcontribution.c, 0, 0, 1
+        )
         self.recipe._ready = False
         # The chi2 is the same as above, plus 4
         res = self.recipe.residual()
@@ -273,7 +275,7 @@ def testResidual(self):
         self.assertAlmostEqual(chi2, dot(res, res))
 
         # Remove those
-        self.fitcontribution.unrestrain(r1)
+        self.fitcontribution.remove_soft_bounds(r1)
         self.recipe._ready = False
         self.fitcontribution.unconstrain(self.fitcontribution.c)
         self.fitcontribution.c.set_value(0)
@@ -314,8 +316,8 @@ def test_boundsToRestraints():
     restraints = list(recipe._restraints)
     assert len(restraints) == 1
     r = restraints[0]
-    actual_lower_bound = r.lb
-    actual_upper_bound = r.ub
+    actual_lower_bound = r.lower_bound
+    actual_upper_bound = r.upper_bound
     actual_sigma = r.sig
     assert actual_lower_bound == expected_lower_bound
     assert actual_upper_bound == expected_upper_bound
@@ -333,8 +335,8 @@ def test_convert_bounds_to_restraints():
     restraints = list(recipe._restraints)
     assert len(restraints) == 1
     r = restraints[0]
-    assert r.lb == -1
-    assert r.ub == 1
+    assert r.lower_bound == -1
+    assert r.upper_bound == 1
     assert r.sig == 2
     assert r.scaled is True
 
@@ -361,7 +363,7 @@ def testPrintFitHook(capturestdout):
     recipe.addContribution(fitcontribution)
 
     recipe.add_variable(fitcontribution.c)
-    recipe.restrain("c", lb=5)
+    recipe.add_soft_bounds("c", lower_bound=5)
     (pfh,) = recipe.getFitHooks()
     out = capturestdout(recipe.scalar_residual)
     assert "" == out
@@ -437,7 +439,7 @@ def test_add_contribution(capturestdout):
     recipe.add_contribution(fitcontribution)
 
     recipe.add_variable(fitcontribution.c)
-    recipe.restrain("c", lb=5)
+    recipe.add_soft_bounds("c", lower_bound=5)
     (pfh,) = recipe.get_fit_hooks()
     out = capturestdout(recipe.scalar_residual)
     assert "" == out
diff --git a/tests/test_recipeorganizer.py b/tests/test_recipeorganizer.py
index bdd38058..08d46dec 100644
--- a/tests/test_recipeorganizer.py
+++ b/tests/test_recipeorganizer.py
@@ -25,6 +25,7 @@
     RecipeContainer,
     RecipeOrganizer,
     equationFromString,
+    get_equation_from_string,
 )
 
 # ----------------------------------------------------------------------------
@@ -32,8 +33,8 @@
 
 class TestEquationFromString(unittest.TestCase):
 
-    def testEquationFromString(self):
-        """Test the equationFromString method."""
+    def test_get_equation_from_string(self):
+        """Test the get_equation_from_string method."""
 
         p1 = Parameter("p1", 1)
         p2 = Parameter("p2", 2)
@@ -46,7 +47,7 @@ def testEquationFromString(self):
         factory.registerArgument("p2", p2)
 
         # Check usage where all parameters are registered with the factory
-        eq = equationFromString("p1+p2", factory)
+        eq = get_equation_from_string("p1+p2", factory)
 
         self.assertEqual(2, len(eq.args))
         self.assertTrue(p1 in eq.args)
@@ -54,7 +55,9 @@ def testEquationFromString(self):
         self.assertEqual(3, eq())
 
         # Try to use a parameter that is not registered
-        self.assertRaises(ValueError, equationFromString, "p1+p2+p3", factory)
+        self.assertRaises(
+            ValueError, get_equation_from_string, "p1+p2+p3", factory
+        )
 
         # Pass that argument in the ns dictionary
         eq = equationFromString("p1+p2+p3", factory, {"p3": p3})
@@ -69,12 +72,16 @@ def testEquationFromString(self):
 
         # Pass and use an unregistered parameter
         self.assertRaises(
-            ValueError, equationFromString, "p1+p2+p3+p4", factory, {"p3": p3}
+            ValueError,
+            get_equation_from_string,
+            "p1+p2+p3+p4",
+            factory,
+            {"p3": p3},
         )
 
         # Try to overload a registered parameter
         self.assertRaises(
-            ValueError, equationFromString, "p1+p2", factory, {"p2": p4}
+            ValueError, get_equation_from_string, "p1+p2", factory, {"p2": p4}
         )
 
         return
@@ -259,7 +266,7 @@ def testRemoveParameter(self):
         self.assertRaises(ValueError, m._remove_parameter, c)
         return
 
-    def test_constrain(self):
+    def test_constrain_parameter(self):
         """Test the constrain method."""
 
         p1 = self.m._new_parameter("p1", 1)
@@ -268,7 +275,7 @@ def test_constrain(self):
 
         self.assertFalse(p1.constrained)
         self.assertEqual(0, len(self.m._constraints))
-        self.m.constrain(p1, "2*p2")
+        self.m.add_constraint(p1, "2*p2")
 
         actual_constrained_params = self.m.getConstrainedPars()
         actual_constrained_params = [p.name for p in actual_constrained_params]
@@ -296,13 +303,13 @@ def test_constrain(self):
         self.assertRaises(ValueError, self.m.constrain, p1, "2*p2", {"p2": p3})
 
         # Remove the constraint
-        self.m.unconstrain(p1)
+        self.m.remove_constraint(p1)
         self.assertFalse(p1.constrained)
         self.assertEqual(0, len(self.m._constraints))
         self.assertFalse(self.m.isConstrained(p1))
 
         # Try an straight constraint
-        self.m.constrain(p1, p2)
+        self.m.add_constraint(p1, p2)
         p2.set_value(7)
         self.m._constraints[p1].update()
         self.assertEqual(7, p1.getValue())
@@ -314,7 +321,7 @@ def test_constrain(self):
         assert actual_constrained_params == expected_constrained_params
 
         # add constraint back and test the old function name `clearConstraints`
-        self.m.constrain(p1, p2)
+        self.m.add_constraint(p1, p2)
         actual_constrained_params = self.m.get_constrained_parmeters()
         actual_constrained_params = [p.name for p in actual_constrained_params]
         expected_constrained_params = [p1.name]
@@ -328,7 +335,7 @@ def test_constrain(self):
 
         return
 
-    def testRestrain(self):
+    def test_add_restraint(self):
         """Test the restrain method."""
 
         p1 = Parameter("p1", 1)
@@ -338,14 +345,14 @@ def testRestrain(self):
         self.m._eqfactory.registerArgument("p2", p2)
 
         self.assertEqual(0, len(self.m._restraints))
-        r = self.m.restrain("p1+p2", ub=10)
+        r = self.m.add_soft_bounds("p1+p2", upper_bound=10)
         self.assertEqual(1, len(self.m._restraints))
         p2.set_value(10)
         self.assertEqual(1, r.penalty())
-        self.m.unrestrain(r)
+        self.m.remove_soft_bounds(r)
         self.assertEqual(0, len(self.m._restraints))
 
-        r = self.m.restrain(p1, ub=10)
+        r = self.m.restrain(p1, upper_bound=10)
         self.assertEqual(1, len(self.m._restraints))
         p1.set_value(11)
         self.assertEqual(1, r.penalty())
@@ -375,8 +382,8 @@ def testGetConstraints(self):
         m2._add_parameter(p3)
         m2._add_parameter(p4)
 
-        self.m.constrain(p1, "p2")
-        m2.constrain(p3, "p4")
+        self.m.add_constraint(p1, "p2")
+        m2.add_constraint(p3, "p4")
 
         cons = self.m._get_constraints()
         self.assertTrue(p1 in cons)
@@ -402,8 +409,8 @@ def testGetRestraints(self):
         m2._add_parameter(p3)
         m2._add_parameter(p4)
 
-        r1 = self.m.restrain("p1 + p2")
-        r2 = m2.restrain("2*p3 + p4")
+        r1 = self.m.add_soft_bounds("p1 + p2")
+        r2 = m2.add_soft_bounds("2*p3 + p4")
 
         res = self.m._get_restraints()
         self.assertTrue(r1 in res)
@@ -625,14 +632,14 @@ def capture_show(*args, **kwargs):
     assert "Constraints" not in lines1
     assert "Restraints" not in lines1
     organizer._new_parameter("z", 7)
-    organizer.constrain("y", "3 * z")
+    organizer.add_constraint("y", "3 * z")
     out2 = capture_show()
     lines2 = out2.strip().split("\n")
     assert 9 == len(lines2)
     assert "Parameters" in lines2
     assert "Constraints" in lines2
     assert "Restraints" not in lines2
-    organizer.restrain("z", lb=2, ub=3, sig=0.001)
+    organizer.add_soft_bounds("z", lower_bound=2, upper_bound=3, sig=0.001)
     out3 = capture_show()
     lines3 = out3.strip().split("\n")
     assert 13 == len(lines3)
diff --git a/tests/test_restraint.py b/tests/test_restraint.py
index c5ef55c7..befae689 100644
--- a/tests/test_restraint.py
+++ b/tests/test_restraint.py
@@ -18,7 +18,7 @@
 
 from diffpy.srfit.equation.builder import EquationFactory
 from diffpy.srfit.fitbase.parameter import Parameter
-from diffpy.srfit.fitbase.recipeorganizer import equationFromString
+from diffpy.srfit.fitbase.recipeorganizer import get_equation_from_string
 from diffpy.srfit.fitbase.restraint import Restraint
 
 
@@ -36,7 +36,7 @@ def testRestraint(self):
         factory.registerArgument("p2", p2)
 
         # Restrain 1 <  p1 + p2 < 5
-        eq = equationFromString("p1 + p2", factory)
+        eq = get_equation_from_string("p1 + p2", factory)
         r = Restraint(eq, 1, 5)
 
         # This should have no penalty
@@ -63,7 +63,7 @@ def testRestraint(self):
         # Make a really large number to check the upper bound
         import numpy
 
-        r.ub = numpy.inf
+        r.upper_bound = numpy.inf
         p1.set_value(1e100)
         self.assertEqual(0, r.penalty())
 
diff --git a/tests/test_sgconstraints.py b/tests/test_sgconstraints.py
index 11859a9f..211d61be 100644
--- a/tests/test_sgconstraints.py
+++ b/tests/test_sgconstraints.py
@@ -87,13 +87,13 @@ def test_ObjCryst_constrain_space_group(pyobjcryst_available):
 
     # Make sure we can't constrain these
     with pytest.raises(ValueError):
-        mn.constrain(mn.x, "y")
+        mn.add_constraint(mn.x, "y")
 
     with pytest.raises(ValueError):
-        mn.constrain(mn.y, "z")
+        mn.add_constraint(mn.y, "z")
 
     with pytest.raises(ValueError):
-        mn.constrain(mn.z, "x")
+        mn.add_constraint(mn.z, "x")
 
     # Nor can we make them into variables
     from diffpy.srfit.fitbase.fitrecipe import FitRecipe