1 /* ----------------------------------------------------------------------
2     This is the
3 
4     ██╗     ██╗ ██████╗  ██████╗  ██████╗ ██╗  ██╗████████╗███████╗
5     ██║     ██║██╔════╝ ██╔════╝ ██╔════╝ ██║  ██║╚══██╔══╝██╔════╝
6     ██║     ██║██║  ███╗██║  ███╗██║  ███╗███████║   ██║   ███████╗
7     ██║     ██║██║   ██║██║   ██║██║   ██║██╔══██║   ██║   ╚════██║
8     ███████╗██║╚██████╔╝╚██████╔╝╚██████╔╝██║  ██║   ██║   ███████║
9     ╚══════╝╚═╝ ╚═════╝  ╚═════╝  ╚═════╝ ╚═╝  ╚═╝   ╚═╝   ╚══════╝®
10 
11     DEM simulation engine, released by
12     DCS Computing Gmbh, Linz, Austria
13     http://www.dcs-computing.com, office@dcs-computing.com
14 
15     LIGGGHTS® is part of CFDEM®project:
16     http://www.liggghts.com | http://www.cfdem.com
17 
18     Core developer and main author:
19     Christoph Kloss, christoph.kloss@dcs-computing.com
20 
21     LIGGGHTS® is open-source, distributed under the terms of the GNU Public
22     License, version 2 or later. It is distributed in the hope that it will
23     be useful, but WITHOUT ANY WARRANTY; without even the implied warranty
24     of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. You should have
25     received a copy of the GNU General Public License along with LIGGGHTS®.
26     If not, see http://www.gnu.org/licenses . See also top-level README
27     and LICENSE files.
28 
29     LIGGGHTS® and CFDEM® are registered trade marks of DCS Computing GmbH,
30     the producer of the LIGGGHTS® software and the CFDEM®coupling software
31     See http://www.cfdem.com/terms-trademark-policy for details.
32 
33 -------------------------------------------------------------------------
34     Contributing author and copyright for this file:
35     This file is from LAMMPS, but has been modified. Copyright for
36     modification:
37 
38     Copyright 2012-     DCS Computing GmbH, Linz
39     Copyright 2009-2012 JKU Linz
40 
41     Copyright of original file:
42     LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
43     http://lammps.sandia.gov, Sandia National Laboratories
44     Steve Plimpton, sjplimp@sandia.gov
45 
46     Copyright (2003) Sandia Corporation.  Under the terms of Contract
47     DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
48     certain rights in this software.  This software is distributed under
49     the GNU General Public License.
50 ------------------------------------------------------------------------- */
51 
52 #include <cmath>
53 #include <stdlib.h>
54 #include <string.h>
55 #include "create_atoms.h"
56 #include "atom.h"
57 #include "atom_vec.h"
58 #include "comm.h"
59 #include "modify.h"
60 #include "fix.h"
61 #include "domain.h"
62 #include "lattice.h"
63 #include "region.h"
64 #include "random_park.h"
65 #include "error.h"
66 #include "force.h"
67 
68 // include last to ensure correct macros
69 #include "domain_definitions.h"
70 
71 using namespace LAMMPS_NS;
72 
73 #define EPSILON  1.0e-6
74 
75 enum{BOX,REGION,SINGLE,RANDOM};
76 
77 /* ---------------------------------------------------------------------- */
78 
CreateAtoms(LAMMPS * lmp)79 CreateAtoms::CreateAtoms(LAMMPS *lmp) :
80     Pointers(lmp),
81     seed_char(NULL)
82 {}
83 
84 /* ---------------------------------------------------------------------- */
85 
command(int narg,char ** arg)86 void CreateAtoms::command(int narg, char **arg)
87 {
88   if (domain->box_exist == 0)
89     error->all(FLERR,"Create_atoms command before simulation box is defined");
90   if (modify->nfix_restart_peratom)
91     error->all(FLERR,"Cannot create_atoms after "
92                "reading restart file with per-atom info");
93 
94   all_in = false;
95   all_in_dist = 0.;
96 
97   // parse arguments
98 
99   if (narg < 2) error->all(FLERR,"Illegal create_atoms command");
100   itype = force->inumeric(FLERR,arg[0]);
101   if (itype <= 0 || itype > atom->ntypes)
102     error->all(FLERR,"Invalid atom type in create_atoms command");
103 
104   int iarg = 0;
105   if (strcmp(arg[1],"box") == 0) {
106     style = BOX;
107     iarg = 2;
108   } else if (strcmp(arg[1],"region") == 0) {
109     style = REGION;
110     if (narg < 3) error->all(FLERR,"Illegal create_atoms command");
111     nregion = domain->find_region(arg[2]);
112     if (nregion == -1) error->all(FLERR,
113                                   "Create_atoms region ID does not exist");
114     domain->regions[nregion]->init();
115     iarg = 3;;
116   } else if (strcmp(arg[1],"single") == 0) {
117     style = SINGLE;
118     if (narg < 5) error->all(FLERR,"Illegal create_atoms command");
119     xone[0] = force->numeric(FLERR,arg[2]);
120     xone[1] = force->numeric(FLERR,arg[3]);
121     xone[2] = force->numeric(FLERR,arg[4]);
122     iarg = 5;
123   } else if (strcmp(arg[1],"random") == 0) {
124     style = RANDOM;
125     if (narg < 5) error->all(FLERR,"Illegal create_atoms command");
126     nrandom = force->inumeric(FLERR,arg[2]);
127     if (seed_char)
128         delete [] seed_char;
129     seed_char = new char [strlen(arg[3])+1];
130     strcpy(seed_char, arg[3]);
131     if (strcmp(arg[4],"NULL") == 0) nregion = -1;
132     else {
133       nregion = domain->find_region(arg[4]);
134       if (nregion == -1) error->all(FLERR,
135                                     "Create_atoms region ID does not exist");
136       domain->regions[nregion]->init();
137     }
138     iarg = 5;
139   } else error->all(FLERR,"Illegal create_atoms command");
140 
141   // process optional keywords
142 
143   int scaleflag = 0;
144   remapflag = 0;
145 
146   nbasis = domain->lattice->nbasis;
147   basistype = new int[nbasis];
148   for (int i = 0; i < nbasis; i++) basistype[i] = itype;
149 
150   while (iarg < narg) {
151     if (strcmp(arg[iarg],"basis") == 0) {
152       if (iarg+3 > narg) error->all(FLERR,"Illegal create_atoms command");
153       int ibasis = force->inumeric(FLERR,arg[iarg+1]);
154       itype = force->inumeric(FLERR,arg[iarg+2]);
155       if (ibasis <= 0 || ibasis > nbasis ||
156           itype <= 0 || itype > atom->ntypes)
157         error->all(FLERR,"Invalid basis setting in create_atoms command");
158       basistype[ibasis-1] = itype;
159       iarg += 3;
160     } else if (strcmp(arg[iarg],"remap") == 0) {
161       if (iarg+2 > narg) error->all(FLERR,"Illegal create_atoms command");
162       if (strcmp(arg[iarg+1],"yes") == 0) remapflag = 1;
163       else if (strcmp(arg[iarg+1],"no") == 0) remapflag = 0;
164       else error->all(FLERR,"Illegal create_atoms command");
165       iarg += 2;
166     } else if (strcmp(arg[iarg],"all_in") == 0) {
167       if (iarg+2 > narg) error->all(FLERR,"Illegal create_atoms command");
168       all_in_dist = force->numeric(FLERR,arg[iarg+1]);
169       if(all_in_dist <= 0.)
170             error->all(FLERR,"Illegal create_atoms command, 'all_in' > 0 required");
171       all_in = true;
172       iarg += 2;
173     } else if (strcmp(arg[iarg],"units") == 0) {
174       if (iarg+2 > narg) error->all(FLERR,"Illegal create_atoms command");
175       if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0;
176       else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1;
177       else error->all(FLERR,"Illegal create_atoms command");
178       iarg += 2;
179     } else error->all(FLERR,"Illegal create_atoms command");
180   }
181 
182   // error checks
183 
184   if (style == RANDOM) {
185     if (nrandom < 0) error->all(FLERR,"Illegal create_atoms command");
186   }
187 
188   if(all_in && ! (REGION==style))
189     error->all(FLERR,"Illegal create_atoms command, 'all_in' can only be used together with region");
190 
191   // demand non-none lattice be defined for BOX and REGION
192   // else setup scaling for SINGLE and RANDOM
193   // could use domain->lattice->lattice2box() to do conversion of
194   //   lattice to box, but not consistent with other uses of units=lattice
195   // triclinic remapping occurs in add_single()
196 
197   if (style == BOX || style == REGION) {
198     if (nbasis == 0)
199       error->all(FLERR,"Cannot create atoms with undefined lattice");
200   } else if (scaleflag == 1) {
201     xone[0] *= domain->lattice->xlattice;
202     xone[1] *= domain->lattice->ylattice;
203     xone[2] *= domain->lattice->zlattice;
204   }
205 
206   // set bounds for my proc in sublo[3] & subhi[3]
207   // if periodic and style = BOX or REGION, i.e. using lattice:
208   //   should create exactly 1 atom when 2 images are both "on" the boundary
209   //   either image may be slightly inside/outside true box due to round-off
210   //   if I am lo proc, decrement lower bound by EPSILON
211   //     this will insure lo image is created
212   //   if I am hi proc, decrement upper bound by 2.0*EPSILON
213   //     this will insure hi image is not created
214   //   thus insertion box is EPSILON smaller than true box
215   //     and is shifted away from true boundary
216   //     which is where atoms are likely to be generated
217 
218   triclinic = domain->triclinic;
219 
220   double epsilon[3];
221   if (triclinic) epsilon[0] = epsilon[1] = epsilon[2] = EPSILON;
222   else {
223     epsilon[0] = domain->prd[0] * EPSILON;
224     epsilon[1] = domain->prd[1] * EPSILON;
225     epsilon[2] = domain->prd[2] * EPSILON;
226   }
227 
228   if (triclinic == 0) {
229     sublo[0] = domain->sublo[0]; subhi[0] = domain->subhi[0];
230     sublo[1] = domain->sublo[1]; subhi[1] = domain->subhi[1];
231     sublo[2] = domain->sublo[2]; subhi[2] = domain->subhi[2];
232   } else {
233     sublo[0] = domain->sublo_lamda[0]; subhi[0] = domain->subhi_lamda[0];
234     sublo[1] = domain->sublo_lamda[1]; subhi[1] = domain->subhi_lamda[1];
235     sublo[2] = domain->sublo_lamda[2]; subhi[2] = domain->subhi_lamda[2];
236   }
237 
238   if (style == BOX || style == REGION) {
239     if (domain->xperiodic) {
240       if (comm->myloc[0] == 0) sublo[0] -= epsilon[0];
241       if (comm->myloc[0] == comm->procgrid[0]-1) subhi[0] -= 2.0*epsilon[0];
242     }
243     if (domain->yperiodic) {
244       if (comm->myloc[1] == 0) sublo[1] -= epsilon[1];
245       if (comm->myloc[1] == comm->procgrid[1]-1) subhi[1] -= 2.0*epsilon[1];
246     }
247     if (domain->zperiodic) {
248       if (comm->myloc[2] == 0) sublo[2] -= epsilon[2];
249       if (comm->myloc[2] == comm->procgrid[2]-1) subhi[2] -= 2.0*epsilon[2];
250     }
251   }
252 
253   // add atoms in one of 3 ways
254 
255   bigint natoms_previous = atom->natoms;
256   int nlocal_previous = atom->nlocal;
257 
258   if (style == SINGLE) add_single();
259   else if (style == RANDOM) add_random();
260   else add_lattice();
261 
262   // invoke set_arrays() for fixes that need initialization of new atoms
263 
264   int nlocal = atom->nlocal;
265   for (int m = 0; m < modify->nfix; m++) {
266     Fix *fix = modify->fix[m];
267     if (fix->create_attribute)
268     {
269       fix->pre_set_arrays();
270       for (int i = nlocal_previous; i < nlocal; i++)
271         fix->set_arrays(i);
272     }
273   }
274 
275   // clean up
276 
277   if (domain->lattice) delete [] basistype;
278 
279   // new total # of atoms
280 
281   bigint nblocal = atom->nlocal;
282   MPI_Allreduce(&nblocal,&atom->natoms,1,MPI_LMP_BIGINT,MPI_SUM,world);
283   if (atom->natoms < 0 || atom->natoms > MAXBIGINT)
284     error->all(FLERR,"Too many total atoms");
285 
286   // print status
287 
288   if (comm->me == 0) {
289     if (screen)
290       fprintf(screen,"Created " BIGINT_FORMAT " atoms\n",
291               atom->natoms-natoms_previous);
292     if (logfile)
293       fprintf(logfile,"Created " BIGINT_FORMAT " atoms\n",
294               atom->natoms-natoms_previous);
295   }
296 
297   // reset simulation now that more atoms are defined
298   // add tags for newly created atoms if possible
299   // if global map exists, reset it
300   // change these to MAXTAGINT when allow tagint = bigint
301 
302   if (atom->natoms > MAXSMALLINT) {
303     if (comm->me == 0)
304       error->warning(FLERR,"Total atom count exceeds ID limit, "
305                      "atoms will not have individual IDs");
306     atom->tag_enable = 0;
307   }
308   if (atom->natoms <= MAXSMALLINT) atom->tag_extend();
309 
310   if (atom->map_style) {
311     atom->nghost = 0;
312     atom->map_init();
313     atom->map_set();
314   }
315 
316   // if a molecular system, set nspecial to 0 for new atoms
317   // NOTE: 31May12, don't think this is needed, avec->create_atom() does it
318 
319   //if (atom->molecular) {
320   //  int **nspecial = atom->nspecial;
321   //  for (int i = nlocal_previous; i < atom->nlocal; i++) {
322   //    nspecial[i][0] = 0;
323   //    nspecial[i][1] = 0;
324   //    nspecial[i][2] = 0;
325   //  }
326   //}
327 
328   // error checks on coarsegraining
329   if(force->cg_active())
330     error->cg(FLERR,"create_atoms");
331 }
332 
333 /* ----------------------------------------------------------------------
334    add single atom with coords at xone if it's in my sub-box
335    if triclinic, xone is in lamda coords
336 ------------------------------------------------------------------------- */
337 
add_single()338 void CreateAtoms::add_single()
339 {
340 
341   // remap atom if requested
342 
343   if (remapflag) {
344     tagint imagetmp = ((tagint) IMGMAX << IMG2BITS) |
345       ((tagint) IMGMAX << IMGBITS) | IMGMAX;
346     domain->remap(xone,imagetmp);
347   }
348 
349   // if triclinic, convert to lamda coords (0-1)
350 
351   double lamda[3],*coord;
352   if (triclinic) {
353     domain->x2lamda(xone,lamda);
354     coord = lamda;
355   } else coord = xone;
356 
357   // if atom is in my subbox, create it
358 
359   if (coord[0] >= sublo[0] && coord[0] < subhi[0] &&
360       coord[1] >= sublo[1] && coord[1] < subhi[1] &&
361       coord[2] >= sublo[2] && coord[2] < subhi[2])
362   {
363         atom->avec->create_atom(itype,xone);
364 
365           for (int j = 0; j < modify->nfix; j++)
366             if (modify->fix[j]->create_attribute)
367             {
368                 modify->fix[j]->pre_set_arrays();
369                 modify->fix[j]->set_arrays(atom->nlocal-1);
370             }
371    }
372 }
373 
374 /* ----------------------------------------------------------------------
375    add Nrandom atoms at random locations
376 ------------------------------------------------------------------------- */
377 
add_random()378 void CreateAtoms::add_random()
379 {
380   double xlo,ylo,zlo,xhi,yhi,zhi,zmid;
381   double lamda[3],*coord;
382   double *boxlo=NULL,*boxhi=NULL;
383 
384   // random number generator, same for all procs
385 
386   RanPark *random = new RanPark(lmp,seed_char);
387 
388   // bounding box for atom creation
389   // in real units, even if triclinic
390   // only limit bbox by region if its bboxflag is set (interior region)
391 
392   if (triclinic == 0) {
393     xlo = domain->boxlo[0]; xhi = domain->boxhi[0];
394     ylo = domain->boxlo[1]; yhi = domain->boxhi[1];
395     zlo = domain->boxlo[2]; zhi = domain->boxhi[2];
396     zmid = zlo + 0.5*(zhi-zlo);
397   } else {
398     xlo = domain->boxlo_bound[0]; xhi = domain->boxhi_bound[0];
399     ylo = domain->boxlo_bound[1]; yhi = domain->boxhi_bound[1];
400     zlo = domain->boxlo_bound[2]; zhi = domain->boxhi_bound[2];
401     zmid = zlo + 0.5*(zhi-zlo);
402     boxlo = domain->boxlo_lamda;
403     boxhi = domain->boxhi_lamda;
404   }
405 
406   if (nregion >= 0 && domain->regions[nregion]->bboxflag) {
407     xlo = MAX(xlo,domain->regions[nregion]->extent_xlo);
408     xhi = MIN(xhi,domain->regions[nregion]->extent_xhi);
409     ylo = MAX(ylo,domain->regions[nregion]->extent_ylo);
410     yhi = MIN(yhi,domain->regions[nregion]->extent_yhi);
411     zlo = MAX(zlo,domain->regions[nregion]->extent_zlo);
412     zhi = MIN(zhi,domain->regions[nregion]->extent_zhi);
413   }
414 
415   // generate random positions for each new atom within bounding box
416   // iterate until atom is within region and triclinic simulation box
417   // if final atom position is in my subbox, create it
418 
419   if (xlo > xhi || ylo > yhi || zlo > zhi)
420     error->all(FLERR,"No overlap of box and region for create_atoms");
421 
422   int valid;
423   for (int i = 0; i < nrandom; i++) {
424     while (1) {
425       xone[0] = xlo + random->uniform() * (xhi-xlo);
426       xone[1] = ylo + random->uniform() * (yhi-ylo);
427       xone[2] = zlo + random->uniform() * (zhi-zlo);
428       if (domain->dimension == 2) xone[2] = zmid;
429 
430       valid = 1;
431       if (nregion >= 0)
432       {
433         if(!all_in && domain->regions[nregion]->match(xone[0],xone[1],xone[2]) == 0)
434             valid = 0;
435 
436         else if(all_in && domain->regions[nregion]->match_shrinkby_cut(xone,all_in_dist) == 0)
437             valid = 0;
438       }
439 
440       if (triclinic) {
441         domain->x2lamda(xone,lamda);
442         coord = lamda;
443         if (coord[0] < boxlo[0] || coord[0] >= boxhi[0] ||
444             coord[1] < boxlo[1] || coord[1] >= boxhi[1] ||
445             coord[2] < boxlo[2] || coord[2] >= boxhi[2]) valid = 0;
446       } else coord = xone;
447 
448       if (valid) break;
449     }
450 
451     // if triclinic, coord is now in lamda units
452 
453     if (coord[0] >= sublo[0] && coord[0] < subhi[0] &&
454         coord[1] >= sublo[1] && coord[1] < subhi[1] &&
455         coord[2] >= sublo[2] && coord[2] < subhi[2])
456     {
457         atom->avec->create_atom(itype,xone);
458 
459           for (int j = 0; j < modify->nfix; j++)
460             if (modify->fix[j]->create_attribute)
461             {
462                 modify->fix[j]->pre_set_arrays();
463                 modify->fix[j]->set_arrays(atom->nlocal-1);
464             }
465     }
466   }
467 
468   // clean-up
469 
470   delete random;
471 }
472 
473 /* ----------------------------------------------------------------------
474    add many atoms by looping over lattice
475 ------------------------------------------------------------------------- */
476 
add_lattice()477 void CreateAtoms::add_lattice()
478 {
479   // convert 8 corners of my subdomain from box coords to lattice coords
480   // for orthogonal, use corner pts of my subbox
481   // for triclinic, use bounding box of my subbox
482   // xyz min to max = bounding box around the domain corners in lattice space
483 
484   double bboxlo[3],bboxhi[3];
485 
486   if (triclinic == 0) {
487     bboxlo[0] = domain->sublo[0]; bboxhi[0] = domain->subhi[0];
488     bboxlo[1] = domain->sublo[1]; bboxhi[1] = domain->subhi[1];
489     bboxlo[2] = domain->sublo[2]; bboxhi[2] = domain->subhi[2];
490   } else domain->bbox(domain->sublo_lamda,domain->subhi_lamda,bboxlo,bboxhi);
491 
492   double xmin,ymin,zmin,xmax,ymax,zmax;
493   xmin = ymin = zmin = BIG;
494   xmax = ymax = zmax = -BIG;
495 
496   domain->lattice->bbox(1,bboxlo[0],bboxlo[1],bboxlo[2],
497                         xmin,ymin,zmin,xmax,ymax,zmax);
498   domain->lattice->bbox(1,bboxhi[0],bboxlo[1],bboxlo[2],
499                         xmin,ymin,zmin,xmax,ymax,zmax);
500   domain->lattice->bbox(1,bboxlo[0],bboxhi[1],bboxlo[2],
501                         xmin,ymin,zmin,xmax,ymax,zmax);
502   domain->lattice->bbox(1,bboxhi[0],bboxhi[1],bboxlo[2],
503                         xmin,ymin,zmin,xmax,ymax,zmax);
504   domain->lattice->bbox(1,bboxlo[0],bboxlo[1],bboxhi[2],
505                         xmin,ymin,zmin,xmax,ymax,zmax);
506   domain->lattice->bbox(1,bboxhi[0],bboxlo[1],bboxhi[2],
507                         xmin,ymin,zmin,xmax,ymax,zmax);
508   domain->lattice->bbox(1,bboxlo[0],bboxhi[1],bboxhi[2],
509                         xmin,ymin,zmin,xmax,ymax,zmax);
510   domain->lattice->bbox(1,bboxhi[0],bboxhi[1],bboxhi[2],
511                         xmin,ymin,zmin,xmax,ymax,zmax);
512 
513   // ilo:ihi,jlo:jhi,klo:khi = loop bounds for lattice overlap of my subbox
514   // overlap = any part of a unit cell (face,edge,pt) in common with my subbox
515   // in lattice space, subbox is a tilted box
516   // but bbox of subbox is aligned with lattice axes
517   // so ilo:khi unit cells should completely tile bounding box
518   // decrement lo, increment hi to avoid round-off issues in lattice->bbox(),
519   //   which can lead to missing atoms in rare cases
520   // extra decrement of lo if min < 0, since static_cast(-1.5) = -1
521 
522   int ilo,ihi,jlo,jhi,klo,khi;
523   ilo = static_cast<int> (xmin) - 1;
524   jlo = static_cast<int> (ymin) - 1;
525   klo = static_cast<int> (zmin) - 1;
526   ihi = static_cast<int> (xmax) + 1;
527   jhi = static_cast<int> (ymax) + 1;
528   khi = static_cast<int> (zmax) + 1;
529 
530   if (xmin < 0.0) ilo--;
531   if (ymin < 0.0) jlo--;
532   if (zmin < 0.0) klo--;
533 
534   // iterate on 3d periodic lattice of unit cells using loop bounds
535   // iterate on nbasis atoms in each unit cell
536   // convert lattice coords to box coords
537   // add atom if it meets all criteria
538 
539   double **basis = domain->lattice->basis;
540   double x[3],lamda[3];
541   double *coord;
542 
543   int i,j,k,m;
544   for (k = klo; k <= khi; k++)
545     for (j = jlo; j <= jhi; j++)
546       for (i = ilo; i <= ihi; i++)
547         for (m = 0; m < nbasis; m++) {
548 
549           x[0] = i + basis[m][0];
550           x[1] = j + basis[m][1];
551           x[2] = k + basis[m][2];
552 
553           // convert from lattice coords to box coords
554 
555           domain->lattice->lattice2box(x[0],x[1],x[2]);
556 
557           // if a region was specified, test if atom is in it
558 
559           if (style == REGION)
560           {
561             if (!all_in && !domain->regions[nregion]->match(x[0],x[1],x[2])) continue;
562             if (all_in && !domain->regions[nregion]->match_shrinkby_cut(x,all_in_dist)) continue;
563           }
564 
565           // test if atom is in my subbox
566 
567           if (triclinic) {
568             domain->x2lamda(x,lamda);
569             coord = lamda;
570           } else coord = x;
571 
572           if (coord[0] < sublo[0] || coord[0] >= subhi[0] ||
573               coord[1] < sublo[1] || coord[1] >= subhi[1] ||
574               coord[2] < sublo[2] || coord[2] >= subhi[2]) continue;
575 
576           // add the atom to my list of atoms
577 
578           atom->avec->create_atom(basistype[m],x);
579 
580           for (int k = 0; k < modify->nfix; k++)
581             if (modify->fix[k]->create_attribute)
582             {
583                 modify->fix[k]->pre_set_arrays();
584                 modify->fix[k]->set_arrays(atom->nlocal-1);
585             }
586         }
587 }
588