diff --git a/include/macis/asci/pt2.hpp b/include/macis/asci/pt2.hpp index 7b9af4d2..9df5fd81 100644 --- a/include/macis/asci/pt2.hpp +++ b/include/macis/asci/pt2.hpp @@ -129,7 +129,6 @@ double asci_pt2_constraint(wavefunction_iterator_t cdets_begin, asci_pairs.reserve(100000000ul); size_t ic = 0; while(ic < ncon_total) { - // Atomically get the next task ID and increment for other // MPI ranks and threads size_t ntake = ic < 1000 ? 1 : 10; diff --git a/include/macis/hamiltonian_generator/sorted_double_loop.hpp b/include/macis/hamiltonian_generator/sorted_double_loop.hpp index bf6c496f..92333898 100644 --- a/include/macis/hamiltonian_generator/sorted_double_loop.hpp +++ b/include/macis/hamiltonian_generator/sorted_double_loop.hpp @@ -76,140 +76,142 @@ class SortedDoubleLoopHamiltonianGenerator // << std::endl; // Populate COO matrix locally - //sparsexx::coo_matrix coo_mat(nbra_dets, nket_dets, 0, 0); + // sparsexx::coo_matrix coo_mat(nbra_dets, nket_dets, 0, + // 0); std::vector row_ind, col_ind; - std::vector nz_val; + std::vector nz_val; - //size_t skip1 = 0; - //size_t skip2 = 0; + // size_t skip1 = 0; + // size_t skip2 = 0; std::mutex coo_mat_thread_mutex; // Loop over uniq alphas in bra/ket auto pop_st = std::chrono::high_resolution_clock::now(); - #pragma omp parallel +#pragma omp parallel { - std::vector row_ind_loc, col_ind_loc; - std::vector nz_val_loc; - std::vector bra_occ_alpha, bra_occ_beta; - #pragma omp for schedule(dynamic) - for(size_t ia_bra = 0; ia_bra < nuniq_bra; ++ia_bra) { - if(unique_alpha_bra[ia_bra].first.any()) { - // Extract alpha bra - const auto bra_alpha = unique_alpha_bra[ia_bra].first; - const size_t beta_st_bra = unique_alpha_bra_idx[ia_bra]; - const size_t beta_en_bra = unique_alpha_bra_idx[ia_bra + 1]; - spin_wfn_traits::state_to_occ(bra_alpha, bra_occ_alpha); - - const auto ket_lower = is_symm ? ia_bra : 0; - for(size_t ia_ket = ket_lower; ia_ket < nuniq_ket; ++ia_ket) { - if(unique_alpha_ket[ia_ket].first.any()) { - // Extract alpha ket - const auto ket_alpha = unique_alpha_ket[ia_ket].first; - - // Compute alpha excitation - const auto ex_alpha = bra_alpha ^ ket_alpha; - const auto ex_alpha_count = spin_wfn_traits::count(ex_alpha); - - // Early exit - if(ex_alpha_count > 4) { - //skip1++; - continue; - } + std::vector row_ind_loc, col_ind_loc; + std::vector nz_val_loc; + std::vector bra_occ_alpha, bra_occ_beta; +#pragma omp for schedule(dynamic) + for(size_t ia_bra = 0; ia_bra < nuniq_bra; ++ia_bra) { + if(unique_alpha_bra[ia_bra].first.any()) { + // Extract alpha bra + const auto bra_alpha = unique_alpha_bra[ia_bra].first; + const size_t beta_st_bra = unique_alpha_bra_idx[ia_bra]; + const size_t beta_en_bra = unique_alpha_bra_idx[ia_bra + 1]; + spin_wfn_traits::state_to_occ(bra_alpha, bra_occ_alpha); + + const auto ket_lower = is_symm ? ia_bra : 0; + for(size_t ia_ket = ket_lower; ia_ket < nuniq_ket; ++ia_ket) { + if(unique_alpha_ket[ia_ket].first.any()) { + // Extract alpha ket + const auto ket_alpha = unique_alpha_ket[ia_ket].first; + + // Compute alpha excitation + const auto ex_alpha = bra_alpha ^ ket_alpha; + const auto ex_alpha_count = spin_wfn_traits::count(ex_alpha); + + // Early exit + if(ex_alpha_count > 4) { + // skip1++; + continue; + } - // Precompute all-alpha excitation if it will be used - const double mat_el_4_alpha = - (ex_alpha_count == 4) - ? this->matrix_element_4(bra_alpha, ket_alpha, ex_alpha) - : 0.0; - if(ex_alpha_count == 4 and std::abs(mat_el_4_alpha) < H_thresh) { - // The only possible matrix element is too-small, skip everyhing - //skip2++; - continue; - } + // Precompute all-alpha excitation if it will be used + const double mat_el_4_alpha = + (ex_alpha_count == 4) + ? this->matrix_element_4(bra_alpha, ket_alpha, ex_alpha) + : 0.0; + if(ex_alpha_count == 4 and std::abs(mat_el_4_alpha) < H_thresh) { + // The only possible matrix element is too-small, skip everyhing + // skip2++; + continue; + } - const size_t beta_st_ket = unique_alpha_ket_idx[ia_ket]; - const size_t beta_en_ket = unique_alpha_ket_idx[ia_ket + 1]; - - // Loop over local betas according to their global indices - for(size_t ibra = beta_st_bra; ibra < beta_en_bra; ++ibra) { - const auto bra_beta = - wfn_traits::beta_string(*(bra_begin + ibra)); - spin_wfn_traits::state_to_occ(bra_beta, bra_occ_beta); - for(size_t iket = beta_st_ket; iket < beta_en_ket; ++iket) { - if(is_symm and (iket < ibra)) continue; - const auto ket_beta = - wfn_traits::beta_string(*(ket_begin + iket)); - - // Compute beta excitation - const auto ex_beta = bra_beta ^ ket_beta; - const auto ex_beta_count = spin_wfn_traits::count(ex_beta); - - if((ex_alpha_count + ex_beta_count) > 4) continue; - - double h_el = 0.0; - if(ex_alpha_count == 4) { - // Use precomputed value - h_el = mat_el_4_alpha; - } else if(ex_beta_count == 4) { - h_el = this->matrix_element_4(bra_beta, ket_beta, ex_beta); - } else if(ex_alpha_count == 2) { - if(ex_beta_count == 2) { - h_el = - this->matrix_element_22(bra_alpha, ket_alpha, ex_alpha, - bra_beta, ket_beta, ex_beta); + const size_t beta_st_ket = unique_alpha_ket_idx[ia_ket]; + const size_t beta_en_ket = unique_alpha_ket_idx[ia_ket + 1]; + + // Loop over local betas according to their global indices + for(size_t ibra = beta_st_bra; ibra < beta_en_bra; ++ibra) { + const auto bra_beta = + wfn_traits::beta_string(*(bra_begin + ibra)); + spin_wfn_traits::state_to_occ(bra_beta, bra_occ_beta); + for(size_t iket = beta_st_ket; iket < beta_en_ket; ++iket) { + if(is_symm and (iket < ibra)) continue; + const auto ket_beta = + wfn_traits::beta_string(*(ket_begin + iket)); + + // Compute beta excitation + const auto ex_beta = bra_beta ^ ket_beta; + const auto ex_beta_count = spin_wfn_traits::count(ex_beta); + + if((ex_alpha_count + ex_beta_count) > 4) continue; + + double h_el = 0.0; + if(ex_alpha_count == 4) { + // Use precomputed value + h_el = mat_el_4_alpha; + } else if(ex_beta_count == 4) { + h_el = this->matrix_element_4(bra_beta, ket_beta, ex_beta); + } else if(ex_alpha_count == 2) { + if(ex_beta_count == 2) { + h_el = this->matrix_element_22(bra_alpha, ket_alpha, + ex_alpha, bra_beta, + ket_beta, ex_beta); + } else { + h_el = + this->matrix_element_2(bra_alpha, ket_alpha, ex_alpha, + bra_occ_alpha, bra_occ_beta); + } + } else if(ex_beta_count == 2) { + h_el = this->matrix_element_2(bra_beta, ket_beta, ex_beta, + bra_occ_beta, bra_occ_alpha); } else { + // Diagonal matrix element h_el = - this->matrix_element_2(bra_alpha, ket_alpha, ex_alpha, - bra_occ_alpha, bra_occ_beta); - } - } else if(ex_beta_count == 2) { - h_el = this->matrix_element_2(bra_beta, ket_beta, ex_beta, - bra_occ_beta, bra_occ_alpha); - } else { - // Diagonal matrix element - h_el = this->matrix_element_diag(bra_occ_alpha, bra_occ_beta); - } - - // Insert matrix element - if(std::abs(h_el) > H_thresh) { - //coo_mat.template insert(ibra, iket, h_el); - row_ind_loc.emplace_back(ibra); - col_ind_loc.emplace_back(iket); - nz_val_loc. emplace_back(h_el); - if(is_symm and ibra != iket) { - //coo_mat.template insert(iket, ibra, h_el); - row_ind_loc.emplace_back(iket); - col_ind_loc.emplace_back(ibra); - nz_val_loc. emplace_back(h_el); + this->matrix_element_diag(bra_occ_alpha, bra_occ_beta); } - } - } // ket beta - } // bra beta - - } - } // Loop over ket alphas - } - } // Loop over bra alphas - - // Atomically insert into larger matrix arrays - #pragma omp critical - { - row_ind.insert(row_ind.end(), row_ind_loc.begin(), row_ind_loc.end()); - //row_ind_loc.clear(); row_ind_loc.shrink_to_fit(); - col_ind.insert(col_ind.end(), col_ind_loc.begin(), col_ind_loc.end()); - //col_ind_loc.clear(); col_ind_loc.shrink_to_fit(); - nz_val.insert(nz_val.end(), nz_val_loc.begin(), nz_val_loc.end()); - //nz_val_loc.clear(); nz_val_loc.shrink_to_fit(); - } + // Insert matrix element + if(std::abs(h_el) > H_thresh) { + // coo_mat.template insert(ibra, iket, h_el); + row_ind_loc.emplace_back(ibra); + col_ind_loc.emplace_back(iket); + nz_val_loc.emplace_back(h_el); + if(is_symm and ibra != iket) { + // coo_mat.template insert(iket, ibra, h_el); + row_ind_loc.emplace_back(iket); + col_ind_loc.emplace_back(ibra); + nz_val_loc.emplace_back(h_el); + } + } - } // OpenMP + } // ket beta + } // bra beta + } + } // Loop over ket alphas + } + } // Loop over bra alphas + +// Atomically insert into larger matrix arrays +#pragma omp critical + { + row_ind.insert(row_ind.end(), row_ind_loc.begin(), row_ind_loc.end()); + // row_ind_loc.clear(); row_ind_loc.shrink_to_fit(); + col_ind.insert(col_ind.end(), col_ind_loc.begin(), col_ind_loc.end()); + // col_ind_loc.clear(); col_ind_loc.shrink_to_fit(); + nz_val.insert(nz_val.end(), nz_val_loc.begin(), nz_val_loc.end()); + // nz_val_loc.clear(); nz_val_loc.shrink_to_fit(); + } + + } // OpenMP auto pop_en = std::chrono::high_resolution_clock::now(); // Generate Sparse Matrix - sparsexx::coo_matrix coo_mat(nbra_dets, nket_dets, - std::move(col_ind), std::move(row_ind), std::move(nz_val), 0); + sparsexx::coo_matrix coo_mat( + nbra_dets, nket_dets, std::move(col_ind), std::move(row_ind), + std::move(nz_val), 0); // Sort for CSR Conversion auto sort_st = std::chrono::high_resolution_clock::now(); diff --git a/include/macis/util/mpi.hpp b/include/macis/util/mpi.hpp index 0d1a380f..47234b99 100644 --- a/include/macis/util/mpi.hpp +++ b/include/macis/util/mpi.hpp @@ -228,21 +228,17 @@ struct mpi_traits> { } }; - - - template class global_atomic { MPI_Win window_; - T* buffer_; - -public: + T* buffer_; + public: global_atomic() = delete; global_atomic(MPI_Comm comm) { MPI_Win_allocate(sizeof(T), sizeof(T), MPI_INFO_NULL, comm, &buffer_, - &window_); + &window_); if(window_ == MPI_WIN_NULL) { throw std::runtime_error("Window creation failed"); } @@ -257,16 +253,15 @@ class global_atomic { global_atomic(const global_atomic&) = delete; global_atomic(global_atomic&&) noexcept = delete; - + T fetch_and_add(T val) { T next_val; MPI_Fetch_and_op(&val, &next_val, mpi_traits::datatype(), 0, 0, MPI_SUM, - window_); - MPI_Win_flush(0,window_); + window_); + MPI_Win_flush(0, window_); return next_val; } }; - } // namespace macis #endif diff --git a/include/macis/util/trexio.hpp b/include/macis/util/trexio.hpp index 1400ef05..41ad9de5 100644 --- a/include/macis/util/trexio.hpp +++ b/include/macis/util/trexio.hpp @@ -7,31 +7,30 @@ * See LICENSE.txt for details */ #include + #include namespace macis { class trexio_exception : public std::exception { std::string msg_; - const char* what() const noexcept override; + const char* what() const noexcept override; -public: + public: trexio_exception(std::string func_name, trexio_exit_code rc); }; class TREXIOFile { - trexio_t * file_handle_ = nullptr; - -public: + trexio_t* file_handle_ = nullptr; + public: TREXIOFile() noexcept = default; - TREXIOFile(std::string name, char mode, int backend); + TREXIOFile(std::string name, char mode, int backend); ~TREXIOFile() noexcept; TREXIOFile(const TREXIOFile&) = delete; - TREXIOFile(TREXIOFile&& other) noexcept; - + TREXIOFile(TREXIOFile&& other) noexcept; int64_t read_mo_num() const; int64_t read_mo_2e_int_eri_size() const; @@ -43,7 +42,6 @@ class TREXIOFile { void write_nucleus_repulsion(double E); void write_mo_1e_int_core_hamiltonian(const double* h); void write_mo_2e_int_eri(const double* V); - }; -} +} // namespace macis diff --git a/src/macis/trexio.cxx b/src/macis/trexio.cxx index e4596051..3929c16d 100644 --- a/src/macis/trexio.cxx +++ b/src/macis/trexio.cxx @@ -6,26 +6,24 @@ * See LICENSE.txt for details */ +#include +#include #include #include -#include -#include namespace macis { -trexio_exception::trexio_exception(std::string func_name, trexio_exit_code rc) : - msg_("TREXIO Error: " + func_name + "\n Error Message: " + std::string(trexio_string_of_error(rc))) {} +trexio_exception::trexio_exception(std::string func_name, trexio_exit_code rc) + : msg_("TREXIO Error: " + func_name + + "\n Error Message: " + std::string(trexio_string_of_error(rc))) {} -const char* trexio_exception::what() const noexcept { - return msg_.c_str(); -} +const char* trexio_exception::what() const noexcept { return msg_.c_str(); } #define TREXIO_EXCEPTION(rc) throw trexio_exception(__PRETTY_FUNCTION__, rc) - TREXIOFile::TREXIOFile(std::string name, char mode, int backend) { trexio_exit_code rc; - file_handle_ = trexio_open( name.c_str(), mode, backend, &rc); + file_handle_ = trexio_open(name.c_str(), mode, backend, &rc); if(rc != TREXIO_SUCCESS) TREXIO_EXCEPTION(rc); } @@ -38,7 +36,6 @@ TREXIOFile::TREXIOFile(TREXIOFile&& other) noexcept { other.file_handle_ = nullptr; } - int64_t TREXIOFile::read_mo_num() const { int64_t nmo; auto rc = trexio_read_mo_num_64(file_handle_, &nmo); @@ -66,49 +63,44 @@ void TREXIOFile::read_mo_1e_int_core_hamiltonian(double* h) const { } void TREXIOFile::read_mo_2e_int_eri(double* V) const { - const auto norb = read_mo_num(); - const auto neri = read_mo_2e_int_eri_size(); + const auto norb = read_mo_num(); + const auto neri = read_mo_2e_int_eri_size(); const size_t nbatch = 100000; - std::vector> idx_batch(nbatch); + std::vector> idx_batch(nbatch); std::vector V_batch(nbatch); size_t ioff = 0; int64_t icount = nbatch; while(icount == nbatch) { if(ioff < neri) { - trexio_read_mo_2e_int_eri(file_handle_, ioff, &icount, - idx_batch[0].data(), V_batch.data()); + trexio_read_mo_2e_int_eri(file_handle_, ioff, &icount, + idx_batch[0].data(), V_batch.data()); } else { icount = 0; } for(size_t ii = 0; ii < icount; ++ii) { - const auto [i,j,k,l] = idx_batch[ii]; + const auto [i, j, k, l] = idx_batch[ii]; const auto v = V_batch[ii]; - V[i + j*norb + k*norb*norb + l*norb*norb*norb] = v; - V[i + j*norb + l*norb*norb + k*norb*norb*norb] = v; - V[j + i*norb + k*norb*norb + l*norb*norb*norb] = v; - V[j + i*norb + l*norb*norb + k*norb*norb*norb] = v; - V[k + l*norb + i*norb*norb + j*norb*norb*norb] = v; - V[k + l*norb + j*norb*norb + i*norb*norb*norb] = v; - V[l + k*norb + i*norb*norb + j*norb*norb*norb] = v; - V[l + k*norb + j*norb*norb + i*norb*norb*norb] = v; + V[i + j * norb + k * norb * norb + l * norb * norb * norb] = v; + V[i + j * norb + l * norb * norb + k * norb * norb * norb] = v; + V[j + i * norb + k * norb * norb + l * norb * norb * norb] = v; + V[j + i * norb + l * norb * norb + k * norb * norb * norb] = v; + V[k + l * norb + i * norb * norb + j * norb * norb * norb] = v; + V[k + l * norb + j * norb * norb + i * norb * norb * norb] = v; + V[l + k * norb + i * norb * norb + j * norb * norb * norb] = v; + V[l + k * norb + j * norb * norb + i * norb * norb * norb] = v; } ioff += icount; } } - - - - - void TREXIOFile::write_mo_num(int64_t nmo) { auto rc = trexio_write_mo_num_64(file_handle_, nmo); if(rc != TREXIO_SUCCESS) TREXIO_EXCEPTION(rc); -} - +} + void TREXIOFile::write_nucleus_repulsion(double E) { auto rc = trexio_write_nucleus_repulsion(file_handle_, E); if(rc != TREXIO_SUCCESS) TREXIO_EXCEPTION(rc); @@ -120,27 +112,28 @@ void TREXIOFile::write_mo_1e_int_core_hamiltonian(const double* h) { } void TREXIOFile::write_mo_2e_int_eri(const double* V) { - const auto norb = read_mo_num(); - const auto npair = (norb * (norb +1)) / 2; - const auto nquad = (npair * (npair+1)) / 2; + const auto norb = read_mo_num(); + const auto npair = (norb * (norb + 1)) / 2; + const auto nquad = (npair * (npair + 1)) / 2; std::vector V_compress(nquad); - std::vector> idx(nquad); + std::vector> idx(nquad); size_t ijkl = 0; - for(size_t i = 0, ij = 0; i < norb; ++i ) - for(size_t j = 0; j <= i; ++j, ++ij) { - for(size_t k = 0, kl = 0; k < norb; ++k ) - for(size_t l = 0; l <= k; ++l, ++kl) { - if(kl > ij) continue; - V_compress.at(ijkl) = V[i + j*norb + k*norb*norb + l*norb*norb*norb]; - idx.at(ijkl) = {i,j,k,l}; - ijkl++; + for(size_t i = 0, ij = 0; i < norb; ++i) + for(size_t j = 0; j <= i; ++j, ++ij) { + for(size_t k = 0, kl = 0; k < norb; ++k) + for(size_t l = 0; l <= k; ++l, ++kl) { + if(kl > ij) continue; + V_compress.at(ijkl) = + V[i + j * norb + k * norb * norb + l * norb * norb * norb]; + idx.at(ijkl) = {i, j, k, l}; + ijkl++; + } } - } - auto rc = trexio_write_mo_2e_int_eri(file_handle_, 0, nquad, - idx[0].data(), V_compress.data()); + auto rc = trexio_write_mo_2e_int_eri(file_handle_, 0, nquad, idx[0].data(), + V_compress.data()); if(rc != TREXIO_SUCCESS) TREXIO_EXCEPTION(rc); } -} +} // namespace macis diff --git a/tests/fcidump_to_trexio.cxx b/tests/fcidump_to_trexio.cxx index 25722238..aa18fa39 100644 --- a/tests/fcidump_to_trexio.cxx +++ b/tests/fcidump_to_trexio.cxx @@ -1,21 +1,20 @@ -#include -#include #include +#include +#include #include #include int main(int argc, char** argv) { - std::vector opts(argc); for(int i = 0; i < argc; ++i) { opts[i] = std::string(argv[i]); } std::string fcidump_fname = opts.at(1); - std::string trexio_fname = opts.at(2); + std::string trexio_fname = opts.at(2); std::cout << "FCIDUMP FILE = " << fcidump_fname << std::endl; - std::cout << "TREXIO FILE = " << trexio_fname << std::endl; + std::cout << "TREXIO FILE = " << trexio_fname << std::endl; // Read the FCIDUMP file std::cout << "Reading FCIDUMP File..." << std::flush; @@ -34,23 +33,24 @@ int main(int argc, char** argv) { // Write TREXIO file std::cout << "Writing TREXIO file..." << std::flush; { - macis::TREXIOFile trexio_file(trexio_fname, 'w', TREXIO_HDF5); - trexio_file.write_mo_num(norb); - trexio_file.write_nucleus_repulsion(E_core); - trexio_file.write_mo_1e_int_core_hamiltonian(T.data()); - trexio_file.write_mo_2e_int_eri(V.data()); + macis::TREXIOFile trexio_file(trexio_fname, 'w', TREXIO_HDF5); + trexio_file.write_mo_num(norb); + trexio_file.write_nucleus_repulsion(E_core); + trexio_file.write_mo_1e_int_core_hamiltonian(T.data()); + trexio_file.write_mo_2e_int_eri(V.data()); } std::cout << "Done!" << std::endl; - + std::vector T_chk(norb2), V_chk(norb4); - size_t norb_chk; double E_core_chk; + size_t norb_chk; + double E_core_chk; std::cout << "Reading TREXIO file..." << std::flush; { - macis::TREXIOFile trexio_file(trexio_fname, 'r', TREXIO_AUTO); - norb_chk = trexio_file.read_mo_num(); - E_core_chk = trexio_file.read_nucleus_repulsion(); - trexio_file.read_mo_1e_int_core_hamiltonian(T_chk.data()); - trexio_file.read_mo_2e_int_eri(V_chk.data()); + macis::TREXIOFile trexio_file(trexio_fname, 'r', TREXIO_AUTO); + norb_chk = trexio_file.read_mo_num(); + E_core_chk = trexio_file.read_nucleus_repulsion(); + trexio_file.read_mo_1e_int_core_hamiltonian(T_chk.data()); + trexio_file.read_mo_2e_int_eri(V_chk.data()); } std::cout << "Done!" << std::endl; @@ -58,10 +58,11 @@ int main(int argc, char** argv) { std::cout << "ECOR_CHECK = " << std::abs(E_core - E_core_chk) << std::endl; double max_diff = 0; - for(auto i = 0; i < norb2; ++i) max_diff = std::max(max_diff, std::abs(T[i] - T_chk[i])); + for(auto i = 0; i < norb2; ++i) + max_diff = std::max(max_diff, std::abs(T[i] - T_chk[i])); std::cout << "T_CHECK = " << max_diff << std::endl; max_diff = 0; - for(auto i = 0; i < norb2; ++i) max_diff = std::max(max_diff, std::abs(V[i] - V_chk[i])); + for(auto i = 0; i < norb2; ++i) + max_diff = std::max(max_diff, std::abs(V[i] - V_chk[i])); std::cout << "V_CHECK = " << max_diff << std::endl; - } diff --git a/tests/standalone_driver.cxx b/tests/standalone_driver.cxx index d813ccbc..c814bc13 100644 --- a/tests/standalone_driver.cxx +++ b/tests/standalone_driver.cxx @@ -23,11 +23,11 @@ #include #include #include -#include #include #include #include #include +#include #include #include #include @@ -101,9 +101,10 @@ int main(int argc, char** argv) { auto nalpha = input.getData("CI.NALPHA"); auto nbeta = input.getData("CI.NBETA"); - std::string reference_data_format = input.getData("CI.REF_DATA_FORMAT"); - std::string reference_data_file = input.getData("CI.REF_DATA_FILE"); - + std::string reference_data_format = + input.getData("CI.REF_DATA_FORMAT"); + std::string reference_data_file = + input.getData("CI.REF_DATA_FILE"); size_t norb, norb2, norb3, norb4; std::vector T, V; @@ -116,11 +117,12 @@ int main(int argc, char** argv) { norb4 = norb2 * norb2; // XXX: Consider reading this into shared memory to avoid replication - T.resize(norb2); V.resize(norb4); + T.resize(norb2); + V.resize(norb4); E_core = macis::read_fcidump_core(reference_data_file); macis::read_fcidump_1body(reference_data_file, T.data(), norb); macis::read_fcidump_2body(reference_data_file, V.data(), norb); - } else { // TREXIO + } else { // TREXIO macis::TREXIOFile trexio_file(reference_data_file, 'r', TREXIO_AUTO); norb = trexio_file.read_mo_num(); norb2 = norb * norb; @@ -128,7 +130,8 @@ int main(int argc, char** argv) { norb4 = norb2 * norb2; // XXX: Consider reading this into shared memory to avoid replication - T.resize(norb2); V.resize(norb4); + T.resize(norb2); + V.resize(norb4); E_core = trexio_file.read_nucleus_repulsion(); trexio_file.read_mo_1e_int_core_hamiltonian(T.data()); trexio_file.read_mo_2e_int_eri(V.data());