Commit 825facad authored by jrgissing's avatar jrgissing
Browse files

create_atoms subset: improvements based on Steve's suggestions

parent 72bbabce

doc/src/create_atoms.txt

+3 −3
Original line number Diff line number Diff line
@@ -24,14 +24,14 @@ style = {box} or {region} or {single} or {random} :l 
    seed = random # seed (positive integer)

    region-ID = create atoms within this region, use NULL for entire simulation box :pre

zero or more keyword/value pairs may be appended :l

keyword = {mol} or {basis} or {randnpos} or {remap} or {var} or {set} or {units} :l

keyword = {mol} or {basis} or {subset} or {remap} or {var} or {set} or {units} :l

  {mol} value = template-ID seed

    template-ID = ID of molecule template specified in a separate "molecule"_molecule.html command

    seed = random # seed (positive integer)

  {basis} values = M itype

    M = which basis atom

    itype = atom type (1-N) to assign to this basis atom

  {randnpos} value = N = number of particles to add at lattice points

  {subset} value = N = number of particles to add at lattice points

  {remap} value = {yes} or {no}

  {var} value = name = variable name to evaluate for test of atom creation

  {set} values = dim name
@@ -194,7 +194,7 @@ command for specifics on how basis atoms are defined for the unit cell 
of the lattice.  By default, all created atoms are assigned the

argument {type} as their atom type.



The {randnpos} keyword can be used in conjunction with the {box} or

The {subset} keyword can be used in conjunction with the {box} or

{region} styles to limit the total number of particles inserted. The

specified number of particles N are placed randomly on the available

lattice points.

src/create_atoms.cpp

+74 −59
Original line number Diff line number Diff line
@@ -35,6 +35,7 @@ 
#include "math_extra.h"

#include "math_const.h"

#include "error.h"

#include "memory.h"



using namespace LAMMPS_NS;

using namespace MathConst;
@@ -109,12 +110,12 @@ void CreateAtoms::command(int narg, char **arg) 
  remapflag = 0;

  mode = ATOM;

  int molseed;

  int randnseed;

  int subsetseed;

  varflag = 0;

  vstr = xstr = ystr = zstr = NULL;

  quatone[0] = quatone[1] = quatone[2] = 0.0;

  nlattpts = randnpos = randnflag = 0;

  Nmask = NULL;

  nlatt = nsubset = subsetflag = 0;

  flag = NULL;



  nbasis = domain->lattice->nbasis;

  basistype = new int[nbasis];
@@ -197,15 +198,15 @@ void CreateAtoms::command(int narg, char **arg) 
      MathExtra::norm3(axisone);

      MathExtra::axisangle_to_quat(axisone,thetaone,quatone);

      iarg += 5;

    } else if (strcmp(arg[iarg],"randnpos") == 0) {

    } else if (strcmp(arg[iarg],"subset") == 0) {

      if (iarg+3 > narg) error->all(FLERR,"Illegal create_atoms command");

      randnpos = force->inumeric(FLERR,arg[iarg+1]);

      randnseed = force->inumeric(FLERR,arg[iarg+2]);

      if (randnpos > 0) randnflag = 1;

      nsubset = force->inumeric(FLERR,arg[iarg+1]);

      subsetseed = force->inumeric(FLERR,arg[iarg+2]);

      if (nsubset > 0) subsetflag = 1;

      else {

        if (me == 0) error->warning(FLERR,"Specifying an 'randnpos' value of "

                                 "'0' is equivalent to no 'randnpos' keyword");

        randnflag = 0;

        if (me == 0) error->warning(FLERR,"Specifying an 'subset' value of "

                                 "'0' is equivalent to no 'subset' keyword");

        subsetflag = 0;

      }

      iarg += 3;

    } else error->all(FLERR,"Illegal create_atoms command");
