TABLE OF CONTENTS


ABINIT/m_mlwfovlp [ Modules ]

[ Top ] [ Modules ]

NAME

  m_mlwfovlp

FUNCTION

  Interface with Wannier90

COPYRIGHT

  Copyright (C) 2005-2018 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 .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_mlwfovlp
28 
29  use defs_basis
30  use defs_datatypes
31  use defs_abitypes
32  use defs_wannier90
33  use m_abicore
34  use m_errors
35  use m_atomdata
36  use m_xmpi
37  use m_sort
38 #ifdef FC_NAG
39  use f90_unix_dir
40 #endif
41 #ifdef HAVE_NETCDF
42  use netcdf
43 #endif
44  use m_nctk
45  use m_hdr
46 
47  use m_io_tools, only : delete_file, get_unit, open_file
48  use m_hide_lapack,     only : matrginv
49  use m_fstrings,      only : strcat, sjoin
50  use m_numeric_tools, only : uniformrandom, simpson_int, c2r, l2int
51  use m_special_funcs,   only : besjm
52  use m_geometry,  only : xred2xcart, rotmat, wigner_seitz
53  use m_fftcore,  only : sphereboundary
54  use m_crystal,  only : crystal_t
55  use m_crystal_io, only : crystal_ncwrite
56  use m_ebands,   only : ebands_ncwrite
57  use m_pawang,   only : pawang_type
58  use m_pawrad,   only : pawrad_type, simp_gen
59  use m_pawtab,   only : pawtab_type
60  use m_pawcprj,  only : pawcprj_type
61  use m_paw_sphharm, only : ylm_cmplx, initylmr
62  use m_paw_overlap, only : smatrix_pawinit
63  use m_evdw_wannier, only : evdw_wannier
64  use m_fft,            only : fourwf
65 
66  implicit none
67 
68  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

PARENTS

      outscfcv

CHILDREN

