Commit 9e942504 authored by Davison, Jacob's avatar Davison, Jacob
Browse files

fixed spin-squared stuff; normal-ordering SS correctly

parent 9e2b03e1
......@@ -98,13 +98,13 @@ def derivative(t, y, inputs):
E, f, G = ravel(y[0:half], hamiltonian.n_sp_states)
generator.f = f
generator.G = G
dE, df, dG = flow.flow(f, G, generator)
dE, df, dG = flow.flow(E, f, G, generator)
# Spin-squared flow
E_spin, f_spin, G_spin = ravel(y[half::], tspinsq.n_sp_states)
# generator_spin.f = f_spin
# generator_spin.G = G_spin
dE_spin, df_spin, dG_spin = flow_spin.flow(f_spin, G_spin, generator)
dE_spin, df_spin, dG_spin = flow_spin.flow(E_spin, f_spin, G_spin, generator)
dy = np.concatenate([unravel(dE, df, dG), unravel(dE_spin, df_spin, dG_spin)], axis=0)
......
......@@ -742,8 +742,8 @@ class WhiteGenerator(Generator):
f = self.f
G = self.G
eta1B, eta2B = self._wrapper_calc_eta_generalized_ph(bas1B, holes, particles, f, G, ref)
#eta1B, eta2B = self._wrapper_calc_eta(bas1B, holes, particles, f, G, ref)
#eta1B, eta2B = self._wrapper_calc_eta_generalized_ph(bas1B, holes, particles, f, G, ref)
eta1B, eta2B = self._wrapper_calc_eta(bas1B, holes, particles, f, G, ref)
#print(np.array_equal(eta1B1, eta1B), np.array_equal(eta2B1, eta2B))
self._eta1B = eta1B
......
......@@ -7,7 +7,7 @@ tn.set_default_backend("numpy")
class TSpinSq(object):
"""Generate the vacuum S^2 operator."""
def __init__(self, n_hole_states, n_particle_states, ref=[], s=0.5):
def __init__(self, n_hole_states, n_particle_states, rho1b, rho2b, ref=[], s=0.5):
"""Class constructor. Instantiate PairingHamiltonian2B object.
Arguments:
......@@ -35,7 +35,7 @@ class TSpinSq(object):
self._sp_basis = np.append(self.holes, self.particles)
self._SS1B, self._SS2B = self.construct()
self._E, self._f, self._G = self.normal_order()
self._E, self._f, self._G = self.normal_order_slow(rho1b, rho2b)
@property
def s(self):
......@@ -267,4 +267,36 @@ class TSpinSq(object):
G = H2B_t
return (E, f, G)
def normal_order_slow(self, rho1b, rho2b):
bas1B = self.sp_basis # get the single particle basis
H1B_t = self.SS1B.astype(np.float32) # get the 1B tensor
H2B_t = self.SS2B.astype(np.float32) # get the 2B tensor
holes = self.holes # get hole states
particles = self.particles # get particle states
n_states = len(bas1B)
# hme = pyci.matrix(len(holes), len(particles), 0.0, H1B_t, H2B_t, H2B_t, imsrg=True)
# w,v = np.linalg.eigh(hme)
# v0 = v[:, 0]
contract_1b = np.einsum('ij,ij->', rho1b, H1B_t)
# rho_reshape_2b = np.reshape(rho2b, (n_states**2,n_states**2))
# h2b_reshape_2b = np.reshape(H2B_t, (n_states**2,n_states**2))
# contract_2b = np.einsum('ij,ij->', rho_reshape_2b, h2b_reshape_2b)
contract_2b = np.einsum('ijkl,ijkl->', rho2b, H2B_t)
E = contract_1b + 0.25*contract_2b
f = H1B_t + np.einsum('piqj,ij->pq', H2B_t, rho1b) #0.5*(np.einsum('piqj,ij->pq', H2B_t, rho1b) - np.einsum('pijq,ij->pq', H2B_t, rho1b))
G = H2B_t
#print("NOH contract_1b, contract_2b, piqj", contract_1b, contract_2b, np.einsum('piqj,ij->pq', H2B_t, rho1b))
return (E, f, G)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment