Commit 11ef281d authored by hebrewsnabla's avatar hebrewsnabla
Browse files

allow rhf_internal in check_stab

parent f19d501b
......@@ -5,7 +5,7 @@ try:
import gaussian
except:
print('fch2py not found. Interface with fch is disabled. Install MOKIT if you need that.')
from automr import stability, dump_mat
from automr import stability, dump_mat, autocas
import time
import copy
......@@ -130,18 +130,25 @@ def mix(xyz, bas, charge=0, conv='loose', cycle=5, skipstb=False, xc=None):
#mf_mix.kernel(dm)
return mf_mix
def check_stab(mf_mix, newton=False):
def check_stab(mf_mix, newton=False, res=False):
if res:
stab = stability.rhf_internal
is_uhf, mf_mix = autocas.check_uhf(mf_mix)
if is_uhf:
raise ValueError('UHF/UKS wavefunction detected. RHF/RKS is required for res=True')
else:
stab = stability.uhf_internal
mf_mix.verbose = 9
mo, stable = stability.uhf_internal(mf_mix)
mo, stable = stab(mf_mix)
cyc = 0
while(not stable and cyc < 10):
mf_mix.verbose = 4
dm_new = scf.uhf.make_rdm1(mo, mf_mix.mo_occ)
dm_new = mf_mix.make_rdm1(mo, mf_mix.mo_occ)
if newton:
mf_mix=mf_mix.newton()
mf_mix.kernel(dm0=dm_new)
mf_mix.verbose = 9
mo, stable = stability.uhf_internal(mf_mix)
mo, stable = stab(mf_mix)
cyc += 1
if not stable:
raise RuntimeError('Stablility Opt failed after %d attempts.' % cyc)
......
......@@ -8,6 +8,39 @@ from pyscf.lib import logger
from pyscf.soscf import newton_ah
from pyscf.scf.stability import _rotate_mo
def rhf_internal(mf, with_symmetry=True, verbose=None):
log = logger.new_logger(mf, verbose)
log.note('**** checking RHF/RKS internal stability ...')
g, hop, hdiag = newton_ah.gen_g_hop_rhf(mf, mf.mo_coeff, mf.mo_occ,
with_symmetry=with_symmetry)
hdiag *= 2
stable = True
def precond(dx, e, x0):
hdiagd = hdiag - e
hdiagd[abs(hdiagd)<1e-8] = 1e-8
return dx/hdiagd
# The results of hop(x) corresponds to a displacement that reduces
# gradients g. It is the vir-occ block of the matrix vector product
# (Hessian*x). The occ-vir block equals to x2.T.conj(). The overall
# Hessian for internal reotation is x2 + x2.T.conj(). This is
# the reason we apply (.real * 2) below
def hessian_x(x):
return hop(x).real * 2
x0 = numpy.zeros_like(g)
x0[g!=0] = 1. / hdiag[g!=0]
if not with_symmetry: # allow to break point group symmetry
x0[numpy.argmin(hdiag)] = 1
e, v = lib.davidson(hessian_x, x0, precond, tol=1e-4, verbose=log)
if e < -1e-5:
log.note('RHF/RKS wavefunction has an internal instability')
mo = _rotate_mo(mf.mo_coeff, mf.mo_occ, v)
stable = False
else:
log.note('RHF/RKS wavefunction is stable in the internal stability analysis')
mo = mf.mo_coeff
return mo, stable
def uhf_internal(mf, with_symmetry=True, verbose=None):
log = logger.new_logger(mf, verbose)
log.note('**** checking UHF/UKS internal stability ...')
......@@ -37,4 +70,4 @@ def uhf_internal(mf, with_symmetry=True, verbose=None):
log.note('UHF/UKS wavefunction is stable in the internal stability analysis')
mo = mf.mo_coeff
stable = True
return mo, stable
\ No newline at end of file
return mo, stable
converged SCF energy = -2088.28076496514
**** checking RHF/RKS internal stability ...
nelec by numeric integration = 48.000004161259334
CPU time for vxc 1.66 sec, wall time 0.05 sec
tol 0.0001 toloose 0.01
max_cycle 50 max_space 12 max_memory 2000 incore True
davidson 0 1 |r|= 4.11 e= [2.061283] max|de|= 2.06 lindep= 0.997
davidson 1 2 |r|= 0.991 e= [1.749281] max|de|= -0.312 lindep= 0.993
davidson 2 3 |r|= 2.89 e= [1.636928] max|de|= -0.112 lindep= 0.977
davidson 3 4 |r|= 1.2 e= [1.460154] max|de|= -0.177 lindep= 0.944
davidson 4 5 |r|= 0.786 e= [1.229078] max|de|= -0.231 lindep= 0.971
davidson 5 6 |r|= 0.707 e= [1.021425] max|de|= -0.208 lindep= 0.983
davidson 6 7 |r|= 0.294 e= [0.796616] max|de|= -0.225 lindep= 0.973
davidson 7 8 |r|= 0.201 e= [0.77611] max|de|= -0.0205 lindep= 0.882
davidson 8 9 |r|= 0.213 e= [0.773811] max|de|= -0.0023 lindep= 0.836
davidson 9 10 |r|= 0.816 e= [0.385896] max|de|= -0.388 lindep= 0.95
davidson 10 11 |r|= 0.485 e= [0.276048] max|de|= -0.11 lindep= 0.962
davidson 11 12 |r|= 0.174 e= [0.238401] max|de|= -0.0376 lindep= 0.844
davidson 12 1 |r|= 0.174 e= [0.238401] max|de|= 1.39e-16 lindep= 1
davidson 13 2 |r|= 0.14 e= [0.236465] max|de|= -0.00194 lindep= 1
davidson 14 3 |r|= 0.0768 e= [0.234487] max|de|= -0.00198 lindep= 0.863
davidson 15 4 |r|= 0.0476 e= [0.233472] max|de|= -0.00101 lindep= 0.861
davidson 16 5 |r|= 0.0174 e= [0.233168] max|de|= -0.000304 lindep= 0.937
root 0 converged |r|= 0.0077 e= 0.2331109827275325 max|de|= -5.68e-05
converged 17 6 |r|= 0.0077 e= [0.233111] max|de|= -5.68e-05
RHF/RKS wavefunction is stable in the internal stability analysis
from pyscf import gto, dft
from automr import guess
#mf=guess.from_fch_simp("v2.fchk", xc='pbe0')
#mf2.verbose=9
#mf2.stability()
mol = gto.Mole(atom='''Cr 0.0 0.0 0.0; Cr 0.0 0.0 1.4''', basis='def2-tzvp').build()
mf = dft.RKS(mol)
mf.xc = 'pbe0'
mf.kernel()
mf2 = guess.check_stab(mf, newton=True, res=True)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment