-
Notifications
You must be signed in to change notification settings - Fork 7
Expand file tree
/
Copy pathhamiltonian_generator.hpp
More file actions
193 lines (152 loc) · 7.84 KB
/
hamiltonian_generator.hpp
File metadata and controls
193 lines (152 loc) · 7.84 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
/*
* MACIS Copyright (c) 2023, The Regents of the University of California,
* through Lawrence Berkeley National Laboratory (subject to receipt of
* any required approvals from the U.S. Dept. of Energy). All rights reserved.
*
* See LICENSE.txt for details
*/
#pragma once
#include <macis/bitset_operations.hpp>
#include <macis/sd_operations.hpp>
#include <macis/types.hpp>
#include <sparsexx/matrix_types/csr_matrix.hpp>
namespace macis {
template <size_t N = 64>
class HamiltonianGenerator {
static_assert(N % 2 == 0, "N Must Be Even");
public:
constexpr static size_t nbits = N;
using full_det_t = std::bitset<N>;
using spin_det_t = std::bitset<N / 2>;
template <typename index_t>
using sparse_matrix_type = sparsexx::csr_matrix<double, index_t>;
using full_det_iterator = wavefunction_iterator_t<N>;
using matrix_span_t = matrix_span<double>;
using rank3_span_t = rank3_span<double>;
using rank4_span_t = rank4_span<double>;
public:
inline spin_det_t alpha_string(full_det_t str) { return bitset_lo_word(str); }
inline spin_det_t beta_string(full_det_t str) { return bitset_hi_word(str); }
size_t norb_;
size_t norb2_;
size_t norb3_;
matrix_span_t T_pq_;
rank4_span_t V_pqrs_;
// G(i,j,k,l) = (ij|kl) - (il|kj)
std::vector<double> G_pqrs_data_;
rank4_span_t G_pqrs_;
// G_red(i,j,k) = G(i,j,k,k)
std::vector<double> G_red_data_;
rank3_span_t G_red_;
// V_red(i,j,k) = (ij|kk)
std::vector<double> V_red_data_;
rank3_span_t V_red_;
// G2_red(i,j) = 0.5 * G(i,i,j,j)
std::vector<double> G2_red_data_;
matrix_span_t G2_red_;
// V2_red(i,j) = (ii|jj)
std::vector<double> V2_red_data_;
matrix_span_t V2_red_;
virtual sparse_matrix_type<int32_t> make_csr_hamiltonian_block_32bit_(
full_det_iterator, full_det_iterator, full_det_iterator,
full_det_iterator, double) = 0;
virtual sparse_matrix_type<int64_t> make_csr_hamiltonian_block_64bit_(
full_det_iterator, full_det_iterator, full_det_iterator,
full_det_iterator, double) = 0;
public:
HamiltonianGenerator(matrix_span_t T, rank4_span_t V);
virtual ~HamiltonianGenerator() noexcept = default;
void generate_integral_intermediates(rank4_span_t V);
inline auto* T() const { return T_pq_.data_handle(); }
inline auto* G_red() const { return G_red_data_.data(); }
inline auto* V_red() const { return V_red_data_.data(); }
inline auto* G() const { return G_pqrs_data_.data(); }
inline auto* V() const { return V_pqrs_.data_handle(); }
double matrix_element_4(spin_det_t bra, spin_det_t ket, spin_det_t ex) const;
double matrix_element_22(spin_det_t bra_alpha, spin_det_t ket_alpha,
spin_det_t ex_alpha, spin_det_t bra_beta,
spin_det_t ket_beta, spin_det_t ex_beta) const;
double matrix_element_2(spin_det_t bra, spin_det_t ket, spin_det_t ex,
const std::vector<uint32_t>& bra_occ_alpha,
const std::vector<uint32_t>& bra_occ_beta) const;
double matrix_element_diag(const std::vector<uint32_t>& occ_alpha,
const std::vector<uint32_t>& occ_beta) const;
double matrix_element(spin_det_t bra_alpha, spin_det_t ket_alpha,
spin_det_t ex_alpha, spin_det_t bra_beta,
spin_det_t ket_beta, spin_det_t ex_beta,
const std::vector<uint32_t>& bra_occ_alpha,
const std::vector<uint32_t>& bra_occ_beta) const;
double single_orbital_en(uint32_t orb, const std::vector<uint32_t>& ss_occ,
const std::vector<uint32_t>& os_occ) const;
std::vector<double> single_orbital_ens(
size_t norb, const std::vector<uint32_t>& ss_occ,
const std::vector<uint32_t>& os_occ) const;
double fast_diag_single(const std::vector<uint32_t>& ss_occ,
const std::vector<uint32_t>& os_occ, uint32_t orb_hol,
uint32_t orb_par, double orig_det_Hii) const;
double fast_diag_single(double hol_en, double par_en, uint32_t orb_hol,
uint32_t orb_par, double orig_det_Hii) const;
double fast_diag_ss_double(double en_hol1, double en_hol2, double en_par1,
double en_par2, uint32_t orb_hol1,
uint32_t orb_hol2, uint32_t orb_par1,
uint32_t orb_par2, double orig_det_Hii) const;
double fast_diag_ss_double(const std::vector<uint32_t>& ss_occ,
const std::vector<uint32_t>& os_occ,
uint32_t orb_hol1, uint32_t orb_hol2,
uint32_t orb_par1, uint32_t orb_par2,
double orig_det_Hii) const;
double fast_diag_os_double(double en_holu, double en_hold, double en_paru,
double en_pard, uint32_t orb_holu,
uint32_t orb_hold, uint32_t orb_paru,
uint32_t orb_pard, double orig_det_Hii) const;
double fast_diag_os_double(const std::vector<uint32_t>& up_occ,
const std::vector<uint32_t>& do_occ,
uint32_t orb_holu, uint32_t orb_hold,
uint32_t orb_paru, uint32_t orb_pard,
double orig_det_Hii) const;
double matrix_element(full_det_t bra, full_det_t ket) const;
template <typename index_t>
sparse_matrix_type<index_t> make_csr_hamiltonian_block(
full_det_iterator bra_begin, full_det_iterator bra_end,
full_det_iterator ket_begin, full_det_iterator ket_end, double H_thresh) {
if constexpr(std::is_same_v<index_t, int32_t>)
return make_csr_hamiltonian_block_32bit_(bra_begin, bra_end, ket_begin,
ket_end, H_thresh);
else if constexpr(std::is_same_v<index_t, int64_t>)
return make_csr_hamiltonian_block_64bit_(bra_begin, bra_end, ket_begin,
ket_end, H_thresh);
else {
throw std::runtime_error("Unsupported index_t");
abort();
}
}
void rdm_contributions_4(spin_det_t bra, spin_det_t ket, spin_det_t ex,
double val, rank4_span_t trdm);
void rdm_contributions_22(spin_det_t bra_alpha, spin_det_t ket_alpha,
spin_det_t ex_alpha, spin_det_t bra_beta,
spin_det_t ket_beta, spin_det_t ex_beta, double val,
rank4_span_t trdm);
void rdm_contributions_2(spin_det_t bra, spin_det_t ket, spin_det_t ex,
const std::vector<uint32_t>& bra_occ_alpha,
const std::vector<uint32_t>& bra_occ_beta,
double val, matrix_span_t ordm, rank4_span_t trdm);
void rdm_contributions_diag(const std::vector<uint32_t>& occ_alpha,
const std::vector<uint32_t>& occ_beta, double val,
matrix_span_t ordm, rank4_span_t trdm);
void rdm_contributions(spin_det_t bra_alpha, spin_det_t ket_alpha,
spin_det_t ex_alpha, spin_det_t bra_beta,
spin_det_t ket_beta, spin_det_t ex_beta,
const std::vector<uint32_t>& bra_occ_alpha,
const std::vector<uint32_t>& bra_occ_beta, double val,
matrix_span_t ordm, rank4_span_t trdm);
virtual void form_rdms(full_det_iterator, full_det_iterator,
full_det_iterator, full_det_iterator, double* C,
matrix_span_t ordm, rank4_span_t trdm) = 0;
void rotate_hamiltonian_ordm(const double* ordm);
virtual void SetJustSingles(bool /*_js*/) {}
virtual bool GetJustSingles() const { return false; }
virtual size_t GetNimp() const { return N / 2; }
};
} // namespace macis
// Implementation
#include <macis/hamiltonian_generator/impl.hpp>