diff --git a/qiskit/circuit/parameterexpression.py b/qiskit/circuit/parameterexpression.py
index c12449bc9082..641a7404c553 100644
--- a/qiskit/circuit/parameterexpression.py
+++ b/qiskit/circuit/parameterexpression.py
@@ -22,7 +22,6 @@
 from qiskit.circuit.exceptions import CircuitError
 from qiskit.utils import optionals as _optionals
 
-
 # This type is redefined at the bottom to insert the full reference to "ParameterExpression", so it
 # can safely be used by runtime type-checkers like Sphinx.  Mypy does not need this because it
 # handles the references by static analysis.
@@ -522,6 +521,20 @@ def is_real(self):
                 return False
         return True
 
+    def sympify(self):
+        """Return symbolic expression as a raw Sympy or Symengine object.
+
+        Symengine is used preferentially; if both are available, the result will always be a
+        ``symengine`` object.  Symengine is a separate library but has integration with Sympy.
+
+        .. note::
+
+            This is for interoperability only.  Qiskit will not accept or work with raw Sympy or
+            Symegine expressions in its parameters, because they do not contain the tracking
+            information used in circuit-parameter binding and assignment.
+        """
+        return self._symbol_expr
+
 
 # Redefine the type so external imports get an evaluated reference; Sphinx needs this to understand
 # the type hints.
diff --git a/qiskit/transpiler/passes/optimization/template_matching/template_substitution.py b/qiskit/transpiler/passes/optimization/template_matching/template_substitution.py
index 640e6c45e2ec..2b8cee9dff7a 100644
--- a/qiskit/transpiler/passes/optimization/template_matching/template_substitution.py
+++ b/qiskit/transpiler/passes/optimization/template_matching/template_substitution.py
@@ -14,9 +14,11 @@
 Template matching substitution, given a list of maximal matches it substitutes
 them in circuit and creates a new optimized dag version of the circuit.
 """
+import collections
 import copy
+import itertools
 
-from qiskit.circuit import ParameterExpression
+from qiskit.circuit import Parameter, ParameterExpression
 from qiskit.dagcircuit.dagcircuit import DAGCircuit
 from qiskit.dagcircuit.dagdependency import DAGDependency
 from qiskit.converters.dagdependency_to_dag import dagdependency_to_dag
@@ -175,7 +177,6 @@ def _rules(self, circuit_sublist, template_sublist, template_complement):
         Returns:
             bool: True if the match respects the given rule for replacement, False otherwise.
         """
-
         if self._quantum_cost(template_sublist, template_complement):
             for elem in circuit_sublist:
                 for config in self.substitution_list:
@@ -269,11 +270,6 @@ def _remove_impossible(self):
         list_predecessors = []
         remove_list = []
 
-        # First remove any scenarios that have parameters in the template.
-        for scenario in self.substitution_list:
-            if scenario.has_parameters():
-                remove_list.append(scenario)
-
         # Initialize predecessors for each group of matches.
         for scenario in self.substitution_list:
             predecessors = set()
@@ -324,8 +320,7 @@ def _substitution(self):
 
             # Fake bind any parameters in the template
             template = self._attempt_bind(template_sublist, circuit_sublist)
-
-            if template is None:
+            if template is None or self._incr_num_parameters(template):
                 continue
 
             template_list = range(0, self.template_dag_dep.size())
@@ -432,7 +427,6 @@ def run_dag_opt(self):
                         cargs = []
                     node = group.template_dag_dep.get_node(index)
                     inst = node.op.copy()
-
                     dag_dep_opt.add_op_node(inst.inverse(), qargs, cargs)
 
             # Add the unmatched gates.
@@ -486,6 +480,11 @@ def _attempt_bind(self, template_sublist, circuit_sublist):
         solution is found then the match is valid and the parameters
         are assigned. If not, None is returned.
 
+        In order to resolve the conflict of the same parameter names in the
+        circuit and template, each variable in the template sublist is
+        re-assigned to a new dummy parameter with a completely separate name
+        if it clashes with one that exists in an input circuit.
+
         Args:
             template_sublist (list): part of the matched template.
             circuit_sublist (list): part of the matched circuit.
@@ -499,51 +498,127 @@ def _attempt_bind(self, template_sublist, circuit_sublist):
         from sympy.parsing.sympy_parser import parse_expr
 
         circuit_params, template_params = [], []
