From 729fb8296c3b57ca503cabb02c4f1b4a61831930 Mon Sep 17 00:00:00 2001
From: Nicolas Essing <n.essing@fz-juelich.de>
Date: Mon, 11 Oct 2021 16:06:00 +0200
Subject: [PATCH] Added bfield to potential where it should be needed.

Calls to add_bfield() whereever the spin-orbit hamiltonian is added to the potential.
Adjusted arguments to the surrounding subroutines and dependencies in the Makefile.

I had to move two subroutines from NonCollinearMagnetism_mod.F90 to a separate
module in NonCollinearMagnetism_Helpers_mod.F90 to avoid circular dependencies.
(They did not change apart from that move.)
---
 source/KKRnano/source/Makefile                |   5 +-
 .../NonCollinearMagnetism_Helpers_mod.F90     | 129 +++++++++++++
 .../source/NonCollinearMagnetism_mod.F90      | 169 +++++-------------
 .../KKRnano/source/ProcessKKRresults_mod.F90  |   2 +-
 .../source/ScatteringCalculation_mod.F90      |   7 +-
 source/KKRnano/source/bfield/bfield.f90       |   2 +-
 source/KKRnano/source/wrappers_mod.F90        |  14 +-
 7 files changed, 193 insertions(+), 135 deletions(-)
 create mode 100644 source/KKRnano/source/NonCollinearMagnetism_Helpers_mod.F90

diff --git a/source/KKRnano/source/Makefile b/source/KKRnano/source/Makefile
index dd88823bf..20f27ba95 100644
--- a/source/KKRnano/source/Makefile
+++ b/source/KKRnano/source/Makefile
@@ -504,7 +504,7 @@ ConfigReader_mod.o: ConfigReaderDictionary_mod.o
 KKROperator_mod.o: Truncation_mod.o ClusterInfo_mod.o fillKKRMatrix_mod.o SparseMatrixDescription_mod.o bsrmm_mod.o
 Lloyd0_new_mod.o: BasisAtom_mod.o RadialMeshData_mod.o ShapefunData_mod.o EnergyMesh_mod.o LDAUData_mod.o GauntCoefficients_mod.o ValenceDensity_mod.o 
 ClusterInfo_mod.o: Statistics_mod.o TruncationZone_mod.o RefCluster_mod.o two_sided_comm_mod.o ExchangeTable_mod.o
-ScatteringCalculation_mod.o: KKRnanoParallel_mod.o Constants_mod.o SingleSiteRef_mod.o arraytest2_mod.o Logging_mod.o InputParams_mod.o TruncationZone_mod.o wrappers_mod.o Main2Arrays_mod.o KKRnano_Comm_mod.o jij_calc_mod.o BCPOperator_mod.o LDAUData_mod.o one_sided_comm_mod.o GauntCoefficients_mod.o InitialGuess_mod.o CalculationData_mod.o EBalanceHandler_mod.o EnergyMesh_mod.o KKROperator_mod.o BasisAtom_mod.o KKRresults_mod.o IterativeSolver_mod.o RefCluster_mod.o TimerMpi_mod.o DimParams_mod.o kloopz1_mod.o ClusterInfo_mod.o JijData_mod.o Truncation_mod.o two_sided_comm_mod.o ChebMeshData_mod.o NonCollinearMagnetism_mod.o
+ScatteringCalculation_mod.o: KKRnanoParallel_mod.o Constants_mod.o SingleSiteRef_mod.o arraytest2_mod.o Logging_mod.o InputParams_mod.o TruncationZone_mod.o wrappers_mod.o Main2Arrays_mod.o KKRnano_Comm_mod.o jij_calc_mod.o BCPOperator_mod.o LDAUData_mod.o one_sided_comm_mod.o GauntCoefficients_mod.o InitialGuess_mod.o CalculationData_mod.o EBalanceHandler_mod.o EnergyMesh_mod.o KKROperator_mod.o BasisAtom_mod.o KKRresults_mod.o IterativeSolver_mod.o RefCluster_mod.o TimerMpi_mod.o DimParams_mod.o kloopz1_mod.o ClusterInfo_mod.o JijData_mod.o Truncation_mod.o two_sided_comm_mod.o ChebMeshData_mod.o NonCollinearMagnetism_mod.o NonCollinearMagnetism_Helpers_mod.o
 total_energy_mod.o: Quadrature_mod.o BasisAtom_mod.o RadialMeshData_mod.o ShapefunData_mod.o ShapeGauntCoefficients_mod.o
 vxcgga.o: XCFunctionals_mod.o
 KKRzero_mod.o: Exceptions_mod.o BrillouinZoneMesh_mod.o PositionReader_mod.o MadelungCalculator_mod.o Startb1_mod.o EnergyMesh_mod.o EnergyMeshHelpers_mod.o Lattice_mod.o BrillouinZone_mod.o ConfigReader_mod.o ConfigReaderDictionary_mod.o InputParams_mod.o Main2Arrays_mod.o DimParams_mod.o BasisAtom_mod.o ShapefunData_mod.o PotentialData_mod.o AtomicCoreData_mod.o RadialMeshData_mod.o ShapefunData_mod.o read_formatted_mod.o read_formatted_shapefun_mod.o ldauinfo_read.o ldaustart.o
@@ -561,8 +561,9 @@ bsrmm_mod.o: CacheOverlap_mod.o
 CacheOverlap_mod.o:
 ChebMeshData_mod.o: InputParams_mod.o RadialMeshData_mod.o
 Truncation_mod.o: Logging_mod.o Exceptions_mod.o TruncationZone_mod.o
-NonCollinearMagnetism_mod.o: RadialMeshData_mod.o ChebMeshData_mod.o read_formatted_shapefun_mod.o
+NonCollinearMagnetism_mod.o: RadialMeshData_mod.o ChebMeshData_mod.o read_formatted_shapefun_mod.o NonCollinearMagnetism_Helpers_mod.o bfield.o
 NonCollinearMagnetismData_mod.o: Exceptions_mod.o
 vintras_new.o: SingleSiteHelpers_mod.o
+bfield.o: NonCollinearMagnetism_Helpers_mod.o
 
 # DO NOT DELETE
diff --git a/source/KKRnano/source/NonCollinearMagnetism_Helpers_mod.F90 b/source/KKRnano/source/NonCollinearMagnetism_Helpers_mod.F90
new file mode 100644
index 000000000..79b183340
--- /dev/null
+++ b/source/KKRnano/source/NonCollinearMagnetism_Helpers_mod.F90
@@ -0,0 +1,129 @@
+module NonCollinearMagnetism_Helpers_mod
+
+  private
+  public :: rotatematrix, create_Umatrix
+  
+  contains
+  
+  subroutine rotatematrix(mat,theta,phi,lmmax,mode)
+  ! rotates a matrix in the local frame pointing in
+  ! the direction of phi and theta to the global frame
+  implicit none
+  !interface
+  double complex,intent(inout)    ::  mat(2*lmmax,2*lmmax)
+  double precision,intent(in)     :: phi
+  double precision,intent(in)     :: theta
+  integer                         :: lmmax
+  integer                         :: mode
+  !local
+  double complex   :: Umat(2*lmmax,2*lmmax)
+  double complex   :: Udeggamat(2*lmmax,2*lmmax)
+  double complex   :: mattemp(2*lmmax,2*lmmax)
+  !double precision :: matmat_zmzm
+
+  !***********************************************************************
+  ! create the rotation matrix:
+  !     | cos(theta/2) exp(-i/2 phi)   -sin(theta/2) exp(-i/2 phi) |
+  !  U= |                                                          |
+  !     | sin(theta/2) exp( i/2 phi)    cos(theta/2) exp( i/2 phi) |
+  !
+  !  Udegga = transpose(complex conjug ( U ) )
+  !***********************************************************************
+
+
+  call create_Umatrix(theta,phi,lmmax,Umat,Udeggamat)
+  !***********************************************************************
+  ! calculate matrix in the global frame:
+  !
+  !  t_glob = U * t_loc * Udegga
+  !***********************************************************************
+
+
+  if (mode==0) then ! 'loc->glob'
+    call zgemm('N','N',2*lmmax,2*lmmax,2*lmmax,(1d0,0d0),mat,2*lmmax,Udeggamat,2*lmmax,(0d0,0d0),mattemp,2*lmmax)
+    call zgemm('N','N',2*lmmax,2*lmmax,2*lmmax,(1d0,0d0),Umat,2*lmmax,mattemp,2*lmmax,(0d0,0d0),mat,2*lmmax)
+  elseif (mode==1) then !'glob->loc'
+    call zgemm('N','N',2*lmmax,2*lmmax,2*lmmax,(1d0,0d0),mat,2*lmmax,Umat,2*lmmax,(0d0,0d0),mattemp,2*lmmax)
+    call zgemm('N','N',2*lmmax,2*lmmax,2*lmmax,(1d0,0d0),Udeggamat,2*lmmax,mattemp,2*lmmax,(0d0,0d0),mat,2*lmmax)
+  else
+    stop '[rotatematrix] mode not known'
+  end if
+  !  writE(324,'(5000F)') tmat
+  ! stop
+
+  end subroutine rotatematrix
+
+  subroutine create_Umatrix(theta,phi,lmmax,Umat,Udeggamat)
+  implicit none
+  !***********************************************************************
+  ! create the rotation matrix:
+  !     | cos(theta/2) exp(-i/2 phi)   -sin(theta/2) exp(-i/2 phi) |
+  !  U= |                                                          |
+  !     | sin(theta/2) exp( i/2 phi)    cos(theta/2) exp( i/2 phi) |
+  !
+  !  Udegga = transpose(complex conjug ( U ) )
+  !***********************************************************************double
+  !precision :: phi
+  !interface
+  double precision,intent(in)     :: phi 
+  double precision,intent(in)     :: theta
+  integer,intent(in)              :: lmmax
+  double complex,intent(out)      :: Umat(2*lmmax,2*lmmax)
+  double complex,intent(out)      :: Udeggamat(2*lmmax,2*lmmax)
+  !local
+  double complex                  :: Umat11,Umat12,Umat21,Umat22
+  double complex                  :: Udeggamat11,Udeggamat12,Udeggamat21,Udeggamat22
+  integer                         :: ival
+  double complex,parameter        :: ci=(0.0D0,1.0D0)
+  character*25               :: spinmode
+
+  spinmode='kkr'
+  if (spinmode=='regular') then
+    Umat11      =  cos(theta/2.0D0)*exp(-ci/2.0D0*phi)
+    Umat12      = -sin(theta/2.0D0)*exp(-ci/2.0D0*phi)
+    Umat21      =  sin(theta/2.0D0)*exp( ci/2.0D0*phi)
+    Umat22      =  cos(theta/2.0D0)*exp( ci/2.0D0*phi)
+  else if (spinmode=='kkr') then
+    Umat11      =  cos(theta/2.0D0)*exp( ci/2.0D0*phi)
+    Umat12      =  sin(theta/2.0D0)*exp( ci/2.0D0*phi)
+    Umat21      = -sin(theta/2.0D0)*exp(-ci/2.0D0*phi)
+    Umat22      =  cos(theta/2.0D0)*exp(-ci/2.0D0*phi)
+  else 
+    stop '[create_Umatrix] mode not known'
+  end if
+
+  Umat=(0.0D0,0.0D0)
+  do ival=1,lmmax
+    Umat(      ival,      ival) = Umat11
+    Umat(      ival,lmmax+ival) = Umat12
+    Umat(lmmax+ival,ival)       = Umat21
+    Umat(lmmax+ival,lmmax+ival) = Umat22
+  end do
+
+  if (spinmode=='regular') then
+  Udeggamat11 =  cos(theta/2.0D0)*exp( ci/2.0D0*phi)
+  Udeggamat12 =  sin(theta/2.0D0)*exp(-ci/2.0D0*phi)
+  Udeggamat21 = -sin(theta/2.0D0)*exp( ci/2.0D0*phi)
+  Udeggamat22 =  cos(theta/2.0D0)*exp(-ci/2.0D0*phi)
+  else if (spinmode=='kkr') then
+  Udeggamat11 =  cos(theta/2.0D0)*exp(-ci/2.0D0*phi)
+  Udeggamat12 = -sin(theta/2.0D0)*exp( ci/2.0D0*phi)
+  Udeggamat21 =  sin(theta/2.0D0)*exp(-ci/2.0D0*phi)
+  Udeggamat22 =  cos(theta/2.0D0)*exp( ci/2.0D0*phi)
+  else 
+    stop '[create_Umatrix] mode not known'
+  end if
+
+
+
+  Udeggamat=(0.0D0,0.0D0)
+  do ival=1,lmmax
+    Udeggamat(      ival,      ival) = Udeggamat11
+    Udeggamat(      ival,lmmax+ival) = Udeggamat12
+    Udeggamat(lmmax+ival,ival)       = Udeggamat21
+    Udeggamat(lmmax+ival,lmmax+ival) = Udeggamat22
+  end do
+
+  end subroutine create_Umatrix
+
+end module
diff --git a/source/KKRnano/source/NonCollinearMagnetism_mod.F90 b/source/KKRnano/source/NonCollinearMagnetism_mod.F90
index 8fc63b771..d4065f141 100644
--- a/source/KKRnano/source/NonCollinearMagnetism_mod.F90
+++ b/source/KKRnano/source/NonCollinearMagnetism_mod.F90
@@ -13,6 +13,9 @@ module NonCollinearMagnetism_mod
 !-------------------------------------------------------------------------------
 use RadialMeshData_mod!, only:
 use ChebMeshData_mod!, only
