src/init1.f90

   1
   2 ! Copyright (C) 2002-2005 J. K. Dewhurst, S. Sharma and C. Ambrosch-Draxl.
   3 ! This file is distributed under the terms of the GNU General Public License.
   4 ! See the file COPYING for license details.
   5
   6 !BOP
   7 ! !ROUTINE: init1
   8 ! !INTERFACE:
   9 subroutine init1
  10 ! !USES:
  11 use modmain
  12 ! !DESCRIPTION:
  13 !   Generates the $k$-point set and then allocates and initialises global
  14 !   variables which depend on the $k$-point set.
  15 !
  16 ! !REVISION HISTORY:
  17 !   Created January 2004 (JKD)
  18 !EOP
  19 !BOC
  20 implicit none
  21 ! local variables
  22 integer ik,is,ia,ias,io,ilo
  23 integer i1,i2,i3,ispn,iv(3)
  24 integer l1,l2,l3,m1,m2,m3,lm1,lm2,lm3
  25 real(8) vl(3),vc(3),boxl(3,4)
  26 real(8) ts0,ts1
  27 ! external functions
  28 complex(8) gauntyry
  29 external gauntyry
  30
  31 call timesec(ts0)
  32
  33 !---------------------!
  34 !     k-point set     !
  35 !---------------------!
  36 ! check if the system is an isolated molecule
  37 if (molecule) then
  38   ngridk(:)=1
  39   vkloff(:)=0.d0
  40   autokpt=.false.
  41 end if
  42 ! setup the default k-point box
  43 boxl(:,:)=0.d0
  44 boxl(:,1)=vkloff(:)/dble(ngridk(:))
  45 boxl(1,2)=1.d0; boxl(2,3)=1.d0; boxl(3,4)=1.d0
  46 ! k-point set and box for Fermi surface plots
  47 if ((task.eq.100).or.(task.eq.101)) then
  48   ngridk(:)=np3d(:)
  49   boxl(:,:)=vclp3d(:,:)
  50 end if
  51 if ((task.eq.20).or.(task.eq.21)) then
  52 ! for band structure plots generate k-points along a line
  53   call connect(bvec,nvp1d,npp1d,vvlp1d,vplp1d,dvp1d,dpp1d)
  54   nkpt=npp1d
  55   if (allocated(vkl)) deallocate(vkl)
  56   allocate(vkl(3,nkpt))
  57   if (allocated(vkc)) deallocate(vkc)
  58   allocate(vkc(3,nkpt))
  59   do ik=1,nkpt
  60     vkl(:,ik)=vplp1d(:,ik)
  61     call r3mv(bvec,vkl(:,ik),vkc(:,ik))
  62   end do
  63 else if (task.eq.25) then
  64 ! effective mass calculation
  65   nkpt=(2*ndspem+1)**3
  66   if (allocated(ivk)) deallocate(ivk)
  67   allocate(ivk(3,nkpt))
  68   if (allocated(vkl)) deallocate(vkl)
  69   allocate(vkl(3,nkpt))
  70   if (allocated(vkc)) deallocate(vkc)
  71   allocate(vkc(3,nkpt))
  72 ! map vector to [0,1)
  73   call r3frac(epslat,vklem,iv)
  74   ik=0
  75   do i3=-ndspem,ndspem
  76     do i2=-ndspem,ndspem
  77       do i1=-ndspem,ndspem
  78         ik=ik+1
  79         ivk(1,ik)=i1; ivk(2,ik)=i2; ivk(3,ik)=i3
  80         vc(1)=dble(i1); vc(2)=dble(i2); vc(3)=dble(i3)
  81         vc(:)=vc(:)*deltaem
  82         call r3mv(binv,vc,vl)
  83         vkl(:,ik)=vklem(:)+vl(:)
  84         call r3mv(bvec,vkl(:,ik),vkc(:,ik))
  85       end do
  86     end do
  87   end do
  88 else
  89 ! determine the k-point grid automatically from radkpt if required
  90   if (autokpt) then
  91     ngridk(:)=int(radkpt/sqrt(avec(1,:)**2+avec(2,:)**2+avec(3,:)**2))+1
  92   end if
  93 ! allocate the reduced k-point set arrays
  94   if (allocated(ivk)) deallocate(ivk)
  95   allocate(ivk(3,ngridk(1)*ngridk(2)*ngridk(3)))
  96   if (allocated(vkl)) deallocate(vkl)
  97   allocate(vkl(3,ngridk(1)*ngridk(2)*ngridk(3)))
  98   if (allocated(vkc)) deallocate(vkc)
  99   allocate(vkc(3,ngridk(1)*ngridk(2)*ngridk(3)))
 100   if (allocated(wkpt)) deallocate(wkpt)
 101   allocate(wkpt(ngridk(1)*ngridk(2)*ngridk(3)))
 102   if (allocated(ikmap)) deallocate(ikmap)
 103   allocate(ikmap(0:ngridk(1)-1,0:ngridk(2)-1,0:ngridk(3)-1))
 104 ! generate the reduced k-point set
 105   call genppts(reducek,.false.,ngridk,boxl,nkpt,ikmap,ivk,vkl,vkc,wkpt)
 106 ! allocate the non-reduced k-point set arrays
 107   nkptnr=ngridk(1)*ngridk(2)*ngridk(3)
 108   if (allocated(ivknr)) deallocate(ivknr)
 109   allocate(ivknr(3,nkptnr))
 110   if (allocated(vklnr)) deallocate(vklnr)
 111   allocate(vklnr(3,nkptnr))
 112   if (allocated(vkcnr)) deallocate(vkcnr)
 113   allocate(vkcnr(3,nkptnr))
 114   if (allocated(wkptnr)) deallocate(wkptnr)
 115   allocate(wkptnr(nkptnr))
 116   if (allocated(ikmapnr)) deallocate(ikmapnr)
 117   allocate(ikmapnr(0:ngridk(1)-1,0:ngridk(2)-1,0:ngridk(3)-1))
 118 ! generate the non-reduced k-point set
 119   call genppts(.false.,.false.,ngridk,boxl,nkptnr,ikmapnr,ivknr,vklnr,vkcnr, &
 120    wkptnr)
 121 end if
 122
 123 !---------------------!
 124 !     G+k vectors     !
 125 !---------------------!
 126 ! determine gkmax
 127 if ((isgkmax.ge.1).and.(isgkmax.le.nspecies)) then
 128   gkmax=rgkmax/rmt(isgkmax)
 129 else
 130   gkmax=rgkmax/2.d0
 131 end if
 132 if (2.d0*gkmax.gt.gmaxvr+epslat) then
 133   write(*,*)
 134   write(*,'("Error(init1): 2*gkmax > gmaxvr  ",2G18.10)') 2.d0*gkmax,gmaxvr
 135   write(*,*)
 136   stop
 137 end if
 138 ! find the maximum number of G+k-vectors
 139 call getngkmax
 140 ! allocate the G+k-vector arrays
 141 if (allocated(ngk)) deallocate(ngk)
 142 allocate(ngk(nspnfv,nkpt))
 143 if (allocated(igkig)) deallocate(igkig)
 144 allocate(igkig(ngkmax,nspnfv,nkpt))
 145 if (allocated(vgkl)) deallocate(vgkl)
 146 allocate(vgkl(3,ngkmax,nspnfv,nkpt))
 147 if (allocated(vgkc)) deallocate(vgkc)
 148 allocate(vgkc(3,ngkmax,nspnfv,nkpt))
 149 if (allocated(gkc)) deallocate(gkc)
 150 allocate(gkc(ngkmax,nspnfv,nkpt))
 151 if (allocated(tpgkc)) deallocate(tpgkc)
 152 allocate(tpgkc(2,ngkmax,nspnfv,nkpt))
 153 if (allocated(sfacgk)) deallocate(sfacgk)
 154 allocate(sfacgk(ngkmax,natmtot,nspnfv,nkpt))
 155 do ik=1,nkpt
 156   do ispn=1,nspnfv
 157     if (spinsprl) then
 158 ! spin-spiral case
 159       if (ispn.eq.1) then
 160         vl(:)=vkl(:,ik)+0.5d0*vqlss(:)
 161         vc(:)=vkc(:,ik)+0.5d0*vqcss(:)
 162       else
 163         vl(:)=vkl(:,ik)-0.5d0*vqlss(:)
 164         vc(:)=vkc(:,ik)-0.5d0*vqcss(:)
 165       end if
 166     else
 167       vl(:)=vkl(:,ik)
 168       vc(:)=vkc(:,ik)
 169     end if
 170 ! generate the G+k-vectors
 171     call gengpvec(vl,vc,ngk(ispn,ik),igkig(:,ispn,ik),vgkl(:,:,ispn,ik), &
 172      vgkc(:,:,ispn,ik),gkc(:,ispn,ik),tpgkc(:,:,ispn,ik))
 173 ! generate structure factors for G+k-vectors
 174     call gensfacgp(ngk(ispn,ik),vgkc(:,:,ispn,ik),ngkmax,sfacgk(:,:,ispn,ik))
 175   end do
 176 end do
 177
 178 !---------------------------------!
 179 !     APWs and local-orbitals     !
 180 !---------------------------------!
 181 ! allocate linearisation energy arrays
 182 if (allocated(apwe)) deallocate(apwe)
 183 allocate(apwe(maxapword,0:lmaxapw,natmtot))
 184 if (allocated(lorbe)) deallocate(lorbe)
 185 allocate(lorbe(maxlorbord,maxlorb,natmtot))
 186 nlomax=0
 187 lolmax=0
 188 apwordmax=0
 189 do is=1,nspecies
 190 ! find the maximum APW order
 191   do l1=0,lmaxapw
 192     apwordmax=max(apwordmax,apword(l1,is))
 193   end do
 194 ! set the APW linearisation energies to the default
 195   do ia=1,natoms(is)
 196     ias=idxas(ia,is)
 197     do l1=0,lmaxapw
 198       do io=1,apword(l1,is)
 199         apwe(io,l1,ias)=apwe0(io,l1,is)
 200       end do
 201     end do
 202   end do
 203 ! find the maximum number of local-orbitals
 204   nlomax=max(nlomax,nlorb(is))
 205 ! set the local-orbital linearisation energies to the default
 206   do ia=1,natoms(is)
 207     ias=idxas(ia,is)
 208     do ilo=1,nlorb(is)
 209       lolmax=max(lolmax,lorbl(ilo,is))
 210       do io=1,lorbord(ilo,is)
 211         lorbe(io,ilo,ias)=lorbe0(io,ilo,is)
 212       end do
 213     end do
 214   end do
 215 end do
 216 lolmmax=(lolmax+1)**2
 217 ! generate the local-orbital index
 218 call genidxlo
 219 ! allocate radial function arrays
 220 if (allocated(apwfr)) deallocate(apwfr)
 221 allocate(apwfr(nrmtmax,2,apwordmax,0:lmaxapw,natmtot))
 222 if (allocated(apwdfr)) deallocate(apwdfr)
 223 allocate(apwdfr(apwordmax,0:lmaxapw,natmtot))
 224 if (allocated(lofr)) deallocate(lofr)
 225 allocate(lofr(nrmtmax,2,nlomax,natmtot))
 226
 227 !------------------------------------!
 228 !     secular equation variables     !
 229 !------------------------------------!
 230 ! number of first-variational states
 231 nstfv=int(chgval/2.d0)+nempty+1
 232 ! overlap and Hamiltonian matrix sizes
 233 if (allocated(nmat)) deallocate(nmat)
 234 allocate(nmat(nspnfv,nkpt))
 235 if (allocated(npmat)) deallocate(npmat)
 236 allocate(npmat(nspnfv,nkpt))
 237 nmatmax=0
 238 do ik=1,nkpt
 239   do ispn=1,nspnfv
 240     nmat(ispn,ik)=ngk(ispn,ik)+nlotot
 241     nmatmax=max(nmatmax,nmat(ispn,ik))
 242 ! packed matrix sizes
 243     npmat(ispn,ik)=(nmat(ispn,ik)*(nmat(ispn,ik)+1))/2
 244 ! the number of first-variational states should not exceed the matrix size
 245     nstfv=min(nstfv,nmat(ispn,ik))
 246   end do
 247 end do
 248 ! number of second-variational states
 249 nstsv=nstfv*nspinor
 250 ! allocate second-variational arrays
 251 if (allocated(evalsv)) deallocate(evalsv)
 252 allocate(evalsv(nstsv,nkpt))
 253 if (allocated(occsv)) deallocate(occsv)
 254 allocate(occsv(nstsv,nkpt))
 255 occsv(:,:)=0.d0
 256 ! allocate overlap and Hamiltonian integral arrays
 257 if (allocated(oalo)) deallocate(oalo)
 258 allocate(oalo(apwordmax,nlomax,natmtot))
 259 if (allocated(ololo)) deallocate(ololo)
 260 allocate(ololo(nlomax,nlomax,natmtot))
 261 if (allocated(haa)) deallocate(haa)
 262 allocate(haa(apwordmax,0:lmaxmat,apwordmax,0:lmaxapw,lmmaxvr,natmtot))
 263 if (allocated(hloa)) deallocate(hloa)
 264 allocate(hloa(nlomax,apwordmax,0:lmaxmat,lmmaxvr,natmtot))
 265 if (allocated(hlolo)) deallocate(hlolo)
 266 allocate(hlolo(nlomax,nlomax,lmmaxvr,natmtot))
 267 ! allocate and generate complex Gaunt coefficient array
 268 if (allocated(gntyry)) deallocate(gntyry)
 269 allocate(gntyry(lmmaxmat,lmmaxvr,lmmaxapw))
 270 do l1=0,lmaxmat
 271   do m1=-l1,l1
 272     lm1=idxlm(l1,m1)
 273     do l2=0,lmaxvr
 274       do m2=-l2,l2
 275         lm2=idxlm(l2,m2)
 276         do l3=0,lmaxapw
 277           do m3=-l3,l3
 278             lm3=idxlm(l3,m3)
 279             gntyry(lm1,lm2,lm3)=gauntyry(l1,l2,l3,m1,m2,m3)
 280           end do
 281         end do
 282       end do
 283     end do
 284   end do
 285 end do
 286
 287 call timesec(ts1)
 288 timeinit=timeinit+ts1-ts0
 289
 290 return
 291 end subroutine
 292 !EOC
 293