Skip to content

Commit 23268a1

Browse files
committed
remove cc
1 parent e418cd1 commit 23268a1

10 files changed

Lines changed: 93 additions & 53 deletions

File tree

‎examples/dh/02-xccode.py‎

Lines changed: 16 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,7 @@
22

33
from pyscf import gto
44
from pyscf.dh import DFDH
5-
from pyscf.dh.xccode import XCDH, XCType, parse_xc_dh, xc_equal
5+
from pyscf.dh.xccode import parse_xc_dh, xc_equal
66

77
mol = gto.M(atom="O; H 1 0.94; H 1 0.94 2 104.5", basis="cc-pVDZ")
88

@@ -17,17 +17,24 @@
1717
print(f"custom XYG3 (2-tuple): {mf.e_tot:.10f}")
1818

1919
# 3. Inspect what the parser returns for any functional
20-
print("\nFunc xc_scf cc c_os c_ss")
20+
print("\nFunc xc_scf c_os c_ss")
2121
for name in ["B2PLYP", "XYG3", "PBE0-DH", "DSD-PBEP86-D3BJ", "SCS-MP2"]:
2222
xc_list, _ = parse_xc_dh(name)
23-
print(f"{name:10s} {xc_list[0]:35s} {xc_list[2]:5g} {xc_list[3]:5g} {xc_list[4]:5g}")
23+
print(f"{name:10s} {xc_list[0]:35s} {xc_list[2]:5g} {xc_list[3]:5g}")
2424

25-
# 4. xDH: separate SCF and non-consistent energy via XCDH
26-
xcdh = XCDH("XYG3")
27-
print(f"\nXYG3 SCF = {xcdh.xc_scf.token}")
28-
low_rung = xcdh.xc_eng.remove(
29-
xcdh.xc_eng.extract_by_xctype(XCType.MP2 | XCType.DFTD3), inplace=False)
30-
print(f"XYG3 xc_n = {low_rung.token} (non-consistent energy)")
25+
# 4. All components returned by parse_xc_dh
26+
for name in ["MP2", "B2PLYP", "XYG3", "DSD-PBEP86-D3BJ", "B2PLYP-D3BJ"]:
27+
(xc, xc_n, c_os, c_ss), xc_add = parse_xc_dh(name)
28+
d3 = xc_add.get("D3", None)
29+
d4 = xc_add.get("D4", None)
30+
print(f"\n{name}:")
31+
print(f" xc = {xc}")
32+
print(f" xc_n = {xc_n}")
33+
print(f" c_os = {c_os:.6g} c_ss = {c_ss:.6g}")
34+
if d3:
35+
print(f" D3 = {d3}")
36+
if d4:
37+
print(f" D4 = {d4}")
3138

3239
# 5. XC token comparison via xc_equal
3340
print(f"\nxc_equal('HF', 'HF,') = {xc_equal('HF', 'HF,')}")

‎pyscf/dh/dh.py‎