+use NonCollinearMagnetism_Helpers_mod, only: rotatematrix, create_Umatrix
+use mod_bfield, only: bfield_data, add_bfield
+
 implicit none
 private
 
@@ -27,7 +30,9 @@ SUBROUTINE tmat_newsolver(ie,nspin,lmax,zat,socscale,  &
         rpan_intervall,ipan_intervall,  &
         rnew,vinsnew,theta,phi,ipot,  &
        ! lly,        &
-        lmpotd,irmd_new,TmatN,soc,enable_quad_prec) ! new input parameters
+        lmpotd,irmd_new,TmatN,soc,enable_quad_prec, &
+        bfield, imt, iteration_number, itscf0, itscf1, &
+        lbfield, lbfield_constr, lbfield_trans, lbfield_mt) 
  
 ! Code converted using TO_F90 by Alan Miller
 ! Date: 2016-04-18  Time: 14:58:02
@@ -58,6 +63,13 @@ INTEGER, INTENT(IN)                      :: irmd_new
 DOUBLE COMPLEX, INTENT(OUT)              :: TmatN(:,:)
 LOGICAL, INTENT(IN)                      :: soc
 LOGICAL, INTENT(IN)                      :: enable_quad_prec
+type(bfield_data), intent(in)            :: bfield !! Information on the noncollinear magnetic fields
+integer, intent(in)                      :: imt !! MT radius (index in cheb mesh)
+integer, intent(in)                      :: iteration_number !! Current iteration number
+integer, intent(in)                      :: itscf0, itscf1 !! Window of iterations when to apply fields
+logical, intent(in)                      :: lbfield, lbfield_constr !! Wheter to apply the fields
+logical, intent(in)                      :: lbfield_trans ! Apply only transveral
+logical, intent(in)                      :: lbfield_mt ! Apply only up do MT radius
 
 INTEGER :: lmmaxd
 INTEGER :: lmmaxso
@@ -152,6 +164,11 @@ allocate(ull(nsra*lmmaxso,lmmaxso,irmdnew))
       vnspll0,vnspll1,'1',soc)
 
   ! add bfield
+  if (lbfield) then
+    call add_bfield(bfield, vnspll1, theta, phi, imt, iteration_number, &
+                    itscf0, itscf1, lbfield_constr, lbfield_trans, &
+                    lbfield_mt, .false.)
+  end if
 
 !c extend matrix for the SRA treatment
   vnspll=czero
@@ -309,6 +326,8 @@ deallocate(rll)
 deallocate(ull)
 
 END SUBROUTINE tmat_newsolver
+
+
 SUBROUTINE rhovalnew(ldorhoef,ielast,nsra,nspin,lmax,ez,wez,zat,  &
         socscale,cleb,icleb,iend,ifunm,lmsp,ncheb,  &
         npan_tot,npan_log,npan_eq,rmesh,irws,  &
@@ -317,7 +336,9 @@ SUBROUTINE rhovalnew(ldorhoef,ielast,nsra,nspin,lmax,ez,wez,zat,  &
         moment_x,moment_y,moment_z, &
         ipot,  &
         den_out,espv,rho2ns,r2nef,gmatn, muorb,  &
-        lpotd,lmaxd,irmd,irmd_new,iemxd,soc,enable_quad_prec) ! new parameters
+        lpotd,lmaxd,irmd,irmd_new,iemxd,soc,enable_quad_prec, &
+        bfield, imt, iteration_number, itscf0, itscf1, &
+        lbfield, lbfield_constr, lbfield_trans, lbfield_mt)
  
 ! Code converted using TO_F90 by Alan Miller
 ! Date: 2016-04-21  Time: 11:39:57
@@ -369,6 +390,14 @@ INTEGER, INTENT(IN)                      :: irmd_new
 INTEGER, INTENT(IN)                      :: iemxd
 LOGICAL, INTENT(IN)                      :: soc
 LOGICAL, INTENT(IN)                      :: enable_quad_prec
+type(bfield_data), intent(in)            :: bfield !! Information on the noncollinear magnetic fields
+integer, intent(in)                      :: imt !! MT radius (index in cheb mesh)
+integer, intent(in)                      :: iteration_number !! Current iteration number
+integer, intent(in)                      :: itscf0, itscf1 !! Window of iterations when to apply fields
+logical, intent(in)                      :: lbfield, lbfield_constr !! Wheter to apply the fields
+logical, intent(in)                      :: lbfield_trans ! Apply only transveral
+logical, intent(in)                      :: lbfield_mt ! Apply only up do MT radius
+
 
 DOUBLE PRECISION, PARAMETER :: cvlight=274.0720442D0
 DOUBLE COMPLEX, PARAMETER :: czero=(0D0,0D0)
@@ -525,7 +554,12 @@ DO ie=1,ielast
       theta,phi,ipan_intervall,rpan_intervall, npan_tot,ncheb,irmdnew,nrmaxd,  &
       vnspll0,vnspll1,'1',soc)
   
-! Add bfield
+  ! Add bfield
+  if (lbfield) then
+    call add_bfield(bfield, vnspll1, theta, phi, imt, iteration_number, &
+                    itscf0, itscf1, lbfield_constr, lbfield_trans, &
+                    lbfield_mt, .false.)
+  end if
 
 !c extend matrix for the SRA treatment
   vnspll=czero
@@ -586,7 +620,12 @@ DO ie=1,ielast
       ipan_intervall,rpan_intervall,npan_tot,ncheb,  &
       irmdnew,nrmaxd,vnspll0,vnspll1, 'transpose',soc)
   
-! Add bfield
+  ! Add bfield
+  if (lbfield) then
+    call add_bfield(bfield, vnspll1, theta, phi, imt, iteration_number, &
+                    itscf0, itscf1, lbfield_constr, lbfield_trans, &
+                    lbfield_mt, .true.)
+  end if
 
 !c extend matrix for the SRA treatment
   vnspll=czero
@@ -4118,54 +4157,6 @@ n=Ncheb+1
 
 end subroutine
 
