From d304e2e69cf03d54cc4640b0e7a84381af2ac1db Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niklas=20G=C3=B6llmann?= Date: Tue, 30 Apr 2024 15:39:40 +0200 Subject: [PATCH] small bugfix in ExchangePotentialD2.cpp --- src/qmoperators/two_electron/ExchangePotentialD2.cpp | 8 ++++---- src/scf_solver/GroundStateSolver.cpp | 5 +++-- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/src/qmoperators/two_electron/ExchangePotentialD2.cpp b/src/qmoperators/two_electron/ExchangePotentialD2.cpp index e263b7fe3..12798d09c 100644 --- 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 index 9467aa707..3579eb3a3 100644 --- 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 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> sortme; for (unsigned int q = 0; q < p.size(); ++q) {