diff --git a/ast/arrays.py b/ast/arrays.py
index 93b718bf6bb88eb6a4cbb8ff0fe8b706842615ea..6908bbe8b0e340bc4d019e00dc755bd3dca77476 100644
--- a/ast/arrays.py
+++ b/ast/arrays.py
@@ -126,6 +126,7 @@ class ArrayAccess:
         self.array = array
         self.indexes = [as_lit_ast(sim, index)]
         self.index = None
+        self.mutable = True
         self.generated = False
         self.check_and_set_index()
 
@@ -177,6 +178,9 @@ class ArrayAccess:
         return self.array.type()
         # return self.array.type() if self.index is None else Type_Array
 
+    def is_mutable(self):
+        return self.mutable
+
     def scope(self):
         if self.index is None:
             scope = None
diff --git a/ast/block.py b/ast/block.py
index f68f9811b0375698e42148f58c78862e2367c68b..de9e32bd9967ce26244aa4f3dc2f448ad5a85726 100644
--- a/ast/block.py
+++ b/ast/block.py
@@ -68,10 +68,8 @@ class Block:
         return fn(self)
 
     def merge_blocks(block1, block2):
-        assert isinstance(block1, Block), \
-            "First block type is not Block!"
-        assert isinstance(block2, Block), \
-            "Second block type is not Block!"
+        assert isinstance(block1, Block), "First block type is not Block!"
+        assert isinstance(block2, Block), "Second block type is not Block!"
         return Block(block1.sim, block1.statements() + block2.statements())
 
     def from_list(sim, block_list):
@@ -79,8 +77,7 @@ class Block:
         result_block = Block(sim, [])
 
         for block in block_list:
-            assert isinstance(block, Block), \
-                "Element in list is not Block!"
+            assert isinstance(block, Block), "Element in list is not Block!"
             result_block = Block.merge_blocks(result_block, block)
 
         return result_block
diff --git a/ast/cast.py b/ast/cast.py
index 970109e06a25556468d34530f076d30961a0cc3f..d219d961b21908f1350beeeaf8609276bf95659f 100644
--- a/ast/cast.py
+++ b/ast/cast.py
@@ -19,6 +19,9 @@ class Cast:
     def type(self):
         return self.cast_type
 
+    def is_mutable(self):
+        return self.expr.is_mutable()
+
     def scope(self):
         return self.expr.scope()
 
diff --git a/ast/expr.py b/ast/expr.py
index ac889de3b5c27616b250834f34a3b8557495f326..c7ec9f4ba3a429387d9a2624c0fd1053aaa978d0 100644
--- a/ast/expr.py
+++ b/ast/expr.py
@@ -21,6 +21,7 @@ class Expr:
         self.expr_type = Expr.infer_type(self.lhs, self.rhs, self.op)
         self.expr_scope = None
         self.vec_generated = []
+        self.mutable = self.lhs.is_mutable() or self.rhs.is_mutable() # Value can change accross references
         self.generated = False
 
     def __str__(self):
@@ -131,6 +132,9 @@ class Expr:
     def type(self):
         return self.expr_type
 
+    def is_mutable(self):
+        return self.mutable
+
     def scope(self):
         if self.expr_scope is None:
             lhs_scp = self.lhs.scope()
@@ -186,6 +190,7 @@ class ExprVec():
         self.expr = expr
         self.index = index
         self.expr_scope = None
+        self.mutable = self.expr.is_mutable() or self.index.is_mutable()
         self.lhs = (expr.lhs if not isinstance(expr.lhs, Expr)
                     else ExprVec(sim, expr.lhs, index))
         self.rhs = (expr.rhs if not isinstance(expr.rhs, Expr)
@@ -231,6 +236,9 @@ class ExprVec():
     def type(self):
         return Type_Float
 
+    def is_mutable(self):
+        return self.mutable
+
     def scope(self):
         if self.expr_scope is None:
             expr_scp = self.expr.scope()
diff --git a/ast/lit.py b/ast/lit.py
index 693c32d463970ddfa248a03019212b4ecb09ac2c..ff09f2397178eeceb8a0dbc8cf9ef96f2c6d5f38 100644
--- a/ast/lit.py
+++ b/ast/lit.py
@@ -48,6 +48,9 @@ class Lit:
     def type(self):
         return self.lit_type
 
+    def is_mutable(self):
+        return False
+
     def scope(self):
         return self.sim.global_scope
 
diff --git a/ast/loops.py b/ast/loops.py
index acf12d5fde7dee0678018d832f570e17ffc07392..817b6215f08671ed6ab58be5a74ed50b1a885605 100644
--- a/ast/loops.py
+++ b/ast/loops.py
@@ -20,6 +20,9 @@ class Iter():
     def type(self):
         return Type_Int
 
+    def is_mutable(self):
+        return False
+
     def scope(self):
         return self.loop.block
 
diff --git a/ast/math.py b/ast/math.py
index e234ff2b6d2bce07c9f9d680d1f8ec6462a9393c..bee39951a0894197903867afb999b4b86af66c32 100644
--- a/ast/math.py
+++ b/ast/math.py
@@ -12,6 +12,9 @@ class Sqrt:
     def type(self):
         return self.expr.type()
 
+    def is_mutable(self):
+        return self.expr.is_mutable()
+
     def scope(self):
         return self.expr.scope()
 
diff --git a/ast/properties.py b/ast/properties.py
index 64c042bedefa68431a3295aaa25edf60376b1064..3e261beb111a66f93cc82226472953e69fb01137 100644
--- a/ast/properties.py
+++ b/ast/properties.py
@@ -46,6 +46,7 @@ class Property:
         self.prop_layout = layout
         self.default_value = default
         self.volatile = volatile
+        self.mutable = True
         self.flattened = False
 
     def __str__(self):
@@ -63,6 +64,9 @@ class Property:
     def default(self):
         return self.default_value
 
+    def is_mutable(self):
+        return self.mutable
+
     def scope(self):
         return self.sim.global_scope
 
diff --git a/ast/sizeof.py b/ast/sizeof.py
index 917328a1b3863e6270ef8f3f480185fc724666b4..df171a81a50fb615f02d74ee4da6314b8eb670da 100644
--- a/ast/sizeof.py
+++ b/ast/sizeof.py
@@ -13,6 +13,9 @@ class Sizeof:
     def type(self):
         return Type_Int
 
+    def is_mutable(self):
+        return False
+
     def scope(self):
         return self.sim.global_scope
 
diff --git a/ast/transform.py b/ast/transform.py
index bdbb99756e334dc646b410484de7fbeb7a57e636..1dad678fe46eee9210f9c9ac3c508814aadbfd50 100644
--- a/ast/transform.py
+++ b/ast/transform.py
@@ -22,7 +22,8 @@ class Transform:
                 item = [f for f in Transform.flattened_list if
                         f[0] == ast.expr.lhs and
                         f[1] == ast.index and
-                        f[2] == ast.expr.rhs]
+                        f[2] == ast.expr.rhs and
+                        not ast.expr.rhs.is_mutable()]
                 if item:
                     return item[0][3]
 
diff --git a/ast/variables.py b/ast/variables.py
index 501f1c92f31922e9c6c976650fcbb4c41bf31cef..9a075eba90d8eda9ce78be3ec3a422a2b28a902a 100644
--- a/ast/variables.py
+++ b/ast/variables.py
@@ -29,6 +29,7 @@ class Var:
         self.var_name = var_name
         self.var_type = var_type
         self.var_init_value = init_value
+        self.mutable = True
         self.var_bonded_arrays = []
 
     def __str__(self):
@@ -94,6 +95,9 @@ class Var:
     def bonded_arrays(self):
         return self.var_bonded_arrays
 
+    def is_mutable(self):
+        return self.mutable
+
     def scope(self):
         return self.sim.global_scope
 
