git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@2676 f3b2605a-c512-4ea7-a41b-209d697bcdaa

src/compute_displace_atom.cpp

@@ -78,7 +78,7 @@ void ComputeDisplaceAtom::compute_peratom()
 
   // grow local displacement array if necessary
 
-  if (atom->nmax > nmax) {
+  if (atom->nlocal > nmax) {
     memory->destroy_2d_double_array(displace);
     nmax = atom->nmax;
     displace =

src/compute_pe_atom.cpp

@@ -80,6 +80,7 @@ void ComputePEAtom::compute_peratom()
     error->all("Per-atom energy was not tallied on needed timestep");
 
   // grow local energy array if necessary
+  // needs to be atom->nmax in length
 
   if (atom->nmax > nmax) {
     memory->sfree(energy);

src/compute_stress_atom.cpp

@@ -94,6 +94,7 @@ void ComputeStressAtom::compute_peratom()
     error->all("Per-atom virial was not tallied on needed timestep");
 
   // grow local stress array if necessary
+  // needs to be atom->nmax in length
 
   if (atom->nmax > nmax) {
     memory->destroy_2d_double_array(stress);

src/fix_evaporate.cpp

@@ -0,0 +1,191 @@
+/* ----------------------------------------------------------------------
+   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 "random_park.h"
+#include "random_mars.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+#define MIN(A,B) ((A) < (B)) ? (A) : (B)
+#define MAX(A,B) ((A) > (B)) ? (A) : (B)
+
+/* ---------------------------------------------------------------------- */
+
+FixEvaporate::FixEvaporate(LAMMPS *lmp, int narg, char **arg) :
+  Fix(lmp, narg, arg)
+{
+  if (narg != 7) error->all("Illegal fix evaporate command");
+
+  scalar_flag = 1;
+  scalar_vector_freq = 1;
+  extscalar = 0;
+
+  nevery = atoi(arg[3]);
+  nflux = atoi(arg[4]);
+  iregion = domain->find_region(arg[5]);
+  int seed = atoi(arg[6]);
+
+  if (nevery <= 0 || nflux <= 0) error->all("Illegal fix evaporate command");
+  if (iregion == -1) error->all("Fix evaporate region ID does not exist");
+  if (seed <= 0) error->all("Illegal fix evaporate command");
+
+  // random number generator, same for all procs
+
+  random = new RanPark(lmp,seed);
+
+  // set up reneighboring
+
+  force_reneighbor = 1;
+  next_reneighbor = (update->ntimestep/nevery)*nevery + nevery;
+  ndeleted = 0;
+
+  nmax = 0;
+  list = NULL;
+  mark = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixEvaporate::~FixEvaporate()
+{
+  delete random;
+  memory->sfree(list);
+  memory->sfree(mark);
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixEvaporate::setmask()
+{
+  int mask = 0;
+  mask |= PRE_EXCHANGE;
+  return mask;
+}
+
+/* ----------------------------------------------------------------------
+   perform particle deletion
+   done before exchange, borders, reneighbor
+   so that ghost atoms and neighbor lists will be correct
+------------------------------------------------------------------------- */
+
+void FixEvaporate::pre_exchange()
+{
+  int i,iwhichglobal,iwhichlocal;
+
+  if (update->ntimestep != next_reneighbor) return;
+
+  // grow list and mark arrays if necessary
+
+  if (atom->nlocal > nmax) {
+    memory->sfree(list);
+    memory->sfree(mark);
+    nmax = atom->nmax;
+    list = (int *) memory->smalloc(nmax*sizeof(int),"fix/evaporate:list");
+    mark = (int *) memory->smalloc(nmax*sizeof(int),"fix/evaporate:mark");
+  }
+
+  // nall = # of deletable atoms in region
+  // nbefore = # on procs before me
+
+  double **x = atom->x;
+  int *mask = atom->mask;
+  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;
+  
+  // nwhack = total number of atoms to delete
+  // choose atoms randomly across all procs and mark them for deletion
+  // shrink local list of candidates as my atoms get marked for deletion
+
+  int nwhack = MIN(nflux,nall);
+
+  for (i = 0; i < nlocal; i++) mark[i] = 0;
+  
+  for (i = 0; i < nwhack; i++) {
+    iwhichglobal = static_cast<int> (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--;
+    }
+    nall--;
+  }
+  
+  // 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);
+      atom->nlocal--;
+    }
+  }
+
+  // reset global natoms
+  // if global map exists, reset it
+
+  atom->natoms -= nwhack;
+  if (nwhack && atom->map_style) {
+    atom->map_init();
+    atom->map_set();
+  }
+
+  // statistics
+  
+  ndeleted += nwhack;
+  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;
+}

src/fix_evaporate.h

@@ -0,0 +1,42 @@
+/* ----------------------------------------------------------------------
+   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.
+------------------------------------------------------------------------- */
+
+#ifndef FIX_EVAPORATE_H
+#define FIX_EVAPORATE_H
+
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixEvaporate : public Fix {
+ public:
+  FixEvaporate(class LAMMPS *, int, char **);
+  ~FixEvaporate();
+  int setmask();
+  void pre_exchange();
+  double compute_scalar();
+  double memory_usage();
+
+ private:
+  int nevery,nflux,iregion;
+  int ndeleted;
+
+  int nmax;
+  int *list,*mark;
+
+  class RanPark *random;
+};
+
+}
+
+#endif

src/read_data.cpp

@@ -441,7 +441,7 @@ void ReadData::atoms()
   if (natoms != atom->natoms) error->all("Did not assign all atoms correctly");
   
   // if any atom ID < 0, error
-  // if all atom IDs == 0, tag_enable = 0
+  // if all atom IDs = 0, tag_enable = 0
   // if any atom ID > 0, error if any atom ID == 0
   // not checking if atom IDs > natoms or are unique
   

src/style.h

@@ -156,6 +156,7 @@ DumpStyle(xyz,DumpXYZ)
 #include "fix_dt_reset.h"
 #include "fix_efield.h"
 #include "fix_enforce2d.h"
+#include "fix_evaporate.h"
 #include "fix_gravity.h"
 #include "fix_gyration.h"
 #include "fix_heat.h"
@@ -216,6 +217,7 @@ FixStyle(drag,FixDrag)
 FixStyle(dt/reset,FixDtReset)
 FixStyle(efield,FixEfield)
 FixStyle(enforce2d,FixEnforce2D)
+FixStyle(evaporate,FixEvaporate)
 FixStyle(gravity,FixGravity)
 FixStyle(gyration,FixGyration)
 FixStyle(heat,FixHeat)