Molecular Dynamics Thread

[masatotanaka]

/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, [e-mail]

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 "string.h"
#include "dump_atom.h"
#include "domain.h"
#include "atom.h"
#include "update.h"
#include "group.h"
#include "error.h"

using namespace LAMMPS_NS;

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

DumpAtom::DumpAtom(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, arg)
{
if (narg != 5) error->all("Illegal dump atom command");

scale_flag = 1;
image_flag = 0;
format_default = NULL;

// one-time file open

if (multifile == 0) openfile();
}

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

void DumpAtom::init()
{
if (image_flag == 0) size_one = 5;
else size_one = 8;

// default format depends on image flags

delete [] format;
if (format_user) {
int n = strlen(format_user) + 2;
format = new char[n];
strcpy(format,format_user);
strcat(format,"\n");
} else {
char *str;
if (image_flag == 0) str = "%d %d %g %g %g";
else str = "%d %d %g %g %g %d %d %d";
int n = strlen(str) + 2;
format = new char[n];
strcpy(format,str);
strcat(format,"\n");
}

// setup function ptrs

if (binary) header_choice = &DumpAtom::header_binary;
else header_choice = &DumpAtom::header_item;

if (scale_flag == 1 && image_flag == 0)
pack_choice = &DumpAtom::pack_scale_noimage;
else if (scale_flag == 1 && image_flag == 1)
pack_choice = &DumpAtom::pack_scale_image;
else if (scale_flag == 0 && image_flag == 0)
pack_choice = &DumpAtom::pack_noscale_noimage;
else if (scale_flag == 0 && image_flag == 1)
pack_choice = &DumpAtom::pack_noscale_image;

if (binary) write_choice = &DumpAtom::write_binary;
else if (image_flag == 0) write_choice = &DumpAtom::write_noimage;
else if (image_flag == 1) write_choice = &DumpAtom::write_image;
}

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

int DumpAtom::modify_param(int narg, char **arg)
{
if (strcmp(arg[0],"scale") == 0) {
if (narg < 2) error->all("Illegal dump_modify command");
if (strcmp(arg[1],"yes") == 0) scale_flag = 1;
else if (strcmp(arg[1],"no") == 0) scale_flag = 0;
else error->all("Illegal dump_modify command");
return 2;
} else if (strcmp(arg[0],"image") == 0) {
if (narg < 2) error->all("Illegal dump_modify command");
if (strcmp(arg[1],"yes") == 0) image_flag = 1;
else if (strcmp(arg[1],"no") == 0) image_flag = 0;
else error->all("Illegal dump_modify command");
return 2;
}
return 0;
}

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

void DumpAtom::write_header(int ndump)
{
(this->*header_choice)(ndump);
}

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

int DumpAtom::count()
{
if (igroup == 0) return atom->nlocal;

int *mask = atom->mask;
int nlocal = atom->nlocal;

int m = 0;
for (int i = 0; i < nlocal; i++)
if (mask & groupbit) m++;
return m;
}

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

int DumpAtom::pack()
{
return (this->*pack_choice)();
}

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

void DumpAtom::write_data(int n, double *buf)
{
(this->*write_choice)(n,buf);
}

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

void DumpAtom::header_binary(int ndump)
{
fwrite(&update->ntimestep,sizeof(int),1,fp);
fwrite(&ndump,sizeof(int),1,fp);
fwrite(&domain->boxxlo,sizeof(double),1,fp);
fwrite(&domain->boxxhi,sizeof(double),1,fp);
fwrite(&domain->boxylo,sizeof(double),1,fp);
fwrite(&domain->boxyhi,sizeof(double),1,fp);
fwrite(&domain->boxzlo,sizeof(double),1,fp);
fwrite(&domain->boxzhi,sizeof(double),1,fp);
fwrite(&size_one,sizeof(int),1,fp);
if (multiproc) {
int one = 1;
fwrite(&one,sizeof(int),1,fp);
} else fwrite(&nprocs,sizeof(int),1,fp);
}

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

void DumpAtom::header_item(int ndump)
{
fprintf(fp,"ITEM: TIMESTEP\n");
fprintf(fp,"%d\n",update->ntimestep);
fprintf(fp,"ITEM: NUMBER OF ATOMS\n");
fprintf(fp,"%d\n",ndump);
fprintf(fp,"ITEM: BOX BOUNDS\n");
fprintf(fp,"%g %g\n",domain->boxxlo,domain->boxxhi);
fprintf(fp,"%g %g\n",domain->boxylo,domain->boxyhi);
fprintf(fp,"%g %g\n",domain->boxzlo,domain->boxzhi);
fprintf(fp,"ITEM: ATOMS\n");
}

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

int DumpAtom::pack_scale_image()
{
int *tag = atom->tag;
int *type = atom->type;
int *image = atom->image;
int *mask = atom->mask;
double **x = atom->x;
int nlocal = atom->nlocal;

double boxxlo = domain->boxxlo;
double boxylo = domain->boxylo;
double boxzlo = domain->boxzlo;
double invxprd = 1.0/domain->xprd;
double invyprd = 1.0/domain->yprd;
double invzprd = 1.0/domain->zprd;

int m = 0;
for (int i = 0; i < nlocal; i++)
if (mask & groupbit) {
buf[m++] = tag;
buf[m++] = type;
buf[m++] = (x[0] - boxxlo) * invxprd;
buf[m++] = (x[1] - boxylo) * invyprd;
buf[m++] = (x[2] - boxzlo) * invzprd;
buf[m++] = (image & 1023) - 512;
buf[m++] = (image >> 10 & 1023) - 512;
buf[m++] = (image >> 20) - 512;
}

return m;
}

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

int DumpAtom::pack_scale_noimage()
{
int *tag = atom->tag;
int *type = atom->type;
int *mask = atom->mask;
double **x = atom->x;
int nlocal = atom->nlocal;

double boxxlo = domain->boxxlo;
double boxylo = domain->boxylo;
double boxzlo = domain->boxzlo;
double invxprd = 1.0/domain->xprd;
double invyprd = 1.0/domain->yprd;
double invzprd = 1.0/domain->zprd;

int m = 0;
for (int i = 0; i < nlocal; i++)
if (mask & groupbit) {
buf[m++] = tag;
buf[m++] = type;
buf[m++] = (x[0] - boxxlo) * invxprd;
buf[m++] = (x[1] - boxylo) * invyprd;
buf[m++] = (x[2] - boxzlo) * invzprd;
}

return m;
}

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

int DumpAtom::pack_noscale_image()
{
int *tag = atom->tag;
int *type = atom->type;
int *image = atom->image;
int *mask = atom->mask;
double **x = atom->x;
int nlocal = atom->nlocal;

int m = 0;
for (int i = 0; i < nlocal; i++)
if (mask & groupbit) {
buf[m++] = tag;
buf[m++] = type;
buf[m++] = x[0];
buf[m++] = x[1];
buf[m++] = x[2];
buf[m++] = (image & 1023) - 512;
buf[m++] = (image >> 10 & 1023) - 512;
buf[m++] = (image >> 20) - 512;
}

return m;
}

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

int DumpAtom::pack_noscale_noimage()
{
int *tag = atom->tag;
int *type = atom->type;
int *mask = atom->mask;
double **x = atom->x;
int nlocal = atom->nlocal;

int m = 0;
for (int i = 0; i < nlocal; i++)
if (mask & groupbit) {
buf[m++] = tag;
buf[m++] = type;
buf[m++] = x[0];
buf[m++] = x[1];
buf[m++] = x[2];
}

return m;
}

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

void DumpAtom::write_binary(int n, double *buf)
{
n *= size_one;
fwrite(&n,sizeof(int),1,fp);
fwrite(buf,sizeof(double),n,fp);
}

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

void DumpAtom::write_image(int n, double *buf)
{
int m = 0;
for (int i = 0; i < n; i++) {
fprintf(fp,format,
static_cast<int> (buf[m]), static_cast<int> (buf[m+1]),
buf[m+2],buf[m+3],buf[m+4], static_cast<int> (buf[m+5]),
static_cast<int> (buf[m+6]), static_cast<int> (buf[m+7]));
m += size_one;
}
}

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

void DumpAtom::write_noimage(int n, double *buf)
{
int m = 0;
for (int i = 0; i < n; i++) {
fprintf(fp,format,
static_cast<int> (buf[m]), static_cast<int> (buf[m+1]),
buf[m+2],buf[m+3],buf[m+4]);
m += size_one;
}
}

[masatotanaka]

/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, [e-mail]

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 "mpi.h"
#include "stdlib.h"
#include "string.h"
#include "stdio.h"
#include "math.h"
#include "domain.h"
#include "atom.h"
#include "force.h"
#include "update.h"
#include "modify.h"
#include "fix.h"
#include "region.h"
#include "lattice.h"
#include "comm.h"
#include "memory.h"
#include "error.h"

#define RegionInclude
#include "style.h"
#undef RegionInclude

using namespace LAMMPS_NS;

#define BIG 1.0e20
#define SMALL 1.0e-4
#define DELTA 1
#define MIN(a,b) ((a) < (b) ? (a) : (b))
#define MAX(a,b) ((a) > (b) ? (a) : (b))

/* ----------------------------------------------------------------------
default is periodic
------------------------------------------------------------------------- */

Domain::Domain(LAMMPS *lmp) : Pointers(lmp)
{
box_exist = 0;

nonperiodic = 0;
xperiodic = yperiodic = zperiodic = 1;
periodicity[0] = xperiodic;
periodicity[1] = yperiodic;
periodicity[2] = zperiodic;

boundary[0][0] = boundary[0][1] = 0;
boundary[1][0] = boundary[1][1] = 0;
boundary[2][0] = boundary[2][1] = 0;

boxxlo = boxylo = boxzlo = -0.5;
boxxhi = boxyhi = boxzhi = 0.5;

lattice = NULL;
nregion = maxregion = 0;
regions = NULL;
}

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

Domain::~Domain()
{
delete lattice;
for (int i = 0; i < nregion; i++) delete regions;
memory->sfree(regions);
}

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

void Domain::init()
{
// set box_change if box dimensions ever change
// due to shrink-wrapping, fix nph, npt, volume/rescale, uniaxial

box_change = 0;
if (nonperiodic == 2) box_change = 1;
for (int i = 0; i < modify->nfix; i++) {
if (strcmp(modify->fix->style,"nph") == 0) box_change = 1;
if (strcmp(modify->fix->style,"npt") == 0) box_change = 1;
if (strcmp(modify->fix->style,"volume/rescale") == 0) box_change = 1;
if (strcmp(modify->fix->style,"uniaxial") == 0) box_change = 1;
}
}

/* ----------------------------------------------------------------------
set initial global box from boxlo/hi (already set by caller)
adjust for shrink-wrapped dims
------------------------------------------------------------------------- */

void Domain::set_initial_box()
{
if (boundary[0][0] == 2) boxxlo -= SMALL;
else if (boundary[0][0] == 3) minxlo = boxxlo;
if (boundary[0][1] == 2) boxxhi += SMALL;
else if (boundary[0][1] == 3) minxhi = boxxhi;

if (boundary[1][0] == 2) boxylo -= SMALL;
else if (boundary[1][0] == 3) minylo = boxylo;
if (boundary[1][1] == 2) boxyhi += SMALL;
else if (boundary[1][1] == 3) minyhi = boxyhi;

if (boundary[2][0] == 2) boxzlo -= SMALL;
else if (boundary[2][0] == 3) minzlo = boxzlo;
if (boundary[2][1] == 2) boxzhi += SMALL;
else if (boundary[2][1] == 3) minzhi = boxzhi;
}

/* ----------------------------------------------------------------------
set global prd, prd_half, prd[], boxlo/hi[] from boxlo/hi
------------------------------------------------------------------------- */

void Domain::set_global_box()
{
xprd = boxxhi - boxxlo;
yprd = boxyhi - boxylo;
zprd = boxzhi - boxzlo;

xprd_half = 0.5*xprd;
yprd_half = 0.5*yprd;
zprd_half = 0.5*zprd;

prd[0] = xprd; prd[1] = yprd; prd[2] = zprd;
boxlo[0] = boxxlo; boxlo[1] = boxylo; boxlo[2] = boxzlo;
boxhi[0] = boxxhi; boxhi[1] = boxyhi; boxhi[2] = boxzhi;
}

/* ----------------------------------------------------------------------
set local subxlo-subzhi, sublo/hi[] from global boxxlo-boxzhi and proc grid
for uppermost proc, insure subhi = boxhi (in case round-off occurs)
------------------------------------------------------------------------- */

void Domain::set_local_box()
{
subxlo = boxxlo + comm->myloc[0] * xprd / comm->procgrid[0];
if (comm->myloc[0] < comm->procgrid[0]-1)
subxhi = boxxlo + (comm->myloc[0]+1) * xprd / comm->procgrid[0];
else subxhi = boxxhi;

subylo = boxylo + comm->myloc[1] * yprd / comm->procgrid[1];
if (comm->myloc[1] < comm->procgrid[1]-1)
subyhi = boxylo + (comm->myloc[1]+1) * yprd / comm->procgrid[1];
else subyhi = boxyhi;

subzlo = boxzlo + comm->myloc[2] * zprd / comm->procgrid[2];
if (comm->myloc[2] < comm->procgrid[2]-1)
subzhi = boxzlo + (comm->myloc[2]+1) * zprd / comm->procgrid[2];
else subzhi = boxzhi;

sublo[0] = subxlo; sublo[1] = subylo; sublo[2] = subzlo;
subhi[0] = subxhi; subhi[1] = subyhi; subhi[2] = subzhi;
}

/* ----------------------------------------------------------------------
reset global & local boxes due to global box boundary changes
if shrink-wrapped, determine atom extent and reset boxlo/hi
set global & local boxes from new boxlo/hi values
------------------------------------------------------------------------- */

void Domain::reset_box()
{
if (nonperiodic == 2) {

// compute extent of atoms on this proc

double extent[3][2],all[3][2];

extent[2][0] = extent[1][0] = extent[0][0] = BIG;
extent[2][1] = extent[1][1] = extent[0][1] = -BIG;

double **x = atom->x;
int nlocal = atom->nlocal;

for (int i = 0; i < nlocal; i++) {
extent[0][0] = MIN(extent[0][0],x[0]);
extent[0][1] = MAX(extent[0][1],x[0]);
extent[1][0] = MIN(extent[1][0],x[1]);
extent[1][1] = MAX(extent[1][1],x[1]);
extent[2][0] = MIN(extent[2][0],x[2]);
extent[2][1] = MAX(extent[2][1],x[2]);
}

// compute extent across all procs
// flip sign of MIN to do it in one Allreduce MAX

extent[0][0] = -extent[0][0];
extent[1][0] = -extent[1][0];
extent[2][0] = -extent[2][0];

MPI_Allreduce(extent,all,6,MPI_DOUBLE,MPI_MAX,world);

// in shrink-wrapped dims, set box by atom extent

if (xperiodic == 0) {
if (boundary[0][0] == 2) boxxlo = -all[0][0] - SMALL;
else if (boundary[0][0] == 3) boxxlo = MIN(-all[0][0]-SMALL,minxlo);
if (boundary[0][1] == 2) boxxhi = all[0][1] + SMALL;
else if (boundary[0][1] == 3) boxxhi = MAX(all[0][1]+SMALL,minxhi);
if (boxxlo > boxxhi) error->all("Illegal simulation box");
}
if (yperiodic == 0) {
if (boundary[1][0] == 2) boxylo = -all[1][0] - SMALL;
else if (boundary[1][0] == 3) boxylo = MIN(-all[1][0]-SMALL,minylo);
if (boundary[1][1] == 2) boxyhi = all[1][1] + SMALL;
else if (boundary[1][1] == 3) boxyhi = MAX(all[1][1]+SMALL,minyhi);
if (boxylo > boxyhi) error->all("Illegal simulation box");
}
if (zperiodic == 0) {
if (boundary[2][0] == 2) boxzlo = -all[2][0] - SMALL;
else if (boundary[2][0] == 3) boxzlo = MIN(-all[2][0]-SMALL,minzlo);
if (boundary[2][1] == 2) boxzhi = all[2][1] + SMALL;
else if (boundary[2][1] == 3) boxzhi = MAX(all[2][1]+SMALL,minzhi);
if (boxzlo > boxzhi) error->all("Illegal simulation box");
}
}

set_global_box();
set_local_box();
}

/* ----------------------------------------------------------------------
enforce PBC and modify box image flags for each atom
called every reneighboring
resulting coord must satisfy lo <= coord < hi
MAX is important since coord - prd < lo can happen when coord = hi
image = 10 bits for each dimension
increment/decrement in wrap-around fashion
------------------------------------------------------------------------- */

void Domain::pbc()
{
int i,idim,otherdims;
int nlocal = atom->nlocal;
double **x = atom->x;
int *image = atom->image;

if (xperiodic) {
for (i = 0; i < nlocal; i++) {
if (x[0] < boxxlo) {
x[0] += xprd;
idim = image & 1023;
otherdims = image ^ idim;
idim--;
idim &= 1023;
image = otherdims | idim;
}
if (x[0] >= boxxhi) {
x[0] -= xprd;
x[0] = MAX(x[0],boxxlo);
idim = image & 1023;
otherdims = image ^ idim;
idim++;
idim &= 1023;
image = otherdims | idim;
}
}
}

if (yperiodic) {
for (i = 0; i < nlocal; i++) {
if (x[1] < boxylo) {
x[1] += yprd;
idim = (image >> 10) & 1023;
otherdims = image ^ (idim << 10);
idim--;
idim &= 1023;
image = otherdims | (idim << 10);
}
if (x[1] >= boxyhi) {
x[1] -= yprd;
x[1] = MAX(x[1],boxylo);
idim = (image >> 10) & 1023;
otherdims = image ^ (idim << 10);
idim++;
idim &= 1023;
image = otherdims | (idim << 10);
}
}
}

if (zperiodic) {
for (i = 0; i < nlocal; i++) {
if (x[2] < boxzlo) {
x[2] += zprd;
idim = image >> 20;
otherdims = image ^ (idim << 20);
idim--;
idim &= 1023;
image = otherdims | (idim << 20);
}
if (x[2] >= boxzhi) {
x[2] -= zprd;
x[2] = MAX(x[2],boxzlo);
idim = image >> 20;
otherdims = image ^ (idim << 20);
idim++;
idim &= 1023;
image = otherdims | (idim << 20);
}
}
}
}

/* ----------------------------------------------------------------------
minimum image convention
use 1/2 of box size as test
------------------------------------------------------------------------- */

void Domain::minimum_image(double *dx, double *dy, double *dz)
{
if (xperiodic) {
if (fabs(*dx) > xprd_half) {
if (*dx < 0.0) *dx += xprd;
else *dx -= xprd;
}
}
if (yperiodic) {
if (fabs(*dy) > yprd_half) {
if (*dy < 0.0) *dy += yprd;
else *dy -= yprd;
}
}
if (zperiodic) {
if (fabs(*dz) > zprd_half) {
if (*dz < 0.0) *dz += zprd;
else *dz -= zprd;
}
}
}

/* ----------------------------------------------------------------------
minimum image convention
use 1/2 of box size as test
------------------------------------------------------------------------- */

void Domain::minimum_image(double *x)
{
if (xperiodic) {
if (fabs(x[0]) > xprd_half) {
if (x[0] < 0.0) x[0] += xprd;
else x[0] -= xprd;
}
}
if (yperiodic) {
if (fabs(x[1]) > yprd_half) {
if (x[1] < 0.0) x[1] += yprd;
else x[1] -= yprd;
}
}
if (zperiodic) {
if (fabs(x[2]) > zprd_half) {
if (x[2] < 0.0) x[2] += zprd;
else x[2] -= zprd;
}
}
}

/* ----------------------------------------------------------------------
remap the point into the periodic box no matter how far away
adjust image accordingly
resulting coord must satisfy lo <= coord < hi
MAX is important since coord - prd < lo can happen when coord = hi
image = 10 bits for each dimension
increment/decrement in wrap-around fashion
------------------------------------------------------------------------- */

void Domain::remap(double &x, double &y, double &z, int &image)
{
if (xperiodic) {
while (x < boxxlo) {
x += xprd;
int idim = image & 1023;
int otherdims = image ^ idim;
idim--;
idim &= 1023;
image = otherdims | idim;
}
while (x >= boxxhi) {
x -= xprd;
int idim = image & 1023;
int otherdims = image ^ idim;
idim++;
idim &= 1023;
image = otherdims | idim;
}
x = MAX(x,boxxlo);
}

if (yperiodic) {
while (y < boxylo) {
y += yprd;
int idim = (image >> 10) & 1023;
int otherdims = image ^ (idim << 10);
idim--;
idim &= 1023;
image = otherdims | (idim << 10);
}
while (y >= boxyhi) {
y -= yprd;
int idim = (image >> 10) & 1023;
int otherdims = image ^ (idim << 10);
idim++;
idim &= 1023;
image = otherdims | (idim << 10);
}
y = MAX(y,boxylo);
}

if (zperiodic) {
while (z < boxzlo) {
z += zprd;
int idim = image >> 20;
int otherdims = image ^ (idim << 20);
idim--;
idim &= 1023;
image = otherdims | (idim << 20);
}
while (z >= boxzhi) {
z -= zprd;
int idim = image >> 20;
int otherdims = image ^ (idim << 20);
idim++;
idim &= 1023;
image = otherdims | (idim << 20);
}
z = MAX(z,boxzlo);
}
}

/* ----------------------------------------------------------------------
unmap the point via image flags
don't reset image flag
------------------------------------------------------------------------- */

void Domain::unmap(double &x, double &y, double &z, int image)
{
int xbox = (image & 1023) - 512;
int ybox = (image >> 10 & 1023) - 512;
int zbox = (image >> 20) - 512;

x = x + xbox*xprd;
y = y + ybox*yprd;
z = z + zbox*zprd;
}

/* ----------------------------------------------------------------------
create a lattice
delete it if style = none
------------------------------------------------------------------------- */

void Domain::set_lattice(int narg, char **arg)
{
delete lattice;
lattice = new Lattice(lmp,narg,arg);
if (lattice->style == 0) {
delete lattice;
lattice = NULL;
}
}

/* ----------------------------------------------------------------------
create a new region
------------------------------------------------------------------------- */

