From 86978a320e326df3bfb28314d69d50b70fe3eaa9 Mon Sep 17 00:00:00 2001
From: Nicolas Essing <n.essing@fz-juelich.de>
Date: Fri, 8 Oct 2021 13:05:51 +0200
Subject: [PATCH] Fixed a bug concerning the initialization of the magnetic
fields.
In CalculationData, the bfield is read for every atom and then initialized
only for the locally treated atoms. Thus, the bfield has to be accessed
differently from the other quantities that are only allocated for the
local atoms.
Also, I restructured and commented lots of code while looking for the bug.
---
source/KKRnano/source/CalculationData_mod.F90 | 8 +-
source/KKRnano/source/bfield/bfield.f90 | 223 ++++++++----------
2 files changed, 102 insertions(+), 129 deletions(-)
diff --git a/source/KKRnano/source/CalculationData_mod.F90 b/source/KKRnano/source/CalculationData_mod.F90
index 9b7d382f8..325be4c22 100644
--- a/source/KKRnano/source/CalculationData_mod.F90
+++ b/source/KKRnano/source/CalculationData_mod.F90
@@ -387,12 +387,12 @@ module CalculationData_mod
! For the local atoms, initialize some fields
do ila = 1, self%num_local_atoms
atom_id = self%atom_ids(ila)
+ ! Beware: self%bfields is allocated and saved for all atoms
call init_bfield(params%noncobfield, params%constr_field, self%bfields(atom_id), dims%lmaxd, &
- self%cheb_mesh_a(atom_id)%npan_tot, self%cheb_mesh_a(atom_id)%npan_lognew, &
- self%cheb_mesh_a(atom_id)%npan_eqnew, self%cheb_mesh_a(atom_id)%ncheb, &
- self%cheb_mesh_a(atom_id)%ipan_intervall, self%cheb_mesh_a(atom_id)%thetasnew, &
+ self%cheb_mesh_a(ila)%npan_lognew, self%cheb_mesh_a(ila)%npan_eqnew, &
+ self%cheb_mesh_a(ila)%ipan_intervall, self%cheb_mesh_a(ila)%thetasnew, &
self%gaunts%iend, self%gaunts%icleb, self%gaunts%cleb(:,1), &
- self%cell_a(atom_id)%ifunm)
+ self%cell_a(ila)%ifunm)
end do
endsubroutine ! constructEverything
diff --git a/source/KKRnano/source/bfield/bfield.f90 b/source/KKRnano/source/bfield/bfield.f90
index 48382a9a7..5c62f2273 100644
--- a/source/KKRnano/source/bfield/bfield.f90
+++ b/source/KKRnano/source/bfield/bfield.f90
@@ -7,9 +7,8 @@
!------------------------------------------------------------------------------------
!> Summary: Module storing the run options and the paramters for bfields and constraining fields
!>
-!> Written by Sascha Brinker, ported from KKRhost to KKRnano by the other two.
-!>
!> Author: Sascha Brinker, Eduardo Mendive, Nicolas Essing
+!> Written by Sascha Brinker, ported from KKRhost to KKRnano by the other two.
!------------------------------------------------------------------------------------
module mod_bfield
@@ -58,16 +57,14 @@ contains
end subroutine
subroutine init_bfield(lbfield, lbfield_constr, bfield, lmax, &
- npan_tot, npan_log, npan_eq, ncheb, ipan_intervall, thetasnew, &
+ npan_log, npan_eq, ipan_intervall, thetasnew, &
iend, icleb, cleb, ifunm)
logical, intent(in) :: lbfield !! Whether to use external nonco bfields
logical, intent(in) :: lbfield_constr !! Whether to use constraint bfields
type(bfield_data), intent(inout) :: bfield !! The bfield data
integer, intent(in) :: lmax !! Angular momentum cutoff
- integer, intent(in) :: npan_tot !! Chebyshev mesh resolution
integer, intent(in) :: npan_log !! Chebyshev mesh resolution
integer, intent(in) :: npan_eq !! Chebyshev mesh resolution
- integer, intent(in) :: ncheb !! Chebyshev mesh resolution
integer, dimension(:), intent(in) :: ipan_intervall !! Indices for important radial points in the mesh
double precision, dimension(:,:), intent(in) :: thetasnew !! Interpolated shape function in Chebychev mesh
integer, intent (in) :: iend !! Number of nonzero gaunt coefficients
@@ -79,9 +76,7 @@ contains
integer :: i_stat
ncleb = size(cleb)
- irmdnew = npan_tot*(ncheb+1)
-
- !TODO? assert size(theasnew) == irmdnew
+ irmdnew = size(thetasnew, 1) ! npan_tot*(ncheb+1)
! Calculate the LL' expansion of the shape function in the new mesh which
! is needed to convolute the magnetic field (only done once and stored to
@@ -90,7 +85,7 @@ contains
if(lbfield .or. lbfield_constr) then
allocate (bfield%thetallmat((lmax+1)**2,(lmax+1)**2,irmdnew), stat=i_stat)
call calc_thetallmat(bfield%thetallmat, lmax, ipan_intervall(npan_log+npan_eq) + 1, &
- iend, irmdnew, thetasnew, ifunm, icleb, cleb)
+ iend, irmdnew, thetasnew, ifunm, icleb, cleb)
end if
end subroutine init_bfield
@@ -99,11 +94,11 @@ contains
!-------------------------------------------------------------------------------
!> Summary: Reads the atom-wise constraining field from bconstr.dat
- !> Author: MdSD
- !> Category: KKRhost, bfield
- !> Deprecated: False
- !> the file has the format:
- !> bz by bz (in Ry), mspin (in muB)
+ !> Author: MdSD, Nicolas Essing
+ !>
+ !> The file should contain three floating point values in each line and a line
+ !> for each atom. The first line is treated as a header and skipped.
+ !> Intepreted as initial constraining bfield in cartesian coordinates, in Ry.
!-------------------------------------------------------------------------------
subroutine read_bconstr(bfields, number_of_atoms)
type(bfield_data), dimension(:), intent(inout) :: bfields
@@ -111,24 +106,23 @@ contains
integer :: iatom
integer :: iostat
-
- open(unit=57493215, file='bconstr_in.dat', iostat=iostat)
- if (iostat == 0) then
- read(57493215,*) ! skip header
+ logical :: file_exists
+
+ inquire(file='bconstr_in.dat', exist=file_exists)
+
+ if (file_exists) then
+ open(unit=57493215, file='bconstr_in.dat', iostat=iostat)
+ read(57493215, *) ! skip header
do iatom = 1, number_of_atoms
- read(57493215,*,iostat=iostat) bfields(iatom)%bfield_constr(:)
- if (iostat > 0) then
- ! iostat should be zero normally and negative in last line, positive is error
+ read(57493215, *, iostat=iostat) bfields(iatom)%bfield_constr(:)
+ if (iostat /= 0) then
write(*,*) "Error reading bconstr_in.dat"
stop
end if
end do
- if (iostat >= 0) then
- ! still zero means there are more lines
- write(*,*) "Error reading bconstr_in.dat or file too long"
- end if
close(57493215)
else
+ ! No 'bconstr_in.dat' given, use default 0
do iatom = 1, number_of_atoms
bfields(iatom)%bfield_constr(:) = 0.
end do
@@ -146,11 +140,12 @@ contains
!-------------------------------------------------------------------------------
!> Summary: Reads the atom-wise magnetic field from bfield.dat
- !> Author: Sascha Brinker
- !> Category: KKRhost, bfield
- !> Deprecated: False
- !> the file has the format:
- !> theta phi bfield_strength (in Tesla)
+ !> Author: Nicolas Essing
+ !>
+ !> The file should contain three numbers per line and one line per atom.
+ !> Read as 'theta phi bfield_strength' with the angles in degrees and the
+ !> strength in Ry ( ! might get changed to Tesla ! )
+ !> Lines can be commented out with a # as first character.
!-------------------------------------------------------------------------------
subroutine read_bfield(bfields,number_of_atoms)
type(bfield_data), dimension(:), intent(inout) :: bfields
@@ -166,22 +161,23 @@ contains
stop
end if
+ ! manual loop over iatom because comments are allowed
iatom = 1
- do while ((iostat == 0) .and. (iatom <= number_of_atoms))
- read(57493215, *, iostat=iostat) linebuffer
- if (linebuffer(1:1) == '#') continue ! input line commented out, get next line
+ do while (iatom <= number_of_atoms)
+ read(57493215, '(A)', iostat=iostat) linebuffer
+ if (iostat /= 0) then
+ write(*,*) "I/O Error while reading bfield.dat"
+ stop
+ end if
+ if (linebuffer(1:1) == '#') cycle ! input line commented out
read(linebuffer, *, iostat=iostat) bfields(iatom)%theta, bfields(iatom)%phi, bfields(iatom)%bfield_strength
if (iostat /= 0) then
- write(*,*) "Malformatted line in bfield.dat."
+ write(*,*) "Error parsing a line in bfield.dat"
stop
end if
iatom = iatom + 1
end do
close(57493215)
- if (iatom <= number_of_atoms + 1) then
- write(*,*) "Error reading bfield.dat or file too short"
- stop
- end if
write(*,*) ' ###############################################'
write(*,*) ' external non-collinear magnetic fields'
@@ -194,7 +190,7 @@ contains
bfields(iatom)%bfield(1) = bfields(iatom)%bfield_strength*cos(bfields(iatom)%phi)*sin(bfields(iatom)%theta)
bfields(iatom)%bfield(2) = bfields(iatom)%bfield_strength*sin(bfields(iatom)%phi)*sin(bfields(iatom)%theta)
bfields(iatom)%bfield(3) = bfields(iatom)%bfield_strength*cos(bfields(iatom)%theta)
- write(*,'(" ",i4,3es16.8)') iatom, bfields(iatom)%theta, bfields(iatom)%phi, bfields(iatom)%bfield_strength
+ write(*,'(2X,I4,3(E16.8))') iatom, bfields(iatom)%theta, bfields(iatom)%phi, bfields(iatom)%bfield_strength
end do
end subroutine read_bfield
@@ -211,99 +207,76 @@ contains
!> @endnote
!------------------------------------------------------------------------------------
subroutine calc_thetallmat(thetansll, lmax, imt1, iend, irmdnew, thetasnew, ifunm, icleb, cleb)
- integer , intent (in) :: lmax !! Angular momentum cut-off
- integer , intent (in) :: imt1 !! index muffin-tin radius
- integer , intent (in) :: iend !! Number of nonzero gaunt coefficients
- integer , intent (in) :: irmdnew !! number of radials point on the Cheby mesh
- integer, dimension(:,:), intent (in) :: icleb !! Pointer array
- double precision, dimension(:), intent (in) :: cleb !! GAUNT coefficients (GAUNT)
- double precision, dimension(:, :), intent (in) :: thetasnew !! shapefun on the Cheby mesh
- integer, dimension(1:(2*lmax+1)**2) , intent (in) :: ifunm !! pointer array for shapefun ! Check index and dimensions!!!!!
- double precision, dimension((lmax+1)**2,(lmax+1)**2,irmdnew), intent (out) :: thetansll !! LL' expansion of the shapefunction
- !------------------------------------------------------------------------------------
- double precision,parameter :: rfpi=3.5449077018110318
- double precision,dimension(1:irmdnew,1:(2*lmax+1)**2) :: shapefun_mod
- double precision :: c0ll
- integer :: lmmax
- integer :: lmmax2
- integer :: ilm
- integer :: ifun
- integer :: ir
- integer :: lm1
- integer :: lm2
- integer :: lm3
- integer :: j
-
- c0ll = 1d0/rfpi
- lmmax = (lmax+1)**2
- lmmax2 = (2*lmax+1)**2
- shapefun_mod(:,:) = 0.d0
- shapefun_mod(1:imt1,1) = rfpi ! is multipled by C_LL^0
- shapefun_mod(imt1+1:irmdnew,1) = thetasnew(imt1+1:irmdnew,1)
- ! convert from pointer to real ilmdo ilm=2,lmmax2
- do ilm=2,lmmax2
- ifun = ifunm(ilm)
- if(.not. ifun == 0) then !shapefun%lmused(ilm)==1) then
- !if(myrank==master .and. debug) write(17*iatom**2,'(" converted ifun= ",i4," to ilm= ",i4)') ifun, ilm
- shapefun_mod(imt1+1:irmdnew,ilm) = thetasnew(imt1+1:irmdnew,ifun)
+ double precision, dimension(:,:,:), intent (out) :: thetansll !! LL' expansion of the shapefunction
+ integer, intent (in) :: lmax !! Angular momentum cut-off
+ integer, intent (in) :: imt1 !! 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
+ double precision, dimension(:, :), intent (in) :: thetasnew !! shapefun on the Cheby mesh
+ integer , dimension(:), intent (in) :: ifunm !! pointer array for shapefun
+
+ double precision, parameter :: rfpi = 3.5449077018110318
+ double precision, parameter :: c0ll = 1.d0 / rfpi
+ integer :: irmdnew ! number of radials point on the Chebyshev mesh
+ integer :: lmmax ! number of angular momentum entries
+ integer :: lmmax2
+ integer :: ifun
+ integer :: lm1, lm2, lm3 ! angular momentum indices
+ integer :: ir ! radial index
+ integer :: j ! index for loop over Gaunt coeffs
+ double precision, dimension(:,:), allocatable :: shapefun_mod
+
+ lmmax = (lmax+1)**2
+ lmmax2 = (2*lmax+1)**2
+ irmdnew = size(thetasnew, 1) ! npan_tot*(ncheb+1)
+
+ allocate(shapefun_mod(irmdnew, lmmax2))
+
+ shapefun_mod(:,:) = 0.d0
+ shapefun_mod(1:imt1,1) = rfpi ! is multipled by C_LL^0
+ shapefun_mod(imt1+1:, 1) = thetasnew(imt1+1:irmdnew,1)
+
+ ! convert from pointer to real l,m
+ do lm1 = 2, lmmax2
+ ifun = ifunm(lm1)
+ if(.not. ifun == 0) then !shapefun%lmused(lm1)==1) then
+ shapefun_mod(imt1+1:, lm1) = thetasnew(imt1+1:, ifun)
end if
end do
- !if(myrank==master .and. debug) write(17*iatom**2,'(" cellnew%shapefun = ")')
- !do ilm= 1,lmmax2
- ! if(sum(abs(thetasnew(:,ilm)))>1e-8) then
- ! if(myrank==master .and. debug) write(17*iatom**2,'(i4,1000es16.8)') ilm, thetasnew(:,ilm)
- ! end if
- !end do
- !if(myrank==master .and. debug) write(17*iatom**2,'(" shapefun_mod = ")')
- !do ilm= 1,lmmax2
- ! if(sum(abs(shapefun_mod(:,ilm)))>1e-8) then
- ! if(myrank==master .and. debug) write(17*iatom**2,'(i4,1000es16.8)') ilm, shapefun_mod(:,ilm)
- ! end if
- !end do
- thetansll(:,:,:)=0.d0
+ thetansll(:,:,:) = 0.d0
+
! diagonal part (not contained in gaunt-coeff)
- do ilm = 1,lmmax
- thetansll(ilm,ilm,:) = shapefun_mod(:,1)*c0ll
+ do lm1 = 1, lmmax
+ thetansll(lm1,lm1,:) = shapefun_mod(:,1) * c0ll
end do
- !write(*,'("ifunm=",1000i4)') ifunm(:)
- do j = 1,iend !gauntcoeff%iend
- lm1 = icleb(j, 1)! gauntcoeff%icleb(j,1) ! lmax
- lm2 = icleb(j, 2)! gauntcoeff%icleb(j,2) ! lmax
- lm3 = icleb(j, 3)! gauntcoeff%icleb(j,3) ! 2*lmax
- !write(17*iatom**2,'("lm1,lm2,lm3 =",3i4,2es16.8)') lm1,lm2,lm3,cleb(j)
- if(lm1<= lmmax .and. lm2 <= lmmax .and. lm3<= lmmax2) then
- ifun = ifunm(lm3)
- if(.not. ifun == 0) then !shapefun%lmused(ilm)==1) then
- !write(*,'("lm1,lm2,lm3,cleb",3i4,10e16.8)') lm1,lm2,lm3,cleb(j)
- do ir = 1,irmdnew!cellnew%nrmaxnew
- thetansll(lm1,lm2,ir) = thetansll(lm1,lm2,ir)+cleb(j)*shapefun_mod(ir,lm3)
- thetansll(lm2,lm1,ir) = thetansll(lm2,lm1,ir)+cleb(j)*shapefun_mod(ir,lm3)
- end do
- end if
+
+ ! This is effectivly a summation over one angular momentum index lm3 for each
+ ! combination of lm1, lm2. Iterate instead over the flattened array cleb
+ ! containing the Gaunt coefficients and extract the angular momentum
+ ! indices for each j.
+ do j = 1, iend
+ lm1 = icleb(j, 1)
+ lm2 = icleb(j, 2)
+ lm3 = icleb(j, 3)
+ if (lm1 > lmmax .or. lm2 > lmmax .or. lm3 > lmmax2) then
+ ! I think this should not happen, as icleb should only map to valid values
+ ! for lm1,lm2,lm3.
+ !TODO check that and remove this warning
+ write(*,*) "Warning from calc_thetallmat in bfield.F90: Maybe invalid values in icleb:"
+ write(*,*) j, lm1, lm2, lm3, lmmax, lmmax2
+ cycle ! Gaunt coeffs zero
+ end if
+ if (ifunm(lm3) == 0) then
+ cycle ! Some symmetry related to shapefunction, I think
end if
+ do ir = 1, irmdnew
+ thetansll(lm1,lm2,ir) = thetansll(lm1,lm2,ir) + cleb(j) * shapefun_mod(ir,lm3)
+ thetansll(lm2,lm1,ir) = thetansll(lm2,lm1,ir) + cleb(j) * shapefun_mod(ir,lm3)
+ end do
end do
- !do lm1 = 1,lmmax
- ! do lm2 = 1,lm1-1
- ! do ir = 1, irmdnew
- ! thetansll(ir,lm2,lm1) = thetansll(ir,lm1,lm2)
- ! end do
- ! end do
- !end do
- !if(myrank==master .and. debug) write(17*iatom**2,'(" thetansll = ")')
- !do lm1=1,lmmax
- ! do lm2=1,lmmax
- ! if(sum(abs(thetansll(:,lm1,lm2)))>1e-8) then
- ! if(myrank==master .and. debug) write(17*iatom**2,'(2i4,1000es16.8)') lm1, lm2, thetansll(:,lm1,lm2)
- ! end if
- ! end do
- !end do
-
+
end subroutine calc_thetallmat
end module mod_bfield
-
-
-
-
-
--
GitLab