adc-connect/adcc | Build 25782308614 | adcc/adc_pp/state_diffdm_2p.py | Coveralls

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

## ---------------------------------------------------------------------

32

    check_doubles_amplitudes, check_singles_amplitudes, check_triples_amplitudes

33

34

35

39

43

45

46

    # one-particle ISR(0) diffdm

49

51

        # N^4: O^4 / N^4: O^4

52

        - 4.0 * einsum("il,jk->ijkl", p1_oo, d_oo)

53

    ).antisymmetrise(0, 1).antisymmetrise(2, 3)

55

        # N^4: O^2V^2 / N^4: O^2V^2

56

        - 1.0 * einsum("ja,ib->iajb", u1, u1)

57

        # N^4: O^2V^2 / N^4: O^2V^2

58

        + 1.0 * einsum("ab,ij->iajb", p1_vv, d_oo)

59

61

62

66

70

72

    # TODO move to intermediates!

73

    # one-particle ISR(0) diffdm

75

76

    # ISR(2) ISR intermediate (TODO Move to intermediates)

78

    # new ones

80

82

        # N^4: O^2V^2 / N^4: O^2V^2

83

        + 4.0 * (

84

            einsum("ib,ja->ijab", ru1, u1)

85

        ).antisymmetrise(0, 1).antisymmetrise(2, 3)

86

        # N^5: O^3V^2 / N^4: O^2V^2

87

        + 2.0 * (

88

            einsum("ijka,kb->ijab", ru1_ooov, u1)

89

        ).antisymmetrise(2, 3)

90

        # N^5: O^3V^2 / N^4: O^2V^2

91

        - 2.0 * (

92

            einsum("jk,ikab->ijab", p1_oo, t2)

93

        ).antisymmetrise(0, 1)

94

95

97

98

101

103

        # ISR(1)-d

106

110

112

            # N^5: O^3V^2 / N^4: O^2V^2

113

            - 2.0 * einsum("kb,ijab->ijka", u1, u2)

114

            # N^4: O^2V^2 / N^4: O^2V^2

115

            - 4.0 * einsum(

116

                "ja,ik->ijka", einsum("lb,jlab->ja", u1, u2), d_oo

117

            ).antisymmetrise(0, 1)

118

119

121

            # N^5: O^2V^3 / N^4: O^1V^3

122

            - 2.0 * einsum("ja,ijbc->iabc", u1, u2)

123

125

        # no doubles contribution

128

129

137

141

143

        # N^6: O^4V^2 / N^4: O^2V^2

144

        + 2.0 * einsum("ijab,klab->ijkl", u2, u2)

145

        # N^6: O^5V^1 / N^4: O^2V^2

146

        - 1.0 * einsum("klmc,ijmc->ijkl", einsum("mb,klbc->klmc", u1, t2),

147

                       einsum("ma,ijac->ijmc", u1, t2))

148

+ (

149

            # N^5: O^3V^2 / N^4: O^2V^2

150

            - 8.0 * einsum("ik,jl->ijkl", einsum("imab,kmab->ik", u2, u2), d_oo)

151

            # N^4: O^4 / N^4: O^4

152

            + 4.0 * einsum("il,jk->ijkl", einsum("ia,la->il", u1, u1), p0.oo)

153

            # N^6: O^5V^1 / N^4: O^2V^2

154

            + 4.0 * einsum("jlmc,ikmc->ijkl", einsum("lb,jmbc->jlmc", u1, t2),

155

                           einsum("ia,kmac->ikmc", u1, t2))

156

            # N^4: O^2V^2 / N^4: O^2V^2

157

            + 4.0 * einsum("il,jk->ijkl",

158

                           einsum("lc,ic->il", einsum("nb,lnbc->lc", u1, t2),

159

                                  einsum("ma,imac->ic", u1, t2)), d_oo)

160

            # N^5: O^3V^2 / N^4: O^2V^2

161

            + 4.0 * einsum("ik,jl->ijkl", einsum("kmnc,imnc->ik",

162

                                                 einsum("nb,kmbc->kmnc", u1, t2),

163

                                                 einsum("na,imac->imnc", u1, t2)),

164

                           d_oo)

165

            # N^5: O^3V^2 / N^4: O^2V^2

166

            + 2.0 * einsum("ik,jl->ijkl",

167

                           einsum("inbc,knbc->ik",

168

                                  einsum("mn,imbc->inbc",

169

                                         einsum("ia,la->il", u1, u1), t2), t2),

170

                           d_oo)

171

        ).antisymmetrise(0, 1).antisymmetrise(2, 3)

172

+ (

173

            # N^5: O^3V^2 / N^4: O^2V^2

174

            + 4.0 * einsum("iklb,jb->ijkl",

175

                           einsum("ic,klbc->iklb", einsum("ma,imac->ic", u1, t2),

176

                                  t2), u1)

177

            # N^6: O^4V^2 / N^4: O^2V^2

178

            + 2.0 * einsum("jklm,im->ijkl",

179

                           einsum("jmbc,klbc->jklm", t2, t2),

180

                           einsum("ia,la->il", u1, u1))

181

        ).antisymmetrise(0, 1).symmetrise([(0, 2), (1, 3)])

182

+ (

183

            # N^4: O^4 / N^4: O^4

184

            + 4.0 * einsum("il,jk->ijkl", einsum("im,lm->il",

185

                                                 einsum("ia,la->il", u1, u1),

186

                                                 p0.oo), d_oo)

187

            # N^4: O^2V^2 / N^4: O^2V^2

188

            + 4.0 * einsum("ik,jl->ijkl",

189

                           einsum("kb,ib->ik",

190

                                  einsum("mc,kmbc->kb",

191

                                         einsum("na,mnac->mc", u1, t2), t2), u1),

192

                           d_oo)

193

        ).antisymmetrise(0, 1).antisymmetrise(2, 3).symmetrise([(0, 2), (1, 3)])

194

196

        # N^6: O^4V^2 / N^4: O^2V^2

197

        + 1.0 * einsum("ijkl,la->ijka", einsum("klbc,ijbc->ijkl", u2, t2), u1)

198

        # N^5: O^3V^2 / N^4: O^2V^2

199

        + 2.0 * einsum("kb,ijab->ijka", einsum("lc,klbc->kb", u1, u2), t2)

200

+ (

201

            # N^5: O^3V^2 / N^4: O^2V^2

202

            + 2.0 * einsum("jk,ia->ijka", einsum("klbc,jlbc->jk", u2, t2), u1)

203

            # N^6: O^4V^2 / N^4: O^2V^2

204

            - 4.0 * einsum("jklb,ilab->ijka", einsum("jc,klbc->jklb", u1, u2), t2)

205

            # N^5: O^3V^2 / N^4: O^2V^2

206

            + 2.0 * einsum("ja,ik->ijka",

207

                           einsum("jlmb,lmab->ja",

208

                                  einsum("jc,lmbc->jlmb", u1, u2), t2), d_oo)

209

            # N^4: O^2V^2 / N^4: O^2V^2

210

            - 4.0 * einsum("ia,jk->ijka",

211

                           einsum("lb,ilab->ia",

212

                                  einsum("mc,lmbc->lb", u1, u2), t2), d_oo)

213

            # N^5: O^3V^2 / N^4: O^2V^2

214

            + 2.0 * einsum("ia,jk->ijka",

215

                           einsum("il,la->ia",

216

                                  einsum("lmbc,imbc->il", u2, t2), u1), d_oo)

217

            # N^4: O^3V^1 / N^4: O^3V^1

218

            + 2.0 * einsum("jk,ia->ijka", einsum("kb,jb->jk", u1, p0.ov), u1)

219

            # N^4: O^3V^1 / N^4: O^3V^1

220

            + 2.0 * einsum("jk,ia->ijka", einsum("jb,kb->jk", u1, u1), p0.ov)

221

            # N^4: O^3V^1 / N^4: O^3V^1

222

            + 2.0 * einsum("ia,jk->ijka",

223

                           einsum("il,la->ia",

224

                                  einsum("lb,ib->il", u1, p0.ov), u1), d_oo)

225

            # N^4: O^3V^1 / N^4: O^3V^1

226

            + 2.0 * einsum("ia,jk->ijka",

227

                           einsum("il,la->ia", einsum("ia,la->il", u1, u1),

228

                                  p0.ov), d_oo)

229

        ).antisymmetrise(0, 1)

230

232

        # N^4: O^2V^2 / N^4: O^2V^2

233

        + 4.0 * einsum("ib,ja->ijab",

234

                       einsum("kc,ikbc->ib", u1, td2),

235

                       u1).antisymmetrise(0, 1).antisymmetrise(2, 3)

236

        # N^5: O^3V^2 / N^4: O^2V^2

237

        + 2.0 * einsum("ijka,kb->ijab",

238

                       einsum("kc,ijac->ijka", u1, td2), u1).antisymmetrise(2, 3)

239

        # N^5: O^3V^2 / N^4: O^2V^2

240

        + 2.0 * einsum("jk,ikab->ijab", einsum("jc,kc->jk", u1, u1), td2)

241

243

        # N^6: O^3V^3 / N^4: O^2V^2

244

        - 4.0 * einsum("jkac,ikbc->iajb", u2, u2)

245

        # N^5: O^2V^3 / N^4: O^2V^2

246

        + 2.0 * einsum("ab,ij->iajb", einsum("klac,klbc->ab", u2, u2), d_oo)

247

        # N^6: O^4V^2 / N^4: O^2V^2

248

        + 1.0 * einsum("iklb,jkla->iajb", einsum("ld,ikbd->iklb", u1, t2),

249

                       einsum("lc,jkac->jkla", u1, t2))

250

        # N^4: O^2V^2 / N^4: O^2V^2

251

        + 1.0 * einsum("ab,ij->iajb", einsum("ka,kb->ab", u1, u1), p0.oo)

252

        # N^4: O^2V^2 / N^4: O^2V^2

253

        - 1.0 * einsum("ij,ab->iajb", einsum("ic,jc->ij", u1, u1), p0.vv)

254

        # N^6: O^3V^3 / N^4: O^2V^2

255

        + 1.0 * einsum("jlac,ilbc->iajb",

256

                       einsum("kl,jkac->jlac", einsum("kd,ld->kl", u1, u1), t2), t2)

257

        # N^6: O^4V^2 / N^4: O^2V^2

258

        + 0.5 * einsum("ijla,lb->iajb",

259

                       einsum("ijkl,ka->ijla",

260

                              einsum("ikcd,jlcd->ijkl", t2, t2), u1), u1)

261

        # N^6: O^4V^2 / N^4: O^2V^2

262

        + 0.5 * einsum("jklb,ikla->iajb",

263

                       einsum("jd,klbd->jklb", u1, t2),

264

                       einsum("ic,klac->ikla", u1, t2))

265

        # N^4: O^2V^2 / N^4: O^2V^2

266

        - 1.0 * einsum("ib,ja->iajb",

267

                       einsum("kd,ikbd->ib", u1, t2), einsum("lc,jlac->ja", u1, t2))

268

        # N^5: O^2V^3 / N^4: O^2V^2

269

        + 1.0 * einsum("ab,ij->iajb",

270

                       einsum("lmad,lmbd->ab",

271

                              einsum("km,klad->lmad",

272

                                     einsum("kc,mc->km", u1, u1), t2), t2), d_oo)

273

        # N^4: O^2V^2 / N^4: O^2V^2

274

        - 1.0 * einsum("ab,ij->iajb",

275

                       einsum("lb,la->ab",

276

                              einsum("kd,klbd->lb", u1, t2),

277

                              einsum("mc,lmac->la", u1, t2)), d_oo)

278

        # N^5: O^3V^2 / N^4: O^2V^2

279

        - 0.5 * einsum("ab,ij->iajb",

280

                       einsum("klmb,klma->ab",

281

                              einsum("md,klbd->klmb", u1, t2),

282

                              einsum("mc,klac->klma", u1, t2)), d_oo)

283

+ (

284

            # N^4: O^2V^2 / N^4: O^2V^2

285

            + 1.0 * einsum("ib,ja->iajb", einsum("ic,bc->ib", u1, p0.vv), u1)

286

            # N^4: O^2V^2 / N^4: O^2V^2

287

            - 1.0 * einsum("ib,ja->iajb", einsum("kb,ik->ib", u1, p0.oo), u1)

288

            # N^5: O^3V^2 / N^4: O^2V^2

289

            + 2.0 * einsum("ijlb,la->iajb",

290

                           einsum("jc,ilbc->ijlb",

291

                                  einsum("kd,jkcd->jc", u1, t2), t2), u1)

292

            # N^5: O^3V^2 / N^4: O^2V^2

293

            + 2.0 * einsum("kb,ijka->iajb",

294

                           einsum("ld,klbd->kb", u1, t2),

295

                           einsum("ic,jkac->ijka", u1, t2))

296

            # N^6: O^4V^2 / N^4: O^2V^2

297

            - 2.0 * einsum("ijkb,ka->iajb",

298

                           einsum("jklc,ilbc->ijkb",

299

                                  einsum("kd,jlcd->jklc", u1, t2), t2), u1)

300

            # N^6: O^4V^2 / N^4: O^2V^2

301

            - 2.0 * einsum("ijlb,la->iajb",

302

                           einsum("ijkc,klbc->ijlb",

303

                                  einsum("id,jkcd->ijkc", u1, t2), t2), u1)

304

            # N^6: O^3V^3 / N^4: O^2V^2

305

            - 2.0 * einsum("ikbc,jkac->iajb",

306

                           einsum("il,klbc->ikbc",

307

                                  einsum("id,ld->il", u1, u1), t2), t2)

308

            # N^4: O^2V^2 / N^4: O^2V^2

309

            - 1.0 * einsum("ib,ja->iajb",

310

                           einsum("kc,ikbc->ib",

311

                                  einsum("ld,klcd->kc", u1, t2), t2), u1)

312

            # N^4: O^2V^2 / N^4: O^2V^2

313

            - 1.0 * einsum("ab,ij->iajb",

314

                           einsum("kb,ka->ab",

315

                                  einsum("kc,bc->kb", u1, p0.vv), u1), d_oo)

316

            # N^4: O^2V^2 / N^4: O^2V^2

317

            - 1.0 * einsum("ab,ij->iajb",

318

                           einsum("mb,ma->ab",

319

                                  einsum("ld,lmbd->mb",

320

                                         einsum("kc,klcd->ld", u1, t2),

321

                                         t2), u1), d_oo)

322

        ).symmetrise([(0, 2), (1, 3)])

323

325

        # N^6: O^3V^3 / N^4: O^1V^3

326

        + 1.0 * einsum("ijka,jkbc->iabc", einsum("id,jkad->ijka", u1, u2), t2)

327

        # N^5: O^2V^3 / N^4: O^1V^3

328

        + 2.0 * einsum("ja,ijbc->iabc", einsum("kd,jkad->ja", u1, u2), t2)

329

+ (

330

            # N^5: O^2V^3 / N^4: O^1V^3

331

            + 2.0 * einsum("ac,ib->iabc", einsum("jkad,jkcd->ac", u2, t2), u1)

332

            # N^6: O^3V^3 / N^4: O^1V^3

333

            - 4.0 * einsum("ikab,kc->iabc", einsum("jkad,ijbd->ikab", u2, t2), u1)

334

            # N^4: O^1V^3 / N^4: O^1V^3

335

            + 2.0 * einsum("ac,ib->iabc", einsum("ja,jc->ac", u1, p0.ov), u1)

336

            # N^4: O^1V^3 / N^4: O^1V^3

337

            + 2.0 * einsum("ac,ib->iabc", einsum("ja,jc->ac", u1, u1), p0.ov)

338

        ).antisymmetrise(2, 3)

339

341

        # N^6: O^2V^4 / N^4: V^4

342

        + 2 * einsum("ijab,ijcd->abcd", u2, u2)

343

        # N^6: O^2V^4 / N^4: V^4

344

        - 1 * einsum("jkab,jkcd->abcd",

345

                     einsum("ij,ikab->jkab", einsum("ie,je->ij", u1, u1), t2), t2)

346

+ (

347

            # N^4: V^4 / N^4: V^4

348

            + 4.0 * einsum("ac,bd->abcd", einsum("ia,ic->ac", u1, u1), p0.vv)

349

            # N^6: O^3V^3 / N^4: V^4

350

            + 4.0 * einsum("jabc,jd->abcd",

351

                           einsum("ijbc,ia->jabc",

352

                                  einsum("jkbe,ikce->ijbc", t2, t2), u1), u1)

353

        ).antisymmetrise(0, 1).antisymmetrise(2, 3)

354

+ (

355

            # N^5: O^1V^4 / N^4: V^4

356

            + 4.0 * einsum("ka,kbcd->abcd",

357

                           einsum("ie,ikae->ka", u1, t2),

358

                           einsum("jb,jkcd->kbcd", u1, t2))

359

            # N^6: O^3V^3 / N^4: V^4

360

            + 2.0 * einsum("ibcd,ia->abcd",

361

                           einsum("ijkb,jkcd->ibcd",

362

                                  einsum("ie,jkbe->ijkb", u1, t2), t2), u1)

363

        ).antisymmetrise(0, 1).symmetrise([(0, 2), (1, 3)])

364

366

367

369

                   intermediates: Intermediates) -> TwoParticleDensity:

371

    # evaluate additional contributions from the S-T block

372

    # if the vector has a triples component

378

    # N^6: O^3V^3 / N^6: O^3V^3

381

382

383

# dict controlling the dispatch of the state_diffdm function

385

    "isr0": diffdm_isr0_2p,

386

    "isr1s": diffdm_isr1s_2p,

387

    "isr1": diffdm_isr1_2p,

388

    "isr2d": diffdm_isr2d_2p,

389

    "isr2": diffdm_isr2_2p,

390

391

392

394

"""

395

    Compute the two-particle difference density matrix of an excited state

396

    in the MO basis.

397

398

    Parameters

399

    ----------

400

    method : str, IsrMethod

401

        The method to use for the computation (e.g. "isr2")

402

    ground_state : LazyMp

403

        The ground state upon which the excitation was based

404

    amplitude : AmplitudeVector

405

        The amplitude vector

406

    intermediates : adcc.Intermediates

407

        Intermediates from the ADC calculation to reuse

408

"""

417

420

                                  f"for {method.name}.")

421

    else:

