1 /* { dg-do compile } */
2 /* { dg-options "-w -O3 -mcpu=cortex-a53" } */
3 typedef struct __sFILE __FILE;
4 typedef __FILE FILE;
5 typedef int atom_id;
6 typedef float real;
7 typedef real rvec[3];
8 typedef real matrix[3][3];
9 enum {
10 ebCGS,ebMOLS,ebSBLOCKS,ebNR
11 };
12 enum {
13 efepNO, efepYES, efepNR
14 };
15 enum {
16 esolNO, esolMNO, esolWATER, esolWATERWATER, esolNR
17 };
18 typedef struct {
19 int nr;
20 atom_id *index;
21 atom_id *a;
22 } t_block;
23 enum {
24 F_LJ,
25 F_LJLR,
26 F_SR,
27 F_LR,
28 F_DVDL,
29 };
30 typedef struct {
31 t_block excl;
32 } t_atoms;
33 typedef struct {
34 t_atoms atoms;
35 t_block blocks[ebNR];
36 } t_topology;
37 typedef struct {
38 } t_nsborder;
39 extern FILE *debug;
40 typedef struct {
41 } t_nrnb;
42 typedef struct {
43 int nri,maxnri;
44 int nrj,maxnrj;
45 int maxlen;
46 int solvent;
47 int *gid;
48 int *jindex;
49 atom_id *jjnr;
50 int *nsatoms;
51 } t_nblist;
52 typedef struct {
53 int nrx,nry,nrz;
54 } t_grid;
55 typedef struct {
56 } t_commrec;
57 enum { eNL_VDWQQ, eNL_VDW, eNL_QQ,
58 eNL_VDWQQ_FREE, eNL_VDW_FREE, eNL_QQ_FREE,
59 eNL_VDWQQ_SOLMNO, eNL_VDW_SOLMNO, eNL_QQ_SOLMNO,
60 eNL_VDWQQ_WATER, eNL_QQ_WATER,
61 eNL_VDWQQ_WATERWATER, eNL_QQ_WATERWATER,
62 eNL_NR };
63 typedef struct {
64 real rlist,rlistlong;
65 real rcoulomb_switch,rcoulomb;
66 real rvdw_switch,rvdw;
67 int efep;
68 int cg0,hcg;
69 int *solvent_type;
70 int *mno_index;
71 rvec *cg_cm;
72 t_nblist nlist_sr[eNL_NR];
73 t_nblist nlist_lr[eNL_NR];
74 int bTwinRange;
75 rvec *f_twin;
76 int *eg_excl;
77 } t_forcerec;
78 typedef struct {
79 real *chargeA,*chargeB,*chargeT;
80 int *bPerturbed;
81 int *typeA,*typeB;
82 unsigned short *cTC,*cENER,*cACC,*cFREEZE,*cXTC,*cVCM;
83 } t_mdatoms;
84 enum { egCOUL, egLJ, egBHAM, egLR, egLJLR, egCOUL14, egLJ14, egNR };
85 typedef struct {
86 real *ee[egNR];
87 } t_grp_ener;
88 typedef struct {
89 t_grp_ener estat;
90 } t_groups;
91 typedef unsigned long t_excl;
reset_nblist(t_nblist * nl)92 static void reset_nblist(t_nblist *nl)
93 {
94 nl->nri = 0;
95 nl->nrj = 0;
96 nl->maxlen = 0;
97 if (nl->maxnri > 0) {
98 nl->gid[0] = -1;
99 if (nl->maxnrj > 1) {
100 nl->jindex[0] = 0;
101 nl->jindex[1] = 0;
102 }
103 }
104 }
reset_neighbor_list(t_forcerec * fr,int bLR,int eNL)105 static void reset_neighbor_list(t_forcerec *fr,int bLR,int eNL)
106 {
107 reset_nblist(&(fr->nlist_lr[eNL]));
108 }
close_i_nblist(t_nblist * nlist)109 static void close_i_nblist(t_nblist *nlist)
110 {
111 int nri = nlist->nri;
112 int len;
113 nlist->jindex[nri+1] = nlist->nrj;
114 len=nlist->nrj - nlist->jindex[nri];
115 if (nlist->solvent==esolMNO)
116 len *= nlist->nsatoms[3*nri];
117 if(len > nlist->maxlen)
118 nlist->maxlen = len;
119 }
close_nblist(t_nblist * nlist)120 static void close_nblist(t_nblist *nlist)
121 {
122 if (nlist->maxnri > 0) {
123 int nri = nlist->nri;
124 if ((nlist->jindex[nri+1] > nlist->jindex[nri]) &&
125 (nlist->gid[nri] != -1)) {
126 nlist->nri++;
127 nlist->jindex[nri+2] = nlist->nrj;
128 }
129 }
130 }
close_neighbor_list(t_forcerec * fr,int bLR,int eNL)131 static void close_neighbor_list(t_forcerec *fr,int bLR,int eNL)
132 {
133 close_nblist(&(fr->nlist_lr[eNL]));
134 }
add_j_to_nblist(t_nblist * nlist,atom_id j_atom)135 static void add_j_to_nblist(t_nblist *nlist,atom_id j_atom)
136 {
137 int nrj=nlist->nrj;
138 nlist->jjnr[nrj] = j_atom;
139 nlist->nrj ++;
140 }
put_in_list(int bHaveLJ[],int ngid,t_mdatoms * md,int icg,int jgid,int nj,atom_id jjcg[],atom_id index[],t_excl bExcl[],int shift,t_forcerec * fr,int bLR,int bVDWOnly,int bCoulOnly)141 static void put_in_list(int bHaveLJ[],
142 int ngid,t_mdatoms *md,
143 int icg,int jgid,int nj,atom_id jjcg[],
144 atom_id index[],
145 t_excl bExcl[],int shift,
146 t_forcerec *fr,int bLR,
147 int bVDWOnly,int bCoulOnly)
148 {
149 t_nblist *vdwc,*vdw,*coul;
150 t_nblist *vdwc_ww=((void *)0),*coul_ww=((void *)0);
151 t_nblist *vdwc_free=((void *)0),*vdw_free=((void *)0),*coul_free=((void *)0);
152 int i,j,jcg,igid,gid,ind_ij;
153 atom_id jj,jj0,jj1,i_atom;
154 int i0,nicg,len;
155 int *type,*typeB;
156 unsigned short *cENER;
157 real *charge,*chargeB;
158 real qi,qiB,qq,rlj;
159 int bWater,bMNO,bFree,bFreeJ,bNotEx,*bPert;
160 charge = md->chargeA;
161 chargeB = md->chargeB;
162 type = md->typeA;
163 typeB = md->typeB;
164 cENER = md->cENER;
165 bPert = md->bPerturbed;
166 i0 = index[icg];
167 nicg = index[icg+1]-i0;
168 bMNO = (fr->solvent_type[icg] == esolMNO);
169 if (bLR) {
170 if (bWater) {
171 vdw = &fr->nlist_lr[eNL_VDW];
172 coul = &fr->nlist_lr[eNL_QQ_WATER];
173 vdwc_ww = &fr->nlist_lr[eNL_VDWQQ_WATERWATER];
174 } else if(bMNO) {
175 vdwc = &fr->nlist_lr[eNL_VDWQQ_SOLMNO];
176 }
177 if (fr->efep != efepNO) {
178 vdw_free = &fr->nlist_lr[eNL_VDW_FREE];
179 coul_free = &fr->nlist_lr[eNL_QQ_FREE];
180 }
181 }
182 else {
183 if (bWater) {
184 } else if(bMNO) {
185 vdwc = &fr->nlist_sr[eNL_VDWQQ_SOLMNO];
186 }
187 if (fr->efep != efepNO) {
188 vdwc_free = &fr->nlist_sr[eNL_VDWQQ_FREE];
189 }
190 }
191 if (fr->efep==efepNO) {
192 if (bWater) {
193 igid = cENER[i_atom];
194 gid = ((igid < jgid) ? (igid*ngid+jgid) : (jgid*ngid+igid));
195 if (!bCoulOnly && !bVDWOnly) {
196 new_i_nblist(vdwc,bLR ? F_LJLR : F_LJ,i_atom,shift,gid,((void *)0));
197 new_i_nblist(vdwc_ww,bLR ? F_LJLR : F_LJ,i_atom,shift,gid,((void *)0));
198 }
199 if (!bCoulOnly)
200 new_i_nblist(vdw,bLR ? F_LJLR : F_LJ,i_atom,shift,gid,((void *)0));
201 if (!bVDWOnly) {
202 new_i_nblist(coul,bLR ? F_LR : F_SR,i_atom,shift,gid,((void *)0));
203 new_i_nblist(coul_ww,bLR ? F_LR : F_SR,i_atom,shift,gid,((void *)0));
204 }
205 for(j=0; (j<nj); j++) {
206 jcg=jjcg[j];
207 if (jcg==icg)
208 jj0 = index[jcg];
209 if (bWater && (fr->solvent_type[jcg] == esolWATER)) {
210 if (bVDWOnly)
211 add_j_to_nblist(vdw,jj0);
212 else {
213 add_j_to_nblist(coul_ww,jj0);
214 add_j_to_nblist(vdwc_ww,jj0);
215 }
216 } else {
217 jj1 = index[jcg+1];
218 if (bCoulOnly) {
219 for(jj=jj0; (jj<jj1); jj++) {
220 if (fabs(charge[jj]) > 1.2e-38)
221 add_j_to_nblist(coul,jj);
222 }
223 } else if (bVDWOnly) {
224 for(jj=jj0; (jj<jj1); jj++)
225 if (bHaveLJ[type[jj]])
226 add_j_to_nblist(vdw,jj);
227 } else {
228 for(jj=jj0; (jj<jj1); jj++) {
229 if (bHaveLJ[type[jj]]) {
230 if (fabs(charge[jj]) > 1.2e-38)
231 add_j_to_nblist(vdwc,jj);
232 add_j_to_nblist(vdw,jj);
233 } else if (fabs(charge[jj]) > 1.2e-38)
234 add_j_to_nblist(coul,jj);
235 }
236 }
237 }
238 }
239 close_i_nblist(vdw);
240 close_i_nblist(coul);
241 close_i_nblist(vdwc);
242 close_i_nblist(coul_ww);
243 close_i_nblist(vdwc_ww);
244 } else if (bMNO) {
245 igid = cENER[i_atom];
246 gid = ((igid < jgid) ? (igid*ngid+jgid) : (jgid*ngid+igid));
247 if (!bCoulOnly && !bVDWOnly)
248 new_i_nblist(vdwc,bLR ? F_LJLR : F_LJ,i_atom,shift,gid,
249 &(fr->mno_index[icg*3]));
250 if (!bCoulOnly)
251 new_i_nblist(vdw,bLR ? F_LR : F_SR,i_atom,shift,gid,
252 &(fr->mno_index[icg*3]));
253 if (!bVDWOnly)
254 new_i_nblist(coul,bLR ? F_LR : F_SR,i_atom,shift,gid,
255 &(fr->mno_index[icg*3]));
256 for(j=0; (j<nj); j++) {
257 jcg=jjcg[j];
258 if (jcg == icg)
259 jj0 = index[jcg];
260 jj1=index[jcg+1];
261 for(jj=jj0; (jj<jj1); jj++) {
262 if (bCoulOnly) {
263 if (fabs(charge[jj]) > 1.2e-38)
264 add_j_to_nblist(coul,jj);
265 } else if (bVDWOnly) {
266 if (bHaveLJ[type[jj]])
267 add_j_to_nblist(vdw,jj);
268 } else {
269 if (bHaveLJ[type[jj]]) {
270 if (fabs(charge[jj]) > 1.2e-38)
271 add_j_to_nblist(vdwc,jj);
272 add_j_to_nblist(vdw,jj);
273 } else if (fabs(charge[jj]) > 1.2e-38)
274 add_j_to_nblist(coul,jj);
275 }
276 }
277 close_i_nblist(vdw);
278 close_i_nblist(coul);
279 close_i_nblist(vdwc);
280 }
281 } else {
282 for(i=0; i<nicg; i++) {
283 igid = cENER[i_atom];
284 gid = ((igid < jgid) ? (igid*ngid+jgid) : (jgid*ngid+igid));
285 qi = charge[i_atom];
286 if (!bCoulOnly && !bVDWOnly)
287 new_i_nblist(vdwc,bLR ? F_LJLR : F_LJ,i_atom,shift,gid,((void *)0));
288 if (!bCoulOnly)
289 new_i_nblist(vdw,bLR ? F_LR : F_SR,i_atom,shift,gid,((void *)0));
290 if (!bVDWOnly)
291 new_i_nblist(coul,bLR ? F_LR : F_SR,i_atom,shift,gid,((void *)0));
292 if (!(bVDWOnly || fabs(qi)<1.2e-38) || !(bCoulOnly || !bHaveLJ[type[i_atom]])) {
293 for(j=0; (j<nj); j++) {
294 jcg=jjcg[j];
295 if (jcg == icg)
296 jj0 = i0 + i + 1;
297 else
298 jj0 = index[jcg];
299 jj1=index[jcg+1];
300 for(jj=jj0; jj<jj1; jj++) {
301 bNotEx = !((int) ((bExcl)[((atom_id) (jj))] & (1<<((atom_id) (i)))));
302 if (bNotEx) {
303 if (bCoulOnly) {
304 if (fabs(charge[jj]) > 1.2e-38)
305 add_j_to_nblist(coul,jj);
306 } else if (bVDWOnly) {
307 if (bHaveLJ[type[jj]])
308 add_j_to_nblist(vdw,jj);
309 } else {
310 if (bHaveLJ[type[jj]]) {
311 if (fabs(qi) > 1.2e-38 && (fabs(charge[jj]) > 1.2e-38))
312 add_j_to_nblist(vdwc,jj);
313 add_j_to_nblist(vdw,jj);
314 } else if (fabs(qi) > 1.2e-38 && (fabs(charge[jj]) > 1.2e-38))
315 add_j_to_nblist(coul,jj);
316 }
317 }
318 }
319 }
320 }
321 close_i_nblist(vdw);
322 close_i_nblist(coul);
323 close_i_nblist(vdwc);
324 }
325 }
326 } else {
327 for(i=0; i<nicg; i++) {
328 igid = cENER[i_atom];
329 gid = ((igid < jgid) ? (igid*ngid+jgid) : (jgid*ngid+igid));
330 qi = charge[i_atom];
331 qiB = chargeB[i_atom];
332 if (!bCoulOnly && !bVDWOnly)
333 new_i_nblist(vdwc,bLR ? F_LJLR : F_LJ,i_atom,shift,gid,
334 bMNO ? &(fr->mno_index[icg*3]) : ((void *)0));
335 if (!bCoulOnly)
336 new_i_nblist(vdw,bLR ? F_LR : F_SR,i_atom,shift,gid,
337 bMNO ? &(fr->mno_index[icg*3]) : ((void *)0));
338 new_i_nblist(coul,bLR ? F_LR : F_SR,i_atom,shift,gid,
339 bMNO ? &(fr->mno_index[icg*3]) : ((void *)0));
340 new_i_nblist(vdw_free,F_DVDL,i_atom,shift,gid,((void *)0));
341 new_i_nblist(coul_free,F_DVDL,i_atom,shift,gid,((void *)0));
342 new_i_nblist(vdwc_free,F_DVDL,i_atom,shift,gid,((void *)0));
343 if (!(bVDWOnly || (fabs(qi)<1.2e-38 && fabs(qiB)<1.2e-38)) ||
344 !(bCoulOnly || (!bHaveLJ[type[i_atom]] && !bHaveLJ[typeB[i_atom]]))) {
345 for(j=0; (j<nj); j++) {
346 jcg=jjcg[j];
347 if (jcg == icg)
348 jj0 = i0 + i + 1;
349 else
350 jj0 = index[jcg];
351 jj1=index[jcg+1];
352 bFree = bPert[i_atom];
353 for(jj=jj0; (jj<jj1); jj++) {
354 bFreeJ = bFree || bPert[jj];
355 if ((!bWater && !bMNO) || i==0 || bFreeJ) {
356 bNotEx = !((int) ((bExcl)[((atom_id) (jj))] & (1<<((atom_id) (i)))));
357 if (bNotEx) {
358 if (bFreeJ) {
359 if (bCoulOnly)
360 add_j_to_nblist(coul_free,jj);
361 else if (bVDWOnly)
362 add_j_to_nblist(vdw_free,jj);
363 add_j_to_nblist(vdwc_free,jj);
364 } else if (bCoulOnly) {
365 add_j_to_nblist(coul,jj);
366 } else if (bVDWOnly) {
367 if (bHaveLJ[type[jj]])
368 add_j_to_nblist(vdw,jj);
369 } else {
370 if (bHaveLJ[type[jj]]) {
371 if (fabs(qi) > 1.2e-38 && (fabs(charge[jj]) > 1.2e-38))
372 add_j_to_nblist(vdwc,jj);
373 add_j_to_nblist(vdw,jj);
374 } else if (fabs(qi) > 1.2e-38 && (fabs(charge[jj]) > 1.2e-38))
375 add_j_to_nblist(coul,jj);
376 }
377 }
378 }
379 }
380 }
381 }
382 close_i_nblist(vdw);
383 close_i_nblist(coul);
384 close_i_nblist(vdwc);
385 if (bWater && (i==0)) {
386 close_i_nblist(coul_ww);
387 close_i_nblist(vdwc_ww);
388 }
389 close_i_nblist(vdw_free);
390 close_i_nblist(coul_free);
391 close_i_nblist(vdwc_free);
392 }
393 }
394 }
setexcl(atom_id start,atom_id end,t_block * excl,int b,t_excl bexcl[])395 static void setexcl(atom_id start,atom_id end,t_block *excl,int b,
396 t_excl bexcl[])
397 {
398 atom_id i,k;
399 if (b) {
400 for(i=start; i<end; i++) {
401 for(k=excl->index[i]; k<excl->index[i+1]; k++) {
402 (bexcl)[((atom_id) (excl->a[k]))] |= (1<<((atom_id) (i-start)));
403 }
404 }
405 }
406 }
calc_naaj(int icg,int cgtot)407 int calc_naaj(int icg,int cgtot)
408 {
409 int naaj;
410 if ((cgtot % 2) == 1) {
411 naaj = 1+(cgtot/2);
412 }
413 else if ((cgtot % 4) == 0) {
414 if (icg < cgtot/2) {
415 if ((icg % 2) == 0)
416 naaj=1+(cgtot/2);
417 }
418 else {
419 if ((icg % 2) == 1)
420 naaj=1+(cgtot/2);
421 }
422 }
423 else {
424 if ((icg % 2) == 0)
425 naaj=1+(cgtot/2);
426 else
427 naaj=cgtot/2;
428 }
429 return naaj;
430 }
get_dx(int Nx,real gridx,real grid_x,real rc2,real x,int * dx0,int * dx1,real * dcx2)431 static void get_dx(int Nx,real gridx,real grid_x,real rc2,real x,
432 int *dx0,int *dx1,real *dcx2)
433 {
434 real dcx,tmp;
435 int xgi,xgi0,xgi1,i;
436 xgi = (int)(Nx+x*grid_x)-Nx;
437 if (xgi < 0) {
438 *dx0 = 0;
439 *dx1 = -1;
440 } else if (xgi >= Nx) {
441 *dx0 = Nx;
442 *dx1 = Nx-1;
443 } else {
444 dcx2[xgi] = 0;
445 *dx0 = xgi;
446 xgi0 = xgi-1;
447 *dx1 = xgi;
448 xgi1 = xgi+1;
449 }
450 for(i=xgi0; i>=0; i--) {
451 dcx = (i+1)*gridx-x;
452 tmp = dcx*dcx;
453 if (tmp >= rc2)
454 *dx0 = i;
455 dcx2[i] = tmp;
456 }
457 for(i=xgi1; i<Nx; i++) {
458 dcx = i*gridx-x;
459 tmp = dcx*dcx;
460 if (tmp >= rc2)
461 *dx1 = i;
462 dcx2[i] = tmp;
463 }
464 }
do_longrange(FILE * log,t_commrec * cr,t_topology * top,t_forcerec * fr,int ngid,t_mdatoms * md,int icg,int jgid,int nlr,atom_id lr[],t_excl bexcl[],int shift,rvec x[],rvec box_size,t_nrnb * nrnb,real lambda,real * dvdlambda,t_groups * grps,int bVDWOnly,int bCoulOnly,int bDoForces,int bHaveLJ[])465 static void do_longrange(FILE *log,t_commrec *cr,t_topology *top,t_forcerec *fr,
466 int ngid,t_mdatoms *md,int icg,
467 int jgid,int nlr,
468 atom_id lr[],t_excl bexcl[],int shift,
469 rvec x[],rvec box_size,t_nrnb *nrnb,
470 real lambda,real *dvdlambda,
471 t_groups *grps,int bVDWOnly,int bCoulOnly,
472 int bDoForces,int bHaveLJ[])
473 {
474 int i;
475 for(i=0; (i<eNL_NR); i++) {
476 if ((fr->nlist_lr[i].nri > fr->nlist_lr[i].maxnri-32) || bDoForces) {
477 close_neighbor_list(fr,1,i);
478 do_fnbf(log,cr,fr,x,fr->f_twin,md,
479 grps->estat.ee[egLJLR],grps->estat.ee[egLR],box_size,
480 nrnb,lambda,dvdlambda,1,i);
481 reset_neighbor_list(fr,1,i);
482 }
483 }
484 if (!bDoForces) {
485 put_in_list(bHaveLJ,ngid,md,icg,jgid,nlr,lr,top->blocks[ebCGS].index,
486 bexcl,shift,fr,
487 1,bVDWOnly,bCoulOnly);
488 }
489 }
ns5_core(FILE * log,t_commrec * cr,t_forcerec * fr,int cg_index[],matrix box,rvec box_size,int ngid,t_topology * top,t_groups * grps,t_grid * grid,rvec x[],t_excl bexcl[],int * bExcludeAlleg,t_nrnb * nrnb,t_mdatoms * md,real lambda,real * dvdlambda,int bHaveLJ[])490 static int ns5_core(FILE *log,t_commrec *cr,t_forcerec *fr,int cg_index[],
491 matrix box,rvec box_size,int ngid,
492 t_topology *top,t_groups *grps,
493 t_grid *grid,rvec x[],t_excl bexcl[],int *bExcludeAlleg,
494 t_nrnb *nrnb,t_mdatoms *md,
495 real lambda,real *dvdlambda,
496 int bHaveLJ[])
497 {
498 static atom_id **nl_lr_ljc,**nl_lr_one,**nl_sr=((void *)0);
499 static int *nlr_ljc,*nlr_one,*nsr;
500 static real *dcx2=((void *)0),*dcy2=((void *)0),*dcz2=((void *)0);
501 t_block *cgs=&(top->blocks[ebCGS]);
502 unsigned short *gid=md->cENER;
503 int tx,ty,tz,dx,dy,dz,cj;
504 int dx0,dx1,dy0,dy1,dz0,dz1;
505 int Nx,Ny,Nz,shift=-1,j,nrj,nns,nn=-1;
506 real gridx,gridy,gridz,grid_x,grid_y,grid_z;
507 int icg=-1,iicg,cgsnr,i0,nri,naaj,min_icg,icg_naaj,jjcg,cgj0,jgid;
508 int bVDWOnly,bCoulOnly;
509 rvec xi,*cgcm;
510 real r2,rs2,rvdw2,rcoul2,rm2,rl2,XI,YI,ZI,dcx,dcy,dcz,tmp1,tmp2;
511 int *i_eg_excl;
512 int use_twinrange,use_two_cutoffs;
513 cgsnr = cgs->nr;
514 rs2 = ((fr->rlist)*(fr->rlist));
515 if (fr->bTwinRange) {
516 rvdw2 = ((fr->rvdw)*(fr->rvdw));
517 rcoul2 = ((fr->rcoulomb)*(fr->rcoulomb));
518 } else {
519 }
520 rm2 = (((rvdw2) < (rcoul2)) ? (rvdw2) : (rcoul2) );
521 rl2 = (((rvdw2) > (rcoul2)) ? (rvdw2) : (rcoul2) );
522 use_twinrange = (rs2 < rm2);
523 use_two_cutoffs = (rm2 < rl2);
524 bVDWOnly = (rvdw2 > rcoul2);
525 bCoulOnly = !bVDWOnly;
526 if (nl_sr == ((void *)0)) {
527 (nl_sr)=save_calloc("nl_sr","ns.c",1341, (ngid),sizeof(*(nl_sr)));
528 (nsr)=save_calloc("nsr","ns.c",1343, (ngid),sizeof(*(nsr)));
529 (nlr_ljc)=save_calloc("nlr_ljc","ns.c",1344, (ngid),sizeof(*(nlr_ljc)));
530 (nlr_one)=save_calloc("nlr_one","ns.c",1345, (ngid),sizeof(*(nlr_one)));
531 if (use_twinrange)
532 (nl_lr_ljc)=save_calloc("nl_lr_ljc","ns.c",1349, (ngid),sizeof(*(nl_lr_ljc)));
533 if (use_two_cutoffs)
534 (nl_lr_one)=save_calloc("nl_lr_one","ns.c",1353, (ngid),sizeof(*(nl_lr_one)));
535 for(j=0; (j<ngid); j++) {
536 (nl_sr[j])=save_calloc("nl_sr[j]","ns.c",1356, (1024),sizeof(*(nl_sr[j])));
537 if (use_twinrange)
538 (nl_lr_ljc[j])=save_calloc("nl_lr_ljc[j]","ns.c",1358, (1024),sizeof(*(nl_lr_ljc[j])));
539 if (use_two_cutoffs)
540 (nl_lr_one[j])=save_calloc("nl_lr_one[j]","ns.c",1360, (1024),sizeof(*(nl_lr_one[j])));
541 }
542 if (debug)
543 fprintf(debug,"ns5_core: rs2 = %g, rvdw2 = %g, rcoul2 = %g (nm^2)\n",
544 rs2,rvdw2,rcoul2);
545 }
546 cgcm = fr->cg_cm;
547 Nx = grid->nrx;
548 Ny = grid->nry;
549 if (dcx2 == ((void *)0)) {
550 (dcx2)=save_calloc("dcx2","ns.c",1379, (Nx*2),sizeof(*(dcx2)));
551 (dcy2)=save_calloc("dcy2","ns.c",1380, (Ny*2),sizeof(*(dcy2)));
552 (dcz2)=save_calloc("dcz2","ns.c",1381, (Nz*2),sizeof(*(dcz2)));
553 }
554 gridx = box[0][0]/grid->nrx;
555 gridy = box[1][1]/grid->nry;
556 gridz = box[2][2]/grid->nrz;
557 grid_x = 1/gridx;
558 grid_y = 1/gridy;
559 grid_z = 1/gridz;
560 for(iicg=fr->cg0; (iicg < fr->hcg); iicg++) {
561 icg = cg_index[iicg];
562 if (icg != iicg)
563 fatal_error(0,"icg = %d, iicg = %d, file %s, line %d",icg,iicg,"ns.c",
564 1408);
565 if(bExcludeAlleg[icg])
566 i_eg_excl = fr->eg_excl + ngid*gid[cgs->index[icg]];
567 setexcl(cgs->index[icg],cgs->index[icg+1],&top->atoms.excl,1,bexcl);
568 naaj = calc_naaj(icg,cgsnr);
569 icg_naaj = icg+naaj;
570 for (tz=-1; tz<=1; tz++) {
571 ZI = cgcm[icg][2]+tz*box[2][2];
572 get_dx(Nz,gridz,grid_z,rcoul2,ZI,&dz0,&dz1,dcz2);
573 if (dz0 > dz1)
574 for (ty=-1; ty<=1; ty++) {
575 YI = cgcm[icg][1]+ty*box[1][1]+tz*box[2][1];
576 get_dx(Ny,gridy,grid_y,rcoul2,YI,&dy0,&dy1,dcy2);
577 for (tx=-1; tx<=1; tx++) {
578 get_dx(Nx,gridx,grid_x,rcoul2,XI,&dx0,&dx1,dcx2);
579 shift=((2*1 +1)*((2*1 +1)*((tz)+1)+(ty)+1)+(tx)+1);
580 for (dx=dx0; (dx<=dx1); dx++) {
581 for (dy=dy0; (dy<=dy1); dy++) {
582 for (dz=dz0; (dz<=dz1); dz++) {
583 if (tmp2 > dcz2[dz]) {
584 for (j=0; (j<nrj); j++) {
585 if (((jjcg >= icg) && (jjcg < icg_naaj)) ||
586 ((jjcg < min_icg))) {
587 if (r2 < rl2) {
588 if (!i_eg_excl[jgid]) {
589 if (r2 < rs2) {
590 if (nsr[jgid] >= 1024) {
591 put_in_list(bHaveLJ,ngid,md,icg,jgid,
592 nsr[jgid],nl_sr[jgid],
593 cgs->index, bexcl,
594 shift,fr,0,0,0);
595 }
596 } else if (r2 < rm2) {
597 } else if (use_two_cutoffs) {
598 if (nlr_one[jgid] >= 1024) {
599 do_longrange(log,cr,top,fr,ngid,md,icg,jgid,
600 nlr_one[jgid],
601 nl_lr_one[jgid],bexcl,shift,x,
602 box_size,nrnb,
603 lambda,dvdlambda,grps,
604 bVDWOnly,bCoulOnly,0,
605 bHaveLJ);
606 }
607 }
608 }
609 }
610 }
611 }
612 }
613 }
614 }
615 }
616 }
617 }
618 }
619 }
620 }
search_neighbours(FILE * log,t_forcerec * fr,rvec x[],matrix box,t_topology * top,t_groups * grps,t_commrec * cr,t_nsborder * nsb,t_nrnb * nrnb,t_mdatoms * md,real lambda,real * dvdlambda)621 int search_neighbours(FILE *log,t_forcerec *fr,
622 rvec x[],matrix box,
623 t_topology *top,t_groups *grps,
624 t_commrec *cr,t_nsborder *nsb,
625 t_nrnb *nrnb,t_mdatoms *md,
626 real lambda,real *dvdlambda)
627 {
628 static t_grid *grid=((void *)0);
629 static t_excl *bexcl;
630 static int *bHaveLJ;
631 static int *cg_index=((void *)0),*slab_index=((void *)0);
632 static int *bExcludeAlleg;
633 rvec box_size;
634 int i,j,m,ngid;
635 int nsearch;
636 nsearch = ns5_core(log,cr,fr,cg_index,box,box_size,ngid,top,grps,
637 grid,x,bexcl,bExcludeAlleg,nrnb,md,lambda,dvdlambda,bHaveLJ);
638 }
639