Skip to content

Commit 5579b8e

Browse files
ethanvoethanvo
committed
Working CC2?
1 parent f87d60b commit 5579b8e

File tree

1 file changed

+172
-88
lines changed

1 file changed

+172
-88
lines changed

‎pyscf/cc/dfrcc2.py

Lines changed: 172 additions & 88 deletions
Original file line numberDiff line numberDiff line change
@@ -9,109 +9,126 @@
99
BLKMIN = getattr(__config__, 'cc_ccsd_blkmin', 4)
1010
MEMORYMIN = getattr(__config__, 'cc_ccsd_memorymin', 2000)
1111

12-
def update_t1(cc, t1, t2, eris):
12+
def energy(mycc, t1=None, eris=None):
13+
'''CCSD correlation energy'''
14+
if t1 is None: t1 = mycc.t1
15+
if eris is None: eris = mycc.ao2mo()
16+
17+
nocc, nvir = t1.shape
18+
fock = eris.fock
19+
mo_e_o = eris.mo_energy[:nocc]
20+
mo_e_v = eris.mo_energy[nocc:] + mycc.level_shift
21+
eia = mo_e_o[:,None] - mo_e_v
22+
Lov = eris.Lov
23+
Bov = make_Bov(mycc, t1, eris)
24+
e = np.einsum('ia,ia', fock[:nocc,nocc:], t1) * 2
25+
for i in range(nocc):
26+
t2i = lib.einsum('La,Ljb->jab', Bov[:, i, :], Bov)
27+
ejab = lib.direct_sum('a,jb->jab', eia[i, :], eia)
28+
t2i /= ejab
29+
taui = t2i + lib.einsum('a,jb->jab', t1[i, :], t1)
30+
e += 2*lib.einsum('jab,La,Ljb', taui, Lov[:, i, :], Lov)
31+
e -= lib.einsum('iab,Lia,Lb', taui, Lov, Lov[:, i, :])
32+
33+
if abs(e.imag) > 1e-4:
34+
logger.warn(mycc, 'Non-zero imaginary part found in CCSD energy %s', e)
35+
return e.real
36+
37+
def update_t1(cc, t1, eris):
1338
# Ref: Hirata et al., J. Chem. Phys. 120, 2581 (2004) Eqs.(35)-(36)
1439
assert(isinstance(eris, ccsd._ChemistsERIs))
1540
nocc, nvir = t1.shape
1641
fock = eris.fock
1742
mo_e_o = eris.mo_energy[:nocc]
1843
mo_e_v = eris.mo_energy[nocc:] + cc.level_shift
1944

45+
eia = mo_e_o[:,None] - mo_e_v
46+
47+
Loo = eris.Loo
48+
Lov = eris.Lov
49+
Lvv = eris.Lvv
50+
51+
Bov = make_Bov(cc, t1, eris)
52+
53+
foo = fock[:nocc,:nocc].copy()
2054
fov = fock[:nocc,nocc:].copy()
55+
fvv = fock[nocc:,nocc:].copy()
2156

22-
Foo = imd.cc_Foo(t1,t2,eris)
23-
Fvv = imd.cc_Fvv(t1,t2,eris)
24-
Fov = imd.cc_Fov(t1,t2,eris)
57+
Fki = np.zeros_like(foo)
58+
Fac = np.zeros_like(fvv)
59+
Fkc = np.zeros_like(fov)
60+
for i in range(nocc):
61+
t2i = lib.einsum('La,Ljb->jab', Bov[:, i, :], Bov)
62+
ejab = lib.direct_sum('a,jb->jab', eia[i, :], eia)
63+
t2i /= ejab
64+
Fki[:, i] += 2*lib.einsum('Lkc,Lld,lcd->k', Lov, Lov, t2i)
65+
Fki[:, i] -= lib.einsum('Lkd,Llc,lcd->k', Lov, Lov, t2i)
66+
Fac +=-2*lib.einsum('Lc,Lld,lad->ac', Lov[:, i, :], Lov, t2i)
67+
Fac += lib.einsum('Ld,Llc,lad->ac', Lov[:, i, :], Lov, t2i)
68+
Fki += 2*lib.einsum('Lkc,Lld,ic,ld->ki', Lov, Lov, t1, t1)
69+
Fki -= lib.einsum('Lkd,Llc,ic,ld->ki', Lov, Lov, t1, t1)
70+
Fac -= 2*lib.einsum('Lkc,Lld,ka,ld->ac', Lov, Lov, t1, t1)
71+
Fac += lib.einsum('Lkd,Llc,ka,ld->ac', Lov, Lov, t1, t1)
72+
Fkc = 2*np.einsum('Lkc,Lld,ld->kc', Lov, Lov, t1)
73+
Fkc -= np.einsum('Lkd,Llc,ld->kc', Lov, Lov, t1)
74+
Fki += foo
75+
Fac += fvv
76+
Fkc += fov
77+
Foo = Fki
78+
Fvv = Fac
79+
Fov = Fkc
2580

