born_charges.F90

Go to the documentation of this file.

!! Copyright (C) 2009 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 born_charges_oct_m
  use atom_oct_m
  use debug_oct_m
  use global_oct_m
  use io_oct_m
  use messages_oct_m
  use mpi_oct_m
  use namespace_oct_m
  use parser_oct_m
  use profiling_oct_m
  use species_oct_m
  use unit_system_oct_m
  use utils_oct_m
 
  implicit none
 
  private
  public ::                &
    born_charges_t,        &
    born_charges_init,     &
    born_charges_end,      &
    born_output_charges
 
  type born_charges_t
    private
    complex(real64), allocatable, public :: charge(:, :, :)
    complex(real64), allocatable  :: sum_ideal(:,:)
    complex(real64), allocatable  :: delta(:,:)
    logical :: correct
  end type born_charges_t
 
contains
 
  ! ---------------------------------------------------------
  subroutine born_charges_init(this, namespace, natoms, val_charge, qtot, dim)
    type(born_charges_t), intent(out) :: this
    type(namespace_t),    intent(in)  :: namespace
    integer,              intent(in)  :: natoms
    real(real64),         intent(in)  :: val_charge
    real(real64),         intent(in)  :: qtot
    integer,              intent(in)  :: dim
 
    integer :: idir
 
    push_sub(born_charges_init)
 
    safe_allocate(this%charge(1:dim, 1:dim, 1:natoms))
    safe_allocate(this%sum_ideal(1:dim, 1:dim))
    safe_allocate(this%delta(1:dim, 1:dim))
    this%charge(1:dim, 1:dim, 1:natoms) = m_zero
    this%delta(1:dim, 1:dim) = m_zero
 
    this%sum_ideal = m_zero
    do idir = 1, dim
      this%sum_ideal(idir, idir) = -(val_charge + qtot) ! total charge
    end do
 
    !%Variable BornChargeSumRuleCorrection
    !%Type logical
    !%Default true
    !%Section Linear Response::Polarizabilities
    !%Description
    !% Enforce the acoustic sum rule by distributing the excess sum of Born charges equally among the atoms.
    !% Sum rule: <math>\sum_{\alpha} Z^{*}_{\alpha, i, j} = Z_{\rm tot} \delta_{ij}</math>.
    !% Violation of the sum rule may be caused by inadequate spacing, box size (in finite directions),
    !% or <i>k</i>-point sampling (in periodic directions).
    !%End
 
    call parse_variable(namespace, 'BornChargeSumRuleCorrection', .true., this%correct)
 
    pop_sub(born_charges_init)
  end subroutine born_charges_init
 
  ! ---------------------------------------------------------
  subroutine born_charges_end(this)
    type(born_charges_t), intent(inout) :: this
 
    push_sub(born_charges_end)
 
    safe_deallocate_a(this%charge)
    safe_deallocate_a(this%sum_ideal)
    safe_deallocate_a(this%delta)
 
    pop_sub(born_charges_end)
  end subroutine born_charges_end
 
  ! ---------------------------------------------------------
  subroutine correct_born_charges(this, natoms, dim)
    type(born_charges_t), intent(inout) :: this
    integer,              intent(in)    :: natoms
    integer,              intent(in)    :: dim
 
    complex(real64) :: born_sum(dim, dim)        ! the sum of born charges from the calculation
    integer :: iatom
 
    push_sub(correct_born_charges)
 
    born_sum = m_zero
 
    do iatom = 1, natoms
      born_sum = born_sum + this%charge(:, :, iatom)
    end do
 
    this%delta = (born_sum - this%sum_ideal) / natoms
 
    if (this%correct) then
      do iatom = 1, natoms
        this%charge(:, :, iatom) = this%charge(:, :, iatom) - this%delta
      end do
    end if
 
    pop_sub(correct_born_charges)
  end subroutine correct_born_charges
 
  ! ---------------------------------------------------------
  subroutine born_output_charges(this, atom, charge, natoms, namespace, dim, dirname, write_real)
    type(born_charges_t), intent(inout) :: this
    type(atom_t),         intent(in)    :: atom(:)
    real(real64),         intent(in)    :: charge(:)
    integer,              intent(in)    :: natoms
    type(namespace_t),    intent(in)    :: namespace
    integer,              intent(in)    :: dim
    character(len=*),     intent(in)    :: dirname
    logical,              intent(in)    :: write_real
    
 
    integer iatom, iunit
    real(real64) :: phase(dim, dim)
 
    push_sub(born_output_charges)
 
    call correct_born_charges(this, natoms, dim)
 
    if (mpi_grp_is_root(mpi_world)) then ! only first node outputs
      iunit = io_open(trim(dirname)//'/born_charges', namespace, action='write')
      write(iunit,'(a)') '# (Frequency-dependent) Born effective charge tensors'
      if (.not. write_real) write(iunit,'(a)') '# Real and imaginary parts'
      do iatom = 1, natoms
        write(iunit,'(a,i5,a,a5,a,f10.4)') 'Index: ', iatom, '   Label: ', trim(atom(iatom)%species%get_label()), &
          '   Ionic charge: ', charge(iatom)
 
        if (.not. write_real) write(iunit,'(a)') 'Real:'
        call output_tensor(real(this%charge(:, :, iatom), real64), dim, unit_one, iunit=iunit)
 
        if (.not. write_real) then
          write(iunit,'(a)') 'Imaginary:'
          call output_tensor(aimag(this%charge(:, :, iatom)), dim, unit_one, iunit=iunit)
        end if
 
        write(iunit,'(a)')
      end do
 
      if (.not. write_real) then
        write(iunit,'(a)') '# Magnitude and phase'
        do iatom = 1, natoms
          write(iunit,'(a,i5,a,a5,a,f10.4)') 'Index: ', iatom, '   Label: ', trim(atom(iatom)%species%get_label()), &
            '   Ionic charge: ', charge(iatom)
 
          write(iunit,'(a)') 'Magnitude:'
          call output_tensor(real(abs(this%charge(:, :, iatom)), real64), dim, unit_one, iunit=iunit)
 
          write(iunit,'(a)') 'Phase:'
 
          where (abs(this%charge(:, :, iatom)) > m_epsilon)
            phase = atan2(aimag(this%charge(:, :, iatom)), real(this%charge(:, :, iatom), real64))
          else where
            phase = m_zero
          end where
          call output_tensor(phase, dim, unit_one, write_average = .false., iunit=iunit)
          write(iunit,'(a)')
        end do
      end if
 
      write(iunit,'(a)') '# Discrepancy of Born effective charges from acoustic sum rule before correction, per atom'
      if (.not. write_real) write(iunit,'(a)') 'Real:'
      call output_tensor(real(this%delta, real64) , dim, unit_one, iunit=iunit)
      if (.not. write_real) then
        write(iunit,'(a)') 'Imaginary:'
        call output_tensor(aimag(this%delta), dim, unit_one, iunit=iunit)
      end if
 
      call io_close(iunit)
    end if
 
    pop_sub(born_output_charges)
  end subroutine born_output_charges
 
end module born_charges_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: