Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fix: DeePKS MD/Relax/Cell-relax calculation convergence #4575

Merged
merged 3 commits into from
Jul 7, 2024
Merged

Fix: DeePKS MD/Relax/Cell-relax calculation convergence #4575

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
c31372a
Fix: DeePKS MD/Relax/Cell-relax calculation convergence
Jul 1, 2024
16eb176
Fix: compiling error
Jul 1, 2024
01d9579
Merge branch 'develop' of github.com:deepmodeling/abacus-develop into…
Jul 4, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
19 changes: 14 additions & 5 deletions source/module_esolver/esolver_ks_lcao.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -244,6 +244,8 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(Input& inp, UnitCell& ucell) {
if (GlobalV::deepks_scf) {
// load the DeePKS model from deep neural network
GlobalC::ld.load_model(INPUT.deepks_model);
// read pdm from file for NSCF or SCF-restart, do it only once in whole calculation
GlobalC::ld.read_projected_DM((GlobalV::init_chg == "file"), GlobalV::deepks_equiv, *GlobalC::ORB.Alpha);
}
#endif

Expand Down Expand Up @@ -1294,6 +1296,18 @@ template <typename TK, typename TR>
bool ESolver_KS_LCAO<TK, TR>::do_after_converge(int& iter) {
ModuleBase::TITLE("ESolver_KS_LCAO", "do_after_converge");

if (GlobalV::dft_plus_u)
{
// use the converged occupation matrix for next MD/Relax SCF calculation
GlobalC::dftu.initialed_locale = true;
}

#ifdef __DEEPKS
if (GlobalV::deepks_scf)
{
GlobalC::ld.set_init_pdm(true);
}
#endif
#ifdef __EXX
if (GlobalC::exx_info.info_ri.real_number) {
return this->exd->exx_after_converge(
Expand All @@ -1314,11 +1328,6 @@ bool ESolver_KS_LCAO<TK, TR>::do_after_converge(int& iter) {
}
#endif // __EXX

if (GlobalV::dft_plus_u) {
// use the converged occupation matrix for next MD/Relax SCF calculation
GlobalC::dftu.initialed_locale = true;
}

return true;
}

Expand Down
3 changes: 2 additions & 1 deletion source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -277,7 +277,8 @@ void Force_LCAO<double>::ftable(const bool isforce,
{
const std::vector<std::vector<double>>& dm_gamma = dm->get_DMK_vector();

GlobalC::ld.cal_projected_DM(dm, ucell, GlobalC::ORB, GlobalC::GridD);
// when deepks_scf is on, the init pdm should be same as the out pdm, so we should not recalculate the pdm
//GlobalC::ld.cal_projected_DM(dm, ucell, GlobalC::ORB, GlobalC::GridD);

GlobalC::ld.cal_descriptor(ucell.nat);

Expand Down
3 changes: 2 additions & 1 deletion source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -355,7 +355,8 @@ void Force_LCAO<std::complex<double>>::ftable(const bool isforce,
{
const std::vector<std::vector<std::complex<double>>>& dm_k = dm->get_DMK_vector();

GlobalC::ld.cal_projected_DM_k(dm, ucell, GlobalC::ORB, GlobalC::GridD);
// when deepks_scf is on, the init pdm should be same as the out pdm, so we should not recalculate the pdm
//GlobalC::ld.cal_projected_DM_k(dm, ucell, GlobalC::ORB, GlobalC::GridD);

GlobalC::ld.cal_descriptor(ucell.nat);

Expand Down
22 changes: 20 additions & 2 deletions source/module_hamilt_lcao/module_deepks/LCAO_deepks.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -278,8 +278,14 @@ class LCAO_Deepks
// 6. check_gdmx, which prints gdmx to a series of .dat files

public:
/// calculate projected density matrix:
/// pdm = sum_i,occ <phi_i|alpha1><alpha2|phi_k>
/**
* @brief calculate projected density matrix:
* pdm = sum_i,occ <phi_i|alpha1><alpha2|phi_k>
* 3 cases to skip calculation of pdm:
* 1. NSCF calculation of DeePKS, init_chg = file and pdm has been read
* 2. SCF calculation of DeePKS with init_chg = file and pdm has been read for restarting SCF
* 3. Relax/Cell-Relax/MD calculation, non-first step will use the convergence pdm from the last step as initial pdm
*/
void cal_projected_DM(const elecstate::DensityMatrix<double, double>* dm,
const UnitCell& ucell,
const LCAO_Orbitals& orb,
Expand Down Expand Up @@ -317,6 +323,18 @@ class LCAO_Deepks
const bool isstress);
void check_gdmx(const int nat);

/**
* @brief set init_pdm to skip the calculation of pdm in SCF iteration
*/
void set_init_pdm(bool ipdm)
{
this->init_pdm = ipdm;
}
/**
* @brief read pdm from file, do it only once in whole calculation
*/
void read_projected_DM(bool read_pdm_file, bool is_equiv, const Numerical_Orbital& alpha);

//-------------------
// LCAO_deepks_vdelta.cpp
//-------------------
Expand Down
6 changes: 4 additions & 2 deletions source/module_hamilt_lcao/module_deepks/LCAO_deepks_interface.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -86,7 +86,8 @@ void LCAO_Deepks_Interface::out_deepks_labels(double etot,
if (GlobalV::deepks_out_labels || GlobalV::deepks_scf)
{
// this part is for integrated test of deepks
ld->cal_projected_DM(dm, ucell, orb, GridD);
// when deepks_scf is on, the init pdm should be same as the out pdm, so we should not recalculate the pdm
if(!GlobalV::deepks_scf) ld->cal_projected_DM(dm, ucell, orb, GridD);
ld->check_projected_dm(); // print out the projected dm for NSCF calculaiton
ld->cal_descriptor(nat); // final descriptor
ld->check_descriptor(ucell);
Expand Down Expand Up @@ -188,7 +189,8 @@ void LCAO_Deepks_Interface::out_deepks_labels(double etot,
{
// this part is for integrated test of deepks
// so it is printed no matter even if deepks_out_labels is not used
ld->cal_projected_DM_k(dm, ucell, orb, GridD);
// when deepks_scf is on, the init pdm should be same as the out pdm, so we should not recalculate the pdm
if(!GlobalV::deepks_scf) ld->cal_projected_DM_k(dm, ucell, orb, GridD);
ld->check_projected_dm(); // print out the projected dm for NSCF calculaiton
ld->cal_descriptor(nat); // final descriptor
ld->check_descriptor(ucell);
Expand Down
87 changes: 48 additions & 39 deletions source/module_hamilt_lcao/module_deepks/LCAO_deepks_pdm.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,42 @@
#include "module_hamilt_lcao/module_hcontainer/atom_pair.h"
#include "module_base/libm/libm.h"

void LCAO_Deepks::read_projected_DM(bool read_pdm_file, bool is_equiv, const Numerical_Orbital& alpha)
{
if (read_pdm_file && !this->init_pdm) //for DeePKS NSCF calculation
{
int pdm_size = 0;
if(!is_equiv)
{
pdm_size = (this->lmaxd * 2 + 1) * (this->lmaxd * 2 + 1);
}
else
{
int nproj = 0;
for(int il = 0; il < this->lmaxd + 1; il++)
{
nproj += (2 * il + 1) * alpha.getNchi(il);
}
pdm_size = nproj * nproj;
}
std::ifstream ifs("pdm.dat");
if (!ifs)
{
ModuleBase::WARNING_QUIT("LCAO_Deepks::cal_projected_DM", "Can not find the file pdm.dat . Please do DeePKS SCF calculation first.");
}
for(int inl=0;inl<this->inlmax;inl++)
{
for(int ind=0;ind<pdm_size;ind++)
{
double c;
ifs >> c;
pdm[inl][ind] = c;
}
}
this->init_pdm = true;
}
}

//this subroutine performs the calculation of projected density matrices
//pdm_m,m'=\sum_{mu,nu} rho_{mu,nu} <chi_mu|alpha_m><alpha_m'|chi_nu>
void LCAO_Deepks::cal_projected_DM(const elecstate::DensityMatrix<double, double>* dm,
Expand All @@ -34,6 +70,12 @@ void LCAO_Deepks::cal_projected_DM(const elecstate::DensityMatrix<double, double
{
ModuleBase::TITLE("LCAO_Deepks", "cal_projected_DM");

// if pdm has been initialized, skip the calculation
if(this->init_pdm)
{
this->init_pdm = false;
return;
}
int pdm_size = 0;
if(!GlobalV::deepks_equiv)
{
Expand All @@ -48,25 +90,6 @@ void LCAO_Deepks::cal_projected_DM(const elecstate::DensityMatrix<double, double
}
pdm_size = nproj * nproj;
}
if (GlobalV::init_chg == "file" && !this->init_pdm) //for DeePKS NSCF calculation
{
std::ifstream ifs("pdm.dat");
if (!ifs)
{
ModuleBase::WARNING_QUIT("LCAO_Deepks::cal_projected_DM", "Can not find the file pdm.dat . Please do DeePKS SCF calculation first.");
}
for(int inl=0;inl<this->inlmax;inl++)
{
for(int ind=0;ind<pdm_size;ind++)
{
double c;
ifs >> c;
pdm[inl][ind] = c;
}
}
this->init_pdm = true;
return;
}

//if(dm.size() == 0 || dm[0].size() == 0)
//{
Expand Down Expand Up @@ -289,6 +312,12 @@ void LCAO_Deepks::cal_projected_DM_k(const elecstate::DensityMatrix<std::complex
const LCAO_Orbitals &orb,
Grid_Driver& GridD)
{
// if pdm has been initialized, skip the calculation
if(this->init_pdm)
{
this->init_pdm = false;
return;
}

int pdm_size = 0;
if(!GlobalV::deepks_equiv)
Expand All @@ -305,26 +334,6 @@ void LCAO_Deepks::cal_projected_DM_k(const elecstate::DensityMatrix<std::complex
pdm_size = nproj * nproj;
}

if (GlobalV::init_chg == "file" && !this->init_pdm) //for DeePKS NSCF calculation
{
std::ifstream ifs("pdm.dat");
if (!ifs)
{
ModuleBase::WARNING_QUIT("LCAO_Deepks::cal_projected_DM_k","Can not find the file pdm.dat . Please do DeePKS SCF calculation first.");
}
for(int inl=0;inl<this->inlmax;inl++)
{
for(int ind=0;ind<pdm_size;ind++)
{
double c;
ifs >> c;
pdm[inl][ind] = c;
}
}
this->init_pdm = true;
return;
}

//check for skipping
//if(dm.size() == 0 || dm[0].size() == 0)
//{
Expand Down
Loading