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