It is possible to run excited states calculations with the quantum package.
To do this, set the n_states variable in the Determinants section to the
number of requested states.
The selection criterion will be the maximum of the selection criteria for each state.
To produce state-average natural orbitals, run
qp_run save_natorb file.ezfio
The MOs will be replaced, so the two-electron integrals and the wave function are invalidated as well.
If you are computing multiple states, it is often a good idea to have the s2_eig flag set to true. This will force the Davidson algorithm to choose only vectors with a value of S^2 equal to the expected_s2. Otherwise, you can have different spin states coming out in the diagonalization.
Also, you need to have a number of states in the Davidson algorithm (n_states_diag) which is greater than n_states.
Usually, it is good practice to use state-averaged MOs so that all states have MOs of comparable quality. For example, if you search for a singly excited state, you can use state-average natural orbitals of a preliminary CIS calculation.
The QP doesn't take account of the symmetry. For reasons due to numerical noise, excited states of different symmetries may enter in the calculation. Note that it is possible to make state-average calculation of states with different symmetries and/or different spin multiplicities.