Logo Search packages:      
Sourcecode: abinit version File versions  Download package

localorb_S.F90

!{\src2tex{textfont=tt}}
!!****f* ABINIT/localorb_S
!! NAME
!! localorb_S
!!
!! FUNCTION
!! Construction of Wannier type localized orbitals (WLO)
!! Ref:cond-mat/0509571
!!
!! COPYRIGHT
!! Copyright (C) 2005-2007 ABINIT group (JB)
!! this file is distributed under the terms of the
!! gnu general public license, see ~abinit/infos/copyright
!! or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!! Needs an unformatted wave function from abinit.
!! chr_inputfname = name of the input file for WLO construction
!! exchn2n3=if 1, n2 and n3 are exchanged
!! headform=format of the wf file
!!
!! natom = number of atoms in the unit cell
!! nr1,nr2,nr3 = grid size (nr1 x nr2 x nr3 = filrho dimension)
!! ntypat = number of atom type
!! ucvol = unit cell volume (> 0)
!! densfileformat = flag for the format of the density file:
!!       0 = ASCII
!!       1 = binary
!! denval = density value exported by interpol, to be wrote in the output file
!! filrho = name of the density file (ASCII or binary)
!! filtau = name of the atomic position file (Xmol format)
!! tau = cartesian coordinates
!! acell = unit cell parameters
!! rprim = orientation of the unit cell axes
!!
!! OUTPUT
!! xcrysden plottable WLO file,
!! corresponding complex data file
!! other analysis and runtime check files

!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!      cut3d
!!
!! CHILDREN
!!      instrng,intagm,inupper,kptindex,metric,overlap_wf,wfread,xredxcart
!!      zgesvd
!!
!! SOURCE

#if defined HAVE_CONFIG_H
#include "config.h"
#endif

subroutine localorb_s(chr_inputfname,ecut,exchn2n3,headform,istwfk,kpt,natom,nband,nkpt,&
                       npwarr,nr1,nr2,nr3,nspinor,nsppol,ntypat,rprimd,tau,typat,znucl)

 use defs_basis
 use defs_datatypes

!This section has been created automatically by the script Abilint (TD). Do not modify these by hand.
#ifdef HAVE_FORTRAN_INTERFACES
 use interfaces_11util
 use interfaces_12geometry
 use interfaces_12parser
 use interfaces_19cut3d, except_this_one => localorb_s
#endif
!End of the abilint section

 implicit none

!Arguments -----------------------------------
!scalars
 integer,intent(in) :: exchn2n3,headform,natom,nkpt,nr1,nr2,nr3,nspinor,nsppol
 integer,intent(in) :: ntypat
 real(dp),intent(in) :: ecut
 character(len=50),intent(in) :: chr_inputfname
!arrays
 integer,intent(in) :: istwfk(nkpt),nband(nkpt),npwarr(nkpt),typat(natom)
 real(dp),intent(in) :: kpt(3,nkpt),rprimd(3,3),znucl(ntypat)
 real(dp),intent(inout) :: tau(3,natom)

!Local variables-------------------------------
!,unout=12
       character(*), parameter :: inputfile='cut.in'
       character*(fnlen) :: ylmnam=" "
! for input file reading :
! for wf :
!scalars
 integer,save :: tim_fourwf=0,tim_rwwf=0
 integer :: cband,cgshift,ckpt,ckpt1,ckpt2,cplex,cspinor,csppol,formeig
 integer :: gridshift1,gridshift2,gridshift3,i,i1,i10,i11,i12,i2,i3,i4,i5,i6,i7
 integer :: i8,i9,ia,iatom,iband,iband1,iband2,icellplot,icellplotx,icellploty
 integer :: icellplotz,icellx,icelly,icellz,icg,ichoice,icountorb,ifile,igeom
 integer :: igx,igx1,igy,igy1,igz,igz1,ii,ii1,ii2,ii3,ikpt,ikpt_x,ikpt_x1
 integer :: ikpt_x2,ikpt_y,ikpt_y1,ikpt_y2,ikpt_z,ikpt_z1,ikpt_z2,ilang,info
 integer :: init_prefact,insmet,ioffkg,iorb,ios,iout,iprompt,ipw,ir1,ir2,ir3
 integer :: ispden,isppol,istop1,ivar,ivect,iwan,iwrite,ix,ixfh,ixint,iy,iz,j
 integer :: j1,j2,k,k1,k2,l,lenline,lenstr,m,m1,marr,matom,mbess,mcg,mgfft
 integer :: mkmem,mlang,mpw,munit,munitx,munity,munitz,n,n4,n5,n6,nband_disk
 integer :: nbands,ncellx,ncellx1,ncelly,ncelly1,ncellz,ncellz1,ndimcsvd,ndisk1
 integer :: ndisk2,nfit,ngx,ngy,ngz,njob,nkpts,nkx,nky,nkz,nline,nlinemax,norb
 integer :: nplwv,npw_k,npwout,nradint,nspden,nstartx,nstartx1,nstarty,nstarty1
 integer :: nstartz,nstartz1,nstopx,nstopx1,nstopy,nstopy1,nstopz,nstopz1,nvar
 integer :: nwancons,nwanmax1,nwanmax2,nwanmax3,nwantot,nxfh,nxhw,nxmax,nxmin
 integer :: nxtot,nxw,nyhw,nymax,nymin,nytot,nyw,nzhw,nzmax,nzmin,nztot,nzw
 integer :: oldcband,oldckpt,oldcspinor,oldcsppol,option,prtsphere,select_exit
 integer :: tread,unkg,unylm=0
 real(dp) :: alpha,arg,bessargmax,bessint_delta,du,efermi,energy,envelop,hx,hy
 real(dp) :: hz,kpgmax,normtot,phy,power,r1,r2,r3,ratsph,rcellx,rcelly,rcellz
 real(dp) :: rcoordx,rcoordy,rcoordz,rcos1,rcos2,rcos3,rcosphy,rcostheta,re1
 real(dp) :: re2,re3,re4,re5,re6,re7,re8,rhx,rhy,rhz,rk,rkpt,rkpt1,rkpt2,rkpt3
 real(dp) :: rkptx,rkpty,rkptz,rkx,rkx1,rkx2,rkx3,rky,rky1,rky2,rky3,rkz,rkz1
 real(dp) :: rkz2,rkz3,rlx,rlx1,rlx1_large,rlx2,rlx2_large,rlx3,rlx3_large
 real(dp) :: rlx_large,rly,rly1,rly1_large,rly2,rly2_large,rly3,rly3_large
 real(dp) :: rly_large,rlz,rlz1,rlz1_large,rlz2,rlz2_large,rlz3,rlz3_large
 real(dp) :: rlz_large,rn1,rn2,rn3,rnorm,rnorm2,roccup,rpi,rrcosphy2,rrsinphy2
 real(dp) :: rsin1,rsin2,rsin3,rsinphy,rsintheta,rwanmax,rx,rxcm,ry,rycm,rz
 real(dp) :: rzcm,shiftx,shifty,shiftz,thetax,thetay,thetaz,tmpi,tmpr,tpi
 real(dp) :: tsmear,ucvol,weight,x,x1,xnow,y,y1,ynow,z,z1,znow
 character(len=1) :: chr2,crhsvd
 character(len=15) :: chr1
 character(len=150) :: chr_string,chr_string1
 character(len=2) :: chr
 character(len=30) :: token
 character(len=4) :: mode_paral
 character(len=50) :: chr_input,chr_wanindex
 character(len=500) :: message
 character(len=fnlen) :: kgnam,output,output1
 character(len=strlen) :: string
 character :: outputchar
 type(hdr_type) :: hdr
 type(mpi_type) :: mpi_enreg
!arrays
 integer :: ngfft(8)
 integer,allocatable :: intarr(:),ltypeorb(:,:),numorb(:,:),nwanindex(:)
 real(dp) :: acell(3),cmax(natom),gmet(3,3),gprimd(3,3),oldkpt(3)
 real(dp) :: phkxred(2,natom),rmet(3,3),shift_tau(3),tau2(3,natom),vec(3)
 real(dp) :: vec1(3),xcart(3,natom),xmat(3,3),xred(3,natom),ymat(3,3),zmat(3,3)
 real(dp),allocatable :: ai(:,:),ar(:,:),dprarr(:),e_kpt(:),fv1(:),fv2(:)
 real(dp),allocatable :: fv3(:),rcoord(:,:),rcoordall(:,:,:),rprojx(:)
 real(dp),allocatable :: rprojy(:),rprojz(:),rtheta(:,:,:),rwork(:),s(:),wi(:)
 real(dp),allocatable :: wr(:),zi(:,:),zr(:,:)
 character(len=15),allocatable :: chrhead(:)
 character(len=fnlen),allocatable :: filename(:)
!no_abirules
 complex(dp) :: cmp1,cmp2,cmp3,cmp4,csum
 complex(dp) :: cf1,cf2,cf3,cf4,cf5,cf6,cf7,cf8,cf9,cf10,cf11,cf12,cf13,cf14,cf15
 complex(dp) :: cwfs1,cwfs2,cwfs3,cwfs4,cwfs5,cwfs6
 complex(kind =8)::orbital,corbital
 complex(dp),allocatable,dimension(:,:)::      cm1,cm2,cm4,cm5,coverlap
 complex(dp),allocatable,dimension(:,:)::      cmat1,cmat2
 complex(dp),allocatable,dimension(:,:) ::     cu,cvt
 complex(dp),allocatable:: csvd(:)
 complex(dp),allocatable ::cwannier(:,:,:) 
 complex(dp),allocatable ::crotl(:,:,:,:,:)
 complex(dp),allocatable :: cbfns(:,:,:),cbfns1(:,:,:),corb(:,:,:)

! *************************************************************************

#ifdef VMS
!DEC$ ATTRIBUTES ALIAS:'ZGESVD' :: zgesvd
#endif

 nbands = maxval(nband)
 ngz = nr3
 ngy = nr2
 ngx = nr1
 tpi=two_pi
 ndimcsvd = 3*nbands + 10
 allocate(csvd(ndimcsvd))

!Finding out nkx,nky,nkz from the kpt array:

 nkx = 1
 nky = 1
 nkz = 1

 do i = 1,nkpt

  re1 = kpt(1,i)
  re2 = kpt(2,i)
  re3 = kpt(3,i)

  if(i.gt.1)then
   do j=1,i-1
    if(re1.eq.kpt(1,j))go to 100
   end do
   nkx = nkx + 1
100 continue
  end if

  if(i.gt.1)then
   do j =1,i-1
    if(re2.eq.kpt(2,j))go to 200
   end do
   nky = nky + 1
200 continue
  end if

  if(i.gt.1)then
   do j =1,i-1
    if(re3.eq.kpt(3,j))go to 300
   end do
   nkz = nkz + 1
300 continue
  end if
 end do ! i=1,nkpt

!==================================================================
!reading the input file for wlo :
!default:
 matom=0
 insmet = 1
 iwrite = 0
 alpha = 1.0_dp
 ncellx1 = 1
 ncelly1 = 1
 ncellz1 = 1
 norb = nbands
 rcellx = 0.0_dp
 rcelly = 0.0_dp
 rcellz = 0.0_dp
 njob = 0
 tsmear=zero
 efermi=zero

 nvar = 13 ! maximum number of variales to be read from the input files
 allocate(chrhead(nvar))
 chrhead(1)='ncenter'
 chrhead(2)='norb'
 chrhead(3)='supercell'
 chrhead(4)='alpha'
 chrhead(5)='lofwrite'
 chrhead(6)='shiftk'
 chrhead(7)='insmet'
 chrhead(8)='efermi'
 chrhead(9)='tsmear'
 chrhead(10)='orbdetails'
 chrhead(11)='wfconstrctn'
 chrhead(12)='wfcell'
 chrhead(13)='jobtype'

 call instrng(chr_inputfname,lenstr,1,strlen,string)
!To make case-insensitive, map characters of string to upper case:
 call inupper(string(1:lenstr))

 marr=3
 allocate(intarr(3),dprarr(3))

 call intagm(dprarr,intarr,0,marr,1,string(1:lenstr),chrhead(1),tread,'INT')
 if(tread==1) matom=intarr(1)

 call intagm(dprarr,intarr,0,marr,1,string(1:lenstr),chrhead(2),tread,'INT')
 if(tread==1) norb=intarr(1)

 if(norb.gt.nbands)then
  print*,'chosen number of bands :',norb
  print*,'nbands in scf :',nbands
  print*,'the former cannot be greater than the latter.'
  stop
 end if
 nbands = norb

 call intagm(dprarr,intarr,0,marr,3,string(1:lenstr),chrhead(3),tread,'INT')
 if(tread==1)then
  ncellx1=intarr(1)
  ncelly1=intarr(2)
  ncellz1=intarr(3)
 end if

 call intagm(dprarr,intarr,0,marr,1,string(1:lenstr),chrhead(4),tread,'DPR')
 if(tread==1) alpha=dprarr(1)

 call intagm(dprarr,intarr,0,marr,1,string(1:lenstr),chrhead(5),tread,'INT')
 if(tread==1) iwrite=intarr(1)

 call intagm(dprarr,intarr,0,marr,3,string(1:lenstr),chrhead(6),tread,'DPR')
 if(tread==1)then
  shiftx=dprarr(1)
  shifty=dprarr(2)
  shiftz=dprarr(3)
 end if

 call intagm(dprarr,intarr,0,marr,1,string(1:lenstr),chrhead(7),tread,'INT')
 if(tread==1) insmet=intarr(1)

 call intagm(dprarr,intarr,0,marr,1,string(1:lenstr),chrhead(8),tread,'DPR')
 if(tread==1) efermi=dprarr(1)

 call intagm(dprarr,intarr,0,marr,1,string(1:lenstr),chrhead(9),tread,'DPR')
 if(tread==1) tsmear=dprarr(1)

 call intagm(dprarr,intarr,0,marr,3,string(1:lenstr),chrhead(12),tread,'DPR')
 if(tread==1)then
  rcellx=intarr(1)
  rcelly=intarr(2)
  rcellz=intarr(3)
 end if

 call intagm(dprarr,intarr,0,marr,1,string(1:lenstr),chrhead(13),tread,'INT')
 if(tread==1) njob=intarr(1)

 deallocate(intarr,dprarr)
 marr=3*nsppol*matom
 allocate(intarr(marr),dprarr(marr))

 allocate(numorb(nsppol,matom),rcoord(3,matom),ltypeorb(nbands,nsppol))
 allocate(e_kpt(nbands),rcoordall(3,nbands,nsppol))

 token = 'rcoord'
 call intagm(dprarr,intarr,0,marr,3*matom,string(1:lenstr),token,tread,'DPR')
 if(tread==1)rcoord(:,:)=reshape(dprarr(1:3*matom),(/3,matom/))

 token = 'numorb'
 call intagm(dprarr,intarr,0,marr,nsppol*matom,string(1:lenstr),token,tread,'INT')
 if(tread==1)numorb(:,:)=reshape(intarr(1:nsppol*matom),(/nsppol,matom/))

 allocate(rtheta(nbands,3,nsppol))
 ii1=sum(numorb(1:nsppol,1:matom))
 if(nbands/=ii1)then
  write(6,'(a,a,a,i4,a,i4,a,a)')' localorb_S : ERROR ',ch10,&
&  ' The sum of the number of orbitals ',ii1,', is not equal to the number of bands =',nbands,ch10,&
&  ' Stop '
  stop
 end if

 do csppol = 1,nsppol
  ii1=0 ; ii2=1
  do i = 1,matom
   j1=numorb(csppol,i) ; ii1 = ii1 + j1
   do k = ii2,ii1
    rcoordall(1:3,k,csppol)=rcoord(1:3,i)
   end do
   ii2 = ii2 +j1
  end do
 end do

 deallocate(intarr,dprarr)
 marr=nbands*3*nsppol
 allocate(intarr(marr),dprarr(marr))

 token = 'ltypeorb'
 call intagm(dprarr,intarr,0,marr,nbands*nsppol,string(1:lenstr),token,tread,'INT')
 if(tread==1)ltypeorb(:,:)=reshape(intarr(1:nbands*nsppol),(/nbands,nsppol/))

 token = 'rtheta'
 call intagm(dprarr,intarr,0,marr,nbands*3*nsppol,string(1:lenstr),token,tread,'DPR')
 if(tread==1)rtheta(:,:,:)=reshape(dprarr(1:nbands*3*nsppol),(/nbands,3,nsppol/))

 nwantot = nbands
 token = 'nwantot'
 call intagm(dprarr,intarr,0,marr,1,string(1:lenstr),token,tread,'INT')
 if(tread==1)nwantot=intarr(1)
 allocate(nwanindex(nwantot))

 if(nwantot == nbands)then
  do i1 = 1,nwantot
   nwanindex(i1) = i1
  end do
 else
  deallocate(intarr,dprarr)
  marr=nwantot
  allocate(intarr(marr),dprarr(marr))
  token = 'nwanindex'
  call intagm(dprarr,intarr,0,marr,nwantot,string(1:lenstr),token,tread,'INT')
  if(tread==1)nwanindex(1:nwantot)=intarr(1:nwantot)
 end if

 deallocate(intarr,dprarr)
