diff --git a/gromacs/force_lj.c b/gromacs/force_lj.c
index fd69675..46cdd94 100644
--- a/gromacs/force_lj.c
+++ b/gromacs/force_lj.c
@@ -20,8 +20,9 @@
* with MD-Bench. If not, see .
* =======================================================================================
*/
-#include
+#include
+#include
#include
#include
#include
@@ -29,6 +30,7 @@
#include
double computeForceLJ(Parameter *param, Atom *atom, Neighbor *neighbor, Stats *stats) {
+ fprintf(stdout, "computeForceLJ begin\n");
int Nlocal = atom->Nlocal;
int* neighs;
MD_FLOAT cutforcesq = param->cutforce * param->cutforce;
@@ -47,7 +49,7 @@ double computeForceLJ(Parameter *param, Atom *atom, Neighbor *neighbor, Stats *s
double S = getTimeStamp();
LIKWID_MARKER_START("force");
-#pragma omp parallel for
+ #pragma omp parallel for
for(int ci = 0; ci < atom->Nclusters_local; ci++) {
MD_FLOAT *ciptr = cluster_pos_ptr(ci);
MD_FLOAT *cifptr = cluster_force_ptr(ci);
@@ -57,7 +59,7 @@ double computeForceLJ(Parameter *param, Atom *atom, Neighbor *neighbor, Stats *s
for(int k = 0; k < numneighs; k++) {
int cj = neighs[k];
MD_FLOAT *cjptr = cluster_pos_ptr(cj);
- for(int cii = 0; cii < CLUSTER_DIM_M; cii++) {
+ for(int cii = 0; cii < atom->clusters[ci].natoms; cii++) {
MD_FLOAT xtmp = cluster_x(ciptr, cii);
MD_FLOAT ytmp = cluster_y(ciptr, cii);
MD_FLOAT ztmp = cluster_z(ciptr, cii);
@@ -66,17 +68,19 @@ double computeForceLJ(Parameter *param, Atom *atom, Neighbor *neighbor, Stats *s
MD_FLOAT fiz = 0;
for(int cjj = 0; cjj < CLUSTER_DIM_N; cjj++) {
- MD_FLOAT delx = xtmp - cluster_x(cjptr, cjj);
- MD_FLOAT dely = ytmp - cluster_y(cjptr, cjj);
- MD_FLOAT delz = ztmp - cluster_z(cjptr, cjj);
- MD_FLOAT rsq = delx * delx + dely * dely + delz * delz;
- if(rsq < cutforcesq) {
- MD_FLOAT sr2 = 1.0 / rsq;
- MD_FLOAT sr6 = sr2 * sr2 * sr2 * sigma6;
- MD_FLOAT force = 48.0 * sr6 * (sr6 - 0.5) * sr2 * epsilon;
- fix += delx * force;
- fiy += dely * force;
- fiz += delz * force;
+ if(ci != cj || cii != cjj) {
+ MD_FLOAT delx = xtmp - cluster_x(cjptr, cjj);
+ MD_FLOAT dely = ytmp - cluster_y(cjptr, cjj);
+ MD_FLOAT delz = ztmp - cluster_z(cjptr, cjj);
+ MD_FLOAT rsq = delx * delx + dely * dely + delz * delz;
+ if(rsq < cutforcesq) {
+ MD_FLOAT sr2 = 1.0 / rsq;
+ MD_FLOAT sr6 = sr2 * sr2 * sr2 * sigma6;
+ MD_FLOAT force = 48.0 * sr6 * (sr6 - 0.5) * sr2 * epsilon;
+ fix += delx * force;
+ fiy += dely * force;
+ fiz += delz * force;
+ }
}
}
@@ -92,5 +96,6 @@ double computeForceLJ(Parameter *param, Atom *atom, Neighbor *neighbor, Stats *s
LIKWID_MARKER_STOP("force");
double E = getTimeStamp();
+ fprintf(stdout, "computeForceLJ end\n");
return E-S;
}
diff --git a/gromacs/includes/pbc.h b/gromacs/includes/pbc.h
index 727ac80..648efd2 100644
--- a/gromacs/includes/pbc.h
+++ b/gromacs/includes/pbc.h
@@ -26,7 +26,7 @@
#ifndef __PBC_H_
#define __PBC_H_
extern void initPbc();
-extern void updatePbc(Atom*, Parameter*);
+extern void updatePbc(Atom*, Parameter*, int);
extern void updateAtomsPbc(Atom*, Parameter*);
extern void setupPbc(Atom*, Parameter*);
#endif
diff --git a/gromacs/main.c b/gromacs/main.c
index 8a05002..b4c1bd9 100644
--- a/gromacs/main.c
+++ b/gromacs/main.c
@@ -48,8 +48,7 @@
extern double computeForceLJ(Parameter*, Atom*, Neighbor*, Stats*);
extern double computeForceEam(Eam*, Parameter*, Atom*, Neighbor*, Stats*);
-void init(Parameter *param)
-{
+void init(Parameter *param) {
param->input_file = NULL;
param->vtk_file = NULL;
param->force_field = FF_LJ;
@@ -72,13 +71,7 @@ void init(Parameter *param)
param->proc_freq = 2.4;
}
-double setup(
- Parameter *param,
- Eam *eam,
- Atom *atom,
- Neighbor *neighbor,
- Stats *stats)
-{
+double setup(Parameter *param, Eam *eam, Atom *atom, Neighbor *neighbor, Stats *stats) {
if(param->force_field == FF_EAM) { initEam(eam, param); }
double S, E;
param->lattice = pow((4.0 / param->rho), (1.0 / 3.0));
@@ -91,11 +84,13 @@ double setup(
initPbc(atom);
initStats(stats);
initNeighbor(neighbor, param);
+
if(param->input_file == NULL) {
createAtom(atom, param);
} else {
readAtom(atom, param);
}
+
setupNeighbor(param, atom);
setupThermo(param, atom->Natoms);
if(param->input_file == NULL) { adjustThermo(param, atom); }
@@ -104,17 +99,11 @@ double setup(
binClusters(atom);
buildNeighbor(atom, neighbor);
E = getTimeStamp();
-
return E-S;
}
-double reneighbour(
- Parameter *param,
- Atom *atom,
- Neighbor *neighbor)
-{
+double reneighbour(Parameter *param, Atom *atom, Neighbor *neighbor) {
double S, E;
-
S = getTimeStamp();
LIKWID_MARKER_START("reneighbour");
updateSingleAtoms(atom);
@@ -125,11 +114,11 @@ double reneighbour(
buildNeighbor(atom, neighbor);
LIKWID_MARKER_STOP("reneighbour");
E = getTimeStamp();
-
return E-S;
}
void initialIntegrate(Parameter *param, Atom *atom) {
+ fprintf(stdout, "initialIntegrate start\n");
for(int ci = 0; ci < atom->Nclusters_local; ci++) {
MD_FLOAT *ciptr = cluster_pos_ptr(ci);
MD_FLOAT *civptr = cluster_velocity_ptr(ci);
@@ -144,9 +133,11 @@ void initialIntegrate(Parameter *param, Atom *atom) {
cluster_z(ciptr, cii) += param->dt * cluster_z(civptr, cii);
}
}
+ fprintf(stdout, "initialIntegrate end\n");
}
void finalIntegrate(Parameter *param, Atom *atom) {
+ fprintf(stdout, "finalIntegrate start\n");
for(int ci = 0; ci < atom->Nclusters_local; ci++) {
MD_FLOAT *civptr = cluster_velocity_ptr(ci);
MD_FLOAT *cifptr = cluster_force_ptr(ci);
@@ -157,10 +148,10 @@ void finalIntegrate(Parameter *param, Atom *atom) {
cluster_z(civptr, cii) += param->dtforce * cluster_z(cifptr, cii);
}
}
+ fprintf(stdout, "finalIntegrate end\n");
}
-void printAtomState(Atom *atom)
-{
+void printAtomState(Atom *atom) {
printf("Atom counts: Natoms=%d Nlocal=%d Nghost=%d Nmax=%d\n",
atom->Natoms, atom->Nlocal, atom->Nghost, atom->Nmax);
@@ -171,8 +162,7 @@ void printAtomState(Atom *atom)
/* } */
}
-int main(int argc, char** argv)
-{
+int main(int argc, char** argv) {
double timer[NUMTIMER];
Eam eam;
Atom atom;
@@ -277,7 +267,7 @@ int main(int argc, char** argv)
initialIntegrate(¶m, &atom);
if((n + 1) % param.every) {
- updatePbc(&atom, ¶m);
+ updatePbc(&atom, ¶m, 0);
} else {
timer[NEIGH] += reneighbour(¶m, &atom, &neighbor);
}
diff --git a/gromacs/neighbor.c b/gromacs/neighbor.c
index 1ae533e..1e65b33 100644
--- a/gromacs/neighbor.c
+++ b/gromacs/neighbor.c
@@ -202,6 +202,7 @@ int atomDistanceInRange(Atom *atom, int ci, int cj, MD_FLOAT rsq) {
}
void buildNeighbor(Atom *atom, Neighbor *neighbor) {
+ printf("buildNeighbor start\n");
int nall = atom->Nclusters_local + atom->Nclusters_ghost;
/* extend atom arrays if necessary */
@@ -213,7 +214,10 @@ void buildNeighbor(Atom *atom, Neighbor *neighbor) {
neighbor->neighbors = (int*) malloc(nmax * neighbor->maxneighs * sizeof(int*));
}
- const MD_FLOAT rBB = cutneighsq / 2.0; // TODO: change this
+ MD_FLOAT bbx = 0.5 * (binsizex + binsizex);
+ MD_FLOAT bby = 0.5 * (binsizey + binsizey);
+ MD_FLOAT rbb_sq = MAX(0.0, cutneigh - 0.5 * sqrt(bbx * bbx + bby * bby));
+ rbb_sq = rbb_sq * rbb_sq;
int resize = 1;
/* loop over each atom, storing neighbors */
@@ -250,7 +254,7 @@ void buildNeighbor(Atom *atom, Neighbor *neighbor) {
double d_bb_sq = getBoundingBoxDistanceSq(atom, ci, cj);
if(d_bb_sq < cutneighsq) {
- if(d_bb_sq < rBB || atomDistanceInRange(atom, ci, cj, cutneighsq)) {
+ if(d_bb_sq < rbb_sq || atomDistanceInRange(atom, ci, cj, cutneighsq)) {
neighptr[n++] = cj;
}
}
@@ -282,6 +286,30 @@ void buildNeighbor(Atom *atom, Neighbor *neighbor) {
neighbor->neighbors = (int*) malloc(atom->Nmax * neighbor->maxneighs * sizeof(int));
}
}
+
+ /*
+ printf("\ncutneighsq = %f, rbb_sq = %f\n", cutneighsq, rbb_sq);
+ for(int ci = 4; ci < 6; ci++) {
+ MD_FLOAT *ciptr = cluster_pos_ptr(ci);
+ int* neighptr = &(neighbor->neighbors[ci * neighbor->maxneighs]);
+ printf("Cluster %d, bbx = {%f, %f}, bby = {%f, %f}, bbz = {%f, %f}\n", ci, atom->clusters[ci].bbminx, atom->clusters[ci].bbmaxx, atom->clusters[ci].bbminy, atom->clusters[ci].bbmaxy, atom->clusters[ci].bbminz, atom->clusters[ci].bbmaxz);
+ for(int cii = 0; cii < CLUSTER_DIM_N; cii++) {
+ fprintf(stdout, "%f, %f, %f\n", cluster_x(ciptr, cii), cluster_y(ciptr, cii), cluster_z(ciptr, cii));
+ }
+
+ printf("Neighbors:\n");
+ for(int k = 0; k < neighbor->numneigh[ci]; k++) {
+ const int cj = neighptr[k];
+ MD_FLOAT *cjptr = cluster_pos_ptr(cj);
+ printf(" Cluster %d, bbx = {%f, %f}, bby = {%f, %f}, bbz = {%f, %f}\n", cj, atom->clusters[cj].bbminx, atom->clusters[cj].bbmaxx, atom->clusters[cj].bbminy, atom->clusters[cj].bbmaxy, atom->clusters[cj].bbminz, atom->clusters[cj].bbmaxz);
+ for(int cjj = 0; cjj < CLUSTER_DIM_N; cjj++) {
+ fprintf(stdout, " %f, %f, %f\n", cluster_x(cjptr, cjj), cluster_y(cjptr, cjj), cluster_z(cjptr, cjj));
+ }
+ }
+ }
+ */
+
+ printf("buildNeighbor end\n");
}
/* internal subroutines */
@@ -348,6 +376,7 @@ void coord2bin2D(MD_FLOAT xin, MD_FLOAT yin, int *ix, int *iy) {
}
void binAtoms(Atom *atom) {
+ printf("binAtoms start\n");
int nall = atom->Nlocal + atom->Nghost;
int resize = 1;
@@ -360,6 +389,7 @@ void binAtoms(Atom *atom) {
for(int i = 0; i < nall; i++) {
int ibin = coord2bin(atom_x(i), atom_y(i));
+ if(ibin < 0 || ibin > mbins) { fprintf(stderr, "%d: %f, %f\n", i, atom_x(i), atom_y(i)); }
if(bincount[ibin] < atoms_per_bin) {
int ac = bincount[ibin]++;
bins[ibin * atoms_per_bin + ac] = i;
@@ -374,10 +404,12 @@ void binAtoms(Atom *atom) {
bins = (int*) malloc(mbins * atoms_per_bin * sizeof(int));
}
}
+ printf("binAtoms end\n");
}
// TODO: Use pigeonhole sorting
void sortAtomsByZCoord(Atom *atom) {
+ printf("sortAtomsByZCoord start\n");
for(int bin = 0; bin < mbins; bin++) {
int c = bincount[bin];
int *bin_ptr = &bins[bin * atoms_per_bin];
@@ -392,7 +424,7 @@ void sortAtomsByZCoord(Atom *atom) {
int j = bin_ptr[ac_j];
MD_FLOAT zj = atom_z(j);
if(zj < min_z) {
- min_ac = zj;
+ min_ac = ac_j;
min_idx = j;
min_z = zj;
}
@@ -402,9 +434,11 @@ void sortAtomsByZCoord(Atom *atom) {
bin_ptr[min_ac] = i;
}
}
+ printf("sortAtomsByZCoord end\n");
}
void buildClusters(Atom *atom) {
+ printf("buildClusters start\n");
atom->Nclusters_local = 0;
/* bin local atoms */
@@ -414,7 +448,7 @@ void buildClusters(Atom *atom) {
for(int bin = 0; bin < mbins; bin++) {
int c = bincount[bin];
int ac = 0;
- const int nclusters = ((c + CLUSTER_DIM_N - 1) / CLUSTER_DIM_N);
+ const int nclusters = ((c + CLUSTER_DIM_M - 1) / CLUSTER_DIM_M);
for(int cl = 0; cl < nclusters; cl++) {
const int ci = atom->Nclusters_local;
@@ -423,6 +457,7 @@ void buildClusters(Atom *atom) {
}
MD_FLOAT *cptr = cluster_pos_ptr(ci);
+ MD_FLOAT *cvptr = cluster_velocity_ptr(ci);
MD_FLOAT bbminx = INFINITY, bbmaxx = -INFINITY;
MD_FLOAT bbminy = INFINITY, bbmaxy = -INFINITY;
MD_FLOAT bbminz = INFINITY, bbmaxz = -INFINITY;
@@ -438,14 +473,17 @@ void buildClusters(Atom *atom) {
cluster_x(cptr, cii) = xtmp;
cluster_y(cptr, cii) = ytmp;
cluster_z(cptr, cii) = ztmp;
+ cluster_x(cvptr, cii) = atom->vx[i];
+ cluster_y(cvptr, cii) = atom->vy[i];
+ cluster_z(cvptr, cii) = atom->vz[i];
// TODO: To create the bounding boxes faster, we can use SIMD operations
if(bbminx > xtmp) { bbminx = xtmp; }
if(bbmaxx < xtmp) { bbmaxx = xtmp; }
if(bbminy > ytmp) { bbminy = ytmp; }
if(bbmaxy < ytmp) { bbmaxy = ytmp; }
- if(bbminz > ytmp) { bbminz = ytmp; }
- if(bbmaxz < ytmp) { bbmaxz = ytmp; }
+ if(bbminz > ztmp) { bbminz = ztmp; }
+ if(bbmaxz < ztmp) { bbmaxz = ztmp; }
atom->clusters[ci].type[cii] = atom->type[i];
atom->clusters[ci].natoms++;
@@ -468,9 +506,23 @@ void buildClusters(Atom *atom) {
atom->Nclusters_local++;
}
}
+
+ /*
+ for(int ci = 4; ci < 9; ci++) {
+ MD_FLOAT *cptr = cluster_pos_ptr(ci);
+ fprintf(stdout, "Cluster %d:\n", ci);
+ fprintf(stdout, "bin=%d, Natoms=%d, bbox={%f,%f},{%f,%f},{%f,%f}\n", atom->clusters[ci].bin, atom->clusters[ci].natoms, atom->clusters[ci].bbminx, atom->clusters[ci].bbmaxx, atom->clusters[ci].bbminy, atom->clusters[ci].bbmaxy, atom->clusters[ci].bbminz, atom->clusters[ci].bbmaxz);
+ for(int cii = 0; cii < CLUSTER_DIM_N; cii++) {
+ fprintf(stdout, "%f, %f, %f\n", cluster_x(cptr, cii), cluster_y(cptr, cii), cluster_z(cptr, cii));
+ }
+ }
+ */
+
+ printf("buildClusters end\n");
}
void binClusters(Atom *atom) {
+ printf("binClusters start\n");
const int nlocal = atom->Nclusters_local;
int resize = 1;
@@ -492,35 +544,37 @@ void binClusters(Atom *atom) {
}
}
- for(int ci = 0; ci < atom->Nclusters_ghost && !resize; ci++) {
- MD_FLOAT *cptr = cluster_pos_ptr(nlocal + ci);
+ for(int cg = 0; cg < atom->Nclusters_ghost && !resize; cg++) {
+ const int ci = nlocal + cg;
+ MD_FLOAT *cptr = cluster_pos_ptr(ci);
MD_FLOAT xtmp, ytmp;
int ix = -1, iy = -1;
xtmp = cluster_x(cptr, 0);
ytmp = cluster_y(cptr, 0);
coord2bin2D(xtmp, ytmp, &ix, &iy);
- ix = MAX(MIN(ix, nbinx - 1), 0);
- iy = MAX(MIN(iy, nbiny - 1), 0);
+ ix = MAX(MIN(ix, mbinx - 1), 0);
+ iy = MAX(MIN(iy, mbiny - 1), 0);
for(int cii = 1; cii < atom->clusters[ci].natoms; cii++) {
int nix, niy;
xtmp = cluster_x(cptr, cii);
ytmp = cluster_y(cptr, cii);
coord2bin2D(xtmp, ytmp, &nix, &niy);
- nix = MAX(MIN(nix, nbinx - 1), 0);
- niy = MAX(MIN(niy, nbiny - 1), 0);
+ nix = MAX(MIN(nix, mbinx - 1), 0);
+ niy = MAX(MIN(niy, mbiny - 1), 0);
// Always put the cluster on the bin of its innermost atom so
// the cluster should be closer to local clusters
- if(atom->PBCx[ci] > 0 && ix > nix) { ix = nix; }
- if(atom->PBCx[ci] < 0 && ix < nix) { ix = nix; }
- if(atom->PBCy[ci] > 0 && iy > niy) { iy = niy; }
- if(atom->PBCy[ci] < 0 && iy < niy) { iy = niy; }
+ if(atom->PBCx[cg] > 0 && ix > nix) { ix = nix; }
+ if(atom->PBCx[cg] < 0 && ix < nix) { ix = nix; }
+ if(atom->PBCy[cg] > 0 && iy > niy) { iy = niy; }
+ if(atom->PBCy[cg] < 0 && iy < niy) { iy = niy; }
}
int bin = iy * mbinx + ix + 1;
int c = cluster_bincount[bin];
if(c < clusters_per_bin) {
- cluster_bins[bin * clusters_per_bin + c] = nlocal + ci;
+ atom->clusters[ci].bin = bin;
+ cluster_bins[bin * clusters_per_bin + c] = ci;
cluster_bincount[bin]++;
} else {
resize = 1;
@@ -533,18 +587,31 @@ void binClusters(Atom *atom) {
cluster_bins = (int*) malloc(mbins * clusters_per_bin * sizeof(int));
}
}
+
+ /*
+ fprintf(stdout, "cluster_bincount\n");
+ for(int i = 0; i < mbins; i++) { fprintf(stdout, "%d, ", cluster_bincount[i]); }
+ fprintf(stdout, "\n");
+ */
+
+ printf("binClusters stop\n");
}
void updateSingleAtoms(Atom *atom) {
+ printf("updateSingleAtoms start\n");
int Natom = 0;
for(int ci = 0; ci < atom->Nclusters_local; ci++) {
MD_FLOAT *cptr = cluster_pos_ptr(ci);
+ MD_FLOAT *cvptr = cluster_velocity_ptr(ci);
for(int cii = 0; cii < atom->clusters[ci].natoms; cii++) {
atom_x(Natom) = cluster_x(cptr, cii);
atom_y(Natom) = cluster_y(cptr, cii);
atom_z(Natom) = cluster_z(cptr, cii);
+ atom->vx[Natom] = cluster_x(cvptr, cii);
+ atom->vy[Natom] = cluster_y(cvptr, cii);
+ atom->vz[Natom] = cluster_z(cvptr, cii);
Natom++;
}
}
@@ -552,4 +619,5 @@ void updateSingleAtoms(Atom *atom) {
if(Natom != atom->Nlocal) {
fprintf(stderr, "updateSingleAtoms(): Number of atoms changed!\n");
}
+ printf("updateSingleAtoms stop\n");
}
diff --git a/gromacs/pbc.c b/gromacs/pbc.c
index 31d8031..54c0bf8 100644
--- a/gromacs/pbc.c
+++ b/gromacs/pbc.c
@@ -22,6 +22,7 @@
*/
#include
#include
+#include
#include
#include
@@ -43,21 +44,48 @@ void initPbc(Atom* atom) {
/* update coordinates of ghost atoms */
/* uses mapping created in setupPbc */
-void updatePbc(Atom *atom, Parameter *param) {
+void updatePbc(Atom *atom, Parameter *param, int updateBoundingBoxes) {
int *border_map = atom->border_map;
int nlocal = atom->Nclusters_local;
MD_FLOAT xprd = param->xprd;
MD_FLOAT yprd = param->yprd;
MD_FLOAT zprd = param->zprd;
- for(int ci = 0; ci < atom->Nclusters_ghost; ci++) {
- MD_FLOAT *cptr = cluster_pos_ptr(nlocal + ci);
- MD_FLOAT *bmap_cptr = cluster_pos_ptr(border_map[ci]);
+ for(int cg = 0; cg < atom->Nclusters_ghost; cg++) {
+ const int ci = nlocal + cg;
+ MD_FLOAT *cptr = cluster_pos_ptr(ci);
+ MD_FLOAT *bmap_cptr = cluster_pos_ptr(border_map[cg]);
+ MD_FLOAT bbminx = INFINITY, bbmaxx = -INFINITY;
+ MD_FLOAT bbminy = INFINITY, bbmaxy = -INFINITY;
+ MD_FLOAT bbminz = INFINITY, bbmaxz = -INFINITY;
for(int cii = 0; cii < atom->clusters[ci].natoms; cii++) {
- cluster_x(cptr, cii) = cluster_x(bmap_cptr, cii) + atom->PBCx[ci] * xprd;
- cluster_y(cptr, cii) = cluster_y(bmap_cptr, cii) + atom->PBCy[ci] * yprd;
- cluster_z(cptr, cii) = cluster_z(bmap_cptr, cii) + atom->PBCz[ci] * zprd;
+ MD_FLOAT xtmp = cluster_x(bmap_cptr, cii) + atom->PBCx[cg] * xprd;
+ MD_FLOAT ytmp = cluster_y(bmap_cptr, cii) + atom->PBCy[cg] * yprd;
+ MD_FLOAT ztmp = cluster_z(bmap_cptr, cii) + atom->PBCz[cg] * zprd;
+
+ cluster_x(cptr, cii) = xtmp;
+ cluster_y(cptr, cii) = ytmp;
+ cluster_z(cptr, cii) = ztmp;
+
+ if(updateBoundingBoxes) {
+ // TODO: To create the bounding boxes faster, we can use SIMD operations
+ if(bbminx > xtmp) { bbminx = xtmp; }
+ if(bbmaxx < xtmp) { bbmaxx = xtmp; }
+ if(bbminy > ytmp) { bbminy = ytmp; }
+ if(bbmaxy < ytmp) { bbmaxy = ytmp; }
+ if(bbminz > ztmp) { bbminz = ztmp; }
+ if(bbmaxz < ztmp) { bbmaxz = ztmp; }
+ }
+ }
+
+ if(updateBoundingBoxes) {
+ atom->clusters[ci].bbminx = bbminx;
+ atom->clusters[ci].bbmaxx = bbmaxx;
+ atom->clusters[ci].bbminy = bbminy;
+ atom->clusters[ci].bbmaxy = bbmaxy;
+ atom->clusters[ci].bbminz = bbminz;
+ atom->clusters[ci].bbmaxz = bbmaxz;
}
}
}
@@ -181,5 +209,5 @@ void setupPbc(Atom *atom, Parameter *param) {
atom->Nclusters = atom->Nclusters_local + Nghost + 1;
// Update created ghost clusters positions
- updatePbc(atom, param);
+ updatePbc(atom, param, 1);
}
diff --git a/gromacs/vtk.c b/gromacs/vtk.c
index b3ff0e6..f339ccd 100644
--- a/gromacs/vtk.c
+++ b/gromacs/vtk.c
@@ -18,8 +18,11 @@ int write_atoms_to_vtk_file(const char* filename, Atom* atom, int timestep) {
fprintf(fp, "ASCII\n");
fprintf(fp, "DATASET UNSTRUCTURED_GRID\n");
fprintf(fp, "POINTS %d double\n", atom->Nlocal);
- for(int i = 0; i < atom->Nlocal; ++i) {
- fprintf(fp, "%.4f %.4f %.4f\n", atom_x(i), atom_y(i), atom_z(i));
+ for(int ci = 0; ci < atom->Nclusters_local; ++ci) {
+ MD_FLOAT *cptr = cluster_pos_ptr(ci);
+ for(int cii = 0; cii < atom->clusters[ci].natoms; ++cii) {
+ fprintf(fp, "%.4f %.4f %.4f\n", cluster_x(cptr, cii), cluster_y(cptr, cii), cluster_z(cptr, cii));
+ }
}
fprintf(fp, "\n\n");
fprintf(fp, "CELLS %d %d\n", atom->Nlocal, atom->Nlocal * 2);