Comparing EE and EA for EOM-XX-CCSD and open-shell doublets

Questions
#1

Hello,
My natural inclination for looking at (geometry optimizations and excitation energies) excited states of small molecule clusters (eg CN-N2) was to use EOM-EE-CCSD. After a closer reading of the manual and some of Dr. Krylov’s papers, I wanted to compare EE to EA, as the latter appears to perform well for Rydberg excited states of open-shell doublets.

In attempting to educate myself, I began with the example in the manual using CN and a couple questions came up. When choosing the occupation, is the EOM_FAKE_IPEA orbital easiest to find by looking through the MOs from the initial HF and finding the position of the 0 eV one (ie 14 in this case)?

Secondly, when comparing the two methods, all of the excitation energies for EA are negative. Am I correct in stating that the first state listed is X (-13.9355 eV) and the first excited state is the second one (-7.5946 eV)?

While the above diatomic has a high level of symmetry, most of the systems I am interested in would likely fall in the C1 point group.
Thanks,
Benj

#2

Benj,

  1. You do not need to use EOM-FAKE-IPEA – this is an old feature which is now obsolete because we have proper EOM-IP and EOM-EA codes. They could be deployed by using
    IP_STATES and EA_STATES, along with the METHOD=EOM-CCSD keyword. This should be executed with ccman2=true (default value).
  2. Doublet CN radical is a perfect EOM-IP system (described by EOM-IP from the anionic reference) – you can describe both 2Pi and 2Sigma states very nicely.
  3. However, if you are interested in higher excited states, including excitations from N2, you may need to consider another EOM method. EOM-EA (from the CN+ reference) can be performed, but you will not be able to describe both 2Pi and 2Sigma states of CN, only one of them, depending which cationic state you take as a reference.
  4. In EOM-IP/EA calculations, what is reported as ‘excitation energies’ are energy differences wrt to the reference state. So your ground state would be your lowest EOM-EA state and your excitation energies should be computed relative to it. Negative values of EOM-EA transition energies mean that your EOM-EA state is below your reference state, which makes perfect sense in this example.
  5. You may find this tutorial useful:
    http://iopenshell.usc.edu/pubs/pdf/revcc-30.4-151.pdf

Anna.

#3

Thank you very much Dr. Krylov!

Both your answers above the the referenced article were helpful in understanding the EA/IP versions of EOM-XX-CCSD in Q-Chem. I am only interested in the first few excited states, which should be localized on the CN (or NO below).

Just to make sure I am on the right track, I also looked at NO. From the info you gave me, I think this would be a good candidate for EOM-EA b/c the pi_u 2p orbitals are filled for NO+, leaving the pi_g 2p orbitals completely empty. EOM-IP would not be suitable in this case because adding an electron to begin with NO- would result in half filled pi_g orbitals which is not ideal.
Thanks,
Benj

#4

Sorry for the extra post Dr. Krylov, another question popped up as I began evaluating triples corrections today. It appears (fT) is not quite working for EOM-EA/IP at the moment. What are your thoughts on N2-CN and N2-NO using EOM-EE-CCSD(fT) instead of EA/IP so I can include non-iterative triples?

#5

Well, I think I answered my own questions by running some optimizations of the A state (CC_STATE_TO_OPT [1,2]) and the NO bond length grew to 1.36 angstrom instead of shrinking. This necessarily threw off the adiabatic excitation energy by more than 1 eV.

#6

I am running 5.4, and (fT) calculations work. Just in case, here are the two alternative ways to run and EOM-IP calculations (and as Anna indicated above, using IP_STATES is the preferred way).

$molecule
0 2
O
H 1 1.2
$end

$rem
method = eom-ccsd(ft)
basis = cc-pvdz
cc_symmetry = false
ip_beta = [1]
eom_nguess_singles = 4
cc_backend = incore
$end

@@@

$molecule
READ
$end

