Octopus
states_elec_io.F90
Go to the documentation of this file.
1!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
2!!
3!! This program is free software; you can redistribute it and/or modify
4!! it under the terms of the GNU General Public License as published by
5!! the Free Software Foundation; either version 2, or (at your option)
6!! any later version.
7!!
8!! This program is distributed in the hope that it will be useful,
9!! but WITHOUT ANY WARRANTY; without even the implied warranty of
10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11!! GNU General Public License for more details.
12!!
13!! You should have received a copy of the GNU General Public License
14!! along with this program; if not, write to the Free Software
15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
16!! 02110-1301, USA.
17!!
18
19#include "global.h"
20
22
24 use accel_oct_m
26 use box_oct_m
27 use comm_oct_m
28 use debug_oct_m
30 use global_oct_m
31 use grid_oct_m
32 use io_oct_m
33 use ions_oct_m
34 use, intrinsic :: iso_fortran_env
36 use math_oct_m
37 use mesh_oct_m
40 use mpi_oct_m
44 use parser_oct_m
45 use phase_oct_m
47 use smear_oct_m
48 use sort_oct_m
49 use space_oct_m
55 use unit_oct_m
57 use utils_oct_m
58 use types_oct_m
60
61 implicit none
62
63 private
64
65 public :: &
70
71contains
72
73 ! ---------------------------------------------------------
75 !
76 subroutine states_elec_write_eigenvalues(nst, st, space, kpoints, error, st_start, compact, iunit, namespace)
77 integer, intent(in) :: nst
78 type(states_elec_t), intent(in) :: st
79 class(space_t), intent(in) :: space
80 type(kpoints_t), intent(in) :: kpoints
81 real(real64), optional, intent(in) :: error(:,:)
82 integer, optional, intent(in) :: st_start
83 logical, optional, intent(in) :: compact
84 integer, optional, intent(in) :: iunit
85 type(namespace_t), optional, intent(in) :: namespace
86
87 integer :: ik, ikk, ist, ns, is, idir, st_start_, iflat, iqn, homo_index, not_printed
88 logical :: print_eigenval
89 real(real64) :: kpoint(space%dim), max_error
90 character(len=120) :: tmp_str(max(space%dim, 3)), cspin
91
92 real(real64), allocatable :: flat_eigenval(:)
93 integer, allocatable :: flat_indices(:, :)
94 integer, parameter :: print_range = 8
95
97
98 st_start_ = 1
99 if (present(st_start)) st_start_ = st_start
100 assert(nst <= st%nst)
101
102
103 if (.not. mpi_grp_is_root(mpi_world)) then
105 return
106 end if
107
108 if (.not. optional_default(compact, .false.)) then
109
110 ns = 1
111 if (st%d%nspin == 2) ns = 2
112
113 message(1) = 'Eigenvalues [' // trim(units_abbrev(units_out%energy)) // ']'
114 call messages_info(1, iunit=iunit, namespace=namespace)
115
116 if (st%d%ispin == spinors) then
117 write(message(1), '(a4,1x,a5,1x,a12,1x,a12,2x,a4,4x,a4,4x,a4)') &
118 '#st',' Spin',' Eigenvalue', 'Occupation ', '<Sx>', '<Sy>', '<Sz>'
119 else
120 write(message(1), '(a4,1x,a5,1x,a12,4x,a12)') &
121 '#st',' Spin',' Eigenvalue', 'Occupation'
122 end if
123 if (present(error)) then
124 write(message(1),'(a,a10)') trim(message(1)), ' Error'
125 end if
126 call messages_info(1, iunit=iunit, namespace=namespace)
127
128 do ik = 1, st%nik, ns
129 if (space%is_periodic()) then
130 ikk = st%d%get_kpoint_index(ik)
131 kpoint(1:space%dim) = kpoints%get_point(ikk, absolute_coordinates = .false.)
132 write(message(1), '(a,i8,a)') '#k =', ikk, ', k = ('
133 do idir = 1, space%dim
134 write(tmp_str(1), '(f10.6)') kpoint(idir)
135 message(1) = trim(message(1))//trim(tmp_str(1))
136 if (idir < space%dim) message(1) = trim(message(1))//','
137 end do
138 message(1) = trim(message(1))//')'
139 call messages_info(1, iunit=iunit, namespace=namespace)
140 end if
141
142 do ist = st_start_, nst
143 do is = 0, ns-1
144 if (is == 0) cspin = 'up'
145 if (is == 1) cspin = 'dn'
146 if (st%d%ispin == unpolarized .or. st%d%ispin == spinors) cspin = '--'
147
148 write(tmp_str(1), '(i4,3x,a2)') ist, trim(cspin)
149 if (st%d%ispin == spinors) then
150 write(tmp_str(2), '(1x,f12.6,5x,f5.2,3x,3f8.4)') &
151 units_from_atomic(units_out%energy, st%eigenval(ist, ik)), st%occ(ist, ik), st%spin(1:3, ist, ik)
152 if (present(error)) write(tmp_str(3), '(a3,es8.1,a1)')' (', error(ist, ik), ')'
153 else
154 write(tmp_str(2), '(1x,f12.6,3x,f12.6)') &
155 units_from_atomic(units_out%energy, st%eigenval(ist, ik+is)), st%occ(ist, ik+is)
156 if (present(error)) write(tmp_str(3), '(a7,es8.1,a1)')' (', error(ist, ik+is), ')'
157 end if
158 if (present(error)) then
159 message(1) = trim(tmp_str(1))//trim(tmp_str(2))//trim(tmp_str(3))
160 else
161 message(1) = trim(tmp_str(1))//trim(tmp_str(2))
162 end if
163 call messages_info(1, iunit=iunit, namespace=namespace)
164 end do
165 end do
166 end do
167
168 else
170 call messages_info(1, iunit=iunit, namespace=namespace)
171
172 safe_allocate(flat_eigenval(1:st%nik*nst))
173 safe_allocate(flat_indices(1:2, 1:st%nik*nst))
174
175 iflat = 1
176 do iqn = 1, st%nik
177 do ist = 1, nst
178
179 flat_eigenval(iflat) = st%eigenval(ist, iqn)
180 flat_indices(1:2, iflat) = (/iqn, ist/)
181
182 iflat = iflat + 1
183 end do
184 end do
185
186 call sort(flat_eigenval, flat_indices(:, :))
187
188 homo_index = st%nik*nst
189 do iflat = 1, st%nik*nst
190 iqn = flat_indices(1, iflat)
191 ist = flat_indices(2, iflat)
192 if (abs(st%occ(ist, iqn)) < 0.1_real64) then
193 homo_index = iflat - 1
194 exit
195 end if
196 end do
197
198 tmp_str(1) = '# State'
199
200 if (space%is_periodic()) tmp_str(1) = trim(tmp_str(1))//' KPoint'
201
202 if (st%d%ispin == spin_polarized) tmp_str(1) = trim(tmp_str(1))//' Spin'
203
204 tmp_str(1) = trim(tmp_str(1))//' Eigenvalue ['// trim(units_abbrev(units_out%energy)) // ']'
205
206 tmp_str(1) = trim(tmp_str(1))//' Occupation'
207
208 if (st%d%ispin == spinors) then
209 tmp_str(1) = trim(tmp_str(1))//' <Sx> <Sy> <Sz>'
210 end if
211
212 if (present(error)) tmp_str(1) = trim(tmp_str(1))//' Error'
213
214 call messages_write(tmp_str(1))
215 call messages_info(iunit=iunit, namespace=namespace)
216
217 not_printed = 0
218 max_error = 0.0_real64
219 do iflat = 1, st%nik*nst
220 iqn = flat_indices(1, iflat)
221 ist = flat_indices(2, iflat)
222 ik = st%d%get_kpoint_index(iqn)
223 is = st%d%get_spin_index(iqn)
224
225 print_eigenval = iflat <= print_range
226 print_eigenval = print_eigenval .or. st%nik*nst - iflat < print_range
227 print_eigenval = print_eigenval .or. abs(iflat - homo_index) <= print_range
228
229 if (print_eigenval) then
230
231 if (not_printed > 0) then
232 call messages_write('')
233 call messages_new_line()
234 call messages_write(' [output of ')
235 call messages_write(not_printed)
236 call messages_write(' eigenvalues skipped')
237 if (present(error)) then
238 call messages_write(': maximum error =')
239 call messages_write(max_error, fmt = '(es7.1)', align_left = .true.)
240 end if
241 call messages_write(']')
242 call messages_new_line()
243 call messages_write('')
244 call messages_info(iunit=iunit, namespace=namespace)
245
246 not_printed = 0
247 max_error = 0.0_real64
248
249 end if
250
251 write(tmp_str(1), '(i7)') ist
252
253 if (space%is_periodic()) then
254 write(tmp_str(1), '(2a,i7)') trim(tmp_str(1)), ' ', ik
255 end if
256
257 if (st%d%ispin == spin_polarized) then
258 if (is == 1) cspin = ' up'
259 if (is == 2) cspin = ' dn'
260 write(tmp_str(1), '(2a,a5)') trim(tmp_str(1)), ' ', cspin
261 end if
262
263
264 if (len(units_abbrev(units_out%energy)) == 1) then
265 write(tmp_str(1), '(2a,f14.6)') trim(tmp_str(1)), ' ', units_from_atomic(units_out%energy, st%eigenval(ist, iqn))
266 else
267 write(tmp_str(1), '(2a,f15.6)') trim(tmp_str(1)), ' ', units_from_atomic(units_out%energy, st%eigenval(ist, iqn))
268 end if
269
270 write(tmp_str(1), '(2a,f11.6)') trim(tmp_str(1)), ' ', st%occ(ist, iqn)
271
272 if (st%d%ispin == spinors) then
273 write(tmp_str(1), '(2a,3f9.4)') trim(tmp_str(1)), ' ', st%spin(1:3, ist, iqn)
274 end if
275
276 if (present(error)) then
277 write(tmp_str(1), '(2a,es8.1,a)') trim(tmp_str(1)), ' (', error(ist, iqn), ')'
278 end if
279
280 call messages_write(tmp_str(1))
281 call messages_info(iunit=iunit, namespace=namespace)
282
283 else
284
285 not_printed = not_printed + 1
286
287 if (present(error)) then
288 max_error = max(max_error, error(ist, iqn))
289 end if
290
291 end if
292
293 end do
294
295 if (nst*st%nik > 1) call print_dos()
296
297 safe_deallocate_a(flat_indices)
298 safe_deallocate_a(flat_eigenval)
299
300 end if
301
302 call smear_write_info(st%smear, namespace, iunit)
303
305
306 contains
307
308 subroutine print_dos()
309
310 integer, parameter :: ndiv = 70, height = 10
311 integer :: histogram(1:ndiv), iev, ien, iline, ife
312 real(real64) :: dhistogram(1:ndiv)
313 character(len=ndiv) :: line
314 real(real64) :: emin, emax, de, maxhist
315
317
318 emin = flat_eigenval(1)
319 emax = flat_eigenval(st%nik*nst)
320 de = (emax - emin)/(ndiv - m_one)
321
322 if (de < m_epsilon) then
324 return
325 end if
326
327 ife = nint((st%smear%e_fermi - emin)/de) + 1
328
329 histogram = 0
330 dhistogram = m_zero
331 do iev = 1, st%nik*nst
332 ien = nint((flat_eigenval(iev) - emin)/de) + 1
333 assert(ien >= 1)
334 assert(ien <= ndiv)
335 dhistogram(ien) = dhistogram(ien) + st%kweights(flat_indices(1, iev))*kpoints%full%npoints
336 end do
337
338 !normalize
339 maxhist = maxval(dhistogram)
340 do ien = 1, ndiv
341 histogram(ien) = nint(dhistogram(ien)*height/maxhist)
342 end do
343
344 call messages_new_line()
345 call messages_write('Density of states:')
346 call messages_new_line()
347 call messages_info(iunit=iunit, namespace=namespace)
348
349 !print histogram
350 do iline = height, 1, -1
351 do ien = 1, ndiv
352 if (histogram(ien) >= iline) then
353 call messages_write('%')
354 else
355 call messages_write('-')
356 end if
357 end do
358 call messages_info(iunit=iunit, namespace=namespace)
359 end do
360
361 line(1:ndiv) = ' '
362 line(ife:ife) = '^'
363 call messages_write(line)
364 call messages_new_line()
365 call messages_info(iunit=iunit, namespace=namespace)
366
368 end subroutine print_dos
369
370 end subroutine states_elec_write_eigenvalues
371
372 ! ---------------------------------------------------------
374 !
375 subroutine states_elec_write_gaps(iunit, st, space)
376 integer, intent(in) :: iunit
377 type(states_elec_t), intent(in) :: st
378 class(space_t), intent(in) :: space
379
380 integer :: ik, ist
381
382 real(real64) :: homo, lumo, egdir, egindir
383 integer :: homok, lumok, egdirk
384
385 push_sub(states_elec_write_gaps)
386
387 if (.not. mpi_grp_is_root(mpi_world) .or. .not. space%is_periodic() &
388 .or. st%smear%method /= smear_semiconductor) then
390 return
391 end if
392
393 !We do not want to compute the bandgap in case there is no extra states
394 if (ceiling(st%qtot/st%smear%el_per_state) == st%nst) then
396 return
397 end if
398
399 homo = -m_huge
400 homok = -1
401 lumo = m_huge
402 lumok = -1
403 egdir = m_huge
404 egdirk = -1
405
406 !TODO: This definition depends on the type of smearing
407 do ik = 1, st%nik
408 if (abs(st%kweights(ik)) < m_epsilon) cycle
409 do ist = 1,st%nst-1
410 if (st%occ(ist, ik) > m_epsilon .and. st%eigenval(ist, ik) > homo) then
411 homo = st%eigenval(ist, ik)
412 homok = ik
413 end if
414
415 if (st%occ(ist+1, ik) <= m_epsilon .and. st%eigenval(ist+1, ik) < lumo) then
416 lumo = st%eigenval(ist+1, ik)
417 lumok = ik
418 end if
419
420 if (st%occ(ist, ik) > m_epsilon .and. st%occ(ist+1, ik) <= m_epsilon &
421 .and. (st%eigenval(ist+1, ik) - st%eigenval(ist, ik)) < egdir) then
422 egdir = (st%eigenval(ist+1, ik) - st%eigenval(ist, ik))
423 egdirk = ik
424 end if
425 end do
426 end do
427
428 egindir = lumo-homo
429
430 if (lumok == -1 .or. egdir <= m_epsilon) then
431 write(message(1),'(a)') 'The system seems to have no gap.'
432 call messages_info(1, iunit)
433 else
434 write(message(1),'(a,i5,a,f7.4,a)') 'Direct gap at ik=', egdirk, ' of ', &
435 units_from_atomic(units_out%energy, egdir), ' ' // trim(units_abbrev(units_out%energy))
436 write(message(2),'(a,i5,a,i5,a,f7.4,a)') 'Indirect gap between ik=', homok, ' and ik=', lumok, &
437 ' of ', units_from_atomic(units_out%energy, egindir), ' ' // trim(units_abbrev(units_out%energy))
438 call messages_info(2, iunit)
439 end if
440
441
443 end subroutine states_elec_write_gaps
444
445
446 ! ---------------------------------------------------------
450 !
451 subroutine states_elec_write_tpa(dir, namespace, space, mesh, st)
452 character(len=*), intent(in) :: dir
453 type(namespace_t), intent(in) :: namespace
454 class(space_t), intent(in) :: space
455 class(mesh_t), intent(in) :: mesh
456 type(states_elec_t), intent(in) :: st
457
458 type(block_t) :: blk
459 integer :: ncols, icoord, ist, ik, tpa_initialst, tpa_initialk
460 integer :: iunit
461
462 real(real64), allocatable :: ff(:)
463 complex(real64), allocatable :: cff(:)
464 real(real64), allocatable :: osc(:)
465 real(real64) :: transition_energy, osc_strength, dsf
466
467 logical :: use_qvector
468 real(real64), allocatable :: qvector(:), psi_initial(:, :), psi_ist(:, :)
469
470 push_sub(states_elec_write_tpa)
471
472 use_qvector = .false.
473
474 ! find the orbital with half-occupation
475 tpa_initialst = -1
476 do ist = 1, st%nst
477 do ik = 1, st%nik
478 if (abs(st%occ(ist,ik)-m_half) < m_min_occ) then
479 tpa_initialst = ist
480 tpa_initialk = ik
481 end if
482 end do
483 end do
484
485 ! make sure that half-occupancy was found
486 if (tpa_initialst == -1) then
487 if (mpi_grp_is_root(mpi_world)) then
488
489 call messages_write('No orbital with half-occupancy found. TPA output is not written.')
490 call messages_warning(namespace=namespace)
491
492 end if
493 pop_sub(states_elec_write_tpa)
494 return
495 end if
496
497 !%Variable MomentumTransfer
498 !%Type block
499 !%Section Output
500 !%Description
501 !% Momentum-transfer vector <math>\vec{q}</math> to be used when calculating matrix elements
502 !% <math>\left< f \left| e^{i \vec{q} \cdot \vec{r}} \right| i \right></math>.
503 !% This enables the calculation of the dynamical structure factor,
504 !% which is closely related to generalized oscillator strengths.
505 !% If the vector is not given, but TPA output is requested (<tt>Output = TPA</tt>),
506 !% only the oscillator strengths are written in the output file.
507 !% For example, to use <math>\vec{q}</math> = (0.1, 0.2, 0.3), set
508 !%
509 !% <tt>%MomentumTransfer
510 !% <br>&nbsp;&nbsp; 0.1 | 0.2 | 0.3
511 !% <br>%</tt>
512 !%End
513 if (parse_block(namespace, 'MomentumTransfer', blk) == 0) then
514
515 ! check if input makes sense
516 ncols = parse_block_cols(blk, 0)
517
518 if (ncols /= space%dim) then ! wrong size
519
520 if (mpi_grp_is_root(mpi_world)) then
521 call messages_write('Inconsistent size of momentum-transfer vector. It will not be used in the TPA calculation.')
522 call messages_warning(namespace=namespace)
523 end if
524
525 else ! correct size
526
527 use_qvector = .true.
528 safe_allocate(qvector(1:space%dim))
529
530 do icoord = 1, space%dim !for x,y,z
531 call parse_block_float(blk, 0, icoord-1, qvector(icoord))
532 qvector(icoord) = units_to_atomic(unit_one / units_inp%length, qvector(icoord))
533 end do
534
535 end if
536
537 end if
538
539 ! calculate the matrix elements
540
541 safe_allocate(ff(1:mesh%np))
542 if (use_qvector) then
543 safe_allocate(cff(1:mesh%np))
544 end if
545 safe_allocate(osc(1:space%dim))
546
547 ! root writes output to file
548
549 if (mpi_grp_is_root(mpi_world)) then
550
551 iunit = io_open(trim(dir)//'/'//trim('tpa_xas'), namespace, action='write')
552
553 ! header
554 if (use_qvector) then
555 write (message(1),'(a1,a30,3(es14.5,1x),a1)') '#', ' momentum-transfer vector : (', &
556 (units_from_atomic(unit_one / units_out%length, qvector(icoord)), icoord=1, space%dim),')'
557 select case (space%dim)
558 case (1)
559 write(message(2), '(a1,4(a15,1x))') '#', 'E' , '<x>', '<f>', 'S(q,omega)'
560 case (2)
561 write(message(2), '(a1,5(a15,1x))') '#', 'E' , '<x>', '<y>', '<f>', 'S(q,omega)'
562 case (3)
563 write(message(2), '(a1,6(a15,1x))') '#', 'E' , '<x>', '<y>', '<z>', '<f>', 'S(q,omega)'
564 end select
565 call messages_info(2, iunit, namespace=namespace)
566 else
567 select case (space%dim)
568 case (1)
569 write(message(1), '(a1,3(a15,1x))') '#', 'E' , '<x>', '<f>'
570 case (2)
571 write(message(1), '(a1,4(a15,1x))') '#', 'E' , '<x>', '<y>', '<f>'
572 case (3)
573 write(message(1), '(a1,5(a15,1x))') '#', 'E' , '<x>', '<y>', '<z>', '<f>'
574 end select
575 call messages_info(1, iunit, namespace=namespace)
576 end if
577
578 end if
579
580 safe_allocate(psi_initial(1:mesh%np, 1:st%d%dim))
581 safe_allocate(psi_ist(1:mesh%np, 1:st%d%dim))
582
583 call states_elec_get_state(st, mesh, tpa_initialst, tpa_initialk, psi_initial)
584
585 ! loop through every state
586 do ist = 1, st%nst
587
588 call states_elec_get_state(st, mesh, ist, tpa_initialk, psi_ist)
589
590 ! final states are the unoccupied ones
591 if (abs(st%occ(ist,tpa_initialk)) < m_min_occ) then
592
593 osc_strength = m_zero
594 transition_energy = st%eigenval(ist, tpa_initialk) - st%eigenval(tpa_initialst, tpa_initialk)
595
596 ! dipole matrix elements <f|x|i> etc. -> oscillator strengths
597 do icoord = 1, space%dim ! for x,y,z
598
599 ff(1:mesh%np) = psi_initial(1:mesh%np, 1)*mesh%x(1:mesh%np, icoord)* &
600 psi_ist(1:mesh%np, 1)
601 osc(icoord) = dmf_integrate(mesh, ff)
602 osc_strength = osc_strength + m_two/real(space%dim, real64) *transition_energy*abs(osc(icoord))**2
603
604 end do
605
606 ! matrix elements <f|exp(iq.r)|i> -> dynamic structure factor
607 if (use_qvector) then
608
609 cff(1:mesh%np) = psi_initial(1:mesh%np, 1)*psi_ist(1:mesh%np, 1)
610
611 do icoord = 1, space%dim ! for x,y,z
612 cff(1:mesh%np) = cff(1:mesh%np)*exp(m_zi*mesh%x(1:mesh%np, icoord)*qvector(icoord))
613 end do
614
615 dsf = abs(zmf_integrate(mesh, cff))**2
616 end if
617
618 ! write oscillator strengths (+ dynamic structure factor if qvector is given) into file
619 if (mpi_grp_is_root(mpi_world)) then
620
621 if (use_qvector) then
622 write(message(1), '(1x,6(es15.8,1x))') units_from_atomic(units_out%energy, transition_energy), osc(:), osc_strength, &
624 else
625 write(message(1), '(1x,6(es15.8,1x))') units_from_atomic(units_out%energy, transition_energy), osc(:), osc_strength
626 end if
627
628 call messages_info(1, iunit, namespace=namespace)
629
630 end if
631
632 end if
633
634 end do
635
636 ! finally close the file
637 if (mpi_grp_is_root(mpi_world)) then
638 call io_close(iunit)
639 end if
640
641 safe_deallocate_a(psi_initial)
642 safe_deallocate_a(psi_ist)
643
644 safe_deallocate_a(ff)
645 if (use_qvector) then
646 safe_deallocate_a(cff)
647 end if
648 safe_deallocate_a(osc)
649 if (use_qvector) then
650 safe_deallocate_a(qvector)
651 end if
652
653 pop_sub(states_elec_write_tpa)
654
655 end subroutine states_elec_write_tpa
656
657
658 ! ---------------------------------------------------------
660 !
661 subroutine states_elec_write_bandstructure(dir, namespace, nst, st, ions, mesh, kpoints, &
662 phase, vec_pot, vec_pot_var)
663 character(len=*), intent(in) :: dir
664 type(namespace_t), intent(in) :: namespace
665 integer, intent(in) :: nst
666 type(states_elec_t), target, intent(in) :: st
667 type(ions_t), target, intent(in) :: ions
668 class(mesh_t), intent(in) :: mesh
669 type(kpoints_t), intent(in) :: kpoints
670 type(phase_t), intent(in) :: phase
671 real(real64), optional, allocatable, intent(in) :: vec_pot(:)
672 ! !! dimension (1:mesh\%box\%dim)
673 real(real64), optional, allocatable, intent(in) :: vec_pot_var(:, :)
674 ! !! dimensions (1:mesh\%box\%dim, 1:os\%sphere\%ns)
675
676 integer :: idir, ist, ik, ns, is,npath, ib, ind
677 integer, allocatable :: iunit(:)
678 real(real64) :: red_kpoint(mesh%box%dim)
679 character(len=80) :: filename
680
681 logical :: projection
682 integer :: ii, ll, mm, nn, work, norb, work2
683 integer :: ia, iorb, idim, maxnorb
684 real(real64), allocatable :: ddot(:,:,:)
685 complex(real64), allocatable :: zdot(:,:,:)
686 real(real64), allocatable :: weight(:,:,:,:,:)
687 type(orbitalset_t) :: os
688 type(wfs_elec_t), pointer :: epsib
689
690
692
693 ! shortcuts
694 ns = 1
695 if (st%d%nspin == 2) ns = 2
696
697 !%Variable BandStructureComputeProjections
698 !%Type logical
699 !%Default false
700 !%Section Output
701 !%Description
702 !% Determines if projections of wavefunctions on the atomic orbitals
703 !% are computed or not for obtaining the orbital resolved band-structure.
704 !%End
705 call parse_variable(namespace, 'BandStructureComputeProjections', .false., projection)
706
707
708 if (mpi_grp_is_root(mpi_world)) then
709 safe_allocate(iunit(0:ns-1))
710
711 do is = 0, ns-1
712 if (ns > 1) then
713 write(filename, '(a,i1.1,a)') 'bandstructure-sp', is+1
714 else
715 write(filename, '(a)') 'bandstructure'
716 end if
717 iunit(is) = io_open(trim(dir)//'/'//trim(filename), namespace, action='write')
718
719 ! write header
720 write(iunit(is),'(a)',advance='no') '# coord. '
721 do idir = 1, mesh%box%dim
722 write(iunit(is),'(3a)',advance='no') 'k', index2axis(idir), ' '
723 end do
724 if (.not. projection) then
725 write(iunit(is),'(a,i6,3a)') '(red. coord.), bands:', nst, ' [', trim(units_abbrev(units_out%energy)), ']'
726 else
727 write(iunit(is),'(a,i6,3a)',advance='no') '(red. coord.), bands:', nst, ' [', trim(units_abbrev(units_out%energy)), '] '
728 do ia = 1, ions%natoms
729 work = orbitalset_utils_count(ions%atom(ia)%species)
730 do norb = 1, work
731 work2 = orbitalset_utils_count(ions%atom(ia)%species, norb)
732 write(iunit(is),'(a, i3.3,a,i1.1,a)',advance='no') 'w(at=',ia,',os=',norb,') '
733 end do
734 end do
735 write(iunit(is),'(a)') ''
736 end if
737 end do
738 end if
739
740 npath = kpoints%nkpt_in_path()*ns
741
742
743 !We need to compute the projections of each wavefunctions on the localized basis
744 if (projection) then
745
746 maxnorb = 0
747 do ia = 1, ions%natoms
748 maxnorb = max(maxnorb, orbitalset_utils_count(ions%atom(ia)%species))
749 end do
750
751 safe_allocate(weight(1:st%nik,1:st%nst, 1:(2*max_l+1), 1:maxnorb, 1:ions%natoms))
752 weight(1:st%nik,1:st%nst, 1:(2*max_l+1), 1:maxnorb, 1:ions%natoms) = m_zero
753
754 do ia = 1, ions%natoms
755
756 !We first count how many orbital set we have
757 work = orbitalset_utils_count(ions%atom(ia)%species)
758
759 !We loop over the orbital sets of the atom ia
760 do norb = 1, work
761 call orbitalset_init(os)
762
763 !We count the orbitals
764 work2 = 0
765 do iorb = 1, ions%atom(ia)%species%get_niwfs()
766 call ions%atom(ia)%species%get_iwf_ilm(iorb, 1, ii, ll, mm)
767 call ions%atom(ia)%species%get_iwf_n(iorb, 1, nn)
768 if (ii == norb) then
769 os%ll = ll
770 os%nn = nn
771 os%ii = ii
772 os%radius = atomic_orbital_get_radius(ions%atom(ia)%species, mesh, iorb, 1, &
773 option__aotruncation__ao_full, 0.01_real64)
774 work2 = work2 + 1
775 end if
776 end do
777 os%norbs = work2
778 os%ndim = 1
779 os%use_submesh = .false.
780 os%spec => ions%atom(ia)%species
781
782 do iorb = 1, os%norbs
783 ! We obtain the orbital
784 if (states_are_real(st)) then
785 call dget_atomic_orbital(namespace, ions%space, ions%latt, ions%pos(:,ia), &
786 ions%atom(ia)%species, mesh, os%sphere, os%ii, os%ll, os%jj, &
787 os, iorb, os%radius, os%ndim, use_mesh=.not.os%use_submesh, &
788 normalize = .true.)
789 else
790 call zget_atomic_orbital(namespace, ions%space, ions%latt, ions%pos(:,ia), &
791 ions%atom(ia)%species, mesh, os%sphere, os%ii, os%ll, os%jj, &
792 os, iorb, os%radius, os%ndim, &
793 use_mesh = .not. phase%is_allocated() .and. .not. os%use_submesh, &
794 normalize = .true.)
795 end if
796 end do !iorb
797
798 if (phase%is_allocated()) then
799 ! In case of complex wavefunction, we allocate the array for the phase correction
800 safe_allocate(os%phase(1:os%sphere%np, st%d%kpt%start:st%d%kpt%end))
801 os%phase(:,:) = m_zero
802 if (.not. os%use_submesh) then
803 safe_allocate(os%eorb_mesh(1:mesh%np, 1:os%norbs, 1:os%ndim, st%d%kpt%start:st%d%kpt%end))
804 os%eorb_mesh(:,:,:,:) = m_zero
805 else
806 safe_allocate(os%eorb_submesh(1:os%sphere%np, 1:os%ndim, 1:os%norbs, st%d%kpt%start:st%d%kpt%end))
807 os%eorb_submesh(:,:,:,:) = m_zero
808 end if
809
810 if (accel_is_enabled() .and. st%d%dim == 1) then
811 os%ldorbs = max(pad_pow2(os%sphere%np), 1)
812 if(.not. os%use_submesh) os%ldorbs = max(pad_pow2(os%sphere%mesh%np), 1)
813
814 safe_allocate(os%buff_eorb(st%d%kpt%start:st%d%kpt%end))
815 do ik= st%d%kpt%start, st%d%kpt%end
816 call accel_create_buffer(os%buff_eorb(ik), accel_mem_read_only, type_cmplx, os%ldorbs*os%norbs)
817 end do
818 end if
819
820 call orbitalset_update_phase(os, mesh%box%dim, st%d%kpt, kpoints, (st%d%ispin==spin_polarized), &
821 vec_pot, vec_pot_var)
822 end if
823
824 if (states_are_real(st)) then
825 safe_allocate(ddot(1:st%d%dim,1:os%norbs, 1:st%block_size))
826 else
827 safe_allocate(zdot(1:st%d%dim,1:os%norbs, 1:st%block_size))
828 end if
829
830 do ik = st%d%kpt%start, st%d%kpt%end
831 if (ik < st%nik - npath + 1) cycle ! We only want points inside the k-point path
832
833 do ib = st%group%block_start, st%group%block_end
834
835 if (phase%is_allocated()) then
836 safe_allocate(epsib)
837 call st%group%psib(ib, ik)%copy_to(epsib)
838 call phase%apply_to(mesh, mesh%np, .false., epsib, src = st%group%psib(ib, ik))
839 else
840 epsib => st%group%psib(ib, ik)
841 end if
842
843 if (states_are_real(st)) then
844 call dorbitalset_get_coeff_batch(os, st%d%dim, epsib, ddot(:, :, :))
845 do ist = 1, st%group%psib(ib, ik)%nst
846 ind = st%group%psib(ib, ik)%ist(ist)
847 do iorb = 1, os%norbs
848 do idim = 1, st%d%dim
849 weight(ist, ind, iorb, norb, ia) = weight(ik, ind, iorb, norb, ia) + abs(ddot(idim, iorb, ist))**2
850 end do
851 end do
852 end do
853 else
854 call zorbitalset_get_coeff_batch(os, st%d%dim, epsib, zdot(:, :, :))
855 do ist = 1, st%group%psib(ib, ik)%nst
856 ind = st%group%psib(ib, ik)%ist(ist)
857 do iorb = 1, os%norbs
858 do idim = 1, st%d%dim
859 weight(ik, ind, iorb, norb, ia) = weight(ik, ind, iorb, norb, ia) + abs(zdot(idim, iorb, ist))**2
860 end do
861 end do
862 end do
863 end if
864
865 if (phase%is_allocated()) then
866 call epsib%end(copy=.false.)
867 safe_deallocate_p(epsib)
868 end if
869
870 end do
871 end do
872
873 safe_deallocate_a(ddot)
874 safe_deallocate_a(zdot)
875
876 call orbitalset_end(os)
877 end do !norb
878
879 if (st%parallel_in_states .or. st%d%kpt%parallel) then
880 call comm_allreduce(st%st_kpt_mpi_grp, weight(:,:,:,:,ia))
881 end if
882 end do !ia
883
884 end if !projection
885
886 if (mpi_grp_is_root(mpi_world)) then
887 ! output bands
888 do ik = st%nik-npath+1, st%nik, ns
889 do is = 0, ns - 1
890 red_kpoint(1:mesh%box%dim) = kpoints%get_point(st%d%get_kpoint_index(ik + is), absolute_coordinates=.false.)
891 write(iunit(is),'(1x)',advance='no')
892 if (st%nik > npath) then
893 write(iunit(is),'(f14.8)',advance='no') kpoints_get_path_coord(kpoints, &
894 st%d%get_kpoint_index(ik + is)-st%d%get_kpoint_index(st%nik -npath))
895 else
896 write(iunit(is),'(f14.8)',advance='no') kpoints_get_path_coord(kpoints, &
897 st%d%get_kpoint_index(ik + is))
898 end if
899 do idir = 1, mesh%box%dim
900 write(iunit(is),'(f14.8)',advance='no') red_kpoint(idir)
901 end do
902 do ist = 1, nst
903 write(iunit(is),'(1x,f14.8)',advance='no') units_from_atomic(units_out%energy, st%eigenval(ist, ik + is))
904 end do
905 if (projection) then
906 do ia = 1, ions%natoms
907 work = orbitalset_utils_count(ions%atom(ia)%species)
908 do norb = 1, work
909 work2 = orbitalset_utils_count(ions%atom(ia)%species, norb)
910 do iorb = 1, work2
911 do ist = 1, nst
912 write(iunit(is),'(es15.8)',advance='no') weight(ik+is,ist,iorb,norb,ia)
913 end do
914 end do
915 end do
916 end do
917 end if
918 end do
919
920 do is = 0, ns-1
921 write(iunit(is), '(a)')
922 end do
923 end do
924
925 do is = 0, ns-1
926 call io_close(iunit(is))
927 end do
928 end if
929
930 if (projection) then
931 safe_deallocate_a(weight)
932 end if
933
934 safe_deallocate_a(iunit)
935
936
938
940
941end module states_elec_io_oct_m
942
943
944!! Local Variables:
945!! mode: f90
946!! coding: utf-8
947!! End:
This is the common interface to a sorting routine. It performs the shell algorithm,...
Definition: sort.F90:149
double exp(double __x) __attribute__((__nothrow__
pure logical function, public accel_is_enabled()
Definition: accel.F90:402
integer, parameter, public accel_mem_read_only
Definition: accel.F90:183
integer, parameter, public max_l
subroutine, public zget_atomic_orbital(namespace, space, latt, pos, species, mesh, sm, ii, ll, jj, os, orbind, radius, d_dim, use_mesh, normalize, index_shift)
This routine returns the atomic orbital basis – provided by the pseudopotential structure in geo.
real(real64) function, public atomic_orbital_get_radius(species, mesh, iorb, ispin, truncation, threshold)
subroutine, public dget_atomic_orbital(namespace, space, latt, pos, species, mesh, sm, ii, ll, jj, os, orbind, radius, d_dim, use_mesh, normalize, index_shift)
This routine returns the atomic orbital basis – provided by the pseudopotential structure in geo.
integer, parameter, public unpolarized
Parameters...
integer, parameter, public spinors
integer, parameter, public spin_polarized
real(real64), parameter, public m_two
Definition: global.F90:190
real(real64), parameter, public m_huge
Definition: global.F90:206
real(real64), parameter, public m_zero
Definition: global.F90:188
complex(real64), parameter, public m_zi
Definition: global.F90:202
real(real64), parameter, public m_epsilon
Definition: global.F90:204
real(real64), parameter, public m_half
Definition: global.F90:194
real(real64), parameter, public m_one
Definition: global.F90:189
real(real64), parameter, public m_min_occ
Minimal occupation that is considered to be non-zero.
Definition: global.F90:211
This module implements the underlying real-space grid.
Definition: grid.F90:117
Definition: io.F90:114
subroutine, public io_close(iunit, grp)
Definition: io.F90:418
integer function, public io_open(file, namespace, action, status, form, position, die, recl, grp)
Definition: io.F90:352
real(real64) pure function, public kpoints_get_path_coord(this, ind)
Definition: kpoints.F90:1656
This module is intended to contain "only mathematical" functions and procedures.
Definition: math.F90:115
integer pure function, public pad_pow2(size)
create array size, which is padded to powers of 2
Definition: math.F90:889
This module defines various routines, operating on mesh functions.
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:118
subroutine, public messages_warning(no_lines, all_nodes, namespace)
Definition: messages.F90:537
subroutine, public messages_new_line()
Definition: messages.F90:1134
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:160
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:616
logical function mpi_grp_is_root(grp)
Is the current MPI process of grpcomm, root.
Definition: mpi.F90:434
type(mpi_grp_t), public mpi_world
Definition: mpi.F90:270
subroutine, public orbitalset_init(this)
Definition: orbitalset.F90:203
subroutine, public orbitalset_end(this)
Definition: orbitalset.F90:229
subroutine, public orbitalset_update_phase(os, dim, kpt, kpoints, spin_polarized, vec_pot, vec_pot_var, kpt_max)
Build the phase correction to the global phase in case the orbital crosses the border of the simulato...
Definition: orbitalset.F90:279
subroutine, public dorbitalset_get_coeff_batch(os, ndim, psib, dot, reduce)
Definition: orbitalset.F90:595
subroutine, public zorbitalset_get_coeff_batch(os, ndim, psib, dot, reduce)
integer function, public orbitalset_utils_count(species, iselect)
Count the number of orbital sets we have for a given atom.
integer function, public parse_block(namespace, name, blk, check_varinfo_)
Definition: parser.F90:618
subroutine, public smear_write_info(this, namespace, iunit)
Definition: smear.F90:840
integer, parameter, public smear_semiconductor
Definition: smear.F90:171
This module is intended to contain "only mathematical" functions and procedures.
Definition: sort.F90:117
pure logical function, public states_are_real(st)
This module handles spin dimensions of the states and the k-point distribution.
This module defines routines to write information about states.
subroutine, public states_elec_write_tpa(dir, namespace, space, mesh, st)
calculate the transition potential and write to a file
subroutine, public states_elec_write_eigenvalues(nst, st, space, kpoints, error, st_start, compact, iunit, namespace)
write the eigenvalues for some states to a file.
subroutine, public states_elec_write_gaps(iunit, st, space)
calculate gaps and write to a file.
subroutine, public states_elec_write_bandstructure(dir, namespace, nst, st, ions, mesh, kpoints, phase, vec_pot, vec_pot_var)
calculate and write the bandstructure
type(type_t), public type_cmplx
Definition: types.F90:134
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
Definition: unit.F90:132
character(len=20) pure function, public units_abbrev(this)
Definition: unit.F90:223
This module defines the unit system, used for input and output.
type(unit_system_t), public units_out
type(unit_system_t), public units_inp
the units systems for reading and writing
type(unit_t), public unit_one
some special units required for particular quantities
This module is intended to contain simple general-purpose utility functions and procedures.
Definition: utils.F90:118
character pure function, public index2axis(idir)
Definition: utils.F90:202
subroutine print_dos()
Describes mesh distribution to nodes.
Definition: mesh.F90:186
A container for the phase.
Definition: phase.F90:177
The states_elec_t class contains all electronic wave functions.
batches of electronic states
Definition: wfs_elec.F90:139
int true(void)