void Domain::add_region(int narg, char **arg)
{
if (narg < 2) error->all("Illegal region command");

// error checks

for (int iregion = 0; iregion < nregion; iregion++)
if (strcmp(arg[0],regions[iregion]->id) == 0)
error->all("Reuse of region ID");

// extend Region list if necessary

if (nregion == maxregion) {
maxregion += DELTA;
regions = (Region **)
memory->srealloc(regions,maxregion*sizeof(Region *),"domain:regions");
}

// create the Region

if (strcmp(arg[1],"none") == 0) error->all("Invalid region style");

#define RegionClass
#define RegionStyle(key,Class) \
else if (strcmp(arg[1],#key) == 0) \
regions[nregion] = new Class(lmp,narg,arg);
#include "style.h"
#undef RegionClass

else error->all("Invalid region style");

nregion++;
}

/* ----------------------------------------------------------------------
boundary settings from the input script
------------------------------------------------------------------------- */

void Domain::set_boundary(int narg, char **arg)
{
if (narg != 3) error->all("Illegal boundary command");

char c;
for (int idim = 0; idim < 3; idim++)
for (int iside = 0; iside < 2; iside++) {
if (iside == 0) c = arg[idim][0];
else if (iside == 1 && strlen(arg[idim]) == 1) c = arg[idim][0];
else c = arg[idim][1];

if (c == 'p') boundary[idim][iside] = 0;
else if (c == 'f') boundary[idim][iside] = 1;
else if (c == 's') boundary[idim][iside] = 2;
else if (c == 'm') boundary[idim][iside] = 3;
else error->all("Illegal boundary command");
}

for (int idim = 0; idim < 3; idim++)
if ((boundary[idim][0] == 0 && boundary[idim][1]) ||
(boundary[idim][0] && boundary[idim][1] == 0))
error->all("Both sides of boundary must be periodic");

if (boundary[0][0] == 0) xperiodic = 1;
else xperiodic = 0;
if (boundary[1][0] == 0) yperiodic = 1;
else yperiodic = 0;
if (boundary[2][0] == 0) zperiodic = 1;
else zperiodic = 0;

periodicity[0] = xperiodic;
periodicity[1] = yperiodic;
periodicity[2] = zperiodic;

nonperiodic = 0;
if (xperiodic == 0 || yperiodic == 0 || zperiodic == 0) {
nonperiodic = 1;
if (boundary[0][0] >= 2 || boundary[0][1] >= 2 ||
boundary[1][0] >= 2 || boundary[1][1] >= 2 ||
boundary[2][0] >= 2 || boundary[2][1] >= 2) nonperiodic = 2;
}
}

[masatotanaka]

/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, [e-mail]

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:
James Fischer (High Performance Technologies, Inc)
Vincent Natoli (Stone Ridge Technology)
David Richie (Stone Ridge Technology)
------------------------------------------------------------------------- */

#include "math.h"
#include "string.h"
#include "ctype.h"
#include "pair_hybrid.h"
#include "atom.h"
#include "force.h"
#include "pair.h"
#include "neighbor.h"
#include "update.h"
#include "comm.h"
#include "memory.h"
#include "error.h"

using namespace LAMMPS_NS;

#define NEIGHEXTRA 10000

#define MIN(a,b) ((a) < (b) ? (a) : (b))
#define MAX(a,b) ((a) > (b) ? (a) : (b))

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

PairHybrid::PairHybrid(LAMMPS *lmp) : Pair(lmp)
{
nstyles = 0;
}

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

PairHybrid::~PairHybrid()
{
if (nstyles) {
for (int m = 0; m < nstyles; m++) delete styles[m];
delete [] styles;
for (int m = 0; m < nstyles; m++) delete [] keywords[m];
delete [] keywords;
}

if (allocated) {
memory->destroy_2d_int_array(setflag);
memory->destroy_2d_int_array(map);
memory->destroy_2d_double_array(cutsq);

delete [] nnlist;
delete [] maxneigh;
for (int m = 0; m < nstyles; m++) memory->sfree(nlist[m]);
delete [] nlist;

delete [] nnlist_full;
delete [] maxneigh_full;
for (int m = 0; m < nstyles; m++) memory->sfree(nlist_full[m]);
delete [] nlist_full;

for (int m = 0; m < nstyles; m++) delete [] firstneigh[m];
delete [] firstneigh;
for (int m = 0; m < nstyles; m++) delete [] numneigh[m];
delete [] numneigh;

for (int m = 0; m < nstyles; m++) delete [] firstneigh_full[m];
delete [] firstneigh_full;
for (int m = 0; m < nstyles; m++) delete [] numneigh_full[m];
delete [] numneigh_full;
}
}

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

void PairHybrid::compute(int eflag, int vflag)
{
int i,j,k,m,n,jfull,nneigh;
int *neighs,*mapi;
double **f_original;

// save ptrs to original neighbor lists

int **firstneigh_original = neighbor->firstneigh;
int *numneigh_original = neighbor->numneigh;
int **firstneigh_full_original = neighbor->firstneigh_full;
int *numneigh_full_original = neighbor->numneigh_full;

// if this is re-neighbor step, create sub-style lists

if (neighbor->ago == 0) {
int *type = atom->type;
int nlocal = atom->nlocal;
int nall = atom->nlocal + atom->nghost;

// realloc per-atom per-style firstneigh/numneigh half/full if necessary

if (nlocal > maxlocal) {
maxlocal = nlocal;
if (neigh_half_every) {
for (m = 0; m < nstyles; m++) {
delete [] firstneigh[m];
delete [] numneigh[m];
}
for (m = 0; m < nstyles; m++) {
firstneigh[m] = new int*[maxlocal];
numneigh[m] = new int[maxlocal];
}
}
if (neigh_full_every) {
for (m = 0; m < nstyles; m++) {
delete [] firstneigh_full[m];
delete [] numneigh_full[m];
}
for (m = 0; m < nstyles; m++) {
firstneigh_full[m] = new int*[maxlocal];
numneigh_full[m] = new int[maxlocal];
}
}
}

// nnlist[] = length of each sub-style half list
// nnlist_full[] = length of each sub-style full list
// count from half and/or full list depending on what sub-styles use

for (m = 0; m < nstyles; m++) nnlist[m] = nnlist_full[m] = 0;

if (neigh_half_every && neigh_full_every) {
for (i = 0; i < nlocal; i++) {
mapi = map[type];
neighs = firstneigh_original;
nneigh = numneigh_original;
for (k = 0; k < nneigh; k++) {
j = neighs[k];
if (j >= nall) j %= nall;
m = mapi[type[j]];
if (styles[m] && styles[m]->neigh_half_every) nnlist[m]++;
}
neighs = firstneigh_full_original;
nneigh = numneigh_full_original;
for (k = 0; k < nneigh; k++) {
j = neighs[k];
if (j >= nall) j %= nall;
m = mapi[type[j]];
if (styles[m] && styles[m]->neigh_full_every) nnlist_full[m]++;
}
}

} else if (neigh_half_every) {
for (i = 0; i < nlocal; i++) {
mapi = map[type];
neighs = firstneigh_original;
nneigh = numneigh_original;
for (k = 0; k < nneigh; k++) {
j = neighs[k];
if (j >= nall) j %= nall;
nnlist[mapi[type[j]]]++;
}
}

} else if (neigh_full_every) {
for (i = 0; i < nlocal; i++) {
mapi = map[type];
neighs = firstneigh_full_original;
nneigh = numneigh_full_original;
for (k = 0; k < nneigh; k++) {
j = neighs[k];
if (j >= nall) j %= nall;
nnlist_full[mapi[type[j]]]++;
}
}
}

// realloc sub-style nlist and nlist_full if necessary

if (neigh_half_every) {
for (m = 0; m < nstyles; m++) {
if (nnlist[m] > maxneigh[m]) {
memory->sfree(nlist[m]);
maxneigh[m] = nnlist[m] + NEIGHEXTRA;
nlist[m] = (int *)
memory->smalloc(maxneigh[m]*sizeof(int),"pair_hybrid:nlist");
}
nnlist[m] = 0;
}
}
if (neigh_full_every) {
for (m = 0; m < nstyles; m++) {
if (nnlist_full[m] > maxneigh_full[m]) {
memory->sfree(nlist_full[m]);
maxneigh_full[m] = nnlist_full[m] + NEIGHEXTRA;
nlist_full[m] = (int *)
memory->smalloc(maxneigh_full[m]*sizeof(int),
"pair_hybrid:nlist_full");
}
nnlist_full[m] = 0;
}
}

// load sub-style half/full list with values from original lists
// load from half and/or full list depending on what sub-styles use

if (neigh_half_every && neigh_full_every) {
for (i = 0; i < nlocal; i++) {
for (m = 0; m < nstyles; m++) {
firstneigh[m] = &nlist[m][nnlist[m]];
numneigh[m] = nnlist[m];
firstneigh_full[m] = &nlist_full[m][nnlist_full[m]];
numneigh_full[m] = nnlist_full[m];
}
mapi = map[type];
neighs = firstneigh_original;
nneigh = numneigh_original;
for (k = 0; k < nneigh; k++) {
j = jfull = neighs[k];
if (j >= nall) j %= nall;
m = mapi[type[j]];
if (styles[m] && styles[m]->neigh_half_every)
nlist[m][nnlist[m]++] = jfull;
}
neighs = firstneigh_full_original;
nneigh = numneigh_full_original;
for (k = 0; k < nneigh; k++) {
j = jfull = neighs[k];
if (j >= nall) j %= nall;
m = mapi[type[j]];
if (styles[m] && styles[m]->neigh_full_every)
nlist_full[m][nnlist_full[m]++] = jfull;
}
for (m = 0; m < nstyles; m++) {
numneigh[m] = nnlist[m] - numneigh[m];
numneigh_full[m] = nnlist_full[m] - numneigh_full[m];
}
}

} else if (neigh_half_every) {
for (i = 0; i < nlocal; i++) {
for (m = 0; m < nstyles; m++) {
firstneigh[m] = &nlist[m][nnlist[m]];
numneigh[m] = nnlist[m];
}
mapi = map[type];
neighs = firstneigh_original;
nneigh = numneigh_original;
for (k = 0; k < nneigh; k++) {
j = jfull = neighs[k];
if (j >= nall) j %= nall;
m = mapi[type[j]];
nlist[m][nnlist[m]++] = jfull;
}
for (m = 0; m < nstyles; m++)
numneigh[m] = nnlist[m] - numneigh[m];
}

} else if (neigh_full_every) {
for (i = 0; i < nlocal; i++) {
for (m = 0; m < nstyles; m++) {
firstneigh_full[m] = &nlist_full[m][nnlist_full[m]];
numneigh_full[m] = nnlist_full[m];
}
mapi = map[type];
neighs = firstneigh_full_original;
nneigh = numneigh_full_original;
for (k = 0; k < nneigh; k++) {
j = jfull = neighs[k];
if (j >= nall) j %= nall;
m = mapi[type[j]];
nlist_full[m][nnlist_full[m]++] = jfull;
}
for (m = 0; m < nstyles; m++)
numneigh_full[m] = nnlist_full[m] - numneigh_full[m];
}
}
}

// call each sub-style's compute function
// set neighbor->firstneigh/numneigh to sub-style lists before call
// set half or full or both depending on what sub-style uses
// for vflag = 1:
// sub-style accumulates in its virial[6]
// sum sub-style virial[6] to hybrid's virial[6]
// for vflag = 2:
// set atom->f to update->f_pair so sub-style will sum its f to f_pair
// call sub-style compute() with vflag % 2 to prevent sub-style
// from calling virial_compute()
// reset atom->f to stored f_original
// call hybrid virial_compute() which will use update->f_pair
// accumulate sub-style energy,virial in hybrid's energy,virial

eng_vdwl = eng_coul = 0.0;
if (vflag) for (n = 0; n < 6; n++) virial[n] = 0.0;

if (vflag == 2) {
f_original = atom->f;
atom->f = update->f_pair;
}

for (m = 0; m < nstyles; m++) {
if (styles[m] == NULL) continue;
if (styles[m]->neigh_half_every) {
neighbor->firstneigh = firstneigh[m];
neighbor->numneigh = numneigh[m];
}
if (styles[m]->neigh_full_every) {
neighbor->firstneigh_full = firstneigh_full[m];
neighbor->numneigh_full = numneigh_full[m];
}
styles[m]->compute(eflag,vflag % 2);
if (eflag) {
eng_vdwl += styles[m]->eng_vdwl;
eng_coul += styles[m]->eng_coul;
}
if (vflag == 1) for (n = 0; n < 6; n++) virial[n] += styles[m]->virial[n];
}

if (vflag == 2) {
atom->f = f_original;
virial_compute();
}

// restore ptrs to original neighbor lists

neighbor->firstneigh = firstneigh_original;
neighbor->numneigh = numneigh_original;
neighbor->firstneigh_full = firstneigh_full_original;
neighbor->numneigh_full = numneigh_full_original;
}

/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */

void PairHybrid::allocate()
{
allocated = 1;
int n = atom->ntypes;

map = memory->create_2d_int_array(n+1,n+1,"pair:map");
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
map[j] = -1;

setflag = memory->create_2d_int_array(n+1,n+1,"pair:setflag");
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[j] = 0;

cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq");

nnlist = new int[nstyles];
maxneigh = new int[nstyles];
nlist = new int*[nstyles];
for (int m = 0; m < nstyles; m++) maxneigh[m] = 0;
for (int m = 0; m < nstyles; m++) nlist[m] = NULL;

nnlist_full = new int[nstyles];
maxneigh_full = new int[nstyles];
nlist_full = new int*[nstyles];
for (int m = 0; m < nstyles; m++) maxneigh_full[m] = 0;
for (int m = 0; m < nstyles; m++) nlist_full[m] = NULL;

maxlocal = 0;
firstneigh = new int**[nstyles];
numneigh = new int*[nstyles];
for (int m = 0; m < nstyles; m++) firstneigh[m] = NULL;
for (int m = 0; m < nstyles; m++) numneigh[m] = NULL;
firstneigh_full = new int**[nstyles];
numneigh_full = new int*[nstyles];
for (int m = 0; m < nstyles; m++) firstneigh_full[m] = NULL;
for (int m = 0; m < nstyles; m++) numneigh_full[m] = NULL;
}

/* ----------------------------------------------------------------------
create one pair style for each arg in list
------------------------------------------------------------------------- */

void PairHybrid::settings(int narg, char **arg)
{
int i,m,istyle;

if (narg < 1) error->all("Illegal pair_style command");

// delete old lists, since cannot just change settings

if (nstyles) {
for (m = 0; m < nstyles; m++) delete styles[m];
delete [] styles;
for (m = 0; m < nstyles; m++) delete [] keywords[m];
delete [] keywords;
}

if (allocated) {
memory->destroy_2d_int_array(setflag);
memory->destroy_2d_int_array(map);
memory->destroy_2d_double_array(cutsq);
}
allocated = 0;

// count sub-styles by skipping numeric args
// one exception is 1st arg of style "table", which is non-numeric word

nstyles = i = 0;
while (i < narg) {
if (strcmp(arg,"table") == 0) i++;
i++;
while (i < narg && !isalpha(arg[0])) i++;
nstyles++;
}

// allocate list of sub-styles

styles = new Pair*[nstyles];
keywords = new char*[nstyles];

// allocate each sub-style and call its settings() with subset of args
// define subset of sub-styles by skipping numeric args
// one exception is 1st arg of style "table", which is non-numeric word

nstyles = i = 0;
while (i < narg) {
for (m = 0; m < nstyles; m++)
if (strcmp(arg,keywords[m]) == 0)
error->all("Pair style hybrid cannot use same pair style twice");
if (strcmp(arg,"hybrid") == 0)
error->all("Pair style hybrid cannot have hybrid as an argument");
styles[nstyles] = force->new_pair(arg);
keywords[nstyles] = new char[strlen(arg)+1];
strcpy(keywords[nstyles],arg);
istyle = i;
if (strcmp(arg,"table") == 0) i++;
i++;
while (i < narg && !isalpha(arg[0])) i++;
if (styles[nstyles]) styles[nstyles]->settings(i-istyle-1,&arg[istyle+1]);
nstyles++;
}

// set comm_forward, comm_reverse to max of any sub-style

for (m = 0; m < nstyles; m++) {
if (styles[m]) comm_forward = MAX(comm_forward,styles[m]->comm_forward);
if (styles[m]) comm_reverse = MAX(comm_reverse,styles[m]->comm_reverse);
}

// neigh_every = 1 if any sub-style = 1

neigh_half_every = neigh_full_every = 0;
for (m = 0; m < nstyles; m++) {
if (styles[m] && styles[m]->neigh_half_every) neigh_half_every = 1;
if (styles[m] && styles[m]->neigh_full_every) neigh_full_every = 1;
}

// single_enable = 0 if any sub-style = 0

for (m = 0; m < nstyles; m++)
if (styles[m] && styles[m]->single_enable == 0) single_enable = 0;
}

/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */

void PairHybrid::coeff(int narg, char **arg)
{
if (narg < 3) error->all("Incorrect args for pair coefficients");
if (!allocated) allocate();

int ilo,ihi,jlo,jhi;
force->bounds(arg[0],atom->ntypes,ilo,ihi);
force->bounds(arg[1],atom->ntypes,jlo,jhi);

// 3rd arg = pair style name

int m;
for (m = 0; m < nstyles; m++)
if (strcmp(arg[2],keywords[m]) == 0) break;
if (m == nstyles) error->all("Pair coeff for hybrid has invalid style");

// move 1st/2nd args to 2nd/3rd args

sprintf(arg[2],"%s",arg[1]);
sprintf(arg[1],"%s",arg[0]);

// invoke sub-style coeff() starting with 1st arg

if (styles[m]) styles[m]->coeff(narg-1,&arg[1]);

// if pair style only allows one pair coeff call (with * * and type mapping)
// then unset any old setflag/map assigned to that style first
// in case pair coeff for this sub-style is being called for 2nd time

if (styles[m] && styles[m]->one_coeff)
for (int i = 1; i <= atom->ntypes; i++)
for (int j = i; j <= atom->ntypes; j++)
if (map[j] == m) {
map[j] = -1;
setflag[j] = 0;
}

// set hybrid map & setflag only if substyle set its setflag
// if sub-style is NULL for "none", still set setflag

int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
if (styles[m] == NULL || styles[m]->setflag[j]) {
map[j] = m;
setflag[j] = 1;
count++;
}
}
}

if (count == 0) error->all("Incorrect args for pair coefficients");
}

