anchor.diff:              if (ios == iostat_end) exit ! reached end of file
               if (ios /= 0) call stop_all('read_input','Problem reading input.')
   
  -            if (.not.uhf) then
  -                ! Working in spin-orbitals but FCIDUMP is in spatial orbitals.
  -                ! Need to only store integrals in one spin-channel in RHF, so
  -                ! can just with (e.g.) beta orbitals.
  -                i = i*2
  -                j = j*2
  -                a = a*2
  -                b = b*2
  -            end if
  -
  -            if (i > 0 .and. j == 0 .and. a == 0 .and. b == 0) then
  -                ! \epsilon_i --- temporarily store for all basis functions,
  -                ! including inactive (frozen) orbitals.
  -                if (uhf) then
  -                    all_basis_fns(i)%sp_eigv = x
  -                else
  -                    all_basis_fns(i-1)%sp_eigv = x
  -                    all_basis_fns(i)%sp_eigv = x
  -                end if
  -            end if
  -
  -            ! Adjust indices to take into account frozen core orbitals.
  +            ! Working in spin orbitals but FCIDUMP is in spatial orbitals in RHF
  +            ! calculations and spin orbitals in UHF calculations, and te basis
  +            ! is ! not necessarily ordered by energy.  We wish to work in spin
  +            ! orbitals ordered by energy.
  +            ! Need to only store integrals in one spin-channel in RHF, so
  +            ! can just get away with referring (e.g.) beta orbitals, hence the
  +            ! factor of 2 in RHF.  The integral store routines convert the
  +            ! indices further to compress the integral stores as much as
  +            ! possible.
  +            i = rhf_fac*sp_fcidump_rank(i)
  +            j = rhf_fac*sp_fcidump_rank(j)
  +            a = rhf_fac*sp_fcidump_rank(a)
  +            b = rhf_fac*sp_fcidump_rank(b)
  +
  +            ! Adjust indices to take into account frozen core orbitals and to
  +            ! convert to an energy ordering.
               ii = i - active_basis_offset
               jj = j - active_basis_offset
               aa = a - active_basis_offset
  
positive.code:                     ! Nuclear energy.
                      Ecore = Ecore + x
  
                  else if (ii > 0 .and. j == 0 .and. a == 0 .and. b == 0) then
  
                      ! \epsilon_i
                      if (uhf) then
                          basis_fns(ii)%sp_eigv = x
                      else
                          basis_fns(ii-1)%sp_eigv = x
                          basis_fns(ii)%sp_eigv = x
                      end if
  
                  else if (j == 0 .and. b == 0) then
  
  
negative.code:         deallocate(all_basis_fns, stat=ierr)
          call check_deallocate('all_basis_fns',ierr)
  
          if (parent) then
              call write_basis_fn_header()
              do i = 1, nbasis
                  call write_basis_fn(basis_fns(i), ind=i, new_line=.true.)
              end do
              write (6,'(/,1X,a8,f18.12)') 'E_core =', Ecore
          end if
  
          if (.not.t_store) then
              ! Should tidy up and deallocate everything we allocated.
  