From 28a774a70b3acf9062a0d620380655a57fd62948 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 25 Jul 2013 22:30:22 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@10367 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/MISC/fix_deposit.cpp | 508 +++++++++++++++ src/MISC/fix_deposit.h | 98 +++ src/MISC/fix_efield.cpp | 406 ++++++++++++ src/MISC/fix_efield.h | 101 +++ src/MISC/fix_evaporate.cpp | 386 ++++++++++++ src/{compute_ti.h => MISC/fix_evaporate.h} | 61 +- src/MISC/fix_orient_fcc.cpp | 593 ++++++++++++++++++ src/MISC/fix_orient_fcc.h | 115 ++++ src/MISC/fix_thermal_conductivity.cpp | 330 ++++++++++ src/MISC/fix_thermal_conductivity.h | 75 +++ src/MISC/fix_ttm.cpp | 685 +++++++++++++++++++++ src/MISC/fix_ttm.h | 156 +++++ src/MISC/fix_viscosity.cpp | 309 ++++++++++ src/MISC/fix_viscosity.h | 80 +++ src/compute_msd_nongauss.cpp | 106 ---- src/compute_msd_nongauss.h | 51 -- src/compute_ti.cpp | 241 -------- src/fix_adapt.cpp | 389 ------------ src/fix_adapt.h | 107 ---- 19 files changed, 3871 insertions(+), 926 deletions(-) create mode 100644 src/MISC/fix_deposit.cpp create mode 100644 src/MISC/fix_deposit.h create mode 100644 src/MISC/fix_efield.cpp create mode 100644 src/MISC/fix_efield.h create mode 100644 src/MISC/fix_evaporate.cpp rename src/{compute_ti.h => MISC/fix_evaporate.h} (53%) create mode 100644 src/MISC/fix_orient_fcc.cpp create mode 100644 src/MISC/fix_orient_fcc.h create mode 100644 src/MISC/fix_thermal_conductivity.cpp create mode 100644 src/MISC/fix_thermal_conductivity.h create mode 100644 src/MISC/fix_ttm.cpp create mode 100644 src/MISC/fix_ttm.h create mode 100644 src/MISC/fix_viscosity.cpp create mode 100644 src/MISC/fix_viscosity.h delete mode 100644 src/compute_msd_nongauss.cpp delete mode 100644 src/compute_msd_nongauss.h delete mode 100644 src/compute_ti.cpp delete mode 100644 src/fix_adapt.cpp delete mode 100644 src/fix_adapt.h diff --git a/src/MISC/fix_deposit.cpp b/src/MISC/fix_deposit.cpp new file mode 100644 index 0000000000..9ee902e0c2 --- /dev/null +++ b/src/MISC/fix_deposit.cpp @@ -0,0 +1,508 @@ +/* ---------------------------------------------------------------------- + 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/MISC/fix_deposit.h b/src/MISC/fix_deposit.h new file mode 100644 index 0000000000..8b2e9f1052 --- /dev/null +++ b/src/MISC/fix_deposit.h @@ -0,0 +1,98 @@ +/* -*- 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/MISC/fix_efield.cpp b/src/MISC/fix_efield.cpp new file mode 100644 index 0000000000..f328a3f294 --- /dev/null +++ b/src/MISC/fix_efield.cpp @@ -0,0 +1,406 @@ +/* ---------------------------------------------------------------------- + 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/MISC/fix_efield.h b/src/MISC/fix_efield.h new file mode 100644 index 0000000000..48bd75174c --- /dev/null +++ b/src/MISC/fix_efield.h @@ -0,0 +1,101 @@ +/* ---------------------------------------------------------------------- + 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/MISC/fix_evaporate.cpp b/src/MISC/fix_evaporate.cpp new file mode 100644 index 0000000000..5fe5c1412f --- /dev/null +++ b/src/MISC/fix_evaporate.cpp @@ -0,0 +1,386 @@ +/* ---------------------------------------------------------------------- + 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/compute_ti.h b/src/MISC/fix_evaporate.h similarity index 53% rename from src/compute_ti.h rename to src/MISC/fix_evaporate.h index 23d3e3a4e9..f59558103a 100644 --- a/src/compute_ti.h +++ b/src/MISC/fix_evaporate.h @@ -11,34 +11,39 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#ifdef COMPUTE_CLASS +#ifdef FIX_CLASS -ComputeStyle(ti,ComputeTI) +FixStyle(evaporate,FixEvaporate) #else -#ifndef COMPUTE_TI_H -#define COMPUTE_TI_H +#ifndef LMP_FIX_EVAPORATE_H +#define LMP_FIX_EVAPORATE_H -#include "compute.h" +#include "fix.h" namespace LAMMPS_NS { -class ComputeTI : public Compute { +class FixEvaporate : public Fix { public: - ComputeTI(class LAMMPS *, int, char **); - ~ComputeTI(); + FixEvaporate(class LAMMPS *, int, char **); + ~FixEvaporate(); + int setmask(); void init(); + void pre_exchange(); double compute_scalar(); + double memory_usage(); private: - int nterms; - int *which; - int *ivar1,*ivar2; - int *ilo, *ihi; - char **var1,**var2; - class Pair **pptr; - char **pstyle; + int nevery,nflux,iregion; + int molflag; + int ndeleted; + char *idregion; + + int nmax; + int *list,*mark; + + class RanPark *random; }; } @@ -54,30 +59,22 @@ 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: Variable name for compute ti does not exist +E: Region ID for fix evaporate does not exist Self-explanatory. -E: Variable for compute ti is invalid style +E: Cannot evaporate atoms in atom_modify first group -Self-explanatory. +This is a restriction due to the way atoms are organized in +a list to enable the atom_modify first command. -E: Compute ti pair style does not exist +W: Fix evaporate may delete atom with non-zero molecule ID -Self-explanatory. +This is probably an error, since you should not delete only one atom +of a molecule. -E: Compute ti tail when pair style does not compute tail corrections +E: Fix evaporate molecule requires atom attribute molecule -Self-explanatory. - -E: Compute ti kspace style does not exist - -Self-explanatory. - -E: Energy was not tallied on needed timestep - -You are using a thermo keyword that requires potentials to -have tallied energy, but they didn't on this timestep. See the -variable doc page for ideas on how to make this work. +The atom style being used does not define a molecule ID. */ diff --git a/src/MISC/fix_orient_fcc.cpp b/src/MISC/fix_orient_fcc.cpp new file mode 100644 index 0000000000..3611478308 --- /dev/null +++ b/src/MISC/fix_orient_fcc.cpp @@ -0,0 +1,593 @@ +/* ---------------------------------------------------------------------- + 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: Koenraad Janssens and David Olmsted (SNL) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "string.h" +#include "stdlib.h" +#include "mpi.h" +#include "fix_orient_fcc.h" +#include "atom.h" +#include "update.h" +#include "respa.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "comm.h" +#include "output.h" +#include "math_const.h" +#include "memory.h" +#include "error.h" +#include "force.h" + +using namespace LAMMPS_NS; +using namespace FixConst; +using namespace MathConst; + +#define BIG 1000000000 + +/* ---------------------------------------------------------------------- */ + +FixOrientFCC::FixOrientFCC(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg) +{ + MPI_Comm_rank(world,&me); + + if (narg != 11) error->all(FLERR,"Illegal fix orient/fcc command"); + + scalar_flag = 1; + global_freq = 1; + extscalar = 1; + + peratom_flag = 1; + size_peratom_cols = 2; + peratom_freq = 1; + + nstats = force->inumeric(FLERR,arg[3]); + direction_of_motion = force->inumeric(FLERR,arg[4]); + a = force->numeric(FLERR,arg[5]); + Vxi = force->numeric(FLERR,arg[6]); + uxif_low = force->numeric(FLERR,arg[7]); + uxif_high = force->numeric(FLERR,arg[8]); + + if (direction_of_motion == 0) { + int n = strlen(arg[9]) + 1; + chifilename = new char[n]; + strcpy(chifilename,arg[9]); + n = strlen(arg[10]) + 1; + xifilename = new char[n]; + strcpy(xifilename,arg[10]); + } else if (direction_of_motion == 1) { + int n = strlen(arg[9]) + 1; + xifilename = new char[n]; + strcpy(xifilename,arg[9]); + n = strlen(arg[10]) + 1; + chifilename = new char[n]; + strcpy(chifilename,arg[10]); + } else error->all(FLERR,"Illegal fix orient/fcc command"); + + // initializations + + half_fcc_nn = 6; + use_xismooth = false; + double xicutoff = 1.57; + xicutoffsq = xicutoff * xicutoff; + cutsq = 0.5 * a*a*xicutoffsq; + nmax = 0; + + // read xi and chi reference orientations from files + + if (me == 0) { + char line[IMGMAX]; + char *result; + int count; + + FILE *infile = fopen(xifilename,"r"); + if (infile == NULL) error->one(FLERR,"Fix orient/fcc file open failed"); + for (int i = 0; i < 6; i++) { + result = fgets(line,IMGMAX,infile); + if (!result) error->one(FLERR,"Fix orient/fcc file read failed"); + count = sscanf(line,"%lg %lg %lg",&Rxi[i][0],&Rxi[i][1],&Rxi[i][2]); + if (count != 3) error->one(FLERR,"Fix orient/fcc file read failed"); + } + fclose(infile); + + infile = fopen(chifilename,"r"); + if (infile == NULL) error->one(FLERR,"Fix orient/fcc file open failed"); + for (int i = 0; i < 6; i++) { + result = fgets(line,IMGMAX,infile); + if (!result) error->one(FLERR,"Fix orient/fcc file read failed"); + count = sscanf(line,"%lg %lg %lg",&Rchi[i][0],&Rchi[i][1],&Rchi[i][2]); + if (count != 3) error->one(FLERR,"Fix orient/fcc file read failed"); + } + fclose(infile); + } + + MPI_Bcast(&Rxi[0][0],18,MPI_DOUBLE,0,world); + MPI_Bcast(&Rchi[0][0],18,MPI_DOUBLE,0,world); + + // make copy of the reference vectors + + for (int i = 0; i < 6; i++) + for (int j = 0; j < 3; j++) { + half_xi_chi_vec[0][i][j] = Rxi[i][j]; + half_xi_chi_vec[1][i][j] = Rchi[i][j]; + } + + // compute xiid,xi0,xi1 for all 12 neighbors + // xi is the favored crystal + // want order parameter when actual is Rchi + + double xi_sq,dxi[3],rchi[3]; + + xiid = 0.0; + for (int i = 0; i < 6; i++) { + rchi[0] = Rchi[i][0]; + rchi[1] = Rchi[i][1]; + rchi[2] = Rchi[i][2]; + find_best_ref(rchi,0,xi_sq,dxi); + xiid += sqrt(xi_sq); + for (int j = 0; j < 3; j++) rchi[j] = -rchi[j]; + find_best_ref(rchi,0,xi_sq,dxi); + xiid += sqrt(xi_sq); + } + + xiid /= 12.0; + xi0 = uxif_low * xiid; + xi1 = uxif_high * xiid; + + // set comm size needed by this Fix + // NOTE: doesn't seem that use_xismooth is ever true + + if (use_xismooth) comm_forward = 62; + else comm_forward = 50; + + added_energy = 0.0; + + nmax = atom->nmax; + nbr = (Nbr *) memory->smalloc(nmax*sizeof(Nbr),"orient/fcc:nbr"); + memory->create(order,nmax,2,"orient/fcc:order"); + array_atom = order; + + // zero the array since a variable may access it before first run + + int nlocal = atom->nlocal; + for (int i = 0; i < nlocal; i++) order[i][0] = order[i][1] = 0.0; +} + +/* ---------------------------------------------------------------------- */ + +FixOrientFCC::~FixOrientFCC() +{ + delete [] xifilename; + delete [] chifilename; + memory->sfree(nbr); + memory->destroy(order); +} + +/* ---------------------------------------------------------------------- */ + +int FixOrientFCC::setmask() +{ + int mask = 0; + mask |= POST_FORCE; + mask |= THERMO_ENERGY; + mask |= POST_FORCE_RESPA; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixOrientFCC::init() +{ + if (strstr(update->integrate_style,"respa")) + nlevels_respa = ((Respa *) update->integrate)->nlevels; + + // need a full neighbor list, built whenever re-neighboring occurs + + int irequest = neighbor->request((void *) this); + neighbor->requests[irequest]->pair = 0; + neighbor->requests[irequest]->fix = 1; + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->full = 1; +} + +/* ---------------------------------------------------------------------- */ + +void FixOrientFCC::init_list(int id, NeighList *ptr) +{ + list = ptr; +} + +/* ---------------------------------------------------------------------- */ + +void FixOrientFCC::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 FixOrientFCC::post_force(int vflag) +{ + int i,j,k,ii,jj,inum,jnum,m,n,nn,nsort,id_self; + int *ilist,*jlist,*numneigh,**firstneigh; + double edelta,omega; + double dx,dy,dz,rsq,xismooth,xi_sq,duxi,duxi_other; + double dxi[3]; + double *dxiptr; + bool found_myself; + + // set local ptrs + + double **x = atom->x; + double **f = atom->f; + int *mask = atom->mask; + int *tag = atom->tag; + int nlocal = atom->nlocal; + int nall = atom->nlocal + atom->nghost; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // insure nbr and order data structures are adequate size + + if (nall > nmax) { + nmax = nall; + memory->destroy(nbr); + memory->destroy(order); + nbr = (Nbr *) memory->smalloc(nmax*sizeof(Nbr),"orient/fcc:nbr"); + memory->create(order,nmax,2,"orient/fcc:order"); + array_atom = order; + } + + // loop over owned atoms and build Nbr data structure of neighbors + // use full neighbor list + + added_energy = 0.0; + int count = 0; + int mincount = BIG; + int maxcount = 0; + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + if (jnum < mincount) mincount = jnum; + if (jnum > maxcount) { + if (maxcount) delete [] sort; + sort = new Sort[jnum]; + maxcount = jnum; + } + + // loop over all neighbors of atom i + // for those within cutsq, build sort data structure + // store local id, rsq, delta vector, xismooth (if included) + + nsort = 0; + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + count++; + + dx = x[i][0] - x[j][0]; + dy = x[i][1] - x[j][1]; + dz = x[i][2] - x[j][2]; + rsq = dx*dx + dy*dy + dz*dz; + + if (rsq < cutsq) { + sort[nsort].id = j; + sort[nsort].rsq = rsq; + sort[nsort].delta[0] = dx; + sort[nsort].delta[1] = dy; + sort[nsort].delta[2] = dz; + if (use_xismooth) { + xismooth = (xicutoffsq - 2.0*rsq/(a*a)) / (xicutoffsq - 1.0); + sort[nsort].xismooth = 1.0 - fabs(1.0-xismooth); + } + nsort++; + } + } + + // sort neighbors by rsq distance + // no need to sort if nsort <= 12 + + if (nsort > 12) qsort(sort,nsort,sizeof(Sort),compare); + + // copy up to 12 nearest neighbors into nbr data structure + // operate on delta vector via find_best_ref() to compute dxi + + n = MIN(12,nsort); + nbr[i].n = n; + if (n == 0) continue; + + double xi_total = 0.0; + for (j = 0; j < n; j++) { + find_best_ref(sort[j].delta,0,xi_sq,dxi); + xi_total += sqrt(xi_sq); + nbr[i].id[j] = sort[j].id; + nbr[i].dxi[j][0] = dxi[0]/n; + nbr[i].dxi[j][1] = dxi[1]/n; + nbr[i].dxi[j][2] = dxi[2]/n; + if (use_xismooth) nbr[i].xismooth[j] = sort[j].xismooth; + } + xi_total /= n; + order[i][0] = xi_total; + + // compute potential derivative to xi + + if (xi_total < xi0) { + nbr[i].duxi = 0.0; + edelta = 0.0; + order[i][1] = 0.0; + } else if (xi_total > xi1) { + nbr[i].duxi = 0.0; + edelta = Vxi; + order[i][1] = 1.0; + } else { + omega = MY_PI2*(xi_total-xi0) / (xi1-xi0); + nbr[i].duxi = MY_PI*Vxi*sin(2.0*omega) / (2.0*(xi1-xi0)); + edelta = Vxi*(1 - cos(2.0*omega)) / 2.0; + order[i][1] = omega / MY_PI2; + } + added_energy += edelta; + } + + if (maxcount) delete [] sort; + + // communicate to acquire nbr data for ghost atoms + + comm->forward_comm_fix(this); + + // compute grain boundary force on each owned atom + // skip atoms not in group + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + if (!(mask[i] & groupbit)) continue; + n = nbr[i].n; + duxi = nbr[i].duxi; + + for (j = 0; j < n; j++) { + dxiptr = &nbr[i].dxi[j][0]; + if (use_xismooth) { + xismooth = nbr[i].xismooth[j]; + f[i][0] += duxi * dxiptr[0] * xismooth; + f[i][1] += duxi * dxiptr[1] * xismooth; + f[i][2] += duxi * dxiptr[2] * xismooth; + } else { + f[i][0] += duxi * dxiptr[0]; + f[i][1] += duxi * dxiptr[1]; + f[i][2] += duxi * dxiptr[2]; + } + + // m = local index of neighbor + // id_self = ID for atom I in atom M's neighbor list + // if M is local atom, id_self will be local ID of atom I + // if M is ghost atom, id_self will be global ID of atom I + + m = nbr[i].id[j]; + if (m < nlocal) id_self = i; + else id_self = tag[i]; + found_myself = false; + nn = nbr[m].n; + + for (k = 0; k < nn; k++) { + if (id_self == nbr[m].id[k]) { + if (found_myself) error->one(FLERR,"Fix orient/fcc found self twice"); + found_myself = true; + duxi_other = nbr[m].duxi; + dxiptr = &nbr[m].dxi[k][0]; + if (use_xismooth) { + xismooth = nbr[m].xismooth[k]; + f[i][0] -= duxi_other * dxiptr[0] * xismooth; + f[i][1] -= duxi_other * dxiptr[1] * xismooth; + f[i][2] -= duxi_other * dxiptr[2] * xismooth; + } else { + f[i][0] -= duxi_other * dxiptr[0]; + f[i][1] -= duxi_other * dxiptr[1]; + f[i][2] -= duxi_other * dxiptr[2]; + } + } + } + } + } + + // print statistics every nstats timesteps + + if (nstats && update->ntimestep % nstats == 0) { + int total; + MPI_Allreduce(&count,&total,1,MPI_INT,MPI_SUM,world); + double ave = total/atom->natoms; + + int min,max; + MPI_Allreduce(&mincount,&min,1,MPI_INT,MPI_MIN,world); + MPI_Allreduce(&maxcount,&max,1,MPI_INT,MPI_MAX,world); + + if (me == 0) { + if (screen) fprintf(screen, + "orient step " BIGINT_FORMAT ": " BIGINT_FORMAT + " atoms have %d neighbors\n", + update->ntimestep,atom->natoms,total); + if (logfile) fprintf(logfile, + "orient step " BIGINT_FORMAT ": " BIGINT_FORMAT + " atoms have %d neighbors\n", + update->ntimestep,atom->natoms,total); + if (screen) + fprintf(screen," neighs: min = %d, max = %d, ave = %g\n", + min,max,ave); + if (logfile) + fprintf(logfile," neighs: min = %d, max = %d, ave = %g\n", + min,max,ave); + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixOrientFCC::post_force_respa(int vflag, int ilevel, int iloop) +{ + if (ilevel == nlevels_respa-1) post_force(vflag); +} + +/* ---------------------------------------------------------------------- */ + +double FixOrientFCC::compute_scalar() +{ + double added_energy_total; + MPI_Allreduce(&added_energy,&added_energy_total,1,MPI_DOUBLE,MPI_SUM,world); + return added_energy_total; +} + +/* ---------------------------------------------------------------------- */ + +int FixOrientFCC::pack_comm(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,k,id,num; + int *tag = atom->tag; + int nlocal = atom->nlocal; + int m = 0; + + for (i = 0; i < n; i++) { + k = list[i]; + num = nbr[k].n; + buf[m++] = num; + buf[m++] = nbr[k].duxi; + + for (j = 0; j < num; j++) { + if (use_xismooth) buf[m++] = nbr[m].xismooth[j]; + buf[m++] = nbr[k].dxi[j][0]; + buf[m++] = nbr[k].dxi[j][1]; + buf[m++] = nbr[k].dxi[j][2]; + + // id stored in buf needs to be global ID + // if k is a local atom, it stores local IDs, so convert to global + // if k is a ghost atom (already comm'd), its IDs are already global + + id = nbr[k].id[j]; + if (k < nlocal) id = tag[id]; + buf[m++] = id; + } + + m += (12-num) * 3; + if (use_xismooth) m += 12-num; + } + + if (use_xismooth) return 62; + return 50; +} + +/* ---------------------------------------------------------------------- */ + +void FixOrientFCC::unpack_comm(int n, int first, double *buf) +{ + int i,j,num; + int last = first + n; + int m = 0; + + for (i = first; i < last; i++) { + nbr[i].n = num = static_cast (buf[m++]); + nbr[i].duxi = buf[m++]; + + for (j = 0; j < num; j++) { + if (use_xismooth) nbr[i].xismooth[j] = buf[m++]; + nbr[i].dxi[j][0] = buf[m++]; + nbr[i].dxi[j][1] = buf[m++]; + nbr[i].dxi[j][2] = buf[m++]; + nbr[i].id[j] = static_cast (buf[m++]); + } + + m += (12-num) * 3; + if (use_xismooth) m += 12-num; + } +} + +/* ---------------------------------------------------------------------- */ + +void FixOrientFCC::find_best_ref(double *displs, int which_crystal, + double &xi_sq, double *dxi) +{ + int i; + double dot,tmp; + + double best_dot = -1.0; // best is biggest (smallest angle) + int best_i = -1; + int best_sign = 0; + + for (i = 0; i < half_fcc_nn; i++) { + dot = displs[0] * half_xi_chi_vec[which_crystal][i][0] + + displs[1] * half_xi_chi_vec[which_crystal][i][1] + + displs[2] * half_xi_chi_vec[which_crystal][i][2]; + if (fabs(dot) > best_dot) { + best_dot = fabs(dot); + best_i = i; + if (dot < 0.0) best_sign = -1; + else best_sign = 1; + } + } + + xi_sq = 0.0; + for (i = 0; i < 3; i++) { + tmp = displs[i] - best_sign * half_xi_chi_vec[which_crystal][best_i][i]; + xi_sq += tmp*tmp; + } + + if (xi_sq > 0.0) { + double xi = sqrt(xi_sq); + for (i = 0; i < 3; i++) + dxi[i] = (best_sign * half_xi_chi_vec[which_crystal][best_i][i] - + displs[i]) / xi; + } else dxi[0] = dxi[1] = dxi[2] = 0.0; +} + +/* ---------------------------------------------------------------------- + compare two neighbors I and J in sort data structure + called via qsort in post_force() method + is a static method so can't access sort data structure directly + return -1 if I < J, 0 if I = J, 1 if I > J + do comparison based on rsq distance +------------------------------------------------------------------------- */ + +int FixOrientFCC::compare(const void *pi, const void *pj) +{ + FixOrientFCC::Sort *ineigh = (FixOrientFCC::Sort *) pi; + FixOrientFCC::Sort *jneigh = (FixOrientFCC::Sort *) pj; + + if (ineigh->rsq < jneigh->rsq) return -1; + else if (ineigh->rsq > jneigh->rsq) return 1; + return 0; +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based arrays +------------------------------------------------------------------------- */ + +double FixOrientFCC::memory_usage() +{ + double bytes = nmax * sizeof(Nbr); + bytes += 2*nmax * sizeof(double); + return bytes; +} diff --git a/src/MISC/fix_orient_fcc.h b/src/MISC/fix_orient_fcc.h new file mode 100644 index 0000000000..727be00539 --- /dev/null +++ b/src/MISC/fix_orient_fcc.h @@ -0,0 +1,115 @@ +/* ---------------------------------------------------------------------- + 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(orient/fcc,FixOrientFCC) + +#else + +#ifndef LMP_FIX_ORIENT_FCC_H +#define LMP_FIX_ORIENT_FCC_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixOrientFCC : public Fix { + public: + struct Nbr { // neighbor info for each owned and ghost atom + int n; // # of closest neighbors (up to 12) + int id[12]; // IDs of each neighbor + // if center atom is owned, these are local IDs + // if center atom is ghost, these are global IDs + double xismooth[12]; // distance weighting factor for each neighbors + double dxi[12][3]; // d order-parameter / dx for each neighbor + double duxi; // d Energy / d order-parameter for atom + }; + + struct Sort { // data structure for sorting to find 12 closest + int id; // ID of neighbor atom + double rsq; // distance between center and neighbor atom + double delta[3]; // displacement between center and neighbor atom + double xismooth; // distance weighting factor + }; + + FixOrientFCC(class LAMMPS *, int, char **); + ~FixOrientFCC(); + int setmask(); + void init(); + void init_list(int, class NeighList *); + void setup(int); + void post_force(int); + void post_force_respa(int, int, int); + double compute_scalar(); + int pack_comm(int, int *, double *, int, int *); + void unpack_comm(int, int, double *); + double memory_usage(); + + private: + int me; + int nlevels_respa; + + int direction_of_motion; // 1 = center shrinks, 0 = center grows + int nstats; // stats output every this many steps + double a; // lattice parameter + double Vxi; // potential value + double uxif_low; // cut-off fraction, low order parameter + double uxif_high; // cut-off fraction, high order parameter + char *xifilename, *chifilename; // file names for 2 crystal orientations + + bool use_xismooth; + double Rxi[12][3],Rchi[12][3],half_xi_chi_vec[2][6][3]; + double xiid,xi0,xi1,xicutoffsq,cutsq,added_energy; + int half_fcc_nn; + + int nmax; // expose 2 per-atom quantities + double **order; // order param and normalized order param + + Nbr *nbr; + Sort *sort; + class NeighList *list; + + void find_best_ref(double *, int, double &, double *); + static int compare(const void *, const void *); +}; + +} + +#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 orient/fcc file open failed + +The fix orient/fcc command could not open a specified file. + +E: Fix orient/fcc file read failed + +The fix orient/fcc command could not read the needed parameters from a +specified file. + +E: Fix orient/fcc found self twice + +The neighbor lists used by fix orient/fcc are messed up. If this +error occurs, it is likely a bug, so send an email to the +"developers"_http://lammps.sandia.gov/authors.html. + +*/ diff --git a/src/MISC/fix_thermal_conductivity.cpp b/src/MISC/fix_thermal_conductivity.cpp new file mode 100644 index 0000000000..b653e34864 --- /dev/null +++ b/src/MISC/fix_thermal_conductivity.cpp @@ -0,0 +1,330 @@ +/* ---------------------------------------------------------------------- + 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/MISC/fix_thermal_conductivity.h b/src/MISC/fix_thermal_conductivity.h new file mode 100644 index 0000000000..22688bb7c0 --- /dev/null +++ b/src/MISC/fix_thermal_conductivity.h @@ -0,0 +1,75 @@ +/* ---------------------------------------------------------------------- + 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/MISC/fix_ttm.cpp b/src/MISC/fix_ttm.cpp new file mode 100644 index 0000000000..abd145a092 --- /dev/null +++ b/src/MISC/fix_ttm.cpp @@ -0,0 +1,685 @@ +/* ---------------------------------------------------------------------- + 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/MISC/fix_ttm.h b/src/MISC/fix_ttm.h new file mode 100644 index 0000000000..f8ed2e626b --- /dev/null +++ b/src/MISC/fix_ttm.h @@ -0,0 +1,156 @@ +/* ---------------------------------------------------------------------- + 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/MISC/fix_viscosity.cpp b/src/MISC/fix_viscosity.cpp new file mode 100644 index 0000000000..ee94a6e6a3 --- /dev/null +++ b/src/MISC/fix_viscosity.cpp @@ -0,0 +1,309 @@ +/* ---------------------------------------------------------------------- + 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/MISC/fix_viscosity.h b/src/MISC/fix_viscosity.h new file mode 100644 index 0000000000..6f694d4f94 --- /dev/null +++ b/src/MISC/fix_viscosity.h @@ -0,0 +1,80 @@ +/* ---------------------------------------------------------------------- + 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. + +*/ diff --git a/src/compute_msd_nongauss.cpp b/src/compute_msd_nongauss.cpp deleted file mode 100644 index eb038cc356..0000000000 --- a/src/compute_msd_nongauss.cpp +++ /dev/null @@ -1,106 +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: Rob Hoy -------------------------------------------------------------------------- */ - -#include "string.h" -#include "compute_msd_nongauss.h" -#include "atom.h" -#include "update.h" -#include "group.h" -#include "domain.h" -#include "fix_store.h" -#include "error.h" - -using namespace LAMMPS_NS; - -/* ---------------------------------------------------------------------- */ - -ComputeMSDNonGauss::ComputeMSDNonGauss(LAMMPS *lmp, int narg, char **arg) : - ComputeMSD(lmp, narg, arg) -{ - size_vector = 3; -} - -/* ---------------------------------------------------------------------- */ - -void ComputeMSDNonGauss::compute_vector() -{ - invoked_vector = update->ntimestep; - - // cm = current center of mass - - double cm[3]; - if (comflag) group->xcm(igroup,masstotal,cm); - else cm[0] = cm[1] = cm[2] = 0.0; - - // dx,dy,dz = displacement of atom from original position - // original unwrapped position is stored by fix - // relative to center of mass if comflag is set - // for triclinic, need to unwrap current atom coord via h matrix - - double **xoriginal = fix->astore; - - double **x = atom->x; - int *mask = atom->mask; - tagint *image = atom->image; - int nlocal = atom->nlocal; - - double *h = domain->h; - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; - - double dx,dy,dz; - int xbox,ybox,zbox; - - double msd[2]; - msd[0] = msd[1] = 0.0; - - if (domain->triclinic == 0) { - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - dx = x[i][0] + xbox*xprd - cm[0] - xoriginal[i][0]; - dy = x[i][1] + ybox*yprd - cm[1] - xoriginal[i][1]; - dz = x[i][2] + zbox*zprd - cm[2] - xoriginal[i][2]; - msd[0] += dx*dx + dy*dy + dz*dz; - msd[1] += (dx*dx + dy*dy + dz*dz)*(dx*dx + dy*dy + dz*dz); - } - - } else { - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - xbox = (image[i] & IMGMASK) - IMGMAX; - ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (image[i] >> IMG2BITS) - IMGMAX; - dx = x[i][0] + h[0]*xbox + h[5]*ybox + h[4]*zbox - - cm[0] - xoriginal[i][0]; - dy = x[i][1] + h[1]*ybox + h[3]*zbox - cm[1] - xoriginal[i][1]; - dz = x[i][2] + h[2]*zbox - cm[2] - xoriginal[i][2]; - msd[0] += dx*dx + dy*dy + dz*dz; - msd[1] += (dx*dx + dy*dy + dz*dz)*(dx*dx + dy*dy + dz*dz); - } - } - - MPI_Allreduce(msd,vector,2,MPI_DOUBLE,MPI_SUM,world); - if (nmsd) { - vector[0] /= nmsd; - vector[1] /= nmsd; - vector[2] = (3*vector[1])/(5*vector[0]*vector[0]) - 1; - } -} diff --git a/src/compute_msd_nongauss.h b/src/compute_msd_nongauss.h deleted file mode 100644 index 86bb653c56..0000000000 --- a/src/compute_msd_nongauss.h +++ /dev/null @@ -1,51 +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 COMPUTE_CLASS - -ComputeStyle(msd/nongauss,ComputeMSDNonGauss) - -#else - -#ifndef LMP_COMPUTE_MSD_NONGAUSS_H -#define LMP_COMPUTE_MSD_NONGAUSS_H - -#include "compute_msd.h" - -namespace LAMMPS_NS { - -class ComputeMSDNonGauss : public ComputeMSD { - public: - ComputeMSDNonGauss(class LAMMPS *, int, char **); - ~ComputeMSDNonGauss() {} - void compute_vector(); -}; - -} - -#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: Could not find compute msd fix ID - -Self-explanatory. - -*/ diff --git a/src/compute_ti.cpp b/src/compute_ti.cpp deleted file mode 100644 index e47bbda556..0000000000 --- a/src/compute_ti.cpp +++ /dev/null @@ -1,241 +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: Sai Jayaraman (University of Notre Dame) -------------------------------------------------------------------------- */ - -#include "mpi.h" -#include "atom.h" -#include "string.h" -#include "compute_ti.h" -#include "update.h" -#include "modify.h" -#include "domain.h" -#include "force.h" -#include "pair.h" -#include "kspace.h" -#include "input.h" -#include "variable.h" -#include "error.h" - -using namespace LAMMPS_NS; - -enum{PAIR,TAIL,KSPACE}; - -/* ---------------------------------------------------------------------- */ - -ComputeTI::ComputeTI(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg) -{ - if (narg < 4) error->all(FLERR,"Illegal compute ti command"); - - peflag = 1; - peratom_flag = 1; - peatomflag = 1; - scalar_flag = 1; - extscalar = 1; - timeflag = 1; - - // terms come in triplets - // changed to quadruplets to include atom type - - nterms = (narg-3) / 4; - if (narg != 4*nterms + 3) error->all(FLERR,"Illegal compute ti command"); - - which = new int[nterms]; - ivar1 = new int[nterms]; - ivar2 = new int[nterms]; - ilo = new int[nterms]; - ihi = new int[nterms]; - var1 = new char*[nterms]; - var2 = new char*[nterms]; - pptr = new Pair*[nterms]; - pstyle = new char*[nterms]; - - for (int m = 0; m < nterms; m++) pstyle[m] = NULL; - - // parse keywords - - nterms = 0; - - int iarg = 3; - while (iarg < narg) { - if (iarg+4 > narg) error->all(FLERR,"Illegal compute ti command"); - if (strcmp(arg[iarg],"kspace") == 0) which[nterms] = KSPACE; - else if (strcmp(arg[iarg],"tail") == 0) which[nterms] = TAIL; - else which[nterms] = PAIR; - - int n = strlen(arg[iarg]) + 1; - pstyle[nterms] = new char[n]; - strcpy(pstyle[nterms],arg[iarg]); - force->bounds(arg[iarg+1],atom->ntypes,ilo[nterms],ihi[nterms]); - iarg += 1; - - if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) { - int n = strlen(&arg[iarg+1][2]) + 1; - var1[nterms] = new char[n]; - strcpy(var1[nterms],&arg[iarg+1][2]); - } else error->all(FLERR,"Illegal compute ti command"); - if (strstr(arg[iarg+2],"v_") == arg[iarg+2]) { - int n = strlen(&arg[iarg+2][2]) + 1; - var2[nterms] = new char[n]; - strcpy(var2[nterms],&arg[iarg+2][2]); - } else error->all(FLERR,"Illegal compute ti command"); - - nterms++; - iarg += 3; - } -} - -/* --------------------------------------------------------------------- */ - -ComputeTI::~ComputeTI() -{ - for (int m = 0; m < nterms; m++) { - delete [] var1[m]; - delete [] var2[m]; - delete [] pstyle[m]; - } - delete [] which; - delete [] ivar1; - delete [] ivar2; - delete [] var1; - delete [] var2; - delete [] ilo; - delete [] ihi; - delete [] pptr; - delete [] pstyle; -} - -/* --------------------------------------------------------------------- */ - -void ComputeTI::init() -{ - // setup and error checks - - for (int m = 0; m < nterms; m++) { - ivar1[m] = input->variable->find(var1[m]); - ivar2[m] = input->variable->find(var2[m]); - if (ivar1[m] < 0 || ivar2 < 0) - error->all(FLERR,"Variable name for compute ti does not exist"); - if (!input->variable->equalstyle(ivar1[m]) || - !input->variable->equalstyle(ivar2[m])) - error->all(FLERR,"Variable for compute ti is invalid style"); - - if (which[m] == PAIR) { - pptr[m] = force->pair_match(pstyle[m],1); - if (pptr[m] == NULL) - error->all(FLERR,"Compute ti pair style does not exist"); - - } else if (which[m] == TAIL) { - if (force->pair == NULL || force->pair->tail_flag == 0) - error->all(FLERR,"Compute ti tail when pair style does not " - "compute tail corrections"); - - } else if (which[m] == KSPACE) { - if (force->kspace == NULL) - error->all(FLERR,"Compute ti kspace style does not exist"); - } - } -} - -/* --------------------------------------------------------------------- */ - -double ComputeTI::compute_scalar() -{ - double eng,engall,value1,value2; - - invoked_scalar = update->ntimestep; - if (update->eflag_global != invoked_scalar) - error->all(FLERR,"Energy was not tallied on needed timestep"); - - double dUdl = 0.0; - - for (int m = 0; m < nterms; m++) { - int total_flag = 0; - if ((ihi[m]-ilo[m])==atom->ntypes) total_flag = 1; - eng = 0.0; - value1 = input->variable->compute_equal(ivar1[m]); - value2 = input->variable->compute_equal(ivar2[m]); - if (value1 == 0.0) continue; - - if (which[m] == PAIR) { - int ntypes = atom->ntypes; - int *mask = atom->mask; - if (total_flag) { - eng = pptr[m]->eng_vdwl + pptr[m]->eng_coul; - MPI_Allreduce(&eng,&engall,1,MPI_DOUBLE,MPI_SUM,world); - } - else { - int nlocal = atom->nlocal; - int npair = nlocal; - int *mask = atom->mask; - int *type = atom->type; - - double *eatom = pptr[m]->eatom; - - if (force->newton) npair += atom->nghost; - for (int i = 0; i < npair; i++) - if ((ilo[m]<=type[i])&(ihi[m]>=type[i])) eng += eatom[i]; - MPI_Allreduce(&eng,&engall,1,MPI_DOUBLE,MPI_SUM,world); - } - dUdl += engall/value1 * value2; - - } else if (which[m] == TAIL) { - double vol = domain->xprd*domain->yprd*domain->zprd; - if (total_flag) - eng = force->pair->etail / vol; - else { - eng = 0; - for (int it = 1; it <= atom->ntypes; it++) { - int jt; - if ((it >= ilo[m])&&(it <=ihi[m])) jt = it; - else jt = ilo[m]; - for (; jt <=ihi[m];jt++) { - if ((force->pair->tail_flag)&&(force->pair->setflag[it][jt])) { - double cut = force->pair->init_one(it,jt); - eng += force->pair->etail_ij; - } - if (it !=jt) eng += force->pair->etail_ij; - } - } - eng /= vol; - } - dUdl += eng/value1 * value2; - - } else if (which[m] == KSPACE) { - int ntypes = atom->ntypes; - int *mask = atom->mask; - if (total_flag) - eng = force->kspace->energy; - else { - int nlocal = atom->nlocal; - int npair = nlocal; - int *mask = atom->mask; - int *type = atom->type; - double *eatom = force->kspace->eatom; - eng = 0; - for(int i = 0; i < nlocal; i++) - if ((ilo[m]<=type[i])&(ihi[m]>=type[i])) - eng += eatom[i]; - MPI_Allreduce(&eng,&engall,1,MPI_DOUBLE,MPI_SUM,world); - eng = engall; - } - dUdl += eng/value1 * value2; - } - } - - scalar = dUdl; - return scalar; -} diff --git a/src/fix_adapt.cpp b/src/fix_adapt.cpp deleted file mode 100644 index 2f1793dfba..0000000000 --- a/src/fix_adapt.cpp +++ /dev/null @@ -1,389 +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 "string.h" -#include "stdlib.h" -#include "fix_adapt.h" -#include "atom.h" -#include "update.h" -#include "modify.h" -#include "force.h" -#include "pair.h" -#include "pair_hybrid.h" -#include "kspace.h" -#include "input.h" -#include "variable.h" -#include "math_const.h" -#include "memory.h" -#include "error.h" - -using namespace LAMMPS_NS; -using namespace FixConst; -using namespace MathConst; - -enum{PAIR,KSPACE,ATOM}; -enum{DIAMETER,CHARGE}; - -/* ---------------------------------------------------------------------- */ - -FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) -{ - if (narg < 5) error->all(FLERR,"Illegal fix adapt command"); - nevery = force->inumeric(FLERR,arg[3]); - if (nevery < 0) error->all(FLERR,"Illegal fix adapt command"); - - // count # of adaptations - - nadapt = 0; - - int iarg = 4; - while (iarg < narg) { - if (strcmp(arg[iarg],"pair") == 0) { - if (iarg+6 > narg) error->all(FLERR,"Illegal fix adapt command"); - nadapt++; - iarg += 6; - } else if (strcmp(arg[iarg],"kspace") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt command"); - nadapt++; - iarg += 2; - } else if (strcmp(arg[iarg],"atom") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal fix adapt command"); - nadapt++; - iarg += 3; - } else break; - } - - if (nadapt == 0) error->all(FLERR,"Illegal fix adapt command"); - adapt = new Adapt[nadapt]; - - // parse keywords - - nadapt = 0; - diamflag = 0; - chgflag = 0; - - iarg = 4; - while (iarg < narg) { - if (strcmp(arg[iarg],"pair") == 0) { - if (iarg+6 > narg) error->all(FLERR,"Illegal fix adapt command"); - adapt[nadapt].which = PAIR; - int n = strlen(arg[iarg+1]) + 1; - adapt[nadapt].pstyle = new char[n]; - strcpy(adapt[nadapt].pstyle,arg[iarg+1]); - n = strlen(arg[iarg+2]) + 1; - adapt[nadapt].pparam = new char[n]; - strcpy(adapt[nadapt].pparam,arg[iarg+2]); - force->bounds(arg[iarg+3],atom->ntypes, - adapt[nadapt].ilo,adapt[nadapt].ihi); - force->bounds(arg[iarg+4],atom->ntypes, - adapt[nadapt].jlo,adapt[nadapt].jhi); - if (strstr(arg[iarg+5],"v_") == arg[iarg+5]) { - n = strlen(&arg[iarg+5][2]) + 1; - adapt[nadapt].var = new char[n]; - strcpy(adapt[nadapt].var,&arg[iarg+5][2]); - } else error->all(FLERR,"Illegal fix adapt command"); - nadapt++; - iarg += 6; - } else if (strcmp(arg[iarg],"kspace") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt command"); - adapt[nadapt].which = KSPACE; - if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) { - int n = strlen(&arg[iarg+1][2]) + 1; - adapt[nadapt].var = new char[n]; - strcpy(adapt[nadapt].var,&arg[iarg+1][2]); - } else error->all(FLERR,"Illegal fix adapt command"); - nadapt++; - iarg += 2; - } else if (strcmp(arg[iarg],"atom") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal fix adapt command"); - adapt[nadapt].which = ATOM; - if (strcmp(arg[iarg+1],"diameter") == 0) { - adapt[nadapt].aparam = DIAMETER; - diamflag = 1; - } else if (strcmp(arg[iarg+1],"charge") == 0) { - adapt[nadapt].aparam = CHARGE; - chgflag = 1; - } else error->all(FLERR,"Illegal fix adapt command"); - if (strstr(arg[iarg+2],"v_") == arg[iarg+2]) { - int n = strlen(&arg[iarg+2][2]) + 1; - adapt[nadapt].var = new char[n]; - strcpy(adapt[nadapt].var,&arg[iarg+2][2]); - } else error->all(FLERR,"Illegal fix adapt command"); - nadapt++; - iarg += 3; - } else break; - } - - // optional keywords - - resetflag = 0; - scaleflag = 0; - - while (iarg < narg) { - if (strcmp(arg[iarg],"reset") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt command"); - if (strcmp(arg[iarg+1],"no") == 0) resetflag = 0; - else if (strcmp(arg[iarg+1],"yes") == 0) resetflag = 1; - else error->all(FLERR,"Illegal fix adapt command"); - iarg += 2; - } else if (strcmp(arg[iarg],"scale") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix adapt command"); - if (strcmp(arg[iarg+1],"no") == 0) scaleflag = 0; - else if (strcmp(arg[iarg+1],"yes") == 0) scaleflag = 1; - else error->all(FLERR,"Illegal fix adapt command"); - iarg += 2; - } else error->all(FLERR,"Illegal fix adapt command"); - } - - // allocate pair style arrays - - int n = atom->ntypes; - for (int m = 0; m < nadapt; m++) - if (adapt[m].which == PAIR) - memory->create(adapt[m].array_orig,n+1,n+1,"adapt:array_orig"); -} - -/* ---------------------------------------------------------------------- */ - -FixAdapt::~FixAdapt() -{ - for (int m = 0; m < nadapt; m++) { - delete [] adapt[m].var; - if (adapt[m].which == PAIR) { - delete [] adapt[m].pstyle; - delete [] adapt[m].pparam; - memory->destroy(adapt[m].array_orig); - } - } - delete [] adapt; -} - -/* ---------------------------------------------------------------------- */ - -int FixAdapt::setmask() -{ - int mask = 0; - mask |= PRE_FORCE; - mask |= POST_RUN; - return mask; -} - -/* ---------------------------------------------------------------------- */ - -void FixAdapt::init() -{ - int i,j; - - // setup and error checks - - anypair = 0; - - for (int m = 0; m < nadapt; m++) { - Adapt *ad = &adapt[m]; - - ad->ivar = input->variable->find(ad->var); - if (ad->ivar < 0) - error->all(FLERR,"Variable name for fix adapt does not exist"); - if (!input->variable->equalstyle(ad->ivar)) - error->all(FLERR,"Variable for fix adapt is invalid style"); - - if (ad->which == PAIR) { - anypair = 1; - - Pair *pair = force->pair_match(ad->pstyle,1); - if (pair == NULL) error->all(FLERR,"Fix adapt pair style does not exist"); - void *ptr = pair->extract(ad->pparam,ad->pdim); - if (ptr == NULL) - error->all(FLERR,"Fix adapt pair style param not supported"); - - ad->pdim = 2; - if (ad->pdim == 0) ad->scalar = (double *) ptr; - if (ad->pdim == 2) ad->array = (double **) ptr; - - // if pair hybrid, test that ilo,ihi,jlo,jhi are valid for sub-style - - if (ad->pdim == 2 && (strcmp(force->pair_style,"hybrid") == 0 || - strcmp(force->pair_style,"hybrid/overlay") == 0)) { - PairHybrid *pair = (PairHybrid *) force->pair; - for (i = ad->ilo; i <= ad->ihi; i++) - for (j = MAX(ad->jlo,i); j <= ad->jhi; j++) - if (!pair->check_ijtype(i,j,ad->pstyle)) - error->all(FLERR,"Fix adapt type pair range is not valid for " - "pair hybrid sub-style"); - } - - } else if (ad->which == KSPACE) { - if (force->kspace == NULL) - error->all(FLERR,"Fix adapt kspace style does not exist"); - kspace_scale = (double *) force->kspace->extract("scale"); - - } else if (ad->which == ATOM) { - if (ad->aparam == DIAMETER) { - if (!atom->radius_flag) - error->all(FLERR,"Fix adapt requires atom attribute diameter"); - } - if (ad->aparam == CHARGE) { - if (!atom->q_flag) - error->all(FLERR,"Fix adapt requires atom attribute charge"); - } - } - } - - // make copy of original pair array values - - for (int m = 0; m < nadapt; m++) { - Adapt *ad = &adapt[m]; - if (ad->which == PAIR && ad->pdim == 2) { - for (i = ad->ilo; i <= ad->ihi; i++) - for (j = MAX(ad->jlo,i); j <= ad->jhi; j++) - ad->array_orig[i][j] = ad->array[i][j]; - } - } -} - -/* ---------------------------------------------------------------------- */ - -void FixAdapt::setup_pre_force(int vflag) -{ - change_settings(); -} - -/* ---------------------------------------------------------------------- */ - -void FixAdapt::pre_force(int vflag) -{ - if (nevery == 0) return; - if (update->ntimestep % nevery) return; - change_settings(); -} - -/* ---------------------------------------------------------------------- */ - -void FixAdapt::post_run() -{ - if (resetflag) restore_settings(); -} - -/* ---------------------------------------------------------------------- - change pair,kspace,atom parameters based on variable evaluation -------------------------------------------------------------------------- */ - -void FixAdapt::change_settings() -{ - int i,j; - - // variable evaluation may invoke computes so wrap with clear/add - - modify->clearstep_compute(); - - for (int m = 0; m < nadapt; m++) { - Adapt *ad = &adapt[m]; - double value = input->variable->compute_equal(ad->ivar); - - // set global scalar or type pair array values - - if (ad->which == PAIR) { - if (ad->pdim == 0) { - if (scaleflag) *ad->scalar = value * ad->scalar_orig; - else *ad->scalar = value; - } else if (ad->pdim == 2) { - if (scaleflag) - for (i = ad->ilo; i <= ad->ihi; i++) - for (j = MAX(ad->jlo,i); j <= ad->jhi; j++) - ad->array[i][j] = value*ad->array_orig[i][j]; - else - for (i = ad->ilo; i <= ad->ihi; i++) - for (j = MAX(ad->jlo,i); j <= ad->jhi; j++) - ad->array[i][j] = value; - } - - // set kspace scale factor - - } else if (ad->which == KSPACE) { - *kspace_scale = value; - - } else if (ad->which == ATOM) { - - // set radius from diameter - // also scale rmass to new value - - if (ad->aparam == DIAMETER) { - int mflag = 0; - if (atom->rmass_flag) mflag = 1; - double density; - - double *radius = atom->radius; - double *rmass = atom->rmass; - int *mask = atom->mask; - int nlocal = atom->nlocal; - - if (mflag == 0) { - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) - radius[i] = 0.5*value; - } else { - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - density = rmass[i] / (4.0*MY_PI/3.0 * - radius[i]*radius[i]*radius[i]); - radius[i] = 0.5*value; - rmass[i] = 4.0*MY_PI/3.0 * - radius[i]*radius[i]*radius[i] * density; - } - } - } else if (ad->aparam == CHARGE) { - double *q = atom->q; - int *mask = atom->mask; - int nlocal = atom->nlocal; - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) q[i] = value; - } - } - } - - modify->addstep_compute(update->ntimestep + nevery); - - // re-initialize pair styles if any PAIR settings were changed - // this resets other coeffs that may depend on changed values, - // and also offset and tail corrections - - if (anypair) force->pair->reinit(); -} - -/* ---------------------------------------------------------------------- - restore pair,kspace.atom parameters to original values -------------------------------------------------------------------------- */ - -void FixAdapt::restore_settings() -{ - for (int m = 0; m < nadapt; m++) { - Adapt *ad = &adapt[m]; - if (ad->which == PAIR) { - if (ad->pdim == 0) *ad->scalar = ad->scalar_orig; - else if (ad->pdim == 2) { - for (int i = ad->ilo; i <= ad->ihi; i++) - for (int j = MAX(ad->jlo,i); j <= ad->jhi; j++) - ad->array[i][j] = ad->array_orig[i][j]; - } - - } else if (ad->which == KSPACE) { - *kspace_scale = 1.0; - - } else if (ad->which == ATOM) { - - } - } - - if (anypair) force->pair->reinit(); -} diff --git a/src/fix_adapt.h b/src/fix_adapt.h deleted file mode 100644 index 88d94d169a..0000000000 --- a/src/fix_adapt.h +++ /dev/null @@ -1,107 +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(adapt,FixAdapt) - -#else - -#ifndef LMP_FIX_ADAPT_H -#define LMP_FIX_ADAPT_H - -#include "fix.h" - -namespace LAMMPS_NS { - -class FixAdapt : public Fix { - public: - int diamflag; // 1 if atom diameters will vary, for AtomVecGranular - int chgflag; - - FixAdapt(class LAMMPS *, int, char **); - ~FixAdapt(); - int setmask(); - void init(); - void setup_pre_force(int); - void pre_force(int); - void post_run(); - - private: - int nadapt,resetflag,scaleflag; - int anypair; - - struct Adapt { - int which,ivar; - char *var; - char *pstyle,*pparam; - int ilo,ihi,jlo,jhi; - int pdim; - double *scalar,scalar_orig; - double **array,**array_orig; - int aparam; - }; - - Adapt *adapt; - double *kspace_scale; - - void change_settings(); - void restore_settings(); -}; - -} - -#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: Variable name for fix adapt does not exist - -Self-explanatory. - -E: Variable for fix adapt is invalid style - -Only equal-style variables can be used. - -E: Fix adapt pair style does not exist - -Self-explanatory - -E: Fix adapt pair style param not supported - -The pair style does not know about the parameter you specified. - -E: Fix adapt type pair range is not valid for pair hybrid sub-style - -Self-explanatory. - -E: Fix adapt kspace style does not exist - -Self-explanatory. - -E: Fix adapt requires atom attribute diameter - -The atom style being used does not specify an atom diameter. - -E: Fix adapt requires atom attribute charge - -The atom style being used does not specify an atom charge. - -*/