Unexpected ROKS difference density

Hello everyone!

I have been trying to compare how ROKS and TDDFT describe the lowest singlet excited state of a Cu(I) complex. It is well known that the S1 state of such a complex is of MLCT nature. I wanted to compare the difference densities obtained by TDDFT and ROKS methods, and while the TDDFT difference density looks reasonable, the ROKS difference density is somewhat bizarre - both in terms of the shape (the density around copper doesn’t even resemble d orbitals), and the sign (blue colour corresponds to increase in electron density, and red - decrease. Since S1 is MLCT state, the electron density should shift from Cu to the diimine ligand - ROKS difference density suggests the opposite direction).

What could be the reason for such results?

S1 TDDFT difference density:

S1 ROKS difference density:

It would be better to show your Q-Chem input files of two jobs. Especially for S1 ROKS, this is tricky and thus it is important to know whether your calculation is reasonable.

Oh, right, sorry. Here is the input for ROKS job (I’m omitting the $basis section at the end):

$molecule
1 1
C -0.0861094 0.2614760 4.4284010
C 0.5096588 0.3757065 3.1816150
C 1.0496294 1.5837268 2.7359403
C 1.0015825 2.6871167 3.5861737
C 0.4124758 2.5873670 4.8404721
C -0.1333138 1.3783079 5.2546878
O 0.6338685 -0.6873111 2.3114184
C -0.3459637 -1.6608830 2.2800812
C -1.5413836 -1.4062650 1.6026722
C -2.4754908 -2.4388502 1.5126745
C -2.2186179 -3.6821212 2.0770633
C -1.0196214 -3.9120354 2.7408295
C -0.0762404 -2.8967669 2.8442780
P -1.7757750 0.2056633 0.7680290
Cu -0.0099175 1.0374774 -0.4300330
N -0.4956207 2.3714615 -1.9491262
N -0.8717911 3.5984554 -2.0622932
N -0.9896102 3.8322542 -3.3630775
C -0.6887778 2.7479908 -4.1100858
C -0.3648591 1.7919690 -3.1788684
H -1.2783050 4.7547966 -3.6728069
C 0.0612041 0.3934239 -3.2684041
N 0.3547329 -0.1762412 -2.1124075
C 0.7603755 -1.4765272 -2.0871350
C 0.8556498 -2.2408204 -3.2793213
C 0.5337349 -1.6016071 -4.4985454
C 0.1432029 -0.2921021 -4.4993900
C 1.2637142 -3.5940535 -3.1972013
C 1.5636959 -4.1527215 -1.9842185
C 1.4784912 -3.3819513 -0.8021984
C 1.0905607 -2.0705839 -0.8484187
P 1.6202624 1.6680347 1.0064964
H 1.0327142 -1.4617953 0.0513896
H 1.7249879 -3.8376432 0.1562885
H 1.8749591 -5.1949944 -1.9250534
H 1.3340135 -4.1802782 -4.1136924
H 0.6002212 -2.1621707 -5.4315054
H -0.1095659 0.2179378 -5.4271557
H -0.7256958 2.7520097 -5.1932227
H -3.0269040 -0.0460659 0.1565302
H -2.2500413 1.0089958 1.8356931
H 0.8744906 -3.0470020 3.3546538
H -0.8161283 -4.8867822 3.1831365
H -2.9599008 -4.4759748 1.9944196
H -3.4158473 -2.2722399 0.9870773
H 1.4191255 3.6397364 3.2594399
H 0.3785236 3.4572959 5.4948016
H -0.5972867 1.2965012 6.2371540
H -0.5036939 -0.6913384 4.7499609
H 2.8470419 0.9599768 1.0487696
H 2.1497972 2.9800872 0.9947822
$end

$rem
jobtype sp
method camb3lyp
dft_d d3_bj
basis gen
purecart 11
PRINT_GENERAL_BASIS false
gen_scfman false
MAX_SCF_CYCLES 500
SCF_GUESS read
SCF_CONVERGENCE 9
scf_final_print true
MEM_STATIC 2000
MEM_TOTAL 80000
ROKS true
MOLDEN_FORMAT true
IQMOL_FCHK true
plots true
make_cube_files true
$end

