Skip to content

Restart fix symm #146

New issue

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

Sign up for GitHub

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

Already on GitHub? Sign in to your account

Merged
shivupa merged 5 commits into from
Apr 24, 2025
Merged

Restart fix symm #146

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
500a6d1
cleaning
shivupa Mar 3, 2025
3d57b36
bug fixes
shivupa Mar 3, 2025
f09bd55
Let restart populate symm orbitals
shivupa Mar 3, 2025
fd5b91c
save symm info to bin file
shivupa Mar 5, 2025
e08f437
read symm info
shivupa Mar 5, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 4 additions & 3 deletions src/calculation/calculation.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -504,8 +504,8 @@ void POLYQUANT_CALCULATION::dump_mf_for_qmcpack(std::string &filename) {
Polyquant_cout("Dumping HDF5 to filename: " + particle_filename);
POLYQUANT_HDF5 hdf5_f(particle_filename);
hdf5_f.dump_mf_to_hdf5_for_QMCPACK(pbc, ecp, complex_vals, restricted, num_ao, num_mo, bohr_unit, num_part_alpha, num_part_beta, num_part_total, multiplicity, num_atom, num_species,
quantum_part_name, scf_calc->E_orbitals_combined[quantum_part_idx], scf_calc->C_combined[quantum_part_idx], atomic_species_ids, atomic_number, atomic_charge,
core_elec, atomic_names, atomic_centers, unique_shells);
quantum_part_name, scf_calc->E_orbitals_combined[quantum_part_idx], scf_calc->C_combined[quantum_part_idx], scf_calc->symm_label_idxs[quantum_part_idx],
scf_calc->symm_labels[quantum_part_idx], atomic_species_ids, atomic_number, atomic_charge, core_elec, atomic_names, atomic_centers, unique_shells);
quantum_part_idx++;
}
}
Expand Down Expand Up @@ -612,7 +612,8 @@ void POLYQUANT_CALCULATION::dump_post_mf_NOs_for_qmcpack(std::string &filename)
Polyquant_cout("Dumping HDF5 to filename: " + particle_filename.str());
POLYQUANT_HDF5 hdf5_f(particle_filename.str());
hdf5_f.dump_mf_to_hdf5_for_QMCPACK(pbc, ecp, complex_vals, restricted, num_ao, num_mo, bohr_unit, num_part_alpha, num_part_beta, num_part_total, multiplicity, num_atom, num_species,
quantum_part_name, ci_calc->occ_nso[state_vec_idx][quantum_part_idx], ci_calc->C_nso[state_vec_idx][quantum_part_idx], atomic_species_ids, atomic_number,
quantum_part_name, ci_calc->occ_nso[state_vec_idx][quantum_part_idx], ci_calc->C_nso[state_vec_idx][quantum_part_idx],
ci_calc->symm_label_idxs[state_vec_idx][quantum_part_idx], ci_calc->symm_labels[state_vec_idx][quantum_part_idx], atomic_species_ids, atomic_number,
atomic_charge, core_elec, atomic_names, atomic_centers, unique_shells);
}
quantum_part_idx++;
Expand Down
15 changes: 11 additions & 4 deletions src/io/hdf5_utilities.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -157,7 +157,8 @@ void POLYQUANT_HDF5::dump_generalparameters(bool complex_vals, bool ecp, bool re
}

void POLYQUANT_HDF5::dump_MOs(std::string quantum_part_name, int num_ao, int num_mo, std::vector<Eigen::Matrix<double, Eigen::Dynamic, 1>> E_orb,
std::vector<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>> mo_coeff) {
std::vector<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>> mo_coeff, std::vector<std::vector<int>> mo_symm_label_idxs,
std::vector<std::vector<std::string>> mo_symm_labels) {
// write MO parameters
Polyquant_cout("dumping MOs");
std::string super_twist_group = "/Super_Twist";
Expand Down Expand Up @@ -190,6 +191,11 @@ void POLYQUANT_HDF5::dump_MOs(std::string quantum_part_name, int num_ao, int num
path = super_twist_group + "/eigenset_" + std::to_string(spin_idx);
// H5Easy::dump(*hdf5_file, path, flattened_mo_coeff, H5Easy::DumpMode::Overwrite);
H5Easy::dump(*hdf5_file, path, mo_coeff[spin_idx].transpose(), H5Easy::DumpMode::Overwrite);

path = super_twist_group + "/eigensymmlab_" + std::to_string(spin_idx);
H5Easy::dump(*hdf5_file, path, mo_symm_labels[spin_idx], H5Easy::DumpMode::Overwrite);
path = super_twist_group + "/eigensymmint_" + std::to_string(spin_idx);
H5Easy::dump(*hdf5_file, path, mo_symm_label_idxs[spin_idx], H5Easy::DumpMode::Overwrite);
}
}

