diff --git a/src/dump_custom.cpp b/src/dump_custom.cpp index 4e4b4224e..547dcfaa6 100644 --- a/src/dump_custom.cpp +++ b/src/dump_custom.cpp @@ -1,2409 +1,2410 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ #include "math.h" #include "stdlib.h" #include "string.h" #include "dump_custom.h" #include "atom.h" #include "force.h" #include "domain.h" #include "region.h" #include "group.h" #include "input.h" #include "variable.h" #include "update.h" #include "modify.h" #include "compute.h" #include "fix.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; // customize by adding keyword // also customize compute_atom_property.cpp enum{ID,MOL,TYPE,ELEMENT,MASS, X,Y,Z,XS,YS,ZS,XSTRI,YSTRI,ZSTRI,XU,YU,ZU,XUTRI,YUTRI,ZUTRI, XSU,YSU,ZSU,XSUTRI,YSUTRI,ZSUTRI, IX,IY,IZ, VX,VY,VZ,FX,FY,FZ, Q,MUX,MUY,MUZ,MU,RADIUS,DIAMETER, OMEGAX,OMEGAY,OMEGAZ,ANGMOMX,ANGMOMY,ANGMOMZ, TQX,TQY,TQZ,SPIN,ERADIUS,ERVEL,ERFORCE, COMPUTE,FIX,VARIABLE}; enum{LT,LE,GT,GE,EQ,NEQ}; enum{INT,DOUBLE,STRING}; #define INVOKED_PERATOM 8 /* ---------------------------------------------------------------------- */ DumpCustom::DumpCustom(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, arg) { if (narg == 5) error->all(FLERR,"No dump custom arguments specified"); clearstep = 1; nevery = force->inumeric(FLERR,arg[3]); // size_one may be shrunk below if additional optional args exist size_one = nfield = narg - 5; pack_choice = new FnPtrPack[nfield]; vtype = new int[nfield]; iregion = -1; idregion = NULL; nthresh = 0; thresh_array = NULL; thresh_op = NULL; thresh_value = NULL; // computes, fixes, variables which the dump accesses memory->create(field2index,nfield,"dump:field2index"); memory->create(argindex,nfield,"dump:argindex"); ncompute = 0; id_compute = NULL; compute = NULL; nfix = 0; id_fix = NULL; fix = NULL; nvariable = 0; id_variable = NULL; variable = NULL; vbuf = NULL; // process attributes // ioptional = start of additional optional args - // only dump image style processes optional args + // only dump image and dump movie styles processes optional args ioptional = parse_fields(narg,arg); - if (ioptional < narg && strcmp(style,"image") != 0) + if ((ioptional < narg) && + (strcmp(style,"image") != 0) && (strcmp(style,"movie") != 0)) error->all(FLERR,"Invalid attribute in dump custom command"); size_one = nfield = ioptional - 5; // atom selection arrays maxlocal = 0; choose = NULL; dchoose = NULL; clist = NULL; // element names ntypes = atom->ntypes; typenames = NULL; // setup format strings vformat = new char*[size_one]; format_default = new char[3*size_one+1]; format_default[0] = '\0'; for (int i = 0; i < size_one; i++) { if (vtype[i] == INT) strcat(format_default,"%d "); else if (vtype[i] == DOUBLE) strcat(format_default,"%g "); else if (vtype[i] == STRING) strcat(format_default,"%s "); vformat[i] = NULL; } // setup column string int n = 0; for (int iarg = 5; iarg < narg; iarg++) n += strlen(arg[iarg]) + 2; columns = new char[n]; columns[0] = '\0'; for (int iarg = 5; iarg < narg; iarg++) { strcat(columns,arg[iarg]); strcat(columns," "); } } /* ---------------------------------------------------------------------- */ DumpCustom::~DumpCustom() { delete [] pack_choice; delete [] vtype; memory->destroy(field2index); memory->destroy(argindex); delete [] idregion; memory->destroy(thresh_array); memory->destroy(thresh_op); memory->destroy(thresh_value); for (int i = 0; i < ncompute; i++) delete [] id_compute[i]; memory->sfree(id_compute); delete [] compute; for (int i = 0; i < nfix; i++) delete [] id_fix[i]; memory->sfree(id_fix); delete [] fix; for (int i = 0; i < nvariable; i++) delete [] id_variable[i]; memory->sfree(id_variable); delete [] variable; for (int i = 0; i < nvariable; i++) memory->destroy(vbuf[i]); delete [] vbuf; memory->destroy(choose); memory->destroy(dchoose); memory->destroy(clist); if (typenames) { for (int i = 1; i <= ntypes; i++) delete [] typenames[i]; delete [] typenames; } for (int i = 0; i < size_one; i++) delete [] vformat[i]; delete [] vformat; delete [] columns; } /* ---------------------------------------------------------------------- */ void DumpCustom::init_style() { delete [] format; char *str; if (format_user) str = format_user; else str = format_default; int n = strlen(str) + 1; format = new char[n]; strcpy(format,str); // default for element names = C if (typenames == NULL) { typenames = new char*[ntypes+1]; for (int itype = 1; itype <= ntypes; itype++) { typenames[itype] = new char[2]; strcpy(typenames[itype],"C"); } } // tokenize the format string and add space at end of each format element char *ptr; for (int i = 0; i < size_one; i++) { if (i == 0) ptr = strtok(format," \0"); else ptr = strtok(NULL," \0"); if (ptr == NULL) error->all(FLERR,"Dump_modify format string is too short"); delete [] vformat[i]; vformat[i] = new char[strlen(ptr) + 2]; strcpy(vformat[i],ptr); vformat[i] = strcat(vformat[i]," "); } // setup boundary string domain->boundary_string(boundstr); // setup function ptrs if (binary && domain->triclinic == 0) header_choice = &DumpCustom::header_binary; else if (binary && domain->triclinic == 1) header_choice = &DumpCustom::header_binary_triclinic; else if (!binary && domain->triclinic == 0) header_choice = &DumpCustom::header_item; else if (!binary && domain->triclinic == 1) header_choice = &DumpCustom::header_item_triclinic; if (binary) write_choice = &DumpCustom::write_binary; else write_choice = &DumpCustom::write_text; // find current ptr for each compute,fix,variable // check that fix frequency is acceptable int icompute; for (int i = 0; i < ncompute; i++) { icompute = modify->find_compute(id_compute[i]); if (icompute < 0) error->all(FLERR,"Could not find dump custom compute ID"); compute[i] = modify->compute[icompute]; } int ifix; for (int i = 0; i < nfix; i++) { ifix = modify->find_fix(id_fix[i]); if (ifix < 0) error->all(FLERR,"Could not find dump custom fix ID"); fix[i] = modify->fix[ifix]; if (nevery % modify->fix[ifix]->peratom_freq) error->all(FLERR,"Dump custom and fix not computed at compatible times"); } int ivariable; for (int i = 0; i < nvariable; i++) { ivariable = input->variable->find(id_variable[i]); if (ivariable < 0) error->all(FLERR,"Could not find dump custom variable name"); variable[i] = ivariable; } // set index and check validity of region if (iregion >= 0) { iregion = domain->find_region(idregion); if (iregion == -1) error->all(FLERR,"Region ID for dump custom does not exist"); } // open single file, one time only if (multifile == 0) openfile(); } /* ---------------------------------------------------------------------- */ void DumpCustom::write_header(bigint ndump) { if (multiproc) (this->*header_choice)(ndump); else if (me == 0) (this->*header_choice)(ndump); } /* ---------------------------------------------------------------------- */ void DumpCustom::header_binary(bigint ndump) { fwrite(&update->ntimestep,sizeof(bigint),1,fp); fwrite(&ndump,sizeof(bigint),1,fp); fwrite(&domain->triclinic,sizeof(int),1,fp); fwrite(&domain->boundary[0][0],6*sizeof(int),1,fp); fwrite(&boxxlo,sizeof(double),1,fp); fwrite(&boxxhi,sizeof(double),1,fp); fwrite(&boxylo,sizeof(double),1,fp); fwrite(&boxyhi,sizeof(double),1,fp); fwrite(&boxzlo,sizeof(double),1,fp); fwrite(&boxzhi,sizeof(double),1,fp); fwrite(&size_one,sizeof(int),1,fp); if (multiproc) fwrite(&nclusterprocs,sizeof(int),1,fp); else fwrite(&nprocs,sizeof(int),1,fp); } /* ---------------------------------------------------------------------- */ void DumpCustom::header_binary_triclinic(bigint ndump) { fwrite(&update->ntimestep,sizeof(bigint),1,fp); fwrite(&ndump,sizeof(bigint),1,fp); fwrite(&domain->triclinic,sizeof(int),1,fp); fwrite(&domain->boundary[0][0],6*sizeof(int),1,fp); fwrite(&boxxlo,sizeof(double),1,fp); fwrite(&boxxhi,sizeof(double),1,fp); fwrite(&boxylo,sizeof(double),1,fp); fwrite(&boxyhi,sizeof(double),1,fp); fwrite(&boxzlo,sizeof(double),1,fp); fwrite(&boxzhi,sizeof(double),1,fp); fwrite(&boxxy,sizeof(double),1,fp); fwrite(&boxxz,sizeof(double),1,fp); fwrite(&boxyz,sizeof(double),1,fp); fwrite(&size_one,sizeof(int),1,fp); if (multiproc) fwrite(&nclusterprocs,sizeof(int),1,fp); else fwrite(&nprocs,sizeof(int),1,fp); } /* ---------------------------------------------------------------------- */ void DumpCustom::header_item(bigint ndump) { fprintf(fp,"ITEM: TIMESTEP\n"); fprintf(fp,BIGINT_FORMAT "\n",update->ntimestep); fprintf(fp,"ITEM: NUMBER OF ATOMS\n"); fprintf(fp,BIGINT_FORMAT "\n",ndump); fprintf(fp,"ITEM: BOX BOUNDS %s\n",boundstr); fprintf(fp,"%g %g\n",boxxlo,boxxhi); fprintf(fp,"%g %g\n",boxylo,boxyhi); fprintf(fp,"%g %g\n",boxzlo,boxzhi); fprintf(fp,"ITEM: ATOMS %s\n",columns); } /* ---------------------------------------------------------------------- */ void DumpCustom::header_item_triclinic(bigint ndump) { fprintf(fp,"ITEM: TIMESTEP\n"); fprintf(fp,BIGINT_FORMAT "\n",update->ntimestep); fprintf(fp,"ITEM: NUMBER OF ATOMS\n"); fprintf(fp,BIGINT_FORMAT "\n",ndump); fprintf(fp,"ITEM: BOX BOUNDS xy xz yz %s\n",boundstr); fprintf(fp,"%g %g %g\n",boxxlo,boxxhi,boxxy); fprintf(fp,"%g %g %g\n",boxylo,boxyhi,boxxz); fprintf(fp,"%g %g %g\n",boxzlo,boxzhi,boxyz); fprintf(fp,"ITEM: ATOMS %s\n",columns); } /* ---------------------------------------------------------------------- */ int DumpCustom::count() { int i; // grow choose and variable vbuf arrays if needed int nlocal = atom->nlocal; if (nlocal > maxlocal) { maxlocal = atom->nmax; memory->destroy(choose); memory->destroy(dchoose); memory->destroy(clist); memory->create(choose,maxlocal,"dump:choose"); memory->create(dchoose,maxlocal,"dump:dchoose"); memory->create(clist,maxlocal,"dump:clist"); for (i = 0; i < nvariable; i++) { memory->destroy(vbuf[i]); memory->create(vbuf[i],maxlocal,"dump:vbuf"); } } // invoke Computes for per-atom quantities if (ncompute) { for (i = 0; i < ncompute; i++) if (!(compute[i]->invoked_flag & INVOKED_PERATOM)) { compute[i]->compute_peratom(); compute[i]->invoked_flag |= INVOKED_PERATOM; } } // evaluate atom-style Variables for per-atom quantities if (nvariable) for (i = 0; i < nvariable; i++) input->variable->compute_atom(variable[i],igroup,vbuf[i],1,0); // choose all local atoms for output for (i = 0; i < nlocal; i++) choose[i] = 1; // un-choose if not in group if (igroup) { int *mask = atom->mask; for (i = 0; i < nlocal; i++) if (!(mask[i] & groupbit)) choose[i] = 0; } // un-choose if not in region if (iregion >= 0) { Region *region = domain->regions[iregion]; double **x = atom->x; for (i = 0; i < nlocal; i++) if (choose[i] && region->match(x[i][0],x[i][1],x[i][2]) == 0) choose[i] = 0; } // un-choose if any threshhold criterion isn't met if (nthresh) { double *ptr; double value; int nstride; int nlocal = atom->nlocal; for (int ithresh = 0; ithresh < nthresh; ithresh++) { // customize by adding to if statement if (thresh_array[ithresh] == ID) { int *tag = atom->tag; for (i = 0; i < nlocal; i++) dchoose[i] = tag[i]; ptr = dchoose; nstride = 1; } else if (thresh_array[ithresh] == MOL) { if (!atom->molecule_flag) error->all(FLERR, "Threshhold for an atom property that isn't allocated"); int *molecule = atom->molecule; for (i = 0; i < nlocal; i++) dchoose[i] = molecule[i]; ptr = dchoose; nstride = 1; } else if (thresh_array[ithresh] == TYPE) { int *type = atom->type; for (i = 0; i < nlocal; i++) dchoose[i] = type[i]; ptr = dchoose; nstride = 1; } else if (thresh_array[ithresh] == ELEMENT) { int *type = atom->type; for (i = 0; i < nlocal; i++) dchoose[i] = type[i]; ptr = dchoose; nstride = 1; } else if (thresh_array[ithresh] == MASS) { if (atom->rmass) { ptr = atom->rmass; nstride = 1; } else { double *mass = atom->mass; int *type = atom->type; for (i = 0; i < nlocal; i++) dchoose[i] = mass[type[i]]; ptr = dchoose; nstride = 1; } } else if (thresh_array[ithresh] == X) { ptr = &atom->x[0][0]; nstride = 3; } else if (thresh_array[ithresh] == Y) { ptr = &atom->x[0][1]; nstride = 3; } else if (thresh_array[ithresh] == Z) { ptr = &atom->x[0][2]; nstride = 3; } else if (thresh_array[ithresh] == XS) { double **x = atom->x; double boxxlo = domain->boxlo[0]; double invxprd = 1.0/domain->xprd; for (i = 0; i < nlocal; i++) dchoose[i] = (x[i][0] - boxxlo) * invxprd; ptr = dchoose; nstride = 1; } else if (thresh_array[ithresh] == YS) { double **x = atom->x; double boxylo = domain->boxlo[1]; double invyprd = 1.0/domain->yprd; for (i = 0; i < nlocal; i++) dchoose[i] = (x[i][1] - boxylo) * invyprd; ptr = dchoose; nstride = 1; } else if (thresh_array[ithresh] == ZS) { double **x = atom->x; double boxzlo = domain->boxlo[2]; double invzprd = 1.0/domain->zprd; for (i = 0; i < nlocal; i++) dchoose[i] = (x[i][2] - boxzlo) * invzprd; ptr = dchoose; nstride = 1; } else if (thresh_array[ithresh] == XSTRI) { double **x = atom->x; double *boxlo = domain->boxlo; double *h_inv = domain->h_inv; for (i = 0; i < nlocal; i++) dchoose[i] = h_inv[0]*(x[i][0]-boxlo[0]) + h_inv[5]*(x[i][1]-boxlo[1]) + h_inv[4]*(x[i][2]-boxlo[2]); ptr = dchoose; nstride = 1; } else if (thresh_array[ithresh] == YSTRI) { double **x = atom->x; double *boxlo = domain->boxlo; double *h_inv = domain->h_inv; for (i = 0; i < nlocal; i++) dchoose[i] = h_inv[1]*(x[i][1]-boxlo[1]) + h_inv[3]*(x[i][2]-boxlo[2]); ptr = dchoose; nstride = 1; } else if (thresh_array[ithresh] == ZSTRI) { double **x = atom->x; double *boxlo = domain->boxlo; double *h_inv = domain->h_inv; for (i = 0; i < nlocal; i++) dchoose[i] = h_inv[2]*(x[i][2]-boxlo[2]); ptr = dchoose; nstride = 1; } else if (thresh_array[ithresh] == XU) { double **x = atom->x; tagint *image = atom->image; double xprd = domain->xprd; for (i = 0; i < nlocal; i++) dchoose[i] = x[i][0] + ((image[i] & IMGMASK) - IMGMAX) * xprd; ptr = dchoose; nstride = 1; } else if (thresh_array[ithresh] == YU) { double **x = atom->x; tagint *image = atom->image; double yprd = domain->yprd; for (i = 0; i < nlocal; i++) dchoose[i] = x[i][1] + ((image[i] >> IMGBITS & IMGMASK) - IMGMAX) * yprd; ptr = dchoose; nstride = 1; } else if (thresh_array[ithresh] == ZU) { double **x = atom->x; tagint *image = atom->image; double zprd = domain->zprd; for (i = 0; i < nlocal; i++) dchoose[i] = x[i][2] + ((image[i] >> IMG2BITS) - IMGMAX) * zprd; ptr = dchoose; nstride = 1; } else if (thresh_array[ithresh] == XUTRI) { double **x = atom->x; tagint *image = atom->image; double *h = domain->h; int xbox,ybox,zbox; for (i = 0; i < nlocal; i++) { xbox = (image[i] & IMGMASK) - IMGMAX; ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; zbox = (image[i] >> IMG2BITS) - IMGMAX; dchoose[i] = x[i][0] + h[0]*xbox + h[5]*ybox + h[4]*zbox; } ptr = dchoose; nstride = 1; } else if (thresh_array[ithresh] == YUTRI) { double **x = atom->x; tagint *image = atom->image; double *h = domain->h; int ybox,zbox; for (i = 0; i < nlocal; i++) { ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; zbox = (image[i] >> IMG2BITS) - IMGMAX; dchoose[i] = x[i][1] + h[1]*ybox + h[3]*zbox; } ptr = dchoose; nstride = 1; } else if (thresh_array[ithresh] == ZUTRI) { double **x = atom->x; tagint *image = atom->image; double *h = domain->h; int zbox; for (i = 0; i < nlocal; i++) { zbox = (image[i] >> IMG2BITS) - IMGMAX; dchoose[i] = x[i][2] + h[2]*zbox; } ptr = dchoose; nstride = 1; } else if (thresh_array[ithresh] == XSU) { double **x = atom->x; tagint *image = atom->image; double boxxlo = domain->boxlo[0]; double invxprd = 1.0/domain->xprd; for (i = 0; i < nlocal; i++) dchoose[i] = (x[i][0] - boxxlo) * invxprd + (image[i] & IMGMASK) - IMGMAX; ptr = dchoose; nstride = 1; } else if (thresh_array[ithresh] == YSU) { double **x = atom->x; tagint *image = atom->image; double boxylo = domain->boxlo[1]; double invyprd = 1.0/domain->yprd; for (i = 0; i < nlocal; i++) dchoose[i] = (x[i][1] - boxylo) * invyprd + (image[i] >> IMGBITS & IMGMASK) - IMGMAX; ptr = dchoose; nstride = 1; } else if (thresh_array[ithresh] == ZSU) { double **x = atom->x; tagint *image = atom->image; double boxzlo = domain->boxlo[2]; double invzprd = 1.0/domain->zprd; for (i = 0; i < nlocal; i++) dchoose[i] = (x[i][2] - boxzlo) * invzprd + (image[i] >> IMG2BITS) - IMGMAX; ptr = dchoose; nstride = 1; } else if (thresh_array[ithresh] == XSUTRI) { double **x = atom->x; tagint *image = atom->image; double *boxlo = domain->boxlo; double *h_inv = domain->h_inv; for (i = 0; i < nlocal; i++) dchoose[i] = h_inv[0]*(x[i][0]-boxlo[0]) + h_inv[5]*(x[i][1]-boxlo[1]) + h_inv[4]*(x[i][2]-boxlo[2]) + (image[i] & IMGMASK) - IMGMAX; ptr = dchoose; nstride = 1; } else if (thresh_array[ithresh] == YSUTRI) { double **x = atom->x; tagint *image = atom->image; double *boxlo = domain->boxlo; double *h_inv = domain->h_inv; for (i = 0; i < nlocal; i++) dchoose[i] = h_inv[1]*(x[i][1]-boxlo[1]) + h_inv[3]*(x[i][2]-boxlo[2]) + (image[i] >> IMGBITS & IMGMASK) - IMGMAX; ptr = dchoose; nstride = 1; } else if (thresh_array[ithresh] == ZSUTRI) { double **x = atom->x; tagint *image = atom->image; double *boxlo = domain->boxlo; double *h_inv = domain->h_inv; for (i = 0; i < nlocal; i++) dchoose[i] = h_inv[2]*(x[i][2]-boxlo[2]) + (image[i] >> IMG2BITS) - IMGMAX; ptr = dchoose; nstride = 1; } else if (thresh_array[ithresh] == IX) { tagint *image = atom->image; for (i = 0; i < nlocal; i++) dchoose[i] = (image[i] & IMGMASK) - IMGMAX; ptr = dchoose; nstride = 1; } else if (thresh_array[ithresh] == IY) { tagint *image = atom->image; for (i = 0; i < nlocal; i++) dchoose[i] = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; ptr = dchoose; nstride = 1; } else if (thresh_array[ithresh] == IZ) { tagint *image = atom->image; for (i = 0; i < nlocal; i++) dchoose[i] = (image[i] >> IMG2BITS) - IMGMAX; ptr = dchoose; nstride = 1; } else if (thresh_array[ithresh] == VX) { ptr = &atom->v[0][0]; nstride = 3; } else if (thresh_array[ithresh] == VY) { ptr = &atom->v[0][1]; nstride = 3; } else if (thresh_array[ithresh] == VZ) { ptr = &atom->v[0][2]; nstride = 3; } else if (thresh_array[ithresh] == FX) { ptr = &atom->f[0][0]; nstride = 3; } else if (thresh_array[ithresh] == FY) { ptr = &atom->f[0][1]; nstride = 3; } else if (thresh_array[ithresh] == FZ) { ptr = &atom->f[0][2]; nstride = 3; } else if (thresh_array[ithresh] == Q) { if (!atom->q_flag) error->all(FLERR, "Threshhold for an atom property that isn't allocated"); ptr = atom->q; nstride = 1; } else if (thresh_array[ithresh] == MUX) { if (!atom->mu_flag) error->all(FLERR, "Threshhold for an atom property that isn't allocated"); ptr = &atom->mu[0][0]; nstride = 4; } else if (thresh_array[ithresh] == MUY) { if (!atom->mu_flag) error->all(FLERR, "Threshhold for an atom property that isn't allocated"); ptr = &atom->mu[0][1]; nstride = 4; } else if (thresh_array[ithresh] == MUZ) { if (!atom->mu_flag) error->all(FLERR, "Threshhold for an atom property that isn't allocated"); ptr = &atom->mu[0][2]; nstride = 4; } else if (thresh_array[ithresh] == MU) { if (!atom->mu_flag) error->all(FLERR, "Threshhold for an atom property that isn't allocated"); ptr = &atom->mu[0][3]; nstride = 4; } else if (thresh_array[ithresh] == RADIUS) { if (!atom->radius_flag) error->all(FLERR, "Threshhold for an atom property that isn't allocated"); ptr = atom->radius; nstride = 1; } else if (thresh_array[ithresh] == DIAMETER) { if (!atom->radius_flag) error->all(FLERR, "Threshhold for an atom property that isn't allocated"); double *radius = atom->radius; for (i = 0; i < nlocal; i++) dchoose[i] = 2.0*radius[i]; ptr = dchoose; nstride = 1; } else if (thresh_array[ithresh] == OMEGAX) { if (!atom->omega_flag) error->all(FLERR, "Threshhold for an atom property that isn't allocated"); ptr = &atom->omega[0][0]; nstride = 3; } else if (thresh_array[ithresh] == OMEGAY) { if (!atom->omega_flag) error->all(FLERR, "Threshhold for an atom property that isn't allocated"); ptr = &atom->omega[0][1]; nstride = 3; } else if (thresh_array[ithresh] == OMEGAZ) { if (!atom->omega_flag) error->all(FLERR, "Threshhold for an atom property that isn't allocated"); ptr = &atom->omega[0][2]; nstride = 3; } else if (thresh_array[ithresh] == ANGMOMX) { if (!atom->angmom_flag) error->all(FLERR, "Threshhold for an atom property that isn't allocated"); ptr = &atom->angmom[0][0]; nstride = 3; } else if (thresh_array[ithresh] == ANGMOMY) { if (!atom->angmom_flag) error->all(FLERR, "Threshhold for an atom property that isn't allocated"); ptr = &atom->angmom[0][1]; nstride = 3; } else if (thresh_array[ithresh] == ANGMOMZ) { if (!atom->angmom_flag) error->all(FLERR, "Threshhold for an atom property that isn't allocated"); ptr = &atom->angmom[0][2]; nstride = 3; } else if (thresh_array[ithresh] == TQX) { if (!atom->torque_flag) error->all(FLERR, "Threshhold for an atom property that isn't allocated"); ptr = &atom->torque[0][0]; nstride = 3; } else if (thresh_array[ithresh] == TQY) { if (!atom->torque_flag) error->all(FLERR, "Threshhold for an atom property that isn't allocated"); ptr = &atom->torque[0][1]; nstride = 3; } else if (thresh_array[ithresh] == TQZ) { if (!atom->torque_flag) error->all(FLERR, "Threshhold for an atom property that isn't allocated"); ptr = &atom->torque[0][2]; nstride = 3; } else if (thresh_array[ithresh] == SPIN) { if (!atom->spin_flag) error->all(FLERR, "Threshhold for an atom property that isn't allocated"); int *spin = atom->spin; for (i = 0; i < nlocal; i++) dchoose[i] = spin[i]; ptr = dchoose; nstride = 1; } else if (thresh_array[ithresh] == ERADIUS) { if (!atom->eradius_flag) error->all(FLERR, "Threshhold for an atom property that isn't allocated"); ptr = atom->eradius; nstride = 1; } else if (thresh_array[ithresh] == ERVEL) { if (!atom->ervel_flag) error->all(FLERR, "Threshhold for an atom property that isn't allocated"); ptr = atom->ervel; nstride = 1; } else if (thresh_array[ithresh] == ERFORCE) { if (!atom->erforce_flag) error->all(FLERR, "Threshhold for an atom property that isn't allocated"); ptr = atom->erforce; nstride = 1; } else if (thresh_array[ithresh] == COMPUTE) { i = nfield + ithresh; if (argindex[i] == 0) { ptr = compute[field2index[i]]->vector_atom; nstride = 1; } else { ptr = &compute[field2index[i]]->array_atom[0][argindex[i]-1]; nstride = compute[field2index[i]]->size_peratom_cols; } } else if (thresh_array[ithresh] == FIX) { i = nfield + ithresh; if (argindex[i] == 0) { ptr = fix[field2index[i]]->vector_atom; nstride = 1; } else { ptr = &fix[field2index[i]]->array_atom[0][argindex[i]-1]; nstride = fix[field2index[i]]->size_peratom_cols; } } else if (thresh_array[ithresh] == VARIABLE) { i = nfield + ithresh; ptr = vbuf[field2index[i]]; nstride = 1; } // unselect atoms that don't meet threshhold criterion value = thresh_value[ithresh]; if (thresh_op[ithresh] == LT) { for (i = 0; i < nlocal; i++, ptr += nstride) if (choose[i] && *ptr >= value) choose[i] = 0; } else if (thresh_op[ithresh] == LE) { for (i = 0; i < nlocal; i++, ptr += nstride) if (choose[i] && *ptr > value) choose[i] = 0; } else if (thresh_op[ithresh] == GT) { for (i = 0; i < nlocal; i++, ptr += nstride) if (choose[i] && *ptr <= value) choose[i] = 0; } else if (thresh_op[ithresh] == GE) { for (i = 0; i < nlocal; i++, ptr += nstride) if (choose[i] && *ptr < value) choose[i] = 0; } else if (thresh_op[ithresh] == EQ) { for (i = 0; i < nlocal; i++, ptr += nstride) if (choose[i] && *ptr != value) choose[i] = 0; } else if (thresh_op[ithresh] == NEQ) { for (i = 0; i < nlocal; i++, ptr += nstride) if (choose[i] && *ptr == value) choose[i] = 0; } } } // compress choose flags into clist // nchoose = # of selected atoms // clist[i] = local index of each selected atom nchoose = 0; for (i = 0; i < nlocal; i++) if (choose[i]) clist[nchoose++] = i; return nchoose; } /* ---------------------------------------------------------------------- */ void DumpCustom::pack(int *ids) { for (int n = 0; n < size_one; n++) (this->*pack_choice[n])(n); if (ids) { int *tag = atom->tag; for (int i = 0; i < nchoose; i++) ids[i] = tag[clist[i]]; } } /* ---------------------------------------------------------------------- */ void DumpCustom::write_data(int n, double *mybuf) { (this->*write_choice)(n,mybuf); } /* ---------------------------------------------------------------------- */ void DumpCustom::write_binary(int n, double *mybuf) { n *= size_one; fwrite(&n,sizeof(int),1,fp); fwrite(mybuf,sizeof(double),n,fp); } /* ---------------------------------------------------------------------- */ void DumpCustom::write_text(int n, double *mybuf) { int i,j; int m = 0; for (i = 0; i < n; i++) { for (j = 0; j < size_one; j++) { if (vtype[j] == INT) fprintf(fp,vformat[j],static_cast<int> (mybuf[m])); else if (vtype[j] == DOUBLE) fprintf(fp,vformat[j],mybuf[m]); else if (vtype[j] == STRING) fprintf(fp,vformat[j],typenames[(int) mybuf[m]]); m++; } fprintf(fp,"\n"); } } /* ---------------------------------------------------------------------- */ int DumpCustom::parse_fields(int narg, char **arg) { // customize by adding to if statement int i; for (int iarg = 5; iarg < narg; iarg++) { i = iarg-5; if (strcmp(arg[iarg],"id") == 0) { pack_choice[i] = &DumpCustom::pack_id; vtype[i] = INT; } else if (strcmp(arg[iarg],"mol") == 0) { if (!atom->molecule_flag) error->all(FLERR,"Dumping an atom property that isn't allocated"); pack_choice[i] = &DumpCustom::pack_molecule; vtype[i] = INT; } else if (strcmp(arg[iarg],"type") == 0) { pack_choice[i] = &DumpCustom::pack_type; vtype[i] = INT; } else if (strcmp(arg[iarg],"element") == 0) { pack_choice[i] = &DumpCustom::pack_type; vtype[i] = STRING; } else if (strcmp(arg[iarg],"mass") == 0) { pack_choice[i] = &DumpCustom::pack_mass; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"x") == 0) { pack_choice[i] = &DumpCustom::pack_x; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"y") == 0) { pack_choice[i] = &DumpCustom::pack_y; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"z") == 0) { pack_choice[i] = &DumpCustom::pack_z; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"xs") == 0) { if (domain->triclinic) pack_choice[i] = &DumpCustom::pack_xs_triclinic; else pack_choice[i] = &DumpCustom::pack_xs; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"ys") == 0) { if (domain->triclinic) pack_choice[i] = &DumpCustom::pack_ys_triclinic; else pack_choice[i] = &DumpCustom::pack_ys; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"zs") == 0) { if (domain->triclinic) pack_choice[i] = &DumpCustom::pack_zs_triclinic; else pack_choice[i] = &DumpCustom::pack_zs; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"xu") == 0) { if (domain->triclinic) pack_choice[i] = &DumpCustom::pack_xu_triclinic; else pack_choice[i] = &DumpCustom::pack_xu; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"yu") == 0) { if (domain->triclinic) pack_choice[i] = &DumpCustom::pack_yu_triclinic; else pack_choice[i] = &DumpCustom::pack_yu; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"zu") == 0) { if (domain->triclinic) pack_choice[i] = &DumpCustom::pack_zu_triclinic; else pack_choice[i] = &DumpCustom::pack_zu; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"xsu") == 0) { if (domain->triclinic) pack_choice[i] = &DumpCustom::pack_xsu_triclinic; else pack_choice[i] = &DumpCustom::pack_xsu; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"ysu") == 0) { if (domain->triclinic) pack_choice[i] = &DumpCustom::pack_ysu_triclinic; else pack_choice[i] = &DumpCustom::pack_ysu; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"zsu") == 0) { if (domain->triclinic) pack_choice[i] = &DumpCustom::pack_zsu_triclinic; else pack_choice[i] = &DumpCustom::pack_zsu; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"ix") == 0) { pack_choice[i] = &DumpCustom::pack_ix; vtype[i] = INT; } else if (strcmp(arg[iarg],"iy") == 0) { pack_choice[i] = &DumpCustom::pack_iy; vtype[i] = INT; } else if (strcmp(arg[iarg],"iz") == 0) { pack_choice[i] = &DumpCustom::pack_iz; vtype[i] = INT; } else if (strcmp(arg[iarg],"vx") == 0) { pack_choice[i] = &DumpCustom::pack_vx; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"vy") == 0) { pack_choice[i] = &DumpCustom::pack_vy; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"vz") == 0) { pack_choice[i] = &DumpCustom::pack_vz; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"fx") == 0) { pack_choice[i] = &DumpCustom::pack_fx; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"fy") == 0) { pack_choice[i] = &DumpCustom::pack_fy; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"fz") == 0) { pack_choice[i] = &DumpCustom::pack_fz; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"q") == 0) { if (!atom->q_flag) error->all(FLERR,"Dumping an atom property that isn't allocated"); pack_choice[i] = &DumpCustom::pack_q; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"mux") == 0) { if (!atom->mu_flag) error->all(FLERR,"Dumping an atom property that isn't allocated"); pack_choice[i] = &DumpCustom::pack_mux; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"muy") == 0) { if (!atom->mu_flag) error->all(FLERR,"Dumping an atom property that isn't allocated"); pack_choice[i] = &DumpCustom::pack_muy; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"muz") == 0) { if (!atom->mu_flag) error->all(FLERR,"Dumping an atom property that isn't allocated"); pack_choice[i] = &DumpCustom::pack_muz; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"mu") == 0) { if (!atom->mu_flag) error->all(FLERR,"Dumping an atom property that isn't allocated"); pack_choice[i] = &DumpCustom::pack_mu; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"radius") == 0) { if (!atom->radius_flag) error->all(FLERR,"Dumping an atom property that isn't allocated"); pack_choice[i] = &DumpCustom::pack_radius; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"diameter") == 0) { if (!atom->radius_flag) error->all(FLERR,"Dumping an atom property that isn't allocated"); pack_choice[i] = &DumpCustom::pack_diameter; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"omegax") == 0) { if (!atom->omega_flag) error->all(FLERR,"Dumping an atom property that isn't allocated"); pack_choice[i] = &DumpCustom::pack_omegax; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"omegay") == 0) { if (!atom->omega_flag) error->all(FLERR,"Dumping an atom property that isn't allocated"); pack_choice[i] = &DumpCustom::pack_omegay; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"omegaz") == 0) { if (!atom->omega_flag) error->all(FLERR,"Dumping an atom property that isn't allocated"); pack_choice[i] = &DumpCustom::pack_omegaz; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"angmomx") == 0) { if (!atom->angmom_flag) error->all(FLERR,"Dumping an atom property that isn't allocated"); pack_choice[i] = &DumpCustom::pack_angmomx; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"angmomy") == 0) { if (!atom->angmom_flag) error->all(FLERR,"Dumping an atom property that isn't allocated"); pack_choice[i] = &DumpCustom::pack_angmomy; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"angmomz") == 0) { if (!atom->angmom_flag) error->all(FLERR,"Dumping an atom property that isn't allocated"); pack_choice[i] = &DumpCustom::pack_angmomz; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"tqx") == 0) { if (!atom->torque_flag) error->all(FLERR,"Dumping an atom property that isn't allocated"); pack_choice[i] = &DumpCustom::pack_tqx; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"tqy") == 0) { if (!atom->torque_flag) error->all(FLERR,"Dumping an atom property that isn't allocated"); pack_choice[i] = &DumpCustom::pack_tqy; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"tqz") == 0) { if (!atom->torque_flag) error->all(FLERR,"Dumping an atom property that isn't allocated"); pack_choice[i] = &DumpCustom::pack_tqz; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"spin") == 0) { if (!atom->spin_flag) error->all(FLERR,"Dumping an atom quantity that isn't allocated"); pack_choice[i] = &DumpCustom::pack_spin; vtype[i] = INT; } else if (strcmp(arg[iarg],"eradius") == 0) { if (!atom->eradius_flag) error->all(FLERR,"Dumping an atom quantity that isn't allocated"); pack_choice[i] = &DumpCustom::pack_eradius; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"ervel") == 0) { if (!atom->ervel_flag) error->all(FLERR,"Dumping an atom quantity that isn't allocated"); pack_choice[i] = &DumpCustom::pack_ervel; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"erforce") == 0) { if (!atom->erforce_flag) error->all(FLERR,"Dumping an atom quantity that isn't allocated"); pack_choice[i] = &DumpCustom::pack_erforce; vtype[i] = DOUBLE; // compute value = c_ID // if no trailing [], then arg is set to 0, else arg is int between [] } else if (strncmp(arg[iarg],"c_",2) == 0) { pack_choice[i] = &DumpCustom::pack_compute; vtype[i] = DOUBLE; int n = strlen(arg[iarg]); char *suffix = new char[n]; strcpy(suffix,&arg[iarg][2]); char *ptr = strchr(suffix,'['); if (ptr) { if (suffix[strlen(suffix)-1] != ']') error->all(FLERR,"Invalid attribute in dump custom command"); argindex[i] = atoi(ptr+1); *ptr = '\0'; } else argindex[i] = 0; n = modify->find_compute(suffix); if (n < 0) error->all(FLERR,"Could not find dump custom compute ID"); if (modify->compute[n]->peratom_flag == 0) error->all(FLERR,"Dump custom compute does not compute per-atom info"); if (argindex[i] == 0 && modify->compute[n]->size_peratom_cols > 0) error->all(FLERR,"Dump custom compute does not calculate per-atom vector"); if (argindex[i] > 0 && modify->compute[n]->size_peratom_cols == 0) error->all(FLERR,"Dump custom compute does not calculate per-atom array"); if (argindex[i] > 0 && argindex[i] > modify->compute[n]->size_peratom_cols) error->all(FLERR,"Dump custom compute vector is accessed out-of-range"); field2index[i] = add_compute(suffix); delete [] suffix; // fix value = f_ID // if no trailing [], then arg is set to 0, else arg is between [] } else if (strncmp(arg[iarg],"f_",2) == 0) { pack_choice[i] = &DumpCustom::pack_fix; vtype[i] = DOUBLE; int n = strlen(arg[iarg]); char *suffix = new char[n]; strcpy(suffix,&arg[iarg][2]); char *ptr = strchr(suffix,'['); if (ptr) { if (suffix[strlen(suffix)-1] != ']') error->all(FLERR,"Invalid attribute in dump custom command"); argindex[i] = atoi(ptr+1); *ptr = '\0'; } else argindex[i] = 0; n = modify->find_fix(suffix); if (n < 0) error->all(FLERR,"Could not find dump custom fix ID"); if (modify->fix[n]->peratom_flag == 0) error->all(FLERR,"Dump custom fix does not compute per-atom info"); if (argindex[i] == 0 && modify->fix[n]->size_peratom_cols > 0) error->all(FLERR,"Dump custom fix does not compute per-atom vector"); if (argindex[i] > 0 && modify->fix[n]->size_peratom_cols == 0) error->all(FLERR,"Dump custom fix does not compute per-atom array"); if (argindex[i] > 0 && argindex[i] > modify->fix[n]->size_peratom_cols) error->all(FLERR,"Dump custom fix vector is accessed out-of-range"); field2index[i] = add_fix(suffix); delete [] suffix; // variable value = v_name } else if (strncmp(arg[iarg],"v_",2) == 0) { pack_choice[i] = &DumpCustom::pack_variable; vtype[i] = DOUBLE; int n = strlen(arg[iarg]); char *suffix = new char[n]; strcpy(suffix,&arg[iarg][2]); argindex[i] = 0; n = input->variable->find(suffix); if (n < 0) error->all(FLERR,"Could not find dump custom variable name"); if (input->variable->atomstyle(n) == 0) error->all(FLERR,"Dump custom variable is not atom-style variable"); field2index[i] = add_variable(suffix); delete [] suffix; } else return iarg; } return narg; } /* ---------------------------------------------------------------------- add Compute to list of Compute objects used by dump return index of where this Compute is in list if already in list, do not add, just return index, else add to list ------------------------------------------------------------------------- */ int DumpCustom::add_compute(char *id) { int icompute; for (icompute = 0; icompute < ncompute; icompute++) if (strcmp(id,id_compute[icompute]) == 0) break; if (icompute < ncompute) return icompute; id_compute = (char **) memory->srealloc(id_compute,(ncompute+1)*sizeof(char *),"dump:id_compute"); delete [] compute; compute = new Compute*[ncompute+1]; int n = strlen(id) + 1; id_compute[ncompute] = new char[n]; strcpy(id_compute[ncompute],id); ncompute++; return ncompute-1; } /* ---------------------------------------------------------------------- add Fix to list of Fix objects used by dump return index of where this Fix is in list if already in list, do not add, just return index, else add to list ------------------------------------------------------------------------- */ int DumpCustom::add_fix(char *id) { int ifix; for (ifix = 0; ifix < nfix; ifix++) if (strcmp(id,id_fix[ifix]) == 0) break; if (ifix < nfix) return ifix; id_fix = (char **) memory->srealloc(id_fix,(nfix+1)*sizeof(char *),"dump:id_fix"); delete [] fix; fix = new Fix*[nfix+1]; int n = strlen(id) + 1; id_fix[nfix] = new char[n]; strcpy(id_fix[nfix],id); nfix++; return nfix-1; } /* ---------------------------------------------------------------------- add Variable to list of Variables used by dump return index of where this Variable is in list if already in list, do not add, just return index, else add to list ------------------------------------------------------------------------- */ int DumpCustom::add_variable(char *id) { int ivariable; for (ivariable = 0; ivariable < nvariable; ivariable++) if (strcmp(id,id_variable[ivariable]) == 0) break; if (ivariable < nvariable) return ivariable; id_variable = (char **) memory->srealloc(id_variable,(nvariable+1)*sizeof(char *), "dump:id_variable"); delete [] variable; variable = new int[nvariable+1]; delete [] vbuf; vbuf = new double*[nvariable+1]; for (int i = 0; i <= nvariable; i++) vbuf[i] = NULL; int n = strlen(id) + 1; id_variable[nvariable] = new char[n]; strcpy(id_variable[nvariable],id); nvariable++; return nvariable-1; } /* ---------------------------------------------------------------------- */ int DumpCustom::modify_param(int narg, char **arg) { if (strcmp(arg[0],"region") == 0) { if (narg < 2) error->all(FLERR,"Illegal dump_modify command"); if (strcmp(arg[1],"none") == 0) iregion = -1; else { iregion = domain->find_region(arg[1]); if (iregion == -1) error->all(FLERR,"Dump_modify region ID does not exist"); int n = strlen(arg[1]) + 1; idregion = new char[n]; strcpy(idregion,arg[1]); } return 2; } if (strcmp(arg[0],"element") == 0) { if (narg < ntypes+1) error->all(FLERR,"Dump modify element names do not match atom types"); if (typenames) { for (int i = 1; i <= ntypes; i++) delete [] typenames[i]; delete [] typenames; typenames = NULL; } typenames = new char*[ntypes+1]; for (int itype = 1; itype <= ntypes; itype++) { int n = strlen(arg[itype]) + 1; typenames[itype] = new char[n]; strcpy(typenames[itype],arg[itype]); } return ntypes+1; } if (strcmp(arg[0],"thresh") == 0) { if (narg < 2) error->all(FLERR,"Illegal dump_modify command"); if (strcmp(arg[1],"none") == 0) { if (nthresh) { memory->destroy(thresh_array); memory->destroy(thresh_op); memory->destroy(thresh_value); thresh_array = NULL; thresh_op = NULL; thresh_value = NULL; } nthresh = 0; return 2; } if (narg < 4) error->all(FLERR,"Illegal dump_modify command"); // grow threshhold arrays memory->grow(thresh_array,nthresh+1,"dump:thresh_array"); memory->grow(thresh_op,(nthresh+1),"dump:thresh_op"); memory->grow(thresh_value,(nthresh+1),"dump:thresh_value"); // set attribute type of threshhold // customize by adding to if statement if (strcmp(arg[1],"id") == 0) thresh_array[nthresh] = ID; else if (strcmp(arg[1],"mol") == 0) thresh_array[nthresh] = MOL; else if (strcmp(arg[1],"type") == 0) thresh_array[nthresh] = TYPE; else if (strcmp(arg[1],"mass") == 0) thresh_array[nthresh] = MASS; else if (strcmp(arg[1],"x") == 0) thresh_array[nthresh] = X; else if (strcmp(arg[1],"y") == 0) thresh_array[nthresh] = Y; else if (strcmp(arg[1],"z") == 0) thresh_array[nthresh] = Z; else if (strcmp(arg[1],"xs") == 0 && domain->triclinic == 0) thresh_array[nthresh] = XS; else if (strcmp(arg[1],"xs") == 0 && domain->triclinic == 1) thresh_array[nthresh] = XSTRI; else if (strcmp(arg[1],"ys") == 0 && domain->triclinic == 0) thresh_array[nthresh] = YS; else if (strcmp(arg[1],"ys") == 0 && domain->triclinic == 1) thresh_array[nthresh] = YSTRI; else if (strcmp(arg[1],"zs") == 0 && domain->triclinic == 0) thresh_array[nthresh] = ZS; else if (strcmp(arg[1],"zs") == 0 && domain->triclinic == 1) thresh_array[nthresh] = ZSTRI; else if (strcmp(arg[1],"xu") == 0 && domain->triclinic == 0) thresh_array[nthresh] = XU; else if (strcmp(arg[1],"xu") == 0 && domain->triclinic == 1) thresh_array[nthresh] = XUTRI; else if (strcmp(arg[1],"yu") == 0 && domain->triclinic == 0) thresh_array[nthresh] = YU; else if (strcmp(arg[1],"yu") == 0 && domain->triclinic == 1) thresh_array[nthresh] = YUTRI; else if (strcmp(arg[1],"zu") == 0 && domain->triclinic == 0) thresh_array[nthresh] = ZU; else if (strcmp(arg[1],"zu") == 0 && domain->triclinic == 1) thresh_array[nthresh] = ZUTRI; else if (strcmp(arg[1],"xsu") == 0 && domain->triclinic == 0) thresh_array[nthresh] = XSU; else if (strcmp(arg[1],"xsu") == 0 && domain->triclinic == 1) thresh_array[nthresh] = XSUTRI; else if (strcmp(arg[1],"ysu") == 0 && domain->triclinic == 0) thresh_array[nthresh] = YSU; else if (strcmp(arg[1],"ysu") == 0 && domain->triclinic == 1) thresh_array[nthresh] = YSUTRI; else if (strcmp(arg[1],"zsu") == 0 && domain->triclinic == 0) thresh_array[nthresh] = ZSU; else if (strcmp(arg[1],"zsu") == 0 && domain->triclinic == 1) thresh_array[nthresh] = ZSUTRI; else if (strcmp(arg[1],"ix") == 0) thresh_array[nthresh] = IX; else if (strcmp(arg[1],"iy") == 0) thresh_array[nthresh] = IY; else if (strcmp(arg[1],"iz") == 0) thresh_array[nthresh] = IZ; else if (strcmp(arg[1],"vx") == 0) thresh_array[nthresh] = VX; else if (strcmp(arg[1],"vy") == 0) thresh_array[nthresh] = VY; else if (strcmp(arg[1],"vz") == 0) thresh_array[nthresh] = VZ; else if (strcmp(arg[1],"fx") == 0) thresh_array[nthresh] = FX; else if (strcmp(arg[1],"fy") == 0) thresh_array[nthresh] = FY; else if (strcmp(arg[1],"fz") == 0) thresh_array[nthresh] = FZ; else if (strcmp(arg[1],"q") == 0) thresh_array[nthresh] = Q; else if (strcmp(arg[1],"mux") == 0) thresh_array[nthresh] = MUX; else if (strcmp(arg[1],"muy") == 0) thresh_array[nthresh] = MUY; else if (strcmp(arg[1],"muz") == 0) thresh_array[nthresh] = MUZ; else if (strcmp(arg[1],"mu") == 0) thresh_array[nthresh] = MU; else if (strcmp(arg[1],"radius") == 0) thresh_array[nthresh] = RADIUS; else if (strcmp(arg[1],"diameter") == 0) thresh_array[nthresh] = DIAMETER; else if (strcmp(arg[1],"omegax") == 0) thresh_array[nthresh] = OMEGAX; else if (strcmp(arg[1],"omegay") == 0) thresh_array[nthresh] = OMEGAY; else if (strcmp(arg[1],"omegaz") == 0) thresh_array[nthresh] = OMEGAZ; else if (strcmp(arg[1],"angmomx") == 0) thresh_array[nthresh] = ANGMOMX; else if (strcmp(arg[1],"angmomy") == 0) thresh_array[nthresh] = ANGMOMY; else if (strcmp(arg[1],"angmomz") == 0) thresh_array[nthresh] = ANGMOMZ; else if (strcmp(arg[1],"tqx") == 0) thresh_array[nthresh] = TQX; else if (strcmp(arg[1],"tqy") == 0) thresh_array[nthresh] = TQY; else if (strcmp(arg[1],"tqz") == 0) thresh_array[nthresh] = TQZ; else if (strcmp(arg[1],"spin") == 0) thresh_array[nthresh] = SPIN; else if (strcmp(arg[1],"eradius") == 0) thresh_array[nthresh] = ERADIUS; else if (strcmp(arg[1],"ervel") == 0) thresh_array[nthresh] = ERVEL; else if (strcmp(arg[1],"erforce") == 0) thresh_array[nthresh] = ERFORCE; // compute value = c_ID // if no trailing [], then arg is set to 0, else arg is between [] // must grow field2index and argindex arrays, since access is beyond nfield else if (strncmp(arg[1],"c_",2) == 0) { thresh_array[nthresh] = COMPUTE; memory->grow(field2index,nfield+nthresh+1,"dump:field2index"); memory->grow(argindex,nfield+nthresh+1,"dump:argindex"); int n = strlen(arg[1]); char *suffix = new char[n]; strcpy(suffix,&arg[1][2]); char *ptr = strchr(suffix,'['); if (ptr) { if (suffix[strlen(suffix)-1] != ']') error->all(FLERR,"Invalid attribute in dump modify command"); argindex[nfield+nthresh] = atoi(ptr+1); *ptr = '\0'; } else argindex[nfield+nthresh] = 0; n = modify->find_compute(suffix); if (n < 0) error->all(FLERR,"Could not find dump modify compute ID"); if (modify->compute[n]->peratom_flag == 0) error->all(FLERR, "Dump modify compute ID does not compute per-atom info"); if (argindex[nfield+nthresh] == 0 && modify->compute[n]->size_peratom_cols > 0) error->all(FLERR, "Dump modify compute ID does not compute per-atom vector"); if (argindex[nfield+nthresh] > 0 && modify->compute[n]->size_peratom_cols == 0) error->all(FLERR, "Dump modify compute ID does not compute per-atom array"); if (argindex[nfield+nthresh] > 0 && argindex[nfield+nthresh] > modify->compute[n]->size_peratom_cols) error->all(FLERR,"Dump modify compute ID vector is not large enough"); field2index[nfield+nthresh] = add_compute(suffix); delete [] suffix; // fix value = f_ID // if no trailing [], then arg is set to 0, else arg is between [] // must grow field2index and argindex arrays, since access is beyond nfield } else if (strncmp(arg[1],"f_",2) == 0) { thresh_array[nthresh] = FIX; memory->grow(field2index,nfield+nthresh+1,"dump:field2index"); memory->grow(argindex,nfield+nthresh+1,"dump:argindex"); int n = strlen(arg[1]); char *suffix = new char[n]; strcpy(suffix,&arg[1][2]); char *ptr = strchr(suffix,'['); if (ptr) { if (suffix[strlen(suffix)-1] != ']') error->all(FLERR,"Invalid attribute in dump modify command"); argindex[nfield+nthresh] = atoi(ptr+1); *ptr = '\0'; } else argindex[nfield+nthresh] = 0; n = modify->find_fix(suffix); if (n < 0) error->all(FLERR,"Could not find dump modify fix ID"); if (modify->fix[n]->peratom_flag == 0) error->all(FLERR,"Dump modify fix ID does not compute per-atom info"); if (argindex[nfield+nthresh] == 0 && modify->fix[n]->size_peratom_cols > 0) error->all(FLERR,"Dump modify fix ID does not compute per-atom vector"); if (argindex[nfield+nthresh] > 0 && modify->fix[n]->size_peratom_cols == 0) error->all(FLERR,"Dump modify fix ID does not compute per-atom array"); if (argindex[nfield+nthresh] > 0 && argindex[nfield+nthresh] > modify->fix[n]->size_peratom_cols) error->all(FLERR,"Dump modify fix ID vector is not large enough"); field2index[nfield+nthresh] = add_fix(suffix); delete [] suffix; // variable value = v_ID // must grow field2index and argindex arrays, since access is beyond nfield } else if (strncmp(arg[1],"v_",2) == 0) { thresh_array[nthresh] = VARIABLE; memory->grow(field2index,nfield+nthresh+1,"dump:field2index"); memory->grow(argindex,nfield+nthresh+1,"dump:argindex"); int n = strlen(arg[1]); char *suffix = new char[n]; strcpy(suffix,&arg[1][2]); argindex[nfield+nthresh] = 0; n = input->variable->find(suffix); if (n < 0) error->all(FLERR,"Could not find dump modify variable name"); if (input->variable->atomstyle(n) == 0) error->all(FLERR,"Dump modify variable is not atom-style variable"); field2index[nfield+nthresh] = add_variable(suffix); delete [] suffix; } else error->all(FLERR,"Invalid dump_modify threshhold operator"); // set operation type of threshhold if (strcmp(arg[2],"<") == 0) thresh_op[nthresh] = LT; else if (strcmp(arg[2],"<=") == 0) thresh_op[nthresh] = LE; else if (strcmp(arg[2],">") == 0) thresh_op[nthresh] = GT; else if (strcmp(arg[2],">=") == 0) thresh_op[nthresh] = GE; else if (strcmp(arg[2],"==") == 0) thresh_op[nthresh] = EQ; else if (strcmp(arg[2],"!=") == 0) thresh_op[nthresh] = NEQ; else error->all(FLERR,"Invalid dump_modify threshhold operator"); // set threshhold value thresh_value[nthresh] = force->numeric(FLERR,arg[3]); nthresh++; return 4; } return 0; } /* ---------------------------------------------------------------------- return # of bytes of allocated memory in buf, choose, variable arrays ------------------------------------------------------------------------- */ bigint DumpCustom::memory_usage() { bigint bytes = Dump::memory_usage(); bytes += memory->usage(choose,maxlocal); bytes += memory->usage(dchoose,maxlocal); bytes += memory->usage(clist,maxlocal); bytes += memory->usage(vbuf,nvariable,maxlocal); return bytes; } /* ---------------------------------------------------------------------- extraction of Compute, Fix, Variable results ------------------------------------------------------------------------- */ void DumpCustom::pack_compute(int n) { double *vector = compute[field2index[n]]->vector_atom; double **array = compute[field2index[n]]->array_atom; int index = argindex[n]; if (index == 0) { for (int i = 0; i < nchoose; i++) { buf[n] = vector[clist[i]]; n += size_one; } } else { index--; for (int i = 0; i < nchoose; i++) { buf[n] = array[clist[i]][index]; n += size_one; } } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_fix(int n) { double *vector = fix[field2index[n]]->vector_atom; double **array = fix[field2index[n]]->array_atom; int index = argindex[n]; if (index == 0) { for (int i = 0; i < nchoose; i++) { buf[n] = vector[clist[i]]; n += size_one; } } else { index--; for (int i = 0; i < nchoose; i++) { buf[n] = array[clist[i]][index]; n += size_one; } } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_variable(int n) { double *vector = vbuf[field2index[n]]; for (int i = 0; i < nchoose; i++) { buf[n] = vector[clist[i]]; n += size_one; } } /* ---------------------------------------------------------------------- one method for every attribute dump custom can output the atom property is packed into buf starting at n with stride size_one customize a new attribute by adding a method ------------------------------------------------------------------------- */ void DumpCustom::pack_id(int n) { int *tag = atom->tag; for (int i = 0; i < nchoose; i++) { buf[n] = tag[clist[i]]; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_molecule(int n) { int *molecule = atom->molecule; for (int i = 0; i < nchoose; i++) { buf[n] = molecule[clist[i]]; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_type(int n) { int *type = atom->type; for (int i = 0; i < nchoose; i++) { buf[n] = type[clist[i]]; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_mass(int n) { int *type = atom->type; double *mass = atom->mass; double *rmass = atom->rmass; if (rmass) { for (int i = 0; i < nchoose; i++) { buf[n] = rmass[clist[i]]; n += size_one; } } else { for (int i = 0; i < nchoose; i++) { buf[n] = mass[type[clist[i]]]; n += size_one; } } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_x(int n) { double **x = atom->x; for (int i = 0; i < nchoose; i++) { buf[n] = x[clist[i]][0]; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_y(int n) { double **x = atom->x; for (int i = 0; i < nchoose; i++) { buf[n] = x[clist[i]][1]; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_z(int n) { double **x = atom->x; for (int i = 0; i < nchoose; i++) { buf[n] = x[clist[i]][2]; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_xs(int n) { double **x = atom->x; double boxxlo = domain->boxlo[0]; double invxprd = 1.0/domain->xprd; for (int i = 0; i < nchoose; i++) { buf[n] = (x[clist[i]][0] - boxxlo) * invxprd; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_ys(int n) { double **x = atom->x; double boxylo = domain->boxlo[1]; double invyprd = 1.0/domain->yprd; for (int i = 0; i < nchoose; i++) { buf[n] = (x[clist[i]][1] - boxylo) * invyprd; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_zs(int n) { double **x = atom->x; double boxzlo = domain->boxlo[2]; double invzprd = 1.0/domain->zprd; for (int i = 0; i < nchoose; i++) { buf[n] = (x[clist[i]][2] - boxzlo) * invzprd; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_xs_triclinic(int n) { int j; double **x = atom->x; double *boxlo = domain->boxlo; double *h_inv = domain->h_inv; for (int i = 0; i < nchoose; i++) { j = clist[i]; buf[n] = h_inv[0]*(x[j][0]-boxlo[0]) + h_inv[5]*(x[j][1]-boxlo[1]) + h_inv[4]*(x[j][2]-boxlo[2]); n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_ys_triclinic(int n) { int j; double **x = atom->x; double *boxlo = domain->boxlo; double *h_inv = domain->h_inv; for (int i = 0; i < nchoose; i++) { j = clist[i]; buf[n] = h_inv[1]*(x[j][1]-boxlo[1]) + h_inv[3]*(x[j][2]-boxlo[2]); n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_zs_triclinic(int n) { double **x = atom->x; double *boxlo = domain->boxlo; double *h_inv = domain->h_inv; for (int i = 0; i < nchoose; i++) { buf[n] = h_inv[2]*(x[clist[i]][2]-boxlo[2]); n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_xu(int n) { int j; double **x = atom->x; tagint *image = atom->image; double xprd = domain->xprd; for (int i = 0; i < nchoose; i++) { j = clist[i]; buf[n] = x[j][0] + ((image[j] & IMGMASK) - IMGMAX) * xprd; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_yu(int n) { int j; double **x = atom->x; tagint *image = atom->image; double yprd = domain->yprd; for (int i = 0; i < nchoose; i++) { j = clist[i]; buf[n] = x[j][1] + ((image[j] >> IMGBITS & IMGMASK) - IMGMAX) * yprd; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_zu(int n) { int j; double **x = atom->x; tagint *image = atom->image; double zprd = domain->zprd; for (int i = 0; i < nchoose; i++) { j = clist[i]; buf[n] = x[j][2] + ((image[j] >> IMG2BITS) - IMGMAX) * zprd; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_xu_triclinic(int n) { int j; double **x = atom->x; tagint *image = atom->image; double *h = domain->h; int xbox,ybox,zbox; for (int i = 0; i < nchoose; i++) { j = clist[i]; xbox = (image[j] & IMGMASK) - IMGMAX; ybox = (image[j] >> IMGBITS & IMGMASK) - IMGMAX; zbox = (image[j] >> IMG2BITS) - IMGMAX; buf[n] = x[j][0] + h[0]*xbox + h[5]*ybox + h[4]*zbox; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_yu_triclinic(int n) { int j; double **x = atom->x; tagint *image = atom->image; double *h = domain->h; int ybox,zbox; for (int i = 0; i < nchoose; i++) { j = clist[i]; ybox = (image[j] >> IMGBITS & IMGMASK) - IMGMAX; zbox = (image[j] >> IMG2BITS) - IMGMAX; buf[n] = x[j][1] + h[1]*ybox + h[3]*zbox; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_zu_triclinic(int n) { int j; double **x = atom->x; tagint *image = atom->image; double *h = domain->h; int zbox; for (int i = 0; i < nchoose; i++) { j = clist[i]; zbox = (image[j] >> IMG2BITS) - IMGMAX; buf[n] = x[j][2] + h[2]*zbox; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_xsu(int n) { int j; double **x = atom->x; tagint *image = atom->image; double boxxlo = domain->boxlo[0]; double invxprd = 1.0/domain->xprd; for (int i = 0; i < nchoose; i++) { j = clist[i]; buf[n] = (x[j][0] - boxxlo) * invxprd + (image[j] & IMGMASK) - IMGMAX; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_ysu(int n) { int j; double **x = atom->x; tagint *image = atom->image; double boxylo = domain->boxlo[1]; double invyprd = 1.0/domain->yprd; for (int i = 0; i < nchoose; i++) { j = clist[i]; buf[n] = (x[j][1] - boxylo) * invyprd + (image[j] >> IMGBITS & IMGMASK) - IMGMAX; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_zsu(int n) { int j; double **x = atom->x; tagint *image = atom->image; double boxzlo = domain->boxlo[2]; double invzprd = 1.0/domain->zprd; for (int i = 0; i < nchoose; i++) { j = clist[i]; buf[n] = (x[j][2] - boxzlo) * invzprd + (image[j] >> IMG2BITS) - IMGMAX; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_xsu_triclinic(int n) { int j; double **x = atom->x; tagint *image = atom->image; double *boxlo = domain->boxlo; double *h_inv = domain->h_inv; for (int i = 0; i < nchoose; i++) { j = clist[i]; buf[n] = h_inv[0]*(x[j][0]-boxlo[0]) + h_inv[5]*(x[j][1]-boxlo[1]) + h_inv[4]*(x[j][2]-boxlo[2]) + (image[j] & IMGMASK) - IMGMAX; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_ysu_triclinic(int n) { int j; double **x = atom->x; tagint *image = atom->image; double *boxlo = domain->boxlo; double *h_inv = domain->h_inv; for (int i = 0; i < nchoose; i++) { j = clist[i]; buf[n] = h_inv[1]*(x[j][1]-boxlo[1]) + h_inv[3]*(x[j][2]-boxlo[2]) + (image[j] >> IMGBITS & IMGMASK) - IMGMAX; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_zsu_triclinic(int n) { int j; double **x = atom->x; tagint *image = atom->image; double *boxlo = domain->boxlo; double *h_inv = domain->h_inv; for (int i = 0; i < nchoose; i++) { j = clist[i]; buf[n] = h_inv[2]*(x[j][2]-boxlo[2]) + (image[j] >> IMG2BITS) - IMGMAX; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_ix(int n) { tagint *image = atom->image; for (int i = 0; i < nchoose; i++) { buf[n] = (image[clist[i]] & IMGMASK) - IMGMAX; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_iy(int n) { tagint *image = atom->image; for (int i = 0; i < nchoose; i++) { buf[n] = (image[clist[i]] >> IMGBITS & IMGMASK) - IMGMAX; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_iz(int n) { tagint *image = atom->image; for (int i = 0; i < nchoose; i++) { buf[n] = (image[clist[i]] >> IMG2BITS) - IMGMAX; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_vx(int n) { double **v = atom->v; for (int i = 0; i < nchoose; i++) { buf[n] = v[clist[i]][0]; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_vy(int n) { double **v = atom->v; for (int i = 0; i < nchoose; i++) { buf[n] = v[clist[i]][1]; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_vz(int n) { double **v = atom->v; for (int i = 0; i < nchoose; i++) { buf[n] = v[clist[i]][2]; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_fx(int n) { double **f = atom->f; for (int i = 0; i < nchoose; i++) { buf[n] = f[clist[i]][0]; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_fy(int n) { double **f = atom->f; for (int i = 0; i < nchoose; i++) { buf[n] = f[clist[i]][1]; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_fz(int n) { double **f = atom->f; for (int i = 0; i < nchoose; i++) { buf[n] = f[clist[i]][2]; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_q(int n) { double *q = atom->q; for (int i = 0; i < nchoose; i++) { buf[n] = q[clist[i]]; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_mux(int n) { double **mu = atom->mu; for (int i = 0; i < nchoose; i++) { buf[n] = mu[clist[i]][0]; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_muy(int n) { double **mu = atom->mu; for (int i = 0; i < nchoose; i++) { buf[n] = mu[clist[i]][1]; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_muz(int n) { double **mu = atom->mu; for (int i = 0; i < nchoose; i++) { buf[n] = mu[clist[i]][2]; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_mu(int n) { double **mu = atom->mu; for (int i = 0; i < nchoose; i++) { buf[n] = mu[clist[i]][3]; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_radius(int n) { double *radius = atom->radius; for (int i = 0; i < nchoose; i++) { buf[n] = radius[clist[i]]; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_diameter(int n) { double *radius = atom->radius; for (int i = 0; i < nchoose; i++) { buf[n] = 2.0*radius[clist[i]]; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_omegax(int n) { double **omega = atom->omega; for (int i = 0; i < nchoose; i++) { buf[n] = omega[clist[i]][0]; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_omegay(int n) { double **omega = atom->omega; for (int i = 0; i < nchoose; i++) { buf[n] = omega[clist[i]][1]; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_omegaz(int n) { double **omega = atom->omega; for (int i = 0; i < nchoose; i++) { buf[n] = omega[clist[i]][2]; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_angmomx(int n) { double **angmom = atom->angmom; for (int i = 0; i < nchoose; i++) { buf[n] = angmom[clist[i]][0]; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_angmomy(int n) { double **angmom = atom->angmom; for (int i = 0; i < nchoose; i++) { buf[n] = angmom[clist[i]][1]; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_angmomz(int n) { double **angmom = atom->angmom; for (int i = 0; i < nchoose; i++) { buf[n] = angmom[clist[i]][2]; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_tqx(int n) { double **torque = atom->torque; for (int i = 0; i < nchoose; i++) { buf[n] = torque[clist[i]][0]; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_tqy(int n) { double **torque = atom->torque; for (int i = 0; i < nchoose; i++) { buf[n] = torque[clist[i]][1]; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_tqz(int n) { double **torque = atom->torque; for (int i = 0; i < nchoose; i++) { buf[n] = torque[clist[i]][2]; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_spin(int n) { int *spin = atom->spin; for (int i = 0; i < nchoose; i++) { buf[n] = spin[clist[i]]; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_eradius(int n) { double *eradius = atom->eradius; for (int i = 0; i < nchoose; i++) { buf[n] = eradius[clist[i]]; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_ervel(int n) { double *ervel = atom->ervel; for (int i = 0; i < nchoose; i++) { buf[n] = ervel[clist[i]]; n += size_one; } } /* ---------------------------------------------------------------------- */ void DumpCustom::pack_erforce(int n) { double *erforce = atom->erforce; for (int i = 0; i < nchoose; i++) { buf[n] = erforce[clist[i]]; n += size_one; } } diff --git a/src/dump_image.cpp b/src/dump_image.cpp index a0063ff30..51c77158e 100644 --- a/src/dump_image.cpp +++ b/src/dump_image.cpp @@ -1,1068 +1,1070 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ #include "math.h" #include "ctype.h" #include "stdlib.h" #include "string.h" #include "dump_image.h" #include "image.h" #include "atom.h" #include "domain.h" #include "group.h" #include "force.h" #include "comm.h" #include "input.h" #include "variable.h" #include "math_const.h" #include "error.h" #include "memory.h" using namespace LAMMPS_NS; using namespace MathConst; #define BIG 1.0e20 -enum{PPM,JPG,PNG}; enum{NUMERIC,ATOM,TYPE,ELEMENT,ATTRIBUTE}; enum{STATIC,DYNAMIC}; enum{NO,YES}; /* ---------------------------------------------------------------------- */ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) : DumpCustom(lmp, narg, arg) { if (binary || multiproc) error->all(FLERR,"Invalid dump image filename"); // force binary flag on to avoid corrupted output on Windows binary = 1; // set filetype based on filename suffix int n = strlen(filename); if (strlen(filename) > 4 && strcmp(&filename[n-4],".jpg") == 0) filetype = JPG; else if (strlen(filename) > 4 && strcmp(&filename[n-4],".JPG") == 0) filetype = JPG; else if (strlen(filename) > 5 && strcmp(&filename[n-5],".jpeg") == 0) filetype = JPG; else if (strlen(filename) > 5 && strcmp(&filename[n-5],".JPEG") == 0) filetype = JPG; else if (strlen(filename) > 4 && strcmp(&filename[n-4],".png") == 0) filetype = PNG; else if (strlen(filename) > 4 && strcmp(&filename[n-4],".PNG") == 0) filetype = PNG; else filetype = PPM; #ifndef LAMMPS_JPEG if (filetype == JPG) error->all(FLERR,"Support for writing images in JPEG format not included"); #endif #ifndef LAMMPS_PNG if (filetype == PNG) error->all(FLERR,"Support for writing images in PNG format not included"); #endif // atom color,diameter settings if (nfield != 2) error->all(FLERR,"Illegal dump image command"); acolor = ATTRIBUTE; if (strcmp(arg[5],"type") == 0) acolor = TYPE; else if (strcmp(arg[5],"element") == 0) acolor = ELEMENT; adiam = ATTRIBUTE; if (strcmp(arg[6],"type") == 0) adiam = TYPE; else if (strcmp(arg[6],"element") == 0) adiam = ELEMENT; // create Image class with single colormap for atoms // change defaults for 2d image = new Image(lmp,1); if (domain->dimension == 2) { image->theta = 0.0; image->phi = 0.0; image->up[0] = 0.0; image->up[1] = 1.0; image->up[2] = 0.0; } // set defaults for optional args atomflag = YES; if (atom->nbondtypes == 0) bondflag = NO; else { bondflag = YES; bcolor = ATOM; bdiam = NUMERIC; bdiamvalue = 0.5; } thetastr = phistr = NULL; cflag = STATIC; cx = cy = cz = 0.5; cxstr = cystr = czstr = NULL; upxstr = upystr = upzstr = NULL; zoomstr = NULL; perspstr = NULL; boxflag = YES; boxdiam = 0.02; axesflag = NO; // parse optional args int iarg = ioptional; while (iarg < narg) { if (strcmp(arg[iarg],"adiam") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal dump image command"); adiam = NUMERIC; adiamvalue = force->numeric(FLERR,arg[iarg+1]); if (adiamvalue <= 0.0) error->all(FLERR,"Illegal dump image command"); iarg += 2; } else if (strcmp(arg[iarg],"atom") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal dump image command"); if (strcmp(arg[iarg+1],"yes") == 0) atomflag = YES; else if (strcmp(arg[iarg+1],"no") == 0) atomflag = NO; else error->all(FLERR,"Illegal dump image command"); iarg += 2; } else if (strcmp(arg[iarg],"bond") == 0) { if (iarg+3 > narg) error->all(FLERR,"Illegal dump image command"); if (atom->nbondtypes == 0) error->all(FLERR,"Dump image bond not allowed with no bond types"); bondflag = YES; if (strcmp(arg[iarg+1],"none") == 0) bondflag = NO; else if (strcmp(arg[iarg+1],"atom") == 0) bcolor = ATOM; else if (strcmp(arg[iarg+1],"type") == 0) bcolor = TYPE; else error->all(FLERR,"Illegal dump image command"); if (!islower(arg[iarg+2][0])) { bdiam = NUMERIC; bdiamvalue = force->numeric(FLERR,arg[iarg+2]); if (bdiamvalue <= 0.0) error->all(FLERR,"Illegal dump image command"); } else if (strcmp(arg[iarg+2],"atom") == 0) bdiam = ATOM; else if (strcmp(arg[iarg+2],"type") == 0) bdiam = TYPE; else if (strcmp(arg[iarg+2],"none") == 0) bondflag = NO; else error->all(FLERR,"Illegal dump image command"); iarg += 3; } else if (strcmp(arg[iarg],"size") == 0) { if (iarg+3 > narg) error->all(FLERR,"Illegal dump image command"); int width = force->inumeric(FLERR,arg[iarg+1]); int height = force->inumeric(FLERR,arg[iarg+2]); if (width <= 0 || height <= 0) error->all(FLERR,"Illegal dump image command"); image->width = width; image->height = height; iarg += 3; } else if (strcmp(arg[iarg],"view") == 0) { if (iarg+3 > narg) error->all(FLERR,"Illegal dump image command"); if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) { int n = strlen(&arg[iarg+1][2]) + 1; thetastr = new char[n]; strcpy(thetastr,&arg[iarg+1][2]); } else { double theta = force->numeric(FLERR,arg[iarg+1]); if (theta < 0.0 || theta > 180.0) error->all(FLERR,"Invalid dump image theta value"); theta *= MY_PI/180.0; image->theta = theta; } if (strstr(arg[iarg+2],"v_") == arg[iarg+2]) { int n = strlen(&arg[iarg+2][2]) + 1; phistr = new char[n]; strcpy(phistr,&arg[iarg+2][2]); } else { double phi = force->numeric(FLERR,arg[iarg+2]); phi *= MY_PI/180.0; image->phi = phi; } iarg += 3; } else if (strcmp(arg[iarg],"center") == 0) { if (iarg+5 > narg) error->all(FLERR,"Illegal dump image command"); if (strcmp(arg[iarg+1],"s") == 0) cflag = STATIC; else if (strcmp(arg[iarg+1],"d") == 0) cflag = DYNAMIC; else error->all(FLERR,"Illegal dump image command"); if (strstr(arg[iarg+2],"v_") == arg[iarg+2]) { int n = strlen(&arg[iarg+2][2]) + 1; cxstr = new char[n]; strcpy(cxstr,&arg[iarg+2][2]); cflag = DYNAMIC; } else cx = force->numeric(FLERR,arg[iarg+2]); if (strstr(arg[iarg+3],"v_") == arg[iarg+3]) { int n = strlen(&arg[iarg+3][2]) + 1; cystr = new char[n]; strcpy(cystr,&arg[iarg+3][2]); cflag = DYNAMIC; } else cy = force->numeric(FLERR,arg[iarg+3]); if (strstr(arg[iarg+4],"v_") == arg[iarg+4]) { int n = strlen(&arg[iarg+4][2]) + 1; czstr = new char[n]; strcpy(czstr,&arg[iarg+4][2]); cflag = DYNAMIC; } else cz = force->numeric(FLERR,arg[iarg+4]); iarg += 5; } else if (strcmp(arg[iarg],"up") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal dump image command"); if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) { int n = strlen(&arg[iarg+1][2]) + 1; upxstr = new char[n]; strcpy(upxstr,&arg[iarg+1][2]); } else image->up[0] = force->numeric(FLERR,arg[iarg+1]); if (strstr(arg[iarg+2],"v_") == arg[iarg+2]) { int n = strlen(&arg[iarg+2][2]) + 1; upystr = new char[n]; strcpy(upystr,&arg[iarg+2][2]); } else image->up[1] = force->numeric(FLERR,arg[iarg+2]); if (strstr(arg[iarg+3],"v_") == arg[iarg+3]) { int n = strlen(&arg[iarg+3][2]) + 1; upzstr = new char[n]; strcpy(upzstr,&arg[iarg+3][2]); } else image->up[2] = force->numeric(FLERR,arg[iarg+3]); iarg += 4; } else if (strcmp(arg[iarg],"zoom") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal dump image command"); if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) { int n = strlen(&arg[iarg+1][2]) + 1; zoomstr = new char[n]; strcpy(zoomstr,&arg[iarg+1][2]); } else { double zoom = force->numeric(FLERR,arg[iarg+1]); if (zoom <= 0.0) error->all(FLERR,"Illegal dump image command"); image->zoom = zoom; } iarg += 2; } else if (strcmp(arg[iarg],"persp") == 0) { error->all(FLERR,"Dump image persp option is not yet supported"); if (iarg+2 > narg) error->all(FLERR,"Illegal dump image command"); if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) { int n = strlen(&arg[iarg+1][2]) + 1; perspstr = new char[n]; strcpy(perspstr,&arg[iarg+1][2]); } else { double persp = force->numeric(FLERR,arg[iarg+1]); if (persp < 0.0) error->all(FLERR,"Illegal dump image command"); image->persp = persp; } iarg += 2; } else if (strcmp(arg[iarg],"box") == 0) { if (iarg+3 > narg) error->all(FLERR,"Illegal dump image command"); if (strcmp(arg[iarg+1],"yes") == 0) boxflag = YES; else if (strcmp(arg[iarg+1],"no") == 0) boxflag = NO; else error->all(FLERR,"Illegal dump image command"); boxdiam = force->numeric(FLERR,arg[iarg+2]); if (boxdiam < 0.0) error->all(FLERR,"Illegal dump image command"); iarg += 3; } else if (strcmp(arg[iarg],"axes") == 0) { if (iarg+3 > narg) error->all(FLERR,"Illegal dump image command"); if (strcmp(arg[iarg+1],"yes") == 0) axesflag = YES; else if (strcmp(arg[iarg+1],"no") == 0) axesflag = NO; else error->all(FLERR,"Illegal dump image command"); axeslen = force->numeric(FLERR,arg[iarg+2]); axesdiam = force->numeric(FLERR,arg[iarg+3]); if (axeslen < 0.0 || axesdiam < 0.0) error->all(FLERR,"Illegal dump image command"); iarg += 4; } else if (strcmp(arg[iarg],"shiny") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal dump image command"); double shiny = force->numeric(FLERR,arg[iarg+1]); if (shiny < 0.0 || shiny > 1.0) error->all(FLERR,"Illegal dump image command"); image->shiny = shiny; iarg += 2; } else if (strcmp(arg[iarg],"ssao") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal dump image command"); if (strcmp(arg[iarg+1],"yes") == 0) image->ssao = YES; else if (strcmp(arg[iarg+1],"no") == 0) image->ssao = NO; else error->all(FLERR,"Illegal dump image command"); int seed = force->inumeric(FLERR,arg[iarg+2]); if (seed <= 0) error->all(FLERR,"Illegal dump image command"); image->seed = seed; double ssaoint = force->numeric(FLERR,arg[iarg+3]); if (ssaoint < 0.0 || ssaoint > 1.0) error->all(FLERR,"Illegal dump image command"); image->ssaoint = ssaoint; iarg += 4; } else error->all(FLERR,"Illegal dump image command"); } // allocate image buffer now that image size is known image->buffers(); // communication neede for bonds colored by atoms if (bondflag) { if (bcolor == ATOM || bdiam == ATOM) comm_forward = 3; else comm_forward = 1; } // additional defaults for dump_modify options diamtype = new double[ntypes+1]; diamelement = new double[ntypes+1]; colortype = new double*[ntypes+1]; colorelement = new double*[ntypes+1]; for (int i = 1; i <= ntypes; i++) { diamtype[i] = 1.0; if (i % 6 == 1) colortype[i] = image->color2rgb("red"); else if (i % 6 == 2) colortype[i] = image->color2rgb("green"); else if (i % 6 == 3) colortype[i] = image->color2rgb("blue"); else if (i % 6 == 4) colortype[i] = image->color2rgb("yellow"); else if (i % 6 == 5) colortype[i] = image->color2rgb("aqua"); else if (i % 6 == 0) colortype[i] = image->color2rgb("cyan"); } if (bondflag) { bdiamtype = new double[atom->nbondtypes+1]; bcolortype = new double*[atom->nbondtypes+1]; for (int i = 1; i <= atom->nbondtypes; i++) { bdiamtype[i] = 0.5; if (i % 6 == 1) bcolortype[i] = image->color2rgb("red"); else if (i % 6 == 2) bcolortype[i] = image->color2rgb("green"); else if (i % 6 == 3) bcolortype[i] = image->color2rgb("blue"); else if (i % 6 == 4) bcolortype[i] = image->color2rgb("yellow"); else if (i % 6 == 5) bcolortype[i] = image->color2rgb("aqua"); else if (i % 6 == 0) bcolortype[i] = image->color2rgb("cyan"); } } else { bdiamtype = NULL; bcolortype = NULL; } // viewflag = DYNAMIC if any view parameter is dynamic viewflag = STATIC; if (thetastr || phistr || cflag == DYNAMIC || upxstr || upystr || upzstr || zoomstr || perspstr) viewflag = DYNAMIC; box_bounds(); if (cflag == STATIC) box_center(); if (viewflag == STATIC) view_params(); // local data maxbufcopy = 0; chooseghost = NULL; bufcopy = NULL; } /* ---------------------------------------------------------------------- */ DumpImage::~DumpImage() { delete image; delete [] diamtype; delete [] diamelement; delete [] colortype; delete [] colorelement; delete [] bdiamtype; delete [] bcolortype; memory->destroy(chooseghost); memory->destroy(bufcopy); } /* ---------------------------------------------------------------------- */ void DumpImage::init_style() { if (multifile == 0) error->all(FLERR,"Dump image requires one snapshot per file"); if (sort_flag) error->all(FLERR,"Dump image cannot perform sorting"); DumpCustom::init_style(); // check variables if (thetastr) { thetavar = input->variable->find(thetastr); if (thetavar < 0) error->all(FLERR,"Variable name for dump image theta does not exist"); if (!input->variable->equalstyle(thetavar)) error->all(FLERR,"Variable for dump image theta is invalid style"); } if (phistr) { phivar = input->variable->find(phistr); if (phivar < 0) error->all(FLERR,"Variable name for dump image phi does not exist"); if (!input->variable->equalstyle(phivar)) error->all(FLERR,"Variable for dump image phi is invalid style"); } if (cxstr) { cxvar = input->variable->find(cxstr); if (cxvar < 0) error->all(FLERR,"Variable name for dump image center does not exist"); if (!input->variable->equalstyle(cxvar)) error->all(FLERR,"Variable for dump image center is invalid style"); } if (cystr) { cyvar = input->variable->find(cystr); if (cyvar < 0) error->all(FLERR,"Variable name for dump image center does not exist"); if (!input->variable->equalstyle(cyvar)) error->all(FLERR,"Variable for dump image center is invalid style"); } if (czstr) { czvar = input->variable->find(czstr); if (czvar < 0) error->all(FLERR,"Variable name for dump image center does not exist"); if (!input->variable->equalstyle(czvar)) error->all(FLERR,"Variable for dump image center is invalid style"); } if (upxstr) { upxvar = input->variable->find(upxstr); if (upxvar < 0) error->all(FLERR,"Variable name for dump image center does not exist"); if (!input->variable->equalstyle(upxvar)) error->all(FLERR,"Variable for dump image center is invalid style"); } if (upystr) { upyvar = input->variable->find(upystr); if (upyvar < 0) error->all(FLERR,"Variable name for dump image center does not exist"); if (!input->variable->equalstyle(upyvar)) error->all(FLERR,"Variable for dump image center is invalid style"); } if (upzstr) { upzvar = input->variable->find(upzstr); if (upzvar < 0) error->all(FLERR,"Variable name for dump image center does not exist"); if (!input->variable->equalstyle(upzvar)) error->all(FLERR,"Variable for dump image center is invalid style"); } if (zoomstr) { zoomvar = input->variable->find(zoomstr); if (zoomvar < 0) error->all(FLERR,"Variable name for dump image zoom does not exist"); if (!input->variable->equalstyle(zoomvar)) error->all(FLERR,"Variable for dump image zoom is invalid style"); } if (perspstr) { perspvar = input->variable->find(perspstr); if (perspvar < 0) error->all(FLERR,"Variable name for dump image persp does not exist"); if (!input->variable->equalstyle(perspvar)) error->all(FLERR,"Variable for dump image persp is invalid style"); } // set up type -> element mapping if (atomflag && acolor == ELEMENT) { for (int i = 1; i <= ntypes; i++) { colorelement[i] = image->element2color(typenames[i]); if (colorelement[i] == NULL) error->all(FLERR,"Invalid dump image element name"); } } if (atomflag && adiam == ELEMENT) { for (int i = 1; i <= ntypes; i++) { diamelement[i] = image->element2diam(typenames[i]); if (diamelement[i] == 0.0) error->all(FLERR,"Invalid dump image element name"); } } } /* ---------------------------------------------------------------------- */ void DumpImage::write() { // open new file openfile(); // reset box center and view parameters if dynamic box_bounds(); if (cflag == DYNAMIC) box_center(); if (viewflag == DYNAMIC) view_params(); // nme = # of atoms this proc will contribute to dump nme = count(); if (nme > maxbuf) { maxbuf = nme; memory->destroy(buf); memory->create(buf,maxbuf*size_one,"dump:buf"); } // pack buf with color & diameter pack(NULL); // set minmax color range if using dynamic atom color map if (acolor == ATTRIBUTE && image->map_dynamic(0)) { double two[2],twoall[2]; double lo = BIG; double hi = -BIG; int m = 0; for (int i = 0; i < nchoose; i++) { lo = MIN(lo,buf[m]); hi = MAX(hi,buf[m]); m += size_one; } two[0] = -lo; two[1] = hi; MPI_Allreduce(two,twoall,2,MPI_DOUBLE,MPI_MAX,world); int flag = image->map_minmax(0,-twoall[0],twoall[1]); if (flag) error->all(FLERR,"Invalid color map min/max values"); } // create image on each proc, then merge them image->clear(); create_image(); image->merge(); // write image file if (me == 0) { if (filetype == JPG) image->write_JPG(fp); else if (filetype == PNG) image->write_PNG(fp); else image->write_PPM(fp); - fclose(fp); + if (multifile) { + fclose(fp); + fp = NULL; + } } } /* ---------------------------------------------------------------------- simulation box bounds ------------------------------------------------------------------------- */ void DumpImage::box_bounds() { if (domain->triclinic == 0) { boxxlo = domain->boxlo[0]; boxxhi = domain->boxhi[0]; boxylo = domain->boxlo[1]; boxyhi = domain->boxhi[1]; boxzlo = domain->boxlo[2]; boxzhi = domain->boxhi[2]; } else { boxxlo = domain->boxlo_bound[0]; boxxhi = domain->boxhi_bound[0]; boxylo = domain->boxlo_bound[1]; boxyhi = domain->boxhi_bound[1]; boxzlo = domain->boxlo_bound[2]; boxzhi = domain->boxhi_bound[2]; boxxy = domain->xy; boxxz = domain->xz; boxyz = domain->yz; } } /* ---------------------------------------------------------------------- reset view parameters called once from constructor if view is STATIC called every snapshot from write() if view is DYNAMIC ------------------------------------------------------------------------- */ void DumpImage::box_center() { if (cxstr) cx = input->variable->compute_equal(cxvar); if (cystr) cy = input->variable->compute_equal(cyvar); if (czstr) cz = input->variable->compute_equal(czvar); image->xctr = boxxlo + cx*(boxxhi-boxxlo); image->yctr = boxylo + cy*(boxyhi-boxylo); image->zctr = boxzlo + cz*(boxzhi-boxzlo); } /* ---------------------------------------------------------------------- reset view parameters in Image class called once from constructor if view is STATIC called every snapshot from write() if view is DYNAMIC ------------------------------------------------------------------------- */ void DumpImage::view_params() { // view direction theta and phi if (thetastr) { double theta = input->variable->compute_equal(thetavar); if (theta < 0.0 || theta > 180.0) error->all(FLERR,"Invalid dump image theta value"); theta *= MY_PI/180.0; image->theta = theta; } if (phistr) { double phi = input->variable->compute_equal(phivar); phi *= MY_PI/180.0; image->phi = phi; } // up vector if (upxstr) image->up[0] = input->variable->compute_equal(upxvar); if (upystr) image->up[1] = input->variable->compute_equal(upyvar); if (upzstr) image->up[2] = input->variable->compute_equal(upzvar); // zoom and perspective if (zoomstr) image->zoom = input->variable->compute_equal(zoomvar); if (image->zoom <= 0.0) error->all(FLERR,"Invalid dump image zoom value"); if (perspstr) image->persp = input->variable->compute_equal(perspvar); if (image->persp < 0.0) error->all(FLERR,"Invalid dump image persp value"); // remainder of view setup is internal to Image class image->view_params(boxxlo,boxxhi,boxylo,boxyhi,boxzlo,boxzhi); } /* ---------------------------------------------------------------------- create image for atoms on this proc every pixel has depth ------------------------------------------------------------------------- */ void DumpImage::create_image() { int i,j,m,itype,atom1,atom2; double diameter,delx,dely,delz; double *color,*color1,*color2; double xmid[3]; // render my atoms if (atomflag) { double **x = atom->x; m = 0; for (i = 0; i < nchoose; i++) { j = clist[i]; if (acolor == TYPE) { itype = static_cast<int> (buf[m]); color = colortype[itype]; } else if (acolor == ELEMENT) { itype = static_cast<int> (buf[m]); color = colorelement[itype]; } else if (acolor == ATTRIBUTE) { color = image->map_value2color(0,buf[m]); } if (adiam == NUMERIC) { diameter = adiamvalue; } else if (adiam == TYPE) { itype = static_cast<int> (buf[m+1]); diameter = diamtype[itype]; } else if (adiam == ELEMENT) { itype = static_cast<int> (buf[m+1]); diameter = diamelement[itype]; } else if (adiam == ATTRIBUTE) { diameter = buf[m+1]; } image->draw_sphere(x[j],color,diameter); m += size_one; } } // render bonds for my atoms // both atoms in bond must be selected for bond to be rendered // if newton_bond is off, only render bond once // render bond in 2 pieces if crosses periodic boundary // if bond is deleted (type = 0), do not render // if bond is turned off (type < 0), still render if (bondflag) { double **x = atom->x; int *tag = atom->tag; int **bond_atom = atom->bond_atom; int **bond_type = atom->bond_type; int *num_bond = atom->num_bond; int *type = atom->type; int nlocal = atom->nlocal; int nall = atom->nlocal + atom->nghost; int newton_bond = force->newton_bond; // communicate choose flag for ghost atoms to know if they are selected // if bcolor/bdiam = ATOM, setup bufcopy to comm atom color/diam attributes if (nall > maxbufcopy) { maxbufcopy = atom->nmax; memory->destroy(chooseghost); memory->create(chooseghost,maxbufcopy,"dump:chooseghost"); if (comm_forward == 3) { memory->destroy(bufcopy); memory->create(bufcopy,maxbufcopy,2,"dump:bufcopy"); } } for (i = 0; i < nlocal; i++) chooseghost[i] = choose[i]; if (comm_forward == 3) { for (i = 0; i < nlocal; i++) bufcopy[i][0] = bufcopy[i][1] = 0.0; m = 0; for (i = 0; i < nchoose; i++) { j = clist[i]; bufcopy[j][0] = buf[m]; bufcopy[j][1] = buf[m+1]; m += size_one; } } comm->forward_comm_dump(this); for (i = 0; i < nchoose; i++) { atom1 = clist[i]; for (m = 0; m < num_bond[atom1]; m++) { atom2 = atom->map(bond_atom[atom1][m]); if (atom2 < 0 || !chooseghost[atom2]) continue; if (newton_bond == 0 && tag[atom1] > tag[atom2]) continue; if (bond_type[atom1][m] == 0) continue; if (bcolor == ATOM) { if (acolor == TYPE) { color1 = colortype[type[atom1]]; color2 = colortype[type[atom2]]; } else if (acolor == ELEMENT) { color1 = colorelement[type[atom1]]; color2 = colorelement[type[atom2]]; } else if (acolor == ATTRIBUTE) { color1 = image->map_value2color(0,bufcopy[atom1][0]); color2 = image->map_value2color(0,bufcopy[atom2][0]); } } else if (bcolor == TYPE) { itype = bond_type[atom1][m]; if (itype < 0) itype = -itype; color = bcolortype[itype]; } if (bdiam == NUMERIC) { diameter = bdiamvalue; } else if (bdiam == ATOM) { if (adiam == NUMERIC) { diameter = adiamvalue; } else if (adiam == TYPE) { diameter = MIN(diamtype[type[atom1]],diamtype[type[atom1]]); } else if (adiam == ELEMENT) { diameter = MIN(diamelement[type[atom1]],diamelement[type[atom1]]); } else if (adiam == ATTRIBUTE) { diameter = MIN(bufcopy[atom1][1],bufcopy[atom2][1]); } } else if (bdiam == TYPE) { itype = bond_type[atom1][m]; if (itype < 0) itype = -itype; diameter = bdiamtype[itype]; } // draw cylinder in 2 pieces if bcolor = ATOM // or bond crosses periodic boundary delx = x[atom2][0] - x[atom1][0]; dely = x[atom2][1] - x[atom1][1]; delz = x[atom2][2] - x[atom1][2]; if (bcolor == ATOM || domain->minimum_image_check(delx,dely,delz)) { domain->minimum_image(delx,dely,delz); xmid[0] = x[atom1][0] + 0.5*delx; xmid[1] = x[atom1][1] + 0.5*dely; xmid[2] = x[atom1][2] + 0.5*delz; if (bcolor == ATOM) image->draw_cylinder(x[atom1],xmid,color1,diameter,3); else image->draw_cylinder(x[atom1],xmid,color,diameter,3); xmid[0] = x[atom2][0] - 0.5*delx; xmid[1] = x[atom2][1] - 0.5*dely; xmid[2] = x[atom2][2] - 0.5*delz; if (bcolor == ATOM) image->draw_cylinder(xmid,x[atom2],color2,diameter,3); else image->draw_cylinder(xmid,x[atom2],color,diameter,3); } else image->draw_cylinder(x[atom1],x[atom2],color,diameter,3); } } } // render outline of simulation box, orthogonal or triclinic if (boxflag) { double diameter = MIN(boxxhi-boxxlo,boxyhi-boxylo); if (domain->dimension == 3) diameter = MIN(diameter,boxzhi-boxzlo); diameter *= boxdiam; double (*boxcorners)[3]; double box[8][3]; if (domain->triclinic == 0) { box[0][0] = boxxlo; box[0][1] = boxylo; box[0][2] = boxzlo; box[1][0] = boxxhi; box[1][1] = boxylo; box[1][2] = boxzlo; box[2][0] = boxxlo; box[2][1] = boxyhi; box[2][2] = boxzlo; box[3][0] = boxxhi; box[3][1] = boxyhi; box[3][2] = boxzlo; box[4][0] = boxxlo; box[4][1] = boxylo; box[4][2] = boxzhi; box[5][0] = boxxhi; box[5][1] = boxylo; box[5][2] = boxzhi; box[6][0] = boxxlo; box[6][1] = boxyhi; box[6][2] = boxzhi; box[7][0] = boxxhi; box[7][1] = boxyhi; box[7][2] = boxzhi; boxcorners = box; } else { domain->box_corners(); boxcorners = domain->corners; } image->draw_box(boxcorners,diameter); } // render XYZ axes in red/green/blue // offset by 10% of box size and scale by axeslen if (axesflag) { double diameter = MIN(boxxhi-boxxlo,boxyhi-boxylo); if (domain->dimension == 3) diameter = MIN(diameter,boxzhi-boxzlo); diameter *= axesdiam; double (*boxcorners)[3]; double axes[4][3]; if (domain->triclinic == 0) { axes[0][0] = boxxlo; axes[0][1] = boxylo; axes[0][2] = boxzlo; axes[1][0] = boxxhi; axes[1][1] = boxylo; axes[1][2] = boxzlo; axes[2][0] = boxxlo; axes[2][1] = boxyhi; axes[2][2] = boxzlo; axes[3][0] = boxxlo; axes[3][1] = boxylo; axes[3][2] = boxzhi; } else { domain->box_corners(); boxcorners = domain->corners; axes[0][0] = boxcorners[0][0]; axes[0][1] = boxcorners[0][1]; axes[0][2] = boxcorners[0][2]; axes[1][0] = boxcorners[1][0]; axes[1][1] = boxcorners[1][1]; axes[1][2] = boxcorners[1][2]; axes[2][0] = boxcorners[2][0]; axes[2][1] = boxcorners[2][1]; axes[2][2] = boxcorners[2][2]; axes[3][0] = boxcorners[4][0]; axes[3][1] = boxcorners[4][1]; axes[3][2] = boxcorners[4][2]; } double offset = MAX(boxxhi-boxxlo,boxyhi-boxylo); if (domain->dimension == 3) offset = MAX(offset,boxzhi-boxzlo); offset *= 0.1; axes[0][0] -= offset; axes[0][1] -= offset; axes[0][2] -= offset; axes[1][0] -= offset; axes[1][1] -= offset; axes[1][2] -= offset; axes[2][0] -= offset; axes[2][1] -= offset; axes[2][2] -= offset; axes[3][0] -= offset; axes[3][1] -= offset; axes[3][2] -= offset; axes[1][0] = axes[0][0] + axeslen*(axes[1][0]-axes[0][0]); axes[1][1] = axes[0][1] + axeslen*(axes[1][1]-axes[0][1]); axes[1][2] = axes[0][2] + axeslen*(axes[1][2]-axes[0][2]); axes[2][0] = axes[0][0] + axeslen*(axes[2][0]-axes[0][0]); axes[2][1] = axes[0][1] + axeslen*(axes[2][1]-axes[0][1]); axes[2][2] = axes[0][2] + axeslen*(axes[2][2]-axes[0][2]); axes[3][0] = axes[0][0] + axeslen*(axes[3][0]-axes[0][0]); axes[3][1] = axes[0][1] + axeslen*(axes[3][1]-axes[0][1]); axes[3][2] = axes[0][2] + axeslen*(axes[3][2]-axes[0][2]); image->draw_axes(axes,diameter); } } /* ---------------------------------------------------------------------- */ int DumpImage::pack_comm(int n, int *list, double *buf, int pbc_flag, int *pbc) { int i,j,m; m = 0; if (comm_forward == 1) { for (i = 0; i < n; i++) { j = list[i]; buf[m++] = chooseghost[j]; } } else { for (i = 0; i < n; i++) { j = list[i]; buf[m++] = chooseghost[j]; buf[m++] = bufcopy[j][0]; buf[m++] = bufcopy[j][1]; } } return comm_forward; } /* ---------------------------------------------------------------------- */ void DumpImage::unpack_comm(int n, int first, double *buf) { int i,m,last; m = 0; last = first + n; if (comm_forward == 1) for (i = first; i < last; i++) chooseghost[i] = static_cast<int> (buf[m++]); else { for (i = first; i < last; i++) { chooseghost[i] = static_cast<int> (buf[m++]); bufcopy[i][0] = buf[m++]; bufcopy[i][1] = buf[m++]; } } } /* ---------------------------------------------------------------------- */ int DumpImage::modify_param(int narg, char **arg) { int n = DumpCustom::modify_param(narg,arg); if (n) return n; if (strcmp(arg[0],"acolor") == 0) { if (narg < 3) error->all(FLERR,"Illegal dump_modify command"); int nlo,nhi; force->bounds(arg[1],atom->ntypes,nlo,nhi); // ptrs = list of ncount colornames separated by '/' int ncount = 1; char *nextptr; char *ptr = arg[2]; while (nextptr = strchr(ptr,'/')) { ptr = nextptr + 1; ncount++; } char **ptrs = new char*[ncount+1]; ncount = 0; ptrs[ncount++] = strtok(arg[2],"/"); while (ptrs[ncount++] = strtok(NULL,"/")); ncount--; // assign each of ncount colors in round-robin fashion to types int m = 0; for (int i = nlo; i <= nhi; i++) { colortype[i] = image->color2rgb(ptrs[m%ncount]); if (colortype[i] == NULL) error->all(FLERR,"Invalid color in dump_modify command"); m++; } delete [] ptrs; return 3; } if (strcmp(arg[0],"adiam") == 0) { if (narg < 3) error->all(FLERR,"Illegal dump_modify command"); int nlo,nhi; force->bounds(arg[1],atom->ntypes,nlo,nhi); double diam = force->numeric(FLERR,arg[2]); if (diam <= 0.0) error->all(FLERR,"Illegal dump_modify command"); for (int i = nlo; i <= nhi; i++) diamtype[i] = diam; return 3; } if (strcmp(arg[0],"amap") == 0) { if (narg < 6) error->all(FLERR,"Illegal dump_modify command"); if (strlen(arg[3]) != 2) error->all(FLERR,"Illegal dump_modify command"); int factor = 2; if (arg[3][0] == 's') factor = 1; int nentry = force->inumeric(FLERR,arg[5]); if (nentry < 1) error->all(FLERR,"Illegal dump_modify command"); int n = 6 + factor*nentry; if (narg < n) error->all(FLERR,"Illegal dump_modify command"); int flag = image->map_reset(0,n-1,&arg[1]); if (flag) error->all(FLERR,"Illegal dump_modify command"); return n; } if (strcmp(arg[0],"bcolor") == 0) { if (narg < 3) error->all(FLERR,"Illegal dump_modify command"); if (atom->nbondtypes == 0) error->all(FLERR,"Dump modify bcolor not allowed with no bond types"); int nlo,nhi; force->bounds(arg[1],atom->nbondtypes,nlo,nhi); // ptrs = list of ncount colornames separated by '/' int ncount = 1; char *nextptr; char *ptr = arg[2]; while (nextptr = strchr(ptr,'/')) { ptr = nextptr + 1; ncount++; } char **ptrs = new char*[ncount+1]; ncount = 0; ptrs[ncount++] = strtok(arg[2],"/"); while (ptrs[ncount++] = strtok(NULL,"/")); ncount--; // assign each of ncount colors in round-robin fashion to types int m = 0; for (int i = nlo; i <= nhi; i++) { bcolortype[i] = image->color2rgb(ptrs[m%ncount]); if (bcolortype[i] == NULL) error->all(FLERR,"Invalid color in dump_modify command"); m++; } delete [] ptrs; return 3; } if (strcmp(arg[0],"bdiam") == 0) { if (narg < 3) error->all(FLERR,"Illegal dump_modify command"); if (atom->nbondtypes == 0) error->all(FLERR,"Dump modify bdiam not allowed with no bond types"); int nlo,nhi; force->bounds(arg[1],atom->ntypes,nlo,nhi); double diam = force->numeric(FLERR,arg[2]); if (diam <= 0.0) error->all(FLERR,"Illegal dump_modify command"); for (int i = nlo; i <= nhi; i++) bdiamtype[i] = diam; return 3; } if (strcmp(arg[0],"backcolor") == 0) { if (narg < 2) error->all(FLERR,"Illegal dump_modify command"); double *color = image->color2rgb(arg[1]); if (color == NULL) error->all(FLERR,"Invalid color in dump_modify command"); image->background[0] = static_cast<int> (color[0]*255.0); image->background[1] = static_cast<int> (color[1]*255.0); image->background[2] = static_cast<int> (color[2]*255.0); return 2; } if (strcmp(arg[0],"boxcolor") == 0) { if (narg < 2) error->all(FLERR,"Illegal dump_modify command"); image->boxcolor = image->color2rgb(arg[1]); if (image->boxcolor == NULL) error->all(FLERR,"Invalid color in dump_modify command"); return 2; } if (strcmp(arg[0],"color") == 0) { if (narg < 5) error->all(FLERR,"Illegal dump_modify command"); int flag = image->addcolor(arg[1],force->numeric(FLERR,arg[2]),force->numeric(FLERR,arg[3]),force->numeric(FLERR,arg[4])); if (flag) error->all(FLERR,"Illegal dump_modify command"); return 5; } return 0; } diff --git a/src/dump_image.h b/src/dump_image.h index ab97d7f4a..9c42c584c 100644 --- a/src/dump_image.h +++ b/src/dump_image.h @@ -1,188 +1,190 @@ /* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ #ifdef DUMP_CLASS DumpStyle(image,DumpImage) #else #ifndef LMP_DUMP_IMAGE_H #define LMP_DUMP_IMAGE_H #include "dump_custom.h" namespace LAMMPS_NS { class DumpImage : public DumpCustom { public: DumpImage(class LAMMPS *, int, char**); - ~DumpImage(); + virtual ~DumpImage(); int pack_comm(int, int *, double *, int, int *); void unpack_comm(int, int, double *); - private: + protected: int filetype; + enum{PPM,JPG,PNG}; + int acolor,adiam; // what determines color/diam of atoms double adiamvalue; // atom diameter value int atomflag,bondflag; // 0/1 for draw atoms,bonds int bcolor,bdiam; // what determines color/diam of bonds double bdiamvalue; // bond diameter value char *thetastr,*phistr; // variables for view theta,phi int thetavar,phivar; // index to theta,phi vars int cflag; // static/dynamic box center double cx,cy,cz; // fractional box center char *cxstr,*cystr,*czstr; // variables for box center int cxvar,cyvar,czvar; // index to box center vars char *upxstr,*upystr,*upzstr; // view up vector variables int upxvar,upyvar,upzvar; // index to up vector vars char *zoomstr,*perspstr; // view zoom and perspective variables int zoomvar,perspvar; // index to zoom,persp vars int boxflag,axesflag; // 0/1 for draw box and axes double boxdiam,axeslen,axesdiam; // params for drawing box and axes int viewflag; // overall view is static or dynamic double *diamtype,*diamelement,*bdiamtype; // per-type diameters double **colortype,**colorelement,**bcolortype; // per-type colors class Image *image; // class that renders each image int *chooseghost; // extended choose array for comm double **bufcopy; // buffer for communicating bond/atom info int maxbufcopy; - void init_style(); + virtual void init_style(); int modify_param(int, char **); void write(); void box_center(); void view_params(); void box_bounds(); void create_image(); }; } #endif #endif /* ERROR/WARNING messages: E: Invalid dump image filename The file produced by dump image cannot be binary and must be for a single processor. E: Support for writing images in JPEG format not included LAMMPS was not built with the -DLAMMPS_JPEG switch in the Makefile. E: Support for writing images in PNG format not included LAMMPS was not built with the -DLAMMPS_PNG switch in the Makefile. E: Illegal ... command Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. E: Dump image bond not allowed with no bond types Self-explanatory. E: Invalid dump image theta value Theta must be between 0.0 and 180.0 inclusive. E: Dump image persp option is not yet supported Self-explanatory. E: Dump image requires one snapshot per file Use a "*" in the filename. E: Dump image cannot perform sorting Self-explanatory. E: Variable name for dump image theta does not exist Self-explanatory. E: Variable for dump image theta is invalid style Must be an equal-style variable. E: Variable name for dump image phi does not exist Self-explanatory. E: Variable for dump image phi is invalid style Must be an equal-style variable. E: Variable name for dump image center does not exist Self-explanatory. E: Variable for dump image center is invalid style Must be an equal-style variable. E: Variable name for dump image zoom does not exist Self-explanatory. E: Variable for dump image zoom is invalid style Must be an equal-style variable. E: Variable name for dump image persp does not exist Self-explanatory. E: Variable for dump image persp is invalid style Must be an equal-style variable. E: Invalid dump image element name The specified element name was not in the standard list of elements. See the dump_modify doc page. E: Invalid dump image zoom value Zoom value must be > 0.0. E: Invalid dump image persp value Persp value must be >= 0.0. E: Invalid color in dump_modify command The specified color name was not in the list of recognized colors. See the dump_modify doc page. E: Dump modify bcolor not allowed with no bond types Self-explanatory. E: Dump modify bdiam not allowed with no bond types Self-explanatory. */ diff --git a/src/dump_movie.cpp b/src/dump_movie.cpp new file mode 100644 index 000000000..0e6fc1952 --- /dev/null +++ b/src/dump_movie.cpp @@ -0,0 +1,59 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include <stdio.h> +#include <stdlib.h> +#include <string.h> + +#include "dump_movie.h" +#include "comm.h" +#include "error.h" +#include "memory.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +DumpMovie::DumpMovie(LAMMPS *lmp, int narg, char **arg) : + DumpImage(lmp, narg, arg) +{ + if (multiproc || compressed || multifile) + error->all(FLERR,"Invalid dump movie filename"); + + filetype = PPM; + bitrate = 2000; + framerate = 24; + fp = NULL; +} + +/* ---------------------------------------------------------------------- */ + +void DumpMovie::openfile() +{ + char moviecmd[1024]; + + if ((comm->me == 0) && (fp == NULL)) { + sprintf(moviecmd,"ffmpeg -y -r %d -f image2pipe -c:v ppm -i - " + "-r 24 -b:v %dk %s 2> /dev/null ", framerate, bitrate, filename); + fp = popen(moviecmd,"w"); + } +} +/* ---------------------------------------------------------------------- */ + +void DumpMovie::init_style() +{ + // initialize image style circumventing multifile check + multifile = 1; + DumpImage::init_style(); + multifile = 0; +} diff --git a/src/dump_movie.h b/src/dump_movie.h new file mode 100644 index 000000000..0ed5c577b --- /dev/null +++ b/src/dump_movie.h @@ -0,0 +1,47 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef DUMP_CLASS + +DumpStyle(movie,DumpMovie) + +#else + +#ifndef LMP_DUMP_MOVIE_H +#define LMP_DUMP_MOVIE_H + +#include "dump_image.h" + +namespace LAMMPS_NS { + +class DumpMovie : public DumpImage { + public: + DumpMovie(LAMMPS *, int, char**); + + virtual void openfile(); + virtual void init_style(); + + protected: + int bitrate; // bitrate of video file in kbps + int framerate; // frame rate of animation +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + + +*/