SOURCE

 138  subroutine mlwfovlp(crystal, ebands, hdr, atindx1,cg,cprj,dtset,dtfil,eigen,gprimd,kg,&
 139 & mband,mcg,mcprj,mgfftc,mkmem,mpi_enreg,mpw,natom,&
 140 & nattyp,nfft,ngfft,nkpt,npwarr,nsppol,ntypat,occ,&
 141 & pawang,pawrad,pawtab,prtvol,psps,rprimd,ucvol,xred)
 142 
 143 
 144 !This section has been created automatically by the script Abilint (TD).
 145 !Do not modify the following lines by hand.
 146 #undef ABI_FUNC
 147 #define ABI_FUNC 'mlwfovlp'
 148 !End of the abilint section
 149 
 150  implicit none
 151 
 152 !Arguments ------------------------------------
 153 !scalars
 154  integer,intent(in) :: mband,mcg,mcprj,mgfftc,mkmem,mpw,natom,nfft,nkpt
 155  integer,intent(in) :: nsppol,ntypat,prtvol
 156  real(dp),intent(in) :: ucvol
 157  type(crystal_t),intent(in) :: crystal
 158  type(ebands_t),intent(in) :: ebands
 159  type(hdr_type),intent(in) :: hdr
 160  type(MPI_type),intent(in) :: mpi_enreg
 161  type(dataset_type),intent(in) :: dtset
 162  type(datafiles_type),intent(in) :: dtfil
 163  type(pawang_type),intent(in) :: pawang
 164  type(pseudopotential_type),intent(in) :: psps
 165 !arrays
 166  integer,intent(in) :: atindx1(natom)
 167  integer :: kg(3,mpw*mkmem),nattyp(ntypat),ngfft(18),npwarr(nkpt)
 168  real(dp),intent(in) :: cg(2,mcg)
 169  real(dp),intent(in) :: eigen(mband*nkpt*nsppol),gprimd(3,3),rprimd(3,3)
 170  real(dp),intent(in) :: occ(mband*nkpt*nsppol)
 171  real(dp),intent(in) :: xred(3,natom)
 172  type(pawcprj_type) :: cprj(natom,mcprj)
 173  type(pawrad_type),intent(in) :: pawrad(psps%ntypat*psps%usepaw)
 174  type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
 175 
 176 !Local variables-------------------------------
 177 !scalars
 178  integer :: band_index,cplex,i,iatom,iband,iband1,iband2,icgtemp
 179  integer :: ig,ii,ikg,ierr
 180  integer :: ikpt,ikpt1,ikpt2,ilmn,intot,isppol,itypat
 181  integer :: iun(nsppol),iun_plot,iwan,jband,jband1,jband2,jj,jj1,jj2,jj3
 182  integer :: lmn_size,lproj,lwanniersetup,mwan,mgfft,n1
 183 #if defined HAVE_WANNIER90
 184  integer :: kk
 185 #endif
 186 #ifdef HAVE_NETCDF
 187  integer :: ncid, ncerr, nrpts
 188  character(len=fnlen) :: abiwan_fname
 189  integer :: have_disentangled_spin(nsppol)
 190  integer,allocatable :: irvec(:,:),ndegen(:)
 191 #endif
 192  integer :: n1tmp,n2,n2tmp,n3,n3tmp,n4,n5,n6,nband_k
 193  integer :: nntot,npw_k,num_nnmax,spacing
 194  integer :: tim_fourwf
 195  integer :: master,max_num_bands,nprocs,spaceComm,spin,rank
 196 !integer :: j,k,l
 197  integer  :: nwan(nsppol),nband_inc(nsppol),num_bands(nsppol)
 198  real(dp) :: weight
 199 #if defined HAVE_WANNIER90
 200  real(dp) :: corrvdw
 201  complex(dpc) :: caux,caux2,caux3
 202 #endif
 203  logical :: gamma_only,leig,lmmn,lwannierrun,spinors !,have_disentangled
 204  character(len=fnlen) :: wfnname
 205  character(len=1000) :: message
 206  character(len=fnlen) :: seed_name(nsppol)
 207  character(len=fnlen) :: fname,filew90_win(nsppol),filew90_wout(nsppol),filew90_amn(nsppol),filew90_ramn(nsppol)
 208  character(len=fnlen) :: filew90_mmn(nsppol),filew90_eig(nsppol)
 209 !arrays
 210  integer :: g1temp(3),ngkpt(3)
 211  integer,allocatable :: g1(:,:,:),gbound(:,:),icg(:,:)
 212  integer,allocatable:: iwav(:,:,:),kg_k(:,:),ovikp(:,:)
 213  integer,allocatable :: proj_l(:,:),proj_m(:,:),proj_radial(:,:)
 214  integer,allocatable :: proj_s_loc(:)
 215  real(dp) :: real_lattice(3,3)
 216  real(dp) :: recip_lattice(3,3)
 217  real(dp),allocatable :: cm1(:,:,:,:,:,:),cm2_paw(:,:,:),cwavef(:,:)
 218  real(dp),allocatable :: denpot(:,:,:)
 219  real(dp),allocatable :: eigenvalues_w(:,:,:),fofgout(:,:),fofr(:,:,:,:)
 220  real(dp),allocatable :: proj_site(:,:,:),proj_x(:,:,:),proj_z(:,:,:),proj_zona(:,:)
 221  real(dp),allocatable :: wann_centres(:,:,:),wann_spreads(:,:),xcart(:,:)
 222  real(dp),allocatable :: proj_s_qaxis_loc(:,:)
 223  complex(dpc),allocatable :: A_paw(:,:,:,:)
 224  complex(dpc),allocatable :: M_matrix(:,:,:,:,:),U_matrix(:,:,:,:)
 225  complex(dpc),allocatable :: U_matrix_opt(:,:,:,:)
 226  complex(dpc),pointer :: A_matrix(:,:,:,:)
 227  logical,allocatable :: band_in(:,:),lwindow(:,:,:)
 228  character(len=3),allocatable :: atom_symbols(:)
 229  logical,allocatable::just_augmentation(:,:)
 230 #if defined HAVE_WANNIER90
 231  real(dp) :: spreadw(3,nsppol)
 232  real(dp),allocatable :: csix(:,:,:,:)
 233  real(dpc),allocatable :: occ_arr(:,:,:),occ_wan(:,:,:)
 234  real(dp),allocatable :: tdocc_wan(:,:)
 235 #endif
 236 
 237 !************************************************************************
 238 
 239 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 240 !1) Initialize variables and allocations
 241 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 242 !
 243 !Some initialization and checks
 244 !
 245  lwanniersetup=1 ! 1 is mandatory ( 0 is for debug)
 246 !to use lwanniersetup=0, one would need
 247 !to define which bands to exclude.
 248  lwannierrun=.true.   ! .false. and .true. are possible
 249  lmmn=.true.          ! .false. and .true. are possible
 250  leig=.true.          ! .false. and .true. are possible
 251 !
 252  gamma_only=.false. !not yet implemented
 253  spinors=.false. !not yet implemented
 254  spin=dtset%useria !spin polarization, just used in case of nsppol==2
 255 !0=> both (default)
 256 !1=> spin up
 257 !2=> spin down
 258 !
 259 !mpi initialization
 260 !
 261  spaceComm=MPI_enreg%comm_cell
 262  nprocs=xmpi_comm_size(spaceComm)
 263  rank=MPI_enreg%me_kpt
 264  master=0
 265 
 266 !write(std_out,'("master ",i3," rank ",i3," nprocs ",i3)') master,rank,nprocs
 267 !
 268 !Generate seed names for wannier90 files, and file names
 269 !
 270  call mlwfovlp_seedname(dtfil%fnameabo_w90,filew90_win,filew90_wout,filew90_amn,&
 271 & filew90_ramn,filew90_mmn,filew90_eig,nsppol,seed_name)
 272 !
 273 !Check the validity of input variables
 274 !FIXME: this is not a check, and prints a warning even if the input is fine!
 275 !must be changed to not print anything if kptopt 3 and istwfk 1 (the latter is easier to check)
 276 !
 277  if(rank==master) then
 278    write(message, '(a,a,a,a)' ) ch10,&
 279 &   '   mlwfovlp:  you should give k-point in the full brillouin zone ',ch10,&
 280 &   '   with explicit k-points (or kptopt=3) and istwfk 1'
 281    call wrtout(ab_out,message,'COLL')
 282    call wrtout(std_out,  message,'COLL')
 283  end if
 284 !
 285  if(MPI_enreg%paral_spinor==1) then
 286    message = ' Parallelization over spinorial components not yet available !'
 287    MSG_ERROR(message)
 288  end if
 289 
 290  if(nsppol==2) then
 291    if(spin==1) then
 292      write(message, '(3a)' ) ch10,&
 293 &     '   mlwfovlp:  Calculating matrices for spin-up polarization  ',ch10
 294    elseif(spin==2) then
 295      write(message, '(3a)' ) ch10,&
 296 &     '   mlwfovlp:  Calculating matrices for spin-down polarization  ',ch10
 297    end if
 298    call wrtout(ab_out,message,'COLL')
 299    call wrtout(std_out,  message,'COLL')
 300    if(spin <0 .or. spin>2) then
 301      message = '  mlwfovlp:  spin variable should be equal to 0, 1 or 2 '
 302      MSG_ERROR(message)
 303    end if
 304  end if
 305 !
 306 !get lattice parameters in wannier90 format
 307 !
 308  do i=1, 3
 309    real_lattice(:,i)=Bohr_Ang*rprimd(i,:)
 310    recip_lattice(:,i)=two_pi*gprimd(i,:)/Bohr_Ang
 311  end do
 312 !
 313  if(psps%npsp/=psps%ntypat) then
 314    MSG_ERROR("prb npsp")
 315  end if
 316 !
 317 !Allocations.
 318 !
 319  num_nnmax=12 !limit fixed for compact structure in wannier_setup.
 320  ABI_ALLOCATE(g1,(3,nkpt,num_nnmax))
 321  ABI_ALLOCATE(ovikp,(nkpt,num_nnmax))
 322  ABI_ALLOCATE(atom_symbols,(natom))
 323  ABI_ALLOCATE(xcart,(3,natom))
 324  ABI_ALLOCATE(band_in,(mband,nsppol))
 325  ABI_ALLOCATE(proj_site,(3,mband,nsppol))
 326  ABI_ALLOCATE(proj_l,(mband,nsppol))
 327  ABI_ALLOCATE(proj_m,(mband,nsppol))
 328  ABI_ALLOCATE(proj_radial,(mband,nsppol))
 329  ABI_ALLOCATE(proj_x,(3,mband,nsppol))
 330  ABI_ALLOCATE(proj_s_loc,(mband))
 331  ABI_ALLOCATE(proj_s_qaxis_loc,(3,mband))
 332  ABI_ALLOCATE(proj_z,(3,mband,nsppol))
 333  ABI_ALLOCATE(proj_zona,(mband,nsppol))
 334 !
 335 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 336 !2) Call to  Wannier setup
 337 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 338 !
 339  nullify(A_matrix)
 340 
 341  call mlwfovlp_setup(atom_symbols,band_in,dtset,filew90_win,gamma_only,&
 342 & g1,lwanniersetup,mband,natom,nband_inc,nkpt,&
 343 & nntot,num_bands,num_nnmax,nsppol,nwan,ovikp,&
 344 & proj_l,proj_m,proj_radial,proj_site,proj_s_loc, proj_s_qaxis_loc, proj_x,proj_z,proj_zona,&
 345 & real_lattice,recip_lattice,rprimd,seed_name,spin,spinors,xcart,xred)
 346 
 347  do isppol=1, nsppol
 348    if(spin.ne.0 .and. spin.ne.isppol) cycle
 349    write(message, '(6a)' ) ch10,&
 350 &   '   mlwfovlp :  mlwfovlp_setup done -',ch10,&
 351 &   '-  see ',trim(filew90_wout(isppol)),' for details.'
 352    call wrtout(ab_out,message,'COLL')
 353    call wrtout(std_out,  message,'COLL')
 354  end do
 355 !
 356 !some allocations after wannier90 setup
 357 !
 358  max_num_bands=maxval(num_bands(:))
 359  mwan=maxval(nwan(:))
 360  ABI_ALLOCATE(eigenvalues_w,(max_num_bands,nkpt,nsppol))
 361  ABI_ALLOCATE(M_matrix,(max_num_bands,max_num_bands,nntot,nkpt,nsppol))
 362  ABI_ALLOCATE(A_matrix,(max_num_bands,mwan,nkpt,nsppol))
 363  ABI_ALLOCATE(iwav,(mband,nkpt,nsppol))
 364 
 365 !
 366 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 367 !3) Write Eigenvalues (file seed_name.eig)
 368 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 369 !
 370  if(leig) then
 371 !
 372 !  Assign file unit numbers
 373    if(rank==master) then
 374      iun(1)=444
 375      if(nsppol==2) iun(2)=455
 376      do isppol=1,nsppol
 377        if(spin.ne.0 .and. spin.ne.isppol) cycle
 378        open(unit=iun(isppol),file=filew90_eig(isppol),form='formatted',status='unknown')
 379      end do
 380    end if !rank==master
 381 !  Loop to write eigenvalues
 382    band_index=0
 383    do isppol=1,nsppol
 384      do ikpt=1,nkpt
 385        nband_k=dtset%nband(ikpt+(isppol-1)*nkpt)
 386        jband=0
 387        do iband=1,mband
 388          if(band_in(iband,isppol)) then
 389            jband=jband+1
 390            if(spin.eq.0 .or. spin.eq.isppol)then
 391 !            Writing data
 392              if(rank==master) write(iun(isppol), '(2i6,4x,f10.5)' ) jband,ikpt,Ha_eV*eigen(iband+band_index)
 393 !            Finish writing, now save eigenvalues
 394              eigenvalues_w(jband,ikpt,isppol)=Ha_eV*eigen(iband+band_index)
 395            end if !spin
 396          end if
 397        end do !iband
 398        band_index=band_index+nband_k
 399      end do !ikpt
 400    end do  !nsppol
 401    if(rank==master) then
 402      do isppol=1,nsppol
 403        if(spin.ne.0 .and. spin.ne.isppol) cycle
 404        close(iun(isppol))
 405      end do
 406      write(message, '(a,a)' ) ch10,&
 407 &     '   mlwfovlp :  eigenvalues written'
 408      call wrtout(std_out,  message,'COLL')
 409    end if !master
 410  end if !leig
 411 !else if( leig . and. lwannierun) then
 412 !read .eig file
 413 
 414 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 415 !4) Calculate overlaps (file seed_name.mmn)
 416 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 417 !
 418 !First calculate indices and shift
 419 !
 420 !write(std_out,*) "Computes shift for cg"
 421  write(message, '(a,a)' ) ch10,&
 422 & '   mlwfovlp : compute shifts for g-points '
 423  call wrtout(std_out,  message,'COLL')
 424 !----------------------------------------------------------------------
 425 !Compute shifts for g points (icg,iwav)
 426 !(here mband is not used, because shifts are internal variables of abinit)
 427 !----------------------------------------------------------------------
 428 !write(std_out,*) mpw*dtset%nspinor*mband*mkmem*nsppol
 429  ABI_ALLOCATE(icg,(nsppol,nkpt))
 430  icg=0
 431  icgtemp=0
 432  iwav(:,:,:)=0
 433  do isppol=1,nsppol
 434    do ikpt=1,nkpt
 435 !
 436 !    MPI:cycle over k-points not treated by this node
 437 !
 438 
 439      if (nprocs>1 ) then !sometimes we can have just one processor
 440        if ( ABS(MPI_enreg%proc_distrb(ikpt,1,isppol)-rank)  /=0) CYCLE
 441      end if
 442 
 443 !    write(std_out,*)'rank',rank,'ikpt',ikpt,'isppol',isppol
 444      nband_k=dtset%nband(ikpt+(isppol-1)*nkpt)
 445 !    write(std_out,*) ikpt+(isppol-1)*nkpt,nkpt
 446      npw_k=npwarr(ikpt)
 447      do iband=1,nband_k
 448        if(iband.gt.mband) then
 449          write(message,'(a,3i0)')" mband",iband,mband,nband_k
 450          MSG_ERROR(message)
 451        end if
 452        iwav(iband,ikpt,isppol)= &
 453 &       (iband-1)*npw_k*dtset%nspinor+icgtemp
 454      end do ! iband
 455      icgtemp=icgtemp+ npw_k*dtset%nspinor*nband_k
 456 !    icg(isppol,ikpt)=icgtemp
 457 !    write(std_out,*) "icg", isppol,ikpt,icg(isppol,ikpt)
 458    end do  ! ikpt
 459  end do   ! isppol
 460 !write(std_out,*) "shift for cg computed"
 461  ABI_DEALLOCATE(icg)
 462 !
 463 !Shifts computed.
 464 !
 465  if( lmmn) then
 466 !
 467 !  In case of parallelization write out cg for all k-points
 468 !
 469    if (nprocs > 1) then
 470 !
 471      if(prtvol>0) then
 472        write(message, '(3a)' ) ch10,&
 473 &       '   mlwfovlp :  Creating temporary files with cg and cprj (PAW)',ch10
 474        call wrtout(ab_out,message,'COLL')
 475        call wrtout(std_out,  message,'COLL')
 476      end if
 477 !
 478      do isppol=1,nsppol
 479        if(spin.ne.0 .and. spin.ne.isppol) cycle
 480        do ikpt=1,nkpt
 481 !
 482 !        MPI:cycle over k-points not treated by this node
 483 !
 484          if (nprocs>1 ) then !sometimes we can have just one processor
 485            if ( ABS(MPI_enreg%proc_distrb(ikpt,1,isppol)-rank)  /=0) CYCLE
 486          end if
 487 
 488 !        write(std_out,*)'writing kpt ',ikpt,'isppol',isppol,' by node ', rank
 489          write(wfnname,'(a,I5.5,".",I1)') trim(dtfil%fnametmp_cg),ikpt,isppol
 490          iun_plot=1000+ikpt+ikpt*(isppol-1)
 491 
 492          open (unit=iun_plot, file=wfnname,form='unformatted')
 493          npw_k=npwarr(ikpt)
 494          do iband=1,mband
 495            do ig=1,npw_k*dtset%nspinor
 496              write(iun_plot) (cg(i,ig+iwav(iband,ikpt,isppol)),i=1,2)
 497            end do
 498          end do
 499          close(iun_plot)
 500        end do !ikpt
 501      end do !isppol
 502 !
 503 !    In the PAW case we also need to write out cprj into files
 504 !
 505      if(psps%usepaw==1) then
 506 !
 507 !      big loop on atoms, kpts, bands and lmn
 508 !
 509        ikpt2=0
 510        do isppol=1,nsppol
 511          if(spin.ne.0 .and. spin.ne.isppol) cycle
 512          do ikpt=1,nkpt
 513 !
 514 !          MPI:cycle over k-points not treated by this node
 515 !
 516            if (nprocs>1 ) then !sometimes we can have just one processor
 517              if ( ABS(MPI_enreg%proc_distrb(ikpt,1,isppol)-MPI_enreg%me)  /=0) CYCLE
 518            end if
 519 
 520            ikpt2=ikpt2+1 !sums just on the k-points treated by this node
 521 !
 522            write(wfnname,'(a,I5.5,".",I1)') trim(dtfil%fnametmp_cprj),ikpt,isppol
 523            iun_plot=1000+ikpt
 524            open (unit=iun_plot, file=wfnname,form='unformatted')
 525 !
 526            do iband=1,mband*dtset%nspinor
 527              ig=iband+(ikpt2-1)*mband*dtset%nspinor +(isppol-1)*nkpt*mband*dtset%nspinor !index for cprj(:,ig)
 528 !
 529              do iatom=1,natom
 530                itypat=dtset%typat(iatom)
 531                lmn_size=pawtab(itypat)%lmn_size
 532 !
 533                do ilmn=1,lmn_size
 534                  write(iun_plot) (( cprj(iatom,ig)%cp(i,ilmn)),i=1,2)
 535                end do !ilmn
 536              end do !iatom
 537            end do !iband
 538 
 539            close(iun_plot)
 540          end do !ikpt
 541        end do !isppol
 542      end if !usepaw==1
 543 
 544 !
 545 !
 546    end if !MPI nprocs>1
 547 !
 548 !  End of MPI preliminarities
 549 !  Calculate PW contribution of overlaps
 550 !
 551    ABI_ALLOCATE(cm1,(2,mband,mband,nntot,nkpt,nsppol))
 552    call mlwfovlp_pw(cg,cm1,g1,iwav,kg,mband,&
 553 &   mkmem,mpi_enreg,mpw,nfft,ngfft,nkpt,nntot,&
 554 &   npwarr,dtset%nspinor,nsppol,ovikp,dtfil%fnametmp_cg,spin)
 555    write(message, '(a,a)' ) ch10,&
 556 &   '   mlwfovlp : PW part of overlap computed   '
 557    call wrtout(std_out,  message,'COLL')
 558 !
 559 !  compute PAW Contribution and add it to PW contribution
 560 !
 561    if(psps%usepaw==1) then
 562      write(message, '(a,a)' ) ch10,&
 563 &     '** smatrix_pawinit : PAW part of overlap  '
 564      call wrtout(std_out,  message,'COLL')
 565      ABI_ALLOCATE(cm2_paw,(2,mband,mband))
 566      do isppol=1,nsppol
 567        if(spin.ne.0 .and. spin.ne.isppol) cycle
 568        do ikpt1=1,nkpt
 569 !
 570 !        MPI:cycle over k-points not treated by this node
 571 !
 572          if (nprocs>1 ) then !sometimes we can have just one processor
 573            if ( ABS(MPI_enreg%proc_distrb(ikpt1,1,isppol)-rank)  /=0) CYCLE
 574          end if
 575 
 576          write(message, '(a,i6,a,2i6)' ) &
 577 &         '   processor',rank,' computes PAW part for kpt and spin',ikpt1,isppol
 578          call wrtout(std_out,  message,'COLL')
 579 
 580          do intot=1,nntot
 581            ikpt2= ovikp(ikpt1,intot)
 582            g1temp(:)=g1(:,ikpt1,intot)
 583            call smatrix_pawinit(atindx1,cm2_paw,cprj,ikpt1,ikpt2,isppol,&
 584 &           g1temp,gprimd,dtset%kpt,mband,mband,mkmem,mpi_enreg,&
 585 &           natom,dtset%nband,nkpt,dtset%nspinor,nsppol,dtset%ntypat,pawang,pawrad,pawtab,rprimd,&
 586 &           dtfil%fnametmp_cprj,dtset%typat,xred)
 587 !          cm1(:,:,:,intot,ikpt1,isppol)=four_pi*cm2_paw(:,:,:)
 588 !           write(6,*) "ikpt1=",ikpt1
 589 !           do iband=1,mband
 590 !             write(6,*) "iband=",iband
 591 !             write(6,*) "Wannier PW       overlap",cm1(:,iband,iband,intot,ikpt1,isppol)
 592 !             write(6,*) "Wannier PAW      overlap",four_pi*cm2_paw(:,iband,iband)
 593 !             write(6,*) "Wannier PW+PAW   overlap",cm1(:,iband,iband,intot,ikpt1,isppol)+four_pi*cm2_paw(:,iband,iband)
 594 !           enddo
 595            cm1(:,:,:,intot,ikpt1,isppol)=cm1(:,:,:,intot,ikpt1,isppol)+four_pi*cm2_paw(:,:,:)
 596          end do ! intot
 597        end do ! ikpt1
 598      end do ! isppol
 599      ABI_DEALLOCATE(cm2_paw)
 600      write(message, '(a,a)' ) ch10,&
 601 &     '   mlwfovlp : PAW part of overlap computed '
 602      call wrtout(std_out,  message,'COLL')
 603    end if ! usepaw
 604 !
 605    call xmpi_barrier(spaceComm)
 606    call xmpi_sum(cm1,spaceComm,ierr)
 607 !
 608 !  write overlap for separate calculation of wannier functions
 609 !
 610    if(rank==master) then
 611      do isppol=1,nsppol !we write separate output files for each isppol
 612        if(spin.ne.0 .and. spin.ne.isppol) cycle
 613        iun(isppol)=220+isppol
 614        open(unit=iun(isppol),file=filew90_mmn(isppol),form='formatted',status='unknown')
 615        write(iun(isppol),*) "nnkp version 90"
 616        write(iun(isppol),*) num_bands(isppol),nkpt,nntot
 617      end do
 618    end if ! rank==master
 619 
 620    do isppol=1,nsppol
 621      if(spin.ne.0 .and. spin.ne.isppol) cycle
 622      do ikpt1=1,nkpt
 623        do intot=1,nntot
 624          if( rank==master) write(iun(isppol),'(2i6,3x,3x,3i5)') ikpt1,ovikp(ikpt1,intot),(g1(jj,ikpt1,intot),jj=1,3)
 625          jband2=0
 626          do iband2=1,mband ! the first index is faster
 627            if(band_in(iband2,isppol)) then
 628              jband2=jband2+1
 629              jband1=0
 630              do iband1=1,mband
 631                if(band_in(iband1,isppol)) then
 632                  jband1=jband1+1
 633                  if(rank==master) write(iun(isppol),*) &
 634 &                 cm1(1,iband1,iband2,intot,ikpt1,isppol),cm1(2,iband1,iband2,intot,ikpt1,isppol)
 635                  M_matrix(jband1,jband2,intot,ikpt1,isppol)=&
 636 &                 cmplx(cm1(1,iband1,iband2,intot,ikpt1,isppol),cm1(2,iband1,iband2,intot,ikpt1,isppol))
 637 !                write(2211,*) ikpt1,intot,iband1,iband2
 638 !                write(2211,*) cm1(1,iband1,iband2,intot,ikpt1,isppol),cm1(2,iband1,iband2,intot,ikpt1,isppol)
 639                end if ! band_in(iband1)
 640              end do ! iband1
 641            end if ! band_in(iband2)
 642          end do ! iband2
 643        end do !intot
 644      end do !ikpt
 645      if( rank==master ) then
 646        close(iun(isppol))
 647        write(message, '(3a)' )  '   ',trim(filew90_mmn(isppol)),' written'
 648        call wrtout(std_out,  message,'COLL')
 649      end if !rank==master
 650    end do !isppol
 651 !
 652    ABI_DEALLOCATE(cm1)
 653 !
 654 !  Write down part of the matrix to the output file
 655 !  This is for the automatic tests
 656 !
 657    if(rank==master) then
 658      write(message, '(4a)' ) ch10,&
 659 &     '   Writing top of the overlap matrix: M_mn(ikb,ik)',ch10,&
 660 &     '   m=n=1:3, ikb=1, ik=1'
 661      call wrtout(ab_out,message,'COLL')
 662      call wrtout(std_out,  message,'COLL')
 663 !
 664 !    just write down the first 3 elements
 665 !
 666      do isppol=1,nsppol
 667        if(spin.ne.0 .and. spin.ne.isppol) cycle
 668        write(message, '( " " )')
 669        if (nsppol>1 ) then
 670          if (isppol==1) write(message,'(2a)')trim(message),'   spin up:'
 671          if (isppol==2) write(message,'(2a)')trim(message),'   spin down:'
 672        end if
 673        do ii=1,3
 674          if(ii>num_bands(isppol)) cycle
 675          write(message,'(3a)') trim(message),ch10,';   ( '
 676          do jj=1,3
 677            if(jj>num_bands(isppol))cycle
 678            write(message, '(a,2f11.6,a)') trim(message),&
 679 &           M_matrix(ii,jj,1,1,isppol),' , '
 680          end do
 681          write(message,'(2a)') trim(message),'    ) '
 682        end do
 683        call wrtout(ab_out,message,'COLL')
 684        call wrtout(std_out,  message,'COLL')
 685      end do
 686 !
 687 !    Now write down bottom of the matrix
 688 !
 689      write(message, '(4a)' ) ch10,&
 690 &     '   Writing bottom of the overlap matrix: M_mn(ikb,ik)',ch10,&
 691 &     '   m=n=num_bands-2:num_bands, ikb=nntot, ik=nkpt'
 692      call wrtout(ab_out,message,'COLL')
 693      call wrtout(std_out,  message,'COLL')
 694 !
 695      do isppol=1,nsppol
 696        if(spin.ne.0 .and. spin.ne.isppol) cycle
 697        write(message, '( " " )')
 698        if (nsppol>1 ) then
 699          if (isppol==1) write(message,'(2a)')trim(message),'   spin up:'
 700          if (isppol==2) write(message,'(2a)')trim(message),'   spin down:'
 701        end if
 702        do ii=num_bands(isppol)-2,num_bands(isppol)
 703          if(ii<1) cycle
 704          write(message,'(3a)') trim(message),ch10,';   ( '
 705          do jj=num_bands(isppol)-2,num_bands(isppol)
 706            if(jj<1)cycle
 707            write(message, '(a,2f11.6,a)') trim(message),&
 708 &           M_matrix(ii,jj,nntot,nkpt,isppol),' , '
 709          end do !j
 710          write(message,'(2a)') trim(message),'    ) '
 711        end do !ii
 712        call wrtout(ab_out,message,'COLL')
 713        call wrtout(std_out,  message,'COLL')
 714      end do !isppol
 715    end if !rank==master
 716 !
 717 !  erase temporary files created for parallel runs
 718 !
 719    if (nprocs > 1) then
 720 !
 721      if(prtvol>0) then
 722        write(message, '(3a)' ) ch10,&
 723 &       '   mlwfovlp :  Removing temporary files with cg and cprj (PAW)',ch10
 724        call wrtout(ab_out,message,'COLL')
 725        call wrtout(std_out,  message,'COLL')
 726      end if
 727 !
 728 !    Just master  node will remove the files
 729 !
 730      if(rank==master) then
 731        do isppol=1,nsppol
 732          if(spin.ne.0 .and. spin.ne.isppol) cycle
 733          do ikpt=1,nkpt
 734            write(wfnname,'(a,I5.5,".",I1)') trim(dtfil%fnametmp_cg),ikpt,isppol
 735            call delete_file(wfnname,ierr)
 736            if(psps%usepaw==1) then
 737              write(wfnname,'(a,I5.5,".",I1)') trim(dtfil%fnametmp_cprj),ikpt,isppol
 738              call delete_file(wfnname,ierr)
 739            end if
 740          end do !ikpt
 741        end do !isppol
 742      end if
 743    end if !MPI nprocs>1
 744 !
 745  end if !lmmn
 746 !if ( lmmn== .false. .and. lwannierun ) the
 747 !read .mmn file
 748 !
 749 !Deallocate arrays  no longer used
 750 !
 751  ABI_DEALLOCATE(ovikp)
 752  ABI_DEALLOCATE(g1)
 753 
 754 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 755 !5) Calculate initial projections
 756 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 757 
 758  if(dtset%w90iniprj/=0 )  then
 759 
 760 
 761 !
 762 !  Set value for lproj (type of projections to be computed)
 763 !  In PAW, options 5 and 6 are not in use.
 764 !  5 means that there will be a contribution from inside the spheres and another from the PW part
 765 !  6 means that we take into account just the inside-spheres contribution
 766 !  2 means that PW part will be calculated
 767 
 768 !
 769    lproj=dtset%w90iniprj
 770    if(dtset%w90iniprj == 5 ) lproj=2 ! Necessary to calculate PW contribution
 771 !
 772    ABI_ALLOCATE(just_augmentation,(mwan,nsppol))
 773    just_augmentation(:,:)=.false.
 774 
 775    if( psps%usepaw==1 .and. (dtset%w90iniprj==2 .or. dtset%w90iniprj>4)) then
 776      if (dtset%w90iniprj==6) just_augmentation(:,:)=.true.
 777      if (dtset%w90iniprj==5) then
 778        do isppol=1,nsppol
 779          if(spin.ne.0 .and. spin.ne.isppol) cycle
 780          do iwan=1,nwan(isppol)
 781 !
 782 !          Trick to skip the planewave contribution for some Wannier functions
 783 !          (Not in production).
 784 !
 785            if(proj_radial(iwan,isppol) > 4) then
 786              just_augmentation(iwan,isppol)=.true.
 787              proj_radial(iwan,isppol)=proj_radial(iwan,isppol)-3
 788              write(message, '(2a,2i4)' ) &
 789 &             '   ','Skiping planewave contribution for iwan, ispin=',iwan,isppol
 790              call wrtout(std_out,  message,'COLL')
 791            end if !proj_radial>4
 792          end do !iwan
 793        end do !isppol
 794      end if !w90iniprj == 5
 795    end if !paw
 796 !
 797 !  Call mlwfovlp_proj (plane waves part of projections)
 798 !
 799    if (dtset%w90iniprj/=6) then ! option 6 not yet in use
 800      call mlwfovlp_proj(A_matrix,band_in,cg,cprj,dtset,gprimd,just_augmentation,kg,&
 801 &     lproj,max_num_bands,mband,mkmem,mpi_enreg,mpw,mwan,natom,&
 802 &     nattyp,nkpt,npwarr,&
 803 &     dtset%nspinor,nsppol,ntypat,num_bands,nwan,pawtab,proj_l,proj_m,&
 804 &     proj_radial,proj_site,proj_x,proj_z,proj_zona,psps,spin,ucvol)
 805      write(message, '(a,a,a,a)' ) ch10,&
 806 &     '   mlwfovlp:  mlwfovlp_proj done -',ch10,&
 807 &     '   Projectors computed.'
 808      call wrtout(std_out,  message,'COLL')
 809    end if !w90proj/=6
 810 !
 811 !  Calculate inside-sphere part of projections (PAW)
 812 !
 813    if (psps%usepaw ==1 .and. ( dtset%w90iniprj>4)) then
 814      ABI_ALLOCATE(A_paw,(max_num_bands,mwan,nkpt,nsppol))
 815      call mlwfovlp_projpaw(A_paw,band_in,cprj,just_augmentation,max_num_bands,mband,mkmem,&
 816 &     mwan,natom,dtset%nband,nkpt,&
 817 &     dtset%nspinor,nsppol,dtset%ntypat,nwan,pawrad,pawtab,&
 818 &     proj_l,proj_m,proj_radial,proj_site,proj_x,proj_z,proj_zona,psps,&
 819 &     rprimd,spin,dtset%typat,xred)
 820 !
 821      write(message, '(a,a,a,a)' ) ch10,&
 822 &     '   mlwfovlp:  mlwfovlp_proj_paw done -',ch10,&
 823 &     '   Inside-spheres part of projectors computed.'
 824      call wrtout(std_out,  message,'COLL')
 825 !
 826 !    Add in-sphere contribution to A_matrix
 827 !
 828 !
 829 !    w90iniprj==5. Plane waves + augmentation contributions
 830 !
 831      if(dtset%w90iniprj==5) A_matrix(:,:,:,:)=A_matrix(:,:,:,:)+A_paw(:,:,:,:)
 832 !
 833 !    w90iniprj==6. Just augmentation contribution
 834 !
 835      if(dtset%w90iniprj==6) A_matrix(:,:,:,:)=A_paw(:,:,:,:)
 836 !
 837 !    deallocations
 838 !
 839      ABI_DEALLOCATE(A_paw)
 840    end if !usepaw==1
 841 
 842    ABI_DEALLOCATE(just_augmentation)
 843 !
 844    call xmpi_barrier(spaceComm)
 845    call xmpi_sum(A_matrix,spaceComm,ierr)
 846 
 847 !
 848 !  write      projections  to a file
 849 !
 850    if(rank==master) then
 851      if(dtset%w90iniprj==1) then
 852        do isppol=1,nsppol
 853          if(spin.ne.0 .and. spin.ne.isppol) cycle
 854          iun(isppol)=219+isppol
 855          open(unit=iun(isppol),file=trim(filew90_ramn(isppol)),form='formatted',status='unknown')
 856          write(iun(isppol),*) 'Projections from Abinit : mband,nkpt,nwan. indices: iband1,iwan,ikpt'
 857          write(iun(isppol),*) num_bands(isppol),nkpt,nwan(isppol)
 858        end do
 859      else
 860        do isppol=1,nsppol
 861          if(spin.ne.0 .and. spin.ne.isppol) cycle
 862          iun(isppol)=220+isppol
 863          open(unit=iun(isppol),file=trim(filew90_amn(isppol)),form='formatted',status='unknown')
 864          write(iun(isppol),*) 'Projections from Abinit : mband,nkpt,nwan. indices: iband1,iwan,ikpt'
 865          write(iun(isppol),*) num_bands(isppol),nkpt,nwan(isppol)
 866        end do
 867      end if
 868 !
 869      do isppol=1,nsppol
 870        if(spin.ne.0 .and. spin.ne.isppol) cycle
 871        do ikpt=1,nkpt
 872          do iwan=1,nwan(isppol)
 873            jband=0
 874            do iband=1,mband
 875              if(band_in(iband,isppol)) then
 876                jband=jband+1
 877                write(iun(isppol),'(3i6,13x,3x,2f18.14)')jband,iwan,ikpt,A_matrix(jband,iwan,ikpt,isppol)
 878              end if !band_in
 879            end do !iband
 880          end do !iwan
 881        end do !ikpt
 882      end do !isppol
 883 !
 884      if(dtset%w90iniprj==1) then
 885        do isppol=1,nsppol
 886          if(spin.ne.0 .and. spin.ne.isppol) cycle
 887          close(iun(isppol))
 888          write(message, '(3a)' ) &
 889 &         '   ',trim(filew90_ramn(isppol)),' written'
 890          call wrtout(std_out,  message,'COLL')
 891        end do
 892      else
 893        do isppol=1,nsppol
 894          if(spin.ne.0 .and. spin.ne.isppol) cycle
 895          close(iun(isppol))
 896          write(message, '(3a)' ) &
 897 &         '   ',trim(filew90_amn(isppol)),' written'
 898          call wrtout(std_out,  message,'COLL')
 899        end do
 900      end if
 901    end if !rank==master
 902 
 903 !
 904 !
 905 !  Write down part of the matrix to the output file
 906 !  This is for the automatic tests
 907 !
 908    if(rank==master) then
 909      write(message, '(4a)' ) ch10,&
 910 &     '   Writing top of the initial projections matrix: A_mn(ik)',ch10,&
 911 &     '   m=1:3, n=1:3, ik=1'
 912      call wrtout(ab_out,message,'COLL')
 913      call wrtout(std_out,  message,'COLL')
 914 !
 915 !    just write down the first 3 elements
 916 !
 917      do isppol=1,nsppol
 918        if(spin.ne.0 .and. spin.ne.isppol) cycle
 919        write(message, '( " " )')
 920        if (nsppol>1 ) then
 921          if (isppol==1) write(message,'(2a)')trim(message),'   spin up:'
 922          if (isppol==2) write(message,'(2a)')trim(message),'   spin down:'
 923        end if
 924        do ii=1,3
 925          if(ii>num_bands(isppol)) cycle
 926          write(message,'(3a)') trim(message),ch10,';   ( '
 927          do jj=1,3
 928            if(jj>nwan(isppol))cycle
 929            write(message, '(a,2f11.6,a)') trim(message),&
 930 &           A_matrix(ii,jj,1,isppol),' , '
 931          end do
 932          write(message,'(2a)') trim(message),'    ) '
 933        end do
 934        call wrtout(ab_out,message,'COLL')
 935        call wrtout(std_out,  message,'COLL')
 936      end do
 937 !
 938 !    Now write down bottom of the matrix
 939 !
 940      write(message, '(4a)' ) ch10,&
 941 &     '   Writing bottom of the initial projections matrix: A_mn(ik)',ch10,&
 942 &     '   m=num_bands-2:num_bands, n=nwan-2:nwan, ik=nkpt'
 943      call wrtout(ab_out,message,'COLL')
 944      call wrtout(std_out,  message,'COLL')
 945 !
 946      do isppol=1,nsppol
 947        if(spin.ne.0 .and. spin.ne.isppol) cycle
 948        write(message, '( " " )')
 949        if (nsppol>1 ) then
 950          if (isppol==1) write(message,'(2a)')trim(message),'   spin up:'
 951          if (isppol==2) write(message,'(2a)')trim(message),'   spin down:'
 952        end if
 953        do ii=num_bands(isppol)-2,num_bands(isppol)
 954          if(ii<1) cycle
 955          write(message,'(3a)') trim(message),ch10,';   ( '
 956          do jj=nwan(isppol)-2,nwan(isppol)
 957            if(jj<1)cycle
 958            write(message, '(a,2f11.6,a)') trim(message),&
 959 &           A_matrix(ii,jj,nkpt,isppol),' , '
 960          end do
 961          write(message,'(2a)') trim(message),'    ) '
 962        end do
 963        call wrtout(ab_out,message,'COLL')
 964        call wrtout(std_out,  message,'COLL')
 965      end do !isppol
 966    end if !rank==master
 967  end if !dtset%w90iniprj/=0
 968 !
 969 !Deallocations
 970 !
 971  ABI_DEALLOCATE(proj_site)
 972  ABI_DEALLOCATE(proj_l)
 973  ABI_DEALLOCATE(proj_m)
 974  ABI_DEALLOCATE(proj_radial)
 975  ABI_DEALLOCATE(proj_x)
 976  ABI_DEALLOCATE(proj_z)
 977  ABI_DEALLOCATE(proj_zona)
 978  ABI_DEALLOCATE(proj_s_loc)
 979  ABI_DEALLOCATE(proj_s_qaxis_loc)
 980 
 981 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 982 !6) write files for wannier function plot
 983 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 984  if( dtset%w90prtunk>0) then
 985    if(psps%usepaw==1) then
 986      write(message, '( a,a,a,a,a,a,a,a,a)')ch10,&
 987 &     "   WARNING: The UNK matrices will not contain the correct wavefunctions ",ch10,&
 988 &     "   since we are just writing the plane wave contribution.",ch10,&
 989 &     "   The contribution from inside the spheres is missing. ",ch10,&
 990 &     "   However, these files can be used for plotting purposes",ch10
 991      call wrtout(std_out,  message,'COLL')
 992    end if
 993 !
 994    spacing = dtset%w90prtunk
 995    write(message, '( 8a,i3,2a)')ch10,&
 996 &   "   UNK files will be written.",ch10,&
 997 &   "   According to the chosen value of w90prtunk",ch10,&
 998 &   "   the wavefunctions are to be written ",ch10, &
 999 &   "   at every ", spacing," records.",ch10
