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