Skip to content

Commit

Permalink
Merging PR #234 and #236 to ROMFPMD (#237)
Browse files Browse the repository at this point in the history 
* Clean up class Rho (#234)

* Clean up class MVPSolver (#236)

---------

Co-authored-by: Jean-Luc Fattebert <fattebertj@ornl.gov>
  • Loading branch information
@dreamer2368 @jeanlucf22
dreamer2368 and jeanlucf22 authored May 22, 2024
1 parent ec0cb30 commit 29c2502
Show file tree
Hide file tree
Showing 6 changed files with 67 additions and 438 deletions.
2 changes: 0 additions & 2 deletions src/DavidsonSolver.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -598,7 +598,6 @@ int DavidsonSolver<OrbitalsType, MatrixType>::solve(

rho_->computeRho(
orbitals, work_orbitals, dm11, dm12, dm21, dm22);
rho_->rescaleTotalCharge();

mgmol_strategy_->update_pot(vh_init, ions_);

Expand Down Expand Up @@ -663,7 +662,6 @@ int DavidsonSolver<OrbitalsType, MatrixType>::solve(

rho_->computeRho(
orbitals, work_orbitals, dm11, dm12, dm21, dm22);
rho_->rescaleTotalCharge();
}

} // inner iterations
Expand Down
1 change: 1 addition & 0 deletions src/MGmol.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -305,6 +305,7 @@ class MGmol : public MGmolInterface
}

OrbitalsType* loadOrbitalFromRestartFile(const std::string filename);

#ifdef MGMOL_HAS_LIBROM
int save_orbital_snapshot(std::string snapshot_dir, OrbitalsType& orbitals);
#endif
Expand Down
87 changes: 31 additions & 56 deletions src/MVPSolver.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,18 +52,22 @@ MVPSolver<OrbitalsType, MatrixType>::MVPSolver(MPI_Comm comm, std::ostream& os,
os_(os),
n_inner_steps_(n_inner_steps),
use_old_dm_(use_old_dm),
ions_(ions)
ions_(ions),
numst_(numst)
{
history_length_ = 2;
eks_history_.resize(history_length_, 100000.);
Control& ct = *(Control::instance());
if (onpe0 && ct.verbose > 0)
{
os_ << "MVPSolver..." << std::endl;
if (use_old_dm_) os_ << "MVPSolver uses old DM..." << std::endl;
}

rho_ = rho;
energy_ = energy;
electrostat_ = electrostat;
mgmol_strategy_ = mgmol_strategy;

numst_ = numst;
work_ = new MatrixType("workMVP", numst_, numst_);
work_ = new MatrixType("workMVP", numst_, numst_);

proj_mat_work_ = new ProjectedMatrices<MatrixType>(numst_, false, kbT);
proj_mat_work_->setup(global_indexes);
Expand Down Expand Up @@ -123,9 +127,9 @@ void MVPSolver<OrbitalsType, MatrixType>::buildTarget_MVP(
{
target_tm_.start();

if (onpe0) std::cout << "buildTarget_MVP()..." << std::endl;

Control& ct = *(Control::instance());
if (onpe0 && ct.verbose > 1)
std::cout << "buildTarget_MVP()..." << std::endl;

proj_mat_work_->assignH(h11);

Expand All @@ -135,10 +139,8 @@ void MVPSolver<OrbitalsType, MatrixType>::buildTarget_MVP(
//
// compute target density matrix, with occupations>0 in numst only
//
// if( onpe0 )os_<<"Compute XN..."<<endl;
proj_mat_work_->setHB2H();

// if( onpe0 )os_<<"Compute DM..."<<endl;
proj_mat_work_->updateDM(orbitals_index);

target = proj_mat_work_->dm();
Expand Down Expand Up @@ -168,18 +170,12 @@ int MVPSolver<OrbitalsType, MatrixType>::solve(OrbitalsType& orbitals)
os_ << "---------------------------------------------------------------"
"-"
<< std::endl;
os_ << "Update DM functions using MVP Solver..." << std::endl;
os_ << "Update DM using MVP Solver..." << std::endl;
os_ << "---------------------------------------------------------------"
"-"
<< std::endl;
}

// step for numerical derivative

// double nel=orbitals.projMatrices()->getNel();
// if( onpe0 )
// os_<<"NEW algorithm: Number of electrons at start = "<<nel<<endl;

KBPsiMatrixSparse kbpsi(nullptr);
kbpsi.setup(ions_);

Expand Down Expand Up @@ -227,61 +223,45 @@ int MVPSolver<OrbitalsType, MatrixType>::solve(OrbitalsType& orbitals)
energy_->saveVofRho();
}

ProjectedMatricesInterface* current_proj_mat
= (inner_it == 0) ? orbitals.getProjMatrices() : proj_mat_work_;
ProjectedMatrices<MatrixType>* current_proj_mat
= (inner_it == 0)
? dynamic_cast<ProjectedMatrices<MatrixType>*>(
orbitals.getProjMatrices())
: proj_mat_work_;

double ts0;
double e0 = 0.;
const int printE = (ct.verbose > 1) ? 1 : 0;

// compute h11 for the current potential by adding local part to
// nonlocal components
MatrixType h11(h11_nl);
mgmol_strategy_->addHlocal2matrix(orbitals, orbitals, h11);

current_proj_mat->assignH(h11);
current_proj_mat->setHB2H();

if (inner_it == 0)
{
// orbitals are new, so a few things need to recomputed
ProjectedMatrices<MatrixType>* projmatrices
= dynamic_cast<ProjectedMatrices<MatrixType>*>(
orbitals.getProjMatrices());

projmatrices->assignH(h11);
projmatrices->setHB2H();

if (use_old_dm_)
{
dmInit = projmatrices->dm();
dmInit = current_proj_mat->dm();
}

ts0 = evalEntropyMVP(projmatrices, true, os_);
e0 = energy_->evaluateTotal(
ts0, projmatrices, orbitals, printE, os_);
}
else
{

dmInit = proj_mat_work_->dm();

proj_mat_work_->assignH(h11);
proj_mat_work_->setHB2H();

ts0 = evalEntropyMVP(proj_mat_work_, true, os_);
e0 = energy_->evaluateTotal(
ts0, proj_mat_work_, orbitals, printE, os_);
}

// N x N target...
// proj_mat_work_->setHiterativeIndex(orbitals.getIterativeIndex(),
// pot.getIterativeIndex());
const double ts0 = evalEntropyMVP(current_proj_mat, true, os_);
const double e0 = energy_->evaluateTotal(
ts0, current_proj_mat, orbitals, printE, os_);

MatrixType target("target", numst_, numst_);
MatrixType delta_dm("delta_dm", numst_, numst_);

buildTarget_MVP(h11, s11, target);

if (use_old_dm_ || inner_it > 0)
{
MatrixType delta_dm("delta_dm", numst_, numst_);
delta_dm = target;
delta_dm -= dmInit;

Expand All @@ -291,18 +271,16 @@ int MVPSolver<OrbitalsType, MatrixType>::solve(OrbitalsType& orbitals)
// evaluate free energy at beta=1
//
if (onpe0 && ct.verbose > 2)
std::cout << "Target energy..." << std::endl;
std::cout << "MVP --- Target energy..." << std::endl;
proj_mat_work_->setDM(target, orbitals.getIterativeIndex());
proj_mat_work_->computeOccupationsFromDM();
if (ct.verbose > 2) current_proj_mat->printOccupations(os_);
double nel = proj_mat_work_->getNel();
const double nel = proj_mat_work_->getNel();
if (onpe0 && ct.verbose > 1)
os_ << "Number of electrons at beta=1 : " << nel
os_ << "MVP --- Number of electrons at beta=1 : " << nel
<< std::endl;

// if( onpe0 )os_<<"Rho..."<<endl;
rho_->computeRho(orbitals, target);
rho_->rescaleTotalCharge();

mgmol_strategy_->update_pot(vh_init, ions_);

Expand All @@ -323,13 +301,13 @@ int MVPSolver<OrbitalsType, MatrixType>::solve(OrbitalsType& orbitals)
ts1, proj_mat_work_, orbitals, ct.verbose - 1, os_);

// line minimization
double beta
const double beta
= minQuadPolynomial(e0, e1, de0, (ct.verbose > 2), os_);

if (onpe0 && ct.verbose > 0)
{
os_ << std::setprecision(12);
os_ << std::fixed << "Inner iteration " << inner_it
os_ << std::fixed << "MVP inner iteration " << inner_it
<< ", E0=" << e0 << ", E1=" << e1;
os_ << std::scientific << ", E0'=" << de0
<< " -> beta=" << beta;
Expand All @@ -355,14 +333,13 @@ int MVPSolver<OrbitalsType, MatrixType>::solve(OrbitalsType& orbitals)
os_ << "Number of electrons for interpolated DM = "
<< pnel << std::endl;
}

// if( onpe0 )os_<<"Rho..."<<endl;
rho_->computeRho(orbitals, *work_);
rho_->rescaleTotalCharge();
}

} // inner iterations

// set DM to latest iteration value
ProjectedMatrices<MatrixType>* projmatrices
= dynamic_cast<ProjectedMatrices<MatrixType>*>(
orbitals.getProjMatrices());
Expand All @@ -375,11 +352,9 @@ int MVPSolver<OrbitalsType, MatrixType>::solve(OrbitalsType& orbitals)
if (onpe0 && ct.verbose > 1)
{
os_ << "---------------------------------------------------------------"
"-"
<< std::endl;
os_ << "End MVP Solver..." << std::endl;
os_ << "---------------------------------------------------------------"
"-"
<< std::endl;
}
solve_tm_.stop();
Expand Down
15 changes: 5 additions & 10 deletions src/MVPSolver.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,27 +24,22 @@ template <class OrbitalsType, class MatrixType>
class MVPSolver
{
private:
MPI_Comm comm_;
const MPI_Comm comm_;
std::ostream& os_;

short n_inner_steps_;
const short n_inner_steps_;

bool use_old_dm_;
const bool use_old_dm_;
Ions& ions_;

int numst_;

Rho<OrbitalsType>* rho_;
Energy<OrbitalsType>* energy_;
Electrostatic* electrostat_;

int history_length_;
std::vector<double> eks_history_;

MGmol<OrbitalsType>* mgmol_strategy_;

double de_old_;
double de_;

int numst_;
MatrixType* work_;
ProjectedMatrices<MatrixType>* proj_mat_work_;

Expand Down
Loading

0 comments on commit 29c2502

Please sign in to comment.