magnetic.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!!
!! 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 magnetic_oct_m
  use boundaries_oct_m
  use comm_oct_m
  use debug_oct_m
  use derivatives_oct_m
  use electron_space_oct_m
  use global_oct_m
  use grid_oct_m
  use ions_oct_m
  use kpoints_oct_m
  use mesh_function_oct_m
  use mesh_oct_m
  use messages_oct_m
  use mpi_oct_m
  use namespace_oct_m
  use poisson_oct_m
  use profiling_oct_m
  use species_oct_m
  use states_abst_oct_m
  use states_elec_oct_m
  use states_elec_dim_oct_m
  use submesh_oct_m
  use unit_oct_m
  use unit_system_oct_m
 
  implicit none
 
  private
  public ::                    &
    magnetic_density,          &
    magnetic_moment,           &
    write_magnetic_moments,    &
    magnetic_local_moments,    &
    magnetic_induced,          &
    magnetic_total_magnetization, &
    write_total_xc_torque     ,&
    calc_xc_torque
 
contains
 
  ! ---------------------------------------------------------
  subroutine magnetic_density(mesh, std, rho, md)
    class(mesh_t),           intent(in)  :: mesh
    type(states_elec_dim_t), intent(in)  :: std
    real(real64),            intent(in)  :: rho(:,:)
    real(real64),            intent(out) :: md(:,:)
 
    push_sub(magnetic_density)
 
    select case (std%ispin)
    case (unpolarized)
      md = m_zero
 
    case (spin_polarized)
      md = m_zero
      md(1:mesh%np, 3) = rho(1:mesh%np, 1) - rho(1:mesh%np, 2)
 
    case (spinors)
      md(1:mesh%np, 1) =  m_two*rho(1:mesh%np, 3)
      md(1:mesh%np, 2) = -m_two*rho(1:mesh%np, 4)
      md(1:mesh%np, 3) = rho(1:mesh%np, 1) - rho(1:mesh%np, 2)
    end select
 
    pop_sub(magnetic_density)
  end subroutine magnetic_density
 
 
  ! ---------------------------------------------------------
  subroutine magnetic_moment(mesh, st, rho, mm)
    class(mesh_t),       intent(in)  :: mesh
    type(states_elec_t), intent(in)  :: st
    real(real64),        intent(in)  :: rho(:,:)
    real(real64),        intent(out) :: mm(3)
 
    real(real64), allocatable :: md(:,:)
 
    push_sub(states_elec_magnetic_moment)
 
    safe_allocate(md(1:mesh%np, 1:3))
    call magnetic_density(mesh, st%d, rho, md)
 
    mm = dmf_integrate(mesh, 3, md)
 
    safe_deallocate_a(md)
 
    pop_sub(states_elec_magnetic_moment)
  end subroutine magnetic_moment
 
 
  ! ---------------------------------------------------------
  subroutine write_magnetic_moments(mesh, st, ions, boundaries, lmm_r, iunit, namespace)
    class(mesh_t),               intent(in) :: mesh
    type(states_elec_t),         intent(in) :: st
    type(ions_t),                intent(in) :: ions
    type(boundaries_t),          intent(in) :: boundaries
    real(real64),                intent(in) :: lmm_r
    integer,           optional, intent(in) :: iunit
    type(namespace_t), optional, intent(in) :: namespace
 
    integer :: ia
    real(real64) :: mm(max(mesh%box%dim, 3))
    real(real64), allocatable :: lmm(:,:)
 
    push_sub(write_magnetic_moments)
 
    call magnetic_moment(mesh, st, st%rho, mm)
    safe_allocate(lmm(1:max(mesh%box%dim, 3), 1:ions%natoms))
    call magnetic_local_moments(mesh, st, ions, boundaries, st%rho, lmm_r, lmm)
 
    message(1) = 'Total Magnetic Moment:'
    call messages_info(1, iunit=iunit, namespace=namespace)
    if (st%d%ispin == spin_polarized) then ! collinear spin
      write(message(1), '(a,f10.6)') ' mz = ', mm(3)
      call messages_info(1, iunit=iunit, namespace=namespace)
    else if (st%d%ispin == spinors) then ! non-collinear
      write(message(1), '(1x,3(a,f10.6,3x))') 'mx = ', mm(1),'my = ', mm(2),'mz = ', mm(3)
      call messages_info(1, iunit=iunit, namespace=namespace)
    end if
 
    write(message(1), '(a,a,a,f7.3,a)') 'Local Magnetic Moments (sphere radius [', &
      trim(units_abbrev(units_out%length)),'] = ', units_from_atomic(units_out%length, lmm_r), '):'
    call messages_info(1, iunit=iunit, namespace=namespace)
    if (st%d%ispin == spin_polarized) then ! collinear spin
      write(message(1),'(a,6x,14x,a)') ' Ion','mz'
      call messages_info(1, iunit=iunit, namespace=namespace)
      do ia = 1, ions%natoms
        write(message(1),'(i4,a10,f15.6)') ia, trim(ions%atom(ia)%species%get_label()), lmm(3, ia)
        call messages_info(1, iunit=iunit, namespace=namespace)
      end do
    else if (st%d%ispin == spinors) then ! non-collinear
      write(message(1),'(a,8x,13x,a,13x,a,13x,a)') ' Ion','mx','my','mz'
      call messages_info(1, iunit=iunit, namespace=namespace)
      do ia = 1, ions%natoms
        write(message(1),'(i4,a10,9f15.6)') ia, trim(ions%atom(ia)%species%get_label()), lmm(:, ia)
        call messages_info(1, iunit=iunit, namespace=namespace)
      end do
    end if
 
    safe_deallocate_a(lmm)
 
    pop_sub(write_magnetic_moments)
  end subroutine write_magnetic_moments
 
  ! ---------------------------------------------------------
  subroutine magnetic_local_moments(mesh, st, ions, boundaries, rho, rr, lmm)
    class(mesh_t),        intent(in)  :: mesh
    type(states_elec_t),  intent(in)  :: st
    type(ions_t),         intent(in)  :: ions
    type(boundaries_t),   intent(in)   :: boundaries
    real(real64),         intent(in)  :: rho(:,:)
    real(real64),         intent(in)  :: rr
    real(real64),         intent(out) :: lmm(max(mesh%box%dim, 3), ions%natoms)
 
    integer :: ia, idir, is
    real(real64), allocatable :: md(:, :)
    type(submesh_t) :: sphere
    real(real64) :: cosqr, sinqr
    complex(real64), allocatable :: phase_spiral(:)
 
    push_sub(magnetic_local_moments)
 
    safe_allocate(md(1:mesh%np, 1:max(mesh%box%dim, 3)))
 
    call magnetic_density(mesh, st%d, rho, md)
    lmm = m_zero
    do ia = 1, ions%natoms
      call submesh_init(sphere, ions%space, mesh, ions%latt, ions%pos(:, ia), rr)
 
      if (boundaries%spiral) then
        safe_allocate(phase_spiral(1:sphere%np))
        do is = 1, sphere%np
          phase_spiral(is) = exp(+m_zi*sum((sphere%rel_x(:,is) + sphere%center - mesh%x(sphere%map(is),:)) &
            *boundaries%spiral_q(1:mesh%box%dim)))
        end do
 
        if (mesh%box%dim>= 3) then
          lmm(1,ia) = m_zero
          lmm(2,ia) = m_zero
 
          do is = 1, sphere%np
            !There is a factor of 1/2 in phase_spiral
            cosqr = real(phase_spiral(is), real64)
            sinqr = aimag(phase_spiral(is))
            lmm(1,ia) = lmm(1,ia)+md(sphere%map(is),1)*cosqr - md(sphere%map(is),2)*sinqr
            lmm(2,ia) = lmm(2,ia)+md(sphere%map(is),1)*sinqr + md(sphere%map(is),2)*cosqr
          end do
          lmm(1,ia) = lmm(1,ia)*mesh%volume_element
          lmm(2,ia) = lmm(2,ia)*mesh%volume_element
          lmm(3,ia) = dsm_integrate_frommesh(mesh, sphere, md(1:mesh%np,3), reduce = .false.)
        else
          assert(.not. boundaries%spiral)
        end if
 
        safe_deallocate_a(phase_spiral)
      else
        do idir = 1, max(mesh%box%dim, 3)
          lmm(idir, ia) = dsm_integrate_frommesh(mesh, sphere, md(1:mesh%np,idir), reduce = .false.)
        end do
      end if
 
      call submesh_end(sphere)
    end do
 
    call mesh%allreduce(lmm)
 
    safe_deallocate_a(md)
 
    pop_sub(magnetic_local_moments)
  end subroutine magnetic_local_moments
 
  ! ---------------------------------------------------------
  subroutine magnetic_total_magnetization(mesh, st, qq, trans_mag)
    class(mesh_t),       intent(in)  :: mesh
    type(states_elec_t), intent(in)  :: st
    real(real64),        intent(in)  :: qq(:)
    complex(real64),     intent(out) :: trans_mag(6)
 
    integer :: ip
    complex(real64), allocatable :: tmp(:,:)
    real(real64), allocatable :: md(:, :)
    complex(real64) :: expqr
 
    push_sub(magnetic_total_magnetization)
 
    call profiling_in("TOTAL_MAGNETIZATION")
 
    safe_allocate(tmp(1:mesh%np, 1:6))
    safe_allocate(md(1:mesh%np, 1:max(mesh%box%dim, 3)))
 
    call magnetic_density(mesh, st%d, st%rho, md)
    do ip = 1, mesh%np
      expqr = exp(-m_zi*sum(mesh%x(ip, 1:mesh%box%dim)*qq(1:mesh%box%dim)))
      tmp(ip,1) = expqr*md(ip,1)
      tmp(ip,2) = expqr*md(ip,2)
      tmp(ip,3) = expqr*md(ip,3)
      tmp(ip,4) = conjg(expqr)*md(ip,1)
      tmp(ip,5) = conjg(expqr)*md(ip,2)
      tmp(ip,6) = conjg(expqr)*md(ip,3)
    end do
 
    trans_mag = zmf_integrate(mesh, 6, tmp)
 
    safe_deallocate_a(md)
    safe_deallocate_a(tmp)
 
    call profiling_out("TOTAL_MAGNETIZATION")
 
    pop_sub(magnetic_total_magnetization)
  end subroutine magnetic_total_magnetization
 
  ! ---------------------------------------------------------
  subroutine magnetic_induced(namespace, gr, st, psolver, kpoints, a_ind, b_ind)
    type(namespace_t),    intent(in)    :: namespace
    type(grid_t),         intent(in)    :: gr
    type(states_elec_t),  intent(inout) :: st
    type(poisson_t),      intent(in)    :: psolver
    type(kpoints_t),      intent(in)    :: kpoints
    real(real64), contiguous,    intent(out)   :: a_ind(:, :)
    real(real64), contiguous,    intent(out)   :: b_ind(:, :)
    
 
    integer :: idir
    real(real64), allocatable :: jj(:, :, :)
 
    push_sub(magnetic_induced)
 
    ! If the states are real, we should never have reached here, but
    ! just in case we return zero.
    if (states_are_real(st)) then
      a_ind = m_zero
      b_ind = m_zero
      pop_sub(magnetic_induced)
      return
    end if
 
    safe_allocate(jj(1:gr%np_part, 1:gr%der%dim, 1:st%d%nspin))
    call states_elec_calc_quantities(gr, st, kpoints, .false., paramagnetic_current = jj)
 
    !We sum the current for up and down, valid for collinear and noncollinear spins
    if (st%d%nspin > 1) then
      do idir = 1, gr%der%dim
        jj(:, idir, 1) = jj(:, idir, 1) + jj(:, idir, 2)
      end do
    end if
 
    a_ind = m_zero
    do idir = 1, gr%der%dim
      call dpoisson_solve(psolver, namespace, a_ind(:, idir), jj(:, idir, 1))
    end do
    ! This minus sign is introduced here because the current that has been used
    ! before is the "number-current density", and not the "charge-current density",
    ! and therefore there is a minus sign missing (electrons are negative charges...)
    a_ind(1:gr%np, 1:gr%der%dim) = - a_ind(1:gr%np, 1:gr%der%dim) / p_c
 
    call dderivatives_curl(gr%der, a_ind, b_ind)
 
    safe_deallocate_a(jj)
    pop_sub(magnetic_induced)
  end subroutine magnetic_induced
 
  subroutine write_total_xc_torque(iunit, mesh, vxc, st)
    integer,                  intent(in) :: iunit
    class(mesh_t),            intent(in) :: mesh
    real(real64),             intent(in) :: vxc(:,:)
    type(states_elec_t),      intent(in) :: st
 
    real(real64), allocatable :: torque(:,:)
    real(real64) :: tt(3)
 
    push_sub(write_total_xc_torque)
 
    safe_allocate(torque(1:mesh%np, 1:3))
 
    call calc_xc_torque(mesh, vxc, st, torque)
 
    tt = dmf_integrate(mesh, 3, torque)
 
    if (mpi_grp_is_root(mpi_world)) then
      write(iunit, '(a)') 'Total xc torque:'
      write(iunit, '(1x,3(a,es10.3,3x))') 'Tx = ', tt(1),'Ty = ', tt(2),'Tz = ', tt(3)
    end if
 
    safe_deallocate_a(torque)
 
    pop_sub(write_total_xc_torque)
  end subroutine write_total_xc_torque
 
  ! ---------------------------------------------------------
  subroutine calc_xc_torque(mesh, vxc, st, torque)
    class(mesh_t),            intent(in) :: mesh
    real(real64),             intent(in) :: vxc(:,:)
    type(states_elec_t),      intent(in) :: st
    real(real64),          intent(inout) :: torque(:,:)
 
    real(real64) :: mag(3), bxc(3)
    integer :: ip
 
    push_sub(calc_xc_torque)
 
    assert(st%d%ispin == spinors)
 
    do ip = 1, mesh%np
      mag(1) =  m_two * st%rho(ip, 3)
      mag(2) = -m_two * st%rho(ip, 4)
      mag(3) = st%rho(ip, 1) - st%rho(ip, 2)
      bxc(1) = -m_two * vxc(ip, 3)
      bxc(2) =  m_two * vxc(ip, 4)
      bxc(3) = -(vxc(ip, 1) - vxc(ip, 2))
      torque(ip, 1) = mag(2) * bxc(3) - mag(3) * bxc(2)
      torque(ip, 2) = mag(3) * bxc(1) - mag(1) * bxc(3)
      torque(ip, 3) = mag(1) * bxc(2) - mag(2) * bxc(1)
    end do
 
    pop_sub(calc_xc_torque)
  end subroutine calc_xc_torque
 
 
 
end module magnetic_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: