From d7d062d4a4ebbf77b848668fcd0f67106fe61a06 Mon Sep 17 00:00:00 2001 From: Qiming Sun Date: Sat, 31 Aug 2024 21:20:12 -0700 Subject: [PATCH] Optimize nuclear hessian (#2241) * Remove chkfile dependency in hessian * Update krylov solver * Update cphf * Update hessian code * Update eph hessian * Fix a shape bug in ucphf * Remove duplicated code * clean up * Adjust tests * Update docstring of lib.linalg_helper.solve --- examples/eph/00-simple_eph.py | 24 ++++ pyscf/df/hessian/rhf.py | 12 +- pyscf/df/hessian/rks.py | 8 +- pyscf/df/hessian/uhf.py | 16 +-- pyscf/df/hessian/uks.py | 9 +- pyscf/eph/rhf.py | 42 +------ pyscf/eph/rks.py | 4 - pyscf/eph/uhf.py | 8 +- pyscf/eph/uks.py | 8 +- pyscf/hessian/rhf.py | 129 +++------------------- pyscf/hessian/rks.py | 7 +- pyscf/hessian/test/test_rhf.py | 62 +++++++++++ pyscf/hessian/test/test_rks.py | 6 +- pyscf/hessian/test/test_uhf.py | 66 +++++++++++ pyscf/hessian/test/test_uks.py | 7 +- pyscf/hessian/uhf.py | 157 ++++----------------------- pyscf/hessian/uks.py | 9 +- pyscf/lib/linalg_helper.py | 124 ++++++++------------- pyscf/lib/test/test_linalg_helper.py | 33 ++++-- pyscf/scf/cphf.py | 64 +++-------- pyscf/scf/ucphf.py | 35 +++--- 21 files changed, 311 insertions(+), 519 deletions(-) diff --git a/examples/eph/00-simple_eph.py b/examples/eph/00-simple_eph.py index 7e629d30c1..3fd8098ebe 100644 --- a/examples/eph/00-simple_eph.py +++ b/examples/eph/00-simple_eph.py @@ -3,7 +3,31 @@ ''' A simple example to run EPH calculation. ''' + from pyscf import gto, dft, eph + +mol = gto.M() +mol.atom = [['O', [0.000000000000, -0.000000000775, 0.923671924285]], + ['H', [-0.000000000000, -1.432564848017, 2.125164039823]], + ['H', [0.000000000000, 1.432564848792, 2.125164035930]]] +mol.unit = 'Bohr' +mol.basis = 'sto3g' +mol.verbose = 4 +mol.build() # this is a pre-computed relaxed geometry + +mf = mol.RHF() +mf.conv_tol = 1e-16 +mf.conv_tol_grad = 1e-10 +mf.kernel() + +myeph = eph.EPH(mf) +grad = mf.nuc_grad_method().kernel() +print("Force on the atoms/au:") +print(grad) +eph, omega = myeph.kernel(mo_rep=True) +print(np.amax(eph)) + + mol = gto.M(atom='N 0 0 0; N 0 0 2.100825', basis='def2-svp', verbose=4, unit="bohr") # this is a pre-computed relaxed molecule # for geometry relaxation, refer to pyscf/example/geomopt diff --git a/pyscf/df/hessian/rhf.py b/pyscf/df/hessian/rhf.py index e55cc993c6..aab57dc190 100644 --- a/pyscf/df/hessian/rhf.py +++ b/pyscf/df/hessian/rhf.py @@ -368,16 +368,8 @@ def make_h1(hessobj, mo_coeff, mo_occ, chkfile=None, atmlst=None, verbose=None): for ia, h1, vj1, vk1 in _gen_jk(hessobj, mo_coeff, mo_occ, chkfile, atmlst, verbose, True): h1 += vj1 - vk1 * .5 - - if chkfile is None: - h1ao[ia] = h1 - else: - key = 'scf_f1ao/%d' % ia - lib.chkfile.save(chkfile, key, h1) - if chkfile is None: - return h1ao - else: - return chkfile + h1ao[ia] = h1 + return h1ao def _gen_jk(hessobj, mo_coeff, mo_occ, chkfile=None, atmlst=None, verbose=None, with_k=True): diff --git a/pyscf/df/hessian/rks.py b/pyscf/df/hessian/rks.py index d526a8e035..e8c83fd87b 100644 --- a/pyscf/df/hessian/rks.py +++ b/pyscf/df/hessian/rks.py @@ -109,13 +109,7 @@ def make_h1(hessobj, mo_coeff, mo_occ, chkfile=None, atmlst=None, verbose=None): for ia, h1, vj1, vk1 in df_rhf_hess._gen_jk( hessobj, mo_coeff, mo_occ, chkfile, atmlst, verbose): h1ao[ia] -= .5 * (alpha - hyb) * vk1 - - if chkfile is None: - return h1ao - else: - for ia in atmlst: - lib.chkfile.save(chkfile, 'scf_f1ao/%d'%ia, h1ao[ia]) - return chkfile + return h1ao class Hessian(rks_hess.Hessian): diff --git a/pyscf/df/hessian/uhf.py b/pyscf/df/hessian/uhf.py index dc264350f6..42560799d9 100644 --- a/pyscf/df/hessian/uhf.py +++ b/pyscf/df/hessian/uhf.py @@ -398,19 +398,9 @@ def make_h1(hessobj, mo_coeff, mo_occ, chkfile=None, atmlst=None, verbose=None): for ia, h1, vj1, vk1 in _gen_jk(hessobj, mo_coeff, mo_occ, chkfile, atmlst, verbose, True): vk1a, vk1b = vk1 - h1a = h1 + vj1 - vk1a - h1b = h1 + vj1 - vk1b - - if chkfile is None: - h1aoa[ia] = h1a - h1aob[ia] = h1b - else: - lib.chkfile.save(chkfile, 'scf_f1ao/0/%d' % ia, h1a) - lib.chkfile.save(chkfile, 'scf_f1ao/1/%d' % ia, h1b) - if chkfile is None: - return (h1aoa,h1aob) - else: - return chkfile + h1aoa[ia] = h1 + vj1 - vk1a + h1aob[ia] = h1 + vj1 - vk1b + return (h1aoa,h1aob) def _gen_jk(hessobj, mo_coeff, mo_occ, chkfile=None, atmlst=None, verbose=None, with_k=True): diff --git a/pyscf/df/hessian/uks.py b/pyscf/df/hessian/uks.py index 853bf84526..0b7b905ce6 100644 --- a/pyscf/df/hessian/uks.py +++ b/pyscf/df/hessian/uks.py @@ -121,14 +121,7 @@ def make_h1(hessobj, mo_coeff, mo_occ, chkfile=None, atmlst=None, verbose=None): vk1a, vk1b = vk1 h1aoa[ia] -= (alpha - hyb) * vk1a h1aob[ia] -= (alpha - hyb) * vk1b - - if chkfile is None: - return h1aoa, h1aob - else: - for ia in atmlst: - lib.chkfile.save(chkfile, 'scf_f1ao/0/%d'%ia, h1aoa[ia]) - lib.chkfile.save(chkfile, 'scf_f1ao/1/%d'%ia, h1aob[ia]) - return chkfile + return h1aoa, h1aob class Hessian(uks_hess.Hessian): diff --git a/pyscf/eph/rhf.py b/pyscf/eph/rhf.py index c78e86b3fb..d6817561fb 100644 --- a/pyscf/eph/rhf.py +++ b/pyscf/eph/rhf.py @@ -23,7 +23,7 @@ import numpy as np import scipy.linalg from pyscf.hessian import rhf -from pyscf.lib import logger, chkfile +from pyscf.lib import logger from pyscf.scf._response_functions import _gen_rhf_response from pyscf import __config__ from pyscf.data.nist import HARTREE2WAVENUMBER, MP_ME @@ -37,16 +37,16 @@ def kernel(ephobj, mo_energy=None, mo_coeff=None, mo_occ=None, mo_rep=False): if mo_coeff is None: mo_coeff = ephobj.base.mo_coeff if mo_occ is None: mo_occ = ephobj.base.mo_occ - # chkfile is used to pass first orbitals from hessian methods to eph methods - # TODO: Remove the dependence to chfile and return first orbitals from a function - assert ephobj.chkfile is not None, 'chkfile is required to save first order orbitals' + h1ao = ephobj.make_h1(mo_coeff, mo_occ) + mo1, mo_e1 = ephobj.solve_mo1(mo_energy, mo_coeff, mo_occ, h1ao) - de = ephobj.hess_elec(mo_energy, mo_coeff, mo_occ) + de = ephobj.hess_elec(mo_energy, mo_coeff, mo_occ, + mo1=mo1, mo_e1=mo_e1, h1ao=h1ao) ephobj.de = de + ephobj.hess_nuc(ephobj.mol) omega, vec = ephobj.get_mode(ephobj.mol, ephobj.de) ephobj.omega, ephobj.vec = omega, vec - ephobj.eph = ephobj.get_eph(ephobj.chkfile, omega, vec, mo_rep) + ephobj.eph = ephobj.get_eph(mo1, omega, vec, mo_rep) return ephobj.eph, ephobj.omega def solve_hmat(mol, hmat, cutoff_frequency=CUTOFF_FREQUENCY, @@ -143,10 +143,6 @@ def _freq_mass_weighted_vec(vec, omega, mass): return vec def get_eph(ephobj, mo1, omega, vec, mo_rep): - if isinstance(mo1, str): - mo1 = chkfile.load(mo1, 'scf_mo1') - mo1 = {int(k): mo1[k] for k in mo1} - mol = ephobj.mol mf = ephobj.base vnuc_deriv = ephobj.vnuc_generator(mol) @@ -209,29 +205,3 @@ def __init__(self, scf_method, cutoff_frequency=CUTOFF_FREQUENCY, get_eph = get_eph vnuc_generator = vnuc_generator kernel = kernel - -if __name__ == '__main__': - from pyscf import gto, scf - - mol = gto.M() - mol.atom = [['O', [0.000000000000, -0.000000000775, 0.923671924285]], - ['H', [-0.000000000000, -1.432564848017, 2.125164039823]], - ['H', [0.000000000000, 1.432564848792, 2.125164035930]]] - mol.unit = 'Bohr' - mol.basis = 'sto3g' - mol.verbose=4 - mol.build() # this is a pre-computed relaxed geometry - - mf = scf.RHF(mol) - mf.conv_tol = 1e-16 - mf.conv_tol_grad = 1e-10 - mf.kernel() - - myeph = EPH(mf) - - grad = mf.nuc_grad_method().kernel() - print("Force on the atoms/au:") - print(grad) - - eph, omega = myeph.kernel(mo_rep=True) - print(np.amax(eph)) diff --git a/pyscf/eph/rks.py b/pyscf/eph/rks.py index 997d649331..6ab473f5af 100644 --- a/pyscf/eph/rks.py +++ b/pyscf/eph/rks.py @@ -100,10 +100,6 @@ def _get_vxc_deriv1(hessobj, mo_coeff, mo_occ, max_memory): return vmat def get_eph(ephobj, mo1, omega, vec, mo_rep): - if isinstance(mo1, str): - mo1 = lib.chkfile.load(mo1, 'scf_mo1') - mo1 = {int(k): mo1[k] for k in mo1} - mol = ephobj.mol mf = ephobj.base ni = mf._numint diff --git a/pyscf/eph/uhf.py b/pyscf/eph/uhf.py index ae4109a1a3..94cbe76594 100644 --- a/pyscf/eph/uhf.py +++ b/pyscf/eph/uhf.py @@ -54,18 +54,12 @@ def fx(mo1): return fx def get_eph(ephobj, mo1, omega, vec, mo_rep): - if isinstance(mo1, str): - mo1 = lib.chkfile.load(mo1, 'scf_mo1') - mo1a = mo1['0'] - mo1b = mo1['1'] - mo1a = {int(k): mo1a[k] for k in mo1a} - mo1b = {int(k): mo1b[k] for k in mo1b} - mol = ephobj.mol mf = ephobj.base vnuc_deriv = ephobj.vnuc_generator(mol) aoslices = mol.aoslice_by_atom() + mo1a, mo1b = mo1 mo_coeff, mo_occ = mf.mo_coeff, mf.mo_occ vind = uhf_deriv_generator(mf, mf.mo_coeff, mf.mo_occ) nao, nmo = mo_coeff[0].shape diff --git a/pyscf/eph/uks.py b/pyscf/eph/uks.py index 6256432e19..7065bab670 100644 --- a/pyscf/eph/uks.py +++ b/pyscf/eph/uks.py @@ -124,13 +124,6 @@ def _get_vxc_deriv1(hessobj, mo_coeff, mo_occ, max_memory): return vmata, vmatb def get_eph(ephobj, mo1, omega, vec, mo_rep): - if isinstance(mo1, str): - mo1 = lib.chkfile.load(mo1, 'scf_mo1') - mo1a = mo1['0'] - mo1b = mo1['1'] - mo1a = {int(k): mo1a[k] for k in mo1a} - mo1b = {int(k): mo1b[k] for k in mo1b} - mol = ephobj.mol mf = ephobj.base ni = mf._numint @@ -141,6 +134,7 @@ def get_eph(ephobj, mo1, omega, vec, mo_rep): vnuc_deriv = ephobj.vnuc_generator(mol) aoslices = mol.aoslice_by_atom() + mo1a, mo1b = mo1 mo_coeff, mo_occ = mf.mo_coeff, mf.mo_occ vind = uhf_deriv_generator(mf, mf.mo_coeff, mf.mo_occ) mem_now = lib.current_memory()[0] diff --git a/pyscf/hessian/rhf.py b/pyscf/hessian/rhf.py index c81bc4829a..175eade2db 100644 --- a/pyscf/hessian/rhf.py +++ b/pyscf/hessian/rhf.py @@ -53,20 +53,13 @@ def hess_elec(hessobj, mo_energy=None, mo_coeff=None, mo_occ=None, max_memory, log) if h1ao is None: - h1ao = hessobj.make_h1(mo_coeff, mo_occ, hessobj.chkfile, atmlst, log) + h1ao = hessobj.make_h1(mo_coeff, mo_occ, None, atmlst, log) t1 = log.timer_debug1('making H1', *time0) if mo1 is None or mo_e1 is None: mo1, mo_e1 = hessobj.solve_mo1(mo_energy, mo_coeff, mo_occ, h1ao, None, atmlst, max_memory, log) t1 = log.timer_debug1('solving MO1', *t1) - if isinstance(h1ao, str): - h1ao = lib.chkfile.load(h1ao, 'scf_f1ao') - h1ao = {int(k): h1ao[k] for k in h1ao} - if isinstance(mo1, str): - mo1 = lib.chkfile.load(mo1, 'scf_mo1') - mo1 = {int(k): mo1[k] for k in mo1} - nao, nmo = mo_coeff.shape mocc = mo_coeff[:,mo_occ>0] s1a = -mol.intor('int1e_ipovlp', comp=3) @@ -240,16 +233,8 @@ def make_h1(hessobj, mo_coeff, mo_occ, chkfile=None, atmlst=None, verbose=None): vhf[:,p0:p1] += vj2 - vk2*.5 h1 = vhf + vhf.transpose(0,2,1) h1 += hcore_deriv(ia) - - if chkfile is None: - h1ao[ia] = h1 - else: - key = 'scf_f1ao/%d' % ia - lib.chkfile.save(chkfile, key, h1) - if chkfile is None: - return h1ao - else: - return chkfile + h1ao[ia] = h1 + return h1ao def get_hcore(mol): '''Part of the second derivatives of core Hamiltonian''' @@ -298,7 +283,7 @@ def _get_jk(mol, intor, comp, aosym, script_dms, lib.hermi_triu(vs[k], hermi=hermi, inplace=True) return vs -def solve_mo1(mf, mo_energy, mo_coeff, mo_occ, h1ao_or_chkfile, +def solve_mo1(mf, mo_energy, mo_coeff, mo_occ, h1ao, fx=None, atmlst=None, max_memory=4000, verbose=None, max_cycle=50, level_shift=0): '''Solve the first order equation @@ -338,12 +323,7 @@ def _ao2mo(mat): s1ao[:,p0:p1] += s1a[:,p0:p1] s1ao[:,:,p0:p1] += s1a[:,p0:p1].transpose(0,2,1) s1vo.append(_ao2mo(s1ao)) - if isinstance(h1ao_or_chkfile, str): - key = 'scf_f1ao/%d' % ia - h1ao = lib.chkfile.load(h1ao_or_chkfile, key) - else: - h1ao = h1ao_or_chkfile[ia] - h1vo.append(_ao2mo(h1ao)) + h1vo.append(_ao2mo(h1ao[ia])) h1vo = numpy.vstack(h1vo) s1vo = numpy.vstack(s1vo) @@ -355,18 +335,10 @@ def _ao2mo(mat): for k in range(ia1-ia0): ia = atmlst[k+ia0] - if isinstance(h1ao_or_chkfile, str): - key = 'scf_mo1/%d' % ia - lib.chkfile.save(h1ao_or_chkfile, key, mo1[k]) - else: - mo1s[ia] = mo1[k] + mo1s[ia] = mo1[k] e1s[ia] = e1[k].reshape(3,nocc,nocc) mo1 = e1 = None - - if isinstance(h1ao_or_chkfile, str): - return h1ao_or_chkfile, e1s - else: - return mo1s, e1s + return mo1s, e1s def gen_vind(mf, mo_coeff, mo_occ): nao, nmo = mo_coeff.shape @@ -431,8 +403,7 @@ def gen_hop(hobj, mo_energy=None, mo_coeff=None, mo_occ=None, verbose=None): max_memory, log) de2 += hobj.hess_nuc() - # Compute H1 integrals and store in hobj.chkfile - hobj.make_h1(mo_coeff, mo_occ, hobj.chkfile, atmlst, log) + h1ao_cache = hobj.make_h1(mo_coeff, mo_occ, None, atmlst, log) aoslices = mol.aoslice_by_atom() s1a = -mol.intor('int1e_ipovlp', comp=3) @@ -445,8 +416,7 @@ def h_op(x): s1ao = 0 for ia in range(natm): shl0, shl1, p0, p1 = aoslices[ia] - h1ao_i = lib.chkfile.load(hobj.chkfile, 'scf_f1ao/%d' % ia) - h1ao += numpy.einsum('x,xij->ij', x[ia], h1ao_i) + h1ao += numpy.einsum('x,xij->ij', x[ia], h1ao_cache[ia]) s1ao_i = numpy.zeros((3,nao,nao)) s1ao_i[:,p0:p1] += s1a[:,p0:p1] s1ao_i[:,:,p0:p1] += s1a[:,p0:p1].transpose(0,2,1) @@ -463,8 +433,7 @@ def h_op(x): for ja in range(natm): q0, q1 = aoslices[ja][2:] - h1ao = lib.chkfile.load(hobj.chkfile, 'scf_f1ao/%s'%ja) - hx[ja] += numpy.einsum('xpq,pq->x', h1ao, dm1) * 4 + hx[ja] += numpy.einsum('xpq,pq->x', h1ao_cache[ja], dm1) * 4 hx[ja] -= numpy.einsum('xpq,pq->x', s1a[:,q0:q1], dme1[q0:q1]) * 2 hx[ja] -= numpy.einsum('xpq,qp->x', s1a[:,q0:q1], dme1[:,q0:q1]) * 2 return hx.ravel() @@ -484,16 +453,15 @@ class HessianBase(lib.StreamObject): level_shift = 0 _keys = { - 'mol', 'base', 'chkfile', 'atmlst', 'de', 'max_cycle', 'level_shift' + 'mol', 'base', 'atmlst', 'de', 'max_cycle', 'level_shift' } def __init__(self, scf_method): self.verbose = scf_method.verbose self.stdout = scf_method.stdout self.mol = scf_method.mol - self.chkfile = scf_method.chkfile - self.max_memory = self.mol.max_memory self.base = scf_method + self.max_memory = self.mol.max_memory self.atmlst = range(self.mol.natm) self.de = numpy.zeros((0,0,3,3)) # (A,B,dR_A,dR_B) @@ -575,9 +543,9 @@ def get_hcore(iatm, jatm): return hcore + hcore.conj().transpose(0,1,3,2) return get_hcore - def solve_mo1(self, mo_energy, mo_coeff, mo_occ, h1ao_or_chkfile, + def solve_mo1(self, mo_energy, mo_coeff, mo_occ, h1ao, fx=None, atmlst=None, max_memory=4000, verbose=None): - return solve_mo1(self.base, mo_energy, mo_coeff, mo_occ, h1ao_or_chkfile, + return solve_mo1(self.base, mo_energy, mo_coeff, mo_occ, h1ao, fx, atmlst, max_memory, verbose, max_cycle=self.max_cycle, level_shift=self.level_shift) @@ -622,72 +590,3 @@ class Hessian(HessianBase): from pyscf import scf scf.hf.RHF.Hessian = lib.class_as_method(Hessian) scf.rohf.ROHF.Hessian = lib.invalid_method('Hessian') - - -if __name__ == '__main__': - from pyscf import scf - - mol = gto.Mole() - mol.verbose = 0 - mol.output = None - mol.atom = [ - [1 , (1. , 0. , 0.000)], - [1 , (0. , 1. , 0.000)], - [1 , (0. , -1.517 , 1.177)], - [1 , (0. , 1.517 , 1.177)] ] - mol.basis = '631g' - mol.unit = 'B' - mol.build() - mf = scf.RHF(mol) - mf.conv_tol = 1e-14 - mf.scf() - n3 = mol.natm * 3 - hobj = mf.Hessian() - e2 = hobj.kernel().transpose(0,2,1,3).reshape(n3,n3) - print(lib.fp(e2) - -0.50693144355876429) - #from hessian import rhf_o0 - #e2ref = rhf_o0.Hessian(mf).kernel().transpose(0,2,1,3).reshape(n3,n3) - #print numpy.linalg.norm(e2-e2ref) - #print numpy.allclose(e2,e2ref) - - def grad_full(ia, inc): - coord = mol.atom_coord(ia).copy() - ptr = mol._atm[ia,gto.PTR_COORD] - de = [] - for i in range(3): - mol._env[ptr+i] = coord[i] + inc - mf = scf.RHF(mol).run(conv_tol=1e-14) - e1a = mf.nuc_grad_method().kernel() - mol._env[ptr+i] = coord[i] - inc - mf = scf.RHF(mol).run(conv_tol=1e-14) - e1b = mf.nuc_grad_method().kernel() - mol._env[ptr+i] = coord[i] - de.append((e1a-e1b)/(2*inc)) - return de - e2ref = [grad_full(ia, .5e-4) for ia in range(mol.natm)] - e2ref = numpy.asarray(e2ref).reshape(n3,n3) - print(numpy.linalg.norm(e2-e2ref)) - print(abs(e2-e2ref).max()) - print(numpy.allclose(e2,e2ref,atol=1e-6)) - -# \partial^2 E / \partial R \partial R' - e2 = hobj.partial_hess_elec(mf.mo_energy, mf.mo_coeff, mf.mo_occ) - e2 += hobj.hess_nuc(mol) - e2 = e2.transpose(0,2,1,3).reshape(n3,n3) - def grad_partial_R(ia, inc): - coord = mol.atom_coord(ia).copy() - ptr = mol._atm[ia,gto.PTR_COORD] - de = [] - for i in range(3): - mol._env[ptr+i] = coord[i] + inc - e1a = mf.nuc_grad_method().kernel() - mol._env[ptr+i] = coord[i] - inc - e1b = mf.nuc_grad_method().kernel() - mol._env[ptr+i] = coord[i] - de.append((e1a-e1b)/(2*inc)) - return de - e2ref = [grad_partial_R(ia, .5e-4) for ia in range(mol.natm)] - e2ref = numpy.asarray(e2ref).reshape(n3,n3) - print(numpy.linalg.norm(e2-e2ref)) - print(abs(e2-e2ref).max()) - print(numpy.allclose(e2,e2ref,atol=1e-8)) diff --git a/pyscf/hessian/rks.py b/pyscf/hessian/rks.py index a76afaccb9..82dbdba731 100644 --- a/pyscf/hessian/rks.py +++ b/pyscf/hessian/rks.py @@ -161,12 +161,7 @@ def make_h1(hessobj, mo_coeff, mo_occ, chkfile=None, atmlst=None, verbose=None): h1ao[ia] += veff + veff.transpose(0,2,1) h1ao[ia] += hcore_deriv(ia) - if chkfile is None: - return h1ao - else: - for ia in atmlst: - lib.chkfile.save(chkfile, 'scf_f1ao/%d'%ia, h1ao[ia]) - return chkfile + return h1ao XX, XY, XZ = 4, 5, 6 YX, YY, YZ = 5, 7, 8 diff --git a/pyscf/hessian/test/test_rhf.py b/pyscf/hessian/test/test_rhf.py index 8d33874c09..fc99b668ca 100644 --- a/pyscf/hessian/test/test_rhf.py +++ b/pyscf/hessian/test/test_rhf.py @@ -100,6 +100,68 @@ def test_finite_diff_rhf_hess(self): e2 = g_scanner(pmol.set_geom_('O 0. 0. -.0001; 1 0. -0.757 0.587; 1 0. 0.757 0.587'))[1] self.assertAlmostEqual(abs(hess[0,:,2] - (e1-e2)/2e-4*lib.param.BOHR).max(), 0, 4) + def test_finite_diff_rhf_hess1(self): + mol = gto.Mole() + mol.verbose = 0 + mol.output = None + mol.atom = [ + [1 , (1. , 0. , 0.000)], + [1 , (0. , 1. , 0.000)], + [1 , (0. , -1.517 , 1.177)], + [1 , (0. , 1.517 , 1.177)] ] + mol.basis = '631g' + mol.unit = 'B' + mol.build() + mf = scf.RHF(mol) + mf.conv_tol = 1e-14 + mf.scf() + n3 = mol.natm * 3 + hobj = mf.Hessian() + e2 = hobj.kernel().transpose(0,2,1,3).reshape(n3,n3) + self.assertAlmostEqual(lib.fp(e2), -0.50693144355876429, 6) + #from hessian import rhf_o0 + #e2ref = rhf_o0.Hessian(mf).kernel().transpose(0,2,1,3).reshape(n3,n3) + #print numpy.linalg.norm(e2-e2ref) + #print numpy.allclose(e2,e2ref) + + def grad_full(ia, inc): + coord = mol.atom_coord(ia).copy() + ptr = mol._atm[ia,gto.PTR_COORD] + de = [] + for i in range(3): + mol._env[ptr+i] = coord[i] + inc + mf = scf.RHF(mol).run(conv_tol=1e-14) + e1a = mf.nuc_grad_method().kernel() + mol._env[ptr+i] = coord[i] - inc + mf = scf.RHF(mol).run(conv_tol=1e-14) + e1b = mf.nuc_grad_method().kernel() + mol._env[ptr+i] = coord[i] + de.append((e1a-e1b)/(2*inc)) + return de + e2ref = [grad_full(ia, .5e-4) for ia in range(mol.natm)] + e2ref = numpy.asarray(e2ref).reshape(n3,n3) + self.assertAlmostEqual(abs(e2-e2ref).max(), 0, 6) + + # \partial^2 E / \partial R \partial R' + e2 = hobj.partial_hess_elec(mf.mo_energy, mf.mo_coeff, mf.mo_occ) + e2 += hobj.hess_nuc(mol) + e2 = e2.transpose(0,2,1,3).reshape(n3,n3) + def grad_partial_R(ia, inc): + coord = mol.atom_coord(ia).copy() + ptr = mol._atm[ia,gto.PTR_COORD] + de = [] + for i in range(3): + mol._env[ptr+i] = coord[i] + inc + e1a = mf.nuc_grad_method().kernel() + mol._env[ptr+i] = coord[i] - inc + e1b = mf.nuc_grad_method().kernel() + mol._env[ptr+i] = coord[i] + de.append((e1a-e1b)/(2*inc)) + return de + e2ref = [grad_partial_R(ia, .5e-4) for ia in range(mol.natm)] + e2ref = numpy.asarray(e2ref).reshape(n3,n3) + self.assertAlmostEqual(abs(e2-e2ref).max(), 0, 8) + @unittest.skipIf(dftd3 is None, "requires the dftd3 library") def test_finite_diff_rhf_d3_hess(self): mf = scf.RHF(mol) diff --git a/pyscf/hessian/test/test_rks.py b/pyscf/hessian/test/test_rks.py index 0f375f7863..a7f376764c 100644 --- a/pyscf/hessian/test/test_rks.py +++ b/pyscf/hessian/test/test_rks.py @@ -109,7 +109,7 @@ def test_finite_diff_lda_hess(self): mf.conv_tol = 1e-14 e0 = mf.kernel() hess = mf.Hessian().kernel() - self.assertAlmostEqual(lib.fp(hess), -0.7828771346902333, 6) + self.assertAlmostEqual(lib.fp(hess), -0.7828771346902333, 4) g_scanner = mf.nuc_grad_method().as_scanner() pmol = mol.copy() @@ -124,7 +124,7 @@ def test_finite_diff_b3lyp_hess(self): mf.xc = 'b3lyp5' e0 = mf.kernel() hess = mf.Hessian().kernel() - self.assertAlmostEqual(lib.fp(hess), -0.7590878171493624, 6) + self.assertAlmostEqual(lib.fp(hess), -0.7590878171493624, 4) g_scanner = mf.nuc_grad_method().as_scanner() pmol = mol.copy() @@ -187,7 +187,7 @@ def test_finite_diff_wb97x_hess(self): mf.xc = 'wb97x' e0 = mf.kernel() hess = mf.Hessian().kernel() - self.assertAlmostEqual(lib.fp(hess), -0.7637876979690904, 6) + self.assertAlmostEqual(lib.fp(hess), -0.7637876979690904, 4) g_scanner = mf.nuc_grad_method().as_scanner() pmol = mol.copy() diff --git a/pyscf/hessian/test/test_uhf.py b/pyscf/hessian/test/test_uhf.py index 85ca8942fc..8b65361b7b 100644 --- a/pyscf/hessian/test/test_uhf.py +++ b/pyscf/hessian/test/test_uhf.py @@ -73,6 +73,72 @@ def test_finite_diff_uhf_hess(self): e2 = g_scanner(pmol.set_geom_('O 0. 0. -.0001; 1 0. -0.757 0.587; 1 0. 0.757 0.587'))[1] self.assertAlmostEqual(abs(hess[0,:,2] - (e1-e2)/2e-4*lib.param.BOHR).max(), 0, 4) + def test_finite_diff_uhf_hess1(self): + mol = gto.Mole() + mol.verbose = 0 + mol.output = None + mol.atom = [ + [1 , (1. , 0. , 0.000)], + [1 , (0. , 1. , 0.000)], + [1 , (0. , -1.517 , 1.177)], + [1 , (0. , 1.517 , 1.177)] ] + mol.basis = '631g' + mol.unit = 'B' + mol.build() + mf = scf.UHF(mol) + mf.conv_tol = 1e-14 + mf.scf() + n3 = mol.natm * 3 + hobj = mf.Hessian() + e2 = hobj.kernel().transpose(0,2,1,3).reshape(n3,n3) + self.assertAlmostEqual(lib.fp(e2), -0.50693144355876429, 6) + + mol.spin = 2 + mf = scf.UHF(mol) + mf.conv_tol = 1e-14 + mf.scf() + n3 = mol.natm * 3 + hobj = mf.Hessian() + e2 = hobj.kernel().transpose(0,2,1,3).reshape(n3,n3) + + def grad_full(ia, inc): + coord = mol.atom_coord(ia).copy() + ptr = mol._atm[ia,gto.PTR_COORD] + de = [] + for i in range(3): + mol._env[ptr+i] = coord[i] + inc + mf = scf.UHF(mol).run(conv_tol=1e-14) + e1a = mf.nuc_grad_method().kernel() + mol._env[ptr+i] = coord[i] - inc + mf = scf.UHF(mol).run(conv_tol=1e-14) + e1b = mf.nuc_grad_method().kernel() + mol._env[ptr+i] = coord[i] + de.append((e1a-e1b)/(2*inc)) + return de + e2ref = [grad_full(ia, .5e-3) for ia in range(mol.natm)] + e2ref = numpy.asarray(e2ref).reshape(n3,n3) + self.assertAlmostEqual(abs(e2-e2ref).max(), 0, 6) + +# \partial^2 E / \partial R \partial R' + e2 = hobj.partial_hess_elec(mf.mo_energy, mf.mo_coeff, mf.mo_occ) + e2 += hobj.hess_nuc(mol) + e2 = e2.transpose(0,2,1,3).reshape(n3,n3) + def grad_partial_R(ia, inc): + coord = mol.atom_coord(ia).copy() + ptr = mol._atm[ia,gto.PTR_COORD] + de = [] + for i in range(3): + mol._env[ptr+i] = coord[i] + inc + e1a = mf.nuc_grad_method().kernel() + mol._env[ptr+i] = coord[i] - inc + e1b = mf.nuc_grad_method().kernel() + mol._env[ptr+i] = coord[i] + de.append((e1a-e1b)/(2*inc)) + return de + e2ref = [grad_partial_R(ia, .5e-4) for ia in range(mol.natm)] + e2ref = numpy.asarray(e2ref).reshape(n3,n3) + self.assertAlmostEqual(abs(e2-e2ref).max(), 0, 8) + @unittest.skipIf(dftd3 is None, "requires the dftd3 library") def test_finite_diff_uhf_d3_hess(self): mf = scf.UHF(mol) diff --git a/pyscf/hessian/test/test_uks.py b/pyscf/hessian/test/test_uks.py index f6b5995dea..1ca25b3277 100644 --- a/pyscf/hessian/test/test_uks.py +++ b/pyscf/hessian/test/test_uks.py @@ -112,7 +112,7 @@ def test_finite_diff_lda_hess(self): mf.conv_tol = 1e-14 e0 = mf.kernel() hess = mf.Hessian().kernel() - self.assertAlmostEqual(lib.fp(hess), -0.8503072107510495, 6) + self.assertAlmostEqual(lib.fp(hess), -0.8503072107510495, 4) g_scanner = mf.nuc_grad_method().as_scanner() pmol = mol.copy() @@ -126,8 +126,9 @@ def test_finite_diff_b3lyp_hess(self): mf.conv_tol = 1e-14 mf.xc = 'b3lyp5' e0 = mf.kernel() + mf.conv_tol_cpscf = 1e-9 hess = mf.Hessian().kernel() - self.assertAlmostEqual(lib.fp(hess), -0.8208641727673912, 6) + self.assertAlmostEqual(lib.fp(hess), -0.8208641727673912, 5) g_scanner = mf.nuc_grad_method().as_scanner() pmol = mol.copy() @@ -174,7 +175,7 @@ def test_finite_diff_wb97x_hess(self): mf.xc = 'wb97x' e0 = mf.kernel() hess = mf.Hessian().kernel() - self.assertAlmostEqual(lib.fp(hess), -0.8207572641132195, 6) + self.assertAlmostEqual(lib.fp(hess), -0.8207572641132195, 4) g_scanner = mf.nuc_grad_method().as_scanner() pmol = mol.copy() diff --git a/pyscf/hessian/uhf.py b/pyscf/hessian/uhf.py index fc90246728..e93795ecd5 100644 --- a/pyscf/hessian/uhf.py +++ b/pyscf/hessian/uhf.py @@ -53,29 +53,15 @@ def hess_elec(hessobj, mo_energy=None, mo_coeff=None, mo_occ=None, max_memory, log) if h1ao is None: - h1ao = hessobj.make_h1(mo_coeff, mo_occ, hessobj.chkfile, atmlst, log) + h1ao = hessobj.make_h1(mo_coeff, mo_occ, None, atmlst, log) t1 = log.timer_debug1('making H1', *time0) if mo1 is None or mo_e1 is None: mo1, mo_e1 = hessobj.solve_mo1(mo_energy, mo_coeff, mo_occ, h1ao, None, atmlst, max_memory, log) t1 = log.timer_debug1('solving MO1', *t1) - if isinstance(h1ao, str): - h1ao = lib.chkfile.load(h1ao, 'scf_f1ao') - h1aoa = h1ao['0'] - h1aob = h1ao['1'] - h1aoa = {int(k): h1aoa[k] for k in h1aoa} - h1aob = {int(k): h1aob[k] for k in h1aob} - else: - h1aoa, h1aob = h1ao - if isinstance(mo1, str): - mo1 = lib.chkfile.load(mo1, 'scf_mo1') - mo1a = mo1['0'] - mo1b = mo1['1'] - mo1a = {int(k): mo1a[k] for k in mo1a} - mo1b = {int(k): mo1b[k] for k in mo1b} - else: - mo1a, mo1b = mo1 + h1aoa, h1aob = h1ao + mo1a, mo1b = mo1 mo_e1a, mo_e1b = mo_e1 nao, nmo = mo_coeff[0].shape @@ -245,21 +231,11 @@ def make_h1(hessobj, mo_coeff, mo_occ, chkfile=None, atmlst=None, verbose=None): vhfa[:,p0:p1] += vj2 - vk2a vhfb[:,p0:p1] += vj2 - vk2b h1 = hcore_deriv(ia) - h1a = h1 + vhfa + vhfa.transpose(0,2,1) - h1b = h1 + vhfb + vhfb.transpose(0,2,1) - - if chkfile is None: - h1aoa[ia] = h1a - h1aob[ia] = h1b - else: - lib.chkfile.save(chkfile, 'scf_f1ao/0/%d' % ia, h1a) - lib.chkfile.save(chkfile, 'scf_f1ao/1/%d' % ia, h1b) - if chkfile is None: - return (h1aoa,h1aob) - else: - return chkfile - -def solve_mo1(mf, mo_energy, mo_coeff, mo_occ, h1ao_or_chkfile, + h1aoa[ia] = h1 + vhfa + vhfa.transpose(0,2,1) + h1aob[ia] = h1 + vhfb + vhfb.transpose(0,2,1) + return (h1aoa,h1aob) + +def solve_mo1(mf, mo_energy, mo_coeff, mo_occ, h1ao, fx=None, atmlst=None, max_memory=4000, verbose=None, max_cycle=50, level_shift=0): mol = mf.mol @@ -274,6 +250,7 @@ def solve_mo1(mf, mo_energy, mo_coeff, mo_occ, h1ao_or_chkfile, if fx is None: fx = gen_vind(mf, mo_coeff, mo_occ) s1a = -mol.intor('int1e_ipovlp', comp=3) + h1aoa, h1aob = h1ao def _ao2mo(mat, mo_coeff, mocc): return numpy.asarray([reduce(numpy.dot, (mo_coeff.T, x, mocc)) for x in mat]) @@ -299,14 +276,8 @@ def _ao2mo(mat, mo_coeff, mocc): s1ao[:,:,p0:p1] += s1a[:,p0:p1].transpose(0,2,1) s1voa.append(_ao2mo(s1ao, mo_coeff[0], mocca)) s1vob.append(_ao2mo(s1ao, mo_coeff[1], moccb)) - if isinstance(h1ao_or_chkfile, str): - h1aoa = lib.chkfile.load(h1ao_or_chkfile, 'scf_f1ao/0/%d'%ia) - h1aob = lib.chkfile.load(h1ao_or_chkfile, 'scf_f1ao/1/%d'%ia) - else: - h1aoa = h1ao_or_chkfile[0][ia] - h1aob = h1ao_or_chkfile[1][ia] - h1voa.append(_ao2mo(h1aoa, mo_coeff[0], mocca)) - h1vob.append(_ao2mo(h1aob, mo_coeff[1], moccb)) + h1voa.append(_ao2mo(h1aoa[ia], mo_coeff[0], mocca)) + h1vob.append(_ao2mo(h1aob[ia], mo_coeff[1], moccb)) h1vo = (numpy.vstack(h1voa), numpy.vstack(h1vob)) s1vo = (numpy.vstack(s1voa), numpy.vstack(s1vob)) @@ -320,20 +291,13 @@ def _ao2mo(mat, mo_coeff, mocc): for k in range(ia1-ia0): ia = atmlst[k+ia0] - if isinstance(h1ao_or_chkfile, str): - lib.chkfile.save(h1ao_or_chkfile, 'scf_mo1/0/%d'%ia, mo1a[k]) - lib.chkfile.save(h1ao_or_chkfile, 'scf_mo1/1/%d'%ia, mo1b[k]) - else: - mo1sa[ia] = mo1a[k] - mo1sb[ia] = mo1b[k] + mo1sa[ia] = mo1a[k] + mo1sb[ia] = mo1b[k] e1sa[ia] = e1a[k].reshape(3,nocca,nocca) e1sb[ia] = e1b[k].reshape(3,noccb,noccb) mo1 = e1 = mo1a = mo1b = e1a = e1b = None - if isinstance(h1ao_or_chkfile, str): - return h1ao_or_chkfile, (e1sa,e1sb) - else: - return (mo1sa,mo1sb), (e1sa,e1sb) + return (mo1sa,mo1sb), (e1sa,e1sb) def gen_vind(mf, mo_coeff, mo_occ): nao, nmoa = mo_coeff[0].shape @@ -388,8 +352,8 @@ def gen_hop(hobj, mo_energy=None, mo_coeff=None, mo_occ=None, verbose=None): max_memory, log) de2 += hobj.hess_nuc() - # Compute H1 integrals and store in hobj.chkfile - hobj.make_h1(mo_coeff, mo_occ, hobj.chkfile, atmlst, log) + h1ao_cache = hobj.make_h1(mo_coeff, mo_occ, None, atmlst, log) + h1aoa_cache, h1aob_cache = h1ao_cache aoslices = mol.aoslice_by_atom() s1a = -mol.intor('int1e_ipovlp', comp=3) @@ -403,10 +367,8 @@ def h_op(x): s1ao = 0 for ia in range(natm): shl0, shl1, p0, p1 = aoslices[ia] - h1ao_i = lib.chkfile.load(hobj.chkfile, 'scf_f1ao/0/%d' % ia) - h1aoa += numpy.einsum('x,xij->ij', x[ia], h1ao_i) - h1ao_i = lib.chkfile.load(hobj.chkfile, 'scf_f1ao/1/%d' % ia) - h1aob += numpy.einsum('x,xij->ij', x[ia], h1ao_i) + h1aoa += numpy.einsum('x,xij->ij', x[ia], h1aoa_cache[ia]) + h1aob += numpy.einsum('x,xij->ij', x[ia], h1aob_cache[ia]) s1ao_i = numpy.zeros((3,nao,nao)) s1ao_i[:,p0:p1] += s1a[:,p0:p1] s1ao_i[:,:,p0:p1] += s1a[:,p0:p1].transpose(0,2,1) @@ -432,10 +394,8 @@ def h_op(x): for ja in range(natm): q0, q1 = aoslices[ja][2:] - h1aoa = lib.chkfile.load(hobj.chkfile, 'scf_f1ao/0/%d' % ja) - h1aob = lib.chkfile.load(hobj.chkfile, 'scf_f1ao/1/%d' % ja) - hx[ja] += numpy.einsum('xpq,pq->x', h1aoa, dm1a) * 2 - hx[ja] += numpy.einsum('xpq,pq->x', h1aob, dm1b) * 2 + hx[ja] += numpy.einsum('xpq,pq->x', h1aoa_cache[ja], dm1a) * 2 + hx[ja] += numpy.einsum('xpq,pq->x', h1aob_cache[ja], dm1b) * 2 hx[ja] -= numpy.einsum('xpq,pq->x', s1a[:,q0:q1], dme1[q0:q1]) hx[ja] -= numpy.einsum('xpq,qp->x', s1a[:,q0:q1], dme1[:,q0:q1]) return hx.ravel() @@ -452,84 +412,11 @@ class Hessian(rhf_hess.HessianBase): make_h1 = make_h1 gen_hop = gen_hop - def solve_mo1(self, mo_energy, mo_coeff, mo_occ, h1ao_or_chkfile, + def solve_mo1(self, mo_energy, mo_coeff, mo_occ, h1ao, fx=None, atmlst=None, max_memory=4000, verbose=None): - return solve_mo1(self.base, mo_energy, mo_coeff, mo_occ, h1ao_or_chkfile, + return solve_mo1(self.base, mo_energy, mo_coeff, mo_occ, h1ao, fx, atmlst, max_memory, verbose, max_cycle=self.max_cycle, level_shift=self.level_shift) from pyscf import scf scf.uhf.UHF.Hessian = lib.class_as_method(Hessian) - -if __name__ == '__main__': - from pyscf import gto - from pyscf import scf - - mol = gto.Mole() - mol.verbose = 0 - mol.output = None - mol.atom = [ - [1 , (1. , 0. , 0.000)], - [1 , (0. , 1. , 0.000)], - [1 , (0. , -1.517 , 1.177)], - [1 , (0. , 1.517 , 1.177)] ] - mol.basis = '631g' - mol.unit = 'B' - mol.build() - mf = scf.UHF(mol) - mf.conv_tol = 1e-14 - mf.scf() - n3 = mol.natm * 3 - hobj = mf.Hessian() - e2 = hobj.kernel().transpose(0,2,1,3).reshape(n3,n3) - print(lib.fp(e2) - -0.50693144355876429) - - mol.spin = 2 - mf = scf.UHF(mol) - mf.conv_tol = 1e-14 - mf.scf() - n3 = mol.natm * 3 - hobj = Hessian(mf) - e2 = hobj.kernel().transpose(0,2,1,3).reshape(n3,n3) - - def grad_full(ia, inc): - coord = mol.atom_coord(ia).copy() - ptr = mol._atm[ia,gto.PTR_COORD] - de = [] - for i in range(3): - mol._env[ptr+i] = coord[i] + inc - mf = scf.UHF(mol).run(conv_tol=1e-14) - e1a = mf.nuc_grad_method().kernel() - mol._env[ptr+i] = coord[i] - inc - mf = scf.UHF(mol).run(conv_tol=1e-14) - e1b = mf.nuc_grad_method().kernel() - mol._env[ptr+i] = coord[i] - de.append((e1a-e1b)/(2*inc)) - return de - e2ref = [grad_full(ia, .5e-3) for ia in range(mol.natm)] - e2ref = numpy.asarray(e2ref).reshape(n3,n3) - print(numpy.linalg.norm(e2-e2ref)) - print(abs(e2-e2ref).max()) - print(numpy.allclose(e2,e2ref,atol=1e-4)) - -# \partial^2 E / \partial R \partial R' - e2 = hobj.partial_hess_elec(mf.mo_energy, mf.mo_coeff, mf.mo_occ) - e2 += hobj.hess_nuc(mol) - e2 = e2.transpose(0,2,1,3).reshape(n3,n3) - def grad_partial_R(ia, inc): - coord = mol.atom_coord(ia).copy() - ptr = mol._atm[ia,gto.PTR_COORD] - de = [] - for i in range(3): - mol._env[ptr+i] = coord[i] + inc - e1a = mf.nuc_grad_method().kernel() - mol._env[ptr+i] = coord[i] - inc - e1b = mf.nuc_grad_method().kernel() - mol._env[ptr+i] = coord[i] - de.append((e1a-e1b)/(2*inc)) - return de - e2ref = [grad_partial_R(ia, .5e-4) for ia in range(mol.natm)] - e2ref = numpy.asarray(e2ref).reshape(n3,n3) - print(numpy.linalg.norm(e2-e2ref)) - print(abs(e2-e2ref).max()) - print(numpy.allclose(e2,e2ref,atol=1e-8)) diff --git a/pyscf/hessian/uks.py b/pyscf/hessian/uks.py index 525436ce43..74bdb85bac 100644 --- a/pyscf/hessian/uks.py +++ b/pyscf/hessian/uks.py @@ -192,14 +192,7 @@ def make_h1(hessobj, mo_coeff, mo_occ, chkfile=None, atmlst=None, verbose=None): h1aoa[ia] += h1 + veffa + veffa.transpose(0,2,1) h1aob[ia] += h1 + veffb + veffb.transpose(0,2,1) - if chkfile is None: - return h1aoa, h1aob - else: - for ia in atmlst: - lib.chkfile.save(chkfile, 'scf_f1ao/0/%d'%ia, h1aoa[ia]) - lib.chkfile.save(chkfile, 'scf_f1ao/1/%d'%ia, h1aob[ia]) - return chkfile - + return h1aoa, h1aob XX, XY, XZ = 4, 5, 6 YX, YY, YZ = 5, 7, 8 diff --git a/pyscf/lib/linalg_helper.py b/pyscf/lib/linalg_helper.py index 4a750c4ea5..59a78b5c68 100644 --- a/pyscf/lib/linalg_helper.py +++ b/pyscf/lib/linalg_helper.py @@ -34,7 +34,7 @@ SAFE_EIGH_LINDEP = getattr(__config__, 'lib_linalg_helper_safe_eigh_lindep', 1e-15) DAVIDSON_LINDEP = getattr(__config__, 'lib_linalg_helper_davidson_lindep', 1e-14) DSOLVE_LINDEP = getattr(__config__, 'lib_linalg_helper_dsolve_lindep', 1e-13) -MAX_MEMORY = getattr(__config__, 'lib_linalg_helper_davidson_max_memory', 2000) # 2GB +MAX_MEMORY = getattr(__config__, 'lib_linalg_helper_davidson_max_memory', 4000) # 4GB # sort by similarity has problem which flips the ordering of eigenvalues when # the initial guess is closed to excited state. In this situation, function @@ -1298,33 +1298,28 @@ def krylov(aop, b, x0=None, tol=1e-10, max_cycle=30, dot=numpy.dot, if isinstance(aop, numpy.ndarray) and aop.ndim == 2: return numpy.linalg.solve(aop+numpy.eye(aop.shape[0]), b) - if isinstance(verbose, logger.Logger): - log = verbose - else: - log = logger.Logger(sys.stdout, verbose) + log = logger.new_logger(verbose=verbose) if not (isinstance(b, numpy.ndarray) and b.ndim == 1): b = numpy.asarray(b) - if x0 is None: - x1 = b - else: + if x0 is not None: b = b - (x0 + aop(x0)) - x1 = b + x1 = b if x1.ndim == 1: x1 = x1.reshape(1, x1.size) nroots, ndim = x1.shape - # Not exactly QR, vectors are orthogonal but not normalized - x1, rmat = _qr(x1, dot, lindep) - for i in range(len(x1)): - x1[i] *= rmat[i,i] - - innerprod = [dot(xi.conj(), xi).real for xi in x1] - if innerprod: + if nroots > 1: + # Not exactly QR, vectors are orthogonal but not normalized + x1, rmat = _qr(x1, dot) + x1 *= rmat.diagonal()[:,None] + innerprod = (rmat.diagonal().real ** 2).tolist() max_innerprod = max(innerprod) else: - max_innerprod = 0 + max_innerprod = dot(x1[0].conj(), x1[0]).real + innerprod = [max_innerprod] + if max_innerprod < lindep or max_innerprod < tol**2: if x0 is None: return numpy.zeros_like(b) @@ -1353,27 +1348,34 @@ def krylov(aop, b, x0=None, tol=1e-10, max_cycle=30, dot=numpy.dot, x1 = axt.copy() for i in range(len(xs)): xsi = numpy.asarray(xs[i]) - for j, axj in enumerate(axt): - x1[j] -= xsi * (dot(xsi.conj(), axj) / innerprod[i]) + w = dot(axt, xsi.conj()) / innerprod[i] + x1 -= xsi * w[:,None] axt = None - max_innerprod = 0 - idx = [] - for i, xi in enumerate(x1): - innerprod1 = dot(xi.conj(), xi).real - max_innerprod = max(max_innerprod, innerprod1) - if innerprod1 > lindep and innerprod1 > tol**2: - idx.append(i) - innerprod.append(innerprod1) + if nroots > 1: + x1, rmat = _qr(x1, dot) + x1 *= rmat.diagonal()[:,None] + innerprod1 = rmat.diagonal().real ** 2 + else: + innerprod1 = [dot(x1[0].conj(), x1[0]).real] + max_innerprod = max(innerprod1, default=0.) + log.debug('krylov cycle %d r = %g', cycle, max_innerprod**.5) if max_innerprod < lindep or max_innerprod < tol**2: break - x1 = x1[idx] + if nroots > 1: + mask = (innerprod1 > lindep) & (innerprod1 > tol**2) + x1 = x1[mask] + innerprod.extend(innerprod1[mask]) + else: + innerprod.append(innerprod1[0]) - nd = cycle + 1 - h = numpy.empty((nd,nd), dtype=x1.dtype) + else: + raise RuntimeError("Krylov solver failed to converge.") + nd = len(xs) + h = numpy.empty((nd,nd), dtype=x1.dtype) for i in range(nd): xi = numpy.asarray(xs[i]) if hermi: @@ -1384,72 +1386,37 @@ def krylov(aop, b, x0=None, tol=1e-10, max_cycle=30, dot=numpy.dot, for j in range(nd): h[i,j] = dot(xi.conj(), ax[j]) xi = None - # Add the contribution of I in (1+a) for i in range(nd): h[i,i] += innerprod[i] g = numpy.zeros((nd,nroots), dtype=x1.dtype) - if b.ndim == 1: - g[0] = innerprod[0] - else: - # Restore the first nroots vectors, which are array b or b-(1+a)x0 - for i in range(min(nd, nroots)): - xsi = numpy.asarray(xs[i]) - for j in range(nroots): - g[i,j] = dot(xsi.conj(), b[j]) + for i in range(nd): + g[i] = dot(b, numpy.asarray(xs[i]).conj()) c = numpy.linalg.solve(h, g) x = _gen_x0(c, xs) if b.ndim == 1: x = x[0] + # Restore the first nroots vectors, which are array b or b-(1+a)x0 if x0 is not None: x += x0 return x -def solve(aop, b, precond, tol=1e-12, max_cycle=30, dot=numpy.dot, +def solve(aop, b, precond=None, tol=1e-12, max_cycle=60, dot=numpy.dot, lindep=DSOLVE_LINDEP, verbose=0, tol_residual=None): - '''Davidson iteration to solve linear equation. + '''Solve a linear equation using the GMRES solver from the scipy.sparse module. ''' - msg = ('linalg_helper.solve is a bad solver for linear equations. ' - 'You should not use this solver in product.') - warnings.warn(msg) - if tol_residual is None: - toloose = numpy.sqrt(tol) - else: - toloose = tol_residual - - assert callable(precond) - - xs = [precond(b)] - ax = [aop(xs[-1])] - - dtype = numpy.result_type(ax[0], xs[0]) - aeff = numpy.zeros((max_cycle,max_cycle), dtype=dtype) - beff = numpy.zeros((max_cycle), dtype=dtype) - for istep in range(max_cycle): - beff[istep] = dot(xs[istep].conj(), b) - for i in range(istep+1): - aeff[istep,i] = dot(xs[istep].conj(), ax[i]) - aeff[i,istep] = dot(xs[i].conj(), ax[istep]) - - v = scipy.linalg.solve(aeff[:istep+1,:istep+1], beff[:istep+1]) - xtrial = _outprod_to_subspace(v, xs) - dx = b - _outprod_to_subspace(v, ax) - rr = numpy_helper.norm(dx) - if verbose: - print('davidson', istep, rr) - if rr < toloose: - break - xs.append(precond(dx)) - ax.append(aop(xs[-1])) - - if verbose: - print(istep) - - return xtrial + from scipy.sparse.linalg import gmres, LinearOperator + assert b.ndim == 1 + n = b.size + c, info = gmres(LinearOperator((n,n), matvec=aop), b, maxiter=max_cycle, + atol=max(lindep, tol)**.5) + if info != 0: + raise RuntimeError(f'scipy.sparse fails to solve linear equation info={info}') + return c dsolve = solve @@ -1471,7 +1438,6 @@ def cho_solve(a, b, strict_sym_pos=True): (fname, lineno)) return scipy.linalg.solve(a, b) - def _qr(xs, dot, lindep=1e-14): '''QR decomposition for a list of vectors (for linearly independent vectors only). xs = (r.T).dot(qs) diff --git a/pyscf/lib/test/test_linalg_helper.py b/pyscf/lib/test/test_linalg_helper.py index efea3cb2c6..158938f606 100644 --- a/pyscf/lib/test/test_linalg_helper.py +++ b/pyscf/lib/test/test_linalg_helper.py @@ -93,12 +93,10 @@ def aop(x): self.assertEqual(e0.size, 1) self.assertAlmostEqual(e0, 2, 8) - @unittest.skip('bad solver. experimental only') def test_solve(self): numpy.random.seed(12) n = 100 a = numpy.random.rand(n,n) - a = a + a.conj().T a += numpy.diag(numpy.random.random(n))* 10 b = numpy.random.random(n) def aop(x): @@ -106,15 +104,8 @@ def aop(x): def precond(x, *args): return x / a.diagonal() xref = numpy.linalg.solve(a, b) - x1 = linalg_helper.dsolve(aop, b, precond, max_cycle=50) - self.assertAlmostEqual(abs(xref - x1).max(), 0, 6) - a_diag = a.diagonal() - aop = lambda x: numpy.dot(a-numpy.diag(a_diag), x.ravel())/a_diag - x1 = linalg_helper.krylov(aop, b/a_diag, max_cycle=50) - self.assertAlmostEqual(abs(xref - x1).max(), 0, 6) - x1 = linalg_helper.krylov(aop, b/a_diag, None, max_cycle=10) - x1 = linalg_helper.krylov(aop, b/a_diag, x1, max_cycle=30) - self.assertAlmostEqual(abs(xref - x1).max(), 0, 6) + x1 = linalg_helper.dsolve(aop, b, precond, max_cycle=80) + self.assertAlmostEqual(abs(xref - x1).max(), 0, 3) def test_krylov_with_level_shift(self): numpy.random.seed(10) @@ -136,6 +127,25 @@ def test_krylov_with_level_shift(self): c = linalg_helper.krylov(aop, b/a_diag, max_cycle=17) self.assertAlmostEqual(abs(ref - c).max(), 0, 8) + def test_krylov_multiple_roots(self): + numpy.random.seed(10) + n = 100 + a = numpy.random.rand(n,n) * .1 + b = numpy.random.rand(4, n) + ref = numpy.linalg.solve(numpy.eye(n) + a, b.T).T + + aop = lambda x: x.dot(a.T) + c = linalg_helper.krylov(aop, b, lindep=1e-14) + self.assertAlmostEqual(abs(ref - c).max(), 0, 7) + + a = numpy.random.rand(n,n) * .1 + numpy.random.rand(n,n) * .1j + b = numpy.random.rand(4, n) + numpy.random.rand(4, n) * .5j + ref = numpy.linalg.solve(numpy.eye(n) + a, b.T).T + + aop = lambda x: x.dot(a.T) + c = linalg_helper.krylov(aop, b, lindep=1e-14) + self.assertAlmostEqual(abs(ref - c).max(), 0, 7) + def test_dgeev(self): numpy.random.seed(12) n = 100 @@ -214,7 +224,6 @@ def abop(x): print((abs(vl[1]) - abs(ul[:,1])).max()) print((abs(vl[2]) - abs(ul[:,2])).max()) - @unittest.skip('difficult to converge') def test_eig_difficult_problem(self): N = 40 neig = 4 diff --git a/pyscf/scf/cphf.py b/pyscf/scf/cphf.py index cda7b94d56..1c91b7d288 100644 --- a/pyscf/scf/cphf.py +++ b/pyscf/scf/cphf.py @@ -68,15 +68,16 @@ def solve_nos1(fvind, mo_energy, mo_occ, h1, e_i = mo_energy[mo_occ>0] e_ai = 1 / (e_a[:,None] + level_shift - e_i) mo1base = h1 * -e_ai + nvir, nocc = e_ai.shape def vind_vo(mo1): - mo1 = mo1.reshape(h1.shape) - v = fvind(mo1).reshape(h1.shape) + mo1 = mo1.reshape(-1, nvir, nocc) + v = fvind(mo1).reshape(-1, nvir, nocc) if level_shift != 0: v -= mo1 * level_shift v *= e_ai - return v.ravel() - mo1 = lib.krylov(vind_vo, mo1base.ravel(), + return v.reshape(-1, nvir*nocc) + mo1 = lib.krylov(vind_vo, mo1base.reshape(-1, nvir*nocc), tol=tol, max_cycle=max_cycle, hermi=hermi, verbose=log) log.timer('krylov solver in CPHF', *t0) return mo1.reshape(h1.shape), None @@ -119,20 +120,20 @@ def solve_withs1(fvind, mo_energy, mo_occ, h1, s1, mo1base[:,occidx] = -s1[:,occidx] * .5 def vind_vo(mo1): - mo1 = mo1.reshape(mo1base.shape) - v = fvind(mo1).reshape(mo1base.shape) + mo1 = mo1.reshape(-1, nmo, nocc) + v = fvind(mo1).reshape(-1, nmo, nocc) 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(), + return v.reshape(-1, nmo*nocc) + mo1 = lib.krylov(vind_vo, mo1base.reshape(-1, nmo*nocc), tol=tol, max_cycle=max_cycle, hermi=hermi, verbose=log) - mo1 = mo1.reshape(mo1base.shape) + mo1 = mo1.reshape(-1, nmo, nocc) mo1[:,occidx] = mo1base[:,occidx] log.timer('krylov solver in CPHF', *t0) - hs += fvind(mo1).reshape(mo1base.shape) + hs += fvind(mo1).reshape(-1, nmo, nocc) mo1[:,viridx] = hs[:,viridx] / (e_i - e_a[:,None]) # mo_e1 has the same symmetry as the first order Fock matrix (hermitian or @@ -143,44 +144,5 @@ def vind_vo(mo1): if h1.ndim == 3: return mo1, mo_e1 else: - return mo1.reshape(h1.shape), mo_e1.reshape(nocc,nocc) - -if __name__ == '__main__': - numpy.random.seed(1) - nd = 3 - nocc = 5 - nmo = 12 - nvir = nmo - nocc - a = numpy.random.random((nocc*nvir,nocc*nvir)) - a = a + a.T - def fvind(x): - v = numpy.dot(a,x[:,nocc:].reshape(-1,nocc*nvir).T) - v1 = numpy.zeros((nd,nmo,nocc)) - v1[:,nocc:] = v.T.reshape(nd,nvir,nocc) - return v1 - mo_energy = numpy.sort(numpy.random.random(nmo)) * 10 - mo_occ = numpy.zeros(nmo) - mo_occ[:nocc] = 2 - e_i = mo_energy[mo_occ>0] - e_a = mo_energy[mo_occ==0] - e_ai = 1 / lib.direct_sum('a-i->ai', e_a, e_i) - h1 = numpy.random.random((nd,nmo,nocc)) - h1[:,:nocc,:nocc] = h1[:,:nocc,:nocc] + h1[:,:nocc,:nocc].transpose(0,2,1) - s1 = numpy.random.random((nd,nmo,nocc)) - s1[:,:nocc,:nocc] = s1[:,:nocc,:nocc] + s1[:,:nocc,:nocc].transpose(0,2,1) - - x = solve(fvind, mo_energy, mo_occ, h1, s1, max_cycle=30)[0] - print(numpy.linalg.norm(x)-6.272581531366389) - hs = h1.reshape(-1,nmo,nocc) - s1.reshape(-1,nmo,nocc)*e_i - print(abs(hs[:,nocc:] + fvind(x)[:,nocc:]+x[:,nocc:]/e_ai).sum()) - -################ - xref = solve(fvind, mo_energy, mo_occ, h1, s1*0, max_cycle=30)[0][:,mo_occ==0] - def fvind(x): - return numpy.dot(a,x.reshape(nd,nocc*nvir).T).T.reshape(nd,nvir,nocc) - h1 = h1[:,nocc:] - x0 = numpy.linalg.solve(numpy.diag(1/e_ai.ravel())+a, -h1.reshape(nd,-1).T).T.reshape(nd,nvir,nocc) - x1 = solve(fvind, mo_energy, mo_occ, h1, max_cycle=30)[0] - print(abs(x0-x1).sum()) - print(abs(xref-x1).sum()) - print(abs(h1 + fvind(x1)+x1/e_ai).sum()) + assert h1.ndim == 2 + return mo1[0], mo_e1[0] diff --git a/pyscf/scf/ucphf.py b/pyscf/scf/ucphf.py index e225dff496..e1149048d5 100644 --- a/pyscf/scf/ucphf.py +++ b/pyscf/scf/ucphf.py @@ -67,25 +67,26 @@ def solve_nos1(fvind, mo_energy, mo_occ, h1, mo1base = numpy.hstack((h1[0].reshape(-1,nvira*nocca), h1[1].reshape(-1,nvirb*noccb))) mo1base *= -e_ai + nov = e_ai.size def vind_vo(mo1): - mo1 = mo1.reshape(mo1base.shape) - v = fvind(mo1).reshape(mo1base.shape) + nd = mo1.shape[0] + v = fvind(mo1).reshape(nd, nov) if level_shift != 0: v -= mo1 * level_shift v *= e_ai - return v.ravel() - mo1 = lib.krylov(vind_vo, mo1base.ravel(), + return v.reshape(-1, nov) + mo1 = lib.krylov(vind_vo, mo1base.reshape(-1, nov), tol=tol, max_cycle=max_cycle, hermi=hermi, verbose=log) log.timer('krylov solver in CPHF', *t0) + mo1 = mo1.reshape(mo1base.shape) + mo1_a = mo1[:,:nvira*nocca].reshape(-1,nvira,nocca) + mo1_b = mo1[:,nvira*nocca:].reshape(-1,nvirb,noccb) if isinstance(h1[0], numpy.ndarray) and h1[0].ndim == 2: - mo1 = (mo1[:nocca*nvira].reshape(nvira,nocca), - mo1[nocca*nvira:].reshape(nvirb,noccb)) + mo1 = (mo1_a[0], mo1_b[0]) else: - mo1 = mo1.reshape(mo1base.shape) - mo1_a = mo1[:,:nvira*nocca].reshape(-1,nvira,nocca) - mo1_b = mo1[:,nvira*nocca:].reshape(-1,nvirb,noccb) + assert h1[0].ndim == 3 mo1 = (mo1_a, mo1_b) return mo1, None @@ -128,20 +129,22 @@ def solve_withs1(fvind, mo_energy, mo_occ, h1, s1, mo1base_a[:,occidxa] = -s1_a[:,occidxa] * .5 mo1base_b[:,occidxb] = -s1_b[:,occidxb] * .5 mo1base = numpy.hstack((mo1base_a.reshape(nset,-1), mo1base_b.reshape(nset,-1))) + nov = nocca * nmoa + noccb * nmob def vind_vo(mo1): - mo1 = mo1.reshape(mo1base.shape) - v = fvind(mo1).reshape(mo1base.shape) + nd = mo1.shape[0] + mo1 = mo1.reshape(nd, nov) + v = fvind(mo1).reshape(nd, nov) 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 = v[:,:nmoa*nocca].reshape(nd,nmoa,nocca) + v1b = v[:,nmoa*nocca:].reshape(nd,nmob,noccb) v1a[:,viridxa] *= eai_a v1b[:,viridxb] *= eai_b v1a[:,occidxa] = 0 v1b[:,occidxb] = 0 - return v.ravel() - mo1 = lib.krylov(vind_vo, mo1base.ravel(), + return v.reshape(nd, nov) + mo1 = lib.krylov(vind_vo, mo1base.reshape(-1, nov), tol=tol, max_cycle=max_cycle, hermi=hermi, verbose=log) mo1 = mo1.reshape(mo1base.shape) mo1_a = mo1[:,:nmoa*nocca].reshape(nset,nmoa,nocca) @@ -164,4 +167,6 @@ def vind_vo(mo1): if isinstance(h1[0], numpy.ndarray) and h1[0].ndim == 2: mo1_a, mo1_b = mo1_a[0], mo1_b[0] mo_e1_a, mo_e1_b = mo_e1_a[0], mo_e1_b[0] + else: + assert h1[0].ndim == 3 return (mo1_a, mo1_b), (mo_e1_a, mo_e1_b)