1000    call wrtout(std_out,  message,'COLL')
1001 !
1002    ABI_ALLOCATE(kg_k,(3,mpw))
1003    n1=ngfft(1)
1004    n2=ngfft(2)
1005    n3=ngfft(3)
1006    n4=ngfft(4)
1007    n5=ngfft(5)
1008    n6=ngfft(6)
1009    cplex=1
1010    mgfft=mgfftc ! error
1011    do isppol=1,nsppol
1012      ikg=0
1013      do ikpt=1,nkpt
1014 !
1015 !      MPI:cycle over k-points not treated by this node
1016 !
1017        if (nprocs>1 ) then !sometimes we can have just one processor
1018          if ( ABS(MPI_enreg%proc_distrb(ikpt,1,isppol)-rank)  /=0) CYCLE
1019        end if
1020 !
1021        if(spin.eq.0 .or. spin.eq.isppol) then
1022          npw_k=npwarr(ikpt)
1023          kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
1024          ABI_ALLOCATE(denpot,(cplex*n4,n5,n6))
1025          ABI_ALLOCATE(cwavef,(2,npw_k))
1026          ABI_ALLOCATE(fofr,(2,n4,n5,n6))
1027          ABI_ALLOCATE(gbound,(2*mgfft+8,2))
1028          ABI_ALLOCATE(fofgout,(2,npw_k))
1029          iun_plot=1000+ikpt+ikpt*(isppol-1)
1030          write(wfnname,'("UNK",I5.5,".",I1)') ikpt, isppol
1031 !        open (unit=iun_plot, file=wfnname,form='formatted')
1032          open(unit=iun_plot, file=wfnname,form='unformatted')
1033 !        optimizing grid for UNK files
1034          n1tmp = n1/spacing
1035          n2tmp = n2/spacing
1036          n3tmp = n3/spacing
1037          if( mod(n1,spacing) /= 0) then
1038            n1tmp = n1tmp + 1
1039          end if
1040          if( mod(n2,spacing) /= 0) then
1041            n2tmp = n2tmp + 1
1042          end if
1043          if( mod(n3,spacing) /= 0) then
1044            n3tmp = n3tmp + 1
1045          end if
1046 !        write(iun_plot,*) n1tmp,n2tmp,n3tmp,ikpt,nband_inc
1047          write(iun_plot) n1tmp,n2tmp,n3tmp,ikpt,nband_inc(isppol)
1048 !        gbound=zero
1049          call sphereboundary(gbound,dtset%istwfk(ikpt),kg_k,mgfft,npw_k)
1050          write(std_out,*) "  writes UNK file for ikpt, spin=",ikpt,isppol
1051          denpot(:,:,:)=zero
1052          weight = one
1053          do iband=1,mband
1054            if(band_in(iband,isppol)) then
1055              do ig=1,npw_k*dtset%nspinor
1056                cwavef(1,ig)=cg(1,ig+iwav(iband,ikpt,isppol))
1057                cwavef(2,ig)=cg(2,ig+iwav(iband,ikpt,isppol))
1058              end do
1059              tim_fourwf=0
1060              call fourwf(cplex,denpot,cwavef,fofgout,fofr,&
1061 &             gbound,gbound,dtset%istwfk(ikpt),kg_k,kg_k,mgfft,&
1062 &             mpi_enreg,1,ngfft,npw_k,npw_k,n4,n5,n6,0,dtset%paral_kgb,&
1063 &             tim_fourwf,weight,weight,use_gpu_cuda=dtset%use_gpu_cuda)
1064 !            do jj3=1,n3,spacing
1065 !            do jj2=1,n2,spacing
1066 !            do jj1=1,n1,spacing
1067 !            write(iun_plot,*) fofr(1,jj1,jj2,jj3),&
1068 !            & fofr(2,jj1,jj2,jj3)
1069 !            end do !jj1
1070 !            end do !jj2
1071 !            end do !jj3
1072 !            unformatted (must be one record)
1073              write(iun_plot) (((fofr(1,jj1,jj2,jj3),fofr(2,jj1,jj2,jj3),&
1074 &             jj1=1,n1,spacing),jj2=1,n2,spacing),jj3=1,n3,spacing)
1075            end if !iband
1076          end do ! iband
1077          ABI_DEALLOCATE(cwavef)
1078          ABI_DEALLOCATE(fofr)
1079          ABI_DEALLOCATE(gbound)
1080          ABI_DEALLOCATE(denpot)
1081          ABI_DEALLOCATE(fofgout)
1082        end  if !spin
1083        ikg=ikg+npw_k
1084        if(spin.eq.0 .or. spin.eq.isppol) close(iun_plot)
1085      end do  ! ikpt
1086    end do  ! nsppol
1087    ABI_DEALLOCATE(kg_k)
1088 !
1089    write(message, '(4a)' )ch10, &
1090 &   '   ','UNK files written',ch10
1091    call wrtout(std_out,  message,'COLL')
1092  end if !dtset%w90prtunk
1093 !
1094  ABI_DEALLOCATE(iwav)
1095 
1096 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1097 !7) Call to  Wannier90
1098 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1099  if(lwannierrun) then
1100    if(lwanniersetup.ne.1) MSG_ERROR("lwanniersetup.ne.1")
1101    ABI_ALLOCATE(U_matrix,(mwan,mwan,nkpt,nsppol))
1102    ABI_ALLOCATE(U_matrix_opt,(max_num_bands,mwan,nkpt,nsppol))
1103    ABI_ALLOCATE(lwindow,(max_num_bands,nkpt,nsppol))
1104    ABI_ALLOCATE(wann_centres,(3,mwan,nsppol))
1105    ABI_ALLOCATE(wann_spreads,(mwan,nsppol))
1106 !  Initialize
1107    U_matrix(:,:,:,:)=czero
1108    U_matrix_opt(:,:,:,:)=czero
1109    lwindow(:,:,:)=.false.
1110    wann_centres(:,:,:)=zero
1111    wann_spreads(:,:)=zero
1112 !
1113 !  write(std_out,*) seed_name
1114 !  write(std_out,*) ngkpt
1115    ngkpt(1)=dtset%kptrlatt(1,1)
1116    ngkpt(2)=dtset%kptrlatt(2,2) ! ajouter test de verif que kptrlatt est bien diagonal
1117    ngkpt(3)=dtset%kptrlatt(3,3)
1118 !  write(std_out,*) nkpt
1119 !  write(std_out,*) rprimd*Bohr_Ang
1120 !  write(std_out,*) two_pi*gprimd/Bohr_Ang
1121 !  write(std_out,*) mband
1122 !  write(std_out,*) "nwan",nwan
1123 !  write(std_out,*) nntot
1124 !  write(std_out,*) natom
1125 !  write(std_out,*) atom_symbols
1126 !  write(std_out,*) xcart
1127 !  write(std_out,*) num_bands,num_bands,nntot,nkpt
1128 !  write(std_out,*) wann_spreads
1129 !  wann_spreads=2
1130 !  do i=1, nkpt
1131 !  do j=1, nntot
1132 !  write(std_out,*) i,j
1133 !  do k=1, num_bands
1134 !  do l=1, num_bands
1135 !  write(std_out,*) "m",M_matrix(l,k,j,i,1)
1136 !  enddo
1137 !  enddo
1138 !  enddo
1139 !  enddo
1140 
1141 
1142 
1143 #if defined HAVE_WANNIER90
1144    do isppol=1,nsppol
1145      if(spin.ne.0 .and. spin.ne.isppol) cycle
1146 !    when nsppol>1, master runs isppol 1 and rank==1 runs isppol 2
1147      if(nprocs>1 .and. isppol==1.and.rank.ne.master) cycle
1148      if(nprocs>1 .and. isppol==2.and.rank.ne.1) cycle
1149 
1150      write(message, '(8a)' ) ch10,&
1151 &     '** mlwfovlp :   call wannier90 library subroutine wannier_run ',ch10,&
1152 &     '   Calculation is running         ',ch10,&
1153 &     '-  see ',trim(filew90_wout(isppol)),' for details.'
1154      call wrtout(std_out,  message,'COLL')
1155 !
1156      call wannier_run(trim(seed_name(isppol)),ngkpt,nkpt,&            !input
1157 &    real_lattice,recip_lattice,dtset%kpt,num_bands(isppol),& !input
1158 &    nwan(isppol),nntot,natom,atom_symbols,&                  !input
1159 &    xcart*Bohr_Ang,gamma_only,M_matrix(:,:,:,:,isppol),A_matrix(:,:,:,isppol),eigenvalues_w(:,:,isppol),& !input
1160 &    U_matrix(1:nwan(isppol),1:nwan(isppol),:,isppol),& !output
1161 &    U_matrix_opt(1:num_bands(isppol),1:nwan(isppol),:,isppol),& !output
1162 &    lwindow_loc=lwindow(1:num_bands(isppol),:,isppol),& !output
1163 &    wann_centres_loc=wann_centres(:,1:nwan(isppol),isppol),&     !output
1164 &    wann_spreads_loc=wann_spreads(1:nwan(isppol),isppol),spread_loc=spreadw(:,isppol))                            !output
1165 
1166 !    ----------------------------------------------------------------------------------------------
1167 
1168      write(message, '(7a)' ) ch10,&
1169 &     '   mlwfovlp :  mlwfovlp_run completed -',ch10,&
1170 &     '-  see ',trim(filew90_wout(isppol)),' for details.',ch10
1171      call wrtout(ab_out,message,'COLL')
1172      call wrtout(std_out,  message,'COLL')
1173 
1174    end do !isppol
1175 
1176 !  collect output of  wannier90 from different processors
1177    call xmpi_barrier(spaceComm)
1178 
1179    call xmpi_sum(U_matrix,spaceComm,ierr)
1180    call xmpi_sum(U_matrix_opt,spaceComm,ierr)
1181    call xmpi_lor(lwindow,spaceComm)
1182    call xmpi_sum(wann_centres,spaceComm,ierr)
1183    call xmpi_sum(wann_spreads,spaceComm,ierr)
1184 
1185    ! Output ABIWAN.nc file
1186 #ifdef HAVE_NETCDF
1187    if (dtset%kptopt == 0) then
1188      MSG_WARNING("Output of ABIWAN.nc requires kptopt /= 0. ABIWAN.nc file won't be produced!")
1189      ! Need kptrlatt in wigner_seitz and client code need to know the k-grid.
1190    end if
1191    if (rank == master .and. dtset%kptopt /= 0) then
1192      abiwan_fname = strcat(dtfil%filnam_ds(4), "_ABIWAN.nc")
1193      call wrtout(std_out, sjoin("Saving wannier90 ouput results in:", abiwan_fname))
1194      call wigner_seitz([zero, zero, zero], [2, 2, 2], dtset%kptrlatt, crystal%rmet, nrpts, irvec, ndegen)
1195      ! We know if disentanglement has been done by looking at the output values of lwindow
1196      ! Not elegant but it is the only way to avoid the parsing of the wannier input.
1197      ! In wannier_run lwindow is set to True if not disentanglement
1198      have_disentangled_spin = 0
1199      do isppol=1,nsppol
1200        !if nwan(isppol) < num_bands(isppol)
1201        if (.not. all(lwindow(:,:,isppol))) have_disentangled_spin(isppol) = 1
1202      end do
1203 
1204      NCF_CHECK(nctk_open_create(ncid, abiwan_fname, xmpi_comm_self))
1205      NCF_CHECK(hdr_ncwrite(hdr, ncid, fform_from_ext("ABIWAN"), nc_define=.True.))
1206      NCF_CHECK(crystal_ncwrite(crystal, ncid))
1207      NCF_CHECK(ebands_ncwrite(ebands, ncid))
1208 
1209      ncerr = nctk_def_dims(ncid, [ &
1210        nctkdim_t("mwan", mwan), &
1211        nctkdim_t("max_num_bands", max_num_bands), &
1212        nctkdim_t("nrpts", nrpts) &
1213      ], defmode=.True.)
1214      NCF_CHECK(ncerr)
1215 
1216      ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "nntot"])
1217      NCF_CHECK(ncerr)
1218      !ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: "fermi_energy", "smearing_width"])
1219      !NCF_CHECK(ncerr)
1220 
1221      ncerr = nctk_def_arrays(ncid, [ &
1222        nctkarr_t("nwan", "int", "number_of_spins"), &
1223        nctkarr_t("num_bands", "int", "number_of_spins"), &
1224        nctkarr_t("band_in_int", "int", "max_number_of_states, number_of_spins"), &
1225        nctkarr_t("lwindow_int", "int", "max_num_bands, number_of_kpoints, number_of_spins"), &
1226        !nctkarr_t("exclude_bands", "int", "max_number_of_states, number_of_spins"), &
1227        !nctkarr_t("eigenvalues_w", "int", "max_num_bands, number_of_kpoints, number_of_spins"), &
1228        nctkarr_t("spread", "dp", "three, number_of_spins"), &
1229        !nctkarr_t("A_matrix", "dp", "two, max_num_bands, mwan, number_of_kpoints, number_of_spins"), &
1230        nctkarr_t("irvec", "int", "three, nrpts"), &
1231        nctkarr_t("ndegen", "int", "nrpts"), &
1232        nctkarr_t("have_disentangled_spin", "int", "number_of_spins"), &
1233        nctkarr_t("U_matrix", "dp", "two, mwan, mwan, number_of_kpoints, number_of_spins"), &
1234        nctkarr_t("U_matrix_opt", "dp", "two, max_num_bands, mwan, number_of_kpoints, number_of_spins"), &
1235        nctkarr_t("wann_centres", "dp", "three, mwan, number_of_spins"), &
1236        nctkarr_t("wann_spreads", "dp", "mwan, number_of_spins") &
1237        ])
1238      NCF_CHECK(ncerr)
1239 
1240      ! Write data.
1241      NCF_CHECK(nctk_set_datamode(ncid))
1242      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "nntot"), nntot))
1243      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "nwan"), nwan))
1244      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "num_bands"), num_bands))
1245      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "band_in_int"), l2int(band_in)))
1246      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "lwindow_int"), l2int(lwindow)))
1247      !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "exclude_bands"), exclude_bands))
1248      !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "eigenvalues_w"), eigenvalues_w))
1249      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "spread"), spreadw))
1250      !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "A_matrix"), c2r(A_matrix)))
1251      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "irvec"), irvec))
1252      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ndegen"), ndegen))
1253      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "have_disentangled_spin"), have_disentangled_spin))
1254      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "U_matrix"), c2r(U_matrix)))
1255      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "U_matrix_opt"), c2r(U_matrix_opt)))
1256      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "wann_centres"), wann_centres))
1257      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "wann_spreads"), wann_spreads))
1258      NCF_CHECK(nf90_close(ncid))
1259 
1260      ABI_FREE(irvec)
1261      ABI_FREE(ndegen)
1262    end if
1263 #endif
1264 
1265 !  CALL SILVESTRELLI'S APPROACH TO EVALUATE vdW INTERACTION ENERGY USING MLWF!!
1266 !  ----------------------------------------------------------------------------------------------
1267    if (dtset%vdw_xc==10.or.dtset%vdw_xc==11.or.dtset%vdw_xc==12.or.dtset%vdw_xc==14.and.rank==master) then
1268 !    vdw_xc==10,11,12,14 starts the
1269 !    vdW interaction using MLWFs
1270      write(std_out,*) 'nwan(nsppol)=',ch10
1271      do ii=1,nsppol
1272        write(std_out,*) 'nsppol=',ii, 'nwan(nsppol)=',nwan(ii),ch10
1273      end do
1274      write(std_out,*) 'mwan=', mwan, ch10
1275 
1276      ABI_ALLOCATE(occ_arr,(mband,nkpt,isppol))
1277      ABI_ALLOCATE(occ_wan,(mwan,nkpt,nsppol))
1278      ABI_ALLOCATE(tdocc_wan,(mwan,nsppol))
1279 
1280      occ_arr(:,:,:)=zero
1281      occ_wan(:,:,:)=zero
1282      tdocc_wan(:,:)=zero
1283      jj = 0
1284      do isppol=1,nsppol
1285        do ikpt=1,nkpt
1286          do iband=1,num_bands(isppol)
1287            jj = jj + 1
1288            occ_arr(iband,ikpt,isppol) = occ(jj)
1289          end do
1290        end do
1291      end do
1292 
1293      do isppol=1,nsppol
1294        do ikpt=1,nkpt
1295          do iwan=1,nwan(isppol)
1296            caux=czero
1297            caux2=czero
1298            caux3=czero
1299            do iband=1,num_bands(isppol) !nband_inc(isppol) !nwan(isppol)
1300              do ii=1,nwan(isppol)
1301                caux=U_matrix(ii,iwan,ikpt,isppol)*U_matrix_opt(iband,ii,ikpt,isppol)
1302 !              DEBUG
1303 !              if(ISNAN(dble(caux))) then
1304 !              write(std_out,*) 'NaN: caux(ikpt,iwan,iband,ii):',ikpt,iwan,iband,ii,ch10
1305 !              end if
1306 !              END DEBUG
1307                do kk=1,nwan(isppol)
1308                  caux2=conjg(U_matrix(kk,iwan,ikpt,isppol))*conjg(U_matrix_opt(iband,kk,ikpt,isppol))
1309                  caux3= caux3+caux*caux2*occ_arr(iband,ikpt,isppol) !take care here as exclude_bands case is not well
1310 !                DEBUG
1311 !                if(ISNAN(dble(caux2))) then
1312 !                write(std_out,*) 'NaN: caux2(ikpt,iwan,iband,kk):',ikpt,iwan,iband,kk,ch10
1313 !                end if
1314 !                if(ISNAN(dble(caux3))) then
1315 !                write(std_out,*) 'NaN: caux3(ikpt,iwan,iband,kk,jj):',ikpt,iwan,iband,kk,jj
1316 !                end if
1317 !                END DEBUG
1318                end do
1319              end do
1320            end do
1321            occ_wan(iwan,ikpt,isppol) = dble(caux3)
1322 !          DEBUG
1323 !          write(std_out,*) occ_wan(iwan,ikpt,isppol)
1324 !          END DEBUG
1325 !          end do
1326          end do
1327        end do
1328      end do
1329 
1330      write(std_out,*) ch10,'MLWFs Occupation Matrix diagonal terms:',ch10
1331 
1332      do jj=1,nsppol
1333        forall(iwan=1:nwan(jj)) tdocc_wan(iwan,jj) = sum(occ_wan(iwan,1:nkpt,jj)) / real(nkpt,dp)
1334        write(std_out,*) 'tdocc_wan(iwan),isppol:',ch10
1335        write(std_out,*) (tdocc_wan(iwan,jj),iwan=1,nwan(jj)),jj
1336      end do
1337 
1338      ABI_ALLOCATE(csix,(mwan,mwan,nsppol,nsppol))
1339 
1340      call evdw_wannier(csix,corrvdw,mwan,natom,nsppol,nwan,tdocc_wan,dtset%vdw_nfrag,&
1341 &     dtset%vdw_supercell,dtset%vdw_typfrag,dtset%vdw_xc,rprimd,wann_centres,wann_spreads,xcart)
1342 
1343      ABI_DEALLOCATE(csix)
1344      ABI_DEALLOCATE(occ_arr)
1345      ABI_DEALLOCATE(occ_wan)
1346      ABI_DEALLOCATE(tdocc_wan)
1347    end if
1348 #else
1349    ABI_UNUSED(occ)
1350 #endif
1351 !  FIXME: looks like there is no automatic test which goes through here: g95 bot did not catch
1352 !  the missing deallocations
1353    ABI_DEALLOCATE(wann_centres)
1354    ABI_DEALLOCATE(wann_spreads)
1355    ABI_DEALLOCATE(U_matrix)
1356    ABI_DEALLOCATE(U_matrix_opt)
1357    ABI_DEALLOCATE(lwindow)
1358 
1359  end if !lwannierrun
1360 !
1361 !deallocation
1362 !
1363  ABI_DEALLOCATE(band_in)
1364  ABI_DEALLOCATE(atom_symbols)
1365  ABI_DEALLOCATE(xcart)
1366  ABI_DEALLOCATE(eigenvalues_w)
1367  ABI_DEALLOCATE(M_matrix)
1368  ABI_DEALLOCATE(A_matrix)
1369 
1370 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
  spin = just used for nsppol>1 ; 0 both, 1 just spin up, 2 just spin down

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

PARENTS

      mlwfovlp

CHILDREN

      mlwfovlp_radial,mlwfovlp_ylmfac,wrtout,ylm_cmplx

SOURCE

