Calculating transition dipole moment between different CI excited states

1

Hello, world! Hi everyone, I am new to psi4 and meet with some problems this week. Thanks for reviewing my topic.

I am trying to implement the Boys diabatization method of excited states with CIS using the psi4numpy interface. To do this, I need to get the transition dipole moments between different CI excited states and the energy of all excited states.

Here are my questions:

  1. is there a way to get CI root energy in the psi4numpy interface?

  2. It seems that I can only get the transition dipole moment between Root 0 and other Roots. I am curious if I can get the transition dipole moment of different CI excited states, i.e., excited state 1 and excited state 2.

Here is my input file:

import psi4
import numpy as np

psi4.set_output_file("output.dat", True)
psi4.set_memory(int(15e9))
numpy_memory = 15

# set molecule
mol = psi4.geometry("""
C    -0.0109277294   -0.0003084931   -0.0000000168
C     1.3195668096   -0.0334284831    0.0000000469
H     1.9161697726    0.8757547346    0.0000016635
H     1.8706155430   -0.9709447058   -0.0000015270
H    -0.5620633765    0.9371525300    0.0000014862
H    -0.6076652628   -0.9094033754   -0.0000016528
C   0.00747483   5.12734788   0.00000047
C   1.30103304   4.83888510   0.00000575
F   2.24935463   5.76457895   0.00000226
F   1.76651200   3.59760567   0.00000333
F  -0.45783697   6.36873346  -0.00000783
F  -0.94085552   4.20168837  -0.00000397
""")

# set option
    psi4.set_options({'basis': '6-31g*',
                      'scf_type': 'pk',
                      'FCI':'false',
                      'opdm':'true',
                      'num_roots':6,
                      'tdm':'true',
                      'dipmom':'true',
                      'ex_level':1,
                      'e_convergence': 1e-8,
                      'd_convergence': 1e-8})


#get excitation energy
ci_e, wfn = psi4.energy('detci', return_wfn=True)
ener = psi4.get_variable('CI ROOT 0 -> ROOT 2 EXCITATION ENERGY')
print(ci_e)
print(ener)


#transition density matrix
opdm_0 = wfn.get_opdm(1,4,"SUM", True)
dm = opdm_0.to_array()
print(dm)


#two electron integral
mints = psi4.core.MintsHelper(wfn.basisset())
I = mints.ao_eri()
D = I.to_array()
print(I.np)

#get transition dipole moment
props = ['DIPOLE', 'TRANSITION_DIPOLE']
psi4.properties('detci',properties=props)
M = ['X','Y','Z']
for i in range(6):
        for j in M:
                dipole=psi4.get_variable('CI ROOT %d -> ROOT %d DIPOLE %s'%(i,i+1,j))
                print('CI ROOT %d -> ROOT %d DIPOLE %s = %s'%(i,i+1,j,dipole))

Some relevant topics and pages: calculating transition and state density matrices, Is there a way to calculate transition dipole moments?, Transition dipole moment using CI, Excited state pes; ci-property

2
  1. Yes. You can see the list of DETCI variables in the documentation, and “CI Root n Total Energy” is in the list.
  2. This feature was added in Psi 1.4. A release candidate is in the works, so you’ll need to use the nightly build to get this. See here for a Python example of the interface.
3

Thank you so much for your reply! It’s really helpful!
I am still curious about question 2. As you said, I downloaded the nightly build and successfully get the transition density matrices. However, I cannot get the transition dipole moment between ci roots that don’t start from root 0. In the output file, it only calculates the transition dipole moment between CI root 0 and other roots. Here is part of my print output.output.dat (90.3 KB)

CI ROOT 0 -> ROOT 1 DIPOLE X = -8.675700964631533e-06
CI ROOT 0 -> ROOT 1 DIPOLE Y = 9.209535034382956e-06
CI ROOT 0 -> ROOT 1 DIPOLE Z = -0.007876256777787869
CI ROOT 1 -> ROOT 2 DIPOLE X = 0.0
CI ROOT 1 -> ROOT 2 DIPOLE Y = 0.0
CI ROOT 1 -> ROOT 2 DIPOLE Z = 0.0
CI ROOT 2 -> ROOT 3 DIPOLE X = 0.0
CI ROOT 2 -> ROOT 3 DIPOLE Y = 0.0
CI ROOT 2 -> ROOT 3 DIPOLE Z = 0.0
CI ROOT 3 -> ROOT 4 DIPOLE X = 0.0
CI ROOT 3 -> ROOT 4 DIPOLE Y = 0.0
CI ROOT 3 -> ROOT 4 DIPOLE Z = 0.0
CI ROOT 4 -> ROOT 5 DIPOLE X = 0.0
CI ROOT 4 -> ROOT 5 DIPOLE Y = 0.0
CI ROOT 4 -> ROOT 5 DIPOLE Z = 0.0
CI ROOT 5 -> ROOT 6 DIPOLE X = 0.0
CI ROOT 5 -> ROOT 6 DIPOLE Y = 0.0
CI ROOT 5 -> ROOT 6 DIPOLE Z = 0.0
CI ROOT 0 -> ROOT 2 DIPOLE X = 3.8604125133899414
CI ROOT 0 -> ROOT 2 DIPOLE Y = 0.0008397200658493508
CI ROOT 0 -> ROOT 2 DIPOLE Z = -1.9931410534163673e-06s
4

Ah, you’re right. Looking over the code, MultipolePropCalc::compute_dipole in oeprop.cc doesn’t access the TDMs computed by detci but just uses the last density stored on it, regardless of whether that’s a transition density or not. That’s annoying.

While there is probably a workaround involving creating an OEProp object and manually setting the density, it certainly won’t be user-friendly.

Unfortunately, I have other bugs to worry about and can’t commit to fixing this one. I do recommend you open an issue about this.

1 Like
5

Thanks for your information! I am not familiar with C++. I guess I can calculate the transition dipole moment from ao_dipole, transition density matrix, and nuclear dipole moment. I have opened an issue in GitHub (#2090).