1 /* ----------------------------------------------------------------------
2 This is the
3
4 ██╗ ██╗ ██████╗ ██████╗ ██████╗ ██╗ ██╗████████╗███████╗
5 ██║ ██║██╔════╝ ██╔════╝ ██╔════╝ ██║ ██║╚══██╔══╝██╔════╝
6 ██║ ██║██║ ███╗██║ ███╗██║ ███╗███████║ ██║ ███████╗
7 ██║ ██║██║ ██║██║ ██║██║ ██║██╔══██║ ██║ ╚════██║
8 ███████╗██║╚██████╔╝╚██████╔╝╚██████╔╝██║ ██║ ██║ ███████║
9 ╚══════╝╚═╝ ╚═════╝ ╚═════╝ ╚═════╝ ╚═╝ ╚═╝ ╚═╝ ╚══════╝®
10
11 DEM simulation engine, released by
12 DCS Computing Gmbh, Linz, Austria
13 http://www.dcs-computing.com, office@dcs-computing.com
14
15 LIGGGHTS® is part of CFDEM®project:
16 http://www.liggghts.com | http://www.cfdem.com
17
18 Core developer and main author:
19 Christoph Kloss, christoph.kloss@dcs-computing.com
20
21 LIGGGHTS® is open-source, distributed under the terms of the GNU Public
22 License, version 2 or later. It is distributed in the hope that it will
23 be useful, but WITHOUT ANY WARRANTY; without even the implied warranty
24 of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. You should have
25 received a copy of the GNU General Public License along with LIGGGHTS®.
26 If not, see http://www.gnu.org/licenses . See also top-level README
27 and LICENSE files.
28
29 LIGGGHTS® and CFDEM® are registered trade marks of DCS Computing GmbH,
30 the producer of the LIGGGHTS® software and the CFDEM®coupling software
31 See http://www.cfdem.com/terms-trademark-policy for details.
32
33 -------------------------------------------------------------------------
34 Contributing author and copyright for this file:
35 This file is from LAMMPS
36 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
37 http://lammps.sandia.gov, Sandia National Laboratories
38 Steve Plimpton, sjplimp@sandia.gov
39
40 Copyright (2003) Sandia Corporation. Under the terms of Contract
41 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
42 certain rights in this software. This software is distributed under
43 the GNU General Public License.
44 ------------------------------------------------------------------------- */
45
46 #include <cmath>
47 #include "angle.h"
48 #include "atom.h"
49 #include "comm.h"
50 #include "force.h"
51 #include "math_const.h"
52 #include "suffix.h"
53 #include "atom_masks.h"
54 #include "memory.h"
55 #include "error.h"
56
57 using namespace LAMMPS_NS;
58 using namespace MathConst;
59
60 /* ---------------------------------------------------------------------- */
61
Angle(LAMMPS * lmp)62 Angle::Angle(LAMMPS *lmp) : Pointers(lmp)
63 {
64 energy = 0.0;
65 writedata = 1;
66
67 allocated = 0;
68 suffix_flag = Suffix::NONE;
69
70 maxeatom = maxvatom = 0;
71 eatom = NULL;
72 vatom = NULL;
73 setflag = NULL;
74
75 datamask = ALL_MASK;
76 datamask_ext = ALL_MASK;
77 }
78
79 /* ---------------------------------------------------------------------- */
80
~Angle()81 Angle::~Angle()
82 {
83 memory->destroy(eatom);
84 memory->destroy(vatom);
85 }
86
87 /* ----------------------------------------------------------------------
88 check if all coeffs are set
89 ------------------------------------------------------------------------- */
90
init()91 void Angle::init()
92 {
93 if (!allocated && atom->nangletypes)
94 error->all(FLERR,"Angle coeffs are not set");
95 for (int i = 1; i <= atom->nangletypes; i++)
96 if (setflag[i] == 0) error->all(FLERR,"All angle coeffs are not set");
97
98 init_style();
99 }
100
101 /* ----------------------------------------------------------------------
102 setup for energy, virial computation
103 see integrate::ev_set() for values of eflag (0-3) and vflag (0-6)
104 ------------------------------------------------------------------------- */
105
ev_setup(int eflag,int vflag)106 void Angle::ev_setup(int eflag, int vflag)
107 {
108 int i,n;
109
110 evflag = 1;
111
112 eflag_either = eflag;
113 eflag_global = eflag % 2;
114 eflag_atom = eflag / 2;
115
116 vflag_either = vflag;
117 vflag_global = vflag % 4;
118 vflag_atom = vflag / 4;
119
120 // reallocate per-atom arrays if necessary
121
122 if (eflag_atom && atom->nmax > maxeatom) {
123 maxeatom = atom->nmax;
124 memory->destroy(eatom);
125 memory->create(eatom,comm->nthreads*maxeatom,"angle:eatom");
126 }
127 if (vflag_atom && atom->nmax > maxvatom) {
128 maxvatom = atom->nmax;
129 memory->destroy(vatom);
130 memory->create(vatom,comm->nthreads*maxvatom,6,"angle:vatom");
131 }
132
133 // zero accumulators
134
135 if (eflag_global) energy = 0.0;
136 if (vflag_global) for (i = 0; i < 6; i++) virial[i] = 0.0;
137 if (eflag_atom) {
138 n = atom->nlocal;
139 if (force->newton_bond) n += atom->nghost;
140 for (i = 0; i < n; i++) eatom[i] = 0.0;
141 }
142 if (vflag_atom) {
143 n = atom->nlocal;
144 if (force->newton_bond) n += atom->nghost;
145 for (i = 0; i < n; i++) {
146 vatom[i][0] = 0.0;
147 vatom[i][1] = 0.0;
148 vatom[i][2] = 0.0;
149 vatom[i][3] = 0.0;
150 vatom[i][4] = 0.0;
151 vatom[i][5] = 0.0;
152 }
153 }
154 }
155
156 /* ----------------------------------------------------------------------
157 tally energy and virial into global and per-atom accumulators
158 virial = r1F1 + r2F2 + r3F3 = (r1-r2) F1 + (r3-r2) F3 = del1*f1 + del2*f3
159 ------------------------------------------------------------------------- */
160
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)161 void Angle::ev_tally(int i, int j, int k, int nlocal, int newton_bond,
162 double eangle, double *f1, double *f3,
163 double delx1, double dely1, double delz1,
164 double delx2, double dely2, double delz2)
165 {
166 double eanglethird,v[6];
167
168 if (eflag_either) {
169 if (eflag_global) {
170 if (newton_bond) energy += eangle;
171 else {
172 eanglethird = THIRD*eangle;
173 if (i < nlocal) energy += eanglethird;
174 if (j < nlocal) energy += eanglethird;
175 if (k < nlocal) energy += eanglethird;
176 }
177 }
178 if (eflag_atom) {
179 eanglethird = THIRD*eangle;
180 if (newton_bond || i < nlocal) eatom[i] += eanglethird;
181 if (newton_bond || j < nlocal) eatom[j] += eanglethird;
182 if (newton_bond || k < nlocal) eatom[k] += eanglethird;
183 }
184 }
185
186 if (vflag_either) {
187 v[0] = delx1*f1[0] + delx2*f3[0];
188 v[1] = dely1*f1[1] + dely2*f3[1];
189 v[2] = delz1*f1[2] + delz2*f3[2];
190 v[3] = delx1*f1[1] + delx2*f3[1];
191 v[4] = delx1*f1[2] + delx2*f3[2];
192 v[5] = dely1*f1[2] + dely2*f3[2];
193
194 if (vflag_global) {
195 if (newton_bond) {
196 virial[0] += v[0];
197 virial[1] += v[1];
198 virial[2] += v[2];
199 virial[3] += v[3];
200 virial[4] += v[4];
201 virial[5] += v[5];
202 } else {
203 if (i < nlocal) {
204 virial[0] += THIRD*v[0];
205 virial[1] += THIRD*v[1];
206 virial[2] += THIRD*v[2];
207 virial[3] += THIRD*v[3];
208 virial[4] += THIRD*v[4];
209 virial[5] += THIRD*v[5];
210 }
211 if (j < nlocal) {
212 virial[0] += THIRD*v[0];
213 virial[1] += THIRD*v[1];
214 virial[2] += THIRD*v[2];
215 virial[3] += THIRD*v[3];
216 virial[4] += THIRD*v[4];
217 virial[5] += THIRD*v[5];
218 }
219 if (k < nlocal) {
220 virial[0] += THIRD*v[0];
221 virial[1] += THIRD*v[1];
222 virial[2] += THIRD*v[2];
223 virial[3] += THIRD*v[3];
224 virial[4] += THIRD*v[4];
225 virial[5] += THIRD*v[5];
226 }
227 }
228 }
229
230 if (vflag_atom) {
231 if (newton_bond || i < nlocal) {
232 vatom[i][0] += THIRD*v[0];
233 vatom[i][1] += THIRD*v[1];
234 vatom[i][2] += THIRD*v[2];
235 vatom[i][3] += THIRD*v[3];
236 vatom[i][4] += THIRD*v[4];
237 vatom[i][5] += THIRD*v[5];
238 }
239 if (newton_bond || j < nlocal) {
240 vatom[j][0] += THIRD*v[0];
241 vatom[j][1] += THIRD*v[1];
242 vatom[j][2] += THIRD*v[2];
243 vatom[j][3] += THIRD*v[3];
244 vatom[j][4] += THIRD*v[4];
245 vatom[j][5] += THIRD*v[5];
246 }
247 if (newton_bond || k < nlocal) {
248 vatom[k][0] += THIRD*v[0];
249 vatom[k][1] += THIRD*v[1];
250 vatom[k][2] += THIRD*v[2];
251 vatom[k][3] += THIRD*v[3];
252 vatom[k][4] += THIRD*v[4];
253 vatom[k][5] += THIRD*v[5];
254 }
255 }
256 }
257 }
258
259 /* ---------------------------------------------------------------------- */
260
memory_usage()261 double Angle::memory_usage()
262 {
263 double bytes = comm->nthreads*maxeatom * sizeof(double);
264 bytes += comm->nthreads*maxvatom*6 * sizeof(double);
265 return bytes;
266 }
267