From efc3c1c1d093ed8172453c4d841dba1e94d38053 Mon Sep 17 00:00:00 2001 From: Matthew R Hermes Date: Mon, 9 Dec 2024 09:48:40 -0600 Subject: [PATCH 1/2] debug min_2frag_tps.py --- my_pyscf/lassi/min_2frag_tps.py | 89 +++++++++++++++++++++------------ 1 file changed, 56 insertions(+), 33 deletions(-) diff --git a/my_pyscf/lassi/min_2frag_tps.py b/my_pyscf/lassi/min_2frag_tps.py index c33b01ba..1462ec4e 100644 --- a/my_pyscf/lassi/min_2frag_tps.py +++ b/my_pyscf/lassi/min_2frag_tps.py @@ -3,17 +3,20 @@ from scipy import linalg from mrh.my_pyscf.lassi.citools import get_lroots -def hess_ss (ham_pq, si): +def get_hess_ss (ham_pq, si, lroots, nroots, proj): + nstates = ham_pq.shape[1] + p = np.prod (lroots) hs = np.dot (ham_pq, si) e = np.dot (si.conj (), hs) - hess = ham_pq - (np.multiply.outer (si.conj (), si) * e) + hess = ham_pq - e*np.eye (ham_pq.shape[0]) hs -= e*si hess -= np.multiply.outer (si.conj (), hs) hess -= np.multiply.outer (hs.conj (), si) hess += hess.conj ().T + hess = proj.conj ().T @ hess @ proj return hess -def hess_us (ham_pq, si, lroots, nroots): +def get_hess_us (ham_pq, si, lroots, nroots, proj): nstates = ham_pq.shape[1] p = np.prod (lroots) hs = np.dot (ham_pq, si) @@ -31,25 +34,26 @@ def hess_us (ham_pq, si, lroots, nroots): idx_u = np.repeat (list (range (lroots[1])), lroots[0]) idx_v = np.tile (list (range (lroots[0])), lroots[1]) for i in range (nroots): - delta_u = np.where (idx_u==i)[0] - hess_us[i,:,delta_u] += hs_p[:,idx_v[delta_u]] - hess_us[:,i,delta_u] -= hs_p[:,idx_v[delta_u]] + delta_u = (idx_u==i) + hess_us[i,:,:p][:,delta_u] += hs_p[:,idx_v[delta_u]] + hess_us[:,i,:p][:,delta_u] -= hs_p[:,idx_v[delta_u]] hess_us = hess_us[nroots:,:nroots,:].reshape (-1,nstates) for i in range (nroots): - delta_v = np.where (idx_v==i)[0] - hess_vs[i,:,delta_v] += hs_p[idx_u[delta_v],:].T - hess_vs[:,i,delta_v] -= hs_p[idx_u[delta_v],:].T + delta_v = (idx_v==i) + hess_vs[i,:,:p][:,delta_v] += hs_p[idx_u[delta_v],:].T + hess_vs[:,i,:p][:,delta_v] -= hs_p[idx_u[delta_v],:].T hess_vs = hess_vs[nroots:,:nroots,:].reshape (-1,nstates) hess = np.append (hess_us, hess_vs, axis=0) + hess = np.dot (hess, proj) return hess + hess.conj () -def hess_uu (ham_pq, si, lroots, nroots): +def get_hess_uu (ham_pq, si, lroots, nroots): p = np.prod (lroots) hs = np.dot (ham_pq, si)[:p].reshape (lroots[1],lroots[0]) h_pp = ham_pq[:p,:p].reshape (lroots[1],lroots[0],lroots[1],lroots[0]) - si_p = si[:p].reshape (lroots[1],lroots[0]) + si = si[:p].reshape (lroots[1],lroots[0]) hess = lib.einsum ('im,jmln,kn->ijkl', si.conj (), h_pp, si) - fu = np.dot (si_p.conj (), hs.T) + fu = np.dot (si.conj (), hs.T) for i in range (lroots[1]): hess[:,i,i,:] += fu hess -= hess.transpose (1,0,2,3) @@ -58,14 +62,14 @@ def hess_uu (ham_pq, si, lroots, nroots): hess = hess[nroots:,:nroots,nroots:,:nroots].reshape (nel, nel) return hess + hess.conj () -def hess_vv (ham_pq, si, lroots, nroots): +def get_hess_vv (ham_pq, si, lroots, nroots): p = np.prod (lroots) hs = np.dot (ham_pq, si)[:p].reshape (lroots[1],lroots[0]) h_pp = ham_pq[:p,:p].reshape (lroots[1],lroots[0],lroots[1],lroots[0]) - si_p = si[:p].reshape (lroots[1],lroots[0]) + si = si[:p].reshape (lroots[1],lroots[0]) hess = lib.einsum ('mi,mjnl,nk->ijkl', si.conj (), h_pp, si) - fv = np.dot (si_p.conj ().T, hs) - for i in range (lroots[1]): + fv = np.dot (si.conj ().T, hs) + for i in range (lroots[0]): hess[:,i,i,:] += fv hess -= hess.transpose (1,0,2,3) hess -= hess.transpose (0,1,3,2) @@ -73,11 +77,11 @@ def hess_vv (ham_pq, si, lroots, nroots): hess = hess[nroots:,:nroots,nroots:,:nroots].reshape (nel, nel) return hess + hess.conj () -def hess_uv (ham_pq, si, lroots, nroots): +def get_hess_uv (ham_pq, si, lroots, nroots): p = np.prod (lroots) hs = np.dot (ham_pq, si)[:p].reshape (lroots[1],lroots[0]) h_pp = ham_pq[:p,:p].reshape (lroots[1],lroots[0],lroots[1],lroots[0]) - si_p = si[:p].reshape (lroots[1],lroots[0]) + si = si[:p].reshape (lroots[1],lroots[0]) hess = lib.einsum ('im,jmnl,nk->ijkl', si.conj (), h_pp, si) hess += np.multiply.outer (si.conj (), hs).transpose (0,2,1,3) hess -= hess.transpose (1,0,2,3) @@ -87,25 +91,26 @@ def hess_uv (ham_pq, si, lroots, nroots): hess = hess[nroots:,:nroots,nroots:,:nroots].reshape (nelu, nelv) return hess + hess.conj () -def hess (ham_pq, si, lroots, nroots): - hess_uu = hess_uu (ham_pq, si, lroots, nroots) - hess_uv = hess_uv (ham_pq, si, lroots, nroots) - hess_vv = hess_vv (ham_pq, si, lroots, nroots) - hess_uu = np.append ([np.append (hess_uu, hess_uv, axis=1), - np.append (hess_uv.T, hess_vv, axis=1)], - axis=0) - hess_us = hess_us (ham_pq, si, lroots, nroots) - hess_ss = hess_ss (ham_pq, si) - hess = np.append ([np.append (hess_uu, hess_us, axis=1), - np.append (hess_us.T, hess_ss, axis=1)], +def get_hess (ham_pq, si, lroots, nroots, proj): + huu = get_hess_uu (ham_pq, si, lroots, nroots) + huv = get_hess_uv (ham_pq, si, lroots, nroots) + hvv = get_hess_vv (ham_pq, si, lroots, nroots) + huu = np.append (np.append (huu, huv, axis=1), + np.append (huv.T, hvv, axis=1), + axis=0) + hus = get_hess_us (ham_pq, si, lroots, nroots, proj) + hss = get_hess_ss (ham_pq, si, lroots, nroots, proj) + hess = np.append (np.append (huu, hus, axis=1), + np.append (hus.T, hss, axis=1), axis=0) return hess -def grad (ham_pq, si, lroots, nroots): +def get_grad (ham_pq, si, lroots, nroots, proj): p = np.prod (lroots) hs = np.dot (ham_pq, si) e = np.dot (si.conj (), hs) grad_s = hs - e*si + grad_s = np.dot (grad_s, proj) hs = hs[:p].reshape (lroots[1], lroots[0]) si = si[:p].reshape (lroots[1], lroots[0]) fu = np.dot (si.conj (), hs.T) @@ -116,14 +121,32 @@ def grad (ham_pq, si, lroots, nroots): grad_v = grad_v[nroots:,:nroots].ravel () return np.concatenate ([grad_u, grad_v, grad_s]) +def get_proj (si, lroots, nroots): + nstates = len (si) + p = np.prod (lroots) + idx_p = np.zeros (lroots, dtype=bool) + idx_p[:nroots,:nroots] = True + idx_q = np.ones (nstates-p, dtype=bool) + idx = np.append (idx_p.ravel (), idx_q) + x = np.eye (nstates) + x = x[:,idx] + Q, R = linalg.qr (x.conj ().T @ si[:,None]) + proj = x @ Q[:,1:] + return proj + def quadratic_step (ham_pq, ci1, si_p, si_q): nroots = len (si_p) lroots = get_lroots (ci1) si = np.zeros (lroots[::-1]) si[:nroots,:nroots] = np.diag (si_p) si = np.append (si.ravel (), si_q) - x = linalg.solve (hess (ham_pq, si, lroots, nroots), - -grad (ham_pq, si, lroots, nroots)) + ham_pq = ham_pq.copy () + e = np.dot (si.conj (), np.dot (ham_pq, si)) + ham_pq -= np.eye (ham_pq.shape[0])*e + proj = get_proj (si, lroots, nroots) + hess = get_hess (ham_pq, si, lroots, nroots, proj) + grad = get_grad (ham_pq, si, lroots, nroots, proj) + x = linalg.solve (hess, -grad) nelu = (lroots[1]-nroots)*nroots nelv = (lroots[0]-nroots)*nroots kappa_u = np.zeros ((lroots[1],lroots[1]), dtype=x.dtype) @@ -132,7 +155,7 @@ def quadratic_step (ham_pq, ci1, si_p, si_q): kappa_u -= kappa_u.T u = linalg.expm (kappa_u)[:,:nroots] kappa_v = np.zeros ((lroots[0],lroots[0]), dtype=x.dtype) - kappa_u[nroots:,:nroots] = x[:nelv].reshape ( + kappa_v[nroots:,:nroots] = x[:nelv].reshape ( lroots[0]-nroots, nroots) kappa_v -= kappa_v.T vh = linalg.expm (kappa_v)[:,:nroots].conj ().T From 18b093226e5ec3882345a824416066f74f677a65 Mon Sep 17 00:00:00 2001 From: Matthew R Hermes Date: Mon, 9 Dec 2024 10:56:56 -0600 Subject: [PATCH 2/2] test min_2frag_tps --- my_pyscf/lassi/excitations.py | 4 ++-- my_pyscf/lassi/min_2frag_tps.py | 25 +++++++++++++++++++++++-- 2 files changed, 25 insertions(+), 4 deletions(-) diff --git a/my_pyscf/lassi/excitations.py b/my_pyscf/lassi/excitations.py index eb50aea8..79e31ff1 100644 --- a/my_pyscf/lassi/excitations.py +++ b/my_pyscf/lassi/excitations.py @@ -626,11 +626,11 @@ def truncrot_ham_pq (self, ham_pq, u, vh): p = np.prod (lroots) nstates = ham_pq.shape[1] h_pr = ham_pq[:p,:].reshape (lroots[1],lroots[0],nstates) - h_pr = np.dot (uh, np.dot (vh, h_pr)).reshape (nroots*nroots,nstates) + h_pr = np.dot (uh, np.dot (vh, h_pr)).reshape (u.shape[1]*v.shape[1],nstates) ham_pq = np.append (h_pr, ham_pq[p:,:], axis=0) nstates = ham_pq.shape[0] h_rp = ham_pq[:,:p].reshape (nstates,lroots[1],lroots[0]) - h_rp = lib.einsum ('rij,ia,jb->rab', h_rp, u, v).reshape (nstates,nroots*nroots) + h_rp = lib.einsum ('rij,ia,jb->rab', h_rp, u, v).reshape (nstates,u.shape[1]*v.shape[1]) ham_pq = np.append (h_rp, ham_pq[:,p:], axis=1) return ham_pq diff --git a/my_pyscf/lassi/min_2frag_tps.py b/my_pyscf/lassi/min_2frag_tps.py index 1462ec4e..71ad71e1 100644 --- a/my_pyscf/lassi/min_2frag_tps.py +++ b/my_pyscf/lassi/min_2frag_tps.py @@ -153,14 +153,35 @@ def quadratic_step (ham_pq, ci1, si_p, si_q): kappa_u[nroots:,:nroots] = x[:nelu].reshape ( lroots[1]-nroots, nroots) kappa_u -= kappa_u.T - u = linalg.expm (kappa_u)[:,:nroots] + u = linalg.expm (-kappa_u/2)[:,:nroots] kappa_v = np.zeros ((lroots[0],lroots[0]), dtype=x.dtype) kappa_v[nroots:,:nroots] = x[:nelv].reshape ( lroots[0]-nroots, nroots) kappa_v -= kappa_v.T - vh = linalg.expm (kappa_v)[:,:nroots].conj ().T + vh = linalg.expm (-kappa_v/2)[:,:nroots].conj ().T return u, vh +def idx_down_ham_pq (ham_pq, lroots, nroots): + p = np.prod (lroots) + q = ham_pq.shape[0] - p + h_pr = ham_pq[:p,:].reshape (lroots[1],lroots[0],p+q) + h_pr = h_pr[:nroots,:nroots,:].reshape (nroots*nroots,p+q) + ham_pq = np.append (h_pr, ham_pq[p:,:], axis=0) + h_rp = ham_pq[:,:p].reshape (nroots*nroots+q,lroots[1],lroots[0]) + h_rp = h_rp[:,:nroots,:nroots].reshape (nroots*nroots+q,nroots*nroots) + ham_pq = np.append (h_rp, ham_pq[:,p:], axis=1) + return ham_pq + +def subspace_eig (ham_pq, lroots, nroots): + p = np.prod (lroots) + q = ham_pq.shape[0] - p + p0 = nroots*nroots + evals, evecs = linalg.eigh (idx_down_ham_pq (ham_pq, lroots, nroots)) + si_p = np.zeros ((lroots[1],lroots[0],p0+q), dtype=evecs.dtype) + si_p[:nroots,:nroots,:] = evecs[:p0,:].reshape (nroots,nroots,p0+q) + si = np.append (si_p.reshape (p,p0+q), evecs[p0:,:], axis=0) + return evals, si +