Fix : update code for S(t+dt/2) (#2753)

* update code for S(t+dt/2) and fix restart * update autotest * update autotest
deepmodeling · Jul 25, 2023 · 4bf23f6 · 4bf23f6
1 parent c2d273b
commit 4bf23f6
Show file tree

Hide file tree

Showing 26 changed files with 162 additions and 80 deletions.
diff --git a/source/module_esolver/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp
@@ -49,6 +49,14 @@ ESolver_KS_LCAO_TDDFT::~ESolver_KS_LCAO_TDDFT()
         }
         delete Hk_laststep;
     }
+    if (Sk_laststep != nullptr)
+    {
+        for (int ik = 0; ik < kv.nks; ++ik)
+        {
+            delete Sk_laststep[ik];
+        }
+        delete Sk_laststep;
+    }
 }
 
 void ESolver_KS_LCAO_TDDFT::Init(Input& inp, UnitCell& ucell)
@@ -131,6 +139,7 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density(int istep, int iter, double ethr)
                                                  this->psi,
                                                  this->psi_laststep,
                                                  this->Hk_laststep,
+                                                 this->Sk_laststep,
                                                  this->pelec_td->ekb,
                                                  td_htype,
                                                  INPUT.propagator,
@@ -149,6 +158,7 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density(int istep, int iter, double ethr)
                                              this->psi,
                                              this->psi_laststep,
                                              this->Hk_laststep,
+                                             this->Sk_laststep,
                                              this->pelec_td->ekb,
                                              td_htype,
                                              INPUT.propagator,
@@ -335,6 +345,15 @@ void ESolver_KS_LCAO_TDDFT::updatepot(const int istep, const int iter)
                     ModuleBase::GlobalFunc::ZEROS(Hk_laststep[ik], this->LOC.ParaV->nloc);
                 }
             }
+            if (this->Sk_laststep == nullptr)
+            {
+                this->Sk_laststep = new std::complex<double>*[kv.nks];
+                for (int ik = 0; ik < kv.nks; ++ik)
+                {
+                    this->Sk_laststep[ik] = new std::complex<double>[this->LOC.ParaV->nloc];
+                    ModuleBase::GlobalFunc::ZEROS(Sk_laststep[ik], this->LOC.ParaV->nloc);
+                }
+            }
         }
 
         for (int ik = 0; ik < kv.nks; ++ik)
@@ -352,11 +371,12 @@ void ESolver_KS_LCAO_TDDFT::updatepot(const int istep, const int iter)
                 hamilt::MatrixBlock<complex<double>> h_mat, s_mat;
                 this->p_hamilt->matrix(h_mat, s_mat);
                 BlasConnector::copy(this->LOC.ParaV->nloc, h_mat.p, 1, Hk_laststep[ik], 1);
+                BlasConnector::copy(this->LOC.ParaV->nloc, s_mat.p, 1, Sk_laststep[ik], 1);
             }
         }
 
         // calculate energy density matrix for tddft
-        if (istep > 1 && module_tddft::Evolve_elec::td_edm == 0)
+        if (istep >= (wf.init_wfc == "file" ? 0 : 2) && module_tddft::Evolve_elec::td_edm == 0)
             this->cal_edm_tddft();
     }
 

diff --git a/source/module_esolver/esolver_ks_lcao_tddft.h b/source/module_esolver/esolver_ks_lcao_tddft.h
@@ -23,6 +23,7 @@ class ESolver_KS_LCAO_TDDFT : public ESolver_KS_LCAO
 
     psi::Psi<std::complex<double>>* psi_laststep = nullptr;
     std::complex<double>** Hk_laststep = nullptr;
+    std::complex<double>** Sk_laststep = nullptr;
     //same as pelec
     elecstate::ElecStateLCAO_TDDFT* pelec_td = nullptr;
     int td_htype = 1;

