1 /* { dg-do compile } */
2 /* { dg-options "-w -O3 -mcpu=cortex-a53" } */
3 typedef struct __sFILE __FILE
;
8 typedef real matrix
[3][3];
10 ebCGS
,ebMOLS
,ebSBLOCKS
,ebNR
13 efepNO
, efepYES
, efepNR
16 esolNO
, esolMNO
, esolWATER
, esolWATERWATER
, esolNR
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
,
65 real rcoulomb_switch
,rcoulomb
;
66 real rvdw_switch
,rvdw
;
72 t_nblist nlist_sr
[eNL_NR
];
73 t_nblist nlist_lr
[eNL_NR
];
79 real
*chargeA
,*chargeB
,*chargeT
;
82 unsigned short *cTC
,*cENER
,*cACC
,*cFREEZE
,*cXTC
,*cVCM
;
84 enum { egCOUL
, egLJ
, egBHAM
, egLR
, egLJLR
, egCOUL14
, egLJ14
, egNR
};
91 typedef unsigned long t_excl
;
92 static void reset_nblist(t_nblist
*nl
)
105 static void reset_neighbor_list(t_forcerec
*fr
,int bLR
,int eNL
)
107 reset_nblist(&(fr
->nlist_lr
[eNL
]));
109 static void close_i_nblist(t_nblist
*nlist
)
111 int nri
= nlist
->nri
;
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
)
120 static void close_nblist(t_nblist
*nlist
)
122 if (nlist
->maxnri
> 0) {
123 int nri
= nlist
->nri
;
124 if ((nlist
->jindex
[nri
+1] > nlist
->jindex
[nri
]) &&
125 (nlist
->gid
[nri
] != -1)) {
127 nlist
->jindex
[nri
+2] = nlist
->nrj
;
131 static void close_neighbor_list(t_forcerec
*fr
,int bLR
,int eNL
)
133 close_nblist(&(fr
->nlist_lr
[eNL
]));
135 static void add_j_to_nblist(t_nblist
*nlist
,atom_id j_atom
)
138 nlist
->jjnr
[nrj
] = j_atom
;
141 static void put_in_list(int bHaveLJ
[],
142 int ngid
,t_mdatoms
*md
,
143 int icg
,int jgid
,int nj
,atom_id jjcg
[],
145 t_excl bExcl
[],int shift
,
146 t_forcerec
*fr
,int bLR
,
147 int bVDWOnly
,int bCoulOnly
)
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
;
156 unsigned short *cENER
;
157 real
*charge
,*chargeB
;
159 int bWater
,bMNO
,bFree
,bFreeJ
,bNotEx
,*bPert
;
160 charge
= md
->chargeA
;
161 chargeB
= md
->chargeB
;
165 bPert
= md
->bPerturbed
;
167 nicg
= index
[icg
+1]-i0
;
168 bMNO
= (fr
->solvent_type
[icg
] == esolMNO
);
171 vdw
= &fr
->nlist_lr
[eNL_VDW
];
172 coul
= &fr
->nlist_lr
[eNL_QQ_WATER
];
173 vdwc_ww
= &fr
->nlist_lr
[eNL_VDWQQ_WATERWATER
];
175 vdwc
= &fr
->nlist_lr
[eNL_VDWQQ_SOLMNO
];
177 if (fr
->efep
!= efepNO
) {
178 vdw_free
= &fr
->nlist_lr
[eNL_VDW_FREE
];
179 coul_free
= &fr
->nlist_lr
[eNL_QQ_FREE
];
185 vdwc
= &fr
->nlist_sr
[eNL_VDWQQ_SOLMNO
];
187 if (fr
->efep
!= efepNO
) {
188 vdwc_free
= &fr
->nlist_sr
[eNL_VDWQQ_FREE
];
191 if (fr
->efep
==efepNO
) {
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));
200 new_i_nblist(vdw
,bLR
? F_LJLR
: F_LJ
,i_atom
,shift
,gid
,((void *)0));
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));
205 for(j
=0; (j
<nj
); j
++) {
209 if (bWater
&& (fr
->solvent_type
[jcg
] == esolWATER
)) {
211 add_j_to_nblist(vdw
,jj0
);
213 add_j_to_nblist(coul_ww
,jj0
);
214 add_j_to_nblist(vdwc_ww
,jj0
);
219 for(jj
=jj0
; (jj
<jj1
); jj
++) {
220 if (fabs(charge
[jj
]) > 1.2e-38)
221 add_j_to_nblist(coul
,jj
);
223 } else if (bVDWOnly
) {
224 for(jj
=jj0
; (jj
<jj1
); jj
++)
225 if (bHaveLJ
[type
[jj
]])
226 add_j_to_nblist(vdw
,jj
);
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
);
240 close_i_nblist(coul
);
241 close_i_nblist(vdwc
);
242 close_i_nblist(coul_ww
);
243 close_i_nblist(vdwc_ww
);
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]));
251 new_i_nblist(vdw
,bLR
? F_LR
: F_SR
,i_atom
,shift
,gid
,
252 &(fr
->mno_index
[icg
*3]));
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
++) {
261 for(jj
=jj0
; (jj
<jj1
); jj
++) {
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
);
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
);
278 close_i_nblist(coul
);
279 close_i_nblist(vdwc
);
282 for(i
=0; i
<nicg
; i
++) {
283 igid
= cENER
[i_atom
];
284 gid
= ((igid
< jgid
) ? (igid
*ngid
+jgid
) : (jgid
*ngid
+igid
));
286 if (!bCoulOnly
&& !bVDWOnly
)
287 new_i_nblist(vdwc
,bLR
? F_LJLR
: F_LJ
,i_atom
,shift
,gid
,((void *)0));
289 new_i_nblist(vdw
,bLR
? F_LR
: F_SR
,i_atom
,shift
,gid
,((void *)0));
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
++) {
300 for(jj
=jj0
; jj
<jj1
; jj
++) {
301 bNotEx
= !((int) ((bExcl
)[((atom_id
) (jj
))] & (1<<((atom_id
) (i
)))));
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
);
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
);
322 close_i_nblist(coul
);
323 close_i_nblist(vdwc
);
327 for(i
=0; i
<nicg
; i
++) {
328 igid
= cENER
[i_atom
];
329 gid
= ((igid
< jgid
) ? (igid
*ngid
+jgid
) : (jgid
*ngid
+igid
));
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));
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
++) {
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
)))));
360 add_j_to_nblist(coul_free
,jj
);
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
);
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
);
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
);
389 close_i_nblist(vdw_free
);
390 close_i_nblist(coul_free
);
391 close_i_nblist(vdwc_free
);
395 static void setexcl(atom_id start
,atom_id end
,t_block
*excl
,int 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
)));
407 int calc_naaj(int icg
,int cgtot
)
410 if ((cgtot
% 2) == 1) {
413 else if ((cgtot
% 4) == 0) {
431 static void get_dx(int Nx
,real gridx
,real grid_x
,real rc2
,real x
,
432 int *dx0
,int *dx1
,real
*dcx2
)
436 xgi
= (int)(Nx
+x
*grid_x
)-Nx
;
440 } else if (xgi
>= Nx
) {
450 for(i
=xgi0
; i
>=0; i
--) {
457 for(i
=xgi1
; i
<Nx
; i
++) {
465 static void do_longrange(FILE *log
,t_commrec
*cr
,t_topology
*top
,t_forcerec
*fr
,
466 int ngid
,t_mdatoms
*md
,int icg
,
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
[])
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
);
485 put_in_list(bHaveLJ
,ngid
,md
,icg
,jgid
,nlr
,lr
,top
->blocks
[ebCGS
].index
,
487 1,bVDWOnly
,bCoulOnly
);
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
,
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
;
510 real r2
,rs2
,rvdw2
,rcoul2
,rm2
,rl2
,XI
,YI
,ZI
,dcx
,dcy
,dcz
,tmp1
,tmp2
;
512 int use_twinrange
,use_two_cutoffs
;
514 rs2
= ((fr
->rlist
)*(fr
->rlist
));
515 if (fr
->bTwinRange
) {
516 rvdw2
= ((fr
->rvdw
)*(fr
->rvdw
));
517 rcoul2
= ((fr
->rcoulomb
)*(fr
->rcoulomb
));
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
)));
532 (nl_lr_ljc
)=save_calloc("nl_lr_ljc","ns.c",1349, (ngid
),sizeof(*(nl_lr_ljc
)));
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
])));
538 (nl_lr_ljc
[j
])=save_calloc("nl_lr_ljc[j]","ns.c",1358, (1024),sizeof(*(nl_lr_ljc
[j
])));
540 (nl_lr_one
[j
])=save_calloc("nl_lr_one[j]","ns.c",1360, (1024),sizeof(*(nl_lr_one
[j
])));
543 fprintf(debug
,"ns5_core: rs2 = %g, rvdw2 = %g, rcoul2 = %g (nm^2)\n",
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
)));
554 gridx
= box
[0][0]/grid
->nrx
;
555 gridy
= box
[1][1]/grid
->nry
;
556 gridz
= box
[2][2]/grid
->nrz
;
560 for(iicg
=fr
->cg0
; (iicg
< fr
->hcg
); iicg
++) {
561 icg
= cg_index
[iicg
];
563 fatal_error(0,"icg = %d, iicg = %d, file %s, line %d",icg
,iicg
,"ns.c",
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
);
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
);
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
))) {
588 if (!i_eg_excl
[jgid
]) {
590 if (nsr
[jgid
] >= 1024) {
591 put_in_list(bHaveLJ
,ngid
,md
,icg
,jgid
,
592 nsr
[jgid
],nl_sr
[jgid
],
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
,
601 nl_lr_one
[jgid
],bexcl
,shift
,x
,
603 lambda
,dvdlambda
,grps
,
604 bVDWOnly
,bCoulOnly
,0,
621 int search_neighbors(FILE *log
,t_forcerec
*fr
,
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
)
628 static t_grid
*grid
=((void *)0);
629 static t_excl
*bexcl
;
631 static int *cg_index
=((void *)0),*slab_index
=((void *)0);
632 static int *bExcludeAlleg
;
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
);