diff --git a/pyscf/dft/xcfun.py b/pyscf/dft/xcfun.py index 8fd498b338..3d937f3f92 100644 --- a/pyscf/dft/xcfun.py +++ b/pyscf/dft/xcfun.py @@ -252,7 +252,7 @@ 'B97XC' , 'B97_1XC' , 'B97_2XC' , 'M05XC' , 'TPSSH' , 'HFLYP'} RSH_XC = {'CAMB3LYP'} -MAX_DERIV_ORDER = ctypes.c_int.in_dll(_itrf, 'XCFUN_max_deriv_order').value +MAX_DERIV_ORDER = 3 VV10_XC = { 'B97M_V' : (6.0, 0.01), @@ -338,8 +338,12 @@ def rsh_coeff(xc_code): beta = hyb - alpha return omega, alpha, beta +@lru_cache(100) def max_deriv_order(xc_code): - return MAX_DERIV_ORDER + try: + return ctypes.c_int.in_dll(_itrf, 'XCFUN_max_deriv_order').value + except ValueError: # make it compatible with old xcfun_itrf.so + return MAX_DERIV_ORDER def test_deriv_order(xc_code, deriv, raise_error=False): support = deriv <= max_deriv_order(xc_code) diff --git a/pyscf/hessian/rhf.py b/pyscf/hessian/rhf.py index 9b1d69fc00..e56d3b46a6 100644 --- a/pyscf/hessian/rhf.py +++ b/pyscf/hessian/rhf.py @@ -296,7 +296,8 @@ def _get_jk(mol, intor, comp, aosym, script_dms, return vs def solve_mo1(mf, mo_energy, mo_coeff, mo_occ, h1ao_or_chkfile, - fx=None, atmlst=None, max_memory=4000, verbose=None, max_cycle=50): + fx=None, atmlst=None, max_memory=4000, verbose=None, + max_cycle=50, level_shift=0): '''Solve the first order equation Kwargs: @@ -343,7 +344,8 @@ def _ao2mo(mat): h1vo = numpy.vstack(h1vo) s1vo = numpy.vstack(s1vo) - mo1, e1 = cphf.solve(fx, mo_energy, mo_occ, h1vo, s1vo, max_cycle=max_cycle) + mo1, e1 = cphf.solve(fx, mo_energy, mo_occ, h1vo, s1vo, + max_cycle=max_cycle, level_shift=level_shift) mo1 = numpy.einsum('pq,xqi->xpi', mo_coeff, mo1).reshape(-1,3,nao,nocc) e1 = e1.reshape(-1,3,nocc,nocc) @@ -470,8 +472,15 @@ def h_op(x): class HessianBase(lib.StreamObject): '''Non-relativistic restricted Hartree-Fock hessian''' + # Max. number of iterations for Krylov solver + max_cycle = 50 + # Shift virtual orbitals to slightly improve the convergence speed of Krylov solver + # A small level_shift ~ 0.1 is often helpful to decrease 2 - 3 iterations + # while the error of cphf solver may be increased by one magnitude. + level_shift = 0 + _keys = { - 'mol', 'base', 'chkfile', 'atmlst', 'de', 'max_cycle' + 'mol', 'base', 'chkfile', 'atmlst', 'de', 'max_cycle', 'level_shift' } def __init__(self, scf_method): @@ -481,7 +490,6 @@ def __init__(self, scf_method): self.base = scf_method self.chkfile = scf_method.chkfile self.max_memory = self.mol.max_memory - self.max_cycle = 50 self.atmlst = range(self.mol.natm) self.de = numpy.zeros((0,0,3,3)) # (A,B,dR_A,dR_B) @@ -566,7 +574,8 @@ def get_hcore(iatm, jatm): def solve_mo1(self, mo_energy, mo_coeff, mo_occ, h1ao_or_chkfile, fx=None, atmlst=None, max_memory=4000, verbose=None): return solve_mo1(self.base, mo_energy, mo_coeff, mo_occ, h1ao_or_chkfile, - fx, atmlst, max_memory, verbose, max_cycle=self.max_cycle) + fx, atmlst, max_memory, verbose, + max_cycle=self.max_cycle, level_shift=self.level_shift) def hess_nuc(self, mol=None, atmlst=None): if mol is None: mol = self.mol diff --git a/pyscf/hessian/test/test_rhf.py b/pyscf/hessian/test/test_rhf.py index e1ae1f7087..7f3bfdb2f9 100644 --- a/pyscf/hessian/test/test_rhf.py +++ b/pyscf/hessian/test/test_rhf.py @@ -35,6 +35,18 @@ def tearDownModule(): del mol class KnownValues(unittest.TestCase): + def test_rhf_hess(self): + mf = scf.RHF(mol) + e0 = mf.kernel() + hess = hessian.RHF(mf).kernel() + self.assertAlmostEqual(lib.fp(hess), -0.7816352153153946, 6) + + hobj = hessian.RHF(mf) + hobj.max_cycle = 10 + hobj.level_shift = .1 + hess = hobj.kernel() + self.assertAlmostEqual(lib.fp(hess), -0.7816352153153946, 6) + def test_finite_diff_x2c_rhf_hess(self): mf = scf.RHF(mol).x2c() mf.conv_tol = 1e-14 diff --git a/pyscf/hessian/test/test_uhf.py b/pyscf/hessian/test/test_uhf.py index eb2622d10b..06d32b38ad 100644 --- a/pyscf/hessian/test/test_uhf.py +++ b/pyscf/hessian/test/test_uhf.py @@ -36,6 +36,15 @@ def tearDownModule(): del mol class KnownValues(unittest.TestCase): + def test_uhf_hess(self): + mf = scf.UHF(mol) + mf.conv_tol = 1e-14 + e0 = mf.kernel() + hobj = mf.Hessian() + hobj.level_shift = .05 + hess = hobj.kernel() + self.assertAlmostEqual(lib.fp(hess), -0.20243405976628576, 5) + def test_finite_diff_rhf_hess(self): mf = scf.UHF(mol) mf.conv_tol = 1e-14 diff --git a/pyscf/hessian/uhf.py b/pyscf/hessian/uhf.py index 1be5ccc587..1b30e264ee 100644 --- a/pyscf/hessian/uhf.py +++ b/pyscf/hessian/uhf.py @@ -257,7 +257,8 @@ def make_h1(hessobj, mo_coeff, mo_occ, chkfile=None, atmlst=None, verbose=None): return chkfile def solve_mo1(mf, mo_energy, mo_coeff, mo_occ, h1ao_or_chkfile, - fx=None, atmlst=None, max_memory=4000, verbose=None): + fx=None, atmlst=None, max_memory=4000, verbose=None, + max_cycle=50, level_shift=0): mol = mf.mol if atmlst is None: atmlst = range(mol.natm) @@ -306,7 +307,8 @@ def _ao2mo(mat, mo_coeff, mocc): h1vo = (numpy.vstack(h1voa), numpy.vstack(h1vob)) s1vo = (numpy.vstack(s1voa), numpy.vstack(s1vob)) - mo1, e1 = ucphf.solve(fx, mo_energy, mo_occ, h1vo, s1vo) + mo1, e1 = ucphf.solve(fx, mo_energy, mo_occ, h1vo, s1vo, + max_cycle=max_cycle, level_shift=level_shift) mo1a = numpy.einsum('pq,xqi->xpi', mo_coeff[0], mo1[0]).reshape(-1,3,nao,nocca) mo1b = numpy.einsum('pq,xqi->xpi', mo_coeff[1], mo1[1]).reshape(-1,3,nao,noccb) e1a = e1[0].reshape(-1,3,nocca,nocca) @@ -449,7 +451,8 @@ class Hessian(rhf_hess.HessianBase): def solve_mo1(self, mo_energy, mo_coeff, mo_occ, h1ao_or_chkfile, fx=None, atmlst=None, max_memory=4000, verbose=None): return solve_mo1(self.base, mo_energy, mo_coeff, mo_occ, h1ao_or_chkfile, - fx, atmlst, max_memory, verbose) + fx, atmlst, max_memory, verbose, + max_cycle=self.max_cycle, level_shift=self.level_shift) def to_gpu(self): raise NotImplementedError diff --git a/pyscf/lib/linalg_helper.py b/pyscf/lib/linalg_helper.py index fe2892364c..1fcc1a9265 100644 --- a/pyscf/lib/linalg_helper.py +++ b/pyscf/lib/linalg_helper.py @@ -1290,9 +1290,9 @@ def krylov(aop, b, x0=None, tol=1e-10, max_cycle=30, dot=numpy.dot, >>> from pyscf import lib >>> a = numpy.random.random((10,10)) * 1e-2 >>> b = numpy.random.random(10) - >>> aop = lambda x: numpy.dot(a,x) + >>> aop = lambda x: a.dot(x.T).T >>> x = lib.krylov(aop, b) - >>> numpy.allclose(numpy.dot(a,x)+x, b) + >>> numpy.allclose(aop(x)+x, b) True ''' if isinstance(aop, numpy.ndarray) and aop.ndim == 2: diff --git a/pyscf/lib/test/test_linalg_helper.py b/pyscf/lib/test/test_linalg_helper.py index 6cc7e463c4..6c8c57f56f 100644 --- a/pyscf/lib/test/test_linalg_helper.py +++ b/pyscf/lib/test/test_linalg_helper.py @@ -116,6 +116,24 @@ def precond(x, *args): x1 = linalg_helper.krylov(aop, b/a_diag, x1, max_cycle=30) self.assertAlmostEqual(abs(xref - x1).max(), 0, 6) + def test_krylov_with_level_shift(self): + numpy.random.seed(10) + n = 100 + a = numpy.random.rand(n,n) * .1 + a = a.dot(a.T) + a_diag = numpy.random.rand(n) + b = numpy.random.rand(n) + ref = numpy.linalg.solve(numpy.diag(a_diag) + a, b) + + #((diag+shift) + (a-shift)) x = b + shift = .1 + a_diag += shift + a -= numpy.eye(n)*shift + + aop = lambda x: (a.dot(x.T).T/a_diag) + c = linalg_helper.krylov(aop, b/a_diag, max_cycle=18) + self.assertAlmostEqual(abs(ref - c).max(), 0, 9) + def test_dgeev(self): numpy.random.seed(12) n = 100 diff --git a/pyscf/scf/cphf.py b/pyscf/scf/cphf.py index 19e2d2d10b..73cbdc2010 100644 --- a/pyscf/scf/cphf.py +++ b/pyscf/scf/cphf.py @@ -27,7 +27,8 @@ def solve(fvind, mo_energy, mo_occ, h1, s1=None, - max_cycle=20, tol=1e-9, hermi=False, verbose=logger.WARN): + max_cycle=50, tol=1e-9, hermi=False, verbose=logger.WARN, + level_shift=0): ''' Args: fvind : function @@ -36,29 +37,43 @@ def solve(fvind, mo_energy, mo_occ, h1, s1=None, Kwargs: hermi : boolean Whether the matrix defined by fvind is Hermitian or not. + level_shift : float + Add to diagonal terms to slightly improve the convergence speed of + Krylov solver ''' if s1 is None: return solve_nos1(fvind, mo_energy, mo_occ, h1, - max_cycle, tol, hermi, verbose) + max_cycle, tol, hermi, verbose, level_shift) else: return solve_withs1(fvind, mo_energy, mo_occ, h1, s1, - max_cycle, tol, hermi, verbose) + max_cycle, tol, hermi, verbose, level_shift) kernel = solve # h1 shape is (:,nvir,nocc) def solve_nos1(fvind, mo_energy, mo_occ, h1, - max_cycle=20, tol=1e-9, hermi=False, verbose=logger.WARN): - '''For field independent basis. First order overlap matrix is zero''' + max_cycle=50, tol=1e-9, hermi=False, verbose=logger.WARN, + level_shift=0): + '''For field independent basis. First order overlap matrix is zero + + Kwargs: + level_shift : float + Add to diagonal terms to slightly improve the convergence speed of + Krylov solver + ''' + assert not hermi log = logger.new_logger(verbose=verbose) t0 = (logger.process_clock(), logger.perf_counter()) e_a = mo_energy[mo_occ==0] e_i = mo_energy[mo_occ>0] - e_ai = 1 / lib.direct_sum('a-i->ai', e_a, e_i) + e_ai = 1 / (e_a[:,None] + level_shift - e_i) mo1base = h1 * -e_ai def vind_vo(mo1): - v = fvind(mo1.reshape(h1.shape)).reshape(h1.shape) + mo1 = mo1.reshape(h1.shape) + v = fvind(mo1).reshape(h1.shape) + if level_shift != 0: + v -= mo1 * level_shift v *= e_ai return v.ravel() mo1 = lib.krylov(vind_vo, mo1base.ravel(), @@ -68,20 +83,23 @@ def vind_vo(mo1): # h1 shape is (:,nocc+nvir,nocc) def solve_withs1(fvind, mo_energy, mo_occ, h1, s1, - max_cycle=20, tol=1e-9, hermi=False, verbose=logger.WARN): + max_cycle=50, tol=1e-9, hermi=False, verbose=logger.WARN, + level_shift=0): '''For field dependent basis. First order overlap matrix is non-zero. The first order orbitals are set to C^1_{ij} = -1/2 S1 e1 = h1 - s1*e0 + (e0_j-e0_i)*c1 + vhf[c1] Kwargs: - hermi : boolean - Whether the matrix defined by fvind is Hermitian or not. + level_shift : float + Add to diagonal terms to slightly improve the convergence speed of + Krylov solver Returns: First order orbital coefficients (in MO basis) and first order orbital energy matrix ''' + assert not hermi log = logger.new_logger(verbose=verbose) t0 = (logger.process_clock(), logger.perf_counter()) @@ -89,34 +107,38 @@ def solve_withs1(fvind, mo_energy, mo_occ, h1, s1, viridx = mo_occ == 0 e_a = mo_energy[viridx] e_i = mo_energy[occidx] - e_ai = 1 / lib.direct_sum('a-i->ai', e_a, e_i) + e_ai = 1 / (e_a[:,None] + level_shift - e_i) nvir, nocc = e_ai.shape nmo = nocc + nvir s1 = s1.reshape(-1,nmo,nocc) hs = mo1base = h1.reshape(-1,nmo,nocc) - s1*e_i - mo_e1 = hs[:,occidx,:].copy() + mo1base = hs.copy() mo1base[:,viridx] *= -e_ai mo1base[:,occidx] = -s1[:,occidx] * .5 def vind_vo(mo1): - v = fvind(mo1.reshape(h1.shape)).reshape(-1,nmo,nocc) + mo1 = mo1.reshape(mo1base.shape) + v = fvind(mo1).reshape(mo1base.shape) + if level_shift != 0: + v -= mo1 * level_shift v[:,viridx,:] *= e_ai v[:,occidx,:] = 0 return v.ravel() mo1 = lib.krylov(vind_vo, mo1base.ravel(), tol=tol, max_cycle=max_cycle, hermi=hermi, verbose=log) mo1 = mo1.reshape(mo1base.shape) + mo1[:,occidx] = mo1base[:,occidx] log.timer('krylov solver in CPHF', *t0) - v1mo = fvind(mo1.reshape(h1.shape)).reshape(-1,nmo,nocc) - mo1[:,viridx] = mo1base[:,viridx] - v1mo[:,viridx]*e_ai + hs += fvind(mo1).reshape(mo1base.shape) + mo1[:,viridx] = hs[:,viridx] / (e_i - e_a[:,None]) # mo_e1 has the same symmetry as the first order Fock matrix (hermitian or # anti-hermitian). mo_e1 = v1mo - s1*lib.direct_sum('i+j->ij',e_i,e_i) - mo_e1 += mo1[:,occidx] * lib.direct_sum('i-j->ij', e_i, e_i) - mo_e1 += v1mo[:,occidx,:] + mo_e1 = hs[:,occidx,:] + mo_e1 += mo1[:,occidx] * (e_i[:,None] - e_i) if h1.ndim == 3: return mo1, mo_e1 diff --git a/pyscf/scf/ucphf.py b/pyscf/scf/ucphf.py index cf4a7e641a..1b794b3d0e 100644 --- a/pyscf/scf/ucphf.py +++ b/pyscf/scf/ucphf.py @@ -27,7 +27,8 @@ def solve(fvind, mo_energy, mo_occ, h1, s1=None, - max_cycle=20, tol=1e-9, hermi=False, verbose=logger.WARN): + max_cycle=50, tol=1e-9, hermi=False, verbose=logger.WARN, + level_shift=0): ''' Args: fvind : function @@ -35,16 +36,18 @@ def solve(fvind, mo_energy, mo_occ, h1, s1=None, ''' if s1 is None: return solve_nos1(fvind, mo_energy, mo_occ, h1, - max_cycle, tol, hermi, verbose) + max_cycle, tol, hermi, verbose, level_shift) else: return solve_withs1(fvind, mo_energy, mo_occ, h1, s1, - max_cycle, tol, hermi, verbose) + max_cycle, tol, hermi, verbose, level_shift) kernel = solve # h1 shape is (:,nvir,nocc) def solve_nos1(fvind, mo_energy, mo_occ, h1, - max_cycle=20, tol=1e-9, hermi=False, verbose=logger.WARN): + max_cycle=20, tol=1e-9, hermi=False, verbose=logger.WARN, + level_shift=0): '''For field independent basis. First order overlap matrix is zero''' + assert not hermi log = logger.new_logger(verbose=verbose) t0 = (logger.process_clock(), logger.perf_counter()) @@ -56,15 +59,20 @@ def solve_nos1(fvind, mo_energy, mo_occ, h1, noccb = numpy.count_nonzero(occidxb) nvira = mo_occ[0].size - nocca nvirb = mo_occ[1].size - noccb - e_ai = numpy.hstack(((mo_energy[0][viridxa,None]-mo_energy[0][occidxa]).ravel(), - (mo_energy[1][viridxb,None]-mo_energy[1][occidxb]).ravel())) + mo_ea, mo_eb = mo_energy + e_ai = numpy.hstack( + ((mo_ea[viridxa,None]+level_shift - mo_ea[occidxa]).ravel(), + (mo_eb[viridxb,None]+level_shift - mo_eb[occidxb]).ravel())) e_ai = 1 / e_ai mo1base = numpy.hstack((h1[0].reshape(-1,nvira*nocca), h1[1].reshape(-1,nvirb*noccb))) mo1base *= -e_ai def vind_vo(mo1): - v = fvind(mo1.reshape(mo1base.shape)).reshape(mo1base.shape) + mo1 = mo1.reshape(mo1base.shape) + v = fvind(mo1).reshape(mo1base.shape) + if level_shift != 0: + v -= mo1 * level_shift v *= e_ai return v.ravel() mo1 = lib.krylov(vind_vo, mo1base.ravel(), @@ -83,12 +91,14 @@ def vind_vo(mo1): # h1 shape is (:,nvir+nocc,nocc) def solve_withs1(fvind, mo_energy, mo_occ, h1, s1, - max_cycle=20, tol=1e-9, hermi=False, verbose=logger.WARN): + max_cycle=20, tol=1e-9, hermi=False, verbose=logger.WARN, + level_shift=0): '''For field dependent basis. First order overlap matrix is non-zero. The first order orbitals are set to C^1_{ij} = -1/2 S1 e1 = h1 - s1*e0 + (e0_j-e0_i)*c1 + vhf[c1] ''' + assert not hermi log = logger.new_logger(verbose=verbose) t0 = (logger.process_clock(), logger.perf_counter()) @@ -99,27 +109,31 @@ def solve_withs1(fvind, mo_energy, mo_occ, h1, s1, nocca = numpy.count_nonzero(occidxa) noccb = numpy.count_nonzero(occidxb) nmoa, nmob = mo_occ[0].size, mo_occ[1].size - eai_a = mo_energy[0][viridxa,None] - mo_energy[0][occidxa] - eai_b = mo_energy[1][viridxb,None] - mo_energy[1][occidxb] + ei_a = mo_energy[0][occidxa] + ei_b = mo_energy[1][occidxb] + ea_a = mo_energy[0][viridxa] + ea_b = mo_energy[1][viridxb] + eai_a = 1. / (ea_a[:,None] + level_shift - ei_a) + eai_b = 1. / (ea_b[:,None] + level_shift - ei_b) s1_a = s1[0].reshape(-1,nmoa,nocca) nset = s1_a.shape[0] s1_b = s1[1].reshape(nset,nmob,noccb) - hs_a = mo1base_a = h1[0].reshape(nset,nmoa,nocca) - s1_a * mo_energy[0][occidxa] - hs_b = mo1base_b = h1[1].reshape(nset,nmob,noccb) - s1_b * mo_energy[1][occidxb] - mo_e1_a = hs_a[:,occidxa].copy() - mo_e1_b = hs_b[:,occidxb].copy() + hs_a = h1[0].reshape(nset,nmoa,nocca) - s1_a * ei_a + hs_b = h1[1].reshape(nset,nmob,noccb) - s1_b * ei_b - mo1base_a[:,viridxa]/= -eai_a - mo1base_b[:,viridxb]/= -eai_b + mo1base_a = hs_a.copy() + mo1base_b = hs_b.copy() + mo1base_a[:,viridxa] *= -eai_a + mo1base_b[:,viridxb] *= -eai_b mo1base_a[:,occidxa] = -s1_a[:,occidxa] * .5 mo1base_b[:,occidxb] = -s1_b[:,occidxb] * .5 - - eai_a = 1. / eai_a - eai_b = 1. / eai_b mo1base = numpy.hstack((mo1base_a.reshape(nset,-1), mo1base_b.reshape(nset,-1))) def vind_vo(mo1): + mo1 = mo1.reshape(mo1base.shape) v = fvind(mo1).reshape(mo1base.shape) + if level_shift != 0: + v -= mo1 * level_shift v1a = v[:,:nmoa*nocca].reshape(nset,nmoa,nocca) v1b = v[:,nmoa*nocca:].reshape(nset,nmob,noccb) v1a[:,viridxa] *= eai_a @@ -129,21 +143,27 @@ def vind_vo(mo1): return v.ravel() mo1 = lib.krylov(vind_vo, mo1base.ravel(), tol=tol, max_cycle=max_cycle, hermi=hermi, verbose=log) + mo1 = mo1.reshape(mo1base.shape) + mo1_a = mo1[:,:nmoa*nocca].reshape(nset,nmoa,nocca) + mo1_b = mo1[:,nmoa*nocca:].reshape(nset,nmob,noccb) + mo1_a[:,occidxa] = mo1base_a[:,occidxa] + mo1_b[:,occidxb] = mo1base_b[:,occidxb] log.timer('krylov solver in CPHF', *t0) v1mo = fvind(mo1).reshape(mo1base.shape) v1a = v1mo[:,:nmoa*nocca].reshape(nset,nmoa,nocca) v1b = v1mo[:,nmoa*nocca:].reshape(nset,nmob,noccb) - mo1 = mo1.reshape(mo1base.shape) - mo1_a = mo1[:,:nmoa*nocca].reshape(nset,nmoa,nocca) - mo1_b = mo1[:,nmoa*nocca:].reshape(nset,nmob,noccb) - mo1_a[:,viridxa] = mo1base_a[:,viridxa] - v1a[:,viridxa] * eai_a - mo1_b[:,viridxb] = mo1base_b[:,viridxb] - v1b[:,viridxb] * eai_b - mo_e1_a += mo1_a[:,occidxa] * (mo_energy[0][occidxa,None] - mo_energy[0][occidxa]) - mo_e1_b += mo1_b[:,occidxb] * (mo_energy[1][occidxb,None] - mo_energy[1][occidxb]) - mo_e1_a += v1mo[:,:nmoa*nocca].reshape(nset,nmoa,nocca)[:,occidxa] - mo_e1_b += v1mo[:,nmoa*nocca:].reshape(nset,nmob,noccb)[:,occidxb] + v1mo = fvind(mo1).reshape(mo1base.shape) + hs_a += v1mo[:,:nmoa*nocca].reshape(nset,nmoa,nocca) + hs_b += v1mo[:,nmoa*nocca:].reshape(nset,nmob,noccb) + mo1_a[:,viridxa] = hs_a[:,viridxa] / (ei_a - ea_a[:,None]) + mo1_b[:,viridxb] = hs_b[:,viridxb] / (ei_b - ea_b[:,None]) + + mo_e1_a = hs_a[:,occidxa] + mo_e1_b = hs_b[:,occidxb] + mo_e1_a += mo1_a[:,occidxa] * (ei_a[:,None] - ei_a) + mo_e1_b += mo1_b[:,occidxb] * (ei_b[:,None] - ei_b) if isinstance(h1[0], numpy.ndarray) and h1[0].ndim == 2: mo1_a, mo1_b = mo1_a[0], mo1_b[0]