diff --git a/src/dump_cfg.cpp b/src/dump_cfg.cpp index 471b5d3bb0..de7790bb89 100755 --- a/src/dump_cfg.cpp +++ b/src/dump_cfg.cpp @@ -30,6 +30,8 @@ #include "memory.h" #include "error.h" +#define UNWRAPEXPAND 10.0 + using namespace LAMMPS_NS; enum{INT,DOUBLE}; // same as in dump_custom.cpp @@ -41,10 +43,20 @@ DumpCFG::DumpCFG(LAMMPS *lmp, int narg, char **arg) : { if (narg < 10 || strcmp(arg[5],"id") != 0 || strcmp(arg[6],"type") != 0 || - strcmp(arg[7],"xs") != 0 || strcmp(arg[8],"ys") != 0 || - strcmp(arg[9],"zs") != 0) - error->all("Dump cfg arguments must start with 'id type xs ys zs'"); + (strcmp(arg[7],"xs") != 0 && strcmp(arg[7],"xsu") != 0) || + (strcmp(arg[8],"ys") != 0 && strcmp(arg[8],"ysu") != 0) || + (strcmp(arg[9],"zs") != 0 && strcmp(arg[9],"zsu") != 0) + ) + error->all("Dump cfg arguments must start with 'id type xs ys zs' or 'id type xsu ysu zsu'"); + if (strcmp(arg[7],"xs") == 0) + if (strcmp(arg[8],"ysu") == 0 || strcmp(arg[9],"zsu") == 0) + error->all("Dump cfg arguments can not mix xs|ys|zs with xsu|ysu|zsu"); + else unwrapflag = 0; + else if (strcmp(arg[8],"ys") == 0 || strcmp(arg[9],"zs") == 0) + error->all("Dump cfg arguments can not mix xs|ys|zs with xsu|ysu|zsu"); + else unwrapflag = 1; + ntypes = atom->ntypes; typenames = NULL; @@ -189,7 +201,9 @@ void DumpCFG::write_header(bigint n) // special handling for atom style peri // use average volume of particles to scale particles to mimic C atoms // scale box dimension to sc lattice for C with sigma = 1.44 Angstroms - + + // Special handling for unwrapped coordinates + double scale; if (atom->peri_flag) { int nlocal = atom->nlocal; @@ -199,9 +213,9 @@ void DumpCFG::write_header(bigint n) MPI_Allreduce(&vone,&vave,1,MPI_DOUBLE,MPI_SUM,world); if (atom->natoms) vave /= atom->natoms; if (vave > 0.0) scale = 1.44 / pow(vave,1.0/3.0); - else scale = 1.0; - } else scale = 1.0; - + } else if (unwrapflag == 1) scale = UNWRAPEXPAND; + else scale = 1.0; + if (me == 0 || multiproc) { char str[64]; sprintf(str,"Number of particles = %s\n",BIGINT_FORMAT); @@ -261,6 +275,8 @@ void DumpCFG::write_data(int n, double *mybuf) // write data lines in rbuf to file after transfer is done + double unwrap_coord; + if (nlines == nchosen) { for (itype = 1; itype <= ntypes; itype++) { for (i = 0; i < nchosen; i++) @@ -271,11 +287,30 @@ void DumpCFG::write_data(int n, double *mybuf) fprintf(fp,"%s\n",typenames[itype]); for (; i < nchosen; i++) { if (rbuf[i][1] == itype) { - for (j = 2; j < size_one; j++) { - if (vtype[j] == INT) - fprintf(fp,vformat[j],static_cast (rbuf[i][j])); - else fprintf(fp,vformat[j],rbuf[i][j]); - } + if (unwrapflag == 0) + for (j = 2; j < size_one; j++) { + if (vtype[j] == INT) + fprintf(fp,vformat[j],static_cast (rbuf[i][j])); + else fprintf(fp,vformat[j],rbuf[i][j]); + } + else + + // Unwrapped scaled coordinates are shifted to + // center of expanded box, to prevent + // rewrapping by AtomEye. Dividing by + // expansion factor restores correct + // interatomic distances. + + for (j = 2; j < 5; j++) { + unwrap_coord = (rbuf[i][j] - 0.5)/UNWRAPEXPAND + 0.5; + fprintf(fp,vformat[j],unwrap_coord); + } + for (j = 5; j < size_one; j++) { + if (vtype[j] == INT) + fprintf(fp,vformat[j],static_cast (rbuf[i][j])); + else fprintf(fp,vformat[j],rbuf[i][j]); + } + fprintf(fp,"\n"); } } diff --git a/src/dump_cfg.h b/src/dump_cfg.h index 0467983f76..3ddfc8fa44 100755 --- a/src/dump_cfg.h +++ b/src/dump_cfg.h @@ -36,6 +36,7 @@ class DumpCFG : public DumpCustom { int nchosen; // # of lines to be written on a writing proc int nlines; // # of lines transferred from buf to rbuf double **rbuf; // buf of data lines for data lines rearrangement + int unwrapflag; // 1 if unwrapped coordinates are requested void init_style(); void write_header(bigint); diff --git a/src/dump_custom.cpp b/src/dump_custom.cpp index 823c5c0dbf..b0dcc982a2 100644 --- a/src/dump_custom.cpp +++ b/src/dump_custom.cpp @@ -35,7 +35,9 @@ using namespace LAMMPS_NS; // same list as in compute_property.cpp, also customize that command enum{ID,MOL,TYPE,MASS, - X,Y,Z,XS,YS,ZS,XSTRI,YSTRI,ZSTRI,XU,YU,ZU,XUTRI,YUTRI,ZUTRI,IX,IY,IZ, + 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,OMEGAX,OMEGAY,OMEGAZ,ANGMOMX,ANGMOMY,ANGMOMZ, TQX,TQY,TQZ,SPIN,ERADIUS,ERVEL,ERFORCE, @@ -563,6 +565,70 @@ int DumpCustom::count() ptr = dchoose; nstride = 1; + } else if (thresh_array[ithresh] == XSU) { + double **x = atom->x; + int *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] & 1023) - 512; + ptr = dchoose; + nstride = 1; + + } else if (thresh_array[ithresh] == YSU) { + double **x = atom->x; + int *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] >> 10 & 1023) - 512; + ptr = dchoose; + nstride = 1; + + } else if (thresh_array[ithresh] == ZSU) { + double **x = atom->x; + int *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] >> 20) - 512; + ptr = dchoose; + nstride = 1; + + } else if (thresh_array[ithresh] == XSUTRI) { + double **x = atom->x; + int *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] & 1023) - 512; + ptr = dchoose; + nstride = 1; + } else if (thresh_array[ithresh] == YSUTRI) { + double **x = atom->x; + int *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] >> 10 & 1023) - 512; + ptr = dchoose; + nstride = 1; + } else if (thresh_array[ithresh] == ZSUTRI) { + double **x = atom->x; + int *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] >> 20) - 512; + ptr = dchoose; + nstride = 1; + } else if (thresh_array[ithresh] == IX) { int *image = atom->image; for (i = 0; i < nlocal; i++) @@ -879,6 +945,18 @@ void DumpCustom::parse_fields(int narg, char **arg) 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; @@ -1254,6 +1332,19 @@ int DumpCustom::modify_param(int narg, char **arg) 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; @@ -1821,6 +1912,120 @@ void DumpCustom::pack_zu_triclinic(int n) /* ---------------------------------------------------------------------- */ +void DumpCustom::pack_xsu(int n) +{ + double **x = atom->x; + int *image = atom->image; + int nlocal = atom->nlocal; + + double boxxlo = domain->boxlo[0]; + double invxprd = 1.0/domain->xprd; + + for (int i = 0; i < nlocal; i++) + if (choose[i]) { + buf[n] = (x[i][0] - boxxlo) * invxprd + (image[i] & 1023) - 512; + n += size_one; + } +} + +/* ---------------------------------------------------------------------- */ + +void DumpCustom::pack_ysu(int n) +{ + double **x = atom->x; + int *image = atom->image; + int nlocal = atom->nlocal; + + double boxylo = domain->boxlo[1]; + double invyprd = 1.0/domain->yprd; + + for (int i = 0; i < nlocal; i++) + if (choose[i]) { + buf[n] = (x[i][1] - boxylo) * invyprd + (image[i] >> 10 & 1023) - 512; + n += size_one; + } +} + +/* ---------------------------------------------------------------------- */ + +void DumpCustom::pack_zsu(int n) +{ + double **x = atom->x; + int *image = atom->image; + int nlocal = atom->nlocal; + + double boxzlo = domain->boxlo[2]; + double invzprd = 1.0/domain->zprd; + + for (int i = 0; i < nlocal; i++) + if (choose[i]) { + buf[n] = (x[i][2] - boxzlo) * invzprd + (image[i] >> 20) - 512; + n += size_one; + } +} + +/* ---------------------------------------------------------------------- */ + +void DumpCustom::pack_xsu_triclinic(int n) +{ + double **x = atom->x; + int *image = atom->image; + int nlocal = atom->nlocal; + + double *boxlo = domain->boxlo; + double *h_inv = domain->h_inv; + + for (int i = 0; i < nlocal; i++) + if (choose[i]) { + buf[n] = 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] & 1023) - 512; + n += size_one; + } +} + +/* ---------------------------------------------------------------------- */ + +void DumpCustom::pack_ysu_triclinic(int n) +{ + double **x = atom->x; + int *image = atom->image; + int nlocal = atom->nlocal; + + double *boxlo = domain->boxlo; + double *h_inv = domain->h_inv; + + for (int i = 0; i < nlocal; i++) + if (choose[i]) { + buf[n] = h_inv[1]*(x[i][1]-boxlo[1]) + + h_inv[3]*(x[i][2]-boxlo[2]) + + (image[i] >> 10 & 1023) - 512; + n += size_one; + } +} + +/* ---------------------------------------------------------------------- */ + +void DumpCustom::pack_zsu_triclinic(int n) +{ + double **x = atom->x; + int *image = atom->image; + int nlocal = atom->nlocal; + + double *boxlo = domain->boxlo; + double *h_inv = domain->h_inv; + + for (int i = 0; i < nlocal; i++) + if (choose[i]) { + buf[n] = h_inv[2]*(x[i][2]-boxlo[2]) + + (image[i] >> 20) - 512; + n += size_one; + } +} + +/* ---------------------------------------------------------------------- */ + void DumpCustom::pack_ix(int n) { int *image = atom->image;