← Back to branch summary

~siesta-pseudos-bases/siesta/trunk-psml

« back to all changes in this revision

Viewing changes to Util/Eig2DOS/Eig2DOS.f90

  • Committer: Alberto Garcia
  • Date: 2019-09-02 14:09:43 UTC
  • mfrom: (427.6.323 trunk)
  • Revision ID: albertog@icmab.es-20190902140943-mzmbe1jacgefpgxw
Sync to trunk-776 (notably nc/soc wavefunction support)


Show diffs side-by-side

added added

removed removed

Lines of Context:
Util/Eig2DOS/Eig2DOS.f90
26
26
!  -e/-m, -E/-M - energy window: Emin and Emax
27
27
!  -F/-f  - shift E_F to zero
28
28
!
29
 
! Density of states in standard output, for both spins.
30
 
! If nspin = 1, the first spin component is multiplied by 2.
31
 
!
32
29
!   NOT printed in this version:
33
30
!   the number of electrons (states) in the energy window, 
34
31
!   with and without broadening (number of eigenvalues and 
44
41
  integer :: nargs, iostat, n_opts, nlabels
45
42
 
46
43
  integer :: nk, nspin, nband, ie, ik, ika, is, ib, nel
47
 
  integer :: nk_kpoints, ik_read
 
44
  integer :: nk_kpoints, ik_read, nspin_blocks
48
45
 
49
46
  real(dp) :: e, eincr, ef, pi, x, sta, norm
50
47
  real(dp) :: emin_file, emax_file
56
53
  character(len=256) :: eig_file, kpoint_file
57
54
 
58
55
  integer  :: npts_energy = 200
59
 
  real(dp) :: emin        = huge(1.0_dp)
60
 
  real(dp) :: emax        = -huge(1.0_dp)
 
56
  real(dp) :: emin        = -huge(1.0_dp)
 
57
  real(dp) :: emax        =  huge(1.0_dp)
61
58
  real(dp) :: smear       = 0.2_dp
62
59
  logical  :: loren       = .false.
63
60
  logical  :: emin_given  = .false.
64
61
  logical  :: emax_given  = .false.
65
62
  logical  :: using_weights = .false.
66
63
  logical  :: shift_efermi = .false.
 
64
  logical  :: non_coll
67
65
  integer  :: min_band = 0
68
66
  integer  :: max_band = 0
69
67
 
157
155
       action="read")
158
156
  read(1,*) Ef
159
157
  read(1,*) nband, nspin, nk
 
158
  non_coll = (nspin == 4)
 
159
 
 
160
  if (non_coll) then
 
161
     nspin_blocks = 1
 
162
  else
 
163
     nspin_blocks = nspin
 
164
  endif
 
165
  
160
166
  write(*,"(a)") "# Eigenvalues read from " // trim(eig_file)
161
167
  if ( debug ) print *, "Ef, nband, nspin, nk:", Ef, nband, nspin, nk
162
168
 
171
177
     write(*,"(a)") "# Kpoint weights read from " // trim(kpoint_file)
172
178
  endif
173
179
 
174
 
  allocate(eig(nband,nspin))
175
 
  allocate(DOS(npts_energy,nspin))
 
180
  allocate(eig(nband,nspin_blocks))
 
181
  allocate(DOS(npts_energy,nspin_blocks))
176
182
 
177
183
  ! Sanity checks
178
184
 
198
204
  emin_file = huge(1.0_dp)
199
205
  emax_file = -huge(1.0_dp)
200
206
  do ik = 1, nk
201
 
     read(1,*) ika, ((eig(ib,is), ib = 1, nband), is = 1, nspin)
202
 
     emin_file = min(emin_file,minval(eig(min_band:max_band,1:nspin)))
203
 
     emax_file = max(emax_file,maxval(eig(min_band:max_band,1:nspin)))
 
207
     read(1,*) ika, ((eig(ib,is), ib = 1, nband), is = 1, nspin_blocks)
 
208
     emin_file = min(emin_file,minval(eig(min_band:max_band,1:nspin_blocks)))
 
