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 "angle.h"
16 #include "atom.h"
17 #include "comm.h"
18 #include "force.h"
19 #include "math_const.h"
20 #include "suffix.h"
21 #include "atom_masks.h"
22 #include "memory.h"
23 #include "error.h"
24
25 using namespace LAMMPS_NS;
26 using namespace MathConst;
27
28 /* ---------------------------------------------------------------------- */
29
Angle(LAMMPS * lmp)30 Angle::Angle(LAMMPS *lmp) : Pointers(lmp)
31 {
32 energy = 0.0;
33 virial[0] = virial[1] = virial[2] = virial[3] = virial[4] = virial[5] = 0.0;
34 writedata = 1;
35
36 allocated = 0;
37 suffix_flag = Suffix::NONE;
38
39 maxeatom = maxvatom = maxcvatom = 0;
40 eatom = nullptr;
41 vatom = nullptr;
42 cvatom = nullptr;
43 setflag = nullptr;
44 centroidstressflag = CENTROID_AVAIL;
45
46 execution_space = Host;
47 datamask_read = ALL_MASK;
48 datamask_modify = ALL_MASK;
49
50 copymode = 0;
51 }
52
53 /* ---------------------------------------------------------------------- */
54
~Angle()55 Angle::~Angle()
56 {
57 if (copymode) return;
58
59 memory->destroy(eatom);
60 memory->destroy(vatom);
61 memory->destroy(cvatom);
62 }
63
64 /* ----------------------------------------------------------------------
65 check if all coeffs are set
66 ------------------------------------------------------------------------- */
67
init()68 void Angle::init()
69 {
70 if (!allocated && atom->nangletypes)
71 error->all(FLERR,"Angle coeffs are not set");
72 for (int i = 1; i <= atom->nangletypes; i++)
73 if (setflag[i] == 0) error->all(FLERR,"All angle coeffs are not set");
74
75 init_style();
76 }
77
78 /* ----------------------------------------------------------------------
79 setup for energy, virial computation
80 see integrate::ev_set() for bitwise settings of eflag/vflag
81 set the following flags, values are otherwise set to 0:
82 evflag != 0 if any bits of eflag or vflag are set
83 eflag_global != 0 if ENERGY_GLOBAL bit of eflag set
84 eflag_atom != 0 if ENERGY_ATOM bit of eflag set
85 eflag_either != 0 if eflag_global or eflag_atom is set
86 vflag_global != 0 if VIRIAL_PAIR or VIRIAL_FDOTR bit of vflag set
87 vflag_atom != 0 if VIRIAL_ATOM bit of vflag set
88 vflag_atom != 0 if VIRIAL_CENTROID bit of vflag set
89 and centroidstressflag != CENTROID_AVAIL
90 cvflag_atom != 0 if VIRIAL_CENTROID bit of vflag set
91 and centroidstressflag = CENTROID_AVAIL
92 vflag_either != 0 if any of vflag_global, vflag_atom, cvflag_atom is set
93 ------------------------------------------------------------------------- */
94
ev_setup(int eflag,int vflag,int alloc)95 void Angle::ev_setup(int eflag, int vflag, int alloc)
96 {
97 int i,n;
98
99 evflag = 1;
100
101 eflag_either = eflag;
102 eflag_global = eflag & ENERGY_GLOBAL;
103 eflag_atom = eflag & ENERGY_ATOM;
104
105 vflag_global = vflag & (VIRIAL_PAIR | VIRIAL_FDOTR);
106 vflag_atom = vflag & VIRIAL_ATOM;
107 if (vflag & VIRIAL_CENTROID && centroidstressflag != CENTROID_AVAIL)
108 vflag_atom = 1;
109 cvflag_atom = 0;
110 if (vflag & VIRIAL_CENTROID && centroidstressflag == CENTROID_AVAIL)
111 cvflag_atom = 1;
112 vflag_either = vflag_global || vflag_atom || cvflag_atom;
113
114 // reallocate per-atom arrays if necessary
115
116 if (eflag_atom && atom->nmax > maxeatom) {
117 maxeatom = atom->nmax;
118 if (alloc) {
119 memory->destroy(eatom);
120 memory->create(eatom,comm->nthreads*maxeatom,"angle:eatom");
121 }
122 }
123 if (vflag_atom && atom->nmax > maxvatom) {
124 maxvatom = atom->nmax;
125 if (alloc) {
126 memory->destroy(vatom);
127 memory->create(vatom,comm->nthreads*maxvatom,6,"angle:vatom");
128 }
129 }
130 if (cvflag_atom && atom->nmax > maxcvatom) {
131 maxcvatom = atom->nmax;
132 if (alloc) {
133 memory->destroy(cvatom);
134 memory->create(cvatom,comm->nthreads*maxcvatom,9,"angle:cvatom");
135 }
136 }
137
138 // zero accumulators
139
140 if (eflag_global) energy = 0.0;
141 if (vflag_global) for (i = 0; i < 6; i++) virial[i] = 0.0;
142 if (eflag_atom && alloc) {
143 n = atom->nlocal;
144 if (force->newton_bond) n += atom->nghost;
145 for (i = 0; i < n; i++) eatom[i] = 0.0;
146 }
147 if (vflag_atom && alloc) {
148 n = atom->nlocal;
149 if (force->newton_bond) n += atom->nghost;
150 for (i = 0; i < n; i++) {
151 vatom[i][0] = 0.0;
152 vatom[i][1] = 0.0;
153 vatom[i][2] = 0.0;
154 vatom[i][3] = 0.0;
155 vatom[i][4] = 0.0;
156 vatom[i][5] = 0.0;
157 }
158 }
159 if (cvflag_atom && alloc) {
160 n = atom->nlocal;
161 if (force->newton_bond) n += atom->nghost;
162 for (i = 0; i < n; i++) {
163 cvatom[i][0] = 0.0;
164 cvatom[i][1] = 0.0;
165 cvatom[i][2] = 0.0;
166 cvatom[i][3] = 0.0;
167 cvatom[i][4] = 0.0;
168 cvatom[i][5] = 0.0;
169 cvatom[i][6] = 0.0;
170 cvatom[i][7] = 0.0;
171 cvatom[i][8] = 0.0;
172 cvatom[i][9] = 0.0;
173 }
174 }
175 }
176
177 /* ----------------------------------------------------------------------
178 tally energy and virial into global and per-atom accumulators
179 virial = r1F1 + r2F2 + r3F3 = (r1-r2) F1 + (r3-r2) F3 = del1*f1 + del2*f3
180 ------------------------------------------------------------------------- */
181
ev_tally(int i,int j,int k,int nlocal,int newton_bond,double eangle,double * f1,double * f3,double delx1,double dely1,double delz1,double delx2,double dely2,double delz2)182 void Angle::ev_tally(int i, int j, int k, int nlocal, int newton_bond,
183 double eangle, double *f1, double *f3,
184 double delx1, double dely1, double delz1,
185 double delx2, double dely2, double delz2)
186 {
187 double eanglethird,v[6];
188
189 if (eflag_either) {
190 if (eflag_global) {
191 if (newton_bond) energy += eangle;
192 else {
193 eanglethird = THIRD*eangle;
194 if (i < nlocal) energy += eanglethird;
195 if (j < nlocal) energy += eanglethird;
196 if (k < nlocal) energy += eanglethird;
197 }
198 }
199 if (eflag_atom) {
200 eanglethird = THIRD*eangle;
201 if (newton_bond || i < nlocal) eatom[i] += eanglethird;
202 if (newton_bond || j < nlocal) eatom[j] += eanglethird;
203 if (newton_bond || k < nlocal) eatom[k] += eanglethird;
204 }
205 }
206
207 if (vflag_either) {
208 v[0] = delx1*f1[0] + delx2*f3[0];
209 v[1] = dely1*f1[1] + dely2*f3[1];
210 v[2] = delz1*f1[2] + delz2*f3[2];
211 v[3] = delx1*f1[1] + delx2*f3[1];
212 v[4] = delx1*f1[2] + delx2*f3[2];
213 v[5] = dely1*f1[2] + dely2*f3[2];
214
215 if (vflag_global) {
216 if (newton_bond) {
217 virial[0] += v[0];
218 virial[1] += v[1];
219 virial[2] += v[2];
220 virial[3] += v[3];
221 virial[4] += v[4];
222 virial[5] += v[5];
223 } else {
224 if (i < nlocal) {
225 virial[0] += THIRD*v[0];
226 virial[1] += THIRD*v[1];
227 virial[2] += THIRD*v[2];
228 virial[3] += THIRD*v[3];
229 virial[4] += THIRD*v[4];
230 virial[5] += THIRD*v[5];
231 }
232 if (j < nlocal) {
233 virial[0] += THIRD*v[0];
234 virial[1] += THIRD*v[1];
235 virial[2] += THIRD*v[2];
236 virial[3] += THIRD*v[3];
237 virial[4] += THIRD*v[4];
238 virial[5] += THIRD*v[5];
239 }
240 if (k < nlocal) {
241 virial[0] += THIRD*v[0];
242 virial[1] += THIRD*v[1];
243 virial[2] += THIRD*v[2];
244 virial[3] += THIRD*v[3];
245 virial[4] += THIRD*v[4];
246 virial[5] += THIRD*v[5];
247 }
248 }
249 }
250
251 if (vflag_atom) {
252 if (newton_bond || i < nlocal) {
253 vatom[i][0] += THIRD*v[0];
254 vatom[i][1] += THIRD*v[1];
255 vatom[i][2] += THIRD*v[2];
256 vatom[i][3] += THIRD*v[3];
257 vatom[i][4] += THIRD*v[4];
258 vatom[i][5] += THIRD*v[5];
259 }
260 if (newton_bond || j < nlocal) {
261 vatom[j][0] += THIRD*v[0];
262 vatom[j][1] += THIRD*v[1];
263 vatom[j][2] += THIRD*v[2];
264 vatom[j][3] += THIRD*v[3];
265 vatom[j][4] += THIRD*v[4];
266 vatom[j][5] += THIRD*v[5];
267 }
268 if (newton_bond || k < nlocal) {
269 vatom[k][0] += THIRD*v[0];
270 vatom[k][1] += THIRD*v[1];
271 vatom[k][2] += THIRD*v[2];
272 vatom[k][3] += THIRD*v[3];
273 vatom[k][4] += THIRD*v[4];
274 vatom[k][5] += THIRD*v[5];
275 }
276 }
277 }
278
279 // per-atom centroid virial
280 if (cvflag_atom) {
281
282 // r0 = (r1+r2+r3)/3
283 // rij = ri-rj
284 // total virial = r10*f1 + r20*f2 + r30*f3
285 // del1: r12
286 // del2: r32
287
288 if (newton_bond || i < nlocal) {
289 double a1[3];
290
291 // a1 = r10 = (2*r12 - r32)/3
292 a1[0] = THIRD*(2*delx1-delx2);
293 a1[1] = THIRD*(2*dely1-dely2);
294 a1[2] = THIRD*(2*delz1-delz2);
295
296 cvatom[i][0] += a1[0]*f1[0];
297 cvatom[i][1] += a1[1]*f1[1];
298 cvatom[i][2] += a1[2]*f1[2];
299 cvatom[i][3] += a1[0]*f1[1];
300 cvatom[i][4] += a1[0]*f1[2];
301 cvatom[i][5] += a1[1]*f1[2];
302 cvatom[i][6] += a1[1]*f1[0];
303 cvatom[i][7] += a1[2]*f1[0];
304 cvatom[i][8] += a1[2]*f1[1];
305 }
306 if (newton_bond || j < nlocal) {
307 double a2[3];
308 double f2[3];
309
310 // a2 = r20 = ( -r12 - r32)/3
311 a2[0] = THIRD*(-delx1-delx2);
312 a2[1] = THIRD*(-dely1-dely2);
313 a2[2] = THIRD*(-delz1-delz2);
314
315 f2[0] = - f1[0] - f3[0];
316 f2[1] = - f1[1] - f3[1];
317 f2[2] = - f1[2] - f3[2];
318
319 cvatom[j][0] += a2[0]*f2[0];
320 cvatom[j][1] += a2[1]*f2[1];
321 cvatom[j][2] += a2[2]*f2[2];
322 cvatom[j][3] += a2[0]*f2[1];
323 cvatom[j][4] += a2[0]*f2[2];
324 cvatom[j][5] += a2[1]*f2[2];
325 cvatom[j][6] += a2[1]*f2[0];
326 cvatom[j][7] += a2[2]*f2[0];
327 cvatom[j][8] += a2[2]*f2[1];
328 }
329 if (newton_bond || k < nlocal) {
330 double a3[3];
331
332 // a3 = r30 = ( -r12 + 2*r32)/3
333 a3[0] = THIRD*(-delx1+2*delx2);
334 a3[1] = THIRD*(-dely1+2*dely2);
335 a3[2] = THIRD*(-delz1+2*delz2);
336
337 cvatom[k][0] += a3[0]*f3[0];
338 cvatom[k][1] += a3[1]*f3[1];
339 cvatom[k][2] += a3[2]*f3[2];
340 cvatom[k][3] += a3[0]*f3[1];
341 cvatom[k][4] += a3[0]*f3[2];
342 cvatom[k][5] += a3[1]*f3[2];
343 cvatom[k][6] += a3[1]*f3[0];
344 cvatom[k][7] += a3[2]*f3[0];
345 cvatom[k][8] += a3[2]*f3[1];
346 }
347 }
348 }
349
350 /* ---------------------------------------------------------------------- */
351
memory_usage()352 double Angle::memory_usage()
353 {
354 double bytes = (double)comm->nthreads*maxeatom * sizeof(double);
355 bytes += (double)comm->nthreads*maxvatom*6 * sizeof(double);
356 bytes += (double)comm->nthreads*maxcvatom*9 * sizeof(double);
357 return bytes;
358 }
359