Skip to content

Commit 7198ec1

Browse files
JGHan7PeizeLinpre-commit-ci-lite[bot]mohanchenmaki49
authored
add RDMFT(Reduced Density Matrix Functional Theory) code (#5325)
* in last commit, 2 processor has 'malloc(): memory corruption:' because I do three for() to realize gemm and it doesn't support many processors... * first push with rdmft writed in abacus-v3.4.4 * update README.md * update README.md again * find all parameters used in Veff_rdmft (for hartree and local, XC still has two parameters don't pass), and pass in rdmft_cal() successfully * use class Veff_rdmft get correct Ehartree and Eexx, but Etotal is wrong * all energy terms correct, ETV_RDMFT = E_one_elec, Ehartree_RDMFT = E_Hartree, EXC_RDMFT = E_exx, but my Etotal doesn't include ion energy terms, etc. * z/dgemm_() is used correctly, next step to use pz/dgemm_() * just one processor is correct, 2 processors have some problem and maybe the error is occurred in processor 2 * 2 processors is wrong but 6 processors is correct... * check 2 processors error * update, realize print Dmn in different processors * get successful result in different numbers of processors, and the printout has not been deleted yet * get the sum of different processors, that is, get the true total energy of each term * delete some printout and next step is replace the pz/dgemm_() used in psiDotPsi() * add some comments * do some little modifications to make the code more concise * just save * just save * last save before realizing gamma_only calculation * gamma_only calculation with errors * save with errors in gamma_only * gamma_only realization with errors in E_T+nonlocal * gamma_only successful with one processors, but multi-processor errors * support g(wg)*H_wfc and g(wg)*wfcHwfc with different g(wg) which is provided by HF_XCfunc or power_XCfunc * get the special DM used in constructing V_XC, and can be used for split_m2D_ktoR() * get the special DM used in constructing V_XC and run successfully (not in numerica meaning) * modify the namespace of rdmft_test files from hamilt to rdmft * error with using spilt_m2D_ktoR() * successfully get Ds_XC, next step is update abacus to have the new interface * successfully run rdmft with power functional in multi-k calculation, have tested alpha=1.0, 0.95, 0.5 * do some little modifications to make the code more concise * fix interface provided by class OperatorEXX but having errors when do V_XC.contributeHk(ik) * print data of tensor from Ds_XC * just save * check the length and data of LM->Hexxc and Vxc_fromRI_c.Hexxs * have modified op_exx_lcao.h/.cpp, RI_2D_Comm.h/.cpp files, delete the modifying in the future * save * successfullu add power XC_functional and run correctly with Si soilds in 6 nodes and 6 processors * successfully realize power functional with testing Si solids, and some printout need to be deleted * Print: last push before delete many printout * delete modifying in source code * Refactor: first add class RDMFT * save * get DM used for updating charge/rho * try to initial DMR of class DensityMatrix * save before modifying gint.cpp * update charge using wg&wfc which provided by me successfully * delete some comment * last save before transfer wg -> natural occupation numbers (occupation_num) * successfully transfer wg -> occ_number and get correct result in Muller_XC * modifying some variable names, function names and comments (for wg -> occNum) * add module_rdmft and copy code here working successfully * delete rdmft_test.h/.cpp files in module_hamilt_lcao/... and modify corresponding CMakeLists.txt * add rdmft.h, rdmft.cpp, rdmft_tools.cpp files to create class RDMFT * just save. In the process of factor function rdmft_cal() -> class RDMFT * modify rdmft_solver in esolver_ks_lcao.h/.cpp * push with compile error in cal_dm.h( psiMulPsiMpi() function with the same name is wrong ) * add class rdmft and do rdmft_solver.init() successfully, next step is init V_TV, V_XC... * successfully move functions of class RDMFT to rdmft.cpp * just save * just save * factor function rdmft_cal() -> class RDMFT and with running error * just save * Factor: transfer function rdmft_cal() -> class RDMFT and run in afterscf() successfully * just save * now get V_local in eiter, next step move it to istep, i.e. together with T_nonlocal * move the obtaining of V_local from elec-step to ion-step * successfully find where to do update_ion() * class RDMFT get correct E_TV and E_hartree, but E_XC=0.0 * Factor: complete func rdmft_cal() -> class RDMFT and find where to use the RDMFT class object rdmft_solver appropriately * save * save for debug * Fix: fix bug in HkPsi(), transfer 'T'->'C' * save * Fix: fix nk_total in class RDMFT depend on the Symmetry::symm_flag * save * test different branch * add Dell_abaInstall.sh for Dell_server * modify module_rdmft to latest version * just save for updating branch rdmft * update module_xc for rdmft(WP22,CWP22 and so on) * Fix memory bug in XC_Functional_Libxc::cal_gdr() * Refactor XC_Functional_Libxc::convert_vtxc_v() * merge to abacus-v3.8.0 but don't debug(many, because many functions and interfaces have changed). The scaling_factor_xc issue is not handled * merge to abacus-v3.8.0 and add scaling_factor in libxc_cpp files * merge to abacus-v3.8.0 and debug is not finished * just save * just save * after merging into abacus-v3.8.0, the first compilation passed, but the program has not been run yet * after merging to abacus-v3.8.0, the program ran successfully for the first time, but the call of V_exx was not correct. * just save * just save * just save * the Vexx may have error when symmetry==1, skip it and continue do my own code * add parameters used in rdmft in input_parameters.h and so on * merge all rdmft-code to abacus-v3.8.0, next step is debug and test * debug now ,some energy is little different from abacus * merge to the latest aABACUS * Fix the bug caused by merging the latest ABAUCS, need to test other bugs such as wp22, cwp22 and the difference of energy * just save * add XC-functional BLYP_LIBXC and BLYP_LR to test * stop tracking abaInstall_HZWpara.sh * just save * just save * save * just save * modify cmake/FindELPA.cmake * save * add RDMFT(Reduced Density Matrix Functional Theory) code * [pre-commit.ci lite] apply automatic fixes * delete lib64/ and include/ * delete some production code and modify some variable/function names * save * merge to the latest abacus * just save * add tests for rdmft * recover accidentally deleted files * move rdmft related functions and their calls from esolver_ks to esolver_ks_lcao, and so on * save * save * [pre-commit.ci lite] apply automatic fixes * add '#ifdef __MPI', delete run_rdmft() and so on * just save * save * debug * add if(ENABLE_LCAO) in rdmft/CMakeLists.txt * save * save * fix compile problem * save * modify the error in module_io/test/read_input_ptest.cpp * prepare for PR * save * add symmetry calculation of exx in rdmft * add test in esolver_ks_lcao.cpp, function after_scf(): update exx to test the exx_energy get by rdmft is exactly equal to exx_energy by abacus * delete rdmft_test.h/.cpp * delete the test for exx in esolver_ks_lcao.cpp, after_scf( * update test with symmetry, tests/integrate/1001_NO_Si2_dzp_rdmft * save * [pre-commit.ci lite] apply automatic fixes * save * delete the new ord_file in tests/ * last save befor move cal_V_ TV/hartree/XC from rdmft.cpp to cal_V_rdmft.cpp * last save before move the functions about update_state from rdmft.cpp to update_state_rdmft.cpp * save * save * [pre-commit.ci lite] apply automatic fixes * save * modify ab_initio_type to bool rdmft * save * Fix: let nk_total *= nspin, in rdmft * modified according to Mr. Chen's suggestion * save * save * merge the latest abacus * [pre-commit.ci lite] apply automatic fixes * save * add some tests for rdmft * save * save * [pre-commit.ci lite] apply automatic fixes * reduce calculation time of tests: rdmft * modify some comments * modify tests * save * modify class RDMFT * save * modify Comments * save * save * save * save * save * modify * save * [pre-commit.ci lite] apply automatic fixes * add a comment in rdmft_cal_hamilt/V.cpp * rename rdmft_cal_hamilt.cpp to rdmft_pot_mat.cpp * Refactore psiDotPsi() * save * save * save * save --------- Co-authored-by: linpz <[email protected]> Co-authored-by: pre-commit-ci-lite[bot] <117423508+pre-commit-ci-lite[bot]@users.noreply.github.com> Co-authored-by: Mohan Chen <[email protected]> Co-authored-by: maki49 <[email protected]> Co-authored-by: wqzhou <[email protected]>
1 parent f084d42 commit 7198ec1

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

56 files changed

+2728
-36
lines changed

‎CMakeLists.txt‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -725,7 +725,8 @@ if(ENABLE_LCAO)
725725
hcontainer
726726
deltaspin
727727
numerical_atomic_orbitals
728-
lr)
728+
lr
729+
rdmft)
729730
if(USE_ELPA)
730731
target_link_libraries(${ABACUS_BIN_NAME} genelpa)
731732
endif()

‎docs/advanced/input_files/input-main.md‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -435,6 +435,9 @@
435435
- [abs\_broadening](#abs_broadening)
436436
- [ri\_hartree\_benchmark](#ri_hartree_benchmark)
437437
- [aims\_nbasis](#aims_nbasis)
438+
- [Reduced Density Matrix Functional Theory](#Reduced-Density-Matrix-Functional-Theory)
439+
- [rdmft](#rdmft)
440+
- [rdmft\_power\_alpha](#rdmft_power_alpha)
438441

439442
[back to top](#full-list-of-input-keywords)
440443
## System variables
@@ -4031,4 +4034,21 @@ The output files are `OUT.${suffix}/Excitation_Energy.dat` and `OUT.${suffix}/Ex
40314034
- **Description**: Atomic basis set size for each atom type (with the same order as in `STRU`) in FHI-aims.
40324035
- **Default**: {} (empty list, where ABACUS use its own basis set size)
40334036

4037+
## Reduced Density Matrix Functional Theory
4038+
4039+
ab-initio methods and the xc-functional parameters used in RDMFT.
4040+
The physical quantities that RDMFT temporarily expects to output are the kinetic energy, total energy, and 1-RDM of the system in the ground state, etc.
4041+
4042+
### rdmft
4043+
4044+
- **Type**: Boolean
4045+
- **Description**: Whether to perform rdmft calculation (reduced density matrix funcional theory)
4046+
- **Default**: false
4047+
4048+
### rdmft_power_alpha
4049+
4050+
- **Type**: Real
4051+
- **Description**: The alpha parameter of power-functional(or other exx-type/hybrid functionals) which used in RDMFT, g(occ_number) = occ_number^alpha
4052+
- **Default**: 0.656
4053+
40344054
[back to top](#full-list-of-input-keywords)

‎source/CMakeLists.txt‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,9 @@ add_subdirectory(module_ri)
1818
add_subdirectory(module_parameter)
1919
add_subdirectory(module_lr)
2020

21+
# add by jghan
22+
add_subdirectory(module_rdmft)
23+
2124
add_library(
2225
driver
2326
OBJECT

‎source/Makefile.Objects‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -74,6 +74,7 @@ VPATH=./src_global:\
7474
./module_lr/operator_casida:\
7575
./module_lr/potentials:\
7676
./module_lr/utils:\
77+
./module_rdmft:\
7778
./\
7879

7980
OBJS_ABACUS_PW=${OBJS_MAIN}\
@@ -115,6 +116,7 @@ ${OBJS_DELTASPIN}\
115116
${OBJS_TENSOR}\
116117
${OBJS_HSOLVER_PEXSI}\
117118
${OBJS_LR}\
119+
${OBJS_RDMFT}
118120

119121
OBJS_MAIN=main.o\
120122
driver.o\
@@ -726,3 +728,8 @@ OBJS_TENSOR=tensor.o\
726728
lr_spectrum.o\
727729
hamilt_casida.o\
728730
esolver_lrtd_lcao.o\
731+
732+
OBJS_RDMFT=rdmft.o\
733+
rdmft_tools.o\
734+
rdmft_pot.o\
735+
update_state_rdmft.o\

‎source/module_esolver/esolver_ks_lcao.cpp‎

Lines changed: 30 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -63,6 +63,10 @@
6363
// #include "module_elecstate/cal_dm.h"
6464
//---------------------------------------------------
6565

66+
// test RDMFT
67+
#include "module_rdmft/rdmft.h"
68+
#include <iostream>
69+
6670
namespace ModuleESolver
6771
{
6872

@@ -250,6 +254,13 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
250254
}
251255
}
252256

257+
// 14) initialize rdmft, added by jghan
258+
if( PARAM.inp.rdmft == true )
259+
{
260+
rdmft_solver.init( this->GG, this->GK, this->pv, ucell, this->kv, *(this->pelec),
261+
this->orb_, two_center_bundle_, PARAM.inp.dft_functional, PARAM.inp.rdmft_power_alpha);
262+
}
263+
253264
ModuleBase::timer::tick("ESolver_KS_LCAO", "before_all_runners");
254265
return;
255266
}
@@ -839,6 +850,7 @@ void ESolver_KS_LCAO<TK, TR>::update_pot(UnitCell& ucell, const int istep, const
839850
{
840851
this->pelec->cal_converged();
841852
}
853+
842854
}
843855

844856
//------------------------------------------------------------------------------
@@ -1042,6 +1054,23 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(UnitCell& ucell, const int istep)
10421054
ModuleBase::timer::tick("ESolver_KS_LCAO", "out_deepks_labels");
10431055
#endif
10441056

1057+
/******** test RDMFT *********/
1058+
if ( PARAM.inp.rdmft == true ) // rdmft, added by jghan, 2024-10-17
1059+
{
1060+
ModuleBase::matrix occ_number_ks(this->pelec->wg);
1061+
for(int ik=0; ik < occ_number_ks.nr; ++ik) { for(int inb=0; inb < occ_number_ks.nc; ++inb) occ_number_ks(ik, inb) /= this->kv.wk[ik]; }
1062+
this->rdmft_solver.update_elec(occ_number_ks, *(this->psi));
1063+
1064+
//initialize the gradients of Etotal on occupation numbers and wfc, and set all elements to 0.
1065+
ModuleBase::matrix dE_dOccNum(this->pelec->wg.nr, this->pelec->wg.nc, true);
1066+
psi::Psi<TK> dE_dWfc(this->psi->get_nk(), this->psi->get_nbands(), this->psi->get_nbasis());
1067+
dE_dWfc.zero_out();
1068+
1069+
double Etotal_RDMFT = this->rdmft_solver.run(dE_dOccNum, dE_dWfc);
1070+
}
1071+
/******** test RDMFT *********/
1072+
1073+
10451074
#ifdef __EXX
10461075
// 11) write rpa information
10471076
if (PARAM.inp.rpa)
@@ -1228,6 +1257,7 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(UnitCell& ucell, const int istep)
12281257
}
12291258
}
12301259

1260+
12311261
//------------------------------------------------------------------------------
12321262
//! the 20th,21th,22th functions of ESolver_KS_LCAO
12331263
//! mohan add 2024-05-11

‎source/module_esolver/esolver_ks_lcao.h‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,9 @@
1212
#include "module_basis/module_nao/two_center_bundle.h"
1313
#include "module_io/output_mat_sparse.h"
1414

15+
// added by jghan for rdmft calculation
16+
#include "module_rdmft/rdmft.h"
17+
1518
#include <memory>
1619

1720
namespace LR
@@ -68,6 +71,8 @@ class ESolver_KS_LCAO : public ESolver_KS<TK> {
6871

6972
TwoCenterBundle two_center_bundle_;
7073

74+
rdmft::RDMFT<TK, TR> rdmft_solver; // added by jghan for rdmft calculation, 2024-03-16
75+
7176
// temporary introduced during removing GlobalC::ORB
7277
LCAO_Orbitals orb_;
7378

‎source/module_esolver/lcao_before_scf.cpp‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -370,6 +370,14 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
370370
}
371371

372372
this->p_hamilt->non_first_scf = istep;
373+
374+
// update in ion-step
375+
if( PARAM.inp.rdmft == true )
376+
{
377+
// necessary operation of these parameters have be done with p_esolver->Init() in source/driver_run.cpp
378+
rdmft_solver.update_ion(ucell, *(this->pw_rho), GlobalC::ppcell.vloc, this->sf.strucFac); // add by jghan, 2024-03-16/2024-10-08
379+
}
380+
373381
return;
374382
}
375383

‎source/module_hamilt_general/module_xc/xc_functional.cpp‎

Lines changed: 51 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,7 @@ std::vector<int> XC_Functional::func_id(1);
1818
int XC_Functional::func_type = 0;
1919
bool XC_Functional::use_libxc = true;
2020
double XC_Functional::hybrid_alpha = 0.25;
21+
std::map<int, double> XC_Functional::scaling_factor_xc = { {1, 1.0} }; // added by jghan, 2024-10-10
2122

2223
void XC_Functional::set_hybrid_alpha(const double alpha_in)
2324
{
@@ -61,6 +62,7 @@ void XC_Functional::set_xc_type(const std::string xc_func_in)
6162
// func_id.push_back(XC_GGA_C_PBE);
6263

6364
func_id.clear();
65+
scaling_factor_xc.clear(); // added by jghan, 2024-07-07
6466
std::string xc_func = xc_func_in;
6567
std::transform(xc_func.begin(), xc_func.end(), xc_func.begin(), (::toupper));
6668
if( xc_func == "LDA" || xc_func == "PZ" || xc_func == "SLAPZNOGXNOGC") //SLA+PZ
@@ -196,13 +198,60 @@ void XC_Functional::set_xc_type(const std::string xc_func_in)
196198
{
197199
// not doing anything
198200
}
201+
else if( xc_func == "MULLER" || xc_func == "POWER" ) // added by jghan, 2024-07-06
202+
{
203+
func_type = 4;
204+
use_libxc = false;
205+
}
199206
#ifdef USE_LIBXC
200207
else if( xc_func == "HSE")
201208
{
202209
func_id.push_back(XC_HYB_GGA_XC_HSE06);
203210
func_type = 4;
204211
use_libxc = true;
205212
}
213+
// added by jghan, 2024-07-06
214+
else if( xc_func == "WP22")
215+
{
216+
func_id.push_back(XC_GGA_X_ITYH); // short-range of B88_X, id=529
217+
func_id.push_back(XC_GGA_C_LYPR); // short-range of LYP_C, id=624
218+
func_type = 4;
219+
use_libxc = true;
220+
}
221+
else if( xc_func == "CWP22")
222+
{
223+
// BLYP_XC_lr = -BLYP_XC_sr + BLYP_XC, the realization of it is in v_xc_libxc() function, xc_functional_libxc_vxc.cpp
224+
func_id.push_back(XC_GGA_X_ITYH); // short-range of B88_X, id=529
225+
func_id.push_back(XC_GGA_C_LYPR); // short-range of LYP_C, id=624
226+
func_id.push_back(XC_GGA_X_B88); // complete B88_X, id=106
227+
func_id.push_back(XC_GGA_C_LYP); // complete LYP_C, id=131
228+
229+
// the scaling factor of CWP22-functionals
230+
scaling_factor_xc[XC_GGA_X_ITYH] = -1.0;
231+
scaling_factor_xc[XC_GGA_C_LYPR] = -1.0;
232+
scaling_factor_xc[XC_GGA_X_B88] = 1.0;
233+
scaling_factor_xc[XC_GGA_X_B88] = 1.0;
234+
235+
func_type = 4;
236+
use_libxc = true;
237+
}
238+
else if( xc_func == "BLYP_LR")
239+
{
240+
// BLYP_XC_lr = -BLYP_XC_sr + BLYP_XC, the realization of it is in v_xc_libxc() function, xc_functional_libxc_vxc.cpp
241+
func_id.push_back(XC_GGA_X_ITYH); // short-range of B88_X, id=529
242+
func_id.push_back(XC_GGA_C_LYPR); // short-range of LYP_C, id=624
243+
func_id.push_back(XC_GGA_X_B88); // complete B88_X, id=106
244+
func_id.push_back(XC_GGA_C_LYP); // complete LYP_C, id=131
245+
246+
// the scaling factor of BLYP_LR-functionals
247+
scaling_factor_xc[XC_GGA_X_ITYH] = -1.0;
248+
scaling_factor_xc[XC_GGA_C_LYPR] = -1.0;
249+
scaling_factor_xc[XC_GGA_X_B88] = 1.0;
250+
scaling_factor_xc[XC_GGA_X_B88] = 1.0;
251+
252+
func_type = 2;
253+
use_libxc = true;
254+
}
206255
#endif
207256
else
208257
{
@@ -243,7 +292,8 @@ void XC_Functional::set_xc_type(const std::string xc_func_in)
243292
#endif
244293

245294
#ifndef USE_LIBXC
246-
if(xc_func == "SCAN" || xc_func == "HSE" || xc_func == "SCAN0")
295+
if(xc_func == "SCAN" || xc_func == "HSE" || xc_func == "SCAN0"
296+
|| xc_func == "MULLER" || xc_func == "POWER" || xc_func == "WP22" || xc_func == "CWP22")
247297
{
248298
ModuleBase::WARNING_QUIT("set_xc_type","to use SCAN, SCAN0, or HSE, LIBXC is required");
249299
}

‎source/module_hamilt_general/module_xc/xc_functional.h‎

Lines changed: 7 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -20,6 +20,8 @@
2020
#include "module_elecstate/module_charge/charge.h"
2121
#include "module_cell/unitcell.h"
2222

23+
#include <map> // added by jghan, 2024-10-10
24+
2325
class XC_Functional
2426
{
2527
public:
@@ -80,6 +82,10 @@ class XC_Functional
8082
//exx_hybrid_alpha for mixing exx in hybrid functional:
8183
static double hybrid_alpha;
8284

85+
// added by jghan, 2024-07-07
86+
// as a scaling factor for different xc-functionals
87+
static std::map<int, double> scaling_factor_xc;
88+
8389
public:
8490
static std::vector<int> get_func_id() { return func_id; }
8591

@@ -162,7 +168,7 @@ class XC_Functional
162168
ModulePW::PW_Basis* rhopw,
163169
const UnitCell* ucell,
164170
std::vector<double>& stress_gga,
165-
const bool is_stress = 0);
171+
const bool is_stress = false);
166172
template <typename T, typename Device,
167173
typename Real = typename GetTypeReal<T>::type>
168174
static void grad_wfc(

‎source/module_hamilt_general/module_xc/xc_functional_libxc.cpp‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -135,6 +135,20 @@ std::vector<xc_func_type> XC_Functional_Libxc::init_func(const std::vector<int>
135135
GlobalC::exx_info.info_global.hse_omega };
136136
xc_func_set_ext_params(&funcs.back(), parameter_hse);
137137
}
138+
// added by jghan, 2024-07-06
139+
else if( id == XC_GGA_X_ITYH ) // short-range of B88_X
140+
{
141+
add_func( XC_GGA_X_ITYH );
142+
double parameter_omega[1] = {PARAM.inp.exx_hse_omega}; // GlobalC::exx_info.info_global.hse_omega
143+
xc_func_set_ext_params(&funcs.back(), parameter_omega);
144+
}
145+
else if( id == XC_GGA_C_LYPR ) // short-range of LYP_C
146+
{
147+
add_func( XC_GGA_C_LYPR );
148+
// the first six parameters come from libxc, and may need to be modified in some cases
149+
double parameter_lypr[7] = {0.04918, 0.132, 0.2533, 0.349, 0.35/2.29, 2.0/2.29, PARAM.inp.exx_hse_omega};
150+
xc_func_set_ext_params(&funcs.back(), parameter_lypr);
151+
}
138152
#endif
139153
else
140154
{

0 commit comments

Comments
 (0)