/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */

double PairHybrid::init_one(int i, int j)
{
// if i,j is set explicity, call its sub-style
// if i,j is not set and i,i sub-style = j,j sub-style
// then set map[j] to this sub-style and call sub-style for init/mixing
// else i,j has not been set by user
// check for special case = style none

double cut = 0.0;
if (setflag[j]) {
if (styles[map[j]]) {
cut = styles[map[j]]->init_one(i,j);
styles[map[j]]->cutsq[j] =
styles[map[j]]->cutsq[j] = cut*cut;
if (tail_flag) {
etail_ij = styles[map[j]]->etail_ij;
ptail_ij = styles[map[j]]->ptail_ij;
}
}
} else if (map == map[j][j]) {
map[j] = map;
if (styles[map[j]]) {
cut = styles[map[j]]->init_one(i,j);
styles[map[j]]->cutsq[j] =
styles[map[j]]->cutsq[j] = cut*cut;
if (tail_flag) {
etail_ij = styles[map[j]]->etail_ij;
ptail_ij = styles[map[j]]->ptail_ij;
}
}
} else error->one("All pair coeffs are not set");

map[j] = map[j];

return cut;
}

/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */

void PairHybrid::init_style()
{
for (int m = 0; m < nstyles; m++)
if (styles[m]) styles[m]->init_style();
}

/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */

void PairHybrid::write_restart(FILE *fp)
{
fwrite(&nstyles,sizeof(int),1,fp);

// each sub-style writes its settings

int n;
for (int m = 0; m < nstyles; m++) {
n = strlen(keywords[m]) + 1;
fwrite(&n,sizeof(int),1,fp);
fwrite(keywords[m],sizeof(char),n,fp);
if (styles[m]) styles[m]->write_restart_settings(fp);
}
}

/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */

void PairHybrid::read_restart(FILE *fp)
{
allocate();

int me = comm->me;
if (me == 0) fread(&nstyles,sizeof(int),1,fp);
MPI_Bcast(&nstyles,1,MPI_INT,0,world);

styles = new Pair*[nstyles];
keywords = new char*[nstyles];

// each sub-style is created via new_pair() and reads its settings

int n;
for (int m = 0; m < nstyles; m++) {
if (me == 0) fread(&n,sizeof(int),1,fp);
MPI_Bcast(&n,1,MPI_INT,0,world);
keywords[m] = new char[n];
if (me == 0) fread(keywords[m],sizeof(char),n,fp);
MPI_Bcast(keywords[m],n,MPI_CHAR,0,world);
styles[m] = force->new_pair(keywords[m]);
if (styles[m]) styles[m]->read_restart_settings(fp);
}
}

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

void PairHybrid::single(int i, int j, int itype, int jtype,
double rsq, double factor_coul, double factor_lj,
int eflag, One &one)
{
if (map[itype][jtype] == -1)
error->one("Invoked pair single on pair style none");

styles[map[itype][jtype]]->
single(i,j,itype,jtype,rsq,factor_coul,factor_lj,eflag,one);
}

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

void PairHybrid::single_embed(int i, int itype, double &phi)
{
if (map[itype][itype] == -1)
error->one("Invoked pair single on pair style none");

styles[map[itype][itype]]->single_embed(i,itype,phi);
}

/* ----------------------------------------------------------------------
modify parameters of the pair style
simply pass command args to each sub-style of hybrid
------------------------------------------------------------------------- */

void PairHybrid::modify_params(int narg, char **arg)
{
for (int m = 0; m < nstyles; m++)
if (styles[m]) styles[m]->modify_params(narg,arg);
}

/* ----------------------------------------------------------------------
memory usage of sub-style firstneigh, numneigh, neighbor list
add in memory usage of each sub-style itself
------------------------------------------------------------------------- */

int PairHybrid::memory_usage()
{
int bytes = nstyles*maxlocal * (sizeof(int *) + sizeof(int));
for (int m = 0; m < nstyles; m++) bytes += maxneigh[m] * sizeof(int);
for (int m = 0; m < nstyles; m++)
if (styles[m]) bytes += styles[m]->memory_usage();
return bytes;
}

[masatotanaka]

/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, [e-mail]

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 "mpi.h"
#include "string.h"
#include "stdlib.h"
#include "sys/types.h"
#include "dirent.h"
#include "read_restart.h"
#include "atom.h"
#include "atom_vec.h"
#include "domain.h"
#include "comm.h"
#include "update.h"
#include "modify.h"
#include "group.h"
#include "force.h"
#include "pair.h"
#include "bond.h"
#include "angle.h"
#include "dihedral.h"
#include "improper.h"
#include "special.h"
#include "universe.h"
#include "memory.h"
#include "error.h"

using namespace LAMMPS_NS;

#define LB_FACTOR 1.1

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

ReadRestart::ReadRestart(LAMMPS *lmp) : Pointers(lmp) {}

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

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

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

MPI_Comm_rank(world,&me);

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

char *file = new char[strlen(arg[0]) + 16];
if (strchr(arg[0],'*')) file_search(arg[0],file);
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

if (me == 0) {
if (screen) fprintf(screen,"Reading restart file ...\n");
char *hfile;
if (multiproc) {
char *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(str);
}
if (multiproc) delete [] hfile;
}

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

header();
domain->box_exist = 1;

// problem setup using info from header

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

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

domain->set_initial_box();
domain->set_global_box();
comm->set_procs();
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);
if (atom->mass) mass();
if (atom->dipole) dipole();
force_fields();

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

// if single file:
// proc 0 reads atoms from file, one chunk per proc (nprocs_file)
// else if one file per proc:
// proc 0 reads chunks from series of files (nprocs_file)
// proc 0 bcasts each chunk to other procs
// each proc unpacks the atoms, saving ones in it's sub-domain

AtomVec *avec = atom->avec;

int maxbuf = 0;
double *buf = NULL;

double subxlo = domain->subxlo;
double subxhi = domain->subxhi;
double subylo = domain->subylo;
double subyhi = domain->subyhi;
double subzlo = domain->subzlo;
double subzhi = domain->subzhi;

int m;
double xtmp,ytmp,ztmp;
char *perproc = new char[strlen(file) + 16];
char *ptr = strchr(file,'%');

for (int iproc = 0; iproc < nprocs_file; iproc++) {
if (me == 0) {
if (multiproc) {
fclose(fp);
*ptr = '\0';
sprintf(perproc,"%s%d%s",file,iproc,ptr+1);
*ptr = '%';
fp = fopen(perproc,"rb");
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open restart file %s",perproc);
error->one(str);
}
}
fread(&n,sizeof(int),1,fp);
}

MPI_Bcast(&n,1,MPI_INT,0,world);
if (n > maxbuf) {
maxbuf = n;
delete [] buf;
buf = new double[maxbuf];
}

if (me == 0) fread(buf,sizeof(double),n,fp);
MPI_Bcast(buf,n,MPI_DOUBLE,0,world);

m = 0;
while (m < n) {
xtmp = buf[m+1];
ytmp = buf[m+2];
ztmp = buf[m+3];
if (xtmp >= subxlo && xtmp < subxhi &&
ytmp >= subylo && ytmp < subyhi &&
ztmp >= subzlo && ztmp < subzhi)
m += avec->unpack_restart(&buf[m]);
else m += static_cast<int> (buf[m]);
}
}