$plots
grid_points 200 200 200
total_density 0
$end

TDDFT input:

$molecule
read
$end

$rem
jobtype sp
method camb3lyp
dft_d d3_bj
basis gen
purecart 11
PRINT_GENERAL_BASIS false
gen_scfman false
MAX_SCF_CYCLES 500
SCF_GUESS read
SCF_CONVERGENCE 9
scf_final_print true
MEM_STATIC 2000
MEM_TOTAL 80000
CIS_N_ROOTS 10
CIS_TRIPLETS false
CIS_SINGLETS true
CIS_RELAXED_DENSITY true
RPA false
MOLDEN_FORMAT true
IQMOL_FCHK true
state_analysis true
make_cube_files true
plots true
$end

$plots
grid_points 200 200 200
$end

This is second part if the S1 tddft calculation, in the first part I only generate the GS orbitals in my basis set. The basis set is the same in all my calculations.

Have you tried looking at the overlap with the ground state? You accomplish this by using the TRANS_MOM_SAVE = 1 $rem variable for a previous ground state calculation and then TRANS_MOM_READ = 2 for the ROKS calculation. Large overlaps are bad news. I believe the GS calculation needs to be unrestricted for TRANS_MOM_SAVE to work…

I re-ran the calculation with this option, and the overlap was zero.

Could you post the output for that ROKS calculation?

Sorry, I mixed up two geometries. I get zero for the perpendicular geometry of the system, and here I posted data for the 95 degree dihedral. I’m sorry for the confusion.

So in the proper case the overlap was 0.020555

Output of the ROKS calculation:

Standard Electronic Orientation quadrupole field applied
 Nucleus-field energy     =     0.0000000649 hartrees
 Guess MOs from SCF MO coefficient file
 Reading MOs from coefficient file
 Reading MOs from coefficient file
 Long-range K will be added via erf
 Coulomb attenuation parameter = 0.33 bohr**(-1)
 A restricted hybrid HF-DFT SCF calculation will be
 performed using Pulay DIIS extrapolation
 Exchange:     0.1900 Hartree-Fock + 0.3500 B88 + 0.4600 muB88 + LR-HF
 Correlation:  0.1900 VWN5 + 0.8100 LYP
 Using SG-1 standard quadrature grid
 Dispersion:   Grimme D3
 SCF converges when DIIS error is below 1.0E-09
 Using Q-Chem read-in guess as SCF_GUESS READ specified.

ROKS singly occupied orbitals: (125,126)
 using 24 threads for integral computing
 -------------------------------------------------------
 OpenMP Integral computing Module                
 Release: version 1.0, May 2013, Q-Chem Inc. Pittsburgh 
 -------------------------------------------------------
 Building long-range K (w = 3.300e-01)
 Building long-range K (w = 3.300e-01)
Timing for roks_project_singlet(): 3.95 sec.
Timing for roks_build(): 349.88 sec.
 ---------------------------------------
  Cycle       Energy         DIIS Error
 ---------------------------------------
    1   -3505.0703927334      1.25E-02
  Det[Ca.S.C0] = -0.760754
  Det[Cb.S.C0] = -0.005371
  <S1|S0>  = 0.005778
 Building long-range K (w = 3.300e-01)
 Building long-range K (w = 3.300e-01)
Timing for roks_project_singlet(): 3.90 sec.
Timing for roks_build(): 376.18 sec.
    2   -3502.1006273839      9.13E-03
  Det[Ca.S.C0] = -0.879849
  Det[Cb.S.C0] = -0.005003
  <S1|S0>  = 0.006226
 Building long-range K (w = 3.300e-01)
 Building long-range K (w = 3.300e-01)
Timing for roks_project_singlet(): 3.92 sec.
Timing for roks_build(): 364.24 sec.
    3   -3504.0822708645      4.59E-03
  Det[Ca.S.C0] = -0.936623
  Det[Cb.S.C0] = -0.004648
  <S1|S0>  = 0.006157
 Building long-range K (w = 3.300e-01)
 Building long-range K (w = 3.300e-01)
Timing for roks_project_singlet(): 3.90 sec.
Timing for roks_build(): 358.29 sec.
    4   -3504.7295149752      2.56E-03
  Det[Ca.S.C0] = 0.967282
  Det[Cb.S.C0] = 0.006658
  <S1|S0>  = 0.009108
 Building long-range K (w = 3.300e-01)
 Building long-range K (w = 3.300e-01)
Timing for roks_project_singlet(): 3.92 sec.
Timing for roks_build(): 350.97 sec.
    5   -3504.9883217290      1.32E-03
  Det[Ca.S.C0] = 0.974097
  Det[Cb.S.C0] = 0.005653
  <S1|S0>  = 0.007787
 Building long-range K (w = 3.300e-01)
 Building long-range K (w = 3.300e-01)
Timing for roks_project_singlet(): 4.00 sec.
Timing for roks_build(): 350.68 sec.
    6   -3505.0346733375      9.74E-04
  Det[Ca.S.C0] = -0.978704
  Det[Cb.S.C0] = -0.006548
  <S1|S0>  = 0.009064
 Building long-range K (w = 3.300e-01)
 Building long-range K (w = 3.300e-01)
Timing for roks_project_singlet(): 3.97 sec.
Timing for roks_build(): 348.45 sec.
    7   -3505.0633345492      6.93E-04
  Det[Ca.S.C0] = 0.981559
  Det[Cb.S.C0] = 0.007090
  <S1|S0>  = 0.009841
 Building long-range K (w = 3.300e-01)
 Building long-range K (w = 3.300e-01)
Timing for roks_project_singlet(): 3.90 sec.
Timing for roks_build(): 348.70 sec.
    8   -3505.0793772928      4.80E-04
  Det[Ca.S.C0] = 0.983893
  Det[Cb.S.C0] = 0.007813
  <S1|S0>  = 0.010872
 Building long-range K (w = 3.300e-01)
 Building long-range K (w = 3.300e-01)
Timing for roks_project_singlet(): 3.90 sec.
Timing for roks_build(): 349.69 sec.
    9   -3505.0916994333      2.14E-04
  Det[Ca.S.C0] = 0.984519
  Det[Cb.S.C0] = 0.008355
  <S1|S0>  = 0.011633
 Building long-range K (w = 3.300e-01)
 Building long-range K (w = 3.300e-01)
Timing for roks_project_singlet(): 3.93 sec.
Timing for roks_build(): 348.59 sec.
   10   -3505.0945412366      6.19E-05
  Det[Ca.S.C0] = 0.984593
  Det[Cb.S.C0] = 0.008613
  <S1|S0>  = 0.011993
 Building long-range K (w = 3.300e-01)
 Building long-range K (w = 3.300e-01)
Timing for roks_project_singlet(): 3.93 sec.
Timing for roks_build(): 348.22 sec.
   11   -3505.0947841246      1.54E-05
  Det[Ca.S.C0] = -0.984601
  Det[Cb.S.C0] = -0.008776
  <S1|S0>  = 0.012220
 Building long-range K (w = 3.300e-01)
 Building long-range K (w = 3.300e-01)
Timing for roks_project_singlet(): 3.93 sec.
Timing for roks_build(): 347.81 sec.
   12   -3505.0948024863      2.22E-06
  Det[Ca.S.C0] = -0.984597
  Det[Cb.S.C0] = -0.008908
  <S1|S0>  = 0.012404
 Building long-range K (w = 3.300e-01)
 Building long-range K (w = 3.300e-01)
Timing for roks_project_singlet(): 3.93 sec.
Timing for roks_build(): 349.90 sec.
   13   -3505.0948046481      1.49E-06
  Det[Ca.S.C0] = -0.984592
  Det[Cb.S.C0] = -0.009072
  <S1|S0>  = 0.012632
 Building long-range K (w = 3.300e-01)
 Building long-range K (w = 3.300e-01)
