Skip to content

Commit 1c54579

Browse files
lyb9812sunliang98
lyb9812
and
sunliang98
authored
[Feature] TD-OFDFT (#6538)
* 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 --------- Co-authored-by: Liang Sun <[email protected]>
1 parent 73fc8a9 commit 1c54579

File tree

10 files changed

+419
-4
lines changed

10 files changed

+419
-4
lines changed

‎source/Makefile.Objects‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -266,6 +266,7 @@ OBJS_ESOLVER=esolver.o\
266266
esolver_lj.o\
267267
esolver_dp.o\
268268
esolver_of.o\
269+
esolver_of_tddft.o\
269270
esolver_of_tool.o\
270271
esolver_of_interface.o\
271272
pw_others.o\
@@ -360,6 +361,7 @@ OBJS_HAMILT_OF=kedf_tf.o\
360361
kedf_xwm.o\
361362
kedf_lkt.o\
362363
kedf_manager.o\
364+
evolve_ofdft.o\
363365

364366
OBJS_HAMILT_LCAO=hamilt_lcao.o\
365367
operator_lcao.o\

‎source/source_esolver/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@ list(APPEND objects
88
esolver_lj.cpp
99
esolver_dp.cpp
1010
esolver_of.cpp
11+
esolver_of_tddft.cpp
1112
esolver_of_interface.cpp
1213
esolver_of_tool.cpp
1314
pw_others.cpp

‎source/source_esolver/esolver.cpp‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -20,6 +20,7 @@ extern "C"
2020
#include "esolver_dp.h"
2121
#include "esolver_lj.h"
2222
#include "esolver_of.h"
23+
#include "esolver_of_tddft.h"
2324
#include "source_io/module_parameter/md_parameter.h"
2425

2526
#include <stdexcept>
@@ -40,6 +41,10 @@ std::string determine_type()
4041
{
4142
esolver_type = "ofdft";
4243
}
44+
else if (PARAM.inp.esolver_type == "tdofdft")
45+
{
46+
esolver_type = "tdofdft";
47+
}
4348
else if (PARAM.inp.esolver_type == "ksdft")
4449
{
4550
esolver_type = "ksdft_pw";
@@ -321,6 +326,10 @@ ESolver* init_esolver(const Input_para& inp, UnitCell& ucell)
321326
{
322327
return new ESolver_OF();
323328
}
329+
else if (esolver_type == "tdofdft")
330+
{
331+
return new ESolver_OF_TDDFT();
332+
}
324333
else if (esolver_type == "lj_pot")
325334
{
326335
return new ESolver_LJ();

‎source/source_esolver/esolver_of.h‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -27,7 +27,7 @@ class ESolver_OF : public ESolver_FP
2727

2828
virtual void cal_stress(UnitCell& ucell, ModuleBase::matrix& stress) override;
2929

30-
private:
30+
protected:
3131
// ======================= variables ==========================
3232
// ---------- the kinetic energy density functionals ----------
3333
KEDF_Manager* kedf_manager_ = nullptr; // KEDF manager, which will be initialized in before_all_runners
@@ -83,7 +83,7 @@ class ESolver_OF : public ESolver_FP
8383

8484
// ============================ tools ===============================
8585
// --------------------- initialize ---------------------------------
86-
void init_elecstate(UnitCell& ucell);
86+
void init_elecstate(UnitCell& ucell);
8787
void allocate_array();
8888

8989
// --------------------- calculate physical qualities ---------------

‎source/source_esolver/esolver_of_tddft.cpp‎

Lines changed: 129 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,129 @@
1+
#include "esolver_of_tddft.h"
2+
3+
#include "source_io/module_parameter/parameter.h"
4+
#include "source_io/cube_io.h"
5+
#include "source_io/output_log.h"
6+
#include "source_io/write_elecstat_pot.h"
7+
//-----------temporary-------------------------
8+
#include "source_base/global_function.h"
9+
#include "source_estate/module_charge/symmetry_rho.h"
10+
#include "source_hamilt/module_ewald/H_Ewald_pw.h"
11+
#include "source_pw/module_pwdft/global.h"
12+
#include "source_io/print_info.h"
13+
#include "source_estate/cal_ux.h"
14+
//-----force-------------------
15+
#include "source_pw/module_pwdft/forces.h"
16+
//-----stress------------------
17+
#include "source_pw/module_ofdft/of_stress_pw.h"
18+
19+
namespace ModuleESolver
20+
{
21+
22+
ESolver_OF_TDDFT::ESolver_OF_TDDFT()
23+
{
24+
this->classname = "ESolver_OF_TDDFT";
25+
this->evolve_ofdft=new Evolve_OFDFT();
26+
}
27+
28+
ESolver_OF_TDDFT::~ESolver_OF_TDDFT()
29+
{
30+
delete this->evolve_ofdft;
31+
}
32+
33+
34+
void ESolver_OF_TDDFT::runner(UnitCell& ucell, const int istep)
35+
{
36+
ModuleBase::timer::tick("ESolver_OF_TDDFT", "runner");
37+
// get Ewald energy, initial rho and phi if necessary
38+
this->before_opt(istep, ucell);
39+
this->iter_ = 0;
40+
41+
bool conv_esolver = false; // this conv_esolver is added by mohan 20250302
42+
#ifdef __MPI
43+
this->iter_time = MPI_Wtime();
44+
#else
45+
this->iter_time = std::chrono::system_clock::now();
46+
#endif
47+
48+
if (istep==0)
49+
{
50+
this->phi_td.resize(PARAM.inp.nspin*this->pw_rho->nrxx);
51+
}
52+
53+
if ((istep<1) && PARAM.inp.init_chg != "file")
54+
{
55+
while (true)
56+
{
57+
// once we get a new rho and phi, update potential
58+
this->update_potential(ucell);
59+
60+
// calculate the energy of new rho and phi
61+
this->energy_llast_ = this->energy_last_;
62+
this->energy_last_ = this->energy_current_;
63+
this->energy_current_ = this->cal_energy();
64+
65+
66+
// check if the job is done
67+
if (this->check_exit(conv_esolver))
68+
{
69+
break;
70+
}
71+
72+
// find the optimization direction and step lenghth theta according to the potential
73+
this->optimize(ucell);
74+
75+
// update the rho and phi based on the direction and theta
76+
this->update_rho();
77+
78+
this->iter_++;
79+
80+
ESolver_FP::iter_finish(ucell, istep, this->iter_, conv_esolver);
81+
}
82+
83+
#ifdef _OPENMP
84+
#pragma omp parallel for collapse(2)
85+
#endif
86+
for (int is = 0; is < PARAM.inp.nspin; ++is)
87+
{
88+
for (int ir = 0; ir < this->pw_rho->nrxx; ++ir)
89+
{
90+
phi_td[is*this->pw_rho->nrxx+ir]=pphi_[is][ir];
91+
}
92+
}
93+
}
94+
else if ((istep<1) && PARAM.inp.init_chg == "file")
95+
{
96+
#ifdef _OPENMP
97+
#pragma omp parallel for collapse(2)
98+
#endif
99+
for (int is = 0; is < PARAM.inp.nspin; ++is)
100+
{
101+
for (int ir = 0; ir < this->pw_rho->nrxx; ++ir)
102+
{
103+
phi_td[is*this->pw_rho->nrxx+ir]=pphi_[is][ir];
104+
}
105+
}
106+
conv_esolver=true;
107+
}
108+
else
109+
{
110+
this->evolve_ofdft->propagate_psi(this->pelec, this->chr, ucell, this->phi_td, this->pw_rho);
111+
#ifdef _OPENMP
112+
#pragma omp parallel for collapse(2)
113+
#endif
114+
for (int is = 0; is < PARAM.inp.nspin; ++is)
115+
{
116+
for (int ir = 0; ir < this->pw_rho->nrxx; ++ir)
117+
{
118+
pphi_[is][ir]=std::abs(phi_td[is*this->pw_rho->nrxx+ir]);
119+
}
120+
}
121+
conv_esolver=true;
122+
}
123+
124+
this->after_opt(istep, ucell, conv_esolver);
125+
126+
ModuleBase::timer::tick("ESolver_OF_TDDFT", "runner");
127+
}
128+
129+
} // namespace ModuleESolver

‎source/source_esolver/esolver_of_tddft.h‎

Lines changed: 23 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,23 @@
1+
#ifndef ESOLVER_OF_TDDFT_H
2+
#define ESOLVER_OF_TDDFT_H
3+
4+
#include "esolver_of.h"
5+
#include "source_pw/module_ofdft/evolve_ofdft.h"
6+
7+
namespace ModuleESolver
8+
{
9+
class ESolver_OF_TDDFT : public ESolver_OF
10+
{
11+
public:
12+
ESolver_OF_TDDFT();
13+
~ESolver_OF_TDDFT();
14+
15+
virtual void runner(UnitCell& ucell, const int istep) override;
16+
17+
protected:
18+
std::vector<std::complex<double>> phi_td; // pphi[i] = ppsi.get_pointer(i), which will be freed in ~Psi().
19+
Evolve_OFDFT* evolve_ofdft=nullptr;
20+
};
21+
} // namespace ModuleESolver
22+
23+
#endif

‎source/source_io/read_input_item_system.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -108,10 +108,10 @@ void ReadInput::item_system()
108108
}
109109
{
110110
Input_Item item("esolver_type");
111-
item.annotation = "the energy solver: ksdft, sdft, ofdft, tddft, lj, dp, ks-lr, lr";
111+
item.annotation = "the energy solver: ksdft, sdft, ofdft, tdofdft, tddft, lj, dp, ks-lr, lr";
112112
read_sync_string(input.esolver_type);
113113
item.check_value = [](const Input_Item& item, const Parameter& para) {
114-
const std::vector<std::string> esolver_types = { "ksdft", "sdft", "ofdft", "tddft", "lj", "dp", "lr", "ks-lr" };
114+
const std::vector<std::string> esolver_types = { "ksdft", "sdft", "ofdft", "tdofdft", "tddft", "lj", "dp", "lr", "ks-lr" };
115115
if (std::find(esolver_types.begin(), esolver_types.end(), para.input.esolver_type) == esolver_types.end())
116116
{
117117
const std::string warningstr = nofound_str(esolver_types, "esolver_type");

‎source/source_pw/module_ofdft/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@ list(APPEND hamilt_ofdft_srcs
66
kedf_lkt.cpp
77
kedf_manager.cpp
88
of_stress_pw.cpp
9+
evolve_ofdft.cpp
910
)
1011

1112
add_library(

0 commit comments

Comments
 (0)