Source code for pennylane.qchem.hamiltonian

# 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.
"""
from functools import singledispatch

# pylint: disable= too-many-branches, too-many-arguments, too-many-locals, too-many-nested-blocks
# pylint: disable=consider-using-generator, protected-access
import pennylane as qml
from pennylane.operation import active_new_opmath
from pennylane.pauli.utils import simplify

from .basis_data import atomic_numbers
from .hartree_fock import nuclear_energy, scf
from .molecule import Molecule
from .observable_hf import fermionic_observable, qubit_observable

# Bohr-Angstrom correlation coefficient (https://physics.nist.gov/cgi-bin/cuu/Value?bohrrada0)
bohr_angs = 0.529177210903


[docs]def 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 (~qchem.molecule.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 = qml.qchem.Molecule(symbols, geometry, alpha=alpha) >>> args = [alpha] >>> 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 = scf(mol)(*args) one = qml.math.einsum("qr,rs,st->qt", coeffs.T, h_core, coeffs) two = qml.math.swapaxes( qml.math.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 = qml.math.zeros(one.shape) o[p, q] = 1.0 one = one + (2 * two[i][p][q][i] - two[i][p][i][q]) * o one = one[qml.math.ix_(active, active)] two = two[qml.math.ix_(active, active, active, active)] return core_constant, one, two return _electron_integrals
[docs]def fermionic_hamiltonian(mol, cutoff=1.0e-12, core=None, active=None): r"""Return a function that computes the fermionic Hamiltonian. Args: mol (~qchem.molecule.Molecule): the molecule object cutoff (float): cutoff value for discarding the negligible electronic integrals core (list[int]): indices of the core orbitals active (list[int]): indices of the active orbitals 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 = qml.qchem.Molecule(symbols, geometry, alpha=alpha) >>> args = [alpha] >>> h = 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: FermiSentence: fermionic Hamiltonian """ core_constant, one, two = electron_integrals(mol, core, active)(*args) return fermionic_observable(core_constant, one, two, cutoff) return _fermionic_hamiltonian
[docs]def diff_hamiltonian(mol, cutoff=1.0e-12, core=None, active=None, mapping="jordan_wigner"): r"""Return a function that computes the qubit Hamiltonian. Args: mol (~qchem.molecule.Molecule): the molecule object cutoff (float): cutoff value for discarding the negligible electronic integrals core (list[int]): indices of the core orbitals active (list[int]): indices of the active orbitals mapping (str): Specifies the fermion-to-qubit mapping. Input values can be ``'jordan_wigner'``, ``'parity'`` or ``'bravyi_kitaev'``. Returns: function: function that computes the qubit hamiltonian **Example** >>> from pennylane import numpy as pnp >>> symbols = ['H', 'H'] >>> geometry = pnp.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], requires_grad = False) >>> alpha = pnp.array([[3.42525091, 0.62391373, 0.1688554], >>> [3.42525091, 0.62391373, 0.1688554]], requires_grad=True) >>> mol = qml.qchem.Molecule(symbols, geometry, alpha=alpha) >>> args = [alpha] >>> h = qml.qchem.diff_hamiltonian(mol)(*args) >>> h.terms()[0] [tensor(0.29817878, requires_grad=True), tensor(0.20813366, requires_grad=True), tensor(-0.34724872, requires_grad=True), tensor(0.13290292, requires_grad=True), tensor(0.20813366, requires_grad=True), tensor(0.17860977, requires_grad=True), tensor(0.04256036, requires_grad=True), tensor(-0.04256036, requires_grad=True), tensor(-0.04256036, requires_grad=True), tensor(0.04256036, requires_grad=True), tensor(-0.34724872, requires_grad=True), tensor(0.17546328, requires_grad=True), tensor(0.13290292, requires_grad=True), tensor(0.17546328, requires_grad=True), tensor(0.18470917, requires_grad=True)] """ def _molecular_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 = fermionic_hamiltonian(mol, cutoff, core, active)(*args) return qubit_observable(h_ferm, mapping=mapping) return _molecular_hamiltonian
[docs]def molecular_hamiltonian(*args, **kwargs): """molecular_hamiltonian(molecule, method="dhf", active_electrons=None, active_orbitals=None,\ mapping="jordan_wigner", outpath=".", wires=None, args=None, convert_tol=1e12) Generate the qubit Hamiltonian of a molecule. This function drives the construction of the second-quantized electronic Hamiltonian of a molecule and its transformation to the basis of Pauli matrices. The net charge of the molecule can be given to simulate cationic/anionic systems. Also, the spin multiplicity can be input to determine the number of unpaired electrons occupying the HF orbitals as illustrated in the left panel of the figure below. The basis of Gaussian-type *atomic* orbitals used to represent the *molecular* orbitals can be specified to go beyond the minimum basis approximation. An active space can be defined for a given number of *active electrons* occupying a reduced set of *active orbitals* as sketched in the right panel of the figure below. | .. figure:: ../../_static/qchem/fig_mult_active_space.png :align: center :width: 90% | Args: molecule (~qchem.molecule.Molecule): the molecule object method (str): Quantum chemistry method used to solve the mean field electronic structure problem. Available options are ``method="dhf"`` to specify the built-in differentiable Hartree-Fock solver, ``method="pyscf"`` to use the PySCF package (requires ``pyscf`` to be installed), or ``method="openfermion"`` to use the OpenFermion-PySCF plugin (this requires ``openfermionpyscf`` to be installed). active_electrons (int): Number of active electrons. If not specified, all electrons are considered to be active. active_orbitals (int): Number of active orbitals. If not specified, all orbitals are considered to be active. mapping (str): transformation used to map the fermionic Hamiltonian to the qubit Hamiltonian outpath (str): path to the directory containing output files wires (Wires, list, tuple, dict): Custom wire mapping for connecting to Pennylane ansatz. For types ``Wires``/``list``/``tuple``, each item in the iterable represents a wire label corresponding to the qubit number equal to its index. For type dict, only int-keyed dict (for qubit-to-wire conversion) is accepted for partial mapping. If None, will use identity map. args (array[array[float]]): initial values of the differentiable parameters convert_tol (float): Tolerance in `machine epsilon <https://numpy.org/doc/stable/reference/generated/numpy.real_if_close.html>`_ for the imaginary part of the Hamiltonian coefficients created by openfermion. Coefficients with imaginary part less than 2.22e-16*tol are considered to be real. Returns: tuple[pennylane.Hamiltonian, int]: the fermionic-to-qubit transformed Hamiltonian and the number of qubits .. note:: The ``molecular_hamiltonian`` function accepts a ``Molecule`` object as its first argument. Look at the `Usage Details` for more details on the old interface. **Example** >>> symbols = ['H', 'H'] >>> coordinates = np.array([[0., 0., -0.66140414], [0., 0., 0.66140414]]) >>> molecule = qml.qchem.Molecule(symbols, coordinates) >>> H, qubits = qml.qchem.molecular_hamiltonian(molecule) >>> print(qubits) 4 >>> print(H) (-0.04207897647782188) [I0] + (0.17771287465139934) [Z0] + (0.1777128746513993) [Z1] + (-0.24274280513140484) [Z2] + (-0.24274280513140484) [Z3] + (0.17059738328801055) [Z0 Z1] + (0.04475014401535161) [Y0 X1 X2 Y3] + (-0.04475014401535161) [Y0 Y1 X2 X3] + (-0.04475014401535161) [X0 X1 Y2 Y3] + (0.04475014401535161) [X0 Y1 Y2 X3] + (0.12293305056183801) [Z0 Z2] + (0.1676831945771896) [Z0 Z3] + (0.1676831945771896) [Z1 Z2] + (0.12293305056183801) [Z1 Z3] + (0.176276408043196) [Z2 Z3] .. details:: :title: Usage Details The old interface for this method involved passing molecular information as separate arguments: ``molecular_hamiltonian``\\ (`symbols, coordinates, name='molecule', charge=0, mult=1, basis='sto-3g',` `method='dhf', active_electrons=None, active_orbitals=None, mapping='jordan_wigner', outpath='.',` `wires=None, alpha=None, coeff=None, args=None, load_data=False, convert_tol=1e12`) Molecule-based Arguments: - **symbols** (list[str]): symbols of the atomic species in the molecule - **coordinates** (array[float]): atomic positions in Cartesian coordinates. The atomic coordinates must be in atomic units and can be given as either a 1D array of size ``3*N``, or a 2D array of shape ``(N, 3)`` where ``N`` is the number of atoms. name (str): name of the molecule - **charge** (int): Net charge of the molecule. If not specified a neutral system is assumed. - **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`. If not specified, a closed-shell HF state is assumed. - **basis** (str): atomic basis set used to represent the molecular orbitals - **alpha** (array[float]): exponents of the primitive Gaussian functions - **coeff** (array[float]): coefficients of the contracted Gaussian functions Therefore, a molecular Hamiltonian had to be constructed in the following manner: .. code-block:: python from pennylane import qchem symbols = ["H", "H"] geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]]) H, qubit = qchem.molecular_hamiltonian(symbols, geometry, charge=0) As part of the new interface, we are shifting towards extracting all the molecular information from the :class:`~.qchem.molecule.Molecule` within the ``molecular_hamiltonian`` method. """ if len(args) != 0: return _molecular_hamiltonian_dispatch(*args, **kwargs) method = kwargs.pop("symbols", None) or kwargs.pop("molecule", None) if method is not None: return _molecular_hamiltonian_dispatch(method, **kwargs) raise NotImplementedError( "The provided arguments do not contain information about symbols in the molecule. " "Please provide that information in the form of a molecule object or as a list of symbols." )
@singledispatch def _molecular_hamiltonian_dispatch(*args, **kwargs): r"""Generate the qubit Hamiltonian of a molecule.""" raise NotImplementedError( "molecular_hamiltonian supports only list or molecule object types. " "Please provide one of them." ) @_molecular_hamiltonian_dispatch.register(Molecule) def _( molecule, method="dhf", active_electrons=None, active_orbitals=None, mapping="jordan_wigner", outpath=".", wires=None, args=None, convert_tol=1e12, ): return _molecular_hamiltonian( molecule.symbols, molecule.coordinates, molecule.name, molecule.charge, molecule.mult, molecule.basis_name, method, active_electrons, active_orbitals, mapping, outpath, wires, molecule.alpha, molecule.coeff, args, molecule.load_data, convert_tol, ) @_molecular_hamiltonian_dispatch.register(list) def _( symbols, coordinates, unit="bohr", name="molecule", charge=0, mult=1, basis="sto-3g", method="dhf", active_electrons=None, active_orbitals=None, mapping="jordan_wigner", outpath=".", wires=None, alpha=None, coeff=None, args=None, load_data=False, convert_tol=1e12, ): if (coord_unit := unit.strip().lower()) not in ("angstrom", "bohr"): raise ValueError( f"The provided unit '{unit}' is not supported. " f"Please set 'unit' to 'bohr' or 'angstrom'." ) if coord_unit == "angstrom": coordinates = coordinates / bohr_angs return _molecular_hamiltonian( symbols, coordinates=coordinates, name=name, charge=charge, mult=mult, basis=basis, method=method, active_electrons=active_electrons, active_orbitals=active_orbitals, mapping=mapping, outpath=outpath, wires=wires, alpha=alpha, coeff=coeff, args=args, load_data=load_data, convert_tol=convert_tol, ) def _molecular_hamiltonian( symbols=None, coordinates=None, name="molecule", charge=0, mult=1, basis="sto-3g", method="dhf", active_electrons=None, active_orbitals=None, mapping="jordan_wigner", outpath=".", wires=None, alpha=None, coeff=None, args=None, load_data=False, convert_tol=1e12, ): # pylint:disable=too-many-arguments, too-many-statements r"""Generate the qubit Hamiltonian of a molecule.""" if method not in ["dhf", "pyscf", "openfermion"]: raise ValueError("Only 'dhf', 'pyscf' and 'openfermion' backends are supported.") if mapping.strip().lower() not in ["jordan_wigner", "parity", "bravyi_kitaev"]: raise ValueError( f"'{mapping}' is not supported." f"Please set the mapping to 'jordan_wigner', 'parity' or 'bravyi_kitaev'." ) if len(coordinates) == len(symbols) * 3: geometry_dhf = qml.numpy.array(coordinates.reshape(len(symbols), 3)) geometry_hf = coordinates elif len(coordinates) == len(symbols): geometry_dhf = qml.numpy.array(coordinates) geometry_hf = coordinates.flatten() wires_map = None if wires: wires_new = qml.qchem.convert._process_wires(wires) wires_map = dict(zip(range(len(wires_new)), list(wires_new.labels))) if method in ("dhf", "pyscf"): n_electrons = sum([atomic_numbers[s] for s in symbols]) - charge if n_electrons % 2 == 1 or mult != 1: raise ValueError( "Open-shell systems are not supported for the requested backend. Use " "method = 'openfermion' or change the charge or spin multiplicity of the molecule." ) if method == "dhf": if args is None and isinstance(geometry_dhf, qml.numpy.tensor): geometry_dhf.requires_grad = False mol = qml.qchem.Molecule( symbols, geometry_dhf, charge=charge, mult=mult, basis_name=basis, load_data=load_data, alpha=alpha, coeff=coeff, ) core, active = qml.qchem.active_space( mol.n_electrons, mol.n_orbitals, mult, active_electrons, active_orbitals ) requires_grad = args is not None h = ( qml.qchem.diff_hamiltonian(mol, core=core, active=active, mapping=mapping)(*args) if requires_grad else qml.qchem.diff_hamiltonian(mol, core=core, active=active, mapping=mapping)() ) if active_new_opmath(): h_as_ps = qml.pauli.pauli_sentence(h) coeffs = qml.numpy.real(list(h_as_ps.values()), requires_grad=requires_grad) h_as_ps = qml.pauli.PauliSentence(dict(zip(h_as_ps.keys(), coeffs))) h = ( qml.s_prod(0, qml.Identity(h.wires[0])) if len(h_as_ps) == 0 else h_as_ps.operation() ) else: coeffs = qml.numpy.real(h.coeffs, requires_grad=requires_grad) h = qml.Hamiltonian(coeffs, h.ops) if wires: h = qml.map_wires(h, wires_map) return h, 2 * len(active) if method == "pyscf": core_constant, one_mo, two_mo = qml.qchem.openfermion_pyscf._pyscf_integrals( symbols, geometry_hf, charge, mult, basis, active_electrons, active_orbitals ) hf = qml.qchem.fermionic_observable(core_constant, one_mo, two_mo) mapping = mapping.strip().lower() qubits = len(hf.wires) if active_new_opmath(): if mapping == "jordan_wigner": h_pl = qml.jordan_wigner(hf, wire_map=wires_map, tol=1.0e-10) elif mapping == "parity": h_pl = qml.parity_transform(hf, qubits, wire_map=wires_map, tol=1.0e-10) elif mapping == "bravyi_kitaev": h_pl = qml.bravyi_kitaev(hf, qubits, wire_map=wires_map, tol=1.0e-10) h_pl.simplify() else: if mapping == "jordan_wigner": h_pl = qml.jordan_wigner(hf, ps=True, wire_map=wires_map, tol=1.0e-10) elif mapping == "parity": h_pl = qml.parity_transform(hf, qubits, ps=True, wire_map=wires_map, tol=1.0e-10) elif mapping == "bravyi_kitaev": h_pl = qml.bravyi_kitaev(hf, qubits, ps=True, wire_map=wires_map, tol=1.0e-10) h_pl = simplify(h_pl.hamiltonian()) return h_pl, len(h_pl.wires) h_pl = qml.qchem.openfermion_pyscf._openfermion_hamiltonian( symbols, geometry_hf, name, charge, mult, basis, active_electrons, active_orbitals, mapping, outpath, wires, convert_tol, ) return h_pl, len(h_pl.wires)