-subroutine rotatematrix(mat,theta,phi,lmmax,mode)
-! rotates a matrix in the local frame pointing in
-! the direction of phi and theta to the global frame
-implicit none
-!interface
-double complex,intent(inout)    ::  mat(2*lmmax,2*lmmax)
-double precision,intent(in)     :: phi
-double precision,intent(in)     :: theta
-integer                         :: lmmax
-integer                         :: mode
-!local
-double complex   :: Umat(2*lmmax,2*lmmax)
-double complex   :: Udeggamat(2*lmmax,2*lmmax)
-double complex   :: mattemp(2*lmmax,2*lmmax)
-!double precision :: matmat_zmzm
-
-!***********************************************************************
-! create the rotation matrix:
-!     | cos(theta/2) exp(-i/2 phi)   -sin(theta/2) exp(-i/2 phi) |
-!  U= |                                                          |
-!     | sin(theta/2) exp( i/2 phi)    cos(theta/2) exp( i/2 phi) |
-!
-!  Udegga = transpose(complex conjug ( U ) )
-!***********************************************************************
-
-
-call create_Umatrix(theta,phi,lmmax,Umat,Udeggamat)
-!***********************************************************************
-! calculate matrix in the global frame:
-!
-!  t_glob = U * t_loc * Udegga
-!***********************************************************************
-
-
-if (mode==0) then ! 'loc->glob'
-  call zgemm('N','N',2*lmmax,2*lmmax,2*lmmax,(1d0,0d0),mat,2*lmmax,Udeggamat,2*lmmax,(0d0,0d0),mattemp,2*lmmax)
-  call zgemm('N','N',2*lmmax,2*lmmax,2*lmmax,(1d0,0d0),Umat,2*lmmax,mattemp,2*lmmax,(0d0,0d0),mat,2*lmmax)
-elseif (mode==1) then !'glob->loc'
-  call zgemm('N','N',2*lmmax,2*lmmax,2*lmmax,(1d0,0d0),mat,2*lmmax,Umat,2*lmmax,(0d0,0d0),mattemp,2*lmmax)
-  call zgemm('N','N',2*lmmax,2*lmmax,2*lmmax,(1d0,0d0),Udeggamat,2*lmmax,mattemp,2*lmmax,(0d0,0d0),mat,2*lmmax)
-else
-  stop '[rotatematrix] mode not known'
-end if
-!  writE(324,'(5000F)') tmat
-! stop
-
-end subroutine rotatematrix
-
 
 SUBROUTINE spin_orbit_compl(lmax,lmmaxd,l_s)
 
@@ -4520,80 +4511,6 @@ end do
 end subroutine getCmatrix
 
 