!The reading of data is finished
!--------------------------------------------------------------------------------------

 if(rcellx+rcelly+rcellz.ne.0.0)then
  if( anint(rcellx).gt.ncellx1)rcellx=real(ncellx1)
  if( anint(rcelly).gt.ncelly1)rcelly=real(ncelly1)
  if( anint(rcellz).gt.ncellz1)rcellz=real(ncellz1)
 end if

 write(6,'(3a)')ch10,' Echoing out the read input variables for WLO calculation:',ch10

 write(6,'(a,a12,i4)')' ',chrhead(1),matom
 write(6,'(a,a12,i4)')' ',chrhead(2),nbands
 write(6,'(a,a12,3i4)')' ',chrhead(3),ncellx1,ncelly1,ncellz1
 if(rcellx+rcelly+rcellz.eq.0.0)then
  write(6,'(a,a12,3i4)')' ',chrhead(12),ncellx1,ncelly1,ncellz1
 else
  write(6,'(a,a12,3i4)')' ',chrhead(12),rcellx,rcelly,rcellz
 end if
 write(6,'(a,a12,es16.6)')' ',chrhead(4),alpha
 write(6,'(a,a12,i4)')' ',chrhead(5),iwrite
 write(6,'(a,a12,3es16.6)')' ',chrhead(6),shiftx,shifty,shiftz
 write(6,'(a,a12,i4)')' ',chrhead(7),insmet
 write(6,'(a,a12,es16.6)')' ',chrhead(8),efermi
 write(6,'(a,a12,es16.6)')' ',chrhead(9),tsmear
 write(6,'(a)')' For each atom :        rcoord(1:3)                                 numorb'
 do i=1,matom
  write(6,'(8x,i4,3es16.6,8x,2i3)')i,rcoord(:,i),numorb(1:nsppol,i)
 end do
 write(6,'(a)')' For each orbital : spin, ltypeorb     and      rtheta(1:3) '
 do csppol=1,nsppol
  do ii2=1,nbands
   write(6,'(8x,i4,2i10,4x,3es16.6)')ii2,csppol,ltypeorb(ii2,csppol),rtheta(ii2,:,csppol)
  end do
 end do
 write(6,'(1x,a12,i4)')chrhead(11),nwantot
 write(6,'(a)')' They are :'
 write(6,'(a,20i4)')' ',(nwanindex(i1),i1 =1,nwantot)

!-----------------------------------------------------------------------------
 allocate( corb(ngx,ngy,ngz),cbfns(ngx,ngy,ngz))
 allocate( cbfns1(ngx,ngy,ngz))

!svd :
 allocate(coverlap(nbands,nbands),cm1(nbands,nbands),cm2(nbands,nbands))
 allocate(cm4(nbands,nbands))
 allocate(cm5(nbands,nbands),cmat1(nbands,nbands),cmat2(nbands,nbands))
 allocate(s(nbands),cu(nbands,nbands),cvt(nbands,nbands),rwork(5*nbands))
 allocate(ar(nbands,nbands),ai(nbands,nbands),zr(nbands,nbands))
 allocate(zi(nbands,nbands),fv1(nbands))
 allocate(fv2(nbands),fv3(nbands),wr(nbands),wi(nbands))

 allocate( cwannier(ncellx1*ngx,ncelly1*ngy,ncellz1*ngz))
 allocate( crotl(nkx,nky,nkz,nbands,nbands))
!-------------------------------------------------------------------------------
!begin executable section
!ios = 0
 mpi_enreg%paralbd=0

 formeig=0
 oldckpt=0
 oldcband=0
 oldcsppol=0
 oldcspinor=0

 iout=-1
 call metric(gmet,gprimd,iout,rmet,rprimd,ucvol)

!get xred
 call xredxcart(natom,-1,rprimd,tau,xred)

 do i=1,natom
  xcart(:,i)=rprimd(:,1)*xred(1,i)+&
&            rprimd(:,2)*xred(2,i)+&
&            rprimd(:,3)*xred(3,i)
 end do

 tpi = 8.0*atan(1.0)
 rpi  = 4.0*atan(1.0)

 ngz = nr3
 ngy = nr2
 ngx = nr1

 if(ncellx1.ge.2)icellplotx = 2
 if(ncellx1.lt.2)icellplotx = 1

 if(ncelly1.ge.2)icellploty = 2
 if(ncelly1.lt.2)icellploty = 1

 if(ncellz1.ge.2)icellplotz = 2
 if(ncellz1.lt.2)icellplotz = 1

 rlx1_large=rprimd(1,1)*real(ncellx1)
 rly1_large=rprimd(2,1)*real(ncellx1)
 rlz1_large=rprimd(3,1)*real(ncellx1)

 rlx2_large=rprimd(1,2)*real(ncelly1)
 rly2_large=rprimd(2,2)*real(ncelly1)
 rlz2_large=rprimd(3,2)*real(ncelly1)

 rlx3_large=rprimd(1,3)*real(ncellz1)
 rly3_large=rprimd(2,3)*real(ncellz1)
 rlz3_large=rprimd(3,3)*real(ncellz1)

 rlx_large = rlx1_large + rlx2_large + rlx3_large
 rly_large = rly1_large + rly2_large + rly3_large
 rlz_large = rlz1_large + rlz2_large + rlz3_large

 rlx1=rprimd(1,1)
 rly1=rprimd(2,1)
 rlz1=rprimd(3,1)

 rlx2=rprimd(1,2)
 rly2=rprimd(2,2)
 rlz2=rprimd(3,2)

 rlx3=rprimd(1,3)
 rly3=rprimd(2,3)
 rlz3=rprimd(3,3)

 rlx = rlx1 + rlx2 + rlx3
 rly = rly1 + rly2 + rly3
 rlz = rlz1 + rlz2 + rlz3

 write(*, '(a,3i6)' )' real space grids:',ngx,ngy,ngz
 write(*, '(a,3i6)' )' k space grids:',nkx,nky,nkz
 write(*, '(a,3es16.6)' )' rlx rly rlz',rlx,rly,rlz

!=============================================================================
! reciprocal space vectors :
!----------------------------
        rkx1 = rly2*rlz3 - rly3*rlz2
        rky1 = rlz2*rlx3 - rlz3*rlx2
        rkz1 = rlx2*rly3 - rlx3*rly2
        rk   = rkx1*rlx1 + rky1*rly1 + rkz1*rlz1
        rkx1 = tpi*rkx1/rk
        rky1 = tpi*rky1/rk
        rkz1 = tpi*rkz1/rk

        rkx2 = rly3*rlz1 - rly1*rlz3
        rky2 = rlz3*rlx1 - rlz1*rlx3
        rkz2 = rlx3*rly1 - rlx1*rly3
        rk   = rkx2*rlx2 + rky2*rly2 + rkz2*rlz2
        rkx2 = tpi*rkx2/rk
        rky2 = tpi*rky2/rk
        rkz2 = tpi*rkz2/rk

        rkx3 = rly1*rlz2 - rly2*rlz1
        rky3 = rlz1*rlx2 - rlz2*rlx1
        rkz3 = rlx1*rly2 - rlx2*rly1
        rk   = rkx3*rlx3 + rky3*rly3 + rkz3*rlz3
        rkx3 = tpi*rkx3/rk
        rky3 = tpi*rky3/rk
        rkz3 = tpi*rkz3/rk
!=============================================================================
        hx=rprimd(1,1)/real(ngx)
        hy=rprimd(2,2)/real(ngy)
        hz=rprimd(3,3)/real(ngz)

      open(26,file="crot.dat",form='unformatted')
      open(29,file="overlap_abs.dat",form='formatted')
      open(30,file="singularvalues.dat",form='formatted')
      open(31,file="lo.status",form='formatted')
      open(6661,file="lo_center.dat",form='formatted')
      open(6662,file="lo_loclength.dat",form='formatted')

       rewind(26)
       rewind(29)
       rewind(30)
       rewind(6661)
       rewind(6662)

       write(29,*)' in the overlap matrix each row correspond to each orbital'
       write(29,*)' described below orbdetails variable in the input file.'
       write(29,*)' '

! diverting for only WLO construction :
       if(njob .ne. 0) go to 8765

!=====================================================
! Box in which the Gaussian has appreciable magnitude
         nxmax = -10000000
         nymax = -10000000
         nzmax = -10000000
         nxmin = 10000000
         nymin = 10000000
         nzmin = 10000000
         do k = 1,ngz
         do j = 1,ngy
         do i = 1,ngx

           r1 = (real(i)-(1.0-mod(real(ngx),2.0)*0.5) )/real(ngx)
           r2 = (real(j)-(1.0-mod(real(ngy),2.0)*0.5) )/real(ngy)
           r3 = (real(k)-(1.0-mod(real(ngz),2.0)*0.5) )/real(ngz)

           x=rlx*(-0.5)+rlx1*r1+rlx2*r2+rlx3*r3
           y=rly*(-0.5)+rly1*r1+rly2*r2+rly3*r3
           z=rlz*(-0.5)+rlz1*r1+rlz2*r2+rlz3*r3
           envelop=exp(-alpha*x*x)*exp(-alpha*y*y)*exp(-alpha*z*z)
           if(envelop.gt.0.000001)then
             if(i.gt.nxmax)nxmax = i
             if(i.lt.nxmin)nxmin = i
             if(j.gt.nymax)nymax = j
             if(j.lt.nymin)nymin = j
             if(k.gt.nzmax)nzmax = k
             if(k.lt.nzmin)nzmin = k
           end if

        end do
        end do
        end do
        nxw = nxmax - nxmin + 1
        nyw = nymax - nymin + 1
        nzw = nzmax - nzmin + 1

        nxhw = (nxw + anint(mod(real(nxw),2.0)) )/2
        nyhw = (nyw + anint(mod(real(nyw),2.0)) )/2
        nzhw = (nzw + anint(mod(real(nzw),2.0)) )/2

        print*,'Box width outside which localized template is set to zero:'
        print*,nxw,' x ',nyw,' x ',nzw
    !stop
!##########################################################################
        do csppol = 1,nsppol
         cspinor = 1 ! only scalar wave functions.

       do ikpt_z= 1,nkz
       do ikpt_y= 1,nky
       do ikpt_x= 1,nkx

        rewind(31)
! ---------------------------------------------
! k index for ascending order of wave vector k
!----------------------------------------------
        call kptindex(shiftx,shifty,shiftz,ikpt_x,ikpt_y,ikpt_z,nkx,nky,nkz,&
     &                     i,j,k)
        ckpt=i + (j - 1)*nkx + (k - 1)*nkx*nky
        write(31, '(a)' )" Doing : Overlap matrix and its singular value decomposition "
        write(31, '(a,3i4)' )" kpoint index: ",ikpt_x,ikpt_y,ikpt_z
        write(31, '(a,3i4)' )' folded to:',i,j,k

        write(*, '(a)' )" Doing : Overlap matrix and its singular value decomposition "
        write(*, '(a,3i4)' )" kpoint index: ",ikpt_x,ikpt_y,ikpt_z
        write(*, '(a,3i4)' )' folded to:',i,j,k

        rkpt1 = kpt(1,ckpt)
        rkpt2 = kpt(2,ckpt)
        rkpt3 = kpt(3,ckpt)

        if(ikpt_x .eq. nkx+1)rkpt1 = rkpt1 + 1.0
        write(31,*)' kpt:',rkpt1,rkpt2,rkpt3

        rkx = rkpt1*rkx1 + rkpt2*rkx2 + rkpt3*rkx3
        rky = rkpt1*rky1 + rkpt2*rky2 + rkpt3*rky3
        rkz = rkpt1*rkz1 + rkpt2*rkz2 + rkpt3*rkz3

        ncellx = 3
        ncelly = 3  ! we fix 3 unitcells only for the inverse f transform.
        ncellz = 3

        if(nkx .eq.1)ncellx = 1
        if(nky .eq.1)ncelly = 1  ! for molecules or low dimensional systems
        if(nkz .eq.1)ncellz = 1

!##############################################################
       do iband2 =  1, nbands  ! FIRST BAND INDEX STARTS HERE
! ---------------------------------------------
!        print*,'iband:',iband2
! real space grid :

         rnorm = 0.0
         corb = 0.0
        do icellx = 1,ncellx
        do icelly = 1,ncelly
        do icellz = 1,ncellz

          rn1 = real(icellx-(ncellx + anint(mod(real(ncellx),2.0)) )/2)
          rn2 = real(icelly-(ncelly + anint(mod(real(ncelly),2.0)) )/2)
          rn3 = real(icellz-(ncellz + anint(mod(real(ncellz),2.0)) )/2)

          munitx = (ngx+anint(mod(real(ngx),2.0)))/2 +
                              anint(mod(real(ngx)+1.0,2.0))
          munity = (ngy+anint(mod(real(ngy),2.0)))/2 +
                              anint(mod(real(ngy)+1.0,2.0))
          munitz = (ngz+anint(mod(real(ngz),2.0)))/2 +
                              anint(mod(real(ngz)+1.0,2.0))

          re1 = rcoordall(1,iband2,csppol)
          re2 = rcoordall(2,iband2,csppol)
          re3 = rcoordall(3,iband2,csppol)
          munitx = munitx + anint(re1*real(ngx))
          munity = munity + anint(re2*real(ngy))
          munitz = munitz + anint(re3*real(ngz))

          nstartx = munitx + anint(rn1)*ngx - nxhw
          nstopx  = munitx + anint(rn1)*ngx + nxhw

          nstarty = munity + anint(rn2)*ngy - nyhw
          nstopy  = munity + anint(rn2)*ngy + nyhw

          nstartz = munitz + anint(rn3)*ngz - nzhw
          nstopz  = munitz + anint(rn3)*ngz + nzhw

          if(nstartx.le.1)nstartx = 1
          if(nstartx.gt.ngx)nstartx = ngx
          if(nstarty.le.1)nstarty = 1
          if(nstarty.gt.ngy)nstarty = ngy
          if(nstartz.le.1)nstartz = 1
          if(nstartz.gt.ngz)nstartz = ngz

          if(nstopx.le.1)nstopx = 1
          if(nstopx.gt.ngx)nstopx = ngx
          if(nstopy.le.1)nstopy = 1
          if(nstopy.gt.ngy)nstopy = ngy
          if(nstopz.le.1)nstopz = 1
          if(nstopz.gt.ngz)nstopz = ngz

         do k = 1,ngz
         do j = 1,ngy
         do i = 1,ngx

           if(i.ge.nstartx.and.i.le.nstopx.and.&
      &       j.ge.nstarty.and.j.le.nstopy.and.&
      &       k.ge.nstartz.and.k.le.nstopz)then

           r1 = (real(i)-(1.0-mod(real(ngx),2.0)*0.5) )/real(ngx)
           r2 = (real(j)-(1.0-mod(real(ngy),2.0)*0.5) )/real(ngy)
           r3 = (real(k)-(1.0-mod(real(ngz),2.0)*0.5) )/real(ngz)

        x1=rlx*(-0.5)+rlx1*r1+rlx2*r2+rlx3*r3
        y1=rly*(-0.5)+rly1*r1+rly2*r2+rly3*r3
        z1=rlz*(-0.5)+rlz1*r1+rlz2*r2+rlz3*r3
! ---------------------------------------------
! for each atom :
          icountorb = 0
          do iatom = 1,matom

        rcoordx=rcoord(1,iatom)*rlx1+rcoord(2,iatom)*rlx2+rcoord(3,iatom)*rlx3
        rcoordy=rcoord(1,iatom)*rly1+rcoord(2,iatom)*rly2+rcoord(3,iatom)*rly3
        rcoordz=rcoord(1,iatom)*rlz1+rcoord(2,iatom)*rlz2+rcoord(3,iatom)*rlz3
! ---------------------------------------------
! for different orbitals at each atom :
           do iorb = 1,numorb(csppol,iatom)
            icountorb = icountorb + 1

!          *************************************************
            if(iband2 .eq. icountorb ) then
!          *************************************************

            i1 = ltypeorb(icountorb,csppol)
            thetax = rpi*rtheta(icountorb,1,csppol)/180.0
            thetay = rpi*rtheta(icountorb,2,csppol)/180.0
            thetaz = rpi*rtheta(icountorb,3,csppol)/180.0
! ---------------------------------------------
        corbital = 0.0

          x=rlx*(-0.5)+rlx1*r1+rlx2*r2+rlx3*r3-rcoordx-rn1*rlx1-rn2*rlx2-rn3*rlx3
          y=rly*(-0.5)+rly1*r1+rly2*r2+rly3*r3-rcoordy-rn1*rly1-rn2*rly2-rn3*rly3
          z=rlz*(-0.5)+rlz1*r1+rlz2*r2+rlz3*r3-rcoordz-rn1*rlz1-rn2*rlz2-rn3*rlz3