Expand Down Expand Up @@ -376,14 +382,15 @@ void POLYQUANT_HDF5::dump_basis(std::string quantum_part_name, std::vector<std::
void POLYQUANT_HDF5::dump_mf_to_hdf5_for_QMCPACK(bool pbc, bool complex_vals, bool ecp, bool restricted, int num_ao, int num_mo, bool bohr_unit, int num_part_alpha, int num_part_beta,
int num_part_total, int multiplicity, int num_atom, int num_species, std::string quantum_part_name,
std::vector<Eigen::Matrix<double, Eigen::Dynamic, 1>> E_orb, std::vector<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>> mo_coeff,
std::vector<int> atomic_species_ids, std::vector<int> atomic_number, std::vector<int> atomic_charge, std::vector<int> core_elec,
std::vector<std::string> atomic_names, std::vector<std::vector<double>> atomic_centers, std::vector<std::vector<libint2::Shell>> unique_shells) {
std::vector<std::vector<int>> mo_symm_label_idxs, std::vector<std::vector<std::string>> mo_symm_labels, std::vector<int> atomic_species_ids,
std::vector<int> atomic_number, std::vector<int> atomic_charge, std::vector<int> core_elec, std::vector<std::string> atomic_names,
std::vector<std::vector<double>> atomic_centers, std::vector<std::vector<libint2::Shell>> unique_shells) {
// create file
Polyquant_cout("dumping file");
this->dump_application();
this->dump_PBC(pbc);
this->dump_generalparameters(complex_vals, ecp, restricted, num_ao, num_mo, bohr_unit, num_part_alpha, num_part_beta, num_part_total, multiplicity);
this->dump_MOs(quantum_part_name, num_ao, num_mo, E_orb, mo_coeff);
this->dump_MOs(quantum_part_name, num_ao, num_mo, E_orb, mo_coeff, mo_symm_label_idxs, mo_symm_labels);
this->dump_atoms(num_atom, num_species, atomic_species_ids, atomic_number, atomic_charge, core_elec, atomic_names, atomic_centers);
this->dump_basis(quantum_part_name, atomic_names, unique_shells);
}
Expand Down
7 changes: 4 additions & 3 deletions src/io/hdf5_utilities.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -70,12 +70,13 @@ class POLYQUANT_HDF5 {
std::vector<std::string> atomic_names, std::vector<std::vector<double>> atomic_centers);
void dump_generalparameters(bool complex_vals, bool ecp, bool restricted, int num_ao, int num_mo, bool bohr_unit, int num_part_alpha, int num_part_beta, int num_part_total, int multiplicity);
void dump_MOs(std::string quantum_part_name, int num_ao, int num_mo, std::vector<Eigen::Matrix<double, Eigen::Dynamic, 1>> E_orb,
std::vector<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>> mo_coeff);
std::vector<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>> mo_coeff, std::vector<std::vector<int>> mo_symm_label_idxs, std::vector<std::vector<std::string>> mo_symm_labels);
void dump_basis(std::string quantum_part_name, std::vector<std::string> atomic_names, std::vector<std::vector<libint2::Shell>> unique_shells);
void dump_mf_to_hdf5_for_QMCPACK(bool pbc, bool complex_vals, bool ecp, bool restricted, int num_ao, int num_mo, bool bohr_unit, int num_part_alpha, int num_part_beta, int num_part_total,
int multiplicity, int num_atom, int num_species, std::string quantum_part_name, std::vector<Eigen::Matrix<double, Eigen::Dynamic, 1>> E_orb,
std::vector<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>> mo_coeff, std::vector<int> atomic_species_ids, std::vector<int> atomic_number,
std::vector<int> atomic_charge, std::vector<int> core_elec, std::vector<std::string> atomic_names, std::vector<std::vector<double>> atomic_centers,
std::vector<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>> mo_coeff, std::vector<std::vector<int>> mo_symm_label_idxs,
std::vector<std::vector<std::string>> mo_symm_labels, std::vector<int> atomic_species_ids, std::vector<int> atomic_number, std::vector<int> atomic_charge,
std::vector<int> core_elec, std::vector<std::string> atomic_names, std::vector<std::vector<double>> atomic_centers,
std::vector<std::vector<libint2::Shell>> unique_shells);
void dump_post_mf_to_hdf5_for_QMCPACK(std::vector<std::vector<std::vector<std::vector<uint64_t>>>> dets, Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> C, int N_dets, int N_states, int N_mo);
};
Expand Down
93 changes: 65 additions & 28 deletions src/scf/epscf.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1075,7 +1075,6 @@ void POLYQUANT_EPSCF::print_iteration() {

void POLYQUANT_EPSCF::form_occ_helper_aufbau(Eigen::Matrix<double, Eigen::Dynamic, 1> &part_occ, const int quantum_part_idx, const int quantum_part_spin_idx, const int quantum_part_irrep_idx,
const int num_parts, const double occval) {

part_occ.setZero(part_occ.size());
for (auto i = 0; i < num_parts; i++) {
part_occ[i] = occval;
Expand Down Expand Up @@ -1626,30 +1625,15 @@ void POLYQUANT_EPSCF::setup_from_file(std::string &filename) {
POLYQUANT_HDF5 hdf5_file(file_to_load);
Polyquant_cout("Reading coefficients from file : " + file_to_load);

// if (!root_group.exists("Super_Twist")) {
// APP_ABORT("Reading coefficients failed. No Super_Twist group in HDF5 file.");
// }
// auto Super_Twist_group = root_group.get_group("Super_Twist");

auto num_basis = this->input_basis->num_basis[quantum_part_idx];
auto quantum_part_irrep_idx = 0;
auto num_mo = this->num_mo[quantum_part_idx];

Polyquant_cout("Reading eigenset_" + std::to_string(idx));
// std::vector<double> data;
// data.resize(num_mo * num_basis);
std::string hpath = "/Super_Twist/eigenset_" + std::to_string(idx);
// hdf5_file.load_data(data, hpath);
hdf5_file.load_data(this->C_combined[quantum_part_idx][0], hpath);
this->C_combined[quantum_part_idx][0].transposeInPlace();

// #pragma omp parallel for
// for (auto i = 0; i < num_mo; i++) {
// for (auto j = 0; j < num_basis; j++) {
// this->C_combined[quantum_part_idx][0](j, i) = data[i * num_basis + j];
// }
// }

reorthogonalize_MOs(this->C_combined[quantum_part_idx][0], quantum_part_idx);

if (this->input_symmetry->do_symmetry == false) {
Expand All @@ -1659,22 +1643,12 @@ void POLYQUANT_EPSCF::setup_from_file(std::string &filename) {
if (quantum_part.num_parts > 1 && quantum_part.restricted == false) {
hpath = "/Super_Twist/eigenset_" + std::to_string(idx + 1);
Polyquant_cout("Reading eigenset_" + std::to_string(idx + 1));
// data.clear();
// data.resize(num_mo * num_basis);
// hdf5_file.load_data(data, hpath);
// #pragma omp parallel for
// for (auto i = 0; i < num_mo; i++) {
// for (auto j = 0; j < num_basis; j++) {
// this->C_combined[quantum_part_idx][1](j, i) = data[i * num_basis + j];
// }
// }
hdf5_file.load_data(this->C_combined[quantum_part_idx][1], hpath);
this->C_combined[quantum_part_idx][1].transposeInPlace();

reorthogonalize_MOs(this->C_combined[quantum_part_idx][1], quantum_part_idx);
if (this->input_symmetry->do_symmetry == false) {
auto quantum_part_irrep_idx = 0;
this->C[quantum_part_idx][1][quantum_part_irrep_idx] = this->C_combined[quantum_part_idx][0];
this->C[quantum_part_idx][1][quantum_part_irrep_idx] = this->C_combined[quantum_part_idx][1];
}
}
idx += 2;
Expand All @@ -1684,8 +1658,70 @@ void POLYQUANT_EPSCF::setup_from_file(std::string &filename) {
permute_initial_MOs();
}
// this->form_combined_orbitals();
// todo. What if we have symm info saved to bin?
bool symm_info_not_populated = true;

if (this->input_symmetry->do_symmetry == true and symm_info_not_populated) {
try {
Polyquant_cout("Attempting to read symmetry info from h5 file...");
quantum_part_idx = 0ul;
for (auto const &[quantum_part_key, quantum_part] : this->input_molecule->quantum_particles) {
std::filesystem::path path(filename);
std::string dir = path.parent_path().string();
std::string file = path.filename().string();
std::string file_to_load = dir + quantum_part_key + "_" + file;

POLYQUANT_HDF5 hdf5_file(file_to_load);

auto spin_idx = 0;
// read if we have this in bin
std::string hpath = "/Super_Twist/eigensymmlab_" + std::to_string(spin_idx);
hdf5_file.load_data(this->symm_labels[quantum_part_idx][spin_idx], hpath);
hpath = "/Super_Twist/eigensymmint_" + std::to_string(spin_idx);
hdf5_file.load_data(this->symm_label_idxs[quantum_part_idx][spin_idx], hpath);
if (quantum_part.num_parts > 1 && quantum_part.restricted == false) {
spin_idx = 1;
hpath = "/Super_Twist/eigensymmlab_" + std::to_string(spin_idx);
hdf5_file.load_data(this->symm_labels[quantum_part_idx][spin_idx], hpath);
hpath = "/Super_Twist/eigensymmint_" + std::to_string(spin_idx);
hdf5_file.load_data(this->symm_label_idxs[quantum_part_idx][spin_idx], hpath);
}
quantum_part_idx++;
}
} catch (...) { // failed to load symm info
APP_WARN(
"We are restarting a calculation with symmetry but symmetry info was not found in the bin file. Attempting to symmetrize orbitals... This is a numerically tenuous process, and may fail...");
this->symmetrize_orbitals(this->C_combined, this->symm_label_idxs, this->symm_labels);
}
}

// need to fill C into C_irrep
if (this->input_symmetry->do_symmetry == true) {
this->symmetrize_orbitals(this->C_combined, this->symm_label_idxs, this->symm_labels);
auto qpidx = 0;
for (auto const &[quantum_part_key, quantum_part] : this->input_molecule->quantum_particles) {
auto num_basis = this->input_basis->num_basis[quantum_part_idx];
auto spinidx = 0;
auto nmo = this->symm_label_idxs[qpidx][spinidx].size();
auto nirrep = this->C[qpidx][spinidx].size();
std::vector<int> populated_mos_per_irrep(nirrep, 0);
for (auto mo_idx = 0; mo_idx < nmo; mo_idx++) {
auto irrep_idx = this->symm_label_idxs[qpidx][0][mo_idx];
this->C[qpidx][spinidx][irrep_idx].col(populated_mos_per_irrep[irrep_idx]) = this->C_combined[qpidx][spinidx].col(mo_idx);
populated_mos_per_irrep[irrep_idx]++;
}
if (quantum_part.num_parts > 1 && quantum_part.restricted == false) {
auto num_basis = this->input_basis->num_basis[quantum_part_idx];
auto spinidx = 1;
auto nmo = this->symm_label_idxs[qpidx][spinidx].size();
auto nirrep = this->C[qpidx][spinidx].size();
std::vector<int> populated_mos_per_irrep(nirrep, 0);
for (auto mo_idx = 0; mo_idx < nmo; mo_idx++) {
auto irrep_idx = this->symm_label_idxs[qpidx][0][mo_idx];
this->C[qpidx][spinidx][irrep_idx].col(populated_mos_per_irrep[irrep_idx]) = this->C_combined[qpidx][spinidx].col(mo_idx);
populated_mos_per_irrep[irrep_idx]++;
}
}
}
}

this->form_occ();
Expand All @@ -1694,6 +1730,7 @@ void POLYQUANT_EPSCF::setup_from_file(std::string &filename) {
// Polyquant_cout("Running a single SCF iteration");
// this->run_iteration();
this->print_iteration();
this->calculate_E_elec();
this->calculate_E_total();
Polyquant_cout(this->E_total);
this->print_combined_orbitals("GUESS ORBITALS FROM FILE");
Expand Down