-
-subroutine create_Umatrix(theta,phi,lmmax,Umat,Udeggamat)
-implicit none
-!***********************************************************************
-! create the rotation matrix:
-!     | cos(theta/2) exp(-i/2 phi)   -sin(theta/2) exp(-i/2 phi) |
-!  U= |                                                          |
-!     | sin(theta/2) exp( i/2 phi)    cos(theta/2) exp( i/2 phi) |
-!
-!  Udegga = transpose(complex conjug ( U ) )
-!***********************************************************************double
-!precision :: phi
-!interface
-double precision,intent(in)     :: phi 
-double precision,intent(in)     :: theta
-integer,intent(in)              :: lmmax
-double complex,intent(out)      :: Umat(2*lmmax,2*lmmax)
-double complex,intent(out)      :: Udeggamat(2*lmmax,2*lmmax)
-!local
-double complex                  :: Umat11,Umat12,Umat21,Umat22
-double complex                  :: Udeggamat11,Udeggamat12,Udeggamat21,Udeggamat22
-integer                         :: ival
-double complex,parameter        :: ci=(0.0D0,1.0D0)
-character*25               :: spinmode
-
-spinmode='kkr'
-if (spinmode=='regular') then
-  Umat11      =  cos(theta/2.0D0)*exp(-ci/2.0D0*phi)
-  Umat12      = -sin(theta/2.0D0)*exp(-ci/2.0D0*phi)
-  Umat21      =  sin(theta/2.0D0)*exp( ci/2.0D0*phi)
-  Umat22      =  cos(theta/2.0D0)*exp( ci/2.0D0*phi)
-else if (spinmode=='kkr') then
-  Umat11      =  cos(theta/2.0D0)*exp( ci/2.0D0*phi)
-  Umat12      =  sin(theta/2.0D0)*exp( ci/2.0D0*phi)
-  Umat21      = -sin(theta/2.0D0)*exp(-ci/2.0D0*phi)
-  Umat22      =  cos(theta/2.0D0)*exp(-ci/2.0D0*phi)
-else 
-  stop '[create_Umatrix] mode not known'
-end if
-
-Umat=(0.0D0,0.0D0)
-do ival=1,lmmax
-  Umat(      ival,      ival) = Umat11
-  Umat(      ival,lmmax+ival) = Umat12
-  Umat(lmmax+ival,ival)       = Umat21
-  Umat(lmmax+ival,lmmax+ival) = Umat22
-end do
-
-if (spinmode=='regular') then
-Udeggamat11 =  cos(theta/2.0D0)*exp( ci/2.0D0*phi)
-Udeggamat12 =  sin(theta/2.0D0)*exp(-ci/2.0D0*phi)
-Udeggamat21 = -sin(theta/2.0D0)*exp( ci/2.0D0*phi)
-Udeggamat22 =  cos(theta/2.0D0)*exp(-ci/2.0D0*phi)
-else if (spinmode=='kkr') then
-Udeggamat11 =  cos(theta/2.0D0)*exp(-ci/2.0D0*phi)
-Udeggamat12 = -sin(theta/2.0D0)*exp( ci/2.0D0*phi)
-Udeggamat21 =  sin(theta/2.0D0)*exp(-ci/2.0D0*phi)
-Udeggamat22 =  cos(theta/2.0D0)*exp( ci/2.0D0*phi)
-else 
-  stop '[create_Umatrix] mode not known'
-end if
-
-
-
-Udeggamat=(0.0D0,0.0D0)
-do ival=1,lmmax
-  Udeggamat(      ival,      ival) = Udeggamat11
-  Udeggamat(      ival,lmmax+ival) = Udeggamat12
-  Udeggamat(lmmax+ival,ival)       = Udeggamat21
-  Udeggamat(lmmax+ival,lmmax+ival) = Udeggamat22
-end do
-
-end subroutine create_Umatrix
-
 SUBROUTINE spin_orbit_one_l(lmax,l_s)
 
 IMPLICIT NONE
diff --git a/source/KKRnano/source/ProcessKKRresults_mod.F90 b/source/KKRnano/source/ProcessKKRresults_mod.F90
index 43a67fb69..e54d45425 100644
--- a/source/KKRnano/source/ProcessKKRresults_mod.F90
+++ b/source/KKRnano/source/ProcessKKRresults_mod.F90
@@ -521,7 +521,7 @@ module ProcessKKRresults_mod
                           calc%noco_data%theta_noco_old(atom_id), calc%noco_data%phi_noco_old(atom_id), &
                           calc%noco_data%angle_fixed(atom_id), & 
                           calc%noco_data%moment_x(atom_id),calc%noco_data%moment_y(atom_id), calc%noco_data%moment_z(atom_id), &
-                          densities%muorb, densities%iemxd, params)
+                          densities%muorb, densities%iemxd, params, calc%bfields(atom_id), calc%mesh_a(ila)%imt, iter)
 
       ! LDAU
       if (ldau_data%LDAU .and. ldau_data%NLDAU >= 1) then
diff --git a/source/KKRnano/source/ScatteringCalculation_mod.F90 b/source/KKRnano/source/ScatteringCalculation_mod.F90
index dbdf267b4..101835862 100644
--- a/source/KKRnano/source/ScatteringCalculation_mod.F90
+++ b/source/KKRnano/source/ScatteringCalculation_mod.F90
@@ -75,7 +75,8 @@ implicit none
     use two_sided_commD_mod, only: distribute
     
     use ChebMeshData_mod, only: interpolate_poten  ! NOCO
-    use NonCollinearMagnetism_mod, only: tmat_newsolver, rotatematrix  ! NOCO
+    use NonCollinearMagnetism_mod, only: tmat_newsolver  ! NOCO
+    use NonCollinearMagnetism_Helpers_mod, only: rotatematrix  ! NOCO
 
     integer, intent(in) :: iter
     type(CalculationData), intent(inout) :: calc
