File: gndstate.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 (427 lines) | stat: -rw-r--r-- 12,920 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
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427

! Copyright (C) 2002-2013 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: gndstate
! !INTERFACE:
subroutine gndstate
! !USES:
use modmain
use moddftu
use modulr
use modmpi
use modomp
! !DESCRIPTION:
!   Computes the self-consistent Kohn-Sham ground-state. General information is
!   written to the file {\tt INFO.OUT}. First- and second-variational
!   eigenvalues, eigenvectors and occupancies are written to the unformatted
!   files {\tt EVALFV.OUT}, {\tt EVALSV.OUT}, {\tt EVECFV.OUT}, {\tt EVECSV.OUT}
!   and {\tt OCCSV.OUT}. The density, magnetisation, Kohn-Sham potential and
!   magnetic field are written to {\tt STATE.OUT}.
!
! !REVISION HISTORY:
!   Created October 2002 (JKD)
!   Added MPI, August 2010 (JKD)
!EOP
!BOC
implicit none
! local variables
logical exist
integer ik,nwork,n
real(8) dv,etp,de,timetot
! allocatable arrays
real(8), allocatable :: v(:),work(:)
! initialise global variables
call init0
call init1
! initialise q-vector-dependent variables if required
if (xctype(1).lt.0) call init2
! apply strain to the G, k, G+k and q-vectors if required
call straingkq
if (task.eq.0) trdstate=.false.
if (task.eq.1) trdstate=.true.
! only the MPI master process should write files
if (mp_mpi) then
! write the real and reciprocal lattice vectors to file
  call writelat
! write symmetry matrices to file
  call writesym
! output the k-point set to file
  call writekpts
