Only ever compute gamma, keep update schemes for future reference

src/scf_solver/ExcitedStatesSolver.cpp

@@ -81,6 +81,9 @@ json ExcitedStatesSolver::optimize(Molecule &mol, FockBuilder &F_0, FockBuilder
 
     RankZeroOperator V_0 = F_0.potential();
     RankZeroOperator V_1 = F_1.potential();
+
+    bool use_harrison = false;  // use harrison´s update scheme for the excitation energies
+    bool update_omega = false;  // use Kottmann´s update scheme for the excitation energies
 
     double err_o = 1.0;
     double err_t = 1.0;
@@ -154,10 +157,11 @@ json ExcitedStatesSolver::optimize(Molecule &mol, FockBuilder &F_0, FockBuilder
 
             // Apply Helmholtz operators
             OrbitalVector X_np1 = H_x.apply(V_0, X_n, Psi);
-            Psi.clear();
+
             // Projecting (1 - rho_0)X
             mrcpp::print::header(2, "Projecting occupied space");
             t_lap.start();
+
             orbital::orthogonalize(this->orth_prec, X_np1, Phi_0);
             orbital::orthogonalize(this->orth_prec, X_np1);
 
@@ -168,23 +172,40 @@ json ExcitedStatesSolver::optimize(Molecule &mol, FockBuilder &F_0, FockBuilder
             // orbital::orthogonalize(this->orth_prec, X_np1); // X_n orbitals should be orthogonal wrt. each other
             // before next iteration, maybe do this before kain acceleration
             // Compute update and errors
+
+            if (update_omega) {
+                domega_n = updateOmega(X_n, X_np1);
+            }
+
             OrbitalVector dX_n = orbital::add(1.0, X_np1, -1.0, X_n);
-            domega_n = updateOmega(X_n, X_np1);
             errors_x = orbital::get_norms(dX_n);
 
+
+
             // Compute KAIN update:
             kain_x.accelerate(orb_prec, X_n, dX_n);
+
+
+            if (use_harrison) {
+                auto V_0_x = V_0(X_n);
+                auto left_hand = orbital::add(1.0, V_0_x, 1.0, Psi);
+                X_np1 = orbital::add(1.0, X_n, 1.0, dX_n);
+                domega_n = - orbital::dot(left_hand, dX_n).sum().real()/orbital::dot(X_np1, X_np1).sum().real();
+            }
+
+            Psi.clear();
 
             // Prepare for next iteration
             X_n = orbital::add(1.0, X_n, 1.0, dX_n);
             // orbital::orthogonalize(this->orth_prec, X_n);
 
             // Setup perturbed Fock operator (including V_1)
-            F_1.setup(orb_prec); // do the x orbitals being uodated change this? it obviously should, but check
+            V_1.setup(orb_prec); 
 
             // Compute omega
             mrcpp::print::header(2, "Computing frequency update");
             t_lap.start();
+
             auto omega_np1 = computeOmega(Phi_0, X_n, F_0, V_1); /*  computeOmega(Phi_0, X_n, Y_n, F_0, F_1, F_mat_0); */
             domega_n = omega_np1 - omega_n;                      // maybe I should do this before normalization
             omega_n += domega_n;