TABLE OF CONTENTS
ABINIT/dfpt_mkrho [ Functions ]
NAME
dfpt_mkrho
FUNCTION
Compute RF charge density rho1(r) and rho1(G) in electrons/bohr**3 from input RF and GS wavefunctions, band occupations, and k point weights.
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,mpw*nspinor*mband*mkmem*nsppol)=wf in G space cg1(2,mpw1*nspinor*mband*mk1mem*nsppol)=first-order wf in G space cplex=1 if rhor1 is real, 2 if rhor1 is complex gprimd(3,3)=dimensional reciprocal space primitive translations irrzon(nfft**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data istwfk_rbz(nkpt_rbz)=input option parameter that describes the storage of wfs kg(3,mpw*mkmem)=reduced planewave coordinates, GS data. kg1(3,mpw1*mkmem1)=reduced planewave coordinates, RF data. mband=maximum number of bands mgfft=maximum size of 1D FFTs mkmem=Number of k points treated by this node (GS data) mk1mem=Number of k points treated by this node (RF data) mpi_enreg=information about MPI parallelization mpw=maximum allowed value for npw (GS wfs) mpw1=maximum allowed value for npw1 (RF data) nband_rbz(nkpt_rbz*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_rbz=number of k points in the reduced Brillouin zone npwarr(nkpt_rbz)=number of planewaves and boundary planewaves at k points npwar1(nkpt_rbz)=number of planewaves and boundary planewaves at k+q points nspden=number of spin-density components nspinor=number of spinorial components of the wavefunctions nsppol=1 for unpolarized, 2 for spin-polarized nsym=number of symmetry elements in group (at least 1 for identity) occ_rbz(mband*nkpt_rbz*nsppol)=occupation numbers for each band (usually 2.0) at each k point of the reduced Brillouin zone phnons(2,nfft**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases rprimd(3,3)=dimensional real space primitive translations symafm(nsym)=(anti)ferromagnetic part of symmetry operations symrel(3,3,nsym)=symmetry matrices in real space (integers) ucvol=unit cell volume (Bohr**3) wtk_rbz(nkpt_rbz)=k point weights (they sum to 1.0).
OUTPUT
rhog1(2,nfft)=total electron density in G space rhor1(cplex*nfft,nspden)=electron density in r space (if spin polarized, array contains total density in first half and spin-up density in second half)
PARENTS
dfpt_looppert
CHILDREN
cg_zcopy,fftpac,fourwf,sphereboundary,symrhg,timab,wrtout,xmpi_sum
SOURCE
69 #if defined HAVE_CONFIG_H 70 #include "config.h" 71 #endif 72 73 #include "abi_common.h" 74 75 76 subroutine dfpt_mkrho(cg,cg1,cplex,gprimd,irrzon,istwfk_rbz,& 77 & kg,kg1,mband,mgfft,mkmem,mk1mem,mpi_enreg,mpw,mpw1,nband_rbz,& 78 & nfft,ngfft,nkpt_rbz,npwarr,npwar1,nspden,nspinor,nsppol,nsym,& 79 & occ_rbz,paral_kgb,phnons,rhog1,rhor1,rprimd,symafm,symrel,ucvol,wtk_rbz) 80 81 use defs_basis 82 use defs_abitypes 83 use m_profiling_abi 84 use m_errors 85 use m_cgtools 86 use m_xmpi 87 88 use m_io_tools, only : get_unit, iomode_from_fname 89 90 !This section has been created automatically by the script Abilint (TD). 91 !Do not modify the following lines by hand. 92 #undef ABI_FUNC 93 #define ABI_FUNC 'dfpt_mkrho' 94 use interfaces_14_hidewrite 95 use interfaces_18_timing 96 use interfaces_32_util 97 use interfaces_52_fft_mpi_noabirule 98 use interfaces_53_ffts 99 use interfaces_67_common 100 !End of the abilint section 101 102 implicit none 103 104 !Arguments ------------------------------------ 105 !scalars 106 integer,intent(in) :: cplex,mband,mgfft,mk1mem,mkmem,mpw,mpw1,nfft,nkpt_rbz 107 integer,intent(in) :: nspden,nspinor,nsppol,nsym,paral_kgb 108 real(dp),intent(in) :: ucvol 109 type(MPI_type),intent(in) :: mpi_enreg 110 !arrays 111 integer,intent(in) :: irrzon(nfft**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4)) 112 integer,intent(in) :: istwfk_rbz(nkpt_rbz),kg(3,mpw*mkmem),kg1(3,mpw1*mk1mem) 113 integer,intent(in) :: nband_rbz(nkpt_rbz*nsppol),ngfft(18),npwar1(nkpt_rbz) 114 integer,intent(in) :: npwarr(nkpt_rbz),symafm(nsym),symrel(3,3,nsym) 115 real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol) 116 real(dp),intent(in) :: cg1(2,mpw1*nspinor*mband*mk1mem*nsppol),gprimd(3,3) 117 real(dp),intent(in) :: occ_rbz(mband*nkpt_rbz*nsppol) 118 real(dp),intent(in) :: phnons(2,nfft**(1-1/nsym),(nspden/nsppol)-3*(nspden/4)) 119 real(dp),intent(in) :: rprimd(3,3),wtk_rbz(nkpt_rbz) 120 real(dp),intent(out) :: rhog1(2,nfft),rhor1(cplex*nfft,nspden) 121 122 !Local variables------------------------------- 123 !scalars 124 integer,parameter :: tim_fourwf7=7,tim_rwwf15=15 125 integer,save :: nskip=0 126 integer :: bdtot_index,i1,i2,i3,iband,icg,icg1,ierr,ifft,ikg,ptr 127 integer :: ikg1,ikpt,ispden,ispinor,isppol,istwf_k,ptr1,ptr2 128 integer :: me,n1,n2,n3,n4,n5,n6,nband_k,npw1_k 129 integer :: npw_k,spaceworld 130 real(dp) :: im0,im1,re0,re1,weight 131 real(dp) :: im0_up,im1_up,re0_up,re1_up,im0_down,im1_down,re0_down,re1_down 132 character(len=500) :: message 133 !arrays 134 integer,allocatable :: gbound(:,:),gbound1(:,:),kg1_k(:,:) 135 integer,allocatable :: kg_k(:,:) 136 real(dp) :: tsec(2) 137 real(dp),allocatable,target :: cwavef(:,:),cwavef1(:,:) 138 real(dp),allocatable :: dummy(:,:),rhoaug(:,:,:,:) 139 real(dp),allocatable :: rhoaug1(:,:,:,:),wfraug(:,:,:,:),wfraug1(:,:,:,:) 140 real(dp),allocatable :: wfraug1_up(:,:,:,:),wfraug1_down(:,:,:,:) 141 real(dp),allocatable :: wfraug_up(:,:,:,:),wfraug_down(:,:,:,:) 142 real(dp),allocatable :: cwave0_up(:,:),cwave0_down(:,:),cwave1_up(:,:),cwave1_down(:,:) 143 144 ! ************************************************************************* 145 146 !DBG_ENTER("COLL") 147 148 if(nspden==4)then 149 ! NOTE : see mkrho for the modifications needed for non-collinear treatment 150 write(message, '(a,a,a,a,a,a,a,a)' ) ch10,& 151 & ' dfpt_mkrho : WARNING -',ch10,& 152 & ' Linear-response calculations are under construction with nspden=4',ch10,& 153 & ' Action : modify value of nspden in input file unless you know what you are doing.' 154 ! call wrtout(ab_out,message,'COLL') 155 call wrtout(std_out,message,'COLL') 156 end if 157 158 !Init spaceworld 159 spaceworld=mpi_enreg%comm_cell 160 me=mpi_enreg%me_kpt 161 162 !zero the charge density array in real space 163 !$OMP PARALLEL DO 164 do ispden=1,nspden 165 do ifft=1,cplex*nfft 166 rhor1(ifft,ispden)=zero 167 end do 168 end do 169 170 !start loop over spin and k points 171 bdtot_index=0; icg=0; icg1=0 172 173 n1=ngfft(1); n2=ngfft(2); n3=ngfft(3) 174 n4=ngfft(4); n5=ngfft(5); n6=ngfft(6) !n4,n5,n6 are FFT dimensions, modified to avoid cache trashing 175 176 !Note that the dimensioning of cwavef and cwavef1 does not include nspinor 177 ABI_ALLOCATE(cwavef,(2,mpw)) 178 ABI_ALLOCATE(cwavef1,(2,mpw1)) 179 !Actually, rhoaug is not needed, except for strong dimensioning requirement 180 ABI_ALLOCATE(dummy,(2,1)) 181 ABI_ALLOCATE(rhoaug,(n4,n5,n6,nspinor**2)) 182 ABI_ALLOCATE(rhoaug1,(cplex*n4,n5,n6,nspinor**2)) 183 ABI_ALLOCATE(wfraug,(2,n4,n5,n6)) 184 ABI_ALLOCATE(wfraug1,(2,n4,n5,n6)) 185 186 ! EB FR Separate collinear and non-collinear magnetism 187 if (nspden /= 4) then ! EB FR nspden check 188 do isppol=1,nsppol 189 190 ikg=0; ikg1=0 191 192 rhoaug1(:,:,:,:)=zero 193 194 do ikpt=1,nkpt_rbz 195 196 nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz) 197 istwf_k=istwfk_rbz(ikpt) 198 npw_k=npwarr(ikpt) 199 npw1_k=npwar1(ikpt) 200 201 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then 202 bdtot_index=bdtot_index+nband_k 203 cycle 204 end if 205 206 ABI_ALLOCATE(gbound,(2*mgfft+8,2)) 207 ABI_ALLOCATE(kg_k,(3,npw_k)) 208 ABI_ALLOCATE(gbound1,(2*mgfft+8,2)) 209 ABI_ALLOCATE(kg1_k,(3,npw1_k)) 210 211 kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 212 call sphereboundary(gbound,istwf_k,kg_k,mgfft,npw_k) 213 214 kg1_k(:,1:npw1_k)=kg1(:,1+ikg1:npw1_k+ikg1) 215 call sphereboundary(gbound1,istwf_k,kg1_k,mgfft,npw1_k) 216 217 ! Loop over bands to fft and square for rho(r) 218 do iband=1,nband_k 219 if (mpi_enreg%proc_distrb(ikpt,iband,isppol)/=me) cycle 220 ! Only treat occupied states 221 if (abs(occ_rbz(iband+bdtot_index))>tol8) then 222 ! Treat separately the two spinor components 223 do ispinor=1,nspinor 224 ! Obtain Fourier transform in fft box and accumulate the density 225 ptr = 1 + (ispinor-1)*npw_k + (iband-1)*npw_k*nspinor + icg 226 call cg_zcopy(npw_k, cg(1,ptr), cwavef) 227 228 ! In these two calls, rhoaug, rhoaug1 and weight are dummy variables, and are not modified 229 call fourwf(1,rhoaug,cwavef,dummy,wfraug,gbound,gbound,& 230 & istwf_k,kg_k,kg_k,mgfft,mpi_enreg,1,ngfft,npw_k,1,n4,n5,n6,0,paral_kgb,tim_fourwf7,weight,weight) 231 232 ! TODO: here ispinor should be ispinorp to get full matrix and nspden 4 233 ptr = 1 + (ispinor-1)*npw1_k + (iband-1)*npw1_k*nspinor + icg1 234 call cg_zcopy(npw1_k, cg1(1,ptr), cwavef1) 235 236 call fourwf(cplex,rhoaug1,cwavef1,dummy,wfraug1,gbound1,gbound1,& 237 & istwf_k,kg1_k,kg1_k,mgfft,mpi_enreg,1,ngfft,npw1_k,1,n4,n5,n6,0,& 238 & paral_kgb,tim_fourwf7,weight,weight) 239 240 ! Compute the weight, note that the factor 2 is 241 ! not the spin factor (see Eq.44 of PRB55,10337 (1997)) 242 weight=two*occ_rbz(iband+bdtot_index)*wtk_rbz(ikpt)/ucvol 243 244 ! Accumulate density 245 if(cplex==2)then 246 !$OMP PARALLEL DO PRIVATE(im0,im1,re0,re1) 247 do i3=1,n3 248 do i2=1,n2 249 do i1=1,n1 250 re0=wfraug(1,i1,i2,i3) ; im0=wfraug(2,i1,i2,i3) 251 re1=wfraug1(1,i1,i2,i3); im1=wfraug1(2,i1,i2,i3) 252 rhoaug1(2*i1-1,i2,i3,1)=rhoaug1(2*i1-1,i2,i3,1)+weight*(re0*re1+im0*im1) 253 rhoaug1(2*i1 ,i2,i3,1)=rhoaug1(2*i1 ,i2,i3,1)+weight*(re0*im1-im0*re1) 254 end do 255 end do 256 end do 257 else 258 !$OMP PARALLEL DO 259 do i3=1,n3 260 do i2=1,n2 261 do i1=1,n1 262 rhoaug1(i1,i2,i3,1)=rhoaug1(i1,i2,i3,1)+& 263 & weight*( wfraug(1,i1,i2,i3)*wfraug1(1,i1,i2,i3) + wfraug(2,i1,i2,i3)*wfraug1(2,i1,i2,i3) ) 264 end do 265 end do 266 end do 267 end if ! cplex 268 end do ! ispinor 269 else !abs(occ_rbz(iband+bdtot_index))>tol8 270 nskip=nskip+1 ! if the state is not occupied. Accumulate the number of one-way 3D ffts skipped 271 end if ! abs(occ_rbz(iband+bdtot_index))>tol8 272 end do ! iband 273 274 ABI_DEALLOCATE(gbound) 275 ABI_DEALLOCATE(kg_k) 276 ABI_DEALLOCATE(gbound1) 277 ABI_DEALLOCATE(kg1_k) 278 279 bdtot_index=bdtot_index+nband_k 280 281 icg=icg+npw_k*nband_k 282 ikg=ikg+npw_k 283 icg1=icg1+npw1_k*nband_k 284 ikg1=ikg1+npw1_k 285 286 end do ! ikpt 287 288 if (xmpi_paral==0) then ! Write the number of one-way 3D ffts skipped until now 289 write(message,'(a,i8)')' mkrho3 : number of one-way 3D ffts skipped in mkrho3 until now =',nskip 290 call wrtout(std_out,message,'PERS') 291 end if 292 293 ! TODO: if n+magnetization basis is used for the density, need to rotate rhoaug1 to that now, before packing into rhor1 294 295 ! Transfer density on augmented fft grid to normal fft grid in real space 296 ! Take also into account the spin, to place it correctly in rhor1. 297 ! Note the use of cplex 298 call fftpac(isppol,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,ngfft,rhor1,rhoaug1,1) 299 300 end do ! loop over isppol spins 301 302 else ! nspden = 4 303 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 304 ! Part added for the non collinear magnetism 305 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 306 ! The same lines of code are in 72_response/accrho3.F90 307 ! TODO: merge these lines in a single routine??!! 308 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 309 310 ikg=0; ikg1=0 311 312 rhoaug1(:,:,:,:)=zero 313 314 do ikpt=1,nkpt_rbz 315 316 nband_k=nband_rbz(ikpt) 317 istwf_k=istwfk_rbz(ikpt) 318 npw_k=npwarr(ikpt) 319 npw1_k=npwar1(ikpt) 320 321 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,1,me)) then 322 bdtot_index=bdtot_index+nband_k 323 cycle 324 end if 325 326 ABI_ALLOCATE(gbound,(2*mgfft+8,2)) 327 ABI_ALLOCATE(kg_k,(3,npw_k)) 328 ABI_ALLOCATE(gbound1,(2*mgfft+8,2)) 329 ABI_ALLOCATE(kg1_k,(3,npw1_k)) 330 331 kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 332 call sphereboundary(gbound,istwf_k,kg_k,mgfft,npw_k) 333 334 kg1_k(:,1:npw1_k)=kg1(:,1+ikg1:npw1_k+ikg1) 335 call sphereboundary(gbound1,istwf_k,kg1_k,mgfft,npw1_k) 336 337 ! Loop over bands to fft and square for rho(r) 338 do iband=1,nband_k 339 340 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband,1,me)) cycle 341 342 ! Only treat occupied states 343 if (abs(occ_rbz(iband+bdtot_index))>tol8) then 344 345 ! Build the four components of rho. We use only norm quantities and, so fourwf. 346 347 ABI_ALLOCATE(wfraug_up,(2,n4,n5,n6)) 348 ABI_ALLOCATE(wfraug_down,(2,n4,n5,n6)) 349 ABI_ALLOCATE(wfraug1_up,(2,n4,n5,n6)) 350 ABI_ALLOCATE(wfraug1_down,(2,n4,n5,n6)) 351 ABI_ALLOCATE(cwave0_up,(2,mpw)) 352 ABI_ALLOCATE(cwave0_down,(2,mpw)) 353 ABI_ALLOCATE(cwave1_up,(2,mpw1)) 354 ABI_ALLOCATE(cwave1_down,(2,mpw1)) 355 356 ! EB FR build spinorial wavefunctions 357 ! Obtain Fourier transform in fft box and accumulate the density 358 ! EB FR How do we manage the following lines for non-collinear???? 359 ! zero order up and down spins 360 ptr1 = 1 + (iband-1)*npw_k*nspinor + icg 361 call cg_zcopy(npw_k, cg(1,ptr1), cwave0_up) 362 ptr2 = 1 + npw_k + (iband-1)*npw_k*nspinor + icg 363 call cg_zcopy(npw_k, cg(1,ptr2), cwave0_down) 364 ! first order up and down spins 365 ptr1 = 1 + (iband-1)*npw1_k*nspinor + icg1 366 call cg_zcopy(npw1_k, cg1(1,ptr1), cwave1_up) 367 ptr2 = 1 + npw1_k + (iband-1)*npw1_k*nspinor + icg1 368 call cg_zcopy(npw1_k, cg1(1,ptr2), cwave1_down) 369 370 ! TODO: here ispinor should be ispinorp to get full matrix and nspden 4 371 ! ptr = 1 + (ispinor-1)*npw1_k + (iband-1)*npw1_k*nspinor + icg1 372 ! call cg_zcopy(npw1_k, cg1(1,ptr), cwavef1) 373 374 ! EB FR lines to be managed (?????) 375 376 ! zero order 377 ! cwave0_up => cwavef(:,1:npw_k) 378 ! cwave0_down => cwavef(:,1+npw_k:2*npw_k) 379 ! first order 380 ! cwave1_up => cwavef1(:,1:npw_k) 381 ! cwave1_down => cwavef1(:,1+npw_k:2*npw_k) 382 383 ! The factor 2 is not the spin factor (see Eq.44 of PRB55,10337 (1997)??) 384 weight=two*occ_rbz(iband+bdtot_index)*wtk_rbz(ikpt)/ucvol 385 !density components 386 !GS wfk Fourrier Tranform 387 ! EB FR in the fourwf calls rhoaug(:,:,:,2) is a dummy argument 388 call fourwf(1,rhoaug(:,:,:,2),cwave0_up,dummy,wfraug_up,gbound,gbound,istwf_k,kg_k,kg_k,& 389 & mgfft,mpi_enreg,1,ngfft,npw_k,1,n4,n5,n6,& 390 & 0,paral_kgb,tim_fourwf7,weight,weight) 391 call fourwf(1,rhoaug(:,:,:,2),cwave0_down,dummy,wfraug_down,gbound,gbound,istwf_k,kg_k,kg_k,& 392 & mgfft,mpi_enreg,1,ngfft,npw_k,1,n4,n5,n6,& 393 & 0,paral_kgb,tim_fourwf7,weight,weight) 394 !1st order wfk Fourrier Transform 395 call fourwf(1,rhoaug1(:,:,:,2),cwave1_up,dummy,wfraug1_up,gbound,gbound,istwf_k,kg_k,kg_k,& 396 & mgfft,mpi_enreg,1,ngfft,npw_k,1,n4,n5,n6,& 397 & 0,paral_kgb,tim_fourwf7,weight,weight) 398 call fourwf(1,rhoaug1(:,:,:,2),cwave1_down,dummy,wfraug1_down,gbound,gbound,istwf_k,kg_k,kg_k,& 399 & mgfft,mpi_enreg,1,ngfft,npw_k,1,n4,n5,n6,& 400 & 0,paral_kgb,tim_fourwf7,weight,weight) 401 402 ! Accumulate 1st-order density (x component) 403 if (cplex==2) then 404 do i3=1,n3 405 do i2=1,n2 406 do i1=1,n1 407 re0_up=wfraug_up(1,i1,i2,i3) ; im0_up=wfraug_up(2,i1,i2,i3) 408 re1_up=wfraug1_up(1,i1,i2,i3) ; im1_up=wfraug1_up(2,i1,i2,i3) 409 re0_down=wfraug_down(1,i1,i2,i3) ; im0_down=wfraug_down(2,i1,i2,i3) 410 re1_down=wfraug1_down(1,i1,i2,i3) ; im1_down=wfraug1_down(2,i1,i2,i3) 411 rhoaug1(2*i1-1,i2,i3,1)=rhoaug1(2*i1-1,i2,i3,1)+weight*(re0_up*re1_up+im0_up*im1_up) !n_upup 412 rhoaug1(2*i1 ,i2,i3,1)=zero ! imag part of n_upup at k 413 rhoaug1(2*i1-1,i2,i3,4)=rhoaug1(2*i1-1,i2,i3,4)+weight*(re0_down*re1_down+im0_down*im1_down) ! n_dndn 414 rhoaug1(2*i1 ,i2,i3,4)=zero ! imag part of n_dndn at k 415 rhoaug1(2*i1-1,i2,i3,2)=rhoaug1(2*i1-1,i2,i3,2)+weight*(re1_up*re0_down+re0_up*re1_down & 416 & +im0_up*im1_down+im0_down*im1_up) ! mx; the factor two is inside weight 417 rhoaug1(2*i1 ,i2,i3,2)=zero ! imag part of mx 418 rhoaug1(2*i1-1,i2,i3,3)=rhoaug1(2*i1-1,i2,i3,3)+weight*(re1_up*im0_down-im1_up*re0_down & 419 & +re0_up*im1_down-im0_up*re1_down) ! my; the factor two is inside weight 420 rhoaug1(2*i1 ,i2,i3,3)=zero ! imag part of my at k 421 end do 422 end do 423 end do 424 else 425 do i3=1,n3 426 do i2=1,n2 427 do i1=1,n1 428 re0_up=wfraug_up(1,i1,i2,i3) ; im0_up=wfraug_up(2,i1,i2,i3) 429 re1_up=wfraug1_up(1,i1,i2,i3) ; im1_up=wfraug1_up(2,i1,i2,i3) 430 re0_down=wfraug_down(1,i1,i2,i3) ; im0_down=wfraug_down(2,i1,i2,i3) 431 re1_down=wfraug1_down(1,i1,i2,i3) ; im1_down=wfraug1_down(2,i1,i2,i3) 432 rhoaug1(i1,i2,i3,1)=rhoaug1(i1,i2,i3,1)+weight*(re0_up*re1_up+im0_up*im1_up) ! n_upup 433 rhoaug1(i1,i2,i3,4)=rhoaug1(i1,i2,i3,4)+weight*(re0_down*re1_down+im0_down*im1_down) ! n_dndn 434 rhoaug1(i1,i2,i3,2)=rhoaug1(i1,i2,i3,2)+weight*(re1_up*re0_down+re0_up*re1_down & 435 & +im0_up*im1_down+im0_down*im1_up) !mx; the factor two is inside weight 436 rhoaug1(i1,i2,i3,3)=rhoaug1(i1,i2,i3,3)+weight*(re1_up*im0_down-im1_up*re0_down & 437 & +re0_up*im1_down-im0_up*re1_down) !my; the factor two is inside weight 438 end do 439 end do 440 end do 441 end if 442 ABI_DEALLOCATE(wfraug_up) 443 ABI_DEALLOCATE(wfraug_down) 444 ABI_DEALLOCATE(wfraug1_up) 445 ABI_DEALLOCATE(wfraug1_down) 446 ABI_DEALLOCATE(cwave0_up) 447 ABI_DEALLOCATE(cwave0_down) 448 ABI_DEALLOCATE(cwave1_up) 449 ABI_DEALLOCATE(cwave1_down) 450 451 end if ! occupied states 452 end do ! End loop on iband 453 454 ABI_DEALLOCATE(gbound) 455 ABI_DEALLOCATE(kg_k) 456 ABI_DEALLOCATE(gbound1) 457 ABI_DEALLOCATE(kg1_k) 458 459 bdtot_index=bdtot_index+nband_k 460 461 icg=icg+npw_k*nband_k 462 ikg=ikg+npw_k 463 464 icg1=icg1+npw1_k*nband_k 465 ikg1=ikg1+npw1_k 466 467 end do ! End loop on ikpt 468 469 470 if (xmpi_paral==0) then ! Write the number of one-way 3D ffts skipped until now 471 write(message,'(a,i8)')' dfpt_mkrho : number of one-way 3D ffts skipped in mkrho3 until now =',nskip 472 call wrtout(std_out,message,'PERS') 473 end if 474 475 ! TODO: if n+magnetization basis is used for the density, need to rotate rhoaug1 to that now, before packing into rhor1 476 477 ! Transfer density on augmented fft grid to normal fft grid in real space 478 ! Take also into account the spin, to place it correctly in rhor1. 479 ! Note the use of cplex 480 call fftpac(1,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,ngfft,rhor1,rhoaug1,1) 481 end if ! nspden /= 4 482 483 !if (xmpi_paral==1) then 484 !call timab(63,1,tsec) 485 !call wrtout(std_out,'dfpt_mkrho: loop on k-points and spins done in parallel','COLL') 486 !call xmpi_barrier(spaceworld) 487 !call timab(63,2,tsec) 488 !end if 489 490 ABI_DEALLOCATE(cwavef) 491 ABI_DEALLOCATE(cwavef1) 492 ABI_DEALLOCATE(dummy) 493 ABI_DEALLOCATE(rhoaug) 494 ABI_DEALLOCATE(rhoaug1) 495 ABI_DEALLOCATE(wfraug) 496 ABI_DEALLOCATE(wfraug1) 497 498 !Recreate full rhor1 on all proc. 499 call timab(48,1,tsec) 500 call timab(71,1,tsec) 501 call xmpi_sum(rhor1,spaceworld,ierr) 502 call timab(71,2,tsec) 503 call timab(48,2,tsec) 504 505 call symrhg(cplex,gprimd,irrzon,mpi_enreg,nfft,nfft,ngfft,nspden,nsppol,nsym,paral_kgb,phnons,& 506 & rhog1,rhor1,rprimd,symafm,symrel) 507 508 !We now have both rho(r) and rho(G), symmetrized, and if nsppol=2 509 !we also have the spin-up density, symmetrized, in rhor1(:,2). 510 511 !DBG_EXIT("COLL") 512 513 end subroutine dfpt_mkrho