[Work in progress] LNOCC #115

hongzhouye · 2025-02-25T20:22:43Z

This PR adds the molecular and periodic LNO-CCSD/CCSD(T). The LNO method was originally introduced by Kállay and co-workers:

J. Chem. Phys. 139, 094105 (2013) (link)
J. Chem. Phys. 135, 104111 (2011) (link)

And this PR is particularly based on a recent extension:

J. Chem. Theory Comput. 2024, 20, 20, 8948–8959 (link)

for spin-restricted molecular & periodic LNO-CC and an unpublished work for spin-unrestricted references (@fishjojo, @MoleOrbitalHybridAnalyst).

The PR is still work in progress with the following to-do's:

add examples and more unit tests (me and @tberkel)
make the C library pyscf-forge/pyscf/lib/lno compatible with regular pyscf (@fishjojo)
remove redundant functions in lno/lnoccsd.py for precomputation (perhaps need to slightly refactor pyscf's CCSD code) (@fishjojo )
add support for spin-unrestricted references (@fishjojo, @MoleOrbitalHybridAnalyst)

tberkel · 2025-04-25T21:14:39Z

I'm working on adding examples, but frozen seems broken with unrestricted LNO. @fishjojo or @MoleOrbitalHybridAnalyst have you tested frozen orbitals? Skimming the code suggests that you did, but I get a correlation energy of 0 when I freeze orbitals. Example is below. Works fine when frozen = 0.

from pyscf import gto, scf, mp, lo
from pyscf.data import elements
from pyscf import lno

mol = gto.Mole()
mol.atom = '''
    O   -1.485163346097   -0.114724564047    0.000000000000
    H   -1.868415346097    0.762298435953    0.000000000000
    H   -0.533833346097    0.040507435953    0.000000000000
    O    1.416468653903    0.111264435953    0.000000000000
    H    1.746241653903   -0.373945564047   -0.758561000000
    H    1.746241653903   -0.373945564047    0.758561000000
'''
mol.basis = 'ccpvdz'
mol.verbose = 7
mol.spin = 0
mol.build()

mf = scf.UHF(mol).density_fit()
mf.kernel()

frozen = elements.chemcore(mol)
# frozen = 0

mymp = mp.MP2(mf, frozen=frozen)
mymp.kernel(with_t2=False)

thresh = [1e-4] * 2

nocca, noccb = mymp.get_nocc()
moidxa, moidxb = mymp.get_frozen_mask()
mo_coeff = mf.mo_coeff[0][:, moidxa], mf.mo_coeff[1][:, moidxb]
orbocca = mo_coeff[0][:, :nocca]
lo_coeff_a = lo.PipekMezey(mf.mol, orbocca).kernel()
orboccb = mo_coeff[1][:, :noccb]
lo_coeff_b = lo.PipekMezey(mf.mol, orboccb).kernel()
lo_coeff = [lo_coeff_a, lo_coeff_b]

frag_lolist_a = [[[i], []] for i in range(nocca)]
frag_lolist_b = [[[], [i]] for i in range(noccb)]
frag_lolist = frag_lolist_a + frag_lolist_b

mlcc = lno.ULNOCCSD(mf, lo_coeff, frag_lolist, frozen=frozen,
                    lno_thresh=thresh)
mlcc.kernel()
# E(LNOMP2) = -152.062490646908  E_corr = 0
# E(LNOCCSD) = -152.062490646908  E_corr = 0

tberkel · 2025-05-06T13:13:54Z

I fixed the ULNO frozen orbitals problem (help from @hczhai) and cleaned up the unittest for ULNO.

tberkel · 2025-05-06T13:18:33Z

Another question for @fishjojo and/or @MoleOrbitalHybridAnalyst - I notice that the results for ULNO and (R)LNO are different for the same thresholds, even when the reference doesn't break spin symmetry. Did you observe the same? I guess the issue is that the orbital-specific DMs treat up and down differently (e.g., commonly a different number of up and down orbitals are included in the active space). I wonder if there's a way to spin average, so that the results agree in the limit of no spin symmetry breaking.

fishjojo · 2025-05-06T21:24:16Z

@tberkel, thanks for the fix for frozen orbitals. The ULNO results match those of RLNO from the old implementation. I notice there are some changes to the RLNO code for determining the LNOs, and I’m not sure the default is still the original ‘ie’ flavor for building the LNOs. Currently, ULNO only works for ‘ie’, as there is difficulty matching the RLNO results using other flavors. Maybe Hongzhou could provide some insights of this. Timothy Berkelbach ***@***.***>于2025年5月6日周二上午6:18写道：

…

*tberkel* left a comment (pyscf/pyscf-forge#115) <#115 (comment)> Another question for @fishjojo <https://github.com/fishjojo> and/or @MoleOrbitalHybridAnalyst <https://github.com/MoleOrbitalHybridAnalyst> - I notice that the results for ULNO and (R)LNO are different for the same thresholds, even when the reference doesn't break spin symmetry. Did you observe the same? I guess the issue is that the orbital-specific DMs treat up and down differently (e.g., commonly a different number of up and down orbitals are included in the active space). I wonder if there's a way to spin average, so that the results agree in the limit of no spin symmetry breaking. — Reply to this email directly, view it on GitHub <#115 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACCZ6FRLAEPZ6D7G2VSY74L25CZD7AVCNFSM6AAAAABX3R2U7SVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDQNJUGU2TKMBSGU> . You are receiving this because you were mentioned.Message ID: ***@***.***>

hongzhouye · 2025-05-08T01:54:09Z

@fishjojo @tberkel Could it be related to this line? I had updated the local MP2 DM definition at some point to better align with the canonical MP2 DM definition, which is twice the value of what the earlier version produced. I’m okay with defaulting to either convention.

- update k2s_scf to include hcore from kpts calculation

…maginary parts in `k2s_iao`.

…n results with lno.py

Update to unrestricted lno code

Improve k2s utilities

…sh for large molecules, I've corrected the contractions to fix this.

Bug Fix: Density Matrix Calculation for ulno

tberkel · 2025-11-17T18:17:39Z

I think this is ready to merge and meets the bar for forge. We have some more small updates, but they'll be easier to add in future PRs once this is merged.

MatthewRHermes · 2025-11-17T18:30:27Z

pyscf/lno/test/test_lnoccsd.py

+        orbocc = mf.mo_coeff[:,frozen:np.count_nonzero(mf.mo_occ)]
+        mlo = lo.PipekMezey(mol, orbocc)
+        lo_coeff = mlo.kernel()
+        while True: # always performing jacobi sweep to avoid trapping in local minimum/saddle point


Could you implement this as a for loop with a maxiter? It should behave for now but this is scary to me.

MatthewRHermes · 2025-11-17T18:32:33Z

pyscf/lno/test/test_lnoccsd.py

+            lo_coeff1 = mlo.stability_jacobi()[1]
+            if lo_coeff1 is lo_coeff:


The discarded first return argument literally should evaluate to the same value as the conditional. It's more readable and "sane" if you just use that, since who knows how stability_jacobi may be refactored in the future?

MatthewRHermes · 2025-11-17T18:34:06Z

pyscf/lno/test/test_ulnoccsd.py

+            while True: # always performing jacobi sweep to avoid trapping in local minimum/saddle point
+                lo_coeff1_s = mlo.stability_jacobi()[1]
+                if lo_coeff1_s is lo_coeff_s:


Same comments as in test_lnoccsd.py

MatthewRHermes · 2025-11-17T18:43:38Z

pyscf/lno/lno.py

+    ex = ex.real
+    ess = ed*0.5 + ex
+    eos = ed*0.5
+    return lib.tag_array(ess+eos, spin_comp=np.array((ess, eos)))


Aren't we trying to move away from lib.tag_array?

MatthewRHermes · 2025-11-17T18:45:20Z

pyscf/lno/lno.py

+    # def mo_splitter(self, kind='mask'):
+    #     r''' Return index arrays that split MOs into
+    #         - frozen occupieds
+    #         - active occupieds
+    #         - active virtuals
+    #         - frozen virtuals
+    #
+    #     Args:
+    #         kind (str):
+    #             'mask'  : return masks each of length nmo
+    #             'index' : return index arrays
+    #             'idx'   : same as 'index'
+    #     '''
+    #     maskact = self.get_frozen_mask()
+    #     maskocc = self.mo_occ > 1e-10
+    #     return mo_splitter(maskact, maskocc, kind=kind)
+    #
+    # def split_mo_coeff(self):
+    #     r''' Return the four components of MOs specified in :func:`mo_splitter`
+    #     '''
+    #     mo = self.mo_coeff
+    #     masks = self.mo_splitter()
+    #     return [mo[:,m] for m in masks]
+    #
+    # def split_mo_energy(self):
+    #     moe = self.mo_energy
+    #     masks = self.mo_splitter()
+    #     return [moe[m] for m in masks]


cleanup on aisle 3

MatthewRHermes · 2025-11-17T18:59:30Z

pyscf/lno/make_lno_rdm1.py

It's probably a good idea to have unittest coverage for a file like this, in which many functions which are different variations on a single theme are implemented. That way, the next time numpy or scipy changes something and breaks all of these functions simultaneously, you immediately get a checklist of all the places you need to update whatever it is.

MatthewRHermes · 2025-11-17T19:03:03Z

pyscf/pbc/lno/test/test_klnoccsd.py

+        del cls.cell, cls.kmf, cls.frozen
+        del cls.scell, cls.smf
+
+    # def test_lno_pm_by_thresh(self):


Use @unittest.skip ('reason') if you intend to enable this test in the future; or delete it if not.

MatthewRHermes · 2025-11-17T19:04:18Z

pyscf/pbc/lno/klnoccsd.py

+    # def contract1(LbraR, LbraI, LketR, LketI):
+    #     out_shape = (LbraR.shape[1], LketR.shape[1])
+    #     out = np.zeros(out_shape, dtype=REAL)
+    #     for q1,q2 in enumerate(k2sdf.qconserv):
+    #         i0,i1 = k2sdf.get_auxslice(q1)
+    #         j0,j1 = k2sdf.get_auxslice(q2)
+    #         zdotNNtoR(LbraR[i0:i1].T, LbraI[i0:i1].T, LketR[j0:j1], LketI[j0:j1], 1, out, 1)
+    #     return out
+    #
+    # eris.oooo[:] = contract1(LooR, LooI, LooR, LooI).reshape(nocc,nocc,nocc,nocc)
+    # eris.ovoo[:] = contract1(LovR, LovI, LooR, LooI).reshape(nocc,nvir,nocc,nocc)
+    # eris.oovv[:] = lib.unpack_tril(
+    #                contract1(LooR, LooI, LvvR, LvvI)).reshape(nocc,nocc,nvir,nvir)
+    # eris.ovvo[:] = contract1(LovR, LovI, LvoR, LvoI).reshape(nocc,nvir,nvir,nocc)
+    # eris.ovov[:] = contract1(LovR, LovI, LovR, LovI).reshape(nocc,nvir,nocc,nvir)
+    # eris.ovvv[:] = contract1(LovR, LovI, LvvR, LvvI).reshape(nocc,nvir,nvir_pair)
+    # eris.vvvv[:] = contract1(LvvR, LvvI, LvvR, LvvI).reshape(nvir_pair,nvir_pair)


MatthewRHermes · 2025-11-17T19:05:06Z

pyscf/pbc/lno/make_lno_rdm1.py

same comment as for the OBC version

MatthewRHermes

Could we get a quick and dirty example input script or two?

jeanwsr · 2025-11-18T13:13:19Z

pyscf/lno/lnoccsd.py

+
+    @property
+    def e_corr_pt2_scs(self):
+        e_corr = 0.333*self.e_corr_pt2_ss + 1.2*self.e_corr_pt2_os


maybe the ss, os factor can be configurable arguments, like the DF SCS-MP2 does

https://github.com/pyscf/pyscf/blob/e155c487f6f44dc249089c7d6aafded58ce53e0e/pyscf/mp/dfmp2_native.py#L138

https://github.com/pyscf/pyscf/blob/e155c487f6f44dc249089c7d6aafded58ce53e0e/pyscf/mp/dfmp2_native.py#L262

hongzhouye and others added 6 commits February 25, 2025 15:02

add spin restricted LNOCC

379c552

fixing style

49a42fc

Merge remote-tracking branch 'hongzhou/lnocc' into lnocc

554d072

add ulnoccsd

adddbba

update C backend

a30f1cf

Merge remote-tracking branch 'upstream/master' into lnocc

41ce90e

Fix ULNO for frozen orbs, rewrite test

b70d77a

bug fix for outcore lno

138fab9

yangjunjie0320 and others added 14 commits May 22, 2025 22:44

- rename k2sdf to _k2sdf in MODIFIED_K2SCCSD to avoid warning

39f507f

- update k2s_scf to include hcore from kpts calculation

- adds stream object inheritance to K2SDF and a warning for large i…

6c53521

…maginary parts in `k2s_iao`.

- addk2sdf to _keys in MODIFIED_K2SCCSD

ef9c144

Removed rs_density_fit in k2s_scf for broader DF compatibility

40432f2

Add convergence check to impurity_solve

f86190d

Fix logic for iao coeff real part

5c176fb

Optimize ulno.py for memory and include multiplicative factor to alig…

27ce793

…n results with lno.py

Update ulno unit test using new dm convention

348659c

Merge pull request #2 from afarahva/lnocc

ea020a2

Update to unrestricted lno code

Merge pull request #1 from yangjunjie0320/lnocc

75e0ee0

Improve k2s utilities

Bug Fix: Density Matrix Calculation for ulno was bugged and would cra…

d56122d

…sh for large molecules, I've corrected the contractions to fix this.

removed print statements

b815a61

Merge pull request #3 from afarahva/lnocc

2ff1359

Bug Fix: Density Matrix Calculation for ulno

Merge branch 'master' into lnocc

3f5c7c8

tberkel requested a review from MatthewRHermes November 17, 2025 18:16

MatthewRHermes reviewed Nov 17, 2025

View reviewed changes

pyscf/pbc/lno/make_lno_rdm1.py

Copy link

Collaborator

MatthewRHermes Nov 17, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same comment as for the OBC version

MatthewRHermes reviewed Nov 17, 2025

View reviewed changes

jeanwsr reviewed Nov 18, 2025

View reviewed changes

		lo_coeff1 = mlo.stability_jacobi()[1]
		if lo_coeff1 is lo_coeff:

[Work in progress] LNOCC #115

Are you sure you want to change the base?

[Work in progress] LNOCC #115

Uh oh!

Conversation

hongzhouye commented Feb 25, 2025 • edited by fishjojo Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

tberkel commented Apr 25, 2025

Uh oh!

tberkel commented May 6, 2025

Uh oh!

tberkel commented May 6, 2025

Uh oh!

fishjojo commented May 6, 2025 via email

Uh oh!

hongzhouye commented May 8, 2025

Uh oh!

tberkel commented Nov 17, 2025

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

MatthewRHermes left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

7 participants

hongzhouye commented Feb 25, 2025 •

edited by fishjojo

Loading