Feature: Input and output vector potential in tddft calculation (#4173)

* outHR * update HR output * update * fix out_mat_R * remove GlobalC and GlobalV int td operators * input and output of vector potential * update * Delete abacus_cmake.sh * Update input-main.md * Update H_TDDFT_pw.h * Update LCAO_matrix.cpp * Update td_velocity.h * Update Makefile.Objects * Update for_testing_input_conv.h * Update td_velocity.cpp * Update td_velocity.cpp * Update td_velocity.cpp * Update td_velocity.cpp * Update td_velocity.cpp * clang-format new file * some update * Update td_velocity.h * update --------- Co-authored-by: Mohan Chen <[email protected]>
deepmodeling · May 25, 2024 · 7843f32 · 7843f32
1 parent 3b7c631
commit 7843f32
Show file tree

Hide file tree

Showing 14 changed files with 216 additions and 33 deletions.
diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md
@@ -343,6 +343,8 @@
     - [td\_heavi\_amp](#td_heavi_amp)
     - [out\_dipole](#out_dipole)
     - [out\_efield](#out_efield)
+    - [out\_vecpot](#out_vecpot)
+    - [init\_vecpot\_file](#init_vecpot_file)
     - [ocp](#ocp)
     - [ocp\_set](#ocp_set)
   - [Variables useful for debugging](#variables-useful-for-debugging)
@@ -3262,6 +3264,22 @@ These variables are used to control berry phase and wannier90 interface paramete
   - False: do not output efield.
 - **Default**: False
 
+### out_vecpot
+
+- **Type**: Boolean
+- **Description**: output TDDFT Vector potential or not(a.u.)
+  - True: output Vector potential in file "OUT.suffix/At.dat"
+  - False: do not output Vector potential.
+- **Default**: False
+
+### init_vecpot_file
+
+- **Type**: Boolean
+- **Description**: Init vector potential through file or not
+  - True: init vector potential from file "At.dat".(a.u.) It consists of four columns, representing istep and vector potential on each direction.
+  - False: calculate vector potential by integral of Efield
+- **Default**: False
+
 ### ocp
 
 - **Type**: Boolean

diff --git a/source/module_elecstate/potentials/H_TDDFT_pw.cpp b/source/module_elecstate/potentials/H_TDDFT_pw.cpp
@@ -38,7 +38,7 @@ double H_TDDFT_pw::lcut1;
 double H_TDDFT_pw::lcut2;
 
 //velocity gauge
-double H_TDDFT_pw::At[3]={0.0,0.0,0.0};
+ModuleBase::Vector3<double> H_TDDFT_pw::At;
 
 // time domain parameters
 

diff --git a/source/module_elecstate/potentials/H_TDDFT_pw.h b/source/module_elecstate/potentials/H_TDDFT_pw.h
@@ -54,7 +54,7 @@ class H_TDDFT_pw : public PotBase
     static double lcut2;
 
     //velocity gauge, vector magnetic potential
-    static double At[3];
+    static ModuleBase::Vector3<double> At;
 
     // time domain parameters
 

diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.cpp
@@ -60,7 +60,7 @@ void TDEkinetic<OperatorLCAO<std::complex<double>, double>>::td_ekinetic_scalar(
 template <typename TK, typename TR>
 void TDEkinetic<OperatorLCAO<TK, TR>>::td_ekinetic_grad(std::complex<double>* Hloc, int nnr, ModuleBase::Vector3<double> grad_overlap){
 	std::complex<double> tmp= {0, grad_overlap*cart_At};
-	Hloc[nnr] += tmp;
+	Hloc[nnr] -= tmp;
 	return;
 }
 
@@ -193,10 +193,9 @@ void TDEkinetic<OperatorLCAO<TK, TR>>::init_td(void)
 {
     TD_Velocity::td_vel_op = &td_velocity;
     //calculate At in cartesian coorinates.
-	double l_norm[3]={this->ucell->a1.norm() ,this->ucell->a2.norm() ,this->ucell->a3.norm()};
-    double (&A)[3] = elecstate::H_TDDFT_pw::At;
-	cart_At = this->ucell->a1*A[0]/l_norm[0] + this->ucell->a2*A[1]/l_norm[1] + this->ucell->a3*A[2]/l_norm[2];
-    std::cout << "cart_At: " << cart_At[0] << " " <<cart_At[1]<< " " << cart_At[2] << std::endl;
+    td_velocity.cal_cart_At(this->ucell->a1, this->ucell->a2, this->ucell->a3, elecstate::H_TDDFT_pw::At);
+    this->cart_At = td_velocity.cart_At;
+    std::cout<<"cart_At: "<<cart_At[0]<< " "<<cart_At[1]<< " "<<cart_At[2]<<std::endl;
 
     //init MOT,MGT
     const LCAO_Orbitals& orb = LCAO_Orbitals::get_const_instance();

diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_nonlocal_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_nonlocal_lcao.cpp
@@ -47,10 +47,7 @@ template <typename TK, typename TR>
 void hamilt::TDNonlocal<hamilt::OperatorLCAO<TK, TR>>::init_td(void)
 {
     //calculate At in cartesian coorinates.
-	double l_norm[3]={this->ucell->a1.norm() ,this->ucell->a2.norm() ,this->ucell->a3.norm()};
-    double (&A)[3] = elecstate::H_TDDFT_pw::At;
-	cart_At = -(this->ucell->a1*A[0]/l_norm[0] + this->ucell->a2*A[1]/l_norm[1] + this->ucell->a3*A[2]/l_norm[2]);
-    std::cout << "cart_At: " << cart_At[0] << " " <<cart_At[1]<< " " << cart_At[2] << std::endl;
+	this->cart_At=TD_Velocity::td_vel_op->cart_At;
 }
 // initialize_HR()
 template <typename TK, typename TR>
@@ -207,7 +204,7 @@ void hamilt::TDNonlocal<hamilt::OperatorLCAO<TK, TR>>::calculate_HR()
                                             atom1->iw2n[iw1],
                                             tau0 * this->ucell->lat0,
                                             T0,
-                                            -cart_At,
+                                            -cart_At/2.0,
                                             0);
 #else
                 uot.snap_psibeta_half_tddft(orb,
@@ -220,7 +217,7 @@ void hamilt::TDNonlocal<hamilt::OperatorLCAO<TK, TR>>::calculate_HR()
                                             atom1->iw2n[iw1],
                                             tau0 * this->ucell->lat0,
                                             T0,
-                                            -cart_At,
+                                            -cart_At/2.0,
                                             0);
 #endif
                 nlm_tot[ad].insert({all_indexes[iw1l], nlm[0]});

diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.cpp
@@ -376,4 +376,4 @@ void sparse_format::clear_zero_elements(
     }
 
     return;
-}
+}
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.h b/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.h
@@ -46,4 +46,4 @@ namespace sparse_format
 
 }
 
-#endif
+#endif
diff --git a/source/module_hamilt_lcao/module_tddft/td_velocity.cpp b/source/module_hamilt_lcao/module_tddft/td_velocity.cpp
@@ -1,24 +1,145 @@
-#include "module_base/timer.h"
 #include "td_velocity.h"
 
+#include "module_base/timer.h"
+#include "module_elecstate/potentials/H_TDDFT_pw.h"
 
 bool TD_Velocity::tddft_velocity = false;
 bool TD_Velocity::out_mat_R = false;
+bool TD_Velocity::out_vecpot = false;
+bool TD_Velocity::init_vecpot_file = false;
+
 TD_Velocity* TD_Velocity::td_vel_op = nullptr;
 
+int TD_Velocity::istep = -1;
+int TD_Velocity::max_istep = -1;
+std::vector<ModuleBase::Vector3<double>> TD_Velocity::At_from_file;
+
 TD_Velocity::TD_Velocity()
 {
-  return;
+    if (init_vecpot_file && istep == -1)
+    {
+        this->read_cart_At();
+    }
+    return;
 }
 TD_Velocity::~TD_Velocity()
 {
-  this->destroy_HS_R_td_sparse();
+    this->destroy_HS_R_td_sparse();
+}
+
+void TD_Velocity::output_cart_At(const std::string& out_dir)
+{
+    if (GlobalV::MY_RANK == 0)
+    {
+        std::string out_file;
+        // generate the output file name
+        out_file = out_dir + "At.dat";
+        std::ofstream ofs;
+        // output title
+        if (istep == 0)
+        {
+            ofs.open(out_file.c_str(), std::ofstream::out);
+            ofs << std::left << std::setw(8) << "#istep" << std::setw(15) << "A_x" << std::setw(15) << "A_y"
+                << std::setw(15) << "A_z" << std::endl;
+        }
+        else
+        {
+            ofs.open(out_file.c_str(), std::ofstream::app);
+        }
+        // output the vector potential
+        ofs << std::left << std::setw(8) << istep;
+        // divide by 2.0 to get the atomic unit
+        for (int i = 0; i < 3; i++)
+        {
+            ofs << std::scientific << std::setprecision(4) << std::setw(15) << cart_At[i] / 2.0;
+        }
+        ofs << std::endl;
+        ofs.close();
+    }
+    return;
 }
 
+void TD_Velocity::cal_cart_At(const ModuleBase::Vector3<double>& a0,
+                              const ModuleBase::Vector3<double>& a1,
+                              const ModuleBase::Vector3<double>& a2,
+                              const ModuleBase::Vector3<double>& At)
+{
+    istep++;
+    if (init_vecpot_file)
+    {
+        this->cart_At = At_from_file[istep > max_istep ? max_istep : istep];
+    }
+    else
+    {
+        const double l_norm[3] = {a0.norm(), a1.norm(), a2.norm()};
+        this->cart_At = a0 * At[0] / l_norm[0] + a1 * At[1] / l_norm[1] + a2 * At[2] / l_norm[2];
+    }
+    // output the vector potential if needed
+    if (out_vecpot == true)
+    {
+        this->output_cart_At(GlobalV::global_out_dir);
+    }
+}
+
+void TD_Velocity::read_cart_At(void)
+{
+    std::string in_file;
+    // generate the input file name
+    in_file = "At.dat";
+    std::ifstream ifs(in_file.c_str());
+    // check if the file is exist
+    if (!ifs)
+    {
+        ModuleBase::WARNING_QUIT("TD_Velocity::read_cart_At", "Cannot open Vector potential file!");
+    }
+    std::string line;
+    std::vector<std::string> str_vec;
+    // use tmp to skip the istep number
+    int tmp = 0;
+    while (std::getline(ifs, line))
+    {
+        // A tmporary vector3 to store the data of this line
+        ModuleBase::Vector3<double> At;
+        if (line[0] == '#')
+        {
+            continue;
+        }
+        std::istringstream iss(line);
+        // skip the istep number
+        if (!(iss >> tmp))
+        {
+            ModuleBase::WARNING_QUIT("TD_Velocity::read_cart_At", "Error reading istep!");
+        }
+        // read the vector potential
+        double component = 0;
+        // Read three components
+        for (int i = 0; i < 3; i++)
+        {
+            if (!(iss >> component))
+            {
+                ModuleBase::WARNING_QUIT("TD_Velocity::read_cart_At",
+                                         "Error reading component " + std::to_string(i + 1) + " for istep "
+                                             + std::to_string(tmp) + "!");
+            }
+            At[i] = component;
+        }
+        // unit transform ,change unit into Ry/bohr/e*t_a.u.
+        At *= 2.0;
+        // add the tmporary vector3 to the vector potential vector
+        At_from_file.push_back(At);
+    }
+    // set the max_istep
+    max_istep = At_from_file.size() - 1;
+    ifs.close();
+
+    return;
+}
 void TD_Velocity::destroy_HS_R_td_sparse(void)
 {
-  std::map<Abfs::Vector3_Order<int>, std::map<size_t, std::map<size_t, std::complex<double>>>> empty_HR_sparse_td_vel_up;
-  std::map<Abfs::Vector3_Order<int>, std::map<size_t, std::map<size_t, std::complex<double>>>> empty_HR_sparse_td_vel_down;
-  HR_sparse_td_vel[0].swap(empty_HR_sparse_td_vel_up);
-  HR_sparse_td_vel[1].swap(empty_HR_sparse_td_vel_down);
-}
+    std::map<Abfs::Vector3_Order<int>, std::map<size_t, std::map<size_t, std::complex<double>>>>
+        empty_HR_sparse_td_vel_up;
+    std::map<Abfs::Vector3_Order<int>, std::map<size_t, std::map<size_t, std::complex<double>>>>
+        empty_HR_sparse_td_vel_down;
+    HR_sparse_td_vel[0].swap(empty_HR_sparse_td_vel_up);
+    HR_sparse_td_vel[1].swap(empty_HR_sparse_td_vel_down);
+}
diff --git a/source/module_hamilt_lcao/module_tddft/td_velocity.h b/source/module_hamilt_lcao/module_tddft/td_velocity.h
@@ -1,17 +1,19 @@
 #ifndef TD_VELOCITY_H
 #define TD_VELOCITY_H
-#include <map>
-
-#include "module_base/timer.h"
 #include "module_base/abfs-vector3_order.h"
-//Class to store TDDFT velocity gague infos.
+#include "module_base/timer.h"
+
+#include <map>
+// Class to store TDDFT velocity gague infos.
 class TD_Velocity
 {
   public:
     TD_Velocity();
     ~TD_Velocity();
 
-    /// @brief Judge if in tddft calculation or not 
+    void init();
+
+    /// @brief Judge if in tddft calculation or not
     static bool tddft_velocity;
 
     /// @brief switch to control the output of HR
@@ -20,15 +22,42 @@ class TD_Velocity
     /// @brief pointer to the only TD_Velocity object itself
     static TD_Velocity* td_vel_op;
 
-    //For TDDFT velocity gague, to fix the output of HR
-    std::map<Abfs::Vector3_Order<int>, std::map<size_t, std::map<size_t, std::complex<double>>>> HR_sparse_td_vel[2];
+    /// @brief switch to control the output of At
+    static bool out_vecpot;
+
+    /// @brief switch to control the source of At
+    static bool init_vecpot_file;
 
+    /// @brief Store the vector potential for tddft calculation
+    ModuleBase::Vector3<double> cart_At;
+
+    /// @brief calculate the At in cartesian coordinate
+    void cal_cart_At(const ModuleBase::Vector3<double>& a0,
+                     const ModuleBase::Vector3<double>& a1,
+                     const ModuleBase::Vector3<double>& a2,
+                     const ModuleBase::Vector3<double>& At);
+
+    // For TDDFT velocity gague, to fix the output of HR
+    std::map<Abfs::Vector3_Order<int>, std::map<size_t, std::map<size_t, std::complex<double>>>> HR_sparse_td_vel[2];
 
   private:
+    /// @brief read At from output file
+    void read_cart_At(void);
+
+    /// @brief output cart_At to output file
+    void output_cart_At(const std::string& out_dir);
+
+    /// @brief store isteps now
+    static int istep;
+
+    /// @brief total steps of read in At
+    static int max_istep;
+
+    /// @brief store the read in At_data
+    static std::vector<ModuleBase::Vector3<double>> At_from_file;
 
     /// @brief destory HSR data stored
     void destroy_HS_R_td_sparse(void);
-
 };
 
-#endif
+#endif
diff --git a/source/module_io/input.cpp b/source/module_io/input.cpp
@@ -470,6 +470,8 @@ void Input::Default(void)
     out_dipole = false;
     out_efield = false;
     out_current = false;
+    out_vecpot = false;
+    init_vecpot_file = false;
 
     td_print_eij = -1.0;
     td_edm = 0;
@@ -1813,6 +1815,14 @@ bool Input::Read(const std::string& fn)
         {
             read_value(ifs, out_efield);
         }
+        else if (strcmp("out_vecpot", word) == 0)
+        {
+            read_value(ifs, out_vecpot);
+        }
+        else if (strcmp("init_vecpot_file", word) == 0)
+        {
+            read_value(ifs, init_vecpot_file);
+        }
         else if (strcmp("td_print_eij", word) == 0)
         {
             read_value(ifs, td_print_eij);
@@ -3752,6 +3762,8 @@ void Input::Bcast()
     Parallel_Common::bcast_bool(out_dipole);
     Parallel_Common::bcast_bool(out_efield);
     Parallel_Common::bcast_bool(out_current);
+    Parallel_Common::bcast_bool(out_vecpot);
+    Parallel_Common::bcast_bool(init_vecpot_file);
     Parallel_Common::bcast_double(td_print_eij);
     Parallel_Common::bcast_int(td_edm);
     Parallel_Common::bcast_bool(test_skip_ewald);

diff --git a/source/module_io/input.h b/source/module_io/input.h
@@ -416,6 +416,8 @@ class Input
     bool out_dipole; // output the dipole or not
     bool out_efield; // output the efield or not
     bool out_current; //output the current or not
+    bool out_vecpot; // output the vector potential or not
+    bool init_vecpot_file; // initialize the vector potential, though file or integral
 
     double td_print_eij; // threshold to output Eij elements
     int td_edm; //0: new edm method   1: old edm method

diff --git a/source/module_io/input_conv.cpp b/source/module_io/input_conv.cpp
@@ -550,6 +550,8 @@ void Input_Conv::Convert(void)
     module_tddft::Evolve_elec::out_current = INPUT.out_current;
     module_tddft::Evolve_elec::td_print_eij = INPUT.td_print_eij;
     module_tddft::Evolve_elec::td_edm = INPUT.td_edm;
+    TD_Velocity::out_vecpot = INPUT.out_vecpot;
+    TD_Velocity::init_vecpot_file = INPUT.init_vecpot_file;
     read_td_efield();
 #endif