• Home
  • Features
  • Pricing
  • Docs
  • Announcements
  • Sign In

adc-connect / adcc / 22232189935

20 Feb 2026 04:30PM UTC coverage: 74.509% (+0.3%) from 74.248%
22232189935

Pull #202

github

web-flow
Merge b55489e7f into b49d31427
Pull Request #202: <S^2> calculation for UHF

1292 of 2028 branches covered (63.71%)

Branch coverage included in aggregate %.

166 of 205 new or added lines in 8 files covered. (80.98%)

1 existing line in 1 file now uncovered.

8076 of 10545 relevant lines covered (76.59%)

167693.6 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

87.84
/adcc/adc_pp/state_diffdm_2p.py
1
#!/usr/bin/env python3
2
## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab
3
## ---------------------------------------------------------------------
4
##
5
## Copyright (C) 2026 by the adcc authors
6
##
7
## This file is part of adcc.
8
##
9
## adcc is free software: you can redistribute it and/or modify
10
## it under the terms of the GNU General Public License as published
11
## by the Free Software Foundation, either version 3 of the License, or
12
## (at your option) any later version.
13
##
14
## adcc is distributed in the hope that it will be useful,
15
## but WITHOUT ANY WARRANTY; without even the implied warranty of
16
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
17
## GNU General Public License for more details.
18
##
19
## You should have received a copy of the GNU General Public License
20
## along with adcc. If not, see <http://www.gnu.org/licenses/>.
21
##
22
## ---------------------------------------------------------------------
23
from adcc import block as b
2✔
24
from adcc.LazyMp import LazyMp
2✔
25
from adcc.AdcMethod import AdcMethod
2✔
26
from adcc.functions import einsum, zeros_like
2✔
27
from adcc.Intermediates import Intermediates
2✔
28
from adcc.AmplitudeVector import AmplitudeVector
2✔
29
from adcc.TwoParticleDensity import TwoParticleDensity
2✔
30
from adcc.NParticleOperator import OperatorSymmetry
2✔
31

32
from .util import check_doubles_amplitudes, check_singles_amplitudes
2✔
33

34

35
def diffdm_adc0_2p(mp, amplitude, intermediates):
2✔
36
    check_singles_amplitudes([b.o, b.v], amplitude)
2✔
37
    u1 = amplitude.ph
2✔
38

39
    hf = mp.reference_state
2✔
40
    d_oo = zeros_like(hf.foo)
2✔
41
    d_oo.set_mask("ii", 1)
2✔
42

43
    dm = TwoParticleDensity(mp, symmetry=OperatorSymmetry.HERMITIAN)
2✔
44

45
    # ADC(1) diffdm
46
    p1_oo = -einsum("ia,la->il", u1, u1).evaluate()
2✔
47
    p1_vv = einsum("ka,kb->ab", u1, u1).evaluate()
2✔
48

49
    dm.oooo = (
2✔
50
        # N^4: O^4 / N^4: O^4
51
        - 4.0 * einsum("il,jk->ijkl", p1_oo, d_oo)
52
    ).antisymmetrise(0, 1).antisymmetrise(2, 3)
53
    dm.ovov = (
2✔
54
        # N^4: O^2V^2 / N^4: O^2V^2
55
        - 1.0 * einsum("ja,ib->iajb", u1, u1)
56
        # N^4: O^2V^2 / N^4: O^2V^2
57
        + 1.0 * einsum("ab,ij->iajb", p1_vv, d_oo)
58
    )
59
    return dm
2✔
60

61

62
def diffdm_adc1_2p(mp, amplitude, intermediates):
2✔
63
    dm = diffdm_adc0_2p(mp, amplitude, intermediates)  # Get ADC(0) result
2✔
64
    u1 = amplitude.ph
2✔
65

66
    hf = mp.reference_state
2✔
67
    d_oo = zeros_like(hf.foo)
2✔
68
    d_oo.set_mask("ii", 1)
2✔
69

70
    t2 = mp.t2(b.oovv)
2✔
71
    # TODO move to intermediates!
72
    # ADC(1) diffdm
73
    p1_oo = -einsum("ia,la->il", u1, u1).evaluate()
2✔
74

75
    # ADC(2) ISR intermediate (TODO Move to intermediates)
76
    ru1 = einsum("ijab,jb->ia", t2, u1).evaluate()
2✔
77
    # new ones
78
    ru1_ooov = einsum("kc,ijac->ijka", u1, t2).evaluate()
2✔
79

80
    dm.oovv += (
2✔
81
        # N^4: O^2V^2 / N^4: O^2V^2
82
        + 4.0 * (
83
            einsum("ib,ja->ijab", ru1, u1)
84
        ).antisymmetrise(0, 1).antisymmetrise(2, 3)
85
        # N^5: O^3V^2 / N^4: O^2V^2
86
        + 2.0 * (
87
            einsum("ijka,kb->ijab", ru1_ooov, u1)
88
        ).antisymmetrise(2, 3)
89
        # N^5: O^3V^2 / N^4: O^2V^2
90
        - 2.0 * (
91
            einsum("jk,ikab->ijab", p1_oo, t2)
92
        ).antisymmetrise(0, 1)
93
    )
94
    return dm
2✔
95

96

97
def diffdm_adc2_2p(mp, amplitude, intermediates):
2✔
98
    dm = diffdm_adc1_2p(mp, amplitude, intermediates)  # Get ADC(1) result
2✔
99
    check_doubles_amplitudes([b.o, b.o, b.v, b.v], amplitude)
2✔
100
    u1, u2 = amplitude.ph, amplitude.pphh
2✔
101
    hf = mp.reference_state
2✔
102
    d_oo = zeros_like(hf.foo)
2✔
103
    d_oo.set_mask("ii", 1)
2✔
104

105
    t2 = mp.t2(b.oovv)
2✔
106
    td2 = mp.td2(b.oovv)
2✔
107
    p0 = mp.mp2_diffdm
2✔
108

109
    dm.oooo += (
2✔
110
        # N^6: O^4V^2 / N^4: O^2V^2
111
        + 2.0 * einsum("ijab,klab->ijkl", u2, u2)
112
        # N^6: O^5V^1 / N^4: O^2V^2
113
        - 1.0 * einsum("klmc,ijmc->ijkl", einsum("mb,klbc->klmc", u1, t2),
114
                       einsum("ma,ijac->ijmc", u1, t2))
115
        + (
116
            # N^5: O^3V^2 / N^4: O^2V^2
117
            - 8.0 * einsum("ik,jl->ijkl", einsum("imab,kmab->ik", u2, u2), d_oo)
118
            # N^4: O^4 / N^4: O^4
119
            + 4.0 * einsum("il,jk->ijkl", einsum("ia,la->il", u1, u1), p0.oo)
120
            # N^6: O^5V^1 / N^4: O^2V^2
121
            + 4.0 * einsum("jlmc,ikmc->ijkl", einsum("lb,jmbc->jlmc", u1, t2),
122
                           einsum("ia,kmac->ikmc", u1, t2))
123
            # N^4: O^2V^2 / N^4: O^2V^2
124
            + 4.0 * einsum("il,jk->ijkl",
125
                           einsum("lc,ic->il", einsum("nb,lnbc->lc", u1, t2),
126
                                  einsum("ma,imac->ic", u1, t2)), d_oo)
127
            # N^5: O^3V^2 / N^4: O^2V^2
128
            + 4.0 * einsum("ik,jl->ijkl", einsum("kmnc,imnc->ik",
129
                                                 einsum("nb,kmbc->kmnc", u1, t2),
130
                                                 einsum("na,imac->imnc", u1, t2)),
131
                           d_oo)
132
            # N^5: O^3V^2 / N^4: O^2V^2
133
            + 2.0 * einsum("ik,jl->ijkl",
134
                           einsum("inbc,knbc->ik",
135
                                  einsum("mn,imbc->inbc",
136
                                         einsum("ia,la->il", u1, u1), t2), t2),
137
                           d_oo)
138
        ).antisymmetrise(0, 1).antisymmetrise(2, 3)
139
        + (
140
            # N^5: O^3V^2 / N^4: O^2V^2
141
            + 4.0 * einsum("iklb,jb->ijkl",
142
                           einsum("ic,klbc->iklb", einsum("ma,imac->ic", u1, t2),
143
                                  t2), u1)
144
            # N^6: O^4V^2 / N^4: O^2V^2
145
            + 2.0 * einsum("jklm,im->ijkl",
146
                           einsum("jmbc,klbc->jklm", t2, t2),
147
                           einsum("ia,la->il", u1, u1))
148
        ).antisymmetrise(0, 1).symmetrise([(0, 2), (1, 3)])
149
        + (
150
            # N^4: O^4 / N^4: O^4
151
            + 4.0 * einsum("il,jk->ijkl", einsum("im,lm->il",
152
                                                 einsum("ia,la->il", u1, u1),
153
                                                 p0.oo), d_oo)
154
            # N^4: O^2V^2 / N^4: O^2V^2
155
            + 4.0 * einsum("ik,jl->ijkl",
156
                           einsum("kb,ib->ik",
157
                                  einsum("mc,kmbc->kb",
158
                                         einsum("na,mnac->mc", u1, t2), t2), u1),
159
                           d_oo)
160
        ).antisymmetrise(0, 1).antisymmetrise(2, 3).symmetrise([(0, 2), (1, 3)])
161
    )
162
    dm.ooov += (
2✔
163
        # N^5: O^3V^2 / N^4: O^2V^2
164
        - 2.0 * einsum("kb,ijab->ijka", u1, u2)
165
        # N^6: O^4V^2 / N^4: O^2V^2
166
        + 1.0 * einsum("ijkl,la->ijka", einsum("klbc,ijbc->ijkl", u2, t2), u1)
167
        # N^5: O^3V^2 / N^4: O^2V^2
168
        + 2.0 * einsum("kb,ijab->ijka", einsum("lc,klbc->kb", u1, u2), t2)
169
        + (
170
            # N^4: O^2V^2 / N^4: O^2V^2
171
            - 4.0 * einsum("ja,ik->ijka", einsum("lb,jlab->ja", u1, u2), d_oo)
172
            # N^5: O^3V^2 / N^4: O^2V^2
173
            + 2.0 * einsum("jk,ia->ijka", einsum("klbc,jlbc->jk", u2, t2), u1)
174
            # N^6: O^4V^2 / N^4: O^2V^2
175
            - 4.0 * einsum("jklb,ilab->ijka", einsum("jc,klbc->jklb", u1, u2), t2)
176
            # N^5: O^3V^2 / N^4: O^2V^2
177
            + 2.0 * einsum("ja,ik->ijka",
178
                           einsum("jlmb,lmab->ja",
179
                                  einsum("jc,lmbc->jlmb", u1, u2), t2), d_oo)
180
            # N^4: O^2V^2 / N^4: O^2V^2
181
            - 4.0 * einsum("ia,jk->ijka",
182
                           einsum("lb,ilab->ia",
183
                                  einsum("mc,lmbc->lb", u1, u2), t2), d_oo)
184
            # N^5: O^3V^2 / N^4: O^2V^2
185
            + 2.0 * einsum("ia,jk->ijka",
186
                           einsum("il,la->ia",
187
                                  einsum("lmbc,imbc->il", u2, t2), u1), d_oo)
188
            # N^4: O^3V^1 / N^4: O^3V^1
189
            + 2.0 * einsum("jk,ia->ijka", einsum("kb,jb->jk", u1, p0.ov), u1)
190
            # N^4: O^3V^1 / N^4: O^3V^1
191
            + 2.0 * einsum("jk,ia->ijka", einsum("jb,kb->jk", u1, u1), p0.ov)
192
            # N^4: O^3V^1 / N^4: O^3V^1
193
            + 2.0 * einsum("ia,jk->ijka",
194
                           einsum("il,la->ia",
195
                                  einsum("lb,ib->il", u1, p0.ov), u1), d_oo)
196
            # N^4: O^3V^1 / N^4: O^3V^1
197
            + 2.0 * einsum("ia,jk->ijka",
198
                           einsum("il,la->ia", einsum("ia,la->il", u1, u1),
199
                                  p0.ov), d_oo)
200
        ).antisymmetrise(0, 1)
201
    )
202
    dm.oovv += (
2✔
203
        # N^4: O^2V^2 / N^4: O^2V^2
204
        + 4.0 * einsum("ib,ja->ijab",
205
                       einsum("kc,ikbc->ib", u1, td2),
206
                       u1).antisymmetrise(0, 1).antisymmetrise(2, 3)
207
        # N^5: O^3V^2 / N^4: O^2V^2
208
        + 2.0 * einsum("ijka,kb->ijab",
209
                       einsum("kc,ijac->ijka", u1, td2), u1).antisymmetrise(2, 3)
210
        # N^5: O^3V^2 / N^4: O^2V^2
211
        + 2.0 * einsum("jk,ikab->ijab", einsum("jc,kc->jk", u1, u1), td2)
212
    )
213
    dm.ovov += (
2✔
214
        # N^6: O^3V^3 / N^4: O^2V^2
215
        - 4.0 * einsum("jkac,ikbc->iajb", u2, u2)
216
        # N^5: O^2V^3 / N^4: O^2V^2
217
        + 2.0 * einsum("ab,ij->iajb", einsum("klac,klbc->ab", u2, u2), d_oo)
218
        # N^6: O^4V^2 / N^4: O^2V^2
219
        + 1.0 * einsum("iklb,jkla->iajb", einsum("ld,ikbd->iklb", u1, t2),
220
                       einsum("lc,jkac->jkla", u1, t2))
221
        # N^4: O^2V^2 / N^4: O^2V^2
222
        + 1.0 * einsum("ab,ij->iajb", einsum("ka,kb->ab", u1, u1), p0.oo)
223
        # N^4: O^2V^2 / N^4: O^2V^2
224
        - 1.0 * einsum("ij,ab->iajb", einsum("ic,jc->ij", u1, u1), p0.vv)
225
        # N^6: O^3V^3 / N^4: O^2V^2
226
        + 1.0 * einsum("jlac,ilbc->iajb",
227
                       einsum("kl,jkac->jlac", einsum("kd,ld->kl", u1, u1), t2), t2)
228
        # N^6: O^4V^2 / N^4: O^2V^2
229
        + 0.5 * einsum("ijla,lb->iajb",
230
                       einsum("ijkl,ka->ijla",
231
                              einsum("ikcd,jlcd->ijkl", t2, t2), u1), u1)
232
        # N^6: O^4V^2 / N^4: O^2V^2
233
        + 0.5 * einsum("jklb,ikla->iajb",
234
                       einsum("jd,klbd->jklb", u1, t2),
235
                       einsum("ic,klac->ikla", u1, t2))
236
        # N^4: O^2V^2 / N^4: O^2V^2
237
        - 1.0 * einsum("ib,ja->iajb",
238
                       einsum("kd,ikbd->ib", u1, t2), einsum("lc,jlac->ja", u1, t2))
239
        # N^5: O^2V^3 / N^4: O^2V^2
240
        + 1.0 * einsum("ab,ij->iajb",
241
                       einsum("lmad,lmbd->ab",
242
                              einsum("km,klad->lmad",
243
                                     einsum("kc,mc->km", u1, u1), t2), t2), d_oo)
244
        # N^4: O^2V^2 / N^4: O^2V^2
245
        - 1.0 * einsum("ab,ij->iajb",
246
                       einsum("lb,la->ab",
247
                              einsum("kd,klbd->lb", u1, t2),
248
                              einsum("mc,lmac->la", u1, t2)), d_oo)
249
        # N^5: O^3V^2 / N^4: O^2V^2
250
        - 0.5 * einsum("ab,ij->iajb",
251
                       einsum("klmb,klma->ab",
252
                              einsum("md,klbd->klmb", u1, t2),
253
                              einsum("mc,klac->klma", u1, t2)), d_oo)
254
        + (
255
            # N^4: O^2V^2 / N^4: O^2V^2
256
            + 1.0 * einsum("ib,ja->iajb", einsum("ic,bc->ib", u1, p0.vv), u1)
257
            # N^4: O^2V^2 / N^4: O^2V^2
258
            - 1.0 * einsum("ib,ja->iajb", einsum("kb,ik->ib", u1, p0.oo), u1)
259
            # N^5: O^3V^2 / N^4: O^2V^2
260
            + 2.0 * einsum("ijlb,la->iajb",
261
                           einsum("jc,ilbc->ijlb",
262
                                  einsum("kd,jkcd->jc", u1, t2), t2), u1)
263
            # N^5: O^3V^2 / N^4: O^2V^2
264
            + 2.0 * einsum("kb,ijka->iajb",
265
                           einsum("ld,klbd->kb", u1, t2),
266
                           einsum("ic,jkac->ijka", u1, t2))
267
            # N^6: O^4V^2 / N^4: O^2V^2
268
            - 2.0 * einsum("ijkb,ka->iajb",
269
                           einsum("jklc,ilbc->ijkb",
270
                                  einsum("kd,jlcd->jklc", u1, t2), t2), u1)
271
            # N^6: O^4V^2 / N^4: O^2V^2
272
            - 2.0 * einsum("ijlb,la->iajb",
273
                           einsum("ijkc,klbc->ijlb",
274
                                  einsum("id,jkcd->ijkc", u1, t2), t2), u1)
275
            # N^6: O^3V^3 / N^4: O^2V^2
276
            - 2.0 * einsum("ikbc,jkac->iajb",
277
                           einsum("il,klbc->ikbc",
278
                                  einsum("id,ld->il", u1, u1), t2), t2)
279
            # N^4: O^2V^2 / N^4: O^2V^2
280
            - 1.0 * einsum("ib,ja->iajb",
281
                           einsum("kc,ikbc->ib",
282
                                  einsum("ld,klcd->kc", u1, t2), t2), u1)
283
            # N^4: O^2V^2 / N^4: O^2V^2
284
            - 1.0 * einsum("ab,ij->iajb",
285
                           einsum("kb,ka->ab",
286
                                  einsum("kc,bc->kb", u1, p0.vv), u1), d_oo)
287
            # N^4: O^2V^2 / N^4: O^2V^2
288
            - 1.0 * einsum("ab,ij->iajb",
289
                           einsum("mb,ma->ab",
290
                                  einsum("ld,lmbd->mb",
291
                                         einsum("kc,klcd->ld", u1, t2),
292
                                         t2), u1), d_oo)
293
        ).symmetrise([(0, 2), (1, 3)])
294
    )
295
    dm.ovvv += (
2✔
296
        # N^5: O^2V^3 / N^4: O^1V^3
297
        - 2.0 * einsum("ja,ijbc->iabc", u1, u2)
298
        # N^6: O^3V^3 / N^4: O^1V^3
299
        + 1.0 * einsum("ijka,jkbc->iabc", einsum("id,jkad->ijka", u1, u2), t2)
300
        # N^5: O^2V^3 / N^4: O^1V^3
301
        + 2.0 * einsum("ja,ijbc->iabc", einsum("kd,jkad->ja", u1, u2), t2)
302
        + (
303
            # N^5: O^2V^3 / N^4: O^1V^3
304
            + 2.0 * einsum("ac,ib->iabc", einsum("jkad,jkcd->ac", u2, t2), u1)
305
            # N^6: O^3V^3 / N^4: O^1V^3
306
            - 4.0 * einsum("ikab,kc->iabc", einsum("jkad,ijbd->ikab", u2, t2), u1)
307
            # N^4: O^1V^3 / N^4: O^1V^3
308
            + 2.0 * einsum("ac,ib->iabc", einsum("ja,jc->ac", u1, p0.ov), u1)
309
            # N^4: O^1V^3 / N^4: O^1V^3
310
            + 2.0 * einsum("ac,ib->iabc", einsum("ja,jc->ac", u1, u1), p0.ov)
311
        ).antisymmetrise(2, 3)
312
    )
