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