Dear All,
Could you help me in the following - I’m interested in excited states of diatomic molecules and for that I need to translate molecular term (which I will usually know) into Psi4 input.

I have tested successfully 3 electronic states of NH molecule just specifying multiplicity (3Sigma - ground state, and 5Sigma and 1Delta are excited states). In general there will be more states with the same multiplicity - I guess I should know how to translate term symbol into Psi4 input.

My code was:

molecule dimer {units bohr
0 5 % (multiplicity which was the only thing I’ve changed)
N
H 1 R
}
Rvals=[2 (more points in general)]
set basis aug-cc-pVQZ
set reference uhf
set maxiter 500
ecp={}
for R in Rvals:
dimer.R=R
ecp[R]=energy(‘ccsd(t)’)

for R in Rvals:
e=ecp[R]
psi4.print_out(" %3.1f %10.12f\n" % (R,e))

Yes, I do. I have full term symbol and I want to make sure I’m calculating this exact state not any other. I think I should also specify dissociation limit of molecule because there can be more then 1 state with exactly the same term symbol.

Most electronic structure methods, CCSD(T) among them, are adapted for computing the ground state of a given electronic symmetry and cannot compute excited states (for a given symmetry) effectively. There are electronic structure methods that can compute excited states, such as EOM-CC, equation-of-motion coupled cluster. Psi4 also has some ADC and TD-DFT capabilities if you want excitation energies at given geometries. I can talk more about those if you tell me what you need.

Now, once you’ve decided on a term symbol, how do you compute its lowest energy state? Psi4, along with most other programs, supports real orbitals and abelian symmetry groups only. That means that you’ll need to correlate orbital occupations from the full C-infinity-v symmetry down to C2v symmetry. From there, you can specify occupations with the DOCC and SOCC keywords. See here.

Let me know if I need to go into more detail on anything.

I have tested 5Sigma state with EOM-CC3. My input was:
olecule dimer {units bohr
N
H 1 R
}
Rvals=[2]
set basis aug-cc-pVQZ
set reference uhf
set maxiter 100
set DOCC [2,0,0,0]
set SOCC [2,0,1,1]
set roots_per_irrep [1,1,1,1]
ecp={}
for R in Rvals:
dimer.R=R
ecp[R]=energy(‘eom-cc3’)

for R in Rvals:
e=ecp[R]
psi4.print_out(" %3.1f %10.12f\n" % (R,e))
On the basis of molecular orbitals I have established DOCC and SOCC and set roots_per_irrep [1,1,1,1] (lets say I don’t want to bother about resulting symmetry). I also did not set multiplicity (5) and I got correct CC3 energy -54.874 (in agreement with my previous results without EOM; atomic units).

Now my question is how to interpret the output - there is CC3 value, and there are also EOM energies:
A1 -54.609
A2 -54.790
B1 -54.775
B2 -53.775
What is the relation of those energies? I expected that one of symmetries (A1,A2,B1,B2) will be the value I expected. (I noticed that they are all higher as if they were further excitations). Does DOCC and SOCC values are for all those energies.

To sum up, I am confused because I expected that there will be no cc3 energy but only the eom-cc3 energies. Please explain explain me the logic of those energies.

I am attaching also output file. NHe.txt (53.4 KB)

Final Energetic Summary for Converged Roots of Irrep A1
Excitation Energy Total Energy
(eV) (cm^-1) (au) (au)
EOM State 2 7.218 58217.9 0.2652601280 -54.609369851622

Psi reports both the excitation energy and the total energy obtained by adding the excitation energy to the ground state’s energy.

Do I understand correctly that I have to specify the ground state (by DOCC, SOCC and spin multiplicity) and get energies of states with various symmerties (according to root_per_irrep)?

In the meantime I run calculation like (input file specified the ground state) that and got for 5Sigma eom-cc3 A2 energy -54.863 vs. -54.873 from cc3 or ccsd(t). But I understood your initial explanation that I have to specify the exited state I want with DOCC and SOCC.

So my input file should specify the ground state and I get excited states energies (now my calculations supports that) OR it should specify excited state I am aiming at? (then most probably I have a problem; I am assuming in this reasoning that my previous energy calculated just by setting multiplicity 5 is correct)