small bugfix in ExchangePotentialD2.cpp

MRChemSoft · May 2, 2024 · 5c4b85f · 5c4b85f
1 parent ca71fce
commit 5c4b85f
Show file tree

Hide file tree

Showing 2 changed files with 7 additions and 6 deletions.
diff --git a/src/qmoperators/two_electron/ExchangePotentialD2.cpp b/src/qmoperators/two_electron/ExchangePotentialD2.cpp
@@ -130,8 +130,8 @@ Orbital ExchangePotentialD2::apply(Orbital phi_p) {
             calcExchange_kij(precf, phi_i, y_i, phi_p, ex_iyp);
             func_vec.push_back(ex_xip);
             func_vec.push_back(ex_iyp);
-            coef_vec.push_back(spin_fac / phi_i.squaredNorm());
-            coef_vec.push_back(spin_fac / phi_i.squaredNorm());
+            coef_vec.push_back(spin_fac * phi_i.occ() / phi_i.squaredNorm());
+            coef_vec.push_back(spin_fac * phi_i.occ() / phi_i.squaredNorm());
         }
         if (not mrcpp::mpi::my_orb(phi_i)) phi_i.free(NUMBER::Total);
         if (not mrcpp::mpi::my_orb(x_i)) x_i.free(NUMBER::Total);
@@ -189,8 +189,8 @@ Orbital ExchangePotentialD2::dagger(Orbital phi_p) {
             calcExchange_kij(precf, y_i, phi_i, phi_p, ex_yip);
             func_vec.push_back(ex_ixp);
             func_vec.push_back(ex_yip);
-            coef_vec.push_back(spin_fac / phi_i.squaredNorm());
-            coef_vec.push_back(spin_fac / phi_i.squaredNorm());
+            coef_vec.push_back(spin_fac * phi_i.occ() / phi_i.squaredNorm());
+            coef_vec.push_back(spin_fac * phi_i.occ() / phi_i.squaredNorm());
         }
         if (not mrcpp::mpi::my_orb(phi_i)) phi_i.free(NUMBER::Total);
         if (not mrcpp::mpi::my_orb(x_i)) x_i.free(NUMBER::Total);

diff --git a/src/scf_solver/GroundStateSolver.cpp b/src/scf_solver/GroundStateSolver.cpp
@@ -472,7 +472,8 @@ bool GroundStateSolver::needDiagonalization(int nIter, bool converged) const {
  */
 DoubleVector GroundStateSolver::getNewOccupations(OrbitalVector &Phi_n, OrbitalVector &Phi_mom) {
     ComplexMatrix overlap = orbital::calc_overlap_matrix(Phi_mom, Phi_n);
-    DoubleVector occ = orbital::get_occupations(Phi_mom);
+    std::cout << "overlap: " << std::endl << overlap << std::endl;
+    DoubleVector occ = orbital::get_occupations(Phi_mom);// get occupation numbers of the orbitals of the first iteration
     // get all unique occupation numbers
     std::set<double> occupationNumbers(occ.begin(), occ.end());
     // for each unique occupation number, determine which orbitals should be assigned this occupation number
@@ -501,7 +502,7 @@ DoubleVector GroundStateSolver::getNewOccupations(OrbitalVector &Phi_n, OrbitalV
     // only consider overlap with orbitals with the current unique occupation number
     ComplexMatrix occOverlap = currOcc.asDiagonal() * overlap;
     ComplexVector p = occOverlap.colwise().norm();
-    std::cout << "p: " << p << std::endl;
+    std::cout << "p: " << std::endl << p << std::endl;
     // sort by highest overlap
     std::vector<std::pair<double, unsigned>> sortme;
     for (unsigned int q = 0; q < p.size(); ++q) {