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 batch_ops_oct_m
  use boundaries_oct_m
  use comm_oct_m
  use debug_oct_m
  use derivatives_oct_m
  use electron_space_oct_m
  use epot_oct_m
  use global_oct_m
  use lalg_basic_oct_m
  use grid_oct_m
  use ions_oct_m
  use kpoints_oct_m
  use math_oct_m
  use mesh_function_oct_m
  use mesh_oct_m
  use messages_oct_m
  use mpi_oct_m
  use namespace_oct_m
  use phase_oct_m
  use poisson_oct_m
  use projector_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
  use wfs_elec_oct_m
 
  implicit none
 
  private
  public ::                    &
    magnetic_density,          &
    magnetic_moment,           &
    compute_and_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), contiguous, intent(in) :: rho(:,:)
    real(real64), contiguous, intent(out) :: md(:,:)
 
    integer :: ip
 
    push_sub(magnetic_density)
 
    select case (std%ispin)
    case (unpolarized)
      md = m_zero
 
    case (spin_polarized)
      md = m_zero
      !$omp parallel do
      do ip = 1, mesh%np
        md(ip, 3) = rho(ip, 1) - rho(ip, 2)
      end do
      !$omp end parallel do
 
    case (spinors)
      !$omp parallel do
      do ip = 1, mesh%np
        md(ip, 1) =  m_two*rho(ip, 3)
        md(ip, 2) = -m_two*rho(ip, 4)
        md(ip, 3) = rho(ip, 1) - rho(ip, 2)
      end do
      !$omp end parallel do
    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), contiguous, 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 compute_and_write_magnetic_moments(gr, st, phase, ep, ions, lmm_r, calc_orb_moments, iunit, namespace)
    type(grid_t),                intent(in) :: gr
    type(states_elec_t),         intent(in) :: st
    type(phase_t),               intent(in) :: phase
    type(epot_t),         intent(in)  :: ep
    type(ions_t),                intent(in) :: ions
    real(real64),                intent(in) :: lmm_r
    logical,           optional, intent(in) :: calc_orb_moments
    integer,           optional, intent(in) :: iunit
    type(namespace_t), optional, intent(in) :: namespace
 
    integer :: ia
    real(real64) :: mm(max(gr%box%dim, 3))
    real(real64), allocatable :: lmm(:,:)
 
    push_sub(compute_and_write_magnetic_moments)
 
    call magnetic_moment(gr, st, st%rho, mm)
    safe_allocate(lmm(1:max(gr%box%dim, 3), 1:ions%natoms))
    call magnetic_local_moments(gr, st, ions, gr%der%boundaries, st%rho, lmm_r, lmm)
 
    message(1) = 'Total Spin 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 Spin 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
 
    if (optional_default(calc_orb_moments, .false.)) then
      assert(st%d%ispin == spinors)
 
      if (st%system_grp%is_root()) write(iunit, '(1x)')
 
      call magnetic_orbital_moments(gr, st, phase, ep, ions, st%rho, lmm_r, lmm)
 
      write(message(1), '(a,a,a,f7.3,a)') 'Local Orbital 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)
      write(message(1),'(a,8x,13x,a,13x,a,13x,a)') ' Ion','Lx','Ly','Lz'
      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(compute_and_write_magnetic_moments)
  end subroutine compute_and_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), contiguous, 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, ww
    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))
            ! meshes with per-point volume elements carry the weight per point
            ww = mesh%volume_element
            if (mesh%use_vol_pp) ww = mesh%vol_pp(sphere%map(is))
            lmm(1,ia) = lmm(1,ia)+(md(sphere%map(is),1)*cosqr - md(sphere%map(is),2)*sinqr)*ww
            lmm(2,ia) = lmm(2,ia)+(md(sphere%map(is),1)*sinqr + md(sphere%map(is),2)*cosqr)*ww
          end do
          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_orbital_moments(gr, st, phase, ep, ions, rho, rr, lom)
    type(grid_t), intent(in) :: gr
    type(states_elec_t), intent(in) :: st
    type(phase_t), intent(in) :: phase
    type(epot_t), intent(in) :: ep
    type(ions_t), intent(in) :: ions
    real(real64), contiguous, intent(in) :: rho(:,:)
    real(real64), intent(in) :: rr
    real(real64), intent(out) :: lom(3, ions%natoms)
 
    integer :: ia, idir, is, map_ip, ibatch, idim, ist, ib, ik
    type(submesh_t) :: sphere
    type(wfs_elec_t) :: epsib
    type(wfs_elec_t), allocatable :: gpsib(:)
 
    complex(real64), allocatable :: psi(:,:), gf(:,:,:), gf_nl(:,:,:)
    complex(real64) :: acc(3), acc_nl(3)
    real(real64) :: x1, x2, x3, factor
 
    push_sub(magnetic_orbital_moments)
 
    assert(st%d%ispin == spinors)
    assert(ions%space%dim == 3)
    if (gr%use_vol_pp) then
      call messages_not_implemented("Orbital magnetic moments with per-point volume elements")
    end if
 
    lom = m_zero
 
    safe_allocate(psi(1:gr%np_part, 1:st%d%dim))
    safe_allocate(gf(1:gr%np, 1:st%d%dim, 1:gr%der%dim))
    safe_allocate(gf_nl(1:gr%np, 1:st%d%dim, 1:gr%der%dim))
    safe_allocate_type_array(wfs_elec_t, gpsib, (1:gr%der%dim))
 
    do ik = st%d%kpt%start, st%d%kpt%end
      do ib = st%group%block_start, st%group%block_end
        ! copy st%group%psib(ib, ik) to epsib and set the phase
        call phase%copy_and_set_phase(gr, st%d%kpt, st%group%psib(ib, ik), epsib)
 
        ! this now takes non-orthogonal axis into account
        call zderivatives_batch_grad(gr%der, epsib, gpsib, set_bc=.false.)
 
        do ibatch = 1, epsib%nst
          ist = st%group%psib(ib, ik)%ist(ibatch)
          factor = st%kweights(ik)*st%occ(ist, ik)
          if (abs(factor) < m_epsilon) cycle
 
          ! Note: this is not efficient at all, but we do not have ready-to-use batch functions for this yet
          ! The regeneration of many submeshes is also very inefficient
          call batch_get_state(epsib, ibatch, gr%np, psi)
          call batch_get_state(gpsib(1), ibatch, gr%np, gf(:,:,1))
          call batch_get_state(gpsib(2), ibatch, gr%np, gf(:,:,2))
          call batch_get_state(gpsib(3), ibatch, gr%np, gf(:,:,3))
 
          ! TODO: precompute \psi^* (-i\nabla) \psi as Im[\psi^* (\nabla) \psi]
 
          do ia = 1, ions%natoms
            call submesh_init(sphere, ions%space, gr, ions%latt, ions%pos(:, ia), rr)
 
            gf_nl = m_zero
            call zprojector_commute_r(ep%proj(ia), gr, gr%der%boundaries, st%d%dim, 1, ik, psi, gf_nl(:, :, 1))
            call zprojector_commute_r(ep%proj(ia), gr, gr%der%boundaries, st%d%dim, 2, ik, psi, gf_nl(:, :, 2))
            call zprojector_commute_r(ep%proj(ia), gr, gr%der%boundaries, st%d%dim, 3, ik, psi, gf_nl(:, :, 3))
 
            acc = m_zero
            acc_nl = m_zero
            do idim = 1, st%d%dim
              do is = 1, sphere%np
                map_ip = sphere%map(is)
                x1 = sphere%rel_x(1, is)
                x2 = sphere%rel_x(2, is)
                x3 = sphere%rel_x(3, is)
 
                acc(1) = acc(1) + conjg(psi(map_ip, idim)) * (x2 * gf(map_ip, idim, 3) - x3 * gf(map_ip, idim, 2))
                acc(2) = acc(2) + conjg(psi(map_ip, idim)) * (x3 * gf(map_ip, idim, 1) - x1 * gf(map_ip, idim, 3))
                acc(3) = acc(3) + conjg(psi(map_ip, idim)) * (x1 * gf(map_ip, idim, 2) - x2 * gf(map_ip, idim, 1))
 
                acc_nl(:) = acc_nl(:) + conjg(psi(map_ip, idim)) * gf_nl(map_ip, idim, :)
              end do
            end do
            call submesh_end(sphere)
 
            ! The -i is absorbed in the imaginary part
            lom(:, ia) = lom(:, ia) + factor * gr%volume_element * (aimag(acc) + dcross_product(ions%pos(:, ia), aimag(acc_nl)))
          end do
        end do
 
        do idir = 1, gr%der%dim
          call gpsib(idir)%end()
        end do
 
        call epsib%end()
      end do
    end do
 
    call gr%allreduce(lom)
    if (st%parallel_in_states .or. st%d%kpt%parallel) then
      call comm_allreduce(st%st_kpt_mpi_grp, lom)
    end if
 
    safe_deallocate_a(psi)
    safe_deallocate_a(gf)
    safe_deallocate_a(gf_nl)
    safe_deallocate_a(gpsib)
 
    pop_sub(magnetic_orbital_moments)
  end subroutine magnetic_orbital_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), contiguous, 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(1:mesh%box%dim, ip)*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...)
    call lalg_scal(gr%np, gr%der%dim, -m_one / p_c, a_ind)
 
    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), contiguous, 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 (st%system_grp%is_root()) 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), contiguous, intent(in) :: vxc(:,:)
    type(states_elec_t), intent(in) :: st
    real(real64), contiguous, intent(out) :: 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: