diff --git a/source/KKRnano/source/Makefile b/source/KKRnano/source/Makefile
index 5895e07c0828a3b0a56d1ef40becf3b2d6957713..9bc6e4fe624d983cd41ac190456f5463c914c85b 100644
--- a/source/KKRnano/source/Makefile
+++ b/source/KKRnano/source/Makefile
@@ -561,9 +561,10 @@ 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_Helpers_mod.o bfield.o
+NonCollinearMagnetism_mod.o: RadialMeshData_mod.o ChebMeshData_mod.o read_formatted_shapefun_mod.o NonCollinearMagnetism_Helpers_mod.o bfield.o torque.o
NonCollinearMagnetismData_mod.o: Exceptions_mod.o
vintras_new.o: SingleSiteHelpers_mod.o
bfield.o: NonCollinearMagnetism_Helpers_mod.o
+torque.o: NonCollinearMagnetism_Helpers_mod.o bfield.o
# DO NOT DELETE
diff --git a/source/KKRnano/source/NonCollinearMagnetism_Helpers_mod.F90 b/source/KKRnano/source/NonCollinearMagnetism_Helpers_mod.F90
index 79b183340edd38c01bf81d3b1c71d41514c88c76..446b77292842237bf73d5c89921c5d350bca7e06 100644
--- a/source/KKRnano/source/NonCollinearMagnetism_Helpers_mod.F90
+++ b/source/KKRnano/source/NonCollinearMagnetism_Helpers_mod.F90
@@ -1,7 +1,7 @@
module NonCollinearMagnetism_Helpers_mod
private
- public :: rotatematrix, create_Umatrix
+ public :: rotatematrix, create_Umatrix, intcheb_cell
contains
@@ -126,4 +126,59 @@ module NonCollinearMagnetism_Helpers_mod
end subroutine create_Umatrix
+ subroutine intcheb_cell(cden,den,rpan_intervall,ipan_intervall, &
+ npan_tot,ncheb,irmdnew)
+ !***********************************************************************
+ ! integrate the complex density of states for LM=1
+ ! gives the total complex charge which is then
+ ! transformed to the xyz component of the magnetic
+ ! moment
+ !***********************************************************************
+ implicit none
+
+ integer :: ncheb,npan_tot,irmdnew
+ integer :: ipan_intervall(0:npan_tot)
+ double precision :: rpan_intervall(0:npan_tot)
+ double complex :: cden(irmdnew),den
+ integer :: irstart,irstop,ipan
+ double precision :: widthfac
+ double complex :: int1
+
+ den=(0.0D0,0.0D0)
+
+ do ipan=1,npan_tot
+ irstart=ipan_intervall(ipan-1)+1
+ irstop = ipan_intervall(ipan)
+ widthfac = 0.5D0*(rpan_intervall(ipan)-rpan_intervall(ipan-1))
+ call intcheb_complex(ncheb,cden(irstart:irstop),int1)
+ den=den+int1*widthfac
+ end do
+
+ end subroutine
+
+ subroutine intcheb_complex(ncheb,arr1,result1)
+ implicit none
+ integer, intent(in) :: ncheb
+ double complex, intent(in) :: arr1(0:ncheb)
+ double complex, intent(out) :: result1
+ double precision :: pi
+ double precision :: intweight(0:ncheb)
+ integer :: icheb1,icheb2
+
+ pi=4d0*datan(1d0)
+ intweight=1.0D0
+ do icheb1=0,ncheb
+ do icheb2=2,ncheb,2
+ intweight(icheb1)=intweight(icheb1)+(-2.0D0/(icheb2**2-1.0D0))*dcos(icheb2*pi*(icheb1+0.5D0)/(Ncheb+1))
+ end do
+ intweight(icheb1)=intweight(icheb1)*2.0D0/(Ncheb+1)
+ end do
+
+ result1=(0.0D0,0.0D0)
+ do icheb1=0,ncheb
+ result1=result1+intweight(icheb1)*arr1(icheb1)
+ end do
+
+ end subroutine
+
end module
diff --git a/source/KKRnano/source/NonCollinearMagnetism_mod.F90 b/source/KKRnano/source/NonCollinearMagnetism_mod.F90
index d4065f1414f046a99ab2bfca05310728a938c297..07329ec08fb1fb1c3c1db0e12fedba68034a5630 100644
--- a/source/KKRnano/source/NonCollinearMagnetism_mod.F90
+++ b/source/KKRnano/source/NonCollinearMagnetism_mod.F90
@@ -13,8 +13,9 @@ module NonCollinearMagnetism_mod
!-------------------------------------------------------------------------------
use RadialMeshData_mod!, only:
use ChebMeshData_mod!, only
-use NonCollinearMagnetism_Helpers_mod, only: rotatematrix, create_Umatrix
+use NonCollinearMagnetism_Helpers_mod, only: rotatematrix, create_Umatrix, intcheb_cell
use mod_bfield, only: bfield_data, add_bfield
+use mod_torque, only: calc_torque
implicit none
private
@@ -338,7 +339,7 @@ SUBROUTINE rhovalnew(ldorhoef,ielast,nsra,nspin,lmax,ez,wez,zat, &
den_out,espv,rho2ns,r2nef,gmatn, muorb, &
lpotd,lmaxd,irmd,irmd_new,iemxd,soc,enable_quad_prec, &
bfield, imt, iteration_number, itscf0, itscf1, &
- lbfield, lbfield_constr, lbfield_trans, lbfield_mt)
+ lbfield, lbfield_constr, ltorque, lbfield_trans, lbfield_mt)
! Code converted using TO_F90 by Alan Miller
! Date: 2016-04-21 Time: 11:39:57
@@ -390,11 +391,12 @@ 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
+type(bfield_data), intent(inout) :: 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) :: ltorque !! Whether to calculate torques
logical, intent(in) :: lbfield_trans ! Apply only transveral
logical, intent(in) :: lbfield_mt ! Apply only up do MT radius
@@ -819,6 +821,14 @@ DO ir=1,irmdnew
END DO
r2nefc(:,:,:) = r2nefc(:,:,:) + r2nefc_loop(:,:,:)
+! Calculate torque
+! Compare this position to KKhost: calling after the rho2nsc_loop, before cheb2oldgrid
+if (ltorque) then
+ call calc_torque(bfield, vins, rho2nsc, theta, phi, lmax, rpan_intervall, ipan_intervall, &
+ npan_tot, ncheb, imt1, iend, icleb, cleb, ifunm, thetasnew, &
+ lbfield_constr, lbfield_mt, itscf0, itscf1, iteration_number, angle_fixed)
+end if
+
allocate(rhotemp(irmdnew,lmpotd))
allocate(rhonewtemp(irws,lmpotd))
DO jspin=1,4
@@ -856,7 +866,8 @@ deallocate(rhotemp)
deallocate(rhonewtemp)
! calculate new THETA and PHI for non-colinear
!IF (.NOT.test('FIXMOM ')) THEN
-if (angle_fixed == 0) then ! angle not fixed
+if (angle_fixed == 0 .or. & ! angle not fixed
+ angle_fixed == 4 .or. angle_fixed == 5 ) then ! use constraining fiels without fixing here
rho2ns_temp(1,1)=rho2int(1)
rho2ns_temp(2,2)=rho2int(2)
rho2ns_temp(1,2)=rho2int(3)
@@ -963,7 +974,8 @@ DO lm1=0,lmaxd1
END DO
! UPDATE ANGLES
-if (angle_fixed == 0) then
+if (angle_fixed == 0 .or. & ! angle not fixed
+ angle_fixed == 4 .or. angle_fixed == 5 ) then ! use constraining fiels without fixing here
phi = phinew
theta = thetanew
endif
@@ -4668,61 +4680,6 @@ deallocate(l_up)
END SUBROUTINE spin_orbit_one_l
-subroutine intcheb_cell(cden,den,rpan_intervall,ipan_intervall, &
- npan_tot,ncheb,irmdnew)
-!***********************************************************************
-! integrate the complex density of states for LM=1
-! gives the total complex charge which is then
-! transformed to the xyz component of the magnetic
-! moment
-!***********************************************************************
-implicit none
-
-integer :: ncheb,npan_tot,irmdnew
-integer :: ipan_intervall(0:npan_tot)
-double precision :: rpan_intervall(0:npan_tot)
-double complex :: cden(irmdnew),den
-integer :: irstart,irstop,ipan
-double precision :: widthfac
-double complex :: int1
-
-den=(0.0D0,0.0D0)
-
- do ipan=1,npan_tot
- irstart=ipan_intervall(ipan-1)+1
- irstop = ipan_intervall(ipan)
- widthfac = 0.5D0*(rpan_intervall(ipan)-rpan_intervall(ipan-1))
- call intcheb_complex(ncheb,cden(irstart:irstop),int1)
- den=den+int1*widthfac
- end do
-
-end subroutine intcheb_cell
-
-subroutine intcheb_complex(ncheb,arr1,result1)
-implicit none
-integer, intent(in) :: ncheb
-double complex, intent(in) :: arr1(0:ncheb)
-double complex, intent(out) :: result1
-double precision :: pi
-double precision :: intweight(0:ncheb)
-integer :: icheb1,icheb2
-
-pi=4d0*datan(1d0)
-intweight=1.0D0
- do icheb1=0,ncheb
- do icheb2=2,ncheb,2
- intweight(icheb1)=intweight(icheb1)+(-2.0D0/(icheb2**2-1.0D0))*dcos(icheb2*pi*(icheb1+0.5D0)/(Ncheb+1))
- end do
- intweight(icheb1)=intweight(icheb1)*2.0D0/(Ncheb+1)
- end do
-
-result1=(0.0D0,0.0D0)
-do icheb1=0,ncheb
- result1=result1+intweight(icheb1)*arr1(icheb1)
-end do
-
-end subroutine
-
subroutine cheb2oldgrid(nrmax,nrmaxnew,lmmaxpot,rmesh,ncheb, &
npan_tot,rpan_intervall,ipan_intervall, &
arrayin,arrayout,irmd)
diff --git a/source/KKRnano/source/bfield/torque.f90 b/source/KKRnano/source/bfield/torque.f90
new file mode 100644
index 0000000000000000000000000000000000000000..5cb0cd20b11a9e8c0fcf7a6efe553166ea5ad18b
--- /dev/null
+++ b/source/KKRnano/source/bfield/torque.f90
@@ -0,0 +1,263 @@
+!-----------------------------------------------------------------------------------------!
+! Copyright (c) 2018 Peter Grünberg Institut, Forschungszentrum Jülich, Germany !
+! This file is part of Jülich KKR code and available as free software under the conditions!
+! of the MIT license as expressed in the LICENSE.md file in more detail. !
+!-----------------------------------------------------------------------------------------!
+
+!------------------------------------------------------------------------------------
+!> Summary: Calculation of the magnetic torque methods
+!>
+!> Author: Sascha Brinker, Nicolas Essing
+!------------------------------------------------------------------------------------
+module mod_torque
+ use mod_bfield, only: bfield_data
+ use NonCollinearMagnetism_Helpers_mod, only: rotatematrix, intcheb_cell
+
+ implicit none
+
+contains
+
+ !-------------------------------------------------------------------------------
+ !> Summary: Calculation of the density for the new solver
+ !> Author: Sascha Brinker, Nicolas Essing
+ !-------------------------------------------------------------------------------
+ subroutine calc_torque(bfield, vpot, rho2nsc, theta, phi, lmax, rpan_intervall, ipan_intervall, &
+ npan_tot, ncheb, imt, iend, icleb, cleb, ifunm, thetasnew, &
+ lbfield_constr, lbfield_mt, itscf0, itscf1, iteration, constr_mode)
+ type(bfield_data), intent(inout) :: bfield !! Information on the magnetic field
+ double precision, dimension(:,:,:), intent(in) :: vpot !! The potential
+ double complex, dimension(:,:,:), intent(in) :: rho2nsc !! complex density matrix
+ double precision, intent(in) :: theta, phi !! Angles of the (old) magnetic moment
+ integer, intent(in) :: lmax !! Angular momentum cutoff
+ double precision, dimension(:), intent(in) :: rpan_intervall
+ integer, dimension(:), intent(in) :: ipan_intervall
+ integer, intent(in) :: npan_tot
+ integer, intent(in) :: ncheb
+ integer, intent(in) :: imt !! Index of the muffin tin radius
+ integer, intent(in) :: iend !! Number of nonzero gaunt coefficients
+ integer, dimension(:,:), intent(in) :: icleb !! l and m values for the Gaunt coefficients
+ double precision, dimension(:), intent(in) :: cleb !! Gaunt coefficients
+ integer , dimension(:), intent(in) :: ifunm !! pointer array for shapefun
+ double precision, dimension(:, :), intent(in) :: thetasnew !! shapefun on the Cheby mesh
+ logical, intent(in) :: lbfield_constr !! Use constraint fields
+ logical, intent(in) :: lbfield_mt !! Use magnetic fields only inside the muffin tin
+ integer, intent(in) :: itscf0, itscf1 !! Apply magnetic fields between these iterations
+ integer, intent(in) :: iteration !! Current iteration
+ integer, intent(in) :: constr_mode !! Mode of the constraining field self-consistency
+
+ double complex, parameter :: cone = (1., 0.) ! complex one
+ double precision, parameter :: rfpi = 3.5449077018110318 ! sqrt(4*pi)
+
+ double precision, parameter :: constraint_bfields_mixing_parameter = 0.03
+
+ integer :: irmd, lmmax, lmpotd ! radial and angular momentum sizes
+ integer :: i, ir, ilm ! generic, radial and angular momentum loop indices
+ double complex :: integrate_result ! to output the result of some integrations
+ double precision :: mag_mom_len ! absolute value of magnetic moment
+ double precision, dimension(3) :: dir, mag_mom_dir ! unit vectors of old and new magnetic moment
+ double precision, dimension(3) :: torque, torque_mt ! torque calculated over whole cell or mt only
+ double precision, dimension(3) :: old_b_constr ! temporary storage for the old constraint field
+ double complex, dimension(2,2) :: rho2ns_temp ! temporary matrix to rotate rho2nsc
+ double complex, dimension(:), allocatable :: integrand ! temporary array to integrate stuff
+ double precision, dimension(:,:), allocatable :: bxc ! xc magnetic field in collinear form
+ double precision, dimension(:,:,:), allocatable :: mag_den ! magnetization density
+ double precision, dimension(:,:,:), allocatable :: mag_den_conv ! mag. den. convoluted with shapefun
+ double precision, dimension(:,:), allocatable :: mag ! integrated lm-resolved magnetic moment
+ double precision, dimension(:,:), allocatable :: mag_mt ! integrated lm-resolved magnetic moment up to mt
+
+ ! Get some dimensions
+ lmmax = (lmax+1)**2
+ irmd = size(vpot,1)
+ lmpotd = size(vpot,2)
+
+ ! Allocate temporary arrays
+ allocate(integrand(irmd))
+ allocate(bxc(irmd, lmpotd))
+ allocate(mag_den(irmd, lmpotd, 3))
+ allocate(mag(lmmax, 3))
+ allocate(mag_mt(lmmax, 3))
+
+ ! Get the xc magnetic field. In the rigid spin approximation, the direction
+ ! of the field is fixed by the local frame, but it is radial dependent.
+ ! Can be calculated as the difference of spin up and spin down part of the
+ ! potential.
+ bxc(:,:) = (vpot(:,:,1) - vpot(:,:,2)) / 2.
+
+ ! Get magnetization density (in the global frame). For each angular momentum
+ ! and radial index, rotate to the global frame.
+ do ilm = 1, lmpotd
+ do ir = 1, irmd
+ rho2ns_temp(1,1) = rho2nsc(ir,ilm,1)
+ rho2ns_temp(2,2) = rho2nsc(ir,ilm,2)
+ rho2ns_temp(1,2) = rho2nsc(ir,ilm,3)
+ rho2ns_temp(2,1) = rho2nsc(ir,ilm,4)
+
+ call rotatematrix(rho2ns_temp, theta, phi, 1, 0)
+
+ mag_den(ir,ilm,1) = aimag( rho2ns_temp(1,2) + rho2ns_temp(2,1) )
+ mag_den(ir,ilm,2) = real( rho2ns_temp(2,1) - rho2ns_temp(1,2) )
+ mag_den(ir,ilm,3) = aimag( rho2ns_temp(2,2) - rho2ns_temp(1,1) )
+ end do
+ end do
+
+ ! Integrate magnetization density (to mt or end depending on parameter)
+ ! to get magnetic moment. First convolute with the shapefun.
+ call calc_mag_mom(lmax, imt, iend, icleb, cleb, ifunm, thetasnew, mag_den, mag_den_conv)
+ do i = 1, 3
+ do ilm = 1, lmmax
+ integrand(:) = mag_den_conv(:, ilm, i) * cone
+ call intcheb_cell(integrand,integrate_result,rpan_intervall,ipan_intervall,npan_tot,ncheb,irmd)
+ mag(ilm,i) = real(integrate_result * rfpi)
+ ! Same for only the muffin tin: Set to zero outside, integrate again
+ integrand(imt:) = (0., 0.)
+ call intcheb_cell(integrand(:),integrate_result,rpan_intervall,ipan_intervall,npan_tot,ncheb,irmd)
+ mag_mt(ilm,i) = real(integrate_result * rfpi)
+ end do
+ end do
+
+ !TODO Output magnetic moments
+
+ ! Calculate direction and absolute value of the magnetic moment
+ if (lbfield_mt) then
+ mag_mom_len = sqrt(dot_product(mag_mt(1,:), mag_mt(1,:)))
+ mag_mom_dir = mag_mt(1,:) / mag_mom_len
+ else
+ mag_mom_len = sqrt(dot_product(mag(1,:), mag(1,:)))
+ mag_mom_dir = mag(1,:) / mag_mom_len
+ end if
+
+ ! Calculate the torque (also mt or end)
+ ! Integrate the xc bfield times the magnetization density.
+ torque(:) = 0.
+ torque_mt(:) = 0.
+ do i = 1, 3
+ do ilm = 1, lmmax
+ integrand(:) = bxc(:, ilm) * mag_den(:, ilm, i) * cone
+ call intcheb_cell(integrand(:),integrate_result,rpan_intervall,ipan_intervall,npan_tot,ncheb,irmd)
+ torque(i) = torque(i) + real(integrate_result * rfpi) ! rfpi only contained in convoluted quantities
+ ! Same for mt
+ integrand(imt:) = 0.
+ call intcheb_cell(integrand(:),integrate_result,rpan_intervall,ipan_intervall,npan_tot,ncheb,irmd)
+ torque_mt(i) = torque_mt(i) + real(integrate_result * rfpi)
+ end do
+ end do
+
+ ! Project torque to perpendicular of magnetic moment
+ dir(1) = sin(theta) * cos(phi)
+ dir(2) = sin(theta) * sin(phi)
+ dir(3) = cos(theta)
+ torque(:) = torque(:) - dir(:) * dot_product(dir(:), torque(:))
+ torque_mt(:) = torque_mt(:) - dir(:) * dot_product(dir(:), torque_mt(:))
+
+ !TODO Output torque
+
+ ! Save torque in bfield_data
+ if (lbfield_mt) then
+ bfield%mag_torque(:) = torque_mt
+ else
+ bfield%mag_torque(:) = torque
+ end if
+
+ ! Scf-cycle for constraint fields, based either on torque or on fields alone
+ if (lbfield_constr .and. itscf0 <= iteration .and. iteration <= itscf1) then
+ if (constr_mode == 3 .or. constr_mode == 5) then
+ bfield%bfield_constr(:) = bfield%bfield_constr(:) - torque(:) / mag_mom_len
+ else if (constr_mode == 2 .or. constr_mode == 4) 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)*constraint_bfields_mixing_parameter
+ else
+ ! There might be other modes that are calculated somewhere else. Do nothing.
+ end if
+ end if
+
+ ! Deallocate
+ deallocate(bxc, integrand, mag_den, mag_den_conv, mag, mag_mt)
+
+ end subroutine calc_torque
+
+ !------------------------------------------------------------------------------------
+ !> Summary: Convolute the magnetic moments with the shapefunction.
+ !>
+ !> Author: Sascha Brinker, Nicolas Essing
+ !------------------------------------------------------------------------------------
+ subroutine calc_mag_mom(lmax, imt, iend, icleb, cleb, ifunm, thetasnew, mag_den, mag_den_conv)
+ integer, intent(in) :: lmax !! Angular momentum cut-off
+ integer, intent(in) :: imt !! index muffin-tin radius
+ integer, intent(in) :: iend !! Number of nonzero gaunt coefficients
+ integer, dimension(:,:), intent(in) :: icleb !! l and m values for the Gaunt coefficients
+ double precision, dimension(:), intent(in) :: cleb !! Gaunt coefficients
+ integer , dimension(:), intent(in) :: ifunm !! pointer array for shapefun
+ double precision, dimension(:, :), intent(in) :: thetasnew !! shapefun on the Cheby mesh
+ double precision, dimension(:,:,:), intent(in) :: mag_den !! uncovoluted magnetization density
+ double precision, dimension(:,:,:), allocatable, intent(out) :: mag_den_conv !! result
+
+ double precision, parameter :: rfpi = 3.5449077018110318 ! sqrt(4*pi)
+ double precision, parameter :: c0ll = 1.d0 / rfpi
+ integer :: lmmax, lmpotd ! Angular momentum sizes
+ integer :: irmdnew ! Radial size
+ integer :: i, j ! Generic loop indices
+ integer :: lm1, lm2, lm3 ! Angular momentum indices
+ integer :: ifun ! Pointer index for thetasnew
+ double precision, dimension(:,:), allocatable :: shapefun_mod ! Complete shapefun
+
+ ! Get some dimensions
+ lmmax = (lmax+1)**2
+ lmpotd = size(mag_den, 2)
+ irmdnew = size(thetasnew, 1) ! npan_tot*(ncheb+1)
+
+ ! Allocate temporary shapefun and the result
+ allocate(shapefun_mod(irmdnew, lmpotd))
+ allocate(mag_den_conv(irmdnew, lmmax, 3))
+
+ ! Build the shapefun array. Start with all zero. Inside muffin tin only l=0
+ ! component is /= 0 (and constant), copy l=0 component for r outside of mt
+ ! from thetasnew.
+ shapefun_mod(:,:) = 0.d0
+ shapefun_mod(1:imt, 1) = rfpi ! is multipled by C_LL^0
+ shapefun_mod(imt+1:, 1) = thetasnew(imt+1:, 1)
+
+ ! Copy other components from thetasnew. Convert from pointer indices to
+ ! normal (l,m)-index
+ do lm1 = 2, lmpotd
+ ifun = ifunm(lm1)
+ if(ifun /= 0) then !shapefun%lmused(lm1)==1) then
+ shapefun_mod(imt+1:, lm1) = thetasnew(imt+1:, ifun)
+ end if
+ end do
+
+ ! Initialize convoluted magnetization density.
+ mag_den_conv(:,:,:) = 0.
+
+ ! Diagonal part (not contained in gaunt-coeff)
+ do i = 1, 3
+ do lm1 = 1, lmpotd
+ mag_den_conv(:,1,i) = mag_den_conv(:,1,i) + mag_den(:,lm1,i) * shapefun_mod(:,lm1) * c0ll
+ end do
+ end do
+
+ ! Offdiagonal part. This is effectively a loop over angular momentum indices
+ ! lm1,lm2,lm3. Iterate instead over the flattened array cleb
+ ! containing the Gaunt coefficients and extract the angular momentum
+ ! indices for each j.
+ do i = 1, 3
+ do j = 1, iend
+ lm1 = icleb(j, 1)
+ lm2 = icleb(j, 2)
+ lm3 = icleb(j, 3)
+ if (lm1 > lmmax .or. lm2 > lmmax .or. lm3 > lmmax) then
+ !TODO I think this should not happen, as icleb should only contain
+ ! valid values. It might be defined for higher lm3. Check that.
+ cycle
+ end if
+ mag_den_conv(:,lm3,i) = mag_den_conv(:,lm3,i) + cleb(j) * mag_den(:,lm1,i) * shapefun_mod(:,lm2)
+ mag_den_conv(:,lm3,i) = mag_den_conv(:,lm3,i) + cleb(j) * mag_den(:,lm2,i) * shapefun_mod(:,lm1)
+ end do
+ end do
+
+ ! Deallocate temporary shapefun
+ deallocate(shapefun_mod)
+
+ end subroutine
+
+end module mod_torque
diff --git a/source/KKRnano/source/wrappers_mod.F90 b/source/KKRnano/source/wrappers_mod.F90
index 0535615d5390115d6309ca9cd9889abdada6da06..cd698824c58abbf831ec7c7fa951205055815b89 100644
--- a/source/KKRnano/source/wrappers_mod.F90
+++ b/source/KKRnano/source/wrappers_mod.F90
@@ -67,7 +67,7 @@ module wrappers_mod
double precision, intent(out) :: muorb(0:,:) ! NOCO
integer, intent(in) :: iemxd ! NOCO
type(InputParams), intent(in) :: params ! NOCO
- type(bfield_data), intent(in) :: bfield
+ type(bfield_data), intent(inout):: bfield
integer, intent(in) :: imt
integer, intent(in) :: iteration_number
@@ -99,7 +99,7 @@ module wrappers_mod
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, &
bfield, imt, iteration_number, params%itbfield0, params%itbfield1, &
- params%noncobfield, params%constr_field, params%trans_bfield, params%mt_bfield)
+ params%noncobfield, params%constr_field, params%torque, params%trans_bfield, params%mt_bfield)
! calculate correct orbital moment
do ispin=1,nspind