Dear all,
I want to compute the spin-orbit corrected EOM-CCSD-IP energies at the equilibrium geometry of the neutral ground state. To my knowledge, with CALC_SOC keyword, qchem prints the SOCs between all possible pairs of IP states. I have two questions:
Can we make qchem to compute SOCs between all states simultaneously ? It requires the diagonalization of the whole BP Hamiltonian consisting all IP states as diagonal and SOCs as off-diagonal elements.
Can we print the energies of the SO EOM-CCSD-IP states? Something like SO-CASSCF energies printed in Molpro or Molcas output !
My knowledge maybe slightly outdated (I coded EOM SOC a while ago and the current code might be slightly modified), but there is no unique way of doing so. The reason is because EOM-CC formalism is non-Hermitian, and A->B and B-> A transitions have numerically different density matrices. The implemented arithmetic averaging does a good job at predicting the matrix elements, but this can be debated. As usual, the user should carefully monitor this.
I actually advise not to print the resulting eigenenergies in Molpro style due to the following. In any state-interaction scheme one should monitor the rate of convergence of the interested eigenergies or eigenenergy differences with respect to the number of interacting states. Only this monitoring can tell whether the selected manifold of SO-coupled non-relativistic states is sufficient to produce the SO splittings and related properties. Not checking this convergence can lead to incorrect estimates.
You mentioned a great point in the second. Indeed, the number of states/configurations are very important for state interacting calculations whereas only important states are calculated with EOM- CCSD. EOM-CCSD are great for non-relativistic energies. Things get tricky for heavy elements, when multiple states ( >2) contribute to one SO state.
In fact, even for light elements several states may be required to produce a converged estimate of SO splittings. In our SOC EOM-CC paper (link) there is an iron complex requiring at least 3 multiplets.
Developers are aware that printing an SO summary would be helpful, so, I think, it will be available in the future. In the mean time, I can recommend either to write a simple script or to use a trivial python script here (link).