iOpenShell » Q-Chem specific questions » Problem in calculation of Dyson orbitals for EOM-IP spin-orbit states

Problem in calculation of Dyson orbitals for EOM-IP spin-orbit states

Page: 1

Author Post
Registered: Apr 2017
Posts: 2
Dear Q-Chem users and developers,
I have started using Q-Chem recently. I am trying to compute Dyson orbitals between the ground singlet (neutral molecule ) and cationic doublet states using the EOM-IP-CCSD + SOC formalism of Q-Chem 5.0. However the job crashes after the spin-orbit calculations between the EOM-IP states are finished, before going to the Dyson orbital calculation step. Again the calculation finishes normally if spin-orbit coupling (SOC) is deactivated. Following are the $rem part of my input and part of the output containing error message -


basis general
ecp general
cc_symmetry false
ccman2 true
method eom-ccsd
purecart 111
ip_states [6] !Target EOM-IP states
calc_soc true
cc_state_to_opt [1,1]
cc_do_dyson true ! calculate Dyson orbitals
cc_trans_prop true
state_analysis true
ianlty 200
make_cube_files true



1/A <--> 6/A
Transition dipole moment (a.u.):
A->B: 0.243062 (X -0.150342, Y -0.190988, Z 0.000000)
B->A: 0.261530 (X -0.161709, Y -0.205544, Z 0.000000)
Oscillator strength (a.u.): 0.009063
Norm of one-particle transition density matrix:
A->B: 0.824515; B->A: 0.981015
One-electron spin-orbit coupling (cm-1):
A->B: (Re 0.000000, Im 0.000000)
B->A: (Re 0.000000, Im 0.000000)
Avg: (Re 0.000000, Im 0.000000)
Mean-field spin-orbit coupling (cm-1):
A->B: (Re 0.000000, Im 0.000000)
B->A: (Re 0.000000, Im 0.000000)
Avg: (Re 0.000000, Im 0.000000)
Atomic CT numbers
omega = 1.6345

Exciton analysis of the transition density matrix
<r_h> [Ang]: [ 0.000000, 0.000000, -0.170445]
<r_e> [Ang]: [ -0.000000, -0.000000, -0.019511]
|<r_e - r_h>| [Ang]: 0.150935
Hole size [Ang]: 1.628268
Cartesian components [Ang]: [ 0.540238, 0.539057, 1.438339]
Electron size [Ang]: 1.727301
Cartesian components [Ang]: [ 0.759768, 0.854831, 1.294443]
RMS electron-hole separation [Ang]: 2.293645
Cartesian components [Ang]: [ 0.927491, 1.004201, 1.841778]
Covariance(r_h, r_e) [Ang^2]: 0.198399
Correlation coefficient: 0.070542
Center-of-mass size [Ang]: 1.227968
Cartesian components [Ang]: [ 0.468500, 0.508482, 1.014820]

EOMIP-CCSD calculation: CPU 541.74 s wall 541.63 s


contents of the context



p0_58468: p4_error: interrupt SIGx: 6

In this regard I have few questions to ask-

1. Is it possible to calculate Dyson orbitals for the EOM-IP spin-orbit states ? If so, how ? I have been able to calculate dyson orbital for EOM-IP states. But for my system SOC plays an important role.

2. To calculate properties (like SOC) between ip-states, the keyword 'cc_state_to_opt' specifies the initial eom-state . In case of calculating properties between two different pairs of eom states how would one use the 'cc_state_to_op' keyword ?

3.How to analyze the spin-orbit part of the output eg., how to extract the spin-orbit splitting between two interacting states ? I looked into the manual but could not get enough details about it.

Any kind of help/suggestions in this regard will be highly appreciated. Sorry for the length of the question.


Registered: Oct 2017
Posts: 17
1. SOCs for EOM-CC in Q-Chem are computed PERTURBATIVELY. This means that the states are computed without SO interaction first, and then SOC is computed as an interstate property without changing these electronic states. If I understand your request correctly, you want to have Dyson orbitals for the states of relativistic nature, which is not implemented (it would be possible to do this with variational treatment of SO). If your calculation would not crash, you would see the same Dyson orbitals as without SO. For details see

