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: Pablo Piaggi (EPFL Lausanne)
17 ------------------------------------------------------------------------- */
18
19 #include "compute_entropy_atom.h"
20 #include <cmath>
21 #include <cstring>
22 #include "atom.h"
23 #include "update.h"
24 #include "modify.h"
25 #include "neighbor.h"
26 #include "neigh_list.h"
27 #include "neigh_request.h"
28 #include "force.h"
29 #include "pair.h"
30 #include "comm.h"
31 #include "math_const.h"
32 #include "memory.h"
33 #include "error.h"
34 #include "domain.h"
35
36 using namespace LAMMPS_NS;
37 using namespace MathConst;
38
39 /* ---------------------------------------------------------------------- */
40
41 ComputeEntropyAtom::
ComputeEntropyAtom(LAMMPS * lmp,int narg,char ** arg)42 ComputeEntropyAtom(LAMMPS *lmp, int narg, char **arg) :
43 Compute(lmp, narg, arg),
44 pair_entropy(nullptr), pair_entropy_avg(nullptr)
45 {
46 if (narg < 5 || narg > 10)
47 error->all(FLERR,"Illegal compute entropy/atom command; wrong number"
48 " of arguments");
49
50 // Arguments are: sigma cutoff avg yes/no cutoff2 local yes/no
51 // sigma is the gaussian width
52 // cutoff is the cutoff for the calculation of g(r)
53 // avg is optional and allows averaging the pair entropy over neighbors
54 // the next argument should be yes or no
55 // cutoff2 is the cutoff for the averaging
56 // local is optional and allows using the local density to normalize
57 // the g(r)
58
59 sigma = utils::numeric(FLERR,arg[3],false,lmp);
60 if (sigma <= 0.0) error->all(FLERR,"Illegal compute entropy/atom"
61 " command; sigma must be positive");
62 cutoff = utils::numeric(FLERR,arg[4],false,lmp);
63 if (cutoff <= 0.0) error->all(FLERR,"Illegal compute entropy/atom"
64 " command; cutoff must be positive");
65
66 cutoff2 = 0.;
67 avg_flag = 0;
68 local_flag = 0;
69
70 // optional keywords
71 int iarg = 5;
72 while (iarg < narg) {
73 if (strcmp(arg[iarg],"avg") == 0) {
74 if (iarg+2 > narg)
75 error->all(FLERR,"Illegal compute entropy/atom;"
76 " missing arguments after avg");
77 if (strcmp(arg[iarg+1],"yes") == 0) avg_flag = 1;
78 else if (strcmp(arg[iarg+1],"no") == 0) avg_flag = 0;
79 else error->all(FLERR,"Illegal compute entropy/atom;"
80 " argument after avg should be yes or no");
81 cutoff2 = utils::numeric(FLERR,arg[iarg+2],false,lmp);
82 if (cutoff2 < 0.0) error->all(FLERR,"Illegal compute entropy/atom"
83 " command; negative cutoff2");
84 cutsq2 = cutoff2*cutoff2;
85 iarg += 3;
86 } else if (strcmp(arg[iarg],"local") == 0) {
87 if (strcmp(arg[iarg+1],"yes") == 0) local_flag = 1;
88 else if (strcmp(arg[iarg+1],"no") == 0) local_flag = 0;
89 else error->all(FLERR,"Illegal compute entropy/atom;"
90 " argument after local should be yes or no");
91 iarg += 2;
92 } else error->all(FLERR,"Illegal compute entropy/atom; argument after"
93 " sigma and cutoff should be avg or local");
94 }
95
96
97 cutsq = cutoff*cutoff;
98 nbin = static_cast<int>(cutoff / sigma) + 1;
99 nmax = 0;
100 maxneigh = 0;
101 // Number of bins above and below the central one that will be
102 // considered as affected by the gaussian kernel
103 // 3 seems a good compromise between speed and good mollification
104 deltabin = 3;
105 deltar = sigma;
106 peratom_flag = 1;
107 size_peratom_cols = 0;
108 }
109
110 /* ---------------------------------------------------------------------- */
111
~ComputeEntropyAtom()112 ComputeEntropyAtom::~ComputeEntropyAtom()
113 {
114 memory->destroy(pair_entropy);
115 if (avg_flag) memory->destroy(pair_entropy_avg);
116 }
117
118 /* ---------------------------------------------------------------------- */
119
init()120 void ComputeEntropyAtom::init()
121 {
122 if (force->pair == nullptr)
123 error->all(FLERR,"Compute entropy/atom requires a pair style be"
124 " defined");
125
126 if ((cutoff+cutoff2) > (force->pair->cutforce + neighbor->skin))
127 {
128 error->all(FLERR,"Compute entropy/atom cutoff is longer than the"
129 " pairwise cutoff. Increase the neighbor list skin"
130 " distance.");
131 }
132
133 int count = 0;
134 for (int i = 0; i < modify->ncompute; i++)
135 if (strcmp(modify->compute[i]->style,"entropy/atom") == 0) count++;
136 if (count > 1 && comm->me == 0)
137 error->warning(FLERR,"More than one compute entropy/atom");
138
139 // Request neighbor list
140 int irequest = neighbor->request(this,instance_me);
141 neighbor->requests[irequest]->pair = 0;
142 neighbor->requests[irequest]->compute = 1;
143 neighbor->requests[irequest]->half = 0;
144 neighbor->requests[irequest]->full = 1;
145 if (avg_flag) {
146 // need a full neighbor list with neighbors of the ghost atoms
147 neighbor->requests[irequest]->occasional = 0;
148 neighbor->requests[irequest]->ghost = 1;
149 } else {
150 // need a regular full neighbor list
151 // can build it occasionally
152 neighbor->requests[irequest]->occasional = 1;
153 neighbor->requests[irequest]->ghost = 0;
154 }
155
156 }
157
158 /* ---------------------------------------------------------------------- */
159
init_list(int,NeighList * ptr)160 void ComputeEntropyAtom::init_list(int /*id*/, NeighList *ptr)
161 {
162 list = ptr;
163 }
164
165 /* ---------------------------------------------------------------------- */
166
compute_peratom()167 void ComputeEntropyAtom::compute_peratom()
168 {
169 int i,j,ii,jj,inum,jnum;
170 double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
171 int *ilist,*jlist,*numneigh,**firstneigh;
172 double *rbin = new double[nbin];
173 double *rbinsq = new double[nbin];
174
175 invoked_peratom = update->ntimestep;
176
177 // Initialize distance vectors
178 for (int i = 0; i < nbin; i++) {
179 rbin[i] = i*deltar;
180 rbinsq[i] = rbin[i]*rbin[i];
181 }
182
183 // grow pair_entropy and pair_entropy_avg array if necessary
184
185 if (atom->nmax > nmax) {
186 if (!avg_flag) {
187 memory->destroy(pair_entropy);
188 nmax = atom->nmax;
189 memory->create(pair_entropy,nmax,"entropy/atom:pair_entropy");
190 vector_atom = pair_entropy;
191 } else {
192 memory->destroy(pair_entropy);
193 memory->destroy(pair_entropy_avg);
194 nmax = atom->nmax;
195 memory->create(pair_entropy,nmax,"entropy/atom:pair_entropy");
196 memory->create(pair_entropy_avg,nmax,
197 "entropy/atom:pair_entropy_avg");
198 vector_atom = pair_entropy_avg;
199 }
200 }
201
202 // invoke occasional neighbor list build (if not perpetual)
203
204 if (!avg_flag) neighbor->build_one(list);
205
206 inum = list->inum + list->gnum;
207 ilist = list->ilist;
208 numneigh = list->numneigh;
209 firstneigh = list->firstneigh;
210
211 // Compute some constants
212 double sigmasq2=2*sigma*sigma;
213 double volume = domain->xprd * domain->yprd * domain->zprd;
214 double density = atom->natoms / volume;
215
216 // compute pair entropy for each atom in group
217 // use full neighbor list
218
219 double **x = atom->x;
220 int *mask = atom->mask;
221 double *gofr = new double[nbin];
222 double *integrand = new double[nbin];
223
224 for (ii = 0; ii < inum; ii++) {
225 i = ilist[ii];
226 if (mask[i] & groupbit) {
227 xtmp = x[i][0];
228 ytmp = x[i][1];
229 ztmp = x[i][2];
230 jlist = firstneigh[i];
231 jnum = numneigh[i];
232
233 // If local density is used, calculate it
234 if (local_flag) {
235 double neigh_cutoff = force->pair->cutforce + neighbor->skin;
236 double volume =
237 (4./3.)*MY_PI*neigh_cutoff*neigh_cutoff*neigh_cutoff;
238 density = jnum / volume;
239 }
240
241 // calculate kernel normalization
242 // Normalization of g(r)
243 double normConstantBase = 4*MY_PI*density;
244 // Normalization of gaussian
245 normConstantBase *= sqrt(2.*MY_PI)*sigma;
246 double invNormConstantBase = 1./normConstantBase;
247
248 // loop over list of all neighbors within force cutoff
249
250 // initialize gofr
251 for (int k=0;k<nbin;++k) gofr[k]=0.;
252
253 for (jj = 0; jj < jnum; jj++) {
254 j = jlist[jj];
255 j &= NEIGHMASK;
256
257 delx = xtmp - x[j][0];
258 dely = ytmp - x[j][1];
259 delz = ztmp - x[j][2];
260 rsq = delx*delx + dely*dely + delz*delz;
261 if (rsq < cutsq) {
262 // contribute to gofr
263 double r=sqrt(rsq);
264 int bin=floor(r/deltar);
265 int minbin, maxbin;
266 minbin=bin - deltabin;
267 if (minbin < 0) minbin=0;
268 if (minbin > (nbin-1)) minbin=nbin-1;
269 maxbin=bin + deltabin;
270 if (maxbin > (nbin-1)) maxbin=nbin-1;
271 for (int k=minbin;k<maxbin+1;k++) {
272 double invNormKernel=invNormConstantBase/rbinsq[k];
273 double distance = r - rbin[k];
274 gofr[k] += invNormKernel*exp(-distance*distance/sigmasq2);
275 }
276 }
277 }
278
279 // Calculate integrand
280 for (int k=0;k<nbin;++k) {
281 if (gofr[k]<1.e-10) {
282 integrand[k] = rbinsq[k];
283 } else {
284 integrand[k] = (gofr[k]*log(gofr[k])-gofr[k]+1)*rbinsq[k];
285 }
286 }
287
288 // Integrate with trapezoid rule
289 double value = 0.;
290 for (int k=1;k<nbin-1;++k) {
291 value += integrand[k];
292 }
293 value += 0.5*integrand[0];
294 value += 0.5*integrand[nbin-1];
295 value *= deltar;
296
297 pair_entropy[i] = -2*MY_PI*density*value;
298 }
299 }
300 delete [] gofr;
301 delete [] integrand;
302
303 if (avg_flag) {
304 for (ii = 0; ii < inum; ii++) {
305 i = ilist[ii];
306 if (mask[i] & groupbit) {
307 xtmp = x[i][0];
308 ytmp = x[i][1];
309 ztmp = x[i][2];
310 jlist = firstneigh[i];
311 jnum = numneigh[i];
312
313 pair_entropy_avg[i] = pair_entropy[i];
314 double counter = 1;
315
316 // loop over list of all neighbors within force cutoff
317 for (jj = 0; jj < jnum; jj++) {
318 j = jlist[jj];
319 j &= NEIGHMASK;
320
321 delx = xtmp - x[j][0];
322 dely = ytmp - x[j][1];
323 delz = ztmp - x[j][2];
324 rsq = delx*delx + dely*dely + delz*delz;
325 if (rsq < cutsq2) {
326 pair_entropy_avg[i] += pair_entropy[j];
327 counter += 1;
328 }
329 }
330 pair_entropy_avg[i] /= counter;
331 }
332 }
333 }
334 delete [] rbin;
335 delete [] rbinsq;
336 }
337
338
339 /* ----------------------------------------------------------------------
340 memory usage of local atom-based array
341 ------------------------------------------------------------------------- */
342
memory_usage()343 double ComputeEntropyAtom::memory_usage()
344 {
345 double bytes;
346 if (!avg_flag) {
347 bytes = nmax * sizeof(double);
348 } else {
349 bytes = 2 * nmax * sizeof(double);
350 }
351 return bytes;
352 }
353