diff --git a/ABACUS.develop/source/src_global/sph_bessel_recursive.h b/ABACUS.develop/source/src_global/sph_bessel_recursive.h index 3eb51b10fc0..219a638fe0a 100644 --- a/ABACUS.develop/source/src_global/sph_bessel_recursive.h +++ b/ABACUS.develop/source/src_global/sph_bessel_recursive.h @@ -7,6 +7,7 @@ #define SPH_BESSEL_RECURSIVE_H #include +#include"stddef.h" class Sph_Bessel_Recursive { diff --git a/ABACUS.develop/source/src_pw/diago_cg.cpp b/ABACUS.develop/source/src_pw/diago_cg.cpp index ed6fb16dcc8..6d69f022ecf 100644 --- a/ABACUS.develop/source/src_pw/diago_cg.cpp +++ b/ABACUS.develop/source/src_pw/diago_cg.cpp @@ -444,7 +444,7 @@ void Diago_CG::schmit_orth //qianrui replace 2021-3-15 char trans2='N'; zgemv_(&trans2,&dim,&m,&NEG_ONE,psi.c,&dmx,lagrange,&inc,&ONE,psi_m,&inc); - psi_norm -= ddot_real(m,lagrange,lagrange); + psi_norm -= ddot_real(m,lagrange,lagrange,false); //====================================================================== /*for (int j = 0; j < m; j++) { @@ -486,23 +486,8 @@ double Diago_CG::ddot_real ( const int &dim, const complex* psi_L, - const complex* psi_R -) -{ - complex result(0,0); - for (int i=0;i Diago_CG::ddot -( - const int & dim, - const complex * psi_L, - const complex * psi_R + const complex* psi_R, + const bool reduce ) { //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< @@ -513,7 +498,7 @@ complex Diago_CG::ddot pL=(double *)psi_L; pR=(double *)psi_R; double result=LapackConnector::dot(dim2,pL,1,pR,1); - Parallel_Reduce::reduce_double_pool( result ); + if(reduce) Parallel_Reduce::reduce_double_pool( result ); return result; //====================================================================== /*complex result(0,0); @@ -524,8 +509,23 @@ complex Diago_CG::ddot Parallel_Reduce::reduce_complex_double_pool( result ); return result.real();*/ //>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> -} // end of ddot +} +complex Diago_CG::ddot +( + const int & dim, + const complex * psi_L, + const complex * psi_R +) +{ + complex result(0, 0); + for (int i = 0; i < dim ; i++) + { + result += conj(psi_L[i]) * psi_R[i] ; + } + Parallel_Reduce::reduce_complex_double_pool( result ); + return result; +} // end of ddot // this return complex Diago_CG::ddot diff --git a/ABACUS.develop/source/src_pw/diago_cg.h b/ABACUS.develop/source/src_pw/diago_cg.h index d389d904225..7065ade9dec 100644 --- a/ABACUS.develop/source/src_pw/diago_cg.h +++ b/ABACUS.develop/source/src_pw/diago_cg.h @@ -15,13 +15,15 @@ class Diago_CG static double ddot_real( const int & dim, const complex* psi_L, - const complex* psi_R) ; + const complex* psi_R, + const bool reduce = true) ; static complex ddot( const int & dim, const complex* psi_L, const complex* psi_R ) ; + static complex ddot( const int & dim, const ComplexMatrix &psi, diff --git a/ABACUS.develop/source/src_pw/hamilt_pw.cpp b/ABACUS.develop/source/src_pw/hamilt_pw.cpp index ef47f13b815..b0de0211df6 100644 --- a/ABACUS.develop/source/src_pw/hamilt_pw.cpp +++ b/ABACUS.develop/source/src_pw/hamilt_pw.cpp @@ -13,7 +13,6 @@ Hamilt_PW::Hamilt_PW() spsi = new complex[1]; GR_index = new int[1]; Bec = new complex[1]; - Ps = new complex[1]; } Hamilt_PW::~Hamilt_PW() @@ -22,7 +21,6 @@ Hamilt_PW::~Hamilt_PW() delete[] spsi; delete[] GR_index; delete[] Bec; - delete[] Ps; } @@ -43,13 +41,11 @@ void Hamilt_PW::allocate( delete[] spsi; delete[] GR_index; delete[] Bec; - delete[] Ps; this->hpsi = new complex [npwx * npol]; this->spsi = new complex [npwx * npol]; this->GR_index = new int[nrxx]; this->Bec = new complex [nkb]; - this->Ps = new complex [nkb * npol]; ZEROS(this->hpsi, npwx * npol); ZEROS(this->spsi, npwx * npol); @@ -135,16 +131,19 @@ void Hamilt_PW::cinitcgg( complex *aux=new complex [dmax*nstart]; complex *paux = aux; complex *ppsi = psi.c; - for(int m=0;mh_psi(ppsi, paux); - paux += dmax; - ppsi += dmax; - } + //qianrui replace it + this->h_psi(psi.c, aux, nstart); + //for(int m=0;mh_psi(ppsi, paux); + // paux += dmax; + // ppsi += dmax; + //} char trans1 = 'C'; char trans2 = 'N'; zgemm_(&trans1,&trans2,&nstart,&nstart,&dmin,&ONE,psi.c,&dmax,aux,&dmax,&ZERO,hc.c,&nstart); hc=transpose(hc,false); + zgemm_(&trans1,&trans2,&nstart,&nstart,&dmin,&ONE,psi.c,&dmax,psi.c,&dmax,&ZERO,sc.c,&nstart); sc=transpose(sc,false); //After psis are strictly normalized, we should use this part. @@ -462,103 +461,157 @@ void Hamilt_PW::s_1psi } -void Hamilt_PW::h_psi(const complex *psi_in, complex *hpsi) +void Hamilt_PW::h_psi(const complex *psi_in, complex *hpsi, const int m) { timer::tick("Hamilt_PW","h_psi",'H'); int i = 0; int j = 0; int ig= 0; - if(NSPIN!=4) ZEROS(hpsi, wf.npw); - else ZEROS(hpsi, wf.npwx * NPOL);//added by zhengdy-soc + //if(NSPIN!=4) ZEROS(hpsi, wf.npw); + //else ZEROS(hpsi, wf.npwx * NPOL);//added by zhengdy-soc + int dmax = wf.npwx * NPOL; //------------------------------------ //(1) the kinetical energy. //------------------------------------ + complex *tmhpsi; + const complex *tmpsi_in; if(T_IN_H) { - for (ig = 0;ig < wf.npw;ig++) + tmhpsi = hpsi; + tmpsi_in = psi_in; + for(int ib = 0 ; ib < m; ++ib) { - hpsi[ig] = wf.g2kin[ig] * psi_in[ig]; - } - //added by zhengdy-soc - if(NSPIN==4) - { - for (ig = wf.npwx;ig < wf.npw + wf.npwx;ig++) + for(ig = 0;ig < wf.npw; ++ig) { - hpsi[ig] = wf.g2kin[ig - wf.npwx] * psi_in[ig]; + tmhpsi[ig] = wf.g2kin[ig] * tmpsi_in[ig]; } + if(NSPIN==4){ + for(ig=wf.npw; ig < wf.npwx; ++ig) + { + tmhpsi[ig] = 0; + } + tmhpsi +=wf.npwx; + tmpsi_in += wf.npwx; + for (ig = 0;ig < wf.npw ;++ig) + { + tmhpsi[ig] = wf.g2kin[ig] * tmpsi_in[ig]; + } + for(ig=wf.npw; ig < wf.npwx; ++ig) + { + tmhpsi[ig] =0; + } + } + tmhpsi += wf.npwx; + tmpsi_in += wf.npwx; } } //------------------------------------ //(2) the local potential. - //------------------------------------ + //----------------------------------- + timer::tick("Hamilt_PW","vloc",'H'); if(VL_IN_H) { - if(NSPIN!=4) + tmhpsi = hpsi; + tmpsi_in = psi_in; + for(int ib = 0 ; ib < m; ++ib) { - ZEROS( UFFT.porter, pw.nrxx); - UFFT.RoundTrip( psi_in, pot.vr_eff1, GR_index, UFFT.porter ); - - for (j = 0;j < wf.npw;j++) - { - hpsi[j] += UFFT.porter[ GR_index[j] ]; - } - } - else - { - complex* porter1 = new complex[pw.nrxx]; - ZEROS( UFFT.porter, pw.nrxx); - ZEROS( porter1, pw.nrxx); - for (int ig=0; ig< wf.npw; ig++) - { - UFFT.porter[ GR_index[ig] ] = psi_in[ig]; - porter1[ GR_index[ig] ] = psi_in[ig + wf.npwx]; + if(NSPIN!=4){ + ZEROS( UFFT.porter, pw.nrxx); + UFFT.RoundTrip( tmpsi_in, pot.vr_eff1, GR_index, UFFT.porter ); + for (j = 0;j < wf.npw;j++) + { + tmhpsi[j] += UFFT.porter[ GR_index[j] ]; + } } - // (2) fft to real space and doing things. - pw.FFT_wfc.FFT3D( UFFT.porter, 1); - pw.FFT_wfc.FFT3D( porter1, 1); - complex sup,sdown; - for (int ir=0; ir< pw.nrxx; ir++) + else { - sup = UFFT.porter[ir] * (pot.vr_eff(0,ir) + pot.vr_eff(3,ir)) + - porter1[ir] * (pot.vr_eff(1,ir) - complex(0.0,1.0) * pot.vr_eff(2,ir)); - sdown = porter1[ir] * (pot.vr_eff(0,ir) - pot.vr_eff(3,ir)) + - UFFT.porter[ir] * (pot.vr_eff(1,ir) + complex(0.0,1.0) * pot.vr_eff(2,ir)); - UFFT.porter[ir] = sup; - porter1[ir] = sdown; - } - // (3) fft back to G space. - pw.FFT_wfc.FFT3D( UFFT.porter, -1); - pw.FFT_wfc.FFT3D( porter1, -1); + complex* porter1 = new complex[pw.nrxx]; + ZEROS( UFFT.porter, pw.nrxx); + ZEROS( porter1, pw.nrxx); + for (int ig=0; ig< wf.npw; ig++) + { + UFFT.porter[ GR_index[ig] ] = tmpsi_in[ig]; + porter1[ GR_index[ig] ] = tmpsi_in[ig + wf.npwx]; + } + // (2) fft to real space and doing things. + pw.FFT_wfc.FFT3D( UFFT.porter, 1); + pw.FFT_wfc.FFT3D( porter1, 1); + complex sup,sdown; + for (int ir=0; ir< pw.nrxx; ir++) + { + sup = UFFT.porter[ir] * (pot.vr_eff(0,ir) + pot.vr_eff(3,ir)) + + porter1[ir] * (pot.vr_eff(1,ir) - complex(0.0,1.0) * pot.vr_eff(2,ir)); + sdown = porter1[ir] * (pot.vr_eff(0,ir) - pot.vr_eff(3,ir)) + + UFFT.porter[ir] * (pot.vr_eff(1,ir) + complex(0.0,1.0) * pot.vr_eff(2,ir)); + UFFT.porter[ir] = sup; + porter1[ir] = sdown; + } + // (3) fft back to G space. + pw.FFT_wfc.FFT3D( UFFT.porter, -1); + pw.FFT_wfc.FFT3D( porter1, -1); - for (j = 0;j < wf.npw;j++) - { - hpsi[j] += UFFT.porter[ GR_index[j] ]; - hpsi[j+wf.npwx] += porter1[ GR_index[j] ]; + for (j = 0;j < wf.npw;j++) + { + tmhpsi[j] += UFFT.porter[ GR_index[j] ]; + } + for (j = 0;j < wf.npw;j++ ) + { + tmhpsi[j+wf.npwx] += porter1[ GR_index[j] ]; + } + delete[] porter1; } - delete[] porter1; + tmhpsi += dmax; + tmpsi_in += dmax; } } + timer::tick("Hamilt_PW","vloc",'H'); //------------------------------------ // (3) the nonlocal pseudopotential. //------------------------------------ + timer::tick("Hamilt_PW","vnl",'H'); if(VNL_IN_H) { if ( ppcell.nkb > 0) { //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< - //qianrui improve 2021-3-16 + //qianrui optimize 2021-3-31 int nkb=ppcell.nkb; - ComplexMatrix becp(NPOL,nkb,false); + ComplexMatrix becp(NPOL * m, nkb, false); char transa = 'C'; char transb = 'N'; - zgemm_(&transa,&transb,&nkb,&NPOL,&wf.npw,&ONE,ppcell.vkb.c,&wf.npwx,psi_in,&wf.npwx,&ZERO,becp.c,&nkb); - becp=transpose(becp,false); - Parallel_Reduce::reduce_complex_double_pool( becp.c, ppcell.nkb * NPOL); - this->add_vuspsi(hpsi, becp.c); + if(m==1 && NPOL==1) + { + int inc = 1; + zgemv_(&transa, &wf.npw, &nkb, &ONE, ppcell.vkb.c, &wf.npwx, psi_in, &inc, &ZERO, becp.c, &inc); + } + else + { + int npm = NPOL * m; + zgemm_(&transa,&transb,&nkb,&npm,&wf.npw,&ONE,ppcell.vkb.c,&wf.npwx,psi_in,&wf.npwx,&ZERO,becp.c,&nkb); + //add_vuspsi is moddified, thus tranpose not needed here. + //if(NONCOLIN) + //{ + // ComplexMatrix partbecp(NPOL, nkb ,false); + // for(int ib = 0; ib < m; ++ib) + // { +// + // for ( i = 0;i < NPOL;i++) + // for (j = 0;j < nkb;j++) + // partbecp(i, j) = tmbecp[i*nkb+j]; + // for (j = 0; j < nkb; j++) + // for (i = 0;i < NPOL;i++) + // tmbecp[j*NPOL+i] = partbecp(i, j); + // tmbecp += NPOL * nkb; + // } + //} + } + + Parallel_Reduce::reduce_complex_double_pool( becp.c, nkb * NPOL * m); + this->add_vuspsi(hpsi, becp.c, m); //====================================================================== /*complex *becp = new complex[ ppcell.nkb * NPOL ]; ZEROS(becp,ppcell.nkb * NPOL); @@ -582,21 +635,25 @@ void Hamilt_PW::h_psi(const complex *psi_in, complex *hpsi) //>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> } } - + timer::tick("Hamilt_PW","vnl",'H'); timer::tick("Hamilt_PW","h_psi",'H'); return; } -void Hamilt_PW::add_vuspsi(complex *hpsi_in,const complex *becp) +void Hamilt_PW::add_vuspsi(complex *hpsi_in,const complex *becp, const int m) { timer::tick("Hamilt_PW","add_vuspsi",'I'); - ZEROS( Ps, ppcell.nkb * NPOL ); + int nkb = ppcell.nkb; + complex *Ps = new complex [nkb * NPOL * m]; + ZEROS( Ps, NPOL * m * nkb); int sum = 0; int iat = 0; // this function sum up each non-local pseudopotential located in each atom, // all we need to do is put the right Dij coefficient to each becp, which // is calculated before. + if(NSPIN!=4) + { for (int it=0; it *hpsi_in,const complex *becp) { for (int ip2=0; ip2Ps[sum+ip2] += ppcell.deeq(CURRENT_SPIN, iat, ip, ip2) * becp[sum+ip]; - else + for(int ib = 0; ib < m ; ++ib) { - this->Ps[sum+ ip2*2] += ppcell.deeq_nc(0, iat, ip2, ip) * becp[sum+ip*2] - +ppcell.deeq_nc(1, iat, ip2, ip) * becp[sum+ip*2+1]; - this->Ps[sum+ ip2*2+1] += ppcell.deeq_nc(2, iat, ip2, ip) * becp[sum+ip*2] - +ppcell.deeq_nc(3, iat, ip2, ip) * becp[sum+ip*2+1]; - } - }// end ih + Ps[(sum + ip2) * m + ib] += ppcell.deeq(CURRENT_SPIN, iat, ip, ip2) * becp[ib * nkb + sum + ip]; + }//end ib + }// end ih + }//end jh + sum += Nprojs; + ++iat; + } //end na + } //end nt + } + else + { + for (int it=0; it becp1,becp2; + const int Nprojs = ucell.atoms[it].nh; + for (int ia=0; iaPs[i]*ppcell.vkb(i,ig); + hpsi_in[ig]+=Ps[i]*ppcell.vkb(i,ig); } } */ @@ -639,10 +725,19 @@ void Hamilt_PW::add_vuspsi(complex *hpsi_in,const complex *becp) // use simple method. //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< - //qianrui improve 2021-3-16 + //qianrui optimize 2021-3-31 char transa = 'N'; char transb = 'T'; - zgemm_(&transa,&transb,&wf.npw,&NPOL,&ppcell.nkb,&ONE,ppcell.vkb.c,&wf.npwx,Ps,&NPOL,&ONE,hpsi_in,&wf.npwx); + if(NPOL==1 && m==1) + { + int inc = 1; + zgemv_(&transa, &wf.npw, &ppcell.nkb, &ONE, ppcell.vkb.c, &wf.npwx, Ps, &inc, &ONE, hpsi_in, &inc); + } + else + { + int npm = NPOL*m; + zgemm_(&transa,&transb,&wf.npw,&npm,&ppcell.nkb,&ONE,ppcell.vkb.c,&wf.npwx,Ps,&npm,&ONE,hpsi_in,&wf.npwx); + } //====================================================================== /*if(!NONCOLIN) for(int i=0; i *hpsi_in,const complex *becp) } }*/ //>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> - + delete []Ps; timer::tick("Hamilt_PW","add_vuspsi",'I'); return; } diff --git a/ABACUS.develop/source/src_pw/hamilt_pw.h b/ABACUS.develop/source/src_pw/hamilt_pw.h index 6113c5ae32e..b79d094e2a4 100644 --- a/ABACUS.develop/source/src_pw/hamilt_pw.h +++ b/ABACUS.develop/source/src_pw/hamilt_pw.h @@ -50,7 +50,7 @@ class Hamilt_PW complex *hpsi, complex *spsi); - void h_psi( const complex *psi, complex *hpsi); + void h_psi( const complex *psi, complex *hpsi, const int m = 1); // qianrui add a default parameter 2021-3-31 void s_1psi( const int npw, @@ -66,11 +66,9 @@ class Hamilt_PW // hpsi , spsi complex *hpsi; complex *spsi; - complex *Bec; - complex *Ps; - void add_vuspsi(complex *hpsi, const complex *becp); + void add_vuspsi(complex *hpsi, const complex *becp, const int m); private: