From e7869286d7fd915057312a34305ac1596fa78c4f Mon Sep 17 00:00:00 2001 From: Jan Eitzinger Date: Tue, 18 Aug 2020 14:27:28 +0200 Subject: [PATCH] Introduce modular version. --- include_CLANG.mk | 4 +- src/allocate.c | 20 +- src/atom.c | 186 +++++++ src/includes/allocate.h | 4 +- src/includes/atom.h | 41 ++ src/includes/neighbor.h | 45 ++ src/includes/parameter.h | 46 ++ src/includes/pbc.h | 36 ++ src/includes/thermo.h | 35 ++ src/includes/util.h | 41 ++ src/main.c | 929 +++-------------------------------- src/neighbor.c | 345 +++++++++++++ src/pbc.c | 161 ++++++ src/thermo.c | 132 +++++ src/util.c | 46 ++ util/mdBench.c | 1000 ++++++++++++++++++++++++++++++++++++++ 16 files changed, 2208 insertions(+), 863 deletions(-) create mode 100644 src/atom.c create mode 100644 src/includes/atom.h create mode 100644 src/includes/neighbor.h create mode 100644 src/includes/parameter.h create mode 100644 src/includes/pbc.h create mode 100644 src/includes/thermo.h create mode 100644 src/includes/util.h create mode 100644 src/neighbor.c create mode 100644 src/pbc.c create mode 100644 src/thermo.c create mode 100644 src/util.c create mode 100644 util/mdBench.c diff --git a/include_CLANG.mk b/include_CLANG.mk index c5faaac..d2992b7 100644 --- a/include_CLANG.mk +++ b/include_CLANG.mk @@ -7,11 +7,11 @@ ANSI_CFLAGS += -std=c99 ANSI_CFLAGS += -pedantic ANSI_CFLAGS += -Wextra -CFLAGS = -O2 -g #$(ANSI_CFLAGS) +CFLAGS = -O3 #-g $(ANSI_CFLAGS) ASFLAGS = -masm=intel CXXFLAGS = $(CFLAGS) FCFLAGS = LFLAGS = -DEFINES = -D_GNU_SOURCE +DEFINES = -D_GNU_SOURCE -DALIGNMENT=64 INCLUDES = LIBS = diff --git a/src/allocate.c b/src/allocate.c index 2ad9236..dca0c97 100644 --- a/src/allocate.c +++ b/src/allocate.c @@ -2,7 +2,7 @@ * ======================================================================================= * * Author: Jan Eitzinger (je), jan.eitzinger@fau.de - * Copyright (c) 2019 RRZE, University Erlangen-Nuremberg + * Copyright (c) 2020 RRZE, University Erlangen-Nuremberg * * Permission is hereby granted, free of charge, to any person obtaining a copy * of this software and associated documentation files (the "Software"), to deal @@ -27,6 +27,7 @@ #include #include +#include #include void* allocate (int alignment, size_t bytesize) @@ -56,3 +57,20 @@ void* allocate (int alignment, size_t bytesize) return ptr; } + +void* reallocate ( + void* ptr, + int alignment, + size_t newBytesize, + size_t oldBytesize) +{ + void* newarray; + newarray = allocate(alignment, newBytesize); + + if(ptr != NULL) { + memcpy(newarray, ptr, oldBytesize); + free(ptr); + } + + return newarray; +} diff --git a/src/atom.c b/src/atom.c new file mode 100644 index 0000000..b61060c --- /dev/null +++ b/src/atom.c @@ -0,0 +1,186 @@ +/* + * ======================================================================================= + * + * Author: Jan Eitzinger (je), jan.eitzinger@fau.de + * Copyright (c) 2020 RRZE, University Erlangen-Nuremberg + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in all + * copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + * + * ======================================================================================= + */ +#include + +#include +#include +#include +#include + +#define DELTA 20000 + +void createAtom(Atom *atom, Parameter *param) +{ + double xlo = 0.0; double xhi = param->xprd; + double ylo = 0.0; double yhi = param->yprd; + double zlo = 0.0; double zhi = param->zprd; + atom->Natoms = 4 * param->nx * param->ny * param->nz; + atom->Nlocal = 0; + double alat = pow((4.0 / param->rho), (1.0 / 3.0)); + int ilo = (int) (xlo / (0.5 * alat) - 1); + int ihi = (int) (xhi / (0.5 * alat) + 1); + int jlo = (int) (ylo / (0.5 * alat) - 1); + int jhi = (int) (yhi / (0.5 * alat) + 1); + int klo = (int) (zlo / (0.5 * alat) - 1); + int khi = (int) (zhi / (0.5 * alat) + 1); + + ilo = MAX(ilo, 0); + ihi = MIN(ihi, 2 * param->nx - 1); + jlo = MAX(jlo, 0); + jhi = MIN(jhi, 2 * param->ny - 1); + klo = MAX(klo, 0); + khi = MIN(khi, 2 * param->nz - 1); + + double xtmp, ytmp, ztmp, vxtmp, vytmp, vztmp; + int i, j, k, m, n; + int sx = 0; int sy = 0; int sz = 0; + int ox = 0; int oy = 0; int oz = 0; + int subboxdim = 8; + + while(oz * subboxdim <= khi) { + + k = oz * subboxdim + sz; + j = oy * subboxdim + sy; + i = ox * subboxdim + sx; + + if(((i + j + k) % 2 == 0) && + (i >= ilo) && (i <= ihi) && + (j >= jlo) && (j <= jhi) && + (k >= klo) && (k <= khi)) { + + xtmp = 0.5 * alat * i; + ytmp = 0.5 * alat * j; + ztmp = 0.5 * alat * k; + + if( xtmp >= xlo && xtmp < xhi && + ytmp >= ylo && ytmp < yhi && + ztmp >= zlo && ztmp < zhi ) { + + n = k * (2 * param->ny) * (2 * param->nx) + + j * (2 * param->nx) + + i + 1; + + for(m = 0; m < 5; m++) { + myrandom(&n); + } + vxtmp = myrandom(&n); + + for(m = 0; m < 5; m++){ + myrandom(&n); + } + vytmp = myrandom(&n); + + for(m = 0; m < 5; m++) { + myrandom(&n); + } + vztmp = myrandom(&n); + + if(atom->Nlocal == atom->Nmax) { + growAtom(atom); + } + + atom->x[atom->Nlocal] = xtmp; + atom->y[atom->Nlocal] = ytmp; + atom->z[atom->Nlocal] = ztmp; + atom->vx[atom->Nlocal] = vxtmp; + atom->vy[atom->Nlocal] = vytmp; + atom->vz[atom->Nlocal] = vztmp; + atom->Nlocal++; + } + } + + sx++; + + if(sx == subboxdim) { sx = 0; sy++; } + if(sy == subboxdim) { sy = 0; sz++; } + if(sz == subboxdim) { sz = 0; ox++; } + if(ox * subboxdim > ihi) { ox = 0; oy++; } + if(oy * subboxdim > jhi) { oy = 0; oz++; } + } +} + +void growAtom(Atom *atom) +{ + int nold = atom->Nmax; + atom->Nmax += DELTA; + + atom->x = (double*) reallocate(atom->x, ALIGNMENT, atom->Nmax * sizeof(double), nold * sizeof(double)); + atom->y = (double*) reallocate(atom->y, ALIGNMENT, atom->Nmax * sizeof(double), nold * sizeof(double)); + atom->z = (double*) reallocate(atom->z, ALIGNMENT, atom->Nmax * sizeof(double), nold * sizeof(double)); + atom->vx = (double*) reallocate(atom->vx, ALIGNMENT, atom->Nmax * sizeof(double), nold * sizeof(double)); + atom->vy = (double*) reallocate(atom->vy, ALIGNMENT, atom->Nmax * sizeof(double), nold * sizeof(double)); + atom->vz = (double*) reallocate(atom->vz, ALIGNMENT, atom->Nmax * sizeof(double), nold * sizeof(double)); + atom->fx = (double*) reallocate(atom->fx, ALIGNMENT, atom->Nmax * sizeof(double), nold * sizeof(double)); + atom->fy = (double*) reallocate(atom->fy, ALIGNMENT, atom->Nmax * sizeof(double), nold * sizeof(double)); + atom->fz = (double*) reallocate(atom->fz, ALIGNMENT, atom->Nmax * sizeof(double), nold * sizeof(double)); +} + + +/* void sortAtom() */ +/* { */ +/* binatoms(neighbor); */ +/* int* binpos = neighbor->bincount; */ +/* int* bins = neighbor->bins; */ + +/* int mbins = neighbor->mbins; */ +/* int atoms_per_bin = neighbor->atoms_per_bin; */ + +/* for(int i=1; i0?binpos[mybin-1]:0; */ +/* int count = binpos[mybin] - start; */ + +/* for(int k=0; k #ifndef __ALLOCATE_H_ #define __ALLOCATE_H_ - extern void* allocate (int alignment, size_t bytesize); - +extern void* reallocate (void* ptr, int alignment, size_t newBytesize, size_t oldBytesize); #endif diff --git a/src/includes/atom.h b/src/includes/atom.h new file mode 100644 index 0000000..168e556 --- /dev/null +++ b/src/includes/atom.h @@ -0,0 +1,41 @@ +/* + * ======================================================================================= + * + * Author: Jan Eitzinger (je), jan.eitzinger@fau.de + * Copyright (c) 2020 RRZE, University Erlangen-Nuremberg + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in all + * copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + * + * ======================================================================================= + */ +#include + +#ifndef __ATOM_H_ +#define __ATOM_H_ + +typedef struct { + int Natoms, Nlocal, Nghost, Nmax; + double *x, *y, *z; + double *vx, *vy, *vz; + double *fx, *fy, *fz; +} Atom; + +extern void createAtom(Atom*, Parameter*); +extern void growAtom(Atom*); +#endif diff --git a/src/includes/neighbor.h b/src/includes/neighbor.h new file mode 100644 index 0000000..8dc1016 --- /dev/null +++ b/src/includes/neighbor.h @@ -0,0 +1,45 @@ +/* + * ======================================================================================= + * + * Author: Jan Eitzinger (je), jan.eitzinger@fau.de + * Copyright (c) 2020 RRZE, University Erlangen-Nuremberg + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in all + * copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + * + * ======================================================================================= + */ +#include +#include + +#ifndef __NEIGHBOR_H_ +#define __NEIGHBOR_H_ +typedef struct { + /* double cutneigh; // neighbor cutoff */ + int every; + int ncalls; + int* neighbors; + int maxneighs; + int* numneigh; +} Neighbor; + +extern void initNeighbor(Neighbor*, Parameter*); +extern void setupNeighbor(Parameter*); +extern void binatoms(Atom*); +extern void buildNeighbor(Atom*, Neighbor*); +#endif diff --git a/src/includes/parameter.h b/src/includes/parameter.h new file mode 100644 index 0000000..152f2ed --- /dev/null +++ b/src/includes/parameter.h @@ -0,0 +1,46 @@ +/* + * ======================================================================================= + * + * Author: Jan Eitzinger (je), jan.eitzinger@fau.de + * Copyright (c) 2020 RRZE, University Erlangen-Nuremberg + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in all + * copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + * + * ======================================================================================= + */ +#ifndef __PARAMETER_H_ +#define __PARAMETER_H_ + +typedef struct { + double epsilon; + double sigma6; + double temp; + double rho; + double mass; + int ntimes; + int nstat; + int every; + double dt; + double dtforce; + double cutforce; + double cutneigh; + int nx, ny, nz; + double xprd, yprd, zprd; +} Parameter; +#endif diff --git a/src/includes/pbc.h b/src/includes/pbc.h new file mode 100644 index 0000000..6bbccd8 --- /dev/null +++ b/src/includes/pbc.h @@ -0,0 +1,36 @@ +/* + * ======================================================================================= + * + * Author: Jan Eitzinger (je), jan.eitzinger@fau.de + * Copyright (c) 2020 RRZE, University Erlangen-Nuremberg + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in all + * copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + * + * ======================================================================================= + */ +#include +#include + +#ifndef __PBC_H_ +#define __PBC_H_ +extern void initPbc(); +extern void updatePbc(Atom*, Parameter*); +extern void updateAtomsPbc(Atom*, Parameter*); +extern void setupPbc(Atom*, Parameter*); +#endif diff --git a/src/includes/thermo.h b/src/includes/thermo.h new file mode 100644 index 0000000..9321969 --- /dev/null +++ b/src/includes/thermo.h @@ -0,0 +1,35 @@ +/* + * ======================================================================================= + * + * Author: Jan Eitzinger (je), jan.eitzinger@fau.de + * Copyright (c) 2020 RRZE, University Erlangen-Nuremberg + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in all + * copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + * + * ======================================================================================= + */ +#include +#include + +#ifndef __THERMO_H_ +#define __THERMO_H_ +extern void setupThermo(Parameter*, int); +extern void computeThermo(int, Parameter*, Atom*); +extern void adjustThermo(Parameter*, Atom*); +#endif diff --git a/src/includes/util.h b/src/includes/util.h new file mode 100644 index 0000000..91ef64d --- /dev/null +++ b/src/includes/util.h @@ -0,0 +1,41 @@ +/* + * ======================================================================================= + * + * Author: Jan Eitzinger (je), jan.eitzinger@fau.de + * Copyright (c) 2020 RRZE, University Erlangen-Nuremberg + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in all + * copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + * + * ======================================================================================= + */ +#ifndef __UTIL_H_ +#define __UTIL_H_ + +#ifndef MIN +#define MIN(x,y) ((x)<(y)?(x):(y)) +#endif +#ifndef MAX +#define MAX(x,y) ((x)>(y)?(x):(y)) +#endif +#ifndef ABS +#define ABS(a) ((a) >= 0 ? (a) : -(a)) +#endif + +extern double myrandom(int*); +#endif diff --git a/src/main.c b/src/main.c index 8138ef6..841bf4f 100644 --- a/src/main.c +++ b/src/main.c @@ -2,7 +2,7 @@ * ======================================================================================= * * Author: Jan Eitzinger (je), jan.eitzinger@fau.de - * Copyright (c) 2019 RRZE, University Erlangen-Nuremberg + * Copyright (c) 2020 RRZE, University Erlangen-Nuremberg * * Permission is hereby granted, free of charge, to any person obtaining a copy * of this software and associated documentation files (the "Software"), to deal @@ -33,290 +33,24 @@ #include #include +#include +#include +#include +#include +#include +#include + #define HLINE "----------------------------------------------------------------------------\n" #define FACTOR 0.999 #define SMALL 1.0e-6 #define DELTA 20000 -#ifndef MIN -#define MIN(x,y) ((x)<(y)?(x):(y)) -#endif -#ifndef MAX -#define MAX(x,y) ((x)>(y)?(x):(y)) -#endif -#ifndef ABS -#define ABS(a) ((a) >= 0 ? (a) : -(a)) -#endif - -static int Natoms, Nlocal, Nghost, Nmax; -static double Cutneigh; // neighbor cutoff -static double xprd, yprd, zprd; -static double xlo, xhi; -static double ylo, yhi; -static double zlo, zhi; -static double *x, *y, *z; -static double *vx, *vy, *vz; -static double *fx, *fy, *fz; - -static int NmaxGhost; -static int *BorderMap; -static int *PBCx, *PBCy, *PBCz; - -typedef struct { - int* numneigh; - int* neighbors; - int maxneighs; - int nbinx, nbiny, nbinz; - /* double cutneigh; // neighbor cutoff */ - double cutneighsq; // neighbor cutoff squared - int every; - int ncalls; - int max_totalneigh; - int *bincount; - int *bins; - int nmax; - int nstencil; // # of bins in stencil - int* stencil; // stencil list of bin offsets - int mbins; //total number of bins - int atoms_per_bin; // max atoms per bin - int mbinx, mbiny, mbinz; // n bins in x, y, z - int mbinxlo, mbinylo, mbinzlo; - double binsizex, binsizey, binsizez; - double bininvx, bininvy, bininvz; -} Neighbor; - -typedef struct { - double epsilon; - double sigma6; - double temp; - double rho; - double mass; - int ntimes; - int nstat; - double dt; - double dtforce; - double cutforce; - int nx, ny, nz; -} Parameter; - -typedef struct { - int *steparr; - double *tmparr; - double *engarr; - double *prsarr; - double mvv2e; - int dof_boltz; - double t_scale; - double p_scale; - double e_scale; - double t_act; - double p_act; - double e_act; - int mstat; -} Thermo; - -/* Park/Miller RNG w/out MASKING, so as to be like f90s version */ -#define IA 16807 -#define IM 2147483647 -#define AM (1.0/IM) -#define IQ 127773 -#define IR 2836 -#define MASK 123459876 - -double myrandom(int* idum) +void init(Atom *atom, Parameter *param) { - int k= (*idum) / IQ; - double ans; - - *idum = IA * (*idum - k * IQ) - IR * k; - if(*idum < 0) *idum += IM; - ans = AM * (*idum); - return ans; -} - -int coord2bin(Neighbor* neighbor, double xin, double yin, double zin) -{ - int ix, iy, iz; - double bininvx = neighbor->bininvx; - double bininvy = neighbor->bininvy; - double bininvz = neighbor->bininvz; - int mbinxlo = neighbor->mbinxlo; - int mbinylo = neighbor->mbinylo; - int mbinzlo = neighbor->mbinzlo; - - if(xin >= xprd) { - ix = (int)((xin - xprd) * bininvx) + neighbor->nbinx - mbinxlo; - } else if(xin >= 0.0) { - ix = (int)(xin * bininvx) - mbinxlo; - } else { - ix = (int)(xin * bininvx) - mbinxlo - 1; - } - - if(yin >= yprd) { - iy = (int)((yin - yprd) * bininvy) + neighbor->nbiny - mbinylo; - } else if(yin >= 0.0) { - iy = (int)(yin * bininvy) - mbinylo; - } else { - iy = (int)(yin * bininvy) - mbinylo - 1; - } - - if(zin >= zprd) { - iz = (int)((zin - zprd) * bininvz) + neighbor->nbinz - mbinzlo; - } else if(zin >= 0.0) { - iz = (int)(zin * bininvz) - mbinzlo; - } else { - iz = (int)(zin * bininvz) - mbinzlo - 1; - } - - return (iz * neighbor->mbiny * neighbor->mbinx + iy * neighbor->mbinx + ix + 1); -} - -void binatoms(Neighbor *neighbor) -{ - int* bincount = neighbor->bincount; - int mbins = neighbor->mbins; - int nall = Nlocal + Nghost; - int resize = 1; - - while(resize > 0) { - resize = 0; - - for(int i = 0; i < mbins; i++) { - bincount[i] = 0; - } - - for(int i = 0; i < nall; i++) { - int ibin = coord2bin(neighbor, x[i], y[i], z[i]); - - if(bincount[ibin] < neighbor->atoms_per_bin) { - int ac = neighbor->bincount[ibin]++; - neighbor->bins[ibin * neighbor->atoms_per_bin + ac] = i; - } else { - resize = 1; - } - } - - if(resize) { - free(neighbor->bins); - neighbor->atoms_per_bin *= 2; - neighbor->bins = (int*) malloc(mbins * neighbor->atoms_per_bin * sizeof(int)); - } - } -} - -double bindist(Neighbor *neighbor, int i, int j, int k) -{ - double delx, dely, delz; - - if(i > 0) { - delx = (i - 1) * neighbor->binsizex; - } else if(i == 0) { - delx = 0.0; - } else { - delx = (i + 1) * neighbor->binsizex; - } - - if(j > 0) { - dely = (j - 1) * neighbor->binsizey; - } else if(j == 0) { - dely = 0.0; - } else { - dely = (j + 1) * neighbor->binsizey; - } - - if(k > 0) { - delz = (k - 1) * neighbor->binsizez; - } else if(k == 0) { - delz = 0.0; - } else { - delz = (k + 1) * neighbor->binsizez; - } - - return (delx * delx + dely * dely + delz * delz); -} - -void buildNeighborlist(Neighbor *neighbor) -{ - neighbor->ncalls++; - int nall = Nlocal + Nghost; - - /* extend atom arrays if necessary */ - if(nall > neighbor->nmax) { - neighbor->nmax = nall; - if(neighbor->numneigh) free(neighbor->numneigh); - if(neighbor->neighbors) free(neighbor->neighbors); - neighbor->numneigh = (int*) malloc(neighbor->nmax * sizeof(int)); - neighbor->neighbors = (int*) malloc(neighbor->nmax * neighbor->maxneighs * sizeof(int*)); - } - - /* bin local & ghost atoms */ - binatoms(neighbor); - int resize = 1; - - /* loop over each atom, storing neighbors */ - while(resize) { - int new_maxneighs = neighbor->maxneighs; - resize = 0; - - for(int i = 0; i < Nlocal; i++) { - int* neighptr = &neighbor->neighbors[i * neighbor->maxneighs]; - int n = 0; - double xtmp = x[i]; - double ytmp = y[i]; - double ztmp = z[i]; - int ibin = coord2bin(neighbor, xtmp, ytmp, ztmp); - - for(int k = 0; k < neighbor->nstencil; k++) { - int jbin = ibin + neighbor->stencil[k]; - int* loc_bin = &neighbor->bins[jbin * neighbor->atoms_per_bin]; - - for(int m = 0; m < neighbor->bincount[jbin]; m++) { - int j = loc_bin[m]; - - if ( j == i ){ - continue; - } - - double delx = xtmp - x[j]; - double dely = ytmp - y[j]; - double delz = ztmp - z[j]; - double rsq = delx * delx + dely * dely + delz * delz; - - if( rsq <= neighbor->cutneighsq ) { - neighptr[n++] = j; - } - } - } - - neighbor->numneigh[i] = n; - - if(n >= neighbor->maxneighs) { - resize = 1; - - if(n >= new_maxneighs) { - new_maxneighs = n; - } - } - } - - if(resize) { - neighbor->maxneighs = new_maxneighs * 1.2; - free(neighbor->neighbors); - neighbor->neighbors = (int*) malloc(Nmax* neighbor->maxneighs * sizeof(int)); - } - } -} - -void init(Neighbor *neighbor, Parameter *param) -{ - x = NULL; y = NULL; z = NULL; - vx = NULL; vy = NULL; vz = NULL; - fx = NULL; fy = NULL; fz = NULL; - - NmaxGhost = 0; - BorderMap = NULL; - PBCx = NULL; PBCy = NULL; PBCz = NULL; + atom->x = NULL; atom->y = NULL; atom->z = NULL; + atom->vx = NULL; atom->vy = NULL; atom->vz = NULL; + atom->fx = NULL; atom->fy = NULL; atom->fz = NULL; param->epsilon = 1.0; param->sigma6 = 1.0; @@ -327,569 +61,25 @@ void init(Neighbor *neighbor, Parameter *param) param->ny = 32; param->nz = 32; param->cutforce = 2.5; + param->cutneigh = param->cutforce + 0.30; param->temp = 1.44; param->nstat = 100; param->mass = 1.0; param->dtforce = 0.5 * param->dt; - - Cutneigh = param->cutforce + 0.30; - double neighscale = 5.0 / 6.0; - neighbor->nbinx = neighscale * param->nx; - neighbor->nbiny = neighscale * param->ny; - neighbor->nbinz = neighscale * param->nz; - neighbor->every = 20; - neighbor->ncalls = 0; - neighbor->nmax = 0; - neighbor->atoms_per_bin = 8; - neighbor->maxneighs = 100; - /* neighbor->cutneigh = param->cutforce + 0.30; */ - neighbor->numneigh = NULL; - neighbor->neighbors = NULL; - neighbor->stencil = NULL; - neighbor->bins = NULL; - neighbor->bincount = NULL; -} - -void setup(Neighbor *neighbor, Parameter *param) -{ + param->every = 20; double lattice = pow((4.0 / param->rho), (1.0 / 3.0)); - double coord; - int mbinxhi, mbinyhi, mbinzhi; - int nextx, nexty, nextz; - - xprd = param->nx * lattice; - yprd = param->ny * lattice; - zprd = param->nz * lattice; - - xlo = 0.0; xhi = xprd; - ylo = 0.0; yhi = yprd; - zlo = 0.0; zhi = zprd; - - neighbor->cutneighsq = Cutneigh * Cutneigh; - neighbor->binsizex = xprd / neighbor->nbinx; - neighbor->binsizey = yprd / neighbor->nbiny; - neighbor->binsizez = zprd / neighbor->nbinz; - - neighbor->bininvx = 1.0 / neighbor->binsizex; - neighbor->bininvy = 1.0 / neighbor->binsizey; - neighbor->bininvz = 1.0 / neighbor->binsizez; - - coord = xlo - Cutneigh - SMALL * xprd; - neighbor->mbinxlo = (int) (coord * neighbor->bininvx); - if (coord < 0.0) { - neighbor->mbinxlo = neighbor->mbinxlo - 1; - } - coord = xhi + Cutneigh + SMALL * xprd; - mbinxhi = (int) (coord * neighbor->bininvx); - - coord = ylo - Cutneigh - SMALL * yprd; - neighbor->mbinylo = (int) (coord * neighbor->bininvy); - if (coord < 0.0) { - neighbor->mbinylo = neighbor->mbinylo - 1; - } - coord = yhi + Cutneigh + SMALL * yprd; - mbinyhi = (int) (coord * neighbor->bininvy); - - coord = zlo - Cutneigh - SMALL * zprd; - neighbor->mbinzlo = (int) (coord * neighbor->bininvz); - if (coord < 0.0) { - neighbor->mbinzlo = neighbor->mbinzlo - 1; - } - coord = zhi + Cutneigh + SMALL * zprd; - mbinzhi = (int) (coord * neighbor->bininvz); - - neighbor->mbinxlo = neighbor->mbinxlo - 1; - mbinxhi = mbinxhi + 1; - neighbor->mbinx = mbinxhi - neighbor->mbinxlo + 1; - - neighbor->mbinylo = neighbor->mbinylo - 1; - mbinyhi = mbinyhi + 1; - neighbor->mbiny = mbinyhi - neighbor->mbinylo + 1; - - neighbor->mbinzlo = neighbor->mbinzlo - 1; - mbinzhi = mbinzhi + 1; - neighbor->mbinz = mbinzhi - neighbor->mbinzlo + 1; - - nextx = (int) (Cutneigh * neighbor->bininvx); - if(nextx * neighbor->binsizex < FACTOR * Cutneigh) nextx++; - - nexty = (int) (Cutneigh * neighbor->bininvy); - if(nexty * neighbor->binsizey < FACTOR * Cutneigh) nexty++; - - nextz = (int) (Cutneigh * neighbor->bininvz); - if(nextz * neighbor->binsizez < FACTOR * Cutneigh) nextz++; - - if (neighbor->stencil) { - free(neighbor->stencil); - } - - neighbor->stencil = (int*) malloc( - (2 * nextz + 1) * (2 * nexty + 1) * (2 * nextx + 1) * sizeof(int)); - - neighbor->nstencil = 0; - int kstart = -nextz; - - for(int k = kstart; k <= nextz; k++) { - for(int j = -nexty; j <= nexty; j++) { - for(int i = -nextx; i <= nextx; i++) { - if(bindist(neighbor, i, j, k) < neighbor->cutneighsq) { - neighbor->stencil[neighbor->nstencil++] = - k * neighbor->mbiny * neighbor->mbinx + j * neighbor->mbinx + i; - } - } - } - } - - neighbor->mbins = neighbor->mbinx * neighbor->mbiny * neighbor->mbinz; - - if (neighbor->bincount) { - free(neighbor->bincount); - } - neighbor->bincount = (int*) malloc(neighbor->mbins * sizeof(int)); - - if (neighbor->bins) { - free(neighbor->bins); - } - neighbor->bins = (int*) malloc(neighbor->mbins * neighbor->atoms_per_bin * sizeof(int)); + param->xprd = param->nx * lattice; + param->yprd = param->ny * lattice; + param->zprd = param->nz * lattice; } -double* myrealloc(double *ptr, int n, int nold) { - - double* newarray; - newarray = (double*) malloc(n * sizeof(double)); - - if(nold) { - memcpy(newarray, ptr, nold * sizeof(double)); - } - if(ptr) { - free(ptr); - } - - return newarray; -} - -int* myreallocInt(int *ptr, int n, int nold) { - - int* newarray; - - newarray = (int*) malloc(n * sizeof(int)); - - if(nold) { - memcpy(newarray, ptr, nold * sizeof(int)); - } - if(ptr) { - free(ptr); - } - - return newarray; -} - -void growBoundary() +void initialIntegrate(Parameter *param, Atom *atom) { - int nold = NmaxGhost; - NmaxGhost += DELTA; + double* x = atom->x; double* y = atom->y; double* z = atom->z; + double* fx = atom->fx; double* fy = atom->fy; double* fz = atom->fz; + double* vx = atom->vx; double* vy = atom->vy; double* vz = atom->vz; - BorderMap = myreallocInt(BorderMap, NmaxGhost, nold); - PBCx = myreallocInt(PBCx, NmaxGhost, nold); - PBCy = myreallocInt(PBCy, NmaxGhost, nold); - PBCz = myreallocInt(PBCz, NmaxGhost, nold); - - if(BorderMap == NULL || PBCx == NULL || PBCy == NULL || PBCz == NULL ) { - printf("ERROR: No memory for Boundary\n"); - } -} - -void growarray() -{ - int nold = Nmax; - Nmax += DELTA; - - x = myrealloc(x, Nmax, nold); y = myrealloc(y, Nmax, nold); z = myrealloc(z, Nmax, nold); - vx = myrealloc(vx, Nmax, nold); vy = myrealloc(vy, Nmax, nold); vz = myrealloc(vz, Nmax, nold); - fx = myrealloc(fx, Nmax, nold); fy = myrealloc(fy, Nmax, nold); fz = myrealloc(fz, Nmax, nold); - - if(x == NULL || y == NULL || z == NULL || - vx == NULL || vy == NULL || vz == NULL || - fx == NULL || fy == NULL || fz == NULL ) { - printf("ERROR: No memory for atoms\n"); - } -} - -void updateBorders() -{ - for(int i = 0; i < Nghost; i++) { - x[Nlocal + i] = x[BorderMap[i]] + PBCx[i] * xprd; - y[Nlocal + i] = y[BorderMap[i]] + PBCy[i] * yprd; - z[Nlocal + i] = z[BorderMap[i]] + PBCz[i] * zprd; - } -} - -void updateAtomLocations() -{ - for(int i = 0; i < Nlocal; i++) { - - if(x[i] < 0.0) { - x[i] += xprd; - } else if(x[i] >= xprd) { - x[i] -= xprd; - } - - if(y[i] < 0.0) { - y[i] += yprd; - } else if(y[i] >= yprd) { - y[i] -= yprd; - } - - if(z[i] < 0.0) { - z[i] += zprd; - } else if(z[i] >= zprd) { - z[i] -= zprd; - } - } -} - -void setupBordersNew() -{ - int lastidx = 0; - int nghostprev = 0; - Nghost = 0; - - for (int i = 0; i < Nlocal; i++) { - - if (Nlocal + Nghost + 1 >= Nmax) { - growarray(); - } - - if (x[i] < Cutneigh) { - Nghost++; - x[i+lastidx] = x[i] + xprd; - y[i+lastidx] = y[i]; - z[i+lastidx] = z[i]; - lastidx++; - } else if (x[i] >= xprd - Cutneigh) { - Nghost++; - x[i+lastidx] = x[i] - xprd; - y[i+lastidx] = y[i]; - z[i+lastidx] = z[i]; - lastidx++; - } - } - - nghostprev = Nghost+1; - - for (int i = 0; i < Nlocal + nghostprev ; i++) { - - if (Nlocal + Nghost + 1 >= Nmax) { - growarray(); - } - - if (y[i] < Cutneigh) { - Nghost++; - x[i+lastidx] = x[i]; - y[i+lastidx] = y[i] + yprd; - z[i+lastidx] = z[i]; - lastidx++; - } else if (y[i] >= yprd - Cutneigh) { - Nghost++; - x[i+lastidx] = x[i]; - y[i+lastidx] = y[i] - yprd; - z[i+lastidx] = z[i]; - lastidx++; - } - } - - nghostprev = Nghost+1; - - for (int i = 0; i < Nlocal + nghostprev; i++) { - - if (Nlocal + Nghost + 1 >= Nmax) { - growarray(); - } - - if (z[i] < Cutneigh) { - Nghost++; - x[i+lastidx] = x[i]; - y[i+lastidx] = y[i]; - z[i+lastidx] = z[i] + zprd; - lastidx++; - } else if(z[i] >= zprd - Cutneigh) { - Nghost++; - x[i+lastidx] = x[i]; - y[i+lastidx] = y[i]; - z[i+lastidx] = z[i] - zprd; - lastidx++; - } - } - - Nghost++; -} - -#define ADDGHOST(dx,dy,dz) Nghost++; BorderMap[Nghost] = i; PBCx[Nghost] = dx; PBCy[Nghost] = dy; PBCz[Nghost] = dz; -void setupBorders() -{ - Nghost = -1; - - for(int i = 0; i < Nlocal; i++) { - - if (Nlocal + Nghost + 7 >= Nmax) { - growarray(); - } - if (Nghost + 7 >= NmaxGhost) { - growBoundary(); - } - - /* Setup ghost atoms */ - /* 6 planes */ - if (x[i] < Cutneigh) { ADDGHOST(+1,0,0); } - if (x[i] >= (xprd-Cutneigh)) { ADDGHOST(-1,0,0); } - if (y[i] < Cutneigh) { ADDGHOST(0,+1,0); } - if (y[i] >= (yprd-Cutneigh)) { ADDGHOST(0,-1,0); } - if (z[i] < Cutneigh) { ADDGHOST(0,0,+1); } - if (z[i] >= (zprd-Cutneigh)) { ADDGHOST(0,0,-1); } - /* 8 corners */ - if (x[i] < Cutneigh && y[i] < Cutneigh && z[i] < Cutneigh) { ADDGHOST(+1,+1,+1); } - if (x[i] < Cutneigh && y[i] >= (yprd-Cutneigh) && z[i] < Cutneigh) { ADDGHOST(+1,-1,+1); } - if (x[i] < Cutneigh && y[i] >= Cutneigh && z[i] >= (zprd-Cutneigh)) { ADDGHOST(+1,+1,-1); } - if (x[i] < Cutneigh && y[i] >= (yprd-Cutneigh) && z[i] >= (zprd-Cutneigh)) { ADDGHOST(+1,-1,-1); } - if (x[i] >= (xprd-Cutneigh) && y[i] < Cutneigh && z[i] < Cutneigh) { ADDGHOST(-1,+1,+1); } - if (x[i] >= (xprd-Cutneigh) && y[i] >= (yprd-Cutneigh) && z[i] < Cutneigh) { ADDGHOST(-1,-1,+1); } - if (x[i] >= (xprd-Cutneigh) && y[i] < Cutneigh && z[i] >= (zprd-Cutneigh)) { ADDGHOST(-1,+1,-1); } - if (x[i] >= (xprd-Cutneigh) && y[i] >= (yprd-Cutneigh) && z[i] >= (zprd-Cutneigh)) { ADDGHOST(-1,-1,-1); } - /* 12 edges */ - if (x[i] < Cutneigh && z[i] < Cutneigh) { ADDGHOST(+1,0,+1); } - if (x[i] < Cutneigh && z[i] >= (zprd-Cutneigh)) { ADDGHOST(+1,0,-1); } - if (x[i] >= (xprd-Cutneigh) && z[i] < Cutneigh) { ADDGHOST(-1,0,+1); } - if (x[i] >= (xprd-Cutneigh) && z[i] >= (zprd-Cutneigh)) { ADDGHOST(-1,0,-1); } - if (y[i] < Cutneigh && z[i] < Cutneigh) { ADDGHOST(0,+1,+1); } - if (y[i] < Cutneigh && z[i] >= (zprd-Cutneigh)) { ADDGHOST(0,+1,-1); } - if (y[i] >= (yprd-Cutneigh) && z[i] < Cutneigh) { ADDGHOST(0,-1,+1); } - if (y[i] >= (yprd-Cutneigh) && z[i] >= (zprd-Cutneigh)) { ADDGHOST(0,-1,-1); } - if (y[i] < Cutneigh && x[i] < Cutneigh) { ADDGHOST(+1,+1,0); } - if (y[i] < Cutneigh && x[i] >= (xprd-Cutneigh)) { ADDGHOST(-1,+1,0); } - if (y[i] >= (yprd-Cutneigh) && x[i] < Cutneigh) { ADDGHOST(+1,-1,0); } - if (y[i] >= (yprd-Cutneigh) && x[i] >= (xprd-Cutneigh)) { ADDGHOST(-1,-1,0); } - } - // increase by one to make it the ghost atom count - Nghost++; -} - -void sortAtoms(Neighbor *neighbor) -{ - binatoms(neighbor); - int* binpos = neighbor->bincount; - int* bins = neighbor->bins; - - int mbins = neighbor->mbins; - int atoms_per_bin = neighbor->atoms_per_bin; - - for(int i=1; i0?binpos[mybin-1]:0; - int count = binpos[mybin] - start; - - for(int k=0; knx * param->ny * param->nz; - Nlocal = 0; - double alat = pow((4.0 / param->rho), (1.0 / 3.0)); - int ilo = (int) (xlo / (0.5 * alat) - 1); - int ihi = (int) (xhi / (0.5 * alat) + 1); - int jlo = (int) (ylo / (0.5 * alat) - 1); - int jhi = (int) (yhi / (0.5 * alat) + 1); - int klo = (int) (zlo / (0.5 * alat) - 1); - int khi = (int) (zhi / (0.5 * alat) + 1); - - ilo = MAX(ilo, 0); - ihi = MIN(ihi, 2 * param->nx - 1); - jlo = MAX(jlo, 0); - jhi = MIN(jhi, 2 * param->ny - 1); - klo = MAX(klo, 0); - khi = MIN(khi, 2 * param->nz - 1); - - double xtmp, ytmp, ztmp, vxtmp, vytmp, vztmp; - int i, j, k, m, n; - int sx = 0; int sy = 0; int sz = 0; - int ox = 0; int oy = 0; int oz = 0; - int subboxdim = 8; - - while(oz * subboxdim <= khi) { - - k = oz * subboxdim + sz; - j = oy * subboxdim + sy; - i = ox * subboxdim + sx; - - if(((i + j + k) % 2 == 0) && - (i >= ilo) && (i <= ihi) && - (j >= jlo) && (j <= jhi) && - (k >= klo) && (k <= khi)) { - - xtmp = 0.5 * alat * i; - ytmp = 0.5 * alat * j; - ztmp = 0.5 * alat * k; - - if( xtmp >= xlo && xtmp < xhi && - ytmp >= ylo && ytmp < yhi && - ztmp >= zlo && ztmp < zhi ) { - - n = k * (2 * param->ny) * (2 * param->nx) + - j * (2 * param->nx) + - i + 1; - - for(m = 0; m < 5; m++) { - myrandom(&n); - } - vxtmp = myrandom(&n); - - for(m = 0; m < 5; m++){ - myrandom(&n); - } - vytmp = myrandom(&n); - - for(m = 0; m < 5; m++) { - myrandom(&n); - } - vztmp = myrandom(&n); - - if(Nlocal == Nmax) { - growarray(); - } - - x[Nlocal] = xtmp; y[Nlocal] = ytmp; z[Nlocal] = ztmp; - vx[Nlocal] = vxtmp; vy[Nlocal] = vytmp; vz[Nlocal] = vztmp; - Nlocal++; - } - } - - sx++; - - if(sx == subboxdim) { sx = 0; sy++; } - if(sy == subboxdim) { sy = 0; sz++; } - if(sz == subboxdim) { sz = 0; ox++; } - if(ox * subboxdim > ihi) { ox = 0; oy++; } - if(oy * subboxdim > jhi) { oy = 0; oz++; } - } -} - -void adjustVelocity(Parameter *param, Thermo *thermo) -{ - /* zero center-of-mass motion */ - double vxtot = 0.0; - double vytot = 0.0; - double vztot = 0.0; - - for(int i = 0; i < Nlocal; i++) { - vxtot += vx[i]; - vytot += vy[i]; - vztot += vz[i]; - } - - vxtot = vxtot / Natoms; - vytot = vytot / Natoms; - vztot = vztot / Natoms; - - for(int i = 0; i < Nlocal; i++) { - vx[i] -= vxtot; - vy[i] -= vytot; - vz[i] -= vztot; - } - - thermo->t_act = 0; - double t = 0.0; - - for(int i = 0; i < Nlocal; i++) { - t += (vx[i] * vx[i] + vy[i] * vy[i] + vz[i] * vz[i]) * param->mass; - } - - t *= thermo->t_scale; - double factor = sqrt(param->temp / t); - - for(int i = 0; i < Nlocal; i++) { - vx[i] *= factor; - vy[i] *= factor; - vz[i] *= factor; - } -} - -void thermoSetup(Parameter *param, Thermo *thermo) -{ - int maxstat = param->ntimes / param->nstat + 2; - - thermo->steparr = (int*) malloc(maxstat * sizeof(int)); - thermo->tmparr = (double*) malloc(maxstat * sizeof(double)); - thermo->engarr = (double*) malloc(maxstat * sizeof(double)); - thermo->prsarr = (double*) malloc(maxstat * sizeof(double)); - - thermo->mvv2e = 1.0; - thermo->dof_boltz = (Natoms * 3 - 3); - thermo->t_scale = thermo->mvv2e / thermo->dof_boltz; - thermo->p_scale = 1.0 / 3 / xprd / yprd / zprd; - thermo->e_scale = 0.5; - - printf("step\ttemp\t\tpressure\n"); -} - - -void thermoCompute(int iflag, Parameter *param, Thermo *thermo) -{ - double t = 0.0, p; - - for(int i = 0; i < Nlocal; i++) { - t += (vx[i] * vx[i] + vy[i] * vy[i] + vz[i] * vz[i]) * param->mass; - } - - t = t * thermo->t_scale; - p = (t * thermo->dof_boltz) * thermo->p_scale; - - int istep = iflag; - - if(iflag == -1){ - istep = param->ntimes; - } - if(iflag == 0){ - thermo->mstat = 0; - } - - thermo->steparr[thermo->mstat] = istep; - thermo->tmparr[thermo->mstat] = t; - thermo->prsarr[thermo->mstat] = p; - thermo->mstat++; - fprintf(stdout, "%i\t%e\t%e\n", istep, t, p); -} - -void initialIntegrate(Parameter *param) -{ - for(int i = 0; i < Nlocal; i++) { + for(int i = 0; i < atom->Nlocal; i++) { vx[i] += param->dtforce * fx[i]; vy[i] += param->dtforce * fy[i]; vz[i] += param->dtforce * fz[i]; @@ -899,21 +89,27 @@ void initialIntegrate(Parameter *param) } } -void finalIntegrate(Parameter *param) +void finalIntegrate(Parameter *param, Atom *atom) { - for(int i = 0; i < Nlocal; i++) { + double* fx = atom->fx; double* fy = atom->fy; double* fz = atom->fz; + double* vx = atom->vx; double* vy = atom->vy; double* vz = atom->vz; + + for(int i = 0; i < atom->Nlocal; i++) { vx[i] += param->dtforce * fx[i]; vy[i] += param->dtforce * fy[i]; vz[i] += param->dtforce * fz[i]; } } -void computeForce(Neighbor *neighbor, Parameter *param) +void computeForce(Neighbor *neighbor, Atom *atom, Parameter *param) { + int Nlocal = atom->Natoms; int* neighs; double cutforcesq = param->cutforce * param->cutforce; double sigma6 = param->sigma6; double epsilon = param->epsilon; + double* x = atom->x; double* y = atom->y; double* z = atom->z; + double* fx = atom->fx; double* fy = atom->fy; double* fz = atom->fz; for(int i = 0; i < Nlocal; i++) { fx[i] = 0.0; @@ -955,46 +151,63 @@ void computeForce(Neighbor *neighbor, Parameter *param) } } + +void printAtomState(Atom *atom) +{ + printf("Atom counts: Natoms=%d Nlocal=%d Nghost=%d Nmax=%d\n", + atom->Natoms, atom->Nlocal, atom->Nghost, atom->Nmax); + + int nall = atom->Nlocal + atom->Nghost; + + for (int i=0; ix[i], atom->y[i], atom->z[i]); + } +} + + int main (int argc, char** argv) { + Atom atom; Neighbor neighbor; Parameter param; - Thermo thermo; - init(&neighbor, ¶m); - setup(&neighbor, ¶m); - create_atoms(¶m); - thermoSetup(¶m, &thermo); - adjustVelocity(¶m, &thermo); - setupBorders(); - updateBorders(); - buildNeighborlist(&neighbor); - thermoCompute(0, ¶m, &thermo); - computeForce(&neighbor, ¶m); + init(&atom, ¶m); + initNeighbor(&neighbor, ¶m); + initPbc(); + setupNeighbor(¶m); + createAtom(&atom, ¶m); + setupThermo(¶m, atom.Natoms); + adjustThermo(¶m, &atom); + /* printAtomState(&atom); */ + setupPbc(&atom, ¶m); + updatePbc(&atom, ¶m); + buildNeighbor(&atom, &neighbor); + computeThermo(0, ¶m, &atom); + computeForce(&neighbor, &atom, ¶m); for(int n = 0; n < param.ntimes; n++) { - initialIntegrate(¶m); + initialIntegrate(¶m, &atom); - if((n + 1) % neighbor.every) { - updateBorders(); + if((n + 1) % param.every) { + updatePbc(&atom, ¶m); } else { - updateAtomLocations(); - setupBorders(); - updateBorders(); - /* sortAtoms(&neighbor); */ - buildNeighborlist(&neighbor); + updateAtomsPbc(&atom, ¶m); + setupPbc(&atom, ¶m); + updatePbc(&atom, ¶m); + /* sortAtom(); */ + buildNeighbor(&atom, &neighbor); } - computeForce(&neighbor, ¶m); - finalIntegrate(¶m); + computeForce(&neighbor, &atom, ¶m); + finalIntegrate(¶m, &atom); if(!((n + 1) % param.nstat) && (n+1) < param.ntimes) { - thermoCompute(n + 1, ¶m, &thermo); + computeThermo(n + 1, ¶m, &atom); } } - thermoCompute(-1, ¶m, &thermo); + computeThermo(-1, ¶m, &atom); return EXIT_SUCCESS; } diff --git a/src/neighbor.c b/src/neighbor.c new file mode 100644 index 0000000..a39cb48 --- /dev/null +++ b/src/neighbor.c @@ -0,0 +1,345 @@ +/* + * ======================================================================================= + * + * Author: Jan Eitzinger (je), jan.eitzinger@fau.de + * Copyright (c) 2020 RRZE, University Erlangen-Nuremberg + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in all + * copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + * + * ======================================================================================= + */ +#include +#include +#include + +#include +#include +#include + +#define SMALL 1.0e-6 +#define FACTOR 0.999 + +static double xprd, yprd, zprd; +static double bininvx, bininvy, bininvz; +static int mbinxlo, mbinylo, mbinzlo; +static int nbinx, nbiny, nbinz; +static int mbinx, mbiny, mbinz; // n bins in x, y, z +static int *bincount; +static int *bins; +static int mbins; //total number of bins +static int atoms_per_bin; // max atoms per bin +static double cutneigh; +static double cutneighsq; // neighbor cutoff squared +static int nmax; +static int nstencil; // # of bins in stencil +static int* stencil; // stencil list of bin offsets +static double binsizex, binsizey, binsizez; + +static int coord2bin(double, double , double); +static double bindist(int, int, int); + +/* exported subroutines */ +void initNeighbor(Neighbor *neighbor, Parameter *param) +{ + double lattice = pow((4.0 / param->rho), (1.0 / 3.0)); + double neighscale = 5.0 / 6.0; + xprd = param->nx * lattice; + yprd = param->ny * lattice; + zprd = param->nz * lattice; + cutneigh = param->cutneigh; + nbinx = neighscale * param->nx; + nbiny = neighscale * param->ny; + nbinz = neighscale * param->nz; + nmax = 0; + atoms_per_bin = 8; + stencil = NULL; + bins = NULL; + bincount = NULL; + neighbor->maxneighs = 100; + neighbor->numneigh = NULL; + neighbor->neighbors = NULL; +} + +void setupNeighbor(Parameter *param) +{ + double coord; + int mbinxhi, mbinyhi, mbinzhi; + int nextx, nexty, nextz; + double xlo = 0.0; double xhi = xprd; + double ylo = 0.0; double yhi = yprd; + double zlo = 0.0; double 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; + + coord = xlo - cutneigh - SMALL * xprd; + mbinxlo = (int) (coord * bininvx); + if (coord < 0.0) { + mbinxlo = mbinxlo - 1; + } + coord = xhi + cutneigh + SMALL * xprd; + mbinxhi = (int) (coord * bininvx); + + coord = ylo - cutneigh - SMALL * yprd; + mbinylo = (int) (coord * bininvy); + if (coord < 0.0) { + mbinylo = mbinylo - 1; + } + coord = yhi + cutneigh + SMALL * yprd; + mbinyhi = (int) (coord * bininvy); + + coord = zlo - cutneigh - SMALL * zprd; + mbinzlo = (int) (coord * bininvz); + if (coord < 0.0) { + mbinzlo = mbinzlo - 1; + } + coord = zhi + cutneigh + SMALL * zprd; + mbinzhi = (int) (coord * bininvz); + + mbinxlo = mbinxlo - 1; + mbinxhi = mbinxhi + 1; + mbinx = mbinxhi - mbinxlo + 1; + + mbinylo = mbinylo - 1; + mbinyhi = mbinyhi + 1; + mbiny = mbinyhi - mbinylo + 1; + + mbinzlo = mbinzlo - 1; + mbinzhi = mbinzhi + 1; + mbinz = mbinzhi - mbinzlo + 1; + + nextx = (int) (cutneigh * bininvx); + if(nextx * binsizex < FACTOR * cutneigh) nextx++; + + nexty = (int) (cutneigh * bininvy); + if(nexty * binsizey < FACTOR * cutneigh) nexty++; + + nextz = (int) (cutneigh * bininvz); + if(nextz * binsizez < FACTOR * cutneigh) nextz++; + + if (stencil) { + free(stencil); + } + + stencil = (int*) malloc( + (2 * nextz + 1) * (2 * nexty + 1) * (2 * nextx + 1) * sizeof(int)); + + nstencil = 0; + int kstart = -nextz; + + for(int k = kstart; k <= nextz; k++) { + for(int j = -nexty; j <= nexty; j++) { + for(int i = -nextx; i <= nextx; i++) { + if(bindist(i, j, k) < cutneighsq) { + stencil[nstencil++] = + k * mbiny * mbinx + j * mbinx + i; + } + } + } + } + + mbins = mbinx * mbiny * mbinz; + + if (bincount) { + free(bincount); + } + bincount = (int*) malloc(mbins * sizeof(int)); + + if (bins) { + free(bins); + } + bins = (int*) malloc(mbins * atoms_per_bin * sizeof(int)); +} + +void buildNeighbor(Atom *atom, Neighbor *neighbor) +{ + int nall = atom->Nlocal + atom->Nghost; + + /* extend atom arrays if necessary */ + if(nall > nmax) { + nmax = nall; + if(neighbor->numneigh) free(neighbor->numneigh); + if(neighbor->neighbors) free(neighbor->neighbors); + neighbor->numneigh = (int*) malloc(nmax * sizeof(int)); + neighbor->neighbors = (int*) malloc(nmax * neighbor->maxneighs * sizeof(int*)); + } + + /* bin local & ghost atoms */ + binatoms(atom); + int resize = 1; + double* x = atom->x; + double* y = atom->y; + double* z = atom->z; + + /* loop over each atom, storing neighbors */ + while(resize) { + int new_maxneighs = neighbor->maxneighs; + resize = 0; + + for(int i = 0; i < atom->Nlocal; i++) { + int* neighptr = &(neighbor->neighbors[i * neighbor->maxneighs]); + int n = 0; + double xtmp = x[i]; + double ytmp = y[i]; + double ztmp = z[i]; + int ibin = coord2bin(xtmp, ytmp, ztmp); + + for(int k = 0; k < nstencil; k++) { + int jbin = ibin + stencil[k]; + int* loc_bin = &bins[jbin * atoms_per_bin]; + + for(int m = 0; m < bincount[jbin]; m++) { + int j = loc_bin[m]; + + if ( j == i ){ + continue; + } + + double delx = xtmp - x[j]; + double dely = ytmp - y[j]; + double delz = ztmp - z[j]; + double rsq = delx * delx + dely * dely + delz * delz; + + if( rsq <= cutneighsq ) { + neighptr[n++] = j; + } + } + } + + neighbor->numneigh[i] = n; + + if(n >= neighbor->maxneighs) { + resize = 1; + + if(n >= new_maxneighs) { + new_maxneighs = n; + } + } + } + + if(resize) { + printf("RESIZE %d\n", neighbor->maxneighs); + neighbor->maxneighs = new_maxneighs * 1.2; + free(neighbor->neighbors); + neighbor->neighbors = (int*) malloc(atom->Nmax * neighbor->maxneighs * sizeof(int)); + } + } +} + +/* internal subroutines */ +double bindist(int i, int j, int k) +{ + double delx, dely, delz; + + if(i > 0) { + delx = (i - 1) * binsizex; + } else if(i == 0) { + delx = 0.0; + } else { + delx = (i + 1) * binsizex; + } + + if(j > 0) { + dely = (j - 1) * binsizey; + } else if(j == 0) { + dely = 0.0; + } else { + dely = (j + 1) * binsizey; + } + + if(k > 0) { + delz = (k - 1) * binsizez; + } else if(k == 0) { + delz = 0.0; + } else { + delz = (k + 1) * binsizez; + } + + return (delx * delx + dely * dely + delz * delz); +} + +int coord2bin(double xin, double yin, double zin) +{ + int ix, iy, iz; + + if(xin >= xprd) { + ix = (int)((xin - xprd) * bininvx) + nbinx - mbinxlo; + } else if(xin >= 0.0) { + ix = (int)(xin * bininvx) - mbinxlo; + } else { + ix = (int)(xin * bininvx) - mbinxlo - 1; + } + + if(yin >= yprd) { + iy = (int)((yin - yprd) * bininvy) + nbiny - mbinylo; + } else if(yin >= 0.0) { + iy = (int)(yin * bininvy) - mbinylo; + } else { + iy = (int)(yin * bininvy) - mbinylo - 1; + } + + if(zin >= zprd) { + iz = (int)((zin - zprd) * bininvz) + nbinz - mbinzlo; + } else if(zin >= 0.0) { + iz = (int)(zin * bininvz) - mbinzlo; + } else { + iz = (int)(zin * bininvz) - mbinzlo - 1; + } + + return (iz * mbiny * mbinx + iy * mbinx + ix + 1); +} + +void binatoms(Atom *atom) +{ + int nall = atom->Nlocal + atom->Nghost; + double* x = atom->x; + double* y = atom->y; + double* z = atom->z; + int resize = 1; + + while(resize > 0) { + resize = 0; + + for(int i = 0; i < mbins; i++) { + bincount[i] = 0; + } + + for(int i = 0; i < nall; i++) { + int ibin = coord2bin(x[i], y[i], z[i]); + + if(bincount[ibin] < atoms_per_bin) { + int ac = bincount[ibin]++; + bins[ibin * atoms_per_bin + ac] = i; + } else { + resize = 1; + } + } + + if(resize) { + free(bins); + atoms_per_bin *= 2; + bins = (int*) malloc(mbins * atoms_per_bin * sizeof(int)); + } + } +} diff --git a/src/pbc.c b/src/pbc.c new file mode 100644 index 0000000..ba4a0ae --- /dev/null +++ b/src/pbc.c @@ -0,0 +1,161 @@ +/* + * ======================================================================================= + * + * Author: Jan Eitzinger (je), jan.eitzinger@fau.de + * Copyright (c) 2020 RRZE, University Erlangen-Nuremberg + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in all + * copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + * + * ======================================================================================= + */ +#include +#include + +#include +#include +#include + +#define DELTA 20000 + +static int NmaxGhost; +static int *BorderMap; +static int *PBCx, *PBCy, *PBCz; + +static void growPbc(); + +void growPbc() +{ + int nold = NmaxGhost; + NmaxGhost += DELTA; + + BorderMap = (int*) reallocate(BorderMap, ALIGNMENT, NmaxGhost * sizeof(int), nold * sizeof(int)); + PBCx = (int*) reallocate(PBCx, ALIGNMENT, NmaxGhost * sizeof(int), nold * sizeof(int)); + PBCy = (int*) reallocate(PBCy, ALIGNMENT, NmaxGhost * sizeof(int), nold * sizeof(int)); + PBCz = (int*) reallocate(PBCz, ALIGNMENT, NmaxGhost * sizeof(int), nold * sizeof(int)); +} + +void initPbc() +{ + NmaxGhost = 0; + BorderMap = NULL; + PBCx = NULL; PBCy = NULL; PBCz = NULL; +} + +void updatePbc(Atom *atom, Parameter *param) +{ + int nlocal = atom->Nlocal; + double* x = atom->x; + double* y = atom->y; + double* z = atom->z; + double xprd = param->xprd; + double yprd = param->yprd; + double zprd = param->zprd; + + for(int i = 0; i < atom->Nghost; i++) { + x[nlocal + i] = x[BorderMap[i]] + PBCx[i] * xprd; + y[nlocal + i] = y[BorderMap[i]] + PBCy[i] * yprd; + z[nlocal + i] = z[BorderMap[i]] + PBCz[i] * zprd; + } +} + +void updateAtomsPbc(Atom *atom, Parameter *param) +{ + double* x = atom->x; + double* y = atom->y; + double* z = atom->z; + double xprd = param->xprd; + double yprd = param->yprd; + double zprd = param->zprd; + + for(int i = 0; i < atom->Nlocal; i++) { + + if(x[i] < 0.0) { + x[i] += xprd; + } else if(x[i] >= xprd) { + x[i] -= xprd; + } + + if(y[i] < 0.0) { + y[i] += yprd; + } else if(y[i] >= yprd) { + y[i] -= yprd; + } + + if(z[i] < 0.0) { + z[i] += zprd; + } else if(z[i] >= zprd) { + z[i] -= zprd; + } + } +} + +#define ADDGHOST(dx,dy,dz) Nghost++; BorderMap[Nghost] = i; PBCx[Nghost] = dx; PBCy[Nghost] = dy; PBCz[Nghost] = dz; +void setupPbc(Atom *atom, Parameter *param) +{ + double* x = atom->x; double* y = atom->y; double* z = atom->z; + double xprd = param->xprd; + double yprd = param->yprd; + double zprd = param->zprd; + double Cutneigh = param->cutneigh; + int Nghost = -1; + + for(int i = 0; i < atom->Nlocal; i++) { + + if (atom->Nlocal + Nghost + 7 >= atom->Nmax) { + growAtom(atom); + x = atom->x; y = atom->y; z = atom->z; + } + if (Nghost + 7 >= NmaxGhost) { + growPbc(); + } + + /* Setup ghost atoms */ + /* 6 planes */ + if (x[i] < Cutneigh) { ADDGHOST(+1,0,0); } + if (x[i] >= (xprd-Cutneigh)) { ADDGHOST(-1,0,0); } + if (y[i] < Cutneigh) { ADDGHOST(0,+1,0); } + if (y[i] >= (yprd-Cutneigh)) { ADDGHOST(0,-1,0); } + if (z[i] < Cutneigh) { ADDGHOST(0,0,+1); } + if (z[i] >= (zprd-Cutneigh)) { ADDGHOST(0,0,-1); } + /* 8 corners */ + if (x[i] < Cutneigh && y[i] < Cutneigh && z[i] < Cutneigh) { ADDGHOST(+1,+1,+1); } + if (x[i] < Cutneigh && y[i] >= (yprd-Cutneigh) && z[i] < Cutneigh) { ADDGHOST(+1,-1,+1); } + if (x[i] < Cutneigh && y[i] >= Cutneigh && z[i] >= (zprd-Cutneigh)) { ADDGHOST(+1,+1,-1); } + if (x[i] < Cutneigh && y[i] >= (yprd-Cutneigh) && z[i] >= (zprd-Cutneigh)) { ADDGHOST(+1,-1,-1); } + if (x[i] >= (xprd-Cutneigh) && y[i] < Cutneigh && z[i] < Cutneigh) { ADDGHOST(-1,+1,+1); } + if (x[i] >= (xprd-Cutneigh) && y[i] >= (yprd-Cutneigh) && z[i] < Cutneigh) { ADDGHOST(-1,-1,+1); } + if (x[i] >= (xprd-Cutneigh) && y[i] < Cutneigh && z[i] >= (zprd-Cutneigh)) { ADDGHOST(-1,+1,-1); } + if (x[i] >= (xprd-Cutneigh) && y[i] >= (yprd-Cutneigh) && z[i] >= (zprd-Cutneigh)) { ADDGHOST(-1,-1,-1); } + /* 12 edges */ + if (x[i] < Cutneigh && z[i] < Cutneigh) { ADDGHOST(+1,0,+1); } + if (x[i] < Cutneigh && z[i] >= (zprd-Cutneigh)) { ADDGHOST(+1,0,-1); } + if (x[i] >= (xprd-Cutneigh) && z[i] < Cutneigh) { ADDGHOST(-1,0,+1); } + if (x[i] >= (xprd-Cutneigh) && z[i] >= (zprd-Cutneigh)) { ADDGHOST(-1,0,-1); } + if (y[i] < Cutneigh && z[i] < Cutneigh) { ADDGHOST(0,+1,+1); } + if (y[i] < Cutneigh && z[i] >= (zprd-Cutneigh)) { ADDGHOST(0,+1,-1); } + if (y[i] >= (yprd-Cutneigh) && z[i] < Cutneigh) { ADDGHOST(0,-1,+1); } + if (y[i] >= (yprd-Cutneigh) && z[i] >= (zprd-Cutneigh)) { ADDGHOST(0,-1,-1); } + if (y[i] < Cutneigh && x[i] < Cutneigh) { ADDGHOST(+1,+1,0); } + if (y[i] < Cutneigh && x[i] >= (xprd-Cutneigh)) { ADDGHOST(-1,+1,0); } + if (y[i] >= (yprd-Cutneigh) && x[i] < Cutneigh) { ADDGHOST(+1,-1,0); } + if (y[i] >= (yprd-Cutneigh) && x[i] >= (xprd-Cutneigh)) { ADDGHOST(-1,-1,0); } + } + // increase by one to make it the ghost atom count + atom->Nghost = Nghost + 1; +} diff --git a/src/thermo.c b/src/thermo.c new file mode 100644 index 0000000..b3bd651 --- /dev/null +++ b/src/thermo.c @@ -0,0 +1,132 @@ +/* + * ======================================================================================= + * + * Author: Jan Eitzinger (je), jan.eitzinger@fau.de + * Copyright (c) 2020 RRZE, University Erlangen-Nuremberg + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in all + * copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + * + * ======================================================================================= + */ +#include +#include +#include + +#include + +static int *steparr; +static double *tmparr; +static double *engarr; +static double *prsarr; +static double mvv2e; +static int dof_boltz; +static double t_scale; +static double p_scale; +static double e_scale; +static double t_act; +static double p_act; +static double e_act; +static int mstat; + +void setupThermo(Parameter *param, int natoms) +{ + int maxstat = param->ntimes / param->nstat + 2; + + steparr = (int*) malloc(maxstat * sizeof(int)); + tmparr = (double*) malloc(maxstat * sizeof(double)); + engarr = (double*) malloc(maxstat * sizeof(double)); + prsarr = (double*) malloc(maxstat * sizeof(double)); + + mvv2e = 1.0; + dof_boltz = (natoms * 3 - 3); + t_scale = mvv2e / dof_boltz; + p_scale = 1.0 / 3 / param->xprd / param->yprd / param->zprd; + e_scale = 0.5; + + printf("step\ttemp\t\tpressure\n"); +} + +void computeThermo(int iflag, Parameter *param, Atom *atom) +{ + double t = 0.0, p; + double* vx = atom->vx; + double* vy = atom->vy; + double* vz = atom->vz; + + for(int i = 0; i < atom->Nlocal; i++) { + t += (vx[i] * vx[i] + vy[i] * vy[i] + vz[i] * vz[i]) * param->mass; + } + + t = t * t_scale; + p = (t * dof_boltz) * p_scale; + + int istep = iflag; + + if(iflag == -1){ + istep = param->ntimes; + } + if(iflag == 0){ + mstat = 0; + } + + steparr[mstat] = istep; + tmparr[mstat] = t; + prsarr[mstat] = p; + mstat++; + fprintf(stdout, "%i\t%e\t%e\n", istep, t, p); +} + +void adjustThermo(Parameter *param, Atom *atom) +{ + /* zero center-of-mass motion */ + double vxtot = 0.0; double vytot = 0.0; double vztot = 0.0; + double* vx = atom->vx; double* vy = atom->vy; double* vz = atom->vz; + + for(int i = 0; i < atom->Nlocal; i++) { + vxtot += vx[i]; + vytot += vy[i]; + vztot += vz[i]; + } + + vxtot = vxtot / atom->Natoms; + vytot = vytot / atom->Natoms; + vztot = vztot / atom->Natoms; + + for(int i = 0; i < atom->Nlocal; i++) { + vx[i] -= vxtot; + vy[i] -= vytot; + vz[i] -= vztot; + } + + t_act = 0; + double t = 0.0; + + for(int i = 0; i < atom->Nlocal; i++) { + t += (vx[i] * vx[i] + vy[i] * vy[i] + vz[i] * vz[i]) * param->mass; + } + + t *= t_scale; + double factor = sqrt(param->temp / t); + + for(int i = 0; i < atom->Nlocal; i++) { + vx[i] *= factor; + vy[i] *= factor; + vz[i] *= factor; + } +} diff --git a/src/util.c b/src/util.c new file mode 100644 index 0000000..6a0aa58 --- /dev/null +++ b/src/util.c @@ -0,0 +1,46 @@ +/* + * ======================================================================================= + * + * Author: Jan Eitzinger (je), jan.eitzinger@fau.de + * Copyright (c) 2020 RRZE, University Erlangen-Nuremberg + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in all + * copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + * + * ======================================================================================= + */ +#include + +/* Park/Miller RNG w/out MASKING, so as to be like f90s version */ +#define IA 16807 +#define IM 2147483647 +#define AM (1.0/IM) +#define IQ 127773 +#define IR 2836 +#define MASK 123459876 + +double myrandom(int* idum) +{ + int k= (*idum) / IQ; + double ans; + + *idum = IA * (*idum - k * IQ) - IR * k; + if(*idum < 0) *idum += IM; + ans = AM * (*idum); + return ans; +} diff --git a/util/mdBench.c b/util/mdBench.c new file mode 100644 index 0000000..8138ef6 --- /dev/null +++ b/util/mdBench.c @@ -0,0 +1,1000 @@ +/* + * ======================================================================================= + * + * Author: Jan Eitzinger (je), jan.eitzinger@fau.de + * Copyright (c) 2019 RRZE, University Erlangen-Nuremberg + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in all + * copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + * + * ======================================================================================= + */ + +#include +#include +#include +#include +#include +#include +#include + +#define HLINE "----------------------------------------------------------------------------\n" + +#define FACTOR 0.999 +#define SMALL 1.0e-6 +#define DELTA 20000 + +#ifndef MIN +#define MIN(x,y) ((x)<(y)?(x):(y)) +#endif +#ifndef MAX +#define MAX(x,y) ((x)>(y)?(x):(y)) +#endif +#ifndef ABS +#define ABS(a) ((a) >= 0 ? (a) : -(a)) +#endif + +static int Natoms, Nlocal, Nghost, Nmax; +static double Cutneigh; // neighbor cutoff +static double xprd, yprd, zprd; +static double xlo, xhi; +static double ylo, yhi; +static double zlo, zhi; +static double *x, *y, *z; +static double *vx, *vy, *vz; +static double *fx, *fy, *fz; + +static int NmaxGhost; +static int *BorderMap; +static int *PBCx, *PBCy, *PBCz; + +typedef struct { + int* numneigh; + int* neighbors; + int maxneighs; + int nbinx, nbiny, nbinz; + /* double cutneigh; // neighbor cutoff */ + double cutneighsq; // neighbor cutoff squared + int every; + int ncalls; + int max_totalneigh; + int *bincount; + int *bins; + int nmax; + int nstencil; // # of bins in stencil + int* stencil; // stencil list of bin offsets + int mbins; //total number of bins + int atoms_per_bin; // max atoms per bin + int mbinx, mbiny, mbinz; // n bins in x, y, z + int mbinxlo, mbinylo, mbinzlo; + double binsizex, binsizey, binsizez; + double bininvx, bininvy, bininvz; +} Neighbor; + +typedef struct { + double epsilon; + double sigma6; + double temp; + double rho; + double mass; + int ntimes; + int nstat; + double dt; + double dtforce; + double cutforce; + int nx, ny, nz; +} Parameter; + +typedef struct { + int *steparr; + double *tmparr; + double *engarr; + double *prsarr; + double mvv2e; + int dof_boltz; + double t_scale; + double p_scale; + double e_scale; + double t_act; + double p_act; + double e_act; + int mstat; +} Thermo; + +/* Park/Miller RNG w/out MASKING, so as to be like f90s version */ +#define IA 16807 +#define IM 2147483647 +#define AM (1.0/IM) +#define IQ 127773 +#define IR 2836 +#define MASK 123459876 + +double myrandom(int* idum) +{ + int k= (*idum) / IQ; + double ans; + + *idum = IA * (*idum - k * IQ) - IR * k; + if(*idum < 0) *idum += IM; + ans = AM * (*idum); + return ans; +} + +int coord2bin(Neighbor* neighbor, double xin, double yin, double zin) +{ + int ix, iy, iz; + double bininvx = neighbor->bininvx; + double bininvy = neighbor->bininvy; + double bininvz = neighbor->bininvz; + int mbinxlo = neighbor->mbinxlo; + int mbinylo = neighbor->mbinylo; + int mbinzlo = neighbor->mbinzlo; + + if(xin >= xprd) { + ix = (int)((xin - xprd) * bininvx) + neighbor->nbinx - mbinxlo; + } else if(xin >= 0.0) { + ix = (int)(xin * bininvx) - mbinxlo; + } else { + ix = (int)(xin * bininvx) - mbinxlo - 1; + } + + if(yin >= yprd) { + iy = (int)((yin - yprd) * bininvy) + neighbor->nbiny - mbinylo; + } else if(yin >= 0.0) { + iy = (int)(yin * bininvy) - mbinylo; + } else { + iy = (int)(yin * bininvy) - mbinylo - 1; + } + + if(zin >= zprd) { + iz = (int)((zin - zprd) * bininvz) + neighbor->nbinz - mbinzlo; + } else if(zin >= 0.0) { + iz = (int)(zin * bininvz) - mbinzlo; + } else { + iz = (int)(zin * bininvz) - mbinzlo - 1; + } + + return (iz * neighbor->mbiny * neighbor->mbinx + iy * neighbor->mbinx + ix + 1); +} + +void binatoms(Neighbor *neighbor) +{ + int* bincount = neighbor->bincount; + int mbins = neighbor->mbins; + int nall = Nlocal + Nghost; + int resize = 1; + + while(resize > 0) { + resize = 0; + + for(int i = 0; i < mbins; i++) { + bincount[i] = 0; + } + + for(int i = 0; i < nall; i++) { + int ibin = coord2bin(neighbor, x[i], y[i], z[i]); + + if(bincount[ibin] < neighbor->atoms_per_bin) { + int ac = neighbor->bincount[ibin]++; + neighbor->bins[ibin * neighbor->atoms_per_bin + ac] = i; + } else { + resize = 1; + } + } + + if(resize) { + free(neighbor->bins); + neighbor->atoms_per_bin *= 2; + neighbor->bins = (int*) malloc(mbins * neighbor->atoms_per_bin * sizeof(int)); + } + } +} + +double bindist(Neighbor *neighbor, int i, int j, int k) +{ + double delx, dely, delz; + + if(i > 0) { + delx = (i - 1) * neighbor->binsizex; + } else if(i == 0) { + delx = 0.0; + } else { + delx = (i + 1) * neighbor->binsizex; + } + + if(j > 0) { + dely = (j - 1) * neighbor->binsizey; + } else if(j == 0) { + dely = 0.0; + } else { + dely = (j + 1) * neighbor->binsizey; + } + + if(k > 0) { + delz = (k - 1) * neighbor->binsizez; + } else if(k == 0) { + delz = 0.0; + } else { + delz = (k + 1) * neighbor->binsizez; + } + + return (delx * delx + dely * dely + delz * delz); +} + +void buildNeighborlist(Neighbor *neighbor) +{ + neighbor->ncalls++; + int nall = Nlocal + Nghost; + + /* extend atom arrays if necessary */ + if(nall > neighbor->nmax) { + neighbor->nmax = nall; + if(neighbor->numneigh) free(neighbor->numneigh); + if(neighbor->neighbors) free(neighbor->neighbors); + neighbor->numneigh = (int*) malloc(neighbor->nmax * sizeof(int)); + neighbor->neighbors = (int*) malloc(neighbor->nmax * neighbor->maxneighs * sizeof(int*)); + } + + /* bin local & ghost atoms */ + binatoms(neighbor); + int resize = 1; + + /* loop over each atom, storing neighbors */ + while(resize) { + int new_maxneighs = neighbor->maxneighs; + resize = 0; + + for(int i = 0; i < Nlocal; i++) { + int* neighptr = &neighbor->neighbors[i * neighbor->maxneighs]; + int n = 0; + double xtmp = x[i]; + double ytmp = y[i]; + double ztmp = z[i]; + int ibin = coord2bin(neighbor, xtmp, ytmp, ztmp); + + for(int k = 0; k < neighbor->nstencil; k++) { + int jbin = ibin + neighbor->stencil[k]; + int* loc_bin = &neighbor->bins[jbin * neighbor->atoms_per_bin]; + + for(int m = 0; m < neighbor->bincount[jbin]; m++) { + int j = loc_bin[m]; + + if ( j == i ){ + continue; + } + + double delx = xtmp - x[j]; + double dely = ytmp - y[j]; + double delz = ztmp - z[j]; + double rsq = delx * delx + dely * dely + delz * delz; + + if( rsq <= neighbor->cutneighsq ) { + neighptr[n++] = j; + } + } + } + + neighbor->numneigh[i] = n; + + if(n >= neighbor->maxneighs) { + resize = 1; + + if(n >= new_maxneighs) { + new_maxneighs = n; + } + } + } + + if(resize) { + neighbor->maxneighs = new_maxneighs * 1.2; + free(neighbor->neighbors); + neighbor->neighbors = (int*) malloc(Nmax* neighbor->maxneighs * sizeof(int)); + } + } +} + +void init(Neighbor *neighbor, Parameter *param) +{ + x = NULL; y = NULL; z = NULL; + vx = NULL; vy = NULL; vz = NULL; + fx = NULL; fy = NULL; fz = NULL; + + NmaxGhost = 0; + BorderMap = NULL; + PBCx = NULL; PBCy = NULL; PBCz = NULL; + + param->epsilon = 1.0; + param->sigma6 = 1.0; + param->rho = 0.8442; + param->ntimes = 200; + param->dt = 0.005; + param->nx = 32; + param->ny = 32; + param->nz = 32; + param->cutforce = 2.5; + param->temp = 1.44; + param->nstat = 100; + param->mass = 1.0; + param->dtforce = 0.5 * param->dt; + + Cutneigh = param->cutforce + 0.30; + double neighscale = 5.0 / 6.0; + neighbor->nbinx = neighscale * param->nx; + neighbor->nbiny = neighscale * param->ny; + neighbor->nbinz = neighscale * param->nz; + neighbor->every = 20; + neighbor->ncalls = 0; + neighbor->nmax = 0; + neighbor->atoms_per_bin = 8; + neighbor->maxneighs = 100; + /* neighbor->cutneigh = param->cutforce + 0.30; */ + neighbor->numneigh = NULL; + neighbor->neighbors = NULL; + neighbor->stencil = NULL; + neighbor->bins = NULL; + neighbor->bincount = NULL; +} + +void setup(Neighbor *neighbor, Parameter *param) +{ + double lattice = pow((4.0 / param->rho), (1.0 / 3.0)); + double coord; + int mbinxhi, mbinyhi, mbinzhi; + int nextx, nexty, nextz; + + xprd = param->nx * lattice; + yprd = param->ny * lattice; + zprd = param->nz * lattice; + + xlo = 0.0; xhi = xprd; + ylo = 0.0; yhi = yprd; + zlo = 0.0; zhi = zprd; + + neighbor->cutneighsq = Cutneigh * Cutneigh; + neighbor->binsizex = xprd / neighbor->nbinx; + neighbor->binsizey = yprd / neighbor->nbiny; + neighbor->binsizez = zprd / neighbor->nbinz; + + neighbor->bininvx = 1.0 / neighbor->binsizex; + neighbor->bininvy = 1.0 / neighbor->binsizey; + neighbor->bininvz = 1.0 / neighbor->binsizez; + + coord = xlo - Cutneigh - SMALL * xprd; + neighbor->mbinxlo = (int) (coord * neighbor->bininvx); + if (coord < 0.0) { + neighbor->mbinxlo = neighbor->mbinxlo - 1; + } + coord = xhi + Cutneigh + SMALL * xprd; + mbinxhi = (int) (coord * neighbor->bininvx); + + coord = ylo - Cutneigh - SMALL * yprd; + neighbor->mbinylo = (int) (coord * neighbor->bininvy); + if (coord < 0.0) { + neighbor->mbinylo = neighbor->mbinylo - 1; + } + coord = yhi + Cutneigh + SMALL * yprd; + mbinyhi = (int) (coord * neighbor->bininvy); + + coord = zlo - Cutneigh - SMALL * zprd; + neighbor->mbinzlo = (int) (coord * neighbor->bininvz); + if (coord < 0.0) { + neighbor->mbinzlo = neighbor->mbinzlo - 1; + } + coord = zhi + Cutneigh + SMALL * zprd; + mbinzhi = (int) (coord * neighbor->bininvz); + + neighbor->mbinxlo = neighbor->mbinxlo - 1; + mbinxhi = mbinxhi + 1; + neighbor->mbinx = mbinxhi - neighbor->mbinxlo + 1; + + neighbor->mbinylo = neighbor->mbinylo - 1; + mbinyhi = mbinyhi + 1; + neighbor->mbiny = mbinyhi - neighbor->mbinylo + 1; + + neighbor->mbinzlo = neighbor->mbinzlo - 1; + mbinzhi = mbinzhi + 1; + neighbor->mbinz = mbinzhi - neighbor->mbinzlo + 1; + + nextx = (int) (Cutneigh * neighbor->bininvx); + if(nextx * neighbor->binsizex < FACTOR * Cutneigh) nextx++; + + nexty = (int) (Cutneigh * neighbor->bininvy); + if(nexty * neighbor->binsizey < FACTOR * Cutneigh) nexty++; + + nextz = (int) (Cutneigh * neighbor->bininvz); + if(nextz * neighbor->binsizez < FACTOR * Cutneigh) nextz++; + + if (neighbor->stencil) { + free(neighbor->stencil); + } + + neighbor->stencil = (int*) malloc( + (2 * nextz + 1) * (2 * nexty + 1) * (2 * nextx + 1) * sizeof(int)); + + neighbor->nstencil = 0; + int kstart = -nextz; + + for(int k = kstart; k <= nextz; k++) { + for(int j = -nexty; j <= nexty; j++) { + for(int i = -nextx; i <= nextx; i++) { + if(bindist(neighbor, i, j, k) < neighbor->cutneighsq) { + neighbor->stencil[neighbor->nstencil++] = + k * neighbor->mbiny * neighbor->mbinx + j * neighbor->mbinx + i; + } + } + } + } + + neighbor->mbins = neighbor->mbinx * neighbor->mbiny * neighbor->mbinz; + + if (neighbor->bincount) { + free(neighbor->bincount); + } + neighbor->bincount = (int*) malloc(neighbor->mbins * sizeof(int)); + + if (neighbor->bins) { + free(neighbor->bins); + } + neighbor->bins = (int*) malloc(neighbor->mbins * neighbor->atoms_per_bin * sizeof(int)); +} + +double* myrealloc(double *ptr, int n, int nold) { + + double* newarray; + newarray = (double*) malloc(n * sizeof(double)); + + if(nold) { + memcpy(newarray, ptr, nold * sizeof(double)); + } + if(ptr) { + free(ptr); + } + + return newarray; +} + +int* myreallocInt(int *ptr, int n, int nold) { + + int* newarray; + + newarray = (int*) malloc(n * sizeof(int)); + + if(nold) { + memcpy(newarray, ptr, nold * sizeof(int)); + } + if(ptr) { + free(ptr); + } + + return newarray; +} + +void growBoundary() +{ + int nold = NmaxGhost; + NmaxGhost += DELTA; + + BorderMap = myreallocInt(BorderMap, NmaxGhost, nold); + PBCx = myreallocInt(PBCx, NmaxGhost, nold); + PBCy = myreallocInt(PBCy, NmaxGhost, nold); + PBCz = myreallocInt(PBCz, NmaxGhost, nold); + + if(BorderMap == NULL || PBCx == NULL || PBCy == NULL || PBCz == NULL ) { + printf("ERROR: No memory for Boundary\n"); + } +} + +void growarray() +{ + int nold = Nmax; + Nmax += DELTA; + + x = myrealloc(x, Nmax, nold); y = myrealloc(y, Nmax, nold); z = myrealloc(z, Nmax, nold); + vx = myrealloc(vx, Nmax, nold); vy = myrealloc(vy, Nmax, nold); vz = myrealloc(vz, Nmax, nold); + fx = myrealloc(fx, Nmax, nold); fy = myrealloc(fy, Nmax, nold); fz = myrealloc(fz, Nmax, nold); + + if(x == NULL || y == NULL || z == NULL || + vx == NULL || vy == NULL || vz == NULL || + fx == NULL || fy == NULL || fz == NULL ) { + printf("ERROR: No memory for atoms\n"); + } +} + +void updateBorders() +{ + for(int i = 0; i < Nghost; i++) { + x[Nlocal + i] = x[BorderMap[i]] + PBCx[i] * xprd; + y[Nlocal + i] = y[BorderMap[i]] + PBCy[i] * yprd; + z[Nlocal + i] = z[BorderMap[i]] + PBCz[i] * zprd; + } +} + +void updateAtomLocations() +{ + for(int i = 0; i < Nlocal; i++) { + + if(x[i] < 0.0) { + x[i] += xprd; + } else if(x[i] >= xprd) { + x[i] -= xprd; + } + + if(y[i] < 0.0) { + y[i] += yprd; + } else if(y[i] >= yprd) { + y[i] -= yprd; + } + + if(z[i] < 0.0) { + z[i] += zprd; + } else if(z[i] >= zprd) { + z[i] -= zprd; + } + } +} + +void setupBordersNew() +{ + int lastidx = 0; + int nghostprev = 0; + Nghost = 0; + + for (int i = 0; i < Nlocal; i++) { + + if (Nlocal + Nghost + 1 >= Nmax) { + growarray(); + } + + if (x[i] < Cutneigh) { + Nghost++; + x[i+lastidx] = x[i] + xprd; + y[i+lastidx] = y[i]; + z[i+lastidx] = z[i]; + lastidx++; + } else if (x[i] >= xprd - Cutneigh) { + Nghost++; + x[i+lastidx] = x[i] - xprd; + y[i+lastidx] = y[i]; + z[i+lastidx] = z[i]; + lastidx++; + } + } + + nghostprev = Nghost+1; + + for (int i = 0; i < Nlocal + nghostprev ; i++) { + + if (Nlocal + Nghost + 1 >= Nmax) { + growarray(); + } + + if (y[i] < Cutneigh) { + Nghost++; + x[i+lastidx] = x[i]; + y[i+lastidx] = y[i] + yprd; + z[i+lastidx] = z[i]; + lastidx++; + } else if (y[i] >= yprd - Cutneigh) { + Nghost++; + x[i+lastidx] = x[i]; + y[i+lastidx] = y[i] - yprd; + z[i+lastidx] = z[i]; + lastidx++; + } + } + + nghostprev = Nghost+1; + + for (int i = 0; i < Nlocal + nghostprev; i++) { + + if (Nlocal + Nghost + 1 >= Nmax) { + growarray(); + } + + if (z[i] < Cutneigh) { + Nghost++; + x[i+lastidx] = x[i]; + y[i+lastidx] = y[i]; + z[i+lastidx] = z[i] + zprd; + lastidx++; + } else if(z[i] >= zprd - Cutneigh) { + Nghost++; + x[i+lastidx] = x[i]; + y[i+lastidx] = y[i]; + z[i+lastidx] = z[i] - zprd; + lastidx++; + } + } + + Nghost++; +} + +#define ADDGHOST(dx,dy,dz) Nghost++; BorderMap[Nghost] = i; PBCx[Nghost] = dx; PBCy[Nghost] = dy; PBCz[Nghost] = dz; +void setupBorders() +{ + Nghost = -1; + + for(int i = 0; i < Nlocal; i++) { + + if (Nlocal + Nghost + 7 >= Nmax) { + growarray(); + } + if (Nghost + 7 >= NmaxGhost) { + growBoundary(); + } + + /* Setup ghost atoms */ + /* 6 planes */ + if (x[i] < Cutneigh) { ADDGHOST(+1,0,0); } + if (x[i] >= (xprd-Cutneigh)) { ADDGHOST(-1,0,0); } + if (y[i] < Cutneigh) { ADDGHOST(0,+1,0); } + if (y[i] >= (yprd-Cutneigh)) { ADDGHOST(0,-1,0); } + if (z[i] < Cutneigh) { ADDGHOST(0,0,+1); } + if (z[i] >= (zprd-Cutneigh)) { ADDGHOST(0,0,-1); } + /* 8 corners */ + if (x[i] < Cutneigh && y[i] < Cutneigh && z[i] < Cutneigh) { ADDGHOST(+1,+1,+1); } + if (x[i] < Cutneigh && y[i] >= (yprd-Cutneigh) && z[i] < Cutneigh) { ADDGHOST(+1,-1,+1); } + if (x[i] < Cutneigh && y[i] >= Cutneigh && z[i] >= (zprd-Cutneigh)) { ADDGHOST(+1,+1,-1); } + if (x[i] < Cutneigh && y[i] >= (yprd-Cutneigh) && z[i] >= (zprd-Cutneigh)) { ADDGHOST(+1,-1,-1); } + if (x[i] >= (xprd-Cutneigh) && y[i] < Cutneigh && z[i] < Cutneigh) { ADDGHOST(-1,+1,+1); } + if (x[i] >= (xprd-Cutneigh) && y[i] >= (yprd-Cutneigh) && z[i] < Cutneigh) { ADDGHOST(-1,-1,+1); } + if (x[i] >= (xprd-Cutneigh) && y[i] < Cutneigh && z[i] >= (zprd-Cutneigh)) { ADDGHOST(-1,+1,-1); } + if (x[i] >= (xprd-Cutneigh) && y[i] >= (yprd-Cutneigh) && z[i] >= (zprd-Cutneigh)) { ADDGHOST(-1,-1,-1); } + /* 12 edges */ + if (x[i] < Cutneigh && z[i] < Cutneigh) { ADDGHOST(+1,0,+1); } + if (x[i] < Cutneigh && z[i] >= (zprd-Cutneigh)) { ADDGHOST(+1,0,-1); } + if (x[i] >= (xprd-Cutneigh) && z[i] < Cutneigh) { ADDGHOST(-1,0,+1); } + if (x[i] >= (xprd-Cutneigh) && z[i] >= (zprd-Cutneigh)) { ADDGHOST(-1,0,-1); } + if (y[i] < Cutneigh && z[i] < Cutneigh) { ADDGHOST(0,+1,+1); } + if (y[i] < Cutneigh && z[i] >= (zprd-Cutneigh)) { ADDGHOST(0,+1,-1); } + if (y[i] >= (yprd-Cutneigh) && z[i] < Cutneigh) { ADDGHOST(0,-1,+1); } + if (y[i] >= (yprd-Cutneigh) && z[i] >= (zprd-Cutneigh)) { ADDGHOST(0,-1,-1); } + if (y[i] < Cutneigh && x[i] < Cutneigh) { ADDGHOST(+1,+1,0); } + if (y[i] < Cutneigh && x[i] >= (xprd-Cutneigh)) { ADDGHOST(-1,+1,0); } + if (y[i] >= (yprd-Cutneigh) && x[i] < Cutneigh) { ADDGHOST(+1,-1,0); } + if (y[i] >= (yprd-Cutneigh) && x[i] >= (xprd-Cutneigh)) { ADDGHOST(-1,-1,0); } + } + // increase by one to make it the ghost atom count + Nghost++; +} + +void sortAtoms(Neighbor *neighbor) +{ + binatoms(neighbor); + int* binpos = neighbor->bincount; + int* bins = neighbor->bins; + + int mbins = neighbor->mbins; + int atoms_per_bin = neighbor->atoms_per_bin; + + for(int i=1; i0?binpos[mybin-1]:0; + int count = binpos[mybin] - start; + + for(int k=0; knx * param->ny * param->nz; + Nlocal = 0; + double alat = pow((4.0 / param->rho), (1.0 / 3.0)); + int ilo = (int) (xlo / (0.5 * alat) - 1); + int ihi = (int) (xhi / (0.5 * alat) + 1); + int jlo = (int) (ylo / (0.5 * alat) - 1); + int jhi = (int) (yhi / (0.5 * alat) + 1); + int klo = (int) (zlo / (0.5 * alat) - 1); + int khi = (int) (zhi / (0.5 * alat) + 1); + + ilo = MAX(ilo, 0); + ihi = MIN(ihi, 2 * param->nx - 1); + jlo = MAX(jlo, 0); + jhi = MIN(jhi, 2 * param->ny - 1); + klo = MAX(klo, 0); + khi = MIN(khi, 2 * param->nz - 1); + + double xtmp, ytmp, ztmp, vxtmp, vytmp, vztmp; + int i, j, k, m, n; + int sx = 0; int sy = 0; int sz = 0; + int ox = 0; int oy = 0; int oz = 0; + int subboxdim = 8; + + while(oz * subboxdim <= khi) { + + k = oz * subboxdim + sz; + j = oy * subboxdim + sy; + i = ox * subboxdim + sx; + + if(((i + j + k) % 2 == 0) && + (i >= ilo) && (i <= ihi) && + (j >= jlo) && (j <= jhi) && + (k >= klo) && (k <= khi)) { + + xtmp = 0.5 * alat * i; + ytmp = 0.5 * alat * j; + ztmp = 0.5 * alat * k; + + if( xtmp >= xlo && xtmp < xhi && + ytmp >= ylo && ytmp < yhi && + ztmp >= zlo && ztmp < zhi ) { + + n = k * (2 * param->ny) * (2 * param->nx) + + j * (2 * param->nx) + + i + 1; + + for(m = 0; m < 5; m++) { + myrandom(&n); + } + vxtmp = myrandom(&n); + + for(m = 0; m < 5; m++){ + myrandom(&n); + } + vytmp = myrandom(&n); + + for(m = 0; m < 5; m++) { + myrandom(&n); + } + vztmp = myrandom(&n); + + if(Nlocal == Nmax) { + growarray(); + } + + x[Nlocal] = xtmp; y[Nlocal] = ytmp; z[Nlocal] = ztmp; + vx[Nlocal] = vxtmp; vy[Nlocal] = vytmp; vz[Nlocal] = vztmp; + Nlocal++; + } + } + + sx++; + + if(sx == subboxdim) { sx = 0; sy++; } + if(sy == subboxdim) { sy = 0; sz++; } + if(sz == subboxdim) { sz = 0; ox++; } + if(ox * subboxdim > ihi) { ox = 0; oy++; } + if(oy * subboxdim > jhi) { oy = 0; oz++; } + } +} + +void adjustVelocity(Parameter *param, Thermo *thermo) +{ + /* zero center-of-mass motion */ + double vxtot = 0.0; + double vytot = 0.0; + double vztot = 0.0; + + for(int i = 0; i < Nlocal; i++) { + vxtot += vx[i]; + vytot += vy[i]; + vztot += vz[i]; + } + + vxtot = vxtot / Natoms; + vytot = vytot / Natoms; + vztot = vztot / Natoms; + + for(int i = 0; i < Nlocal; i++) { + vx[i] -= vxtot; + vy[i] -= vytot; + vz[i] -= vztot; + } + + thermo->t_act = 0; + double t = 0.0; + + for(int i = 0; i < Nlocal; i++) { + t += (vx[i] * vx[i] + vy[i] * vy[i] + vz[i] * vz[i]) * param->mass; + } + + t *= thermo->t_scale; + double factor = sqrt(param->temp / t); + + for(int i = 0; i < Nlocal; i++) { + vx[i] *= factor; + vy[i] *= factor; + vz[i] *= factor; + } +} + +void thermoSetup(Parameter *param, Thermo *thermo) +{ + int maxstat = param->ntimes / param->nstat + 2; + + thermo->steparr = (int*) malloc(maxstat * sizeof(int)); + thermo->tmparr = (double*) malloc(maxstat * sizeof(double)); + thermo->engarr = (double*) malloc(maxstat * sizeof(double)); + thermo->prsarr = (double*) malloc(maxstat * sizeof(double)); + + thermo->mvv2e = 1.0; + thermo->dof_boltz = (Natoms * 3 - 3); + thermo->t_scale = thermo->mvv2e / thermo->dof_boltz; + thermo->p_scale = 1.0 / 3 / xprd / yprd / zprd; + thermo->e_scale = 0.5; + + printf("step\ttemp\t\tpressure\n"); +} + + +void thermoCompute(int iflag, Parameter *param, Thermo *thermo) +{ + double t = 0.0, p; + + for(int i = 0; i < Nlocal; i++) { + t += (vx[i] * vx[i] + vy[i] * vy[i] + vz[i] * vz[i]) * param->mass; + } + + t = t * thermo->t_scale; + p = (t * thermo->dof_boltz) * thermo->p_scale; + + int istep = iflag; + + if(iflag == -1){ + istep = param->ntimes; + } + if(iflag == 0){ + thermo->mstat = 0; + } + + thermo->steparr[thermo->mstat] = istep; + thermo->tmparr[thermo->mstat] = t; + thermo->prsarr[thermo->mstat] = p; + thermo->mstat++; + fprintf(stdout, "%i\t%e\t%e\n", istep, t, p); +} + +void initialIntegrate(Parameter *param) +{ + for(int i = 0; i < Nlocal; i++) { + vx[i] += param->dtforce * fx[i]; + vy[i] += param->dtforce * fy[i]; + vz[i] += param->dtforce * fz[i]; + x[i] += param->dt * vx[i]; + y[i] += param->dt * vy[i]; + z[i] += param->dt * vz[i]; + } +} + +void finalIntegrate(Parameter *param) +{ + for(int i = 0; i < Nlocal; i++) { + vx[i] += param->dtforce * fx[i]; + vy[i] += param->dtforce * fy[i]; + vz[i] += param->dtforce * fz[i]; + } +} + +void computeForce(Neighbor *neighbor, Parameter *param) +{ + int* neighs; + double cutforcesq = param->cutforce * param->cutforce; + double sigma6 = param->sigma6; + double epsilon = param->epsilon; + + for(int i = 0; i < Nlocal; i++) { + fx[i] = 0.0; + fy[i] = 0.0; + fz[i] = 0.0; + } + + for(int i = 0; i < Nlocal; i++) { + neighs = &neighbor->neighbors[i * neighbor->maxneighs]; + int numneighs = neighbor->numneigh[i]; + double xtmp = x[i]; + double ytmp = y[i]; + double ztmp = z[i]; + + double fix = 0; + double fiy = 0; + double fiz = 0; + + for(int k = 0; k < numneighs; k++) { + int j = neighs[k]; + double delx = xtmp - x[j]; + double dely = ytmp - y[j]; + double delz = ztmp - z[j]; + double rsq = delx * delx + dely * dely + delz * delz; + + if(rsq < cutforcesq) { + double sr2 = 1.0 / rsq; + double sr6 = sr2 * sr2 * sr2 * sigma6; + double force = 48.0 * sr6 * (sr6 - 0.5) * sr2 * epsilon; + fix += delx * force; + fiy += dely * force; + fiz += delz * force; + } + } + + fx[i] += fix; + fy[i] += fiy; + fz[i] += fiz; + } +} + +int main (int argc, char** argv) +{ + Neighbor neighbor; + Parameter param; + Thermo thermo; + + init(&neighbor, ¶m); + setup(&neighbor, ¶m); + create_atoms(¶m); + thermoSetup(¶m, &thermo); + adjustVelocity(¶m, &thermo); + setupBorders(); + updateBorders(); + buildNeighborlist(&neighbor); + thermoCompute(0, ¶m, &thermo); + computeForce(&neighbor, ¶m); + + for(int n = 0; n < param.ntimes; n++) { + + initialIntegrate(¶m); + + if((n + 1) % neighbor.every) { + updateBorders(); + } else { + updateAtomLocations(); + setupBorders(); + updateBorders(); + /* sortAtoms(&neighbor); */ + buildNeighborlist(&neighbor); + } + + computeForce(&neighbor, ¶m); + finalIntegrate(¶m); + + if(!((n + 1) % param.nstat) && (n+1) < param.ntimes) { + thermoCompute(n + 1, ¶m, &thermo); + } + } + + thermoCompute(-1, ¶m, &thermo); + + return EXIT_SUCCESS; +}