From 7df9447c5667bdf14b9979e1b63d9966a946f3d7 Mon Sep 17 00:00:00 2001
From: fbischoff <florian.bischoff@hu-berlin.de>
Date: Fri, 23 Aug 2024 12:19:00 -0400
Subject: [PATCH] analyze MOs in complex calculations

---
 src/madness/chem/SCF.cc           | 14 +++++++++-----
 src/madness/chem/SCF.h            |  7 +++++--
 src/madness/chem/molecularbasis.h |  2 +-
 src/madness/chem/oep.cc           |  5 ++++-
 src/madness/chem/znemo.cc         | 11 +++++++++++
 5 files changed, 30 insertions(+), 9 deletions(-)

diff --git a/src/madness/chem/SCF.cc b/src/madness/chem/SCF.cc
index 7572602fa06..01be3b90f46 100644
--- a/src/madness/chem/SCF.cc
+++ b/src/madness/chem/SCF.cc
@@ -592,8 +592,11 @@ vecfuncT SCF::project_ao_basis_only(World& world, const AtomicBasisSet& aobasis,
     return ao;
 }
 
-void SCF::analyze_vectors(World& world, const vecfuncT& mo, const tensorT& occ,
-                          const tensorT& energy, const std::vector<int>& set) {
+void SCF::analyze_vectors(World& world, const vecfuncT& mo,
+            const vecfuncT& ao, double vtol,
+            const Molecule& molecule, const int print_level,
+            const AtomicBasisSet& aobasis, const tensorT& occ,
+            const tensorT& energy, const std::vector<int>& set) {
     START_TIMER(world);
     PROFILE_MEMBER_FUNC(SCF);
     tensorT Saomo = matrix_inner(world, ao, mo);
@@ -624,7 +627,8 @@ void SCF::analyze_vectors(World& world, const vecfuncT& mo, const tensorT& occ,
     for (long i = 0; i < nmo; ++i) {
         size_t ncoeffi = mo[i].size();
         ncoeff += ncoeffi;
-        if (world.rank() == 0 and (param.print_level() > 1)) {
+        // if (world.rank() == 0 and (param.print_level() > 1)) {
+        if (world.rank() == 0 and (print_level > 1)) {
             printf("  MO%4ld : ", i);
             if (set.size())
                 printf("set=%d : ", set[i]);
@@ -2355,12 +2359,12 @@ void SCF::solve(World& world) {
     }
 
     if (param.nwfile() == "none") {
-        analyze_vectors(world, amo, aocc, aeps);
+        analyze_vectors(world, amo, ao, vtol, molecule, param.print_level(), aobasis, aocc, aeps);
         if (param.nbeta() != 0 && !param.spin_restricted()) {
             if (world.rank() == 0 and (param.print_level() > 1))
                 print("Analysis of beta MO vectors");
 
-            analyze_vectors(world, bmo, bocc, beps);
+            analyze_vectors(world, bmo, ao, vtol, molecule, param.print_level(), aobasis, bocc, beps);
         }
     }
 
diff --git a/src/madness/chem/SCF.h b/src/madness/chem/SCF.h
index 2e284fdd593..b99fc21b4eb 100644
--- a/src/madness/chem/SCF.h
+++ b/src/madness/chem/SCF.h
@@ -352,8 +352,11 @@ class SCF {
     std::vector<int> group_orbital_sets(World& world, const tensorT& eps,
                                         const tensorT& occ, const int nmo) const;
 
-    void analyze_vectors(World& world, const vecfuncT& mo, const tensorT& occ = tensorT(),
-                         const tensorT& energy = tensorT(), const std::vector<int>& set = std::vector<int>());
+    static void analyze_vectors(World& world, const vecfuncT& mo,
+            const vecfuncT& ao, double vtol,
+            const Molecule& molecule, const int print_level,
+            const AtomicBasisSet& aobasis, const tensorT& occ = tensorT(),
+            const tensorT& energy = tensorT(), const std::vector<int>& set = std::vector<int>());
 
     distmatT kinetic_energy_matrix(World& world, const vecfuncT& v) const;
 
diff --git a/src/madness/chem/molecularbasis.h b/src/madness/chem/molecularbasis.h
index d201e519124..51205321656 100644
--- a/src/madness/chem/molecularbasis.h
+++ b/src/madness/chem/molecularbasis.h
@@ -683,7 +683,7 @@ class AtomicBasisSet {
     /// - basis function number
     /// - MO coeff
     template <typename T>
-    void print_anal(const Molecule& molecule, const Tensor<T>& v) {
+    void print_anal(const Molecule& molecule, const Tensor<T>& v) const {
         const double thresh = 0.2*v.normf();
         if (thresh == 0.0) {
             printf("    zero vector\n");
diff --git a/src/madness/chem/oep.cc b/src/madness/chem/oep.cc
index f9adb8d0b38..f9509963654 100644
--- a/src/madness/chem/oep.cc
+++ b/src/madness/chem/oep.cc
@@ -244,7 +244,10 @@ double OEP::iterate(const std::string model, const vecfuncT& HF_nemo, const tens
 	    if (param.do_localize()) {
 	    	for (size_t i=0; i<KS_nemo.size(); ++i) calc->aeps(i)=KS_Fock(i,i);
 	    	KS_nemo=localize(KS_nemo,param.econv(),iter==0);
-	    	if (param.print_level()>=10) calc->analyze_vectors(world,KS_nemo,calc->aocc,tensorT(),calc->aset);
+	    	if (param.print_level()>=10)
+	    		SCF::analyze_vectors(world, KS_nemo, calc->ao, calc->vtol, calc->molecule, param.print_level(),
+	    			calc->aobasis, calc->aocc, tensorT(), calc->aset );
+	    	// calc->analyze_vectors(world,KS_nemo,calc->aocc,tensorT(),calc->aset);
 	    }
 	    if (do_symmetry()) {
 		    std::vector<std::string> str_irreps;
diff --git a/src/madness/chem/znemo.cc b/src/madness/chem/znemo.cc
index 3c68b8c8a8a..701d31a3455 100644
--- a/src/madness/chem/znemo.cc
+++ b/src/madness/chem/znemo.cc
@@ -426,6 +426,17 @@ double Znemo::analyze() {
 	save_orbitals("plot");
 	save(density,"density");
 	save(spindensity,"spindensity");
+
+	auto components=std::vector<std::vector<real_function_3d>>({real(amo),imag(amo),real(bmo),imag(bmo)});
+	auto component_names=std::vector<std::string>({"real_amo","imag_amo","real_bmo","imag_bmo"});
+	std::vector<real_function_3d> real_aos=SCF::project_ao_basis_only(world, aobasis, mol);
+	for (size_t i=0; i<components.size(); ++i) {
+		if (world.rank()==0) print("analysis of MO component ",component_names[i]);
+		if (components[i].size()>0) SCF::analyze_vectors(world, components[i], real_aos,
+			FunctionDefaults<3>::get_thresh()*0.1, molecule(), cparam.print_level(), aobasis);
+
+	}
+
     return energy;
 }