+        # Set of all parameter names that are present in the circuits to be optimised.
+        circuit_params_set = set()
 
         template_dag_dep = copy.deepcopy(self.template_dag_dep)
 
-        for idx, t_idx in enumerate(template_sublist):
+        # add parameters from circuit to circuit_params
+        for idx, _ in enumerate(template_sublist):
             qc_idx = circuit_sublist[idx]
-            circuit_params += self.circuit_dag_dep.get_node(qc_idx).op.params
-            template_params += template_dag_dep.get_node(t_idx).op.params
+            parameters = self.circuit_dag_dep.get_node(qc_idx).op.params
+            circuit_params += parameters
+            for parameter in parameters:
+                if isinstance(parameter, ParameterExpression):
+                    circuit_params_set.update(x.name for x in parameter.parameters)
+
+        _dummy_counter = itertools.count()
+
+        def dummy_parameter():
+            # Strictly not _guaranteed_ to avoid naming clashes, but if someone's calling their
+            # parameters this then that's their own fault.
+            return Parameter(f"_qiskit_template_dummy_{next(_dummy_counter)}")
+
+        # Substitutions for parameters that have clashing names between the input circuits and the
+        # defined templates.
+        template_clash_substitutions = collections.defaultdict(dummy_parameter)
+
+        # add parameters from template to template_params, replacing parameters with names that
+        # clash with those in the circuit.
+        for t_idx in template_sublist:
+            node = template_dag_dep.get_node(t_idx)
+            sub_node_params = []
+            for t_param_exp in node.op.params:
+                if isinstance(t_param_exp, ParameterExpression):
+                    for t_param in t_param_exp.parameters:
+                        if t_param.name in circuit_params_set:
+                            new_param = template_clash_substitutions[t_param.name]
+                            t_param_exp = t_param_exp.assign(t_param, new_param)
+                sub_node_params.append(t_param_exp)
+                template_params.append(t_param_exp)
+            node.op.params = sub_node_params
+
+        for node in template_dag_dep.get_nodes():
+            sub_node_params = []
+            for param_exp in node.op.params:
+                if isinstance(param_exp, ParameterExpression):
+                    for param in param_exp.parameters:
+                        if param.name in template_clash_substitutions:
+                            param_exp = param_exp.assign(
+                                param, template_clash_substitutions[param.name]
+                            )
+                sub_node_params.append(param_exp)
+
+            node.op.params = sub_node_params
 
         # Create the fake binding dict and check
-        equations, symbols, sol, fake_bind = [], set(), {}, {}
-        for t_idx, params in enumerate(template_params):
-            if isinstance(params, ParameterExpression):
-                equations.append(sym.Eq(parse_expr(str(params)), circuit_params[t_idx]))
-                for param in params.parameters:
-                    symbols.add(param)
-
-        if not symbols:
+        equations, circ_dict, temp_symbols, sol, fake_bind = [], {}, {}, {}, {}
+        for circuit_param, template_param in zip(circuit_params, template_params):
+            if isinstance(template_param, ParameterExpression):
+                if isinstance(circuit_param, ParameterExpression):
+                    circ_param_sym = circuit_param.sympify()
+                else:
+                    circ_param_sym = parse_expr(str(circuit_param))
+                equations.append(sym.Eq(template_param.sympify(), circ_param_sym))
+
+                for param in template_param.parameters:
+                    temp_symbols[param] = param.sympify()
+
+                if isinstance(circuit_param, ParameterExpression):
+                    for param in circuit_param.parameters:
+                        circ_dict[param] = param.sympify()
+            elif template_param != circuit_param:
+                # Both are numeric parameters, but aren't equal.
+                return None
+
+        if not temp_symbols:
             return template_dag_dep
 
         # Check compatibility by solving the resulting equation
-        sym_sol = sym.solve(equations)
+        sym_sol = sym.solve(equations, set(temp_symbols.values()))
         for key in sym_sol:
             try:
-                sol[str(key)] = float(sym_sol[key])
+                sol[str(key)] = ParameterExpression(circ_dict, sym_sol[key])
             except TypeError:
                 return None
 
         if not sol:
             return None
 
-        for param in symbols:
-            fake_bind[param] = sol[str(param)]
+        for key in temp_symbols:
+            fake_bind[key] = sol[str(key)]
 
         for node in template_dag_dep.get_nodes():
             bound_params = []
