1 // clang-format off
2 /* ----------------------------------------------------------------------
3    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
4    https://www.lammps.org/, Sandia National Laboratories
5    Steve Plimpton, sjplimp@sandia.gov
6 
7    Copyright (2003) Sandia Corporation.  Under the terms of Contract
8    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
9    certain rights in this software.  This software is distributed under
10    the GNU General Public License.
11 
12    See the README file in the top-level LAMMPS directory.
13 ------------------------------------------------------------------------- */
14 
15 /* ----------------------------------------------------------------------
16    Contributing author: Naveen Michaud-Agrawal (Johns Hopkins U)
17 ------------------------------------------------------------------------- */
18 
19 #include "fix_recenter.h"
20 
21 #include "atom.h"
22 #include "comm.h"
23 #include "domain.h"
24 #include "error.h"
25 #include "group.h"
26 #include "lattice.h"
27 #include "modify.h"
28 #include "respa.h"
29 #include "update.h"
30 
31 #include <cmath>
32 #include <cstring>
33 
34 using namespace LAMMPS_NS;
35 using namespace FixConst;
36 
37 enum{BOX,LATTICE,FRACTION};
38 
39 /* ---------------------------------------------------------------------- */
40 
FixRecenter(LAMMPS * lmp,int narg,char ** arg)41 FixRecenter::FixRecenter(LAMMPS *lmp, int narg, char **arg) :
42   Fix(lmp, narg, arg)
43 {
44   if (narg < 6) error->all(FLERR,"Illegal fix recenter command");
45 
46   xcom = ycom = zcom = 0.0;
47   xflag = yflag = zflag = 1;
48   xinitflag = yinitflag = zinitflag = 0;
49   shift[0] = shift[1] = shift[2] = 0.0;
50   distance = 0.0;
51   scalar_flag = 1;
52   vector_flag = 1;
53   size_vector = 3;
54   extscalar = 1;
55   extvector = 1;
56   global_freq = 1;
57   dynamic_group_allow = 1;
58 
59   if (strcmp(arg[3],"NULL") == 0) xflag = 0;
60   else if (strcmp(arg[3],"INIT") == 0) xinitflag = 1;
61   else xcom = utils::numeric(FLERR,arg[3],false,lmp);
62   if (strcmp(arg[4],"NULL") == 0) yflag = 0;
63   else if (strcmp(arg[4],"INIT") == 0) yinitflag = 1;
64   else ycom = utils::numeric(FLERR,arg[4],false,lmp);
65   if (strcmp(arg[5],"NULL") == 0) zflag = 0;
66   else if (strcmp(arg[5],"INIT") == 0) zinitflag = 1;
67   else zcom = utils::numeric(FLERR,arg[5],false,lmp);
68 
69   // optional args
70 
71   group2bit = groupbit;
72   scaleflag = LATTICE;
73 
74   int iarg = 6;
75   while (iarg < narg) {
76     if (strcmp(arg[iarg],"shift") == 0) {
77       int igroup2 = group->find(arg[iarg+1]);
78       if (igroup2 < 0) error->all(FLERR,"Could not find fix recenter group ID");
79       group2bit = group->bitmask[igroup2];
80       iarg += 2;
81     } else if (strcmp(arg[iarg],"units") == 0) {
82       if (strcmp(arg[iarg+1],"box") == 0) scaleflag = BOX;
83       else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = LATTICE;
84       else if (strcmp(arg[iarg+1],"fraction") == 0) scaleflag = FRACTION;
85       else error->all(FLERR,"Illegal fix recenter command");
86       iarg += 2;
87     } else error->all(FLERR,"Illegal fix recenter command");
88   }
89 
90   // scale xcom,ycom,zcom
91 
92   double xscale,yscale,zscale;
93   if (scaleflag == LATTICE) {
94     xscale = domain->lattice->xlattice;
95     yscale = domain->lattice->ylattice;
96     zscale = domain->lattice->zlattice;
97   }
98   else xscale = yscale = zscale = 1.0;
99 
100   xcom *= xscale;
101   ycom *= yscale;
102   zcom *= zscale;
103 
104   // cannot have 0 atoms in group
105 
106   if (group->count(igroup) == 0)
107     error->all(FLERR,"Fix recenter group has no atoms");
108 }
109 
110 /* ---------------------------------------------------------------------- */
111 
setmask()112 int FixRecenter::setmask()
113 {
114   int mask = 0;
115   mask |= INITIAL_INTEGRATE;
116   mask |= INITIAL_INTEGRATE_RESPA;
117   return mask;
118 }
119 
120 /* ---------------------------------------------------------------------- */
121 
init()122 void FixRecenter::init()
123 {
124   // warn if any integrate fix comes after this one
125 
126   int after = 0;
127   int flag = 0;
128   for (int i = 0; i < modify->nfix; i++) {
129     if (strcmp(id,modify->fix[i]->id) == 0) after = 1;
130     else if ((modify->fmask[i] & INITIAL_INTEGRATE) && after) flag = 1;
131   }
132   if (flag && comm->me == 0)
133     error->warning(FLERR,"Fix recenter should come after all other "
134                    "integration fixes");
135 
136   masstotal = group->mass(igroup);
137 
138   // if any components of requested COM were INIT, store initial COM
139 
140   if (xinitflag || yinitflag || zinitflag) {
141     double xcm[3];
142     group->xcm(igroup,masstotal,xcm);
143     xinit = xcm[0];
144     yinit = xcm[1];
145     zinit = xcm[2];
146   }
147 
148   if (utils::strmatch(update->integrate_style,"^respa"))
149     nlevels_respa = ((Respa *) update->integrate)->nlevels;
150 }
151 
152 /* ---------------------------------------------------------------------- */
153 
initial_integrate(int)154 void FixRecenter::initial_integrate(int /*vflag*/)
155 {
156   // target COM
157   // bounding box around domain works for both orthogonal and triclinic
158 
159   double xtarget,ytarget,ztarget;
160   double *bboxlo,*bboxhi;
161 
162   if (scaleflag == FRACTION) {
163     if (domain->triclinic == 0) {
164       bboxlo = domain->boxlo;
165       bboxhi = domain->boxhi;
166     } else {
167       bboxlo = domain->boxlo_bound;
168       bboxhi = domain->boxhi_bound;
169     }
170   }
171 
172   if (xinitflag) xtarget = xinit;
173   else if (scaleflag == FRACTION)
174     xtarget = bboxlo[0] + xcom*(bboxhi[0] - bboxlo[0]);
175   else xtarget = xcom;
176 
177   if (yinitflag) ytarget = yinit;
178   else if (scaleflag == FRACTION)
179     ytarget = bboxlo[1] + ycom*(bboxhi[1] - bboxlo[1]);
180   else ytarget = ycom;
181 
182   if (zinitflag) ztarget = zinit;
183   else if (scaleflag == FRACTION)
184     ztarget = bboxlo[2] + zcom*(bboxhi[2] - bboxlo[2]);
185   else ztarget = zcom;
186 
187   // current COM
188 
189   double xcm[3];
190   if (group->dynamic[igroup])
191     masstotal = group->mass(igroup);
192 
193   group->xcm(igroup,masstotal,xcm);
194 
195   // shift coords by difference between actual COM and requested COM
196 
197   double **x = atom->x;
198   int *mask = atom->mask;
199   int nlocal = atom->nlocal;
200 
201   shift[0] = xflag ? (xtarget - xcm[0]) : 0.0;
202   shift[1] = yflag ? (ytarget - xcm[1]) : 0.0;
203   shift[2] = zflag ? (ztarget - xcm[2]) : 0.0;
204   distance = sqrt(shift[0]*shift[0] + shift[1]*shift[1] + shift[2]*shift[2]);
205 
206   for (int i = 0; i < nlocal; i++)
207     if (mask[i] & group2bit) {
208       x[i][0] += shift[0];
209       x[i][1] += shift[1];
210       x[i][2] += shift[2];
211     }
212 }
213 
214 /* ---------------------------------------------------------------------- */
215 
initial_integrate_respa(int vflag,int ilevel,int)216 void FixRecenter::initial_integrate_respa(int vflag, int ilevel, int /*iloop*/)
217 {
218   // outermost level - operate recenter
219   // all other levels - nothing
220 
221   if (ilevel == nlevels_respa-1) initial_integrate(vflag);
222 }
223 
224 /* ---------------------------------------------------------------------- */
225 
compute_scalar()226 double FixRecenter::compute_scalar()
227 {
228   return distance;
229 }
230 
231 /* ---------------------------------------------------------------------- */
232 
compute_vector(int n)233 double FixRecenter::compute_vector(int n)
234 {
235   return shift[n];
236 }
237