diff --git a/LICENSE b/LICENSE
new file mode 100644
index 0000000000000000000000000000000000000000..5f4049087480f2bbf810dc15d4f8fa4dbe0ee247
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2021 Rafael Ravedutti Lucio Machado
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
diff --git a/examples/kernels b/examples/kernels
new file mode 100644
index 0000000000000000000000000000000000000000..6fcb7d8905a05538e27c395fed908b180b2071bf
--- /dev/null
+++ b/examples/kernels
@@ -0,0 +1,834 @@
+digraph AST {
+	node [color=lightblue2 style=filled]
+	size="6,6"
+	n139866028692368 [label=Block]
+	n139866028640632 [label=ParticleFor]
+	n139866028692368 -> n139866028640632
+	n139866028694664 [label=ParticleFor]
+	n139866028692368 -> n139866028694664
+	n139866028640632 [label=ParticleFor]
+	n139866028694720 [label="Iter(0)"]
+	n139866028640632 -> n139866028694720
+	n139866028694888 [label=Block]
+	n139866028640632 -> n139866028694888
+	n139866028694776 [label=0]
+	n139866028640632 -> n139866028694776
+	n139866028694832 [label=0]
+	n139866028640632 -> n139866028694832
+	n139866028694720 [label="Iter(0)"]
+	n139866028694888 [label=Block]
+	n139866028695056 [label=For]
+	n139866028694888 -> n139866028695056
+	n139866028695336 [label="*"]
+	n139866028694720 [label="Iter(0)"]
+	n139866028695336 -> n139866028694720
+	n139866028691808 [label=neighborlist_capacity]
+	n139866028695336 -> n139866028691808
+	n139866028691808 [label=neighborlist_capacity]
+	n139866028692480 [label=PropertyAccess]
+	n139866028639904 [label=position]
+	n139866028692480 -> n139866028639904
+	n139866028694720 [label="Iter(0)"]
+	n139866028692480 -> n139866028694720
+	n139866028388816 [label="*"]
+	n139866028692480 -> n139866028388816
+	n139866028390552 [label="+"]
+	n139866028692480 -> n139866028390552
+	n139866028392120 [label="+"]
+	n139866028692480 -> n139866028392120
+	n139866028639904 [label=position]
+	n139866028388816 [label="*"]
+	n139866028694720 [label="Iter(0)"]
+	n139866028388816 -> n139866028694720
+	n139866028388872 [label=3]
+	n139866028388816 -> n139866028388872
+	n139866028388872 [label=3]
+	n139866028390552 [label="+"]
+	n139866028390384 [label="*"]
+	n139866028390552 -> n139866028390384
+	n139866028390608 [label=1]
+	n139866028390552 -> n139866028390608
+	n139866028390384 [label="*"]
+	n139866028694720 [label="Iter(0)"]
+	n139866028390384 -> n139866028694720
+	n139866028390440 [label=3]
+	n139866028390384 -> n139866028390440
+	n139866028390440 [label=3]
+	n139866028390608 [label=1]
+	n139866028392120 [label="+"]
+	n139866028391952 [label="*"]
+	n139866028392120 -> n139866028391952
+	n139866028392176 [label=2]
+	n139866028392120 -> n139866028392176
+	n139866028391952 [label="*"]
+	n139866028694720 [label="Iter(0)"]
+	n139866028391952 -> n139866028694720
+	n139866028392008 [label=3]
+	n139866028391952 -> n139866028392008
+	n139866028392008 [label=3]
+	n139866028392176 [label=2]
+	n139866028216728 [label="*"]
+	n139866028694720 [label="Iter(0)"]
+	n139866028216728 -> n139866028694720
+	n139866028216784 [label=3]
+	n139866028216728 -> n139866028216784
+	n139866028216784 [label=3]
+	n139866028217512 [label="*"]
+	n139866028694720 [label="Iter(0)"]
+	n139866028217512 -> n139866028694720
+	n139866028217568 [label=3]
+	n139866028217512 -> n139866028217568
+	n139866028217568 [label=3]
+	n139866028218408 [label="*"]
+	n139866028694720 [label="Iter(0)"]
+	n139866028218408 -> n139866028694720
+	n139866028218464 [label=3]
+	n139866028218408 -> n139866028218464
+	n139866028218464 [label=3]
+	n139866028218576 [label="+"]
+	n139866028218408 [label="*"]
+	n139866028218576 -> n139866028218408
+	n139866028218632 [label=1]
+	n139866028218576 -> n139866028218632
+	n139866028218632 [label=1]
+	n139866028219192 [label="*"]
+	n139866028694720 [label="Iter(0)"]
+	n139866028219192 -> n139866028694720
+	n139866028219248 [label=3]
+	n139866028219192 -> n139866028219248
+	n139866028219248 [label=3]
+	n139866028219360 [label="+"]
+	n139866028219192 [label="*"]
+	n139866028219360 -> n139866028219192
+	n139866028219416 [label=1]
+	n139866028219360 -> n139866028219416
+	n139866028219416 [label=1]
+	n139866028220200 [label="*"]
+	n139866028694720 [label="Iter(0)"]
+	n139866028220200 -> n139866028694720
+	n139866028220256 [label=3]
+	n139866028220200 -> n139866028220256
+	n139866028220256 [label=3]
+	n139866028220368 [label="+"]
+	n139866028220200 [label="*"]
+	n139866028220368 -> n139866028220200
+	n139866028355656 [label=2]
+	n139866028220368 -> n139866028355656
+	n139866028355656 [label=2]
+	n139866028356216 [label="*"]
+	n139866028694720 [label="Iter(0)"]
+	n139866028356216 -> n139866028694720
+	n139866028356272 [label=3]
+	n139866028356216 -> n139866028356272
+	n139866028356272 [label=3]
+	n139866028356384 [label="+"]
+	n139866028356216 [label="*"]
+	n139866028356384 -> n139866028356216
+	n139866028356440 [label=2]
+	n139866028356384 -> n139866028356440
+	n139866028356440 [label=2]
+	n139866028359408 [label="*"]
+	n139866028694720 [label="Iter(0)"]
+	n139866028359408 -> n139866028694720
+	n139866028359464 [label=3]
+	n139866028359408 -> n139866028359464
+	n139866028359464 [label=3]
+	n139866028388536 [label="*"]
+	n139866028694720 [label="Iter(0)"]
+	n139866028388536 -> n139866028694720
+	n139866028388592 [label=3]
+	n139866028388536 -> n139866028388592
+	n139866028388592 [label=3]
+	n139866028389712 [label="*"]
+	n139866028694720 [label="Iter(0)"]
+	n139866028389712 -> n139866028694720
+	n139866028389768 [label=3]
+	n139866028389712 -> n139866028389768
+	n139866028389768 [label=3]
+	n139866028389880 [label="+"]
+	n139866028389712 [label="*"]
+	n139866028389880 -> n139866028389712
+	n139866028389936 [label=1]
+	n139866028389880 -> n139866028389936
+	n139866028389936 [label=1]
+	n139866028390104 [label="*"]
+	n139866028694720 [label="Iter(0)"]
+	n139866028390104 -> n139866028694720
+	n139866028390160 [label=3]
+	n139866028390104 -> n139866028390160
+	n139866028390160 [label=3]
+	n139866028390272 [label="+"]
+	n139866028390104 [label="*"]
+	n139866028390272 -> n139866028390104
+	n139866028390328 [label=1]
+	n139866028390272 -> n139866028390328
+	n139866028390328 [label=1]
+	n139866028391280 [label="*"]
+	n139866028694720 [label="Iter(0)"]
+	n139866028391280 -> n139866028694720
+	n139866028391336 [label=3]
+	n139866028391280 -> n139866028391336
+	n139866028391336 [label=3]
+	n139866028391448 [label="+"]
+	n139866028391280 [label="*"]
+	n139866028391448 -> n139866028391280
+	n139866028391504 [label=2]
+	n139866028391448 -> n139866028391504
+	n139866028391504 [label=2]
+	n139866028391672 [label="*"]
+	n139866028694720 [label="Iter(0)"]
+	n139866028391672 -> n139866028694720
+	n139866028391728 [label=3]
+	n139866028391672 -> n139866028391728
+	n139866028391728 [label=3]
+	n139866028391840 [label="+"]
+	n139866028391672 [label="*"]
+	n139866028391840 -> n139866028391672
+	n139866028391896 [label=2]
+	n139866028391840 -> n139866028391896
+	n139866028391896 [label=2]
+	n139866028695056 [label=For]
+	n139866028695112 [label="Iter(1)"]
+	n139866028695056 -> n139866028695112
+	n139866028695224 [label=Block]
+	n139866028695056 -> n139866028695224
+	n139866028695168 [label=0]
+	n139866028695056 -> n139866028695168
+	n139866028695000 [label=ArrayAccess]
+	n139866028695056 -> n139866028695000
+	n139866028695112 [label="Iter(1)"]
+	n139866028695224 [label=Block]
+	n139866028357336 [label=Filter]
+	n139866028695224 -> n139866028357336
+	n139866028695448 [label="+"]
+	n139866028695336 [label="*"]
+	n139866028695448 -> n139866028695336
+	n139866028695112 [label="Iter(1)"]
+	n139866028695448 -> n139866028695112
+	n139866028216448 [label=PropertyAccess]
+	n139866028639904 [label=position]
+	n139866028216448 -> n139866028639904
+	n139866028695280 [label=ArrayAccess]
+	n139866028216448 -> n139866028695280
+	n139866028389152 [label="*"]
+	n139866028216448 -> n139866028389152
+	n139866028390888 [label="+"]
+	n139866028216448 -> n139866028390888
+	n139866028425288 [label="+"]
+	n139866028216448 -> n139866028425288
+	n139866028695280 [label=ArrayAccess]
+	n139866028691864 [label=neighborlists]
+	n139866028695280 -> n139866028691864
+	n139866028695448 [label="+"]
+	n139866028695280 -> n139866028695448
+	n139866028691864 [label=neighborlists]
+	n139866028389152 [label="*"]
+	n139866028695280 [label=ArrayAccess]
+	n139866028389152 -> n139866028695280
+	n139866028389208 [label=3]
+	n139866028389152 -> n139866028389208
+	n139866028389208 [label=3]
+	n139866028390888 [label="+"]
+	n139866028390720 [label="*"]
+	n139866028390888 -> n139866028390720
+	n139866028390944 [label=1]
+	n139866028390888 -> n139866028390944
+	n139866028390720 [label="*"]
+	n139866028695280 [label=ArrayAccess]
+	n139866028390720 -> n139866028695280
+	n139866028390776 [label=3]
+	n139866028390720 -> n139866028390776
+	n139866028390776 [label=3]
+	n139866028390944 [label=1]
+	n139866028425288 [label="+"]
+	n139866028392288 [label="*"]
+	n139866028425288 -> n139866028392288
+	n139866028425344 [label=2]
+	n139866028425288 -> n139866028425344
+	n139866028392288 [label="*"]
+	n139866028695280 [label=ArrayAccess]
+	n139866028392288 -> n139866028695280
+	n139866028392344 [label=3]
+	n139866028392288 -> n139866028392344
+	n139866028392344 [label=3]
+	n139866028425344 [label=2]
+	n139866028216560 [label="-"]
+	n139866028692480 [label=PropertyAccess]
+	n139866028216560 -> n139866028692480
+	n139866028216448 [label=PropertyAccess]
+	n139866028216560 -> n139866028216448
+	n139866028217064 [label="*"]
+	n139866028695280 [label=ArrayAccess]
+	n139866028217064 -> n139866028695280
+	n139866028217120 [label=3]
+	n139866028217064 -> n139866028217120
+	n139866028217120 [label=3]
+	n139866028217848 [label="*"]
+	n139866028695280 [label=ArrayAccess]
+	n139866028217848 -> n139866028695280
+	n139866028217904 [label=3]
+	n139866028217848 -> n139866028217904
+	n139866028217904 [label=3]
+	n139866028218240 [label="*"]
+	n139866028217400 [label=VectorAccess]
+	n139866028218240 -> n139866028217400
+	n139866028218184 [label=VectorAccess]
+	n139866028218240 -> n139866028218184
+	n139866028217400 [label=VectorAccess]
+	n139866028216560 [label="-"]
+	n139866028217400 -> n139866028216560
+	n139866028218184 [label=VectorAccess]
+	n139866028216560 [label="-"]
+	n139866028218184 -> n139866028216560
+	n139866028218744 [label="*"]
+	n139866028695280 [label=ArrayAccess]
+	n139866028218744 -> n139866028695280
+	n139866028218800 [label=3]
+	n139866028218744 -> n139866028218800
+	n139866028218800 [label=3]
+	n139866028218912 [label="+"]
+	n139866028218744 [label="*"]
+	n139866028218912 -> n139866028218744
+	n139866028218968 [label=1]
+	n139866028218912 -> n139866028218968
+	n139866028218968 [label=1]
+	n139866028219528 [label="*"]
+	n139866028695280 [label=ArrayAccess]
+	n139866028219528 -> n139866028695280
+	n139866028219584 [label=3]
+	n139866028219528 -> n139866028219584
+	n139866028219584 [label=3]
+	n139866028219696 [label="+"]
+	n139866028219528 [label="*"]
+	n139866028219696 -> n139866028219528
+	n139866028219752 [label=1]
+	n139866028219696 -> n139866028219752
+	n139866028219752 [label=1]
+	n139866028219920 [label="*"]
+	n139866028219080 [label=VectorAccess]
+	n139866028219920 -> n139866028219080
+	n139866028219864 [label=VectorAccess]
+	n139866028219920 -> n139866028219864
+	n139866028219080 [label=VectorAccess]
+	n139866028216560 [label="-"]
+	n139866028219080 -> n139866028216560
+	n139866028219864 [label=VectorAccess]
+	n139866028216560 [label="-"]
+	n139866028219864 -> n139866028216560
+	n139866028220032 [label="+"]
+	n139866028218240 [label="*"]
+	n139866028220032 -> n139866028218240
+	n139866028219920 [label="*"]
+	n139866028220032 -> n139866028219920
+	n139866028355768 [label="*"]
+	n139866028695280 [label=ArrayAccess]
+	n139866028355768 -> n139866028695280
+	n139866028355824 [label=3]
+	n139866028355768 -> n139866028355824
+	n139866028355824 [label=3]
+	n139866028355936 [label="+"]
+	n139866028355768 [label="*"]
+	n139866028355936 -> n139866028355768
+	n139866028355992 [label=2]
+	n139866028355936 -> n139866028355992
+	n139866028355992 [label=2]
+	n139866028356552 [label="*"]
+	n139866028695280 [label=ArrayAccess]
+	n139866028356552 -> n139866028695280
+	n139866028356608 [label=3]
+	n139866028356552 -> n139866028356608
+	n139866028356608 [label=3]
+	n139866028356720 [label="+"]
+	n139866028356552 [label="*"]
+	n139866028356720 -> n139866028356552
+	n139866028356776 [label=2]
+	n139866028356720 -> n139866028356776
+	n139866028356776 [label=2]
+	n139866028356944 [label="*"]
+	n139866028356104 [label=VectorAccess]
+	n139866028356944 -> n139866028356104
+	n139866028356888 [label=VectorAccess]
+	n139866028356944 -> n139866028356888
+	n139866028356104 [label=VectorAccess]
+	n139866028216560 [label="-"]
+	n139866028356104 -> n139866028216560
+	n139866028356888 [label=VectorAccess]
+	n139866028216560 [label="-"]
+	n139866028356888 -> n139866028216560
+	n139866028357056 [label="+"]
+	n139866028220032 [label="+"]
+	n139866028357056 -> n139866028220032
+	n139866028356944 [label="*"]
+	n139866028357056 -> n139866028356944
+	n139866028357168 [label="<"]
+	n139866028357056 [label="+"]
+	n139866028357168 -> n139866028357056
+	n139866028357224 [label=2.5]
+	n139866028357168 -> n139866028357224
+	n139866028357224 [label=2.5]
+	n139866028357336 [label=Filter]
+	n139866028357168 [label="<"]
+	n139866028357336 -> n139866028357168
+	n139866028357392 [label=Block]
+	n139866028357336 -> n139866028357392
+	n139866028357392 [label=Block]
+	n139866028359296 [label=Assign]
+	n139866028357392 -> n139866028359296
+	n139866028357728 [label="/"]
+	n139866028357784 [label=1.0]
+	n139866028357728 -> n139866028357784
+	n139866028357056 [label="+"]
+	n139866028357728 -> n139866028357056
+	n139866028357784 [label=1.0]
+	n139866028357896 [label="*"]
+	n139866028357728 [label="/"]
+	n139866028357896 -> n139866028357728
+	n139866028357728 [label="/"]
+	n139866028357896 -> n139866028357728
+	n139866028358008 [label="*"]
+	n139866028357896 [label="*"]
+	n139866028358008 -> n139866028357896
+	n139866028357728 [label="/"]
+	n139866028358008 -> n139866028357728
+	n139866028357672 [label=PropertyAccess]
+	n139866028640016 [label=force]
+	n139866028357672 -> n139866028640016
+	n139866028694720 [label="Iter(0)"]
+	n139866028357672 -> n139866028694720
+	n139866028388536 [label="*"]
+	n139866028357672 -> n139866028388536
+	n139866028390272 [label="+"]
+	n139866028357672 -> n139866028390272
+	n139866028391840 [label="+"]
+	n139866028357672 -> n139866028391840
+	n139866028640016 [label=force]
+	n139866028358344 [label="*"]
+	n139866028216560 [label="-"]
+	n139866028358344 -> n139866028216560
+	n139866028358400 [label=48.0]
+	n139866028358344 -> n139866028358400
+	n139866028358400 [label=48.0]
+	n139866028358512 [label="*"]
+	n139866028216560 [label="-"]
+	n139866028358512 -> n139866028216560
+	n139866027164896 [label="*"]
+	n139866028358512 -> n139866027164896
+	n139866027164896 [label="*"]
+	n139866028358400 [label=48.0]
+	n139866027164896 -> n139866028358400
+	n139866028358008 [label="*"]
+	n139866027164896 -> n139866028358008
+	n139866028358624 [label="-"]
+	n139866028358008 [label="*"]
+	n139866028358624 -> n139866028358008
+	n139866028358680 [label=0.5]
+	n139866028358624 -> n139866028358680
+	n139866028358680 [label=0.5]
+	n139866028358792 [label="*"]
+	n139866028216560 [label="-"]
+	n139866028358792 -> n139866028216560
+	n139866027165008 [label="*"]
+	n139866028358792 -> n139866027165008
+	n139866027165008 [label="*"]
+	n139866027164896 [label="*"]
+	n139866027165008 -> n139866027164896
+	n139866028358624 [label="-"]
+	n139866027165008 -> n139866028358624
+	n139866028358904 [label="*"]
+	n139866028216560 [label="-"]
+	n139866028358904 -> n139866028216560
+	n139866027165120 [label="*"]
+	n139866028358904 -> n139866027165120
+	n139866027165120 [label="*"]
+	n139866027165008 [label="*"]
+	n139866027165120 -> n139866027165008
+	n139866028357728 [label="/"]
+	n139866027165120 -> n139866028357728
+	n139866028359016 [label="*"]
+	n139866028216560 [label="-"]
+	n139866028359016 -> n139866028216560
+	n139866027165120 [label="*"]
+	n139866028359016 -> n139866027165120
+	n139866028359184 [label="+"]
+	n139866028357672 [label=PropertyAccess]
+	n139866028359184 -> n139866028357672
+	n139866028359016 [label="*"]
+	n139866028359184 -> n139866028359016
+	n139866028359296 [label=Assign]
+	n139866028389600 [label=VectorAccess]
+	n139866028359296 -> n139866028389600
+	n139866028388480 [label=VectorAccess]
+	n139866028359296 -> n139866028388480
+	n139866028391168 [label=VectorAccess]
+	n139866028359296 -> n139866028391168
+	n139866028390048 [label=VectorAccess]
+	n139866028359296 -> n139866028390048
+	n139866028425568 [label=VectorAccess]
+	n139866028359296 -> n139866028425568
+	n139866028391616 [label=VectorAccess]
+	n139866028359296 -> n139866028391616
+	n139866028389600 [label=VectorAccess]
+	n139866028357672 [label=PropertyAccess]
+	n139866028389600 -> n139866028357672
+	n139866028388480 [label=VectorAccess]
+	n139866028359184 [label="+"]
+	n139866028388480 -> n139866028359184
+	n139866028391168 [label=VectorAccess]
+	n139866028357672 [label=PropertyAccess]
+	n139866028391168 -> n139866028357672
+	n139866028390048 [label=VectorAccess]
+	n139866028359184 [label="+"]
+	n139866028390048 -> n139866028359184
+	n139866028425568 [label=VectorAccess]
+	n139866028357672 [label=PropertyAccess]
+	n139866028425568 -> n139866028357672
+	n139866028391616 [label=VectorAccess]
+	n139866028359184 [label="+"]
+	n139866028391616 -> n139866028359184
+	n139866028695168 [label=0]
+	n139866028695000 [label=ArrayAccess]
+	n139866028691920 [label=numneighs]
+	n139866028695000 -> n139866028691920
+	n139866028694720 [label="Iter(0)"]
+	n139866028695000 -> n139866028694720
+	n139866028691920 [label=numneighs]
+	n139866028694776 [label=0]
+	n139866028694832 [label=0]
+	n139866028694664 [label=ParticleFor]
+	n139866028693656 [label="Iter(2)"]
+	n139866028694664 -> n139866028693656
+	n139866028692592 [label=Block]
+	n139866028694664 -> n139866028692592
+	n139866028693544 [label=0]
+	n139866028694664 -> n139866028693544
+	n139866028693320 [label=0]
+	n139866028694664 -> n139866028693320
+	n139866028693656 [label="Iter(2)"]
+	n139866028692592 [label=Block]
+	n139866028426240 [label=Assign]
+	n139866028692592 -> n139866028426240
+	n139866027930848 [label=Assign]
+	n139866028692592 -> n139866027930848
+	n139866028357560 [label=PropertyAccess]
+	n139866028639960 [label=velocity]
+	n139866028357560 -> n139866028639960
+	n139866028693656 [label="Iter(2)"]
+	n139866028357560 -> n139866028693656
+	n139866028426744 [label="*"]
+	n139866028357560 -> n139866028426744
+	n139866028428480 [label="+"]
+	n139866028357560 -> n139866028428480
+	n139866027930064 [label="+"]
+	n139866028357560 -> n139866027930064
+	n139866028639960 [label=velocity]
+	n139866028426744 [label="*"]
+	n139866028693656 [label="Iter(2)"]
+	n139866028426744 -> n139866028693656
+	n139866028427136 [label=3]
+	n139866028426744 -> n139866028427136
+	n139866028427136 [label=3]
+	n139866028428480 [label="+"]
+	n139866028427976 [label="*"]
+	n139866028428480 -> n139866028427976
+	n139866028428536 [label=1]
+	n139866028428480 -> n139866028428536
+	n139866028427976 [label="*"]
+	n139866028693656 [label="Iter(2)"]
+	n139866028427976 -> n139866028693656
+	n139866028428368 [label=3]
+	n139866028427976 -> n139866028428368
+	n139866028428368 [label=3]
+	n139866028428536 [label=1]
+	n139866027930064 [label="+"]
+	n139866028429208 [label="*"]
+	n139866027930064 -> n139866028429208
+	n139866027930120 [label=2]
+	n139866027930064 -> n139866027930120
+	n139866028429208 [label="*"]
+	n139866028693656 [label="Iter(2)"]
+	n139866028429208 -> n139866028693656
+	n139866027929952 [label=3]
+	n139866028429208 -> n139866027929952
+	n139866027929952 [label=3]
+	n139866027930120 [label=2]
+	n139866028425624 [label=PropertyAccess]
+	n139866028640016 [label=force]
+	n139866028425624 -> n139866028640016
+	n139866028693656 [label="Iter(2)"]
+	n139866028425624 -> n139866028693656
+	n139866028426800 [label="*"]
+	n139866028425624 -> n139866028426800
+	n139866028428200 [label="+"]
+	n139866028425624 -> n139866028428200
+	n139866027929784 [label="+"]
+	n139866028425624 -> n139866027929784
+	n139866028426800 [label="*"]
+	n139866028693656 [label="Iter(2)"]
+	n139866028426800 -> n139866028693656
+	n139866028426856 [label=3]
+	n139866028426800 -> n139866028426856
+	n139866028426856 [label=3]
+	n139866028428200 [label="+"]
+	n139866028428032 [label="*"]
+	n139866028428200 -> n139866028428032
+	n139866028428256 [label=1]
+	n139866028428200 -> n139866028428256
+	n139866028428032 [label="*"]
+	n139866028693656 [label="Iter(2)"]
+	n139866028428032 -> n139866028693656
+	n139866028428088 [label=3]
+	n139866028428032 -> n139866028428088
+	n139866028428088 [label=3]
+	n139866028428256 [label=1]
+	n139866027929784 [label="+"]
+	n139866028429264 [label="*"]
+	n139866027929784 -> n139866028429264
+	n139866027929840 [label=2]
+	n139866027929784 -> n139866027929840
+	n139866028429264 [label="*"]
+	n139866028693656 [label="Iter(2)"]
+	n139866028429264 -> n139866028693656
+	n139866027929672 [label=3]
+	n139866028429264 -> n139866027929672
+	n139866027929672 [label=3]
+	n139866027929840 [label=2]
+	n139866028425736 [label="*"]
+	n139866028425792 [label=0.005]
+	n139866028425736 -> n139866028425792
+	n139866028425624 [label=PropertyAccess]
+	n139866028425736 -> n139866028425624
+	n139866028425792 [label=0.005]
+	n139866028425904 [label=PropertyAccess]
+	n139866028639848 [label=mass]
+	n139866028425904 -> n139866028639848
+	n139866028693656 [label="Iter(2)"]
+	n139866028425904 -> n139866028693656
+	n139866028639848 [label=mass]
+	n139866028426016 [label="/"]
+	n139866028425736 [label="*"]
+	n139866028426016 -> n139866028425736
+	n139866028425904 [label=PropertyAccess]
+	n139866028426016 -> n139866028425904
+	n139866028426128 [label="+"]
+	n139866028357560 [label=PropertyAccess]
+	n139866028426128 -> n139866028357560
+	n139866028426016 [label="/"]
+	n139866028426128 -> n139866028426016
+	n139866028426352 [label="*"]
+	n139866028693656 [label="Iter(2)"]
+	n139866028426352 -> n139866028693656
+	n139866028426408 [label=3]
+	n139866028426352 -> n139866028426408
+	n139866028426408 [label=3]
+	n139866028427584 [label="*"]
+	n139866028693656 [label="Iter(2)"]
+	n139866028427584 -> n139866028693656
+	n139866028427640 [label=3]
+	n139866028427584 -> n139866028427640
+	n139866028427640 [label=3]
+	n139866028427752 [label="+"]
+	n139866028427584 [label="*"]
+	n139866028427752 -> n139866028427584
+	n139866028427808 [label=1]
+	n139866028427752 -> n139866028427808
+	n139866028427808 [label=1]
+	n139866028428816 [label="*"]
+	n139866028693656 [label="Iter(2)"]
+	n139866028428816 -> n139866028693656
+	n139866028428872 [label=3]
+	n139866028428816 -> n139866028428872
+	n139866028428872 [label=3]
+	n139866028428984 [label="+"]
+	n139866028428816 [label="*"]
+	n139866028428984 -> n139866028428816
+	n139866028429040 [label=2]
+	n139866028428984 -> n139866028429040
+	n139866028429040 [label=2]
+	n139866028426240 [label=Assign]
+	n139866028427472 [label=VectorAccess]
+	n139866028426240 -> n139866028427472
+	n139866028426688 [label=VectorAccess]
+	n139866028426240 -> n139866028426688
+	n139866028428704 [label=VectorAccess]
+	n139866028426240 -> n139866028428704
+	n139866028427920 [label=VectorAccess]
+	n139866028426240 -> n139866028427920
+	n139866027930288 [label=VectorAccess]
+	n139866028426240 -> n139866027930288
+	n139866028429152 [label=VectorAccess]
+	n139866028426240 -> n139866028429152
+	n139866028427472 [label=VectorAccess]
+	n139866028357560 [label=PropertyAccess]
+	n139866028427472 -> n139866028357560
+	n139866028426688 [label=VectorAccess]
+	n139866028426128 [label="+"]
+	n139866028426688 -> n139866028426128
+	n139866028428704 [label=VectorAccess]
+	n139866028357560 [label=PropertyAccess]
+	n139866028428704 -> n139866028357560
+	n139866028427920 [label=VectorAccess]
+	n139866028426128 [label="+"]
+	n139866028427920 -> n139866028426128
+	n139866027930288 [label=VectorAccess]
+	n139866028357560 [label=PropertyAccess]
+	n139866027930288 -> n139866028357560
+	n139866028429152 [label=VectorAccess]
+	n139866028426128 [label="+"]
+	n139866028429152 -> n139866028426128
+	n139866027930344 [label=PropertyAccess]
+	n139866028639904 [label=position]
+	n139866027930344 -> n139866028639904
+	n139866028693656 [label="Iter(2)"]
+	n139866027930344 -> n139866028693656
+	n139866027931688 [label="*"]
+	n139866027930344 -> n139866027931688
+	n139866027933088 [label="+"]
+	n139866027930344 -> n139866027933088
+	n139866027967152 [label="+"]
+	n139866027930344 -> n139866027967152
+	n139866027931688 [label="*"]
+	n139866028693656 [label="Iter(2)"]
+	n139866027931688 -> n139866028693656
+	n139866027931744 [label=3]
+	n139866027931688 -> n139866027931744
+	n139866027931744 [label=3]
+	n139866027933088 [label="+"]
+	n139866027932920 [label="*"]
+	n139866027933088 -> n139866027932920
+	n139866027933144 [label=1]
+	n139866027933088 -> n139866027933144
+	n139866027932920 [label="*"]
+	n139866028693656 [label="Iter(2)"]
+	n139866027932920 -> n139866028693656
+	n139866027932976 [label=3]
+	n139866027932920 -> n139866027932976
+	n139866027932976 [label=3]
+	n139866027933144 [label=1]
+	n139866027967152 [label="+"]
+	n139866027966984 [label="*"]
+	n139866027967152 -> n139866027966984
+	n139866027967208 [label=2]
+	n139866027967152 -> n139866027967208
+	n139866027966984 [label="*"]
+	n139866028693656 [label="Iter(2)"]
+	n139866027966984 -> n139866028693656
+	n139866027967040 [label=3]
+	n139866027966984 -> n139866027967040
+	n139866027967040 [label=3]
+	n139866027967208 [label=2]
+	n139866027930456 [label=PropertyAccess]
+	n139866028639960 [label=velocity]
+	n139866027930456 -> n139866028639960
+	n139866028693656 [label="Iter(2)"]
+	n139866027930456 -> n139866028693656
+	n139866027931352 [label="*"]
+	n139866027930456 -> n139866027931352
+	n139866027932752 [label="+"]
+	n139866027930456 -> n139866027932752
+	n139866027966816 [label="+"]
+	n139866027930456 -> n139866027966816
+	n139866027931352 [label="*"]
+	n139866028693656 [label="Iter(2)"]
+	n139866027931352 -> n139866028693656
+	n139866027931408 [label=3]
+	n139866027931352 -> n139866027931408
+	n139866027931408 [label=3]
+	n139866027932752 [label="+"]
+	n139866027932584 [label="*"]
+	n139866027932752 -> n139866027932584
+	n139866027932808 [label=1]
+	n139866027932752 -> n139866027932808
+	n139866027932584 [label="*"]
+	n139866028693656 [label="Iter(2)"]
+	n139866027932584 -> n139866028693656
+	n139866027932640 [label=3]
+	n139866027932584 -> n139866027932640
+	n139866027932640 [label=3]
+	n139866027932808 [label=1]
+	n139866027966816 [label="+"]
+	n139866027966648 [label="*"]
+	n139866027966816 -> n139866027966648
+	n139866027966872 [label=2]
+	n139866027966816 -> n139866027966872
+	n139866027966648 [label="*"]
+	n139866028693656 [label="Iter(2)"]
+	n139866027966648 -> n139866028693656
+	n139866027966704 [label=3]
+	n139866027966648 -> n139866027966704
+	n139866027966704 [label=3]
+	n139866027966872 [label=2]
+	n139866027930568 [label="*"]
+	n139866027930624 [label=0.005]
+	n139866027930568 -> n139866027930624
+	n139866027930456 [label=PropertyAccess]
+	n139866027930568 -> n139866027930456
+	n139866027930624 [label=0.005]
+	n139866027930736 [label="+"]
+	n139866027930344 [label=PropertyAccess]
+	n139866027930736 -> n139866027930344
+	n139866027930568 [label="*"]
+	n139866027930736 -> n139866027930568
+	n139866027930960 [label="*"]
+	n139866028693656 [label="Iter(2)"]
+	n139866027930960 -> n139866028693656
+	n139866027931016 [label=3]
+	n139866027930960 -> n139866027931016
+	n139866027931016 [label=3]
+	n139866027932192 [label="*"]
+	n139866028693656 [label="Iter(2)"]
+	n139866027932192 -> n139866028693656
+	n139866027932248 [label=3]
+	n139866027932192 -> n139866027932248
+	n139866027932248 [label=3]
+	n139866027932360 [label="+"]
+	n139866027932192 [label="*"]
+	n139866027932360 -> n139866027932192
+	n139866027932416 [label=1]
+	n139866027932360 -> n139866027932416
+	n139866027932416 [label=1]
+	n139866027933424 [label="*"]
+	n139866028693656 [label="Iter(2)"]
+	n139866027933424 -> n139866028693656
+	n139866027933480 [label=3]
+	n139866027933424 -> n139866027933480
+	n139866027933480 [label=3]
+	n139866027933592 [label="+"]
+	n139866027933424 [label="*"]
+	n139866027933592 -> n139866027933424
+	n139866027933648 [label=2]
+	n139866027933592 -> n139866027933648
+	n139866027933648 [label=2]
+	n139866027930848 [label=Assign]
+	n139866027932080 [label=VectorAccess]
+	n139866027930848 -> n139866027932080
+	n139866027931296 [label=VectorAccess]
+	n139866027930848 -> n139866027931296
+	n139866027933312 [label=VectorAccess]
+	n139866027930848 -> n139866027933312
+	n139866027932528 [label=VectorAccess]
+	n139866027930848 -> n139866027932528
+	n139866027967376 [label=VectorAccess]
+	n139866027930848 -> n139866027967376
+	n139866027966592 [label=VectorAccess]
+	n139866027930848 -> n139866027966592
+	n139866027932080 [label=VectorAccess]
+	n139866027930344 [label=PropertyAccess]
+	n139866027932080 -> n139866027930344
+	n139866027931296 [label=VectorAccess]
+	n139866027930736 [label="+"]
+	n139866027931296 -> n139866027930736
+	n139866027933312 [label=VectorAccess]
+	n139866027930344 [label=PropertyAccess]
+	n139866027933312 -> n139866027930344
+	n139866027932528 [label=VectorAccess]
+	n139866027930736 [label="+"]
+	n139866027932528 -> n139866027930736
+	n139866027967376 [label=VectorAccess]
+	n139866027930344 [label=PropertyAccess]
+	n139866027967376 -> n139866027930344
+	n139866027966592 [label=VectorAccess]
+	n139866027930736 [label="+"]
+	n139866027966592 -> n139866027930736
+	n139866028693544 [label=0]
+	n139866028693320 [label=0]
+}
diff --git a/examples/kernels.pdf b/examples/kernels.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..54c0d45a7732838899d0d081e18c622211358f87
Binary files /dev/null and b/examples/kernels.pdf differ
diff --git a/lift.py b/examples/lift.py
similarity index 100%
rename from lift.py
rename to examples/lift.py
diff --git a/particle.py b/examples/lj_embedded.py
similarity index 65%
rename from particle.py
rename to examples/lj_embedded.py
index 8680d2b46ecfbfa920437be69b4b3e12acd5123a..160c6081bacc9e0a3cd6e8561ebff9ea2eb7bf15 100644
--- a/particle.py
+++ b/examples/lj_embedded.py
@@ -1,5 +1,5 @@
 import pairs
-from ir.layouts import Layout_SoA
+
 
 dt = 0.005
 cutoff_radius = 2.5
@@ -13,21 +13,6 @@ mass = psim.add_real_property('mass', 1.0)
 position = psim.add_vector_property('position')
 velocity = psim.add_vector_property('velocity')
 force = psim.add_vector_property('force', vol=True)
-
-#def lj(i, j):
-#    sr2 = 1.0 / rsq(i, j)
-#    sr6 = sr2 * sr2 * sr2 * sigma6
-#    force[i] += delta(i, j) * 48.0 * sr6 * (sr6 - 0.5) * sr2 * epsilon
-
-#def euler(i):
-#    velocity[i] += dt * force[i] / mass[i]
-#    position[i] += dt * velocity[i]
-
-# psim.compute(lj)
-# psim.compute(euler)
-
-#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.create_neighbor_lists()
diff --git a/examples/lj_func.py b/examples/lj_func.py
new file mode 100644
index 0000000000000000000000000000000000000000..5b2a7d4862ed49bd27b182f8192d718e6c35f9d6
--- /dev/null
+++ b/examples/lj_func.py
@@ -0,0 +1,45 @@
+import pairs
+
+
+def delta(i, j):
+    return position[i] - position[j]
+
+
+def rsq(i, j):
+    dp = delta(i, j)
+    return dp.x() * dp.x() + dp.y() * dp.y() + dp.z() * dp.z()
+
+
+def lj(i, j):
+    sr2 = 1.0 / rsq
+    sr6 = sr2 * sr2 * sr2 * sigma6
+    #f = 48.0 * sr6 * (sr6 - 0.5) * sr2 * epsilon
+    #force[i] += delta * f
+    force[i] += delta * 48.0 * sr6 * (sr6 - 0.5) * sr2 * epsilon
+
+
+def euler(i):
+    velocity[i] += dt * force[i] / mass[i]
+    position[i] += dt * velocity[i]
+
+
+dt = 0.005
+cutoff_radius = 2.5
+skin = 0.3
+sigma = 1.0
+epsilon = 1.0
+sigma6 = sigma ** 6
+
+psim = pairs.simulation("lj_ns")
+psim.add_real_property('mass', 1.0)
+psim.add_vector_property('position')
+psim.add_vector_property('velocity')
+psim.add_vector_property('force', vol=True)
+psim.from_file("data/minimd_setup_4x4x4.input", ['mass', 'position', 'velocity'])
+psim.create_cell_lists(2.8, 2.8)
+psim.create_neighbor_lists()
+psim.periodic(2.8)
+psim.vtk_output("output/test")
+psim.compute(psim, lj, cutoff_radius, 'position', {'sigma6': sigma6, 'epsilon': epsilon})
+psim.compute(psim, euler, symbols={'dt': dt})
+psim.generate()
diff --git a/examples/lj_ns.cpp b/examples/lj_ns.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..bc536ce388b0def6feb02d9711eff97ae4137561
--- /dev/null
+++ b/examples/lj_ns.cpp
@@ -0,0 +1,825 @@
+#include <math.h>
+#include <stdbool.h>
+#include <stdio.h>
+#include <stdlib.h>
+//---
+#include "runtime/pairs.hpp"
+#include "runtime/read_from_file.hpp"
+#include "runtime/vtk.hpp"
+
+using namespace pairs;
+
+int main() {
+    PairsSim *ps = new PairsSim();
+    int particle_capacity = 10000;
+    int nlocal = 0;
+    int nghost = 0;
+    double grid0_d0_min = 0;
+    double grid0_d0_max = 0;
+    double grid0_d1_min = 0;
+    double grid0_d1_max = 0;
+    double grid0_d2_min = 0;
+    double grid0_d2_max = 0;
+    int nstencil = 0;
+    int ncells = 1;
+    int ncells_capacity = 100;
+    int cell_capacity = 20;
+    int neighborlist_capacity = 32;
+    int npbc = 0;
+    int pbc_capacity = 100;
+    int resize = 0;
+    double grid_buffer[6];
+    int dim_cells[3];
+    int *cell_particles = (int *) malloc((sizeof(int) * (ncells_capacity * cell_capacity)));
+    int *cell_sizes = (int *) malloc((sizeof(int) * ncells_capacity));
+    int *stencil = (int *) malloc((sizeof(int) * 27));
+    int *particle_cell = (int *) malloc((sizeof(int) * particle_capacity));
+    int *neighborlists = (int *) malloc((sizeof(int) * (particle_capacity * neighborlist_capacity)));
+    int *numneighs = (int *) malloc((sizeof(int) * particle_capacity));
+    int *pbc_map = (int *) malloc((sizeof(int) * pbc_capacity));
+    int *pbc_mult = (int *) malloc((sizeof(int) * (pbc_capacity * 3)));
+    double *mass = (double *) malloc((sizeof(double) * (particle_capacity + pbc_capacity)));
+    ps->addProperty(Property(0, "mass", mass, Prop_Float));
+    double *position = (double *) malloc((sizeof(double) * ((particle_capacity + pbc_capacity) * 3)));
+    ps->addProperty(Property(1, "position", position, Prop_Vector, AoS, (particle_capacity + pbc_capacity), 3));
+    double *velocity = (double *) malloc((sizeof(double) * ((particle_capacity + pbc_capacity) * 3)));
+    ps->addProperty(Property(2, "velocity", velocity, Prop_Vector, AoS, (particle_capacity + pbc_capacity), 3));
+    double *force = (double *) malloc((sizeof(double) * ((particle_capacity + pbc_capacity) * 3)));
+    ps->addProperty(Property(3, "force", force, Prop_Vector, AoS, (particle_capacity + pbc_capacity), 3));
+    const int prop_list_0[] = {0, 1, 2};
+    nlocal = pairs::read_particle_data(ps, "data/minimd_setup_4x4x4.input", grid_buffer, prop_list_0, 3);
+    const double a0 = grid_buffer[0];
+    grid0_d0_min = a0;
+    const double a1 = grid_buffer[1];
+    grid0_d0_max = a1;
+    const double a2 = grid_buffer[2];
+    grid0_d1_min = a2;
+    const double a3 = grid_buffer[3];
+    grid0_d1_max = a3;
+    const double a4 = grid_buffer[4];
+    grid0_d2_min = a4;
+    const double a5 = grid_buffer[5];
+    grid0_d2_max = a5;
+    fprintf(stdout, "CellListsStencilBuild\n");
+    fflush(stdout);
+    const double e462 = grid0_d0_max - grid0_d0_min;
+    const double e463 = e462 / 2.8;
+    const int e464 = ceil(e463) + 2;
+    dim_cells[0] = e464;
+    const double e466 = grid0_d1_max - grid0_d1_min;
+    const double e467 = e466 / 2.8;
+    const int e468 = ceil(e467) + 2;
+    dim_cells[1] = e468;
+    const int a54 = dim_cells[0];
+    const int a56 = dim_cells[1];
+    const int e469 = a54 * a56;
+    const double e470 = grid0_d2_max - grid0_d2_min;
+    const double e471 = e470 / 2.8;
+    const int e472 = ceil(e471) + 2;
+    dim_cells[2] = e472;
+    const int a58 = dim_cells[2];
+    const int e473 = e469 * a58;
+    ncells = e473;
+    resize = 1;
+    while((resize > 0)) {
+        resize = 0;
+        const bool e475 = ncells >= ncells_capacity;
+        if(e475) {
+            resize = ncells;
+        }
+        const bool e476 = resize > 0;
+        if(e476) {
+            fprintf(stdout, "Resize ncells_capacity\n");
+            fflush(stdout);
+            const int e477 = resize * 2;
+            ncells_capacity = e477;
+            cell_particles = (int *) realloc(cell_particles, (sizeof(int) * (ncells_capacity * cell_capacity)));
+            cell_sizes = (int *) realloc(cell_sizes, (sizeof(int) * ncells_capacity));
+        }
+    }
+    nstencil = 0;
+    for(int i15 = -1; i15 < 2; i15++) {
+        const int a59 = dim_cells[0];
+        const int e481 = i15 * a59;
+        for(int i16 = -1; i16 < 2; i16++) {
+            const int e482 = e481 + i16;
+            const int a60 = dim_cells[1];
+            const int e483 = e482 * a60;
+            for(int i17 = -1; i17 < 2; i17++) {
+                const int e484 = e483 + i17;
+                stencil[nstencil] = e484;
+                const int e485 = nstencil + 1;
+                nstencil = e485;
+            }
+        }
+    }
+    const int e486 = nlocal + npbc;
+    pairs::vtk_write_data(ps, "output/test_local", 0, nlocal, 0);
+    pairs::vtk_write_data(ps, "output/test_pbc", nlocal, e486, 0);
+    for(int i14 = 0; i14 < 101; i14++) {
+        const int e452 = i14 % 20;
+        const bool e453 = e452 == 0;
+        if(e453) {
+            fprintf(stdout, "EnforcePBC\n");
+            fflush(stdout);
+            const double e115 = grid0_d0_max - grid0_d0_min;
+            const double e122 = grid0_d0_max - grid0_d0_min;
+            const double e129 = grid0_d1_max - grid0_d1_min;
+            const double e136 = grid0_d1_max - grid0_d1_min;
+            const double e143 = grid0_d2_max - grid0_d2_min;
+            const double e150 = grid0_d2_max - grid0_d2_min;
+            for(int i3 = 0; i3 < nlocal; i3++) {
+                const int e110 = i3 * 3;
+                const double p8_0 = position[e110];
+                const bool e112 = p8_0 < grid0_d0_min;
+                if(e112) {
+                    const int e113 = i3 * 3;
+                    const double p9_0 = position[e113];
+                    const double e116 = p9_0 + e115;
+                    position[e113] = e116;
+                }
+                const int e117 = i3 * 3;
+                const double p10_0 = position[e117];
+                const bool e119 = p10_0 > grid0_d0_max;
+                if(e119) {
+                    const int e120 = i3 * 3;
+                    const double p11_0 = position[e120];
+                    const double e123 = p11_0 - e122;
+                    position[e120] = e123;
+                }
+                const int e124 = i3 * 3;
+                const int e125 = e124 + 1;
+                const double p12_1 = position[e125];
+                const bool e126 = p12_1 < grid0_d1_min;
+                if(e126) {
+                    const int e127 = i3 * 3;
+                    const int e128 = e127 + 1;
+                    const double p13_1 = position[e128];
+                    const double e130 = p13_1 + e129;
+                    position[e128] = e130;
+                }
+                const int e131 = i3 * 3;
+                const int e132 = e131 + 1;
+                const double p14_1 = position[e132];
+                const bool e133 = p14_1 > grid0_d1_max;
+                if(e133) {
+                    const int e134 = i3 * 3;
+                    const int e135 = e134 + 1;
+                    const double p15_1 = position[e135];
+                    const double e137 = p15_1 - e136;
+                    position[e135] = e137;
+                }
+                const int e138 = i3 * 3;
+                const int e139 = e138 + 2;
+                const double p16_2 = position[e139];
+                const bool e140 = p16_2 < grid0_d2_min;
+                if(e140) {
+                    const int e141 = i3 * 3;
+                    const int e142 = e141 + 2;
+                    const double p17_2 = position[e142];
+                    const double e144 = p17_2 + e143;
+                    position[e142] = e144;
+                }
+                const int e145 = i3 * 3;
+                const int e146 = e145 + 2;
+                const double p18_2 = position[e146];
+                const bool e147 = p18_2 > grid0_d2_max;
+                if(e147) {
+                    const int e148 = i3 * 3;
+                    const int e149 = e148 + 2;
+                    const double p19_2 = position[e149];
+                    const double e151 = p19_2 - e150;
+                    position[e149] = e151;
+                }
+            }
+        }
+        const int e454 = i14 % 20;
+        const bool e455 = e454 == 0;
+        if(e455) {
+            fprintf(stdout, "SetupPBC\n");
+            fflush(stdout);
+            resize = 1;
+            while((resize > 0)) {
+                resize = 0;
+                npbc = 0;
+                const int e153 = nlocal + npbc;
+                const double e157 = grid0_d0_min + 2.8;
+                const double e168 = grid0_d0_max - grid0_d0_min;
+                const double e185 = grid0_d0_max - 2.8;
+                const double e196 = grid0_d0_max - grid0_d0_min;
+                for(int i4 = 0; i4 < e153; i4++) {
+                    const int e154 = nlocal + npbc;
+                    const int e155 = i4 * 3;
+                    const double p20_0 = position[e155];
+                    const bool e158 = p20_0 < e157;
+                    if(e158) {
+                        const bool e159 = npbc >= pbc_capacity;
+                        if(e159) {
+                            const bool e160 = resize > npbc;
+                            resize = (e160) ? ((resize + 1)) : (npbc);
+                        } else {
+                            pbc_map[npbc] = i4;
+                            const int e162 = npbc * 3;
+                            pbc_mult[e162] = 1;
+                            const int e164 = e154 * 3;
+                            const double p21_0 = position[e164];
+                            const int e166 = i4 * 3;
+                            const double p22_0 = position[e166];
+                            const double e169 = p22_0 + e168;
+                            position[e164] = e169;
+                            const int e170 = npbc * 3;
+                            const int e171 = e170 + 1;
+                            pbc_mult[e171] = 0;
+                            const int e172 = e154 * 3;
+                            const int e173 = e172 + 1;
+                            const double p23_1 = position[e173];
+                            const int e174 = i4 * 3;
+                            const int e175 = e174 + 1;
+                            const double p24_1 = position[e175];
+                            position[e173] = p24_1;
+                            const int e176 = npbc * 3;
+                            const int e177 = e176 + 2;
+                            pbc_mult[e177] = 0;
+                            const int e178 = e154 * 3;
+                            const int e179 = e178 + 2;
+                            const double p25_2 = position[e179];
+                            const int e180 = i4 * 3;
+                            const int e181 = e180 + 2;
+                            const double p26_2 = position[e181];
+                            position[e179] = p26_2;
+                            const int e182 = npbc + 1;
+                            npbc = e182;
+                        }
+                    }
+                    const int e183 = i4 * 3;
+                    const double p27_0 = position[e183];
+                    const bool e186 = p27_0 > e185;
+                    if(e186) {
+                        const bool e187 = npbc >= pbc_capacity;
+                        if(e187) {
+                            const bool e188 = resize > npbc;
+                            resize = (e188) ? ((resize + 1)) : (npbc);
+                        } else {
+                            pbc_map[npbc] = i4;
+                            const int e190 = npbc * 3;
+                            pbc_mult[e190] = -1;
+                            const int e192 = e154 * 3;
+                            const double p28_0 = position[e192];
+                            const int e194 = i4 * 3;
+                            const double p29_0 = position[e194];
+                            const double e197 = p29_0 - e196;
+                            position[e192] = e197;
+                            const int e198 = npbc * 3;
+                            const int e199 = e198 + 1;
+                            pbc_mult[e199] = 0;
+                            const int e200 = e154 * 3;
+                            const int e201 = e200 + 1;
+                            const double p30_1 = position[e201];
+                            const int e202 = i4 * 3;
+                            const int e203 = e202 + 1;
+                            const double p31_1 = position[e203];
+                            position[e201] = p31_1;
+                            const int e204 = npbc * 3;
+                            const int e205 = e204 + 2;
+                            pbc_mult[e205] = 0;
+                            const int e206 = e154 * 3;
+                            const int e207 = e206 + 2;
+                            const double p32_2 = position[e207];
+                            const int e208 = i4 * 3;
+                            const int e209 = e208 + 2;
+                            const double p33_2 = position[e209];
+                            position[e207] = p33_2;
+                            const int e210 = npbc + 1;
+                            npbc = e210;
+                        }
+                    }
+                }
+                const int e211 = nlocal + npbc;
+                const double e215 = grid0_d1_min + 2.8;
+                const double e226 = grid0_d1_max - grid0_d1_min;
+                const double e243 = grid0_d1_max - 2.8;
+                const double e254 = grid0_d1_max - grid0_d1_min;
+                for(int i5 = 0; i5 < e211; i5++) {
+                    const int e212 = nlocal + npbc;
+                    const int e213 = i5 * 3;
+                    const int e214 = e213 + 1;
+                    const double p34_1 = position[e214];
+                    const bool e216 = p34_1 < e215;
+                    if(e216) {
+                        const bool e217 = npbc >= pbc_capacity;
+                        if(e217) {
+                            const bool e218 = resize > npbc;
+                            resize = (e218) ? ((resize + 1)) : (npbc);
+                        } else {
+                            pbc_map[npbc] = i5;
+                            const int e220 = npbc * 3;
+                            const int e221 = e220 + 1;
+                            pbc_mult[e221] = 1;
+                            const int e222 = e212 * 3;
+                            const int e223 = e222 + 1;
+                            const double p35_1 = position[e223];
+                            const int e224 = i5 * 3;
+                            const int e225 = e224 + 1;
+                            const double p36_1 = position[e225];
+                            const double e227 = p36_1 + e226;
+                            position[e223] = e227;
+                            const int e228 = npbc * 3;
+                            pbc_mult[e228] = 0;
+                            const int e230 = e212 * 3;
+                            const double p37_0 = position[e230];
+                            const int e232 = i5 * 3;
+                            const double p38_0 = position[e232];
+                            position[e230] = p38_0;
+                            const int e234 = npbc * 3;
+                            const int e235 = e234 + 2;
+                            pbc_mult[e235] = 0;
+                            const int e236 = e212 * 3;
+                            const int e237 = e236 + 2;
+                            const double p39_2 = position[e237];
+                            const int e238 = i5 * 3;
+                            const int e239 = e238 + 2;
+                            const double p40_2 = position[e239];
+                            position[e237] = p40_2;
+                            const int e240 = npbc + 1;
+                            npbc = e240;
+                        }
+                    }
+                    const int e241 = i5 * 3;
+                    const int e242 = e241 + 1;
+                    const double p41_1 = position[e242];
+                    const bool e244 = p41_1 > e243;
+                    if(e244) {
+                        const bool e245 = npbc >= pbc_capacity;
+                        if(e245) {
+                            const bool e246 = resize > npbc;
+                            resize = (e246) ? ((resize + 1)) : (npbc);
+                        } else {
+                            pbc_map[npbc] = i5;
+                            const int e248 = npbc * 3;
+                            const int e249 = e248 + 1;
+                            pbc_mult[e249] = -1;
+                            const int e250 = e212 * 3;
+                            const int e251 = e250 + 1;
+                            const double p42_1 = position[e251];
+                            const int e252 = i5 * 3;
+                            const int e253 = e252 + 1;
+                            const double p43_1 = position[e253];
+                            const double e255 = p43_1 - e254;
+                            position[e251] = e255;
+                            const int e256 = npbc * 3;
+                            pbc_mult[e256] = 0;
+                            const int e258 = e212 * 3;
+                            const double p44_0 = position[e258];
+                            const int e260 = i5 * 3;
+                            const double p45_0 = position[e260];
+                            position[e258] = p45_0;
+                            const int e262 = npbc * 3;
+                            const int e263 = e262 + 2;
+                            pbc_mult[e263] = 0;
+                            const int e264 = e212 * 3;
+                            const int e265 = e264 + 2;
+                            const double p46_2 = position[e265];
+                            const int e266 = i5 * 3;
+                            const int e267 = e266 + 2;
+                            const double p47_2 = position[e267];
+                            position[e265] = p47_2;
+                            const int e268 = npbc + 1;
+                            npbc = e268;
+                        }
+                    }
+                }
+                const int e269 = nlocal + npbc;
+                const double e273 = grid0_d2_min + 2.8;
+                const double e284 = grid0_d2_max - grid0_d2_min;
+                const double e301 = grid0_d2_max - 2.8;
+                const double e312 = grid0_d2_max - grid0_d2_min;
+                for(int i6 = 0; i6 < e269; i6++) {
+                    const int e270 = nlocal + npbc;
+                    const int e271 = i6 * 3;
+                    const int e272 = e271 + 2;
+                    const double p48_2 = position[e272];
+                    const bool e274 = p48_2 < e273;
+                    if(e274) {
+                        const bool e275 = npbc >= pbc_capacity;
+                        if(e275) {
+                            const bool e276 = resize > npbc;
+                            resize = (e276) ? ((resize + 1)) : (npbc);
+                        } else {
+                            pbc_map[npbc] = i6;
+                            const int e278 = npbc * 3;
+                            const int e279 = e278 + 2;
+                            pbc_mult[e279] = 1;
+                            const int e280 = e270 * 3;
+                            const int e281 = e280 + 2;
+                            const double p49_2 = position[e281];
+                            const int e282 = i6 * 3;
+                            const int e283 = e282 + 2;
+                            const double p50_2 = position[e283];
+                            const double e285 = p50_2 + e284;
+                            position[e281] = e285;
+                            const int e286 = npbc * 3;
+                            pbc_mult[e286] = 0;
+                            const int e288 = e270 * 3;
+                            const double p51_0 = position[e288];
+                            const int e290 = i6 * 3;
+                            const double p52_0 = position[e290];
+                            position[e288] = p52_0;
+                            const int e292 = npbc * 3;
+                            const int e293 = e292 + 1;
+                            pbc_mult[e293] = 0;
+                            const int e294 = e270 * 3;
+                            const int e295 = e294 + 1;
+                            const double p53_1 = position[e295];
+                            const int e296 = i6 * 3;
+                            const int e297 = e296 + 1;
+                            const double p54_1 = position[e297];
+                            position[e295] = p54_1;
+                            const int e298 = npbc + 1;
+                            npbc = e298;
+                        }
+                    }
+                    const int e299 = i6 * 3;
+                    const int e300 = e299 + 2;
+                    const double p55_2 = position[e300];
+                    const bool e302 = p55_2 > e301;
+                    if(e302) {
+                        const bool e303 = npbc >= pbc_capacity;
+                        if(e303) {
+                            const bool e304 = resize > npbc;
+                            resize = (e304) ? ((resize + 1)) : (npbc);
+                        } else {
+                            pbc_map[npbc] = i6;
+                            const int e306 = npbc * 3;
+                            const int e307 = e306 + 2;
+                            pbc_mult[e307] = -1;
+                            const int e308 = e270 * 3;
+                            const int e309 = e308 + 2;
+                            const double p56_2 = position[e309];
+                            const int e310 = i6 * 3;
+                            const int e311 = e310 + 2;
+                            const double p57_2 = position[e311];
+                            const double e313 = p57_2 - e312;
+                            position[e309] = e313;
+                            const int e314 = npbc * 3;
+                            pbc_mult[e314] = 0;
+                            const int e316 = e270 * 3;
+                            const double p58_0 = position[e316];
+                            const int e318 = i6 * 3;
+                            const double p59_0 = position[e318];
+                            position[e316] = p59_0;
+                            const int e320 = npbc * 3;
+                            const int e321 = e320 + 1;
+                            pbc_mult[e321] = 0;
+                            const int e322 = e270 * 3;
+                            const int e323 = e322 + 1;
+                            const double p60_1 = position[e323];
+                            const int e324 = i6 * 3;
+                            const int e325 = e324 + 1;
+                            const double p61_1 = position[e325];
+                            position[e323] = p61_1;
+                            const int e326 = npbc + 1;
+                            npbc = e326;
+                        }
+                    }
+                }
+                const bool e327 = resize > 0;
+                if(e327) {
+                    fprintf(stdout, "Resize pbc_capacity\n");
+                    fflush(stdout);
+                    const int e328 = resize * 2;
+                    pbc_capacity = e328;
+                    pbc_map = (int *) realloc(pbc_map, (sizeof(int) * pbc_capacity));
+                    pbc_mult = (int *) realloc(pbc_mult, (sizeof(int) * (pbc_capacity * 3)));
+                    mass = (double *) realloc(mass, (sizeof(double) * (particle_capacity + pbc_capacity)));
+                    ps->updateProperty(0, mass);
+                    position = (double *) realloc(position, (sizeof(double) * ((particle_capacity + pbc_capacity) * 3)));
+                    ps->updateProperty(1, position, (particle_capacity + pbc_capacity), 3);
+                    velocity = (double *) realloc(velocity, (sizeof(double) * ((particle_capacity + pbc_capacity) * 3)));
+                    ps->updateProperty(2, velocity, (particle_capacity + pbc_capacity), 3);
+                    force = (double *) realloc(force, (sizeof(double) * ((particle_capacity + pbc_capacity) * 3)));
+                    ps->updateProperty(3, force, (particle_capacity + pbc_capacity), 3);
+                }
+            }
+        } else {
+            fprintf(stdout, "UpdatePBC\n");
+            fflush(stdout);
+            const double e348 = grid0_d0_max - grid0_d0_min;
+            const double e358 = grid0_d1_max - grid0_d1_min;
+            const double e368 = grid0_d2_max - grid0_d2_min;
+            for(int i7 = 0; i7 < npbc; i7++) {
+                const int e341 = nlocal + i7;
+                const int e342 = e341 * 3;
+                const double p62_0 = position[e342];
+                const int a32 = pbc_map[i7];
+                const int e344 = a32 * 3;
+                const double p63_0 = position[e344];
+                const int e346 = i7 * 3;
+                const int a33 = pbc_mult[e346];
+                const double e349 = a33 * e348;
+                const double e350 = p63_0 + e349;
+                position[e342] = e350;
+                const int e351 = nlocal + i7;
+                const int e352 = e351 * 3;
+                const int e353 = e352 + 1;
+                const double p64_1 = position[e353];
+                const int a34 = pbc_map[i7];
+                const int e354 = a34 * 3;
+                const int e355 = e354 + 1;
+                const double p65_1 = position[e355];
+                const int e356 = i7 * 3;
+                const int e357 = e356 + 1;
+                const int a35 = pbc_mult[e357];
+                const double e359 = a35 * e358;
+                const double e360 = p65_1 + e359;
+                position[e353] = e360;
+                const int e361 = nlocal + i7;
+                const int e362 = e361 * 3;
+                const int e363 = e362 + 2;
+                const double p66_2 = position[e363];
+                const int a36 = pbc_map[i7];
+                const int e364 = a36 * 3;
+                const int e365 = e364 + 2;
+                const double p67_2 = position[e365];
+                const int e366 = i7 * 3;
+                const int e367 = e366 + 2;
+                const int a37 = pbc_mult[e367];
+                const double e369 = a37 * e368;
+                const double e370 = p67_2 + e369;
+                position[e363] = e370;
+            }
+        }
+        const int e456 = i14 % 20;
+        const bool e457 = e456 == 0;
+        if(e457) {
+            fprintf(stdout, "CellListsBuild\n");
+            fflush(stdout);
+            resize = 1;
+            while((resize > 0)) {
+                resize = 0;
+                for(int i8 = 0; i8 < ncells; i8++) {
+                    cell_sizes[i8] = 0;
+                }
+                const int e511 = nlocal + npbc;
+                for(int i9 = 0; i9 < e511; i9++) {
+                    const int e372 = i9 * 3;
+                    const double p68_0 = position[e372];
+                    const double e374 = p68_0 - grid0_d0_min;
+                    const double e375 = e374 / 2.8;
+                    const int e376 = i9 * 3;
+                    const int e377 = e376 + 1;
+                    const double p69_1 = position[e377];
+                    const double e378 = p69_1 - grid0_d1_min;
+                    const double e379 = e378 / 2.8;
+                    const int e380 = i9 * 3;
+                    const int e381 = e380 + 2;
+                    const double p70_2 = position[e381];
+                    const double e382 = p70_2 - grid0_d2_min;
+                    const double e383 = e382 / 2.8;
+                    const int a39 = dim_cells[1];
+                    const int e384 = (int)(e375) * a39;
+                    const int e385 = e384 + (int)(e379);
+                    const int a40 = dim_cells[2];
+                    const int e386 = e385 * a40;
+                    const int e387 = e386 + (int)(e383);
+                    const bool e388 = e387 >= 0;
+                    const bool e389 = e387 <= ncells;
+                    const bool e390 = e388 && e389;
+                    if(e390) {
+                        const int a41 = cell_sizes[e387];
+                        const bool e391 = a41 >= cell_capacity;
+                        if(e391) {
+                            resize = a41;
+                        } else {
+                            const int e392 = e387 * cell_capacity;
+                            const int e393 = e392 + a41;
+                            cell_particles[e393] = i9;
+                            particle_cell[i9] = e387;
+                        }
+                        const int e394 = a41 + 1;
+                        cell_sizes[e387] = e394;
+                    }
+                }
+                const bool e395 = resize > 0;
+                if(e395) {
+                    fprintf(stdout, "Resize cell_capacity\n");
+                    fflush(stdout);
+                    const int e396 = resize * 2;
+                    cell_capacity = e396;
+                    cell_particles = (int *) realloc(cell_particles, (sizeof(int) * (ncells_capacity * cell_capacity)));
+                }
+            }
+        }
+        const int e458 = i14 % 20;
+        const bool e459 = e458 == 0;
+        if(e459) {
+            fprintf(stdout, "NeighborListsBuild\n");
+            fflush(stdout);
+            resize = 1;
+            while((resize > 0)) {
+                resize = 0;
+                for(int i10 = 0; i10 < nlocal; i10++) {
+                    numneighs[i10] = 0;
+                    for(int i11 = 0; i11 < nstencil; i11++) {
+                        const int a46 = particle_cell[i10];
+                        const int a47 = stencil[i11];
+                        const int e400 = a46 + a47;
+                        const bool e401 = e400 >= 0;
+                        const bool e402 = e400 <= ncells;
+                        const bool e403 = e401 && e402;
+                        if(e403) {
+                            const int e404 = e400 * cell_capacity;
+                            const int e412 = i10 * 3;
+                            const double p71_0 = position[e412];
+                            const int e421 = i10 * 3;
+                            const int e422 = e421 + 1;
+                            const double p71_1 = position[e422];
+                            const int e431 = i10 * 3;
+                            const int e432 = e431 + 2;
+                            const double p71_2 = position[e432];
+                            const int e439 = i10 * neighborlist_capacity;
+                            const int a48 = cell_sizes[e400];
+                            for(int i12 = 0; i12 < a48; i12++) {
+                                const int e405 = e404 + i12;
+                                const int a49 = cell_particles[e405];
+                                const bool e406 = a49 != i10;
+                                if(e406) {
+                                    const int e414 = a49 * 3;
+                                    const double p72_0 = position[e414];
+                                    const int e423 = a49 * 3;
+                                    const int e424 = e423 + 1;
+                                    const double p72_1 = position[e424];
+                                    const int e433 = a49 * 3;
+                                    const int e434 = e433 + 2;
+                                    const double p72_2 = position[e434];
+                                    const double e407_0 = p71_0 - p72_0;
+                                    const double e407_1 = p71_1 - p72_1;
+                                    const double e407_2 = p71_2 - p72_2;
+                                    const double e416 = e407_0 * e407_0;
+                                    const double e425 = e407_1 * e407_1;
+                                    const double e426 = e416 + e425;
+                                    const double e435 = e407_2 * e407_2;
+                                    const double e436 = e426 + e435;
+                                    const bool e437 = e436 < 2.8;
+                                    if(e437) {
+                                        const int a50 = numneighs[i10];
+                                        const bool e438 = a50 >= neighborlist_capacity;
+                                        if(e438) {
+                                            resize = a50;
+                                        } else {
+                                            const int e440 = e439 + a50;
+                                            neighborlists[e440] = a49;
+                                            const int e441 = a50 + 1;
+                                            numneighs[i10] = e441;
+                                        }
+                                    }
+                                }
+                            }
+                        }
+                    }
+                }
+                const bool e442 = resize > 0;
+                if(e442) {
+                    fprintf(stdout, "Resize neighborlist_capacity\n");
+                    fflush(stdout);
+                    const int e443 = resize * 2;
+                    neighborlist_capacity = e443;
+                    neighborlists = (int *) realloc(neighborlists, (sizeof(int) * (particle_capacity * neighborlist_capacity)));
+                }
+            }
+        }
+        fprintf(stdout, "PropertiesResetVolatile\n");
+        fflush(stdout);
+        for(int i13 = 0; i13 < nlocal; i13++) {
+            const int e446 = i13 * 3;
+            const double p73_0 = force[e446];
+            const int e448 = i13 * 3;
+            const int e449 = e448 + 1;
+            const double p73_1 = force[e449];
+            const int e450 = i13 * 3;
+            const int e451 = e450 + 2;
+            const double p73_2 = force[e451];
+            force[e446] = 0.0;
+            force[e449] = 0.0;
+            force[e451] = 0.0;
+        }
+        for(int i0 = 0; i0 < nlocal; i0++) {
+            const int e1 = i0 * neighborlist_capacity;
+            const int e47 = i0 * 3;
+            const double p0_0 = position[e47];
+            const int e55 = i0 * 3;
+            const int e56 = e55 + 1;
+            const double p0_1 = position[e56];
+            const int e63 = i0 * 3;
+            const int e64 = e63 + 2;
+            const double p0_2 = position[e64];
+            const int e51 = i0 * 3;
+            const int e59 = i0 * 3;
+            const int e60 = e59 + 1;
+            const int e67 = i0 * 3;
+            const int e68 = e67 + 2;
+            const int a6 = numneighs[i0];
+            for(int i1 = 0; i1 < a6; i1++) {
+                const int e2 = e1 + i1;
+                const int a7 = neighborlists[e2];
+                const int e49 = a7 * 3;
+                const double p1_0 = position[e49];
+                const int e57 = a7 * 3;
+                const int e58 = e57 + 1;
+                const double p1_1 = position[e58];
+                const int e65 = a7 * 3;
+                const int e66 = e65 + 2;
+                const double p1_2 = position[e66];
+                const double e3_0 = p0_0 - p1_0;
+                const double e3_1 = p0_1 - p1_1;
+                const double e3_2 = p0_2 - p1_2;
+                const double e12 = e3_0 * e3_0;
+                const double e21 = e3_1 * e3_1;
+                const double e22 = e12 + e21;
+                const double e31 = e3_2 * e3_2;
+                const double e32 = e22 + e31;
+                const bool e33 = e32 < 2.5;
+                if(e33) {
+                    const double e34 = 1.0 / e32;
+                    const double e35 = e34 * e34;
+                    const double e36 = e35 * e34;
+                    const double p2_0 = force[e51];
+                    const double p2_1 = force[e60];
+                    const double p2_2 = force[e68];
+                    const double e40 = e36 - 0.5;
+                    const double e507 = 48.0 * e36;
+                    const double e508 = e507 * e40;
+                    const double e509 = e508 * e34;
+                    const double e43_0 = e3_0 * e509;
+                    const double e43_1 = e3_1 * e509;
+                    const double e43_2 = e3_2 * e509;
+                    const double e44_0 = p2_0 + e43_0;
+                    const double e44_1 = p2_1 + e43_1;
+                    const double e44_2 = p2_2 + e43_2;
+                    force[e51] = e44_0;
+                    force[e60] = e44_1;
+                    force[e68] = e44_2;
+                }
+            }
+        }
+        for(int i2 = 0; i2 < nlocal; i2++) {
+            const int e76 = i2 * 3;
+            const double p3_0 = velocity[e76];
+            const int e82 = i2 * 3;
+            const int e83 = e82 + 1;
+            const double p3_1 = velocity[e83];
+            const int e88 = i2 * 3;
+            const int e89 = e88 + 2;
+            const double p3_2 = velocity[e89];
+            const int e74 = i2 * 3;
+            const double p4_0 = force[e74];
+            const int e80 = i2 * 3;
+            const int e81 = e80 + 1;
+            const double p4_1 = force[e81];
+            const int e86 = i2 * 3;
+            const int e87 = e86 + 2;
+            const double p4_2 = force[e87];
+            const double e69_0 = 0.005 * p4_0;
+            const double e69_1 = 0.005 * p4_1;
+            const double e69_2 = 0.005 * p4_2;
+            const double p5 = mass[i2];
+            const double e70_0 = e69_0 / p5;
+            const double e70_1 = e69_1 / p5;
+            const double e70_2 = e69_2 / p5;
+            const double e71_0 = p3_0 + e70_0;
+            const double e71_1 = p3_1 + e70_1;
+            const double e71_2 = p3_2 + e70_2;
+            velocity[e76] = e71_0;
+            velocity[e83] = e71_1;
+            velocity[e89] = e71_2;
+            const int e96 = i2 * 3;
+            const double p6_0 = position[e96];
+            const int e102 = i2 * 3;
+            const int e103 = e102 + 1;
+            const double p6_1 = position[e103];
+            const int e108 = i2 * 3;
+            const int e109 = e108 + 2;
+            const double p6_2 = position[e109];
+            const int e94 = i2 * 3;
+            const double p7_0 = velocity[e94];
+            const int e100 = i2 * 3;
+            const int e101 = e100 + 1;
+            const double p7_1 = velocity[e101];
+            const int e106 = i2 * 3;
+            const int e107 = e106 + 2;
+            const double p7_2 = velocity[e107];
+            const double e90_0 = 0.005 * p7_0;
+            const double e90_1 = 0.005 * p7_1;
+            const double e90_2 = 0.005 * p7_2;
+            const double e91_0 = p6_0 + e90_0;
+            const double e91_1 = p6_1 + e90_1;
+            const double e91_2 = p6_2 + e90_2;
+            position[e96] = e91_0;
+            position[e103] = e91_1;
+            position[e109] = e91_2;
+        }
+        const int e461 = nlocal + npbc;
+        const int e460 = i14 + 1;
+        pairs::vtk_write_data(ps, "output/test_local", 0, nlocal, e460);
+        pairs::vtk_write_data(ps, "output/test_pbc", nlocal, e461, e460);
+    }
+}
diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 0000000000000000000000000000000000000000..374b58cbf4636f1e28bacf987ac2fe89ed27ccba
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,6 @@
+[build-system]
+requires = [
+    "setuptools>=42",
+    "wheel"
+]
+build-backend = "setuptools.build_meta"
diff --git a/setup.py b/setup.py
new file mode 100644
index 0000000000000000000000000000000000000000..c9a1f5f637fa85a78e6db81a96684a698b1e69d1
--- /dev/null
+++ b/setup.py
@@ -0,0 +1,46 @@
+import setuptools
+
+
+def readme():
+    with open('README.md') as f:
+        return f.read()
+
+
+modules = [
+    'code_gen',
+    'coupling',
+    'graph',
+    'ir',
+    'mapping',
+    'runtime',
+    'sim',
+    'transformations',
+]
+
+setuptools.setup(name='pairs',
+    description="A code generator for particle simulations",
+    version="0.0.1",
+    long_description=readme(),
+    long_description_content_type="text/markdown",
+    author="Rafael Ravedutti Lucio Machado",
+    license='MIT',
+    author_email="rafael.r.ravedutti@fau.de",
+    url="https://github.com/rafaelravedutti/pairs",
+    install_requires=[],
+    packages=['pairs'] + [f"pairs.{mod}" for mod in modules],
+    package_dir={'pairs': 'src/pairs'},
+    package_data={'pairs.runtime': ['runtime/*.hpp']},
+    classifiers=[
+        "Programming Language :: Python :: 3",
+        "License :: OSI Approved :: MIT License",
+        "Operating System :: OS Independent",
+    ],
+    project_urls={
+        "Bug Tracker": "https://github.com/rafaelravedutti/pairs",
+        "Documentation": "https://github.com/rafaelravedutti/pairs",
+        "Source Code": "https://github.com/rafaelravedutti/pairs",
+    },
+    extras_require={},
+    tests_require=[],
+    python_requires=">=3.6",
+)
diff --git a/code_gen/__init__.py b/src/__init__.py
similarity index 100%
rename from code_gen/__init__.py
rename to src/__init__.py
diff --git a/pairs.py b/src/pairs/__init__.py
similarity index 53%
rename from pairs.py
rename to src/pairs/__init__.py
index d816433b91e4b1ca5260685c9bd7c971eac7b0b3..5a71e5818a272123daa5491b5a079eaf127dfef5 100644
--- a/pairs.py
+++ b/src/pairs/__init__.py
@@ -1,5 +1,5 @@
-from code_gen.cgen import CGen
-from sim.particle_simulation import ParticleSimulation
+from pairs.code_gen.cgen import CGen
+from pairs.sim.particle_simulation import ParticleSimulation
 
 
 def simulation(ref, dims=3, timesteps=100):
diff --git a/coupling/__init__.py b/src/pairs/code_gen/__init__.py
similarity index 100%
rename from coupling/__init__.py
rename to src/pairs/code_gen/__init__.py
diff --git a/code_gen/cgen.py b/src/pairs/code_gen/cgen.py
similarity index 93%
rename from code_gen/cgen.py
rename to src/pairs/code_gen/cgen.py
index 93ea927c2a9aff76d6a56663fc2e761c420a6a9e..219e9e950597390529636c53f7eaaed917d281e0 100644
--- a/code_gen/cgen.py
+++ b/src/pairs/code_gen/cgen.py
@@ -1,23 +1,23 @@
-from ir.assign import Assign
-from ir.arrays import Array, ArrayAccess, ArrayDecl
-from ir.block import Block
-from ir.branches import Branch
-from ir.cast import Cast
-from ir.bin_op import BinOp, Decl, VectorAccess
-from ir.data_types import Type_Int, Type_Float, Type_String, Type_Vector
-from ir.functions import Call
-from ir.layouts import Layout_AoS, Layout_SoA, Layout_Invalid
-from ir.lit import Lit
-from ir.loops import For, Iter, ParticleFor, While
-from ir.math import Ceil, Sqrt
-from ir.memory import Malloc, Realloc
-from ir.properties import Property, PropertyAccess, PropertyList, RegisterProperty, UpdateProperty
-from ir.select import Select
-from ir.sizeof import Sizeof
-from ir.utils import Print
-from ir.variables import Var, VarDecl
-from sim.timestep import Timestep
-from code_gen.printer import Printer
+from pairs.ir.assign import Assign
+from pairs.ir.arrays import Array, ArrayAccess, ArrayDecl
+from pairs.ir.block import Block
+from pairs.ir.branches import Branch
+from pairs.ir.cast import Cast
+from pairs.ir.bin_op import BinOp, Decl, VectorAccess
+from pairs.ir.data_types import Type_Int, Type_Float, Type_String, Type_Vector
+from pairs.ir.functions import Call
+from pairs.ir.layouts import Layout_AoS, Layout_SoA, Layout_Invalid
+from pairs.ir.lit import Lit
+from pairs.ir.loops import For, Iter, ParticleFor, While
+from pairs.ir.math import Ceil, Sqrt
+from pairs.ir.memory import Malloc, Realloc
+from pairs.ir.properties import Property, PropertyAccess, PropertyList, RegisterProperty, UpdateProperty
+from pairs.ir.select import Select
+from pairs.ir.sizeof import Sizeof
+from pairs.ir.utils import Print
+from pairs.ir.variables import Var, VarDecl
+from pairs.sim.timestep import Timestep
+from pairs.code_gen.printer import Printer
 
 
 class CGen:
diff --git a/code_gen/printer.py b/src/pairs/code_gen/printer.py
similarity index 100%
rename from code_gen/printer.py
rename to src/pairs/code_gen/printer.py
diff --git a/ir/__init__.py b/src/pairs/coupling/__init__.py
similarity index 100%
rename from ir/__init__.py
rename to src/pairs/coupling/__init__.py
diff --git a/coupling/parse_cpp.py b/src/pairs/coupling/parse_cpp.py
similarity index 97%
rename from coupling/parse_cpp.py
rename to src/pairs/coupling/parse_cpp.py
index 726475ef663cea2cff65e04384bddaef0ad04a01..8c94702a2e8140f32c02a012258eaa072355d946 100644
--- a/coupling/parse_cpp.py
+++ b/src/pairs/coupling/parse_cpp.py
@@ -1,10 +1,10 @@
-from ast.block import Block
-from ast.branches import Branch
-from ast.data_types import Type_Float, Type_Vector
-from ast.math import Sqrt
-from ast.select import Select
 import clang.cindex
 from clang.cindex import CursorKind as kind
+from pairs.ast.block import Block
+from pairs.ast.branches import Branch
+from pairs.ast.data_types import Type_Float, Type_Vector
+from pairs.ast.math import Sqrt
+from pairs.ast.select import Select
 
 
 def print_tree(node, indent=0):
diff --git a/sim/__init__.py b/src/pairs/graph/__init__.py
similarity index 100%
rename from sim/__init__.py
rename to src/pairs/graph/__init__.py
diff --git a/graph/graphviz.py b/src/pairs/graph/graphviz.py
similarity index 86%
rename from graph/graphviz.py
rename to src/pairs/graph/graphviz.py
index e078b421b5875b6fbcf8dc924aaf478332785c4e..535abe4840a2f793dcd8e4c294f5dbf684939e96 100644
--- a/graph/graphviz.py
+++ b/src/pairs/graph/graphviz.py
@@ -1,11 +1,11 @@
-from ir.arrays import Array
-from ir.bin_op import BinOp, Decl
-from ir.lit import Lit
-from ir.loops import Iter
-from ir.properties import Property
-from ir.variables import Var
-from ir.visitor import Visitor
 from graphviz import Digraph
+from pairs.ir.arrays import Array
+from pairs.ir.bin_op import BinOp, Decl
+from pairs.ir.lit import Lit
+from pairs.ir.loops import Iter
+from pairs.ir.properties import Property
+from pairs.ir.variables import Var
+from pairs.ir.visitor import Visitor
 
 
 class ASTGraph:
diff --git a/ir/operators.py b/src/pairs/ir/__init__.py
similarity index 100%
rename from ir/operators.py
rename to src/pairs/ir/__init__.py
diff --git a/ir/arrays.py b/src/pairs/ir/arrays.py
similarity index 94%
rename from ir/arrays.py
rename to src/pairs/ir/arrays.py
index a49e39a952247f6adefac8530a79eaf21dbd3faf..49aeb88718aad1ff486d5209a89e18e4424b3ab0 100644
--- a/ir/arrays.py
+++ b/src/pairs/ir/arrays.py
@@ -1,12 +1,12 @@
-from ir.assign import Assign
-from ir.ast_node import ASTNode
-from ir.bin_op import BinOp, ASTTerm
-from ir.data_types import Type_Array
-from ir.layouts import Layout_AoS, Layout_SoA
-from ir.lit import as_lit_ast
-from ir.memory import Realloc
-from ir.variables import Var
 from functools import reduce
+from pairs.ir.assign import Assign
+from pairs.ir.ast_node import ASTNode
+from pairs.ir.bin_op import BinOp, ASTTerm
+from pairs.ir.data_types import Type_Array
+from pairs.ir.layouts import Layout_AoS, Layout_SoA
+from pairs.ir.lit import as_lit_ast
+from pairs.ir.memory import Realloc
+from pairs.ir.variables import Var
 
 
 class Arrays:
diff --git a/ir/assign.py b/src/pairs/ir/assign.py
similarity index 83%
rename from ir/assign.py
rename to src/pairs/ir/assign.py
index 4465699badf5f75b9916bf2953aac56ef2eff0e4..97af6ea5045ff50c9c8bec28e2f61e84e53de207 100644
--- a/ir/assign.py
+++ b/src/pairs/ir/assign.py
@@ -1,8 +1,8 @@
-from ir.ast_node import ASTNode
-from ir.data_types import Type_Vector
-from ir.lit import as_lit_ast
-from ir.vector_expr import VectorExpression
 from functools import reduce
+from pairs.ir.ast_node import ASTNode
+from pairs.ir.data_types import Type_Vector
+from pairs.ir.lit import as_lit_ast
+from pairs.ir.vector_expr import VectorExpression
 
 
 class Assign(ASTNode):
diff --git a/ir/ast_node.py b/src/pairs/ir/ast_node.py
similarity index 88%
rename from ir/ast_node.py
rename to src/pairs/ir/ast_node.py
index cbfc0f44dfa7360ad0dbcd998ad6794183b24b48..94e75fa1b260d4b0d6456d1b1c8c147102398873 100644
--- a/ir/ast_node.py
+++ b/src/pairs/ir/ast_node.py
@@ -1,4 +1,4 @@
-from ir.data_types import Type_Invalid
+from pairs.ir.data_types import Type_Invalid
 
 
 class ASTNode:
diff --git a/ir/bin_op.py b/src/pairs/ir/bin_op.py
similarity index 96%
rename from ir/bin_op.py
rename to src/pairs/ir/bin_op.py
index 46c3e4deb8caf619236447a5b798ecc755f6839d..95f48dc000aa6b6129b8ba1c420be61bba1d5317 100644
--- a/ir/bin_op.py
+++ b/src/pairs/ir/bin_op.py
@@ -1,8 +1,8 @@
-from ir.ast_node import ASTNode
-from ir.assign import Assign
-from ir.data_types import Type_Float, Type_Bool, Type_Vector
-from ir.lit import as_lit_ast
-from ir.vector_expr import VectorExpression
+from pairs.ir.ast_node import ASTNode
+from pairs.ir.assign import Assign
+from pairs.ir.data_types import Type_Float, Type_Bool, Type_Vector
+from pairs.ir.lit import as_lit_ast
+from pairs.ir.vector_expr import VectorExpression
 
 
 class Decl(ASTNode):
diff --git a/ir/block.py b/src/pairs/ir/block.py
similarity index 97%
rename from ir/block.py
rename to src/pairs/ir/block.py
index b136db8b0c3afb826cf92d3c6824428f94cb97db..853de45289b9680abe943d5619fa4f45d3611bc8 100644
--- a/ir/block.py
+++ b/src/pairs/ir/block.py
@@ -1,4 +1,4 @@
-from ir.ast_node import ASTNode
+from pairs.ir.ast_node import ASTNode
 
 
 class Block(ASTNode):
diff --git a/ir/branches.py b/src/pairs/ir/branches.py
similarity index 92%
rename from ir/branches.py
rename to src/pairs/ir/branches.py
index 224e67146f48d7f5367ae82567d25f0ad1616249..93962ada8da15565684d5bc7ac43d1c1c1193f50 100644
--- a/ir/branches.py
+++ b/src/pairs/ir/branches.py
@@ -1,6 +1,6 @@
-from ir.ast_node import ASTNode
-from ir.block import Block
-from ir.lit import as_lit_ast
+from pairs.ir.ast_node import ASTNode
+from pairs.ir.block import Block
+from pairs.ir.lit import as_lit_ast
 
 
 class Branch(ASTNode):
diff --git a/ir/cast.py b/src/pairs/ir/cast.py
similarity index 86%
rename from ir/cast.py
rename to src/pairs/ir/cast.py
index 7207b10e0396ed09810bad0bef1c9121a1b35cdd..061eb47e88f80381cd0cd7030b210381c4b741d4 100644
--- a/ir/cast.py
+++ b/src/pairs/ir/cast.py
@@ -1,5 +1,5 @@
-from ir.bin_op import ASTTerm
-from ir.data_types import Type_Int, Type_Float
+from pairs.ir.bin_op import ASTTerm
+from pairs.ir.data_types import Type_Int, Type_Float
 
 
 class Cast(ASTTerm):
diff --git a/ir/data_types.py b/src/pairs/ir/data_types.py
similarity index 100%
rename from ir/data_types.py
rename to src/pairs/ir/data_types.py
diff --git a/ir/functions.py b/src/pairs/ir/functions.py
similarity index 85%
rename from ir/functions.py
rename to src/pairs/ir/functions.py
index 046cd77ba4a898db0c1e3ec4239db54f62f8e4ae..8b01b90856327a0423484bb9b5305d10a637de07 100644
--- a/ir/functions.py
+++ b/src/pairs/ir/functions.py
@@ -1,6 +1,6 @@
-from ir.bin_op import ASTTerm
-from ir.data_types import Type_Int, Type_Invalid
-from ir.lit import as_lit_ast
+from pairs.ir.bin_op import ASTTerm
+from pairs.ir.data_types import Type_Int, Type_Invalid
+from pairs.ir.lit import as_lit_ast
 
 
 class Call(ASTTerm):
diff --git a/ir/layouts.py b/src/pairs/ir/layouts.py
similarity index 100%
rename from ir/layouts.py
rename to src/pairs/ir/layouts.py
diff --git a/ir/lit.py b/src/pairs/ir/lit.py
similarity index 88%
rename from ir/lit.py
rename to src/pairs/ir/lit.py
index ab6215714fa3e2f287d36c7c48a7398e52fe4b91..ff738b5a1dc0d8fba8069cca8f48097a94fd7cf3 100644
--- a/ir/lit.py
+++ b/src/pairs/ir/lit.py
@@ -1,5 +1,5 @@
-from ir.ast_node import ASTNode
-from ir.data_types import Type_Invalid, Type_Int, Type_Float, Type_Bool, Type_String, Type_Vector
+from pairs.ir.ast_node import ASTNode
+from pairs.ir.data_types import Type_Invalid, Type_Int, Type_Float, Type_Bool, Type_String, Type_Vector
 
 
 def is_literal(a):
