diff --git a/src/compute_rdf.cpp b/src/compute_rdf.cpp index 7c5f433bb5..fe4be0cd3b 100644 --- a/src/compute_rdf.cpp +++ b/src/compute_rdf.cpp @@ -95,6 +95,7 @@ ComputeRDF::ComputeRDF(LAMMPS *lmp, int narg, char **arg) : typecount = new int[ntypes+1]; icount = new int[npairs]; jcount = new int[npairs]; + duplicates = new int[npairs]; } /* ---------------------------------------------------------------------- */ @@ -113,13 +114,14 @@ ComputeRDF::~ComputeRDF() delete [] typecount; delete [] icount; delete [] jcount; + delete [] duplicates; } /* ---------------------------------------------------------------------- */ void ComputeRDF::init() { - int i,m; + int i,j,m; if (force->pair) delr = force->pair->cutforce / nbin; else error->all(FLERR,"Compute rdf requires a pair style be defined"); @@ -143,12 +145,17 @@ void ComputeRDF::init() // 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]; @@ -156,10 +163,11 @@ void ComputeRDF::init() 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; // need an occasional half neighbor list - int irequest = neighbor->request(this,instance_me); neighbor->requests[irequest]->pair = 0; neighbor->requests[irequest]->compute = 1; @@ -265,10 +273,12 @@ void ComputeRDF::compute_array() MPI_Allreduce(hist[0],histall[0],npairs*nbin,MPI_DOUBLE,MPI_SUM,world); // convert counts to g(r) and coord(r) and copy into output array - // nideal = # of J atoms surrounding single I atom in a single bin - // assuming J atoms are at uniform density + // vfrac = fraction of volume in shell m + // npairs = number of pairs, corrected for duplicates + // duplicates = pairs in which both atoms are the same - double constant,nideal,gr,ncoord,rlower,rupper; + double constant,vfrac,gr,ncoord,rlower,rupper; + int pairtot; if (domain->dimension == 3) { constant = 4.0*MY_PI / (3.0*domain->xprd*domain->yprd*domain->zprd); @@ -278,12 +288,13 @@ void ComputeRDF::compute_array() for (ibin = 0; ibin < nbin; ibin++) { rlower = ibin*delr; rupper = (ibin+1)*delr; - nideal = constant * - (rupper*rupper*rupper - rlower*rlower*rlower) * jcount[m]; - if (icount[m]*nideal != 0.0) - gr = histall[m][ibin] / (icount[m]*nideal); + vfrac = constant * (rupper*rupper*rupper - rlower*rlower*rlower); + pairtot = icount[m] * jcount[m] - duplicates[m]; + if (vfrac * pairtot != 0.0) + gr = histall[m][ibin] / (vfrac * pairtot); else gr = 0.0; - ncoord += gr*nideal; + if (icount[m] != 0) + ncoord += gr * vfrac * pairtot / icount[m]; array[ibin][1+2*m] = gr; array[ibin][2+2*m] = ncoord; } @@ -297,11 +308,13 @@ void ComputeRDF::compute_array() for (ibin = 0; ibin < nbin; ibin++) { rlower = ibin*delr; rupper = (ibin+1)*delr; - nideal = constant * (rupper*rupper - rlower*rlower) * jcount[m]; - if (icount[m]*nideal != 0.0) - gr = histall[m][ibin] / (icount[m]*nideal); + vfrac = constant * (rupper*rupper - rlower*rlower); + pairtot = icount[m] * jcount[m] - duplicates[m]; + if (vfrac * pairtot != 0.0) + gr = histall[m][ibin] / (vfrac * pairtot); else gr = 0.0; - ncoord += gr*nideal; + if (icount[m] != 0) + ncoord += gr * vfrac * pairtot / icount[m]; array[ibin][1+2*m] = gr; array[ibin][2+2*m] = ncoord; } diff --git a/src/compute_rdf.h b/src/compute_rdf.h index 794fcc9ee6..0da56609e4 100644 --- a/src/compute_rdf.h +++ b/src/compute_rdf.h @@ -44,7 +44,7 @@ class ComputeRDF : public Compute { double **histall; // summed histogram bins across all procs int *typecount; - int *icount,*jcount; + int *icount,*jcount,*duplicates; class NeighList *list; // half neighbor list }; diff --git a/src/finish.cpp b/src/finish.cpp index 36a7973724..b2ca44945a 100644 --- a/src/finish.cpp +++ b/src/finish.cpp @@ -121,9 +121,9 @@ void Finish::end(int flag) const char fmt1[] = "Loop time of %g on %d procs " "for %d steps with " BIGINT_FORMAT " atoms\n\n"; if (screen) fprintf(screen,fmt1,time_loop,ntasks,update->nsteps, - atom->natoms,cpu_loop); + atom->natoms); if (logfile) fprintf(logfile,fmt1,time_loop,ntasks,update->nsteps, - atom->natoms,cpu_loop); + atom->natoms); // Gromacs/NAMD-style performance metric for suitable unit settings @@ -749,8 +749,10 @@ void Finish::end(int flag) nspec_all/atom->natoms); fprintf(screen,"Neighbor list builds = " BIGINT_FORMAT "\n", neighbor->ncalls); - fprintf(screen,"Dangerous builds = " BIGINT_FORMAT "\n", - neighbor->ndanger); + if (neighbor->dist_check) + fprintf(screen,"Dangerous builds = " BIGINT_FORMAT "\n", + neighbor->ndanger); + else fprintf(screen,"Dangerous builds not checked\n"); } if (logfile) { if (nall < 2.0e9) @@ -764,8 +766,10 @@ void Finish::end(int flag) nspec_all/atom->natoms); fprintf(logfile,"Neighbor list builds = " BIGINT_FORMAT "\n", neighbor->ncalls); - fprintf(logfile,"Dangerous builds = " BIGINT_FORMAT "\n", - neighbor->ndanger); + if (neighbor->dist_check) + fprintf(logfile,"Dangerous builds = " BIGINT_FORMAT "\n", + neighbor->ndanger); + else fprintf(logfile,"Dangerous builds not checked\n"); } } } diff --git a/src/info.cpp b/src/info.cpp index c5a663ee6b..f36cee7b83 100644 --- a/src/info.cpp +++ b/src/info.cpp @@ -161,7 +161,8 @@ void Info::command(int narg, char **arg) if (flags & CONFIG) { - fprintf(out,"\nLAMMPS version: %s\n", universe->version); + fprintf(out,"\nLAMMPS version: %s / %s\n", + universe->version, universe->num_ver); fprintf(out,"sizeof(smallint): %3d-bit\n",(int)sizeof(smallint)*8); fprintf(out,"sizeof(imageint): %3d-bit\n",(int)sizeof(imageint)*8); fprintf(out,"sizeof(tagint): %3d-bit\n",(int)sizeof(tagint)*8); diff --git a/src/main.cpp b/src/main.cpp index 15c736119a..d85a79d831 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -26,16 +26,10 @@ int main(int argc, char **argv) { MPI_Init(&argc,&argv); - //try { - LAMMPS *lammps = new LAMMPS(argc,argv,MPI_COMM_WORLD); - lammps->input->file(); - delete lammps; - //} - //catch(...) { - //int me; - //MPI_Comm_rank(MPI_COMM_WORLD,&me); - //fprintf(stderr,"Unknown exception caught on rank %d\n",me); - //} + LAMMPS *lammps = new LAMMPS(argc,argv,MPI_COMM_WORLD); + + lammps->input->file(); + delete lammps; MPI_Barrier(MPI_COMM_WORLD); MPI_Finalize(); diff --git a/src/neighbor.cpp b/src/neighbor.cpp index 1f959baccd..4e93ae9cd9 100644 --- a/src/neighbor.cpp +++ b/src/neighbor.cpp @@ -811,6 +811,8 @@ void Neighbor::init() fprintf(logfile," %d neighbor list requests\n",nrequest); fprintf(logfile," update every %d steps, delay %d steps, check %s\n", every,delay,dist_check ? "yes" : "no"); + fprintf(logfile," max neighbors/atom: %d, page size: %d\n", + oneatom, pgsize); fprintf(logfile," master list distance cutoff = %g\n",cutneighmax); fprintf(logfile," ghost atom cutoff = %g\n",cutghost); if (style != NSQ) @@ -823,6 +825,8 @@ void Neighbor::init() fprintf(screen," %d neighbor list requests\n",nrequest); fprintf(screen," update every %d steps, delay %d steps, check %s\n", every,delay,dist_check ? "yes" : "no"); + fprintf(screen," max neighbors/atom: %d, page size: %d\n", + oneatom, pgsize); fprintf(screen," master list distance cutoff = %g\n",cutneighmax); fprintf(screen," ghost atom cutoff = %g\n",cutghost); if (style != NSQ)