Timing for roks_project_singlet(): 3.93 sec.
Timing for roks_build(): 349.12 sec.
   14   -3505.0948067508      1.32E-06
  Det[Ca.S.C0] = 0.984585
  Det[Cb.S.C0] = 0.009397
  <S1|S0>  = 0.013085
 Building long-range K (w = 3.300e-01)
 Building long-range K (w = 3.300e-01)
Timing for roks_project_singlet(): 3.95 sec.
Timing for roks_build(): 349.83 sec.
   15   -3505.0948090291      1.17E-06
  Det[Ca.S.C0] = -0.984573
  Det[Cb.S.C0] = -0.010028
  <S1|S0>  = 0.013963
 Building long-range K (w = 3.300e-01)
 Building long-range K (w = 3.300e-01)
Timing for roks_project_singlet(): 3.94 sec.
Timing for roks_build(): 349.03 sec.
   16   -3505.0948100492      1.04E-06
  Det[Ca.S.C0] = -0.984556
  Det[Cb.S.C0] = -0.010862
  <S1|S0>  = 0.015124
 Building long-range K (w = 3.300e-01)
 Building long-range K (w = 3.300e-01)
Timing for roks_project_singlet(): 3.93 sec.
Timing for roks_build(): 348.79 sec.
   17   -3505.0948052955      8.77E-07
  Det[Ca.S.C0] = -0.984532
  Det[Cb.S.C0] = -0.011941
  <S1|S0>  = 0.016625
 Building long-range K (w = 3.300e-01)
 Building long-range K (w = 3.300e-01)
Timing for roks_project_singlet(): 3.88 sec.
Timing for roks_build(): 349.13 sec.
   18   -3505.0947908689      6.96E-07
  Det[Ca.S.C0] = -0.984505
  Det[Cb.S.C0] = -0.012932
  <S1|S0>  = 0.018005
 Building long-range K (w = 3.300e-01)
 Building long-range K (w = 3.300e-01)
Timing for roks_project_singlet(): 3.87 sec.
Timing for roks_build(): 347.66 sec.
   19   -3505.0947703687      5.13E-07
  Det[Ca.S.C0] = -0.984481
  Det[Cb.S.C0] = -0.013673
  <S1|S0>  = 0.019036
 Building long-range K (w = 3.300e-01)
 Building long-range K (w = 3.300e-01)
Timing for roks_project_singlet(): 3.91 sec.
Timing for roks_build(): 348.02 sec.
   20   -3505.0947509411      3.70E-07
  Det[Ca.S.C0] = -0.984466
  Det[Cb.S.C0] = -0.014069
  <S1|S0>  = 0.019587
 Building long-range K (w = 3.300e-01)
 Building long-range K (w = 3.300e-01)
Timing for roks_project_singlet(): 3.90 sec.
Timing for roks_build(): 347.81 sec.
   21   -3505.0947394175      2.76E-07
  Det[Ca.S.C0] = -0.984456
  Det[Cb.S.C0] = -0.014326
  <S1|S0>  = 0.019945
 Building long-range K (w = 3.300e-01)
 Building long-range K (w = 3.300e-01)
Timing for roks_project_singlet(): 3.91 sec.
Timing for roks_build(): 348.50 sec.
   22   -3505.0947315107      2.04E-07
  Det[Ca.S.C0] = 0.984452
  Det[Cb.S.C0] = 0.014459
  <S1|S0>  = 0.020131
 Building long-range K (w = 3.300e-01)
 Building long-range K (w = 3.300e-01)
Timing for roks_project_singlet(): 3.92 sec.
Timing for roks_build(): 349.20 sec.
   23   -3505.0947270909      1.66E-07
  Det[Ca.S.C0] = 0.984449
  Det[Cb.S.C0] = 0.014548
  <S1|S0>  = 0.020254
 Building long-range K (w = 3.300e-01)
 Building long-range K (w = 3.300e-01)
Timing for roks_project_singlet(): 3.90 sec.
Timing for roks_build(): 347.91 sec.
   24   -3505.0947243592      1.79E-07
  Det[Ca.S.C0] = 0.984446
  Det[Cb.S.C0] = 0.014630
  <S1|S0>  = 0.020369
 Building long-range K (w = 3.300e-01)
 Building long-range K (w = 3.300e-01)
Timing for roks_project_singlet(): 3.93 sec.
Timing for roks_build(): 348.17 sec.
   25   -3505.0947216356      9.30E-08
  Det[Ca.S.C0] = 0.984445
  Det[Cb.S.C0] = 0.014687
  <S1|S0>  = 0.020448
 Building long-range K (w = 3.300e-01)
 Building long-range K (w = 3.300e-01)
Timing for roks_project_singlet(): 3.89 sec.
Timing for roks_build(): 348.73 sec.
   26   -3505.0947196915      8.32E-08
  Det[Ca.S.C0] = 0.984443
  Det[Cb.S.C0] = 0.014734
  <S1|S0>  = 0.020512
 Building long-range K (w = 3.300e-01)
 Building long-range K (w = 3.300e-01)
Timing for roks_project_singlet(): 3.91 sec.
Timing for roks_build(): 349.04 sec.
   27   -3505.0947180776      5.98E-08
  Det[Ca.S.C0] = 0.984443
  Det[Cb.S.C0] = 0.014750
  <S1|S0>  = 0.020535
 Building long-range K (w = 3.300e-01)
 Building long-range K (w = 3.300e-01)
Timing for roks_project_singlet(): 3.88 sec.
Timing for roks_build(): 348.78 sec.
   28   -3505.0947174973      3.46E-08
  Det[Ca.S.C0] = 0.984442
  Det[Cb.S.C0] = 0.014760
  <S1|S0>  = 0.020549
 Building long-range K (w = 3.300e-01)
 Building long-range K (w = 3.300e-01)
Timing for roks_project_singlet(): 3.91 sec.
Timing for roks_build(): 348.29 sec.
   29   -3505.0947171109      2.37E-08
  Det[Ca.S.C0] = 0.984442
  Det[Cb.S.C0] = 0.014767
  <S1|S0>  = 0.020558
 Building long-range K (w = 3.300e-01)
 Building long-range K (w = 3.300e-01)
Timing for roks_project_singlet(): 3.90 sec.
Timing for roks_build(): 348.00 sec.
   30   -3505.0947168584      1.60E-08
  Det[Ca.S.C0] = 0.984442
  Det[Cb.S.C0] = 0.014768
  <S1|S0>  = 0.020560
 Building long-range K (w = 3.300e-01)
 Building long-range K (w = 3.300e-01)
Timing for roks_project_singlet(): 3.90 sec.
Timing for roks_build(): 348.52 sec.
   31   -3505.0947168073      7.73E-09
  Det[Ca.S.C0] = 0.984442
  Det[Cb.S.C0] = 0.014766
  <S1|S0>  = 0.020557
 Building long-range K (w = 3.300e-01)
 Building long-range K (w = 3.300e-01)
Timing for roks_project_singlet(): 3.69 sec.
Timing for roks_build(): 349.31 sec.
   32   -3505.0947168674      3.58E-09
  Det[Ca.S.C0] = 0.984442
  Det[Cb.S.C0] = 0.014766
  <S1|S0>  = 0.020557
 Building long-range K (w = 3.300e-01)
 Building long-range K (w = 3.300e-01)
Timing for roks_project_singlet(): 3.90 sec.
Timing for roks_build(): 348.58 sec.
   33   -3505.0947168716      2.22E-09
  Det[Ca.S.C0] = 0.984442
  Det[Cb.S.C0] = 0.014765
  <S1|S0>  = 0.020556
 Building long-range K (w = 3.300e-01)
 Building long-range K (w = 3.300e-01)
Timing for roks_project_singlet(): 3.90 sec.
Timing for roks_build(): 349.67 sec.
   34   -3505.0947168959      1.35E-09
  Det[Ca.S.C0] = 0.984442
  Det[Cb.S.C0] = 0.014765
  <S1|S0>  = 0.020556
 Building long-range K (w = 3.300e-01)
 Building long-range K (w = 3.300e-01)
Timing for roks_project_singlet(): 3.90 sec.
Timing for roks_build(): 347.50 sec.
   35   -3505.0947168974      1.09E-09
  Det[Ca.S.C0] = 0.984442
  Det[Cb.S.C0] = 0.014765
  <S1|S0>  = 0.020556
 Building long-range K (w = 3.300e-01)
 Building long-range K (w = 3.300e-01)