2351  subroutine mlwfovlp_proj(A_matrix,band_in,cg,cprj,dtset,gprimd,just_augmentation,kg,&
2352 &lproj,max_num_bands,mband,mkmem,mpi_enreg,mpw,mwan,natom,nattyp,&
2353 &nkpt,npwarr,nspinor,&
2354 &nsppol,ntypat,num_bands,nwan,pawtab,proj_l,proj_m,proj_radial,&
2355 &proj_site,proj_x,proj_z,proj_zona,psps,spin,ucvol)
2356 
2357 
2358 !This section has been created automatically by the script Abilint (TD).
2359 !Do not modify the following lines by hand.
2360 #undef ABI_FUNC
2361 #define ABI_FUNC 'mlwfovlp_proj'
2362 !End of the abilint section
2363 
2364  implicit none
2365 
2366 !Arguments ------------------------------------
2367 !scalars
2368  complex(dpc),parameter :: c1=(1._dp,0._dp)
2369  integer,intent(in) :: lproj,max_num_bands,mband,mkmem,mpw,mwan,natom,nkpt,nspinor,nsppol
2370  integer,intent(in) :: ntypat,spin
2371  type(MPI_type),intent(in) :: mpi_enreg
2372  type(dataset_type),intent(in) :: dtset
2373  type(pseudopotential_type),intent(in) :: psps
2374 !arrays
2375  integer ::nattyp(ntypat)
2376  integer,intent(in) :: kg(3,mpw*mkmem),npwarr(nkpt),num_bands(nsppol),nwan(nsppol),proj_l(mband,nsppol)
2377  integer,intent(in) :: proj_m(mband,nsppol)
2378  integer,intent(inout)::proj_radial(mband,nsppol)
2379  real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol)
2380  real(dp),intent(in) :: gprimd(3,3),proj_site(3,mband,nsppol)
2381  real(dp),intent(in) :: proj_x(3,mband,nsppol),proj_z(3,mband,nsppol),proj_zona(mband,nsppol)
2382  complex(dpc),intent(out) :: A_matrix(max_num_bands,mwan,nkpt,nsppol)
2383 !character(len=fnlen),intent(in) :: filew90_win(nsppol)
2384  logical,intent(in) :: band_in(mband,nsppol)
2385  logical,intent(in)::just_augmentation(mwan,nsppol)
2386  type(pawcprj_type) :: cprj(natom,nspinor*mband*mkmem*nsppol)
2387  type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
2388 
2389 !Local variables-------------------------------
2390 !scalars
2391  integer :: iatom,iatprjn,iband,iband1,iband2,ibg,icat,icg,icg_shift
2392  integer :: idum,idx,ikg,ikpt,ilmn,ipw,iproj
2393  integer :: ispinor,isppol,itypat,iwan,jband,jj1,libprjn
2394  integer :: lmn_size,natprjn,nband_k,nbprjn,npw_k
2395  integer :: sumtmp
2396  integer :: max_lmax,max_lmax2,mproj,nprocs,spaceComm,rank
2397  real(dp),parameter :: qtol=2.0d-8
2398  real(dp) :: arg,norm_error,norm_error_bar
2399  real(dp) :: ucvol,x1,x2,xnorm,xnormb,xx,yy,zz
2400  complex(dpc) :: amn_tmp(nspinor)
2401  complex(dpc) :: cstr_fact
2402  character(len=500) :: message
2403 !arrays
2404  integer :: kg_k(3,mpw),lmax(nsppol),lmax2(nsppol),nproj(nsppol)
2405  integer,allocatable :: lprjn(:),npprjn(:)
2406  real(dp) :: kpg(3),kpt(3)
2407  real(dp),allocatable :: amn(:,:,:,:,:),amn2(:,:,:,:,:,:,:)
2408  real(dp),allocatable :: gsum2(:),kpg2(:),radial(:)
2409  complex(dpc),allocatable :: gf(:,:),gft_lm(:)
2410  complex(dpc),allocatable :: ylmc_fac(:,:,:),ylmcp(:)
2411 
2412 !no_abirules
2413 !Tables 3.1 & 3.2, User guide
2414  integer,save :: orb_l_defs(-5:3)=(/2,2,1,1,1,0,1,2,3/)
2415 ! integer,parameter :: mtransfo(0:3,7)=&
2416 !&  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/))
2417 
2418 !************************************************************************
2419 
2420 !mpi initialization
2421  spaceComm=MPI_enreg%comm_cell
2422  nprocs=xmpi_comm_size(spaceComm)
2423  rank=MPI_enreg%me_kpt
2424 
2425 !Check input variables
2426  if ((lproj/=1).and.(lproj/=2).and.(lproj/=5)) then
2427    write(message, '(3a)' )&
2428 &   ' Value of lproj no allowed ',ch10,&
2429 &   ' Action : change lproj.'
2430    MSG_ERROR(message)
2431  end if
2432 
2433  write(message, '(a,a)' )ch10,&
2434 & '** mlwfovlp_proj:  compute A_matrix of initial guess for wannier functions'
2435  call wrtout(std_out,message,'COLL')
2436 
2437 !Initialize to 0.d0
2438  A_matrix(:,:,:,:)=cmplx(0.d0,0.d0)
2439 
2440 !End of preliminarities
2441 
2442 !********************* Write Random projectors
2443  if(lproj==1) then
2444    idum=123456
2445 !  Compute random projections
2446    ABI_ALLOCATE(amn,(2,mband,mwan,nkpt,nsppol))
2447    amn=zero
2448    do isppol=1,nsppol
2449      if(spin.ne.0 .and. spin.ne.isppol) cycle
2450      do ikpt=1,nkpt
2451 !
2452 !      MPI: cycle over kpts not treated by this node
2453 !
2454        if (ABS(MPI_enreg%proc_distrb(ikpt,1,isppol)-rank)/=0) CYCLE
2455 !      write(std_out,'("kpt loop2: ikpt",i3," rank ",i3)') ikpt,rank
2456 
2457 !
2458        do iband1=1,mband
2459          xnormb=0.d0
2460          do iband2=1,nwan(isppol)
2461            x1=uniformrandom(idum)
2462            x2=uniformrandom(idum)
2463            xnorm=sqrt(x1**2+x2**2)
2464            xnormb=xnormb+xnorm
2465            amn(1,iband1,iband2,ikpt,isppol)=x1
2466            amn(2,iband1,iband2,ikpt,isppol)=x2
2467          end do
2468          do iband2=1,nwan(isppol)
2469            amn(1,iband1,iband2,ikpt,isppol)=amn(1,iband1,iband2,ikpt,isppol)/xnormb
2470            amn(2,iband1,iband2,ikpt,isppol)=amn(2,iband1,iband2,ikpt,isppol)/xnormb
2471          end do !iband2
2472        end do !iband1
2473      end do !ikpt
2474    end do !isppol
2475    do isppol=1,nsppol
2476      if(spin.ne.0 .and. spin.ne.isppol) cycle
2477      do ikpt=1,nkpt
2478 !
2479 !      MPI: cycle over kpts not treated by this node
2480 !
2481        if (ABS(MPI_enreg%proc_distrb(ikpt,1,isppol)-rank)/=0) CYCLE
2482 !
2483        do iband2=1,nwan(isppol)
2484          jband=0
2485          do iband1=1,mband
2486            if(band_in(iband1,isppol)) then
2487              jband=jband+1
2488              if(jband.gt.num_bands(isppol)) then
2489                write(message, '(3a)' )&
2490 &               'Value of jband is above num_bands ',ch10,&
2491 &               'Action: contact Abinit group'
2492                MSG_ERROR(message)
2493              end if
2494              A_matrix(jband,iband2,ikpt,isppol)=cmplx(amn(1,iband1,iband2,ikpt,isppol),amn(2,iband1,iband2,ikpt,isppol))
2495            end if
2496          end do !iband1
2497        end do !iband2
2498      end do !ikpt
2499    end do !isppol
2500    ABI_DEALLOCATE(amn)
2501  end if
2502 
2503 !********************* Projection on atomic orbitals based on .win file
2504  if( lproj==2) then !based on .win file
2505    nproj(:)=nwan(:)/nspinor !if spinors, then the number of projections are
2506    mproj=maxval(nproj(:))
2507 !  half the total of wannier functions
2508 !
2509 !  obtain lmax and lmax2
2510    lmax(:)=0
2511    lmax2(:)=0
2512 !
2513    do isppol=1,nsppol
2514      if(spin.ne.0 .and. spin.ne.isppol) cycle
2515      do iproj=1,nproj(isppol)
2516        lmax(isppol)=max(lmax(isppol),orb_l_defs(proj_l(iproj,isppol)))
2517      end do !iproj
2518      lmax2(isppol)=(lmax(isppol)+1)**2
2519    end do !isppol
2520    max_lmax=maxval(lmax(:))
2521    max_lmax2=maxval(lmax2(:))
2522 !  Allocate arrays
2523    ABI_ALLOCATE(ylmc_fac,(max_lmax2,mproj,nsppol))
2524 !
2525 !  get ylmfac, factor used for rotations and hybrid orbitals
2526    do isppol=1,nsppol
2527      if(spin.ne.0 .and. spin.ne.isppol) cycle
2528      call mlwfovlp_ylmfac(ylmc_fac(1:lmax2(isppol),1:nproj(isppol),isppol),lmax(isppol),lmax2(isppol),&
2529 &     mband,nproj(isppol),proj_l(:,isppol),proj_m(:,isppol),proj_x(:,:,isppol),proj_z(:,:,isppol))
2530    end do
2531 !
2532    norm_error=zero
2533    norm_error_bar=zero
2534    icg=0
2535 !
2536    do isppol=1,nsppol
2537 !    Allocate arrays
2538      if(spin.eq.0 .or. spin.eq.isppol) then
2539 !      this has to be done this way because the variable icg changes at the end of the
2540 !      cycle. We cannot just skip the hole cycle.
2541        ABI_ALLOCATE(gf,(mpw,nproj(isppol)))
2542        ABI_ALLOCATE(gft_lm,(lmax2(isppol)))
2543        ABI_ALLOCATE(gsum2,(nproj(isppol)))
2544        ABI_ALLOCATE(kpg2,(mpw))
2545        ABI_ALLOCATE(radial,(lmax2(isppol)))
2546        ABI_ALLOCATE(ylmcp,(lmax2(isppol)))
2547      end if
2548 !
2549      ikg=0
2550      do ikpt=1, nkpt
2551 !
2552 !      MPI: cycle over kpts not treated by this node
2553 !
2554        if (ABS(MPI_enreg%proc_distrb(ikpt,1,isppol)-rank)/=0) CYCLE
2555 !
2556        if(spin.eq.0 .or. spin.eq.isppol) then
2557          write(message, '(a,i6,a,2i6)' ) &
2558 &         '   processor',rank,' will compute k-point,spin=',ikpt,isppol
2559          call wrtout(std_out,  message,'COLL')
2560        end if
2561 !
2562 !      Initialize variables
2563        npw_k=npwarr(ikpt)
2564        gsum2(:)=0.d0
2565        gf(:,:) = (0.d0,0.d0)
2566        kpt(:)=dtset%kpt(:,ikpt)
2567        kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
2568 
2569        do ipw=1, npw_k
2570          kpg(1)= (kpt(1) + real(kg_k(1,ipw),dp))     !k+G
2571          kpg(2)= (kpt(2) + real(kg_k(2,ipw),dp))
2572          kpg(3)= (kpt(3) + real(kg_k(3,ipw),dp))
2573 !
2574 !        Calculate modulus of k+G
2575          xx=gprimd(1,1)*kpg(1)+gprimd(1,2)*kpg(2)+gprimd(1,3)*kpg(3)
2576          yy=gprimd(2,1)*kpg(1)+gprimd(2,2)*kpg(2)+gprimd(2,3)*kpg(3)
2577          zz=gprimd(3,1)*kpg(1)+gprimd(3,2)*kpg(2)+gprimd(3,3)*kpg(3)
2578          kpg2(ipw)= two_pi*sqrt(xx**2+yy**2+zz**2)
2579 !
2580 !        Complex Y_lm for k+G
2581          if(lmax(isppol)==0) then
2582            ylmcp(1)=c1/sqrt(four_pi)
2583          else
2584            call ylm_cmplx(lmax(isppol),ylmcp,xx,yy,zz)
2585          end if
2586 !
2587          if(spin.eq.0 .or. spin.eq.isppol) then
2588 !          !
2589            do iproj=1,nproj(isppol)
2590 !
2591 !            In PAW, we can use proj_radial > 4 to indicate that we just
2592 !            want the in-sphere contribution
2593 !
2594              if( psps%usepaw==1) then
2595                if( just_augmentation(iproj,isppol)) cycle
2596              end if
2597 !
2598 !            obtain radial part
2599              call mlwfovlp_radial(proj_zona(iproj,isppol),lmax(isppol),lmax2(isppol)&
2600 &             ,radial,proj_radial(iproj,isppol),kpg2(ipw))
2601 !
2602 !            scale complex representation of projector orbital with radial functions
2603 !            of appropriate l
2604              gft_lm(:)=radial(:)*ylmc_fac(1:lmax2(isppol),iproj,isppol)
2605 !
2606 !            complex structure factor for projector orbital position
2607              arg = ( kpg(1)*proj_site(1,iproj,isppol) + &
2608 &             kpg(2)*proj_site(2,iproj,isppol) + &
2609 &             kpg(3)*proj_site(3,iproj,isppol) ) * 2*pi
2610              cstr_fact = cmplx(cos(arg), -sin(arg) )
2611 !
2612 !            obtain guiding functions
2613              gf(ipw,iproj)=cstr_fact*dot_product(ylmcp,gft_lm)
2614 !
2615              gsum2(iproj)=gsum2(iproj)+real(gf(ipw,iproj))**2+aimag(gf(ipw,iproj))**2
2616            end do !iproj
2617          end if !spin
2618        end do !ipw
2619 !
2620        if(spin.eq.0 .or. spin.eq.isppol) then
2621          do iproj=1,nproj(isppol)
2622 !
2623 !          In PAW, we can use proj_radial > 4 to indicate that we just
2624 !          want the in-sphere contribution
2625 !
2626            if(psps%usepaw==1 ) then
2627              if (just_augmentation(iproj,isppol)) cycle
2628            end if
2629 !
2630            gsum2(iproj)=16._dp*pi**2*gsum2(iproj)/ucvol
2631            gf(:,iproj)=gf(:,iproj)/sqrt(gsum2(iproj))
2632            norm_error=max(abs(gsum2(iproj)-one),norm_error)
2633            norm_error_bar=norm_error_bar+(gsum2(iproj)-one)**2
2634          end do !iproj
2635 !
2636 !        Guiding functions are computed.
2637 !        compute overlaps of gaussian projectors and wave functions
2638          do iproj=1,nproj(isppol)
2639 !
2640 !          In PAW, we can use proj_radial > 4 to indicate that we just
2641 !          want the in-sphere contribution
2642 !
2643            if(psps%usepaw==1 ) then
2644              if ( just_augmentation(iproj,isppol)) cycle
2645            end if
2646 !
2647            jband=0
2648            do iband=1,mband
2649              if(band_in(iband,isppol)) then
2650                icg_shift=npw_k*nspinor*(iband-1)+icg
2651                jband=jband+1
2652                amn_tmp(:)=cmplx(0.d0,0.d0)
2653                do ispinor=1,nspinor
2654                  do ipw=1,npw_k
2655 !
2656 !                  The case of spinors is tricky, we have nproj =  nwan/2
2657 !                  so we project to spin up and spin down separately, to have at
2658 !                  the end an amn matrix with nwan projections.
2659 !
2660 !                  idx=ipw*nspinor - (nspinor-ispinor)
2661                    idx=ipw+(ispinor-1)*npw_k
2662                    amn_tmp(ispinor)=amn_tmp(ispinor)+gf(ipw,iproj)*cmplx(cg(1,idx+icg_shift),-cg(2,idx+icg_shift))
2663                  end do !ipw
2664                end do !ispinor
2665                do ispinor=1,nspinor
2666                  iwan=(iproj*nspinor)- (nspinor-ispinor)
2667                  A_matrix(jband,iwan,ikpt,isppol)=amn_tmp(ispinor)
2668                end do
2669              end if !band_in
2670            end do !iband
2671          end do !iproj
2672        end if !spin==isppol
2673        icg=icg+npw_k*nspinor*mband
2674        ikg=ikg+npw_k
2675      end do !ikpt
2676 !    Deallocations
2677      if(spin.eq.0 .or. spin.eq.isppol) then
2678        ABI_DEALLOCATE(gf)
2679        ABI_DEALLOCATE(gft_lm)
2680        ABI_DEALLOCATE(gsum2)
2681        ABI_DEALLOCATE(kpg2)
2682        ABI_DEALLOCATE(radial)
2683        ABI_DEALLOCATE(ylmcp)
2684      end if
2685    end do !isppol
2686 !
2687 !  if(isppol==1) then
2688 !  norm_error_bar=sqrt(norm_error_bar/real(nkpt*(nwan(1)),dp))
2689 !  else
2690 !  if(spin==0)    norm_error_bar=sqrt(norm_error_bar/real(nkpt*(nwan(1)+nwan(2)),dp))
2691 !  if(spin==1)    norm_error_bar=sqrt(norm_error_bar/real(nkpt*nwan(1),dp))
2692 !  if(spin==2)    norm_error_bar=sqrt(norm_error_bar/real(nkpt*nwan(2),dp))
2693 !  end if
2694 !  if(norm_error>0.05_dp) then
2695 !  write(message, '(6a,f6.3,a,f6.3,12a)' )ch10,&
2696 !  &     ' mlwfovlp_proj : WARNING',ch10,&
2697 !  &     '  normalization error for wannier projectors',ch10,&
2698 !  &     '  is',norm_error_bar,' (average) and',norm_error,' (max).',ch10,&
2699 !  &     '  this may indicate more cell-to-cell overlap of the radial functions',ch10,&
2700 !  &     '  than you want.',ch10,&
2701 !  &     '  Action : modify zona (inverse range of radial functions)',ch10,&
2702 !  '  under "begin projectors" in ',trim(filew90_win),' file',ch10
2703 !  call wrtout(std_out,message,'COLL')
2704 !  end if
2705 !
2706 !  !Deallocate
2707 !  deallocate(ylmc_fac)
2708 !
2709    ABI_DEALLOCATE(ylmc_fac)
2710  end if !lproj==2
2711 
2712 
2713 !*************** computes projection  from PROJECTORS ********************
2714  if(lproj==3) then  !! if LPROJPRJ
2715 !  ----- set values for projections --------------------- ! INPUT
2716 !  nbprjn:number of  different l-values for projectors
2717 !  lprjn: value of l for each projectors par ordre croissant
2718 !  npprjn: number of projectors for each lprjn
2719    natprjn=1  ! atoms with wannier functions are first
2720    if(natprjn/=1) then ! in this case lprjn should depend on iatprjn
2721      MSG_ERROR("natprjn/=1")
2722    end if
2723    nbprjn=2
2724    ABI_ALLOCATE(lprjn,(nbprjn))
2725    lprjn(1)=0
2726    lprjn(2)=1
2727    ABI_ALLOCATE(npprjn,(0:lprjn(nbprjn)))
2728    npprjn(0)=1
2729    npprjn(1)=1
2730 !  --- test coherence of nbprjn and nwan
2731    sumtmp=0
2732    do iatprjn=1,natprjn
2733      do libprjn=0,lprjn(nbprjn)
2734        sumtmp=sumtmp+(2*libprjn+1)*npprjn(libprjn)
2735      end do
2736    end do
2737    if(sumtmp/=nwan(1)) then
2738      write(std_out,*) "Number of Wannier orbitals is not equal to number of projections"
2739      write(std_out,*) "Action: check values of lprjn,npprjn % nwan"
2740      write(std_out,*) "nwan, sumtmp=",nwan,sumtmp
2741      MSG_ERROR("Aborting now")
2742    end if
2743 !  --- end test of coherence
2744    ABI_ALLOCATE(amn2,(2,natom,nsppol,nkpt,mband,nspinor,nwan(1)))
2745    if(psps%usepaw==1) then
2746      amn2=zero
2747      ibg=0
2748      do isppol=1,nsppol
2749        do ikpt=1,nkpt
2750          nband_k=dtset%nband(ikpt+(isppol-1)*nkpt)
2751          do iband=1,nband_k
2752 !          write(std_out,*)"amn2",iband,ibg,ikpt
2753            do ispinor=1,nspinor
2754              icat=1
2755              do itypat=1,dtset%ntypat
2756                lmn_size=pawtab(itypat)%lmn_size
2757                do iatom=icat,icat+nattyp(itypat)-1
2758                  jj1=0
2759                  do ilmn=1,lmn_size
2760                    if(iatom.le.natprjn) then
2761 !                    do iwan=1,nwan
2762                      do libprjn=0,lprjn(nbprjn)
2763 !                      if (psps%indlmn(1,ilmn,itypat)==proj_l(iwan)) then
2764 !                      if (psps%indlmn(2,ilmn,itypat)==mtransfo(proj_l(iwan),proj_m(iwan))) then
2765                        if (psps%indlmn(1,ilmn,itypat)==libprjn) then
2766                          if (psps%indlmn(3,ilmn,itypat)<=npprjn(libprjn)) then
2767                            if(band_in(iband,isppol)) then
2768                              jj1=jj1+1
2769                              if(jj1>nwan(isppol)) then
2770                                write(std_out,*) "number of wannier orbitals is lower than lmn_size"
2771                                write(std_out,*) jj1,nwan(isppol)
2772                                MSG_ERROR("Aborting now")
2773                              end if
2774                              amn2(1,iatom,isppol,ikpt,iband,ispinor,jj1)=cprj(iatom,iband+ibg)%cp(1,ilmn)
2775                              amn2(2,iatom,isppol,ikpt,iband,ispinor,jj1)=cprj(iatom,iband+ibg)%cp(2,ilmn)
2776                            end if
2777                          end if
2778                        end if
2779                      end do ! libprjn
2780 !                    endif
2781 !                    endif
2782 !                    enddo ! iwan
2783                    end if ! natprjn
2784                  end do !ilmn
2785                end do ! iatom
2786                icat=icat+nattyp(itypat)
2787              end do ! itypat
2788            end do ! ispinor
2789          end do !iband
2790          ibg=ibg+nband_k*nspinor
2791 !        write(std_out,*)'amn2b',iband,ibg,ikpt
2792        end do !ikpt
2793      end do ! isppol
2794 
2795 !    -----------------------  Save Amn   --------------------
2796      do isppol=1,nsppol
2797        do ikpt=1,nkpt
2798          do iband2=1,nwan(isppol)
2799            jband=0
2800            do iband1=1,mband
2801              if(band_in(iband1,isppol)) then
2802                jband=jband+1
2803                A_matrix(jband,iband2,ikpt,isppol)=&
2804 &               cmplx(amn2(1,1,1,ikpt,iband1,1,iband2),amn2(2,1,1,ikpt,iband1,1,iband2))
2805              end if
2806            end do
2807          end do
2808        end do
2809      end do
2810    end if !usepaw
2811    ABI_DEALLOCATE(amn2)
2812    ABI_DEALLOCATE(npprjn)
2813    ABI_DEALLOCATE(lprjn)
2814 
2815  end if ! lproj==3
2816 
2817 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.
  spin = just used for nsppol>1 ; 0 both, 1 just spin up, 2 just spin down
  typat(natom)= atom type
  xred(3,natom)=reduced dimensionless atomic coordinates