diff --git a/ir/loops.py b/src/pairs/ir/loops.py
similarity index 93%
rename from ir/loops.py
rename to src/pairs/ir/loops.py
index 0243dcf4134035d9bb543aee7ba6f420cc0c456b..40d73d228fee13ef80286e14cbd7e24a6637cf7c 100644
--- a/ir/loops.py
+++ b/src/pairs/ir/loops.py
@@ -1,9 +1,9 @@
-from ir.ast_node import ASTNode
-from ir.bin_op import BinOp, ASTTerm
-from ir.block import Block
-from ir.branches import Filter
-from ir.data_types import Type_Int
-from ir.lit import as_lit_ast
+from pairs.ir.ast_node import ASTNode
+from pairs.ir.bin_op import BinOp, ASTTerm
+from pairs.ir.block import Block
+from pairs.ir.branches import Filter
+from pairs.ir.data_types import Type_Int
+from pairs.ir.lit import as_lit_ast
 
 
 class Iter(ASTTerm):
diff --git a/ir/math.py b/src/pairs/ir/math.py
similarity index 89%
rename from ir/math.py
rename to src/pairs/ir/math.py
index 93d6269f5a8dd3427d7fd1e48ddfbd19b23225af..9ef762f1f75d4071baf651d9b5f288dbfe0447a7 100644
--- a/ir/math.py
+++ b/src/pairs/ir/math.py
@@ -1,5 +1,5 @@
-from ir.bin_op import ASTTerm
-from ir.data_types import Type_Int, Type_Float
+from pairs.ir.bin_op import ASTTerm
+from pairs.ir.data_types import Type_Int, Type_Float
 
 
 class Sqrt(ASTTerm):
diff --git a/ir/memory.py b/src/pairs/ir/memory.py
similarity index 88%
rename from ir/memory.py
rename to src/pairs/ir/memory.py
index 21b0ace6c30b70c22b38e5158ade203b39fe5128..72cd1f0617af41f6eb994c73a28746294a0fd2b1 100644
--- a/ir/memory.py
+++ b/src/pairs/ir/memory.py
@@ -1,8 +1,8 @@
-from ir.ast_node import ASTNode
-from ir.bin_op import BinOp
-from ir.sizeof import Sizeof
 from functools import reduce
 import operator
+from pairs.ir.ast_node import ASTNode
+from pairs.ir.bin_op import BinOp
+from pairs.ir.sizeof import Sizeof
 
 
 class Malloc(ASTNode):
diff --git a/ir/mutator.py b/src/pairs/ir/mutator.py
similarity index 100%
rename from ir/mutator.py
rename to src/pairs/ir/mutator.py
diff --git a/src/pairs/ir/operators.py b/src/pairs/ir/operators.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/ir/properties.py b/src/pairs/ir/properties.py
similarity index 94%
rename from ir/properties.py
rename to src/pairs/ir/properties.py
index 2d5be9225b79cca1edd6e835b099d48aa3d6106f..c17c78a59198f94831281f879565e10c6e3dad3b 100644
--- a/ir/properties.py
+++ b/src/pairs/ir/properties.py
@@ -1,10 +1,10 @@
-from ir.ast_node import ASTNode
-from ir.assign import Assign
-from ir.bin_op import BinOp, Decl, ASTTerm, VectorAccess
-from ir.data_types import Type_Vector
-from ir.layouts import Layout_AoS
-from ir.lit import as_lit_ast
-from ir.vector_expr import VectorExpression
+from pairs.ir.ast_node import ASTNode
+from pairs.ir.assign import Assign
+from pairs.ir.bin_op import BinOp, Decl, ASTTerm, VectorAccess
+from pairs.ir.data_types import Type_Vector
+from pairs.ir.layouts import Layout_AoS
+from pairs.ir.lit import as_lit_ast
+from pairs.ir.vector_expr import VectorExpression
 
 
 class Properties:
diff --git a/ir/select.py b/src/pairs/ir/select.py
similarity index 77%
rename from ir/select.py
rename to src/pairs/ir/select.py
index 0079d2f5a7cd68e681f3f554e3658d9612bc72d1..8e5e3aebbb22f02eae1adaf68ac4d4dba4a09222 100644
--- a/ir/select.py
+++ b/src/pairs/ir/select.py
@@ -1,6 +1,6 @@
-from ir.ast_node import ASTNode
-from ir.bin_op import BinOp
-from ir.lit import as_lit_ast
+from pairs.ir.ast_node import ASTNode
+from pairs.ir.bin_op import BinOp
+from pairs.ir.lit import as_lit_ast
 
 
 class Select(ASTNode):
diff --git a/ir/sizeof.py b/src/pairs/ir/sizeof.py
similarity index 69%
rename from ir/sizeof.py
rename to src/pairs/ir/sizeof.py
index 9d042bb110a53d5629d4fc8efd7980e42bbcf7d0..407231390893e624b09e88623f0da9617bac689c 100644
--- a/ir/sizeof.py
+++ b/src/pairs/ir/sizeof.py
@@ -1,5 +1,5 @@
-from ir.bin_op import ASTTerm
-from ir.data_types import Type_Int
+from pairs.ir.bin_op import ASTTerm
+from pairs.ir.data_types import Type_Int
 
 
 class Sizeof(ASTTerm):
diff --git a/ir/transform.py b/src/pairs/ir/transform.py
similarity index 92%
rename from ir/transform.py
rename to src/pairs/ir/transform.py
index ba738b6070bcbad0ed4090fa5a285c4e6f75725e..4705552594e6ddf4276f84bad8f67a9ad46a16a2 100644
--- a/ir/transform.py
+++ b/src/pairs/ir/transform.py
@@ -1,10 +1,10 @@
-from ir.arrays import ArrayAccess
-from ir.bin_op import BinOp
-from ir.data_types import Type_Int, Type_Vector
-from ir.layouts import Layout_AoS, Layout_SoA
-from ir.lit import Lit
-from ir.loops import Iter
-from ir.properties import Property
+from pairs.ir.arrays import ArrayAccess
+from pairs.ir.bin_op import BinOp
+from pairs.ir.data_types import Type_Int, Type_Vector
+from pairs.ir.layouts import Layout_AoS, Layout_SoA
+from pairs.ir.lit import Lit
+from pairs.ir.loops import Iter
+from pairs.ir.properties import Property
 
 
 class Transform:
diff --git a/ir/utils.py b/src/pairs/ir/utils.py
similarity index 82%
rename from ir/utils.py
rename to src/pairs/ir/utils.py
index a4446f17274d6dded8a43e17335d4627846e4651..8151d68497d8884637f3e8f1657c445a71a1f476 100644
--- a/ir/utils.py
+++ b/src/pairs/ir/utils.py
@@ -1,4 +1,4 @@
-from ir.ast_node import ASTNode
+from pairs.ir.ast_node import ASTNode
 
 
 class Print(ASTNode):
diff --git a/ir/variables.py b/src/pairs/ir/variables.py
similarity index 93%
rename from ir/variables.py
rename to src/pairs/ir/variables.py
index b4b40dacf9b3f45d2a805250f5d74bfdddf607b5..9f456e872ae91e517288f7e006e1bf067e65dd3d 100644
--- a/ir/variables.py
+++ b/src/pairs/ir/variables.py
@@ -1,6 +1,6 @@
-from ir.ast_node import ASTNode
-from ir.assign import Assign
-from ir.bin_op import ASTTerm 
+from pairs.ir.ast_node import ASTNode
+from pairs.ir.assign import Assign
+from pairs.ir.bin_op import ASTTerm 
 
 
 class Variables:
diff --git a/ir/vector_expr.py b/src/pairs/ir/vector_expr.py
similarity index 92%
rename from ir/vector_expr.py
rename to src/pairs/ir/vector_expr.py
index bd3bb0beb9dfd9c7258ad280af34aea9fa7ec6fa..706a07d75c1cae451d78ea7cae6a6b4536569e97 100644
--- a/ir/vector_expr.py
+++ b/src/pairs/ir/vector_expr.py
@@ -1,6 +1,6 @@
-from ir.ast_node import ASTNode
-from ir.data_types import Type_Vector
-from ir.lit import Lit
+from pairs.ir.ast_node import ASTNode
+from pairs.ir.data_types import Type_Vector
+from pairs.ir.lit import Lit
 
 
 class VectorExpression(ASTNode):
diff --git a/ir/visitor.py b/src/pairs/ir/visitor.py
similarity index 100%
rename from ir/visitor.py
rename to src/pairs/ir/visitor.py
diff --git a/src/pairs/mapping/.funcs.py.swp b/src/pairs/mapping/.funcs.py.swp
new file mode 100644
index 0000000000000000000000000000000000000000..a5728d920b85bc0845eb39ccdc3a0ad5fbdf7362
Binary files /dev/null and b/src/pairs/mapping/.funcs.py.swp differ
diff --git a/src/pairs/mapping/__init__.py b/src/pairs/mapping/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/new_syntax.py b/src/pairs/mapping/funcs.py
similarity index 75%
rename from new_syntax.py
rename to src/pairs/mapping/funcs.py
index bdbeb0d5ef215601f87773d4bbbaef6e8b0ef164..7b47ec66e54073f197a78b80edda14a0589e214e 100644
--- a/new_syntax.py
+++ b/src/pairs/mapping/funcs.py
@@ -1,30 +1,7 @@
-import pairs
-from ir.assign import Assign
-from ir.bin_op import BinOp
 import ast
 import inspect
-
-
-def delta(i, j):
-    return position[i] - position[j]
-
-
-def rsq(i, j):
-    dp = delta(i, j)
-    return dp.x() * dp.x() + dp.y() * dp.y() + dp.z() * dp.z()
-
-
-def lj(i, j):
-    sr2 = 1.0 / rsq
-    sr6 = sr2 * sr2 * sr2 * sigma6
-    #f = 48.0 * sr6 * (sr6 - 0.5) * sr2 * epsilon
-    #force[i] += delta * f
-    force[i] += delta * 48.0 * sr6 * (sr6 - 0.5) * sr2 * epsilon
-
-
-def euler(i):
-    velocity[i] += dt * force[i] / mass[i]
-    position[i] += dt * velocity[i]
+from pairs.ir.assign import Assign
+from pairs.ir.bin_op import BinOp
 
 
 class UndefinedSymbol():
@@ -142,7 +119,7 @@ class BuildParticleIR(ast.NodeVisitor):
         return self.visit(node.value)[self.visit(node.slice)]
 
 
-def add_kernel(sim, func, cutoff_radius=None, position=None, symbols={}):
+def compute(sim, func, cutoff_radius=None, position=None, symbols={}):
     src = inspect.getsource(func)
     tree = ast.parse(src, mode='exec')
     #print(ast.dump(ast.parse(src, mode='exec')))
@@ -162,33 +139,9 @@ def add_kernel(sim, func, cutoff_radius=None, position=None, symbols={}):
             ir.visit(tree)
 
     elif nparams == 2:
-        for i, j, delta, rsq in psim.particle_pairs(cutoff_radius, sim.property(position)):
+        for i, j, delta, rsq in sim.particle_pairs(cutoff_radius, sim.property(position)):
             ir.add_symbols({params[0]: i, params[1]: j, 'delta': delta, 'rsq': rsq})
             ir.visit(tree)
 
     else:
         raise Exception(f"Invalid number of parameters: {nparams}")
-
-
-dt = 0.005
-cutoff_radius = 2.5
-skin = 0.3
-sigma = 1.0
-epsilon = 1.0
-sigma6 = sigma ** 6
-
-psim = pairs.simulation("lj_ns")
-psim.add_real_property('mass', 1.0)
-psim.add_vector_property('position')
-psim.add_vector_property('velocity')
-psim.add_vector_property('force', vol=True)
-psim.from_file("data/minimd_setup_4x4x4.input", ['mass', 'position', 'velocity'])
-psim.create_cell_lists(2.8, 2.8)
-psim.create_neighbor_lists()
-psim.periodic(2.8)
-psim.vtk_output("output/test")
-
-add_kernel(psim, lj, cutoff_radius, 'position', {'sigma6': sigma6, 'epsilon': epsilon})
-add_kernel(psim, euler, symbols={'dt': dt})
-
-psim.generate()
diff --git a/src/pairs/runtime/__init__.py b/src/pairs/runtime/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/runtime/pairs.hpp b/src/pairs/runtime/pairs.hpp
similarity index 100%
rename from runtime/pairs.hpp
rename to src/pairs/runtime/pairs.hpp
diff --git a/runtime/read_from_file.hpp b/src/pairs/runtime/read_from_file.hpp
similarity index 100%
rename from runtime/read_from_file.hpp
rename to src/pairs/runtime/read_from_file.hpp
diff --git a/runtime/vtk.hpp b/src/pairs/runtime/vtk.hpp
similarity index 100%
rename from runtime/vtk.hpp
rename to src/pairs/runtime/vtk.hpp
diff --git a/src/pairs/sim/__init__.py b/src/pairs/sim/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/sim/arrays.py b/src/pairs/sim/arrays.py
similarity index 82%
rename from sim/arrays.py
rename to src/pairs/sim/arrays.py
index 7406019480a38d11a4c0d73ed9119abff07220f1..232964a9aea286747119b06e5c608590f22aa599 100644
--- a/sim/arrays.py
+++ b/src/pairs/sim/arrays.py
@@ -1,5 +1,5 @@
-from ir.memory import Malloc
-from ir.arrays import ArrayDecl
+from pairs.ir.memory import Malloc
+from pairs.ir.arrays import ArrayDecl
 
 
 class ArraysDecl:
diff --git a/sim/cell_lists.py b/src/pairs/sim/cell_lists.py
similarity index 93%
rename from sim/cell_lists.py
rename to src/pairs/sim/cell_lists.py
index 40ec5188e77268a0543a4efbf12bdfebbfaca432..6e94efde7f7c7e1147c8b70ce374bab0fdc757e5 100644
--- a/sim/cell_lists.py
+++ b/src/pairs/sim/cell_lists.py
@@ -1,13 +1,13 @@
-from ir.bin_op import BinOp
-from ir.branches import Branch, Filter
-from ir.cast import Cast
-from ir.data_types import Type_Int
-from ir.math import Ceil
-from ir.loops import For, ParticleFor
-from ir.utils import Print
 from functools import reduce
-from sim.resize import Resize
 import math
+from pairs.ir.bin_op import BinOp
+from pairs.ir.branches import Branch, Filter
+from pairs.ir.cast import Cast
+from pairs.ir.data_types import Type_Int
+from pairs.ir.math import Ceil
+from pairs.ir.loops import For, ParticleFor
+from pairs.ir.utils import Print
+from pairs.sim.resize import Resize
 
 
 class CellLists:
diff --git a/sim/grid.py b/src/pairs/sim/grid.py
similarity index 97%
rename from sim/grid.py
rename to src/pairs/sim/grid.py
index 0eea7ed01ac93abaa86715b0fe107ed980479193..6d71d2f3b9026bdb53de25fe2b40b2aaf9f29fa0 100644
--- a/sim/grid.py
+++ b/src/pairs/sim/grid.py
@@ -1,4 +1,4 @@
-from ir.data_types import Type_Float
+from pairs.ir.data_types import Type_Float
 
 
 class Grid:
diff --git a/sim/kernel_wrapper.py b/src/pairs/sim/kernel_wrapper.py
similarity index 87%
rename from sim/kernel_wrapper.py
rename to src/pairs/sim/kernel_wrapper.py
index 7ac32e36d74894bd5b5928883c8381440eba9fde..ccf7f5adfbe3568048590616af819fc05af77530 100644
--- a/sim/kernel_wrapper.py
+++ b/src/pairs/sim/kernel_wrapper.py
@@ -1,4 +1,4 @@
-from ir.block import Block
+from pairs.ir.block import Block
 
 
 class KernelWrapper():
diff --git a/sim/lattice.py b/src/pairs/sim/lattice.py
similarity index 95%
rename from sim/lattice.py
rename to src/pairs/sim/lattice.py
index b9985e670a729e35a235afe2cf7ddec362c3b47f..d2332c7ebe6e900ffb638f370ad07e96fab4db99 100644
--- a/sim/lattice.py
+++ b/src/pairs/sim/lattice.py
@@ -1,5 +1,5 @@
-from ir.data_types import Type_Vector
-from ir.loops import For
+from pairs.ir.data_types import Type_Vector
+from pairs.ir.loops import For
 
 
 class ParticleLattice():
diff --git a/sim/neighbor_lists.py b/src/pairs/sim/neighbor_lists.py
similarity index 89%
rename from sim/neighbor_lists.py
rename to src/pairs/sim/neighbor_lists.py
index f2fa90b3de80137605a1fdc625f062732792e9a7..f134c1db4f3dc99e0b5d9f03ace5a9600c2e0e51 100644
--- a/sim/neighbor_lists.py
+++ b/src/pairs/sim/neighbor_lists.py
@@ -1,8 +1,8 @@
-from ir.branches import Branch, Filter
-from ir.data_types import Type_Int
-from ir.loops import For, ParticleFor, NeighborFor
-from ir.utils import Print
-from sim.resize import Resize
+from pairs.ir.branches import Branch, Filter
+from pairs.ir.data_types import Type_Int
+from pairs.ir.loops import For, ParticleFor, NeighborFor
+from pairs.ir.utils import Print
+from pairs.sim.resize import Resize
 
 
 class NeighborLists:
diff --git a/sim/particle_simulation.py b/src/pairs/sim/particle_simulation.py
similarity index 82%
rename from sim/particle_simulation.py
rename to src/pairs/sim/particle_simulation.py
index 9bf20a5d28ac07f1277a533e6f55f3395fdc7276..2eb5c89d455f4a80f3d2fb46e369d3783b223af0 100644
--- a/sim/particle_simulation.py
+++ b/src/pairs/sim/particle_simulation.py
@@ -1,29 +1,30 @@
-from ir.arrays import Arrays
-from ir.block import Block
-from ir.branches import Filter
-from ir.data_types import Type_Int, Type_Float, Type_Vector
-from ir.layouts import Layout_AoS
-from ir.loops import ParticleFor, NeighborFor
-from ir.properties import Properties
-from ir.variables import Variables
-from graph.graphviz import ASTGraph
-from sim.arrays import ArraysDecl
-from sim.cell_lists import CellLists, CellListsBuild, CellListsStencilBuild
-from sim.grid import Grid2D, Grid3D
-from sim.kernel_wrapper import KernelWrapper
-from sim.lattice import ParticleLattice
-from sim.neighbor_lists import NeighborLists, NeighborListsBuild
-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
-from sim.vtk import VTKWrite
-from transformations.prioritize_scalar_ops import prioritaze_scalar_ops
-from transformations.set_used_bin_ops import set_used_bin_ops
-from transformations.simplify import simplify_expressions
-from transformations.LICM import move_loop_invariant_code
+from pairs.ir.arrays import Arrays
+from pairs.ir.block import Block
+from pairs.ir.branches import Filter
+from pairs.ir.data_types import Type_Int, Type_Float, Type_Vector
+from pairs.ir.layouts import Layout_AoS
+from pairs.ir.loops import ParticleFor, NeighborFor
+from pairs.ir.properties import Properties
+from pairs.ir.variables import Variables
+from pairs.graph.graphviz import ASTGraph
+from pairs.mapping.funcs import compute
+from pairs.sim.arrays import ArraysDecl
+from pairs.sim.cell_lists import CellLists, CellListsBuild, CellListsStencilBuild
+from pairs.sim.grid import Grid2D, Grid3D
+from pairs.sim.kernel_wrapper import KernelWrapper
+from pairs.sim.lattice import ParticleLattice
+from pairs.sim.neighbor_lists import NeighborLists, NeighborListsBuild
+from pairs.sim.pbc import PBC, UpdatePBC, EnforcePBC, SetupPBC
+from pairs.sim.properties import PropertiesAlloc, PropertiesResetVolatile
+from pairs.sim.read_from_file import ReadFromFile
+from pairs.sim.setup_wrapper import SetupWrapper
+from pairs.sim.timestep import Timestep
+from pairs.sim.variables import VariablesDecl
+from pairs.sim.vtk import VTKWrite
+from pairs.transformations.prioritize_scalar_ops import prioritaze_scalar_ops
+from pairs.transformations.set_used_bin_ops import set_used_bin_ops
+from pairs.transformations.simplify import simplify_expressions
+from pairs.transformations.LICM import move_loop_invariant_code
 
 
 class ParticleSimulation:
@@ -128,6 +129,9 @@ class ParticleSimulation:
         self.properties.add_capacity(self.pbc.pbc_capacity)
         return self.pbc
 
+    def compute(self, sim, func, cutoff_radius=None, position=None, symbols={}):
+        return compute(sim, func, cutoff_radius, position, symbols)
+
     def particle_pairs(self, cutoff_radius=None, position=None):
         self.clear_block()
         for i in ParticleFor(self):
diff --git a/sim/pbc.py b/src/pairs/sim/pbc.py
similarity index 94%
rename from sim/pbc.py
rename to src/pairs/sim/pbc.py
index 1f066be079ec89809b1ec3c8f4eb370823f6506c..c480c4d3a5ebb58468c7ca2a43fdaa9abf7ad2c5 100644
--- a/sim/pbc.py
+++ b/src/pairs/sim/pbc.py
@@ -1,9 +1,9 @@
-from ir.branches import Branch, Filter
-from ir.data_types import Type_Int
-from ir.loops import For, ParticleFor
-from ir.utils import Print
-from ir.select import Select
-from sim.resize import Resize
+from pairs.ir.branches import Branch, Filter
+from pairs.ir.data_types import Type_Int
+from pairs.ir.loops import For, ParticleFor
+from pairs.ir.utils import Print
+from pairs.ir.select import Select
+from pairs.sim.resize import Resize
 
 
 class PBC:
diff --git a/sim/properties.py b/src/pairs/sim/properties.py
similarity index 84%
rename from sim/properties.py
rename to src/pairs/sim/properties.py
index eca36175845b14fe9cffd6f8fbb2a93b7cda35db..7839eed72d759edf9cefcd6f8d708678b716d477 100644
--- a/sim/properties.py
+++ b/src/pairs/sim/properties.py
@@ -1,9 +1,9 @@
+from pairs.ir.data_types import Type_Float, Type_Vector
+from pairs.ir.loops import ParticleFor
+from pairs.ir.memory import Malloc, Realloc
+from pairs.ir.properties import RegisterProperty, UpdateProperty
+from pairs.ir.utils import Print
 from functools import reduce
-from ir.data_types import Type_Float, Type_Vector
-from ir.loops import ParticleFor
-from ir.memory import Malloc, Realloc
-from ir.properties import RegisterProperty, UpdateProperty
-from ir.utils import Print
 import operator
 
 
diff --git a/sim/read_from_file.py b/src/pairs/sim/read_from_file.py
similarity index 81%
rename from sim/read_from_file.py
rename to src/pairs/sim/read_from_file.py
index 6a597a2a50175f3bb2bd0468843e9bafe5bdb345..dc3355bb083a99094079a3aee17f0ab28f086651 100644
--- a/sim/read_from_file.py
+++ b/src/pairs/sim/read_from_file.py
@@ -1,7 +1,7 @@
-from ir.data_types import Type_Float
-from ir.functions import Call_Int
-from ir.properties import PropertyList
-from sim.grid import MutableGrid
+from pairs.ir.data_types import Type_Float
+from pairs.ir.functions import Call_Int
+from pairs.ir.properties import PropertyList
+from pairs.sim.grid import MutableGrid
 
 
 class ReadFromFile():
diff --git a/sim/resize.py b/src/pairs/sim/resize.py
similarity index 84%
rename from sim/resize.py
rename to src/pairs/sim/resize.py
index f862578410e845f0d7d4da522b0e13ba8aeb923a..308ceaa9efad698856e58d8ec94c5edf1249e394 100644
--- a/sim/resize.py
+++ b/src/pairs/sim/resize.py
@@ -1,10 +1,10 @@
+from pairs.ir.branches import Filter
+from pairs.ir.data_types import Type_Int, Type_Float, Type_Vector
+from pairs.ir.loops import While
+from pairs.ir.memory import Realloc
+from pairs.ir.properties import UpdateProperty
+from pairs.ir.utils import Print
 from functools import reduce
-from ir.branches import Filter
-from ir.data_types import Type_Int, Type_Float, Type_Vector
-from ir.loops import While
-from ir.memory import Realloc
-from ir.properties import UpdateProperty
-from ir.utils import Print
 import operator
 
 
diff --git a/sim/setup_wrapper.py b/src/pairs/sim/setup_wrapper.py
similarity index 87%
rename from sim/setup_wrapper.py
rename to src/pairs/sim/setup_wrapper.py
index bc0d5355cb3ee38595831e5a87d280ae3e61d211..629766722bb5546eda98a917b7c6487d048943a0 100644
--- a/sim/setup_wrapper.py
+++ b/src/pairs/sim/setup_wrapper.py
@@ -1,4 +1,4 @@
-from ir.block import Block
+from pairs.ir.block import Block
 
 
 class SetupWrapper():
diff --git a/sim/timestep.py b/src/pairs/sim/timestep.py
similarity index 91%
rename from sim/timestep.py
rename to src/pairs/sim/timestep.py
index f15ee7c32ef6cd3bb1e83df037abd633c83d572e..3e312f281a1657d17bd4244e7860adb1f7a03bba 100644
--- a/sim/timestep.py
+++ b/src/pairs/sim/timestep.py
@@ -1,7 +1,7 @@
-from ir.bin_op import BinOp
-from ir.block import Block
-from ir.branches import Branch
-from ir.loops import For
+from pairs.ir.bin_op import BinOp
+from pairs.ir.block import Block
+from pairs.ir.branches import Branch
+from pairs.ir.loops import For
 
 
 class Timestep:
diff --git a/sim/variables.py b/src/pairs/sim/variables.py
similarity index 85%
rename from sim/variables.py
rename to src/pairs/sim/variables.py
index 4d37e32ca17bc2b2bc385f99437649b67fb1c91a..6869ed1578bf0f93a7e18af5993e34d6c3d96727 100644
--- a/sim/variables.py
+++ b/src/pairs/sim/variables.py
@@ -1,4 +1,4 @@
-from ir.variables import VarDecl
+from pairs.ir.variables import VarDecl
 
 
 class VariablesDecl:
diff --git a/sim/vtk.py b/src/pairs/sim/vtk.py
similarity index 84%
rename from sim/vtk.py
rename to src/pairs/sim/vtk.py
index c499726c43ca20571ae4c18e22f9fa2c090ff2b9..86de55122918cf7d4df1c264861bcc4f0f1ff34f 100644
--- a/sim/vtk.py
+++ b/src/pairs/sim/vtk.py
@@ -1,6 +1,6 @@
-from ir.ast_node import ASTNode
-from ir.functions import Call_Void
-from ir.lit import as_lit_ast
+from pairs.ir.ast_node import ASTNode
+from pairs.ir.functions import Call_Void
+from pairs.ir.lit import as_lit_ast
 
 
 class VTKWrite(ASTNode):
diff --git a/transformations/LICM.py b/src/pairs/transformations/LICM.py
similarity index 97%
rename from transformations/LICM.py
rename to src/pairs/transformations/LICM.py
index 0ac4d37d15d1018d115c1f17a5edd26d37fb258b..264aa37fea89a79fa0716511e45bd1bbb6acd1e5 100644
--- a/transformations/LICM.py
+++ b/src/pairs/transformations/LICM.py
@@ -1,8 +1,8 @@
-from ir.bin_op import BinOp
-from ir.loops import For, While
-from ir.mutator import Mutator
-from ir.properties import PropertyAccess
-from ir.visitor import Visitor
+from pairs.ir.bin_op import BinOp
+from pairs.ir.loops import For, While
+from pairs.ir.mutator import Mutator
+from pairs.ir.properties import PropertyAccess
+from pairs.ir.visitor import Visitor
 
 
 class SetBlockVariants(Mutator):
diff --git a/src/pairs/transformations/__init__.py b/src/pairs/transformations/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/transformations/prioritize_scalar_ops.py b/src/pairs/transformations/prioritize_scalar_ops.py
similarity index 93%
rename from transformations/prioritize_scalar_ops.py
rename to src/pairs/transformations/prioritize_scalar_ops.py
index cf099a8e4497f80f9747d1a5816aa09ffa1a6345..47957f64c78da0edf214af51fadc5b4f9e7da94e 100644
--- a/transformations/prioritize_scalar_ops.py
+++ b/src/pairs/transformations/prioritize_scalar_ops.py
@@ -1,6 +1,6 @@
-from ir.bin_op import BinOp
-from ir.data_types import Type_Float, Type_Vector
-from ir.mutator import Mutator
+from pairs.ir.bin_op import BinOp
+from pairs.ir.data_types import Type_Float, Type_Vector
+from pairs.ir.mutator import Mutator
 
 
 class PrioritazeScalarOps(Mutator):
diff --git a/transformations/set_used_bin_ops.py b/src/pairs/transformations/set_used_bin_ops.py
similarity index 93%
rename from transformations/set_used_bin_ops.py
rename to src/pairs/transformations/set_used_bin_ops.py
index 404cb65d11c6365685499380f459fef0d87cc4f2..ac26bf13f7f2b272b2558817dcfe7d345b3625a4 100644
--- a/transformations/set_used_bin_ops.py
+++ b/src/pairs/transformations/set_used_bin_ops.py
@@ -1,4 +1,4 @@
-from ir.visitor import Visitor
+from pairs.ir.visitor import Visitor
 
 
 class SetUsedBinOps(Visitor):
diff --git a/transformations/simplify.py b/src/pairs/transformations/simplify.py
similarity index 90%
rename from transformations/simplify.py
rename to src/pairs/transformations/simplify.py
index 69803d816b6758ca26bd01f0fa1d7012df69b819..ecba0d6580f239849bf1b813c5073a55bca9427e 100644
--- a/transformations/simplify.py
+++ b/src/pairs/transformations/simplify.py
@@ -1,6 +1,6 @@
-from ir.data_types import Type_Int
-from ir.lit import Lit
-from ir.mutator import Mutator
+from pairs.ir.data_types import Type_Int
+from pairs.ir.lit import Lit
+from pairs.ir.mutator import Mutator
 
 
 class SimplifyExpressions(Mutator):