@@ -270,7 +271,9 @@ implicit none
                                     noco%theta_noco(i1),noco%phi_noco(i1),1,  & !ipot=1 because potential has only one or two entries (spin polarized case)
                                     !dims%lly,        &    
                                     atomdata%potential%lmpot,atomdata%chebmesh_ptr%irmd_new, &
-                                    kkr(ila)%TmatN(:,:,ispin),params%soc,params%enable_quad_prec)
+                                    kkr(ila)%TmatN(:,:,ispin),params%soc,params%enable_quad_prec, &
+                                    calc%bfields(i1), calc%mesh_a(ila)%imt, iter, params%itbfield0, params%itbfield1, &
+                                    params%noncobfield, params%constr_field, params%trans_bfield, params%mt_bfield)
                
                 call rotatematrix(kkr(ila)%TmatN(:,:,ispin),noco%theta_noco(i1),noco%phi_noco(i1),lmmaxd,0)
               else
diff --git a/source/KKRnano/source/bfield/bfield.f90 b/source/KKRnano/source/bfield/bfield.f90
index 57797d628..7dcb6fecb 100644
--- a/source/KKRnano/source/bfield/bfield.f90
+++ b/source/KKRnano/source/bfield/bfield.f90
@@ -12,7 +12,7 @@
 !------------------------------------------------------------------------------------
 module mod_bfield
   
-  use :: NonCollinearMagnetism_mod, only : rotatematrix
+  use :: NonCollinearMagnetism_Helpers_mod, only : rotatematrix
 
   implicit none
 
diff --git a/source/KKRnano/source/wrappers_mod.F90 b/source/KKRnano/source/wrappers_mod.F90
index befda5ce8..0535615d5 100644
--- a/source/KKRnano/source/wrappers_mod.F90
+++ b/source/KKRnano/source/wrappers_mod.F90
@@ -6,7 +6,10 @@
 !> @todo check intents
 module wrappers_mod
 #include "macros.h"
-use Warnings_mod, only: launch_warning
+  
+  use Warnings_mod, only: launch_warning
+  use mod_bfield, only: bfield_data
+  
   implicit none
   private
   
@@ -27,7 +30,7 @@ use Warnings_mod, only: launch_warning
                             korbit, theta_noco, phi_noco, theta_noco_old, &                          
                             phi_noco_old, angle_fixed, &
                             moment_x, moment_y, moment_z, &
-                            muorb, iemxd, params) ! NOCO/SOC
+                            muorb, iemxd, params, bfield, imt, iteration_number) ! NOCO/SOC
     use BasisAtom_mod, only: BasisAtom
     use GauntCoefficients_mod, only: GauntCoefficients
     use EnergyMesh_mod, only: EnergyMesh
@@ -64,6 +67,9 @@ use Warnings_mod, only: launch_warning
     double precision, intent(out)   :: muorb(0:,:)     ! NOCO
     integer, intent(in)             :: iemxd           ! NOCO
     type(InputParams), intent(in)   :: params          ! NOCO
+    type(bfield_data), intent(in)   :: bfield
+    integer, intent(in)             :: imt
+    integer, intent(in)             :: iteration_number
     
     integer :: ispin, nspind, irmind, irnsd, lmaxd, l
 
@@ -91,7 +97,9 @@ use Warnings_mod, only: launch_warning
                       theta_noco,phi_noco,angle_fixed,moment_x,moment_y,moment_z,&
                       1,  &  ! ipot=1
                       den,espv,rho2ns,r2nef, gmatn(:,:,:,1), muorb,  & ! just one spin component of gmatn needed
-                      atomdata%potential%lpot,lmaxd,mesh%irmd,chebmesh%irmd_new,iemxd, params%soc,params%enable_quad_prec)
+                      atomdata%potential%lpot,lmaxd,mesh%irmd,chebmesh%irmd_new,iemxd, params%soc,params%enable_quad_prec, &
+                      bfield, imt, iteration_number, params%itbfield0, params%itbfield1, &
+                      params%noncobfield, params%constr_field, params%trans_bfield, params%mt_bfield)
  
        ! calculate correct orbital moment
        do ispin=1,nspind
-- 
GitLab