Skip to content

Commit c87639c

Browse files
lyb9812lyb9812
committed
move memeber function of esolver_of_tddft to module_ofdft
1 parent 96dd8cc commit c87639c

File tree

6 files changed

+213
-294
lines changed

6 files changed

+213
-294
lines changed

‎source/Makefile.Objects‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -361,6 +361,7 @@ OBJS_HAMILT_OF=kedf_tf.o\
361361
kedf_xwm.o\
362362
kedf_lkt.o\
363363
kedf_manager.o\
364+
evolve_phi.o\
364365

365366
OBJS_HAMILT_LCAO=hamilt_lcao.o\
366367
operator_lcao.o\

‎source/source_esolver/esolver_of_tddft.cpp‎

Lines changed: 3 additions & 285 deletions
Original file line numberDiff line numberDiff line change
@@ -22,17 +22,12 @@ namespace ModuleESolver
2222
ESolver_OF_TDDFT::ESolver_OF_TDDFT()
2323
{
2424
this->classname = "ESolver_OF_TDDFT";
25-
/*this->pphi_td= new std::complex<double>*[PARAM.inp.nspin];
26-
27-
for (int is = 0; is < PARAM.inp.nspin; ++is)
28-
{
29-
this->pphi_td[is]= new std::complex<double>[pw_rho->nrxx];
30-
}*/
25+
this->evolve_psi=new EVOLVE_PSI();
3126
}
3227

