← Back to branch summary

~albertog/siesta/efield-2.5-jjunquera

« back to all changes in this revision

Viewing changes to Util/COOP/mprop.f90

  • Committer: Alberto Garcia
  • Date: 2007-11-22 21:31:50 UTC
  • mfrom: (unknown (missing))
  • Revision ID: Arch-1:siesta@uam.es--2006%siesta-devel--reference--2.5--patch-1
Merged COOP/COHP functionality
Merged the COOP/COHP/PDOS code from the 2.1 branch complex.

��-- Functionality changes:

* The wavefunction file is now in WFSX format (well-packed, single-precision, and thus
significantly smaller). The Util/wfsx2wfs program can be used to re-generate the old format.

* The xijo array (for relative positions between interacting orbitals) is always produced.

-- New features:

Option COOP.Write triggers the creation of a wavefunction file (for the SCF-k-point set)
and an HSX file (enhanced HS format) that can later be processed by Util/COOP/mprop
for the off-line generation of COOP/COHP/(P)DOS.

(replayed patch-1 of outlier branch ref-2.4)


Show diffs side-by-side

added added

removed removed

Lines of Context:
Util/COOP/mprop.f90
 
1
!==================================================================
 
2
program mprop
 
3
 
 
4
  use main_vars
 
5
  use orbital_set, only: get_orbital_set
 
6
  use io_hs, only: read_hs_file
 
7
  use read_curves, only: read_curve_information, mask_to_arrays
 
8
 
 
9
  implicit none
 
10
 
 
11
  logical :: gamma_wfsx
 
12
 
 
13
  !
 
14
  !     Process options
 
15
  !
 
16
  n_opts = 0
 
17
  do
 
18
     call getopts('dhls:n:m:M:R:',opt_name,opt_arg,n_opts,iostat)
 
19
     if (iostat /= 0) exit
 
20
     select case(opt_name)
 
21
     case ('d')
 
22
        debug = .true.
 
23
     case ('+d')
 
24
        debug = .false.
 
25
     case ('l')
 
26
        energies_only = .true.
 
27
     case ('s')
 
28
        read(opt_arg,*) smear
 
29
     case ('n')
 
30
        read(opt_arg,*) npts_energy
 
31
     case ('m')
 
32
        read(opt_arg,*) minimum_spec_energy
 
33
     case ('M')
 
34
        read(opt_arg,*) maximum_spec_energy
 
35
     case ('R')
 
36
        ref_line_given = .true.
 
37
        ref_line = opt_arg
 
38
     case ('h')
 
39
        call manual()
 
40
     case ('?',':')
 
41
        write(0,*) "Invalid option: ", opt_arg(1:1)
 
42
        write(0,*) "Usage: mprop [ -d ] [ -h ] MPROP_FILE_ROOT"
 
43
        write(0,*) "Use -h option for manual"
 
44
        STOP
 
45
     end select
 
46
  enddo
 
47
 
 
48
  nargs = command_argument_count()
 
49
  nlabels = nargs - n_opts + 1
 
50
  if (nlabels /= 1)  then
 
51
     write(0,*) "Usage: mprop [ -d ] [ -h ] MPROP_FILE_ROOT"
 
52
     write(0,*) "Use -h option for manual"
 
53
     STOP
 
54
  endif
 
55
 
 
56
  call get_command_argument(n_opts,value=mflnm,status=iostat)
 
57
  if (iostat /= 0) then
 
58
     STOP "Cannot get .mprop file root"
 
59
  endif
 
60
 
 
61
  print "(a,f7.3)", "Using smearing parameter: ", smear
 
62
  print "(a,i6,a)", "Using ", npts_energy, " points in energy range"
 
63
 
 
64
  !==================================================
 
65
 
 
66
  ierr=0
 
67
 
 
68
  ! * LECTURA TIPUS CALCUL
 
