From 3a37ec68fb4f243fcbc714a7c4f0bacb3920e1a7 Mon Sep 17 00:00:00 2001
From: Nicolas Essing <n.essing@fz-juelich.de>
Date: Mon, 18 Oct 2021 16:07:56 +0200
Subject: [PATCH] Implemented torque.
The subroutine calc_torque() in bfield/torque.f90 calculates the
magnetic torque and iterates the selfconsistency loop for the
constraint magnetic fields, based on torque or on the fields alone.
The routine is called from rhovalnew() as in KKRhost. Adapted the
inputs and included the call. Adapted wrappers_mod, where this
routine is called.
I had to move another subroutine from NonCollinearMagnetism_mod.F90
to the helpers module to avoid circular dependencies. Did not change
the routine otherwise.
Included the new dependencies in the Makefile.
Compiles and did not crash in a short test run,
but not completely tested.
---
source/KKRnano/source/Makefile | 3 +-
.../NonCollinearMagnetism_Helpers_mod.F90 | 57 +++-
.../source/NonCollinearMagnetism_mod.F90 | 77 ++---
source/KKRnano/source/bfield/torque.f90 | 263 ++++++++++++++++++
source/KKRnano/source/wrappers_mod.F90 | 4 +-
5 files changed, 340 insertions(+), 64 deletions(-)
create mode 100644 source/KKRnano/source/bfield/torque.f90
diff --git a/source/KKRnano/source/Makefile b/source/KKRnano/source/Makefile
index 5895e07c0..9bc6e4fe6 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 79b183340..446b77292 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 d4065f141..07329ec08 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 000000000..5cb0cd20b
--- /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 0535615d5..cd698824c 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
--
GitLab