diff --git a/gromacs/includes/parameter.h b/gromacs/includes/parameter.h
index 69d8f3d..4b31aca 100644
--- a/gromacs/includes/parameter.h
+++ b/gromacs/includes/parameter.h
@@ -34,7 +34,7 @@ typedef struct {
     char* param_file;
     char* input_file;
     char* vtk_file;
-    char *xtc_file;
+    char* xtc_file;
     MD_FLOAT epsilon;
     MD_FLOAT sigma;
     MD_FLOAT sigma6;
diff --git a/gromacs/includes/simd/avx512_float.h b/gromacs/includes/simd/avx512_float.h
index ac5ea8c..8043bbe 100644
--- a/gromacs/includes/simd/avx512_float.h
+++ b/gromacs/includes/simd/avx512_float.h
@@ -49,6 +49,13 @@ static inline MD_FLOAT simd_h_reduce_sum(MD_SIMD_FLOAT a) {
     return 0.0;
 }
 
+static inline MD_FLOAT simd_incr_reduced_sum(MD_FLOAT *m, MD_SIMD_FLOAT v0, MD_SIMD_FLOAT v1, MD_SIMD_FLOAT v2, MD_SIMD_FLOAT v3) {
+    // This would only be called in a Mx16 configuration, which is not valid in GROMACS
+    fprintf(stderr, "simd_h_reduce_sum(): Called with AVX512 intrinsics and single-precision which is not valid!\n");
+    exit(-1);
+    return 0.0;
+}
+
 static inline MD_SIMD_FLOAT simd_load_h_duplicate(const float* m) {
     return _mm512_castpd_ps(_mm512_broadcast_f64x4(_mm256_load_pd((const double *)(m))));
 }
@@ -57,7 +64,7 @@ static inline MD_SIMD_FLOAT simd_load_h_dual(const float* m) {
     return _mm512_shuffle_f32x4(_mm512_broadcastss_ps(_mm_load_ss(m)), _mm512_broadcastss_ps(_mm_load_ss(m + 1)), 0x44);
 }
 
-static inline MD_FLOAT simd_h_dual_reduce_sum(float* m, MD_SIMD_FLOAT v0, MD_SIMD_FLOAT v1) {
+static inline MD_FLOAT simd_h_dual_incr_reduced_sum(float* m, MD_SIMD_FLOAT v0, MD_SIMD_FLOAT v1) {
     __m512 t0, t1;
     __m128 t2, t3;
 
diff --git a/gromacs/includes/util.h b/gromacs/includes/util.h
index 4e123ea..5adabe2 100644
--- a/gromacs/includes/util.h
+++ b/gromacs/includes/util.h
@@ -42,10 +42,6 @@
 #define MAXLINE 4096
 #endif
 
-#ifndef MAX
-#define MAX(a,b)    ((a) > (b) ? (a) : (b))
-#endif
-
 #define FF_LJ   0
 #define FF_EAM  1
 
diff --git a/lammps/atom.c b/lammps/atom.c
index e46d9ee..79b8f64 100644
--- a/lammps/atom.c
+++ b/lammps/atom.c
@@ -171,8 +171,9 @@ int type_str2int(const char *type) {
 int readAtom(Atom* atom, Parameter* param) {
     int len = strlen(param->input_file);
     if(strncmp(&param->input_file[len - 4], ".pdb", 4) == 0) { return readAtom_pdb(atom, param); }
+    if(strncmp(&param->input_file[len - 4], ".gro", 4) == 0) { return readAtom_gro(atom, param); }
     if(strncmp(&param->input_file[len - 4], ".dmp", 4) == 0) { return readAtom_dmp(atom, param); }
-    fprintf(stderr, "Invalid input file extension: %s\nValid choices are: pdb, dmp\n", param->input_file);
+    fprintf(stderr, "Invalid input file extension: %s\nValid choices are: pdb, gro, dmp\n", param->input_file);
     exit(-1);
     return -1;
 }
@@ -262,8 +263,86 @@ int readAtom_pdb(Atom* atom, Parameter* param) {
     return read_atoms;
 }
 
-int readAtom_dmp(Atom* atom, Parameter* param)
-{
+int readAtom_gro(Atom* atom, Parameter* param) {
+    FILE *fp = fopen(param->input_file, "r");
+    char line[MAXLINE];
+    char desc[MAXLINE];
+    int read_atoms = 0;
+    int atoms_to_read = 0;
+    int i = 0;
+
+    if(!fp) {
+        fprintf(stderr, "Could not open input file: %s\n", param->input_file);
+        exit(-1);
+        return -1;
+    }
+
+    fgets(desc, MAXLINE, fp);
+    for(i = 0; desc[i] != '\n'; i++);
+    desc[i] = '\0';
+    fgets(line, MAXLINE, fp);
+    atoms_to_read = atoi(strtok(line, " "));
+    fprintf(stdout, "System: %s with %d atoms\n", desc, atoms_to_read);
+
+    while(!feof(fp) && read_atoms < atoms_to_read) {
+        fgets(line, MAXLINE, fp);
+        char *label = strtok(line, " ");
+        int type = type_str2int(strtok(NULL, " "));
+        int atom_id = atoi(strtok(NULL, " ")) - 1;
+        atom_id = read_atoms;
+        while(atom_id + 1 >= atom->Nmax) {
+            growAtom(atom);
+        }
+
+        atom->type[atom_id] = type;
+        atom_x(atom_id) = atof(strtok(NULL, " "));
+        atom_y(atom_id) = atof(strtok(NULL, " "));
+        atom_z(atom_id) = atof(strtok(NULL, " "));
+        atom->vx[atom_id] = atof(strtok(NULL, " "));
+        atom->vy[atom_id] = atof(strtok(NULL, " "));
+        atom->vz[atom_id] = atof(strtok(NULL, " "));
+        atom->ntypes = MAX(atom->type[atom_id] + 1, atom->ntypes);
+        atom->Natoms++;
+        atom->Nlocal++;
+        read_atoms++;
+    }
+
+    if(!feof(fp)) {
+        fgets(line, MAXLINE, fp);
+        param->xlo = 0.0;
+        param->xhi = atof(strtok(line, " "));
+        param->ylo = 0.0;
+        param->yhi = atof(strtok(NULL, " "));
+        param->zlo = 0.0;
+        param->zhi = atof(strtok(NULL, " "));
+        param->xprd = param->xhi - param->xlo;
+        param->yprd = param->yhi - param->ylo;
+        param->zprd = param->zhi - param->zlo;
+    }
+
+    if(read_atoms != atoms_to_read) {
+        fprintf(stderr, "Input error: Number of atoms read do not match (%d/%d).\n", read_atoms, atoms_to_read);
+        exit(-1);
+        return -1;
+    }
+
+    atom->epsilon = allocate(ALIGNMENT, atom->ntypes * atom->ntypes * sizeof(MD_FLOAT));
+    atom->sigma6 = allocate(ALIGNMENT, atom->ntypes * atom->ntypes * sizeof(MD_FLOAT));
+    atom->cutforcesq = allocate(ALIGNMENT, atom->ntypes * atom->ntypes * sizeof(MD_FLOAT));
+    atom->cutneighsq = allocate(ALIGNMENT, atom->ntypes * atom->ntypes * sizeof(MD_FLOAT));
+    for(int i = 0; i < atom->ntypes * atom->ntypes; i++) {
+        atom->epsilon[i] = param->epsilon;
+        atom->sigma6[i] = param->sigma6;
+        atom->cutneighsq[i] = param->cutneigh * param->cutneigh;
+        atom->cutforcesq[i] = param->cutforce * param->cutforce;
+    }
+
+    fprintf(stdout, "Read %d atoms from %s\n", read_atoms, param->input_file);
+    fclose(fp);
+    return read_atoms;
+}
+
+int readAtom_dmp(Atom* atom, Parameter* param) {
     FILE *fp = fopen(param->input_file, "r");
     char line[MAXLINE];
     int natoms = 0;
@@ -355,8 +434,7 @@ int readAtom_dmp(Atom* atom, Parameter* param)
     return natoms;
 }
 
-void growAtom(Atom *atom)
-{
+void growAtom(Atom *atom) {
     int nold = atom->Nmax;
     atom->Nmax += DELTA;
 
diff --git a/lammps/includes/atom.h b/lammps/includes/atom.h
index 1f21df8..159a497 100644
--- a/lammps/includes/atom.h
+++ b/lammps/includes/atom.h
@@ -43,6 +43,7 @@ extern void initAtom(Atom*);
 extern void createAtom(Atom*, Parameter*);
 extern int readAtom(Atom*, Parameter*);
 extern int readAtom_pdb(Atom*, Parameter*);
+extern int readAtom_gro(Atom*, Parameter*);
 extern int readAtom_dmp(Atom*, Parameter*);
 extern void growAtom(Atom*);
 
diff --git a/lammps/includes/parameter.h b/lammps/includes/parameter.h
index 1da21ff..5c8eba7 100644
--- a/lammps/includes/parameter.h
+++ b/lammps/includes/parameter.h
@@ -31,9 +31,11 @@
 
 typedef struct {
     int force_field;
+    char* param_file;
     char* input_file;
     char* vtk_file;
     MD_FLOAT epsilon;
+    MD_FLOAT sigma;
     MD_FLOAT sigma6;
     MD_FLOAT temp;
     MD_FLOAT rho;
@@ -41,8 +43,10 @@ typedef struct {
     int ntypes;
     int ntimes;
     int nstat;
-    int every;
-    int halfneigh;
+    int reneigh_every;
+    int x_out_every;
+    int v_out_every;
+    int half_neigh;
     MD_FLOAT dt;
     MD_FLOAT dtforce;
     MD_FLOAT skin;
@@ -55,4 +59,9 @@ typedef struct {
     double proc_freq;
     char* eam_file;
 } Parameter;
+
+void initParameter(Parameter*);
+void readParameter(Parameter*, const char*);
+void printParameter(Parameter*);
+
 #endif
diff --git a/lammps/includes/util.h b/lammps/includes/util.h
index 20d8806..5adabe2 100644
--- a/lammps/includes/util.h
+++ b/lammps/includes/util.h
@@ -32,6 +32,15 @@
 #ifndef ABS
 #define ABS(a) ((a) >= 0 ? (a) : -(a))
 #endif
+#ifdef DEBUG
+#define DEBUG_MESSAGE   printf
+#else
+#define DEBUG_MESSAGE
+#endif
+
+#ifndef MAXLINE
+#define MAXLINE 4096
+#endif
 
 #define FF_LJ   0
 #define FF_EAM  1
diff --git a/lammps/main.c b/lammps/main.c
index bde9abb..9c7f188 100644
--- a/lammps/main.c
+++ b/lammps/main.c
@@ -49,39 +49,7 @@ extern double computeForceLJFullNeigh(Parameter*, Atom*, Neighbor*, Stats*);
 extern double computeForceLJHalfNeigh(Parameter*, Atom*, Neighbor*, Stats*);
 extern double computeForceEam(Eam*, Parameter*, Atom*, Neighbor*, Stats*);
 
-void init(Parameter *param)
-{
-    param->input_file = NULL;
-    param->vtk_file = NULL;
-    param->force_field = FF_LJ;
-    param->epsilon = 1.0;
-    param->sigma6 = 1.0;
-    param->rho = 0.8442;
-    param->ntypes = 4;
-    param->ntimes = 200;
-    param->dt = 0.005;
-    param->nx = 32;
-    param->ny = 32;
-    param->nz = 32;
-    param->cutforce = 2.5;
-    param->skin = 0.3;
-    param->cutneigh = param->cutforce + param->skin;
-    param->temp = 1.44;
-    param->nstat = 100;
-    param->mass = 1.0;
-    param->dtforce = 0.5 * param->dt;
-    param->every = 20;
-    param->proc_freq = 2.4;
-    param->halfneigh = 1;
-}
-
-double setup(
-        Parameter *param,
-        Eam *eam,
-        Atom *atom,
-        Neighbor *neighbor,
-        Stats *stats)
-{
+double setup(Parameter *param, Eam *eam, Atom *atom, Neighbor *neighbor, Stats *stats) {
     if(param->force_field == FF_EAM) { initEam(eam, param); }
     double S, E;
     param->lattice = pow((4.0 / param->rho), (1.0 / 3.0));
@@ -99,6 +67,7 @@ double setup(
     } else {
         readAtom(atom, param);
     }
+
     setupNeighbor(param);
     setupThermo(param, atom->Natoms);
     if(param->input_file == NULL) { adjustThermo(param, atom); }
@@ -106,17 +75,11 @@ double setup(
     updatePbc(atom, param);
     buildNeighbor(atom, neighbor);
     E = getTimeStamp();
-
     return E-S;
 }
 
-double reneighbour(
-        Parameter *param,
-        Atom *atom,
-        Neighbor *neighbor)
-{
+double reneighbour(Parameter *param, Atom *atom, Neighbor *neighbor) {
     double S, E;
-
     S = getTimeStamp();
     LIKWID_MARKER_START("reneighbour");
     updateAtomsPbc(atom, param);
@@ -126,12 +89,10 @@ double reneighbour(
     buildNeighbor(atom, neighbor);
     LIKWID_MARKER_STOP("reneighbour");
     E = getTimeStamp();
-
     return E-S;
 }
 
-void initialIntegrate(Parameter *param, Atom *atom)
-{
+void initialIntegrate(Parameter *param, Atom *atom) {
     MD_FLOAT* fx = atom->fx; MD_FLOAT* fy = atom->fy; MD_FLOAT* fz = atom->fz;
     MD_FLOAT* vx = atom->vx; MD_FLOAT* vy = atom->vy; MD_FLOAT* vz = atom->vz;
 
@@ -145,8 +106,7 @@ void initialIntegrate(Parameter *param, Atom *atom)
     }
 }
 
-void finalIntegrate(Parameter *param, Atom *atom)
-{
+void finalIntegrate(Parameter *param, Atom *atom) {
     MD_FLOAT* fx = atom->fx; MD_FLOAT* fy = atom->fy; MD_FLOAT* fz = atom->fz;
     MD_FLOAT* vx = atom->vx; MD_FLOAT* vy = atom->vy; MD_FLOAT* vz = atom->vz;
 
@@ -157,8 +117,7 @@ void finalIntegrate(Parameter *param, Atom *atom)
     }
 }
 
-void printAtomState(Atom *atom)
-{
+void printAtomState(Atom *atom) {
     printf("Atom counts: Natoms=%d Nlocal=%d Nghost=%d Nmax=%d\n",
             atom->Natoms, atom->Nlocal, atom->Nghost, atom->Nmax);
 
@@ -169,8 +128,7 @@ void printAtomState(Atom *atom)
     /*     } */
 }
 
-int main(int argc, char** argv)
-{
+int main(int argc, char** argv) {
     double timer[NUMTIMER];
     Eam eam;
     Atom atom;
@@ -185,10 +143,13 @@ int main(int argc, char** argv)
         //LIKWID_MARKER_REGISTER("reneighbour");
         //LIKWID_MARKER_REGISTER("pbc");
     }
-    init(&param);
 
-    for(int i = 0; i < argc; i++)
-    {
+    initParameter(&param);
+    for(int i = 0; i < argc; i++) {
+        if((strcmp(argv[i], "-p") == 0)) {
+            readParameter(&param, argv[++i]);
+            continue;
+        }
         if((strcmp(argv[i], "-f") == 0)) {
             if((param.force_field = str2ff(argv[++i])) < 0) {
                 fprintf(stderr, "Invalid force field!\n");
@@ -239,6 +200,7 @@ int main(int argc, char** argv)
         if((strcmp(argv[i], "-h") == 0) || (strcmp(argv[i], "--help") == 0)) {
             printf("MD Bench: A minimalistic re-implementation of miniMD\n");
             printf(HLINE);
+            printf("-p <string>:          file to read parameters from (can be specified more than once)\n");
             printf("-f <string>:          force field (lj or eam), default lj\n");
             printf("-i <string>:          input file with atom positions (dump)\n");
             printf("-e <string>:          input file for EAM\n");
@@ -255,6 +217,9 @@ int main(int argc, char** argv)
 
     param.cutneigh = param.cutforce + param.skin;
     setup(&param, &eam, &atom, &neighbor, &stats);
+    printParameter(&param);
+
+    printf("step\ttemp\t\tpressure\n");
     computeThermo(0, &param, &atom);
 #if defined(MEM_TRACER) || defined(INDEX_TRACER)
     traceAddresses(&param, &atom, &neighbor, n + 1);
@@ -262,7 +227,7 @@ int main(int argc, char** argv)
     if(param.force_field == FF_EAM) {
         timer[FORCE] = computeForceEam(&eam, &param, &atom, &neighbor, &stats);
     } else {
-        if( param.halfneigh ) {
+        if(param.half_neigh) {
             timer[FORCE] = computeForceLJHalfNeigh(&param, &atom, &neighbor, &stats);
         } else {
             timer[FORCE] = computeForceLJFullNeigh(&param, &atom, &neighbor, &stats);
@@ -279,7 +244,7 @@ int main(int argc, char** argv)
     for(int n = 0; n < param.ntimes; n++) {
         initialIntegrate(&param, &atom);
 
-        if((n + 1) % param.every) {
+        if((n + 1) % param.reneigh_every) {
             updatePbc(&atom, &param);
         } else {
             timer[NEIGH] += reneighbour(&param, &atom, &neighbor);
@@ -292,7 +257,7 @@ int main(int argc, char** argv)
         if(param.force_field == FF_EAM) {
             timer[FORCE] += computeForceEam(&eam, &param, &atom, &neighbor, &stats);
         } else {
-            if( param.halfneigh ) {
+            if(param.half_neigh) {
                 timer[FORCE] = computeForceLJHalfNeigh(&param, &atom, &neighbor, &stats);
             } else {
                 timer[FORCE] = computeForceLJFullNeigh(&param, &atom, &neighbor, &stats);
diff --git a/lammps/neighbor.c b/lammps/neighbor.c
index cac2b1f..a7c684e 100644
--- a/lammps/neighbor.c
+++ b/lammps/neighbor.c
@@ -70,7 +70,7 @@ void initNeighbor(Neighbor *neighbor, Parameter *param)
     neighbor->maxneighs = 100;
     neighbor->numneigh = NULL;
     neighbor->neighbors = NULL;
-    neighbor->halfneigh = param->halfneigh;
+    neighbor->halfneigh = param->half_neigh;
 }
 
 void setupNeighbor(Parameter* param)
@@ -377,17 +377,17 @@ void sortAtom(Atom* atom) {
     }
 
 #ifdef AOS
-    double* new_x = (double*) malloc(Nmax * sizeof(MD_FLOAT) * 3);
+    MD_FLOAT* new_x = (MD_FLOAT*) malloc(Nmax * sizeof(MD_FLOAT) * 3);
 #else
-    double* new_x = (double*) malloc(Nmax * sizeof(MD_FLOAT));
-    double* new_y = (double*) malloc(Nmax * sizeof(MD_FLOAT));
-    double* new_z = (double*) malloc(Nmax * sizeof(MD_FLOAT));
+    MD_FLOAT* new_x = (MD_FLOAT*) malloc(Nmax * sizeof(MD_FLOAT));
+    MD_FLOAT* new_y = (MD_FLOAT*) malloc(Nmax * sizeof(MD_FLOAT));
+    MD_FLOAT* new_z = (MD_FLOAT*) malloc(Nmax * sizeof(MD_FLOAT));
 #endif
-    double* new_vx = (double*) malloc(Nmax * sizeof(MD_FLOAT));
-    double* new_vy = (double*) malloc(Nmax * sizeof(MD_FLOAT));
-    double* new_vz = (double*) malloc(Nmax * sizeof(MD_FLOAT));
-    double* old_x = atom->x; double* old_y = atom->y; double* old_z = atom->z;
-    double* old_vx = atom->vx; double* old_vy = atom->vy; double* old_vz = atom->vz;
+    MD_FLOAT* new_vx = (MD_FLOAT*) malloc(Nmax * sizeof(MD_FLOAT));
+    MD_FLOAT* new_vy = (MD_FLOAT*) malloc(Nmax * sizeof(MD_FLOAT));
+    MD_FLOAT* new_vz = (MD_FLOAT*) malloc(Nmax * sizeof(MD_FLOAT));
+    MD_FLOAT* old_x = atom->x; MD_FLOAT* old_y = atom->y; MD_FLOAT* old_z = atom->z;
+    MD_FLOAT* old_vx = atom->vx; MD_FLOAT* old_vy = atom->vy; MD_FLOAT* old_vz = atom->vz;
 
     for(int mybin = 0; mybin<mbins; mybin++) {
         int start = mybin>0?binpos[mybin-1]:0;
diff --git a/lammps/parameter.c b/lammps/parameter.c
new file mode 100644
index 0000000..5bb9632
--- /dev/null
+++ b/lammps/parameter.c
@@ -0,0 +1,150 @@
+/*
+ * =======================================================================================
+ *
+ *   Author:   Jan Eitzinger (je), jan.eitzinger@fau.de
+ *   Copyright (c) 2020 RRZE, University Erlangen-Nuremberg
+ *
+ *   This file is part of MD-Bench.
+ *
+ *   MD-Bench is free software: you can redistribute it and/or modify it
+ *   under the terms of the GNU Lesser General Public License as published
+ *   by the Free Software Foundation, either version 3 of the License, or
+ *   (at your option) any later version.
+ *
+ *   MD-Bench is distributed in the hope that it will be useful, but WITHOUT ANY
+ *   WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
+ *   PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
+ *   details.
+ *
+ *   You should have received a copy of the GNU Lesser General Public License along
+ *   with MD-Bench.  If not, see <https://www.gnu.org/licenses/>.
+ * =======================================================================================
+ */
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+//---
+#include <parameter.h>
+#include <util.h>
+
+void initParameter(Parameter *param) {
+    param->input_file = NULL;
+    param->vtk_file = NULL;
+    param->eam_file = NULL;
+    param->force_field = FF_LJ;
+    param->epsilon = 1.0;
+    param->sigma = 1.0;
+    param->sigma6 = 1.0;
+    param->rho = 0.8442;
+    param->ntypes = 4;
+    param->ntimes = 200;
+    param->dt = 0.005;
+    param->nx = 32;
+    param->ny = 32;
+    param->nz = 32;
+    param->cutforce = 2.5;
+    param->skin = 0.3;
+    param->cutneigh = param->cutforce + param->skin;
+    param->temp = 1.44;
+    param->nstat = 100;
+    param->mass = 1.0;
+    param->dtforce = 0.5 * param->dt;
+    param->reneigh_every = 20;
+    param->x_out_every = 20;
+    param->v_out_every = 5;
+    param->half_neigh = 1;
+    param->proc_freq = 2.4;
+}
+
+void readParameter(Parameter *param, const char *filename) {
+    FILE *fp = fopen(filename, "r");
+    char line[MAXLINE];
+    int i;
+
+    if(!fp) {
+        fprintf(stderr, "Could not open parameter file: %s\n", filename);
+        exit(-1);
+    }
+
+    while(!feof(fp)) {
+        line[0] = '\0';
+        fgets(line, MAXLINE, fp);
+        for(i = 0; line[i] != '\0' && line[i] != '#'; i++);
+        line[i] = '\0';
+
+        char *tok = strtok(line, " ");
+        char *val = strtok(NULL, " ");
+
+        #define PARSE_PARAM(p,f)   if(strncmp(tok, #p, sizeof(#p) / sizeof(#p[0]) - 1) == 0) { param->p = f(val); }
+        #define PARSE_STRING(p)    PARSE_PARAM(p, strdup)
+        #define PARSE_INT(p)       PARSE_PARAM(p, atoi)
+        #define PARSE_REAL(p)      PARSE_PARAM(p, atof)
+
+        if(tok != NULL && val != NULL) {
+            PARSE_PARAM(force_field, str2ff);
+            PARSE_STRING(input_file);
+            PARSE_STRING(eam_file);
+            PARSE_STRING(vtk_file);
+            PARSE_REAL(epsilon);
+            PARSE_REAL(sigma);
+            PARSE_REAL(rho);
+            PARSE_REAL(dt);
+            PARSE_REAL(cutforce);
+            PARSE_REAL(skin);
+            PARSE_REAL(temp);
+            PARSE_REAL(mass);
+            PARSE_REAL(proc_freq);
+            PARSE_INT(ntypes);
+            PARSE_INT(ntimes);
+            PARSE_INT(nx);
+            PARSE_INT(ny);
+            PARSE_INT(nz);
+            PARSE_INT(nstat);
+            PARSE_INT(reneigh_every);
+            PARSE_INT(x_out_every);
+            PARSE_INT(v_out_every);
+            PARSE_INT(half_neigh);
+        }
+    }
+
+    // Update sigma6 parameter
+    MD_FLOAT s2 = param->sigma * param->sigma;
+    param->sigma6 = s2 * s2 * s2;
+    fclose(fp);
+}
+
+void printParameter(Parameter *param) {
+    printf("Parameters:\n");
+    if(param->input_file != NULL) {
+        printf("Input file: %s\n", param->input_file);
+    }
+
+    if(param->vtk_file != NULL) {
+        printf("VTK file: %s\n", param->vtk_file);
+    }
+
+    if(param->eam_file != NULL) {
+        printf("EAM file: %s\n", param->eam_file);
+    }
+
+    printf("\tForce field: %s\n", ff2str(param->force_field));
+    printf("\tUnit cells (nx, ny, nz): %d, %d, %d\n", param->nx, param->ny, param->nz);
+    printf("\tDomain box sizes (x, y, z): %e, %e, %e\n", param->xprd, param->yprd, param->zprd);
+    printf("\tLattice size: %e\n", param->lattice);
+    printf("\tEpsilon: %e\n", param->epsilon);
+    printf("\tSigma: %e\n", param->sigma);
+    printf("\tTemperature: %e\n", param->temp);
+    printf("\tRHO: %e\n", param->rho);
+    printf("\tMass: %e\n", param->mass);
+    printf("\tNumber of types: %d\n", param->ntypes);
+    printf("\tNumber of timesteps: %d\n", param->ntimes);
+    printf("\tReport stats every (timesteps): %d\n", param->nstat);
+    printf("\tReneighbor every (timesteps): %d\n", param->reneigh_every);
+    printf("\tOutput positions every (timesteps): %d\n", param->x_out_every);
+    printf("\tOutput velocities every (timesteps): %d\n", param->v_out_every);
+    printf("\tDelta time (dt): %e\n", param->dt);
+    printf("\tCutoff radius: %e\n", param->cutforce);
+    printf("\tSkin: %e\n", param->skin);
+    printf("\tHalf neighbor lists: %d\n", param->half_neigh);
+    printf("\tProcessor frequency (GHz): %.4f\n\n", param->proc_freq);
+}
diff --git a/lammps/thermo.c b/lammps/thermo.c
index 6f83663..73e94ae 100644
--- a/lammps/thermo.c
+++ b/lammps/thermo.c
@@ -65,8 +65,6 @@ void setupThermo(Parameter *param, int natoms)
         e_scale = 524287.985533;//16.0;
         param->dtforce /= mvv2e;
     }
-
-    printf("step\ttemp\t\tpressure\n");
 }
 
 void computeThermo(int iflag, Parameter *param, Atom *atom)