Lines changed: 12 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -38,11 +38,20 @@ def __init__(self,
3838
elif isinstance(xc, (tuple, list)) and len(xc) == 2 and isinstance(xc[0], str):
3939
xc_list, xc_add = parse_xc_dh(xc)
4040
elif len(xc) == 5:
41+
import warnings
42+
warnings.warn(
43+
"5-tuple XC format is deprecated. "
44+
"Use a string name, a code string, or a 2-tuple (code_scf, code_eng) for xDH.",
45+
FutureWarning
46+
)
47+
xc_list = xc[0], xc[1], xc[3], xc[4]
48+
xc_add = {}
49+
elif len(xc) == 4:
4150
xc_list = xc
4251
xc_add = {}
4352
else:
4453
xc_list, xc_add = xc
45-
self.xc, self.xc_n, self.cc, self.c_os, self.c_ss = xc_list
54+
self.xc, self.xc_n, self.c_os, self.c_ss = xc_list
4655
self.xc_add = xc_add
4756
if self._scf is not None:
4857
if not hasattr(self._scf, 'xc'):
@@ -97,11 +106,11 @@ def calc_batch_size(self, unit_flop, pre_flop=0, fixed_mem=None):
97106

98107
@property
99108
def eval_ss(self):
100-
return abs(self.cc * self.c_ss) > 1e-7
109+
return abs(self.c_ss) > 1e-7
101110

102111
@property
103112
def eval_os(self):
104-
return abs(self.cc * self.c_os) > 1e-7
113+
return abs(self.c_os) > 1e-7
105114

106115
@property
107116
def eval_pt2(self):

‎pyscf/dh/dhutil.py‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -191,9 +191,9 @@ def tot_size(*args):
191191
return size
192192

193193

194-
def restricted_biorthogonalize(t_ijab, cc, c_os, c_ss):
195-
coef_0 = cc * (c_os + c_ss)
196-
coef_1 = - cc * c_ss
194+
def restricted_biorthogonalize(t_ijab, c_os, c_ss):
195+
coef_0 = c_os + c_ss
196+
coef_1 = -c_ss
197197
if abs(coef_1) < 1e-7:
198198
return coef_0 * t_ijab
199199
t_shape = t_ijab.shape

‎pyscf/dh/polar/rdfdh.py‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -201,11 +201,11 @@ def prepare_pt2_deriv(self):
201201
pdA_t_ijab += einsum("Aca, ijcb -> Aijab", pdA_F_0_mo[:, sv, sv], t_ijab)
202202
pdA_t_ijab /= D_ijab
203203

204-
cc, c_os, c_ss = self.cc, self.c_os, self.c_ss
205-
# T_ijab = cc * ((c_os + c_ss) * t_ijab - c_ss * t_ijab.swapaxes(-1, -2))
206-
# pdA_T_ijab = cc * ((c_os + c_ss) * pdA_t_ijab - c_ss * pdA_t_ijab.swapaxes(-1, -2))
207-
T_ijab = restricted_biorthogonalize(t_ijab, cc, c_os, c_ss)
208-
pdA_T_ijab = restricted_biorthogonalize(pdA_t_ijab, cc, c_os, c_ss)
204+
c_os, c_ss = self.c_os, self.c_ss
205+
# T_ijab = (c_os + c_ss) * t_ijab - c_ss * t_ijab.swapaxes(-1, -2)
206+
# pdA_T_ijab = (c_os + c_ss) * pdA_t_ijab - c_ss * pdA_t_ijab.swapaxes(-1, -2)
207+
T_ijab = restricted_biorthogonalize(t_ijab, c_os, c_ss)
208+
pdA_T_ijab = restricted_biorthogonalize(pdA_t_ijab, c_os, c_ss)
209209

210210
pdA_G_ia_ri[:, :, sI] += einsum("Aijab, Pjb -> APia", pdA_T_ijab, Y_ia_ri)
211211
pdA_G_ia_ri[:, :, sI] += einsum("ijab, APjb -> APia", T_ijab, pdA_Y_ia_ri)

‎pyscf/dh/polar/udfdh.py‎

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -209,7 +209,7 @@ def prepare_pdA_Y_ia_ri(self):
209209

210210
def prepare_pt2_deriv(self):
211211
tensors = self.tensors
212-
cc, c_os, c_ss = self.cc, self.c_os, self.c_ss
212+
c_os, c_ss = self.c_os, self.c_ss
213213
nocc, nvir, nmo, naux = self.nocc, self.nvir, self.nmo, self.df_ri.get_naoaux()
214214
mocc, mvir = max(nocc), max(nvir)
215215
so, sv = self.so, self.sv
@@ -248,21 +248,21 @@ def prepare_pt2_deriv(self):
248248
if σς in (αα, ββ):
249249
# T_ijab = cc * 0.5 * c_ss * (t_ijab - t_ijab.swapaxes(-1, -2))
250250
# pdA_T_ijab = cc * 0.5 * c_ss * (pdA_t_ijab - pdA_t_ijab.swapaxes(-1, -2))
251-
T_ijab = cc * 0.5 * c_ss * hermi_sum_last2dim(t_ijab, hermi=ANTIHERMI, inplace=False)
252-
pdA_T_ijab = cc * 0.5 * c_ss * hermi_sum_last2dim(pdA_t_ijab, hermi=ANTIHERMI, inplace=False)
251+
T_ijab = 0.5 * c_ss * hermi_sum_last2dim(t_ijab, hermi=ANTIHERMI, inplace=False)
252+
pdA_T_ijab = 0.5 * c_ss * hermi_sum_last2dim(pdA_t_ijab, hermi=ANTIHERMI, inplace=False)
253253
pdA_D_rdm1[σ][:, so[σ], so[σ]] -= 2 * einsum("kiba, Akjba -> Aij", T_ijab, pdA_t_ijab)
254254
pdA_D_rdm1[σ][:, sv[σ], sv[σ]] += 2 * einsum("ijac, Aijbc -> Aab", T_ijab, pdA_t_ijab)
255255
pdA_G_ia_ri[σ][:, :, sI] += 4 * einsum("ijab, APjb -> APia", T_ijab, pdA_Y_ia_ri[σ])
256256
pdA_G_ia_ri[σ][:, :, sI] += 4 * einsum("Aijab, Pjb -> APia", pdA_T_ijab, Y_ia_ri[σ])
257257
else:
258-
T_ijab = cc * c_os * t_ijab
259-
pdA_T_ijab = cc * c_os * pdA_t_ijab
258+
T_ijab = c_os * t_ijab
259+
pdA_T_ijab = c_os * pdA_t_ijab
260260
for sJ in gen_batch(0, nocc[α], nbatch):
261261
if sI == sJ:
262262
T_jkab = T_ijab
263263
else:
264264
t_jkab = tensors["t_ijab" + str(αβ)][sJ]
265-
T_jkab = cc * c_os * t_jkab
265+
T_jkab = c_os * t_jkab
266266
pdA_D_rdm1[α][:, sI, sJ] -= einsum("jkba, Aikba -> Aij", T_jkab, pdA_t_ijab)
267267
# pdA_D_rdm1[α][:, so[α], so[α]] -= einsum("ikab, Ajkab -> Aij", T_ijab, pdA_t_ijab)
268268
pdA_D_rdm1[β][:, so[β], so[β]] -= einsum("kiba, Akjba -> Aij", T_ijab, pdA_t_ijab)

‎pyscf/dh/rdfdh.py‎

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -59,13 +59,13 @@ def energy_elec_nc(mf, mo_coeff=None, h1e=None, vhf=None, **_):
5959
def energy_elec_pt2(mf, params=None, eng_bi=None, **kwargs):
6060
if not mf.eval_pt2:
6161
return 0, 0, 0
62-
cc, c_os, c_ss = params if params else mf.cc, mf.c_os, mf.c_ss
62+
c_os, c_ss = params if params else mf.c_os, mf.c_ss
6363
emp2_0, eng_bi1, eng_bi2 = eng_bi if eng_bi else mf.energy_elec_mp2(eval_ss=mf.eval_ss, **kwargs)
6464
if getattr(mf, 'mp2_backend', None) == "dfmp2_native":
65-
return cc * emp2_0, None, None
65+
return emp2_0, None, None
6666
if getattr(mf, 'mp2_backend', None) == "dfmp2":
67-
return cc * (c_os * eng_bi1 + c_ss * eng_bi2), eng_bi1, eng_bi2
68-
return (cc * ((c_os + c_ss) * eng_bi1 - c_ss * eng_bi2),
67+
return c_os * eng_bi1 + c_ss * eng_bi2, eng_bi1, eng_bi2
68+
return (((c_os + c_ss) * eng_bi1 - c_ss * eng_bi2),
6969
eng_bi1,
7070
eng_bi1 - eng_bi2)
7171

@@ -344,7 +344,7 @@ def prepare_pt2(self, dump_t_ijab=True):
344344
e = self.mo_energy
345345
naux = self.df_ri.get_naoaux()
346346
so, sv = self.so, self.sv
347-
cc, c_os, c_ss = self.cc, self.c_os, self.c_ss
347+
c_os, c_ss = self.c_os, self.c_ss
348348

349349
D_rdm1 = np.zeros((nmo, nmo))
350350

@@ -371,7 +371,7 @@ def build(sI, t_ijab, g_ijab):
371371
eng_bi2[0] += einsum("ijab, ijba ->", t_ijab, g_ijab)
372372
if dump_t_ijab:
373373
tensors["t_ijab"][sI] = t_ijab
374-
T_ijab = restricted_biorthogonalize(t_ijab, cc, c_os, c_ss)
374+
T_ijab = restricted_biorthogonalize(t_ijab, c_os, c_ss)
375375
D_rdm1[sv, sv] += 2 * einsum("ijac, ijbc -> ab", T_ijab, t_ijab)
376376
D_rdm1[so, so] -= 2 * einsum("ijab, ikab -> jk", T_ijab, t_ijab)
377377
G_ia_ri[:, sI] = einsum("ijab, Pjb -> Pia", T_ijab, Y_ia_ri)

‎pyscf/dh/test/test_rdfdh.py‎

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -41,3 +41,24 @@ def test_XYG3(self):
4141
mol = gto.Mole(atom="O; H 1 2; H 1 2 2 104.2458898548", unit="AU", basis="cc-pVTZ").build()
4242
mf = dh.DFDH(mol, xc="XYG3", mp2_backend="dfmp2", frozen=1).run()
4343
self.assertAlmostEqual(mf.e_tot, REF_ETOT, places=4)
44+
45+
def test_dsdpbep86_d3bj(self):
46+
"""Cross-validation against REST (Rust Electronic-Structure Toolkit)."""
47+
from pyscf.dh.xccode import parse_xc_dh
48+
(xc, xc_n, c_os, c_ss), xc_add = parse_xc_dh("DSD-PBEP86-D3BJ")
49+
self.assertEqual(xc, "0.69*HF + 0.31*PBE, 0.44*P86")
50+
self.assertAlmostEqual(c_os, 0.52, places=6)
51+
self.assertAlmostEqual(c_ss, 0.22, places=6)
52+
self.assertEqual(xc_add["D3"]["version"], "d3bj")
53+
54+
REF_SCF = -56.3349732703
55+
REF_ETOT = -56.4305347110
56+
mol = gto.Mole(atom="""
57+
N -2.1988391019 1.8973746268 0.0000000000
58+
H -1.1788391019 1.8973746268 0.0000000000
59+
H -2.5388353987 1.0925460144 -0.5263586446
60+
H -2.5388400276 2.7556271745 -0.4338224694
61+
""", basis="cc-pVDZ", unit="Angstrom", verbose=0).build()
62+
mf = dh.DFDH(mol, xc="DSD-PBEP86-D3BJ", auxbasis_jk="def2-universal-jkfit", auxbasis_ri="def2-universal-jkfit").run()
63+
self.assertAlmostEqual(mf.mf_s.e_tot, REF_SCF, places=3)
64+
self.assertAlmostEqual(mf.e_tot, REF_ETOT, places=6)

‎pyscf/dh/udfdh.py‎

Lines changed: 8 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -42,24 +42,22 @@ def energy_elec_nc(mf, mo_coeff=None, h1e=None, vhf=None, **_):
4242

4343

4444
def energy_elec_pt2(mf, params=None, eng_bi=None, **kwargs):
45-
cc, c_os, c_ss = params if params else mf.cc, mf.c_os, mf.c_ss
46-
if not mf.eval_pt2:
47-
return 0, 0, 0
45+
c_os, c_ss = params if params else mf.c_os, mf.c_ss
4846
emp2_0, eng_bi1, eng_bi2 = eng_bi if eng_bi else mf.energy_elec_mp2(eval_ss=mf.eval_ss, **kwargs)
4947
if getattr(mf, 'mp2_backend', None) == "dfmp2_native":
50-
return cc * emp2_0, None, None
48+
return emp2_0, None, None
5149
if getattr(mf, 'mp2_backend', None) == "dfmp2":
52-
return cc * (c_os * eng_bi1 + c_ss * eng_bi2), eng_bi1, eng_bi2
50+
return c_os * eng_bi1 + c_ss * eng_bi2, eng_bi1, eng_bi2
5351
eng_os = eng_bi1[αβ]
5452
eng_ss = 0.5 * (eng_bi1[αα] + eng_bi1[ββ] - eng_bi2[αα] - eng_bi2[ββ])
55-
eng_pt2 = cc * (c_os * eng_os + c_ss * eng_ss)
53+
eng_pt2 = c_os * eng_os + c_ss * eng_ss
5654
return eng_pt2, eng_os, eng_ss
5755

5856

5957
def energy_elec(mf, params=None, **kwargs):
6058
eng_nc = mf.energy_elec_nc(**kwargs)[0]
6159
nocc, nvir = mf.nocc, mf.nvir
62-
_, _, c_ss = params if params else mf.cc, mf.c_os, mf.c_ss
60+
_, c_ss = params if params else mf.c_os, mf.c_ss
6361
eval_ss = True if abs(c_ss) > 1e-7 else False
6462
t_ijab_blk = None
6563
if mf.with_t_ijab:
@@ -337,7 +335,7 @@ def prepare_pt2(self, dump_t_ijab=True, fast_trans=True):
337335
eo, ev = self.eo, self.ev
338336
naux = self.df_ri.get_naoaux()
339337
so, sv = self.so, self.sv
340-
cc, c_os, c_ss = self.cc, self.c_os, self.c_ss
338+
c_os, c_ss = self.c_os, self.c_ss
341339
eval_ss = True if abs(c_ss) > 1e-7 else False
342340

343341
D_rdm1 = np.zeros((2, nmo, nmo))
@@ -377,12 +375,12 @@ def prepare_pt2(self, dump_t_ijab=True, fast_trans=True):
377375
t_ijab = tensors["t_ijab" + str(σς)][sI]
378376
if σς in (αα, ββ):
379377
# T_ijab = cc * 0.5 * c_ss * (t_ijab - t_ijab.swapaxes(-1, -2))
380-
T_ijab = cc * 0.5 * c_ss * hermi_sum_last2dim(t_ijab, hermi=ANTIHERMI, inplace=False)
378+
T_ijab = 0.5 * c_ss * hermi_sum_last2dim(t_ijab, hermi=ANTIHERMI, inplace=False)
381379
D_rdm1[σ, so[σ], so[σ]] -= 2 * einsum("kiab, kjab -> ij", T_ijab, t_ijab)
382380
D_rdm1[σ, sv[σ], sv[σ]] += 2 * einsum("ijac, ijbc -> ab", T_ijab, t_ijab)
383381
G_ia_ri[σ][:, sI] += 4 * einsum("ijab, Pjb -> Pia", T_ijab, Y_ia_ri[σ])
384382
else: # σς == αβ
385-
T_ijab = cc * c_os * t_ijab
383+
T_ijab = c_os * t_ijab
386384
# D_rdm1[α, so[α], so[α]] -= einsum("ikab, jkab -> ij", T_ijab, t_ijab)
387385
# D_rdm1[β, so[β], so[β]] -= einsum("kiba, kjba -> ij", T_ijab, t_ijab)
388386
# D_rdm1[α, sv[α], sv[α]] += einsum("ijac, ijbc -> ab", T_ijab, t_ijab)

‎pyscf/dh/xccode/functionals/family_bdh.json‎

Lines changed: 9 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -45,14 +45,13 @@
4545
"tested": "QChem 5.1.1"
4646
},
4747
"DSD_PBEP86": {
48-
"code": "0.72*HF + 0.28*PBE, 0.44*P86_FT + 0.51*MP2_OS + 0.36*MP2_SS",
48+
"code": "0.72*HF + 0.28*PBE, 0.44*P86 + 0.51*MP2_OS + 0.36*MP2_SS",
4949
"ref": "10.1002/jcc.23391"
5050
},
5151
"DSD_PBEP86_D3BJ": {
52-
"code": "0.69*HF + 0.31*PBE, 0.44*P86_FT + 0.52*MP2_OS + 0.22*MP2_SS + DFTD3(BJ, s6=0.48, a1=0, s8=0, a2=5.6)",
52+
"code": "0.69*HF + 0.31*PBE, 0.44*P86 + 0.52*MP2_OS + 0.22*MP2_SS + DFTD3(BJ, xc=dsdpbep86)",
5353
"ref": "10.1002/jcc.23391",
54-
"tested": "QChem 5.1.1",
55-
"notes": "P86_FT used in Gaussian, QChem, etc instead of P86. Also see https://github.com/pyscf/pyscf/issues/1004."
54+
"tested": "MRCC 2022-03-18"
5655
},
5756
"DSD_PBEPBE_D3BJ": {
5857
"code": "0.68*HF + 0.32*PBE, 0.49*PBE + 0.55*MP2_OS + 0.13*MP2_SS + DFTD3(BJ, s6=0.78, a1=0, s8=0, a2=6.1)",
@@ -113,6 +112,12 @@
113112
"PWPB95"
114113
]
115114
},
115+
"QCHEM_DSD_PBEP86_D3BJ": {
116+
"code": "0.69*HF + 0.31*PBE, 0.44*P86_FT + 0.52*MP2_OS + 0.22*MP2_SS + DFTD3(BJ, xc=dsdpbep86)",
117+
"ref": "10.1002/jcc.23391",
118+
"tested": "QChem 5.1.1",
119+
"notes": "P86_FT used in Gaussian/QChem instead of P86; see https://github.com/pyscf/pyscf/issues/1004"
120+
},
116121
"REV_DSD_PBEP86_D3BJ": {
117122
"code": "0.69*HF + 0.31*PBE, 0.4296*P86_FT + 0.5785*MP2_OS + 0.0799*MP2_SS + DFTD3(BJ, S6=0.4377, A1=0, S8=0, A2=5.5)",
118123
"ref": "10.1002/jcc.23391",

‎pyscf/dh/xccode/parse.py‎

Lines changed: 7 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -42,9 +42,9 @@ def _extract_components(xc_eng, xc_scf, name):
4242
mp2_list = xc_eng.extract_by_xctype(XCType.MP2)
4343
if len(mp2_list) > 0:
4444
mp2 = mp2_list[0]
45-
cc, c_os, c_ss = mp2.fac, mp2.parameters[0], mp2.parameters[1]
45+
c_os, c_ss = mp2.fac * mp2.parameters[0], mp2.fac * mp2.parameters[1]
4646
else:
47-
cc, c_os, c_ss = 0, 0, 0
47+
c_os, c_ss = 0, 0
4848
xc_add = {}
4949
d3_list = xc_eng.extract_by_xctype(XCType.DFTD3)
5050
if len(d3_list) > 0:
@@ -62,19 +62,19 @@ def _extract_components(xc_eng, xc_scf, name):
6262
d4 = d4_list[0]
6363
d4_xc = d4.additional.get("XC", _strip_d3_suffix(name))
6464
xc_add["D4"] = {"xc": d4_xc, "version": "d4"}
65-
return xc, xc_n, cc, c_os, c_ss, xc_add
65+
return xc, xc_n, c_os, c_ss, xc_add
6666

6767

6868
def parse_xc_dh(xc_dh: str):
6969
if isinstance(xc_dh, tuple) and len(xc_dh) == 2:
7070
xc_scf = XCList(xc_dh[0], code_scf=True).token
7171
xc_eng = XCList(xc_dh[1], code_scf=False)
72-
xc, xc_n, cc, c_os, c_ss, xc_add = _extract_components(xc_eng, xc_scf, "")
73-
return (xc, xc_n, cc, c_os, c_ss), xc_add
72+
xc, xc_n, c_os, c_ss, xc_add = _extract_components(xc_eng, xc_scf, "")
73+
return (xc, xc_n, c_os, c_ss), xc_add
7474

7575
xcdh = XCDH(xc_dh)
7676
_check_unsupported(xc_dh)
7777
xc_scf = xcdh.xc_scf.token
7878
xc_eng = xcdh.xc_eng
79-
xc, xc_n, cc, c_os, c_ss, xc_add = _extract_components(xc_eng, xc_scf, xc_dh)
80-
return (xc, xc_n, cc, c_os, c_ss), xc_add
79+
xc, xc_n, c_os, c_ss, xc_add = _extract_components(xc_eng, xc_scf, xc_dh)
80+
return (xc, xc_n, c_os, c_ss), xc_add

0 commit comments

Comments
 (0)