2681
# Move energy terms to the other side
2782
Foo[np.diag_indices(nocc)] -= mo_e_o
2883
Fvv[np.diag_indices(nvir)] -= mo_e_v
2984

30-
Loo = eris.Loo
31-
Lov = eris.Lov
32-
Lvv = eris.Lvv
33-
34-
eia = mo_e_o[:,None] - mo_e_v
35-
3685
# T1 equation
3786
t1new =-2*np.einsum('kc,ka,ic->ia', fov, t1, t1)
3887
t1new += np.einsum('ac,ic->ia', Fvv, t1)
3988
t1new += -np.einsum('ki,ka->ia', Foo, t1)
40-
t1new += 2*np.einsum('kc,kica->ia', Fov, t2)
41-
t1new += -np.einsum('kc,ikca->ia', Fov, t2)
89+
for i in range(nocc):
90+
t2i = lib.einsum('La,Ljb->jab', Bov[:, i, :], Bov)
91+
ejab = lib.direct_sum('a,jb->jab', eia[i, :], eia)
92+
t2i /= ejab
93+
t1new += 2*lib.einsum('c,ica->ia', Fov[i, :], t2i)
94+
t1new[i, :] += -lib.einsum('kc,kca->a', Fov, t2i)
95+
t1new[i, :] += 2*lib.einsum('Lkd,Lac,kcd->a', Lov, Lvv, t2i)
96+
t1new[i, :] += -lib.einsum('Lkc,Lad,kcd->a', Lov, Lvv, t2i)
97+
t1new +=-2*lib.einsum('Llc,Li,lac->ia', Lov, Loo[:, i, :], t2i)
98+
t1new += lib.einsum('Lc,Lli,lac->ia', Lov[:, i, :], Loo, t2i)
99+
42100
t1new += np.einsum('kc,ic,ka->ia', Fov, t1, t1)
43101
t1new += fov.conj()
44102
t1new += 2*np.einsum('Lkc,Lia,kc->ia', Lov, Lov, t1)
45103
t1new += -np.einsum('Lki,Lac,kc->ia', Loo, Lvv, t1)
46-
t1new += 2*lib.einsum('Lkd,Lac,ikcd->ia', Lov, Lvv, t2)
47-
t1new += -lib.einsum('Lkc,Lad,ikcd->ia', Lov, Lvv, t2)
48104
t1new += 2*lib.einsum('Lkd,Lac,kd,ic->ia', Lov, Lvv, t1, t1)
49105
t1new += -lib.einsum('Lkc,Lad,kd,ic->ia', Lov, Lvv, t1, t1)
50-
t1new +=-2*lib.einsum('Llc,Lki,klac->ia', Lov, Loo, t2)
51-
t1new += lib.einsum('Lkc,Lli,klac->ia', Lov, Loo, t2)
52106
t1new +=-2*lib.einsum('Llc,Lki,lc,ka->ia', Lov, Loo, t1, t1)
53107
t1new += lib.einsum('Lkc,Lli,lc,ka->ia', Lov, Loo, t1, t1)
54108

55109
t1new /= eia
56110

57111
return t1new
58112

