From 32836eebcb36eda06dadd55d2d6ed21ffed7437a Mon Sep 17 00:00:00 2001 From: Rafael Ravedutti Date: Thu, 7 Jul 2022 02:11:50 +0200 Subject: [PATCH] Setup first DEM example with input file from lecture Signed-off-by: Rafael Ravedutti --- data/dem/collision.in | 3 ++ data/dem/mdbench_params.conf | 5 +++ lammps/atom.c | 64 +++++++++++++++++++++++++++++++++++- lammps/includes/atom.h | 1 + lammps/includes/parameter.h | 1 + lammps/parameter.c | 13 ++++++++ 6 files changed, 86 insertions(+), 1 deletion(-) create mode 100644 data/dem/collision.in create mode 100644 data/dem/mdbench_params.conf diff --git a/data/dem/collision.in b/data/dem/collision.in new file mode 100644 index 0000000..0dedb0a --- /dev/null +++ b/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 diff --git a/data/dem/mdbench_params.conf b/data/dem/mdbench_params.conf new file mode 100644 index 0000000..622c19f --- /dev/null +++ b/data/dem/mdbench_params.conf @@ -0,0 +1,5 @@ +ntimes 200 +dt 0.001 +reneigh_every 20 +k_s 100 +k_dn 10 diff --git a/lammps/atom.c b/lammps/atom.c index 46e4110..27199fb 100644 --- a/lammps/atom.c +++ b/lammps/atom.c @@ -174,7 +174,8 @@ int readAtom(Atom* atom, Parameter* param) { if(strncmp(¶m->input_file[len - 4], ".pdb", 4) == 0) { return readAtom_pdb(atom, param); } if(strncmp(¶m->input_file[len - 4], ".gro", 4) == 0) { return readAtom_gro(atom, param); } if(strncmp(¶m->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(¶m->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; diff --git a/lammps/includes/atom.h b/lammps/includes/atom.h index 4c473d8..8c9a534 100644 --- a/lammps/includes/atom.h +++ b/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 diff --git a/lammps/includes/parameter.h b/lammps/includes/parameter.h index 497f4ce..30c03bd 100644 --- a/lammps/includes/parameter.h +++ b/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*); diff --git a/lammps/parameter.c b/lammps/parameter.c index e16f538..51c0ade 100644 --- a/lammps/parameter.c +++ b/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);