Setup first DEM example with input file from lecture

Signed-off-by: Rafael Ravedutti <rafaelravedutti@gmail.com>

data/dem/collision.in

@@ -0,0 +1,3 @@
+2
+1 0.1 0.166 0.55 0 1 0 0
+1 0.1 0.833 0.45 0 -1 0 0

data/dem/mdbench_params.conf

@@ -0,0 +1,5 @@
+ntimes 200
+dt 0.001
+reneigh_every 20
+k_s 100
+k_dn 10

lammps/atom.c

@@ -174,7 +174,8 @@ int readAtom(Atom* atom, Parameter* param) {
     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, gro, dmp\n", param->input_file);
+    if(strncmp(&param->input_file[len - 3], ".in",  3) == 0) { return readAtom_in(atom, param); }
+    fprintf(stderr, "Invalid input file extension: %s\nValid choices are: pdb, gro, dmp, in\n", param->input_file);
     exit(-1);
     return -1;
 }
@@ -435,6 +436,67 @@ int readAtom_dmp(Atom* atom, Parameter* param) {
     return natoms;
 }
 
+int readAtom_in(Atom* atom, Parameter* param) {
+    FILE *fp = fopen(param->input_file, "r");
+    char line[MAXLINE];
+    int natoms = 0;
+    int atom_id = 0;
+
+    if(!fp) {
+        fprintf(stderr, "Could not open input file: %s\n", param->input_file);
+        exit(-1);
+        return -1;
+    }
+
+    atom->ntypes = 1;
+    while(!feof(fp)) {
+        fgets(line, MAXLINE, fp);
+        natoms = atoi(line);
+        for(int i = 0; i < natoms; i++) {
+            fgets(line, MAXLINE, fp);
+
+            // TODO: store mass per atom
+            char *s_mass = strtok(line, " ");
+            if(strncmp(s_mass, "inf", 3) == 0) {
+                // Set atom's mass to INFINITY
+            } else {
+                param->mass = atof(s_mass);
+            }
+
+            atom->radius[atom_id] = atof(strtok(NULL, " "));
+            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->type[atom_id] = 0;
+            atom->ntypes = MAX(atom->type[atom_id], atom->ntypes);
+            atom_id++;
+        }
+    }
+
+    if(!natoms) {
+        fprintf(stderr, "Input error: atom data was not read!\n");
+        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", natoms, param->input_file);
+    return natoms;
+}
+
 void growAtom(Atom *atom) {
     int nold = atom->Nmax;
     atom->Nmax += DELTA;

lammps/includes/atom.h

@@ -49,6 +49,7 @@ extern int readAtom(Atom*, Parameter*);
 extern int readAtom_pdb(Atom*, Parameter*);
 extern int readAtom_gro(Atom*, Parameter*);
 extern int readAtom_dmp(Atom*, Parameter*);
+extern int readAtom_in(Atom*, Parameter*);
 extern void growAtom(Atom*);
 
 #ifdef AOS

lammps/includes/parameter.h

@@ -62,6 +62,7 @@ typedef struct {
     MD_FLOAT k_s;
     MD_FLOAT k_dn;
     MD_FLOAT gx, gy, gz;
+    MD_FLOAT reflect_x, reflect_y, reflect_z;
 } Parameter;
 
 void initParameter(Parameter*);

lammps/parameter.c

@@ -60,6 +60,9 @@ void initParameter(Parameter *param) {
     param->gx = 0.0;
     param->gy = 0.0;
     param->gz = 0.0;
+    param->reflect_x = 0.0;
+    param->reflect_y = 0.0;
+    param->reflect_z = 0.0;
 }
 
 void readParameter(Parameter *param, const char *filename) {
@@ -93,6 +96,14 @@ void readParameter(Parameter *param, const char *filename) {
             PARSE_STRING(vtk_file);
             PARSE_REAL(epsilon);
             PARSE_REAL(sigma);
+            PARSE_REAL(k_s);
+            PARSE_REAL(k_dn);
+            PARSE_REAL(reflect_x);
+            PARSE_REAL(reflect_y);
+            PARSE_REAL(reflect_z);
+            PARSE_REAL(gx);
+            PARSE_REAL(gy);
+            PARSE_REAL(gz);
             PARSE_REAL(rho);
             PARSE_REAL(dt);
             PARSE_REAL(cutforce);
@@ -139,6 +150,8 @@ void printParameter(Parameter *param) {
     printf("\tLattice size: %e\n", param->lattice);
     printf("\tEpsilon: %e\n", param->epsilon);
     printf("\tSigma: %e\n", param->sigma);
+    printf("\tSpring constant: %e\n", param->k_s);
+    printf("\tDamping constant: %e\n", param->k_dn);
     printf("\tTemperature: %e\n", param->temp);
     printf("\tRHO: %e\n", param->rho);
     printf("\tMass: %e\n", param->mass);