OUTPUT

  A_paw(max_num_bands,nwan,nkpt) = A matrix containing initial guess for MLWFs
                          (augmentation part of the matrix)

SIDE EFFECTS

NOTES

 This routine is still under developement

PARENTS

      mlwfovlp

CHILDREN

      mlwfovlp_ylmfar,simp_gen,simpson_int,wrtout,xred2xcart

SOURCE

2877 subroutine mlwfovlp_projpaw(A_paw,band_in,cprj,just_augmentation,max_num_bands,mband,mkmem,&
2878 &mwan,natom,nband,nkpt,&
2879 &nspinor,nsppol,ntypat,nwan,pawrad,pawtab,&
2880 &proj_l,proj_m,proj_radial,proj_site,proj_x,proj_z,proj_zona,psps,&
2881 &rprimd,spin,typat,xred)
2882 
2883 
2884 !This section has been created automatically by the script Abilint (TD).
2885 !Do not modify the following lines by hand.
2886 #undef ABI_FUNC
2887 #define ABI_FUNC 'mlwfovlp_projpaw'
2888 !End of the abilint section
2889 
2890  implicit none
2891 
2892 !Arguments ------------------------------------
2893  integer,intent(in) :: max_num_bands,mband,mkmem,mwan,natom,nkpt
2894  integer,intent(in) :: nspinor,nsppol,ntypat,spin
2895  !arrays
2896  integer,intent(in) :: nband(nsppol*nkpt),nwan(nsppol)
2897  integer,intent(in) :: proj_l(mband,nsppol),proj_m(mband,nsppol),proj_radial(mband,nsppol)
2898  integer,intent(in) :: typat(natom)
2899  real(dp),intent(in):: proj_site(3,mband,nsppol)
2900  real(dp),intent(in) :: proj_x(3,mband,nsppol),proj_z(3,mband,nsppol),proj_zona(mband,nsppol)
2901  real(dp),intent(in) :: rprimd(3,3),xred(3,natom)
2902  complex(dpc),intent(out) :: A_paw(max_num_bands,mwan,nkpt,nsppol)
2903  logical,intent(in) :: band_in(mband,nsppol)
2904  logical,intent(in)::just_augmentation(mwan,nsppol)
2905  type(pawcprj_type) :: cprj(natom,nspinor*mband*mkmem*nsppol)
2906  type(pawrad_type),intent(in) :: pawrad(ntypat)
2907  type(pawtab_type),intent(in) :: pawtab(ntypat)
2908  type(pseudopotential_type),intent(in) :: psps
2909 
2910 !Local variables-------------------------------
2911  !local variables
2912  integer :: basis_size,iatom,iband,ii
2913  integer :: ikpt,ir,isppol,itypat,iwan,jband
2914  integer :: ll,lm,ln,mm,ilmn
2915  integer :: lmn_size,max_lmax2, mesh_size,nn
2916  integer :: lmax(nsppol),lmax2(nsppol)
2917  real(dp):: aa,int_rad2,prod_real,prod_imag
2918  real(dp),parameter :: dx=0.015d0,rmax=10.d0,xmin=0.d0
2919  real(dp):: sum,wan_lm_fac,x
2920  complex(dpc)::prod
2921  character(len=500) :: message
2922  !arrays
2923  integer :: index(mband,nkpt,nsppol)
2924  real(dp):: dist,norm(mwan,nsppol)
2925  real(dp):: proj_cart(3,mwan,nsppol),proj_site_unit(3,mwan,nsppol)
2926  real(dp):: xcart_unit(3,natom),xred_unit(3,natom)
2927  real(dp),allocatable ::aux(:),ff(:),r(:),int_rad(:),rad_int(:)
2928  real(dp),allocatable::ylmr_fac(:,:,:)
2929 
2930 !no_abirules
2931 !Tables 3.1 & 3.2, User guide
2932  integer,save :: orb_l_defs(-5:3)=(/2,2,1,1,1,0,1,2,3/)
2933 !real(dp),allocatable :: ylm(:,:)
2934 
2935 
2936 ! *************************************************************************
2937 
2938  write(message, '(a,a)' )ch10,'** mlwfovlp_proj:  compute in-sphere part of A_matrix'
2939  call wrtout(std_out,message,'COLL')
2940 
2941 !Check input variables
2942  do isppol=1,nsppol
2943    if(spin.ne.0 .and. spin.ne.isppol) cycle
2944    do iwan=1,nwan(nsppol)
2945      if(proj_radial(iwan,isppol)<1 .or. proj_radial(iwan,isppol)>4)then
2946        write(message,'(a,a,a,i6)')&
2947 &       '  proj_radial should be between 1 and 4,',ch10,&
2948 &       '  however, proj_radial=',proj_radial(iwan,isppol)
2949        MSG_BUG(message)
2950      end if
2951    end do
2952  end do
2953 
2954 !Initialize
2955  A_paw(:,:,:,:)=cmplx(0.d0,0.d0)
2956 
2957 !Get index for cprj
2958  ii=0
2959  do isppol=1,nsppol
2960    do ikpt=1,nkpt
2961      do iband=1,nband(ikpt)
2962        ii=ii+1
2963        index(iband,ikpt,isppol)=ii
2964      end do
2965    end do
2966  end do
2967 
2968 !obtain lmax and lmax2
2969  lmax(:)=0
2970  lmax2(:)=0
2971  do isppol=1,nsppol
2972    if(spin.ne.0 .and. spin.ne.isppol) cycle
2973    do iwan=1,nwan(isppol)
2974      lmax(isppol)=max(lmax(isppol),orb_l_defs(proj_l(iwan,isppol)))
2975    end do !iwan
2976    lmax2(isppol)=(lmax(isppol)+1)**2
2977  end do
2978  max_lmax2=maxval(lmax2(:))
2979 !
2980 !get ylmfac, factor used for rotations and hybrid orbitals
2981 !
2982  ABI_ALLOCATE(ylmr_fac,(max_lmax2,mwan,nsppol))
2983 
2984 
2985  do isppol=1,nsppol
2986    if(spin.ne.0 .and. spin.ne.isppol) cycle
2987    call mlwfovlp_ylmfar(ylmr_fac(1:lmax2(isppol),1:nwan(isppol),isppol),&
2988 &   lmax(isppol),lmax2(isppol),mband,nwan(isppol),proj_l(:,isppol),proj_m(:,isppol),&
2989 &   proj_x(:,:,isppol),proj_z(:,:,isppol))
2990 !
2991 !  Shift projection centers and atom centers to the primitive cell
2992 !  This will be useful after, when we check if the Wannier function
2993 !  lies on one specific atom
2994 !
2995    proj_site_unit(:,:,:)=0.d0
2996    do iwan=1,nwan(isppol)
2997      do ii=1,3
2998        proj_site_unit(ii,iwan,isppol)=ABS(proj_site(ii,iwan,isppol)-AINT(proj_site(ii,iwan,isppol)) )
2999      end do
3000    end do
3001    do iatom=1,natom
3002      do ii=1,3
3003        xred_unit(ii,iatom)=ABS(xred(ii,iatom)-AINT(xred(ii,iatom)) )
3004      end do
3005    end do
3006    call xred2xcart(natom,rprimd,xcart_unit,xred_unit)
3007    call xred2xcart(mwan,rprimd,proj_cart(:,:,isppol),proj_site_unit(:,:,isppol))
3008 !
3009 !  Normalize the Wannier functions
3010 !
3011 !  Radial part
3012    mesh_size= nint((rmax - xmin ) / dx + 1)
3013    ABI_ALLOCATE( ff,(mesh_size))
3014    ABI_ALLOCATE(r,(mesh_size))
3015    ABI_ALLOCATE(rad_int,(mesh_size))
3016    ABI_ALLOCATE(aux,(mesh_size))
3017    do ir=1, mesh_size
3018      x=xmin+DBLE(ir-1)*dx
3019      r(ir)=x
3020    end do   !ir
3021    do iwan=1,nwan(isppol)
3022 !    write(std_out,*)'iwan',iwan
3023 !    radial functions shown in table 3.3 of wannier90 manual
3024      if(proj_radial(iwan,isppol)==1) ff(:) = 2.d0 * proj_zona(iwan,isppol)**(1.5d0) * exp(-proj_zona(iwan,isppol)*r(:))
3025      if(proj_radial(iwan,isppol)==2) ff(:) = 1.d0/(2.d0*sqrt(2.d0))*proj_zona(iwan,isppol)**(1.5d0) *&
3026 &     (2.d0 - proj_zona(iwan,isppol)*r(:))*exp(-proj_zona(iwan,isppol)*r(:)/2.d0)
3027      if(proj_radial(iwan,isppol)==3) ff(:) = sqrt(4.d0/27.d0)*proj_zona(iwan,isppol)**(1.5d0)&
3028 &     * (1.d0 - 2.d0*proj_zona(iwan,isppol)*r(:)/3.d0 + 2.d0*proj_zona(iwan,isppol)**2*r(:)**2/27.d0)&
3029 &     * exp(-proj_zona(iwan,isppol) * r(:)/3.d0)
3030 
3031      if(proj_radial(iwan,isppol)/=4) then
3032        aux(:)=ff(:)**2*r(:)**2
3033        call simpson_int(mesh_size,dx,aux,rad_int)
3034        sum=0.d0
3035        do ir=1,mesh_size
3036          sum=sum+rad_int(ir)
3037        end do
3038        int_rad2=sum/real(mesh_size,dp)
3039 !
3040 !      do ir=1,mesh_size
3041 !      if(iwan==1) write(400,*)r(ir),aux(ir),rad_int(ir)
3042 !      end do
3043      else
3044 !
3045 !      ==4: gaussian function
3046 !      f(x)=\exp(-1/4(x/aa)**2)
3047 !      \int f(x)f(x) dx = \int \exp(-1/2(x/aa)**2) = aa*sqrt(2pi)
3048 !
3049        int_rad2=sqrt(2.d0*pi)*proj_zona(iwan,isppol)
3050      end if
3051 
3052 !
3053 !    Now angular part
3054 !
3055      prod_real=0.d0
3056      do lm=1,lmax2(isppol)
3057        wan_lm_fac=ylmr_fac(lm,iwan,isppol)
3058 !      write(std_out,*)'wan_lm_fac',wan_lm_fac
3059 !      write(std_out,*)'int_rad2',int_rad2
3060        prod_real= prod_real + wan_lm_fac**2 * int_rad2
3061      end do
3062      norm(iwan,isppol)=sqrt(prod_real)
3063    end do !iwan
3064    ABI_DEALLOCATE(ff)
3065    ABI_DEALLOCATE(r)
3066    ABI_DEALLOCATE(rad_int)
3067    ABI_DEALLOCATE(aux)
3068 !
3069 !  Now that we found our guiding functions
3070 !  We proceed with the internal product of
3071 !  our guiding functions and the wave function
3072 !  Amn=<G_m|\Psi_n> inside the sphere.
3073 !  The term <G_m|\Psi_n> inside the sphere is:
3074 !  = \sum_i <G_n | \phi_i - \tphi_i> <p_im|\Psi_m>
3075 !
3076 !
3077 !  G_n \phi and \tphi can be decomposed in
3078 !  a radial function times an angular function.
3079 !
3080 !
3081 !  Big loop on iwan and iatom
3082 !
3083    do iwan=1,nwan(isppol)
3084      do iatom=1,natom
3085 !
3086 !      check if center of wannier function coincides
3087 !      with the center of the atom
3088 !
3089        dist=((proj_cart(1,iwan,isppol)-xcart_unit(1,iatom))**2 + &
3090 &       (proj_cart(2,iwan,isppol)-xcart_unit(2,iatom))**2 + &
3091 &       (proj_cart(3,iwan,isppol)-xcart_unit(3,iatom))**2)**0.5
3092 !
3093 !      if the distance between the centers is major than 0.1 angstroms skip
3094 !
3095        if( dist > 0.188972613) cycle
3096 !
3097        write(message, '(2a,i4,a,i4,2a)')ch10, '   Wannier function center',iwan,' is on top of atom',&
3098 &       iatom,ch10,'      Calculating in-sphere contribution'
3099        call wrtout(ab_out,message,'COLL')
3100        call wrtout(std_out,  message,'COLL')
3101 !
3102 !      Get useful quantities
3103 !
3104        itypat=typat(iatom)
3105        lmn_size=pawtab(itypat)%lmn_size
3106        basis_size=pawtab(itypat)%basis_size
3107        mesh_size=pawtab(itypat)%mesh_size
3108        ABI_ALLOCATE(int_rad,(basis_size))
3109        ABI_ALLOCATE(ff,(mesh_size))
3110        ABI_ALLOCATE(aux,(mesh_size))
3111 
3112 !
3113 !      Integrate first the radial part
3114 !      and save it into an array
3115 !
3116 !
3117 !      radial functions shown in table 3.3 of wannier90 manual
3118 !
3119        if(proj_radial(iwan,isppol)==1) aux(1:mesh_size) = 2.d0 * proj_zona(iwan,isppol)**(1.5d0) *&
3120 &       exp(-proj_zona(iwan,isppol)*pawrad(itypat)%rad(1:mesh_size))
3121        if(proj_radial(iwan,isppol)==2) aux(1:mesh_size) = 1.d0/(2.d0*sqrt(2.d0))*proj_zona(iwan,isppol)**(1.5d0) *&
3122 &       (2.d0 - proj_zona(iwan,isppol)*pawrad(itypat)%rad(1:mesh_size)) &
3123 &       * exp(-proj_zona(iwan,isppol)*pawrad(itypat)%rad(1:mesh_size)/2.d0)
3124        if(proj_radial(iwan,isppol)==3) aux(1:mesh_size) = sqrt(4.d0/27.d0)*proj_zona(iwan,isppol)**(1.5d0)&
3125 &       * (1.d0 - 2.d0*proj_zona(iwan,isppol)*pawrad(itypat)%rad(1:mesh_size)/3.d0 &
3126 &       + 2.d0*proj_zona(iwan,isppol)**2 *pawrad(itypat)%rad(1:mesh_size)**2/27.d0)&
3127 &       * exp(-proj_zona(iwan,isppol) * pawrad(itypat)%rad(1:mesh_size)/3.d0)
3128 !
3129 !      ==4: gaussian function
3130 !      f(x)=\exp(-1/4(x/aa)**2)
3131 !
3132        if(proj_radial(iwan,isppol)==4) then
3133          aa=1.d0/proj_zona(iwan,isppol)
3134          aux(1:mesh_size)= exp(-0.25d0*(pawrad(itypat)%rad(1:mesh_size)*aa)**2)
3135        end if
3136 !
3137 !      Normalize aux
3138        aux(:)=aux(:)/norm(iwan,isppol)
3139 !
3140        do ln=1,basis_size
3141          if(just_augmentation(iwan,isppol)) then
3142 !
3143 !          just augmentation region contribution
3144 !          In this case there is no need to use \tphi
3145 !          ff= \int R_wan(r) (R_phi(ln;r)/r ) r^2 dr
3146 !
3147            ff(1:mesh_size)= aux(1:mesh_size) * pawtab(itypat)%phi(1:mesh_size,ln) &
3148 &           * pawrad(itypat)%rad(1:mesh_size)
3149          else
3150 !          Inside sphere contribution = \phi - \tphi
3151 !          ff= \int R_wan(r) (R_phi(ln;r)/r - R_tphi(ln;r)/r) r^2 dr
3152            ff(1:mesh_size)= aux(1:mesh_size) * (pawtab(itypat)%phi(1:mesh_size,ln)-pawtab(itypat)%tphi(1:mesh_size,ln)) &
3153 &           * pawrad(itypat)%rad(1:mesh_size)
3154          end if
3155 !
3156 !        Integration with simpson routine
3157 !
3158          call simp_gen(int_rad(ln),ff,pawrad(itypat))
3159 !        do ii=1,mesh_size
3160 !        unit_ln=400+ln
3161 !        if( iwan==1 ) write(unit_ln,*)pawrad(itypat)%rad(ii),ff(ii),int_rad(ln)
3162 !        end do
3163        end do !ln
3164        ABI_DEALLOCATE(ff)
3165        ABI_DEALLOCATE(aux)
3166 !
3167 !      Now integrate the angular part
3168 !      Cycle on i indices
3169 !
3170 !      prod_real=0.d0
3171        do ilmn=1, lmn_size
3172          ll=Psps%indlmn(1,ilmn,itypat)
3173          mm=Psps%indlmn(2,ilmn,itypat)
3174          nn=Psps%indlmn(3,ilmn,itypat)
3175          lm=Psps%indlmn(4,ilmn,itypat)
3176          ln=Psps%indlmn(5,ilmn,itypat)
3177 !        write(std_out,*)'ll ',ll,' mm ',mm,'nn',nn,"lm",lm,"ln",ln
3178 !
3179 !        Get wannier factor for that lm component
3180          if(lm <=lmax2(isppol)) then
3181            wan_lm_fac=ylmr_fac(lm,iwan,isppol)
3182 !          Make delta product
3183 !          Here we integrate the angular part
3184 !          Since the integral of the product of two spherical harmonics
3185 !          is a delta function
3186            if( abs(wan_lm_fac) > 0.0d0) then
3187 !            write(std_out,*) 'll',ll,'mm',mm,'lm',lm,'ln',ln,'factor',wan_lm_fac !lm index for wannier function
3188 !
3189 !            Calculate Amn_paw, now that the radial and angular integrations are done
3190 !
3191              prod=cmplx(0.d0,0.d0)
3192              do ikpt=1,nkpt
3193                jband=0
3194                do iband=1,nband(ikpt)
3195                  if(band_in(iband,isppol)) then
3196                    jband=jband+1
3197 
3198                    prod_real= cprj(iatom,index(iband,ikpt,isppol))%cp(1,ilmn) * int_rad(ln) * wan_lm_fac
3199                    prod_imag= cprj(iatom,index(iband,ikpt,isppol))%cp(2,ilmn) * int_rad(ln) * wan_lm_fac
3200                    prod=cmplx(prod_real,prod_imag)
3201 
3202                    A_paw(jband,iwan,ikpt,isppol)=A_paw(jband,iwan,ikpt,isppol)+prod
3203                  end if !band_in
3204                end do !iband
3205              end do !ikpt
3206 !
3207            end if !lm<=lmax2
3208          end if  ! abs(wan_lm_fac) > 0.0d0
3209        end do !ilmn=1, lmn_size
3210        ABI_DEALLOCATE(int_rad)
3211      end do !iatom
3212    end do !iwan
3213  end do !isppol
3214 
3215 !Deallocate quantities
3216  ABI_DEALLOCATE(ylmr_fac)
3217 
3218 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=informations 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
  spin = just used for nsppol>1 ; 0 both, 1 just spin up, 2 just spin down

OUTPUT

  cm1(2,mband,mband,nntot,nkpt,nsppol): overlap <u_(nk1)|u_(mk1+b)>.

SIDE EFFECTS

  (only writing, printing)

NOTES

PARENTS

      mlwfovlp

CHILDREN

      wrtout,xmpi_barrier,xmpi_sum

SOURCE

1984 subroutine mlwfovlp_pw(cg,cm1,g1,iwav,kg,mband,mkmem,mpi_enreg,mpw,nfft,ngfft,nkpt,nntot,&
1985 &  npwarr,nspinor,nsppol,ovikp,seed_name,spin)
1986 
1987 
1988 !This section has been created automatically by the script Abilint (TD).
1989 !Do not modify the following lines by hand.
1990 #undef ABI_FUNC
1991 #define ABI_FUNC 'mlwfovlp_pw'
1992 !End of the abilint section
1993 
1994  implicit none
1995 
1996 !Arguments ------------------------------------
1997 !scalars
1998  integer,intent(in) :: mband,mkmem,mpw,nfft,nkpt,nntot
1999  integer,intent(in) :: nspinor,nsppol,spin
2000  character(len=fnlen) ::  seed_name  !seed names of files containing cg info used in case of MPI
2001  type(MPI_type),intent(in) :: mpi_enreg
2002 !arrays
2003  integer,intent(in) :: g1(3,nkpt,nntot),kg(3,mpw*mkmem),ngfft(18),npwarr(nkpt)
2004  integer,intent(in) :: iwav(mband,nkpt,nsppol)
2005  integer,intent(in) :: ovikp(nkpt,nntot)
2006  real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol)
2007  real(dp),intent(out) :: cm1(2,mband,mband,nntot,nkpt,nsppol)
2008 
2009 !Local variables-------------------------------
2010 !scalars
2011  integer :: iband1,iband2,ierr,ig,ig1,ig1b,ig2,ig2b
2012  integer :: ig3,ig3b,igk1,igk2,igks1,igks2,ii,ikg,ikpt,ikpt1,ikpt2,imntot,index,intot,ios,ipw
2013  integer :: ispinor,isppol,iunit,me,n1,n2,n3,npoint,npoint2,npw_k,npw_k2
2014  integer :: nprocs,spaceComm
2015  integer,allocatable::indpwk(:,:),kg_k(:,:)
2016  integer,allocatable :: invpwk(:,:)
2017  character(len=500) :: message
2018  character(len=fnlen) ::  cg_file  !file containing cg info used in case of MPI
2019  logical::lfile
2020  real(dp),allocatable :: cg_read(:,:) !to be used in case of MPI
2021 
2022 !************************************************************************
2023 
2024  write(message, '(a,a)' ) ch10,&
2025 & '** mlwfovlp_pw : compute pw part of overlap'
2026  call wrtout(std_out,  message,'COLL')
2027 
2028 !initialize flags
2029  lfile=.false.
2030 !mpi initialization
2031  spaceComm=MPI_enreg%comm_cell
2032  nprocs=xmpi_comm_size(spaceComm)
2033  me=MPI_enreg%me_kpt
2034 
2035  if(nprocs>1) then
2036    ABI_ALLOCATE(cg_read,(2,nspinor*mpw*mband))
2037  end if
2038 
2039 
2040 !****************compute intermediate quantities  (index, shifts) ******
2041 !------------compute index for g points--------------------------------
2042 !ig is a plane waves which belongs to the sphere ecut for ikpt (they
2043 !are npwarr(ikpt))
2044 !npoint is the position in the grid of planes waves
2045 !(they are nfft)
2046 !indpwk is a application ig-> npoint
2047 !invpwk is not an application (some npoint have no ig corresponding)
2048 !cg are ordered with npw_k !
2049 !----------------------------------------------------------------------
2050 !------------compute index for g points--------------------------------
2051 !----------------------------------------------------------------------
2052  write(message, '(a,a)' ) ch10,&
2053 & '   first compute index for g-points'
2054  call wrtout(std_out,  message,'COLL')
2055 !
2056 !Allocations
2057  ABI_ALLOCATE(kg_k,(3,mpw))
2058  ABI_ALLOCATE(indpwk,(nkpt,mpw))
2059  ABI_ALLOCATE(invpwk,(nkpt,nfft))
2060 !
2061  n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
2062  invpwk=0
2063  indpwk=0
2064  kg_k=0
2065  do isppol=1,1  !invpwk is not spin dependent
2066 !  so we just do it once
2067    ikg=0
2068    do ikpt=1,nkpt
2069 !
2070 !    MPI:cycle over k-points not treated by this node
2071 !
2072      if ( ABS(MPI_enreg%proc_distrb(ikpt,1,isppol)-me)  /=0) CYCLE
2073 
2074 !
2075 !    write(std_out,*)'me',me,'ikpt',ikpt,'isppol',isppol
2076      do npoint=1,nfft
2077        if(invpwk(ikpt,npoint)/=0 )then
2078          write(std_out,*) "error0 , invpwk is overwritten"
2079          write(std_out,*) ikpt,npoint
2080          MSG_ERROR("Aborting now")
2081        end if
2082      end do
2083      npw_k=npwarr(ikpt)
2084 !    write(std_out,*) ikpt,npw_k,nfft
2085      kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
2086      do ig=1,npw_k
2087        if(ig.gt.mpw) then
2088          write(std_out,*)"error ig",ig,"greater than mpw ",mpw
2089          MSG_ERROR("Aborting now")
2090        end if
2091        if(indpwk(ikpt,ig)/=0) then
2092          write(std_out,*) "error, indpwk is overwritten"
2093          write(std_out,*) ikpt,ig,indpwk(ikpt,ig)
2094          MSG_ERROR("Aborting now")
2095        end if
2096        ig1=modulo(kg_k(1,ig),n1)
2097        ig2=modulo(kg_k(2,ig),n2)
2098        ig3=modulo(kg_k(3,ig),n3)
2099        indpwk(ikpt,ig)=ig1+1+n1*(ig2+n2*ig3)
2100        npoint=indpwk(ikpt,ig)
2101        if(npoint.gt.nfft) then
2102          MSG_ERROR("error npoint")
2103        end if
2104 !      write(std_out,*) ikpt,ig,npoint,invpwk(ikpt,npoint)
2105        if(invpwk(ikpt,npoint)/=0) then
2106          write(std_out,*) "error, invpwk is overwritten"
2107          write(std_out,*) ikpt,ig,npoint,invpwk(ikpt,npoint)
2108          MSG_ERROR("Aborting now")
2109        end if
2110        invpwk(ikpt,npoint)=ig
2111 !      write(std_out,*)'ikpt,npoint,invpwk',ikpt,npoint,invpwk(ikpt,npoint)
2112 !      if(ikpt.eq.1) write(std_out,*) "ig npoint",ig, npoint
2113 !      write(std_out,*) "ikpt ig npoint",ikpt,ig, npoint
2114      end do
2115      ikg=ikg+npw_k
2116 
2117    end do !ikpt
2118  end do !isppol
2119 !write(std_out,*) "index for g points has been computed"
2120 
2121  call xmpi_barrier(spaceComm)
2122  call xmpi_sum(invpwk,spaceComm,ierr)
2123 
2124 !----------------------------------------------------------------------
2125 !------------test invpwk-----------------------------------------------
2126 !----------------------------------------------------------------------
2127 !write(std_out,*) "TEST INVPWK"
2128 !ikpt=3
2129 !isppol=1
2130 !do ig=1,npwarr(ikpt)
2131 !npoint=indpwk(ikpt,ig)
2132 !write(std_out,*) "ig npoint    ",ig, npoint
2133 !write(std_out,*) "ig npoint inv",invpwk(ikpt,npoint),npoint
2134 !end do
2135 !do ig3=1,n3
2136 !do ig2=1,n2
2137 !do ig1=1,n1
2138 !npoint=ig1+(ig2-1)*n1+(ig3-1)*n2*n1
2139 !ig=invpwk(ikpt,npoint)
2140 !!   if(ig/=0)  write(std_out,*) "ig npoint",ig, npoint
2141 !end do
2142 !end do
2143 !end do
2144 
2145 !Deallocate unused variables
2146  ABI_DEALLOCATE(kg_k)
2147  ABI_DEALLOCATE(indpwk)
2148 
2149 !***********************************************************************
2150 !**calculate overlap M_{mn}(k,b)=<\Psi_{k,m}|e^{-ibr}|\Psi_{k+b,n}>*****
2151 !***********************************************************************
2152  write(message, '(a,a)' ) ch10,&
2153 & '   mlwfovlp_pw : compute overlaps '
2154  call wrtout(std_out,  message,'COLL')
2155  write(message, '(a,a)' ) ch10,&
2156 & "     nkpt  nntot  mband "
2157  call wrtout(std_out,  message,'COLL')
2158  write(message, '(i6,2x,i6,2x,i6,2x,i6)' ) &
2159 & nkpt,nntot,mband
2160  call wrtout(std_out,  message,'COLL')
2161  cm1=zero
2162  write(message, '(a)' )  '  '
2163  call wrtout(std_out,  message,'COLL')
2164  do isppol=1,nsppol
2165    if(spin.ne.0 .and. spin.ne.isppol) cycle
2166    imntot=0
2167    do ikpt1=1,nkpt
2168 !
2169 !    MPI:cycle over k-points not treated by this node
2170 !
2171      if ( ABS(MPI_enreg%proc_distrb(ikpt1,1,isppol)-me)  /=0) CYCLE
2172 !
2173      write(message, '(a,i6,a,i6,a,i6)' ) &
2174 &     '     Processor',me,' computes k-point',ikpt1,' and spin=',isppol
2175      call wrtout(std_out,  message,'COLL')
2176 !    write(std_out,*)trim(message)
2177 
2178      do intot=1,nntot
2179        lfile=.false. !flag to know if this kpt will be read from a file, see below
2180        imntot=imntot+1
2181        ikpt2= ovikp(ikpt1,intot)
2182 !      write(std_out,*)'me',me,'ikpt1',ikpt1,'ikpt2',ikpt2,'intot',intot,'isppol',isppol
2183 
2184 !
2185 !      MPI: if ikpt2 not found in this processor then
2186 !      read info from an unformatted file
2187 !
2188        if ( ABS(MPI_enreg%proc_distrb(ikpt2,1,isppol)-me)  /=0) then
2189          lfile=.true.
2190          write(cg_file,'(a,I5.5,".",I1)') trim(seed_name),ikpt2,isppol
2191          iunit=1000+ikpt2+ikpt2*(isppol-1)
2192          npw_k2=npwarr(ikpt2)
2193          open (unit=iunit, file=cg_file,form='unformatted',status='old',iostat=ios)
2194          if(ios /= 0) then
2195            write(message,*) " mlwfovlp_pw: file",trim(cg_file), "not found"
2196            MSG_ERROR(message)
2197          end if
2198 !
2199          do iband2=1,mband
2200            do ipw=1,npw_k2*nspinor
2201              index=ipw+(iband2-1)*npw_k2*nspinor
2202              read(iunit) (cg_read(ii,index),ii=1,2)
2203 !            if(me==0 .and. ikpt2==4)write(300,*)'ipw,iband2,index',ipw,iband2,index,cg_read(:,index)
2204 !            if(me==1 .and. ikpt2==4)write(301,*)'ipw,iband2,index',ipw,iband2,index,cg_read(:,index)
2205            end do
2206          end do
2207          close(iunit)
2208        end if
2209 !
2210        npw_k=npwarr(ikpt1)
2211        npw_k2=npwarr(ikpt2)
2212        do ig3=1,n3
2213          do ig2=1,n2
2214            do ig1=1,n1
2215 !            write(std_out,*) isppol,ikpt1,iband1,iband2,intot
2216              npoint=ig1+(ig2-1)*n1+(ig3-1)*n2*n1
2217              if(npoint.gt.nfft) then
2218                write(std_out,*) "error npoint  b"
2219                MSG_ERROR("Aborting now")
2220              end if
2221              ig1b=ig1+g1(1,ikpt1,intot)
2222              ig2b=ig2+g1(2,ikpt1,intot)
2223              ig3b=ig3+g1(3,ikpt1,intot)
2224 !            write(std_out,*) ig1,ig2,ig3
2225 !            write(std_out,*) ig1b,ig2b,ig3b
2226              if(ig1b.lt.1) ig1b=ig1b+n1
2227              if(ig2b.lt.1) ig2b=ig2b+n2
2228              if(ig3b.lt.1) ig3b=ig3b+n3
2229              if(ig1b.gt.n1) ig1b=ig1b-n1
2230              if(ig2b.gt.n2) ig2b=ig2b-n2
2231              if(ig3b.gt.n3) ig3b=ig3b-n3
2232              npoint2=ig1b+(ig2b-1)*n1+(ig3b-1)*n2*n1
2233              if(npoint2.gt.nfft) then
2234                write(std_out,*)"error npoint  c"
2235                MSG_ERROR("Aborting now")
2236              end if
2237              igk1=invpwk(ikpt1,npoint)
2238              igk2=invpwk(ikpt2,npoint2)
2239 
2240 !            if(intot==10) write(std_out,*)'Before igk1 and igk2',ikpt1,ikpt2,isppol
2241 
2242              if(igk1/=0.and.igk2/=0) then
2243                do iband2=1,mband
2244                  do iband1=1,mband
2245                    do ispinor=1,nspinor
2246 !                    igks1= (igk1*nspinor)-(nspinor-ispinor)
2247 !                    igks2= (igk2*nspinor)-(nspinor-ispinor)
2248                      igks1= igk1+ (ispinor-1)*npw_k
2249                      igks2= igk2+ (ispinor-1)*npw_k2
2250 
2251 !                    Here the igks is to include the spinor component missing in igk
2252                      if(lfile) index=igks2+npw_k2*nspinor*(iband2-1) !In case of MPI, see below
2253 !
2254 !                    If MPI sometimes the info was read from an unformatted file
2255 !                    If that is the case lfile==.true.
2256 !
2257                      if(lfile) then
2258                        cm1(1,iband1,iband2,intot,ikpt1,isppol)=cm1(1,iband1,iband2,intot,ikpt1,isppol)+ &
2259 &                       cg(1,igks1+iwav(iband1,ikpt1,isppol))*cg_read(1,index)&
2260 &                       + cg(2,igks1+iwav(iband1,ikpt1,isppol))*cg_read(2,index)
2261                        cm1(2,iband1,iband2,intot,ikpt1,isppol)=cm1(2,iband1,iband2,intot,ikpt1,isppol)+ &
2262 &                       cg(1,igks1+iwav(iband1,ikpt1,isppol))*cg_read(2,index)&
2263 &                       - cg(2,igks1+iwav(iband1,ikpt1,isppol))*cg_read(1,index)
2264 !
2265                      else
2266                        cm1(1,iband1,iband2,intot,ikpt1,isppol)=cm1(1,iband1,iband2,intot,ikpt1,isppol)+ &
2267 &                       cg(1,igks1+iwav(iband1,ikpt1,isppol))*cg(1,igks2+iwav(iband2,ikpt2,isppol))&
2268 &                       + cg(2,igks1+iwav(iband1,ikpt1,isppol))*cg(2,igks2+iwav(iband2,ikpt2,isppol))
2269                        cm1(2,iband1,iband2,intot,ikpt1,isppol)=cm1(2,iband1,iband2,intot,ikpt1,isppol)+ &
2270 &                       cg(1,igks1+iwav(iband1,ikpt1,isppol))*cg(2,igks2+iwav(iband2,ikpt2,isppol))&
2271 &                       - cg(2,igks1+iwav(iband1,ikpt1,isppol))*cg(1,igks2+iwav(iband2,ikpt2,isppol))
2272                      end if
2273                    end do !ispinor
2274                  end do ! iband1
2275                end do ! iband2
2276              end if
2277            end do
2278          end do
2279        end do
2280      end do ! intot
2281    end do ! ikpt1
2282  end do ! isppol
2283 !
2284 !Deallocations
2285 !
2286  ABI_DEALLOCATE(invpwk)
2287  if(nprocs>1)  then
2288    ABI_DEALLOCATE(cg_read)
2289  end if
2290 
2291  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.

PARENTS

      mlwfovlp_proj

CHILDREN

      besjm,simpson_int

SOURCE

3258 subroutine mlwfovlp_radial(alpha,lmax,lmax2,radial,rvalue,xx)
3259 
3260 
3261 !This section has been created automatically by the script Abilint (TD).
3262 !Do not modify the following lines by hand.
3263 #undef ABI_FUNC
3264 #define ABI_FUNC 'mlwfovlp_radial'
3265 !End of the abilint section
3266 
3267  implicit none
3268 
3269 !Arguments ------------------------------------
3270 !scalars
3271  integer,intent(in) :: lmax,lmax2,rvalue
3272  real(dp),intent(in) :: alpha,xx
3273 !arrays
3274  real(dp),intent(out) :: radial(lmax2)
3275 
3276 !Local variables
3277 !scalars
3278  integer :: ir,ll,lm,mesh,mm
3279  real(dp),parameter :: dx=0.015d0,rmax=10.d0,xmin=0.d0
3280  real(dp) :: aa,ftmp,gauss,rtmp,x
3281  character(len=500) :: message
3282 !arrays
3283  real(dp),save :: dblefact(4)=(/1_dp,3_dp,15_dp,105_dp/)
3284  real(dp),allocatable :: aux(:),bes(:),cosr(:),func_r(:),r(:),rad_int(:)
3285  real(dp),allocatable :: sinr(:)
3286 
3287 ! *************************************************************************
3288 
3289 !Radial functions in the form of hydrogenic orbitals as defined in the
3290 !wannier90 manual.
3291  if(( rvalue > 0 ).and.(rvalue < 4)) then
3292 
3293 !  mesh
3294    mesh= nint((rmax - xmin ) / dx + 1)
3295    ABI_ALLOCATE( bes,(mesh))
3296    ABI_ALLOCATE(func_r,(mesh))
3297    ABI_ALLOCATE(r,(mesh))
3298    ABI_ALLOCATE(rad_int,(mesh))
3299    ABI_ALLOCATE( aux,(mesh))
3300    ABI_ALLOCATE(cosr,(mesh))
3301    ABI_ALLOCATE(sinr,(mesh))
3302    do ir=1, mesh
3303      x=xmin+DBLE(ir-1)*dx
3304      r(ir)=x
3305    end do   !ir
3306 
3307 !  radial functions shown in table 3.3 of wannier90 manual
3308    if (rvalue==1) func_r(:) = 2.d0 * alpha**(3.d0/2.d0) * exp(-alpha*r(:))
3309    if (rvalue==2) func_r(:) = 1.d0/(2.d0*sqrt(2.d0))*alpha**(3.d0/2.d0) *&
3310 &   (2.d0 - alpha*r(:))*exp(-alpha*r(:)/2.d0)
3311    if (rvalue==3) func_r(:) = sqrt(4.d0/27.d0)*alpha**(3.d0/2.d0)&
3312 &   * (1.d0 - 2.d0*alpha*r(:)/3.d0 + 2.d0*alpha**2*r(:)**2/27.d0)&
3313 &   * exp(-alpha * r(:)/3.d0)
3314 
3315 !  compute spherical bessel functions
3316    cosr(:)=cos(xx*r(:))
3317    sinr(:)=sin(xx*r(:))
3318    lm=0
3319    do ll=0,lmax
3320      call besjm(xx,bes,cosr,ll,mesh,sinr,r)
3321      aux(:)=bes(:)*func_r(:)*r(:)
3322 !    do ir=1,mesh
3323 !    write(310,*) r(ir),bes(ir)
3324 !    end do
3325      call simpson_int(mesh,dx,aux,rad_int)
3326      rtmp=rad_int(mesh)/mesh
3327      do mm=-ll,ll
3328        lm=lm+1
3329        radial(lm)=rtmp
3330      end do !mm
3331    end do !ll
3332    ABI_DEALLOCATE(bes)
3333    ABI_DEALLOCATE(func_r)
3334    ABI_DEALLOCATE(r)
3335    ABI_DEALLOCATE(aux)
3336    ABI_DEALLOCATE(rad_int)
3337    ABI_DEALLOCATE(cosr)
3338    ABI_DEALLOCATE(sinr)
3339 
3340 !  Radial part in the form of Gaussian functions of a given width
3341 !  Taken by code of made by drh.
3342  elseif ( rvalue == 4) then
3343    aa=1._dp/alpha
3344    gauss=exp(-0.25_dp*(aa*xx)**2)
3345    lm=0
3346    do ll=0,lmax
3347      ftmp=(0.5_dp*pi)**(0.25_dp)*aa*sqrt(aa/dblefact(ll+1))*(aa*xx)**ll*gauss
3348      do mm=-ll,ll
3349        lm=lm+1
3350        radial(lm)=ftmp
3351      end do
3352    end do
3353  else ! rvalue < 0 of rvalue > 4
3354    write(message,'(a,i6,5a)')&
3355 &   '  Radial function r=',rvalue,ch10,&
3356 &   '  is not defined',ch10,&
3357 &   '  Modify .win file',ch10
3358    MSG_BUG(message)
3359  end if !rvalue
3360 
3361 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

PARENTS

      mlwfovlp

CHILDREN

      wrtout

SOURCE

1484 subroutine mlwfovlp_seedname(fname_w90,filew90_win,filew90_wout,filew90_amn,&
1485 & filew90_ramn,filew90_mmn,filew90_eig,nsppol,seed_name)
1486 
1487 
1488 !This section has been created automatically by the script Abilint (TD).
1489 !Do not modify the following lines by hand.
1490 #undef ABI_FUNC
1491 #define ABI_FUNC 'mlwfovlp_seedname'
1492 !End of the abilint section
1493 
1494  implicit none
1495 
1496 !Arguments ------------------------------------
1497  integer,intent(in) :: nsppol
1498  character(len=fnlen),intent(out) :: filew90_win(nsppol),filew90_wout(nsppol),filew90_amn(nsppol),filew90_ramn(nsppol)
1499  character(len=fnlen),intent(out) :: filew90_mmn(nsppol),filew90_eig(nsppol),seed_name(nsppol)
1500  character(len=fnlen),intent(in) :: fname_w90
1501 
1502 !Local variables-------------------------------
1503  integer::isppol
1504  character(len=fnlen) :: test_win1,test_win2,test_win3
1505  logical :: lfile
1506  character(len=2000) :: message   
1507  character(len=10)::postfix
1508 ! *************************************************************************
1509 
1510  seed_name(:)=trim(fname_w90)
1511  do isppol=1,nsppol
1512    if(nsppol==1)postfix='.win'
1513    if(nsppol==2 .and. isppol==1)postfix='_up.win'
1514    if(nsppol==2 .and. isppol==2)postfix='_down.win'
1515 !
1516    filew90_win(isppol)=trim(seed_name(isppol))//trim(postfix)
1517    test_win1=filew90_win(isppol)
1518    inquire(file=filew90_win(isppol),exist=lfile)
1519 
1520    if(.not.lfile) then
1521      seed_name(isppol)='wannier90'
1522      filew90_win(isppol)=trim(seed_name(isppol))//trim(postfix)
1523      test_win2=filew90_win(isppol)
1524      inquire(file=filew90_win(isppol),exist=lfile)
1525    end if
1526 
1527    if(.not.lfile) then
1528      seed_name(isppol)='w90'
1529      filew90_win=trim(seed_name(isppol))//trim(postfix)
1530      test_win3=filew90_win(isppol)
1531      inquire(file=filew90_win(isppol),exist=lfile)
1532    end if
1533 
1534    if(.not.lfile) then
1535      write(message,'(17a)')ch10,&
1536 &     ' mlwfovlp_seedname : ERROR - ',ch10,&
1537 &     ' wannier90 interface needs one of the following files:',ch10,&
1538 &     '      ',trim(test_win1),ch10,&
1539 &     '      ',trim(test_win2),ch10,&
1540 &     '      ',trim(test_win3),ch10,&
1541 &     ' Action: read wannier90 tutorial and/or user manual',ch10,&
1542 &     '  and supply proper *.win file'
1543      MSG_ERROR(message)
1544    end if
1545  end do !isppol
1546 
1547 
1548 !Files having different names for
1549 !different spin polarizations
1550  if(nsppol==1) then
1551    filew90_win(1) =trim(seed_name(1))//'.win'
1552    filew90_wout(1)=trim(seed_name(1))//'.wout'
1553    filew90_ramn(1)=trim(seed_name(1))//'random.amn'
1554    filew90_amn(1) =trim(seed_name(1))//'.amn'
1555    filew90_mmn(1) =trim(seed_name(1))//'.mmn'
1556    filew90_eig(1) =trim(seed_name(1))//'.eig'
1557  elseif(nsppol==2) then
1558    filew90_win(1) =trim(seed_name(1))//'_up.win'
1559    filew90_win(2) =trim(seed_name(2))//'_down.win'
1560 !
1561    filew90_wout(1)=trim(seed_name(1))//'_up.wout'
1562    filew90_wout(2)=trim(seed_name(2))//'_down.wout'
1563 !
1564    filew90_ramn(1)=trim(seed_name(1))//'random_up.amn'
1565    filew90_ramn(2)=trim(seed_name(2))//'random_down.amn'
1566 !
1567    filew90_amn(1)=trim(seed_name(1))//'_up.amn'
1568    filew90_amn(2)=trim(seed_name(2))//'_down.amn'
1569 !
1570    filew90_mmn(1)=trim(seed_name(1))//'_up.mmn'
1571    filew90_mmn(2)=trim(seed_name(2))//'_down.mmn'
1572 !
1573    filew90_eig(1)=trim(seed_name(1))//'_up.eig'
1574    filew90_eig(2)=trim(seed_name(2))//'_down.eig'
1575  end if
1576 !change also seed_name for nsppol=2
1577  if(nsppol==2) then
1578    seed_name(1)=trim(seed_name(1))//'_up'
1579    seed_name(2)=trim(seed_name(2))//'_down'
1580  end if
1581 !End file-name section
1582 
1583  write(message, '(a,a)' ) ch10,&
1584 & '---------------------------------------------------------------'
1585  call wrtout(ab_out,message,'COLL')
1586  call wrtout(std_out,  message,'COLL')
1587  write(message, '(5a)' ) ch10,&
1588 & '  Calculation of overlap and call to wannier90 library ',ch10,&
1589 & '  to obtain maximally localized wannier functions ',ch10
1590 
1591  call wrtout(std_out,  message,'COLL')
1592  call wrtout(ab_out,message,'COLL')
1593 
1594  if(nsppol==1) then
1595    write(message, '(23a)' ) &
1596 &   '  - ',trim(filew90_win(1)),' is a mandatory secondary input',ch10,&
1597 &   '  - ',trim(filew90_wout(1)),' is the output for the library',ch10,&
1598 &   '  - ',trim(filew90_ramn(1)),' contains random projections',ch10,&
1599 &   '  - ',trim(filew90_amn(1)),' contains projections',ch10,&
1600 &   '  - ',trim(filew90_mmn(1)),' contains the overlap',ch10,&
1601 &   '  - ',trim(filew90_eig(1)),' contains the eigenvalues'
1602  elseif(nsppol==2) then
1603    write(message, '(41a)' ) &
1604 &   '  - ',trim(filew90_win(1)),&
1605 &   ' and ',trim(filew90_win(2)),ch10,'are mandatory secondary input',ch10,&
1606 &   '  - ',trim(filew90_wout(1)),&
1607 &   ' and ',trim(filew90_wout(2)),ch10,' are the output for the library',ch10,&
1608 &   '  - ',trim(filew90_ramn(1)),&
1609 &   ' and ',trim(filew90_ramn(2)),ch10,' contain random projections',ch10,&
1610 &   '  - ',trim(filew90_amn(1)),&
1611 &   ' and ',trim(filew90_amn(2)),ch10,' contain projections',ch10,&
1612 &   '  - ',trim(filew90_mmn(1)),&
1613 &   ' and ',trim(filew90_mmn(2)),ch10,' contain the overlap',ch10,&
1614 &   '  - ',trim(filew90_eig(1)),&
1615 &   ' and ',trim(filew90_eig(2)),ch10,' contain the eigenvalues'
1616  end if
1617  call wrtout(std_out,  message,'COLL')
1618  call wrtout(ab_out,message,'COLL')
1619 
1620  write(message, '(a,a)' ) ch10,&
1621 & '---------------------------------------------------------------'
1622  call wrtout(ab_out,message,'COLL')
1623  call wrtout(std_out,  message,'COLL')
1624 
1625 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
  spin = just used for nsppol>1 ; 0 both, 1 just spin up, 2 just spin down
  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

PARENTS

      mlwfovlp

CHILDREN

      atomdata_from_znucl,wannier_setup,wrtout

SOURCE

1680  subroutine mlwfovlp_setup(atom_symbols,band_in,dtset,filew90_win,gamma_only,&
1681 & g1,lwanniersetup,mband,natom,nband_inc,nkpt,&
1682 & nntot,num_bands,num_nnmax,nsppol,nwan,ovikp,&
1683 & proj_l,proj_m,proj_radial,proj_site,proj_s_loc,proj_s_qaxis_loc,proj_x,proj_z,proj_zona,&
1684 & real_lattice,recip_lattice,rprimd,seed_name,spin,spinors,xcart,xred)
1685 
1686 
1687 !This section has been created automatically by the script Abilint (TD).
1688 !Do not modify the following lines by hand.
1689 #undef ABI_FUNC
1690 #define ABI_FUNC 'mlwfovlp_setup'
1691 !End of the abilint section
1692 
1693  implicit none
1694 
1695 !Arguments---------------------------
1696 ! scalars
1697 !scalars
1698  integer,intent(in) :: lwanniersetup,mband,natom,nkpt,nsppol
1699  integer,intent(in) :: num_nnmax,spin
1700  integer,intent(out) :: nband_inc(nsppol),nntot,num_bands(nsppol),nwan(nsppol)
1701  logical,intent(in) :: gamma_only,spinors
1702  type(dataset_type),intent(in) :: dtset
1703 !arrays
1704  integer,intent(out) :: g1(3,nkpt,num_nnmax),ovikp(nkpt,num_nnmax)
1705  integer,intent(out) :: proj_l(mband,nsppol),proj_m(mband,nsppol),proj_radial(mband,nsppol)
1706  real(dp),intent(in) :: real_lattice(3,3)
1707  real(dp),intent(in) :: recip_lattice(3,3),rprimd(3,3),xred(3,natom)
1708  real(dp),intent(out) :: proj_site(3,mband,nsppol),proj_x(3,mband,nsppol),proj_z(3,mband,nsppol)
1709  real(dp),intent(out) :: proj_zona(mband,nsppol),xcart(3,natom)
1710  logical,intent(out) :: band_in(mband,nsppol)
1711  character(len=3),intent(out) :: atom_symbols(natom)
1712  character(len=fnlen),intent(in) :: seed_name(nsppol),filew90_win(nsppol)
1713 
1714  integer, optional, intent(out) :: proj_s_loc(mband)
1715  real(dp), optional, intent(out) :: proj_s_qaxis_loc(3,mband)
1716 
1717 !Local variables---------------------------
1718 !scalars
1719  integer :: iatom,icb,ikpt,ikpt1,intot,isppol,itypat,jj,mband_,unt
1720  real(dp) :: znucl1
1721  character(len=2) :: symbol
1722  character(len=500) :: message
1723  character(len=fnlen) :: filew90_nnkp
1724  type(atomdata_t) :: atom
1725 !arrays
1726  integer :: exclude_bands(mband,nsppol),ngkpt(3)
1727 
1728 ! *************************************************************************
1729 
1730 !^^^^^^^^^^^^^^^^read wannier90.nnkp^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1731  if(lwanniersetup==0) then  !this part is not coded for nsppol>1
1732    isppol=1
1733    filew90_nnkp=trim(seed_name(isppol))//'.nnkp'
1734    if (open_file(filew90_nnkp,message,newunit=unt,form='formatted',status='old') /= 0) then
1735      MSG_ERROR(message)
1736    end if
1737    read(unt,*)
1738    read(unt,*) nntot , mband_, nwan(1)
1739    write(message, '(a,a,i6,i6,i6)' )ch10,&
1740 &   ' mlwfovlp_setup nntot,mband,nwan ', nntot,mband_,nwan(1)
1741    call wrtout(std_out,message,'COLL')
1742    if(mband_.ne.mband) then
1743      write(message, '(4a)' )&
1744 &     'mband_ is not equal to mband ',ch10,&
1745 &     'Action: check ',trim(filew90_nnkp)
1746      MSG_ERROR(message)
1747    end if
1748    if(nwan(1)>mband) then
1749      write(message, '(4a)' )&
1750 &     'nwan > mband ',ch10,&
1751 &     'Action: check ',trim(filew90_nnkp)
1752      MSG_ERROR(message)
1753    end if
1754    if(nwan(1)==0) then
1755      write(message, '(4a)' )&
1756 &     'nwan = 0 ',ch10,&
1757 &     'Action: check ',trim(filew90_nnkp)
1758      MSG_ERROR(message)
1759    end if
1760    do ikpt=1,nkpt
1761      do intot=1,nntot
1762 !      ikpt1: k point  (ikpt=ikpt1)
1763 !      ovikp(intot,ikpt): neighbour number intot for ikpt
1764 !      g1(1:3,intot,ikpt): non reciprocal space vector between the 2 k-points
1765        read(unt,*)  &
1766 &       ikpt1,ovikp(ikpt,intot),(g1(jj,ikpt,intot),jj=1,3)
1767        if(ikpt1.ne.ikpt) write(std_out,*) "warning: ikpt1 .ne ikpt : ?"
1768      end do
1769    end do
1770    close(unt)
1771    write(message, '(3a)' )ch10,&
1772 &   trim(filew90_nnkp),'wannier90.nnkp has been read !'
1773    call wrtout(std_out,message,'COLL')
1774 
1775    message = ' exclude bands is not given in this case (not implemented) '
1776    MSG_ERROR(message)
1777 
1778 !  ^^^^^^^^^^^^^^^^^^^^^^^ call wannier_setup begin^^^^^^^^^^^^^^^^^^^^^^^^
1779  else if (lwanniersetup==1) then
1780    num_bands(:)=mband
1781 !  num_nnmax=12 !limit fixed for compact structure in wannier_setup.
1782    ovikp=0.d0
1783 !  "When nshiftk=1, kptrlatt is initialized as a diagonal (3x3) matrix, whose diagonal
1784 !  elements are the three values ngkpt(1:3)"
1785    ngkpt(1)=dtset%kptrlatt(1,1)
1786    ngkpt(2)=dtset%kptrlatt(2,2) !  have to verif kptrlatt is diagonal
1787    ngkpt(3)=dtset%kptrlatt(3,3)
1788    do iatom=1,natom
1789      itypat=dtset%typat(iatom)
1790      znucl1=dtset%znucl(itypat)
1791      call atomdata_from_znucl(atom, znucl1)
1792      symbol=trim(adjustl(atom%symbol))
1793 !    write(309,*) symbol
1794      atom_symbols(iatom)=symbol
1795      xcart(:,iatom)=rprimd(:,1)*xred(1,iatom)+&
1796 &     rprimd(:,2)*xred(2,iatom)+&
1797 &     rprimd(:,3)*xred(3,iatom)
1798    end do ! iatom
1799 !  write(std_out,*) xcart
1800 !  write(std_out,*) Bohr_Ang
1801 !  write(std_out,*) rprimd*Bohr_Ang
1802 !  write(std_out,*) seed_name
1803 !  write(std_out,*) ngkpt
1804 !  write(std_out,*) nkpt
1805 !  write(std_out,*) mband
1806 !  write(std_out,*) natom
1807 !  write(std_out,*) atom_symbols
1808    write(message, '(a,a)' )ch10,&
1809 &   '** mlwfovlp_setup:  call wannier90 library subroutine wannier_setup'
1810    call wrtout(std_out,message,'COLL')
1811 #if defined HAVE_WANNIER90
1812    nwan(:)=0
1813    num_bands(:)=0
1814    do isppol=1,nsppol
1815      if(spin.ne.0 .and. spin.ne.isppol) cycle
1816 #ifndef HAVE_WANNIER90_V1
1817 !WANNIER90_V2 has the 2 optional arguments
1818      if (present(proj_s_loc)) then
1819        call wannier_setup(seed_name(isppol),ngkpt,nkpt&            !input
1820 &      ,real_lattice,recip_lattice,dtset%kpt&                      !input
1821 &      ,mband,natom,atom_symbols,xcart*Bohr_Ang&                   !input
1822 &      ,gamma_only,spinors&                                        !input
1823 &      ,nntot,ovikp,g1,num_bands(isppol),nwan(isppol)&             !output
1824 &      ,proj_site(:,:,isppol),proj_l(:,isppol)&                    !output
1825 &      ,proj_m(:,isppol),proj_radial(:,isppol)&                    !output
1826 &      ,proj_z(:,:,isppol),proj_x(:,:,isppol)&                     !output
1827 &      ,proj_zona(:,isppol),exclude_bands(:,isppol)&               !output
1828 &      ,proj_s_loc,proj_s_qaxis_loc)                               !output
1829      else
1830 #else
1831 !WANNIER90_V1 or no proj_s_loc provided
1832        call wannier_setup(seed_name(isppol),ngkpt,nkpt&            !input
1833 &      ,real_lattice,recip_lattice,dtset%kpt&                      !input
1834 &      ,mband,natom,atom_symbols,xcart*Bohr_Ang&                   !input
1835 &      ,gamma_only,spinors&                                        !input
1836 &      ,nntot,ovikp,g1,num_bands(isppol),nwan(isppol)&             !output
1837 &      ,proj_site(:,:,isppol),proj_l(:,isppol)&                    !output
1838 &      ,proj_m(:,isppol),proj_radial(:,isppol)&                    !output
1839 &      ,proj_z(:,:,isppol),proj_x(:,:,isppol)&                     !output
1840 &      ,proj_zona(:,isppol),exclude_bands(:,isppol))               !output
1841 #endif
1842 #ifndef HAVE_WANNIER90_V1
1843      end if
1844 #endif
1845    end do !isppol
1846 #else
1847    ABI_UNUSED(gamma_only)
1848    ABI_UNUSED(real_lattice)
1849    ABI_UNUSED(recip_lattice)
1850    ABI_UNUSED(spinors)
1851 #endif
1852 !  do isppol=1,nsppol
1853 !  if(spin.ne.0 .and. spin.ne.isppol) cycle
1854 !  write(std_out,*)  "1", nntot,nwan(isppol)
1855 !  write(std_out,*)  "2",num_bands(isppol)  ! states on which wannier functions are computed
1856 !  write(std_out,*)  "3", proj_site(:,1:nwan(isppol),isppol)
1857 !  write(std_out,*)  "4",proj_l(1:nwan(isppol),isppol)
1858 !  write(std_out,*)  "5",proj_m(1:nwan(isppol),isppol)
1859 !  write(std_out,*)  "6", proj_radial(1:nwan(isppol),isppol)
1860 !  write(std_out,*)  "7", proj_z(:,1:nwan(isppol),isppol)
1861 !  write(std_out,*)  "8", proj_x(:,1:nwan(isppol),isppol)
1862 !  write(std_out,*)  "9",proj_zona(1:nwan(isppol),isppol)
1863 !  write(std_out,*)  "10",exclude_bands(:,isppol)
1864 !  end do!isppol
1865 !  testdebug:  ovikp(1,1)=1
1866 !  ^^^^^^^^^^^^^^^^^^^^^^^ end ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1867  end if  ! lwanniersetup
1868  do isppol=1,nsppol
1869    if(spin.ne.0 .and. spin.ne.isppol) cycle
1870    band_in(:,isppol)=.true.
1871    do icb=1,mband
1872      if(exclude_bands(icb,isppol).ne.0)  band_in(exclude_bands(icb,isppol),isppol)=.false.
1873    end do
1874    nband_inc(isppol)=0
1875    do icb=1, mband
1876      if (band_in(icb,isppol)) then
1877        nband_inc(isppol)=nband_inc(isppol)+1
1878      end if
1879    end do
1880  end do !isppol
1881  if(any(mband.gt.num_bands(:))) then
1882    write(message, '(a,a)' )ch10,&
1883 &   '   Following bands are excluded from the calculation of wannier functions:'
1884    call wrtout(std_out,message,'COLL')
1885 
1886    do isppol=1,nsppol
1887      if(spin.ne.0 .and. spin.ne.isppol) cycle
1888      if(nsppol==2) then
1889        write(message,'("For spin",i4)')isppol
1890 !      write(message,'(a,i)')'For spin=',isppol
1891        call wrtout(std_out,message,'COLL')
1892      end if !nsppol
1893      do jj=1,mband-num_bands(isppol),10
1894        write(message,'(10i7)') exclude_bands(jj:min(jj+9,mband-num_bands(isppol)),isppol)
1895        call wrtout(std_out,message,'COLL')
1896      end do
1897    end do !isppol
1898  end if
1899 
1900  do isppol=1,nsppol
1901    if(spin.ne.0 .and. spin.ne.isppol) cycle
1902    if(nsppol==2) then
1903      write(message,'("For spin",i4)')isppol
1904      call wrtout(std_out,message,'COLL')
1905    end if !nsppol
1906    write(message, '(a,i6,3a)' )ch10,&
1907 &   nwan(isppol),' wannier functions will be computed (see ',trim(filew90_win(isppol)),')'
1908    call wrtout(std_out,message,'COLL')
1909 !  write(std_out,*) exclude_bands(icb),band_in(icb)
1910 !  ^^^^^^^^^^^^^^^END OF READING
1911    write(message, '(a,i6,a)' )ch10,&
1912 &   num_bands(isppol),' bands will be used to extract wannier functions'
1913    call wrtout(std_out,message,'COLL')
1914    if(num_bands(isppol).lt.nwan(isppol)) then
1915      write(message, '(4a)' )&
1916 &     ' number of bands is lower than the number of wannier functions',ch10,&
1917 &     ' Action : check input file and ',trim(filew90_win(isppol))
1918      MSG_ERROR(message)
1919    else if (num_bands(isppol)==nwan(isppol)) then
1920      write(message, '(a,a,a,a)' )ch10,&
1921 &     '   Number of bands is equal to the number of wannier functions',ch10,&
1922 &     '   Disentanglement will not be necessary'
1923      call wrtout(std_out,message,'COLL')
1924    else if  (num_bands(isppol).gt.nwan(isppol)) then
1925      write(message, '(a,a,a,a)' )ch10,&
1926 &     '   Number of bands is larger than the number of wannier functions',ch10,&
1927 &     '   Disentanglement will be necessary'
1928      call wrtout(std_out,message,'COLL')
1929    end if
1930    write(message, '(2x,a,a,i3,1x,a)' )ch10,&
1931 &   '   Each k-point has', nntot,'neighbours'
1932    call wrtout(std_out,message,'COLL')
1933 
1934  end do !isppol
1935 
1936 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

PARENTS

      mlwfovlp_proj

CHILDREN

      rotmat,ylm_cmplx,zgesv

SOURCE

3409 subroutine mlwfovlp_ylmfac(ylmc_fac,lmax,lmax2,mband,nwan,proj_l,proj_m,proj_x,proj_z)
3410 
3411 
3412 !This section has been created automatically by the script Abilint (TD).
3413 !Do not modify the following lines by hand.
3414 #undef ABI_FUNC
3415 #define ABI_FUNC 'mlwfovlp_ylmfac'
3416 !End of the abilint section
3417 
3418  implicit none
3419 
3420 !Arguments ------------------------------------
3421  integer, intent(in):: lmax,lmax2,nwan,mband
3422 ! arrays
3423  integer,intent(in) :: proj_l(mband),proj_m(mband)
3424  real(dp),intent(in) :: proj_x(3,mband),proj_z(3,mband)
3425  complex(dp),intent(out)::ylmc_fac(lmax2,nwan)
3426 !
3427 !Local variables-------------------------------
3428 !
3429  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
3430  integer :: idum,ii,info,inversion_flag
3431  integer :: ir,iwan,jj,ll,lm,lmc,mm,mr
3432  real(dp):: onem,test
3433 ! arrays
3434  integer:: ipiv(lmax2)
3435  real(dp)::r(3,lmax2),rp(3,lmax2)
3436  real(dp)::rs2,rs3,rs6,rs12,umat(3,3)
3437  complex(dp)::crot(lmax2,lmax2),ctor(lmax2,lmax2),orb_lm(lmax2,-5:3,7)
3438  complex(dp):: ylmcp(lmax2)
3439  complex(dp):: ylmc_rr(lmax2,lmax2),ylmc_rr_save(lmax2,lmax2)
3440  complex(dp):: ylmc_rrinv(lmax2,lmax2),ylmc_rp(lmax2,lmax2)
3441  complex(dp),parameter :: c0=(0._dp,0._dp),c1=(1._dp,0._dp),ci=(0._dp,1._dp)
3442  character(len=500) :: message                   ! to be uncommented, if needed
3443 
3444 ! *************************************************************************
3445 
3446 
3447 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3448 !DEBUG
3449 !write(std_out,*)'lmax ',lmax,'lmax2 ',lmax2
3450 !write(std_out,*)'mband ',mband,'nwan ',nwan
3451 !
3452 !do iwan=1,nwan
3453 !write(std_out,*)'iwan,proj_l, proj_m',proj_l(iwan),proj_m(iwan)
3454 !write(std_out,*)'iwan,proj_x, proj_z',iwan,proj_x(:,iwan),proj_z(:,iwan)
3455 !end do
3456 !!END DEBUG
3457 
3458 !constants for linear combinations of ylm's
3459  rs2=1._dp/sqrt(2._dp)
3460  rs3=1._dp/sqrt(3._dp)
3461  rs6=1._dp/sqrt(6._dp)
3462  rs12=1._dp/sqrt(12._dp)
3463 
3464 !complex lm coefficients for real spherical harmonics in conventional order
3465 !s, py,pz,px, dxy,dyz,dz2,dxz,dx2-y2, fy(3x2-y2),fxyz,fyz2,fz3,fxz2,
3466 !fz(x2-y2),fx(x2-3y2)
3467  ctor(:,:)=c0
3468  do ll=0,lmax
3469    mm=0
3470    lm= ll**2+ll+mm+1
3471    ctor(lm,lm)=c1
3472    if(ll>0) then
3473      onem=one
3474      do mm=1,ll
3475        onem=-onem !(-1^mm)
3476        lm= ll**2+ll+mm+1
3477        lmc=ll**2+ll-mm+1
3478        ctor(lm ,lm )=rs2*c1
3479        ctor(lmc,lm )=onem*rs2*c1
3480        ctor(lm ,lmc)=rs2*ci
3481        ctor(lmc,lmc)=-onem*rs2*ci
3482      end do
3483    end if
3484  end do
3485 
3486  lm=0
3487  do ll=0,lmax
3488    do mm=-ll,ll
3489      lm=lm+1
3490      ctor(:,lm)=ctor(:,lm)*conjg(ci)**ll
3491    end do !mm
3492  end do !ll
3493 
3494 
3495 !coefficients for basic wannier orbitals in Table 3.1 order
3496  orb_lm(:,:,:)=c0
3497  ii=0
3498  do ll=0,lmax
3499    do mr=1,2*ll+1
3500      ii=ii+1
3501      orb_lm(:,ll,mr)=ctor(:,orb_idx(ii))
3502    end do
3503  end do
3504 
3505 
3506 
3507 !coefficients for linear combinations in table 3.2 order
3508  if(lmax>=1) then
3509 !  s            px
3510    orb_lm(:,-1,1)=rs2*ctor(:,1)+rs2*ctor(:,4)
3511    orb_lm(:,-1,2)=rs2*ctor(:,1)-rs2*ctor(:,4)
3512 !  s            px            py
3513    orb_lm(:,-2,1)=rs3*ctor(:,1)-rs6*ctor(:,4)+rs2*ctor(:,2)
3514    orb_lm(:,-2,2)=rs3*ctor(:,1)-rs6*ctor(:,4)-rs2*ctor(:,2)
3515    orb_lm(:,-2,3)=rs3*ctor(:,1)+2._dp*rs6*ctor(:,4)
3516 !  s        px        py        pz
3517    orb_lm(:,-3,1)=half*(ctor(:,1)+ctor(:,4)+ctor(:,2)+ctor(:,3))
3518    orb_lm(:,-3,2)=half*(ctor(:,1)+ctor(:,4)-ctor(:,2)-ctor(:,3))
3519    orb_lm(:,-3,3)=half*(ctor(:,1)-ctor(:,4)+ctor(:,2)-ctor(:,3))
3520    orb_lm(:,-3,4)=half*(ctor(:,1)-ctor(:,4)-ctor(:,2)+ctor(:,3))
3521  end if
3522  if(lmax>=2) then
3523 !  s            px            py
3524    orb_lm(:,-4,1)=rs3*ctor(:,1)-rs6*ctor(:,4)+rs2*ctor(:,2)
3525    orb_lm(:,-4,2)=rs3*ctor(:,1)-rs6*ctor(:,4)-rs2*ctor(:,2)
3526    orb_lm(:,-4,3)=rs3*ctor(:,1)+2._dp*rs6*ctor(:,4)
3527 !  pz           dz2
3528    orb_lm(:,-4,4)= rs2*ctor(:,3)+rs2*ctor(:,7)
3529    orb_lm(:,-4,5)=-rs2*ctor(:,3)+rs2*ctor(:,7)
3530 !  s            px            dz2         dx2-y2
3531    orb_lm(:,-5,1)=rs6*ctor(:,1)-rs2*ctor(:,4)-rs12*ctor(:,7)+half*ctor(:,9)
3532    orb_lm(:,-5,2)=rs6*ctor(:,1)+rs2*ctor(:,4)-rs12*ctor(:,7)+half*ctor(:,9)
3533 !  s            py            dz2         dx2-y2
3534    orb_lm(:,-5,3)=rs6*ctor(:,1)-rs2*ctor(:,2)-rs12*ctor(:,7)-half*ctor(:,9)
3535    orb_lm(:,-5,4)=rs6*ctor(:,1)+rs2*ctor(:,2)-rs12*ctor(:,7)-half*ctor(:,9)
3536 !  s            pz           dz2
3537    orb_lm(:,-5,5)=rs6*ctor(:,1)-rs2*ctor(:,3)+rs3*ctor(:,7)
3538    orb_lm(:,-5,6)=rs6*ctor(:,1)+rs2*ctor(:,3)+rs3*ctor(:,7)
3539  end if
3540 
3541 !stuff complex wannier orbital coefficient array
3542  do iwan=1,nwan
3543    ylmc_fac(:,iwan)=orb_lm(:,proj_l(iwan),proj_m(iwan))
3544  end do
3545 
3546 
3547 !setup to rotate ylmc_fac to new axes if called for
3548 !skip if only s projectors are used
3549  if ( lmax>0 ) then
3550 !  generate a set of nr=lmax2 random vectors
3551 !  idum=123456
3552    do ir=1,lmax2
3553      do ii=1,3
3554        r(ii,ir) = uniformrandom(idum)-0.5d0
3555      end do !ii
3556      call ylm_cmplx(lmax,ylmcp,r(1,ir),r(2,ir),r(3,ir))
3557      ylmc_rr(ir,:)=conjg(ylmcp(:))
3558      ylmc_rr_save(ir,:)=conjg(ylmcp(:))
3559    end do !ir
3560 
3561    ylmc_rrinv(:,:)=c0
3562    do ii=1,lmax2
3563      ylmc_rrinv(ii,ii)=c1
3564    end do !ii
3565 !  calculate inverse of ylmc(ir,lm) matrix
3566    call ZGESV(lmax2,lmax2,ylmc_rr,lmax2,ipiv,ylmc_rrinv,lmax2,info)
3567 
3568 !  check that r points are independent (ie., that matrix inversion wasn't
3569 !  too close to singular)
3570    ylmc_rr=matmul(ylmc_rrinv,ylmc_rr_save)
3571    test=zero
3572    do ii=1,lmax2
3573      ylmc_rr(ii,ii)=ylmc_rr(ii,ii)-c1
3574      do jj=1,lmax2
3575        test=max(abs(ylmc_rr(ii,jj)),test)
3576      end do !ii
3577    end do !jj
3578    if(test>tol8) then
3579      write(message, '(5a)' )&
3580 &     '  matrix inversion error for wannier rotations',ch10,&
3581 &     '  random vectors r(j,1:nr) are not all independent !! ',ch10,&
3582 &     '  Action : re-seed uniformrandom or maybe just try again'
3583      MSG_ERROR(message)
3584    end if !test>tol8
3585 
3586 !  end of the preliminaries, now to the rotations of the wannier orbitals
3587    do iwan=1,nwan
3588 !    don't bother for s orbitals
3589      if(proj_l(iwan)==0) cycle
3590 !    check for default axes and cycle if found
3591      if(proj_z(1,iwan)==zero .and. proj_z(2,iwan)==zero .and.&
3592 &     proj_z(3,iwan)== one .and. proj_x(1,iwan)==one .and.&
3593 &     proj_x(2,iwan)==zero .and. proj_x(3,iwan)==zero) cycle
3594 
3595 !    get the u matrix that rotates the reference frame
3596      call rotmat(proj_x(:,iwan),proj_z(:,iwan),inversion_flag,umat)
3597 
3598 !    find rotated r-vectors. Optional inversion
3599 !    operation is an extension of the wannier90 axis-setting options
3600 !    which only allow for proper axis rotations
3601      if(inversion_flag==1) then
3602        rp(:,:)= -matmul ( umat(:,:),  r(:,:) )
3603      else
3604        rp(:,:) = matmul ( umat(:,:) , r(:,:) )
3605      end if !inversion_flag
3606 
3607      do ir=1,lmax2
3608 !      get the ylm representation of the rotated vectors
3609        call ylm_cmplx(lmax,ylmcp,rp(1,ir),rp(2,ir),rp(3,ir))
3610        ylmc_rp(ir,:)=conjg(ylmcp(:))
3611      end do !ir
3612 !    the matrix product sum(ir) ylmc_rrinv(lm,ir)*ylmc_rp(ir,lm') gives the
3613 !    the complex lmXlm matrix representation of the coordinate rotation
3614      crot(:,:)=matmul(ylmc_rrinv(:,:),ylmc_rp(:,:))
3615 
3616 !    now rotate the current wannier orbital
3617      ylmcp(:)=matmul(crot(:,:),ylmc_fac(:,iwan))
3618      ylmc_fac(:,iwan)=ylmcp(:)
3619 
3620 !    write(std_out,*)'ylmc_fac',ylmc_fac(:,iwan)
3621    end do !iwan
3622  end if !lmax>0
3623 
3624 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

PARENTS

      mlwfovlp_projpaw

CHILDREN

      initylmr,matrginv,rotmat

SOURCE

3673 subroutine mlwfovlp_ylmfar(ylmr_fac,lmax,lmax2,mband,nwan,proj_l,proj_m,proj_x,proj_z)
3674 
3675 
3676 !This section has been created automatically by the script Abilint (TD).
3677 !Do not modify the following lines by hand.
3678 #undef ABI_FUNC
3679 #define ABI_FUNC 'mlwfovlp_ylmfar'
3680 !End of the abilint section
3681 
3682  implicit none
3683 
3684 !Arguments ------------------------------------
3685  integer, intent(in):: lmax,lmax2,nwan,mband
3686 ! arrays
3687  integer,intent(in) :: proj_l(mband),proj_m(mband)
3688  real(dp),intent(in) :: proj_x(3,mband),proj_z(3,mband)
3689  real(dp),intent(out)::ylmr_fac(lmax2,nwan)
3690 !
3691 !Local variables-------------------------------
3692 !
3693  integer :: idum,ii,inversion_flag
3694  integer :: ir,iwan,jj,ll,lm,mm,mr
3695  real(dp):: onem,test
3696 ! arrays
3697  real(dp),allocatable::dummy(:,:),nrm(:)
3698  real(dp)::r(3,lmax2),rp(3,lmax2)
3699  real(dp)::rs2,rs3,rs6,rs12,umat(3,3)
3700  real(dp)::rot(lmax2,lmax2),tor(lmax2,lmax2),orb_lm(lmax2,-5:3,7)
3701  real(dp):: ylmrp(lmax2)
3702  real(dp):: ylmr_rr(lmax2,lmax2),ylmr_rr_save(lmax2,lmax2)
3703  real(dp):: ylmr_rrinv(lmax2,lmax2),ylmr_rp(lmax2,lmax2)
3704  character(len=500) :: message                   ! to be uncommented, if needed
3705 !no_abirules
3706 !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
3707 
3708 ! *************************************************************************
3709 
3710 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3711 !DEBUG
3712 !write(std_out,*)'lmax ',lmax,'lmax2 ',lmax2
3713 !write(std_out,*)'mband ',mband,'nwan ',nwan
3714 !
3715 !do iwan=1,nwan
3716 !write(std_out,*)'iwan,proj_l, proj_m',proj_l(iwan),proj_m(iwan)
3717 !write(std_out,*)'iwan,proj_x, proj_z',iwan,proj_x(:,iwan),proj_z(:,iwan)
3718 !end do
3719 !!END DEBUG
3720 
3721 !constants for linear combinations of ylm's
3722  rs2=1._dp/sqrt(2._dp)
3723  rs3=1._dp/sqrt(3._dp)
3724  rs6=1._dp/sqrt(6._dp)
3725  rs12=1._dp/sqrt(12._dp)
3726 
3727 !
3728 !mapping lm coefficients for real spherical harmonics
3729 !table 3.1 of Wannier90 user guide with real spherical harmonics in routine initylmr
3730 !s, py,pz,px, dxy,dyz,dz2,dxz,dx2-y2, fy(3x2-y2),fxyz,fyz2,fz3,fxz2,
3731 !fz(x2-y2),fx(x2-3y2)
3732 !note: check ordering of f orbitals, it might be wrong
3733 
3734  tor(:,:)=0.d0
3735  lm=0
3736  do ll=0,lmax
3737    do mm=-ll,ll
3738      onem=(-1.d0)**mm
3739      lm=lm+1
3740      if(ll == 0) then
3741        tor(lm,lm)=1.d0
3742      else
3743        tor(lm,lm)=onem*1.d0
3744      end if
3745    end do !mm
3746  end do !ll
3747 !do lm=1,16
3748 !write(std_out,*)'tor lm=',lm,tor(:,lm)
3749 !end do
3750 
3751 !coefficients for basic wannier orbitals in Table 3.1 order
3752  orb_lm(:,:,:)=0.d0
3753  ii=0
3754  do ll=0,lmax
3755    do mr=1,2*ll+1
3756      ii=ii+1
3757      orb_lm(:,ll,mr)= tor(:,ii)
3758 !    write(std_out,*)'ii',ii,'orb_lm',orb_lm(:,ll,mr)
3759    end do
3760  end do
3761 
3762 
3763 
3764 !coefficients for linear combinations in table 3.2 order
3765  if(lmax>=1) then
3766 !  s            px
3767    orb_lm(:,-1,1)=rs2*tor(:,1)+rs2*tor(:,4)
3768    orb_lm(:,-1,2)=rs2*tor(:,1)-rs2*tor(:,4)
3769 !  s            px            py
3770    orb_lm(:,-2,1)=rs3*tor(:,1)-rs6*tor(:,4)+rs2*tor(:,2)
3771    orb_lm(:,-2,2)=rs3*tor(:,1)-rs6*tor(:,4)-rs2*tor(:,2)
3772    orb_lm(:,-2,3)=rs3*tor(:,1)+2._dp*rs6*tor(:,4)
3773 !  s        px        py        pz
3774    orb_lm(:,-3,1)=half*(tor(:,1)+tor(:,4)+tor(:,2)+tor(:,3))
3775    orb_lm(:,-3,2)=half*(tor(:,1)+tor(:,4)-tor(:,2)-tor(:,3))
3776    orb_lm(:,-3,3)=half*(tor(:,1)-tor(:,4)+tor(:,2)-tor(:,3))
3777    orb_lm(:,-3,4)=half*(tor(:,1)-tor(:,4)-tor(:,2)+tor(:,3))
3778  end if
3779  if(lmax>=2) then
3780 !  s            px            py
3781    orb_lm(:,-4,1)=rs3*tor(:,1)-rs6*tor(:,4)+rs2*tor(:,2)
3782    orb_lm(:,-4,2)=rs3*tor(:,1)-rs6*tor(:,4)-rs2*tor(:,2)
3783    orb_lm(:,-4,3)=rs3*tor(:,1)+2._dp*rs6*tor(:,4)
3784 !  pz           dz2
3785    orb_lm(:,-4,4)= rs2*tor(:,3)+rs2*tor(:,7)
3786    orb_lm(:,-4,5)=-rs2*tor(:,3)+rs2*tor(:,7)
3787 !  s            px            dz2         dx2-y2
3788    orb_lm(:,-5,1)=rs6*tor(:,1)-rs2*tor(:,4)-rs12*tor(:,7)+half*tor(:,9)
3789    orb_lm(:,-5,2)=rs6*tor(:,1)+rs2*tor(:,4)-rs12*tor(:,7)+half*tor(:,9)
3790 !  s            py            dz2         dx2-y2
3791    orb_lm(:,-5,3)=rs6*tor(:,1)-rs2*tor(:,2)-rs12*tor(:,7)-half*tor(:,9)
3792    orb_lm(:,-5,4)=rs6*tor(:,1)+rs2*tor(:,2)-rs12*tor(:,7)-half*tor(:,9)
3793 !  s            pz           dz2
3794    orb_lm(:,-5,5)=rs6*tor(:,1)-rs2*tor(:,3)+rs3*tor(:,7)
3795    orb_lm(:,-5,6)=rs6*tor(:,1)+rs2*tor(:,3)+rs3*tor(:,7)
3796  end if
3797 
3798 !real wannier orbital coefficient array
3799  do iwan=1,nwan
3800    ylmr_fac(:,iwan)=orb_lm(:,proj_l(iwan),proj_m(iwan))
3801  end do
3802 
3803 
3804 !setup to rotate ylmr_fac to new axes if called for
3805 !skip if only s projetors are used
3806  if ( lmax>0 ) then
3807 !  generate a set of nr=lmax2 random vetors
3808    idum=123456
3809    do ir=1,lmax2
3810      do ii=1,3
3811        r(ii,ir) = uniformrandom(idum)-0.5d0
3812      end do !ii
3813    end do !ir
3814    ABI_ALLOCATE(nrm,(lmax2))
3815    nrm(:)=sqrt(r(1,:)**2+r(2,:)**2+r(3,:)**2)**0.5
3816    call initylmr(lmax+1,1,lmax2,nrm,1,r(:,:),ylmr_rr_save(:,:),dummy)
3817    ylmr_rr(:,:)=ylmr_rr_save(:,:)
3818    do ir=1,lmax2
3819      ylmr_rr_save(ir,:)=ylmr_rr(:,ir)
3820    end do
3821    ABI_DEALLOCATE(nrm)
3822 
3823    ylmr_rrinv(:,:)=0.d0
3824    do ii=1,lmax2
3825      ylmr_rrinv(ii,ii)=1.d0
3826    end do !ii
3827 !  calculate inverse of ylmr(ir,lm) matrix
3828    ylmr_rrinv(:,:)=ylmr_rr_save(:,:)
3829    call matrginv(ylmr_rrinv,lmax2,lmax2)
3830 
3831 !  check that r points are independent (ie., that matrix inversion wasn't
3832 !  too close to singular)
3833    ylmr_rr=matmul(ylmr_rrinv,ylmr_rr_save)
3834    test=0.d0
3835    do ii=1,lmax2
3836      ylmr_rr(ii,ii)=ylmr_rr(ii,ii)-1.d0
3837      do jj=1,lmax2
3838        test=max(abs(ylmr_rr(ii,jj)),test)
3839      end do !ii
3840    end do !jj
3841    if(test>tol8) then
3842      write(message, '(5a)' )&
3843 &     '  matrix inversion error for wannier rotations',ch10,&
3844 &     '  random vetors r(j,1:nr) are not all independent !! ',ch10,&
3845 &     '  Action : re-seed uniformrandom or maybe just try again'
3846      MSG_ERROR(message)
3847    end if !test>tol8
3848 
3849 !  end of the preliminaries, now to the rotations of the wannier orbitals
3850    do iwan=1,nwan
3851 !    don't bother for s orbitals
3852      if(proj_l(iwan)==0) cycle
3853 !    check for default axes and cycle if found
3854      if(proj_z(1,iwan)==0.d0 .and. proj_z(2,iwan)==0.d0 .and.&
3855 &     proj_z(3,iwan)== 1.d0 .and. proj_x(1,iwan)==1.d0 .and.&
3856 &     proj_x(2,iwan)==0.d0 .and. proj_x(3,iwan)==0.d0) cycle
3857 
3858 !    get the u matrix that rotates the reference frame
3859      call rotmat(proj_x(:,iwan),proj_z(:,iwan),inversion_flag,umat)
3860 !
3861 !    find rotated r-vetors. Optional inversion
3862 !    operation is an extension of the wannier90 axis-setting options
3863 !    which only allow for proper axis rotations
3864      if(inversion_flag==1) then
3865        rp(:,:)= -matmul ( umat(:,:),  r(:,:) )
3866      else
3867        rp(:,:) = matmul ( umat(:,:) , r(:,:) )
3868      end if !inversion_flag
3869 
3870 !    get the ylm representation of the rotated vetors
3871      ABI_ALLOCATE(nrm,(lmax2))
3872      nrm(:)=sqrt(rp(1,:)**2+rp(2,:)**2+rp(3,:)**2)**0.5
3873      call initylmr(lmax+1,1,lmax2,nrm,1,rp(:,:),ylmr_rp(:,:),dummy)
3874      ylmr_rr(:,:)=ylmr_rp(:,:)
3875      do ir=1,lmax2
3876        ylmr_rp(ir,:)=ylmr_rr(:,ir)
3877      end do
3878      ABI_DEALLOCATE(nrm)
3879 !    the matrix product sum(ir) ylmr_rrinv(lm,ir)*ylmr_rp(ir,lm') gives the
3880 !    the  lmXlm matrix representation of the coordinate rotation
3881 
3882      rot(:,:)=matmul(ylmr_rrinv(:,:),ylmr_rp(:,:))
3883 !
3884 !    now rotate the current wannier orbital
3885      ylmrp(:)=matmul(rot(:,:),ylmr_fac(:,iwan))
3886      ylmr_fac(:,iwan)=ylmrp(:)
3887    end do !iwan
3888  end if !lmax>0
3889 
3890 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

PARENTS

CHILDREN

SOURCE

1391  subroutine read_chkunit(seed_name,nkpt,ndimwin,ierr)
1392 
1393 
1394 !This section has been created automatically by the script Abilint (TD).
1395 !Do not modify the following lines by hand.
1396 #undef ABI_FUNC
1397 #define ABI_FUNC 'read_chkunit'
1398 !End of the abilint section
1399 
1400  implicit none
1401 
1402 !Arguments ------------------------------------
1403 !scalars
1404  integer,intent(in) :: nkpt
1405  character(len=*),intent(in) :: seed_name
1406  integer,intent(out) :: ierr
1407 !arrays
1408  integer,intent(out) :: ndimwin(nkpt)
1409 
1410 !Local variables-------------------------------
1411 !scalars
1412  integer :: chk_unit,ios,ikpt
1413  logical :: have_disentangled
1414 
1415 !************************************************************************
1416 
1417    chk_unit=get_unit()
1418    fname=TRIM(seed_name)//'.chk'
1419    open(unit=chk_unit,file=fname,form='unformatted',status='old',iostat=ios)
1420 
1421    ierr=0
1422    read(chk_unit) ! header                                   ! Date and time
1423    read(chk_unit) ! ((real_lattice(i,j),i=1,3),j=1,3)        ! Real lattice
1424    read(chk_unit) ! ((recip_lattice(i,j),i=1,3),j=1,3)       ! Reciprocal lattice
1425    read(chk_unit) ! num_kpts
1426    read(chk_unit) ! ((kpt_latt(i,nkp),i=1,3),nkp=1,num_kpts) ! K-points
1427    read(chk_unit) ! nntot                  ! Number of nearest k-point neighbours
1428    read(chk_unit) ! num_wann               ! Number of wannier functions
1429    read(chk_unit) ! chkpt1                 ! Position of checkpoint
1430    read(chk_unit) have_disentangled        ! Whether a disentanglement has been performed
1431    if (have_disentangled) then
1432 !    read(chk_unit) ! omega_invariant     ! Omega invariant
1433 !    read(chk_unit) ((lwindow(i,nkp),i=1,num_bands),nkp=1,num_kpts)
1434      read(chk_unit) (ndimwin(ikpt),ikpt=1,nkpt)
1435 !    read(chk_unit) (((u_matrix_opt(i,j,nkp),i=1,num_bands),j=1,num_wann),nkp=1,num_kpts)
1436    else
1437 !    this is not expected. we should have disentanglement. Report the error.
1438      ierr=-1
1439    end if
1440 !  read(chk_unit)  (((u_matrix(i,j,k),i=1,num_wann),j=1,num_wann),k=1,num_kpts)               ! U_matrix
1441 !  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
1442 !  read(chk_unit)  ((wannier_centres(i,j),i=1,3),j=1,num_wann)
1443    close(chk_unit)
1444 
1445 end subroutine read_chkunit