subroutine write_spat_doub_occ_stats(initial)
! routine to write out the double occupancy data
logical, intent(in), optional :: initial
character(*), parameter :: this_routine = "write_spat_doub_occ_stats"
type(write_state_t), save :: state
logical, save :: inited = .false.
integer :: i
character(12) :: num
def_default(state%init, initial, .false.)
! If the output file hasn't been opened yet, then create it.
if (iProcIndex == Root .and. .not. inited) then
state%funit = get_free_unit()
call init_spat_doub_occ_stats(state%funit)
inited = .true.
end if
if (iProcIndex == root) then
if (state%init .or. state%prepend) then
write(state%funit, '("#")', advance='no')
state%prepend = state%init
else if (.not. state%prepend) then
write(state%funit, '(" ")', advance='no')
end if
state%cols = 0
state%cols_mc = 0
state%mc_out = tMCOutput
call stats_out(state, .false., iter + PreviousCycles, 'Iter.')
do i = 1, nBasis / 2
write(num, '(i12)') i
call stats_out(state, .false., all_inst_spatial_doub_occ(i) / &
(sum(all_norm_psi_squared) / real(inum_runs, dp) * real(StepsSft, dp)), &
'orbital '//trim(adjustl(num)))
end do
! And we are done
write(state%funit, *)
call neci_flush(state%funit)
end if
end subroutine write_spat_doub_occ_stats