@@ -78,86 +78,116 @@ void Gint_k::transfer_pvpR(hamilt::HContainer<std::complex<double>>* hR,
7878 ModuleBase::TITLE (" Gint_k" " transfer_pvpR"
7979 ModuleBase::timer::tick (" Gint_k" " transfer_pvpR"
8080
81- this ->hRGintCd ->set_zero ();
82-
83- for (int iap = 0 ; iap < this ->hRGintCd ->size_atom_pairs (); iap++)
84- {
85- auto * ap = &this ->hRGintCd ->get_atom_pair (iap);
86- const int iat1 = ap->get_atom_i ();
87- const int iat2 = ap->get_atom_j ();
88- if (iat1 <= iat2)
89- {
90- hamilt::AtomPair<std::complex <double >>* upper_ap = ap;
91- hamilt::AtomPair<std::complex <double >>* lower_ap = this ->hRGintCd ->find_pair (iat2, iat1);
92- const hamilt::AtomPair<double >* ap_nspin_0 = this ->hRGint_tmp [0 ]->find_pair (iat1, iat2);
93- const hamilt::AtomPair<double >* ap_nspin_3 = this ->hRGint_tmp [3 ]->find_pair (iat1, iat2);
94- for (int ir = 0 ; ir < upper_ap->get_R_size (); ir++)
95- {
96- const auto R_index = upper_ap->get_R_index (ir);
97- auto upper_mat = upper_ap->find_matrix (R_index);
98- auto mat_nspin_0 = ap_nspin_0->find_matrix (R_index);
99- auto mat_nspin_3 = ap_nspin_3->find_matrix (R_index);
81+ auto ijr_info = hR->get_ijr_info ();
10082
101- // The row size and the col size of upper_matrix is double that of matrix_nspin_0
102- for (int irow = 0 ; irow < mat_nspin_0->get_row_size (); ++irow)
103- {
104- for (int icol = 0 ; icol < mat_nspin_0->get_col_size (); ++icol)
105- {
106- upper_mat->get_value (2 *irow, 2 *icol) = mat_nspin_0->get_value (irow, icol) + mat_nspin_3->get_value (irow, icol);
107- upper_mat->get_value (2 *irow+1 , 2 *icol+1 ) = mat_nspin_0->get_value (irow, icol) - mat_nspin_3->get_value (irow, icol);
108- }
109- }
83+ #ifdef __MPI
84+ int mg = hR->get_paraV ()->get_global_row_size ()/2 ;
85+ int ng = hR->get_paraV ()->get_global_col_size ()/2 ;
86+ int nb = hR->get_paraV ()->get_block_size ()/2 ;
87+ int blacs_ctxt = hR->get_paraV ()->blacs_ctxt ;
88+ std::vector<int > iat2iwt (ucell_in->nat );
89+ for (int iat = 0 ; iat < ucell_in->nat ; iat++) {
90+ iat2iwt[iat] = ucell_in->get_iat2iwt ()[iat]/2 ;
91+ }
92+ Parallel_Orbitals *pv = new Parallel_Orbitals ();
93+ pv->set (mg, ng, nb, blacs_ctxt);
94+ pv->set_atomic_trace (iat2iwt.data (), ucell_in->nat , mg);
95+ this ->hR_tmp = new hamilt::HContainer<std::complex <double >>(pv, nullptr , &ijr_info);
96+ #endif
97+
98+ ModuleBase::Memory::record (" Gint::hRGintCd" this ->hR_tmp ->get_memory_size ());
11099
111- if (PARAM.globalv .domag )
112- {
113- const hamilt::AtomPair<double >* ap_nspin_1 = this ->hRGint_tmp [1 ]->find_pair (iat1, iat2);
114- const hamilt::AtomPair<double >* ap_nspin_2 = this ->hRGint_tmp [2 ]->find_pair (iat1, iat2);
115- const auto mat_nspin_1 = ap_nspin_1->find_matrix (R_index);
116- const auto mat_nspin_2 = ap_nspin_2->find_matrix (R_index);
117- for (int irow = 0 ; irow < mat_nspin_1->get_row_size (); ++irow)
100+ // select hRGint_tmp
101+ std::vector<int > first = {0 , 1 , 1 , 0 };
102+ std::vector<int > second= {3 , 2 , 2 , 3 };
103+ // select position in the big matrix
104+ std::vector<int > row_set = {0 , 0 , 1 , 1 };
105+ std::vector<int > col_set = {0 , 1 , 0 , 1 };
106+ // construct complex matrix
107+ std::vector<int > clx_i = {1 , 0 , 0 , -1 };
108+ std::vector<int > clx_j = {0 , 1 , -1 , 0 };
109+ for (int is = 0 ; is < 4 ; is++){
110+ if (!PARAM.globalv .domag && (is==1 || is==2 )) continue ;
111+ this ->hR_tmp ->set_zero ();
112+ hamilt::HContainer<std::complex <double >>* hRGint_tmpCd = new hamilt::HContainer<std::complex <double >>(this ->ucell ->nat );
113+ hRGint_tmpCd->insert_ijrs (this ->gridt ->get_ijr_info (), *(this ->ucell ));
114+ hRGint_tmpCd->allocate (nullptr , true );
115+ hRGint_tmpCd->set_zero ();
116+ for (int iap = 0 ; iap < hRGint_tmpCd->size_atom_pairs (); iap++)
117+ {
118+ auto * ap = &hRGint_tmpCd->get_atom_pair (iap);
119+ const int iat1 = ap->get_atom_i ();
120+ const int iat2 = ap->get_atom_j ();
121+ if (iat1 <= iat2)
122+ {
123+ hamilt::AtomPair<std::complex <double >>* upper_ap = ap;
124+ hamilt::AtomPair<std::complex <double >>* lower_ap = hRGint_tmpCd->find_pair (iat2, iat1);
125+ const hamilt::AtomPair<double >* ap_nspin1 = this ->hRGint_tmp [first[is]] ->find_pair (iat1, iat2);
126+ const hamilt::AtomPair<double >* ap_nspin2 = this ->hRGint_tmp [second[is]] ->find_pair (iat1, iat2);
127+ for (int ir = 0 ; ir < upper_ap->get_R_size (); ir++)
128+ {
129+ const auto R_index = upper_ap->get_R_index (ir);
130+ auto upper_mat = upper_ap->find_matrix (R_index);
131+ auto mat_nspin1 = ap_nspin1->find_matrix (R_index);
132+ auto mat_nspin2 = ap_nspin2->find_matrix (R_index);
133+ // The row size and the col size of upper_matrix is double that of matrix_nspin_0
134+ for (int irow = 0 ; irow < mat_nspin1->get_row_size (); ++irow)
118135 {
119- for (int icol = 0 ; icol < mat_nspin_1 ->get_col_size (); ++icol)
136+ for (int icol = 0 ; icol < mat_nspin1 ->get_col_size (); ++icol)
120137 {
121- 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);
122- 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);
138+ upper_mat->get_value (irow, icol) = mat_nspin1 ->get_value (irow, icol)
139+ + std::complex <double >(clx_i[is], clx_j[is] ) * mat_nspin2 ->get_value (irow, icol);
123140 }
124141 }
125- }
126-
127- // fill the lower triangle matrix
128- if (iat1 < iat2)
129- {
130- auto lower_mat = lower_ap->find_matrix (-R_index);
131- for (int irow = 0 ; irow < upper_mat->get_row_size (); ++irow)
142+ // fill the lower triangle matrix
143+ // When is=0 or 3, the real part does not need conjugation;
144+ // when is=1 or 2, the small matrix is not Hermitian, so conjugation is not needed
145+ if (iat1 < iat2)
132146 {
133- for (int icol = 0 ; icol < upper_mat->get_col_size (); ++icol)
147+ auto lower_mat = lower_ap->find_matrix (-R_index);
148+ for (int irow = 0 ; irow < upper_mat->get_row_size (); ++irow)
134149 {
135- lower_mat->get_value (icol, irow) = conj (upper_mat->get_value (irow, icol));
150+ for (int icol = 0 ; icol < upper_mat->get_col_size (); ++icol)
151+ {
152+ lower_mat->get_value (icol, irow) = upper_mat->get_value (irow, icol);
153+ }
136154 }
137155 }
156+
138157 }
139158 }
140159 }
141- }
142-
143- // ===================================
144- // transfer HR from Gint to Veff<OperatorLCAO<std::complex<double>, std::complex<double>>>
145- // ===================================
146160#ifdef __MPI
147- int size;
148- MPI_Comm_size (MPI_COMM_WORLD, &size);
149- if (size == 1 )
150- {
151- hR->add (*this ->hRGintCd );
152- }
153- else
154- {
155- hamilt::transferSerials2Parallels<std::complex <double >>(*this ->hRGintCd , hR);
156- }
161+ // transfer hRGint_tmpCd to parallel hR_tmp
162+ hamilt::transferSerials2Parallels ( *hRGint_tmpCd, this ->hR_tmp );
157163#else
158- hR-> add (* this ->hRGintCd ) ;
164+ this ->hR_tmp = hRGint_tmpCd ;
159165#endif
160-
166+ // merge hR_tmp to hR
167+ for (int iap = 0 ; iap < hR->size_atom_pairs (); iap++)
168+ {
169+ auto * ap = &hR->get_atom_pair (iap);
170+ const int iat1 = ap->get_atom_i ();
171+ const int iat2 = ap->get_atom_j ();
172+ auto * ap_nspin = this ->hR_tmp ->find_pair (iat1, iat2);
173+ for (int ir = 0 ; ir < ap->get_R_size (); ir++)
174+ {
175+ const auto R_index = ap->get_R_index (ir);
176+ auto upper_mat = ap->find_matrix (R_index);
177+ auto mat_nspin = ap_nspin->find_matrix (R_index);
178+ // The row size and the col size of upper_matrix is double that of matrix_nspin_0
179+ for (int irow = 0 ; irow < mat_nspin->get_row_size (); ++irow)
180+ {
181+ for (int icol = 0 ; icol < mat_nspin->get_col_size (); ++icol)
182+ {
183+ upper_mat->get_value (2 *irow+row_set[is], 2 *icol+col_set[is]) =
184+ mat_nspin->get_value (irow, icol);
185+ }
186+ }
187+ }
188+ }
189+ delete hRGint_tmpCd;
190+ }
161191 ModuleBase::timer::tick (" Gint_k" " transfer_pvpR"
162192 return ;
163193}
0 commit comments