diff --git a/src/fix.h b/src/fix.h index 551358dc49..9d1fc93cfd 100644 --- a/src/fix.h +++ b/src/fix.h @@ -93,7 +93,7 @@ class Fix : protected Pointers { virtual double thermo(int) {return 0.0;} virtual int dof(int) {return 0;} - virtual void dilate(int, double, double, double, double) {} + virtual void deform(int) {} virtual int modify_param(int, char **) {return 0;} diff --git a/src/fix_nph.h b/src/fix_nph.h index 5b47bef6bb..60b79846fb 100644 --- a/src/fix_nph.h +++ b/src/fix_nph.h @@ -40,7 +40,7 @@ class FixNPH : public Fix { double boltz,nktv2p; double vol0,nkt; - int press_couple,dilate_partial; + int press_couple,allremap; int p_flag[3]; // 1 if control P on this dim, 0 if not double p_start[3],p_stop[3]; double p_freq[3],p_target[3]; @@ -60,7 +60,7 @@ class FixNPH : public Fix { int tflag,pflag; void couple(); - void box_dilate(int); + void remap(int); }; } diff --git a/src/fix_npt.cpp b/src/fix_npt.cpp index 048b2aa425..6b2c6f2b99 100644 --- a/src/fix_npt.cpp +++ b/src/fix_npt.cpp @@ -22,7 +22,6 @@ #include "atom.h" #include "force.h" #include "comm.h" -#include "output.h" #include "modify.h" #include "compute.h" #include "kspace.h" @@ -117,7 +116,7 @@ FixNPT::FixNPT(LAMMPS *lmp, int narg, char **arg) : // process extra keywords drag = 0.0; - dilate_partial = 0; + allremap = 1; int iarg; if (press_couple == 0) iarg = 10; @@ -130,8 +129,8 @@ FixNPT::FixNPT(LAMMPS *lmp, int narg, char **arg) : iarg += 2; } else if (strcmp(arg[iarg],"dilate") == 0) { if (iarg+2 > narg) error->all("Illegal fix npt command"); - if (strcmp(arg[iarg+1],"all") == 0) dilate_partial = 0; - else if (strcmp(arg[iarg+1],"partial") == 0) dilate_partial = 1; + if (strcmp(arg[iarg+1],"all") == 0) allremap = 1; + else if (strcmp(arg[iarg+1],"partial") == 0) allremap = 0; else error->all("Illegal fix npt command"); iarg += 2; } else error->all("Illegal fix npt command"); @@ -284,7 +283,7 @@ void FixNPT::init() step_respa = ((Respa *) update->integrate)->step; } - // detect if any rigid fixes exist so rigid bodies move when box is dilated + // detect if any rigid fixes exist so rigid bodies move when box is remapped // rfix[] = indices to each fix rigid delete [] rfix; @@ -381,9 +380,9 @@ void FixNPT::initial_integrate() } } - // rescale simulation box and all owned atoms by 1/2 step + // remap simulation box and all owned atoms by 1/2 step - box_dilate(0); + remap(0); // x update by full step only for atoms in NPT group @@ -395,10 +394,10 @@ void FixNPT::initial_integrate() } } - // rescale simulation box and all owned atoms by 1/2 step + // remap simulation box and all owned atoms by 1/2 step // redo KSpace coeffs since volume has changed - box_dilate(0); + remap(0); if (kspace_flag) force->kspace->setup(); } @@ -468,11 +467,11 @@ void FixNPT::final_integrate() void FixNPT::initial_integrate_respa(int ilevel, int flag) { // if flag = 1, then is 2nd call at outermost level from rRESPA - // perform 2nd half of box dilate on own + ghost atoms and return + // perform 2nd half of box remap on own + ghost atoms and return // redo KSpace coeffs since volume has changed if (flag == 1) { - box_dilate(1); + remap(1); if (kspace_flag) force->kspace->setup(); return; } @@ -494,7 +493,7 @@ void FixNPT::initial_integrate_respa(int ilevel, int flag) int *mask = atom->mask; int nlocal = atom->nlocal; - // outermost level - update eta_dot and omega_dot, apply to v, dilate box + // outermost level - update eta_dot and omega_dot, apply to v, remap box // all other levels - NVE update of v // x,v updates only performed for atoms in NPT group @@ -540,9 +539,9 @@ void FixNPT::initial_integrate_respa(int ilevel, int flag) } } - // rescale simulation box and all owned atoms by 1/2 step + // remap simulation box and all owned atoms by 1/2 step - box_dilate(0); + remap(0); } else { @@ -638,89 +637,62 @@ void FixNPT::couple() } /* ---------------------------------------------------------------------- - dilate the box around center of box + change box size + remap owned or owned+ghost atoms depending on flag + remap all atoms or fix group atoms depending on allremap flag + if rigid bodies exist, scale rigid body centers-of-mass ------------------------------------------------------------------------- */ -void FixNPT::box_dilate(int flag) +void FixNPT::remap(int flag) { int i,n; - - // ctr = geometric center of box in a dimension - // scale owned or owned+ghost atoms depending on flag - // re-define simulation box via xprd/yprd/zprd - // scale atom coords for all atoms or only for fix group atoms - // if fix rigid exists, scale rigid body centers-of-mass - // don't do anything if non-periodic or press style is constant volume + double oldlo,oldhi,ctr; double **x = atom->x; int *mask = atom->mask; if (flag) n = atom->nlocal + atom->nghost; else n = atom->nlocal; - double oldlo,oldhi,ctr; + // convert pertinent atoms and rigid bodies to lamda coords - if (domain->xperiodic && p_flag[0]) { - oldlo = domain->boxlo[0]; - oldhi = domain->boxhi[0]; - ctr = 0.5 * (oldlo + oldhi); - domain->boxlo[0] = (oldlo-ctr)*dilation[0] + ctr; - domain->boxhi[0] = (oldhi-ctr)*dilation[0] + ctr; - domain->prd[0] = domain->xprd = domain->boxhi[0] - domain->boxlo[0]; - if (dilate_partial) { - for (i = 0; i < n; i++) - if (mask[i] & groupbit) - x[i][0] = ctr + (x[i][0]-ctr)*dilation[0]; - } else { - for (i = 0; i < n; i++) - x[i][0] = ctr + (x[i][0]-ctr)*dilation[0]; - } - if (nrigid) - for (i = 0; i < nrigid; i++) - modify->fix[rfix[i]]-> - dilate(0,oldlo,oldhi,domain->boxlo[0],domain->boxhi[0]); + if (allremap) domain->x2lamda(n); + else { + for (i = 0; i < n; i++) + if (mask[i] & groupbit) + domain->x2lamda(x[i],x[i]); } - if (domain->yperiodic && p_flag[1]) { - oldlo = domain->boxlo[1]; - oldhi = domain->boxhi[1]; - ctr = 0.5 * (oldlo + oldhi); - domain->boxlo[1] = (oldlo-ctr)*dilation[1] + ctr; - domain->boxhi[1] = (oldhi-ctr)*dilation[1] + ctr; - domain->prd[1] = domain->yprd = domain->boxhi[1] - domain->boxlo[1]; - if (dilate_partial) { - for (i = 0; i < n; i++) - if (mask[i] & groupbit) - x[i][1] = ctr + (x[i][1]-ctr)*dilation[1]; - } else { - for (i = 0; i < n; i++) - x[i][1] = ctr + (x[i][1]-ctr)*dilation[1]; + if (nrigid) + for (i = 0; i < nrigid; i++) + modify->fix[rfix[i]]->deform(0); + + // reset global and local box to new size/shape + + for (i = 0; i < 3; i++) { + if (p_flag[i]) { + oldlo = domain->boxlo[i]; + oldhi = domain->boxhi[i]; + ctr = 0.5 * (oldlo + oldhi); + domain->boxlo[i] = (oldlo-ctr)*dilation[i] + ctr; + domain->boxhi[i] = (oldhi-ctr)*dilation[i] + ctr; } - if (nrigid) - for (i = 0; i < nrigid; i++) - modify->fix[rfix[i]]-> - dilate(1,oldlo,oldhi,domain->boxlo[1],domain->boxhi[1]); } - if (domain->zperiodic && p_flag[2]) { - oldlo = domain->boxlo[2]; - oldhi = domain->boxhi[2]; - ctr = 0.5 * (oldlo + oldhi); - domain->boxlo[2] = (oldlo-ctr)*dilation[2] + ctr; - domain->boxhi[2] = (oldhi-ctr)*dilation[2] + ctr; - domain->prd[2] = domain->zprd = domain->boxhi[2] - domain->boxlo[2]; - if (dilate_partial) { - for (i = 0; i < n; i++) - if (mask[i] & groupbit) - x[i][2] = ctr + (x[i][2]-ctr)*dilation[2]; - } else { - for (i = 0; i < n; i++) - x[i][2] = ctr + (x[i][2]-ctr)*dilation[2]; - } - if (nrigid) - for (i = 0; i < nrigid; i++) - modify->fix[rfix[i]]-> - dilate(2,oldlo,oldhi,domain->boxlo[2],domain->boxhi[2]); + domain->set_global_box(); + domain->set_local_box(); + + // convert pertinent atoms and rigid bodies back to box coords + + if (allremap) domain->lamda2x(n); + else { + for (i = 0; i < n; i++) + if (mask[i] & groupbit) + domain->lamda2x(x[i],x[i]); } + + if (nrigid) + for (i = 0; i < nrigid; i++) + modify->fix[rfix[i]]->deform(1); } /* ---------------------------------------------------------------------- diff --git a/src/fix_npt.h b/src/fix_npt.h index 15f36d55dc..0a32120ba4 100644 --- a/src/fix_npt.h +++ b/src/fix_npt.h @@ -45,7 +45,7 @@ class FixNPT : public Fix { double t_freq; double f_eta,eta_dot,eta; - int press_couple,dilate_partial; + int press_couple,allremap; int p_flag[3]; // 1 if control P on this dim, 0 if not double p_start[3],p_stop[3]; double p_freq[3],p_target[3]; @@ -65,7 +65,7 @@ class FixNPT : public Fix { int tflag,pflag; void couple(); - void box_dilate(int); + void remap(int); }; } diff --git a/src/fix_rigid.cpp b/src/fix_rigid.cpp index 8a5c4d749c..7b767bdfcb 100644 --- a/src/fix_rigid.cpp +++ b/src/fix_rigid.cpp @@ -847,18 +847,19 @@ int FixRigid::dof(int igroup) } /* ---------------------------------------------------------------------- - adjust xcm of each rigid body due to box dilation in idim - called by various fixes that change box + adjust xcm of each rigid body due to box deformation + called by various fixes that change box size/shape + flag = 0/1 means map from box to lamda coords or vice versa ------------------------------------------------------------------------- */ -void FixRigid::dilate(int idim, - double oldlo, double oldhi, double newlo, double newhi) +void FixRigid::deform(int flag) { - double ratio; - for (int ibody = 0; ibody < nbody; ibody++) { - ratio = (xcm[ibody][idim] - oldlo) / (oldhi - oldlo); - xcm[ibody][idim] = newlo + ratio*(newhi - newlo); - } + if (flag == 0) + for (int ibody = 0; ibody < nbody; ibody++) + domain->x2lamda(xcm[ibody],xcm[ibody]); + else + for (int ibody = 0; ibody < nbody; ibody++) + domain->lamda2x(xcm[ibody],xcm[ibody]); } /* ----------------------------------------------------------------------