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: 9subroutine init1 10! !USES: 11use modmain 12use moddftu 13use modulr 14use modtddft 15use modtest 16use modvars 17use modstore 18! !DESCRIPTION: 19! Generates the $k$-point set and then allocates and initialises global 20! variables which depend on the $k$-point set. 21! 22! !REVISION HISTORY: 23! Created January 2004 (JKD) 24!EOP 25!BOC 26implicit none 27! local variables 28logical lsym(48) 29integer is,ia,ias,nppt 30integer io,ilo,i1,i2,i3 31integer ik,isym,jspn 32integer l1,l2,l3,m1,m2,m3 33integer lm1,lm2,lm3,n 34real(8) vl(3),vc(3),t1 35real(8) boxl(3,0:3) 36real(8) ts0,ts1 37! external functions 38complex(8), external :: gauntyry 39 40call timesec(ts0) 41 42!---------------------! 43! k-point set ! 44!---------------------! 45! check if the system is an isolated molecule 46if (molecule) then 47 ngridk(:)=1 48 vkloff(:)=0.d0 49 autokpt=.false. 50end if 51! store the point group symmetries for reducing the k-point set 52if (reducek.eq.0) then 53 nsymkpt=1 54 symkpt(:,:,1)=symlat(:,:,1) 55else 56 lsym(:)=.false. 57 do isym=1,nsymcrys 58 if (reducek.eq.2) then 59! check symmetry is symmorphic 60 if (.not.tv0symc(isym)) goto 10 61! check also that the spin rotation is the same as the spatial rotation 62 if (spinpol) then 63 if (lspnsymc(isym).ne.lsplsymc(isym)) goto 10 64 end if 65 end if 66 lsym(lsplsymc(isym))=.true. 6710 continue 68 end do 69 nsymkpt=0 70 do isym=1,nsymlat 71 if (lsym(isym)) then 72 nsymkpt=nsymkpt+1 73 symkpt(:,:,nsymkpt)=symlat(:,:,isym) 74 end if 75 end do 76end if 77if (any(task.eq.[20,21,22,23])) then 78! generate k-points along a path for band structure plots 79 call plotpt1d(bvec,nvp1d,npp1d,vvlp1d,vplp1d,dvp1d,dpp1d) 80 nkpt=npp1d 81 if (allocated(vkl)) deallocate(vkl) 82 allocate(vkl(3,nkpt)) 83 if (allocated(vkc)) deallocate(vkc) 84 allocate(vkc(3,nkpt)) 85 do ik=1,nkpt 86 vkl(:,ik)=vplp1d(:,ik) 87 call r3mv(bvec,vkl(:,ik),vkc(:,ik)) 88 end do 89 nkptnr=nkpt 90else if (task.eq.25) then 91! effective mass calculation 92 nkpt=(2*ndspem+1)**3 93 if (allocated(ivk)) deallocate(ivk) 94 allocate(ivk(3,nkpt)) 95 if (allocated(vkl)) deallocate(vkl) 96 allocate(vkl(3,nkpt)) 97 if (allocated(vkc)) deallocate(vkc) 98 allocate(vkc(3,nkpt)) 99! map vector to [0,1) 100 call r3frac(epslat,vklem) 101 ik=0 102 do i3=-ndspem,ndspem 103 do i2=-ndspem,ndspem 104 do i1=-ndspem,ndspem 105 ik=ik+1 106 ivk(1,ik)=i1; ivk(2,ik)=i2; ivk(3,ik)=i3 107 vc(1)=dble(i1); vc(2)=dble(i2); vc(3)=dble(i3) 108 vc(:)=vc(:)*deltaem 109 call r3mv(binv,vc,vl) 110 vkl(:,ik)=vklem(:)+vl(:) 111 call r3mv(bvec,vkl(:,ik),vkc(:,ik)) 112 end do 113 end do 114 end do 115 nkptnr=nkpt 116else 117! determine the k-point grid automatically from radkpt if required 118 if (autokpt) then 119 t1=radkpt/twopi 120 ngridk(:)=int(t1*sqrt(bvec(1,:)**2+bvec(2,:)**2+bvec(3,:)**2))+1 121 end if 122! set up the default k-point box 123 boxl(:,0)=vkloff(:)/dble(ngridk(:)) 124 if (task.eq.102) boxl(:,0)=0.d0 125 boxl(:,1)=boxl(:,0) 126 boxl(:,2)=boxl(:,0) 127 boxl(:,3)=boxl(:,0) 128 boxl(1,1)=boxl(1,1)+1.d0 129 boxl(2,2)=boxl(2,2)+1.d0 130 boxl(3,3)=boxl(3,3)+1.d0 131! k-point set and box for Fermi surface plots 132 if (any(task.eq.[100,101,102])) then 133 ngridk(:)=np3d(:) 134 if (task.ne.102) boxl(:,:)=vclp3d(:,:) 135 end if 136! allocate the k-point set arrays 137 if (allocated(ivkik)) deallocate(ivkik) 138 allocate(ivkik(0:ngridk(1)-1,0:ngridk(2)-1,0:ngridk(3)-1)) 139 if (allocated(ivkiknr)) deallocate(ivkiknr) 140 allocate(ivkiknr(0:ngridk(1)-1,0:ngridk(2)-1,0:ngridk(3)-1)) 141 nkptnr=ngridk(1)*ngridk(2)*ngridk(3) 142 if (allocated(ivk)) deallocate(ivk) 143 allocate(ivk(3,nkptnr)) 144 if (allocated(vkl)) deallocate(vkl) 145 allocate(vkl(3,nkptnr)) 146 if (allocated(vkc)) deallocate(vkc) 147 allocate(vkc(3,nkptnr)) 148 if (allocated(wkpt)) deallocate(wkpt) 149 allocate(wkpt(nkptnr)) 150! generate the k-point set 151 call genppts(.false.,nsymkpt,symkpt,ngridk,nkptnr,epslat,bvec,boxl,nkpt, & 152 ivkik,ivkiknr,ivk,vkl,vkc,wkpt,wkptnr) 153! write to VARIABLES.OUT 154 call writevars('nsymkpt',iv=nsymkpt) 155 call writevars('symkpt',nv=9*nsymkpt,iva=symkpt) 156 call writevars('ngridk',nv=3,iva=ngridk) 157 call writevars('vkloff',nv=3,rva=vkloff) 158 call writevars('nkpt',iv=nkpt) 159 call writevars('ivkik',nv=nkptnr,iva=ivkik) 160 call writevars('ivk',nv=3*nkptnr,iva=ivk) 161 call writevars('vkl',nv=3*nkptnr,rva=vkl) 162 call writevars('wkpt',nv=nkpt,rva=wkpt) 163end if 164if (any(task.eq.[700,701,731,732,733,741,742,743,771,772,773])) then 165! generate ultracell reciprocal lattice vectors if required 166 call reciplat(avecu,bvecu,omegau,omegabzu) 167! generate the kappa, k+kappa and Q-points if required 168 call genkpakq 169end if 170! write the k-points to test file 171call writetest(910,'k-points (Cartesian)',nv=3*nkpt,tol=1.d-8,rva=vkc) 172 173!---------------------! 174! G+k-vectors ! 175!---------------------! 176if ((xctype(1).lt.0).or.tddos.or. & 177 (any(task.eq.[5,10,205,300,600,620,630,800,801]))) then 178 nppt=nkptnr 179else 180 nppt=nkpt 181end if 182! find the maximum number of G+k-vectors 183call findngkmax(nkpt,vkc,nspnfv,vqcss,ngvc,vgc,gkmax,ngkmax) 184! allocate the G+k-vector arrays 185if (allocated(ngk)) deallocate(ngk) 186allocate(ngk(nspnfv,nppt)) 187if (allocated(igkig)) deallocate(igkig) 188allocate(igkig(ngkmax,nspnfv,nppt)) 189if (allocated(vgkl)) deallocate(vgkl) 190allocate(vgkl(3,ngkmax,nspnfv,nppt)) 191if (allocated(vgkc)) deallocate(vgkc) 192allocate(vgkc(3,ngkmax,nspnfv,nppt)) 193if (allocated(gkc)) deallocate(gkc) 194allocate(gkc(ngkmax,nspnfv,nppt)) 195if (allocated(sfacgk)) deallocate(sfacgk) 196allocate(sfacgk(ngkmax,natmtot,nspnfv,nppt)) 197do ik=1,nppt 198 do jspn=1,nspnfv 199 vl(:)=vkl(:,ik) 200 vc(:)=vkc(:,ik) 201! spin-spiral case 202 if (spinsprl) then 203 if (jspn.eq.1) then 204 vl(:)=vl(:)+0.5d0*vqlss(:) 205 vc(:)=vc(:)+0.5d0*vqcss(:) 206 else 207 vl(:)=vl(:)-0.5d0*vqlss(:) 208 vc(:)=vc(:)-0.5d0*vqcss(:) 209 end if 210 end if 211! generate the G+k-vectors 212 call gengkvec(ngvc,ivg,vgc,vl,vc,gkmax,ngkmax,ngk(jspn,ik), & 213 igkig(:,jspn,ik),vgkl(:,:,jspn,ik),vgkc(:,:,jspn,ik),gkc(:,jspn,ik)) 214! generate structure factors for G+k-vectors 215 call gensfacgp(ngk(jspn,ik),vgkc(:,:,jspn,ik),ngkmax,sfacgk(:,:,jspn,ik)) 216 end do 217end do 218! write to VARIABLES.OUT 219call writevars('nspnfv',iv=nspnfv) 220call writevars('ngk',nv=nspnfv*nkpt,iva=ngk) 221do ik=1,nkpt 222 do jspn=1,nspnfv 223 call writevars('igkig',l=jspn,m=ik,nv=ngk(jspn,ik),iva=igkig(:,jspn,ik)) 224 end do 225end do 226 227!---------------------------------! 228! APWs and local-orbitals ! 229!---------------------------------! 230apwordmax=0 231lorbordmax=0 232nlomax=0 233lolmax=0 234do is=1,nspecies 235 lmoapw(is)=0 236 do l1=0,lmaxapw 237! find the maximum APW order 238 apwordmax=max(apwordmax,apword(l1,is)) 239! find total number of APW coefficients (l, m and order) 240 lmoapw(is)=lmoapw(is)+(2*l1+1)*apword(l1,is) 241 end do 242! find the maximum number of local-orbitals 243 nlomax=max(nlomax,nlorb(is)) 244! find the maximum local-orbital order and angular momentum 245 do ilo=1,nlorb(is) 246 lolmax=max(lolmax,lorbl(ilo,is)) 247 lorbordmax=max(lorbordmax,lorbord(ilo,is)) 248 end do 249end do 250lolmmax=(lolmax+1)**2 251! polynomial order used for APW and local-orbital radial derivatives 252npapw=max(apwordmax+1,4) 253nplorb=max(lorbordmax+1,4) 254! set the APW and local-orbital linearisation energies to the default 255if (allocated(apwe)) deallocate(apwe) 256allocate(apwe(apwordmax,0:lmaxapw,natmtot)) 257if (allocated(lorbe)) deallocate(lorbe) 258allocate(lorbe(lorbordmax,maxlorb,natmtot)) 259do is=1,nspecies 260 do l1=0,lmaxapw 261 do io=1,apword(l1,is) 262 do ia=1,natoms(is) 263 ias=idxas(ia,is) 264 apwe(io,l1,ias)=apwe0(io,l1,is) 265 end do 266 end do 267 end do 268 do ilo=1,nlorb(is) 269 do io=1,lorbord(ilo,is) 270 do ia=1,natoms(is) 271 ias=idxas(ia,is) 272 lorbe(io,ilo,ias)=lorbe0(io,ilo,is) 273 end do 274 end do 275 end do 276end do 277! generate the local-orbital index 278call genidxlo 279! allocate radial function arrays 280if (allocated(apwfr)) deallocate(apwfr) 281allocate(apwfr(nrmtmax,2,apwordmax,0:lmaxapw,natmtot)) 282if (allocated(apwdfr)) deallocate(apwdfr) 283allocate(apwdfr(apwordmax,0:lmaxapw,natmtot)) 284if (allocated(lofr)) deallocate(lofr) 285allocate(lofr(nrmtmax,2,nlomax,natmtot)) 286 287!-------------------------! 288! DFT+U variables ! 289!-------------------------! 290if (dftu.ne.0) then 291! allocate energy arrays to calculate Slater integrals with Yukawa potential 292 if (allocated(fdue)) deallocate(fdue) 293 allocate(fdue(0:lmaxdm,natmtot)) 294! allocate radial functions to calculate Slater integrals with Yukawa potential 295 if (allocated(fdufr)) deallocate(fdufr) 296 allocate(fdufr(nrmtmax,0:lmaxdm,natmtot)) 297end if 298 299!---------------------------------------! 300! eigenvalue equation variables ! 301!---------------------------------------! 302! total number of empty states (M. Meinert) 303nempty=nint(nempty0*max(natmtot,1)) 304if (nempty.lt.1) nempty=1 305! number of first-variational states 306nstfv=int(chgval/2.d0)+nempty+1 307! overlap and Hamiltonian matrix sizes 308if (allocated(nmat)) deallocate(nmat) 309allocate(nmat(nspnfv,nkpt)) 310nmatmax=0 311do ik=1,nkpt 312 do jspn=1,nspnfv 313 n=ngk(jspn,ik)+nlotot 314 if (nstfv.gt.n) then 315 write(*,*) 316 write(*,'("Error(init1): number of first-variational states larger than & 317 &matrix size")') 318 write(*,'("Increase rgkmax or decrease nempty")') 319 write(*,*) 320 stop 321 end if 322 nmat(jspn,ik)=n 323 nmatmax=max(nmatmax,n) 324 end do 325end do 326! number of second-variational states 327nstsv=nstfv*nspinor 328! allocate second-variational arrays 329if (allocated(evalsv)) deallocate(evalsv) 330allocate(evalsv(nstsv,nkpt)) 331if (allocated(occsv)) deallocate(occsv) 332allocate(occsv(nstsv,nkpt)) 333occsv(:,:)=0.d0 334! allocate overlap and Hamiltonian integral arrays 335if (allocated(oalo)) deallocate(oalo) 336allocate(oalo(apwordmax,nlomax,natmtot)) 337if (allocated(ololo)) deallocate(ololo) 338allocate(ololo(nlomax,nlomax,natmtot)) 339if (allocated(haa)) deallocate(haa) 340allocate(haa(lmmaxo,apwordmax,0:lmaxapw,apwordmax,0:lmaxapw,natmtot)) 341if (allocated(hloa)) deallocate(hloa) 342allocate(hloa(lmmaxo,apwordmax,0:lmaxapw,nlomax,natmtot)) 343if (allocated(hlolo)) deallocate(hlolo) 344allocate(hlolo(lmmaxo,nlomax,nlomax,natmtot)) 345! allocate and generate complex Gaunt coefficient array 346if (allocated(gntyry)) deallocate(gntyry) 347allocate(gntyry(lmmaxo,lmmaxapw,lmmaxapw)) 348do l1=0,lmaxapw 349 do m1=-l1,l1 350 lm1=l1*(l1+1)+m1+1 351 do l3=0,lmaxapw 352 do m3=-l3,l3 353 lm3=l3*(l3+1)+m3+1 354 do l2=0,lmaxo 355 do m2=-l2,l2 356 lm2=l2*(l2+1)+m2+1 357 gntyry(lm2,lm3,lm1)=gauntyry(l1,l2,l3,m1,m2,m3) 358 end do 359 end do 360 end do 361 end do 362 end do 363end do 364! write to VARIABLES.OUT 365call writevars('nempty',iv=nempty) 366call writevars('nstfv',iv=nstfv) 367call writevars('nlotot',iv=nlotot) 368call writevars('nstsv',iv=nstsv) 369 370call timesec(ts1) 371timeinit=timeinit+ts1-ts0 372 373end subroutine 374!EOC 375 376