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: C.D. Barrett, cdb333@cavs.msstate.edu
17 Copyright (C) 2013
18 ------------------------------------------------------------------------- */
19
20 #include <cmath>
21 #include <cstring>
22 #include "compute_basal_atom.h"
23 #include "atom.h"
24 #include "update.h"
25 #include "modify.h"
26 #include "neighbor.h"
27 #include "neigh_list.h"
28 #include "neigh_request.h"
29 #include "force.h"
30 #include "pair.h"
31 #include "comm.h"
32 #include "memory.h"
33 #include "error.h"
34
35 using namespace LAMMPS_NS;
36
37 /* ---------------------------------------------------------------------- */
38
ComputeBasalAtom(LAMMPS * lmp,int narg,char ** arg)39 ComputeBasalAtom::ComputeBasalAtom(LAMMPS *lmp, int narg, char **arg) :
40 Compute(lmp, narg, arg)
41 {
42 if (narg != 3) error->all(FLERR,"Illegal compute basal/atom command");
43
44 peratom_flag = 1;
45 size_peratom_cols = 3;
46
47 nmax = 0;
48 BPV = nullptr;
49 maxneigh = 0;
50 distsq = nullptr;
51 nearest = nullptr;
52 nearest_n0 = nullptr;
53 nearest_n1 = nullptr;
54 }
55
56 /* ---------------------------------------------------------------------- */
57
~ComputeBasalAtom()58 ComputeBasalAtom::~ComputeBasalAtom()
59 {
60 memory->destroy(BPV);
61 memory->destroy(distsq);
62 memory->destroy(nearest);
63 memory->destroy(nearest_n0);
64 memory->destroy(nearest_n1);
65 }
66
67 /* ---------------------------------------------------------------------- */
68
init()69 void ComputeBasalAtom::init()
70 {
71 // need an occasional full neighbor list
72
73 int irequest = neighbor->request(this,instance_me);
74 neighbor->requests[irequest]->pair = 0;
75 neighbor->requests[irequest]->compute = 1;
76 neighbor->requests[irequest]->half = 0;
77 neighbor->requests[irequest]->full = 1;
78 neighbor->requests[irequest]->occasional = 1;
79
80 int count1 = 0;
81 for (int i = 0; i < modify->ncompute; i++)
82 if (strcmp(modify->compute[i]->style,"basal/atom") == 0) count1++;
83 if (count1 > 1 && comm->me == 0)
84 error->warning(FLERR,"More than one compute basal/atom");
85 }
86
87 /* ---------------------------------------------------------------------- */
88
init_list(int,NeighList * ptr)89 void ComputeBasalAtom::init_list(int /*id*/, NeighList *ptr)
90 {
91 list = ptr;
92 }
93
94 /* ---------------------------------------------------------------------- */
95
compute_peratom()96 void ComputeBasalAtom::compute_peratom()
97 {
98 int i,j,ii,jj,k,n,inum,jnum;
99 double xtmp,ytmp,ztmp,delx,dely,delz,rsq,var5,var6,var7;
100 int *ilist,*jlist,*numneigh,**firstneigh;
101 int chi[8];
102 int value;
103 int count;
104 int k2[3];
105 int j1[3];
106 double x4[3],y4[3],z4[3],x5[3],y5[3],z5[3],x6[3],y6[3],z6[3];
107 double x7[3],y7[3],z7[3];
108
109 invoked_peratom = update->ntimestep;
110
111 // grow structure array if necessary
112
113 if (atom->nmax > nmax) {
114 memory->destroy(BPV);
115 nmax = atom->nmax;
116 memory->create(BPV,nmax,3,"basal/atom:basal");
117 array_atom = BPV;
118 }
119
120 // invoke full neighbor list (will copy or build if necessary)
121
122 neighbor->build_one(list);
123
124 inum = list->inum;
125 ilist = list->ilist;
126 numneigh = list->numneigh;
127 firstneigh = list->firstneigh;
128
129 // compute structure parameter for each atom in group
130 // use full neighbor list
131
132 double **x = atom->x;
133 int *mask = atom->mask;
134 double cutsq = force->pair->cutforce * force->pair->cutforce;
135
136 for (ii = 0; ii < inum; ii++) {
137 i = ilist[ii];
138 if (mask[i] & groupbit) {
139 xtmp = x[i][0];
140 ytmp = x[i][1];
141 ztmp = x[i][2];
142 jlist = firstneigh[i];
143 jnum = numneigh[i];
144
145 // ensure distsq and nearest arrays are long enough
146
147 if (jnum > maxneigh) {
148 memory->destroy(distsq);
149 memory->destroy(nearest);
150 memory->destroy(nearest_n0);
151 memory->destroy(nearest_n1);
152 maxneigh = jnum;
153 memory->create(distsq,maxneigh,"compute/basal/atom:distsq");
154 memory->create(nearest,maxneigh,"compute/basal/atom:nearest");
155 memory->create(nearest_n0,maxneigh,"compute/basal/atom:nearest_n0");
156 memory->create(nearest_n1,maxneigh,"compute/basal/atom:nearest_n1");
157 }
158 // neighbor selection is identical to ackland/atom algorithm
159
160 // loop over list of all neighbors within force cutoff
161 // distsq[] = distance sq to each
162 // nearest[] = atom indices of neighbors
163
164 n = 0;
165 for (jj = 0; jj < jnum; jj++) {
166 j = jlist[jj];
167 j &= NEIGHMASK;
168
169 delx = xtmp - x[j][0];
170 dely = ytmp - x[j][1];
171 delz = ztmp - x[j][2];
172 rsq = delx*delx + dely*dely + delz*delz;
173 if (rsq < cutsq) {
174 distsq[n] = rsq;
175 nearest[n++] = j;
176 }
177 }
178
179 // Select 6 nearest neighbors
180
181 select2(6,n,distsq,nearest);
182
183 // Mean squared separation
184
185 double r0_sq = 0.0;
186 for (j = 0; j < 6; j++) r0_sq += distsq[j];
187 r0_sq /= 6.0;
188
189 // n0 near neighbors with: distsq<1.45*r0_sq
190 // n1 near neighbors with: distsq<1.55*r0_sq
191
192 double n0_dist_sq = 1.45*r0_sq,
193 n1_dist_sq = 1.55*r0_sq;
194 int n0 = 0, n1 = 0;
195 for (j = 0; j < n; j++) {
196 if (distsq[j] < n1_dist_sq) {
197 nearest_n1[n1++] = nearest[j];
198 if (distsq[j] < n0_dist_sq) {
199 nearest_n0[n0++] = nearest[j];
200 }
201 }
202 }
203
204 // Evaluate all angles <(r_ij,rik) forall n0 particles with: distsq<1.45*r0_sq
205 double bond_angle;
206 double norm_j, norm_k;
207 chi[0] = chi[1] = chi[2] = chi[3] = chi[4] = chi[5] = chi[6] = chi[7] = 0;
208 double x_ij, y_ij, z_ij, x_ik, y_ik, z_ik, xmean5, ymean5, zmean5,
209 xmean6, ymean6, zmean6, xmean7, ymean7, zmean7;
210 double *x3 = new double[n0];
211 double *y3 = new double[n0];
212 double *z3 = new double[n0];
213 for (j = 0; j < n0; j++) {
214 x_ij = x[i][0]-x[nearest_n0[j]][0];
215 y_ij = x[i][1]-x[nearest_n0[j]][1];
216 z_ij = x[i][2]-x[nearest_n0[j]][2];
217 norm_j = sqrt (x_ij*x_ij + y_ij*y_ij + z_ij*z_ij);
218 if (norm_j <= 0.) {continue;}
219 for (k = j+1; k < n0; k++) {
220 x_ik = x[i][0]-x[nearest_n0[k]][0];
221 y_ik = x[i][1]-x[nearest_n0[k]][1];
222 z_ik = x[i][2]-x[nearest_n0[k]][2];
223 norm_k = sqrt (x_ik*x_ik + y_ik*y_ik + z_ik*z_ik);
224 if (norm_k <= 0.) {continue;}
225 bond_angle = (x_ij*x_ik + y_ij*y_ik + z_ij*z_ik) / (norm_j*norm_k);
226 //find all bond angles that are about 180 degrees
227 if (-1. <= bond_angle && bond_angle < -0.945) {
228 x3[chi[0]] = x_ik - x_ij;
229 y3[chi[0]] = y_ik - y_ij;
230 z3[chi[0]] = z_ik - z_ij;
231 chi[0]++;
232 }
233 }
234 }
235 // for atoms that have 2 or 3 ~180 bond angles:
236 if (2 == chi[0] || 3 == chi[0]) {
237 count = value = 0;
238 if (chi[0] == 2) {
239 k2[0] = 0;
240 j1[0] = 1;
241 }
242 else {
243 k2[0] = 0;
244 k2[1] = 0;
245 k2[2] = 1;
246 j1[0]=1;
247 j1[1]=2;
248 j1[2]=2;
249 }
250 xmean5 = ymean5 = zmean5 = xmean6 = ymean6 = zmean6 = xmean7 = ymean7 = zmean7 = 0.0;
251 for (j = 0; j < chi[0]; j++) {
252 for (k = j+1; k < chi[0]; k++) {
253 //get cross products
254 x4[count] = y3[j1[count]]*z3[k2[count]]-y3[k2[count]]*z3[j1[count]];
255 y4[count] = z3[j1[count]]*x3[k2[count]]-z3[k2[count]]*x3[j1[count]];
256 z4[count] = x3[j1[count]]*y3[k2[count]]-x3[k2[count]]*y3[j1[count]];
257 //get all sign combinations of cross products
258 x5[count] = x4[count]*copysign(1.0,x4[count]);
259 y5[count] = y4[count]*copysign(1.0,x4[count]);
260 z5[count] = z4[count]*copysign(1.0,x4[count]);
261 x6[count] = x4[count]*copysign(1.0,y4[count]);
262 y6[count] = y4[count]*copysign(1.0,y4[count]);
263 z6[count] = z4[count]*copysign(1.0,y4[count]);
264 x7[count] = x4[count]*copysign(1.0,z4[count]);
265 y7[count] = y4[count]*copysign(1.0,z4[count]);
266 z7[count] = z4[count]*copysign(1.0,z4[count]);
267 //get average cross products
268 xmean5 += x5[count];
269 ymean5 += y5[count];
270 zmean5 += z5[count];
271 xmean6 += x6[count];
272 ymean6 += y6[count];
273 zmean6 += z6[count];
274 xmean7 += x7[count];
275 ymean7 += y7[count];
276 zmean6 += z7[count];
277 count++;
278 }
279 }
280 if (count > 0) {
281 xmean5 /= count;
282 xmean6 /= count;
283 xmean7 /= count;
284 ymean5 /= count;
285 ymean6 /= count;
286 ymean7 /= count;
287 zmean5 /= count;
288 zmean6 /= count;
289 zmean7 /= count;
290 }
291 var5 = var6 = var7 = 0.0;
292 //find standard deviations
293 for (j=0;j<count;j++) {
294 var5 = var5 + x5[j]*x5[j]-2*x5[j]*xmean5+xmean5*xmean5+y5[j]*y5[j]-2*y5[j]*ymean5+ymean5*ymean5+z5[j]*z5[j]-2*z5[j]*zmean5+zmean5*zmean5;
295 var6 = var6 + x6[j]*x6[j]-2*x6[j]*xmean6+xmean6*xmean6+y6[j]*y6[j]-2*y6[j]*ymean6+ymean6*ymean6+z6[j]*z6[j]-2*z6[j]*zmean6+zmean6*zmean6;
296 var7 = var7 + x7[j]*x7[j]-2*x7[j]*xmean7+xmean7*xmean7+y7[j]*y7[j]-2*y7[j]*ymean7+ymean7*ymean7+z7[j]*z7[j]-2*z7[j]*zmean7+zmean7*zmean7;
297 }
298 //select sign combination with minimum standard deviation
299 if (var5 < var6) {
300 if (var5 < var7) { value = 0;}
301 else {value = 2;}
302 }
303 else if (var6 < var7) {value = 1;}
304 else {value = 2;}
305 //BPV is average of cross products of all neighbor vectors which are part of 180 degree angles
306 BPV[i][0] = 0;
307 BPV[i][1] = 0;
308 BPV[i][2] = 0;
309 for (k=0;k<count;k++) {
310 if (value == 0) {
311 BPV[i][0] = BPV[i][0]+x5[k];
312 BPV[i][1] = BPV[i][1]+y5[k];
313 BPV[i][2] = BPV[i][2]+z5[k];
314 }
315 else if (value == 1) {
316 BPV[i][0] = BPV[i][0]+x6[k];
317 BPV[i][1] = BPV[i][1]+y6[k];
318 BPV[i][2] = BPV[i][2]+z6[k];
319 }
320 else {
321 BPV[i][0] = BPV[i][0]+x7[k];
322 BPV[i][1] = BPV[i][1]+y7[k];
323 BPV[i][2] = BPV[i][2]+z7[k];
324 }
325 }
326 }
327 //for atoms with more than three 180 degree bond angles:
328 else if (chi[0] > 3) {
329 double x44[3], y44[3], z44[3], S0;
330 int l, m;
331 count = value = 0;
332 S0 = 100000;
333 k2[0] = 0;
334 k2[1] = 0;
335 k2[2] = 1;
336 j1[0]=1;
337 j1[1]=2;
338 j1[2]=2;
339 //algorithm is as above, but now all combinations of three 180 degree angles are compared, and the combination with minimum standard deviation is chosen
340 for (j=0; j<chi[0]; j++) {
341 for (k=j+1; k<chi[0]; k++) {
342 for (l=k+1; l<chi[0]; l++) {
343 if (k >= chi[0] || l >= chi[0]) continue;
344 //get unique combination of three neighbor vectors
345 x4[0] = x3[j];
346 x4[1] = x3[k];
347 x4[2] = x3[l];
348 y4[0] = y3[j];
349 y4[1] = y3[k];
350 y4[2] = y3[l];
351 z4[0] = z3[j];
352 z4[1] = z3[k];
353 z4[2] = z3[l];
354 xmean5 = ymean5 = zmean5 = xmean6 = ymean6 = zmean6 = xmean7 = ymean7 = zmean7 = 0;
355 for (m=0;m<3;m++) {
356 //get cross products
357 x44[m] = y4[j1[m]]*z4[k2[m]]-y4[k2[m]]*z4[j1[m]];
358 y44[m] = z4[j1[m]]*x4[k2[m]]-z4[k2[m]]*x4[j1[m]];
359 z44[m] = x4[j1[m]]*y4[k2[m]]-x4[k2[m]]*y4[j1[m]];
360 x5[m] = x44[m]*copysign(1.0,x44[m]);
361 y5[m] = y44[m]*copysign(1.0,x44[m]);
362 z5[m] = z44[m]*copysign(1.0,x44[m]);
363 x6[m] = x44[m]*copysign(1.0,y44[m]);
364 y6[m] = y44[m]*copysign(1.0,y44[m]);
365 z6[m] = z44[m]*copysign(1.0,y44[m]);
366 x7[m] = x44[m]*copysign(1.0,z44[m]);
367 y7[m] = y44[m]*copysign(1.0,z44[m]);
368 z7[m] = z44[m]*copysign(1.0,z44[m]);
369 //get average cross products
370 xmean5 = xmean5 + x5[m];
371 ymean5 = ymean5 + y5[m];
372 zmean5 = zmean5 + z5[m];
373 xmean6 = xmean6 + x6[m];
374 ymean6 = ymean6 + y6[m];
375 zmean6 = zmean6 + z6[m];
376 xmean7 = xmean7 + x7[m];
377 ymean7 = ymean7 + y7[m];
378 zmean6 = zmean6 + z7[m];
379 }
380 xmean5 = xmean5/3;
381 xmean6 = xmean6/3;
382 xmean7 = xmean7/3;
383 ymean5 = ymean5/3;
384 ymean6 = ymean6/3;
385 ymean7 = ymean7/3;
386 zmean5 = zmean5/3;
387 zmean6 = zmean6/3;
388 zmean7 = zmean7/3;
389 var5 = var6 = var7 = 0;
390 //get standard deviations
391 for (m=0;m<3;m++) {
392 var5 = var5 + x5[m]*x5[m]-2*x5[m]*xmean5+xmean5*xmean5+y5[m]*y5[m]-2*y5[m]*ymean5+ymean5*ymean5+z5[m]*z5[m]-2*z5[m]*zmean5+zmean5*zmean5;
393 var6 = var6 + x6[m]*x6[m]-2*x6[m]*xmean6+xmean6*xmean6+y6[m]*y6[m]-2*y6[m]*ymean6+ymean6*ymean6+z6[m]*z6[m]-2*z6[m]*zmean6+zmean6*zmean6;
394 var7 = var7 + x7[m]*x7[m]-2*x7[m]*xmean7+xmean7*xmean7+y7[m]*y7[m]-2*y7[m]*ymean7+ymean7*ymean7+z7[m]*z7[m]-2*z7[m]*zmean7+zmean7*zmean7;
395 }
396 //choose minimum standard deviation
397 if (var5 < S0) {
398 S0 = var5;
399 BPV[i][0] = (x5[0]+x5[1]+x5[2])/3;
400 BPV[i][1] = (y5[0]+y5[1]+x5[2])/3;
401 BPV[i][2] = (z5[0]+z5[1]+z5[2])/3;
402 }
403 if (var6 < S0) {
404 S0 = var6;
405 BPV[i][0] = (x6[0]+x6[1]+x6[2])/3;
406 BPV[i][1] = (y6[0]+y6[1]+x6[2])/3;
407 BPV[i][2] = (z6[0]+z6[1]+z6[2])/3;
408 }
409 if (var7 < S0) {
410 S0 = var7;
411 BPV[i][0] = (x7[0]+x7[1]+x7[2])/3;
412 BPV[i][1] = (y7[0]+y7[1]+x7[2])/3;
413 BPV[i][2] = (z7[0]+z7[1]+z7[2])/3;
414 }
415 }
416 }
417 }
418 //if there are less than two ~180 degree bond angles, the algorithm returns null
419 } else BPV[i][0] = BPV[i][1] = BPV[i][2] = 0.0;
420
421 delete[] x3;
422 delete[] y3;
423 delete[] z3;
424
425 //normalize BPV:
426 double Mag = sqrt(BPV[i][0]*BPV[i][0] +
427 BPV[i][1]*BPV[i][1] + BPV[i][2]*BPV[i][2]);
428 if (Mag > 0) {
429 BPV[i][0] = BPV[i][0]/Mag;
430 BPV[i][1] = BPV[i][1]/Mag;
431 BPV[i][2] = BPV[i][2]/Mag;
432 }
433 } else BPV[i][0] = BPV[i][1] = BPV[i][2] = 0.0;
434 }
435 }
436 /* ----------------------------------------------------------------------
437 2 select routines from Numerical Recipes (slightly modified)
438 find k smallest values in array of length n
439 2nd routine sorts auxiliary array at same time
440 ------------------------------------------------------------------------- */
441
442 #define SWAP(a,b) tmp = a; a = b; b = tmp;
443 #define ISWAP(a,b) itmp = a; a = b; b = itmp;
444
select(int k,int n,double * arr)445 void ComputeBasalAtom::select(int k, int n, double *arr)
446 {
447 int i,ir,j,l,mid;
448 double a,tmp;
449
450 arr--;
451 l = 1;
452 ir = n;
453 for (;;) {
454 if (ir <= l+1) {
455 if (ir == l+1 && arr[ir] < arr[l]) {
456 SWAP(arr[l],arr[ir])
457 }
458 return;
459 } else {
460 mid=(l+ir) >> 1;
461 SWAP(arr[mid],arr[l+1])
462 if (arr[l] > arr[ir]) {
463 SWAP(arr[l],arr[ir])
464 }
465 if (arr[l+1] > arr[ir]) {
466 SWAP(arr[l+1],arr[ir])
467 }
468 if (arr[l] > arr[l+1]) {
469 SWAP(arr[l],arr[l+1])
470 }
471 i = l+1;
472 j = ir;
473 a = arr[l+1];
474 for (;;) {
475 do i++; while (arr[i] < a);
476 do j--; while (arr[j] > a);
477 if (j < i) break;
478 SWAP(arr[i],arr[j])
479 }
480 arr[l+1] = arr[j];
481 arr[j] = a;
482 if (j >= k) ir = j-1;
483 if (j <= k) l = i;
484 }
485 }
486 }
487
488 /* ---------------------------------------------------------------------- */
489
select2(int k,int n,double * arr,int * iarr)490 void ComputeBasalAtom::select2(int k, int n, double *arr, int *iarr)
491 {
492 int i,ir,j,l,mid,ia,itmp;
493 double a,tmp;
494
495 arr--;
496 iarr--;
497 l = 1;
498 ir = n;
499 for (;;) {
500 if (ir <= l+1) {
501 if (ir == l+1 && arr[ir] < arr[l]) {
502 SWAP(arr[l],arr[ir])
503 ISWAP(iarr[l],iarr[ir])
504 }
505 return;
506 } else {
507 mid=(l+ir) >> 1;
508 SWAP(arr[mid],arr[l+1])
509 ISWAP(iarr[mid],iarr[l+1])
510 if (arr[l] > arr[ir]) {
511 SWAP(arr[l],arr[ir])
512 ISWAP(iarr[l],iarr[ir])
513 }
514 if (arr[l+1] > arr[ir]) {
515 SWAP(arr[l+1],arr[ir])
516 ISWAP(iarr[l+1],iarr[ir])
517 }
518 if (arr[l] > arr[l+1]) {
519 SWAP(arr[l],arr[l+1])
520 ISWAP(iarr[l],iarr[l+1])
521 }
522 i = l+1;
523 j = ir;
524 a = arr[l+1];
525 ia = iarr[l+1];
526 for (;;) {
527 do i++; while (arr[i] < a);
528 do j--; while (arr[j] > a);
529 if (j < i) break;
530 SWAP(arr[i],arr[j])
531 ISWAP(iarr[i],iarr[j])
532 }
533 arr[l+1] = arr[j];
534 arr[j] = a;
535 iarr[l+1] = iarr[j];
536 iarr[j] = ia;
537 if (j >= k) ir = j-1;
538 if (j <= k) l = i;
539 }
540 }
541 }
542
543 /* ----------------------------------------------------------------------
544 memory usage of local atom-based array
545 ------------------------------------------------------------------------- */
546
memory_usage()547 double ComputeBasalAtom::memory_usage()
548 {
549 double bytes = 3*nmax * sizeof(double);
550 return bytes;
551 }
552