adc-connect / adcc / 25782308614

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

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 .. import block as b	2✔
24	from ..LazyMp import LazyMp	2✔
25	from ..AdcMethod import IsrMethod	2✔
26	from ..functions import einsum, zeros_like	2✔
27	from ..Intermediates import Intermediates	2✔
28	from ..AmplitudeVector import AmplitudeVector	2✔
29	from ..TwoParticleDensity import TwoParticleDensity	2✔
30	from ..NParticleOperator import OperatorSymmetry	2✔
31	from .util import (	2✔
32	check_doubles_amplitudes, check_singles_amplitudes, check_triples_amplitudes
33	)
34
35
36	def diffdm_isr0_2p(mp, amplitude, intermediates):	2✔
37	check_singles_amplitudes([b.o, b.v], amplitude)	2✔
38	u1 = amplitude.ph	2✔
39
40	hf = mp.reference_state	2✔
41	d_oo = zeros_like(hf.foo)	2✔
42	d_oo.set_mask("ii", 1)	2✔
43
44	dm = TwoParticleDensity(mp, symmetry=OperatorSymmetry.HERMITIAN)	2✔
45
46	# one-particle ISR(0) diffdm
47	p1_oo = -einsum("ia,la->il", u1, u1).evaluate()	2✔
48	p1_vv = einsum("ka,kb->ab", u1, u1).evaluate()	2✔
49
50	dm.oooo = (	2✔
51	# N^4: O^4 / N^4: O^4
52	- 4.0 * einsum("il,jk->ijkl", p1_oo, d_oo)
53	).antisymmetrise(0, 1).antisymmetrise(2, 3)
54	dm.ovov = (	2✔
55	# N^4: O^2V^2 / N^4: O^2V^2
56	- 1.0 * einsum("ja,ib->iajb", u1, u1)
57	# N^4: O^2V^2 / N^4: O^2V^2
58	+ 1.0 * einsum("ab,ij->iajb", p1_vv, d_oo)
59	)
60	return dm	2✔
61
62
63	def diffdm_isr1s_2p(mp, amplitude, intermediates):	2✔
64	dm = diffdm_isr0_2p(mp, amplitude, intermediates) # Get ISR(0) result	2✔
65	u1 = amplitude.ph	2✔
66
67	hf = mp.reference_state	2✔
68	d_oo = zeros_like(hf.foo)	2✔
69	d_oo.set_mask("ii", 1)	2✔
70
71	t2 = mp.t2(b.oovv)	2✔
72	# TODO move to intermediates!
73	# one-particle ISR(0) diffdm
74	p1_oo = -einsum("ia,la->il", u1, u1).evaluate()	2✔
75
76	# ISR(2) ISR intermediate (TODO Move to intermediates)
77	ru1 = einsum("ijab,jb->ia", t2, u1).evaluate()	2✔
78	# new ones
79	ru1_ooov = einsum("kc,ijac->ijka", u1, t2).evaluate()	2✔
80
81	dm.oovv += (	2✔
82	# N^4: O^2V^2 / N^4: O^2V^2
83	+ 4.0 * (
84	einsum("ib,ja->ijab", ru1, u1)
85	).antisymmetrise(0, 1).antisymmetrise(2, 3)
86	# N^5: O^3V^2 / N^4: O^2V^2
87	+ 2.0 * (
88	einsum("ijka,kb->ijab", ru1_ooov, u1)
89	).antisymmetrise(2, 3)
90	# N^5: O^3V^2 / N^4: O^2V^2
91	- 2.0 * (
92	einsum("jk,ikab->ijab", p1_oo, t2)
93	).antisymmetrise(0, 1)
94	)
95
96	return dm	2✔
97
98
99	def diffdm_isr1_2p(mp, amplitude, intermediates):	2✔
100	dm = diffdm_isr1s_2p(mp, amplitude, intermediates) # Get ISR(1)-s result	2✔
101
102	try:	2✔
103	# ISR(1)-d
104	check_doubles_amplitudes([b.o, b.o, b.v, b.v], amplitude)	2✔
105	u1, u2 = amplitude.ph, amplitude.pphh	2✔
106
107	hf = mp.reference_state	2✔
108	d_oo = zeros_like(hf.foo)	2✔
109	d_oo.set_mask("ii", 1)	2✔
110
111	dm.ooov += (	2✔
112	# N^5: O^3V^2 / N^4: O^2V^2
113	- 2.0 * einsum("kb,ijab->ijka", u1, u2)
114	# N^4: O^2V^2 / N^4: O^2V^2
115	- 4.0 * einsum(
116	"ja,ik->ijka", einsum("lb,jlab->ja", u1, u2), d_oo
117	).antisymmetrise(0, 1)
118	)
119
120	dm.ovvv += (	2✔
121	# N^5: O^2V^3 / N^4: O^1V^3
122	- 2.0 * einsum("ja,ijbc->iabc", u1, u2)
123	)
124	except ValueError:	×
125	# no doubles contribution
126	pass	×
127	return dm	2✔
128
129
130	def diffdm_isr2d_2p(mp, amplitude, intermediates):	2✔
131	dm = diffdm_isr1_2p(mp, amplitude, intermediates) # Get ISR(1) result	2✔
132	check_doubles_amplitudes([b.o, b.o, b.v, b.v], amplitude)	2✔
133	u1, u2 = amplitude.ph, amplitude.pphh	2✔
134	hf = mp.reference_state	2✔
135	d_oo = zeros_like(hf.foo)	2✔
136	d_oo.set_mask("ii", 1)	2✔
137
138	t2 = mp.t2(b.oovv)	2✔
139	td2 = mp.td2(b.oovv)	2✔
140	p0 = mp.mp2_diffdm	2✔
141
142	dm.oooo += (	2✔
143	# N^6: O^4V^2 / N^4: O^2V^2
144	+ 2.0 * einsum("ijab,klab->ijkl", u2, u2)
145	# N^6: O^5V^1 / N^4: O^2V^2
146	- 1.0 * einsum("klmc,ijmc->ijkl", einsum("mb,klbc->klmc", u1, t2),
147	einsum("ma,ijac->ijmc", u1, t2))
148	+ (
149	# N^5: O^3V^2 / N^4: O^2V^2
150	- 8.0 * einsum("ik,jl->ijkl", einsum("imab,kmab->ik", u2, u2), d_oo)
151	# N^4: O^4 / N^4: O^4
152	+ 4.0 * einsum("il,jk->ijkl", einsum("ia,la->il", u1, u1), p0.oo)
153	# N^6: O^5V^1 / N^4: O^2V^2
154	+ 4.0 * einsum("jlmc,ikmc->ijkl", einsum("lb,jmbc->jlmc", u1, t2),
155	einsum("ia,kmac->ikmc", u1, t2))
156	# N^4: O^2V^2 / N^4: O^2V^2
157	+ 4.0 * einsum("il,jk->ijkl",
158	einsum("lc,ic->il", einsum("nb,lnbc->lc", u1, t2),
159	einsum("ma,imac->ic", u1, t2)), d_oo)
160	# N^5: O^3V^2 / N^4: O^2V^2
161	+ 4.0 * einsum("ik,jl->ijkl", einsum("kmnc,imnc->ik",
162	einsum("nb,kmbc->kmnc", u1, t2),
163	einsum("na,imac->imnc", u1, t2)),
164	d_oo)
165	# N^5: O^3V^2 / N^4: O^2V^2
166	+ 2.0 * einsum("ik,jl->ijkl",
167	einsum("inbc,knbc->ik",
168	einsum("mn,imbc->inbc",
169	einsum("ia,la->il", u1, u1), t2), t2),
170	d_oo)
171	).antisymmetrise(0, 1).antisymmetrise(2, 3)
172	+ (
173	# N^5: O^3V^2 / N^4: O^2V^2
174	+ 4.0 * einsum("iklb,jb->ijkl",
175	einsum("ic,klbc->iklb", einsum("ma,imac->ic", u1, t2),
176	t2), u1)
177	# N^6: O^4V^2 / N^4: O^2V^2
178	+ 2.0 * einsum("jklm,im->ijkl",
179	einsum("jmbc,klbc->jklm", t2, t2),
180	einsum("ia,la->il", u1, u1))
181	).antisymmetrise(0, 1).symmetrise([(0, 2), (1, 3)])
182	+ (
183	# N^4: O^4 / N^4: O^4
184	+ 4.0 * einsum("il,jk->ijkl", einsum("im,lm->il",
185	einsum("ia,la->il", u1, u1),
186	p0.oo), d_oo)
187	# N^4: O^2V^2 / N^4: O^2V^2
188	+ 4.0 * einsum("ik,jl->ijkl",
189	einsum("kb,ib->ik",
190	einsum("mc,kmbc->kb",
191	einsum("na,mnac->mc", u1, t2), t2), u1),
192	d_oo)
193	).antisymmetrise(0, 1).antisymmetrise(2, 3).symmetrise([(0, 2), (1, 3)])
194	)
195	dm.ooov += (	2✔
196	# N^6: O^4V^2 / N^4: O^2V^2
197	+ 1.0 * einsum("ijkl,la->ijka", einsum("klbc,ijbc->ijkl", u2, t2), u1)
198	# N^5: O^3V^2 / N^4: O^2V^2
199	+ 2.0 * einsum("kb,ijab->ijka", einsum("lc,klbc->kb", u1, u2), t2)
200	+ (
201	# N^5: O^3V^2 / N^4: O^2V^2
202	+ 2.0 * einsum("jk,ia->ijka", einsum("klbc,jlbc->jk", u2, t2), u1)
203	# N^6: O^4V^2 / N^4: O^2V^2
204	- 4.0 * einsum("jklb,ilab->ijka", einsum("jc,klbc->jklb", u1, u2), t2)
205	# N^5: O^3V^2 / N^4: O^2V^2
206	+ 2.0 * einsum("ja,ik->ijka",
207	einsum("jlmb,lmab->ja",
208	einsum("jc,lmbc->jlmb", u1, u2), t2), d_oo)
209	# N^4: O^2V^2 / N^4: O^2V^2
210	- 4.0 * einsum("ia,jk->ijka",
211	einsum("lb,ilab->ia",
212	einsum("mc,lmbc->lb", u1, u2), t2), d_oo)
213	# N^5: O^3V^2 / N^4: O^2V^2
214	+ 2.0 * einsum("ia,jk->ijka",
215	einsum("il,la->ia",
216	einsum("lmbc,imbc->il", u2, t2), u1), d_oo)
217	# N^4: O^3V^1 / N^4: O^3V^1
218	+ 2.0 * einsum("jk,ia->ijka", einsum("kb,jb->jk", u1, p0.ov), u1)
219	# N^4: O^3V^1 / N^4: O^3V^1
220	+ 2.0 * einsum("jk,ia->ijka", einsum("jb,kb->jk", u1, u1), p0.ov)
221	# N^4: O^3V^1 / N^4: O^3V^1
222	+ 2.0 * einsum("ia,jk->ijka",
223	einsum("il,la->ia",
224	einsum("lb,ib->il", u1, p0.ov), u1), d_oo)
225	# N^4: O^3V^1 / N^4: O^3V^1
226	+ 2.0 * einsum("ia,jk->ijka",
227	einsum("il,la->ia", einsum("ia,la->il", u1, u1),
228	p0.ov), d_oo)
229	).antisymmetrise(0, 1)
230	)
231	dm.oovv += (	2✔
232	# N^4: O^2V^2 / N^4: O^2V^2
233	+ 4.0 * einsum("ib,ja->ijab",
234	einsum("kc,ikbc->ib", u1, td2),
235	u1).antisymmetrise(0, 1).antisymmetrise(2, 3)
236	# N^5: O^3V^2 / N^4: O^2V^2
237	+ 2.0 * einsum("ijka,kb->ijab",
238	einsum("kc,ijac->ijka", u1, td2), u1).antisymmetrise(2, 3)
239	# N^5: O^3V^2 / N^4: O^2V^2
240	+ 2.0 * einsum("jk,ikab->ijab", einsum("jc,kc->jk", u1, u1), td2)
241	)
242	dm.ovov += (	2✔
243	# N^6: O^3V^3 / N^4: O^2V^2
244	- 4.0 * einsum("jkac,ikbc->iajb", u2, u2)
245	# N^5: O^2V^3 / N^4: O^2V^2
246	+ 2.0 * einsum("ab,ij->iajb", einsum("klac,klbc->ab", u2, u2), d_oo)
247	# N^6: O^4V^2 / N^4: O^2V^2
248	+ 1.0 * einsum("iklb,jkla->iajb", einsum("ld,ikbd->iklb", u1, t2),
249	einsum("lc,jkac->jkla", u1, t2))
250	# N^4: O^2V^2 / N^4: O^2V^2
251	+ 1.0 * einsum("ab,ij->iajb", einsum("ka,kb->ab", u1, u1), p0.oo)
252	# N^4: O^2V^2 / N^4: O^2V^2
253	- 1.0 * einsum("ij,ab->iajb", einsum("ic,jc->ij", u1, u1), p0.vv)
254	# N^6: O^3V^3 / N^4: O^2V^2
255	+ 1.0 * einsum("jlac,ilbc->iajb",
256	einsum("kl,jkac->jlac", einsum("kd,ld->kl", u1, u1), t2), t2)
257	# N^6: O^4V^2 / N^4: O^2V^2
258	+ 0.5 * einsum("ijla,lb->iajb",
259	einsum("ijkl,ka->ijla",
260	einsum("ikcd,jlcd->ijkl", t2, t2), u1), u1)
261	# N^6: O^4V^2 / N^4: O^2V^2
262	+ 0.5 * einsum("jklb,ikla->iajb",
263	einsum("jd,klbd->jklb", u1, t2),
264	einsum("ic,klac->ikla", u1, t2))
265	# N^4: O^2V^2 / N^4: O^2V^2
266	- 1.0 * einsum("ib,ja->iajb",
267	einsum("kd,ikbd->ib", u1, t2), einsum("lc,jlac->ja", u1, t2))
268	# N^5: O^2V^3 / N^4: O^2V^2
269	+ 1.0 * einsum("ab,ij->iajb",
270	einsum("lmad,lmbd->ab",
271	einsum("km,klad->lmad",
272	einsum("kc,mc->km", u1, u1), t2), t2), d_oo)
273	# N^4: O^2V^2 / N^4: O^2V^2
274	- 1.0 * einsum("ab,ij->iajb",
275	einsum("lb,la->ab",
276	einsum("kd,klbd->lb", u1, t2),
277	einsum("mc,lmac->la", u1, t2)), d_oo)
278	# N^5: O^3V^2 / N^4: O^2V^2
279	- 0.5 * einsum("ab,ij->iajb",
280	einsum("klmb,klma->ab",
281	einsum("md,klbd->klmb", u1, t2),
282	einsum("mc,klac->klma", u1, t2)), d_oo)
283	+ (
284	# N^4: O^2V^2 / N^4: O^2V^2
285	+ 1.0 * einsum("ib,ja->iajb", einsum("ic,bc->ib", u1, p0.vv), u1)
286	# N^4: O^2V^2 / N^4: O^2V^2
287	- 1.0 * einsum("ib,ja->iajb", einsum("kb,ik->ib", u1, p0.oo), u1)
288	# N^5: O^3V^2 / N^4: O^2V^2
289	+ 2.0 * einsum("ijlb,la->iajb",
290	einsum("jc,ilbc->ijlb",
291	einsum("kd,jkcd->jc", u1, t2), t2), u1)
292	# N^5: O^3V^2 / N^4: O^2V^2
293	+ 2.0 * einsum("kb,ijka->iajb",
294	einsum("ld,klbd->kb", u1, t2),
295	einsum("ic,jkac->ijka", u1, t2))
296	# N^6: O^4V^2 / N^4: O^2V^2
297	- 2.0 * einsum("ijkb,ka->iajb",
298	einsum("jklc,ilbc->ijkb",
299	einsum("kd,jlcd->jklc", u1, t2), t2), u1)
300	# N^6: O^4V^2 / N^4: O^2V^2
301	- 2.0 * einsum("ijlb,la->iajb",
302	einsum("ijkc,klbc->ijlb",
303	einsum("id,jkcd->ijkc", u1, t2), t2), u1)
304	# N^6: O^3V^3 / N^4: O^2V^2
305	- 2.0 * einsum("ikbc,jkac->iajb",
306	einsum("il,klbc->ikbc",
307	einsum("id,ld->il", u1, u1), t2), t2)
308	# N^4: O^2V^2 / N^4: O^2V^2
309	- 1.0 * einsum("ib,ja->iajb",
310	einsum("kc,ikbc->ib",
311	einsum("ld,klcd->kc", u1, t2), t2), u1)
312	# N^4: O^2V^2 / N^4: O^2V^2
313	- 1.0 * einsum("ab,ij->iajb",
314	einsum("kb,ka->ab",
315	einsum("kc,bc->kb", u1, p0.vv), u1), d_oo)
316	# N^4: O^2V^2 / N^4: O^2V^2
317	- 1.0 * einsum("ab,ij->iajb",
318	einsum("mb,ma->ab",
319	einsum("ld,lmbd->mb",
320	einsum("kc,klcd->ld", u1, t2),
321	t2), u1), d_oo)
322	).symmetrise([(0, 2), (1, 3)])
323	)
324	dm.ovvv += (	2✔
325	# N^6: O^3V^3 / N^4: O^1V^3
326	+ 1.0 * einsum("ijka,jkbc->iabc", einsum("id,jkad->ijka", u1, u2), t2)
327	# N^5: O^2V^3 / N^4: O^1V^3
328	+ 2.0 * einsum("ja,ijbc->iabc", einsum("kd,jkad->ja", u1, u2), t2)
329	+ (
330	# N^5: O^2V^3 / N^4: O^1V^3
331	+ 2.0 * einsum("ac,ib->iabc", einsum("jkad,jkcd->ac", u2, t2), u1)
332	# N^6: O^3V^3 / N^4: O^1V^3
333	- 4.0 * einsum("ikab,kc->iabc", einsum("jkad,ijbd->ikab", u2, t2), u1)
334	# N^4: O^1V^3 / N^4: O^1V^3
335	+ 2.0 * einsum("ac,ib->iabc", einsum("ja,jc->ac", u1, p0.ov), u1)
336	# N^4: O^1V^3 / N^4: O^1V^3
337	+ 2.0 * einsum("ac,ib->iabc", einsum("ja,jc->ac", u1, u1), p0.ov)
338	).antisymmetrise(2, 3)
339	)
340	dm.vvvv += (	2✔
341	# N^6: O^2V^4 / N^4: V^4
342	+ 2 * einsum("ijab,ijcd->abcd", u2, u2)
343	# N^6: O^2V^4 / N^4: V^4
344	- 1 * einsum("jkab,jkcd->abcd",
345	einsum("ij,ikab->jkab", einsum("ie,je->ij", u1, u1), t2), t2)
346	+ (
347	# N^4: V^4 / N^4: V^4
348	+ 4.0 * einsum("ac,bd->abcd", einsum("ia,ic->ac", u1, u1), p0.vv)
349	# N^6: O^3V^3 / N^4: V^4
350	+ 4.0 * einsum("jabc,jd->abcd",
351	einsum("ijbc,ia->jabc",
352	einsum("jkbe,ikce->ijbc", t2, t2), u1), u1)
353	).antisymmetrise(0, 1).antisymmetrise(2, 3)
354	+ (
355	# N^5: O^1V^4 / N^4: V^4
356	+ 4.0 * einsum("ka,kbcd->abcd",
357	einsum("ie,ikae->ka", u1, t2),
358	einsum("jb,jkcd->kbcd", u1, t2))
359	# N^6: O^3V^3 / N^4: V^4
360	+ 2.0 * einsum("ibcd,ia->abcd",
361	einsum("ijkb,jkcd->ibcd",
362	einsum("ie,jkbe->ijkb", u1, t2), t2), u1)
363	).antisymmetrise(0, 1).symmetrise([(0, 2), (1, 3)])
364	)
365	return dm	2✔
366
367
368	def diffdm_isr2_2p(mp: LazyMp, amplitude: AmplitudeVector,	2✔
369	intermediates: Intermediates) -> TwoParticleDensity:
370	dm = diffdm_isr2d_2p(mp, amplitude, intermediates) # Get ISR(1) result	2✔
371	# evaluate additional contributions from the S-T block
372	# if the vector has a triples component
373	try:	2✔
374	check_triples_amplitudes([b.o, b.o, b.o, b.v, b.v, b.v], amplitude)	2✔
375	except ValueError:	2✔
376	return dm	2✔
NEW 377	u1, u3 = amplitude.ph, amplitude.ppphhh	×
378	# N^6: O^3V^3 / N^6: O^3V^3
NEW 379	dm.oovv += - 6 * einsum("kc,ijkabc->ijab", u1, u3)	×
NEW 380	return dm	×
381
382
383	# dict controlling the dispatch of the state_diffdm function
384	DISPATCH = {	2✔
385	"isr0": diffdm_isr0_2p,
386	"isr1s": diffdm_isr1s_2p,
387	"isr1": diffdm_isr1_2p,
388	"isr2d": diffdm_isr2d_2p,
389	"isr2": diffdm_isr2_2p,
390	}
391
392
393	def state_diffdm_2p(method, ground_state, amplitude, intermediates=None):	2✔
394	"""
395	Compute the two-particle difference density matrix of an excited state
396	in the MO basis.
397
398	Parameters
399	----------
400	method : str, IsrMethod
401	The method to use for the computation (e.g. "isr2")
402	ground_state : LazyMp
403	The ground state upon which the excitation was based
404	amplitude : AmplitudeVector
405	The amplitude vector
406	intermediates : adcc.Intermediates
407	Intermediates from the ADC calculation to reuse
408	"""
409	if not isinstance(method, IsrMethod):	2!
410	method = IsrMethod(method)	×
411	if not isinstance(ground_state, LazyMp):	2!
412	raise TypeError("ground_state should be a LazyMp object.")	×
413	if not isinstance(amplitude, AmplitudeVector):	2!
414	raise TypeError("amplitude should be an AmplitudeVector object.")	×
415	if intermediates is None:	2✔
416	intermediates = Intermediates(ground_state)	2✔
417
418	if method.name not in DISPATCH:	2!
419	raise NotImplementedError("state_diffdm_2p is not implemented "	×
420	f"for {method.name}.")
421	else:
422	ret = DISPATCH[method.name](ground_state, amplitude, intermediates)	2✔
423	return ret.evaluate()	2✔

adc-connect / adcc / 25782308614

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

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