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