TABLE OF CONTENTS
ABINIT/test_charge [ Functions ]
NAME
test_charge
FUNCTION
Reports info on the electronic charge as well as Drude plasma frequency. Mainly used in the GW part.
INPUTS
nelectron_exp=Expected total number of electrons (used to normalize the charge)
OUTPUT
PARENTS
bethe_salpeter,mrgscr,screening,sigma
CHILDREN
wrtout
SOURCE
314 subroutine test_charge(nfftf,nelectron_exp,nspden,rhor,ucvol,& 315 & usepaw,usexcnhat,usefinegrid,compch_sph,compch_fft,omegaplasma) 316 317 use defs_basis 318 use m_profiling_abi 319 320 !This section has been created automatically by the script Abilint (TD). 321 !Do not modify the following lines by hand. 322 #undef ABI_FUNC 323 #define ABI_FUNC 'test_charge' 324 use interfaces_14_hidewrite 325 !End of the abilint section 326 327 implicit none 328 329 !Arguments ------------------------------------ 330 !scalars 331 integer,intent(in) :: nfftf,nspden,usefinegrid,usepaw,usexcnhat 332 real(dp),intent(in) :: compch_fft,compch_sph,ucvol,nelectron_exp 333 real(dp),intent(out) :: omegaplasma 334 !arrays 335 real(dp),intent(inout) :: rhor(nfftf,nspden) 336 337 !Local variables ------------------------------ 338 !scalars 339 real(dp) :: nelectron_tot,nelectron_fft 340 real(dp) :: nelectron_pw,nelectron_sph,rhoav,rs,nratio 341 character(len=500) :: msg 342 343 !************************************************************************* 344 345 ! ABI_UNUSED(usexcnhat) 346 if (usexcnhat==0)then 347 end if 348 349 ! === For PAW output of compensation charges === 350 if (usepaw==1) then 351 !if (usepaw==1.and.usexcnhat>0) then ! TODO I still dont understand this if! 352 write(msg,'(4a)')ch10,' PAW TEST:',ch10,' ==== Compensation charge inside spheres ============' 353 if (compch_sph<greatest_real.and.compch_fft<greatest_real) & 354 & write(msg,'(3a)')TRIM(msg),ch10,' The following values must be close...' 355 if (compch_sph<greatest_real) & 356 & write(msg,'(3a,f22.15)')TRIM(msg),ch10,' Compensation charge over spherical meshes = ',compch_sph 357 if (compch_fft<greatest_real) then 358 if (usefinegrid==1) then 359 write(msg,'(3a,f22.15)')TRIM(msg),ch10,' Compensation charge over fine fft grid = ',compch_fft 360 else 361 write(msg,'(3a,f22.15)')TRIM(msg),ch10,' Compensation charge over fft grid = ',compch_fft 362 end if 363 end if 364 call wrtout(ab_out,msg,'COLL') 365 call wrtout(std_out,msg,'COLL') 366 write(msg,'(a)')ch10 367 call wrtout(ab_out,msg,'COLL') 368 call wrtout(std_out,msg,'COLL') 369 end if !PAW 370 371 nelectron_pw =SUM(rhor(:,1))*ucvol/nfftf 372 nelectron_tot=nelectron_pw 373 nratio =nelectron_exp/nelectron_tot 374 375 if (usepaw==1) then 376 nelectron_sph=nelectron_pw+compch_sph 377 nelectron_fft=nelectron_pw+compch_fft 378 nelectron_tot=nelectron_sph 379 nratio=(nelectron_exp-nelectron_sph)/nelectron_pw 380 end if 381 382 rhoav=nelectron_tot/ucvol ; rs=(three/(four_pi*rhoav))**third 383 if (usepaw==0) then 384 write(msg,'(2(a,f9.4))')& 385 & ' Number of electrons calculated from density = ',nelectron_tot,'; Expected = ',nelectron_exp 386 else 387 write(msg,'(2(a,f9.4),a)')& 388 & ' Total number of electrons per unit cell = ',nelectron_sph,' (Spherical mesh), ',nelectron_fft,' (FFT mesh)' 389 end if 390 call wrtout(std_out,msg,'COLL') 391 call wrtout(ab_out,msg,'COLL') 392 393 !$write(msg,'(a,f9.4)')' Renormalizing smooth charge density using nratio = ',nratio 394 !! rhor(:,:)=nratio*rhor(:,:) 395 396 write(msg,'(a,f9.6)')' average of density, n = ',rhoav 397 call wrtout(std_out,msg,'COLL') 398 call wrtout(ab_out,msg,'COLL') 399 write(msg,'(a,f9.4)')' r_s = ',rs 400 call wrtout(std_out,msg,'COLL') 401 call wrtout(ab_out,msg,'COLL') 402 omegaplasma=SQRT(four_pi*rhoav) 403 write(msg,'(a,f9.4,2a)')' omega_plasma = ',omegaplasma*Ha_eV,' [eV]',ch10 404 call wrtout(std_out,msg,'COLL') 405 call wrtout(ab_out,msg,'COLL') 406 407 end subroutine test_charge
ABINIT/wfd_mkrho [ Functions ]
NAME
wfd_mkrho
FUNCTION
Calculate the charge density on the fine FFT grid in real space.
COPYRIGHT
Copyright (C) 1999-2018 ABINIT group (MG,GMR, VO, LR, RWG, RShaltaf) 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
ngfftf(18)=array containing all the information for the "fine" FFT. Cryst<crystal_t> Info on the crystalline structure optcalc=option for calculation. If =0 (default value) perform calculation of electronic density. If =1, perform calculation of kinetic energy density. In both cases, the result is returned in rhor. Psps<type(pseudopotential_type)>=variables related to pseudopotentials nfftf=Total number of points on the fine FFT grid (for this processor) Kmesh<kmesh_t>= Info on the k-sampling: Wfd<wfd_t)=datatype gathering info on the wavefunctions. [optcalc]=Optional option used to calculate the kinetic energy density. Defaults to 0.
OUTPUT
rhor(nfftf,nspden)=The density in the real space on the fine FFT grid. If nsppol==2, total charge in first half, spin-up component in second half. (for non-collinear magnetism, first element: total density, 3 next ones: mx,my,mz in units of hbar/2) If optcalc==1 (optional argument, default value is 0), then rhor will actually contains kinetic energy density (taur) instead of electronic density.
NOTES
In the case of PAW calculations:
All computations are done on the fine FFT grid.
All variables (nfftf,ngfftf,mgfftf) refer to this fine FFT grid.
All arrays (densities/potentials...) are computed on this fine FFT grid.
Developers have to be careful when introducing others arrays:
they have to be stored on the fine FFT grid.
In the case of norm-conserving calculations:
The mesh is the usual augmented FFT grid to treat correctly the convolution.
PARENTS
bethe_salpeter,screening,sigma
CHILDREN
wrtout
SOURCE
53 #if defined HAVE_CONFIG_H 54 #include "config.h" 55 #endif 56 57 #include "abi_common.h" 58 59 subroutine wfd_mkrho(Wfd,Cryst,Psps,Kmesh,Bands,ngfftf,nfftf,rhor,& 60 & optcalc) ! optional arguments 61 62 63 use defs_basis 64 use defs_datatypes 65 use m_profiling_abi 66 use m_xmpi 67 use m_errors 68 69 use m_fstrings, only : sjoin, itoa 70 use m_iterators, only : iter2_t, iter_yield, iter_len, iter_free 71 use m_crystal, only : crystal_t 72 use m_bz_mesh, only : kmesh_t 73 use m_fft, only : fft_ug 74 use m_wfd, only : wfd_t, wfd_get_ur, wfd_update_bkstab, wfd_iterator_bks, wfd_change_ngfft 75 76 !This section has been created automatically by the script Abilint (TD). 77 !Do not modify the following lines by hand. 78 #undef ABI_FUNC 79 #define ABI_FUNC 'wfd_mkrho' 80 use interfaces_14_hidewrite 81 use interfaces_56_recipspace 82 use interfaces_67_common 83 !End of the abilint section 84 85 implicit none 86 87 !Arguments ------------------------------------ 88 !scalars 89 integer,intent(in) :: nfftf 90 integer,intent(in),optional :: optcalc 91 type(ebands_t),intent(in) :: Bands 92 type(kmesh_t),intent(in) :: Kmesh 93 type(crystal_t),intent(in) :: Cryst 94 type(Pseudopotential_type),intent(in) :: Psps 95 type(wfd_t),intent(inout) :: Wfd 96 !arrays 97 integer,intent(in) :: ngfftf(18) 98 real(dp),intent(out) :: rhor(nfftf,Wfd%nspden) 99 100 !Local variables ------------------------------ 101 !scalars 102 integer,parameter :: ndat1=1 103 integer :: cplex,ib,ib_iter,ierr,ik,ir,is,n1,n2,n3,nfftotf 104 integer :: alpha,nalpha,ipw,myoptcalc 105 real(dp) :: kpt_cart,kg_k_cart,gp2pi1,gp2pi2,gp2pi3,cwftmp,bks_weight 106 character(len=500) :: msg 107 !arrays 108 integer,allocatable :: irrzon(:,:,:) 109 real(dp),allocatable :: phnons(:,:,:),rhog(:,:),rhor_down(:),rhor_mx(:),rhor_my(:),cwavef(:,:) 110 complex(dpc),allocatable :: wfr_x(:),wfr_y(:) 111 complex(gwpc),allocatable :: gradug(:),work(:) 112 complex(gwpc),allocatable,target :: wfr(:) 113 complex(gwpc), ABI_CONTIGUOUS pointer :: cwavef1(:),cwavef2(:) 114 type(iter2_t) :: Iter_bks 115 116 !************************************************************************* 117 118 DBG_ENTER("COLL") 119 120 ! Consistency check. 121 ABI_CHECK(Wfd%nsppol == Bands%nsppol, "Mismatch in nsppol") 122 123 if (ANY(ngfftf(1:3) /= Wfd%ngfft(1:3))) call wfd_change_ngfft(Wfd,Cryst,Psps,ngfftf) 124 125 ! Calculate IBZ contribution to the charge density. 126 ABI_MALLOC(wfr, (nfftf*Wfd%nspinor)) 127 128 if (wfd%nspden == 4) then 129 ABI_MALLOC(wfr_x, (nfftf)) 130 ABI_MALLOC(wfr_y, (nfftf)) 131 ABI_MALLOC(rhor_down, (nfftf)) 132 ABI_MALLOC(rhor_mx, (nfftf)) 133 ABI_MALLOC(rhor_my, (nfftf)) 134 rhor_down = zero; rhor_mx = zero; rhor_my = zero 135 end if 136 137 ! Update the (b,k,s) distribution table. 138 call wfd_update_bkstab(Wfd) 139 140 ! Calculate the unsymmetrized density. 141 rhor=zero 142 myoptcalc=0; if (present(optcalc)) myoptcalc=optcalc 143 nalpha=1; if (myoptcalc==1) nalpha=3 144 if (myoptcalc == 1 .and. wfd%nspinor == 2) then 145 MSG_ERROR("kinetic energy density with nspinor == 2 not implemented") 146 end if 147 148 ! Build the iterator that will distribute the states in an automated way. 149 Iter_bks = wfd_iterator_bks(Wfd,bks_mask=ABS(Bands%occ)>=tol8) 150 151 do alpha=1,nalpha 152 do is=1,Wfd%nsppol 153 do ik=1,Wfd%nkibz 154 do ib_iter=1,iter_len(Iter_bks,ik,is) 155 ib = iter_yield(Iter_bks,ib_iter,ik,is) 156 bks_weight = Bands%occ(ib,ik,is) * Kmesh%wt(ik) / Cryst%ucvol 157 158 call wfd_get_ur(Wfd,ib,ik,is,wfr) 159 160 cwavef1 => wfr(1:nfftf) 161 if (myoptcalc == 1) then 162 ABI_MALLOC(gradug,(Wfd%Kdata(ik)%npw)) 163 ABI_MALLOC(cwavef,(2,Wfd%Kdata(ik)%npw)) 164 ABI_MALLOC(work,(nfftf)) 165 cwavef(1,:)= REAL(Wfd%Wave(ib,ik,is)%ug(:)) 166 cwavef(2,:)=AIMAG(Wfd%Wave(ib,ik,is)%ug(:)) 167 ! Multiplication by 2pi i (k+G)_alpha 168 gp2pi1=Cryst%gprimd(alpha,1)*two_pi 169 gp2pi2=Cryst%gprimd(alpha,2)*two_pi 170 gp2pi3=Cryst%gprimd(alpha,3)*two_pi 171 kpt_cart=gp2pi1*Wfd%kibz(1,ik)+gp2pi2*Wfd%kibz(2,ik)+gp2pi3*Wfd%kibz(3,ik) 172 do ipw=1,Wfd%Kdata(ik)%npw 173 kg_k_cart= gp2pi1*Wfd%Kdata(ik)%kg_k(1,ipw) + & 174 & gp2pi2*Wfd%Kdata(ik)%kg_k(2,ipw) + & 175 & gp2pi3*Wfd%Kdata(ik)%kg_k(3,ipw)+kpt_cart 176 ! ipwsp=ipw!+(ispinor-1)*Wfd%Kdata(ik)%npw 177 cwftmp=-cwavef(2,ipw)*kg_k_cart 178 cwavef(2,ipw)=cwavef(1,ipw)*kg_k_cart 179 cwavef(1,ipw)=cwftmp 180 end do 181 gradug(:)=CMPLX(cwavef(1,:),cwavef(2,:),gwpc) 182 call fft_ug(Wfd%npwarr(ik),nfftf,Wfd%nspinor,ndat1,Wfd%mgfft,Wfd%ngfft,& 183 & Wfd%istwfk(ik),Wfd%Kdata(ik)%kg_k,Wfd%Kdata(ik)%gbound,gradug,work) 184 cwavef1(:)=work(:) 185 ABI_FREE(work) 186 ABI_FREE(cwavef) 187 ABI_FREE(gradug) 188 end if 189 190 !$OMP PARALLEL DO 191 do ir=1,nfftf 192 rhor(ir,is) = rhor(ir,is) + CONJG(cwavef1(ir)) * cwavef1(ir) * bks_weight 193 end do 194 !call cplx_addtorho(n1,n2,n3,n4,n5,n6,ndat,weight_r,ur,rho) 195 196 if (wfd%nspinor == 2 .and. wfd%nspden == 1) then 197 cwavef2 => wfr(1+nfftf:2*nfftf) 198 do ir=1,nfftf 199 rhor(ir, 1) = rhor(ir, 1) + CONJG(cwavef2(ir)) * cwavef2(ir) * bks_weight 200 end do 201 end if 202 203 if (wfd%nspinor == 2 .and. wfd%nspden == 4) then 204 cwavef2 => wfr(1+nfftf:2*nfftf) 205 wfr_x(:) = cwavef1(:) + cwavef2(:) ! $(\Psi^{1}+\Psi^{2})$ 206 wfr_y(:) = cwavef1(:) -j_dpc*cwavef2(:) ! $(\Psi^{1}-i\Psi^{2})$ 207 !$OMP PARALLEL DO 208 do ir=1,nfftf 209 rhor_down(ir) = rhor_down(ir) + CONJG(cwavef2(ir)) * cwavef2(ir) * bks_weight 210 rhor_mx(ir) = rhor_mx(ir) + CONJG(wfr_x(ir)) * wfr_x(ir) * bks_weight 211 rhor_my(ir) = rhor_my(ir) + CONJG(wfr_y(ir)) * wfr_y(ir) * bks_weight 212 end do 213 end if 214 215 end do 216 end do 217 end do 218 219 end do ! enddo alpha 220 221 call iter_free(Iter_bks) 222 223 select case (myoptcalc) 224 case (0) 225 ! density 226 if (wfd%nspden == 4) then 227 rhor(:, 2) = rhor_mx 228 rhor(:, 3) = rhor_my 229 rhor(:, 4) = rhor_down 230 end if 231 case (1) 232 ! convention for taur = 1/2 Sum_i |grad phi_i|^2 233 rhor(:,:)=half*rhor(:,:) 234 235 case default 236 MSG_ERROR(sjoin("Wrong myoptcalc:", itoa(myoptcalc))) 237 end select 238 239 call xmpi_sum(rhor,Wfd%comm,ierr) 240 241 ! Symmetrization in G-space implementing also the AFM case 242 n1=ngfftf(1); n2=ngfftf(2); n3=ngfftf(3); nfftotf=n1*n2*n3 243 244 ABI_MALLOC(irrzon,(nfftotf**(1-1/Cryst%nsym),2,(Wfd%nspden/Wfd%nsppol)-3*(Wfd%nspden/4))) 245 ABI_MALLOC(phnons,(2,nfftotf,(Wfd%nspden/Wfd%nsppol)-3*(Wfd%nspden/4))) 246 247 if (Cryst%nsym/=1) then 248 call irrzg(irrzon,Wfd%nspden,Wfd%nsppol,Cryst%nsym,n1,n2,n3,phnons,Cryst%symafm,Cryst%symrel,Cryst%tnons) 249 end if 250 251 ! Symmetrize rho(r), and pack nspden components following abinit conventions. 252 cplex=1 253 ABI_MALLOC(rhog,(2,cplex*nfftf)) 254 255 call symrhg(cplex,Cryst%gprimd,irrzon,Wfd%MPI_enreg,nfftf,nfftotf,ngfftf,Wfd%nspden,Wfd%nsppol,& 256 & Cryst%nsym,Wfd%paral_kgb,phnons,rhog,rhor,Cryst%rprimd,Cryst%symafm,Cryst%symrel) 257 258 ABI_FREE(rhog) 259 ABI_FREE(phnons) 260 ABI_FREE(irrzon) 261 262 ! Find and print minimum and maximum total electron density 263 ! (or total kinetic energy density, or total element of kinetic energy density tensor) and locations 264 !call wrtout(std_out,'mkrho: echo density (plane-wave part only)','COLL') 265 !call prtrhomxmn(std_out,wfd%mpi_enreg,nfftf,ngfftf,wfd%nspden,1,rhor,optrhor=optcalc,ucvol=crystl%ucvol) 266 267 write(msg,'(a,f9.4)')' planewave contribution to nelect: ',SUM(rhor(:,1))*Cryst%ucvol/nfftf 268 call wrtout(std_out,msg,'COLL') 269 270 if (Wfd%nspden==4) then 271 write(msg,'(a,3f9.4)')& 272 ' mx, my, mz: ',SUM(rhor(:,2))*Cryst%ucvol/nfftf,SUM(rhor(:,3))*Cryst%ucvol/nfftf,SUM(rhor(:,4))*Cryst%ucvol/nfftf 273 call wrtout(std_out,msg,'COLL') 274 end if 275 276 ABI_FREE(wfr) 277 278 if (Wfd%nspden == 4) then 279 ABI_FREE(wfr_x) 280 ABI_FREE(wfr_y) 281 ABI_FREE(rhor_down) 282 ABI_FREE(rhor_mx) 283 ABI_FREE(rhor_my) 284 end if 285 286 DBG_EXIT("COLL") 287 288 end subroutine wfd_mkrho