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