From e64c3345bc0cd5900b8ad4148c89d825f6f36867 Mon Sep 17 00:00:00 2001 From: Rafael Ravedutti Date: Mon, 31 Jan 2022 23:46:20 +0100 Subject: [PATCH] Fix a few more bugs on gromacs variant Signed-off-by: Rafael Ravedutti --- gromacs/force_lj.c | 11 +++++++++- gromacs/includes/atom.h | 2 +- gromacs/includes/vtk.h | 3 ++- gromacs/main.c | 6 ++++-- gromacs/neighbor.c | 19 +++++++++-------- gromacs/pbc.c | 25 +++++++++++----------- gromacs/vtk.c | 47 +++++++++++++++++++++++++++++++++++++++-- 7 files changed, 85 insertions(+), 28 deletions(-) diff --git a/gromacs/force_lj.c b/gromacs/force_lj.c index 46cdd94..7ff61cb 100644 --- a/gromacs/force_lj.c +++ b/gromacs/force_lj.c @@ -67,7 +67,7 @@ double computeForceLJ(Parameter *param, Atom *atom, Neighbor *neighbor, Stats *s MD_FLOAT fiy = 0; MD_FLOAT fiz = 0; - for(int cjj = 0; cjj < CLUSTER_DIM_N; cjj++) { + for(int cjj = 0; cjj < atom->clusters[cj].natoms; cjj++) { if(ci != cj || cii != cjj) { MD_FLOAT delx = xtmp - cluster_x(cjptr, cjj); MD_FLOAT dely = ytmp - cluster_y(cjptr, cjj); @@ -80,6 +80,11 @@ double computeForceLJ(Parameter *param, Atom *atom, Neighbor *neighbor, Stats *s fix += delx * force; fiy += dely * force; fiz += delz * force; + /* + if(force < -50.0 || force > 50.0) { + fprintf(stdout, "%d-%d/%d-%d: %f, %f, %f ---- %f, %f, %f, %f\n", ci, cii, cj, cjj, xtmp, ytmp, ztmp, fix, fiy, fiz, force); + } + */ } } } @@ -87,6 +92,10 @@ double computeForceLJ(Parameter *param, Atom *atom, Neighbor *neighbor, Stats *s cluster_x(cifptr, cii) += fix; cluster_y(cifptr, cii) += fiy; cluster_z(cifptr, cii) += fiz; + /* + if(fix < -100.0 || fix > 100.0 || fiy < -100.0 || fiy > 100.0 || fiz < -100.0 || fiz > 100.0) { + fprintf(stdout, "%d-%d: %f, %f, %f ---- %f, %f, %f\n", ci, cii, xtmp, ytmp, ztmp, fix, fiy, fiz); + }*/ } } diff --git a/gromacs/includes/atom.h b/gromacs/includes/atom.h index 17065cf..4b04fda 100644 --- a/gromacs/includes/atom.h +++ b/gromacs/includes/atom.h @@ -25,8 +25,8 @@ #ifndef __ATOM_H_ #define __ATOM_H_ -#define CLUSTER_DIM_N 4 #define CLUSTER_DIM_M 4 +#define CLUSTER_DIM_N 4 typedef struct { int bin; diff --git a/gromacs/includes/vtk.h b/gromacs/includes/vtk.h index 0f03a74..66b3b32 100644 --- a/gromacs/includes/vtk.h +++ b/gromacs/includes/vtk.h @@ -24,5 +24,6 @@ #ifndef __VTK_H_ #define __VTK_H_ -extern int write_atoms_to_vtk_file(const char* filename, Atom* atom, int timestep); +extern int write_local_atoms_to_vtk_file(const char* filename, Atom* atom, int timestep); +extern int write_ghost_atoms_to_vtk_file(const char* filename, Atom* atom, int timestep); #endif diff --git a/gromacs/main.c b/gromacs/main.c index b4c1bd9..3556c6a 100644 --- a/gromacs/main.c +++ b/gromacs/main.c @@ -260,7 +260,8 @@ int main(int argc, char** argv) { timer[TOTAL] = getTimeStamp(); if(param.vtk_file != NULL) { - write_atoms_to_vtk_file(param.vtk_file, &atom, 0); + write_local_atoms_to_vtk_file(param.vtk_file, &atom, 0); + write_ghost_atoms_to_vtk_file(param.vtk_file, &atom, 0); } for(int n = 0; n < param.ntimes; n++) { @@ -289,7 +290,8 @@ int main(int argc, char** argv) { } if(param.vtk_file != NULL) { - write_atoms_to_vtk_file(param.vtk_file, &atom, n + 1); + write_local_atoms_to_vtk_file(param.vtk_file, &atom, n + 1); + write_ghost_atoms_to_vtk_file(param.vtk_file, &atom, n + 1); } } diff --git a/gromacs/neighbor.c b/gromacs/neighbor.c index 1e65b33..ad7455a 100644 --- a/gromacs/neighbor.c +++ b/gromacs/neighbor.c @@ -289,7 +289,7 @@ void buildNeighbor(Atom *atom, Neighbor *neighbor) { /* printf("\ncutneighsq = %f, rbb_sq = %f\n", cutneighsq, rbb_sq); - for(int ci = 4; ci < 6; ci++) { + for(int ci = 0; 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); @@ -377,7 +377,6 @@ 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; while(resize > 0) { @@ -387,7 +386,7 @@ void binAtoms(Atom *atom) { bincount[i] = 0; } - for(int i = 0; i < nall; i++) { + for(int i = 0; i < atom->Nlocal; 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) { @@ -507,8 +506,15 @@ void buildClusters(Atom *atom) { } } + printf("buildClusters end\n"); +} + +void binClusters(Atom *atom) { + printf("binClusters start\n"); + /* - for(int ci = 4; ci < 9; ci++) { + fprintf(stdout, "Nghost = %d\n", atom->Nclusters_ghost); + for(int ci = atom->Nclusters_local; ci < atom->Nclusters_local + 4; 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); @@ -518,11 +524,6 @@ void buildClusters(Atom *atom) { } */ - printf("buildClusters end\n"); -} - -void binClusters(Atom *atom) { - printf("binClusters start\n"); const int nlocal = atom->Nclusters_local; int resize = 1; diff --git a/gromacs/pbc.c b/gromacs/pbc.c index 54c0bf8..3d93c3b 100644 --- a/gromacs/pbc.c +++ b/gromacs/pbc.c @@ -122,19 +122,18 @@ void updateAtomsPbc(Atom *atom, Parameter *param) { * defining ghost atoms around domain * only creates mapping and coordinate corrections * that are then enforced in updatePbc */ -#define ADDGHOST(dx,dy,dz) \ - Nghost++; \ - border_map[Nghost] = ci; \ - atom->PBCx[Nghost] = dx; \ - atom->PBCy[Nghost] = dy; \ - atom->PBCz[Nghost] = dz; \ - copy_cluster_types(atom, atom->Nclusters_local + Nghost, ci) - -void copy_cluster_types(Atom *atom, int dest, int src) { - for(int cii = 0; cii < atom->clusters[src].natoms; cii++) { - atom->clusters[dest].type[cii] = atom->clusters[src].type[cii]; +#define ADDGHOST(dx,dy,dz); \ + Nghost++; \ + const int g_atom_idx = atom->Nclusters_local + Nghost; \ + border_map[Nghost] = ci; \ + atom->PBCx[Nghost] = dx; \ + atom->PBCy[Nghost] = dy; \ + atom->PBCz[Nghost] = dz; \ + atom->clusters[g_atom_idx].natoms = atom->clusters[ci].natoms; \ + Nghost_atoms += atom->clusters[g_atom_idx].natoms; \ + for(int cii = 0; cii < atom->clusters[ci].natoms; cii++) { \ + atom->clusters[g_atom_idx].type[cii] = atom->clusters[ci].type[cii]; \ } -} /* internal subroutines */ void growPbc(Atom* atom) { @@ -154,6 +153,7 @@ void setupPbc(Atom *atom, Parameter *param) { MD_FLOAT zprd = param->zprd; MD_FLOAT Cutneigh = param->cutneigh; int Nghost = -1; + int Nghost_atoms = 0; for(int ci = 0; ci < atom->Nclusters_local; ci++) { if (atom->Nclusters_local + Nghost + 7 >= atom->Nclusters_max) { @@ -205,6 +205,7 @@ void setupPbc(Atom *atom, Parameter *param) { } // increase by one to make it the ghost atom count + atom->Nghost = Nghost_atoms; atom->Nclusters_ghost = Nghost + 1; atom->Nclusters = atom->Nclusters_local + Nghost + 1; diff --git a/gromacs/vtk.c b/gromacs/vtk.c index f339ccd..4b788ea 100644 --- a/gromacs/vtk.c +++ b/gromacs/vtk.c @@ -3,9 +3,9 @@ #include -int write_atoms_to_vtk_file(const char* filename, Atom* atom, int timestep) { +int write_local_atoms_to_vtk_file(const char* filename, Atom* atom, int timestep) { char timestep_filename[128]; - snprintf(timestep_filename, sizeof timestep_filename, "%s_%d.vtk", filename, timestep); + snprintf(timestep_filename, sizeof timestep_filename, "%s_local_%d.vtk", filename, timestep); FILE* fp = fopen(timestep_filename, "wb"); if(fp == NULL) { @@ -45,3 +45,46 @@ int write_atoms_to_vtk_file(const char* filename, Atom* atom, int timestep) { fclose(fp); return 0; } + +int write_ghost_atoms_to_vtk_file(const char* filename, Atom* atom, int timestep) { + char timestep_filename[128]; + snprintf(timestep_filename, sizeof timestep_filename, "%s_ghost_%d.vtk", filename, timestep); + FILE* fp = fopen(timestep_filename, "wb"); + + if(fp == NULL) { + fprintf(stderr, "Could not open VTK file for writing!\n"); + return -1; + } + + fprintf(fp, "# vtk DataFile Version 2.0\n"); + fprintf(fp, "Particle data\n"); + fprintf(fp, "ASCII\n"); + fprintf(fp, "DATASET UNSTRUCTURED_GRID\n"); + fprintf(fp, "POINTS %d double\n", atom->Nghost); + for(int ci = atom->Nclusters_local; ci < atom->Nclusters_local + atom->Nclusters_ghost; ++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->Nghost, atom->Nghost * 2); + for(int i = 0; i < atom->Nghost; ++i) { + fprintf(fp, "1 %d\n", i); + } + fprintf(fp, "\n\n"); + fprintf(fp, "CELL_TYPES %d\n", atom->Nghost); + for(int i = 0; i < atom->Nghost; ++i) { + fprintf(fp, "1\n"); + } + fprintf(fp, "\n\n"); + fprintf(fp, "POINT_DATA %d\n", atom->Nghost); + fprintf(fp, "SCALARS mass double\n"); + fprintf(fp, "LOOKUP_TABLE default\n"); + for(int i = 0; i < atom->Nghost; i++) { + fprintf(fp, "1.0\n"); + } + fprintf(fp, "\n\n"); + fclose(fp); + return 0; +}