!--------------------------------------------------------------------------------! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany! This file is part of FLEUR and available as free software under the conditions! of the MIT license as expressed in the LICENSE file in more detail.!--------------------------------------------------------------------------------

#include "cpp_double.h"#ifdef CPP_MPIinclude'mpif.h'integerierr(3)integercpu_indexintegerstt(MPI_STATUS_SIZE)#endifreal,allocatable::eigg(:)realkpoints(nkptd)!!! a and b coeff. constructed for each k-pointcomplex,allocatable::acof(:,:,:)complex,allocatable::bcof(:,:,:)complex,allocatable::ccof(:,:,:,:)!!! the parameters for the number of wfsinteger::nwfs!!! the potential in the spheres and the vacuumreal,allocatable::vr(:,:,:),vz(:,:,:),vrf(:,:,:,:)!!! auxiliary potentialscomplex,allocatable::vpw(:,:),vzxy(:,:,:,:)!!! bkpts dataintegernntot,ikpt_helpinteger,allocatable::gb(:,:,:),bpt(:,:)!!! radial wavefunctions in the muffin-tins and more ...real,allocatable::flo(:,:,:,:)real,allocatable::ff(:,:,:,:),gg(:,:,:,:)

ccyclebyspinsstarts! NO NON-COLLINEAR WHATSOEVER!do110jspin=1,jspd! cycle by spinsprint*,"spin=",jspinjsp_start=jspin;jsp_end=jspinc**************************************************************cforbzsym=.true.:determinemappingbetweenkptsandw90kptsc**************************************************************if(l_bzsym)thenl_file=.false.inquire(file='w90kpts',exist=l_file)IF(.NOT.l_file)CALLjuDFT_error+("w90kpts not found, needed if bzsym",calledby+="wann_plotw90")open(412,file='w90kpts',form='formatted')read(412,*)fullnkptscif(l_matrixmmn)thencallocate(kpoi(3,fullnkpts))cdok=1,fullnkptscread(412,*)kpoi(:,k)cenddocendifclose(412)print*,"fullnkpts=",fullnkptsIF(fullnkpts<=nkpts)CALLjuDFT_error("fullnkpts.le.nkpts"+,calledby="wann_plotw90")allocate(irreduc(fullnkpts),mapkoper(fullnkpts))l_file=.false.inquire(file='kptsmap',exist=l_file)IF(.NOT.l_file)CALLjuDFT_error+("kptsmap not found, needed if bzsym",calledby+="wann_plotw90")open(713,file='kptsmap')doi=1,fullnkptsread(713,*)kpt,irreduc(i),mapkoper(i)IF(kpt/=i)CALLjuDFT_error("kpt.ne.i",calledby+="wann_plotw90")print*,i,irreduc(i),mapkoper(i)enddoclose(713)IF(MAXVAL(irreduc(:))/=nkpts)CALLjuDFT_error+("max(irreduc(:))/=nkpts",calledby="wann_plotw90")elsefullnkpts=nkptscif(l_matrixmmn)thencallocate(kpoi(3,fullnkpts))copen(412,file='kpts',form='formatted')cread(412,*)cdok=1,fullnkptscread(412,*)kpoi(:,k)cenddocendifendifcif(l_matrixmmn)thencprint*,kpoi(:,:)cendifc**************************************************************if(l_byindex)num_bands=band_max-band_min+1if(l_byenergy)then!determine number of bandsl_file=.false.inquire(file=spin12(jspin)//'.eig',exist=l_file)

<ac,bc,u,ue)endif!preparations for vacuumc**************************************************************************c**************************************************************************nbmin=1nbmax=nslibdband:DOnbn=nbmin,nbmaxDOiz=0,grid(3)-1DOiy=0,grid(2)-1xloop:DOix=0,grid(1)-1posi=ix+1+iy*(grid(1))+iz*(grid(1))*(grid(2))point=zero+vec1*(ix+0.0)/(grid(1))+vec2*(iy+0.0)$/(grid(2))IF(.NOT.twodim)point=point+vec3*(iz+0.0)/(grid(3))!Check if the point is in MT-sphereif(film)thenfas=-bkpt(3)*bmat(3,3)*z1factor=cmplx(cos(fas),sin(fas))elsefactor=cmplx(1.0,0.0)endifii1=3ii2=3ii3=3IF(film.AND..NOT.odi%d1)ii3=0IF(odi%d1)THENii1=0;ii2=0ENDIFDOi1=-ii1,ii1DOi2=-ii2,ii2DOi3=-ii3,ii3pt=point+MATMUL(amat,(/i1,i2,i3/))na=0DOnt=1,ntypeDOnq=1,neq(nt)na=na+1s=SQRT(dot_PRODUCT(pos(:,na)-pt,pos(:,na)-pt))IF(s<rmsh(jri(nt),nt))THENCALLwann_real(

<xdnout)wannierfunc(posi,nbn,ikpt)=xdnout*factorCYCLExloopENDIFENDDOENDDO!ntENDDOENDDOENDDO!i1!Check for point in vacuumIF(film.AND..NOT.odi%d1.AND.ABS(point(3))>=z1)THENivac=1if(point(3).lt.0.0)ivac=2jvac=ivacif(nvac==1)jvac=1callwann_plot_vac(point,z1,nmzd,nv2d,n3d,nvac,>nmz,delz,bmat,bbmat,evac(:,jspin),bkpt,vz,jspin,>k1(:,jspin),k2(:,jspin),k3(:,jspin),nvd,

>ac(:,nbn,ivac),&bc(:,nbn,ivac),&u(:,:,jvac),ue(:,:,jvac),xdnout)xdnout=xdnout*factorif(real(xdnout).gt.9.0.or.real(xdnout).lt.-9.0&.or.imag(xdnout).gt.9.0.or.imag(xdnout).lt.-9.0)thenxdnout=cmplx(0.0,0.0)print*,"vac-problem at z=",point(3)endifcCALLwann_real(c>point,0,0,1,0,bkpt,c>n3d,nmzxyd,n2d,ntypsd,lmaxd,jmtd,c>natd,ntypd,nmzd,nop,nop2,mrot,tau,invtab,c>nq3,nvac,invs,z1,delz,nmz,nmzxy,nq2,c>lmax,rmsh,jri,pos,ngopr,ntypsy,nvd,c>omtil,amat,bmat,odi,ods,nlod,llod,nlo,llo,c>ff,gg,flo,acof(nbn,:,:),bcof(nbn,:,:),c>ccof(:,nbn,:,:),z(:,nbn),c>nv,k1,k2,k3,lmd,nbasfcn,c<xdnout)wannierfunc(posi,nbn,ikpt)=xdnout*factorCYCLExloopENDIFIF(odi%d1)THENIF(SQRT((pt(1))**2+(pt(2))**2)>=z1)THENCALLwann_real(