diff --git a/src/GPU/pair_coul_long_gpu.cpp b/src/GPU/pair_coul_long_gpu.cpp index a1864a4002..465d82232f 100644 --- a/src/GPU/pair_coul_long_gpu.cpp +++ b/src/GPU/pair_coul_long_gpu.cpp @@ -144,7 +144,14 @@ void PairCoulLongGPU::init_style() if (force->newton_pair) error->all(FLERR,"Cannot use newton pair with coul/long/gpu pair style"); - // Repeat cutsq calculation because done after call to init_style + // Call init_one calculation make sure scale is correct + for (int i = 1; i <= atom->ntypes; i++) { + for (int j = i; j <= atom->ntypes; j++) { + if (setflag[i][j] != 0 || (setflag[i][i] != 0 && setflag[j][j] != 0)) { + double cut = init_one(i,j); + } + } + } double cell_size = cut_coul + neighbor->skin; cut_coulsq = cut_coul * cut_coul; diff --git a/src/MOLECULE/bond_table.cpp b/src/MOLECULE/bond_table.cpp index fc42faeeea..3a9f7cc4d0 100644 --- a/src/MOLECULE/bond_table.cpp +++ b/src/MOLECULE/bond_table.cpp @@ -29,9 +29,10 @@ using namespace LAMMPS_NS; -enum{LINEAR,SPLINE}; +enum{NONE,LINEAR,SPLINE}; #define MAXLINE 1024 +#define BIGNUM 1.0e300 /* ---------------------------------------------------------------------- */ @@ -134,6 +135,7 @@ void BondTable::settings(int narg, char **arg) { if (narg != 2) error->all(FLERR,"Illegal bond_style command"); + tabstyle = NONE; if (strcmp(arg[0],"linear") == 0) tabstyle = LINEAR; else if (strcmp(arg[0],"spline") == 0) tabstyle = SPLINE; else error->all(FLERR,"Unknown table style in bond style table"); @@ -289,6 +291,7 @@ void BondTable::free_table(Table *tb) void BondTable::read_table(Table *tb, char *file, char *keyword) { char line[MAXLINE]; + double emin = BIGNUM; // open file @@ -326,14 +329,65 @@ void BondTable::read_table(Table *tb, char *file, char *keyword) // read r,e,f table values from file int itmp; + int cerror = 0; + int r0idx = -1; + fgets(line,MAXLINE,fp); for (int i = 0; i < tb->ninput; i++) { - fgets(line,MAXLINE,fp); - sscanf(line,"%d %lg %lg %lg", - &itmp,&tb->rfile[i],&tb->efile[i],&tb->ffile[i]); + if (NULL == fgets(line,MAXLINE,fp)) + error->one(FLERR,"Premature end of file in bond table"); + if (4 != sscanf(line,"%d %lg %lg %lg", + &itmp,&tb->rfile[i],&tb->efile[i],&tb->ffile[i])) ++cerror; + if (tb->efile[i] < emin) { + emin = tb->efile[i]; + r0idx = i; + } + } + fclose(fp); + + // infer r0 from minimum of potential, if not given explicitly + + if ((tb->r0 == 0.0) && (r0idx >= 0)) tb->r0 = tb->rfile[r0idx]; + + // warn if force != dE/dr at any point that is not an inflection point + // check via secant approximation to dE/dr + // skip two end points since do not have surrounding secants + // inflection point is where curvature changes sign + + double r,e,f,rprev,rnext,eprev,enext,fleft,fright; + + int ferror = 0; + for (int i = 1; i < tb->ninput-1; i++) { + r = tb->rfile[i]; + rprev = tb->rfile[i-1]; + rnext = tb->rfile[i+1]; + e = tb->efile[i]; + eprev = tb->efile[i-1]; + enext = tb->efile[i+1]; + f = tb->ffile[i]; + fleft = - (e-eprev) / (r-rprev); + fright = - (enext-e) / (rnext-r); + if (f < fleft && f < fright) ferror++; + if (f > fleft && f > fright) ferror++; + //printf("Values %d: %g %g %g\n",i,r,e,f); + //printf(" secant %d %d %g: %g %g %g\n",i,ferror,r,fleft,fright,f); } - fclose(fp); + if (ferror) { + char str[128]; + sprintf(str,"%d of %d force values in table are inconsistent with -dE/dr.\n" + " Should only be flagged at inflection points",ferror,tb->ninput); + error->warning(FLERR,str); + } + + // warn if data was read incompletely, e.g. columns were missing + + if (cerror) { + char str[128]; + sprintf(str,"%d of %d lines in table were incomplete or could not be" + " parsed completely",cerror,tb->ninput); + error->warning(FLERR,str); + } } /* ---------------------------------------------------------------------- diff --git a/src/USER-FEP/README b/src/USER-FEP/README index 967f8e6efe..8c942a90d8 100644 --- a/src/USER-FEP/README +++ b/src/USER-FEP/README @@ -11,3 +11,8 @@ There are auxiliary tools for using this package in tools/fep. The person who created this package is Agilio Padua at Université Blaise Pascal Clermont-Ferrand (agilio.padua at univ-bpclermont.fr) Contact him directly if you have questions. + +Pair style morse/soft was contributed by Stefan Paquay, s.paquay at tue.nl +Applied Physics/Theory of Polymers and Soft Matter, +Eindhoven University of Technology (TU/e), The Netherlands +Contact him in case of problems with this pair style. diff --git a/src/compute_dipole_chunk.cpp b/src/compute_dipole_chunk.cpp new file mode 100644 index 0000000000..449fb002b6 --- /dev/null +++ b/src/compute_dipole_chunk.cpp @@ -0,0 +1,295 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "string.h" +#include "compute_dipole_chunk.h" +#include "atom.h" +#include "update.h" +#include "modify.h" +#include "compute_chunk_atom.h" +#include "domain.h" +#include "memory.h" +#include "error.h" +#include "math_special.h" + +using namespace LAMMPS_NS; +using namespace MathSpecial; + +enum { MASSCENTER, GEOMCENTER }; + +/* ---------------------------------------------------------------------- */ + +ComputeDipoleChunk::ComputeDipoleChunk(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg) +{ + if ((narg != 4) && (narg != 5)) error->all(FLERR,"Illegal compute dipole/chunk command"); + + array_flag = 1; + size_array_cols = 4; + size_array_rows = 0; + size_array_rows_variable = 1; + extarray = 0; + + // ID of compute chunk/atom + + int n = strlen(arg[3]) + 1; + idchunk = new char[n]; + strcpy(idchunk,arg[3]); + + usecenter = MASSCENTER; + + if (narg == 5) { + if (strncmp(arg[4],"geom",4) == 0) usecenter = GEOMCENTER; + else if (strcmp(arg[4],"mass") == 0) usecenter = MASSCENTER; + else error->all(FLERR,"Illegal compute dipole/chunk command"); + } + + init(); + + // chunk-based data + + nchunk = 1; + maxchunk = 0; + massproc = masstotal = NULL; + chrgproc = chrgtotal = NULL; + com = comall = NULL; + dipole = dipoleall = NULL; + allocate(); +} + +/* ---------------------------------------------------------------------- */ + +ComputeDipoleChunk::~ComputeDipoleChunk() +{ + delete [] idchunk; + memory->destroy(massproc); + memory->destroy(masstotal); + memory->destroy(chrgproc); + memory->destroy(chrgtotal); + memory->destroy(com); + memory->destroy(comall); + memory->destroy(dipole); + memory->destroy(dipoleall); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeDipoleChunk::init() +{ + int icompute = modify->find_compute(idchunk); + if (icompute < 0) + error->all(FLERR,"Chunk/atom compute does not exist for " + "compute dipole/chunk"); + cchunk = (ComputeChunkAtom *) modify->compute[icompute]; + if (strcmp(cchunk->style,"chunk/atom") != 0) + error->all(FLERR,"Compute dipole/chunk does not use chunk/atom compute"); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeDipoleChunk::compute_array() +{ + int i,index; + double massone; + double unwrap[3]; + + invoked_array = update->ntimestep; + + // compute chunk/atom assigns atoms to chunk IDs + // extract ichunk index vector from compute + // ichunk = 1 to Nchunk for included atoms, 0 for excluded atoms + + nchunk = cchunk->setup_chunks(); + cchunk->compute_ichunk(); + int *ichunk = cchunk->ichunk; + + if (nchunk > maxchunk) allocate(); + size_array_rows = nchunk; + + // zero local per-chunk values + + for (int i = 0; i < nchunk; i++) { + massproc[i] = chrgproc[i] = 0.0; + com[i][0] = com[i][1] = com[i][2] = 0.0; + dipole[i][0] = dipole[i][1] = dipole[i][2] = dipole[i][3] = 0.0; + } + + // compute COM for each chunk + + double **x = atom->x; + int *mask = atom->mask; + int *type = atom->type; + imageint *image = atom->image; + double *mass = atom->mass; + double *rmass = atom->rmass; + double *q = atom->q; + double **mu = atom->mu; + + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + index = ichunk[i]-1; + if (index < 0) continue; + if (usecenter == MASSCENTER) { + if (rmass) massone = rmass[i]; + else massone = mass[type[i]]; + } else massone = 1.0; // usecenter == GEOMCENTER + + domain->unmap(x[i],image[i],unwrap); + massproc[index] += massone; + if (atom->q_flag) chrgproc[index] += atom->q[i]; + com[index][0] += unwrap[0] * massone; + com[index][1] += unwrap[1] * massone; + com[index][2] += unwrap[2] * massone; + } + + MPI_Allreduce(massproc,masstotal,nchunk,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(chrgproc,chrgtotal,nchunk,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&com[0][0],&comall[0][0],3*nchunk,MPI_DOUBLE,MPI_SUM,world); + + for (int i = 0; i < nchunk; i++) { + if (masstotal[i] == 0.0) masstotal[i] = 1.0; + comall[i][0] /= masstotal[i]; + comall[i][1] /= masstotal[i]; + comall[i][2] /= masstotal[i]; + } + + // compute dipole for each chunk + + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + index = ichunk[i]-1; + if (index < 0) continue; + domain->unmap(x[i],image[i],unwrap); + if (atom->q_flag) { + dipole[index][0] += q[i]*unwrap[0]; + dipole[index][1] += q[i]*unwrap[1]; + dipole[index][2] += q[i]*unwrap[2]; + } + if (atom->mu_flag) { + dipole[index][0] += mu[i][0]; + dipole[index][1] += mu[i][1]; + dipole[index][2] += mu[i][2]; + } + } + } + + MPI_Allreduce(&dipole[0][0],&dipoleall[0][0],4*nchunk, + MPI_DOUBLE,MPI_SUM,world); + + for (i = 0; i < nchunk; i++) { + // correct for position dependence with charged chunks + dipoleall[i][0] -= chrgtotal[i]*comall[i][0]; + dipoleall[i][1] -= chrgtotal[i]*comall[i][1]; + dipoleall[i][2] -= chrgtotal[i]*comall[i][2]; + // compute total dipole moment + dipoleall[i][3] = sqrt(square(dipoleall[i][0]) + + square(dipoleall[i][1]) + + square(dipoleall[i][2])); + } +} + +/* ---------------------------------------------------------------------- + lock methods: called by fix ave/time + these methods insure vector/array size is locked for Nfreq epoch + by passing lock info along to compute chunk/atom +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + increment lock counter +------------------------------------------------------------------------- */ + +void ComputeDipoleChunk::lock_enable() +{ + cchunk->lockcount++; +} + +/* ---------------------------------------------------------------------- + decrement lock counter in compute chunk/atom, it if still exists +------------------------------------------------------------------------- */ + +void ComputeDipoleChunk::lock_disable() +{ + int icompute = modify->find_compute(idchunk); + if (icompute >= 0) { + cchunk = (ComputeChunkAtom *) modify->compute[icompute]; + cchunk->lockcount--; + } +} + +/* ---------------------------------------------------------------------- + calculate and return # of chunks = length of vector/array +------------------------------------------------------------------------- */ + +int ComputeDipoleChunk::lock_length() +{ + nchunk = cchunk->setup_chunks(); + return nchunk; +} + +/* ---------------------------------------------------------------------- + set the lock from startstep to stopstep +------------------------------------------------------------------------- */ + +void ComputeDipoleChunk::lock(Fix *fixptr, bigint startstep, bigint stopstep) +{ + cchunk->lock(fixptr,startstep,stopstep); +} + +/* ---------------------------------------------------------------------- + unset the lock +------------------------------------------------------------------------- */ + +void ComputeDipoleChunk::unlock(Fix *fixptr) +{ + cchunk->unlock(fixptr); +} + +/* ---------------------------------------------------------------------- + free and reallocate per-chunk arrays +------------------------------------------------------------------------- */ + +void ComputeDipoleChunk::allocate() +{ + memory->destroy(massproc); + memory->destroy(masstotal); + memory->destroy(chrgproc); + memory->destroy(chrgtotal); + memory->destroy(com); + memory->destroy(comall); + memory->destroy(dipole); + memory->destroy(dipoleall); + maxchunk = nchunk; + memory->create(massproc,maxchunk,"dipole/chunk:massproc"); + memory->create(masstotal,maxchunk,"dipole/chunk:masstotal"); + memory->create(chrgproc,maxchunk,"dipole/chunk:chrgproc"); + memory->create(chrgtotal,maxchunk,"dipole/chunk:chrgtotal"); + memory->create(com,maxchunk,3,"dipole/chunk:com"); + memory->create(comall,maxchunk,3,"dipole/chunk:comall"); + memory->create(dipole,maxchunk,4,"dipole/chunk:dipole"); + memory->create(dipoleall,maxchunk,4,"dipole/chunk:dipoleall"); + array = dipoleall; +} + +/* ---------------------------------------------------------------------- + memory usage of local data +------------------------------------------------------------------------- */ + +double ComputeDipoleChunk::memory_usage() +{ + double bytes = (bigint) maxchunk * 2 * sizeof(double); + bytes += (bigint) maxchunk * 2*3 * sizeof(double); + bytes += (bigint) maxchunk * 2*4 * sizeof(double); + return bytes; +} diff --git a/src/compute_dipole_chunk.h b/src/compute_dipole_chunk.h new file mode 100644 index 0000000000..b0605b1c35 --- /dev/null +++ b/src/compute_dipole_chunk.h @@ -0,0 +1,77 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef COMPUTE_CLASS + +ComputeStyle(dipole/chunk,ComputeDipoleChunk) + +#else + +#ifndef LMP_COMPUTE_DIPOLE_CHUNK_H +#define LMP_COMPUTE_DIPOLE_CHUNK_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeDipoleChunk : public Compute { + public: + ComputeDipoleChunk(class LAMMPS *, int, char **); + ~ComputeDipoleChunk(); + void init(); + void compute_array(); + + void lock_enable(); + void lock_disable(); + int lock_length(); + void lock(class Fix *, bigint, bigint); + void unlock(class Fix *); + + double memory_usage(); + + private: + int nchunk,maxchunk; + char *idchunk; + class ComputeChunkAtom *cchunk; + + double *massproc,*masstotal; + double *chrgproc,*chrgtotal; + double **com,**comall; + double **dipole,**dipoleall; + int usecenter; + + void allocate(); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Chunk/atom compute does not exist for compute dipole/chunk + +Self-explanatory. + +E: Compute dipole/chunk does not use chunk/atom compute + +The style of the specified compute is not chunk/atom. + +*/ diff --git a/src/info.cpp b/src/info.cpp index fe6f354fcd..9db3dd01de 100644 --- a/src/info.cpp +++ b/src/info.cpp @@ -27,6 +27,7 @@ #include "fix.h" #include "force.h" #include "pair.h" +#include "pair_hybrid.h" #include "group.h" #include "input.h" #include "modify.h" @@ -281,6 +282,13 @@ void Info::command(int narg, char **arg) fprintf(out,"Atoms = " BIGINT_FORMAT ", types = %d, style = %s\n", atom->natoms, atom->ntypes, force->pair_style); + if (force->pair && strstr(force->pair_style,"hybrid")) { + PairHybrid *hybrid = (PairHybrid *)force->pair; + fprintf(out,"Hybrid sub-styles:"); + for (int i=0; i < hybrid->nstyles; ++i) + fprintf(out," %s", hybrid->keywords[i]); + fputc('\n',out); + } if (atom->molecular > 0) { const char *msg; msg = force->bond_style ? force->bond_style : "none"; diff --git a/src/library.cpp b/src/library.cpp index da81194e3a..e51813413d 100644 --- a/src/library.cpp +++ b/src/library.cpp @@ -417,8 +417,9 @@ void lammps_gather_atoms(void *ptr, char *name, int flag = 0; if (lmp->atom->tag_enable == 0 || lmp->atom->tag_consecutive() == 0) flag = 1; if (lmp->atom->natoms > MAXSMALLINT) flag = 1; - if (flag && lmp->comm->me == 0) { - lmp->error->warning(FLERR,"Library error in lammps_gather_atoms"); + if (flag) { + if (lmp->comm->me == 0) + lmp->error->warning(FLERR,"Library error in lammps_gather_atoms"); return; } @@ -506,8 +507,9 @@ void lammps_scatter_atoms(void *ptr, char *name, if (lmp->atom->tag_enable == 0 || lmp->atom->tag_consecutive() == 0) flag = 1; if (lmp->atom->natoms > MAXSMALLINT) flag = 1; if (lmp->atom->map_style == 0) flag = 1; - if (flag && lmp->comm->me == 0) { - lmp->error->warning(FLERR,"Library error in lammps_scatter_atoms"); + if (flag) { + if (lmp->comm->me == 0) + lmp->error->warning(FLERR,"Library error in lammps_scatter_atoms"); return; } diff --git a/src/pair.cpp b/src/pair.cpp index e3c4a1aa36..c169177ea9 100644 --- a/src/pair.cpp +++ b/src/pair.cpp @@ -38,11 +38,11 @@ #include "suffix.h" #include "atom_masks.h" #include "memory.h" +#include "math_const.h" #include "error.h" using namespace LAMMPS_NS; - -#define EWALD_F 1.12837917 +using namespace MathConst; enum{NONE,RLINEAR,RSQ,BMP}; @@ -56,8 +56,6 @@ Pair::Pair(LAMMPS *lmp) : Pointers(lmp) { instance_me = instance_total++; - THIRD = 1.0/3.0; - eng_vdwl = eng_coul = 0.0; comm_forward = comm_reverse = comm_reverse_off = 0; @@ -385,7 +383,7 @@ void Pair::init_tables(double cut_coul, double *cut_respa) ftable[i] = qqrd2e/r * fgamma; etable[i] = qqrd2e/r * egamma; } else { - ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2); + ftable[i] = qqrd2e/r * (derfc + MY_ISPI4*grij*expm2); etable[i] = qqrd2e/r * derfc; } } else { @@ -397,9 +395,9 @@ void Pair::init_tables(double cut_coul, double *cut_respa) etable[i] = qqrd2e/r * egamma; vtable[i] = qqrd2e/r * fgamma; } else { - ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2 - 1.0); + ftable[i] = qqrd2e/r * (derfc + MY_ISPI4*grij*expm2 - 1.0); etable[i] = qqrd2e/r * derfc; - vtable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2); + vtable[i] = qqrd2e/r * (derfc + MY_ISPI4*grij*expm2); } if (rsq_lookup.f > cut_respa[2]*cut_respa[2]) { if (rsq_lookup.f < cut_respa[3]*cut_respa[3]) { @@ -408,7 +406,7 @@ void Pair::init_tables(double cut_coul, double *cut_respa) ctable[i] = qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw); } else { if (msmflag) ftable[i] = qqrd2e/r * fgamma; - else ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2); + else ftable[i] = qqrd2e/r * (derfc + MY_ISPI4*grij*expm2); ctable[i] = qqrd2e/r; } } @@ -481,7 +479,7 @@ void Pair::init_tables(double cut_coul, double *cut_respa) f_tmp = qqrd2e/r * fgamma; e_tmp = qqrd2e/r * egamma; } else { - f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2); + f_tmp = qqrd2e/r * (derfc + MY_ISPI4*grij*expm2); e_tmp = qqrd2e/r * derfc; } } else { @@ -492,9 +490,9 @@ void Pair::init_tables(double cut_coul, double *cut_respa) e_tmp = qqrd2e/r * egamma; v_tmp = qqrd2e/r * fgamma; } else { - f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2 - 1.0); + f_tmp = qqrd2e/r * (derfc + MY_ISPI4*grij*expm2 - 1.0); e_tmp = qqrd2e/r * derfc; - v_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2); + v_tmp = qqrd2e/r * (derfc + MY_ISPI4*grij*expm2); } if (rsq_lookup.f > cut_respa[2]*cut_respa[2]) { if (rsq_lookup.f < cut_respa[3]*cut_respa[3]) { @@ -503,7 +501,7 @@ void Pair::init_tables(double cut_coul, double *cut_respa) c_tmp = qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw); } else { if (msmflag) f_tmp = qqrd2e/r * fgamma; - else f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2); + else f_tmp = qqrd2e/r * (derfc + MY_ISPI4*grij*expm2); c_tmp = qqrd2e/r; } } diff --git a/src/pair.h b/src/pair.h index 5f33f04e46..b2f0ece46e 100644 --- a/src/pair.h +++ b/src/pair.h @@ -219,8 +219,6 @@ class Pair : protected Pointers { typedef union {int i; float f;} union_int_float_t; - double THIRD; - int vflag_fdotr; int maxeatom,maxvatom; diff --git a/src/pair_hybrid.h b/src/pair_hybrid.h index 4962fd4fbe..4d224dafc3 100644 --- a/src/pair_hybrid.h +++ b/src/pair_hybrid.h @@ -31,6 +31,7 @@ class PairHybrid : public Pair { friend class FixOMP; friend class Force; friend class Respa; + friend class Info; public: PairHybrid(class LAMMPS *); virtual ~PairHybrid(); diff --git a/src/random_park.h b/src/random_park.h index 2f300db089..dfa0b24312 100644 --- a/src/random_park.h +++ b/src/random_park.h @@ -19,7 +19,6 @@ namespace LAMMPS_NS { class RanPark : protected Pointers { - friend class Set; public: RanPark(class LAMMPS *, int); double uniform(); diff --git a/src/write_restart.cpp b/src/write_restart.cpp index 9305c34ba6..7b3e341b94 100644 --- a/src/write_restart.cpp +++ b/src/write_restart.cpp @@ -111,28 +111,30 @@ void WriteRestart::command(int narg, char **arg) // init entire system since comm->exchange is done // comm::init needs neighbor::init needs pair::init needs kspace::init, etc - if (comm->me == 0 && screen) - fprintf(screen,"System init for write_restart ...\n"); - lmp->init(); + if (noinit == 0) { + if (comm->me == 0 && screen) + fprintf(screen,"System init for write_restart ...\n"); + lmp->init(); - // move atoms to new processors before writing file - // enforce PBC in case atoms are outside box - // call borders() to rebuild atom map since exchange() destroys map - // NOTE: removed call to setup_pre_exchange - // used to be needed by fixShearHistory for granular - // to move history info from neigh list to atoms between runs - // but now that is done via FIx::post_run() - // don't think any other fix needs this or should do it - // e.g. fix evaporate should not delete more atoms + // move atoms to new processors before writing file + // enforce PBC in case atoms are outside box + // call borders() to rebuild atom map since exchange() destroys map + // NOTE: removed call to setup_pre_exchange + // used to be needed by fixShearHistory for granular + // to move history info from neigh list to atoms between runs + // but now that is done via FIx::post_run() + // don't think any other fix needs this or should do it + // e.g. fix evaporate should not delete more atoms - // modify->setup_pre_exchange(); - if (domain->triclinic) domain->x2lamda(atom->nlocal); - domain->pbc(); - domain->reset_box(); - comm->setup(); - comm->exchange(); - comm->borders(); - if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost); + // modify->setup_pre_exchange(); + if (domain->triclinic) domain->x2lamda(atom->nlocal); + domain->pbc(); + domain->reset_box(); + comm->setup(); + comm->exchange(); + comm->borders(); + if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost); + } // write single restart file @@ -220,6 +222,9 @@ void WriteRestart::multiproc_options(int multiproc_caller, int mpiioflag_caller, else filewriter = 0; iarg += 2; + } else if (strcmp(arg[iarg],"noinit") == 0) { + noinit = 1; + iarg++; } else error->all(FLERR,"Illegal write_restart command"); } } diff --git a/src/write_restart.h b/src/write_restart.h index df710759ab..2253977bda 100644 --- a/src/write_restart.h +++ b/src/write_restart.h @@ -36,6 +36,7 @@ class WriteRestart : protected Pointers { int me,nprocs; FILE *fp; bigint natoms; // natoms (sum of nlocal) to write into file + int noinit; int multiproc; // 0 = proc 0 writes for all // else # of procs writing files