|
| 1 | +#!/usr/bin/env python |
| 2 | + |
| 3 | +import numpy |
| 4 | +import numpy.linalg |
| 5 | +from pyscf import gto, scf, mcscf |
| 6 | + |
| 7 | +mol = gto.M(atom=['H 0 0 %f'%i for i in range(10)], unit='Bohr', |
| 8 | + basis='ccpvtz') |
| 9 | + |
| 10 | +# |
| 11 | +# A regular SCF calculation for this sytem will raise a warning message |
| 12 | +# |
| 13 | +# Warn: Singularity detected in overlap matrix (condition number = 5.47e+09). SCF may be inaccurate and hard to converge. |
| 14 | +# |
| 15 | +# The linear dependency can cause HF, MCSCF etc methods converging to wrong |
| 16 | +# answer. This example shows how to remove linear dependency from overlap |
| 17 | +# matrix and use the linearly independent basis in the HF, MCSCF calculations. |
| 18 | +# |
| 19 | +# There is a shortcut function to remove linear-dependency, eg |
| 20 | +# |
| 21 | +# mf = scf.RHF(mol).apply(scf.addons.remove_linear_dep_) |
| 22 | +# |
| 23 | +# This example demonstrated how the linear dependency is removed in our |
| 24 | +# implementation. |
| 25 | +# |
| 26 | + |
| 27 | +# |
| 28 | +# The smallest eigenvalue of overlap matrix is 10^{-9} |
| 29 | +# |
| 30 | +s = mol.intor('cint1e_ovlp_sph') |
| 31 | +print(numpy.linalg.eigh(s)[0][:8]) |
| 32 | +#[ 1.96568587e-09 8.58358923e-08 7.86870520e-07 1.89728026e-06 |
| 33 | +# 2.14355169e-06 8.96267338e-06 2.46812168e-05 3.26534277e-05] |
| 34 | + |
| 35 | +def eig(h, s): |
| 36 | + d, t = numpy.linalg.eigh(s) |
| 37 | +# Removing the eigenvectors assoicated to the smallest eigenvalue, the new |
| 38 | +# basis defined by x matrix has 139 vectors. |
| 39 | + x = t[:,d>1e-8] / numpy.sqrt(d[d>1e-8]) |
| 40 | + xhx = reduce(numpy.dot, (x.T, h, x)) |
| 41 | + e, c = numpy.linalg.eigh(xhx) |
| 42 | + c = numpy.dot(x, c) |
| 43 | +# Return 139 eigenvalues and 139 eigenvectors. |
| 44 | + return e, c |
| 45 | +# |
| 46 | +# Replacing the default eig function with the above one, the HF solver |
| 47 | +# generate only 139 canonical orbitals |
| 48 | +# |
| 49 | +mf = scf.RHF(mol) |
| 50 | +mf.eig = eig |
| 51 | +mf.verbose = 4 |
| 52 | +mf.kernel() |
| 53 | + |
| 54 | +# |
| 55 | +# The CASSCF solver takes the HF orbital as initial guess. The MCSCF problem |
| 56 | +# size is (0 core, 10 active, 129 external) orbitals. This information can be |
| 57 | +# found in the output. |
| 58 | +# |
| 59 | +mc = mcscf.CASSCF(mf, 10, 10) |
| 60 | +mc.verbose = 4 |
| 61 | +mc.kernel() |
| 62 | + |
| 63 | + |
| 64 | + |
| 65 | +# |
| 66 | +# For symmetry adapted calculation, similar treatments can be applied. |
| 67 | +# |
| 68 | +# Here by assigning symmetry=1, mol.irrep_name, mol.irrep_id and mol.symm_orb |
| 69 | +# (see pyscf/gto/mole.py) are initialized in the mol object. They are the |
| 70 | +# irrep symbols, IDs, and symmetry-adapted-basis. |
| 71 | +# |
| 72 | +mol = gto.M(atom=['H 0 0 %f'%i for i in range(10)], unit='Bohr', |
| 73 | + basis='ccpvtz', symmetry=1) |
| 74 | + |
| 75 | +# |
| 76 | +# The smallest eigenvalue is associated to A1u irrep. Removing the relevant |
| 77 | +# basis will not break the symmetry |
| 78 | +# |
| 79 | +s = mol.intor('cint1e_ovlp_sph') |
| 80 | +for i, c in enumerate(mol.symm_orb): |
| 81 | + s1 = reduce(numpy.dot, (c.T, s, c)) |
| 82 | + print(mol.irrep_name[i], numpy.linalg.eigh(s1)[0]) |
| 83 | +#A1g [ 8.58358928e-08 2.14355169e-06 2.46812168e-05 3.26534277e-05 |
| 84 | +#... |
| 85 | +#E1gx [ 1.67409011e-04 2.38132838e-03 4.51022127e-03 9.89429994e-03 |
| 86 | +#... |
| 87 | +#E1gy [ 1.67409011e-04 2.38132838e-03 4.51022127e-03 9.89429994e-03 |
| 88 | +#... |
| 89 | +#A1u [ 1.96568605e-09 7.86870519e-07 1.89728026e-06 8.96267338e-06 |
| 90 | +#... |
| 91 | + |
| 92 | +# pyscf/scf/hf_symm.py |
| 93 | +def eig(h, s): |
| 94 | + from pyscf import symm |
| 95 | + nirrep = len(mol.symm_orb) |
| 96 | + h = symm.symmetrize_matrix(h, mol.symm_orb) |
| 97 | + s = symm.symmetrize_matrix(s, mol.symm_orb) |
| 98 | + cs = [] |
| 99 | + es = [] |
| 100 | +# |
| 101 | +# Linear dependency are removed by looping over different symmetry irreps. |
| 102 | +# |
| 103 | + for ir in range(nirrep): |
| 104 | + d, t = numpy.linalg.eigh(s[ir]) |
| 105 | + x = t[:,d>1e-8] / numpy.sqrt(d[d>1e-8]) |
| 106 | + xhx = reduce(numpy.dot, (x.T, h[ir], x)) |
| 107 | + e, c = numpy.linalg.eigh(xhx) |
| 108 | + cs.append(reduce(numpy.dot, (mol.symm_orb[ir], x, c))) |
| 109 | + es.append(e) |
| 110 | + e = numpy.hstack(es) |
| 111 | + c = numpy.hstack(cs) |
| 112 | + return e, c |
| 113 | +mf = scf.RHF(mol) |
| 114 | +mf.eig = eig |
| 115 | +mf.verbose = 4 |
| 116 | +mf.kernel() |
| 117 | + |
| 118 | +mc = mcscf.CASSCF(mf, 10, 10) |
| 119 | +mc.verbose = 4 |
| 120 | +mc.kernel() |
0 commit comments