From 46dc8d62c6c99b483c00543b0feb057b227adb96 Mon Sep 17 00:00:00 2001
From: Nicolas Essing <n.essing@fz-juelich.de>
Date: Wed, 3 Nov 2021 16:58:46 +0100
Subject: [PATCH] Added the mean xc bfield as a scale for updating the
 constraint fields.

The mean xc bfield is calculated together with the torques as a quantity
that is needed to update the constraining fields. It is stored in the
bfield_data type.
The calculation is not a physically meaningfull thing as it ignores
the underlying mesh details, but good enough for this application.

This fixes the suprisingly small mixing needed for the constraint method
based only on the moments. As pointed out by Eduardo, comparing with the
definition in literature (e.g. Journal of Applied Physics 85, 4824 (1999)),
there was some magnetic field missing to get the update loop to the
right units (i.e. the right order of magnitude).
Tested only for bcc Fe so far, this constraint mode now works well
with a mixing factor of around 1 (0.9 is now default).
---
 source/KKRnano/source/bfield/bfield.f90                  | 9 +++++++--
 source/KKRnano/source/bfield/torque.f90                  | 6 +++++-
 source/KKRnano/source/datastructures/InputParamsNew.txt  | 2 +-
 source/KKRnano/source/datastructures/InputParams_mod.F90 | 4 ++--
 4 files changed, 15 insertions(+), 6 deletions(-)

diff --git a/source/KKRnano/source/bfield/bfield.f90 b/source/KKRnano/source/bfield/bfield.f90
index 24df00425..732d339cc 100644
--- a/source/KKRnano/source/bfield/bfield.f90
+++ b/source/KKRnano/source/bfield/bfield.f90
@@ -29,11 +29,16 @@ module mod_bfield
   !> Summary: A type storing information on magnetic fields for a single atom
   !-------------------------------------------------------------------------------
   type :: bfield_data
+    ! Noncollinear magnetic fields
     double precision, dimension(3) :: bfield_ext    !! external magnetic field in cartesian coordinates
     double precision, dimension(3) :: bfield_constr !! constraining field in cartesian coordinates
-    double precision, dimension(3) :: mag_torque    !! Magnetic torque
-    double precision, dimension(3) :: mag_mom       !! Magnetic moment
 
+    ! Torque and other information used to update the constraints fields
+    double precision, dimension(3) :: mag_torque     !! Magnetic torque
+    double precision, dimension(3) :: mag_mom        !! Magnetic moment
+    double precision               :: mean_xc_bfield !! Mean magnitude of xc bfield
+
+    ! Precalculated intermediate results
     double precision, dimension(:,:,:), allocatable :: thetallmat !! shapefun in the ll' expansion
   end type
 
diff --git a/source/KKRnano/source/bfield/torque.f90 b/source/KKRnano/source/bfield/torque.f90
index 660bf8ef6..d61174350 100644
--- a/source/KKRnano/source/bfield/torque.f90
+++ b/source/KKRnano/source/bfield/torque.f90
@@ -77,6 +77,9 @@ contains
     ! potential.
     bxc(:,:) = (vpot(:,:,1) - vpot(:,:,2)) / 2.
 
+    ! Calculate and save mean xc bfield strength
+    bfield%mean_xc_bfield = sum(abs(bxc(:,1))) / irmd
+
     ! Get magnetization density (in the global frame). For each angular momentum
     ! and radial index, rotate to the global frame.
     do ilm = 1, lmpotd
@@ -184,7 +187,8 @@ contains
     else if (constr_mode == 2) then
       old_b_constr = bfield%bfield_constr(:)
       bfield%bfield_constr(:) = old_b_constr - dot_product(old_b_constr,dir)*dir - &
-              (mag_mom_dir - dot_product(mag_mom_dir,dir)*dir)*constr_bfield_mixing
+              ( (mag_mom_dir - dot_product(mag_mom_dir,dir)*dir) * &
+                bfield%mean_xc_bfield * constr_bfield_mixing )
     else
       ! There might be other modes that are calculated somewhere else
       ! (e.g. mode 1, which only fixes the direction by not changing the local
diff --git a/source/KKRnano/source/datastructures/InputParamsNew.txt b/source/KKRnano/source/datastructures/InputParamsNew.txt
index c204040f0..270644bc1 100644
--- a/source/KKRnano/source/datastructures/InputParamsNew.txt
+++ b/source/KKRnano/source/datastructures/InputParamsNew.txt
@@ -133,7 +133,7 @@ l trans_bfield .FALSE.
 ### [bfields] apply magnetic field only in the muffin-tin (called lbfield_mt in KKRhost code)
 l mt_bfield .FALSE.
 ### [bfields] mixing parameter for the selfconsistency of the constraint magnetic fields
-d constr_bfield_mixing 0.03
+d constr_bfield_mixing 0.9
 ### [bfields] start applying magnetic field at iteration (called ibfield_itscf0 in KKRhost code)
 i itbfield0 0
 ### [bfields] stop applying magnetic field after iteration itscf1 (called ibfield_itscf1 in KKRhost code)
diff --git a/source/KKRnano/source/datastructures/InputParams_mod.F90 b/source/KKRnano/source/datastructures/InputParams_mod.F90
index 84abbf7d3..1cf06ab5e 100644
--- a/source/KKRnano/source/datastructures/InputParams_mod.F90
+++ b/source/KKRnano/source/datastructures/InputParams_mod.F90
@@ -682,9 +682,9 @@ integer function getValues(filename, self) result(ierror)
     destroy_and_return
   endif
 
-  ierror = getValue(cr, "constr_bfield_mixing", self%constr_bfield_mixing , def=0.03)
+  ierror = getValue(cr, "constr_bfield_mixing", self%constr_bfield_mixing , def=0.9)
   if (ierror == use_default) then
-    write(*,*) "WARNING: Bad/no value given for constr_bfield_mixing. Set constr_bfield_mixing to 0.03"
+    write(*,*) "WARNING: Bad/no value given for constr_bfield_mixing. Set constr_bfield_mixing to 0.9"
     ierror = 0 ! ok, no error
   elseif (ierror /= 0) then
     write(*,*) "Bad/no value given for constr_bfield_mixing."
-- 
GitLab