diff --git a/src/nest/deltascf/fr.py b/src/nest/deltascf/fr.py index 709e918..171f58c 100644 --- a/src/nest/deltascf/fr.py +++ b/src/nest/deltascf/fr.py @@ -22,6 +22,8 @@ finite projector shifts in the paper's unrestricted spin channel. """ +import weakref + import numpy as np from scipy.linalg import eigh from pyscf import lib, scf, symm @@ -30,6 +32,8 @@ def _set_imom(mf, mo_ref, setocc): + # Both occupation closures need access to mf, but must not own it. + mf = weakref.proxy(mf) scf.addons.mom_occ(mf, mo_ref, setocc) imom_occ = mf.get_occ if mf.mol.symmetry: @@ -40,9 +44,7 @@ def get_occ(mo_energy=None, mo_coeff=None): if not mf.mol.symmetry: occupations = imom_occ(mo_energy, mo_coeff) else: - # PySCF mom_occ ranks all orbitals together. Here retain the - # reference alpha/beta electron counts within each irrep, so - # overlap ranking cannot move electrons between symmetry blocks. + # Preserve the reference spin occupations within each irrep. if mo_coeff is None: mo_coeff = mf.mo_coeff orbital_symmetry = mf.get_orbsym(mo_coeff) @@ -68,25 +70,39 @@ def get_occ(mo_energy=None, mo_coeff=None): class FR(lib.StreamObject): - """Freeze both singly occupied orbitals, relax the others, then run IMOM. + """Freeze both SOMOs, then release all orbitals with IMOM (ROHF/ROKS). + + Requires spin=2 and real S-orthonormal input orbitals. With symmetry, + each orbital must have a definite irrep in an Abelian point group; + multiple orbitals may share an irrep. Dooh, Coov and SO3 are unsupported. + + Attributes: + freeze_tol : float + Frozen-stage orbital-gradient threshold. Default is 1e-5. + freeze_max_cycle : int + Maximum frozen-stage iterations. Default is 100. + Failure to converge raises before release starts. + conv_tol, conv_tol_grad : float + Energy and gradient thresholds for release; inherited from mf. + max_cycle : int + Maximum release iterations; inherited from mf. + callback : callable + Called after each iteration with the PySCF environment dict. + env['fr_stage'] is 'freeze' or 'release'. Default is None. - Use ``mf.FR().kernel(mo_coeff, mo_occ)`` on a plain ROHF/ROKS object. - Requires real S-orthonormal orbitals and spin=2. With symmetry enabled, - input orbitals must belong to individual irreps of an Abelian point - group (e.g. Cs, C2v, D2h). Both stages retain the initial alpha/beta - electron counts in each irrep. For linear molecules or atoms, select - an Abelian subgroup explicitly; Dooh, Coov and SO3 are not supported. - Both stages use ordinary SCF and fresh CDIIS histories, with no Hessian. - Energies always come from the physical Hamiltonian. + Saved results: + freeze_converged : bool + Whether the frozen stage converged. + freeze_cycles : int + Number of frozen-stage iterations. + converged, cycles + Convergence status and iteration count of the release stage. + e_tot, mo_energy, mo_coeff, mo_occ + Final energy, orbital energies, coefficients and occupations. - freeze_tol (default 1e-5) bounds PySCF's projected orbital-gradient norm; - freeze_max_cycle (100) bounds the frozen stage. Failure raises before - release. freeze_cycles/freeze_converged describe that stage; cycles and - converged describe release, governed by the usual SCF tolerances. - callback receives both stages, identified by env['fr_stage']. - PySCF's extra relaxed convergence check is disabled in both stages. - At INFO verbosity (verbose=4), stage starts and the final SOMO-subspace - overlap singular values relative to the input orbitals are logged. + Examples: + >>> excited = mf.FR() + >>> excited.kernel(mo_coeff, mo_occ) """ __name_mixin__ = 'FR' @@ -107,10 +123,35 @@ def __init__(self, mf): if mf is not self: self.__dict__.update(mf.__dict__) + def check_sanity(self): + if self.mol.spin != 2: + raise ValueError('FR requires spin=2') + if self.mol.symmetry and self.mol.groupname in ('Dooh', 'Coov', 'SO3'): + raise NotImplementedError('FR supports Abelian point groups; select an Abelian subgroup') + if not np.isfinite(self.freeze_tol) or self.freeze_tol <= 0: + raise ValueError('freeze_tol must be finite and positive') + if not isinstance(self.freeze_max_cycle, (int, np.integer)) or self.freeze_max_cycle < 1: + raise ValueError('freeze_max_cycle must be a positive integer') + super().check_sanity() + return self + + def dump_flags(self, verbose=None): + super().dump_flags(verbose) + log = logger.new_logger(self, verbose) + log.info('FR freeze gradient tolerance = %g', self.freeze_tol) + log.info('FR freeze max_cycle = %d', self.freeze_max_cycle) + log.info('FR release occupation method = IMOM') + return self + def kernel(self, mo_coeff=None, mo_occ=None): + self.check_sanity() + if self.verbose >= logger.INFO: + self.dump_flags() + log = logger.new_logger(self) + c = np.array(self.mo_coeff if mo_coeff is None else mo_coeff, copy=True) occ = np.array(self.mo_occ if mo_occ is None else mo_occ, copy=True) - if (self.mol.spin != 2 or occ.ndim != 1 or + if (occ.ndim != 1 or not np.all(np.isin(occ, [0, 1, 2])) or np.count_nonzero(occ == 1) != 2 or occ.sum() != self.mol.nelectron): raise ValueError('FR requires a spin=2 occupation vector with exactly two SOMOs') @@ -121,8 +162,6 @@ def kernel(self, mo_coeff=None, mo_occ=None): raise ValueError('FR requires S-orthonormal input orbitals') orbital_symmetry = np.zeros(occ.size, dtype=int) if self.mol.symmetry: - if self.mol.groupname not in ('C1', 'Ci', 'Cs', 'C2', 'C2v', 'C2h', 'D2', 'D2h'): - raise NotImplementedError('FR supports Abelian point groups; select an Abelian subgroup') orbital_symmetry = symm.label_orb_symm( self.mol, self.mol.irrep_id, self.mol.symm_orb, c, s=s, check=True) c = lib.tag_array(c, orbsym=orbital_symmetry) @@ -135,8 +174,6 @@ def kernel(self, mo_coeff=None, mo_occ=None): else tuple(requested) == (alpha, beta)) if not matches: raise ValueError('Input occupations conflict with irrep_nelec for ' + name) - if self.freeze_tol <= 0 or self.freeze_max_cycle < 1: - raise ValueError('freeze_tol and freeze_max_cycle must be positive') self.converged = False self.freeze_converged = False @@ -156,8 +193,7 @@ def kernel(self, mo_coeff=None, mo_occ=None): mf.diis.rollback = self.diis_space_rollback mf.diis.damp = self.diis_damp else: - # Let the ordinary kernel build its orthogonalized DIIS - # residual, without carrying over the frozen history. + # Start release with a fresh DIIS history. mf.diis = True mf.DIIS = scf.diis.CDIIS def stage_callback(env, stage=stage): @@ -165,23 +201,18 @@ def stage_callback(env, stage=stage): callback(dict(env, fr_stage=stage)) mf.callback = stage_callback - # The two singly occupied columns never change during the freeze - # iterations. All other orbitals are linear combinations of the - # INITIAL doubly occupied and virtual columns, an S-orthonormal basis - # for the allowed space. Put occupied indices first: eigh returns - # increasing energies, so the lowest eigenvectors fill these slots. + # Exclude SOMOs; place doubly occupied slots before virtual slots. doubly_occupied_indices = np.flatnonzero(occ == 2) virtual_indices = np.flatnonzero(occ == 0) free_indices = np.concatenate((doubly_occupied_indices, virtual_indices)) free_orbitals = c[:, free_indices] free_blocks = [free_indices[orbital_symmetry[free_indices] == irrep] for irrep in np.unique(orbital_symmetry[free_indices])] - original_get_grad = frozen.get_grad + # Bind to the unmodified base, not the object that owns this closure. + original_get_grad = base.get_grad def frozen_eig(fock, overlap, **kwargs): - # In each free block B, B.T @ S @ B = I. Thus this is an ordinary - # eigenproblem, not an AO generalized eigenproblem. Separate - # irreps must not mix, even if their eigenvalues are degenerate. + # B.T @ S @ B = I; diagonalize each irrep separately. coeff = c.copy() energies = np.einsum('pi,pi->i', c, fock @ c) for indices in free_blocks: @@ -194,16 +225,12 @@ def frozen_eig(fock, overlap, **kwargs): return energies, coeff def frozen_get_occ(mo_energy=None, mo_coeff=None): - # frozen_eig already placed the lowest free eigenvectors into - # the doubly occupied slots; keep the initial occupation labels. + # frozen_eig already orders the free orbitals by energy. return occ.copy() def frozen_get_grad(coeff, occupations, fock=None): - # PySCF packs all nonredundant ROHF rotations into a 1-D vector: - # alpha: occupied -> virtual; beta: doubly -> singly/virtual. - # Select only doubly occupied -> virtual entries of THAT vector. - # Derive masks from the supplied occupations because symmetry - # SCF can reorder orbitals during its finalization step. + # Select doubly occupied -> virtual entries in PySCF's packed gradient. + # Rebuild the mask because symmetry finalization can reorder orbitals. alpha_rotations = (occupations == 0)[:, None] & (occupations > 0)[None, :] beta_rotations = (occupations != 2)[:, None] & (occupations == 2)[None, :] all_rotations = alpha_rotations | beta_rotations @@ -212,25 +239,17 @@ def frozen_get_grad(coeff, occupations, fock=None): gradient = original_get_grad(coeff, occupations, fock) return gradient[free_gradient_entries] - # These three hooks implement the constrained eigenproblem, filling, - # and convergence test. Density/Fock/energy updates stay in PySCF. frozen.eig = frozen_eig frozen.get_occ = frozen_get_occ frozen.get_grad = frozen_get_grad if frozen.diis: - # CDIIS must ignore rotations involving fixed singly occupied - # orbitals too. For the physical ROHF Fock, its virtual/doubly - # occupied block is (F_alpha + F_beta)/2. Multiplying by the - # occupation difference 2 gives exactly PySCF's orbital gradient: - # ||B.T @ (F D S - S D F) @ B||_F = sqrt(2) * ||g_free||_2, - # where D is the total (alpha + beta) density. The same-irrep - # restriction is applied by CDIIS using the orbsym tag below. + # Project the DIIS residual onto the free orbital space. if self.mol.symmetry: free_orbitals = lib.tag_array(free_orbitals, orbsym=orbital_symmetry[free_indices]) frozen.diis.Corth = free_orbitals frozen.conv_tol_grad = self.freeze_tol frozen.max_cycle = self.freeze_max_cycle - logger.info(self, 'FR Freeze stage begins: fix both singly occupied orbitals; ' + log.info('FR Freeze stage begins: fix both singly occupied orbitals; ' 'relax doubly occupied/virtual orbitals') frozen.kernel(dm0=frozen.make_rdm1(c, occ)) self.freeze_cycles = frozen.cycles @@ -238,26 +257,18 @@ def frozen_get_grad(coeff, occupations, fock=None): if not self.freeze_converged: raise RuntimeError('FR frozen SCF did not converge; release was not started') - # PySCF mom_occ captures these occupied reference spaces once (IMOM). - # Resetting the reference after freeze incorporates spectator relaxation. - logger.info(self, 'FR Release stage begins: relax all orbitals with the post-freeze IMOM reference') - # Symmetry SCF finalization sorts columns AND occupations. Use both - # returned arrays together, rather than reusing initial column labels. + log.info('FR Release stage begins: relax all orbitals with the post-freeze IMOM reference') + # Use the post-freeze reference, including symmetry reordering. setocc = np.asarray((frozen.mo_occ > 0, frozen.mo_occ == 2), dtype=float) _set_imom(release, frozen.mo_coeff, setocc) release.kernel(dm0=frozen.make_rdm1()) - # With conv_check=False the ordinary kernel already checks energy - # AND orbital-gradient convergence. Honor its result, including any - # user-supplied check_convergence callback, without a second criterion. - # Retain the ordinary MF results and caches, but not stage-local hooks. + # Copy results without retaining stage-local hooks. for key, value in release.__dict__.items(): if key not in ('callback', 'get_occ', 'diis', 'DIIS'): self.__dict__[key] = value _set_imom(self, frozen.mo_coeff, setocc) - # Use the final occupation labels: IMOM and symmetry finalization - # can move the singly occupied orbitals to different column indices. somo_overlap = c[:, occ == 1].T @ s @ self.mo_coeff[:, self.mo_occ == 1] - logger.info(self, 'FR target SOMO overlap singular values: %s', + log.note('FR target SOMO overlap singular values: %s', np.linalg.svd(somo_overlap, compute_uv=False)) return self.e_tot diff --git a/src/nest/deltascf/sgm.py b/src/nest/deltascf/sgm.py index d88ee2b..0d5d505 100644 --- a/src/nest/deltascf/sgm.py +++ b/src/nest/deltascf/sgm.py @@ -1,5 +1,6 @@ #!/usr/bin/env python # Copyright 2026 The NEST Developers. All Rights Reserved. +# Copyright 2014-2021 The PySCF Developers. All Rights Reserved. # # Licensed under the Apache License, Version 2.0 (the "License"); # you may not use this file except in compliance with the License. @@ -23,7 +24,7 @@ from pyscf import __config__, lib from pyscf.lib import logger -from pyscf.scf import hf, hf_symm +from pyscf.scf import addons, hf, hf_symm _ARMIJO_C1 = 1e-4 @@ -33,84 +34,90 @@ _MAX_PRECOND = 1e4 +# Adapted from PySCF PR #3448, commit e45d0fe9c201825dd6b7069ade7d6532b9952617. +# https://github.com/pyscf/pyscf/pull/3448 def gen_g_hop_rohf(mf, mo_coeff, mo_occ, fock_ao=None, h1e=None, with_symmetry=True): - '''ROHF orbital gradient and full-K Hessian-vector product.''' - mol = mf.mol + '''ROHF orbital gradient and symmetric orbital Hessian action.''' mo_occ = numpy.asarray(mo_occ) - if h1e is None: - h1e = mf.get_hcore(mol) - if fock_ao is None or getattr(fock_ao, 'focka', None) is None: + mol = mf.mol + mo_coeff0 = mo_coeff + if getattr(mf, '_scf', None) and mf._scf.mol != mol: + # TODO: construct vind with dual-basis treatment + mo_coeff = addons.project_mo_nr2nr(mf._scf.mol, mo_coeff, mol) + + if getattr(fock_ao, 'focka', None) is None: + if getattr(mf, '_scf', None) and mf._scf.mol != mol: + h1e = mf.get_hcore(mol) dm0 = mf.make_rdm1(mo_coeff, mo_occ) - vhf = mf.get_veff(mol, dm0) - focka_ao = h1e + vhf[0] - fockb_ao = h1e + vhf[1] + fock_ao = mf.get_fock(h1e, dm=dm0) + focka = reduce(numpy.dot, (mo_coeff.conj().T, fock_ao.focka, mo_coeff)) + fockb = reduce(numpy.dot, (mo_coeff.conj().T, fock_ao.fockb, mo_coeff)) else: - focka_ao, fockb_ao = fock_ao.focka, fock_ao.fockb - - focka = reduce(numpy.dot, (mo_coeff.conj().T, focka_ao, mo_coeff)) - fockb = reduce(numpy.dot, (mo_coeff.conj().T, fockb_ao, mo_coeff)) - occidxa = mo_occ > 0 - occidxb = mo_occ == 2 + focka = reduce(numpy.dot, (mo_coeff0.conj().T, fock_ao.focka, mo_coeff0)) + fockb = reduce(numpy.dot, (mo_coeff0.conj().T, fock_ao.fockb, mo_coeff0)) + mo_occa = occidxa = mo_occ > 0 + mo_occb = occidxb = mo_occ == 2 viridxa = ~occidxa viridxb = ~occidxb - uniq_var_a = viridxa[:, None] & occidxa - uniq_var_b = viridxb[:, None] & occidxb - uniq_var = uniq_var_a | uniq_var_b - - gmat = numpy.zeros_like(focka) - gmat[uniq_var_a] = focka[uniq_var_a] - gmat[uniq_var_b] += fockb[uniq_var_b] - g = gmat[uniq_var] - - focka_diag = focka.diagonal().real - fockb_diag = fockb.diagonal().real - h_diag_mat = numpy.zeros_like(focka_diag[:, None] - focka_diag) - h_diag_mat[uniq_var_a] = ( - focka_diag[:, None] - focka_diag)[uniq_var_a] - h_diag_mat[uniq_var_b] += ( - fockb_diag[:, None] - fockb_diag)[uniq_var_b] - h_diag = h_diag_mat[uniq_var] + uniq_var_a = viridxa[:,None] & occidxa + uniq_var_b = viridxb[:,None] & occidxb + uniq_ab = uniq_var_a | uniq_var_b + nmo = mo_coeff.shape[-1] + orboa = mo_coeff[:,occidxa] + orbob = mo_coeff[:,occidxb] + orbva = mo_coeff[:,viridxa] + orbvb = mo_coeff[:,viridxb] + + g = numpy.zeros_like(focka) + g[uniq_var_a] = focka[uniq_var_a] + g[uniq_var_b] += fockb[uniq_var_b] + g = g[uniq_ab] + ea = focka.diagonal().real + eb = fockb.diagonal().real + h_diag = numpy.zeros_like(focka.real) + h_diag[uniq_var_a] = (ea[:,None] - ea)[uniq_var_a] + h_diag[uniq_var_b] += (eb[:,None] - eb)[uniq_var_b] + h_diag = h_diag[uniq_ab] if with_symmetry and mol.symmetry: orbsym = hf_symm.get_orbsym(mol, mo_coeff) - sym_forbid = (orbsym[:, None] != orbsym)[uniq_var] - g = g.copy() - h_diag = h_diag.copy() + sym_forbid = (orbsym[:,None] != orbsym)[uniq_ab] g[sym_forbid] = 0 h_diag[sym_forbid] = 0 - mo_occ_a = occidxa.astype(numpy.double) - mo_occ_b = occidxb.astype(numpy.double) - vind = mf.gen_response((mo_coeff, mo_coeff), (mo_occ_a, mo_occ_b), + gmat = hf.unpack_uniq_var(g, mo_occ) + vind = mf.gen_response((mo_coeff,)*2, (mo_occa, mo_occb), hermi=1, with_nlc=False) def h_op(x): if with_symmetry and mol.symmetry: x = x.copy() x[sym_forbid] = 0 - kappa = hf.unpack_uniq_var(x, mo_occ) - - dm1a_mo = kappa * (mo_occ_a[None, :] - mo_occ_a[:, None]) - dm1b_mo = kappa * (mo_occ_b[None, :] - mo_occ_b[:, None]) - dm1 = numpy.asarray(( - reduce(numpy.dot, (mo_coeff, dm1a_mo, mo_coeff.conj().T)), - reduce(numpy.dot, (mo_coeff, dm1b_mo, mo_coeff.conj().T)))) + x1 = numpy.zeros((nmo,nmo), dtype=x.dtype) + x1[uniq_ab] = x + x1a = x1[uniq_var_a].reshape(orbva.shape[1], orboa.shape[1]) + x1b = x1[uniq_var_b].reshape(orbvb.shape[1], orbob.shape[1]) + d1a = reduce(numpy.dot, (orbva, x1a, orboa.conj().T)) + d1b = reduce(numpy.dot, (orbvb, x1b, orbob.conj().T)) + dm1 = numpy.array((d1a+d1a.conj().T, d1b+d1b.conj().T)) v1a, v1b = vind(dm1) + # Keep the shared rotation's oo/vv contributions for each spin. + kappa = hf.unpack_uniq_var(x, mo_occ) hmat_a = focka.dot(kappa) - kappa.dot(focka) hmat_b = fockb.dot(kappa) - kappa.dot(fockb) hmat_a += reduce(numpy.dot, (mo_coeff.conj().T, v1a, mo_coeff)) hmat_b += reduce(numpy.dot, (mo_coeff.conj().T, v1b, mo_coeff)) - - out = numpy.zeros_like(focka) - out[uniq_var_a] = hmat_a[uniq_var_a] - out[uniq_var_b] += hmat_b[uniq_var_b] - out = out[uniq_var] + hx = numpy.zeros_like(hmat_a) + hx[uniq_var_a] = hmat_a[uniq_var_a] + hx[uniq_var_b] += hmat_b[uniq_var_b] + # Transport the moving-frame gradient back to the reference frame. + hx += .5 * (kappa.dot(gmat) - gmat.dot(kappa)) + hx = hx[uniq_ab] if with_symmetry and mol.symmetry: - out = out.copy() - out[sym_forbid] = 0 - return out + hx[sym_forbid] = 0 + return hx return g, h_op, h_diag @@ -222,8 +229,7 @@ class SGM(lib.StreamObject): def __new__(cls, mf): if isinstance(mf, SGM): return mf - assert isinstance(mf, hf.SCF) - if not mf.istype('ROHF'): + if not isinstance(mf, hf.SCF) or not mf.istype('ROHF'): raise NotImplementedError('SGM currently supports ROHF/ROKS objects') obj = object.__new__(cls) return lib.set_class(obj, (cls, mf.__class__)) @@ -234,6 +240,15 @@ def __init__(self, mf): self.__dict__.update(mf.__dict__) self._scf = mf + def check_sanity(self): + if not numpy.isfinite(self.gradient_scale) or self.gradient_scale <= 0: + raise ValueError('gradient_scale must be finite and positive') + tol = self.tol if self.conv_tol_grad is None else self.conv_tol_grad + if not numpy.isfinite(tol) or tol <= 0: + raise ValueError('SGM gradient tolerance must be finite and positive') + super().check_sanity() + return self + def dump_flags(self, verbose=None): super().dump_flags(verbose) log = logger.new_logger(self, verbose) @@ -307,7 +322,10 @@ def _descent_direction(self, history, grad_delta, opt_grad_delta, h0_inv, return direction def kernel(self, mo_coeff=None, mo_occ=None): - log = logger.new_logger(self, self.verbose) + self.check_sanity() + if self.verbose >= logger.INFO: + self.dump_flags() + log = logger.new_logger(self) if mo_coeff is None: mo_coeff = self.mo_coeff if mo_occ is None: @@ -317,9 +335,6 @@ def kernel(self, mo_coeff=None, mo_occ=None): if mo_coeff is None: raise RuntimeError('mo_coeff must be specified for SGM') tol = self.tol if self.conv_tol_grad is None else self.conv_tol_grad - if not numpy.isfinite(self.gradient_scale) or self.gradient_scale <= 0: - raise ValueError('gradient_scale must be finite and positive') - self.dump_flags() mo_guess = mo_coeff.copy() h1e = self.get_hcore() diff --git a/src/nest/deltascf/tests/test_fr.py b/src/nest/deltascf/tests/test_fr.py index 29b98e3..bf6f98f 100644 --- a/src/nest/deltascf/tests/test_fr.py +++ b/src/nest/deltascf/tests/test_fr.py @@ -14,174 +14,103 @@ # limitations under the License. import unittest -from unittest import mock - -import numpy as np -from scipy.linalg import expm -from pyscf import gto, scf, symm +from pyscf import gto from nest.deltascf import fr class KnownValues(unittest.TestCase): @classmethod def setUpClass(cls): - cls.mol = gto.M(atom='O 0 0 0; H 0 -.757 .587; H 0 .757 .587', - basis='6-31g', spin=2, verbose=0) - - def test_freeze_release(self): - for xc, use_symmetry in ((None, False), ('BHandHLYP', False), - (None, True), ('BHandHLYP', True)): - with self.subTest(xc=xc, symmetry=use_symmetry): - mol = gto.M(atom=self.mol.atom, basis=self.mol.basis, spin=2, - symmetry=use_symmetry, verbose=0) - ground_mol = mol.copy() - ground_mol.spin = 0 - ground = ground_mol.RHF() if xc is None else ground_mol.RKS(xc=xc) - mf = mol.ROHF() if xc is None else mol.ROKS(xc=xc) - if xc is not None: - ground.grids.level = mf.grids.level = 0 - ground.kernel() - occ = ground.mo_occ.copy() - occ[[4, 5]] = 1 - c = ground.mo_coeff.copy() - s = mf.get_ovlp() - stages = [] - - def callback(env): - stages.append(env['fr_stage']) - coeff = env['mo_coeff'] - np.testing.assert_allclose(coeff.T @ s @ coeff, np.eye(occ.size), atol=1e-10) - if env['fr_stage'] == 'freeze': - np.testing.assert_array_equal(coeff[:, occ == 1], c[:, occ == 1]) - free_basis = env['mf_diis'].Corth - np.testing.assert_allclose(free_basis.T @ s @ c[:, occ == 1], 0, atol=1e-10) - # Independently verify that projected CDIIS residual - # and the physical frozen gradient have the same zero. - dm_total = env['dm'].sum(axis=0) - fock = env['fock'] # Physical Fock at CURRENT density. - residual = free_basis.T @ (fock @ dm_total @ s - s @ dm_total @ fock) @ free_basis - full_gradient = mf.get_grad(coeff, occ, fock) - all_rotations = ((occ == 0)[:, None] & (occ > 0)) | ((occ != 2)[:, None] & (occ == 2)) - free_entries = ((occ == 0)[:, None] & (occ == 2))[all_rotations] - frozen_gradient = full_gradient[free_entries] - self.assertAlmostEqual(np.linalg.norm(residual), - np.sqrt(2) * np.linalg.norm(frozen_gradient), places=10) - np.testing.assert_allclose(env['mf'].get_grad(coeff, occ, fock), frozen_gradient, atol=1e-12) - if stages.count('freeze') == 1: - # Energy finite differences are independent of - # the gradient packing and the DIIS construction. - direction = np.zeros_like(full_gradient) - direction[free_entries] = frozen_gradient - direction /= np.linalg.norm(direction) - rotation = scf.hf.unpack_uniq_var(direction, occ) - eps = 1e-4 - plus = mf.energy_tot(dm=mf.make_rdm1(coeff @ expm(eps * rotation), occ)) - minus = mf.energy_tot(dm=mf.make_rdm1(coeff @ expm(-eps * rotation), occ)) - self.assertAlmostEqual((plus-minus)/(2*eps), 2 * full_gradient @ direction, places=6) - else: - self.assertEqual(np.count_nonzero(env['mo_occ'] == 1), 2) - - opt = mf.FR().set(conv_tol=1e-10, conv_tol_grad=1e-6, max_cycle=100, - callback=callback) - self.assertIsInstance(opt, fr.FR) - self.assertTrue(opt.istype('ROHF')) - # First-order implementation must not require response kernels. - with mock.patch.object(type(mf), 'gen_response', side_effect=AssertionError('Hessian used')): - opt.kernel(c, occ) - self.assertTrue(opt.freeze_converged) - self.assertTrue(opt.converged) - self.assertIn('freeze', stages) - self.assertIn('release', stages) - self.assertLess(np.linalg.norm(opt.get_grad(opt.mo_coeff, opt.mo_occ)), 1e-6) - self.assertAlmostEqual(opt.e_tot, opt.energy_tot(dm=opt.make_rdm1()), places=10) - np.testing.assert_array_equal(opt.get_occ(opt.mo_energy, opt.mo_coeff), opt.mo_occ) - # Independent conventional triplet SCF reference for this valence state. - ref = mf.copy().set(conv_tol=1e-11, conv_tol_grad=1e-7) - ref.kernel() - self.assertTrue(ref.converged) - self.assertAlmostEqual(opt.e_tot, ref.e_tot, places=7) - self.assertIsNone(mf.mo_coeff) - self.assertIs(opt.callback, callback) - self.assertIs(opt.FR(), opt) - if use_symmetry: - initial_symmetry = symm.label_orb_symm(mol, mol.irrep_id, mol.symm_orb, c) - final_symmetry = symm.label_orb_symm(mol, mol.irrep_id, mol.symm_orb, opt.mo_coeff) - np.testing.assert_array_equal(opt.mo_coeff.orbsym, final_symmetry) - for irrep in mol.irrep_id: - for occupied in (lambda o: o > 0, lambda o: o == 2): - self.assertEqual(np.count_nonzero(occupied(occ) & (initial_symmetry == irrep)), - np.count_nonzero(occupied(opt.mo_occ) & (final_symmetry == irrep))) - - def test_failed_freeze_does_not_release(self): - ground_mol = self.mol.copy() - ground_mol.spin = 0 - ground = ground_mol.RHF().run() - occ = ground.mo_occ.copy() - occ[[0, 5]] = 1 - stages = [] - opt = fr.FR(self.mol.ROHF()).set(freeze_max_cycle=1, freeze_tol=1e-12, - callback=lambda env: stages.append(env['fr_stage'])) - with self.assertRaisesRegex(RuntimeError, 'release was not started'): - opt.kernel(ground.mo_coeff, occ) - self.assertFalse(opt.converged) - self.assertFalse(opt.freeze_converged) - self.assertEqual(stages, ['freeze']) - - def test_reject_invalid_reference(self): - mf = self.mol.ROHF().run() - with self.assertRaisesRegex(ValueError, 'S-orthonormal'): - mf.FR().kernel(mf.mo_coeff * 2, mf.mo_occ) - with self.assertRaisesRegex(ValueError, 'two SOMOs'): - mf.FR().kernel(mf.mo_coeff, np.zeros_like(mf.mo_occ)) - with self.assertRaises(NotImplementedError): - fr.FR(self.mol.UHF()) - - def test_symmetry_finalization_reorders_columns(self): - mol = gto.M(atom=self.mol.atom, basis=self.mol.basis, spin=2, symmetry=True, verbose=0) - mf = mol.ROHF().set(conv_tol=1e-11, conv_tol_grad=1e-8).run() - order = np.arange(mf.mo_occ.size)[::-1] - # In particular, singly occupied columns are no longer immediately - # after doubly occupied ones. PySCF finalization will reorder them. - coeff = np.asarray(mf.mo_coeff[:, order]) - occupations = mf.mo_occ[order] - opt = mf.FR().set(conv_tol=1e-10, conv_tol_grad=1e-7) - opt.kernel(coeff, occupations) + atom = ''' + O 0.000 0.000 0.000 + H 0.000 -0.757 0.587 + H 0.000 0.757 0.587 + ''' + cls.mol0 = gto.M(atom=atom, basis='6-31g', spin=0, verbose=0) + cls.mol1 = gto.M(atom=atom, basis='6-31g', spin=2, verbose=0) + + def test_hcho_rydberg_triplet(self): + # Same Rydberg excitation as test_sgm.py; Q-Chem SOMOs 8 and 15. + # Oracle: vsQChem/HCHO/QChem/hcho.log, triplet energy -113.131417 Ha. + atom = ''' + C 0.00000000 0.00000000 -1.13947666 + O 0.00000000 0.00000000 1.14402883 + H 0.00000000 1.76627623 -2.23398653 + H 0.00000000 -1.76627623 -2.23398653 + ''' + mol = gto.M(atom=atom, basis='cc-pvdz', spin=0, verbose=0) + ground = mol.RKS(xc='BHandHLYP') + ground.grids.atom_grid = (75, 302) + ground.kernel() + occupations = ground.mo_occ.copy() + occupations[7] = 1 + occupations[14] = 1 + + triplet = mol.copy() + triplet.spin = 2 + mf = triplet.ROKS(xc='BHandHLYP') + mf.grids.atom_grid = (75, 302) + opt = fr.FR(mf) + opt.max_cycle = 100 + e_tot = opt.kernel(ground.mo_coeff, occupations) + + self.assertTrue(opt.freeze_converged) + self.assertTrue(opt.converged) + self.assertAlmostEqual(e_tot, -113.131417, delta=1e-4) + # PySCF SGM reference for the same state and numerical settings. + self.assertAlmostEqual(e_tot, -113.131420102713, delta=1e-6) + + def test_h2o_triplet_rohf(self): + ground = self.mol0.RHF().run() + occupations = ground.mo_occ.copy() + occupations[4] = 1 # HOMO -> LUMO: 2/0 becomes 1/1. + occupations[5] = 1 + + opt = fr.FR(self.mol1.ROHF()) + opt.conv_tol = 1e-10 + opt.conv_tol_grad = 1e-6 + opt.max_cycle = 100 + e_tot = opt.kernel(ground.mo_coeff, occupations) + + # This valence excitation reaches the ordinary triplet SCF solution. + reference = self.mol1.ROHF() + reference.conv_tol = 1e-11 + reference.conv_tol_grad = 1e-7 + reference.kernel() + + self.assertTrue(opt.freeze_converged) + self.assertTrue(opt.converged) + self.assertTrue(reference.converged) + self.assertAlmostEqual(e_tot, reference.e_tot, places=7) + + def test_h2o_triplet_roks(self): + ground = self.mol0.RKS(xc='BHandHLYP') + ground.grids.level = 0 + ground.kernel() + occupations = ground.mo_occ.copy() + occupations[4] = 1 # HOMO -> LUMO: 2/0 becomes 1/1. + occupations[5] = 1 + + mf = self.mol1.ROKS(xc='BHandHLYP') + mf.grids.level = 0 + opt = fr.FR(mf) + opt.conv_tol = 1e-10 + opt.conv_tol_grad = 1e-6 + opt.max_cycle = 100 + e_tot = opt.kernel(ground.mo_coeff, occupations) + + reference = self.mol1.ROKS(xc='BHandHLYP') + reference.grids.level = 0 + reference.conv_tol = 1e-11 + reference.conv_tol_grad = 1e-7 + reference.kernel() + + self.assertTrue(opt.freeze_converged) self.assertTrue(opt.converged) - self.assertAlmostEqual(opt.e_tot, mf.e_tot, places=9) - np.testing.assert_allclose(opt.make_rdm1(), mf.make_rdm1(), atol=1e-6) - np.testing.assert_array_equal(opt.get_occ(opt.mo_energy, opt.mo_coeff), opt.mo_occ) - - def test_symmetry_rejects_mixed_orbitals_and_conflicting_counts(self): - mol = gto.M(atom=self.mol.atom, basis=self.mol.basis, spin=2, symmetry=True, verbose=0) - mf = mol.ROHF().run() - labels = mf.mo_coeff.orbsym - first = 0 - second = np.flatnonzero(labels != labels[first])[0] - mixed = mf.mo_coeff.copy() - mixed[:, first] = (mf.mo_coeff[:, first] + mf.mo_coeff[:, second]) / np.sqrt(2) - mixed[:, second] = (mf.mo_coeff[:, first] - mf.mo_coeff[:, second]) / np.sqrt(2) - with self.assertRaisesRegex(ValueError, 'not symmetrized'): - mf.FR().kernel(mixed, mf.mo_occ) - name = mol.irrep_name[mol.irrep_id.index(int(labels[first]))] - mf.irrep_nelec = {name: (0, 0)} - with self.assertRaisesRegex(ValueError, 'conflict'): - mf.FR().kernel(mf.mo_coeff, mf.mo_occ) - - def test_imom_rejects_nonnested_spaces(self): - mf = self.mol.ROHF().run() - setocc = np.asarray((mf.mo_occ > 0, mf.mo_occ == 2), dtype=float) - - def invalid_imom(mf, coeff, occupations): - invalid = np.zeros_like(mf.mo_occ) - invalid[:3] = 2 - invalid[3:7] = 1 # Same electron count, but four SOMOs. - mf.get_occ = lambda *args: invalid - - with mock.patch.object(scf.addons, 'mom_occ', side_effect=invalid_imom): - fr._set_imom(mf, mf.mo_coeff, setocc) - with self.assertRaisesRegex(RuntimeError, 'nonnested'): - mf.get_occ(mf.mo_energy, mf.mo_coeff) + self.assertTrue(reference.converged) + self.assertAlmostEqual(e_tot, reference.e_tot, places=7) if __name__ == '__main__': + print('Full tests for deltascf.fr') unittest.main() diff --git a/src/nest/deltascf/tests/test_sgm.py b/src/nest/deltascf/tests/test_sgm.py index 386f183..baa05ea 100644 --- a/src/nest/deltascf/tests/test_sgm.py +++ b/src/nest/deltascf/tests/test_sgm.py @@ -15,7 +15,6 @@ import unittest import numpy - from pyscf import gto from nest.deltascf import sgm