One last question: In case I would like to use Orca, should I generate the .molden file (using orca_2mkl) from the .gbw of a tddft calculation and then using Multiwfn (options 18,13) generate the natural orbitals file?

Thank you again,
Lorenzo

In your way the wavefunction of excited state cannot be directly passed to Multiwfn. Although density matrix of excited state was written into .chk file due to "density=current", Multiwfn doesn't directly utilize it because by default Multiwfn only loads orbital information from the inputted wavefunction file.

The most convenient way of performing this analysis for an excited state in Multiwfn is using "out=wfn" keyword and specify the .wfn file to be exported at the end of Gaussian input file. After TDDFT calculation, the natural orbitals of the excited state chosen by "root" (default is 1) will be generated and exported to the .wfn file, which can be used as input file of Multiwfn.

By the way, using "nstates=1" is deprecated. Due to the TDDFT solving algorithm reason, when you are interest in excited state i, "nstates" should be set to at least i+2 to guarantee enough solving accuracy for state i.

Everything runs fine and I obtain a value for the coefficent quite similar to Amber GAFF (using Hirshfeld-I partitioning and atmrad folder) using the .fchk file for the molecule and the .wfn files for carbon and hydrogen. Then, I wanted to check if excited states could also be dealt with. My input file for the excited state is:

Following the exact same procedure (substituting the .fchk file with the excited state one and keeping the same .wfn files for the isolated atoms as in the ground state), I obtain identical parameters for the C6 coefficients as if the program was reading the ground state density from the .fchk instead of the excited state density.

Could you please give me some insights abouth this?
Thank you very much for your help,
Lorenzo