-
-            for param in node.op.params:
-                if isinstance(param, ParameterExpression):
-                    try:
-                        bound_params.append(float(param.bind(fake_bind)))
-                    except KeyError:
-                        return None
+            for param_exp in node.op.params:
+                if isinstance(param_exp, ParameterExpression):
+                    for param in param_exp.parameters:
+                        if param in fake_bind:
+                            if fake_bind[param] not in bound_params:
+                                param_exp = param_exp.assign(param, fake_bind[param])
                 else:
-                    bound_params.append(param)
+                    param_exp = float(param_exp)
+                bound_params.append(param_exp)
 
             node.op.params = bound_params
 
         return template_dag_dep
+
+    def _incr_num_parameters(self, template):
+        """
+        Checks if template substitution would increase the number of
+        parameters in the circuit.
+        """
+        template_params = set()
+        for param_list in (node.op.params for node in template.get_nodes()):
+            for param_exp in param_list:
+                if isinstance(param_exp, ParameterExpression):
+                    template_params.update(param_exp.parameters)
+
+        circuit_params = set()
+        for param_list in (node.op.params for node in self.circuit_dag_dep.get_nodes()):
+            for param_exp in param_list:
+                if isinstance(param_exp, ParameterExpression):
+                    circuit_params.update(param_exp.parameters)
+
+        return len(template_params) > len(circuit_params)
diff --git a/releasenotes/notes/add-parameters-to-template-substitution-a1379cdbfcc10b5c.yaml b/releasenotes/notes/add-parameters-to-template-substitution-a1379cdbfcc10b5c.yaml
new file mode 100644
index 000000000000..c93051834649
--- /dev/null
+++ b/releasenotes/notes/add-parameters-to-template-substitution-a1379cdbfcc10b5c.yaml
@@ -0,0 +1,74 @@
+---
+features:
+  - |
+    The :class:`.ParameterExpression` class is now allowed in the 
+    template optimization transpiler pass. An illustrative example 
+    of using :class:`.Parameter`\s is the following:
+
+    .. code-block:: 
+
+      from qiskit import QuantumCircuit, transpile, schedule
+      from qiskit.circuit import Parameter
+      
+      from qiskit.transpiler import PassManager
+      from qiskit.transpiler.passes import TemplateOptimization
+      
+      # New contributions to the template optimization
+      from qiskit.transpiler.passes.calibration import RZXCalibrationBuilder, rzx_templates
+      
+      from qiskit.test.mock import FakeCasablanca
+      backend = FakeCasablanca()
+      
+      phi = Parameter('φ')
+      
+      qc = QuantumCircuit(2)
+      qc.cx(0,1)
+      qc.p(2*phi, 1)
+      qc.cx(0,1)
+      print('Original circuit:')
+      print(qc)
+      
+      pass_ = TemplateOptimization(**rzx_templates.rzx_templates(['zz2']))
+      qc_cz = PassManager(pass_).run(qc)
+      print('ZX based circuit:')
+      print(qc_cz)
+      
+      # Add the calibrations
+      pass_ = RZXCalibrationBuilder(backend)
+      cal_qc = PassManager(pass_).run(qc_cz.bind_parameters({phi: 0.12}))
+      
+      # Transpile to the backend basis gates
+      cal_qct = transpile(cal_qc, backend)
+      qct = transpile(qc.bind_parameters({phi: 0.12}), backend)
+      
+      # Compare the schedule durations
+      print('Duration of schedule with the calibration:')
+      print(schedule(cal_qct, backend).duration)
+      print('Duration of standard with two CNOT gates:')
+      print(schedule(qct, backend).duration)
+          
+    outputs
+
+    .. parsed-literal::
+
+      Original circuit:
+      			 
+      q_0: ──■──────────────■──
+           ┌─┴─┐┌────────┐┌─┴─┐
+      q_1: ┤ X ├┤ P(2*φ) ├┤ X ├
+           └───┘└────────┘└───┘
+      ZX based circuit:
+      					 ┌─────────────┐            »
+      q_0: ────────────────────────────────────┤0            ├────────────»
+           ┌──────────┐┌──────────┐┌──────────┐│  Rzx(2.0*φ) │┌──────────┐»
+      q_1: ┤ Rz(-π/2) ├┤ Rx(-π/2) ├┤ Rz(-π/2) ├┤1            ├┤ Rx(-2*φ) ├»
+           └──────────┘└──────────┘└──────────┘└─────────────┘└──────────┘»
+      «                                                     
+      «q_0: ────────────────────────────────────────────────
+      «     ┌──────────┐┌──────────┐┌──────────┐┌──────────┐
+      «q_1: ┤ Rz(-π/2) ├┤ Rx(-π/2) ├┤ Rz(-π/2) ├┤ P(2.0*φ) ├
+      «     └──────────┘└──────────┘└──────────┘└──────────┘
+      Duration of schedule with the calibration:
+      1600
+      Duration of standard with two CNOT gates:
+      6848
diff --git a/test/python/transpiler/test_template_matching.py b/test/python/transpiler/test_template_matching.py
index 2c042dc57ae6..803ff50bbded 100644
--- a/test/python/transpiler/test_template_matching.py
+++ b/test/python/transpiler/test_template_matching.py
@@ -16,18 +16,37 @@
 import unittest
 import numpy as np
 from qiskit import QuantumRegister, QuantumCircuit
