diff --git a/src/fix_deposit.cpp b/src/fix_deposit.cpp deleted file mode 100644 index 9ee902e0c2..0000000000 --- a/src/fix_deposit.cpp +++ /dev/null @@ -1,508 +0,0 @@ -/* ---------------------------------------------------------------------- - 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 "math.h" -#include "stdlib.h" -#include "string.h" -#include "fix_deposit.h" -#include "atom.h" -#include "atom_vec.h" -#include "force.h" -#include "update.h" -#include "modify.h" -#include "fix.h" -#include "comm.h" -#include "domain.h" -#include "lattice.h" -#include "region.h" -#include "random_park.h" -#include "memory.h" -#include "error.h" - -using namespace LAMMPS_NS; -using namespace FixConst; - -/* ---------------------------------------------------------------------- */ - -FixDeposit::FixDeposit(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg) -{ - if (narg < 7) error->all(FLERR,"Illegal fix deposit command"); - - restart_global = 1; - time_depend = 1; - - // required args - - ninsert = force->inumeric(FLERR,arg[3]); - ntype = force->inumeric(FLERR,arg[4]); - nfreq = force->inumeric(FLERR,arg[5]); - seed = force->inumeric(FLERR,arg[6]); - - if (seed <= 0) error->all(FLERR,"Illegal fix deposit command"); - - // set defaults - - iregion = -1; - idregion = NULL; - idnext = 0; - 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; - targetflag = 0; - - // read options from end of input line - - options(narg-7,&arg[7]); - - // error checks on region and its extent being inside simulation box - - if (iregion == -1) error->all(FLERR,"Must specify a region in fix deposit"); - if (domain->regions[iregion]->bboxflag == 0) - error->all(FLERR,"Fix deposit region does not support a bounding box"); - if (domain->regions[iregion]->dynamic_check()) - error->all(FLERR,"Fix deposit region cannot be dynamic"); - - 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; - - if (domain->triclinic == 0) { - if (xlo < domain->boxlo[0] || xhi > domain->boxhi[0] || - ylo < domain->boxlo[1] || yhi > domain->boxhi[1] || - zlo < domain->boxlo[2] || zhi > domain->boxhi[2]) - error->all(FLERR,"Deposition region extends outside simulation box"); - } else { - if (xlo < domain->boxlo_bound[0] || xhi > domain->boxhi_bound[0] || - ylo < domain->boxlo_bound[1] || yhi > domain->boxhi_bound[1] || - zlo < domain->boxlo_bound[2] || zhi > domain->boxhi_bound[2]) - error->all(FLERR,"Deposition region extends outside simulation box"); - } - - // setup scaling - - 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 (domain->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; - tx *= xscale; - ty *= yscale; - tz *= zscale; - - // maxtag_all = current max tag for all atoms - - if (idnext) { - int *tag = atom->tag; - int nlocal = atom->nlocal; - - int maxtag = 0; - for (int i = 0; i < nlocal; i++) maxtag = MAX(maxtag,tag[i]); - MPI_Allreduce(&maxtag,&maxtag_all,1,MPI_INT,MPI_MAX,world); - } - - // random number generator, same for all procs - - random = new RanPark(lmp,seed); - - // set up reneighboring - - force_reneighbor = 1; - next_reneighbor = update->ntimestep + 1; - nfirst = next_reneighbor; - ninserted = 0; -} - -/* ---------------------------------------------------------------------- */ - -FixDeposit::~FixDeposit() -{ - delete random; - delete [] idregion; -} - -/* ---------------------------------------------------------------------- */ - -int FixDeposit::setmask() -{ - int mask = 0; - mask |= PRE_EXCHANGE; - return mask; -} - -/* ---------------------------------------------------------------------- */ - -void FixDeposit::init() -{ - // set index and check validity of region - - iregion = domain->find_region(idregion); - if (iregion == -1) - error->all(FLERR,"Region ID for fix deposit does not exist"); -} - -/* ---------------------------------------------------------------------- - perform particle insertion -------------------------------------------------------------------------- */ - -void FixDeposit::pre_exchange() -{ - int i,j; - int flag,flagall; - double coord[3],lamda[3],delx,dely,delz,rsq; - double *newcoord; - - // 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; - - double *sublo,*subhi; - if (domain->triclinic == 0) { - sublo = domain->sublo; - subhi = domain->subhi; - } else { - sublo = domain->sublo_lamda; - subhi = domain->subhi_lamda; - } - - // attempt an insertion until successful - - int nfix = modify->nfix; - Fix **fix = modify->fix; - - int success = 0; - int attempt = 0; - while (attempt < maxattempt) { - attempt++; - - // choose random position for new atom within region - - coord[0] = xlo + random->uniform() * (xhi-xlo); - coord[1] = ylo + random->uniform() * (yhi-ylo); - coord[2] = zlo + random->uniform() * (zhi-zlo); - while (domain->regions[iregion]->match(coord[0],coord[1],coord[2]) == 0) { - coord[0] = xlo + random->uniform() * (xhi-xlo); - coord[1] = ylo + random->uniform() * (yhi-ylo); - coord[2] = zlo + random->uniform() * (zhi-zlo); - } - - // adjust vertical coord by offset - - if (domain->dimension == 2) coord[1] += offset; - else coord[2] += 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 (domain->dimension == 2) { - dim = 1; - max = domain->boxlo[1]; - } else { - dim = 2; - max = domain->boxlo[2]; - } - - double **x = atom->x; - int nlocal = atom->nlocal; - for (i = 0; i < nlocal; i++) { - if (localflag) { - delx = coord[0] - x[i][0]; - dely = coord[1] - x[i][1]; - delz = 0.0; - domain->minimum_image(delx,dely,delz); - if (domain->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 (domain->dimension == 2) - coord[1] = maxall + lo + random->uniform()*(hi-lo); - else - coord[2] = maxall + lo + random->uniform()*(hi-lo); - } - - // now have final coord - // if distance to any atom is less than near, try again - - double **x = atom->x; - int nlocal = atom->nlocal; - - flag = 0; - for (i = 0; i < nlocal; i++) { - delx = coord[0] - x[i][0]; - dely = coord[1] - x[i][1]; - delz = coord[2] - 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); - - // if we have a sputter target change velocity vector accordingly - if (targetflag) { - double vel = sqrt(vxtmp*vxtmp + vytmp*vytmp + vztmp*vztmp); - delx = tx - coord[0]; - dely = ty - coord[1]; - delz = tz - coord[2]; - double rsq = delx*delx + dely*dely + delz*delz; - if (rsq > 0.0) { - double rinv = sqrt(1.0/rsq); - vxtmp = delx*rinv*vel; - vytmp = dely*rinv*vel; - vztmp = delz*rinv*vel; - } - } - - // 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_atom() - // initialize info about the atoms - // set group mask to "all" plus fix group - - if (domain->triclinic) { - domain->x2lamda(coord,lamda); - newcoord = lamda; - } else newcoord = coord; - - flag = 0; - if (newcoord[0] >= sublo[0] && newcoord[0] < subhi[0] && - newcoord[1] >= sublo[1] && newcoord[1] < subhi[1] && - newcoord[2] >= sublo[2] && newcoord[2] < subhi[2]) flag = 1; - else if (domain->dimension == 3 && newcoord[2] >= domain->boxhi[2] && - comm->myloc[2] == comm->procgrid[2]-1 && - newcoord[0] >= sublo[0] && newcoord[0] < subhi[0] && - newcoord[1] >= sublo[1] && newcoord[1] < subhi[1]) flag = 1; - else if (domain->dimension == 2 && newcoord[1] >= domain->boxhi[1] && - comm->myloc[1] == comm->procgrid[1]-1 && - newcoord[0] >= sublo[0] && newcoord[0] < subhi[0]) flag = 1; - - if (flag) { - atom->avec->create_atom(ntype,coord); - 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; - for (j = 0; j < nfix; j++) - if (fix[j]->create_attribute) fix[j]->set_arrays(m); - } - 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 (!success && comm->me == 0) - error->warning(FLERR,"Particle deposition was unsuccessful",0); - - // reset global natoms - // if idnext, set new atom ID to incremented maxtag_all - // else set new atom ID to value beyond all current atoms - // if global map exists, reset it now instead of waiting for comm - // since adding an atom messes up ghosts - - if (success) { - atom->natoms += 1; - if (atom->tag_enable) { - if (idnext) { - maxtag_all++; - if (atom->nlocal && atom->tag[atom->nlocal-1] == 0) - atom->tag[atom->nlocal-1] = maxtag_all; - } else atom->tag_extend(); - if (atom->map_style) { - atom->nghost = 0; - atom->map_init(); - atom->map_set(); - } - } - } - - // next timestep to insert - // next_reneighbor = 0 if done - - if (success) ninserted++; - 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(FLERR,"Illegal fix indent command"); - - int iarg = 0; - while (iarg < narg) { - if (strcmp(arg[iarg],"region") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command"); - iregion = domain->find_region(arg[iarg+1]); - if (iregion == -1) - error->all(FLERR,"Region ID for fix deposit does not exist"); - int n = strlen(arg[iarg+1]) + 1; - idregion = new char[n]; - strcpy(idregion,arg[iarg+1]); - iarg += 2; - } else if (strcmp(arg[iarg],"id") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command"); - if (strcmp(arg[iarg+1],"max") == 0) idnext = 0; - else if (strcmp(arg[iarg+1],"next") == 0) idnext = 1; - else error->all(FLERR,"Illegal fix deposit command"); - iarg += 2; - } else if (strcmp(arg[iarg],"global") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal fix deposit command"); - globalflag = 1; - localflag = 0; - lo = force->numeric(FLERR,arg[iarg+1]); - hi = force->numeric(FLERR,arg[iarg+2]); - iarg += 3; - } else if (strcmp(arg[iarg],"local") == 0) { - if (iarg+4 > narg) error->all(FLERR,"Illegal fix deposit command"); - localflag = 1; - globalflag = 0; - lo = force->numeric(FLERR,arg[iarg+1]); - hi = force->numeric(FLERR,arg[iarg+2]); - deltasq = force->numeric(FLERR,arg[iarg+3])*force->numeric(FLERR,arg[iarg+3]); - iarg += 4; - } else if (strcmp(arg[iarg],"near") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command"); - nearsq = force->numeric(FLERR,arg[iarg+1])*force->numeric(FLERR,arg[iarg+1]); - iarg += 2; - } else if (strcmp(arg[iarg],"attempt") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command"); - maxattempt = force->inumeric(FLERR,arg[iarg+1]); - iarg += 2; - } else if (strcmp(arg[iarg],"rate") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command"); - rateflag = 1; - rate = force->numeric(FLERR,arg[iarg+1]); - iarg += 2; - } else if (strcmp(arg[iarg],"vx") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal fix deposit command"); - vxlo = force->numeric(FLERR,arg[iarg+1]); - vxhi = force->numeric(FLERR,arg[iarg+2]); - iarg += 3; - } else if (strcmp(arg[iarg],"vy") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal fix deposit command"); - vylo = force->numeric(FLERR,arg[iarg+1]); - vyhi = force->numeric(FLERR,arg[iarg+2]); - iarg += 3; - } else if (strcmp(arg[iarg],"vz") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal fix deposit command"); - vzlo = force->numeric(FLERR,arg[iarg+1]); - vzhi = force->numeric(FLERR,arg[iarg+2]); - iarg += 3; - } else if (strcmp(arg[iarg],"units") == 0) { - if (iarg+2 > narg) error->all(FLERR,"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(FLERR,"Illegal fix deposit command"); - iarg += 2; - } else if (strcmp(arg[iarg],"target") == 0) { - if (iarg+4 > narg) error->all(FLERR,"Illegal fix deposit command"); - tx = force->numeric(FLERR,arg[iarg+1]); - ty = force->numeric(FLERR,arg[iarg+2]); - tz = force->numeric(FLERR,arg[iarg+3]); - targetflag = 1; - iarg += 4; - } else error->all(FLERR,"Illegal fix deposit command"); - } -} - -/* ---------------------------------------------------------------------- - pack entire state of Fix into one write -------------------------------------------------------------------------- */ - -void FixDeposit::write_restart(FILE *fp) -{ - int n = 0; - double list[4]; - list[n++] = random->state(); - list[n++] = ninserted; - list[n++] = nfirst; - list[n++] = next_reneighbor; - - if (comm->me == 0) { - int size = n * sizeof(double); - fwrite(&size,sizeof(int),1,fp); - fwrite(list,sizeof(double),n,fp); - } -} - -/* ---------------------------------------------------------------------- - use state info from restart file to restart the Fix -------------------------------------------------------------------------- */ - -void FixDeposit::restart(char *buf) -{ - int n = 0; - double *list = (double *) buf; - - seed = static_cast (list[n++]); - ninserted = static_cast (list[n++]); - nfirst = static_cast (list[n++]); - next_reneighbor = static_cast (list[n++]); - - random->reset(seed); -} diff --git a/src/fix_deposit.h b/src/fix_deposit.h deleted file mode 100644 index 8b2e9f1052..0000000000 --- a/src/fix_deposit.h +++ /dev/null @@ -1,98 +0,0 @@ -/* -*- 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 FIX_CLASS - -FixStyle(deposit,FixDeposit) - -#else - -#ifndef LMP_FIX_DEPOSIT_H -#define LMP_FIX_DEPOSIT_H - -#include "stdio.h" -#include "fix.h" - -namespace LAMMPS_NS { - -class FixDeposit : public Fix { - public: - FixDeposit(class LAMMPS *, int, char **); - ~FixDeposit(); - int setmask(); - void init(); - void pre_exchange(); - void write_restart(FILE *); - void restart(char *); - - private: - int ninsert,ntype,nfreq,seed; - int iregion,globalflag,localflag,maxattempt,rateflag,scaleflag,targetflag; - char *idregion; - double lo,hi,deltasq,nearsq,rate; - double vxlo,vxhi,vylo,vyhi,vzlo,vzhi; - double xlo,xhi,ylo,yhi,zlo,zhi; - double tx,ty,tz; - int nfirst,ninserted; - int idnext,maxtag_all; - class RanPark *random; - - void options(int, char **); -}; - -} - -#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: Must specify a region in fix deposit - -The region keyword must be specified with this fix. - -E: Fix deposit region does not support a bounding box - -Not all regions represent bounded volumes. You cannot use -such a region with the fix deposit command. - -E: Fix deposit region cannot be dynamic - -Only static regions can be used with fix deposit. - -E: Deposition region extends outside simulation box - -Self-explanatory. - -E: Region ID for fix deposit does not exist - -Self-explanatory. - -W: Particle deposition was unsuccessful - -The fix deposit command was not able to insert as many atoms as -needed. The requested volume fraction may be too high, or other atoms -may be in the insertion region. - -U: Use of fix deposit with undefined lattice - -Must use lattice command with compute fix deposit command if units -option is set to lattice. - -*/ diff --git a/src/fix_efield.cpp b/src/fix_efield.cpp deleted file mode 100644 index f328a3f294..0000000000 --- a/src/fix_efield.cpp +++ /dev/null @@ -1,406 +0,0 @@ -/* ---------------------------------------------------------------------- - 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. -------------------------------------------------------------------------- */ - -/* ---------------------------------------------------------------------- - Contributing author: Christina Payne (Vanderbilt U) - Stan Moore (Sandia) for dipole terms -------------------------------------------------------------------------- */ - -#include "math.h" -#include "string.h" -#include "stdlib.h" -#include "fix_efield.h" -#include "atom.h" -#include "update.h" -#include "domain.h" -#include "comm.h" -#include "modify.h" -#include "force.h" -#include "respa.h" -#include "input.h" -#include "variable.h" -#include "memory.h" -#include "error.h" - -using namespace LAMMPS_NS; -using namespace FixConst; - -enum{NONE,CONSTANT,EQUAL,ATOM}; - -/* ---------------------------------------------------------------------- */ - -FixEfield::FixEfield(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg) -{ - if (narg != 6) error->all(FLERR,"Illegal fix efield command"); - - vector_flag = 1; - scalar_flag = 1; - size_vector = 3; - global_freq = 1; - extvector = 1; - extscalar = 1; - - qe2f = force->qe2f; - xstr = ystr = zstr = NULL; - - if (strstr(arg[3],"v_") == arg[3]) { - int n = strlen(&arg[3][2]) + 1; - xstr = new char[n]; - strcpy(xstr,&arg[3][2]); - } else { - ex = qe2f * force->numeric(FLERR,arg[3]); - xstyle = CONSTANT; - } - - if (strstr(arg[4],"v_") == arg[4]) { - int n = strlen(&arg[4][2]) + 1; - ystr = new char[n]; - strcpy(ystr,&arg[4][2]); - } else { - ey = qe2f * force->numeric(FLERR,arg[4]); - ystyle = CONSTANT; - } - - if (strstr(arg[5],"v_") == arg[5]) { - int n = strlen(&arg[5][2]) + 1; - zstr = new char[n]; - strcpy(zstr,&arg[5][2]); - } else { - ez = qe2f * force->numeric(FLERR,arg[5]); - zstyle = CONSTANT; - } - - // optional args - - estr = NULL; - - int iarg = 6; - while (iarg < narg) { - if (strcmp(arg[iarg],"energy") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix efield command"); - if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) { - int n = strlen(&arg[iarg+1][2]) + 1; - estr = new char[n]; - strcpy(estr,&arg[iarg+1][2]); - } else error->all(FLERR,"Illegal fix efield command"); - iarg += 2; - } else error->all(FLERR,"Illegal fix efield command"); - } - - force_flag = 0; - fsum[0] = fsum[1] = fsum[2] = fsum[3] = 0.0; - - maxatom = 0; - efield = NULL; -} - -/* ---------------------------------------------------------------------- */ - -FixEfield::~FixEfield() -{ - delete [] xstr; - delete [] ystr; - delete [] zstr; - delete [] estr; - memory->destroy(efield); -} - -/* ---------------------------------------------------------------------- */ - -int FixEfield::setmask() -{ - int mask = 0; - mask |= THERMO_ENERGY; - mask |= POST_FORCE; - mask |= POST_FORCE_RESPA; - mask |= MIN_POST_FORCE; - return mask; -} - -/* ---------------------------------------------------------------------- */ - -void FixEfield::init() -{ - qflag = muflag = 0; - if (atom->q_flag) qflag = 1; - if (atom->mu_flag && atom->torque_flag) muflag = 1; - if (!qflag && !muflag) - error->all(FLERR,"Fix efield requires atom attribute q or mu"); - - // check variables - - if (xstr) { - xvar = input->variable->find(xstr); - if (xvar < 0) - error->all(FLERR,"Variable name for fix efield does not exist"); - if (input->variable->equalstyle(xvar)) xstyle = EQUAL; - else if (input->variable->atomstyle(xvar)) xstyle = ATOM; - else error->all(FLERR,"Variable for fix efield is invalid style"); - } - if (ystr) { - yvar = input->variable->find(ystr); - if (yvar < 0) - error->all(FLERR,"Variable name for fix efield does not exist"); - if (input->variable->equalstyle(yvar)) ystyle = EQUAL; - else if (input->variable->atomstyle(yvar)) ystyle = ATOM; - else error->all(FLERR,"Variable for fix efield is invalid style"); - } - if (zstr) { - zvar = input->variable->find(zstr); - if (zvar < 0) - error->all(FLERR,"Variable name for fix efield does not exist"); - if (input->variable->equalstyle(zvar)) zstyle = EQUAL; - else if (input->variable->atomstyle(zvar)) zstyle = ATOM; - else error->all(FLERR,"Variable for fix efield is invalid style"); - } - if (estr) { - evar = input->variable->find(estr); - if (evar < 0) - error->all(FLERR,"Variable name for fix efield does not exist"); - if (input->variable->atomstyle(evar)) estyle = ATOM; - else error->all(FLERR,"Variable for fix efield is invalid style"); - } else estyle = NONE; - - if (xstyle == ATOM || ystyle == ATOM || zstyle == ATOM) - varflag = ATOM; - else if (xstyle == EQUAL || ystyle == EQUAL || zstyle == EQUAL) - varflag = EQUAL; - else varflag = CONSTANT; - - if (muflag && varflag == ATOM) - error->all(FLERR,"Fix efield with dipoles cannot use atom-style variables"); - - if (muflag && update->whichflag == 2 && comm->me == 0) - error->warning(FLERR, - "The minimizer does not re-orient dipoles " - "when using fix efield"); - - if (varflag == CONSTANT && estyle != NONE) - error->all(FLERR,"Cannot use variable energy with " - "constant efield in fix efield"); - if ((varflag == EQUAL || varflag == ATOM) && - update->whichflag == 2 && estyle == NONE) - error->all(FLERR,"Must use variable energy with fix efield"); - - if (strstr(update->integrate_style,"respa")) - nlevels_respa = ((Respa *) update->integrate)->nlevels; -} - -/* ---------------------------------------------------------------------- */ - -void FixEfield::setup(int vflag) -{ - if (strstr(update->integrate_style,"verlet")) - post_force(vflag); - else { - ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1); - post_force_respa(vflag,nlevels_respa-1,0); - ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1); - } -} - -/* ---------------------------------------------------------------------- */ - -void FixEfield::min_setup(int vflag) -{ - post_force(vflag); -} - -/* ---------------------------------------------------------------------- - apply F = qE -------------------------------------------------------------------------- */ - -void FixEfield::post_force(int vflag) -{ - double **f = atom->f; - double *q = atom->q; - int *mask = atom->mask; - tagint *image = atom->image; - int nlocal = atom->nlocal; - - // reallocate efield array if necessary - - if (varflag == ATOM && nlocal > maxatom) { - maxatom = atom->nmax; - memory->destroy(efield); - memory->create(efield,maxatom,4,"efield:efield"); - } - - // fsum[0] = "potential energy" for added force - // fsum[123] = extra force added to atoms - - fsum[0] = fsum[1] = fsum[2] = fsum[3] = 0.0; - force_flag = 0; - - double **x = atom->x; - double fx,fy,fz; - - // constant efield - - if (varflag == CONSTANT) { - double unwrap[3]; - - // charge interactions - // force = qE, potential energy = F dot x in unwrapped coords - - if (qflag) { - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - fx = q[i]*ex; - fy = q[i]*ey; - fz = q[i]*ez; - f[i][0] += fx; - f[i][1] += fy; - f[i][2] += fz; - - domain->unmap(x[i],image[i],unwrap); - fsum[0] -= fx*unwrap[0]+fy*unwrap[1]+fz*unwrap[2]; - fsum[1] += fx; - fsum[2] += fy; - fsum[3] += fz; - } - } - - // dipole interactions - // no force, torque = mu cross E, potential energy = -mu dot E - - if (muflag) { - double **mu = atom->mu; - double **t = atom->torque; - double tx,ty,tz; - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - tx = ez*mu[i][1] - ey*mu[i][2]; - ty = ex*mu[i][2] - ez*mu[i][0]; - tz = ey*mu[i][0] - ex*mu[i][1]; - t[i][0] += tx; - t[i][1] += ty; - t[i][2] += tz; - fsum[0] -= mu[i][0]*ex + mu[i][1]*ey + mu[i][2]*ez; - } - } - - // variable efield, wrap with clear/add - // potential energy = evar if defined, else 0.0 - - } else { - - modify->clearstep_compute(); - - if (xstyle == EQUAL) ex = qe2f * input->variable->compute_equal(xvar); - else if (xstyle == ATOM && efield) - input->variable->compute_atom(xvar,igroup,&efield[0][0],3,0); - if (ystyle == EQUAL) ey = qe2f * input->variable->compute_equal(yvar); - else if (ystyle == ATOM && efield) - input->variable->compute_atom(yvar,igroup,&efield[0][1],3,0); - if (zstyle == EQUAL) ez = qe2f * input->variable->compute_equal(zvar); - else if (zstyle == ATOM && efield) - input->variable->compute_atom(zvar,igroup,&efield[0][2],3,0); - if (estyle == ATOM && efield) - input->variable->compute_atom(evar,igroup,&efield[0][3],4,0); - - modify->addstep_compute(update->ntimestep + 1); - - // charge interactions - // force = qE - - if (qflag) { - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - if (xstyle == ATOM) fx = qe2f * q[i]*efield[i][0]; - else fx = q[i]*ex; - f[i][0] += fx; - fsum[1] += fx; - if (ystyle == ATOM) fy = qe2f * q[i]*efield[i][1]; - else fy = q[i]*ey; - f[i][1] += fy; - fsum[2] += fy; - if (zstyle == ATOM) fz = qe2f * q[i]*efield[i][2]; - else fz = q[i]*ez; - f[i][2] += fz; - fsum[3] += fz; - if (estyle == ATOM) fsum[0] += efield[0][3]; - } - } - - // dipole interactions - // no force, torque = mu cross E - - if (muflag) { - double **mu = atom->mu; - double **t = atom->torque; - double tx,ty,tz; - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - tx = ez*mu[i][1] - ey*mu[i][2]; - ty = ex*mu[i][2] - ez*mu[i][0]; - tz = ey*mu[i][0] - ex*mu[i][1]; - t[i][0] += tx; - t[i][1] += ty; - t[i][2] += tz; - } - } - } -} - -/* ---------------------------------------------------------------------- */ - -void FixEfield::post_force_respa(int vflag, int ilevel, int iloop) -{ - if (ilevel == nlevels_respa-1) post_force(vflag); -} - -/* ---------------------------------------------------------------------- */ - -void FixEfield::min_post_force(int vflag) -{ - post_force(vflag); -} - -/* ---------------------------------------------------------------------- - memory usage of local atom-based array -------------------------------------------------------------------------- */ - -double FixEfield::memory_usage() -{ - double bytes = 0.0; - if (varflag == ATOM) bytes = atom->nmax*4 * sizeof(double); - return bytes; -} - -/* ---------------------------------------------------------------------- - return energy added by fix -------------------------------------------------------------------------- */ - -double FixEfield::compute_scalar(void) -{ - if (force_flag == 0) { - MPI_Allreduce(fsum,fsum_all,4,MPI_DOUBLE,MPI_SUM,world); - force_flag = 1; - } - return fsum_all[0]; -} - -/* ---------------------------------------------------------------------- - return total extra force due to fix -------------------------------------------------------------------------- */ - -double FixEfield::compute_vector(int n) -{ - if (force_flag == 0) { - MPI_Allreduce(fsum,fsum_all,4,MPI_DOUBLE,MPI_SUM,world); - force_flag = 1; - } - return fsum_all[n+1]; -} diff --git a/src/fix_efield.h b/src/fix_efield.h deleted file mode 100644 index 48bd75174c..0000000000 --- a/src/fix_efield.h +++ /dev/null @@ -1,101 +0,0 @@ -/* ---------------------------------------------------------------------- - 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 FIX_CLASS - -FixStyle(efield,FixEfield) - -#else - -#ifndef LMP_FIX_EFIELD_H -#define LMP_FIX_EFIELD_H - -#include "fix.h" - -namespace LAMMPS_NS { - -class FixEfield : public Fix { - public: - FixEfield(class LAMMPS *, int, char **); - ~FixEfield(); - int setmask(); - void init(); - void setup(int); - void min_setup(int); - void post_force(int); - void post_force_respa(int, int, int); - void min_post_force(int); - double memory_usage(); - double compute_scalar(); - double compute_vector(int); - - private: - double ex,ey,ez; - int varflag; - char *xstr,*ystr,*zstr,*estr; - int xvar,yvar,zvar,evar,xstyle,ystyle,zstyle,estyle; - int nlevels_respa; - double qe2f; - double fdotx; - int qflag,muflag; - - int maxatom; - double **efield; - - int force_flag; - double fsum[4],fsum_all[4]; -}; - -} - -#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: Fix efield requires atom attribute q or mu - -Self-explanatory. - -E: Variable name for fix efield does not exist - -Self-explanatory. - -E: Variable for fix efield is invalid style - -Only equal-style or atom-style variables can be used. - -E: Fix efield with dipoles cannot use atom-style variables - -This feature is not yet supported. - -W: The minimizer does not re-orient dipoles when using fix efield - -Self-explanatory. - -E: Cannot use variable energy with constant efield in fix efield - -Self-explanatory. - -E: Must use variable energy with fix efield - -One or more variables are defined for fix efield, which require -variable energy when using the minimizer. - -*/ diff --git a/src/fix_evaporate.cpp b/src/fix_evaporate.cpp deleted file mode 100644 index 5fe5c1412f..0000000000 --- a/src/fix_evaporate.cpp +++ /dev/null @@ -1,386 +0,0 @@ -/* ---------------------------------------------------------------------- - 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 "math.h" -#include "stdlib.h" -#include "string.h" -#include "fix_evaporate.h" -#include "atom.h" -#include "atom_vec.h" -#include "update.h" -#include "domain.h" -#include "region.h" -#include "comm.h" -#include "force.h" -#include "group.h" -#include "random_park.h" -#include "random_mars.h" -#include "memory.h" -#include "error.h" - -using namespace LAMMPS_NS; -using namespace FixConst; - -/* ---------------------------------------------------------------------- */ - -FixEvaporate::FixEvaporate(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg) -{ - if (narg < 7) error->all(FLERR,"Illegal fix evaporate command"); - - scalar_flag = 1; - global_freq = 1; - extscalar = 0; - - nevery = force->inumeric(FLERR,arg[3]); - nflux = force->inumeric(FLERR,arg[4]); - iregion = domain->find_region(arg[5]); - int n = strlen(arg[5]) + 1; - idregion = new char[n]; - strcpy(idregion,arg[5]); - int seed = force->inumeric(FLERR,arg[6]); - - if (nevery <= 0 || nflux <= 0) - error->all(FLERR,"Illegal fix evaporate command"); - if (iregion == -1) - error->all(FLERR,"Region ID for fix evaporate does not exist"); - if (seed <= 0) error->all(FLERR,"Illegal fix evaporate command"); - - // random number generator, same for all procs - - random = new RanPark(lmp,seed); - - // optional args - - molflag = 0; - - int iarg = 7; - while (iarg < narg) { - if (strcmp(arg[iarg],"molecule") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix evaporate command"); - if (strcmp(arg[iarg+1],"no") == 0) molflag = 0; - else if (strcmp(arg[iarg+1],"yes") == 0) molflag = 1; - else error->all(FLERR,"Illegal fix evaporate command"); - iarg += 2; - } else error->all(FLERR,"Illegal fix evaporate command"); - } - - // set up reneighboring - - force_reneighbor = 1; - next_reneighbor = (update->ntimestep/nevery)*nevery + nevery; - ndeleted = 0; - - nmax = 0; - list = NULL; - mark = NULL; -} - -/* ---------------------------------------------------------------------- */ - -FixEvaporate::~FixEvaporate() -{ - delete [] idregion; - delete random; - memory->destroy(list); - memory->destroy(mark); -} - -/* ---------------------------------------------------------------------- */ - -int FixEvaporate::setmask() -{ - int mask = 0; - mask |= PRE_EXCHANGE; - return mask; -} - -/* ---------------------------------------------------------------------- */ - -void FixEvaporate::init() -{ - // set index and check validity of region - - iregion = domain->find_region(idregion); - if (iregion == -1) - error->all(FLERR,"Region ID for fix evaporate does not exist"); - - // check that no deletable atoms are in atom->firstgroup - // deleting such an atom would not leave firstgroup atoms first - - if (atom->firstgroup >= 0) { - int *mask = atom->mask; - int nlocal = atom->nlocal; - int firstgroupbit = group->bitmask[atom->firstgroup]; - - int flag = 0; - for (int i = 0; i < nlocal; i++) - if ((mask[i] & groupbit) && (mask[i] && firstgroupbit)) flag = 1; - - int flagall; - MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world); - - if (flagall) - error->all(FLERR,"Cannot evaporate atoms in atom_modify first group"); - } - - // if molflag not set, warn if any deletable atom has a mol ID - - if (molflag == 0 && atom->molecule_flag) { - int *molecule = atom->molecule; - int *mask = atom->mask; - int nlocal = atom->nlocal; - int flag = 0; - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) - if (molecule[i]) flag = 1; - int flagall; - MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world); - if (flagall && comm->me == 0) - error->warning(FLERR, - "Fix evaporate may delete atom with non-zero molecule ID"); - } - - if (molflag && atom->molecule_flag == 0) - error->all(FLERR, - "Fix evaporate molecule requires atom attribute molecule"); -} - -/* ---------------------------------------------------------------------- - perform particle deletion - done before exchange, borders, reneighbor - so that ghost atoms and neighbor lists will be correct -------------------------------------------------------------------------- */ - -void FixEvaporate::pre_exchange() -{ - int i,j,m,iwhichglobal,iwhichlocal; - int ndel,ndeltopo[4]; - - if (update->ntimestep != next_reneighbor) return; - - // grow list and mark arrays if necessary - - if (atom->nlocal > nmax) { - memory->destroy(list); - memory->destroy(mark); - nmax = atom->nmax; - memory->create(list,nmax,"evaporate:list"); - memory->create(mark,nmax,"evaporate:mark"); - } - - // ncount = # of deletable atoms in region that I own - // nall = # on all procs - // nbefore = # on procs before me - // list[ncount] = list of local indices of atoms I can delete - - double **x = atom->x; - int *mask = atom->mask; - int *tag = atom->tag; - int nlocal = atom->nlocal; - - int ncount = 0; - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) - if (domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2])) - list[ncount++] = i; - - int nall,nbefore; - MPI_Allreduce(&ncount,&nall,1,MPI_INT,MPI_SUM,world); - MPI_Scan(&ncount,&nbefore,1,MPI_INT,MPI_SUM,world); - nbefore -= ncount; - - // ndel = total # of atom deletions, in or out of region - // ndeltopo[1,2,3,4] = ditto for bonds, angles, dihedrals, impropers - // mark[] = 1 if deleted - - ndel = 0; - for (i = 0; i < nlocal; i++) mark[i] = 0; - - // atomic deletions - // choose atoms randomly across all procs and mark them for deletion - // shrink eligible list as my atoms get marked - // keep ndel,ncount,nall,nbefore current after each atom deletion - - if (molflag == 0) { - while (nall && ndel < nflux) { - iwhichglobal = static_cast (nall*random->uniform()); - if (iwhichglobal < nbefore) nbefore--; - else if (iwhichglobal < nbefore + ncount) { - iwhichlocal = iwhichglobal - nbefore; - mark[list[iwhichlocal]] = 1; - list[iwhichlocal] = list[ncount-1]; - ncount--; - } - ndel++; - nall--; - } - - // molecule deletions - // choose one atom in one molecule randomly across all procs - // bcast mol ID and delete all atoms in that molecule on any proc - // update deletion count by total # of atoms in molecule - // shrink list of eligible candidates as any of my atoms get marked - // keep ndel,ndeltopo,ncount,nall,nbefore current after each mol deletion - - } else { - int me,proc,iatom,imolecule,ndelone,ndelall; - int *molecule = atom->molecule; - - ndeltopo[0] = ndeltopo[1] = ndeltopo[2] = ndeltopo[3] = 0; - - while (nall && ndel < nflux) { - - // pick an iatom,imolecule on proc me to delete - - iwhichglobal = static_cast (nall*random->uniform()); - if (iwhichglobal >= nbefore && iwhichglobal < nbefore + ncount) { - iwhichlocal = iwhichglobal - nbefore; - iatom = list[iwhichlocal]; - imolecule = molecule[iatom]; - me = comm->me; - } else me = -1; - - // bcast mol ID to delete all atoms from - // if mol ID > 0, delete any atom in molecule and decrement counters - // if mol ID == 0, delete single iatom - // be careful to delete correct # of bond,angle,etc for newton on or off - - MPI_Allreduce(&me,&proc,1,MPI_INT,MPI_MAX,world); - MPI_Bcast(&imolecule,1,MPI_INT,proc,world); - ndelone = 0; - for (i = 0; i < nlocal; i++) { - if (imolecule && molecule[i] == imolecule) { - mark[i] = 1; - ndelone++; - - if (atom->avec->bonds_allow) { - if (force->newton_bond) ndeltopo[0] += atom->num_bond[i]; - else { - for (j = 0; j < atom->num_bond[i]; j++) { - if (tag[i] < atom->bond_atom[i][j]) ndeltopo[0]++; - } - } - } - if (atom->avec->angles_allow) { - if (force->newton_bond) ndeltopo[1] += atom->num_angle[i]; - else { - for (j = 0; j < atom->num_angle[i]; j++) { - m = atom->map(atom->angle_atom2[i][j]); - if (m >= 0 && m < nlocal) ndeltopo[1]++; - } - } - } - if (atom->avec->dihedrals_allow) { - if (force->newton_bond) ndeltopo[2] += atom->num_dihedral[i]; - else { - for (j = 0; j < atom->num_dihedral[i]; j++) { - m = atom->map(atom->dihedral_atom2[i][j]); - if (m >= 0 && m < nlocal) ndeltopo[2]++; - } - } - } - if (atom->avec->impropers_allow) { - if (force->newton_bond) ndeltopo[3] += atom->num_improper[i]; - else { - for (j = 0; j < atom->num_improper[i]; j++) { - m = atom->map(atom->improper_atom2[i][j]); - if (m >= 0 && m < nlocal) ndeltopo[3]++; - } - } - } - - } else if (me == proc && i == iatom) { - mark[i] = 1; - ndelone++; - } - } - - // remove any atoms marked for deletion from my eligible list - - i = 0; - while (i < ncount) { - if (mark[list[i]]) { - list[i] = list[ncount-1]; - ncount--; - } else i++; - } - - // update ndel,ncount,nall,nbefore - // ndelall is total atoms deleted on this iteration - // ncount is already correct, so resum to get nall and nbefore - - MPI_Allreduce(&ndelone,&ndelall,1,MPI_INT,MPI_SUM,world); - ndel += ndelall; - MPI_Allreduce(&ncount,&nall,1,MPI_INT,MPI_SUM,world); - MPI_Scan(&ncount,&nbefore,1,MPI_INT,MPI_SUM,world); - nbefore -= ncount; - } - } - - // delete my marked atoms - // loop in reverse order to avoid copying marked atoms - - AtomVec *avec = atom->avec; - - for (i = nlocal-1; i >= 0; i--) { - if (mark[i]) { - avec->copy(atom->nlocal-1,i,1); - atom->nlocal--; - } - } - - // reset global natoms and bonds, angles, etc - // if global map exists, reset it now instead of waiting for comm - // since deleting atoms messes up ghosts - - atom->natoms -= ndel; - if (molflag) { - int all[4]; - MPI_Allreduce(ndeltopo,all,4,MPI_INT,MPI_SUM,world); - atom->nbonds -= all[0]; - atom->nangles -= all[1]; - atom->ndihedrals -= all[2]; - atom->nimpropers -= all[3]; - } - - if (ndel && atom->map_style) { - atom->nghost = 0; - atom->map_init(); - atom->map_set(); - } - - // statistics - - ndeleted += ndel; - next_reneighbor = update->ntimestep + nevery; -} - -/* ---------------------------------------------------------------------- - return number of deleted particles -------------------------------------------------------------------------- */ - -double FixEvaporate::compute_scalar() -{ - return 1.0*ndeleted; -} - -/* ---------------------------------------------------------------------- - memory usage of local atom-based arrays -------------------------------------------------------------------------- */ - -double FixEvaporate::memory_usage() -{ - double bytes = 2*nmax * sizeof(int); - return bytes; -} diff --git a/src/fix_evaporate.h b/src/fix_evaporate.h deleted file mode 100644 index f59558103a..0000000000 --- a/src/fix_evaporate.h +++ /dev/null @@ -1,80 +0,0 @@ -/* ---------------------------------------------------------------------- - 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 FIX_CLASS - -FixStyle(evaporate,FixEvaporate) - -#else - -#ifndef LMP_FIX_EVAPORATE_H -#define LMP_FIX_EVAPORATE_H - -#include "fix.h" - -namespace LAMMPS_NS { - -class FixEvaporate : public Fix { - public: - FixEvaporate(class LAMMPS *, int, char **); - ~FixEvaporate(); - int setmask(); - void init(); - void pre_exchange(); - double compute_scalar(); - double memory_usage(); - - private: - int nevery,nflux,iregion; - int molflag; - int ndeleted; - char *idregion; - - int nmax; - int *list,*mark; - - class RanPark *random; -}; - -} - -#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: Region ID for fix evaporate does not exist - -Self-explanatory. - -E: Cannot evaporate atoms in atom_modify first group - -This is a restriction due to the way atoms are organized in -a list to enable the atom_modify first command. - -W: Fix evaporate may delete atom with non-zero molecule ID - -This is probably an error, since you should not delete only one atom -of a molecule. - -E: Fix evaporate molecule requires atom attribute molecule - -The atom style being used does not define a molecule ID. - -*/ diff --git a/src/fix_thermal_conductivity.cpp b/src/fix_thermal_conductivity.cpp deleted file mode 100644 index b653e34864..0000000000 --- a/src/fix_thermal_conductivity.cpp +++ /dev/null @@ -1,330 +0,0 @@ -/* ---------------------------------------------------------------------- - 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. -------------------------------------------------------------------------- */ - -/* ---------------------------------------------------------------------- - Contributing author: Craig Tenney (UND) added support - for swapping atoms of different masses -------------------------------------------------------------------------- */ - -#include "math.h" -#include "mpi.h" -#include "string.h" -#include "stdlib.h" -#include "fix_thermal_conductivity.h" -#include "atom.h" -#include "force.h" -#include "domain.h" -#include "modify.h" -#include "error.h" - -using namespace LAMMPS_NS; -using namespace FixConst; - -#define BIG 1.0e10 - -/* ---------------------------------------------------------------------- */ - -FixThermalConductivity::FixThermalConductivity(LAMMPS *lmp, - int narg, char **arg) : - Fix(lmp, narg, arg) -{ - if (narg < 6) error->all(FLERR,"Illegal fix thermal/conductivity command"); - - MPI_Comm_rank(world,&me); - - nevery = force->inumeric(FLERR,arg[3]); - if (nevery <= 0) error->all(FLERR,"Illegal fix thermal/conductivity command"); - - scalar_flag = 1; - global_freq = nevery; - extscalar = 0; - - if (strcmp(arg[4],"x") == 0) edim = 0; - else if (strcmp(arg[4],"y") == 0) edim = 1; - else if (strcmp(arg[4],"z") == 0) edim = 2; - else error->all(FLERR,"Illegal fix thermal/conductivity command"); - - nbin = force->inumeric(FLERR,arg[5]); - if (nbin % 2 || nbin <= 2) - error->all(FLERR,"Illegal fix thermal/conductivity command"); - - // optional keywords - - nswap = 1; - - int iarg = 6; - while (iarg < narg) { - if (strcmp(arg[iarg],"swap") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal fix thermal/conductivity command"); - nswap = force->inumeric(FLERR,arg[iarg+1]); - if (nswap <= 0) - error->all(FLERR, - "Fix thermal/conductivity swap value must be positive"); - iarg += 2; - } else error->all(FLERR,"Illegal fix thermal/conductivity command"); - } - - // initialize array sizes to nswap+1 so have space to shift values down - - index_lo = new int[nswap+1]; - index_hi = new int[nswap+1]; - ke_lo = new double[nswap+1]; - ke_hi = new double[nswap+1]; - - e_exchange = 0.0; -} - -/* ---------------------------------------------------------------------- */ - -FixThermalConductivity::~FixThermalConductivity() -{ - delete [] index_lo; - delete [] index_hi; - delete [] ke_lo; - delete [] ke_hi; -} - -/* ---------------------------------------------------------------------- */ - -int FixThermalConductivity::setmask() -{ - int mask = 0; - mask |= END_OF_STEP; - return mask; -} - -/* ---------------------------------------------------------------------- */ - -void FixThermalConductivity::init() -{ - // warn if any fix ave/spatial comes after this fix - // can cause glitch in averaging since ave will happen after swap - - int foundme = 0; - for (int i = 0; i < modify->nfix; i++) { - if (modify->fix[i] == this) foundme = 1; - if (foundme && strcmp(modify->fix[i]->style,"ave/spatial") == 0 && me == 0) - error->warning(FLERR, - "Fix thermal/conductivity comes before fix ave/spatial"); - } - - // set bounds of 2 slabs in edim - // only necessary for static box, else re-computed in end_of_step() - // lo bin is always bottom bin - // hi bin is just above half height - - if (domain->box_change == 0) { - prd = domain->prd[edim]; - boxlo = domain->boxlo[edim]; - boxhi = domain->boxhi[edim]; - double binsize = (boxhi-boxlo) / nbin; - slablo_lo = boxlo; - slablo_hi = boxlo + binsize; - slabhi_lo = boxlo + (nbin/2)*binsize; - slabhi_hi = boxlo + (nbin/2+1)*binsize; - } - - periodicity = domain->periodicity[edim]; -} - -/* ---------------------------------------------------------------------- */ - -void FixThermalConductivity::end_of_step() -{ - int i,j,m,insert; - double coord,ke; - MPI_Status status; - struct { - double value; - int proc; - } mine[2],all[2]; - - // if box changes, recompute bounds of 2 slabs in edim - - if (domain->box_change) { - prd = domain->prd[edim]; - boxlo = domain->boxlo[edim]; - boxhi = domain->boxhi[edim]; - double binsize = (boxhi-boxlo) / nbin; - slablo_lo = boxlo; - slablo_hi = boxlo + binsize; - slabhi_lo = boxlo + (nbin/2)*binsize; - slabhi_hi = boxlo + (nbin/2+1)*binsize; - } - - // make 2 lists of up to nswap atoms - // hottest atoms in lo slab, coldest atoms in hi slab (really mid slab) - // lo slab list is sorted by hottest, hi slab is sorted by coldest - // map atoms back into periodic box if necessary - // insert = location in list to insert new atom - - 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; - - nhi = nlo = 0; - - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - coord = x[i][edim]; - if (coord < boxlo && periodicity) coord += prd; - else if (coord >= boxhi && periodicity) coord -= prd; - - if (coord >= slablo_lo && coord < slablo_hi) { - ke = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]; - if (rmass) ke *= 0.5*rmass[i]; - else ke *= 0.5*mass[type[i]]; - if (nlo < nswap || ke > ke_lo[nswap-1]) { - for (insert = nlo-1; insert >= 0; insert--) - if (ke < ke_lo[insert]) break; - insert++; - for (m = nlo-1; m >= insert; m--) { - ke_lo[m+1] = ke_lo[m]; - index_lo[m+1] = index_lo[m]; - } - ke_lo[insert] = ke; - index_lo[insert] = i; - if (nlo < nswap) nlo++; - } - } - - if (coord >= slabhi_lo && coord < slabhi_hi) { - ke = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]; - if (rmass) ke *= 0.5*rmass[i]; - else ke *= 0.5*mass[type[i]]; - if (nhi < nswap || ke < ke_hi[nswap-1]) { - for (insert = nhi-1; insert >= 0; insert--) - if (ke > ke_hi[insert]) break; - insert++; - for (m = nhi-1; m >= insert; m--) { - ke_hi[m+1] = ke_hi[m]; - index_hi[m+1] = index_hi[m]; - } - ke_hi[insert] = ke; - index_hi[insert] = i; - if (nhi < nswap) nhi++; - } - } - } - - // loop over nswap pairs - // pair 2 global atoms at beginning of sorted lo/hi slab lists via Allreduce - // BIG values are for procs with no atom to contribute - // use negative of hottest KE since is doing a MINLOC - // MINLOC also communicates which procs own them - // exchange kinetic energy between the 2 particles - // if I own both particles just swap, else point2point comm of velocities - - double sbuf[4],rbuf[4],vcm[3]; - double eswap = 0.0; - - mine[0].proc = mine[1].proc = me; - int ilo = 0; - int ihi = 0; - - for (m = 0; m < nswap; m++) { - if (ilo < nlo) mine[0].value = -ke_lo[ilo]; - else mine[0].value = BIG; - if (ihi < nhi) mine[1].value = ke_hi[ihi]; - else mine[1].value = BIG; - - MPI_Allreduce(mine,all,2,MPI_DOUBLE_INT,MPI_MINLOC,world); - if (all[0].value == BIG || all[1].value == BIG) continue; - - if (me == all[0].proc && me == all[1].proc) { - i = index_lo[ilo++]; - j = index_hi[ihi++]; - sbuf[0] = v[j][0]; - sbuf[1] = v[j][1]; - sbuf[2] = v[j][2]; - if (rmass) sbuf[3] = rmass[j]; - else sbuf[3] = mass[type[j]]; - rbuf[0] = v[i][0]; - rbuf[1] = v[i][1]; - rbuf[2] = v[i][2]; - if (rmass) rbuf[3] = rmass[i]; - else rbuf[3] = mass[type[i]]; - vcm[0] = (sbuf[3]*sbuf[0] + rbuf[3]*rbuf[0]) / (sbuf[3] + rbuf[3]); - vcm[1] = (sbuf[3]*sbuf[1] + rbuf[3]*rbuf[1]) / (sbuf[3] + rbuf[3]); - vcm[2] = (sbuf[3]*sbuf[2] + rbuf[3]*rbuf[2]) / (sbuf[3] + rbuf[3]); - v[j][0] = 2.0 * vcm[0] - sbuf[0]; - v[j][1] = 2.0 * vcm[1] - sbuf[1]; - v[j][2] = 2.0 * vcm[2] - sbuf[2]; - eswap += sbuf[3] * (vcm[0] * (vcm[0] - sbuf[0]) + - vcm[1] * (vcm[1] - sbuf[1]) + - vcm[2] * (vcm[2] - sbuf[2])); - v[i][0] = 2.0 * vcm[0] - rbuf[0]; - v[i][1] = 2.0 * vcm[1] - rbuf[1]; - v[i][2] = 2.0 * vcm[2] - rbuf[2]; - eswap -= rbuf[3] * (vcm[0] * (vcm[0] - rbuf[0]) + - vcm[1] * (vcm[1] - rbuf[1]) + - vcm[2] * (vcm[2] - rbuf[2])); - - } else if (me == all[0].proc) { - j = index_lo[ilo++]; - sbuf[0] = v[j][0]; - sbuf[1] = v[j][1]; - sbuf[2] = v[j][2]; - if (rmass) sbuf[3] = rmass[j]; - else sbuf[3] = mass[type[j]]; - MPI_Sendrecv(sbuf,4,MPI_DOUBLE,all[1].proc,0, - rbuf,4,MPI_DOUBLE,all[1].proc,0,world,&status); - vcm[0] = (sbuf[3]*sbuf[0] + rbuf[3]*rbuf[0]) / (sbuf[3] + rbuf[3]); - vcm[1] = (sbuf[3]*sbuf[1] + rbuf[3]*rbuf[1]) / (sbuf[3] + rbuf[3]); - vcm[2] = (sbuf[3]*sbuf[2] + rbuf[3]*rbuf[2]) / (sbuf[3] + rbuf[3]); - v[j][0] = 2.0 * vcm[0] - sbuf[0]; - v[j][1] = 2.0 * vcm[1] - sbuf[1]; - v[j][2] = 2.0 * vcm[2] - sbuf[2]; - eswap -= sbuf[3] * (vcm[0] * (vcm[0] - sbuf[0]) + - vcm[1] * (vcm[1] - sbuf[1]) + - vcm[2] * (vcm[2] - sbuf[2])); - - } else if (me == all[1].proc) { - j = index_hi[ihi++]; - sbuf[0] = v[j][0]; - sbuf[1] = v[j][1]; - sbuf[2] = v[j][2]; - if (rmass) sbuf[3] = rmass[j]; - else sbuf[3] = mass[type[j]]; - MPI_Sendrecv(sbuf,4,MPI_DOUBLE,all[0].proc,0, - rbuf,4,MPI_DOUBLE,all[0].proc,0,world,&status); - vcm[0] = (sbuf[3]*sbuf[0] + rbuf[3]*rbuf[0]) / (sbuf[3] + rbuf[3]); - vcm[1] = (sbuf[3]*sbuf[1] + rbuf[3]*rbuf[1]) / (sbuf[3] + rbuf[3]); - vcm[2] = (sbuf[3]*sbuf[2] + rbuf[3]*rbuf[2]) / (sbuf[3] + rbuf[3]); - v[j][0] = 2.0 * vcm[0] - sbuf[0]; - v[j][1] = 2.0 * vcm[1] - sbuf[1]; - v[j][2] = 2.0 * vcm[2] - sbuf[2]; - eswap += sbuf[3] * (vcm[0] * (vcm[0] - sbuf[0]) + - vcm[1] * (vcm[1] - sbuf[1]) + - vcm[2] * (vcm[2] - sbuf[2])); - } - } - - // tally energy exchange from all swaps - - double eswap_all; - MPI_Allreduce(&eswap,&eswap_all,1,MPI_DOUBLE,MPI_SUM,world); - e_exchange += force->mvv2e * eswap_all; -} - -/* ---------------------------------------------------------------------- */ - -double FixThermalConductivity::compute_scalar() -{ - return e_exchange; -} diff --git a/src/fix_thermal_conductivity.h b/src/fix_thermal_conductivity.h deleted file mode 100644 index 22688bb7c0..0000000000 --- a/src/fix_thermal_conductivity.h +++ /dev/null @@ -1,75 +0,0 @@ -/* ---------------------------------------------------------------------- - 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 FIX_CLASS - -FixStyle(thermal/conductivity,FixThermalConductivity) - -#else - -#ifndef LMP_FIX_THERMAL_CONDUCTIVITY_H -#define LMP_FIX_THERMAL_CONDUCTIVITY_H - -#include "fix.h" - -namespace LAMMPS_NS { - -class FixThermalConductivity : public Fix { - public: - FixThermalConductivity(class LAMMPS *, int, char **); - ~FixThermalConductivity(); - int setmask(); - void init(); - void end_of_step(); - double compute_scalar(); - - private: - int me; - int edim,nbin,periodicity; - int nswap; - double prd,boxlo,boxhi; - double slablo_lo,slablo_hi,slabhi_lo,slabhi_hi; - double e_exchange; - - int nlo,nhi; - int *index_lo,*index_hi; - double *ke_lo,*ke_hi; -}; - -} - -#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: Fix thermal/conductivity swap value must be positive - -Self-explanatory. - -W: Fix thermal/conductivity comes before fix ave/spatial - -The order of these 2 fixes in your input script is such that fix -thermal/conductivity comes first. If you are using fix ave/spatial to -measure the temperature profile induced by fix viscosity, then this -may cause a glitch in the profile since you are averaging immediately -after swaps have occurred. Flipping the order of the 2 fixes -typically helps. - -*/ diff --git a/src/fix_ttm.cpp b/src/fix_ttm.cpp deleted file mode 100644 index abd145a092..0000000000 --- a/src/fix_ttm.cpp +++ /dev/null @@ -1,685 +0,0 @@ -/* ---------------------------------------------------------------------- - 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. -------------------------------------------------------------------------- */ - -/* ---------------------------------------------------------------------- - Contributing authors: Paul Crozier (SNL) - Carolyn Phillips (University of Michigan) -------------------------------------------------------------------------- */ - -#include "lmptype.h" -#include "mpi.h" -#include "math.h" -#include "string.h" -#include "stdlib.h" -#include "fix_ttm.h" -#include "atom.h" -#include "force.h" -#include "update.h" -#include "domain.h" -#include "region.h" -#include "respa.h" -#include "comm.h" -#include "random_mars.h" -#include "memory.h" -#include "error.h" - -using namespace LAMMPS_NS; -using namespace FixConst; - -#define MAXLINE 1024 - -/* ---------------------------------------------------------------------- */ - -FixTTM::FixTTM(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg) -{ - if (narg < 15) error->all(FLERR,"Illegal fix ttm command"); - - vector_flag = 1; - size_vector = 2; - global_freq = 1; - extvector = 1; - nevery = 1; - restart_peratom = 1; - restart_global = 1; - - seed = force->inumeric(FLERR,arg[3]); - electronic_specific_heat = force->numeric(FLERR,arg[4]); - electronic_density = force->numeric(FLERR,arg[5]); - electronic_thermal_conductivity = force->numeric(FLERR,arg[6]); - gamma_p = force->numeric(FLERR,arg[7]); - gamma_s = force->numeric(FLERR,arg[8]); - v_0 = force->numeric(FLERR,arg[9]); - nxnodes = force->inumeric(FLERR,arg[10]); - nynodes = force->inumeric(FLERR,arg[11]); - nznodes = force->inumeric(FLERR,arg[12]); - - fpr = fopen(arg[13],"r"); - if (fpr == NULL) { - char str[128]; - sprintf(str,"Cannot open file %s",arg[13]); - error->one(FLERR,str); - } - - nfileevery = force->inumeric(FLERR,arg[14]); - - if (nfileevery) { - if (narg != 16) error->all(FLERR,"Illegal fix ttm command"); - MPI_Comm_rank(world,&me); - if (me == 0) { - fp = fopen(arg[15],"w"); - if (fp == NULL) { - char str[128]; - sprintf(str,"Cannot open fix ttm file %s",arg[15]); - error->one(FLERR,str); - } - } - } - - // error check - - if (seed <= 0) error->all(FLERR,"Invalid random number seed in fix ttm command"); - if (electronic_specific_heat <= 0.0) - error->all(FLERR,"Fix ttm electronic_specific_heat must be > 0.0"); - if (electronic_density <= 0.0) - error->all(FLERR,"Fix ttm electronic_density must be > 0.0"); - if (electronic_thermal_conductivity < 0.0) - error->all(FLERR,"Fix ttm electronic_thermal_conductivity must be >= 0.0"); - if (gamma_p <= 0.0) error->all(FLERR,"Fix ttm gamma_p must be > 0.0"); - if (gamma_s < 0.0) error->all(FLERR,"Fix ttm gamma_s must be >= 0.0"); - if (v_0 < 0.0) error->all(FLERR,"Fix ttm v_0 must be >= 0.0"); - if (nxnodes <= 0 || nynodes <= 0 || nznodes <= 0) - error->all(FLERR,"Fix ttm number of nodes must be > 0"); - - v_0_sq = v_0*v_0; - - // initialize Marsaglia RNG with processor-unique seed - - random = new RanMars(lmp,seed + comm->me); - - // allocate per-type arrays for force prefactors - - gfactor1 = new double[atom->ntypes+1]; - gfactor2 = new double[atom->ntypes+1]; - - // allocate 3d grid variables - - total_nnodes = nxnodes*nynodes*nznodes; - - memory->create(nsum,nxnodes,nynodes,nznodes,"ttm:nsum"); - memory->create(nsum_all,nxnodes,nynodes,nznodes,"ttm:nsum_all"); - memory->create(T_initial_set,nxnodes,nynodes,nznodes,"ttm:T_initial_set"); - memory->create(sum_vsq,nxnodes,nynodes,nznodes,"ttm:sum_vsq"); - memory->create(sum_mass_vsq,nxnodes,nynodes,nznodes,"ttm:sum_mass_vsq"); - memory->create(sum_vsq_all,nxnodes,nynodes,nznodes,"ttm:sum_vsq_all"); - memory->create(sum_mass_vsq_all,nxnodes,nynodes,nznodes, - "ttm:sum_mass_vsq_all"); - memory->create(T_electron_old,nxnodes,nynodes,nznodes,"ttm:T_electron_old"); - memory->create(T_electron,nxnodes,nynodes,nznodes,"ttm:T_electron"); - memory->create(net_energy_transfer,nxnodes,nynodes,nznodes, - "TTM:net_energy_transfer"); - memory->create(net_energy_transfer_all,nxnodes,nynodes,nznodes, - "TTM:net_energy_transfer_all"); - - flangevin = NULL; - grow_arrays(atom->nmax); - - // zero out the flangevin array - - for (int i = 0; i < atom->nmax; i++) { - flangevin[i][0] = 0; - flangevin[i][1] = 0; - flangevin[i][2] = 0; - } - - atom->add_callback(0); - atom->add_callback(1); - - // set initial electron temperatures from user input file - - if (me == 0) read_initial_electron_temperatures(); - MPI_Bcast(&T_electron[0][0][0],total_nnodes,MPI_DOUBLE,0,world); -} - -/* ---------------------------------------------------------------------- */ - -FixTTM::~FixTTM() -{ - if (nfileevery && me == 0) fclose(fp); - - delete random; - - delete [] gfactor1; - delete [] gfactor2; - - memory->destroy(nsum); - memory->destroy(nsum_all); - memory->destroy(T_initial_set); - memory->destroy(sum_vsq); - memory->destroy(sum_mass_vsq); - memory->destroy(sum_vsq_all); - memory->destroy(sum_mass_vsq_all); - memory->destroy(T_electron_old); - memory->destroy(T_electron); - memory->destroy(flangevin); - memory->destroy(net_energy_transfer); - memory->destroy(net_energy_transfer_all); -} - -/* ---------------------------------------------------------------------- */ - -int FixTTM::setmask() -{ - int mask = 0; - mask |= POST_FORCE; - mask |= POST_FORCE_RESPA; - mask |= END_OF_STEP; - return mask; -} - -/* ---------------------------------------------------------------------- */ - -void FixTTM::init() -{ - if (domain->dimension == 2) - error->all(FLERR,"Cannot use fix ttm with 2d simulation"); - if (domain->nonperiodic != 0) - error->all(FLERR,"Cannot use nonperiodic boundares with fix ttm"); - if (domain->triclinic) - error->all(FLERR,"Cannot use fix ttm with triclinic box"); - - // set force prefactors - - for (int i = 1; i <= atom->ntypes; i++) { - gfactor1[i] = - gamma_p / force->ftm2v; - gfactor2[i] = - sqrt(24.0*force->boltz*gamma_p/update->dt/force->mvv2e) / force->ftm2v; - } - - for (int ixnode = 0; ixnode < nxnodes; ixnode++) - for (int iynode = 0; iynode < nynodes; iynode++) - for (int iznode = 0; iznode < nznodes; iznode++) - net_energy_transfer_all[ixnode][iynode][iznode] = 0; - - if (strstr(update->integrate_style,"respa")) - nlevels_respa = ((Respa *) update->integrate)->nlevels; -} - -/* ---------------------------------------------------------------------- */ - -void FixTTM::setup(int vflag) -{ - if (strstr(update->integrate_style,"verlet")) - post_force_setup(vflag); - else { - ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1); - post_force_respa_setup(vflag,nlevels_respa-1,0); - ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1); - } -} - -/* ---------------------------------------------------------------------- */ - -void FixTTM::post_force(int vflag) -{ - double **x = atom->x; - double **v = atom->v; - double **f = atom->f; - int *type = atom->type; - int *mask = atom->mask; - int nlocal = atom->nlocal; - - double gamma1,gamma2; - - // apply damping and thermostat to all atoms in fix group - - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - - double xscale = (x[i][0] - domain->boxlo[0])/domain->xprd; - double yscale = (x[i][1] - domain->boxlo[1])/domain->yprd; - double zscale = (x[i][2] - domain->boxlo[2])/domain->zprd; - int ixnode = static_cast(xscale*nxnodes); - int iynode = static_cast(yscale*nynodes); - int iznode = static_cast(zscale*nznodes); - while (ixnode > nxnodes-1) ixnode -= nxnodes; - while (iynode > nynodes-1) iynode -= nynodes; - while (iznode > nznodes-1) iznode -= nznodes; - while (ixnode < 0) ixnode += nxnodes; - while (iynode < 0) iynode += nynodes; - while (iznode < 0) iznode += nznodes; - - if (T_electron[ixnode][iynode][iznode] < 0) - error->all(FLERR,"Electronic temperature dropped below zero"); - - double tsqrt = sqrt(T_electron[ixnode][iynode][iznode]); - - gamma1 = gfactor1[type[i]]; - double vsq = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]; - if (vsq > v_0_sq) gamma1 *= (gamma_p + gamma_s)/gamma_p; - gamma2 = gfactor2[type[i]] * tsqrt; - - flangevin[i][0] = gamma1*v[i][0] + gamma2*(random->uniform()-0.5); - flangevin[i][1] = gamma1*v[i][1] + gamma2*(random->uniform()-0.5); - flangevin[i][2] = gamma1*v[i][2] + gamma2*(random->uniform()-0.5); - - f[i][0] += flangevin[i][0]; - f[i][1] += flangevin[i][1]; - f[i][2] += flangevin[i][2]; - } - } -} - -/* ---------------------------------------------------------------------- */ - -void FixTTM::post_force_setup(int vflag) -{ - double **f = atom->f; - int *mask = atom->mask; - int nlocal = atom->nlocal; - - // apply langevin forces that have been stored from previous run - - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - f[i][0] += flangevin[i][0]; - f[i][1] += flangevin[i][1]; - f[i][2] += flangevin[i][2]; - } - } -} - -/* ---------------------------------------------------------------------- */ - -void FixTTM::post_force_respa(int vflag, int ilevel, int iloop) -{ - if (ilevel == nlevels_respa-1) post_force(vflag); -} - -/* ---------------------------------------------------------------------- */ - -void FixTTM::post_force_respa_setup(int vflag, int ilevel, int iloop) -{ - if (ilevel == nlevels_respa-1) post_force_setup(vflag); -} - -/* ---------------------------------------------------------------------- */ - -void FixTTM::reset_dt() -{ - for (int i = 1; i <= atom->ntypes; i++) - gfactor2[i] = - sqrt(24.0*force->boltz*gamma_p/update->dt/force->mvv2e) / force->ftm2v; -} - -/* ---------------------------------------------------------------------- - read in initial electron temperatures from a user-specified file - only called by proc 0 -------------------------------------------------------------------------- */ - -void FixTTM::read_initial_electron_temperatures() -{ - char line[MAXLINE]; - - for (int ixnode = 0; ixnode < nxnodes; ixnode++) - for (int iynode = 0; iynode < nynodes; iynode++) - for (int iznode = 0; iznode < nznodes; iznode++) - T_initial_set[ixnode][iynode][iznode] = 0; - - // read initial electron temperature values from file - - int ixnode,iynode,iznode; - double T_tmp; - while (1) { - if (fgets(line,MAXLINE,fpr) == NULL) break; - sscanf(line,"%d %d %d %lg",&ixnode,&iynode,&iznode,&T_tmp); - if (T_tmp < 0.0) error->one(FLERR,"Fix ttm electron temperatures must be > 0.0"); - T_electron[ixnode][iynode][iznode] = T_tmp; - T_initial_set[ixnode][iynode][iznode] = 1; - } - - for (int ixnode = 0; ixnode < nxnodes; ixnode++) - for (int iynode = 0; iynode < nynodes; iynode++) - for (int iznode = 0; iznode < nznodes; iznode++) - if (T_initial_set[ixnode][iynode][iznode] == 0) - error->one(FLERR,"Initial temperatures not all set in fix ttm"); - - // close file - - fclose(fpr); -} - -/* ---------------------------------------------------------------------- */ - -void FixTTM::end_of_step() -{ - 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; - - for (int ixnode = 0; ixnode < nxnodes; ixnode++) - for (int iynode = 0; iynode < nynodes; iynode++) - for (int iznode = 0; iznode < nznodes; iznode++) - net_energy_transfer[ixnode][iynode][iznode] = 0; - - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - double xscale = (x[i][0] - domain->boxlo[0])/domain->xprd; - double yscale = (x[i][1] - domain->boxlo[1])/domain->yprd; - double zscale = (x[i][2] - domain->boxlo[2])/domain->zprd; - int ixnode = static_cast(xscale*nxnodes); - int iynode = static_cast(yscale*nynodes); - int iznode = static_cast(zscale*nznodes); - while (ixnode > nxnodes-1) ixnode -= nxnodes; - while (iynode > nynodes-1) iynode -= nynodes; - while (iznode > nznodes-1) iznode -= nznodes; - while (ixnode < 0) ixnode += nxnodes; - while (iynode < 0) iynode += nynodes; - while (iznode < 0) iznode += nznodes; - net_energy_transfer[ixnode][iynode][iznode] += - (flangevin[i][0]*v[i][0] + flangevin[i][1]*v[i][1] + - flangevin[i][2]*v[i][2]); - } - - MPI_Allreduce(&net_energy_transfer[0][0][0], - &net_energy_transfer_all[0][0][0], - total_nnodes,MPI_DOUBLE,MPI_SUM,world); - - double dx = domain->xprd/nxnodes; - double dy = domain->yprd/nynodes; - double dz = domain->zprd/nznodes; - double del_vol = dx*dy*dz; - - // num_inner_timesteps = # of inner steps (thermal solves) - // required this MD step to maintain a stable explicit solve - - int num_inner_timesteps = 1; - double inner_dt = update->dt; - double stability_criterion = 1.0 - - 2.0*inner_dt/(electronic_specific_heat*electronic_density) * - (electronic_thermal_conductivity*(1.0/dx/dx + 1.0/dy/dy + 1.0/dz/dz)); - if (stability_criterion < 0.0) { - inner_dt = 0.5*(electronic_specific_heat*electronic_density) / - (electronic_thermal_conductivity*(1.0/dx/dx + 1.0/dy/dy + 1.0/dz/dz)); - num_inner_timesteps = static_cast(update->dt/inner_dt) + 1; - inner_dt = update->dt/double(num_inner_timesteps); - if (num_inner_timesteps > 1000000) - error->warning(FLERR,"Too many inner timesteps in fix ttm",0); - } - - for (int ith_inner_timestep = 0; ith_inner_timestep < num_inner_timesteps; - ith_inner_timestep++) { - - for (int ixnode = 0; ixnode < nxnodes; ixnode++) - for (int iynode = 0; iynode < nynodes; iynode++) - for (int iznode = 0; iznode < nznodes; iznode++) - T_electron_old[ixnode][iynode][iznode] = - T_electron[ixnode][iynode][iznode]; - - // compute new electron T profile - - for (int ixnode = 0; ixnode < nxnodes; ixnode++) - for (int iynode = 0; iynode < nynodes; iynode++) - for (int iznode = 0; iznode < nznodes; iznode++) { - int right_xnode = ixnode + 1; - int right_ynode = iynode + 1; - int right_znode = iznode + 1; - if (right_xnode == nxnodes) right_xnode = 0; - if (right_ynode == nynodes) right_ynode = 0; - if (right_znode == nznodes) right_znode = 0; - int left_xnode = ixnode - 1; - int left_ynode = iynode - 1; - int left_znode = iznode - 1; - if (left_xnode == -1) left_xnode = nxnodes - 1; - if (left_ynode == -1) left_ynode = nynodes - 1; - if (left_znode == -1) left_znode = nznodes - 1; - T_electron[ixnode][iynode][iznode] = - T_electron_old[ixnode][iynode][iznode] + - inner_dt/(electronic_specific_heat*electronic_density) * - (electronic_thermal_conductivity * - ((T_electron_old[right_xnode][iynode][iznode] + - T_electron_old[left_xnode][iynode][iznode] - - 2*T_electron_old[ixnode][iynode][iznode])/dx/dx + - (T_electron_old[ixnode][right_ynode][iznode] + - T_electron_old[ixnode][left_ynode][iznode] - - 2*T_electron_old[ixnode][iynode][iznode])/dy/dy + - (T_electron_old[ixnode][iynode][right_znode] + - T_electron_old[ixnode][iynode][left_znode] - - 2*T_electron_old[ixnode][iynode][iznode])/dz/dz) - - (net_energy_transfer_all[ixnode][iynode][iznode])/del_vol); - } - } - - // output nodal temperatures for current timestep - - if ((nfileevery) && !(update->ntimestep % nfileevery)) { - - // compute atomic Ta for each grid point - - for (int ixnode = 0; ixnode < nxnodes; ixnode++) - for (int iynode = 0; iynode < nynodes; iynode++) - for (int iznode = 0; iznode < nznodes; iznode++) { - nsum[ixnode][iynode][iznode] = 0; - nsum_all[ixnode][iynode][iznode] = 0; - sum_vsq[ixnode][iynode][iznode] = 0.0; - sum_mass_vsq[ixnode][iynode][iznode] = 0.0; - sum_vsq_all[ixnode][iynode][iznode] = 0.0; - sum_mass_vsq_all[ixnode][iynode][iznode] = 0.0; - } - - double massone; - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - if (rmass) massone = rmass[i]; - else massone = mass[type[i]]; - double xscale = (x[i][0] - domain->boxlo[0])/domain->xprd; - double yscale = (x[i][1] - domain->boxlo[1])/domain->yprd; - double zscale = (x[i][2] - domain->boxlo[2])/domain->zprd; - int ixnode = static_cast(xscale*nxnodes); - int iynode = static_cast(yscale*nynodes); - int iznode = static_cast(zscale*nznodes); - while (ixnode > nxnodes-1) ixnode -= nxnodes; - while (iynode > nynodes-1) iynode -= nynodes; - while (iznode > nznodes-1) iznode -= nznodes; - while (ixnode < 0) ixnode += nxnodes; - while (iynode < 0) iynode += nynodes; - while (iznode < 0) iznode += nznodes; - double vsq = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]; - nsum[ixnode][iynode][iznode] += 1; - sum_vsq[ixnode][iynode][iznode] += vsq; - sum_mass_vsq[ixnode][iynode][iznode] += massone*vsq; - } - - MPI_Allreduce(&nsum[0][0][0],&nsum_all[0][0][0],total_nnodes, - MPI_INT,MPI_SUM,world); - MPI_Allreduce(&sum_vsq[0][0][0],&sum_vsq_all[0][0][0],total_nnodes, - MPI_DOUBLE,MPI_SUM,world); - MPI_Allreduce(&sum_mass_vsq[0][0][0],&sum_mass_vsq_all[0][0][0], - total_nnodes,MPI_DOUBLE,MPI_SUM,world); - - if (me == 0) { - fprintf(fp,BIGINT_FORMAT,update->ntimestep); - - double T_a; - for (int ixnode = 0; ixnode < nxnodes; ixnode++) - for (int iynode = 0; iynode < nynodes; iynode++) - for (int iznode = 0; iznode < nznodes; iznode++) { - T_a = 0; - if (nsum_all[ixnode][iynode][iznode] > 0) - T_a = sum_mass_vsq_all[ixnode][iynode][iznode]/ - (3.0*force->boltz*nsum_all[ixnode][iynode][iznode]/force->mvv2e); - fprintf(fp," %f",T_a); - } - - fprintf(fp,"\t"); - for (int ixnode = 0; ixnode < nxnodes; ixnode++) - for (int iynode = 0; iynode < nynodes; iynode++) - for (int iznode = 0; iznode < nznodes; iznode++) - fprintf(fp,"%f ",T_electron[ixnode][iynode][iznode]); - fprintf(fp,"\n"); - } - } -} - -/* ---------------------------------------------------------------------- - memory usage of 3d grid -------------------------------------------------------------------------- */ - -double FixTTM::memory_usage() -{ - double bytes = 0.0; - bytes += 5*total_nnodes * sizeof(int); - bytes += 14*total_nnodes * sizeof(double); - return bytes; -} - -/* ---------------------------------------------------------------------- */ - -void FixTTM::grow_arrays(int ngrow) -{ - - memory->grow(flangevin,ngrow,3,"TTM:flangevin"); - -} - -/* ---------------------------------------------------------------------- - return the energy of the electronic subsystem or the net_energy transfer - between the subsystems -------------------------------------------------------------------------- */ - -double FixTTM::compute_vector(int n) -{ - double e_energy = 0.0; - double transfer_energy = 0.0; - - double dx = domain->xprd/nxnodes; - double dy = domain->yprd/nynodes; - double dz = domain->zprd/nznodes; - double del_vol = dx*dy*dz; - - for (int ixnode = 0; ixnode < nxnodes; ixnode++) - for (int iynode = 0; iynode < nynodes; iynode++) - for (int iznode = 0; iznode < nznodes; iznode++) { - e_energy += - T_electron[ixnode][iynode][iznode]*electronic_specific_heat* - electronic_density*del_vol; - transfer_energy += - net_energy_transfer_all[ixnode][iynode][iznode]*update->dt; - } - - if (n == 0) return e_energy; - if (n == 1) return transfer_energy; - return 0.0; -} - -/* ---------------------------------------------------------------------- - pack entire state of Fix into one write -------------------------------------------------------------------------- */ - -void FixTTM::write_restart(FILE *fp) -{ - double *rlist; - memory->create(rlist,nxnodes*nynodes*nznodes+1,"TTM:rlist"); - - int n = 0; - rlist[n++] = seed; - - for (int ixnode = 0; ixnode < nxnodes; ixnode++) - for (int iynode = 0; iynode < nynodes; iynode++) - for (int iznode = 0; iznode < nznodes; iznode++) - rlist[n++] = T_electron[ixnode][iynode][iznode]; - - if (comm->me == 0) { - int size = n * sizeof(double); - fwrite(&size,sizeof(int),1,fp); - fwrite(rlist,sizeof(double),n,fp); - } - - memory->destroy(rlist); -} - -/* ---------------------------------------------------------------------- - use state info from restart file to restart the Fix -------------------------------------------------------------------------- */ - -void FixTTM::restart(char *buf) -{ - int n = 0; - double *rlist = (double *) buf; - - // the seed must be changed from the initial seed - - seed = static_cast (0.5*rlist[n++]); - - for (int ixnode = 0; ixnode < nxnodes; ixnode++) - for (int iynode = 0; iynode < nynodes; iynode++) - for (int iznode = 0; iznode < nznodes; iznode++) - T_electron[ixnode][iynode][iznode] = rlist[n++]; - - delete random; - random = new RanMars(lmp,seed+comm->me); -} - -/* ---------------------------------------------------------------------- - pack values in local atom-based arrays for restart file -------------------------------------------------------------------------- */ - -int FixTTM::pack_restart(int i, double *buf) -{ - buf[0] = 4; - buf[1] = flangevin[i][0]; - buf[2] = flangevin[i][1]; - buf[3] = flangevin[i][2]; - return 4; -} - -/* ---------------------------------------------------------------------- - unpack values from atom->extra array to restart the fix -------------------------------------------------------------------------- */ - -void FixTTM::unpack_restart(int nlocal, int nth) -{ - double **extra = atom->extra; - - // skip to Nth set of extra values - - int m = 0; - for (int i = 0; i < nth; i++) m += static_cast (extra[nlocal][m]); - m++; - - flangevin[nlocal][0] = extra[nlocal][m++]; - flangevin[nlocal][1] = extra[nlocal][m++]; - flangevin[nlocal][2] = extra[nlocal][m++]; -} - -/* ---------------------------------------------------------------------- - maxsize of any atom's restart data -------------------------------------------------------------------------- */ - -int FixTTM::maxsize_restart() -{ - return 4; -} - -/* ---------------------------------------------------------------------- - size of atom nlocal's restart data -------------------------------------------------------------------------- */ - -int FixTTM::size_restart(int nlocal) -{ - return 4; -} diff --git a/src/fix_ttm.h b/src/fix_ttm.h deleted file mode 100644 index f8ed2e626b..0000000000 --- a/src/fix_ttm.h +++ /dev/null @@ -1,156 +0,0 @@ -/* ---------------------------------------------------------------------- - 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 FIX_CLASS - -FixStyle(ttm,FixTTM) - -#else - -#ifndef LMP_FIX_TTM_H -#define LMP_FIX_TTM_H - -#include "fix.h" - -namespace LAMMPS_NS { - -class FixTTM : public Fix { - public: - FixTTM(class LAMMPS *, int, char **); - ~FixTTM(); - int setmask(); - void init(); - void setup(int); - void post_force(int); - void post_force_respa(int, int, int); - void post_force_setup(int); - void post_force_respa_setup(int, int, int); - void end_of_step(); - void reset_dt(); - void write_restart(FILE *); - void restart(char *); - int pack_restart(int, double *); - void unpack_restart(int, int); - int size_restart(int); - int maxsize_restart(); - double memory_usage(); - void grow_arrays(int); - double compute_vector(int); - - private: - int me; - int nfileevery; - int nlevels_respa; - int seed; - class RanMars *random; - FILE *fp,*fpr; - int nxnodes,nynodes,nznodes,total_nnodes; - int ***nsum; - int ***nsum_all,***T_initial_set; - double *gfactor1,*gfactor2,*ratio; - double **flangevin; - double ***T_electron,***T_electron_old; - double ***sum_vsq,***sum_mass_vsq; - double ***sum_vsq_all,***sum_mass_vsq_all; - double ***net_energy_transfer,***net_energy_transfer_all; - double electronic_specific_heat,electronic_density; - double electronic_thermal_conductivity; - double gamma_p,gamma_s,v_0,v_0_sq; - - void read_initial_electron_temperatures(); -}; - -} - -#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: Cannot open file %s - -The specified file cannot be opened. Check that the path and name are -correct. - -E: Cannot open fix ttm file %s - -The output file for the fix ttm command cannot be opened. Check that -the path and name are correct. - -E: Invalid random number seed in fix ttm command - -Random number seed must be > 0. - -E: Fix ttm electronic_specific_heat must be > 0.0 - -Self-explanatory. - -E: Fix ttm electronic_density must be > 0.0 - -Self-explanatory. - -E: Fix ttm electronic_thermal_conductivity must be >= 0.0 - -Self-explanatory. - -E: Fix ttm gamma_p must be > 0.0 - -Self-explanatory. - -E: Fix ttm gamma_s must be >= 0.0 - -Self-explanatory. - -E: Fix ttm v_0 must be >= 0.0 - -Self-explanatory. - -E: Fix ttm number of nodes must be > 0 - -Self-explanatory. - -E: Cannot use fix ttm with 2d simulation - -This is a current restriction of this fix due to the grid it creates. - -E: Cannot use nonperiodic boundares with fix ttm - -This fix requires a fully periodic simulation box. - -E: Cannot use fix ttm with triclinic box - -This is a current restriction of this fix due to the grid it creates. - -E: Electronic temperature dropped below zero - -Something has gone wrong with the fix ttm electron temperature model. - -E: Fix ttm electron temperatures must be > 0.0 - -Self-explanatory. - -E: Initial temperatures not all set in fix ttm - -Self-explantory. - -W: Too many inner timesteps in fix ttm - -Self-explanatory. - -*/ diff --git a/src/fix_viscosity.cpp b/src/fix_viscosity.cpp deleted file mode 100644 index ee94a6e6a3..0000000000 --- a/src/fix_viscosity.cpp +++ /dev/null @@ -1,309 +0,0 @@ -/* ---------------------------------------------------------------------- - 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. -------------------------------------------------------------------------- */ - -/* ---------------------------------------------------------------------- - Contributing author: Craig Tenney (UND) added support - for swapping atoms of different masses -------------------------------------------------------------------------- */ - -#include "math.h" -#include "mpi.h" -#include "string.h" -#include "stdlib.h" -#include "fix_viscosity.h" -#include "atom.h" -#include "domain.h" -#include "modify.h" -#include "error.h" -#include "force.h" - -using namespace LAMMPS_NS; -using namespace FixConst; - -// needs to be big, but not so big that lose precision when subtract velocity - -#define BIG 1.0e10 - -/* ---------------------------------------------------------------------- */ - -FixViscosity::FixViscosity(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg) -{ - if (narg < 7) error->all(FLERR,"Illegal fix viscosity command"); - - MPI_Comm_rank(world,&me); - - nevery = force->inumeric(FLERR,arg[3]); - if (nevery <= 0) error->all(FLERR,"Illegal fix viscosity command"); - - scalar_flag = 1; - global_freq = nevery; - extscalar = 0; - - if (strcmp(arg[4],"x") == 0) vdim = 0; - else if (strcmp(arg[4],"y") == 0) vdim = 1; - else if (strcmp(arg[4],"z") == 0) vdim = 2; - else error->all(FLERR,"Illegal fix viscosity command"); - - if (strcmp(arg[5],"x") == 0) pdim = 0; - else if (strcmp(arg[5],"y") == 0) pdim = 1; - else if (strcmp(arg[5],"z") == 0) pdim = 2; - else error->all(FLERR,"Illegal fix viscosity command"); - - nbin = force->inumeric(FLERR,arg[6]); - if (nbin % 2 || nbin <= 2) error->all(FLERR,"Illegal fix viscosity command"); - - // optional keywords - - nswap = 1; - vtarget = BIG; - - int iarg = 7; - while (iarg < narg) { - if (strcmp(arg[iarg],"swap") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix viscosity command"); - nswap = force->inumeric(FLERR,arg[iarg+1]); - if (nswap <= 0) - error->all(FLERR,"Fix viscosity swap value must be positive"); - iarg += 2; - } else if (strcmp(arg[iarg],"vtarget") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix viscosity command"); - if (strcmp(arg[iarg+1],"INF") == 0) vtarget = BIG; - else vtarget = force->numeric(FLERR,arg[iarg+1]); - if (vtarget <= 0.0) - error->all(FLERR,"Fix viscosity vtarget value must be positive"); - iarg += 2; - } else error->all(FLERR,"Illegal fix viscosity command"); - } - - // initialize array sizes to nswap+1 so have space to shift values down - - pos_index = new int[nswap+1]; - neg_index = new int[nswap+1]; - pos_delta = new double[nswap+1]; - neg_delta = new double[nswap+1]; - - p_exchange = 0.0; -} - -/* ---------------------------------------------------------------------- */ - -FixViscosity::~FixViscosity() -{ - delete [] pos_index; - delete [] neg_index; - delete [] pos_delta; - delete [] neg_delta; -} - -/* ---------------------------------------------------------------------- */ - -int FixViscosity::setmask() -{ - int mask = 0; - mask |= END_OF_STEP; - return mask; -} - -/* ---------------------------------------------------------------------- */ - -void FixViscosity::init() -{ - // warn if any fix ave/spatial comes after this fix - // can cause glitch in averaging since ave will happen after swap - - int foundme = 0; - for (int i = 0; i < modify->nfix; i++) { - if (modify->fix[i] == this) foundme = 1; - if (foundme && strcmp(modify->fix[i]->style,"ave/spatial") == 0 && me == 0) - error->warning(FLERR,"Fix viscosity comes before fix ave/spatial"); - } - - // set bounds of 2 slabs in pdim - // only necessary for static box, else re-computed in end_of_step() - // lo bin is always bottom bin - // hi bin is just above half height - - if (domain->box_change == 0) { - prd = domain->prd[pdim]; - boxlo = domain->boxlo[pdim]; - boxhi = domain->boxhi[pdim]; - double binsize = (boxhi-boxlo) / nbin; - slablo_lo = boxlo; - slablo_hi = boxlo + binsize; - slabhi_lo = boxlo + (nbin/2)*binsize; - slabhi_hi = boxlo + (nbin/2+1)*binsize; - } - - periodicity = domain->periodicity[pdim]; -} - -/* ---------------------------------------------------------------------- */ - -void FixViscosity::end_of_step() -{ - int i,m,insert; - double coord,delta; - MPI_Status status; - struct { - double value; - int proc; - } mine[2],all[2]; - - // if box changes, recompute bounds of 2 slabs in pdim - - if (domain->box_change) { - prd = domain->prd[pdim]; - boxlo = domain->boxlo[pdim]; - boxhi = domain->boxhi[pdim]; - double binsize = (boxhi-boxlo) / nbin; - slablo_lo = boxlo; - slablo_hi = boxlo + binsize; - slabhi_lo = boxlo + (nbin/2)*binsize; - slabhi_hi = boxlo + (nbin/2+1)*binsize; - } - - // make 2 lists of up to nswap atoms with velocity closest to +/- vtarget - // lists are sorted by closeness to vtarget - // only consider atoms in the bottom/middle slabs - // map atoms back into periodic box if necessary - // insert = location in list to insert new atom - - double **x = atom->x; - double **v = atom->v; - int *type = atom->type; - int *mask = atom->mask; - int nlocal = atom->nlocal; - - npositive = nnegative = 0; - - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - coord = x[i][pdim]; - if (coord < boxlo && periodicity) coord += prd; - else if (coord >= boxhi && periodicity) coord -= prd; - - if (coord >= slablo_lo && coord < slablo_hi) { - if (v[i][vdim] < 0.0) continue; - delta = fabs(v[i][vdim] - vtarget); - if (npositive < nswap || delta < pos_delta[nswap-1]) { - for (insert = npositive-1; insert >= 0; insert--) - if (delta > pos_delta[insert]) break; - insert++; - for (m = npositive-1; m >= insert; m--) { - pos_delta[m+1] = pos_delta[m]; - pos_index[m+1] = pos_index[m]; - } - pos_delta[insert] = delta; - pos_index[insert] = i; - if (npositive < nswap) npositive++; - } - } - - if (coord >= slabhi_lo && coord < slabhi_hi) { - if (v[i][vdim] > 0.0) continue; - delta = fabs(v[i][vdim] + vtarget); - if (nnegative < nswap || delta < neg_delta[nswap-1]) { - for (insert = nnegative-1; insert >= 0; insert--) - if (delta > neg_delta[insert]) break; - insert++; - for (m = nnegative-1; m >= insert; m--) { - neg_delta[m+1] = neg_delta[m]; - neg_index[m+1] = neg_index[m]; - } - neg_delta[insert] = delta; - neg_index[insert] = i; - if (nnegative < nswap) nnegative++; - } - } - } - - // loop over nswap pairs - // find 2 global atoms with smallest delta in bottom/top slabs - // BIG values are for procs with no atom to contribute - // MINLOC also communicates which procs own them - // exchange momenta between the 2 particles - // if I own both particles just swap, else point2point comm of vel,mass - - double *mass = atom->mass; - double *rmass = atom->rmass; - - int ipos,ineg; - double sbuf[2],rbuf[2],vcm; - - double pswap = 0.0; - mine[0].proc = mine[1].proc = me; - int ipositive = 0; - int inegative = 0; - - for (m = 0; m < nswap; m++) { - if (ipositive < npositive) mine[0].value = pos_delta[ipositive]; - else mine[0].value = BIG; - if (inegative < nnegative) mine[1].value = neg_delta[inegative]; - else mine[1].value = BIG; - - MPI_Allreduce(mine,all,2,MPI_DOUBLE_INT,MPI_MINLOC,world); - - if (all[0].value == BIG || all[1].value == BIG) continue; - - if (me == all[0].proc && me == all[1].proc) { - ipos = pos_index[ipositive++]; - ineg = neg_index[inegative++]; - rbuf[0] = v[ipos][vdim]; - if (rmass) rbuf[1] = rmass[ipos]; - else rbuf[1] = mass[type[ipos]]; - sbuf[0] = v[ineg][vdim]; - if (rmass) sbuf[1] = rmass[ineg]; - else sbuf[1] = mass[type[ineg]]; - vcm = (sbuf[1]*sbuf[0] + rbuf[1]*rbuf[0]) / (sbuf[1] + rbuf[1]); - v[ineg][vdim] = 2.0 * vcm - sbuf[0]; - v[ipos][vdim] = 2.0 * vcm - rbuf[0]; - pswap += rbuf[1] * (vcm - rbuf[0]) - sbuf[1] * (vcm - sbuf[0]); - - } else if (me == all[0].proc) { - ipos = pos_index[ipositive++]; - sbuf[0] = v[ipos][vdim]; - if (rmass) sbuf[1] = rmass[ipos]; - else sbuf[1] = mass[type[ipos]]; - MPI_Sendrecv(sbuf,2,MPI_DOUBLE,all[1].proc,0, - rbuf,2,MPI_DOUBLE,all[1].proc,0,world,&status); - vcm = (sbuf[1]*sbuf[0] + rbuf[1]*rbuf[0]) / (sbuf[1] + rbuf[1]); - v[ipos][vdim] = 2.0 * vcm - sbuf[0]; - pswap += sbuf[1] * (vcm - sbuf[0]); - - } else if (me == all[1].proc) { - ineg = neg_index[inegative++]; - sbuf[0] = v[ineg][vdim]; - if (rmass) sbuf[1] = rmass[ineg]; - else sbuf[1] = mass[type[ineg]]; - MPI_Sendrecv(sbuf,2,MPI_DOUBLE,all[0].proc,0, - rbuf,2,MPI_DOUBLE,all[0].proc,0,world,&status); - vcm = (sbuf[1]*sbuf[0] + rbuf[1]*rbuf[0]) / (sbuf[1] + rbuf[1]); - v[ineg][vdim] = 2.0 * vcm - sbuf[0]; - pswap -= sbuf[1] * (vcm - sbuf[0]); - } - } - - // tally momentum exchange from all swaps - - double pswap_all; - MPI_Allreduce(&pswap,&pswap_all,1,MPI_DOUBLE,MPI_SUM,world); - p_exchange += pswap_all; -} - -/* ---------------------------------------------------------------------- */ - -double FixViscosity::compute_scalar() -{ - return p_exchange; -} diff --git a/src/fix_viscosity.h b/src/fix_viscosity.h deleted file mode 100644 index 6f694d4f94..0000000000 --- a/src/fix_viscosity.h +++ /dev/null @@ -1,80 +0,0 @@ -/* ---------------------------------------------------------------------- - 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 FIX_CLASS - -FixStyle(viscosity,FixViscosity) - -#else - -#ifndef LMP_FIX_VISCOSITY_H -#define LMP_FIX_VISCOSITY_H - -#include "fix.h" - -namespace LAMMPS_NS { - -class FixViscosity : public Fix { - public: - FixViscosity(class LAMMPS *, int, char **); - ~FixViscosity(); - int setmask(); - void init(); - void end_of_step(); - double compute_scalar(); - - private: - int me; - int vdim,pdim,nbin,periodicity; - int nswap; - double vtarget; - double prd,boxlo,boxhi; - double slablo_lo,slablo_hi,slabhi_lo,slabhi_hi; - double p_exchange; - - int npositive,nnegative; - int *pos_index,*neg_index; - double *pos_delta,*neg_delta; -}; - -} - -#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: Fix viscosity swap value must be positive - -Self-explanatory. - -E: Fix viscosity vtarget value must be positive - -Self-explanatory. - -W: Fix viscosity comes before fix ave/spatial - -The order of these 2 fixes in your input script is such that -fix viscosity comes first. If you are using fix ave/spatial -to measure the velocity profile induced by fix viscosity, then -this may cause a glitch in the profile since you are averaging -immediately after swaps have occurred. Flipping the order -of the 2 fixes typically helps. - -*/