diff --git a/src/Hamiltonian.cc b/src/Hamiltonian.cc index 5d3683ec..c7efa2f6 100644 --- a/src/Hamiltonian.cc +++ b/src/Hamiltonian.cc @@ -141,7 +141,7 @@ void Hamiltonian::applyLocal(const int ncolors, T& phi, T& hphi) for (int i = 0; i < ncolors; i++) { using memory_space_type = typename T::memory_space_type; - ORBDTYPE* ihphi = hphi.getPsi(i); + auto ihphi = hphi.getPsi(i); unsigned int const size = hphi.getNumpt(); ORBDTYPE* ihphi_host_view = MemorySpace::Memory::allocate_host_view(size); @@ -221,6 +221,8 @@ template <> void Hamiltonian::addHlocal2matrix(LocGridOrbitals& phi1, LocGridOrbitals& phi2, ReplicatedMatrix& hij, const bool force) { + (void)hij; + applyLocal(phi2, force); // phi1.addDotWithNcol2Matrix(*hlphi_, hij); diff --git a/src/Hamiltonian.h b/src/Hamiltonian.h index 08194f87..0c87ceda 100644 --- a/src/Hamiltonian.h +++ b/src/Hamiltonian.h @@ -44,7 +44,7 @@ class Hamiltonian template void addHlocal2matrix(OrbitalsType& orbitals1, OrbitalsType& orbitals2, - MatrixType& mat, const bool force = false); + MatrixType& mat, const bool force); void addHlocalij(OrbitalsType& orbitals1, OrbitalsType& orbitals2, ProjectedMatricesInterface*); void addHlocalij(OrbitalsType& orbitals1, ProjectedMatricesInterface*); diff --git a/src/HamiltonianMVPSolver.cc b/src/HamiltonianMVPSolver.cc index 715f52ba..17f1f66b 100644 --- a/src/HamiltonianMVPSolver.cc +++ b/src/HamiltonianMVPSolver.cc @@ -151,7 +151,7 @@ int HamiltonianMVPSolver::solve( // compute new h11 for the current potential by adding local part to // nonlocal components h11 = h11nl; - hamiltonian_->addHlocal2matrix(orbitals, orbitals, h11); + hamiltonian_->addHlocal2matrix(orbitals, orbitals, h11, false); projmatrices->assignH(h11); projmatrices->setHB2H(); @@ -179,7 +179,7 @@ int HamiltonianMVPSolver::solve( // update H and compute energy at midpoint h11 = h11nl; - hamiltonian_->addHlocal2matrix(orbitals, orbitals, h11); + hamiltonian_->addHlocal2matrix(orbitals, orbitals, h11, false); projmatrices->assignH(h11); projmatrices->setHB2H(); @@ -214,7 +214,7 @@ int HamiltonianMVPSolver::solve( // update H with new potential h11 = h11nl; - hamiltonian_->addHlocal2matrix(orbitals, orbitals, h11); + hamiltonian_->addHlocal2matrix(orbitals, orbitals, h11, false); projmatrices->assignH(h11); projmatrices->setHB2H(); @@ -270,7 +270,7 @@ int HamiltonianMVPSolver::solve( // update H h11 = h11nl; - hamiltonian_->addHlocal2matrix(orbitals, orbitals, h11); + hamiltonian_->addHlocal2matrix(orbitals, orbitals, h11, false); projmatrices->assignH(h11); projmatrices->setHB2H(); diff --git a/src/MGmol.cc b/src/MGmol.cc index c20fef32..2c61447d 100644 --- a/src/MGmol.cc +++ b/src/MGmol.cc @@ -1171,33 +1171,28 @@ void MGmol::projectOutKernel(OrbitalsType& phi) } template -void MGmol::setGamma( - const pb::Lap& lapOper, const Potentials& pot) +void MGmol::precond_mg(OrbitalsType& phi) { assert(orbitals_precond_); Control& ct = *(Control::instance()); - orbitals_precond_->setGamma( - lapOper, pot, ct.getMGlevels(), proj_matrices_.get()); -} + Potentials& pot = hamiltonian_->potential(); + pb::Lap* lapOper = hamiltonian_->lapOper(); -template -void MGmol::precond_mg(OrbitalsType& phi) -{ - assert(orbitals_precond_); + orbitals_precond_->setGamma( + *lapOper, pot, ct.getMGlevels(), proj_matrices_.get()); orbitals_precond_->precond_mg(phi); } template -double MGmol::computeResidual(OrbitalsType& orbitals, - OrbitalsType& work_orbitals, Ions& ions, OrbitalsType& res, - const bool print_residual, const bool norm_res) +double MGmol::computeResidual(OrbitalsType& phi, + OrbitalsType& hphi, Ions& ions, OrbitalsType& res, + const KBPsiMatrixSparse* const kbpsi, const bool print_residual, + const bool norm_res) { - assert(orbitals.getIterativeIndex() >= 0); - comp_res_tm_.start(); // os_<<"computeResidual()"<::computeResidual(OrbitalsType& orbitals, proj_matrices_->computeInvB(); - Potentials& pot = hamiltonian_->potential(); - pb::Lap* lapop = hamiltonian_->lapOper(); - - setGamma(*lapop, pot); - - // get H*psi stored in work_orbitals.psi + // get H*phi stored in hphi // and psi^T H psi in Hij - getHpsiAndTheta(ions, orbitals, work_orbitals); + getHpsiAndTheta(ions, phi, hphi, kbpsi); - double norm2Res = computeConstraintResidual( - orbitals, work_orbitals, res, print_residual, norm_res); + double norm2Res + = computeConstraintResidual(phi, hphi, res, print_residual, norm_res); - if (ct.isSpreadFunctionalEnergy()) addResidualSpreadPenalty(orbitals, res); + if (ct.isSpreadFunctionalEnergy()) addResidualSpreadPenalty(phi, res); comp_res_tm_.stop(); @@ -1344,19 +1334,8 @@ double MGmol::computePrecondResidual(OrbitalsType& phi, { Control& ct = *(Control::instance()); - proj_matrices_->computeInvB(); - - Potentials& pot = hamiltonian_->potential(); - pb::Lap* lapop = hamiltonian_->lapOper(); - - setGamma(*lapop, pot); - - // get H*psi stored in hphi - // and psi^T H psi in Hij - getHpsiAndTheta(ions, phi, hphi, kbpsi); - - double norm2Res - = computeConstraintResidual(phi, hphi, res, print_residual, norm_res); + double norm2Res = computeResidual( + phi, hphi, ions, res, kbpsi, print_residual, norm_res); if (ct.withPreconditioner()) { @@ -1366,8 +1345,6 @@ double MGmol::computePrecondResidual(OrbitalsType& phi, orbitals_precond_->precond_mg(res); } - // if( ct.isSpreadFunctionalActive() )addResidualSpreadPenalty(phi,res); - return norm2Res; } diff --git a/src/MGmol.h b/src/MGmol.h index 9c400b54..9cdc345d 100644 --- a/src/MGmol.h +++ b/src/MGmol.h @@ -321,10 +321,17 @@ class MGmol : public MGmolInterface void projectOutKernel(OrbitalsType& phi); void precond_mg(OrbitalsType& orbitals); - void setGamma(const pb::Lap& lapOper, const Potentials& pot); + double computeResidual(OrbitalsType& orbitals, OrbitalsType& work_orbitals, + Ions& ions, OrbitalsType& res, const KBPsiMatrixSparse* const kbpsi, + const bool print_residual, const bool norm_res); double computeResidual(OrbitalsType& orbitals, OrbitalsType& work_orbitals, Ions& ions, OrbitalsType& res, const bool print_residual, - const bool norm_res); + const bool norm_res) + { + return computeResidual(orbitals, work_orbitals, ions, res, + g_kbpsi_.get(), print_residual, norm_res); + } + void applyAOMMprojection(OrbitalsType&); void force(OrbitalsType& orbitals, Ions& ions) { diff --git a/src/MVPSolver.cc b/src/MVPSolver.cc index 29021a95..787c9bdb 100644 --- a/src/MVPSolver.cc +++ b/src/MVPSolver.cc @@ -239,7 +239,7 @@ int MVPSolver::solve(OrbitalsType& orbitals) // compute h11 for the current potential by adding local part to // nonlocal components MatrixType h11(h11_nl); - hamiltonian_->addHlocal2matrix(orbitals, orbitals, h11); + hamiltonian_->addHlocal2matrix(orbitals, orbitals, h11, false); current_proj_mat->assignH(h11); current_proj_mat->setHB2H(); @@ -317,10 +317,9 @@ int MVPSolver::solve(OrbitalsType& orbitals) energy_->saveVofRho(); // update h11 - { - h11 = h11_nl; - hamiltonian_->addHlocal2matrix(orbitals, orbitals, h11); - } + h11 = h11_nl; + hamiltonian_->addHlocal2matrix( + orbitals, orbitals, h11, false); proj_mat_work_->assignH(h11); proj_mat_work_->setHB2H(); diff --git a/src/computeHij.cc b/src/computeHij.cc index d229c89d..2786d4c6 100644 --- a/src/computeHij.cc +++ b/src/computeHij.cc @@ -197,7 +197,7 @@ void MGmol::computeHij_private(OrbitalsType& orbitals_i, ss2dm->accumulate(submat, hij, 0.); // add local Hamiltonian part to phi^T*H*phi - hamiltonian_->addHlocal2matrix(orbitals_i, orbitals_j, hij); + hamiltonian_->addHlocal2matrix(orbitals_i, orbitals_j, hij, false); } template @@ -324,13 +324,6 @@ void MGmol::computeHnlPhiAndAdd2HPhi(Ions& ions, hphi.setIterativeIndex(phi.getIterativeIndex()); } -template -void MGmol::getHpsiAndTheta( - Ions& ions, OrbitalsType& phi, OrbitalsType& hphi) -{ - getHpsiAndTheta(ions, phi, hphi, g_kbpsi_.get()); -} - template void MGmol::getHpsiAndTheta(Ions& ions, OrbitalsType& phi, OrbitalsType& hphi, const KBPsiMatrixSparse* const kbpsi) @@ -345,7 +338,7 @@ void MGmol::getHpsiAndTheta(Ions& ions, OrbitalsType& phi, os_ << " getHpsiAndTheta" << std::endl; #endif - hphi.assign(hamiltonian_->applyLocal(phi)); + hamiltonian_->applyLocal(phi.chromatic_number(), phi, hphi); // Compute "nstates" columns of matrix // Hij = phi**T * H_loc * phi and save in sh @@ -370,7 +363,13 @@ void MGmol::getHpsiAndTheta(Ions& ions, OrbitalsType& phi, kbpsi->computeHvnlMatrix(ions, proj_matrices_.get()); // add local part of H to sh - hamiltonian_->addHlocalij(phi, proj_matrices_.get()); + SquareLocalMatrices slm( + phi.subdivx(), phi.chromatic_number()); + + phi.computeLocalProduct(hphi, slm); + proj_matrices_->setLocalMatrixElementsHl(slm); + + proj_matrices_->consolidateH(); energy_->saveVofRho();