$rem
method = eom-ccsd(ft)
basis = cc-pvdz
cc_symmetry = false
ee_states = [1]
eom_fake_ipea = 1
eom_nguess_singles = 4
cc_backend = incore
$end

In the output:

 EOMIP transition 1/A
 Total energy = -75.10163267 a.u.  Excitation energy = 11.6723 eV.

The EOM-IP-CCSD(fT) energy is: 0.0030922528

 EOMEE transition 1/A
 Total energy = -75.10163268 a.u.  Excitation energy = 11.6723 eV.

 EOM-EE-CCSD(fT) energy = 0.0030922585
#7

Thanks for the info, will something similar to your suggestion work for EA (which was failing with NO)?

#8

EA is a bit trickier than IP to reproduce using via EOM-EE, but here’s an input for NO doublet reference:

$molecule
0 2
N
O 1 1.2
$end

$rem
method = eom-ccsd(ft)
basis = cc-pvdz
scf_algorithm = gdm
cc_symmetry = false
ea_beta = [1]
eom_nguess_singles = 4
cc_backend = incore
$end

@@@

$molecule
0 2
N
O 1 1.2
$end

$rem
method = hf
basis = cc-pvdz
scf_algorithm = gdm
cc_symmetry = false
eom_fake_ipea = 1
$end

@@@

$molecule
-1 1
N
O 1 1.2
$end

$rem
method = eom-ccsd(ft)
basis = cc-pvdz
unrestricted = true
scf_guess = read
scf_max_cycles = 0
cc_symmetry = false
ee_states = [1]
eom_fake_ipea = 1
cc_backend = incore
$end

$occupied
1 2 3 4 5 6 7 8
1 2 3 4 5 6 7 8
$end

Output:

 EOMEA transition 1/A
 Total energy = -129.50980189 a.u.  Excitation energy = 2.0151 eV.

The EOM-EA-CCSD(fT) energy is: -0.0051501946

 EOMEE transition 1/A
 Total energy = -129.50980170 a.u.  Excitation energy = 2.0151 eV.

 EOM-EE-CCSD(fT) energy = -0.0051502254

What kind of failure are you experiencing? Which version of Q-Chem?

#9

Thank you very much for the help with EA as well! For both of these systems I am interested in the first state above the (mostly) degenerate ground states.

The problem my EA calculations looks similar to one mentioned in another thread. I am using 5.3.0 and 5.4.0 and the error occurs with both.

Input file:
Using geom from no-A-opt-eomeaccsd-avdz-a

$molecule
+1 1
O
N 1 bond

bond 1.060669
$end

$rem
CC_SYMMETRY False
JOBTYPE sp
BASIS aug-cc-pvdz
SCF_ALGORITHM gdm
METHOD EOM-CCSD(FT)
EOM_NGUESS_SINGLES 4
EA_STATES [5]
CC_BACKEND incore
$end

Output file:

 Davidson procedure converged

biorthogonolize_left_to_right(): degenerate block detected, ndeg=2

biorthogonolize_left_to_right(): degenerate block detected, ndeg=2

EOMEA-CCSD/MP2 transition 1/A
EOM-EA-CCSD(fT) energy will be calculated without batching

Q-Chem fatal error occurred in module /scratch/qcdevops/jenkins/workspace/build_qchem_linux_distrib/tags/qc540/qchem/ccman2/qchem/ccman2_main.C, line 26:

libtensor::gen_block_tensor<N, BtTraits>::get_block(const index&, bool), /scratch/qcdevops/jenkins/workspace/build_qchem_linux_distrib/tags/qc540/qchem/libtensor/libtensor/gen_block_tensor/impl/gen_bl (189), symmetry_violation
Index does not correspond to a canonical block.

Please submit a crash report at q-chem.com/reporter

#10