-from qiskit.circuit import Parameter, Gate
-from qiskit.extensions import UnitaryGate
+from qiskit.circuit import Parameter
 from qiskit.quantum_info import Operator
 from qiskit.circuit.library.templates import template_nct_2a_2, template_nct_5a_3
 from qiskit.converters.circuit_to_dag import circuit_to_dag
 from qiskit.converters.circuit_to_dagdependency import circuit_to_dagdependency
 from qiskit.transpiler import PassManager
 from qiskit.transpiler.passes import TemplateOptimization
+from qiskit.transpiler.passes.calibration.rzx_templates import rzx_templates
 from qiskit.test import QiskitTestCase
 from qiskit.transpiler.exceptions import TranspilerError
 
 
+def _ry_to_rz_template_pass(parameter: Parameter = None, extra_costs=None):
+    """Create a simple pass manager that runs a template optimisation with a single transformation.
+    It turns ``RX(pi/2).RY(parameter).RX(-pi/2)`` into the equivalent virtual ``RZ`` rotation, where
+    if ``parameter`` is given, it will be the instance used in the template."""
+    if parameter is None:
+        parameter = Parameter("_ry_rz_template_inner")
+    template = QuantumCircuit(1)
+    template.rx(-np.pi / 2, 0)
+    template.ry(parameter, 0)
+    template.rx(np.pi / 2, 0)
+    template.rz(-parameter, 0)  # pylint: disable=invalid-unary-operand-type
+
+    costs = {"rx": 16, "ry": 16, "rz": 0}
+    if extra_costs is not None:
+        costs.update(extra_costs)
+
+    return PassManager(TemplateOptimization([template], user_cost_dict=costs))
+
+
 class TestTemplateMatching(QiskitTestCase):
     """Test the TemplateOptimization pass."""
 
@@ -221,14 +240,18 @@ def test_accept_dagdependency(self):
         pass_ = TemplateOptimization(template_list=templates)
         circuit_out = PassManager(pass_).run(circuit_in)
 
-        self.assertEqual(circuit_out.count_ops().get("cx", 0), 0)
+        # these are NOT equal if template optimization works
+        self.assertNotEqual(circuit_in, circuit_out)
+
+        # however these are equivalent if the operators are the same
+        self.assertTrue(Operator(circuit_in).equiv(circuit_out))
 
     def test_parametric_template(self):
         """
         Check matching where template has parameters.
              ┌───────────┐                  ┌────────┐
         q_0: ┤ P(-1.0*β) ├──■────────────■──┤0       ├
-             ├───────────┤┌─┴─┐┌──────┐┌─┴─┐│  CZ(β) │
+             ├───────────┤┌─┴─┐┌──────┐┌─┴─┐│  CU(2β)│
         q_1: ┤ P(-1.0*β) ├┤ X ├┤ P(β) ├┤ X ├┤1       ├
              └───────────┘└───┘└──────┘└───┘└────────┘
         First test try match on
@@ -249,28 +272,14 @@ def test_parametric_template(self):
               └──────┘                          └───┘└──────┘└───┘
         """
 
