magnetic_constrain.F90

Go to the documentation of this file.

!! Copyright (C) 2021 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 magnetic_constrain_oct_m
  use debug_oct_m
  use electron_space_oct_m
  use global_oct_m
  use, intrinsic :: iso_fortran_env
  use lattice_vectors_oct_m
  use magnetic_oct_m
  use messages_oct_m
  use mesh_oct_m
  use namespace_oct_m
  use parser_oct_m
  use profiling_oct_m
  use space_oct_m
  use states_elec_dim_oct_m
  use submesh_oct_m
 
  implicit none
 
  private
  public ::                          &
    magnetic_constrain_t,            &
    magnetic_constrain_init,         &
    magnetic_constrain_update,       &
    magnetic_constrain_end
 
  integer, public, parameter ::        &
    CONSTRAIN_NONE                = 0, &
    constrain_dir                 = 1, &
    constrain_full                = 2
 
 
  type magnetic_constrain_t
    private
    integer, public           :: level  = constrain_none 
    real(real64)              :: lambda = m_zero         
    real(real64), allocatable :: constrain(:,:)
    real(real64), allocatable, public :: pot(:,:)
    real(real64), public      :: energy
 
    real(real64)       :: lmm_r
  end type magnetic_constrain_t
 
  ! Threshold for considering zero magnetic moments
  real(real64), parameter :: tol_mag_norm = 1.0e-6_real64
 
contains
 
  ! ---------------------------------------------------------
  subroutine magnetic_constrain_init(this, namespace, mesh, std, natoms, min_dist)
    type(magnetic_constrain_t), intent(inout) :: this
    type(namespace_t),         intent(in)    :: namespace
    class(mesh_t),             intent(in)    :: mesh
    type(states_elec_dim_t),   intent(in)    :: std
    integer,                   intent(in)    :: natoms
    real(real64),              intent(in)    :: min_dist
 
    integer :: ia, idir
    real(real64)   :: rmin, norm
    type(block_t) :: blk
 
    push_sub(magnetic_constrain_init)
 
    !%Variable MagneticConstrain
    !%Type integer
    !%Default no
    !%Section Hamiltonian
    !%Description
    !% This variable selects which magnetic constrain expression is added to the Hamiltonian.
    !% The <tt>AtomsMagnetDirection</tt> block is used for determining the constrained direction.
    !%
    !% This is ignored for independent particles, as in this case, there is nothing that
    !% determines the magnetic direction.
    !%
    !%Option constrain_none 0
    !% No constrain is added to the Hamiltonian.
    !%Option constrain_dir 1
    !% We are adding a constrain for the direction and sign of the magnetic moments only,
    !% see Ma and Dudarev, PRB 91, 054420 (2015).
    !%Option constrain_full 2
    !% We are adding a constrain for the direction and norm of the magnetic moments only,
    !% see Ma and Dudarev, PRB 91, 054420 (2015).
    !%End
    call parse_variable(namespace, 'MagneticConstrain', constrain_none, this%level)
    call messages_print_var_value('MagneticConstrain', this%level, namespace=namespace)
    if (this%level == constrain_none) then
      pop_sub(magnetic_constrain_init)
      return
    end if
 
    call messages_experimental('MagneticConstrain', namespace=namespace)
 
    !%Variable MagneticConstrainStrength
    !%Type float
    !%Default 0.01
    !%Section Hamiltonian
    !%Description
    !% This variable determines the value of the Lagrange multiplier used for the constrain term.
    !%End
    call parse_variable(namespace, 'MagneticConstrainStrength', 0.01_real64, this%lambda)
    call messages_print_var_value('MagneticConstrainStrength', this%lambda, namespace=namespace)
 
 
    !We are using the AtomsMagnetDirection block for the contrain
    if (parse_block(namespace, 'AtomsMagnetDirection', blk) < 0) then
      message(1) = "AtomsMagnetDirection block is not defined."
      message(2) = "Magnetic constrained DFT cannot be used without it."
      call messages_fatal(2, namespace=namespace)
    end if
 
    if (parse_block_n(blk) /= natoms) then
      message(1) = "AtomsMagnetDirection block has the wrong number of rows."
      call messages_fatal(1, namespace=namespace)
    end if
 
    if(std%ispin == unpolarized) then
      message(1) = "Magnetic constrains can only be used for spin-polized and spinor calculations."
      call messages_fatal(1, namespace=namespace)
    end if
 
    safe_allocate(this%constrain(1:3, natoms))
    do ia = 1, natoms
      !Read from AtomsMagnetDirection block
      if (std%ispin == spin_polarized) then
        call parse_block_float(blk, ia-1, 0, this%constrain(3, ia))
        this%constrain(1:2, ia) = m_zero
      elseif (std%ispin == spinors) then
        do idir = 1, 3
          call parse_block_float(blk, ia-1, idir-1, this%constrain(idir, ia))
          if (abs(this%constrain(idir, ia)) < tol_mag_norm) this%constrain(idir, ia) = m_zero
        end do
      end if
 
      if(this%level == constrain_dir) then
        !Renormalization for constrained directions
        norm = norm2(this%constrain(1:3, ia))
        this%constrain(1:3, ia) = this%constrain(1:3, ia) / norm
      end if
    end do
    call parse_block_end(blk)
 
 
    safe_allocate(this%pot(1:mesh%np, 1:std%nspin))
    this%pot = m_zero
 
    rmin = min_dist
    call parse_variable(namespace, 'LocalMagneticMomentsSphereRadius', min(m_half*rmin, 100.0_real64), &
      this%lmm_r)
 
    pop_sub(magnetic_constrain_init)
  end subroutine magnetic_constrain_init
 
  ! ---------------------------------------------------------
  subroutine magnetic_constrain_end(this)
    type(magnetic_constrain_t), intent(inout) :: this
 
    push_sub(magnetic_constrain_end)
 
    safe_deallocate_a(this%constrain)
    safe_deallocate_a(this%pot)
 
    pop_sub(magnetic_constrain_end)
  end subroutine magnetic_constrain_end
 
  ! ---------------------------------------------------------
  subroutine magnetic_constrain_update(this, mesh, std, space, latt, pos, rho)
    type(magnetic_constrain_t), intent(inout) :: this
    class(mesh_t),              intent(in)    :: mesh
    type(states_elec_dim_t),    intent(in)    :: std
    class(space_t),             intent(in)    :: space
    type(lattice_vectors_t),    intent(in)    :: latt
    real(real64),               intent(in)    :: pos(:,:)
    real(real64),               intent(in)    :: rho(:,:)
 
    integer :: ia, idir, ip
    real(real64) :: bb(3), b2, lmm(3), dotp, norm, xx, max_r
    real(real64), allocatable :: md(:,:), mdf(:), mask(:)
    type(submesh_t) :: sphere
    integer :: natoms
    ! Threshold for computes sinc(x) with a Taylor expansion
    ! The error in the Taylor expansion is below 8e-11 for this value
    real(real64), parameter :: tol_sinc = 0.01_real64
 
    if (this%level == constrain_none) return
 
    push_sub(magnetic_constrain_update)
 
    call profiling_in(tostring(magnetic_constrain))
 
    this%pot = m_zero
    this%energy = m_zero
 
    safe_allocate(md(1:mesh%np, 1:max(space%dim, 3)))
    call magnetic_density(mesh, std, rho, md)
 
    natoms = size(pos, dim=2)
 
    do ia = 1, natoms
      call submesh_init(sphere, space, mesh, latt, pos(:, ia), this%lmm_r)
 
      ! Create a mask to be applied to the magnetic moment, see
      ! Ma and Dudarev, PRB 91, 054420 (2015)
      ! This helps the convergence in avoiding jitter at the edge of the sphere,
      ! and produces a continuous potential
      max_r = maxval(sphere%r)
      safe_allocate(mask(1:sphere%np))
      mask = m_zero
      !$omp parallel do private(xx)
      do ip = 1, sphere%np
        xx = sphere%r(ip) * m_pi / max_r
        if(xx < tol_sinc) then
          mask(ip) = m_one-xx*xx/6.0_real64
        else
          mask(ip) = sin(xx)/xx
        end if
      end do
 
      ! We compute the local moment here
      ! We multiply here by a function that decreases monotonically to zero
      ! at the boundary of the atomic sphere.
      lmm = m_zero
      safe_allocate(mdf(1:sphere%np))
      do idir = 1, max(space%dim, 3)
        !$omp parallel do
        do ip = 1, sphere%np
          mdf(ip) = md(sphere%map(ip), idir) * mask(ip)
        end do
        lmm(idir) = dsm_integrate(mesh, sphere, mdf)
      end do
      safe_deallocate_a(mdf)
 
      ! See Ma and Dudarev, PRB 91, 054420 (2015)
      ! See for instance https://www.vasp.at/wiki/index.php/I_CONSTRAINED_M
      select case(this%level)
      case(constrain_dir)
        norm = norm2(lmm)
        if (norm > tol_mag_norm) then
          dotp = dot_product(lmm, this%constrain(1:3, ia))
          this%energy = this%energy + this%lambda * (norm - dotp)
          bb = lmm/norm - this%constrain(:,ia)
        else
          cycle
        end if
      case(constrain_full)
        this%energy = this%energy + this%lambda * sum((lmm - this%constrain(:,ia))**2)
        bb = m_two * (lmm - this%constrain(:,ia))
      end select
      bb = bb * this%lambda
 
      ! We are adding an effective Zeeman term within the atomic sphere
      select case(std%ispin)
      case (spin_polarized)
        b2 = norm2(bb(1:max(space%dim, 3)))
        do ip = 1, sphere%np
          this%pot(sphere%map(ip), 1) = this%pot(sphere%map(ip), 1) + b2 * mask(ip)
          this%pot(sphere%map(ip), 2) = this%pot(sphere%map(ip), 2) - b2 * mask(ip)
        end do
 
      case (spinors)
        do ip = 1, sphere%np
          this%pot(sphere%map(ip), 1) = this%pot(sphere%map(ip), 1) + bb(3) * mask(ip)
          this%pot(sphere%map(ip), 2) = this%pot(sphere%map(ip), 2) - bb(3) * mask(ip)
          this%pot(sphere%map(ip), 3) = this%pot(sphere%map(ip), 3) + bb(1) * mask(ip)
          this%pot(sphere%map(ip), 4) = this%pot(sphere%map(ip), 4) - bb(2) * mask(ip)
        end do
 
      end select
 
      safe_deallocate_a(mask)
      call submesh_end(sphere)
    end do
 
    safe_deallocate_a(md)
    safe_deallocate_a(mdf)
 
    call profiling_out(tostring(magnetic_constrain))
 
    pop_sub(magnetic_constrain_update)
  end subroutine magnetic_constrain_update
 
end module magnetic_constrain_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: