Repository: PennyLaneAI/pennylane
Pull request #1610 (2021-08-30 21:29:26)

Create qchem hf feature branch
**Context:**
This PR creates a feature branch for adding the differentiable Hartree-Fock solver functions. The approved PRs related to the solver will be merged into this feature branch, instead of the pennylane master branch, until we have the minimum requirements to perform a full differentiable Hartree-Fock calculation. Then this feature branch will be merged into master.

**Description of the Change:**
`hf` sub-folder is created inside `pennylane` to store all the new differentiable Hartree-Fock  additions. 
`hf` sub-folder is created inside `tests` to store all the differentiable Hartree-Fock unit tests. 

**Benefits:**
This feature branch collects all the approved HF additions and prevents adding an incomplete HF feature to PL.

**Possible Drawbacks:**
No specific drawbacks.

**Related GitHub Issues:**
No related issues.


github-actions[bot]: Hello. You may have forgotten to update the changelog!
Please edit [.github/CHANGELOG.md](/PennyLaneAI/pennylane/blob/master/.github/CHANGELOG.md) with:
* A one-to-two sentence description of the change. You may include a small working example for new features.
* A link back to this PR.
* Your name (or GitHub username) in the contributors section.
codecov[bot]: # [Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1610?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
> Merging [#1610](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1610?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (a895650) into [master](https://codecov.io/gh/PennyLaneAI/pennylane/commit/f147f22d8d99ba5891edd45a6a1f7dd679c8a23c?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (f147f22) will **not change** coverage.
> The diff coverage is `n/a`.

[![Impacted file tree graph](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1610/graphs/tree.svg?width=650&height=150&src=pr&token=9S3YzlUgGk&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI)](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1610?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI)

```diff
@@           Coverage Diff           @@
##           master    #1610   +/-   ##
=======================================
  Coverage   99.13%   99.13%           
=======================================
  Files         195      195           
  Lines       14114    14114           
=======================================
  Hits        13992    13992           
  Misses        122      122           
```



------

[Continue to review full report at Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1610?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1610?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI). Last update [f147f22...a895650](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1610?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).


Patch:
diff --git a/doc/code/qml_hf.rst b/doc/code/qml_hf.rst
new file mode 100644
index 000000000..33db7274b
--- /dev/null
+++ b/doc/code/qml_hf.rst
@@ -0,0 +1,7 @@
+qml.hf
+======
+
+.. currentmodule:: pennylane.hf
+
+.. automodapi:: pennylane.hf
+
diff --git a/doc/index.rst b/doc/index.rst
index 0989077e0..d9c579cc3 100644
--- a/doc/index.rst
+++ b/doc/index.rst
@@ -211,6 +211,7 @@ PennyLane is **free** and **open source**, released under the Apache License, Ve
    code/qml_kernels
    code/qml_math
    code/qml_qaoa
+   code/qml_hf
    code/qml_qchem
    code/qml_qnn
    code/qml_tape
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 2e1a0be60..88daccd80 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -4,6 +4,101 @@
 
 <h3>New features since last release</h3>
 
+* A differentiable Hartree-Fock (HF) solver has been added. It can be used to construct molecular Hamiltonians 
+  that can be differentiated with respect to nuclear coordinates and basis-set parameters.
+  [(#1610)](https://github.com/PennyLaneAI/pennylane/pull/1610)
+
+  The HF solver computes the integrals over basis functions, constructs the relevant matrices, and
+  performs self-consistent-field iterations to obtain a set of optimized molecular orbital
+  coefficients. These coefficients and the computed integrals over basis functions are used to
+  construct the one- and two-body electron integrals in the molecular orbital basis which can be
+  used to generate a differentiable second-quantized Hamiltonian in the fermionic and qubit basis.
+
+  The following code shows the construction of the Hamiltonian for the hydrogen molecule where the
+  geometry of the molecule and the basis set parameters are all differentiable.
+
+  ```python
+  import pennylane as qml
+  from pennylane import numpy as np
+  
+  symbols = ["H", "H"]
+  geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 2.0]], requires_grad=True)
+  alpha = np.array([[3.42525091, 0.62391373, 0.1688554],
+                    [3.42525091, 0.62391373, 0.1688554]], requires_grad = True)
+  coeff = np.array([[0.15432897, 0.53532814, 0.44463454],
+                    [0.15432897, 0.53532814, 0.44463454]], requires_grad = True)
+
+  # we create a molecule object with differentiable atomic coordinates and basis set parameters
+  # alpha and coeff are the exponentents and contraction coefficients of the Gaussian functions
+  mol = qml.hf.Molecule(symbols, geometry, alpha=alpha, coeff=coeff)
+  args = [geometry, alpha, coeff] # initial values of the differentiable parameters
+  
+  hamiltonian = qml.hf.generate_hamiltonian(mol)(*args)
+  ```
+
+  The generated Hamiltonian can be used in a circuit where the molecular geometry, the basis set
+  parameters and the circuit parameters are optimized simultaneously.
+
+  ```python
+  import autograd
+  
+  params = [np.array([0.0], requires_grad=True)]
+  dev = qml.device("default.qubit", wires=4)
+  hf_state = np.array([1, 1, 0, 0])
+  
+  def generate_circuit(mol):
+      @qml.qnode(dev)
+      def circuit(*args):
+          qml.BasisState(hf_state, wires=[0, 1, 2, 3])
+          qml.DoubleExcitation(*args[0][0], wires=[0, 1, 2, 3])
+          return qml.expval(hf.generate_hamiltonian(mol)(*args[1:]))
+      return circuit
+  
+  for n in range(10): # geometry and parameter optimization loop
+  
+      # we create a molecule object with differentiable atomic coordinates and basis set parameters
+      # alpha and coeff are the exponentents and contraction coefficients of the Gaussian functions 
+      mol = hf.Molecule(symbols, geometry, alpha=alpha, coeff=coeff)
+      args_ = [params, *args] # initial values of the differentiable parameters
+  
+      # compute gradients with respect to the circuit parameters and update the parameters
+      g_params = autograd.grad(generate_circuit(mol), argnum = 0)(*args_)
+      params = params - 0.1 * g_params[0]
+  
+      # compute gradients with respect to the nuclear coordinates and update geometry
+      forces = autograd.grad(generate_circuit(mol), argnum = 1)(*args_)
+      geometry = geometry - 0.5 * forces
+  
+      # compute gradients with respect to the Gaussian exponents and update the exponents
+      g_alpha = autograd.grad(generate_circuit(mol), argnum = 2)(*args_)
+      alpha = alpha - 0.1 * g_alpha
+  
+      # compute gradients with respect to the Gaussian contraction coefficients and update them
+      g_coeff = autograd.grad(generate_circuit(mol), argnum = 3)(*args_)
+      coeff = coeff - 0.1 * g_coeff
+  ```
+
+  The components of the HF solver can also be differentiated individually. For instance, the overlap
+  integral can be differentiated with respect to the basis set parameters as
+
+  ```python
+  symbols = ["H", "H"]
+  geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 2.0]], requires_grad=False)
+  alpha = np.array([[3.42525091, 0.62391373, 0.1688554],
+                    [3.42525091, 0.62391373, 0.1688554]], requires_grad = True)
+  coeff = np.array([[0.15432897, 0.53532814, 0.44463454],
+                    [0.15432897, 0.53532814, 0.44463454]], requires_grad = True)
+
+  mol = qml.hf.Molecule(symbols, geometry, alpha=alpha, coeff=coeff)
+  args = [alpha, coeff]
+  
+  a = mol.basis_set[0]
+  b = mol.basis_set[1]
+
+  g_alpha = autograd.grad(qml.hf.generate_overlap(a, b), argnum = 0)(*args)
+  g_coeff = autograd.grad(qml.hf.generate_overlap(a, b), argnum = 1)(*args)
+  ```
+
 * The `insert` transform has now been added, providing a way to insert single-qubit operations into
   a quantum circuit. The transform can apply to quantum functions, tapes, and devices.
   [(#1795)](https://github.com/PennyLaneAI/pennylane/pull/1795)
@@ -932,6 +1027,8 @@
 
 This release contains contributions from (in alphabetical order):
 
-Utkarsh Azad, Akash Narayanan B, Sam Banning, Thomas Bromley, Olivia Di Matteo, Andrew Gardhouse, David Ittah, 
-Josh Izaac, Christina Lee, Romain Moyard, Carrie-Anne Rubidge, Maria Schuld, Rishabh Singh, Jay Soni, Ingrid Strandberg, 
-Antal Száva, Cody Wang, David Wierichs, Moritz Willmann.
+Juan Miguel Arrazola, Utkarsh Azad, Akash Narayanan B, Sam Banning, Thomas Bromley, Jack Ceroni,
+Alain Delgado, Olivia Di Matteo, Andrew Gardhouse, David Ittah, Josh Izaac, Soran Jahangiri,
+Christina Lee, Romain Moyard, Carrie-Anne Rubidge, Maria Schuld, Rishabh Singh, Jay Soni,
+Ingrid Strandberg, Antal Száva, Teresa Tamayo-Mendoza, Rodrigo Vargas, Cody Wang, David Wierichs,
+Moritz Willmann.
diff --git a/pennylane/__init__.py b/pennylane/__init__.py
index 309faceb6..3b0bd6d11 100644
--- a/pennylane/__init__.py
+++ b/pennylane/__init__.py
@@ -31,6 +31,7 @@ import pennylane.math
 import pennylane.operation
 import pennylane.qnn
 import pennylane.templates
+import pennylane.hf
 from pennylane._device import Device, DeviceError
 from pennylane._grad import grad, jacobian, finite_diff
 from pennylane._qubit_device import QubitDevice
diff --git a/pennylane/hf/__init__.py b/pennylane/hf/__init__.py
new file mode 100644
index 000000000..bb0353074
--- /dev/null
+++ b/pennylane/hf/__init__.py
@@ -0,0 +1,46 @@
+# Copyright 2018-2021 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""
+This subpackage provides the functionality to perform differentiable Hartree-Fock calculations.
+"""
+from .basis_data import STO3G, atomic_numbers
+from .basis_set import BasisFunction, atom_basis_data, mol_basis_data
+from .hamiltonian import (
+    generate_electron_integrals,
+    generate_fermionic_hamiltonian,
+    generate_hamiltonian,
+)
+from .hartree_fock import generate_scf, hf_energy, nuclear_energy
+from .integrals import (
+    contracted_norm,
+    electron_repulsion,
+    expansion,
+    gaussian_kinetic,
+    gaussian_overlap,
+    generate_attraction,
+    generate_kinetic,
+    generate_overlap,
+    generate_repulsion,
+    nuclear_attraction,
+    primitive_norm,
+)
+from .matrices import (
+    generate_attraction_matrix,
+    generate_core_matrix,
+    generate_kinetic_matrix,
+    generate_overlap_matrix,
+    generate_repulsion_tensor,
+    molecular_density_matrix,
+)
+from .molecule import Molecule
diff --git a/pennylane/hf/basis_data.py b/pennylane/hf/basis_data.py
new file mode 100644
index 000000000..5ef70f2b5
--- /dev/null
+++ b/pennylane/hf/basis_data.py
@@ -0,0 +1,140 @@
+# Copyright 2018-2021 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""
+This module contains basis set parameters defining Gaussian-type orbitals for a selected number of
+atoms. The data are taken from the Basis Set Exchange `library <https://www.basissetexchange.org>`_.
+The current data include the STO-3G basis set for atoms with atomic numbers 1-10.
+"""
+
+atomic_numbers = {
+    "H": 1,
+    "He": 2,
+    "Li": 3,
+    "Be": 4,
+    "B": 5,
+    "C": 6,
+    "N": 7,
+    "O": 8,
+    "F": 9,
+    "Ne": 10,
+}
+
+STO3G = {
+    "H": {
+        "orbitals": ["S"],
+        "exponents": [[0.3425250914e01, 0.6239137298e00, 0.1688554040e00]],
+        "coefficients": [[0.1543289673e00, 0.5353281423e00, 0.4446345422e00]],
+    },
+    "He": {
+        "orbitals": ["S"],
+        "exponents": [[0.6362421394e01, 0.1158922999e01, 0.3136497915e00]],
+        "coefficients": [[0.1543289673e00, 0.5353281423e00, 0.4446345422e00]],
+    },
+    "Li": {
+        "orbitals": ["S", "SP"],
+        "exponents": [
+            [0.1611957475e02, 0.2936200663e01, 0.7946504870e00],
+            [0.6362897469e00, 0.1478600533e00, 0.4808867840e-01],
+        ],
+        "coefficients": [
+            [0.1543289673e00, 0.5353281423e00, 0.4446345422e00],
+            [-0.9996722919e-01, 0.3995128261e00, 0.7001154689e00],
+            [0.1559162750e00, 0.6076837186e00, 0.3919573931e00],
+        ],
+    },
+    "Be": {
+        "orbitals": ["S", "SP"],
+        "exponents": [
+            [0.3016787069e02, 0.5495115306e01, 0.1487192653e01],
+            [0.1314833110e01, 0.3055389383e00, 0.9937074560e-01],
+        ],
+        "coefficients": [
+            [0.1543289673e00, 0.5353281423e00, 0.4446345422e00],
+            [-0.9996722919e-01, 0.3995128261e00, 0.7001154689e00],
+            [0.1559162750e00, 0.6076837186e00, 0.3919573931e00],
+        ],
+    },
+    "B": {
+        "orbitals": ["S", "SP"],
+        "exponents": [
+            [0.4879111318e02, 0.8887362172e01, 0.2405267040e01],
+            [0.2236956142e01, 0.5198204999e00, 0.1690617600e00],
+        ],
+        "coefficients": [
+            [0.1543289673e00, 0.5353281423e00, 0.4446345422e00],
+            [-0.9996722919e-01, 0.3995128261e00, 0.7001154689e00],
+            [0.1559162750e00, 0.6076837186e00, 0.3919573931e00],
+        ],
+    },
+    "C": {
+        "orbitals": ["S", "SP"],
+        "exponents": [
+            [0.7161683735e02, 0.1304509632e02, 0.3530512160e01],
+            [0.2941249355e01, 0.6834830964e00, 0.2222899159e00],
+        ],
+        "coefficients": [
+            [0.1543289673e00, 0.5353281423e00, 0.4446345422e00],
+            [-0.9996722919e-01, 0.3995128261e00, 0.7001154689e00],
+            [0.1559162750e00, 0.6076837186e00, 0.3919573931e00],
+        ],
+    },
+    "N": {
+        "orbitals": ["S", "SP"],
+        "exponents": [
+            [0.9910616896e02, 0.1805231239e02, 0.4885660238e01],
+            [0.3780455879e01, 0.8784966449e00, 0.2857143744e00],
+        ],
+        "coefficients": [
+            [0.1543289673e00, 0.5353281423e00, 0.4446345422e00],
+            [-0.9996722919e-01, 0.3995128261e00, 0.7001154689e00],
+            [0.1559162750e00, 0.6076837186e00, 0.3919573931e00],
+        ],
+    },
+    "O": {
+        "orbitals": ["S", "SP"],
+        "exponents": [
+            [0.1307093214e03, 0.2380886605e02, 0.6443608313e01],
+            [0.5033151319e01, 0.1169596125e01, 0.3803889600e00],
+        ],
+        "coefficients": [
+            [0.1543289673e00, 0.5353281423e00, 0.4446345422e00],
+            [-0.9996722919e-01, 0.3995128261e00, 0.7001154689e00],
+            [0.1559162750e00, 0.6076837186e00, 0.3919573931e00],
+        ],
+    },
+    "F": {
+        "orbitals": ["S", "SP"],
+        "exponents": [
+            [0.1666791340e03, 0.3036081233e02, 0.8216820672e01],
+            [0.6464803249e01, 0.1502281245e01, 0.4885884864e00],
+        ],
+        "coefficients": [
+            [0.1543289673e00, 0.5353281423e00, 0.4446345422e00],
+            [-0.9996722919e-01, 0.3995128261e00, 0.7001154689e00],
+            [0.1559162750e00, 0.6076837186e00, 0.3919573931e00],
+        ],
+    },
+    "Ne": {
+        "orbitals": ["S", "SP"],
+        "exponents": [
+            [0.2070156070e03, 0.3770815124e02, 0.1020529731e02],
+            [0.8246315120e01, 0.1916266291e01, 0.6232292721e00],
+        ],
+        "coefficients": [
+            [0.1543289673e00, 0.5353281423e00, 0.4446345422e00],
+            [-0.9996722919e-01, 0.3995128261e00, 0.7001154689e00],
+            [0.1559162750e00, 0.6076837186e00, 0.3919573931e00],
+        ],
+    },
+}
diff --git a/pennylane/hf/basis_set.py b/pennylane/hf/basis_set.py
new file mode 100644
index 000000000..74ad62464
--- /dev/null
+++ b/pennylane/hf/basis_set.py
@@ -0,0 +1,139 @@
+# Copyright 2018-2021 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""
+This module contains functions and classes to create a BasisFunction object from standard basis sets
+such as STO-3G.
+"""
+# pylint: disable=too-few-public-methods
+from pennylane.hf.basis_data import STO3G
+
+
+class BasisFunction:
+    r"""Create a basis function object.
+
+    A basis set is composed of a set of basis functions that are typically constructed as a linear
+    combination of primitive Gaussian functions. For instance, a basis function in the STO-3G basis
+    set is formed as
+
+    .. math::
+
+        \psi = a_1 G_1 + a_2 G_2 + a_3 G_3,
+
+    where :math:`a` denotes the contraction coefficients and :math:`G` is a Gaussian function
+    defined as
+
+    .. math::
+
+        G = x^l y^m z^n e^{-\alpha r^2}.
+
+    Each Gaussian function is characterized by the angular momentum numbers :math:`(l, m, n)` that
+    determine the type of the orbital, the exponent :math:`\alpha` and the position vector
+    :math:`r = (x, y, z)`. These parameters and the contraction coefficients :math:`a` define
+    atomic basis functions. Predefined values of the exponents and contraction coefficients for
+    each atomic orbital of a given chemical element can be obtained from reference libraries such as
+    the Basis Set Exchange `library <https://www.basissetexchange.org>`_.
+
+    The basis function object created by the BasisFunction class stores all the basis set parameters
+    including the angular momentum, exponents, positions and coefficients of the Gaussian functions.
+
+    The basis function object can be easily passed to the functions that compute various types of
+    integrals over such functions, e.g., overlap integrals, which are essential for Hartree-Fock
+    calculations.
+
+    Args:
+        l (tuple[int]): angular momentum numbers of the basis function.
+        alpha (array(float)): exponents of the primitive Gaussian functions
+        coeff (array(float)): coefficients of the contracted Gaussian functions
+        r (array(float)): positions of the Gaussian functions
+    """
+
+    def __init__(self, l, alpha, coeff, r):
+        self.l = l
+        self.alpha = alpha
+        self.coeff = coeff
+        self.r = r
+        self.params = [self.alpha, self.coeff, self.r]
+
+
+def atom_basis_data(name, atom):
+    r"""Generate default basis set parameters for an atom.
+
+    This function extracts the angular momentum, exponents, and contraction coefficients of
+    Gaussian functions forming atomic orbitals for a given atom. These values are taken from the
+    basis set data provided in ``basis_data.py``.
+
+    Args:
+        name (str): name of the basis set
+        atom (str): atomic symbol of the chemical element
+
+    Returns:
+        list(tuple): tuple containing the angular momentum, the exponents and contraction
+        coefficients of a basis function
+
+    **Example**
+
+    >>> params = atom_basis_data('sto-3g', 'H')
+    >>> print(params)
+    [((0, 0, 0), [3.425250914, 0.6239137298, 0.168855404], [0.1543289673, 0.5353281423, 0.4446345422])]
+    """
+    basis_sets = {"sto-3g": STO3G}
+
+    s = [(0, 0, 0)]
+    p = [(1, 0, 0), (0, 1, 0), (0, 0, 1)]  # for px, py, pz, respectively
+
+    basis = basis_sets[name][atom]
+    params = []
+    for i, j in enumerate(basis["orbitals"]):
+        if j == "S":
+            params.append((s[0], basis["exponents"][i], basis["coefficients"][i]))
+        elif j == "SP":
+            for term in j:
+                if term == "S":
+                    params.append((s[0], basis["exponents"][i], basis["coefficients"][i]))
+                if term == "P":
+                    for l in p:
+                        params.append((l, basis["exponents"][i], basis["coefficients"][i + 1]))
+    return params
+
+
+def mol_basis_data(name, symbols):
+    r"""Generates default basis set parameters for a molecule.
+
+    This function generates the default basis set parameters for a list of atomic symbols and
+    computes the total number of basis functions for each atom.
+
+    Args:
+        name (str): name of the basis set
+        symbols (list[str]): symbols of the atomic species in the molecule
+
+    Returns:
+        tuple(list, tuple): the number of atomic basis functions and the basis set parameters for
+        each atom in the molecule
+
+    **Example**
+
+    >>> n_basis, params = mol_basis_data('sto-3g', ['H', 'H'])
+    >>> print(n_basis)
+    [1, 1]
+    >>> print(params)
+    (((0, 0, 0), [3.425250914, 0.6239137298, 0.168855404], [0.1543289673, 0.5353281423, 0.4446345422]),
+     ((0, 0, 0), [3.425250914, 0.6239137298, 0.168855404], [0.1543289673, 0.5353281423, 0.4446345422]))
+    """
+    n_basis = []
+    basis_set = []
+    for s in symbols:
+        basis = atom_basis_data(name, s)
+        n_basis += [len(basis)]
+        basis_set += basis
+    return n_basis, tuple(basis_set)
diff --git a/pennylane/hf/hamiltonian.py b/pennylane/hf/hamiltonian.py
new file mode 100644
index 000000000..a00d6016d
--- /dev/null
+++ b/pennylane/hf/hamiltonian.py
@@ -0,0 +1,459 @@
+# Copyright 2018-2021 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""
+This module contains the functions needed for computing the molecular Hamiltonian.
+"""
+# pylint: disable= too-many-branches, too-many-arguments, too-many-locals, too-many-nested-blocks
+import autograd.numpy as anp
+import pennylane as qml
+from pennylane import numpy as np
+from pennylane.hf.hartree_fock import generate_scf, nuclear_energy
+
+
+def generate_electron_integrals(mol, core=None, active=None):
+    r"""Return a function that computes the one- and two-electron integrals in the molecular orbital
+    basis.
+
+    The one- and two-electron integrals are required to construct a molecular Hamiltonian in the
+    second-quantized form
+
+    .. math::
+
+        H = \sum_{pq} h_{pq} c_p^{\dagger} c_q + \frac{1}{2} \sum_{pqrs} h_{pqrs} c_p^{\dagger} c_q^{\dagger} c_r c_s,
+
+    where :math:`c^{\dagger}` and :math:`c` are the creation and annihilation operators,
+    respectively, and :math:`h_{pq}` and :math:`h_{pqrs}` are the one- and two-electron integrals.
+    These integrals can be computed by integrating over molecular orbitals :math:`\phi` as
+
+    .. math::
+
+        h_{pq} = \int \phi_p(r)^* \left ( -\frac{\nabla_r^2}{2} - \sum_i \frac{Z_i}{|r-R_i|} \right )  \phi_q(r) dr,
+
+    and
+
+    .. math::
+
+        h_{pqrs} = \int \frac{\phi_p(r_1)^* \phi_q(r_2)^* \phi_r(r_2) \phi_s(r_1)}{|r_1 - r_2|} dr_1 dr_2.
+
+    The molecular orbitals are constructed as a linear combination of atomic orbitals as
+
+    .. math::
+
+        \phi_i = \sum_{\nu}c_{\nu}^i \chi_{\nu}.
+
+    The one- and two-electron integrals can be written in the molecular orbital basis as
+
+    .. math::
+
+        h_{pq} = \sum_{\mu \nu} C_{p \mu} h_{\mu \nu} C_{\nu q},
+
+    and
+
+    .. math::
+
+        h_{pqrs} = \sum_{\mu \nu \rho \sigma} C_{p \mu} C_{q \nu} h_{\mu \nu \rho \sigma} C_{\rho r} C_{\sigma s}.
+
+    The :math:`h_{\mu \nu}` and :math:`h_{\mu \nu \rho \sigma}` terms refer to the elements of the
+    core matrix and the electron repulsion tensor, respectively, and :math:`C` is the molecular
+    orbital expansion coefficient matrix.
+
+    Args:
+        mol (Molecule): the molecule object
+        core (list[int]): indices of the core orbitals
+        active (list[int]): indices of the active orbitals
+
+    Returns:
+        function: function that computes the core constant, and the one- and two-electron integrals
+
+    **Example**
+
+    >>> symbols  = ['H', 'H']
+    >>> geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad = False)
+    >>> alpha = np.array([[3.42525091, 0.62391373, 0.1688554],
+    >>>                   [3.42525091, 0.62391373, 0.1688554]], requires_grad=True)
+    >>> mol = Molecule(symbols, geometry, alpha=alpha)
+    >>> args = [alpha]
+    >>> generate_electron_integrals(mol)(*args)
+    (1.0,
+     array([[-1.3902192695e+00,  0.0000000000e+00],
+            [-4.4408920985e-16, -2.9165331336e-01]]),
+     array([[[[ 7.1443907755e-01, -2.7755575616e-17],
+              [ 5.5511151231e-17,  1.7024144301e-01]],
+             [[ 5.5511151231e-17,  1.7024144301e-01],
+              [ 7.0185315353e-01,  6.6613381478e-16]]],
+            [[[-1.3877787808e-16,  7.0185315353e-01],
+              [ 1.7024144301e-01,  2.2204460493e-16]],
+             [[ 1.7024144301e-01, -4.4408920985e-16],
+              [ 6.6613381478e-16,  7.3883668974e-01]]]]))
+    """
+
+    def electron_integrals(*args):
+        r"""Compute the one- and two-electron integrals in the molecular orbital basis.
+
+        Args:
+            args (array[array[float]]): initial values of the differentiable parameters
+
+        Returns:
+            tuple[array[float]]: 1D tuple containing core constant, one- and two-electron integrals
+        """
+        _, coeffs, _, h_core, repulsion_tensor = generate_scf(mol)(*args)
+        one = anp.einsum("qr,rs,st->qt", coeffs.T, h_core, coeffs)
+        two = anp.swapaxes(
+            anp.einsum(
+                "ab,cd,bdeg,ef,gh->acfh", coeffs.T, coeffs.T, repulsion_tensor, coeffs, coeffs
+            ),
+            1,
+            3,
+        )
+        core_constant = nuclear_energy(mol.nuclear_charges, mol.coordinates)(*args)
+
+        if core is None and active is None:
+            return core_constant, one, two
+
+        for i in core:
+            core_constant = core_constant + 2 * one[i][i]
+            for j in core:
+                core_constant = core_constant + 2 * two[i][j][j][i] - two[i][j][i][j]
+
+        for p in active:
+            for q in active:
+                for i in core:
+                    o = anp.zeros(one.shape)
+                    o[p, q] = 1.0
+                    one = one + (2 * two[i][p][q][i] - two[i][p][i][q]) * o
+
+        one = one[anp.ix_(active, active)]
+        two = two[anp.ix_(active, active, active, active)]
+
+        return core_constant, one, two
+
+    return electron_integrals
+
+
+def generate_fermionic_hamiltonian(mol, cutoff=1.0e-12, core=None, active=None):
+    r"""Return a function that computes the fermionic hamiltonian.
+
+    Args:
+        mol (Molecule): the molecule object
+        cutoff (float): cutoff value for discarding the negligible electronic integrals
+
+    Returns:
+        function: function that computes the fermionic hamiltonian
+
+    **Example**
+
+    >>> symbols  = ['H', 'H']
+    >>> geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad = False)
+    >>> alpha = np.array([[3.42525091, 0.62391373, 0.1688554],
+    >>>                   [3.42525091, 0.62391373, 0.1688554]], requires_grad=True)
+    >>> mol = Molecule(symbols, geometry, alpha=alpha)
+    >>> args = [alpha]
+    >>> h = generate_fermionic_hamiltonian(mol)(*args)
+    """
+
+    def fermionic_hamiltonian(*args):
+        r"""Compute the fermionic hamiltonian.
+
+        Args:
+            args (array[array[float]]): initial values of the differentiable parameters
+
+        Returns:
+            tuple(array[float], list[list[int]]): the Hamiltonian coefficients and operators
+        """
+        core_constant, one, two = generate_electron_integrals(mol, core, active)(*args)
+
+        core_constant = anp.array([core_constant])
+
+        indices_one = anp.argwhere(abs(one) >= cutoff)
+        operators_one = (indices_one * 2).tolist() + (
+            indices_one * 2 + 1
+        ).tolist()  # up-up + down-down terms
+        coeffs_one = anp.tile(one[abs(one) >= cutoff], 2)
+
+        indices_two = anp.argwhere(abs(two) >= cutoff)
+        n = len(indices_two)
+        operators_two = (
+            [(indices_two[i] * 2).tolist() for i in range(n)]  # up-up-up-up term
+            + [
+                (indices_two[i] * 2 + [0, 1, 1, 0]).tolist() for i in range(n)
+            ]  # up-down-down-up term
+            + [
+                (indices_two[i] * 2 + [1, 0, 0, 1]).tolist() for i in range(n)
+            ]  # down-up-up-down term
+            + [(indices_two[i] * 2 + 1).tolist() for i in range(n)]  # down-down-down-down term
+        )
+        coeffs_two = anp.tile(two[abs(two) >= cutoff], 4) / 2
+
+        coeffs = anp.concatenate((core_constant, coeffs_one, coeffs_two))
+        operators = [[]] + operators_one + operators_two
+        indices_sort = [operators.index(i) for i in sorted(operators)]
+
+        return coeffs[indices_sort], sorted(operators)
+
+    return fermionic_hamiltonian
+
+
+def generate_hamiltonian(mol, cutoff=1.0e-12, core=None, active=None):
+    r"""Return a function that computes the qubit hamiltonian.
+
+    Args:
+        mol (Molecule): the molecule object
+        cutoff (float): cutoff value for discarding the negligible electronic integrals
+
+    Returns:
+        function: function that computes the qubit hamiltonian
+
+    **Example**
+
+    >>> symbols  = ['H', 'H']
+    >>> geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad = False)
+    >>> alpha = np.array([[3.42525091, 0.62391373, 0.1688554],
+    >>>                   [3.42525091, 0.62391373, 0.1688554]], requires_grad=True)
+    >>> mol = Molecule(symbols, geometry, alpha=alpha)
+    >>> args = [alpha]
+    >>> h = generate_hamiltonian(mol)(*args)
+    >>> h.terms[0]
+    tensor([ 0.29817879+0.j,  0.20813365+0.j,  0.20813365+0.j,
+             0.17860977+0.j,  0.04256036+0.j, -0.04256036+0.j,
+            -0.04256036+0.j,  0.04256036+0.j, -0.34724873+0.j,
+             0.13290293+0.j, -0.34724873+0.j,  0.17546329+0.j,
+             0.17546329+0.j,  0.13290293+0.j,  0.18470917+0.j], requires_grad=True)
+    """
+
+    def hamiltonian(*args):
+        r"""Compute the qubit hamiltonian.
+
+        Args:
+            args (array[array[float]]): initial values of the differentiable parameters
+
+        Returns:
+            Hamiltonian: the qubit Hamiltonian
+        """
+        h_ferm = generate_fermionic_hamiltonian(mol, cutoff, core, active)(*args)
+
+        ops = []
+
+        for n, t in enumerate(h_ferm[1]):
+
+            if len(t) == 0:
+                coeffs = np.array([h_ferm[0][n]])
+                ops = ops + [qml.Identity(0)]
+
+            elif len(t) == 2:
+                op = _generate_qubit_operator(t)
+                if op != 0:
+                    for i, o in enumerate(op[1]):
+                        if len(o) == 0:
+                            op[1][i] = qml.Identity(0)
+                        if len(o) == 1:
+                            op[1][i] = _return_pauli(o[0][1])(o[0][0])
+                        if len(o) > 1:
+                            k = qml.Identity(0)
+                            for o_ in o:
+                                k = k @ _return_pauli(o_[1])(o_[0])
+                            op[1][i] = k
+                    coeffs = np.concatenate([coeffs, np.array(op[0]) * h_ferm[0][n]])
+                    ops = ops + op[1]
+
+            elif len(t) == 4:
+                op = _generate_qubit_operator(t)
+                if op != 0:
+                    for i, o in enumerate(op[1]):
+                        if len(o) == 0:
+                            op[1][i] = qml.Identity(0)
+                        if len(o) == 1:
+                            op[1][i] = _return_pauli(o[0][1])(o[0][0])
+                        if len(o) > 1:
+                            k = qml.Identity(0)
+                            for o_ in o:
+                                k = k @ _return_pauli(o_[1])(o_[0])
+                            op[1][i] = k
+                    coeffs = np.concatenate([coeffs, np.array(op[0]) * h_ferm[0][n]])
+                    ops = ops + op[1]
+
+        h = qml.Hamiltonian(coeffs, ops, simplify=True)
+
+        return h
+
+    return hamiltonian
+
+
+def _generate_qubit_operator(op):
+    r"""Convert a fermionic operator to a qubit operator using the Jordan-Wigner mapping.
+
+    The one-body fermionic operator :math:`a_2^\dagger a_0` is constructed as [2, 0] and the
+    two-body operator :math:`a_4^\dagger a_3^\dagger a_2 a_1` is constructed as [4, 3, 2, 1].
+
+    Args:
+        op (list[int]): the fermionic operator
+
+    Returns
+        tuple(list[complex], list[list[int, str]]): list of coefficients and the qubit-operator terms
+
+    **Example**
+
+    >>> f  = [0, 0]
+    >>> q = _generate_qubit_operator(f)
+    >>> q
+    ([(0.5+0j), (-0.5+0j)], [[], [(0, 'Z')]]) # corresponds to :math:`\frac{1}{2}(I_0 - Z_0)`
+    """
+    if len(op) == 2:
+        op = [((op[0], 1), (op[1], 0))]
+
+    if len(op) == 4:
+        op = [((op[0], 1), (op[1], 1), (op[2], 0), (op[3], 0))]
+
+        if op[0][0][0] == op[0][1][0] or op[0][2][0] == op[0][3][0]:
+            return 0
+
+    for t in op:
+        for l in t:
+            z = [(index, "Z") for index in range(l[0])]
+            x = z + [(l[0], "X"), 0.5]
+
+            if l[1]:
+                y = z + [(l[0], "Y"), -0.5j]
+
+            else:
+                y = z + [(l[0], "Y"), 0.5j]
+
+            if t.index(l) == 0:
+                q = [x, y]
+            else:
+                m = []
+                for t1 in q:
+                    for t2 in [x, y]:
+                        q1, c1 = _pauli_mult(t1[:-1], t2[:-1])
+                        m.append(q1 + [c1 * t1[-1] * t2[-1]])
+                q = m
+
+    c = [p[-1] for p in q]
+    o = [p[:-1] for p in q]
+
+    for item in o:
+        k = [i for i, x in enumerate(o) if x == item]
+        if len(k) >= 2:
+            for j in k[::-1][:-1]:
+                del o[j]
+                c[k[0]] = c[k[0]] + c[j]
+                del c[j]
+
+    return c, o
+
+
+def _pauli_mult(p1, p2):
+    r"""Return the result of multiplication between two tensor products of Pauli operators.
+
+    The Pauli operator :math:`(P_0)` is denoted by [(0, 'P')], where :math:`P` represents
+    :math:`X`, :math:`Y` or :math:`Z`.
+
+    Args:
+        p1 (list[tuple[int, str]]): the first tensor product of Pauli operators
+        p2 (list[tuple[int, str]]): the second tensor product of Pauli operators
+
+    Returns
+        tuple(list[tuple[int, str]], complex): list of the Pauli operators and the coefficient
+
+    **Example**
+
+    >>> p1 = [(0, "X"), (1, "Y")],  # X_0 @ Y_1
+    >>> p2 = [(0, "X"), (2, "Y")],  # X_0 @ Y_2
+    >>> _pauli_mult(p1, p2)
+    ([(2, "Y"), (1, "Y")], 1.0) # p1 @ p2 = X_0 @ Y_1 @ X_0 @ Y_2
+    """
+    c = 1.0
+
+    t1 = [t[0] for t in p1]
+    t2 = [t[0] for t in p2]
+
+    k = []
+
+    for i in p1:
+        if i[0] in t1 and i[0] not in t2:
+            k.append((i[0], pauli_mult[i[1]]))
+        for j in p2:
+            if j[0] in t2 and j[0] not in t1:
+                k.append((j[0], pauli_mult[j[1]]))
+
+            if i[0] == j[0]:
+                if i[1] + j[1] in pauli_coeff:
+                    k.append((i[0], pauli_mult[i[1] + j[1]]))
+                    c = c * pauli_coeff[i[1] + j[1]]
+                else:
+                    k.append((i[0], pauli_mult[i[1] + j[1]]))
+
+    k = [i for i in k if "I" not in i[1]]
+
+    for item in k:
+        k_ = [i for i, x in enumerate(k) if x == item]
+        if len(k_) >= 2:
+            for j in k_[::-1][:-1]:
+                del k[j]
+
+    return k, c
+
+
+def _return_pauli(p):
+    r"""Return the PennyLane Pauli operator.
+
+    Args:
+        args (str): symbol representing the Pauli operator
+
+    Returns:
+        pennylane.ops: the PennyLane Pauli operator
+
+    **Example**
+
+    >>> _return_pauli('X')
+    qml.PauliX
+    """
+    if p == "X":
+        return qml.PauliX
+    if p == "Y":
+        return qml.PauliY
+
+    return qml.PauliZ
+
+
+pauli_mult = {
+    "XX": "I",
+    "YY": "I",
+    "ZZ": "I",
+    "ZX": "Y",
+    "XZ": "Y",
+    "ZY": "X",
+    "YZ": "X",
+    "XY": "Z",
+    "YX": "Z",
+    "IX": "X",
+    "IY": "Y",
+    "IZ": "Z",
+    "XI": "X",
+    "YI": "Y",
+    "ZI": "Z",
+    "I": "I",
+    "II": "I",
+    "X": "X",
+    "Y": "Y",
+    "Z": "Z",
+}
+
+pauli_coeff = {
+    "ZX": 1.0j,
+    "XZ": -1.0j,
+    "ZY": -1.0j,
+    "YZ": 1.0j,
+    "XY": 1.0j,
+    "YX": -1.0j,
+}
diff --git a/pennylane/hf/hartree_fock.py b/pennylane/hf/hartree_fock.py
new file mode 100644
index 000000000..ddf8a004c
--- /dev/null
+++ b/pennylane/hf/hartree_fock.py
@@ -0,0 +1,261 @@
+# Copyright 2018-2021 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""
+This module contains the functions needed for performing the self-consistent-field calculations.
+"""
+
+import itertools
+
+import autograd.numpy as anp
+from pennylane.hf.matrices import (
+    generate_core_matrix,
+    generate_overlap_matrix,
+    generate_repulsion_tensor,
+    molecular_density_matrix,
+)
+
+
+def generate_scf(mol, n_steps=50, tol=1e-8):
+    r"""Return a function that performs the self-consistent-field calculations.
+
+    In the Hartree-Fock method, molecular orbitals are typically constructed as a linear combination
+    of atomic orbitals
+
+    .. math::
+
+        \phi_i(r) = \sum_{\mu} C_{\mu i} \chi_{\mu}(r),
+
+    with coefficients :math:`C_{\mu i}` that are initially unknown. The self-consistent-field
+    iterations are performed to find a converged set of molecular orbital coefficients that minimize
+    the total energy of the molecular system. This optimization problem can be reduced to solving a
+    linear system of equations which are usually written as
+
+    .. math::
+
+        FC = SCE,
+
+    where :math:`E` is a diagonal matrix of eigenvalues, representing the molecular orbital
+    energies, :math:`C` is the matrix of molecular orbital coefficients, :math:`S` is the overlap
+    matrix and :math:`F` is the Fock matrix, which also depends on the coefficients. Fixing an
+    initial guess :math:`C_0`, the corresponding :math:`F_0` is built and the system
+    :math:`F_0C_0 = SC_0E` is solved to obtain a solution :math:`C_1`. This process is iteratively
+    repeated until the coefficients are converged.
+
+    The key step in in this process is constructing the Fock matrix which is defined as
+
+    .. math::
+
+        F = H + \frac{1}{2} J - K,
+
+    where :math:`H`, :math:`J` and :math:`K` are the core Hamiltonian matrix, Coulomb matrix and
+    exchange matrix, respectively. The entries of :math:`H` are computed from the electronic kinetic
+    energy and the electron-nuclear attraction integrals, which are integrals over atomic basis
+    functions. The elements of the :math:`J` and :math:`K` matrices are obtained from the Coulomb
+    and exchange integrals over the basis functions.
+
+    Following the procedure in
+    [`Lehtola et al. Molecules 2020, 25, 1218 <https://www.mdpi.com/1420-3049/25/5/1218>`_], we
+    express the molecular orbital coefficients in terms of a matrix :math:`X` as
+    :math:`C = X \tilde{C}` which gives the following transformed equation
+
+    .. math::
+
+         \tilde{F} \tilde{C} = \tilde{S} \tilde{C} E,
+
+    where :math:`\tilde{F} = X^T F X`, :math:`\tilde{S} = X^T S X` and :math:`S` is the overlap
+    matrix. We chose :math:`X` such that :math:`\tilde{S} = 1` as
+
+    .. math::
+
+        X = V \Lambda^{-1/2} V^T,
+
+    where :math:`V` and :math:`\Lambda` are the eigenvectors and eigenvalues of :math:`S`,
+    respectively. This gives the eigenvalue equation
+
+    .. math::
+
+         \tilde{F}\tilde{C} = \tilde{C}E,
+
+    which is solved with conventional methods iteratively.
+
+    Args:
+        mol (Molecule): the molecule object
+        n_steps (int): the number of iterations
+        tol (float): convergence tolerance
+
+    Returns:
+        function: function that performs the self-consistent-field calculations
+
+    **Example**
+
+    >>> symbols  = ['H', 'H']
+    >>> geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad = False)
+    >>> alpha = np.array([[3.42525091, 0.62391373, 0.1688554],
+    >>>                   [3.42525091, 0.62391373, 0.1688554]], requires_grad=True),
+    >>> mol = Molecule(symbols, geometry, alpha=alpha)
+    >>> args = [alpha]
+    >>> v_fock, coeffs, fock_matrix, h_core, repulsion_tensor = generate_hartree_fock(mol)(*args)
+    >>> v_fock
+    array([-0.67578019,  0.94181155])
+    """
+
+    def scf(*args):
+        r"""Perform the self-consistent-field iterations.
+
+        Args:
+            args (array[array[float]]): initial values of the differentiable parameters
+
+        Returns:
+            tuple(array[float]): eigenvalues of the Fock matrix, molecular orbital coefficients,
+            Fock matrix, core matrix
+        """
+        basis_functions = mol.basis_set
+        charges = mol.nuclear_charges
+        r = mol.coordinates
+        n_electron = mol.n_electrons
+
+        if r.requires_grad:
+            args_r = [[args[0][i]] * mol.n_basis[i] for i in range(len(mol.n_basis))]
+            args_ = [*args] + [anp.vstack(list(itertools.chain(*args_r)))]
+            repulsion_tensor = generate_repulsion_tensor(basis_functions)(*args_[1:])
+            s = generate_overlap_matrix(basis_functions)(*args_[1:])
+            h_core = generate_core_matrix(basis_functions, charges, r)(*args_)
+        else:
+            repulsion_tensor = generate_repulsion_tensor(basis_functions)(*args)
+            s = generate_overlap_matrix(basis_functions)(*args)
+            h_core = generate_core_matrix(basis_functions, charges, r)(*args)
+
+        s = s + anp.diag(anp.random.rand(len(s)) * 1.0e-12)
+
+        w, v = anp.linalg.eigh(s)
+        x = v @ anp.diag(anp.array([1 / anp.sqrt(i) for i in w])) @ v.T
+
+        eigvals, w_fock = anp.linalg.eigh(x.T @ h_core @ x)  # initial guess for the scf problem
+        coeffs = x @ w_fock
+
+        p = molecular_density_matrix(n_electron, coeffs)
+
+        for _ in range(n_steps):
+
+            j = anp.einsum("pqrs,rs->pq", repulsion_tensor, p)
+            k = anp.einsum("psqr,rs->pq", repulsion_tensor, p)
+
+            fock_matrix = h_core + 2 * j - k
+
+            eigvals, w_fock = anp.linalg.eigh(x.T @ fock_matrix @ x)
+
+            coeffs = x @ w_fock
+
+            p_update = molecular_density_matrix(n_electron, coeffs)
+
+            if anp.linalg.norm(p_update - p) <= tol:
+                break
+
+            p = p_update
+
+        return eigvals, coeffs, fock_matrix, h_core, repulsion_tensor
+
+    return scf
+
+
+def nuclear_energy(charges, r):
+    r"""Return a function that computes the nuclear-repulsion energy.
+
+    The nuclear-repulsion energy is computed as
+
+    .. math::
+
+        \sum_{i>j}^n \frac{q_i q_j}{r_{ij}},
+
+    where :math:`q`, :math:`r` and :math:`n` denote the nuclear charges (atomic numbers), nuclear
+    positions and the number of nuclei, respectively.
+
+    Args:
+        charges (list[int]): nuclear charges in atomic units
+        r (array[float]): nuclear positions
+
+    Returns:
+        function: function that computes the nuclear-repulsion energy
+
+    **Example**
+
+    >>> symbols  = ['H', 'F']
+    >>> geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 2.0]], requires_grad = True)
+    >>> mol = Molecule(symbols, geometry)
+    >>> args = [mol.coordinates]
+    >>> e = nuclear_energy(mol.nuclear_charges, mol.coordinates)(*args)
+    >>> print(e)
+    4.5
+    """
+
+    def nuclear(*args):
+        r"""Compute the nuclear-repulsion energy.
+
+        Args:
+            args (array[array[float]]): initial values of the differentiable parameters
+
+        Returns:
+            array[float]: nuclear-repulsion energy
+        """
+        if r.requires_grad:
+            coor = args[0]
+        else:
+            coor = r
+        e = 0
+        for i, r1 in enumerate(coor):
+            for j, r2 in enumerate(coor[i + 1 :]):
+                e = e + (charges[i] * charges[i + j + 1] / anp.linalg.norm(r1 - r2))
+        return e
+
+    return nuclear
+
+
+def hf_energy(mol):
+    r"""Return a function that computes the Hartree-Fock energy.
+
+    Args:
+        mol (Molecule): the molecule object
+
+    Returns:
+        function: function that computes the Hartree-Fock energy
+
+    **Example**
+
+    >>> symbols  = ['H', 'H']
+    >>> geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad = False)
+    >>> alpha = np.array([[3.42525091, 0.62391373, 0.1688554],
+    >>>                   [3.42525091, 0.62391373, 0.1688554]], requires_grad=True),
+    >>> mol = Molecule(symbols, geometry, alpha=alpha)
+    >>> args = [alpha]
+    >>> hf_energy(mol)(*args)
+    -1.065999461545263
+    """
+
+    def energy(*args):
+        r"""Compute the Hartree-Fock energy.
+
+        Args:
+            args (array[array[float]]): initial values of the differentiable parameters
+
+        Returns:
+            float: the Hartree-Fock energy
+        """
+        _, coeffs, fock_matrix, h_core, _ = generate_scf(mol)(*args)
+        e_rep = nuclear_energy(mol.nuclear_charges, mol.coordinates)(*args)
+        e_elec = anp.einsum(
+            "pq,qp", fock_matrix + h_core, molecular_density_matrix(mol.n_electrons, coeffs)
+        )
+        return e_elec + e_rep
+
+    return energy
diff --git a/pennylane/hf/integrals.py b/pennylane/hf/integrals.py
new file mode 100644
index 000000000..82b1457ff
--- /dev/null
+++ b/pennylane/hf/integrals.py
@@ -0,0 +1,853 @@
+# Copyright 2018-2021 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""
+This module contains the functions needed for computing integrals over basis functions.
+"""
+# pylint: disable= unbalanced-tuple-unpacking, too-many-arguments
+import autograd.numpy as anp
+import autograd.scipy as asp
+from scipy.special import factorial2 as fac2
+
+
+def primitive_norm(l, alpha):
+    r"""Compute the normalization constant for a primitive Gaussian function.
+
+    A Gaussian function centred at the position :math:`r = (x, y, z)` is defined as
+
+    .. math::
+
+        G = x^{l_x} y^{l_y} z^{l_z} e^{-\alpha r^2},
+
+    where :math:`l = (l_x, l_y, l_z)` defines the angular momentum quantum number and :math:`\alpha`
+    is the Gaussian function exponent. The normalization constant for this function is computed as
+
+    .. math::
+
+        N(l, \alpha) = (\frac{2\alpha}{\pi})^{3/4} \frac{(4 \alpha)^{(l_x + l_y + l_z)/2}}
+        {(2l_x-1)!! (2l_y-1)!! (2l_z-1)!!)^{1/2}}.
+
+    Args:
+        l (tuple[int]): angular momentum quantum number of the basis function
+        alpha (array[float]): exponent of the primitive Gaussian function
+
+    Returns:
+        array[float]: normalization coefficient
+
+    **Example**
+
+    >>> l = (0, 0, 0)
+    >>> alpha = np.array([3.425250914])
+    >>> n = gaussian_norm(l, alpha)
+    >>> print(n)
+    array([1.79444183])
+    """
+    lx, ly, lz = l
+    n = (
+        (2 * alpha / anp.pi) ** 0.75
+        * (4 * alpha) ** (sum(l) / 2)
+        / anp.sqrt(fac2(2 * lx - 1) * fac2(2 * ly - 1) * fac2(2 * lz - 1))
+    )
+    return n
+
+
+def contracted_norm(l, alpha, a):
+    r"""Compute the normalization constant for a contracted Gaussian function.
+
+    A contracted Gaussian function is defined as
+
+    .. math::
+
+        \psi = a_1 G_1 + a_2 G_2 + a_3 G_3,
+
+    where :math:`a` denotes the contraction coefficients and :math:`G` is a primitive Gaussian function. The
+    normalization constant for this function is computed as
+
+    .. math::
+
+        N(l, \alpha, a) = [\frac{\pi^{3/2}(2l_x-1)!! (2l_y-1)!! (2l_z-1)!!}{2^{l_x + l_y + l_z}}
+        \sum_{i,j} \frac{a_i a_j}{(\alpha_i + \alpha_j)^{{l_x + l_y + l_z+3/2}}}]^{-1/2}
+
+    where :math:`l` and :math:`\alpha` denote the angular momentum quantum number and the exponent
+    of the Gaussian function, respectively.
+
+    Args:
+        l (tuple[int]): angular momentum quantum number of the primitive Gaussian functions
+        alpha (array[float]): exponents of the primitive Gaussian functions
+        a (array[float]): coefficients of the contracted Gaussian functions
+
+    Returns:
+        array[float]: normalization coefficient
+
+    **Example**
+
+    >>> l = (0, 0, 0)
+    >>> alpha = np.array([3.425250914, 0.6239137298, 0.168855404])
+    >>> a = np.array([1.79444183, 0.50032649, 0.18773546])
+    >>> n = contracted_norm(l, alpha, a)
+    >>> print(n)
+    0.39969026908800853
+    """
+    lx, ly, lz = l
+    c = anp.pi ** 1.5 / 2 ** sum(l) * fac2(2 * lx - 1) * fac2(2 * ly - 1) * fac2(2 * lz - 1)
+    s = (
+        (a.reshape(len(a), 1) * a) / ((alpha.reshape(len(alpha), 1) + alpha) ** (sum(l) + 1.5))
+    ).sum()
+    n = 1 / anp.sqrt(c * s)
+    return n
+
+
+def _generate_params(params, args):
+    """Generate basis set parameters. The default values are used for the non-differentiable
+    parameters and the user-defined values are used for the differentiable ones.
+
+    Args:
+        params (list(array[float])): default values of the basis set parameters
+        args (list(array[float])): initial values of the differentiable basis set parameters
+
+    Returns:
+        list(array[float]): basis set parameters
+    """
+    basis_params = []
+    c = 0
+    for p in params:
+        if p.requires_grad:
+            basis_params.append(args[c])
+            c += 1
+        else:
+            basis_params.append(p)
+    return basis_params
+
+
+def expansion(la, lb, ra, rb, alpha, beta, t):
+    r"""Compute Hermite Gaussian expansion coefficients recursively for two Gaussian functions.
+
+    An overlap distribution, which defines the product of two Gaussians, can be written as a Hermite
+    expansion as [`Helgaker (1995) p798 <https://www.worldscientific.com/doi/abs/10.1142/9789812832115_0001>`_]
+
+    .. math::
+
+        \Omega_{ij} = \sum_{t=0}^{i+j} E_t^{ij} \Lambda_t,
+
+    where :math:`\Lambda` is a Hermite polynomial of degree :math:`t`, :math:`E` denotes the expansion
+    coefficients, :math:`\Omega_{ij} = G_i G_j`, and :math:`G` is a Gaussian function. The overlap
+    integral between two Gaussian functions can be simply computed by integrating over the overlap
+    distribution which requires obtaining the expansion coefficients. This can be done recursively
+    as [`Helgaker (1995) p799 <https://www.worldscientific.com/doi/abs/10.1142/9789812832115_0001>`_]
+
+    .. math::
+
+        E_t^{i+1,j} = \frac{1}{2p} E_{t-1}^{ij} - \frac{qr}{\alpha} E_{t}^{ij} + (t+1) E_{t+1}^{ij},
+
+    and
+
+    .. math::
+
+        E_t^{i,j+1} = \frac{1}{2p} E_{t-1}^{ij} + \frac{qr}{\beta} E_{t}^{ij} + (t+1) E_{t+1}^{ij},
+
+    where :math:`p = \alpha + \beta` and :math:`q = \alpha \beta / (\alpha + \beta)` are computed
+    from the Gaussian exponents :math:`\alpha, \beta` and the position :math:`r` is computed as
+    :math:`r = r_\alpha - r_\beta`. The starting coefficient is
+
+    .. math::
+
+        E_0^{00} = e^{-qr^2},
+
+    and :math:`E_t^{ij} = 0` is :math:`t < 0` or :math:`t > (i+j)`.
+
+    Args:
+        la (integer): angular momentum component for the first Gaussian function
+        lb (integer): angular momentum component for the second Gaussian function
+        ra (float): position component of the the first Gaussian function
+        rb (float): position component of the the second Gaussian function
+        alpha (array[float]): exponent of the first Gaussian function
+        beta (array[float]): exponent of the second Gaussian function
+        t (integer): number of nodes in the Hermite Gaussian
+
+    Returns:
+        array[float]: expansion coefficients for each Gaussian combination
+
+    **Example**
+
+    >>> la, lb = 0, 0
+    >>> ra, rb = 0.0, 0.0
+    >>> alpha = np.array([3.42525091])
+    >>> beta =  np.array([3.42525091])
+    >>> t = 0
+    >>> c = expansion(la, lb, ra, rb, alpha, beta, t)
+    >>> c
+    array([1.])
+    """
+    p = anp.array(alpha + beta)
+    q = anp.array(alpha * beta / p)
+    r = anp.array(ra - rb)
+
+    if la == lb == t == 0:
+        return anp.exp(-q * r ** 2)
+
+    if t < 0 or t > (la + lb):
+        return 0.0
+
+    if lb == 0:
+        return (
+            (1 / (2 * p)) * expansion(la - 1, lb, ra, rb, alpha, beta, t - 1)
+            - (q * r / alpha) * expansion(la - 1, lb, ra, rb, alpha, beta, t)
+            + (t + 1) * expansion(la - 1, lb, ra, rb, alpha, beta, t + 1)
+        )
+
+    return (
+        (1 / (2 * p)) * expansion(la, lb - 1, ra, rb, alpha, beta, t - 1)
+        + (q * r / beta) * expansion(la, lb - 1, ra, rb, alpha, beta, t)
+        + (t + 1) * expansion(la, lb - 1, ra, rb, alpha, beta, t + 1)
+    )
+
+
+def gaussian_overlap(la, lb, ra, rb, alpha, beta):
+    r"""Compute overlap integral for two primitive Gaussian functions.
+
+    The overlap integral between two Gaussian functions denoted by :math:`a` and :math:`b` can be
+    computed as [`Helgaker (1995) p803 <https://www.worldscientific.com/doi/abs/10.1142/9789812832115_0001>`_]:
+
+    .. math::
+
+        S_{ab} = E^{ij} E^{kl} E^{mn} \left (\frac{\pi}{p}  \right )^{3/2},
+
+    where :math:`E` is a coefficient that can be computed recursively, :math:`i-n` are the angular
+    momentum quantum numbers corresponding to different Cartesian components and :math:`p` is
+    computed from the exponents of the two Gaussian functions as :math:`p = \alpha + \beta`.
+
+    Args:
+        la (integer): angular momentum for the first Gaussian function
+        lb (integer): angular momentum for the second Gaussian function
+        ra (float): position vector of the the first Gaussian function
+        rb (float): position vector of the the second Gaussian function
+        alpha (array[float]): exponent of the first Gaussian function
+        beta (array[float]): exponent of the second Gaussian function
+
+    Returns:
+        array[float]: overlap integral between primitive Gaussian functions
+
+    **Example**
+
+    >>> la, lb = (0, 0, 0), (0, 0, 0)
+    >>> ra, rb = np.array(([0., 0., 0.]), np.array(([0., 0., 0.])
+    >>> alpha = np.array([np.pi/2])
+    >>> beta = np.array([np.pi/2])
+    >>> o = gaussian_overlap(la, lb, ra, rb, alpha, beta)
+    >>> o
+    array([1.])
+    """
+    p = alpha + beta
+    s = 1.0
+    for i in range(3):
+        s = s * anp.sqrt(anp.pi / p) * expansion(la[i], lb[i], ra[i], rb[i], alpha, beta, 0)
+    return s
+
+
+def generate_overlap(basis_a, basis_b):
+    r"""Return a function that computes the overlap integral for two contracted Gaussian functions.
+
+    Args:
+        basis_a (BasisFunction): first basis function
+        basis_b (BasisFunction): second basis function
+
+    Returns:
+        function: function that computes the overlap integral
+
+    **Example**
+
+    >>> symbols  = ['H', 'H']
+    >>> geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad = False)
+    >>> mol = Molecule(symbols, geometry)
+    >>> args = []
+    >>> generate_overlap(mol.basis_set[0], mol.basis_set[0])(*args)
+    1.0
+    """
+
+    def overlap_integral(*args):
+        r"""Normalize and compute the overlap integral for two contracted Gaussian functions.
+
+        Args:
+            args (array[float]): initial values of the differentiable parameters
+
+        Returns:
+            array[float]: the overlap integral between two contracted Gaussian orbitals
+        """
+
+        args_a = [i[0] for i in args]
+        args_b = [i[1] for i in args]
+
+        alpha, ca, ra = _generate_params(basis_a.params, args_a)
+        beta, cb, rb = _generate_params(basis_b.params, args_b)
+
+        ca = ca * primitive_norm(basis_a.l, alpha)
+        cb = cb * primitive_norm(basis_b.l, beta)
+
+        na = contracted_norm(basis_a.l, alpha, ca)
+        nb = contracted_norm(basis_b.l, beta, cb)
+
+        return (
+            na
+            * nb
+            * (
+                (ca[:, anp.newaxis] * cb)
+                * gaussian_overlap(basis_a.l, basis_b.l, ra, rb, alpha[:, anp.newaxis], beta)
+            ).sum()
+        )
+
+    return overlap_integral
+
+
+def _diff2(i, j, ri, rj, alpha, beta):
+    r"""Compute the second order differentiated integral needed for evaluating a kinetic integral.
+
+    The second-order integral :math:`D_{ij}^2`, where :math:`i` and :math:`j` denote angular
+    momentum components of Gaussian functions, is computed from overlap integrals :math:`S` and the
+    Gaussian exponent :math:`\beta` as
+    [`Helgaker (1995) p804 <https://www.worldscientific.com/doi/abs/10.1142/9789812832115_0001>`_]:
+
+    .. math::
+
+        D_{ij}^2 = j(j-1)S_{i,j-2}^0 - 2\beta(2j+1)S_{i,j}^0 + 4\beta^2 S_{i,j+2}^0.
+
+    Args:
+        i (integer): angular momentum component for the first Gaussian function
+        j (integer): angular momentum component for the second Gaussian function
+        ri (array[float]): position component of the the first Gaussian function
+        rj (array[float]): position component of the the second Gaussian function
+        alpha (array[float]): exponent of the first Gaussian function
+        beta (array[float]): exponent of the second Gaussian function
+
+    Returns:
+        array[float]: second-order differentiated integral between two Gaussian functions
+    """
+    p = alpha + beta
+
+    d1 = j * (j - 1) * anp.sqrt(anp.pi / p) * expansion(i, j - 2, ri, rj, alpha, beta, 0)
+    d2 = -2 * beta * (2 * j + 1) * anp.sqrt(anp.pi / p) * expansion(i, j, ri, rj, alpha, beta, 0)
+    d3 = 4 * beta ** 2 * anp.sqrt(anp.pi / p) * expansion(i, j + 2, ri, rj, alpha, beta, 0)
+
+    return d1 + d2 + d3
+
+
+def gaussian_kinetic(la, lb, ra, rb, alpha, beta):
+    r"""Compute the kinetic integral for two primitive Gaussian functions.
+
+    The kinetic integral between two Gaussian functions denoted by :math:`a` and :math:`b` is
+    computed as
+    [`Helgaker (1995) p805 <https://www.worldscientific.com/doi/abs/10.1142/9789812832115_0001>`_]:
+
+    .. math::
+
+        T_{ab} = -\frac{1}{2} \left ( D_{ij}^2 D_{kl}^0 D_{mn}^0 + D_{ij}^0 D_{kl}^2 D_{mn}^0 + D_{ij}^0 D_{kl}^0 D_{mn}^2\right ),
+
+    where :math:`D_{ij}^0 = S_{ij}^0` is an overlap integral and :math:`D_{ij}^2` is computed from
+    overlap integrals :math:`S` and the Gaussian exponent :math:`\beta` as
+
+    .. math::
+
+        D_{ij}^2 = j(j-1)S_{i,j-2}^0 - 2\beta(2j+1)S_{i,j}^0 + 4\beta^2 S_{i,j+2}^0.
+
+    Args:
+        la (tuple[int]): angular momentum for the first Gaussian function
+        lb (tuple[int]): angular momentum for the second Gaussian function
+        ra (array[float]): position vector of the the first Gaussian function
+        rb (array[float]): position vector of the the second Gaussian function
+        alpha (array[float]): exponent of the first Gaussian function
+        beta (array[float]): exponent of the second Gaussian function
+
+    Returns:
+        array[float]: kinetic integral between two Gaussian functions
+
+    **Example**
+
+    >>> la, lb = (0, 0, 0), (0, 0, 0)
+    >>> ra, rb = np.array(([0., 0., 0.]), np.array(([0., 0., 0.])
+    >>> alpha = np.array([np.pi/2])
+    >>> beta = np.array([np.pi/2])
+    >>> t = gaussian_kinetic(la, lb, ra, rb, alpha, beta)
+    >>> t
+    array([2.35619449])
+    """
+
+    p = alpha + beta
+
+    t1 = (
+        _diff2(la[0], lb[0], ra[0], rb[0], alpha, beta)
+        * anp.sqrt(anp.pi / p)
+        * expansion(la[1], lb[1], ra[1], rb[1], alpha, beta, 0)
+        * anp.sqrt(anp.pi / p)
+        * expansion(la[2], lb[2], ra[2], rb[2], alpha, beta, 0)
+    )
+
+    t2 = (
+        anp.sqrt(anp.pi / p)
+        * expansion(la[0], lb[0], ra[0], rb[0], alpha, beta, 0)
+        * _diff2(la[1], lb[1], ra[1], rb[1], alpha, beta)
+        * anp.sqrt(anp.pi / p)
+        * expansion(la[2], lb[2], ra[2], rb[2], alpha, beta, 0)
+    )
+
+    t3 = (
+        anp.sqrt(anp.pi / p)
+        * expansion(la[0], lb[0], ra[0], rb[0], alpha, beta, 0)
+        * anp.sqrt(anp.pi / p)
+        * expansion(la[1], lb[1], ra[1], rb[1], alpha, beta, 0)
+        * _diff2(la[2], lb[2], ra[2], rb[2], alpha, beta)
+    )
+
+    return -0.5 * (t1 + t2 + t3)
+
+
+def generate_kinetic(basis_a, basis_b):
+    r"""Return a function that computes the kinetic integral for two contracted Gaussian functions.
+
+    Args:
+        basis_a (BasisFunction): first basis function
+        basis_b (BasisFunction): second basis function
+
+    Returns:
+        function: function that computes the kinetic integral
+
+    **Example**
+
+    >>> symbols  = ['H', 'H']
+    >>> geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad = False)
+    >>> alpha = np.array([[3.425250914, 0.6239137298, 0.168855404],
+    >>>                   [3.425250914, 0.6239137298, 0.168855404]], requires_grad = True)
+    >>> mol = hf.Molecule(symbols, geometry, alpha=alpha)
+    >>> args = [mol.alpha]
+    >>> generate_kinetic(mol.basis_set[0], mol.basis_set[1])(*args)
+    0.38325367405312843
+    """
+
+    def kinetic_integral(*args):
+        r"""Compute the kinetic integral for two contracted Gaussian functions.
+
+        Args:
+            args (array[float]): initial values of the differentiable parameters
+
+        Returns:
+            array[float]: the kinetic integral between two contracted Gaussian orbitals
+        """
+        args_a = [i[0] for i in args]
+        args_b = [i[1] for i in args]
+
+        alpha, ca, ra = _generate_params(basis_a.params, args_a)
+        beta, cb, rb = _generate_params(basis_b.params, args_b)
+
+        ca = ca * primitive_norm(basis_a.l, alpha)
+        cb = cb * primitive_norm(basis_b.l, beta)
+
+        na = contracted_norm(basis_a.l, alpha, ca)
+        nb = contracted_norm(basis_b.l, beta, cb)
+
+        return (
+            na
+            * nb
+            * (
+                (ca[:, anp.newaxis] * cb)
+                * gaussian_kinetic(basis_a.l, basis_b.l, ra, rb, alpha[:, anp.newaxis], beta)
+            ).sum()
+        )
+
+    return kinetic_integral
+
+
+def _boys(n, t):
+    r"""Evaluate the Boys function.
+
+    The :math:`n`-th order `Boys function <https://arxiv.org/abs/2107.01488>`_ is defined as
+
+    .. math::
+
+        F_n(t) = \int_{0}^{1}x^{2n} e^{-tx^2}dx.
+
+    The Boys function is related to the lower incomplete Gamma
+    `function <https://en.wikipedia.org/wiki/Incomplete_gamma_function>`_, :math:`\gamma`, as
+
+    .. math::
+
+        F_n(t) = \frac{1}{2t^{n + 0.5}} \gamma(n + 0.5, t),
+
+    where
+
+    .. math::
+
+        \gamma(m, t) = \int_{0}^{t} x^{m-1} e^{-x} dx.
+
+    Args:
+        n (float): order of the Boys function
+        t (float): exponent of the Boys function
+
+    Returns:
+        float: value of the Boys function
+    """
+    if t == 0.0:
+        return 1 / (2 * n + 1)
+
+    return asp.special.gammainc(n + 0.5, t) * asp.special.gamma(n + 0.5) / (2 * t ** (n + 0.5))
+
+
+def _hermite_coulomb(t, u, v, n, p, dr):
+    """Evaluate the Hermite integral needed to compute the nuclear attraction and electron repulsion
+    integrals.
+
+    These integrals are computed recursively starting from the Boys function
+    [`Helgaker (1995) p817 <https://www.worldscientific.com/doi/abs/10.1142/9789812832115_0001>`_]:
+
+
+    .. math::
+
+        R_{000}^n = (-2p)^n F_n(pR_{CP}^2),
+
+    where :math:`F_n` is the Boys function, :math:`p` is computed from the exponents of the two
+    Gaussian functions as :math:`p = \alpha + \beta`, and :math:`R_{CP}` is the distance between the
+    center of the composite Gaussian centered at :math:`P` and the electrostatic potential generated
+    by a nucleus at :math:`C`. The following recursive equations are used to compute the
+    higher-order Hermite integrals
+
+    .. math::
+
+        R_{t+1, u, v}^n = t R_{t-1, u, v}^{n+1} + x R_{t, u, v}^{n+1},
+
+        R_{t, u+1, v}^n = u R_{t, u-1, v}^{n+1} + y R_{t, u, v}^{n+1},
+
+        R_{t, u, v+1}^n = v R_{t, u, v-1}^{n+1} + z R_{t, u, v}^{n+1},
+
+    where :math:`x`, :math:`y` and :math:`z` are the Cartesian components of :math:`R_{CP}`.
+
+    Args:
+        t (integer): order of Hermite derivative in x
+        u (integer): order of Hermite derivative in y
+        v (float): order of Hermite derivative in z
+        n (integer): order of the Boys function
+        p (float): sum of the Gaussian exponents
+        dr (array[float]): distance between the center of the composite Gaussian and the nucleus
+
+    Returns:
+        array[float]: value of the Hermite integral
+    """
+    x, y, z = dr[0], dr[1], dr[2]
+    T = p * (dr ** 2).sum(axis=0)
+    r = 0
+
+    if t == u == v == 0:
+        f = []
+        for term in T.flatten():
+            f.append(_boys(n, term))
+        return ((-2 * p) ** n) * anp.array(f).reshape(T.shape)
+
+    if t == u == 0:
+        if v > 1:
+            r = r + (v - 1) * _hermite_coulomb(t, u, v - 2, n + 1, p, dr)
+        r = r + z * _hermite_coulomb(t, u, v - 1, n + 1, p, dr)
+        return r
+
+    if t == 0:
+        if u > 1:
+            r = r + (u - 1) * _hermite_coulomb(t, u - 2, v, n + 1, p, dr)
+        r = r + y * _hermite_coulomb(t, u - 1, v, n + 1, p, dr)
+        return r
+
+    if t > 1:
+        r = r + (t - 1) * _hermite_coulomb(t - 2, u, v, n + 1, p, dr)
+    r = r + x * _hermite_coulomb(t - 1, u, v, n + 1, p, dr)
+    return r
+
+
+def nuclear_attraction(la, lb, ra, rb, alpha, beta, r):
+    r"""Compute nuclear attraction integral between primitive Gaussian functions.
+
+    The nuclear attraction integral between two Gaussian functions denoted by :math:`a` and
+    :math:`b` can be computed as
+    [`Helgaker (1995) p820 <https://www.worldscientific.com/doi/abs/10.1142/9789812832115_0001>`_]
+
+    .. math::
+
+        V_{ab} = \frac{2\pi}{p} \sum_{tuv} E_t^{ij} E_u^{kl} E_v^{mn} R_{tuv},
+
+    where :math:`E` and :math:`R` represent the Hermite Gaussian expansion coefficients and the
+    Hermite Coulomb integral, respectively. The sum goes over :math:`i + j + 1`, :math:`k + l + 1`
+    and :math:`m + n + 1` for :math:`t`, :math:`u` and :math:`v`, respectively, and :math:`p` is
+    computed from the exponents of the two Gaussian functions as :math:`p = \alpha + \beta`.
+
+    Args:
+        la (tuple[int]): angular momentum for the first Gaussian function
+        lb (tuple[int]): angular momentum for the second Gaussian function
+        ra (array[float]): position vector of the the first Gaussian function
+        rb (array[float]): position vector of the the second Gaussian function
+        alpha (array[float]): exponent of the first Gaussian function
+        beta (array[float]): exponent of the second Gaussian function
+        r (array[float]): position vector of nucleus
+
+    Returns:
+        array[float]: nuclear attraction integral between two Gaussian functions
+    """
+    l1, m1, n1 = la
+    l2, m2, n2 = lb
+    p = alpha + beta
+    rgp = (alpha * ra[:, anp.newaxis, anp.newaxis] + beta * rb[:, anp.newaxis, anp.newaxis]) / (
+        alpha + beta
+    )
+    dr = rgp - anp.array(r)[:, anp.newaxis, anp.newaxis]
+
+    a = 0.0
+    for t in range(l1 + l2 + 1):
+        for u in range(m1 + m2 + 1):
+            for v in range(n1 + n2 + 1):
+                a = a + expansion(l1, l2, ra[0], rb[0], alpha, beta, t) * expansion(
+                    m1, m2, ra[1], rb[1], alpha, beta, u
+                ) * expansion(n1, n2, ra[2], rb[2], alpha, beta, v) * _hermite_coulomb(
+                    t, u, v, 0, p, dr
+                )
+    a = a * 2 * anp.pi / p
+    return a
+
+
+def generate_attraction(r, basis_a, basis_b):
+    r"""Return a function that computes the nuclear attraction integral for two contracted Gaussian
+    functions.
+
+    Args:
+        r (array[float]): position vector of nucleus
+        basis_a (BasisFunction): first basis function
+        basis_b (BasisFunction): second basis function
+
+    Returns:
+        function: function that computes the electron-nuclear attraction integral
+
+    **Example**
+
+    >>> symbols  = ['H', 'H']
+    >>> geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad = False)
+    >>> alpha = np.array([[3.425250914, 0.6239137298, 0.168855404],
+    >>>                   [3.425250914, 0.6239137298, 0.168855404]], requires_grad = True)
+    >>> mol = Molecule(symbols, geometry, alpha=alpha)
+    >>> basis_a = mol.basis_set[0]
+    >>> basis_b = mol.basis_set[1]
+    >>> args = [mol.alpha]
+    >>> generate_attraction(geometry[0], basis_a, basis_b)(*args)
+    0.801208332328965
+    """
+
+    def attraction_integral(*args):
+        r"""Compute the electron-nuclear attraction integral for two contracted Gaussian functions.
+
+        Args:
+            args (array[float]): initial values of the differentiable parameters
+
+        Returns:
+            array[float]: the electron-nuclear attraction integral
+        """
+        if r.requires_grad:
+            coor = args[0]
+            args_a = [i[0] for i in args[1:]]
+            args_b = [i[1] for i in args[1:]]
+        else:
+            coor = r
+            args_a = [i[0] for i in args]
+            args_b = [i[1] for i in args]
+
+        alpha, ca, ra = _generate_params(basis_a.params, args_a)
+        beta, cb, rb = _generate_params(basis_b.params, args_b)
+
+        ca = ca * primitive_norm(basis_a.l, alpha)
+        cb = cb * primitive_norm(basis_b.l, beta)
+
+        na = contracted_norm(basis_a.l, alpha, ca)
+        nb = contracted_norm(basis_b.l, beta, cb)
+
+        v = (
+            na
+            * nb
+            * (
+                (ca * cb[:, anp.newaxis])
+                * nuclear_attraction(
+                    basis_a.l, basis_b.l, ra, rb, alpha, beta[:, anp.newaxis], coor
+                )
+            ).sum()
+        )
+        return v
+
+    return attraction_integral
+
+
+def electron_repulsion(la, lb, lc, ld, ra, rb, rc, rd, alpha, beta, gamma, delta):
+    r"""Compute the electron-electron repulsion integral between four primitive Gaussian functions.
+
+    The electron repulsion integral between four Gaussian functions denoted by :math:`a`, :math:`b`
+    , :math:`c` and :math:`d` can be computed as
+    [`Helgaker (1995) p820 <https://www.worldscientific.com/doi/abs/10.1142/9789812832115_0001>`_]
+
+    .. math::
+
+        g_{abcd} = \frac{2\pi^{5/2}}{pq\sqrt{p+q}} \sum_{tuv} E_t^{o_a o_b} E_u^{m_a m_b}
+        E_v^{n_a n_b} \sum_{rsw} (-1)^{r+s+w} E_r^{o_c o_d} E_s^{m_c m_d} E_w^{n_c n_d}
+        R_{t+r, u+s, v+w},
+
+    where :math:`E` and :math:`R` are the Hermite Gaussian expansion coefficients and the
+    Hermite Coulomb integral, respectively. The sums go over the angular momentum quantum numbers
+    :math:`o_i + o_j + 1`, :math:`m_i + m_j + 1` and :math:`n_i + n_j + 1` respectively for
+    :math:`t, u, v` and :math:`r, s, w`. The exponents of the Gaussian functions are used to compute
+    :math:`p` and :math:`q` as :math:`p = \alpha + \beta` and :math:`q = \gamma + \delta`.
+
+    Args:
+        la (tuple[int]): angular momentum for the first Gaussian function
+        lb (tuple[int]): angular momentum for the second Gaussian function
+        lc (tuple[int]): angular momentum for the third Gaussian function
+        ld (tuple[int]): angular momentum for the forth Gaussian function
+        ra (array[float]): position vector of the the first Gaussian function
+        rb (array[float]): position vector of the the second Gaussian function
+        rc (array[float]): position vector of the the third Gaussian function
+        rd (array[float]): position vector of the the forth Gaussian function
+        alpha (array[float]): exponent of the first Gaussian function
+        beta (array[float]): exponent of the second Gaussian function
+        gamma (array[float]): exponent of the third Gaussian function
+        delta (array[float]): exponent of the forth Gaussian function
+
+    Returns:
+        array[float]: electron-electron repulsion integral between four Gaussian functions
+    """
+    l1, m1, n1 = la
+    l2, m2, n2 = lb
+    l3, m3, n3 = lc
+    l4, m4, n4 = ld
+
+    p = alpha + beta
+    q = gamma + delta
+
+    p_ab = (
+        alpha * ra[:, anp.newaxis, anp.newaxis, anp.newaxis, anp.newaxis]
+        + beta * rb[:, anp.newaxis, anp.newaxis, anp.newaxis, anp.newaxis]
+    ) / (alpha + beta)
+
+    p_cd = (
+        gamma * rc[:, anp.newaxis, anp.newaxis, anp.newaxis, anp.newaxis]
+        + delta * rd[:, anp.newaxis, anp.newaxis, anp.newaxis, anp.newaxis]
+    ) / (gamma + delta)
+
+    g = 0.0
+    for t in range(l1 + l2 + 1):
+        for u in range(m1 + m2 + 1):
+            for v in range(n1 + n2 + 1):
+                for r in range(l3 + l4 + 1):
+                    for s in range(m3 + m4 + 1):
+                        for w in range(n3 + n4 + 1):
+                            g = g + expansion(l1, l2, ra[0], rb[0], alpha, beta, t) * expansion(
+                                m1, m2, ra[1], rb[1], alpha, beta, u
+                            ) * expansion(n1, n2, ra[2], rb[2], alpha, beta, v) * expansion(
+                                l3, l4, rc[0], rd[0], gamma, delta, r
+                            ) * expansion(
+                                m3, m4, rc[1], rd[1], gamma, delta, s
+                            ) * expansion(
+                                n3, n4, rc[2], rd[2], gamma, delta, w
+                            ) * (
+                                (-1) ** (r + s + w)
+                            ) * _hermite_coulomb(
+                                t + r, u + s, v + w, 0, (p * q) / (p + q), p_ab - p_cd
+                            )
+
+    g = g * 2 * (anp.pi ** 2.5) / (p * q * anp.sqrt(p + q))
+
+    return g
+
+
+def generate_repulsion(basis_a, basis_b, basis_c, basis_d):
+    r"""Return a function that computes the electron-electron repulsion integral for four contracted
+    Gaussian functions.
+
+    Args:
+        basis_a (BasisFunction): first basis function
+        basis_b (BasisFunction): second basis function
+        basis_c (BasisFunction): third basis function
+        basis_d (BasisFunction): fourth basis function
+    Returns:
+        function: function that computes the electron repulsion integral
+
+    **Example**
+
+    >>> symbols  = ['H', 'H']
+    >>> geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad = False)
+    >>> alpha = np.array([[3.425250914, 0.6239137298, 0.168855404],
+    >>>                   [3.425250914, 0.6239137298, 0.168855404],
+    >>>                   [3.425250914, 0.6239137298, 0.168855404],
+    >>>                   [3.425250914, 0.6239137298, 0.168855404]], requires_grad = True)
+    >>> mol = Molecule(symbols, geometry, alpha=alpha)
+    >>> basis_a = mol.basis_set[0]
+    >>> basis_b = mol.basis_set[1]
+    >>> args = [mol.alpha]
+    >>> generate_repulsion(basis_a, basis_b, basis_a, basis_b)(*args)
+    0.45590152106593573
+    """
+
+    def repulsion_integral(*args):
+        r"""Compute the electron-electron repulsion integral for four contracted Gaussian functions.
+
+        Args:
+            args (array[float]): initial values of the differentiable parameters
+
+        Returns:
+            array[float]: the electron repulsion integral between four contracted Gaussian functions
+        """
+        args_a = [i[0] for i in args]
+        args_b = [i[1] for i in args]
+        args_c = [i[2] for i in args]
+        args_d = [i[3] for i in args]
+
+        alpha, ca, ra = _generate_params(basis_a.params, args_a)
+        beta, cb, rb = _generate_params(basis_b.params, args_b)
+        gamma, cc, rc = _generate_params(basis_c.params, args_c)
+        delta, cd, rd = _generate_params(basis_d.params, args_d)
+
+        ca = ca * primitive_norm(basis_a.l, alpha)
+        cb = cb * primitive_norm(basis_b.l, beta)
+        cc = cc * primitive_norm(basis_c.l, gamma)
+        cd = cd * primitive_norm(basis_d.l, delta)
+
+        n1 = contracted_norm(basis_a.l, alpha, ca)
+        n2 = contracted_norm(basis_b.l, beta, cb)
+        n3 = contracted_norm(basis_c.l, gamma, cc)
+        n4 = contracted_norm(basis_d.l, delta, cd)
+
+        e = (
+            n1
+            * n2
+            * n3
+            * n4
+            * (
+                (
+                    ca
+                    * cb[:, anp.newaxis]
+                    * cc[:, anp.newaxis, anp.newaxis]
+                    * cd[:, anp.newaxis, anp.newaxis, anp.newaxis]
+                )
+                * electron_repulsion(
+                    basis_a.l,
+                    basis_b.l,
+                    basis_c.l,
+                    basis_d.l,
+                    ra,
+                    rb,
+                    rc,
+                    rd,
+                    alpha,
+                    beta[:, anp.newaxis],
+                    gamma[:, anp.newaxis, anp.newaxis],
+                    delta[:, anp.newaxis, anp.newaxis, anp.newaxis],
+                )
+            ).sum()
+        )
+        return e
+
+    return repulsion_integral
diff --git a/pennylane/hf/matrices.py b/pennylane/hf/matrices.py
new file mode 100644
index 000000000..614c15844
--- /dev/null
+++ b/pennylane/hf/matrices.py
@@ -0,0 +1,346 @@
+# Copyright 2018-2021 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""
+This module contains the functions needed for computing matrices.
+"""
+# pylint: disable= too-many-branches
+import autograd.numpy as anp
+from pennylane.hf.integrals import (
+    generate_attraction,
+    generate_kinetic,
+    generate_overlap,
+    generate_repulsion,
+)
+
+
+def molecular_density_matrix(n_electron, c):
+    r"""Compute the molecular density matrix.
+
+    The density matrix :math:`P` is computed from the molecular orbital coefficients :math:`C` as
+
+    .. math::
+
+        P_{\mu \nu} = \sum_{i=1}^{N} C_{\mu i} C_{\nu i},
+
+    where :math:`N = N_{electrons} / 2` is the number of occupied orbitals. Note that the total
+    density matrix is the sum of the :math:`\alpha` and :math:`\beta` density
+    matrices, :math:`P = P^{\alpha} + P^{\beta}`.
+
+    Args:
+        n_electron (integer): number of electrons
+        c (array[array[float]]): molecular orbital coefficients
+
+    Returns:
+        array[array[float]]: density matrix
+
+    **Example**
+
+    >>> c = np.array([[-0.54828771,  1.21848441], [-0.54828771, -1.21848441]])
+    >>> n_electron = 2
+    >>> density_matrix(n_electron, c)
+    array([[0.30061941, 0.30061941], [0.30061941, 0.30061941]])
+    """
+    p = anp.dot(c[:, : n_electron // 2], anp.conjugate(c[:, : n_electron // 2]).T)
+    return p
+
+
+def generate_overlap_matrix(basis_functions):
+    r"""Return a function that computes the overlap matrix for a given set of basis functions.
+
+    Args:
+        basis_functions (list[BasisFunction]): basis functions
+
+    Returns:
+        function: function that computes the overlap matrix
+
+    **Example**
+
+    >>> symbols  = ['H', 'H']
+    >>> geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad = False)
+    >>> alpha = np.array([[3.42525091, 0.62391373, 0.1688554],
+    >>>                   [3.42525091, 0.62391373, 0.1688554]], requires_grad=True),
+    >>> mol = Molecule(symbols, geometry, alpha=alpha)
+    >>> args = [alpha]
+    >>> overlap_matrix(mol.basis_set)(*args)
+    array([[1.0, 0.7965883009074122], [0.7965883009074122, 1.0]])
+    """
+
+    def overlap(*args):
+        r"""Construct the overlap matrix for a given set of basis functions.
+
+        Args:
+            args (array[array[float]]): initial values of the differentiable parameters
+
+        Returns:
+            array[array[float]]: the overlap matrix
+        """
+        n = len(basis_functions)
+        overlap_matrix = anp.eye(len(basis_functions))
+        for i, a in enumerate(basis_functions):
+            for j, b in enumerate(basis_functions):
+                if i < j:
+                    if args:
+                        args_ab = []
+                        for arg in args:
+                            args_ab.append(arg[[i, j]])
+                        overlap_integral = generate_overlap(a, b)(*args_ab)
+                    else:
+                        overlap_integral = generate_overlap(a, b)()
+                    o = anp.zeros((n, n))
+                    o[i, j] = o[j, i] = 1.0
+                    overlap_matrix = overlap_matrix + overlap_integral * o
+        return overlap_matrix
+
+    return overlap
+
+
+def generate_kinetic_matrix(basis_functions):
+    r"""Return a function that computes the kinetic matrix for a given set of basis functions.
+
+    Args:
+        basis_functions (list[BasisFunction]): basis functions
+
+    Returns:
+        function: function that computes the kinetic matrix
+
+    **Example**
+
+    >>> symbols  = ['H', 'H']
+    >>> geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad = False)
+    >>> alpha = np.array([[3.42525091, 0.62391373, 0.1688554],
+    >>>                   [3.42525091, 0.62391373, 0.1688554]], requires_grad=True),
+    >>> mol = Molecule(symbols, geometry, alpha=alpha)
+    >>> args = [alpha]
+    >>> kinetic_matrix(mol.basis_set)(*args)
+    array([[0.76003189, 0.38325367], [0.38325367, 0.76003189]])
+    """
+
+    def kinetic(*args):
+        r"""Construct the kinetic matrix for a given set of basis functions.
+
+        Args:
+            args (array[array[float]]): initial values of the differentiable parameters
+
+        Returns:
+            array[array[float]]: the kinetic matrix
+        """
+        n = len(basis_functions)
+        kinetic_matrix = anp.zeros((n, n))
+        for i, a in enumerate(basis_functions):
+            for j, b in enumerate(basis_functions):
+                if i <= j:
+                    if args:
+                        args_ab = []
+                        for arg in args:
+                            args_ab.append(arg[[i, j]])
+                        kinetic_integral = generate_kinetic(a, b)(*args_ab)
+                    else:
+                        kinetic_integral = generate_kinetic(a, b)()
+                    o = anp.zeros((n, n))
+                    o[i, j] = o[j, i] = 1.0
+                    kinetic_matrix = kinetic_matrix + kinetic_integral * o
+        return kinetic_matrix
+
+    return kinetic
+
+
+def generate_attraction_matrix(basis_functions, charges, r):
+    r"""Return a function that computes the electron-nuclear attraction matrix for a given set of
+    basis functions.
+
+    Args:
+        basis_functions (list[BasisFunction]): basis functions
+        charges (list[int]): nuclear charges
+        r (array[float]): nuclear positions
+
+    Returns:
+        function: function that computes the electron-nuclear attraction matrix
+
+    **Example**
+
+    >>> symbols  = ['H', 'H']
+    >>> geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad = False)
+    >>> alpha = np.array([[3.42525091, 0.62391373, 0.1688554],
+    >>>                   [3.42525091, 0.62391373, 0.1688554]], requires_grad=True),
+    >>> mol = Molecule(symbols, geometry, alpha=alpha)
+    >>> args = [alpha]
+    >>> attraction_matrix(mol.basis_set, mol.nuclear_charges, mol.coordinates)(*args)
+    array([[-2.03852057, -1.60241667], [-1.60241667, -2.03852057]])
+    """
+
+    def attraction(*args):
+        r"""Construct the electron-nuclear attraction matrix for a given set of basis functions.
+
+        Args:
+            args (array[array[float]]): initial values of the differentiable parameters
+
+        Returns:
+            array[array[float]]: the electron-nuclear attraction matrix
+        """
+        n = len(basis_functions)
+        attraction_matrix = anp.zeros((n, n))
+        for i, a in enumerate(basis_functions):
+            for j, b in enumerate(basis_functions):
+                attraction_integral = 0
+                if i <= j:
+                    if args:
+                        args_ab = []
+
+                        if r.requires_grad:
+                            for l in range(len(args) - 1):
+                                args_ab.append(args[l + 1][[i, j]])
+                        else:
+                            for arg in args:
+                                args_ab.append(arg[[i, j]])
+
+                        for k, c in enumerate(r):
+                            if c.requires_grad:
+                                args_ab = [args[0][k]] + args_ab
+                            attraction_integral = attraction_integral - charges[
+                                k
+                            ] * generate_attraction(c, a, b)(*args_ab)
+                            if c.requires_grad:
+                                args_ab = args_ab[1:]
+                    else:
+                        for k, c in enumerate(r):
+                            attraction_integral = (
+                                attraction_integral - charges[k] * generate_attraction(c, a, b)()
+                            )
+
+                    o = anp.zeros((n, n))
+                    o[i, j] = o[j, i] = 1.0
+                    attraction_matrix = attraction_matrix + attraction_integral * o
+        return attraction_matrix
+
+    return attraction
+
+
+def generate_repulsion_tensor(basis_functions):
+    r"""Return a function that computes the electron repulsion tensor for a given set of basis
+    functions.
+
+    Args:
+        basis_functions (list[BasisFunction]): basis functions
+
+    Returns:
+        function: function that computes the electron repulsion tensor
+
+    **Example**
+
+    >>> symbols  = ['H', 'H']
+    >>> geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad = False)
+    >>> alpha = np.array([[3.42525091, 0.62391373, 0.1688554],
+    >>>                   [3.42525091, 0.62391373, 0.1688554]], requires_grad=True),
+    >>> mol = Molecule(symbols, geometry, alpha=alpha)
+    >>> args = [alpha]
+    >>> repulsion_tensor(mol.basis_set)(*args)
+    array([[[[0.77460595, 0.56886144], [0.56886144, 0.65017747]],
+            [[0.56886144, 0.45590152], [0.45590152, 0.56886144]]],
+           [[[0.56886144, 0.45590152], [0.45590152, 0.56886144]],
+            [[0.65017747, 0.56886144],[0.56886144, 0.77460595]]]])
+    """
+
+    def repulsion(*args):
+        r"""Construct the electron repulsion tensor for a given set of basis functions.
+
+        Permutational symmetries are taken from [D.F. Brailsford and G.G. Hall, International
+        Journal of Quantum Chemistry, 1971, 5, 657-668].
+
+        Args:
+            args (array[array[float]]): initial values of the differentiable parameters
+
+        Returns:
+            array[array[float]]: the electron repulsion tensor
+        """
+        n = len(basis_functions)
+        repulsion_tensor = anp.zeros((n, n, n, n))
+        e_calc = []
+
+        for i, a in enumerate(basis_functions):
+            for j, b in enumerate(basis_functions):
+                for k, c in enumerate(basis_functions):
+                    for l, d in enumerate(basis_functions):
+
+                        if [i, j, k, l] not in e_calc:
+                            if args:
+                                args_abcd = []
+                                for arg in args:
+                                    args_abcd.append(arg[[i, j, k, l]])
+                                repulsion_integral = generate_repulsion(a, b, c, d)(*args_abcd)
+                            else:
+                                repulsion_integral = generate_repulsion(a, b, c, d)()
+
+                            o = anp.zeros((n, n, n, n))
+                            o[i, j, k, l] = o[k, l, i, j] = o[j, i, l, k] = o[l, k, j, i] = 1.0
+                            o[j, i, k, l] = o[l, k, i, j] = o[i, j, l, k] = o[k, l, j, i] = 1.0
+                            repulsion_tensor = repulsion_tensor + repulsion_integral * o
+                            e_calc = e_calc + [
+                                [i, j, k, l],
+                                [k, l, i, j],
+                                [j, i, l, k],
+                                [l, k, j, i],
+                                [j, i, k, l],
+                                [l, k, i, j],
+                                [i, j, l, k],
+                                [k, l, j, i],
+                            ]
+
+        return repulsion_tensor
+
+    return repulsion
+
+
+def generate_core_matrix(basis_functions, charges, r):
+    r"""Return a function that computes the core matrix for a given set of basis functions.
+
+    The core matrix is computed as a sum of the kinetic and electron-nuclear attraction matrices.
+
+    Args:
+        basis_functions (list[BasisFunction]): basis functions
+        charges (list[int]): nuclear charges
+        r (array[float]): nuclear positions
+
+    Returns:
+        function: function that computes the core matrix
+
+    **Example**
+
+    >>> symbols  = ['H', 'H']
+    >>> geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad = False)
+    >>> alpha = np.array([[3.42525091, 0.62391373, 0.1688554],
+    >>>                   [3.42525091, 0.62391373, 0.1688554]], requires_grad=True),
+    >>> mol = Molecule(symbols, geometry, alpha=alpha)
+    >>> args = [alpha]
+    >>> core_matrix(mol.basis_set, mol.nuclear_charges, mol.coordinates)(*args)
+    array([[-1.27848869, -1.21916299], [-1.21916299, -1.27848869]])
+    """
+
+    def core(*args):
+        r"""Construct the core matrix for a given set of basis functions.
+
+        Args:
+            args (array[array[float]]): initial values of the differentiable parameters
+
+        Returns:
+            array[array[float]]: the core matrix
+        """
+        if r.requires_grad:
+            t = generate_kinetic_matrix(basis_functions)(*args[1:])
+        else:
+            t = generate_kinetic_matrix(basis_functions)(*args)
+
+        a = generate_attraction_matrix(basis_functions, charges, r)(*args)
+        return t + a
+
+    return core
diff --git a/pennylane/hf/molecule.py b/pennylane/hf/molecule.py
new file mode 100644
index 000000000..5db87ea74
--- /dev/null
+++ b/pennylane/hf/molecule.py
@@ -0,0 +1,112 @@
+# Copyright 2018-2021 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""
+This module contains functions and classes to create a Molecule object. This object stores all
+the necessary information to perform a Hartree-Fock calculation for a given molecule.
+"""
+import itertools
+
+# pylint: disable=too-few-public-methods, too-many-arguments, too-many-instance-attributes
+import numpy as np
+from pennylane import numpy as pnp
+from pennylane.hf.basis_data import atomic_numbers
+from pennylane.hf.basis_set import BasisFunction, mol_basis_data
+
+
+class Molecule:
+    r"""Create a molecule object that stores molecular information and default basis set parameters.
+
+    The molecule object can be passed to functions that perform a Hartree-Fock calculation.
+
+    Args:
+        symbols (list[str]): Symbols of the atomic species in the molecule. Currently, atoms with
+            atomic numbers 1-10 are supported.
+        coordinates (array[float]): 1D array with the atomic positions in Cartesian coordinates. The
+            coordinates must be given in atomic units and the size of the array should be ``3*N``
+            where ``N`` is the number of atoms.
+        charge (int): net charge of the molecule
+        mult (int): Spin multiplicity :math:`\mathrm{mult}=N_\mathrm{unpaired} + 1` for
+            :math:`N_\mathrm{unpaired}` unpaired electrons occupying the HF orbitals. Possible
+            values of ``mult`` are :math:`1, 2, 3, \ldots`.
+        basis_name (str): Atomic basis set used to represent the molecular orbitals. Currently, the
+            only supported basis set is 'sto-3g'.
+        l (tuple[int]): angular momentum quantum numbers of the basis function
+        alpha (array[float]): exponents of the primitive Gaussian functions
+        coeff (array[float]): coefficients of the contracted Gaussian functions
+        r (array[float]): positions of the Gaussian functions
+
+    **Example**
+
+    >>> symbols  = ['H', 'H']
+    >>> geometry = pnp.array([[0.0, 0.0, -0.694349],
+    >>>                       [0.0, 0.0,  0.694349]], requires_grad = True)
+    >>> mol = Molecule(symbols, geometry)
+    >>> print(mol.n_electrons)
+    2
+    """
+
+    def __init__(
+        self,
+        symbols,
+        coordinates,
+        charge=0,
+        mult=1,
+        basis_name="sto-3g",
+        l=None,
+        alpha=None,
+        coeff=None,
+    ):
+
+        if basis_name not in ["sto-3g", "STO-3G"]:
+            raise ValueError("Currently, the only supported basis set is 'sto-3g'.")
+
+        if set(symbols) - set(atomic_numbers):
+            raise ValueError(f"Atoms in {set(symbols) - set(atomic_numbers)} are not supported.")
+
+        self.symbols = symbols
+        self.coordinates = coordinates
+        self.charge = charge
+        self.mult = mult
+        self.basis_name = basis_name
+
+        self.n_basis, self.basis_data = mol_basis_data(self.basis_name, self.symbols)
+
+        if l is None:
+            l = [i[0] for i in self.basis_data]
+
+        if alpha is None:
+            alpha = [pnp.array(i[1], requires_grad=False) for i in self.basis_data]
+
+        if coeff is None:
+            coeff = [pnp.array(i[2], requires_grad=False) for i in self.basis_data]
+
+        r = list(
+            itertools.chain(
+                *[[self.coordinates[i]] * self.n_basis[i] for i in range(len(self.n_basis))]
+            )
+        )
+
+        self.l = l
+        self.alpha = alpha
+        self.coeff = coeff
+        self.r = r
+
+        self.basis_set = [
+            BasisFunction(self.l[i], self.alpha[i], self.coeff[i], self.r[i]) for i in range(len(l))
+        ]
+        self.n_orbitals = len(self.l)
+
+        self.nuclear_charges = [atomic_numbers[s] for s in self.symbols]
+
+        self.n_electrons = sum(np.array(self.nuclear_charges)) - self.charge


Test patch:
diff --git a/tests/hf/test_basis_set.py b/tests/hf/test_basis_set.py
new file mode 100644
index 000000000..7b60cb610
--- /dev/null
+++ b/tests/hf/test_basis_set.py
@@ -0,0 +1,112 @@
+# Copyright 2018-2021 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""
+Unit tests for generating basis set default parameters.
+"""
+# pylint: disable=no-self-use
+import pytest
+from pennylane import numpy as np
+from pennylane.hf.basis_set import BasisFunction, atom_basis_data, mol_basis_data
+
+basis_data_H = [
+    (
+        "sto-3g",
+        "H",
+        1,
+        [(0, 0, 0)],
+        [0.3425250914e01, 0.6239137298e00, 0.1688554040e00],
+        [0.1543289673e00, 0.5353281423e00, 0.4446345422e00],
+        [0.0000000000e00, 0.0000000000e00, -0.694349000e00],
+    )
+]
+
+
+class TestBasis:
+    """Tests for generating basis set default parameters"""
+
+    @pytest.mark.parametrize("basis_data", basis_data_H)
+    def test_basisfunction(self, basis_data):
+        """Test that BasisFunction class creates basis function objects correctly."""
+        _, _, _, l, alpha, coeff, r = basis_data
+
+        basis_function = BasisFunction(l, alpha, coeff, r)
+
+        assert np.allclose(np.array(basis_function.alpha), np.array(alpha))
+        assert np.allclose(np.array(basis_function.coeff), np.array(coeff))
+        assert np.allclose(np.array(basis_function.r), np.array(r))
+
+        assert np.allclose(basis_function.params[0], alpha)
+        assert np.allclose(basis_function.params[1], coeff)
+        assert np.allclose(basis_function.params[2], r)
+
+    @pytest.mark.parametrize("basis_data", basis_data_H)
+    def test_atom_basis_data(self, basis_data):
+        """Test that correct basis set parameters are generated for a given atom."""
+        basis, symbol, _, l, alpha, coeff, _ = basis_data
+        params = atom_basis_data(basis, symbol)[0]
+
+        assert np.allclose(params[0], l)
+        assert np.allclose(params[1], alpha)
+        assert np.allclose(params[2], coeff)
+
+    basis_data_HF = [
+        (
+            "sto-3g",
+            ["H", "F"],
+            [1, 5],
+            (
+                (
+                    (0, 0, 0),
+                    [0.3425250914e01, 0.6239137298e00, 0.1688554040e00],
+                    [0.1543289673e00, 0.5353281423e00, 0.4446345422e00],
+                ),
+                (
+                    (0, 0, 0),
+                    [0.1666791340e03, 0.3036081233e02, 0.8216820672e01],
+                    [0.1543289673e00, 0.5353281423e00, 0.4446345422e00],
+                ),
+                (
+                    (0, 0, 0),
+                    [0.6464803249e01, 0.1502281245e01, 0.4885884864e00],
+                    [-0.9996722919e-01, 0.3995128261e00, 0.7001154689e00],
+                ),
+                (
+                    (1, 0, 0),
+                    [0.6464803249e01, 0.1502281245e01, 0.4885884864e00],
+                    [0.1559162750e00, 0.6076837186e00, 0.3919573931e00],
+                ),
+                (
+                    (0, 1, 0),
+                    [0.6464803249e01, 0.1502281245e01, 0.4885884864e00],
+                    [0.1559162750e00, 0.6076837186e00, 0.3919573931e00],
+                ),
+                (
+                    (0, 0, 1),
+                    [0.6464803249e01, 0.1502281245e01, 0.4885884864e00],
+                    [0.1559162750e00, 0.6076837186e00, 0.3919573931e00],
+                ),
+            ),
+        )
+    ]
+
+    @pytest.mark.parametrize("basis_data", basis_data_HF)
+    def test_mol_basis_data(self, basis_data):
+        """Test that correct basis set parameters are generated for a given molecule represented as
+        a list of atoms."""
+        basis, symbols, n_ref, params_ref = basis_data
+        n_basis, params = mol_basis_data(basis, symbols)
+
+        assert n_basis == n_ref
+
+        assert np.allclose(params, params_ref)
diff --git a/tests/hf/test_hamiltonians.py b/tests/hf/test_hamiltonians.py
new file mode 100644
index 000000000..a0ece4bda
--- /dev/null
+++ b/tests/hf/test_hamiltonians.py
@@ -0,0 +1,381 @@
+# Copyright 2018-2021 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""
+Unit tests for functions needed for computing the Hamiltonian.
+"""
+import autograd
+import pennylane as qml
+import pytest
+from pennylane import Identity, PauliX, PauliY, PauliZ
+from pennylane import numpy as np
+from pennylane.hf.hamiltonian import (
+    _generate_qubit_operator,
+    _pauli_mult,
+    _return_pauli,
+    generate_electron_integrals,
+    generate_fermionic_hamiltonian,
+    generate_hamiltonian,
+)
+from pennylane.hf.molecule import Molecule
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "core", "active", "e_core", "one_ref", "two_ref"),
+    [
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad=False),
+            None,
+            None,
+            np.array([1.0000000000321256]),
+            # computed with OpenFermion using molecule.one_body_integrals
+            np.array([[-1.39021927e00, -1.28555566e-16], [-3.52805508e-16, -2.91653305e-01]]),
+            # computed with OpenFermion using molecule.two_body_integrals
+            np.array(
+                [
+                    [
+                        [[7.14439079e-01, 6.62555256e-17], [2.45552260e-16, 1.70241444e-01]],
+                        [[2.45552260e-16, 1.70241444e-01], [7.01853156e-01, 6.51416091e-16]],
+                    ],
+                    [
+                        [[6.62555256e-17, 7.01853156e-01], [1.70241444e-01, 2.72068603e-16]],
+                        [[1.70241444e-01, 2.72068603e-16], [6.51416091e-16, 7.38836693e-01]],
+                    ],
+                ]
+            ),
+        ),
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad=False),
+            [],
+            [0, 1],
+            np.array([1.0000000000321256]),
+            # computed with OpenFermion using molecule.one_body_integrals
+            np.array([[-1.39021927e00, -1.28555566e-16], [-3.52805508e-16, -2.91653305e-01]]),
+            # computed with OpenFermion using molecule.two_body_integrals
+            np.array(
+                [
+                    [
+                        [[7.14439079e-01, 6.62555256e-17], [2.45552260e-16, 1.70241444e-01]],
+                        [[2.45552260e-16, 1.70241444e-01], [7.01853156e-01, 6.51416091e-16]],
+                    ],
+                    [
+                        [[6.62555256e-17, 7.01853156e-01], [1.70241444e-01, 2.72068603e-16]],
+                        [[1.70241444e-01, 2.72068603e-16], [6.51416091e-16, 7.38836693e-01]],
+                    ],
+                ]
+            ),
+        ),
+        (
+            ["Li", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad=False),
+            [0, 1, 2, 3],
+            [4, 5],
+            # reference values of e_core, one and two are computed with our initial prototype code
+            np.array([-5.141222763432437]),
+            np.array([[1.17563204e00, -5.75186616e-18], [-5.75186616e-18, 1.78830226e00]]),
+            np.array(
+                [
+                    [
+                        [[3.12945511e-01, 4.79898448e-19], [4.79898448e-19, 9.78191587e-03]],
+                        [[4.79898448e-19, 9.78191587e-03], [3.00580620e-01, 4.28570365e-18]],
+                    ],
+                    [
+                        [[4.79898448e-19, 3.00580620e-01], [9.78191587e-03, 4.28570365e-18]],
+                        [[9.78191587e-03, 4.28570365e-18], [4.28570365e-18, 5.10996835e-01]],
+                    ],
+                ]
+            ),
+        ),
+    ],
+)
+def test_generate_electron_integrals(symbols, geometry, core, active, e_core, one_ref, two_ref):
+    r"""Test that generate_electron_integrals returns the correct values."""
+    mol = Molecule(symbols, geometry)
+    args = []
+
+    e, one, two = generate_electron_integrals(mol, core=core, active=active)(*args)
+
+    assert np.allclose(e, e_core)
+    assert np.allclose(one, one_ref)
+    assert np.allclose(two, two_ref)
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "alpha", "coeffs_h_ref", "ops_h_ref"),
+    [
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad=False),
+            np.array(
+                [[3.42525091, 0.62391373, 0.1688554], [3.42525091, 0.62391373, 0.1688554]],
+                requires_grad=True,
+            ),
+            # Hamiltonian coefficients and operators computed with OpenFermion using
+            # openfermion.transforms.get_fermion_operator(molecule.get_molecular_hamiltonian())
+            # The "^" symbols in the operators are removed and "," is added for consistency
+            np.array(
+                [
+                    1.0000000000321256,
+                    -1.3902192706002598,
+                    0.35721953951840535,
+                    0.08512072192002007,
+                    0.35721953951840535,
+                    0.08512072192002007,
+                    0.08512072192002007,
+                    0.35092657803574406,
+                    0.08512072192002007,
+                    0.35092657803574406,
+                    0.35721953951840535,
+                    0.08512072192002007,
+                    -1.3902192706002598,
+                    0.35721953951840535,
+                    0.08512072192002007,
+                    0.08512072192002007,
+                    0.35092657803574406,
+                    0.08512072192002007,
+                    0.35092657803574406,
+                    0.35092657803574495,
+                    0.08512072192002007,
+                    0.35092657803574495,
+                    0.08512072192002007,
+                    -0.2916533049477536,
+                    0.08512072192002007,
+                    0.3694183466586136,
+                    0.08512072192002007,
+                    0.3694183466586136,
+                    0.35092657803574495,
+                    0.08512072192002007,
+                    0.35092657803574495,
+                    0.08512072192002007,
+                    0.08512072192002007,
+                    0.3694183466586136,
+                    -0.2916533049477536,
+                    0.08512072192002007,
+                    0.3694183466586136,
+                ]
+            ),
+            [
+                [],
+                [0, 0],
+                [0, 0, 0, 0],
+                [0, 0, 2, 2],
+                [0, 1, 1, 0],
+                [0, 1, 3, 2],
+                [0, 2, 0, 2],
+                [0, 2, 2, 0],
+                [0, 3, 1, 2],
+                [0, 3, 3, 0],
+                [1, 0, 0, 1],
+                [1, 0, 2, 3],
+                [1, 1],
+                [1, 1, 1, 1],
+                [1, 1, 3, 3],
+                [1, 2, 0, 3],
+                [1, 2, 2, 1],
+                [1, 3, 1, 3],
+                [1, 3, 3, 1],
+                [2, 0, 0, 2],
+                [2, 0, 2, 0],
+                [2, 1, 1, 2],
+                [2, 1, 3, 0],
+                [2, 2],
+                [2, 2, 0, 0],
+                [2, 2, 2, 2],
+                [2, 3, 1, 0],
+                [2, 3, 3, 2],
+                [3, 0, 0, 3],
+                [3, 0, 2, 1],
+                [3, 1, 1, 3],
+                [3, 1, 3, 1],
+                [3, 2, 0, 1],
+                [3, 2, 2, 3],
+                [3, 3],
+                [3, 3, 1, 1],
+                [3, 3, 3, 3],
+            ],
+        )
+    ],
+)
+def test_generate_fermionic_hamiltonian(symbols, geometry, alpha, coeffs_h_ref, ops_h_ref):
+    r"""Test that generate_fermionic_hamiltonian returns the correct Hamiltonian."""
+    mol = Molecule(symbols, geometry, alpha=alpha)
+    args = [alpha]
+    h = generate_fermionic_hamiltonian(mol)(*args)
+
+    assert np.allclose(h[0], coeffs_h_ref)
+    assert h[1] == ops_h_ref
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "h_ref_data"),
+    [
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad=False),
+            # computed with qchem.convert_observable and an OpenFermion Hamiltonian; data reordered
+            # h_mol = molecule.get_molecular_hamiltonian()
+            # h_f = openfermion.transforms.get_fermion_operator(h_mol)
+            # h_q = openfermion.transforms.jordan_wigner(h_f)
+            # h_pl = qchem.convert_observable(h_q, wires=[0, 1, 2, 3], tol=(5e-5))
+            (
+                np.array(
+                    [
+                        0.2981788017,
+                        0.2081336485,
+                        0.2081336485,
+                        0.1786097698,
+                        0.042560361,
+                        -0.042560361,
+                        -0.042560361,
+                        0.042560361,
+                        -0.3472487379,
+                        0.1329029281,
+                        -0.3472487379,
+                        0.175463289,
+                        0.175463289,
+                        0.1329029281,
+                        0.1847091733,
+                    ]
+                ),
+                [
+                    Identity(wires=[0]),
+                    PauliZ(wires=[0]),
+                    PauliZ(wires=[1]),
+                    PauliZ(wires=[0]) @ PauliZ(wires=[1]),
+                    PauliY(wires=[0]) @ PauliX(wires=[1]) @ PauliX(wires=[2]) @ PauliY(wires=[3]),
+                    PauliY(wires=[0]) @ PauliY(wires=[1]) @ PauliX(wires=[2]) @ PauliX(wires=[3]),
+                    PauliX(wires=[0]) @ PauliX(wires=[1]) @ PauliY(wires=[2]) @ PauliY(wires=[3]),
+                    PauliX(wires=[0]) @ PauliY(wires=[1]) @ PauliY(wires=[2]) @ PauliX(wires=[3]),
+                    PauliZ(wires=[2]),
+                    PauliZ(wires=[0]) @ PauliZ(wires=[2]),
+                    PauliZ(wires=[3]),
+                    PauliZ(wires=[0]) @ PauliZ(wires=[3]),
+                    PauliZ(wires=[1]) @ PauliZ(wires=[2]),
+                    PauliZ(wires=[1]) @ PauliZ(wires=[3]),
+                    PauliZ(wires=[2]) @ PauliZ(wires=[3]),
+                ],
+            ),
+        )
+    ],
+)
+def test_generate_hamiltonian(symbols, geometry, h_ref_data):
+    r"""Test that generate_hamiltonian returns the correct Hamiltonian."""
+
+    mol = Molecule(symbols, geometry)
+    args = []
+    h = generate_hamiltonian(mol)(*args)
+    h_ref = qml.Hamiltonian(h_ref_data[0], h_ref_data[1])
+
+    assert np.allclose(h.terms[0], h_ref.terms[0])
+    assert qml.Hamiltonian(np.ones(len(h.terms[0])), h.terms[1]).compare(
+        qml.Hamiltonian(np.ones(len(h_ref.terms[0])), h_ref.terms[1])
+    )
+
+
+@pytest.mark.parametrize(
+    ("f_operator", "q_operator"),
+    [
+        (
+            [0, 0],
+            # obtained with openfermion using: jordan_wigner(FermionOperator('0^ 0', 1))
+            # reformatted the original openfermion output: (0.5+0j) [] + (-0.5+0j) [Z0]
+            ([(0.5 + 0j), (-0.5 + 0j)], [[], [(0, "Z")]]),
+        ),
+    ],
+)
+def test_generate_qubit_operator(f_operator, q_operator):
+    r"""Test that _generate_qubit_operator returns the correct operator."""
+    result = _generate_qubit_operator(f_operator)
+
+    assert result == q_operator
+
+
+@pytest.mark.parametrize(
+    ("p1", "p2", "p_ref"),
+    [
+        (
+            [(0, "X"), (1, "Y")],  # X_0 @ Y_1
+            [(0, "X"), (2, "Y")],  # X_0 @ Y_2
+            ([(2, "Y"), (1, "Y")], 1.0),  # X_0 @ Y_1 @ X_0 @ Y_2
+        ),
+    ],
+)
+def test_pauli_mult(p1, p2, p_ref):
+    r"""Test that _generate_qubit_operator returns the correct operator."""
+    result = _pauli_mult(p1, p2)
+
+    assert result == p_ref
+
+
+@pytest.mark.parametrize(
+    ("symbol", "operator"),
+    [
+        ("X", qml.PauliX),
+        ("Y", qml.PauliY),
+        ("Z", qml.PauliZ),
+    ],
+)
+def test_return_pauli(symbol, operator):
+    r"""Test that_return_pauli returns the correct operator."""
+    p = _return_pauli(symbol)
+    assert p is operator
+
+
+def test_gradient_expvalH():
+    r"""Test that the gradient of expval(H) computed with ``autograd.grad`` is equal to the value
+    obtained with the finite difference method."""
+    symbols = ["H", "H"]
+    geometry = (
+        np.array([[0.0, 0.0, -0.3674625962], [0.0, 0.0, 0.3674625962]], requires_grad=False)
+        / 0.529177210903
+    )
+    alpha = np.array(
+        [[3.42525091, 0.62391373, 0.1688554], [3.42525091, 0.62391373, 0.1688554]],
+        requires_grad=True,
+    )
+
+    mol = Molecule(symbols, geometry, alpha=alpha)
+    args = [alpha]
+    dev = qml.device("default.qubit", wires=4)
+
+    def energy(mol):
+        @qml.qnode(dev)
+        def circuit(*args):
+            qml.PauliX(0)
+            qml.PauliX(1)
+            qml.DoubleExcitation(0.22350048111151138, wires=[0, 1, 2, 3])
+            h_qubit = generate_hamiltonian(mol)(*args)
+            return qml.expval(h_qubit)
+
+        return circuit
+
+    grad_autograd = autograd.grad(energy(mol), argnum=0)(*args)
+
+    alpha_1 = np.array(
+        [[3.42515091, 0.62391373, 0.1688554], [3.42525091, 0.62391373, 0.1688554]],
+        requires_grad=False,
+    )  # alpha[0][0] -= 0.0001
+
+    alpha_2 = np.array(
+        [[3.42535091, 0.62391373, 0.1688554], [3.42525091, 0.62391373, 0.1688554]],
+        requires_grad=False,
+    )  # alpha[0][0] += 0.0001
+
+    e_1 = energy(mol)(*[alpha_1])
+    e_2 = energy(mol)(*[alpha_2])
+
+    grad_finitediff = (e_2 - e_1) / 0.0002
+
+    assert np.allclose(grad_autograd[0][0], grad_finitediff)
diff --git a/tests/hf/test_hartree_fock.py b/tests/hf/test_hartree_fock.py
new file mode 100644
index 000000000..c1eedb924
--- /dev/null
+++ b/tests/hf/test_hartree_fock.py
@@ -0,0 +1,171 @@
+# Copyright 2018-2021 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""
+Unit tests for for Hartree-Fock functions.
+"""
+import autograd
+import pytest
+from pennylane import numpy as np
+from pennylane.hf.hartree_fock import generate_scf, hf_energy, nuclear_energy
+from pennylane.hf.molecule import Molecule
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "v_fock", "coeffs", "fock_matrix", "h_core", "repulsion_tensor"),
+    [
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad=False),
+            np.array([-0.67578019, 0.94181155]),
+            np.array([[-0.52754647, -1.56782303], [-0.52754647, 1.56782303]]),
+            np.array([[-0.51126165, -0.70283714], [-0.70283714, -0.51126165]]),
+            np.array([[-1.27848869, -1.21916299], [-1.21916299, -1.27848869]]),
+            np.array(
+                [
+                    [
+                        [[0.77460595, 0.56886144], [0.56886144, 0.65017747]],
+                        [[0.56886144, 0.45590152], [0.45590152, 0.56886144]],
+                    ],
+                    [
+                        [[0.56886144, 0.45590152], [0.45590152, 0.56886144]],
+                        [[0.65017747, 0.56886144], [0.56886144, 0.77460595]],
+                    ],
+                ]
+            ),
+        )
+    ],
+)
+def test_scf(symbols, geometry, v_fock, coeffs, fock_matrix, h_core, repulsion_tensor):
+    r"""Test that generate_scf returns the correct values."""
+    mol = Molecule(symbols, geometry)
+    v, c, f, h, e = generate_scf(mol)()
+
+    assert np.allclose(v, v_fock)
+    assert np.allclose(c, coeffs)
+    assert np.allclose(f, fock_matrix)
+    assert np.allclose(h, h_core)
+    assert np.allclose(e, repulsion_tensor)
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "charge", "e_ref"),
+    [
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad=False),
+            0,
+            # HF energy computed with pyscf using scf.hf.SCF(mol).kernel(numpy.eye(mol.nao_nr()))
+            np.array([-1.06599931664376]),
+        ),
+        (
+            ["H", "H", "H"],
+            np.array([[0.0, 1.0, 0.0], [0.0, 0.0, 0.0], [1.0, 0.0, 0.0]], requires_grad=False),
+            1,
+            # HF energy computed with pyscf using scf.hf.SCF(mol).kernel(numpy.eye(mol.nao_nr()))
+            np.array([-0.948179228995941]),
+        ),
+        (
+            ["H", "F"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad=False),
+            0,
+            # HF energy computed with pyscf using scf.hf.SCF(mol).kernel(numpy.eye(mol.nao_nr()))
+            np.array([-97.8884541671664]),
+        ),
+    ],
+)
+def test_hf_energy(symbols, geometry, charge, e_ref):
+    r"""Test that hf_energy returns the correct energy."""
+    mol = Molecule(symbols, geometry, charge=charge)
+    e = hf_energy(mol)()
+    assert np.allclose(e, e_ref)
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "g_ref"),
+    [
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad=True),
+            # HF gradient computed with pyscf using rhf.nuc_grad_method().kernel()
+            np.array([[0.0, 0.0, 0.3650435], [0.0, 0.0, -0.3650435]]),
+        ),
+        (
+            ["H", "Li"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 2.0]], requires_grad=True),
+            # HF gradient computed with pyscf using rhf.nuc_grad_method().kernel()
+            np.array([[0.0, 0.0, 0.21034957], [0.0, 0.0, -0.21034957]]),
+        ),
+    ],
+)
+def test_hf_energy_gradient(symbols, geometry, g_ref):
+    r"""Test that the gradient of the Hartree-Fock energy wrt differentiable parameters is
+    correct."""
+    mol = Molecule(symbols, geometry)
+    args = [mol.coordinates]
+    g = autograd.grad(hf_energy(mol))(*args)
+
+    assert np.allclose(g, g_ref)
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "e_ref"),
+    [
+        # e_repulsion = \sum_{ij} (q_i * q_j / r_{ij})
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad=False),
+            np.array([1.0]),
+        ),
+        (
+            ["H", "F"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 2.0]], requires_grad=True),
+            np.array([4.5]),
+        ),
+        (
+            ["H", "O", "H"],
+            np.array([[0.0, 1.0, 0.0], [0.0, 0.0, 0.0], [1.0, 0.0, 0.0]], requires_grad=True),
+            np.array([16.707106781186546]),
+        ),
+    ],
+)
+def test_nuclear_energy(symbols, geometry, e_ref):
+    r"""Test that nuclear_energy returns the correct energy."""
+    mol = Molecule(symbols, geometry)
+    args = [mol.coordinates]
+    e = nuclear_energy(mol.nuclear_charges, mol.coordinates)(*args)
+    assert np.allclose(e, e_ref)
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "g_ref"),
+    [
+        # gradient = d(q_i * q_j / (xi - xj)) / dxi, ...
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad=True),
+            np.array([[0.0, 0.0, 1.0], [0.0, 0.0, -1.0]]),
+        ),
+        (
+            ["H", "F"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 2.0]], requires_grad=True),
+            np.array([[0.0, 0.0, 2.25], [0.0, 0.0, -2.25]]),
+        ),
+    ],
+)
+def test_nuclear_energy_gradient(symbols, geometry, g_ref):
+    r"""Test that nuclear energy gradients are correct."""
+    mol = Molecule(symbols, geometry)
+    args = [mol.coordinates]
+    g = autograd.grad(nuclear_energy(mol.nuclear_charges, mol.coordinates))(*args)
+    assert np.allclose(g, g_ref)
diff --git a/tests/hf/test_integrals.py b/tests/hf/test_integrals.py
new file mode 100644
index 000000000..abd99bcec
--- /dev/null
+++ b/tests/hf/test_integrals.py
@@ -0,0 +1,690 @@
+# Copyright 2018-2021 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""
+Unit tests for functions needed to computing integrals over basis functions.
+"""
+import autograd
+import pytest
+from pennylane import numpy as np
+from pennylane.hf.integrals import (
+    _boys,
+    _diff2,
+    _generate_params,
+    _hermite_coulomb,
+    contracted_norm,
+    expansion,
+    gaussian_kinetic,
+    gaussian_overlap,
+    generate_attraction,
+    generate_kinetic,
+    generate_overlap,
+    generate_repulsion,
+    primitive_norm,
+)
+from pennylane.hf.molecule import Molecule
+
+
+@pytest.mark.parametrize(
+    ("l", "alpha", "n"),
+    [
+        # normalization constant for an s orbital is :math:`(\frac {2 \alpha}{\pi})^{{3/4}}`.
+        ((0, 0, 0), np.array([3.425250914]), np.array([1.79444183])),
+    ],
+)
+def test_gaussian_norm(l, alpha, n):
+    r"""Test that the computed normalization constant of a Gaussian function is correct."""
+    assert np.allclose(primitive_norm(l, alpha), n)
+
+
+@pytest.mark.parametrize(
+    ("l", "alpha", "a", "n"),
+    [
+        # normalization constant for a contracted Gaussian function composed of three normalized
+        # s orbital is :math:`1/3`.
+        (
+            (0, 0, 0),
+            np.array([3.425250914, 3.425250914, 3.425250914]),
+            np.array([1.79444183, 1.79444183, 1.79444183]),
+            np.array([0.33333333]),
+        )
+    ],
+)
+def test_contraction_norm(l, alpha, a, n):
+    r"""Test that the computed normalization constant of a contracted Gaussian function is correct."""
+    assert np.allclose(contracted_norm(l, alpha, a), n)
+
+
+@pytest.mark.parametrize(
+    ("alpha", "coeff", "r"),
+    [
+        (
+            np.array([3.42525091, 0.62391373, 0.1688554], requires_grad=True),
+            np.array([0.15432897, 0.53532814, 0.44463454], requires_grad=True),
+            np.array([0.0, 0.0, 0.0], requires_grad=False),
+        ),
+        (
+            np.array([3.42525091, 0.62391373, 0.1688554], requires_grad=False),
+            np.array([0.15432897, 0.53532814, 0.44463454], requires_grad=False),
+            np.array([0.0, 0.0, 0.0], requires_grad=True),
+        ),
+    ],
+)
+def test_generate_params(alpha, coeff, r):
+    r"""Test that test_generate_params returns correct basis set parameters."""
+    params = [alpha, coeff, r]
+    args = [p for p in [alpha, coeff, r] if p.requires_grad]
+    basis_params = _generate_params(params, args)
+
+    assert np.allclose(basis_params, (alpha, coeff, r))
+
+
+@pytest.mark.parametrize(
+    ("la", "lb", "ra", "rb", "alpha", "beta", "t", "c"),
+    [
+        (
+            0,
+            0,
+            np.array([1.2]),
+            np.array([1.2]),
+            np.array([3.42525091]),
+            np.array([3.42525091]),
+            0,
+            np.array([1.0]),
+        ),
+        (
+            1,
+            0,
+            np.array([0.0]),
+            np.array([0.0]),
+            np.array([3.42525091]),
+            np.array([3.42525091]),
+            0,
+            np.array([0.0]),
+        ),
+        (
+            1,
+            1,
+            np.array([0.0]),
+            np.array([10.0]),
+            np.array([3.42525091]),
+            np.array([3.42525091]),
+            0,
+            np.array([0.0]),
+        ),
+    ],
+)
+def test_expansion(la, lb, ra, rb, alpha, beta, t, c):
+    r"""Test that expansion function returns correct value."""
+    assert np.allclose(expansion(la, lb, ra, rb, alpha, beta, t), c)
+    assert np.allclose(expansion(la, lb, ra, rb, alpha, beta, -1), np.array([0.0]))
+    assert np.allclose(expansion(0, 1, ra, rb, alpha, beta, 2), np.array([0.0]))
+
+
+@pytest.mark.parametrize(
+    ("la", "lb", "ra", "rb", "alpha", "beta", "o"),
+    [
+        # two normalized s orbitals
+        (
+            (0, 0, 0),
+            (0, 0, 0),
+            np.array([0.0, 0.0, 0.0]),
+            np.array([0.0, 0.0, 0.0]),
+            np.array([np.pi / 2]),
+            np.array([np.pi / 2]),
+            np.array([1.0]),
+        ),
+        (
+            (0, 0, 0),
+            (0, 0, 0),
+            np.array([0.0, 0.0, 0.0]),
+            np.array([20.0, 0.0, 0.0]),
+            np.array([3.42525091]),
+            np.array([3.42525091]),
+            np.array([0.0]),
+        ),
+        (
+            (1, 0, 0),
+            (0, 0, 1),
+            np.array([0.0, 0.0, 0.0]),
+            np.array([0.0, 0.0, 0.0]),
+            np.array([6.46480325]),
+            np.array([6.46480325]),
+            np.array([0.0]),
+        ),
+    ],
+)
+def test_gaussian_overlap(la, lb, ra, rb, alpha, beta, o):
+    r"""Test that gaussian overlap function returns a correct value."""
+    assert np.allclose(gaussian_overlap(la, lb, ra, rb, alpha, beta), o)
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "alpha", "coef", "r", "o_ref"),
+    [
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 20.0]], requires_grad=False),
+            np.array(
+                [[3.42525091, 0.62391373, 0.1688554], [3.42525091, 0.62391373, 0.1688554]],
+                requires_grad=False,
+            ),
+            np.array(
+                [[0.15432897, 0.53532814, 0.44463454], [0.15432897, 0.53532814, 0.44463454]],
+                requires_grad=True,
+            ),
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 20.0]], requires_grad=True),
+            np.array([0.0]),
+        ),
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]], requires_grad=False),
+            np.array(
+                [[3.42525091, 0.62391373, 0.1688554], [3.42525091, 0.62391373, 0.1688554]],
+                requires_grad=False,
+            ),
+            np.array(
+                [[0.15432897, 0.53532814, 0.44463454], [0.15432897, 0.53532814, 0.44463454]],
+                requires_grad=False,
+            ),
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]], requires_grad=True),
+            np.array([1.0]),
+        ),
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]], requires_grad=True),
+            np.array(
+                [[3.42525091, 0.62391373, 0.1688554], [3.42525091, 0.62391373, 0.1688554]],
+                requires_grad=True,
+            ),
+            np.array(
+                [[0.15432897, 0.53532814, 0.44463454], [0.15432897, 0.53532814, 0.44463454]],
+                requires_grad=True,
+            ),
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]], requires_grad=True),
+            np.array([1.0]),
+        ),
+    ],
+)
+def test_generate_overlap(symbols, geometry, alpha, coef, r, o_ref):
+    r"""Test that generate_overlap function returns a correct value for the overlap integral."""
+    mol = Molecule(symbols, geometry)
+    basis_a = mol.basis_set[0]
+    basis_b = mol.basis_set[1]
+    args = [p for p in [alpha, coef, r] if p.requires_grad]
+
+    o = generate_overlap(basis_a, basis_b)(*args)
+    assert np.allclose(o, o_ref)
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "alpha", "coeff"),
+    [
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad=False),
+            np.array(
+                [[3.42525091, 0.62391373, 0.1688554], [3.42525091, 0.62391373, 0.1688554]],
+                requires_grad=True,
+            ),
+            np.array(
+                [[0.15432897, 0.53532814, 0.44463454], [0.15432897, 0.53532814, 0.44463454]],
+                requires_grad=True,
+            ),
+        ),
+    ],
+)
+def test_gradient(symbols, geometry, alpha, coeff):
+    r"""Test that the overlap gradient computed with respect to the basis parameters is correct."""
+    mol = Molecule(symbols, geometry, alpha=alpha, coeff=coeff)
+    basis_a = mol.basis_set[0]
+    basis_b = mol.basis_set[1]
+    args = [mol.alpha, mol.coeff]
+
+    g_alpha = autograd.grad(generate_overlap(basis_a, basis_b), argnum=0)(*args)
+    g_coeff = autograd.grad(generate_overlap(basis_a, basis_b), argnum=1)(*args)
+
+    # compute overlap gradients with respect to alpha and coeff using finite diff
+    delta = 0.0001
+    g_ref_alpha = np.zeros(6).reshape(alpha.shape)
+    g_ref_coeff = np.zeros(6).reshape(coeff.shape)
+
+    for i in range(len(alpha)):
+        for j in range(len(alpha[0])):
+
+            alpha_minus = alpha.copy()
+            alpha_plus = alpha.copy()
+            alpha_minus[i][j] = alpha_minus[i][j] - delta
+            alpha_plus[i][j] = alpha_plus[i][j] + delta
+            o_minus = generate_overlap(basis_a, basis_b)(*[alpha_minus, coeff])
+            o_plus = generate_overlap(basis_a, basis_b)(*[alpha_plus, coeff])
+            g_ref_alpha[i][j] = (o_plus - o_minus) / (2 * delta)
+
+            coeff_minus = coeff.copy()
+            coeff_plus = coeff.copy()
+            coeff_minus[i][j] = coeff_minus[i][j] - delta
+            coeff_plus[i][j] = coeff_plus[i][j] + delta
+            o_minus = generate_overlap(basis_a, basis_b)(*[alpha, coeff_minus])
+            o_plus = generate_overlap(basis_a, basis_b)(*[alpha, coeff_plus])
+            g_ref_coeff[i][j] = (o_plus - o_minus) / (2 * delta)
+
+    assert np.allclose(g_alpha, g_ref_alpha)
+    assert np.allclose(g_coeff, g_ref_coeff)
+
+
+@pytest.mark.parametrize(
+    ("i", "j", "ri", "rj", "alpha", "beta", "d"),
+    [
+        # _diff2 must return 0.0 for two Gaussians centered far apart at 0.0 and 20.0
+        (
+            0,
+            1,
+            np.array([0.0]),
+            np.array([20.0]),
+            np.array([3.42525091]),
+            np.array([3.42525091]),
+            np.array([0.0]),
+        ),
+        # computed manually
+        (
+            0,
+            0,
+            np.array([0.0]),
+            np.array([1.0]),
+            np.array([3.42525091]),
+            np.array([3.42525091]),
+            np.array([1.01479665]),
+        ),
+    ],
+)
+def test_diff2(i, j, ri, rj, alpha, beta, d):
+    r"""Test that _diff2 function returns a correct value."""
+    assert np.allclose(_diff2(i, j, ri, rj, alpha, beta), d)
+
+
+@pytest.mark.parametrize(
+    ("la", "lb", "ra", "rb", "alpha", "beta", "t"),
+    [
+        # gaussian_kinetic must return 0.0 for two Gaussians centered far apart
+        (
+            (0, 0, 0),
+            (0, 0, 0),
+            np.array([0.0, 0.0, 0.0]),
+            np.array([20.0, 0.0, 0.0]),
+            np.array([3.42525091]),
+            np.array([3.42525091]),
+            np.array([0.0]),
+        ),
+    ],
+)
+def test_gaussian_kinetic(la, lb, ra, rb, alpha, beta, t):
+    r"""Test that gaussian_kinetic function returns a correct value."""
+    assert np.allclose(gaussian_kinetic(la, lb, ra, rb, alpha, beta), t)
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "alpha", "coeff", "t_ref"),
+    [
+        # generate_kinetic must return 0.0 for two Gaussians centered far apart
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 20.0]], requires_grad=False),
+            np.array(
+                [[3.42525091, 0.62391373, 0.1688554], [3.42525091, 0.62391373, 0.1688554]],
+                requires_grad=False,
+            ),
+            np.array(
+                [[0.15432897, 0.53532814, 0.44463454], [0.15432897, 0.53532814, 0.44463454]],
+                requires_grad=True,
+            ),
+            np.array([0.0]),
+        ),
+        # kinetic integral obtained from pyscf using mol.intor('int1e_kin')
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad=False),
+            np.array(
+                [[3.42525091, 0.62391373, 0.1688554], [3.42525091, 0.62391373, 0.1688554]],
+                requires_grad=False,
+            ),
+            np.array(
+                [[0.15432897, 0.53532814, 0.44463454], [0.15432897, 0.53532814, 0.44463454]],
+                requires_grad=True,
+            ),
+            np.array([0.38325384]),
+        ),
+    ],
+)
+def test_generate_kinetic(symbols, geometry, alpha, coeff, t_ref):
+    r"""Test that generate_kinetic function returns a correct value for the kinetic integral."""
+    mol = Molecule(symbols, geometry, alpha=alpha, coeff=coeff)
+    basis_a = mol.basis_set[0]
+    basis_b = mol.basis_set[1]
+    args = [p for p in [alpha, coeff] if p.requires_grad]
+
+    t = generate_kinetic(basis_a, basis_b)(*args)
+    assert np.allclose(t, t_ref)
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "alpha", "coeff"),
+    [
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad=False),
+            np.array(
+                [[3.42525091, 0.62391373, 0.1688554], [3.42525091, 0.62391373, 0.1688554]],
+                requires_grad=True,
+            ),
+            np.array(
+                [[0.15432897, 0.53532814, 0.44463454], [0.15432897, 0.53532814, 0.44463454]],
+                requires_grad=True,
+            ),
+        ),
+    ],
+)
+def test_gradient_kinetic(symbols, geometry, alpha, coeff):
+    r"""Test that the kinetic gradient computed with respect to the basis parameters is correct."""
+    mol = Molecule(symbols, geometry, alpha=alpha, coeff=coeff)
+    basis_a = mol.basis_set[0]
+    basis_b = mol.basis_set[1]
+    args = [mol.alpha, mol.coeff]
+
+    g_alpha = autograd.grad(generate_kinetic(basis_a, basis_b), argnum=0)(*args)
+    g_coeff = autograd.grad(generate_kinetic(basis_a, basis_b), argnum=1)(*args)
+
+    # compute kinetic gradients with respect to alpha, coeff and r using finite diff
+    delta = 0.0001
+    g_ref_alpha = np.zeros(6).reshape(alpha.shape)
+    g_ref_coeff = np.zeros(6).reshape(coeff.shape)
+
+    for i in range(len(alpha)):
+        for j in range(len(alpha[0])):
+
+            alpha_minus = alpha.copy()
+            alpha_plus = alpha.copy()
+            alpha_minus[i][j] = alpha_minus[i][j] - delta
+            alpha_plus[i][j] = alpha_plus[i][j] + delta
+            t_minus = generate_kinetic(basis_a, basis_b)(*[alpha_minus, coeff])
+            t_plus = generate_kinetic(basis_a, basis_b)(*[alpha_plus, coeff])
+            g_ref_alpha[i][j] = (t_plus - t_minus) / (2 * delta)
+
+            coeff_minus = coeff.copy()
+            coeff_plus = coeff.copy()
+            coeff_minus[i][j] = coeff_minus[i][j] - delta
+            coeff_plus[i][j] = coeff_plus[i][j] + delta
+            t_minus = generate_kinetic(basis_a, basis_b)(*[alpha, coeff_minus])
+            t_plus = generate_kinetic(basis_a, basis_b)(*[alpha, coeff_plus])
+            g_ref_coeff[i][j] = (t_plus - t_minus) / (2 * delta)
+
+    assert np.allclose(g_alpha, g_ref_alpha)
+    assert np.allclose(g_coeff, g_ref_coeff)
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "alpha", "coeff", "a_ref"),
+    [
+        # trivial case: integral should be zero since atoms are located very far apart
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 20.0]], requires_grad=False),
+            np.array(
+                [[3.42525091, 0.62391373, 0.1688554], [3.42525091, 0.62391373, 0.1688554]],
+                requires_grad=True,
+            ),
+            np.array(
+                [[0.15432897, 0.53532814, 0.44463454], [0.15432897, 0.53532814, 0.44463454]],
+                requires_grad=True,
+            ),
+            np.array([0.0]),
+        ),
+        # nuclear attraction integral obtained from pyscf using mol.intor('int1e_nuc')
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad=True),
+            np.array(
+                [[3.42525091, 0.62391373, 0.1688554], [3.42525091, 0.62391373, 0.1688554]],
+                requires_grad=True,
+            ),
+            np.array(
+                [[0.15432897, 0.53532814, 0.44463454], [0.15432897, 0.53532814, 0.44463454]],
+                requires_grad=True,
+            ),
+            np.array([0.80120855]),
+        ),
+    ],
+)
+def test_generate_attraction(symbols, geometry, alpha, coeff, a_ref):
+    r"""Test that generate_attraction function returns a correct value for the kinetic integral."""
+    mol = Molecule(symbols, geometry, alpha=alpha, coeff=coeff)
+    basis_a = mol.basis_set[0]
+    basis_b = mol.basis_set[1]
+    args = [p for p in [alpha, coeff] if p.requires_grad]
+
+    if geometry.requires_grad:
+        args = [geometry[0]] + args + [geometry]
+
+    a = generate_attraction(geometry[0], basis_a, basis_b)(*args)
+    assert np.allclose(a, a_ref)
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "alpha", "coeff"),
+    [
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad=False),
+            np.array(
+                [[3.42525091, 0.62391373, 0.1688554], [3.42525091, 0.62391373, 0.1688554]],
+                requires_grad=True,
+            ),
+            np.array(
+                [[0.15432897, 0.53532814, 0.44463454], [0.15432897, 0.53532814, 0.44463454]],
+                requires_grad=True,
+            ),
+        ),
+    ],
+)
+def test_gradient_attraction(symbols, geometry, alpha, coeff):
+    r"""Test that the attraction gradient computed with respect to the basis parameters is correct."""
+    mol = Molecule(symbols, geometry, alpha=alpha, coeff=coeff)
+    basis_a = mol.basis_set[0]
+    basis_b = mol.basis_set[1]
+    args = [mol.alpha, mol.coeff]
+    r_nuc = geometry[0]
+
+    g_alpha = autograd.grad(generate_attraction(r_nuc, basis_a, basis_b), argnum=0)(*args)
+    g_coeff = autograd.grad(generate_attraction(r_nuc, basis_a, basis_b), argnum=1)(*args)
+
+    # compute attraction gradients with respect to alpha and coeff using finite diff
+    delta = 0.0001
+    g_ref_alpha = np.zeros(6).reshape(alpha.shape)
+    g_ref_coeff = np.zeros(6).reshape(coeff.shape)
+
+    for i in range(len(alpha)):
+        for j in range(len(alpha[0])):
+
+            alpha_minus = alpha.copy()
+            alpha_plus = alpha.copy()
+            alpha_minus[i][j] = alpha_minus[i][j] - delta
+            alpha_plus[i][j] = alpha_plus[i][j] + delta
+            a_minus = generate_attraction(r_nuc, basis_a, basis_b)(*[alpha_minus, coeff])
+            a_plus = generate_attraction(r_nuc, basis_a, basis_b)(*[alpha_plus, coeff])
+            g_ref_alpha[i][j] = (a_plus - a_minus) / (2 * delta)
+
+            coeff_minus = coeff.copy()
+            coeff_plus = coeff.copy()
+            coeff_minus[i][j] = coeff_minus[i][j] - delta
+            coeff_plus[i][j] = coeff_plus[i][j] + delta
+            a_minus = generate_attraction(r_nuc, basis_a, basis_b)(*[alpha, coeff_minus])
+            a_plus = generate_attraction(r_nuc, basis_a, basis_b)(*[alpha, coeff_plus])
+            g_ref_coeff[i][j] = (a_plus - a_minus) / (2 * delta)
+
+    assert np.allclose(g_alpha, g_ref_alpha)
+    assert np.allclose(g_coeff, g_ref_coeff)
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "alpha", "coeff", "e_ref"),
+    [
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 20.0]], requires_grad=False),
+            np.array(
+                [
+                    [3.42525091, 0.62391373, 0.1688554],
+                    [3.42525091, 0.62391373, 0.1688554],
+                    [3.42525091, 0.62391373, 0.1688554],
+                    [3.42525091, 0.62391373, 0.1688554],
+                ],
+                requires_grad=False,
+            ),
+            np.array(
+                [
+                    [0.15432897, 0.53532814, 0.44463454],
+                    [0.15432897, 0.53532814, 0.44463454],
+                    [0.15432897, 0.53532814, 0.44463454],
+                    [0.15432897, 0.53532814, 0.44463454],
+                ],
+                requires_grad=True,
+            ),
+            np.array([0.0]),
+        ),
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad=False),
+            np.array(
+                [
+                    [3.42525091, 0.62391373, 0.1688554],
+                    [3.42525091, 0.62391373, 0.1688554],
+                    [3.42525091, 0.62391373, 0.1688554],
+                    [3.42525091, 0.62391373, 0.1688554],
+                ],
+                requires_grad=False,
+            ),
+            np.array(
+                [
+                    [0.15432897, 0.53532814, 0.44463454],
+                    [0.15432897, 0.53532814, 0.44463454],
+                    [0.15432897, 0.53532814, 0.44463454],
+                    [0.15432897, 0.53532814, 0.44463454],
+                ],
+                requires_grad=True,
+            ),
+            np.array([0.45590169]),
+        ),
+    ],
+)
+def test_generate_repulsion(symbols, geometry, alpha, coeff, e_ref):
+    r"""Test that generate_repulsion function returns a correct value for the repulsion integral."""
+    mol = Molecule(symbols, geometry, alpha=alpha, coeff=coeff)
+    basis_a = mol.basis_set[0]
+    basis_b = mol.basis_set[1]
+    args = [p for p in [alpha, coeff] if p.requires_grad]
+
+    a = generate_repulsion(basis_a, basis_b, basis_a, basis_b)(*args)
+
+    assert np.allclose(a, e_ref)
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "alpha", "coeff"),
+    [
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad=False),
+            np.array(
+                [
+                    [3.42525091, 0.62391373, 0.1688554],
+                    [3.42525091, 0.62391373, 0.1688554],
+                    [3.42525091, 0.62391373, 0.1688554],
+                    [3.42525091, 0.62391373, 0.1688554],
+                ],
+                requires_grad=True,
+            ),
+            np.array(
+                [
+                    [0.15432897, 0.53532814, 0.44463454],
+                    [0.15432897, 0.53532814, 0.44463454],
+                    [0.15432897, 0.53532814, 0.44463454],
+                    [0.15432897, 0.53532814, 0.44463454],
+                ],
+                requires_grad=True,
+            ),
+        ),
+    ],
+)
+def test_gradient_repulsion(symbols, geometry, alpha, coeff):
+    r"""Test that the repulsion gradient computed with respect to the basis parameters is correct."""
+    mol = Molecule(symbols, geometry, alpha=alpha, coeff=coeff)
+    basis_a = mol.basis_set[0]
+    basis_b = mol.basis_set[1]
+    args = [mol.alpha, mol.coeff]
+
+    g_alpha = autograd.grad(generate_repulsion(basis_a, basis_b, basis_a, basis_b), argnum=0)(*args)
+    g_coeff = autograd.grad(generate_repulsion(basis_a, basis_b, basis_a, basis_b), argnum=1)(*args)
+
+    # compute repulsion gradients with respect to alpha and coeff using finite diff
+    delta = 0.0001
+    g_ref_alpha = np.zeros(12).reshape(alpha.shape)
+    g_ref_coeff = np.zeros(12).reshape(coeff.shape)
+
+    for i in range(len(alpha)):
+        for j in range(len(alpha[0])):
+
+            alpha_minus = alpha.copy()
+            alpha_plus = alpha.copy()
+            alpha_minus[i][j] = alpha_minus[i][j] - delta
+            alpha_plus[i][j] = alpha_plus[i][j] + delta
+            e_minus = generate_repulsion(basis_a, basis_b, basis_a, basis_b)(*[alpha_minus, coeff])
+            e_plus = generate_repulsion(basis_a, basis_b, basis_a, basis_b)(*[alpha_plus, coeff])
+            g_ref_alpha[i][j] = (e_plus - e_minus) / (2 * delta)
+
+            coeff_minus = coeff.copy()
+            coeff_plus = coeff.copy()
+            coeff_minus[i][j] = coeff_minus[i][j] - delta
+            coeff_plus[i][j] = coeff_plus[i][j] + delta
+            e_minus = generate_repulsion(basis_a, basis_b, basis_a, basis_b)(*[alpha, coeff_minus])
+            e_plus = generate_repulsion(basis_a, basis_b, basis_a, basis_b)(*[alpha, coeff_plus])
+            g_ref_coeff[i][j] = (e_plus - e_minus) / (2 * delta)
+
+    assert np.allclose(g_alpha, g_ref_alpha)
+    assert np.allclose(g_coeff, g_ref_coeff)
+
+
+@pytest.mark.parametrize(
+    ("n", "t", "f_ref"),
+    [(1.25, 0.00, 0.2857142857142857), (2.75, 1.23, 0.061750771828252976)],
+)
+def test_boys(n, t, f_ref):
+    r"""Test that the Boys function is evaluated correctly."""
+    f = _boys(n, t)
+    assert np.allclose(f, f_ref)
+
+
+@pytest.mark.parametrize(
+    ("t", "u", "v", "n", "p", "dr", "h_ref"),
+    [
+        (0, 0, 0, 0, 6.85050183, np.array([0.0, 0.0, 0.0]), 1.0),
+        (0, 0, 1, 0, 6.85050183, np.array([0.0, 0.0, 0.0]), 0.0),
+        (0, 0, 2, 0, 6.85050183, np.array([0.0, 0.0, 0.0]), -4.56700122),
+        (0, 2, 0, 0, 6.85050183, np.array([0.0, 0.0, 0.0]), -4.56700122),
+        (0, 1, 0, 0, 6.85050183, np.array([0.0, 0.0, 0.0]), 0.0),
+        (2, 1, 0, 0, 6.85050183, np.array([0.0, 0.0, 0.0]), 0.0),
+        (1, 1, 0, 0, 6.85050183, np.array([0.0, 0.0, 0.0]), 0.0),
+    ],
+)
+def test_hermite_coulomb(t, u, v, n, p, dr, h_ref):
+    r"""Test that the _hermite_coulomb function returns a correct value."""
+    h = _hermite_coulomb(t, u, v, n, p, dr)
+    assert np.allclose(h, h_ref)
diff --git a/tests/hf/test_matrices.py b/tests/hf/test_matrices.py
new file mode 100644
index 000000000..357354bda
--- /dev/null
+++ b/tests/hf/test_matrices.py
@@ -0,0 +1,535 @@
+# Copyright 2018-2021 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""
+Unit tests for functions needed for computing matrices.
+"""
+import autograd
+import pytest
+from pennylane import numpy as np
+from pennylane.hf.matrices import (
+    generate_attraction_matrix,
+    generate_core_matrix,
+    generate_kinetic_matrix,
+    generate_overlap_matrix,
+    generate_repulsion_tensor,
+    molecular_density_matrix,
+)
+from pennylane.hf.molecule import Molecule
+
+
+@pytest.mark.parametrize(
+    ("n_electron", "c", "p_ref"),
+    [
+        (
+            2,
+            np.array([[-0.54828771, 1.21848441], [-0.54828771, -1.21848441]]),
+            # all P elements are computed as 0.54828771**2 = 0.3006194129370441
+            np.array([[0.30061941, 0.30061941], [0.30061941, 0.30061941]]),
+        ),
+    ],
+)
+def test_molecular_density_matrix(n_electron, c, p_ref):
+    r"""Test that molecular_density_matrix returns the correct matrix."""
+    p = molecular_density_matrix(n_electron, c)
+    assert np.allclose(p, p_ref)
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "alpha", "s_ref"),
+    [
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad=False),
+            np.array(
+                [[3.42525091, 0.62391373, 0.1688554], [3.42525091, 0.62391373, 0.1688554]],
+                requires_grad=True,
+            ),
+            np.array([[1.0, 0.7965883009074122], [0.7965883009074122, 1.0]]),
+        )
+    ],
+)
+def test_overlap_matrix(symbols, geometry, alpha, s_ref):
+    r"""Test that overlap_matrix returns the correct matrix."""
+    mol = Molecule(symbols, geometry, alpha=alpha)
+    args = [alpha]
+    s = generate_overlap_matrix(mol.basis_set)(*args)
+    assert np.allclose(s, s_ref)
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "s_ref"),
+    [
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad=False),
+            np.array([[1.0, 0.7965883009074122], [0.7965883009074122, 1.0]]),
+        )
+    ],
+)
+def test_overlap_matrix_nodiff(symbols, geometry, s_ref):
+    r"""Test that overlap_matrix returns the correct matrix when no differentiable parameter is
+    used."""
+    mol = Molecule(symbols, geometry)
+    s = generate_overlap_matrix(mol.basis_set)()
+    assert np.allclose(s, s_ref)
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "alpha", "coeff", "g_alpha_ref", "g_coeff_ref"),
+    [
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad=False),
+            np.array(
+                [[3.42525091, 0.62391373, 0.1688554], [3.42525091, 0.62391373, 0.1688554]],
+                requires_grad=True,
+            ),
+            np.array(
+                [[0.15432897, 0.53532814, 0.44463454], [0.15432897, 0.53532814, 0.44463454]],
+                requires_grad=True,
+            ),
+            # Jacobian matrix contains gradient of S11, S12, S21, S22 wrt arg_1, arg_2.
+            np.array(
+                [
+                    [
+                        [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]],
+                        [
+                            [-0.00043783, -0.09917143, -0.11600206],
+                            [-0.00043783, -0.09917143, -0.11600206],
+                        ],
+                    ],
+                    [
+                        [
+                            [-0.00043783, -0.09917143, -0.11600206],
+                            [-0.00043783, -0.09917143, -0.11600206],
+                        ],
+                        [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]],
+                    ],
+                ]
+            ),
+            np.array(
+                [
+                    [
+                        [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]],
+                        [
+                            [-0.15627637, -0.02812029, 0.08809831],
+                            [-0.15627637, -0.02812029, 0.08809831],
+                        ],
+                    ],
+                    [
+                        [
+                            [-0.15627637, -0.02812029, 0.08809831],
+                            [-0.15627637, -0.02812029, 0.08809831],
+                        ],
+                        [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]],
+                    ],
+                ]
+            ),
+        )
+    ],
+)
+def test_gradient_overlap_matrix(symbols, geometry, alpha, coeff, g_alpha_ref, g_coeff_ref):
+    r"""Test that the overlap gradients are correct."""
+    mol = Molecule(symbols, geometry, alpha=alpha, coeff=coeff)
+    args = [mol.alpha, mol.coeff]
+    g_alpha = autograd.jacobian(generate_overlap_matrix(mol.basis_set), argnum=0)(*args)
+    g_coeff = autograd.jacobian(generate_overlap_matrix(mol.basis_set), argnum=1)(*args)
+    assert np.allclose(g_alpha, g_alpha_ref)
+    assert np.allclose(g_coeff, g_coeff_ref)
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "alpha", "t_ref"),
+    [
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad=False),
+            np.array(
+                [[3.42525091, 0.62391373, 0.1688554], [3.42525091, 0.62391373, 0.1688554]],
+                requires_grad=True,
+            ),
+            np.array(
+                [
+                    [0.7600318862777408, 0.38325367405372557],
+                    [0.38325367405372557, 0.7600318862777408],
+                ]
+            ),
+        )
+    ],
+)
+def test_kinetic_matrix(symbols, geometry, alpha, t_ref):
+    r"""Test that kinetic_matrix returns the correct matrix."""
+    mol = Molecule(symbols, geometry, alpha=alpha)
+    args = [alpha]
+    t = generate_kinetic_matrix(mol.basis_set)(*args)
+    assert np.allclose(t, t_ref)
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "t_ref"),
+    [
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad=False),
+            np.array(
+                [
+                    [0.7600318862777408, 0.38325367405372557],
+                    [0.38325367405372557, 0.7600318862777408],
+                ]
+            ),
+        )
+    ],
+)
+def test_kinetic_matrix_nodiff(symbols, geometry, t_ref):
+    r"""Test that kinetic_matrix returns the correct matrix when no differentiable parameter is
+    used."""
+    mol = Molecule(symbols, geometry)
+    t = generate_kinetic_matrix(mol.basis_set)()
+    assert np.allclose(t, t_ref)
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "alpha", "coeff", "g_alpha_ref", "g_coeff_ref"),
+    [
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad=False),
+            np.array(
+                [[3.42525091, 0.62391373, 0.1688554], [3.42525091, 0.62391373, 0.1688554]],
+                requires_grad=True,
+            ),
+            np.array(
+                [[0.15432897, 0.53532814, 0.44463454], [0.15432897, 0.53532814, 0.44463454]],
+                requires_grad=True,
+            ),
+            # Jacobian matrix contains gradient of T11, T12, T21, T22 wrt arg_1, arg_2.
+            np.array(
+                [
+                    [
+                        [[0.03263157, 0.85287851, 0.68779528], [0.0, 0.0, 0.0]],
+                        [
+                            [-0.00502729, 0.08211579, 0.3090185],
+                            [-0.00502729, 0.08211579, 0.3090185],
+                        ],
+                    ],
+                    [
+                        [
+                            [-0.00502729, 0.08211579, 0.3090185],
+                            [-0.00502729, 0.08211579, 0.3090185],
+                        ],
+                        [[0.0, 0.0, 0.0], [0.03263157, 0.85287851, 0.68779528]],
+                    ],
+                ]
+            ),
+            np.array(
+                [
+                    [
+                        [[1.824217, 0.10606991, -0.76087597], [0.0, 0.0, 0.0]],
+                        [
+                            [-0.00846016, 0.08488012, -0.09925695],
+                            [-0.00846016, 0.08488012, -0.09925695],
+                        ],
+                    ],
+                    [
+                        [
+                            [-0.00846016, 0.08488012, -0.09925695],
+                            [-0.00846016, 0.08488012, -0.09925695],
+                        ],
+                        [[0.0, 0.0, 0.0], [1.824217, 0.10606991, -0.76087597]],
+                    ],
+                ]
+            ),
+        )
+    ],
+)
+def test_gradient_kinetic_matrix(symbols, geometry, alpha, coeff, g_alpha_ref, g_coeff_ref):
+    r"""Test that the kinetic gradients are correct."""
+    mol = Molecule(symbols, geometry, alpha=alpha, coeff=coeff)
+    args = [mol.alpha, mol.coeff]
+    g_alpha = autograd.jacobian(generate_kinetic_matrix(mol.basis_set), argnum=0)(*args)
+    g_coeff = autograd.jacobian(generate_kinetic_matrix(mol.basis_set), argnum=1)(*args)
+    assert np.allclose(g_alpha, g_alpha_ref)
+    assert np.allclose(g_coeff, g_coeff_ref)
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "alpha", "v_ref"),
+    [
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad=False),
+            np.array(
+                [[3.42525091, 0.62391373, 0.1688554], [3.42525091, 0.62391373, 0.1688554]],
+                requires_grad=True,
+            ),
+            # attraction matrix obtained from pyscf using mol.intor('int1e_nuc')
+            np.array(
+                [
+                    [-2.03852075, -1.6024171],
+                    [-1.6024171, -2.03852075],
+                ]
+            ),
+        )
+    ],
+)
+def test_attraction_matrix(symbols, geometry, alpha, v_ref):
+    r"""Test that attraction_matrix returns the correct matrix."""
+    mol = Molecule(symbols, geometry, alpha=alpha)
+    args = [mol.alpha]
+    v = generate_attraction_matrix(mol.basis_set, mol.nuclear_charges, mol.coordinates)(*args)
+    assert np.allclose(v, v_ref)
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "alpha", "v_ref"),
+    [
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad=True),
+            np.array(
+                [[3.42525091, 0.62391373, 0.1688554], [3.42525091, 0.62391373, 0.1688554]],
+                requires_grad=True,
+            ),
+            # attraction matrix obtained from pyscf using mol.intor('int1e_nuc')
+            np.array(
+                [
+                    [-2.03852075, -1.6024171],
+                    [-1.6024171, -2.03852075],
+                ]
+            ),
+        )
+    ],
+)
+def test_attraction_matrix_diffR(symbols, geometry, alpha, v_ref):
+    r"""Test that attraction_matrix returns the correct matrix when positions are differentiable."""
+    mol = Molecule(symbols, geometry, alpha=alpha)
+    r_basis = mol.coordinates
+    args = [mol.coordinates, mol.alpha, r_basis]
+    v = generate_attraction_matrix(mol.basis_set, mol.nuclear_charges, mol.coordinates)(*args)
+    assert np.allclose(v, v_ref)
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "v_ref"),
+    [
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad=False),
+            # attraction matrix obtained from pyscf using mol.intor('int1e_nuc')
+            np.array(
+                [
+                    [-2.03852075, -1.6024171],
+                    [-1.6024171, -2.03852075],
+                ]
+            ),
+        )
+    ],
+)
+def test_attraction_matrix_nodiff(symbols, geometry, v_ref):
+    r"""Test that attraction_matrix returns the correct matrix."""
+    mol = Molecule(symbols, geometry)
+    v = generate_attraction_matrix(mol.basis_set, mol.nuclear_charges, mol.coordinates)()
+    assert np.allclose(v, v_ref)
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "alpha", "coeff", "g_r_ref"),
+    [
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad=True),
+            np.array(
+                [[3.42525091, 0.62391373, 0.1688554], [3.42525091, 0.62391373, 0.1688554]],
+                requires_grad=True,
+            ),
+            np.array(
+                [[0.15432897, 0.53532814, 0.44463454], [0.15432897, 0.53532814, 0.44463454]],
+                requires_grad=True,
+            ),
+            np.array(
+                [
+                    [
+                        [[0.0, 0.0, 0.0], [0.0, 0.0, 0.44900112]],
+                        [
+                            [0.0, 0.0, -0.26468668],
+                            [0.0, 0.0, 0.26468668],
+                        ],
+                    ],
+                    [
+                        [
+                            [0.0, 0.0, -0.26468668],
+                            [0.0, 0.0, 0.26468668],
+                        ],
+                        [[0.0, 0.0, -0.44900112], [0.0, 0.0, 0.0]],
+                    ],
+                ]
+            ),
+        )
+    ],
+)
+def test_gradient_attraction_matrix(symbols, geometry, alpha, coeff, g_r_ref):
+    r"""Test that the attraction gradients are correct."""
+    mol = Molecule(symbols, geometry, alpha=alpha, coeff=coeff)
+    r_basis = mol.coordinates
+    args = [mol.coordinates, mol.alpha, mol.coeff, r_basis]
+
+    g_r = autograd.jacobian(
+        generate_attraction_matrix(mol.basis_set, mol.nuclear_charges, mol.coordinates), argnum=0
+    )(*args)
+    assert np.allclose(g_r, g_r_ref)
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "alpha", "e_ref"),
+    [
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad=False),
+            np.array(
+                [[3.42525091, 0.62391373, 0.1688554], [3.42525091, 0.62391373, 0.1688554]],
+                requires_grad=True,
+            ),
+            # electron repulsion tensor obtained from pyscf with mol.intor('int2e')
+            np.array(
+                [
+                    [
+                        [[0.77460594, 0.56886157], [0.56886157, 0.65017755]],
+                        [[0.56886157, 0.45590169], [0.45590169, 0.56886157]],
+                    ],
+                    [
+                        [[0.56886157, 0.45590169], [0.45590169, 0.56886157]],
+                        [[0.65017755, 0.56886157], [0.56886157, 0.77460594]],
+                    ],
+                ]
+            ),
+        )
+    ],
+)
+def test_repulsion_tensor(symbols, geometry, alpha, e_ref):
+    r"""Test that repulsion_tensor returns the correct matrix."""
+    mol = Molecule(symbols, geometry, alpha=alpha)
+    args = [mol.alpha]
+    e = generate_repulsion_tensor(mol.basis_set)(*args)
+    assert np.allclose(e, e_ref)
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "e_ref"),
+    [
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad=False),
+            # electron repulsion tensor obtained from pyscf with mol.intor('int2e')
+            np.array(
+                [
+                    [
+                        [[0.77460594, 0.56886157], [0.56886157, 0.65017755]],
+                        [[0.56886157, 0.45590169], [0.45590169, 0.56886157]],
+                    ],
+                    [
+                        [[0.56886157, 0.45590169], [0.45590169, 0.56886157]],
+                        [[0.65017755, 0.56886157], [0.56886157, 0.77460594]],
+                    ],
+                ]
+            ),
+        )
+    ],
+)
+def test_repulsion_tensor_nodiff(symbols, geometry, e_ref):
+    r"""Test that repulsion_tensor returns the correct matrix when no differentiable parameter is
+    used."""
+    mol = Molecule(symbols, geometry)
+    e = generate_repulsion_tensor(mol.basis_set)()
+    assert np.allclose(e, e_ref)
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "alpha", "c_ref"),
+    [
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad=False),
+            np.array(
+                [[3.42525091, 0.62391373, 0.1688554], [3.42525091, 0.62391373, 0.1688554]],
+                requires_grad=True,
+            ),
+            # core matrix obtained from pyscf using scf.RHF(mol).get_hcore()
+            np.array(
+                [
+                    [-1.27848886, -1.21916326],
+                    [-1.21916326, -1.27848886],
+                ]
+            ),
+        )
+    ],
+)
+def test_core_matrix(symbols, geometry, alpha, c_ref):
+    r"""Test that core_matrix returns the correct matrix."""
+    mol = Molecule(symbols, geometry, alpha=alpha)
+    args = [mol.alpha]
+    c = generate_core_matrix(mol.basis_set, mol.nuclear_charges, mol.coordinates)(*args)
+    assert np.allclose(c, c_ref)
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "c_ref"),
+    [
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad=False),
+            # core matrix obtained from pyscf using scf.RHF(mol).get_hcore()
+            np.array(
+                [
+                    [-1.27848886, -1.21916326],
+                    [-1.21916326, -1.27848886],
+                ]
+            ),
+        )
+    ],
+)
+def test_core_matrix_nodiff(symbols, geometry, c_ref):
+    r"""Test that core_matrix returns the correct matrix when no differentiable parameter is
+    used."""
+    mol = Molecule(symbols, geometry)
+    c = generate_core_matrix(mol.basis_set, mol.nuclear_charges, mol.coordinates)()
+    assert np.allclose(c, c_ref)
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "alpha", "c_ref"),
+    [
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad=True),
+            np.array(
+                [[3.42525091, 0.62391373, 0.1688554], [3.42525091, 0.62391373, 0.1688554]],
+                requires_grad=True,
+            ),
+            # core matrix obtained from pyscf using scf.RHF(mol).get_hcore()
+            np.array(
+                [
+                    [-1.27848886, -1.21916326],
+                    [-1.21916326, -1.27848886],
+                ]
+            ),
+        )
+    ],
+)
+def test_core_matrix_diff_positions(symbols, geometry, alpha, c_ref):
+    r"""Test that core_matrix returns the correct matrix when positions are differentiable."""
+    mol = Molecule(symbols, geometry, alpha=alpha)
+    r_basis = mol.coordinates
+    args = [mol.coordinates, mol.alpha, r_basis]
+    c = generate_core_matrix(mol.basis_set, mol.nuclear_charges, mol.coordinates)(*args)
+    assert np.allclose(c, c_ref)
diff --git a/tests/hf/test_molecule.py b/tests/hf/test_molecule.py
new file mode 100644
index 000000000..4218d96d6
--- /dev/null
+++ b/tests/hf/test_molecule.py
@@ -0,0 +1,179 @@
+# Copyright 2018-2021 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""
+Unit tests for the molecule object.
+"""
+# pylint: disable=no-self-use
+import pytest
+from pennylane import numpy as np
+from pennylane.hf.molecule import Molecule
+from pennylane.hf.basis_set import BasisFunction
+
+
+class TestMolecule:
+    """Tests for generating a molecule object."""
+
+    @pytest.mark.parametrize(
+        ("symbols", "geometry"),
+        [
+            (["H", "F"], np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]])),
+        ],
+    )
+    def test_build_molecule(self, symbols, geometry):
+        r"""Test that the generated molecule object has the correct type."""
+        mol = Molecule(symbols, geometry)
+        assert isinstance(mol, Molecule)
+
+    @pytest.mark.parametrize(
+        ("symbols", "geometry"),
+        [
+            (["H", "F"], np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]])),
+        ],
+    )
+    def test_basis_error(self, symbols, geometry):
+        r"""Test that an error is raised if a wrong basis set name is entered."""
+        with pytest.raises(ValueError, match="Currently, the only supported basis set is"):
+            Molecule(symbols, geometry, basis_name="6-31g")
+
+    @pytest.mark.parametrize(
+        ("symbols", "geometry"),
+        [
+            (["H", "Og"], np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]])),
+        ],
+    )
+    def test_symbol_error(self, symbols, geometry):
+        r"""Test that an error is raised if a wrong/not-supported atomic symbol is entered."""
+        with pytest.raises(ValueError, match="are not supported"):
+            Molecule(symbols, geometry)
+
+    @pytest.mark.parametrize(
+        ("symbols", "geometry", "charge", "mult", "basis_name"),
+        [
+            (["H", "F"], np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]]), 0, 1, "sto-3g"),
+        ],
+    )
+    def test_default_inputs(self, symbols, geometry, charge, mult, basis_name):
+        r"""Test that the molecule object contains correct default molecular input data."""
+        mol = Molecule(symbols, geometry)
+
+        assert mol.symbols == symbols
+        assert np.allclose(mol.coordinates, geometry)
+        assert mol.charge == charge
+        assert mol.mult == mult
+        assert mol.basis_name == basis_name
+
+    @pytest.mark.parametrize(
+        ("symbols", "geometry", "n_electrons", "n_orbitals", "nuclear_charges"),
+        [
+            (
+                ["H", "F"],
+                np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]]),
+                10,
+                6,
+                [1, 9],
+            ),
+        ],
+    )
+    def test_molecule_data(self, symbols, geometry, n_electrons, n_orbitals, nuclear_charges):
+        r"""Test that the molecule object contains correct molecular data."""
+        mol = Molecule(symbols, geometry)
+
+        assert mol.n_electrons == n_electrons
+        assert mol.n_orbitals == n_orbitals
+        assert mol.nuclear_charges == nuclear_charges
+
+    @pytest.mark.parametrize(
+        ("symbols", "geometry", "n_basis", "basis_data"),
+        [
+            (
+                ["H", "H"],
+                np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]]),
+                [1, 1],
+                (
+                    (
+                        (0, 0, 0),
+                        [3.42525091, 0.62391373, 0.1688554],
+                        [0.15432897, 0.53532814, 0.44463454],
+                    ),
+                    (
+                        (0, 0, 0),
+                        [3.42525091, 0.62391373, 0.1688554],
+                        [0.15432897, 0.53532814, 0.44463454],
+                    ),
+                ),
+            ),
+        ],
+    )
+    def test_default_basisdata(self, symbols, geometry, n_basis, basis_data):
+        r"""Test that the molecule object contains correct default basis data for a given molecule."""
+        mol = Molecule(symbols, geometry)
+
+        assert mol.n_basis == n_basis
+        assert np.allclose(mol.basis_data, basis_data)
+
+    @pytest.mark.parametrize(
+        ("symbols", "geometry", "l", "alpha", "coeff", "r"),
+        [
+            (
+                ["H", "H"],
+                np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]]),
+                [(0, 0, 0), (0, 0, 0)],
+                [[3.42525091, 0.62391373, 0.1688554], [3.42525091, 0.62391373, 0.1688554]],
+                [[0.15432897, 0.53532814, 0.44463454], [0.15432897, 0.53532814, 0.44463454]],
+                np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]]),
+            ),
+            (
+                ["H", "F"],
+                np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]]),
+                [(0, 0, 0), (0, 0, 0), (0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1)],
+                [
+                    [0.3425250914e01, 0.6239137298e00, 0.1688554040e00],
+                    [0.1666791340e03, 0.3036081233e02, 0.8216820672e01],
+                    [0.6464803249e01, 0.1502281245e01, 0.4885884864e00],
+                    [0.6464803249e01, 0.1502281245e01, 0.4885884864e00],
+                    [0.6464803249e01, 0.1502281245e01, 0.4885884864e00],
+                    [0.6464803249e01, 0.1502281245e01, 0.4885884864e00],
+                ],
+                [
+                    [0.1543289673e00, 0.5353281423e00, 0.4446345422e00],
+                    [0.1543289673e00, 0.5353281423e00, 0.4446345422e00],
+                    [-0.9996722919e-01, 0.3995128261e00, 0.7001154689e00],
+                    [0.1559162750e00, 0.6076837186e00, 0.3919573931e00],
+                    [0.1559162750e00, 0.6076837186e00, 0.3919573931e00],
+                    [0.1559162750e00, 0.6076837186e00, 0.3919573931e00],
+                ],
+                np.array(
+                    [
+                        [0.0, 0.0, 0.0],
+                        [0.0, 0.0, 1.0],
+                        [0.0, 0.0, 1.0],
+                        [0.0, 0.0, 1.0],
+                        [0.0, 0.0, 1.0],
+                        [0.0, 0.0, 1.0],
+                    ]
+                ),
+            ),
+        ],
+    )
+    def test_basisset(self, symbols, geometry, l, alpha, coeff, r):
+        r"""Test that the molecule object contains the correct basis set and non-default basis data
+        for a given molecule.
+        """
+        mol = Molecule(symbols, geometry)
+
+        assert set(map(type, mol.basis_set)) == {BasisFunction}
+        assert mol.l == l
+        assert np.allclose(mol.alpha, alpha)
+        assert np.allclose(mol.coeff, coeff)
+        assert np.allclose(mol.r, r)