This routine prints the band energies (and possibly the band derivatives)
This routine is parallel, even if the scaling is very bad since at the moment everything must be written by the root node (we need that the output is sorted). But at least if works independently of the number of processors. I think that a way to write in parallel to the output would help a lot, so that we don't have to send all eigenvalues to the root node.
subroutine geninterp_main()
!! This routine prints the band energies (and possibly the band derivatives)
!!
!! This routine is parallel, even if ***the scaling is very bad*** since at the moment
!! everything must be written by the root node (we need that the output is sorted).
!! But at least if works independently of the number of processors.
!! I think that a way to write in parallel to the output would help a lot,
!! so that we don't have to send all eigenvalues to the root node.
integer :: kpt_unit, outdat_unit, num_kpts, ierr, i, j, k, enidx
character(len=500) :: commentline
character(len=50) :: cdum
integer, dimension(:), allocatable :: kpointidx, localkpointidx
real(kind=dp), dimension(:, :), allocatable :: kpoints, localkpoints
complex(kind=dp), dimension(:, :), allocatable :: HH
complex(kind=dp), dimension(:, :), allocatable :: UU
complex(kind=dp), dimension(:, :, :), allocatable :: delHH
real(kind=dp), dimension(3) :: kpt, frac
real(kind=dp), dimension(:, :, :), allocatable :: localdeleig
real(kind=dp), dimension(:, :, :), allocatable :: globaldeleig
real(kind=dp), dimension(:, :), allocatable :: localeig
real(kind=dp), dimension(:, :), allocatable :: globaleig
logical :: absoluteCoords
character(len=200) :: outdat_filename
integer, dimension(0:num_nodes - 1) :: counts
integer, dimension(0:num_nodes - 1) :: displs
if (on_root .and. (timing_level > 0)) call io_stopwatch('geninterp_main', 1)
if (on_root) then
write (stdout, *)
write (stdout, '(1x,a)') '*---------------------------------------------------------------------------*'
write (stdout, '(1x,a)') '| Generic Band Interpolation routines |'
write (stdout, '(1x,a)') '*---------------------------------------------------------------------------*'
kpt_unit = io_file_unit()
open (unit=kpt_unit, file=trim(seedname)//'_geninterp.kpt', form='formatted', status='old', err=105)
! First line: comment (e.g. creation date, author, ...)
read (kpt_unit, '(A500)', err=106, end=106) commentline
read (kpt_unit, *, err=106, end=106) cdum
if (index(cdum, 'crystal') > 0) then
absoluteCoords = .false.
elseif (index(cdum, 'frac') > 0) then
absoluteCoords = .false.
elseif (index(cdum, 'cart') > 0) then
absoluteCoords = .true.
elseif (index(cdum, 'abs') > 0) then
absoluteCoords = .true.
else
call io_error('Error on second line of file '//trim(seedname)//'_geninterp.kpt: '// &
'unable to recognize keyword')
end if
! Third line: number of following kpoints
read (kpt_unit, *, err=106, end=106) num_kpts
end if
call comms_bcast(num_kpts, 1)
allocate (HH(num_wann, num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating HH in calcTDF')
allocate (UU(num_wann, num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating UU in calcTDF')
if (geninterp_alsofirstder) then
allocate (delHH(num_wann, num_wann, 3), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating delHH in calcTDF')
end if
! I call once the routine to calculate the Hamiltonian in real-space <0n|H|Rm>
call get_HH_R
if (on_root) then
allocate (kpointidx(num_kpts), stat=ierr)
if (ierr /= 0) call io_error('Error allocating kpointidx in geinterp_main.')
allocate (kpoints(3, num_kpts), stat=ierr)
if (ierr /= 0) call io_error('Error allocating kpoints in geinterp_main.')
if (geninterp_single_file) then
allocate (globaleig(num_wann, num_kpts), stat=ierr)
if (ierr /= 0) call io_error('Error allocating globaleig in geinterp_main.')
allocate (globaldeleig(num_wann, 3, num_kpts), stat=ierr)
if (ierr /= 0) call io_error('Error allocating globaldeleig in geinterp_main.')
end if
else
! On the other nodes, I still allocate them with size 1 to avoid
! that some compilers still try to access the memory
allocate (kpointidx(1), stat=ierr)
if (ierr /= 0) call io_error('Error allocating kpointidx in geinterp_main.')
allocate (kpoints(1, 1), stat=ierr)
if (ierr /= 0) call io_error('Error allocating kpoints in geinterp_main.')
if (geninterp_single_file) then
allocate (globaleig(num_wann, 1), stat=ierr)
if (ierr /= 0) call io_error('Error allocating globaleig in geinterp_main.')
allocate (globaldeleig(num_wann, 3, 1), stat=ierr)
if (ierr /= 0) call io_error('Error allocating globaldeleig in geinterp_main.')
end if
end if
! I precalculate how to split on different nodes
call comms_array_split(num_kpts, counts, displs)
allocate (localkpoints(3, max(1, counts(my_node_id))), stat=ierr)
if (ierr /= 0) call io_error('Error allocating localkpoints in geinterp_main.')
allocate (localeig(num_wann, max(1, counts(my_node_id))), stat=ierr)
if (ierr /= 0) call io_error('Error allocating localeig in geinterp_main.')
allocate (localdeleig(num_wann, 3, max(1, counts(my_node_id))), stat=ierr)
if (ierr /= 0) call io_error('Error allocating localdeleig in geinterp_main.')
! On root, I read numpoints_thischunk points
if (on_root) then
! Lines with integer identifier and three coordinates
! (in crystallographic coordinates relative to the reciprocal lattice vectors)
do i = 1, num_kpts
read (kpt_unit, *, err=106, end=106) kpointidx(i), kpt
! Internally, I need the relative (fractional) coordinates in units of the reciprocal-lattice vectors
if (absoluteCoords .eqv. .false.) then
kpoints(:, i) = kpt
else
kpoints(:, i) = 0._dp
! I use the real_lattice (transposed) and a factor of 2pi instead of inverting again recip_lattice
do j = 1, 3
kpoints(j, i) = real_lattice(j, 1)*kpt(1) + real_lattice(j, 2)*kpt(2) + &
real_lattice(j, 3)*kpt(3)
end do
kpoints(:, i) = kpoints(:, i)/(2._dp*pi)
end if
end do
close (kpt_unit)
end if
! Now, I distribute the kpoints; 3* because I send kx, ky, kz
call comms_scatterv(localkpoints, 3*counts(my_node_id), kpoints, 3*counts, 3*displs)
if (.not. geninterp_single_file) then
! Allocate at least one entry, even if we don't use it
allocate (localkpointidx(max(1, counts(my_node_id))), stat=ierr)
if (ierr /= 0) call io_error('Error allocating localkpointidx in geinterp_main.')
call comms_scatterv(localkpointidx, counts(my_node_id), kpointidx, counts, displs)
end if
! I open the output file(s)
if (geninterp_single_file) then
if (on_root) then
outdat_filename = trim(seedname)//'_geninterp.dat'
outdat_unit = io_file_unit()
open (unit=outdat_unit, file=trim(outdat_filename), form='formatted', err=107)
call internal_write_header(outdat_unit, commentline)
end if
else
if (num_nodes > 99999) then
write (outdat_filename, '(a,a,I0,a)') trim(seedname), '_geninterp_', my_node_id, '.dat'
else
write (outdat_filename, '(a,a,I5.5,a)') trim(seedname), '_geninterp_', my_node_id, '.dat'
endif
outdat_unit = io_file_unit()
open (unit=outdat_unit, file=trim(outdat_filename), form='formatted', err=107)
call comms_bcast(commentline, len(commentline))
call internal_write_header(outdat_unit, commentline)
end if
! And now, each node calculates its own k points
do i = 1, counts(my_node_id)
kpt = localkpoints(:, i)
! Here I get the band energies and the velocities (if required)
if (geninterp_alsofirstder) then
call wham_get_eig_deleig(kpt, localeig(:, i), localdeleig(:, :, i), HH, delHH, UU)
else
call pw90common_fourier_R_to_k(kpt, HH_R, HH, 0)
call utility_diagonalize(HH, num_wann, localeig(:, i), UU)
end if
end do
if (geninterp_single_file) then
! Now, I get the results from the different nodes
call comms_gatherv(localeig, num_wann*counts(my_node_id), globaleig, &
num_wann*counts, num_wann*displs)
if (geninterp_alsofirstder) then
call comms_gatherv(localdeleig, 3*num_wann*counts(my_node_id), globaldeleig, &
3*num_wann*counts, 3*num_wann*displs)
end if
! Now the printing, only on root node
if (on_root) then
do i = 1, num_kpts
kpt = kpoints(:, i)
! First calculate the absolute coordinates for printing
frac = 0._dp
do j = 1, 3
frac(j) = recip_lattice(1, j)*kpt(1) + recip_lattice(2, j)*kpt(2) + recip_lattice(3, j)*kpt(3)
end do
! I print each line
if (geninterp_alsofirstder) then
do enidx = 1, num_wann
write (outdat_unit, '(I10,7G18.10)') kpointidx(i), frac, &
globaleig(enidx, i), globaldeleig(enidx, :, i)
end do
else
do enidx = 1, num_wann
write (outdat_unit, '(I10,4G18.10)') kpointidx(i), frac, globaleig(enidx, i)
end do
end if
end do
close (outdat_unit)
end if
else
! Each node simply writes to its own file
do i = 1, counts(my_node_id)
kpt = localkpoints(:, i)
! First calculate the absolute coordinates for printing
frac = 0._dp
do j = 1, 3
frac(j) = recip_lattice(1, j)*kpt(1) + recip_lattice(2, j)*kpt(2) + recip_lattice(3, j)*kpt(3)
end do
! I print each line
if (geninterp_alsofirstder) then
do enidx = 1, num_wann
write (outdat_unit, '(I10,7G18.10)') localkpointidx(i), frac, &
localeig(enidx, i), localdeleig(enidx, :, i)
end do
else
do enidx = 1, num_wann
write (outdat_unit, '(I10,4G18.10)') localkpointidx(i), frac, localeig(enidx, i)
end do
end if
end do
close (outdat_unit)
end if
! All k points processed: Final processing/deallocations
if (on_root) then
write (stdout, '(1x,a)') '| All done. |'
write (stdout, '(1x,a)') '*---------------------------------------------------------------------------*'
end if
if (allocated(kpointidx)) deallocate (kpointidx)
if (allocated(kpoints)) deallocate (kpoints)
if (allocated(localkpoints)) deallocate (localkpoints)
if (allocated(HH)) deallocate (HH)
if (allocated(UU)) deallocate (UU)
if (allocated(delHH)) deallocate (delHH)
if (allocated(localeig)) deallocate (localeig)
if (allocated(localdeleig)) deallocate (localdeleig)
if (allocated(globaleig)) deallocate (globaleig)
if (allocated(globaldeleig)) deallocate (globaldeleig)
if (on_root .and. (timing_level > 0)) call io_stopwatch('geninterp_main', 2)
return
105 call io_error('Error: Problem opening k-point file '//trim(seedname)//'_geninterp.kpt')
106 call io_error('Error: Problem reading k-point file '//trim(seedname)//'_geninterp.kpt')
107 call io_error('Error: Problem opening output file '//trim(outdat_filename))
end subroutine geninterp_main