@@ -620,10 +620,13 @@ void LCAO_Hamilt::calculate_STN_R_sparse(const int ¤t_spin, const double &
620620
621621 if (GlobalV::NSPIN!=4 )
622622 {
623- temp_value_double = GlobalC::LM.SlocR [index];
624- if (std::abs (temp_value_double) > sparse_threshold)
623+ if (current_spin == 0 )
625624 {
626- GlobalC::LM.SR_sparse [dR][iw1_all][iw2_all] = temp_value_double;
625+ temp_value_double = GlobalC::LM.SlocR [index];
626+ if (std::abs (temp_value_double) > sparse_threshold)
627+ {
628+ GlobalC::LM.SR_sparse [dR][iw1_all][iw2_all] = temp_value_double;
629+ }
627630 }
628631
629632 temp_value_double = GlobalC::LM.Hloc_fixedR [index];
@@ -659,6 +662,113 @@ void LCAO_Hamilt::calculate_STN_R_sparse(const int ¤t_spin, const double &
659662}
660663
661664
665+ void LCAO_Hamilt::calculate_STN_R_sparse_for_S (const double &sparse_threshold)
666+ {
667+ ModuleBase::TITLE (" LCAO_Hamilt" ," calculate_STN_R_sparse_for_S" );
668+
669+ int index = 0 ;
670+ ModuleBase::Vector3<double > dtau, tau1, tau2;
671+ ModuleBase::Vector3<double > dtau1, dtau2, tau0;
672+
673+ double temp_value_double;
674+ std::complex <double > temp_value_complex;
675+
676+ for (int T1 = 0 ; T1 < GlobalC::ucell.ntype ; ++T1)
677+ {
678+ Atom* atom1 = &GlobalC::ucell.atoms [T1];
679+ for (int I1 = 0 ; I1 < atom1->na ; ++I1)
680+ {
681+ tau1 = atom1->tau [I1];
682+ GlobalC::GridD.Find_atom (GlobalC::ucell, tau1, T1, I1);
683+ Atom* atom1 = &GlobalC::ucell.atoms [T1];
684+ const int start = GlobalC::ucell.itiaiw2iwt (T1,I1,0 );
685+
686+ for (int ad = 0 ; ad < GlobalC::GridD.getAdjacentNum ()+1 ; ++ad)
687+ {
688+ const int T2 = GlobalC::GridD.getType (ad);
689+ const int I2 = GlobalC::GridD.getNatom (ad);
690+ Atom* atom2 = &GlobalC::ucell.atoms [T2];
691+
692+ tau2 = GlobalC::GridD.getAdjacentTau (ad);
693+ dtau = tau2 - tau1;
694+ double distance = dtau.norm () * GlobalC::ucell.lat0 ;
695+ double rcut = GlobalC::ORB.Phi [T1].getRcut () + GlobalC::ORB.Phi [T2].getRcut ();
696+
697+ bool adj = false ;
698+
699+ if (distance < rcut) adj = true ;
700+ else if (distance >= rcut)
701+ {
702+ for (int ad0 = 0 ; ad0 < GlobalC::GridD.getAdjacentNum ()+1 ; ++ad0)
703+ {
704+ const int T0 = GlobalC::GridD.getType (ad0);
705+
706+ tau0 = GlobalC::GridD.getAdjacentTau (ad0);
707+ dtau1 = tau0 - tau1;
708+ dtau2 = tau0 - tau2;
709+
710+ double distance1 = dtau1.norm () * GlobalC::ucell.lat0 ;
711+ double distance2 = dtau2.norm () * GlobalC::ucell.lat0 ;
712+
713+ double rcut1 = GlobalC::ORB.Phi [T1].getRcut () + GlobalC::ucell.infoNL .Beta [T0].get_rcut_max ();
714+ double rcut2 = GlobalC::ORB.Phi [T2].getRcut () + GlobalC::ucell.infoNL .Beta [T0].get_rcut_max ();
715+
716+ if ( distance1 < rcut1 && distance2 < rcut2 )
717+ {
718+ adj = true ;
719+ break ;
720+ }
721+ }
722+ }
723+
724+ if (adj)
725+ {
726+ const int start2 = GlobalC::ucell.itiaiw2iwt (T2,I2,0 );
727+
728+ Abfs::Vector3_Order<int > dR (GlobalC::GridD.getBox (ad).x , GlobalC::GridD.getBox (ad).y , GlobalC::GridD.getBox (ad).z );
729+
730+ for (int ii=0 ; ii<atom1->nw *GlobalV::NPOL; ii++)
731+ {
732+ const int iw1_all = start + ii;
733+ const int mu = GlobalC::ParaO.trace_loc_row [iw1_all];
734+
735+ if (mu<0 )continue ;
736+
737+ for (int jj=0 ; jj<atom2->nw *GlobalV::NPOL; jj++)
738+ {
739+ int iw2_all = start2 + jj;
740+ const int nu = GlobalC::ParaO.trace_loc_col [iw2_all];
741+
742+ if (nu<0 )continue ;
743+
744+ if (GlobalV::NSPIN!=4 )
745+ {
746+ temp_value_double = GlobalC::LM.SlocR [index];
747+ if (std::abs (temp_value_double) > sparse_threshold)
748+ {
749+ GlobalC::LM.SR_sparse [dR][iw1_all][iw2_all] = temp_value_double;
750+ }
751+ }
752+ else
753+ {
754+ temp_value_complex = GlobalC::LM.SlocR_soc [index];
755+ if (std::abs (temp_value_complex) > sparse_threshold)
756+ {
757+ GlobalC::LM.SR_soc_sparse [dR][iw1_all][iw2_all] = temp_value_complex;
758+ }
759+ }
760+
761+ ++index;
762+ }
763+ }
764+ }
765+ }
766+ }
767+ }
768+
769+ return ;
770+ }
771+
662772void LCAO_Hamilt::calculate_HSR_sparse (const int ¤t_spin, const double &sparse_threshold)
663773{
664774 ModuleBase::TITLE (" LCAO_Hamilt" ," calculate_HSR_sparse" );
@@ -692,6 +802,13 @@ void LCAO_Hamilt::calculate_HSR_sparse(const int ¤t_spin, const double &sp
692802
693803}
694804
805+ void LCAO_Hamilt::calculate_SR_sparse (const double &sparse_threshold)
806+ {
807+ ModuleBase::TITLE (" LCAO_Hamilt" ," calculate_SR_sparse" );
808+ set_R_range_sparse ();
809+ calculate_STN_R_sparse_for_S (sparse_threshold);
810+ }
811+
695812void LCAO_Hamilt::calculat_HR_dftu_sparse (const int ¤t_spin, const double &sparse_threshold)
696813{
697814 ModuleBase::TITLE (" LCAO_Hamilt" ," calculat_HR_dftu_sparse" );
0 commit comments