add GHF-based pprpa and gradient for rpprpa

examples/grad/05-rpprpa_grad.py

@@ -0,0 +1,95 @@
+r"""
+    This is an example to calculate the molecular geometrical gradient of the ppRPA energy.
+"""
+from pyscf import gto, scf, dft
+from pyscf.pprpa import rpprpa_davidson, rpprpa_direct
+from pyscf.grad import rpprpa
+import numpy as np
+
+mult = "s"
+istate = 0
+nocc = 5
+nvir = 10
+
+def mfobj(dx):
+    mol = gto.Mole()
+    mol.verbose = 0
+    mol.atom = [["O", (0.0, 0.0, 0.0)],
+                ["H", (0.0, 1.0, 1.0+dx)],
+                ["H", (0.0, -1.0, 1.0)]]
+    mol.basis = "cc-pvdz"
+    mol.charge = 0
+    mol.unit = "Bohr"
+    mol.build()
+
+    # density fitting is required to ensure the numerical and analytical gradients are the same
+    # mf = scf.RHF(mol).density_fit()
+
+    # DFT grid response not implemented in pprpa
+    # Higher grid level for better accuracy
+    mf = dft.RKS(mol, xc="b3lyp").density_fit()
+    mf.grids.level = 7
+    mf.conv_tol = 1e-12
+    return mf
+
+def pprpaobj(mf, nocc_act, nvir_act):
+    mp = rpprpa_davidson.RppRPADavidson(mf, nocc_act, nvir_act, channel="pp", nroot=3)
+    mp.mu = 0.0
+    return mp
+
+
+
+
+dx = 0.001
+mf1 = mfobj(dx)
+mf2 = mfobj(-dx)
+
+e1 = mf1.kernel()
+# f1 = mf1.get_fock(dm=dm0_hf)
+h1 = mf1.get_hcore()
+
+e2 = mf2.kernel()
+# f2 = mf2.get_fock(dm=dm0_hf)
+h2 = mf2.get_hcore()
+
+# dfock = (f1 - f2)/(2.0*dx)
+dh = (h1 - h2)/(2.0*dx)
+e_hf = (e1 - e2)/dx/2.0
+
+mp1 = pprpaobj(mf1, nocc, nvir)
+mp1.kernel(mult)
+mp1.analyze()
+if mult == "s":
+    e1 += mp1.exci_s[istate] if mp1.channel == "pp" else -mp1.exci_s[istate]
+else:
+    e1 += mp1.exci_t[istate] if mp1.channel == "pp" else -mp1.exci_t[istate]
+mp2 = pprpaobj(mf2, nocc, nvir)
+mp2.kernel(mult)
+mp2.analyze()
+if mult == "s":
+    e2 += mp2.exci_s[istate] if mp2.channel == "pp" else -mp2.exci_s[istate]
+else:
+    e2 += mp2.exci_t[istate] if mp2.channel == "pp" else -mp2.exci_t[istate]
+
+e_rpa = (e1 - e2)/dx/2.0
+
+mf = mfobj(0.0)
+mf.kernel()
+if mult == "s":
+    oo_dim = nocc * (nocc + 1) // 2
+else:
+    oo_dim = nocc * (nocc - 1) // 2
+mp = pprpaobj(mf, nocc, nvir)
+mp.kernel(mult)
+mp.analyze()
+
+xy = mp.xy_s[istate] if mult == "s" else mp.xy_t[istate]
+from lib_pprpa.grad import pprpa
+mpg = mp.Gradients(mult,istate)
+mpg.kernel()
+
+print("analytical (Total):  ", mpg.de)
+print("numerical (Total):  ", e_rpa)
+print("diff: ", e_rpa - mpg.de[1,2])
+# print(e1,e2)
+print(mpg.de[0]+mpg.de[1]+mpg.de[2]) # Should be zero

examples/pprpa/06-gpprpa_excitation_energy.py

@@ -0,0 +1,48 @@
+from pyscf import gto, dft
+
+# For ppRPA excitation energy of N-electron system in particle-particle channel
+# mean field is (N-2)-electron
+
+mol = gto.Mole()
+mol.verbose = 4
+mol.atom = [
+    ["O", (0.0, 0.0, 0.0)],
+    ["H", (0.0, -0.7571, 0.5861)],
+    ["H", (0.0, 0.7571, 0.5861)]]
+mol.basis = 'def2-svp'
+# create a (N-2)-electron system for charged-neutral H2O
+mol.charge = 2
+mol.build()
+
+# =====> Part I. Restricted and its real-valued generalized ppRPA <=====
+# restricted KS-DFT calculation as starting point of RppRPA
+rmf = dft.RKS(mol)
+rmf.xc = "b3lyp"
+rmf.kernel()
+from pyscf.pprpa.rpprpa_davidson import RppRPADavidson
+pp = RppRPADavidson(rmf, nocc_act=3, nvir_act=None, nroot=1)
+# solve for singlet states
+pp.kernel("s")
+# solve for triplet states
+pp.kernel("t")
+pp.analyze()
+
+# Davidson algorithm for GKS, equivalent to the above RKS case
+from pyscf.pprpa.gpprpa_davidson import GppRPADavidson
+gmf = rmf.to_gks()
+pp = GppRPADavidson(gmf, nocc_act=6, nvir_act=None, nroot=4)
+# solve for singlet and triplet states
+pp.kernel()
+pp.analyze()
+
+
+
+# =====> Part II. Complex-valued generalized ppRPA <=====
+gmf = dft.GKS(mol).x2c1e()
+gmf.xc = "b3lyp"
+gmf.kernel()
+pp = GppRPADavidson(gmf, nocc_act=6, nvir_act=None, nroot=4)
+# solve for "singlet" and "triplet" states
+# triplets are no longer degenerate with SOC
+pp.kernel()
+pp.analyze()