Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Openshell DMET #208

Merged
Merged
Show file tree
Hide file tree
Changes from 26 commits
Commits
Show all changes
29 commits
Select commit Hold shift + click to select a range
45514a1
Open-shell DMET.
alexfleury-sb Apr 21, 2022
4778089
Fix for get_rdm CCSD.
alexfleury-sb Apr 21, 2022
63849e2
Merge remote-tracking branch 'upstream/develop' into openshell_dmet.
alexfleury-sb Sep 9, 2022
08f7bfd
Added NAO localization.
alexfleury-sb Sep 9, 2022
21ef023
Added LiO2 spin=1 DMET test.
alexfleury-sb Sep 9, 2022
8bf770a
Merge remote-tracking branch 'upstream/develop' into openshell_dmet.
alexfleury-sb Nov 29, 2022
4f14881
Get_rdm based on develop.
alexfleury-sb Dec 6, 2022
32d3a92
Reinstated the fixes for CCSD RDMS.
alexfleury-sb Dec 6, 2022
fc2ca02
Added UHF MF for DMET. New get_rdm for VQESolver.
alexfleury-sb Dec 14, 2022
e0f2112
Fix syntax error.
alexfleury-sb Dec 15, 2022
6fe0f33
Fixing tests.
alexfleury-sb Dec 15, 2022
5373000
Conformance tests.
alexfleury-sb Dec 15, 2022
1a867da
Docs.
alexfleury-sb Dec 16, 2022
f15c66b
Fix doc type.
alexfleury-sb Dec 16, 2022
41d16d7
Merge remote-tracking branch 'upstream/develop' into openshell_dmet.
alexfleury-sb Dec 22, 2022
8c8a30c
Fixes for first review.
alexfleury-sb Dec 22, 2022
5e4c10a
Variable name in verbose print.
alexfleury-sb Dec 22, 2022
a2b7f38
Merge remote-tracking branch 'upstream/develop' into openshell_dmet.
alexfleury-sb Dec 23, 2022
f68018b
Second round of fixes.
alexfleury-sb Dec 23, 2022
b020781
Removed reduce, dict in function for alpha/beta, vqe changes.
alexfleury-sb Dec 23, 2022
1a26d17
Working sate with \.
alexfleury-sb Dec 23, 2022
ddc6656
Fix failing test because of removed param.
alexfleury-sb Dec 24, 2022
3706d7f
No numpy for diagonal substraction, doc and variable names.
alexfleury-sb Jan 3, 2023
8b1bc53
Conformance test for indent + too long.
alexfleury-sb Jan 3, 2023
cf1eb80
Removed int.
alexfleury-sb Jan 3, 2023
9729dde
Typos in active.
alexfleury-sb Jan 3, 2023
4fa47a6
List to tuple.
alexfleury-sb Jan 3, 2023
6633b99
Breaking lines and doc.
alexfleury-sb Jan 3, 2023
6b008cf
Type hint.
alexfleury-sb Jan 3, 2023
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
34 changes: 23 additions & 11 deletions tangelo/algorithms/classical/ccsd_solver.py
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@
"""Class performing electronic structure calculation employing the CCSD method.
"""