!-------------------------------------------
!rotating the orbitals as desired :
!-----------------------------------------
          rsin1=sin(thetax)
          rcos1=cos(thetax)

          rsin2=sin(thetay)
          rcos2=cos(thetay)

          rsin3=sin(thetaz)
          rcos3=cos(thetaz)
!-------------------------------------------
!        xmat  -->  clockwise rotation wrt x
!-------------------------------------------
      if(abs(thetax) .ne. 0.0)then
          xmat(1,1)=1.0
          xmat(1,2)=0.0
          xmat(1,3)=0.0

          xmat(2,1)=0.0
          xmat(2,2)=rcos1
          xmat(2,3)=-rsin1

          xmat(3,1)=0.0
          xmat(3,2)=rsin1
          xmat(3,3)=rcos1

          vec(1) = x ; vec(2) = y ; vec(3) = z
          vec1 = matmul(xmat,vec)
          x = vec1(1) ; y = vec1(2) ; z = vec1(3)
      end if
!-------------------------------------------
!        ymat  -->  clockwise rotation wrt y
!-------------------------------------------
      if(abs(thetay) .ne. 0.0)then
          ymat(1,1)=rcos2
          ymat(1,2)=0.0
          ymat(1,3)=-rsin2

          ymat(2,1)=0.0
          ymat(2,2)=1.0
          ymat(2,3)=0.0

          ymat(3,1)=rsin2
          ymat(3,2)=0.0
          ymat(3,3)=rcos2

          vec(1) = x ; vec(2) = y ; vec(3) = z
          vec1 = matmul(ymat,vec)
          x = vec1(1) ; y = vec1(2) ; z = vec1(3)
      end if
!-------------------------------------------
!        zmat  -->  clockwise rotation wrt z
!-------------------------------------------
      if(abs(thetaz) .ne. 0.0)then
          zmat(1,1)=rcos3
          zmat(1,2)=-rsin3
          zmat(1,3)=0.0

          zmat(2,1)=rsin3
          zmat(2,2)=rcos3
          zmat(2,3)=0.0

          zmat(3,1)=0.0
          zmat(3,2)=0.0
          zmat(3,3)=1.0

          vec(1) = x ; vec(2) = y ; vec(3) = z
          vec1 = matmul(zmat,vec)

          x = vec1(1) ; y = vec1(2) ; z = vec1(3)
      end if
!-------------------------------------------
! gaussian envelop :
      envelop=exp(-alpha*x*x)*exp(-alpha*y*y)*exp(-alpha*z*z)
      if((x*x + y*y + z*z).eq.0.0)then
        rsintheta=0.0
        rcostheta=0.0
        else
        rsintheta=sqrt(z*z + y*y)/sqrt(x*x + y*y + z*z)
        rcostheta=x/sqrt(x*x + y*y + z*z)
      end if

      if((y*y + z*z).eq.0.0)then
        rsinphy  =0.0
        rcosphy  =0.0
        else
        rsinphy  =z/sqrt(z*z + y*y)
        rcosphy  =y/sqrt(z*z + y*y)
      end if

      rrsinphy2 =2.0*rsinphy*rcosphy
      rrcosphy2 =rcosphy*rcosphy-rsinphy*rsinphy
!------------------------------------------------------
! sherical harmonics :
!------------------------------------------------------
! s
      if(i1.eq.1)orbital= envelop*sqrt(1.0/(4.0*rpi))

      if(i1.eq.2.or.i1.eq.3)then
      cmp1=envelop*rsintheta*(rcosphy+cmplx(0.0,1.0)*rsinphy)*sqrt(3.0/(8.0*rpi))
      cmp2=envelop*(-1.0)*rsintheta*(rcosphy-cmplx(0.0,1.0)*rsinphy)*sqrt(3.0/(8.0*rpi))
! py
      if(i1.eq.2)orbital= (1.0/sqrt(2.0))*(cmp1-cmp2)
! pz

      if(i1.eq.3)orbital= (-cmplx(0.0,1.0)/sqrt(2.0))*(cmp1+ cmp2)
      end if
! px
     if(i1.eq.4)then
      cmp1=envelop*rcostheta*sqrt(3.0/(4.0*rpi))
      orbital=  cmp1
     end if

     if(i1.eq.5.or.i1.eq.6)then
      cmp1=envelop*(-1.0)*sqrt(15.0/(8.0*rpi))*rsintheta*rcostheta*(rcosphy+cmplx(0.0,1.0)*rsinphy)
      cmp2=envelop*sqrt(15.0/(8.0*rpi))*rsintheta*rcostheta*(rcosphy-cmplx(0.0,1.0)*rsinphy)
! dyz
      if(i1.eq.5)orbital= (cmplx(0.0,1.0)/sqrt(2.0))*(cmp1+cmp2)
! dxz
      if(i1.eq.6)orbital= (-1.0/sqrt(2.0))*(cmp1-cmp2)
     end if

     if(i1.eq.7.or.i1.eq.8)then
      cmp1=envelop*sqrt(15.0/(32.0*rpi))*rsintheta*rsintheta*(rrcosphy2+cmplx(0.0,1.0)*rrsinphy2)
      cmp2=envelop*sqrt(15.0/(32.0*rpi))*rsintheta*rsintheta*(rrcosphy2-cmplx(0.0,1.0)*rrsinphy2)
! dxy
      if(i1.eq.7)orbital= (-cmplx(0.0,1.0)/sqrt(2.0))*(cmp1 - cmp2)
! dx2-y2
      if(i1.eq.8)orbital= (1.0/sqrt(2.0))*(cmp1 + cmp2)
      end if

! dz2
      if(i1.eq.9)then
      cmp1=envelop*sqrt(5.0/(16.0*rpi))*(3.0*rcostheta*rcostheta-1.0)
      orbital= cmp1
      end if
!--------------------------------------------------------------------------
         corbital = corbital  +  orbital*&
&           exp(cmplx(0.0,1.0)*(rkx*(rn1*rlx1+rn2*rlx2+rn3*rlx3 - x1)+&
&                               rky*(rn1*rly1+rn2*rly2+rn3*rly3 - y1)+&
&                               rkz*(rn1*rlz1+rn2*rlz2+rn3*rlz3 - z1) ))
        !   corb(i,j,k)= corbital
           corb(i,j,k)= corb(i,j,k) + corbital
        !   rnorm = rnorm + conjg(corbital)*corbital

 !          *************************************************
             end if
!          *************************************************

           end do  ! orbital type
          end do ! atom
       if(icountorb.ne.nbands)then
         write(*,*)icountorb,' -vs- ',nbands
         write(*,*)'total number of orbitals and nband do not match'
         stop
       end if
         !-----------new
         end if
           !-------

! normalization :
          if(icellx.eq.ncellx.and.icelly.eq.ncelly.and.icellz.eq.ncellz)then
              rnorm = rnorm + conjg(corb(i,j,k))*corb(i,j,k)
           end if

        end do  ! real space grid
        end do
        end do

         end do
         end do  ! sum over neighbouring unitcell ends
         end do

! normalization :
       ! print*,'norm:',rnorm
        corb = corb / sqrt(rnorm)

!#########################################

       do iband1 =  1, nbands
       ! print*,'   overlap beqween', iband1,' and ',iband2
! reading wave function from "_wfk" file and calculating the overlap :

       call overlap_wf(corb,e_kpt,exchn2n3,csppol,iband1,ckpt,&
     & ecut,headform,istwfk,kpt,nband,nbands,nkpt,npwarr,&
     & nr1,nr2,nr3,nspinor,nsppol,rprimd,cmp2)

       coverlap(iband2,iband1) = cmp2

       end do  ! iband1
!#########################################
       end do  ! iband2
!#####################################################

