From faa0629e8f405617f43dada9c4862cf9f41e062c Mon Sep 17 00:00:00 2001
From: Nicolas Essing <n.essing@fz-juelich.de>
Date: Thu, 14 Oct 2021 13:48:24 +0200
Subject: [PATCH] Fixed a bug in constraint bfield calculation related to
 muffin tins.

The muffin tin index was taken from the old mesh, which is not correct
and can even give errors if it is larger than the total size in the
new mesh.

Introduced a function in the Chebyshev mesh module that calculates
the index of the muffin tin radius. Formula for this is taken from
KKRhost.
Used this function where before the muffin tin index of the old mesh
was used. Added the dependency on this function to the Makefile.
---
 source/KKRnano/source/Makefile                           | 2 +-
 source/KKRnano/source/ProcessKKRresults_mod.F90          | 8 +++++---
 source/KKRnano/source/ScatteringCalculation_mod.F90      | 5 +++--
 .../KKRnano/source/datastructures/ChebMeshData_mod.F90   | 9 ++++++++-
 4 files changed, 17 insertions(+), 7 deletions(-)

diff --git a/source/KKRnano/source/Makefile b/source/KKRnano/source/Makefile
index 20f27ba95..5895e07c0 100644
--- a/source/KKRnano/source/Makefile
+++ b/source/KKRnano/source/Makefile
@@ -466,7 +466,7 @@ endif
 KKRnano.o: Logging_mod.o KKRzero_mod.o PotentialConverter_mod.o KKRnanoParallel_mod.o BasisAtom_mod.o AtomicCore_mod.o RadialMeshData_mod.o main2_aux_mod.o ScatteringCalculation_mod.o Main2Arrays_mod.o KKRnano_Comm_mod.o ProcessKKRresults_mod.o InputParams_mod.o EBalanceHandler_mod.o LDAUData_mod.o TimerMpi_mod.o EnergyMesh_mod.o DimParams_mod.o CalculationData_mod.o
 TruncationZone_mod.o: Main2Arrays_mod.o
 KKRmat_mod.o: ExchangeTable_mod.o two_sided_comm_mod.o arraytest2_mod.o Logging_mod.o TimerMpi_mod.o KKROperator_mod.o one_sided_comm_mod.o ClusterInfo_mod.o IterativeSolver_mod.o jij_calc_mod.o Truncation_mod.o SparseMatrixDescription_mod.o BCPOperator_mod.o SolverStats_mod.o InitialGuess_mod.o fillKKRMatrix_mod.o DirectSolver_mod.o
-ProcessKKRresults_mod.o: arraytest2_mod.o KKRnanoParallel_mod.o AtomicForce_mod.o Logging_mod.o Lloyds_formula_mod.o wrappers_mod.o brydbm_new_com_mod.o Main2Arrays_mod.o ShapefunData_mod.o EnergyResults_mod.o LDAUData_mod.o ShapeGauntCoefficients_mod.o GauntCoefficients_mod.o CalculationData_mod.o muffin_tin_zero_mod.o RadialMeshData_mod.o MadelungPotential_mod.o broyden_kkr_mod.o EnergyMesh_mod.o BasisAtom_mod.o DensityResults_mod.o Lloyd0_new_mod.o KKRresults_mod.o InputParams_mod.o TimerMpi_mod.o DimParams_mod.o BroydenData_mod.o NearField_calc_mod.o debug_morgan_mod.o total_energy_mod.o
+ProcessKKRresults_mod.o: arraytest2_mod.o KKRnanoParallel_mod.o AtomicForce_mod.o Logging_mod.o Lloyds_formula_mod.o wrappers_mod.o brydbm_new_com_mod.o Main2Arrays_mod.o ShapefunData_mod.o EnergyResults_mod.o LDAUData_mod.o ShapeGauntCoefficients_mod.o GauntCoefficients_mod.o CalculationData_mod.o muffin_tin_zero_mod.o RadialMeshData_mod.o MadelungPotential_mod.o broyden_kkr_mod.o EnergyMesh_mod.o BasisAtom_mod.o DensityResults_mod.o Lloyd0_new_mod.o KKRresults_mod.o InputParams_mod.o TimerMpi_mod.o DimParams_mod.o BroydenData_mod.o NearField_calc_mod.o debug_morgan_mod.o total_energy_mod.o ChebMeshData_mod.o
 NearField_kkr_mod.o: Constants_mod.o
 NearField_calc_mod.o: RadialMeshData_mod.o NearField_com_mod.o BasisAtom_mod.o CalculationData_mod.o
 NearField_mod.o: Harmonics_mod.o MadelungCalculator_mod.o Constants_mod.o NearField_kkr_mod.o
diff --git a/source/KKRnano/source/ProcessKKRresults_mod.F90 b/source/KKRnano/source/ProcessKKRresults_mod.F90
index e54d45425..2ca697c40 100644
--- a/source/KKRnano/source/ProcessKKRresults_mod.F90
+++ b/source/KKRnano/source/ProcessKKRresults_mod.F90
@@ -420,6 +420,7 @@ module ProcessKKRresults_mod
     use KKRresults_mod, only: KKRresults
     use DensityResults_mod, only: DensityResults
     use EnergyResults_mod, only: EnergyResults
+    use ChebMeshData_mod, only: get_muffin_tin_index
     
     use lloyd0_new_mod, only: lloyd0_wrapper_com
     use wrappers_mod, only: rhoval_wrapper, RHOTOTB_wrapper, RHOMOM_NEW_wrapper
@@ -519,9 +520,10 @@ module ProcessKKRresults_mod
                           calc%gaunts, emesh, ldau_data, params%Volterra, &
                           dims%korbit, calc%noco_data%theta_noco(atom_id), calc%noco_data%phi_noco(atom_id), &
                           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, calc%bfields(atom_id), calc%mesh_a(ila)%imt, iter)
+                          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, calc%bfields(atom_id), &
+                          get_muffin_tin_index(atomdata%chebmesh_ptr), 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 101835862..4f53b7396 100644
--- a/source/KKRnano/source/ScatteringCalculation_mod.F90
+++ b/source/KKRnano/source/ScatteringCalculation_mod.F90
@@ -74,7 +74,7 @@ implicit none
     
     use two_sided_commD_mod, only: distribute
     
-    use ChebMeshData_mod, only: interpolate_poten  ! NOCO
+    use ChebMeshData_mod, only: interpolate_poten, get_muffin_tin_index  ! NOCO
     use NonCollinearMagnetism_mod, only: tmat_newsolver  ! NOCO
     use NonCollinearMagnetism_Helpers_mod, only: rotatematrix  ! NOCO
 
@@ -272,7 +272,8 @@ implicit none
                                     !dims%lly,        &    
                                     atomdata%potential%lmpot,atomdata%chebmesh_ptr%irmd_new, &
                                     kkr(ila)%TmatN(:,:,ispin),params%soc,params%enable_quad_prec, &
-                                    calc%bfields(i1), calc%mesh_a(ila)%imt, iter, params%itbfield0, params%itbfield1, &
+                                    calc%bfields(i1), get_muffin_tin_index(atomdata%chebmesh_ptr), &
+                                    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)
diff --git a/source/KKRnano/source/datastructures/ChebMeshData_mod.F90 b/source/KKRnano/source/datastructures/ChebMeshData_mod.F90
index e16f623bb..e6d9071e5 100644
--- a/source/KKRnano/source/datastructures/ChebMeshData_mod.F90
+++ b/source/KKRnano/source/datastructures/ChebMeshData_mod.F90
@@ -17,7 +17,7 @@ module ChebMeshData_mod
   public :: writeChebMeshDataDA           
   public :: readChebMeshDataDA, openChebMeshDataDAFile, closeChebMeshDataDAFile, writeChebMeshDataIndexDA      
   public :: readChebMeshDataIndexDA, openChebMeshDataIndexDAFile, closeChebMeshDataIndexDAFile 
-  public :: interpolate_poten  
+  public :: interpolate_poten, get_muffin_tin_index
  ! public :: readChebMeshDataHeader    
   
   type ChebMeshData
@@ -533,6 +533,13 @@ ipot=1
   end do ! ispin
 end subroutine interpolate_poten
 
+!> Get the index of the muffin tin
+integer function get_muffin_tin_index(chebmesh) result(imt)
+  type(ChebMeshData), intent(in) :: chebmesh
+
+  imt = chebmesh%ipan_intervall(chebmesh%npan_lognew + chebmesh%npan_eqnew) + 1
+end function
+
   !----------------------------------------------------------------------------
   !> Helping routine for ConstructChebMesh 
 subroutine chebmesh(npan,ncheb,ri,ro)
-- 
GitLab