File: init1.f90

package info (click to toggle)
elkcode 5.4.24-2
  • links: PTS, VCS
  • area: main
  • in suites: buster, sid
  • size: 12,840 kB
  • sloc: f90: 48,415; fortran: 22,457; perl: 965; makefile: 384; sh: 369; python: 105; ansic: 67
file content (379 lines) | stat: -rw-r--r-- 11,142 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379

! Copyright (C) 2002-2005 J. K. Dewhurst, S. Sharma and C. Ambrosch-Draxl.
! This file is distributed under the terms of the GNU General Public License.
! See the file COPYING for license details.

!BOP
! !ROUTINE: init1
! !INTERFACE:
subroutine init1
! !USES:
use modmain
use moddftu
use modulr
use modtddft
use modtest
use modvars
use modstore
! !DESCRIPTION:
!   Generates the $k$-point set and then allocates and initialises global
!   variables which depend on the $k$-point set.
!
! !REVISION HISTORY:
!   Created January 2004 (JKD)
!EOP
!BOC
implicit none
! local variables
logical lsym(48)
integer is,ia,ias,nppt
integer io,ilo,i1,i2,i3
integer ik,igk,isym,jspn
integer l1,l2,l3,m1,m2,m3
integer lm1,lm2,lm3
real(8) vl(3),vc(3),t1
real(8) ts0,ts1
! external functions
complex(8) gauntyry
external gauntyry

call timesec(ts0)

!---------------------!
!     k-point set     !
!---------------------!
! check if the system is an isolated molecule
if (molecule) then
  ngridk(:)=1
  vkloff(:)=0.d0
  autokpt=.false.
end if
! store the point group symmetries for reducing the k-point set
if (reducek.eq.0) then
  nsymkpt=1
  symkpt(:,:,1)=symlat(:,:,1)
else
  lsym(:)=.false.
  do isym=1,nsymcrys
    if (reducek.eq.2) then
! check symmetry is symmorphic
      t1=abs(vtlsymc(1,isym))+abs(vtlsymc(2,isym))+abs(vtlsymc(3,isym))
      if (t1.gt.epslat) goto 10
! check also that the spin rotation is the same as the spatial rotation
      if (spinpol) then
        if (lspnsymc(isym).ne.lsplsymc(isym)) goto 10
      end if
    end if
    lsym(lsplsymc(isym))=.true.
10 continue
  end do
  nsymkpt=0
  do isym=1,nsymlat
    if (lsym(isym)) then
      nsymkpt=nsymkpt+1
      symkpt(:,:,nsymkpt)=symlat(:,:,isym)
    end if
  end do
end if
if ((task.eq.20).or.(task.eq.21)) then
! generate k-points along a path for band structure plots
  call plotpt1d(bvec,nvp1d,npp1d,vvlp1d,vplp1d,dvp1d,dpp1d)
  nkpt=npp1d
  if (allocated(vkl)) deallocate(vkl)
  allocate(vkl(3,nkpt))
  if (allocated(vkc)) deallocate(vkc)
  allocate(vkc(3,nkpt))
  do ik=1,nkpt
    vkl(:,ik)=vplp1d(:,ik)
    call r3mv(bvec,vkl(:,ik),vkc(:,ik))
  end do
  nkptnr=nkpt
else if (task.eq.25) then
! effective mass calculation
  nkpt=(2*ndspem+1)**3
  if (allocated(ivk)) deallocate(ivk)
  allocate(ivk(3,nkpt))
  if (allocated(vkl)) deallocate(vkl)
  allocate(vkl(3,nkpt))
  if (allocated(vkc)) deallocate(vkc)
  allocate(vkc(3,nkpt))
! map vector to [0,1)
  call r3frac(epslat,vklem)
  ik=0
  do i3=-ndspem,ndspem
    do i2=-ndspem,ndspem
      do i1=-ndspem,ndspem
        ik=ik+1
        ivk(1,ik)=i1; ivk(2,ik)=i2; ivk(3,ik)=i3
        vc(1)=dble(i1); vc(2)=dble(i2); vc(3)=dble(i3)
        vc(:)=vc(:)*deltaem
        call r3mv(binv,vc,vl)
        vkl(:,ik)=vklem(:)+vl(:)
        call r3mv(bvec,vkl(:,ik),vkc(:,ik))
      end do
    end do
  end do
  nkptnr=nkpt
