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 #include "compute_displace_atom.h"
16
17 #include "atom.h"
18 #include "domain.h"
19 #include "error.h"
20 #include "fix_store.h"
21 #include "group.h"
22 #include "input.h"
23 #include "memory.h"
24 #include "modify.h"
25 #include "update.h"
26 #include "variable.h"
27
28 #include <cmath>
29 #include <cstring>
30
31 using namespace LAMMPS_NS;
32
33 /* ---------------------------------------------------------------------- */
34
ComputeDisplaceAtom(LAMMPS * lmp,int narg,char ** arg)35 ComputeDisplaceAtom::ComputeDisplaceAtom(LAMMPS *lmp, int narg, char **arg) :
36 Compute(lmp, narg, arg),
37 displace(nullptr), id_fix(nullptr)
38 {
39 if (narg < 3) error->all(FLERR,"Illegal compute displace/atom command");
40
41 peratom_flag = 1;
42 size_peratom_cols = 4;
43 create_attribute = 1;
44
45 // optional args
46
47 refreshflag = 0;
48 rvar = nullptr;
49
50 int iarg = 3;
51 while (iarg < narg) {
52 if (strcmp(arg[iarg],"refresh") == 0) {
53 if (iarg+2 > narg)
54 error->all(FLERR,"Illegal compute displace/atom command");
55 refreshflag = 1;
56 delete [] rvar;
57 rvar = utils::strdup(arg[iarg+1]);
58 iarg += 2;
59 } else error->all(FLERR,"Illegal compute displace/atom command");
60 }
61
62 // error check
63
64 if (refreshflag) {
65 ivar = input->variable->find(rvar);
66 if (ivar < 0)
67 error->all(FLERR,"Variable name for compute displace/atom does not exist");
68 if (input->variable->atomstyle(ivar) == 0)
69 error->all(FLERR,"Compute displace/atom variable "
70 "is not atom-style variable");
71 }
72
73 // create a new fix STORE style
74 // id = compute-ID + COMPUTE_STORE, fix group = compute group
75
76 id_fix = utils::strdup(std::string(id) + "_COMPUTE_STORE");
77 fix = (FixStore *) modify->add_fix(fmt::format("{} {} STORE peratom 1 3",
78 id_fix, group->names[igroup]));
79
80 // calculate xu,yu,zu for fix store array
81 // skip if reset from restart file
82
83 if (fix->restart_reset) fix->restart_reset = 0;
84 else {
85 double **xoriginal = fix->astore;
86
87 double **x = atom->x;
88 int *mask = atom->mask;
89 imageint *image = atom->image;
90 int nlocal = atom->nlocal;
91
92 for (int i = 0; i < nlocal; i++)
93 if (mask[i] & groupbit) domain->unmap(x[i],image[i],xoriginal[i]);
94 else xoriginal[i][0] = xoriginal[i][1] = xoriginal[i][2] = 0.0;
95 }
96
97 // per-atom displacement array
98
99 nmax = nvmax = 0;
100 varatom = nullptr;
101 }
102
103 /* ---------------------------------------------------------------------- */
104
~ComputeDisplaceAtom()105 ComputeDisplaceAtom::~ComputeDisplaceAtom()
106 {
107 // check nfix in case all fixes have already been deleted
108
109 if (modify->nfix) modify->delete_fix(id_fix);
110
111 delete [] id_fix;
112 memory->destroy(displace);
113 delete [] rvar;
114 memory->destroy(varatom);
115 }
116
117 /* ---------------------------------------------------------------------- */
118
init()119 void ComputeDisplaceAtom::init()
120 {
121 // set fix which stores original atom coords
122
123 int ifix = modify->find_fix(id_fix);
124 if (ifix < 0) error->all(FLERR,"Could not find compute displace/atom fix ID");
125 fix = (FixStore *) modify->fix[ifix];
126
127 if (refreshflag) {
128 ivar = input->variable->find(rvar);
129 if (ivar < 0)
130 error->all(FLERR,"Variable name for compute displace/atom does not exist");
131 }
132 }
133
134 /* ---------------------------------------------------------------------- */
135
compute_peratom()136 void ComputeDisplaceAtom::compute_peratom()
137 {
138 invoked_peratom = update->ntimestep;
139
140 // grow local displacement array if necessary
141
142 if (atom->nmax > nmax) {
143 memory->destroy(displace);
144 nmax = atom->nmax;
145 memory->create(displace,nmax,4,"displace/atom:displace");
146 array_atom = displace;
147 }
148
149 // dx,dy,dz = displacement of atom from original position
150 // original unwrapped position is stored by fix
151 // for triclinic, need to unwrap current atom coord via h matrix
152
153 double **xoriginal = fix->astore;
154
155 double **x = atom->x;
156 int *mask = atom->mask;
157 imageint *image = atom->image;
158 int nlocal = atom->nlocal;
159
160 double *h = domain->h;
161 double xprd = domain->xprd;
162 double yprd = domain->yprd;
163 double zprd = domain->zprd;
164
165 int xbox,ybox,zbox;
166 double dx,dy,dz;
167
168 if (domain->triclinic == 0) {
169 for (int i = 0; i < nlocal; i++)
170 if (mask[i] & groupbit) {
171 xbox = (image[i] & IMGMASK) - IMGMAX;
172 ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
173 zbox = (image[i] >> IMG2BITS) - IMGMAX;
174 dx = x[i][0] + xbox*xprd - xoriginal[i][0];
175 dy = x[i][1] + ybox*yprd - xoriginal[i][1];
176 dz = x[i][2] + zbox*zprd - xoriginal[i][2];
177 displace[i][0] = dx;
178 displace[i][1] = dy;
179 displace[i][2] = dz;
180 displace[i][3] = sqrt(dx*dx + dy*dy + dz*dz);
181 } else displace[i][0] = displace[i][1] =
182 displace[i][2] = displace[i][3] = 0.0;
183
184 } else {
185 for (int i = 0; i < nlocal; i++)
186 if (mask[i] & groupbit) {
187 xbox = (image[i] & IMGMASK) - IMGMAX;
188 ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
189 zbox = (image[i] >> IMG2BITS) - IMGMAX;
190 dx = x[i][0] + h[0]*xbox + h[5]*ybox + h[4]*zbox - xoriginal[i][0];
191 dy = x[i][1] + h[1]*ybox + h[3]*zbox - xoriginal[i][1];
192 dz = x[i][2] + h[2]*zbox - xoriginal[i][2];
193 displace[i][0] = dx;
194 displace[i][1] = dy;
195 displace[i][2] = dz;
196 displace[i][3] = sqrt(dx*dx + dy*dy + dz*dz);
197 } else displace[i][0] = displace[i][1] =
198 displace[i][2] = displace[i][3] = 0.0;
199 }
200 }
201
202 /* ----------------------------------------------------------------------
203 initialize one atom's storage values, called when atom is created
204 ------------------------------------------------------------------------- */
205
set_arrays(int i)206 void ComputeDisplaceAtom::set_arrays(int i)
207 {
208 double **xoriginal = fix->astore;
209 double **x = atom->x;
210 xoriginal[i][0] = x[i][0];
211 xoriginal[i][1] = x[i][1];
212 xoriginal[i][2] = x[i][2];
213 }
214
215 /* ----------------------------------------------------------------------
216 reset per-atom storage values, based on atom-style variable evaluation
217 called by dump when dump_modify refresh is set
218 ------------------------------------------------------------------------- */
219
refresh()220 void ComputeDisplaceAtom::refresh()
221 {
222 if (atom->nmax > nvmax) {
223 nvmax = atom->nmax;
224 memory->destroy(varatom);
225 memory->create(varatom,nvmax,"displace/atom:varatom");
226 }
227
228 input->variable->compute_atom(ivar,igroup,varatom,1,0);
229
230 double **xoriginal = fix->astore;
231 double **x = atom->x;
232 imageint *image = atom->image;
233 int nlocal = atom->nlocal;
234
235 for (int i = 0; i < nlocal; i++)
236 if (varatom[i]) domain->unmap(x[i],image[i],xoriginal[i]);
237 }
238
239 /* ----------------------------------------------------------------------
240 memory usage of local atom-based array
241 ------------------------------------------------------------------------- */
242
memory_usage()243 double ComputeDisplaceAtom::memory_usage()
244 {
245 double bytes = (double)nmax*4 * sizeof(double);
246 bytes += (double)nvmax * sizeof(double);
247 return bytes;
248 }
249