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 /* ----------------------------------------------------------------------
16 Contributing author: Axel Kohlmeyer (Temple U)
17 ------------------------------------------------------------------------- */
18
19 #include "improper_cossq_omp.h"
20
21 #include "atom.h"
22 #include "comm.h"
23 #include "force.h"
24 #include "neighbor.h"
25
26 #include <cmath>
27
28 #include "omp_compat.h"
29 #include "suffix.h"
30 using namespace LAMMPS_NS;
31
32 #define TOLERANCE 0.05
33 #define SMALL 0.001
34
35 /* ---------------------------------------------------------------------- */
36
ImproperCossqOMP(class LAMMPS * lmp)37 ImproperCossqOMP::ImproperCossqOMP(class LAMMPS *lmp)
38 : ImproperCossq(lmp), ThrOMP(lmp,THR_IMPROPER)
39 {
40 suffix_flag |= Suffix::OMP;
41 }
42
43 /* ---------------------------------------------------------------------- */
44
compute(int eflag,int vflag)45 void ImproperCossqOMP::compute(int eflag, int vflag)
46 {
47 ev_init(eflag,vflag);
48
49 const int nall = atom->nlocal + atom->nghost;
50 const int nthreads = comm->nthreads;
51 const int inum = neighbor->nimproperlist;
52
53 #if defined(_OPENMP)
54 #pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(eflag,vflag)
55 #endif
56 {
57 int ifrom, ito, tid;
58
59 loop_setup_thr(ifrom, ito, tid, inum, nthreads);
60 ThrData *thr = fix->get_thr(tid);
61 thr->timer(Timer::START);
62 ev_setup_thr(eflag, vflag, nall, eatom, vatom, cvatom, thr);
63
64 if (inum > 0) {
65 if (evflag) {
66 if (eflag) {
67 if (force->newton_bond) eval<1,1,1>(ifrom, ito, thr);
68 else eval<1,1,0>(ifrom, ito, thr);
69 } else {
70 if (force->newton_bond) eval<1,0,1>(ifrom, ito, thr);
71 else eval<1,0,0>(ifrom, ito, thr);
72 }
73 } else {
74 if (force->newton_bond) eval<0,0,1>(ifrom, ito, thr);
75 else eval<0,0,0>(ifrom, ito, thr);
76 }
77 }
78 thr->timer(Timer::BOND);
79 reduce_thr(this, eflag, vflag, thr);
80 } // end of omp parallel region
81 }
82
83 template <int EVFLAG, int EFLAG, int NEWTON_BOND>
eval(int nfrom,int nto,ThrData * const thr)84 void ImproperCossqOMP::eval(int nfrom, int nto, ThrData * const thr)
85 {
86 int i1,i2,i3,i4,n,type;
87 double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z;
88 double eimproper,f1[3],f2[3],f3[3],f4[3];
89 double rjisq, rji, rlksq, rlk, cosphi, angfac;
90 double cjiji, clkji, clklk, cfact1, cfact2, cfact3;
91
92 eimproper = 0.0;
93
94 const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0];
95 dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0];
96 const int5_t * _noalias const improperlist = (int5_t *) neighbor->improperlist[0];
97 const int nlocal = atom->nlocal;
98
99 for (n = nfrom; n < nto; n++) {
100 i1 = improperlist[n].a;
101 i2 = improperlist[n].b;
102 i3 = improperlist[n].c;
103 i4 = improperlist[n].d;
104 type = improperlist[n].t;
105
106 /* separation vector between i1 and i2, (i2-i1) */
107 vb1x = x[i2].x - x[i1].x;
108 vb1y = x[i2].y - x[i1].y;
109 vb1z = x[i2].z - x[i1].z;
110 rjisq = vb1x*vb1x + vb1y*vb1y + vb1z*vb1z ;
111 rji = sqrt(rjisq);
112
113 /* separation vector between i2 and i3 (i3-i2) */
114 vb2x = x[i3].x - x[i2].x;
115 vb2y = x[i3].y - x[i2].y;
116 vb2z = x[i3].z - x[i2].z;
117
118 /* separation vector between i3 and i4, (i4-i3) */
119 vb3x = x[i4].x - x[i3].x;
120 vb3y = x[i4].y - x[i3].y;
121 vb3z = x[i4].z - x[i3].z;
122 rlksq = vb3x*vb3x + vb3y*vb3y + vb3z*vb3z ;
123 rlk = sqrt(rlksq);
124
125 cosphi = (vb3x*vb1x + vb3y*vb1y + vb3z*vb1z)/(rji * rlk);
126
127 /* Check that cos(phi) is in the correct limits. */
128 if (cosphi > 1.0 + TOLERANCE || cosphi < (-1.0 - TOLERANCE))
129 problem(FLERR, i1, i2, i3, i4);
130
131 /* Apply corrections to round-off errors. */
132 if (cosphi > 1.0) cosphi -= SMALL;
133 if (cosphi < -1.0) cosphi += SMALL;
134
135 /* Calculate the angle: */
136 double torangle = acos(cosphi);
137 cosphi = cos(torangle - chi[type]);
138
139 if (EFLAG) eimproper = 0.5 * k[type] * cosphi * cosphi;
140
141 /*
142 printf("The tags: %d-%d-%d-%d, of type %d .\n",atom->tag[i1],atom->tag[i2],atom->tag[i3],atom->tag[i4],type);
143 printf("The ji vector: %f, %f, %f.\nThe lk vector: %f, %f, %f.\n", vb1x,vb1y,vb1z,vb3x,vb3y,vb3z);
144 printf("The cosine of the angle: %-1.16e.\n", cosphi);
145 printf("The energy of the improper: %-1.16e with prefactor %-1.16e.\n", eimproper, 0.5*k[type]);
146 */
147
148 /* Work out forces. */
149 angfac = - k[type] * cosphi;
150
151 cjiji = rjisq;
152 clklk = rlksq;
153 /*CLKJI = RXLK * RXJI + RYLK * RYJI + RZLK * RZJI */
154 clkji = vb3x*vb1x + vb3y*vb1y + vb3z*vb1z;
155
156 /*CFACT1 = CLKLK * CJIJI
157 CFACT1 = SQRT(CFACT1)
158 CFACT1 = ANGFAC / CFACT1*/
159 cfact1 = angfac / sqrt(clklk * cjiji);
160 /*CFACT2 = CLKJI / CLKLK*/
161 cfact2 = clkji / clklk;
162 /*CFACT3 = CLKJI / CJIJI*/
163 cfact3 = clkji / cjiji;
164
165 /*FIX = -RXLK + CFACT3 * RXJI
166 FIY = -RYLK + CFACT3 * RYJI
167 FIZ = -RZLK + CFACT3 * RZJI*/
168 f1[0] = - vb3x + cfact3 * vb1x;
169 f1[1] = - vb3y + cfact3 * vb1y;
170 f1[2] = - vb3z + cfact3 * vb1z;
171
172 /*FJX = -FIX
173 FJY = -FIY
174 FJZ = -FIZ*/
175 f2[0] = - f1[0];
176 f2[1] = - f1[1];
177 f2[2] = - f1[2];
178
179 /*FKX = CFACT2 * RXLK - RXJI
180 FKY = CFACT2 * RYLK - RYJI
181 FKZ = CFACT2 * RZLK - RZJI*/
182 f3[0] = cfact2 * vb3x - vb1x;
183 f3[1] = cfact2 * vb3y - vb1y;
184 f3[2] = cfact2 * vb3z - vb1z;
185
186 /*FLX = -FKX
187 FLY = -FKY
188 FLZ = -FKZ*/
189 f4[0] = - f3[0];
190 f4[1] = - f3[1];
191 f4[2] = - f3[2];
192
193 /*FIX = FIX * CFACT1
194 FIY = FIY * CFACT1
195 FIZ = FIZ * CFACT1*/
196 f1[0] *= cfact1;
197 f1[1] *= cfact1;
198 f1[2] *= cfact1;
199
200 /*FJX = FJX * CFACT1
201 FJY = FJY * CFACT1
202 FJZ = FJZ * CFACT1*/
203 f2[0] *= cfact1;
204 f2[1] *= cfact1;
205 f2[2] *= cfact1;
206
207 /*FKX = FKX * CFACT1
208 FKY = FKY * CFACT1
209 FKZ = FKZ * CFACT1*/
210 f3[0] *= cfact1;
211 f3[1] *= cfact1;
212 f3[2] *= cfact1;
213
214 /*FLX = FLX * CFACT1
215 FLY = FLY * CFACT1
216 FLZ = FLZ * CFACT1*/
217 f4[0] *= cfact1;
218 f4[1] *= cfact1;
219 f4[2] *= cfact1;
220
221 /* Apply force to each of 4 atoms */
222 if (NEWTON_BOND || i1 < nlocal) {
223 f[i1].x += f1[0];
224 f[i1].y += f1[1];
225 f[i1].z += f1[2];
226 }
227
228 if (NEWTON_BOND || i2 < nlocal) {
229 f[i2].x += f2[0];
230 f[i2].y += f2[1];
231 f[i2].z += f2[2];
232 }
233
234 if (NEWTON_BOND || i3 < nlocal) {
235 f[i3].x += f3[0];
236 f[i3].y += f3[1];
237 f[i3].z += f3[2];
238 }
239
240 if (NEWTON_BOND || i4 < nlocal) {
241 f[i4].x += f4[0];
242 f[i4].y += f4[1];
243 f[i4].z += f4[2];
244 }
245
246 if (EVFLAG)
247 ev_tally_thr(this,i1,i2,i3,i4,nlocal,NEWTON_BOND,eimproper,f1,f3,f4,
248 -vb1x,-vb1y,-vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,thr);
249 }
250 }
251