```
subroutine FindNatOrbs()
! Fed into this routine will be the wavefunction and its amplitudes
! within the given excitation level.
! First need to set up the orbital labels and symmetries etc. This is
! done slightly differently for spin and spatial and whether or not we
! are truncating the virtual space when writing out the final ROFCIDUMP
! file.
! Allocate the matrix used to find the natural orbitals.
integer :: ierr
character(len=*), parameter :: t_r = 'FindNatOrbs'
allocate(NatOrbMat(NoOrbs, NoOrbs), stat=ierr)
call LogMemAlloc('NatOrbMat', NoOrbs**2, 8, t_r, NatOrbMatTag, ierr)
if (ierr /= 0) call stop_all(t_r, "Mem allocation for NatOrbMat failed.")
NatOrbMat(:, :) = 0.0_dp
allocate(Evalues(NoOrbs), stat=ierr)
call LogMemAlloc('Evalues', NoOrbs, 8, t_r, EvaluesTag, ierr)
if (ierr /= 0) call stop_all(t_r, "Mem allocation for Evalues failed.")
Evalues(:) = 0.0_dp
! First need to fill the relevant matrix for calculating the type of
! natural orbitals we want.
if (tFindCINatOrbs) then
! For the CISD, CISDT etc natural orbitals, the relevant matrix is
! the one electron reduced density matrix from the previous
! spawning calculation (trucated at a certain excitation).
call FillOneRDM()
else if (tUseMP2VarDenMat) then
! For the MP2 natural orbitals, the natural orbital matrix is the
! MP2 variational density matrix.
call FillMP2VDM()
end if
! Then need to diagonalise this, maintaining the various symmetries
! (spin and spatial).
call DiagNatOrbMat()
! We then need to put the resulting eigenvectors back into the ordering
! we want, and copy these over to CoeffT1.
call OrderCoeffT1()
end subroutine FindNatOrbs
```

NECI was developed by Alavi Group

© 2023

Documentation generated by FORD