Source code for qiskit_nature.second_q.formats.qcschema_translator

# This code is part of a Qiskit project.
#
# (C) Copyright IBM 2022, 2023.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

"""Translator methods for the QCSchema."""

from __future__ import annotations

import logging
from typing import cast

import numpy as np

from qiskit_nature.units import DistanceUnit
from qiskit_nature.second_q.problems import ElectronicBasis, ElectronicStructureProblem
from qiskit_nature.second_q.hamiltonians import ElectronicEnergy
from qiskit_nature.second_q.operators import ElectronicIntegrals
from qiskit_nature.second_q.operators.symmetric_two_body import (
    S1Integrals,
    S4Integrals,
    S8Integrals,
)
from qiskit_nature.second_q.properties import (
    AngularMomentum,
    ElectronicDipoleMoment,
    Magnetization,
    ParticleNumber,
)
from qiskit_nature.second_q.transformers import BasisTransformer

from .molecule_info import MoleculeInfo
from .qcschema import QCSchema

LOGGER = logging.getLogger(__name__)


[docs]def qcschema_to_problem( qcschema: QCSchema, *, basis: ElectronicBasis = ElectronicBasis.MO, include_dipole: bool = True, ) -> ElectronicStructureProblem: """Builds out an :class:`.ElectronicStructureProblem` from a :class:`.QCSchema` instance. This method centralizes the construction of an :class:`.ElectronicStructureProblem` from a :class:`.QCSchema`. Args: qcschema: the :class:`.QCSchema` object from which to build the problem. basis: the :class:`.ElectronicBasis` of the generated problem. include_dipole: whether or not to include an :class:`.ElectronicDipoleMoment` property in the generated problem (if the data is available). Raises: QiskitNatureError: if either of the required 1- or 2-body electronic integrals are missing. NotImplementedError: if an unsupported :class:`.ElectronicBasis` is requested. Returns: An :class:`.ElectronicStructureProblem` instance. """ norb = qcschema.properties.calcinfo_nmo hamiltonian: ElectronicEnergy = None dipole_x: ElectronicIntegrals | None = None dipole_y: ElectronicIntegrals | None = None dipole_z: ElectronicIntegrals | None = None overlap: np.ndarray | None = None if basis == ElectronicBasis.AO: hamiltonian = _get_ao_hamiltonian(qcschema) if qcschema.wavefunction.scf_overlap is not None: nao = hamiltonian.register_length overlap = _reshape_2(qcschema.wavefunction.scf_overlap, nao, nao) if include_dipole: dipole_x = _get_ao_dipole(qcschema, "x") dipole_y = _get_ao_dipole(qcschema, "y") dipole_z = _get_ao_dipole(qcschema, "z") elif basis == ElectronicBasis.MO: hamiltonian = _get_mo_hamiltonian(qcschema) if qcschema.wavefunction.scf_overlap is not None: try: overlap = get_overlap_ab_from_qcschema(qcschema) except AttributeError: if not hamiltonian.electronic_integrals.beta.is_empty() and not np.allclose( hamiltonian.electronic_integrals.alpha["+-"], hamiltonian.electronic_integrals.beta["+-"], ): LOGGER.warning( "Your MO coefficients appear to differ for the alpha- and beta-spin " "orbitals but you did not provide any orbital overlap data. This may lead " "to faulty expectation values of observables which mix alpha- and beta-spin" " components (for example the AngularMomentum)." ) if include_dipole: dipole_x = _get_mo_dipole(qcschema, "x") dipole_y = _get_mo_dipole(qcschema, "y") dipole_z = _get_mo_dipole(qcschema, "z") else: raise NotImplementedError(f"The basis {basis} is not supported by the translation method.") hamiltonian.nuclear_repulsion_energy = qcschema.properties.nuclear_repulsion_energy natm = len(qcschema.molecule.symbols) geo = qcschema.molecule.geometry molecule = MoleculeInfo( symbols=qcschema.molecule.symbols, coords=[(geo[3 * i], geo[3 * i + 1], geo[3 * i + 2]) for i in range(natm)], multiplicity=qcschema.molecule.molecular_multiplicity or 1, charge=qcschema.molecule.molecular_charge or 0, units=DistanceUnit.BOHR, masses=qcschema.molecule.masses, ) num_particles = (qcschema.properties.calcinfo_nalpha, qcschema.properties.calcinfo_nbeta) problem = ElectronicStructureProblem(hamiltonian) problem.basis = basis problem.molecule = molecule problem.num_particles = num_particles problem.num_spatial_orbitals = norb problem.reference_energy = qcschema.properties.return_energy problem.properties.angular_momentum = AngularMomentum(norb, overlap) problem.properties.magnetization = Magnetization(norb) problem.properties.particle_number = ParticleNumber(norb) if qcschema.wavefunction.scf_eigenvalues_a is not None: problem.orbital_energies = np.asarray(qcschema.wavefunction.scf_eigenvalues_a) if qcschema.wavefunction.scf_eigenvalues_b is not None: problem.orbital_energies_b = np.asarray(qcschema.wavefunction.scf_eigenvalues_b) if include_dipole and dipole_x is not None and dipole_y is not None and dipole_z is not None: dipole = ElectronicDipoleMoment(dipole_x, dipole_y, dipole_z) if qcschema.properties.nuclear_dipole_moment is not None: dipole.nuclear_dipole_moment = qcschema.properties.nuclear_dipole_moment problem.properties.electronic_dipole_moment = dipole return problem
def _reshape_2(arr, dim, dim_2=None): return np.asarray(arr).reshape((dim, dim_2 if dim_2 is not None else dim)) def _reshape_4(arr, dim): npair = dim * (dim + 1) // 2 if len(arr) == npair * (npair + 1) // 2: return S8Integrals(np.asarray(arr)) if len(arr) == npair**2: return S4Integrals(np.asarray(arr).reshape((npair,) * 2)) if len(arr) == dim**4: return S1Integrals(np.asarray(arr).reshape((dim,) * 4)) return arr def _get_ao_hamiltonian(qcschema: QCSchema) -> ElectronicEnergy: nao = int(np.sqrt(len(qcschema.wavefunction.scf_fock_a))) hcore = _reshape_2(qcschema.wavefunction.scf_fock_a, nao) hcore_b = None if qcschema.wavefunction.scf_fock_b is not None: hcore_b = _reshape_2(qcschema.wavefunction.scf_fock_b, nao) eri = _reshape_4(qcschema.wavefunction.scf_eri, nao) hamiltonian = ElectronicEnergy.from_raw_integrals(hcore, eri, hcore_b) return hamiltonian def _get_mo_hamiltonian(qcschema: QCSchema) -> ElectronicEnergy: if qcschema.wavefunction.scf_fock_mo_a is not None: return _get_mo_hamiltonian_direct(qcschema) hamiltonian = _get_ao_hamiltonian(qcschema) transformer = get_ao_to_mo_from_qcschema(qcschema) return cast(ElectronicEnergy, transformer.transform_hamiltonian(hamiltonian)) def _get_mo_hamiltonian_direct(qcschema: QCSchema) -> ElectronicEnergy: norb = qcschema.properties.calcinfo_nmo hij = _reshape_2(qcschema.wavefunction.scf_fock_mo_a, norb) hijkl = _reshape_4(qcschema.wavefunction.scf_eri_mo_aa, norb) hij_b = None hijkl_bb = None hijkl_ba = None if qcschema.wavefunction.scf_fock_mo_b is not None: hij_b = _reshape_2(qcschema.wavefunction.scf_fock_mo_b, norb) if qcschema.wavefunction.scf_eri_mo_bb is not None: hijkl_bb = _reshape_4(qcschema.wavefunction.scf_eri_mo_bb, norb) if qcschema.wavefunction.scf_eri_mo_ba is not None: hijkl_ba = _reshape_4(qcschema.wavefunction.scf_eri_mo_ba, norb) if qcschema.wavefunction.scf_eri_mo_ab is not None and hijkl_ba is None: hijkl_ba = np.transpose(_reshape_4(qcschema.wavefunction.scf_eri_mo_ab, norb)) hamiltonian = ElectronicEnergy.from_raw_integrals(hij, hijkl, hij_b, hijkl_bb, hijkl_ba) return hamiltonian def _get_ao_dipole(qcschema, axis) -> ElectronicIntegrals | None: name = f"scf_dipole_{axis}" integrals = getattr(qcschema.wavefunction, name, None) if integrals is None: return None nao = int(np.sqrt(len(integrals))) ints = _reshape_2(integrals, nao) return ElectronicIntegrals.from_raw_integrals(ints) def _get_mo_dipole(qcschema, axis) -> ElectronicIntegrals | None: name_a = f"scf_dipole_mo_{axis}_a" if getattr(qcschema.wavefunction, name_a, None) is not None: return _get_mo_dipole_direct(qcschema, axis) dipole = _get_ao_dipole(qcschema, axis) transformer = get_ao_to_mo_from_qcschema(qcschema) return transformer.transform_electronic_integrals(dipole) def _get_mo_dipole_direct(qcschema, axis) -> ElectronicIntegrals | None: norb = qcschema.properties.calcinfo_nmo name_a = f"scf_dipole_mo_{axis}_a" integrals_a = getattr(qcschema.wavefunction, name_a, None) if integrals_a is None: return None ints_a = _reshape_2(integrals_a, norb) name_b = f"scf_dipole_mo_{axis}_b" integrals_b = getattr(qcschema.wavefunction, name_b, None) ints_b = None if integrals_b is not None: ints_b = _reshape_2(integrals_b, norb) return ElectronicIntegrals.from_raw_integrals(ints_a, h1_b=ints_b)
[docs]def get_ao_to_mo_from_qcschema(qcschema: QCSchema) -> BasisTransformer: """Builds out a :class:`.BasisTransformer` from a :class:`.QCSchema` instance. This utility extracts the AO-to-MO conversion coefficients, contained in a :class:`.QCSchema` object. Args: qcschema: the :class:`.QCSchema` object from which to build the problem. Returns: A :class:`.BasisTransformer` instance. """ nmo = qcschema.properties.calcinfo_nmo nao = len(qcschema.wavefunction.scf_orbitals_a) // nmo coeff_a = _reshape_2(qcschema.wavefunction.scf_orbitals_a, nao, nmo) coeff_b = None if qcschema.wavefunction.scf_orbitals_b is not None: coeff_b = _reshape_2(qcschema.wavefunction.scf_orbitals_b, nao, nmo) transformer = BasisTransformer( ElectronicBasis.AO, ElectronicBasis.MO, ElectronicIntegrals.from_raw_integrals(coeff_a, h1_b=coeff_b), ) return transformer
[docs]def get_overlap_ab_from_qcschema(qcschema: QCSchema) -> np.ndarray | None: """Builds the alpha-beta spin orbital overlap matrix from a :class:`.QCSchema` instance. Args: qcschema: the :class:`.QCSchema` object from which to build the problem. Raises: AttributeError: when the overlap or MO coefficients are missing. Returns: The overlap matrix as a 2-dimensional numpy array or `None` if the overlap is equal to the identity matrix. """ if qcschema.wavefunction.scf_orbitals_a is None or qcschema.wavefunction.scf_overlap is None: raise AttributeError nmo = qcschema.properties.calcinfo_nmo nao = len(qcschema.wavefunction.scf_orbitals_a) // nmo coeff_a = _reshape_2(qcschema.wavefunction.scf_orbitals_a, nao, nmo) coeff_b = None if qcschema.wavefunction.scf_orbitals_b is not None: coeff_b = _reshape_2(qcschema.wavefunction.scf_orbitals_b, nao, nmo) if coeff_b is None: # when coeff_b is None, that means it is equal to coeff_a, which in turn means that the # overlap is equal to the identity return None overlap = _reshape_2(qcschema.wavefunction.scf_overlap, nao, nao) return coeff_a.T @ overlap @ coeff_b