C++ AtomVec类代码示例

void writefile(AtomVec& atoms, OriginBox& obox){
    // The .xyz format is simple:
    // Line 1: [Number of atoms]
    // Line 2: [Comment line, whatever you want, left blank here]
    // Line 3: [element type, here C for carbon]\t[x]\t[y]\t[z]
    // Line 4: [element type, here C for carbon]\t[x]\t[y]\t[z]
    // ...
    // Line N+1: [element type, here C for carbon]\t[x]\t[y]\t[z]
    // Line N+2: [element type, here C for carbon]\t[x]\t[y]\t[z]
    // Line N+3: [Number of atoms]
    // Line N+4: [Comment line, whatever you want, left blank here]
    // Line N+5: [element type, here C for carbon]\t[x]\t[y]\t[z]
    // ...
    // And so on. each set of N atoms/coordinates corresponds to a "frame",
    // which then make a movie. There must be the same N in each frame for VMD.
    // Note that if this is compiled as a 2D simulation, it will leave out
    // the z coordinate, and VMD can't handle that.
    
    ofstream outf;
    outf.open("packing.xyz", ios::out);
    outf.precision(24);
    
    outf << atoms.size() << endl;
    outf << "L=" << obox.L() << endl; // blank line for comment
    for(uint i=0; i<atoms.size(); i++){
        if(i < Ns){
            outf << "C";
        } else {
            outf << "O";
        }
        Vec normloc = obox.diff(Vec::Zero(), atoms[i].x);
        for(uint j=0; j<NDIM; j++){
            outf << "\t" << normloc[j];
        }
        outf << endl;
    };

    // Unnecessary extra:
    // Write a "tcl" file with the box boundaries
    // the "tcl" file is made specifically for VMD
    ofstream pbcfile;
    pbcfile.open("packing.tcl", ios::out);
    pbcfile << "set cell [pbc set {";
    for(uint j=0; j<NDIM; j++){
        pbcfile << obox.box_shape()[j] << " ";
    }
    pbcfile << "} -all];\n";
    pbcfile << "pbc box -toggle -center origin -color red;\n";

    pbcfile << "set natoms [atomselect 0 \"name C\";];\n"
            << "$natoms set radius " << (sigma/2.0) << ";\n"
            << "set natoms [atomselect 0 \"name O\";];\n"
            << "$natoms set radius " << (sigmal/2.0) << ";\n";


    // Now you should be able to run "vmd -e LJatoms-pbc.tcl LJatoms.xyz"
    // and it will show you the movie and also the bounding box
    // if you have .vmdrc in that same directory, you should also be able
    // to toggle the box with the "b" button
};

void writefile(ofstream& outf, AtomVec& atoms, Box& bx){
    // The .xyz format is simple:
    // Line 1: [Number of atoms]
    // Line 2: [Comment line, whatever you want, left blank here]
    // Line 3: [element type, here C for carbon]\t[x]\t[y]\t[z]
    // Line 4: [element type, here C for carbon]\t[x]\t[y]\t[z]
    // ...
    // Line N+1: [element type, here C for carbon]\t[x]\t[y]\t[z]
    // Line N+2: [element type, here C for carbon]\t[x]\t[y]\t[z]
    // Line N+3: [Number of atoms]
    // Line N+4: [Comment line, whatever you want, left blank here]
    // Line N+5: [element type, here C for carbon]\t[x]\t[y]\t[z]
    // ...
    // And so on. each set of N atoms/coordinates corresponds to a "frame", 
    // which then make a movie. There must be the same N in each frame for VMD.
    // Note that if this is compiled as a 2D simulation, it will leave out
    // the z coordinate, and VMD can't handle that.
    
    outf << atoms.size() << endl;
    outf << endl; // blank line for comment
    for(uint i=0; i<atoms.size(); i++){
        outf << "C";
        Vec normloc = bx.diff(Vec::Zero(), atoms[i].x);
        for(uint j=0; j<NDIM; j++){
            outf << "\t" << normloc[j];
        }
        outf << endl;
    };
};

void CommBrick::reverse_comm()
{
  int n;
  MPI_Request request;
  AtomVec *avec = atom->avec;
  double **f = atom->f;
  double *buf;

  // exchange data with another proc
  // if other proc is self, just copy
  // if comm_f_only set, exchange or copy directly from f, don't pack

  for (int iswap = nswap-1; iswap >= 0; iswap--) {
    if (sendproc[iswap] != me) {
      if (comm_f_only) {
        if (size_reverse_recv[iswap])
          MPI_Irecv(buf_recv,size_reverse_recv[iswap],MPI_DOUBLE,
                    sendproc[iswap],0,world,&request);
        if (size_reverse_send[iswap]) {
          if (size_reverse_send[iswap]) buf = f[firstrecv[iswap]];
          else buf = NULL;
          MPI_Send(buf,size_reverse_send[iswap],MPI_DOUBLE,
                   recvproc[iswap],0,world);
        }
        if (size_reverse_recv[iswap]) MPI_Wait(&request,MPI_STATUS_IGNORE);
      } else {
        if (size_reverse_recv[iswap])
          MPI_Irecv(buf_recv,size_reverse_recv[iswap],MPI_DOUBLE,
                    sendproc[iswap],0,world,&request);
        n = avec->pack_reverse(recvnum[iswap],firstrecv[iswap],buf_send);
        if (n) MPI_Send(buf_send,n,MPI_DOUBLE,recvproc[iswap],0,world);
        if (size_reverse_recv[iswap]) MPI_Wait(&request,MPI_STATUS_IGNORE);
      }
      avec->unpack_reverse(sendnum[iswap],sendlist[iswap],buf_recv);

    } else {
      if (comm_f_only) {
        if (sendnum[iswap])
          avec->unpack_reverse(sendnum[iswap],sendlist[iswap],
                               f[firstrecv[iswap]]);
      } else {
        avec->pack_reverse(recvnum[iswap],firstrecv[iswap],buf_send);
        avec->unpack_reverse(sendnum[iswap],sendlist[iswap],buf_send);
      }
    }
  }
}

//.........这里部分代码省略.........
      rn = random->uniform();
      h = (hi_current-shift_randompos(radtmp)) - rn * (hi_current-lo_current-2.*shift_randompos(radtmp));
      attempt++;
      xyz_random(h,coord,radtmp);
      for (i = 0; i < nnear; i++) {
        if(overlaps_xnear_i(coord,radtmp,xnear,i)) break; 
      }
      if (i == nnear) {
	    success = 1;
	    break;
      }
    }
    if (success) {
        nnear=insert_in_xnear(xnear,nnear,coord,radtmp); 
    } else break;
  }

  // warn if not all insertions were performed

  ninserted += nnear-nprevious;
  if (nnear - nprevious < nnew && me == 0)
    error->warning("Less insertions than requested");

  // 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 atom
  // type, diameter, density set from fix parameters
  // group mask set to "all" plus fix group
  // z velocity set to what velocity would be if particle
  //   had fallen from top of insertion region
  // this gives continuous stream of atoms
  // set npartner for new atom to 0 (assume not touching any others)

  AtomVec *avec = atom->avec;
  int j,m,flag;
  double denstmp,vxtmp,vytmp,vztmp;
  double g = grav; 
  //double g = 1.; //originally
  double *sublo = domain->sublo;
  double *subhi = domain->subhi;
  int b_id;  

  int nfix = modify->nfix;
  Fix **fix = modify->fix;

  for (i = nprevious; i < nnear; i++) {
    coord[0] = xnear[i][0];
    coord[1] = xnear[i][1];
    coord[2] = xnear[i][2];
    radtmp = xnear[i][3];
    b_id = static_cast<int>(xnear[i][4]);
    denstmp = rand_pour(density_lo,density_hi,density_ran_style)*density_scaling();

    calc_insert_velocities(i,g,xnear,vxtmp,vytmp,vztmp);

    flag = 0;
    if (coord[0] >= sublo[0] && coord[0] < subhi[0] &&
	coord[1] >= sublo[1] && coord[1] < subhi[1] &&
	coord[2] >= sublo[2] && coord[2] < subhi[2]) flag = 1;
    else if (domain->dimension == 3 && coord[2] >= domain->boxhi[2] &&
	     comm->myloc[2] == comm->procgrid[2]-1 &&
	     coord[0] >= sublo[0] && coord[0] < subhi[0] &&
	     coord[1] >= sublo[1] && coord[1] < subhi[1]) flag = 1;
    else if (domain->dimension == 2 && coord[1] >= domain->boxhi[1] &&
	     comm->myloc[1] == comm->procgrid[1]-1 &&
	     coord[0] >= sublo[0] && coord[0] < subhi[0]) flag = 1;

void WriteRestart::write(char *file)
{
  // special case where reneighboring is not done in integrator
  //   on timestep restart file is written (due to build_once being set)
  // if box is changing, must be reset, else restart file will have
  //   wrong box size and atoms will be lost when restart file is read
  // other calls to pbc and domain and comm are not made,
  //   b/c they only make sense if reneighboring is actually performed

  if (neighbor->build_once) domain->reset_box();

  // natoms = sum of nlocal = value to write into restart file
  // if unequal and thermo lostflag is "error", don't write restart file

  bigint nblocal = atom->nlocal;

  if(region)
  {
      nblocal = 0.;
      for (int i = 0; i < atom->nlocal; i++)
          if(region->match(atom->x[i][0],atom->x[i][1],atom->x[i][2]))
             nblocal += 1;
  }
  
  MPI_Allreduce(&nblocal,&natoms,1,MPI_LMP_BIGINT,MPI_SUM,world);
  if (natoms != atom->natoms && output->thermo->lostflag == ERROR && !region) 
    error->all(FLERR,"Atom count is inconsistent, cannot write restart file");

  // check if filename contains "%"

  int multiproc;
  if (strchr(file,'%')) multiproc = 1;
  else multiproc = 0;

  // open single restart file or base file for multiproc case

  if (me == 0) {
    char *hfile;
    if (multiproc) {
      hfile = new char[strlen(file) + 16];
      char *ptr = strchr(file,'%');
      *ptr = '\0';
      sprintf(hfile,"%s%s%s",file,"base",ptr+1);
      *ptr = '%';
    } else hfile = file;
    fp = fopen(hfile,"wb");
    if (fp == NULL) {
      char str[512];
      sprintf(str,"Cannot open restart file %s",hfile);
      error->one(FLERR,str);
    }
    if (multiproc) delete [] hfile;
  }

  // proc 0 writes header, groups, ntype-length arrays, force field
  // all procs write fix info

  if (me == 0) {
    header();
    group->write_restart(fp);
    type_arrays();
    force_fields();
  }

  modify->write_restart(fp);

  // communication buffer for all my atom's info
  // max_size = largest buffer needed by any proc

  int max_size;
  int send_size = atom->avec->size_restart();
  MPI_Allreduce(&send_size,&max_size,1,MPI_INT,MPI_MAX,world);

  double *buf;
  if (me == 0) memory->create(buf,max_size,"write_restart:buf");
  else memory->create(buf,send_size,"write_restart:buf");

  // pack my atom data into buf

  AtomVec *avec = atom->avec;
  int n = 0;

  if(!region) 
  {
      for (int i = 0; i < atom->nlocal; i++) n += avec->pack_restart(i,&buf[n]);
  }
  else
  {
      for (int i = 0; i < atom->nlocal; i++)
          if(region->match(atom->x[i][0],atom->x[i][1],atom->x[i][2]))
             n += avec->pack_restart(i,&buf[n]);

      send_size = n;
  }

  // if any fix requires it, remap each atom's coords via PBC
  // is because fix changes atom coords (excepting an integrate fix)
  // just remap in buffer, not actual atoms

  if (modify->restart_pbc_any) {
//.........这里部分代码省略.........

void DeleteAtoms::command(int narg, char **arg)
{
  if (domain->box_exist == 0)
    error->all(FLERR,"Delete_atoms command before simulation box is defined");
  if (narg < 1) error->all(FLERR,"Illegal delete_atoms command");
  if (atom->tag_enable == 0)
    error->all(FLERR,"Cannot use delete_atoms unless atoms have IDs");

  // store state before delete

  bigint natoms_previous = atom->natoms;
  bigint nbonds_previous = atom->nbonds;
  bigint nangles_previous = atom->nangles;
  bigint ndihedrals_previous = atom->ndihedrals;
  bigint nimpropers_previous = atom->nimpropers;

  // delete the atoms

  if (strcmp(arg[0],"group") == 0) delete_group(narg,arg);
  else if (strcmp(arg[0],"region") == 0) delete_region(narg,arg);
  else if (strcmp(arg[0],"overlap") == 0) delete_overlap(narg,arg);
  else if (strcmp(arg[0],"porosity") == 0) delete_porosity(narg,arg);
  else error->all(FLERR,"Illegal delete_atoms command");

  // optionally delete additional bonds or atoms in molecules

  if (bond_flag) delete_bond();
  if (mol_flag) delete_molecule();

  // delete local atoms flagged in dlist
  // reset nlocal

  AtomVec *avec = atom->avec;
  int nlocal = atom->nlocal;

  int i = 0;
  while (i < nlocal) {
    if (dlist[i]) {
      avec->copy(nlocal-1,i,1);
      dlist[i] = dlist[nlocal-1];
      nlocal--;
    } else i++;
  }

  atom->nlocal = nlocal;
  memory->destroy(dlist);

  // if non-molecular system and compress flag set,
  // reset atom tags to be contiguous
  // set all atom IDs to 0, call tag_extend()

  if (atom->molecular == 0 && compress_flag) {
    tagint *tag = atom->tag;
    for (i = 0; i < nlocal; i++) tag[i] = 0;
    atom->tag_extend();
  }

  // reset atom->natoms and also topology counts
  // reset atom->map if it exists
  // set nghost to 0 so old ghosts of deleted atoms won't be mapped

  bigint nblocal = atom->nlocal;
  MPI_Allreduce(&nblocal,&atom->natoms,1,MPI_LMP_BIGINT,MPI_SUM,world);
  if (atom->map_style) {
    atom->nghost = 0;
    atom->map_init();
    atom->map_set();
  }

  recount_topology();

  // print before and after atom and topology counts

  bigint ndelete = natoms_previous - atom->natoms;
  bigint ndelete_bonds = nbonds_previous - atom->nbonds;
  bigint ndelete_angles = nangles_previous - atom->nangles;
  bigint ndelete_dihedrals = ndihedrals_previous - atom->ndihedrals;
  bigint ndelete_impropers = nimpropers_previous - atom->nimpropers;

  if (comm->me == 0) {
    if (screen) {
      fprintf(screen,"Deleted " BIGINT_FORMAT
              " atoms, new total = " BIGINT_FORMAT "\n",
              ndelete,atom->natoms);
      if (bond_flag || mol_flag) {
        if (nbonds_previous) 
          fprintf(screen,"Deleted " BIGINT_FORMAT
                  " bonds, new total = " BIGINT_FORMAT "\n",
                  ndelete_bonds,atom->nbonds);
        if (nangles_previous) 
          fprintf(screen,"Deleted " BIGINT_FORMAT
                  " angles, new total = " BIGINT_FORMAT "\n",
                  ndelete_angles,atom->nangles);
        if (ndihedrals_previous) 
          fprintf(screen,"Deleted " BIGINT_FORMAT
                  " dihedrals, new total = " BIGINT_FORMAT "\n",
                  ndelete_dihedrals,atom->ndihedrals);
        if (nimpropers_previous) 
          fprintf(screen,"Deleted " BIGINT_FORMAT
                  " impropers, new total = " BIGINT_FORMAT "\n",
//.........这里部分代码省略.........

void FixSRP::post_run()
{
  // all bond particles are removed after each run
  // useful for write_data and write_restart commands
  // since those commands occur between runs

  bigint natoms_previous = atom->natoms;
  int nlocal = atom->nlocal;
  int* dlist;
  memory->create(dlist,nlocal,"fix_srp:dlist");

  for (int i = 0; i < nlocal; i++){
    if(atom->type[i] == bptype)
      dlist[i] = 1;
    else
      dlist[i] = 0;
  }

  // delete local atoms flagged in dlist
  // reset nlocal

  AtomVec *avec = atom->avec;

  int i = 0;
  while (i < nlocal) {
    if (dlist[i]) {
      avec->copy(nlocal-1,i,1);
      dlist[i] = dlist[nlocal-1];
      nlocal--;
    } else i++;
  }

  atom->nlocal = nlocal;
  memory->destroy(dlist);

  // reset atom->natoms
  // reset atom->map if it exists
  // set nghost to 0 so old ghosts won't be mapped

  bigint nblocal = atom->nlocal;
  MPI_Allreduce(&nblocal,&atom->natoms,1,MPI_LMP_BIGINT,MPI_SUM,world);
  if (atom->map_style) {
    atom->nghost = 0;
    atom->map_init();
    atom->map_set();
  }

  // print before and after atom count

  bigint ndelete = natoms_previous - atom->natoms;

  if (comm->me == 0) {
    if (screen) fprintf(screen,"Deleted " BIGINT_FORMAT
                        " atoms, new total = " BIGINT_FORMAT "\n",
                        ndelete,atom->natoms);
    if (logfile) fprintf(logfile,"Deleted " BIGINT_FORMAT
                         " atoms, new total = " BIGINT_FORMAT "\n",
                         ndelete,atom->natoms);
  }

  // verlet calls box_too_small_check() in post_run
  // this check maps all bond partners
  // therefore need ghosts

  // need to convert to lambda coords before apply pbc
  if (domain->triclinic) domain->x2lamda(atom->nlocal);
  domain->pbc();
  comm->setup();
  comm->exchange();
  if (atom->sortfreq > 0) atom->sort();
  comm->borders();
  // change back to box coordinates
  if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
}

void CommCuda::forward_comm_pack_cuda()
{
	static int count=0;
	static double kerneltime=0.0;
	static double copytime=0.0;
    timespec time1,time2,time3;
  int n;  // initialize comm buffers & exchange memory

  MPI_Request request;
  MPI_Status status;
  AtomVec *avec = atom->avec;
  double **x = atom->x;

  cuda->shared_data.domain.xy=domain->xy;
  cuda->shared_data.domain.xz=domain->xz;
  cuda->shared_data.domain.yz=domain->yz;
  cuda->shared_data.domain.prd[0]=domain->prd[0];
  cuda->shared_data.domain.prd[1]=domain->prd[1];
  cuda->shared_data.domain.prd[2]=domain->prd[2];
  cuda->shared_data.domain.triclinic=domain->triclinic;
  if(not comm_x_only && not avec->cudable) cuda->downloadAll();  //if not comm_x_only the communication routine of the atom_vec style class is used 

  // exchange data with another proc
  // if other proc is self, just copy
  // if comm_x_only set, exchange or copy directly to x, don't unpack

  for (int iswap = 0; iswap < nswap; iswap++) {
    if (sendproc[iswap] != me) 
    {   
      if (comm_x_only) 
      {

		
clock_gettime(CLOCK_REALTIME,&time1);
        
      //  n = Cuda_CommCuda_PackComm(&cuda->shared_data,sendnum[iswap],iswap,(void*) cuda->shared_data.comm.buf_send[iswap],pbc[iswap],pbc_flag[iswap]);
		  n = Cuda_CommCuda_PackComm(&cuda->shared_data,sendnum[iswap],iswap,(void*)buf_send,pbc[iswap],pbc_flag[iswap]);
		
clock_gettime(CLOCK_REALTIME,&time2);
		
        if((sizeof(X_FLOAT)!=sizeof(double)) && n) //some complicated way to safe some transfer size if single precision is used
          n=(n+1)*sizeof(X_FLOAT)/sizeof(double);
		cuda->shared_data.comm.send_size[iswap]=n;
      } 
      else if (ghost_velocity) 
      {
clock_gettime(CLOCK_REALTIME,&time1);

       // n = Cuda_CommCuda_PackComm_Vel(&cuda->shared_data,sendnum[iswap],iswap,(void*) &buf_send[iswap*maxsend],pbc[iswap],pbc_flag[iswap]);
		
clock_gettime(CLOCK_REALTIME,&time2);
		
        if((sizeof(X_FLOAT)!=sizeof(double)) && n) //some complicated way to safe some transfer size if single precision is used
          n=(n+1)*sizeof(X_FLOAT)/sizeof(double);
		cuda->shared_data.comm.send_size[iswap]=n;
       } 
      else 
      {
	    MPI_Irecv(buf_recv,size_forward_recv[iswap],MPI_DOUBLE,
		  recvproc[iswap],0,world,&request);
        
        if(avec->cudable)
          n = avec->pack_comm(sendnum[iswap],&iswap,
                           cuda->shared_data.comm.buf_send[iswap],pbc_flag[iswap],pbc[iswap]);
        else
	      n = avec->pack_comm(sendnum[iswap],sendlist[iswap],
			    cuda->shared_data.comm.buf_send[iswap],pbc_flag[iswap],pbc[iswap]);
	
	    MPI_Send(buf_send,n,MPI_DOUBLE,sendproc[iswap],0,world);
	    MPI_Wait(&request,&status);
	    avec->unpack_comm(recvnum[iswap],firstrecv[iswap],buf_recv);
      }

    }
    else  //sendproc == me
    {
      if (comm_x_only) 
      {     	
	    if (sendnum[iswap])
		{
          n = Cuda_CommCuda_PackComm_Self(&cuda->shared_data,sendnum[iswap],iswap,firstrecv[iswap],pbc[iswap],pbc_flag[iswap]);
          if(n<0) error->all(FLERR," # CUDA ERRROR on PackComm_Self");
          if((sizeof(X_FLOAT)!=sizeof(double)) && n)
            n=(n+1)*sizeof(X_FLOAT)/sizeof(double);
		}
      } 
      else if (ghost_velocity) 
      {
		n = avec->pack_comm_vel(sendnum[iswap],sendlist[iswap],
				buf_send,pbc_flag[iswap],pbc[iswap]);
	    avec->unpack_comm_vel(recvnum[iswap],firstrecv[iswap],buf_send);
      } 
      else 
      {
		n = avec->pack_comm(sendnum[iswap],sendlist[iswap],
			    buf_send,pbc_flag[iswap],pbc[iswap]);
		avec->unpack_comm(recvnum[iswap],firstrecv[iswap],buf_send);
      }
    }
  }
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	radsum = radtmp + xnear[i][3];
	if (rsq <= radsum*radsum) break;
      }
      if (i == nnear) {
	success = 1;
	break;
      }
    }
    if (success) {
      xnear[nnear][0] = coord[0];
      xnear[nnear][1] = coord[1];
      xnear[nnear][2] = coord[2];
      xnear[nnear][3] = radtmp;
      nnear++;
    } else break;
  }

  // warn if not all insertions were performed

  ninserted += nnear-nprevious;
  if (nnear - nprevious < nnew && me == 0)
    error->warning("Less insertions than requested");

  // 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 atom
  // type, diameter, density set from fix parameters
  // group mask set to "all" plus fix group
  // z velocity set to what velocity would be if particle
  //   had fallen from top of insertion region
  // this gives continuous stream of atoms
  // set npartner for new atom to 0 (assume not touching any others)

  AtomVec *avec = atom->avec;
  int m,flag;
  double denstmp,vxtmp,vytmp,vztmp;
  double g = 1.0;
  double *sublo = domain->sublo;
  double *subhi = domain->subhi;

  for (i = nprevious; i < nnear; i++) {
    coord[0] = xnear[i][0];
    coord[1] = xnear[i][1];
    coord[2] = xnear[i][2];
    radtmp = xnear[i][3];
    denstmp = density_lo + random->uniform() * (density_hi-density_lo);
    if (domain->dimension == 3) {
      vxtmp = vxlo + random->uniform() * (vxhi-vxlo);
      vytmp = vylo + random->uniform() * (vyhi-vylo);
      vztmp = vz - sqrt(2.0*g*(hi_current-coord[2]));
    } else {
      vxtmp = vxlo + random->uniform() * (vxhi-vxlo);
      vytmp = vy - sqrt(2.0*g*(hi_current-coord[1]));
      vztmp = 0.0;
    }

    flag = 0;
    if (coord[0] >= sublo[0] && coord[0] < subhi[0] &&
	coord[1] >= sublo[1] && coord[1] < subhi[1] &&
	coord[2] >= sublo[2] && coord[2] < subhi[2]) flag = 1;
    else if (domain->dimension == 3 && coord[2] >= domain->boxhi[2] &&
	     comm->myloc[2] == comm->procgrid[2]-1 &&
	     coord[0] >= sublo[0] && coord[0] < subhi[0] &&
	     coord[1] >= sublo[1] && coord[1] < subhi[1]) flag = 1;
    else if (domain->dimension == 2 && coord[1] >= domain->boxhi[1] &&
	     comm->myloc[1] == comm->procgrid[1]-1 &&

void Replicate::command(int narg, char **arg)
{
  int i,j,m,n;

  if (domain->box_exist == 0)
    error->all("Replicate command before simulation box is defined");
  if (narg != 3) error->all("Illegal replicate command");

  int me = comm->me;
  int nprocs = comm->nprocs;

  if (me == 0 && screen) fprintf(screen,"Replicating atoms ...\n");

  // nrep = total # of replications

  int nx = atoi(arg[0]);
  int ny = atoi(arg[1]);
  int nz = atoi(arg[2]);
  int nrep = nx*ny*nz;

  // error and warning checks

  if (nx <= 0 || ny <= 0 || nz <= 0) error->all("Illegal replicate command");
  if (domain->dimension == 2 && nz != 1)
    error->all("Cannot replicate 2d simulation in z dimension");
  if ((nx > 1 && domain->xperiodic == 0) || 
      (ny > 1 && domain->yperiodic == 0) ||
      (nz > 1 && domain->zperiodic == 0)) 
    error->warning("Replicating in a non-periodic dimension");

  if (atom->nextra_grow || atom->nextra_restart || atom->nextra_store)
    error->all("Cannot replicate with fixes that store atom quantities");

  // maxtag = largest atom tag across all existing atoms

  int maxtag = 0;
  for (i = 0; i < atom->nlocal; i++) maxtag = MAX(atom->tag[i],maxtag);
  int maxtag_all;
  MPI_Allreduce(&maxtag,&maxtag_all,1,MPI_INT,MPI_MAX,world);
  maxtag = maxtag_all;

  // maxmol = largest molecule tag across all existing atoms

  int maxmol = 0;
  if (atom->molecular) {
    for (i = 0; i < atom->nlocal; i++) maxmol = MAX(atom->molecule[i],maxmol);
    int maxmol_all;
    MPI_Allreduce(&maxmol,&maxmol_all,1,MPI_INT,MPI_MAX,world);
    maxmol = maxmol_all;
  }

  // unmap existing atoms via image flags

  for (i = 0; i < atom->nlocal; i++)
    domain->unmap(atom->x[i],atom->image[i]);

  // communication buffer for all my atom's info
  // max_size = largest buffer needed by any proc
  // must do before new Atom class created,
  //   since size_restart() uses atom->nlocal

  int max_size;
  int send_size = atom->avec->size_restart();
  MPI_Allreduce(&send_size,&max_size,1,MPI_INT,MPI_MAX,world);

  double *buf = 
    (double *) memory->smalloc(max_size*sizeof(double),"replicate:buf");

  // old = original atom class
  // atom = new replicated atom class
  // if old atom style was hybrid, pass sub-style names to create_avec

  Atom *old = atom;
  atom = new Atom(lmp);

  int nstyles = 0;
  char **keywords = NULL;
  if (strcmp(old->atom_style,"hybrid") == 0) {
    AtomVecHybrid *avec_hybrid = (AtomVecHybrid *) old->avec;
    nstyles = avec_hybrid->nstyles;
    keywords = avec_hybrid->keywords;
  }
  atom->create_avec(old->atom_style,nstyles,keywords);

  // check that new problem size will not be too large
  // if N > 2^31, turn off tags
  // if molecular, N/Nbonds/etc cannot be > 2^31 else tags/counts invalid

  double rep = nrep;
  if (rep*old->natoms > MAXATOMS) atom->tag_enable = 0;

  if (atom->molecular) {
    if (rep*old->natoms > MAXATOMS || rep*old->nbonds > MAXATOMS ||
	rep*old->nangles > MAXATOMS || rep*old->ndihedrals > MAXATOMS ||
	rep*old->nimpropers > MAXATOMS)
      error->all("Too big a problem to replicate with molecular atom style");
  }

  // assign atom and topology counts in new class from old one

//.........这里部分代码省略.........

void CommBrick::borders()
{
  int i,n,itype,iswap,dim,ineed,twoneed;
  int nsend,nrecv,sendflag,nfirst,nlast,ngroup;
  double lo,hi;
  int *type;
  double **x;
  double *buf,*mlo,*mhi;
  MPI_Request request;
  AtomVec *avec = atom->avec;

  // do swaps over all 3 dimensions

  iswap = 0;
  smax = rmax = 0;

  for (dim = 0; dim < 3; dim++) {
    nlast = 0;
    twoneed = 2*maxneed[dim];
    for (ineed = 0; ineed < twoneed; ineed++) {

      // find atoms within slab boundaries lo/hi using <= and >=
      // check atoms between nfirst and nlast
      //   for first swaps in a dim, check owned and ghost
      //   for later swaps in a dim, only check newly arrived ghosts
      // store sent atom indices in sendlist for use in future timesteps

      x = atom->x;
      if (mode == SINGLE) {
        lo = slablo[iswap];
        hi = slabhi[iswap];
      } else {
        type = atom->type;
        mlo = multilo[iswap];
        mhi = multihi[iswap];
      }
      if (ineed % 2 == 0) {
        nfirst = nlast;
        nlast = atom->nlocal + atom->nghost;
      }

      nsend = 0;

      // sendflag = 0 if I do not send on this swap
      // sendneed test indicates receiver no longer requires data
      // e.g. due to non-PBC or non-uniform sub-domains

      if (ineed/2 >= sendneed[dim][ineed % 2]) sendflag = 0;
      else sendflag = 1;

      // find send atoms according to SINGLE vs MULTI
      // all atoms eligible versus only atoms in bordergroup
      // can only limit loop to bordergroup for first sends (ineed < 2)
      // on these sends, break loop in two: owned (in group) and ghost

      if (sendflag) {
        if (!bordergroup || ineed >= 2) {
          if (mode == SINGLE) {
            for (i = nfirst; i < nlast; i++)
              if (x[i][dim] >= lo && x[i][dim] <= hi) {
                if (nsend == maxsendlist[iswap]) grow_list(iswap,nsend);
                sendlist[iswap][nsend++] = i;
              }
          } else {
            for (i = nfirst; i < nlast; i++) {
              itype = type[i];
              if (x[i][dim] >= mlo[itype] && x[i][dim] <= mhi[itype]) {
                if (nsend == maxsendlist[iswap]) grow_list(iswap,nsend);
                sendlist[iswap][nsend++] = i;
              }
            }
          }

        } else {
          if (mode == SINGLE) {
            ngroup = atom->nfirst;
            for (i = 0; i < ngroup; i++)
              if (x[i][dim] >= lo && x[i][dim] <= hi) {
                if (nsend == maxsendlist[iswap]) grow_list(iswap,nsend);
                sendlist[iswap][nsend++] = i;
              }
            for (i = atom->nlocal; i < nlast; i++)
              if (x[i][dim] >= lo && x[i][dim] <= hi) {
                if (nsend == maxsendlist[iswap]) grow_list(iswap,nsend);
                sendlist[iswap][nsend++] = i;
              }
          } else {
            ngroup = atom->nfirst;
            for (i = 0; i < ngroup; i++) {
              itype = type[i];
              if (x[i][dim] >= mlo[itype] && x[i][dim] <= mhi[itype]) {
                if (nsend == maxsendlist[iswap]) grow_list(iswap,nsend);
                sendlist[iswap][nsend++] = i;
              }
            }
            for (i = atom->nlocal; i < nlast; i++) {
              itype = type[i];
              if (x[i][dim] >= mlo[itype] && x[i][dim] <= mhi[itype]) {
                if (nsend == maxsendlist[iswap]) grow_list(iswap,nsend);
                sendlist[iswap][nsend++] = i;
//.........这里部分代码省略.........

void CommBrick::exchange()
{
  int i,m,nsend,nrecv,nrecv1,nrecv2,nlocal;
  double lo,hi,value;
  double **x;
  double *sublo,*subhi;
  MPI_Request request;
  AtomVec *avec = atom->avec;

  // clear global->local map for owned and ghost atoms
  // b/c atoms migrate to new procs in exchange() and
  //   new ghosts are created in borders()
  // map_set() is done at end of borders()
  // clear ghost count and any ghost bonus data internal to AtomVec

  if (map_style) atom->map_clear();
  atom->nghost = 0;
  atom->avec->clear_bonus();

  // insure send buf is large enough for single atom
  // bufextra = max size of one atom = allowed overflow of sendbuf
  // fixes can change per-atom size requirement on-the-fly

  int bufextra_old = bufextra;
  maxexchange = maxexchange_atom + maxexchange_fix;
  bufextra = maxexchange + BUFEXTRA;
  if (bufextra > bufextra_old)
    memory->grow(buf_send,maxsend+bufextra,"comm:buf_send");

  // subbox bounds for orthogonal or triclinic

  if (triclinic == 0) {
    sublo = domain->sublo;
    subhi = domain->subhi;
  } else {
    sublo = domain->sublo_lamda;
    subhi = domain->subhi_lamda;
  }

  // loop over dimensions

  int dimension = domain->dimension;

  for (int dim = 0; dim < dimension; dim++) {

    // fill buffer with atoms leaving my box, using < and >=
    // when atom is deleted, fill it in with last atom

    x = atom->x;
    lo = sublo[dim];
    hi = subhi[dim];
    nlocal = atom->nlocal;
    i = nsend = 0;

    while (i < nlocal) {
      if (x[i][dim] < lo || x[i][dim] >= hi) {
        if (nsend > maxsend) grow_send(nsend,1);
        nsend += avec->pack_exchange(i,&buf_send[nsend]);
        avec->copy(nlocal-1,i,1);
        nlocal--;
      } else i++;
    }
    atom->nlocal = nlocal;

    // send/recv atoms in both directions
    // send size of message first so receiver can realloc buf_recv if needed
    // if 1 proc in dimension, no send/recv
    //   set nrecv = 0 so buf_send atoms will be lost
    // if 2 procs in dimension, single send/recv
    // if more than 2 procs in dimension, send/recv to both neighbors

    if (procgrid[dim] == 1) nrecv = 0;
    else {
      MPI_Sendrecv(&nsend,1,MPI_INT,procneigh[dim][0],0,
                   &nrecv1,1,MPI_INT,procneigh[dim][1],0,world,MPI_STATUS_IGNORE);
      nrecv = nrecv1;
      if (procgrid[dim] > 2) {
        MPI_Sendrecv(&nsend,1,MPI_INT,procneigh[dim][1],0,
                     &nrecv2,1,MPI_INT,procneigh[dim][0],0,world,MPI_STATUS_IGNORE);
        nrecv += nrecv2;
      }
      if (nrecv > maxrecv) grow_recv(nrecv);

      MPI_Irecv(buf_recv,nrecv1,MPI_DOUBLE,procneigh[dim][1],0,
                world,&request);
      MPI_Send(buf_send,nsend,MPI_DOUBLE,procneigh[dim][0],0,world);
      MPI_Wait(&request,MPI_STATUS_IGNORE);

      if (procgrid[dim] > 2) {
        MPI_Irecv(&buf_recv[nrecv1],nrecv2,MPI_DOUBLE,procneigh[dim][0],0,
                  world,&request);
        MPI_Send(buf_send,nsend,MPI_DOUBLE,procneigh[dim][1],0,world);
        MPI_Wait(&request,MPI_STATUS_IGNORE);
      }
    }

    // check incoming atoms to see if they are in my box
    // if so, add to my list
    // box check is only for this dimension,
    //   atom may be passed to another proc in later dims
//.........这里部分代码省略.........

void CommBrick::forward_comm(int dummy)
{
  int n;
  MPI_Request request;
  AtomVec *avec = atom->avec;
  double **x = atom->x;
  double *buf;

  // exchange data with another proc
  // if other proc is self, just copy
  // if comm_x_only set, exchange or copy directly to x, don't unpack

  for (int iswap = 0; iswap < nswap; iswap++) {
    if (sendproc[iswap] != me) {
      if (comm_x_only) {
        if (size_forward_recv[iswap]) {
          if (size_forward_recv[iswap]) buf = x[firstrecv[iswap]];
          else buf = NULL;
          MPI_Irecv(buf,size_forward_recv[iswap],MPI_DOUBLE,
                    recvproc[iswap],0,world,&request);
        }
        n = avec->pack_comm(sendnum[iswap],sendlist[iswap],
                            buf_send,pbc_flag[iswap],pbc[iswap]);
        if (n) MPI_Send(buf_send,n,MPI_DOUBLE,sendproc[iswap],0,world);
        if (size_forward_recv[iswap]) MPI_Wait(&request,MPI_STATUS_IGNORE);
      } else if (ghost_velocity) {
        if (size_forward_recv[iswap])
          MPI_Irecv(buf_recv,size_forward_recv[iswap],MPI_DOUBLE,
                    recvproc[iswap],0,world,&request);
        n = avec->pack_comm_vel(sendnum[iswap],sendlist[iswap],
                                buf_send,pbc_flag[iswap],pbc[iswap]);
        if (n) MPI_Send(buf_send,n,MPI_DOUBLE,sendproc[iswap],0,world);
        if (size_forward_recv[iswap]) MPI_Wait(&request,MPI_STATUS_IGNORE);
        avec->unpack_comm_vel(recvnum[iswap],firstrecv[iswap],buf_recv);
      } else {
        if (size_forward_recv[iswap])
          MPI_Irecv(buf_recv,size_forward_recv[iswap],MPI_DOUBLE,
                    recvproc[iswap],0,world,&request);
        n = avec->pack_comm(sendnum[iswap],sendlist[iswap],
                            buf_send,pbc_flag[iswap],pbc[iswap]);
        if (n) MPI_Send(buf_send,n,MPI_DOUBLE,sendproc[iswap],0,world);
        if (size_forward_recv[iswap]) MPI_Wait(&request,MPI_STATUS_IGNORE);
        avec->unpack_comm(recvnum[iswap],firstrecv[iswap],buf_recv);
      }

    } else {
      if (comm_x_only) {
        if (sendnum[iswap])
          avec->pack_comm(sendnum[iswap],sendlist[iswap],
                          x[firstrecv[iswap]],pbc_flag[iswap],pbc[iswap]);
      } else if (ghost_velocity) {
        avec->pack_comm_vel(sendnum[iswap],sendlist[iswap],
                            buf_send,pbc_flag[iswap],pbc[iswap]);
        avec->unpack_comm_vel(recvnum[iswap],firstrecv[iswap],buf_send);
      } else {
        avec->pack_comm(sendnum[iswap],sendlist[iswap],
                        buf_send,pbc_flag[iswap],pbc[iswap]);
        avec->unpack_comm(recvnum[iswap],firstrecv[iswap],buf_send);
      }
    }
  }
}

void ReadRestart::command(int narg, char **arg)
{
  if (narg != 1) error->all(FLERR,"Illegal read_restart command");

  if (domain->box_exist)
    error->all(FLERR,"Cannot read_restart after simulation box is defined");

  MPI_Comm_rank(world,&me);
  MPI_Comm_size(world,&nprocs);

  // if filename contains "*", search dir for latest restart file

  char *file = new char[strlen(arg[0]) + 16];
  if (strchr(arg[0],'*')) {
    int n;
    if (me == 0) {
      file_search(arg[0],file);
      n = strlen(file) + 1;
    }
    MPI_Bcast(&n,1,MPI_INT,0,world);
    MPI_Bcast(file,n,MPI_CHAR,0,world);
  } else strcpy(file,arg[0]);

  // check if filename contains "%"

  int multiproc;
  if (strchr(file,'%')) multiproc = 1;
  else multiproc = 0;

  // open single restart file or base file for multiproc case
  // auto-detect whether byte swapping needs to be done as file is read

  if (me == 0) {
    if (screen) fprintf(screen,"Reading restart file ...\n");
    char *hfile;
    if (multiproc) {
      hfile = new char[strlen(file) + 16];
      char *ptr = strchr(file,'%');
      *ptr = '\0';
      sprintf(hfile,"%s%s%s",file,"base",ptr+1);
      *ptr = '%';
    } else hfile = file;
    fp = fopen(hfile,"rb");
    if (fp == NULL) {
      char str[128];
      sprintf(str,"Cannot open restart file %s",hfile);
      error->one(FLERR,str);
    }
    swapflag = autodetect(&fp,hfile);
    if (multiproc) delete [] hfile;
  }

  MPI_Bcast(&swapflag,1,MPI_INT,0,world);

  // read header info and create atom style and simulation box

  header();
  domain->box_exist = 1;

  // problem setup using info from header

  int n;
  if (nprocs == 1) n = static_cast<int> (atom->natoms);
  else n = static_cast<int> (LB_FACTOR * atom->natoms / nprocs);

  atom->allocate_type_arrays();
  atom->avec->grow(n);
  n = atom->nmax;

  domain->print_box("  ");
  domain->set_initial_box();
  domain->set_global_box();
  comm->set_proc_grid();
  domain->set_local_box();

  // read groups, ntype-length arrays, force field, fix info from file
  // nextra = max # of extra quantities stored with each atom

  group->read_restart(fp);
  type_arrays();
  force_fields();

  int nextra = modify->read_restart(fp);
  atom->nextra_store = nextra;
  memory->create(atom->extra,n,nextra,"atom:extra");

  // single file:
  // nprocs_file = # of chunks in file
  // proc 0 reads chunks one at a time and bcasts it to other procs
  // each proc unpacks the atoms, saving ones in it's sub-domain
  // check for atom in sub-domain differs for orthogonal vs triclinic box
  // close restart file when done

  AtomVec *avec = atom->avec;

  int maxbuf = 0;
  double *buf = NULL;
  int m;

  if (multiproc == 0) {
//.........这里部分代码省略.........

//.........这里部分代码省略.........
                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 (molecular == 2) {
            if (molatom[i] == 0) {
              index = molindex[i];
              ndeltopo[0] += onemols[index]->nbonds;
              ndeltopo[1] += onemols[index]->nangles;
              ndeltopo[2] += onemols[index]->ndihedrals;
              ndeltopo[3] += onemols[index]->nimpropers;
            }
          }

        } 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;
}