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 "dihedral.h"
16 #include "atom.h"
17 #include "comm.h"
18 #include "force.h"
19 #include "atom_masks.h"
20 #include "memory.h"
21 #include "error.h"
22 #include "suffix.h"
23 #include "update.h"
24
25 using namespace LAMMPS_NS;
26
27 /* ----------------------------------------------------------------------
28 set dihedral contribution to Vdwl and Coulombic energy to 0.0
29 DihedralCharmm will override this
30 ------------------------------------------------------------------------- */
31
Dihedral(LAMMPS * lmp)32 Dihedral::Dihedral(LAMMPS *lmp) : Pointers(lmp)
33 {
34 energy = 0.0;
35 writedata = 0;
36
37 allocated = 0;
38 suffix_flag = Suffix::NONE;
39
40 maxeatom = maxvatom = maxcvatom = 0;
41 eatom = nullptr;
42 vatom = nullptr;
43 cvatom = nullptr;
44 setflag = nullptr;
45 centroidstressflag = CENTROID_AVAIL;
46
47 execution_space = Host;
48 datamask_read = ALL_MASK;
49 datamask_modify = ALL_MASK;
50
51 copymode = 0;
52 }
53
54 /* ---------------------------------------------------------------------- */
55
~Dihedral()56 Dihedral::~Dihedral()
57 {
58 if (copymode) return;
59
60 memory->destroy(eatom);
61 memory->destroy(vatom);
62 memory->destroy(cvatom);
63 }
64
65 /* ----------------------------------------------------------------------
66 check if all coeffs are set
67 ------------------------------------------------------------------------- */
68
init()69 void Dihedral::init()
70 {
71 if (!allocated && atom->ndihedraltypes)
72 error->all(FLERR,"Dihedral coeffs are not set");
73 for (int i = 1; i <= atom->ndihedraltypes; i++)
74 if (setflag[i] == 0) error->all(FLERR,"All dihedral coeffs are not set");
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 Dihedral::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,"dihedral: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,"dihedral: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,"dihedral: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 + r4F4 = (r1-r2) F1 + (r3-r2) F3 + (r4-r2) F4
180 = (r1-r2) F1 + (r3-r2) F3 + (r4-r3 + r3-r2) F4
181 = vb1*f1 + vb2*f3 + (vb3+vb2)*f4
182 ------------------------------------------------------------------------- */
183
ev_tally(int i1,int i2,int i3,int i4,int nlocal,int newton_bond,double edihedral,double * f1,double * f3,double * f4,double vb1x,double vb1y,double vb1z,double vb2x,double vb2y,double vb2z,double vb3x,double vb3y,double vb3z)184 void Dihedral::ev_tally(int i1, int i2, int i3, int i4,
185 int nlocal, int newton_bond,
186 double edihedral, double *f1, double *f3, double *f4,
187 double vb1x, double vb1y, double vb1z,
188 double vb2x, double vb2y, double vb2z,
189 double vb3x, double vb3y, double vb3z)
190 {
191 double edihedralquarter,v[6];
192
193 if (eflag_either) {
194 if (eflag_global) {
195 if (newton_bond) energy += edihedral;
196 else {
197 edihedralquarter = 0.25*edihedral;
198 if (i1 < nlocal) energy += edihedralquarter;
199 if (i2 < nlocal) energy += edihedralquarter;
200 if (i3 < nlocal) energy += edihedralquarter;
201 if (i4 < nlocal) energy += edihedralquarter;
202 }
203 }
204 if (eflag_atom) {
205 edihedralquarter = 0.25*edihedral;
206 if (newton_bond || i1 < nlocal) eatom[i1] += edihedralquarter;
207 if (newton_bond || i2 < nlocal) eatom[i2] += edihedralquarter;
208 if (newton_bond || i3 < nlocal) eatom[i3] += edihedralquarter;
209 if (newton_bond || i4 < nlocal) eatom[i4] += edihedralquarter;
210 }
211 }
212
213 if (vflag_either) {
214 v[0] = vb1x*f1[0] + vb2x*f3[0] + (vb3x+vb2x)*f4[0];
215 v[1] = vb1y*f1[1] + vb2y*f3[1] + (vb3y+vb2y)*f4[1];
216 v[2] = vb1z*f1[2] + vb2z*f3[2] + (vb3z+vb2z)*f4[2];
217 v[3] = vb1x*f1[1] + vb2x*f3[1] + (vb3x+vb2x)*f4[1];
218 v[4] = vb1x*f1[2] + vb2x*f3[2] + (vb3x+vb2x)*f4[2];
219 v[5] = vb1y*f1[2] + vb2y*f3[2] + (vb3y+vb2y)*f4[2];
220
221 if (vflag_global) {
222 if (newton_bond) {
223 virial[0] += v[0];
224 virial[1] += v[1];
225 virial[2] += v[2];
226 virial[3] += v[3];
227 virial[4] += v[4];
228 virial[5] += v[5];
229 } else {
230 if (i1 < nlocal) {
231 virial[0] += 0.25*v[0];
232 virial[1] += 0.25*v[1];
233 virial[2] += 0.25*v[2];
234 virial[3] += 0.25*v[3];
235 virial[4] += 0.25*v[4];
236 virial[5] += 0.25*v[5];
237 }
238 if (i2 < nlocal) {
239 virial[0] += 0.25*v[0];
240 virial[1] += 0.25*v[1];
241 virial[2] += 0.25*v[2];
242 virial[3] += 0.25*v[3];
243 virial[4] += 0.25*v[4];
244 virial[5] += 0.25*v[5];
245 }
246 if (i3 < nlocal) {
247 virial[0] += 0.25*v[0];
248 virial[1] += 0.25*v[1];
249 virial[2] += 0.25*v[2];
250 virial[3] += 0.25*v[3];
251 virial[4] += 0.25*v[4];
252 virial[5] += 0.25*v[5];
253 }
254 if (i4 < nlocal) {
255 virial[0] += 0.25*v[0];
256 virial[1] += 0.25*v[1];
257 virial[2] += 0.25*v[2];
258 virial[3] += 0.25*v[3];
259 virial[4] += 0.25*v[4];
260 virial[5] += 0.25*v[5];
261 }
262 }
263 }
264
265 if (vflag_atom) {
266 if (newton_bond || i1 < nlocal) {
267 vatom[i1][0] += 0.25*v[0];
268 vatom[i1][1] += 0.25*v[1];
269 vatom[i1][2] += 0.25*v[2];
270 vatom[i1][3] += 0.25*v[3];
271 vatom[i1][4] += 0.25*v[4];
272 vatom[i1][5] += 0.25*v[5];
273 }
274 if (newton_bond || i2 < nlocal) {
275 vatom[i2][0] += 0.25*v[0];
276 vatom[i2][1] += 0.25*v[1];
277 vatom[i2][2] += 0.25*v[2];
278 vatom[i2][3] += 0.25*v[3];
279 vatom[i2][4] += 0.25*v[4];
280 vatom[i2][5] += 0.25*v[5];
281 }
282 if (newton_bond || i3 < nlocal) {
283 vatom[i3][0] += 0.25*v[0];
284 vatom[i3][1] += 0.25*v[1];
285 vatom[i3][2] += 0.25*v[2];
286 vatom[i3][3] += 0.25*v[3];
287 vatom[i3][4] += 0.25*v[4];
288 vatom[i3][5] += 0.25*v[5];
289 }
290 if (newton_bond || i4 < nlocal) {
291 vatom[i4][0] += 0.25*v[0];
292 vatom[i4][1] += 0.25*v[1];
293 vatom[i4][2] += 0.25*v[2];
294 vatom[i4][3] += 0.25*v[3];
295 vatom[i4][4] += 0.25*v[4];
296 vatom[i4][5] += 0.25*v[5];
297 }
298 }
299 }
300
301 // per-atom centroid virial
302 if (cvflag_atom) {
303
304 // r0 = (r1+r2+r3+r4)/4
305 // rij = ri-rj
306 // total virial = r10*f1 + r20*f2 + r30*f3 + r40*f4
307 // vb1: r12
308 // vb2: r32
309 // vb3: r43
310
311 if (newton_bond || i1 < nlocal) {
312 double a1[3];
313
314 // a1 = r10 = (3*r12 - 2*r32 - r43)/4
315 a1[0] = 0.25*(3*vb1x - 2*vb2x - vb3x);
316 a1[1] = 0.25*(3*vb1y - 2*vb2y - vb3y);
317 a1[2] = 0.25*(3*vb1z - 2*vb2z - vb3z);
318
319 cvatom[i1][0] += a1[0]*f1[0];
320 cvatom[i1][1] += a1[1]*f1[1];
321 cvatom[i1][2] += a1[2]*f1[2];
322 cvatom[i1][3] += a1[0]*f1[1];
323 cvatom[i1][4] += a1[0]*f1[2];
324 cvatom[i1][5] += a1[1]*f1[2];
325 cvatom[i1][6] += a1[1]*f1[0];
326 cvatom[i1][7] += a1[2]*f1[0];
327 cvatom[i1][8] += a1[2]*f1[1];
328 }
329 if (newton_bond || i2 < nlocal) {
330 double a2[3];
331 double f2[3];
332
333 // a2 = r20 = ( -r12 - 2*r32 - r43)/4
334 a2[0] = 0.25*(-vb1x - 2*vb2x - vb3x);
335 a2[1] = 0.25*(-vb1y - 2*vb2y - vb3y);
336 a2[2] = 0.25*(-vb1z - 2*vb2z - vb3z);
337
338 f2[0] = - f1[0] - f3[0] - f4[0];
339 f2[1] = - f1[1] - f3[1] - f4[1];
340 f2[2] = - f1[2] - f3[2] - f4[2];
341
342 cvatom[i2][0] += a2[0]*f2[0];
343 cvatom[i2][1] += a2[1]*f2[1];
344 cvatom[i2][2] += a2[2]*f2[2];
345 cvatom[i2][3] += a2[0]*f2[1];
346 cvatom[i2][4] += a2[0]*f2[2];
347 cvatom[i2][5] += a2[1]*f2[2];
348 cvatom[i2][6] += a2[1]*f2[0];
349 cvatom[i2][7] += a2[2]*f2[0];
350 cvatom[i2][8] += a2[2]*f2[1];
351 }
352 if (newton_bond || i3 < nlocal) {
353 double a3[3];
354
355 // a3 = r30 = ( -r12 + 2*r32 - r43)/4
356 a3[0] = 0.25*(-vb1x + 2*vb2x - vb3x);
357 a3[1] = 0.25*(-vb1y + 2*vb2y - vb3y);
358 a3[2] = 0.25*(-vb1z + 2*vb2z - vb3z);
359
360 cvatom[i3][0] += a3[0]*f3[0];
361 cvatom[i3][1] += a3[1]*f3[1];
362 cvatom[i3][2] += a3[2]*f3[2];
363 cvatom[i3][3] += a3[0]*f3[1];
364 cvatom[i3][4] += a3[0]*f3[2];
365 cvatom[i3][5] += a3[1]*f3[2];
366 cvatom[i3][6] += a3[1]*f3[0];
367 cvatom[i3][7] += a3[2]*f3[0];
368 cvatom[i3][8] += a3[2]*f3[1];
369 }
370 if (newton_bond || i4 < nlocal) {
371 double a4[3];
372
373 // a4 = r40 = ( -r12 + 2*r32 + 3*r43)/4
374 a4[0] = 0.25*(-vb1x + 2*vb2x + 3*vb3x);
375 a4[1] = 0.25*(-vb1y + 2*vb2y + 3*vb3y);
376 a4[2] = 0.25*(-vb1z + 2*vb2z + 3*vb3z);
377
378 cvatom[i4][0] += a4[0]*f4[0];
379 cvatom[i4][1] += a4[1]*f4[1];
380 cvatom[i4][2] += a4[2]*f4[2];
381 cvatom[i4][3] += a4[0]*f4[1];
382 cvatom[i4][4] += a4[0]*f4[2];
383 cvatom[i4][5] += a4[1]*f4[2];
384 cvatom[i4][6] += a4[1]*f4[0];
385 cvatom[i4][7] += a4[2]*f4[0];
386 cvatom[i4][8] += a4[2]*f4[1];
387 }
388 }
389 }
390
391 /* ---------------------------------------------------------------------- */
392
problem(const char * filename,int lineno,int i1,int i2,int i3,int i4)393 void Dihedral::problem(const char *filename, int lineno,
394 int i1, int i2, int i3, int i4)
395 {
396 const auto x = atom->x;
397 auto warn = fmt::format("Dihedral problem: {} {} {} {} {} {}\n",
398 comm->me, update->ntimestep, atom->tag[i1],
399 atom->tag[i2], atom->tag[i3], atom->tag[i4]);
400 warn += fmt::format("WARNING: 1st atom: {} {:.8} {:.8} {:.8}\n",
401 comm->me,x[i1][0],x[i1][1],x[i1][2]);
402 warn += fmt::format("WARNING: 2nd atom: {} {:.8} {:.8} {:.8}\n",
403 comm->me,x[i2][0],x[i2][1],x[i2][2]);
404 warn += fmt::format("WARNING: 3rd atom: {} {:.8} {:.8} {:.8}\n",
405 comm->me,x[i3][0],x[i3][1],x[i3][2]);
406 warn += fmt::format("WARNING: 4th atom: {} {:.8} {:.8} {:.8}",
407 comm->me,x[i4][0],x[i4][1],x[i4][2]);
408 error->warning(filename, lineno, warn);
409 }
410
411 /* ---------------------------------------------------------------------- */
412
memory_usage()413 double Dihedral::memory_usage()
414 {
415 double bytes = (double)comm->nthreads*maxeatom * sizeof(double);
416 bytes += (double)comm->nthreads*maxvatom*6 * sizeof(double);
417 bytes += (double)comm->nthreads*maxcvatom*9 * sizeof(double);
418 return bytes;
419 }
420