Timing for roks_project_singlet(): 3.84 sec.
Timing for roks_build(): 347.10 sec.
   36   -3505.0947169094      9.12E-10 Convergence criterion met
 ---------------------------------------
 One-Electron    Energy = -12934.6127500560
 Total Coulomb   Energy =  5629.0257291864
 Alpha Exchange  Energy =   -35.3734631443
 Beta  Exchange  Energy =   -36.1129889774
 DFT   Exchange  Energy =  -199.3215623933
 DFT Correlation Energy =   -10.8076683049
 Nuclear Repu.   Energy =  4082.1079867803
 Nuclear Attr.   Energy = -16428.0189075641
 Kinetic         Energy =  3493.4061575081
  Det[Ca.S.C0] = 0.984442
  Det[Cb.S.C0] = 0.014764
  <S1|S0>  = 0.020555
 SCF time:  CPU 12699.87 s  wall 626.31 s

That overlap does not look bad at all. Is there any particular reason you’re using GEN_SCFMAN = FALSE? You could try with GEN_SCFMAN and use different SCF solvers (DIIS + MOM, SGM, I even use GDM some times for these…) Also, are you sure the guess you’re providing is correct? In other words, are you sure that the HOMO and LUMO from the guess look like the NTOs of the state you’re trying to target?

I ran a test on a small basis with GEN_SCFMAN = TRUE and using GDM as the SCF solver, plus a small trick I implemented for ROKS (will be out in 6.2). I cant show you a difference density because I dont know how to plot those but the two open-shell orbitals of the ROKS calculations seem consistent with the difference density you show for TDDFT. Isovalues of 0.100.

OS 1 (Cu d-like):

OS 2 (ligand pi*):

I’d suggest looking at your guess orbitals to see they are correct and maybe try different solvers.

Hi, could I just follow up on this?

I tried to find the trans_mom_read option in the manual but did not see this. Would that allow me to compute ROKS transition moments?

Cheers
Felix

So long as the orbitals for the initial state have been saved in a previous calc with the trans_mom_save, yes.

1 Like

I finally got back to the ROKS calculations.
Thank you for your advice! Changing the SCF solver to GDM and GEN_SCFMAN to true worked, the difference density looks similarly to TD-DFT one.
I was using gen_scfman = false because in the first job I was generating ground state orbitals in a basis set that’s not built in in Q-Chem (DZP).
However, I have a couple more questions.
I am doing a PES scan of S1 state along the rotation of N-Cu-N and P-Cu-P planes. Using the two SCF solvers I got the following picture:

Also, the overlap between S1 and S0 state is visibly different (and higher) when using GDM:
obraz
(absolute values are depicted)

Are such differences, both in energy (up to 0,1 eV) and overlap to be expected when using two different solvers?

Glad it worked. The unfortunate answer is yes. Unlike ROHF for high-spin open-shell systems, the ROKS energy is not invariant to rotations among the open shells. As a consequence the orbital space may have several local minima that look wonky (weird difference densities, large overlaps with the ground state, bad energies…) and different solvers will sometimes land you in different ones. In general, I like the ROKS solutions that have the smallest overlap with the ground state and that “look” correct (MOs / orbital differences make sense). In my experience they have often gone hand-in-hand but yours is the first time I see when a smaller overlap (DIIS solution) produces wonky orbitals. Could it be that the wonky DIIS solution is a legit solution that just doesnt look like the one you’re interested on? (i.e. a nearby excited state)

As I mentioned, I implemented a trick that supressess the open-shell mixing in ROKS and it does wonders with GDM and SGM. So, if you provide it a guess with sensible open-shell orbitals (ground state orbitals) it’ll try to prevent the open-shells from getting scrambled up and producing weird results. It is an unpublished idea since I didnt have the time to write it up so use at your own discretion! It’ll be out in 6.2 so you can try it then. I’d also be happy to run a test one for you to see what happens but I’d need the exact input file.