analyze MOs in complex calculations

m-a-d-n-e-s-s · Aug 23, 2024 · 7df9447 · 7df9447
1 parent 0e709ee
commit 7df9447
Show file tree

Hide file tree

Showing 5 changed files with 30 additions and 9 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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;
 }