From bea8025b9345e2a1fa084fcc1244f5d0a135caa6 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Wed, 13 Dec 2006 00:34:21 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@186 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/GRANULAR/Install.csh | 8 +- src/GRANULAR/{fix_insert.cpp => fix_pour.cpp} | 95 +++-- src/GRANULAR/{fix_insert.h => fix_pour.h} | 10 +- src/GRANULAR/pair_gran_history.cpp | 8 +- src/GRANULAR/pair_gran_history.h | 1 - src/GRANULAR/style_granular.h | 4 +- src/KSPACE/pair_lj_cut_coul_long.cpp | 2 +- src/fix_deposit.cpp | 369 ++++++++++++++++++ src/fix_deposit.h | 40 ++ src/fix_gravity.h | 2 +- src/fix_indent.cpp | 8 +- src/fix_shear_history.h | 2 +- src/region_block.h | 2 +- src/region_cylinder.h | 2 +- src/region_prism.h | 2 - src/style.h | 2 + src/style_granular.h | 4 +- src/temp_full.cpp | 55 +-- src/temp_full.h | 3 + src/temp_partial.cpp | 52 ++- src/temp_partial.h | 2 + src/temp_ramp.cpp | 45 ++- src/temp_ramp.h | 2 + src/temp_region.cpp | 46 ++- src/temperature.cpp | 67 ++-- src/temperature.h | 12 +- 26 files changed, 658 insertions(+), 187 deletions(-) rename src/GRANULAR/{fix_insert.cpp => fix_pour.cpp} (85%) rename src/GRANULAR/{fix_insert.h => fix_pour.h} (92%) create mode 100644 src/fix_deposit.cpp create mode 100644 src/fix_deposit.h diff --git a/src/GRANULAR/Install.csh b/src/GRANULAR/Install.csh index 104514fed4..ccc65ca90f 100644 --- a/src/GRANULAR/Install.csh +++ b/src/GRANULAR/Install.csh @@ -9,8 +9,8 @@ if ($1 == 1) then cp atom_granular.cpp .. cp fix_freeze.cpp .. cp fix_gran_diag.cpp .. - cp fix_insert.cpp .. cp fix_nve_gran.cpp .. + cp fix_pour.cpp .. cp fix_shear_history.cpp .. cp fix_wall_gran.cpp .. cp pair_gran_hertzian.cpp .. @@ -20,8 +20,8 @@ if ($1 == 1) then cp atom_granular.h .. cp fix_freeze.h .. cp fix_gran_diag.h .. - cp fix_insert.h .. cp fix_nve_gran.h .. + cp fix_pour.h .. # cp fix_shear_history.h .. cp fix_wall_gran.h .. cp pair_gran_hertzian.h .. @@ -36,8 +36,8 @@ else if ($1 == 0) then rm ../atom_granular.cpp rm ../fix_freeze.cpp rm ../fix_gran_diag.cpp - rm ../fix_insert.cpp rm ../fix_nve_gran.cpp + rm ../fix_pour.cpp rm ../fix_shear_history.cpp rm ../fix_wall_gran.cpp rm ../pair_gran_hertzian.cpp @@ -47,8 +47,8 @@ else if ($1 == 0) then rm ../atom_granular.h rm ../fix_freeze.h rm ../fix_gran_diag.h - rm ../fix_insert.h rm ../fix_nve_gran.h + rm ../fix_pour.h # rm ../fix_shear_history.h rm ../fix_wall_gran.h rm ../pair_gran_hertzian.h diff --git a/src/GRANULAR/fix_insert.cpp b/src/GRANULAR/fix_pour.cpp similarity index 85% rename from src/GRANULAR/fix_insert.cpp rename to src/GRANULAR/fix_pour.cpp index 94c6e3ec33..775b58005a 100644 --- a/src/GRANULAR/fix_insert.cpp +++ b/src/GRANULAR/fix_pour.cpp @@ -14,14 +14,14 @@ #include "math.h" #include "stdlib.h" #include "string.h" -#include "fix_insert.h" +#include "fix_pour.h" #include "atom.h" #include "force.h" #include "update.h" +#include "comm.h" #include "modify.h" #include "fix_gravity.h" #include "fix_shear_history.h" -#include "neighbor.h" #include "domain.h" #include "region.h" #include "region_block.h" @@ -34,12 +34,12 @@ /* ---------------------------------------------------------------------- */ -FixInsert::FixInsert(int narg, char **arg) : Fix(narg, arg) +FixPour::FixPour(int narg, char **arg) : Fix(narg, arg) { - if (narg < 6) error->all("Illegal fix insert command"); + if (narg < 6) error->all("Illegal fix pour command"); if (atom->check_style("granular") == 0) - error->all("Must use fix insert with atom style granular"); + error->all("Must use fix pour with atom style granular"); // required args @@ -64,34 +64,34 @@ FixInsert::FixInsert(int narg, char **arg) : Fix(narg, arg) int iarg = 6; while (iarg < narg) { if (strcmp(arg[iarg],"region") == 0) { - if (iarg+2 > narg) error->all("Illegal fix insert command"); + if (iarg+2 > narg) error->all("Illegal fix pour command"); for (iregion = 0; iregion < domain->nregion; iregion++) if (strcmp(arg[iarg+1],domain->regions[iregion]->id) == 0) break; if (iregion == domain->nregion) - error->all("Fix insert region ID does not exist"); + error->all("Fix pour region ID does not exist"); iarg += 2; } else if (strcmp(arg[iarg],"diam") == 0) { - if (iarg+3 > narg) error->all("Illegal fix insert command"); + if (iarg+3 > narg) error->all("Illegal fix pour command"); radius_lo = 0.5 * atof(arg[iarg+1]); radius_hi = 0.5 * atof(arg[iarg+2]); iarg += 3; } else if (strcmp(arg[iarg],"dens") == 0) { - if (iarg+3 > narg) error->all("Illegal fix insert command"); + if (iarg+3 > narg) error->all("Illegal fix pour command"); density_lo = atof(arg[iarg+1]); density_hi = atof(arg[iarg+2]); iarg += 3; } else if (strcmp(arg[iarg],"vol") == 0) { - if (iarg+3 > narg) error->all("Illegal fix insert command"); + if (iarg+3 > narg) error->all("Illegal fix pour command"); volfrac = atof(arg[iarg+1]); maxattempt = atoi(arg[iarg+2]); iarg += 3; } else if (strcmp(arg[iarg],"rate") == 0) { - if (iarg+2 > narg) error->all("Illegal fix insert command"); + if (iarg+2 > narg) error->all("Illegal fix pour command"); rate = atof(arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"vel") == 0) { if (force->dimension == 3) { - if (iarg+6 > narg) error->all("Illegal fix insert command"); + if (iarg+6 > narg) error->all("Illegal fix pour command"); vxlo = atof(arg[iarg+1]); vxhi = atof(arg[iarg+2]); vylo = atof(arg[iarg+3]); @@ -99,24 +99,24 @@ FixInsert::FixInsert(int narg, char **arg) : Fix(narg, arg) vz = atof(arg[iarg+5]); iarg += 6; } else { - if (iarg+4 > narg) error->all("Illegal fix insert command"); + if (iarg+4 > narg) error->all("Illegal fix pour command"); vxlo = atof(arg[iarg+1]); vxhi = atof(arg[iarg+2]); vy = atof(arg[iarg+3]); vz = 0.0; iarg += 4; } - } else error->all("Illegal fix insert command"); + } else error->all("Illegal fix pour command"); } // error check that a valid region was specified - if (iregion == -1) error->all("Must specify a region in fix insert"); + if (iregion == -1) error->all("Must specify a region in fix pour"); // error checks on region if (domain->regions[iregion]->interior == 0) - error->all("Must use region with side = in with fix insert"); + error->all("Must use region with side = in with fix pour"); if (strcmp(domain->regions[iregion]->style,"block") == 0) { region_style = 1; @@ -139,15 +139,15 @@ FixInsert::FixInsert(int narg, char **arg) : Fix(narg, arg) zlo = ((RegCylinder *) domain->regions[iregion])->lo; zhi = ((RegCylinder *) domain->regions[iregion])->hi; if (axis != 'z') - error->all("Must use a z-axis cylinder with fix insert"); + error->all("Must use a z-axis cylinder with fix pour"); if (xc-rc < domain->boxxlo || xc+rc > domain->boxxhi || yc-rc < domain->boxylo || yc+rc > domain->boxyhi || zlo < domain->boxzlo || zhi > domain->boxzhi) error->all("Insertion region extends outside simulation box"); - } else error->all("Must use a block or cylinder region with fix insert"); + } else error->all("Must use a block or cylinder region with fix pour"); if (region_style == 2 && force->dimension == 2) - error->all("Must use a block region with fix insert for 2d simulations"); + error->all("Must use a block region with fix pour for 2d simulations"); // random number generator, same for all procs @@ -219,7 +219,7 @@ FixInsert::FixInsert(int narg, char **arg) : Fix(narg, arg) /* ---------------------------------------------------------------------- */ -FixInsert::~FixInsert() +FixPour::~FixPour() { delete random; delete [] recvcounts; @@ -228,7 +228,7 @@ FixInsert::~FixInsert() /* ---------------------------------------------------------------------- */ -int FixInsert::setmask() +int FixPour::setmask() { int mask = 0; mask |= PRE_EXCHANGE; @@ -237,7 +237,7 @@ int FixInsert::setmask() /* ---------------------------------------------------------------------- */ -void FixInsert::init() +void FixPour::init() { // insure gravity fix exists // for 3d must point in -z, for 2d must point in -y @@ -247,7 +247,7 @@ void FixInsert::init() for (ifix = 0; ifix < modify->nfix; ifix++) if (strcmp(modify->fix[ifix]->style,"gravity") == 0) break; if (ifix == modify->nfix) - error->all("Must use fix gravity with fix insert"); + error->all("Must use fix gravity with fix pour"); double phi = ((FixGravity *) modify->fix[ifix])->phi; double theta = ((FixGravity *) modify->fix[ifix])->theta; @@ -260,11 +260,11 @@ void FixInsert::init() if (force->dimension == 3) { if (fabs(xgrav) > EPSILON || fabs(ygrav) > EPSILON || fabs(zgrav+1.0) > EPSILON) - error->all("Gravity must point in -z to use with fix insert in 3d"); + error->all("Gravity must point in -z to use with fix pour in 3d"); } else { if (fabs(xgrav) > EPSILON || fabs(ygrav+1.0) > EPSILON || fabs(zgrav) > EPSILON) - error->all("Gravity must point in -y to use with fix insert in 2d"); + error->all("Gravity must point in -y to use with fix pour in 2d"); } // check if a shear history fix exists @@ -279,7 +279,7 @@ void FixInsert::init() perform particle insertion ------------------------------------------------------------------------- */ -void FixInsert::pre_exchange() +void FixPour::pre_exchange() { int i; @@ -316,9 +316,9 @@ void FixInsert::pre_exchange() // xnear is for atoms from all procs + atoms to be inserted double **xmine = - memory->create_2d_double_array(ncount,4,"fix_insert:xmine"); + memory->create_2d_double_array(ncount,4,"fix_pour:xmine"); double **xnear = - memory->create_2d_double_array(nprevious+nnew,4,"fix_insert:xnear"); + memory->create_2d_double_array(nprevious+nnew,4,"fix_pour:xnear"); int nnear = nprevious; // setup for allgatherv @@ -402,16 +402,17 @@ void FixInsert::pre_exchange() if (nnear - nprevious < nnew && me == 0) error->warning("Less insertions than requested"); - // add new atoms in my sub-box to my arrays - // initialize info about the atoms + // check if new atom is in my sub-box or above it if I'm highest proc + // if so, add to my list via create_one() + // initialize info about the atom // type, diameter, density set from fix parameters // group mask set to "all" plus fix group // z velocity set to what velocity would be if particle // had fallen from top of insertion region // this gives continuous stream of atoms - // set npartner for new atoms to 0 (assume not touching any others) + // set npartner for new atom to 0 (assume not touching any others) - int m; + int m,flag; double denstmp,vxtmp,vytmp,vztmp; double g = 1.0; @@ -431,9 +432,19 @@ void FixInsert::pre_exchange() vztmp = 0.0; } + flag = 0; if (xtmp >= domain->subxlo && xtmp < domain->subxhi && ytmp >= domain->subylo && ytmp < domain->subyhi && - ztmp >= domain->subzlo && ztmp < domain->subzhi) { + ztmp >= domain->subzlo && ztmp < domain->subzhi) flag = 1; + else if (force->dimension == 3 && ztmp >= domain->boxzhi && + comm->myloc[2] == comm->procgrid[2]-1 && + xtmp >= domain->subxlo && xtmp < domain->subxhi && + ytmp >= domain->subylo && ytmp < domain->subyhi) flag = 1; + else if (force->dimension == 2 && ytmp >= domain->boxyhi && + comm->myloc[1] == comm->procgrid[1]-1 && + xtmp >= domain->subxlo && xtmp < domain->subxhi) flag = 1; + + if (flag) { atom->create_one(ntype,xtmp,ytmp,ztmp); m = atom->nlocal - 1; atom->type[m] = ntype; @@ -452,15 +463,17 @@ void FixInsert::pre_exchange() } } - // tag # of new particles grow beyond all previous atoms + // set tag # of new particles beyond all previous atoms // reset global natoms // if global map exists, reset it - atom->tag_extend(); - atom->natoms += nnear - nprevious; - if (atom->map_style) { - atom->map_init(); - atom->map_set(); + if (atom->tag_enable) { + atom->tag_extend(); + atom->natoms += nnear - nprevious; + if (atom->map_style) { + atom->map_init(); + atom->map_set(); + } } // free local memory @@ -480,7 +493,7 @@ void FixInsert::pre_exchange() use maximum diameter for inserted particle ------------------------------------------------------------------------- */ -int FixInsert::overlap(int i) +int FixPour::overlap(int i) { double delta = radius_hi + atom->radius[i]; double **x = atom->x; @@ -508,7 +521,7 @@ int FixInsert::overlap(int i) /* ---------------------------------------------------------------------- */ -void FixInsert::xyz_random(double h, double &x, double &y, double &z) +void FixPour::xyz_random(double h, double &x, double &y, double &z) { if (force->dimension == 3) { if (region_style == 1) { diff --git a/src/GRANULAR/fix_insert.h b/src/GRANULAR/fix_pour.h similarity index 92% rename from src/GRANULAR/fix_insert.h rename to src/GRANULAR/fix_pour.h index f8b89157a1..2fea9b3576 100644 --- a/src/GRANULAR/fix_insert.h +++ b/src/GRANULAR/fix_pour.h @@ -11,21 +11,21 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#ifndef FIX_INSERT_H -#define FIX_INSERT_H +#ifndef FIX_POUR_H +#define FIX_POUR_H #include "fix.h" class RanPark; -class FixInsert : public Fix { +class FixPour : public Fix { friend class PairGranHistory; friend class PairGranHertzian; friend class PairGranNoHistory; public: - FixInsert(int, char **); - ~FixInsert(); + FixPour(int, char **); + ~FixPour(); int setmask(); void init(); void pre_exchange(); diff --git a/src/GRANULAR/pair_gran_history.cpp b/src/GRANULAR/pair_gran_history.cpp index f3a44b7c05..8c34eb6dc4 100644 --- a/src/GRANULAR/pair_gran_history.cpp +++ b/src/GRANULAR/pair_gran_history.cpp @@ -26,7 +26,7 @@ #include "update.h" #include "modify.h" #include "fix.h" -#include "fix_insert.h" +#include "fix_pour.h" #include "comm.h" #include "memory.h" #include "neighbor.h" @@ -390,10 +390,10 @@ void PairGranHistory::init_style() MPI_Allreduce(&mine,&maxrad_dynamic,1,MPI_DOUBLE,MPI_MAX,world); for (i = 0; i < modify->nfix; i++) - if (strcmp(modify->fix[i]->style,"insert") == 0) + if (strcmp(modify->fix[i]->style,"pour") == 0) maxrad_dynamic = - MAX(maxrad_dynamic,((FixInsert *) modify->fix[i])->radius_hi); - + MAX(maxrad_dynamic,((FixPour *) modify->fix[i])->radius_hi); + double maxrad_frozen = 0.0; for (i = 0; i < nlocal; i++) if (mask[i] & freeze_group_bit) diff --git a/src/GRANULAR/pair_gran_history.h b/src/GRANULAR/pair_gran_history.h index d8be75d258..60f2329082 100644 --- a/src/GRANULAR/pair_gran_history.h +++ b/src/GRANULAR/pair_gran_history.h @@ -20,7 +20,6 @@ class PairGranHistory : public Pair { friend class Neighbor; friend class FixWallGran; friend class FixGranDiag; - friend class FixInsert; public: PairGranHistory(); diff --git a/src/GRANULAR/style_granular.h b/src/GRANULAR/style_granular.h index 5b571b79a6..b221551953 100644 --- a/src/GRANULAR/style_granular.h +++ b/src/GRANULAR/style_granular.h @@ -22,8 +22,8 @@ AtomStyle(granular,AtomGranular) #ifdef FixInclude #include "fix_freeze.h" #include "fix_gran_diag.h" -#include "fix_insert.h" #include "fix_nve_gran.h" +#include "fix_pour.h" #include "fix_shear_history.h" #include "fix_wall_gran.h" #endif @@ -31,8 +31,8 @@ AtomStyle(granular,AtomGranular) #ifdef FixClass FixStyle(freeze,FixFreeze) FixStyle(gran/diag,FixGranDiag) -FixStyle(insert,FixInsert) FixStyle(nve/gran,FixNVEGran) +FixStyle(pour,FixPour) FixStyle(SHEAR_HISTORY,FixShearHistory) FixStyle(wall/gran,FixWallGran) #endif diff --git a/src/KSPACE/pair_lj_cut_coul_long.cpp b/src/KSPACE/pair_lj_cut_coul_long.cpp index d306df3a4d..ced25d88e6 100644 --- a/src/KSPACE/pair_lj_cut_coul_long.cpp +++ b/src/KSPACE/pair_lj_cut_coul_long.cpp @@ -686,7 +686,7 @@ double PairLJCutCoulLong::init_one(int i, int j) lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0); lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0); lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0); - + if (offset_flag) { double ratio = sigma[i][j] / cut_lj[i][j]; offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0)); diff --git a/src/fix_deposit.cpp b/src/fix_deposit.cpp new file mode 100644 index 0000000000..2b96796d7f --- /dev/null +++ b/src/fix_deposit.cpp @@ -0,0 +1,369 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + www.cs.sandia.gov/~sjplimp/lammps.html + Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories + + 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 "math.h" +#include "stdlib.h" +#include "string.h" +#include "fix_deposit.h" +#include "atom.h" +#include "force.h" +#include "update.h" +#include "comm.h" +#include "domain.h" +#include "lattice.h" +#include "region.h" +#include "random_park.h" +#include "memory.h" +#include "error.h" + +/* ---------------------------------------------------------------------- */ + +FixDeposit::FixDeposit(int narg, char **arg) : Fix(narg, arg) +{ + if (narg < 7) error->all("Illegal fix deposit command"); + + // required args + + ninsert = atoi(arg[3]); + ntype = atoi(arg[4]); + nfreq = atoi(arg[5]); + seed = atoi(arg[6]); + + // set defaults + + iregion = -1; + globalflag = localflag = 0; + lo = hi = deltasq = 0.0; + nearsq = 0.0; + maxattempt = 10; + rateflag = 0; + vxlo = vxhi = vylo = vyhi = vzlo = vzhi = 0.0; + scaleflag = 1; + + // read options from end of input line + + options(narg-7,&arg[7]); + + // error check on region + + if (iregion == -1) error->all("Must specify a region in fix deposit"); + if (domain->regions[iregion]->interior == 0) + error->all("Must use region with side = in with fix deposit"); + + // store extent of region + + xlo = domain->regions[iregion]->extent_xlo; + xhi = domain->regions[iregion]->extent_xhi; + ylo = domain->regions[iregion]->extent_ylo; + yhi = domain->regions[iregion]->extent_yhi; + zlo = domain->regions[iregion]->extent_zlo; + zhi = domain->regions[iregion]->extent_zhi; + + // setup scaling + + if (scaleflag && domain->lattice == NULL) + error->all("Use of fix deposit with undefined lattice"); + + double xscale,yscale,zscale; + if (scaleflag) { + xscale = domain->lattice->xlattice; + yscale = domain->lattice->ylattice; + zscale = domain->lattice->zlattice; + } + else xscale = yscale = zscale = 1.0; + + // apply scaling to all input parameters with dist/vel units + + if (force->dimension == 2) { + lo *= yscale; + hi *= yscale; + rate *= yscale; + } else { + lo *= zscale; + hi *= zscale; + rate *= zscale; + } + deltasq *= xscale*xscale; + nearsq *= xscale*xscale; + vxlo *= xscale; + vxhi *= xscale; + vylo *= yscale; + vyhi *= yscale; + vzlo *= zscale; + vzhi *= zscale; + + // store extent of region + + xlo = domain->regions[iregion]->extent_xlo; + xhi = domain->regions[iregion]->extent_xhi; + ylo = domain->regions[iregion]->extent_ylo; + yhi = domain->regions[iregion]->extent_yhi; + zlo = domain->regions[iregion]->extent_zlo; + zhi = domain->regions[iregion]->extent_zhi; + + // random number generator, same for all procs + + random = new RanPark(seed); + + // set up reneighboring + + force_reneighbor = 1; + next_reneighbor = update->ntimestep + 1; + nfirst = next_reneighbor; + ninserted = 0; +} + +/* ---------------------------------------------------------------------- */ + +FixDeposit::~FixDeposit() +{ + delete random; +} + +/* ---------------------------------------------------------------------- */ + +int FixDeposit::setmask() +{ + int mask = 0; + mask |= PRE_EXCHANGE; + return mask; +} + +/* ---------------------------------------------------------------------- + perform particle insertion +------------------------------------------------------------------------- */ + +void FixDeposit::pre_exchange() +{ + int flag,flagall; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + + // just return if should not be called on this timestep + + if (next_reneighbor != update->ntimestep) return; + + // compute current offset = bottom of insertion volume + + double offset = 0.0; + if (rateflag) offset = (update->ntimestep - nfirst) * update->dt * rate; + + // attempt an insertion until successful + + int success = 0; + int attempt = 0; + while (attempt < maxattempt) { + attempt++; + + // choose random position for new atom within region + + xtmp = xlo + random->uniform() * (xhi-xlo); + ytmp = ylo + random->uniform() * (yhi-ylo); + ztmp = zlo + random->uniform() * (zhi-zlo); + while (domain->regions[iregion]->match(xtmp,ytmp,ztmp) == 0) { + xtmp = xlo + random->uniform() * (xhi-xlo); + ytmp = ylo + random->uniform() * (yhi-ylo); + ztmp = zlo + random->uniform() * (zhi-zlo); + } + + // adjust vertical coord by offset + + if (force->dimension == 2) ytmp += offset; + else ztmp += offset; + + // if global, reset vertical coord to be lo-hi above highest atom + // if local, reset vertical coord to be lo-hi above highest "nearby" atom + // local computation computes lateral distance between 2 particles w/ PBC + + if (globalflag || localflag) { + int dim; + double max,maxall,delx,dely,delz,rsq; + + if (force->dimension == 2) { + dim = 1; + max = domain->boxylo; + } else { + dim = 2; + max = domain->boxzlo; + } + + double **x = atom->x; + int nlocal = atom->nlocal; + for (int i = 0; i < nlocal; i++) { + if (localflag) { + delx = xtmp - x[i][0]; + dely = ytmp - x[i][1]; + delz = 0.0; + domain->minimum_image(&delx,&dely,&delz); + if (force->dimension == 2) rsq = delx*delx; + else rsq = delx*delx + dely*dely; + if (rsq > deltasq) continue; + } + if (x[i][dim] > max) max = x[i][dim]; + } + + MPI_Allreduce(&max,&maxall,1,MPI_DOUBLE,MPI_MAX,world); + if (force->dimension == 2) + ytmp = maxall + lo + random->uniform()*(hi-lo); + else + ztmp = maxall + lo + random->uniform()*(hi-lo); + } + + // now have final xtmp,ytmp,ztmp + // if distance to any atom is less than near, try again + + double **x = atom->x; + int nlocal = atom->nlocal; + + flag = 0; + for (int i = 0; i < nlocal; i++) { + delx = xtmp - x[i][0]; + dely = ytmp - x[i][1]; + delz = ztmp - x[i][2]; + domain->minimum_image(&delx,&dely,&delz); + rsq = delx*delx + dely*dely + delz*delz; + if (rsq < nearsq) flag = 1; + } + MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_MAX,world); + if (flagall) continue; + + // insertion will proceed + // choose random velocity for new atom + + double vxtmp = vxlo + random->uniform() * (vxhi-vxlo); + double vytmp = vylo + random->uniform() * (vyhi-vylo); + double vztmp = vzlo + random->uniform() * (vzhi-vzlo); + + // check if new atom is in my sub-box or above it if I'm highest proc + // if so, add to my list via create_one() + // initialize info about the atoms + // set group mask to "all" plus fix group + + flag = 0; + if (xtmp >= domain->subxlo && xtmp < domain->subxhi && + ytmp >= domain->subylo && ytmp < domain->subyhi && + ztmp >= domain->subzlo && ztmp < domain->subzhi) flag = 1; + else if (force->dimension == 3 && ztmp >= domain->boxzhi && + comm->myloc[2] == comm->procgrid[2]-1 && + xtmp >= domain->subxlo && xtmp < domain->subxhi && + ytmp >= domain->subylo && ytmp < domain->subyhi) flag = 1; + else if (force->dimension == 2 && ytmp >= domain->boxyhi && + comm->myloc[1] == comm->procgrid[1]-1 && + xtmp >= domain->subxlo && xtmp < domain->subxhi) flag = 1; + + if (flag) { + atom->create_one(ntype,xtmp,ytmp,ztmp); + int m = atom->nlocal - 1; + atom->type[m] = ntype; + atom->mask[m] = 1 | groupbit; + atom->v[m][0] = vxtmp; + atom->v[m][1] = vytmp; + atom->v[m][2] = vztmp; + } + MPI_Allreduce(&flag,&success,1,MPI_INT,MPI_MAX,world); + break; + } + + // warn if not successful b/c too many attempts or no proc owned particle + + if (comm->me == 0) + if (success == 0) + error->warning("Particle deposition was unsuccessful"); + + // set tag # of new particle beyond all previous atoms + // reset global natoms + // if global map exists, reset it + + if (success && atom->tag_enable) { + atom->tag_extend(); + atom->natoms += 1; + if (atom->map_style) { + atom->map_init(); + atom->map_set(); + } + } + + // next timestep to insert + + if (ninserted < ninsert) next_reneighbor += nfreq; + else next_reneighbor = 0; +} + +/* ---------------------------------------------------------------------- + parse optional parameters at end of input line +------------------------------------------------------------------------- */ + +void FixDeposit::options(int narg, char **arg) +{ + if (narg < 0) error->all("Illegal fix indent command"); + + int iarg = 0; + while (iarg < narg) { + if (strcmp(arg[iarg],"region") == 0) { + if (iarg+2 > narg) error->all("Illegal fix deposit command"); + for (iregion = 0; iregion < domain->nregion; iregion++) + if (strcmp(arg[iarg+1],domain->regions[iregion]->id) == 0) break; + if (iregion == domain->nregion) + error->all("Fix deposit region ID does not exist"); + iarg += 2; + } else if (strcmp(arg[iarg],"global") == 0) { + if (iarg+3 > narg) error->all("Illegal fix deposit command"); + globalflag = 1; + localflag = 0; + lo = atof(arg[iarg+1]); + hi = atof(arg[iarg+2]); + iarg += 3; + } else if (strcmp(arg[iarg],"local") == 0) { + if (iarg+4 > narg) error->all("Illegal fix deposit command"); + localflag = 1; + globalflag = 0; + lo = atof(arg[iarg+1]); + hi = atof(arg[iarg+2]); + deltasq = atof(arg[iarg+3])*atof(arg[iarg+3]); + iarg += 4; + } else if (strcmp(arg[iarg],"near") == 0) { + if (iarg+2 > narg) error->all("Illegal fix deposit command"); + nearsq = atof(arg[iarg+1])*atof(arg[iarg+1]); + iarg += 2; + } else if (strcmp(arg[iarg],"attempt") == 0) { + if (iarg+2 > narg) error->all("Illegal fix deposit command"); + maxattempt = atoi(arg[iarg+1]); + iarg += 2; + } else if (strcmp(arg[iarg],"rate") == 0) { + if (iarg+2 > narg) error->all("Illegal fix deposit command"); + rateflag = 1; + rate = atof(arg[iarg+1]); + iarg += 2; + } else if (strcmp(arg[iarg],"vx") == 0) { + if (iarg+3 > narg) error->all("Illegal fix deposit command"); + vxlo = atof(arg[iarg+1]); + vxhi = atof(arg[iarg+2]); + iarg += 3; + } else if (strcmp(arg[iarg],"vy") == 0) { + if (iarg+3 > narg) error->all("Illegal fix deposit command"); + vylo = atof(arg[iarg+1]); + vyhi = atof(arg[iarg+2]); + iarg += 3; + } else if (strcmp(arg[iarg],"vz") == 0) { + if (iarg+3 > narg) error->all("Illegal fix deposit command"); + vzlo = atof(arg[iarg+1]); + vzhi = atof(arg[iarg+2]); + iarg += 3; + } else if (strcmp(arg[iarg],"units") == 0) { + if (iarg+2 > narg) error->all("Illegal fix deposit command"); + if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0; + else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1; + else error->all("Illegal fix deposit command"); + iarg += 2; + } else error->all("Illegal fix deposit command"); + } +} diff --git a/src/fix_deposit.h b/src/fix_deposit.h new file mode 100644 index 0000000000..a0fee52ff5 --- /dev/null +++ b/src/fix_deposit.h @@ -0,0 +1,40 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + www.cs.sandia.gov/~sjplimp/lammps.html + Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories + + 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. +------------------------------------------------------------------------- */ + +#ifndef FIX_DEPOSIT_H +#define FIX_DEPOSIT_H + +#include "fix.h" + +class RanPark; + +class FixDeposit : public Fix { + public: + FixDeposit(int, char **); + ~FixDeposit(); + int setmask(); + void pre_exchange(); + + private: + int ninsert,ntype,nfreq,seed; + int iregion,globalflag,localflag,maxattempt,rateflag,scaleflag; + double lo,hi,deltasq,nearsq,rate; + double vxlo,vxhi,vylo,vyhi,vzlo,vzhi; + double xlo,xhi,ylo,yhi,zlo,zhi; + int nfirst,ninserted; + RanPark *random; + + void options(int, char **); +}; + +#endif diff --git a/src/fix_gravity.h b/src/fix_gravity.h index 1f597415fd..5f5301e703 100644 --- a/src/fix_gravity.h +++ b/src/fix_gravity.h @@ -17,7 +17,7 @@ #include "fix.h" class FixGravity : public Fix { - friend class FixInsert; + friend class FixPour; public: FixGravity(int, char **); diff --git a/src/fix_indent.cpp b/src/fix_indent.cpp index 511f43cfa4..ff491f2ef5 100644 --- a/src/fix_indent.cpp +++ b/src/fix_indent.cpp @@ -38,13 +38,13 @@ FixIndent::FixIndent(int narg, char **arg) : Fix(narg, arg) if (narg < 4) error->all("Illegal fix indent command"); k = atof(arg[3]); - // set input line defaults + // set defaults istyle = NONE; vx = vy = vz = 0.0; - scaleflag = 1; radflag = 0; r0_start = 0.0; + scaleflag = 1; // read options from end of input line @@ -63,7 +63,7 @@ FixIndent::FixIndent(int narg, char **arg) : Fix(narg, arg) } else xscale = yscale = zscale = 1.0; - // apply scaling to indenter force constant, geometry, and velocity + // apply scaling to force constant, velocity, and geometry k /= xscale; k3 = k/3.0; @@ -290,7 +290,7 @@ int FixIndent::thermo_compute(double *values) } /* ---------------------------------------------------------------------- - parse optional parameters at end of fix indent input line + parse optional parameters at end of input line ------------------------------------------------------------------------- */ void FixIndent::options(int narg, char **arg) diff --git a/src/fix_shear_history.h b/src/fix_shear_history.h index 17a375027c..f9e9be0fd2 100644 --- a/src/fix_shear_history.h +++ b/src/fix_shear_history.h @@ -18,7 +18,7 @@ class FixShearHistory : public Fix { friend class Neighbor; - friend class FixInsert; + friend class FixPour; public: FixShearHistory(int, char **); diff --git a/src/region_block.h b/src/region_block.h index 056a2e65d6..5feae6267f 100644 --- a/src/region_block.h +++ b/src/region_block.h @@ -17,7 +17,7 @@ #include "region.h" class RegBlock : public Region { - friend class FixInsert; + friend class FixPour; public: RegBlock(int, char **); diff --git a/src/region_cylinder.h b/src/region_cylinder.h index 0355f4c998..2ae2220e84 100644 --- a/src/region_cylinder.h +++ b/src/region_cylinder.h @@ -17,7 +17,7 @@ #include "region.h" class RegCylinder : public Region { - friend class FixInsert; + friend class FixPour; public: RegCylinder(int, char **); diff --git a/src/region_prism.h b/src/region_prism.h index 39f3c4565f..e758295946 100644 --- a/src/region_prism.h +++ b/src/region_prism.h @@ -17,8 +17,6 @@ #include "region.h" class RegPrism : public Region { - friend class FixInsert; - public: RegPrism(int, char **); ~RegPrism() {} diff --git a/src/style.h b/src/style.h index d75feff684..b1988d4760 100644 --- a/src/style.h +++ b/src/style.h @@ -97,6 +97,7 @@ DumpStyle(xyz,DumpXYZ) #include "fix_centro.h" #include "fix_com.h" #include "fix_drag.h" +#include "fix_deposit.h" #include "fix_efield.h" #include "fix_energy.h" #include "fix_enforce2d.h" @@ -142,6 +143,7 @@ FixStyle(aveforce,FixAveForce) FixStyle(CENTRO,FixCentro) FixStyle(com,FixCOM) FixStyle(drag,FixDrag) +FixStyle(deposit,FixDeposit) FixStyle(efield,FixEfield) FixStyle(ENERGY,FixEnergy) FixStyle(enforce2d,FixEnforce2D) diff --git a/src/style_granular.h b/src/style_granular.h index 5b571b79a6..b221551953 100644 --- a/src/style_granular.h +++ b/src/style_granular.h @@ -22,8 +22,8 @@ AtomStyle(granular,AtomGranular) #ifdef FixInclude #include "fix_freeze.h" #include "fix_gran_diag.h" -#include "fix_insert.h" #include "fix_nve_gran.h" +#include "fix_pour.h" #include "fix_shear_history.h" #include "fix_wall_gran.h" #endif @@ -31,8 +31,8 @@ AtomStyle(granular,AtomGranular) #ifdef FixClass FixStyle(freeze,FixFreeze) FixStyle(gran/diag,FixGranDiag) -FixStyle(insert,FixInsert) FixStyle(nve/gran,FixNVEGran) +FixStyle(pour,FixPour) FixStyle(SHEAR_HISTORY,FixShearHistory) FixStyle(wall/gran,FixWallGran) #endif diff --git a/src/temp_full.cpp b/src/temp_full.cpp index d55400d204..10b6516149 100644 --- a/src/temp_full.cpp +++ b/src/temp_full.cpp @@ -15,21 +15,32 @@ #include "temp_full.h" #include "atom.h" #include "force.h" +#include "group.h" #include "error.h" /* ---------------------------------------------------------------------- */ -TempFull::TempFull(int narg, char **arg) : Temperature(narg, arg) {} +TempFull::TempFull(int narg, char **arg) : Temperature(narg, arg) +{ + options(narg-3,&arg[3]); +} /* ---------------------------------------------------------------------- */ void TempFull::init() { - count_atoms(); count_fix(); - dof = force->dimension * ncount; + recount(); +} + +/* ---------------------------------------------------------------------- */ + +void TempFull::recount() +{ + double natoms = group->count(igroup); + dof = force->dimension * natoms; dof -= extra_dof + fix_dof; - if (ncount > 0) tfactor = force->mvv2e / (dof * force->boltz); + if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } @@ -58,6 +69,7 @@ double TempFull::compute() } MPI_Allreduce(&t,&t_total,1,MPI_DOUBLE,MPI_SUM,world); + if (dynamic) recount(); t_total *= tfactor; return t_total; } @@ -74,33 +86,22 @@ void TempFull::tensor() int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; + int mass_require = atom->mass_require; double massone,t[6]; for (i = 0; i < 6; i++) t[i] = 0.0; - if (atom->mass_require) { - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - massone = mass[type[i]]; - t[0] += massone * v[i][0]*v[i][0]; - t[1] += massone * v[i][1]*v[i][1]; - t[2] += massone * v[i][2]*v[i][2]; - t[3] += massone * v[i][0]*v[i][1]; - t[4] += massone * v[i][0]*v[i][2]; - t[5] += massone * v[i][1]*v[i][2]; - } - } else { - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - massone = rmass[i]; - t[0] += massone * v[i][0]*v[i][0]; - t[1] += massone * v[i][1]*v[i][1]; - t[2] += massone * v[i][2]*v[i][2]; - t[3] += massone * v[i][0]*v[i][1]; - t[4] += massone * v[i][0]*v[i][2]; - t[5] += massone * v[i][1]*v[i][2]; - } - } + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + if (mass_require) massone = mass[type[i]]; + else massone = rmass[i]; + t[0] += massone * v[i][0]*v[i][0]; + t[1] += massone * v[i][1]*v[i][1]; + t[2] += massone * v[i][2]*v[i][2]; + t[3] += massone * v[i][0]*v[i][1]; + t[4] += massone * v[i][0]*v[i][2]; + t[5] += massone * v[i][1]*v[i][2]; + } MPI_Allreduce(&t,&ke_tensor,6,MPI_DOUBLE,MPI_SUM,world); for (i = 0; i < 6; i++) ke_tensor[i] *= force->mvv2e; diff --git a/src/temp_full.h b/src/temp_full.h index 355258626f..d413f1ec8c 100644 --- a/src/temp_full.h +++ b/src/temp_full.h @@ -23,6 +23,9 @@ class TempFull : public Temperature { void init(); double compute(); void tensor(); + + private: + void recount(); }; #endif diff --git a/src/temp_partial.cpp b/src/temp_partial.cpp index f8236d19d1..219880642c 100644 --- a/src/temp_partial.cpp +++ b/src/temp_partial.cpp @@ -23,6 +23,8 @@ TempPartial::TempPartial(int narg, char **arg) : Temperature(narg, arg) { + options(narg-6,&arg[6]); + xflag = atoi(arg[3]); yflag = atoi(arg[4]); zflag = atoi(arg[5]); @@ -32,11 +34,18 @@ TempPartial::TempPartial(int narg, char **arg) : Temperature(narg, arg) void TempPartial::init() { - count_atoms(); count_fix(); - dof = (xflag+yflag+zflag) * ncount; + recount(); +} + +/* ---------------------------------------------------------------------- */ + +void TempPartial::recount() +{ + double natoms = group->count(igroup); + dof = (xflag+yflag+zflag) * natoms; dof -= extra_dof + fix_dof; - if (ncount > 0) tfactor = force->mvv2e / (dof * force->boltz); + if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } @@ -46,17 +55,27 @@ double TempPartial::compute() { double **v = atom->v; double *mass = atom->mass; + double *rmass = atom->rmass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double t = 0.0; - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) - t += (xflag*v[i][0]*v[i][0] + yflag*v[i][1]*v[i][1] + - zflag*v[i][2]*v[i][2]) * mass[type[i]]; + if (atom->mass_require) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + t += (xflag*v[i][0]*v[i][0] + yflag*v[i][1]*v[i][1] + + zflag*v[i][2]*v[i][2]) * mass[type[i]]; + } else { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + t += (xflag*v[i][0]*v[i][0] + yflag*v[i][1]*v[i][1] + + zflag*v[i][2]*v[i][2]) * rmass[i]; + } + MPI_Allreduce(&t,&t_total,1,MPI_DOUBLE,MPI_SUM,world); + if (dynamic) recount(); t_total *= tfactor; return t_total; } @@ -69,22 +88,25 @@ void TempPartial::tensor() double **v = atom->v; double *mass = atom->mass; + double *rmass = atom->rmass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; + int mass_require = atom->mass_require; - double rmass,t[6]; + double massone,t[6]; for (i = 0; i < 6; i++) t[i] = 0.0; for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - rmass = mass[type[i]]; - t[0] += rmass * xflag*v[i][0]*v[i][0]; - t[1] += rmass * yflag*v[i][1]*v[i][1]; - t[2] += rmass * zflag*v[i][2]*v[i][2]; - t[3] += rmass * xflag*yflag*v[i][0]*v[i][1]; - t[4] += rmass * xflag*zflag*v[i][0]*v[i][2]; - t[5] += rmass * yflag*zflag*v[i][1]*v[i][2]; + if (mass_require) massone = mass[type[i]]; + else massone = rmass[i]; + t[0] += massone * xflag*v[i][0]*v[i][0]; + t[1] += massone * yflag*v[i][1]*v[i][1]; + t[2] += massone * zflag*v[i][2]*v[i][2]; + t[3] += massone * xflag*yflag*v[i][0]*v[i][1]; + t[4] += massone * xflag*zflag*v[i][0]*v[i][2]; + t[5] += massone * yflag*zflag*v[i][1]*v[i][2]; } MPI_Allreduce(&t,&ke_tensor,6,MPI_DOUBLE,MPI_SUM,world); diff --git a/src/temp_partial.h b/src/temp_partial.h index 8916a09c5f..3d834f4bfe 100644 --- a/src/temp_partial.h +++ b/src/temp_partial.h @@ -26,6 +26,8 @@ class TempPartial : public Temperature { private: int xflag,yflag,zflag; + + void recount(); }; #endif diff --git a/src/temp_ramp.cpp b/src/temp_ramp.cpp index 7964068d46..d544bb3d58 100644 --- a/src/temp_ramp.cpp +++ b/src/temp_ramp.cpp @@ -17,6 +17,7 @@ #include "temp_ramp.h" #include "atom.h" #include "force.h" +#include "group.h" #include "error.h" #define MIN(A,B) ((A) < (B)) ? (A) : (B) @@ -26,6 +27,8 @@ TempRamp::TempRamp(int narg, char **arg) : Temperature(narg, arg) { + options(narg-9,&arg[9]); + if (strcmp(arg[3],"vx") == 0) v_dim = 0; else if (strcmp(arg[3],"vy") == 0) v_dim = 1; else if (strcmp(arg[3],"vz") == 0) v_dim = 2; @@ -63,11 +66,18 @@ TempRamp::TempRamp(int narg, char **arg) : Temperature(narg, arg) void TempRamp::init() { - count_atoms(); count_fix(); - dof = force->dimension * ncount; + recount(); +} + +/* ---------------------------------------------------------------------- */ + +void TempRamp::recount() +{ + double natoms = group->count(igroup); + dof = force->dimension * natoms; dof -= extra_dof + fix_dof; - if (ncount > 0) tfactor = force->mvv2e / (dof * force->boltz); + if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } @@ -80,9 +90,11 @@ double TempRamp::compute() double **x = atom->x; double **v = atom->v; double *mass = atom->mass; + double *rmass = atom->rmass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; + int mass_require = atom->mass_require; double t = 0.0; for (int i = 0; i < nlocal; i++) @@ -95,11 +107,15 @@ double TempRamp::compute() vtmp[1] = v[i][1]; vtmp[2] = v[i][2]; vtmp[v_dim] -= vramp; - t += (vtmp[0]*vtmp[0] + vtmp[1]*vtmp[1] + vtmp[2]*vtmp[2]) * - mass[type[i]]; + if (mass_require) + t += (vtmp[0]*vtmp[0] + vtmp[1]*vtmp[1] + vtmp[2]*vtmp[2]) * + mass[type[i]]; + else + t += (vtmp[0]*vtmp[0] + vtmp[1]*vtmp[1] + vtmp[2]*vtmp[2]) * rmass[i]; } MPI_Allreduce(&t,&t_total,1,MPI_DOUBLE,MPI_SUM,world); + if (dynamic) recount(); t_total *= tfactor; return t_total; } @@ -114,11 +130,13 @@ void TempRamp::tensor() double **x = atom->x; double **v = atom->v; double *mass = atom->mass; + double *rmass = atom->rmass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; + int mass_require = atom->mass_require; - double rmass,t[6]; + double massone,t[6]; for (i = 0; i < 6; i++) t[i] = 0.0; for (i = 0; i < nlocal; i++) @@ -132,13 +150,14 @@ void TempRamp::tensor() vtmp[2] = v[i][2]; vtmp[v_dim] -= vramp; - rmass = mass[type[i]]; - t[0] += rmass * vtmp[0]*vtmp[0]; - t[1] += rmass * vtmp[1]*vtmp[1]; - t[2] += rmass * vtmp[2]*vtmp[2]; - t[3] += rmass * vtmp[0]*vtmp[1]; - t[4] += rmass * vtmp[0]*vtmp[2]; - t[5] += rmass * vtmp[1]*vtmp[2]; + if (mass_require) massone = mass[type[i]]; + else massone = rmass[i]; + t[0] += massone * vtmp[0]*vtmp[0]; + t[1] += massone * vtmp[1]*vtmp[1]; + t[2] += massone * vtmp[2]*vtmp[2]; + t[3] += massone * vtmp[0]*vtmp[1]; + t[4] += massone * vtmp[0]*vtmp[2]; + t[5] += massone * vtmp[1]*vtmp[2]; } MPI_Allreduce(&t,&ke_tensor,6,MPI_DOUBLE,MPI_SUM,world); diff --git a/src/temp_ramp.h b/src/temp_ramp.h index 5cf7f6201b..ab8674d20d 100644 --- a/src/temp_ramp.h +++ b/src/temp_ramp.h @@ -29,6 +29,8 @@ class TempRamp : public Temperature { double coord_lo,coord_hi; int v_dim; double v_lo,v_hi; + + void recount(); }; #endif diff --git a/src/temp_region.cpp b/src/temp_region.cpp index 344b0ee116..15766e45f3 100644 --- a/src/temp_region.cpp +++ b/src/temp_region.cpp @@ -24,6 +24,8 @@ TempRegion::TempRegion(int narg, char **arg) : Temperature(narg, arg) { + options(narg-4,&arg[4]); + for (iregion = 0; iregion < domain->nregion; iregion++) if (strcmp(arg[3],domain->regions[iregion]->id) == 0) break; if (iregion == domain->nregion) @@ -44,19 +46,30 @@ double TempRegion::compute() double **x = atom->x; double **v = atom->v; double *mass = atom->mass; + double *rmass = atom->rmass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; int count = 0; double t = 0.0; - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit && - domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2])) { - count++; - t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * - mass[type[i]]; - } + + if (atom->mass_require) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit && + domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2])) { + count++; + t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * + mass[type[i]]; + } + } else { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit && + domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2])) { + count++; + t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i]; + } + } double tarray[2],tarray_all[2]; tarray[0] = count; @@ -77,23 +90,26 @@ void TempRegion::tensor() double **x = atom->x; double **v = atom->v; double *mass = atom->mass; + double *rmass = atom->rmass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; + int mass_require = atom->mass_require; - double rmass,t[6]; + double massone,t[6]; for (i = 0; i < 6; i++) t[i] = 0.0; for (i = 0; i < nlocal; i++) if (mask[i] & groupbit && domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2])) { - rmass = mass[type[i]]; - t[0] += rmass * v[i][0]*v[i][0]; - t[1] += rmass * v[i][1]*v[i][1]; - t[2] += rmass * v[i][2]*v[i][2]; - t[3] += rmass * v[i][0]*v[i][1]; - t[4] += rmass * v[i][0]*v[i][2]; - t[5] += rmass * v[i][1]*v[i][2]; + if (mass_require) massone = mass[type[i]]; + else massone = rmass[i]; + t[0] += massone * v[i][0]*v[i][0]; + t[1] += massone * v[i][1]*v[i][1]; + t[2] += massone * v[i][2]*v[i][2]; + t[3] += massone * v[i][0]*v[i][1]; + t[4] += massone * v[i][0]*v[i][2]; + t[5] += massone * v[i][1]*v[i][2]; } MPI_Allreduce(&t,&ke_tensor,6,MPI_DOUBLE,MPI_SUM,world); diff --git a/src/temperature.cpp b/src/temperature.cpp index 86fc68fa30..8af6a24d64 100644 --- a/src/temperature.cpp +++ b/src/temperature.cpp @@ -45,32 +45,7 @@ Temperature::Temperature(int narg, char **arg) // set modify defaults extra_dof = 3; - - // set input line defaults - - scaleflag = 1; - - // read options from end of input line - - if (strcmp(style,"full") == 0) options(narg-3,&arg[3]); - else if (strcmp(style,"partial") == 0) options(narg-6,&arg[6]); - else if (strcmp(style,"ramp") == 0) options(narg-9,&arg[9]); - else if (strcmp(style,"region") == 0) options(narg-4,&arg[4]); - - // set scaling for RAMP style - - if (strcmp(style,"ramp") == 0) { - - if (scaleflag && domain->lattice == NULL) - error->all("Use of temperature ramp with undefined lattice"); - - if (scaleflag) { - xscale = domain->lattice->xlattice; - yscale = domain->lattice->ylattice; - zscale = domain->lattice->zlattice; - } - else xscale = yscale = zscale = 1.0; - } + dynamic = 0; } /* ---------------------------------------------------------------------- */ @@ -93,25 +68,16 @@ void Temperature::modify_params(int narg, char **arg) if (iarg+2 > narg) error->all("Illegal temp_modify command"); extra_dof = atoi(arg[iarg+1]); iarg += 2; + } else if (strcmp(arg[iarg],"dynamic") == 0) { + if (iarg+2 > narg) error->all("Illegal temp_modify command"); + if (strcmp(arg[iarg+1],"no") == 0) dynamic = 0; + else if (strcmp(arg[iarg+1],"yes") == 0) dynamic = 1; + else error->all("Illegal temp_modify command"); + iarg += 2; } else error->all("Illegal temp_modify command"); } } -/* ---------------------------------------------------------------------- - ncount = atoms in Temperature group -------------------------------------------------------------------------- */ - -void Temperature::count_atoms() -{ - int *mask = atom->mask; - int nlocal = atom->nlocal; - - int icount = 0; - for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) icount++; - double rcount = icount; - MPI_Allreduce(&rcount,&ncount,1,MPI_DOUBLE,MPI_SUM,world); -} - /* ---------------------------------------------------------------------- count degrees of freedom subtracted by fixes ------------------------------------------------------------------------- */ @@ -131,6 +97,10 @@ void Temperature::options(int narg, char **arg) { if (narg < 0) error->all("Illegal temperature command"); + // option defaults + + scaleflag = 1; + int iarg = 0; while (iarg < narg) { if (strcmp(arg[iarg],"units") == 0) { @@ -141,4 +111,19 @@ void Temperature::options(int narg, char **arg) iarg += 2; } else error->all("Illegal temperature command"); } + + // set scaling for RAMP style + + if (strcmp(style,"ramp") == 0) { + + if (scaleflag && domain->lattice == NULL) + error->all("Use of temperature ramp with undefined lattice"); + + if (scaleflag) { + xscale = domain->lattice->xlattice; + yscale = domain->lattice->ylattice; + zscale = domain->lattice->zlattice; + } + else xscale = yscale = zscale = 1.0; + } } diff --git a/src/temperature.h b/src/temperature.h index f33aa5de84..396ee40cc4 100644 --- a/src/temperature.h +++ b/src/temperature.h @@ -22,21 +22,21 @@ class Temperature : public LAMMPS { int igroup,groupbit; double t_total,dof; double ke_tensor[6]; - double tfactor,ncount; - int extra_dof,fix_dof; - int scaleflag; - double xscale,yscale,zscale; Temperature(int, char **); virtual ~Temperature(); void modify_params(int, char **); - void count_atoms(); void count_fix(); virtual void init() = 0; virtual double compute() = 0; virtual void tensor() = 0; - private: + protected: + double tfactor; + int extra_dof,fix_dof; + int scaleflag,dynamic; + double xscale,yscale,zscale; + void options(int, char **); };