fix the spectra strength of singlet

Author: maki49

source/module_lr/esolver_lrtd_lcao.cpp

@@ -564,12 +564,15 @@ void LR::ESolver_LR<T, TR>::after_all_runners(UnitCell& ucell)
             &this->pelec->ekb.c[is * nstates], this->X[is].template data<T>(), nstates, openshell,
             LR_Util::tolower(input.abs_gauge));
         spectrum.transition_analysis(spin_types[is]);
-        spectrum.optical_absorption_method1(freq, input.abs_broadening, spin_types[is]);
-        // =============================================== for test ====================================================
-        // spectrum.optical_absorption_method2(freq, input.abs_broadening, spin_types[is]);
-        // spectrum.test_transition_dipoles_velocity_ks(eig_ks.c);
-        // spectrum.write_transition_dipole(PARAM.globalv.global_out_dir + "dipole_velocity_ks_" + spin_types[is] + ".dat");
-        // =============================================== for test ====================================================
+        if (spin_types[is] != "triplet")        // triplets has no transition dipole and no contribution to the spectrum
+        {
+            spectrum.optical_absorption_method1(freq, input.abs_broadening);
+            // =============================================== for test ====================================================
+            // spectrum.optical_absorption_method2(freq, input.abs_broadening, spin_types[is]);
+            // 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

@@ -81,6 +81,7 @@ ModuleBase::Vector3<double> LR::LR_Spectrum<double>::cal_transition_dipole_istat
     }
     trans_dipole *= (ucell.omega / static_cast<double>(gint->get_ncxyz()));   // dv
     trans_dipole *= static_cast<double>(this->nk);  // nk is divided inside DM_trans, now recover it
+    if (this->nspin_x == 1) { trans_dipole *= sqrt(2.0); } // *2 for 2 spins, /sqrt(2) for the halfed dimension of X in the normalizaiton
     Parallel_Reduce::reduce_all(trans_dipole.x);
     Parallel_Reduce::reduce_all(trans_dipole.y);
     Parallel_Reduce::reduce_all(trans_dipole.z);
@@ -134,6 +135,7 @@ ModuleBase::Vector3<std::complex<double>> LR::LR_Spectrum<std::complex<double>>:
     }
     trans_dipole *= (ucell.omega / static_cast<double>(gint->get_ncxyz()));   // dv
     trans_dipole *= static_cast<double>(this->nk);  // nk is divided inside DM_trans, now recover it
+    if (this->nspin_x == 1) { trans_dipole *= sqrt(2.0); } // *2 for 2 spins, /sqrt(2) for the halfed dimension of X in the normalizaiton
     Parallel_Reduce::reduce_all(trans_dipole.x);
     Parallel_Reduce::reduce_all(trans_dipole.y);
     Parallel_Reduce::reduce_all(trans_dipole.z);
@@ -177,21 +179,21 @@ void LR::LR_Spectrum<T>::oscillator_strength()
 }
 
 template<typename T>
-void LR::LR_Spectrum<T>::optical_absorption_method1(const std::vector<double>& freq, const double eta, const std::string& spintype)
+void LR::LR_Spectrum<T>::optical_absorption_method1(const std::vector<double>& freq, const double eta)
 {
     // ============test dipole================
     // this->cal_transition_dipoles_length();
-    // this->write_transition_dipole(PARAM.globalv.global_out_dir + "dipole_length_" + spintype + ".dat");
+    // this->write_transition_dipole(PARAM.globalv.global_out_dir + "dipole_length.dat");
     // this->cal_transition_dipoles_velocity();
-    // this->write_transition_dipole(PARAM.globalv.global_out_dir + "dipole_velocity_" + spintype + ".dat");
+    // this->write_transition_dipole(PARAM.globalv.global_out_dir + "dipole_velocity.dat");
     // exit(0);
     // ============test dipole================
     ModuleBase::TITLE("LR::LR_Spectrum", "optical_absorption");
     // 4*pi^2/V * mean_squared_dipole *delta(w-Omega_S)
     // = -8*pi*Omega_S/V * mean_squared_dipole * Im[1/[(w+i\eta)^2-\Omega_S^2]]
     // = -4*pi/V * oscilator_strength * Im[1/[(w+i\eta)^2-\Omega_S^2]]
     std::vector<double>& osc = this->oscillator_strength_;
-    std::ofstream ofs(PARAM.globalv.global_out_dir + "absorption_" + spintype + ".dat");
+    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; }
     double FourPI_div_c = ModuleBase::FOUR_PI / 137.036;
     double fac = 4 * M_PI / ucell.omega * ModuleBase::e2 / this->nk;   // e2 for Ry to Hartree in the denominator

source/module_lr/lr_spectrum.h

@@ -32,9 +32,9 @@ namespace LR
             this->oscillator_strength();
         };
         /// @brief calculate the optical absorption spectrum with $Im[1/[(w+i\eta)^2-\Omega_S^2]]$
-        void optical_absorption_method1(const std::vector<double>& freq, const double eta, const std::string& spintype);
+        void optical_absorption_method1(const std::vector<double>& freq, const double eta);
         /// @brief calculate the optical absorption spectrum with lorentzian delta function
-        void optical_absorption_method2(const std::vector<double>& freq, const double eta, const std::string& spintype);
+        void optical_absorption_method2(const std::vector<double>& freq, const double eta);
         /// @brief print out the transition dipole moment and the main contributions to the transition amplitude
         void transition_analysis(const std::string& spintype);
 

source/module_lr/lr_spectrum_velocity.cpp

@@ -56,6 +56,7 @@ namespace LR
                 trans_dipole[i] += LR_Util::dot_R_matrix(*vR.get_current_term_pointer(i), *DM_trans.get_DMR_pointer(is + 1), ucell.nat) * fac;
             }   // end for spin_x, only matter in open-shell system
             trans_dipole[i] *= static_cast<double>(this->nk);  // nk is divided inside DM_trans, now recover it
+            if (this->nspin_x == 1) { trans_dipole[i] *= sqrt(2.0); } // *2 for 2 spins, /sqrt(2) for the halfed dimension of X in the normalizaiton
             Parallel_Reduce::reduce_all(trans_dipole[i]);
         }   // end for direction
         return convert_vector_to_vector3<T>(trans_dipole);
@@ -82,6 +83,7 @@ namespace LR
                 }
             }   // end for spin_x, only matter in open-shell system
             trans_dipole[i] *= static_cast<double>(this->nk);  // nk is divided inside DM_trans, now recover it
+            if (this->nspin_x == 1) { trans_dipole[i] *= sqrt(2.0); } // *2 for 2 spins, /sqrt(2) for the halfed dimension of X in the normalizaiton
             Parallel_Reduce::reduce_all(trans_dipole[i]);
         }   // end for direction
         return convert_vector_to_vector3<T>(trans_dipole);
@@ -101,11 +103,11 @@ namespace LR
     }
 
     template<typename T>
-    void LR::LR_Spectrum<T>::optical_absorption_method2(const std::vector<double>& freq, const double eta, const std::string& spintype)
+    void LR::LR_Spectrum<T>::optical_absorption_method2(const std::vector<double>& freq, const double eta)
     {
         ModuleBase::TITLE("LR::LR_Spectrum", "optical_absorption_velocity");
         // 4*pi^2/V * mean_squared_dipole *delta(w-Omega_S)
-        std::ofstream ofs(PARAM.globalv.global_out_dir + "absorption_" + spintype + ".dat");
+        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
         for (int f = 0;f < freq.size();++f)