TABLE OF CONTENTS
ABINIT/mlwfovlp_projpaw [ Functions ]
NAME
mlwfovlp_projpaw
FUNCTION
Calculates the functions that are given to Wannier90 as an starting guess. Here we project them inside the PAW spheres
COPYRIGHT
Copyright (C) 2009-2018 ABINIT group (T. Rangel) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
band_in(mband)= logical array which indicates the bands to be excluded from the calculation
cprj(natom,nspinor*mband*mkmem*nsppol)= <p_lmn|Cnk> coefficients for each WF |Cnk>
and each |p_lmn> non-local projector
just_augmentation= flag used to indicate that we are just going
to compute augmentation part of the matrix
and we are excluding the plane wave part.
mband= maximum number of bands
mkmem= number of k points which can fit in memory; set to 0 if use disk
natom= number of atoms in cell.
nband(nkpt*nsppol)= array cointaining number of bands at each k-point and isppol
nkpt=number of k points.
num_bands=number of bands actually used to construct the wannier function (NOT USED IN 6.7.1 SO WAS TEMPORARILY REMOVED)
nspinor=number of spinorial components of the wavefunctions
nsppol=1 for unpolarized, 2 for spin-polarized
ntypat=number of types of atoms in unit cell.
nwan= number of wannier fonctions (read in wannier90.win).
pawrad(ntypat)= type(pawrad_type) radial information of paw objects
pawtab(ntypat)= For PAW, TABulated data initialized at start
proj_l(mband)= angular part of the projection function (quantum number l)
proj_m(mband)= angular part of the projection function (quantum number m)
proj_radial(mband)= radial part of the projection.
proj_site(3,mband)= site of the projection.
proj_x(3,mband)= x axis for the projection.
proj_z(3,mband)= z axis for the projection.
proj_zona(mband)= extension of the radial part.
psps <type(pseudopotential_type)>=variables related to pseudopotentials
rprimd(3,3)= Direct lattice vectors, Bohr units.
spin = just used for nsppol>1 ; 0 both, 1 just spin up, 2 just spin down
typat(natom)= atom type
xred(3,natom)=reduced dimensionless atomic coordinates
OUTPUT
A_paw(max_num_bands,nwan,nkpt) = A matrix containing initial guess for MLWFs
(augmentation part of the matrix)
SIDE EFFECTS
NOTES
This routine is still under developement
PARENTS
mlwfovlp
CHILDREN
mlwfovlp_ylmfar,simp_gen,simpson_int,wrtout,xred2xcart
SOURCE
66 #if defined HAVE_CONFIG_H 67 #include "config.h" 68 #endif 69 70 #include "abi_common.h" 71 72 subroutine mlwfovlp_projpaw(A_paw,band_in,cprj,just_augmentation,max_num_bands,mband,mkmem,& 73 &mwan,natom,nband,nkpt,& 74 &nspinor,nsppol,ntypat,nwan,pawrad,pawtab,& 75 &proj_l,proj_m,proj_radial,proj_site,proj_x,proj_z,proj_zona,psps,& 76 &rprimd,spin,typat,xred) 77 78 use defs_basis 79 use defs_datatypes 80 use defs_wannier90 81 use m_errors 82 use m_profiling_abi 83 84 use m_numeric_tools, only : simpson_int 85 use m_geometry, only : xred2xcart 86 use m_pawrad, only : pawrad_type, simp_gen 87 use m_pawtab, only : pawtab_type 88 use m_pawcprj, only : pawcprj_type 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 'mlwfovlp_projpaw' 94 use interfaces_14_hidewrite 95 use interfaces_67_common, except_this_one => mlwfovlp_projpaw 96 !End of the abilint section 97 98 implicit none 99 100 !Arguments ------------------------------------ 101 integer,intent(in) :: max_num_bands,mband,mkmem,mwan,natom,nkpt 102 integer,intent(in) :: nspinor,nsppol,ntypat,spin 103 !arrays 104 integer,intent(in) :: nband(nsppol*nkpt),nwan(nsppol) 105 integer,intent(in) :: proj_l(mband,nsppol),proj_m(mband,nsppol),proj_radial(mband,nsppol) 106 integer,intent(in) :: typat(natom) 107 real(dp),intent(in):: proj_site(3,mband,nsppol) 108 real(dp),intent(in) :: proj_x(3,mband,nsppol),proj_z(3,mband,nsppol),proj_zona(mband,nsppol) 109 real(dp),intent(in) :: rprimd(3,3),xred(3,natom) 110 complex(dpc),intent(out) :: A_paw(max_num_bands,mwan,nkpt,nsppol) 111 logical,intent(in) :: band_in(mband,nsppol) 112 logical,intent(in)::just_augmentation(mwan,nsppol) 113 type(pawcprj_type) :: cprj(natom,nspinor*mband*mkmem*nsppol) 114 type(pawrad_type),intent(in) :: pawrad(ntypat) 115 type(pawtab_type),intent(in) :: pawtab(ntypat) 116 type(pseudopotential_type),intent(in) :: psps 117 118 !Local variables------------------------------- 119 !local variables 120 integer :: basis_size,iatom,iband,ii 121 integer :: ikpt,ir,isppol,itypat,iwan,jband 122 integer :: ll,lm,ln,mm,ilmn 123 integer :: lmn_size,max_lmax2, mesh_size,nn 124 integer :: lmax(nsppol),lmax2(nsppol) 125 real(dp):: aa,int_rad2,prod_real,prod_imag 126 real(dp),parameter :: dx=0.015d0,rmax=10.d0,xmin=0.d0 127 real(dp):: sum,wan_lm_fac,x 128 complex(dpc)::prod 129 character(len=500) :: message 130 !arrays 131 integer :: index(mband,nkpt,nsppol) 132 real(dp):: dist,norm(mwan,nsppol) 133 real(dp):: proj_cart(3,mwan,nsppol),proj_site_unit(3,mwan,nsppol) 134 real(dp):: xcart_unit(3,natom),xred_unit(3,natom) 135 real(dp),allocatable ::aux(:),ff(:),r(:),int_rad(:),rad_int(:) 136 real(dp),allocatable::ylmr_fac(:,:,:) 137 138 !no_abirules 139 !Tables 3.1 & 3.2, User guide 140 integer,save :: orb_l_defs(-5:3)=(/2,2,1,1,1,0,1,2,3/) 141 !real(dp),allocatable :: ylm(:,:) 142 143 144 ! ************************************************************************* 145 146 !DEBUG 147 !write (std_out,*) ' mlwfovlp_projpaw : enter' 148 !ENDDEBUG 149 150 !DEBUG ! to be uncommented, if needed 151 152 write(message, '(a,a)' )ch10,& 153 & '** mlwfovlp_proj: compute in-sphere part of A_matrix' 154 call wrtout(std_out,message,'COLL') 155 156 ! 157 !Check input variables 158 ! 159 do isppol=1,nsppol 160 if(spin.ne.0 .and. spin.ne.isppol) cycle 161 do iwan=1,nwan(nsppol) 162 if(proj_radial(iwan,isppol)<1 .or. proj_radial(iwan,isppol)>4)then 163 write(message,'(a,a,a,i6)')& 164 & ' proj_radial should be between 1 and 4,',ch10,& 165 & ' however, proj_radial=',proj_radial(iwan,isppol) 166 MSG_BUG(message) 167 end if 168 end do 169 end do 170 171 ! 172 !Initialize 173 ! 174 A_paw(:,:,:,:)=cmplx(0.d0,0.d0) 175 ! 176 !Get index for cprj 177 ! 178 ii=0 179 do isppol=1,nsppol 180 do ikpt=1,nkpt 181 do iband=1,nband(ikpt) 182 ii=ii+1 183 index(iband,ikpt,isppol)=ii 184 end do 185 end do 186 end do 187 ! 188 !obtain lmax and lmax2 189 ! 190 lmax(:)=0 191 lmax2(:)=0 192 do isppol=1,nsppol 193 if(spin.ne.0 .and. spin.ne.isppol) cycle 194 do iwan=1,nwan(isppol) 195 lmax(isppol)=max(lmax(isppol),orb_l_defs(proj_l(iwan,isppol))) 196 end do !iwan 197 lmax2(isppol)=(lmax(isppol)+1)**2 198 end do 199 max_lmax2=maxval(lmax2(:)) 200 ! 201 !get ylmfac, factor used for rotations and hybrid orbitals 202 ! 203 ABI_ALLOCATE(ylmr_fac,(max_lmax2,mwan,nsppol)) 204 205 206 do isppol=1,nsppol 207 if(spin.ne.0 .and. spin.ne.isppol) cycle 208 call mlwfovlp_ylmfar(ylmr_fac(1:lmax2(isppol),1:nwan(isppol),isppol),& 209 & lmax(isppol),lmax2(isppol),mband,nwan(isppol),proj_l(:,isppol),proj_m(:,isppol),& 210 & proj_x(:,:,isppol),proj_z(:,:,isppol)) 211 ! 212 ! Shift projection centers and atom centers to the primitive cell 213 ! This will be useful after, when we check if the Wannier function 214 ! lies on one specific atom 215 ! 216 proj_site_unit(:,:,:)=0.d0 217 do iwan=1,nwan(isppol) 218 do ii=1,3 219 proj_site_unit(ii,iwan,isppol)=ABS(proj_site(ii,iwan,isppol)-AINT(proj_site(ii,iwan,isppol)) ) 220 end do 221 end do 222 do iatom=1,natom 223 do ii=1,3 224 xred_unit(ii,iatom)=ABS(xred(ii,iatom)-AINT(xred(ii,iatom)) ) 225 end do 226 end do 227 call xred2xcart(natom,rprimd,xcart_unit,xred_unit) 228 call xred2xcart(mwan,rprimd,proj_cart(:,:,isppol),proj_site_unit(:,:,isppol)) 229 ! 230 ! Normalize the Wannier functions 231 ! 232 ! Radial part 233 mesh_size= nint((rmax - xmin ) / dx + 1) 234 ABI_ALLOCATE( ff,(mesh_size)) 235 ABI_ALLOCATE(r,(mesh_size)) 236 ABI_ALLOCATE(rad_int,(mesh_size)) 237 ABI_ALLOCATE(aux,(mesh_size)) 238 do ir=1, mesh_size 239 x=xmin+DBLE(ir-1)*dx 240 r(ir)=x 241 end do !ir 242 do iwan=1,nwan(isppol) 243 ! write(std_out,*)'iwan',iwan 244 ! radial functions shown in table 3.3 of wannier90 manual 245 if(proj_radial(iwan,isppol)==1) ff(:) = 2.d0 * proj_zona(iwan,isppol)**(1.5d0) * exp(-proj_zona(iwan,isppol)*r(:)) 246 if(proj_radial(iwan,isppol)==2) ff(:) = 1.d0/(2.d0*sqrt(2.d0))*proj_zona(iwan,isppol)**(1.5d0) *& 247 & (2.d0 - proj_zona(iwan,isppol)*r(:))*exp(-proj_zona(iwan,isppol)*r(:)/2.d0) 248 if(proj_radial(iwan,isppol)==3) ff(:) = sqrt(4.d0/27.d0)*proj_zona(iwan,isppol)**(1.5d0)& 249 & * (1.d0 - 2.d0*proj_zona(iwan,isppol)*r(:)/3.d0 + 2.d0*proj_zona(iwan,isppol)**2*r(:)**2/27.d0)& 250 & * exp(-proj_zona(iwan,isppol) * r(:)/3.d0) 251 252 if(proj_radial(iwan,isppol)/=4) then 253 aux(:)=ff(:)**2*r(:)**2 254 call simpson_int(mesh_size,dx,aux,rad_int) 255 sum=0.d0 256 do ir=1,mesh_size 257 sum=sum+rad_int(ir) 258 end do 259 int_rad2=sum/real(mesh_size,dp) 260 ! 261 ! do ir=1,mesh_size 262 ! if(iwan==1) write(400,*)r(ir),aux(ir),rad_int(ir) 263 ! end do 264 else 265 ! 266 ! ==4: gaussian function 267 ! f(x)=\exp(-1/4(x/aa)**2) 268 ! \int f(x)f(x) dx = \int \exp(-1/2(x/aa)**2) = aa*sqrt(2pi) 269 ! 270 int_rad2=sqrt(2.d0*pi)*proj_zona(iwan,isppol) 271 end if 272 273 ! 274 ! Now angular part 275 ! 276 prod_real=0.d0 277 do lm=1,lmax2(isppol) 278 wan_lm_fac=ylmr_fac(lm,iwan,isppol) 279 ! write(std_out,*)'wan_lm_fac',wan_lm_fac 280 ! write(std_out,*)'int_rad2',int_rad2 281 prod_real= prod_real + wan_lm_fac**2 * int_rad2 282 end do 283 norm(iwan,isppol)=sqrt(prod_real) 284 end do !iwan 285 ABI_DEALLOCATE(ff) 286 ABI_DEALLOCATE(r) 287 ABI_DEALLOCATE(rad_int) 288 ABI_DEALLOCATE(aux) 289 ! 290 ! Now that we found our guiding functions 291 ! We proceed with the internal product of 292 ! our guiding functions and the wave function 293 ! Amn=<G_m|\Psi_n> inside the sphere. 294 ! The term <G_m|\Psi_n> inside the sphere is: 295 ! = \sum_i <G_n | \phi_i - \tphi_i> <p_im|\Psi_m> 296 ! 297 ! 298 ! G_n \phi and \tphi can be decomposed in 299 ! a radial function times an angular function. 300 ! 301 ! 302 ! Big loop on iwan and iatom 303 ! 304 do iwan=1,nwan(isppol) 305 do iatom=1,natom 306 ! 307 ! check if center of wannier function coincides 308 ! with the center of the atom 309 ! 310 dist=((proj_cart(1,iwan,isppol)-xcart_unit(1,iatom))**2 + & 311 & (proj_cart(2,iwan,isppol)-xcart_unit(2,iatom))**2 + & 312 & (proj_cart(3,iwan,isppol)-xcart_unit(3,iatom))**2)**0.5 313 ! 314 ! if the distance between the centers is major than 0.1 angstroms skip 315 ! 316 if( dist > 0.188972613) cycle 317 ! 318 write(message, '(2a,i4,a,i4,2a)')ch10, ' Wannier function center',iwan,' is on top of atom',& 319 & iatom,ch10,' Calculating in-sphere contribution' 320 call wrtout(ab_out,message,'COLL') 321 call wrtout(std_out, message,'COLL') 322 ! 323 ! Get useful quantities 324 ! 325 itypat=typat(iatom) 326 lmn_size=pawtab(itypat)%lmn_size 327 basis_size=pawtab(itypat)%basis_size 328 mesh_size=pawtab(itypat)%mesh_size 329 ABI_ALLOCATE(int_rad,(basis_size)) 330 ABI_ALLOCATE(ff,(mesh_size)) 331 ABI_ALLOCATE(aux,(mesh_size)) 332 333 ! 334 ! Integrate first the radial part 335 ! and save it into an array 336 ! 337 ! 338 ! radial functions shown in table 3.3 of wannier90 manual 339 ! 340 if(proj_radial(iwan,isppol)==1) aux(1:mesh_size) = 2.d0 * proj_zona(iwan,isppol)**(1.5d0) *& 341 & exp(-proj_zona(iwan,isppol)*pawrad(itypat)%rad(1:mesh_size)) 342 if(proj_radial(iwan,isppol)==2) aux(1:mesh_size) = 1.d0/(2.d0*sqrt(2.d0))*proj_zona(iwan,isppol)**(1.5d0) *& 343 & (2.d0 - proj_zona(iwan,isppol)*pawrad(itypat)%rad(1:mesh_size)) & 344 & * exp(-proj_zona(iwan,isppol)*pawrad(itypat)%rad(1:mesh_size)/2.d0) 345 if(proj_radial(iwan,isppol)==3) aux(1:mesh_size) = sqrt(4.d0/27.d0)*proj_zona(iwan,isppol)**(1.5d0)& 346 & * (1.d0 - 2.d0*proj_zona(iwan,isppol)*pawrad(itypat)%rad(1:mesh_size)/3.d0 & 347 & + 2.d0*proj_zona(iwan,isppol)**2 *pawrad(itypat)%rad(1:mesh_size)**2/27.d0)& 348 & * exp(-proj_zona(iwan,isppol) * pawrad(itypat)%rad(1:mesh_size)/3.d0) 349 ! 350 ! ==4: gaussian function 351 ! f(x)=\exp(-1/4(x/aa)**2) 352 ! 353 if(proj_radial(iwan,isppol)==4) then 354 aa=1.d0/proj_zona(iwan,isppol) 355 aux(1:mesh_size)= exp(-0.25d0*(pawrad(itypat)%rad(1:mesh_size)*aa)**2) 356 end if 357 ! 358 ! Normalize aux 359 aux(:)=aux(:)/norm(iwan,isppol) 360 ! 361 do ln=1,basis_size 362 if(just_augmentation(iwan,isppol)) then 363 ! 364 ! just augmentation region contribution 365 ! In this case there is no need to use \tphi 366 ! ff= \int R_wan(r) (R_phi(ln;r)/r ) r^2 dr 367 ! 368 ff(1:mesh_size)= aux(1:mesh_size) * pawtab(itypat)%phi(1:mesh_size,ln) & 369 & * pawrad(itypat)%rad(1:mesh_size) 370 else 371 ! Inside sphere contribution = \phi - \tphi 372 ! ff= \int R_wan(r) (R_phi(ln;r)/r - R_tphi(ln;r)/r) r^2 dr 373 ff(1:mesh_size)= aux(1:mesh_size) * (pawtab(itypat)%phi(1:mesh_size,ln)-pawtab(itypat)%tphi(1:mesh_size,ln)) & 374 & * pawrad(itypat)%rad(1:mesh_size) 375 end if 376 ! 377 ! Integration with simpson routine 378 ! 379 call simp_gen(int_rad(ln),ff,pawrad(itypat)) 380 ! do ii=1,mesh_size 381 ! unit_ln=400+ln 382 ! if( iwan==1 ) write(unit_ln,*)pawrad(itypat)%rad(ii),ff(ii),int_rad(ln) 383 ! end do 384 end do !ln 385 ABI_DEALLOCATE(ff) 386 ABI_DEALLOCATE(aux) 387 ! 388 ! Now integrate the angular part 389 ! Cycle on i indices 390 ! 391 ! prod_real=0.d0 392 do ilmn=1, lmn_size 393 ll=Psps%indlmn(1,ilmn,itypat) 394 mm=Psps%indlmn(2,ilmn,itypat) 395 nn=Psps%indlmn(3,ilmn,itypat) 396 lm=Psps%indlmn(4,ilmn,itypat) 397 ln=Psps%indlmn(5,ilmn,itypat) 398 ! write(std_out,*)'ll ',ll,' mm ',mm,'nn',nn,"lm",lm,"ln",ln 399 ! 400 ! Get wannier factor for that lm component 401 if(lm <=lmax2(isppol)) then 402 wan_lm_fac=ylmr_fac(lm,iwan,isppol) 403 ! Make delta product 404 ! Here we integrate the angular part 405 ! Since the integral of the product of two spherical harmonics 406 ! is a delta function 407 if( abs(wan_lm_fac) > 0.0d0) then 408 ! write(std_out,*) 'll',ll,'mm',mm,'lm',lm,'ln',ln,'factor',wan_lm_fac !lm index for wannier function 409 ! 410 ! Calculate Amn_paw, now that the radial and angular integrations are done 411 ! 412 prod=cmplx(0.d0,0.d0) 413 do ikpt=1,nkpt 414 jband=0 415 do iband=1,nband(ikpt) 416 if(band_in(iband,isppol)) then 417 jband=jband+1 418 419 prod_real= cprj(iatom,index(iband,ikpt,isppol))%cp(1,ilmn) * int_rad(ln) * wan_lm_fac 420 prod_imag= cprj(iatom,index(iband,ikpt,isppol))%cp(2,ilmn) * int_rad(ln) * wan_lm_fac 421 prod=cmplx(prod_real,prod_imag) 422 423 A_paw(jband,iwan,ikpt,isppol)=A_paw(jband,iwan,ikpt,isppol)+prod 424 end if !band_in 425 end do !iband 426 end do !ikpt 427 ! 428 end if !lm<=lmax2 429 end if ! abs(wan_lm_fac) > 0.0d0 430 end do !ilmn=1, lmn_size 431 ABI_DEALLOCATE(int_rad) 432 end do !iatom 433 end do !iwan 434 end do !isppol 435 ! 436 !Deallocate quantities 437 ! 438 ABI_DEALLOCATE(ylmr_fac) 439 440 441 !DEBUG 442 !write (std_out,*) ' mlwfovlp_projpaw : exit' 443 !stop 444 !ENDDEBUG 445 446 end subroutine mlwfovlp_projpaw