Excited state normal mode frequency doesn't match for NH3

Question

I am trying to optimize the geometry of excited state NH3 and get the corresponding normal modes for excited states. I am trying to get the frequencies as mentioned in this reference. They claim to use CC2/TZVPP for calculations. I am trying to re-create their results with my script but there is a big mismatch between the values they report and I am getting.

Their reported values:

I am able to get their ground state frequencies but not the excited state. My GAMESS script for excited state calculation:

!   File created for NH3 Excited states
 $BASIS GBASIS=TZV $END
 $CONTRL SCFTYP=RHF RUNTYP=OPTIMIZE CCTYP=CCD $END
 $SCF DIRSCF=.TRUE. $END
 $STATPT  OPTTOL=0.0001 NSTEP=20 HSSEND=.T. $END
 $FORCE METHOD=FULLNUM VIBSIZ=0.001 VIBANL=.TRUE. PURIFY=.TRUE. $END
 $eominp nstate(1)=1,0,0,0,0,0,0,0 iroot(1)=2 ccprpe=.true.
         minit=2 noact=3 nuact=5 $end
 $DATA
Title
C1
N     7.0    -0.00000    -0.00000     0.05500
H     1.0    -0.01280    -0.96010    -0.25500
H     1.0     0.83790     0.46900    -0.25500
H     1.0    -0.82510     0.49110    -0.25500
 $END

The resultant frequencies that I get after running this is: {364.17, 1680.03, 1686.99, 3658.94, 3866.18, 3867.96} Could you help in fixing the issue?

Edit:

Based on one of the answers, I updated the input file. Resulting frequencies with EOM-CCSD and TZVPP basis functions are, {237.63, 875.84, 982.98, 3508.41, 3524.99, 3533.73} Still way off than what is being reported. But it looks like moving in the right direction.

!   File created for NH3 Excited states
 $BASIS GBASIS=TZVPP $END
 $CONTRL SCFTYP=RHF RUNTYP=OPTIMIZE CCTYP=EOM-CCSD $END
 $SCF DIRSCF=.TRUE. $END
 $STATPT  OPTTOL=0.0001 NSTEP=20 HSSEND=.T. $END
 $FORCE METHOD=FULLNUM VIBSIZ=0.001 VIBANL=.TRUE. PURIFY=.TRUE. $END
 $eominp nstate(1)=1,0,0,0,0,0,0,0 iroot(1)=2 ccprpe=.true.
         minit=2 noact=3 nuact=5 $end
 $DATA
Title
C1
N     7.0    -0.00000    -0.00000     0.05500
H     1.0    -0.01280    -0.96010    -0.25500
H     1.0     0.83790     0.46900    -0.25500
H     1.0    -0.82510     0.49110    -0.25500
 $END

Answer 1

This is not a full answer, because it does not solve the problem. But I hope to shed some light on why you are not getting the correct frequencies.

Short answer: I suspect there is a bug in GAMESS

Long answer:

First of all I would advise you to go through the manual of GAMESS. There are lots of problems in the input file, and most of them can be easily fixed by looking at the manual.

1. You are using EOM-CCSD, so you have to use CCTYP=EOM-CCSD.
2. This method does not have analytic gradients, so you must add NUMGRD=.t. in $CONTRL. Your input file as of now, calculates the analytic gradient of the RHF part only, and produces the wrong result. (EOM-CCSD calculation has three parts—ground state RHF, ground state CCSD, excited state EOM-CCSD)
3. You cannot use PURIFY=.t. without internal coordinates! Since you are using Cartesian coordinates, the correct keyword is PROJCT=.t.
4. The basis set and method problem has been already mentioned in the answer by @Tyberius. The basis set you are looking for is def2-TZVPP, which is not available in GAMESS. (GBASIS=KTZVPP is not the same basis) . You will have to input it from outside. (As a blatant self-plug, I have written a blog post about using outside basis sets here. It's too long to write here anyways)
5. The nstate(1)=1,0,0,0,0,0,0,0 keyword request for only one symmetric excited state. You should use iroot(1)=1,1. The first number in iroot(1)=indicates the symmetry (1 for completely symmetric i.e. A). The second number indicates the state (1 for first excited state). (Your current input file requests a state that does not exist)
6. The ccprpe=.true. keyword calculates the EOM-CCSD density matrix for properties. Unless you need properties (like dipole moment etc.) I would suggest you turn this off. Remember, the density matrix is calculated at every stage of geometry optimization and numerical differentiaion for gradient, so it's clearly wasting cpu resources.

Finally, GAMESS cannot run numerical optimizations and numerical hessians in the same run (i.e. you cannot use HSSEND=.t. with NUMGRD=.t.). Therefore, you have to optimize first, with NUMGRD=.t., then take that optimized structure and then calculate its hessian with RUNTYP=HESSIAN and METHOD=FULLNUM.

Even after fixing all of this, when I try to calculate the hessian of the excited $\ce{NH3}$, I get 3 imaginary frequencies. Which is odd, because I just optimized it! I believe the numeric hessian code of GAMESS is not using the EOM-CCSD first excited state energy values, but is instead using the ground state CCSD or RHF energy values for calculating the hessian. I am not completely sure if this is a bug, or whether some keyword needs to be changed. The good news is that at least the gradient calculation run correctly (I am getting a planar optimized structure for the excited state).

Answer 2

I can see two major problems:

1. CCD != CC2 and I'm not certain if CC2 is available in GAMESS. CC2 is an approximation of CCSD, which is available in GAMESS. You may be able to use this since $\ce{NH3}$ is a fairly small molecule, though I don't know how time consuming this will be relative to your initial CCD calculation.

2. TZV != TZVPP. The def2-TZVPP basis set has two sets of polarization functions added to it relative to def2-TZV (TZV denotes an entirely different basis set in GAMESS). The def2 basis sets are actually not available by default in GAMESS and must be pulled from an external source (the KTZV and KTZVPP in GAMESS are the older def basis sets). Using a different basis set will significantly affect the results.

I haven't used GAMESS, so there could be other issues with your input.