@@ -78,86 +78,113 @@ 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+ int mg = hR->get_paraV ()->get_global_row_size ()/2 ;
82+ int ng = hR->get_paraV ()->get_global_col_size ()/2 ;
83+ int nb = hR->get_paraV ()->get_block_size ()/2 ;
84+ #ifdef __MPI
85+ int blacs_ctxt = hR->get_paraV ()->blacs_ctxt ;
86+ std::vector<int > iat2iwt (ucell_in->nat );
87+ for (int iat = 0 ; iat < ucell_in->nat ; iat++) {
88+ iat2iwt[iat] = ucell_in->get_iat2iwt ()[iat]/2 ;
89+ }
90+ Parallel_Orbitals *pv = new Parallel_Orbitals ();
91+ pv->set (mg, ng, nb, blacs_ctxt);
92+ pv->set_atomic_trace (iat2iwt.data (), ucell_in->nat , mg);
93+ auto ijr_info = hR->get_ijr_info ();
10094
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)
95+ this ->hR_tmp = new hamilt::HContainer<std::complex <double >>(pv, nullptr , &ijr_info);
96+ ModuleBase::Memory::record (" Gint::hRGintCd" this ->hR_tmp ->get_memory_size ());
97+
98+ // select hRGint_tmp
99+ std::vector<int > first = {0 , 1 , 1 , 0 };
100+ std::vector<int > second= {3 , 2 , 2 , 3 };
101+ // select position in the big matrix
102+ std::vector<int > row_set = {0 , 0 , 1 , 1 };
103+ std::vector<int > col_set = {0 , 1 , 0 , 1 };
104+ // construct complex matrix
105+ std::vector<int > clx_i = {1 , 0 , 0 , -1 };
106+ std::vector<int > clx_j = {0 , 1 , -1 , 0 };
107+ for (int is = 0 ; is < 4 ; is++){
108+ if (!PARAM.globalv .domag && (is==1 || is==2 )) continue ;
109+ this ->hR_tmp ->set_zero ();
110+ hamilt::HContainer<std::complex <double >>* hRGint_tmpCd = new hamilt::HContainer<std::complex <double >>(this ->ucell ->nat );
111+ hRGint_tmpCd->insert_ijrs (this ->gridt ->get_ijr_info (), *(this ->ucell ));
112+ hRGint_tmpCd->allocate (nullptr , true );
113+ hRGint_tmpCd->set_zero ();
114+ for (int iap = 0 ; iap < hRGint_tmpCd->size_atom_pairs (); iap++)
115+ {
116+ auto * ap = &hRGint_tmpCd->get_atom_pair (iap);
117+ const int iat1 = ap->get_atom_i ();
118+ const int iat2 = ap->get_atom_j ();
119+ if (iat1 <= iat2)
120+ {
121+ hamilt::AtomPair<std::complex <double >>* upper_ap = ap;
122+ hamilt::AtomPair<std::complex <double >>* lower_ap = hRGint_tmpCd->find_pair (iat2, iat1);
123+ const hamilt::AtomPair<double >* ap_nspin1 = this ->hRGint_tmp [first[is]] ->find_pair (iat1, iat2);
124+ const hamilt::AtomPair<double >* ap_nspin2 = this ->hRGint_tmp [second[is]] ->find_pair (iat1, iat2);
125+ for (int ir = 0 ; ir < upper_ap->get_R_size (); ir++)
126+ {
127+ const auto R_index = upper_ap->get_R_index (ir);
128+ auto upper_mat = upper_ap->find_matrix (R_index);
129+ auto mat_nspin1 = ap_nspin1->find_matrix (R_index);
130+ auto mat_nspin2 = ap_nspin2->find_matrix (R_index);
131+ // The row size and the col size of upper_matrix is double that of matrix_nspin_0
132+ for (int irow = 0 ; irow < mat_nspin1->get_row_size (); ++irow)
105133 {
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);
134+ for (int icol = 0 ; icol < mat_nspin1->get_col_size (); ++icol)
135+ {
136+ upper_mat->get_value (irow, icol) = mat_nspin1->get_value (irow, icol)
137+ + std::complex <double >(clx_i[is], clx_j[is]) * mat_nspin2->get_value (irow, icol);
138+ }
108139 }
109- }
110-
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)
140+ // fill the lower triangle matrix
141+ // When is=0 or 3, the real part does not need conjugation;
142+ // when is=1 or 2, the small matrix is not Hermitian, so conjugation is not needed
143+ if (iat1 < iat2)
118144 {
119- for (int icol = 0 ; icol < mat_nspin_1->get_col_size (); ++icol)
145+ auto lower_mat = lower_ap->find_matrix (-R_index);
146+ for (int irow = 0 ; irow < upper_mat->get_row_size (); ++irow)
120147 {
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);
148+ for (int icol = 0 ; icol < upper_mat->get_col_size (); ++icol)
149+ {
150+ lower_mat->get_value (icol, irow) = upper_mat->get_value (irow, icol);
151+ }
123152 }
124153 }
125- }
126154
127- // fill the lower triangle matrix
128- if (iat1 < iat2)
155+ }
156+ }
157+ }
158+ // transfer hRGint_tmpCd to parallel hR_tmp
159+ hamilt::transferSerials2Parallels ( *hRGint_tmpCd, this ->hR_tmp );
160+ // merge hR_tmp to hR
161+ for (int iap = 0 ; iap < hR->size_atom_pairs (); iap++)
162+ {
163+ auto * ap = &hR->get_atom_pair (iap);
164+ const int iat1 = ap->get_atom_i ();
165+ const int iat2 = ap->get_atom_j ();
166+ auto * ap_nspin = this ->hR_tmp ->find_pair (iat1, iat2);
167+ for (int ir = 0 ; ir < ap->get_R_size (); ir++)
168+ {
169+ const auto R_index = ap->get_R_index (ir);
170+ auto upper_mat = ap->find_matrix (R_index);
171+ auto mat_nspin = ap_nspin->find_matrix (R_index);
172+ // The row size and the col size of upper_matrix is double that of matrix_nspin_0
173+ for (int irow = 0 ; irow < mat_nspin->get_row_size (); ++irow)
129174 {
130- auto lower_mat = lower_ap->find_matrix (-R_index);
131- for (int irow = 0 ; irow < upper_mat->get_row_size (); ++irow)
175+ for (int icol = 0 ; icol < mat_nspin->get_col_size (); ++icol)
132176 {
133- for (int icol = 0 ; icol < upper_mat->get_col_size (); ++icol)
134- {
135- lower_mat->get_value (icol, irow) = conj (upper_mat->get_value (irow, icol));
136- }
177+ upper_mat->get_value (2 *irow+row_set[is], 2 *icol+col_set[is]) =
178+ mat_nspin->get_value (irow, icol);
137179 }
138180 }
139181 }
140182 }
141- }
142-
143- // ===================================
144- // transfer HR from Gint to Veff<OperatorLCAO<std::complex<double>, std::complex<double>>>
145- // ===================================
146- #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);
183+ delete hRGint_tmpCd;
156184 }
157185#else
158- hR->add (*this ->hRGintCd );
159- #endif
160186
187+ #endif
161188 ModuleBase::timer::tick (" Gint_k" " transfer_pvpR"
162189 return ;
163190}
0 commit comments