diff --git a/src/USER-COLVARS/fix_colvars.cpp b/src/USER-COLVARS/fix_colvars.cpp index 2f8ca857cb..9adde33f6d 100644 --- a/src/USER-COLVARS/fix_colvars.cpp +++ b/src/USER-COLVARS/fix_colvars.cpp @@ -300,7 +300,7 @@ FixColvars::FixColvars(LAMMPS *lmp, int narg, char **arg) : if (atom->rmass_flag) error->all(FLERR,"Cannot use fix colvars for atoms with rmass attribute"); - if (instances) + if (instances > 0) error->all(FLERR,"Only one fix colvars can be active at a time"); ++instances; @@ -313,6 +313,7 @@ FixColvars::FixColvars(LAMMPS *lmp, int narg, char **arg) : conf_file = strdup(arg[3]); rng_seed = 1966; + unwrap_flag = 0; inp_name = NULL; out_name = NULL; @@ -327,6 +328,14 @@ FixColvars::FixColvars(LAMMPS *lmp, int narg, char **arg) : out_name = strdup(arg[argsdone+1]); } else if (0 == strcmp(arg[argsdone], "seed")) { rng_seed = atoi(arg[argsdone+1]); + } else if (0 == strcmp(arg[argsdone], "unwrap")) { + if (0 == strcmp(arg[argsdone+1], "yes")) { + unwrap_flag = 1; + } else if (0 == strcmp(arg[argsdone+1], "no")) { + unwrap_flag = 0; + } else { + error->all(FLERR,"Incorrect fix colvars unwrap flag"); + } } else if (0 == strcmp(arg[argsdone], "tstat")) { tmp_name = strdup(arg[argsdone+1]); } else { @@ -538,6 +547,14 @@ void FixColvars::init() memory->create(comm_buf,nmax,"colvars:comm_buf"); const double * const * const x = atom->x; + const tagint * const image = atom->image; + + const double xprd = domain->xprd; + const double yprd = domain->yprd; + const double zprd = domain->zprd; + const double xy = domain->xy; + const double xz = domain->xz; + const double yz = domain->yz; if (me == 0) { @@ -549,12 +566,24 @@ void FixColvars::init() for (i=0; imap(taglist[i]); if ((k >= 0) && (k < nlocal)) { + of[i].tag = cd[i].tag = tag[k]; of[i].type = cd[i].type = type[k]; - cd[i].x = x[k][0]; - cd[i].y = x[k][1]; - cd[i].z = x[k][2]; of[i].x = of[i].y = of[i].z = 0.0; + + if (unwrap_flag) { + const int ix = (image[k] & IMGMASK) - IMGMAX; + const int iy = (image[k] >> IMGBITS & IMGMASK) - IMGMAX; + const int iz = (image[k] >> IMG2BITS) - IMGMAX; + + cd[i].x = x[k][0] + ix * xprd + iy * xy + iz * xz; + cd[i].y = x[k][1] + iy * yprd + iz * yz; + cd[i].z = x[k][2] + iz * zprd; + } else { + cd[i].x = x[k][0]; + cd[i].y = x[k][1]; + cd[i].z = x[k][2]; + } } } @@ -588,11 +617,24 @@ void FixColvars::init() for (i=0; imap(taglist[i]); if ((k >= 0) && (k < nlocal)) { + comm_buf[nme].tag = tag[k]; comm_buf[nme].type = type[k]; - comm_buf[nme].x = x[k][0]; - comm_buf[nme].y = x[k][1]; - comm_buf[nme].z = x[k][2]; + + if (unwrap_flag) { + const int ix = (image[k] & IMGMASK) - IMGMAX; + const int iy = (image[k] >> IMGBITS & IMGMASK) - IMGMAX; + const int iz = (image[k] >> IMG2BITS) - IMGMAX; + + comm_buf[nme].x = x[k][0] + ix * xprd + iy * xy + iz * xz; + comm_buf[nme].y = x[k][1] + iy * yprd + iz * yz; + comm_buf[nme].z = x[k][2] + iz * zprd; + } else { + comm_buf[nme].x = x[k][0]; + comm_buf[nme].y = x[k][1]; + comm_buf[nme].z = x[k][2]; + } + ++nme; } } @@ -643,6 +685,14 @@ void FixColvars::post_force(int vflag) const int * const tag = atom->tag; const double * const * const x = atom->x; double * const * const f = atom->f; + const tagint * const image = atom->image; + + const double xprd = domain->xprd; + const double yprd = domain->yprd; + const double zprd = domain->zprd; + const double xy = domain->xy; + const double xz = domain->xz; + const double yz = domain->yz; const int nlocal = atom->nlocal; /* check and potentially grow local communication buffers. */ @@ -671,9 +721,20 @@ void FixColvars::post_force(int vflag) for (i=0; imap(taglist[i]); if ((k >= 0) && (k < nlocal)) { - cd[i].x = x[k][0]; - cd[i].y = x[k][1]; - cd[i].z = x[k][2]; + + if (unwrap_flag) { + const int ix = (image[k] & IMGMASK) - IMGMAX; + const int iy = (image[k] >> IMGBITS & IMGMASK) - IMGMAX; + const int iz = (image[k] >> IMG2BITS) - IMGMAX; + + cd[i].x = x[k][0] + ix * xprd + iy * xy + iz * xz; + cd[i].y = x[k][1] + iy * yprd + iz * yz; + cd[i].z = x[k][2] + iz * zprd; + } else { + cd[i].x = x[k][0]; + cd[i].y = x[k][1]; + cd[i].z = x[k][2]; + } } } @@ -703,9 +764,21 @@ void FixColvars::post_force(int vflag) const int k = atom->map(taglist[i]); if ((k >= 0) && (k < nlocal)) { comm_buf[nme].tag = tag[k]; - comm_buf[nme].x = x[k][0]; - comm_buf[nme].y = x[k][1]; - comm_buf[nme].z = x[k][2]; + + if (unwrap_flag) { + const int ix = (image[k] & IMGMASK) - IMGMAX; + const int iy = (image[k] >> IMGBITS & IMGMASK) - IMGMAX; + const int iz = (image[k] >> IMG2BITS) - IMGMAX; + + comm_buf[nme].x = x[k][0] + ix * xprd + iy * xy + iz * xz; + comm_buf[nme].y = x[k][1] + iy * yprd + iz * yz; + comm_buf[nme].z = x[k][2] + iz * zprd; + } else { + comm_buf[nme].x = x[k][0]; + comm_buf[nme].y = x[k][1]; + comm_buf[nme].z = x[k][2]; + } + ++nme; } } diff --git a/src/USER-COLVARS/fix_colvars.h b/src/USER-COLVARS/fix_colvars.h index de144edc58..6e1a5f8993 100644 --- a/src/USER-COLVARS/fix_colvars.h +++ b/src/USER-COLVARS/fix_colvars.h @@ -80,6 +80,7 @@ class FixColvars : public Fix { int nlevels_respa; // flag to determine respa levels. int store_forces; // flag to determine whether to store total forces + int unwrap_flag; // 1 if atom coords are unwrapped, 0 if not static int instances; // count fix instances, since colvars currently // only supports one instance at a time };