// close restart file and clean-up memory

if (me == 0) fclose(fp);
delete [] buf;
delete [] file;
delete [] perproc;

// check that all atoms were assigned to procs

double natoms;
double rlocal = atom->nlocal;
MPI_Allreduce(&rlocal,&natoms,1,MPI_DOUBLE,MPI_SUM,world);

if (me == 0) {
if (screen) fprintf(screen," %.15g atoms\n",natoms);
if (logfile) fprintf(logfile," %.15g atoms\n",natoms);
}

if (natoms != atom->natoms) error->all("Did not assign all atoms correctly");

if (me == 0) {
if (atom->nbonds) {
if (screen) fprintf(screen," %d bonds\n",atom->nbonds);
if (logfile) fprintf(logfile," %d bonds\n",atom->nbonds);
}
if (atom->nangles) {
if (screen) fprintf(screen," %d angles\n",atom->nangles);
if (logfile) fprintf(logfile," %d angles\n",atom->nangles);
}
if (atom->ndihedrals) {
if (screen) fprintf(screen," %d dihedrals\n",atom->ndihedrals);
if (logfile) fprintf(logfile," %d dihedrals\n",atom->ndihedrals);
}
if (atom->nimpropers) {
if (screen) fprintf(screen," %d impropers\n",atom->nimpropers);
if (logfile) fprintf(logfile," %d impropers\n",atom->nimpropers);
}
}

// check if tags are being used
// create global mapping and bond topology now that system is defined

int flag = 0;
for (int i = 0; i < atom->nlocal; i++)
if (atom->tag > 0) flag = 1;
int flag_all;
MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_MAX,world);
if (flag_all == 0) atom->tag_enable = 0;

if (atom->map_style) {
atom->map_init();
atom->map_set();
}
if (atom->molecular) {
Special special(lmp);
special.build();
}
}

/* ----------------------------------------------------------------------
search for all files in dir matching infile which contains "*"
replace "*" with latest timestep value to create outfile name
if infile also contains "%", need to use "base" when search directory
------------------------------------------------------------------------- */

void ReadRestart::file_search(char *infile, char *outfile)
{
// if filename contains "%" replace "%" with "base"

char *pattern = new char[strlen(infile) + 16];
char *ptr;

if (ptr = strchr(infile,'%')) {
*ptr = '\0';
sprintf(pattern,"%s%s%s",infile,"base",ptr+1);
*ptr = '%';
} else strcpy(pattern,infile);

// scan all files in directory, searching for files that match pattern
// maxnum = largest int that matches "*"

int n = strlen(pattern) + 16;
char *begin = new char[n];
char *middle = new char[n];
char *end = new char[n];

ptr = strchr(pattern,'*');
*ptr = '\0';
strcpy(begin,pattern);
strcpy(end,ptr+1);
int nbegin = strlen(begin);
int maxnum = -1;

if (me == 0) {
struct dirent *ep;
DIR *dp = opendir("./");
if (dp == NULL)
error->one("Cannot open dir to search for restart file");
while (ep = readdir(dp)) {
if (strstr(ep->d_name,begin) != ep->d_name) continue;
if ((ptr = strstr(&ep->d_name[nbegin],end)) == NULL) continue;
if (strlen(end) == 0) ptr = ep->d_name + strlen(ep->d_name);
*ptr = '\0';
if (strlen(&ep->d_name[nbegin]) < n) {
strcpy(middle,&ep->d_name[nbegin]);
if (atoi(middle) > maxnum) maxnum = atoi(middle);
}
}
closedir(dp);
if (maxnum < 0) error->one("Found no restart file matching pattern");
}

delete [] pattern;
delete [] begin;
delete [] middle;
delete [] end;

// create outfile with maxint substituted for "*"
// use original infile, not pattern, since need to retain "%" in filename

ptr = strchr(infile,'*');
*ptr = '\0';
sprintf(outfile,"%s%d%s",infile,maxnum,ptr+1);
*ptr = '*';
}

/* ----------------------------------------------------------------------
read header of restart file
------------------------------------------------------------------------- */

void ReadRestart::header()
{
int px,py,pz;
int xperiodic,yperiodic,zperiodic;
int boundary[3][2];

// read flags and values until flag = -1

int flag = read_int();
while (flag >= 0) {

// check restart file version, compare to LAMMPS version

if (flag == 0) {
char *version = read_char();
if (strcmp(version,universe->version) != 0) {
error->warning("Restart file version does not match LAMMPS version");
if (screen) fprintf(screen," restart file = %s, LAMMPS = %s\n",
version,universe->version);
}
delete [] version;

// reset unit_style only if different
// so that timestep,neighbor-skin are not changed

} else if (flag == 1) {
char *style = read_char();
if (strcmp(style,update->unit_style) != 0 && me == 0)
error->warning("Resetting unit_style to restart file value");
if (strcmp(style,update->unit_style) != 0) update->set_units(style);
delete [] style;

} else if (flag == 2) {
update->ntimestep = read_int();

// set dimension from restart file, warn if different

} else if (flag == 3) {
int dimension = read_int();
if (dimension != force->dimension && me == 0)
error->warning("Resetting dimension to restart file value");
force->dimension = dimension;
if (force->dimension == 2 && domain->zperiodic == 0)
error->all("Cannot run 2d simulation with nonperiodic Z dimension");

// read nprocs_file from restart file, warn if different

} else if (flag == 4) {
nprocs_file = read_int();
if (nprocs_file != comm->nprocs && me == 0)
error->warning("Restart file used different # of processors");

// don't set procgrid, warn if different

} else if (flag == 5) {
px = read_int();
} else if (flag == 6) {
py = read_int();
} else if (flag == 7) {
pz = read_int();
if (comm->user_procgrid[0] != 0 &&
(px != comm->user_procgrid[0] || py != comm->user_procgrid[1] ||
pz != comm->user_procgrid[2]) && me == 0)
error->warning("Restart file used different 3d processor grid");

// don't set newton_pair, warn if different
// set newton_bond from restart file, warn if different

} else if (flag == 8) {
int newton_pair = read_int();
if (newton_pair != force->newton_pair && me == 0)
error->warning("Restart file used different newton pair setting");
} else if (flag == 9) {
int newton_bond = read_int();
if (newton_bond != force->newton_bond && me == 0)
error->warning("Resetting newton bond to restart file value");
force->newton_bond = newton_bond;
if (force->newton_pair || force->newton_bond) force->newton = 1;
else force->newton = 0;

// set boundary settings from restart file, warn if different

} else if (flag == 10) {
xperiodic = read_int();
} else if (flag == 11) {
yperiodic = read_int();
} else if (flag == 12) {
zperiodic = read_int();
} else if (flag == 13) {
boundary[0][0] = read_int();
} else if (flag == 14) {
boundary[0][1] = read_int();
} else if (flag == 15) {
boundary[1][0] = read_int();
} else if (flag == 16) {
boundary[1][1] = read_int();
} else if (flag == 17) {
boundary[2][0] = read_int();
} else if (flag == 18) {
boundary[2][1] = read_int();

int flag = 0;
if ((xperiodic != domain->xperiodic || yperiodic != domain->yperiodic ||
zperiodic != domain->zperiodic)) flag = 1;
if (boundary[0][0] != domain->boundary[0][0] ||
boundary[0][1] != domain->boundary[0][1] ||
boundary[1][0] != domain->boundary[1][0] ||
boundary[1][1] != domain->boundary[1][1] ||
boundary[2][0] != domain->boundary[2][0] ||
boundary[2][1] != domain->boundary[2][1]) flag = 1;
if (flag && me == 0)
error->warning("Resetting boundary settings to restart file values");

domain->xperiodic = xperiodic;
domain->yperiodic = yperiodic;
domain->zperiodic = zperiodic;
domain->boundary[0][0] = boundary[0][0];
domain->boundary[0][1] = boundary[0][1];
domain->boundary[1][0] = boundary[1][0];
domain->boundary[1][1] = boundary[1][1];
domain->boundary[2][0] = boundary[2][0];
domain->boundary[2][1] = boundary[2][1];

domain->nonperiodic = 0;
if (xperiodic == 0 || yperiodic == 0 || zperiodic == 0) {
domain->nonperiodic = 1;
if (boundary[0][0] >= 2 || boundary[0][1] >= 2 ||
boundary[1][0] >= 2 || boundary[1][1] >= 2 ||
boundary[2][0] >= 2 || boundary[2][1] >= 2)
domain->nonperiodic = 2;
}

// create new AtomVec class
// if style = hybrid, read additional sub-class arguments

} else if (flag == 19) {
char *style = read_char();

int nwords = 0;
char **words = NULL;

if (strcmp(style,"hybrid") == 0) {
nwords = read_int();
words = new char*[nwords];
for (int i = 0; i < nwords; i++) words = read_char();
}

atom->create_avec(style,nwords,words);
for (int i = 0; i < nwords; i++) delete [] words;
delete [] words;
delete [] style;

} else if (flag == 20) {
atom->natoms = read_double();
} else if (flag == 21) {
atom->ntypes = read_int();
} else if (flag == 22) {
atom->nbonds = read_int();
} else if (flag == 23) {
atom->nbondtypes = read_int();
} else if (flag == 24) {
atom->bond_per_atom = read_int();
} else if (flag == 25) {
atom->nangles = read_int();
} else if (flag == 26) {
atom->nangletypes = read_int();
} else if (flag == 27) {
atom->angle_per_atom = read_int();
} else if (flag == 28) {
atom->ndihedrals = read_int();
} else if (flag == 29) {
atom->ndihedraltypes = read_int();
} else if (flag == 30) {
atom->dihedral_per_atom = read_int();
} else if (flag == 31) {
atom->nimpropers = read_int();
} else if (flag == 32) {
atom->nimpropertypes = read_int();
} else if (flag == 33) {
atom->improper_per_atom = read_int();

} else if (flag == 34) {
domain->boxxlo = read_double();
} else if (flag == 35) {
domain->boxxhi = read_double();
} else if (flag == 36) {
domain->boxylo = read_double();
} else if (flag == 37) {
domain->boxyhi = read_double();
} else if (flag == 38) {
domain->boxzlo = read_double();
} else if (flag == 39) {
domain->boxzhi = read_double();

} else if (flag == 40) {
force->special_lj[1] = read_double();
} else if (flag == 41) {
force->special_lj[2] = read_double();
} else if (flag == 42) {
force->special_lj[3] = read_double();
} else if (flag == 43) {
force->special_coul[1] = read_double();
} else if (flag == 44) {
force->special_coul[2] = read_double();
} else if (flag == 45) {
force->special_coul[3] = read_double();

} else error->all("Invalid flag in header of restart file");

flag = read_int();
}
}

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

