Fix: optimize lr_spectrum (#5805)

* update the broadening function in lr_spectrum

* optimize transition analysis and fix norm bug

Author: maki49

source/module_lr/esolver_lrtd_lcao.cpp

@@ -568,7 +568,7 @@ void LR::ESolver_LR<T, TR>::after_all_runners(UnitCell& ucell)
         {
             spectrum.optical_absorption_method1(freq, input.abs_broadening);
             // =============================================== for test ====================================================
-            // spectrum.optical_absorption_method2(freq, input.abs_broadening, spin_types[is]);
+            // spectrum.optical_absorption_method2(freq, input.abs_broadening);
             // spectrum.test_transition_dipoles_velocity_ks(eig_ks.c);
             // spectrum.write_transition_dipole(PARAM.globalv.global_out_dir + "dipole_velocity_ks.dat");
             // =============================================== for test ====================================================

source/module_lr/lr_spectrum.cpp

@@ -148,7 +148,8 @@ template<> double LR::LR_Spectrum<double>::cal_mean_squared_dipole(ModuleBase::V
 }
 template<> double LR::LR_Spectrum<std::complex<double>>::cal_mean_squared_dipole(ModuleBase::Vector3<std::complex<double>> dipole)
 {
-    return dipole.norm2().real() / 3.;
+    // return dipole.norm2().real() / 3.;       // ModuleBase::Vector3::norm2 calculates x*x + y*y + z*z, but here we need x*x.conj() + y*y.conj() + z*z.conj()
+    return (std::norm(dipole.x) + std::norm(dipole.y) + std::norm(dipole.z)) / 3.;
 }
 
 template<typename T>
@@ -217,12 +218,12 @@ void LR::LR_Spectrum<T>::transition_analysis(const std::string& spintype)
     ofs << std::setw(40) << spintype << std::endl;
     ofs << "==================================================================== " << std::endl;
     ofs << std::setw(8) << "State" << std::setw(30) << "Excitation Energy (Ry, eV)" <<
-        std::setw(45) << "Transition dipole x, y, z (a.u.)" << std::setw(30) << "Oscillator strength(a.u.)" << std::endl;
+        std::setw(90) << "Transition dipole x, y, z (a.u.)" << std::setw(30) << "Oscillator strength(a.u.)" << std::endl;
     ofs << "------------------------------------------------------------------------------------ " << std::endl;
     for (int istate = 0;istate < nstate;++istate)
         ofs << std::setw(8) << istate << std::setw(15) << std::setprecision(6) << eig[istate] << std::setw(15) << eig[istate] * ModuleBase::Ry_to_eV
-        << std::setw(15) << transition_dipole_[istate].x << std::setw(15) << transition_dipole_[istate].y << std::setw(15) << transition_dipole_[istate].z
-        << std::setw(30) << oscillator_strength_[istate] << std::endl;
+        << std::setprecision(4) << std::setw(30) << transition_dipole_[istate].x << std::setw(30) << transition_dipole_[istate].y << std::setw(30) << transition_dipole_[istate].z
+        << std::setprecision(6) << std::setw(30) << oscillator_strength_[istate] << std::endl;
     ofs << "------------------------------------------------------------------------------------ " << std::endl;
     ofs << std::setw(8) << "State" << std::setw(20) << "Occupied orbital"
         << std::setw(20) << "Virtual orbital" << std::setw(30) << "Excitation amplitude"

source/module_lr/lr_spectrum_velocity.cpp

@@ -23,9 +23,14 @@ namespace LR
         return vR;
     }
 
-    inline double lorentz_delta(const double freq_diff, const double eta)
+    inline double lorentz_delta(const double dfreq_au, const double eta_au)
     {
-        return eta / (freq_diff * freq_diff + eta * eta) / M_PI;
+        return eta_au / (dfreq_au * dfreq_au + eta_au * eta_au) / M_PI;
+    }
+    inline double gauss_delta(const double dfreq_au, const double eta_au)
+    {
+        const double c = eta_au / std::sqrt(2. * std::log(2.));
+        return std::exp(-dfreq_au * dfreq_au / (2 * c * c)) / (std::sqrt(2 * M_PI) * c);
     }
 
     template<typename T> inline ModuleBase::Vector3<T> convert_vector_to_vector3(const std::vector<std::complex<double>>& vec);
@@ -109,13 +114,13 @@ namespace LR
         // 4*pi^2/V * mean_squared_dipole *delta(w-Omega_S)
         std::ofstream ofs(PARAM.globalv.global_out_dir + "absorption.dat");
         if (GlobalV::MY_RANK == 0) { ofs << "Frequency (eV) | wave length(nm) | Absorption (a.u.)" << std::endl; }
-        const double fac = 4 * M_PI * M_PI / ucell.omega * ModuleBase::e2 / this->nk;  // e2: Ry to Hartree in the denominator
+        const double fac = 4 * M_PI * M_PI / ucell.omega / this->nk;
         for (int f = 0;f < freq.size();++f)
         {
             double abs_value = 0.0;
             for (int i = 0;i < nstate;++i)
             {
-                abs_value += this->mean_squared_transition_dipole_[i] * lorentz_delta((freq[f] - eig[i]), eta);
+                abs_value += this->mean_squared_transition_dipole_[i] * lorentz_delta((freq[f] - eig[i]) / ModuleBase::e2, eta / ModuleBase::e2); // e2: Ry to Hartree 
             }
             abs_value *= fac;
             if (GlobalV::MY_RANK == 0) { ofs << freq[f] * ModuleBase::Ry_to_eV << "\t" << 91.126664 / freq[f] << "\t" << abs_value << std::endl; }