diff --git a/src/osdi/osdi.h b/src/osdi/osdi.h
index d35d58aef..6f633b037 100644
--- a/src/osdi/osdi.h
+++ b/src/osdi/osdi.h
@@ -189,7 +189,7 @@ typedef struct OsdiDescriptor {
                                      double temperature, uint32_t num_terminals,
                                      OsdiSimParas *sim_params, OsdiInitInfo *res);
 
-  uint32_t (*eval)(void *handle, void *inst, void *model, OsdiSimInfo *info);
+  uint32_t (*eval)(void *handle, void *inst, void *model, const OsdiSimInfo *info);
   void (*load_noise)(void *inst, void *model, double freq, double *noise_dens,
                   double *ln_noise_dens);
   void (*load_residual_resist)(void *inst, void* model, double *dst);
diff --git a/src/osdi/osdidefs.h b/src/osdi/osdidefs.h
index 81f8d3e0f..45d117336 100644
--- a/src/osdi/osdidefs.h
+++ b/src/osdi/osdidefs.h
@@ -59,7 +59,7 @@ typedef struct OsdiExtraInstData {
   double temp;
   bool temp_given;
   bool dt_given;
-  bool finish;
+  uint32_t eval_flags;
 
 } ALIGN(MAX_ALIGN) OsdiExtraInstData;
 
diff --git a/src/osdi/osdiload.c b/src/osdi/osdiload.c
index 56c352c9b..e00c40efb 100644
--- a/src/osdi/osdiload.c
+++ b/src/osdi/osdiload.c
@@ -28,11 +28,10 @@ char *sim_params[NUM_SIM_PARAMS + 1] = {
     "gdev", "gmin", "tnom", "simulatorVersion", "sourceScaleFactor", NULL};
 char *sim_params_str[1] = {NULL};
 
-double sim_param_vals[NUM_SIM_PARAMS] = {0,0,0,0,0};
+double sim_param_vals[NUM_SIM_PARAMS] = {0, 0, 0, 0, 0};
 
 /* values returned by $simparam*/
 OsdiSimParas get_simparams(const CKTcircuit *ckt) {
-
   double simulatorVersion = strtod(PACKAGE_VERSION, NULL);
   double gdev = ckt->CKTgmin;
   double sourceScaleFactor = ckt->CKTsrcFact;
@@ -48,6 +47,71 @@ OsdiSimParas get_simparams(const CKTcircuit *ckt) {
   return sim_params_;
 }
 
+static void eval(const OsdiDescriptor *descr, const GENinstance *gen_inst,
+                 void *inst, OsdiExtraInstData *extra_inst_data,
+                 const void *model, const OsdiSimInfo *sim_info) {
+
+  OsdiNgspiceHandle handle =
+      (OsdiNgspiceHandle){.kind = 3, .name = gen_inst->GENname};
+  /* TODO initial conditions? */
+  extra_inst_data->eval_flags = descr->eval(&handle, inst, model, sim_info);
+}
+
+static void load(CKTcircuit *ckt, const GENinstance *gen_inst, void *model,
+                 void *inst, OsdiExtraInstData *extra_inst_data, bool is_tran,
+                 bool is_init_tran, const OsdiDescriptor *descr) {
+  double dump;
+  if (is_tran) {
+    /* load dc matrix and capacitances (charge derivative multiplied with
+     * CKTag[0]) */
+    descr->load_jacobian_tran(inst, model, ckt->CKTag[0]);
+
+    /* load static rhs and dynamic linearized rhs (SUM Vb * dIa/dVb)*/
+    descr->load_spice_rhs_tran(inst, model, ckt->CKTrhs, ckt->CKTrhsOld,
+                               ckt->CKTag[0]);
+
+    uint32_t *node_mapping =
+        (uint32_t *)(((char *)inst) + descr->node_mapping_offset);
+
+    /* use numeric integration to obtain the remainer of the RHS*/
+    int state = gen_inst->GENstate + (int)descr->num_states;
+    for (uint32_t i = 0; i < descr->num_nodes; i++) {
+      if (descr->nodes[i].react_residual_off != UINT32_MAX) {
+
+        double residual_react =
+            *((double *)(((char *)inst) + descr->nodes[i].react_residual_off));
+
+        /* store charges in state vector*/
+        ckt->CKTstate0[state] = residual_react;
+        if (is_init_tran) {
+          ckt->CKTstate1[state] = residual_react;
+        }
+
+        /* we only care about the numeric integration itself not ceq/geq
+        because those are already calculated by load_jacobian_tran and
+        load_spice_rhs_tran*/
+        NIintegrate(ckt, &dump, &dump, 0, state);
+
+        /* add the numeric derivative to the rhs */
+        ckt->CKTrhs[node_mapping[i]] -= ckt->CKTstate0[state + 1];
+
+        if (is_init_tran) {
+          ckt->CKTstate1[state + 1] = ckt->CKTstate0[state + 1];
+        }
+
+        state += 2;
+      }
+    }
+  } else {
+    /* copy internal derivatives into global matrix */
+    descr->load_jacobian_resist(inst, model);
+
+    /* calculate spice RHS from internal currents and store into global RHS
+     */
+    descr->load_spice_rhs_dc(inst, model, ckt->CKTrhs, ckt->CKTrhsOld);
+  }
+}
+
 extern int OSDIload(GENmodel *inModel, CKTcircuit *ckt) {
   OsdiNgspiceHandle handle;
   GENmodel *gen_model;
@@ -113,6 +177,35 @@ extern int OSDIload(GENmodel *inModel, CKTcircuit *ckt) {
 
   OsdiRegistryEntry *entry = osdi_reg_entry_model(inModel);
   const OsdiDescriptor *descr = entry->descriptor;
+  uint32_t eval_flags = 0;
+
+#ifdef USE_OMP
+  /* use openmp 3.0 tasks to parallelize linked list transveral */
+#pragma omp parallel
+#pragma omp single
+  {
+    for (gen_model = inModel; gen_model; gen_model = gen_model->GENnextModel) {
+      void *model = osdi_model_data(gen_model);
+
+      for (gen_inst = gen_model->GENinstances; gen_inst;
+           gen_inst = gen_inst->GENnextInstance) {
+
+        void *inst = osdi_instance_data(entry, gen_inst);
+
+        OsdiExtraInstData *extra_inst_data =
+            osdi_extra_instance_data(entry, gen_inst);
+
+#pragma omp task firstprivate(gen_inst, inst, extra_inst_data, model)
+        eval(descr, gen_inst, inst, extra_inst_data, model, &sim_info);
+      }
+    }
+  }
+
+  /* init small signal analysis does not require loading values into
+   * matrix/rhs*/
+  if (is_init_smsig) {
+    return ret;
+  }
 
   for (gen_model = inModel; gen_model; gen_model = gen_model->GENnextModel) {
     void *model = osdi_model_data(gen_model);
@@ -120,92 +213,49 @@ extern int OSDIload(GENmodel *inModel, CKTcircuit *ckt) {
     for (gen_inst = gen_model->GENinstances; gen_inst;
          gen_inst = gen_inst->GENnextInstance) {
       void *inst = osdi_instance_data(entry, gen_inst);
+      OsdiExtraInstData *extra_inst_data =
+          osdi_extra_instance_data(entry, gen_inst);
+      load(ckt, gen_inst, model, inst, extra_inst_data, is_tran, is_init_tran,
+           descr);
+      eval_flags |= extra_inst_data->eval_flags;
+    }
+  }
+#else
+  for (gen_model = inModel; gen_model; gen_model = gen_model->GENnextModel) {
+    void *model = osdi_model_data(gen_model);
 
-      /* hpyothetically this could run in parallel we do not write any shared
-       data here*/
-      handle = (OsdiNgspiceHandle){.kind = 3, .name = gen_inst->GENname};
-      /* TODO initial conditions? */
-      uint32_t ret_flags = descr->eval(&handle, inst, model, &sim_info);
+    for (gen_inst = gen_model->GENinstances; gen_inst;
+         gen_inst = gen_inst->GENnextInstance) {
+      void *inst = osdi_instance_data(entry, gen_inst);
 
-      /* call to $fatal in Verilog-A abort!*/
-      if (ret_flags & EVAL_RET_FLAG_FATAL) {
-        return E_PANIC;
-      }
+      OsdiExtraInstData *extra_inst_data =
+          osdi_extra_instance_data(entry, gen_inst);
+      eval(descr, gen_inst, inst, extra_inst_data, model, &sim_info);
 
       /* init small signal analysis does not require loading values into
        * matrix/rhs*/
-      if (is_init_smsig) {
-        continue;
-      }
-
-      /* handle calls to $finish, $limit, $stop
-       * TODO actually do something with extra_inst_data->finish  and
-       * extra_inst_data->limt
-       * */
-      OsdiExtraInstData *extra_inst_data =
-          osdi_extra_instance_data(entry, gen_inst);
-      if (ret_flags & EVAL_RET_FLAG_FINISH) {
-        extra_inst_data->finish = true;
-      }
-      if (ret_flags & EVAL_RET_FLAG_LIM) {
-        ckt->CKTnoncon++;
-        ckt->CKTtroubleElt = gen_inst;
-
-      }
-      if (ret_flags & EVAL_RET_FLAG_STOP) {
-        ret = (E_PAUSE);
-      }
-
-      if (is_tran) {
-        /* load dc matrix and capacitances (charge derivative multiplied with
-         * CKTag[0]) */
-        descr->load_jacobian_tran(inst, model, ckt->CKTag[0]);
-
-        /* load static rhs and dynamic linearized rhs (SUM Vb * dIa/dVb)*/
-        descr->load_spice_rhs_tran(inst, model, ckt->CKTrhs, ckt->CKTrhsOld,
-                                   ckt->CKTag[0]);
-
-        uint32_t *node_mapping =
-            (uint32_t *)(((char *)inst) + descr->node_mapping_offset);
-
-        /* use numeric integration to obtain the remainer of the RHS*/
-        int state = gen_inst->GENstate + (int) descr->num_states;
-        for (uint32_t i = 0; i < descr->num_nodes; i++) {
-          if (descr->nodes[i].react_residual_off != UINT32_MAX) {
-
-            double residual_react = *((
-                double *)(((char *)inst) + descr->nodes[i].react_residual_off));
-
-            /* store charges in state vector*/
-            ckt->CKTstate0[state] = residual_react;
-            if (is_init_tran) {
-              ckt->CKTstate1[state] = residual_react;
-            }
-
-            /* we only care about the numeric integration itself not ceq/geq
-            because those are already calculated by load_jacobian_tran and
-            load_spice_rhs_tran*/
-            NIintegrate(ckt, &dump, &dump, 0, state);
-
-            /* add the numeric derivative to the rhs */
-            ckt->CKTrhs[node_mapping[i]] -= ckt->CKTstate0[state + 1];
-
-            if (is_init_tran) {
-              ckt->CKTstate1[state + 1] = ckt->CKTstate0[state + 1];
-            }
-
-            state += 2;
-          }
-        }
-      } else {
-        /* copy internal derivatives into global matrix */
-        descr->load_jacobian_resist(inst, model);
-
-        /* calculate spice RHS from internal currents and store into global RHS
-         */
-        descr->load_spice_rhs_dc(inst, model, ckt->CKTrhs, ckt->CKTrhsOld);
+      if (!is_init_smsig) {
+        load(ckt, gen_inst, model, inst, extra_inst_data, is_tran, is_init_tran,
+             descr);
+        eval_flags |= extra_inst_data->eval_flags;
       }
     }
   }
-  return ret;
+#endif
+
+  /* call to $fatal in Verilog-A abort simulation!*/
+  if (eval_flags & EVAL_RET_FLAG_FATAL) {
+    return E_PANIC;
+  }
+
+  if (eval_flags & EVAL_RET_FLAG_LIM) {
+    ckt->CKTnoncon++;
+    ckt->CKTtroubleElt = gen_inst;
+  }
+
+  if (eval_flags & EVAL_RET_FLAG_STOP) {
+    return E_PAUSE;
+  }
+
+  return OK;
 }