3328
ESolver_OF_TDDFT::~ESolver_OF_TDDFT()
3429
{
35-
//delete psi_td;
30+
delete this->evolve_psi;
3631
for (int is = 0; is < PARAM.inp.nspin; ++is)
3732
{
3833
delete[] this->pphi_td[is];
@@ -107,7 +102,7 @@ void ESolver_OF_TDDFT::runner(UnitCell& ucell, const int istep)
107102
}
108103
else
109104
{
110-
this->propagate_psi(ucell,this->pphi_td,this->pw_rho);
105+
this->evolve_psi->propagate_psi(this->pelec, this->chr, ucell, this->pphi_td, this->pw_rho);
111106
for (int is = 0; is < PARAM.inp.nspin; ++is)
112107
{
113108
for (int ir = 0; ir < pw_rho->nrxx; ++ir)
@@ -123,281 +118,4 @@ void ESolver_OF_TDDFT::runner(UnitCell& ucell, const int istep)
123118
ModuleBase::timer::tick("ESolver_OF_TDDFT", "runner");
124119
}
125120

126-
/**
127-
* @brief Prepare to optimize the charge density,
128-
* update elecstate, kedf, and opts if needed
129-
* calculate ewald energy, initialize the rho, phi, theta
130-
*
131-
* @param istep
132-
* @param ucell
133-
*/
134-
void ESolver_OF_TDDFT::before_opt(const int istep, UnitCell& ucell)
135-
{
136-
ModuleBase::TITLE("ESolver_OF", "before_opt");
137-
ModuleBase::timer::tick("ESolver_OF", "before_opt");
138-
139-
//! 1) call before_scf() of ESolver_FP
140-
ESolver_FP::before_scf(ucell, istep);
141-
142-
if (ucell.cell_parameter_updated)
143-
{
144-
this->dV_ = ucell.omega / this->pw_rho->nxyz;
145-
146-
// initialize elecstate, including potential
147-
this->init_elecstate(ucell);
148-
149-
// Initialize KEDF
150-
this->kedf_manager_->init(PARAM.inp, this->pw_rho, this->dV_, this->nelec_[0]);
151-
152-
// Initialize optimization methods
153-
this->init_opt();
154-
155-
// Refresh the arrays
156-
delete this->psi_;
157-
delete this->psi_td;
158-
this->psi_ = new psi::Psi<double>(1,
159-
PARAM.inp.nspin,
160-
this->pw_rho->nrxx,
161-
this->pw_rho->nrxx,
162-
true);
163-
/*this->psi_td = new psi::Psi<std::complex<double>>(1,
164-
PARAM.inp.nspin,
165-
this->pw_rho->nrxx,
166-
this->pw_rho->nrxx,
167-
true);*/
168-
//this->pphi_td= new std::complex<double>*[PARAM.inp.nspin];
169-
170-
for (int is = 0; is < PARAM.inp.nspin; ++is)
171-
{
172-
this->pphi_[is] = this->psi_->get_pointer(is);
173-
//this->pphi_td[is] = this->psi_td->get_pointer(is);
174-
//this->pphi_td[is]= new std::complex<double>[pw_rho->nrxx];
175-
}
176-
177-
delete this->ptemp_rho_;
178-
this->ptemp_rho_ = new Charge();
179-
this->ptemp_rho_->set_rhopw(this->pw_rho);
180-
this->ptemp_rho_->allocate(PARAM.inp.nspin);
181-
182-
for (int is = 0; is < PARAM.inp.nspin; ++is)
183-
{
184-
delete[] this->pdLdphi_[is];
185-
delete[] this->pdEdphi_[is];
186-
delete[] this->pdirect_[is];
187-
delete[] this->precip_dir_[is];
188-
this->pdLdphi_[is] = new double[this->pw_rho->nrxx];
189-
this->pdEdphi_[is] = new double[this->pw_rho->nrxx];
190-
this->pdirect_[is] = new double[this->pw_rho->nrxx];
191-
this->precip_dir_[is] = new std::complex<double>[pw_rho->npw];
192-
}
193-
}
194-
195-
this->pelec->init_scf(istep, ucell, Pgrid, sf.strucFac, locpp.numeric, ucell.symm);
196-
197-
Symmetry_rho srho;
198-
for (int is = 0; is < PARAM.inp.nspin; is++)
199-
{
200-
srho.begin(is, this->chr, this->pw_rho, ucell.symm);
201-
}
202-
203-
for (int is = 0; is < PARAM.inp.nspin; ++is)
204-
{
205-
if (PARAM.inp.init_chg != "file")
206-
{
207-
for (int ibs = 0; ibs < this->pw_rho->nrxx; ++ibs)
208-
{
209-
// Here we initialize rho to be uniform,
210-
// because the rho got by pot.init_pot -> Charge::atomic_rho may contain minus elements.
211-
this->chr.rho[is][ibs] = this->nelec_[is] / this->pelec->omega;
212-
this->pphi_[is][ibs] = sqrt(this->chr.rho[is][ibs]);
213-
}
214-
}
215-
else
216-
{
217-
for (int ibs = 0; ibs < this->pw_rho->nrxx; ++ibs)
218-
{
219-
this->pphi_[is][ibs] = sqrt(this->chr.rho[is][ibs]);
220-
}
221-
}
222-
}
223-
224-
for (int is = 0; is < PARAM.inp.nspin; ++is)
225-
{
226-
this->pelec->eferm.set_efval(is, 0);
227-
this->theta_[is] = 0.;
228-
ModuleBase::GlobalFunc::ZEROS(this->pdLdphi_[is], this->pw_rho->nrxx);
229-
ModuleBase::GlobalFunc::ZEROS(this->pdEdphi_[is], this->pw_rho->nrxx);
230-
ModuleBase::GlobalFunc::ZEROS(this->pdirect_[is], this->pw_rho->nrxx);
231-
}
232-
if (PARAM.inp.nspin == 1)
233-
{
234-
this->theta_[0] = 0.2;
235-
}
236-
237-
ModuleBase::timer::tick("ESolver_OF", "before_opt");
238-
}
239-
240-
void ESolver_OF_TDDFT::get_Hpsi(UnitCell& ucell, const std::complex<double>* const* psi_, ModulePW::PW_Basis* pw_rho, std::complex<double>** Hpsi)
241-
{
242-
// update rho
243-
for (int is = 0; is < PARAM.inp.nspin; ++is)
244-
{
245-
for (int ir = 0; ir < pw_rho->nrxx; ++ir)
246-
{
247-
this->chr.rho[is][ir] = abs(psi_[is][ir])*abs(psi_[is][ir]);
248-
}
249-
}
250-
251-
this->pelec->pot->update_from_charge(&this->chr, &ucell); // Hartree + XC + external
252-
this->get_tf_potential(this->chr.rho,pw_rho ,this->pelec->pot->get_effective_v()); // TF potential
253-
for (int is = 0; is < PARAM.inp.nspin; ++is)
254-
{
255-
const double* vr_eff = this->pelec->pot->get_effective_v(is);
256-
for (int ir = 0; ir < pw_rho->nrxx; ++ir)
257-
{
258-
Hpsi[is][ir] = vr_eff[ir]*psi_[is][ir];
259-
}
260-
}
261-
this->get_vw_potential_phi(psi_, pw_rho, Hpsi);
262-
}
263-
264-
void ESolver_OF_TDDFT::get_tf_potential(const double* const* prho, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpot)
265-
{
266-
const double c_tf_
267-
= 3.0 / 10.0 * std::pow(3 * std::pow(M_PI, 2.0), 2.0 / 3.0)
268-
* 2;
269-
if (PARAM.inp.nspin == 1)
270-
{
271-
for (int ir = 0; ir < pw_rho->nrxx; ++ir)
272-
{
273-
rpot(0, ir) += 5.0 / 3.0 * c_tf_ * std::pow(prho[0][ir], 2. / 3.);
274-
}
275-
}
276-
else if (PARAM.inp.nspin == 2)
277-
{
278-
for (int is = 0; is < PARAM.inp.nspin; ++is)
279-
{
280-
for (int ir = 0; ir < pw_rho->nrxx; ++ir)
281-
{
282-
rpot(is, ir) += 5.0 / 3.0 * c_tf_ * std::pow(2. * prho[is][ir], 2. / 3.);
283-
}
284-
}
285-
}
286-
}
287-
288-
void ESolver_OF_TDDFT::get_vw_potential_phi(const std::complex<double>* const* pphi, ModulePW::PW_Basis* pw_rho, std::complex<double>** Hpsi)
289-
{
290-
std::complex<double>** rLapPhi = new std::complex<double>*[PARAM.inp.nspin];
291-
for (int is = 0; is < PARAM.inp.nspin; ++is) {
292-
rLapPhi[is] = new std::complex<double>[pw_rho->nrxx];
293-
}
294-
std::complex<double>** recipPhi = new std::complex<double>*[PARAM.inp.nspin];
295-
for (int is = 0; is < PARAM.inp.nspin; ++is)
296-
{
297-
recipPhi[is] = new std::complex<double>[pw_rho->npw];
298-
299-
pw_rho->real2recip(pphi[is], recipPhi[is]);
300-
for (int ik = 0; ik < pw_rho->npw; ++ik)
301-
{
302-
recipPhi[is][ik] *= pw_rho->gg[ik] * pw_rho->tpiba2;
303-
}
304-
pw_rho->recip2real(recipPhi[is], rLapPhi[is]);
305-
for (int ik = 0; ik < pw_rho->npw; ++ik)
306-
{
307-
Hpsi[is][ik]+=rLapPhi[is][ik];
308-
}
309-
}
310-
311-
for (int is = 0; is < PARAM.inp.nspin; ++is)
312-
{
313-
delete[] recipPhi[is];
314-
delete[] rLapPhi[is];
315-
}
316-
delete[] recipPhi;
317-
delete[] rLapPhi;
318-
}
319-
320-
void ESolver_OF_TDDFT::get_CD_potential(const std::complex<double>* const* psi_, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpot)
321-
{
322-
for (int is = 0; is < PARAM.inp.nspin; ++is)
323-
{
324-
//recipCurrent = new std::complex<double>[pw_rho->npw];
325-
//delete[] recipCurrent;
326-
}
327-
}
328-
329-
void ESolver_OF_TDDFT::propagate_psi(UnitCell& ucell, std::complex<double>** pphi_, ModulePW::PW_Basis* pw_rho)
330-
{
331-
ModuleBase::timer::tick("ESolver_OF_TDDFT", "propagte_psi");
332-
333-
std::complex<double> imag(0.0,1.0);
334-
double dt=PARAM.inp.mdp.md_dt;
335-
std::complex<double>** K1 = new std::complex<double>*[PARAM.inp.nspin];
336-
std::complex<double>** K2 = new std::complex<double>*[PARAM.inp.nspin];
337-
std::complex<double>** K3 = new std::complex<double>*[PARAM.inp.nspin];
338-
std::complex<double>** K4 = new std::complex<double>*[PARAM.inp.nspin];
339-
std::complex<double>** psi1 = new std::complex<double>*[PARAM.inp.nspin];
340-
std::complex<double>** psi2 = new std::complex<double>*[PARAM.inp.nspin];
341-
std::complex<double>** psi3 = new std::complex<double>*[PARAM.inp.nspin];
342-
for (int is = 0; is < PARAM.inp.nspin; ++is) {
343-
K1[is] = new std::complex<double>[pw_rho->nrxx];
344-
K2[is] = new std::complex<double>[pw_rho->nrxx];
345-
K3[is] = new std::complex<double>[pw_rho->nrxx];
346-
K4[is] = new std::complex<double>[pw_rho->nrxx];
347-
psi1[is] = new std::complex<double>[pw_rho->nrxx];
348-
psi2[is] = new std::complex<double>[pw_rho->nrxx];
349-
psi3[is] = new std::complex<double>[pw_rho->nrxx];
350-
}
351-
get_Hpsi(ucell,pphi_,pw_rho,K1);
352-
for (int is = 0; is < PARAM.inp.nspin; ++is){
353-
for (int ir = 0; ir < pw_rho->nrxx; ++ir)
354-
{
355-
K1[is][ir]=-1.0*K1[is][ir]*dt*imag;
356-
psi1[is][ir]=pphi_[is][ir]+0.5*K1[is][ir];
357-
}
358-
}
359-
get_Hpsi(ucell,psi1,pw_rho,K2);
360-
for (int is = 0; is < PARAM.inp.nspin; ++is){
361-
for (int ir = 0; ir < pw_rho->nrxx; ++ir)
362-
{
363-
K2[is][ir]=-1.0*K2[is][ir]*dt*imag;
364-
psi2[is][ir]=pphi_[is][ir]+0.5*K2[is][ir];
365-
}
366-
}
367-
get_Hpsi(ucell,psi2,pw_rho,K3);
368-
for (int is = 0; is < PARAM.inp.nspin; ++is){
369-
for (int ir = 0; ir < pw_rho->nrxx; ++ir)
370-
{
371-
K3[is][ir]=-1.0*K3[is][ir]*dt*imag;
372-
psi3[is][ir]=pphi_[is][ir]+K3[is][ir];
373-
}
374-
}
375-
get_Hpsi(ucell,psi3,pw_rho,K4);
376-
for (int is = 0; is < PARAM.inp.nspin; ++is){
377-
for (int ir = 0; ir < pw_rho->nrxx; ++ir)
378-
{
379-
K4[is][ir]=-1.0*K4[is][ir]*dt*imag;
380-
pphi_[is][ir]+=1.0/6.0*(K1[is][ir]+2.0*K2[is][ir]+2.0*K3[is][ir]+K4[is][ir]);
381-
}
382-
}
383-
384-
for (int is = 0; is < PARAM.inp.nspin; ++is) {
385-
delete[] K1[is];
386-
delete[] K2[is];
387-
delete[] K3[is];
388-
delete[] K4[is];
389-
delete[] psi1[is];
390-
delete[] psi2[is];
391-
delete[] psi3[is];
392-
}
393-
delete[] K1;
394-
delete[] K2;
395-
delete[] K3;
396-
delete[] K4;
397-
delete[] psi1;
398-
delete[] psi2;
399-
delete[] psi3;
400-
ModuleBase::timer::tick("ESolver_OF_TDDFT", "propagte_psi");
401-
}
402-
403121
} // namespace ModuleESolver

‎source/source_esolver/esolver_of_tddft.h‎

Lines changed: 2 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,7 @@
22
#define ESOLVER_OF_TDDFT_H
33

44
#include "esolver_of.h"
5+
#include "source_pw/module_ofdft/evolve_psi.h"
56

67
namespace ModuleESolver
78
{
@@ -15,15 +16,7 @@ class ESolver_OF_TDDFT : public ESolver_OF
1516

1617
protected:
1718
std::complex<double>** pphi_td = nullptr; // pphi[i] = ppsi.get_pointer(i), which will be freed in ~Psi().
18-
psi::Psi<std::complex<double>>* psi_td = nullptr;
19-
//std::complex<double>** pdEdphi_phi_ = nullptr; // (dE/dphi)*phi
20-
// ==================== main process of TDOFDFT ======================
21-
void before_opt(const int istep, UnitCell& ucell);
22-
void get_Hpsi(UnitCell& ucell, const std::complex<double>* const* psi_, ModulePW::PW_Basis* pw_rho, std::complex<double>** Hpsi);
23-
void get_tf_potential(const double* const* prho, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpot);
24-
void get_vw_potential_phi(const std::complex<double>* const* pphi, ModulePW::PW_Basis* pw_rho, std::complex<double>** Hpsi); // -1/2 \nabla^2 \phi
25-
void get_CD_potential(const std::complex<double>* const* psi_, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpot);
26-
void propagate_psi(UnitCell& ucell, std::complex<double>** pphi_, ModulePW::PW_Basis* pw_rho);
19+
EVOLVE_PSI* evolve_psi=nullptr;
2720
};
2821
} // namespace ModuleESolver
2922

‎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_phi.cpp
910
)
1011

1112
add_library(

0 commit comments

Comments
 (0)