@@ -245,8 +246,8 @@ void CreateAtoms::command(int narg, char **arg) 
    ranmol = new RanMars(lmp,molseed+me);

  }



  if (randnflag) {

    ranbox = new RanMars(lmp,randnseed+me);

  if (subsetflag) {

    ranlatt = new RanMars(lmp,subsetseed+me);

  }



  // error check and further setup for variable test
@@ -545,7 +546,7 @@ void CreateAtoms::command(int narg, char **arg) 
  delete [] xstr;

  delete [] ystr;

  delete [] zstr;

  delete [] Nmask;

  memory->destroy(flag);



  // print status
@@ -779,14 +780,14 @@ void CreateAtoms::add_lattice() 


  int i,j,k,m;



  // one pass for default mode, two passes for randnpos mode:

  // one pass for default mode, two passes for subset mode:

  // first pass: count how many particles will be inserted

  // second pass: filter to N number of particles (and insert)

  int maskcntr = 0;

  int iflag = 0;

  int npass = 1;

  if (randnflag) npass = 2;

  if (subsetflag) npass = 2;

  for (int pass = 0; pass < npass; pass++) {

    if (pass == 1) lattice_mask();

    if (pass == 1) get_subset();

    for (k = klo; k <= khi; k++)

      for (j = jlo; j <= jhi; j++)

        for (i = ilo; i <= ihi; i++)
@@ -821,9 +822,9 @@ void CreateAtoms::add_lattice() 
                coord[2] < sublo[2] || coord[2] >= subhi[2]) continue;



            // add the atom or entire molecule to my list of atoms

            if (randnflag && pass == 0) nlattpts++;

            if (subsetflag && pass == 0) nlatt++;

            else {

              if (!randnflag || Nmask[maskcntr++] == 1)

              if (!subsetflag || flag[iflag++] == 1)

                if (mode == ATOM) atom->avec->create_atom(basistype[m],x);

                else add_molecule(x);

            }
@@ -835,64 +836,78 @@ void CreateAtoms::add_lattice() 
   define a random mask to insert only N particles on lattice

------------------------------------------------------------------------- */



void CreateAtoms::lattice_mask()

void CreateAtoms::get_subset()

{

  int i,j,temp,irand,nboxme;

  enum{ATOMS,HOLES};

  int i,j,temp,irand,mysubset,npicks,pickmode;

  double myrand;

  if (nprocs > 1) {

    int allnlattpts[nprocs];

    int allnboxmes[nprocs];

    int allnlatts[nprocs];

    int allsubsets[nprocs];

    double proc_sects[nprocs];



   MPI_Allgather(&nlattpts, 1, MPI_INT, &allnlattpts, 1, MPI_INT, world);

    MPI_Allgather(&nlatt, 1, MPI_INT, &allnlatts, 1, MPI_INT, world);



    if (me == 0) {

     int total_lattpts = 0;

      int ntotal = 0;

      for (i = 0; i < nprocs; i++)

       total_lattpts += allnlattpts[i];

        ntotal += allnlatts[i];



     if (randnpos > total_lattpts)

      if (nsubset > ntotal)

         error->one(FLERR,"Attempting to insert more particles than available lattice points");



      // using proc 0, let's randomly choose how many lattice points filled for each proc

      int *allNmask = new int[total_lattpts]();

      // define regions of unity based on a proc's fraction of total lattice points

      proc_sects[0] = allnlatts[0]/ntotal;

      for (i = 1; i < nprocs; i++)

        proc_sects[i] = proc_sects[i-1] + (double) allnlatts[i] / (double) ntotal;



      if (nsubset > ntotal/2) {

        pickmode = HOLES;

        npicks = ntotal - nsubset;

      } else {

        pickmode = ATOMS;

        npicks = nsubset;

      }



      int ilatt = 0;

      for (i = 0; i < nprocs; i++)

        for (j = 0; j < allnlattpts[i]; j++)

          allNmask[ilatt++] = i;

        allsubsets[i] = 0;



      // shuffle lattice points, labelled by proc

      for (i = total_lattpts-1; i > 0; --i) {

        irand = floor(ranbox->uniform()*(i+1));

        temp = allNmask[i];

        allNmask[i] = allNmask[irand];

        allNmask[irand] = temp;

      for (i = 0; i < npicks; i++) {

        myrand = ranlatt->uniform();

        for (j = 0; j < nprocs; j++)

          if (myrand < proc_sects[j]) {

            allsubsets[j]++;

            break;

          }

      }



      if (pickmode == HOLES)

        for (i = 0; i < nprocs; i++)

        allnboxmes[i] = 0;



      for (i = 0; i < randnpos; i++)

        allnboxmes[allNmask[i]]++;

          allsubsets[i] = allnlatts[i] - allsubsets[i];

    }



    MPI_Scatter(&allnboxmes, 1, MPI_INT, &nboxme, 1, MPI_INT, 0, world);

    MPI_Scatter(&allsubsets, 1, MPI_INT, &mysubset, 1, MPI_INT, 0, world);

  } else {

    if (randnpos > nlattpts)

    if (nsubset > nlatt)

       error->one(FLERR,"Attempting to insert more particles than available lattice points");

    nboxme = randnpos;

    mysubset = nsubset;

  }



  // each processor chooses its random lattice points

  Nmask = new int[nlattpts]();

  for (i = 0; i < nboxme; i++)

    Nmask[i] = 1;

  memory->create(flag,nlatt,"create_atoms:flag");



  for (i = 0; i < nlatt; i++)

    if (i < mysubset)

      flag[i] = 1;

    else

      flag[i] = 0;



  // shuffle filled lattice points

  for (i = nlattpts-1; i > 0; --i) {

    irand = floor(ranbox->uniform()*(i+1));

    temp = Nmask[i];

    Nmask[i] = Nmask[irand];

    Nmask[irand] = temp;

  for (i = nlatt-1; i > 0; --i) {

    irand = floor(ranlatt->uniform()*(i+1));

    temp = flag[i];

    flag[i] = flag[irand];

    flag[irand] = temp;

  }

}

src/create_atoms.h

+5 −5
Original line number Diff line number Diff line
@@ -40,12 +40,12 @@ class CreateAtoms : protected Pointers { 
  char *vstr,*xstr,*ystr,*zstr;

  char *xstr_copy,*ystr_copy,*zstr_copy;



  int nlattpts,randnpos,randnflag;

  int *Nmask;                 // used to insert N number of particles on lattice

  int nlatt,nsubset,subsetflag;

  int *flag;                  // used to insert N number of particles on lattice



  class Molecule *onemol;

  class RanMars *ranmol;

  class RanMars *ranbox;

  class RanMars *ranlatt;



  int triclinic;

  double sublo[3],subhi[3];   // epsilon-extended proc sub-box for adding atoms
@@ -53,7 +53,7 @@ class CreateAtoms : protected Pointers { 
  void add_single();

  void add_random();

  void add_lattice();

  void lattice_mask();

  void get_subset();

  void add_molecule(double *, double * = NULL);

  int vartest(double *);      // evaluate a variable with new atom position

};

CRA Git | Maintained and supported by SUSTech CRA and CCSE