floquet.F90

Go to the documentation of this file.

!! Copyright (C) 2016 H. Huebener
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
program oct_floquet
  use blas_oct_m
  use calc_mode_par_oct_m
  use comm_oct_m
  use debug_oct_m
  use density_oct_m
  use fft_oct_m
  use gauge_field_oct_m
  use global_oct_m
  use grid_oct_m
  use hamiltonian_elec_oct_m
  use io_oct_m
  use lalg_adv_oct_m
  use mesh_oct_m
  use messages_oct_m
  use mpi_oct_m
  use multicomm_oct_m
  use namespace_oct_m
  use parser_oct_m
  use profiling_oct_m
  use restart_oct_m
  use space_oct_m
  use states_elec_oct_m
  use states_elec_restart_oct_m
  use electrons_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use utils_oct_m
  use v_ks_oct_m
  use xc_oct_m
 
  implicit none
 
  integer :: ierr
 
  type(electrons_t), pointer :: sys
  type(states_elec_t) :: st
  complex(real64), allocatable :: hmss(:,:), psi(:,:,:), hpsi(:,:,:), temp_state1(:,:)
  complex(real64), allocatable :: HFloquet(:,:,:), HFloq_eff(:,:), temp(:,:)
  real(real64), allocatable :: eigenval(:), bands(:,:)
  character(len=80) :: filename
  integer :: it, nT, ik, ist, in, im, file, idim, nik, ik_count
  integer :: Forder, Fdim, m0, n0, n1, nst, ii, jj, lim_nst
  real(real64) :: dt, Tcycle,omega
  logical :: downfolding = .false.
  class(mesh_t), pointer :: mesh
  type(restart_t) :: restart
 
  call global_init()
 
  call parser_init()
 
  call messages_init()
 
  call io_init()
  call profiling_init(global_namespace)
 
  call print_header()
  call messages_print_with_emphasis(msg="Non-interacting Floquet", namespace=global_namespace)
  call messages_print_with_emphasis(namespace=global_namespace)
 
  call messages_experimental("oct-floquet utility")
  call fft_all_init(global_namespace)
  call unit_system_init(global_namespace)
  call restart_module_init(global_namespace)
 
  call calc_mode_par%set_parallelization(p_strategy_states, default = .false.)
  sys => electrons_t(global_namespace)
  call sys%init_parallelization(mpi_world)
  ! make shortcut copies
  st = sys%st
 
  ! generate the full hamiltonian following the sequence in td_init
  call hamiltonian_elec_epot_generate(sys%hm, global_namespace, sys%space, sys%gr, sys%ions, &
    sys%ext_partners, st, time=m_zero)
  call sys%hm%update(sys%gr, global_namespace, sys%space, sys%ext_partners, time = m_zero)
 
  call states_elec_allocate_wfns(st, sys%gr)
 
  call restart_init(restart, global_namespace, restart_gs, restart_type_load, sys%mc, ierr, mesh=sys%gr, exact=.true.)
  if (ierr == 0) then
    call states_elec_load(restart, global_namespace, sys%space, st, sys%gr, sys%kpoints, ierr, label = ": gs")
  end if
  if (ierr /= 0) then
    message(1) = 'Unable to read ground-state wavefunctions.'
    call messages_fatal(1)
  end if
 
  call density_calc(st, sys%gr, st%rho)
  call v_ks_calc(sys%ks, global_namespace, sys%space, sys%hm, st, sys%ions, sys%ext_partners, &
    calc_eigenval=.true., time = m_zero)
  call sys%hm%update(sys%gr, global_namespace, sys%space, sys%ext_partners, time = m_zero)
 
  call floquet_init()
 
  call floquet_solve_non_interacting()
 
 
  ! wait for all processors to finish
  call mpi_world%barrier()
 
  call fft_all_end()
  safe_deallocate_p(sys)
  call profiling_end(global_namespace)
  call io_end()
  call print_date("Calculation ended on ")
  call messages_end()
 
  call parser_end()
  call global_end()
 