else
! determine the k-point grid automatically from radkpt if required
  if (autokpt) then
    t1=radkpt/twopi
    ngridk(:)=int(t1*sqrt(bvec(1,:)**2+bvec(2,:)**2+bvec(3,:)**2))+1
  end if
! set up the default k-point box
  kptboxl(:,0)=vkloff(:)/dble(ngridk(:))
  if (task.eq.102) kptboxl(:,0)=0.d0
  kptboxl(:,1)=kptboxl(:,0)
  kptboxl(:,2)=kptboxl(:,0)
  kptboxl(:,3)=kptboxl(:,0)
  kptboxl(1,1)=kptboxl(1,1)+1.d0
  kptboxl(2,2)=kptboxl(2,2)+1.d0
  kptboxl(3,3)=kptboxl(3,3)+1.d0
! k-point set and box for Fermi surface plots
  if ((task.eq.100).or.(task.eq.101).or.(task.eq.102)) then
    ngridk(:)=np3d(:)
    if (task.ne.102) kptboxl(:,:)=vclp3d(:,:)
  end if
! allocate the k-point set arrays
  if (allocated(ivkik)) deallocate(ivkik)
  allocate(ivkik(0:ngridk(1)-1,0:ngridk(2)-1,0:ngridk(3)-1))
  if (allocated(ivkiknr)) deallocate(ivkiknr)
  allocate(ivkiknr(0:ngridk(1)-1,0:ngridk(2)-1,0:ngridk(3)-1))
  nkptnr=ngridk(1)*ngridk(2)*ngridk(3)
  if (allocated(ivk)) deallocate(ivk)
  allocate(ivk(3,nkptnr))
  if (allocated(vkl)) deallocate(vkl)
  allocate(vkl(3,nkptnr))
  if (allocated(vkc)) deallocate(vkc)
  allocate(vkc(3,nkptnr))
  if (allocated(wkpt)) deallocate(wkpt)
  allocate(wkpt(nkptnr))
! generate the k-point set
  call genppts(.false.,nsymkpt,symkpt,ngridk,nkptnr,epslat,bvec,kptboxl,nkpt, &
   ivkik,ivkiknr,ivk,vkl,vkc,wkpt,wkptnr)
! write to VARIABLES.OUT
  call writevars('nsymkpt',iv=nsymkpt)
  call writevars('symkpt',nv=9*nsymkpt,iva=symkpt)
  call writevars('ngridk',nv=3,iva=ngridk)
  call writevars('vkloff',nv=3,rva=vkloff)
  call writevars('nkpt',iv=nkpt)
  call writevars('ivkik',nv=nkptnr,iva=ivkik)
  call writevars('ivk',nv=3*nkptnr,iva=ivk)
  call writevars('vkl',nv=3*nkptnr,rva=vkl)
  call writevars('wkpt',nv=nkpt,rva=wkpt)
end if
if (task.eq.700) then
! generate ultracell reciprocal lattice vectors if required
  call reciplat(avecu,bvecu,omegau,omegabzu)
! generate the kappa, k+kappa and Q-points if required
  call genkpakq
end if
! write the k-points to test file
call writetest(910,'k-points (Cartesian)',nv=3*nkpt,tol=1.d-8,rva=vkc)

!---------------------!
!     G+k-vectors     !
!---------------------!
if ((xctype(1).lt.0).or.(task.eq.5).or.(task.eq.10).or.(task.eq.205).or. &
 (task.eq.300).or.(task.eq.600).or.(task.eq.620).or.(task.eq.630).or.tddos) then
  nppt=nkptnr
else
  nppt=nkpt
