diff --git a/source/KKRnano/source/CalculationData_mod.F90 b/source/KKRnano/source/CalculationData_mod.F90
index e50f08335424ad85d0cf1df2dad958c4f6e77a47..abd9c75a9d44c948197cd80f2ab5227cbd37e50f 100644
--- a/source/KKRnano/source/CalculationData_mod.F90
+++ b/source/KKRnano/source/CalculationData_mod.F90
@@ -139,10 +139,9 @@ module CalculationData_mod
     allocate(self%jij_data_a(num_local_atoms)) 
     if(num_local_atoms > 1 .and. params%Jij) warn(6, "Jij work with max. 1 atom so far!") 
     
-    ! Allocate bfields for all atoms, makes it easier to read from disk.
-    ! The types itself are very small, the larger components will only be allocated for local atoms.
+    ! Allocate bfields only if used
     if (params%noncobfield) then
-      allocate(self%bfields(dims%naez))
+      allocate(self%bfields(num_local_atoms))
     end if
 
     allocate(self%atom_ids(num_local_atoms))
@@ -388,12 +387,10 @@ module CalculationData_mod
 
     if (params%noncobfield) then
       ! Initialize the noncolinear magnetic field. If present, read from disk
-      call load_bfields_from_disk(self%bfields, params%constr_field)
-      ! For the local atoms, initialize some fields
+      call load_bfields_from_disk(self%bfields, params%constr_field, dims%naez, self%atom_ids)
+      ! Initialize the 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(self%bfields(atom_id), dims%lmaxd, &
+        call init_bfield(self%bfields(ila), dims%lmaxd, &
                          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), &
diff --git a/source/KKRnano/source/ProcessKKRresults_mod.F90 b/source/KKRnano/source/ProcessKKRresults_mod.F90
index 2ca697c40f93733e1e25b7db853a85cba4ea1df6..7f82afb631ac225772923c9c2fe1611180148a92 100644
--- a/source/KKRnano/source/ProcessKKRresults_mod.F90
+++ b/source/KKRnano/source/ProcessKKRresults_mod.F90
@@ -522,7 +522,7 @@ module ProcessKKRresults_mod
                           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), &
+                          densities%muorb, densities%iemxd, params, calc%bfields(ila), &
                           get_muffin_tin_index(atomdata%chebmesh_ptr), iter)
 
       ! LDAU
diff --git a/source/KKRnano/source/ScatteringCalculation_mod.F90 b/source/KKRnano/source/ScatteringCalculation_mod.F90
index 4f53b7396a9d1e1ed94423900617dd5643eadf14..7e33e5c347278c7a211125f91e4ad6c190eafc05 100644
--- a/source/KKRnano/source/ScatteringCalculation_mod.F90
+++ b/source/KKRnano/source/ScatteringCalculation_mod.F90
@@ -272,7 +272,7 @@ 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), get_muffin_tin_index(atomdata%chebmesh_ptr), &
+                                    calc%bfields(ila), get_muffin_tin_index(atomdata%chebmesh_ptr), &
                                     iter, params%itbfield0, params%itbfield1, &
                                     params%noncobfield, params%constr_field, params%trans_bfield, params%mt_bfield)
                
diff --git a/source/KKRnano/source/bfield/bfield.f90 b/source/KKRnano/source/bfield/bfield.f90
index 297a62a61cf3e55089c908c98cb77082a75fc380..0e139ecbd95826daabca386a64a25f80eb1c0b43 100644
--- a/source/KKRnano/source/bfield/bfield.f90
+++ b/source/KKRnano/source/bfield/bfield.f90
@@ -43,20 +43,38 @@ contains
   !> Load the external noncollinear magnetic field (if a file is present) and
   !> the initial guess for the constraining fields (if constraint magnetism is
   !> used and a file is present) from disk.
-  !> This subroutine loads information for all atoms into an array of structs.
-  subroutine load_bfields_from_disk(bfields, lbfield_constr)
-    type(bfield_data), intent(inout) :: bfields(:)
-    logical,           intent(in)    :: lbfield_constr
+  !> This subroutine loads information for the locally treated atoms into an
+  !> array of structs.
+  subroutine load_bfields_from_disk(bfields, lbfield_constr, naez, atom_ids)
+    type(bfield_data), intent(inout) :: bfields(:)        !! Array for the read fields
+    logical,           intent(in)    :: lbfield_constr    !! Wheter constraint magnetism is used
+    integer,                       intent(in) :: naez     !! Number of atoms in the unit cell
+    integer(kind=4), dimension(:), intent(in) :: atom_ids !! Locally treated atoms
     
-    integer :: number_of_atoms
+    integer :: num_local_atoms ! number of atoms treated in this MPI rank
+    integer :: ila             ! atom index
+    type(bfield_data), dimension(:), allocatable :: all_bfields ! array of all bfields to read
 
-    number_of_atoms = size(bfields)
-    
-    call read_bfield(bfields, number_of_atoms)
+    ! Get the number of locally treated atoms form the size of the output array
+    num_local_atoms = size(bfields)
+
+    ! Allocate an array for all atoms in the unit cell to read the files (which of course
+    ! contain the fields of every atom)
+    allocate(all_bfields(naez))
+
+    ! Read the files
+    call read_bfield(all_bfields)
     if (lbfield_constr) then
-      call read_bconstr(bfields, number_of_atoms)
+      call read_bconstr(all_bfields)
     end if
-    
+
+    ! Copy the locally treated atoms to the output
+    do ila = 1, num_local_atoms
+      bfields(ila) = all_bfields(atom_ids(ila))
+    end do
+
+    ! Forget the other atoms
+    deallocate(all_bfields)
   end subroutine
 
   !> Initialize a bfield_data struct. Allocates and calculates the shapefunction
@@ -99,14 +117,16 @@ contains
   !> 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)
+  subroutine read_bconstr(bfields)
     type(bfield_data), dimension(:), intent(inout) :: bfields
-    integer                        , intent(in)  :: number_of_atoms
     
+    integer :: number_of_atoms
     integer :: iatom
     integer :: iostat
     logical :: file_exists
 
+    number_of_atoms = size(bfields)
+
     inquire(file='bconstr_in.dat', exist=file_exists)
 
     if (file_exists) then
@@ -146,14 +166,16 @@ contains
   !> 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)
+  subroutine read_bfield(bfields)
     type(bfield_data), dimension(:), intent(inout) :: bfields
-    integer                        , intent(in)   :: number_of_atoms
-    
+
+    integer        :: number_of_atoms
     integer        :: iatom, iostat
     character(256) :: linebuffer
     logical        :: file_exists
 
+    number_of_atoms = size(bfields)
+
     inquire(file='bfield.dat', exist=file_exists)
    
     if (file_exists) then
@@ -200,6 +222,7 @@ contains
    end do
   end subroutine read_bfield
 
+
   !>------------------------------------------------------------------------------
   !> Summary: Adds the magnetic field to the potential
   !> Author: Sascha Brinker, Nicolas Essing