Fix: Hefei-NAMD interface of file syns_nao.csr and parameter cal_syns

docs/advanced/input_files/input-main.md

@@ -3472,8 +3472,9 @@ These variables are used to control molecular dynamics calculations. For more in
 ### cal_syns
 
 - **Type**: Boolean
-- **Description**: Whether the asynchronous overlap matrix is calculated for Hefei-NAMD.
+- **Description**: Whether to calculate and output asynchronous overlap matrix for Hefei-NAMD interface. When enabled, calculates `<phi(t-1)|phi(t)>` by computing overlap between basis functions at atomic positions from previous time step and current time step. The overlap is calculated by shifting atom positions backward by `velocity × md_dt`. Output file: `OUT.*/syns_nao.csr` in CSR format.
 - **Default**: False
+- **Note**: Only works with LCAO basis and molecular dynamics calculations. Requires atomic velocities. Output starts from the second MD step (istep > 0).
 
 ### dmax
 

source/source_esolver/esolver_ks_lcao.cpp

@@ -12,10 +12,10 @@
 #include "../source_lcao/module_ri/exx_opt_orb.h"
 #endif
 #include "source_lcao/module_rdmft/rdmft.h"
-#include "source_estate/module_charge/chgmixing.h" // use charge mixing, mohan add 20251006 
+#include "source_estate/module_charge/chgmixing.h" // use charge mixing, mohan add 20251006
 #include "source_estate/module_dm/init_dm.h" // init dm from electronic wave functions
-#include "source_io/ctrl_runner_lcao.h" // use ctrl_runner_lcao() 
-#include "source_io/ctrl_iter_lcao.h" // use ctrl_iter_lcao() 
+#include "source_io/ctrl_runner_lcao.h" // use ctrl_runner_lcao()
+#include "source_io/ctrl_iter_lcao.h" // use ctrl_iter_lcao()
 #include "source_io/ctrl_scf_lcao.h" // use ctrl_scf_lcao()
 #include "source_io/print_info.h"
 #include "source_lcao/rho_tau_lcao.h" // mohan add 20251024
@@ -554,7 +554,7 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(UnitCell& ucell, const int istep, const
             PARAM.inp, this->kv, this->pelec, this->dmat.dm, this->pv,
             this->gd, this->psi, hamilt_lcao, this->dftu, this->two_center_bundle_,
             this->orb_, this->pw_wfc, this->pw_rho, this->pw_big, this->sf,
-            this->rdmft_solver, this->deepks, this->exx_nao, 
+            this->rdmft_solver, this->deepks, this->exx_nao,
             this->conv_esolver, this->scf_nmax_flag, istep);
 
     //! 3) Clean up RA, which is used to serach for adjacent atoms

source/source_io/ctrl_scf_lcao.cpp

@@ -31,6 +31,7 @@
 #include "source_io/to_qo.h"                // use toQO
 #include "source_lcao/module_rdmft/rdmft.h" // use RDMFT codes
 #include "source_lcao/rho_tau_lcao.h"       // mohan add 2025-10-24
+#include "source_lcao/module_operator_lcao/overlap_new.h" // use hamilt::OverlapNew for NAMD
 
 template <typename TK, typename TR>
 void ModuleIO::ctrl_scf_lcao(UnitCell& ucell,
@@ -459,6 +460,40 @@ void ModuleIO::ctrl_scf_lcao(UnitCell& ucell,
         tqo.calculate();
     }
 
+    //------------------------------------------------------------------
+    //! 18) Calculate and output asynchronous overlap matrix for Hefei-NAMD
+    //------------------------------------------------------------------
+    if (inp.cal_syns && (istep > 0 || inp.init_vel))
+    {
+        ModuleBase::TITLE("ModuleIO", "output_namd_async_overlap");
+        ModuleBase::timer::tick("ModuleIO", "output_namd_async_overlap");
+
+        // Create a new OverlapNew instance specifically for SR_async calculation
+        // This allows SR_async to be initialized with velocity-shifted dtau
+        hamilt::OverlapNew<hamilt::OperatorLCAO<TK, TR>>* overlap_async =
+            new hamilt::OverlapNew<hamilt::OperatorLCAO<TK, TR>>(
+                nullptr,  // hsk_in: not needed for SR_async calculation
+                kv.kvec_d,
+                nullptr,  // hR_in: not needed for SR_async calculation
+                nullptr,  // SR_in: not needed for SR_async calculation
+                &ucell,
+                orb.cutoffs(),
+                &gd,
+                two_center_bundle.overlap_orb.get());
+
+        // Use same precision as DMR output (default 8 if not specified)
+        const int precision = inp.out_dmr[0] > 0 ? inp.out_dmr[1] : 8;
+        const Parallel_Orbitals* paraV = p_hamilt->getSR()->get_paraV();
+        hamilt::HContainer<TR>* SR_async = overlap_async->calculate_SR_async(ucell, PARAM.mdp.md_dt, paraV);
+        overlap_async->output_SR_async_csr(istep, SR_async, precision);
+
+        // Clean up
+        delete SR_async;
+        delete overlap_async;
+
+        ModuleBase::timer::tick("ModuleIO", "output_namd_async_overlap");
+    }
+
     ModuleBase::timer::tick("ModuleIO", "ctrl_scf_lcao");
 }
 

source/source_lcao/FORCE_gamma.cpp

@@ -109,14 +109,6 @@ void Force_LCAO<double>::allocate(const UnitCell& ucell,
                               &gd,
                               nullptr);
 
-    // calculate asynchronous S matrix to output for Hefei-NAMD
-    if (PARAM.inp.cal_syns)
-    {
-        cal_deri = false;
-        ModuleBase::timer::tick("Forces", "allocate");
-        ModuleBase::WARNING_QUIT("cal_syns", "this function has been broken and will be fixed later.");
-    }
-
     ModuleBase::timer::tick("Forces", "allocate");
     return;
 }

source/source_lcao/FORCE_k.cpp

@@ -132,33 +132,6 @@ void Force_LCAO<std::complex<double>>::allocate(const UnitCell& ucell,
                               &gd,
                               nullptr); // delete lm.Hloc_fixedR
 
-    // calculate asynchronous S matrix to output for Hefei-NAMD
-    if (PARAM.inp.cal_syns)
-    {
-        cal_deri = false;
-
-        ModuleBase::WARNING_QUIT("cal_syns", "This function has been broken and will be fixed later.");
-
-        LCAO_domain::build_ST_new(fsr,
-                                  'S',
-                                  cal_deri,
-                                  PARAM.inp.cal_stress,
-                                  ucell,
-                                  orb,
-                                  pv,
-                                  two_center_bundle,
-                                  &(gd),
-                                  nullptr, // delete lm.SlocR
-                                  PARAM.inp.cal_syns,
-                                  PARAM.inp.dmax);
-
-        for (int ik = 0; ik < nks; ik++)
-        {
-
-            bool bit = false; // LiuXh, 2017-03-21
-        }
-    }
-
     ModuleBase::timer::tick("Forces", "allocate");
     return;
 }