anchor.diff:                                      seen_ijij(tri_ind_reorder(i,j)) = seen_ijij(tri_ind_reorder(i,j)) + 2
                                   end if
                               end if
  -                        else if (min(ii,jj,aa,bb) <= 0) then
  -                            ! Have an integral involving some core and some
  -                            ! active orbitals.
  -                            if (ii == aa .and. min(jj,bb) >= 1) then
  -                                ! Update <j|h|b> with contribution <ij|ib>.
  -                                if (mod(seen_iaib(i, tri_ind_reorder(jj,bb)),2) == 0) then
  -                                    x = get_one_body_int_mol(one_e_h_integrals, jj, bb) + x*rhf_fac
  -                                    call store_one_body_int_mol(jj, bb, x, one_e_h_integrals)
  -                                    seen_iaib(i, tri_ind_reorder(jj,bb)) = seen_iaib(i, tri_ind_reorder(jj,bb)) + 1
  +                        case(2)
  +                            ! Have an integral involving two core and two active orbitals.
  +                            ic = 1
  +                            ia = 1
  +                            do iorb = 1, 4
  +                                if (orbs(iorb) > 0) then
  +                                    active(ia) = orbs(iorb)
  +                                    ia = ia +1
  +                                else
  +                                    core(ic) = orbs(iorb)
  +                                    ic = ic + 1
                                   end if
  -                            else if (jj == bb .and. min(ii,aa) >= 1) then
  -                                ! Update <i|h|a> with contribution <ij|aj>
  -                                if (mod(seen_iaib(j, tri_ind_reorder(ii,aa)),2) == 0) then
  -                                    x = get_one_body_int_mol(one_e_h_integrals, ii, aa) + x*rhf_fac
  -                                    call store_one_body_int_mol(ii, aa, x, one_e_h_integrals)
  -                                    seen_iaib(j, tri_ind_reorder(ii,aa)) = seen_iaib(j, tri_ind_reorder(ii,aa)) + 1
  -                                end if
  -                            else if (ii == bb .and. min(jj,aa) >= 1) then
  -                                ! Update <j|h|a> with contribution <ij|ai>.
  -                                if (seen_iaib(i, tri_ind_reorder(jj,aa)) < 2) then
  -                                    x = get_one_body_int_mol(one_e_h_integrals, jj, aa) - x
  -                                    call store_one_body_int_mol(jj, aa, x, one_e_h_integrals)
  -                                    seen_iaib(i, tri_ind_reorder(jj,aa)) = seen_iaib(j, tri_ind_reorder(jj,aa)) + 2
  -                                end if
  -                            else if (jj == aa .and. min(ii,bb) >= 1) then
  -                                ! Update <i|h|b> with contribution <ij|jb>
  -                                if (seen_iaib(j, tri_ind_reorder(ii,bb)) < 2) then
  -                                    x = get_one_body_int_mol(one_e_h_integrals, ii, bb) - x
  -                                    call store_one_body_int_mol(ii, bb, x, one_e_h_integrals)
  -                                    seen_iaib(j, tri_ind_reorder(ii,bb)) = seen_iaib(j, tri_ind_reorder(ii,bb)) + 2
  +                            end do
  +                            if (core(1) == core(2)) then
  +                                ! Have integral of type < i a | i b > or < i a | b i >,
  +                                ! where i is a core orbital and a and b are
  +                                ! active orbitals.
  +                                if ((ii == core(1) .and. aa == core(1)) .or. (jj == core(1) .and. bb == core(1))) then
  +                                    ! < i a | i b >
  +                                    if (mod(seen_iaib(core(1), tri_ind_reorder(active(1),active(2))),2) == 0) then
  +                                        ! Update <a|h|b> with contribution <ia|ib>.
  +                                        x = get_one_body_int_mol(one_e_h_integrals, active(1), active(2)) + x*rhf_fac
  +                                        call store_one_body_int_mol(active(1), active(2), x, one_e_h_integrals)
  +                                        seen_iaib(core(1), tri_ind_reorder(active(1),active(2))) = &
  +                                            seen_iaib(core(1), tri_ind_reorder(active(1),active(2))) + 1
  +                                    end if
  +                                else
  +                                    ! < i a | b i > (or a permutation thereof)
  +                                    if (seen_iaib(core(1), tri_ind_reorder(active(1),active(2))) < 2) then
  +                                        ! Update <j|h|a> with contribution <ij|ai>.
  +                                        x = get_one_body_int_mol(one_e_h_integrals, active(1), active(2)) - x
  +                                        call store_one_body_int_mol(active(1), active(2), x, one_e_h_integrals)
  +                                        seen_iaib(core(1), tri_ind_reorder(active(1),active(2))) = &
  +                                            seen_iaib(core(1), tri_ind_reorder(active(1),active(2))) + 2
  +                                    end if
                                   end if
                               end if
  -                        else if (all( (/ ii,jj,aa,bb /) > 0)) then
  +                        case(4)
                               ! Have <ij|ab> involving active orbitals.
                               call store_two_body_int_mol(ii, jj, aa, bb, x, coulomb_integrals)
  -                        end if
  +                        end select
                       end if
   
                   end if
  
positive.code: 
      end subroutine read_in_fcidump
  
      subroutine init_basis_fns_read_in(norb, uhf, orbsym, basis_arr)
  
          ! Initialise list of basis functions based upon the FCI namelist`data
          ! read in from an FCIDUMP file.
  
          ! In:
          !    norb: number of orbitals/functions in the FCIDUMP file.
          !    uhf: true if FCIDUMP file was produced from an unrestricted (UHF)
          !         calculation.  norb refers to spatial orbitals in a RHF-based
          !         FCIDUMP file and spin orbitals in a UHF-based FCIDUMP file.
          !    orbsym: list of symmetries of each orbital/function.
          ! In/Out:
          !    basis_arr: array of basis functions.  On entry it is allocated to
          !    (at least) the size of the basis.  On exit the individual basis_fn
  
negative.code:             rhf_fac = 2  ! need to double count some integrals.
          end if
  
          if (present(cas_info)) then
              if (any(cas_info < 1)) then
                  cas = (/ nel, nbasis/2 /)
  