checkpoint: He lrcc2 works

m-a-d-n-e-s-s · Jul 19, 2024 · 41393a8 · 41393a8
1 parent 58b69bd
commit 41393a8
Show file tree

Hide file tree

Showing 4 changed files with 14 additions and 11 deletions.
diff --git a/src/madness/chem/CC2.cc b/src/madness/chem/CC2.cc
@@ -859,9 +859,9 @@ CC2::solve_lrcc2(Pairs<CCPair>& gs_doubles, const CC_vecfunction& gs_singles, co
     const double omega_cis = ex_singles.omega;
 
     for (size_t iter = 0; iter < parameters.iter_max(); iter++) {
-        print_header2("Macroiteration " + std::to_string(int(iter)) + " of LRCC2 for energy "+std::to_string(ex_singles.omega));
-        bool sconv = iterate_lrcc2_pairs(world, gs_singles, ex_singles, ex_doubles, info);
-        bool dconv = iterate_lrcc2_singles(world, gs_singles, gs_doubles, ex_singles, ex_doubles, info);
+        print_header2("Macroiteration " + std::to_string(int(iter)) + " of LRCC2 for excitation energy "+std::to_string(ex_singles.omega));
+        bool dconv = iterate_lrcc2_pairs(world, gs_singles, ex_singles, ex_doubles, info);
+        bool sconv = iterate_lrcc2_singles(world, gs_singles, gs_doubles, ex_singles, ex_doubles, info);
         // update_reg_residues_ex(world, gs_singles, ex_singles, ex_doubles, info);
         if (sconv and dconv) break;
     }
@@ -1090,7 +1090,7 @@ CC2::initialize_pairs(Pairs<CCPair>& pairs, const CCState ftype, const CalcType
 //                tmp.excitation = excitation;
                 tmp.constant_part = const_part;
                 // tmp.constant_part = vconst_part[0];
-                print_header1("loading constant part");
+                // print_header1("loading constant part");
                 pairs.insert(i, j, tmp);
                 CCPotentials::compute_excited_pair_energy(world, pairs(i, j), x, info);
             } else error("Unknown pairtype");

diff --git a/src/madness/chem/CC2.h b/src/madness/chem/CC2.h
@@ -307,9 +307,10 @@ class CC2 : public OptimizationTargetInterface, public QCPropertyInterface {
 
             // Normalize Singles if it is excited state
             if (ctype == CT_LRCCS or ctype == CT_LRCC2 or ctype == CT_ADC2) {
-                output("Normalizing new singles");
-                const double norm=inner(GV,info.R_square*GV);
-                scale(world, GV, 1.0 / norm);
+                Nemo::normalize(GV, info.R);
+                // output("Normalizing new singles");
+                // const double norm=inner(GV,info.R_square*GV);
+                // scale(world, GV, 1.0 / norm);
             } else output("Singles not normalized");
 
             // residual
@@ -353,9 +354,10 @@ class CC2 : public OptimizationTargetInterface, public QCPropertyInterface {
             singles.omega = omega;
             vector_real_function_3d new_singles = truncate(GV);
             if (info.parameters.kain()) new_singles = solver.update(singles.get_vecfunction(), residual);
-            if (info.parameters.debug()) {
-                print_size(world, new_singles, "new_singles");
-            }
+            if (info.parameters.debug()) print_size(world, new_singles, "new_singles");
+            if (ctype == CT_LRCCS or ctype == CT_LRCC2 or ctype == CT_ADC2) Nemo::normalize(new_singles, info.R);
+            if (info.parameters.debug()) print_size(world, new_singles, "new_singles normalized");
+
             for (size_t i = 0; i < GV.size(); i++) {
                 singles(i + info.parameters.freeze()).function = copy(new_singles[i]);
             }

diff --git a/src/madness/chem/CCPotentials.h b/src/madness/chem/CCPotentials.h
@@ -35,6 +35,7 @@ class CCPotentials {
         info.molecular_coordinates = nemo->get_calc()->molecule.get_all_coords_vec();
         info.parameters = parameters;
         info.R_square = nemo->R_square;
+        info.R = nemo->R;
         info.U1 = nemo->ncf->U1vec();
         info.U2 = nemo->ncf->U2();
         info.intermediate_potentials = get_potentials;

diff --git a/src/madness/chem/CCStructures.h b/src/madness/chem/CCStructures.h
@@ -1209,7 +1209,7 @@ struct Info {
     CCParameters parameters;
     std::vector<double> orbital_energies;
     CCIntermediatePotentials intermediate_potentials;
-    Function<double,3> R_square, U2;
+    Function<double,3> R_square, U2, R;;
     std::vector<Function<double,3>> U1;
 
     vector_real_function_3d get_active_mo_ket() const {