2. cc_state_to_opt specifies from which excited/IP/EA state the transition density matrix will be computed. The syntax is the following (taken from the manual):
Specifies which state to optimize (or from which state compute EOM-EOM inter-state
[i,j] optimize the jth state of the ith irrep.

If you want to calculate couplings between, for example, the second ionized state (without symmetry) and all other states, put CC_STATE_TO_OPT=[1,2]

3. This not easy, but not impossible. It requires some knowledge of quantum mechanics. At first, you need to form a SO matrix. The actual SOCs are not invariant with respect to rotations, but SO operator is a scalar operator, and its form is covariant (it does not change). Therefore, you need to consider how your states change when you do rotations of the coordinate system. It depends on multiplicity of the states, their spin projection, and point group irrep (for non-abelian point groups of symmetry the states will rotate in the multidimensional irrep, changing irreps of abelian point subgroups). In the simplest case you will need only Wigner matrices, coming from spin. Having them, rotate the system (do not forget to specify symmetry = false and sym_ignore = true), and substitute the couplings for the rotated states into a system of equations with Wigner matrices. If you solve the system of equations, you will have the whole SOC matrix for the initial molecular orientation. If you add energies on the diagonal, and diagonalize the matrix, you will get the energies of the states with SO splitting included.

4. You do not specify the system that you study, making hard to make suggestions. Just be careful with ECP - not all of them are good for SO calculations. Q-Chem by default does not compute full SOC, only one-electron and two-electron from two-electron density matrix (SO mean field). You may need full SO operator (which is more expensive) - in this case request two-electron properties.
« Last edit by pavel on Sat Oct 14, 2017 6:04 am. »
Registered: Apr 2017
Posts: 2
Thank you very much for the help and information. The system I am working on is ICN. I am using aug-cc-pVTZ basis for C, N from Q-Chem library and the following basis for Iodine -

I 0
I-ECP 3 36
2 2.00659990 -2.60678697
2 5.41200018 -15.30154037
2 20.88789940 -73.32185364
1 63.82109833 -23.43061638
2 1.53009999 -144.29801941
2 1.85660005 252.49822998
2 2.58430004 -200.15640259
2 3.75419998 85.98131561
1 1.13870001 32.24239349
0 12.47010040 6.45168209
2 0.91109997 32.43734360
2 1.06579995 -101.57602692
2 1.38989997 169.25859070
2 1.79579997 -105.27653503
1 2.54789996 12.12060356
0 1.17250001 6.45733500
2 0.86530000 3.41730189
2 0.89960003 -3.97894001
2 4.58720016 -37.31042862
2 5.52820015 77.44996643
1 22.20800018 18.01509857
0 19.14010048 6.86909294

I 0
S 1 1.00
0.6803000 1.00000000
S 1 1.00
0.5208000 1.00000000
S 1 1.00
0.1143000 1.00000000
P 1 1.00
1.1250000 1.00000000
P 1 1.00
0.3508000 1.00000000
P 1 1.00
0.1166000 1.00000000
D 1 1.00
3.9502000 1.00000000
D 1 1.00
2.5742000 1.00000000
D 1 1.00
1.3071000 1.00000000
D 1 1.00
0.5116000 1.00000000
D 1 1.00
0.2660000 1.00000000
SP 1 1.00
0.08000 1.00000000 1.00000000
SP 1 1.00
0.03680 1.00000000 1.00000000
SP 1 1.00
0.007360 1.00000000 1.00000000
SP 1 1.00
0.001470 1.00000000 1.00000000
SP 1 1.00
0.000294 1.00000000 1.00000000
SP 1 1.00
0.0000588 1.00000000 1.00000000
SP 1 1.00
0.0000118 1.00000000 1.00000000
SP 1 1.00
0.00000235 1.00000000 1.00000000

The basis gives excellent values for ionization energies of the IP states. However the spin-orbit coupling values between the IP states is very less. For example, the SOC between the first two IP states are ~ 15 cm-1 (for one electron SOC) and ~ 5 cm-1 (SOMF). However the literature value is around 4000 cm-1.

Does the discrepancy in SOC comes from the basis for Iodine being used ? I know little about the above basis since I got it from a Q-Chem sample input. I do not really know the source of the basis or any references regarding the basis. I also took the contribution from two electron SOC. But that gave values similar like SOMF.

Thanks for your time.

« Last edit by soumitra on Fri Oct 20, 2017 8:39 am. »
Registered: Oct 2017
Posts: 17
The magnitude of several thousands cm-1 should not surprise, it is typical for iodine. What are the numbers 5 cm-1 and 15 cm-1? Just matrix elements or rotationary invariant constants? If they are matrix elements, than as I wrote before, they depend on the orientation of the molecule.

SOC is sensitive to what is happening with the core orbitals, especially for heavy elements (see the paper about SOC). What is this ECP? Does it provide a correct nodal structure? SO calculations require special integrals, and I do not know if they were implemented for ECP case. Can you run your calculation without ECP (very crude, but it will elucidate if there is a problem with integrals)?

Page: 1

iOpenShell » Q-Chem specific questions » Problem in calculation of Dyson orbitals for EOM-IP spin-orbit states