adc-connect/adcc | Build 24572346840 | 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

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

31

33

34

38

42

44

45

    # one-particle ISR(0) diffdm

48

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)

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

60

61

65

69

71

    # TODO move to intermediates!

72

    # one-particle ISR(0) diffdm

74

75

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

77

    # new ones

79

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

96

97

101

105

107

    # TODO move to intermediates!

108

    # one-particle ISR(0) diffdm

110

111

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

113

    # new ones

115

117

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

118

        + 4.0 * (

119

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

120

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

121

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

122

        + 2.0 * (

123

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

124

        ).antisymmetrise(2, 3)

125

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

126

        - 2.0 * (

127

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

128

        ).antisymmetrise(0, 1)

129

130

132

        # ISR(1)-d

135

137

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

138

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

139

140

142

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

143

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

144

146

        # no doubles contribution

149

150

158

162

164

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

165

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

166

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

167

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

168

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

169

+ (

170

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

171

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

172

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

173

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

174

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

175

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

176

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

177

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

178

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

179

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

180

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

181

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

182

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

183

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

184

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

185

                           d_oo)

186

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

187

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

188

                           einsum("inbc,knbc->ik",

189

                                  einsum("mn,imbc->inbc",

190

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

191

                           d_oo)

192

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

193

+ (

194

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

195

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

196

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

197

                                  t2), u1)

198

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

199

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

200

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

201

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

202

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

203

+ (

204

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

205

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

206

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

207

                                                 p0.oo), d_oo)

208

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

209

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

210

                           einsum("kb,ib->ik",

211

                                  einsum("mc,kmbc->kb",

212

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

213

                           d_oo)

214

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

215

217

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

218

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

219

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

220

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

221

+ (

222

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

223

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

224

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

225

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

226

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

227

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

228

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

229

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

230

                           einsum("jlmb,lmab->ja",

231

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

232

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

233

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

234

                           einsum("lb,ilab->ia",

235

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

236

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

237

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

238

                           einsum("il,la->ia",

239

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

240

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

241

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

242

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

243

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

244

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

245

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

246

                           einsum("il,la->ia",

247

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

248

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

249

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

250

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

251

                                  p0.ov), d_oo)

252

        ).antisymmetrise(0, 1)

253

255

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

256

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

257

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

258

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

259

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

260

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

261

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

262

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

263

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

264

266

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

267

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

268

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

269

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

270

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

271

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

272

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

273

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

274

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

275

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

276

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

277

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

278

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

279

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

280

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

281

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

282

                       einsum("ijkl,ka->ijla",

283

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

284

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

285

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

286

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

287

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

288

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

289

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

290

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

291

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

292

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

293

                       einsum("lmad,lmbd->ab",

294

                              einsum("km,klad->lmad",

295

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

296

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

297

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

298

                       einsum("lb,la->ab",

299

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

300

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

301

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

302

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

303

                       einsum("klmb,klma->ab",

304

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

305

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

306

+ (

307

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

308

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

309

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

310

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

311

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

312

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

313

                           einsum("jc,ilbc->ijlb",

314

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

315

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

316

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

317

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

318

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

319

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

320

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

321

                           einsum("jklc,ilbc->ijkb",

322

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

323

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

324

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

325

                           einsum("ijkc,klbc->ijlb",

326

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

327

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

328

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

329

                           einsum("il,klbc->ikbc",

330

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

331

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

332

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

333

                           einsum("kc,ikbc->ib",

334

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

335

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

336

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

337

                           einsum("kb,ka->ab",

338

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

339

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

340

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

341

                           einsum("mb,ma->ab",

342

                                  einsum("ld,lmbd->mb",

343

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

344

                                         t2), u1), d_oo)

345

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

346

348

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

349

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

350

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

351

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

352

+ (

353

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

354

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

355

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

356

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

357

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

358

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

359

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

360

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

361

        ).antisymmetrise(2, 3)

362

364

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

365

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

366

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

367

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

368

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

369

+ (

370

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

371

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

372

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

373

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

374

                           einsum("ijbc,ia->jabc",

375

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

376

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

377

+ (

378

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

379

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

380

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

381

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

382

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

383

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

384

                           einsum("ijkb,jkcd->ibcd",

385

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

386

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

387

389

390

391

# dict controlling the dispatch of the state_diffdm function

393

    "isr0": diffdm_isr0_2p,

394

    "isr1s": diffdm_isr1s_2p,

395

    "isr1": diffdm_isr1_2p,

396

    "isr2": diffdm_isr2_2p,

397

398

399

401

"""

402

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

403

    in the MO basis.

404

405

    Parameters

406

    ----------

407

    method : str, IsrMethod

408

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

409

    ground_state : LazyMp

410

        The ground state upon which the excitation was based

411

    amplitude : AmplitudeVector

412

        The amplitude vector

413

    intermediates : adcc.Intermediates

414

        Intermediates from the ADC calculation to reuse

415

"""

424

427

                                  f"for {method.name}.")

428

    else:

adc-connect / adcc / 24572346840

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

adc-connect / adcc / 24572346840

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