Digging into the differences some more, your molecule sections above have 0 2 while mine have +1 1. From the EA operator having one extra creation operator I thought the input file needed to begin with one fewer electron when trying to compute the excited states of the open-shell doublet, which matches up with the EA input files in the manual. When I use the above (+1 1) and eom-ccsd and EA_states = [5] I get adiabatic excitation energies close to 5.31 eV using bond distances that also agree well with the experiments (X: 1.158 A, A: 1.073 A).

Did I miss something regarding the initial charge and spin multiplicity for EA?

#11

The examples above are for illustration only, I wasn’t aiming for the doublet excited states. I can also reproduce the symmetry issue with the latest version of Q-Chem and submitted an internal bug report regarding this error. The EOM-EA in this case can be reproduced with EOM_FAKE_IPEA using the following input:

$molecule
+1 1
O
N 1 bond

bond 1.060669
$end

$rem
CC_SYMMETRY False
JOBTYPE sp
BASIS aug-cc-pvdz
SCF_ALGORITHM gdm
METHOD EOM-CCSD
EOM_NGUESS_SINGLES 4
EA_STATES [5]
CC_BACKEND incore
$end

@@@

$molecule
+1 1
O
N 1 bond

bond 1.060669
$end

$rem
CC_SYMMETRY False
JOBTYPE sp
BASIS aug-cc-pvdz
SCF_ALGORITHM gdm
EOM_FAKE_IPEA 1
METHOD HF
UNRESTRICTED TRUE
$end

@@@

$molecule
0 2
O
N 1 bond

bond 1.060669
$end

$rem
CC_SYMMETRY False
JOBTYPE sp
BASIS aug-cc-pvdz
SCF_GUESS READ
SCF_MAX_CYCLES 0
EOM_FAKE_IPEA 1
METHOD EOM-CCSD(fT)
EE_STATES [5]
CC_BACKEND incore
$end

$occupied
1 2 3 4 5 6 7 14
1 2 3 4 5 6 7
$end

The first section computes EOM-EA energies (for reference) and the second & third sections compute EOM-EA(fT) energies via EOM-EE. You will notice that EOM-EA excitation energies match between the two ways of obtaining EA states.

1 Like
#12

Thank you very much! Those results look great and your explanation makes perfect sense. Thank you as well for submitting the bug report.

Out of curiosity, will the above input template work for EOM-EA-CC(2,3) as well (my EA jobs do not produce any output, but IP and EE seem to work fine)?

#13

I don’t see why not, except that EOM(2,3) will go through the old CC codes in Q-Chem with corresponding limitations. Just use caution and verify results at every step.

1 Like
#14

Thanks! I will give that a try and be sure to look closely at all of the results.

One last question to make sure I use EOM_FAKE_IPEA correctly. If I change basis sets I need to look at the orbital energies and find the position of the -0.000 eV orbital, which I then set in the occupation block. Is this correct?

#15

Correct, the zero-energy orbital must be made occupied when doing EOM-EA via EOM-EE. The number of electrons increases by one and multiplicity must also be adjusted accordingly.

1 Like
#16

EOMEE transition 1/A
Total energy = -129.50980170 a.u. Excitation energy = 2.0151 eV.
EOM-EE-CCSD(fT) energy = -0.0051502254

I saw such an output message in the above posts. I am also doing eom-ccsd(fT) calculation. I was wondering what unit the energy value of EOM-EE-CCSD(fT) is in here and if the energy of eom-ccsd should be added to it. It’s probably explained somewhere, but I missed it. I would be glad if you can help.

#17

The fT correction is in hartrees, and you would add it to the CCSD energy (-129.50980170) to get the total.

#18

Have these issues with EOM-EA-CCSD(fT) and EOM-EA-CC(2,3) been fixed in 6.1 (or earlier)?

#19

We’re pleased to report that we have added a fix for the EOM-IP-CCSD(fT) and EOM-EA-CCSD(fT) calculations; this fix will be in the next Q-Chem bugfix release (6.1.1).