! singular value decomposition :
!-----------------------------------
       write(29,*)" kpoint:",ikpt_x,ikpt_y,ikpt_z
       write(30,*)" kpoint:",ikpt_x,ikpt_y,ikpt_z
       write(29,*)"  "
       write(30,*)"  "

       do i = 1,nbands
       write(29,"(20(f12.5))")(abs(coverlap(i,j)),j=1,nbands)
       end do
       crhsvd = 'a'
       call zgesvd(crhsvd,crhsvd,nbands,nbands,coverlap,nbands,s,cu,nbands,&
                          cvt,nbands,csvd, ndimcsvd, rwork, info )

               cm5 = 0.0
               do i = 1,nbands
                 if(abs(s(i)).gt.0.0000001)cm5(i,i)= 1.0
               end do
               write(30,"(10f8.4)")(s(j),j=1,nbands)
               write(31,*)' singular values:'
               write(31,"(10f8.4)")(s(j),j=1,nbands)
               write(31,*)' note: if you find some of the singular values consitently'
               write(31,*)'       very close to zero (0.00xxx) please stop the'
               write(31,*)'       calculation and look at the overlap_abs.dat file.'
               write(31,*)'       exclude orbitals (from under the orbdetails variable '
               write(31,*)'       in the input file) corresponding to predominantly'
               write(31,*)'       zero rows of overlap matrix. '

               cm1=matmul(cu,cvt)
               cm2=conjg(cm1)

             !  crotl(ikpt_x,ikpt_y,ikpt_z,1:nbands,1:nbands) = cm2(1:nbands,1:nbands)

          write(26)((cm2(i1,i2),i1=1,nbands),i2 = 1,nbands)

!--------------------------------------------------------------------------
       end do  ! ikpt_x
       end do  ! ikpt_y  ! ikpt loop ends here
       end do  ! ikpt_z

       end do  ! csppol
       close(26)

8765    continue  ! if jobtype is only WLO construction

! ##########################################################################
! #                                                                        #
! #                   Wannier function construction                        #
! #                                                                        #
! ##########################################################################
      normtot = 0.0
      do  iwan = 1,nwantot  ! BAND INDEX, ie, WANNIER INDEX STARTS HERE

       open(26,file="crot.dat",form='unformatted')
       rewind(26)

        iband1 = nwanindex(iwan)
        cwannier = 0.0
        rwanmax = -1000.0
! ======================================================
      do csppol = 1,nsppol
         cspinor = 1 ! only scalar wave functions.

        write(*,*)' reading'
        do ikpt_z = 1,nkz
        do ikpt_y = 1,nky
        do ikpt_x = 1,nkx
           rewind(31)
           write(31,*)' Doing : WLO construction.'
           write(31,*)' WLO index:',iband1

           write(*,*)' Doing : WLO construction.'
           write(*,*)' WLO index:',iband1

           read(26)((cm2(i1,i2),i1=1,nbands),i2 = 1,nbands)

! ---------------------------------------------
! k index for ascending order of wave vector k
!----------------------------------------------
        call kptindex(shiftx,shifty,shiftz,ikpt_x,ikpt_y,ikpt_z,nkx,nky,nkz,i,j,k)
        ckpt=i + (j - 1)*nkx + (k - 1)*nkx*nky
        write(31, '(a,3i4)' )" kpoint index : ",ikpt_x,ikpt_y,ikpt_z
        write(31, '(a,3i4)' )" kpoint index : ",ikpt_x,ikpt_y,ikpt_z
        write(31, '(a,3i4)' )' folded to:',i,j,k

        write(*,'(a,3i4)')" kpoint index : ",ikpt_x,ikpt_y,ikpt_z
        rkpt1 = kpt(1,ckpt)
        rkpt2 = kpt(2,ckpt)
        rkpt3 = kpt(3,ckpt)

        write(31,*)'kpt:',rkpt1,rkpt2,rkpt3

        rkx = rkpt1*rkx1 + rkpt2*rkx2 + rkpt3*rkx3
        rky = rkpt1*rky1 + rkpt2*rky2 + rkpt3*rky3
        rkz = rkpt1*rkz1 + rkpt2*rkz2 + rkpt3*rkz3
!----------------------------------
! rotating the energy eigenstates :
!----------------------------------
      cbfns1  = cmplx(0.0_dp,0.0_dp)
      do iband2 = 1,nbands

       call wfread(cbfns,e_kpt,exchn2n3,csppol,iband2,ckpt,&
     & ecut,headform,istwfk,kpt,nband,nbands,nkpt,npwarr,&
     & nr1,nr2,nr3,nspinor,nsppol,rprimd)

       if(insmet.eq.2)then
         roccup=(e_kpt(iband2)-efermi)/tsmear ! rkt
         roccup = 1.0_dp/(1.0_dp+exp(roccup))
        else
          roccup = 1.0
       end if
         do igz = 1,ngz
          do igy = 1,ngy
           do igx = 1,ngx

       cbfns1(igx,igy,igz)=cbfns1(igx,igy,igz)&
     &      +cm2(iband1,iband2)*cbfns(igx,igy,igz)*sqrt(roccup)

           end do
          end do
         end do
      end do ! iband2

!------------------
! wf construction:
!------------------
        if(rcellx+rcelly+rcellz .ne. 0.0) then  ! 0 is default

          munitx = (ngx*ncellx1+anint(mod(real(ngx*ncellx1),2.0)))/2 +
                              anint(mod(real(ngx*ncellx1)+1.0,2.0))
          munity = (ngy*ncelly1+anint(mod(real(ngy*ncelly1),2.0)))/2 +
                              anint(mod(real(ngy*ncelly1)+1.0,2.0))
          munitz = (ngz*ncellz1+anint(mod(real(ngz*ncellz1),2.0)))/2 +
                              anint(mod(real(ngz*ncellz1)+1.0,2.0))

          re1 = rcoordall(1,iband1,csppol)
          re2 = rcoordall(2,iband1,csppol)
          re3 = rcoordall(3,iband1,csppol)
          munitx = munitx + anint(re1*real(ngx))
          munity = munity + anint(re2*real(ngy))
          munitz = munitz + anint(re3*real(ngz))

          nstartx = munitx  - anint(rcellx*0.5*real(ngx))
          nstopx  = munitx  + anint(rcellx*0.5*real(ngx))

          nstarty = munity  - anint(rcelly*0.5*real(ngy))
          nstopy  = munity  + anint(rcelly*0.5*real(ngy))

          nstartz = munitz  - anint(rcellz*0.5*real(ngz))
          nstopz  = munitz  + anint(rcellz*0.5*real(ngz))

          if(nstartx.le.1)nstartx = 1
          if(nstartx.gt.ngx*ncellx1)nstartx = ngx*ncellx1
          if(nstarty.le.1)nstarty = 1
          if(nstarty.gt.ngy*ncelly1)nstarty = ngy*ncelly1
          if(nstartz.le.1)nstartz = 1
          if(nstartz.gt.ngz*ncellz1)nstartz = ngz*ncellz1

          if(nstopx.le.1)nstopx = 1
          if(nstopx.gt.ngx*ncellx1)nstopx = ngx*ncellx1
          if(nstopy.le.1)nstopy = 1
          if(nstopy.gt.ngy*ncelly1)nstopy = ngy*ncelly1
          if(nstopz.le.1)nstopz = 1
          if(nstopz.gt.ngz*ncellz1)nstopz = ngz*ncellz1

         else

          nstartx = 1
          nstopx  = ngx*ncellx1

          nstarty = 1
          nstopy  = ngy*ncelly1

          nstartz = 1
          nstopz  = ngz*ncellz1

          end if

       do icellz = 1,ncellz1
         do icelly = 1,ncelly1
          do icellx = 1,ncellx1


         do igz = 1,ngz
         do igy = 1,ngy
         do igx = 1,ngx

           l = ngx*(icellx-1) + igx
           m = ngy*(icelly-1) + igy
           n = ngz*(icellz-1) + igz

           if(l.ge.nstartx.and.l.le.nstopx.and.&
      &       m.ge.nstarty.and.m.le.nstopy.and.&
      &       n.ge.nstartz.and.n.le.nstopz)then

           if(mod(real(ncellx1),2.0).eq.1)then
                   ii1 = igx
           else
                 if(mod(real(ngx),2.0).eq.1)then
                   ii1 = igx+(ngx-1)/2
                   if(ii1.gt.ngx)ii1 = ii1 - ngx
                  else
                   ii1 = igx+ngx/2
                   if(ii1.gt.ngx)ii1 = ii1 - ngx
                 end if
           end if

           if(mod(real(ncelly1),2.0).eq.1)then
                   ii2 = igy
           else
                 if(mod(real(ngy),2.0).eq.1)then
                   ii2 = igy+(ngy-1)/2
                   if(ii2.gt.ngy)ii2 = ii2 - ngy
                  else
                   ii2 = igy+ngy/2
                   if(ii2.gt.ngy)ii2 = ii2 - ngy
                 end if
           end if

           if(mod(real(ncellz1),2.0).eq.1)then
                   ii3 = igz
           else
                 if(mod(real(ngz),2.0).eq.1)then
                   ii3 = igz+(ngz-1)/2
                   if(ii3.gt.ngz)ii3 = ii3 - ngz
                  else
                   ii3 = igz+ngz/2
                   if(ii3.gt.ngz)ii3 = ii3 - ngz
                 end if
           end if

        r1 = (real(l)-(1.0-mod(real(ngx*ncellx1),2.0)*0.5) )/real(ngx)
        r2 = (real(m)-(1.0-mod(real(ngy*ncelly1),2.0)*0.5) )/real(ngy)
        r3 = (real(n)-(1.0-mod(real(ngz*ncellz1),2.0)*0.5) )/real(ngz)

        x = rlx_large*(-0.5)+rlx1*r1+rlx2*r2+rlx3*r3
        y = rly_large*(-0.5)+rly1*r1+rly2*r2+rly3*r3
        z = rlz_large*(-0.5)+rlz1*r1+rlz2*r2+rlz3*r3

           cwannier(l,m,n)=cwannier(l,m,n)+&
