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