-        class CZp(Gate):
-            """CZ gates used for the test."""
-
-            def __init__(self, num_qubits, params):
-                super().__init__("cz", num_qubits, params)
-
-            def inverse(self):
-                inverse = UnitaryGate(np.diag([1.0, 1.0, 1.0, np.exp(-2.0j * self.params[0])]))
-                inverse.name = "icz"
-                return inverse
-
-        def template_czp2():
-            beta = Parameter("β")
-            qc = QuantumCircuit(2)
-            qc.p(-beta, 0)
-            qc.p(-beta, 1)
-            qc.cx(0, 1)
-            qc.p(beta, 1)
-            qc.cx(0, 1)
-            qc.append(CZp(2, [beta]), [0, 1])
-
-            return qc
+        beta = Parameter("β")
+        template = QuantumCircuit(2)
+        template.p(-beta, 0)
+        template.p(-beta, 1)
+        template.cx(0, 1)
+        template.p(beta, 1)
+        template.cx(0, 1)
+        template.cu(0, 2.0 * beta, 0, 0, 0, 1)
 
         def count_cx(qc):
             """Counts the number of CX gates for testing."""
@@ -288,7 +297,10 @@ def count_cx(qc):
         circuit_in.p(3, 2)
         circuit_in.cx(1, 2)
 
-        pass_ = TemplateOptimization(template_list=[template_czp2()])
+        pass_ = TemplateOptimization(
+            template_list=[template],
+            user_cost_dict={"cx": 6, "p": 0, "cu": 8},
+        )
         circuit_out = PassManager(pass_).run(circuit_in)
 
         np.testing.assert_almost_equal(Operator(circuit_out).data[3, 3], np.exp(-4.0j))
@@ -308,54 +320,330 @@ def count_cx(qc):
         circuit_in.p(3, 2)
         circuit_in.cx(1, 2)
 
-        pass_ = TemplateOptimization(template_list=[template_czp2()])
+        pass_ = TemplateOptimization(
+            template_list=[template],
+            user_cost_dict={"cx": 6, "p": 0, "cu": 8},
+        )
         circuit_out = PassManager(pass_).run(circuit_in)
 
-        self.assertEqual(count_cx(circuit_out), 2)  # One match => two CX gates.
-        np.testing.assert_almost_equal(Operator(circuit_in).data, Operator(circuit_out).data)
+        # these are NOT equal if template optimization works
+        self.assertNotEqual(circuit_in, circuit_out)
+
+        # however these are equivalent if the operators are the same
+        self.assertTrue(Operator(circuit_in).equiv(circuit_out))
 
-    def test_unbound_parameters(self):
+    def test_optimizer_does_not_replace_unbound_partial_match(self):
         """
         Test that partial matches with parameters will not raise errors.
         This tests that if parameters are still in the temporary template after
         _attempt_bind then they will not be used.
         """
 
-        class PhaseSwap(Gate):
-            """CZ gates used for the test."""
+        beta = Parameter("β")
+        template = QuantumCircuit(2)
+        template.cx(1, 0)
+        template.cx(1, 0)
+        template.p(beta, 1)
+        template.cu(0, 0, 0, -beta, 0, 1)
 
-            def __init__(self, num_qubits, params):
-                super().__init__("p", num_qubits, params)
+        circuit_in = QuantumCircuit(2)
+        circuit_in.cx(1, 0)
+        circuit_in.cx(1, 0)
+        pass_ = TemplateOptimization(
+            template_list=[template],
+            user_cost_dict={"cx": 6, "p": 0, "cu": 8},
+        )
 
-            def inverse(self):
-                inverse = UnitaryGate(
-                    np.diag(
-                        [1.0, 1.0, np.exp(-1.0j * self.params[0]), np.exp(-1.0j * self.params[0])]
-                    )
-                )
-                inverse.name = "p"
-                return inverse
+        circuit_out = PassManager(pass_).run(circuit_in)
 
-        def template():
-            beta = Parameter("β")
-            qc = QuantumCircuit(2)
-            qc.cx(1, 0)
-            qc.cx(1, 0)
-            qc.p(beta, 1)
-            qc.append(PhaseSwap(2, [beta]), [0, 1])
+        # The template optimisation should not have replaced anything, because
+        # that would require it to leave dummy parameters in place without
+        # binding them.
+        self.assertEqual(circuit_in, circuit_out)
 
-            return qc
+    def test_unbound_parameters_in_rzx_template(self):
+        """
+        Test that rzx template ('zz2') functions correctly for a simple
+        circuit with an unbound ParameterExpression. This uses the same
+        Parameter (theta) as the template, so this also checks that template
+        substitution handle this correctly.
+        """
 
+        theta = Parameter("ϴ")
         circuit_in = QuantumCircuit(2)
