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

Add new setups for Copper melting with LJ and EAM

Browse Source 
Signed-off-by: Rafael Ravedutti <rafaelravedutti@gmail.com>
This commit is contained in:
Rafael Ravedutti
2021-11-30 01:33:55 +01:00
parent f7010113bf
commit bb21a885a1
14 changed files with 96831 additions and 21 deletions
Download Patch File Download Diff File Expand all files Collapse all files

102
src/atom.c
View File

@@ -24,6 +24,7 @@
*/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <atom.h>
@@ -32,6 +33,14 @@
#define DELTA 20000
#ifndef MAXLINE
#define MAXLINE 4096
#endif
#ifndef MAX
#define MAX(a,b) ((a) > (b) ? (a) : (b))
#endif
void initAtom(Atom *atom)
{
atom->x = NULL; atom->y = NULL; atom->z = NULL;
@@ -152,6 +161,99 @@ void createAtom(Atom *atom, Parameter *param)
}
}
int readAtom(Atom* atom, Parameter* param)
{
FILE *fp = fopen(param->input_file, "r");
char line[MAXLINE];
int natoms = 0;
int read_atoms = 0;
int atom_id = -1;
int ts = -1;
if(!fp) {
fprintf(stderr, "Could not open input file: %s\n", param->input_file);
exit(-1);
return -1;
}
while(!feof(fp) && ts < 1 && !read_atoms) {
fgets(line, MAXLINE, fp);
if(strncmp(line, "ITEM: ", 6) == 0) {
char *item = &line[6];
if(strncmp(item, "TIMESTEP", 8) == 0) {
fgets(line, MAXLINE, fp);
ts = atoi(line);
} else if(strncmp(item, "NUMBER OF ATOMS", 15) == 0) {
fgets(line, MAXLINE, fp);
natoms = atoi(line);
atom->Natoms = natoms;
atom->Nlocal = natoms;
while(atom->Nlocal >= atom->Nmax) {
growAtom(atom);
}
} else if(strncmp(item, "BOX BOUNDS pp pp pp", 19) == 0) {
fgets(line, MAXLINE, fp);
param->xlo = atof(strtok(line, " "));
param->xhi = atof(strtok(NULL, " "));
param->xprd = param->xhi - param->xlo;
fgets(line, MAXLINE, fp);
param->ylo = atof(strtok(line, " "));
param->yhi = atof(strtok(NULL, " "));
param->yprd = param->yhi - param->ylo;
fgets(line, MAXLINE, fp);
param->zlo = atof(strtok(line, " "));
param->zhi = atof(strtok(NULL, " "));
param->zprd = param->zhi - param->zlo;
} else if(strncmp(item, "ATOMS id type x y z vx vy vz", 28) == 0) {
for(int i = 0; i < natoms; i++) {
fgets(line, MAXLINE, fp);
atom_id = atoi(strtok(line, " ")) - 1;
atom->type[atom_id] = atoi(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->ntypes = MAX(atom->type[atom_id], atom->ntypes);
read_atoms++;
}
} else {
fprintf(stderr, "Invalid item: %s\n", item);
exit(-1);
return -1;
}
} else {
fprintf(stderr, "Invalid input from file, expected item reference but got:\n%s\n", line);
exit(-1);
return -1;
}
}
if(ts < 0 || !natoms || !read_atoms) {
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;

4
src/eam_utils.c
View File

@@ -43,7 +43,7 @@ void initEam(Eam* eam, Parameter* param) {
}
void coeff(Eam* eam, Parameter* param) {
read_file(&eam->file, param->input_file);
read_eam_file(&eam->file, param->eam_file);
param->mass = eam->file.mass;
param->cutforce = eam->file.cut;
param->cutneigh = param->cutforce + 1.0;
@@ -59,7 +59,7 @@ void init_style(Eam* eam, Parameter* param) {
array2spline(eam, param);
}
void read_file(Funcfl* file, const char* filename) {
void read_eam_file(Funcfl* file, const char* filename) {
FILE* fptr;
char line[MAXLINE];

1
src/includes/atom.h
View File

@@ -41,6 +41,7 @@ typedef struct {
extern void initAtom(Atom*);
extern void createAtom(Atom*, Parameter*);
extern int readAtom(Atom*, Parameter*);
extern void growAtom(Atom*);
#ifdef AOS

2
src/includes/eam.h
View File

@@ -47,7 +47,7 @@ typedef struct {
void initEam(Eam* eam, Parameter* param);
void coeff(Eam* eam, Parameter* param);
void init_style(Eam* eam, Parameter *param);
void read_file(Funcfl* file, const char* filename);
void read_eam_file(Funcfl* file, const char* filename);
void file2array(Eam* eam);
void array2spline(Eam* eam, Parameter* param);
void interpolate(int n, MD_FLOAT delta, MD_FLOAT* f, MD_FLOAT* spline);

2
src/includes/neighbor.h
View File

@@ -34,7 +34,7 @@ typedef struct {
} Neighbor;
extern void initNeighbor(Neighbor*, Parameter*);
extern void setupNeighbor();
extern void setupNeighbor(Parameter*);
extern void binatoms(Atom*);
extern void buildNeighbor(Atom*, Neighbor*);
extern void sortAtom(Atom*);

2
src/includes/parameter.h
View File

@@ -51,7 +51,9 @@ typedef struct {
MD_FLOAT cutneigh;
int nx, ny, nz;
MD_FLOAT lattice;
MD_FLOAT xlo, xhi, ylo, yhi, zlo, zhi;
MD_FLOAT xprd, yprd, zprd;
double proc_freq;
char* eam_file;
} Parameter;
#endif

20
src/main.c
View File

@@ -88,13 +88,17 @@ double setup(
S = getTimeStamp();
initAtom(atom);
initNeighbor(neighbor, param);
initPbc(atom);
initStats(stats);
setupNeighbor();
createAtom(atom, param);
initNeighbor(neighbor, param);
if(param->input_file == NULL) {
createAtom(atom, param);
} else {
readAtom(atom, param);
}
setupNeighbor(param);
setupThermo(param, atom->Natoms);
adjustThermo(param, atom);
if(param->input_file == NULL) { adjustThermo(param, atom); }
setupPbc(atom, param);
updatePbc(atom, param);
buildNeighbor(atom, neighbor);
@@ -209,6 +213,11 @@ int main(int argc, char** argv)
param.input_file = strdup(argv[++i]);
continue;
}
if((strcmp(argv[i], "-e") == 0))
{
param.eam_file = strdup(argv[++i]);
continue;
}
if((strcmp(argv[i], "-n") == 0) || (strcmp(argv[i], "--nsteps") == 0))
{
param.ntimes = atoi(argv[++i]);
@@ -244,7 +253,8 @@ int main(int argc, char** argv)
printf("MD Bench: A minimalistic re-implementation of miniMD\n");
printf(HLINE);
printf("-f <string>: force field (lj or eam), default lj\n");
printf("-i <string>: input file for EAM\n");
printf("-i <string>: input file with atom positions (dump)\n");
printf("-e <string>: input file for EAM\n");
printf("-n / --nsteps <int>: set number of timesteps for simulation\n");
printf("-nx/-ny/-nz <int>: set linear dimension of systembox in x/y/z direction\n");
printf("--freq <real>: processor frequency (GHz)\n");

38
src/neighbor.c
View File

@@ -71,22 +71,46 @@ void initNeighbor(Neighbor *neighbor, Parameter *param)
neighbor->neighbors = NULL;
}
void setupNeighbor()
void setupNeighbor(Parameter* param)
{
MD_FLOAT coord;
int mbinxhi, mbinyhi, mbinzhi;
int nextx, nexty, nextz;
if(param->input_file != NULL) {
xprd = param->xprd;
yprd = param->yprd;
zprd = param->zprd;
}
// TODO: update lo and hi for standard case and use them here instead
MD_FLOAT xlo = 0.0; MD_FLOAT xhi = xprd;
MD_FLOAT ylo = 0.0; MD_FLOAT yhi = yprd;
MD_FLOAT zlo = 0.0; MD_FLOAT zhi = zprd;
cutneighsq = cutneigh * cutneigh;
binsizex = xprd / nbinx;
binsizey = yprd / nbiny;
binsizez = zprd / nbinz;
bininvx = 1.0 / binsizex;
bininvy = 1.0 / binsizey;
bininvz = 1.0 / binsizez;
if(param->input_file != NULL) {
binsizex = cutneigh * 0.5;
binsizey = cutneigh * 0.5;
binsizez = cutneigh * 0.5;
nbinx = (int)((param->xhi - param->xlo) / binsizex);
nbiny = (int)((param->yhi - param->ylo) / binsizey);
nbinz = (int)((param->zhi - param->zlo) / binsizez);
if(nbinx == 0) { nbinx = 1; }
if(nbiny == 0) { nbiny = 1; }
if(nbinz == 0) { nbinz = 1; }
bininvx = nbinx / (param->xhi - param->xlo);
bininvy = nbiny / (param->yhi - param->ylo);
bininvz = nbinz / (param->zhi - param->zlo);
} else {
binsizex = xprd / nbinx;
binsizey = yprd / nbiny;
binsizez = zprd / nbinz;
bininvx = 1.0 / binsizex;
bininvy = 1.0 / binsizey;
bininvz = 1.0 / binsizez;
}
coord = xlo - cutneigh - SMALL * xprd;
mbinxlo = (int) (coord * bininvx);

47
src/util.c
View File

@@ -30,13 +30,50 @@
#define IR 2836
#define MASK 123459876
double myrandom(int* idum)
double myrandom(int* seed)
{
int k= (*idum) / IQ;
int k= (*seed) / IQ;
double ans;
*idum = IA * (*idum - k * IQ) - IR * k;
if(*idum < 0) *idum += IM;
ans = AM * (*idum);
*seed = IA * (*seed - k * IQ) - IR * k;
if(*seed < 0) *seed += IM;
ans = AM * (*seed);
return ans;
}
void random_reset(int *seed, int ibase, double *coord)
{
int i;
char *str = (char *) &ibase;
int n = sizeof(int);
unsigned int hash = 0;
for (i = 0; i < n; i++) {
hash += str[i];
hash += (hash << 10);
hash ^= (hash >> 6);
}
str = (char *) coord;
n = 3 * sizeof(double);
for (i = 0; i < n; i++) {
hash += str[i];
hash += (hash << 10);
hash ^= (hash >> 6);
}
hash += (hash << 3);
hash ^= (hash >> 11);
hash += (hash << 15);
// keep 31 bits of unsigned int as new seed
// do not allow seed = 0, since will cause hang in gaussian()
*seed = hash & 0x7ffffff;
if (!(*seed)) *seed = 1;
// warm up the RNG
for (i = 0; i < 5; i++) myrandom(seed);
//save = 0;
}