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

@@ -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
 

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

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)

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."