59-
def update_t2(cc, t1, t2, eris):
60-
# Ref: Hirata et al., J. Chem. Phys. 120, 2581 (2004) Eqs.(35)-(36)
61-
assert(isinstance(eris, ccsd._ChemistsERIs))
62-
nocc, nvir = t1.shape
63-
fock = eris.fock
64-
mo_e_o = eris.mo_energy[:nocc]
65-
mo_e_v = eris.mo_energy[nocc:] + cc.level_shift
66-
67-
fov = fock[:nocc,nocc:].copy()
68-
foo = fock[:nocc,:nocc].copy()
69-
fvv = fock[nocc:,nocc:].copy()
70-
71-
Foo = imd.cc_Foo(t1,t2,eris)
72-
Fvv = imd.cc_Fvv(t1,t2,eris)
73-
74-
# Move energy terms to the other side
75-
Foo[np.diag_indices(nocc)] -= mo_e_o
76-
Fvv[np.diag_indices(nvir)] -= mo_e_v
77-
78-
Loo = eris.Loo
79-
Lov = eris.Lov
80-
Lvv = eris.Lvv
81-
82-
# T2 equation
83-
tmp2 = lib.einsum('Lki,Lbc,ka->abic', Loo, Lvv, -t1)
84-
tmp2 += lib.einsum('Lia,Lbc->acib', Lov, Lvv).conj()
85-
tmp = lib.einsum('abic,jc->ijab', tmp2, t1)
86-
t2new = tmp + tmp.transpose(1,0,3,2)
87-
tmp2 = lib.einsum('Lkc,Lia,jc->akij', Lov, Lov, t1)
88-
tmp2 += lib.einsum('Lia,Ljk->akij', Lov, Loo).conj()
89-
tmp = lib.einsum('akij,kb->ijab', tmp2, t1)
90-
t2new -= tmp + tmp.transpose(1,0,3,2)
91-
t2new += lib.einsum('Lia,Ljb->ijab', Lov, Lov).conj()
92-
Woooo2 = lib.einsum('Lij,Lkl->ikjl', Loo, Loo)
93-
Woooo2 += lib.einsum('Llc,Lki,jc->klij', Lov, Loo, t1)
94-
Woooo2 += lib.einsum('Lkc,Llj,ic->klij', Lov, Loo, t1)
95-
Woooo2 += lib.einsum('Lkc,Lld,ic,jd->klij', Lov, Lov, t1, t1)
96-
t2new += lib.einsum('klij,ka,lb->ijab', Woooo2, t1, t1)
97-
Wvvvv = lib.einsum('Lkc,Lbd,ka->abcd', Lov, Lvv, -t1)
98-
Wvvvv = Wvvvv + Wvvvv.transpose(1,0,3,2)
99-
Wvvvv += lib.einsum('Lab,Lcd->acbd', Lvv, Lvv)
100-
t2new += lib.einsum('abcd,ic,jd->ijab', Wvvvv, t1, t1)
101-
Lvv2 = fvv - np.einsum('kc,ka->ac', fov, t1)
102-
Lvv2 -= np.diag(np.diag(fvv))
103-
tmp = lib.einsum('ac,ijcb->ijab', Lvv2, t2)
104-
t2new += (tmp + tmp.transpose(1,0,3,2))
105-
Loo2 = foo + np.einsum('kc,ic->ki', fov, t1)
106-
Loo2 -= np.diag(np.diag(foo))
107-
tmp = lib.einsum('ki,kjab->ijab', Loo2, t2)
108-
t2new -= (tmp + tmp.transpose(1,0,3,2))
109-
110-
eia = mo_e_o[:,None] - mo_e_v
111-
eijab = lib.direct_sum('ia,jb->ijab',eia,eia)
112-
t2new /= eijab
113-
114-
return t2new
113+
def make_Bov(cc, t1, eris):
114+
nocc = eris.nocc
115+
nmo = eris.fock.shape[0]
116+
nvir = nmo - nocc
117+
naux = cc._scf.with_df.get_naoaux()
118+
C = eris.mo_coeff.copy()
119+
X = C[:, nocc:] - lib.einsum('ui,ia->ua', C[:, :nocc], t1)
120+
Y = C[:, :nocc] + lib.einsum('ua,ia->ui', C[:, nocc:], t1)
121+
C[:, :nocc] = Y
122+
C[:, nocc:] = X
123+
Bov = np.empty((naux,nocc,nvir))
124+
ijslice = (0, nmo, 0, nmo)
125+
Lpq = None
126+
p1 = 0
127+
for eri1 in cc._scf.with_df.loop():
128+
Lpq = _ao2mo.nr_e2(eri1, C, ijslice, aosym='s2', out=Lpq).reshape(-1,nmo,nmo)
129+
p0, p1 = p1, p1 + Lpq.shape[0]
130+
Bov[p0:p1] = Lpq[:,:nocc,nocc:]
131+
return Bov
115132

