@@ -47,95 +47,149 @@ void compose_hr_gint(HContainer<double>& hr_gint)
4747 ModuleBase::timer::tick (" Gint" , " compose_hr_gint" );
4848}
4949
50- void compose_hr_gint ( const std::vector<HContainer< double >>& hr_gint_part,
51- HContainer<std:: complex < double >>& hr_gint_full )
50+ template < typename T>
51+ void transfer_hr_gint_to_hR ( const HContainer<T>& hr_gint, HContainer<T>& hR )
5252{
53- ModuleBase::TITLE (" Gint" , " compose_hr_gint" );
54- ModuleBase::timer::tick (" Gint" , " compose_hr_gint" );
55- for (int iap = 0 ; iap < hr_gint_full.size_atom_pairs (); iap++)
53+ ModuleBase::TITLE (" Gint" , " transfer_hr_gint_to_hR" );
54+ ModuleBase::timer::tick (" Gint" , " transfer_hr_gint_to_hR" );
55+ #ifdef __MPI
56+ int size = 0 ;
57+ MPI_Comm_size (MPI_COMM_WORLD, &size);
58+ if (size == 1 )
59+ {
60+ hR.add (hr_gint);
61+ }
62+ else
5663 {
57- auto * ap = &(hr_gint_full.get_atom_pair (iap));
58- const int iat1 = ap->get_atom_i ();
59- const int iat2 = ap->get_atom_j ();
60- if (iat1 <= iat2)
64+ hamilt::transferSerials2Parallels (hr_gint, &hR);
65+ }
66+ #else
67+ hR.add (hr_gint);
68+ #endif
69+ ModuleBase::timer::tick (" Gint" , " transfer_hr_gint_to_hR" );
70+ }
71+
72+
73+ void merge_hr_part_to_hR (const std::vector<hamilt::HContainer<double >>& hRGint_tmp,
74+ hamilt::HContainer<std::complex <double >>* hR,
75+ const GintInfo& gint_info){
76+ ModuleBase::TITLE (" Gint_k" , " transfer_pvpR" );
77+ ModuleBase::timer::tick (" Gint_k" , " transfer_pvpR" );
78+
79+ int mg = hR->get_paraV ()->get_global_row_size ()/2 ;
80+ int ng = hR->get_paraV ()->get_global_col_size ()/2 ;
81+ int nb = hR->get_paraV ()->get_block_size ()/2 ;
82+ const UnitCell* ucell_in = gint_info.get_ucell ();
83+ #ifdef __MPI
84+ int blacs_ctxt = hR->get_paraV ()->blacs_ctxt ;
85+ std::vector<int > iat2iwt (ucell_in->nat );
86+ for (int iat = 0 ; iat < ucell_in->nat ; iat++) {
87+ iat2iwt[iat] = ucell_in->get_iat2iwt ()[iat]/2 ;
88+ }
89+ Parallel_Orbitals *pv = new Parallel_Orbitals ();
90+ pv->set (mg, ng, nb, blacs_ctxt);
91+ pv->set_atomic_trace (iat2iwt.data (), ucell_in->nat , mg);
92+ auto ijr_info = hR->get_ijr_info ();
93+
94+ auto * hR_tmp = new hamilt::HContainer<std::complex <double >>(pv, nullptr , &ijr_info);
95+
96+ // select hRGint_tmp
97+ std::vector<int > first = {0 , 1 , 1 , 0 };
98+ std::vector<int > second= {3 , 2 , 2 , 3 };
99+ // select position in the big matrix
100+ std::vector<int > row_set = {0 , 0 , 1 , 1 };
101+ std::vector<int > col_set = {0 , 1 , 0 , 1 };
102+ // construct complex matrix
103+ std::vector<int > clx_i = {1 , 0 , 0 , -1 };
104+ std::vector<int > clx_j = {0 , 1 , -1 , 0 };
105+ for (int is = 0 ; is < 4 ; is++){
106+ if (!PARAM.globalv .domag && (is==1 || is==2 )) continue ;
107+ hR_tmp->set_zero ();
108+ hamilt::HContainer<std::complex <double >>* hRGint_tmpCd = new hamilt::HContainer<std::complex <double >>(ucell_in->nat );
109+ hRGint_tmpCd->insert_ijrs (&ijr_info, *ucell_in);
110+ hRGint_tmpCd->allocate (nullptr , true );
111+ hRGint_tmpCd->set_zero ();
112+ for (int iap = 0 ; iap < hRGint_tmpCd->size_atom_pairs (); iap++)
61113 {
62- hamilt::AtomPair<std::complex <double >>* upper_ap = ap;
63- hamilt::AtomPair<std::complex <double >>* lower_ap = hr_gint_full.find_pair (iat2, iat1);
64- const hamilt::AtomPair<double >* ap_nspin_0 = hr_gint_part[0 ].find_pair (iat1, iat2);
65- const hamilt::AtomPair<double >* ap_nspin_3 = hr_gint_part[3 ].find_pair (iat1, iat2);
66- for (int ir = 0 ; ir < upper_ap->get_R_size (); ir++)
114+ auto * ap = &hRGint_tmpCd->get_atom_pair (iap);
115+ const int iat1 = ap->get_atom_i ();
116+ const int iat2 = ap->get_atom_j ();
117+ if (iat1 <= iat2)
67118 {
68- const auto R_index = upper_ap->get_R_index (ir);
69- auto upper_mat = upper_ap->find_matrix (R_index);
70- auto mat_nspin_0 = ap_nspin_0->find_matrix (R_index);
71- auto mat_nspin_3 = ap_nspin_3->find_matrix (R_index);
72-
73- // The row size and the col size of upper_matrix is double that of matrix_nspin_0
74- for (int irow = 0 ; irow < mat_nspin_0->get_row_size (); ++irow)
75- {
76- for (int icol = 0 ; icol < mat_nspin_0->get_col_size (); ++icol)
119+ hamilt::AtomPair<std::complex <double >>* upper_ap = ap;
120+ hamilt::AtomPair<std::complex <double >>* lower_ap = hRGint_tmpCd->find_pair (iat2, iat1);
121+ const hamilt::AtomPair<double >* ap_nspin1 = hRGint_tmp[first[is]].find_pair (iat1, iat2);
122+ const hamilt::AtomPair<double >* ap_nspin2 = hRGint_tmp[second[is]].find_pair (iat1, iat2);
123+ for (int ir = 0 ; ir < upper_ap->get_R_size (); ir++)
124+ {
125+ const auto R_index = upper_ap->get_R_index (ir);
126+ auto upper_mat = upper_ap->find_matrix (R_index);
127+ auto mat_nspin1 = ap_nspin1->find_matrix (R_index);
128+ auto mat_nspin2 = ap_nspin2->find_matrix (R_index);
129+ // The row size and the col size of upper_matrix is double that of matrix_nspin_0
130+ for (int irow = 0 ; irow < mat_nspin1->get_row_size (); ++irow)
77131 {
78- upper_mat->get_value (2 *irow, 2 *icol) = mat_nspin_0->get_value (irow, icol) + mat_nspin_3->get_value (irow, icol);
79- upper_mat->get_value (2 *irow+1 , 2 *icol+1 ) = mat_nspin_0->get_value (irow, icol) - mat_nspin_3->get_value (irow, icol);
132+ for (int icol = 0 ; icol < mat_nspin1->get_col_size (); ++icol)
133+ {
134+ upper_mat->get_value (irow, icol) = mat_nspin1->get_value (irow, icol)
135+ + std::complex <double >(clx_i[is], clx_j[is]) * mat_nspin2->get_value (irow, icol);
136+ }
80137 }
81- }
82-
83- if (PARAM.globalv .domag )
84- {
85- const hamilt::AtomPair<double >* ap_nspin_1 = hr_gint_part[1 ].find_pair (iat1, iat2);
86- const hamilt::AtomPair<double >* ap_nspin_2 = hr_gint_part[2 ].find_pair (iat1, iat2);
87- const auto mat_nspin_1 = ap_nspin_1->find_matrix (R_index);
88- const auto mat_nspin_2 = ap_nspin_2->find_matrix (R_index);
89- for (int irow = 0 ; irow < mat_nspin_1->get_row_size (); ++irow)
138+ // fill the lower triangle matrix
139+ // When is=0 or 3, the real part does not need conjugation;
140+ // when is=1 or 2, the small matrix is not Hermitian, so conjugation is not needed
141+ if (iat1 < iat2)
90142 {
91- for (int icol = 0 ; icol < mat_nspin_1->get_col_size (); ++icol)
143+ auto lower_mat = lower_ap->find_matrix (-R_index);
144+ for (int irow = 0 ; irow < upper_mat->get_row_size (); ++irow)
92145 {
93- upper_mat->get_value (2 *irow, 2 *icol+1 ) = mat_nspin_1->get_value (irow, icol) + std::complex <double >(0.0 , 1.0 ) * mat_nspin_2->get_value (irow, icol);
94- upper_mat->get_value (2 *irow+1 , 2 *icol) = mat_nspin_1->get_value (irow, icol) - std::complex <double >(0.0 , 1.0 ) * mat_nspin_2->get_value (irow, icol);
146+ for (int icol = 0 ; icol < upper_mat->get_col_size (); ++icol)
147+ {
148+ lower_mat->get_value (icol, irow) = upper_mat->get_value (irow, icol);
149+ }
95150 }
96151 }
152+
97153 }
154+ }
155+ }
98156
99- // fill the lower triangle matrix
100- if (iat1 < iat2)
157+ // transfer hRGint_tmpCd to parallel hR_tmp
158+ hamilt::transferSerials2Parallels ( *hRGint_tmpCd, hR_tmp);
159+ // merge hR_tmp to hR
160+ for (int iap = 0 ; iap < hR->size_atom_pairs (); iap++)
161+ {
162+ auto * ap = &hR->get_atom_pair (iap);
163+ const int iat1 = ap->get_atom_i ();
164+ const int iat2 = ap->get_atom_j ();
165+ auto * ap_nspin = hR_tmp ->find_pair (iat1, iat2);
166+ for (int ir = 0 ; ir < ap->get_R_size (); ir++)
167+ {
168+ const auto R_index = ap->get_R_index (ir);
169+ auto upper_mat = ap->find_matrix (R_index);
170+ auto mat_nspin = ap_nspin->find_matrix (R_index);
171+ // The row size and the col size of upper_matrix is double that of matrix_nspin_0
172+ for (int irow = 0 ; irow < mat_nspin->get_row_size (); ++irow)
101173 {
102- auto lower_mat = lower_ap->find_matrix (-R_index);
103- for (int irow = 0 ; irow < upper_mat->get_row_size (); ++irow)
174+ for (int icol = 0 ; icol < mat_nspin->get_col_size (); ++icol)
104175 {
105- for (int icol = 0 ; icol < upper_mat->get_col_size (); ++icol)
106- {
107- lower_mat->get_value (icol, irow) = conj (upper_mat->get_value (irow, icol));
108- }
176+ upper_mat->get_value (2 *irow+row_set[is], 2 *icol+col_set[is]) =
177+ mat_nspin->get_value (irow, icol);
109178 }
110179 }
111180 }
112181 }
113- }
114- ModuleBase::timer::tick (" Gint" , " compose_hr_gint" );
115- }
116-
117- template <typename T>
118- void transfer_hr_gint_to_hR (const HContainer<T>& hr_gint, HContainer<T>& hR)
119- {
120- ModuleBase::TITLE (" Gint" , " transfer_hr_gint_to_hR" );
121- ModuleBase::timer::tick (" Gint" , " transfer_hr_gint_to_hR" );
122- #ifdef __MPI
123- int size = 0 ;
124- MPI_Comm_size (MPI_COMM_WORLD, &size);
125- if (size == 1 )
126- {
127- hR.add (hr_gint);
128- }
129- else
130- {
131- hamilt::transferSerials2Parallels (hr_gint, &hR);
182+ delete hRGint_tmpCd;
132183 }
133184#else
134- hR. add (hr_gint);
185+
135186#endif
136- ModuleBase::timer::tick (" Gint" , " transfer_hr_gint_to_hR" );
187+ ModuleBase::timer::tick (" Gint_k" , " transfer_pvpR" );
188+ return ;
137189}
138190
191+
192+
139193// gint_info should not have been a parameter, but it was added to initialize dm_gint_full
140194// In the future, we might try to remove the gint_info parameter
141195template <typename T>
0 commit comments