1/* { dg-do compile } */
2/* { dg-options "-w -O3 -mcpu=cortex-a53" } */
3typedef struct __sFILE __FILE;
4typedef __FILE FILE;
5typedef int atom_id;
6typedef float real;
7typedef real rvec[3];
8typedef real matrix[3][3];
9enum {
10  ebCGS,ebMOLS,ebSBLOCKS,ebNR
11};
12enum {
13  efepNO, efepYES, efepNR
14};
15enum {
16  esolNO, esolMNO, esolWATER, esolWATERWATER, esolNR
17};
18typedef struct {
19  int nr;
20  atom_id *index;
21  atom_id *a;
22} t_block;
23enum {
24  F_LJ,
25  F_LJLR,
26  F_SR,
27  F_LR,
28  F_DVDL,
29};
30typedef struct {
31  t_block excl;
32} t_atoms;
33typedef struct {
34  t_atoms atoms;
35  t_block blocks[ebNR];
36} t_topology;
37typedef struct {
38} t_nsborder;
39extern FILE *debug;
40typedef struct {
41} t_nrnb;
42typedef 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;
52typedef struct {
53  int nrx,nry,nrz;
54} t_grid;
55typedef struct {
56} t_commrec;
57enum { 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 };
63typedef 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;
78typedef 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;
84enum { egCOUL, egLJ, egBHAM, egLR, egLJLR, egCOUL14, egLJ14, egNR };
85typedef struct {
86  real *ee[egNR];
87} t_grp_ener;
88typedef struct {
89  t_grp_ener estat;
90} t_groups;
91typedef unsigned long t_excl;
92static 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}
105static void reset_neighbor_list(t_forcerec *fr,int bLR,int eNL)
106{
107    reset_nblist(&(fr->nlist_lr[eNL]));
108}
109static 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}
120static 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}
131static void close_neighbor_list(t_forcerec *fr,int bLR,int eNL)
132{
133    close_nblist(&(fr->nlist_lr[eNL]));
134}
135static 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}
141static 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}
395static 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}
407int 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}
431static 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}
465static 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}
490static 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}
621int 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