berry.F90

Go to the documentation of this file.

!! Copyright (C) 2010 D. Strubbe
!!
!! 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"
 
module berry_oct_m
  use debug_oct_m
  use eigensolver_oct_m
  use electron_space_oct_m
  use global_oct_m
  use grid_oct_m
  use hamiltonian_elec_oct_m
  use interaction_partner_oct_m
  use ions_oct_m
  use, intrinsic :: iso_fortran_env
  use lalg_adv_oct_m
  use lattice_vectors_oct_m
  use mesh_oct_m
  use mesh_function_oct_m
  use messages_oct_m
  use namespace_oct_m
  use profiling_oct_m
  use parser_oct_m
  use smear_oct_m
  use space_oct_m
  use states_elec_oct_m
  use states_elec_dim_oct_m
  use v_ks_oct_m
 
  implicit none
 
  private
  public ::                 &
    berry_t,                &
    berry_init,             &
    berry_perform_internal_scf, &
    calc_dipole,            &
    berry_dipole,           &
    berry_potential,        &
    berry_energy_correction
 
  type berry_t
    private
 
    integer :: max_iter_berry
  end type berry_t
contains
 
  ! ---------------------------------------------------------
  subroutine berry_init(this, namespace)
    type(berry_t),     intent(inout) :: this
    type(namespace_t), intent(in)    :: namespace
 
    push_sub(berry_init)
 
    !%Variable MaximumIterBerry
    !%Type integer
    !%Default 10
    !%Section SCF::Convergence
    !%Description
    !% Maximum number of iterations for the Berry potential, within each SCF iteration.
    !% Only applies if a <tt>StaticElectricField</tt> is applied in a periodic direction.
    !% The code will move on to the next SCF iteration even if convergence
    !% has not been achieved. -1 means unlimited.
    !%End
    call parse_variable(namespace, 'MaximumIterBerry', 10, this%max_iter_berry)
    if (this%max_iter_berry < 0) this%max_iter_berry = huge(this%max_iter_berry)
 
 
    pop_sub(berry_init)
  end subroutine berry_init
 
  ! ---------------------------------------------------------
  subroutine berry_perform_internal_scf(this, namespace, space, eigensolver, gr, st, hm, &
    iter, ks, ions, ext_partners)
    type(berry_t),            intent(in)    :: this
    type(namespace_t),        intent(in)    :: namespace
    type(electron_space_t),   intent(in)    :: space
    type(eigensolver_t),      intent(inout) :: eigensolver
    type(grid_t),             intent(in)    :: gr
    type(states_elec_t),      intent(inout) :: st
    type(hamiltonian_elec_t), intent(inout) :: hm
    integer,                  intent(in)    :: iter
    type(v_ks_t),             intent(inout) :: ks
    type(ions_t),             intent(in)    :: ions
    type(partner_list_t),     intent(in)    :: ext_partners
 
    integer :: iberry, idir
    logical :: berry_conv
    real(real64) :: dipole_prev(space%dim), dipole(space%dim)
    real(real64), parameter :: tol = 1e-5_real64
 
    push_sub(berry_perform_internal_scf)
 
    assert(allocated(hm%vberry))
 
    if (st%parallel_in_states) then
      call messages_not_implemented("Berry phase parallel in states", namespace=namespace)
    end if
 
    call calc_dipole(dipole, space, gr, st, ions)
 
    do iberry = 1, this%max_iter_berry
      eigensolver%converged = 0
      call eigensolver%run(namespace, gr, st, hm, iter)
 
      !Calculation of the Berry potential
      call berry_potential(st, namespace, space, gr, hm%ions%latt, hm%ep%e_field, hm%vberry)
 
      !Calculation of the corresponding energy
      hm%energy%berry = berry_energy_correction(st, space, gr, hm%ions%latt, &
        hm%ep%e_field(1:space%periodic_dim), hm%vberry(1:gr%np, 1:hm%d%nspin))
 
      !We recompute the KS potential
      call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_current=.false.)
 
      dipole_prev = dipole
      call calc_dipole(dipole, space, gr, st, ions)
      write(message(1),'(a,9f12.6)') 'Dipole = ', dipole(1:space%dim)
      call messages_info(1, namespace=namespace)
 
      berry_conv = .true.
      do idir = 1, space%periodic_dim
        if (abs(dipole_prev(idir)) > 1e-10_real64) then
          berry_conv = berry_conv .and. (abs(dipole(idir) - dipole_prev(idir)) < tol &
            .or.(abs((dipole(idir) - dipole_prev(idir)) / dipole_prev(idir)) < tol))
        else
          berry_conv = berry_conv .and. (abs(dipole(idir) - dipole_prev(idir)) < tol)
        end if
      end do
 
      if (berry_conv) exit
    end do
 
    pop_sub(berry_perform_internal_scf)
  end subroutine berry_perform_internal_scf
 
  ! ---------------------------------------------------------
  !TODO: This should be a method of the electronic system directly
  ! that can be exposed
  subroutine calc_dipole(dipole, space, mesh, st, ions)
    real(real64),          intent(out)   :: dipole(:)
    class(space_t),        intent(in)    :: space
    class(mesh_t),         intent(in)    :: mesh
    type(states_elec_t),   intent(in)    :: st
    type(ions_t),          intent(in)    :: ions
 
    integer :: ispin, idir
    real(real64) :: e_dip(space%dim + 1, st%d%nspin), n_dip(space%dim), nquantumpol
 
    push_sub(calc_dipole)
 
    assert(.not. ions%latt%nonorthogonal)
 
    dipole(1:space%dim) = m_zero
 
    do ispin = 1, st%d%nspin
      call dmf_multipoles(mesh, st%rho(:, ispin), 1, e_dip(:, ispin))
    end do
 
    n_dip = ions%dipole()
 
    do idir = 1, space%dim
      ! in periodic directions use single-point Berry`s phase calculation
      if (idir  <=  space%periodic_dim) then
        dipole(idir) = -n_dip(idir) - berry_dipole(st, mesh, ions%latt, space, idir)
 
        ! use quantum of polarization to reduce to smallest possible magnitude
        nquantumpol = nint(dipole(idir)/norm2(ions%latt%rlattice(:, idir)))
        dipole(idir) = dipole(idir) - nquantumpol * norm2(ions%latt%rlattice(:, idir))
        ! in aperiodic directions use normal dipole formula
 
      else
        e_dip(idir + 1, 1) = sum(e_dip(idir + 1, :))
        dipole(idir) = -n_dip(idir) - e_dip(idir + 1, 1)
      end if
    end do
 
    pop_sub(calc_dipole)
  end subroutine calc_dipole
 
 
  ! ---------------------------------------------------------
  real(real64) function berry_dipole(st, mesh, latt, space, dir) result(dipole)
    type(states_elec_t), intent(in)     :: st
    class(mesh_t),       intent(in)     :: mesh
    type(lattice_vectors_t), intent(in) :: latt
    class(space_t),      intent(in)     :: space
    integer,             intent(in)     :: dir
 
    integer :: ik
    complex(real64) :: det
 
    push_sub(berry_dipole)
 
    !There is a missing reduction here.
    !However, as we are limited to a single k-point at the moment, this is not a problem.
    assert(st%nik == 1)
 
    dipole = m_zero
    do ik = st%d%kpt%start, st%d%kpt%end ! determinants for different spins and k-points multiply since matrix is block-diagonal
      det  = berry_phase_det(st, mesh, latt, space, dir, ik)
      dipole = dipole + aimag(log(det))
    end do
 
    dipole = -(norm2(latt%rlattice(:,dir)) / (m_two*m_pi)) * dipole
 
    pop_sub(berry_dipole)
  end function berry_dipole
 
 
  ! ---------------------------------------------------------
  !! \f[
  !! det <\psi_i|exp(-i(2*\pi/L)x)|\psi_j>
  !! \f]
  !! E Yaschenko, L Fu, L Resca, R Resta, Phys. Rev. B 58, 1222-1229 (1998)
  complex(real64) function berry_phase_det(st, mesh, latt, space, dir, ik) result(det)
    type(states_elec_t),     intent(in) :: st
    class(mesh_t),           intent(in) :: mesh
    type(lattice_vectors_t), intent(in) :: latt
    class(space_t),          intent(in) :: space
    integer,                 intent(in) :: dir
    integer,                 intent(in) :: ik
 
    integer :: ist, noccst, gvector(3)
    complex(real64), allocatable :: matrix(:, :), tmp(:), phase(:)
    complex(real64), allocatable :: psi(:, :), psi2(:, :)
 
    push_sub(berry_phase_det)
 
    ! find how many states are occupied on this k-point. Formalism only works if semiconducting anyway.
    noccst = 0
    do ist = 1, st%nst
      if (st%occ(ist, ik) > m_epsilon) noccst = ist
    end do
 
    safe_allocate(matrix(1:noccst, 1:noccst))
 
    gvector(:) = 0
    gvector(dir) = 1
    call berry_phase_matrix(st, mesh, latt, space, noccst, ik, ik, gvector, matrix)
 
    if (noccst > 0) then
      det = lalg_determinant(noccst, matrix, preserve_mat=.false.) ** st%smear%el_per_state
    else
      det = m_one
    end if
 
    safe_deallocate_a(matrix)
    safe_deallocate_a(tmp)
    safe_deallocate_a(phase)
    safe_deallocate_a(psi)
    safe_deallocate_a(psi2)
 
    pop_sub(berry_phase_det)
  end function berry_phase_det
 
 
  ! ---------------------------------------------------------
  subroutine berry_phase_matrix(st, mesh, latt, space, nst, ik, ik2, gvector, matrix)
    type(states_elec_t),     intent(in)  :: st
    class(mesh_t),           intent(in)  :: mesh
    type(lattice_vectors_t), intent(in)  :: latt
    class(space_t),          intent(in)  :: space
    integer,                 intent(in)  :: nst
    integer,                 intent(in)  :: ik
    integer,                 intent(in)  :: ik2
    integer,                 intent(in)  :: gvector(:)
    complex(real64),         intent(out) :: matrix(:,:)
 
    integer :: ist, ist2, idim, ip, idir
    complex(real64), allocatable :: tmp(:), phase(:)
    complex(real64), allocatable :: psi(:, :), psi2(:, :)
    real(real64) :: norm
    ! FIXME: real/cplx versions
 
    push_sub(berry_phase_matrix)
 
    assert(.not. st%parallel_in_states)
    assert(.not. latt%nonorthogonal)
 
    safe_allocate(tmp(1:mesh%np))
    safe_allocate(phase(1:mesh%np))
    safe_allocate(psi(1:mesh%np, 1:st%d%dim))
    safe_allocate(psi2(1:mesh%np, 1:st%d%dim))
 
    phase = m_zero
    do idir = 1, space%dim
      if (gvector(idir) == 0) cycle
      norm = norm2(latt%rlattice(:,idir))
      if (idir > space%periodic_dim) norm = norm * m_two * mesh%box%bounding_box_l(idir)
      do ip = 1, mesh%np
        phase(ip) = phase(ip) + exp(-m_zi*gvector(idir)*(m_two*m_pi/norm)*mesh%x(ip, idir))
      end do
    end do
 
    do ist = 1, nst
      call states_elec_get_state(st, mesh, ist, ik, psi)
      do ist2 = 1, nst
        call states_elec_get_state(st, mesh, ist2, ik2, psi2)
        matrix(ist, ist2) = m_z0
        do idim = 1, st%d%dim ! spinor components
 
          do ip = 1, mesh%np
            tmp(ip) = conjg(psi(ip, idim))*phase(ip)*psi2(ip, idim)
          end do
 
          matrix(ist, ist2) = matrix(ist, ist2) + zmf_integrate(mesh, tmp)
        end do
      end do !ist2
    end do !ist
 
    safe_deallocate_a(tmp)
    safe_deallocate_a(phase)
    safe_deallocate_a(psi)
    safe_deallocate_a(psi2)
 
    pop_sub(berry_phase_matrix)
  end subroutine berry_phase_matrix
 
 
  ! ---------------------------------------------------------
  subroutine berry_potential(st, namespace, space, mesh, latt, e_field, pot)
    type(states_elec_t),     intent(in)  :: st
    type(namespace_t),       intent(in)  :: namespace
    class(space_t),          intent(in)  :: space
    class(mesh_t),           intent(in)  :: mesh
    type(lattice_vectors_t), intent(in)  :: latt
    real(real64),            intent(in)  :: e_field(:)
    real(real64),            intent(out) :: pot(:,:)
 
    integer :: ispin, idir
    complex(real64) :: factor, det
 
    push_sub(berry_potential)
 
    if (latt%nonorthogonal) then
      call messages_not_implemented("Berry phase for non-orthogonal cells", namespace=namespace)
    end if
 
    if (st%nik > 1) then
      call messages_not_implemented("Berry phase with k-points", namespace=namespace)
    end if
 
    pot(1:mesh%np, 1:st%d%nspin) = m_zero
 
    do ispin = 1, st%d%nspin
      do idir = 1, space%periodic_dim
        if (abs(e_field(idir)) > m_epsilon) then
          ! calculate the ip-independent part first
          det = berry_phase_det(st, mesh, latt, space, idir, ispin)
          if (abs(det) > m_epsilon) then
            factor = e_field(idir) * (norm2(latt%rlattice(:,idir)) / (m_two*m_pi)) / det
          else
            ! If det = 0, mu = -infinity, so this condition should never be reached
            ! if things are working properly.
            write(message(1),*) "Divide by zero: dir = ", idir, " Berry-phase determinant = ", det
            call messages_fatal(1, namespace=namespace)
          end if
          pot(1:mesh%np, ispin) = pot(1:mesh%np, ispin) + &
            aimag(factor * exp(m_two * m_pi * m_zi * mesh%x(1:mesh%np, idir) / norm2(latt%rlattice(:,idir))))
        end if
      end do
    end do
 
    pop_sub(berry_potential)
  end subroutine berry_potential
 
 
  ! ---------------------------------------------------------
  real(real64) function berry_energy_correction(st, space, mesh, latt, e_field, vberry) result(delta)
    type(states_elec_t),     intent(in) :: st
    class(space_t),          intent(in) :: space
    class(mesh_t),           intent(in) :: mesh
    type(lattice_vectors_t), intent(in) :: latt
    real(real64),            intent(in) :: e_field(:)
    real(real64),            intent(in) :: vberry(:,:)
 
    integer :: ispin, idir
 
    push_sub(berry_energy_correction)
 
    ! first we calculate expectation value of Berry potential, to subtract off
    delta = m_zero
    do ispin = 1, st%d%nspin
      delta = delta - dmf_dotp(mesh, st%rho(1:mesh%np, ispin), vberry(1:mesh%np, ispin))
    end do
 
    ! the real energy contribution is -mu.E
    do idir = 1, space%periodic_dim
      delta = delta - berry_dipole(st, mesh, latt, space, idir) * e_field(idir)
    end do
 
    pop_sub(berry_energy_correction)
  end function berry_energy_correction
 
end module berry_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: