diff --git a/polyamorphism_optimization/innitialize_sim.py b/polyamorphism_optimization/innitialize_sim.py
new file mode 100644
index 0000000..bd9b470
--- /dev/null
+++ b/polyamorphism_optimization/innitialize_sim.py
@@ -0,0 +1,44 @@
+#!/usr/bin/env python3
+import argparse
+import os
+import subprocess
+from pathlib import Path
+from importlib.resources import files
+import configparser
+
+
+# --- Core function ---
+def innitialize_simulation_dir(directory: Path):
+    os.chdir(directory)
+    
+    params = configparser.ConfigParser()
+    params = params.read('params.ini')
+
+    template_path = files("polyamorphism_optimization.templates") / f"topology_{params['model']}.top"
+    TOPOLOGY = template_path.read_text()
+    template_path = files("polyamorphism_optimization.templates") / "mdp_parameters.mdp"
+    MDP_PARAMETERS = template_path.read_text()
+    
+    # Write files
+    Path("mdp_parameters.mdp").write_text(MDP_PARAMETERS)
+    Path("topology.top").write_text(TOPOLOGY)
+
+    script_dir = Path(__file__).parent
+    subprocess.run([script_dir / "create_emin.sh"], check=True)
+    subprocess.run([script_dir / "create_annealing1.sh"], check=True)
+
+    print(f"finished and simulation directory {directory} created and filled")
+    for f in os.listdir(directory):
+        print(f)
+
+# --- CLI entry point ---
+def main():
+    parser = argparse.ArgumentParser(description="Create simulation input files.")
+    parser.add_argument("-d", "--directory", default=".", help="Base directory")
+    args = parser.parse_args()
+
+    create_simulation_dir(Path(args.directory))
+
+if __name__ == "__main__":
+    main()
+
diff --git a/polyamorphism_optimization/setup_sim.py b/polyamorphism_optimization/setup_sim.py
deleted file mode 100644
index 9f20cc7..0000000
--- a/polyamorphism_optimization/setup_sim.py
+++ /dev/null
@@ -1,110 +0,0 @@
-#!/usr/bin/env python3
-import argparse
-import os
-import subprocess
-from pathlib import Path
-
-# --- Template definitions ---
-MDP_PARAMETERS = """\
-integrator               =  INTEGRATOR
-dt                       = 0.001
-nsteps                   =  NSTEPS
-nstcomm                  = 10
-nstlog                   =  NLOG
-nstcalcenergy            = 10
-nstenergy                =  NENERGY
-nstxout-compressed       =  NSTXOUT
-compressed-x-precision   = 1000
-energygrps               = system
-rlist                    = 1.2
-nstlist                  = 10
-coulombtype              = pme
-rcoulomb                 = 1.2
-coulomb-modifier         = potential-shift-verlet
-vdwtype                  = pme
-rvdw                     = 1.2
-vdw-modifier             = potential-shift-verlet
-fourierspacing           = 0.144
-tcoupl                   = v-rescale
-tau_t                    =  TAUT
-gen_vel                  = yes
-gen_temp                 =  TEMP
-ref_t                    =  TEMP
-tc-grps                  = system
-pcoupl                   =  PCOUPL
-tau_p                    =  TAUP
-ref_p                    =  PRESSURE
-compressibility          = 4.5e-5
-nstpcouple               = 10
-nsttcouple               = 10
-"""
-
-PARAMS_EMIN = """\
-INTEGRATOR = steep
-NSTEPS = 10000
-NLOG = 100
-NENERGY = 100
-NSTXOUT = 100
-TAUT = 1.0
-TEMP = 300
-PCOUPL = no
-TCOUPL = no
-TAUP = 1.0
-"""
-
-PARAMS_AN1 = """\
-INTEGRATOR = md
-NSTEPS = 2000000
-NLOG = 1000
-NENERGY = 100
-NSTXOUT = 1000
-TCOUPL = v-rescale
-TAUT = 0.5
-TEMP = 10.0
-PCOUPL = c-rescale
-TAUP = 1.0
-"""
-
-TOPOLOGY = """\
-[ defaults ]
-; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ
-  1             2               yes              1.0     1.0
-[atomtypes]
-;name     mass      charge   ptype    sigma        epsilon
-A     MASS       0.000       A      SIGMA_A         EPSILON_A
-B     MASS       0.000       A      SIGMA_B         EPSILON_B
-D     0          0.000       D   0.0           0.0
-...
-"""
-
-# --- Core function ---
-def create_simulation_dir(directory: Path):
-    directory.mkdir(parents=True, exist_ok=True)
-    os.chdir(directory)
-
-    # Write files
-    (directory / "mdp_parameters.mdp").write_text(MDP_PARAMETERS)
-    (directory / "params_emin.ini").write_text(PARAMS_EMIN)
-    (directory / "params_an1.ini").write_text(PARAMS_AN1)
-    (directory / "topology.top").write_text(TOPOLOGY)
-    # ... add the rest of your templates here ...
-
-    script_dir = Path(__file__).parent
-    subprocess.run([script_dir / "create_emin.sh"], check=True)
-    subprocess.run([script_dir / "create_annealing1.sh"], check=True)
-
-    print(f"finished and simulation directory {directory} created and filled")
-    for f in os.listdir(directory):
-        print(f)
-
-# --- CLI entry point ---
-def main():
-    parser = argparse.ArgumentParser(description="Create simulation input files.")
-    parser.add_argument("-d", "--directory", default=".", help="Base directory")
-    args = parser.parse_args()
-
-    create_simulation_dir(Path(args.directory))
-
-if __name__ == "__main__":
-    main()
-
diff --git a/polyamorphism_optimization/setup_sim.sh b/polyamorphism_optimization/setup_sim.sh
deleted file mode 100755
index c2ecefc..0000000
--- a/polyamorphism_optimization/setup_sim.sh
+++ /dev/null
@@ -1,182 +0,0 @@
-#!/bin/bash -l
-
-usage() { echo "Usage: $0 [-d basedirectory ]" 1>&2; exit 1; }
-
-DIRECTORY="."
-
-while getopts "d:" opt; do
-  case $opt in
-    d )
-        DIRECTORY=$OPTARG
-        ;;
-    h )
-        usage
-        ;;
-   \? )
-        usage
-        ;;
-  esac
-done
-
-cd $DIRECTORY
-
-#=====================================================================================
-cat > mdp_parameters.mdp << EOF
-integrator               =  INTEGRATOR 
-dt                       = 0.001
-nsteps                   =  NSTEPS 
-nstcomm                  = 10
-
-nstlog                   =  NLOG 
-nstcalcenergy            = 10
-nstenergy                =  NENERGY 
-nstxout-compressed       =  NSTXOUT 
-compressed-x-precision   = 1000
-energygrps               = system
-
-rlist                    = 1.2
-nstlist                  = 10
-
-coulombtype              = pme
-rcoulomb                 = 1.2
-coulomb-modifier         = potential-shift-verlet
-vdwtype                  = pme
-rvdw                     = 1.2
-vdw-modifier             = potential-shift-verlet
-fourierspacing           = 0.144
-
-tcoupl                   = v-rescale
-tau_t                    =  TAUT 
-gen_vel                  = yes
-gen_temp                 =  TEMP 
-ref_t                    =  TEMP 
-tc-grps                  = system
-
-pcoupl                   =  PCOUPL 
-tau_p                    =  TAUP 
-ref_p                    =  PRESSURE 
-compressibility          = 4.5e-5
-
-nstpcouple               = 10
-nsttcouple               = 10
-
-EOF
-
-#=====================================================================================
-cat > params_emin.ini << EOF
-INTEGRATOR = steep
-NSTEPS = 10000
-NLOG = 100
-NENERGY = 100
-NSTXOUT = 100
-TAUT = 1.0
-TEMP = 300
-PCOUPL = no
-TCOUPL = no
-TAUP = 1.0
-EOF
-
-cat > params_an1.ini << EOF
-INTEGRATOR = md
-NSTEPS = 2000000
-NLOG = 1000
-NENERGY = 100
-NSTXOUT = 1000
-TCOUPL = v-rescale
-TAUT = 0.5
-TEMP = 10.0
-PCOUPL = c-rescale
-TAUP = 1.0
-EOF
-
-#=====================================================================================
-
-cat > topology.top << EOF
-[ defaults ]
-; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ
-  1             2               yes              1.0     1.0
-
-[atomtypes]
-;name     mass      charge   ptype    sigma        epsilon
-A     MASS       0.000       A      SIGMA_A         EPSILON_A 
-B     MASS       0.000       A      SIGMA_B         EPSILON_B 
-D     0             0.000       D   0.0           0.0
-
-[moleculetype]
-; name nrexcl
-L  1
-
-[atoms]
-; nr type resnr residu atom cgnr charge
-1     A   1     L  A    1     0        MASS 
-2     B   1     L  B1   1    CHARGE_B      MASS 
-3     B   1     L  B2   1    CHARGE_B      MASS 
-4     D   1     L  D    1    CHARGE_A      0.0
-
-[settles]
-;i      funct   doh     dhh
-1       1       D_AB     D_BB 
-
-
-[dummies3]
-; Dummy from                            funct       a        b
-4       1       2       3       1       DUMMY_AB      DUMMY_AB 
-
-[exclusions]
-1       2       3       4
-2       1       3       4
-3       1       2       4
-4       1       2       3
-
-[system]
- MODEL 
-
-[molecules]
-L                 NMOL 
-EOF
-
-cat > topology_atomistic.top << EOF
-[ defaults ]
-; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ
-  1             2               yes              1.0     1.0
-
-[atomtypes]
-;name     mass      charge   ptype    sigma        epsilon
-A     MASS       0.000       A      SIGMA_A         EPSILON_A 
-B     MASS       0.000       A      SIGMA_B         EPSILON_B 
-
-[moleculetype]
-; name nrexcl
-A  1
-
-[atoms]
-; nr type resnr residu atom cgnr charge
-1     A   1     L  A   1    CHARGE_A      MASS 
-
-[moleculetype]
-; name nrexcl
-B  1
-
-[atoms]
-; nr type resnr residu atom cgnr charge
-1     B   1     L  B   1    CHARGE_B      MASS 
-
-[system]
- MODEL 
-
-[molecules]
-A                1000
-B                2000
-EOF
-
-#=====================================================================================
-echo $(dirname $0)
-bash ".$(dirname $0)/create_emin.sh"
-bash ".$(dirname $0)/create_annealing1.sh"
-
-
-#=====================================================================================
-
-echo "finished and simulation directory $DIRECTORY created and filled"
-ls
-exit 0
diff --git a/polyamorphism_optimization/templates/NPT.ini b/polyamorphism_optimization/templates/NPT.ini
new file mode 100644
index 0000000..7aeb8d1
--- /dev/null
+++ b/polyamorphism_optimization/templates/NPT.ini
@@ -0,0 +1,9 @@
+INTEGRATOR = md
+NSTEPS = 2000000
+NLOG = 1000
+NENERGY = 100
+NSTXOUT = 1000
+TCOUPL = v-rescale
+TAUT = 1.0
+PCOUPL = c-rescale
+TAUP = 1.0
diff --git a/polyamorphism_optimization/templates/NVT.ini b/polyamorphism_optimization/templates/NVT.ini
new file mode 100644
index 0000000..86ed9dd
--- /dev/null
+++ b/polyamorphism_optimization/templates/NVT.ini
@@ -0,0 +1,9 @@
+INTEGRATOR = md
+NSTEPS = 2000000
+NLOG = 1000
+NENERGY = 100
+NSTXOUT = 1000
+TCOUPL = v-rescale
+TAUT = 1.0
+PCOUPL = no
+TAUP = 1.0
diff --git a/polyamorphism_optimization/templates/mdp_parameters.mdp b/polyamorphism_optimization/templates/mdp_parameters.mdp
new file mode 100644
index 0000000..a1607a2
--- /dev/null
+++ b/polyamorphism_optimization/templates/mdp_parameters.mdp
@@ -0,0 +1,31 @@
+integrator               =  INTEGRATOR
+dt                       = 0.001
+nsteps                   =  NSTEPS
+nstcomm                  = 10
+nstlog                   =  NLOG
+nstcalcenergy            = 10
+nstenergy                =  NENERGY
+nstxout-compressed       =  NSTXOUT
+compressed-x-precision   = 1000
+energygrps               = system
+rlist                    = 1.2
+nstlist                  = 10
+coulombtype              = pme
+rcoulomb                 = 1.2
+coulomb-modifier         = potential-shift-verlet
+vdwtype                  = pme
+rvdw                     = 1.2
+vdw-modifier             = potential-shift-verlet
+fourierspacing           = 0.144
+tcoupl                   = v-rescale
+tau_t                    =  TAUT
+gen_vel                  = yes
+gen_temp                 =  TEMP
+ref_t                    =  TEMP
+tc-grps                  = system
+pcoupl                   =  PCOUPL
+tau_p                    =  TAUP
+ref_p                    =  PRESSURE
+compressibility          = 4.5e-5
+nstpcouple               = 10
+nsttcouple               = 10
diff --git a/polyamorphism_optimization/templates/params/annealing1.ini b/polyamorphism_optimization/templates/params/annealing1.ini
new file mode 100644
index 0000000..0d49ba8
--- /dev/null
+++ b/polyamorphism_optimization/templates/params/annealing1.ini
@@ -0,0 +1,10 @@
+INTEGRATOR = md
+NSTEPS = 2000000
+NLOG = 1000
+NENERGY = 100
+NSTXOUT = 1000
+TCOUPL = v-rescale
+TAUT = 0.5
+TEMP = 10.0
+PCOUPL = c-rescale
+TAUP = 1.0
diff --git a/polyamorphism_optimization/templates/params/annealing2.ini b/polyamorphism_optimization/templates/params/annealing2.ini
new file mode 100644
index 0000000..47fee70
--- /dev/null
+++ b/polyamorphism_optimization/templates/params/annealing2.ini
@@ -0,0 +1,10 @@
+INTEGRATOR = md
+NSTEPS = 5000000
+NLOG = 1000
+NENERGY = 100
+NSTXOUT = 1000
+TCOUPL = v-rescale
+TAUT = 0.5
+TEMP = 10.0
+PCOUPL = c-rescale
+TAUP = 1.0
diff --git a/polyamorphism_optimization/templates/params/emin.ini b/polyamorphism_optimization/templates/params/emin.ini
new file mode 100644
index 0000000..e2a8ae4
--- /dev/null
+++ b/polyamorphism_optimization/templates/params/emin.ini
@@ -0,0 +1,10 @@
+INTEGRATOR = steep
+NSTEPS = 10000
+NLOG = 100
+NENERGY = 100
+NSTXOUT = 100
+TAUT = 1.0
+TEMP = 300
+PCOUPL = no
+TCOUPL = no
+TAUP = 1.0
diff --git a/polyamorphism_optimization/templates/top_3point.top b/polyamorphism_optimization/templates/top_3point.top
new file mode 100644
index 0000000..e5cef22
--- /dev/null
+++ b/polyamorphism_optimization/templates/top_3point.top
@@ -0,0 +1,33 @@
+[ defaults ]
+; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ
+  1             2               yes              1.0     1.0
+[atomtypes]
+;name     mass      charge   ptype    sigma        epsilon
+A     MASS       0.000       A      SIGMA_A         EPSILON_A
+B     MASS       0.000       A      SIGMA_B         EPSILON_B
+
+[moleculetype]
+; name nrexcl
+L  1
+
+[atoms]
+; nr type resnr residu atom cgnr charge
+1     A   1     L  A    1    CHARGE_A      MASS 
+2     B   1     L  B1   1    CHARGE_B      MASS 
+3     B   1     L  B2   1    CHARGE_B      MASS 
+
+[settles]
+;i      funct   doh     dhh
+1       1       D_AB     D_BB 
+
+
+[exclusions]
+1       2       3
+2       1       3
+3       1       2
+
+[system]
+ MODEL 
+
+[molecules]
+L                 NMOL 
diff --git a/polyamorphism_optimization/templates/top_4point.top b/polyamorphism_optimization/templates/top_4point.top
new file mode 100644
index 0000000..0bc968b
--- /dev/null
+++ b/polyamorphism_optimization/templates/top_4point.top
@@ -0,0 +1,40 @@
+[ defaults ]
+; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ
+  1             2               yes              1.0     1.0
+[atomtypes]
+;name     mass      charge   ptype    sigma        epsilon
+A     MASS       0.000       A      SIGMA_A         EPSILON_A
+B     MASS       0.000       A      SIGMA_B         EPSILON_B
+D     0          0.000       D   0.0           0.0
+
+[moleculetype]
+; name nrexcl
+L  1
+
+[atoms]
+; nr type resnr residu atom cgnr charge
+1     A   1     L  A    1     0        MASS 
+2     B   1     L  B1   1    CHARGE_B      MASS 
+3     B   1     L  B2   1    CHARGE_B      MASS 
+4     D   1     L  D    1    CHARGE_A      0.0
+
+[settles]
+;i      funct   doh     dhh
+1       1       D_AB     D_BB 
+
+
+[dummies3]
+; Dummy from                            funct       a        b
+4       1       2       3       1       DUMMY_AB      DUMMY_AB 
+
+[exclusions]
+1       2       3       4
+2       1       3       4
+3       1       2       4
+4       1       2       3
+
+[system]
+ MODEL 
+
+[molecules]
+L                 NMOL 
diff --git a/polyamorphism_optimization/templates/top_atomistic.top b/polyamorphism_optimization/templates/top_atomistic.top
new file mode 100644
index 0000000..bb5cd42
--- /dev/null
+++ b/polyamorphism_optimization/templates/top_atomistic.top
@@ -0,0 +1,31 @@
+[ defaults ]
+; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ
+  1             2               yes              1.0     1.0
+
+[atomtypes]
+;name     mass      charge   ptype    sigma        epsilon
+A     MASS       0.000       A      SIGMA_A         EPSILON_A 
+B     MASS       0.000       A      SIGMA_B         EPSILON_B 
+
+[moleculetype]
+; name nrexcl
+A  1
+
+[atoms]
+; nr type resnr residu atom cgnr charge
+1     A   1     L  A   1    CHARGE_A      MASS 
+
+[moleculetype]
+; name nrexcl
+B  1
+
+[atoms]
+; nr type resnr residu atom cgnr charge
+1     B   1     L  B   1    CHARGE_B      MASS 
+
+[system]
+ MODEL 
+
+[molecules]
+A                NMOL
+B                NMOLB