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

@@ -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

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

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)

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)