&             cbfns1(ii1,ii2,ii3)*exp(cmplx(0.0,1.0)*(rkx*x+rky*y+rkz*z))

           if(ikpt_x.eq.nkx.and.ikpt_y.eq.nky.and.ikpt_z.eq.nkz)then
             if(abs(cwannier(l,m,n)).gt.rwanmax)then
               rwanmax = abs(cwannier(l,m,n))
               nwanmax1 = l
               nwanmax2 = m
               nwanmax3 = n
             end if
           end if
!       ---------------new
          end if
    !------------
          end do
          end do
         end do
         end do
        end do
        end do

       end do  ! ikpt_x
       end do  ! ikpt_y  ! ikpt loop ends here
       end do  ! ikpt_z

!--------------------------------------------------------------------------
! processing of wlo :
!-------------------
!normalization  :
      cwannier=cwannier/real(nkx*nky*nkz)

        re1 = 0.0_dp
        re2 = 0.0_dp
        re3 = 0.0_dp
        re4 = 0.0_dp
        re5 = 0.0_dp
        re6 = 0.0_dp
        re7 = 0.0_dp

             do k = nstartz,nstopz
              do j = nstarty,nstopy
               do i = nstartx,nstopx

        r1 = (real(i)-(1.0-mod(real(ngx*ncellx1),2.0)*0.5) )/real(ngx)
        r2 = (real(j)-(1.0-mod(real(ngy*ncelly1),2.0)*0.5) )/real(ngy)
        r3 = (real(k)-(1.0-mod(real(ngz*ncellz1),2.0)*0.5) )/real(ngz)

        x = rlx_large*(-0.5)+rlx1*r1+rlx2*r2+rlx3*r3
        y = rly_large*(-0.5)+rly1*r1+rly2*r2+rly3*r3
        z = rlz_large*(-0.5)+rlz1*r1+rlz2*r2+rlz3*r3

           csum=cwannier(i,j,k)*conjg(cwannier(i,j,k))
        re1 = re1 + x*csum
        re2 = re2 + y*csum
        re3 = re3 + z*csum
        re4 = re4 + x*x*csum
        re5 = re5 + y*y*csum
        re6 = re6 + z*z*csum
        re7 = re7 + csum

            end do
          end do
        end do
    if(re7.gt.0.0000001)then
        re1 = re1/re7  ; rxcm = re1/re7
        re2 = re2/re7  ; rycm = re2/re7
        re3 = re3/re7  ; rzcm = re3/re7
        re4 = re4/re7
        re5 = re5/re7
        re6 = re6/re7
    end if
        normtot = normtot + re7
        write(*,*)" ========================"
        write(*,*)" WLO : ",iband1,'csppol:',csppol

        write(*,*)' normalization:  ',re7
        write(*,*)' qo center x',re1
        write(*,*)' qo center y',re2
        write(*,*)' qo center z',re3
        write(*,*)"-----------------------------"
        write(*,*)' loclength_x*x',(re4-re1*re1)
        write(*,*)' loclength_y*y',(re5-re2*re2)
        write(*,*)' loclength_z*z',(re6-re3*re3)
        write(6662,"(i4,4(1x,f8.5),i2)")iband1,re4-re1*re1,re5-re2*re2,re6-re3*re3,re7,csppol
        write(6661,"(i4,4(1x,f8.5),i2)")iband1,re1,re2,re3,re7,csppol
        if(iwan.eq.nwantot)write(6661,*)' total norm:',normtot
        if(iwan.eq.nwantot)write(6662,*)' total norm:',normtot

!---------------------------------
!   Plotting :
          munitx = (ngx*ncellx1+anint(mod(real(ngx*ncellx1),2.0)))/2 +
                             anint(mod(real(ngx*ncellx1)+1.0,2.0))
          munity = (ngy*ncelly1+anint(mod(real(ngy*ncelly1),2.0)))/2 +
                             anint(mod(real(ngy*ncelly1)+1.0,2.0))
          munitz = (ngz*ncellz1+anint(mod(real(ngz*ncellz1),2.0)))/2 +
                             anint(mod(real(ngz*ncellz1)+1.0,2.0))

          write(6, '(a,3i5)' ) ' mid points :',munitx,munity,munitz
          re1 = rcoordall(1,iband1,csppol)
          re2 = rcoordall(2,iband1,csppol)
          re3 = rcoordall(3,iband1,csppol)
          munitx = munitx + anint(re1*real(ngx))
          munity = munity + anint(re2*real(ngy))
          munitz = munitz + anint(re3*real(ngz))

           write(6, '(a,3i5)' ) ' center sites :',munitx,munity,munitz

           if(icellplotx.gt.1)then
            nstartx1 = munitx - ngx + 1
            nstopx1  = icellplotx*ngx + nstartx1 - 1
            if(nstartx1.le.1)then
              nstartx1 = 1
              nstopx1  = icellplotx*ngx
            end if
            if(nstopx1.ge.ncellx1*ngx)then
               nstartx1 = (ncellx1 -2)*ngx + 1
               nstopx1  = ncellx1*ngx
            end if
           else
            nstartx1 = 1
            nstopx1  = icellplotx*ngx
           end if

           if(icellploty.gt.1)then
            nstarty1 = munity - ngy + 1
            nstopy1  = icellploty*ngy + nstarty1 - 1
            if(nstarty1.le.1)then
              nstarty1 = 1
              nstopy1  = icellploty*ngy
            end if
            if(nstopy1.ge.ncelly1*ngy)then
              nstarty1 = (ncelly1 -2)*ngy + 1
              nstopy1  = ncelly1*ngy
            end if
           else
            nstarty1 = 1
            nstopy1 = icellploty*ngy
          end if

           if(icellplotz.gt.1)then
            nstartz1 = munitz - ngz + 1
            nstopz1  = icellplotz*ngz + nstartz1 - 1
            if(nstartz1.le.1)then
               nstartz1 = 1
               nstopz1  = icellplotz*ngz
            end if
            if(nstopz1.ge.ncellz1*ngz)then
               nstartz1 = (ncellz1 -2)*ngz + 1
               nstopz1  = ncellz1*ngz
            end if
           else
            nstartz1 = 1
            nstopz1  = icellplotz*ngz
          end if

           nstartx = max(nstartx,nstartx1)
           nstarty = max(nstarty,nstarty1)
           nstartz = max(nstartz,nstartz1)

           nstopx = min(nstopx1,nstopx)
           nstopy = min(nstopy1,nstopy)
           nstopz = min(nstopz1,nstopz)

           write(6, '(a)' )' plotting range of LO :'
           write(6, '(a,2i5)' )' plot x :',nstartx,nstopx
           write(6, '(a,2i5)' )' plot y :',nstarty,nstopy
           write(6, '(a,2i5)' )' plot z :',nstartz,nstopz

       write(6, '(a,3i5)' )' ng :',ngx,ngy,ngz

       write(6, '(a,3i5)' )'     ',shape(cwannier)

      rewind(10000+iband1*10+csppol)
      if(iwrite.eq.1)rewind(20000+iband1*10+csppol)

