@@ -76,140 +76,142 @@ class SortedDoubleLoopHamiltonianGenerator
7676 // << std::endl;
7777
7878 // Populate COO matrix locally
79- // sparsexx::coo_matrix<double, index_t> coo_mat(nbra_dets, nket_dets, 0, 0);
79+ // sparsexx::coo_matrix<double, index_t> coo_mat(nbra_dets, nket_dets, 0,
80+ // 0);
8081 std::vector<index_t > row_ind, col_ind;
81- std::vector<double > nz_val;
82+ std::vector<double > nz_val;
8283
83- // size_t skip1 = 0;
84- // size_t skip2 = 0;
84+ // size_t skip1 = 0;
85+ // size_t skip2 = 0;
8586
8687 std::mutex coo_mat_thread_mutex;
8788
8889 // Loop over uniq alphas in bra/ket
8990 auto pop_st = std::chrono::high_resolution_clock::now ();
90- #pragma omp parallel
91+ #pragma omp parallel
9192 {
92- std::vector<index_t > row_ind_loc, col_ind_loc;
93- std::vector<double > nz_val_loc;
94- std::vector<uint32_t > bra_occ_alpha, bra_occ_beta;
95- #pragma omp for schedule(dynamic)
96- for (size_t ia_bra = 0 ; ia_bra < nuniq_bra; ++ia_bra) {
97- if (unique_alpha_bra[ia_bra].first .any ()) {
98- // Extract alpha bra
99- const auto bra_alpha = unique_alpha_bra[ia_bra].first ;
100- const size_t beta_st_bra = unique_alpha_bra_idx[ia_bra];
101- const size_t beta_en_bra = unique_alpha_bra_idx[ia_bra + 1 ];
102- spin_wfn_traits::state_to_occ (bra_alpha, bra_occ_alpha);
103-
104- const auto ket_lower = is_symm ? ia_bra : 0 ;
105- for (size_t ia_ket = ket_lower; ia_ket < nuniq_ket; ++ia_ket) {
106- if (unique_alpha_ket[ia_ket].first .any ()) {
107- // Extract alpha ket
108- const auto ket_alpha = unique_alpha_ket[ia_ket].first ;
109-
110- // Compute alpha excitation
111- const auto ex_alpha = bra_alpha ^ ket_alpha;
112- const auto ex_alpha_count = spin_wfn_traits::count (ex_alpha);
113-
114- // Early exit
115- if (ex_alpha_count > 4 ) {
116- // skip1++;
117- continue ;
118- }
93+ std::vector<index_t > row_ind_loc, col_ind_loc;
94+ std::vector<double > nz_val_loc;
95+ std::vector<uint32_t > bra_occ_alpha, bra_occ_beta;
96+ #pragma omp for schedule(dynamic)
97+ for (size_t ia_bra = 0 ; ia_bra < nuniq_bra; ++ia_bra) {
98+ if (unique_alpha_bra[ia_bra].first .any ()) {
99+ // Extract alpha bra
100+ const auto bra_alpha = unique_alpha_bra[ia_bra].first ;
101+ const size_t beta_st_bra = unique_alpha_bra_idx[ia_bra];
102+ const size_t beta_en_bra = unique_alpha_bra_idx[ia_bra + 1 ];
103+ spin_wfn_traits::state_to_occ (bra_alpha, bra_occ_alpha);
104+
105+ const auto ket_lower = is_symm ? ia_bra : 0 ;
106+ for (size_t ia_ket = ket_lower; ia_ket < nuniq_ket; ++ia_ket) {
107+ if (unique_alpha_ket[ia_ket].first .any ()) {
108+ // Extract alpha ket
109+ const auto ket_alpha = unique_alpha_ket[ia_ket].first ;
110+
111+ // Compute alpha excitation
112+ const auto ex_alpha = bra_alpha ^ ket_alpha;
113+ const auto ex_alpha_count = spin_wfn_traits::count (ex_alpha);
114+
115+ // Early exit
116+ if (ex_alpha_count > 4 ) {
117+ // skip1++;
118+ continue ;
119+ }
119120
120- // Precompute all-alpha excitation if it will be used
121- const double mat_el_4_alpha =
122- (ex_alpha_count == 4 )
123- ? this ->matrix_element_4 (bra_alpha, ket_alpha, ex_alpha)
124- : 0.0 ;
125- if (ex_alpha_count == 4 and std::abs (mat_el_4_alpha) < H_thresh) {
126- // The only possible matrix element is too-small, skip everyhing
127- // skip2++;
128- continue ;
129- }
121+ // Precompute all-alpha excitation if it will be used
122+ const double mat_el_4_alpha =
123+ (ex_alpha_count == 4 )
124+ ? this ->matrix_element_4 (bra_alpha, ket_alpha, ex_alpha)
125+ : 0.0 ;
126+ if (ex_alpha_count == 4 and std::abs (mat_el_4_alpha) < H_thresh) {
127+ // The only possible matrix element is too-small, skip everyhing
128+ // skip2++;
129+ continue ;
130+ }
130131
131- const size_t beta_st_ket = unique_alpha_ket_idx[ia_ket];
132- const size_t beta_en_ket = unique_alpha_ket_idx[ia_ket + 1 ];
133-
134- // Loop over local betas according to their global indices
135- for (size_t ibra = beta_st_bra; ibra < beta_en_bra; ++ibra) {
136- const auto bra_beta =
137- wfn_traits::beta_string (*(bra_begin + ibra));
138- spin_wfn_traits::state_to_occ (bra_beta, bra_occ_beta);
139- for (size_t iket = beta_st_ket; iket < beta_en_ket; ++iket) {
140- if (is_symm and (iket < ibra)) continue ;
141- const auto ket_beta =
142- wfn_traits::beta_string (*(ket_begin + iket));
143-
144- // Compute beta excitation
145- const auto ex_beta = bra_beta ^ ket_beta;
146- const auto ex_beta_count = spin_wfn_traits::count (ex_beta);
147-
148- if ((ex_alpha_count + ex_beta_count) > 4 ) continue ;
149-
150- double h_el = 0.0 ;
151- if (ex_alpha_count == 4 ) {
152- // Use precomputed value
153- h_el = mat_el_4_alpha;
154- } else if (ex_beta_count == 4 ) {
155- h_el = this ->matrix_element_4 (bra_beta, ket_beta, ex_beta);
156- } else if (ex_alpha_count == 2 ) {
157- if (ex_beta_count == 2 ) {
158- h_el =
159- this ->matrix_element_22 (bra_alpha, ket_alpha, ex_alpha,
160- bra_beta, ket_beta, ex_beta);
132+ const size_t beta_st_ket = unique_alpha_ket_idx[ia_ket];
133+ const size_t beta_en_ket = unique_alpha_ket_idx[ia_ket + 1 ];
134+
135+ // Loop over local betas according to their global indices
136+ for (size_t ibra = beta_st_bra; ibra < beta_en_bra; ++ibra) {
137+ const auto bra_beta =
138+ wfn_traits::beta_string (*(bra_begin + ibra));
139+ spin_wfn_traits::state_to_occ (bra_beta, bra_occ_beta);
140+ for (size_t iket = beta_st_ket; iket < beta_en_ket; ++iket) {
141+ if (is_symm and (iket < ibra)) continue ;
142+ const auto ket_beta =
143+ wfn_traits::beta_string (*(ket_begin + iket));
144+
145+ // Compute beta excitation
146+ const auto ex_beta = bra_beta ^ ket_beta;
147+ const auto ex_beta_count = spin_wfn_traits::count (ex_beta);
148+
149+ if ((ex_alpha_count + ex_beta_count) > 4 ) continue ;
150+
151+ double h_el = 0.0 ;
152+ if (ex_alpha_count == 4 ) {
153+ // Use precomputed value
154+ h_el = mat_el_4_alpha;
155+ } else if (ex_beta_count == 4 ) {
156+ h_el = this ->matrix_element_4 (bra_beta, ket_beta, ex_beta);
157+ } else if (ex_alpha_count == 2 ) {
158+ if (ex_beta_count == 2 ) {
159+ h_el = this ->matrix_element_22 (bra_alpha, ket_alpha,
160+ ex_alpha, bra_beta,
161+ ket_beta, ex_beta);
162+ } else {
163+ h_el =
164+ this ->matrix_element_2 (bra_alpha, ket_alpha, ex_alpha,
165+ bra_occ_alpha, bra_occ_beta);
166+ }
167+ } else if (ex_beta_count == 2 ) {
168+ h_el = this ->matrix_element_2 (bra_beta, ket_beta, ex_beta,
169+ bra_occ_beta, bra_occ_alpha);
161170 } else {
171+ // Diagonal matrix element
162172 h_el =
163- this ->matrix_element_2 (bra_alpha, ket_alpha, ex_alpha,
164- bra_occ_alpha, bra_occ_beta);
165- }
166- } else if (ex_beta_count == 2 ) {
167- h_el = this ->matrix_element_2 (bra_beta, ket_beta, ex_beta,
168- bra_occ_beta, bra_occ_alpha);
169- } else {
170- // Diagonal matrix element
171- h_el = this ->matrix_element_diag (bra_occ_alpha, bra_occ_beta);
172- }
173-
174- // Insert matrix element
175- if (std::abs (h_el) > H_thresh) {
176- // coo_mat.template insert<false>(ibra, iket, h_el);
177- row_ind_loc.emplace_back (ibra);
178- col_ind_loc.emplace_back (iket);
179- nz_val_loc. emplace_back (h_el);
180- if (is_symm and ibra != iket) {
181- // coo_mat.template insert<false>(iket, ibra, h_el);
182- row_ind_loc.emplace_back (iket);
183- col_ind_loc.emplace_back (ibra);
184- nz_val_loc. emplace_back (h_el);
173+ this ->matrix_element_diag (bra_occ_alpha, bra_occ_beta);
185174 }
186- }
187175
188- } // ket beta
189- } // bra beta
190-
191- }
192- } // Loop over ket alphas
193- }
194- } // Loop over bra alphas
195-
196- // Atomically insert into larger matrix arrays
197- #pragma omp critical
198- {
199- row_ind.insert (row_ind.end (), row_ind_loc.begin (), row_ind_loc.end ());
200- // row_ind_loc.clear(); row_ind_loc.shrink_to_fit();
201- col_ind.insert (col_ind.end (), col_ind_loc.begin (), col_ind_loc.end ());
202- // col_ind_loc.clear(); col_ind_loc.shrink_to_fit();
203- nz_val.insert (nz_val.end (), nz_val_loc.begin (), nz_val_loc.end ());
204- // nz_val_loc.clear(); nz_val_loc.shrink_to_fit();
205- }
176+ // Insert matrix element
177+ if (std::abs (h_el) > H_thresh) {
178+ // coo_mat.template insert<false>(ibra, iket, h_el);
179+ row_ind_loc.emplace_back (ibra);
180+ col_ind_loc.emplace_back (iket);
181+ nz_val_loc.emplace_back (h_el);
182+ if (is_symm and ibra != iket) {
183+ // coo_mat.template insert<false>(iket, ibra, h_el);
184+ row_ind_loc.emplace_back (iket);
185+ col_ind_loc.emplace_back (ibra);
186+ nz_val_loc.emplace_back (h_el);
187+ }
188+ }
206189
207- } // OpenMP
190+ } // ket beta
191+ } // bra beta
192+ }
193+ } // Loop over ket alphas
194+ }
195+ } // Loop over bra alphas
196+
197+ // Atomically insert into larger matrix arrays
198+ #pragma omp critical
199+ {
200+ row_ind.insert (row_ind.end (), row_ind_loc.begin (), row_ind_loc.end ());
201+ // row_ind_loc.clear(); row_ind_loc.shrink_to_fit();
202+ col_ind.insert (col_ind.end (), col_ind_loc.begin (), col_ind_loc.end ());
203+ // col_ind_loc.clear(); col_ind_loc.shrink_to_fit();
204+ nz_val.insert (nz_val.end (), nz_val_loc.begin (), nz_val_loc.end ());
205+ // nz_val_loc.clear(); nz_val_loc.shrink_to_fit();
206+ }
207+
208+ } // OpenMP
208209 auto pop_en = std::chrono::high_resolution_clock::now ();
209210
210211 // Generate Sparse Matrix
211- sparsexx::coo_matrix<double , index_t > coo_mat (nbra_dets, nket_dets,
212- std::move (col_ind), std::move (row_ind), std::move (nz_val), 0 );
212+ sparsexx::coo_matrix<double , index_t > coo_mat (
213+ nbra_dets, nket_dets, std::move (col_ind), std::move (row_ind),
214+ std::move (nz_val), 0 );
213215
214216 // Sort for CSR Conversion
215217 auto sort_st = std::chrono::high_resolution_clock::now ();
0 commit comments