diff --git a/spin_systems/ClassicalIsingTermalState.py b/spin_systems/ClassicalIsingTermalState.py
index d125e29..992e3bf 100644
--- a/spin_systems/ClassicalIsingTermalState.py
+++ b/spin_systems/ClassicalIsingTermalState.py
@@ -1,5 +1,7 @@
 import numpy as np
+import opt_einsum as oe
 from ..tensors.MatrixProductState import MPS
+from ..tensors.iMatrixProductState import iMPS
 
 def Ising_ThermalState(L, K):
     """
@@ -36,4 +38,13 @@ def Ising_ThermalState(L, K):
         tensors[i] = tensors[i].reshape(shp[0],shp[1]*shp[2],shp[3])
     rho_t = MPS(L,d=4,tensors=tensors)
     rho_t.right_normalize()
-    return rho_t
\ No newline at end of file
+    return rho_t
+
+def iMPSstate(K):
+    B = Ising_ThermalState(3, K).tensors[1]
+    theta = oe.contract('abc,cde->abde',B,B).reshape(8,8)
+    Sb = np.linalg.svd(theta,compute_uv=False)[:2]
+    Sb /= np.linalg.norm(Sb)
+    return iMPS(Sb,B,B,d=4)
+    
+    
\ No newline at end of file
diff --git a/spin_systems/Hamiltonians.py b/spin_systems/Hamiltonians.py
index 8695939..341567b 100644
--- a/spin_systems/Hamiltonians.py
+++ b/spin_systems/Hamiltonians.py
@@ -91,5 +91,4 @@ def Hierarchical_Dyson_model(N,sigma,h,J=1):
             tensors[k][i+1,i] = np.eye(2)
             if k+i<L:
                 tensors[k][-1,i,:,:] = A[k,k+i]*sigma_x
-    return MPO(L,tensors=tensors)
-    
\ No newline at end of file
+    return MPO(L,tensors=tensors)
\ No newline at end of file
diff --git a/spin_systems/__init__.py b/spin_systems/__init__.py
index e9286c9..01e7239 100644
--- a/spin_systems/__init__.py
+++ b/spin_systems/__init__.py
@@ -1,4 +1,5 @@
 from . import Hamiltonians
 from . import observables
 from . import states
-from . import ClassicalIsingTermalState
\ No newline at end of file
+from . import ClassicalIsingTermalState
+from . import iHamiltonians
\ No newline at end of file
diff --git a/spin_systems/__pycache__/ClassicalIsingTermalState.cpython-38.pyc b/spin_systems/__pycache__/ClassicalIsingTermalState.cpython-38.pyc
index d319931..d5104a0 100644
Binary files a/spin_systems/__pycache__/ClassicalIsingTermalState.cpython-38.pyc and b/spin_systems/__pycache__/ClassicalIsingTermalState.cpython-38.pyc differ
diff --git a/spin_systems/iHamiltonians.py b/spin_systems/iHamiltonians.py
new file mode 100644
index 0000000..180bfad
--- /dev/null
+++ b/spin_systems/iHamiltonians.py
@@ -0,0 +1,18 @@
+import numpy as np
+
+def IsingChainWithAncilla(J=1., h_z=0.25, h_x=0.):
+    sz = np.array([[1,0],[0,-1]])
+    sx = np.array([[0,1],[1,0]])
+    Id2 = np.eye(2)
+    Id4 = np.eye(4)
+    
+    sz_sys = np.kron(sz,Id2)
+    sx_sys = np.kron(sx,Id2)
+    
+    sz_anc = np.kron(Id2,sz)
+    sx_anc = np.kron(Id2,sz)
+
+    H  = -J*(np.kron(sx_sys,sx_sys))-h_z*(np.kron(sz_sys,Id4)+np.kron(Id4,sz_sys))-h_x*(np.kron(sx_sys,Id4)+np.kron(Id4,sx_sys))
+    H  -= -J*(np.kron(sx_anc,sx_anc))-h_z*(np.kron(sz_anc,Id4)+np.kron(Id4,sz_anc))-h_x*(np.kron(sx_anc,Id4)+np.kron(Id4,sx_anc))
+    
+    return H
\ No newline at end of file