TABLE OF CONTENTS
- ABINIT/m_mlwfovlp
- m_mlwfovlp/mlwfovlp
- m_mlwfovlp/mlwfovlp_proj
- m_mlwfovlp/mlwfovlp_projpaw
- m_mlwfovlp/mlwfovlp_pw
- m_mlwfovlp/mlwfovlp_radial
- m_mlwfovlp/mlwfovlp_seedname
- m_mlwfovlp/mlwfovlp_setup
- m_mlwfovlp/mlwfovlp_ylmfac
- m_mlwfovlp/mlwfovlp_ylmfar
- mlwfovlp/read_chkunit
ABINIT/m_mlwfovlp [ Modules ]
NAME
m_mlwfovlp
FUNCTION
Interface with Wannier90
COPYRIGHT
Copyright (C) 2005-2024 ABINIT group (BAmadon, CEspejo, FJollet, TRangel, DRH) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_mlwfovlp 23 24 use defs_basis 25 use defs_wannier90 26 use m_abicore 27 use m_errors 28 use m_atomdata 29 use m_xmpi 30 use m_sort 31 #ifdef FC_NAG 32 use f90_unix_dir 33 #endif 34 #ifdef HAVE_NETCDF 35 use netcdf 36 #endif 37 use m_nctk 38 use m_hdr 39 use m_dtset 40 use m_dtfil 41 42 use defs_datatypes, only : pseudopotential_type, ebands_t 43 use defs_abitypes, only : MPI_type 44 use m_io_tools, only : delete_file, get_unit, open_file 45 use m_hide_lapack, only : matrginv 46 use m_fstrings, only : strcat, sjoin 47 use m_numeric_tools, only : uniformrandom, simpson_int, c2r, l2int 48 use m_special_funcs, only : besjm 49 use m_geometry, only : xred2xcart, rotmat, wigner_seitz 50 use m_fftcore, only : sphereboundary 51 use m_crystal, only : crystal_t 52 use m_ebands, only : ebands_ncwrite 53 use m_pawang, only : pawang_type 54 use m_pawrad, only : pawrad_type, simp_gen 55 use m_pawtab, only : pawtab_type 56 use m_pawcprj, only : pawcprj_type 57 use m_paw_sphharm, only : ylm_cmplx, initylmr 58 use m_paw_overlap, only : smatrix_pawinit 59 use m_evdw_wannier, only : evdw_wannier 60 use m_fft, only : fourwf 61 62 implicit none 63 64 private
m_mlwfovlp/mlwfovlp [ Functions ]
[ Top ] [ m_mlwfovlp ] [ Functions ]
NAME
mlwfovlp
FUNCTION
Routine which computes overlap M_{mn}(k,b) and projection A_{mn}(k)
for Wannier code (www.wannier.org f90 version).
Various file are written (wannier90.*) which can be used to run a
separate wannier calculation with the wannier90 code.
INPUTS
crystal<crystal_t>=Info on the crystalline structure. ebands<ebands_t>=The object describing the band structure. hdr <type(hdr_type)>=the header of wf, den and pot files atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f) cg(2,mcg)=planewave coefficients of wavefunctions. cprj(natom,mcprj)= <p_lmn|Cnk> coefficients for each WF |Cnk> and each |p_lmn> non-local projector dtset <type(dataset_type)>=all input variables for this dataset dtfil <type(datafiles_type)>=variables related to files ecut=cut-off energy for plane wave basis sphere (Ha) eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree) gprimd(3,3)=dimensional reciprocal space primitive translations kg(3,mpw*mkmem)=reduced planewave coordinates. mband=maximum number of bands mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol mgfft=maximum size of 1D FFTs mgfftc=maximum size of 1D FFTs (coarse grid) mkmem =number of k points treated by this node. mpi_enreg=information about MPI parallelization mpw=maximum dimensioned size of npw. natom=number of atoms in cell. nattyp(ntypat)= # atoms of each type. nfft=(effective) number of FFT grid points (for this processor) (see NOTES at beginning of scfcv) ngfft(18)=contain all needed information about 3D FFT (see NOTES at beginning of scfcv) nkpt=number of k points. npwarr(nkpt)=number of planewaves in basis at this k point nsppol=1 for unpolarized, 2 for spin-polarized ntypat=number of types of atoms in unit cell. occ(mband*nkpt*nsppol) Occupation number for each band (often 2) for each k point. prtvol=control print volume and debugging output psps <type(pseudopotential_type)>=variables related to pseudopotentials rprimd(3,3)=dimensional primitive translations for real space (bohr) ucvol=unit cell volume (bohr**3) xred(3,natom)=reduced dimensionless atomic coordinates
OUTPUT
(only writing, printing)
SIDE EFFECTS
(only writing, printing)
NOTES
SOURCE
129 subroutine mlwfovlp(crystal, ebands, hdr, atindx1,cg,cprj,dtset,dtfil,eigen,gprimd,kg,& 130 & mband,mcg,mcprj,mgfftc,mkmem,mpi_enreg,mpw,natom,& 131 & nattyp,nfft,ngfft,nkpt,npwarr,nsppol,ntypat,occ,& 132 & pawang,pawrad,pawtab,prtvol,psps,rprimd,ucvol,xred) 133 134 !Arguments ------------------------------------ 135 !scalars 136 integer,intent(in) :: mband,mcg,mcprj,mgfftc,mkmem,mpw,natom,nfft,nkpt 137 integer,intent(in) :: nsppol,ntypat,prtvol 138 real(dp),intent(in) :: ucvol 139 type(crystal_t),intent(in) :: crystal 140 type(ebands_t),intent(in) :: ebands 141 type(hdr_type),intent(in) :: hdr 142 type(MPI_type),intent(in) :: mpi_enreg 143 type(dataset_type),intent(in) :: dtset 144 type(datafiles_type),intent(in) :: dtfil 145 type(pawang_type),intent(in) :: pawang 146 type(pseudopotential_type),intent(in) :: psps 147 !arrays 148 integer,intent(in) :: atindx1(natom) 149 integer :: kg(3,mpw*mkmem),nattyp(ntypat),ngfft(18),npwarr(nkpt) 150 real(dp),intent(in) :: cg(2,mcg) 151 real(dp),intent(in) :: eigen(mband*nkpt*nsppol),gprimd(3,3),rprimd(3,3) 152 real(dp),intent(in) :: occ(mband*nkpt*nsppol) 153 real(dp),intent(in) :: xred(3,natom) 154 type(pawcprj_type) :: cprj(natom,mcprj) 155 type(pawrad_type),intent(in) :: pawrad(psps%ntypat*psps%usepaw) 156 type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw) 157 158 !Local variables------------------------------- 159 !scalars 160 integer :: band_index,cplex,i,iatom,iband,iband1,iband2,icgtemp 161 integer :: ig,ii,ikg,ierr 162 integer :: ikpt,ikpt1,ikpt2,ilmn,intot,isppol,itypat 163 integer :: iun(nsppol),iun_plot,iwan,jband,jband1,jband2,jj,jj1,jj2,jj3 164 integer :: lmn_size,lproj,lwanniersetup,mwan,mgfft,n1 165 #if defined HAVE_WANNIER90 166 integer :: kk 167 #ifdef HAVE_NETCDF 168 integer :: ncid, ncerr, nrpts 169 character(len=fnlen) :: abiwan_fname 170 integer :: have_disentangled_spin(nsppol) 171 integer,allocatable :: irvec(:,:),ndegen(:) 172 #endif 173 #endif 174 integer :: n1tmp,n2,n2tmp,n3,n3tmp,n4,n5,n6,nband_k 175 integer :: nntot,npw_k,num_nnmax,spacing 176 integer :: tim_fourwf 177 integer :: master,max_num_bands,nprocs,spaceComm,rank 178 integer :: nwan(nsppol),nband_inc(nsppol),num_bands(nsppol) 179 real(dp) :: weight 180 #if defined HAVE_WANNIER90 181 real(dp) :: corrvdw 182 complex(dpc) :: caux,caux2,caux3 183 #endif 184 logical :: gamma_only,leig,lmmn,lwannierrun,spinors !,have_disentangled 185 character(len=fnlen) :: wfnname 186 character(len=1000) :: message 187 character(len=fnlen) :: seed_name(nsppol) 188 character(len=fnlen) :: fname,filew90_win(nsppol),filew90_wout(nsppol),filew90_amn(nsppol),filew90_ramn(nsppol) 189 character(len=fnlen) :: filew90_mmn(nsppol),filew90_eig(nsppol) 190 !arrays 191 integer :: g1temp(3),ngkpt(3) 192 integer,allocatable :: g1(:,:,:),gbound(:,:),icg(:,:) 193 integer,allocatable:: iwav(:,:,:),kg_k(:,:),ovikp(:,:) 194 integer,allocatable :: proj_l(:,:),proj_m(:,:),proj_radial(:,:) 195 integer,allocatable :: proj_s_loc(:) 196 real(dp) :: real_lattice(3,3) 197 real(dp) :: recip_lattice(3,3) 198 real(dp),allocatable :: cm1(:,:,:,:,:,:),cm2_paw(:,:,:),cwavef(:,:) 199 real(dp),allocatable :: denpot(:,:,:) 200 real(dp),allocatable :: eigenvalues_w(:,:,:),fofgout(:,:),fofr(:,:,:,:) 201 real(dp),allocatable :: proj_site(:,:,:),proj_x(:,:,:),proj_z(:,:,:),proj_zona(:,:) 202 real(dp),allocatable :: wann_centres(:,:,:),wann_spreads(:,:),xcart(:,:) 203 real(dp),allocatable :: proj_s_qaxis_loc(:,:) 204 complex(dpc),allocatable :: A_paw(:,:,:,:) 205 complex(dpc),allocatable :: M_matrix(:,:,:,:,:),U_matrix(:,:,:,:) 206 complex(dpc),allocatable :: U_matrix_opt(:,:,:,:) 207 complex(dpc),pointer :: A_matrix(:,:,:,:) 208 logical,allocatable :: band_in(:,:),lwindow(:,:,:) 209 character(len=3),allocatable :: atom_symbols(:) 210 logical,allocatable::just_augmentation(:,:) 211 #if defined HAVE_WANNIER90 212 real(dp) :: spreadw(3,nsppol) 213 real(dp),allocatable :: csix(:,:,:,:) 214 real(dpc),allocatable :: occ_arr(:,:,:),occ_wan(:,:,:) 215 real(dp),allocatable :: tdocc_wan(:,:) 216 #endif 217 218 !************************************************************************ 219 220 ABI_UNUSED((/crystal%natom, ebands%nkpt, hdr%nkpt/)) 221 222 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 223 !1) Initialize variables and allocations 224 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 225 ! 226 !Some initialization and checks 227 ! 228 lwanniersetup=1 ! 1 is mandatory ( 0 is for debug) 229 !to use lwanniersetup=0, one would need 230 !to define which bands to exclude. 231 lwannierrun=.true. ! .false. and .true. are possible 232 lmmn=.true. ! .false. and .true. are possible 233 leig=.true. ! .false. and .true. are possible 234 ! 235 gamma_only=.false. !not yet implemented 236 spinors=.false. !not yet implemented 237 ! 238 !mpi initialization 239 ! 240 spaceComm=MPI_enreg%comm_cell 241 nprocs=xmpi_comm_size(spaceComm) 242 rank=MPI_enreg%me_kpt 243 master=0 244 245 !write(std_out,'("master ",i3," rank ",i3," nprocs ",i3)') master,rank,nprocs 246 ! 247 !Generate seed names for wannier90 files, and file names 248 ! 249 call mlwfovlp_seedname(dtfil%fnameabo_w90,filew90_win,filew90_wout,filew90_amn,& 250 & filew90_ramn,filew90_mmn,filew90_eig,nsppol,seed_name) 251 ! 252 !Check the validity of input variables 253 !FIXME: this is not a check, and prints a warning even if the input is fine! 254 !must be changed to not print anything if kptopt 3 and istwfk 1 (the latter is easier to check) 255 ! 256 if(rank==master) then 257 write(message, '(a,a,a,a)' ) ch10,& 258 & ' mlwfovlp: you should give k-point in the full brillouin zone ',ch10,& 259 & ' with explicit k-points (or kptopt=3) and istwfk 1' 260 call wrtout(ab_out,message,'COLL') 261 call wrtout(std_out, message,'COLL') 262 end if 263 ! 264 if(MPI_enreg%paral_spinor==1) then 265 message = ' Parallelization over spinorial components not yet available !' 266 ABI_ERROR(message) 267 end if 268 269 if(nsppol==2) then 270 write(message, '(3a)' ) ch10,& 271 & ' mlwfovlp: Calculating matrices for both spin polarization ',ch10 272 call wrtout(ab_out,message,'COLL') 273 call wrtout(std_out, message,'COLL') 274 end if 275 ! 276 !get lattice parameters in wannier90 format 277 ! 278 do i=1, 3 279 real_lattice(:,i)=Bohr_Ang*rprimd(i,:) 280 recip_lattice(:,i)=two_pi*gprimd(i,:)/Bohr_Ang 281 end do 282 ! 283 if(psps%npsp/=psps%ntypat) then 284 ABI_ERROR("prb npsp") 285 end if 286 ! 287 !Allocations. 288 ! 289 num_nnmax=12 !limit fixed for compact structure in wannier_setup. 290 ABI_MALLOC(g1,(3,nkpt,num_nnmax)) 291 ABI_MALLOC(ovikp,(nkpt,num_nnmax)) 292 ABI_MALLOC(atom_symbols,(natom)) 293 ABI_MALLOC(xcart,(3,natom)) 294 ABI_MALLOC(band_in,(mband,nsppol)) 295 ABI_MALLOC(proj_site,(3,mband,nsppol)) 296 ABI_MALLOC(proj_l,(mband,nsppol)) 297 ABI_MALLOC(proj_m,(mband,nsppol)) 298 ABI_MALLOC(proj_radial,(mband,nsppol)) 299 ABI_MALLOC(proj_x,(3,mband,nsppol)) 300 ABI_MALLOC(proj_s_loc,(mband)) 301 ABI_MALLOC(proj_s_qaxis_loc,(3,mband)) 302 ABI_MALLOC(proj_z,(3,mband,nsppol)) 303 ABI_MALLOC(proj_zona,(mband,nsppol)) 304 ! 305 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 306 !2) Call to Wannier setup 307 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 308 ! 309 nullify(A_matrix) 310 311 ! 312 call mlwfovlp_setup(atom_symbols,band_in,dtset,filew90_win,gamma_only,& 313 & g1,lwanniersetup,mband,natom,nband_inc,nkpt,& 314 & nntot,num_bands,num_nnmax,nsppol,nwan,ovikp,& 315 & proj_l,proj_m,proj_radial,proj_site,proj_s_loc, proj_s_qaxis_loc, proj_x,proj_z,proj_zona,& 316 & real_lattice,recip_lattice,rprimd,seed_name,spinors,xcart,xred) 317 318 do isppol=1, nsppol 319 write(message, '(6a)' ) ch10,& 320 & ' mlwfovlp : mlwfovlp_setup done -',ch10,& 321 & '- see ',trim(filew90_wout(isppol)),' for details.' 322 call wrtout(ab_out,message,'COLL') 323 call wrtout(std_out, message,'COLL') 324 end do 325 326 ! 327 !some allocations after wannier90 setup 328 ! 329 max_num_bands=maxval(num_bands(:)) 330 mwan=maxval(nwan(:)) 331 ABI_MALLOC(eigenvalues_w,(max_num_bands,nkpt,nsppol)) 332 ABI_MALLOC(M_matrix,(max_num_bands,max_num_bands,nntot,nkpt,nsppol)) 333 ABI_MALLOC(A_matrix,(max_num_bands,mwan,nkpt,nsppol)) 334 ABI_MALLOC(iwav,(mband,nkpt,nsppol)) 335 336 ! 337 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 338 !3) Write Eigenvalues (file seed_name.eig) 339 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 340 ! 341 if(leig) then 342 ! 343 ! Assign file unit numbers 344 if(rank==master) then 345 iun(1)=444 346 if(nsppol==2) iun(2)=455 347 do isppol=1,nsppol 348 open(unit=iun(isppol),file=filew90_eig(isppol),form='formatted',status='unknown') 349 end do 350 end if !rank==master 351 ! Loop to write eigenvalues 352 band_index=0 353 do isppol=1,nsppol 354 do ikpt=1,nkpt 355 nband_k=dtset%nband(ikpt+(isppol-1)*nkpt) 356 jband=0 357 do iband=1,mband 358 if(band_in(iband,isppol)) then 359 jband=jband+1 360 ! Writing data 361 if(rank==master) write(iun(isppol), '(2i6,4x,f10.5)' ) jband,ikpt,Ha_eV*eigen(iband+band_index) 362 ! Finish writing, now save eigenvalues 363 eigenvalues_w(jband,ikpt,isppol)=Ha_eV*eigen(iband+band_index) 364 end if 365 end do !iband 366 band_index=band_index+nband_k 367 end do !ikpt 368 end do !nsppol 369 if(rank==master) then 370 do isppol=1,nsppol 371 close(iun(isppol)) 372 end do 373 write(message, '(a,a)' ) ch10,& 374 & ' mlwfovlp : eigenvalues written' 375 call wrtout(std_out, message,'COLL') 376 end if !master 377 end if !leig 378 !else if( leig . and. lwannierun) then 379 !read .eig file 380 381 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 382 !4) Calculate overlaps (file seed_name.mmn) 383 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 384 ! 385 !First calculate indices and shift 386 ! 387 !write(std_out,*) "Computes shift for cg" 388 write(message, '(a,a)' ) ch10,& 389 & ' mlwfovlp : compute shifts for g-points ' 390 call wrtout(std_out, message,'COLL') 391 !---------------------------------------------------------------------- 392 !Compute shifts for g points (icg,iwav) 393 !(here mband is not used, because shifts are internal variables of abinit) 394 !---------------------------------------------------------------------- 395 !write(std_out,*) mpw*dtset%nspinor*mband*mkmem*nsppol 396 ABI_MALLOC(icg,(nsppol,nkpt)) 397 icg=0 398 icgtemp=0 399 iwav(:,:,:)=0 400 do isppol=1,nsppol 401 do ikpt=1,nkpt 402 ! 403 ! MPI:cycle over k-points not treated by this node 404 ! 405 406 if (nprocs>1 ) then !sometimes we can have just one processor 407 if ( ABS(MPI_enreg%proc_distrb(ikpt,1,isppol)-rank) /=0) CYCLE 408 end if 409 410 ! write(std_out,*)'rank',rank,'ikpt',ikpt,'isppol',isppol 411 nband_k=dtset%nband(ikpt+(isppol-1)*nkpt) 412 ! write(std_out,*) ikpt+(isppol-1)*nkpt,nkpt 413 npw_k=npwarr(ikpt) 414 do iband=1,nband_k 415 if(iband.gt.mband) then 416 write(message,'(a,3i0)')" mband",iband,mband,nband_k 417 ABI_ERROR(message) 418 end if 419 iwav(iband,ikpt,isppol)= & 420 & (iband-1)*npw_k*dtset%nspinor+icgtemp 421 end do ! iband 422 icgtemp=icgtemp+ npw_k*dtset%nspinor*nband_k 423 ! icg(isppol,ikpt)=icgtemp 424 ! write(std_out,*) "icg", isppol,ikpt,icg(isppol,ikpt) 425 end do ! ikpt 426 end do ! isppol 427 !write(std_out,*) "shift for cg computed" 428 ABI_FREE(icg) 429 ! 430 !Shifts computed. 431 ! 432 if( lmmn) then 433 ! 434 ! In case of parallelization write out cg for all k-points 435 ! 436 if (nprocs > 1) then 437 ! 438 if(prtvol>0) then 439 write(message, '(3a)' ) ch10,& 440 & ' mlwfovlp : Creating temporary files with cg and cprj (PAW)',ch10 441 call wrtout(ab_out,message,'COLL') 442 call wrtout(std_out, message,'COLL') 443 end if 444 ! 445 do isppol=1,nsppol 446 do ikpt=1,nkpt 447 ! 448 ! MPI:cycle over k-points not treated by this node 449 ! 450 if (nprocs>1 ) then !sometimes we can have just one processor 451 if ( ABS(MPI_enreg%proc_distrb(ikpt,1,isppol)-rank) /=0) CYCLE 452 end if 453 454 ! write(std_out,*)'writing kpt ',ikpt,'isppol',isppol,' by node ', rank 455 write(wfnname,'(a,I5.5,".",I1)') trim(dtfil%fnametmp_cg),ikpt,isppol 456 iun_plot=1000+ikpt+ikpt*(isppol-1) 457 458 open (unit=iun_plot, file=wfnname,form='unformatted') 459 npw_k=npwarr(ikpt) 460 do iband=1,mband 461 do ig=1,npw_k*dtset%nspinor 462 write(iun_plot) (cg(i,ig+iwav(iband,ikpt,isppol)),i=1,2) 463 end do 464 end do 465 close(iun_plot) 466 end do !ikpt 467 end do !isppol 468 ! 469 ! In the PAW case we also need to write out cprj into files 470 ! 471 if(psps%usepaw==1) then 472 ! 473 ! big loop on atoms, kpts, bands and lmn 474 ! 475 ikpt2=0 476 do isppol=1,nsppol 477 do ikpt=1,nkpt 478 ! 479 ! MPI:cycle over k-points not treated by this node 480 ! 481 if (nprocs>1 ) then !sometimes we can have just one processor 482 if ( ABS(MPI_enreg%proc_distrb(ikpt,1,isppol)-MPI_enreg%me) /=0) CYCLE 483 end if 484 485 ikpt2=ikpt2+1 !sums just on the k-points treated by this node 486 ! 487 write(wfnname,'(a,I5.5,".",I1)') trim(dtfil%fnametmp_cprj),ikpt,isppol 488 iun_plot=1000+ikpt 489 open (unit=iun_plot, file=wfnname,form='unformatted') 490 ! 491 do iband=1,mband*dtset%nspinor 492 ig=iband+(ikpt2-1)*mband*dtset%nspinor +(isppol-1)*nkpt*mband*dtset%nspinor !index for cprj(:,ig) 493 ! 494 do iatom=1,natom 495 itypat=dtset%typat(iatom) 496 lmn_size=pawtab(itypat)%lmn_size 497 ! 498 do ilmn=1,lmn_size 499 write(iun_plot) (( cprj(iatom,ig)%cp(i,ilmn)),i=1,2) 500 end do !ilmn 501 end do !iatom 502 end do !iband 503 504 close(iun_plot) 505 end do !ikpt 506 end do !isppol 507 end if !usepaw==1 508 509 ! 510 ! 511 end if !MPI nprocs>1 512 ! 513 ! End of MPI preliminarities 514 ! Calculate PW contribution of overlaps 515 ! 516 ABI_MALLOC(cm1,(2,mband,mband,nntot,nkpt,nsppol)) 517 ! this loops over spin internally 518 call mlwfovlp_pw(cg,cm1,g1,iwav,kg,mband,& 519 & mkmem,mpi_enreg,mpw,nfft,ngfft,nkpt,nntot,& 520 & npwarr,dtset%nspinor,nsppol,ovikp,dtfil%fnametmp_cg) 521 write(message, '(a,a)' ) ch10,& 522 & ' mlwfovlp : PW part of overlap computed ' 523 call wrtout(std_out, message,'COLL') 524 ! 525 ! compute PAW Contribution and add it to PW contribution 526 ! 527 if(psps%usepaw==1) then 528 write(message, '(a,a)' ) ch10,& 529 & '** smatrix_pawinit : PAW part of overlap ' 530 call wrtout(std_out, message,'COLL') 531 ABI_MALLOC(cm2_paw,(2,mband,mband)) 532 do isppol=1,nsppol 533 do ikpt1=1,nkpt 534 ! 535 ! MPI:cycle over k-points not treated by this node 536 ! 537 if (nprocs>1 ) then !sometimes we can have just one processor 538 if ( ABS(MPI_enreg%proc_distrb(ikpt1,1,isppol)-rank) /=0) CYCLE 539 end if 540 541 write(message, '(a,i6,a,2i6)' ) & 542 & ' processor',rank,' computes PAW part for kpt and spin',ikpt1,isppol 543 call wrtout(std_out, message,'COLL') 544 545 do intot=1,nntot 546 ikpt2= ovikp(ikpt1,intot) 547 g1temp(:)=g1(:,ikpt1,intot) 548 call smatrix_pawinit(atindx1,cm2_paw,cprj,ikpt1,ikpt2,isppol,& 549 & g1temp,gprimd,dtset%kpt,mband,mband,mkmem,mpi_enreg,& 550 & natom,dtset%nband,nkpt,dtset%nspinor,nsppol,dtset%ntypat,pawang,pawrad,pawtab,rprimd,& 551 & dtfil%fnametmp_cprj,dtset%typat,xred) 552 ! cm1(:,:,:,intot,ikpt1,isppol)=four_pi*cm2_paw(:,:,:) 553 ! write(6,*) "ikpt1=",ikpt1 554 ! do iband=1,mband 555 ! write(6,*) "iband=",iband 556 ! write(6,*) "Wannier PW overlap",cm1(:,iband,iband,intot,ikpt1,isppol) 557 ! write(6,*) "Wannier PAW overlap",four_pi*cm2_paw(:,iband,iband) 558 ! write(6,*) "Wannier PW+PAW overlap",cm1(:,iband,iband,intot,ikpt1,isppol)+four_pi*cm2_paw(:,iband,iband) 559 ! enddo 560 cm1(:,:,:,intot,ikpt1,isppol)=cm1(:,:,:,intot,ikpt1,isppol)+four_pi*cm2_paw(:,:,:) 561 end do ! intot 562 end do ! ikpt1 563 end do ! isppol 564 ABI_FREE(cm2_paw) 565 write(message, '(a,a)' ) ch10,& 566 & ' mlwfovlp : PAW part of overlap computed ' 567 call wrtout(std_out, message,'COLL') 568 end if ! usepaw 569 ! 570 call xmpi_barrier(spaceComm) 571 call xmpi_sum(cm1,spaceComm,ierr) 572 ! 573 ! write overlap for separate calculation of wannier functions 574 ! 575 if(rank==master) then 576 do isppol=1,nsppol !we write separate output files for each isppol 577 iun(isppol)=220+isppol 578 open(unit=iun(isppol),file=filew90_mmn(isppol),form='formatted',status='unknown') 579 write(iun(isppol),*) "nnkp version 90" 580 write(iun(isppol),*) num_bands(isppol),nkpt,nntot 581 end do 582 end if ! rank==master 583 584 do isppol=1,nsppol 585 do ikpt1=1,nkpt 586 do intot=1,nntot 587 if( rank==master) write(iun(isppol),'(2i6,3x,3x,3i5)') ikpt1,ovikp(ikpt1,intot),(g1(jj,ikpt1,intot),jj=1,3) 588 jband2=0 589 do iband2=1,mband ! the first index is faster 590 if(band_in(iband2,isppol)) then 591 jband2=jband2+1 592 jband1=0 593 do iband1=1,mband 594 if(band_in(iband1,isppol)) then 595 jband1=jband1+1 596 if(rank==master) write(iun(isppol),*) & 597 & cm1(1,iband1,iband2,intot,ikpt1,isppol),cm1(2,iband1,iband2,intot,ikpt1,isppol) 598 M_matrix(jband1,jband2,intot,ikpt1,isppol)=& 599 & cmplx(cm1(1,iband1,iband2,intot,ikpt1,isppol),cm1(2,iband1,iband2,intot,ikpt1,isppol)) 600 ! write(2211,*) ikpt1,intot,iband1,iband2 601 ! write(2211,*) cm1(1,iband1,iband2,intot,ikpt1,isppol),cm1(2,iband1,iband2,intot,ikpt1,isppol) 602 end if ! band_in(iband1) 603 end do ! iband1 604 end if ! band_in(iband2) 605 end do ! iband2 606 end do !intot 607 end do !ikpt 608 if( rank==master ) then 609 close(iun(isppol)) 610 write(message, '(3a)' ) ' ',trim(filew90_mmn(isppol)),' written' 611 call wrtout(std_out, message,'COLL') 612 end if !rank==master 613 end do !isppol 614 ! 615 ABI_FREE(cm1) 616 ! 617 ! Write down part of the matrix to the output file 618 ! This is for the automatic tests 619 ! 620 if(rank==master) then 621 write(message, '(4a)' ) ch10,& 622 & ' Writing top of the overlap matrix: M_mn(ikb,ik)',ch10,& 623 & ' m=n=1:3, ikb=1, ik=1' 624 call wrtout(ab_out,message,'COLL') 625 call wrtout(std_out, message,'COLL') 626 ! 627 ! just write down the first 3 elements 628 ! 629 do isppol=1,nsppol 630 write(message, '( " " )') 631 if (nsppol>1 ) then 632 if (isppol==1) write(message,'(2a)')trim(message),' spin up:' 633 if (isppol==2) write(message,'(2a)')trim(message),' spin down:' 634 end if 635 do ii=1,3 636 if(ii>num_bands(isppol)) cycle 637 write(message,'(3a)') trim(message),ch10,'; ( ' 638 do jj=1,3 639 if(jj>num_bands(isppol))cycle 640 write(message, '(a,2f11.6,a)') trim(message),& 641 & M_matrix(ii,jj,1,1,isppol),' , ' 642 end do 643 write(message,'(2a)') trim(message),' ) ' 644 end do 645 call wrtout(ab_out,message,'COLL') 646 call wrtout(std_out, message,'COLL') 647 end do 648 ! 649 ! Now write down bottom of the matrix 650 ! 651 write(message, '(4a)' ) ch10,& 652 & ' Writing bottom of the overlap matrix: M_mn(ikb,ik)',ch10,& 653 & ' m=n=num_bands-2:num_bands, ikb=nntot, ik=nkpt' 654 call wrtout(ab_out,message,'COLL') 655 call wrtout(std_out, message,'COLL') 656 ! 657 do isppol=1,nsppol 658 write(message, '( " " )') 659 if (nsppol>1 ) then 660 if (isppol==1) write(message,'(2a)')trim(message),' spin up:' 661 if (isppol==2) write(message,'(2a)')trim(message),' spin down:' 662 end if 663 do ii=num_bands(isppol)-2,num_bands(isppol) 664 if(ii<1) cycle 665 write(message,'(3a)') trim(message),ch10,'; ( ' 666 do jj=num_bands(isppol)-2,num_bands(isppol) 667 if(jj<1)cycle 668 write(message, '(a,2f11.6,a)') trim(message),& 669 & M_matrix(ii,jj,nntot,nkpt,isppol),' , ' 670 end do !j 671 write(message,'(2a)') trim(message),' ) ' 672 end do !ii 673 call wrtout(ab_out,message,'COLL') 674 call wrtout(std_out, message,'COLL') 675 end do !isppol 676 end if !rank==master 677 ! 678 ! erase temporary files created for parallel runs 679 ! 680 if (nprocs > 1) then 681 ! 682 if(prtvol>0) then 683 write(message, '(3a)' ) ch10,& 684 & ' mlwfovlp : Removing temporary files with cg and cprj (PAW)',ch10 685 call wrtout(ab_out,message,'COLL') 686 call wrtout(std_out, message,'COLL') 687 end if 688 ! 689 ! Just master node will remove the files 690 ! 691 if(rank==master) then 692 do isppol=1,nsppol 693 do ikpt=1,nkpt 694 write(wfnname,'(a,I5.5,".",I1)') trim(dtfil%fnametmp_cg),ikpt,isppol 695 call delete_file(wfnname,ierr) 696 if(psps%usepaw==1) then 697 write(wfnname,'(a,I5.5,".",I1)') trim(dtfil%fnametmp_cprj),ikpt,isppol 698 call delete_file(wfnname,ierr) 699 end if 700 end do !ikpt 701 end do !isppol 702 end if 703 end if !MPI nprocs>1 704 ! 705 end if !lmmn 706 !if ( lmmn== .false. .and. lwannierun ) the 707 !read .mmn file 708 ! 709 !Deallocate arrays no longer used 710 ! 711 ABI_FREE(ovikp) 712 ABI_FREE(g1) 713 714 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 715 !5) Calculate initial projections 716 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 717 718 if(dtset%w90iniprj/=0 ) then 719 720 721 ! 722 ! Set value for lproj (type of projections to be computed) 723 ! In PAW, options 5 and 6 are not in use. 724 ! 5 means that there will be a contribution from inside the spheres and another from the PW part 725 ! 6 means that we take into account just the inside-spheres contribution 726 ! 2 means that PW part will be calculated 727 728 ! 729 lproj=dtset%w90iniprj 730 if(dtset%w90iniprj == 5 ) lproj=2 ! Necessary to calculate PW contribution 731 ! 732 ABI_MALLOC(just_augmentation,(mwan,nsppol)) 733 just_augmentation(:,:)=.false. 734 735 if( psps%usepaw==1 .and. (dtset%w90iniprj==2 .or. dtset%w90iniprj>4)) then 736 if (dtset%w90iniprj==6) just_augmentation(:,:)=.true. 737 if (dtset%w90iniprj==5) then 738 do isppol=1,nsppol 739 do iwan=1,nwan(isppol) 740 ! 741 ! Trick to skip the planewave contribution for some Wannier functions 742 ! (Not in production). 743 ! 744 if(proj_radial(iwan,isppol) > 4) then 745 just_augmentation(iwan,isppol)=.true. 746 proj_radial(iwan,isppol)=proj_radial(iwan,isppol)-3 747 write(message, '(2a,2i4)' ) & 748 & ' ','Skiping planewave contribution for iwan, ispin=',iwan,isppol 749 call wrtout(std_out, message,'COLL') 750 end if !proj_radial>4 751 end do !iwan 752 end do !isppol 753 end if !w90iniprj == 5 754 end if !paw 755 ! 756 ! Call mlwfovlp_proj (plane waves part of projections) 757 ! 758 if (dtset%w90iniprj/=6) then ! option 6 not yet in use 759 call mlwfovlp_proj(A_matrix,band_in,cg,cprj,dtset,gprimd,just_augmentation,kg,& 760 & lproj,max_num_bands,mband,mkmem,mpi_enreg,mpw,mwan,natom,& 761 & nattyp,nkpt,npwarr,& 762 & dtset%nspinor,nsppol,ntypat,num_bands,nwan,pawtab,proj_l,proj_m,& 763 & proj_radial,proj_site,proj_x,proj_z,proj_zona,psps,ucvol) 764 write(message, '(a,a,a,a)' ) ch10,& 765 & ' mlwfovlp: mlwfovlp_proj done -',ch10,& 766 & ' Projectors computed.' 767 call wrtout(std_out, message,'COLL') 768 end if !w90proj/=6 769 ! 770 ! Calculate inside-sphere part of projections (PAW) 771 ! 772 if (psps%usepaw ==1 .and. ( dtset%w90iniprj>4)) then 773 ABI_MALLOC(A_paw,(max_num_bands,mwan,nkpt,nsppol)) 774 call mlwfovlp_projpaw(A_paw,band_in,cprj,just_augmentation,max_num_bands,mband,mkmem,& 775 & mwan,natom,dtset%nband,nkpt,& 776 & dtset%nspinor,nsppol,dtset%ntypat,nwan,pawrad,pawtab,& 777 & proj_l,proj_m,proj_radial,proj_site,proj_x,proj_z,proj_zona,psps,& 778 & rprimd,dtset%typat,xred) 779 ! 780 write(message, '(a,a,a,a)' ) ch10,& 781 & ' mlwfovlp: mlwfovlp_proj_paw done -',ch10,& 782 & ' Inside-spheres part of projectors computed.' 783 call wrtout(std_out, message,'COLL') 784 ! 785 ! Add in-sphere contribution to A_matrix 786 ! 787 ! 788 ! w90iniprj==5. Plane waves + augmentation contributions 789 ! 790 if(dtset%w90iniprj==5) A_matrix(:,:,:,:)=A_matrix(:,:,:,:)+A_paw(:,:,:,:) 791 ! 792 ! w90iniprj==6. Just augmentation contribution 793 ! 794 if(dtset%w90iniprj==6) A_matrix(:,:,:,:)=A_paw(:,:,:,:) 795 ! 796 ! deallocations 797 ! 798 ABI_FREE(A_paw) 799 end if !usepaw==1 800 801 ABI_FREE(just_augmentation) 802 ! 803 call xmpi_barrier(spaceComm) 804 call xmpi_sum(A_matrix,spaceComm,ierr) 805 806 ! 807 ! write projections to a file 808 ! 809 if(rank==master) then 810 if(dtset%w90iniprj==1) then 811 do isppol=1,nsppol 812 iun(isppol)=219+isppol 813 open(unit=iun(isppol),file=trim(filew90_ramn(isppol)),form='formatted',status='unknown') 814 write(iun(isppol),*) 'Projections from Abinit : mband,nkpt,nwan. indices: iband1,iwan,ikpt' 815 write(iun(isppol),*) num_bands(isppol),nkpt,nwan(isppol) 816 end do 817 else 818 do isppol=1,nsppol 819 iun(isppol)=220+isppol 820 open(unit=iun(isppol),file=trim(filew90_amn(isppol)),form='formatted',status='unknown') 821 write(iun(isppol),*) 'Projections from Abinit : mband,nkpt,nwan. indices: iband1,iwan,ikpt' 822 write(iun(isppol),*) num_bands(isppol),nkpt,nwan(isppol) 823 end do 824 end if 825 ! 826 do isppol=1,nsppol 827 do ikpt=1,nkpt 828 do iwan=1,nwan(isppol) 829 jband=0 830 do iband=1,mband 831 if(band_in(iband,isppol)) then 832 jband=jband+1 833 write(iun(isppol),'(3i6,13x,3x,2f18.14)')jband,iwan,ikpt,A_matrix(jband,iwan,ikpt,isppol) 834 end if !band_in 835 end do !iband 836 end do !iwan 837 end do !ikpt 838 end do !isppol 839 ! 840 if(dtset%w90iniprj==1) then 841 do isppol=1,nsppol 842 close(iun(isppol)) 843 write(message, '(3a)' ) & 844 & ' ',trim(filew90_ramn(isppol)),' written' 845 call wrtout(std_out, message,'COLL') 846 end do 847 else 848 do isppol=1,nsppol 849 close(iun(isppol)) 850 write(message, '(3a)' ) & 851 & ' ',trim(filew90_amn(isppol)),' written' 852 call wrtout(std_out, message,'COLL') 853 end do 854 end if 855 end if !rank==master 856 857 ! 858 ! 859 ! Write down part of the matrix to the output file 860 ! This is for the automatic tests 861 ! 862 if(rank==master) then 863 write(message, '(4a)' ) ch10,& 864 & ' Writing top of the initial projections matrix: A_mn(ik)',ch10,& 865 & ' m=1:3, n=1:3, ik=1' 866 call wrtout(ab_out,message,'COLL') 867 call wrtout(std_out, message,'COLL') 868 ! 869 ! just write down the first 3 elements 870 ! 871 do isppol=1,nsppol 872 write(message, '( " " )') 873 if (nsppol>1 ) then 874 if (isppol==1) write(message,'(2a)')trim(message),' spin up:' 875 if (isppol==2) write(message,'(2a)')trim(message),' spin down:' 876 end if 877 do ii=1,3 878 if(ii>num_bands(isppol)) cycle 879 write(message,'(3a)') trim(message),ch10,'; ( ' 880 do jj=1,3 881 if(jj>nwan(isppol))cycle 882 write(message, '(a,2f11.6,a)') trim(message),& 883 & A_matrix(ii,jj,1,isppol),' , ' 884 end do 885 write(message,'(2a)') trim(message),' ) ' 886 end do 887 call wrtout(ab_out,message,'COLL') 888 call wrtout(std_out, message,'COLL') 889 end do 890 ! 891 ! Now write down bottom of the matrix 892 ! 893 write(message, '(4a)' ) ch10,& 894 & ' Writing bottom of the initial projections matrix: A_mn(ik)',ch10,& 895 & ' m=num_bands-2:num_bands, n=nwan-2:nwan, ik=nkpt' 896 call wrtout(ab_out,message,'COLL') 897 call wrtout(std_out, message,'COLL') 898 ! 899 do isppol=1,nsppol 900 write(message, '( " " )') 901 if (nsppol>1 ) then 902 if (isppol==1) write(message,'(2a)')trim(message),' spin up:' 903 if (isppol==2) write(message,'(2a)')trim(message),' spin down:' 904 end if 905 do ii=num_bands(isppol)-2,num_bands(isppol) 906 if(ii<1) cycle 907 write(message,'(3a)') trim(message),ch10,'; ( ' 908 do jj=nwan(isppol)-2,nwan(isppol) 909 if(jj<1)cycle 910 write(message, '(a,2f11.6,a)') trim(message),& 911 & A_matrix(ii,jj,nkpt,isppol),' , ' 912 end do 913 write(message,'(2a)') trim(message),' ) ' 914 end do 915 call wrtout(ab_out,message,'COLL') 916 call wrtout(std_out, message,'COLL') 917 end do !isppol 918 end if !rank==master 919 end if !dtset%w90iniprj/=0 920 ! 921 !Deallocations 922 ! 923 ABI_FREE(proj_site) 924 ABI_FREE(proj_l) 925 ABI_FREE(proj_m) 926 ABI_FREE(proj_radial) 927 ABI_FREE(proj_x) 928 ABI_FREE(proj_z) 929 ABI_FREE(proj_zona) 930 ABI_FREE(proj_s_loc) 931 ABI_FREE(proj_s_qaxis_loc) 932 933 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 934 !6) write files for wannier function plot 935 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 936 if( dtset%w90prtunk>0) then 937 if(psps%usepaw==1) then 938 write(message, '( a,a,a,a,a,a,a,a,a)')ch10,& 939 & " WARNING: The UNK matrices will not contain the correct wavefunctions ",ch10,& 940 & " since we are just writing the plane wave contribution.",ch10,& 941 & " The contribution from inside the spheres is missing. ",ch10,& 942 & " However, these files can be used for plotting purposes",ch10 943 call wrtout(std_out, message,'COLL') 944 end if 945 ! 946 spacing = dtset%w90prtunk 947 write(message, '( 8a,i3,2a)')ch10,& 948 & " UNK files will be written.",ch10,& 949 & " According to the chosen value of w90prtunk",ch10,& 950 & " the wavefunctions are to be written ",ch10, & 951 & " at every ", spacing," records.",ch10 952 call wrtout(std_out, message,'COLL') 953 ! 954 ABI_MALLOC(kg_k,(3,mpw)) 955 n1=ngfft(1) 956 n2=ngfft(2) 957 n3=ngfft(3) 958 n4=ngfft(4) 959 n5=ngfft(5) 960 n6=ngfft(6) 961 cplex=1 962 mgfft=mgfftc ! error 963 do isppol=1,nsppol 964 ikg=0 965 do ikpt=1,nkpt 966 ! 967 ! MPI:cycle over k-points not treated by this node 968 ! 969 if (nprocs>1 ) then !sometimes we can have just one processor 970 if ( ABS(MPI_enreg%proc_distrb(ikpt,1,isppol)-rank) /=0) CYCLE 971 end if 972 ! 973 npw_k=npwarr(ikpt) 974 kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 975 ABI_MALLOC(denpot,(cplex*n4,n5,n6)) 976 ABI_MALLOC(cwavef,(2,npw_k)) 977 ABI_MALLOC(fofr,(2,n4,n5,n6)) 978 ABI_MALLOC(gbound,(2*mgfft+8,2)) 979 ABI_MALLOC(fofgout,(2,npw_k)) 980 iun_plot=1000+ikpt+ikpt*(isppol-1) 981 write(wfnname,'("UNK",I5.5,".",I1)') ikpt, isppol 982 ! open (unit=iun_plot, file=wfnname,form='formatted') 983 open(unit=iun_plot, file=wfnname,form='unformatted') 984 ! optimizing grid for UNK files 985 n1tmp = n1/spacing 986 n2tmp = n2/spacing 987 n3tmp = n3/spacing 988 if( mod(n1,spacing) /= 0) then 989 n1tmp = n1tmp + 1 990 end if 991 if( mod(n2,spacing) /= 0) then 992 n2tmp = n2tmp + 1 993 end if 994 if( mod(n3,spacing) /= 0) then 995 n3tmp = n3tmp + 1 996 end if 997 ! write(iun_plot,*) n1tmp,n2tmp,n3tmp,ikpt,nband_inc 998 write(iun_plot) n1tmp,n2tmp,n3tmp,ikpt,nband_inc(isppol) 999 ! gbound=zero 1000 call sphereboundary(gbound,dtset%istwfk(ikpt),kg_k,mgfft,npw_k) 1001 write(std_out,*) " writes UNK file for ikpt, spin=",ikpt,isppol 1002 denpot(:,:,:)=zero 1003 weight = one 1004 do iband=1,mband 1005 if(band_in(iband,isppol)) then 1006 do ig=1,npw_k*dtset%nspinor 1007 cwavef(1,ig)=cg(1,ig+iwav(iband,ikpt,isppol)) 1008 cwavef(2,ig)=cg(2,ig+iwav(iband,ikpt,isppol)) 1009 end do 1010 tim_fourwf=0 1011 call fourwf(cplex,denpot,cwavef,fofgout,fofr,& 1012 & gbound,gbound,dtset%istwfk(ikpt),kg_k,kg_k,mgfft,& 1013 & mpi_enreg,1,ngfft,npw_k,npw_k,n4,n5,n6,0,& 1014 & tim_fourwf,weight,weight,gpu_option=dtset%gpu_option) 1015 ! do jj3=1,n3,spacing 1016 ! do jj2=1,n2,spacing 1017 ! do jj1=1,n1,spacing 1018 ! write(iun_plot,*) fofr(1,jj1,jj2,jj3),& 1019 ! & fofr(2,jj1,jj2,jj3) 1020 ! end do !jj1 1021 ! end do !jj2 1022 ! end do !jj3 1023 ! unformatted (must be one record) 1024 write(iun_plot) (((fofr(1,jj1,jj2,jj3),fofr(2,jj1,jj2,jj3),& 1025 & jj1=1,n1,spacing),jj2=1,n2,spacing),jj3=1,n3,spacing) 1026 end if !iband 1027 end do ! iband 1028 ABI_FREE(cwavef) 1029 ABI_FREE(fofr) 1030 ABI_FREE(gbound) 1031 ABI_FREE(denpot) 1032 ABI_FREE(fofgout) 1033 ikg=ikg+npw_k 1034 close(iun_plot) 1035 end do ! ikpt 1036 end do ! nsppol 1037 ABI_FREE(kg_k) 1038 ! 1039 write(message, '(4a)' )ch10, & 1040 & ' ','UNK files written',ch10 1041 call wrtout(std_out, message,'COLL') 1042 end if !dtset%w90prtunk 1043 ! 1044 ABI_FREE(iwav) 1045 1046 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1047 !7) Call to Wannier90 1048 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1049 if(lwannierrun) then 1050 if(lwanniersetup.ne.1) ABI_ERROR("lwanniersetup.ne.1") 1051 ABI_MALLOC(U_matrix,(mwan,mwan,nkpt,nsppol)) 1052 ABI_MALLOC(U_matrix_opt,(max_num_bands,mwan,nkpt,nsppol)) 1053 ABI_MALLOC(lwindow,(max_num_bands,nkpt,nsppol)) 1054 ABI_MALLOC(wann_centres,(3,mwan,nsppol)) 1055 ABI_MALLOC(wann_spreads,(mwan,nsppol)) 1056 ! Initialize 1057 U_matrix(:,:,:,:)=czero 1058 U_matrix_opt(:,:,:,:)=czero 1059 lwindow(:,:,:)=.false. 1060 wann_centres(:,:,:)=zero 1061 wann_spreads(:,:)=zero 1062 ! 1063 ! write(std_out,*) seed_name 1064 ! write(std_out,*) ngkpt 1065 ngkpt(1)=dtset%kptrlatt(1,1) 1066 ngkpt(2)=dtset%kptrlatt(2,2) ! ajouter test de verif que kptrlatt est bien diagonal 1067 ngkpt(3)=dtset%kptrlatt(3,3) 1068 ! write(std_out,*) nkpt 1069 ! write(std_out,*) rprimd*Bohr_Ang 1070 ! write(std_out,*) two_pi*gprimd/Bohr_Ang 1071 ! write(std_out,*) mband 1072 ! write(std_out,*) "nwan",nwan 1073 ! write(std_out,*) nntot 1074 ! write(std_out,*) natom 1075 ! write(std_out,*) atom_symbols 1076 ! write(std_out,*) xcart 1077 ! write(std_out,*) num_bands,num_bands,nntot,nkpt 1078 ! write(std_out,*) wann_spreads 1079 ! wann_spreads=2 1080 ! do i=1, nkpt 1081 ! do j=1, nntot 1082 ! write(std_out,*) i,j 1083 ! do k=1, num_bands 1084 ! do l=1, num_bands 1085 ! write(std_out,*) "m",M_matrix(l,k,j,i,1) 1086 ! enddo 1087 ! enddo 1088 ! enddo 1089 ! enddo 1090 1091 1092 1093 #if defined HAVE_WANNIER90 1094 do isppol=1,nsppol 1095 ! when nsppol>1, master runs isppol 1 and rank==1 runs isppol 2 1096 if(nprocs>1 .and. isppol==1.and.rank.ne.master) cycle 1097 if(nprocs>1 .and. isppol==2.and.rank.ne.1) cycle 1098 1099 write(message, '(8a)' ) ch10,& 1100 & '** mlwfovlp : call wannier90 library subroutine wannier_run ',ch10,& 1101 & ' Calculation is running ',ch10,& 1102 & '- see ',trim(filew90_wout(isppol)),' for details.' 1103 call wrtout(std_out, message,'COLL') 1104 ! 1105 call wannier_run(trim(seed_name(isppol)),ngkpt,nkpt,& !input 1106 & real_lattice,recip_lattice,dtset%kpt,num_bands(isppol),& !input 1107 & nwan(isppol),nntot,natom,atom_symbols,& !input 1108 & xcart*Bohr_Ang,gamma_only,M_matrix(:,:,:,:,isppol),A_matrix(:,:,:,isppol),eigenvalues_w(:,:,isppol),& !input 1109 & U_matrix(1:nwan(isppol),1:nwan(isppol),:,isppol),& !output 1110 & U_matrix_opt(1:num_bands(isppol),1:nwan(isppol),:,isppol),& !output 1111 & lwindow_loc=lwindow(1:num_bands(isppol),:,isppol),& !output 1112 & wann_centres_loc=wann_centres(:,1:nwan(isppol),isppol),& !output 1113 & wann_spreads_loc=wann_spreads(1:nwan(isppol),isppol),spread_loc=spreadw(:,isppol)) !output 1114 1115 ! ---------------------------------------------------------------------------------------------- 1116 1117 write(message, '(7a)' ) ch10,& 1118 & ' mlwfovlp : mlwfovlp_run completed -',ch10,& 1119 & '- see ',trim(filew90_wout(isppol)),' for details.',ch10 1120 call wrtout(ab_out,message,'COLL') 1121 call wrtout(std_out, message,'COLL') 1122 1123 end do !isppol 1124 1125 ! collect output of wannier90 from different processors 1126 call xmpi_barrier(spaceComm) 1127 1128 call xmpi_sum(U_matrix,spaceComm,ierr) 1129 call xmpi_sum(U_matrix_opt,spaceComm,ierr) 1130 call xmpi_lor(lwindow,spaceComm) 1131 call xmpi_sum(wann_centres,spaceComm,ierr) 1132 call xmpi_sum(wann_spreads,spaceComm,ierr) 1133 1134 ! Output ABIWAN.nc file 1135 #ifdef HAVE_NETCDF 1136 if (dtset%kptopt == 0) then 1137 ABI_WARNING("Output of ABIWAN.nc requires kptopt /= 0. ABIWAN.nc file won't be produced!") 1138 ! Need kptrlatt in wigner_seitz and client code need to know the k-grid. 1139 end if 1140 if (rank == master .and. dtset%kptopt /= 0) then 1141 abiwan_fname = strcat(dtfil%filnam_ds(4), "_ABIWAN.nc") 1142 call wrtout(std_out, sjoin("Saving wannier90 ouput results in:", abiwan_fname)) 1143 call wigner_seitz([zero, zero, zero], [2, 2, 2], dtset%kptrlatt, crystal%rmet, nrpts, irvec, ndegen) 1144 ! We know if disentanglement has been done by looking at the output values of lwindow 1145 ! Not elegant but it is the only way to avoid the parsing of the wannier input. 1146 ! In wannier_run lwindow is set to True if not disentanglement 1147 have_disentangled_spin = 0 1148 do isppol=1,nsppol 1149 !if nwan(isppol) < num_bands(isppol) 1150 if (.not. all(lwindow(:,:,isppol))) have_disentangled_spin(isppol) = 1 1151 end do 1152 1153 NCF_CHECK(nctk_open_create(ncid, abiwan_fname, xmpi_comm_self)) 1154 NCF_CHECK(hdr%ncwrite(ncid, fform_from_ext("ABIWAN"), nc_define=.True.)) 1155 NCF_CHECK(crystal%ncwrite(ncid)) 1156 NCF_CHECK(ebands_ncwrite(ebands, ncid)) 1157 1158 ncerr = nctk_def_dims(ncid, [ & 1159 nctkdim_t("mwan", mwan), & 1160 nctkdim_t("max_num_bands", max_num_bands), & 1161 nctkdim_t("nrpts", nrpts) & 1162 ], defmode=.True.) 1163 NCF_CHECK(ncerr) 1164 1165 ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "nntot"]) 1166 NCF_CHECK(ncerr) 1167 !ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: "fermi_energy", "smearing_width"]) 1168 !NCF_CHECK(ncerr) 1169 1170 ncerr = nctk_def_arrays(ncid, [ & 1171 nctkarr_t("nwan", "int", "number_of_spins"), & 1172 nctkarr_t("num_bands", "int", "number_of_spins"), & 1173 nctkarr_t("band_in_int", "int", "max_number_of_states, number_of_spins"), & 1174 nctkarr_t("lwindow_int", "int", "max_num_bands, number_of_kpoints, number_of_spins"), & 1175 !nctkarr_t("exclude_bands", "int", "max_number_of_states, number_of_spins"), & 1176 !nctkarr_t("eigenvalues_w", "int", "max_num_bands, number_of_kpoints, number_of_spins"), & 1177 nctkarr_t("spread", "dp", "three, number_of_spins"), & 1178 !nctkarr_t("A_matrix", "dp", "two, max_num_bands, mwan, number_of_kpoints, number_of_spins"), & 1179 nctkarr_t("irvec", "int", "three, nrpts"), & 1180 nctkarr_t("ndegen", "int", "nrpts"), & 1181 nctkarr_t("have_disentangled_spin", "int", "number_of_spins"), & 1182 nctkarr_t("U_matrix", "dp", "two, mwan, mwan, number_of_kpoints, number_of_spins"), & 1183 nctkarr_t("U_matrix_opt", "dp", "two, max_num_bands, mwan, number_of_kpoints, number_of_spins"), & 1184 nctkarr_t("wann_centres", "dp", "three, mwan, number_of_spins"), & 1185 nctkarr_t("wann_spreads", "dp", "mwan, number_of_spins") & 1186 ]) 1187 NCF_CHECK(ncerr) 1188 1189 ! Write data. 1190 NCF_CHECK(nctk_set_datamode(ncid)) 1191 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "nntot"), nntot)) 1192 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "nwan"), nwan)) 1193 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "num_bands"), num_bands)) 1194 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "band_in_int"), l2int(band_in))) 1195 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "lwindow_int"), l2int(lwindow))) 1196 !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "exclude_bands"), exclude_bands)) 1197 !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "eigenvalues_w"), eigenvalues_w)) 1198 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "spread"), spreadw)) 1199 !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "A_matrix"), c2r(A_matrix))) 1200 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "irvec"), irvec)) 1201 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ndegen"), ndegen)) 1202 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "have_disentangled_spin"), have_disentangled_spin)) 1203 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "U_matrix"), c2r(U_matrix))) 1204 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "U_matrix_opt"), c2r(U_matrix_opt))) 1205 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "wann_centres"), wann_centres)) 1206 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "wann_spreads"), wann_spreads)) 1207 NCF_CHECK(nf90_close(ncid)) 1208 1209 ABI_FREE(irvec) 1210 ABI_FREE(ndegen) 1211 end if 1212 #endif 1213 1214 ! CALL SILVESTRELLI'S APPROACH TO EVALUATE vdW INTERACTION ENERGY USING MLWF!! 1215 ! ---------------------------------------------------------------------------------------------- 1216 if (dtset%vdw_xc==10.or.dtset%vdw_xc==11.or.dtset%vdw_xc==12.or.dtset%vdw_xc==14.and.rank==master) then 1217 ! vdw_xc==10,11,12,14 starts the 1218 ! vdW interaction using MLWFs 1219 write(std_out,*) 'nwan(nsppol)=',ch10 1220 do ii=1,nsppol 1221 write(std_out,*) 'nsppol=',ii, 'nwan(nsppol)=',nwan(ii),ch10 1222 end do 1223 write(std_out,*) 'mwan=', mwan, ch10 1224 1225 ABI_MALLOC(occ_arr,(mband,nkpt,isppol)) 1226 ABI_MALLOC(occ_wan,(mwan,nkpt,nsppol)) 1227 ABI_MALLOC(tdocc_wan,(mwan,nsppol)) 1228 1229 occ_arr(:,:,:)=zero 1230 occ_wan(:,:,:)=zero 1231 tdocc_wan(:,:)=zero 1232 jj = 0 1233 do isppol=1,nsppol 1234 do ikpt=1,nkpt 1235 do iband=1,num_bands(isppol) 1236 jj = jj + 1 1237 occ_arr(iband,ikpt,isppol) = occ(jj) 1238 end do 1239 end do 1240 end do 1241 1242 do isppol=1,nsppol 1243 do ikpt=1,nkpt 1244 do iwan=1,nwan(isppol) 1245 caux=czero 1246 caux2=czero 1247 caux3=czero 1248 do iband=1,num_bands(isppol) !nband_inc(isppol) !nwan(isppol) 1249 do ii=1,nwan(isppol) 1250 caux=U_matrix(ii,iwan,ikpt,isppol)*U_matrix_opt(iband,ii,ikpt,isppol) 1251 ! DEBUG 1252 ! if(ISNAN(dble(caux))) then 1253 ! write(std_out,*) 'NaN: caux(ikpt,iwan,iband,ii):',ikpt,iwan,iband,ii,ch10 1254 ! end if 1255 ! END DEBUG 1256 do kk=1,nwan(isppol) 1257 caux2=conjg(U_matrix(kk,iwan,ikpt,isppol))*conjg(U_matrix_opt(iband,kk,ikpt,isppol)) 1258 caux3= caux3+caux*caux2*occ_arr(iband,ikpt,isppol) !take care here as exclude_bands case is not well 1259 ! DEBUG 1260 ! if(ISNAN(dble(caux2))) then 1261 ! write(std_out,*) 'NaN: caux2(ikpt,iwan,iband,kk):',ikpt,iwan,iband,kk,ch10 1262 ! end if 1263 ! if(ISNAN(dble(caux3))) then 1264 ! write(std_out,*) 'NaN: caux3(ikpt,iwan,iband,kk,jj):',ikpt,iwan,iband,kk,jj 1265 ! end if 1266 ! END DEBUG 1267 end do 1268 end do 1269 end do 1270 occ_wan(iwan,ikpt,isppol) = dble(caux3) 1271 ! DEBUG 1272 ! write(std_out,*) occ_wan(iwan,ikpt,isppol) 1273 ! END DEBUG 1274 ! end do 1275 end do 1276 end do 1277 end do 1278 1279 write(std_out,*) ch10,'MLWFs Occupation Matrix diagonal terms:',ch10 1280 1281 do jj=1,nsppol 1282 forall(iwan=1:nwan(jj)) tdocc_wan(iwan,jj) = sum(occ_wan(iwan,1:nkpt,jj)) / real(nkpt,dp) 1283 write(std_out,*) 'tdocc_wan(iwan),isppol:',ch10 1284 write(std_out,*) (tdocc_wan(iwan,jj),iwan=1,nwan(jj)),jj 1285 end do 1286 1287 ABI_MALLOC(csix,(mwan,mwan,nsppol,nsppol)) 1288 1289 call evdw_wannier(csix,corrvdw,mwan,natom,nsppol,nwan,tdocc_wan,dtset%vdw_nfrag,& 1290 & dtset%vdw_supercell,dtset%vdw_typfrag,dtset%vdw_xc,rprimd,wann_centres,wann_spreads,xcart) 1291 1292 ABI_FREE(csix) 1293 ABI_FREE(occ_arr) 1294 ABI_FREE(occ_wan) 1295 ABI_FREE(tdocc_wan) 1296 end if 1297 #else 1298 ABI_UNUSED(occ) 1299 #endif 1300 ! FIXME: looks like there is no automatic test which goes through here: g95 bot did not catch 1301 ! the missing deallocations 1302 ABI_FREE(wann_centres) 1303 ABI_FREE(wann_spreads) 1304 ABI_FREE(U_matrix) 1305 ABI_FREE(U_matrix_opt) 1306 ABI_FREE(lwindow) 1307 1308 end if !lwannierrun 1309 ! 1310 !deallocation 1311 ! 1312 ABI_FREE(band_in) 1313 ABI_FREE(atom_symbols) 1314 ABI_FREE(xcart) 1315 ABI_FREE(eigenvalues_w) 1316 ABI_FREE(M_matrix) 1317 ABI_FREE(A_matrix) 1318 1319 contains
m_mlwfovlp/mlwfovlp_proj [ Functions ]
[ Top ] [ m_mlwfovlp ] [ Functions ]
NAME
mlwfovlp_proj
FUNCTION
Routine which computes projection A_{mn}(k) for Wannier code (www.wannier.org f90 version).
INPUTS
cg(2,mpw*nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions
cprj(natom,nspinor*mband*mkmem*nsppol)= <p_lmn|Cnk> coefficients for each WF |Cnk>
and each |p_lmn> non-local projector
dtset <type(dataset_type)>=all input variables for this dataset
filew90_win = secondary input file for wannier90 (WAS NOT USED IN v6.7.1 - so has been temporarily removed)
kg(3,mpw*mkmem)=reduced planewave coordinates.
lproj= flag 0: no projections, 1: random projections,
2: projections on atomic orbitals
3: projections on projectors
mband=maximum number of bands
mkmem =number of k points treated by this node.
npwarr(nkpt)=number of planewaves in basis at this k point
mpi_enreg=information about MPI parallelization
mpw=maximum dimensioned size of npw.
natom=number of atoms in cell.
nattyp(ntypat)= # atoms of each type.
nkpt=number of k points.
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.
num_bands=number of bands actually used to construct the wannier function
nwan= number of wannier fonctions (read in wannier90.win).
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
OUTPUT
A_matrix(num_bands,nwan,nkpt,nsppol)= Matrix of projections needed by wannier_run ( also wannier90random.amn is written)
SIDE EFFECTS
(only writing, printing)
NOTES
SOURCE
2236 subroutine mlwfovlp_proj(A_matrix,band_in,cg,cprj,dtset,gprimd,just_augmentation,kg,& 2237 &lproj,max_num_bands,mband,mkmem,mpi_enreg,mpw,mwan,natom,nattyp,& 2238 &nkpt,npwarr,nspinor,& 2239 &nsppol,ntypat,num_bands,nwan,pawtab,proj_l,proj_m,proj_radial,& 2240 &proj_site,proj_x,proj_z,proj_zona,psps,ucvol) 2241 2242 !Arguments ------------------------------------ 2243 !scalars 2244 complex(dpc),parameter :: c1=(1._dp,0._dp) 2245 integer,intent(in) :: lproj,max_num_bands,mband,mkmem,mpw,mwan,natom,nkpt,nspinor,nsppol 2246 integer,intent(in) :: ntypat 2247 type(MPI_type),intent(in) :: mpi_enreg 2248 type(dataset_type),intent(in) :: dtset 2249 type(pseudopotential_type),intent(in) :: psps 2250 !arrays 2251 integer ::nattyp(ntypat) 2252 integer,intent(in) :: kg(3,mpw*mkmem),npwarr(nkpt),num_bands(nsppol),nwan(nsppol),proj_l(mband,nsppol) 2253 integer,intent(in) :: proj_m(mband,nsppol) 2254 integer,intent(inout)::proj_radial(mband,nsppol) 2255 real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol) 2256 real(dp),intent(in) :: gprimd(3,3),proj_site(3,mband,nsppol) 2257 real(dp),intent(in) :: proj_x(3,mband,nsppol),proj_z(3,mband,nsppol),proj_zona(mband,nsppol) 2258 complex(dpc),intent(out) :: A_matrix(max_num_bands,mwan,nkpt,nsppol) 2259 !character(len=fnlen),intent(in) :: filew90_win(nsppol) 2260 logical,intent(in) :: band_in(mband,nsppol) 2261 logical,intent(in)::just_augmentation(mwan,nsppol) 2262 type(pawcprj_type) :: cprj(natom,nspinor*mband*mkmem*nsppol) 2263 type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw) 2264 2265 !Local variables------------------------------- 2266 !scalars 2267 integer :: iatom,iatprjn,iband,iband1,iband2,ibg,icat,icg,icg_shift 2268 integer :: idum,idx,ikg,ikpt,ilmn,ipw,iproj 2269 integer :: ispinor,isppol,itypat,iwan,jband,jj1,libprjn 2270 integer :: lmn_size,natprjn,nband_k,nbprjn,npw_k 2271 integer :: sumtmp 2272 integer :: max_lmax,max_lmax2,mproj,nprocs,spaceComm,rank 2273 real(dp),parameter :: qtol=2.0d-8 2274 real(dp) :: arg,norm_error,norm_error_bar 2275 real(dp) :: ucvol,x1,x2,xnorm,xnormb,xx,yy,zz 2276 complex(dpc) :: amn_tmp(nspinor) 2277 complex(dpc) :: cstr_fact 2278 character(len=500) :: message 2279 !arrays 2280 integer :: kg_k(3,mpw),lmax(nsppol),lmax2(nsppol),nproj(nsppol) 2281 integer,allocatable :: lprjn(:),npprjn(:) 2282 real(dp) :: kpg(3),kpt(3) 2283 real(dp),allocatable :: amn(:,:,:,:,:),amn2(:,:,:,:,:,:,:) 2284 real(dp),allocatable :: gsum2(:),kpg2(:),radial(:) 2285 complex(dpc),allocatable :: gf(:,:),gft_lm(:) 2286 complex(dpc),allocatable :: ylmc_fac(:,:,:),ylmcp(:) 2287 2288 !no_abirules 2289 !Tables 3.1 & 3.2, User guide 2290 integer,save :: orb_l_defs(-5:3)=(/2,2,1,1,1,0,1,2,3/) 2291 ! integer,parameter :: mtransfo(0:3,7)=& 2292 !& reshape((/1,0,0,0,0,0,0,1,1,1,0,0,0,0,0,-2,-1,2,1,0,0,0,-1,1,2,-2,-3,3/),(/4,7/)) 2293 2294 !************************************************************************ 2295 2296 !mpi initialization 2297 spaceComm=MPI_enreg%comm_cell 2298 nprocs=xmpi_comm_size(spaceComm) 2299 rank=MPI_enreg%me_kpt 2300 2301 !Check input variables 2302 if ((lproj/=1).and.(lproj/=2).and.(lproj/=5)) then 2303 write(message, '(3a)' )& 2304 & ' Value of lproj no allowed ',ch10,& 2305 & ' Action : change lproj.' 2306 ABI_ERROR(message) 2307 end if 2308 2309 write(message, '(a,a)' )ch10,& 2310 & '** mlwfovlp_proj: compute A_matrix of initial guess for wannier functions' 2311 call wrtout(std_out,message,'COLL') 2312 2313 !Initialize to 0.d0 2314 A_matrix(:,:,:,:)=cmplx(0.d0,0.d0) 2315 2316 !End of preliminarities 2317 2318 !********************* Write Random projectors 2319 if(lproj==1) then 2320 idum=123456 2321 ! Compute random projections 2322 ABI_MALLOC(amn,(2,mband,mwan,nkpt,nsppol)) 2323 amn=zero 2324 do isppol=1,nsppol 2325 do ikpt=1,nkpt 2326 ! 2327 ! MPI: cycle over kpts not treated by this node 2328 ! 2329 if (ABS(MPI_enreg%proc_distrb(ikpt,1,isppol)-rank)/=0) CYCLE 2330 ! write(std_out,'("kpt loop2: ikpt",i3," rank ",i3)') ikpt,rank 2331 2332 ! 2333 do iband1=1,mband 2334 xnormb=0.d0 2335 do iband2=1,nwan(isppol) 2336 x1=uniformrandom(idum) 2337 x2=uniformrandom(idum) 2338 xnorm=sqrt(x1**2+x2**2) 2339 xnormb=xnormb+xnorm 2340 amn(1,iband1,iband2,ikpt,isppol)=x1 2341 amn(2,iband1,iband2,ikpt,isppol)=x2 2342 end do 2343 do iband2=1,nwan(isppol) 2344 amn(1,iband1,iband2,ikpt,isppol)=amn(1,iband1,iband2,ikpt,isppol)/xnormb 2345 amn(2,iband1,iband2,ikpt,isppol)=amn(2,iband1,iband2,ikpt,isppol)/xnormb 2346 end do !iband2 2347 end do !iband1 2348 end do !ikpt 2349 end do !isppol 2350 do isppol=1,nsppol 2351 do ikpt=1,nkpt 2352 ! 2353 ! MPI: cycle over kpts not treated by this node 2354 ! 2355 if (ABS(MPI_enreg%proc_distrb(ikpt,1,isppol)-rank)/=0) CYCLE 2356 ! 2357 do iband2=1,nwan(isppol) 2358 jband=0 2359 do iband1=1,mband 2360 if(band_in(iband1,isppol)) then 2361 jband=jband+1 2362 if(jband.gt.num_bands(isppol)) then 2363 write(message, '(3a)' )& 2364 & 'Value of jband is above num_bands ',ch10,& 2365 & 'Action: contact Abinit group' 2366 ABI_ERROR(message) 2367 end if 2368 A_matrix(jband,iband2,ikpt,isppol)=cmplx(amn(1,iband1,iband2,ikpt,isppol),amn(2,iband1,iband2,ikpt,isppol)) 2369 end if 2370 end do !iband1 2371 end do !iband2 2372 end do !ikpt 2373 end do !isppol 2374 ABI_FREE(amn) 2375 end if 2376 2377 !********************* Projection on atomic orbitals based on .win file 2378 if( lproj==2) then !based on .win file 2379 nproj(:)=nwan(:)/nspinor !if spinors, then the number of projections are 2380 mproj=maxval(nproj(:)) 2381 ! half the total of wannier functions 2382 ! 2383 ! obtain lmax and lmax2 2384 lmax(:)=0 2385 lmax2(:)=0 2386 ! 2387 do isppol=1,nsppol 2388 do iproj=1,nproj(isppol) 2389 lmax(isppol)=max(lmax(isppol),orb_l_defs(proj_l(iproj,isppol))) 2390 end do !iproj 2391 lmax2(isppol)=(lmax(isppol)+1)**2 2392 end do !isppol 2393 max_lmax=maxval(lmax(:)) 2394 max_lmax2=maxval(lmax2(:)) 2395 ! Allocate arrays 2396 ABI_MALLOC(ylmc_fac,(max_lmax2,mproj,nsppol)) 2397 ! 2398 ! get ylmfac, factor used for rotations and hybrid orbitals 2399 do isppol=1,nsppol 2400 call mlwfovlp_ylmfac(ylmc_fac(1:lmax2(isppol),1:nproj(isppol),isppol),lmax(isppol),lmax2(isppol),& 2401 & mband,nproj(isppol),proj_l(:,isppol),proj_m(:,isppol),proj_x(:,:,isppol),proj_z(:,:,isppol)) 2402 end do 2403 ! 2404 norm_error=zero 2405 norm_error_bar=zero 2406 icg=0 2407 ! 2408 do isppol=1,nsppol 2409 ! Allocate arrays 2410 ! this has to be done this way because the variable icg changes at the end of the 2411 ! cycle. We cannot just skip the hole cycle. 2412 ABI_MALLOC(gf,(mpw,nproj(isppol))) 2413 ABI_MALLOC(gft_lm,(lmax2(isppol))) 2414 ABI_MALLOC(gsum2,(nproj(isppol))) 2415 ABI_MALLOC(kpg2,(mpw)) 2416 ABI_MALLOC(radial,(lmax2(isppol))) 2417 ABI_MALLOC(ylmcp,(lmax2(isppol))) 2418 ! 2419 ikg=0 2420 do ikpt=1, nkpt 2421 ! 2422 ! MPI: cycle over kpts not treated by this node 2423 ! 2424 if (ABS(MPI_enreg%proc_distrb(ikpt,1,isppol)-rank)/=0) CYCLE 2425 ! 2426 write(message, '(a,i6,a,2i6)' ) & 2427 & ' processor',rank,' will compute k-point,spin=',ikpt,isppol 2428 call wrtout(std_out, message,'COLL') 2429 2430 ! Initialize variables 2431 npw_k=npwarr(ikpt) 2432 gsum2(:)=0.d0 2433 gf(:,:) = (0.d0,0.d0) 2434 kpt(:)=dtset%kpt(:,ikpt) 2435 kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 2436 2437 do ipw=1, npw_k 2438 kpg(1)= (kpt(1) + real(kg_k(1,ipw),dp)) !k+G 2439 kpg(2)= (kpt(2) + real(kg_k(2,ipw),dp)) 2440 kpg(3)= (kpt(3) + real(kg_k(3,ipw),dp)) 2441 ! 2442 ! Calculate modulus of k+G 2443 xx=gprimd(1,1)*kpg(1)+gprimd(1,2)*kpg(2)+gprimd(1,3)*kpg(3) 2444 yy=gprimd(2,1)*kpg(1)+gprimd(2,2)*kpg(2)+gprimd(2,3)*kpg(3) 2445 zz=gprimd(3,1)*kpg(1)+gprimd(3,2)*kpg(2)+gprimd(3,3)*kpg(3) 2446 kpg2(ipw)= two_pi*sqrt(xx**2+yy**2+zz**2) 2447 ! 2448 ! Complex Y_lm for k+G 2449 if(lmax(isppol)==0) then 2450 ylmcp(1)=c1/sqrt(four_pi) 2451 else 2452 call ylm_cmplx(lmax(isppol),ylmcp,xx,yy,zz) 2453 end if 2454 ! 2455 ! ! 2456 do iproj=1,nproj(isppol) 2457 ! 2458 ! In PAW, we can use proj_radial > 4 to indicate that we just 2459 ! want the in-sphere contribution 2460 ! 2461 if( psps%usepaw==1) then 2462 if( just_augmentation(iproj,isppol)) cycle 2463 end if 2464 ! 2465 ! obtain radial part 2466 call mlwfovlp_radial(proj_zona(iproj,isppol),lmax(isppol),lmax2(isppol)& 2467 & ,radial,proj_radial(iproj,isppol),kpg2(ipw)) 2468 ! 2469 ! scale complex representation of projector orbital with radial functions 2470 ! of appropriate l 2471 gft_lm(:)=radial(:)*ylmc_fac(1:lmax2(isppol),iproj,isppol) 2472 ! 2473 ! complex structure factor for projector orbital position 2474 arg = ( kpg(1)*proj_site(1,iproj,isppol) + & 2475 & kpg(2)*proj_site(2,iproj,isppol) + & 2476 & kpg(3)*proj_site(3,iproj,isppol) ) * 2*pi 2477 cstr_fact = cmplx(cos(arg), -sin(arg) ) 2478 ! 2479 ! obtain guiding functions 2480 gf(ipw,iproj)=cstr_fact*dot_product(ylmcp,gft_lm) 2481 ! 2482 gsum2(iproj)=gsum2(iproj)+real(gf(ipw,iproj))**2+aimag(gf(ipw,iproj))**2 2483 end do !iproj 2484 end do !ipw 2485 ! 2486 do iproj=1,nproj(isppol) 2487 ! 2488 ! In PAW, we can use proj_radial > 4 to indicate that we just 2489 ! want the in-sphere contribution 2490 ! 2491 if(psps%usepaw==1 ) then 2492 if (just_augmentation(iproj,isppol)) cycle 2493 end if 2494 ! 2495 gsum2(iproj)=16._dp*pi**2*gsum2(iproj)/ucvol 2496 gf(:,iproj)=gf(:,iproj)/sqrt(gsum2(iproj)) 2497 norm_error=max(abs(gsum2(iproj)-one),norm_error) 2498 norm_error_bar=norm_error_bar+(gsum2(iproj)-one)**2 2499 end do !iproj 2500 ! 2501 ! Guiding functions are computed. 2502 ! compute overlaps of gaussian projectors and wave functions 2503 do iproj=1,nproj(isppol) 2504 ! 2505 ! In PAW, we can use proj_radial > 4 to indicate that we just 2506 ! want the in-sphere contribution 2507 ! 2508 if(psps%usepaw==1 ) then 2509 if ( just_augmentation(iproj,isppol)) cycle 2510 end if 2511 ! 2512 jband=0 2513 do iband=1,mband 2514 if(band_in(iband,isppol)) then 2515 icg_shift=npw_k*nspinor*(iband-1)+icg 2516 jband=jband+1 2517 amn_tmp(:)=cmplx(0.d0,0.d0) 2518 do ispinor=1,nspinor 2519 do ipw=1,npw_k 2520 ! 2521 ! The case of spinors is tricky, we have nproj = nwan/2 2522 ! so we project to spin up and spin down separately, to have at 2523 ! the end an amn matrix with nwan projections. 2524 ! 2525 ! idx=ipw*nspinor - (nspinor-ispinor) 2526 idx=ipw+(ispinor-1)*npw_k 2527 amn_tmp(ispinor)=amn_tmp(ispinor)+gf(ipw,iproj)*cmplx(cg(1,idx+icg_shift),-cg(2,idx+icg_shift)) 2528 end do !ipw 2529 end do !ispinor 2530 do ispinor=1,nspinor 2531 iwan=(iproj*nspinor)- (nspinor-ispinor) 2532 A_matrix(jband,iwan,ikpt,isppol)=amn_tmp(ispinor) 2533 end do 2534 end if !band_in 2535 end do !iband 2536 end do !iproj 2537 icg=icg+npw_k*nspinor*mband 2538 ikg=ikg+npw_k 2539 end do !ikpt 2540 ! Deallocations 2541 ABI_FREE(gf) 2542 ABI_FREE(gft_lm) 2543 ABI_FREE(gsum2) 2544 ABI_FREE(kpg2) 2545 ABI_FREE(radial) 2546 ABI_FREE(ylmcp) 2547 end do !isppol 2548 ! 2549 ! if(isppol==1) then 2550 ! norm_error_bar=sqrt(norm_error_bar/real(nkpt*(nwan(1)),dp)) 2551 ! else 2552 ! norm_error_bar=sqrt(norm_error_bar/real(nkpt*(nwan(1)+nwan(2)),dp)) 2553 ! end if 2554 ! if(norm_error>0.05_dp) then 2555 ! write(message, '(6a,f6.3,a,f6.3,12a)' )ch10,& 2556 ! & ' mlwfovlp_proj : WARNING',ch10,& 2557 ! & ' normalization error for wannier projectors',ch10,& 2558 ! & ' is',norm_error_bar,' (average) and',norm_error,' (max).',ch10,& 2559 ! & ' this may indicate more cell-to-cell overlap of the radial functions',ch10,& 2560 ! & ' than you want.',ch10,& 2561 ! & ' Action : modify zona (inverse range of radial functions)',ch10,& 2562 ! ' under "begin projectors" in ',trim(filew90_win),' file',ch10 2563 ! call wrtout(std_out,message,'COLL') 2564 ! end if 2565 ! 2566 ! !Deallocate 2567 ! deallocate(ylmc_fac) 2568 ! 2569 ABI_FREE(ylmc_fac) 2570 end if !lproj==2 2571 2572 2573 !*************** computes projection from PROJECTORS ******************** 2574 if(lproj==3) then !! if LPROJPRJ 2575 ! ----- set values for projections --------------------- ! INPUT 2576 ! nbprjn:number of different l-values for projectors 2577 ! lprjn: value of l for each projectors par ordre croissant 2578 ! npprjn: number of projectors for each lprjn 2579 natprjn=1 ! atoms with wannier functions are first 2580 if(natprjn/=1) then ! in this case lprjn should depend on iatprjn 2581 ABI_ERROR("natprjn/=1") 2582 end if 2583 nbprjn=2 2584 ABI_MALLOC(lprjn,(nbprjn)) 2585 lprjn(1)=0 2586 lprjn(2)=1 2587 ABI_MALLOC(npprjn,(0:lprjn(nbprjn))) 2588 npprjn(0)=1 2589 npprjn(1)=1 2590 ! --- test coherence of nbprjn and nwan 2591 sumtmp=0 2592 do iatprjn=1,natprjn 2593 do libprjn=0,lprjn(nbprjn) 2594 sumtmp=sumtmp+(2*libprjn+1)*npprjn(libprjn) 2595 end do 2596 end do 2597 if(sumtmp/=nwan(1)) then 2598 write(std_out,*) "Number of Wannier orbitals is not equal to number of projections" 2599 write(std_out,*) "Action: check values of lprjn,npprjn % nwan" 2600 write(std_out,*) "nwan, sumtmp=",nwan,sumtmp 2601 ABI_ERROR("Aborting now") 2602 end if 2603 ! --- end test of coherence 2604 ABI_MALLOC(amn2,(2,natom,nsppol,nkpt,mband,nspinor,nwan(1))) 2605 if(psps%usepaw==1) then 2606 amn2=zero 2607 ibg=0 2608 do isppol=1,nsppol 2609 do ikpt=1,nkpt 2610 nband_k=dtset%nband(ikpt+(isppol-1)*nkpt) 2611 do iband=1,nband_k 2612 ! write(std_out,*)"amn2",iband,ibg,ikpt 2613 do ispinor=1,nspinor 2614 icat=1 2615 do itypat=1,dtset%ntypat 2616 lmn_size=pawtab(itypat)%lmn_size 2617 do iatom=icat,icat+nattyp(itypat)-1 2618 jj1=0 2619 do ilmn=1,lmn_size 2620 if(iatom.le.natprjn) then 2621 ! do iwan=1,nwan 2622 do libprjn=0,lprjn(nbprjn) 2623 ! if (psps%indlmn(1,ilmn,itypat)==proj_l(iwan)) then 2624 ! if (psps%indlmn(2,ilmn,itypat)==mtransfo(proj_l(iwan),proj_m(iwan))) then 2625 if (psps%indlmn(1,ilmn,itypat)==libprjn) then 2626 if (psps%indlmn(3,ilmn,itypat)<=npprjn(libprjn)) then 2627 if(band_in(iband,isppol)) then 2628 jj1=jj1+1 2629 if(jj1>nwan(isppol)) then 2630 write(std_out,*) "number of wannier orbitals is lower than lmn_size" 2631 write(std_out,*) jj1,nwan(isppol) 2632 ABI_ERROR("Aborting now") 2633 end if 2634 amn2(1,iatom,isppol,ikpt,iband,ispinor,jj1)=cprj(iatom,iband+ibg)%cp(1,ilmn) 2635 amn2(2,iatom,isppol,ikpt,iband,ispinor,jj1)=cprj(iatom,iband+ibg)%cp(2,ilmn) 2636 end if 2637 end if 2638 end if 2639 end do ! libprjn 2640 ! endif 2641 ! endif 2642 ! enddo ! iwan 2643 end if ! natprjn 2644 end do !ilmn 2645 end do ! iatom 2646 icat=icat+nattyp(itypat) 2647 end do ! itypat 2648 end do ! ispinor 2649 end do !iband 2650 ibg=ibg+nband_k*nspinor 2651 ! write(std_out,*)'amn2b',iband,ibg,ikpt 2652 end do !ikpt 2653 end do ! isppol 2654 2655 ! ----------------------- Save Amn -------------------- 2656 do isppol=1,nsppol 2657 do ikpt=1,nkpt 2658 do iband2=1,nwan(isppol) 2659 jband=0 2660 do iband1=1,mband 2661 if(band_in(iband1,isppol)) then 2662 jband=jband+1 2663 A_matrix(jband,iband2,ikpt,isppol)=& 2664 & cmplx(amn2(1,1,1,ikpt,iband1,1,iband2),amn2(2,1,1,ikpt,iband1,1,iband2)) 2665 end if 2666 end do 2667 end do 2668 end do 2669 end do 2670 end if !usepaw 2671 ABI_FREE(amn2) 2672 ABI_FREE(npprjn) 2673 ABI_FREE(lprjn) 2674 2675 end if ! lproj==3 2676 2677 end subroutine mlwfovlp_proj
m_mlwfovlp/mlwfovlp_projpaw [ Functions ]
[ Top ] [ m_mlwfovlp ] [ 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
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.
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
SOURCE
2730 subroutine mlwfovlp_projpaw(A_paw,band_in,cprj,just_augmentation,max_num_bands,mband,mkmem,& 2731 &mwan,natom,nband,nkpt,& 2732 &nspinor,nsppol,ntypat,nwan,pawrad,pawtab,& 2733 &proj_l,proj_m,proj_radial,proj_site,proj_x,proj_z,proj_zona,psps,& 2734 &rprimd,typat,xred) 2735 2736 !Arguments ------------------------------------ 2737 integer,intent(in) :: max_num_bands,mband,mkmem,mwan,natom,nkpt 2738 integer,intent(in) :: nspinor,nsppol,ntypat 2739 !arrays 2740 integer,intent(in) :: nband(nsppol*nkpt),nwan(nsppol) 2741 integer,intent(in) :: proj_l(mband,nsppol),proj_m(mband,nsppol),proj_radial(mband,nsppol) 2742 integer,intent(in) :: typat(natom) 2743 real(dp),intent(in):: proj_site(3,mband,nsppol) 2744 real(dp),intent(in) :: proj_x(3,mband,nsppol),proj_z(3,mband,nsppol),proj_zona(mband,nsppol) 2745 real(dp),intent(in) :: rprimd(3,3),xred(3,natom) 2746 complex(dpc),intent(out) :: A_paw(max_num_bands,mwan,nkpt,nsppol) 2747 logical,intent(in) :: band_in(mband,nsppol) 2748 logical,intent(in)::just_augmentation(mwan,nsppol) 2749 type(pawcprj_type) :: cprj(natom,nspinor*mband*mkmem*nsppol) 2750 type(pawrad_type),intent(in) :: pawrad(ntypat) 2751 type(pawtab_type),intent(in) :: pawtab(ntypat) 2752 type(pseudopotential_type),intent(in) :: psps 2753 2754 !Local variables------------------------------- 2755 !local variables 2756 integer :: basis_size,iatom,iband,ii 2757 integer :: ikpt,ir,isppol,itypat,iwan,jband 2758 integer :: ll,lm,ln,mm,ilmn 2759 integer :: lmn_size,max_lmax2, mesh_size,nn 2760 integer :: lmax(nsppol),lmax2(nsppol) 2761 real(dp):: aa,int_rad2,prod_real,prod_imag 2762 real(dp),parameter :: dx=0.015d0,rmax=10.d0,xmin=0.d0 2763 real(dp):: sum,wan_lm_fac,x 2764 complex(dpc)::prod 2765 character(len=500) :: message 2766 !arrays 2767 integer :: index(mband,nkpt,nsppol) 2768 real(dp):: dist,norm(mwan,nsppol) 2769 real(dp):: proj_cart(3,mwan,nsppol),proj_site_unit(3,mwan,nsppol) 2770 real(dp):: xcart_unit(3,natom),xred_unit(3,natom) 2771 real(dp),allocatable ::aux(:),ff(:),r(:),int_rad(:),rad_int(:) 2772 real(dp),allocatable::ylmr_fac(:,:,:) 2773 2774 !no_abirules 2775 !Tables 3.1 & 3.2, User guide 2776 integer,save :: orb_l_defs(-5:3)=(/2,2,1,1,1,0,1,2,3/) 2777 !real(dp),allocatable :: ylm(:,:) 2778 2779 2780 ! ************************************************************************* 2781 2782 write(message, '(a,a)' )ch10,'** mlwfovlp_proj: compute in-sphere part of A_matrix' 2783 call wrtout(std_out,message,'COLL') 2784 2785 !Check input variables 2786 do isppol=1,nsppol 2787 do iwan=1,nwan(nsppol) 2788 if(proj_radial(iwan,isppol)<1 .or. proj_radial(iwan,isppol)>4)then 2789 write(message,'(a,a,a,i6)')& 2790 & ' proj_radial should be between 1 and 4,',ch10,& 2791 & ' however, proj_radial=',proj_radial(iwan,isppol) 2792 ABI_BUG(message) 2793 end if 2794 end do 2795 end do 2796 2797 !Initialize 2798 A_paw(:,:,:,:)=cmplx(0.d0,0.d0) 2799 2800 !Get index for cprj 2801 ii=0 2802 do isppol=1,nsppol 2803 do ikpt=1,nkpt 2804 do iband=1,nband(ikpt) 2805 ii=ii+1 2806 index(iband,ikpt,isppol)=ii 2807 end do 2808 end do 2809 end do 2810 2811 !obtain lmax and lmax2 2812 lmax(:)=0 2813 lmax2(:)=0 2814 do isppol=1,nsppol 2815 do iwan=1,nwan(isppol) 2816 lmax(isppol)=max(lmax(isppol),orb_l_defs(proj_l(iwan,isppol))) 2817 end do !iwan 2818 lmax2(isppol)=(lmax(isppol)+1)**2 2819 end do 2820 max_lmax2=maxval(lmax2(:)) 2821 ! 2822 !get ylmfac, factor used for rotations and hybrid orbitals 2823 ! 2824 ABI_MALLOC(ylmr_fac,(max_lmax2,mwan,nsppol)) 2825 2826 2827 do isppol=1,nsppol 2828 call mlwfovlp_ylmfar(ylmr_fac(1:lmax2(isppol),1:nwan(isppol),isppol),& 2829 & lmax(isppol),lmax2(isppol),mband,nwan(isppol),proj_l(:,isppol),proj_m(:,isppol),& 2830 & proj_x(:,:,isppol),proj_z(:,:,isppol)) 2831 ! 2832 ! Shift projection centers and atom centers to the primitive cell 2833 ! This will be useful after, when we check if the Wannier function 2834 ! lies on one specific atom 2835 ! 2836 proj_site_unit(:,:,:)=0.d0 2837 do iwan=1,nwan(isppol) 2838 do ii=1,3 2839 proj_site_unit(ii,iwan,isppol)=ABS(proj_site(ii,iwan,isppol)-AINT(proj_site(ii,iwan,isppol)) ) 2840 end do 2841 end do 2842 do iatom=1,natom 2843 do ii=1,3 2844 xred_unit(ii,iatom)=ABS(xred(ii,iatom)-AINT(xred(ii,iatom)) ) 2845 end do 2846 end do 2847 call xred2xcart(natom,rprimd,xcart_unit,xred_unit) 2848 call xred2xcart(mwan,rprimd,proj_cart(:,:,isppol),proj_site_unit(:,:,isppol)) 2849 ! 2850 ! Normalize the Wannier functions 2851 ! 2852 ! Radial part 2853 mesh_size= nint((rmax - xmin ) / dx + 1) 2854 ABI_MALLOC( ff,(mesh_size)) 2855 ABI_MALLOC(r,(mesh_size)) 2856 ABI_MALLOC(rad_int,(mesh_size)) 2857 ABI_MALLOC(aux,(mesh_size)) 2858 do ir=1, mesh_size 2859 x=xmin+DBLE(ir-1)*dx 2860 r(ir)=x 2861 end do !ir 2862 do iwan=1,nwan(isppol) 2863 ! write(std_out,*)'iwan',iwan 2864 ! radial functions shown in table 3.3 of wannier90 manual 2865 if(proj_radial(iwan,isppol)==1) ff(:) = 2.d0 * proj_zona(iwan,isppol)**(1.5d0) * exp(-proj_zona(iwan,isppol)*r(:)) 2866 if(proj_radial(iwan,isppol)==2) ff(:) = 1.d0/(2.d0*sqrt(2.d0))*proj_zona(iwan,isppol)**(1.5d0) *& 2867 & (2.d0 - proj_zona(iwan,isppol)*r(:))*exp(-proj_zona(iwan,isppol)*r(:)/2.d0) 2868 if(proj_radial(iwan,isppol)==3) ff(:) = sqrt(4.d0/27.d0)*proj_zona(iwan,isppol)**(1.5d0)& 2869 & * (1.d0 - 2.d0*proj_zona(iwan,isppol)*r(:)/3.d0 + 2.d0*proj_zona(iwan,isppol)**2*r(:)**2/27.d0)& 2870 & * exp(-proj_zona(iwan,isppol) * r(:)/3.d0) 2871 2872 if(proj_radial(iwan,isppol)/=4) then 2873 aux(:)=ff(:)**2*r(:)**2 2874 call simpson_int(mesh_size,dx,aux,rad_int) 2875 sum=0.d0 2876 do ir=1,mesh_size 2877 sum=sum+rad_int(ir) 2878 end do 2879 int_rad2=sum/real(mesh_size,dp) 2880 ! 2881 ! do ir=1,mesh_size 2882 ! if(iwan==1) write(400,*)r(ir),aux(ir),rad_int(ir) 2883 ! end do 2884 else 2885 ! 2886 ! ==4: gaussian function 2887 ! f(x)=\exp(-1/4(x/aa)**2) 2888 ! \int f(x)f(x) dx = \int \exp(-1/2(x/aa)**2) = aa*sqrt(2pi) 2889 ! 2890 int_rad2=sqrt(2.d0*pi)*proj_zona(iwan,isppol) 2891 end if 2892 2893 ! 2894 ! Now angular part 2895 ! 2896 prod_real=0.d0 2897 do lm=1,lmax2(isppol) 2898 wan_lm_fac=ylmr_fac(lm,iwan,isppol) 2899 ! write(std_out,*)'wan_lm_fac',wan_lm_fac 2900 ! write(std_out,*)'int_rad2',int_rad2 2901 prod_real= prod_real + wan_lm_fac**2 * int_rad2 2902 end do 2903 norm(iwan,isppol)=sqrt(prod_real) 2904 end do !iwan 2905 ABI_FREE(ff) 2906 ABI_FREE(r) 2907 ABI_FREE(rad_int) 2908 ABI_FREE(aux) 2909 ! 2910 ! Now that we found our guiding functions 2911 ! We proceed with the internal product of 2912 ! our guiding functions and the wave function 2913 ! Amn=<G_m|\Psi_n> inside the sphere. 2914 ! The term <G_m|\Psi_n> inside the sphere is: 2915 ! = \sum_i <G_n | \phi_i - \tphi_i> <p_im|\Psi_m> 2916 ! 2917 ! 2918 ! G_n \phi and \tphi can be decomposed in 2919 ! a radial function times an angular function. 2920 ! 2921 ! 2922 ! Big loop on iwan and iatom 2923 ! 2924 do iwan=1,nwan(isppol) 2925 do iatom=1,natom 2926 ! 2927 ! check if center of wannier function coincides 2928 ! with the center of the atom 2929 ! 2930 dist=((proj_cart(1,iwan,isppol)-xcart_unit(1,iatom))**2 + & 2931 & (proj_cart(2,iwan,isppol)-xcart_unit(2,iatom))**2 + & 2932 & (proj_cart(3,iwan,isppol)-xcart_unit(3,iatom))**2)**0.5 2933 ! 2934 ! if the distance between the centers is major than 0.1 angstroms skip 2935 ! 2936 if( dist > 0.188972613) cycle 2937 ! 2938 write(message, '(2a,i4,a,i4,2a)')ch10, ' Wannier function center',iwan,' is on top of atom',& 2939 & iatom,ch10,' Calculating in-sphere contribution' 2940 call wrtout(ab_out,message,'COLL') 2941 call wrtout(std_out, message,'COLL') 2942 ! 2943 ! Get useful quantities 2944 ! 2945 itypat=typat(iatom) 2946 lmn_size=pawtab(itypat)%lmn_size 2947 basis_size=pawtab(itypat)%basis_size 2948 mesh_size=pawtab(itypat)%mesh_size 2949 ABI_MALLOC(int_rad,(basis_size)) 2950 ABI_MALLOC(ff,(mesh_size)) 2951 ABI_MALLOC(aux,(mesh_size)) 2952 2953 ! 2954 ! Integrate first the radial part 2955 ! and save it into an array 2956 ! 2957 ! 2958 ! radial functions shown in table 3.3 of wannier90 manual 2959 ! 2960 if(proj_radial(iwan,isppol)==1) aux(1:mesh_size) = 2.d0 * proj_zona(iwan,isppol)**(1.5d0) *& 2961 & exp(-proj_zona(iwan,isppol)*pawrad(itypat)%rad(1:mesh_size)) 2962 if(proj_radial(iwan,isppol)==2) aux(1:mesh_size) = 1.d0/(2.d0*sqrt(2.d0))*proj_zona(iwan,isppol)**(1.5d0) *& 2963 & (2.d0 - proj_zona(iwan,isppol)*pawrad(itypat)%rad(1:mesh_size)) & 2964 & * exp(-proj_zona(iwan,isppol)*pawrad(itypat)%rad(1:mesh_size)/2.d0) 2965 if(proj_radial(iwan,isppol)==3) aux(1:mesh_size) = sqrt(4.d0/27.d0)*proj_zona(iwan,isppol)**(1.5d0)& 2966 & * (1.d0 - 2.d0*proj_zona(iwan,isppol)*pawrad(itypat)%rad(1:mesh_size)/3.d0 & 2967 & + 2.d0*proj_zona(iwan,isppol)**2 *pawrad(itypat)%rad(1:mesh_size)**2/27.d0)& 2968 & * exp(-proj_zona(iwan,isppol) * pawrad(itypat)%rad(1:mesh_size)/3.d0) 2969 ! 2970 ! ==4: gaussian function 2971 ! f(x)=\exp(-1/4(x/aa)**2) 2972 ! 2973 if(proj_radial(iwan,isppol)==4) then 2974 aa=1.d0/proj_zona(iwan,isppol) 2975 aux(1:mesh_size)= exp(-0.25d0*(pawrad(itypat)%rad(1:mesh_size)*aa)**2) 2976 end if 2977 ! 2978 ! Normalize aux 2979 aux(:)=aux(:)/norm(iwan,isppol) 2980 ! 2981 do ln=1,basis_size 2982 if(just_augmentation(iwan,isppol)) then 2983 ! 2984 ! just augmentation region contribution 2985 ! In this case there is no need to use \tphi 2986 ! ff= \int R_wan(r) (R_phi(ln;r)/r ) r^2 dr 2987 ! 2988 ff(1:mesh_size)= aux(1:mesh_size) * pawtab(itypat)%phi(1:mesh_size,ln) & 2989 & * pawrad(itypat)%rad(1:mesh_size) 2990 else 2991 ! Inside sphere contribution = \phi - \tphi 2992 ! ff= \int R_wan(r) (R_phi(ln;r)/r - R_tphi(ln;r)/r) r^2 dr 2993 ff(1:mesh_size)= aux(1:mesh_size) * (pawtab(itypat)%phi(1:mesh_size,ln)-pawtab(itypat)%tphi(1:mesh_size,ln)) & 2994 & * pawrad(itypat)%rad(1:mesh_size) 2995 end if 2996 ! 2997 ! Integration with simpson routine 2998 ! 2999 call simp_gen(int_rad(ln),ff,pawrad(itypat)) 3000 ! do ii=1,mesh_size 3001 ! unit_ln=400+ln 3002 ! if( iwan==1 ) write(unit_ln,*)pawrad(itypat)%rad(ii),ff(ii),int_rad(ln) 3003 ! end do 3004 end do !ln 3005 ABI_FREE(ff) 3006 ABI_FREE(aux) 3007 ! 3008 ! Now integrate the angular part 3009 ! Cycle on i indices 3010 ! 3011 ! prod_real=0.d0 3012 do ilmn=1, lmn_size 3013 ll=Psps%indlmn(1,ilmn,itypat) 3014 mm=Psps%indlmn(2,ilmn,itypat) 3015 nn=Psps%indlmn(3,ilmn,itypat) 3016 lm=Psps%indlmn(4,ilmn,itypat) 3017 ln=Psps%indlmn(5,ilmn,itypat) 3018 ! write(std_out,*)'ll ',ll,' mm ',mm,'nn',nn,"lm",lm,"ln",ln 3019 ! 3020 ! Get wannier factor for that lm component 3021 if(lm <=lmax2(isppol)) then 3022 wan_lm_fac=ylmr_fac(lm,iwan,isppol) 3023 ! Make delta product 3024 ! Here we integrate the angular part 3025 ! Since the integral of the product of two spherical harmonics 3026 ! is a delta function 3027 if( abs(wan_lm_fac) > 0.0d0) then 3028 ! write(std_out,*) 'll',ll,'mm',mm,'lm',lm,'ln',ln,'factor',wan_lm_fac !lm index for wannier function 3029 ! 3030 ! Calculate Amn_paw, now that the radial and angular integrations are done 3031 ! 3032 prod=cmplx(0.d0,0.d0) 3033 do ikpt=1,nkpt 3034 jband=0 3035 do iband=1,nband(ikpt) 3036 if(band_in(iband,isppol)) then 3037 jband=jband+1 3038 3039 prod_real= cprj(iatom,index(iband,ikpt,isppol))%cp(1,ilmn) * int_rad(ln) * wan_lm_fac 3040 prod_imag= cprj(iatom,index(iband,ikpt,isppol))%cp(2,ilmn) * int_rad(ln) * wan_lm_fac 3041 prod=cmplx(prod_real,prod_imag) 3042 3043 A_paw(jband,iwan,ikpt,isppol)=A_paw(jband,iwan,ikpt,isppol)+prod 3044 end if !band_in 3045 end do !iband 3046 end do !ikpt 3047 ! 3048 end if !lm<=lmax2 3049 end if ! abs(wan_lm_fac) > 0.0d0 3050 end do !ilmn=1, lmn_size 3051 ABI_FREE(int_rad) 3052 end do !iatom 3053 end do !iwan 3054 end do !isppol 3055 3056 !Deallocate quantities 3057 ABI_FREE(ylmr_fac) 3058 3059 end subroutine mlwfovlp_projpaw
m_mlwfovlp/mlwfovlp_pw [ Functions ]
[ Top ] [ m_mlwfovlp ] [ Functions ]
NAME
mlwfovlp_pw
FUNCTION
Routine which computes PW part of overlap M_{mn}(k,b)
for Wannier code (www.wannier.org f90 version).
INPUTS
cg(2,mpw*nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions. g1(3,nkpt,nntot) = G vector shift which is necessary to obtain k1+b iwav(mband,nkpt,nsppol): shift for pw components in cg. kg(3,mpw*mkmem)=reduced planewave coordinates. mband=maximum number of bands mgfft=maximum size of 1D FFTs mkmem =number of k points treated by this node. mpi_enreg=information about MPI parallelization mpw=maximum dimensioned size of npw. nfft=(effective) number of FFT grid points (for this processor) (see NOTES at beginning of scfcv) ngfft(18)=contain all needed information about 3D FFT (see NOTES at beginning of scfcv) nkpt=number of k points. npwarr(nkpt)=number of planewaves in basis at this k point nspinor=number of spinorial components of the wavefunctions nsppol=1 for unpolarized, 2 for spin-polarized ovikp(nkpt,nntot)= gives nntot value of k2 (in the BZ) for each k1 (k2=k1+b mod(G)) seed_name= seed_name of files containing cg for all k-points to be used with MPI
OUTPUT
cm1(2,mband,mband,nntot,nkpt,nsppol): overlap <u_(nk1)|u_(mk1+b)>.
SIDE EFFECTS
(only writing, printing)
NOTES
SOURCE
1885 subroutine mlwfovlp_pw(cg,cm1,g1,iwav,kg,mband,mkmem,mpi_enreg,mpw,nfft,ngfft,nkpt,nntot,& 1886 & npwarr,nspinor,nsppol,ovikp,seed_name) 1887 1888 !Arguments ------------------------------------ 1889 !scalars 1890 integer,intent(in) :: mband,mkmem,mpw,nfft,nkpt,nntot 1891 integer,intent(in) :: nspinor,nsppol 1892 character(len=fnlen) :: seed_name !seed names of files containing cg info used in case of MPI 1893 type(MPI_type),intent(in) :: mpi_enreg 1894 !arrays 1895 integer,intent(in) :: g1(3,nkpt,nntot),kg(3,mpw*mkmem),ngfft(18),npwarr(nkpt) 1896 integer,intent(in) :: iwav(mband,nkpt,nsppol) 1897 integer,intent(in) :: ovikp(nkpt,nntot) 1898 real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol) 1899 real(dp),intent(out) :: cm1(2,mband,mband,nntot,nkpt,nsppol) 1900 1901 !Local variables------------------------------- 1902 !scalars 1903 integer :: iband1,iband2,ierr,ig,ig1,ig1b,ig2,ig2b 1904 integer :: ig3,ig3b,igk1,igk2,igks1,igks2,ii,ikg,ikpt,ikpt1,ikpt2,imntot,index,intot,ios,ipw 1905 integer :: ispinor,isppol,iunit,me,n1,n2,n3,npoint,npoint2,npw_k,npw_k2 1906 integer :: nprocs,spaceComm 1907 integer,allocatable::indpwk(:,:),kg_k(:,:) 1908 integer,allocatable :: invpwk(:,:) 1909 character(len=500) :: message 1910 character(len=fnlen) :: cg_file !file containing cg info used in case of MPI 1911 logical::lfile 1912 real(dp),allocatable :: cg_read(:,:) !to be used in case of MPI 1913 1914 !************************************************************************ 1915 1916 write(message, '(a,a)' ) ch10,& 1917 & '** mlwfovlp_pw : compute pw part of overlap' 1918 call wrtout(std_out, message,'COLL') 1919 1920 !initialize flags 1921 lfile=.false. 1922 !mpi initialization 1923 spaceComm=MPI_enreg%comm_cell 1924 nprocs=xmpi_comm_size(spaceComm) 1925 me=MPI_enreg%me_kpt 1926 1927 if(nprocs>1) then 1928 ABI_MALLOC(cg_read,(2,nspinor*mpw*mband)) 1929 end if 1930 1931 1932 !****************compute intermediate quantities (index, shifts) ****** 1933 !------------compute index for g points-------------------------------- 1934 !ig is a plane waves which belongs to the sphere ecut for ikpt (they 1935 !are npwarr(ikpt)) 1936 !npoint is the position in the grid of planes waves 1937 !(they are nfft) 1938 !indpwk is a application ig-> npoint 1939 !invpwk is not an application (some npoint have no ig corresponding) 1940 !cg are ordered with npw_k ! 1941 !---------------------------------------------------------------------- 1942 !------------compute index for g points-------------------------------- 1943 !---------------------------------------------------------------------- 1944 write(message, '(a,a)' ) ch10,& 1945 & ' first compute index for g-points' 1946 call wrtout(std_out, message,'COLL') 1947 ! 1948 !Allocations 1949 ABI_MALLOC(kg_k,(3,mpw)) 1950 ABI_MALLOC(indpwk,(nkpt,mpw)) 1951 ABI_MALLOC(invpwk,(nkpt,nfft)) 1952 ! 1953 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 1954 invpwk=0 1955 indpwk=0 1956 kg_k=0 1957 do isppol=1,1 !invpwk is not spin dependent 1958 ! so we just do it once 1959 ikg=0 1960 do ikpt=1,nkpt 1961 ! 1962 ! MPI:cycle over k-points not treated by this node 1963 ! 1964 if ( ABS(MPI_enreg%proc_distrb(ikpt,1,isppol)-me) /=0) CYCLE 1965 1966 ! 1967 ! write(std_out,*)'me',me,'ikpt',ikpt,'isppol',isppol 1968 do npoint=1,nfft 1969 if(invpwk(ikpt,npoint)/=0 )then 1970 write(std_out,*) "error0 , invpwk is overwritten" 1971 write(std_out,*) ikpt,npoint 1972 ABI_ERROR("Aborting now") 1973 end if 1974 end do 1975 npw_k=npwarr(ikpt) 1976 ! write(std_out,*) ikpt,npw_k,nfft 1977 kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 1978 do ig=1,npw_k 1979 if(ig.gt.mpw) then 1980 write(std_out,*)"error ig",ig,"greater than mpw ",mpw 1981 ABI_ERROR("Aborting now") 1982 end if 1983 if(indpwk(ikpt,ig)/=0) then 1984 write(std_out,*) "error, indpwk is overwritten" 1985 write(std_out,*) ikpt,ig,indpwk(ikpt,ig) 1986 ABI_ERROR("Aborting now") 1987 end if 1988 ig1=modulo(kg_k(1,ig),n1) 1989 ig2=modulo(kg_k(2,ig),n2) 1990 ig3=modulo(kg_k(3,ig),n3) 1991 indpwk(ikpt,ig)=ig1+1+n1*(ig2+n2*ig3) 1992 npoint=indpwk(ikpt,ig) 1993 if(npoint.gt.nfft) then 1994 ABI_ERROR("error npoint") 1995 end if 1996 ! write(std_out,*) ikpt,ig,npoint,invpwk(ikpt,npoint) 1997 if(invpwk(ikpt,npoint)/=0) then 1998 write(std_out,*) "error, invpwk is overwritten" 1999 write(std_out,*) ikpt,ig,npoint,invpwk(ikpt,npoint) 2000 ABI_ERROR("Aborting now") 2001 end if 2002 invpwk(ikpt,npoint)=ig 2003 ! write(std_out,*)'ikpt,npoint,invpwk',ikpt,npoint,invpwk(ikpt,npoint) 2004 ! if(ikpt.eq.1) write(std_out,*) "ig npoint",ig, npoint 2005 ! write(std_out,*) "ikpt ig npoint",ikpt,ig, npoint 2006 end do 2007 ikg=ikg+npw_k 2008 2009 end do !ikpt 2010 end do !isppol 2011 !write(std_out,*) "index for g points has been computed" 2012 2013 call xmpi_barrier(spaceComm) 2014 call xmpi_sum(invpwk,spaceComm,ierr) 2015 2016 !---------------------------------------------------------------------- 2017 !------------test invpwk----------------------------------------------- 2018 !---------------------------------------------------------------------- 2019 !write(std_out,*) "TEST INVPWK" 2020 !ikpt=3 2021 !isppol=1 2022 !do ig=1,npwarr(ikpt) 2023 !npoint=indpwk(ikpt,ig) 2024 !write(std_out,*) "ig npoint ",ig, npoint 2025 !write(std_out,*) "ig npoint inv",invpwk(ikpt,npoint),npoint 2026 !end do 2027 !do ig3=1,n3 2028 !do ig2=1,n2 2029 !do ig1=1,n1 2030 !npoint=ig1+(ig2-1)*n1+(ig3-1)*n2*n1 2031 !ig=invpwk(ikpt,npoint) 2032 !! if(ig/=0) write(std_out,*) "ig npoint",ig, npoint 2033 !end do 2034 !end do 2035 !end do 2036 2037 !Deallocate unused variables 2038 ABI_FREE(kg_k) 2039 ABI_FREE(indpwk) 2040 2041 !*********************************************************************** 2042 !**calculate overlap M_{mn}(k,b)=<\Psi_{k,m}|e^{-ibr}|\Psi_{k+b,n}>***** 2043 !*********************************************************************** 2044 write(message, '(a,a)' ) ch10,& 2045 & ' mlwfovlp_pw : compute overlaps ' 2046 call wrtout(std_out, message,'COLL') 2047 write(message, '(a,a)' ) ch10,& 2048 & " nkpt nntot mband " 2049 call wrtout(std_out, message,'COLL') 2050 write(message, '(i6,2x,i6,2x,i6,2x,i6)' ) & 2051 & nkpt,nntot,mband 2052 call wrtout(std_out, message,'COLL') 2053 cm1=zero 2054 write(message, '(a)' ) ' ' 2055 call wrtout(std_out, message,'COLL') 2056 do isppol=1,nsppol 2057 imntot=0 2058 do ikpt1=1,nkpt 2059 ! 2060 ! MPI:cycle over k-points not treated by this node 2061 ! 2062 if ( ABS(MPI_enreg%proc_distrb(ikpt1,1,isppol)-me) /=0) CYCLE 2063 ! 2064 write(message, '(a,i6,a,i6,a,i6)' ) & 2065 & ' Processor',me,' computes k-point',ikpt1,' and spin=',isppol 2066 call wrtout(std_out, message,'COLL') 2067 ! write(std_out,*)trim(message) 2068 2069 do intot=1,nntot 2070 lfile=.false. !flag to know if this kpt will be read from a file, see below 2071 imntot=imntot+1 2072 ikpt2= ovikp(ikpt1,intot) 2073 ! write(std_out,*)'me',me,'ikpt1',ikpt1,'ikpt2',ikpt2,'intot',intot,'isppol',isppol 2074 2075 ! 2076 ! MPI: if ikpt2 not found in this processor then 2077 ! read info from an unformatted file 2078 ! 2079 if ( ABS(MPI_enreg%proc_distrb(ikpt2,1,isppol)-me) /=0) then 2080 lfile=.true. 2081 write(cg_file,'(a,I5.5,".",I1)') trim(seed_name),ikpt2,isppol 2082 iunit=1000+ikpt2+ikpt2*(isppol-1) 2083 npw_k2=npwarr(ikpt2) 2084 open (unit=iunit, file=cg_file,form='unformatted',status='old',iostat=ios) 2085 if(ios /= 0) then 2086 write(message,*) " mlwfovlp_pw: file",trim(cg_file), "not found" 2087 ABI_ERROR(message) 2088 end if 2089 ! 2090 do iband2=1,mband 2091 do ipw=1,npw_k2*nspinor 2092 index=ipw+(iband2-1)*npw_k2*nspinor 2093 read(iunit) (cg_read(ii,index),ii=1,2) 2094 ! if(me==0 .and. ikpt2==4)write(300,*)'ipw,iband2,index',ipw,iband2,index,cg_read(:,index) 2095 ! if(me==1 .and. ikpt2==4)write(301,*)'ipw,iband2,index',ipw,iband2,index,cg_read(:,index) 2096 end do 2097 end do 2098 close(iunit) 2099 end if 2100 ! 2101 npw_k=npwarr(ikpt1) 2102 npw_k2=npwarr(ikpt2) 2103 do ig3=1,n3 2104 do ig2=1,n2 2105 do ig1=1,n1 2106 ! write(std_out,*) isppol,ikpt1,iband1,iband2,intot 2107 npoint=ig1+(ig2-1)*n1+(ig3-1)*n2*n1 2108 if(npoint.gt.nfft) then 2109 write(std_out,*) "error npoint b" 2110 ABI_ERROR("Aborting now") 2111 end if 2112 ig1b=ig1+g1(1,ikpt1,intot) 2113 ig2b=ig2+g1(2,ikpt1,intot) 2114 ig3b=ig3+g1(3,ikpt1,intot) 2115 ! write(std_out,*) ig1,ig2,ig3 2116 ! write(std_out,*) ig1b,ig2b,ig3b 2117 if(ig1b.lt.1) ig1b=ig1b+n1 2118 if(ig2b.lt.1) ig2b=ig2b+n2 2119 if(ig3b.lt.1) ig3b=ig3b+n3 2120 if(ig1b.gt.n1) ig1b=ig1b-n1 2121 if(ig2b.gt.n2) ig2b=ig2b-n2 2122 if(ig3b.gt.n3) ig3b=ig3b-n3 2123 npoint2=ig1b+(ig2b-1)*n1+(ig3b-1)*n2*n1 2124 if(npoint2.gt.nfft) then 2125 write(std_out,*)"error npoint c" 2126 ABI_ERROR("Aborting now") 2127 end if 2128 igk1=invpwk(ikpt1,npoint) 2129 igk2=invpwk(ikpt2,npoint2) 2130 2131 ! if(intot==10) write(std_out,*)'Before igk1 and igk2',ikpt1,ikpt2,isppol 2132 2133 if(igk1/=0.and.igk2/=0) then 2134 do iband2=1,mband 2135 do iband1=1,mband 2136 do ispinor=1,nspinor 2137 ! igks1= (igk1*nspinor)-(nspinor-ispinor) 2138 ! igks2= (igk2*nspinor)-(nspinor-ispinor) 2139 igks1= igk1+ (ispinor-1)*npw_k 2140 igks2= igk2+ (ispinor-1)*npw_k2 2141 2142 ! Here the igks is to include the spinor component missing in igk 2143 if(lfile) index=igks2+npw_k2*nspinor*(iband2-1) !In case of MPI, see below 2144 ! 2145 ! If MPI sometimes the info was read from an unformatted file 2146 ! If that is the case lfile==.true. 2147 ! 2148 ! TODO: this filter should be outside, not inside 1000 loops!!! 2149 if(lfile) then 2150 cm1(1,iband1,iband2,intot,ikpt1,isppol)=cm1(1,iband1,iband2,intot,ikpt1,isppol)+ & 2151 & cg(1,igks1+iwav(iband1,ikpt1,isppol))*cg_read(1,index)& 2152 & + cg(2,igks1+iwav(iband1,ikpt1,isppol))*cg_read(2,index) 2153 cm1(2,iband1,iband2,intot,ikpt1,isppol)=cm1(2,iband1,iband2,intot,ikpt1,isppol)+ & 2154 & cg(1,igks1+iwav(iband1,ikpt1,isppol))*cg_read(2,index)& 2155 & - cg(2,igks1+iwav(iband1,ikpt1,isppol))*cg_read(1,index) 2156 ! 2157 else 2158 cm1(1,iband1,iband2,intot,ikpt1,isppol)=cm1(1,iband1,iband2,intot,ikpt1,isppol)+ & 2159 & cg(1,igks1+iwav(iband1,ikpt1,isppol))*cg(1,igks2+iwav(iband2,ikpt2,isppol))& 2160 & + cg(2,igks1+iwav(iband1,ikpt1,isppol))*cg(2,igks2+iwav(iband2,ikpt2,isppol)) 2161 cm1(2,iband1,iband2,intot,ikpt1,isppol)=cm1(2,iband1,iband2,intot,ikpt1,isppol)+ & 2162 & cg(1,igks1+iwav(iband1,ikpt1,isppol))*cg(2,igks2+iwav(iband2,ikpt2,isppol))& 2163 & - cg(2,igks1+iwav(iband1,ikpt1,isppol))*cg(1,igks2+iwav(iband2,ikpt2,isppol)) 2164 end if 2165 end do !ispinor 2166 end do ! iband1 2167 end do ! iband2 2168 end if 2169 end do ! ig1 2170 end do ! ig2 2171 end do ! ig3 2172 end do ! intot 2173 end do ! ikpt1 2174 end do ! isppol 2175 ! 2176 !Deallocations 2177 ! 2178 ABI_FREE(invpwk) 2179 if(nprocs>1) then 2180 ABI_FREE(cg_read) 2181 end if 2182 2183 end subroutine mlwfovlp_pw
m_mlwfovlp/mlwfovlp_radial [ Functions ]
[ Top ] [ m_mlwfovlp ] [ Functions ]
NAME
mlwfovlp_radial
FUNCTION
Calculates the radial part of the initial functions given to Wannier90 as an starting point for the minimization. The trial functions are a set of solutions to the radial part of the hydrogenic Schrodinger equation as it is explained in Table 3.3 of the Wannier90 user guide.
INPUTS
alpha= Z/a = zona lmax= maximum value of l rvalue= integer defining the choice for radial functions R(r). It can take values from 1-3. It is associted to the radial part of the hydrogenic Schrodinger equation for l=0, See the manual of Wannier90 for more information. (www.wannier.org) xx= scalar number used to calculate the spherical bessel function. J_il(xx)
OUTPUT
mlwfovlp_radial= radial part for initial projections used to construct MLWF
SIDE EFFECTS
None
NOTES
Calculates the radial part of the initial functions given as an initial guess by the user to construct the MLWF.
SOURCE
3093 subroutine mlwfovlp_radial(alpha,lmax,lmax2,radial,rvalue,xx) 3094 3095 !Arguments ------------------------------------ 3096 !scalars 3097 integer,intent(in) :: lmax,lmax2,rvalue 3098 real(dp),intent(in) :: alpha,xx 3099 !arrays 3100 real(dp),intent(out) :: radial(lmax2) 3101 3102 !Local variables 3103 !scalars 3104 integer :: ir,ll,lm,mesh,mm 3105 real(dp),parameter :: dx=0.015d0,rmax=10.d0,xmin=0.d0 3106 real(dp) :: aa,ftmp,gauss,rtmp,x 3107 character(len=500) :: message 3108 !arrays 3109 real(dp),save :: dblefact(4)=(/1_dp,3_dp,15_dp,105_dp/) 3110 real(dp),allocatable :: aux(:),bes(:),cosr(:),func_r(:),r(:),rad_int(:) 3111 real(dp),allocatable :: sinr(:) 3112 3113 ! ************************************************************************* 3114 3115 !Radial functions in the form of hydrogenic orbitals as defined in the 3116 !wannier90 manual. 3117 if(( rvalue > 0 ).and.(rvalue < 4)) then 3118 3119 ! mesh 3120 mesh= nint((rmax - xmin ) / dx + 1) 3121 ABI_MALLOC( bes,(mesh)) 3122 ABI_MALLOC(func_r,(mesh)) 3123 ABI_MALLOC(r,(mesh)) 3124 ABI_MALLOC(rad_int,(mesh)) 3125 ABI_MALLOC( aux,(mesh)) 3126 ABI_MALLOC(cosr,(mesh)) 3127 ABI_MALLOC(sinr,(mesh)) 3128 do ir=1, mesh 3129 x=xmin+DBLE(ir-1)*dx 3130 r(ir)=x 3131 end do !ir 3132 3133 ! radial functions shown in table 3.3 of wannier90 manual 3134 if (rvalue==1) func_r(:) = 2.d0 * alpha**(3.d0/2.d0) * exp(-alpha*r(:)) 3135 if (rvalue==2) func_r(:) = 1.d0/(2.d0*sqrt(2.d0))*alpha**(3.d0/2.d0) *& 3136 & (2.d0 - alpha*r(:))*exp(-alpha*r(:)/2.d0) 3137 if (rvalue==3) func_r(:) = sqrt(4.d0/27.d0)*alpha**(3.d0/2.d0)& 3138 & * (1.d0 - 2.d0*alpha*r(:)/3.d0 + 2.d0*alpha**2*r(:)**2/27.d0)& 3139 & * exp(-alpha * r(:)/3.d0) 3140 3141 ! compute spherical bessel functions 3142 cosr(:)=cos(xx*r(:)) 3143 sinr(:)=sin(xx*r(:)) 3144 lm=0 3145 do ll=0,lmax 3146 call besjm(xx,bes,cosr,ll,mesh,sinr,r) 3147 aux(:)=bes(:)*func_r(:)*r(:) 3148 ! do ir=1,mesh 3149 ! write(310,*) r(ir),bes(ir) 3150 ! end do 3151 call simpson_int(mesh,dx,aux,rad_int) 3152 rtmp=rad_int(mesh)/mesh 3153 do mm=-ll,ll 3154 lm=lm+1 3155 radial(lm)=rtmp 3156 end do !mm 3157 end do !ll 3158 ABI_FREE(bes) 3159 ABI_FREE(func_r) 3160 ABI_FREE(r) 3161 ABI_FREE(aux) 3162 ABI_FREE(rad_int) 3163 ABI_FREE(cosr) 3164 ABI_FREE(sinr) 3165 3166 ! Radial part in the form of Gaussian functions of a given width 3167 ! Taken by code of made by drh. 3168 elseif ( rvalue == 4) then 3169 aa=1._dp/alpha 3170 gauss=exp(-0.25_dp*(aa*xx)**2) 3171 lm=0 3172 do ll=0,lmax 3173 ftmp=(0.5_dp*pi)**(0.25_dp)*aa*sqrt(aa/dblefact(ll+1))*(aa*xx)**ll*gauss 3174 do mm=-ll,ll 3175 lm=lm+1 3176 radial(lm)=ftmp 3177 end do 3178 end do 3179 else ! rvalue < 0 of rvalue > 4 3180 write(message,'(a,i6,5a)')& 3181 & ' Radial function r=',rvalue,ch10,& 3182 & ' is not defined',ch10,& 3183 & ' Modify .win file',ch10 3184 ABI_BUG(message) 3185 end if !rvalue 3186 3187 end subroutine mlwfovlp_radial
m_mlwfovlp/mlwfovlp_seedname [ Functions ]
[ Top ] [ m_mlwfovlp ] [ Functions ]
NAME
mlwfovlp_seedname
FUNCTION
Get seed name and file names of all wannier90 related files
INPUTS
fname_w90=root name of file appended with _w90
OUTPUT
filew90_win= main input file for Wannier90 filew90_wout= main output file for Wannier90 filew90_amn= file containing Amn matrix filew90_ramn= file containing Amn matrix (random initial projections) filew90_mmn= file containing Mmn matrix filew90_eig= file containing eigenvalues nsppol= number of spin polarizations seed_name= common seed name for all wannier90 related files
SIDE EFFECTS
NOTES
SOURCE
1414 subroutine mlwfovlp_seedname(fname_w90,filew90_win,filew90_wout,filew90_amn,& 1415 & filew90_ramn,filew90_mmn,filew90_eig,nsppol,seed_name) 1416 1417 !Arguments ------------------------------------ 1418 integer,intent(in) :: nsppol 1419 character(len=fnlen),intent(out) :: filew90_win(nsppol),filew90_wout(nsppol),filew90_amn(nsppol),filew90_ramn(nsppol) 1420 character(len=fnlen),intent(out) :: filew90_mmn(nsppol),filew90_eig(nsppol),seed_name(nsppol) 1421 character(len=fnlen),intent(in) :: fname_w90 1422 1423 !Local variables------------------------------- 1424 integer::isppol 1425 character(len=fnlen) :: test_win1,test_win2,test_win3 1426 logical :: lfile 1427 character(len=2000) :: message 1428 character(len=10)::postfix 1429 ! ************************************************************************* 1430 1431 seed_name(:)=trim(fname_w90) 1432 do isppol=1,nsppol 1433 if(nsppol==1)postfix='.win' 1434 if(nsppol==2 .and. isppol==1)postfix='_up.win' 1435 if(nsppol==2 .and. isppol==2)postfix='_down.win' 1436 ! 1437 filew90_win(isppol)=trim(seed_name(isppol))//trim(postfix) 1438 test_win1=filew90_win(isppol) 1439 inquire(file=filew90_win(isppol),exist=lfile) 1440 1441 if(.not.lfile) then 1442 seed_name(isppol)='wannier90' 1443 filew90_win(isppol)=trim(seed_name(isppol))//trim(postfix) 1444 test_win2=filew90_win(isppol) 1445 inquire(file=filew90_win(isppol),exist=lfile) 1446 end if 1447 1448 if(.not.lfile) then 1449 seed_name(isppol)='w90' 1450 filew90_win=trim(seed_name(isppol))//trim(postfix) 1451 test_win3=filew90_win(isppol) 1452 inquire(file=filew90_win(isppol),exist=lfile) 1453 end if 1454 1455 if(.not.lfile) then 1456 write(message,'(17a)')ch10,& 1457 & ' mlwfovlp_seedname : ERROR - ',ch10,& 1458 & ' wannier90 interface needs one of the following files:',ch10,& 1459 & ' ',trim(test_win1),ch10,& 1460 & ' ',trim(test_win2),ch10,& 1461 & ' ',trim(test_win3),ch10,& 1462 & ' Action: read wannier90 tutorial and/or user manual',ch10,& 1463 & ' and supply proper *.win file' 1464 ABI_ERROR(message) 1465 end if 1466 end do !isppol 1467 1468 1469 !Files having different names for 1470 !different spin polarizations 1471 if(nsppol==1) then 1472 filew90_win(1) =trim(seed_name(1))//'.win' 1473 filew90_wout(1)=trim(seed_name(1))//'.wout' 1474 filew90_ramn(1)=trim(seed_name(1))//'random.amn' 1475 filew90_amn(1) =trim(seed_name(1))//'.amn' 1476 filew90_mmn(1) =trim(seed_name(1))//'.mmn' 1477 filew90_eig(1) =trim(seed_name(1))//'.eig' 1478 elseif(nsppol==2) then 1479 filew90_win(1) =trim(seed_name(1))//'_up.win' 1480 filew90_win(2) =trim(seed_name(2))//'_down.win' 1481 ! 1482 filew90_wout(1)=trim(seed_name(1))//'_up.wout' 1483 filew90_wout(2)=trim(seed_name(2))//'_down.wout' 1484 ! 1485 filew90_ramn(1)=trim(seed_name(1))//'random_up.amn' 1486 filew90_ramn(2)=trim(seed_name(2))//'random_down.amn' 1487 ! 1488 filew90_amn(1)=trim(seed_name(1))//'_up.amn' 1489 filew90_amn(2)=trim(seed_name(2))//'_down.amn' 1490 ! 1491 filew90_mmn(1)=trim(seed_name(1))//'_up.mmn' 1492 filew90_mmn(2)=trim(seed_name(2))//'_down.mmn' 1493 ! 1494 filew90_eig(1)=trim(seed_name(1))//'_up.eig' 1495 filew90_eig(2)=trim(seed_name(2))//'_down.eig' 1496 end if 1497 !change also seed_name for nsppol=2 1498 if(nsppol==2) then 1499 seed_name(1)=trim(seed_name(1))//'_up' 1500 seed_name(2)=trim(seed_name(2))//'_down' 1501 end if 1502 !End file-name section 1503 1504 write(message, '(a,a)' ) ch10,& 1505 & '---------------------------------------------------------------' 1506 call wrtout(ab_out,message,'COLL') 1507 call wrtout(std_out, message,'COLL') 1508 write(message, '(5a)' ) ch10,& 1509 & ' Calculation of overlap and call to wannier90 library ',ch10,& 1510 & ' to obtain maximally localized wannier functions ',ch10 1511 1512 call wrtout(std_out, message,'COLL') 1513 call wrtout(ab_out,message,'COLL') 1514 1515 if(nsppol==1) then 1516 write(message, '(23a)' ) & 1517 & ' - ',trim(filew90_win(1)),' is a mandatory secondary input',ch10,& 1518 & ' - ',trim(filew90_wout(1)),' is the output for the library',ch10,& 1519 & ' - ',trim(filew90_ramn(1)),' contains random projections',ch10,& 1520 & ' - ',trim(filew90_amn(1)),' contains projections',ch10,& 1521 & ' - ',trim(filew90_mmn(1)),' contains the overlap',ch10,& 1522 & ' - ',trim(filew90_eig(1)),' contains the eigenvalues' 1523 elseif(nsppol==2) then 1524 write(message, '(41a)' ) & 1525 & ' - ',trim(filew90_win(1)),& 1526 & ' and ',trim(filew90_win(2)),ch10,'are mandatory secondary input',ch10,& 1527 & ' - ',trim(filew90_wout(1)),& 1528 & ' and ',trim(filew90_wout(2)),ch10,' are the output for the library',ch10,& 1529 & ' - ',trim(filew90_ramn(1)),& 1530 & ' and ',trim(filew90_ramn(2)),ch10,' contain random projections',ch10,& 1531 & ' - ',trim(filew90_amn(1)),& 1532 & ' and ',trim(filew90_amn(2)),ch10,' contain projections',ch10,& 1533 & ' - ',trim(filew90_mmn(1)),& 1534 & ' and ',trim(filew90_mmn(2)),ch10,' contain the overlap',ch10,& 1535 & ' - ',trim(filew90_eig(1)),& 1536 & ' and ',trim(filew90_eig(2)),ch10,' contain the eigenvalues' 1537 end if 1538 call wrtout(std_out, message,'COLL') 1539 call wrtout(ab_out,message,'COLL') 1540 1541 write(message, '(a,a)' ) ch10,& 1542 & '---------------------------------------------------------------' 1543 call wrtout(ab_out,message,'COLL') 1544 call wrtout(std_out, message,'COLL') 1545 1546 end subroutine mlwfovlp_seedname
m_mlwfovlp/mlwfovlp_setup [ Functions ]
[ Top ] [ m_mlwfovlp ] [ Functions ]
NAME
mlwfovlp_setup
FUNCTION
Routine which creates table g1 and ovikp necessary to compute overlap for Wannier code (www.wannier.org f90 version).
INPUTS
atom_symbols(natom)= table of symbol for each atom
and each |p_lmn> non-local projector
dtset <type(dataset_type)>=all input variables for this dataset
filew90_win(nsppol) secondary input files for w90
lwanniersetup= flag: only 1 is fully working.
natom =number of atoms in cell.
mband=maximum number of bands
natom=number of atoms in cell.
nkpt=number of k points.
num_bands(isppol)=number of bands actually used to construct the wannier function
nwan(isppol)= number of wannier fonctions (read in wannier90.win).
dtset <type(dataset_type)>=all input variables for this dataset
real_lattice(3,3)=dimensional primitive translations for real space
in format required by wannier90
recip_lattice(3,3)=dimensional primitive translations for reciprocal space
in format required by wannier90
rprimd(3,3)=dimensional primitive translations for real space (bohr)
seed_name=character string for generating wannier90 filenames
xcart(3,natom)=atomic coordinates in bohr
xred(3,natom)=reduced dimensionless atomic coordinates
OUTPUT
band_in(mband,nsppol) = band to take into account for wannier calculation g1(3,nkpt,nntot) = G vector shift which is necessary to obtain k1+b from k2 in the case where k1+b does not belong to the 1st BZ. nband_inc(nsppol) = # of included bands nntot = number of k-point neighbour ovikp(nkpt,nntot)= gives nntot value of k2 (in the BZ) for each k1 (k2=k1+b mod(G))
SIDE EFFECTS
(only writing, printing)
NOTES
SOURCE
1594 subroutine mlwfovlp_setup(atom_symbols,band_in,dtset,filew90_win,gamma_only,& 1595 & g1,lwanniersetup,mband,natom,nband_inc,nkpt,& 1596 & nntot,num_bands,num_nnmax,nsppol,nwan,ovikp,& 1597 & proj_l,proj_m,proj_radial,proj_site,proj_s_loc,proj_s_qaxis_loc,proj_x,proj_z,proj_zona,& 1598 & real_lattice,recip_lattice,rprimd,seed_name,spinors,xcart,xred) 1599 1600 !Arguments--------------------------- 1601 ! scalars 1602 !scalars 1603 integer,intent(in) :: lwanniersetup,mband,natom,nkpt,nsppol 1604 integer,intent(in) :: num_nnmax 1605 integer,intent(out) :: nband_inc(nsppol),nntot,num_bands(nsppol),nwan(nsppol) 1606 logical,intent(in) :: gamma_only,spinors 1607 type(dataset_type),intent(in) :: dtset 1608 !arrays 1609 integer,intent(out) :: g1(3,nkpt,num_nnmax),ovikp(nkpt,num_nnmax) 1610 integer,intent(out) :: proj_l(mband,nsppol),proj_m(mband,nsppol),proj_radial(mband,nsppol) 1611 real(dp),intent(in) :: real_lattice(3,3) 1612 real(dp),intent(in) :: recip_lattice(3,3),rprimd(3,3),xred(3,natom) 1613 real(dp),intent(out) :: proj_site(3,mband,nsppol),proj_x(3,mband,nsppol),proj_z(3,mband,nsppol) 1614 real(dp),intent(out) :: proj_zona(mband,nsppol),xcart(3,natom) 1615 logical,intent(out) :: band_in(mband,nsppol) 1616 character(len=3),intent(out) :: atom_symbols(natom) 1617 character(len=fnlen),intent(in) :: seed_name(nsppol),filew90_win(nsppol) 1618 1619 integer, optional, intent(out) :: proj_s_loc(mband) 1620 real(dp), optional, intent(out) :: proj_s_qaxis_loc(3,mband) 1621 1622 !Local variables--------------------------- 1623 !scalars 1624 integer :: iatom,icb,ikpt,ikpt1,intot,isppol,itypat,jj,mband_,unt 1625 real(dp) :: znucl1 1626 character(len=2) :: symbol 1627 character(len=500) :: message 1628 character(len=fnlen) :: filew90_nnkp 1629 type(atomdata_t) :: atom 1630 !arrays 1631 integer :: exclude_bands(mband,nsppol),ngkpt(3) 1632 1633 ! ************************************************************************* 1634 1635 !^^^^^^^^^^^^^^^^read wannier90.nnkp^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 1636 if(lwanniersetup==0) then !this part is not coded for nsppol>1 1637 isppol=1 1638 filew90_nnkp=trim(seed_name(isppol))//'.nnkp' 1639 if (open_file(filew90_nnkp,message,newunit=unt,form='formatted',status='old') /= 0) then 1640 ABI_ERROR(message) 1641 end if 1642 read(unt,*) 1643 read(unt,*) nntot , mband_, nwan(1) 1644 write(message, '(a,a,i6,i6,i6)' )ch10,& 1645 & ' mlwfovlp_setup nntot,mband,nwan ', nntot,mband_,nwan(1) 1646 call wrtout(std_out,message,'COLL') 1647 if(mband_.ne.mband) then 1648 write(message, '(4a)' )& 1649 & 'mband_ is not equal to mband ',ch10,& 1650 & 'Action: check ',trim(filew90_nnkp) 1651 ABI_ERROR(message) 1652 end if 1653 if(nwan(1)>mband) then 1654 write(message, '(4a)' )& 1655 & 'nwan > mband ',ch10,& 1656 & 'Action: check ',trim(filew90_nnkp) 1657 ABI_ERROR(message) 1658 end if 1659 if(nwan(1)==0) then 1660 write(message, '(4a)' )& 1661 & 'nwan = 0 ',ch10,& 1662 & 'Action: check ',trim(filew90_nnkp) 1663 ABI_ERROR(message) 1664 end if 1665 do ikpt=1,nkpt 1666 do intot=1,nntot 1667 ! ikpt1: k point (ikpt=ikpt1) 1668 ! ovikp(intot,ikpt): neighbour number intot for ikpt 1669 ! g1(1:3,intot,ikpt): non reciprocal space vector between the 2 k-points 1670 read(unt,*) & 1671 & ikpt1,ovikp(ikpt,intot),(g1(jj,ikpt,intot),jj=1,3) 1672 if(ikpt1.ne.ikpt) write(std_out,*) "warning: ikpt1 .ne ikpt : ?" 1673 end do 1674 end do 1675 close(unt) 1676 write(message, '(3a)' )ch10,& 1677 & trim(filew90_nnkp),'wannier90.nnkp has been read !' 1678 call wrtout(std_out,message,'COLL') 1679 1680 message = ' exclude bands is not given in this case (not implemented) ' 1681 ABI_ERROR(message) 1682 1683 ! ^^^^^^^^^^^^^^^^^^^^^^^ call wannier_setup begin^^^^^^^^^^^^^^^^^^^^^^^^ 1684 else if (lwanniersetup==1) then 1685 num_bands(:)=mband 1686 ! num_nnmax=12 !limit fixed for compact structure in wannier_setup. 1687 ovikp=0.d0 1688 ! "When nshiftk=1, kptrlatt is initialized as a diagonal (3x3) matrix, whose diagonal 1689 ! elements are the three values ngkpt(1:3)" 1690 ngkpt(1)=dtset%kptrlatt(1,1) 1691 ngkpt(2)=dtset%kptrlatt(2,2) ! have to verif kptrlatt is diagonal 1692 ngkpt(3)=dtset%kptrlatt(3,3) 1693 do iatom=1,natom 1694 itypat=dtset%typat(iatom) 1695 znucl1=dtset%znucl(itypat) 1696 call atomdata_from_znucl(atom, znucl1) 1697 symbol=trim(adjustl(atom%symbol)) 1698 ! write(309,*) symbol 1699 atom_symbols(iatom)=symbol 1700 xcart(:,iatom)=rprimd(:,1)*xred(1,iatom)+& 1701 & rprimd(:,2)*xred(2,iatom)+& 1702 & rprimd(:,3)*xred(3,iatom) 1703 end do ! iatom 1704 ! write(std_out,*) xcart 1705 ! write(std_out,*) Bohr_Ang 1706 ! write(std_out,*) rprimd*Bohr_Ang 1707 ! write(std_out,*) seed_name 1708 ! write(std_out,*) ngkpt 1709 ! write(std_out,*) nkpt 1710 ! write(std_out,*) mband 1711 ! write(std_out,*) natom 1712 ! write(std_out,*) atom_symbols 1713 write(message, '(a,a)' )ch10,& 1714 & '** mlwfovlp_setup: call wannier90 library subroutine wannier_setup' 1715 call wrtout(std_out,message,'COLL') 1716 #if defined HAVE_WANNIER90 1717 nwan(:)=0 1718 num_bands(:)=0 1719 do isppol=1,nsppol 1720 #ifdef HAVE_WANNIER90_V1 1721 call wannier_setup(seed_name(isppol),ngkpt,nkpt& !input 1722 & ,real_lattice,recip_lattice,dtset%kpt& !input 1723 & ,mband,natom,atom_symbols,xcart*Bohr_Ang& !input 1724 & ,gamma_only,spinors& !input 1725 & ,nntot,ovikp,g1,num_bands(isppol),nwan(isppol)& !output 1726 & ,proj_site(:,:,isppol),proj_l(:,isppol)& !output 1727 & ,proj_m(:,isppol),proj_radial(:,isppol)& !output 1728 & ,proj_z(:,:,isppol),proj_x(:,:,isppol)& !output 1729 & ,proj_zona(:,isppol),exclude_bands(:,isppol)) !output 1730 #else 1731 !WANNIER90_V2 has the 2 optional arguments 1732 if (present(proj_s_loc)) then 1733 call wannier_setup(seed_name(isppol),ngkpt,nkpt& !input 1734 & ,real_lattice,recip_lattice,dtset%kpt& !input 1735 & ,mband,natom,atom_symbols,xcart*Bohr_Ang& !input 1736 & ,gamma_only,spinors& !input 1737 & ,nntot,ovikp,g1,num_bands(isppol),nwan(isppol)& !output 1738 & ,proj_site(:,:,isppol),proj_l(:,isppol)& !output 1739 & ,proj_m(:,isppol),proj_radial(:,isppol)& !output 1740 & ,proj_z(:,:,isppol),proj_x(:,:,isppol)& !output 1741 & ,proj_zona(:,isppol),exclude_bands(:,isppol)& !output 1742 & ,proj_s_loc,proj_s_qaxis_loc) !output 1743 else 1744 !no proj_s_loc provided 1745 call wannier_setup(seed_name(isppol),ngkpt,nkpt& !input 1746 & ,real_lattice,recip_lattice,dtset%kpt& !input 1747 & ,mband,natom,atom_symbols,xcart*Bohr_Ang& !input 1748 & ,gamma_only,spinors& !input 1749 & ,nntot,ovikp,g1,num_bands(isppol),nwan(isppol)& !output 1750 & ,proj_site(:,:,isppol),proj_l(:,isppol)& !output 1751 & ,proj_m(:,isppol),proj_radial(:,isppol)& !output 1752 & ,proj_z(:,:,isppol),proj_x(:,:,isppol)& !output 1753 & ,proj_zona(:,isppol),exclude_bands(:,isppol)) !output 1754 end if 1755 #endif 1756 end do !isppol 1757 ! if we do not have w90, avoid complaints about unused input variables 1758 #else 1759 ABI_UNUSED(gamma_only) 1760 ABI_UNUSED(real_lattice) 1761 ABI_UNUSED(recip_lattice) 1762 ABI_UNUSED(spinors) 1763 #endif 1764 ! do isppol=1,nsppol 1765 ! write(std_out,*) "1", nntot,nwan(isppol) 1766 ! write(std_out,*) "2",num_bands(isppol) ! states on which wannier functions are computed 1767 ! write(std_out,*) "3", proj_site(:,1:nwan(isppol),isppol) 1768 ! write(std_out,*) "4",proj_l(1:nwan(isppol),isppol) 1769 ! write(std_out,*) "5",proj_m(1:nwan(isppol),isppol) 1770 ! write(std_out,*) "6", proj_radial(1:nwan(isppol),isppol) 1771 ! write(std_out,*) "7", proj_z(:,1:nwan(isppol),isppol) 1772 ! write(std_out,*) "8", proj_x(:,1:nwan(isppol),isppol) 1773 ! write(std_out,*) "9",proj_zona(1:nwan(isppol),isppol) 1774 ! write(std_out,*) "10",exclude_bands(:,isppol) 1775 ! end do!isppol 1776 ! testdebug: ovikp(1,1)=1 1777 ! ^^^^^^^^^^^^^^^^^^^^^^^ end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 1778 end if ! lwanniersetup 1779 do isppol=1,nsppol 1780 band_in(:,isppol)=.true. 1781 do icb=1,mband 1782 if(exclude_bands(icb,isppol).ne.0) band_in(exclude_bands(icb,isppol),isppol)=.false. 1783 end do 1784 nband_inc(isppol)=0 1785 do icb=1, mband 1786 if (band_in(icb,isppol)) then 1787 nband_inc(isppol)=nband_inc(isppol)+1 1788 end if 1789 end do 1790 end do !isppol 1791 if(any(mband.gt.num_bands(:))) then 1792 write(message, '(a,a)' )ch10,& 1793 & ' Following bands are excluded from the calculation of wannier functions:' 1794 call wrtout(std_out,message,'COLL') 1795 1796 do isppol=1,nsppol 1797 if(nsppol==2) then 1798 write(message,'("For spin",i4)')isppol 1799 ! write(message,'(a,i)')'For spin=',isppol 1800 call wrtout(std_out,message,'COLL') 1801 end if !nsppol 1802 do jj=1,mband-num_bands(isppol),10 1803 write(message,'(10i7)') exclude_bands(jj:min(jj+9,mband-num_bands(isppol)),isppol) 1804 call wrtout(std_out,message,'COLL') 1805 end do 1806 end do !isppol 1807 end if 1808 1809 do isppol=1,nsppol 1810 if(nsppol==2) then 1811 write(message,'("For spin",i4)')isppol 1812 call wrtout(std_out,message,'COLL') 1813 end if !nsppol 1814 write(message, '(a,i6,3a)' )ch10,& 1815 & nwan(isppol),' wannier functions will be computed (see ',trim(filew90_win(isppol)),')' 1816 call wrtout(std_out,message,'COLL') 1817 ! write(std_out,*) exclude_bands(icb),band_in(icb) 1818 ! ^^^^^^^^^^^^^^^END OF READING 1819 write(message, '(a,i6,a)' )ch10,& 1820 & num_bands(isppol),' bands will be used to extract wannier functions' 1821 call wrtout(std_out,message,'COLL') 1822 if(num_bands(isppol).lt.nwan(isppol)) then 1823 write(message, '(4a)' )& 1824 & ' number of bands is lower than the number of wannier functions',ch10,& 1825 & ' Action : check input file and ',trim(filew90_win(isppol)) 1826 ABI_ERROR(message) 1827 else if (num_bands(isppol)==nwan(isppol)) then 1828 write(message, '(a,a,a,a)' )ch10,& 1829 & ' Number of bands is equal to the number of wannier functions',ch10,& 1830 & ' Disentanglement will not be necessary' 1831 call wrtout(std_out,message,'COLL') 1832 else if (num_bands(isppol).gt.nwan(isppol)) then 1833 write(message, '(a,a,a,a)' )ch10,& 1834 & ' Number of bands is larger than the number of wannier functions',ch10,& 1835 & ' Disentanglement will be necessary' 1836 call wrtout(std_out,message,'COLL') 1837 end if 1838 write(message, '(2x,a,a,i3,1x,a)' )ch10,& 1839 & ' Each k-point has', nntot,'neighbours' 1840 call wrtout(std_out,message,'COLL') 1841 1842 end do !isppol 1843 1844 end subroutine mlwfovlp_setup
m_mlwfovlp/mlwfovlp_ylmfac [ Functions ]
[ Top ] [ m_mlwfovlp ] [ Functions ]
NAME
mlwfovlp_ylmfac
FUNCTION
Routine that produces a factor by which the initial guess of functions will be multiplied for the Wannier90 interface. It is just used if there are rotations, or if the functions required are linear combinations of the ylm real functions. Example, For a function G(r)= 1/2 s + 1/3 px - 1/2 pz it would produce a matrix of the following form: [1/2,-1/2,1/3,0,0...0] The real spherical harmonics are given as factors of complex spherical harmonics The real spherical harmonics are given in table 3.1 of Wannier90 user guide.
INPUTS
lmax= maximum l value for spherical harmonics lmax2=number of ylm functions mband=maximum number of bands nwan = number of wannier functions 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_x(3,mband)= x axis for the projection. proj_z(3,mband)= z axis for the projection.
OUTPUT
ylmc_fac(lmax2,nwan)=matrix containig a factor for ylm hybrid orbitals
SIDE EFFECTS
(only writing, printing)
NOTES
SOURCE
3229 subroutine mlwfovlp_ylmfac(ylmc_fac,lmax,lmax2,mband,nwan,proj_l,proj_m,proj_x,proj_z) 3230 3231 !Arguments ------------------------------------ 3232 integer, intent(in):: lmax,lmax2,nwan,mband 3233 ! arrays 3234 integer,intent(in) :: proj_l(mband),proj_m(mband) 3235 real(dp),intent(in) :: proj_x(3,mband),proj_z(3,mband) 3236 complex(dp),intent(out)::ylmc_fac(lmax2,nwan) 3237 ! 3238 !Local variables------------------------------- 3239 ! 3240 integer :: orb_idx(16)=(/1,3,4,2,7,8,6,9,5,13,14,12,15,11,16,10/) !Tab3.1 Wannier90 user guide 3241 integer :: idum,ii,info,inversion_flag 3242 integer :: ir,iwan,jj,ll,lm,lmc,mm,mr 3243 real(dp):: onem,test 3244 ! arrays 3245 integer:: ipiv(lmax2) 3246 real(dp)::r(3,lmax2),rp(3,lmax2) 3247 real(dp)::rs2,rs3,rs6,rs12,umat(3,3) 3248 complex(dp)::crot(lmax2,lmax2),ctor(lmax2,lmax2),orb_lm(lmax2,-5:3,7) 3249 complex(dp):: ylmcp(lmax2) 3250 complex(dp):: ylmc_rr(lmax2,lmax2),ylmc_rr_save(lmax2,lmax2) 3251 complex(dp):: ylmc_rrinv(lmax2,lmax2),ylmc_rp(lmax2,lmax2) 3252 complex(dp),parameter :: c0=(0._dp,0._dp),c1=(1._dp,0._dp),ci=(0._dp,1._dp) 3253 character(len=500) :: message ! to be uncommented, if needed 3254 3255 ! ************************************************************************* 3256 3257 3258 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3259 !DEBUG 3260 !write(std_out,*)'lmax ',lmax,'lmax2 ',lmax2 3261 !write(std_out,*)'mband ',mband,'nwan ',nwan 3262 ! 3263 !do iwan=1,nwan 3264 !write(std_out,*)'iwan,proj_l, proj_m',proj_l(iwan),proj_m(iwan) 3265 !write(std_out,*)'iwan,proj_x, proj_z',iwan,proj_x(:,iwan),proj_z(:,iwan) 3266 !end do 3267 !!END DEBUG 3268 3269 !constants for linear combinations of ylm's 3270 rs2=1._dp/sqrt(2._dp) 3271 rs3=1._dp/sqrt(3._dp) 3272 rs6=1._dp/sqrt(6._dp) 3273 rs12=1._dp/sqrt(12._dp) 3274 3275 !complex lm coefficients for real spherical harmonics in conventional order 3276 !s, py,pz,px, dxy,dyz,dz2,dxz,dx2-y2, fy(3x2-y2),fxyz,fyz2,fz3,fxz2, 3277 !fz(x2-y2),fx(x2-3y2) 3278 ctor(:,:)=c0 3279 do ll=0,lmax 3280 mm=0 3281 lm= ll**2+ll+mm+1 3282 ctor(lm,lm)=c1 3283 if(ll>0) then 3284 onem=one 3285 do mm=1,ll 3286 onem=-onem !(-1^mm) 3287 lm= ll**2+ll+mm+1 3288 lmc=ll**2+ll-mm+1 3289 ctor(lm ,lm )=rs2*c1 3290 ctor(lmc,lm )=onem*rs2*c1 3291 ctor(lm ,lmc)=rs2*ci 3292 ctor(lmc,lmc)=-onem*rs2*ci 3293 end do 3294 end if 3295 end do 3296 3297 lm=0 3298 do ll=0,lmax 3299 do mm=-ll,ll 3300 lm=lm+1 3301 ctor(:,lm)=ctor(:,lm)*conjg(ci)**ll 3302 end do !mm 3303 end do !ll 3304 3305 3306 !coefficients for basic wannier orbitals in Table 3.1 order 3307 orb_lm(:,:,:)=c0 3308 ii=0 3309 do ll=0,lmax 3310 do mr=1,2*ll+1 3311 ii=ii+1 3312 orb_lm(:,ll,mr)=ctor(:,orb_idx(ii)) 3313 end do 3314 end do 3315 3316 3317 3318 !coefficients for linear combinations in table 3.2 order 3319 if(lmax>=1) then 3320 ! s px 3321 orb_lm(:,-1,1)=rs2*ctor(:,1)+rs2*ctor(:,4) 3322 orb_lm(:,-1,2)=rs2*ctor(:,1)-rs2*ctor(:,4) 3323 ! s px py 3324 orb_lm(:,-2,1)=rs3*ctor(:,1)-rs6*ctor(:,4)+rs2*ctor(:,2) 3325 orb_lm(:,-2,2)=rs3*ctor(:,1)-rs6*ctor(:,4)-rs2*ctor(:,2) 3326 orb_lm(:,-2,3)=rs3*ctor(:,1)+2._dp*rs6*ctor(:,4) 3327 ! s px py pz 3328 orb_lm(:,-3,1)=half*(ctor(:,1)+ctor(:,4)+ctor(:,2)+ctor(:,3)) 3329 orb_lm(:,-3,2)=half*(ctor(:,1)+ctor(:,4)-ctor(:,2)-ctor(:,3)) 3330 orb_lm(:,-3,3)=half*(ctor(:,1)-ctor(:,4)+ctor(:,2)-ctor(:,3)) 3331 orb_lm(:,-3,4)=half*(ctor(:,1)-ctor(:,4)-ctor(:,2)+ctor(:,3)) 3332 end if 3333 if(lmax>=2) then 3334 ! s px py 3335 orb_lm(:,-4,1)=rs3*ctor(:,1)-rs6*ctor(:,4)+rs2*ctor(:,2) 3336 orb_lm(:,-4,2)=rs3*ctor(:,1)-rs6*ctor(:,4)-rs2*ctor(:,2) 3337 orb_lm(:,-4,3)=rs3*ctor(:,1)+2._dp*rs6*ctor(:,4) 3338 ! pz dz2 3339 orb_lm(:,-4,4)= rs2*ctor(:,3)+rs2*ctor(:,7) 3340 orb_lm(:,-4,5)=-rs2*ctor(:,3)+rs2*ctor(:,7) 3341 ! s px dz2 dx2-y2 3342 orb_lm(:,-5,1)=rs6*ctor(:,1)-rs2*ctor(:,4)-rs12*ctor(:,7)+half*ctor(:,9) 3343 orb_lm(:,-5,2)=rs6*ctor(:,1)+rs2*ctor(:,4)-rs12*ctor(:,7)+half*ctor(:,9) 3344 ! s py dz2 dx2-y2 3345 orb_lm(:,-5,3)=rs6*ctor(:,1)-rs2*ctor(:,2)-rs12*ctor(:,7)-half*ctor(:,9) 3346 orb_lm(:,-5,4)=rs6*ctor(:,1)+rs2*ctor(:,2)-rs12*ctor(:,7)-half*ctor(:,9) 3347 ! s pz dz2 3348 orb_lm(:,-5,5)=rs6*ctor(:,1)-rs2*ctor(:,3)+rs3*ctor(:,7) 3349 orb_lm(:,-5,6)=rs6*ctor(:,1)+rs2*ctor(:,3)+rs3*ctor(:,7) 3350 end if 3351 3352 !stuff complex wannier orbital coefficient array 3353 do iwan=1,nwan 3354 ylmc_fac(:,iwan)=orb_lm(:,proj_l(iwan),proj_m(iwan)) 3355 end do 3356 3357 3358 !setup to rotate ylmc_fac to new axes if called for 3359 !skip if only s projectors are used 3360 if ( lmax>0 ) then 3361 ! generate a set of nr=lmax2 random vectors 3362 ! idum=123456 3363 do ir=1,lmax2 3364 do ii=1,3 3365 r(ii,ir) = uniformrandom(idum)-0.5d0 3366 end do !ii 3367 call ylm_cmplx(lmax,ylmcp,r(1,ir),r(2,ir),r(3,ir)) 3368 ylmc_rr(ir,:)=conjg(ylmcp(:)) 3369 ylmc_rr_save(ir,:)=conjg(ylmcp(:)) 3370 end do !ir 3371 3372 ylmc_rrinv(:,:)=c0 3373 do ii=1,lmax2 3374 ylmc_rrinv(ii,ii)=c1 3375 end do !ii 3376 ! calculate inverse of ylmc(ir,lm) matrix 3377 call ZGESV(lmax2,lmax2,ylmc_rr,lmax2,ipiv,ylmc_rrinv,lmax2,info) 3378 3379 ! check that r points are independent (ie., that matrix inversion wasn't 3380 ! too close to singular) 3381 ylmc_rr=matmul(ylmc_rrinv,ylmc_rr_save) 3382 test=zero 3383 do ii=1,lmax2 3384 ylmc_rr(ii,ii)=ylmc_rr(ii,ii)-c1 3385 do jj=1,lmax2 3386 test=max(abs(ylmc_rr(ii,jj)),test) 3387 end do !ii 3388 end do !jj 3389 if(test>tol8) then 3390 write(message, '(5a)' )& 3391 & ' matrix inversion error for wannier rotations',ch10,& 3392 & ' random vectors r(j,1:nr) are not all independent !! ',ch10,& 3393 & ' Action : re-seed uniformrandom or maybe just try again' 3394 ABI_ERROR(message) 3395 end if !test>tol8 3396 3397 ! end of the preliminaries, now to the rotations of the wannier orbitals 3398 do iwan=1,nwan 3399 ! don't bother for s orbitals 3400 if(proj_l(iwan)==0) cycle 3401 ! check for default axes and cycle if found 3402 if(proj_z(1,iwan)==zero .and. proj_z(2,iwan)==zero .and.& 3403 & proj_z(3,iwan)== one .and. proj_x(1,iwan)==one .and.& 3404 & proj_x(2,iwan)==zero .and. proj_x(3,iwan)==zero) cycle 3405 3406 ! get the u matrix that rotates the reference frame 3407 call rotmat(proj_x(:,iwan),proj_z(:,iwan),inversion_flag,umat) 3408 3409 ! find rotated r-vectors. Optional inversion 3410 ! operation is an extension of the wannier90 axis-setting options 3411 ! which only allow for proper axis rotations 3412 if(inversion_flag==1) then 3413 rp(:,:)= -matmul ( umat(:,:), r(:,:) ) 3414 else 3415 rp(:,:) = matmul ( umat(:,:) , r(:,:) ) 3416 end if !inversion_flag 3417 3418 do ir=1,lmax2 3419 ! get the ylm representation of the rotated vectors 3420 call ylm_cmplx(lmax,ylmcp,rp(1,ir),rp(2,ir),rp(3,ir)) 3421 ylmc_rp(ir,:)=conjg(ylmcp(:)) 3422 end do !ir 3423 ! the matrix product sum(ir) ylmc_rrinv(lm,ir)*ylmc_rp(ir,lm') gives the 3424 ! the complex lmXlm matrix representation of the coordinate rotation 3425 crot(:,:)=matmul(ylmc_rrinv(:,:),ylmc_rp(:,:)) 3426 3427 ! now rotate the current wannier orbital 3428 ylmcp(:)=matmul(crot(:,:),ylmc_fac(:,iwan)) 3429 ylmc_fac(:,iwan)=ylmcp(:) 3430 3431 ! write(std_out,*)'ylmc_fac',ylmc_fac(:,iwan) 3432 end do !iwan 3433 end if !lmax>0 3434 3435 end subroutine mlwfovlp_ylmfac
m_mlwfovlp/mlwfovlp_ylmfar [ Functions ]
[ Top ] [ m_mlwfovlp ] [ Functions ]
NAME
mlwfovlp_ylmfar
FUNCTION
Routine that produces a fator by which the initial guess of functions will be multiplied for the Wannier90 interface. It is just used if there are rotations, or if the functions required are linear combinations of the ylm real functions. Example, For a function G(r)= 1/2 s + 1/3 px - 1/2 pz it would produce a matrix of the following form: [1/2,-1/2,1/3,0,0...0] This function is similar to mlwfovlp_ylmfac, but the factors it uses real spherical harmonics instead of complex spherical harmonics. Remember that real spherical harmonics are linear combinations of complex spherical harmonics
INPUTS
lmax= maximum l value for spherical harmonics lmax2=number of ylm functions mband=maximum number of bands nwan = number of wannier functions 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_x(3,mband)= x axis for the projection. proj_z(3,mband)= z axis for the projection.
OUTPUT
ylmc_fac(lmax2,nwan)=matrix containig a factor for ylm hybrid orbitals
SIDE EFFECTS
NOTES
SOURCE
3478 subroutine mlwfovlp_ylmfar(ylmr_fac,lmax,lmax2,mband,nwan,proj_l,proj_m,proj_x,proj_z) 3479 3480 !Arguments ------------------------------------ 3481 integer, intent(in):: lmax,lmax2,nwan,mband 3482 ! arrays 3483 integer,intent(in) :: proj_l(mband),proj_m(mband) 3484 real(dp),intent(in) :: proj_x(3,mband),proj_z(3,mband) 3485 real(dp),intent(out)::ylmr_fac(lmax2,nwan) 3486 ! 3487 !Local variables------------------------------- 3488 ! 3489 integer :: idum,ii,inversion_flag 3490 integer :: ir,iwan,jj,ll,lm,mm,mr 3491 real(dp):: onem,test 3492 ! arrays 3493 real(dp),allocatable::dummy(:,:),nrm(:) 3494 real(dp)::r(3,lmax2),rp(3,lmax2) 3495 real(dp)::rs2,rs3,rs6,rs12,umat(3,3) 3496 real(dp)::rot(lmax2,lmax2),tor(lmax2,lmax2),orb_lm(lmax2,-5:3,7) 3497 real(dp):: ylmrp(lmax2) 3498 real(dp):: ylmr_rr(lmax2,lmax2),ylmr_rr_save(lmax2,lmax2) 3499 real(dp):: ylmr_rrinv(lmax2,lmax2),ylmr_rp(lmax2,lmax2) 3500 character(len=500) :: message ! to be uncommented, if needed 3501 !no_abirules 3502 !integer :: orb_idx(16)=(/1,3,4,2,7,8,6,9,5,13,14,12,15,11,16,10/) !Tab3.1 Wannier90 user guide 3503 3504 ! ************************************************************************* 3505 3506 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3507 !DEBUG 3508 !write(std_out,*)'lmax ',lmax,'lmax2 ',lmax2 3509 !write(std_out,*)'mband ',mband,'nwan ',nwan 3510 ! 3511 !do iwan=1,nwan 3512 !write(std_out,*)'iwan,proj_l, proj_m',proj_l(iwan),proj_m(iwan) 3513 !write(std_out,*)'iwan,proj_x, proj_z',iwan,proj_x(:,iwan),proj_z(:,iwan) 3514 !end do 3515 !!END DEBUG 3516 3517 !constants for linear combinations of ylm's 3518 rs2=1._dp/sqrt(2._dp) 3519 rs3=1._dp/sqrt(3._dp) 3520 rs6=1._dp/sqrt(6._dp) 3521 rs12=1._dp/sqrt(12._dp) 3522 3523 ! 3524 !mapping lm coefficients for real spherical harmonics 3525 !table 3.1 of Wannier90 user guide with real spherical harmonics in routine initylmr 3526 !s, py,pz,px, dxy,dyz,dz2,dxz,dx2-y2, fy(3x2-y2),fxyz,fyz2,fz3,fxz2, 3527 !fz(x2-y2),fx(x2-3y2) 3528 !note: check ordering of f orbitals, it might be wrong 3529 3530 tor(:,:)=0.d0 3531 lm=0 3532 do ll=0,lmax 3533 do mm=-ll,ll 3534 onem=(-1.d0)**mm 3535 lm=lm+1 3536 if(ll == 0) then 3537 tor(lm,lm)=1.d0 3538 else 3539 tor(lm,lm)=onem*1.d0 3540 end if 3541 end do !mm 3542 end do !ll 3543 !do lm=1,16 3544 !write(std_out,*)'tor lm=',lm,tor(:,lm) 3545 !end do 3546 3547 !coefficients for basic wannier orbitals in Table 3.1 order 3548 orb_lm(:,:,:)=0.d0 3549 ii=0 3550 do ll=0,lmax 3551 do mr=1,2*ll+1 3552 ii=ii+1 3553 orb_lm(:,ll,mr)= tor(:,ii) 3554 ! write(std_out,*)'ii',ii,'orb_lm',orb_lm(:,ll,mr) 3555 end do 3556 end do 3557 3558 3559 3560 !coefficients for linear combinations in table 3.2 order 3561 if(lmax>=1) then 3562 ! s px 3563 orb_lm(:,-1,1)=rs2*tor(:,1)+rs2*tor(:,4) 3564 orb_lm(:,-1,2)=rs2*tor(:,1)-rs2*tor(:,4) 3565 ! s px py 3566 orb_lm(:,-2,1)=rs3*tor(:,1)-rs6*tor(:,4)+rs2*tor(:,2) 3567 orb_lm(:,-2,2)=rs3*tor(:,1)-rs6*tor(:,4)-rs2*tor(:,2) 3568 orb_lm(:,-2,3)=rs3*tor(:,1)+2._dp*rs6*tor(:,4) 3569 ! s px py pz 3570 orb_lm(:,-3,1)=half*(tor(:,1)+tor(:,4)+tor(:,2)+tor(:,3)) 3571 orb_lm(:,-3,2)=half*(tor(:,1)+tor(:,4)-tor(:,2)-tor(:,3)) 3572 orb_lm(:,-3,3)=half*(tor(:,1)-tor(:,4)+tor(:,2)-tor(:,3)) 3573 orb_lm(:,-3,4)=half*(tor(:,1)-tor(:,4)-tor(:,2)+tor(:,3)) 3574 end if 3575 if(lmax>=2) then 3576 ! s px py 3577 orb_lm(:,-4,1)=rs3*tor(:,1)-rs6*tor(:,4)+rs2*tor(:,2) 3578 orb_lm(:,-4,2)=rs3*tor(:,1)-rs6*tor(:,4)-rs2*tor(:,2) 3579 orb_lm(:,-4,3)=rs3*tor(:,1)+2._dp*rs6*tor(:,4) 3580 ! pz dz2 3581 orb_lm(:,-4,4)= rs2*tor(:,3)+rs2*tor(:,7) 3582 orb_lm(:,-4,5)=-rs2*tor(:,3)+rs2*tor(:,7) 3583 ! s px dz2 dx2-y2 3584 orb_lm(:,-5,1)=rs6*tor(:,1)-rs2*tor(:,4)-rs12*tor(:,7)+half*tor(:,9) 3585 orb_lm(:,-5,2)=rs6*tor(:,1)+rs2*tor(:,4)-rs12*tor(:,7)+half*tor(:,9) 3586 ! s py dz2 dx2-y2 3587 orb_lm(:,-5,3)=rs6*tor(:,1)-rs2*tor(:,2)-rs12*tor(:,7)-half*tor(:,9) 3588 orb_lm(:,-5,4)=rs6*tor(:,1)+rs2*tor(:,2)-rs12*tor(:,7)-half*tor(:,9) 3589 ! s pz dz2 3590 orb_lm(:,-5,5)=rs6*tor(:,1)-rs2*tor(:,3)+rs3*tor(:,7) 3591 orb_lm(:,-5,6)=rs6*tor(:,1)+rs2*tor(:,3)+rs3*tor(:,7) 3592 end if 3593 3594 !real wannier orbital coefficient array 3595 do iwan=1,nwan 3596 ylmr_fac(:,iwan)=orb_lm(:,proj_l(iwan),proj_m(iwan)) 3597 end do 3598 3599 3600 !setup to rotate ylmr_fac to new axes if called for 3601 !skip if only s projetors are used 3602 if ( lmax>0 ) then 3603 ! generate a set of nr=lmax2 random vetors 3604 idum=123456 3605 do ir=1,lmax2 3606 do ii=1,3 3607 r(ii,ir) = uniformrandom(idum)-0.5d0 3608 end do !ii 3609 end do !ir 3610 ABI_MALLOC(nrm,(lmax2)) 3611 nrm(:)=sqrt(r(1,:)**2+r(2,:)**2+r(3,:)**2)**0.5 3612 call initylmr(lmax+1,1,lmax2,nrm,1,r(:,:),ylmr_rr_save(:,:),dummy) 3613 ylmr_rr(:,:)=ylmr_rr_save(:,:) 3614 do ir=1,lmax2 3615 ylmr_rr_save(ir,:)=ylmr_rr(:,ir) 3616 end do 3617 ABI_FREE(nrm) 3618 3619 ylmr_rrinv(:,:)=0.d0 3620 do ii=1,lmax2 3621 ylmr_rrinv(ii,ii)=1.d0 3622 end do !ii 3623 ! calculate inverse of ylmr(ir,lm) matrix 3624 ylmr_rrinv(:,:)=ylmr_rr_save(:,:) 3625 call matrginv(ylmr_rrinv,lmax2,lmax2) 3626 3627 ! check that r points are independent (ie., that matrix inversion wasn't 3628 ! too close to singular) 3629 ylmr_rr=matmul(ylmr_rrinv,ylmr_rr_save) 3630 test=0.d0 3631 do ii=1,lmax2 3632 ylmr_rr(ii,ii)=ylmr_rr(ii,ii)-1.d0 3633 do jj=1,lmax2 3634 test=max(abs(ylmr_rr(ii,jj)),test) 3635 end do !ii 3636 end do !jj 3637 if(test>tol8) then 3638 write(message, '(5a)' )& 3639 & ' matrix inversion error for wannier rotations',ch10,& 3640 & ' random vetors r(j,1:nr) are not all independent !! ',ch10,& 3641 & ' Action : re-seed uniformrandom or maybe just try again' 3642 ABI_ERROR(message) 3643 end if !test>tol8 3644 3645 ! end of the preliminaries, now to the rotations of the wannier orbitals 3646 do iwan=1,nwan 3647 ! don't bother for s orbitals 3648 if(proj_l(iwan)==0) cycle 3649 ! check for default axes and cycle if found 3650 if(proj_z(1,iwan)==0.d0 .and. proj_z(2,iwan)==0.d0 .and.& 3651 & proj_z(3,iwan)== 1.d0 .and. proj_x(1,iwan)==1.d0 .and.& 3652 & proj_x(2,iwan)==0.d0 .and. proj_x(3,iwan)==0.d0) cycle 3653 3654 ! get the u matrix that rotates the reference frame 3655 call rotmat(proj_x(:,iwan),proj_z(:,iwan),inversion_flag,umat) 3656 ! 3657 ! find rotated r-vetors. Optional inversion 3658 ! operation is an extension of the wannier90 axis-setting options 3659 ! which only allow for proper axis rotations 3660 if(inversion_flag==1) then 3661 rp(:,:)= -matmul ( umat(:,:), r(:,:) ) 3662 else 3663 rp(:,:) = matmul ( umat(:,:) , r(:,:) ) 3664 end if !inversion_flag 3665 3666 ! get the ylm representation of the rotated vetors 3667 ABI_MALLOC(nrm,(lmax2)) 3668 nrm(:)=sqrt(rp(1,:)**2+rp(2,:)**2+rp(3,:)**2)**0.5 3669 call initylmr(lmax+1,1,lmax2,nrm,1,rp(:,:),ylmr_rp(:,:),dummy) 3670 ylmr_rr(:,:)=ylmr_rp(:,:) 3671 do ir=1,lmax2 3672 ylmr_rp(ir,:)=ylmr_rr(:,ir) 3673 end do 3674 ABI_FREE(nrm) 3675 ! the matrix product sum(ir) ylmr_rrinv(lm,ir)*ylmr_rp(ir,lm') gives the 3676 ! the lmXlm matrix representation of the coordinate rotation 3677 3678 rot(:,:)=matmul(ylmr_rrinv(:,:),ylmr_rp(:,:)) 3679 ! 3680 ! now rotate the current wannier orbital 3681 ylmrp(:)=matmul(rot(:,:),ylmr_fac(:,iwan)) 3682 ylmr_fac(:,iwan)=ylmrp(:) 3683 end do !iwan 3684 end if !lmax>0 3685 3686 end subroutine mlwfovlp_ylmfar
mlwfovlp/read_chkunit [ Functions ]
[ Top ] [ mlwfovlp ] [ Functions ]
NAME
read_chkunit
FUNCTION
Function which reads the .chk file produced by Wannier90
INPUTS
OUTPUT
SOURCE
1336 subroutine read_chkunit(seed_name,nkpt,ndimwin,ierr) 1337 1338 !Arguments ------------------------------------ 1339 !scalars 1340 integer,intent(in) :: nkpt 1341 character(len=*),intent(in) :: seed_name 1342 integer,intent(out) :: ierr 1343 !arrays 1344 integer,intent(out) :: ndimwin(nkpt) 1345 1346 !Local variables------------------------------- 1347 !scalars 1348 integer :: chk_unit,ios,ikpt 1349 logical :: have_disentangled 1350 1351 !************************************************************************ 1352 1353 chk_unit=get_unit() 1354 fname=TRIM(seed_name)//'.chk' 1355 open(unit=chk_unit,file=fname,form='unformatted',status='old',iostat=ios) 1356 1357 ierr=0 1358 read(chk_unit) ! header ! Date and time 1359 read(chk_unit) ! ((real_lattice(i,j),i=1,3),j=1,3) ! Real lattice 1360 read(chk_unit) ! ((recip_lattice(i,j),i=1,3),j=1,3) ! Reciprocal lattice 1361 read(chk_unit) ! num_kpts 1362 read(chk_unit) ! ((kpt_latt(i,nkp),i=1,3),nkp=1,num_kpts) ! K-points 1363 read(chk_unit) ! nntot ! Number of nearest k-point neighbours 1364 read(chk_unit) ! num_wann ! Number of wannier functions 1365 read(chk_unit) ! chkpt1 ! Position of checkpoint 1366 read(chk_unit) have_disentangled ! Whether a disentanglement has been performed 1367 if (have_disentangled) then 1368 ! read(chk_unit) ! omega_invariant ! Omega invariant 1369 ! read(chk_unit) ((lwindow(i,nkp),i=1,num_bands),nkp=1,num_kpts) 1370 read(chk_unit) (ndimwin(ikpt),ikpt=1,nkpt) 1371 ! read(chk_unit) (((u_matrix_opt(i,j,nkp),i=1,num_bands),j=1,num_wann),nkp=1,num_kpts) 1372 else 1373 ! this is not expected. we should have disentanglement. Report the error. 1374 ierr=-1 1375 end if 1376 ! read(chk_unit) (((u_matrix(i,j,k),i=1,num_wann),j=1,num_wann),k=1,num_kpts) ! U_matrix 1377 ! read(chk_unit) ((((m_matrix(i,j,k,l,1),i=1,num_wann),j=1,num_wann),k=1,nntot),l=1,num_kpts) ! M_matrix 1378 ! read(chk_unit) ((wannier_centres(i,j),i=1,3),j=1,num_wann) 1379 close(chk_unit) 1380 1381 end subroutine read_chkunit