1#!/usr/bin/env python
2#
3# Author: Qiming Sun <osirpt.sun@gmail.com>
4#
5
6import numpy
7from pyscf import gto
8from pyscf import scf
9
10'''
11Mixing decoration, for density fitting, scalar relativistic effects, and
12second order (Newton-Raphson) SCF.
13
14Density fitting and scalar relativistic effects can be applied together,
15regardless to the order you apply the decoration.
16
17NOTE the second order SCF (New in version 1.1) decorating operation are not
18commutable with scf.density_fit operation
19        [scf.density_fit, scf.sfx2c1e    ] == 0
20        [scf.newton     , scf.sfx2c1e    ] != 0
21        [scf.newton     , scf.density_fit] != 0
22* scf.density_fit(scf.newton(scf.RHF(mol))) is the SOSCF for regular 2e
23  integrals, but with density fitting integrals for the Hessian.  It's an
24  approximate SOSCF optimization method;
25* scf.newton(scf.density_fit(scf.RHF(mol))) is the exact second order
26  optimization for the given scf object which is a density-fitted-scf method.
27  The SOSCF is not an approximate scheme.
28* scf.density_fit(scf.newton(scf.density_fit(scf.RHF(mol))), auxbasis='ahlrichs')
29  is an approximate SOSCF scheme for the given density-fitted-scf method.
30  Here we use small density fitting basis (ahlrichs cfit basis) to approximate
31  the Hessian for the large-basis-density-fitted-scf scheme.
32'''
33
34mol = gto.Mole()
35mol.build(
36    verbose = 0,
37    atom = '''8  0  0.     0
38              1  0  -0.757 0.587
39              1  0  0.757  0.587''',
40    basis = 'ccpvdz',
41)
42
43#
44# 1. spin-free X2C-HF with density fitting approximation on 2E integrals
45#
46mf = scf.density_fit(scf.sfx2c1e(scf.RHF(mol)))
47mf = scf.RHF(mol).x2c().density_fit()  # Stream style
48energy = mf.kernel()
49print('E = %.12f, ref = -76.075408156180' % energy)
50
51#
52# 2. spin-free X2C correction for density-fitting HF.  Since X2C correction is
53# commutable with density fitting operation, it is fully equivalent to case 1.
54#
55mf = scf.sfx2c1e(scf.density_fit(scf.RHF(mol)))
56mf = scf.RHF(mol).density_fit().x2c()  # Stream style
57energy = mf.kernel()
58print('E = %.12f, ref = -76.075408156180' % energy)
59
60#
61# 3. The order to apply X2C or newton method matters.  If relativistic effects
62# need to be included in the calculation, you should first call x2c then
63# newton method.
64#
65e1 = scf.RHF(mol).kernel()
66e2 = scf.RHF(mol).x2c().kernel()
67print('E(NR)         = %.12f  E(X2C)        = %.12f' % (e1, e2))
68e1 = scf.RHF(mol).newton().x2c().kernel()
69e2 = scf.RHF(mol).x2c().newton().kernel()
70print('E(Newton,X2C) = %.12f  E(X2C,Newton) = %.12f' % (e1, e2))
71
72#
73# 4. Newton method for non-relativistic HF
74#
75mf = scf.newton(scf.RHF(mol))
76mf = scf.RHF(mol).newton()  # Stream style
77energy = mf.kernel()
78print('E = %.12f, ref = -76.026765673120' % energy)
79
80#
81# 5. Newton method for non-relativistic HF with density fitting for orbital
82# hessian of newton solver.  Note the answer is equal to case 3, but the
83# solver "mf" is different.
84#
85mf = scf.density_fit(scf.newton(scf.RHF(mol)))
86mf = scf.RHF(mol).newton().density_fit()
87energy = mf.kernel()
88print('E = %.12f, ref = -76.026765673120' % energy)
89
90#
91# 6. Newton method to solve the density-fitting approximated HF object.  There
92# is no approximation for newton method (orbital hessian).  Note the density
93# fitting is applied on HF object only.  It does not affect the Newton solver.
94#
95mf = scf.newton(scf.density_fit(scf.RHF(mol)))
96mf = scf.RHF(mol).density_fit().newton()
97energy = mf.kernel()
98print('E = %.12f, ref = -76.026744737357' % energy)
99
100#
101# 7. Newton method for density-fitting HF, and the hessian of Newton solver is
102# also approximated with density fitting.  Note the anwser is equivalent to
103# case 6, but the solver "mf" is different.  Here the fitting basis for HF and
104# Newton solver are different.  HF is approximated with the default density
105# fitting basis (Weigend cfit basis).  Newton solver is approximated with
106# Ahlrichs cfit basis.
107#
108mf = scf.density_fit(scf.newton(scf.density_fit(scf.RHF(mol))), 'ahlrichs')
109mf = scf.RHF(mol).density_fit().newton().density_fit(auxbasis='ahlrichs')
110energy = mf.kernel()
111print('E = %.12f, ref = -76.026744737357' % energy)
112
113