diff --git a/layers/TensorDense.py b/layers/TensorDense.py
index c11ce24..e138fa3 100644
--- a/layers/TensorDense.py
+++ b/layers/TensorDense.py
@@ -42,10 +42,6 @@ def build(self, tt_input_shape):
 
             self.tt_input_shape = tuple([roots] * self.cores_number)
 
-            #print("reshaped", tt_input_shape)
-
-        #print("Input shape", self.tt_input_shape)
-
         self.bias = tf.Variable(tf.zeros(shape=self.tt_output_shape), name="bias", trainable=True)
         
         last_idx = self.cores_number-1        
diff --git a/tensor/__init__.py b/tensor/__init__.py
index 8787b7e..dfdf8fc 100644
--- a/tensor/__init__.py
+++ b/tensor/__init__.py
@@ -1,4 +1,5 @@
 from syngular.tensor.matrix_product_state import MatrixProductState
 from syngular.tensor.matrix_product_operator import MatrixProductOperator
 from syngular.tensor.matrix_product_density_operator import MatrixProductDensityOperator
-from syngular.tensor.differential_matrix_product_operator import DifferentialMatrixProductOperator
\ No newline at end of file
+from syngular.tensor.differential_matrix_product_operator import DifferentialMatrixProductOperator
+from syngular.tensor.generalized_site import GeneralizedSite, Index
\ No newline at end of file
diff --git a/tensor/generalized_product_state.py b/tensor/generalized_product_state.py
new file mode 100644
index 0000000..e69de29
diff --git a/tensor/generalized_site.py b/tensor/generalized_site.py
new file mode 100644
index 0000000..6a76d1d
--- /dev/null
+++ b/tensor/generalized_site.py
@@ -0,0 +1,21 @@
+from __future__ import annotations
+from typing import Tuple, List, Union
+
+import numpy as np
+from scipy import linalg
+
+from opt_einsum import contract
+
+class Index:
+    name: str
+    dim: int
+
+GeneralizedSite = dict[str, list[Index]]
+
+# class GeneralizedSite:
+
+#     def __init__(self, groups: Group) -> None:
+#         self.groups = groups
+
+#     def __getitem__(self, index: str):
+#         return self.groups[]
\ No newline at end of file
diff --git a/tensor/peps.py b/tensor/peps.py
new file mode 100644
index 0000000..2466902
--- /dev/null
+++ b/tensor/peps.py
@@ -0,0 +1,59 @@
+from __future__ import annotations
+from typing import Tuple, List, Union
+
+import numpy as np
+from pandas import array
+from scipy import linalg
+
+from opt_einsum import contract
+
+class PEPState:
+
+    """
+    Example: a 3x3 PEPS
+           |   |   |
+        ---O---O---O---
+         / | / | / |
+        ---O---O---O---
+         / | / | / |
+        ---O---O---O---
+         / | / | / |
+
+    Convention:
+
+               | (t)
+        (l) ---O--- (r)  ==> (l, t, i, b, r)
+        (i)  / | (b)
+
+    """
+
+    def __init__(self, tensor, bond_matrix: List[List[int]]) -> None:
+        self.parameters_number = 0
+        self.real_parameters_number = 0
+
+        if tensor is not None:
+            self.tensor_shape = tensor.shape
+            self.bond_matrix = np.array(bond_matrix)
+
+            self.width = bond_matrix.shape[0]
+            self.height = bond_matrix.shape[1]
+
+            self.order = len(self.tensor_shape)
+            self.real_parameters_number = np.prod(self.tensor_shape)
+            
+            self.sites_number = len(self.bond_shape)+1
+            self.sites = [None] * self.sites_number
+
+            if self.bond_matrix != ():
+                self.shape = []
+                for w in self.width:
+                    for h in self.height:
+                        shp = (0,0,0,0)
+                        if w == 0:
+                            shp[0] = 1
+                        elif w == self.width-1:
+                            shp[3] = 1
+
+    @staticmethod
+    def ground(width: int, height: int) -> PEPState:
+        pass
diff --git a/test/test_dmrg.py b/test/test_dmrg.py
new file mode 100644
index 0000000..2959bab
--- /dev/null
+++ b/test/test_dmrg.py
@@ -0,0 +1,13 @@
+from numpy import outer
+from syngular.variational import DMRG
+from syngular.tensor import MatrixProductOperator
+
+if __name__ == "__main__":
+
+    input_shape = (3,3,3,3,3,3)
+    output_shape = (3,3,3,3,3,3)
+    bond_shape = (2,2,2,2,2)
+
+    mpo = MatrixProductOperator.random(input_shape=input_shape, output_shape=output_shape, bond_shape=bond_shape)
+
+    DMRG.solve(mpo)
\ No newline at end of file
diff --git a/test/test_sites.py b/test/test_sites.py
new file mode 100644
index 0000000..5bb2e99
--- /dev/null
+++ b/test/test_sites.py
@@ -0,0 +1,8 @@
+from syngular.tensor import GeneralizedSite, Index
+
+site_1: GeneralizedSite = {
+    "input": {
+        'i1': 10,
+        'i2': 5
+    }
+}
\ No newline at end of file
diff --git a/variational/__init__.py b/variational/__init__.py
new file mode 100644
index 0000000..8cc4a75
--- /dev/null
+++ b/variational/__init__.py
@@ -0,0 +1,2 @@
+from syngular.variational.dmrg import DMRG
+from syngular.variational.lanczos import Lanczos
\ No newline at end of file
diff --git a/variational/dmrg.py b/variational/dmrg.py
new file mode 100644
index 0000000..1805a91
--- /dev/null
+++ b/variational/dmrg.py
@@ -0,0 +1,76 @@
+from syngular.tensor import MatrixProductOperator, MatrixProductState
+
+from opt_einsum import contract
+
+class DMRG:
+
+
+    # Does not work with operator with input_shape != output_shape!
+    # To remedy use rho = AA^t the d ensity matrix of the operator
+
+    def solve(operator: MatrixProductOperator):
+
+        n = operator.sites_number
+        input_shape = operator.input_shape
+        bond_shape = operator.bond_shape
+        
+        state = MatrixProductState.random(input_shape=input_shape, bond_shape=bond_shape)
+        state.right_orthonormalization()
+
+        print('[DMRG] Initialization')
+
+        left_blocks = DMRG.__left_blocks(operator=operator, state=state)
+        right_blocks = DMRG.__right_blocks(operator=operator, state=state)
+
+        print([r.shape for r in right_blocks])
+        print([r.shape for r in left_blocks])
+
+
+    def __right_blocks(operator: MatrixProductOperator, state: MatrixProductState):
+        right_blocks = []
+        n = operator.sites_number
+
+        print('[DMRG] Right blocks')
+
+        for k in range(n-1,1,-1):
+
+            state_site = state.sites[k]
+            op_site = operator.sites[k]
+
+            struct = [
+                state_site, [1,3,2],
+                op_site, [4,6,3,5],
+                state_site, [7, 6, 8]
+            ]
+            if k != n-1: struct += [right_blocks[n-k-2], [2,5,8]]
+            struct += [[1, 4, 7]]
+            
+            block = contract(*struct)
+            right_blocks.append(block)
+        
+        return right_blocks
+
+
+    def __left_blocks(operator: MatrixProductOperator, state: MatrixProductState):
+        left_blocks = []
+        n = operator.sites_number
+
+        print('[DMRG] Left blocks')
+
+        for k in range(n-2):
+
+            state_site = state.sites[k]
+            op_site = operator.sites[k]
+
+            struct = [
+                state_site, [1,3,2],
+                op_site, [4,6,3,5],
+                state_site, [7, 6, 8]
+            ]
+            if k != 0: struct += [left_blocks[k-1], [1, 4, 7]]
+            struct += [[2, 5, 8]]
+            
+            block = contract(*struct)
+            left_blocks.append(block)
+        
+        return left_blocks
\ No newline at end of file
diff --git a/variational/lanczos.py b/variational/lanczos.py
new file mode 100644
index 0000000..4767f56
--- /dev/null
+++ b/variational/lanczos.py
@@ -0,0 +1,4 @@
+
+
+class Lanczos:
+    pass
\ No newline at end of file