Skip to content

Commit 681a7a7

Browse files
maki49maki49
committed
rename init_fxc as init_xc_kernel
1 parent 7157fd0 commit 681a7a7

File tree

9 files changed

+341
-25
lines changed

9 files changed

+341
-25
lines changed

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

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -3921,14 +3921,14 @@ These parameters are used to solve the excited states using. e.g. LR-TDDFT.
39213921
Currently supported: `RPA`, `LDA`, `PBE`, `HSE`, `HF`.
39223922
- **Default**: LDA
39233923

3924-
### init_fxc
3924+
### init_xc_kernel
39253925

39263926
- **Type**: String
39273927
- **Description**: The method to initalize the xc kernel.
3928-
- "gs": Calculate xc kerenel ($f_\text{xc}$) from the ground-state charge density.
3929-
- "file_fxc": Read the xc kernel $f_\text{xc}$ on grid from the provided files. The following words should be the paths of ".cube" files, where the first 3 (spin-aa, spin-ab and spin-bb) will be read in. The parameter [xc_kernel](#xc_kernel) will be invalid. Now only LDA-type kernel is supproted as the potential will be calculated by directly multiplying the transition density.
3930-
- "file_chg": Calculate fxc from the charge density read from the provided files. The following words should be the paths of ".cube" files, where the first [nspin]($nspin) files will be read in.
3931-
- **Default**: "gs"
3928+
- "from_chg_groundstate": Calculate xc kerenel ($f_\text{xc}$) from the ground-state charge density.
3929+
- "file": Read the xc kernel $f_\text{xc}$ on grid from the provided files. The following words should be the paths of ".cube" files, where the first 3 (spin-aa, spin-ab and spin-bb) will be read in. The parameter [xc_kernel](#xc_kernel) will be invalid. Now only LDA-type kernel is supproted as the potential will be calculated by directly multiplying the transition density.
3930+
- "from_chg_file": Calculate fxc from the charge density read from the provided files. The following words should be the paths of ".cube" files, where the first [nspin]($nspin) files will be read in.
3931+
- **Default**: "from_chg_groundstate"
39323932

39333933
### lr_solver
39343934

‎source/module_io/read_input_item_tddft.cpp‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -310,17 +310,17 @@ void ReadInput::item_lr_tddft()
310310
this->add_item(item);
311311
}
312312
{
313-
Input_Item item("init_fxc");
313+
Input_Item item("init_xc_kernel");
314314
item.annotation = "The method to initalize the xc kernel";
315315
item.read_value = [](const Input_Item& item, Parameter& para) {
316316
size_t count = item.get_size();
317-
auto& ifxc = para.input.init_fxc;
317+
auto& ifxc = para.input.init_xc_kernel;
318318
for (int i = 0; i < count; i++) { ifxc.push_back(item.str_values[i]); }
319319
};
320320
item.reset_value = [](const Input_Item& item, Parameter& para) {
321-
if (para.input.init_fxc.empty()) { para.input.init_fxc.push_back("gs"); }
321+
if (para.input.init_xc_kernel.empty()) { para.input.init_xc_kernel.push_back("from_chg_groundstate"); }
322322
};
323-
sync_stringvec(input.init_fxc, para.input.init_fxc.size(), "gs");
323+
sync_stringvec(input.init_xc_kernel, para.input.init_xc_kernel.size(), "from_chg_groundstate");
324324
this->add_item(item);
325325
}
326326
{

‎source/module_io/test/read_input_ptest.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -417,7 +417,7 @@ TEST_F(InputParaTest, ParaRead)
417417
EXPECT_EQ(param.inp.nocc, param.inp.nbands);
418418
EXPECT_EQ(param.inp.nvirt, 1);
419419
EXPECT_EQ(param.inp.xc_kernel, "LDA");
420-
EXPECT_EQ(param.inp.init_fxc[0], "gs");
420+
EXPECT_EQ(param.inp.init_xc_kernel[0], "from_chg_groundstate");
421421
EXPECT_EQ(param.inp.lr_solver, "dav");
422422
EXPECT_DOUBLE_EQ(param.inp.lr_thr, 1e-2);
423423
EXPECT_FALSE(param.inp.lr_unrestricted);

‎source/module_io/write_libxc_r.cpp‎

Lines changed: 285 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,285 @@
1+
//======================
2+
// AUTHOR : Peize Lin
3+
// DATE : 2024-09-12
4+
//======================
5+
6+
#ifdef USE_LIBXC
7+
8+
#include "write_libxc_r.h"
9+
#include "module_hamilt_general/module_xc/xc_functional.h"
10+
#include "module_hamilt_general/module_xc/xc_functional_libxc.h"
11+
#include "module_elecstate/module_charge/charge.h"
12+
#include "module_esolver/esolver_fp.h"
13+
#include "module_io/cube_io.h"
14+
#include "module_base/global_variable.h"
15+
#include "module_parameter/parameter.h"
16+
#include "module_base/timer.h"
17+
#include "module_base/tool_title.h"
18+
19+
#include <xc.h>
20+
#include <iostream>
21+
#include <vector>
22+
#include <stdexcept>
23+
#include <unistd.h>
24+
25+
void ModuleIO::write_libxc_r(
26+
const int order,
27+
const std::vector<int> &func_id,
28+
const int &nrxx, // number of real-space grid
29+
const double &omega, // volume of cell
30+
const double tpiba,
31+
const Charge*const chr,
32+
const ModuleESolver::ESolver_FP*const esolver)
33+
{
34+
ModuleBase::TITLE("ModuleIO","write_libxc_r");
35+
ModuleBase::timer::tick("ModuleIO","write_libxc_r");
36+
37+
const int nspin =
38+
(PARAM.inp.nspin == 1 || ( PARAM.inp.nspin ==4 && !PARAM.globalv.domag && !PARAM.globalv.domag_z))
39+
? 1 : 2;
40+
41+
//----------------------------------------------------------
42+
// xc_func_type is defined in Libxc package
43+
// to understand the usage of xc_func_type,
44+
// use can check on website, for example:
45+
// https://www.tddft.org/programs/libxc/manual/libxc-5.1.x/
46+
//----------------------------------------------------------
47+
48+
std::vector<xc_func_type> funcs = XC_Functional_Libxc::init_func( func_id, (1==nspin) ? XC_UNPOLARIZED:XC_POLARIZED );
49+
50+
const bool is_gga = [&funcs]()
51+
{
52+
for( xc_func_type &func : funcs )
53+
{
54+
switch( func.info->family )
55+
{
56+
case XC_FAMILY_GGA:
57+
case XC_FAMILY_HYB_GGA:
58+
return true;
59+
}
60+
}
61+
return false;
62+
}();
63+
64+
// converting rho
65+
std::vector<double> rho;
66+
std::vector<double> amag;
67+
if(1==nspin || 2==PARAM.inp.nspin)
68+
{
69+
rho = XC_Functional_Libxc::convert_rho(nspin, nrxx, chr);
70+
}
71+
else
72+
{
73+
std::tuple<std::vector<double>,std::vector<double>> rho_amag = XC_Functional_Libxc::convert_rho_amag_nspin4(nspin, nrxx, chr);
74+
rho = std::get<0>(std::move(rho_amag));
75+
amag = std::get<1>(std::move(rho_amag));
76+
}
77+
78+
std::vector<double> sigma;
79+
if(is_gga)
80+
{
81+
const std::vector<std::vector<ModuleBase::Vector3<double>>> gdr = XC_Functional_Libxc::cal_gdr(nspin, nrxx, rho, tpiba, chr);
82+
sigma = XC_Functional_Libxc::convert_sigma(gdr);
83+
}
84+
85+
std::vector<double> exc;
86+
std::vector<double> vrho;
87+
std::vector<double> vsigma;
88+
std::vector<double> v2rho2;
89+
std::vector<double> v2rhosigma;
90+
std::vector<double> v2sigma2;
91+
std::vector<double> v3rho3;
92+
std::vector<double> v3rho2sigma;
93+
std::vector<double> v3rhosigma2;
94+
std::vector<double> v3sigma3;
95+
std::vector<double> v4rho4;
96+
std::vector<double> v4rho3sigma;
97+
std::vector<double> v4rho2sigma2;
98+
std::vector<double> v4rhosigma3;
99+
std::vector<double> v4sigma4;
100+
// attention: order 4321 don't break
101+
switch( order )
102+
{
103+
case 4: v4rho4.resize( nrxx * ((1==nspin)?1:5) );
104+
case 3: v3rho3.resize( nrxx * ((1==nspin)?1:4) );
105+
case 2: v2rho2.resize( nrxx * ((1==nspin)?1:3) );
106+
case 1: vrho .resize( nrxx * nspin );
107+
case 0: exc .resize( nrxx );
108+
break;
109+
default: throw std::domain_error("order ="+std::to_string(order)
110+
+" unfinished in "+std::string(__FILE__)+" line "+std::to_string(__LINE__));
111+
break;
112+
}
113+
if(is_gga)
114+
{
115+
switch( order )
116+
{
117+
case 4: v4rho3sigma .resize( nrxx * ((1==nspin)?1:12) );
118+
v4rho2sigma2.resize( nrxx * ((1==nspin)?1:15) );
119+
v4rhosigma3 .resize( nrxx * ((1==nspin)?1:20) );
120+
v4sigma4 .resize( nrxx * ((1==nspin)?1:15) );
121+
case 3: v3rho2sigma .resize( nrxx * ((1==nspin)?1:9) );
122+
v3rhosigma2 .resize( nrxx * ((1==nspin)?1:12) );
123+
v3sigma3 .resize( nrxx * ((1==nspin)?1:10) );
124+
case 2: v2rhosigma .resize( nrxx * ((1==nspin)?1:6) );
125+
v2sigma2 .resize( nrxx * ((1==nspin)?1:6) );
126+
case 1: vsigma .resize( nrxx * ((1==nspin)?1:3) );
127+
case 0: break;
128+
default: throw std::domain_error("order ="+std::to_string(order)
129+
+" unfinished in "+std::string(__FILE__)+" line "+std::to_string(__LINE__));
130+
break;
131+
}
132+
}
133+
134+
for( xc_func_type &func : funcs )
135+
{
136+
// jiyy add for threshold
137+
constexpr double rho_threshold = 1E-6;
138+
constexpr double grho_threshold = 1E-10;
139+
140+
xc_func_set_dens_threshold(&func, rho_threshold);
141+
142+
// sgn for threshold mask
143+
const std::vector<double> sgn = XC_Functional_Libxc::cal_sgn(rho_threshold, grho_threshold, func, nspin, nrxx, rho, sigma);
144+
145+
// call libxc function
146+
// attention: order 432 don't break
147+
switch( func.info->family )
148+
{
149+
case XC_FAMILY_LDA:
150+
{
151+
switch( order )
152+
{
153+
case 4: xc_lda_lxc ( &func, nrxx, rho.data(), v4rho4.data() );
154+
case 3: xc_lda_kxc ( &func, nrxx, rho.data(), v3rho3.data() );
155+
case 2: xc_lda_fxc ( &func, nrxx, rho.data(), v2rho2.data() );
156+
case 1: xc_lda_exc_vxc( &func, nrxx, rho.data(), exc.data(), vrho.data() );
157+
break;
158+
case 0: xc_lda_exc ( &func, nrxx, rho.data(), exc.data() );
159+
break;
160+
default: throw std::domain_error("order ="+std::to_string(order)
161+
+" unfinished in "+std::string(__FILE__)+" line "+std::to_string(__LINE__));
162+
break;
163+
}
164+
break;
165+
}
166+
case XC_FAMILY_GGA:
167+
case XC_FAMILY_HYB_GGA:
168+
{
169+
switch( order )
170+
{
171+
case 4: xc_gga_lxc ( &func, nrxx, rho.data(), sigma.data(), v4rho4.data(), v4rho3sigma.data(), v4rho2sigma2.data(), v4rhosigma3.data(), v4sigma4.data() );
172+
case 3: xc_gga_kxc ( &func, nrxx, rho.data(), sigma.data(), v3rho3.data(), v3rho2sigma.data(), v3rhosigma2.data(), v3sigma3.data() );
173+
case 2: xc_gga_fxc ( &func, nrxx, rho.data(), sigma.data(), v2rho2.data(), v2rhosigma.data(), v2sigma2.data() );
174+
case 1: xc_gga_exc_vxc( &func, nrxx, rho.data(), sigma.data(), exc.data(), vrho.data(), vsigma.data() );
175+
break;
176+
case 0: xc_gga_exc ( &func, nrxx, rho.data(), sigma.data(), exc.data() );
177+
break;
178+
default: throw std::domain_error("order ="+std::to_string(order)
179+
+" unfinished in "+std::string(__FILE__)+" line "+std::to_string(__LINE__));
180+
break;
181+
}
182+
break;
183+
}
184+
default:
185+
{
186+
throw std::domain_error("func.info->family ="+std::to_string(func.info->family)
187+
+" unfinished in "+std::string(__FILE__)+" line "+std::to_string(__LINE__));
188+
break;
189+
}
190+
} // end switch( func.info->family )
191+
} // end for( xc_func_type &func : funcs )
192+
193+
auto write_data = [&esolver](
194+
const std::string data_name,
195+
const std::vector<double> &data,
196+
const int number_spin)
197+
{
198+
for(int is=0; is<number_spin; ++is)
199+
{
200+
std::ofstream ofs;
201+
if(GlobalV::MY_RANK==0)
202+
{
203+
const std::string folder_name = PARAM.globalv.global_out_dir + "/xc_r/";
204+
const std::string command0 = "test -d " + folder_name + " || mkdir " + folder_name;
205+
system( command0.c_str() );
206+
const std::string file_name = folder_name + data_name+"_"+std::to_string(is);
207+
ofs.open(file_name);
208+
209+
ofs.unsetf(std::ostream::fixed);
210+
ofs << std::setprecision(PARAM.inp.out_xc_r[1]);
211+
ofs << std::scientific;
212+
}
213+
214+
#ifdef __MPI
215+
ModuleIO::write_cube_core(
216+
ofs,
217+
esolver->pw_big->bz,
218+
esolver->pw_big->nbz,
219+
esolver->pw_rhod->nplane * number_spin,
220+
esolver->pw_rhod->startz_current,
221+
data.data()+is,
222+
esolver->pw_rhod->nx * esolver->pw_rhod->ny,
223+
esolver->pw_rhod->nz,
224+
number_spin,
225+
esolver->pw_rhod->nz);
226+
#else
227+
if(nspin!=1)
228+
{ throw std::invalid_argument("nspin="+std::to_string(nspin)+" is invalid for ModuleIO::write_cube_core without MPI. see "+std:string(__FILE__)+" line "+std::to_string(__LINE__)); }
229+
ModuleIO::write_cube_core(
230+
ofs,
231+
data.data()+is,
232+
esolver->pw_rhod->nx * esolver->pw_rhod->ny,
233+
esolver->pw_rhod->nz,
234+
esolver->pw_rhod->nz);
235+
#endif
236+
}
237+
};
238+
239+
write_data( "rho", rho, nspin );
240+
241+
if(1!=nspin && 2!=PARAM.inp.nspin)
242+
write_data( "amag", amag, 1 );
243+
244+
if(is_gga)
245+
write_data( "sigma", sigma, (1==nspin)?1:3 );
246+
247+
switch( order )
248+
{
249+
case 4: write_data( "v4rho4", v4rho4, (1==nspin)?1:5 );
250+
case 3: write_data( "v3rho3", v3rho3, (1==nspin)?1:4 );
251+
case 2: write_data( "v2rho2", v2rho2, (1==nspin)?1:3 );
252+
case 1: write_data( "vrho" , vrho , nspin );
253+
case 0: write_data( "exc" , exc , 1 );
254+
break;
255+
default: throw std::domain_error("order ="+std::to_string(order)
256+
+" unfinished in "+std::string(__FILE__)+" line "+std::to_string(__LINE__));
257+
break;
258+
}
259+
if(is_gga)
260+
{
261+
switch( order )
262+
{
263+
case 4: write_data( "v4rho3sigma" , v4rho3sigma , (1==nspin)?1:12 );
264+
write_data( "v4rho2sigma2", v4rho2sigma2, (1==nspin)?1:15 );
265+
write_data( "v4rhosigma3" , v4rhosigma3 , (1==nspin)?1:20 );
266+
write_data( "v4sigma4" , v4sigma4 , (1==nspin)?1:15 );
267+
case 3: write_data( "v3rho2sigma" , v3rho2sigma , (1==nspin)?1:9 );
268+
write_data( "v3rhosigma2" , v3rhosigma2 , (1==nspin)?1:12 );
269+
write_data( "v3sigma3" , v3sigma3 , (1==nspin)?1:10 );
270+
case 2: write_data( "v2rhosigma" , v2rhosigma , (1==nspin)?1:6 );
271+
write_data( "v2sigma2" , v2sigma2 , (1==nspin)?1:6 );
272+
case 1: write_data( "vsigma" , vsigma , (1==nspin)?1:3 );
273+
case 0: break;
274+
default: throw std::domain_error("order ="+std::to_string(order)
275+
+" unfinished in "+std::string(__FILE__)+" line "+std::to_string(__LINE__));
276+
break;
277+
}
278+
}
279+
280+
XC_Functional_Libxc::finish_func(funcs);
281+
282+
ModuleBase::timer::tick("ModuleIO","write_libxc_r");
283+
}
284+
285+
#endif // USE_LIBXC

‎source/module_io/write_libxc_r.h‎

Lines changed: 30 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,30 @@
1+
//======================
2+
// AUTHOR : Peize Lin
3+
// DATE : 2024-09-12
4+
//======================
5+
6+
#ifndef WRITE_LIBXC_R_H
7+
#define WRITE_LIBXC_R_H
8+
9+
#ifdef USE_LIBXC
10+
11+
#include <vector>
12+
13+
class Charge;
14+
namespace ModuleESolver{ class ESolver_FP; }
15+
16+
namespace ModuleIO
17+
{
18+
extern void write_libxc_r(
19+
const int order,
20+
const std::vector<int> &func_id,
21+
const int &nrxx, // number of real-space grid
22+
const double &omega, // volume of cell
23+
const double tpiba,
24+
const Charge* const chr,
25+
const ModuleESolver::ESolver_FP*const esolver);
26+
}
27+
28+
#endif // USE_LIBXC
29+
30+
#endif // WRITE_LIBXC_R_H

‎source/module_lr/esolver_lrtd_lcao.cpp‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -87,7 +87,7 @@ template<typename T, typename TR>
8787
void LR::ESolver_LR<T, TR>::parameter_check()const
8888
{
8989
std::set<std::string> lr_solvers = { "dav", "lapack" , "spectrum", "dav_subspace" };
90-
std::set<std::string> xc_kernels = { "rpa", "lda", "pbe", "hf" , "hse" };
90+
std::set<std::string> xc_kernels = { "rpa", "lda", "pwlda", "pbe", "hf" , "hse" };
9191
if (lr_solvers.find(this->input.lr_solver) == lr_solvers.end()) {
9292
throw std::invalid_argument("ESolver_LR: unknown type of lr_solver");
9393
}
@@ -603,11 +603,11 @@ void LR::ESolver_LR<T, TR>::init_pot(const Charge& chg_gs)
603603
{
604604
using ST = PotHxcLR::SpinType;
605605
case 1:
606-
this->pot[0] = std::make_shared<PotHxcLR>(xc_kernel, this->pw_rho, &ucell, &chg_gs, GlobalC::Pgrid, ST::S1, input.init_fxc);
606+
this->pot[0] = std::make_shared<PotHxcLR>(xc_kernel, this->pw_rho, &ucell, &chg_gs, GlobalC::Pgrid, ST::S1, input.init_xc_kernel);
607607
break;
608608
case 2:
609-
this->pot[0] = std::make_shared<PotHxcLR>(xc_kernel, this->pw_rho, &ucell, &chg_gs, GlobalC::Pgrid, openshell ? ST::S2_updown : ST::S2_singlet, input.init_fxc);
610-
this->pot[1] = std::make_shared<PotHxcLR>(xc_kernel, this->pw_rho, &ucell, &chg_gs, GlobalC::Pgrid, openshell ? ST::S2_updown : ST::S2_triplet, input.init_fxc);
609+
this->pot[0] = std::make_shared<PotHxcLR>(xc_kernel, this->pw_rho, &ucell, &chg_gs, GlobalC::Pgrid, openshell ? ST::S2_updown : ST::S2_singlet, input.init_xc_kernel);
610+
this->pot[1] = std::make_shared<PotHxcLR>(xc_kernel, this->pw_rho, &ucell, &chg_gs, GlobalC::Pgrid, openshell ? ST::S2_updown : ST::S2_triplet, input.init_xc_kernel);
611611
break;
612612
default:
613613
throw std::invalid_argument("ESolver_LR: nspin must be 1 or 2");

0 commit comments

Comments
 (0)