Step 1.

Optimization geometry of ground electronic state of your molecule (ICN)

Step 2

Computation of (by TD-DFT) excitation energies assuming the geometry obtained in Step 1 and that molecule has one electron less than initial molecule (in your case it is ICN+). Please remember about the selection of correct multiplicity.

Then you have a list of excited states of ICN+.

For first ionization energy, you need

1) found last line containing 'SCF Done:' in computation for ICN

2) found first line containing 'SCF Done:' in computation for ICN+

The difference of energies given in those lines is ionization energies in Hartree.

Then you can add the energy of excited state for ICN+ for next ionization energies.

You need select proper functional and basis set for this computation. The simple B3LYP could give good results, but it could be not enough. Let consider CAM-B3LYP or another functional (lots of benchmark studies was done). Another issue is basis set. For reliable results, you need big basis set.

Anyway, the vertical approximation isn't the best choice. You can optimize (together with frequencies computations) the ground states of ICN  and ICN+. Then you can extract 0-0 transition energies (differences between energies from lines starting from 'Sum of electronic and zero-point Energies'. Then you can add energies of excited states from TD-DFT.

The ICN is a small molecule. Let consider more reliable methods like CCSD or CASPT2

Dear Marcin, thank you  very much for your comment.  It was really helpful. I appreciate it. I have few further queries.

Since I am studying the electronic spectra of the molecule and  electronic transitions are mostly of Franck-Condon type, vertical Ionzation energies are  better choice  for me over  the adiabetic ionization energies. Now  when I follow the second suggestion of optimizing ground states of both  ICN and ICN+ ,  taking the difference and adding consecutive excited states,  gives excellent matches with experimental VIEs.  But I guess this  method should calculate adiabetic ionization energies since they are associated with optimized geometries.   Could you please correct me where I am wrong ?

I have another question.  While calculating six doublets of ICN+,  I put  td=(nstates=6) in the route section . It calculates the ground state. Now it should have five excited states. Why  does it  calculate  excitation energies for  six excited states ?  Its really confusing to me.  Also among six doublets  1-2 and 4-5  are  degenerate (from molpro calculation). But from  TDDFT  output no two states appear to be degenerate.

For your last concern, yes , I have calculated the CASPT2 energies too ( in  Molpro).  I am yet to start EOM-CCSD calculations (in Gaussian)  which I think is pretty much like TDDFT. 

Thanks again.

About comparisons between vertical and adiabatic approach. Only adiabatic energies can be measured experimentally. In principle, adiabatic energies should be compared with thru 0-0 transition energies. So, if your experimental spectra contain nice vibronic structure, you should rather compare them with results of adiabatic computations. However, in a case of several spectra (especially in solvent), you can have only broad band. For this case, vertical approximation works quite well and you can compare results of computations with the position of the maximum of the band.

About method described previously. The predicted energies are not pure adiabatic energies. Only first ionization energy is calculated in a pure adiabatic way. Then, in practice, you calculated the vertical excitation energies of ICN+, and added them to adiabatic first ionization energy (transition between the ground state of ICN and ground state of ICN+). Anyway, the whole optimization of an excited state of cation can typically reduce energy by no more than few hundreds of meV. It is much less important than the difference between adiabatic and vertical first ionization potential.

About TD(nstates=n). This keyword informs Gaussian that you would like to compute energies of n excited states. The ground state energy was calculated by “pure” DFT method. So, ground state energy is not a part of TD computations. It is earlier (obviously required) step.

About degeneracy. As far, as I know, the Gaussian treats degenerated states as two separated states and print them separately. So, for example, my molecule has following the pattern of states: X (not degenerate), A(no degenerate), B(degenerated). In the output, state A appears as one states (state number 1), state B appears as two separated states (states 2 and 3). States 2 and 3 both have exactly the same energies.

Probably you saw several benchmark studies for TD end EOM computations (I’m not able to give the whole list of publications on this topic). For example, from my recently published benchmark studies for atoms, I showed that EOM-CCSD is much more precise than TD-DFT for each basis set. The best TD-DFT functionals barely reach the precision of CIS(D) computations. Anyway, the logic and Gaussian implementation of EOM and TD are similar.

Dear Marcin, thank you very much for the detailed information. It taught me a lot , provided I am very new to Gaussian. Now I am turning up to EOM-CCSD energies of  doublet states of ICN+. I want a relative study between them(TDDFT and EOM-CCSD). Hope I will get your help this time also, if needed . 

Thanks again.

I am puzzled by the Marcin's statement "About comparisons between vertical and adiabatic approach. Only adiabatic energies can be measured experimentally." My experience is that most experimental IEs (measured by PES of gaseous species...) are vertical energies. Most organic molecules have broad PES spectra and do not show clear vibrational structures. PES spectra of small molecules like HCN have probably clear vibrational structures that help extracting vibrational progressions and adiabatic energies.