From 218e792a37f3040cb336b2b793d5a89937baa901 Mon Sep 17 00:00:00 2001 From: Sean Reiter Date: Fri, 18 Aug 2023 10:59:49 -0400 Subject: [PATCH] Batching of determinants (#69) * Implementation naive algo * pt2 selection now iterates over the full internal space and only a chunk of the connected. selection_step is now an mpi function * fixed doctest * for real this time * for real this time * some initial changes before adding parameter to chunk connected space by input size * spelling * doctest fix * Added functionality for user to specify the size of the working portion of the connected space * Local/global sorts seperate in the selection step + added few more tests * Changes from recent review * Initial PR for chunking the external space (#53) * pt2 selection now iterates over the full internal space and only a chunk of the connected. selection_step is now an mpi function * fixed doctest * for real this time * for real this time * some initial changes before adding parameter to chunk connected space by input size * spelling * doctest fix * Added functionality for user to specify the size of the working portion of the connected space * Local/global sorts seperate in the selection step + added few more tests * Changes from recent review * Constraint based triplet excitations (#7) * Add constraint-based excitation functions (no testing yet) * small name change * now constrained excitation handles opposite spin ones as well; this is a starting point for tests * no way this is working first try? * first test passing * Minimal test added and passing * Constraint-based generation of connected dets from wf + new test * Now constraints working for beta electrons as well * black * black * changes * build fix attempt * attempt 2 * Update test_python.yml * update requirements.txt * Update test_python.yml * remove enumerate * Update test_python.yml * nathan PR (#9) * External chunk into clean (#61) * Implementation naive algo * pt2 selection now iterates over the full internal space and only a chunk of the connected. selection_step is now an mpi function * fixed doctest * for real this time * for real this time * some initial changes before adding parameter to chunk connected space by input size * spelling * doctest fix * Added functionality for user to specify the size of the working portion of the connected space * Local/global sorts seperate in the selection step + added few more tests * Changes from recent review * Initial PR for chunking the external space (#53) * pt2 selection now iterates over the full internal space and only a chunk of the connected. selection_step is now an mpi function * fixed doctest * for real this time * for real this time * some initial changes before adding parameter to chunk connected space by input size * spelling * doctest fix * Added functionality for user to specify the size of the working portion of the connected space * Local/global sorts seperate in the selection step + added few more tests * Changes from recent review * Constraint based triplet excitations (#7) * Add constraint-based excitation functions (no testing yet) * small name change * now constrained excitation handles opposite spin ones as well; this is a starting point for tests * no way this is working first try? * first test passing * Minimal test added and passing * Constraint-based generation of connected dets from wf + new test * Now constraints working for beta electrons as well * black * black * changes * build fix attempt * attempt 2 * Update test_python.yml * update requirements.txt * Update test_python.yml * remove enumerate * Update test_python.yml --------- Co-authored-by: Thomas Applencourt * only master loads wave function and integrals (#62) * only master loads wave function and integrals * added # because it was missing --------- Co-authored-by: Thomas Applencourt Co-authored-by: Nathan W Johnson <115745236+nathanjohnsongithub@users.noreply.github.com> * New pt2 (#11) * remove useless decorator * added two electron parent class * Not all tests will pass (intdriven not implemented yet); detdriven on the fly pt2 is done * remove old tests for chunked connected space (no longer necessary since all pt2 done on the fly) * forgot to remove something * selection step modified to work with on the fly pt2 * remove leftover print statements * very naive integral driven on the fly pt2 * load balancing * Static load balancing, currently very slow. Optimize later * update comment * Integral-driven constraints for category C; before re-factoring * Integral-driven constrained PT2 done for category D * Category E singles, onto doubles * Integral-driven constraints done for PT2, all categories; onto re-formatting * doctest * Some comments to make integral-driven code more readable * load balanced constraints; need to test for multiple ranks still * remove print * commit before testing on other branch * commit * parallel working * merge nathan's changes (#10) * External chunk into clean (#61) * Implementation naive algo * pt2 selection now iterates over the full internal space and only a chunk of the connected. selection_step is now an mpi function * fixed doctest * for real this time * for real this time * some initial changes before adding parameter to chunk connected space by input size * spelling * doctest fix * Added functionality for user to specify the size of the working portion of the connected space * Local/global sorts seperate in the selection step + added few more tests * Changes from recent review * Initial PR for chunking the external space (#53) * pt2 selection now iterates over the full internal space and only a chunk of the connected. selection_step is now an mpi function * fixed doctest * for real this time * for real this time * some initial changes before adding parameter to chunk connected space by input size * spelling * doctest fix * Added functionality for user to specify the size of the working portion of the connected space * Local/global sorts seperate in the selection step + added few more tests * Changes from recent review * Constraint based triplet excitations (#7) * Add constraint-based excitation functions (no testing yet) * small name change * now constrained excitation handles opposite spin ones as well; this is a starting point for tests * no way this is working first try? * first test passing * Minimal test added and passing * Constraint-based generation of connected dets from wf + new test * Now constraints working for beta electrons as well * black * black * changes * build fix attempt * attempt 2 * Update test_python.yml * update requirements.txt * Update test_python.yml * remove enumerate * Update test_python.yml --------- Co-authored-by: Thomas Applencourt * only master loads wave function and integrals (#62) * only master loads wave function and integrals * added # because it was missing * nathan PR (#9) * External chunk into clean (#61) * Implementation naive algo * pt2 selection now iterates over the full internal space and only a chunk of the connected. selection_step is now an mpi function * fixed doctest * for real this time * for real this time * some initial changes before adding parameter to chunk connected space by input size * spelling * doctest fix * Added functionality for user to specify the size of the working portion of the connected space * Local/global sorts seperate in the selection step + added few more tests * Changes from recent review * Initial PR for chunking the external space (#53) * pt2 selection now iterates over the full internal space and only a chunk of the connected. selection_step is now an mpi function * fixed doctest * for real this time * for real this time * some initial changes before adding parameter to chunk connected space by input size * spelling * doctest fix * Added functionality for user to specify the size of the working portion of the connected space * Local/global sorts seperate in the selection step + added few more tests * Changes from recent review * Constraint based triplet excitations (#7) * Add constraint-based excitation functions (no testing yet) * small name change * now constrained excitation handles opposite spin ones as well; this is a starting point for tests * no way this is working first try? * first test passing * Minimal test added and passing * Constraint-based generation of connected dets from wf + new test * Now constraints working for beta electrons as well * black * black * changes * build fix attempt * attempt 2 * Update test_python.yml * update requirements.txt * Update test_python.yml * remove enumerate * Update test_python.yml --------- Co-authored-by: Thomas Applencourt * only master loads wave function and integrals (#62) * only master loads wave function and integrals * added # because it was missing --------- Co-authored-by: Thomas Applencourt Co-authored-by: Nathan W Johnson <115745236+nathanjohnsongithub@users.noreply.github.com> --------- Co-authored-by: Thomas Applencourt Co-authored-by: Nathan W Johnson <115745236+nathanjohnsongithub@users.noreply.github.com> * comments * before pr --------- Co-authored-by: Thomas Applencourt Co-authored-by: Nathan W Johnson <115745236+nathanjohnsongithub@users.noreply.github.com> * remove unittest * remove useless value * Flexible types (#12) * initial flexible typing * basic excitation drivers for tuple dets * fix build * move over constraint + phase drivers * mid re-factoring; commit before changing test files to new types * tests passing for tuple determinant * black * Excitation and phase classes fully moved to determinant * doctest * doctest * no more doctests please * Some re-factoring; generic methods moved to generic determinant class * black * Update test_python.yml Need 3.10 for .bit_count() * new io load_wf to handle bitstring representation * Overloaded spindet logical operations, Determinant NamedTuple class re=implemented to be more flexible * doctest + kwargs option for Determinant type * finish updating __new__() for Determinant * TIL about reverse operators * overloaded operators for bitstring type * re-factored types fully working for tuples * docs * io * first review * Phase changes (#13) Loopless phase algorithm. Now agnostic to spin determinant type * fix indexerror in selection * black after resolve merge conflict * comment * Some re-factoring + removal of dependency on tuple * batch drivers * batched operations --------- Co-authored-by: Thomas Applencourt Co-authored-by: Nathan W Johnson <115745236+nathanjohnsongithub@users.noreply.github.com> --- qe/drivers.py | 572 ++++++++++++-------- qe/fundamental_types.py | 774 ++++++++++++++------------- qe/io.py | 5 +- tests/test_everything_all_at_once.py | 38 +- 4 files changed, 766 insertions(+), 623 deletions(-) diff --git a/qe/drivers.py b/qe/drivers.py index 7eccf65..08e4593 100644 --- a/qe/drivers.py +++ b/qe/drivers.py @@ -79,6 +79,142 @@ def integral_category(i, j, k, l): return "G" +# ______ _ _ _____ _ _ _ _ +# | ___ \ | | | | ___| (_) | | | (_) +# | |_/ / |__ __ _ ___ ___ __ _ _ __ __| | | |____ _____ _| |_ __ _| |_ _ ___ _ __ +# | __/| '_ \ / _` / __|/ _ \ / _` | '_ \ / _` | | __\ \/ / __| | __/ _` | __| |/ _ \| '_ \ +# | | | | | | (_| \__ \ __/ | (_| | | | | (_| | | |___> < (__| | || (_| | |_| | (_) | | | | +# \_| |_| |_|\__,_|___/\___| \__,_|_| |_|\__,_| \____/_/\_\___|_|\__\__,_|\__|_|\___/|_| |_| +# +# + +# Drivers for `batched` operations + + +def batch_apply_single_excitation( + batch: List[Determinant], h: OrbitalIdx, p: OrbitalIdx, spin: str +) -> List[Determinant]: + """Apply single excitation to a batch of determinants at a time + :param b: Batch of Determinants to apply h->p excitation to + :param hash: Hash map to check if excited determinant is present in the wave function + :param h: Hole in excitation + :param p: Particle in excitation + :param spintype: Effectively a bool, is the excitation applied to alpha or beta spin determinants + + >>> batch = [Determinant((0, 1), (0, 1)), + ... Determinant((0, 2), (0, 2))] + >>> batch_apply_single_excitation(batch, 0, 3, "alpha") + [Determinant(alpha=(1, 3), beta=(0, 1)), Determinant(alpha=(2, 3), beta=(0, 2))] + >>> batch = [Determinant((0, 1), (0, 1)), + ... Determinant((0, 2), (0, 2))] + >>> batch_apply_single_excitation(batch, 0, 3, "beta") + [Determinant(alpha=(0, 1), beta=(1, 3)), Determinant(alpha=(0, 2), beta=(2, 3))] + """ + + # TODO: Once in C, replace with allocation of space for new determinants... + # For now, we do list comprehension + + return [det_I.apply_single_excitation(h, p, spin) for det_I in batch] + + +def batch_apply_same_spin_double_excitation( + batch: List[Determinant], + h1: OrbitalIdx, + p1: OrbitalIdx, + h2: OrbitalIdx, + p2: OrbitalIdx, + spin: str, +) -> List[Determinant]: + """Apply double (same-spin) excitation to a batch of determinants at a time + :param b: Batch of Determinants to apply h->p excitation to + :param size: Size of batch b + :param h1, h2: Hole in excitation + :param p1, p2: Particle in excitation + :param spintype: Effectively a bool, is the excitation applied to alpha or beta spin determinants + >>> batch = [Determinant((0, 1, 2), (0, 1, 2)), + ... Determinant((0, 2, 3), (0, 2, 3))] + >>> batch_apply_same_spin_double_excitation(batch, 0, 4, 2, 5, "alpha") + [Determinant(alpha=(1, 4, 5), beta=(0, 1, 2)), Determinant(alpha=(3, 4, 5), beta=(0, 2, 3))] + """ + + # TODO: Once in C, replace with allocation of space for new determinants... + # For now, we do list comprehension + + return [det_I.apply_same_spin_double_excitation(h1, p1, h2, p2, spin) for det_I in batch] + + +def batch_apply_opposite_spin_double_excitation( + batch: List[Determinant], h1: OrbitalIdx, p1: OrbitalIdx, h2: OrbitalIdx, p2: OrbitalIdx +) -> List[Determinant]: + """Apply double (same-spin) excitation to a batch of determinants at a time + :param b: Batch of Determinants to apply h->p excitation to + :param size: Size of batch b + :param h1, h2: Hole in excitation + :param p1, p2: Particle in excitation + :param spintype: Effectively a bool, is the excitation applied to alpha or beta spin determinants + >>> batch = [Determinant((0, 1, 2), (0, 1, 2)), + ... Determinant((0, 2, 3), (0, 2, 3))] + >>> batch_apply_opposite_spin_double_excitation(batch, 0, 4, 2, 5) + [Determinant(alpha=(1, 2, 4), beta=(0, 1, 5)), Determinant(alpha=(2, 3, 4), beta=(0, 3, 5))] + """ + + # TODO: Once in C, replace with allocation of space for new determinants... + # For now, we do list comprehension + + return [det_I.apply_opposite_spin_double_excitation(h1, p1, h2, p2) for det_I in batch] + + +def batch_single_phase( + batch: List[Determinant], h: OrbitalIdx, p: OrbitalIdx, spin: str +) -> List[int]: + """Compute phase for batch of determinants related by excitation h <-> p + :param b: Batch of Determinants to apply h->p excitation to + :param size: Size of batch b + :param h: Hole in excitation + :param p: Particle in excitation + :param spintype: Effectively a bool, is the excitation applied to alpha or beta spin determinants + >>> batch = [Determinant((0, 1, 4, 7, 8), ()), + ... Determinant((0, 1, 7, 11, 13), ())] + >>> batch_single_phase(batch, 1, 12, "alpha") + array([-1, 1]) + """ + return np.array([getattr(det, spin).single_phase(h, p) for det in batch], dtype="int") + + +def batch_double_phase( + batch: List[Determinant], + h1: OrbitalIdx, + p1: OrbitalIdx, + h2: OrbitalIdx, + p2: OrbitalIdx, + spin: str, +) -> List[int]: + """Compute phase for batch of determinants related by double (same-spin) excitation h1, h2 <-> p1, p2 + :param b: Batch of Determinants to apply h->p excitation to + :param size: Size of batch b + :param h1, h2: Holes in excitation + :param p1, p2: Particles in excitation + :param spintype: Effectively a bool, is the excitation applied to alpha or beta spin determinants + >>> batch = [Determinant((0, 1, 2, 3, 4, 5, 6, 7, 8), ()), + ... Determinant((0, 1, 2, 3, 4, 5, 6, 7, 13), ())] + >>> batch_double_phase(batch, 2, 11, 3, 12, "alpha") + array([1, 1]) + >>> batch = [Determinant((0, 1, 4, 5, 6, 7, 8), ()), + ... Determinant((0, 1, 4, 5, 6, 7, 13), ())] + >>> batch_double_phase(batch, 1, 11, 4, 3, "alpha") + array([ 1, -1]) + """ + # batch_single_phase returns np.array; supports element-wise multiplication + phase_of_batch = batch_single_phase(batch, h1, p1, spin) * batch_single_phase( + batch, h2, p2, spin + ) + if h2 < p1: + phase_of_batch *= -1 + if p2 < h1: + phase_of_batch *= -1 + return phase_of_batch + + # _ _ _ _ _ _ # | | | | (_) | | (_) # | |_| | __ _ _ __ ___ _| | |_ ___ _ __ _ __ _ _ __ @@ -119,7 +255,7 @@ def H_ij(self, det_i: Determinant, det_j: Determinant) -> Energy: def H_ij_spindet(sdet_j, spin: str) -> Energy: """, when I and J differ by exactly one orbital.""" # det_i is accessible - phase, m, p = det_i.single_exc(sdet_j, spin) + phase, m, p = getattr(det_i, spin).single_exc(sdet_j) return self.H_ij_orbital(m, p) * phase ed_up, ed_dn = det_i.exc_degree(det_j) @@ -222,7 +358,7 @@ def H_ij_single_indices( spin: str, ) -> Iterator[Two_electron_integral_index_phase]: """, when I and J differ by exactly one orbital""" - phase, h, p = det_i.single_exc(sdet_j, spin) + phase, h, p = getattr(det_i, spin).single_exc(sdet_j) for i in sdet_i: yield (h, i, p, i), phase @@ -235,7 +371,7 @@ def H_ij_doubleAA_indices( ) -> Iterator[Two_electron_integral_index_phase]: """, when I and J differ by exactly two orbitals within the same spin.""" - phase, h1, h2, p1, p2 = det_i.double_exc(sdet_j, spin) + phase, h1, h2, p1, p2 = getattr(det_i, spin).double_exc(sdet_j) yield (h1, h2, p1, p2), phase yield (h1, h2, p2, p1), -phase @@ -246,8 +382,8 @@ def H_ij_doubleAB_2e_index( """, when I and J differ by exactly one alpha spin-orbital and one beta spin-orbital.""" - phaseA, h1, p1 = det_i.single_exc(det_j.alpha, "alpha") - phaseB, h2, p2 = det_i.single_exc(det_j.beta, "beta") + phaseA, h1, p1 = det_i.alpha.single_exc(det_j.alpha) + phaseB, h2, p2 = det_i.beta.single_exc(det_j.beta) yield (h1, h2, p1, p2), phaseA * phaseB @@ -350,6 +486,7 @@ def get_dets_via_orbital_occupancy( ) # if len(det_indices) == 0, `set()` is returned + # TODO: Possible, return Determinants as opposed to indices? return det_indices @staticmethod @@ -383,25 +520,28 @@ def do_single( Called by category functions corresponding to single excitations For use in building the Hamiltonian in the variational step""" - for I in det_indices: - det = psi_internal[I] # Grab |Determinant| I from w.f. - if spin == "alpha": - # Create new |Determinant| via excitation from h -> p - # Hole, particle are alpha spin orbitals - excited_det = det.apply_excitation(((h,), (p,)), ((), ())) - else: # Spin is "beta" - # Hole, particle are beta spin orbitals - excited_det = det.apply_excitation(((), ()), ((h,), (p,))) - # If the excited determinant is in the internal wave function -> - # Compute phase of (single) excitation pair and yield (I, J), phase - if excited_det in det_to_index: - # :param `phasemod` is \pm 1, sign of two-electron integral in definition of SC-rules - - phase = phasemod * det.single_phase(h, p, spin) - - yield (I, det_to_index[excited_det]), phase - else: - pass + # TODO: Just for proof of concept, can do this outside the function later + # det_indices is from itertools, so it is de-allocated after one pass through it. So, yield (I, D_I) pairs + batch_of_dets = [(I, psi_internal[I]) for I in det_indices] + # In order + # 1. Apply simultaneous single h -> p excitation to pre-filtered determinants in `batch_of_dets` + exc_batch = batch_apply_single_excitation( + [index_det_pair[1] for index_det_pair in batch_of_dets], h, p, spin + ) + # 2. Batch filter of excitation pairs based on whether excitation is \in psi_internal + filtered_pairs = [ + (index_det_pair[0], det_to_index[exc_det_I]) + for index_det_pair, exc_det_I in zip(batch_of_dets, exc_batch) + if (exc_det_I in det_to_index) + ] + # 3. Compute phase for filtered pairs + # (phasemod is \pm 1; accounts for whether integral is coulomb or exchange) + phase_of_batch = phasemod * batch_single_phase( + [psi_internal[I] for (I, _) in filtered_pairs], h, p, spin + ) + # Yield (I, J), phase pairs for computing + for (I, J), phase in zip(filtered_pairs, phase_of_batch): + yield (I, J), phase @staticmethod def do_single_pt2( @@ -419,24 +559,25 @@ def do_single_pt2( Called by category functions corresponding to single excitations For use in computing E_pt2 and subsequent selection""" - for I in det_indices: - det = psi_internal[I] # Grab |Determinant| I from w.f. - if spin == "alpha": - # Create new |Determinant| via excitation from h -> p - # Hole, particle are alpha spin orbitals - excited_det = det.apply_excitation(((h,), (p,)), ((), ())) - else: # Spin is "beta" - # Hole, particle are beta spin orbitals - excited_det = det.apply_excitation(((), ()), ((h,), (p,))) - # Assert the excited determinant satisfies the appropriate constraint C - assert check_constraint(excited_det) == C - # param `phasemod` is \pm 1, sign of two-electron integral in definition of SC-rules - - phase = phasemod * det.single_phase(h, p, spin) - - # excited_det is in the connected space by default; yield (I, det_J) phase - # Cannot yield index of excited_det since this requires storing (knowing) connected space a priori - yield (I, excited_det), phase + # det_indices is from itertools, so it is de-allocated after one pass through it. So, yield (I, D_I) pairs + batch_of_dets = [(I, psi_internal[I]) for I in det_indices] + # In order + # 1. Apply simultaneous single h -> p excitation to pre-filtered determinants in `batch_of_dets` + exc_batch = batch_apply_single_excitation( + [index_det_pair[1] for index_det_pair in batch_of_dets], h, p, spin + ) + # 2. Organize excitation pairs of (I, det_J); I is index of det_I \in psi_internal, det_J in connected + int_ext_pairs = [ + (int_index_det_pair[0], exc_det_I) + for int_index_det_pair, exc_det_I in zip(batch_of_dets, exc_batch) + ] + # 3. Compute phase for batch excitation pairs + phase_of_batch = phasemod * batch_single_phase( + [det_I for (_, det_I) in batch_of_dets], h, p, spin + ) + # Yield (I, J), phase pairs for computing + for (I, det_J), phase in zip(int_ext_pairs, phase_of_batch): + yield (I, det_J), phase @staticmethod def do_double_samespin( @@ -455,30 +596,33 @@ def do_double_samespin( # Unpack hole-particle pairs h1, p1 = hp1 h2, p2 = hp2 - # Pre-filter; get list of candidate determinants thaat are possibly related via (same-spin) excitations h1 -> p1, h2 -> p2 + # Pre-filter; get list of candidate determinants that are possibly related via (same-spin) excitations h1 -> p1, h2 -> p2 det_indices_AA = Hamiltonian_two_electrons_integral_driven.get_dets_via_orbital_occupancy( - spindet_occ_i, {}, {"same": {h1, h2}}, {"same": {p1, p1}} + spindet_occ_i, {}, {"same": {h1, h2}}, {"same": {p1, p2}} ) - for I in det_indices_AA: - det = psi_internal[I] - if spin == "alpha": - # Create new |Determinant| via excitation from h1 -> p1, h2 -> p1 - excited_det = det.apply_excitation(((h1, h2), (p1, p2)), ((), ())) - else: - # Holes, particles are beta spin orbitals - excited_det = det.apply_excitation(((), ()), ((h1, h2), (p1, p2))) - # If the excited determinant is in the internal wave function -> - # Compute phase of (double) excitation pair and yield (I, J), phase - if excited_det in det_to_index: - phase = det.double_phase((h1, h2), (p1, p2), spin) - if h2 < h1: - phase = -phase - if p2 < p1: - phase = -phase - - yield (I, det_to_index[excited_det]), phase - else: - pass + # Get (I, D_I) pairs of dets indicaated by det_indices_AA + batch_of_dets = [(I, psi_internal[I]) for I in det_indices_AA] + # In order + # 1. Apply simultaneous single h -> p excitation to pre-filtered determinants in `batch_of_dets` + exc_batch = batch_apply_same_spin_double_excitation( + [index_det_pair[1] for index_det_pair in batch_of_dets], h1, p1, h2, p2, spin + ) + # 2. Batch filter of excitation pairs based on whether excitation is \in psi_internal + filtered_pairs = [ + (index_det_pair[0], det_to_index[exc_det_I]) + for index_det_pair, exc_det_I in zip(batch_of_dets, exc_batch) + if (exc_det_I in det_to_index) + ] + # 3. Compute phase for filtered pairs + phase_of_batch = batch_double_phase( + [psi_internal[I] for (I, _) in filtered_pairs], h1, p1, h2, p2, spin + ) + # For exchange integrals; + if np.sign(h2 - h1) != np.sign(p2 - p1): + phase_of_batch *= -1 + # Yield (I, J), phase pairs for computing + for (I, J), phase in zip(filtered_pairs, phase_of_batch): + yield (I, J), phase @staticmethod def do_double_samespin_pt2( @@ -518,7 +662,7 @@ def do_double_samespin_pt2( # Related pairs must be: # Occupied in: (alpha) h1, h2 C = {a1, a2, a3} (beta) none # Empty in: (alpha) p1, p2 {a1 + 1, a1 + 2, ... N_orb - 1} - {a1, a2, a3, h1, h2} (beta) none - det_indices = ( + det_indices_AA = ( Hamiltonian_two_electrons_integral_driven.get_dets_via_orbital_occupancy( spindet_occ, {}, @@ -536,7 +680,7 @@ def do_double_samespin_pt2( # Related pairs must be: # Occupied in: (alpha) h1, h2 C - {p2} (beta) none # Empty in: (alpha) p1, p2 {a1 + 1, a1 + 2, ... N_orb - 1} - ((C - {p2}) | {h1, h2}) (beta) none - det_indices = ( + det_indices_AA = ( Hamiltonian_two_electrons_integral_driven.get_dets_via_orbital_occupancy( spindet_occ, {}, @@ -553,7 +697,7 @@ def do_double_samespin_pt2( # Related pairs must be: # Occupied in: (alpha) h1, h2 C - {p1} (beta) none # Empty in: (alpha) p1, p2 {a1 + 1, a1 + 2, ... N_orb - 1} - ((C - {p1}) | {h1, h2}) (beta) none - det_indices = ( + det_indices_AA = ( Hamiltonian_two_electrons_integral_driven.get_dets_via_orbital_occupancy( spindet_occ, {}, @@ -569,7 +713,7 @@ def do_double_samespin_pt2( # Related pairs must be: # Occupied in: (alpha) h1, h2 C - {p1, p2} (beta) none # Empty in: (alpha) p1, p2 {a1 + 1, a1 + 2, ... N_orb - 1} - ((C - {p1, p2}) | {h1, h2}) (beta) none - det_indices = ( + det_indices_AA = ( Hamiltonian_two_electrons_integral_driven.get_dets_via_orbital_occupancy( spindet_occ, {}, @@ -583,30 +727,39 @@ def do_double_samespin_pt2( # Related pairs must be: # Occupied in: (alpha) C (beta) h1, h2 # Empty in: (alpha) {a1 + 1, a1 + 2, ... N_orb - 1} - C | {h1, h2}) (beta) p1, p2 - det_indices = Hamiltonian_two_electrons_integral_driven.get_dets_via_orbital_occupancy( - spindet_occ, - oppspindet_occ, - {"same": {h1, h2}, "opposite": set(C)}, - {"same": {p1, p2}, "opposite": unocc_oppspin_orbitals}, - ) - - # Now, loop through the filtered `candidate` determinants; generate the excitation pairs related by h1 -> p1, h2 -> p2 - for I in det_indices: - det = psi[I] - if spin == "alpha": - # Create new |Determinant| via excitation from i -> j, k -> l - excited_det = det.apply_excitation(((h1, h2), (p1, p2)), ((), ())) - else: - # Holes, particles are beta spin orbitals - excited_det = det.apply_excitation(((), ()), ((h1, h2), (p1, p2))) - # Assert the excited determinant satisfies the appropriate constraint C - assert check_constraint(excited_det) == C - phase = det.double_phase((h1, h2), (p1, p2), spin) - if h2 < h1: - phase = -phase - if p2 < p1: - phase = -phase - yield (I, excited_det), phase + det_indices_AA = ( + Hamiltonian_two_electrons_integral_driven.get_dets_via_orbital_occupancy( + spindet_occ, + oppspindet_occ, + {"same": {h1, h2}, "opposite": set(C)}, + {"same": {p1, p2}, "opposite": unocc_oppspin_orbitals}, + ) + ) + + # Organized filtered dets into index, determinant (I, D_I) pairs + batch_of_dets = [(I, psi[I]) for I in det_indices_AA] + # In order + # 1. Apply simultaneous single h1, h2 -> p1, p2 (same-spin) excitation to pre-filtered determinants in `batch_of_dets` + exc_batch = batch_apply_same_spin_double_excitation( + [index_det_pair[1] for index_det_pair in batch_of_dets], h1, p1, h2, p2, spin + ) + # 2. Organize excitation pairs of (I, det_J); + # a. I is index of det_I \in psi + # b. det_J is the determinant connected to det_I via h1, h2 -> p1, p2 + int_ext_pairs = [ + (int_index_det_pair[0], exc_det_I) + for int_index_det_pair, exc_det_I in zip(batch_of_dets, exc_batch) + ] + # 3. Compute phase for filtered pairs + phase_of_batch = batch_double_phase( + [psi[I] for (I, _) in int_ext_pairs], h1, p1, h2, p2, spin + ) + # For exchange integrals; + if np.sign(h2 - h1) != np.sign(p2 - p1): + phase_of_batch *= -1 + # Yield (I, J), phase pairs for computing + for (I, det_J), phase in zip(int_ext_pairs, phase_of_batch): + yield (I, det_J), phase @staticmethod def do_double_oppspin( @@ -632,25 +785,38 @@ def do_double_oppspin( {"same": {h1}, "opposite": {h2}}, {"same": {p1}, "opposite": {p2}}, ) - for I in det_indices_AB: - det = psi_internal[I] - if spin == "alpha": - # h1, p1 -> alpha spin-orbitals, h2, p2 -> beta - excited_det = det.apply_excitation(((h1,), (p1,)), ((h2,), (p2,))) - # Phase is computed from both spindets - phaseA = det.single_phase(h1, p1, spin) - phaseB = det.single_phase(h2, p2, "beta") - else: - # h1, p1 -> beta spin-orbitals, h2, p2 -> alpha - excited_det = det.apply_excitation(((h2,), (p2,)), ((h1,), (p1,))) - phaseA = det.single_phase(h1, p1, spin) - phaseB = det.single_phase(h2, p2, "alpha") - # If the excited determinant is in the internal wave function -> - # Compute phase of (double) excitation pair and yield (I, J), phase - if excited_det in det_to_index: - yield (I, det_to_index[excited_det]), phaseA * phaseB - else: - pass + # Get (I, D_I) pairs of dets indicaated by det_indices_AA + batch_of_dets = [(I, psi_internal[I]) for I in det_indices_AB] + # In order + # 1. Apply simultaneous single ha, hb -> pa, pb excitation to pre-filtered determinants in `batch_of_dets` + if spin == "alpha": + exc_batch = batch_apply_opposite_spin_double_excitation( + [index_det_pair[1] for index_det_pair in batch_of_dets], h1, p1, h2, p2 + ) + oppspin = "beta" + else: # Spin is `beta` + exc_batch = batch_apply_opposite_spin_double_excitation( + [index_det_pair[1] for index_det_pair in batch_of_dets], h2, p2, h1, p1 + ) + oppspin = "alpha" + # 2. Batch filter of excitation pairs based on whether excitation is \in psi_internal + filtered_pairs = [ + (index_det_pair[0], det_to_index[exc_det_I]) + for index_det_pair, exc_det_I in zip(batch_of_dets, exc_batch) + if (exc_det_I in det_to_index) + ] + # 3. Compute phase for filtered pairs + phaseA_of_batch = batch_single_phase( + [psi_internal[I] for (I, _) in filtered_pairs], h1, p1, spin + ) + phaseB_of_batch = batch_single_phase( + [psi_internal[I] for (I, _) in filtered_pairs], h2, p2, oppspin + ) + # Element-wise multiplication of phaseA, phaseB + phase_of_batch = phaseA_of_batch * phaseB_of_batch + # Yield (I, J), phase pairs for computing + for (I, J), phase in zip(filtered_pairs, phase_of_batch): + yield (I, J), phase @staticmethod def do_double_oppspin_pt2( @@ -687,7 +853,7 @@ def do_double_oppspin_pt2( # Related pairs must be: # Occupied in: (alpha) h1, C = {a1, a2, a3} (beta) h2 # Empty in: (alpha) p1, {a1 + 1, a1 + 2, ... N_orb - 1} - {a1, a2, a3, h1} (beta) p2 - det_indices = ( + det_indices_AB = ( Hamiltonian_two_electrons_integral_driven.get_dets_via_orbital_occupancy( spindet_occ, oppspindet_occ, @@ -702,7 +868,7 @@ def do_double_oppspin_pt2( # Related pairs must be: # Occupied in: (alpha) h1, C - {p1} (beta) h2 # Empty in: (alpha) p1, {a1 + 1, a1 + 2, ... N_orb - 1} - ((C - {p1}) | {h1}) (beta) p2 - det_indices = ( + det_indices_AB = ( Hamiltonian_two_electrons_integral_driven.get_dets_via_orbital_occupancy( spindet_occ, oppspindet_occ, @@ -724,7 +890,7 @@ def do_double_oppspin_pt2( # Related pairs must be: # Occupied in: (alpha) h2, C = {a1, a2, a3} (beta) h1 # Empty in: (alpha) p2, {a1 + 1, a1 + 2, ... N_orb - 1} - {a1, a2, a3, h2} (beta) p1 - det_indices = ( + det_indices_AB = ( Hamiltonian_two_electrons_integral_driven.get_dets_via_orbital_occupancy( spindet_occ, oppspindet_occ, @@ -739,7 +905,7 @@ def do_double_oppspin_pt2( # Related pairs must be: # Occupied in: (alpha) h2, C - {p2} (beta) h1 # Empty in: (alpha) p2, {a1 + 1, a1 + 2, ... N_orb - 1} - ((C - {p2}) | {h2}) (beta) p1 - det_indices = ( + det_indices_AB = ( Hamiltonian_two_electrons_integral_driven.get_dets_via_orbital_occupancy( spindet_occ, oppspindet_occ, @@ -747,23 +913,36 @@ def do_double_oppspin_pt2( {"same": {p1}, "opposite": unocc_spin_orbitals}, ) ) - # Now, loop through filtered dets to yield excitation pairs s.to C - for I in det_indices: - det = psi[I] - if spin == "alpha": - # h1, p1 -> alpha spin-orbitals, h2, p2 -> beta - excited_det = det.apply_excitation(((h1,), (p1,)), ((h2,), (p2,))) - # Phase is computed from both spindets - phaseA = det.single_phase(h1, p1, spin) - phaseB = det.single_phase(h2, p2, "beta") - else: - # h1, p1 -> beta spin-orbitals, h2, p2 -> alpha - excited_det = det.apply_excitation(((h2,), (p2,)), ((h1,), (p1,))) - phaseA = det.single_phase(h1, p1, spin) - phaseB = det.single_phase(h2, p2, "alpha") - # Assert excited det satisfies constraint and yield - assert check_constraint(excited_det) == C - yield (I, excited_det), phaseA * phaseB + + # Organized filtered dets into index, determinant (I, D_I) pairs + batch_of_dets = [(I, psi[I]) for I in det_indices_AB] + # In order + # 1. Apply simultaneous single ha, hb -> pa, pb excitation to pre-filtered determinants in `batch_of_dets` + if spin == "alpha": + exc_batch = batch_apply_opposite_spin_double_excitation( + [index_det_pair[1] for index_det_pair in batch_of_dets], h1, p1, h2, p2 + ) + oppspin = "beta" + else: # Spin is `beta` + exc_batch = batch_apply_opposite_spin_double_excitation( + [index_det_pair[1] for index_det_pair in batch_of_dets], h2, p2, h1, p1 + ) + oppspin = "alpha" + # 2. Organize excitation pairs of (I, det_J); + # a. I is index of det_I \in psi + # b. det_J is the determinant connected to det_I via h1, h2 -> p1, p2 + int_ext_pairs = [ + (int_index_det_pair[0], exc_det_I) + for int_index_det_pair, exc_det_I in zip(batch_of_dets, exc_batch) + ] + # 3. Compute phase for filtered pairs + phaseA_of_batch = batch_single_phase([psi[I] for (I, _) in int_ext_pairs], h1, p1, spin) + phaseB_of_batch = batch_single_phase([psi[I] for (I, _) in int_ext_pairs], h2, p2, oppspin) + # Element-wise multiplication of phaseA, phaseB + phase_of_batch = phaseA_of_batch * phaseB_of_batch + # Yield (I, J), phase pairs for computing + for (I, det_J), phase in zip(int_ext_pairs, phase_of_batch): + yield (I, det_J), phase @staticmethod def category_A( @@ -2714,7 +2893,6 @@ def psi_external_pt2( E_pt2_J = nominator_conts_table.values() nominator_conts = np.array(list(E_pt2_J), dtype="float") # TODO: For integral driven, loop over integrals? In general, be more efficient in this area. - # TODO: In general, though, we need a different way to do this. It goes back to storing the connected space. We can even do the denominator conts in a separate array.. psi_connected_C = [det_J for det_J in nominator_conts_table.keys()] denominator_conts = np.divide( 1.0, @@ -2879,77 +3057,6 @@ def global_sort_pt2_energies(comm, local_best_dets: Psi_det, local_best_energies # Next, we have functions that split the connected space -def get_chunk_of_connected_determinants(psi_det: Psi_det, n_orb: int, L=None) -> Iterator[Psi_det]: - """ - MPI function, generates chunks of connected determinants of size L - - Inputs: - :param psi_det: list of determinants - :param L: integer, maximally allowed `chunk' of the conneceted space to yield at a time - default is L = None, if no chunk size specified, a chunk of the connected space is allocated to each rank - - >>> d1 = Determinant((0, 1), (0,) ) ; d2 = Determinant((0, 2), (0,) ) - >>> for psi_chunk in get_chunk_of_connected_determinants( [ d1,d2 ], 4): - ... len(psi_chunk) - 22 - >>> for psi_chunk in get_chunk_of_connected_determinants( [ d1,d2 ], 4, 11 ): - ... len(psi_chunk) - 11 - 11 - >>> for psi_chunk in get_chunk_of_connected_determinants( [ d1,d2 ], 4, 10 ): - ... len(psi_chunk) - 10 - 10 - 2 - """ - - def gen_all_connected_determinant(psi_det: Psi_det, n_orb: int) -> Psi_det: - """ - >>> d1 = Determinant((0, 1), (0,) ) ; d2 = Determinant((0, 2), (0,) ) - >>> len(gen_all_connected_determinant( [ d1,d2 ], 4 )) - 22 - - We remove the connected determinant who are already inside the wave function. Order doesn't matter - """ - # Literal programing - # return set(chain.from_iterable(map(self.gen_all_connected_det_from_det, psi_det)))- set(psi_det) - - # Naive algorithm 13 - l_global = [] - for i, det in enumerate(psi_det): - for det_connected in det.gen_all_connected_det(n_orb): - # Remove determinant who are in psi_det - if det_connected in psi_det: - continue - # If it's was already generated by an old determinant, just drop it - if any(det_connected.is_connected(d) for d in psi_det[:i]): - continue - - l_global.append(det_connected) - - # Return connected space - return l_global - - # Naive: Each ranks generates all connected determinants, and takes what is theirs - psi_connected = gen_all_connected_determinant(psi_det, n_orb) - world_size = MPI.COMM_WORLD.Get_size() - # TODO: len(psi_connected) will not scale, but is needed for this naive representation - full_idx = np.arange(len(psi_connected)) - split_idx = np.array_split(full_idx, world_size) - # Part of connected space available to each rank - full_chunk = [psi_connected[i] for i in split_idx[MPI.COMM_WORLD.Get_rank()]] - - if L is None: # If no argument passed by the user - L = len(full_chunk) - number_of_chunks, leftovers = divmod(len(full_chunk), L) - # Yield chunks of size L one at a time - for i in np.arange(number_of_chunks): - yield full_chunk[i * L : (i + 1) * L] - # Yield `leftover' determinants (size of chunk < L) - if leftovers > 0: # If there are leftover determinants to handle - yield full_chunk[number_of_chunks * L : number_of_chunks * L + leftovers] - - def generate_all_constraints(n_elec: int, n_orb: int, m=3) -> List[Tuple[OrbitalIdx, ...]]: """Generate all `m'-constraints, characterized by m most highly occupied (usually alpha) electrons >>> generate_all_constraints(3, 4) @@ -2989,70 +3096,75 @@ def dispatch_local_constraints( # Initialize array to track workload of each rank W = np.zeros(shape=(comm.Get_size(),), dtype="i") C_loc = [] # Pre-allocate space for local constraints - na = len(getattr(psi[0], "alpha")) # No. of alpha electrons - nb = len(getattr(psi[0], "beta")) # No. of beta electrons + na = getattr(psi[0], "alpha").popcnt() # No. of alpha electrons + nb = getattr(psi[0], "beta").popcnt() # No. of beta electrons # Pass through all triplet constraints to distribute H = [] # Track work dist. for C in generate_all_constraints(na, n_orb): - B_upper = set(range(min(C) + 1, n_orb)) # Upper bitmask - B_lower = set(range(min(C))) # Lower bitmask + B_upper = tuple(range(min(C) + 1, n_orb)) # Upper bitmask + B_lower = tuple(range(min(C))) # Lower bitmask h = 0 # Track work of this constraint for det in psi: det_a = getattr(det, "alpha") - constraint_orbitals_occupied = set(det_a) & set(C) - nonconstrained_orbitals_occupied = (set(det_a) & B_upper) - set(C) + constraint_orbitals_occupied = det_a & C + # Operations below is equivalent to `det_a & B_upper - C` + nonconstrained_orbitals_occupied = (det_a & B_upper).get_holes(C) # n_particles = [np_a, np_b, np_aa, np_bb, np_ab] # Number of particles (or pairs) that (could possibly) involve an excitation satisfying C - if len(constraint_orbitals_occupied) == 0: + if constraint_orbitals_occupied.popcnt() == 0: # No excitations will satisfy C -> Pass n_particles = np.zeros(5, dtype="i") - elif len(constraint_orbitals_occupied) == 1: + elif constraint_orbitals_occupied.popcnt() == 1: # To satisfy C, excitation must be aa (into empty constraint orbitals) n_particles = np.array([0, 0, 1, 0, 0], dtype="i") - elif len(constraint_orbitals_occupied) == 2: + elif constraint_orbitals_occupied.popcnt() == 2: # To satisfy C, only possible single is a, must excite into empty constraint orbital - na_orbs_unocc_lower = B_lower - set(det_a) + # Operations below is equivalent to `B_lower - det_a` in set operations + na_orbs_unocc_lower = det_a.get_particles(B_lower) # No bb; for ab doubles, a must excite into empty constraint orbital, so 1 * (self.n_orb - nb) ab pairs - n_particles = np.array([1, 0, len(na_orbs_unocc_lower), 0, (n_orb - nb)], dtype="i") - elif len(constraint_orbitals_occupied) == 3: + n_particles = np.array( + [1, 0, na_orbs_unocc_lower.popcnt(), 0, (n_orb - nb)], dtype="i" + ) + elif constraint_orbitals_occupied.popcnt() == 3: # To satisfy C, any a or aa excitaion into `lower` unoccupied alpha orbitals # All possible b or bb excitations satisfy C - na_orbs_unocc_lower = B_lower - set(det_a) + # Operations below is equivalent to `B_lower - det_a` in set operations + na_orbs_unocc_lower = det_a.get_particles(B_lower) # Divide some things by 2 to avoid repeats due to permuation (e.g., (p1, p2) = (1, 2) <-> (p1, p2) = (2, 1)) n_particles = np.array( [ - len(na_orbs_unocc_lower), + na_orbs_unocc_lower.popcnt(), n_orb - nb, - len(na_orbs_unocc_lower) * (len(na_orbs_unocc_lower) - 1) / 2, + na_orbs_unocc_lower.popcnt() * (na_orbs_unocc_lower.popcnt() - 1) / 2, (n_orb - nb) * (n_orb - nb - 1) / 2, - (n_orb - nb) * len(na_orbs_unocc_lower), + (n_orb - nb) * na_orbs_unocc_lower.popcnt(), ], dtype="i", ) # n_holes = [nh_a, nh_b, nh_aa, nh_bb, nh_ab] # Number of holes (or pairs) that (could possibly) involve an excitation satisfying C - if len(nonconstrained_orbitals_occupied) > 2: + if nonconstrained_orbitals_occupied.popcnt() > 2: # No excitations will satisfy C -> Pass n_holes = np.zeros(5, dtype="i") - elif len(nonconstrained_orbitals_occupied) == 2: + elif nonconstrained_orbitals_occupied.popcnt() == 2: # To satisfy C, excitation must be aa (out of `higher` non-constraint orbitals) n_holes = np.array([0, 0, 1, 0, 0], dtype="i") - elif len(nonconstrained_orbitals_occupied) == 1: + elif nonconstrained_orbitals_occupied.popcnt() == 1: # To satisfy C, only possible single is a, must excite out of `higher` constraint orbitals - na_orbs_occ_lower = set(det_a) & B_lower + na_orbs_occ_lower = det_a & B_lower # No bb; for ab doubles, a must excite out of `higher` constraint orbital, so 1 * nb ab pairs - n_holes = [1, 0, len(na_orbs_occ_lower), 0, nb] - elif len(nonconstrained_orbitals_occupied) == 0: + n_holes = [1, 0, na_orbs_occ_lower.popcnt(), 0, nb] + elif nonconstrained_orbitals_occupied.popcnt() == 0: # To satisfy C, can excite out of any `lower` occupied alpha orbitals # All possible b or bb excitations satisfy C - na_orbs_occ_lower = set(det_a) & B_lower + na_orbs_occ_lower = det_a & B_lower n_holes = [ - len(na_orbs_occ_lower), + na_orbs_occ_lower.popcnt(), nb, - len(na_orbs_occ_lower) * (len(na_orbs_occ_lower) - 1) / 2, + na_orbs_occ_lower.popcnt() * (na_orbs_occ_lower.popcnt() - 1) / 2, nb * (nb - 1) / 2, - len(na_orbs_occ_lower) * nb, + na_orbs_occ_lower.popcnt() * nb, ] # Number of singly/doubly connected determinants to |det> satisfying constraint C diff --git a/qe/fundamental_types.py b/qe/fundamental_types.py index 2dbe2fc..d257f16 100644 --- a/qe/fundamental_types.py +++ b/qe/fundamental_types.py @@ -1,3 +1,6 @@ +# For forward declaration in type hints +from __future__ import annotations + from typing import Tuple, Dict, NamedTuple, List, NewType, Iterator from itertools import chain, product, combinations, takewhile from functools import partial, cache @@ -15,6 +18,18 @@ # $ = \int \phi_i(r) (-\frac{1}{2} \Delta + V_en ) \phi_k(r) dr$ One_electron_integral = Dict[Tuple[OrbitalIdx, OrbitalIdx], float] +# # _ +# # |_) o _|_ _ _ _|_ +# # |_) | |_ _> (/_ |_ +# # + +# class Spin_determinant_bitset: +# def __init__(self, int): +# # TODO: Define these things +# self.tag = qelib.SPIN_DET_TYPE_BITSET +# npa = np.array(sorted(t), dtype=np.intc) +# npa_size = len(t) +# self.handle = qelib.qe_spin_det_bitset_create(npa, npa_size) # ___ # | ._ | _ @@ -45,7 +60,9 @@ def convert_repr(self, Norb): # Return |Bitstring| representation of this |Determinant| return Spin_determinant_bitstring(int(("".join(bitstr)), 2)) - def __and__(self, s_tuple: Tuple[OrbitalIdx, ...]) -> Tuple[OrbitalIdx, ...]: + # Use iterator methods inherent to |Tuple|; no need to overload as in bitset representation + + def __and__(self, s_tuple: Spin_determinant_tuple) -> Spin_determinant_tuple: """Overload `&` operator to perform set intersection Return type |Spin_determinant_tuple| >>> Spin_determinant_tuple((0, 1)) & Spin_determinant_tuple((0, 2)) @@ -55,7 +72,7 @@ def __and__(self, s_tuple: Tuple[OrbitalIdx, ...]) -> Tuple[OrbitalIdx, ...]: """ return Spin_determinant_tuple(sorted(set(self) & set(s_tuple))) - def __rand__(self, s_tuple: Tuple[OrbitalIdx, ...]) -> Tuple[OrbitalIdx]: + def __rand__(self, s_tuple: Spin_determinant_tuple) -> Spin_determinant_tuple: """Reverse overloaded __and__ >>> (0, 1) & Spin_determinant_tuple((0, 2)) (0,) @@ -64,7 +81,7 @@ def __rand__(self, s_tuple: Tuple[OrbitalIdx, ...]) -> Tuple[OrbitalIdx]: """ return self.__and__(s_tuple) - def __or__(self, s_tuple: Tuple[OrbitalIdx, ...]) -> Tuple[OrbitalIdx, ...]: + def __or__(self, s_tuple: Spin_determinant_tuple) -> Spin_determinant_tuple: """Overload `|` operator to perform set union Return type |Spin_determinant_tuple| >>> Spin_determinant_tuple((0, 1)) | Spin_determinant_tuple((0, 2)) @@ -76,14 +93,14 @@ def __or__(self, s_tuple: Tuple[OrbitalIdx, ...]) -> Tuple[OrbitalIdx, ...]: """ return Spin_determinant_tuple(sorted(set(self) | set(s_tuple))) - def __ror__(self, s_tuple: Tuple[OrbitalIdx, ...]) -> Tuple[OrbitalIdx]: + def __ror__(self, s_tuple: Spin_determinant_tuple) -> Spin_determinant_tuple: """Reverse overloaded __or__ >>> (0, 1) | Spin_determinant_tuple((0, 2)) (0, 1, 2) """ return self.__or__(s_tuple) - def __xor__(self, s_tuple: Tuple[OrbitalIdx, ...]) -> Tuple[OrbitalIdx, ...]: + def __xor__(self, s_tuple: Spin_determinant_tuple) -> Spin_determinant_tuple: """Overload `^` operator to perform symmetric set difference Return type |Spin_determinant_tuple| >>> Spin_determinant_tuple((0, 1)) ^ Spin_determinant_tuple((0, 2)) @@ -95,14 +112,14 @@ def __xor__(self, s_tuple: Tuple[OrbitalIdx, ...]) -> Tuple[OrbitalIdx, ...]: """ return Spin_determinant_tuple(sorted(set(self) ^ set(s_tuple))) - def __rxor__(self, s_tuple: Tuple[OrbitalIdx, ...]) -> Tuple[OrbitalIdx]: + def __rxor__(self, s_tuple: Spin_determinant_tuple) -> Spin_determinant_tuple: """Reverse overloaded __xor__ >>> (0, 1) ^ Spin_determinant_tuple((0, 2)) (1, 2) """ return self.__xor__(s_tuple) - def __sub__(self, s_tuple: Tuple[OrbitalIdx, ...]) -> Tuple[OrbitalIdx, ...]: + def __sub__(self, s_tuple: Spin_determinant_tuple) -> Spin_determinant_tuple: """Overload `-` operator to perform set difference Return type |Spin_determinant_tuple| >>> Spin_determinant_tuple((0, 1)) - Spin_determinant_tuple((0, 2)) @@ -112,9 +129,9 @@ def __sub__(self, s_tuple: Tuple[OrbitalIdx, ...]) -> Tuple[OrbitalIdx, ...]: >>> Spin_determinant_tuple((0, 1)) - Spin_determinant_tuple((0, 1)) () """ - return Spin_determinant_tuple(set(self) - set(s_tuple)) + return Spin_determinant_tuple(sorted(set(self) - set(s_tuple))) - def __rsub__(self, s_tuple: Tuple[OrbitalIdx, ...]) -> Tuple[OrbitalIdx]: + def __rsub__(self, s_tuple: Spin_determinant_tuple) -> Spin_determinant_tuple: """Reverse overloaded __sub__ Convert arg `s_tuple' to |Spin_determinant_tuple|, then perform __sub__ since operation is not communative >>> (0, 1, 2, 3) - Spin_determinant_tuple((0, 2)) @@ -124,9 +141,89 @@ def __rsub__(self, s_tuple: Tuple[OrbitalIdx, ...]) -> Tuple[OrbitalIdx]: def popcnt(self) -> int: """Perform a `popcount'; return length of the tuple""" - return len(self) + return int(len(self)) + + # _ + # |_ _ o _|_ _. _|_ o _ ._ + # |_ >< (_ | |_ (_| |_ | (_) | | + # + + def apply_single_excitation( + self, hole: OrbitalIdx, particle: OrbitalIdx + ) -> Spin_determinant_tuple: + """Apply single hole -> particle excitation to instance of |Spin_determinant| + >>> Spin_determinant_tuple((0, 1)).apply_single_excitation(0, 2) + (1, 2) + >>> Spin_determinant_tuple((1, 4)).apply_single_excitation(4, 0) + (0, 1) + >>> Spin_determinant_tuple((0, 1, 3, 4, 7, 10)).apply_single_excitation(4, 8) + (0, 1, 3, 7, 8, 10) + >>> Spin_determinant_tuple((0, 1, 3, 4, 7, 10)).apply_single_excitation(4, 5) + (0, 1, 3, 5, 7, 10) + >>> Spin_determinant_tuple((0, 1, 3, 4, 7, 10)).apply_single_excitation(7, 2) + (0, 1, 2, 3, 4, 10) + >>> Spin_determinant_tuple((0, 1, 3, 4, 7, 10)).apply_single_excitation(7, 5) + (0, 1, 3, 4, 5, 10) + """ + hole_index = self.index(hole) + # Where to put the particle so things are sorted? + particle_index = 0 + while (particle_index < self.popcnt() - 1) & (self[particle_index] < particle): + particle_index += 1 + # Upon return, either i == to self.popcnt() -1, or self[particle_index] > particle + # print(particle, self[hole_index + 1]) + if hole_index < particle_index: + if self[particle_index] > particle: + # Place partice before self[particle_index] in new tuple + t = ( + self[:hole_index] + + self[hole_index + 1 : particle_index] + + (particle,) + + self[particle_index:] + ) + else: + # Here, self[particle_index] < particle, so place particle after + t = ( + self[:hole_index] + + self[hole_index + 1 : particle_index + 1] + + (particle,) + + self[particle_index + 1 :] + ) + else: + t = ( + self[:particle_index] + + (particle,) + + self[particle_index:hole_index] + + self[hole_index + 1 :] + ) + return Spin_determinant_tuple(t) + + def apply_double_excitation( + self, h1: OrbitalIdx, p1: OrbitalIdx, h2: OrbitalIdx, p2: OrbitalIdx + ) -> Spin_determinant_tuple: + """Apply double h1 -> p1, h2 -> p2 excitation to instance of |Spin_determinant| + >>> Spin_determinant_tuple((0, 1)).apply_double_excitation(0, 2, 1, 3) + (2, 3) + >>> Spin_determinant_tuple((0, 4)).apply_double_excitation(4, 1, 0, 2) + (1, 2) + >>> Spin_determinant_tuple((0, 1, 3, 4, 7, 10)).apply_double_excitation(4, 8, 10, 2) + (0, 1, 2, 3, 7, 8) + """ + temp = self.apply_single_excitation(h1, p1) + return temp.apply_single_excitation(h2, p2) + + def exc_degree_spindet(self, right: Spin_determinant_tuple) -> int: + """Return excitation degree between two |Spin_determinant| + >>> Spin_determinant_tuple((0, 1)).exc_degree_spindet((1, 2)) + 1 + >>> Spin_determinant_tuple((0, 1)).exc_degree_spindet((2, 3)) + 2 + >>> Spin_determinant_tuple((0, 1)).exc_degree_spindet((0, 1)) + 0 + """ + return ((self ^ right).popcnt()) // 2 - def gen_all_connected_spindet(self, ed: int, n_orb: int) -> Iterator[Tuple[OrbitalIdx, ...]]: + def gen_all_connected_spindet(self, ed: int, n_orb: int) -> Iterator[Spin_determinant_tuple]: """Generate all connected spin determinants to self relative to a particular excitation degree :param n_orb: global parameter >>> sorted(Spin_determinant_tuple((0, 1)).gen_all_connected_spindet(1, 4)) @@ -143,6 +240,119 @@ def gen_all_connected_spindet(self, ed: int, n_orb: int) -> Iterator[Tuple[Orbit return [self ^ tuple((set(h) | set(p))) for h, p in l_hp_pairs] + # _ + # |_) |_ _. _ _ |_| _ | _ + # | | | (_| _> (/_ o | | (_) | (/_ + # _ / + # _. ._ _| |_) _. ._ _|_ o _ | _ + # (_| | | (_| | (_| | |_ | (_ | (/_ + + # Driver functions for computing phase, hole and particle between determinant pairs + + def get_holes(self, right: Spin_determinant_tuple) -> Spin_determinant_tuple: + """Get holes involved in excitation between two |Spin_determinant| + >>> Spin_determinant_tuple((0, 1)).get_holes((1, 2)) + (0,) + >>> Spin_determinant_tuple((0, 1)).get_holes((2, 3)) + (0, 1) + >>> Spin_determinant_tuple((0, 1)).get_holes((0, 1)) + () + """ + return (self ^ right) & self + + def get_particles(self, right: Spin_determinant_tuple) -> Spin_determinant_tuple: + """Get particles involved in excitation between two |Spin_determinant| + >>> Spin_determinant_tuple((0, 1)).get_particles((1, 2)) + (2,) + >>> Spin_determinant_tuple((0, 1)).get_particles((2, 3)) + (2, 3) + >>> Spin_determinant_tuple((0, 1)).get_particles((0, 1)) + () + """ + return (self ^ right) & right + + def single_phase( + self, + h: OrbitalIdx, + p: OrbitalIdx, + ): + """Function to compute phase for when I and J differ by exactly one orbital h <-> pd + + >>> Spin_determinant_tuple((0, 4, 6)).single_phase(4, 5) + 1 + >>> Spin_determinant_tuple((0, 1, 8)).single_phase(1, 17) + -1 + >>> Spin_determinant_tuple((0, 1, 4, 8)).single_phase(1, 17) + 1 + >>> Spin_determinant_tuple((0, 1, 4, 7, 8)).single_phase(1, 17) + -1 + """ + # Naive; compute phase for |Spin_determinant| pairs related by excitataion from h <-> p + j, k = min(h, p), max(h, p) + pmask = tuple((i for i in range(j + 1, k))) + parity = (self & pmask).popcnt() % 2 + return -1 if parity else 1 + + def double_phase(self, h1: OrbitalIdx, p1: OrbitalIdx, h2: OrbitalIdx, p2: OrbitalIdx): + """Function to compute phase for when I and J differ by exactly two orbitals h1, h2 <-> p1, p2 + Only for same spin double excitations + h1, h2 is occupied in self, p1, p2 is unoccupied + >>> Spin_determinant_tuple((0, 1, 2, 3, 4, 5, 6, 7, 8)).double_phase(2, 11, 3, 12) + 1 + >>> Spin_determinant_tuple((0, 1, 2, 3, 4, 5, 6, 7, 8)).double_phase(2, 11, 8, 17) + -1 + """ + # TODO: NOTE, elsewhere in qe.drivers code, will have to switch hole particle order in args! + # Compute phase. Loopless as in https://arxiv.org/abs/1311.6244 + phase = self.single_phase(h1, p1) * self.single_phase(h2, p2) + # if max(h1, p1) > min(h2, p2): + # return -phase + # else: + # return phase + if h2 < p1: + phase *= -1 + if p2 < h1: + phase *= -1 + return phase + + def single_exc(self, right: Spin_determinant_tuple) -> Tuple[int, OrbitalIdx, OrbitalIdx]: + """phase, hole, particle of when I and J differ by exactly one orbital + h is occupied only in I + p is occupied only in J + + >>> Spin_determinant_tuple((0, 4, 6)).single_exc((0, 5, 6)) + (1, 4, 5) + >>> Spin_determinant_tuple((0, 4, 6)).single_exc((0, 22, 6)) + (-1, 4, 22) + >>> Spin_determinant_tuple((0, 1, 8)).single_exc((0, 8, 17)) + (-1, 1, 17) + """ + # Get holes, particle in exc + + (h,) = self.get_holes(right) + (p,) = self.get_particles(right) + + return self.single_phase(h, p), h, p + + def double_exc( + self, right: Spin_determinant_tuple + ) -> Tuple[int, OrbitalIdx, OrbitalIdx, OrbitalIdx, OrbitalIdx]: + """phase, holes, particles of when I and J differ by exactly two orbitals + h1, h2 are occupied only in I + p1, p2 are occupied only in J + + >>> Spin_determinant_tuple((0, 1, 2, 3, 4, 5, 6, 7, 8)).double_exc((0, 1, 4, 5, 6, 7, 8, 11, 12)) + (1, 2, 3, 11, 12) + >>> Spin_determinant_tuple((0, 1, 2, 3, 4, 5, 6, 7, 8)).double_exc((0, 1, 3, 4, 5, 6, 7, 11, 17)) + (-1, 2, 8, 11, 17) + """ + # Holes + h1, h2 = self.get_holes(right) + # Particles + p1, p2 = self.get_particles(right) + + return self.double_phase(h1, p1, h2, p2), h1, h2, p1, p2 + # _ # |_) o _|_ _ _|_ ._ o ._ _ @@ -328,50 +538,6 @@ def popcnt(self) -> int: """Perform a `popcount'; number of bits set to True in self""" return self.bit_count() - def gen_all_connected_spindet(self, ed: int, n_orb: int) -> Iterator[Tuple[OrbitalIdx, ...]]: - """Generate all connected spin determinants to self relative to a particular excitation degree - :param n_orb: global parameter, used to pad bitstring with necessary 0s - - >>> for excited_sdet in sorted(Spin_determinant_bitstring(0b11).gen_all_connected_spindet(1, 4)): - ... bin(excited_sdet) - '0b101' - '0b110' - '0b1001' - '0b1010' - >>> for excited_sdet in sorted(Spin_determinant_bitstring(0b1).gen_all_connected_spindet(1, 4)): - ... bin(excited_sdet) - '0b10' - '0b100' - '0b1000' - >>> for excited_sdet in sorted(Spin_determinant_bitstring(0b11).gen_all_connected_spindet(2, 4)): - ... bin(excited_sdet) - '0b1100' - >>> sorted(Spin_determinant_bitstring(0b11).gen_all_connected_spindet(2, 2)) - [] - >>> for excited_sdet in sorted(Spin_determinant_bitstring(0b1000).gen_all_connected_spindet(1, 4)): - ... bin(excited_sdet) - '0b1' - '0b10' - '0b100' - """ - # Run through bitstring to create lists of particles and holes - holes = [] - particles = [] - # Some pre-processing; Create intermediate bitstring that is padded with zeros s.t len(inter) = n_orb - # If n_orb bit is already set -> No affect - inter = format(self, "#0" + f"{n_orb + 2}" + "b") - # One pass right -> left through reflected bitstring ([:-2] skips 'b0' in the flipped reflected bitstring) - for i, bit in enumerate(inter[::-1][:-2]): - # If ith bit is set, append to holes - if bit == "1": - holes.append(i) - # Else ith bit is not set, append to particles - else: - particles.append(i) - l_hp_pairs = product(combinations(tuple(holes), ed), combinations(tuple(particles), ed)) - - return [self ^ tuple(sorted(set(h) | set(p))) for h, p in l_hp_pairs] - # ______ _ _ _ # | _ \ | | (_) | | @@ -383,213 +549,154 @@ def gen_all_connected_spindet(self, ed: int, n_orb: int) -> Iterator[Tuple[Orbit # -class Determinant(tuple): - """Slater determinant: Product of 2 determinants. - One for $\alpha$ electrons and one for \beta electrons. - Abstract |Determinant| class; mimics behaviour of |NamedTuple|""" - - __slots__ = () - - _fields = ("alpha", "beta") - - def __new__(_cls, *args, **kwargs): - """Create new |Determinant| instance - If type of alpha, beta are |int| -> Return as is - ' ' |tuple| -> Convert to |Spin_determinant_tuple|, then return""" - - # Can either... - # Pass arguments in alpha, beta order as unnamed objects; e.g., Determinant(alpha_sdet, beta_sdet) - if len(args) > 0: - if ( - len(args) == 2 - ): # Most often, Determinants will be created via Determinant(alpha_sdet, beta_sdet) - alpha, beta = args - elif len(args) == 1: - # Arg passed might be tuple of sdets; Determinant(((alpha_sdet, beta_sdet)) - # Note: For some reason, mpi4py does this sometime when performing collectives - try: - # Unpack length-1 tuple - (_sdets,) = args - alpha, beta = _sdets - except: - raise TypeError(f"Expected 2 arguments, got {args}") - # Or... - # Pass arguments with keywords; e.g., Determinant(alpha=alpha_sdet, beta=beta_sdet) - elif len(kwargs) > 0: - try: - alpha = kwargs["alpha"] - beta = kwargs["beta"] - except: - raise KeyError(f"Expected keyword arguments for 'alpha', 'beta', got {kwargs}") - else: - raise TypeError(f"Expected two keyword arguments for 'alpha', 'beta', got {kwargs}") +class Determinant: + """Generic Slater determinant claass: Product of 2 determinants. + One for $\alpha$ electrons and one for \beta electrons.""" - # Once arguments are parsed, determine data representation - if isinstance(alpha, tuple): - _alpha = Spin_determinant_tuple(alpha) - elif isinstance(alpha, int): - _alpha = Spin_determinant_bitstring(alpha) - else: - raise TypeError( - f"Expected 'alpha' argument of type or , got {type(alpha)}" - ) - if isinstance(beta, tuple): - _beta = Spin_determinant_tuple(beta) - elif isinstance(beta, int): - _beta = Spin_determinant_bitstring(beta) + def __init__(self, alpha, beta, representation="tuple"): + self.flag = representation + if self.flag == "tuple": + self.alpha = Spin_determinant_tuple(alpha) + self.beta = Spin_determinant_tuple(beta) else: - raise TypeError( - f"Expected 'beta' argument of type or , got {type(beta)}" - ) - return tuple.__new__(_cls, (_alpha, _beta)) - - @classmethod - def _make(cls, iterable): - """Make a new Determinant object from a sequence or iterable - Used in batch `apply_excitation' calls""" - result = tuple.__new__(cls, iterable) - if len(result) != 2: - raise TypeError(f"Expected 2 arguments, got {len(result)}") - return result + raise NotImplementedError def __repr__(self): """Return a nicely formatted representation string""" - return "Determinant(alpha=%r, beta=%r)" % self + return "Determinant(alpha=%r, beta=%r)" % (self.alpha, self.beta) - @property - def alpha(self): - """Return alpha spin determinant""" - return self[0] + def __eq__(self, right: Determinant): + """Comparison operator""" + if (self.alpha == right.alpha) & (self.beta == right.beta): + return True + else: + return False - @property - def beta(self): - """Return beta spin determinant""" - return self[1] + def __hash__(self): + """Hash custom Det; Concatenate tuples, in [\alpha, \beta] order. Return hash of this 2 * N_a * N_b tuple""" + return hash((self.alpha, self.beta)) - def convert_repr(self, Norb=None): - """Conver to bitstring (tuple) representation of |Determinant| if self is tuple (bitstring) - >>> Determinant((0, 2), (1,)).convert_repr(4) - Determinant(alpha=5, beta=2) - >>> Determinant((0, 2), ()).convert_repr(4) - Determinant(alpha=5, beta=0) - - >>> Determinant(0b0101, 0b0001).convert_repr() - Determinant(alpha=(0, 2), beta=(0,)) - >>> Determinant(0b0101, 0b0000).convert_repr() - Determinant(alpha=(0, 2), beta=()) - """ - # Each spin determinant class has member `convert_repr` function - return Determinant(self.alpha.convert_repr(Norb), self.beta.convert_repr(Norb)) + def __iter__(self): + """Unpack determinant""" + return iter((self.alpha, self.beta)) # _ # |_ _ o _|_ _. _|_ o _ ._ # |_ >< (_ | |_ (_| |_ | (_) | | # - def apply_excitation( - self, - alpha_exc: Tuple[Tuple[OrbitalIdx, ...], Tuple[OrbitalIdx, ...]], - beta_exc: Tuple[Tuple[OrbitalIdx, ...], Tuple[OrbitalIdx, ...]], - ) -> NamedTuple: - """Apply excitation to self, produce new |Determinant| - Each type |Determinant_tuple| and |Determinant_bitstring| has own implementation of `apply_excitation_to_spindet` based on type + def apply_single_excitation(self, h: OrbitalIdx, p: OrbitalIdx, spin: str) -> Determinant: + """Apply single excitation to self, return new abstract |Determinant| class + Each fundamental type has own implementation of spindet excitaations Inputs - :param `alpha_exc` (`beta_exc`): Specifies holes, particles involved in excitation of alpha (beta) |Spin_determinant| - - If either argument is empty (), no excitation is applied - >>> Determinant((0, 1), (0, 1)).apply_excitation(((1,), (2,)), ((1,), (2,))) - Determinant(alpha=(0, 2), beta=(0, 2)) - >>> Determinant((0, 1), (0, 1)).apply_excitation(((0, 1), (2, 3)), ((0, 1), (3, 4))) - Determinant(alpha=(2, 3), beta=(3, 4)) - >>> Determinant((0, 1), (0, 1)).apply_excitation(((), ()), ((), ())) - Determinant(alpha=(0, 1), beta=(0, 1)) - - >>> Determinant(0b11, 0b11).apply_excitation(((1,), (2,)), ((1,), (2,))) - Determinant(alpha=5, beta=5) - >>> Determinant(0b11, 0b11).apply_excitation(((0, 1), (2, 3)), ((0, 1), (3, 4))) - Determinant(alpha=12, beta=24) - >>> Determinant(0b11, 0b11).apply_excitation(((), ()), ((), ())) - Determinant(alpha=3, beta=3) - >>> Determinant(0b11, 0b11).apply_excitation(((1,), (2,)), ((), ())) - Determinant(alpha=5, beta=3) + :param h, p: Specifies hole, particle involved in excitation of |Spin_determinant| + :param spin: Spin-type of excitation, "alpha" or "beta" + + >>> Determinant((0, 1), (0, 1), "tuple").apply_single_excitation(1, 2, "alpha") + Determinant(alpha=(0, 2), beta=(0, 1)) + >>> Determinant((0, 1, 2), (0, 1), "tuple").apply_single_excitation(0, 4, "alpha") + Determinant(alpha=(1, 2, 4), beta=(0, 1)) + >>> Determinant((0, 1), (0, 1), "tuple").apply_single_excitation(0, 2, "beta") + Determinant(alpha=(0, 1), beta=(1, 2)) """ - # Unpack alpha, beta holes - lh_a, lp_a = alpha_exc - lh_b, lp_b = beta_exc + if spin == "alpha": + return Determinant(self.alpha.apply_single_excitation(h, p), self.beta, self.flag) + elif spin == "beta": + return Determinant(self.alpha, self.beta.apply_single_excitation(h, p), self.flag) + else: + raise NotImplementedError + + def apply_same_spin_double_excitation( + self, h1: OrbitalIdx, p1: OrbitalIdx, h2: OrbitalIdx, p2: OrbitalIdx, spin: str + ) -> Determinant: + """Apply double excitation to self, same-spin, return new abstract |Determinant| class + Each fundamental type has own implementation of spindet excitaations + Inputs + :param h1, p1, h2, p2: Specifies holes, particles involved in excitation of |Spin_determinant| + :param spin: Spin-type of excitation, "alpha" or "beta" + + + >>> Determinant((0, 1), (0, 1), "tuple").apply_same_spin_double_excitation(0, 3, 1, 2, "alpha") + Determinant(alpha=(2, 3), beta=(0, 1)) + >>> Determinant((0, 1, 3), (0, 1), "tuple").apply_same_spin_double_excitation(3, 2, 1, 4, "alpha") + Determinant(alpha=(0, 2, 4), beta=(0, 1)) + >>> Determinant((0, 1), (0, 1), "tuple").apply_same_spin_double_excitation(0, 2, 1, 3, "beta") + Determinant(alpha=(0, 1), beta=(2, 3)) + """ + if spin == "alpha": + return Determinant( + self.alpha.apply_double_excitation(h1, p1, h2, p2), self.beta, self.flag + ) + elif spin == "beta": + return Determinant( + self.alpha, self.beta.apply_double_excitation(h1, p1, h2, p2), self.flag + ) + else: + raise NotImplementedError + + def apply_opposite_spin_double_excitation( + self, h1: OrbitalIdx, p1: OrbitalIdx, h2: OrbitalIdx, p2: OrbitalIdx + ) -> Determinant: + """Apply double excitation to self, opposite-spin, return new abstract |Determinant| class + Each fundamental type has own implementation of spindet excitaations + Inputs + :param h1, p1, h2, p2: Specifies holes, particles involved in excitation. Assumed that h1, p1 (h2, p2) -> alpha exc (beta exc) - excited_sdet_a = self.alpha ^ (tuple(sorted(set(lh_a) | set(lp_a)))) - excited_sdet_b = self.beta ^ (tuple(sorted(set(lh_b) | set(lp_b)))) - return Determinant(excited_sdet_a, excited_sdet_b) + >>> Determinant((0, 1), (0, 1), "tuple").apply_opposite_spin_double_excitation(0, 3, 1, 2) + Determinant(alpha=(1, 3), beta=(0, 2)) + >>> Determinant((0, 1, 3), (0, 1), "tuple").apply_opposite_spin_double_excitation(3, 2, 1, 4) + Determinant(alpha=(0, 1, 2), beta=(0, 4)) + """ + return Determinant( + self.alpha.apply_single_excitation(h1, p1), + self.beta.apply_single_excitation(h2, p2), + self.flag, + ) - def exc_degree(self, det_J: NamedTuple) -> Tuple[int, int]: + def exc_degree(self, right: Determinant) -> int: """Compute excitation degree; the number of orbitals which differ between two |Determinants| self, det_J - >>> Determinant((0, 1), (0, 1)).exc_degree(Determinant(alpha=(0, 2), beta=(4, 6))) - (1, 2) - >>> Determinant((0, 1), (0, 1)).exc_degree(Determinant(alpha=(0, 1), beta=(4, 6))) - (0, 2) - >>> Determinant(0b11, 0b11).exc_degree(Determinant(0b101, 0b101000)) + >>> Determinant((0, 1), (0, 1), "tuple").exc_degree(Determinant((0, 2), (4, 6), "tuple")) (1, 2) - >>> Determinant(0b11, 0b11).exc_degree(Determinant(0b11, 0b101000)) + >>> Determinant((0, 1), (0, 1), "tuple").exc_degree(Determinant((0, 1), (4, 6), "tuple")) (0, 2) """ - ed_up = (self.alpha ^ det_J.alpha).popcnt() // 2 - ed_dn = (self.beta ^ det_J.beta).popcnt() // 2 + ed_up = (self.alpha ^ right.alpha).popcnt() // 2 + ed_dn = (self.beta ^ right.beta).popcnt() // 2 return (ed_up, ed_dn) - def is_connected(self, det_j) -> Tuple[int, int]: + def is_connected(self, det_j: Determinant) -> Tuple[int, int]: """Compute the excitation degree, the number of orbitals which differ between the two determinants. Return bool; `Is det_j (singley or doubley) connected to instance of self? - >>> Determinant((0, 1), (0, 1)).is_connected(Determinant((0, 1), (0, 2))) + >>> Determinant((0, 1), (0, 1), "tuple").is_connected(Determinant((0, 1), (0, 2), "tuple")) True - >>> Determinant((0, 1), (0, 1)).is_connected(Determinant((0, 2), (0, 2))) + >>> Determinant((0, 1), (0, 1), "tuple").is_connected(Determinant((0, 2), (0, 2), "tuple")) True - >>> Determinant((0, 1), (0, 1)).is_connected(Determinant((2, 3), (0, 1))) + >>> Determinant((0, 1), (0, 1), "tuple").is_connected(Determinant((2, 3), (0, 1), "tuple")) True - >>> Determinant((0, 1), (0, 1)).is_connected(Determinant((2, 3), (0, 2))) + >>> Determinant((0, 1), (0, 1), "tuple").is_connected(Determinant((2, 3), (0, 2), "tuple")) False - >>> Determinant((0, 1), (0, 1)).is_connected(Determinant((0, 1), (0, 1))) - False - - >>> Determinant(0b11, 0b11).is_connected(Determinant(0b11, 0b101)) - True - >>> Determinant(0b11, 0b11).is_connected(Determinant(0b101, 0b101)) - True - >>> Determinant(0b11, 0b11).is_connected(Determinant(0b1100, 0b11)) - True - >>> Determinant(0b11, 0b11).is_connected(Determinant(0b1100, 0b101)) - False - >>> Determinant(0b11, 0b11).is_connected(Determinant(0b11, 0b11)) + >>> Determinant((0, 1), (0, 1), "tuple").is_connected(Determinant((0, 1), (0, 1), "tuple")) False """ return sum(self.exc_degree(det_j)) in [1, 2] - def gen_all_connected_det(self, n_orb: int) -> Iterator[NamedTuple]: + def gen_all_connected_det(self, n_orb: int) -> Iterator[Determinant]: """Generate all determinants that are singly or doubly connected to self :param n_orb: global parameter, needed to cap possible excitations - >>> sorted(Determinant((0, 1), (0,)).gen_all_connected_det(3)) - [Determinant(alpha=(0, 1), beta=(1,)), - Determinant(alpha=(0, 1), beta=(2,)), - Determinant(alpha=(0, 2), beta=(0,)), - Determinant(alpha=(0, 2), beta=(1,)), - Determinant(alpha=(0, 2), beta=(2,)), - Determinant(alpha=(1, 2), beta=(0,)), - Determinant(alpha=(1, 2), beta=(1,)), - Determinant(alpha=(1, 2), beta=(2,))] - - >>> for excited_det in sorted(Determinant(0b11, 0b1).gen_all_connected_det(3)): - ... [bin(excited_det.alpha), bin(excited_det.beta)] - ['0b11', '0b10'] - ['0b11', '0b100'] - ['0b101', '0b1'] - ['0b101', '0b10'] - ['0b101', '0b100'] - ['0b110', '0b1'] - ['0b110', '0b10'] - ['0b110', '0b100'] + >>> [d for d in (Determinant((0, 1), (0,), "tuple").gen_all_connected_det(3))] + [Determinant(alpha=(0, 2), beta=(0,)), + Determinant(alpha=(1, 2), beta=(0,)), + Determinant(alpha=(0, 1), beta=(1,)), + Determinant(alpha=(0, 1), beta=(2,)), + Determinant(alpha=(0, 2), beta=(1,)), + Determinant(alpha=(0, 2), beta=(2,)), + Determinant(alpha=(1, 2), beta=(1,)), + Determinant(alpha=(1, 2), beta=(2,))] """ + # Generate all singles from constituent alpha and beta spin determinants # Then, opposite-spin and same-spin doubles @@ -598,35 +705,44 @@ def gen_all_connected_det(self, n_orb: int) -> Iterator[NamedTuple]: l_double_aa = self.alpha.gen_all_connected_spindet(2, n_orb) # Singles and doubles; alpha spin - exc_a = (Determinant(det_alpha, self.beta) for det_alpha in chain(l_single_a, l_double_aa)) + exc_a = ( + Determinant(det_alpha, self.beta, self.flag) + for det_alpha in chain(l_single_a, l_double_aa) + ) l_single_b = set(self.beta.gen_all_connected_spindet(1, n_orb)) l_double_bb = self.beta.gen_all_connected_spindet(2, n_orb) # Singles and doubles; beta spin - exc_b = (Determinant(self.alpha, det_beta) for det_beta in chain(l_single_b, l_double_bb)) + exc_b = ( + Determinant(self.alpha, det_beta, self.flag) + for det_beta in chain(l_single_b, l_double_bb) + ) l_double_ab = product(l_single_a, l_single_b) # Doubles; opposite-spin - exc_ab = (Determinant(det_alpha, det_beta) for det_alpha, det_beta in l_double_ab) + exc_ab = ( + Determinant(det_alpha, det_beta, self.flag) for det_alpha, det_beta in l_double_ab + ) return chain(exc_a, exc_b, exc_ab) def triplet_constrained_single_excitations_from_det( self, constraint: Tuple[OrbitalIdx, ...], n_orb: int, spin="alpha" - ) -> Iterator[NamedTuple]: - """Called by inherited classes; Generate singlet excitations from constraint""" + ) -> Iterator[Determinant]: + """Called by inherited classes; Generate singlet excitations from constraint + :param spin: Refers to the spin type of `constraint`""" ha, pa, hb, pb = self.get_holes_particles_for_constrained_singles(constraint, n_orb, spin) # Excitations of argument `spin` for h, p in product(ha, pa): if spin == "alpha": # Then, det_a is alpha spindet - excited_det = self.apply_excitation(((h,), (p,)), ((), ())) + excited_det = self.apply_single_excitation(h, p, "alpha") else: # det_a is beta spindet - excited_det = self.apply_excitation(((), ()), ((h,), (p,))) + excited_det = self.apply_single_excitation(h, p, "beta") assert getattr(excited_det, spin)[-3:] == constraint yield excited_det @@ -634,18 +750,19 @@ def triplet_constrained_single_excitations_from_det( for h, p in product(hb, pb): if spin == "alpha": # Then, det_b is beta spindet - excited_det = self.apply_excitation(((), ()), ((h,), (p,))) + excited_det = self.apply_single_excitation(h, p, "beta") else: # det_b is alpha spindet - excited_det = self.apply_excitation(((h,), (p,)), ((), ())) + excited_det = self.apply_single_excitation(h, p, "alpha") # TODO: Assertion for bitstring? Shouldn't need though assert getattr(excited_det, spin)[-3:] == constraint yield excited_det def triplet_constrained_double_excitations_from_det( self, constraint: Tuple[OrbitalIdx, ...], n_orb: int, spin="alpha" - ) -> Iterator[NamedTuple]: - """Called by inherited classes; Generate singlet excitations from constraint""" + ) -> Iterator[Determinant]: + """Called by inherited classes; Generate singlet excitations from constraint + :param spin: Refers to the spin type of `constraint`""" # |Determinant_tuple| and |Determinant_bitstring| each have this method haa, paa, hbb, pbb, hab, pab = self.get_holes_particles_for_constrained_doubles( @@ -653,24 +770,24 @@ def triplet_constrained_double_excitations_from_det( ) # Excitations of argument `spin` # Same-spin doubles, for argument `spin` - for holes, particles in product(haa, paa): + for (h1, h2), (p1, p2) in product(haa, paa): if spin == "alpha": # Then, det_a is alpha spindet - excited_det = self.apply_excitation((holes, particles), ((), ())) + excited_det = self.apply_same_spin_double_excitation(h1, p1, h2, p2, "alpha") else: # det_a is beta spindet - excited_det = self.apply_excitation(((), ()), (holes, particles)) + excited_det = self.apply_same_spin_double_excitation(h1, p1, h2, p2, "beta") assert getattr(excited_det, spin)[-3:] == constraint yield excited_det # Same-spin doubles, for opposite-spin to `spin` - for holes, particles in product(hbb, pbb): + for (h1, h2), (p1, p2) in product(hbb, pbb): if spin == "alpha": # Then, det_b is beta spindet - excited_det = self.apply_excitation(((), ()), (holes, particles)) + excited_det = self.apply_same_spin_double_excitation(h1, p1, h2, p2, "beta") else: # det_b is alpha spindet - excited_det = self.apply_excitation((holes, particles), ((), ())) + excited_det = self.apply_same_spin_double_excitation(h1, p1, h2, p2, "alpha") assert getattr(excited_det, spin)[-3:] == constraint yield excited_det @@ -680,10 +797,10 @@ def triplet_constrained_double_excitations_from_det( pa, pb = particles if spin == "alpha": # det_a is alpha, det_b is beta - excited_det = self.apply_excitation(((ha,), (pa,)), ((hb,), (pb,))) + excited_det = self.apply_opposite_spin_double_excitation(ha, pa, hb, pb) else: # det_a is beta, det_b is beta - excited_det = self.apply_excitation(((hb,), (pb,)), ((ha,), (pa,))) + excited_det = self.apply_opposite_spin_double_excitation(hb, pb, ha, pa) assert getattr(excited_det, spin)[-3:] == constraint yield excited_det @@ -724,40 +841,45 @@ def get_holes_particles_for_constrained_singles( constraint_orbitals_occupied = det_a & constraint # Which `higher-order` (spin) orbitals o >= a1 that are not {a1, a2, a3} are occupied in |D_a>? (If any) # TODO: Different from Tubman paper, which has an error if I reada it correctly - nonconstrained_orbitals_occupied = (det_a & B) - constraint + # Equivalent to `det_a & B - constraint` in set operations + nonconstrained_orbitals_occupied = (det_a & B).get_holes(constraint) # If no double excitation of |D> will produce |J> satisfying constraint - if len(constraint_orbitals_occupied) == 1 or len(nonconstrained_orbitals_occupied) > 1: + if ( + constraint_orbitals_occupied.popcnt() == 1 + or nonconstrained_orbitals_occupied.popcnt() > 1 + ): # No single excitations generated by the inputted |Determinant|: {det} satisfy given constraint: {constraint} # These are empty lists return (ha, pa, hb, pb) # Create list of possible `particles` s.to constraint - if len(constraint_orbitals_occupied) == 2: + if constraint_orbitals_occupied.popcnt() == 2: # (Two constraint orbitals occupied) e.g., a1, a2 \in |D_a> -> A single (a) x_a \in ha to a_unocc is necessary to satisfy the constraint # A single (b) will still not satisfy constraint - (a_unocc,) = ( - (det_a | constraint) - (det_a & constraint) - ) & constraint # The 1 unoccupied constraint orbital + # Operation equivalent to `(det_a ^ constraint) & constraint ` + (a_unocc,) = det_a.get_particles(constraint) # The 1 unoccupied constraint orbital pa.append(a_unocc) - elif len(constraint_orbitals_occupied) == 3: + if constraint_orbitals_occupied.popcnt() == 3: # a1, a2, a3 \in |D_a> -> |D> satisfies constraint! -> # Any single x_a \in ha to w_a where w_a < a1 will satisfy constraint - det_unocc_a_orbs = all_orbs - det_a + # Operation below is equivalent to all_orbs - det_al; get orbs that are unoccupied in alpha spindet + det_unocc_a_orbs = det_a.get_particles(all_orbs) for w_a in takewhile(lambda x: x < a1, det_unocc_a_orbs): pa.append(w_a) # Any single x_b \in hb to w_b - det_unocc_b_orbs = all_orbs - det_b + # Operation below is equivalent to all_orbs - det_b; get orbs that are unoccupied in beta spindet + det_unocc_b_orbs = det_b.get_particles(all_orbs) for w_b in det_unocc_b_orbs: pb.append(w_b) # Create list of possible `holes` s.to constraint - if len(nonconstrained_orbitals_occupied) == 1: + if nonconstrained_orbitals_occupied.popcnt() == 1: # x_a > a1 \in |D_a> with x_a \not\in {a1, a2, a3} -> A single (a) x_a to w_a \in pa is necessary to satisfy constraint # A single (b) will not satisfy (x_a,) = nonconstrained_orbitals_occupied # Has length 1; unpack ha.append(x_a) - elif len(nonconstrained_orbitals_occupied) == 0: + elif nonconstrained_orbitals_occupied.popcnt() == 0: # No `higher` orbitals \not\in {a1, a2, a3} occupied in |D> -> # A single (a) x_a to w_a \in pa, where x_a < a1 (so as not to ruin constraint) for x_a in takewhile(lambda x: x < a1, det_a): @@ -818,31 +940,36 @@ def get_holes_particles_for_constrained_doubles( constraint_orbitals_occupied = det_a & constraint # Which `higher-order`(spin) orbitals o >= a1 that are not {a1, a2, a3} are occupied in |D>? (If any) # TODO: Different from Tubman paper, which has an error if I read it correctly... - nonconstrained_orbitals_occupied = (det_a & B) - constraint + # Equivalent to `det_a & B - constraint` in set operations + nonconstrained_orbitals_occupied = (det_a & B).get_holes(constraint) # If this -> no double excitation of |D> will produce |J> satisfying constraint |T> - if len(constraint_orbitals_occupied) == 0 or len(nonconstrained_orbitals_occupied) > 2: + if ( + constraint_orbitals_occupied.popcnt() == 0 + or nonconstrained_orbitals_occupied.popcnt() > 2 + ): # No double excitations generated by the inputted |Determinant|: {det} satisfy given constraint: {constraint} # These are empty lists return (haa, paa, hbb, pbb, hab, pab) # Create list of possible `particles` s.to constraint - if len(constraint_orbitals_occupied) == 1: + if constraint_orbitals_occupied.popcnt() == 1: # (One constraint orbital occupied) e.g., a1 \in |D_a> -> A same-spin (aa) double to (x_a, y_a) \in haa to (a2, a3) is necessary # No same-spin (bb) or opposite-spin (ab) excitations will satisfy constraint # New particles -> a2, a3 - a_unocc_1, a_unocc_2 = ((det_a | constraint) - (det_a & constraint)) & ( + # Operation equivalent to `(det_a ^ constraint) & constraint ` + a_unocc_1, a_unocc_2 = det_a.get_particles( constraint ) # This set will have length 2; unpack paa.append((a_unocc_1, a_unocc_2)) - elif len(constraint_orbitals_occupied) == 2: + elif constraint_orbitals_occupied.popcnt() == 2: # (Two constraint orbitals occupied) e.g., a1, a2 \in |D_a> -> # A same-spin (aa) double (x_a, y_a) \in haa to (z_a, a_unocc), where z_a\not\in|D_a>, and z_a < a1 (if excited to z_a > a1, constraint ruined!) - (a_unocc,) = ( - (det_a | constraint) - (det_a & constraint) - ) & constraint # The 1 unoccupied constraint orbital - det_unocc_a_orbs = all_orbs - det_a # Unocc orbitals in |D_a> + # Operation equivalent to `(det_a ^ constraint) & constraint ` + (a_unocc,) = det_a.get_particles(constraint) # The 1 unoccupied constraint orbital + # Operation below is equivalent to `all_orbs - det_a` + det_unocc_a_orbs = det_a.get_particles(all_orbs) # Unocc orbitals in |D_a> for z_a in takewhile(lambda x: x < a1, det_unocc_a_orbs): # z < a_unocc trivially, no need to check they are distinct paa.append((z_a, a_unocc)) @@ -852,15 +979,17 @@ def get_holes_particles_for_constrained_doubles( for z_b in det_unocc_b_orbs: pab.append((a_unocc, z_b)) - elif len(constraint_orbitals_occupied) == 3: + elif constraint_orbitals_occupied.popcnt() == 3: # a1, a2, a3 \in |D_a> -> |D> satisfies constraint! -> # Any same-spin (aa) double (x_a, y_a) \in haa to (w_a, z_a), where w_a, z_a < a1 - det_unocc_a_orbs = all_orbs - det_a + # Operations below is equivalent to `all_orbs - det_a` + det_unocc_a_orbs = det_a.get_particles(all_orbs) for w_a in takewhile(lambda x: x < a1, det_unocc_a_orbs): for z_a in takewhile(lambda z: z < w_a, det_unocc_a_orbs): paa.append((w_a, z_a)) # Any same-spin (bb) double (x_b, y_b) \in hbb to (w_b, z_b) - det_unocc_b_orbs = all_orbs - det_b # Unocc orbitals in |D_a> + # Operations below is equivalent to `all_orbs - det_b` + det_unocc_b_orbs = det_b.get_particles(all_orbs) # Unocc orbitals in |D_b> for w_b in det_unocc_b_orbs: for z_b in takewhile(lambda x: x < w_b, det_unocc_b_orbs): pbb.append((w_b, z_b)) @@ -870,13 +999,13 @@ def get_holes_particles_for_constrained_doubles( pab.append((w_a, z_b)) # Create list of possible `holes` s.to constraint - if len(nonconstrained_orbitals_occupied) == 2: + if nonconstrained_orbitals_occupied.popcnt() == 2: # x_a, y_a \in |D_a> with x_a, y_a > a1 and \not\in {a1, a2, a3} -> A same-spin (aa) double (x_a, y_a) to (w_a, z_a) \in paa is necessary # No same-spin (bb) or opposite-spin (ab) excitations will satisfy constraint # New holes -> x, y x_a, y_a = nonconstrained_orbitals_occupied # This set will have length 2; unpack haa.append((x_a, y_a)) - elif len(nonconstrained_orbitals_occupied) == 1: + elif nonconstrained_orbitals_occupied.popcnt() == 1: # x_a > a1 \in |D_a> with x_a \not\in {a1, a2, a3} -> # A same-spin (aa) double (x_a, y_a) to (w_a, z_a) \in paa, where y_a < a1 (exciting y_a < a1 doesn't remove constraint) (x_a,) = nonconstrained_orbitals_occupied # Has length 1; unpack @@ -888,7 +1017,7 @@ def get_holes_particles_for_constrained_doubles( for y_b in det_b: hab.append((x_a, y_b)) - elif len(nonconstrained_orbitals_occupied) == 0: + elif nonconstrained_orbitals_occupied.popcnt() == 0: # No `higher` orbitals \not\in {a1, a2, a3} occupied in |D> -> # A same-spin (aa) double (x_a, y_a) to (w_a, z_a) \in paa, where x_a, y_a < a1 for x_a in takewhile(lambda x: x < a1, det_a): @@ -912,116 +1041,7 @@ def get_holes_particles_for_constrained_doubles( # _. ._ _| |_) _. ._ _|_ o _ | _ # (_| | | (_| | (_| | |_ | (_ | (/_ - # Driver functions for computing phase, hole and particle between determinant pairs - # TODO: These are ONLY implemented for |Spin_determinant_tuple| at the moment; - # So, might just keep as is - - def single_phase( - self, - h: OrbitalIdx, - p: OrbitalIdx, - spin: str, - ): - """Function to compute phase for when I and J differ by exactly one orbital h <-> p - h is occupied in sdet = getattr(self, spin), p is unoccupied - - >>> Determinant((0, 4, 6), ()).single_phase(4, 5, "alpha") - 1 - >>> Determinant((), (0, 1, 8)).single_phase(1, 17, "beta") - -1 - >>> Determinant((0, 1, 4, 8), ()).single_phase(1, 17, "alpha") - 1 - >>> Determinant((0, 1, 4, 7, 8), ()).single_phase(1, 17, "alpha") - -1 - - >>> Determinant(0b1010001, 0b0).single_phase(4, 6, "alpha") - 1 - >>> Determinant(0b0, 0b100000011).single_phase(1, 17, "beta") - -1 - >>> Determinant(0b100010011, 0b0).single_phase(1, 17, "alpha") - 1 - >>> Determinant(0b110010011, 0b0).single_phase(1, 17, "alpha") - -1 - """ - # Naive; compute phase for |Spin_determinant| pairs related by excitataion from h <-> p - sdet = getattr(self, spin) - j, k = min(h, p), max(h, p) - - if isinstance(sdet, tuple): - pmask = tuple((i for i in range(j + 1, k))) - elif isinstance(sdet, int): - # Strings immutable -> Work with lists for ON assignment, convert to |str| -> |int| when done - pmask = ["0", "b"] - pmask.extend(["0"] * k) - for i in range(j + 1, k): - pmask[-(i + 1)] = "1" - pmask = int(("".join(pmask)), 2) - - parity = (sdet & pmask).popcnt() % 2 - return -1 if parity else 1 - - def double_phase(self, holes: Tuple[OrbitalIdx, ...], particles: Tuple[OrbitalIdx, ...], spin): - """Function to compute phase for when I and J differ by exactly two orbitals h1, h2 <-> p1, p2 - Only for same spin double excitations - h1, h2 is occupied in sdet = getattr(self, spin), p1, p2 is unoccupied - >>> Determinant((0, 1, 2, 3, 4, 5, 6, 7, 8), ()).double_phase((2, 3), (11, 12), "alpha") - 1 - >>> Determinant((0, 1, 2, 3, 4, 5, 6, 7, 8), ()).double_phase((2, 8), (11, 17), "alpha") - -1 - """ - # Compute phase. Loopless as in https://arxiv.org/abs/1311.6244 - h1, h2 = holes - p1, p2 = particles - phase = self.single_phase(h1, p1, spin) * self.single_phase(h2, p2, spin) - # if max(h1, p1) > min(h2, p2): - # return -phase - # else: - # return phase - if h2 < p1: - phase *= -1 - if p2 < h1: - phase *= -1 - return phase - - def single_exc(self, sdet_j, spin: str) -> Tuple[int, OrbitalIdx, OrbitalIdx]: - """phase, hole, particle of when I and J differ by exactly one orbital - h is occupied only in I - p is occupied only in J - - >>> Determinant((0, 4, 6), ()).single_exc((0, 5, 6), "alpha") - (1, 4, 5) - >>> Determinant((0, 4, 6), ()).single_exc((0, 22, 6), "alpha") - (-1, 4, 22) - >>> Determinant((), (0, 1, 8)).single_exc((0, 8, 17), "beta") - (-1, 1, 17) - """ - # Get holes, particle in exc - sdet_i = getattr(self, spin) - (h,) = sdet_i - sdet_j - (p,) = sdet_j - sdet_i - - return self.single_phase(h, p, spin), h, p - - def double_exc( - self, sdet_j: Tuple[OrbitalIdx, ...], spin: str - ) -> Tuple[int, OrbitalIdx, OrbitalIdx, OrbitalIdx, OrbitalIdx]: - """phase, holes, particles of when I and J differ by exactly two orbitals - h1, h2 are occupied only in I - p1, p2 are occupied only in J - - >>> Determinant((0, 1, 2, 3, 4, 5, 6, 7, 8), ()).double_exc((0, 1, 4, 5, 6, 7, 8, 11, 12), "alpha") - (1, 2, 3, 11, 12) - >>> Determinant((), (0, 1, 2, 3, 4, 5, 6, 7, 8)).double_exc((0, 1, 3, 4, 5, 6, 7, 11, 17), "beta") - (-1, 2, 8, 11, 17) - """ - sdet_i = getattr(self, spin) - # Holes - h1, h2 = sorted(sdet_i - sdet_j) - # Particles - p1, p2 = sorted(sdet_j - sdet_i) - - return self.double_phase((h1, h2), (p1, p2), spin), h1, h2, p1, p2 - + # These are used for testing; all other phase, hole, particle functions are implemented on the spin determinant level @staticmethod def single_exc_no_phase( sdet_i: Tuple[OrbitalIdx, ...], sdet_j: Tuple[OrbitalIdx, ...] diff --git a/qe/io.py b/qe/io.py index 4b8a7ec..a0b25b4 100644 --- a/qe/io.py +++ b/qe/io.py @@ -159,12 +159,15 @@ def grouper(iterable, n): Determinant( tuple(decode_det(det_i, det_representation)), tuple(decode_det(det_j, det_representation)), + "tuple", ) ) elif det_representation == "bitstring": alpha_str = ["0", "b"] + [bit for bit in decode_det(det_i, det_representation)][::-1] beta_str = ["0", "b"] + [bit for bit in decode_det(det_j, det_representation)][::-1] - det.append(Determinant(int(("".join(alpha_str)), 2), int(("".join(beta_str)), 2))) + det.append( + Determinant(int(("".join(alpha_str)), 2), int(("".join(beta_str)), 2), "bitset") + ) else: raise NotImplementedError diff --git a/tests/test_everything_all_at_once.py b/tests/test_everything_all_at_once.py index a466d1e..8ec8ec1 100755 --- a/tests/test_everything_all_at_once.py +++ b/tests/test_everything_all_at_once.py @@ -268,8 +268,8 @@ def psi_and_integral(self): # 4 Electron in 4 Orbital # I'm stupid so let's do the product n_orb = 4 - psi = [Determinant((0, 1), (0, 1))] - for det in Determinant((0, 1), (0, 1)).gen_all_connected_det(n_orb): + psi = [Determinant((0, 1), (0, 1), "tuple")] + for det in Determinant((0, 1), (0, 1), "tuple").gen_all_connected_det(n_orb): psi.append(det) d_two_e_integral = {} for i, j, k, l in product(range(n_orb), repeat=4): @@ -281,7 +281,10 @@ def psi_and_integral_PT2(self): # minimal psi_and_integral, psi_i != psi_j # Na = 3, Nb = 3, Norb = 6 to account for triplet constraints n_orb = 6 - psi_i = [Determinant((0, 1, 2), (0, 1, 2)), Determinant((1, 2, 3), (1, 2, 3))] + psi_i = [ + Determinant((0, 1, 2), (0, 1, 2), "tuple"), + Determinant((1, 2, 3), (1, 2, 3), "tuple"), + ] psi_j = [] for i, det in enumerate(psi_i): for det_connected in det.gen_all_connected_det(n_orb): @@ -340,8 +343,8 @@ class Test_Constrained_Excitation(Timing, unittest.TestCase): def psi_and_norb_2det(self): # Do 5 e, 10 orb return 10, [ - Determinant((0, 1, 2, 3, 4), (0, 1, 2, 3, 4)), - Determinant((1, 2, 3, 4, 5), (1, 2, 3, 4, 5)), + Determinant((0, 1, 2, 3, 4), (0, 1, 2, 3, 4), "tuple"), + Determinant((1, 2, 3, 4, 5), (1, 2, 3, 4, 5), "tuple"), ] @cached_property @@ -389,16 +392,10 @@ def connected_by_constraint(self): ) for det in self.psi_connected_2det: # What triplet constraint does this det satisfy? - con = self.check_constraint(det, "alpha") + con = check_constraint(det, "alpha") d[con].append(det) return d - def check_constraint(self, det: Determinant, spin="alpha"): - # Give me a determinant. What constraint does it satisfy? (What are three most highly occupied alpha spin orbitas) - spindet = getattr(det, spin) - # Return constraint as |Spin_determinant| - return spindet[-3:] - def test_constrained_excitations(self, spin="alpha"): # 1. For each constraint, pass through wf n_orb, psi = self.psi_and_norb_2det @@ -429,12 +426,23 @@ def test_constrained_excitations(self, spin="alpha"): # This is for ONE determinant.. So there should be no duplicates in the dict per generator key for gen_det_I, det_I_conn_by_C in sd_by_C.items(): # This list contains all determinants connected to det_I - ref_det_I_conn = list(self.connected_by_det[gen_det_I]) + # Generate all singles and doubles of gen_det_I + det_I_sd = [det_J for det_J in gen_det_I.gen_all_connected_det(n_orb)] + ref_det_I_conn = [] # Empty list + for det_J in det_I_sd: + if det_J not in psi: + ref_det_I_conn.append(det_J) + # Filter out those that satisfy constraint C ref_det_I_conn_by_C = [ - det_J for det_J in ref_det_I_conn if (self.check_constraint(det_J) == C) + det_J for det_J in ref_det_I_conn if (check_constraint(det_J) == C) ] - self.assertListEqual(sorted(ref_det_I_conn_by_C), sorted(det_I_conn_by_C)) + ref_det_I_conn_by_C = sorted(ref_det_I_conn_by_C, key=lambda x: x.alpha) + ref_det_I_conn_by_C = sorted(ref_det_I_conn_by_C, key=lambda x: x.beta) + + det_I_conn_by_C = sorted(det_I_conn_by_C, key=lambda x: x.alpha) + det_I_conn_by_C = sorted(det_I_conn_by_C, key=lambda x: x.beta) + self.assertListEqual(ref_det_I_conn_by_C, det_I_conn_by_C) class Test_Integral_Driven_Categories(Test_Minimal):