TABLE OF CONTENTS
ABINIT/mkrho [ Functions ]
NAME
mkrho
FUNCTION
Depending on option argument value: --Compute charge density rho(r) and rho(G) in electrons/bohr**3 from input wavefunctions, band occupations, and k point wts. --Compute kinetic energy density tau(r) and tau(G) in bohr**-5 from input wavefunctions, band occupations, and k point wts. --Compute a given element of the kinetic energy density tensor tau_{alpha,beta}(r) and tau_{alpha,beta}(G) in bohr**-5 from input wavefunctions, band occupations, and k point wts.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, LSI, AR, MB) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
cg(2,mcg)=wf in G space dtset <type(dataset_type)>=all input variables for this dataset | istwfk(nkpt)=input option parameter that describes the storage of wfs | mband=maximum number of bands | mgfft=maximum size of 1D FFTs | mkmem=Number of k points treated by this node | mpw=maximum allowed value for npw | nband(nkpt*nsppol)=number of bands to be included in summation | at each k point for each spin channel | nfft=(effective) number of FFT grid points (for this processor) | ngfft(18)=contain all needed information about 3D FFT, | see ~abinit/doc/variables/vargs.htm#ngfft | nkpt=number of k points | nspden=number of spin-density components | nsppol=1 for unpolarized, 2 for spin-polarized | nsym=number of symmetry elements in group (at least 1 for identity) | symafm(nsym)=(anti)ferromagnetic part of symmetry operations | symrel(3,3,nsym)=symmetry matrices in real space (integers) | wtk(nkpt)=k point weights (they sum to 1.0) gprimd(3,3)=dimensional reciprocal space primitive translations irrzon(nfft**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data kg(3,mpw*mkmem)=reduced planewave coordinates mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol mpi_enreg=informations about MPI parallelization npwarr(nkpt)=number of planewaves and boundary planewaves at each k occ(mband*nkpt*nsppol)= occupation numbers for each band (usually 2.0) at each k point option if 0: compute rhor (electron density) if 1: compute taur (kinetic energy density) (i.e. Trace over the kinetic energy density tensor) if 2: compute taur_{alpha,beta} !!NOT YET IMPLEMENTED (a given element of the kinetic energy density tensor) paw_dmft <type(paw_dmft_type)>= paw+dmft related data phnons(2,nfft**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases rprimd(3,3)=dimensional real space primitive translations tim_mkrho=timing code of the calling routine(can be set to 0 if not attributed) ucvol=unit cell volume (Bohr**3) wvl_den <type(wvl_denspot_type)>=density informations for wavelets wvl_wfs <type(wvl_projector_type)>=wavefunctions informations for wavelets
OUTPUT
rhog(2,nfft)=total electron density in G space rhor(nfft,nspden)=electron density in r space (if spin polarized, array contains total density in first half and spin-up density in second half) (for non-collinear magnetism, first element: total density, 3 next ones: mx,my,mz in units of hbar/2)
PARENTS
afterscfloop,energy,gstate,respfn,scfcv,vtorho
CHILDREN
bandfft_kpt_set_ikpt,fftpac,fourwf,prep_fourwf,prtrhomxmn sphereboundary,symrhg,timab,wrtout,wvl_mkrho,xmpi_sum
SOURCE
79 #if defined HAVE_CONFIG_H 80 #include "config.h" 81 #endif 82 83 #include "abi_common.h" 84 85 subroutine mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,& 86 & rhog,rhor,rprimd,tim_mkrho,ucvol,wvl_den,wvl_wfs,& 87 & option) !optional 88 89 use defs_basis 90 use defs_abitypes 91 use defs_wvltypes 92 use m_profiling_abi 93 use m_xmpi 94 use m_errors 95 96 use m_hamiltonian, only : gs_hamiltonian_type 97 use m_bandfft_kpt, only : bandfft_kpt_set_ikpt 98 use m_paw_dmft, only : paw_dmft_type 99 100 !This section has been created automatically by the script Abilint (TD). 101 !Do not modify the following lines by hand. 102 #undef ABI_FUNC 103 #define ABI_FUNC 'mkrho' 104 use interfaces_14_hidewrite 105 use interfaces_18_timing 106 use interfaces_32_util 107 use interfaces_52_fft_mpi_noabirule 108 use interfaces_53_ffts 109 use interfaces_66_wfs 110 use interfaces_67_common, except_this_one => mkrho 111 !End of the abilint section 112 113 implicit none 114 115 !Arguments ------------------------------------ 116 !scalars 117 integer,intent(in) :: mcg,tim_mkrho 118 integer,intent(in),optional :: option 119 real(dp),intent(in) :: ucvol 120 type(MPI_type),intent(inout) :: mpi_enreg 121 type(dataset_type),intent(in) :: dtset 122 type(paw_dmft_type), intent(in) :: paw_dmft 123 type(wvl_wf_type),intent(inout) :: wvl_wfs 124 type(wvl_denspot_type), intent(inout) :: wvl_den 125 !no_abirules 126 !nfft**(1-1/nsym) is 1 if nsym==1, and nfft otherwise 127 integer, intent(in) :: irrzon(dtset%nfft**(1-1/dtset%nsym),2, & 128 & (dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)) 129 integer, intent(in) :: kg(3,dtset%mpw*dtset%mkmem),npwarr(dtset%nkpt) 130 real(dp), intent(in) :: gprimd(3,3) 131 real(dp), intent(in) :: cg(2,mcg) 132 real(dp), intent(in) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol) 133 !nfft**(1-1/nsym) is 1 if nsym==1, and nfft otherwise 134 real(dp), intent(in) :: phnons(2,(dtset%ngfft(1)*dtset%ngfft(2)*dtset%ngfft(3))**(1-1/dtset%nsym), & 135 & (dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)) 136 real(dp), intent(in) :: rprimd(3,3) 137 real(dp), intent(out) :: rhor(dtset%nfft,dtset%nspden),rhog(2,dtset%nfft) 138 139 !Local variables------------------------------- 140 !scalars 141 integer,save :: nskip=0 142 integer :: alpha,use_nondiag_occup_dmft,bdtot_index,beta,blocksize,iband,iband1,ibandc1,ib,iblock,icg,ierr 143 integer :: ifft,ikg,ikpt,ioption,ipw,ipwsp,ishf,ispden,ispinor,ispinor_index 144 integer :: isppol,istwf_k,jspinor_index 145 integer :: me,my_nspinor,n1,n2,n3,n4,n5,n6,nalpha,nband_k,nbandc1,nbdblock,nbeta 146 integer :: ndat,nfftot,npw_k,spaceComm,tim_fourwf 147 real(dp) :: kpt_cart,kg_k_cart,gp2pi1,gp2pi2,gp2pi3,cwftmp 148 real(dp) :: weight,weight_i 149 character(len=500) :: message 150 !arrays 151 integer,allocatable :: gbound(:,:),kg_k(:,:) 152 logical :: locc_test,nspinor1TreatedByThisProc,nspinor2TreatedByThisProc 153 real(dp) :: dummy(2,1),tsec(2) 154 real(dp),allocatable :: cwavef(:,:,:),cwavefb(:,:,:),cwavef_x(:,:) 155 real(dp),allocatable :: cwavef_y(:,:),cwavefb_x(:,:),cwavefb_y(:,:),kg_k_cart_block(:) 156 real(dp),allocatable :: occ_k(:),rhoaug(:,:,:),rhoaug_down(:,:,:) 157 real(dp),allocatable :: rhoaug_mx(:,:,:),rhoaug_my(:,:,:),rhoaug_up(:,:,:) 158 real(dp),allocatable :: taur_alphabeta(:,:,:,:),wfraug(:,:,:,:) 159 160 ! ************************************************************************* 161 162 DBG_ENTER("COLL") 163 164 call timab(790+tim_mkrho,1,tsec) 165 call timab(799,1,tsec) 166 167 if(.not.(present(option))) then 168 ioption=0 169 else 170 ioption=option 171 end if 172 173 if(ioption/=0.and.paw_dmft%use_sc_dmft==1) then 174 message = ' option argument value of this routines should be 0 if usedmft=1. ' 175 MSG_ERROR(message) 176 end if 177 if(paw_dmft%use_sc_dmft/=0) then 178 nbandc1=(paw_dmft%mbandc-1)*paw_dmft%use_sc_dmft+1 179 else 180 nbandc1=1 181 end if 182 use_nondiag_occup_dmft=0 183 184 !if(dtset%nspinor==2.and.paw_dmft%use_sc_dmft==1) then 185 !write(message, '(a,a,a,a)' )ch10,& 186 !& ' mkrho : ERROR -',ch10,& 187 !& ' nspinor argument value of this routines should be 1 if usedmft=1. ' 188 !call wrtout(std_out,message,'COLL') 189 !end if 190 191 my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor) 192 if (mpi_enreg%paral_spinor==0) then 193 ispinor_index=1;jspinor_index=1 194 nspinor1TreatedByThisProc=.true. 195 nspinor2TreatedByThisProc=(dtset%nspinor==2) 196 else 197 ispinor_index=mpi_enreg%me_spinor+1;jspinor_index=3-ispinor_index 198 nspinor1TreatedByThisProc=(mpi_enreg%me_spinor==0) 199 nspinor2TreatedByThisProc=(mpi_enreg%me_spinor==1) 200 end if 201 202 !Set local variable which depend on option argument 203 204 !nalpha*nbeta is the number of element of the kinetic energy density tensor 205 !to be computed in the irreducible Brillouin Zone (BZ) to get the result in the full BZ. 206 !In case of electron density calculation, nalpha=nbeta=1 207 select case (ioption) 208 case (0) 209 nalpha = 1 210 nbeta = 1 211 case (1) 212 nalpha = 3 213 nbeta = 1 214 ABI_ALLOCATE(taur_alphabeta,(dtset%nfft,dtset%nspden,3,1)) 215 case (2) 216 nalpha = 3 217 nbeta = 3 218 ABI_ALLOCATE(taur_alphabeta,(dtset%nfft,dtset%nspden,3,3)) 219 case default 220 MSG_BUG(' ioption argument value should be 0,1 or 2.') 221 end select 222 223 !Init me 224 me=mpi_enreg%me_kpt 225 !zero the charge density array in real space 226 !$OMP PARALLEL DO COLLAPSE(2) 227 do ispden=1,dtset%nspden 228 do ifft=1,dtset%nfft 229 rhor(ifft,ispden)=zero 230 end do 231 end do 232 233 !WVL - Branching with a separate mkrho procedure 234 !in wavelet. 235 if (dtset%usewvl == 1) then 236 select case(ioption) 237 case (0) 238 call wvl_mkrho(dtset, irrzon, mpi_enreg, phnons, rhor, wvl_wfs, wvl_den) 239 return 240 case (1) 241 ! call wvl_mkrho(dtset, mpi_enreg, occ, rhor, wvl_wfs, wvl_den) 242 message = ' Sorry, kinetic energy density (taur) is not yet implemented in wavelet formalism.' 243 MSG_ERROR(message) 244 case (2) 245 ! call wvl_mkrho(dtset, mpi_enreg, occ, rhor, wvl_wfs, wvl_den) 246 message = ' Sorry, kinetic energy density tensor (taur_(alpha,beta)) is not yet implemented in wavelet formalism.' 247 MSG_BUG(message) 248 end select 249 end if 250 !WVL - Following is done in plane waves. 251 252 !start loop over alpha and beta 253 254 do alpha=1,nalpha 255 do beta=1,nbeta 256 257 ! start loop over spin and k points 258 bdtot_index=0 259 icg=0 260 261 ! n4,n5,n6 are FFT dimensions, modified to avoir cache trashing 262 n1 = dtset%ngfft(1) ; n2 = dtset%ngfft(2) ; n3 = dtset%ngfft(3) 263 n4 = dtset%ngfft(4) ; n5 = dtset%ngfft(5) ; n6 = dtset%ngfft(6) 264 ndat = 1 ; if (mpi_enreg%paral_kgb==1) ndat = mpi_enreg%bandpp 265 ABI_ALLOCATE(cwavef,(2,dtset%mpw,my_nspinor)) 266 ABI_ALLOCATE(rhoaug,(n4,n5,n6)) 267 ABI_ALLOCATE(wfraug,(2,n4,n5,n6*ndat)) 268 ABI_ALLOCATE(cwavefb,(2,dtset%mpw*paw_dmft%use_sc_dmft,my_nspinor)) 269 if(dtset%nspden==4) then 270 ABI_ALLOCATE(rhoaug_up,(n4,n5,n6)) 271 ABI_ALLOCATE(rhoaug_down,(n4,n5,n6)) 272 ABI_ALLOCATE(rhoaug_mx,(n4,n5,n6)) 273 ABI_ALLOCATE(rhoaug_my,(n4,n5,n6)) 274 rhoaug_up(:,:,:)=zero 275 rhoaug_down(:,:,:)=zero 276 rhoaug_mx(:,:,:)=zero 277 rhoaug_my(:,:,:)=zero 278 end if 279 280 do isppol=1,dtset%nsppol 281 ikg=0 282 283 rhoaug(:,:,:)=zero 284 do ikpt=1,dtset%nkpt 285 286 nband_k = dtset%nband(ikpt+(isppol-1)*dtset%nkpt) 287 npw_k=npwarr(ikpt) 288 istwf_k = dtset%istwfk(ikpt) 289 290 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then 291 bdtot_index=bdtot_index+nband_k 292 cycle 293 end if 294 295 ABI_ALLOCATE(gbound,(2*dtset%mgfft+8,2)) 296 ABI_ALLOCATE(kg_k,(3,npw_k)) 297 298 kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 299 call sphereboundary(gbound,istwf_k,kg_k,dtset%mgfft,npw_k) 300 301 ! Loop over bands to fft and square for rho(r) 302 ! Shoulb be changed to treat bands by batch always 303 304 if(mpi_enreg%paral_kgb /= 1) then ! Not yet parallelized on spinors 305 do iband=1,nband_k 306 ! if(paw_dmft%use_sc_dmft==1) then 307 ! write(std_out,*) 'iband ',iband,occ(iband+bdtot_index),paw_dmft%occnd(iband,iband,ikpt,isppol) 308 ! else 309 ! write(std_out,*) 'iband ',iband,occ(iband+bdtot_index) 310 ! endif 311 if(mpi_enreg%paralbd==1)then 312 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband,isppol,me)) cycle 313 end if 314 do ibandc1=1,nbandc1 ! in case of DMFT 315 ! Check if DMFT and only treat occupied states (check on occ.) 316 if(paw_dmft%use_sc_dmft == 1) then 317 iband1 = paw_dmft%include_bands(ibandc1) 318 if(paw_dmft%band_in(iband)) then 319 if(.not. paw_dmft%band_in(iband1)) stop 320 use_nondiag_occup_dmft = 1 321 locc_test = abs(paw_dmft%occnd(1,iband,iband1,ikpt,isppol)) +& 322 & abs(paw_dmft%occnd(2,iband,iband1,ikpt,isppol))>tol8 323 ! write(std_out,*) "mkrho,ikpt,iband,use_occnd",ikpt,iband 324 else 325 use_nondiag_occup_dmft = 0 326 locc_test = abs(occ(iband+bdtot_index))>tol8 327 if(ibandc1 /=1 .and. .not. paw_dmft%band_in(iband)) cycle 328 end if 329 else 330 use_nondiag_occup_dmft = 0 331 locc_test = abs(occ(iband+bdtot_index))>tol8 332 end if 333 334 if (locc_test) then 335 ! Obtain Fourier transform in fft box and accumulate the density or the kinetic energy density 336 ! Not yet parallise on nspinor if paral_kgb non equal to 1 337 ipwsp=(iband-1)*npw_k*my_nspinor +icg 338 cwavef(:,1:npw_k,1)=cg(:,1+ipwsp:ipwsp+npw_k) 339 if (my_nspinor==2) cwavef(:,1:npw_k,2)=cg(:,ipwsp+npw_k+1:ipwsp+2*npw_k) 340 341 if(ioption==1)then 342 ! Multiplication by 2pi i (k+G)_alpha 343 gp2pi1=gprimd(alpha,1)*two_pi ; gp2pi2=gprimd(alpha,2)*two_pi ; gp2pi3=gprimd(alpha,3)*two_pi 344 kpt_cart=gp2pi1*dtset%kptns(1,ikpt)+gp2pi2*dtset%kptns(2,ikpt)+gp2pi3*dtset%kptns(3,ikpt) 345 do ispinor=1,my_nspinor 346 do ipw=1,npw_k 347 kg_k_cart=gp2pi1*kg_k(1,ipw)+gp2pi2*kg_k(2,ipw)+gp2pi3*kg_k(3,ipw)+kpt_cart 348 cwftmp=-cwavef(2,ipw,ispinor)*kg_k_cart 349 cwavef(2,ipw,ispinor)=cwavef(1,ipw,ispinor)*kg_k_cart 350 cwavef(1,ipw,ispinor)=cwftmp 351 end do 352 end do 353 else if(ioption==2)then 354 message = ' Sorry, kinetic energy density tensor (taur_(alpha,beta)) is not yet implemented.' 355 MSG_ERROR(message) 356 end if 357 358 ! Non diag occupation in DMFT. 359 if(use_nondiag_occup_dmft==1) then 360 ipwsp=(iband1-1)*npw_k*my_nspinor +icg 361 cwavefb(:,1:npw_k,1)=cg(:,1+ipwsp:ipwsp+npw_k) 362 if (my_nspinor==2) cwavefb(:,1:npw_k,2)=cg(:,ipwsp+npw_k+1:ipwsp+2*npw_k) 363 weight =paw_dmft%occnd(1,iband,iband1,ikpt,isppol)*dtset%wtk(ikpt)/ucvol 364 if(dtset%nspinor==1) weight_i=zero 365 if(dtset%nspinor==2) weight_i=paw_dmft%occnd(2,iband,iband1,ikpt,isppol)*dtset%wtk(ikpt)/ucvol 366 367 else 368 weight=occ(iband+bdtot_index)*dtset%wtk(ikpt)/ucvol 369 weight_i=weight 370 end if 371 372 if(mpi_enreg%paralbd==0) tim_fourwf=3 373 if(mpi_enreg%paralbd==1) tim_fourwf=6 374 375 ! The same section of code is also found in vtowfk.F90 : should be rationalized ! 376 377 call fourwf(1,rhoaug,cwavef(:,:,1),dummy,wfraug,gbound,gbound,& 378 & istwf_k,kg_k,kg_k,dtset%mgfft,mpi_enreg,1,dtset%ngfft,& 379 & npw_k,1,n4,n5,n6,1,dtset%paral_kgb,tim_fourwf,weight,weight_i,& 380 & use_ndo=use_nondiag_occup_dmft,fofginb=cwavefb(:,:,1),& 381 & use_gpu_cuda=dtset%use_gpu_cuda) 382 383 384 if(dtset%nspinor==2)then 385 ! DEBUG GZ !To obtain a x-directed magnetization(test) 386 ! cwavef1(1,1:npw_k)=-cwavef(2,1:npw_k) 387 ! cwavef1(2,1:npw_k)= cwavef(1,1:npw_k) 388 ! ENDDEBUG 389 390 if(dtset%nspden==1) then 391 ! We need only the total density : accumulation continues on top of rhoaug 392 393 call fourwf(1,rhoaug,cwavef(:,:,2),dummy,wfraug,gbound,gbound,& 394 & istwf_k,kg_k,kg_k,dtset%mgfft,mpi_enreg,1,dtset%ngfft,& 395 & npw_k,1,n4,n5,n6,1,dtset%paral_kgb,tim_fourwf,weight,weight_i,& 396 & use_ndo=use_nondiag_occup_dmft,fofginb=cwavefb(:,:,2),use_gpu_cuda=dtset%use_gpu_cuda) 397 398 399 else if(dtset%nspden==4) then 400 ! Build the four components of rho. We use only norm quantities and, so fourwf. 401 !$\sum_{n} f_n \Psi^{* \alpha}_n \Psi^{\alpha}_n =\rho^{\alpha \alpha}$ 402 !$\sum_{n} f_n (\Psi^{1}+\Psi^{2})^*_n (\Psi^{1}+\Psi^{2})_n=rho+m_x$ 403 !$\sum_{n} f_n (\Psi^{1}-i \Psi^{2})^*_n (\Psi^{1}-i \Psi^{2})_n=rho+m_y$ 404 ABI_ALLOCATE(cwavef_x,(2,npw_k)) 405 ABI_ALLOCATE(cwavef_y,(2,npw_k)) 406 ABI_ALLOCATE(cwavefb_x,(2,npw_k*paw_dmft%use_sc_dmft)) 407 ABI_ALLOCATE(cwavefb_y,(2,npw_k*paw_dmft%use_sc_dmft)) 408 !$(\Psi^{1}+\Psi^{2})$ 409 cwavef_x(:,:)=cwavef(:,1:npw_k,1)+cwavef(:,1:npw_k,2) 410 !$(\Psi^{1}-i \Psi^{2})$ 411 cwavef_y(1,:)=cwavef(1,1:npw_k,1)+cwavef(2,1:npw_k,2) 412 cwavef_y(2,:)=cwavef(2,1:npw_k,1)-cwavef(1,1:npw_k,2) 413 if(use_nondiag_occup_dmft==1) then 414 cwavefb_x(:,:)=cwavefb(:,1:npw_k,1)+cwavefb(:,1:npw_k,2) 415 cwavefb_y(1,:)=cwavefb(1,1:npw_k,1)+cwavefb(2,1:npw_k,2) 416 cwavefb_y(2,:)=cwavefb(2,1:npw_k,1)-cwavefb(1,1:npw_k,2) 417 end if 418 rhoaug_up(:,:,:)=rhoaug(:,:,:) !Already computed 419 420 call fourwf(1,rhoaug_down,cwavef(:,:,2),dummy,wfraug,gbound,gbound,& 421 & istwf_k,kg_k,kg_k,dtset%mgfft,mpi_enreg,1,dtset%ngfft,& 422 & npw_k,1,n4,n5,n6,1,dtset%paral_kgb,tim_fourwf,weight,weight_i,& 423 & use_ndo=use_nondiag_occup_dmft,fofginb=cwavefb(:,:,2),use_gpu_cuda=dtset%use_gpu_cuda) 424 425 call fourwf(1,rhoaug_mx,cwavef_x,dummy,wfraug,gbound,gbound,& 426 & istwf_k,kg_k,kg_k,dtset%mgfft,mpi_enreg,1,dtset%ngfft,& 427 & npw_k,1,n4,n5,n6,1,dtset%paral_kgb,tim_fourwf,weight,weight_i,& 428 & use_ndo=use_nondiag_occup_dmft,fofginb=cwavefb_x,use_gpu_cuda=dtset%use_gpu_cuda) 429 430 call fourwf(1,rhoaug_my,cwavef_y,dummy,wfraug,gbound,gbound,& 431 & istwf_k,kg_k,kg_k,dtset%mgfft,mpi_enreg,1,dtset%ngfft,& 432 & npw_k,1,n4,n5,n6,1,dtset%paral_kgb,tim_fourwf,weight,weight_i,& 433 & use_ndo=use_nondiag_occup_dmft,fofginb=cwavefb_y,use_gpu_cuda=dtset%use_gpu_cuda) 434 435 ABI_DEALLOCATE(cwavef_x) 436 ABI_DEALLOCATE(cwavef_y) 437 ABI_DEALLOCATE(cwavefb_x) 438 ABI_DEALLOCATE(cwavefb_y) 439 440 end if ! dtset%nspden/=4 441 442 443 end if 444 445 else 446 ! Accumulate the number of one-way 3D ffts skipped 447 nskip=nskip+1 448 end if ! abs(occ(iband+bdtot_index))>tol8 449 ! End loop on iband 450 end do ! iband1=1,(nband_k-1)*paw_dmft%use_sc_dmft+1 451 end do ! iband=1,nband_k 452 453 else !paral_kgb==1 454 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) cycle 455 456 call bandfft_kpt_set_ikpt(ikpt,mpi_enreg) 457 458 nbdblock=nband_k/(mpi_enreg%nproc_band * mpi_enreg%bandpp) 459 blocksize=nband_k/nbdblock 460 if(allocated(cwavef)) then 461 ABI_DEALLOCATE(cwavef) 462 end if 463 ABI_ALLOCATE(cwavef,(2,npw_k*blocksize,dtset%nspinor)) 464 if(ioption==1) then 465 ABI_ALLOCATE(kg_k_cart_block,(npw_k)) 466 end if 467 ABI_ALLOCATE(occ_k,(nband_k)) 468 occ_k(:)=occ(bdtot_index+1:bdtot_index+nband_k) 469 470 do iblock=1,nbdblock 471 if (dtset%nspinor==1) then 472 cwavef(:,1:npw_k*blocksize,1)=cg(:,1+(iblock-1)*npw_k*blocksize+icg:iblock*npw_k*blocksize+icg) 473 else 474 if (mpi_enreg%paral_spinor==0) then 475 ishf=(iblock-1)*npw_k*my_nspinor*blocksize+icg 476 do ib=1,blocksize 477 cwavef(:,(ib-1)*npw_k+1:ib*npw_k,1)=cg(:,1+(2*ib-2)*npw_k+ishf:(2*ib-1)*npw_k+ishf) 478 cwavef(:,(ib-1)*npw_k+1:ib*npw_k,2)=cg(:,1+(2*ib-1)*npw_k+ishf:ib*2*npw_k+ishf) 479 end do 480 else 481 ishf=(iblock-1)*npw_k*my_nspinor*blocksize+icg 482 do ib=1,blocksize 483 cwavef(:,(ib-1)*npw_k+1:ib*npw_k,ispinor_index)=& 484 & cg(:,1+(ib-1)*npw_k+ishf:ib*npw_k+ishf) 485 cwavef(:,(ib-1)*npw_k+1:ib*npw_k,jspinor_index)=zero 486 end do 487 call xmpi_sum(cwavef,mpi_enreg%comm_spinor,ierr) 488 end if 489 end if 490 if(ioption==1)then 491 ! Multiplication by 2pi i (k+G)_alpha 492 gp2pi1=gprimd(alpha,1)*two_pi ; gp2pi2=gprimd(alpha,2)*two_pi ; gp2pi3=gprimd(alpha,3)*two_pi 493 kpt_cart=gp2pi1*dtset%kptns(1,ikpt)+gp2pi2*dtset%kptns(2,ikpt)+gp2pi3*dtset%kptns(3,ikpt) 494 kg_k_cart_block(1:npw_k)=gp2pi1*kg_k(1,1:npw_k)+gp2pi2*kg_k(2,1:npw_k)+gp2pi3*kg_k(3,1:npw_k)+kpt_cart 495 do ib=1,blocksize 496 do ipw=1,npw_k 497 cwftmp=-cwavef(2,ipw+(ib-1)*npw_k,1)*kg_k_cart_block(ipw) 498 cwavef(2,ipw,1)=cwavef(1,ipw+(ib-1)*npw_k,1)*kg_k_cart_block(ipw) 499 cwavef(1,ipw,1)=cwftmp 500 if (my_nspinor==2) then 501 cwftmp=-cwavef(2,ipw+(ib-1)*npw_k,2)*kg_k_cart_block(ipw) 502 cwavef(2,ipw,2)=cwavef(1,ipw+(ib-1)*npw_k,2)*kg_k_cart_block(ipw) 503 cwavef(1,ipw,2)=cwftmp 504 end if 505 end do 506 end do 507 else if(ioption==2)then 508 message = ' Sorry, kinetic energy density tensor (taur_(alpha,beta)) is not yet implemented.' 509 MSG_ERROR(message) 510 end if 511 512 call timab(538,1,tsec) 513 if (nspinor1TreatedByThisProc) then 514 call prep_fourwf(rhoaug,blocksize,cwavef(:,:,1),wfraug,iblock,istwf_k,dtset%mgfft,mpi_enreg,& 515 & nband_k,ndat,dtset%ngfft,npw_k,n4,n5,n6,occ_k,1,ucvol,& 516 & dtset%wtk(ikpt),use_gpu_cuda=dtset%use_gpu_cuda) 517 end if 518 call timab(538,2,tsec) 519 if(dtset%nspinor==2)then 520 if (dtset%nspden==1) then 521 if (nspinor2TreatedByThisProc) then 522 call prep_fourwf(rhoaug,blocksize,cwavef(:,:,2),wfraug,& 523 & iblock,istwf_k,dtset%mgfft,mpi_enreg,& 524 & nband_k,ndat,dtset%ngfft,npw_k,n4,n5,n6,occ_k,1,ucvol,& 525 & dtset%wtk(ikpt),use_gpu_cuda=dtset%use_gpu_cuda) 526 end if 527 else if(dtset%nspden==4 ) then 528 ABI_ALLOCATE(cwavef_x,(2,npw_k*blocksize)) 529 ABI_ALLOCATE(cwavef_y,(2,npw_k*blocksize)) 530 cwavef_x(:,:)=cwavef(:,:,1)+cwavef(:,:,2) 531 cwavef_y(1,:)=cwavef(1,:,1)+cwavef(2,:,2) 532 cwavef_y(2,:)=cwavef(2,:,1)-cwavef(1,:,2) 533 call timab(538,1,tsec) 534 if (nspinor1TreatedByThisProc) then 535 call prep_fourwf(rhoaug_down,blocksize,cwavef(:,:,2),wfraug,& 536 & iblock,istwf_k,dtset%mgfft,mpi_enreg,& 537 & nband_k,ndat,dtset%ngfft,npw_k,n4,n5,n6,occ_k,1,ucvol,& 538 & dtset%wtk(ikpt),use_gpu_cuda=dtset%use_gpu_cuda) 539 end if 540 if (nspinor2TreatedByThisProc) then 541 call prep_fourwf(rhoaug_mx,blocksize,cwavef_x,wfraug,& 542 & iblock,istwf_k,dtset%mgfft,mpi_enreg,& 543 & nband_k,ndat,dtset%ngfft,npw_k,n4,n5,n6,occ_k,1,ucvol,& 544 & dtset%wtk(ikpt),use_gpu_cuda=dtset%use_gpu_cuda) 545 call prep_fourwf(rhoaug_my,blocksize,cwavef_y,wfraug,& 546 & iblock,istwf_k,dtset%mgfft,mpi_enreg,& 547 & nband_k,ndat,dtset%ngfft,npw_k,n4,n5,n6,occ_k,1,ucvol,& 548 & dtset%wtk(ikpt),use_gpu_cuda=dtset%use_gpu_cuda) 549 end if 550 call timab(538,2,tsec) 551 ABI_DEALLOCATE(cwavef_x) 552 ABI_DEALLOCATE(cwavef_y) 553 end if 554 end if 555 end do !iblock 556 if(ioption==1) then 557 ABI_DEALLOCATE(kg_k_cart_block) 558 end if 559 if (allocated(cwavef)) then 560 ABI_DEALLOCATE(cwavef) 561 end if 562 ABI_DEALLOCATE(occ_k) 563 end if 564 565 ABI_DEALLOCATE(gbound) 566 ABI_DEALLOCATE(kg_k) 567 568 bdtot_index=bdtot_index+nband_k 569 570 if (dtset%mkmem/=0) then 571 icg=icg+npw_k*my_nspinor*nband_k 572 ikg=ikg+npw_k 573 end if 574 575 ! End loop on ikpt: 576 end do 577 578 579 if(mpi_enreg%paral_kgb == 1) then 580 call bandfft_kpt_set_ikpt(-1,mpi_enreg) 581 if (dtset%nspden==4) then 582 ! Sum the contribution of the band and of the FFT 583 call xmpi_sum(rhoaug ,mpi_enreg%comm_bandspinorfft, ierr) 584 call xmpi_sum(rhoaug_down,mpi_enreg%comm_bandspinorfft, ierr) 585 call xmpi_sum(rhoaug_mx ,mpi_enreg%comm_bandspinorfft, ierr) 586 call xmpi_sum(rhoaug_my ,mpi_enreg%comm_bandspinorfft, ierr) 587 rhoaug_up(:,:,:) = rhoaug(:,:,:) 588 else 589 call xmpi_sum(rhoaug,mpi_enreg%comm_bandspinorfft,ierr) 590 end if 591 end if 592 593 ! Transfer density on augmented fft grid to normal fft grid in real space 594 ! Take also into account the spin, to place it correctly in rhor. 595 if(dtset%nspden==1 .or. dtset%nspden==2) then 596 call fftpac(isppol,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhor,rhoaug,1) 597 else if(dtset%nspden==4) then 598 ispden=1 599 call fftpac(ispden,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhor,rhoaug_up,1) 600 ispden=2 601 call fftpac(ispden,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhor,rhoaug_mx,1) 602 ispden=3 603 call fftpac(ispden,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhor,rhoaug_my,1) 604 ispden=4 605 call fftpac(ispden,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhor,rhoaug_down,1) 606 ABI_DEALLOCATE(rhoaug_up) 607 ABI_DEALLOCATE(rhoaug_down) 608 ABI_DEALLOCATE(rhoaug_mx) 609 ABI_DEALLOCATE(rhoaug_my) 610 end if 611 612 end do ! isppol=1,dtset%nsppol 613 614 if(allocated(cwavef)) then 615 ABI_DEALLOCATE(cwavef) 616 end if 617 ABI_DEALLOCATE(cwavefb) 618 ABI_DEALLOCATE(rhoaug) 619 ABI_DEALLOCATE(wfraug) 620 621 ! Recreate full rhor on all proc. 622 call timab(48,1,tsec) 623 call timab(71,1,tsec) 624 spaceComm=mpi_enreg%comm_cell 625 if (mpi_enreg%paral_hf==1)spaceComm=mpi_enreg%comm_kpt 626 if(mpi_enreg%paral_kgb==1)spaceComm=mpi_enreg%comm_kpt 627 call xmpi_sum(rhor,spaceComm,ierr) 628 call timab(71,2,tsec) 629 call timab(48,2,tsec) 630 631 call timab(799,2,tsec) 632 call timab(549,1,tsec) 633 634 if(ioption==1 .or. ioption==2) then 635 !$OMP PARALLEL DO COLLAPSE(2) 636 do ispden=1,dtset%nspden 637 do ifft=1,dtset%nfft 638 taur_alphabeta(ifft,ispden,alpha,beta) = rhor(ifft,ispden) 639 end do 640 end do 641 end if 642 643 end do ! beta=1,nbeta 644 end do ! alpha=1,nalpha 645 646 !Compute the trace over the kinetic energy density tensor. i.e. Sum of the 3 diagonal elements. 647 if(ioption==1)then 648 ! zero rhor array in real space 649 do ispden=1,dtset%nspden 650 !$OMP PARALLEL DO 651 do ifft=1,dtset%nfft 652 rhor(ifft,ispden)=zero 653 end do 654 end do 655 do alpha = 1, nalpha 656 !$OMP PARALLEL DO COLLAPSE(2) 657 do ispden=1,dtset%nspden 658 do ifft=1,dtset%nfft 659 rhor(ifft,ispden) = rhor(ifft,ispden) + taur_alphabeta(ifft,ispden,alpha,1) 660 end do 661 end do 662 end do 663 end if 664 665 nfftot=dtset%ngfft(1) * dtset%ngfft(2) * dtset%ngfft(3) 666 667 select case (ioption) 668 case(0,1) 669 call symrhg(1,gprimd,irrzon,mpi_enreg,dtset%nfft,nfftot,dtset%ngfft,dtset%nspden,dtset%nsppol,dtset%nsym,& 670 dtset%paral_kgb,phnons,rhog,rhor,rprimd,dtset%symafm,dtset%symrel) 671 if(ioption==1)then 672 !$OMP PARALLEL DO 673 do ifft=1,dtset%nfft 674 do ispden=1,dtset%nspden 675 rhor(ifft,ispden)=1.0d0/2.0d0*rhor(ifft,ispden) 676 end do 677 rhog(:,ifft)=1.0d0/2.0d0*rhog(:,ifft) 678 end do 679 end if 680 case(2) 681 message = ' Sorry, kinetic energy density tensor (taur_(alpha,beta)) is not yet implemented.' 682 MSG_BUG(message) 683 684 ! call symtaug(1,gprimd,irrzon,mpi_enreg,dtset%nfft,nfftot,dtset%ngfft,dtset%nspden,dtset%nsppol,dtset%nsym,& 685 ! dtset%paral_kgb,phnons,rhog,rhor,rprimd,dtset%symafm,dtset%symrel) 686 end select 687 688 call timab(549,2,tsec) 689 690 !We now have both rho(r) and rho(G), symmetrized, and if dtset%nsppol=2 691 !we also have the spin-up density, symmetrized, in rhor(:,2). 692 !In case of non collinear magnetism, we have rho,mx,my,mz. No symmetry is applied 693 694 call timab(799,1,tsec) 695 696 if(ioption==1 .or. ioption==2) then 697 ABI_DEALLOCATE(taur_alphabeta) 698 end if 699 700 !Find and print minimum and maximum total electron density 701 !(or total kinetic energy density, or total element of kinetic energy density tensor) and locations 702 call wrtout(std_out,'mkrho: echo density (plane-wave part only)','COLL') 703 call prtrhomxmn(std_out,mpi_enreg,dtset%nfft,dtset%ngfft,dtset%nspden,1,rhor,optrhor=ioption,ucvol=ucvol) 704 705 call timab(799,2,tsec) 706 call timab(790+tim_mkrho,2,tsec) 707 708 DBG_EXIT("COLL") 709 710 end subroutine mkrho