Skip to content

Commit

Permalink
Merge pull request #116 from UCL-CCS/specify-frozen-core
Browse files Browse the repository at this point in the history 
Specify frozen core
  • Loading branch information
@MIWdlB
MIWdlB committed Nov 8, 2023
2 parents 72f4d9c + b3ef938 commit 72c9b1a
Show file tree
Hide file tree
Showing 4 changed files with 252 additions and 117 deletions.
1 change: 1 addition & 0 deletions CHANGELOG.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
- Location of explainer notebooks moved up a level to `docs/notebooks`
- Update of dependencies
- Python dependency updated to `>=3.8, <3.11` to pull in symmer.
- `HamBuilder.build()` now accepts MO indices for active space reduction. This overwrites `n_qubits` if used.
### Removed
- `savefile` argument removed from driver, as it is not used.

Expand Down
7 changes: 0 additions & 7 deletions docs/notebooks/2. Orbital Localization Functions.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -162,13 +162,6 @@
"see https://pubs.acs.org/doi/10.1021/acs.jctc.8b01112 for further info"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
Expand Down
187 changes: 145 additions & 42 deletions nbed/ham_builder.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,6 +43,7 @@ def __init__(
self.scf_method = scf_method
self.constant_e_shift = constant_e_shift
self.transform = transform
self._restricted = isinstance(scf_method, (scf.rhf.RHF, dft.rks.RKS))

@property
def _one_body_integrals(self) -> np.ndarray:
Expand All @@ -64,7 +65,7 @@ def _one_body_integrals(self) -> np.ndarray:
hcore = self.scf_method.get_hcore()

# one body terms
if isinstance(self.scf_method, (scf.uhf.UHF, dft.uks.UKS)):
if not self._restricted:
logger.info("Calculating unrestricted one body intergrals.")
one_body_integrals_alpha = (
c_matrix_active[0].T @ hcore[0] @ c_matrix_active[0]
Expand Down Expand Up @@ -95,7 +96,7 @@ def _two_body_integrals(self) -> np.ndarray:
logger.debug("Calculating two body integrals.")
c_matrix_active = self.scf_method.mo_coeff

if isinstance(self.scf_method, (scf.uhf.UHF, dft.uks.UKS)):
if not self._restricted:
n_orbs_alpha = c_matrix_active[0].shape[1]
n_orbs_beta = c_matrix_active[1].shape[1]

Expand Down Expand Up @@ -148,51 +149,119 @@ def _two_body_integrals(self) -> np.ndarray:
logger.debug(f"{two_body_integrals.shape}")
return two_body_integrals

def _reduce_active_space(
def _reduced_orbitals(
self,
qubit_reduction: int,
one_body_integrals: np.ndarray,
two_body_integrals: np.ndarray,
) -> Tuple[float, np.ndarray, np.ndarray]:
"""Reduce the active space to accommodate a certain number of qubits."""
logger.debug("Reducing the active space.")
) -> Tuple[np.ndarray, np.ndarray]:
"""Find the orbitals which correspond to the active space and core.
Args:
qubit_reduction (int): Number of qubits to reduce by.
Returns:
active_indices (np.ndarray): Indices of active orbitals.
core_indices (np.ndarray): Indices of core orbitals.
"""
logger.debug("Finding active space.")
logger.debug(f"Reducing by {qubit_reduction} qubits.")

if type(qubit_reduction) is not int:
logger.error("Invalid qubit_reduction of type %s.", type(qubit_reduction))
raise HamiltonianBuilderError("qubit_reduction must be an Intger")
if qubit_reduction == 0:
logger.debug("No active space reduction required.")
return 0, one_body_integrals, two_body_integrals
if self._restricted:
return np.array([]), np.where(self.scf_method.mo_occ >= 0)[0]
else:
return (
np.array([]),
np.where(self.scf_method.mo_occ.sum(axis=0) >= 0)[0],
)

# +1 because each MO is 2 qubits for closed shell
orbital_reduction = (qubit_reduction + 1) // 2

occupation = self.scf_method.mo_occ
if not self._restricted:
occupation = occupation.sum(axis=0)
occupied = np.where(occupation > 0)[0]
virtual = np.where(occupation == 0)[0]

# find where the last occupied level is
occupied = np.where(self.scf_method.mo_occ > 0)[0]
unoccupied = np.where(self.scf_method.mo_occ == 0)[0]
logger.debug("Occupied orbitals %s.", occupied)
logger.debug("virtual orbitals %s.", virtual)

occupied_reduction = (
orbital_reduction * len(occupied)
) // self._one_body_integrals.shape[-1]
virtual_reduction = orbital_reduction - occupied_reduction
logger.debug(f"Reducing occupied by {occupied_reduction} spatial orbitals.")
logger.debug(f"Reducing virtual by {virtual_reduction} spatial orbitals.")

# +1 because each MO is 2 qubits for closed shell.
n_orbitals = (qubit_reduction + 1) // 2
logger.debug(f"Reducing to {n_orbitals}.")
# Again +1 because we want to use odd numbers to reduce
# occupied orbitals
occupied_reduction = (n_orbitals + 1) // 2
unoccupied_reduction = qubit_reduction - occupied_reduction
core_indices = np.array([])
removed_virtual = np.array([])

if occupied_reduction > 0:
core_indices = occupied[:occupied_reduction]
occupied = occupied[occupied_reduction:]

# We want the MOs nearest the fermi level
# unoccupied orbitals go from 0->N and occupied from N->M
active_indices = np.append(
occupied[occupied_reduction:], unoccupied[:unoccupied_reduction]
)
self._active_space_indices = active_indices
logger.debug(f"Active indices {self._active_space_indices}.")
if virtual_reduction > 0:
removed_virtual = virtual[-virtual_reduction:]
virtual = virtual[:-virtual_reduction]

active_indices = np.append(occupied, virtual)
logger.debug(f"Core indices {core_indices}.")
logger.debug(f"Active indices {active_indices}.")
logger.debug(f"Removed virtual indices {removed_virtual}.")
return core_indices, active_indices

def _reduce_active_space(
self,
one_body_integrals: np.ndarray,
two_body_integrals: np.ndarray,
core_indices: np.ndarray,
active_indices: np.ndarray,
) -> Tuple[float, np.ndarray, np.ndarray]:
"""Reduce the active space to accommodate a certain number of qubits.
Args:
one_body_integrals (np.ndarray): One-electron integrals in physicist notation.
two_body_integrals (np.ndarray): Two-electron integrals in physicist notation.
core_indices (np.ndarray): Indices of core orbitals.
active_indices (np.ndarray): Indices of active orbitals.
"""
logger.debug("Reducing the active space.")
logger.debug(f"{core_indices=}")
logger.debug(f"{active_indices=}")

occupied_indices = np.where(self.scf_method.mo_occ > 0)[0]
core_indices = np.array(core_indices)
active_indices = np.array(active_indices)

# Determine core constant
core_constant = 0.0
for i in occupied_indices:
if core_indices.ndim != 1:
logger.error("Core indices given as dimension %s array.", core_indices.ndim)
raise HamiltonianBuilderError("Core indices must be 1D array.")
if active_indices.ndim != 1:
logger.error(
"Active indices given as dimension %s array.", active_indices.ndim
)
raise HamiltonianBuilderError("Active indices must be 1D array.")
if set(core_indices).intersection(set(active_indices)) != set():
logger.error("Core and active indices overlap.")
raise HamiltonianBuilderError("Core and active indices must not overlap.")
if len(core_indices) + len(active_indices) > self._one_body_integrals.shape[-1]:
logger.error("Too many indices given.")
raise HamiltonianBuilderError(
"Number of core and active indices must not exceed number of orbitals."
)

for i in core_indices:
core_constant += one_body_integrals[0, i, i]
core_constant += one_body_integrals[1, i, i]

for j in occupied_indices:
for j in core_indices:
core_constant += (
two_body_integrals[0, i, j, j, i]
- two_body_integrals[0, i, j, i, j]
Expand All @@ -204,9 +273,9 @@ def _reduce_active_space(

# Modified one electron integrals
one_body_integrals_new = np.copy(one_body_integrals)
for u in active_indices:
for v in active_indices:
for i in occupied_indices:
for i in core_indices:
for u in active_indices:
for v in active_indices:
one_body_integrals_new[0, u, v] += (
two_body_integrals[0, i, u, v, i]
- two_body_integrals[0, i, u, i, v]
Expand All @@ -229,7 +298,6 @@ def _reduce_active_space(
)
]

# Restrict integral ranges and change M appropriately
logger.debug("Active space reduced.")
logger.debug(f"{one_body_integrals_new.shape}")
logger.debug(f"{two_body_integrals_new.shape}")
Expand Down Expand Up @@ -360,7 +428,11 @@ def _taper(self, qham: QubitOperator) -> QubitOperator:
return taper_off_qubits(qham, stabilizers)

def build(
self, n_qubits: Optional[int] = None, taper: Optional[bool] = False
self,
n_qubits: Optional[int] = None,
taper: Optional[bool] = False,
core_indices: Optional[List[int]] = None,
active_indices: Optional[List[int]] = None,
) -> QubitOperator:
"""Returns second quantized fermionic molecular Hamiltonian.
Expand All @@ -372,40 +444,67 @@ def build(
make this approximation.
Args:
scf_method (StreamObject): A pyscf self-consistent method.
constant_e_shift (float): constant energy term to add to Hamiltonian
active_indices (list): A list of spatial orbital indices indicating which orbitals should be
considered active.
occupied_indices (list): A list of spatial orbital indices indicating which orbitals should be
considered doubly occupied.
n_qubits (int): Either total number of qubits to use (positive value) or
number of qubits to reduce size by (negative value).
taper (bool): Whether to taper the Hamiltonian.
core_indices (List[int]): Indices of core orbitals.
active_indices (List[int]): Indices of active orbitals.
Returns:
molecular_hamiltonian (InteractionOperator): fermionic molecular Hamiltonian
molecular_hamiltonian (QubitOperator): Qubit Hamiltonian for molecular system.
"""
logger.debug("Building for %s qubits.", n_qubits)
qubit_reduction = 0

while True:
if n_qubits == 0:
logger.error("n_qubits input as 0.")
message = "n_qubits input as 0.\n"
+"Positive integers can be used to define total qubits used.\n"
+"Negative integers can be used to define a reduction."
raise HamiltonianBuilderError(message)
elif n_qubits is None:
logger.debug("No qubit reduction requested.")
elif n_qubits < 0:
logger.debug("Interpreting negative n_qubits as reduction.")
qubit_reduction = -1 * n_qubits
n_qubits = (self._one_body_integrals.shape[-1] * 2) + n_qubits

logger.info("Building Hamiltonian for %s qubits.", n_qubits)

indices_not_set = (core_indices is None) or (active_indices is None)

max_cycles = 5
for i in range(1, max_cycles + 1):
one_body_integrals = self._one_body_integrals
two_body_integrals = self._two_body_integrals

if indices_not_set:
logger.debug("No active space indices given.")
core_indices, active_indices = self._reduced_orbitals(qubit_reduction)

(
core_constant,
one_body_integrals,
two_body_integrals,
) = self._reduce_active_space(
qubit_reduction, one_body_integrals, two_body_integrals
one_body_integrals, two_body_integrals, core_indices, active_indices
)

one_body_coefficients, two_body_coefficients = self._spinorb_from_spatial(
one_body_integrals, two_body_integrals
)
logger.debug(f"{one_body_coefficients.shape=}")
logger.debug(f"{two_body_coefficients.shape=}")

logger.debug("Building interaction operator.")
molecular_hamiltonian = InteractionOperator(
(self.constant_e_shift + core_constant),
one_body_coefficients,
0.5 * two_body_coefficients,
)
logger.debug(
f"{count_qubits(molecular_hamiltonian)} qubits in Hamiltonian."
)

qham = self._qubit_transform(self.transform, molecular_hamiltonian)

Expand All @@ -421,9 +520,13 @@ def build(
# Wanted to do a recursive thing to get the correct number
# from tapering but it takes ages.
final_n_qubits = count_qubits(qham)
logger.debug(f"{final_n_qubits} qubits used in cycle {i} Hamiltonian.")
if final_n_qubits <= n_qubits:
logger.debug("Hamiltonian reduced to %s qubits.", final_n_qubits)
return qham
if i == max_cycles:
logger.info("Maximum number of cycles reached.")
return qham

# Check that we have the right number of qubits.
qubit_reduction += final_n_qubits - n_qubits
Loading

0 comments on commit 72c9b1a

Please sign in to comment.