+        circuit_in.cx(0, 1)
+        circuit_in.p(2 * theta, 1)
+        circuit_in.cx(0, 1)
+
+        pass_ = TemplateOptimization(**rzx_templates(["zz2"]))
+        circuit_out = PassManager(pass_).run(circuit_in)
+
+        # these are NOT equal if template optimization works
+        self.assertNotEqual(circuit_in, circuit_out)
+
+        # however these are equivalent if the operators are the same
+        theta_set = 0.42
+        self.assertTrue(
+            Operator(circuit_in.bind_parameters({theta: theta_set})).equiv(
+                circuit_out.bind_parameters({theta: theta_set})
+            )
+        )
+
+    def test_two_parameter_template(self):
+        """
+        Test a two-Parameter template based on rzx_templates(["zz3"]),
+
+                                ┌───┐┌───────┐┌───┐┌────────────┐»
+        q_0: ──■─────────────■──┤ X ├┤ Rz(φ) ├┤ X ├┤ Rz(-1.0*φ) ├»
+             ┌─┴─┐┌───────┐┌─┴─┐└─┬─┘└───────┘└─┬─┘└────────────┘»
+        q_1: ┤ X ├┤ Rz(θ) ├┤ X ├──■─────────────■────────────────»
+             └───┘└───────┘└───┘
+        «     ┌─────────┐┌─────────┐┌─────────┐┌───────────┐┌──────────────┐»
+        «q_0: ┤ Rz(π/2) ├┤ Rx(π/2) ├┤ Rz(π/2) ├┤ Rx(1.0*φ) ├┤1             ├»
+        «     └─────────┘└─────────┘└─────────┘└───────────┘│  Rzx(-1.0*φ) │»
+        «q_1: ──────────────────────────────────────────────┤0             ├»
+        «                                                   └──────────────┘»
+        «      ┌─────────┐  ┌─────────┐┌─────────┐                        »
+        «q_0: ─┤ Rz(π/2) ├──┤ Rx(π/2) ├┤ Rz(π/2) ├────────────────────────»
+        «     ┌┴─────────┴─┐├─────────┤├─────────┤┌─────────┐┌───────────┐»
+        «q_1: ┤ Rz(-1.0*θ) ├┤ Rz(π/2) ├┤ Rx(π/2) ├┤ Rz(π/2) ├┤ Rx(1.0*θ) ├»
+        «     └────────────┘└─────────┘└─────────┘└─────────┘└───────────┘»
+        «     ┌──────────────┐
+        «q_0: ┤0             ├─────────────────────────────────
+        «     │  Rzx(-1.0*θ) │┌─────────┐┌─────────┐┌─────────┐
+        «q_1: ┤1             ├┤ Rz(π/2) ├┤ Rx(π/2) ├┤ Rz(π/2) ├
+        «     └──────────────┘└─────────┘└─────────┘└─────────┘
+
+        correctly template matches into a unique circuit, but that it is
+        equivalent to the input circuit when the Parameters are bound to floats
+        and checked with Operator equivalence.
+        """
+        theta = Parameter("θ")
+        phi = Parameter("φ")
+
+        template = QuantumCircuit(2)
+        template.cx(0, 1)
+        template.rz(theta, 1)
+        template.cx(0, 1)
+        template.cx(1, 0)
+        template.rz(phi, 0)
+        template.cx(1, 0)
+        template.rz(-phi, 0)
+        template.rz(np.pi / 2, 0)
+        template.rx(np.pi / 2, 0)
+        template.rz(np.pi / 2, 0)
+        template.rx(phi, 0)
+        template.rzx(-phi, 1, 0)
+        template.rz(np.pi / 2, 0)
+        template.rz(-theta, 1)
+        template.rx(np.pi / 2, 0)
+        template.rz(np.pi / 2, 1)
+        template.rz(np.pi / 2, 0)
+        template.rx(np.pi / 2, 1)
+        template.rz(np.pi / 2, 1)
+        template.rx(theta, 1)
+        template.rzx(-theta, 0, 1)
+        template.rz(np.pi / 2, 1)
+        template.rx(np.pi / 2, 1)
+        template.rz(np.pi / 2, 1)
+
+        alpha = Parameter("$\\alpha$")
+        beta = Parameter("$\\beta$")
+
+        circuit_in = QuantumCircuit(2)
+        circuit_in.cx(0, 1)
+        circuit_in.rz(2 * alpha, 1)
+        circuit_in.cx(0, 1)
         circuit_in.cx(1, 0)