diff --git a/source/module_hamilt_lcao/module_tddft/evolve_elec.cpp b/source/module_hamilt_lcao/module_tddft/evolve_elec.cpp
@@ -30,6 +30,7 @@ void Evolve_elec::solve_psi(const int& istep,
                             psi::Psi<std::complex<double>>* psi,
                             psi::Psi<std::complex<double>>* psi_laststep,
                             std::complex<double>** Hk_laststep,
+                            std::complex<double>** Sk_laststep,
                             ModuleBase::matrix& ekb,
                             int htype,
                             int propagator,
@@ -54,6 +55,7 @@ void Evolve_elec::solve_psi(const int& istep,
                        psi[0].get_pointer(),
                        psi_laststep[0].get_pointer(),
                        nullptr,
+                       nullptr,
                        &(ekb(ik, 0)),
                        htype,
                        propagator);
@@ -67,6 +69,7 @@ void Evolve_elec::solve_psi(const int& istep,
                        psi[0].get_pointer(),
                        psi_laststep[0].get_pointer(),
                        Hk_laststep[ik],
+                       Sk_laststep[ik],
                        &(ekb(ik, 0)),
                        htype,
                        propagator);

diff --git a/source/module_hamilt_lcao/module_tddft/evolve_elec.h b/source/module_hamilt_lcao/module_tddft/evolve_elec.h
@@ -47,6 +47,7 @@ class Evolve_elec
                           psi::Psi<std::complex<double>>* psi,
                           psi::Psi<std::complex<double>>* psi_laststep,
                           std::complex<double>** Hk_laststep,
+                          std::complex<double>** Sk_laststep,
                           ModuleBase::matrix& ekb,
                           int htype,
                           int propagator,

diff --git a/source/module_hamilt_lcao/module_tddft/evolve_psi.cpp b/source/module_hamilt_lcao/module_tddft/evolve_psi.cpp
@@ -22,6 +22,7 @@ void evolve_psi(const int nband,
                 std::complex<double>* psi_k,
                 std::complex<double>* psi_k_laststep,
                 std::complex<double>* H_laststep,
+                std::complex<double>* S_laststep,
                 double* ekb,
                 int htype,
                 int propagator)
@@ -58,7 +59,7 @@ void evolve_psi(const int nband,
     /// @output Htmp
     if (htype == 1 && propagator != 2)
     {
-        half_Hmatrix(pv, nband, nlocal, Htmp, H_laststep, print_matrix);
+        half_Hmatrix(pv, nband, nlocal, Htmp, Stmp, H_laststep, S_laststep, print_matrix);
     }
 
     // (2)->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

diff --git a/source/module_hamilt_lcao/module_tddft/evolve_psi.h b/source/module_hamilt_lcao/module_tddft/evolve_psi.h
@@ -18,6 +18,7 @@ void evolve_psi(const int nband,
                 std::complex<double>* psi_k,
                 std::complex<double>* psi_k_laststep,
                 std::complex<double>* H_laststep,
+                std::complex<double>* S_laststep,
                 double* ekb,
                 int htype,
                 int propagator);

diff --git a/source/module_hamilt_lcao/module_tddft/middle_hamilt.cpp b/source/module_hamilt_lcao/module_tddft/middle_hamilt.cpp
@@ -14,7 +14,9 @@ void half_Hmatrix(const Parallel_Orbitals* pv,
                   const int nband,
                   const int nlocal,
                   std::complex<double>* Htmp,
+                  std::complex<double>* Stmp,
                   const std::complex<double>* H_laststep,
+                  const std::complex<double>* S_laststep,
                   const int print_matrix)
 {
     if (print_matrix)
@@ -48,6 +50,7 @@ void half_Hmatrix(const Parallel_Orbitals* pv,
     std::complex<double> alpha = {0.5, 0.0};
     std::complex<double> beta = {0.5, 0.0};
     ScalapackConnector::geadd('N', nlocal, nlocal, alpha, H_laststep, 1, 1, pv->desc, beta, Htmp, 1, 1, pv->desc);
+    ScalapackConnector::geadd('N', nlocal, nlocal, alpha, S_laststep, 1, 1, pv->desc, beta, Stmp, 1, 1, pv->desc);
 
     if (print_matrix)
     {

diff --git a/source/module_hamilt_lcao/module_tddft/middle_hamilt.h b/source/module_hamilt_lcao/module_tddft/middle_hamilt.h
@@ -26,7 +26,9 @@ void half_Hmatrix(const Parallel_Orbitals* pv,
                   const int nband,
                   const int nlocal,
                   std::complex<double>* Htmp,
+                  std::complex<double>* Stmp,
                   const std::complex<double>* H_laststep,
+                  const std::complex<double>* S_laststep,
                   const int print_matrix);
 #endif
 } // namespace module_tddft

diff --git a/source/module_hamilt_lcao/module_tddft/test/middle_hamilt_test.cpp b/source/module_hamilt_lcao/module_tddft/test/middle_hamilt_test.cpp
@@ -31,6 +31,8 @@ TEST(MiddleHamiltTest, testMiddleHamilt)
 {
     std::complex<double>* Htmp;
     std::complex<double>* Hlaststep;
+    std::complex<double>* Stmp;
+    std::complex<double>* Slaststep;
     int nband = 3;
     int nlocal = 4;
     bool print_matrix = false;
@@ -49,6 +51,8 @@ TEST(MiddleHamiltTest, testMiddleHamilt)
     // Initialize data
     Htmp = new std::complex<double>[nlocal * nlocal];
     Hlaststep = new std::complex<double>[nlocal * nlocal];
+    Stmp = new std::complex<double>[nlocal * nlocal];
+    Slaststep = new std::complex<double>[nlocal * nlocal];
 
     for (int i = 0; i < nlocal; ++i)
     {
@@ -58,16 +62,19 @@ TEST(MiddleHamiltTest, testMiddleHamilt)
             {
                 Htmp[i * nlocal + j] = std::complex<double>(1.0, 0.0);
                 Hlaststep[i * nlocal + j] = std::complex<double>(1.0 + 0.2 * (i * nlocal + j), 0.0);
+                Stmp[i * nlocal + j] = std::complex<double>(1.0, 0.0);
+                Slaststep[i * nlocal + j] = std::complex<double>(1.0 + 0.2 * (i * nlocal + j), 0.0);
             }
             else
             {
                 Hlaststep[i * nlocal + j] = std::complex<double>(0.2 * (i * nlocal + j), 0.0);
+                Slaststep[i * nlocal + j] = std::complex<double>(0.2 * (i * nlocal + j), 0.0);
             }
         }
     }
 
     // Call the function
-    module_tddft::half_Hmatrix(pv, nband, nlocal, Htmp, Hlaststep, print_matrix);
+    module_tddft::half_Hmatrix(pv, nband, nlocal, Htmp, Stmp, Hlaststep, Slaststep, print_matrix);
 
     // Check the results
     EXPECT_NEAR(Htmp[0].real(), 1.0, doublethreshold);
@@ -103,6 +110,41 @@ TEST(MiddleHamiltTest, testMiddleHamilt)
     EXPECT_NEAR(Htmp[15].real(), 2.5, doublethreshold);
     EXPECT_NEAR(Htmp[15].imag(), 0.0, doublethreshold);
 
+    EXPECT_NEAR(Stmp[0].real(), 1.0, doublethreshold);
+    EXPECT_NEAR(Stmp[0].imag(), 0.0, doublethreshold);
+    EXPECT_NEAR(Stmp[1].real(), 0.1, doublethreshold);
+    EXPECT_NEAR(Stmp[1].imag(), 0.0, doublethreshold);
+    EXPECT_NEAR(Stmp[2].real(), 0.2, doublethreshold);
+    EXPECT_NEAR(Stmp[2].imag(), 0.0, doublethreshold);
+    EXPECT_NEAR(Stmp[3].real(), 0.3, doublethreshold);
+    EXPECT_NEAR(Stmp[3].imag(), 0.0, doublethreshold);
+    EXPECT_NEAR(Stmp[4].real(), 0.4, doublethreshold);
+    EXPECT_NEAR(Stmp[4].imag(), 0.0, doublethreshold);
+    EXPECT_NEAR(Stmp[5].real(), 1.5, doublethreshold);
+    EXPECT_NEAR(Stmp[5].imag(), 0.0, doublethreshold);
+    EXPECT_NEAR(Stmp[6].real(), 0.6, doublethreshold);
+    EXPECT_NEAR(Stmp[6].imag(), 0.0, doublethreshold);
+    EXPECT_NEAR(Stmp[7].real(), 0.7, doublethreshold);
+    EXPECT_NEAR(Stmp[7].imag(), 0.0, doublethreshold);
+    EXPECT_NEAR(Stmp[8].real(), 0.8, doublethreshold);
+    EXPECT_NEAR(Stmp[8].imag(), 0.0, doublethreshold);
+    EXPECT_NEAR(Stmp[9].real(), 0.9, doublethreshold);
+    EXPECT_NEAR(Stmp[9].imag(), 0.0, doublethreshold);
+    EXPECT_NEAR(Stmp[10].real(), 2.0, doublethreshold);
+    EXPECT_NEAR(Stmp[10].imag(), 0.0, doublethreshold);
+    EXPECT_NEAR(Stmp[11].real(), 1.1, doublethreshold);
+    EXPECT_NEAR(Stmp[11].imag(), 0.0, doublethreshold);
+    EXPECT_NEAR(Stmp[12].real(), 1.2, doublethreshold);
+    EXPECT_NEAR(Stmp[12].imag(), 0.0, doublethreshold);
+    EXPECT_NEAR(Stmp[13].real(), 1.3, doublethreshold);
+    EXPECT_NEAR(Stmp[13].imag(), 0.0, doublethreshold);
+    EXPECT_NEAR(Stmp[14].real(), 1.4, doublethreshold);
+    EXPECT_NEAR(Stmp[14].imag(), 0.0, doublethreshold);
+    EXPECT_NEAR(Stmp[15].real(), 2.5, doublethreshold);
+    EXPECT_NEAR(Stmp[15].imag(), 0.0, doublethreshold);
+
     delete[] Htmp;
     delete[] Hlaststep;
+    delete[] Stmp;
+    delete[] Slaststep;
 }
diff --git a/tests/integrate/601_NO_TDDFT_CO/result.ref b/tests/integrate/601_NO_TDDFT_CO/result.ref
@@ -1,5 +1,5 @@
-etotref -602.9324267719877
-etotperatomref -301.4662133860
-totalforceref 14.600956
-totalstressref 26.941295
-totaltimeref +7.2570
+etotref -602.9325118915726
+etotperatomref -301.4662559458
+totalforceref 14.587060
+totalstressref 26.916557
+totaltimeref +9.1941
diff --git a/tests/integrate/601_NO_TDDFT_CO_occ/result.ref b/tests/integrate/601_NO_TDDFT_CO_occ/result.ref
@@ -1,5 +1,5 @@
-etotref -602.9324267719870
-etotperatomref -301.4662133860
-totalforceref 14.600956
-totalstressref 26.944193
-totaltimeref +7.2778
+etotref -602.9325118915675
+etotperatomref -301.4662559458
+totalforceref 14.587060
+totalstressref 26.917979
+totaltimeref +9.2693
diff --git a/tests/integrate/601_NO_TDDFT_H2/result.ref b/tests/integrate/601_NO_TDDFT_H2/result.ref
@@ -1,5 +1,5 @@
-etotref -18.11864981293731
-etotperatomref -9.0593249065
-totalforceref 44.341410
-totalstressref 78.497930
-totaltimeref +3.9683
+etotref -18.12265885794639
+etotperatomref -9.0613294290
+totalforceref 44.286294
+totalstressref 78.574455
+totaltimeref +5.8103
diff --git a/tests/integrate/601_NO_TDDFT_H2_etrs/result.ref b/tests/integrate/601_NO_TDDFT_H2_etrs/result.ref
@@ -1,5 +1,5 @@
-etotref -18.03711173613717
+etotref -18.03711173613713
 etotperatomref -9.0185558681
 totalforceref 47.878136
 totalstressref 84.826734
-totaltimeref +3.8888
+totaltimeref +5.7771
diff --git a/tests/integrate/601_NO_TDDFT_H2_kpoint/result.ref b/tests/integrate/601_NO_TDDFT_H2_kpoint/result.ref
@@ -1,5 +1,5 @@
-etotref -18.11805148765627
-etotperatomref -9.0590257438
-totalforceref 44.379294
-totalstressref 78.642667
-totaltimeref +3.9332
+etotref -18.12206279865324
+etotperatomref -9.0610313993
+totalforceref 44.324152
+totalstressref 78.718996
+totaltimeref +5.7482
diff --git a/tests/integrate/601_NO_TDDFT_H2_len_gauss/result.ref b/tests/integrate/601_NO_TDDFT_H2_len_gauss/result.ref
@@ -1,5 +1,5 @@
-etotref -31.36573200236431
-etotperatomref -15.6828660012
-totalforceref 0.553188
-totalstressref 5.401078
-totaltimeref +3.9278
+etotref -31.37136608859633
+etotperatomref -15.6856830443
+totalforceref 0.564312
+totalstressref 5.491632
+totaltimeref +5.8870
diff --git a/tests/integrate/601_NO_TDDFT_H2_len_gauss_dire/result.ref b/tests/integrate/601_NO_TDDFT_H2_len_gauss_dire/result.ref
@@ -1,5 +1,5 @@
-etotref -31.36561540267921
-etotperatomref -15.6828077013
-totalforceref 0.551808
-totalstressref 5.401785
-totaltimeref +4.1191
+etotref -31.37124890787322
+etotperatomref -15.6856244539
+totalforceref 0.562932
+totalstressref 5.492325
+totaltimeref +5.9820
diff --git a/tests/integrate/601_NO_TDDFT_H2_len_heavi/result.ref b/tests/integrate/601_NO_TDDFT_H2_len_heavi/result.ref
@@ -1,5 +1,5 @@
-etotref -38.85116687166003
-etotperatomref -19.4255834358
-totalforceref 0.541042
-totalstressref 77.414057
-totaltimeref +3.9307
+etotref -38.85887113911118
+etotperatomref -19.4294355696
+totalforceref 0.554778
+totalstressref 77.526154
+totaltimeref +5.8598
diff --git a/tests/integrate/601_NO_TDDFT_H2_len_hhg/result.ref b/tests/integrate/601_NO_TDDFT_H2_len_hhg/result.ref
@@ -1,5 +1,5 @@
-etotref -82.60939739426242
-etotperatomref -41.3046986971
-totalforceref 0.725730
-totalstressref 501.148069
-totaltimeref +3.9894
+etotref -82.62918679042025
+etotperatomref -41.3145933952
+totalforceref 0.737814
+totalstressref 501.354693
+totaltimeref +5.7690
diff --git a/tests/integrate/601_NO_TDDFT_H2_len_trape/result.ref b/tests/integrate/601_NO_TDDFT_H2_len_trape/result.ref
@@ -1,5 +1,5 @@
-etotref -30.92097376654489
-etotperatomref -15.4604868833
-totalforceref 0.553928
-totalstressref 1.126626
-totaltimeref +3.6563
+etotref -30.92648470928617
+etotperatomref -15.4632423546
+totalforceref 0.564638
+totalstressref 1.215697
+totaltimeref +5.5098
diff --git a/tests/integrate/601_NO_TDDFT_H2_len_trigo/result.ref b/tests/integrate/601_NO_TDDFT_H2_len_trigo/result.ref
@@ -1,5 +1,5 @@
-etotref -30.91150736235756
-etotperatomref -15.4557536812
-totalforceref 0.553934
-totalstressref 1.035650
-totaltimeref +3.6851
+etotref -30.91701568663496
+etotperatomref -15.4585078433
+totalforceref 0.564638
+totalstressref 1.124692
+totaltimeref +5.5228
diff --git a/tests/integrate/601_NO_TDDFT_H2_oldedm/result.ref b/tests/integrate/601_NO_TDDFT_H2_oldedm/result.ref
@@ -1,5 +1,5 @@
-etotref -18.11805148765627
-etotperatomref -9.0590257438
-totalforceref 44.379294
-totalstressref 78.642667
-totaltimeref +3.8941
+etotref -18.12206279865324
+etotperatomref -9.0610313993
+totalforceref 44.324152
+totalstressref 78.718996
+totaltimeref +5.7885
diff --git a/tests/integrate/601_NO_TDDFT_H2_restart/restart/kpoints b/tests/integrate/601_NO_TDDFT_H2_restart/restart/kpoints
@@ -0,0 +1,6 @@
+                               nkstot now = 1
+      KPT             DirectX             DirectY             DirectZ              Weight
+        1                   0                   0                   0                   1
+                                   nkstot = 1                                                            ibzkpt
+      KPT             DirectX             DirectY             DirectZ     IBZ             DirectX             DirectY             DirectZ
+        1                   0                   0                   0       1                   0                   0                   0
diff --git a/tests/integrate/601_NO_TDDFT_H2_restart/result.ref b/tests/integrate/601_NO_TDDFT_H2_restart/result.ref
@@ -1,5 +1,7 @@
-etotref -18.06686250309962
-etotperatomref -9.0334312515
-totalforceref 35.297146
-totalstressref 62.655472
-totaltimeref +8.2247
+etotref -18.06281562534108
+etotperatomref -9.0314078127
+totalforceref 35.324092
+totalstressref 62.632014
+totaltimeref +6.1024
++6.2997
++6.2593
diff --git a/tests/integrate/601_NO_TDDFT_H2_taylor/result.ref b/tests/integrate/601_NO_TDDFT_H2_taylor/result.ref
@@ -1,5 +1,5 @@
-etotref -17.65560405098635
-etotperatomref -8.8278020255
-totalforceref 26.351130
-totalstressref 46.577654
-totaltimeref +4.0140
+etotref -17.66950569471270
+etotperatomref -8.8347528474
+totalforceref 26.535042
+totalstressref 47.003244
+totaltimeref +5.8190