1#!/usr/bin/env python
2#
3# Author: Oliver J. Backhouse <olbackhouse@gmail.com>
4#         George H. Booth <george.booth@kcl.ac.uk>
5#
6
7'''
8An example of restricted AGF2 with density fitting, obtaining the 1RDM and dipole moment
9
10Default AGF2 corresponds to the AGF2(1,0) method outlined in the papers:
11  - O. J. Backhouse, M. Nusspickel and G. H. Booth, J. Chem. Theory Comput., 16, 1090 (2020).
12  - O. J. Backhouse and G. H. Booth, J. Chem. Theory Comput., 16, 6294 (2020).
13'''
14
15from pyscf import gto, scf, agf2
16import numpy as np
17from functools import reduce
18
19mol = gto.M(atom='O 0 0 0; H 0 0 1; H 0 1 0', basis='cc-pvdz')
20
21mf = scf.RHF(mol).density_fit(auxbasis='cc-pv5z-ri')
22mf.conv_tol = 1e-12
23mf.run()
24
25# Run an AGF2 calculation
26gf2 = agf2.AGF2(mf)
27gf2.conv_tol = 1e-7
28gf2.run(verbose=4)
29
30# Print the first 3 ionization potentials
31gf2.ipagf2(nroots=3)
32
33# Print the first 3 electron affinities
34gf2.eaagf2(nroots=3)
35
36# Get the MO-basis density matrix and calculate dipole moments:
37dm = gf2.make_rdm1()
38
39dipole = [0.0, 0.0, 0.0]
40# Transform dipole moment integrals into MO basis
41mol.set_common_origin([0,0,0])
42r_ints_ao = mol.intor('cint1e_r_sph', comp=3)
43r_ints_mo = np.empty_like(r_ints_ao)
44for i in range(3):
45    r_ints_mo[i] = reduce(np.dot,(mf.mo_coeff.T, r_ints_ao[i], mf.mo_coeff))
46    dipole[i] = -np.trace(np.dot(dm, r_ints_mo[i]))
47    # Add nuclear component
48    for j in range(mol.natm):
49        dipole[i] += mol.atom_charge(j) * mol.atom_coord(j)[i]
50
51print('Dipole moment from AGF2: {} {} {}'.format(dipole[0], dipole[1], dipole[2]))
52