209
     emax_file = max(emax_file,maxval(eig(min_band:max_band,1:nspin_blocks)))
204
210
  end do
205
211
 
206
212
  write(*,"(a,2f15.7)") "# Emin, Emax in file for selected band(s):", emin_file, emax_file
222
228
  if (.not. emax_given) emax = emax_file + 6._dp*smear
223
229
 
224
230
  if (npts_energy .lt. 2) npts_energy = 2
225
 
  eincr = (emax - emin) / real((npts_energy - 1),kind=dp)
 
231
  eincr = (emax - emin) / real( npts_energy - 1, dp)
226
232
 
227
233
  nel = 0
228
234
  DOS(:,:) = 0.0_dp
231
237
 
232
238
  do ik = 1, nk
233
239
 
234
 
     read(1,*) ika, ((eig(ib,is), ib = 1, nband), is = 1, nspin)
 
240
     read(1,*) ika, ((eig(ib,is), ib = 1, nband), is = 1, nspin_blocks)
235
241
     if ( using_weights ) then
236
242
        read(2,*) ik_read, k(:), weight
237
243
        if (ik_read /= ik) STOP "ik mismatch"
239
245
        weight = 1.0_dp / nk
240
246
     end if
241
247
 
242
 
     do is = 1, nspin
 
248
     do is = 1, nspin_blocks
243
249
        do ib = min_band, max_band
244
250
           e = eig(ib,is)
245
251
           if (shift_efermi) e = e - ef
262
268
     norm = sqrt(pi) * smear
263
269
  end if
264
270
 
265
 
  do is = 1, nspin
 
271
  do is = 1, nspin_blocks
266
272
     do ie = 1, npts_energy
267
273
        DOS(ie,is) = DOS(ie,is)/norm
268
274
     end do
270
276
 
271
277
! integral, extremely sophisticated -----------------------------------
272
278
 
273
 
  do is = 1, nspin
 
279
  do is = 1, nspin_blocks
274
280
     integral(is) = 0.0_dp
275
281
     do ie = 1, npts_energy
276
282
        integral(is) = integral(is) + dos(ie,is)*eincr
282
288
  if ( nspin == 1 ) then
283
289
     nel = nel * 2
284
290
  end if
285
 
  sta = nel / real(nk,kind=dp)
 
291
  sta = nel / real(nk, dp)
286
292
 
287
293
! output, prepared for gnuplot ----------------------------------------
288
294
 
 
295
! Density of states in standard output:
 
296
  
 
297
! -- If nspin = 2, the first two columns correspond to the two spin components.
 
298
!                  A third column contains their sum.
 
299
! -- If nspin = 1, two identical columns are printed, plus their sum.
 
300
! -- If nspin = 4 (NC/SOC case), a single column is printed with the complete DOS.
 
301
!
 
302
 
289
303
  write(*,"(a,3i6)") '# Nbands, Nspin, Nk   = ', nband, nspin, nk
290
304
  if (shift_efermi) then
291
305
     write(*,"(a,f10.4,a)") '# E_F                 = ', Ef , ' eV --> (shifted to ZERO)'
302
316
     do ie = 1, npts_energy
303
317
        write(*,"(4f14.6)") emin + (ie-1)*eincr, DOS(ie,1), DOS(ie,1), 2._dp*DOS(ie,1)
304
318
     end do
305
 
  else
 
319
  else if ( nspin == 2 ) then
306
320
     write(*,"(a)") '#        E            N(up)        N(down)         Ntot'
307
321
     do ie = 1, npts_energy
308
322
        write(*,"(4f14.6)") emin + (ie-1)*eincr, DOS(ie,:), sum(DOS(ie,:))
309
323
     end do
 
324
  else  ! nspin == 4   !!  Single column with complete DOS
 
325
     write(*,"(a)") '#        E            Ntot  (complete DOS for NC/SOC case)'
 
326
     do ie = 1, npts_energy
 
327
        write(*,"(4f14.6)") emin + (ie-1)*eincr, DOS(ie,:)
 
328
     end do
310
329
  end if
311
330
 
312
331
  deallocate(eig,dos)