moebiusband/MD-Bench
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Add parameter reading for LAMMPS variant

Browse Source 
Signed-off-by: Rafael Ravedutti <rafaelravedutti@gmail.com>
This commit is contained in:
Rafael Ravedutti 2022-03-17 02:44:34 +01:00
parent d4b34e1fa4
commit 887f41871c
11 changed files with 293 additions and 80 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
gromacs/includes/parameter.h
View File

@ -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;

9
gromacs/includes/simd/avx512_float.h
View File

@ -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;

4
gromacs/includes/util.h
View File

@ -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

88
lammps/atom.c
View File

@ -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;

1
lammps/includes/atom.h
View File

@ -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*);

13
lammps/includes/parameter.h
View File

@ -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

9
lammps/includes/util.h
View File

@ -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

75
lammps/main.c
View File

@ -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);

20
lammps/neighbor.c
View File

@ -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;

150
lammps/parameter.c Normal file
View File

@ -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);
}

2
lammps/thermo.c
View File

@ -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)