From deae88d7b48b66f133325444374af99b2c39640e Mon Sep 17 00:00:00 2001 From: haozhihan Date: Sat, 30 Nov 2024 11:37:23 +0800 Subject: [PATCH] format wf_atomic files --- source/module_psi/wavefunc.cpp | 2 +- source/module_psi/wf_atomic.cpp | 781 +++++++++++++++++--------------- source/module_psi/wf_atomic.h | 26 +- 3 files changed, 434 insertions(+), 375 deletions(-) diff --git a/source/module_psi/wavefunc.cpp b/source/module_psi/wavefunc.cpp index 63a168b05a..45a6e7ea18 100644 --- a/source/module_psi/wavefunc.cpp +++ b/source/module_psi/wavefunc.cpp @@ -353,7 +353,7 @@ void diago_PAO_in_pw_k2(const base_device::DEVICE_GPU* ctx, const int& lmaxkb, hamilt::Hamilt, base_device::DEVICE_GPU>* phm_in) { - // TODO float + // TODO float } template <> diff --git a/source/module_psi/wf_atomic.cpp b/source/module_psi/wf_atomic.cpp index c5d2acf21d..fbfc9837b6 100644 --- a/source/module_psi/wf_atomic.cpp +++ b/source/module_psi/wf_atomic.cpp @@ -1,13 +1,14 @@ #include "wf_atomic.h" -#include "module_hamilt_pw/hamilt_pwdft/global.h" -#include "module_parameter/parameter.h" + #include "module_base/math_integral.h" -#include "module_base/math_sphbes.h" #include "module_base/math_polyint.h" +#include "module_base/math_sphbes.h" #include "module_base/math_ylmreal.h" #include "module_base/timer.h" #include "module_base/tool_quit.h" +#include "module_hamilt_pw/hamilt_pwdft/global.h" #include "module_hamilt_pw/hamilt_pwdft/soc.h" +#include "module_parameter/parameter.h" #include @@ -17,84 +18,91 @@ WF_atomic::WF_atomic() WF_atomic::~WF_atomic() { - if(PARAM.inp.test_deconstructor) - { - std::cout << " ~WF_atomic()" << std::endl; - } - if(this->wanf2!= nullptr) + if (PARAM.inp.test_deconstructor) + { + std::cout << " ~WF_atomic()" << std::endl; + } + if (this->wanf2 != nullptr) { delete[] wanf2; } - if(this->psi != nullptr) + if (this->psi != nullptr) { delete psi; } } -void WF_atomic::init_at_1(Structure_Factor *sf_in, ModuleBase::realArray* tab_at) +void WF_atomic::init_at_1(Structure_Factor* sf_in, ModuleBase::realArray* tab_at) { - if(PARAM.inp.use_paw) { return; -} - if (PARAM.inp.test_wf) { ModuleBase::TITLE("WF_atomic","init_at_1"); -} - ModuleBase::timer::tick("WF_atomic","init_at_1"); + if (PARAM.inp.use_paw) + { + return; + } + if (PARAM.inp.test_wf) + { + ModuleBase::TITLE("WF_atomic", "init_at_1"); + } + ModuleBase::timer::tick("WF_atomic", "init_at_1"); this->psf = sf_in; GlobalV::ofs_running << "\n Make real space PAO into reciprocal space." << std::endl; this->print_PAOs(); -//---------------------------------------------------------- -// EXPLAIN : Find the type of atom that has most mesh points. -//---------------------------------------------------------- + //---------------------------------------------------------- + // EXPLAIN : Find the type of atom that has most mesh points. + //---------------------------------------------------------- int ndm = 0; - for (int it=0; it ndm) ? GlobalC::ucell.atoms[it].ncpp.msh : ndm; } - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"max mesh points in Pseudopotential",ndm); + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "max mesh points in Pseudopotential", ndm); // needed to normalize atomic wfcs (not a bad idea in general and // necessary to compute correctly lda+U projections) tab_at->zero_out(); -//---------------------------------------------------------- -// EXPLAIN : If use gauss orbitals to represent aotmic -// orbitals (controlled by parameters) -// -// USE GLOBAL CLASS VARIABLES : -// NAME : GlobalC::ucell.atoms.nchi -// NAME : GlobalC::ucell.atoms.msh(number of mesh points) -// NAME : GlobalC::ucell.atoms.r -//---------------------------------------------------------- + //---------------------------------------------------------- + // EXPLAIN : If use gauss orbitals to represent aotmic + // orbitals (controlled by parameters) + // + // USE GLOBAL CLASS VARIABLES : + // NAME : GlobalC::ucell.atoms.nchi + // NAME : GlobalC::ucell.atoms.msh(number of mesh points) + // NAME : GlobalC::ucell.atoms.r + //---------------------------------------------------------- const int startq = 0; const double pref = ModuleBase::FOUR_PI / sqrt(GlobalC::ucell.omega); - double *aux = new double[ndm]; - double *vchi = new double[ndm]; + double* aux = new double[ndm]; + double* vchi = new double[ndm]; - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"dq(describe PAO in reciprocal space)",PARAM.globalv.dq); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"max q",PARAM.globalv.nqx); + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "dq(describe PAO in reciprocal space)", PARAM.globalv.dq); + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "max q", PARAM.globalv.nqx); - for (int it=0; itlabel << " is " << atom->ncpp.nchi << std::endl; + GlobalV::ofs_running << "\n number of pseudo atomic orbitals for " << atom->label << " is " << atom->ncpp.nchi + << std::endl; - for (int ic=0; icncpp.nchi ;ic++) + for (int ic = 0; ic < atom->ncpp.nchi; ic++) { - //std::cout << "\n T=" << it << " ic=" << ic << std::endl; + // std::cout << "\n T=" << it << " ic=" << ic << std::endl; int nmesh; - if(PARAM.inp.pseudo_mesh) { + if (PARAM.inp.pseudo_mesh) + { nmesh = atom->ncpp.mesh; - } else { + } + else + { nmesh = atom->ncpp.msh; -} + } // check the unit condition - double *inner_part = new double[nmesh]; - for (int ir=0; irncpp.chi(ic,ir) * atom->ncpp.chi(ic,ir); + inner_part[ir] = atom->ncpp.chi(ic, ir) * atom->ncpp.chi(ic, ir); } double unit = 0.0; ModuleBase::Integral::Simpson_Integral(nmesh, inner_part, atom->ncpp.rab.data(), unit); @@ -182,63 +190,66 @@ void WF_atomic::init_at_1(Structure_Factor *sf_in, ModuleBase::realArray* tab_at if (atom->ncpp.oc[ic] >= 0.0) { const int l = atom->ncpp.lchi[ic]; - for (int iq=startq; iqncpp.msh, atom->ncpp.r.data(), q, l, aux); - for (int ir = 0;ir < atom->ncpp.msh;ir++) + for (int ir = 0; ir < atom->ncpp.msh; ir++) { - vchi[ir] = atom->ncpp.chi(ic,ir) * aux[ir] * atom->ncpp.r[ir]; + vchi[ir] = atom->ncpp.chi(ic, ir) * aux[ir] * atom->ncpp.r[ir]; } double vqint = 0.0; ModuleBase::Integral::Simpson_Integral(atom->ncpp.msh, vchi, atom->ncpp.rab.data(), vqint); - tab_at->operator()(it, ic, iq) = vqint * pref; + tab_at->operator()(it, ic, iq) = vqint * pref; // if( it == 0 && ic == 0 ) // { // // for (ir = 0;ir < GlobalC::ucell.atoms[it].msh;ir++) - // GlobalV::ofs_running << std::setprecision(20) << "\n vchi(" << ir << ")=" << vchi[ir]; - // GlobalV::ofs_running<<"\n aux[0] = "<3) { ModuleBase::TITLE("WF_atomic","atomic_wfc"); -} - ModuleBase::timer::tick("WF_atomic","atomic_wfc"); + if (PARAM.inp.test_wf > 3) + { + ModuleBase::TITLE("WF_atomic", "atomic_wfc"); + } + ModuleBase::timer::tick("WF_atomic", "atomic_wfc"); //========================================================= // This routine computes the superposition of atomic // WF_atomictions for a given k-point. //========================================================= const int total_lm = (lmax_wfc + 1) * (lmax_wfc + 1); ModuleBase::matrix ylm(total_lm, np); - std::complex *aux = new std::complex[np]; - double *chiaux = nullptr; + std::complex* aux = new std::complex[np]; + double* chiaux = nullptr; - ModuleBase::Vector3 *gk = new ModuleBase::Vector3 [np]; - for (int ig=0;ig* gk = new ModuleBase::Vector3[np]; + for (int ig = 0; ig < np; ig++) { - gk[ig] = wfc_basis->getgpluskcar(ik,ig); + gk[ig] = wfc_basis->getgpluskcar(ik, ig); } - //ylm = spherical harmonics functions + // ylm = spherical harmonics functions ModuleBase::YlmReal::Ylm_Real(total_lm, np, gk, ylm); int index = 0; //--------------------------------------------------------- // Calculate G space 3D wave functions //--------------------------------------------------------- - double *flq = new double[np]; - for (int it = 0;it < GlobalC::ucell.ntype;it++) + double* flq = new double[np]; + for (int it = 0; it < GlobalC::ucell.ntype; it++) { - for (int ia = 0;ia < GlobalC::ucell.atoms[it].na;ia++) + for (int ia = 0; ia < GlobalC::ucell.atoms[it].na; ia++) { - //std::cout << "\n it = " << it << " ia = " << ia << std::endl; - std::complex *sk = this->psf->get_sk(ik, it, ia, wfc_basis); + // std::cout << "\n it = " << it << " ia = " << ia << std::endl; + std::complex* sk = this->psf->get_sk(ik, it, ia, wfc_basis); //------------------------------------------------------- // Calculate G space 3D wave functions //------------------------------------------------------- - for (int iw = 0;iw < GlobalC::ucell.atoms[it].ncpp.nchi;iw++) + for (int iw = 0; iw < GlobalC::ucell.atoms[it].ncpp.nchi; iw++) { if (GlobalC::ucell.atoms[it].ncpp.oc[iw] >= 0.0) { @@ -301,267 +314,307 @@ void WF_atomic::atomic_wfc(const int& ik, // WF_atomictions for k=0 that are real in real space //----------------------------------------------------- - //--------------------------------------------------------- - // flq = radial fourier transform of atomic orbitals chi - //--------------------------------------------------------- - for (int ig=0; ig1e-8||fabs(fact[1])>1e-8) - { - for(int is=0;is<2;is++) - { - if(fabs(fact[is])>1e-8) - { - const int ind = lmaxkb + soc.sph_ind(l,j,m,is); - ModuleBase::GlobalFunc::ZEROS(aux, np); - for(int n1=0;n1<2*l+1;n1++){ - const int lm = l*l +n1; - if(std::abs(soc.rotylm(n1,ind))>1e-8) { - for(int ig=0; ignpwx*is ) = lphase * fact[is] * sk[ig] * aux[ig] * flq[ig]; -} - } - else { - for(int ig=0; ignpwx*is) = std::complex(0.0 , 0.0); -} -} - }//is - index++; - }//if - }//m - }//if + fact[0] = soc.spinor(l, j, m, 0); + fact[1] = soc.spinor(l, j, m, 1); + if (fabs(fact[0]) > 1e-8 || fabs(fact[1]) > 1e-8) + { + for (int is = 0; is < 2; is++) + { + if (fabs(fact[is]) > 1e-8) + { + const int ind = lmaxkb + soc.sph_ind(l, j, m, is); + ModuleBase::GlobalFunc::ZEROS(aux, np); + for (int n1 = 0; n1 < 2 * l + 1; n1++) + { + const int lm = l * l + n1; + if (std::abs(soc.rotylm(n1, ind)) > 1e-8) + { + for (int ig = 0; ig < np; ig++) + { + aux[ig] += soc.rotylm(n1, ind) * ylm(lm, ig); + } + } + } + for (int ig = 0; ig < np; ig++) + { + wfcatom(index, ig + this->npwx * is) + = lphase * fact[is] * sk[ig] * aux[ig] * flq[ig]; + } + } + else + { + for (int ig = 0; ig < np; ig++) + { + wfcatom(index, ig + this->npwx * is) + = std::complex(0.0, 0.0); + } + } + } // is + index++; + } // if + } // m + } // if else - {//atomic_wfc_so_mag + { // atomic_wfc_so_mag double alpha, gamma; - std::complex fup,fdown; + std::complex fup, fdown; int nc; - //This routine creates two functions only in the case j=l+1/2 or exit in the other case - if(fabs(j-l+0.5)<1e-4) { continue; -} + // This routine creates two functions only in the case j=l+1/2 or exit in the other case + if (fabs(j - l + 0.5) < 1e-4) + { + continue; + } delete[] chiaux; - chiaux = new double [np]; - //Find the functions j= l- 1/2 - if(l==0) { - for(int ig=0;igGlobalC::ucell.natomwfc) { ModuleBase::WARNING_QUIT("GlobalC::wf.atomic_wfc()","error: too many wfcs"); -} - for(int ig = 0;ig GlobalC::ucell.natomwfc) + { + ModuleBase::WARNING_QUIT("GlobalC::wf.atomic_wfc()", "error: too many wfcs"); + } + for (int ig = 0; ig < np; ig++) { - aux[ig] = sk[ig] * ylm(lm,ig) * chiaux[ig]; + aux[ig] = sk[ig] * ylm(lm, ig) * chiaux[ig]; } - //rotate wfc as needed - //first rotation with angle alpha around (OX) - for(int ig = 0;ignpwx) = (cos(0.5 * gamma) - ModuleBase::IMAG_UNIT * sin(0.5*gamma)) * fdown; - //second rotation with angle gamma around(OZ) - fup = cos(0.5 * (alpha + ModuleBase::PI))*aux[ig]; - fdown = ModuleBase::IMAG_UNIT * sin(0.5 * (alpha + ModuleBase::PI))*aux[ig]; - wfcatom(index+2*l+1,ig) = (cos(0.5*gamma) + ModuleBase::IMAG_UNIT*sin(0.5*gamma))*fup; - wfcatom(index+2*l+1,ig+ this->npwx) = (cos(0.5*gamma) - ModuleBase::IMAG_UNIT*sin(0.5*gamma))*fdown; + fdown = ModuleBase::IMAG_UNIT * sin(0.5 * alpha) * aux[ig]; + // build the orthogonal wfc + // first rotation with angle (alpha + ModuleBase::PI) around (OX) + wfcatom(index, ig) + = (cos(0.5 * gamma) + ModuleBase::IMAG_UNIT * sin(0.5 * gamma)) * fup; + wfcatom(index, ig + this->npwx) + = (cos(0.5 * gamma) - ModuleBase::IMAG_UNIT * sin(0.5 * gamma)) * fdown; + // second rotation with angle gamma around(OZ) + fup = cos(0.5 * (alpha + ModuleBase::PI)) * aux[ig]; + fdown = ModuleBase::IMAG_UNIT * sin(0.5 * (alpha + ModuleBase::PI)) * aux[ig]; + wfcatom(index + 2 * l + 1, ig) + = (cos(0.5 * gamma) + ModuleBase::IMAG_UNIT * sin(0.5 * gamma)) * fup; + wfcatom(index + 2 * l + 1, ig + this->npwx) + = (cos(0.5 * gamma) - ModuleBase::IMAG_UNIT * sin(0.5 * gamma)) * fdown; } index++; } - index += 2*l +1; + index += 2 * l + 1; } } else - {//atomic_wfc_nc + { // atomic_wfc_nc double alpha, gamman; std::complex fup, fdown; - //alpha = GlobalC::ucell.magnet.angle1_[it]; - //gamman = -GlobalC::ucell.magnet.angle2_[it] + 0.5*ModuleBase::PI; + // alpha = GlobalC::ucell.magnet.angle1_[it]; + // gamman = -GlobalC::ucell.magnet.angle2_[it] + 0.5*ModuleBase::PI; alpha = GlobalC::ucell.atoms[it].angle1[ia]; gamman = -1 * GlobalC::ucell.atoms[it].angle2[ia] + 0.5 * ModuleBase::PI; - for(int m = 0;m<2*l+1;m++) + for (int m = 0; m < 2 * l + 1; m++) { - const int lm = l*l +m; - if(index+2*l+1>GlobalC::ucell.natomwfc) { ModuleBase::WARNING_QUIT("GlobalC::wf.atomic_wfc()","error: too many wfcs"); -} - for(int ig = 0;ig GlobalC::ucell.natomwfc) + { + ModuleBase::WARNING_QUIT("GlobalC::wf.atomic_wfc()", "error: too many wfcs"); + } + for (int ig = 0; ig < np; ig++) { - aux[ig] = sk[ig] * ylm(lm,ig) * flq[ig]; + aux[ig] = sk[ig] * ylm(lm, ig) * flq[ig]; } - //rotate function - //first, rotation with angle alpha around(OX) - for(int ig = 0;ignpwx) = (cos(0.5 * gamman) - ModuleBase::IMAG_UNIT * sin(0.5*gamman)) * fdown; - //second rotation with angle gamma around(OZ) - fup = cos(0.5 * (alpha + ModuleBase::PI)) * aux[ig]; - fdown = ModuleBase::IMAG_UNIT * sin(0.5 * (alpha + ModuleBase::PI)) * aux[ig]; - wfcatom(index+2*l+1,ig) = (cos(0.5*gamman) + ModuleBase::IMAG_UNIT*sin(0.5*gamman))*fup; - wfcatom(index+2*l+1,ig+ this->npwx) = (cos(0.5*gamman) - ModuleBase::IMAG_UNIT*sin(0.5*gamman))*fdown; + fup = cos(0.5 * alpha) * aux[ig]; + fdown = ModuleBase::IMAG_UNIT * sin(0.5 * alpha) * aux[ig]; + // build the orthogonal wfc + // first rotation with angle(alpha+ModuleBase::PI) around(OX) + wfcatom(index, ig) + = (cos(0.5 * gamman) + ModuleBase::IMAG_UNIT * sin(0.5 * gamman)) * fup; + wfcatom(index, ig + this->npwx) + = (cos(0.5 * gamman) - ModuleBase::IMAG_UNIT * sin(0.5 * gamman)) * fdown; + // second rotation with angle gamma around(OZ) + fup = cos(0.5 * (alpha + ModuleBase::PI)) * aux[ig]; + fdown = ModuleBase::IMAG_UNIT * sin(0.5 * (alpha + ModuleBase::PI)) * aux[ig]; + wfcatom(index + 2 * l + 1, ig) + = (cos(0.5 * gamman) + ModuleBase::IMAG_UNIT * sin(0.5 * gamman)) * fup; + wfcatom(index + 2 * l + 1, ig + this->npwx) + = (cos(0.5 * gamman) - ModuleBase::IMAG_UNIT * sin(0.5 * gamman)) * fdown; } index++; } - index += 2*l+1; + index += 2 * l + 1; } - // index++; + // index++; } - else{//LSDA and nomagnet case + else + { // LSDA and nomagnet case - for (int m=0; m<2*l+1; m++) - { - const int lm = l * l + m; - for (int ig=0;igirindex[ir]; - const int ip = wfc_basis->fftixy2ip[ir]; + const int ip = wfc_basis->fftixy2ip[ir]; const int nz = wfc_basis->nz; - if(ip == 0 && GlobalV::RANK_IN_POOL ==0) - { - for(int iz=0; izirindex[ir]; - const int ip = wfc_basis->fftixy2ip[ir]; + const int ip = wfc_basis->fftixy2ip[ir]; const int nz = wfc_basis->nz; - if(ip == 0 && GlobalV::RANK_IN_POOL ==0) - { - for(int iz=0; iz* psi, srand(unsigned(PARAM.inp.pw_seed + GlobalC::Pkpoints.startk_pool[GlobalV::MY_POOL] + ik)); const int nxy = wfc_basis->fftnxy; const int nz = wfc_basis->nz; - const int nstnz = wfc_basis->nst*nz; + const int nstnz = wfc_basis->nst * nz; - FPTYPE *stickrr = new FPTYPE[nz]; - FPTYPE *stickarg = new FPTYPE[nz]; - FPTYPE *tmprr = new FPTYPE[nstnz]; - FPTYPE *tmparg = new FPTYPE[nstnz]; - for (int iw = iw_start ;iw < iw_end;iw++) + FPTYPE* stickrr = new FPTYPE[nz]; + FPTYPE* stickarg = new FPTYPE[nz]; + FPTYPE* tmprr = new FPTYPE[nstnz]; + FPTYPE* tmparg = new FPTYPE[nstnz]; + for (int iw = iw_start; iw < iw_end; iw++) { std::complex* ppsi = &(psi[iw * this->npwx * PARAM.globalv.npol]); int startig = 0; - for(int ipol = 0 ; ipol < PARAM.globalv.npol ; ++ipol) + for (int ipol = 0; ipol < PARAM.globalv.npol; ++ipol) { - - for(int ir=0; ir < nxy; ir++) - { - if(wfc_basis->fftixy2ip[ir] < 0) { continue; -} - if(GlobalV::RANK_IN_POOL==0) - { - for(int iz=0; izfftixy2ip[ir] < 0) + { + continue; + } + if (GlobalV::RANK_IN_POOL == 0) + { + for (int iz = 0; iz < nz; iz++) + { + stickrr[iz] = std::rand() / FPTYPE(RAND_MAX); + stickarg[iz] = std::rand() / FPTYPE(RAND_MAX); + } + } + stick_to_pool(stickrr, ir, tmprr, wfc_basis); stick_to_pool(stickarg, ir, tmparg, wfc_basis); - } + } - for (int ig = 0;ig < ng;ig++) + for (int ig = 0; ig < ng; ig++) { - const FPTYPE rr = tmprr[wfc_basis->getigl2isz(ik,ig)]; - const FPTYPE arg= ModuleBase::TWO_PI * tmparg[wfc_basis->getigl2isz(ik,ig)]; - const FPTYPE gk2 = wfc_basis->getgk2(ik,ig); - ppsi[ig+startig] = std::complex(rr * cos(arg), rr * sin(arg)) / FPTYPE(gk2 + 1.0); + const FPTYPE rr = tmprr[wfc_basis->getigl2isz(ik, ig)]; + const FPTYPE arg = ModuleBase::TWO_PI * tmparg[wfc_basis->getigl2isz(ik, ig)]; + const FPTYPE gk2 = wfc_basis->getgk2(ik, ig); + ppsi[ig + startig] = std::complex(rr * cos(arg), rr * sin(arg)) / FPTYPE(gk2 + 1.0); } - for(int ig=ng ;ig(0.0, 0.0); + ppsi[ig + startig] = std::complex(0.0, 0.0); } startig += npwx; } @@ -661,27 +716,29 @@ void WF_atomic::random_t(std::complex* psi, srand(unsigned(PARAM.inp.pw_seed + ik)); } #endif // __MPI - for (int iw = iw_start ;iw < iw_end;iw++) + for (int iw = iw_start; iw < iw_end; iw++) { std::complex* ppsi = &(psi[iw * this->npwx * PARAM.globalv.npol]); - for (int ig = 0;ig < ng;ig++) + for (int ig = 0; ig < ng; ig++) { - const FPTYPE rr = std::rand()/FPTYPE(RAND_MAX); //qianrui add RAND_MAX - const FPTYPE arg= ModuleBase::TWO_PI * std::rand()/FPTYPE(RAND_MAX); - const FPTYPE gk2 = wfc_basis->getgk2(ik,ig); + const FPTYPE rr = std::rand() / FPTYPE(RAND_MAX); // qianrui add RAND_MAX + const FPTYPE arg = ModuleBase::TWO_PI * std::rand() / FPTYPE(RAND_MAX); + const FPTYPE gk2 = wfc_basis->getgk2(ik, ig); ppsi[ig] = std::complex(rr * cos(arg), rr * sin(arg)) / FPTYPE(gk2 + 1.0); } - if(PARAM.globalv.npol==2) {for (int ig = this->npwx;ig < this->npwx + ng;ig++) + if (PARAM.globalv.npol == 2) { - const FPTYPE rr = std::rand()/FPTYPE(RAND_MAX); - const FPTYPE arg= ModuleBase::TWO_PI * std::rand()/FPTYPE(RAND_MAX); - const FPTYPE gk2 = wfc_basis->getgk2(ik,ig-this->npwx); - ppsi[ig] = std::complex(rr * cos(arg), rr * sin(arg)) / FPTYPE(gk2 + 1.0); + for (int ig = this->npwx; ig < this->npwx + ng; ig++) + { + const FPTYPE rr = std::rand() / FPTYPE(RAND_MAX); + const FPTYPE arg = ModuleBase::TWO_PI * std::rand() / FPTYPE(RAND_MAX); + const FPTYPE gk2 = wfc_basis->getgk2(ik, ig - this->npwx); + ppsi[ig] = std::complex(rr * cos(arg), rr * sin(arg)) / FPTYPE(gk2 + 1.0); + } } -} } #ifdef __MPI -// #if ((!defined __CUDA) && (!defined __ROCM)) + // #if ((!defined __CUDA) && (!defined __ROCM)) } // #endif // ((!defined __CUDA) && (!defined __ROCM)) #endif // __MPI @@ -702,42 +759,44 @@ void WF_atomic::atomicrandom(ModuleBase::ComplexMatrix& psi, srand(unsigned(PARAM.inp.pw_seed + GlobalC::Pkpoints.startk_pool[GlobalV::MY_POOL] + ik)); const int nxy = wfc_basis->fftnxy; const int nz = wfc_basis->nz; - const int nstnz = wfc_basis->nst*nz; + const int nstnz = wfc_basis->nst * nz; - double *stickrr = new double[nz]; - double *stickarg = new double[nz]; - double *tmprr = new double[nstnz]; - double *tmparg = new double[nstnz]; - for (int iw = iw_start ;iw < iw_end;iw++) + double* stickrr = new double[nz]; + double* stickarg = new double[nz]; + double* tmprr = new double[nstnz]; + double* tmparg = new double[nstnz]; + for (int iw = iw_start; iw < iw_end; iw++) { int startig = 0; - for(int ipol = 0 ; ipol < PARAM.globalv.npol ; ++ipol) + for (int ipol = 0; ipol < PARAM.globalv.npol; ++ipol) { - for(int ir=0; ir < nxy; ir++) - { - if(wfc_basis->fftixy2ip[ir] < 0) { continue; -} - if(GlobalV::RANK_IN_POOL==0) - { - for(int iz=0; izfftixy2ip[ir] < 0) + { + continue; + } + if (GlobalV::RANK_IN_POOL == 0) + { + for (int iz = 0; iz < nz; iz++) + { + stickrr[iz] = std::rand() / double(RAND_MAX); + stickarg[iz] = std::rand() / double(RAND_MAX); + } + } + stick_to_pool(stickrr, ir, tmprr, wfc_basis); stick_to_pool(stickarg, ir, tmparg, wfc_basis); - } + } - for (int ig = 0;ig < ng;ig++) + for (int ig = 0; ig < ng; ig++) { const double rr = tmprr[wfc_basis->ig2isz[ig]]; - const double arg= ModuleBase::TWO_PI * tmparg[wfc_basis->ig2isz[ig]]; - psi(iw,startig+ig) *= (1.0 + 0.05 * std::complex(rr * cos(arg), rr * sin(arg))); + const double arg = ModuleBase::TWO_PI * tmparg[wfc_basis->ig2isz[ig]]; + psi(iw, startig + ig) *= (1.0 + 0.05 * std::complex(rr * cos(arg), rr * sin(arg))); } - for(int ig=ng ;ig(0.0, 0.0); + psi(iw, startig + ig) = std::complex(0.0, 0.0); } startig += npwx; } @@ -752,24 +811,24 @@ void WF_atomic::atomicrandom(ModuleBase::ComplexMatrix& psi, #else if (PARAM.inp.pw_seed > 0) // qianrui add 2021-8-13 { - srand(unsigned(PARAM.inp.pw_seed + GlobalC::Pkpoints.startk_pool[GlobalV::MY_POOL] + ik)); - } + srand(unsigned(PARAM.inp.pw_seed + GlobalC::Pkpoints.startk_pool[GlobalV::MY_POOL] + ik)); + } #endif double rr, arg; - for (int iw = iw_start ;iw < iw_end;iw++) - { - int startig = 0; - for(int ip = 0 ; ip < PARAM.globalv.npol; ++ip) - { - for(int ig = 0 ; ig < npw ; ++ig) - { - rr = rand()/double(RAND_MAX); - arg = ModuleBase::TWO_PI * rand()/double(RAND_MAX); - psi(iw,startig+ig) *= (1.0 + 0.05 * std::complex(rr * cos(arg), rr * sin(arg))); - } - startig += npwx; - } - } + for (int iw = iw_start; iw < iw_end; iw++) + { + int startig = 0; + for (int ip = 0; ip < PARAM.globalv.npol; ++ip) + { + for (int ig = 0; ig < npw; ++ig) + { + rr = rand() / double(RAND_MAX); + arg = ModuleBase::TWO_PI * rand() / double(RAND_MAX); + psi(iw, startig + ig) *= (1.0 + 0.05 * std::complex(rr * cos(arg), rr * sin(arg))); + } + startig += npwx; + } + } #ifdef __MPI } #endif diff --git a/source/module_psi/wf_atomic.h b/source/module_psi/wf_atomic.h index de2bc27efb..2a8fb17af6 100644 --- a/source/module_psi/wf_atomic.h +++ b/source/module_psi/wf_atomic.h @@ -7,13 +7,11 @@ #include "module_base/realarray.h" #include "module_basis/module_pw/pw_basis_k.h" #include "module_hamilt_pw/hamilt_pwdft/structure_factor.h" - #include "module_psi/psi.h" class WF_atomic { - public: - + public: WF_atomic(); ~WF_atomic(); @@ -22,27 +20,27 @@ class WF_atomic // ModuleBase::IntArray igk; #ifdef __CUDA - double *d_g2kin; + double* d_g2kin; #endif - ModuleBase::realArray table_local;//mohan add 2009-09-10 + ModuleBase::realArray table_local; // mohan add 2009-09-10 - //temporary psi for new code + // temporary psi for new code psi::Psi>* psi = nullptr; - ModuleBase::ComplexMatrix *wanf2 = nullptr; // wannier functions in the PW basis + ModuleBase::ComplexMatrix* wanf2 = nullptr; // wannier functions in the PW basis /** * @brief init a table with the radial Fourier transform of the atomic WF_atomictions * @param sf_in [out] the structure factor * @param tab_at [out] atomic table */ - void init_at_1(Structure_Factor *sf_in, ModuleBase::realArray* tab_at); + void init_at_1(Structure_Factor* sf_in, ModuleBase::realArray* tab_at); - void print_PAOs()const; + void print_PAOs() const; - public: //template change to public, will be refactor later. added by zhengdy 20230302 - int *irindex = nullptr; + public: // template change to public, will be refactor later. added by zhengdy 20230302 + int* irindex = nullptr; void atomic_wfc(const int& ik, const int& np, @@ -69,6 +67,7 @@ class WF_atomic const int iw_end, const int ik, const ModulePW::PW_Basis_K* wfc_basis); + void random(std::complex* psi, const int iw_start, const int iw_end, @@ -86,8 +85,9 @@ class WF_atomic void stick_to_pool(double* stick, const int& ir, double* out, const ModulePW::PW_Basis_K* wfc_basis) const; void stick_to_pool(float* stick, const int& ir, float* out, const ModulePW::PW_Basis_K* wfc_basis) const; #endif + private: - Structure_Factor *psf; + Structure_Factor* psf; }; -#endif +#endif