! plotting in xcrysden format :

      write(10000+iband1*10+csppol,*)"CRYSTAL"
      write(10000+iband1*10+csppol,*)"PRIMVEC"
      write(10000+iband1*10+csppol,"(3f10.5)") rlx1  , rly1 ,rlz1
      write(10000+iband1*10+csppol,"(3f10.5)") rlx2,  rly2  , rlz2
      write(10000+iband1*10+csppol,"(3f10.5)") rlx3, rly3,  rlz3
      write(10000+iband1*10+csppol,*)"CONVVEC"
      write(10000+iband1*10+csppol,"(3f10.5)")  rlx1  , rly1 ,rlz1
      write(10000+iband1*10+csppol,"(3f10.5)") rlx2,  rly2  , rlz2
      write(10000+iband1*10+csppol,"(3f10.5)") rlx3, rly3,  rlz3
      write(10000+iband1*10+csppol,*)"PRIMCOORD"
      write(10000+iband1*10+csppol,*)natom, 1

       do i = 1,natom
     write(10000+iband1*10+csppol,"(i3,3f10.5)")int(znucl(typat(i))),(xcart(j,i),j=1,3)
       end do
      write(10000+iband1*10+csppol,*)"BEGIN_BLOCK_DATAGRID_3D"
      write(10000+iband1*10+csppol,*)"WANNIER_CENTER_3D"
      write(10000+iband1*10+csppol,*)"DATAGRID_3D_WANNIER"

      write(10000+iband1*10+csppol,*)nstopx-nstartx+1,nstopy-nstarty+1,nstopz-nstartz+1

       re4 = real(nstopx-nstartx+1-1)/real(ngx)
       re5 = real(nstopy-nstarty+1-1)/real(ngy)
       re6 = real(nstopz-nstartz+1-1)/real(ngz)

       r1 = (real(nstartx)-(1.0-mod(real(ngx*ncellx1),2.0)*0.5) )/real(ngx)
       r2 = (real(nstarty)-(1.0-mod(real(ngy*ncelly1),2.0)*0.5) )/real(ngy)
       r3 = (real(nstartz)-(1.0-mod(real(ngz*ncellz1),2.0)*0.5) )/real(ngz)

       x = rlx_large*(-0.5)+rlx1*r1+rlx2*r2+rlx3*r3
       y = rly_large*(-0.5)+rly1*r1+rly2*r2+rly3*r3
       z = rlz_large*(-0.5)+rlz1*r1+rlz2*r2+rlz3*r3

      write(10000+iband1*10+csppol,"(3f10.5)")x,y,z
      if(iwrite.eq.1)write(20000+iband1*10+csppol,*)x,y,z

      write(10000+iband1*10+csppol,"(3f10.5)")rlx1*re4,rly1*re4,rlz1*re4
      write(10000+iband1*10+csppol,"(3f10.5)")rlx2*re5,rly2*re5,rlz2*re5
      write(10000+iband1*10+csppol,"(3f10.5)")rlx3*re6,rly3*re6,rlz3*re6

      if(iwrite.eq.1)then
         write(20000+iband1*10+csppol,*)nstartx,nstarty,nstartz
         write(20000+iband1*10+csppol,*)nstopx,nstopy,nstopz
         write(20000+iband1*10+csppol,"(3f10.5)")rlx1*re4,rly1*re4,rlz1*re4
         write(20000+iband1*10+csppol,"(3f10.5)")rlx2*re5,rly2*re5,rlz2*re5
         write(20000+iband1*10+csppol,"(3f10.5)")rlx3*re6,rly3*re6,rlz3*re6
      end if

!     Apparently, the determination of nwanmax1, nwanmax2 and nwanmax3 is no portable accross platforms'
      write(6, '(a,3i5)' )'-wanmax:',nwanmax1,nwanmax2,nwanmax3
      cmp1 = cwannier(nwanmax1,nwanmax2,nwanmax3)
      print*,cmp1
      cmp1 = abs(cmp1)/cmp1

             re1 = 0.0
             do igz1 = nstartz  ,  nstopz
             do igy1 = nstarty  ,  nstopy
             do igx1 = nstartx  ,  nstopx
               write(10000+iband1*10+csppol,"(1f10.5)")real(cwannier(igx1,igy1,igz1)*cmp1)
              if(iwrite.eq.1)write(20000+iband1*10+csppol,*)cwannier(igx1,igy1,igz1)
              re1 = re1 + abs(cwannier(igx1,igy1,igz1))**2.0
             end do
             end do
             end do
             write(*,*)'norm enclosed:',re1

      write(10000+iband1*10+csppol,*)"END_DATAGRID_3D"
      write(10000+iband1*10+csppol,*)"END_BLOCK_DATAGRID_3D"

!------------------------------------------------------
! X Y Z projection calculation
!------------------------------------------------------
!  goto 6776
       re1 = abs(rlx1)*real(ncellx1)+abs(rlx2)*real(ncelly1)+abs(rlx3)*real(ncellz1)
       re2 = abs(rly1)*real(ncellx1)+abs(rly2)*real(ncelly1)+abs(rly3)*real(ncellz1)
       re3 = abs(rlz1)*real(ncellx1)+abs(rlz2)*real(ncelly1)+abs(rlz3)*real(ncellz1)

       re4 = max(abs(rlx1)/real(ngx),abs(rlx2)/real(ngy),abs(rlx3)/real(ngz))
       re5 = max(abs(rly1)/real(ngx),abs(rly2)/real(ngy),abs(rly3)/real(ngz))
       re6 = max(abs(rlz1)/real(ngx),abs(rlz2)/real(ngy),abs(rlz3)/real(ngz))

       nxtot = anint(re1/re4)
       nytot = anint(re2/re5)
       nztot = anint(re3/re6)
       allocate(rprojx(nxtot),rprojy(nytot),rprojz(nztot))
       rprojx = 0.0
       rprojy = 0.0
       rprojz = 0.0

       re4 = min(rlx1*real(ncellx1),rlx2*real(ncelly1),rlx3*real(ncellz1))
       if(re4.ge.0.0)re4 = 0.0
       re5 = min(rly1*real(ncellx1),rly2*real(ncelly1),rly3*real(ncellz1))
       if(re5.ge.0.0)re5 = 0.0
       re6 = min(rlz1*real(ncellx1),rlz2*real(ncelly1),rlz3*real(ncellz1))
       if(re6.ge.0.0)re6 = 0.0

             do k = nstartz,nstopz
              do j = nstarty,nstopy
               do i = nstartx,nstopx

        r1 = (real(i)-(1.0-mod(real(ngx*ncellx1),2.0)*0.5) )/real(ngx)
        r2 = (real(j)-(1.0-mod(real(ngy*ncelly1),2.0)*0.5) )/real(ngy)
        r3 = (real(k)-(1.0-mod(real(ngz*ncellz1),2.0)*0.5) )/real(ngz)

        x = rlx1*r1+rlx2*r2+rlx3*r3
        y = rly1*r1+rly2*r2+rly3*r3
        z = rlz1*r1+rlz2*r2+rlz3*r3

        i1 = 1+anint((x - re4 )/(abs(re1)/real(nxtot)))
        i2 = 1+anint((y - re5 )/(abs(re2)/real(nytot)))
        i3 = 1+anint((z - re6 )/(abs(re3)/real(nztot)))

        if(i1.gt.nxtot)print*,'i1=',i1
        if(i2.gt.nytot)print*,'i2=',i2
        if(i3.gt.nztot)print*,'i3=',i3

        rprojx(i1)= rprojx(i1) + cwannier(i,j,k)*conjg(cwannier(i,j,k))
        rprojy(i2)= rprojy(i2) + cwannier(i,j,k)*conjg(cwannier(i,j,k))
        rprojz(i3)= rprojz(i3) + cwannier(i,j,k)*conjg(cwannier(i,j,k))

               end do
              end do
             end do

             rewind(510000+iband1*10+csppol)
             rewind(520000+iband1*10+csppol)
             rewind(530000+iband1*10+csppol)

             do i = 1,nxtot
              x =  -0.5*re1 +(abs(re1)/real(nxtot))*(real(i)-(1.0-mod(real(nxtot),2.0)*0.5))
              write(510000+iband1*10+csppol,"(10f12.6)")x,rprojx(i)
             end do
             do i = 1,nytot
              y =  -0.5*re2+(abs(re2)/real(nytot))*(real(i)-(1.0-mod(real(nytot),2.0)*0.5))
              write(520000+iband1*10+csppol,"(10f12.6)")y,rprojy(i)
             end do
             do i = 1,nztot
              z =  -0.5*re3+(abs(re3)/real(nztot))*(real(i)-(1.0-mod(real(nztot),2.0)*0.5))
              write(530000+iband1*10+csppol,"(10f12.6)")z,rprojz(i)
             end do

             close(510000+iband1*10+csppol)
             close(520000+iband1*10+csppol)
             close(530000+iband1*10+csppol)
        deallocate(rprojx,rprojy,rprojz)

!------------------------------------------------------------
       end do  ! csppol
       close(26)
!-------------------------------------------------------------
      end do ! iwan

       close(29)
       close(30)
       close(6661)
       close(6662)
       close(unit=31,status='delete')

      close(19)
      return
end subroutine localorb_s
!!***

Generated by  Doxygen 1.6.0   Back to index