From 099d6774a4abb3b5fb82545dd5fc9f3306826019 Mon Sep 17 00:00:00 2001 From: jeffschriber Date: Fri, 19 Jan 2024 16:28:42 -0500 Subject: [PATCH 1/6] Generalize post-processing --- forte/CMakeLists.txt | 2 +- forte/api/forte_python_module.cc | 4 +-- .../{spin_corr.cc => post_process.cc} | 31 +++++++++++++------ .../{spin_corr.h => post_process.h} | 13 +++++--- forte/pymodule.py | 15 +++++++-- forte/register_forte_options.py | 2 ++ 6 files changed, 47 insertions(+), 20 deletions(-) rename forte/post_process/{spin_corr.cc => post_process.cc} (94%) rename forte/post_process/{spin_corr.h => post_process.h} (86%) diff --git a/forte/CMakeLists.txt b/forte/CMakeLists.txt index 3c20c0f69..a3a2270cc 100644 --- a/forte/CMakeLists.txt +++ b/forte/CMakeLists.txt @@ -228,7 +228,7 @@ orbital-helpers/semi_canonicalize.cc orbital-helpers/unpaired_density.cc pci/pci.cc pci/pci_sigma.cc -post_process/spin_corr.cc +post_process/post_process.cc sci/aci.cc sci/aci_build_F.cc sci/gasaci_build_F.cc diff --git a/forte/api/forte_python_module.cc b/forte/api/forte_python_module.cc index 4211c5ee5..ba18466f8 100644 --- a/forte/api/forte_python_module.cc +++ b/forte/api/forte_python_module.cc @@ -65,7 +65,7 @@ #include "sci/tdci.h" #include "genci/ci_occupation.h" -#include "post_process/spin_corr.h" +#include "post_process/post_process.h" namespace py = pybind11; using namespace pybind11::literals; @@ -217,7 +217,7 @@ PYBIND11_MODULE(_forte, m) { "Make an object that holds the molecular orbital integrals for the active orbitals"); m.def("make_dynamic_correlation_solver", &make_dynamic_correlation_solver, "Make a dynamical correlation solver"); - m.def("perform_spin_analysis", &perform_spin_analysis, "Do spin analysis"); + m.def("perform_post_processing", &perform_post_processing, "Do wave function post processing"); m.def("make_dsrg_method", &make_dsrg_method, "Make a DSRG method (spin-integrated implementation)"); m.def("make_sadsrg_method", &make_sadsrg_method, diff --git a/forte/post_process/spin_corr.cc b/forte/post_process/post_process.cc similarity index 94% rename from forte/post_process/spin_corr.cc rename to forte/post_process/post_process.cc index 787c0c4e0..93c871539 100644 --- a/forte/post_process/spin_corr.cc +++ b/forte/post_process/post_process.cc @@ -34,23 +34,34 @@ #include "base_classes/forte_options.h" #include "helpers/printing.h" #include "helpers/helpers.h" -#include "post_process/spin_corr.h" +#include "post_process/post_process.h" using namespace psi; namespace forte { -SpinCorr::SpinCorr(std::shared_ptr rdms, std::shared_ptr options, - std::shared_ptr mo_space_info, - std::shared_ptr as_ints) - : rdms_(rdms), options_(options), mo_space_info_(mo_space_info), as_ints_(as_ints) { +PostProcess::PostProcess(const std::string method, std::shared_ptr rdms, std::shared_ptr options, + std::shared_ptr mo_space_info, + std::shared_ptr as_ints) + : method_(method), rdms_(rdms), options_(options), mo_space_info_(mo_space_info), as_ints_(as_ints) { nactpi_ = mo_space_info_->dimension("ACTIVE"); nirrep_ = nactpi_.n(); nact_ = nactpi_.sum(); } -std::pair, std::shared_ptr> SpinCorr::compute_nos() { +void PostProcess::process() { + + if (method_ == "SPIN_CORRELATION"){ + spin_analysis(); + }else if (method_ == "UNPAIRED_DENSITY"){ + // up_density(); + } + +} + + +std::pair, std::shared_ptr> PostProcess::compute_nos() { print_h2("Natural Orbitals"); @@ -108,7 +119,7 @@ std::pair, std::shared_ptr> SpinCorr:: return std::make_pair(Ua, Ub); } -void SpinCorr::spin_analysis() { +void PostProcess::spin_analysis() { size_t nact = static_cast(nact_); size_t nact2 = nact * nact; size_t nact3 = nact * nact2; @@ -361,11 +372,11 @@ void SpinCorr::spin_analysis() { } } -void perform_spin_analysis(std::shared_ptr rdms, std::shared_ptr options, +void perform_post_processing(const std::string method, std::shared_ptr rdms, std::shared_ptr options, std::shared_ptr mo_space_info, std::shared_ptr as_ints) { - SpinCorr spin(rdms, options, mo_space_info, as_ints); - spin.spin_analysis(); + PostProcess proc(method, rdms, options, mo_space_info, as_ints); + proc.process(); } } // namespace forte diff --git a/forte/post_process/spin_corr.h b/forte/post_process/post_process.h similarity index 86% rename from forte/post_process/spin_corr.h rename to forte/post_process/post_process.h index b5d8be573..73b615902 100644 --- a/forte/post_process/spin_corr.h +++ b/forte/post_process/post_process.h @@ -41,21 +41,26 @@ namespace forte { /** - * @brief The SpinCorr class + * @brief The PostProcess class * This class computes the spin correlation from RDMs */ -class SpinCorr { +class PostProcess { public: - SpinCorr(std::shared_ptr rdms, std::shared_ptr options, + PostProcess(const std::string method, std::shared_ptr rdms, std::shared_ptr options, std::shared_ptr mo_space_info, std::shared_ptr as_ints); std::pair, std::shared_ptr> compute_nos(); + void process(); + void spin_analysis(); private: + + const std::string method_; + std::shared_ptr rdms_; std::shared_ptr options_; @@ -71,7 +76,7 @@ class SpinCorr { psi::Dimension nactpi_; }; -void perform_spin_analysis(std::shared_ptr rdms, std::shared_ptr options, +void perform_post_processing(const std::string method, std::shared_ptr rdms, std::shared_ptr options, std::shared_ptr mo_space_info, std::shared_ptr as_ints); diff --git a/forte/pymodule.py b/forte/pymodule.py index 023ce9a13..b8e2e6a5b 100644 --- a/forte/pymodule.py +++ b/forte/pymodule.py @@ -99,9 +99,9 @@ def forte_driver(data: ForteData): data = ActiveSpaceRDMs(max_rdm_level=max_rdm_level).run(data) write_external_rdm_file(data.rdms) - if options.get_bool("SPIN_ANALYSIS"): - data = ActiveSpaceRDMs(max_rdm_level=2, rdms_type=forte.RDMsType.spin_dependent).run(data) - forte.perform_spin_analysis(data.rdms, options, mo_space_info, as_ints) +# if options.get_bool("SPIN_ANALYSIS"): +# data = ActiveSpaceRDMs(max_rdm_level=2, rdms_type=forte.RDMsType.spin_dependent).run(data) +# forte.perform_spin_analysis(data.rdms, options, mo_space_info, as_ints) # solver for dynamical correlation from DSRG correlation_solver_type = options.get_str("CORRELATION_SOLVER") @@ -126,6 +126,15 @@ def forte_driver(data: ForteData): average_energy = forte.compute_average_state_energy(state_energies_list, state_weights_map) return_en = average_energy + + for proc in options.get_list("POST_PROCESS"): + if proc in ['SPIN_CORRELATION', 'UNPAIRED_DENSITY']: + data = ActiveSpaceRDMs(max_rdm_level=2, rdms_type=forte.RDMsType.spin_dependent).run(data) + forte.perform_post_processing(proc, data.rdms, options, mo_space_info, as_ints) + else: + raise Exception(f"Forte: Requestion post processing routine ({proc}) not implemented") + + return return_en diff --git a/forte/register_forte_options.py b/forte/register_forte_options.py index 769f5d76b..77b734053 100644 --- a/forte/register_forte_options.py +++ b/forte/register_forte_options.py @@ -536,6 +536,8 @@ def register_aci_options(options): options.add_bool("SPIN_ANALYSIS", False, "Do spin correlation analysis?") + options.add_list("POST_PROCESS", "List of post-processing items") + options.add_bool("SPIN_TEST", False, "Do test validity of correlation analysis") options.add_bool("ACI_RELAXED_SPIN", False, "Do spin correlation analysis for relaxed wave function?") From ba21c23756ebdf3f9b514f5223b85c2eea3ef12c Mon Sep 17 00:00:00 2001 From: jeffschriber Date: Fri, 19 Jan 2024 16:52:50 -0500 Subject: [PATCH 2/6] Start migrating unpaired density to post process --- forte/post_process/post_process.cc | 17 +++++++++++++++-- forte/post_process/post_process.h | 2 ++ .../unpaired_density.cc | 0 .../unpaired_density.h | 0 4 files changed, 17 insertions(+), 2 deletions(-) rename forte/{orbital-helpers => post_process}/unpaired_density.cc (100%) rename forte/{orbital-helpers => post_process}/unpaired_density.h (100%) diff --git a/forte/post_process/post_process.cc b/forte/post_process/post_process.cc index 93c871539..ff5da7a6b 100644 --- a/forte/post_process/post_process.cc +++ b/forte/post_process/post_process.cc @@ -61,7 +61,7 @@ void PostProcess::process() { } -std::pair, std::shared_ptr> PostProcess::compute_nos() { +std::pair, std::shared_ptr>, std::ppair,std::shared_ptr> PostProcess::compute_nos() { print_h2("Natural Orbitals"); @@ -116,7 +116,7 @@ std::pair, std::shared_ptr> PostProces } } } - return std::make_pair(Ua, Ub); + return std::make_pair>; } void PostProcess::spin_analysis() { @@ -372,6 +372,19 @@ void PostProcess::spin_analysis() { } } +void PostProcess::unpaired_density() { + + size_t nact = static_cast(nact_); + size_t nact2 = nact * nact; + size_t nact3 = nact * nact2; + + auto UA = std::make_shared(nact, nact); + auto UB = std::make_shared(nact, nact); + + // Get natural orbitals +} + + void perform_post_processing(const std::string method, std::shared_ptr rdms, std::shared_ptr options, std::shared_ptr mo_space_info, std::shared_ptr as_ints) { diff --git a/forte/post_process/post_process.h b/forte/post_process/post_process.h index 73b615902..e873adaf8 100644 --- a/forte/post_process/post_process.h +++ b/forte/post_process/post_process.h @@ -57,6 +57,8 @@ class PostProcess { void spin_analysis(); + void unpaired_density(); + private: const std::string method_; diff --git a/forte/orbital-helpers/unpaired_density.cc b/forte/post_process/unpaired_density.cc similarity index 100% rename from forte/orbital-helpers/unpaired_density.cc rename to forte/post_process/unpaired_density.cc diff --git a/forte/orbital-helpers/unpaired_density.h b/forte/post_process/unpaired_density.h similarity index 100% rename from forte/orbital-helpers/unpaired_density.h rename to forte/post_process/unpaired_density.h From c20de31bcc59661cdd810df12e8c95c6544f33ee Mon Sep 17 00:00:00 2001 From: jeffschriber Date: Thu, 25 Jan 2024 11:47:41 -0500 Subject: [PATCH 3/6] Unpaired electron post-processing working --- forte/CMakeLists.txt | 1 - forte/post_process/post_process.cc | 104 ++++---- forte/post_process/post_process.h | 6 +- forte/post_process/unpaired_density.cc | 344 ------------------------- forte/post_process/unpaired_density.h | 58 ----- forte/pymodule.py | 4 +- 6 files changed, 53 insertions(+), 464 deletions(-) delete mode 100644 forte/post_process/unpaired_density.cc delete mode 100644 forte/post_process/unpaired_density.h diff --git a/forte/CMakeLists.txt b/forte/CMakeLists.txt index a3a2270cc..3350d93a7 100644 --- a/forte/CMakeLists.txt +++ b/forte/CMakeLists.txt @@ -225,7 +225,6 @@ orbital-helpers/orbital_embedding.cc orbital-helpers/orbitaloptimizer.cc orbital-helpers/pao_builder.cc orbital-helpers/semi_canonicalize.cc -orbital-helpers/unpaired_density.cc pci/pci.cc pci/pci_sigma.cc post_process/post_process.cc diff --git a/forte/post_process/post_process.cc b/forte/post_process/post_process.cc index ff5da7a6b..1badcc61c 100644 --- a/forte/post_process/post_process.cc +++ b/forte/post_process/post_process.cc @@ -42,12 +42,16 @@ namespace forte { PostProcess::PostProcess(const std::string method, std::shared_ptr rdms, std::shared_ptr options, std::shared_ptr mo_space_info, + std::shared_ptr ints, std::shared_ptr as_ints) - : method_(method), rdms_(rdms), options_(options), mo_space_info_(mo_space_info), as_ints_(as_ints) { + : method_(method), rdms_(rdms), options_(options), mo_space_info_(mo_space_info), ints_(ints), as_ints_(as_ints) { nactpi_ = mo_space_info_->dimension("ACTIVE"); nirrep_ = nactpi_.n(); nact_ = nactpi_.sum(); + + print_h1("Wavefunction Post-processing"); + } void PostProcess::process() { @@ -55,13 +59,13 @@ void PostProcess::process() { if (method_ == "SPIN_CORRELATION"){ spin_analysis(); }else if (method_ == "UNPAIRED_DENSITY"){ - // up_density(); + unpaired_density(); } } -std::pair, std::shared_ptr>, std::ppair,std::shared_ptr> PostProcess::compute_nos() { +std::tuple, std::shared_ptr, std::shared_ptr,std::shared_ptr> PostProcess::compute_active_nos() { print_h2("Natural Orbitals"); @@ -116,7 +120,7 @@ std::pair, std::shared_ptr>, } } } - return std::make_pair>; + return std::make_tuple(Ua, Ub, OCC_A,OCC_B); } void PostProcess::spin_analysis() { @@ -127,60 +131,11 @@ void PostProcess::spin_analysis() { auto UA = std::make_shared(nact, nact); auto UB = std::make_shared(nact, nact); - // if (options_->get_str("SPIN_BASIS") == "IAO") { - // outfile->Printf("\n Computing spin correlation in IAO basis \n"); - // auto Ca = ints_->Ca(); - // std::shared_ptr IAO = - // IAOBuilder::build(reference_wavefunction_->basisset(), - // reference_wavefunction_->get_basisset("MINAO_BASIS"), Ca, - // options_->; - // outfile->Printf("\n Computing IAOs\n"); - // std::map> iao_info = IAO->build_iaos(); - // std::shared_ptr iao_orbs(iao_info["A"]->clone()); - - // std::shared_ptr Cainv(Ca->clone()); - // Cainv->invert(); - // auto iao_coeffs = psi::Matrix::doublet(Cainv, iao_orbs, false, - // false); - - // size_t new_dim = iao_orbs->colspi()[0]; - - // auto labels = IAO->print_IAO(iao_orbs, new_dim, nmo_, reference_wavefunction_); - // std::vector IAO_inds; - // if (options_->get_bool("PI_ACTIVE_SPACE")) { - // for (size_t i = 0, maxi = labels.size(); i < maxi; ++i) { - // std::string label = labels[i]; - // if (label.find("z") != std::string::npos) { - // IAO_inds.push_back(i); - // } - // } - // } else { - // nact = new_dim; - // for (size_t i = 0; i < new_dim; ++i) { - // IAO_inds.push_back(i); - // } - // } - - // std::vector active_mo = mo_space_info_->get_absolute_mo("ACTIVE"); - // for (size_t i = 0; i < nact; ++i) { - // int idx = IAO_inds[i]; - // outfile->Printf("\n Using IAO %d", idx); - // for (size_t j = 0; j < nact; ++j) { - // int mo = active_mo[j]; - // UA->set(j, i, iao_coeffs->get(mo, idx)); - // } - // } - // UB->copy(UA); - // outfile->Printf("\n"); - - //} else if (options_->get_str("SPIN_BASIS") == "NO") { outfile->Printf("\n Computing spin correlation in NO basis \n"); - auto pair = compute_nos(); - - UA = pair.first; - UB = pair.second; + auto no_tuple = compute_active_nos(); + std::tie(UA, UB,std::ignore,std::ignore) = no_tuple; int nmo = mo_space_info_->size("ALL"); @@ -374,21 +329,54 @@ void PostProcess::spin_analysis() { void PostProcess::unpaired_density() { + outfile->Printf("\n Computing unpaired electron numbers/density in NO basis."); size_t nact = static_cast(nact_); - size_t nact2 = nact * nact; - size_t nact3 = nact * nact2; auto UA = std::make_shared(nact, nact); auto UB = std::make_shared(nact, nact); + auto OCC_A = std::make_shared(nact); + auto OCC_B = std::make_shared(nact); // Get natural orbitals + auto no_tuple = compute_active_nos(); + std::tie(UA, UB,OCC_A,OCC_B) = no_tuple; + + std::vector nus(nact); + + outfile->Printf("\n N.O. nu_i "); + outfile->Printf("\n -------- --------"); + + // the total unpaired electrons (alpha only) + double nu_total = 0.0; + for (size_t i = 0; i < nact; ++i){ + double ni = OCC_A->get(i);// + OCC_B->get(i); + double nu_i = (ni*ni)*(1.0-ni)*(1.0-ni); + nus[i] = nu_i; + nu_total += nu_i; + outfile->Printf("\n %2d \t %7.5f ", i, nu_i); + } + + outfile->Printf("\n total: %5.3f", nu_total); + + // Plotting of unpaired densities should be handled externally + // Here, I'll put the unpaired numbers in a file and + // transform the orbitals, so the plots can be made + std::ofstream ofile; + ofile.open("unpaired_electrons_alpha.txt", std::ofstream::out | std::ofstream::trunc); + for (size_t i = 0; i < nact; ++i) { + ofile << std::setprecision(12) << nus[i] << "\n"; + } + ofile.close(); + + ints_->rotate_orbitals(UA,UB); } void perform_post_processing(const std::string method, std::shared_ptr rdms, std::shared_ptr options, std::shared_ptr mo_space_info, + std::shared_ptr ints, std::shared_ptr as_ints) { - PostProcess proc(method, rdms, options, mo_space_info, as_ints); + PostProcess proc(method, rdms, options, mo_space_info,ints, as_ints); proc.process(); } diff --git a/forte/post_process/post_process.h b/forte/post_process/post_process.h index e873adaf8..39a5cb1f6 100644 --- a/forte/post_process/post_process.h +++ b/forte/post_process/post_process.h @@ -49,9 +49,10 @@ class PostProcess { public: PostProcess(const std::string method, std::shared_ptr rdms, std::shared_ptr options, std::shared_ptr mo_space_info, + std::shared_ptr ints, std::shared_ptr as_ints); - std::pair, std::shared_ptr> compute_nos(); + std::tuple, std::shared_ptr, std::shared_ptr, std::shared_ptr> compute_active_nos(); void process(); @@ -69,6 +70,8 @@ class PostProcess { std::shared_ptr mo_space_info_; + std::shared_ptr ints_; + std::shared_ptr as_ints_; size_t nact_; @@ -80,6 +83,7 @@ class PostProcess { void perform_post_processing(const std::string method, std::shared_ptr rdms, std::shared_ptr options, std::shared_ptr mo_space_info, + std::shared_ptr ints, std::shared_ptr as_ints); } // namespace forte diff --git a/forte/post_process/unpaired_density.cc b/forte/post_process/unpaired_density.cc deleted file mode 100644 index b8a26c8c6..000000000 --- a/forte/post_process/unpaired_density.cc +++ /dev/null @@ -1,344 +0,0 @@ -/* - * @BEGIN LICENSE - * - * Forte: an open-source plugin to Psi4 (https://github.com/psi4/psi4) - * that implements a variety of quantum chemistry methods for strongly - * correlated electrons. - * - * Copyright (c) 2012-2024 by its authors (see COPYING, COPYING.LESSER, AUTHORS). - * - * The copyrights for code used from other parties are included in - * the corresponding files. - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU Lesser General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public License - * along with this program. If not, see http://www.gnu.org/licenses/. - * - * @END LICENSE - */ - -#include "psi4/libmints/molecule.h" -#include "psi4/libpsio/psio.h" -#include "psi4/libpsio/psio.hpp" - -#include "psi4/psi4-dec.h" -#include "psi4/libpsi4util/PsiOutStream.h" - -#include "psi4/libmints/local.h" -#include "psi4/libmints/matrix.h" -#include "psi4/libmints/vector.h" - -#include "base_classes/forte_options.h" -#include "unpaired_density.h" -#include "localize.h" - -using namespace psi; - -namespace forte { - -UPDensity::UPDensity(std::shared_ptr ints, - std::shared_ptr mo_space_info, - std::shared_ptr options, std::shared_ptr Ua, - std::shared_ptr Ub) - : options_(options), ints_(ints), mo_space_info_(mo_space_info), Uas_(Ua), Ubs_(Ub) {} - -void UPDensity::compute_unpaired_density(std::vector& oprdm_a, - std::vector& oprdm_b) { - // TODO: re-enable this code - // psi::Dimension nactpi = mo_space_info_->dimension("ACTIVE"); - // psi::Dimension nmopi = mo_space_info_->dimension("ALL"); - // psi::Dimension ncmopi = mo_space_info_->dimension("CORRELATED"); - // size_t nirrep = ints_->nirrep(); - // psi::Dimension rdocc = mo_space_info_->dimension("RESTRICTED_DOCC"); - // psi::Dimension fdocc = mo_space_info_->dimension("FROZEN_DOCC"); - // - // size_t nact = nactpi.sum(); - // - // // First compute natural orbitals - // auto opdm_a = std::make_shared("OPDM_A", nirrep, nactpi, nactpi); - // auto opdm_b = std::make_shared("OPDM_B", nirrep, nactpi, nactpi); - // - // // Put 1-RDM into Shared matrix - // int offset = 0; - // for (size_t h = 0; h < nirrep; ++h) { - // for (int u = 0; u < nactpi[h]; ++u) { - // for (int v = 0; v < nactpi[h]; ++v) { - // opdm_a->set(h, u, v, oprdm_a[(u + offset) * nact + v + offset]); - // opdm_b->set(h, u, v, oprdm_b[(u + offset) * nact + v + offset]); - // } - // } - // offset += nactpi[h]; - // } - // // opdm_a->transform(Uas_); - // // opdm_b->transform(Ubs_); - // - // // Diagonalize the 1-RDMs - // auto OCC_A = std::make_shared("ALPHA NOCC", nirrep, nactpi) - // auto OCC_B = std::make_shared("BETA NOCC", nirrep, nactpi); - // auto NO_A = std::make_shared(nirrep, nactpi, nactpi); - // auto NO_B = std::make_shared(nirrep, nactpi, nactpi); - // - // opdm_a->diagonalize(NO_A, OCC_A, descending); - // opdm_b->diagonalize(NO_B, OCC_B, descending); - // - // // Build the transformation matrix - // // Only build density for active orbitals - // auto Ua = std::make_shared("Ua", nmopi, nmopi); - // auto Ub = std::make_shared("Ub", nmopi, nmopi); - // - // Ua->zero(); - // Ub->zero(); - // - // // This Ua/Ub build will ensure that the density only includes active orbitals - // // If natural orbitals are desired, change the 1.0 to NO_a->get(p,q) - // for (size_t h = 0; h < nirrep; ++h) { - // size_t irrep_offset = fdocc[h] + rdocc[h]; - // for (int p = 0; p < nactpi[h]; ++p) { - // // for (int q = 0; q < nactpi[h]; ++q) { - // Ua->set(h, p + irrep_offset, p + irrep_offset, 1.0); - // Ub->set(h, p + irrep_offset, p + irrep_offset, 1.0); - // // } - // } - // } - // - // /// ** Compute atom-based unpaired contributions - // - // // ** This will be done in a completely localized basis - // // ** This code has only been tested for pz (pi) orbitals, beware! - // - // auto Ua_act = std::make_shared(nact, nact); - // - // // relocalize to atoms - // - // // Grab matrix that takes the transforms from the NO basis to our local basis - // auto loc = std::make_shared(options_, ints_, mo_space_info_); - // loc->full_localize(); - // Ua_act = loc->get_U()->clone(); - // std::shared_ptr Noinv(NO_A->clone()); - // Noinv->invert(); - // auto Ua_act_r = psi::linalg::doublet(Noinv, Ua_act, false, false); - // - // // Compute sum(p,i) n_i * ( 1 - n_i ) * (U_p,i)^2 - // double total = 0.0; - // std::vector scales(nact); - // for (size_t i = 0; i < nact; ++i) { - // double value = 0.0; - // for (size_t h = 0; h < nirrep; ++h) { - // for (int p = 0; p < nactpi[h]; ++p) { - // // double n_p = OCC_A->get(p) + OCC_B->get(p); - // double n_p = OCC_A->get(p); - // double up_el = n_p * (1.0 - n_p); - // - // value += up_el * Ua_act_r->get(p, i) * Ua_act_r->get(p, i); - // } - // } - // scales[i] = value; - // total += value; - // outfile->Printf("\n MO %d: %1.6f", i, value); - // } - // outfile->Printf("\n Total unpaired electrons: %1.4f", total); - // - // // Build the density using scaled columns of C - // - // auto Ca = ints_->Ca(); - // auto Cb = ints_->Cb(); - // - // auto Ca_new = psi::linalg::doublet(Ca->clone(), Ua, false, false); - // auto Cb_new = psi::linalg::doublet(Cb->clone(), Ub, false, false); - // - // for (size_t h = 0; h < nirrep; ++h) { - // int offset = fdocc[h] + rdocc[h]; - // for (int p = 0; p < nactpi[h]; ++p) { - // // double n_p = OCC_A->get(p) + OCC_B->get(p); - // // double up_el = n_p * (2.0 - n_p); - // double up_el = scales[p]; - // // double up_el = n_p * n_p * (2.0 - n_p ) * (2.0 - n_p); - // // outfile->Printf("\n Weight for orbital (%d,%d): %1.5f", h, p, - // up_el); Ca_new->scale_column(h, offset + p, up_el); Cb_new->scale_column(h, offset - // + p, up_el); - // } - // } - - psi::Dimension nactpi = mo_space_info_->dimension("ACTIVE"); - psi::Dimension nmopi = mo_space_info_->dimension("ALL"); - psi::Dimension ncmopi = mo_space_info_->dimension("CORRELATED"); - size_t nirrep = ints_->nirrep(); - psi::Dimension rdocc = mo_space_info_->dimension("RESTRICTED_DOCC"); - psi::Dimension fdocc = mo_space_info_->dimension("FROZEN_DOCC"); - - size_t nact = nactpi.sum(); - - // First compute natural orbitals - auto opdm_a = std::make_shared("OPDM_A", nirrep, nactpi, nactpi); - auto opdm_b = std::make_shared("OPDM_B", nirrep, nactpi, nactpi); - - // Put 1-RDM into Shared matrix - int offset = 0; - for (size_t h = 0; h < nirrep; ++h) { - for (int u = 0; u < nactpi[h]; ++u) { - for (int v = 0; v < nactpi[h]; ++v) { - opdm_a->set(h, u, v, oprdm_a[(u + offset) * nact + v + offset]); - opdm_b->set(h, u, v, oprdm_b[(u + offset) * nact + v + offset]); - } - } - offset += nactpi[h]; - } - // opdm_a->transform(Uas_); - // opdm_b->transform(Ubs_); - - // Diagonalize the 1-RDMs - auto OCC_A = std::make_shared("ALPHA NOCC", nactpi); - auto OCC_B = std::make_shared("BETA NOCC", nactpi); - auto NO_A = std::make_shared(nirrep, nactpi, nactpi); - auto NO_B = std::make_shared(nirrep, nactpi, nactpi); - - opdm_a->diagonalize(NO_A, OCC_A, descending); - opdm_b->diagonalize(NO_B, OCC_B, descending); - - // Build the transformation matrix - // Only build density for active orbitals - auto Ua = std::make_shared("Ua", nmopi, nmopi); - auto Ub = std::make_shared("Ub", nmopi, nmopi); - - Ua->zero(); - Ub->zero(); - - // This Ua/Ub build will ensure that the density only includes active orbitals - // If natural orbitals are desired, change the 1.0 to NO_a->get(p,q) - for (size_t h = 0; h < nirrep; ++h) { - size_t irrep_offset = fdocc[h] + rdocc[h]; - for (int p = 0; p < nactpi[h]; ++p) { - // for (int q = 0; q < nactpi[h]; ++q) { - Ua->set(h, p + irrep_offset, p + irrep_offset, 1.0); - Ub->set(h, p + irrep_offset, p + irrep_offset, 1.0); - // } - } - } - - /// ** Compute atom-based unpaired contributions - - // ** This will be done in a completely localized basis - // ** This code has only been tested for pz (pi) orbitals, beware! - - auto Ua_act = std::make_shared(nact, nact); - - // relocalize to atoms - - // Grab matrix that takes the transforms from the NO basis to our local basis - auto loc = std::make_shared(options_, ints_, mo_space_info_); - - std::vector actmo = mo_space_info_->absolute_mo("ACTIVE"); - std::vector loc_mo(2); - loc_mo[0] = static_cast(actmo[0]); - loc_mo[1] = static_cast(actmo.back()); - loc->set_orbital_space(loc_mo); - - loc->compute_transformation(); - Ua_act = loc->get_Ua()->clone(); - std::shared_ptr Noinv(NO_A->clone()); - Noinv->invert(); - auto Ua_act_r = psi::linalg::doublet(Noinv, Ua_act, false, false); - - // Compute sum(p,i) n_i * ( 1 - n_i ) * (U_p,i)^2 - double total = 0.0; - std::vector scales(nact); - for (size_t i = 0; i < nact; ++i) { - double value = 0.0; - for (size_t h = 0; h < nirrep; ++h) { - for (int p = 0; p < nactpi[h]; ++p) { - // double n_p = OCC_A->get(p) + OCC_B->get(p); - double n_p = OCC_A->get(p); - double up_el = n_p * (1.0 - n_p); - - value += up_el * Ua_act_r->get(p, i) * Ua_act_r->get(p, i); - } - } - scales[i] = value; - total += value; - outfile->Printf("\n MO %d: %1.6f", i, value); - } - outfile->Printf("\n Total unpaired electrons: %1.4f", total); - - // Build the density using scaled columns of C - - auto Ca = ints_->Ca(); - auto Cb = ints_->Cb(); - - auto Ca_new = psi::linalg::doublet(Ca->clone(), Ua, false, false); - auto Cb_new = psi::linalg::doublet(Cb->clone(), Ub, false, false); - - for (size_t h = 0; h < nirrep; ++h) { - int offset = fdocc[h] + rdocc[h]; - for (int p = 0; p < nactpi[h]; ++p) { - // double n_p = OCC_A->get(p) + OCC_B->get(p); - // double up_el = n_p * (2.0 - n_p); - double up_el = scales[p]; - // double up_el = n_p * n_p * (2.0 - n_p ) * (2.0 - n_p); - // outfile->Printf("\n Weight for orbital (%d,%d): %1.5f", h, p, up_el); - Ca_new->scale_column(h, offset + p, up_el); - Cb_new->scale_column(h, offset + p, up_el); - } - } - - auto Da = ints_->wfn()->Da(); - auto Db = ints_->wfn()->Db(); - - // auto Da_new = std::make_shared("Da_new", nmopi, nmopi); - // auto Db_new = std::make_shared("Db_new", nmopi, nmopi); - - // Da_new->gemm(false, true, 1.0, Ca_new, Ca_new, 0.0); - // Db_new->gemm(false, true, 1.0, Cb_new, Cb_new, 0.0); - - // Da->copy(Da_new); - // Db->copy(Db_new); - - // This is for IAOs, don't really need it - // std::shared_ptr IAO = - // IAOBuilder::build(wfn_->basisset(), - // wfn_->get_basisset("MINAO_BASIS"), Ca, options_); - // outfile->Printf("\n Computing IAOs\n"); - // std::map> iao_info = IAO->build_iaos(); - // std::shared_ptr iao_orbs(iao_info["A"]->clone()); - // - // std::shared_ptr Cainv(Ca->clone()); - // Cainv->invert(); - // auto iao_coeffs = psi::linalg::doublet(Cainv, iao_orbs, false, - // false); - // - // size_t new_dim = iao_orbs->colspi()[0]; - // size_t new_dim2 = new_dim * new_dim; - // size_t new_dim3 = new_dim2 * new_dim; - // - // auto labels = IAO->print_IAO(iao_orbs, new_dim, nmo, wfn_); - // - // outfile->Printf("\n label size: %zu", labels.size()); - // - // std::vector IAO_inds; - // for (int i = 0; i < labels.size(); ++i) { - // std::string label = labels[i]; - // if (label.find("z") != std::string::npos) { - // IAO_inds.push_back(i); - // } - // } - // std::vector active_mo = mo_space_info_->absolute_mo("ACTIVE"); - // for (int i = 0; i < nact; ++i) { - // int idx = IAO_inds[i]; - // outfile->Printf("\n Using IAO %d", idx); - // for (int j = 0; j < nact; ++j) { - // int mo = active_mo[j]; - // Ua_act->set(j, i, iao_coeffs->get(mo, idx)); - // } - // } -} - -UPDensity::~UPDensity() {} -} // namespace forte diff --git a/forte/post_process/unpaired_density.h b/forte/post_process/unpaired_density.h deleted file mode 100644 index 8426f109e..000000000 --- a/forte/post_process/unpaired_density.h +++ /dev/null @@ -1,58 +0,0 @@ -/* - * @BEGIN LICENSE - * - * Forte: an open-source plugin to Psi4 (https://github.com/psi4/psi4) - * that implements a variety of quantum chemistry methods for strongly - * correlated electrons. - * - * Copyright (c) 2012-2024 by its authors (see COPYING, COPYING.LESSER, AUTHORS). - * - * The copyrights for code used from other parties are included in - * the corresponding files. - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU Lesser General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public License - * along with this program. If not, see http://www.gnu.org/licenses/. - * - * @END LICENSE - */ - -#pragma once - -#include "psi4/libmints/wavefunction.h" -#include "psi4/libmints/local.h" - -#include "base_classes/mo_space_info.h" -#include "integrals/integrals.h" -#include "base_classes/rdms.h" -#include "iao_builder.h" - -namespace forte { - -class UPDensity { - public: - UPDensity(std::shared_ptr ints, std::shared_ptr mo_space_info, - std::shared_ptr options, std::shared_ptr Ua, - std::shared_ptr Ub); - - ~UPDensity(); - - void compute_unpaired_density(std::vector& ordm_a, std::vector& ordm_b); - - private: - std::shared_ptr options_; - std::shared_ptr ints_; - std::shared_ptr mo_space_info_; - std::shared_ptr Uas_; - std::shared_ptr Ubs_; -}; -} // namespace forte diff --git a/forte/pymodule.py b/forte/pymodule.py index b8e2e6a5b..fb5ac1d4f 100644 --- a/forte/pymodule.py +++ b/forte/pymodule.py @@ -130,9 +130,9 @@ def forte_driver(data: ForteData): for proc in options.get_list("POST_PROCESS"): if proc in ['SPIN_CORRELATION', 'UNPAIRED_DENSITY']: data = ActiveSpaceRDMs(max_rdm_level=2, rdms_type=forte.RDMsType.spin_dependent).run(data) - forte.perform_post_processing(proc, data.rdms, options, mo_space_info, as_ints) + forte.perform_post_processing(proc, data.rdms, options, mo_space_info,ints, as_ints) else: - raise Exception(f"Forte: Requestion post processing routine ({proc}) not implemented") + raise Exception(f"Forte: Requested post processing routine ({proc}) not implemented") return return_en From bd77482489e4aeae7bc1acf5c284ab233279ffe0 Mon Sep 17 00:00:00 2001 From: jeffschriber Date: Thu, 25 Jan 2024 14:55:36 -0500 Subject: [PATCH 4/6] Remove some comments --- forte/pymodule.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/forte/pymodule.py b/forte/pymodule.py index fb5ac1d4f..2bd73fd2b 100644 --- a/forte/pymodule.py +++ b/forte/pymodule.py @@ -99,10 +99,6 @@ def forte_driver(data: ForteData): data = ActiveSpaceRDMs(max_rdm_level=max_rdm_level).run(data) write_external_rdm_file(data.rdms) -# if options.get_bool("SPIN_ANALYSIS"): -# data = ActiveSpaceRDMs(max_rdm_level=2, rdms_type=forte.RDMsType.spin_dependent).run(data) -# forte.perform_spin_analysis(data.rdms, options, mo_space_info, as_ints) - # solver for dynamical correlation from DSRG correlation_solver_type = options.get_str("CORRELATION_SOLVER") if correlation_solver_type != "NONE": From 4a08c32ca40c9bb9eb7600da8c2aa995040da928 Mon Sep 17 00:00:00 2001 From: jeffschriber Date: Thu, 22 Feb 2024 15:15:23 -0500 Subject: [PATCH 5/6] add definition of unpaired electron consistent with our paper --- forte/post_process/post_process.cc | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/forte/post_process/post_process.cc b/forte/post_process/post_process.cc index 1badcc61c..b22119eaf 100644 --- a/forte/post_process/post_process.cc +++ b/forte/post_process/post_process.cc @@ -343,17 +343,18 @@ void PostProcess::unpaired_density() { std::vector nus(nact); - outfile->Printf("\n N.O. nu_i "); - outfile->Printf("\n -------- --------"); + outfile->Printf("\n N.O. occ nu_i "); + outfile->Printf("\n -------- --------- --------"); // the total unpaired electrons (alpha only) double nu_total = 0.0; for (size_t i = 0; i < nact; ++i){ double ni = OCC_A->get(i);// + OCC_B->get(i); - double nu_i = (ni*ni)*(1.0-ni)*(1.0-ni); + //double nu_i = (ni*ni)*(1.0-ni)*(1.0-ni); + double nu_i = (ni)*(1.0-ni); nus[i] = nu_i; nu_total += nu_i; - outfile->Printf("\n %2d \t %7.5f ", i, nu_i); + outfile->Printf("\n %2d \t %7.5f \t %7.5f ", i, ni,nu_i); } outfile->Printf("\n total: %5.3f", nu_total); From df5ceedeb5d6c7709c74f36f07ee348b210117c2 Mon Sep 17 00:00:00 2001 From: jeffschriber Date: Mon, 17 Jun 2024 11:37:42 -0400 Subject: [PATCH 6/6] Unpaired electron computation now high-spin friendly --- forte/post_process/post_process.cc | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/forte/post_process/post_process.cc b/forte/post_process/post_process.cc index b22119eaf..55c69be76 100644 --- a/forte/post_process/post_process.cc +++ b/forte/post_process/post_process.cc @@ -343,6 +343,7 @@ void PostProcess::unpaired_density() { std::vector nus(nact); + outfile->Printf("\n Number of unpaired electrons (nu_i) in alpha NOs"); outfile->Printf("\n N.O. occ nu_i "); outfile->Printf("\n -------- --------- --------"); @@ -369,6 +370,10 @@ void PostProcess::unpaired_density() { } ofile.close(); + if (ints_->spin_restriction() == IntegralSpinRestriction::Restricted){ + UB = UA->clone(); + } + ints_->rotate_orbitals(UA,UB); }