void ReadRestart::mass()
{
double *mass = new double[atom->ntypes+1];
if (me == 0) fread(&mass[1],sizeof(double),atom->ntypes,fp);
MPI_Bcast(&mass[1],atom->ntypes,MPI_DOUBLE,0,world);
atom->set_mass(mass);
delete [] mass;
}

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

void ReadRestart::dipole()
{
double *dipole = new double[atom->ntypes+1];
if (me == 0) fread(&dipole[1],sizeof(double),atom->ntypes,fp);
MPI_Bcast(&dipole[1],atom->ntypes,MPI_DOUBLE,0,world);
atom->set_dipole(dipole);
delete [] dipole;
}

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

void ReadRestart::force_fields()
{
int n;
char *style;

if (me == 0) fread(&n,sizeof(int),1,fp);
MPI_Bcast(&n,1,MPI_INT,0,world);
if (n) {
style = (char *) memory->smalloc(n*sizeof(char),"read_restart:style");
if (me == 0) fread(style,sizeof(char),n,fp);
MPI_Bcast(style,n,MPI_CHAR,0,world);

if (force->pair == NULL || strcmp(style,force->pair_style)) {
if (force->pair) {
if (me == 0)
error->warning("Resetting pair_style to restart file value");
delete force->pair;
}
force->create_pair(style);
}

memory->sfree(style);
force->pair->read_restart(fp);
} else force->create_pair("none");

if (atom->avec->bonds_allow) {
if (me == 0) fread(&n,sizeof(int),1,fp);
MPI_Bcast(&n,1,MPI_INT,0,world);
if (n) {
style = (char *) memory->smalloc(n*sizeof(char),"read_restart:style");
if (me == 0) fread(style,sizeof(char),n,fp);
MPI_Bcast(style,n,MPI_CHAR,0,world);

if (force->bond == NULL || strcmp(style,force->bond_style)) {
if (force->bond) {
if (me == 0)
error->warning("Resetting bond_style to restart file value");
delete force->bond;
}
force->create_bond(style);
}

memory->sfree(style);
force->bond->read_restart(fp);
} else force->create_bond("none");
}

if (atom->avec->angles_allow) {
if (me == 0) fread(&n,sizeof(int),1,fp);
MPI_Bcast(&n,1,MPI_INT,0,world);
if (n) {
style = (char *) memory->smalloc(n*sizeof(char),"read_restart:style");
if (me == 0) fread(style,sizeof(char),n,fp);
MPI_Bcast(style,n,MPI_CHAR,0,world);

if (force->angle == NULL || strcmp(style,force->angle_style)) {
if (force->angle) {
if (me == 0)
error->warning("Resetting angle_style to restart file value");
delete force->angle;
}
force->create_angle(style);
}

memory->sfree(style);
force->angle->read_restart(fp);
} else force->create_angle("none");
}

if (atom->avec->dihedrals_allow) {
if (me == 0) fread(&n,sizeof(int),1,fp);
MPI_Bcast(&n,1,MPI_INT,0,world);
if (n) {
style = (char *) memory->smalloc(n*sizeof(char),"read_restart:style");
if (me == 0) fread(style,sizeof(char),n,fp);
MPI_Bcast(style,n,MPI_CHAR,0,world);

if (force->dihedral == NULL || strcmp(style,force->dihedral_style)) {
if (force->dihedral) {
if (me == 0)
error->warning("Resetting dihedral_style to restart file value");
delete force->dihedral;
}
force->create_dihedral(style);
}

memory->sfree(style);
force->dihedral->read_restart(fp);
} else force->create_dihedral("none");
}

if (atom->avec->impropers_allow) {
if (me == 0) fread(&n,sizeof(int),1,fp);
MPI_Bcast(&n,1,MPI_INT,0,world);
if (n) {
style = (char *) memory->smalloc(n*sizeof(char),"read_restart:style");
if (me == 0) fread(style,sizeof(char),n,fp);
MPI_Bcast(style,n,MPI_CHAR,0,world);

if (force->improper == NULL || strcmp(style,force->improper_style)) {
if (force->improper) {
if (me == 0)
error->warning("Resetting improper_style to restart file value");
delete force->improper;
}
force->create_improper(style);
}

memory->sfree(style);
force->improper->read_restart(fp);
} else force->create_improper("none");
}
}

/* ----------------------------------------------------------------------
read an int from restart file and bcast it
------------------------------------------------------------------------- */

int ReadRestart::read_int()
{
int value;
if (me == 0) fread(&value,sizeof(int),1,fp);
MPI_Bcast(&value,1,MPI_INT,0,world);
return value;
}

/* ----------------------------------------------------------------------
read a double from restart file and bcast it
------------------------------------------------------------------------- */

double ReadRestart::read_double()
{
double value;
if (me == 0) fread(&value,sizeof(double),1,fp);
MPI_Bcast(&value,1,MPI_DOUBLE,0,world);
return value;
}

/* ----------------------------------------------------------------------
read a char str from restart file and bcast it
str is allocated here, ptr is returned, caller must deallocate
------------------------------------------------------------------------- */

char *ReadRestart::read_char()
{
int n;
if (me == 0) fread(&n,sizeof(int),1,fp);
MPI_Bcast(&n,1,MPI_INT,0,world);
char *value = new char[n];
if (me == 0) fread(value,sizeof(char),n,fp);
MPI_Bcast(value,n,MPI_CHAR,0,world);
return value;
}