Skip to content

Commit 4c72902

Browse files
lyb9812sunliang98
lyb9812
and
sunliang98
authored
[Feature] Add current dependent potential for TD-OFDFT (#6593)
* add esolver_of_tddft * esolver_of_tddft cmakelist * move memeber function of esolver_of_tddft to module_ofdft * fix bug (confuse evolve_psi with evolve_phi) * change evolve_phi with evolve_ofdft; change to vector * details * code optimization * add situation of reading charge file for tdofdft * add CD Potential * add INPUT parameters for tdofdft * add integrate test for TDOFDFT * update doc file for tdofdft * solve conflict * type transformation * type transformation * move integration test for TDOFDFT to 07_OFDFT * add judgement for nspin in TDOFDFT * add judgement for nspin in TDOFDFT * TDOFDFT integration test * change type of warning --------- Co-authored-by: Liang Sun <[email protected]>
1 parent 487555c commit 4c72902

File tree

25 files changed

+350
-6
lines changed

25 files changed

+350
-6
lines changed

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

Lines changed: 23 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -256,6 +256,9 @@
256256
- [of\_ml\_chi\_pnl](#of_ml_chi_pnl)
257257
- [of\_ml\_chi\_qnl](#of_ml_chi_qnl)
258258
- [of\_ml\_local\_test](#of_ml_local_test)
259+
- [TD-OFDFT: time dependent orbital free density functional theory](#tdofdft-time-dependent-orbital-free-density-functional-theory)
260+
- [of\_cd](#of_cd)
261+
- [of\_mcd\_alpha](#of_mcd_alpha)
259262
- [Electric Field and Dipole Correction](#electric-field-and-dipole-correction)
260263
- [efield\_flag](#efield_flag)
261264
- [dip\_cor\_flag](#dip_cor_flag)
@@ -528,6 +531,7 @@ These variables are used to control general system parameters.
528531
- **Description**: choose the energy solver.
529532
- ksdft: Kohn-Sham density functional theory
530533
- ofdft: orbital-free density functional theory
534+
- tdofdft: time-dependent orbital-free density functional theory
531535
- sdft: [stochastic density functional theory](#electronic-structure-sdft)
532536
- tddft: real-time time-dependent density functional theory (TDDFT)
533537
- lj: Leonard Jones potential
@@ -2747,6 +2751,25 @@ Warning: this function is not robust enough for the current version. Please try
27472751

27482752
[back to top](#full-list-of-input-keywords)
27492753

2754+
## TDOFDFT: time dependent orbital free density functional theory
2755+
2756+
### of_cd
2757+
2758+
- **Type**: Boolean
2759+
- **Availability**: TDOFDFT
2760+
- **Type**: Boolean
2761+
- **Description**: Added the current dependent(CD) potential. (https://doi.org/10.1103/PhysRevB.98.144302)
2762+
- True: Added the CD potential.
2763+
- False: Not added the CD potential.
2764+
- **Default**: False
2765+
2766+
### of_mcd_alpha
2767+
2768+
- **Type**: Real
2769+
- **Availability**: TDOFDFT
2770+
- **Description**: The value of the parameter alpha in modified CD potential method. mCDPotenial=alpha*CDPotenial(proposed in paper PhysRevB.98.144302)
2771+
- **Default**: 1.0
2772+
27502773
## Electric field and dipole correction
27512774

27522775
These variables are relevant to electric field and dipole correction

‎source/source_io/module_parameter/input_parameter.h‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -687,5 +687,10 @@ struct Input_para
687687
// src/gga_c_pbe.c
688688
std::vector<double> xc_corr_ext = {
689689
130, 0.06672455060314922, 0.031090690869654895034, 1.00000};
690+
691+
// ============== #Parameters (24.td-ofdft) ===========================
692+
bool of_cd = false; ///< add CD potential or not https://doi.org/10.1103/PhysRevB.98.144302
693+
double of_mCD_alpha = 1.0; /// parameter of modified CD Potential
694+
690695
};
691696
#endif

‎source/source_io/read_input.cpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -101,6 +101,7 @@ ReadInput::ReadInput(const int& rank)
101101
this->item_sdft();
102102
this->item_deepks();
103103
this->item_rt_tddft();
104+
this->item_tdofdft();
104105
this->item_lr_tddft();
105106
this->item_output();
106107
this->item_postprocess();

‎source/source_io/read_input.h‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -108,6 +108,8 @@ class ReadInput
108108
// items for real time tddft
109109
void item_rt_tddft();
110110
// items for linear response tddft
111+
void item_tdofdft();
112+
// items for td-ofdft
111113
void item_lr_tddft();
112114
// items for output
113115
void item_output();

‎source/source_io/read_input_item_tddft.cpp‎

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -305,6 +305,22 @@ void ReadInput::item_rt_tddft()
305305
}
306306

307307

308+
}
309+
void ReadInput::item_tdofdft()
310+
{
311+
// TD-OFDFT
312+
{
313+
Input_Item item("of_cd");
314+
item.annotation = "add CD Potential or not";
315+
read_sync_bool(input.of_cd);
316+
this->add_item(item);
317+
}
318+
{
319+
Input_Item item("of_mcd_alpha");
320+
item.annotation = "parameter of modified CD Potential";
321+
read_sync_double(input.of_mCD_alpha);
322+
this->add_item(item);
323+
}
308324
}
309325
void ReadInput::item_lr_tddft()
310326
{

‎source/source_io/test/read_input_ptest.cpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -364,6 +364,8 @@ TEST_F(InputParaTest, ParaRead)
364364
EXPECT_EQ(param.inp.of_full_pw_dim, 0);
365365
EXPECT_FALSE(param.inp.of_read_kernel);
366366
EXPECT_EQ(param.inp.of_kernel_file, "WTkernel.txt");
367+
EXPECT_FALSE(param.inp.of_cd);
368+
EXPECT_DOUBLE_EQ(param.inp.of_mCD_alpha,1.0);
367369
EXPECT_DOUBLE_EQ(param.inp.of_xwm_kappa, 1.);
368370
EXPECT_DOUBLE_EQ(param.inp.of_xwm_rho_ref, 1.);
369371
EXPECT_EQ(param.inp.device, "cpu");

‎source/source_io/test/support/INPUT‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -387,3 +387,7 @@ nsc_min 4 #Minimum number of spin-constrained iteration
387387
sc_scf_nmin 4 #Minimum number of outer scf loop before initializing lambda loop
388388
alpha_trial 0.02 #Initial trial step size for lambda in eV/uB^2
389389
sccut 4 #Maximal step size for lambda in eV/uB
390+
391+
#Parameters (23. Time-dependent orbital-free DFT)
392+
of_cd 0 #0: no CD potential; 1: add CD potential
393+
of_mCD_alpha 1.0 # parameter of modified CD potential

‎source/source_pw/module_ofdft/evolve_ofdft.cpp‎

Lines changed: 94 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -25,7 +25,11 @@ void Evolve_OFDFT::cal_Hpsi(elecstate::ElecState* pelec,
2525
}
2626

2727
pelec->pot->update_from_charge(&chr, &ucell); // Hartree + XC + external
28-
this->cal_tf_potential(chr.rho,pw_rho ,pelec->pot->get_effective_v()); // TF potential
28+
this->cal_tf_potential(chr.rho, pw_rho, pelec->pot->get_effective_v()); // TF potential
29+
if (PARAM.inp.of_cd)
30+
{
31+
this->cal_CD_potential(psi_, pw_rho, pelec->pot->get_effective_v(), PARAM.inp.of_mCD_alpha); // CD potential
32+
}
2933

3034
#ifdef _OPENMP
3135
#pragma omp parallel for
@@ -72,6 +76,9 @@ void Evolve_OFDFT::cal_vw_potential_phi(std::vector<std::complex<double>> pphi,
7276
ModulePW::PW_Basis* pw_rho,
7377
std::vector<std::complex<double>> Hpsi)
7478
{
79+
if (PARAM.inp.nspin <= 0) {
80+
ModuleBase::WARNING_QUIT("Evolve_OFDFT","nspin must be positive");
81+
}
7582
std::complex<double>** rLapPhi = new std::complex<double>*[PARAM.inp.nspin];
7683
#ifdef _OPENMP
7784
#pragma omp parallel for
@@ -118,13 +125,96 @@ void Evolve_OFDFT::cal_vw_potential_phi(std::vector<std::complex<double>> pphi,
118125

119126
void Evolve_OFDFT::cal_CD_potential(std::vector<std::complex<double>> psi_,
120127
ModulePW::PW_Basis* pw_rho,
121-
ModuleBase::matrix& rpot)
128+
ModuleBase::matrix& rpot,
129+
double mCD_para)
122130
{
131+
std::complex<double> imag(0.0,1.0);
132+
133+
if (PARAM.inp.nspin <= 0) {
134+
ModuleBase::WARNING_QUIT("Evolve_OFDFT","nspin must be positive");
135+
}
136+
std::complex<double>** recipPhi = new std::complex<double>*[PARAM.inp.nspin];
137+
std::complex<double>** rPhi = new std::complex<double>*[PARAM.inp.nspin];
138+
#ifdef _OPENMP
139+
#pragma omp parallel for
140+
#endif
141+
for (int is = 0; is < PARAM.inp.nspin; ++is) {
142+
rPhi[is] = new std::complex<double>[pw_rho->nrxx];
143+
for (int ir = 0; ir < pw_rho->nrxx; ++ir)
144+
{
145+
rPhi[is][ir]=psi_[is * pw_rho->nrxx + ir];
146+
}
147+
}
148+
149+
#ifdef _OPENMP
150+
#pragma omp parallel for
151+
#endif
123152
for (int is = 0; is < PARAM.inp.nspin; ++is)
124153
{
125-
//recipCurrent = new std::complex<double>[pw_rho->npw];
126-
//delete[] recipCurrent;
154+
std::complex<double> *recipCurrent_x=new std::complex<double>[pw_rho->npw];
155+
std::complex<double> *recipCurrent_y=new std::complex<double>[pw_rho->npw];
156+
std::complex<double> *recipCurrent_z=new std::complex<double>[pw_rho->npw];
157+
std::complex<double> *recipCDPotential=new std::complex<double>[pw_rho->npw];
158+
std::complex<double> *rCurrent_x=new std::complex<double>[pw_rho->nrxx];
159+
std::complex<double> *rCurrent_y=new std::complex<double>[pw_rho->nrxx];
160+
std::complex<double> *rCurrent_z=new std::complex<double>[pw_rho->nrxx];
161+
std::complex<double> *kF_r=new std::complex<double>[pw_rho->nrxx];
162+
std::complex<double> *rCDPotential=new std::complex<double>[pw_rho->nrxx];
163+
recipPhi[is] = new std::complex<double>[pw_rho->npw];
164+
165+
for (int ir = 0; ir < pw_rho->nrxx; ++ir)
166+
{
167+
kF_r[ir]=std::pow(3*std::pow(ModuleBase::PI*std::abs(rPhi[is][ir]),2),1/3);
168+
}
169+
170+
pw_rho->real2recip(rPhi[is], recipPhi[is]);
171+
for (int ik = 0; ik < pw_rho->npw; ++ik)
172+
{
173+
recipCurrent_x[ik]=imag*pw_rho->gcar[ik].x*recipPhi[is][ik]* pw_rho->tpiba;
174+
recipCurrent_y[ik]=imag*pw_rho->gcar[ik].y*recipPhi[is][ik]* pw_rho->tpiba;
175+
recipCurrent_z[ik]=imag*pw_rho->gcar[ik].z*recipPhi[is][ik]* pw_rho->tpiba;
176+
}
177+
pw_rho->recip2real(recipCurrent_x,rCurrent_x);
178+
pw_rho->recip2real(recipCurrent_y,rCurrent_y);
179+
pw_rho->recip2real(recipCurrent_z,rCurrent_z);
180+
for (int ir = 0; ir < pw_rho->nrxx; ++ir)
181+
{
182+
rCurrent_x[ir]=std::imag(rCurrent_x[ir]*std::conj(rPhi[is][ir]));
183+
rCurrent_y[ir]=std::imag(rCurrent_y[ir]*std::conj(rPhi[is][ir]));
184+
rCurrent_z[ir]=std::imag(rCurrent_z[ir]*std::conj(rPhi[is][ir]));
185+
}
186+
pw_rho->real2recip(rCurrent_x,recipCurrent_x);
187+
pw_rho->real2recip(rCurrent_y,recipCurrent_y);
188+
pw_rho->real2recip(rCurrent_z,recipCurrent_z);
189+
for (int ik = 0; ik < pw_rho->npw; ++ik)
190+
{
191+
recipCDPotential[ik]=recipCurrent_x[ik]*pw_rho->gcar[ik].x+recipCurrent_y[ik]*pw_rho->gcar[ik].y+recipCurrent_z[ik]*pw_rho->gcar[ik].z;
192+
recipCDPotential[ik]*=imag/pw_rho->gg[ik];
193+
}
194+
pw_rho->recip2real(recipCDPotential,rCDPotential);
195+
196+
for (int ir = 0; ir < pw_rho->nrxx; ++ir)
197+
{
198+
rpot(0, ir) -= mCD_para*2.0*std::real(rCDPotential[ir])*std::pow(ModuleBase::PI,3) / (2.0*std::pow(std::real(kF_r[ir]),2));
199+
}
200+
delete[] recipCurrent_x;
201+
delete[] recipCurrent_y;
202+
delete[] recipCurrent_z;
203+
delete[] rCurrent_x;
204+
delete[] rCurrent_y;
205+
delete[] rCurrent_z;
127206
}
207+
208+
#ifdef _OPENMP
209+
#pragma omp parallel for
210+
#endif
211+
for (int is = 0; is < PARAM.inp.nspin; ++is)
212+
{
213+
delete[] recipPhi[is];
214+
delete[] rPhi[is];
215+
}
216+
delete[] recipPhi;
217+
delete[] rPhi;
128218
}
129219

130220
void Evolve_OFDFT::propagate_psi(elecstate::ElecState* pelec,

‎source/source_pw/module_ofdft/evolve_ofdft.h‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -48,7 +48,8 @@ class Evolve_OFDFT
4848
std::vector<std::complex<double>> Hpsi); // -1/2 \nabla^2 \phi
4949
void cal_CD_potential(std::vector<std::complex<double>> psi_,
5050
ModulePW::PW_Basis* pw_rho,
51-
ModuleBase::matrix& rpot);
51+
ModuleBase::matrix& rpot,
52+
double mCD_para);
5253

5354
};
5455
#endif

‎tests/07_OFDFT/30_TDOFDFT_Al/INPUT‎

Lines changed: 36 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,36 @@
1+
INPUT_PARAMETERS
2+
#Parameters (1.General)
3+
suffix autotest
4+
calculation md
5+
esolver_type tdofdft
6+
7+
pseudo_dir ../../PP_ORB
8+
pseudo_rcut 16
9+
cal_force 1
10+
cal_stress 1
11+
12+
#Parameters (2.Iteration)
13+
ecutwfc 17
14+
scf_nmax 100
15+
16+
#OFDFT
17+
of_kinetic tf+
18+
of_method tn
19+
of_full_pw 1
20+
of_full_pw_dim 1
21+
of_cd 0
22+
23+
#Parameters (3.Basis)
24+
basis_type pw
25+
26+
init_vel 1
27+
28+
md_restart 0
29+
md_type nve
30+
md_nstep 2
31+
md_dt 0.25
32+
md_tfirst 58022.52706
33+
md_dumpfreq 10
34+
md_tfreq 1.08
35+
md_tchain 1
36+
nbspline 10 # be sure fft dimension is odd

0 commit comments

Comments
 (0)