@@ -457,6 +457,11 @@ void DensityMatrix<std::complex<double>, double>::cal_DMR()
457457 {
458458 throw std::string (" Atom-pair not belong this process" );
459459 }
460+ std::vector<std::complex <double >> tmp_DMR;
461+ if (PARAM.inp .nspin == 4 )
462+ {
463+ tmp_DMR.resize (tmp_ap.get_size ());
464+ }
460465 for (int ir = 0 ; ir < tmp_ap.get_R_size (); ++ir)
461466 {
462467 const ModuleBase::Vector3<int > r_index = tmp_ap.get_R_index (ir);
@@ -514,9 +519,7 @@ void DensityMatrix<std::complex<double>, double>::cal_DMR()
514519 // treat DMR as pauli matrix when NSPIN=4
515520 if (PARAM.inp .nspin == 4 )
516521 {
517- std::vector<std::complex <double >> tmp_DMR (this ->_paraV ->get_col_size ()
518- * this ->_paraV ->get_row_size (),
519- std::complex <double >(0.0 , 0.0 ));
522+ tmp_DMR.assign (tmp_ap.get_size (), std::complex <double >(0.0 , 0.0 ));
520523 for (int ik = 0 ; ik < this ->_nk ; ++ik)
521524 {
522525 // cal k_phase
@@ -534,35 +537,34 @@ void DensityMatrix<std::complex<double>, double>::cal_DMR()
534537 // jump DMK to fill DMR
535538 // DMR is row-major, DMK is column-major
536539 tmp_DMK_pointer += col_ap * this ->_paraV ->nrow + row_ap;
537- for (int mu = 0 ; mu < this -> _paraV -> get_row_size (iat1 ); ++mu)
540+ for (int mu = 0 ; mu < tmp_ap. get_row_size (); ++mu)
538541 {
539- BlasConnector::axpy (this -> _paraV -> get_col_size (iat2 ),
542+ BlasConnector::axpy (tmp_ap. get_col_size (),
540543 kphase,
541544 tmp_DMK_pointer,
542545 ld_hk,
543546 tmp_DMR_pointer,
544547 1 );
545548 tmp_DMK_pointer += 1 ;
546- tmp_DMR_pointer += this -> _paraV -> get_col_size (iat2 );
549+ tmp_DMR_pointer += tmp_ap. get_col_size ();
547550 }
548551 }
549552 int npol = 2 ;
550553 // step_trace = 0 for NSPIN=1,2; ={0, 1, local_col, local_col+1} for NSPIN=4
551- std::vector< int > step_trace (npol * npol, 0 ) ;
554+ int step_trace[ 4 ] ;
552555 for (int is = 0 ; is < npol; is++)
553556 {
554557 for (int is2 = 0 ; is2 < npol; is2++)
555558 {
556- step_trace[is * npol + is2] = this ->_paraV ->get_col_size (iat2) * is + is2;
557- // step_trace[is + is2 * npol] = this->_paraV->get_col_size(iat2) * is + is2;
559+ step_trace[is * npol + is2] = tmp_ap.get_col_size () * is + is2;
558560 }
559561 }
560562 std::complex <double > tmp[4 ];
561563 double * target_DMR = tmp_matrix->get_pointer ();
562564 std::complex <double >* tmp_DMR_pointer = tmp_DMR.data ();
563- for (int irow = 0 ; irow < this -> _paraV -> get_row_size (iat1 ); irow += 2 )
565+ for (int irow = 0 ; irow < tmp_ap. get_row_size (); irow += 2 )
564566 {
565- for (int icol = 0 ; icol < this -> _paraV -> get_col_size (iat2 ); icol += 2 )
567+ for (int icol = 0 ; icol < tmp_ap. get_col_size (); icol += 2 )
566568 {
567569 // catch the 4 spin component value of one orbital pair
568570 tmp[0 ] = tmp_DMR_pointer[icol + step_trace[0 ]];
@@ -577,8 +579,8 @@ void DensityMatrix<std::complex<double>, double>::cal_DMR()
577579 = -tmp[1 ].imag () + tmp[2 ].imag (); // (i * (rho_updown - rho_downup)).real()
578580 target_DMR[icol + step_trace[3 ]] = tmp[0 ].real () - tmp[3 ].real ();
579581 }
580- tmp_DMR_pointer += this -> _paraV -> get_col_size (iat2 ) * 2 ;
581- target_DMR += this -> _paraV -> get_col_size (iat2 ) * 2 ;
582+ tmp_DMR_pointer += tmp_ap. get_col_size () * 2 ;
583+ target_DMR += tmp_ap. get_col_size () * 2 ;
582584 }
583585 }
584586 }
0 commit comments