Skip to content

Commit e2ce5ae

Browse files
tberkeltberkel
authored
Merge pull request #173 from lijiachen417/master
add documentation for particle-particle Random Phase Approximation
2 parents c5548e2 + 1fcff12 commit e2ce5ae

File tree

10 files changed

+628
-3
lines changed

10 files changed

+628
-3
lines changed

‎examples/pprpa/01-pprpa_total_energy.py‎

Lines changed: 42 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,42 @@
1+
from pyscf import gto, dft
2+
3+
# For ppRPA total energy of N-electron system
4+
# mean field is also N-electron
5+
6+
mol = gto.Mole()
7+
mol.verbose = 5
8+
mol.atom = [
9+
["O", (0.0, 0.0, 0.0)],
10+
["H", (0.0, -0.7571, 0.5861)],
11+
["H", (0.0, 0.7571, 0.5861)]]
12+
mol.basis = 'def2-svp'
13+
mol.build()
14+
15+
# =====> Part I. Restricted ppRPA <=====
16+
mf = dft.RKS(mol)
17+
mf.xc = "b3lyp"
18+
mf.kernel()
19+
20+
from pyscf.pprpa.rpprpa_direct import RppRPADirect
21+
pp = RppRPADirect(mf)
22+
ec = pp.get_correlation()
23+
etot, ehf, ec = pp.energy_tot()
24+
print("H2O Hartree-Fock energy = %.8f" % ehf)
25+
print("H2O ppRPA correlation energy = %.8f" % ec)
26+
print("H2O ppRPA total energy = %.8f" % etot)
27+
28+
# =====> Part II. Unrestricted ppRPA <=====
29+
# unrestricted KS-DFT calculation as starting point of UppRPA
30+
umf = dft.UKS(mol)
31+
umf.xc = "b3lyp"
32+
umf.kernel()
33+
34+
# direct diagonalization, N6 scaling
35+
from pyscf.pprpa.upprpa_direct import UppRPADirect
36+
pp = UppRPADirect(umf)
37+
ec = pp.get_correlation()
38+
etot, ehf, ec = pp.energy_tot()
39+
print("H2O Hartree-Fock energy = %.8f" % ehf)
40+
print("H2O ppRPA correlation energy = %.8f" % ec)
41+
print("H2O ppRPA total energy = %.8f" % etot)
42+

‎examples/pprpa/02-pprpa_excitation_energy.py‎

Lines changed: 64 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,64 @@
1+
from pyscf import gto, dft
2+
3+
# For ppRPA excitation energy of N-electron system in particle-particle channel
4+
# mean field is (N-2)-electron
5+
6+
mol = gto.Mole()
7+
mol.verbose = 5
8+
mol.atom = [
9+
["O", (0.0, 0.0, 0.0)],
10+
["H", (0.0, -0.7571, 0.5861)],
11+
["H", (0.0, 0.7571, 0.5861)]]
12+
mol.basis = 'def2-svp'
13+
# create a (N-2)-electron system for charged-neutral H2O
14+
mol.charge = 2
15+
mol.build()
16+
17+
# =====> Part I. Restricted ppRPA <=====
18+
# restricted KS-DFT calculation as starting point of RppRPA
19+
rmf = dft.RKS(mol)
20+
rmf.xc = "b3lyp"
21+
rmf.kernel()
22+
23+
# direct diagonalization, N6 scaling
24+
from pyscf.pprpa.rpprpa_direct import RppRPADirect
25+
# ppRPA can be solved in an active space
26+
pp = RppRPADirect(rmf, nocc_act=None, nvir_act=10)
27+
# number of two-electron addition states to print
28+
pp.pp_state = 10
29+
# solve for singlet states
30+
pp.kernel("s")
31+
# solve for triplet states
32+
pp.kernel("t")
33+
pp.analyze()
34+
35+
# Davidson algorithm, N4 scaling
36+
from pyscf.pprpa.rpprpa_davidson import RppRPADavidson
37+
# ppRPA can be solved in an active space
38+
pp = RppRPADavidson(rmf, nocc_act=3, nvir_act=None, nroot=10)
39+
# solve for singlet states
40+
pp.kernel("s")
41+
# solve for triplet states
42+
pp.kernel("t")
43+
pp.analyze()
44+
45+
# =====> Part II. Unrestricted ppRPA <=====
46+
# unrestricted KS-DFT calculation as starting point of UppRPA
47+
umf = dft.UKS(mol)
48+
umf.xc = "b3lyp"
49+
umf.kernel()
50+
51+
# direct diagonalization, N6 scaling
52+
from pyscf.pprpa.upprpa_direct import UppRPADirect
53+
# ppRPA can be solved in an active space
54+
pp = UppRPADirect(umf, nocc_act=None, nvir_act=10)
55+
# number of two-electron addition states to print
56+
pp.pp_state = 10
57+
# solve ppRPA in the (alpha alpha, alpha alpha) subspace
58+
pp.kernel(subspace=['aa'])
59+
# solve ppRPA in the (alpha beta, alpha beta) subspace
60+
pp.kernel(subspace=['ab'])
61+
# solve ppRPA in the (beta beta, beta beta) subspace
62+
pp.kernel(subspace=['bb'])
63+
pp.analyze()
64+

‎examples/pprpa/03-hhrpa_excitation_energy.py‎

Lines changed: 63 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,63 @@
1+
from pyscf import gto, dft
2+
3+
# For ppRPA excitation energy of N-electron system in hole-hole channel
4+
# mean field is (N+2)-electron
5+
6+
mol = gto.Mole()
7+
mol.verbose = 5
8+
mol.atom = [
9+
["O", (0.0, 0.0, 0.0)],
10+
["H", (0.0, -0.7571, 0.5861)],
11+
["H", (0.0, 0.7571, 0.5861)]]
12+
mol.basis = 'def2-svp'
13+
# create a (N+2)-electron system for charged-neutral H2O
14+
mol.charge = -2
15+
mol.build()
16+
17+
# =====> Part I. Restricted ppRPA <=====
18+
mf = dft.RKS(mol)
19+
mf.xc = "b3lyp"
20+
mf.kernel()
21+
22+
# direct diagonalization, N6 scaling
23+
# ppRPA can be solved in an active space
24+
from pyscf.pprpa.rpprpa_direct import RppRPADirect
25+
pp = RppRPADirect(mf, nvir_act=10, nelec="n+2")
26+
# number of two-electron removal states to print
27+
pp.hh_state = 10
28+
# solve for singlet states
29+
pp.kernel("s")
30+
# solve for triplet states
31+
pp.kernel("t")
32+
pp.analyze()
33+
34+
# Davidson algorithm, N4 scaling
35+
# ppRPA can be solved in an active space
36+
from pyscf.pprpa.rpprpa_davidson import RppRPADavidson
37+
pp = RppRPADavidson(mf, nvir_act=10, channel="hh")
38+
# solve for singlet states
39+
pp.kernel("s")
40+
# solve for triplet states
41+
pp.kernel("t")
42+
pp.analyze()
43+
44+
# =====> Part II. Unrestricted ppRPA <=====
45+
# unrestricted KS-DFT calculation as starting point of UppRPA
46+
umf = dft.UKS(mol)
47+
umf.xc = "b3lyp"
48+
umf.kernel()
49+
50+
# direct diagonalization, N6 scaling
51+
from pyscf.pprpa.upprpa_direct import UppRPADirect
52+
# ppRPA can be solved in an active space
53+
pp = UppRPADirect(umf, nvir_act=10, nelec='n+2')
54+
# number of two-electron addition states to print
55+
pp.pp_state = 10
56+
# solve ppRPA in the (alpha alpha, alpha alpha) subspace
57+
pp.kernel(subspace=['aa'])
58+
# solve ppRPA in the (alpha beta, alpha beta) subspace
59+
pp.kernel(subspace=['ab'])
60+
# solve ppRPA in the (beta beta, beta beta) subspace
61+
pp.kernel(subspace=['bb'])
62+
pp.analyze()
63+

‎examples/pprpa/04-gamma_pprpa_excitation_energy.py‎

Lines changed: 166 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,166 @@
1+
import os
2+
3+
from pyscf.pbc import df, dft, gto, lib
4+
5+
# For ppRPA excitation energy of N-electron system in particle-particle channel
6+
# mean field is (N-2)-electron
7+
8+
# carbon vacancy in diamond
9+
# see Table.1 in https://doi.org/10.1021/acs.jctc.4c00829
10+
cell = gto.Cell()
11+
cell.build(unit='angstrom',
12+
a=[[7.136000, 0.000000, 0.000000],
13+
[0.000000, 7.136000, 0.000000],
14+
[0.000000, 0.000000, 7.136000]],
15+
atom=[
16+
["C", (0.001233209267, 0.001233209267, 1.777383281725)],
17+
["C", (0.012225089066, 1.794731051862, 3.561871956432)],
18+
["C", (1.782392880032, 1.782392880032, 5.362763822515)],
19+
["C", (1.780552680464, -0.013167298675, -0.008776569968)],
20+
["C", (3.557268948138, 0.012225089066, 1.790128043568)],
21+
["C", (5.339774910934, 1.794731051862, 1.790128043568)],
22+
["C", (-0.013167298675, 1.780552680464, -0.008776569968)],
23+
["C", (1.782392880032, 3.569607119968, -0.010763822515)],
24+
["C", (1.794731051862, 5.339774910934, 1.790128043568)],
25+
["C", (2.671194343578, 0.900046784093, 0.881569117555)],
26+
["C", (0.887735886768, 0.887735886768, 6.244403380954)],
27+
["C", (0.900046784093, 2.680805656422, 4.470430882445)],
28+
["C", (2.680805656422, 4.451953215907, 0.881569117555)],
29+
["C", (0.895786995277, 6.238213500378, 0.896853305635)],
30+
["C", (0.930821705350, 0.930821705350, 2.673641624787)],
31+
["C", (4.421178294650, 0.930821705350, 2.678358375213)],
32+
["C", (0.900046784093, 2.671194343578, 0.881569117555)],
33+
["C", (6.238213500378, 0.895786995277, 0.896853305635)],
34+
["C", (2.680805656422, 0.900046784093, 4.470430882445)],
35+
["C", (4.451953215907, 2.680805656422, 0.881569117555)],
36+
["C", (0.930821705350, 4.421178294650, 2.678358375213)],
37+
["C", (1.794731051862, 0.012225089066, 3.561871956432)],
38+
["C", (0.012225089066, 3.557268948138, 1.790128043568)],
39+
["C", (3.569607119968, 1.782392880032, -0.010763822515)],
40+
["C", (1.736746319267, 1.736746319267, 1.671367479693)],
41+
["C", (5.351404126874, 0.000595873126, 0.004129648157)],
42+
["C", (0.000595873126, 5.351404126874, 0.004129648157)],
43+
["C", (2.676000000000, 2.676000000000, 6.244000000000)],
44+
["C", (6.244000000000, 2.676000000000, 2.676000000000)],
45+
["C", (2.676000000000, 6.244000000000, 2.676000000000)],
46+
["C", (0.000595873126, 0.000595873126, 5.347870351843)],
47+
["C", (0.001233209267, 5.350766790733, 3.574616718275)],
48+
["C", (1.780552680464, 5.365167298675, 5.360776569968)],
49+
["C", (3.571447319536, -0.013167298675, 5.360776569968)],
50+
["C", (5.365167298675, 1.780552680464, 5.360776569968)],
51+
["C", (5.365167298675, 3.571447319536, -0.008776569968)],
52+
["C", (5.350766790733, 5.350766790733, 1.777383281725)],
53+
["C", (4.464264113232, 0.887735886768, 6.243596619046)],
54+
["C", (4.451953215907, 2.671194343578, 4.470430882445)],
55+
["C", (2.671194343578, 4.451953215907, 4.470430882445)],
56+
["C", (0.895786995277, 6.249786499622, 4.455146694365)],
57+
["C", (4.421178294650, 4.421178294650, 2.673641624787)],
58+
["C", (6.249786499622, 4.456213004723, 0.896853305635)],
59+
["C", (0.887735886768, 4.464264113232, 6.243596619046)],
60+
["C", (6.249786499622, 0.895786995277, 4.455146694365)],
61+
["C", (4.456213004723, 6.249786499622, 0.896853305635)],
62+
["C", (5.350766790733, 0.001233209267, 3.574616718275)],
63+
["C", (-0.013167298675, 3.571447319536, 5.360776569968)],
64+
["C", (3.571447319536, 5.365167298675, -0.008776569968)],
65+
["C", (3.615253680733, 1.736746319267, 3.680632520307)],
66+
["C", (3.615253680733, 3.615253680733, 1.671367479693)],
67+
["C", (6.244000000000, 2.676000000000, 6.244000000000)],
68+
["C", (6.244000000000, 6.244000000000, 2.676000000000)],
69+
["C", (1.736746319267, 3.615253680733, 3.680632520307)],
70+
["C", (2.676000000000, 6.244000000000, 6.244000000000)],
71+
["C", (3.557268948138, 5.339774910934, 3.561871956432)],
72+
["C", (5.351404126874, 5.351404126874, 5.347870351843)],
73+
["C", (3.569607119968, 3.569607119968, 5.362763822515)],
74+
["C", (5.339774910934, 3.557268948138, 3.561871956432)],
75+
["C", (4.464264113232, 4.464264113232, 6.244403380954)],
76+
["C", (4.456213004723, 6.238213500378, 4.455146694365)],
77+
["C", (6.238213500378, 4.456213004723, 4.455146694365)],
78+
["C", (6.244000000000, 6.244000000000, 6.244000000000)],
79+
],
80+
dimension=3,
81+
max_memory=90000,
82+
verbose=5,
83+
basis='cc-pvdz',
84+
# create a (N-2)-electron system
85+
charge=2,
86+
precision=1e-12)
87+
88+
gdf = df.RSDF(cell)
89+
gdf.auxbasis = "cc-pvdz-ri"
90+
gdf_fname = 'gdf_ints.h5'
91+
gdf._cderi_to_save = gdf_fname
92+
if not os.path.isfile(gdf_fname):
93+
gdf.build()
94+
95+
# =====> Part I. Restricted ppRPA <=====
96+
# After SCF, PySCF might fail in Makov-Payne correction
97+
# save chkfile to restart
98+
chkfname = 'scf.chk'
99+
if os.path.isfile(chkfname):
100+
kmf = dft.RKS(cell).rs_density_fit()
101+
kmf.xc = "b3lyp"
102+
kmf.exxdiv = None
103+
kmf.with_df = gdf
104+
kmf.with_df._cderi = gdf_fname
105+
data = lib.chkfile.load(chkfname, 'scf')
106+
kmf.__dict__.update(data)
107+
else:
108+
kmf = dft.RKS(cell).rs_density_fit()
109+
kmf.xc = "b3lyp"
110+
kmf.exxdiv = None
111+
kmf.with_df = gdf
112+
kmf.with_df._cderi = gdf_fname
113+
kmf.chkfile = chkfname
114+
kmf.kernel()
115+
116+
# direct diagonalization, N6 scaling
117+
# ppRPA can be solved in an active space
118+
from pyscf.pprpa.rpprpa_direct import RppRPADirect
119+
pp = RppRPADirect(kmf, nocc_act=50, nvir_act=50)
120+
# number of two-electron addition states to print
121+
pp.pp_state = 50
122+
# solve for singlet states
123+
pp.kernel("s")
124+
# solve for triplet states
125+
pp.kernel("t")
126+
pp.analyze()
127+
128+
# Davidson algorithm, N4 scaling
129+
# ppRPA can be solved in an active space
130+
from pyscf.pprpa.rpprpa_davidson import RppRPADavidson
131+
pp = RppRPADavidson(kmf, nocc_act=50, nvir_act=50, nroot=50)
132+
# solve for singlet states
133+
pp.kernel("s")
134+
# solve for triplet states
135+
pp.kernel("t")
136+
pp.analyze()
137+
138+
# =====> Part II. Unrestricted ppRPA <=====
139+
chkfname = 'uscf.chk'
140+
if os.path.isfile(chkfname):
141+
kmf = dft.UKS(cell).rs_density_fit()
142+
kmf.xc = "b3lyp"
143+
kmf.exxdiv = None
144+
kmf.with_df = gdf
145+
kmf.with_df._cderi = gdf_fname
146+
data = lib.chkfile.load(chkfname, 'scf')
147+
kmf.__dict__.update(data)
148+
else:
149+
kmf = dft.UKS(cell).rs_density_fit()
150+
kmf.xc = "b3lyp"
151+
kmf.exxdiv = None
152+
kmf.with_df = gdf
153+
kmf.with_df._cderi = gdf_fname
154+
kmf.chkfile = chkfname
155+
kmf.kernel()
156+
157+
# direct diagonalization, N6 scaling
158+
# ppRPA can be solved in an active space
159+
from pyscf.pprpa.upprpa_direct import UppRPADirect
160+
pp = UppRPADirect(kmf, nocc_act=50, nvir_act=50)
161+
# number of two-electron addition states to print
162+
pp.pp_state = 50
163+
# solve ppRPA
164+
pp.kernel()
165+
pp.analyze()
166+

0 commit comments

Comments
 (0)