geninterp.F90 Source File

This File Depends On

sourcefile~~geninterp.f90~~EfferentGraph sourcefile~geninterp.f90 geninterp.F90 sourcefile~constants.f90 constants.F90 sourcefile~constants.f90->sourcefile~geninterp.f90 sourcefile~postw90_common.f90 postw90_common.F90 sourcefile~constants.f90->sourcefile~postw90_common.f90 sourcefile~utility.f90 utility.F90 sourcefile~constants.f90->sourcefile~utility.f90 sourcefile~comms.f90 comms.F90 sourcefile~constants.f90->sourcefile~comms.f90 sourcefile~wan_ham.f90 wan_ham.F90 sourcefile~constants.f90->sourcefile~wan_ham.f90 sourcefile~parameters.f90 parameters.F90 sourcefile~constants.f90->sourcefile~parameters.f90 sourcefile~get_oper.f90 get_oper.F90 sourcefile~constants.f90->sourcefile~get_oper.f90 sourcefile~io.f90 io.F90 sourcefile~constants.f90->sourcefile~io.f90 sourcefile~ws_distance.f90 ws_distance.F90 sourcefile~constants.f90->sourcefile~ws_distance.f90 sourcefile~postw90_common.f90->sourcefile~geninterp.f90 sourcefile~postw90_common.f90->sourcefile~wan_ham.f90 sourcefile~postw90_common.f90->sourcefile~get_oper.f90 sourcefile~utility.f90->sourcefile~geninterp.f90 sourcefile~utility.f90->sourcefile~postw90_common.f90 sourcefile~utility.f90->sourcefile~wan_ham.f90 sourcefile~utility.f90->sourcefile~parameters.f90 sourcefile~utility.f90->sourcefile~ws_distance.f90 sourcefile~comms.f90->sourcefile~geninterp.f90 sourcefile~comms.f90->sourcefile~postw90_common.f90 sourcefile~comms.f90->sourcefile~get_oper.f90 sourcefile~wan_ham.f90->sourcefile~geninterp.f90 sourcefile~parameters.f90->sourcefile~geninterp.f90 sourcefile~parameters.f90->sourcefile~postw90_common.f90 sourcefile~parameters.f90->sourcefile~wan_ham.f90 sourcefile~parameters.f90->sourcefile~get_oper.f90 sourcefile~parameters.f90->sourcefile~ws_distance.f90 sourcefile~get_oper.f90->sourcefile~geninterp.f90 sourcefile~get_oper.f90->sourcefile~wan_ham.f90 sourcefile~io.f90->sourcefile~geninterp.f90 sourcefile~io.f90->sourcefile~postw90_common.f90 sourcefile~io.f90->sourcefile~utility.f90 sourcefile~io.f90->sourcefile~comms.f90 sourcefile~io.f90->sourcefile~parameters.f90 sourcefile~io.f90->sourcefile~get_oper.f90 sourcefile~io.f90->sourcefile~ws_distance.f90 sourcefile~ws_distance.f90->sourcefile~postw90_common.f90
Help

Files Dependent On This One

sourcefile~~geninterp.f90~~AfferentGraph sourcefile~geninterp.f90 geninterp.F90 sourcefile~postw90.f90 postw90.F90 sourcefile~geninterp.f90->sourcefile~postw90.f90
Help

Source Code


Source Code

!-*- mode: F90 -*-!
!------------------------------------------------------------!
! This file is distributed as part of the Wannier90 code and !
! under the terms of the GNU General Public License. See the !
! file `LICENSE' in the root directory of the Wannier90      !
! distribution, or http://www.gnu.org/copyleft/gpl.txt       !
!                                                            !
! The webpage of the Wannier90 code is www.wannier.org       !
!                                                            !
! The Wannier90 code is hosted on GitHub:                    !
!                                                            !
! https://github.com/wannier-developers/wannier90            !
!------------------------------------------------------------!
!                                                            !
!                                                            !
!------------------------------------------------------------!

module w90_geninterp
  !! Generic Interpolation Routine 
  !!
  !! written by Giovanni Pizzi
  !! THEOS, EPFL, Station 12, 1015 Lausanne (Switzerland)  
  !! June, 2012 

  use w90_constants
  use w90_parameters, only    : geninterp_alsofirstder, num_wann, recip_lattice, real_lattice, &
       timing_level, geninterp_single_file
  use w90_io, only            : io_error,stdout,io_stopwatch,io_file_unit,seedname,io_stopwatch
  use w90_get_oper, only      : get_HH_R, HH_R
  use w90_comms
  use w90_utility, only       : utility_diagonalize
  use w90_postw90_common, only : pw90common_fourier_R_to_k
  use w90_wan_ham, only       : wham_get_eig_deleig
  use w90_io, only            : io_date
  implicit none

  private 
  public :: geninterp_main

contains 

  subroutine internal_write_header(outdat_unit,commentline)
  !! Writes a header for the output file(s).

    integer, intent(in) :: outdat_unit
    !! Integer with the output file unit. The file must be already open.
    character(len=*)    :: commentline !! no intent?
    !! String with the comment taken from the output, to be written on the output

    character(len=9)   :: cdate, ctime

    call io_date(cdate,ctime)
    write(outdat_unit, '(A)') "# Written on "//cdate//" at "//ctime ! Date and time    
    ! I rewrite the comment line on the output
    write(outdat_unit, '(A)') "# Input file comment: " // trim(commentline)

    if (geninterp_alsofirstder) then
       write(outdat_unit, '(A)') "#  Kpt_idx  K_x (1/ang)       K_y (1/ang)        K_z (1/ang)       Energy (eV)" // &
            "      EnergyDer_x       EnergyDer_y       EnergyDer_z"
    else
       write(outdat_unit, '(A)') "#  Kpt_idx  K_x (1/ang)       K_y (1/ang)        K_z (1/ang)       Energy (eV)"
    end if
  end subroutine internal_write_header
    
  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
    end if
       
    ! I precalculate how to split on different nodes
    call comms_array_split(num_kpts,counts,displs)
       
    allocate(localkpoints(3,counts(my_node_id)),stat=ierr)
    if (ierr/=0) call io_error('Error allocating localkpoints in geinterp_main.')       
    
    allocate(localeig(num_wann,counts(my_node_id)),stat=ierr)
    if (ierr/=0) call io_error('Error allocating localeig in geinterp_main.')
    allocate(localdeleig(num_wann,3,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(1,1),3*counts(my_node_id),kpoints(1,1),3*counts, 3*displs)
    if  (.not.geninterp_single_file) then
       allocate(localkpointidx(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(:,1,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(1,1),num_wann*counts(my_node_id),globaleig(1,1), &
            num_wann*counts, num_wann*displs)
       
       if (geninterp_alsofirstder) then
          call comms_gatherv(localdeleig(1,1,1),3*num_wann*counts(my_node_id),globaldeleig(1,1,1), &
               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

end module w90_geninterp