Refactor: tddft code

source/Makefile.Objects

@@ -382,14 +382,19 @@ OBJS_IO_LCAO=cal_r_overlap_R.o\
 
 OBJS_LCAO=DM_gamma.o\
       DM_k.o\
-      ELEC_evolve.o\
+      evolve_elec.o\
+      evolve_psi.o\
+      bandenergy.o\
+      middle_hamilt.o\
+      norm_psi.o\
+      propagator.o\
+      upsi.o\
       FORCE_STRESS.o\
       FORCE_gamma.o\
       FORCE_gamma_edm.o\
       FORCE_gamma_tvnl.o\
       FORCE_gamma_vl.o\
       FORCE_k.o\
-      LCAO_evolve.o\
       LCAO_gen_fixedH.o\
       LCAO_hamilt.o\
       LCAO_matrix.o\

source/module_elecstate/potentials/H_TDDFT_pw.cpp

@@ -2,7 +2,7 @@
 
 #include "module_base/constants.h"
 #include "module_base/timer.h"
-#include "module_hamilt_lcao/module_tddft/ELEC_evolve.h"
+#include "module_hamilt_lcao/module_tddft/evolve_elec.h"
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
 #include "module_io/input.h"
 #include "module_io/input_conv.h"
@@ -11,11 +11,61 @@ namespace elecstate
 {
 
 int H_TDDFT_pw::istep = -1;
-//==========================================================
-// this function aims to add external time-dependent potential
-// (eg: linear potential) used in tddft
-// fuxiang add in 2017-05
-//==========================================================
+
+double H_TDDFT_pw::amp;
+double H_TDDFT_pw::bmod;
+double H_TDDFT_pw::bvec[3];
+
+int H_TDDFT_pw::stype; // 0 : length gauge  1: velocity gauge
+
+std::vector<int> H_TDDFT_pw::ttype;
+//  0  Gauss type function.
+//  1  trapezoid type function.
+//  2  Trigonometric functions, sin^2.
+//  3  heaviside function.
+//  4  HHG function.
+
+int H_TDDFT_pw::tstart;
+int H_TDDFT_pw::tend;
+double H_TDDFT_pw::dt;
+
+// space domain parameters
+
+// length gauge
+double H_TDDFT_pw::lcut1;
+double H_TDDFT_pw::lcut2;
+
+// time domain parameters
+
+// Gauss
+int H_TDDFT_pw::gauss_count;
+std::vector<double> H_TDDFT_pw::gauss_omega; // time(a.u.)^-1
+std::vector<double> H_TDDFT_pw::gauss_phase;
+std::vector<double> H_TDDFT_pw::gauss_sigma; // time(a.u.)
+std::vector<double> H_TDDFT_pw::gauss_t0;
+std::vector<double> H_TDDFT_pw::gauss_amp; // Ry/bohr
+
+// trapezoid
+int H_TDDFT_pw::trape_count;
+std::vector<double> H_TDDFT_pw::trape_omega; // time(a.u.)^-1
+std::vector<double> H_TDDFT_pw::trape_phase;
+std::vector<double> H_TDDFT_pw::trape_t1;
+std::vector<double> H_TDDFT_pw::trape_t2;
+std::vector<double> H_TDDFT_pw::trape_t3;
+std::vector<double> H_TDDFT_pw::trape_amp; // Ry/bohr
+
+// Trigonometric
+int H_TDDFT_pw::trigo_count;
+std::vector<double> H_TDDFT_pw::trigo_omega1; // time(a.u.)^-1
+std::vector<double> H_TDDFT_pw::trigo_omega2; // time(a.u.)^-1
+std::vector<double> H_TDDFT_pw::trigo_phase1;
+std::vector<double> H_TDDFT_pw::trigo_phase2;
+std::vector<double> H_TDDFT_pw::trigo_amp; // Ry/bohr
+
+// Heaviside
+int H_TDDFT_pw::heavi_count;
+std::vector<double> H_TDDFT_pw::heavi_t0;
+std::vector<double> H_TDDFT_pw::heavi_amp; // Ry/bohr
 
 void H_TDDFT_pw::cal_fixed_v(double *vl_pseudo)
 {
@@ -24,10 +74,8 @@ void H_TDDFT_pw::cal_fixed_v(double *vl_pseudo)
     // time evolve
     H_TDDFT_pw::istep++;
 
-    read_parameters(&INPUT);
-
     // judgement to skip vext
-    if (!ELEC_evolve::td_vext || istep > tend || istep < tstart)
+    if (!module_tddft::Evolve_elec::td_vext || istep > tend || istep < tstart)
     {
         return;
     }
@@ -36,13 +84,17 @@ void H_TDDFT_pw::cal_fixed_v(double *vl_pseudo)
     ModuleBase::timer::tick("H_TDDFT_pw", "cal_fixed_v");
 
     int count = 0;
+    gauss_count = 0;
+    trape_count = 0;
+    trigo_count = 0;
+    heavi_count = 0;
 
-    for (auto direc: ELEC_evolve::td_vext_dire_case)
+    for (auto direc: module_tddft::Evolve_elec::td_vext_dire_case)
     {
         std::vector<double> vext_space(this->rho_basis_->nrxx, 0.0);
         double vext_time = cal_v_time(ttype[count]);
 
-        if (ELEC_evolve::out_efield && GlobalV::MY_RANK == 0)
+        if (module_tddft::Evolve_elec::out_efield && GlobalV::MY_RANK == 0)
         {
             std::stringstream as;
             as << GlobalV::global_out_dir << "efield_" << count << ".dat";
@@ -62,80 +114,6 @@ void H_TDDFT_pw::cal_fixed_v(double *vl_pseudo)
     return;
 }
 
-std::vector<double> H_TDDFT_pw::set_parameters(std::string params, double c)
-{
-    std::vector<double> params_ori;
-    std::vector<double> params_out;
-    Input_Conv::parse_expression(params, params_ori);
-    for (auto param: params_ori)
-        params_out.emplace_back(param * c);
-
-    return params_out;
-}
-
-void H_TDDFT_pw::read_parameters(Input *in)
-{
-    stype = in->td_stype;
-
-    Input_Conv::parse_expression(in->td_ttype, ttype);
-
-    tstart = in->td_tstart;
-    tend = in->td_tend;
-
-    dt = in->mdp.md_dt;
-
-    // space domain parameters
-
-    // length gauge
-    lcut1 = in->td_lcut1;
-    lcut2 = in->td_lcut2;
-
-    // time domain parameters
-
-    // Gauss
-    gauss_count = 0;
-    gauss_omega = set_parameters(in->td_gauss_freq, 2 * ModuleBase::PI * ModuleBase::AU_to_FS); // time(a.u.)^-1
-    gauss_phase = set_parameters(in->td_gauss_phase, 1.0);
-    gauss_sigma = set_parameters(in->td_gauss_sigma, 1 / ModuleBase::AU_to_FS);
-    gauss_t0 = set_parameters(in->td_gauss_t0, 1.0);
-    gauss_amp = set_parameters(in->td_gauss_amp, ModuleBase::BOHR_TO_A / ModuleBase::Ry_to_eV); // Ry/bohr
-
-    // trapezoid
-    trape_count = 0;
-    trape_omega = set_parameters(in->td_trape_freq, 2 * ModuleBase::PI * ModuleBase::AU_to_FS); // time(a.u.)^-1
-    trape_phase = set_parameters(in->td_trape_phase, 1.0);
-    trape_t1 = set_parameters(in->td_trape_t1, 1.0);
-    trape_t2 = set_parameters(in->td_trape_t2, 1.0);
-    trape_t3 = set_parameters(in->td_trape_t3, 1.0);
-    trape_amp = set_parameters(in->td_trape_amp, ModuleBase::BOHR_TO_A / ModuleBase::Ry_to_eV); // Ry/bohr
-
-    // Trigonometric
-    trigo_count = 0;
-    trigo_omega1 = set_parameters(in->td_trigo_freq1, 2 * ModuleBase::PI * ModuleBase::AU_to_FS); // time(a.u.)^-1
-    trigo_omega2 = set_parameters(in->td_trigo_freq2, 2 * ModuleBase::PI * ModuleBase::AU_to_FS); // time(a.u.)^-1
-    trigo_phase1 = set_parameters(in->td_trigo_phase1, 1.0);
-    trigo_phase2 = set_parameters(in->td_trigo_phase2, 1.0);
-    trigo_amp = set_parameters(in->td_trigo_amp, ModuleBase::BOHR_TO_A / ModuleBase::Ry_to_eV); // Ry/bohr
-
-    // Heaviside
-    heavi_count = 0;
-    heavi_t0 = set_parameters(in->td_heavi_t0, 1.0);
-    heavi_amp = set_parameters(in->td_heavi_amp, ModuleBase::BOHR_TO_A / ModuleBase::Ry_to_eV); // Ry/bohr
-
-    // HHG
-    // hhg_count = 0;
-    // hhg_amp1 = set_parameters(in->td_hhg_amp1, ModuleBase::BOHR_TO_A / ModuleBase::Ry_to_eV); // Ry/bohr
-    // hhg_amp2 = set_parameters(in->td_hhg_amp2, ModuleBase::BOHR_TO_A / ModuleBase::Ry_to_eV); // Ry/bohr
-    // hhg_omega1 = set_parameters(in->td_hhg_freq1, 2 * ModuleBase::PI * ModuleBase::AU_to_FS); // time(a.u.)^-1
-    // hhg_omega2 = set_parameters(in->td_hhg_freq2, 2 * ModuleBase::PI * ModuleBase::AU_to_FS); // time(a.u.)^-1
-    // hhg_phase1 = set_parameters(in->td_hhg_phase1, 1.0);
-    // hhg_phase2 = set_parameters(in->td_hhg_phase2, 1.0);
-    // hhg_t0 = set_parameters(in->td_hhg_t0, 1.0);
-    // hhg_sigma = set_parameters(in->td_hhg_sigma, 1/ModuleBase::AU_to_FS);
-
-    return;
-}
-
 void H_TDDFT_pw::cal_v_space(std::vector<double> &vext_space, int direc)
 {
     ModuleBase::TITLE("H_TDDFT_pw", "cal_v_space");
@@ -165,11 +143,7 @@ void H_TDDFT_pw::cal_v_space_length(std::vector<double> &vext_space, int direc)
     ModuleBase::TITLE("H_TDDFT_pw", "cal_v_space_length");
     ModuleBase::timer::tick("H_TDDFT_pw", "cal_v_space_length");
 
-    double bmod[3];
-    for (int i = 0; i < 3; i++)
-    {
-        bmod[i] = prepare(GlobalC::ucell, i);
-    }
+    prepare(GlobalC::ucell, direc);
 
     for (int ir = 0; ir < this->rho_basis_->nrxx; ++ir)
     {
@@ -183,15 +157,15 @@ void H_TDDFT_pw::cal_v_space_length(std::vector<double> &vext_space, int direc)
         switch (direc)
         {
         case 1:
-            vext_space[ir] = cal_v_space_length_potential(x) / bmod[0];
+            vext_space[ir] = cal_v_space_length_potential(x) / bmod;
             break;
 
         case 2:
-            vext_space[ir] = cal_v_space_length_potential(y) / bmod[1];
+            vext_space[ir] = cal_v_space_length_potential(y) / bmod;
             break;
 
         case 3:
-            vext_space[ir] = cal_v_space_length_potential(z) / bmod[2];
+            vext_space[ir] = cal_v_space_length_potential(z) / bmod;
             break;
 
         default:
@@ -336,30 +310,8 @@ double H_TDDFT_pw::cal_v_time_heaviside()
     return vext_time;
 }
 
-// double H_TDDFT_pw::cal_v_time_HHG()
-// {
-//     double vext_time = 0.0;
-//     double t0 = *(hhg_t0.begin() + hhg_count);
-//     double omega1 = *(hhg_omega1.begin() + hhg_count);
-//     double phase1 = *(hhg_phase1.begin() + hhg_count);
-//     double omega2 = *(hhg_omega2.begin() + hhg_count);
-//     double phase2 = *(hhg_phase2.begin() + hhg_count);
-//     double amp1 = *(hhg_amp1.begin() + hhg_count);
-//     double amp2 = *(hhg_amp2.begin() + hhg_count);
-//     double sigma = *(trigo_amp2.begin() + trigo_count);
-
-//     double hhg_t = (istep - t0) * dt;
-//     vext_time = (cos(omega1 * hhg_t + phase1) * amp1 + cos(omega2 * hhg_t + phase2) * amp2)
-//                 * exp(-hhg_t * hhg_t * 0.5 / (sigma * sigma));
-//     hhg_count++;
-
-//     return vext_time;
-// }
-
-double H_TDDFT_pw::prepare(const UnitCell &cell, int &dir)
+void H_TDDFT_pw::prepare(const UnitCell& cell, int& dir)
 {
-    double bvec[3] = {0.0};
-    double bmod = 0.0;
     if (dir == 0)
     {
         bvec[0] = cell.G.e11;
@@ -383,7 +335,22 @@ double H_TDDFT_pw::prepare(const UnitCell &cell, int &dir)
         ModuleBase::WARNING_QUIT("H_TDDFT_pw::prepare", "direction is wrong!");
     }
     bmod = sqrt(pow(bvec[0], 2) + pow(bvec[1], 2) + pow(bvec[2], 2));
-    return bmod;
+}
+
+void H_TDDFT_pw ::compute_force(const UnitCell& cell, ModuleBase::matrix& fe)
+{
+    int iat = 0;
+    for (int it = 0; it < cell.ntype; ++it)
+    {
+        for (int ia = 0; ia < cell.atoms[it].na; ++ia)
+        {
+            for (int jj = 0; jj < 3; ++jj)
+            {
+                fe(iat, jj) = ModuleBase::e2 * amp * cell.atoms[it].ncpp.zv * bvec[jj] / bmod;
+            }
+            ++iat;
+        }
+    }
 }
 
 } // namespace elecstate