import numpy as np
from pyscf import cc, lib
from pyscf.cc.ccsd_rdm import _make_rdm1, _make_rdm2, _gamma1_intermediates, _gamma2_outcore
from pyscf.cc.uccsd_rdm import (_make_rdm1 as _umake_rdm1, _make_rdm2 as _umake_rdm2,
Expand All @@ -39,6 +40,8 @@ class CCSDSolver(ElectronicStructureSolver):
def __init__(self, molecule):
self.cc_fragment = None

self.spin = molecule.spin

self.mean_field = molecule.mean_field
self.frozen = molecule.frozen_mos
self.uhf = molecule.uhf
Expand Down Expand Up @@ -77,18 +80,27 @@ def get_rdm(self):
raise RuntimeError("CCSDSolver: Cannot retrieve RDM. Please run the 'simulate' method first")

# Solve the lambda equation and obtain the reduced density matrix from CC calculation
self.cc_fragment.solve_lambda()
t1 = self.cc_fragment.t1
t2 = self.cc_fragment.t2
l1 = self.cc_fragment.l1
l2 = self.cc_fragment.l2

d1 = _gamma1_intermediates(self.cc_fragment, t1, t2, l1, l2) if not self.uhf else _ugamma1_intermediates(self.cc_fragment, t1, t2, l1, l2)
f = lib.H5TmpFile()
d2 = _gamma2_outcore(self.cc_fragment, t1, t2, l1, l2, f, False) if not self.uhf else _ugamma2_outcore(self.cc_fragment, t1, t2, l1, l2, f, False)

one_rdm = _make_rdm1(self.cc_fragment, d1, with_frozen=False) if not self.uhf else _umake_rdm1(self.cc_fragment, d1, with_frozen=False)
two_rdm = (_make_rdm2(self.cc_fragment, d1, d2, with_dm1=True, with_frozen=False) if not self.uhf
else _umake_rdm2(self.cc_fragment, d1, d2, with_dm1=True, with_frozen=False))
l1, l2 = self.cc_fragment.solve_lambda(t1, t2)

if self.spin == 0 and not self.uhf:
d1 = _gamma1_intermediates(self.cc_fragment, t1, t2, l1, l2)
f = lib.H5TmpFile()
d2 = _gamma2_outcore(self.cc_fragment, t1, t2, l1, l2, f, False)

one_rdm = _make_rdm1(self.cc_fragment, d1, with_frozen=False)
two_rdm = _make_rdm2(self.cc_fragment, d1, d2, with_dm1=True, with_frozen=False)
else:
d1 = _ugamma1_intermediates(self.cc_fragment, t1, t2, l1, l2)
f = lib.H5TmpFile()
d2 = _ugamma2_outcore(self.cc_fragment, t1, t2, l1, l2, f, False)

one_rdm = _umake_rdm1(self.cc_fragment, d1, with_frozen=False)
two_rdm = _umake_rdm2(self.cc_fragment, d1, d2, with_dm1=True, with_frozen=False)

if not self.uhf:
one_rdm = np.sum(one_rdm, axis=0)
two_rdm = np.sum((two_rdm[0], 2*two_rdm[1], two_rdm[2]), axis=0)

return one_rdm, two_rdm
130 changes: 130 additions & 0 deletions tangelo/algorithms/variational/vqe_solver.py
Original file line number Diff line number Diff line change
Expand Up @@ -506,6 +506,136 @@ def get_rdm(self, var_params, resample=False, sum_spin=True, ref_state=Circuit()

return rdm1_spin, rdm2_spin

def get_rdm_uhf(self, var_params, resample=False, ref_state=Circuit()):
"""Compute the 1- and 2- RDM matrices using the VQE energy evaluation.
This method allows to combine the DMET problem decomposition technique
with the VQE as an electronic structure solver. The RDMs are computed by
using each fermionic Hamiltonian term, transforming them and computing
the elements one-by-one. Note that the Hamiltonian coefficients will not
be multiplied as in the energy evaluation. The first element of the
Hamiltonian is the nuclear repulsion energy term, not the Hamiltonian
term.

Args:
var_params (numpy.array or list): variational parameters to use for
rdm calculation
resample (bool): Whether to resample saved frequencies. get_rdm with
savefrequencies=True must be called or a dictionary for each
qubit terms' frequencies must be set to self.rdm_freq_dict
ref_state (Circuit): A reference state preparation circuit.

Returns: TODO
(numpy.array, numpy.array): One & two-particle spin summed RDMs if
sumspin=True or the full One & two-Particle RDMs if
sumspin=False.
"""

self.ansatz.update_var_params(var_params)

# Initialize the RDM arrays
n_mol_orbitals = max(self.molecule.n_active_mos)
rdm1_np_a = np.zeros((n_mol_orbitals,) * 2)
rdm1_np_b = np.zeros((n_mol_orbitals,) * 2)
rdm2_np_a = np.zeros((n_mol_orbitals,) * 4)
rdm2_np_b = np.zeros((n_mol_orbitals,) * 4)
rdm2_np_ba = np.zeros((n_mol_orbitals,) * 4)

# If resampling is requested, check that a previous savefrequencies run has been called
if resample:
if hasattr(self, "rdm_freq_dict"):
qb_freq_dict = self.rdm_freq_dict
resampled_expect_dict = dict()
else:
raise AttributeError("Need to run RDM calculation with savefrequencies=True")
else:
qb_freq_dict = dict()
qb_expect_dict = dict()

# Loop over each element of Hamiltonian (non-zero value)
for key in self.molecule.fermionic_hamiltonian.terms:
# Ignore constant / empty term
if not key:
continue

# Assign indices depending on one- or two-body term
length = len(key)
# One-body terms.
if(length == 2):
alexfleury-sb marked this conversation as resolved.
Show resolved Hide resolved
pele, qele = int(key[0][0]), int(key[1][0])
iele, jele = pele // 2, qele // 2
iele_r, jele_r = pele % 2, qele % 2
# Two-body terms.
elif(length == 4):
pele, qele, rele, sele = int(key[0][0]), int(key[1][0]), int(key[2][0]), int(key[3][0])
iele, jele, kele, lele = pele // 2, qele // 2, rele // 2, sele // 2
iele_r, jele_r, kele_r, lele_r = pele % 2, qele % 2, rele % 2, sele % 2

# Create the Hamiltonian with the correct key (Set coefficient to one)
hamiltonian_temp = FermionOperator(key)

# Obtain qubit Hamiltonian
qubit_hamiltonian2 = fermion_to_qubit_mapping(fermion_operator=hamiltonian_temp,
mapping=self.qubit_mapping,
n_spinorbitals=self.molecule.n_active_sos,
n_electrons=self.molecule.n_active_electrons,
up_then_down=self.up_then_down,
spin=self.molecule.spin)
qubit_hamiltonian2.compress()

# Run through each qubit term separately, use previously calculated result for the qubit term or
# calculate and save results for that qubit term
opt_energy2 = 0.
for qb_term, qb_coef in qubit_hamiltonian2.terms.items():
if qb_term:
if qb_term not in qb_freq_dict:
if resample:
warnings.warn(f"Warning: rerunning circuit for missing qubit term {qb_term}")
basis_circuit = Circuit(measurement_basis_gates(qb_term))
full_circuit = ref_state + self.ansatz.circuit + basis_circuit
qb_freq_dict[qb_term], _ = self.backend.simulate(full_circuit)
if resample:
if qb_term not in resampled_expect_dict:
resampled_freq_dict = get_resampled_frequencies(qb_freq_dict[qb_term], self.backend.n_shots)
resampled_expect_dict[qb_term] = self.backend.get_expectation_value_from_frequencies_oneterm(qb_term, resampled_freq_dict)
expectation = resampled_expect_dict[qb_term]
else:
if qb_term not in qb_expect_dict:
qb_expect_dict[qb_term] = self.backend.get_expectation_value_from_frequencies_oneterm(qb_term, qb_freq_dict[qb_term])
expectation = qb_expect_dict[qb_term]
opt_energy2 += qb_coef * expectation
else:
opt_energy2 += qb_coef

# Put the values in np arrays (differentiate 1- and 2-RDM)
if length == 2:
alexfleury-sb marked this conversation as resolved.
Show resolved Hide resolved
if (iele_r, jele_r) == (0, 0):
rdm1_np_a[iele, jele] += opt_energy2
elif (iele_r, jele_r) == (1, 1):
rdm1_np_b[iele, jele] += opt_energy2
elif length == 4:
if ((iele != lele) or (jele != kele)):
if (iele_r, jele_r, kele_r, lele_r) == (0, 0, 0, 0):
rdm2_np_a[lele, iele, kele, jele] += 0.5 * opt_energy2
rdm2_np_a[iele, lele, jele, kele] += 0.5 * opt_energy2
elif (iele_r, jele_r, kele_r, lele_r) == (1, 1, 1, 1):
rdm2_np_b[lele, iele, kele, jele] += 0.5 * opt_energy2
rdm2_np_b[iele, lele, jele, kele] += 0.5 * opt_energy2
elif (iele_r, jele_r, kele_r, lele_r) == (0, 1, 1, 0):
rdm2_np_ba[iele, lele, jele, kele] += 0.5 * opt_energy2
rdm2_np_ba[lele, iele, kele, jele] += 0.5 * opt_energy2
else:
if (iele_r, jele_r, kele_r, lele_r) == (0, 0, 0, 0):
rdm2_np_a[iele, lele, jele, kele] += opt_energy2
elif (iele_r, jele_r, kele_r, lele_r) == (1, 1, 1, 1):
rdm2_np_b[iele, lele, jele, kele] += opt_energy2
elif (iele_r, jele_r, kele_r, lele_r) == (0, 1, 1, 0):
rdm2_np_ba[iele, lele, jele, kele] += opt_energy2

# save rdm frequency dictionary
self.rdm_freq_dict = qb_freq_dict

return (rdm1_np_a, rdm1_np_b), (rdm2_np_a, rdm2_np_ba, rdm2_np_b)

def _default_optimizer(self, func, var_params):
"""Function used as a default optimizer for VQE when user does not
provide one. Can be used as an example for users who wish to provide
Expand Down
12 changes: 9 additions & 3 deletions tangelo/problem_decomposition/dmet/_helpers/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -12,10 +12,16 @@
# See the License for the specific language governing permissions and
# limitations under the License.

# Common helper functions for all mean-field.
from .dmet_orbitals import dmet_orbitals as _orbitals
from .dmet_fragment import dmet_fragment_constructor as _fragment_constructor
from .dmet_onerdm import dmet_low_rdm as _low_rdm
from .dmet_onerdm import dmet_fragment_rdm as _fragment_rdm
from .dmet_bath import dmet_fragment_bath as _fragment_bath
from .dmet_scf_guess import dmet_fragment_guess as _fragment_guess
from .dmet_scf import dmet_fragment_scf as _fragment_scf

# Specific helper functions for restricted / unrestricted mean-field.
from .dmet_onerdm import dmet_low_rdm_rhf as _low_rdm_rhf
from .dmet_onerdm import dmet_low_rdm_rohf_uhf as _low_rdm_rohf_uhf
from .dmet_scf_guess import dmet_fragment_guess_rhf as _fragment_guess_rhf
from .dmet_scf_guess import dmet_fragment_guess_rohf_uhf as _fragment_guess_rohf_uhf
from .dmet_scf import dmet_fragment_scf_rhf as _fragment_scf_rhf
from .dmet_scf import dmet_fragment_scf_rohf_uhf as _fragment_scf_rohf_uhf
15 changes: 7 additions & 8 deletions tangelo/problem_decomposition/dmet/_helpers/dmet_fragment.py
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@
"""


def dmet_fragment_constructor(mol, atom_list, number_fragment):
def dmet_fragment_constructor(mol, atom_list, n_fragment):
"""Construct orbital list.

Make a list of number of orbitals for each fragment while obtaining the list
Expand All @@ -28,8 +28,7 @@ def dmet_fragment_constructor(mol, atom_list, number_fragment):
Args:
mol (pyscf.gto.Mole): The molecule to simulate (The full molecule).
atom_list (list): The atom list for each fragment (int).
number_fragment (list): Number of element in atom list per fragment
(int).
n_fragment (list): Number of atoms per fragment (int).

Returns:
list: The number of orbitals for each fragment (int).
Expand All @@ -39,18 +38,18 @@ def dmet_fragment_constructor(mol, atom_list, number_fragment):
"""

# Make a new atom list based on how many fragments for DMET calculation
if number_fragment == 0:
if n_fragment == 0:
atom_list2 = atom_list
else:
# Calculate the number of DMET calculations
number_new_fragment = int(len(atom_list)/(number_fragment+1)) # number of DMET calulation per loop
n_new_fragment = int(len(atom_list)/(n_fragment+1))
atom_list2 = []

# Define the number of atoms per DMET calculation
for i in range(number_new_fragment):
for i in range(n_new_fragment):
num = 0
for j in range(number_fragment + 1):
k = (number_fragment+1)*i+j
for j in range(n_fragment + 1):
k = (n_fragment + 1) * i + j
num += atom_list[k]
atom_list2.append(num)

Expand Down
47 changes: 37 additions & 10 deletions tangelo/problem_decomposition/dmet/_helpers/dmet_onerdm.py
Original file line number Diff line number Diff line change
Expand Up @@ -18,41 +18,68 @@
"""

import numpy as np
from functools import reduce


def dmet_low_rdm(active_fock, number_active_electrons):
def dmet_low_rdm_rhf(active_fock, n_active_electrons):
alexfleury-sb marked this conversation as resolved.
Show resolved Hide resolved
"""Construct the one-particle RDM from low-level calculation.

Args:
active_fock (numpy.array): Fock matrix from low-level calculation
(float64).
number_active_electrons (int): Number of electrons in the entire system.
n_active_electrons (int): Number of electrons in the entire system.

Returns:
numpy.array: One-particle RDM of the low-level calculation (float64).
"""

# Extract the occupied part of the one-particle RDM
num_occ = number_active_electrons / 2
num_occ = n_active_electrons // 2
e, c = np.linalg.eigh(active_fock)
new_index = e.argsort()
e = e[new_index]
c = c[:, new_index]
onerdm = np.dot(c[:, : int(num_occ)], c[:, : int(num_occ)].T) * 2
onerdm = np.dot(c[:, : num_occ], c[:, : num_occ].T) * 2

return onerdm


def dmet_fragment_rdm(t_list, bath_orb, core_occupied, number_active_electrons):
def dmet_low_rdm_rohf_uhf(active_fock_alpha, active_fock_beta, n_active_alpha, n_active_beta):
"""Construct the one-particle RDM from low-level calculation.

Args:
active_fock_alpha (numpy.array): Fock matrix from low-level calculation
(float64).
active_fock_beta (numpy.array): Fock matrix from low-level calculation
(float64).
n_active_alpha (int): Number of alpha electrons.
n_active_beta (int): Number of beta electrons.

Returns:
onerdm (numpy.array): One-particle RDM of the low-level calculation
(float64).
"""

e, c = np.linalg.eigh(active_fock_alpha)
new_index = e.argsort()
c = c[:, new_index]
onerdm_alpha = np.dot(c[:, :int(n_active_alpha)], c[:, :int(n_active_alpha)].T)

e, c = np.linalg.eigh(active_fock_beta)
new_index = e.argsort()
c = c[:, new_index]
onerdm_beta = np.dot(c[:, :int(n_active_beta)], c[:, :int(n_active_beta)].T)

return onerdm_alpha + onerdm_beta


def dmet_fragment_rdm(t_list, bath_orb, core_occupied, n_active_electrons):
"""Construct the one-particle RDM for the core orbitals.

Args:
t_list (list): Number of [0] fragment & [1] bath orbitals (int).
bath_orb (numpy.array): The bath orbitals (float64).
core_occupied (numpy.array): Core occupied part of the MO coefficients
(float64).
number_active_electrons (int): Number of electrons in the entire system.
n_active_electrons (int): Number of electrons in the entire system.

Returns:
int: Number of orbitals for fragment calculation.
Expand All @@ -72,9 +99,9 @@ def dmet_fragment_rdm(t_list, bath_orb, core_occupied, number_active_electrons):

# Define the number of electrons in the fragment
number_ele_temp = np.sum(core_occupied)
number_electrons = int(round(number_active_electrons - number_ele_temp))
number_electrons = int(round(n_active_electrons - number_ele_temp))

# Obtain the one particle RDM for the fragment (core)
core_occupied_onerdm = reduce(np.dot, (bath_orb, np.diag(core_occupied), bath_orb.T))
core_occupied_onerdm = bath_orb @ np.diag(core_occupied) @ bath_orb.T

return number_orbitals, number_electrons, core_occupied_onerdm
Loading