Skip to content

Commit

Permalink
Getting close!
Browse files Browse the repository at this point in the history
  • Loading branch information
@lothian
lothian committed Jan 2, 2024
1 parent e5540a5 commit 74219d3
Show file tree
Hide file tree
Showing 2 changed files with 77 additions and 31 deletions.
108 changes: 77 additions & 31 deletions pycc/cceom.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,32 +42,77 @@ def solve_eom(self, nstates=1, e_conv=1e-5, r_conv=1e-5, maxiter=100):
H = hbar.ccwfn.H
contract = hbar.contract

L = nstates * 3 # size of guess space (varies)
maxL = L * 10 # max size of subspace (fixed)
M = nstates * 2 # size of guess space (varies)
maxM = M * 10 # max size of subspace (fixed)

# Initialize guess vectors
_, C1 = self.cis(L, o, v, H, contract)
C2 = np.zeros((L, no, no, nv, nv))

# Compute sigma vectors
s1 = np.zeros_like(C1)
s2 = np.zeros_like(C2)
for state in range(L):
s1[state] = self.s1(hbar, C1[state], C2[state])
s2[state] = self.s2(hbar, C1[state], C2[state])

# Build and diagonalize subspace Hamiltonian
G = np.zeros((L,L))

# Build and orthonormalize correction vectors

#
E, C1 = self.guess(M, o, v, H, contract)
C2 = np.zeros((M, no, no, nv, nv))
E = E[:nstates]
print("EOM Iter %3d: M = %3d" % (0, M))
print(" E/state dE norm")
for state in range(nstates):
print("%20.12f --- ---" % (E[state]))

# Build preconditioner (energy denominator)
hbar_occ = np.diag(hbar.Hoo)
hbar_vir = np.diag(hbar.Hvv)
Dia = hbar_occ.reshape(-1,1) - hbar_vir
Dijab = (hbar_occ.reshape(-1,1,1,1) + hbar_occ.reshape(-1,1,1) -
hbar_vir.reshape(-1,1) - hbar_vir)
D = np.hstack((Dia.flatten(), Dijab.flatten()))

maxiter = 20
for niter in range(maxiter):
E_old = E

# Compute sigma vectors
s1 = np.zeros_like(C1)
s2 = np.zeros_like(C2)
for state in range(M):
s1[state] = self.s1(hbar, C1[state], C2[state])
s2[state] = self.s2(hbar, C1[state], C2[state])

# Build and diagonalize subspace Hamiltonian
S = np.hstack((np.reshape(s1, (M, no*nv)), np.reshape(s2, (M, no*no*nv*nv))))
C = np.hstack((np.reshape(C1, (M, no*nv)), np.reshape(C2, (M, no*no*nv*nv))))
G = C @ S.T
l, a = np.linalg.eigh(G)
E = l[:nstates]

# Build correction vectors
# r --> (M, no*nv + no*no*nv*nv)
# a --> (M, M)
# S --> (M, no*nv + no*no*nv*nv)
# C --> (M, no*nv + no*no*nv*nv)
# l --> (M)
r = a @ (S - np.diag(l) @ C)
delta = r/np.subtract.outer(l,D) # element-by-element division

# Add new vectors to guess space and orthonormalize
Q, _ = np.linalg.qr(np.concatenate((C, delta[:nstates])).T)
C = Q.T.copy()
M = C.shape[0]

# Print status and check convergence and print status
dE = E - E_old
r_norm = np.linalg.norm(r[:nstates], axis=1)
print("EOM Iter %3d: M = %3d" % (niter, M))
print(" E/state dE norm")
for state in range(nstates):
print("%20.12f %20.12f %20.12f" % (E[state], dE[state], r_norm[state]))



# Re-shape guess vectors for next iteration
C1 = np.reshape(C[:,:no*nv], (M,no,nv)).copy()
C2 = np.reshape(C[:,no*nv:], (M,no,no,nv,nv)).copy()

print("\nCCEOM converged in %.3f seconds." % (time.time() - time_init))



def cis(self, nstates, o, v, H, contract):
def guess(self, M, o, v, H, contract, method='CCS'):
"""
Compute CIS excited states as guess vectors
"""
Expand All @@ -76,15 +121,16 @@ def cis(self, nstates, o, v, H, contract):
F = H.F
L = H.L

