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_contact_atom.h"
16 #include <cstring>
17 #include "atom.h"
18 #include "update.h"
19 #include "modify.h"
20 #include "neighbor.h"
21 #include "neigh_list.h"
22 #include "neigh_request.h"
23 #include "force.h"
24 #include "comm.h"
25 #include "memory.h"
26 #include "error.h"
27
28 using namespace LAMMPS_NS;
29
30 /* ---------------------------------------------------------------------- */
31
ComputeContactAtom(LAMMPS * lmp,int narg,char ** arg)32 ComputeContactAtom::ComputeContactAtom(LAMMPS *lmp, int narg, char **arg) :
33 Compute(lmp, narg, arg),
34 contact(nullptr)
35 {
36 if (narg != 3) error->all(FLERR,"Illegal compute contact/atom command");
37
38 peratom_flag = 1;
39 size_peratom_cols = 0;
40 comm_reverse = 1;
41
42 nmax = 0;
43
44 // error checks
45
46 if (!atom->sphere_flag)
47 error->all(FLERR,"Compute contact/atom requires atom style sphere");
48 }
49
50 /* ---------------------------------------------------------------------- */
51
~ComputeContactAtom()52 ComputeContactAtom::~ComputeContactAtom()
53 {
54 memory->destroy(contact);
55 }
56
57 /* ---------------------------------------------------------------------- */
58
init()59 void ComputeContactAtom::init()
60 {
61 if (force->pair == nullptr)
62 error->all(FLERR,"Compute contact/atom requires a pair style be defined");
63
64 int count = 0;
65 for (int i = 0; i < modify->ncompute; i++)
66 if (strcmp(modify->compute[i]->style,"contact/atom") == 0) count++;
67 if (count > 1 && comm->me == 0)
68 error->warning(FLERR,"More than one compute contact/atom");
69
70 // need an occasional neighbor list
71
72 int irequest = neighbor->request(this,instance_me);
73 neighbor->requests[irequest]->size = 1;
74 neighbor->requests[irequest]->pair = 0;
75 neighbor->requests[irequest]->compute = 1;
76 neighbor->requests[irequest]->occasional = 1;
77 }
78
79 /* ---------------------------------------------------------------------- */
80
init_list(int,NeighList * ptr)81 void ComputeContactAtom::init_list(int /*id*/, NeighList *ptr)
82 {
83 list = ptr;
84 }
85
86 /* ---------------------------------------------------------------------- */
87
compute_peratom()88 void ComputeContactAtom::compute_peratom()
89 {
90 int i,j,ii,jj,inum,jnum;
91 double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
92 double radi,radsum,radsumsq;
93 int *ilist,*jlist,*numneigh,**firstneigh;
94
95 invoked_peratom = update->ntimestep;
96
97 // grow contact array if necessary
98
99 if (atom->nmax > nmax) {
100 memory->destroy(contact);
101 nmax = atom->nmax;
102 memory->create(contact,nmax,"contact/atom:contact");
103 vector_atom = contact;
104 }
105
106 // invoke neighbor list (will copy or build if necessary)
107
108 neighbor->build_one(list);
109
110 inum = list->inum;
111 ilist = list->ilist;
112 numneigh = list->numneigh;
113 firstneigh = list->firstneigh;
114
115 // compute number of contacts for each atom in group
116 // contact if distance <= sum of radii
117 // tally for both I and J
118
119 double **x = atom->x;
120 double *radius = atom->radius;
121 int *mask = atom->mask;
122 int nlocal = atom->nlocal;
123 int nall = nlocal + atom->nghost;
124
125 for (i = 0; i < nall; i++) contact[i] = 0.0;
126
127 for (ii = 0; ii < inum; ii++) {
128 i = ilist[ii];
129 if (mask[i] & groupbit) {
130 xtmp = x[i][0];
131 ytmp = x[i][1];
132 ztmp = x[i][2];
133 radi = radius[i];
134 jlist = firstneigh[i];
135 jnum = numneigh[i];
136
137 for (jj = 0; jj < jnum; jj++) {
138 j = jlist[jj];
139 j &= NEIGHMASK;
140
141 delx = xtmp - x[j][0];
142 dely = ytmp - x[j][1];
143 delz = ztmp - x[j][2];
144 rsq = delx*delx + dely*dely + delz*delz;
145 radsum = radi + radius[j];
146 radsumsq = radsum*radsum;
147 if (rsq <= radsumsq) {
148 contact[i] += 1.0;
149 contact[j] += 1.0;
150 }
151 }
152 }
153 }
154
155 // communicate ghost atom counts between neighbor procs if necessary
156
157 if (force->newton_pair) comm->reverse_comm_compute(this);
158 }
159
160 /* ---------------------------------------------------------------------- */
161
pack_reverse_comm(int n,int first,double * buf)162 int ComputeContactAtom::pack_reverse_comm(int n, int first, double *buf)
163 {
164 int i,m,last;
165
166 m = 0;
167 last = first + n;
168 for (i = first; i < last; i++)
169 buf[m++] = contact[i];
170 return m;
171 }
172
173 /* ---------------------------------------------------------------------- */
174
unpack_reverse_comm(int n,int * list,double * buf)175 void ComputeContactAtom::unpack_reverse_comm(int n, int *list, double *buf)
176 {
177 int i,j,m;
178
179 m = 0;
180 for (i = 0; i < n; i++) {
181 j = list[i];
182 contact[j] += buf[m++];
183 }
184 }
185
186 /* ----------------------------------------------------------------------
187 memory usage of local atom-based array
188 ------------------------------------------------------------------------- */
189
memory_usage()190 double ComputeContactAtom::memory_usage()
191 {
192 double bytes = (double)nmax * sizeof(double);
193 return bytes;
194 }
195