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