Skip to content

Commit e5148d6

Browse files
dgasmithdgasmith
authored
Merge pull request #49 from tylertak/srimp2
Added sDF-MP2 files
2 parents 965f0c7 + bd411b0 commit e5148d6

File tree

3 files changed

+158
-0
lines changed

3 files changed

+158
-0
lines changed

‎Moller-Plesset/README.md‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,7 @@ ontop of this result.
99
### Included Reference Implementations
1010
- `MP2.py`: Second-order Moller-Plesset perturbation theory
1111
- `DF-MP2.py`: MP2 utilizing density-fitting to reduce its overall cost.
12+
- `sDF-MP2.py`: Stochastic orbital MP2 utilizing density-fitting (sDF-MP2) to further reduce the computational cost of DF-MP2.
1213
- `MP3.py`: Third-order Moller-Plesset perturbation theory
1314
- `MP3-SO.py`: MP3 in spin-orbital formalism that simplifies the equations and provides contrast to the spin-summed version above.
1415
- `MPn.py`: An example on how to automate the contraction of higher order MP theory.
@@ -20,3 +21,5 @@ ontop of this result.
2021
- [[Takasuka:2008:044112](http://aip.scitation.org/doi/10.1063/1.2958921)] A. Takatsuka, T. Siichiro, and W. Hackbusch, *J. Phys. Chem.*, **129**, 044112 (2008)
2122
3) Equations used in reference implementations:
2223
- [[Szabo:1996](https://books.google.com/books?id=KQ3DAgAAQBAJ&printsec=frontcover&dq=szabo+%26+ostlund&hl=en&sa=X&ved=0ahUKEwiYhv6A8YjUAhXLSCYKHdH5AJ4Q6AEIJjAA#v=onepage&q=szabo%20%26%20ostlund&f=false)] A. Szabo and N. S. Ostlund, *Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory.* Courier Corporation, 1996.
24+
4) sDF-MP2 reference:
25+
- [[Takeshita:2017:4605](https://pubs.acs.org/doi/abs/10.1021/acs.jctc.7b00343)] T.Y. Takeshita, W.A. de Jong, D. Neuhauser, R. Baer and E. Rabani

‎Moller-Plesset/sDF-MP2.py‎

Lines changed: 151 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,151 @@
1+
"""
2+
A reference implementation of stochastic orbital resolution of identity
3+
(or density-fitted) MP2 (sRI-MP2) from a RHF reference.
4+
5+
Reference:
6+
Stochastic Formulation of the Resolution of Identity: Application to
7+
Second Order Moller-Plesset Perturbation Theory
8+
J. Chem. Theory Comput., 2017, 13 (10), pp 4605-4610
9+
DOI: 10.1021/acs.jctc.7b00343
10+
11+
Tyler Y. Takeshita
12+
Department of Chemistry, University of California Berkeley
13+
Materials Sciences Division, Lawrence Berkeley National Laboratory
14+
15+
Wibe A. de Jong
16+
Computational Research Division, Lawrence Berkeley National Laboratory
17+
18+
Daniel Neuhauser
19+
Department of Chemistry and Biochemistry, University of California, Los Angeles
20+
21+
Roi Baer
22+
Fritz Harber Center for Molecular Dynamics, Institute of Chemistry, The Hebrew University of Jerusalem
23+
24+
Eran Rabani
25+
Department of Chemistry, University of California Berkeley
26+
Materials Sciences Division, Lawrence Berkeley National Laboratory
27+
The Sackler Center for Computational Molecular Science, Tel Aviv University
28+
"""
29+
30+
__authors__ = ["Tyler Y. Takeshita", "Daniel G. A. Smith"]
31+
__credits__ = ["Tyler Y. Takeshita", "Daniel G. A. Smith"]
32+
33+
__copyright__ = "(c) 2014-2017, The Psi4NumPy Developers"
34+
__license__ = "BSD-3-Clause"
35+
__date__ = "2018-04-14"
36+
37+
import numpy as np
38+
import psi4
39+
import time
40+
41+
# Set numpy defaults
42+
np.set_printoptions(precision=5, linewidth=200, suppress=True)
43+
# psi4.set_output_file("output.dat")
44+
45+
# ==> Geometry <==
46+
# Note: Symmetry was turned off
47+
mol = psi4.geometry("""
48+
O
49+
H 1 0.96
50+
H 1 0.96 2 104
51+
symmetry c1
52+
""")
53+
54+
# How many samples to run?
55+
nsample = 5000
56+
57+
# ==> Basis sets <==
58+
psi4.set_options({
59+
'basis': 'aug-cc-pvdz',
60+
'scf_type': 'df',
61+
'e_convergence': 1e-10,
62+
'd_convergence': 1e-10
63+
})
64+
65+
wfn = psi4.core.Wavefunction.build(mol, psi4.core.get_global_option('basis'))
66+
67+
# Build auxiliary basis set
68+
aux = psi4.core.BasisSet.build(mol, "DF_BASIS_SCF", "", "RIFIT", "aug-cc-pVDZ")
69+
70+
# Get orbital basis & build zero basis
71+
orb = wfn.basisset()
72+
73+
# The zero basis set
74+
zero_bas = psi4.core.BasisSet.zero_ao_basis_set()
75+
76+
# Build instance of MintsHelper
77+
mints = psi4.core.MintsHelper(orb)
78+
79+
# ==> Build Density-Fitted Integrals <==
80+
81+
# Build (P|pq) raw 3-index ERIs, dimension (1, Naux, nbf, nbf)
82+
Ppq = mints.ao_eri(zero_bas, aux, orb, orb)
83+
84+
# Build & invert Coulomb metric, dimension (1, Naux, 1, Naux)
85+
metric = mints.ao_eri(zero_bas, aux, zero_bas, aux)
86+
metric.power(-0.5, 1.e-14)
87+
88+
# Remove excess dimensions of Ppq, & metric
89+
Ppq = np.squeeze(Ppq)
90+
metric = np.squeeze(metric)
91+
92+
# Build the Qso object
93+
Qpq = np.einsum('QP,Ppq->Qpq', metric, Ppq)
94+
95+
# ==> Perform HF <==
96+
energy, scf_wfn = psi4.energy('scf', return_wfn=True)
97+
print("Finished SCF...")
98+
99+
# ==> AO -> MO Transformation of integrals <==
100+
# Get the MO coefficients and energies
101+
evecs = np.array(scf_wfn.epsilon_a())
102+
103+
Co = scf_wfn.Ca_subset("AO", "OCC")
104+
Cv = scf_wfn.Ca_subset("AO", "VIR")
105+
106+
Qia = np.einsum("Qpq,pi,qa->Qia", Qpq, Co, Cv, optimize=True)
107+
108+
nocc = scf_wfn.nalpha()
109+
nvirt = scf_wfn.nmo() - nocc
110+
111+
denom = 1.0 / (evecs[:nocc].reshape(-1, 1, 1, 1) + evecs[:nocc].reshape(-1, 1, 1) -
112+
evecs[nocc:].reshape(-1, 1) - evecs[nocc:])
113+
114+
t = time.time()
115+
e_srimp2 = 0.0
116+
print("Transformed ERIs...")
117+
118+
# Loop over samples to reduce stochastic noise
119+
print("Starting sample loop...")
120+
for x in range(nsample):
121+
122+
# ==> Build Stochastic Integral Matrices <==
123+
# Create two random vector
124+
vec = np.random.choice([-1, 1], size=(Qia.shape[0]))
125+
vecp = np.random.choice([-1, 1], size=(Qia.shape[0]))
126+
127+
# Generate first R matrices
128+
ia = np.einsum("Q,Qia->ia", vec, Qia)
129+
iap = np.einsum("Q,Qia->ia", vecp, Qia)
130+
131+
# ==> Calculate a single stochastic RI-MP2 (sRI-MP2) sample <==
132+
133+
# Caculate sRI-MP2 correlation energy
134+
e_srimp2 += 2.0 * np.einsum('ijab,ia,ia,jb,jb->', denom, ia, iap, ia, iap)
135+
e_srimp2 -= np.einsum('ijab,ia,ib,jb,ja->', denom, ia, iap, ia, iap)
136+
137+
e_srimp2 /= float(nsample)
138+
total_time = time.time() - t
139+
time_per_sample = total_time / float(nsample)
140+
141+
# Print sample energy to output
142+
print("\nNumber of samples: % 16d" % nsample)
143+
print("Total time (s): % 16.2f" % total_time)
144+
print("Time per sample (us): % 16.2f" % (time_per_sample * 1.e6))
145+
print("sRI-MP2 correlation sample energy: % 16.10f" % e_srimp2)
146+
147+
psi_mp2_energy = psi4.energy("MP2")
148+
mp2_correlation_energy = psi4.get_variable("MP2 CORRELATION ENERGY")
149+
150+
print("\nRI-MP2 energy: % 16.10f" % mp2_correlation_energy)
151+
print("Sample error % 16.10f" % (e_srimp2 - mp2_correlation_energy))

‎tests/test_RI_PT.py‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -13,6 +13,10 @@ def test_MP2(workspace):
1313
exe_py(workspace, tdir, 'MP2')
1414

1515

16+
def test_sDF_MP2(workspace):
17+
exe_py(workspace, tdir, 'sDF-MP2')
18+
19+
1620
def test_MP3(workspace):
1721
exe_py(workspace, tdir, 'MP3')
1822

0 commit comments

Comments
 (0)