69
 
 
70
  open(mpr_u,file=trim(mflnm) // ".mpr", status='old')
 
71
  read(mpr_u,*) sflnm
 
72
  read(mpr_u,*) what
 
73
 
 
74
  dos=(trim(what).eq.'DOS')
 
75
  coop=(trim(what).eq.'COOP')
 
76
  
 
77
  !==================================================
 
78
  ! * LECTURA sflnm.WFSX
 
79
 
 
80
  write(6,"(1x,a,'.WFSX ...')") trim(sflnm)
 
81
 
 
82
  open(wfs_u,file=trim(sflnm)//'.WFSX',status='old',form='unformatted')
 
83
  read(wfs_u) nkp, gamma_wfsx
 
84
  allocate (wk(nkp), pk(3,nkp))
 
85
 
 
86
  read(wfs_u) nsp
 
87
  read(wfs_u) nao
 
88
  read(wfs_u)        !! Symbols, etc
 
89
  if (debug) print *, "WFSX read: nkp, nsp, nnao: ", nkp, nsp, nao
 
90
 
 
91
  allocate (ados(npts_energy,nsp))
 
92
 
 
93
  nwfmx = 0
 
94
  min_energy = huge(1.0_dp)
 
95
  max_energy = -huge(1.0_dp)
 
96
 
 
97
  do ik=1,nkp
 
98
     do is=1,nsp
 
99
 
 
100
        read(wfs_u) idummy, pk(1:3,ik), wk(ik)
 
101
        if (idummy /= ik) stop "ik index mismatch in WFS file"
 
102
 
 
103
        read(wfs_u) is0
 
104
        read(wfs_u) number_of_wfns
 
105
        nwfmx = max(nwfmx,number_of_wfns)
 
106
 
 
107
        do iw=1,number_of_wfns
 
108
           read(wfs_u) iw0
 
109
           read(wfs_u) eigval
 
110
           min_energy = min(min_energy,eigval)
 
111
           max_energy = max(max_energy,eigval)
 
112
           read(wfs_u)
 
113
        enddo
 
114
     enddo
 
115
  enddo
 
116
 
 
117
  print *, " Maximum number of wfs per k-point: ", nwfmx
 
118
  print "(a,2f12.4)", "Min_energy, max_energy on WFS file: ",  &
 
119
       min_energy, max_energy
 
120
 
 
121
 
 
122
 
 
123
  low_e = min_energy - 2*smear
 
124
  high_e = max_energy + 2*smear
 
125
  e_step = (high_e-low_e)/(npts_energy-1)
 
126
  ados(:,1:nsp) = 0.0_dp
 
127
 
 
128
  ! skip four records
 
129
 
 
130
  rewind(wfs_u)
 
131
 
 
132
  read(wfs_u) 
 
133
  read(wfs_u) 
 
134
  read(wfs_u) 
 
135
  read(wfs_u) 
 
136
 
 
137
  do ik=1,nkp
 
138
     do is=1,nsp
 
139
        read(wfs_u)
 
140
        read(wfs_u)
 
141
        read(wfs_u)  number_of_wfns
 
142
        do iw=1,number_of_wfns
 
143
           read(wfs_u) 
 
144
           read(wfs_u) eigval
 
145
           do i = 1, npts_energy
 
146
              energy = low_e + e_step*(i-1)
 
147
              weight = delta(energy-eigval)
 
148
              if (weight < 1.0e-6) CYCLE           ! Will not contribute
 
149
              ados(i,is) = ados(i,is) + wk(ik) * weight
 
150
           enddo
 
151
           read(wfs_u)       ! Skip wfn info
 
152
        enddo
 
153
     enddo
 
154
  enddo
 
155
 
 
156
  call io_assign(idos)
 
157
  open(idos,file=trim(sflnm)//".alldos",form="formatted", &
 
158
       status="unknown",action="write",position="rewind")
 
159
  write(idos,*) "#  Energy   LARGE-SCALE DOS"
 
160
  do i = 1, npts_energy
 
161
     energy = low_e + e_step*(i-1)
 
162
     write(idos,*) energy, (ados(i,is), is=1,nsp)
 
163
  enddo
 
164
 
 
165
  call io_close(idos)
 
166
 
 
167
 
 
168
  ! * LECTURA sflnm.HS
 
169
  ! Will pick up atoms, zval, and thus N_electrons
 
170
 
 
171
  call read_hs_file()
 
172
  if (gamma_wfsx .neqv. gamma) STOP "Gamma mismatch"
 
173
 
 
174
  ztot = 0.0_dp
 
175
  do ia = 1, na_u
 
176
     ztot = ztot + zval(isa(ia))
 
177
  enddo
 
178
 
 
179
 
 
180
  !
 
181
  !       Compute integrated total DOS
 
182
 
 
183
  allocate(intdos(npts_energy))
 
184
  call io_assign(intdos_u)
 
185
  open(intdos_u,file=trim(sflnm)//".intdos",form="formatted", &
 
186
       status="unknown",action="write",position="rewind")
 
187
  intdos(1) = 0.0_dp
 
188
  write(intdos_u,*) low_e, intdos(1)
 
189
  do i = 2, npts_energy
 
190
     energy = low_e + e_step*(i-1)
 
191
     intdos(i) = intdos(i-1) + sum(ados(i,:)) * e_step * 2.0_dp /nsp
 
192
     write(intdos_u,*) energy, intdos(i)
 
193
  enddo
 
194
  call io_close(intdos_u)
 
195
 
 
196
  ! Look for Fermi Energy
 
197
  do i = 2, npts_energy
 
198
     if (intdos(i) > ztot) then
 
199
        ! Found fermi energy
 
200
        energy = low_e + e_step*(i-1)
 
201
        efermi = energy - (intdos(i)-ztot)*e_step/(intdos(i)-intdos(i-1))
 
202
        exit
 
203
     endif
 
204
  enddo
 
205
 
 
206
  print *, "Fermi energy: ", efermi
 
207
 
 
208
  if (energies_only) STOP
 
209
  !-------------------------------------------------------------------
 
210
 
 
211
  if (minimum_spec_energy > min_energy) then
 
212
     min_energy = minimum_spec_energy
 
213
  else
 
214
     min_energy = min_energy - 2*smear
 
215
  endif
 
216
  if (maximum_spec_energy < max_energy) then
 
217
     max_energy = maximum_spec_energy
 
218
  else
 
219
     max_energy = max_energy + 2*smear
 
220
  endif
 
221
 
 
222
  print "(a,2f12.4)", "Min_energy, max_energy used: ",  &
 
223
       min_energy, max_energy
 
224
 
 
225
  !====================
 
226
 
 
227
  ! * Orbital list
 
228
 
 
229
  allocate(za(no_u), zc(no_u), zn(no_u), zl(no_u), zx(no_u), zz(no_u))
 
230
  nao = 0
 
231
  do ia=1,na_u
 
232
     it = isa(ia)
 
233
     io = 0
 
234
     do 
 
235
        io = io + 1
 
236
        if (io > no(it)) exit
 
237
        lorb = lquant(it,io)
 
238
        do ko = 1, 2*lorb + 1
 
239
           nao = nao + 1
 
240
           za(nao)=ia
 
241
           zc(nao)=it
 
242
           zn(nao)=nquant(it,io)
 
243
           zl(nao)=lorb
 
244
           zx(nao)=ko
 
245
           zz(nao)=zeta(it,io)
 
246
        enddo
 
247
        io = io + 2*lorb
 
248
     enddo
 
249
  enddo
 
250
  if (nao /= no_u) STOP "nao /= no_u"
 
251
 
 
252
  ! ==================================
 
253
 
 
254
  write(6,"('Writing files: ',a,'.stt ...')",advance='no') trim(mflnm)
 
255
  open(stt_u,file=trim(mflnm)//'.stt')
 
256
  write(stt_u,"(/'UNIT CELL ATOMS:')")
 
257
  write(stt_u,"(3x,i4,2x,i3,2x,a2)") (i, isa(i), label(isa(i)), i=1,na_u)
 
258
  write(stt_u,"(/'BASIS SET:')")
 
259
  write(stt_u,"(5x,a20,3(3x,a1))") 'spec', 'n', 'l', 'z'
 
260
  do it=1,nspecies
 
261
     write(stt_u,"(5x,a20)") trim(label(it))
 
262
     io = 0
 
263
     do 
 
264
        io = io + 1
 
265
        if (io > no(it)) exit
 
266
        write(stt_u,"(3(2x,i2))") nquant(it,io), lquant(it,io), zeta(it,io)
 
267
        io = io + 2*lquant(it,io)
 
268
     enddo
 
269
  enddo
 
270
 
 
271
  write(stt_u,"(/'SPIN: ',i2)") nsp
 
272
  write(stt_u,"(/'FERMI ENERGY: ',f18.6)") efermi
 
273
 
 
274
  write(stt_u,"(/'AO LIST:')")
 
275
  taux=repeat(' ',len(taux))
 
276
  do io=1,no_u
 
277
     taux(1:30)=repeat(' ',30)
 
278
     ik=1
 
279
     if (zl(io).eq.0) then
 
280
        taux(ik:ik)='s'
 
281
        ik=0
 
282
     elseif (zl(io).eq.1) then
 
283
        if (zx(io).eq.1) taux(ik:ik+1)='py'
 
284
        if (zx(io).eq.2) taux(ik:ik+1)='pz'
 
285
        if (zx(io).eq.3) taux(ik:ik+1)='px'
 
286
        ik=0
 
287
     elseif (zl(io).eq.2) then
 
288
        if (zx(io).eq.1) taux(ik:ik+2)='dxy'
 
289
        if (zx(io).eq.2) taux(ik:ik+2)='dyz'
 
290
        if (zx(io).eq.3) taux(ik:ik+2)='dz2'
 
291
        if (zx(io).eq.4) taux(ik:ik+2)='dxz'
 
292
        if (zx(io).eq.5) taux(ik:ik+5)='dx2-y2'
 
293
        ik=0
 
294
     elseif (zl(io).eq.3) then
 
295
        taux(ik:ik)='f'
 
296
     elseif (zl(io).eq.4) then
 
297
        taux(ik:ik)='g'
 
298
     elseif (zl(io).eq.5) then
 
299
        taux(ik:ik)='h'
 
300
     endif
 
301
     write(stt_u,"(3x,i5,2x,i3,2x,a2)",advance='no')  &
 
302
                          io, za(io), trim(label(zc(io)))
 
303
     if (ik.eq.0) then
 
304
        write(stt_u,"(3x,i2,a)") zn(io), trim(taux)
 
305
     else
 
306
        write(stt_u,"(3x,i2,a,i2.2)") zn(io), trim(taux), zx(io)
 
307
     endif
 
308
  enddo
 
309
 
 
310
  !==================================
 
311
 
 
312
  if (ref_line_given) then
 
313
     allocate(ref_mask(no_u))
 
314
     print *, "Orbital set spec: ", trim(ref_line)
 
315
     call get_orbital_set(ref_line,ref_mask)
 
316
     do io=1, no_u
 
317
        if (ref_mask(io)) write(6,fmt="(i5)",advance="no") io
 
318
     enddo
 
319
     deallocate(ref_mask)
 
320
     write(6,*)
 
321
     STOP "bye from ref_line processing"
 
322
  endif
 
323
 
 
324
!
 
325
! Process orbital sets
 
326
!
 
327
  allocate(orb_mask(no_u,2,ncbmx))
 
328
  allocate (koc(ncbmx,2,no_u))
 
329
  call read_curve_information(dos,coop,  &
 
330
                                    mpr_u,no_u,ncbmx,ncb,tit,orb_mask,dtc)
 
331
  if (dos) then
 
332
     orb_mask(:,2,1:ncb) = .true.       ! All orbitals considered
 
333
  endif
 
334
  call mask_to_arrays(ncb,orb_mask(:,1,:),noc(:,1),koc(:,1,:))
 
335
  call mask_to_arrays(ncb,orb_mask(:,2,:),noc(:,2),koc(:,2,:))
 
336
 
 
337
 
 
338
  !=====================
 
339
 
 
340
 if (coop) then
 
341
    allocate (coop_vals(npts_energy,nsp,ncb))
 
342
    allocate (cohp_vals(npts_energy,nsp,ncb))
 
343
    coop_vals(:,:nsp,:ncb)=0.d0
 
344
    cohp_vals(:,:nsp,:ncb)=0.d0
 
345
 endif
 
346
 if (dos) then
 
347
    allocate (pdos_vals(npts_energy,nsp,ncb))
 
348
    pdos_vals(:,:nsp,:ncb) = 0.0_dp
 
349
 endif
 
350
 
 
351
  ados(:,:nsp) = 0.0_dp
 
352
 
 
353
  !================================================================
 
354
 
 
355
  ! * Curves
 
356
 
 
357
     if (gamma) then
 
358
        allocate(wf(1,1:no_u))
 
359
     else
 
360
        allocate(wf(2,1:no_u))
 
361
     endif
 
362
 
 
363
     !Stream over file, without using too much memory
 
364
 
 
365
     rewind(wfs_u)
 
366
 
 
367
     read(wfs_u) 
 
368
     read(wfs_u) 
 
369
     read(wfs_u) 
 
370
     read(wfs_u) 
 
371
     
 
372
     e_step = (max_energy-min_energy)/(npts_energy-1)
 
373
 
 
374
     do ik=1,nkp
 
375
        do is=1,nsp
 
376
           read(wfs_u)
 
377
           read(wfs_u)
 
378
           read(wfs_u)  number_of_wfns
 
379
           do iw=1,number_of_wfns
 
380
              read(wfs_u) 
 
381
              read(wfs_u) eigval
 
382
              ! Early termination of iteration
 
383
              if (eigval < min_energy .or. eigval > max_energy) then
 
384
                 read(wfs_u)   ! Still need to read this
 
385
                 CYCLE
 
386
              endif
 
387
 
 
388
              read(wfs_u) (wf(:,io), io=1,nao)
 
389
 
 
390
              do i = 1, npts_energy
 
391
                 energy = min_energy + e_step*(i-1)
 
392
                 weight = delta(energy-eigval)
 
393
                 if (weight < 1.0e-6) CYCLE           ! Will not contribute
 
394
 
 
395
                 ados(i,is) = ados(i,is) + wk(ik) * weight
 
396
 
 
397
                 do ic=1,ncb
 
398
 
 
399
                    do i1=1,noc(ic,1)
 
400
                      io1=koc(ic,1,i1)              ! AO Set I
 
401
 
 
402
                      ! On site contribution to DOS  
 
403
                      if (dos) then
 
404
                        if (gamma) then
 
405
                          qcos = wf(1,io1)**2
 
406
                        else
 
407
                          qcos=(wf(1,io1)**2+wf(2,io1)**2)
 
408
                        endif
 
409
 
 
410
                        pdos_vals(i,is,ic)= pdos_vals(i,is,ic)  +   &
 
411
                                               qcos * wk(ik) * weight
 
412
                      endif
 
413
 
 
414
                       do i2=1,noc(ic,2)             
 
415
                          io2=koc(ic,2,i2)              ! AO Set II
 
416
 
 
417
                          ! (qcos, qsin) = C_1*conjg(C_2)
 
418
 
 
419
                          if (gamma) then
 
420
                             qcos = wf(1,io1)*wf(1,io2) 
 
421
                             qsin = 0.0_dp
 
422
                          else
 
423
                             qcos= (wf(1,io1)*wf(1,io2) + &
 
424
                                  wf(2,io1)*wf(2,io2))
 
425
                             qsin= (wf(2,io1)*wf(1,io2) - &
 
426
                                  wf(1,io1)*wf(2,io2))
 
427
                          endif
 
428
 
 
429
                          do isr=1,nsr(io1,io2)
 
430
                             !
 
431
                             ! Note that we  want orbitals in
 
432
                             ! *different* atoms (See Dronskowski)
 
433
                             !
 
434
                             if ( dt(io1,io2,isr) < 1.0e-2_dp) CYCLE
 
435
 
 
436
                             ! Discard if overlap is very small
 
437
                             ! (It happens if not neglecting KB overlap)
 
438
                             if ( abs(sr(io1,io2,isr)) < 1.0e-10_dp) CYCLE
 
439
 
 
440
 
 
441
                             ! No distance restrictions for PDOS
 
442
                             ! Note also that there is no factor of 2
 
443
 
 
444
                             ! k*R_12    (r_2-r_1)
 
445
                             alfa=dot_product(pk(1:3,ik),rn(1:3,io1,io2,isr))
 
446
 
 
447
                            if (dos) then
 
448
                              pdos_vals(i,is,ic)=pdos_vals(i,is,ic)  +   & 
 
449
                                  (qcos*cos(alfa)+qsin*sin(alfa))    *   &
 
450
                                  sr(io1,io2,isr)*wk(ik) * weight
 
451
                            endif
 
452
 
 
453
                            if (coop) then
 
454
                             if (dt(io1,io2,isr).ge.dtc(ic,1) .and.  &
 
455
                                  dt(io1,io2,isr).le.dtc(ic,2)) then
 
456
                                !
 
457
                                ! Crb = Real(C_1*conjg(C_2)*exp(-i*alfa)) * S_12
 
458
                                ! or Real(conjg(C_1)*C_2)*exp(+i*alfa)) * S_12
 
459
                                !
 
460
                                coop_vals(i,is,ic)=coop_vals(i,is,ic)  +   & 
 
461
                                     2.0_dp * (qcos*cos(alfa)+qsin*sin(alfa)) * &
 
462
                                     sr(io1,io2,isr)*wk(ik) * weight
 
463
                                cohp_vals(i,is,ic)=cohp_vals(i,is,ic)  +   & 
 
464
                                     2.0_dp * (qcos*cos(alfa)+qsin*sin(alfa)) * &
 
465
                                     hr(io1,io2,isr,is)*wk(ik) * weight
 
466
                             endif
 
467
                            endif  ! coop
 
468
 
 
469
                          enddo ! isr
 
470
                       enddo  ! i2
 
471
                    enddo  ! i1
 
472
                 enddo    ! ic
 
473
              enddo  ! energy
 
474
           enddo   ! iwf
 
475
        enddo      ! is
 
476
     enddo         ! ik
 
477
 
 
478
!--------------------------------------------------------
 
479
 
 
480
     write(stt_u,"(/'KPOINTS:',i7)") nkp
 
481
     do ik=1,nkp
 
482
        write(stt_u,"(3x,3f9.6)") pk(:,ik)
 
483
     enddo
 
484
 
 
485
  if (dos) then
 
486
     write(stt_u,"(/'PDOS CURVES:')")
 
487
     do ic=1,ncb
 
488
        write(stt_u,"(3x,a)") trim(tit(ic))
 
489
        write(stt_u,"(3x,'AO set I: ',1000i5)") (koc(ic,1,j),j=1,noc(ic,1))
 
490
     enddo
 
491
  endif  ! dos
 
492
 
 
493
  if (coop) then
 
494
     write(stt_u,"(/'COOP CURVES:')")
 
495
     do ic=1,ncb
 
496
        write(stt_u,"(3x,a)") trim(tit(ic))
 
497
        write(stt_u,"(3x,'Distance range:',2f9.4)") dtc(ic,:)
 
498
        write(stt_u,"(3x,'AO set I: ',1000i5)") (koc(ic,1,j),j=1,noc(ic,1))
 
499
        write(stt_u,"(3x,'AO set II:',1000i5)") (koc(ic,2,j),j=1,noc(ic,2))
 
500
     enddo
 
501
  endif
 
502
 
 
503
  close(stt_u)
 
504
 
 
505
     !
 
506
     !      Simple DOS output
 
507
     !
 
508
     call io_assign(idos)
 
509
     open(idos,file=trim(sflnm)//".ados",form="formatted", &
 
510
          status="unknown",action="write",position="rewind")
 
511
     write(idos,*) "#  Energy   SIMPLE DOS"
 
512
     e_step = (max_energy-min_energy)/(npts_energy-1)
 
513
     do i = 1, npts_energy
 
514
        energy = min_energy + e_step*(i-1)
 
515
        write(idos,*) energy, (ados(i,is), is=1,nsp)
 
516
     enddo
 
517
     call io_close(idos)
 
518
 
 
519
   if (coop) then
 
520
     !
 
521
     !      COOP output
 
522
     !
 
523
     do ic = 1, ncb
 
524
        open(tab_u,file=trim(mflnm)// "." // trim(tit(ic)) // '.coop')
 
525
        write(tab_u,"(a1,14x,'ENERGY',4(16x,a1,i1))") '#', ("s",m, m=1,nsp)
 
526
        do i=1, npts_energy
 
527
           energy = min_energy + e_step*(i-1)
 
528
           write(tab_u,"(f20.8,10(2f13.8,5x)))")  &           ! Spin in loop
 
529
                energy, (coop_vals(i,is,ic),is=1,nspin)
 
530
        enddo
 
531
        close(tab_u)
 
532
     enddo
 
533
     !
 
534
     !      COHP output
 
535
     !
 
536
     do ic = 1, ncb
 
537
        open(tab_u,file=trim(mflnm)// "." // trim(tit(ic)) // '.cohp')
 
538
        write(tab_u,"(a1,14x,'ENERGY',4(16x,a1,i1))") '#', ("s",m, m=1,nsp)
 
539
        do i=1, npts_energy
 
540
           energy = min_energy + e_step*(i-1)
 
541
           write(tab_u,"(f20.8,10(2f13.8,5x)))")  &           ! Spin in loop
 
542
                energy, (cohp_vals(i,is,ic),is=1,nspin)
 
543
        enddo
 
544
        close(tab_u)
 
545
     enddo
 
546
 
 
547
  endif  ! coop
 
548
 
 
549
     !      PDOS output
 
550
     !
 
551
  if (dos) then
 
552
     do ic = 1, ncb
 
553
        open(tab_u,file=trim(mflnm)// "." // trim(tit(ic)) // '.pdos')
 
554
        write(tab_u,"(a1,14x,'ENERGY',4(16x,a1,i1))") '#', ("s",m, m=1,nsp)
 
555
        do i=1, npts_energy
 
556
           energy = min_energy + e_step*(i-1)
 
557
           write(tab_u,"(f20.8,10(2f13.8,5x)))")  &           ! Spin in loop
 
558
                energy, (pdos_vals(i,is,ic),is=1,nspin)
 
559
        enddo
 
560
        close(tab_u)
 
561
     enddo
 
562
  endif
 
563
 
 
564
CONTAINS
 
565
 
 
566
  function delta(x) result(res)
 
567
    real(dp), intent(in) :: x
 
568
    real(dp)             :: res
 
569
 
 
570
    if ((abs(x) > 8*smear)) then
 
571
       res = 0.0_dp
 
572
       RETURN
 
573
    endif
 
574
 
 
575
    res = exp(-(x/smear)**2) / (smear*sqrt(pi))
 
576
 
 
577
  end function delta
 
578
end program mprop
 
579
 
 
580