! write lattice vectors and atomic positions to file
  open(50,file='GEOMETRY'//trim(filext),form='FORMATTED')
  call writegeom(50)
  close(50)
! write interatomic distances to file
  open(50,file='IADIST'//trim(filext),form='FORMATTED')
  call writeiad(50)
  close(50)
! open INFO.OUT file
  open(60,file='INFO'//trim(filext),form='FORMATTED')
! open TOTENERGY.OUT
  open(61,file='TOTENERGY'//trim(filext),form='FORMATTED')
! open FERMIDOS.OUT
  open(62,file='FERMIDOS'//trim(filext),form='FORMATTED')
! open MOMENT.OUT if required
  if (spinpol) open(63,file='MOMENT'//trim(filext),form='FORMATTED')
! open GAP.OUT
  open(64,file='GAP'//trim(filext),form='FORMATTED')
! open RMSDVS.OUT
  open(65,file='RMSDVS'//trim(filext),form='FORMATTED')
! open DTOTENERGY.OUT
  open(66,file='DTOTENERGY'//trim(filext),form='FORMATTED')
! open TMDFTU.OUT
  if (tmwrite) open(67,file='TMDFTU'//trim(filext),form='FORMATTED')
! open MOMENTM.OUT
  if (spinpol) open(68,file='MOMENTM'//trim(filext),form='FORMATTED')
! open RESIDUAL.OUT
  if (xctype(1).lt.0) open(69,file='RESIDUAL'//trim(filext),form='FORMATTED')
! write out general information to INFO.OUT
  call writeinfo(60)
  write(60,*)
end if
iscl=0
if (trdstate) then
! read the Kohn-Sham potential and fields from file
  call readstate
  if (mp_mpi) then
    write(60,'("Potential read in from STATE.OUT")')
  end if
  if (autolinengy) call readfermi
else
! initialise the density and magnetisation from atomic data
  call rhoinit
  call maginit
! generate the Kohn-Sham potential and magnetic field
  call potks(.true.)
  if (mp_mpi) then
    write(60,'("Kohn-Sham potential initialised from atomic data")')
  end if
end if
if (mp_mpi) flush(60)
call genvsig
! size of mixing vector
n=npmtmax*natmtot+ngtot
if (spinpol) n=n+ndmag*(npcmtmax*natmtot+ngtot)
if (tvmatmt) n=n+2*((lmmaxdm*nspinor)**2)*natmtot
! allocate mixing array
allocate(v(n))
! determine the size of the mixer work array
nwork=-1
call mixerifc(mixtype,n,v,dv,nwork,v)
allocate(work(nwork))
! initialise the mixer
iscl=0
call mixpack(.true.,n,v)
call mixerifc(mixtype,n,v,dv,nwork,work)
! set the stop signal to .false.
tstop=.false.
! set last self-consistent loop flag
tlast=.false.
etp=0.d0
! begin the self-consistent loop
if (mp_mpi) then
  write(60,*)
  write(60,'("+------------------------------+")')
  write(60,'("| Self-consistent loop started |")')
  write(60,'("+------------------------------+")')
end if
do iscl=1,maxscl
  if (mp_mpi) then
    write(60,*)
    write(60,'("+--------------------+")')
    write(60,'("| Loop number : ",I4," |")') iscl
    write(60,'("+--------------------+")')
  end if
  if (iscl.ge.maxscl) then
    if (mp_mpi) then
      write(60,*)
      write(60,'("Reached self-consistent loops maximum")')
    end if
    write(*,*)
    write(*,'("Warning(gndstate): failed to reach self-consistency after ",I4,&
     &" loops")') iscl
    tlast=.true.
  end if
  if (mp_mpi) flush(60)
! reset the OpenMP thread variables
  call omp_reset
! generate the core wavefunctions and densities
  call gencore
! find the new linearisation energies
  call linengy
! write out the linearisation energies
  if (mp_mpi) call writelinen
! generate the APW and local-orbital radial functions and integrals
  call genapwlofr
! generate the spin-orbit coupling radial functions
  call gensocfr
! generate the first- and second-variational eigenvectors and eigenvalues
  call genevfsv
! find the occupation numbers and Fermi energy
  call occupy
  if (mp_mpi.and.autoswidth) then
    write(60,*)
    write(60,'("New smearing width : ",G18.10)') swidth
  end if
  if (mp_mpi) then
! write the occupation numbers to file
    do ik=1,nkpt
      call putoccsv(filext,ik,occsv(:,ik))
    end do
! write eigenvalues to file
    call writeeval
! write the Fermi energy to file
    call writefermi
  end if
! synchronise MPI processes
  call mpi_barrier(mpicom,ierror)
! generate the density and magnetisation
  call rhomag
! DFT+U or fixed tensor moment calculation
  if ((dftu.ne.0).or.(ftmtype.ne.0)) then
! generate the muffin-tin density matrix used for computing the potential matrix
    call gendmatmt
! write the FTM tensor moments to file
    if (ftmtype.ne.0) call writeftm
! generate the DFT+U or FTM muffin-tin potential matrices
    call genvmatmt
  end if
  if (dftu.ne.0) then
    if (mp_mpi) then
! write the DFT+U matrices to file
      call writedftu
! calculate and write tensor moments to file
      if (tmwrite) then
        if (spinorb) then
          call writetm3du(67)
        else
          call writetm2du(67)
        end if
      end if
    end if
  end if
! compute the Kohn-Sham potentials and magnetic fields
  call potks(.true.)
  if (mp_mpi) then
    if ((xcgrad.eq.3).and.(c_tb09.ne.0.d0)) then
      write(60,*)
      write(60,'("Tran-Blaha ''09 constant c : ",G18.10)') c_tb09
    end if
  end if
! pack interstitial and muffin-tin potential and field into one array
  call mixpack(.true.,n,v)
! mix in the old potential and field with the new
  call mixerifc(mixtype,n,v,dv,nwork,work)
! make sure every MPI process has a numerically identical potential
  if (np_mpi.gt.1) then
    call mpi_bcast(v,n,mpi_double_precision,0,mpicom,ierror)
  end if
! unpack potential and field
  call mixpack(.false.,n,v)
! add the fixed spin moment effective field (after mixing)
  call fsmfield
! Fourier transform Kohn-Sham potential to G-space
  call genvsig
! reduce the external magnetic fields if required
  if (reducebf.lt.1.d0) then
    bfieldc(:)=bfieldc(:)*reducebf
    bfcmt(:,:,:)=bfcmt(:,:,:)*reducebf
  end if
! compute the energy components
  call energy
  if (mp_mpi) then
! output energy components
    call writeengy(60)
    write(60,*)
    write(60,'("Density of states at Fermi energy : ",G18.10)') fermidos
    write(60,'(" (states/Hartree/unit cell)")')
    write(60,*)
    write(60,'("Estimated indirect band gap : ",G18.10)') bandgap(1)
    write(60,'(" from k-point ",I6," to k-point ",I6)') ikgap(1),ikgap(2)
    write(60,'("Estimated direct band gap   : ",G18.10)') bandgap(2)
    write(60,'(" at k-point ",I6)') ikgap(3)
! write total energy to TOTENERGY.OUT
    write(61,'(G22.12)') engytot
    flush(61)
! write DOS at Fermi energy to FERMIDOS.OUT
    write(62,'(G18.10)') fermidos
    flush(62)
! output charges and moments
    call writechg(60)
    if (spinpol) then
! write total moment to MOMENT.OUT
      write(63,'(3G18.10)') momtot(1:ndmag)
      flush(63)
! write total moment magnitude to MOMENTM.OUT
      write(68,'(G18.10)') momtotm
      flush(68)
    end if
! write estimated Kohn-Sham indirect band gap
    write(64,'(G22.12)') bandgap(1)
    flush(64)
! output effective fields for fixed spin moment calculations
    if (fsmtype.ne.0) call writefsm(60)
! check for WRITE file
    inquire(file='WRITE',exist=exist)
    if (exist) then
      write(60,*)
      write(60,'("WRITE file exists - writing STATE.OUT")')
      call writestate
      open(50,file='WRITE')
      close(50,status='DELETE')
    end if
! write STATE.OUT file if required
    if (nwrite.ge.1) then
      if (mod(iscl,nwrite).eq.0) then
        call writestate
        write(60,*)
        write(60,'("Wrote STATE.OUT")')
      end if
    end if
! write OEP residual
    if (xctype(1).lt.0) then
      write(60,*)
      write(60,'("Magnitude of OEP residual : ",G18.10)') resoep
      write(69,'(G18.10)') resoep
      flush(69)
    end if
  end if
! exit self-consistent loop if required
  if (tlast) goto 10
! check for convergence
  if (iscl.ge.2) then
    if (mp_mpi) then
      write(60,*)
      write(60,'("RMS change in Kohn-Sham potential (target) : ",G18.10," (",&
       &G18.10,")")') dv,epspot
      write(65,'(G18.10)') dv
      flush(65)
    end if
    de=abs(engytot-etp)
    if (mp_mpi) then
      write(60,'("Absolute change in total energy (target)   : ",G18.10," (",&
       &G18.10,")")') de,epsengy
      write(66,'(G18.10)') de
      flush(66)
    end if
    if ((dv.lt.epspot).and.(de.lt.epsengy)) then
      if (mp_mpi) then
        write(60,*)
        write(60,'("Convergence targets achieved")')
      end if
      tlast=.true.
    end if
  end if
! average the current and previous total energies and store
  if (iscl.gt.1) then
    etp=0.75d0*engytot+0.25d0*etp
  else
    etp=engytot
  end if
! check for STOP file (only master process)
  if (mp_mpi) then
    inquire(file='STOP',exist=exist)
    if (exist) then
      write(60,*)
      write(60,'("STOP file exists - stopping self-consistent loop")')
      open(50,file='STOP')
      close(50,status='DELETE')
      tstop=.true.
      tlast=.true.
    end if
  end if
! broadcast tlast and tstop from master process to all other processes
  call mpi_bcast(tlast,1,mpi_logical,0,mpicom,ierror)
  call mpi_bcast(tstop,1,mpi_logical,0,mpicom,ierror)
! output the current total CPU time
  timetot=timeinit+timemat+timefv+timesv+timerho+timepot+timefor
  if (mp_mpi) then
    write(60,*)
    write(60,'("Time (CPU seconds) : ",F12.2)') timetot
  end if
! end the self-consistent loop
end do
10 continue
! synchronise MPI processes
call mpi_barrier(mpicom,ierror)
if (mp_mpi) then
  write(60,*)
  write(60,'("+------------------------------+")')
  write(60,'("| Self-consistent loop stopped |")')
  write(60,'("+------------------------------+")')
! write density and potentials to file only if maxscl > 1
  if (maxscl.gt.1) then
    call writestate
    write(60,*)
    write(60,'("Wrote STATE.OUT")')
  end if
end if
! compute forces if required
if (tforce) then
  call force
! output forces to INFO.OUT
  if (mp_mpi) call writeforces(60)
end if
! compute the current density and total current if required
if (tcden) then
  call curden(afieldc)
  if (mp_mpi) then
    write(60,*)
    write(60,'("Total current per unit cell")')
    write(60,'(3G18.10)') curtot
    write(60,'(" magnitude : ",G18.10)') curtotm
  end if
end if
! total time used
timetot=timeinit+timemat+timefv+timesv+timerho+timepot+timefor
if (mp_mpi) then
! output timing information
  write(60,*)
  write(60,'("Timings (CPU seconds) :")')
  write(60,'(" initialisation",T40,": ",F12.2)') timeinit
  write(60,'(" Hamiltonian and overlap matrix set up",T40,": ",F12.2)') timemat
  write(60,'(" first-variational eigenvalue equation",T40,": ",F12.2)') timefv
  if (spinpol) then
    write(60,'(" second-variational calculation",T40,": ",F12.2)') timesv
  end if
  write(60,'(" charge density calculation",T40,": ",F12.2)') timerho
  write(60,'(" potential calculation",T40,": ",F12.2)') timepot
  if (tforce) then
    write(60,'(" force calculation",T40,": ",F12.2)') timefor
  end if
  write(60,'(" total",T40,": ",F12.2)') timetot
  write(60,*)
  write(60,'("+----------------------------+")')
  write(60,'("| Elk version ",I1.1,".",I1.1,".",I2.2," stopped |")') version
  write(60,'("+----------------------------+")')
! close the INFO.OUT file
  close(60)
! close the TOTENERGY.OUT file
  close(61)
! close the FERMIDOS.OUT file
  close(62)
! close the MOMENT.OUT and MOMENTM.OUT files
  if (spinpol) then
    close(63); close(68)
  end if
! close the GAP.OUT file
  close(64)
! close the RMSDVS.OUT file
  close(65)
! close the DTOTENERGY.OUT file
  close(66)
! close TMDFTU.OUT file
  if (tmwrite) close(67)
! close the RESIDUAL.OUT file
  if (xctype(1).lt.0) close(69)
end if
deallocate(v,work)
! synchronise MPI processes
call mpi_barrier(mpicom,ierror)
return
end subroutine
!EOC