end if
! find the maximum number of G+k-vectors
call findngkmax(nkpt,vkc,nspnfv,vqcss,ngvec,vgc,gkmax,ngkmax)
! allocate the G+k-vector arrays
if (allocated(ngk)) deallocate(ngk)
allocate(ngk(nspnfv,nppt))
if (allocated(igkig)) deallocate(igkig)
allocate(igkig(ngkmax,nspnfv,nppt))
if (allocated(vgkl)) deallocate(vgkl)
allocate(vgkl(3,ngkmax,nspnfv,nppt))
if (allocated(vgkc)) deallocate(vgkc)
allocate(vgkc(3,ngkmax,nspnfv,nppt))
if (allocated(gkc)) deallocate(gkc)
allocate(gkc(ngkmax,nspnfv,nppt))
if (allocated(tpgkc)) deallocate(tpgkc)
allocate(tpgkc(2,ngkmax,nspnfv,nppt))
if (allocated(sfacgk)) deallocate(sfacgk)
allocate(sfacgk(ngkmax,natmtot,nspnfv,nppt))
do ik=1,nppt
  do jspn=1,nspnfv
    vl(:)=vkl(:,ik)
    vc(:)=vkc(:,ik)
! spin-spiral case
    if (spinsprl) then
      if (jspn.eq.1) then
        vl(:)=vl(:)+0.5d0*vqlss(:)
        vc(:)=vc(:)+0.5d0*vqcss(:)
      else
        vl(:)=vl(:)-0.5d0*vqlss(:)
        vc(:)=vc(:)-0.5d0*vqcss(:)
      end if
    end if
! generate the G+k-vectors
    call gengkvec(ngvec,ivg,vgc,vl,vc,gkmax,ngkmax,ngk(jspn,ik), &
     igkig(:,jspn,ik),vgkl(:,:,jspn,ik),vgkc(:,:,jspn,ik))
! generate the spherical coordinates of the G+k-vectors
    do igk=1,ngk(jspn,ik)
      call sphcrd(vgkc(:,igk,jspn,ik),gkc(igk,jspn,ik),tpgkc(:,igk,jspn,ik))
    end do
! generate structure factors for G+k-vectors
    call gensfacgp(ngk(jspn,ik),vgkc(:,:,jspn,ik),ngkmax,sfacgk(:,:,jspn,ik))
  end do
end do
! write to VARIABLES.OUT
call writevars('nspnfv',iv=nspnfv)
call writevars('ngk',nv=nspnfv*nkpt,iva=ngk)
do ik=1,nkpt
  do jspn=1,nspnfv
    call writevars('igkig',l=jspn,m=ik,nv=ngk(jspn,ik),iva=igkig(:,jspn,ik))
  end do
end do

!---------------------------------!
!     APWs and local-orbitals     !
!---------------------------------!
apwordmax=0
lorbordmax=0
nlomax=0
lolmax=0
do is=1,nspecies
! find the maximum APW order
  do l1=0,lmaxapw
    apwordmax=max(apwordmax,apword(l1,is))
  end do
! find the maximum number of local-orbitals
  nlomax=max(nlomax,nlorb(is))
! find the maximum local-orbital order and angular momentum
  do ilo=1,nlorb(is)
    lolmax=max(lolmax,lorbl(ilo,is))
    lorbordmax=max(lorbordmax,lorbord(ilo,is))
  end do
end do
lolmmax=(lolmax+1)**2
! polynomial order used for APW and local-orbital radial derivatives
npapw=max(apwordmax+1,4)
nplorb=max(lorbordmax+1,4)
! set the APW and local-orbital linearisation energies to the default
if (allocated(apwe)) deallocate(apwe)
allocate(apwe(apwordmax,0:lmaxapw,natmtot))
if (allocated(lorbe)) deallocate(lorbe)
allocate(lorbe(lorbordmax,maxlorb,natmtot))
do is=1,nspecies
  do l1=0,lmaxapw
    do io=1,apword(l1,is)
      do ia=1,natoms(is)
        ias=idxas(ia,is)
        apwe(io,l1,ias)=apwe0(io,l1,is)
      end do
    end do
  end do
  do ilo=1,nlorb(is)
    do io=1,lorbord(ilo,is)
      do ia=1,natoms(is)
        ias=idxas(ia,is)
        lorbe(io,ilo,ias)=lorbe0(io,ilo,is)
      end do
    end do
  end do
end do
! generate the local-orbital index
call genidxlo
! allocate radial function arrays
if (allocated(apwfr)) deallocate(apwfr)
allocate(apwfr(nrmtmax,2,apwordmax,0:lmaxapw,natmtot))
if (allocated(apwdfr)) deallocate(apwdfr)
allocate(apwdfr(apwordmax,0:lmaxapw,natmtot))
if (allocated(lofr)) deallocate(lofr)
allocate(lofr(nrmtmax,2,nlomax,natmtot))

!-------------------------!
!     DFT+U variables     !
!-------------------------!
if (dftu.ne.0) then
! allocate energy arrays to calculate Slater integrals with Yukawa potential
  if (allocated(fdue)) deallocate(fdue)
  allocate(fdue(0:lmaxdm,natmtot))
! allocate radial functions to calculate Slater integrals with Yukawa potential
  if (allocated(fdufr)) deallocate(fdufr)
  allocate(fdufr(nrmtmax,0:lmaxdm,natmtot))
end if

!---------------------------------------!
!     eigenvalue equation variables     !
!---------------------------------------!
! total number of empty states (M. Meinert)
nempty=nint(nempty0*max(natmtot,1))
if (nempty.lt.1) nempty=1
! number of first-variational states
nstfv=int(chgval/2.d0)+nempty+1
! overlap and Hamiltonian matrix sizes
if (allocated(nmat)) deallocate(nmat)
allocate(nmat(nspnfv,nkpt))
nmatmax=0
do ik=1,nkpt
  do jspn=1,nspnfv
    nmat(jspn,ik)=ngk(jspn,ik)+nlotot
    if (nstfv.gt.nmat(jspn,ik)) then
      write(*,*)
      write(*,'("Error(init1): number of first-variational states larger than &
       &matrix size")')
      write(*,'("Increase rgkmax or decrease nempty")')
      write(*,*)
      stop
    end if
    nmatmax=max(nmatmax,nmat(jspn,ik))
  end do
end do
! number of second-variational states
nstsv=nstfv*nspinor
! allocate second-variational arrays
if (allocated(evalsv)) deallocate(evalsv)
allocate(evalsv(nstsv,nkpt))
if (allocated(occsv)) deallocate(occsv)
allocate(occsv(nstsv,nkpt))
occsv(:,:)=0.d0
! allocate overlap and Hamiltonian integral arrays
if (allocated(oalo)) deallocate(oalo)
allocate(oalo(apwordmax,nlomax,natmtot))
if (allocated(ololo)) deallocate(ololo)
allocate(ololo(nlomax,nlomax,natmtot))
if (allocated(haa)) deallocate(haa)
allocate(haa(lmmaxo,apwordmax,0:lmaxapw,apwordmax,0:lmaxapw,natmtot))
if (allocated(hloa)) deallocate(hloa)
allocate(hloa(lmmaxo,apwordmax,0:lmaxapw,nlomax,natmtot))
if (allocated(hlolo)) deallocate(hlolo)
allocate(hlolo(lmmaxo,nlomax,nlomax,natmtot))
! allocate and generate complex Gaunt coefficient array
if (allocated(gntyry)) deallocate(gntyry)
allocate(gntyry(lmmaxapw,lmmaxo,lmmaxapw))
do l1=0,lmaxapw
  do m1=-l1,l1
    lm1=idxlm(l1,m1)
    do l2=0,lmaxo
      do m2=-l2,l2
        lm2=idxlm(l2,m2)
        do l3=0,lmaxapw
          do m3=-l3,l3
            lm3=idxlm(l3,m3)
            gntyry(lm1,lm2,lm3)=gauntyry(l1,l2,l3,m1,m2,m3)
          end do
        end do
      end do
    end do
  end do
end do
! write to VARIABLES.OUT
call writevars('nempty',iv=nempty)
call writevars('nstfv',iv=nstfv)
call writevars('nlotot',iv=nlotot)
call writevars('nstsv',iv=nstsv)

call timesec(ts1)
timeinit=timeinit+ts1-ts0

return
end subroutine
!EOC