1 // **************************************************************************
2 // gayberne.cu
3 // -------------------
4 // W. Michael Brown (ORNL)
5 //
6 // Device code for Gay-Berne potential acceleration
7 //
8 // __________________________________________________________________________
9 // This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
10 // __________________________________________________________________________
11 //
12 // begin :
13 // email : brownw@ornl.gov
14 // ***************************************************************************
15
16 #if defined(NV_KERNEL) || defined(USE_HIP)
17 #include "lal_ellipsoid_extra.h"
18 #endif
19
compute_eta_torque(numtyp m[9],numtyp m2[9],const numtyp4 shape,numtyp ans[9])20 ucl_inline void compute_eta_torque(numtyp m[9],numtyp m2[9], const numtyp4 shape,
21 numtyp ans[9])
22 {
23 numtyp den = m[3]*m[2]*m[7]-m[0]*m[5]*m[7]-
24 m[2]*m[6]*m[4]+m[1]*m[6]*m[5]-
25 m[3]*m[1]*m[8]+m[0]*m[4]*m[8];
26 den = ucl_recip(den);
27
28 ans[0] = shape.x*(m[5]*m[1]*m2[2]+(numtyp)2.0*m[4]*m[8]*m2[0]-
29 m[4]*m2[2]*m[2]-(numtyp)2.0*m[5]*m2[0]*m[7]+
30 m2[1]*m[2]*m[7]-m2[1]*m[1]*m[8]-
31 m[3]*m[8]*m2[1]+m[6]*m[5]*m2[1]+
32 m[3]*m2[2]*m[7]-m2[2]*m[6]*m[4])*den;
33
34 ans[1] = shape.x*(m[2]*m2[0]*m[7]-m[8]*m2[0]*m[1]+
35 (numtyp)2.0*m[0]*m[8]*m2[1]-m[0]*m2[2]*m[5]-
36 (numtyp)2.0*m[6]*m[2]*m2[1]+m2[2]*m[3]*m[2]-
37 m[8]*m[3]*m2[0]+m[6]*m2[0]*m[5]+
38 m[6]*m2[2]*m[1]-m2[2]*m[0]*m[7])*den;
39
40 ans[2] = shape.x*(m[1]*m[5]*m2[0]-m[2]*m2[0]*m[4]-
41 m[0]*m[5]*m2[1]+m[3]*m[2]*m2[1]-
42 m2[1]*m[0]*m[7]-m[6]*m[4]*m2[0]+
43 (numtyp)2.0*m[4]*m[0]*m2[2]-(numtyp)2.0*m[3]*m2[2]*m[1]+
44 m[3]*m[7]*m2[0]+m[6]*m2[1]*m[1])*den;
45
46 ans[3] = shape.y*(-m[4]*m2[5]*m[2]+(numtyp)2.0*m[4]*m[8]*m2[3]+
47 m[5]*m[1]*m2[5]-(numtyp)2.0*m[5]*m2[3]*m[7]+
48 m2[4]*m[2]*m[7]-m2[4]*m[1]*m[8]-
49 m[3]*m[8]*m2[4]+m[6]*m[5]*m2[4]-
50 m2[5]*m[6]*m[4]+m[3]*m2[5]*m[7])*den;
51
52 ans[4] = shape.y*(m[2]*m2[3]*m[7]-m[1]*m[8]*m2[3]+
53 (numtyp)2.0*m[8]*m[0]*m2[4]-m2[5]*m[0]*m[5]-
54 (numtyp)2.0*m[6]*m2[4]*m[2]-m[3]*m[8]*m2[3]+
55 m[6]*m[5]*m2[3]+m[3]*m2[5]*m[2]-
56 m[0]*m2[5]*m[7]+m2[5]*m[1]*m[6])*den;
57
58 ans[5] = shape.y*(m[1]*m[5]*m2[3]-m[2]*m2[3]*m[4]-
59 m[0]*m[5]*m2[4]+m[3]*m[2]*m2[4]+
60 (numtyp)2.0*m[4]*m[0]*m2[5]-m[0]*m2[4]*m[7]+
61 m[1]*m[6]*m2[4]-m2[3]*m[6]*m[4]-
62 (numtyp)2.0*m[3]*m[1]*m2[5]+m[3]*m2[3]*m[7])*den;
63
64 ans[6] = shape.z*(-m[4]*m[2]*m2[8]+m[1]*m[5]*m2[8]+
65 (numtyp)2.0*m[4]*m2[6]*m[8]-m[1]*m2[7]*m[8]+
66 m[2]*m[7]*m2[7]-(numtyp)2.0*m2[6]*m[7]*m[5]-
67 m[3]*m2[7]*m[8]+m[5]*m[6]*m2[7]-
68 m[4]*m[6]*m2[8]+m[7]*m[3]*m2[8])*den;
69
70 ans[7] = shape.z*-(m[1]*m[8]*m2[6]-m[2]*m2[6]*m[7]-
71 (numtyp)2.0*m2[7]*m[0]*m[8]+m[5]*m2[8]*m[0]+
72 (numtyp)2.0*m2[7]*m[2]*m[6]+m[3]*m2[6]*m[8]-
73 m[3]*m[2]*m2[8]-m[5]*m[6]*m2[6]+
74 m[0]*m2[8]*m[7]-m2[8]*m[1]*m[6])*den;
75
76 ans[8] = shape.z*(m[1]*m[5]*m2[6]-m[2]*m2[6]*m[4]-
77 m[0]*m[5]*m2[7]+m[3]*m[2]*m2[7]-
78 m[4]*m[6]*m2[6]-m[7]*m2[7]*m[0]+
79 (numtyp)2.0*m[4]*m2[8]*m[0]+m[7]*m[3]*m2[6]+
80 m[6]*m[1]*m2[7]-(numtyp)2.0*m2[8]*m[3]*m[1])*den;
81 }
82
k_gayberne(const __global numtyp4 * restrict x_,const __global numtyp4 * restrict q,const __global numtyp4 * restrict shape,const __global numtyp4 * restrict well,const __global numtyp * restrict gum,const __global numtyp2 * restrict sig_eps,const int ntypes,const __global numtyp * restrict lshape,const __global int * dev_nbor,const int stride,__global acctyp4 * restrict ans,const int astride,__global acctyp * restrict engv,__global int * restrict err_flag,const int eflag,const int vflag,const int inum,const int t_per_atom)83 __kernel void k_gayberne(const __global numtyp4 *restrict x_,
84 const __global numtyp4 *restrict q,
85 const __global numtyp4 *restrict shape,
86 const __global numtyp4 *restrict well,
87 const __global numtyp *restrict gum,
88 const __global numtyp2 *restrict sig_eps,
89 const int ntypes,
90 const __global numtyp *restrict lshape,
91 const __global int *dev_nbor,
92 const int stride,
93 __global acctyp4 *restrict ans,
94 const int astride,
95 __global acctyp *restrict engv,
96 __global int *restrict err_flag,
97 const int eflag, const int vflag, const int inum,
98 const int t_per_atom) {
99 int tid, ii, offset;
100 atom_info(t_per_atom,ii,tid,offset);
101
102 __local numtyp sp_lj[4];
103 int n_stride;
104 local_allocate_store_ellipse();
105
106 sp_lj[0]=gum[3];
107 sp_lj[1]=gum[4];
108 sp_lj[2]=gum[5];
109 sp_lj[3]=gum[6];
110
111 acctyp4 f, tor;
112 f.x=(acctyp)0; f.y=(acctyp)0; f.z=(acctyp)0;
113 tor.x=(acctyp)0; tor.y=(acctyp)0; tor.z=(acctyp)0;
114 acctyp energy, virial[6];
115 if (EVFLAG) {
116 energy=(acctyp)0;
117 for (int i=0; i<6; i++) virial[i]=(acctyp)0;
118 }
119
120 if (ii<inum) {
121 int nbor, nbor_end;
122 int i, numj;
123 nbor_info_p(dev_nbor,stride,t_per_atom,ii,offset,i,numj,
124 n_stride,nbor_end,nbor);
125
126 numtyp4 ix; fetch4(ix,i,pos_tex);
127 int itype=ix.w;
128 numtyp a1[9], b1[9], g1[9];
129 numtyp4 ishape=shape[itype];
130 {
131 numtyp t[9];
132 gpu_quat_to_mat_trans(q,i,a1);
133 gpu_diag_times3(ishape,a1,t);
134 gpu_transpose_times3(a1,t,g1);
135 gpu_diag_times3(well[itype],a1,t);
136 gpu_transpose_times3(a1,t,b1);
137 }
138
139 numtyp factor_lj;
140 for ( ; nbor<nbor_end; nbor+=n_stride) {
141 int j=dev_nbor[nbor];
142 factor_lj = sp_lj[sbmask(j)];
143 j &= NEIGHMASK;
144
145 numtyp4 jx; fetch4(jx,j,pos_tex);
146 int jtype=jx.w;
147
148 // Compute r12
149 numtyp r12[3];
150 r12[0] = jx.x-ix.x;
151 r12[1] = jx.y-ix.y;
152 r12[2] = jx.z-ix.z;
153 numtyp ir = gpu_dot3(r12,r12);
154
155 ir = ucl_rsqrt(ir);
156 numtyp r = ucl_recip(ir);
157
158 numtyp a2[9];
159 gpu_quat_to_mat_trans(q,j,a2);
160
161 numtyp u_r, dUr[3], tUr[3], eta, teta[3];
162 { // Compute U_r, dUr, eta, and teta
163 // Compute g12
164 numtyp g12[9];
165 {
166 numtyp g2[9];
167 {
168 gpu_diag_times3(shape[jtype],a2,g12);
169 gpu_transpose_times3(a2,g12,g2);
170 gpu_plus3(g1,g2,g12);
171 }
172
173 { // Compute U_r and dUr
174
175 // Compute kappa
176 numtyp kappa[3];
177 gpu_mldivide3(g12,r12,kappa,err_flag);
178
179 // -- replace r12 with r12 hat
180 r12[0]*=ir;
181 r12[1]*=ir;
182 r12[2]*=ir;
183
184 // -- kappa is now / r
185 kappa[0]*=ir;
186 kappa[1]*=ir;
187 kappa[2]*=ir;
188
189 // energy
190
191 // compute u_r and dUr
192 numtyp uslj_rsq;
193 {
194 // Compute distance of closest approach
195 numtyp h12, sigma12;
196 sigma12 = gpu_dot3(r12,kappa);
197 sigma12 = ucl_rsqrt((numtyp)0.5*sigma12);
198 h12 = r-sigma12;
199
200 // -- kappa is now ok
201 kappa[0]*=r;
202 kappa[1]*=r;
203 kappa[2]*=r;
204
205 int mtype=fast_mul(ntypes,itype)+jtype;
206 numtyp sigma = sig_eps[mtype].x;
207 numtyp epsilon = sig_eps[mtype].y;
208 numtyp varrho = sigma/(h12+gum[0]*sigma);
209 numtyp varrho6 = varrho*varrho*varrho;
210 varrho6*=varrho6;
211 numtyp varrho12 = varrho6*varrho6;
212 u_r = (numtyp)4.0*epsilon*(varrho12-varrho6);
213
214 numtyp temp1 = ((numtyp)2.0*varrho12*varrho-varrho6*varrho)/sigma;
215 temp1 = temp1*(numtyp)24.0*epsilon;
216 uslj_rsq = temp1*sigma12*sigma12*sigma12*(numtyp)0.5;
217 numtyp temp2 = gpu_dot3(kappa,r12);
218 uslj_rsq = uslj_rsq*ir*ir;
219
220 dUr[0] = temp1*r12[0]+uslj_rsq*(kappa[0]-temp2*r12[0]);
221 dUr[1] = temp1*r12[1]+uslj_rsq*(kappa[1]-temp2*r12[1]);
222 dUr[2] = temp1*r12[2]+uslj_rsq*(kappa[2]-temp2*r12[2]);
223 }
224
225 // torque for particle 1
226 {
227 numtyp tempv[3], tempv2[3];
228 tempv[0] = -uslj_rsq*kappa[0];
229 tempv[1] = -uslj_rsq*kappa[1];
230 tempv[2] = -uslj_rsq*kappa[2];
231 gpu_row_times3(kappa,g1,tempv2);
232 gpu_cross3(tempv,tempv2,tUr);
233 }
234 }
235 }
236
237 // Compute eta
238 {
239 eta = (numtyp)2.0*lshape[itype]*lshape[jtype];
240 numtyp det_g12 = gpu_det3(g12);
241 eta = ucl_powr(eta/det_g12,gum[1]);
242 }
243
244 // Compute teta
245 numtyp temp[9], tempv[3], tempv2[3];
246 compute_eta_torque(g12,a1,ishape,temp);
247 numtyp temp1 = -eta*gum[1];
248
249 tempv[0] = temp1*temp[0];
250 tempv[1] = temp1*temp[1];
251 tempv[2] = temp1*temp[2];
252 gpu_cross3(a1,tempv,tempv2);
253 teta[0] = tempv2[0];
254 teta[1] = tempv2[1];
255 teta[2] = tempv2[2];
256
257 tempv[0] = temp1*temp[3];
258 tempv[1] = temp1*temp[4];
259 tempv[2] = temp1*temp[5];
260 gpu_cross3(a1+3,tempv,tempv2);
261 teta[0] += tempv2[0];
262 teta[1] += tempv2[1];
263 teta[2] += tempv2[2];
264
265 tempv[0] = temp1*temp[6];
266 tempv[1] = temp1*temp[7];
267 tempv[2] = temp1*temp[8];
268 gpu_cross3(a1+6,tempv,tempv2);
269 teta[0] += tempv2[0];
270 teta[1] += tempv2[1];
271 teta[2] += tempv2[2];
272 }
273
274 numtyp chi, dchi[3], tchi[3];
275 { // Compute chi and dchi
276
277 // Compute b12
278 numtyp b2[9], b12[9];
279 {
280 gpu_diag_times3(well[jtype],a2,b12);
281 gpu_transpose_times3(a2,b12,b2);
282 gpu_plus3(b1,b2,b12);
283 }
284
285 // compute chi_12
286 r12[0]*=r;
287 r12[1]*=r;
288 r12[2]*=r;
289 numtyp iota[3];
290 gpu_mldivide3(b12,r12,iota,err_flag);
291 // -- iota is now iota/r
292 iota[0]*=ir;
293 iota[1]*=ir;
294 iota[2]*=ir;
295 r12[0]*=ir;
296 r12[1]*=ir;
297 r12[2]*=ir;
298 chi = gpu_dot3(r12,iota);
299 chi = ucl_powr(chi*(numtyp)2.0,gum[2]);
300
301 // -- iota is now ok
302 iota[0]*=r;
303 iota[1]*=r;
304 iota[2]*=r;
305
306 numtyp temp1 = gpu_dot3(iota,r12);
307 numtyp temp2 = (numtyp)-4.0*ir*ir*gum[2]*ucl_powr(chi,(gum[2]-(numtyp)1.0)/
308 gum[2]);
309 dchi[0] = temp2*(iota[0]-temp1*r12[0]);
310 dchi[1] = temp2*(iota[1]-temp1*r12[1]);
311 dchi[2] = temp2*(iota[2]-temp1*r12[2]);
312
313 // compute t_chi
314 numtyp tempv[3];
315 gpu_row_times3(iota,b1,tempv);
316 gpu_cross3(tempv,iota,tchi);
317 tchi[0] *= temp2;
318 tchi[1] *= temp2;
319 tchi[2] *= temp2;
320 }
321
322 numtyp temp2 = factor_lj*eta*chi;
323 if (EVFLAG && eflag)
324 energy+=u_r*temp2;
325 numtyp temp1 = -eta*u_r*factor_lj;
326 if (EVFLAG && vflag) {
327 r12[0]*=-r;
328 r12[1]*=-r;
329 r12[2]*=-r;
330 numtyp ft=temp1*dchi[0]-temp2*dUr[0];
331 f.x+=ft;
332 virial[0]+=r12[0]*ft;
333 ft=temp1*dchi[1]-temp2*dUr[1];
334 f.y+=ft;
335 virial[1]+=r12[1]*ft;
336 virial[3]+=r12[0]*ft;
337 ft=temp1*dchi[2]-temp2*dUr[2];
338 f.z+=ft;
339 virial[2]+=r12[2]*ft;
340 virial[4]+=r12[0]*ft;
341 virial[5]+=r12[1]*ft;
342 } else {
343 f.x+=temp1*dchi[0]-temp2*dUr[0];
344 f.y+=temp1*dchi[1]-temp2*dUr[1];
345 f.z+=temp1*dchi[2]-temp2*dUr[2];
346 }
347
348 // Torque on 1
349 temp1 = -u_r*eta*factor_lj;
350 temp2 = -u_r*chi*factor_lj;
351 numtyp temp3 = -chi*eta*factor_lj;
352 tor.x+=temp1*tchi[0]+temp2*teta[0]+temp3*tUr[0];
353 tor.y+=temp1*tchi[1]+temp2*teta[1]+temp3*tUr[1];
354 tor.z+=temp1*tchi[2]+temp2*teta[2]+temp3*tUr[2];
355
356 } // for nbor
357 } // if ii
358 store_answers_t(f,tor,energy,virial,ii,astride,tid,t_per_atom,offset,eflag,
359 vflag,ans,engv,inum);
360 }
361
362