# Build CIS matrix and diagonalize
CIS_H = L[v,o,o,v].swapaxes(0,1).swapaxes(0,2).copy()
CIS_H += contract('ab,ij->iajb', F[v,v], np.eye(no))
CIS_H -= contract('ij,ab->iajb', F[o,o], np.eye(nv))

eps, c = np.linalg.eigh(np.reshape(CIS_H, (no*nv,no*nv)))

c = np.reshape(c.T[slice(0,nstates),:], (nstates, no, nv)).copy()
# Build list of guess vectors (C1 tensors)
guesses = np.reshape(c.T[slice(0,M),:], (M, no, nv)).copy()

return eps, c
return eps[:M], guesses


def s1(self, hbar, C1, C2):
Expand All @@ -93,13 +139,13 @@ def s1(self, hbar, C1, C2):

s1 = contract('ie,ae->ia', C1, hbar.Hvv)
s1 -= contract('mi,ma->ia', hbar.Hoo, C1)
s1 += 2.0 * contract('maei,me->ia', hbar.Hovvo, C1)
s1 += contract('maei,me->ia', hbar.Hovvo, C1) * 2.0
s1 -= contract('maie,me->ia', hbar.Hovov, C1)
s1 += 2.0 * contract('miea,me->ia', C2, hbar.Hov)
s1 += contract('miea,me->ia', C2, hbar.Hov) * 2.0
s1 -= contract('imea,me->ia', C2, hbar.Hov)
s1 += 2.0 * contract('imef,amef->ia', C2, hbar.Hvovv)
s1 += contract('imef,amef->ia', C2, hbar.Hvovv) * 2.0
s1 -= contract('imef,amfe->ia', C2, hbar.Hvovv)
s1 -= 2.0 * contract('mnie,mnae->ia', hbar.Hooov, C2)
s1 -= contract('mnie,mnae->ia', hbar.Hooov, C2) * 2.0
s1 += contract('nmie,mnae->ia', hbar.Hooov, C2)

return s1
Expand All @@ -113,11 +159,11 @@ def s2(self, hbar, C1, C2):
o = hbar.ccwfn.o
v = hbar.ccwfn.v

Zvv = 2.0 * contract('amef,mf->ae', hbar.Hvovv, C1)
Zvv = contract('amef,mf->ae', hbar.Hvovv, C1) * 2.0
Zvv -= contract('amfe,mf->ae', hbar.Hvovv, C1)
Zvv -= contract('nmaf,nmef->ae', C2, L[o,o,v,v])

Zoo = -2.0 * contract('mnie,ne->mi', hbar.Hooov, C1)
Zoo = contract('mnie,ne->mi', hbar.Hooov, C1) * -2.0
Zoo += contract('nmie,ne->mi', hbar.Hooov, C1)
Zoo -= contract('mnef,inef->mi', L[o,o,v,v], C2)

Expand All @@ -127,11 +173,11 @@ def s2(self, hbar, C1, C2):
s2 += contract('mi,mjab->ijab', Zoo, t2)
s2 += contract('ijeb,ae->ijab', C2, hbar.Hvv)
s2 -= contract('mi,mjab->ijab', hbar.Hoo, C2)
s2 += 0.5 * contract('mnij,mnab->ijab', hbar.Hoooo, C2)
s2 += 0.5 * contract('ijef,abef->ijab', C2, hbar.Hvvvv)
s2 += contract('mnij,mnab->ijab', hbar.Hoooo, C2) * 0.5
s2 += contract('ijef,abef->ijab', C2, hbar.Hvvvv) * 0.5
s2 -= contract('imeb,maje->ijab', C2, hbar.Hovov)
s2 -= contract('imea,mbej->ijab', C2, hbar.Hovvo)
s2 += 2.0 * contract('miea,mbej->ijab', C2, hbar.Hovvo)
s2 += contract('miea,mbej->ijab', C2, hbar.Hovvo) * 2.0
s2 -= contract('miea,mbje->ijab', C2, hbar.Hovov)

return s2 + s2.swapaxes(0,1).swapaxes(2,3)
0 pycc/tests/test_034_eomccsd.py → pycc/tests/test_035_eomccsd.py
View file
File renamed without changes.

0 comments on commit 74219d3

Please sign in to comment.