313
    dm.vvvv += (
2✔
314
        # N^6: O^2V^4 / N^4: V^4
315
        + 2 * einsum("ijab,ijcd->abcd", u2, u2)
316
        # N^6: O^2V^4 / N^4: V^4
317
        - 1 * einsum("jkab,jkcd->abcd",
318
                     einsum("ij,ikab->jkab", einsum("ie,je->ij", u1, u1), t2), t2)
319
        + (
320
            # N^4: V^4 / N^4: V^4
321
            + 4.0 * einsum("ac,bd->abcd", einsum("ia,ic->ac", u1, u1), p0.vv)
322
            # N^6: O^3V^3 / N^4: V^4
323
            + 4.0 * einsum("jabc,jd->abcd",
324
                           einsum("ijbc,ia->jabc",
325
                                  einsum("jkbe,ikce->ijbc", t2, t2), u1), u1)
326
        ).antisymmetrise(0, 1).antisymmetrise(2, 3)
327
        + (
328
            # N^5: O^1V^4 / N^4: V^4
329
            + 4.0 * einsum("ka,kbcd->abcd",
330
                           einsum("ie,ikae->ka", u1, t2),
331
                           einsum("jb,jkcd->kbcd", u1, t2))
332
            # N^6: O^3V^3 / N^4: V^4
333
            + 2.0 * einsum("ibcd,ia->abcd",
334
                           einsum("ijkb,jkcd->ibcd",
335
                                  einsum("ie,jkbe->ijkb", u1, t2), t2), u1)
336
        ).antisymmetrise(0, 1).symmetrise([(0, 2), (1, 3)])
337
    )
338
    return dm
2✔
339

340

341
# dict controlling the dispatch of the state_diffdm function
342
DISPATCH = {
2✔
343
    "adc0": diffdm_adc0_2p,
344
    "adc1": diffdm_adc1_2p,
345
    "adc2": diffdm_adc2_2p,
346
    "adc2x": diffdm_adc2_2p,      # same as ADC(2)
347
}
348

349

350
def state_diffdm_2p(method, ground_state, amplitude, intermediates=None):
2✔
351
    """
352
    Compute the two-particle difference density matrix of an excited state
353
    in the MO basis.
354

355
    Parameters
356
    ----------
357
    method : str, AdcMethod
358
        The method to use for the computation (e.g. "adc2")
359
    ground_state : LazyMp
360
        The ground state upon which the excitation was based
361
    amplitude : AmplitudeVector
362
        The amplitude vector
363
    intermediates : adcc.Intermediates
364
        Intermediates from the ADC calculation to reuse
365
    """
366
    if not isinstance(method, AdcMethod):
2!
NEW
367
        method = AdcMethod(method)
×
368
    if not isinstance(ground_state, LazyMp):
2!
NEW
369
        raise TypeError("ground_state should be a LazyMp object.")
×
370
    if not isinstance(amplitude, AmplitudeVector):
2!
NEW
371
        raise TypeError("amplitude should be an AmplitudeVector object.")
×
372
    if intermediates is None:
2!
373
        intermediates = Intermediates(ground_state)
2✔
374

375
    if method.name not in DISPATCH:
2!
NEW
376
        raise NotImplementedError("state_diffdm_2p is not implemented "
×
377
                                  f"for {method.name}.")
378
    else:
379
        ret = DISPATCH[method.name](ground_state, amplitude, intermediates)
2✔
380
        return ret.evaluate()
2✔
STATUS · Troubleshooting · Open an Issue · Sales · Support · CAREERS · ENTERPRISE · START FREE · SCHEDULE DEMO
ANNOUNCEMENTS · TWITTER · TOS & SLA · Supported CI Services · What's a CI service? · Automated Testing

© 2026 Coveralls, Inc