From fe888e46221910931096e959d48943114ef974ee Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 7 Jul 2017 15:39:25 -0400 Subject: [PATCH] add support for recomputing normalization factors and finite size correction during --- doc/src/compute_rdf.txt | 11 ++++- src/compute_rdf.cpp | 96 ++++++++++++++++++++++++++--------------- src/compute_rdf.h | 2 + 3 files changed, 73 insertions(+), 36 deletions(-) diff --git a/doc/src/compute_rdf.txt b/doc/src/compute_rdf.txt index acbc0e4f0c..e462e85fc0 100644 --- a/doc/src/compute_rdf.txt +++ b/doc/src/compute_rdf.txt @@ -180,9 +180,18 @@ will register an arbitrarily large spike at whatever distance they happen to be at, and zero everywhere else. Coord(r) will show a step change from zero to one at the location of the spike in g(r). +NOTE: compute rdf can handle dynamic groups and systems where atoms +are added or removed, but this causes that certain normalization +parameters need to be recomputed in every step and include collective +communication operations. This will reduce performance and limit +parallel efficiency and scaling. For systems, where only the type +of atoms changes (e.g. when using "fix atom/swap"_fix_atom_swap.html), +you need to explicitly request the dynamic normalization updates +via "compute_modify dynamic yes"_compute_modify.html + [Related commands:] -"fix ave/time"_fix_ave_time.html +"fix ave/time"_fix_ave_time.html, "compute_modify"_compute_modify.html [Default:] diff --git a/src/compute_rdf.cpp b/src/compute_rdf.cpp index fc6ad6d8b6..72eb46f6f5 100644 --- a/src/compute_rdf.cpp +++ b/src/compute_rdf.cpp @@ -48,7 +48,6 @@ ComputeRDF::ComputeRDF(LAMMPS *lmp, int narg, char **arg) : array_flag = 1; extarray = 0; - dynamic_group_allow = 0; nbin = force->inumeric(FLERR,arg[3]); if (nbin < 1) error->all(FLERR,"Illegal compute rdf command"); @@ -125,6 +124,9 @@ ComputeRDF::ComputeRDF(LAMMPS *lmp, int narg, char **arg) : icount = new int[npairs]; jcount = new int[npairs]; duplicates = new int[npairs]; + + dynamic = 0; + natoms_old = 0; } /* ---------------------------------------------------------------------- */ @@ -150,7 +152,6 @@ ComputeRDF::~ComputeRDF() void ComputeRDF::init() { - int i,j,m; if (!force->pair && !cutflag) error->all(FLERR,"Compute rdf requires a pair style be defined " @@ -184,40 +185,12 @@ void ComputeRDF::init() for (int i = 0; i < nbin; i++) array[i][0] = (i+0.5) * delr; - // count atoms of each type that are also in group + // initialize normalization, finite size correction, and changing atom counts - int *mask = atom->mask; - int *type = atom->type; - int nlocal = atom->nlocal; - int ntypes = atom->ntypes; - - for (i = 1; i <= ntypes; i++) typecount[i] = 0; - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) typecount[type[i]]++; - - // icount = # of I atoms participating in I,J pairs for each histogram - // jcount = # of J atoms participating in I,J pairs for each histogram - // duplicates = # of atoms in both groups I and J for each histogram - - for (m = 0; m < npairs; m++) { - icount[m] = 0; - for (i = ilo[m]; i <= ihi[m]; i++) icount[m] += typecount[i]; - jcount[m] = 0; - for (i = jlo[m]; i <= jhi[m]; i++) jcount[m] += typecount[i]; - duplicates[m] = 0; - for (i = ilo[m]; i <= ihi[m]; i++) - for (j = jlo[m]; j <= jhi[m]; j++) - if (i == j) duplicates[m] += typecount[i]; - } - - int *scratch = new int[npairs]; - MPI_Allreduce(icount,scratch,npairs,MPI_INT,MPI_SUM,world); - for (i = 0; i < npairs; i++) icount[i] = scratch[i]; - MPI_Allreduce(jcount,scratch,npairs,MPI_INT,MPI_SUM,world); - for (i = 0; i < npairs; i++) jcount[i] = scratch[i]; - MPI_Allreduce(duplicates,scratch,npairs,MPI_INT,MPI_SUM,world); - for (i = 0; i < npairs; i++) duplicates[i] = scratch[i]; - delete [] scratch; + natoms_old = atom->natoms; + dynamic = group->dynamic[igroup]; + if (dynamic_user) dynamic = 1; + init_norm(); // need an occasional half neighbor list // if user specified, request a cutoff = cutoff_user + skin @@ -246,6 +219,48 @@ void ComputeRDF::init_list(int id, NeighList *ptr) /* ---------------------------------------------------------------------- */ +void ComputeRDF::init_norm() +{ + int i,j,m; + + // count atoms of each type that are also in group + + const int nlocal = atom->nlocal; + const int ntypes = atom->ntypes; + const int * const mask = atom->mask; + const int * const type = atom->type; + + for (i = 1; i <= ntypes; i++) typecount[i] = 0; + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) typecount[type[i]]++; + + // icount = # of I atoms participating in I,J pairs for each histogram + // jcount = # of J atoms participating in I,J pairs for each histogram + // duplicates = # of atoms in both groups I and J for each histogram + + for (m = 0; m < npairs; m++) { + icount[m] = 0; + for (i = ilo[m]; i <= ihi[m]; i++) icount[m] += typecount[i]; + jcount[m] = 0; + for (i = jlo[m]; i <= jhi[m]; i++) jcount[m] += typecount[i]; + duplicates[m] = 0; + for (i = ilo[m]; i <= ihi[m]; i++) + for (j = jlo[m]; j <= jhi[m]; j++) + if (i == j) duplicates[m] += typecount[i]; + } + + int *scratch = new int[npairs]; + MPI_Allreduce(icount,scratch,npairs,MPI_INT,MPI_SUM,world); + for (i = 0; i < npairs; i++) icount[i] = scratch[i]; + MPI_Allreduce(jcount,scratch,npairs,MPI_INT,MPI_SUM,world); + for (i = 0; i < npairs; i++) jcount[i] = scratch[i]; + MPI_Allreduce(duplicates,scratch,npairs,MPI_INT,MPI_SUM,world); + for (i = 0; i < npairs; i++) duplicates[i] = scratch[i]; + delete [] scratch; +} + +/* ---------------------------------------------------------------------- */ + void ComputeRDF::compute_array() { int i,j,m,ii,jj,inum,jnum,itype,jtype,ipair,jpair,ibin,ihisto; @@ -253,6 +268,17 @@ void ComputeRDF::compute_array() int *ilist,*jlist,*numneigh,**firstneigh; double factor_lj,factor_coul; + if (natoms_old != atom->natoms) { + dynamic = 1; + natoms_old = atom->natoms; + } + + // if the number of atoms has changed or we have a dynamic group + // or dynamic updates are requested (e.g. when changing atom types) + // we need to recompute some normalization parameters + + if (dynamic) init_norm(); + invoked_array = update->ntimestep; // invoke half neighbor list (will copy or build if necessary) diff --git a/src/compute_rdf.h b/src/compute_rdf.h index 61d99e0881..4188799b74 100644 --- a/src/compute_rdf.h +++ b/src/compute_rdf.h @@ -51,6 +51,8 @@ class ComputeRDF : public Compute { int *duplicates; class NeighList *list; // half neighbor list + void init_norm(); + bigint natoms_old; }; }