diff --git a/pyscf/pbc/scf/ghf.py b/pyscf/pbc/scf/ghf.py index fc5830006e..e88916e99e 100644 --- a/pyscf/pbc/scf/ghf.py +++ b/pyscf/pbc/scf/ghf.py @@ -117,7 +117,7 @@ def __init__(self, cell, kpt=np.zeros(3), _finalize = mol_ghf.GHF._finalize analyze = lib.invalid_method('analyze') mulliken_pop = lib.invalid_method('mulliken_pop') - mulliken_meta = lib.invalid_method('mulliken_meta') + mulliken_meta = mol_ghf.GHF.mulliken_meta spin_square = mol_ghf.GHF.spin_square stability = mol_ghf.GHF.stability diff --git a/pyscf/pbc/scf/kghf.py b/pyscf/pbc/scf/kghf.py index aa0f88c3ca..461a30dabb 100644 --- a/pyscf/pbc/scf/kghf.py +++ b/pyscf/pbc/scf/kghf.py @@ -136,38 +136,41 @@ def get_occ(mf, mo_energy_kpts=None, mo_coeff_kpts=None): return mo_occ_kpts -def mulliken_meta(cell, dm_ao_kpts, verbose=logger.DEBUG, +def _make_rdm1_meta(cell, dm_ao_kpts, kpts, pre_orth_method, s): + from pyscf.lo import orth + from pyscf.pbc.tools import k2gamma + + kmesh = k2gamma.kpts_to_kmesh(cell, kpts-kpts[0]) + nkpts, nso = dm_ao_kpts.shape[:2] + nao = nso // 2 + scell, phase = k2gamma.get_phase(cell, kpts, kmesh) + s_sc = k2gamma.to_supercell_ao_integrals(cell, kpts, s, kmesh=kmesh, force_real=False) + orth_coeff = orth.orth_ao(scell, 'meta_lowdin', pre_orth_method, s=s_sc)[:,:nao] # cell 0 only + c_inv = np.dot(orth_coeff.T.conj(), s_sc) + c_inv = lib.einsum('aRp,Rk->kap', c_inv.reshape(nao,nkpts,nao), phase) + dm_aa = lib.einsum('kap,kpq,kbq->ab', c_inv, dm_ao_kpts[:,:nao,:nao], c_inv.conj()) + dm_bb = lib.einsum('kap,kpq,kbq->ab', c_inv, dm_ao_kpts[:,nao:,nao:], c_inv.conj()) + + return (dm_aa, dm_bb) + +def mulliken_meta(cell, dm_ao_kpts, kpts, verbose=logger.DEBUG, pre_orth_method=PRE_ORTH_METHOD, s=None): '''A modified Mulliken population analysis, based on meta-Lowdin AOs. - - Note this function only computes the Mulliken population for the gamma - point density matrix. + The results are equivalent to the corresponding supercell calculation. ''' - from pyscf.lo import orth - if s is None: - s = khf.get_ovlp(cell) log = logger.new_logger(cell, verbose) - log.note('Analyze output for *gamma point*.') - log.info(' To include the contributions from k-points, transform to a ' - 'supercell then run the population analysis on the supercell\n' - ' from pyscf.pbc.tools import k2gamma\n' - ' k2gamma.k2gamma(mf).mulliken_meta()') - log.note("KGHF mulliken_meta") - dm_ao_gamma = dm_ao_kpts[0,:,:] - nso = dm_ao_gamma.shape[-1] - nao = nso // 2 - dm_ao_gamma_aa = dm_ao_gamma[:nao,:nao] - dm_ao_gamma_bb = dm_ao_gamma[nao:,nao:] - s_gamma = s[0,:,:].real - s_gamma_aa = s_gamma[:nao,:nao] - orth_coeff = orth.orth_ao(cell, 'meta_lowdin', pre_orth_method, s=s_gamma_aa) - c_inv = np.dot(orth_coeff.T, s_gamma_aa) - dm_aa = reduce(np.dot, (c_inv, dm_ao_gamma_aa, c_inv.T.conj())) - dm_bb = reduce(np.dot, (c_inv, dm_ao_gamma_bb, c_inv.T.conj())) + if s is None: + s = khf.get_ovlp(None, cell=cell, kpts=kpts) + if s is not None: + if s[0].shape == dm_ao_kpts[0].shape: # s in SO + nao = dm_ao_kpts[0].shape[0]//2 + s = lib.asarray(s[:,:nao,:nao], order='C') # keep only one spin sector + + dm_aa, dm_bb = _make_rdm1_meta(cell, dm_ao_kpts, kpts, pre_orth_method, s) log.note(' ** Mulliken pop alpha/beta on meta-lowdin orthogonal AOs **') - return mol_uhf.mulliken_pop(cell, (dm_aa,dm_bb), np.eye(orth_coeff.shape[0]), log) + return mol_uhf.mulliken_pop(cell, (dm_aa,dm_bb), np.eye(dm_aa.shape[0]), log) def _cast_mol_init_guess(fn): def fn_init_guess(mf, cell=None, kpts=None): @@ -259,12 +262,13 @@ def get_bands(self, kpts_band, cell=None, dm_kpts=None, kpts=None): raise NotImplementedError @lib.with_doc(mulliken_meta.__doc__) - def mulliken_meta(self, cell=None, dm=None, verbose=logger.DEBUG, + def mulliken_meta(self, cell=None, dm=None, kpts=None, verbose=logger.DEBUG, pre_orth_method=PRE_ORTH_METHOD, s=None): if cell is None: cell = self.cell if dm is None: dm = self.make_rdm1() - if s is None: s = self.get_ovlp(cell) - return mulliken_meta(cell, dm, s=s, verbose=verbose, + if kpts is None: kpts = self.kpts + if s is None: s = khf.get_ovlp(self, cell, kpts) + return mulliken_meta(cell, dm, kpts, s=s, verbose=verbose, pre_orth_method=pre_orth_method) def mulliken_pop(self): diff --git a/pyscf/pbc/scf/khf.py b/pyscf/pbc/scf/khf.py index e0846c73c9..8024a45953 100644 --- a/pyscf/pbc/scf/khf.py +++ b/pyscf/pbc/scf/khf.py @@ -278,31 +278,36 @@ def analyze(mf, verbose=None, with_meta_lowdin=WITH_META_LOWDIN, #return mf.mulliken_pop(mf.cell, dm, s=ovlp_ao, verbose=verbose) -def mulliken_meta(cell, dm_ao_kpts, verbose=logger.DEBUG, +def _make_rdm1_meta(cell, dm_ao_kpts, kpts, pre_orth_method, s): + from pyscf.lo import orth + from pyscf.pbc.tools import k2gamma + + kmesh = k2gamma.kpts_to_kmesh(cell, kpts-kpts[0]) + nkpts, nao = dm_ao_kpts.shape[:2] + scell, phase = k2gamma.get_phase(cell, kpts, kmesh) + s_sc = k2gamma.to_supercell_ao_integrals(cell, kpts, s, kmesh=kmesh, force_real=False) + orth_coeff = orth.orth_ao(scell, 'meta_lowdin', pre_orth_method, s=s_sc)[:,:nao] # cell 0 only + c_inv = np.dot(orth_coeff.T.conj(), s_sc) + c_inv = lib.einsum('aRp,Rk->kap', c_inv.reshape(nao,nkpts,nao), phase) + dm = lib.einsum('kap,kpq,kbq->ab', c_inv, dm_ao_kpts, c_inv.conj()) + + return dm + + +def mulliken_meta(cell, dm_ao_kpts, kpts, verbose=logger.DEBUG, pre_orth_method=PRE_ORTH_METHOD, s=None): '''A modified Mulliken population analysis, based on meta-Lowdin AOs. - - Note this function only computes the Mulliken population for the gamma - point density matrix. + The results are equivalent to the corresponding supercell calculation. ''' - from pyscf.lo import orth - if s is None: - s = get_ovlp(cell) log = logger.new_logger(cell, verbose) - log.note('Analyze output for *gamma point*') - log.info(' To include the contributions from k-points, transform to a ' - 'supercell then run the population analysis on the supercell\n' - ' from pyscf.pbc.tools import k2gamma\n' - ' k2gamma.k2gamma(mf).mulliken_meta()') - log.note("KRHF mulliken_meta") - dm_ao_gamma = dm_ao_kpts[0,:,:].real - s_gamma = s[0,:,:].real - orth_coeff = orth.orth_ao(cell, 'meta_lowdin', pre_orth_method, s=s_gamma) - c_inv = np.dot(orth_coeff.T, s_gamma) - dm = reduce(np.dot, (c_inv, dm_ao_gamma, c_inv.T.conj())) + + if s is None: + s = get_ovlp(None, cell=cell, kpts=kpts) + + dm = _make_rdm1_meta(cell, dm_ao_kpts, kpts, pre_orth_method, s) log.note(' ** Mulliken pop on meta-lowdin orthogonal AOs **') - return mol_hf.mulliken_pop(cell, dm, np.eye(orth_coeff.shape[0]), log) + return mol_hf.mulliken_pop(cell, dm, np.eye(dm.shape[0]), log) def canonicalize(mf, mo_coeff_kpts, mo_occ_kpts, fock=None): @@ -633,7 +638,7 @@ def dump_chk(self, envs_or_file): def analyze(mf, verbose=None, with_meta_lowdin=WITH_META_LOWDIN, **kwargs): raise NotImplementedError - def mulliken_meta(self, cell=None, dm=None, verbose=logger.DEBUG, + def mulliken_meta(self, cell=None, dm=None, kpts=None, verbose=logger.DEBUG, pre_orth_method=PRE_ORTH_METHOD, s=None): raise NotImplementedError @@ -738,12 +743,13 @@ def get_init_guess(self, cell=None, key='minao', s1e=None): return dm_kpts @lib.with_doc(mulliken_meta.__doc__) - def mulliken_meta(self, cell=None, dm=None, verbose=logger.DEBUG, + def mulliken_meta(self, cell=None, dm=None, kpts=None, verbose=logger.DEBUG, pre_orth_method=PRE_ORTH_METHOD, s=None): if cell is None: cell = self.cell if dm is None: dm = self.make_rdm1() + if kpts is None: kpts = self.kpts if s is None: s = self.get_ovlp(cell) - return mulliken_meta(cell, dm, s=s, verbose=verbose, + return mulliken_meta(cell, dm, kpts, s=s, verbose=verbose, pre_orth_method=pre_orth_method) def nuc_grad_method(self): diff --git a/pyscf/pbc/scf/krohf.py b/pyscf/pbc/scf/krohf.py index b6b6da83a3..bb2cefbc93 100644 --- a/pyscf/pbc/scf/krohf.py +++ b/pyscf/pbc/scf/krohf.py @@ -211,15 +211,12 @@ def get_occ(mf, mo_energy_kpts=None, mo_coeff_kpts=None): @lib.with_doc(khf.mulliken_meta.__doc__) -def mulliken_meta(cell, dm_ao_kpts, verbose=logger.DEBUG, +def mulliken_meta(cell, dm_ao_kpts, kpts, verbose=logger.DEBUG, pre_orth_method=PRE_ORTH_METHOD, s=None): '''Mulliken population analysis, based on meta-Lowdin AOs. - - Note this function only computes the Mulliken population for the gamma - point density matrix. ''' dm = dm_ao_kpts[0] + dm_ao_kpts[1] - return khf.mulliken_meta(cell, dm, verbose, pre_orth_method, s) + return khf.mulliken_meta(cell, dm, kpts, verbose, pre_orth_method, s) def canonicalize(mf, mo_coeff_kpts, mo_occ_kpts, fock=None): @@ -362,12 +359,13 @@ def init_guess_by_chkfile(self, chk=None, project=True, kpts=None): if kpts is None: kpts = self.kpts return init_guess_by_chkfile(self.cell, chk, project, kpts) - def mulliken_meta(self, cell=None, dm=None, verbose=logger.DEBUG, + def mulliken_meta(self, cell=None, dm=None, kpts=None, verbose=logger.DEBUG, pre_orth_method=PRE_ORTH_METHOD, s=None): if cell is None: cell = self.cell if dm is None: dm = self.make_rdm1() + if kpts is None: kpts = self.kpts if s is None: s = self.get_ovlp(cell) - return mulliken_meta(cell, dm, s=s, verbose=verbose, + return mulliken_meta(cell, dm, kpts, s=s, verbose=verbose, pre_orth_method=pre_orth_method) def stability(self, diff --git a/pyscf/pbc/scf/kuhf.py b/pyscf/pbc/scf/kuhf.py index 8827cbad6e..259d848d83 100644 --- a/pyscf/pbc/scf/kuhf.py +++ b/pyscf/pbc/scf/kuhf.py @@ -216,32 +216,52 @@ def energy_elec(mf, dm_kpts=None, h1e_kpts=None, vhf_kpts=None): return (e1+e_coul).real, e_coul.real -def mulliken_meta(cell, dm_ao_kpts, verbose=logger.DEBUG, +def _make_rdm1_meta(cell, dm_ao_kpts, kpts, pre_orth_method, s): + from pyscf.lo import orth + from pyscf.pbc.tools import k2gamma + + kmesh = k2gamma.kpts_to_kmesh(cell, kpts-kpts[0]) + nkpts, nao = dm_ao_kpts[0].shape[:2] + scell, phase = k2gamma.get_phase(cell, kpts, kmesh) + s_sc = k2gamma.to_supercell_ao_integrals(cell, kpts, s, kmesh=kmesh, force_real=False) + orth_coeff = orth.orth_ao(scell, 'meta_lowdin', pre_orth_method, s=s_sc)[:,:nao] # cell 0 only + c_inv = np.dot(orth_coeff.T.conj(), s_sc) + c_inv = lib.einsum('aRp,Rk->kap', c_inv.reshape(nao,nkpts,nao), phase) + dm_a = lib.einsum('kap,kpq,kbq->ab', c_inv, dm_ao_kpts[0], c_inv.conj()) + dm_b = lib.einsum('kap,kpq,kbq->ab', c_inv, dm_ao_kpts[1], c_inv.conj()) + + return (dm_a, dm_b) + + +def mulliken_meta(cell, dm_ao_kpts, kpts, verbose=logger.DEBUG, pre_orth_method=PRE_ORTH_METHOD, s=None): '''A modified Mulliken population analysis, based on meta-Lowdin AOs. - - Note this function only computes the Mulliken population for the gamma - point density matrix. + The results are equivalent to the corresponding supercell calculation. ''' - from pyscf.lo import orth + log = logger.new_logger(cell, verbose) + if s is None: - s = khf.get_ovlp(cell) + s = khf.get_ovlp(None, cell=cell, kpts=kpts) + + dm_a, dm_b = _make_rdm1_meta(cell, dm_ao_kpts, kpts, pre_orth_method, s) + + log.note(' ** Mulliken pop alpha/beta on meta-lowdin orthogonal AOs **') + return mol_uhf.mulliken_pop(cell, (dm_a,dm_b), np.eye(dm_a.shape[0]), log) + + +def mulliken_meta_spin(cell, dm_ao_kpts, kpts, verbose=logger.DEBUG, + pre_orth_method=PRE_ORTH_METHOD, s=None): + '''A modified Mulliken population analysis, based on meta-Lowdin AOs. + ''' log = logger.new_logger(cell, verbose) - log.note('Analyze output for *gamma point*.') - log.info(' To include the contributions from k-points, transform to a ' - 'supercell then run the population analysis on the supercell\n' - ' from pyscf.pbc.tools import k2gamma\n' - ' k2gamma.k2gamma(mf).mulliken_meta()') - log.note("KUHF mulliken_meta") - dm_ao_gamma = dm_ao_kpts[:,0,:,:].real - s_gamma = s[0,:,:].real - orth_coeff = orth.orth_ao(cell, 'meta_lowdin', pre_orth_method, s=s_gamma) - c_inv = np.dot(orth_coeff.T, s_gamma) - dm_a = reduce(np.dot, (c_inv, dm_ao_gamma[0], c_inv.T.conj())) - dm_b = reduce(np.dot, (c_inv, dm_ao_gamma[1], c_inv.T.conj())) + + if s is None: + s = khf.get_ovlp(None, cell=cell, kpts=kpts) + + dm_a, dm_b = _make_rdm1_meta(cell, dm_ao_kpts, kpts, pre_orth_method, s) log.note(' ** Mulliken pop alpha/beta on meta-lowdin orthogonal AOs **') - return mol_uhf.mulliken_pop(cell, (dm_a,dm_b), np.eye(orth_coeff.shape[0]), log) + return mol_uhf.mulliken_spin_pop(cell, (dm_a,dm_b), np.eye(dm_a.shape[0]), log) def canonicalize(mf, mo_coeff_kpts, mo_occ_kpts, fock=None): @@ -509,14 +529,25 @@ def init_guess_by_chkfile(self, chk=None, project=True, kpts=None): return init_guess_by_chkfile(self.cell, chk, project, kpts) @lib.with_doc(mulliken_meta.__doc__) - def mulliken_meta(self, cell=None, dm=None, verbose=logger.DEBUG, + def mulliken_meta(self, cell=None, dm=None, kpts=None, verbose=logger.DEBUG, pre_orth_method=PRE_ORTH_METHOD, s=None): if cell is None: cell = self.cell if dm is None: dm = self.make_rdm1() + if kpts is None: kpts = self.kpts if s is None: s = self.get_ovlp(cell) - return mulliken_meta(cell, dm, s=s, verbose=verbose, + return mulliken_meta(cell, dm, kpts, s=s, verbose=verbose, pre_orth_method=pre_orth_method) + @lib.with_doc(mulliken_meta_spin.__doc__) + def mulliken_meta_spin(self, cell=None, dm=None, kpts=None, verbose=logger.DEBUG, + pre_orth_method=PRE_ORTH_METHOD, s=None): + if cell is None: cell = self.cell + if dm is None: dm = self.make_rdm1() + if kpts is None: kpts = self.kpts + if s is None: s = self.get_ovlp(cell) + return mulliken_meta_spin(cell, dm, kpts, s=s, verbose=verbose, + pre_orth_method=pre_orth_method) + def mulliken_pop(self): raise NotImplementedError diff --git a/pyscf/pbc/scf/test/test_khf.py b/pyscf/pbc/scf/test/test_khf.py index 7d39168fee..101d6b12f7 100644 --- a/pyscf/pbc/scf/test/test_khf.py +++ b/pyscf/pbc/scf/test/test_khf.py @@ -62,13 +62,13 @@ def tearDownModule(): class KnownValues(unittest.TestCase): def test_analyze(self): - rpop, rchg = kmf.analyze()[0] # pop at gamma point + rpop, rchg = kmf.analyze()[0] upop, uchg = kumf.analyze()[0] gpop, gchg = kgmf.analyze()[0] self.assertTrue(isinstance(rpop, np.ndarray) and rpop.ndim == 1) self.assertAlmostEqual(abs(upop[0]+upop[1]-rpop).max(), 0, 7) self.assertAlmostEqual(abs(gpop[0]+gpop[1]-rpop).max(), 0, 5) - self.assertAlmostEqual(lib.fp(rpop), 1.697446, 5) + self.assertAlmostEqual(lib.fp(rpop), 1.638430, 5) def test_kpt_vs_supercell_high_cost(self): # For large n, agreement is always achieved diff --git a/pyscf/pbc/scf/test/test_mulliken_meta.py b/pyscf/pbc/scf/test/test_mulliken_meta.py new file mode 100644 index 0000000000..4117e25e0f --- /dev/null +++ b/pyscf/pbc/scf/test/test_mulliken_meta.py @@ -0,0 +1,114 @@ +#!/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.pbc import gto, scf, tools + + +class KPTvsSUPCELL_noshift(unittest.TestCase): + + scaled_center = None + + def make_cell_supcell_kpts(self, kmesh, **kwargs): + cell = gto.Cell() + cell.verbose = 4 + cell.output = '/dev/null' + cell.atom = ''' + O -1.4771410554 -0.0001473417 -0.0820628725 + H -2.1020842744 -0.6186960400 0.3275021172 + H -1.4531946703 -0.1515066183 -1.0398792446 + ''' + cell.a = np.eye(3) * 4 + cell.basis = 'sto-3g' + cell.set(**kwargs) + cell.build() + scell = tools.super_cell(cell, kmesh) + kpts = cell.make_kpts(kmesh, scaled_center=self.scaled_center) + return cell, scell, kpts + + def test_krhf(self): + kmesh = [3,1,1] + cell, scell, kpts = self.make_cell_supcell_kpts(kmesh) + + kmf = scf.KRHF(cell, kpts=kpts).density_fit() + kdm = kmf.init_guess_by_1e() + kchg = kmf.mulliken_meta(cell, kdm, kpts, verbose=0)[1] + + smf = scf.RHF(scell, kpt=kpts[0]).density_fit() + sdm = smf.init_guess_by_1e() + schg = smf.mulliken_meta(scell, sdm, verbose=0)[1][:cell.natm] + + self.assertAlmostEqual(abs(kchg-schg).max(), 0, 3) + + def test_kuhf(self): + spin0 = 2 + kmesh = [3,1,1] + spin = int(np.prod(kmesh)) * spin0 + cell, scell, kpts = self.make_cell_supcell_kpts(kmesh) + + kmf = scf.KUHF(cell, kpts=kpts).density_fit() + kdm = kmf.init_guess_by_1e() + kchg = kmf.mulliken_meta(cell, kdm, kpts, verbose=0)[1] + kchg_spin = kmf.mulliken_meta_spin(cell, kdm, kpts, verbose=0)[1] + + smf = scf.UHF(scell, kpt=kpts[0]).density_fit() + sdm = smf.init_guess_by_1e() + schg = smf.mulliken_meta(scell, sdm, verbose=0)[1][:cell.natm] + schg_spin = smf.mulliken_meta_spin(scell, sdm, verbose=0)[1][:cell.natm] + + self.assertAlmostEqual(abs(kchg-schg).max(), 0, 3) + self.assertAlmostEqual(abs(kchg_spin-schg_spin).max(), 0, 3) + + def test_kghf(self): + spin0 = 2 + kmesh = [3,1,1] + spin = int(np.prod(kmesh)) * spin0 + cell, scell, kpts = self.make_cell_supcell_kpts(kmesh) + + kmf = scf.KGHF(cell, kpts=kpts).density_fit() + kdm = kmf.init_guess_by_1e() + kchg = kmf.mulliken_meta(cell, kdm, kpts, verbose=0)[1] + + smf = scf.GHF(scell, kpt=kpts[0]).density_fit() + sdm = smf.init_guess_by_1e() + schg = smf.mulliken_meta(scell, sdm, verbose=0)[1][:cell.natm] + + self.assertAlmostEqual(abs(kchg-schg).max(), 0, 3) + + def test_krohf(self): + spin0 = 2 + kmesh = [3,1,1] + spin = int(np.prod(kmesh)) * spin0 + cell, scell, kpts = self.make_cell_supcell_kpts(kmesh) + + kmf = scf.KROHF(cell, kpts=kpts).density_fit() + kdm = kmf.init_guess_by_1e() + kchg = kmf.mulliken_meta(cell, kdm, kpts, verbose=0)[1] + + smf = scf.ROHF(scell, kpt=kpts[0]).density_fit() + sdm = smf.init_guess_by_1e() + schg = smf.mulliken_meta(scell, sdm, verbose=0)[1][:cell.natm] + + self.assertAlmostEqual(abs(kchg-schg).max(), 0, 3) + +class KPTvsSUPCELL_shift(KPTvsSUPCELL_noshift): + + scaled_center = [0.3765, 0.7729, 0.1692] + +if __name__ == "__main__": + print("Full Tests for PBC SCF mulliken_meta charges") + unittest.main() diff --git a/pyscf/pbc/scf/test/test_rohf.py b/pyscf/pbc/scf/test/test_rohf.py index a5c43592c5..4daa93e481 100644 --- a/pyscf/pbc/scf/test/test_rohf.py +++ b/pyscf/pbc/scf/test/test_rohf.py @@ -117,9 +117,9 @@ def test_spin_square(self): def test_analyze(self): pop, chg = kmf.analyze()[0] - self.assertAlmostEqual(lib.fp(pop), 1.1120443320325235, 7) + self.assertAlmostEqual(lib.fp(pop), 1.1514919154737624, 7) self.assertAlmostEqual(sum(chg), 0, 7) - self.assertAlmostEqual(lib.fp(chg), 0.002887875601340767, 7) + self.assertAlmostEqual(lib.fp(chg), -0.04683923436982078, 7) def test_small_system(self): # issue #686 diff --git a/pyscf/pbc/scf/uhf.py b/pyscf/pbc/scf/uhf.py index ea48c75efb..6d75689398 100644 --- a/pyscf/pbc/scf/uhf.py +++ b/pyscf/pbc/scf/uhf.py @@ -127,6 +127,7 @@ class UHF(pbchf.SCF): analyze = mol_uhf.UHF.analyze mulliken_pop = mol_uhf.UHF.mulliken_pop mulliken_meta = mol_uhf.UHF.mulliken_meta + mulliken_meta_spin = mol_uhf.UHF.mulliken_meta_spin canonicalize = mol_uhf.UHF.canonicalize spin_square = mol_uhf.UHF.spin_square stability = mol_uhf.UHF.stability diff --git a/pyscf/pbc/tools/k2gamma.py b/pyscf/pbc/tools/k2gamma.py index 46cd3ed85f..0a24d5069d 100644 --- a/pyscf/pbc/tools/k2gamma.py +++ b/pyscf/pbc/tools/k2gamma.py @@ -257,15 +257,18 @@ def transform(mo_energy, mo_coeff, mo_occ): return mf -def to_supercell_ao_integrals(cell, kpts, ao_ints): +def to_supercell_ao_integrals(cell, kpts, ao_ints, kmesh=None, force_real=True): '''Transform from the unitcell k-point AO integrals to the supercell gamma-point AO integrals. ''' - scell, phase = get_phase(cell, kpts) + scell, phase = get_phase(cell, kpts, kmesh=kmesh) NR, Nk = phase.shape nao = cell.nao scell_ints = np.einsum('Rk,kij,Sk->RiSj', phase, ao_ints, phase.conj()) - return scell_ints.reshape(NR*nao,NR*nao).real + if force_real: + return scell_ints.reshape(NR*nao,NR*nao).real + else: + return scell_ints.reshape(NR*nao,NR*nao) def to_supercell_mo_integrals(kmf, mo_ints): diff --git a/pyscf/pbc/tools/test/test_k2gamma.py b/pyscf/pbc/tools/test/test_k2gamma.py index fa4f04c98c..c2cac7590f 100644 --- a/pyscf/pbc/tools/test/test_k2gamma.py +++ b/pyscf/pbc/tools/test/test_k2gamma.py @@ -67,8 +67,8 @@ def test_k2gamma(self): mf.kernel() popa, popb = mf.mulliken_meta()[0] - self.assertAlmostEqual(lib.fp(popa), 1.5403023058, 7) - self.assertAlmostEqual(lib.fp(popb), 1.5403023058, 7) + self.assertAlmostEqual(lib.fp(popa), 1.2700920989, 7) + self.assertAlmostEqual(lib.fp(popb), 1.2700920989, 7) popa, popb = k2gamma.k2gamma(mf).mulliken_meta()[0] self.assertAlmostEqual(lib.fp(popa), 0.8007278745, 7) diff --git a/pyscf/scf/ghf.py b/pyscf/scf/ghf.py index fe37e7a630..fd8dc368bb 100644 --- a/pyscf/scf/ghf.py +++ b/pyscf/scf/ghf.py @@ -316,7 +316,7 @@ def mulliken_pop(mol, dm, s=None, verbose=logger.DEBUG): dmb = dm[nao:,nao:] if s is not None: assert (s.size == nao**2 or numpy.allclose(s[:nao,:nao], s[nao:,nao:])) - s = s[:nao,:nao] + s = lib.asarray(s[:nao,:nao], order='C') return uhf.mulliken_pop(mol, (dma,dmb), s, verbose) def mulliken_meta(mol, dm_ao, verbose=logger.DEBUG, @@ -328,7 +328,7 @@ def mulliken_meta(mol, dm_ao, verbose=logger.DEBUG, dmb = dm_ao[nao:,nao:] if s is not None: assert (s.size == nao**2 or numpy.allclose(s[:nao,:nao], s[nao:,nao:])) - s = s[:nao,:nao] + s = lib.asarray(s[:nao,:nao], order='C') return uhf.mulliken_meta(mol, (dma,dmb), verbose, pre_orth_method, s) def det_ovlp(mo1, mo2, occ1, occ2, ovlp):