1 /**************************************************************************
2 lj_tip4p_long.cpp
3 -------------------
4 V. Nikolskiy (HSE)
5
6 Class for acceleration of the lj/tip4p/long pair style
7
8 __________________________________________________________________________
9 This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
10 __________________________________________________________________________
11
12 begin :
13 email : thevsevak@gmail.com
14 ***************************************************************************/
15
16 #if defined(USE_OPENCL)
17 #include "lj_tip4p_long_cl.h"
18 #elif defined(USE_CUDART)
19 const char *lj_tip4p=0;
20 #else
21 #include "lj_tip4p_long_cubin.h"
22 #endif
23
24 #include "lal_lj_tip4p_long.h"
25 #include <cassert>
26 namespace LAMMPS_AL {
27 #define LJTIP4PLongT LJ_TIP4PLong<numtyp, acctyp>
28
29 extern Device<PRECISION,ACC_PRECISION> device;
30
31 template <class numtyp, class acctyp>
LJ_TIP4PLong()32 LJTIP4PLongT::LJ_TIP4PLong(): BaseCharge<numtyp,acctyp>(), _allocated(false) {
33 }
34
35 template <class numtyp, class acctyp>
~LJ_TIP4PLong()36 LJTIP4PLongT::~LJ_TIP4PLong() {
37 clear();
38 }
39
40 template <class numtyp, class acctyp>
bytes_per_atom(const int max_nbors) const41 int LJTIP4PLongT::bytes_per_atom(const int max_nbors) const {
42 return this->bytes_per_atom_atomic(max_nbors);
43 }
44
45 template <class numtyp, class acctyp>
init(const int ntypes,double ** host_cutsq,double ** host_lj1,double ** host_lj2,double ** host_lj3,double ** host_lj4,double ** host_offset,double * host_special_lj,const int nlocal,const int tH,const int tO,const double a,const double qd,const int nall,const int max_nbors,const int maxspecial,const double cell_size,const double gpu_split,FILE * _screen,double ** host_cut_ljsq,const double host_cut_coulsq,const double host_cut_coulsqplus,double * host_special_coul,const double qqrd2e,const double g_ewald,int map_size,int max_same)46 int LJTIP4PLongT::init(const int ntypes,
47 double **host_cutsq, double **host_lj1,
48 double **host_lj2, double **host_lj3,
49 double **host_lj4, double **host_offset,
50 double *host_special_lj, const int nlocal,
51 const int tH, const int tO,
52 const double a, const double qd,
53 const int nall, const int max_nbors,
54 const int maxspecial, const double cell_size,
55 const double gpu_split, FILE *_screen,
56 double **host_cut_ljsq,
57 const double host_cut_coulsq, const double host_cut_coulsqplus,
58 double *host_special_coul, const double qqrd2e,
59 const double g_ewald, int map_size, int max_same) {
60 int success;
61 success=this->init_atomic(nlocal,nall,max_nbors,maxspecial,cell_size,gpu_split,
62 _screen,lj_tip4p_long,"k_lj_tip4p_long");
63 if (success!=0)
64 return success;
65 k_pair_distrib.set_function(*this->pair_program,"k_lj_tip4p_long_distrib");
66 k_pair_reneigh.set_function(*this->pair_program,"k_lj_tip4p_reneigh");
67 k_pair_newsite.set_function(*this->pair_program,"k_lj_tip4p_newsite");
68 #if defined(LAL_OCL_EV_JIT)
69 k_pair_distrib_noev.set_function(*this->pair_program_noev,
70 "k_lj_tip4p_long_distrib");
71 #else
72 k_pair_dt_sel = &k_pair_distrib;
73 #endif
74
75 TypeH = tH;
76 TypeO = tO;
77 alpha = a;
78 qdist = qd;
79
80 // If atom type constants fit in shared memory use fast kernel
81 int lj_types=ntypes;
82 shared_types=false;
83 int max_shared_types=this->device->max_shared_types();
84 if (lj_types<=max_shared_types && this->_block_size>=max_shared_types) {
85 lj_types=max_shared_types;
86 shared_types=true;
87 }
88 _lj_types=lj_types;
89
90 // Allocate a host write buffer for data initialization
91 UCL_H_Vec<numtyp> host_write(lj_types*lj_types*32,*(this->ucl_device),
92 UCL_WRITE_ONLY);
93
94 for (int i=0; i<lj_types*lj_types; i++)
95 host_write[i]=0.0;
96
97 lj1.alloc(lj_types*lj_types,*(this->ucl_device),UCL_READ_ONLY);
98 this->atom->type_pack4(ntypes,lj_types,lj1,host_write,host_lj1,host_lj2,
99 host_cut_ljsq);
100
101 lj3.alloc(lj_types*lj_types,*(this->ucl_device),UCL_READ_ONLY);
102 this->atom->type_pack4(ntypes,lj_types,lj3,host_write,host_lj3,host_lj4,
103 host_offset);
104
105 cutsq.alloc(lj_types*lj_types,*(this->ucl_device),UCL_READ_ONLY);
106 this->atom->type_pack1(ntypes,lj_types,cutsq,host_write,host_cutsq);
107
108 sp_lj.alloc(8,*(this->ucl_device),UCL_READ_ONLY);
109 for (int i=0; i<4; i++) {
110 host_write[i]=host_special_lj[i];
111 host_write[i+4]=host_special_coul[i];
112 }
113 ucl_copy(sp_lj,host_write,8,false);
114
115 //force_comp.alloc(72*72, *(this->ucl_device), UCL_READ_WRITE);
116
117 _qqrd2e=qqrd2e;
118 _g_ewald=g_ewald;
119 cut_coulsq = host_cut_coulsq;
120 cut_coulsqplus = host_cut_coulsqplus;
121
122 hneight.alloc(nall*4,*(this->ucl_device), UCL_READ_WRITE);
123 m.alloc(nall,*(this->ucl_device), UCL_READ_WRITE);
124 ansO.alloc(nall,*(this->ucl_device), UCL_READ_WRITE);
125
126 this->tag.alloc(nall,*(this->ucl_device), UCL_READ_ONLY);
127 this->atom_sametag.alloc(max_same, *(this->ucl_device), UCL_READ_ONLY);
128 this->map_array.alloc(map_size,*(this->ucl_device), UCL_READ_ONLY);
129
130 _allocated=true;
131 this->_max_bytes=lj1.row_bytes()+lj3.row_bytes()+cutsq.row_bytes()+
132 sp_lj.row_bytes() + hneight.row_bytes()+m.row_bytes()+
133 this->tag.row_bytes()+this->atom_sametag.row_bytes() +
134 this->map_array.row_bytes();
135 return 0;
136 }
137
138
139 template <class numtyp, class acctyp>
clear()140 void LJTIP4PLongT::clear() {
141 if (!_allocated)
142 return;
143 _allocated=false;
144
145 lj1.clear();
146 lj3.clear();
147 sp_lj.clear();
148 cutsq.clear();
149 hneight.clear();
150 m.clear();
151 tag.clear();
152 atom_sametag.clear();
153 map_array.clear();
154 ansO.clear();
155 //force_comp.clear();
156
157 k_pair_distrib.clear();
158 k_pair_reneigh.clear();
159 k_pair_newsite.clear();
160 #if defined(LAL_OCL_EV_JIT)
161 k_pair_distrib_noev.clear();
162 #endif
163
164 this->clear_atomic();
165 }
166
167 template <class numtyp, class acctyp>
host_memory_usage() const168 double LJTIP4PLongT::host_memory_usage() const {
169 return this->host_memory_usage_atomic()+sizeof(LJ_TIP4PLong<numtyp,acctyp>);
170 }
171
172 // ---------------------------------------------------------------------------
173 // Calculate energies, forces, and torques
174 // ---------------------------------------------------------------------------
175 template <class numtyp, class acctyp>
loop(const int eflag,const int vflag)176 int LJTIP4PLongT::loop(const int eflag, const int vflag) {
177 // Compute the block size and grid size to keep all cores busy
178 const int BX=this->block_size();
179
180 int ainum=this->ans->inum();
181 const int nall = this->atom->nall();
182 int nbor_pitch=this->nbor->nbor_pitch();
183 this->time_pair.start();
184 int GX;
185 GX=static_cast<int>(ceil(static_cast<double>(nall)/BX));
186 if (t_ago == 0) {
187 this->k_pair_reneigh.set_size(GX,BX);
188 this->k_pair_reneigh.run(&this->atom->x,
189 &this->nbor->dev_nbor, &this->_nbor_data->begin(),
190 &nall, &ainum,&nbor_pitch, &this->_threads_per_atom,
191 &hneight, &m, &TypeO, &TypeH,
192 &tag, &map_array, &atom_sametag);
193 }
194 this->k_pair_newsite.set_size(GX,BX);
195 this->k_pair_newsite.run(&this->atom->x,
196 &this->nbor->dev_nbor, &this->_nbor_data->begin(),
197 &nall, &ainum,
198 &nbor_pitch, &this->_threads_per_atom,
199 &hneight, &m, &TypeO, &TypeH, &alpha,
200 &this->atom->q);
201
202 GX=static_cast<int>(ceil(static_cast<double>(this->ans->inum())/
203 (BX/this->_threads_per_atom)));
204 if (vflag) {
205 this->ansO.resize_ib(ainum*3);
206 } else {
207 this->ansO.resize_ib(ainum);
208 }
209 this->ansO.zero();
210 this->device->gpu->sync();
211 if(shared_types) {
212 this->k_pair_sel->set_size(GX,BX);
213 this->k_pair_sel->run(&this->atom->x, &lj1, &lj3, &_lj_types, &sp_lj,
214 &this->nbor->dev_nbor, &this->_nbor_data->begin(),
215 &this->ans->force, &this->ans->engv, &eflag, &vflag,
216 &ainum, &nbor_pitch, &this->_threads_per_atom,
217 &hneight, &m, &TypeO, &TypeH, &alpha,
218 &this->atom->q, &cutsq, &_qqrd2e, &_g_ewald,
219 &cut_coulsq, &cut_coulsqplus, &this->ansO);
220 } else {
221 this->k_pair.set_size(GX,BX);
222 this->k_pair.run(&this->atom->x, &lj1, &lj3, &_lj_types, &sp_lj,
223 &this->nbor->dev_nbor, &this->_nbor_data->begin(),
224 &this->ans->force, &this->ans->engv, &eflag, &vflag,
225 &ainum, &nbor_pitch, &this->_threads_per_atom,
226 &hneight, &m, &TypeO, &TypeH, &alpha,
227 &this->atom->q, &cutsq, &_qqrd2e, &_g_ewald,
228 &cut_coulsq, &cut_coulsqplus, &this->ansO);
229 }
230 #if defined(LAL_OCL_EV_JIT)
231 if (eflag || vflag) k_pair_dt_sel = &k_pair_distrib;
232 else k_pair_dt_sel = &k_pair_distrib_noev;
233 #endif
234
235 GX=static_cast<int>(ceil(static_cast<double>(this->ans->inum())/BX));
236 k_pair_dt_sel->set_size(GX,BX);
237 k_pair_dt_sel->run(&this->atom->x, &this->ans->force, &this->ans->engv,
238 &eflag, &vflag, &ainum, &nbor_pitch,
239 &this->_threads_per_atom, &hneight, &m, &TypeO, &TypeH,
240 &alpha,&this->atom->q, &this->ansO);
241 this->time_pair.stop();
242 return GX;
243 }
244
245
246 template <class numtyp, class acctyp>
copy_relations_data(int n,tagint * tag,int * map_array,int map_size,int * sametag,int max_same,int ago)247 void LJTIP4PLongT::copy_relations_data(int n, tagint *tag, int *map_array,
248 int map_size, int *sametag, int max_same, int ago) {
249 int nall = n;
250 const int hn_sz = n*4; // matrix size = col size * col number
251 hneight.resize_ib(hn_sz);
252
253 m.resize_ib(n);
254 m.zero();
255
256 if (ago == 0) {
257 hneight.zero();
258
259
260 {
261 UCL_H_Vec<tagint> host_tag_write;
262 host_tag_write.view(tag, nall, *(this->ucl_device));
263 this->tag.resize_ib(nall);
264 ucl_copy(this->tag, host_tag_write, false);
265 }
266
267 UCL_H_Vec<int> host_write;
268 host_write.view(sametag, max_same, *(this->ucl_device));
269 this->atom_sametag.resize_ib(max_same);
270 ucl_copy(this->atom_sametag, host_write, false);
271
272 this->map_array.resize_ib(map_size);
273 host_write.view(map_array, map_size, *(this->ucl_device));
274 ucl_copy(this->map_array, host_write, false);
275 }
276 }
277
278 // ---------------------------------------------------------------------------
279 // Copy nbor list from host if necessary and then calculate forces, virials,..
280 // ---------------------------------------------------------------------------
281 template <class numtyp, class acctyp>
compute(const int f_ago,const int inum_full,const int nall,double ** host_x,int * host_type,int * ilist,int * numj,int ** firstneigh,const bool eflag_in,const bool vflag_in,const bool eatom,const bool vatom,int & host_start,const double cpu_time,bool & success,double * host_q,const int nlocal,double * boxlo,double * prd)282 void LJTIP4PLongT::compute(const int f_ago, const int inum_full,
283 const int nall, double **host_x, int *host_type,
284 int *ilist, int *numj, int **firstneigh,
285 const bool eflag_in, const bool vflag_in,
286 const bool eatom, const bool vatom,
287 int &host_start, const double cpu_time,
288 bool &success, double *host_q,
289 const int nlocal, double *boxlo, double *prd) {
290 this->acc_timers();
291 int eflag, vflag;
292 if (eflag_in) eflag=2;
293 else eflag=0;
294 if (vflag_in) vflag=2;
295 else vflag=0;
296
297 this->set_kernel(eflag,vflag);
298 if (inum_full==0) {
299 host_start=0;
300 // Make sure textures are correct if realloc by a different hybrid style
301 this->resize_atom(0,nall,success);
302 this->zero_timers();
303 return;
304 }
305
306 int ago=this->hd_balancer.ago_first(f_ago);
307 int inum=this->hd_balancer.balance(ago,inum_full,cpu_time);
308 this->ans->inum(inum);
309 host_start=inum;
310
311 if (ago==0) {
312 this->reset_nbors(nall, inum, ilist, numj, firstneigh, success);
313 if (!success)
314 return;
315 }
316
317 this->atom->cast_x_data(host_x,host_type);
318 this->atom->cast_q_data(host_q);
319 this->hd_balancer.start_timer();
320 this->atom->add_x_data(host_x,host_type);
321 this->atom->add_q_data();
322
323 this->device->precompute(f_ago,nlocal,nall,host_x,host_type,success,host_q,
324 boxlo, prd);
325
326 t_ago = ago;
327 loop(eflag,vflag);
328 this->ans->copy_answers(eflag_in,vflag_in,eatom,vatom,ilist,inum);
329 this->device->add_ans_object(this->ans);
330 this->hd_balancer.stop_timer();
331 }
332
333 // ---------------------------------------------------------------------------
334 // Reneighbor on GPU if necessary and then compute forces, virials, energies
335 // ---------------------------------------------------------------------------
336 template <class numtyp, class acctyp>
compute(const int ago,const int inum_full,const int nall,double ** host_x,int * host_type,double * sublo,double * subhi,tagint * tag,int * map_array,int map_size,int * sametag,int max_same,int ** nspecial,tagint ** special,const bool eflag_in,const bool vflag_in,const bool eatom,const bool vatom,int & host_start,int ** ilist,int ** jnum,const double cpu_time,bool & success,double * host_q,double * boxlo,double * prd)337 int** LJTIP4PLongT::compute(const int ago, const int inum_full,
338 const int nall, double **host_x, int *host_type,
339 double *sublo, double *subhi, tagint *tag,
340 int *map_array, int map_size, int *sametag,
341 int max_same, int **nspecial, tagint **special,
342 const bool eflag_in, const bool vflag_in,
343 const bool eatom, const bool vatom,
344 int &host_start, int **ilist, int **jnum,
345 const double cpu_time, bool &success,
346 double *host_q, double *boxlo, double *prd) {
347 this->acc_timers();
348 int eflag, vflag;
349 if (eflag_in) eflag=2;
350 else eflag=0;
351 if (vflag_in) vflag=2;
352 else vflag=0;
353
354 this->set_kernel(eflag,vflag);
355 if (inum_full==0) {
356 host_start=0;
357 // Make sure textures are correct if realloc by a different hybrid style
358 this->resize_atom(0,nall,success);
359 this->zero_timers();
360 return nullptr;
361 }
362
363 this->hd_balancer.balance(cpu_time);
364 int inum=this->hd_balancer.get_gpu_count(ago,inum_full);
365 this->ans->inum(inum);
366 host_start=inum;
367
368 // Build neighbor list on GPU if necessary
369 if (ago==0) {
370 this->build_nbor_list(inum, inum_full-inum, nall, host_x, host_type,
371 sublo, subhi, tag, nspecial, special, success);
372 if (!success)
373 return nullptr;
374 this->atom->cast_q_data(host_q);
375 this->hd_balancer.start_timer();
376 } else {
377 this->atom->cast_x_data(host_x,host_type);
378 this->atom->cast_q_data(host_q);
379 this->hd_balancer.start_timer();
380 this->atom->add_x_data(host_x,host_type);
381 }
382 this->atom->add_q_data();
383 *ilist=this->nbor->host_ilist.begin();
384 *jnum=this->nbor->host_acc.begin();
385
386 copy_relations_data(nall, tag, map_array, map_size, sametag, max_same, ago);
387
388 this->device->precompute(ago,inum_full,nall,host_x,host_type,success,host_q,
389 boxlo, prd);
390
391 t_ago = ago;
392 loop(eflag,vflag);
393 this->ans->copy_answers(eflag_in,vflag_in,eatom,vatom,inum);
394 this->device->add_ans_object(this->ans);
395 this->hd_balancer.stop_timer();
396
397 return this->nbor->host_jlist.begin()-host_start;
398 }
399
400
401 template class LJ_TIP4PLong<PRECISION,ACC_PRECISION>;
402 }
403