diff --git a/data/minimd_setup_4x4x4.input b/data/minimd_setup_4x4x4.input
new file mode 100644
index 0000000000000000000000000000000000000000..a31ad8693606608786a77d7678f2ef286705d13f
--- /dev/null
+++ b/data/minimd_setup_4x4x4.input
@@ -0,0 +1,257 @@
+-2.8,9.5184,-2.8,9.5184,-2.8,9.5184
+1,0,0,0,-1.13216,0.0489689,-2.03943
+1,1.6796,0,0,0.689294,0.210488,-1.97539
+1,3.3592,0,0,-1.6486,0.372007,-1.91135
+1,5.0388,0,0,0.172855,0.533527,-1.84732
+1,0.8398,0.8398,0,-1.25431,-1.30387,-1.75126
+1,2.5194,0.8398,0,0.567145,-1.14235,-1.68722
+1,4.199,0.8398,0,-1.77075,-0.980834,-1.62318
+1,5.8786,0.8398,0,0.0507054,-0.819315,-1.55914
+1,0,1.6796,0,0.961435,1.34112,-1.52712
+1,1.6796,1.6796,0,-1.37646,1.50264,-1.46308
+1,3.3592,1.6796,0,0.444995,1.66416,-1.39904
+1,5.0388,1.6796,0,-1.8929,1.82568,-1.335
+1,0.8398,2.5194,0,0.839286,-0.011719,-1.23895
+1,2.5194,2.5194,0,-1.49861,0.1498,-1.17491
+1,4.199,2.5194,0,0.322846,0.311319,-1.11087
+1,5.8786,2.5194,0,-2.01505,0.472839,-1.04683
+1,0,3.3592,0,-1.10432,-1.52608,-1.01481
+1,1.6796,3.3592,0,0.717136,-1.36456,-0.95077
+1,3.3592,3.3592,0,-1.62076,-1.20304,-0.886731
+1,5.0388,3.3592,0,0.200697,-1.04152,-0.822692
+1,0.8398,4.199,0,-1.22647,1.28043,-0.726633
+1,2.5194,4.199,0,0.594987,1.44195,-0.662594
+1,4.199,4.199,0,-1.74291,1.60347,-0.598555
+1,5.8786,4.199,0,0.078548,1.76499,-0.534516
+1,0,5.0388,0,0.989277,-0.233926,-0.502497
+1,1.6796,5.0388,0,-1.34862,-0.0724069,-0.438458
+1,3.3592,5.0388,0,0.472838,0.0891123,-0.374419
+1,5.0388,5.0388,0,-1.86506,0.250631,-0.31038
+1,0.8398,5.8786,0,0.867128,-1.58677,-0.214321
+1,2.5194,5.8786,0,-1.47077,-1.42525,-0.150282
+1,4.199,5.8786,0,0.350689,-1.26373,-0.0862431
+1,5.8786,5.8786,0,-1.98721,-1.10221,-0.0222041
+1,0.8398,0,0.8398,-0.16575,-0.940691,0.0418349
+1,2.5194,0,0.8398,1.65571,-0.779172,0.105874
+1,4.199,0,0.8398,-0.68219,-0.617653,0.169913
+1,5.8786,0,0.8398,1.13927,-0.456133,0.233952
+1,0,0.8398,0.8398,2.05,1.7043,0.265971
+1,1.6796,0.8398,0.8398,-0.287899,1.86582,0.33001
+1,3.3592,0.8398,0.8398,1.53356,2.02734,0.394049
+1,5.0388,0.8398,0.8398,-0.804339,-1.97049,0.458088
+1,0.8398,1.6796,0.8398,1.92785,0.351463,0.554147
+1,2.5194,1.6796,0.8398,-0.410049,0.512982,0.618186
+1,4.199,1.6796,0.8398,1.41141,0.674501,0.682225
+1,5.8786,1.6796,0.8398,-0.926488,0.83602,0.746264
+1,0,2.5194,0.8398,-0.0157586,-1.1629,0.778283
+1,1.6796,2.5194,0.8398,1.8057,-1.00138,0.842322
+1,3.3592,2.5194,0.8398,-0.532198,-0.83986,0.906361
+1,5.0388,2.5194,0.8398,1.28926,-0.67834,0.9704
+1,0.8398,3.3592,0.8398,-0.137908,1.64362,1.06646
+1,2.5194,3.3592,0.8398,1.68355,1.80514,1.1305
+1,4.199,3.3592,0.8398,-0.654347,1.96665,1.19454
+1,5.8786,3.3592,0.8398,1.16711,-2.03118,1.25858
+1,0,4.199,0.8398,2.07784,0.129256,1.2906
+1,1.6796,4.199,0.8398,-0.260057,0.290775,1.35463
+1,3.3592,4.199,0.8398,1.5614,0.452294,1.41867
+1,5.0388,4.199,0.8398,-0.776496,0.613813,1.48271
+1,0.8398,5.0388,0.8398,1.95569,-1.22359,1.57877
+1,2.5194,5.0388,0.8398,-0.382206,-1.06207,1.64281
+1,4.199,5.0388,0.8398,1.43925,-0.900548,1.70685
+1,5.8786,5.0388,0.8398,-0.898645,-0.739028,1.77089
+1,0,5.8786,0.8398,0.012084,1.42141,1.80291
+1,1.6796,5.8786,0.8398,1.83354,1.58293,1.86695
+1,3.3592,5.8786,0.8398,-0.504355,1.74445,1.93099
+1,5.0388,5.8786,0.8398,1.3171,1.90597,1.99502
+1,0,0,1.6796,-1.02079,-2.09187,2.05906
+1,1.6796,0,1.6796,0.800664,-1.93035,-2.03625
+1,3.3592,0,1.6796,-1.53723,-1.76883,-1.97221
+1,5.0388,0,1.6796,0.284225,-1.60731,-1.90818
+1,0.8398,0.8398,1.6796,-1.14294,0.714645,-1.81212
+1,2.5194,0.8398,1.6796,0.678515,0.876164,-1.74808
+1,4.199,0.8398,1.6796,-1.65938,1.03768,-1.68404
+1,5.8786,0.8398,1.6796,0.162076,1.1992,-1.62
+1,0,1.6796,1.6796,1.0728,-0.799716,-1.58798
+1,1.6796,1.6796,1.6796,-1.26509,-0.638197,-1.52394
+1,3.3592,1.6796,1.6796,0.556366,-0.476678,-1.4599
+1,5.0388,1.6796,1.6796,-1.78153,-0.315159,-1.39586
+1,0.8398,2.5194,1.6796,0.950656,2.0068,-1.29981
+1,2.5194,2.5194,1.6796,-1.38724,-1.99104,-1.23577
+1,4.199,2.5194,1.6796,0.434217,-1.82952,-1.17173
+1,5.8786,2.5194,1.6796,-1.90368,-1.668,-1.10769
+1,0,3.3592,1.6796,-0.992952,0.492437,-1.07567
+1,1.6796,3.3592,1.6796,0.828507,0.653957,-1.01163
+1,3.3592,3.3592,1.6796,-1.50939,0.815476,-0.947591
+1,5.0388,3.3592,1.6796,0.312067,0.976995,-0.883552
+1,0.8398,4.199,1.6796,-1.1151,-0.860404,-0.787493
+1,2.5194,4.199,1.6796,0.706357,-0.698885,-0.723454
+1,4.199,4.199,1.6796,-1.63154,-0.537366,-0.659415
+1,5.8786,4.199,1.6796,0.189918,-0.375847,-0.595376
+1,0,5.0388,1.6796,1.10065,1.78459,-0.563357
+1,1.6796,5.0388,1.6796,-1.23725,1.94611,-0.499318
+1,3.3592,5.0388,1.6796,0.584208,-2.05173,-0.435279
+1,5.0388,5.0388,1.6796,-1.75369,-1.89021,-0.37124
+1,0.8398,5.8786,1.6796,0.978498,0.43175,-0.275181
+1,2.5194,5.8786,1.6796,-1.3594,0.593269,-0.211142
+1,4.199,5.8786,1.6796,0.462059,0.754788,-0.147103
+1,5.8786,5.8786,1.6796,-1.87584,0.916307,-0.0830645
+1,0.8398,0,2.5194,-0.05438,1.07783,-0.0190255
+1,2.5194,0,2.5194,1.76708,1.23935,0.0450135
+1,4.199,0,2.5194,-0.570819,1.40086,0.109053
+1,5.8786,0,2.5194,1.25064,1.56238,0.173092
+1,0,0.8398,2.5194,-1.99799,-0.436534,0.205111
+1,1.6796,0.8398,2.5194,-0.176529,-0.275015,0.26915
+1,3.3592,0.8398,2.5194,1.64493,-0.113496,0.333189
+1,5.0388,0.8398,2.5194,-0.692968,0.0480232,0.397228
+1,0.8398,1.6796,2.5194,2.03922,-1.78938,0.493287
+1,2.5194,1.6796,2.5194,-0.298678,-1.62786,0.557326
+1,4.199,1.6796,2.5194,1.52278,-1.46634,0.621365
+1,5.8786,1.6796,2.5194,-0.815118,-1.30482,0.685404
+1,0,2.5194,2.5194,0.0956117,0.855619,0.717423
+1,1.6796,2.5194,2.5194,1.91707,1.01714,0.781462
+1,3.3592,2.5194,2.5194,-0.420828,1.17866,0.845501
+1,5.0388,2.5194,2.5194,1.40063,1.34018,0.90954
+1,0.8398,3.3592,2.5194,-0.0265375,-0.497222,1.0056
+1,2.5194,3.3592,2.5194,1.79492,-0.335703,1.06964
+1,4.199,3.3592,2.5194,-0.542977,-0.174184,1.13368
+1,5.8786,3.3592,2.5194,1.27848,-0.0126647,1.19772
+1,0,4.199,2.5194,-1.97015,-2.01158,1.22973
+1,1.6796,4.199,2.5194,-0.148687,-1.85006,1.29377
+1,3.3592,4.199,2.5194,1.67277,-1.68854,1.35781
+1,5.0388,4.199,2.5194,-0.665126,-1.52703,1.42185
+1,0.8398,5.0388,2.5194,2.06706,0.794931,1.51791
+1,2.5194,5.0388,2.5194,-0.270836,0.956451,1.58195
+1,4.199,5.0388,2.5194,1.55062,1.11797,1.64599
+1,5.8786,5.0388,2.5194,-0.787275,1.27949,1.71003
+1,0,5.8786,2.5194,0.123454,-0.719429,1.74205
+1,1.6796,5.8786,2.5194,1.94491,-0.55791,1.80609
+1,3.3592,5.8786,2.5194,-0.392985,-0.396391,1.87012
+1,5.0388,5.8786,2.5194,1.42847,-0.234872,1.93416
+1,0,0,3.3592,-0.909424,-0.0733526,1.9982
+1,1.6796,0,3.3592,0.912034,0.0881666,2.06224
+1,3.3592,0,3.3592,-1.42586,0.249686,-2.03308
+1,5.0388,0,3.3592,0.395595,0.411205,-1.96904
+1,0.8398,0.8398,3.3592,-1.03157,-1.42619,-1.87298
+1,2.5194,0.8398,3.3592,0.789885,-1.26468,-1.80894
+1,4.199,0.8398,3.3592,-1.54801,-1.10316,-1.7449
+1,5.8786,0.8398,3.3592,0.273446,-0.941637,-1.68086
+1,0,1.6796,3.3592,1.18418,1.2188,-1.64884
+1,1.6796,1.6796,3.3592,-1.15372,1.38032,-1.5848
+1,3.3592,1.6796,3.3592,0.667736,1.54184,-1.52076
+1,5.0388,1.6796,3.3592,-1.67016,1.70336,-1.45672
+1,0.8398,2.5194,3.3592,1.06203,-0.134041,-1.36067
+1,2.5194,2.5194,3.3592,-1.27587,0.0274787,-1.29663
+1,4.199,2.5194,3.3592,0.545587,0.188998,-1.23259
+1,5.8786,2.5194,3.3592,-1.79231,0.350517,-1.16855
+1,0,3.3592,3.3592,-0.881582,-1.6484,-1.13653
+1,1.6796,3.3592,3.3592,0.939877,-1.48688,-1.07249
+1,3.3592,3.3592,3.3592,-1.39802,-1.32536,-1.00845
+1,5.0388,3.3592,3.3592,0.423438,-1.16384,-0.944412
+1,0.8398,4.199,3.3592,-1.00373,1.15811,-0.848354
+1,2.5194,4.199,3.3592,0.817728,1.31963,-0.784315
+1,4.199,4.199,3.3592,-1.52017,1.48115,-0.720276
+1,5.8786,4.199,3.3592,0.301289,1.64267,-0.656237
+1,0,5.0388,3.3592,1.21202,-0.356248,-0.624217
+1,1.6796,5.0388,3.3592,-1.12588,-0.194728,-0.560178
+1,3.3592,5.0388,3.3592,0.695579,-0.0332092,-0.496139
+1,5.0388,5.0388,3.3592,-1.64232,0.12831,-0.4321
+1,0.8398,5.8786,3.3592,1.08987,-1.70909,-0.336042
+1,2.5194,5.8786,3.3592,-1.24803,-1.54757,-0.272003
+1,4.199,5.8786,3.3592,0.573429,-1.38605,-0.207964
+1,5.8786,5.8786,3.3592,-1.76447,-1.22453,-0.143925
+1,0.8398,0,4.199,0.0569903,-1.06301,-0.0798859
+1,2.5194,0,4.199,1.87845,-0.901493,-0.0158469
+1,4.199,0,4.199,-0.459449,-0.739974,0.0481921
+1,5.8786,0,4.199,1.36201,-0.578455,0.112231
+1,0,0.8398,4.199,-1.88662,1.58198,0.144251
+1,1.6796,0.8398,4.199,-0.0651589,1.7435,0.20829
+1,3.3592,0.8398,4.199,1.7563,1.90502,0.272329
+1,5.0388,0.8398,4.199,-0.581598,-2.09282,0.336368
+1,0.8398,1.6796,4.199,-2.00877,0.229141,0.432426
+1,2.5194,1.6796,4.199,-0.187308,0.390661,0.496465
+1,4.199,1.6796,4.199,1.63415,0.55218,0.560504
+1,5.8786,1.6796,4.199,-0.703747,0.713699,0.624543
+1,0,2.5194,4.199,0.206982,-1.28522,0.656563
+1,1.6796,2.5194,4.199,2.02844,-1.1237,0.720602
+1,3.3592,2.5194,4.199,-0.309457,-0.962181,0.784641
+1,5.0388,2.5194,4.199,1.512,-0.800662,0.84868
+1,0.8398,3.3592,4.199,0.0848328,1.5213,0.944738
+1,2.5194,3.3592,4.199,1.90629,1.68281,1.00878
+1,4.199,3.3592,4.199,-0.431606,1.84433,1.07282
+1,5.8786,3.3592,4.199,1.38985,2.00585,1.13686
+1,0,4.199,4.199,-1.85877,0.00693419,1.16887
+1,1.6796,4.199,4.199,-0.0373163,0.168453,1.23291
+1,3.3592,4.199,4.199,1.78414,0.329973,1.29695
+1,5.0388,4.199,4.199,-0.553756,0.491492,1.36099
+1,0.8398,5.0388,4.199,-1.98092,-1.34591,1.45705
+1,2.5194,5.0388,4.199,-0.159465,-1.18439,1.52109
+1,4.199,5.0388,4.199,1.66199,-1.02287,1.58513
+1,5.8786,5.0388,4.199,-0.675905,-0.86135,1.64917
+1,0,5.8786,4.199,0.234825,1.29909,1.68119
+1,1.6796,5.8786,4.199,2.05628,1.46061,1.74523
+1,3.3592,5.8786,4.199,-0.281615,1.62213,1.80926
+1,5.0388,5.8786,4.199,1.53984,1.78365,1.8733
+1,0,0,5.0388,-0.798054,1.94516,1.93734
+1,1.6796,0,5.0388,1.0234,-2.05267,2.00138
+1,3.3592,0,5.0388,-1.31449,-1.89115,2.06542
+1,5.0388,0,5.0388,0.506965,-1.72963,-2.0299
+1,0.8398,0.8398,5.0388,-0.920203,0.592323,-1.93384
+1,2.5194,0.8398,5.0388,0.901256,0.753842,-1.8698
+1,4.199,0.8398,5.0388,-1.43664,0.915362,-1.80576
+1,5.8786,0.8398,5.0388,0.384816,1.07688,-1.74172
+1,0,1.6796,5.0388,1.29555,-0.922038,-1.7097
+1,1.6796,1.6796,5.0388,-1.04235,-0.760518,-1.64566
+1,3.3592,1.6796,5.0388,0.779106,-0.598999,-1.58162
+1,5.0388,1.6796,5.0388,-1.55879,-0.43748,-1.51758
+1,0.8398,2.5194,5.0388,1.1734,1.88448,-1.42153
+1,2.5194,2.5194,5.0388,-1.1645,2.046,-1.35749
+1,4.199,2.5194,5.0388,0.656957,-1.95184,-1.29345
+1,5.8786,2.5194,5.0388,-1.68094,-1.79032,-1.22941
+1,0,3.3592,5.0388,-0.770211,0.370116,-1.19739
+1,1.6796,3.3592,5.0388,1.05125,0.531635,-1.13335
+1,3.3592,3.3592,5.0388,-1.28665,0.693154,-1.06931
+1,5.0388,3.3592,5.0388,0.534808,0.854674,-1.00527
+1,0.8398,4.199,5.0388,-0.89236,-0.982726,-0.909214
+1,2.5194,4.199,5.0388,0.929098,-0.821206,-0.845175
+1,4.199,4.199,5.0388,-1.4088,-0.659687,-0.781136
+1,5.8786,4.199,5.0388,0.412659,-0.498168,-0.717097
+1,0,5.0388,5.0388,1.32339,1.66227,-0.685078
+1,1.6796,5.0388,5.0388,-1.01451,1.82379,-0.621039
+1,3.3592,5.0388,5.0388,0.806949,1.98531,-0.557
+1,5.0388,5.0388,5.0388,-1.53095,-2.01253,-0.492961
+1,0.8398,5.8786,5.0388,1.20124,0.309428,-0.396902
+1,2.5194,5.8786,5.0388,-1.13666,0.470947,-0.332863
+1,4.199,5.8786,5.0388,0.6848,0.632467,-0.268824
+1,5.8786,5.8786,5.0388,-1.6531,0.793986,-0.204785
+1,0.8398,0,5.8786,0.168361,0.955505,-0.140746
+1,2.5194,0,5.8786,1.98982,1.11702,-0.0767072
+1,4.199,0,5.8786,-0.348079,1.27854,-0.0126682
+1,5.8786,0,5.8786,1.47338,1.44006,0.0513707
+1,0,0.8398,5.8786,-1.77525,-0.558856,0.0833902
+1,1.6796,0.8398,5.8786,0.0462114,-0.397337,0.147429
+1,3.3592,0.8398,5.8786,1.86767,-0.235817,0.211468
+1,5.0388,0.8398,5.8786,-0.470228,-0.0742982,0.275507
+1,0.8398,1.6796,5.8786,-1.8974,-1.9117,0.371566
+1,2.5194,1.6796,5.8786,-0.0759377,-1.75018,0.435605
+1,4.199,1.6796,5.8786,1.74552,-1.58866,0.499644
+1,5.8786,1.6796,5.8786,-0.592377,-1.42714,0.563683
+1,0,2.5194,5.8786,0.318352,0.733298,0.595702
+1,1.6796,2.5194,5.8786,-2.01955,0.894817,0.659741
+1,3.3592,2.5194,5.8786,-0.198087,1.05634,0.72378
+1,5.0388,2.5194,5.8786,1.62337,1.21786,0.787819
+1,0.8398,3.3592,5.8786,0.196203,-0.619544,0.883878
+1,2.5194,3.3592,5.8786,2.01766,-0.458025,0.947917
+1,4.199,3.3592,5.8786,-0.320236,-0.296505,1.01196
+1,5.8786,3.3592,5.8786,1.50122,-0.134986,1.07599
+1,0,4.199,5.8786,-1.7474,2.02545,1.10801
+1,1.6796,4.199,5.8786,0.074054,-1.97239,1.17205
+1,3.3592,4.199,5.8786,1.89551,-1.81087,1.23609
+1,5.0388,4.199,5.8786,-0.442385,-1.64935,1.30013
+1,0.8398,5.0388,5.8786,-1.86955,0.67261,1.39619
+1,2.5194,5.0388,5.8786,-0.0480952,0.834129,1.46023
+1,4.199,5.0388,5.8786,1.77336,0.995648,1.52427
+1,5.8786,5.0388,5.8786,-0.564534,1.15717,1.58831
+1,0,5.8786,5.8786,0.346195,-0.841751,1.62033
+1,1.6796,5.8786,5.8786,-1.9917,-0.680232,1.68437
+1,3.3592,5.8786,5.8786,-0.170244,-0.518712,1.7484
+1,5.0388,5.8786,5.8786,1.65121,-0.357193,1.81244
diff --git a/output/.gitignore b/output/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..5e7d2734cfc60289debf74293817c0a8f572ff32
--- /dev/null
+++ b/output/.gitignore
@@ -0,0 +1,4 @@
+# Ignore everything in this directory
+*
+# Except this file
+!.gitignore
diff --git a/particle.py b/particle.py
index 7a4377eb39cf9933a8c41cb3f8816a5885059878..d9b95bf4a97ba01b139f976444d85b500c859864 100644
--- a/particle.py
+++ b/particle.py
@@ -14,9 +14,10 @@ position = psim.add_vector_property('position')
 velocity = psim.add_vector_property('velocity')
 force = psim.add_vector_property('force', vol=True)
 
-grid = psim.grid_3d(0.0, 16.0, 0.0, 16.0, 0.0, 16.0)
-psim.create_particle_lattice(grid, spacing=[1.2, 1.2, 1.2])
-psim.create_cell_lists(grid, 2.8, 2.8)
+#grid = psim.grid_3d(0.0, 4.0, 0.0, 4.0, 0.0, 4.0)
+#psim.create_particle_lattice(grid, spacing=[0.82323, 0.82323, 0.82323])
+psim.from_file("data/minimd_setup_4x4x4.input", ['mass', 'position', 'velocity'])
+psim.create_cell_lists(2.8, 2.8)
 psim.periodic(2.8)
 psim.vtk_output("output/test")
 
diff --git a/sim/particle_simulation.py b/sim/particle_simulation.py
index 3f3c485d5fa405b7f626bf558c074976f2d1aec5..ce0ba31f55cf86cd666f7413ae7b16bf054ed3f6 100644
--- a/sim/particle_simulation.py
+++ b/sim/particle_simulation.py
@@ -14,6 +14,7 @@ from sim.kernel_wrapper import KernelWrapper
 from sim.lattice import ParticleLattice
 from sim.pbc import PBC, UpdatePBC, EnforcePBC, SetupPBC
 from sim.properties import PropertiesAlloc, PropertiesResetVolatile
+from sim.read_from_file import ReadFromFile
 from sim.setup_wrapper import SetupWrapper
 from sim.timestep import Timestep
 from sim.variables import VariablesDecl
@@ -103,8 +104,14 @@ class ParticleSimulation:
         lattice = ParticleLattice(self, grid, spacing, props, positions)
         self.setups.add_setup_block(lattice.lower())
 
-    def create_cell_lists(self, grid, spacing, cutoff_radius):
-        self.cell_lists = CellLists(self, grid, spacing, cutoff_radius)
+    def from_file(self, filename, prop_names):
+        props = [self.property(prop_name) for prop_name in prop_names]
+        read_object = ReadFromFile(self, filename, props)
+        self.setups.add_setup_block(read_object.lower())
+        self.grid = read_object.grid
+
+    def create_cell_lists(self, spacing, cutoff_radius):
+        self.cell_lists = CellLists(self, self.grid, spacing, cutoff_radius)
         return self.cell_lists
 
     def periodic(self, cutneigh, flags=[1, 1, 1]):
@@ -168,8 +175,8 @@ class ParticleSimulation:
     def generate(self):
         timestep = Timestep(self, self.ntimesteps, [
             (EnforcePBC(self.pbc).lower(), 20),
-            (SetupPBC(self.pbc).lower(), 20),
-            UpdatePBC(self.pbc).lower(),
+            #(SetupPBC(self.pbc).lower(), 20),
+            #UpdatePBC(self.pbc).lower(),
             (CellListsBuild(self.cell_lists).lower(), 20),
             PropertiesResetVolatile(self).lower(),
             self.kernels.lower()
@@ -180,6 +187,10 @@ class ParticleSimulation:
         body = Block.from_list(self, [
             CellListsStencilBuild(self.cell_lists).lower(),
             self.setups.lower(),
+            EnforcePBC(self.pbc).lower(),
+            #SetupPBC(self.pbc).lower(),
+            #UpdatePBC(self.pbc).lower(),
+            CellListsBuild(self.cell_lists).lower(),
             Block(self, VTKWrite(self, self.vtk_file, 0)),
             timestep.as_block()
         ])
diff --git a/sim/read_from_file.py b/sim/read_from_file.py
new file mode 100644
index 0000000000000000000000000000000000000000..42608a1161e3e4baa3011939ca3d3527e3875036
--- /dev/null
+++ b/sim/read_from_file.py
@@ -0,0 +1,49 @@
+from ast.data_types import Type_Int, Type_Float, Type_Vector
+from ast.loops import For
+from sim.grid import Grid
+
+class ReadFromFile():
+    def __init__(self, sim, filename, props):
+        self.sim = sim
+        self.filename = filename
+        self.props = props
+        self.grid = None
+
+    def lower(self):
+        ndims = None
+        nlocal = self.sim.nlocal
+        line_number = 0
+
+        self.sim.clear_block()
+        with open(self.filename, "r") as fp:
+            for line in fp:
+                current_data = line.split(',')
+                if line_number == 0:
+                    assert len(current_data) % 2 == 0, "Number of dimensions in header must be even!"
+                    ndims = int(len(current_data) / 2)
+                    config = [[float(current_data[d * 2]), float(current_data[d * 2 + 1])] for d in range(0, ndims)]
+                    self.grid = Grid(self.sim, config)
+                else:
+                    i = 0
+                    for p in self.props:
+                        if p.type() == Type_Vector:
+                            for d in range(0, ndims):
+                                p[nlocal][d].set(float(current_data[i + d]))
+
+                            i += ndims
+
+                        else:
+                            if p.type() == Type_Int:
+                                p[nlocal].set(int(current_data[i]))
+                            elif p.type() == Type_Float:
+                                p[nlocal].set(float(current_data[i]))
+                            else:
+                                raise Exception(f"Invalid property type at line {line_number}!")
+
+                            i += 1
+
+                    nlocal.set(nlocal + 1)
+
+                line_number += 1
+
+        return self.sim.block
diff --git a/sim/timestep.py b/sim/timestep.py
index 222fda1cb66ee6d11218e319e7fa1ddd403efe22..3599933bbd2f34a2217418f0f05ebc56dd31b00a 100644
--- a/sim/timestep.py
+++ b/sim/timestep.py
@@ -8,7 +8,7 @@ class Timestep:
     def __init__(self, sim, nsteps, item_list=None):
         self.sim = sim
         self.block = Block(sim, [])
-        self.timestep_loop = For(sim, 0, nsteps, self.block)
+        self.timestep_loop = For(sim, 1, nsteps + 1, self.block)
 
         if item_list is not None:
             for item in item_list:
@@ -24,9 +24,7 @@ class Timestep:
         assert exec_every >= 0, \
             "exec_every parameter must be higher or equal than zero!"
 
-        stmts = (item if not isinstance(item, Block)
-                 else item.statements())
-
+        stmts = item if not isinstance(item, Block) else item.statements()
         ts = self.timestep_loop.iter()
         if exec_every > 0:
             self.block.add_statement(