diff --git a/src/initialisation_module.f90 b/src/initialisation_module.f90
index 66223a82..c99c1624 100644
--- a/src/initialisation_module.f90
+++ b/src/initialisation_module.f90
@@ -335,7 +335,8 @@ subroutine set_up(find_chdens,level)
                                       nspin, min_layer,                &
                                       glob2node, flag_XLBOMD,          &
                                       flag_neutral_atom, flag_diagonalisation, &
-                                      flag_write_occ_mat, occ_mat, occ_mat_glob
+                                      flag_write_occ_mat, occ_mat, occ_mat_glob, &
+                                      flag_DFTplusU
     use memory_module,          only: reg_alloc_mem, reg_dealloc_mem,  &
                                       type_dbl, type_int
     use group_module,           only: parts
@@ -490,8 +491,14 @@ subroutine set_up(find_chdens,level)
     if (inode == ionode .and. iprint_init > 2) &
          write (io_lun,fmt='(4x,a)') trim(prefix)//'Made covering set for matrix multiplications'
 
-    allocate(occ_mat(7,7,bundle%n_prim,nspin))
-    if(flag_write_occ_mat) allocate(occ_mat_glob(7,7,ni_in_cell,nspin))
+    if(flag_DFTplusU) then
+       allocate(occ_mat(7,7,bundle%n_prim,nspin))
+       occ_mat = zero
+       if(flag_write_occ_mat) then
+          allocate(occ_mat_glob(7,7,ni_in_cell,nspin))
+          occ_mat_glob = zero
+       end if
+    end if
     ! Create all of the indexing required to perform matrix multiplications
     ! at a later point. This routine also identifies all the density
     ! matrix range interactions and hamiltonian range interactions
diff --git a/src/move_atoms.module.f90 b/src/move_atoms.module.f90
index 764104b1..aa90f51d 100644
--- a/src/move_atoms.module.f90
+++ b/src/move_atoms.module.f90
@@ -3661,7 +3661,10 @@ subroutine update_H(fixed_potential)
                                       flag_SFcoeffReuse, flag_diagonalisation, &
                                       ne_spin_in_cell,                 &
                                       ne_in_cell, spin_factor,         &
-                                      ni_in_cell, area_moveatoms
+                                      ni_in_cell, area_moveatoms, &
+                                      flag_write_occ_mat, occ_mat, &
+                                      flag_DFTplusU
+    use primary_module,         only: bundle
     use density_module,         only: set_atomic_density,              &
                                       density, set_density_pcc,        &
                                       get_electronic_density
@@ -3685,6 +3688,12 @@ subroutine update_H(fixed_potential)
     integer :: spin_SF, spin, stat
 
     call start_timer(tmr_l_tmp1,WITH_LEVEL)
+    ! Update DFT+U occupation matrix
+    if(flag_DFTplusU) then
+       deallocate(occ_mat)
+       allocate(occ_mat(7,7,bundle%n_prim,nspin))
+       occ_mat = zero
+    end if
     ! (0) Pseudopotentials: choose correct form
     select case (pseudo_type)
     case (OLDPS)
diff --git a/src/pao2blip.f90 b/src/pao2blip.f90
index 57ae0651..00ce0e16 100644
--- a/src/pao2blip.f90
+++ b/src/pao2blip.f90
@@ -806,7 +806,7 @@ subroutine blips_in_star(inode,ionode,sym_type,&
        scal_prod_sym(n) = scal_prod(m)/coeff(n)
     end do
     if((inode == ionode).and.(iprint_basis >= 4)) then
-       write(unit=io_lun,fmt='(/6x" symmetric scalar product:"/)')
+       write(unit=io_lun,fmt='(/6x," symmetric scalar product:"/)')
        do n = 1, n_star_in_sphere
           write(unit=io_lun,fmt='(6x,3x,i5,3x,e15.6)') n, scal_prod_sym(n)
        end do