berry.F90

Go to the documentation of this file.

!! Copyright (C) 2010 D. Strubbe
!! Copyright (C) 2026 N. Tancogne-Dejean
!!
!! 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 batch_ops_oct_m
  use comm_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 io_oct_m, only: io_open, io_close
  use, intrinsic :: iso_fortran_env
  use kpoints_oct_m
  use lalg_adv_oct_m
  use lattice_vectors_oct_m
  use mesh_oct_m
  use mesh_batch_oct_m
  use mesh_function_oct_m
  use messages_oct_m
  use mpi_oct_m
  use namespace_oct_m
  use profiling_oct_m
  use parser_oct_m
  use phase_oct_m, only: phase_build_vector
  use smear_oct_m
  use space_oct_m
  use states_abst_oct_m
  use states_elec_oct_m
  use states_elec_dim_oct_m
  use types_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use utils_oct_m
  use v_ks_oct_m
  use wfs_elec_oct_m
 
  implicit none
 
  private
  public ::                 &
    berry_t,                &
    berry_init,             &
    berry_perform_internal_scf, &
    calc_dipole,            &
    berry_dipole,           &
    berry_write_polarization, &
    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_write_polarization(namespace, space, mesh, st, ions, kpoints, iunit)
    type(namespace_t),    intent(in) :: namespace
    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
    type(kpoints_t),      intent(in) :: kpoints
    integer,              intent(in) :: iunit
 
    integer :: idir, pol_iunit
    real(real64) :: berry_phase(space%dim)
    real(real64) :: electronic_dipole(space%dim), ionic_dipole(space%dim), total_dipole(space%dim)
    real(real64) :: electronic_pol(space%dim), ionic_pol(space%dim), total_pol(space%dim)
    real(real64) :: pol_quantum(space%dim), volume
 
    push_sub(berry_write_polarization)
 
    if (.not. space%is_periodic()) then
      message(1) = "SCFCalculatePolarization requires a periodic system."
      call messages_fatal(1, namespace=namespace)
    end if
 
    if (kpoints%use_symmetries) then
      message(1) = "SCFCalculatePolarization requires the full k-point grid."
      message(2) = "Set KPointsUseSymmetries = no."
      call messages_fatal(2, namespace=namespace)
    end if
 
    if (bitand(kpoints%method, kpoints_monkh_pack) == 0 .or. &
      bitand(kpoints%method, kpoints_user) /= 0 .or. bitand(kpoints%method, kpoints_path) /= 0) then
      message(1) = "SCFCalculatePolarization currently requires a Monkhorst-Pack KPointsGrid."
      call messages_fatal(1, namespace=namespace)
    end if
 
    if (kpoints%reduced%npoints /= kpoints%full%npoints) then
      message(1) = "SCFCalculatePolarization requires unreduced k-points."
      call messages_fatal(1, namespace=namespace)
    end if
 
    ! With more than one shift the grid labels no longer identify a k-point.
    if (kpoints%full%nshifts /= 1) then
      message(1) = "SCFCalculatePolarization requires a KPointsGrid with a single shift."
      call messages_fatal(1, namespace=namespace)
    end if
 
    if (st%parallel_in_states) then
      call messages_not_implemented("SCFCalculatePolarization with states parallelization", namespace=namespace)
    end if
 
    if (.not. smear_is_semiconducting(st%smear) .and. .not. st%smear%integral_occs) then
      message(1) = "SCFCalculatePolarization requires integer occupations."
      call messages_fatal(1, namespace=namespace)
    end if
 
    call berry_polarization_phase(namespace, space, mesh, st, ions%latt, kpoints, berry_phase)
 
    electronic_dipole = m_zero
    do idir = 1, space%periodic_dim
      electronic_dipole(1:space%dim) = electronic_dipole(1:space%dim) &
        + berry_phase(idir) * ions%latt%rlattice(1:space%dim, idir) / (m_two*m_pi)
    end do
 
    ! ions_dipole() returns P_PROTON_CHARGE * sum_I Z_I R_I, that is minus the
    ! ionic dipole, following the same convention as in calc_dipole below.
    ! The result is defined up to a quantum of polarization
    ionic_dipole = ions%latt%fold_into_cell(-ions%dipole())
    total_dipole = ions%latt%fold_into_cell(electronic_dipole + ionic_dipole)
 
    volume = ions%latt%rcell_volume
    electronic_pol = electronic_dipole / volume
    ionic_pol = ionic_dipole / volume
    total_pol = total_dipole / volume
 
    pol_quantum = m_zero
    do idir = 1, space%periodic_dim
      pol_quantum(idir) = norm2(ions%latt%rlattice(1:space%dim, idir)) / volume
    end do
 
    if (st%system_grp%is_root()) then
      pol_iunit = io_open(static_dir//'polarization', namespace, action='write', position='asis')
      call berry_write_polarization_info(space, pol_iunit, berry_phase, electronic_dipole, ionic_dipole, total_dipole, &
        electronic_pol, ionic_pol, total_pol, pol_quantum)
      call io_close(pol_iunit)
 
      call berry_write_total_polarization_info(space, iunit, total_pol)
    end if
 
    pop_sub(berry_write_polarization)
  end subroutine berry_write_polarization
 
  ! ---------------------------------------------------------
  subroutine berry_polarization_phase(namespace, space, mesh, st, latt, kpoints, berry_phase)
    type(namespace_t),          intent(in)  :: namespace
    class(space_t),             intent(in)  :: space
    class(mesh_t),              intent(in)  :: mesh
    type(states_elec_t),        intent(in)  :: st
    type(lattice_vectors_t),    intent(in)  :: latt
    type(kpoints_t),            intent(in)  :: kpoints
    real(real64), contiguous,   intent(out) :: berry_phase(:)
 
    integer :: idir, iq, iq_current, iq_next, nstrings
    logical :: have_reference
    logical, allocatable :: visited(:)
    integer, allocatable :: neighbor(:, :)
    integer, allocatable :: gvectors(:, :, :)
    real(real64) :: string_phase, reference_phase
    complex(real64) :: string_prod, det
 
    push_sub(berry_polarization_phase)
 
    berry_phase = m_zero
    safe_allocate(visited(1:st%nik))
 
    call berry_build_kpoint_strings(namespace, space, kpoints, neighbor, gvectors)
 
    do idir = 1, space%periodic_dim
      visited = .false.
      have_reference = .false.
      reference_phase = m_zero
      nstrings = kpoints%full%npoints / kpoints%nik_axis(idir)
      if (nstrings < 1) then
        message(1) = "Invalid k-point grid for Berry-phase polarization."
        call messages_fatal(1, namespace=namespace)
      end if
 
      do iq = 1, st%nik
        if (visited(iq)) cycle
 
        string_prod = m_z1
        iq_current = iq
        do
          visited(iq_current) = .true.
          call berry_neighbor_kpoint(namespace, space, mesh, latt, st, neighbor(:, idir), &
            gvectors(:, :, idir), iq_current, iq_next, det)
          string_prod = string_prod * det
          iq_current = iq_next
          if (visited(iq_current)) exit
        end do
 
        if (abs(string_prod) > m_tiny) then
          string_phase = aimag(log(string_prod))
 
          ! The phase of a string is only defined modulo 2*pi, and the logarithm
          ! returns it on the principal branch. Averaging those values directly is
          ! wrong as soon as the strings of a direction straddle the branch cut:
          ! the average then jumps by a fraction of the polarization quantum when
          ! the transverse sampling changes. Each string is therefore brought onto
          ! the branch of the first string of the direction before averaging.
          if (have_reference) then
            string_phase = string_phase &
              + m_two*m_pi*anint((reference_phase - string_phase)/(m_two*m_pi))
          else
            reference_phase = string_phase
            have_reference = .true.
          end if
 
          berry_phase(idir) = berry_phase(idir) + string_phase
        end if
      end do
 
      berry_phase(idir) = berry_phase(idir) / nstrings
    end do
 
    safe_deallocate_a(visited)
    safe_deallocate_a(neighbor)
    safe_deallocate_a(gvectors)
 
    pop_sub(berry_polarization_phase)
  end subroutine berry_polarization_phase
 
  ! ---------------------------------------------------------
  subroutine berry_neighbor_kpoint(namespace, space, mesh, latt, st, neighbor, gvectors, iq, iq_next, det)
    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
    type(states_elec_t),     intent(in)  :: st
    integer, contiguous,     intent(in)  :: neighbor(:)
    integer, contiguous,     intent(in)  :: gvectors(:, :)
    integer,                 intent(in)  :: iq
    integer,                 intent(out) :: iq_next
    complex(real64),         intent(out) :: det
 
    integer :: ik, ik_next, spin, gvector(space%dim), nocc
 
    push_sub(berry_neighbor_kpoint)
 
    ik = st%d%get_kpoint_index(iq)
    spin = st%d%get_spin_index(iq)
    ik_next = neighbor(ik)
    gvector(1:space%dim) = gvectors(1:space%dim, ik)
 
    if (st%d%ispin == spin_polarized) then
      iq_next = (ik_next - 1)*2 + spin
    else
      iq_next = ik_next
    end if
 
    nocc = berry_occupied_states(namespace, st, iq, iq_next)
    det = berry_overlap_det(st, namespace, space, mesh, latt, nocc, iq, iq_next, gvector)
 
    pop_sub(berry_neighbor_kpoint)
  end subroutine berry_neighbor_kpoint
 
  ! ---------------------------------------------------------
  subroutine berry_build_kpoint_strings(namespace, space, kpoints, neighbor, gvectors)
    type(namespace_t),    intent(in)  :: namespace
    class(space_t),       intent(in)  :: space
    type(kpoints_t),      intent(in)  :: kpoints
    integer, allocatable, intent(out) :: neighbor(:, :)
    integer, allocatable, intent(out) :: gvectors(:, :, :)
 
    integer :: ik, ik_next, idir, jdir, npoints, key
    integer :: label(space%periodic_dim)
    integer, allocatable :: lookup(:)
    real(real64) :: dk(space%dim), residual
    real(real64), parameter :: tol = 1e-8_real64
 
    push_sub(berry_build_kpoint_strings)
 
    assert(allocated(kpoints%lk123))
 
    npoints = kpoints%full%npoints
 
    ! Invert the labelling: the label of a k-point gives its position on the
    ! grid, and this gives back the k-point sitting at a given position.
    safe_allocate(lookup(1:npoints))
    lookup = 0
    do ik = 1, npoints
      key = berry_kpoint_key(space, kpoints, kpoints%lk123(ik, 1:space%periodic_dim))
      if (key < 1 .or. key > npoints) then
        message(1) = "Unable to construct Berry-phase k-point strings from KPointsGrid."
        call messages_fatal(1, namespace=namespace)
      end if
      if (lookup(key) /= 0) then
        message(1) = "Repeated k-points in the Berry-phase k-point grid."
        call messages_fatal(1, namespace=namespace)
      end if
      lookup(key) = ik
    end do
 
    safe_allocate(neighbor(1:npoints, 1:space%periodic_dim))
    safe_allocate(gvectors(1:space%dim, 1:npoints, 1:space%periodic_dim))
 
    do idir = 1, space%periodic_dim
      dk = m_zero
      dk(idir) = m_one/kpoints%nik_axis(idir)
 
      do ik = 1, npoints
        label(1:space%periodic_dim) = kpoints%lk123(ik, 1:space%periodic_dim)
        label(idir) = label(idir) - 1
        if (label(idir) < 1) label(idir) = label(idir) + kpoints%nik_axis(idir)
 
        ik_next = lookup(berry_kpoint_key(space, kpoints, label))
        neighbor(ik, idir) = ik_next
 
        ! k - dk and the stored neighbor describe the same state, so they differ
        ! by a reciprocal-lattice vector, which is an integer in reduced
        ! coordinates.
        do jdir = 1, space%dim
          residual = kpoints%full%red_point(jdir, ik_next) - kpoints%full%red_point(jdir, ik) + dk(jdir)
          gvectors(jdir, ik, idir) = nint(residual)
          if (abs(residual - gvectors(jdir, ik, idir)) > tol) then
            message(1) = "Unable to construct Berry-phase k-point strings from KPointsGrid."
            call messages_fatal(1, namespace=namespace)
          end if
        end do
      end do
    end do
 
    safe_deallocate_a(lookup)
 
    pop_sub(berry_build_kpoint_strings)
  end subroutine berry_build_kpoint_strings
 
  ! ---------------------------------------------------------
  integer pure function berry_kpoint_key(space, kpoints, label) result(key)
    class(space_t),  intent(in) :: space
    type(kpoints_t), intent(in) :: kpoints
    integer,         intent(in) :: label(:)
 
    integer :: idir, stride
 
    key = 1
    stride = 1
    do idir = 1, space%periodic_dim
      key = key + (label(idir) - 1)*stride
      stride = stride*kpoints%nik_axis(idir)
    end do
  end function berry_kpoint_key
 
  ! ---------------------------------------------------------
  integer function berry_occupied_states(namespace, st, iq, iq_next) result(nocc)
    type(namespace_t),   intent(in) :: namespace
    type(states_elec_t), intent(in) :: st
    integer,             intent(in) :: iq, iq_next
 
    integer :: ist, nocc_next
 
    push_sub(berry_occupied_states)
 
    nocc = 0
    nocc_next = 0
    do ist = 1, st%nst
      if (st%occ(ist, iq) > m_epsilon) nocc = ist
      if (st%occ(ist, iq_next) > m_epsilon) nocc_next = ist
    end do
 
    if (nocc /= nocc_next) then
      message(1) = "Berry-phase polarization requires the same occupied manifold at neighboring k-points."
      call messages_fatal(1, namespace=namespace)
    end if
 
    pop_sub(berry_occupied_states)
  end function berry_occupied_states
 
  ! ---------------------------------------------------------
  complex(real64) function berry_overlap_det(st, namespace, space, mesh, latt, nocc, iq, iq_next, gvector) result(det)
    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
    integer,                 intent(in) :: nocc
    integer,                 intent(in) :: iq
    integer,                 intent(in) :: iq_next
    integer, contiguous,     intent(in) :: gvector(:)
 
    complex(real64), allocatable :: matrix(:, :)
 
    push_sub(berry_overlap_det)
 
    if (nocc == 0) then
      det = m_z1
      pop_sub(berry_overlap_det)
      return
    end if
 
    safe_allocate(matrix(1:nocc, 1:nocc))
    call berry_overlap_matrix(st, namespace, space, mesh, latt, nocc, iq, iq_next, gvector, matrix)
    det = lalg_determinant(nocc, matrix, preserve_mat=.false.) ** st%smear%el_per_state
    safe_deallocate_a(matrix)
 
    pop_sub(berry_overlap_det)
  end function berry_overlap_det
 
  ! ---------------------------------------------------------
  subroutine berry_overlap_matrix(st, namespace, space, mesh, latt, nocc, iq, iq_next, gvector, matrix)
    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
    integer,                     intent(in)  :: nocc
    integer,                     intent(in)  :: iq
    integer,                     intent(in)  :: iq_next
    integer, contiguous,         intent(in)  :: gvector(:)
    complex(real64), contiguous, intent(out) :: matrix(:, :)
 
    integer :: ist, ib, jb, jst_start, jst_end, nst_block, sender
    logical :: local_iq, local_iq_next
    real(real64) :: gred(space%dim), gcart(space%dim)
    complex(real64), allocatable :: psi_next(:, :, :), dot(:, :), phase(:)
    type(wfs_elec_t) :: psib_next, psib_cur
 
    push_sub(berry_overlap_matrix)
 
    assert(.not. st%parallel_in_states)
 
    matrix = m_z0
    local_iq = (iq >= st%d%kpt%start .and. iq <= st%d%kpt%end)
    local_iq_next = (iq_next >= st%d%kpt%start .and. iq_next <= st%d%kpt%end)
 
    if (any(gvector(1:space%dim) /= 0)) then
      safe_allocate(phase(1:mesh%np))
      gred(1:space%dim) = real(gvector(1:space%dim), real64)
      call kpoints_to_absolute(latt, gred, gcart)
      call phase_build_vector(mesh, mesh%np, gcart(1:space%dim), .false., phase)
    end if
 
    if (local_iq) then
      ! The result of the batched dot product is indexed by the global state indices
      ! of the blocks, and the block containing the state nocc may extend beyond it.
      safe_allocate(dot(1:st%nst, 1:nocc))
      dot = m_z0
    end if
 
    do jb = st%group%block_start, st%group%block_end
      jst_start = st%group%block_range(jb, 1)
      if (jst_start > nocc) cycle
      jst_end = min(st%group%block_range(jb, 2), nocc)
      nst_block = jst_end - jst_start + 1
 
      safe_allocate(psi_next(1:mesh%np, 1:st%d%dim, 1:nst_block))
 
      if (local_iq_next) then
        do ist = jst_start, jst_end
          call states_elec_get_state(st, mesh, ist, iq_next, psi_next(:, :, ist - jst_start + 1))
        end do
      end if
      if (st%d%kpt%parallel) then
        sender = st%d%kpt%process(iq_next)
        call st%d%kpt%mpi_grp%bcast(psi_next, mesh%np*st%d%dim*nst_block, mpi_double_complex, sender)
      end if
 
      if (local_iq) then
        ! The batch is labeled with iq, not iq_next, as the batches entering the dot product
        ! below must carry the same k-point index.
        call wfs_elec_init(psib_next, st%d%dim, jst_start, jst_end, psi_next, iq)
        if (st%group%psib(st%group%block_start, iq)%is_packed()) call psib_next%do_pack(copy = .true.)
        if (allocated(phase)) call batch_mul_mf(mesh%np, phase, psib_next, psib_next)
 
        do ib = st%group%block_start, st%group%block_end
          if (st%group%block_range(ib, 1) > nocc) cycle
          if (states_are_real(st)) then
            ! The phase-shifted states are always complex, so real ground-state
            ! wavefunctions must be promoted before taking the overlap.
            call st%group%psib(ib, iq)%copy_to(psib_cur, copy_data = .true., dest_type = type_cmplx)
            call zmesh_batch_dotp_matrix(mesh, psib_cur, psib_next, dot, reduce = .false.)
            call psib_cur%end()
          else
            call zmesh_batch_dotp_matrix(mesh, st%group%psib(ib, iq), psib_next, dot, reduce = .false.)
          end if
        end do
 
        call psib_next%end(copy = .false.)
      end if
 
      safe_deallocate_a(psi_next)
    end do
 
    if (local_iq) then
      matrix(1:nocc, 1:nocc) = dot(1:nocc, 1:nocc)
      safe_deallocate_a(dot)
    end if
 
    call mesh%allreduce(matrix)
    if (st%d%kpt%parallel) call comm_allreduce(st%d%kpt%mpi_grp, matrix)
 
    safe_deallocate_a(phase)
 
    pop_sub(berry_overlap_matrix)
  end subroutine berry_overlap_matrix
 
  ! ---------------------------------------------------------
  subroutine berry_write_polarization_info(space, iunit, berry_phase, electronic_dipole, ionic_dipole, total_dipole, &
    electronic_pol, ionic_pol, total_pol, pol_quantum)
    class(space_t), intent(in) :: space
    integer,        intent(in) :: iunit
    real(real64), contiguous, intent(in) :: berry_phase(:)
    real(real64), contiguous, intent(in) :: electronic_dipole(:)
    real(real64), contiguous, intent(in) :: ionic_dipole(:)
    real(real64), contiguous, intent(in) :: total_dipole(:)
    real(real64), contiguous, intent(in) :: electronic_pol(:)
    real(real64), contiguous, intent(in) :: ionic_pol(:)
    real(real64), contiguous, intent(in) :: total_pol(:)
    real(real64), contiguous, intent(in) :: pol_quantum(:)
 
    type(unit_t) :: unit_pol, unit_c_m2
    integer :: idir
 
    push_sub(berry_write_polarization_info)
 
    call berry_polarization_units(space, unit_pol, unit_c_m2)
 
    write(message(1), '(a)') 'Berry-phase polarization:'
    write(message(2), '(a)') 'Defined modulo the polarization quantum.'
    write(message(3), '(a)') 'Ionic and total values are reduced to the nearest polarization branch.'
    call messages_info(3, iunit=iunit)
 
    write(message(1), '(a)') '  Berry phases [rad]:'
    do idir = 1, space%periodic_dim
      write(message(1+idir), '(6x,3a,es18.10)') '<', index2axis(idir), '> = ', berry_phase(idir)
    end do
    call messages_info(1+space%periodic_dim, iunit=iunit)
 
    write(message(1), '(a21,a20)') '  Electronic dipole:', '[' // trim(units_abbrev(units_out%length)) // ']'
    do idir = 1, space%dim
      write(message(1+idir), '(6x,3a,es18.10)') '<', index2axis(idir), '> = ', &
        units_from_atomic(units_out%length, electronic_dipole(idir))
    end do
    call messages_info(1+space%dim, iunit=iunit)
 
    write(message(1), '(a21,a20)') '  Ionic dipole:', '[' // trim(units_abbrev(units_out%length)) // ']'
    do idir = 1, space%dim
      write(message(1+idir), '(6x,3a,es18.10)') '<', index2axis(idir), '> = ', &
        units_from_atomic(units_out%length, ionic_dipole(idir))
    end do
    call messages_info(1+space%dim, iunit=iunit)
 
    write(message(1), '(a21,a20)') '  Total dipole:', '[' // trim(units_abbrev(units_out%length)) // ']'
    do idir = 1, space%dim
      write(message(1+idir), '(6x,3a,es18.10)') '<', index2axis(idir), '> = ', &
        units_from_atomic(units_out%length, total_dipole(idir))
    end do
    call messages_info(1+space%dim, iunit=iunit)
 
    write(message(1), '(a,t35,a,t52,a)') '  Electronic polarization:', &
      '[' // trim(units_abbrev(unit_pol)) // ']', '[' // trim(units_abbrev(unit_c_m2)) // ']'
    do idir = 1, space%dim
      write(message(1+idir), '(6x,3a,t35,es14.5,t52,es14.5)') '<', index2axis(idir), '> = ', &
        units_from_atomic(unit_pol, electronic_pol(idir)), units_from_atomic(unit_c_m2, electronic_pol(idir))
    end do
    call messages_info(1+space%dim, iunit=iunit)
 
    write(message(1), '(a,t35,a,t52,a)') '  Ionic polarization:', &
      '[' // trim(units_abbrev(unit_pol)) // ']', '[' // trim(units_abbrev(unit_c_m2)) // ']'
    do idir = 1, space%dim
      write(message(1+idir), '(6x,3a,t35,es14.5,t52,es14.5)') '<', index2axis(idir), '> = ', &
        units_from_atomic(unit_pol, ionic_pol(idir)), units_from_atomic(unit_c_m2, ionic_pol(idir))
    end do
    call messages_info(1+space%dim, iunit=iunit)
 
    write(message(1), '(a,t35,a,t52,a)') '  Total polarization:', &
      '[' // trim(units_abbrev(unit_pol)) // ']', '[' // trim(units_abbrev(unit_c_m2)) // ']'
    do idir = 1, space%dim
      write(message(1+idir), '(6x,3a,t35,es14.5,t52,es14.5)') '<', index2axis(idir), '> = ', &
        units_from_atomic(unit_pol, total_pol(idir)), units_from_atomic(unit_c_m2, total_pol(idir))
    end do
    call messages_info(1+space%dim, iunit=iunit)
 
    write(message(1), '(a,t30,a,t52,a)') '  Polarization quantum:', &
      '[' // trim(units_abbrev(unit_pol)) // ']', '[' // trim(units_abbrev(unit_c_m2)) // ']'
    do idir = 1, space%periodic_dim
      write(message(1+idir), '(6x,3a,t30,es14.5,t52,es14.5)') '<', index2axis(idir), '> = ', &
        units_from_atomic(unit_pol, pol_quantum(idir)), units_from_atomic(unit_c_m2, pol_quantum(idir))
    end do
    call messages_info(1+space%periodic_dim, iunit=iunit)
 
    pop_sub(berry_write_polarization_info)
  end subroutine berry_write_polarization_info
 
  ! ---------------------------------------------------------
  subroutine berry_write_total_polarization_info(space, iunit, total_pol)
    class(space_t), intent(in) :: space
    integer,        intent(in) :: iunit
    real(real64), contiguous, intent(in) :: total_pol(:)
 
    type(unit_t) :: unit_pol, unit_c_m2
    integer :: idir
 
    push_sub(berry_write_total_polarization_info)
 
    call berry_polarization_units(space, unit_pol, unit_c_m2)
 
    write(message(1), '(a,t30,a,t52,a)') 'Total polarization:', &
      '[' // trim(units_abbrev(unit_pol)) // ']', '[' // trim(units_abbrev(unit_c_m2)) // ']'
    do idir = 1, space%dim
      write(message(1+idir), '(4x,3a,t30,es14.5,t52,es14.5)') '<', index2axis(idir), '> = ', &
        units_from_atomic(unit_pol, total_pol(idir)), units_from_atomic(unit_c_m2, total_pol(idir))
    end do
    call messages_info(1+space%dim, iunit=iunit)
 
    pop_sub(berry_write_total_polarization_info)
  end subroutine berry_write_total_polarization_info
 
  ! ---------------------------------------------------------
  subroutine berry_polarization_units(space, unit_pol, unit_c_m2)
    class(space_t), intent(in)  :: space
    type(unit_t),   intent(out) :: unit_pol
    type(unit_t),   intent(out) :: unit_c_m2
 
    push_sub(berry_polarization_units)
 
    unit_pol = units_out%length**(1-space%periodic_dim)
    unit_c_m2%factor = m_one / (1.602176634e-19_real64 / (p_a_b*1.0e-10_real64)**2)
    unit_c_m2%abbrev = 'C/m^2'
    unit_c_m2%name = 'coulomb per square meter'
 
    pop_sub(berry_polarization_units)
  end subroutine berry_polarization_units
 
  ! ---------------------------------------------------------
  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
    ! The Berry potential and dipole are implemented for a single unpolarized k-point only
    if (st%d%ispin == spinors) then
      call messages_not_implemented("Berry phase with spinors", namespace=namespace)
    end if
    if (st%d%ispin == spin_polarized) then
      call messages_not_implemented("Berry phase with spin polarization", 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, space, ext_partners, 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, namespace, 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
 
  ! ---------------------------------------------------------
  subroutine calc_dipole(dipole, space, mesh, st, ions)
    real(real64), contiguous, 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%spin_channels), 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%spin_channels
      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, ions%namespace, 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, namespace, mesh, latt, space, dir) result(dipole)
    type(states_elec_t), intent(in)     :: st
    type(namespace_t),   intent(in)     :: namespace
    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, namespace, mesh, latt, space, dir, ik)
      !Check in case det is zero. This should never occur in practical calculations
      if (abs(det) > m_tiny) then
        dipole = dipole + aimag(log(det))
      end if
    end do
 
    dipole = -(norm2(latt%rlattice(:,dir)) / (m_two*m_pi)) * dipole
 
    pop_sub(berry_dipole)
  end function berry_dipole
 
 
  ! ---------------------------------------------------------
  complex(real64) function berry_phase_det(st, namespace, mesh, latt, space, dir, ik) result(det)
    type(states_elec_t),     intent(in) :: st
    type(namespace_t),       intent(in) :: namespace
    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)
 
    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
 
    gvector(:) = 0
    gvector(dir) = 1
    det = berry_overlap_det(st, namespace, space, mesh, latt, noccst, ik, ik, gvector)
 
    pop_sub(berry_phase_det)
  end function berry_phase_det
 
 
  ! ---------------------------------------------------------
  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), contiguous, intent(in) :: e_field(:)
    real(real64), contiguous, intent(out) :: pot(:,:)
 
    integer :: ispin, idir, ip
    complex(real64) :: factor, det, phase_factor
 
    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, namespace, 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
          ! Note: this is only correct for orthogonal cells
          phase_factor = m_two * m_pi * m_zi / norm2(latt%rlattice(:,idir))
          !$omp parallel do
          do ip = 1, mesh%np
            pot(ip, ispin) = pot(ip, ispin) + &
              aimag(factor * exp(phase_factor * mesh%x_t(ip, idir)))
          end do
          !$omp end parallel do
        end if
      end do
    end do
 
    pop_sub(berry_potential)
  end subroutine berry_potential
 
 
  ! ---------------------------------------------------------
  real(real64) function berry_energy_correction(st, namespace, space, mesh, latt, e_field, vberry) result(delta)
    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), contiguous, 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, namespace, 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: