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