116133
def make_t2(cc, t1, eris):
117134
nocc = eris.nocc
@@ -139,15 +156,9 @@ def make_t2(cc, t1, eris):
139156
t2 /= eijab
140157
return t2
141158

142-
def update_amps(cc, t1, t2, eris):
143-
t2 = make_t2(cc, t1, eris)
144-
t1new = update_t1(cc, t1, t2, eris)
145-
t2new = make_t2(cc, t1new, eris)
146-
# Compare t2_t1 with t2
147-
print("T2 COMPARE")
148-
#print(np.allclose(t2_t1, t2))
149-
#t2new = update_t2(cc, t1, t2, eris)
150-
return t1new, t2new
159+
def update_amps(cc, t1, eris):
160+
t1new = update_t1(cc, t1, eris)
161+
return t1new
151162

152163
# t1: ia
153164
# t2: ijab
@@ -156,14 +167,17 @@ def kernel(mycc, eris=None, t1=None, t2=None, max_cycle=50, tol=1e-8,
156167
log = logger.new_logger(mycc, verbose)
157168
if eris is None:
158169
eris = mycc.ao2mo(mycc.mo_coeff)
159-
if t1 is None and t2 is None:
160-
t1, t2 = mycc.get_init_guess(eris)
161-
elif t2 is None:
162-
t2 = mycc.get_init_guess(eris)[1]
170+
if t1 is None:
171+
mo_e = eris.mo_energy
172+
nocc = mycc.nocc
173+
eia = mo_e[:nocc,None] - mo_e[None,nocc:]
174+
t1 = eris.fock[:nocc,nocc:] / eia
175+
if t2 is None:
176+
t2 = make_t2(mycc, t1, eris)
163177

164178
cput1 = cput0 = (logger.process_clock(), logger.perf_counter())
165179
eold = 0
166-
eccsd = mycc.energy(t1, t2, eris)
180+
eccsd = mycc.energy(t1, eris)
167181
log.info('Init E_corr(CCSD) = %.15g', eccsd)
168182

169183
if isinstance(mycc.diis, lib.diis.DIIS):
@@ -176,7 +190,8 @@ def kernel(mycc, eris=None, t1=None, t2=None, max_cycle=50, tol=1e-8,
176190

177191
conv = False
178192
for istep in range(max_cycle):
179-
t1new, t2new = mycc.update_amps(t1, t2, eris)
193+
t1new = mycc.update_amps(t1, eris)
194+
t2new = make_t2(mycc, t1new, eris)
180195
tmpvec = mycc.amplitudes_to_vector(t1new, t2new)
181196
tmpvec -= mycc.amplitudes_to_vector(t1, t2)
182197
normt = np.linalg.norm(tmpvec)
@@ -189,7 +204,7 @@ def kernel(mycc, eris=None, t1=None, t2=None, max_cycle=50, tol=1e-8,
189204
t1, t2 = t1new, t2new
190205
t1new = t2new = None
191206
t1, t2 = mycc.run_diis(t1, t2, istep, normt, eccsd-eold, adiis)
192-
eold, eccsd = eccsd, mycc.energy(t1, t2, eris)
207+
eold, eccsd = eccsd, mycc.energy(t1, eris)
193208
log.info('cycle = %d E_corr(CCSD) = %.15g dE = %.9g norm(t1,t2) = %.6g',
194209
istep+1, eccsd, eccsd - eold, normt)
195210
cput1 = log.timer('CCSD iter', *cput1)
@@ -198,15 +213,84 @@ def kernel(mycc, eris=None, t1=None, t2=None, max_cycle=50, tol=1e-8,
198213
break
199214
log.timer('CCSD', *cput0)
200215
return conv, eccsd, t1, t2
216+
'''
217+
def kernel(mycc, eris=None, t1=None, max_cycle=50, tol=1e-8,
218+
tolnormt=1e-6, verbose=None):
219+
log = logger.new_logger(mycc, verbose)
220+
if eris is None:
221+
eris = mycc.ao2mo(mycc.mo_coeff)
222+
if t1 is None:
223+
mo_e = eris.mo_energy
224+
nocc = mycc.nocc
225+
eia = mo_e[:nocc,None] - mo_e[None,nocc:]
226+
t1 = eris.fock[:nocc,nocc:] / eia
227+
228+
cput1 = cput0 = (logger.process_clock(), logger.perf_counter())
229+
eold = 0
230+
# t2 = make_t2(mycc, t1, eris)
231+
eccsd = mycc.energy(t1, eris)
232+
log.info('Init E_corr(CCSD) = %.15g', eccsd)
233+
234+
if isinstance(mycc.diis, lib.diis.DIIS):
235+
adiis = mycc.diis
236+
elif mycc.diis:
237+
adiis = lib.diis.DIIS(mycc, mycc.diis_file, incore=mycc.incore_complete)
238+
adiis.space = mycc.diis_space
239+
else:
240+
adiis = None
201241
242+
conv = False
243+
for istep in range(max_cycle):
244+
t1new = mycc.update_amps(t1, eris)
245+
tmpvec = mycc.amplitudes_to_vector(t1new)
246+
tmpvec -= mycc.amplitudes_to_vector(t1)
247+
normt = np.linalg.norm(tmpvec)
248+
normt += np.linalg.norm(make_t2(mycc, t1new, eris).ravel() - make_t2(mycc, t1, eris).ravel())
249+
tmpvec = None
250+
if mycc.iterative_damping < 1.0:
251+
alpha = mycc.iterative_damping
252+
t1new = (1-alpha) * t1 + alpha * t1new
253+
t1 = t1new
254+
t1new = None
255+
t1 = mycc.run_diis(t1, istep, normt, eccsd-eold, adiis)
256+
eold, eccsd = eccsd, mycc.energy(t1, eris)
257+
log.info('cycle = %d E_corr(CCSD) = %.15g dE = %.9g norm(t1,t2) = %.6g',
258+
istep+1, eccsd, eccsd - eold, normt)
259+
cput1 = log.timer('CCSD iter', *cput1)
260+
if abs(eccsd-eold) < tol and normt < tolnormt:
261+
conv = True
262+
break
263+
log.timer('CCSD', *cput0)
264+
return conv, eccsd, t1
265+
'''
202266
class DFRCC2(ccsd.CCSD):
203267
'''restricted CCSD with IP-EOM, EA-EOM, EE-EOM, and SF-EOM capabilities
204268
205269
Ground-state CCSD is performed in optimized ccsd.CCSD and EOM is performed here.
206270
'''
271+
energy = energy
207272
kernel = kernel
208273
update_amps = update_amps
274+
'''
275+
def amplitudes_to_vector(self, t1, out=None):
276+
vector = np.ndarray(t1.size, t1.dtype, buffer=out)
277+
vector = t1.ravel()
278+
return vector
279+
280+
def vector_to_amplitudes(self, vector, nmo=None, nocc=None):
281+
if nocc is None: nocc = self.nocc
282+
if nmo is None: nmo = self.nmo
283+
return vector.reshape(nocc,nmo-nocc)
209284
285+
def run_diis(self, t1, istep, normt, de, adiis):
286+
if (adiis and
287+
istep >= self.diis_start_cycle and
288+
abs(de) < self.diis_start_energy_diff):
289+
vec = self.amplitudes_to_vector(t1)
290+
t1 = self.vector_to_amplitudes(adiis.update(vec))
291+
logger.debug1(self, 'DIIS for step %d', istep)
292+
return t1
293+
'''
210294
def ao2mo(self, mo_coeff=None):
211295
cput0 = (logger.process_clock(), logger.perf_counter())
212296
log = logger.Logger(self.stdout, self.verbose)

0 commit comments

Comments
 (0)