Renaming of excited states dipoles (#2541)

* Renaming of excited states dipoles

* Update tests/cc48/input.dat

Co-authored-by: Lori A. Burns <lori.burns@gmail.com>

* UHF -> not same density

* Include ground state in correlation energy count

* Move densities to new field, in EOM

Co-authored-by: Lori A. Burns <lori.burns@gmail.com>

doc/sphinxman/source/glossary_psivariables.rst

@@ -127,13 +127,19 @@ PSI Variables by Alpha
    energy [Eh] and correlation correction components [Eh] for the compound
    method requested through cbs().
 
-.. psivar:: CC ROOT n DIPOLE
+.. psivar:: CCname ROOT n DIPOLE
+   CCname ROOT n (h) DIPOLE
+   CCname ROOT n DIPOLE - h TRANSITION
 
    Dipole array [e a0] for the requested coupled cluster level of theory and root *n* (number starts at GS = 0), (3,).
+   Conventions for root indexing and whether h refers to transition or root irrep are as in :ref:`sec:psivarnotes`.
 
-.. psivar:: CC ROOT n QUADRUPOLE
+.. psivar:: CCname ROOT n QUADRUPOLE
+   CCname ROOT n (h) QUADRUPOLE
+   CCname ROOT n QUADRUPOLE - h TRANSITION
 
    Redundant quadrupole array [e a0^2] for the requested coupled cluster level of theory and root *n* (number starts at GS = 0), (3, 3).
+   Conventions for root indexing and whether h refers to transition or root irrep are as in :ref:`sec:psivarnotes`.
 
 .. psivar:: CC ROOT n TOTAL ENERGY
    CC ROOT n CORRELATION ENERGY
@@ -1559,7 +1565,9 @@ PSI Variables by Alpha
    ADC ROOT 0 -> ROOT m OSCILLATOR STRENGTH (VEL)
    ADC ROOT 0 (h) -> ROOT m (i) OSCILLATOR STRENGTH (VEL)
    ADC ROOT 0 -> ROOT m OSCILLATOR STRENGTH (VEL) - h TRANSITION
-   CC ROOT n (h) -> ROOT m (i) OSCILLATOR STRENGTH (LEN)
+   CCname ROOT n -> ROOT m OSCILLATOR STRENGTH (LEN)
+   CCname ROOT n (h) -> ROOT m (i) OSCILLATOR STRENGTH (LEN)
+   CCname ROOT n -> ROOT m OSCILLATOR STRENGTH (LEN) - h TRANSITION
    TD-fctl ROOT 0 -> ROOT m OSCILLATOR STRENGTH (LEN)
    TD-fctl ROOT 0 (h) -> ROOT m (i) OSCILLATOR STRENGTH (LEN)
    TD-fctl ROOT 0 -> ROOT m OSCILLATOR STRENGTH (LEN) - h TRANSITION
@@ -1573,8 +1581,12 @@ PSI Variables by Alpha
    Conventions for root indexing and whether h refers to transition or root
    irrep are as in :ref:`sec:psivarnotes`.
 
-.. psivar:: CC ROOT n (h) -> ROOT m (i) ROTARY STRENGTH (LEN)
-   CC ROOT n (h) -> ROOT m (i) ROTARY STRENGTH (VEL)
+.. psivar:: CCname ROOT n -> ROOT m ROTARY STRENGTH (LEN)
+   CCname ROOT n (h) -> ROOT m (i) ROTARY STRENGTH (LEN)
+   CCname ROOT n -> ROOT m ROTARY STRENGTH (LEN) - h TRANSITION
+   CCname ROOT n -> ROOT m ROTARY STRENGTH (VEL)
+   CCname ROOT n (h) -> ROOT m (i) ROTARY STRENGTH (VEL)
+   CCname ROOT n -> ROOT m ROTARY STRENGTH (VEL) - h TRANSITION
    TD-fctl ROOT 0 -> ROOT m ROTARY STRENGTH (LEN)
    TD-fctl ROOT 0 (h) -> ROOT m (i) ROTARY STRENGTH (LEN)
    TD-fctl ROOT 0 -> ROOT m ROTARY STRENGTH (LEN) - h TRANSITION

psi4/driver/p4util/python_helpers.py

@@ -678,7 +678,6 @@ def _qcvar_reshape_set(key, val):
     """Reverse `_qcvar_reshape_get` for internal psi4.core.Matrix storage."""
 
     reshaper = None
-
     if key.upper().startswith("MBIS"):
         if key.upper().endswith("CHARGES"):
             return val
@@ -693,12 +692,13 @@ def _qcvar_reshape_set(key, val):
             val = val.reshape(-1, 3, 3, 3)
             val = np.array([_multipole_compressor(val[iat], 3) for iat in range(len(val))])
             return val
-    elif key.upper().endswith("DIPOLE"):
+    elif key.upper().endswith("DIPOLE") or "DIPOLE -" in key.upper():
         reshaper = (1, 3)
     elif "QUADRUPOLE POLARIZABILITY TENSOR" in key.upper():
         reshaper = (3, 3, 3)
-    elif any(key.upper().endswith(p) for p in _multipole_order):
-        val = _multipole_compressor(val, _multipole_order.index(key.upper().split()[-1]))
+    elif any((key.upper().endswith(p) or f"{p} -" in key.upper()) for p in _multipole_order):
+        p = [p for p in _multipole_order if (key.upper().endswith(p) or f"{p} -" in key.upper())]
+        val = _multipole_compressor(val, _multipole_order.index(p[0]))
         reshaper = (1, -1)
     elif key.upper() in ["MULLIKEN_CHARGES", "LOWDIN_CHARGES", "MULLIKEN CHARGES", "LOWDIN CHARGES"]:
         reshaper = (1, -1)
@@ -727,12 +727,13 @@ def _qcvar_reshape_get(key, val):
             val = val.np.reshape(-1, 10)
             val = np.array([_multipole_plumper(val[iat], 3) for iat in range(len(val))])
             return val
-    elif key.upper().endswith("DIPOLE"):
+    elif key.upper().endswith("DIPOLE") or "DIPOLE -" in key.upper():
         reshaper = (3, )
     elif "QUADRUPOLE POLARIZABILITY TENSOR" in key.upper():
         reshaper = (3, 3, 3)
-    elif any(key.upper().endswith(p) for p in _multipole_order):
-        return _multipole_plumper(val.np.reshape((-1, )), _multipole_order.index(key.upper().split()[-1]))
+    elif any((key.upper().endswith(p) or f"{p} -" in key.upper()) for p in _multipole_order):
+        p = [p for p in _multipole_order if (key.upper().endswith(p) or f"{p} -" in key.upper())]
+        return _multipole_plumper(val.np.reshape((-1, )), _multipole_order.index(p[0]))
     elif key.upper() in ["MULLIKEN_CHARGES", "LOWDIN_CHARGES", "MULLIKEN CHARGES", "LOWDIN CHARGES"]:
         reshaper = (-1, )
 

psi4/driver/procrouting/proc.py

@@ -3165,36 +3165,56 @@ def run_cc_property(name, **kwargs):
         core.cclambda(ccwfn)
         core.ccdensity(ccwfn)
 
+    # => Make OEProp calls <=
     if n_one > 0:
-        # call oe prop for GS density
+        # ==> Initialize OEProp  <==
         oe = core.OEProp(ccwfn)
+        for oe_prop_name in one:
+            oe.add(oe_prop_name.upper())
+        # ==> OEProp for the ground state <==
         # TODO: When Psi is Py 3.9+, transition to the removeprefix version.
         title = name.upper().replace("EOM-", "")
         #title = name.upper().removeprefix("EOM-")
         oe.set_title(title)
-        oe.set_names({title + " {}", "CC {}"})
-        for oe_name in one:
-            oe.add(oe_name.upper())
+        set_of_names = {title + " {}", "CC {}"}
+        if name.startswith("eom"):
+            gs_h = 0
+            for h, i in enumerate(ccwfn.soccpi()):
+                if i % 2:
+                    gs_h = gs_h ^ h
+            ct = ccwfn.molecule().point_group().char_table()
+            total_h_lbl = ct.gamma(0).symbol()
+            gs_h_lbl = ct.gamma(gs_h).symbol()
+            set_of_names.update({title + " ROOT 0 {}", "CC ROOT 0 {}",
+                                 f"{title} ROOT 0 {{}} - {total_h_lbl} TRANSITION",
+                                      f"CC ROOT 0 {{}} - {total_h_lbl} TRANSITION",
+                                 f"{title} ROOT 0 ({gs_h_lbl}) {{}}", f"CC ROOT 0 ({gs_h_lbl}) {{}}"})
+        oe.set_names(set_of_names)
         oe.compute()
-        # call oe prop for each ES density
+        print(ccwfn.variables())
+        # ==> OEProp for Excited States <==
         if name.startswith('eom'):
-            # copy GS CC DIP/QUAD ... to CC ROOT 0 DIP/QUAD ... if we are doing multiple roots
-            if 'DIPOLE' in one:
-                core.set_variable("CC ROOT 0 DIPOLE", core.variable("CC DIPOLE"))
-                # core.set_variable("CC ROOT n DIPOLE", core.variable("CC DIPOLE"))  # P::e CCENERGY
-            if 'QUADRUPOLE' in one:
-                core.set_variable("CC ROOT 0 QUADRUPOLE", core.variable("CC QUADRUPOLE"))
-                # core.set_variable("CC ROOT n QUADRUPOLE", core.variable("CC QUADRUPOLE"))  # P::e CCENERGY
-
-            n_root = sum(core.get_global_option("ROOTS_PER_IRREP"))
-            for rn in range(n_root):
-                oe.set_title("CC ROOT {}".format(rn + 1))
-                Da = ccwfn.variable("CC ROOT {} Da".format(rn + 1))
-                oe.set_Da_so(Da)
-                if core.get_global_option("REFERENCE") == "UHF":
-                    Db = ccwfn.variable("CC ROOT {} Db".format(rn + 1))
-                    oe.set_Db_so(Db)
-                oe.compute()
+            n_root_pi = core.get_global_option("ROOTS_PER_IRREP")
+            for h in range(ccwfn.nirrep()):
+                root_h_lbl = ct.gamma(h).symbol()
+                trans_h_lbl = ct.gamma(h ^ gs_h).symbol()
+                # Don't forget to count the ground state!
+                for i in range(n_root_pi[h]):
+                    if h == gs_h: i += 1
+                    root_title = title + f" ROOT {i} ({root_h_lbl})"
+                    oe.set_title(root_title)
+                    total_idx = ccwfn.total_index(i, h)
+                    set_of_names = {f"{title} ROOT {total_idx} {{}}", f"CC ROOT {total_idx} {{}}",
+                                    f"{title} ROOT {total_idx} {{}} - {trans_h_lbl} TRANSITION",
+                                         f"CC ROOT {total_idx} {{}} - {trans_h_lbl} TRANSITION",
+                                    f"{title} ROOT {i} ({root_h_lbl}) {{}}", f"CC ROOT {i} ({root_h_lbl}) {{}}"}
+                    oe.set_names(set_of_names)
+                    Da = ccwfn.get_alpha_density(root_title + " Da")
+                    oe.set_Da_so(Da)
+                    if not ccwfn.same_a_b_dens():
+                        Db = ccwfn.get_beta_density(root_title + " Db")
+                        oe.set_Db_so(Db)
+                    oe.compute()
 
     core.set_global_option('WFN', 'SCF')
     core.revoke_global_option_changed('WFN')

psi4/src/export_wavefunction.cc

@@ -319,7 +319,9 @@ void export_wavefunction(py::module& m) {
         .def("set_PCM", &Wavefunction::set_PCM, "Set the PCM object")
         .def("get_PCM", &Wavefunction::get_PCM, "Get the PCM object")
 #endif
-        .def("PCM_enabled", &Wavefunction::PCM_enabled, "Whether running a PCM calculation");
+        .def("PCM_enabled", &Wavefunction::PCM_enabled, "Whether running a PCM calculation")
+        .def("get_alpha_density", [](Wavefunction& wfn, std::string name) {return wfn.Da_map_[name] ;}, "Experimental!")
+        .def("get_beta_density", [](Wavefunction& wfn, std::string name) {return wfn.Db_map_[name] ;}, "Experimental!");
 
     py::class_<scf::HF, std::shared_ptr<scf::HF>, Wavefunction>(m, "HF", "docstring")
         .def("compute_fvpi", &scf::HF::compute_fvpi, "Update number of frozen virtuals")
@@ -603,6 +605,8 @@ void export_wavefunction(py::module& m) {
     py::class_<ccenergy::CCEnergyWavefunction, std::shared_ptr<ccenergy::CCEnergyWavefunction>, Wavefunction>(
         m, "CCWavefunction", "Specialized Wavefunction used by the ccenergy, cceom, ccgradient, etc. modules.")
         .def(py::init<std::shared_ptr<Wavefunction>, Options&>())
+        .def("total_index", [](ccenergy::CCEnergyWavefunction& wfn, int i, int h) { return wfn.total_indices[{i, h}]; },
+                """Map an index (i) within irrep (h) to its energy-sorted index among all roots.", "i"_a, "h"_a)
         .def("get_amplitudes", &ccenergy::CCEnergyWavefunction::get_amplitudes, R"pbdoc(
                Get dict of converged T amplitudes.
 

psi4/src/psi4/cc/ccdensity/Params.h

@@ -57,6 +57,11 @@ struct Params {
     int transition;
     int dertype;
     std::string wfn;
+    // TODO: Make meaning of variable invariant.
+    //       Initially, this variable means number of states, but later
+    //       it changes to number of excited states. This inconsistency
+    //       should be fixed. Resolving the problem is not safe until the
+    //       ccdensity module is documented.
     int nstates;
     int prop_sym;
     int prop_root;

psi4/src/psi4/cc/ccdensity/ccdensity.cc

@@ -39,6 +39,7 @@
 #include <cstdlib>
 #include <cstring>
 #include <cmath>
+#include "ccdensity.h"
 #include "MOInfo.h"
 #include "Params.h"
 #include "Frozen.h"
@@ -376,6 +377,15 @@ PsiReturnType ccdensity(std::shared_ptr<ccenergy::CCEnergyWavefunction> ref_wfn,
             Pb_x = block_to_matrix(moinfo.opdm_b);
         }
 
+        std::string short_name;
+        if (params.wfn == "EOM_CC2") {
+            short_name = "CC2";
+        } else if (params.wfn == "EOM_CC3") {
+            short_name = "CC3";
+        } else if (params.wfn == "EOM_CCSD") {
+            short_name = "CCSD";
+        }
+
         /* Transform Da/b to so basis and set in wfn */
         // If this becomes a wavefunction subclass someday, just redefine the densities directly.
         if (ref_wfn->same_a_b_dens()) {
@@ -385,8 +395,7 @@ PsiReturnType ccdensity(std::shared_ptr<ccenergy::CCEnergyWavefunction> ref_wfn,
                 auto ref_Da_so = ref_wfn->Da();
                 ref_Da_so->copy(Pa_so);
             } else {
-                auto var_title = "CC ROOT " + std::to_string(i) + " Da";
-                ref_wfn->set_array_variable(var_title, Pa_so);
+                density_saver(*ref_wfn, &(rho_params[i]), "Da", Pa_so);
             }
         } else {
             auto Pa_so = linalg::triplet(ref_wfn->Ca(), Pa_x, ref_wfn->Ca(), false, false, true);
@@ -397,17 +406,19 @@ PsiReturnType ccdensity(std::shared_ptr<ccenergy::CCEnergyWavefunction> ref_wfn,
                 ref_Da_so->copy(Pa_so);
                 ref_Db_so->copy(Pb_so);
             } else {
-                auto var_title = "CC ROOT " + std::to_string(i) + " Da";
-                ref_wfn->set_array_variable(var_title, Pa_so);
-                var_title = "CC ROOT " + std::to_string(i) + " Db";
-                ref_wfn->set_array_variable(var_title, Pb_so);
+                density_saver(*ref_wfn, &(rho_params[i]), "Da", Pa_so);
+                density_saver(*ref_wfn, &(rho_params[i]), "Db", Pb_so);
             }
         }
 
         // For psivar scraper
 
-        // Process::environment.globals["CC ROOT n DIPOLE"]
-        // Process::environment.globals["CC ROOT n QUADRUPOLE"]
+        // Process::environment.globals["CCname ROOT n DIPOLE"]
+        // Process::environment.globals["CCname ROOT n DIPOLE - h TRANSITION"]
+        // Process::environment.globals["CCname ROOT n (h) DIPOLE"]
+        // Process::environment.globals["CCname ROOT n QUADRUPOLE"]
+        // Process::environment.globals["CCname ROOT n QUADRUPOLE - h TRANSITION"]
+        // Process::environment.globals["CCname ROOT n (h) QUADRUPOLE"]
 
         free_block(moinfo.opdm);
         free_block(moinfo.opdm_a);

psi4/src/psi4/cc/ccdensity/ccdensity.h

@@ -45,6 +45,8 @@ namespace ccdensity {
 void scalar_saver_ground(ccenergy::CCEnergyWavefunction& wfn, struct TD_Params *S, const std::string suffix, double val);
 // Save a scalar describing a excited->excited transition.
 void scalar_saver_excited(ccenergy::CCEnergyWavefunction& wfn, struct TD_Params *S, struct TD_Params *U, const std::string suffix, double val);
+// Save a density.
+void density_saver(ccenergy::CCEnergyWavefunction& wfn, struct RHO_Params *S, const std::string suffix, SharedMatrix val);
 
 }
 }

psi4/src/psi4/cc/ccdensity/ex_oscillator_strength.cc

@@ -283,8 +283,14 @@ void ex_oscillator_strength(ccenergy::CCEnergyWavefunction& wfn, struct TD_Param
     outfile->Printf("\tEinstein B Coefficient   %11.8e \n", einstein_b);
 
     // Save variables to wfn.
+    // Process::environment.globals["CCname ROOT n -> ROOT m OSCILLATOR STRENGTH (LEN)"]
+    // Process::environment.globals["CCname ROOT n -> ROOT m OSCILLATOR STRENGTH (LEN) - h TRANSITION"]
     // Process::environment.globals["CCname ROOT n (h) -> ROOT m (i) OSCILLATOR STRENGTH (LEN)"]
+    // Process::environment.globals["CCname ROOT n -> ROOT m EINSTEIN A (LEN)"]
+    // Process::environment.globals["CCname ROOT n -> ROOT m EINSTEIN A (LEN) - h TRANSITION"]
     // Process::environment.globals["CCname ROOT n (h) -> ROOT m (i) EINSTEIN A (LEN)"]
+    // Process::environment.globals["CCname ROOT n -> ROOT m EINSTEIN B (LEN)"]
+    // Process::environment.globals["CCname ROOT n -> ROOT m EINSTEIN B (LEN) - h TRANSITION"]
     // Process::environment.globals["CCname ROOT n (h) -> ROOT m (i) EINSTEIN B (LEN)"]
     scalar_saver_excited(wfn, S, U, "OSCILLATOR STRENGTH (LEN)", f_x + f_y + f_z);
     scalar_saver_excited(wfn, S, U, "EINSTEIN A (LEN)", einstein_a);

psi4/src/psi4/cc/ccdensity/ex_rotational_strength.cc

@@ -153,6 +153,8 @@ void ex_rotational_strength(ccenergy::CCEnergyWavefunction& wfn, struct TD_Param
     outfile->Printf("\tRotational Strength (10^-40 esu^2 cm^2)  %11.8lf\n", rs * _au2cgs);
 
     // Save rotary strength to wfn.
+    // Process::environment.globals["CCname ROOT n -> ROOT m ROTARY STRENGTH (LEN)"]
+    // Process::environment.globals["CCname ROOT n -> ROOT m ROTARY STRENGTH (LEN) - h TRANSITION"]
     // Process::environment.globals["CCname ROOT n (h) -> ROOT m (i) ROTARY STRENGTH (LEN)"]
     scalar_saver_excited(wfn, S, U, "ROTARY STRENGTH (LEN)", rs);
 
@@ -253,6 +255,8 @@ void ex_rotational_strength(ccenergy::CCEnergyWavefunction& wfn, struct TD_Param
     outfile->Printf("\tRotational Strength (10^-40 esu^2 cm^2)  %11.8lf\n", rs * _au2cgs);
 
     // Save rotary strength to wfn.
+    // Process::environment.globals["CCname ROOT n -> ROOT m ROTARY STRENGTH (VEL)"]
+    // Process::environment.globals["CCname ROOT n -> ROOT m ROTARY STRENGTH (VEL) - h TRANSITION"]
     // Process::environment.globals["CCname ROOT n (h) -> ROOT m (i) ROTARY STRENGTH (VEL)"]
     scalar_saver_excited(wfn, S, U, "ROTARY STRENGTH (VEL)", rs);
 

psi4/src/psi4/cc/ccdensity/wfn_saver.cc

@@ -44,7 +44,7 @@ namespace ccdensity {
 void scalar_saver_ground(ccenergy::CCEnergyWavefunction& wfn, struct TD_Params *S, const std::string suffix, double val) {
     auto target_sym = moinfo.sym ^ S->irrep;
     auto idx_num = S->root + static_cast<int>(S->irrep == 0);
-    auto total_idx = wfn.state_idx_to_identifiers[{idx_num, target_sym}];
+    auto total_idx = wfn.total_indices[{idx_num, target_sym}];
     auto trans_irr_lbl = moinfo.labels[moinfo.sym ^ target_sym];
     std::unordered_set<std::string> names {"CC"};
     if (params.wfn == "EOM_CCSD") {
@@ -69,8 +69,8 @@ void scalar_saver_excited(ccenergy::CCEnergyWavefunction& wfn, struct TD_Params
     auto U_sym = moinfo.sym ^ U->irrep;
     auto S_idx = S->root + static_cast<int>(S->irrep == 0);
     auto U_idx = U->root + static_cast<int>(U->irrep == 0);
-    auto S_total_idx = wfn.state_idx_to_identifiers[{S_idx, S_sym}];
-    auto U_total_idx = wfn.state_idx_to_identifiers[{U_idx, U_sym}];
+    auto S_total_idx = wfn.total_indices[{S_idx, S_sym}];
+    auto U_total_idx = wfn.total_indices[{U_idx, U_sym}];
     auto trans_irr_lbl = moinfo.labels[S_sym ^ U_sym];
     std::unordered_set<std::string> names {"CC"};
     if (params.wfn == "EOM_CCSD") {
@@ -83,13 +83,41 @@ void scalar_saver_excited(ccenergy::CCEnergyWavefunction& wfn, struct TD_Params
     for (const auto name: names) {
         auto varname = name + " ROOT " + std::to_string(S_idx) + " (" + moinfo.labels[S_sym] + ") -> ROOT " + std::to_string(U_idx) + " (" + moinfo.labels[U_sym] + ") " + suffix;
         wfn.set_scalar_variable(varname, val);
-        wfn.set_scalar_variable(varname, val);
         varname = name + " ROOT " + std::to_string(S_total_idx) + " -> ROOT " + std::to_string(U_total_idx) + " " + suffix;
         wfn.set_scalar_variable(varname, val);
         varname = name + " ROOT " + std::to_string(S_total_idx) + " -> ROOT " + std::to_string(U_total_idx) + " " + suffix + " - " + trans_irr_lbl + " TRANSITION";
         wfn.set_scalar_variable(varname, val);
     }
 }
 
+void density_saver(ccenergy::CCEnergyWavefunction& wfn, struct RHO_Params *S, const std::string suffix, SharedMatrix val) {
+    auto target_sym = moinfo.sym ^ S->R_irr;
+    auto idx_num = S->R_root + static_cast<int>(S->R_irr == 0);
+    auto total_idx = wfn.total_indices[{idx_num, target_sym}];
+    auto trans_irr_lbl = moinfo.labels[moinfo.sym ^ target_sym];
+    std::unordered_set<std::string> names {"CC"};
+    if (params.wfn == "EOM_CCSD") {
+        names.insert("CCSD");
+    } else if (params.wfn == "EOM_CC2") {
+        names.insert("CC2");
+    } else {
+        throw PSIEXCEPTION("Unknown wfn type");
+    }
+
+    if (suffix != "Da" and suffix != "Db") {
+        throw PSIEXCEPTION("Unknown spin type.");
+    }
+    std::map<std::string, SharedMatrix>& density_map = (suffix == "Da") ? wfn.Da_map_ : wfn.Db_map_;
+
+    for (const auto name: names) {
+        auto varname = name + " ROOT " + std::to_string(idx_num) + " (" + moinfo.labels[target_sym] + ") " + suffix;
+        density_map[varname] = val;
+        varname = name + " ROOT " + std::to_string(total_idx) + " " + suffix;
+        density_map[varname] = val;
+        varname = name + " ROOT " + std::to_string(total_idx) + " " + suffix + " - " + trans_irr_lbl + " TRANSITION";
+        density_map[varname] = val;
+    }
+}
+
 }
 }