contains
 
  !------------------------------------------------
  subroutine floquet_init()
 
    push_sub(floquet_init)
 
    !for now no domain distribution allowed
    assert(sys%gr%np == sys%gr%np_global)
 
    ! variables documented in td/td_write.F90
    call parse_variable(global_namespace, 'TDFloquetFrequency', m_zero, omega, units_inp%energy)
    call messages_print_var_value('Frequency used for Floquet analysis', omega)
    if (abs(omega) <= m_epsilon) then
      message(1) = "Please give a non-zero value for TDFloquetFrequency"
      call messages_fatal(1)
    end if
 
    ! get time of one cycle
    tcycle = m_two * m_pi / omega
 
    call parse_variable(global_namespace, 'TDFloquetSample',20 ,nt)
    call messages_print_var_value('Number of Floquet time-sampling points', nt)
    dt = tcycle/real(nt, real64)
 
    call parse_variable(global_namespace, 'TDFloquetDimension',-1,forder)
    if (forder.ge.0) then
      call messages_print_var_value('Order of multiphoton Floquet-Hamiltonian', forder)
      !Dimension of multiphoton Floquet-Hamiltonian
      fdim = 2 * forder + 1
    else
      message(1) = 'Floquet-Hamiltonian is downfolded'
      call messages_info(1)
      downfolding = .true.
      forder = 1
      fdim = 3
    end if
 
    dt = tcycle/real(nt, real64)
 
    pop_sub(floquet_init)
 
  end subroutine floquet_init
 
  !---------------------------------------------------
  subroutine floquet_solve_non_interacting()
    type(states_elec_t) :: hm_st
 
    push_sub(floquet_solve_non_interacting)
 
    mesh => sys%gr
    nst = st%nst
 
    safe_allocate(hmss(1:nst,1:nst))
    safe_allocate( psi(1:nst,1:st%d%dim,1:mesh%np))
    safe_allocate(hpsi(1:nst,1:st%d%dim,1:mesh%np))
    safe_allocate(temp_state1(1:mesh%np,1:st%d%dim))
 
    ! this is used to initialize the local state object
    call states_elec_copy(hm_st, st)
 
    ! we are only interested for k-point with zero weight
    nik = sys%kpoints%nik_skip
 
    ! multiphoton Floquet Hamiltonian, layout:
    !     (H_{-n,-m} ...  H_{-n,0} ...  H_{-n,m})
    !     (    .      .      .      .      .    )
    ! H = (H_{0,-m}  ...  H_{0,0}  ...  H_{0,m} )
    !     (    .      .      .      .      .    )
    !     (H_{n,-m}  ...  H_{n,0}  ...  H_{n,m} )
    safe_allocate(hfloquet(1:nik,1:nst*fdim, 1:nst*fdim))
    hfloquet(1:nik,1:nst*fdim, 1:nst*fdim) = m_zero
 
    ! perform time-integral over one cycle
    do it = 1, nt
      ! get non-interacting Hamiltonian at time (offset by one cycle to allow for ramp)
      call sys%hm%update(sys%gr, global_namespace, sys%space, sys%ext_partners, &
        time=tcycle+it*dt)
      ! get hpsi
      call zhamiltonian_elec_apply_all(sys%hm, global_namespace, sys%gr, st, hm_st)
 
      ! project Hamiltonian into grounstates for zero weight k-points
      ik_count = 0
 
      do ik = sys%kpoints%reduced%npoints-nik+1, sys%kpoints%reduced%npoints
        ik_count = ik_count + 1
 
        psi(1:nst, 1:st%d%dim, 1:mesh%np)= m_zero
        hpsi(1:nst, 1:st%d%dim, 1:mesh%np)= m_zero
 
        do ist = st%st_start, st%st_end
          if (state_kpt_is_local(st, ist, ik)) then
            call states_elec_get_state(st, mesh, ist, ik,temp_state1)
            do idim = 1, st%d%dim
              psi(ist, idim, 1:mesh%np) = temp_state1(1:mesh%np, idim)
            end do
            call states_elec_get_state(hm_st, mesh, ist, ik,temp_state1)
            do idim = 1, st%d%dim
              hpsi(ist, idim, 1:mesh%np) = temp_state1(1:mesh%np, idim)
            end do
          end if
        end do
        call comm_allreduce(mpi_world, psi)
        call comm_allreduce(mpi_world, hpsi)
        assert(mesh%np_global*st%d%dim < huge(0_int32))
        hmss(1:nst, 1:nst) = m_zero
        call zgemm( 'n',                        &
          'c',                                  &
          nst,                                  &
          nst,                                  &
          i8_to_i4(mesh%np_global*st%d%dim),    &
          cmplx(mesh%volume_element, m_zero, real64) , &
          hpsi(1, 1, 1),                        &
          ubound(hpsi, dim = 1),                &
          psi(1, 1, 1),                         &
          ubound(psi, dim = 1),                 &
          m_z0,                                 &
          hmss(1, 1),                           &
          ubound(hmss, dim = 1))
 
        hmss(1:nst,1:nst) = conjg(hmss(1:nst,1:nst))
 
        ! accumulate the Floquet integrals
        do in = -forder, forder
          do im = -forder, forder
            ii = (in + forder) * nst
            jj = (im + forder) * nst
            hfloquet(ik_count, ii+1:ii+nst, jj+1:jj+nst) =  &
              hfloquet(ik_count, ii+1:ii+nst, jj+1:jj+nst) &
              + hmss(1:nst, 1:nst) * exp(-(in - im)* m_zi * omega * it * dt)
            ! diagonal term
            if (in == im) then
              do ist = 1,nst
                hfloquet(ik_count, ii+ist, ii+ist) = hfloquet(ik_count, ii+ist, ii+ist) + in * omega
              end do
            end if
          end do
        end do
      end do !ik
    end do ! it
 
    hfloquet(:,:,:) = m_one / nt * hfloquet(:,:,:)
 
    ! diagonalize Floquet Hamiltonian
    if (downfolding) then
      ! here perform downfolding
      safe_allocate(hfloq_eff(1:nst,1:nst))
      safe_allocate(eigenval(1:nst))
      safe_allocate(bands(1:nik,1:nst))
 
      hfloq_eff(1:nst,1:nst) = m_zero
      do ik = 1, nik
        ! the HFloquet blocks are copied directly out of the super matrix
        m0 = nst ! the m=0 start position
        n0 = nst ! the n=0 start postion
        n1 = 2 * nst ! the n=+1 start postion
        hfloq_eff(1:nst,1:nst) = hfloquet(ik,n0+1:n0+nst,m0+1:m0+nst) + &
          m_one/omega*(matmul(hfloquet(ik,1:nst,m0+1:m0+nst), hfloquet(ik,n1+1:n1+nst,m0+1:m0+nst))- &
          matmul(hfloquet(ik,n1+1:n1+nst,m0+1:m0+nst), hfloquet(ik,1:nst,m0+1:m0+nst)))
 
        eigenval(1:nst) = m_zero
        call lalg_eigensolve(nst, hfloq_eff, eigenval)
        bands(ik,1:nst) = eigenval(1:nst)
      end do
      safe_deallocate_a(hfloq_eff)
    else
      ! the full Floquet
      safe_allocate(eigenval(1:nst*fdim))
      safe_allocate(bands(1:nik,1:nst*fdim))
      safe_allocate(temp(1:nst*fdim, 1:nst*fdim))
 
      do ik = 1, nik
        temp(1:nst*fdim,1:nst*fdim) = hfloquet(ik,1:nst*fdim,1:nst*fdim)
        call lalg_eigensolve(nst*fdim, temp, eigenval)
        bands(ik,1:nst*fdim) = eigenval(1:nst*fdim)
      end do
    end if
 
    !write bandstructure to file
    if (downfolding) then
      lim_nst = nst
      filename="downfolded_floquet_bands"
    else
      lim_nst = nst*fdim
      filename="floquet_bands"
    end if
    ! write bands (full or downfolded)
    if (mpi_world%rank == 0) then
      file = 987254
      file = io_open(filename, sys%namespace, action = 'write')
      do ik = 1, nik
        do ist = 1, lim_nst
          write(file,'(e13.6, 1x)',advance='no') bands(ik,ist)
        end do
        write(file,'(1x)')
      end do
      call io_close(file)
    end if
 
    if (.not. downfolding) then
      ! for the full Floquet case compute also the trivially shifted
      ! Floquet bands for reference (i.e. setting H_{nm}=0 for n!=m)
      bands(1:nik,1:nst*fdim) = m_zero
      do ik = 1,nik
        temp(1:nst*fdim, 1:nst*fdim) = m_zero
        do jj = 0, fdim - 1
          ii = jj * nst
          temp(ii+1:ii+nst, ii+1:ii+nst) = hfloquet(ik, ii+1:ii+nst, ii+1:ii+nst)
        end do
        call lalg_eigensolve(nst*fdim, temp, eigenval)
        bands(ik, 1:nst*fdim) = eigenval(1:nst*fdim)
      end do
 
      if (mpi_world%rank == 0) then
        filename='trivial_floquet_bands'
        file = io_open(filename, sys%namespace, action = 'write')
        do ik = 1, nik
          do ist = 1, lim_nst
            write(file,'(e13.6, 1x)',advance='no') bands(ik,ist)
          end do
          write(file,'(1x)')
        end do
        call io_close(file)
      end if
    end if
 
    ! reset time in Hamiltonian
    call sys%hm%update(sys%gr, global_namespace, sys%space, sys%ext_partners, time=m_zero)
 
    safe_deallocate_a(hmss)
    safe_deallocate_a(psi)
    safe_deallocate_a(hpsi)
    safe_deallocate_a(temp_state1)
    safe_deallocate_a(hfloquet)
    safe_deallocate_a(eigenval)
    safe_deallocate_a(bands)
    safe_deallocate_a(temp)
 
    pop_sub(floquet_solve_non_interacting)
 
  end subroutine floquet_solve_non_interacting
 
end program oct_floquet
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: