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_fourier_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
ImproperFourierOMP(class LAMMPS * lmp)37 ImproperFourierOMP::ImproperFourierOMP(class LAMMPS *lmp)
38 : ImproperFourier(lmp), ThrOMP(lmp,THR_IMPROPER)
39 {
40 suffix_flag |= Suffix::OMP;
41 }
42
43 /* ---------------------------------------------------------------------- */
44
compute(int eflag,int vflag)45 void ImproperFourierOMP::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 ImproperFourierOMP::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
89 const double * const * const x = atom->x;
90 const int * const * const improperlist = neighbor->improperlist;
91
92 for (n = nfrom; n < nto; n++) {
93 i1 = improperlist[n][0];
94 i2 = improperlist[n][1];
95 i3 = improperlist[n][2];
96 i4 = improperlist[n][3];
97 type = improperlist[n][4];
98
99 // 1st bond
100
101 vb1x = x[i2][0] - x[i1][0];
102 vb1y = x[i2][1] - x[i1][1];
103 vb1z = x[i2][2] - x[i1][2];
104
105 // 2nd bond
106
107 vb2x = x[i3][0] - x[i1][0];
108 vb2y = x[i3][1] - x[i1][1];
109 vb2z = x[i3][2] - x[i1][2];
110
111 // 3rd bond
112
113 vb3x = x[i4][0] - x[i1][0];
114 vb3y = x[i4][1] - x[i1][1];
115 vb3z = x[i4][2] - x[i1][2];
116
117 add1_thr<EVFLAG,EFLAG,NEWTON_BOND>(i1,i2,i3,i4,type,
118 vb1x,vb1y,vb1z,
119 vb2x,vb2y,vb2z,
120 vb3x,vb3y,vb3z,thr);
121 if (all[type]) {
122 add1_thr<EVFLAG,EFLAG,NEWTON_BOND>(i1,i4,i2,i3,type,
123 vb3x,vb3y,vb3z,
124 vb1x,vb1y,vb1z,
125 vb2x,vb2y,vb2z,thr);
126 add1_thr<EVFLAG,EFLAG,NEWTON_BOND>(i1,i3,i4,i2,type,
127 vb2x,vb2y,vb2z,
128 vb3x,vb3y,vb3z,
129 vb1x,vb1y,vb1z,thr);
130 }
131 }
132 }
133
134 template <int EVFLAG, int EFLAG, int NEWTON_BOND>
add1_thr(const int i1,const int i2,const int i3,const int i4,const int type,const double & vb1x,const double & vb1y,const double & vb1z,const double & vb2x,const double & vb2y,const double & vb2z,const double & vb3x,const double & vb3y,const double & vb3z,ThrData * const thr)135 void ImproperFourierOMP::add1_thr(const int i1,const int i2,
136 const int i3,const int i4,
137 const int type,
138 const double &vb1x,
139 const double &vb1y,
140 const double &vb1z,
141 const double &vb2x,
142 const double &vb2y,
143 const double &vb2z,
144 const double &vb3x,
145 const double &vb3y,
146 const double &vb3z,
147 ThrData * const thr)
148 {
149 double eimproper,f1[3],f2[3],f3[3],f4[3];
150 double c,c2,a,s,projhfg,dhax,dhay,dhaz,dahx,dahy,dahz,cotphi;
151 double ax,ay,az,ra2,rh2,ra,rh,rar,rhr,arx,ary,arz,hrx,hry,hrz;
152
153 double * const * const f = thr->get_f();
154 const int nlocal = atom->nlocal;
155
156 eimproper = 0.0;
157
158 // c0 calculation
159 // A = vb1 X vb2 is perpendicular to IJK plane
160
161 ax = vb1y*vb2z-vb1z*vb2y;
162 ay = vb1z*vb2x-vb1x*vb2z;
163 az = vb1x*vb2y-vb1y*vb2x;
164 ra2 = ax*ax+ay*ay+az*az;
165 rh2 = vb3x*vb3x+vb3y*vb3y+vb3z*vb3z;
166 ra = sqrt(ra2);
167 rh = sqrt(rh2);
168 if (ra < SMALL) ra = SMALL;
169 if (rh < SMALL) rh = SMALL;
170
171 rar = 1/ra;
172 rhr = 1/rh;
173 arx = ax*rar;
174 ary = ay*rar;
175 arz = az*rar;
176 hrx = vb3x*rhr;
177 hry = vb3y*rhr;
178 hrz = vb3z*rhr;
179
180 c = arx*hrx+ary*hry+arz*hrz;
181
182 // error check
183
184 if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE))
185 problem(FLERR, i1, i2, i3, i4);
186
187 if (c > 1.0) c = 1.0;
188 if (c < -1.0) c = -1.0;
189
190 s = sqrt(1.0 - c*c);
191 if (s < SMALL) s = SMALL;
192 cotphi = c/s;
193
194 projhfg = (vb3x*vb1x+vb3y*vb1y+vb3z*vb1z) /
195 sqrt(vb1x*vb1x+vb1y*vb1y+vb1z*vb1z);
196 projhfg += (vb3x*vb2x+vb3y*vb2y+vb3z*vb2z) /
197 sqrt(vb2x*vb2x+vb2y*vb2y+vb2z*vb2z);
198 if (projhfg > 0.0) {
199 s *= -1.0;
200 cotphi *= -1.0;
201 }
202
203 // force and energy
204 // E = k ( C0 + C1 cos(w) + C2 cos(w)
205
206 c2 = 2.0*s*s-1.0;
207 if (EFLAG) eimproper = k[type]*(C0[type]+C1[type]*s+C2[type]*c2);
208
209 // dhax = diffrence between H and A in X direction, etc
210
211 a = k[type]*(C1[type]+4.0*C2[type]*s)*cotphi;
212 dhax = hrx-c*arx;
213 dhay = hry-c*ary;
214 dhaz = hrz-c*arz;
215
216 dahx = arx-c*hrx;
217 dahy = ary-c*hry;
218 dahz = arz-c*hrz;
219
220 f2[0] = (dhay*vb1z - dhaz*vb1y)*rar*a;
221 f2[1] = (dhaz*vb1x - dhax*vb1z)*rar*a;
222 f2[2] = (dhax*vb1y - dhay*vb1x)*rar*a;
223
224 f3[0] = (-dhay*vb2z + dhaz*vb2y)*rar*a;
225 f3[1] = (-dhaz*vb2x + dhax*vb2z)*rar*a;
226 f3[2] = (-dhax*vb2y + dhay*vb2x)*rar*a;
227
228 f4[0] = dahx*rhr*a;
229 f4[1] = dahy*rhr*a;
230 f4[2] = dahz*rhr*a;
231
232 f1[0] = -(f2[0] + f3[0] + f4[0]);
233 f1[1] = -(f2[1] + f3[1] + f4[1]);
234 f1[2] = -(f2[2] + f3[2] + f4[2]);
235
236 // apply force to each of 4 atoms
237
238 if (NEWTON_BOND || i1 < nlocal) {
239 f[i1][0] += f1[0];
240 f[i1][1] += f1[1];
241 f[i1][2] += f1[2];
242 }
243
244 if (NEWTON_BOND || i2 < nlocal) {
245 f[i2][0] += f3[0];
246 f[i2][1] += f3[1];
247 f[i2][2] += f3[2];
248 }
249
250 if (NEWTON_BOND || i3 < nlocal) {
251 f[i3][0] += f2[0];
252 f[i3][1] += f2[1];
253 f[i3][2] += f2[2];
254 }
255
256 if (NEWTON_BOND || i4 < nlocal) {
257 f[i4][0] += f4[0];
258 f[i4][1] += f4[1];
259 f[i4][2] += f4[2];
260 }
261
262 if (EVFLAG)
263 ev_tally_thr(this,i1,i2,i3,i4,nlocal,NEWTON_BOND,eimproper,f1,f2,f4,
264 -vb1x,-vb1y,-vb1z,vb2x-vb1x,vb2y-vb1y,vb2z-vb1z,vb3x-vb2x,vb3y-vb2y,vb3z-vb2z,thr);
265
266 }
267