1#!/usr/bin/env python
2
3'''
4Transition density matrix and transition dipole for CASSCF/CASCI wavefunction
5'''
6
7from functools import reduce
8import numpy
9from pyscf import gto, scf, mcscf, fci
10
11mol = gto.M(atom='''
12C   1.20809735,    0.69749533,   0.00000000
13C   0.00000000,    1.39499067,   0.00000000
14C   -1.20809735,   0.69749533,   0.00000000
15C   -1.20809735,  -0.69749533,   0.00000000
16C   -0.00000000,  -1.39499067,   0.00000000
17C   1.20809735,   -0.69749533,   0.00000000
18H   2.16038781,    1.24730049,   0.00000000
19H   0.00000000,    2.49460097,   0.00000000
20H   -2.16038781,   1.24730049,   0.00000000
21H   -2.16038781,  -1.24730049,   0.00000000
22H   0.00000000,   -2.49460097,   0.00000000
23H   2.16038781,   -1.24730049,   0.00000000''',
24verbose = 4,
25basis = 'ccpvdz',
26symmetry = True,
27)
28mf = scf.RHF(mol)
29mf.kernel()
30#mf.analyze()
31
32#
33# 1. State-average CASSCF to get optimal orbitals
34#
35mc = mcscf.CASSCF(mf, 6, 6)
36solver_ag = fci.direct_spin0_symm.FCI(mol)
37solver_b2u = fci.direct_spin0_symm.FCI(mol)
38solver_b2u.wfnsym = 'B2u'
39mc.fcisolver = mcscf.state_average_mix(mc, [solver_ag,solver_b2u], [.5,.5])
40cas_list = [17,20,21,22,23,30]  # 2pz orbitals
41mo = mcscf.sort_mo(mc, mf.mo_coeff, cas_list)
42mc.kernel(mo)
43#mc.analyze()
44mc_mo = mc.mo_coeff
45
46#
47# 2. Ground state wavefunction.  This step can be passed you approximate it
48# with the state-averaged CASSCF wavefunction
49#
50mc = mcscf.CASCI(mf, 6, 6)
51mc.fcisolver.wfnsym = 'Ag'
52mc.kernel(mc_mo)
53ground_state = mc.ci
54
55#
56# 3. Exited states.  In this example, B2u are bright states.
57#
58# Here, mc.ci[0] is the first excited state.
59#
60mc = mcscf.CASCI(mf, 6, 6)
61mc.fcisolver.wfnsym = 'B2u'
62mc.fcisolver.nroots = 8
63mc.kernel(mc_mo)
64
65#
66# 4. transition density matrix and transition dipole
67#
68# Be careful with the gauge origin of the dipole integrals
69#
70charges = mol.atom_charges()
71coords = mol.atom_coords()
72nuc_charge_center = numpy.einsum('z,zx->x', charges, coords) / charges.sum()
73mol.set_common_orig_(nuc_charge_center)
74dip_ints = mol.intor('cint1e_r_sph', comp=3)
75
76def makedip(ci_id):
77    # transform density matrix in MO representation
78    t_dm1 = mc.fcisolver.trans_rdm1(ground_state, mc.ci[ci_id], mc.ncas, mc.nelecas)
79    # transform density matrix to AO representation
80    orbcas = mc_mo[:,mc.ncore:mc.ncore+mc.ncas]
81    t_dm1_ao = reduce(numpy.dot, (orbcas, t_dm1, orbcas.T))
82    # transition dipoles
83    return numpy.einsum('xij,ji->x', dip_ints, t_dm1_ao)
84
85# 1st and 6th excited states are B2u of D6h point group, dark states
86# 3rd and 4th excited states are triplet states, dipole == 0
87for i in range(8):
88    print('Transition dipole between |0> and |%d>'%(i+1), makedip(i))
89
90