+        circuit_in.rz(3 * beta, 0)
         circuit_in.cx(1, 0)
 
-        pass_ = TemplateOptimization(template_list=[template()])
+        pass_ = TemplateOptimization(
+            [template],
+            user_cost_dict={"cx": 6, "rz": 0, "rx": 1, "rzx": 0},
+        )
         circuit_out = PassManager(pass_).run(circuit_in)
 
-        # This template will not fully match as long as gates with parameters do not
-        # commute with any other gates in the DAG dependency.
-        self.assertEqual(circuit_out.count_ops().get("cx", 0), 2)
+        # these are NOT equal if template optimization works
+        self.assertNotEqual(circuit_in, circuit_out)
+
+        # however these are equivalent if the operators are the same
+        alpha_set = 0.37
+        beta_set = 0.42
+        self.assertTrue(
+            Operator(circuit_in.bind_parameters({alpha: alpha_set, beta: beta_set})).equiv(
+                circuit_out.bind_parameters({alpha: alpha_set, beta: beta_set})
+            )
+        )
+
+    def test_exact_substitution_numeric_parameter(self):
+        """Test that a template match produces the expected value for numeric parameters."""
+        circuit_in = QuantumCircuit(1)
+        circuit_in.rx(-np.pi / 2, 0)
+        circuit_in.ry(1.45, 0)
+        circuit_in.rx(np.pi / 2, 0)
+        circuit_out = _ry_to_rz_template_pass().run(circuit_in)
+
+        expected = QuantumCircuit(1)
+        expected.rz(1.45, 0)
+        self.assertEqual(circuit_out, expected)
+
+    def test_exact_substitution_symbolic_parameter(self):
+        """Test that a template match produces the expected value for numeric parameters."""
+        a_circuit = Parameter("a")
+        circuit_in = QuantumCircuit(1)
+        circuit_in.h(0)
+        circuit_in.rx(-np.pi / 2, 0)
+        circuit_in.ry(a_circuit, 0)
+        circuit_in.rx(np.pi / 2, 0)
+        circuit_out = _ry_to_rz_template_pass(extra_costs={"h": 1}).run(circuit_in)
+
+        expected = QuantumCircuit(1)
+        expected.h(0)
+        expected.rz(a_circuit, 0)
+        self.assertEqual(circuit_out, expected)
+
+    def test_naming_clash(self):
+        """Test that the template matching works and correctly replaces a template if there is a
+        naming clash between it and the circuit.  This should include binding a partial match with a
+        parameter."""
+        # Two instances of parameters with the same name---this is how naming clashes might occur.
+        a_template = Parameter("a")
+        a_circuit = Parameter("a")
+        circuit_in = QuantumCircuit(1)
+        circuit_in.h(0)
+        circuit_in.rx(-np.pi / 2, 0)
+        circuit_in.ry(a_circuit, 0)
+        circuit_in.rx(np.pi / 2, 0)
+        circuit_out = _ry_to_rz_template_pass(a_template, extra_costs={"h": 1}).run(circuit_in)
+
+        expected = QuantumCircuit(1)
+        expected.h(0)
+        expected.rz(a_circuit, 0)
+        self.assertEqual(circuit_out, expected)
+        # Ensure that the bound parameter in the output is referentially the same as the one we put
+        # in the input circuit..
+        self.assertEqual(len(circuit_out.parameters), 1)
+        self.assertIs(circuit_in.parameters[0], a_circuit)
+        self.assertIs(circuit_out.parameters[0], a_circuit)
+
+    def test_naming_clash_in_expression(self):
+        """Test that the template matching works and correctly replaces a template if there is a
+        naming clash between it and the circuit.  This should include binding a partial match with a
+        parameter."""
+        a_template = Parameter("a")
+        a_circuit = Parameter("a")
+        circuit_in = QuantumCircuit(1)
+        circuit_in.h(0)
+        circuit_in.rx(-np.pi / 2, 0)
+        circuit_in.ry(2 * a_circuit, 0)
+        circuit_in.rx(np.pi / 2, 0)
+        circuit_out = _ry_to_rz_template_pass(a_template, extra_costs={"h": 1}).run(circuit_in)
+
+        expected = QuantumCircuit(1)
+        expected.h(0)
+        expected.rz(2 * a_circuit, 0)
+        self.assertEqual(circuit_out, expected)
+        # Ensure that the bound parameter in the output is referentially the same as the one we put
+        # in the input circuit..
+        self.assertEqual(len(circuit_out.parameters), 1)
+        self.assertIs(circuit_in.parameters[0], a_circuit)
+        self.assertIs(circuit_out.parameters[0], a_circuit)
+
+    def test_template_match_with_uninvolved_parameter(self):
+        """Test that the template matching algorithm succeeds at matching a circuit that contains an
+        unbound parameter that is not involved in the subcircuit that matches."""
+        b_circuit = Parameter("b")
+        circuit_in = QuantumCircuit(2)
+        circuit_in.rz(b_circuit, 0)
+        circuit_in.rx(-np.pi / 2, 1)
+        circuit_in.ry(1.45, 1)
+        circuit_in.rx(np.pi / 2, 1)
+        circuit_out = _ry_to_rz_template_pass().run(circuit_in)
+
+        expected = QuantumCircuit(2)
+        expected.rz(b_circuit, 0)
+        expected.rz(1.45, 1)
+        self.assertEqual(circuit_out, expected)
+
+    def test_multiple_numeric_matches_same_template(self):
+        """Test that the template matching will change both instances of a partial match within a
+        longer circuit."""
+        circuit_in = QuantumCircuit(2)
+        # Qubit 0
+        circuit_in.rx(-np.pi / 2, 0)
+        circuit_in.ry(1.32, 0)
+        circuit_in.rx(np.pi / 2, 0)
+        # Qubit 1
+        circuit_in.rx(-np.pi / 2, 1)
+        circuit_in.ry(2.54, 1)
+        circuit_in.rx(np.pi / 2, 1)
+        circuit_out = _ry_to_rz_template_pass().run(circuit_in)
+
+        expected = QuantumCircuit(2)
+        expected.rz(1.32, 0)
+        expected.rz(2.54, 1)
+        self.assertEqual(circuit_out, expected)
+
+    def test_multiple_symbolic_matches_same_template(self):
+        """Test that the template matching will change both instances of a partial match within a
+        longer circuit."""
+        a, b = Parameter("a"), Parameter("b")
+        circuit_in = QuantumCircuit(2)
+        # Qubit 0
+        circuit_in.rx(-np.pi / 2, 0)
+        circuit_in.ry(a, 0)
+        circuit_in.rx(np.pi / 2, 0)
+        # Qubit 1
+        circuit_in.rx(-np.pi / 2, 1)
+        circuit_in.ry(b, 1)
+        circuit_in.rx(np.pi / 2, 1)
+        circuit_out = _ry_to_rz_template_pass().run(circuit_in)
+
+        expected = QuantumCircuit(2)
+        expected.rz(a, 0)
+        expected.rz(b, 1)
+        self.assertEqual(circuit_out, expected)
+
+    def test_template_match_multiparameter(self):
+        """Test that the template matching works on instructions that take more than one
+        parameter."""
+        a = Parameter("a")
+        b = Parameter("b")
+        template = QuantumCircuit(1)
+        template.u(0, a, b, 0)
+        template.rz(-a - b, 0)
+
+        circuit_in = QuantumCircuit(1)
+        circuit_in.u(0, 1.23, 2.45, 0)
+        pm = PassManager(TemplateOptimization([template], user_cost_dict={"u": 16, "rz": 0}))
+        circuit_out = pm.run(circuit_in)
+
+        expected = QuantumCircuit(1)
+        expected.rz(1.23 + 2.45, 0)
+
+        self.assertEqual(circuit_out, expected)
+
+    def test_naming_clash_multiparameter(self):
+        """Test that the naming clash prevention mechanism works with instructions that take
+        multiple parameters."""
+        a_template = Parameter("a")
+        b_template = Parameter("b")
+        template = QuantumCircuit(1)
+        template.u(0, a_template, b_template, 0)
+        template.rz(-a_template - b_template, 0)
+
+        a_circuit = Parameter("a")
+        b_circuit = Parameter("b")
+        circuit_in = QuantumCircuit(1)
+        circuit_in.u(0, a_circuit, b_circuit, 0)
+        pm = PassManager(TemplateOptimization([template], user_cost_dict={"u": 16, "rz": 0}))
+        circuit_out = pm.run(circuit_in)
+
+        expected = QuantumCircuit(1)
+        expected.rz(a_circuit + b_circuit, 0)
+
+        self.assertEqual(circuit_out, expected)
 
 
 if __name__ == "__main__":