anchor.diff:              norb = nbasis/2
           end if
   
  +        ! Set up basis functions used in calculation.
           allocate(basis_fns(nbasis), stat=ierr)
           call check_allocate('basis_fns', nbasis, ierr)
  -
  -        ! Set up basis functions.
  -        call init_basis_fns_read_in(norb, uhf, orbsym(active_basis_offset/rhf_fac+1:norb+active_basis_offset/rhf_fac), basis_fns) 
  +        call init_basis_fns_read_in(norb, uhf, orbsym, sp_eigv, sp_eigv_rank(1+active_basis_offset/rhf_fac:), basis_fns) 
  +
  +        deallocate(sp_eigv, stat=ierr)
  +        call check_deallocate('sp_eigv', ierr)
   
           ! Was a symmetry found for all basis functions?  If not, then we must
           ! turn symmetry off.
  -        if (minval(orbsym(:norb+active_basis_offset/rhf_fac)) < 0) then
  -            write (6,'(1X,a62)') 'Unconverged symmetry found.  Turning point group symmetry off.'
  +        if (minval(basis_fns(:)%sym) < 0) then
  +            if (parent) write (6,'(1X,a62)') 'Unconverged symmetry found.  Turning point group symmetry off.'
               forall (i=1:nbasis) basis_fns(i)%sym = 0
           end if
   
  
positive.code:                                     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
                                  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 if
                              end if
                          else if (all( (/ ii,jj,aa,bb /) > 0)) then
                              ! Have <ij|ab> involving active orbitals.
                              call store_two_body_int_mol(ii, jj, aa, bb, x, coulomb_integrals)
                          end if
                      end if
  
                  end if
  
negative.code:         if (size(basis_fns) /= size(all_basis_fns)) then
              ! We froze some orbitals...
              ! Print out entire original basis.
              if (parent) then
                  call write_basis_fn_header()
                  do i = 1, size(all_basis_fns)
                      call write_basis_fn(all_basis_fns(i), ind=i, new_line=.true.)
                  end do
                  write (6,'(/,1X,"Freezing...",/,1X,"Using complete active space: (",' &
                            //int_fmt(cas(1),0)//',",",'//int_fmt(cas(2),0)//',")",/)') cas
              end if
          end if
  
          do i = 1, size(all_basis_fns)
  