diff --git a/pyscf/cc/gccsd.py b/pyscf/cc/gccsd.py index d7e842ffa1..2db369cfd5 100644 --- a/pyscf/cc/gccsd.py +++ b/pyscf/cc/gccsd.py @@ -362,6 +362,7 @@ def _make_eris_incore(mycc, mo_coeff=None, ao2mofn=None): mo = mo_a + mo_b eri = ao2mo.kernel(mycc._scf._eri, mo) if eri.size == nmo**4: # if mycc._scf._eri is a complex array + eri = eri.reshape(nmo**2, nmo**2) sym_forbid = (orbspin[:,None] != orbspin).ravel() else: # 4-fold symmetry sym_forbid = (orbspin[:,None] != orbspin)[np.tril_indices(nmo)] diff --git a/pyscf/dft/rks.py b/pyscf/dft/rks.py index 70f5fdea37..807e32a0f5 100644 --- a/pyscf/dft/rks.py +++ b/pyscf/dft/rks.py @@ -539,5 +539,4 @@ def to_gpu(self): obj = lib.to_gpu(hf.SCF.reset(self.view(RKS))) # Attributes only defined in gpu4pyscf.RKS obj.screen_tol = 1e-14 - obj.disp = None return obj diff --git a/pyscf/dft/uks.py b/pyscf/dft/uks.py index 19491012a7..cf2f7d7486 100644 --- a/pyscf/dft/uks.py +++ b/pyscf/dft/uks.py @@ -203,5 +203,4 @@ def to_gpu(self): obj = lib.to_gpu(SCF.reset(self.view(UKS))) # Attributes only defined in gpu4pyscf.RKS obj.screen_tol = 1e-14 - obj.disp = None return obj diff --git a/pyscf/grad/mp2.py b/pyscf/grad/mp2.py index 58dcc92cf5..65ff7c8834 100644 --- a/pyscf/grad/mp2.py +++ b/pyscf/grad/mp2.py @@ -241,12 +241,16 @@ def converged(self): def _shell_prange(mol, start, stop, blksize): + l = mol._bas[start:stop,gto.ANG_OF] + if mol.cart: + dims = (l+1)*(l+2)//2 * mol._bas[start:stop,gto.NCTR_OF] + else: + dims = (l*2+1) * mol._bas[start:stop,gto.NCTR_OF] nao = 0 ib0 = start - for ib in range(start, stop): - now = (mol.bas_angular(ib)*2+1) * mol.bas_nctr(ib) + for ib, now in zip(range(start, stop), dims): nao += now - if nao > blksize and nao > now: + if nao > blksize: yield (ib0, ib, nao-now) ib0 = ib nao = now diff --git a/pyscf/grad/test/test_mp2.py b/pyscf/grad/test/test_mp2.py index 1d920cb772..24316ea50c 100644 --- a/pyscf/grad/test/test_mp2.py +++ b/pyscf/grad/test/test_mp2.py @@ -153,6 +153,14 @@ def test_symmetrize(self): g = mol.RHF.run().MP2().run().Gradients().kernel() self.assertAlmostEqual(lib.fp(g), 0.049987975650731625, 6) + # issue 1985 + def test_cart_gto(self): + mol1 = mol.copy() + mol1.cart = True + mol1.basis = '6-31g*' + g = mol.RHF.run().MP2().run().Gradients().kernel() + self.assertAlmostEqual(lib.fp(g), -0.03568120792884476, 6) + if __name__ == "__main__": print("Tests for MP2 gradients") diff --git a/pyscf/gto/mole.py b/pyscf/gto/mole.py index 6e642f66a7..e330758b3c 100644 --- a/pyscf/gto/mole.py +++ b/pyscf/gto/mole.py @@ -419,7 +419,7 @@ def str2atm(line): return list(zip(z, c.tolist())) #TODO: sort exponents -def format_basis(basis_tab): +def format_basis(basis_tab, sort_basis=True): '''Convert the input :attr:`Mole.basis` to the internal data format. ``{ atom: [(l, ((-exp, c_1, c_2, ..), @@ -457,10 +457,16 @@ def format_basis(basis_tab): fmt_basis = {} for atom, atom_basis in basis_tab.items(): symb = _atom_symbol(atom) - fmt_basis[symb] = basis_converter(symb, atom_basis) - - if len(fmt_basis[symb]) == 0: + _basis = basis_converter(symb, atom_basis) + if len(_basis) == 0: raise BasisNotFoundError('Basis not found for %s' % symb) + + # Sort basis accroding to angular momentum. This is important for method + # decontract_basis, which assumes that basis functions with the same + # angular momentum are grouped together. Related to issue #1620 #1770 + if sort_basis: + _basis = sorted([b for b in _basis if b], key=lambda b: b[0]) + fmt_basis[symb] = _basis return fmt_basis def _generate_basis_converter(): @@ -503,11 +509,6 @@ def converter(symb, raw_basis): _basis.extend(nparray_to_list(rawb)) else: _basis = nparray_to_list(raw_basis) - - # Sort basis accroding to angular momentum. This is important for method - # decontract_basis, which assumes that basis functions with the same - # angular momentum are grouped together. Related to issue #1620 #1770 - _basis = sorted([b for b in _basis if b], key=lambda b: b[0]) return _basis return converter @@ -654,6 +655,8 @@ def _to_full_contraction(mol, bas_idx): bas_idx = ib0 + numpy.where(mol._bas[ib0:ib1,ANG_OF] == l)[0] if len(bas_idx) == 0: continue + if bas_idx[0] + len(bas_idx) != bas_idx[-1] + 1: + raise NotImplementedError('Discontinuous bases of same angular momentum') mol_exps, b_coeff = _to_full_contraction(mol, bas_idx) nprim, nc = b_coeff.shape @@ -2523,11 +2526,13 @@ def build(self, dump_input=True, parse_arg=ARGPARSE, if self.verbose >= logger.WARN: self.check_sanity() - self._atom = self.format_atom(self.atom, unit=self.unit) + if self.atom: + self._atom = self.format_atom(self.atom, unit=self.unit) uniq_atoms = {a[0] for a in self._atom} - _basis = _parse_default_basis(self.basis, uniq_atoms) - self._basis = self.format_basis(_basis) + if self.basis: + _basis = _parse_default_basis(self.basis, uniq_atoms) + self._basis = self.format_basis(_basis) env = self._env[:PTR_ENV_START] self._atm, self._bas, self._env = \ self.make_env(self._atom, self._basis, env, self.nucmod, @@ -2643,30 +2648,12 @@ def _build_symmetry(self, *args, **kwargs): for ir in self.irrep_id] return self - @lib.with_doc(format_atom.__doc__) - def format_atom(self, atom, origin=0, axes=None, unit='Ang'): - return format_atom(atom, origin, axes, unit) - - @lib.with_doc(format_basis.__doc__) - def format_basis(self, basis_tab): - return format_basis(basis_tab) - - @lib.with_doc(format_pseudo.__doc__) - def format_pseudo(self, pseudo_tab): - return format_pseudo(pseudo_tab) - - @lib.with_doc(format_ecp.__doc__) - def format_ecp(self, ecp_tab): - return format_ecp(ecp_tab) - - @lib.with_doc(expand_etb.__doc__) - def expand_etb(self, l, n, alpha, beta): - return expand_etb(l, n, alpha, beta) - - @lib.with_doc(expand_etbs.__doc__) - def expand_etbs(self, etbs): - return expand_etbs(etbs) - etbs = expand_etbs + format_atom = staticmethod(format_atom) + format_basis = staticmethod(format_basis) + format_pseudo = staticmethod(format_pseudo) + format_ecp = staticmethod(format_ecp) + expand_etb = staticmethod(expand_etb) + expand_etbs = etbs = staticmethod(expand_etbs) @lib.with_doc(make_env.__doc__) def make_env(self, atoms, basis, pre_env=[], nucmod={}, nucprop=None): @@ -3713,7 +3700,7 @@ def __getattr__(self, key): '_repr_mimebundle_'): # https://github.com/mewwts/addict/issues/26 # https://github.com/jupyter/notebook/issues/2014 - raise AttributeError + raise AttributeError(f'Mole object has no attribute {key}') # Import all available modules. Some methods are registered to other # classes/modules when importing modules in __all__. @@ -3734,8 +3721,10 @@ def __getattr__(self, key): if xc in dft.XC: mf.xc = xc key = 'TDDFT' - else: + elif 'CI' in key or 'CC' in key or 'CAS' in key or 'MP' in key: mf = scf.HF(self) + else: + raise AttributeError(f'Mole object has no attribute {key}') method = getattr(mf, key) diff --git a/pyscf/lib/CMakeLists.txt b/pyscf/lib/CMakeLists.txt index d04e6444dc..98fa206365 100644 --- a/pyscf/lib/CMakeLists.txt +++ b/pyscf/lib/CMakeLists.txt @@ -222,7 +222,7 @@ if(ENABLE_XCFUN AND BUILD_XCFUN) GIT_REPOSITORY https://github.com/dftlibs/xcfun.git GIT_TAG a89b783 # This patch adds the xcfun derivative order up to 5 - PATCH_COMMAND git apply ${PROJECT_SOURCE_DIR}/libxcfun.patch + PATCH_COMMAND git apply --reject ${PROJECT_SOURCE_DIR}/libxcfun.patch PREFIX ${PROJECT_BINARY_DIR}/deps INSTALL_DIR ${PROJECT_SOURCE_DIR}/deps CMAKE_ARGS -DCMAKE_BUILD_TYPE=RELEASE -DBUILD_SHARED_LIBS=1 diff --git a/pyscf/lib/misc.py b/pyscf/lib/misc.py index 0b8e4a2dd0..c2f0c71537 100644 --- a/pyscf/lib/misc.py +++ b/pyscf/lib/misc.py @@ -1340,7 +1340,10 @@ def to_gpu(method): gpu4pyscf objects. ''' import cupy + from pyscf import gto for key, val in method.__dict__.items(): + if isinstance(val, gto.MoleBase): + continue if isinstance(val, numpy.ndarray): setattr(method, key, cupy.asarray(val)) elif hasattr(val, 'to_gpu'): diff --git a/pyscf/lo/boys.py b/pyscf/lo/boys.py index 3178458e0b..7c8cb14395 100644 --- a/pyscf/lo/boys.py +++ b/pyscf/lo/boys.py @@ -115,7 +115,10 @@ def dipole_integral(mol, mo_coeff, charge_center=None): return dip def atomic_init_guess(mol, mo_coeff): - s = mol.intor_symmetric('int1e_ovlp') + if getattr(mol, 'pbc_intor', None): + s = mol.pbc_intor('int1e_ovlp', hermi=1) + else: + s = mol.intor_symmetric('int1e_ovlp') c = orth.orth_ao(mol, s=s) mo = reduce(numpy.dot, (c.conj().T, s, mo_coeff)) # Find the AOs which have largest overlap to MOs diff --git a/pyscf/lo/pipek.py b/pyscf/lo/pipek.py index 7f24d752aa..667328a238 100644 --- a/pyscf/lo/pipek.py +++ b/pyscf/lo/pipek.py @@ -29,9 +29,10 @@ from pyscf.lib import logger from pyscf.lo import orth from pyscf.lo import boys +from pyscf.lo import pipek_jacobi from pyscf import __config__ -def atomic_pops(mol, mo_coeff, method='meta_lowdin', mf=None): +def atomic_pops(mol, mo_coeff, method='meta_lowdin', mf=None, s=None): ''' Kwargs: method : string @@ -53,10 +54,11 @@ def atomic_pops(mol, mo_coeff, method='meta_lowdin', mf=None): nmo = mo_coeff.shape[1] proj = numpy.empty((mol.natm,nmo,nmo)) - if getattr(mol, 'pbc_intor', None): # whether mol object is a cell - s = mol.pbc_intor('int1e_ovlp_sph', hermi=1) - else: - s = mol.intor_symmetric('int1e_ovlp') + if s is None: + if getattr(mol, 'pbc_intor', None): # whether mol object is a cell + s = mol.pbc_intor('int1e_ovlp', hermi=1) + else: + s = mol.intor_symmetric('int1e_ovlp') if method == 'becke': from pyscf.dft import gen_grid @@ -184,7 +186,7 @@ class PipekMezey(boys.OrbitalLocalizer): pop_method = getattr(__config__, 'lo_pipek_PM_pop_method', 'meta_lowdin') conv_tol = getattr(__config__, 'lo_pipek_PM_conv_tol', 1e-6) - exponent = getattr(__config__, 'lo_pipek_PM_exponent', 2) # should be 2 or 4 + exponent = getattr(__config__, 'lo_pipek_PM_exponent', 2) # any integer >= 2 _keys = {'pop_method', 'conv_tol', 'exponent'} @@ -199,78 +201,69 @@ def dump_flags(self, verbose=None): logger.info(self, 'pop_method = %s',self.pop_method) def gen_g_hop(self, u): + exponent = self.exponent mo_coeff = lib.dot(self.mo_coeff, u) - pop = self.atomic_pops(self.mol, mo_coeff, self.pop_method) - if self.exponent == 2: - g0 = numpy.einsum('xii,xip->pi', pop, pop) - g = -self.pack_uniq_var(g0-g0.conj().T) * 2 - elif self.exponent == 4: - pop3 = numpy.einsum('xii->xi', pop)**3 - g0 = numpy.einsum('xi,xip->pi', pop3, pop) - g = -self.pack_uniq_var(g0-g0.conj().T) * 4 - else: - raise NotImplementedError('exponent %s' % self.exponent) - - h_diag = numpy.einsum('xii,xpp->pi', pop, pop) * 2 - g_diag = g0.diagonal() - h_diag-= g_diag + g_diag.reshape(-1,1) - h_diag+= numpy.einsum('xip,xip->pi', pop, pop) * 2 - h_diag+= numpy.einsum('xip,xpi->pi', pop, pop) * 2 - h_diag = -self.pack_uniq_var(h_diag) * 2 - - g0 = g0 + g0.conj().T - if self.exponent == 2: - def h_op(x): - x = self.unpack_uniq_var(x) - norb = x.shape[0] - hx = lib.dot(x.T, g0.T).conj() - hx+= numpy.einsum('xip,xi->pi', pop, numpy.einsum('qi,xiq->xi', x, pop)) * 2 - hx-= numpy.einsum('xpp,xip->pi', pop, - lib.dot(pop.reshape(-1,norb), x).reshape(-1,norb,norb)) * 2 - hx-= numpy.einsum('xip,xp->pi', pop, numpy.einsum('qp,xpq->xp', x, pop)) * 2 - return -self.pack_uniq_var(hx-hx.conj().T) - else: - def h_op(x): - x = self.unpack_uniq_var(x) - norb = x.shape[0] - hx = lib.dot(x.T, g0.T).conj() * 2 - pop2 = numpy.einsum('xii->xi', pop)**2 - pop3 = numpy.einsum('xii->xi', pop)**3 - tmp = numpy.einsum('qi,xiq->xi', x, pop) * pop2 - hx+= numpy.einsum('xip,xi->pi', pop, tmp) * 12 - hx-= numpy.einsum('xp,xip->pi', pop3, - lib.dot(pop.reshape(-1,norb), x).reshape(-1,norb,norb)) * 4 - tmp = numpy.einsum('qp,xpq->xp', x, pop) * pop2 - hx-= numpy.einsum('xip,xp->pi', pop, tmp) * 12 - return -self.pack_uniq_var(hx-hx.conj().T) + projR = self.atomic_pops(self.mol, mo_coeff, self.pop_method).real + pop = lib.einsum('xii->xi', projR) + popexp1 = pop**(exponent-1) + popexp2 = pop**(exponent-2) + + # gradient + g = lib.einsum('xj,xij->ij', popexp1, projR) + g = -2 * exponent * self.pack_uniq_var(g - g.T) + + # hessian diagonal + g1 = lib.einsum('xi,xi->i', popexp1, pop) + g2 = lib.einsum('xi,xj->ij', popexp1, pop) + h_diag = -2 * exponent * (g1[:,None] - g2) + g1 = lib.einsum('xi,xij->ij', popexp2, projR**2.) + h_diag += 4 * exponent * (exponent-1) * g1 + h_diag = -self.pack_uniq_var(h_diag + h_diag.T) + + # hessian vector product + G = lib.einsum('xi,xij->ij', popexp1, projR) + G += G.T + + def h_op(x): + xR = self.unpack_uniq_var(x).real + + # contributions from (nabla proj) x (nabla proj) + j0 = popexp2 * lib.einsum('xik,ki->xi', projR, xR) + j1 = lib.einsum('xi,xij->ij', j0, projR) + hx = -4 * exponent * (exponent-1) * j1 + + # contributions from nabla^2 proj: symmetric terms + j1 = numpy.dot(G,xR) + hx += -exponent * j1 + + # contributions from nabla^2 proj: asymmetric terms + j1 = lib.einsum('xj,xik,kj->ij', popexp1, projR, xR) + hx += 2 * exponent * j1 + + hx -= hx.T + + return -self.pack_uniq_var(hx) return g, h_op, h_diag def get_grad(self, u=None): if u is None: u = numpy.eye(self.mo_coeff.shape[1]) + exponent = self.exponent mo_coeff = lib.dot(self.mo_coeff, u) - pop = self.atomic_pops(self.mol, mo_coeff, self.pop_method) - if self.exponent == 2: - g0 = numpy.einsum('xii,xip->pi', pop, pop) - g = -self.pack_uniq_var(g0-g0.conj().T) * 2 - else: - pop3 = numpy.einsum('xii->xi', pop)**3 - g0 = numpy.einsum('xi,xip->pi', pop3, pop) - g = -self.pack_uniq_var(g0-g0.conj().T) * 4 + projR = self.atomic_pops(self.mol, mo_coeff, self.pop_method).real + popexp1 = lib.einsum('xii->xi', projR)**(exponent-1) + g = lib.einsum('xj,xij->ij', popexp1, projR) + g = -2 * exponent * self.pack_uniq_var(g - g.T) return g def cost_function(self, u=None): if u is None: u = numpy.eye(self.mo_coeff.shape[1]) mo_coeff = lib.dot(self.mo_coeff, u) - pop = self.atomic_pops(self.mol, mo_coeff, self.pop_method) - if self.exponent == 2: - return numpy.einsum('xii,xii->', pop, pop) - else: - pop2 = numpy.einsum('xii->xi', pop)**2 - return numpy.einsum('xi,xi', pop2, pop2) + projR = self.atomic_pops(self.mol, mo_coeff, self.pop_method).real + return (lib.einsum('xii->xi', projR)**self.exponent).sum() @lib.with_doc(atomic_pops.__doc__) - def atomic_pops(self, mol, mo_coeff, method=None): + def atomic_pops(self, mol, mo_coeff, method=None, s=None): if method is None: method = self.pop_method @@ -278,7 +271,11 @@ def atomic_pops(self, mol, mo_coeff, method=None): logger.error(self, 'PM with IAO scheme should include an scf ' 'object when creating PM object.\n PM(mol, mf=scf_object)') raise ValueError('PM attribute method is not valid') - return atomic_pops(mol, mo_coeff, method, self._scf) + return atomic_pops(mol, mo_coeff, method, self._scf, s=s) + + def stability_jacobi(self): + return pipek_jacobi.PipekMezey_stability_jacobi(self) + PM = Pipek = PipekMezey @@ -287,12 +284,21 @@ def atomic_pops(self, mol, mo_coeff, method=None): mol = gto.Mole() mol.atom = ''' - O 0. 0. 0.2 - H 0. -0.5 -0.4 - H 0. 0.5 -0.4 - ''' + O 0.00000 0.00000 0.11779 + H 0.00000 0.75545 -0.47116 + H 0.00000 -0.75545 -0.47116 + ''' mol.basis = 'ccpvdz' mol.build() mf = scf.RHF(mol).run() - mo = PM(mol).kernel(mf.mo_coeff[:,5:9], verbose=4) + mlo = PM(mol) + mlo.verbose = 4 + mlo.exponent = 2 # integer >= 2 + mo0 = mf.mo_coeff[:,mf.mo_occ>1e-6] + mo = mlo.kernel(mo0) + isstable, mo1 = mlo.stability_jacobi() + if not isstable: + mo = mlo.kernel(mo1) + isstable, mo1 = mlo.stability_jacobi() + assert( isstable ) diff --git a/pyscf/lo/pipek_jacobi.py b/pyscf/lo/pipek_jacobi.py new file mode 100644 index 0000000000..dce9d79f38 --- /dev/null +++ b/pyscf/lo/pipek_jacobi.py @@ -0,0 +1,161 @@ +#!/usr/bin/env python +# Copyright 2020-2023 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. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +# Author: Hong-Zhou Ye +# + + +import numpy as np + +from pyscf import lib + + +def PipekMezey_stability_jacobi(mlo, mo_coeff=None, conv_tol=None): + ''' Perform a stability analysis using Jacobi sweep. + + The function first identifies all pairs whose mix may lead to a loss increase greater + than `conv_tol`. If no pairs are found, the function returns the input mo_coeff. + Otherwise, it performs one Jacobi sweep for all pairs identified above and return the + new mo_coeff. + + Note: + All analysis inside this function is performed with exponent = 2, which is generally + good enough for hopping out of local minimum. + + Args: + mlo : PipekMezey object + mo_coeff : np.ndarray + conv_tol : float + If mixing a pair of MOs leads to a loss increase greater than this parameter, + the pair of MOs will be sent to a subsequent Jacob sweep. + + Returns: + A tuple with two entries: + + isstable : bool + If the input MOs are stable under Jacob rotation + mo_coeff_new : np.ndarray + The MOs after Jacobi rotation. If isstable is False, mo_coeff_new is the same + as input mo_coeff. + ''' + log = lib.logger.new_logger(mlo) + + exponent = mlo.exponent + if exponent != 2: + log.warn('Jacobi sweep for PipekMezey with exponent = %d is not implemented. ' + 'Exponent = 2 will be used instead.', exponent) + + if mo_coeff is None: mo_coeff = mlo.mo_coeff + if conv_tol is None: conv_tol = mlo.conv_tol + + mol = mlo.mol + nmo = mo_coeff.shape[1] + + if getattr(mol, 'pbc_intor', None): + s = mol.pbc_intor('int1e_ovlp', hermi=1) + else: + s = mol.intor_symmetric('int1e_ovlp') + + proj = mlo.atomic_pops(mol, mo_coeff, s=s) + pop = lib.einsum('xii->xi', proj) + + idx_tril = np.tril_indices(nmo,-1) + def pack_tril(A): + if A.ndim == 2: + return A[idx_tril] + elif A.ndim == 3: + return np.asarray([A[i][idx_tril] for i in range(A.shape[0])]) + else: + raise RuntimeError + s_proj = pack_tril(proj) + d_proj = pack_tril(lib.direct_sum('ai-aj->aij', pop, pop))*0.5 + m_proj = pack_tril(lib.direct_sum('ai+aj->aij', pop, pop))*0.5 + loss0_allpair = pack_tril(lib.direct_sum('ai+aj->ij', pop**2, pop**2)) + A_allpair = (s_proj**2 - d_proj**2.).sum(axis=0) + B_allpair = (s_proj*d_proj).sum(axis=0) + C_allpair = (2*m_proj**2+s_proj**2+d_proj**2).sum(axis=0) + t_allpair = np.arctan(2*B_allpair/A_allpair) + loss1_allpair = C_allpair - A_allpair*np.cos(4*t_allpair) - 2*B_allpair*np.sin(4*t_allpair) + T_allpair = t_allpair + np.pi*0.25 + loss2_allpair = C_allpair - A_allpair*np.cos(4*T_allpair) - 2*B_allpair*np.sin(4*T_allpair) + loss_allpair = np.maximum(loss1_allpair, loss2_allpair) + dloss_allpair = loss_allpair - loss0_allpair + mo_pair_indices = np.where(dloss_allpair > conv_tol)[0] + dloss_pair = dloss_allpair[mo_pair_indices] + mo_pair_indices = mo_pair_indices[np.argsort(dloss_pair)[::-1]] + + rows = np.floor(((1+8*mo_pair_indices)**0.5-1)*0.5).astype(int) + cols = mo_pair_indices - rows*(rows+1) // 2 + rows += 1 + pairs = [(row,col) for row,col in zip(rows,cols)] + + if len(pairs) > 0: + log.info('Jacobi sweep for %d mo pairs: %s', len(pairs), pairs) + mo_coeff = _jacobi_sweep(mlo, mo_coeff, pairs, conv_tol, s=s) + return False, mo_coeff + else: + return True, mo_coeff + +def _jacobi_sweep(mlo, mo_coeff, mo_pairs, conv_tol, s=None): + mol = mlo.mol + mo_coeff = mo_coeff.copy() + + def get_loss_pair(m12s, d12s, s12s, theta): + c2t = np.cos(2*theta) + s2t = np.sin(2*theta) + loss_pair = np.asarray([((m12s+d12s*c2t-s12s*s2t)**2).sum(), + ((m12s-d12s*c2t+s12s*s2t)**2).sum()]) + return loss_pair.sum(), loss_pair + def rotate_pair(mo_coeff, i, j, theta): + ct = np.cos(theta) + st = np.sin(theta) + mo1 = mo_coeff[:,i]*ct - mo_coeff[:,j]*st + mo2 = mo_coeff[:,i]*st + mo_coeff[:,j]*ct + mo_coeff[:,i] = mo1 + mo_coeff[:,j] = mo2 + return mo_coeff + + if s is None: + if getattr(mol, 'pbc_intor', None): + s = mol.pbc_intor('int1e_ovlp', hermi=1) + else: + s = mol.intor_symmetric('int1e_ovlp') + + for (idx,jdx) in mo_pairs: + mo_pair = mo_coeff[:,[idx,jdx]] + + pop = mlo.atomic_pops(mol, mo_pair, mlo.pop_method, s=s) + n1s = pop[:,0,0] + n2s = pop[:,1,1] + m12s = 0.5*(n1s + n2s) + d12s = 0.5*(n1s - n2s) + s12s = pop[:,0,1].real + l0 = (n1s**2 + n2s**2).sum() + + A = (s12s**2 - d12s**2).sum() + B = (s12s*d12s).sum() + + t1 = 0.25 * np.arctan(2*B / A) + t2 = t1 + np.pi*0.25 + l1, l1s = get_loss_pair(m12s, d12s, s12s, t1) + l2, l2s = get_loss_pair(m12s, d12s, s12s, t2) + + dl1 = l1 - l0 + dl2 = l2 - l0 + if dl1 > conv_tol or dl2 > conv_tol: + (topt, ls) = (t1, l1s) if l1 > l2 else (t2, l2s) + mo_coeff = rotate_pair(mo_coeff, idx, jdx, topt) + + return mo_coeff diff --git a/pyscf/lo/test/test_pipek.py b/pyscf/lo/test/test_pipek.py new file mode 100644 index 0000000000..c288564310 --- /dev/null +++ b/pyscf/lo/test/test_pipek.py @@ -0,0 +1,146 @@ +#!/usr/bin/env python +# Copyright 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. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# + +import unittest +import numpy as np +from pyscf import gto, lo + + +class Water(unittest.TestCase): + @classmethod + def setUpClass(cls): + mol = gto.Mole() + mol.verbose = 4 + mol.output = '/dev/null' + mol.atom = ''' + O 0.000000 0.000000 0.000000 + O 0.000000 0.000000 1.480000 + H 0.895669 0.000000 -0.316667 + H -0.895669 0.000000 1.796667 + ''' + mol.basis = 'sto-3g' + mol.precision = 1e-10 + mol.build() + + cls.mol = mol + + @classmethod + def tearDownClass(cls): + cls.mol.stdout.close() + del cls.mol + + def test_cost(self): + ''' Test for `cost_function` + ''' + mol = self.mol + s = mol.intor_symmetric('int1e_ovlp') + mo_coeff = lo.orth.schmidt(s) + norb = mo_coeff.shape[1] + u = np.eye(norb) + + loss_ref = { + 2: 11.0347269392, + 3: 10.5595998735, + 4: 10.1484720537 + } + + mlo = lo.pipek.PipekMezey(mol, mo_coeff) + for exponent in [2,3,4]: + mlo.set(exponent=exponent) + loss = mlo.cost_function(u) + self.assertAlmostEqual(loss, loss_ref[exponent], 6) + + def test_grad_hess(self): + ''' Test for `get_grad` and `gen_g_hop` + ''' + mol = self.mol + s = mol.intor_symmetric('int1e_ovlp') + mo_coeff = lo.orth.schmidt(s) + + mo_idx = [1,6,9] + mo_coeff = mo_coeff[:,mo_idx] + norb = mo_coeff.shape[1] + u0 = np.eye(norb) + + mlo = lo.pipek.PipekMezey(mol, mo_coeff) + x0 = mlo.pack_uniq_var(np.zeros((norb,norb))) + + step_length = 1e-3 + precision = 4 + + for exponent in [2,3,4]: + mlo.set(exponent=exponent) + g = mlo.get_grad(u0) + g1, h_op, h_diag = mlo.gen_g_hop(u0) + + self.assertAlmostEqual(abs(g-g1).max(), 0, 6) + + H = np.zeros((norb,norb)) + for i in range(norb): + xi = np.zeros(norb) + xi[i] = 1 + H[:,i] = h_op(xi) + + self.assertAlmostEqual(abs(np.diagonal(H)-h_diag).max(), 0, 6) + + def func(x): + u = mlo.extract_rotation(x) + return -mlo.cost_function(u) + + num_g = _num_grad(func, x0, step_length) + num_H = _num_hess(func, x0, step_length) + + self.assertAlmostEqual(abs(g-num_g).max(), 0, precision) + self.assertAlmostEqual(abs(H-num_H).max(), 0, precision) + + +def _num_grad(func, x0, step_length): + x0 = np.asarray(x0) + n = x0.size + g = np.zeros_like(x0) + for i in range(n): + dx = np.zeros(n) + dx[i] = step_length + yf = func(x0+dx) + yb = func(x0-dx) + g[i] = (yf-yb) / (2.*step_length) + return g +def _num_hess(func, x0, step_length): + x0 = np.asarray(x0) + n = x0.size + H = np.zeros((n,n)) + y0 = func(x0) + for i in range(n): + dxi = np.zeros(n) + dxi[i] = step_length + for j in range(i+1,n): + dxj = np.zeros(n) + dxj[j] = step_length + yff = func(x0+dxi+dxj) + yfb = func(x0+dxi-dxj) + ybf = func(x0-dxi+dxj) + ybb = func(x0-dxi-dxj) + H[i,j] = H[j,i] = (yff+ybb-yfb-ybf) / (4.*step_length**2.) + yf = func(x0+dxi) + yb = func(x0-dxi) + H[i,i] = (yf+yb-2*y0) / step_length**2. + + return H + + +if __name__ == "__main__": + print("Full Tests for PipekMezey") + unittest.main() diff --git a/pyscf/pbc/df/df_jk.py b/pyscf/pbc/df/df_jk.py index 535db5ea18..88c782b311 100644 --- a/pyscf/pbc/df/df_jk.py +++ b/pyscf/pbc/df/df_jk.py @@ -296,13 +296,15 @@ def get_k_kpts(mydf, dm_kpts, hermi=1, kpts=numpy.zeros((1,3)), kpts_band=None, skmoR = skmo2R = None if not mydf.force_dm_kbuild: if mo_coeff is not None: - if isinstance(mo_coeff[0], (list, tuple)): + if isinstance(mo_coeff[0], (list, tuple)) or (isinstance(mo_coeff[0], numpy.ndarray) + and mo_coeff[0].ndim == 3): mo_coeff = [mo for mo1 in mo_coeff for mo in mo1] if len(mo_coeff) != nset*nkpts: # wrong shape log.warn('mo_coeff from dm tag has wrong shape. ' 'Calculating mo from dm instead.') mo_coeff = None - elif isinstance(mo_occ[0], (list, tuple)): + elif isinstance(mo_occ[0], (list, tuple)) or (isinstance(mo_occ[0], numpy.ndarray) + and mo_occ[0].ndim == 2): mo_occ = [mo for mo1 in mo_occ for mo in mo1] if mo_coeff is not None: skmoR, skmoI = _format_mo(mo_coeff, mo_occ, shape=(nset,nkpts), order='F', @@ -704,13 +706,15 @@ def get_k_kpts_kshift(mydf, dm_kpts, kshift, hermi=0, kpts=numpy.zeros((1,3)), k skmoR = skmo2R = None if not mydf.force_dm_kbuild: if mo_coeff is not None: - if isinstance(mo_coeff[0], (list, tuple)): + if isinstance(mo_coeff[0], (list, tuple)) or (isinstance(mo_coeff[0], numpy.ndarray) + and mo_coeff[0].ndim == 3): mo_coeff = [mo for mo1 in mo_coeff for mo in mo1] if len(mo_coeff) != nset*nkpts: # wrong shape log.warn('mo_coeff from dm tag has wrong shape. ' 'Calculating mo from dm instead.') mo_coeff = None - elif isinstance(mo_occ[0], (list, tuple)): + elif isinstance(mo_occ[0], (list, tuple)) or (isinstance(mo_occ[0], numpy.ndarray) + and mo_occ[0].ndim == 2): mo_occ = [mo for mo1 in mo_occ for mo in mo1] if mo_coeff is not None: skmoR, skmoI = _format_mo(mo_coeff, mo_occ, shape=(nset,nkpts), order='F', diff --git a/pyscf/pbc/gto/cell.py b/pyscf/pbc/gto/cell.py index 33891fd6c7..813cabe258 100644 --- a/pyscf/pbc/gto/cell.py +++ b/pyscf/pbc/gto/cell.py @@ -971,7 +971,7 @@ def __getattr__(self, key): '_repr_mimebundle_'): # https://github.com/mewwts/addict/issues/26 # https://github.com/jupyter/notebook/issues/2014 - raise AttributeError + raise AttributeError(f'Cell object has no attribute {key}') # Import all available modules. Some methods are registered to other # classes/modules when importing modules in __all__. @@ -994,8 +994,10 @@ def __getattr__(self, key): if xc in XC: mf.xc = xc key = 'KTDDFT' - else: + elif 'CI' in key or 'CC' in key or 'MP' in key: mf = scf.KHF(self) + else: + raise AttributeError(f'Cell object has no attribute {key}') # Remove prefix 'K' because methods are registered without the leading 'K' key = key[1:] else: @@ -1008,8 +1010,10 @@ def __getattr__(self, key): if xc in XC: mf.xc = xc key = 'TDDFT' - else: + elif 'CI' in key or 'CC' in key or 'MP' in key: mf = scf.HF(self) + else: + raise AttributeError(f'Cell object has no attribute {key}') method = getattr(mf, key) diff --git a/pyscf/solvent/_attach_solvent.py b/pyscf/solvent/_attach_solvent.py index aeda9da86c..ab5d58dc7c 100644 --- a/pyscf/solvent/_attach_solvent.py +++ b/pyscf/solvent/_attach_solvent.py @@ -383,7 +383,7 @@ def casci_iter_(ci0, log): de = e_tot - e_last if isinstance(e_cas, (float, numpy.number)): - log.info('Sovlent cycle %d E(CASCI+solvent) = %.15g ' + log.info('Solvent cycle %d E(CASCI+solvent) = %.15g ' 'dE = %g', cycle, e_tot, de) else: for i, e in enumerate(e_tot): @@ -523,7 +523,7 @@ def kernel(self, *args, **kwargs): with_solvent.e, with_solvent.v = with_solvent.kernel(dm) de = e_tot - e_last - log.info('Sovlent cycle %d E_tot = %.15g dE = %g', + log.info('Solvent cycle %d E_tot = %.15g dE = %g', cycle, e_tot, de) if abs(e_tot-e_last).max() < with_solvent.conv_tol: diff --git a/pyscf/tools/molden.py b/pyscf/tools/molden.py index 8d0ac4e48f..07862101be 100644 --- a/pyscf/tools/molden.py +++ b/pyscf/tools/molden.py @@ -239,8 +239,16 @@ def read_one_bas(lsym, nb, fac=1): elif dat[0].upper() in 'SPDFGHIJ': basis[symb].append(read_one_bas(*dat)) - mol.basis = envs['basis'] = basis mol.atom = [atoms[i] for i in atom_seq] + uniq_atoms = {a[0] for a in mol.atom} + + # To avoid the mol.build() sort the basis, disable mol.basis and set the + # internal data _basis directly. It is a workaround to solve issue #1961. + # Mole.decontract_basis function should be rewritten to support + # discontinuous bases that have the same angular momentum. + mol.basis = {} + _basis = gto.mole._parse_default_basis(basis, uniq_atoms) + mol._basis = envs['basis'] = gto.format_basis(_basis, sort_basis=False) return mol def _parse_mo(lines, envs): @@ -353,6 +361,8 @@ def load(moldenfile, verbose=0): else: sys.stderr.write('Unknown section %s\n' % sec_title) + mo_energy, mo_coeff, mo_occ, irrep_labels, spins = None, None, None, None, None + for sec_kind in _SEC_ORDER: if sec_kind == 'MO' and 'MO' in sec_kinds: if len(sec_kinds['MO']) == 1: @@ -392,6 +402,8 @@ def load(moldenfile, verbose=0): if isinstance(mo_occ, tuple): mol.spin = int(mo_occ[0].sum() - mo_occ[1].sum()) + if not mol._built: + mol.build(0, 0) return mol, mo_energy, mo_coeff, mo_occ, irrep_labels, spins parse = read = load @@ -503,8 +515,9 @@ def remove_high_l(mol, mo_coeff=None): ''' pmol = mol.copy() pmol.basis = {} + pmol._basis = {} for symb, bas in mol._basis.items(): - pmol.basis[symb] = [b for b in bas if b[0] <= 4] + pmol._basis[symb] = [b for b in bas if b[0] <= 4] pmol.build(0, 0) if mo_coeff is None: return pmol, None diff --git a/pyscf/tools/test/test_molden.py b/pyscf/tools/test/test_molden.py index 9d0cca7f66..944c7c9178 100644 --- a/pyscf/tools/test/test_molden.py +++ b/pyscf/tools/test/test_molden.py @@ -82,6 +82,53 @@ def test_dump_cartesian_gto_symm_orbital(self): mo_coeff = res[2] self.assertAlmostEqual(abs(mf.mo_coeff-mo_coeff).max(), 0, 12) + def test_basis_not_sorted(self): + with tempfile.NamedTemporaryFile('w') as ftmp: + ftmp.write('''\ +[Molden Format] +made by pyscf v[2.4.0] +[Atoms] (AU) +Ne 1 10 0.00000000000000 0.00000000000000 0.00000000000000 +[GTO] +1 0 + s 8 1.00 + 17880 0.00073769817327139 + 2683 0.0056746782244738 + 611.5 0.028871187450674 + 173.5 0.10849560938601 + 56.64 0.29078802505672 + 20.42 0.44814064476114 + 7.81 0.25792047270532 + 1.653 0.015056839544698 + s 8 1.00 + 17880 -0.00033214909943481 + 2683 -0.0026205019065875 + 611.5 -0.013009816761002 + 173.5 -0.053420003125961 + 56.64 -0.14716522424261 + 20.42 -0.33838075724805 + 7.81 -0.20670101921688 + 1.653 1.0950299234565 + s 1 1.00 + 0.4869 1 + p 3 1.00 + 28.39 0.066171986640049 + 6.27 0.34485329752845 + 1.695 0.73045763818875 + p 1 1.00 + 0.4317 1 + d 1 1.00 + 2.202 1 + s 1 1.00 + 1 1 + +[5d] +[7f] +[9g] +''') + ftmp.flush() + mol = molden.load(ftmp.name)[0] + self.assertEqual(mol._bas[:,1].tolist(), [0, 0, 0, 1, 1, 2, 0]) if __name__ == "__main__": print("Full Tests for molden")