?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
changed reading hint
[gromacs/adressmacs.git] / src / mdlib / ns.c
blob0c67d8a5d658322f6070ce7e6cd5b76c4874899d
1 /*
2 * $Id$
3 *
4 * This source code is part of
5 *
6 * G R O M A C S
7 *
8 * GROningen MAchine for Chemical Simulations
9 *
10 * VERSION 2.0
11 *
12 * Copyright (c) 1991-1999
13 * BIOSON Research Institute, Dept. of Biophysical Chemistry
14 * University of Groningen, The Netherlands
15 *
16 * Please refer to:
17 * GROMACS: A message-passing parallel molecular dynamics implementation
18 * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
19 * Comp. Phys. Comm. 91, 43-56 (1995)
20 *
21 * Also check out our WWW page:
22 * http://md.chem.rug.nl/~gmx
23 * or e-mail to:
24 * gromacs@chem.rug.nl
25 *
26 * And Hey:
27 * GRowing Old MAkes el Chrono Sweat
28 */
29 static char *SRCID_ns_c = "$Id$";
30
31 #include <math.h>
32 #include <string.h>
33 #include "sysstuff.h"
34 #include "assert.h"
35 #include "smalloc.h"
36 #include "macros.h"
37 #include "maths.h"
38 #include "vec.h"
39 #include "network.h"
40 #include "nsgrid.h"
41 #include "force.h"
42 #include "ns.h"
43 #include "pbc.h"
44 #include "names.h"
45 #include "fatal.h"
46 #include "nrnb.h"
47 #include "ns.h"
48 #include "force.h"
49 #include "txtdump.h"
50
51 #define MAX_CG 1024
52
53 typedef struct {
54 int ncg;
55 int nj;
56 atom_id jcg[MAX_CG];
57 } t_ns_buf;
58
59 /*
60 * E X C L U S I O N H A N D L I N G
61 */
62 typedef unsigned long t_excl;
63
64 #ifdef DEBUG
65 static void SETEXCL_(t_excl e[],atom_id i,atom_id j)
66 { e[j] = e[j] | (1<<i); }
67 static void RMEXCL_(t_excl e[],atom_id i,atom_id j)
68 { e[j]=e[j] & ~(1<<i); }
69 static bool ISEXCL_(t_excl e[],atom_id i,atom_id j)
70 { return (bool)(e[j] & (1<<i)); }
71 static bool NOTEXCL_(t_excl e[],atom_id i,atom_id j)
72 { return !(ISEXCL(e,i,j)); }
73 #else
74 #define SETEXCL(e,i,j) (e)[((atom_id) j)] |= (1<<((atom_id) i))
75 #define RMEXCL(e,i,j) (e)[((atom_id) j)] &= (~(1<<((atom_id) i)))
76 #define ISEXCL(e,i,j) (bool) ((e)[((atom_id) j)] & (1<<((atom_id) i)))
77 #define NOTEXCL(e,i,j) !(ISEXCL(e,i,j))
78 #endif
79
80 /************************************************
81 *
82 * U T I L I T I E S F O R N S
83 *
84 ************************************************/
85
86 static int NLI_INC = 1000;
87 static int NLJ_INC = 16384;
88
89 static void reallocate_nblist(t_nblist *nl)
90 {
91 if (debug)
92 fprintf(debug,"reallocating neigborlist il_code=%d, maxnri=%d\n",
93 nl->il_code,nl->maxnri);
94 srenew(nl->iinr, nl->maxnri+2);
95 srenew(nl->gid, nl->maxnri+2);
96 srenew(nl->shift, nl->maxnri+2);
97 srenew(nl->jindex, nl->maxnri+2);
98 }
99
100 static void init_nblist(t_nblist *nl,int homenr,int il_code)
101 {
102 nl->il_code = il_code;
103 /* maxnri is influenced by the number of shifts (maximum is 8)
104 * and the number of energy groups.
105 * If it is not enough, nl memory will be reallocated during the run.
106 * 4 seems to be a reasonable factor, which only causes reallocation
107 * during runs with tiny and many energygroups.
108 */
109 nl->maxnri = homenr*4;
110 nl->maxnrj = 0;
111 nl->nri = 0;
112 nl->nrj = 0;
113 nl->iinr = NULL;
114 nl->gid = NULL;
115 nl->shift = NULL;
116 nl->jindex = NULL;
117 reallocate_nblist(nl);
118 nl->jindex[0] = 0;
119 nl->jindex[1] = 0;
120 if (nl->maxnri > 0)
121 nl->iinr[0] = -1;
122 }
123
124 static unsigned int nbf_index(bool bCoul,bool bRF,bool bBham,
125 bool bTab,bool bWater,bool bEwald)
126 {
127 /* lot of redundancy here,
128 * since we cant have RF and EWALD simultaneously...
129 */
130
131 int inloop[64] = {
132 eNR_LJC, eNR_QQ, eNR_BHAM, eNR_QQ,
133 eNR_LJCRF, eNR_QQRF, eNR_BHAMRF, eNR_QQRF,
134 eNR_TAB, eNR_COULTAB, eNR_BHAMTAB, eNR_COULTAB,
135 eNR_TAB, eNR_COULTAB, eNR_BHAMTAB, eNR_COULTAB,
136 eNR_LJC_WAT, eNR_QQ_WAT, eNR_BHAM_WAT, eNR_QQ_WAT,
137 eNR_LJCRF_WAT, eNR_QQRF_WAT, eNR_BHAMRF_WAT, eNR_QQRF_WAT,
138 eNR_TAB_WAT, eNR_COULTAB_WAT,eNR_BHAMTAB_WAT,eNR_COULTAB_WAT,
139 eNR_TAB_WAT, eNR_COULTAB_WAT,eNR_BHAMTAB_WAT,eNR_COULTAB_WAT,
140 eNR_LJC_EW, eNR_QQ_EW, eNR_BHAM_EW, eNR_QQ_EW,
141 eNR_LJC_EW, eNR_QQ_EW, eNR_BHAM_EW, eNR_QQ_EW,
142 eNR_TAB, eNR_COULTAB, eNR_BHAMTAB, eNR_COULTAB,
143 eNR_TAB, eNR_COULTAB, eNR_BHAMTAB, eNR_COULTAB,
144 eNR_LJC_WAT_EW, eNR_QQ_WAT_EW, eNR_BHAM_WAT_EW,eNR_QQ_WAT_EW,
145 eNR_LJC_WAT_EW, eNR_QQ_WAT_EW, eNR_BHAM_WAT_EW,eNR_QQ_WAT_EW,
146 eNR_TAB_WAT, eNR_COULTAB_WAT,eNR_BHAMTAB_WAT,eNR_COULTAB_WAT,
147 eNR_TAB_WAT, eNR_COULTAB_WAT,eNR_BHAMTAB_WAT,eNR_COULTAB_WAT
148 };
149
150 unsigned int ni;
151
152 ni = bCoul | (bBham << 1) | (bRF << 2) | (bTab << 3) | (bWater << 4)
153 | (bEwald << 5);
154
155 return inloop[ni];
156 }
157
158 void init_neighbor_list(FILE *log,t_forcerec *fr,int homenr)
159 {
160 /* Make maxlr tunable! (does not seem to be a big difference though)
161 * This parameter determines the number of i particles in a long range
162 * neighbourlist. Too few means many function calls, too many means
163 * cache trashing.
164 */
165 int maxsr,maxsr_wat,maxlr,maxlr_wat;
166
167 maxsr = homenr-fr->nWatMol*3;
168 if (maxsr < 0)
169 fatal_error(0,"%s, %d: Negative number of short range atoms.\n"
170 "Call your Gromacs dealer for assistance.",__FILE__,__LINE__);
171 maxsr_wat = fr->nWatMol;
172 maxlr = max(maxsr/10,50);
173 maxlr_wat = max(maxsr_wat/10,50);
174
175 init_nblist(&fr->nlist_sr[eNL_VDW],maxsr,
176 nbf_index(FALSE,fr->bRF,fr->bBHAM,fr->bTab,FALSE,fr->bEwald));
177 init_nblist(&fr->nlist_sr[eNL_QQ],maxsr,
178 nbf_index(TRUE,fr->bRF,fr->bBHAM,fr->bTab,FALSE,fr->bEwald));
179 if (fr->bPert)
180 init_nblist(&fr->nlist_sr[eNL_FREE],maxsr,
181 fr->bBHAM ? eNR_BHAM_FREE : eNR_LJC_FREE);
182 if (fr->bWaterOpt) {
183 init_nblist(&fr->nlist_sr[eNL_VDW_WAT],maxsr_wat,
184 nbf_index(FALSE,fr->bRF,fr->bBHAM,fr->bTab,TRUE,fr->bEwald));
185 init_nblist(&fr->nlist_sr[eNL_QQ_WAT],maxsr_wat,
186 nbf_index(TRUE,fr->bRF,fr->bBHAM,fr->bTab,TRUE,fr->bEwald));
187 }
188 if (fr->bTwinRange) {
189 fprintf(log,"Allocating space for long range neighbour list of %d atoms\n",
190 maxlr);
191 init_nblist(&fr->nlist_lr[eNL_VDW],maxlr,
192 nbf_index(FALSE,fr->bRF,fr->bBHAM,fr->bTab,FALSE,fr->bEwald));
193 init_nblist(&fr->nlist_lr[eNL_QQ],maxlr,
194 nbf_index(TRUE,fr->bRF,fr->bBHAM,fr->bTab,FALSE,fr->bEwald));
195 if (fr->bPert)
196 init_nblist(&fr->nlist_lr[eNL_FREE],maxlr,
197 fr->bBHAM ? eNR_BHAM_FREE : eNR_LJC_FREE);
198 if (fr->bWaterOpt) {
199 init_nblist(&fr->nlist_lr[eNL_VDW_WAT],maxlr_wat,
200 nbf_index(FALSE,fr->bRF,fr->bBHAM,fr->bTab,TRUE,fr->bEwald));
201 init_nblist(&fr->nlist_lr[eNL_QQ_WAT],maxlr_wat,
202 nbf_index(TRUE,fr->bRF,fr->bBHAM,fr->bTab,TRUE,fr->bEwald));
203 }
204 }
205 }
206
207 static void reset_nblist(t_nblist *nl)
208 {
209 nl->nri = 0;
210 nl->nrj = 0;
211 if (nl->maxnri > 0) {
212 nl->iinr[0] = -1;
213 if (nl->maxnrj > 1) {
214 nl->jindex[0] = 0;
215 nl->jindex[1] = 0;
216 }
217 }
218 }
219
220 static void reset_neighbor_list(t_forcerec *fr,bool bLR,int eNL)
221 {
222 int i;
223
224 if (bLR)
225 reset_nblist(&(fr->nlist_lr[eNL]));
226 else {
227 for(i=0; (i<eNL_NR); i++)
228 reset_nblist(&(fr->nlist_sr[i]));
229 }
230 }
231
232 static gmx_inline void new_i_nblist(t_nblist *nlist,
233 int ftype,int i_atom,int shift,int gid)
234 {
235 int i,nri,nshift;
236
237 if (nlist->maxnrj <= nlist->nrj + NLJ_INC-1) {
238 if (debug)
239 fprintf(debug,"Adding %5d J particles for nblist %s\n",NLJ_INC,
240 interaction_function[ftype].longname);
241
242 nlist->maxnrj += NLJ_INC;
243 srenew(nlist->jjnr,nlist->maxnrj);
244 }
245
246 nri = nlist->nri;
247
248 /* Check whether we have to increase the i counter */
249 if ((nlist->iinr[nri] != i_atom) ||
250 (nlist->shift[nri] != shift) ||
251 (nlist->gid[nri] != gid)) {
252 /* This is something else. Now see if any entries have
253 * been added in the list of the previous atom.
254 */
255 if ((nlist->jindex[nri+1] > nlist->jindex[nri]) &&
256 (nlist->iinr[nri] != -1)) {
257
258 /* If so increase the counter */
259 nlist->nri++;
260 nri++;
261 if (nlist->nri >= nlist->maxnri) {
262 nlist->maxnri += NLI_INC;
263 reallocate_nblist(nlist);
264 }
265 }
266 /* Set the number of neighbours and the atom number */
267 nlist->jindex[nri+1] = nlist->jindex[nri];
268 nlist->iinr[nri] = i_atom;
269 nlist->gid[nri] = gid;
270 nlist->shift[nri] = shift;
271 }
272 }
273
274 #ifdef SORTNLIST
275 #define aswap(v,i,j) { \
276 atom_id temp; \
277 \
278 temp=v[i]; \
279 v[i]=v[j]; \
280 v[j]=temp; \
281 }
282
283 static void quicksort(atom_id v[], int left, int right)
284 {
285 int i,last;
286
287 if (left >= right) /* Do nothing if array contains */
288 return; /* fewer than two elements */
289 aswap(v,left,(left+right)/2); /* Move partition element */
290 last=left; /* to v[0] */
291 for(i=left+1; (i<=right); i++) /* partition */
292 if (v[i] < v[left]) {
293 last++;
294 aswap(v,last,i); /* watch out for macro trick */
295 }
296 aswap(v,left,last); /* restore partition element */
297 quicksort(v,left,last-1);
298 quicksort(v,last+1,right);
299 }
300 #endif
301
302 static gmx_inline void close_i_nblist(t_nblist *nlist)
303 {
304 int nri = nlist->nri;
305
306 nlist->jindex[nri+1] = nlist->nrj;
307 }
308
309 static gmx_inline void close_nblist(t_nblist *nlist)
310 {
311 if (nlist->maxnri > 0) {
312 int nri = nlist->nri;
313
314 if ((nlist->jindex[nri+1] > nlist->jindex[nri]) &&
315 (nlist->iinr[nri] != -1)) {
316 nlist->nri++;
317 nlist->jindex[nri+2] = nlist->nrj;
318 }
319 }
320 }
321
322 static gmx_inline void close_neighbor_list(t_forcerec *fr,bool bLR,int eNL)
323 {
324 int i;
325
326 if (bLR)
327 close_nblist(&(fr->nlist_lr[eNL]));
328 else {
329 for(i=0; (i<eNL_NR); i++)
330 close_nblist(&(fr->nlist_sr[i]));
331 }
332 }
333
334 static void add_j_to_nblist(t_nblist *nlist,int j_atom)
335 {
336 int nrj=nlist->nrj;
337
338 nlist->jjnr[nrj] = j_atom;
339 nlist->nrj ++;
340 }
341
342 static gmx_inline void put_in_list(bool bHaveLJ[],
343 int nWater,int ngid,t_mdatoms *md,
344 int icg,int jgid,int nj,atom_id jjcg[],
345 atom_id index[],atom_id a[],
346 t_excl bExcl[],int shift,
347 t_forcerec *fr,bool bLR,bool bCoulOnly)
348 {
349 t_nblist *vdw,*coul,*free=NULL;
350
351 int i,j,jcg,igid,gid,ind_ij;
352 atom_id jj,jj0,jj1,i_atom,j_atom;
353 int i0,nicg;
354
355 int *type;
356 unsigned short *cENER;
357 real *charge;
358 real qi,qq,rlj;
359 bool bWater,bFree,bFreeJ,bNotEx,*bPert;
360
361 #ifdef SORTNLIST
362 /* Quicksort the charge groups in the neighbourlist to obtain
363 * better caching properties. We do this only for the short range,
364 * i.e. when we use the nlist more than once
365 */
366 if (!bLR)
367 quicksort(jjcg,0,nj-1);
368 #endif
369 /* Copy some pointers */
370 charge = md->chargeA;
371 type = md->typeA;
372 cENER = md->cENER;
373 bPert = md->bPerturbed;
374
375 /* Check whether this molecule is a water molecule */
376 i0 = index[icg];
377 nicg = index[icg+1]-i0;
378 bWater = (((type[a[i0]] == nWater) && (nicg == 3)) &&
379 (!bPert[a[i0]]) && (!bPert[a[i0+1]]) && (!bPert[a[i0+2]]));
380 if (bWater)
381 nicg = 1;
382
383 if (bLR) {
384 if (bWater) {
385 vdw = &fr->nlist_lr[eNL_VDW_WAT];
386 coul = &fr->nlist_lr[eNL_QQ_WAT];
387 }
388 else {
389 vdw = &fr->nlist_lr[eNL_VDW];
390 coul = &fr->nlist_lr[eNL_QQ];
391 }
392 if (fr->bPert)
393 free = &fr->nlist_lr[eNL_FREE];
394 }
395 else {
396 if (bWater) {
397 vdw = &fr->nlist_sr[eNL_VDW_WAT];
398 coul = &fr->nlist_sr[eNL_QQ_WAT];
399 }
400 else {
401 vdw = &fr->nlist_sr[eNL_VDW];
402 coul = &fr->nlist_sr[eNL_QQ];
403 }
404 if (fr->bPert)
405 free = &fr->nlist_sr[eNL_FREE];
406 }
407
408 /* Loop over the atoms in the i charge group */
409 for(i=0; (i<nicg); i++) {
410 i_atom = a[i0+i];
411 igid = cENER[i_atom];
412 gid = GID(igid,jgid,ngid);
413
414 /* Create new i_atom for each energy group */
415 if (!bCoulOnly)
416 new_i_nblist(vdw,bLR ? F_LJLR : F_LJ,i_atom,shift,gid);
417 new_i_nblist(coul,bLR ? F_LR : F_SR,i_atom,shift,gid);
418 if (fr->bPert)
419 new_i_nblist(free,F_DVDL,i_atom,shift,gid);
420
421 /* Loop over the j charge groups */
422 for(j=0; (j<nj); j++) {
423 jcg=jjcg[j];
424
425 if (bWater && (jcg==icg))
426 continue;
427
428 /* Check for large charge groups */
429 if (jcg == icg)
430 jj0 = i0 + i + 1;
431 else
432 jj0 = index[jcg];
433
434 jj1=index[jcg+1];
435
436 if (bWater && !fr->bPert) {
437 if (bCoulOnly) {
438 for(jj=jj0; (jj<jj1); jj++) {
439 j_atom = a[jj];
440
441 add_j_to_nblist(coul,j_atom);
442 }
443 }
444 else {
445 for(jj=jj0; (jj<jj1); jj++) {
446 j_atom = a[jj];
447
448 if (bHaveLJ[type[j_atom]])
449 add_j_to_nblist(vdw,j_atom);
450 else if (charge[j_atom] != 0)
451 add_j_to_nblist(coul,j_atom);
452 }
453 }
454 }
455 else {
456 /* Finally loop over the atoms in the j-charge group */
457 bFree = bPert[i_atom];
458 qi = charge[i_atom];
459 for(jj=jj0; (jj<jj1); jj++) {
460 j_atom = a[jj];
461 bFreeJ = bFree || bPert[j_atom];
462 bNotEx = NOTEXCL(bExcl,i,j_atom);
463
464 if (bNotEx) {
465 if (bFreeJ)
466 add_j_to_nblist(free,j_atom);
467 else if (bCoulOnly)
468 /* This is done whether or not bWater is set */
469 add_j_to_nblist(coul,j_atom);
470 else {
471 if (bHaveLJ[type[j_atom]])
472 add_j_to_nblist(vdw,j_atom);
473 else if (qi*charge[j_atom] != 0)
474 add_j_to_nblist(coul,j_atom);
475 }
476 }
477 }
478 }
479 }
480 close_i_nblist(coul);
481 close_i_nblist(vdw);
482 if (fr->bPert)
483 close_i_nblist(free);
484 }
485 }
486
487 void setexcl(int nri,atom_id ia[],t_block *excl,bool b,
488 t_excl bexcl[])
489 {
490 atom_id k;
491 int i,inr;
492
493 if (b) {
494 for(i=0; (i<nri); i++) {
495 inr = ia[i];
496 for(k=excl->index[inr]; (k<excl->index[inr+1]); k++) {
497 SETEXCL(bexcl,i,excl->a[k]);
498 }
499 }
500 }
501 else {
502 for(i=0; (i<nri); i++) {
503 inr = ia[i];
504 for(k=excl->index[inr]; (k<excl->index[inr+1]); k++) {
505 RMEXCL(bexcl,i,excl->a[k]);
506 }
507 }
508 }
509 }
510
511 int calc_naaj(int icg,int cgtot)
512 {
513 int naaj;
514
515 if ((cgtot % 2) == 1) {
516 /* Odd number of charge groups, easy */
517 naaj = 1+(cgtot/2);
518 }
519 else if ((cgtot % 4) == 0) {
520 /* Multiple of four is hard */
521 if (icg < cgtot/2) {
522 if ((icg % 2) == 0)
523 naaj=1+(cgtot/2);
524 else
525 naaj=cgtot/2;
526 }
527 else {
528 if ((icg % 2) == 1)
529 naaj=1+(cgtot/2);
530 else
531 naaj=cgtot/2;
532 }
533 }
534 else {
535 /* cgtot/2 = odd */
536 if ((icg % 2) == 0)
537 naaj=1+(cgtot/2);
538 else
539 naaj=cgtot/2;
540 }
541 #ifdef DEBUG
542 fprintf(log,"naaj=%d\n",naaj);
543 #endif
544 return naaj;
545 }
546
547 /************************************************
548 *
549 * S I M P L E C O R E S T U F F
550 *
551 ************************************************/
552
553 static real calc_image_rect(rvec xi,rvec xj,rvec box_size,
554 rvec b_inv,int *shift)
555 {
556 const real h15=1.5;
557 real ddx,ddy,ddz;
558 real dx,dy,dz;
559 real r2;
560 int tx,ty,tz;
561
562 /* Compute diff vector */
563 dx=xj[XX]-xi[XX];
564 dy=xj[YY]-xi[YY];
565 dz=xj[ZZ]-xi[ZZ];
566
567 /* Perform NINT operation, using trunc operation, therefore
568 * we first add 1.5 then subtract 1 again
569 */
570 tx=dx*b_inv[XX]+h15;
571 ty=dy*b_inv[YY]+h15;
572 tz=dz*b_inv[ZZ]+h15;
573 tx--;
574 ty--;
575 tz--;
576
577 /* Correct diff vector for translation */
578 ddx=tx*box_size[XX]-dx;
579 ddy=ty*box_size[YY]-dy;
580 ddz=tz*box_size[ZZ]-dz;
581
582 /* Distance squared */
583 r2=(ddx*ddx)+(ddy*ddy)+(ddz*ddz);
584
585 *shift=XYZ2IS(tx,ty,tz);
586
587 return r2;
588 }
589
590 static void ns_inner_rect(rvec x[],int icg,int njcg,atom_id jcg[],
591 bool bBox,rvec box_size,rvec b_inv,real rcut2,
592 t_block *cgs,t_ns_buf **ns_buf,unsigned short gid[])
593 {
594 int shift;
595 int j,nrj,jgid;
596 atom_id cg_j,*cgindex,*cga;
597 t_ns_buf *nsbuf;
598
599 cgindex = cgs->index;
600 cga = cgs->a;
601 if (bBox) {
602 shift = CENTRAL;
603 for(j=0; (j<njcg); j++) {
604 cg_j = jcg[j];
605 nrj = cgindex[cg_j+1]-cgindex[cg_j];
606 if (calc_image_rect(x[icg],x[cg_j],box_size,b_inv,&shift) < rcut2) {
607 jgid = gid[cga[cgindex[cg_j]]];
608 nsbuf = &ns_buf[jgid][shift];
609 if (nsbuf->ncg >= MAX_CG)
610 fatal_error(0,"Too many charge groups (%d) in buffer",nsbuf->ncg);
611 nsbuf->jcg[nsbuf->ncg++]=cg_j;
612 nsbuf->nj += nrj;
613 }
614 }
615 } else {
616 for(j=0; (j<njcg); j++) {
617 cg_j = jcg[j];
618 nrj = cgindex[cg_j+1]-cgindex[cg_j];
619 if ((rcut2 == 0) || (distance2(x[icg],x[cg_j]) < rcut2)) {
620 jgid = gid[cga[cgindex[cg_j]]];
621 nsbuf = &ns_buf[jgid][CENTRAL];
622 if (nsbuf->ncg >= MAX_CG)
623 fatal_error(0,"Too many charge groups (%d) in buffer",nsbuf->ncg);
624 nsbuf->jcg[nsbuf->ncg++]=cg_j;
625 nsbuf->nj += nrj;
626 }
627 }
628 }
629 }
630
631 static int ns_simple_core(t_forcerec *fr,
632 t_topology *top,
633 t_mdatoms *md,
634 rvec box_size,t_excl bexcl[],
635 int ngid,t_ns_buf **ns_buf,
636 bool bHaveLJ[])
637 {
638 static atom_id *aaj=NULL;
639 int naaj,k;
640 real rlist2;
641 int nsearch,icg,jcg,i0,nri,nn;
642 t_ns_buf *nsbuf;
643 atom_id *i_atoms;
644 t_block *cgs=&(top->blocks[ebCGS]);
645 t_block *excl=&(top->atoms.excl);
646 rvec b_inv;
647 int m;
648 bool bBox;
649
650 if (aaj==NULL) {
651 snew(aaj,2*cgs->nr);
652 for(jcg=0; (jcg<cgs->nr); jcg++) {
653 aaj[jcg]=jcg;
654 aaj[jcg+cgs->nr]=jcg;
655 }
656 }
657 rlist2 = sqr(fr->rlist);
658
659 bBox = (fr->eBox != ebtNONE);
660 if (bBox)
661 for(m=0; (m<DIM); m++)
662 b_inv[m]=divide(1.0,box_size[m]);
663
664 nsearch=0;
665 for (icg=fr->cg0; (icg<fr->hcg); icg++) {
666 i0 = cgs->index[icg];
667 nri = cgs->index[icg+1]-i0;
668 i_atoms = &(cgs->a[i0]);
669 setexcl(nri,i_atoms,excl,TRUE,bexcl);
670
671 naaj=calc_naaj(icg,cgs->nr);
672 ns_inner_rect(fr->cg_cm,icg,naaj,&(aaj[icg]),
673 bBox,box_size,b_inv,rlist2,cgs,ns_buf,md->cENER);
674 nsearch += naaj;
675
676 for(nn=0; (nn<ngid); nn++) {
677 for(k=0; (k<SHIFTS); k++) {
678 nsbuf = &(ns_buf[nn][k]);
679 if (nsbuf->ncg > 0) {
680 put_in_list(bHaveLJ,fr->nWater,
681 ngid,md,icg,nn,nsbuf->ncg,nsbuf->jcg,
682 cgs->index,cgs->a,bexcl,k,fr,FALSE,FALSE);
683 nsbuf->ncg=nsbuf->nj=0;
684 }
685 }
686 }
687 setexcl(nri,i_atoms,excl,FALSE,bexcl);
688 }
689 close_neighbor_list(fr,FALSE,-1);
690
691 return nsearch;
692 }
693
694 /************************************************
695 *
696 * N S 5 G R I D S T U F F
697 *
698 ************************************************/
699 static bool get_dx(int cx,int Nx,int tx,int delta,int *dx0,int *dx1)
700 {
701 if (tx == 1) {
702 if (cx >= delta)
703 return TRUE;
704 else {
705 *dx0=Nx+cx-delta;
706 *dx1=Nx-1;
707 }
708 } else if (tx == 0) {
709 if (cx < delta)
710 *dx0=0;
711 else
712 *dx0=cx-delta;
713 if (cx < Nx-delta)
714 *dx1=cx+delta;
715 else
716 *dx1=Nx-1;
717 }
718 else {
719 if (cx < Nx-delta)
720 return TRUE;
721 else {
722 *dx0=0;
723 *dx1=cx+delta-Nx;
724 }
725 }
726
727 return FALSE;
728 }
729
730 #define sqr(x) ((x)*(x))
731 #define calc_dx2(XI,YI,ZI,y) (sqr(XI-y[XX])+sqr(YI-y[YY])+sqr(ZI-y[ZZ]))
732 #define calc_cyl_dx2(XI,YI,y) (sqr(XI-y[XX])+sqr(YI-y[YY]))
733 /****************************************************
734 *
735 * F A S T N E I G H B O R S E A R C H I N G
736 *
737 * Optimized neighboursearching routine using grid
738 * of at least 5x5x5, see GROMACS manual
739 *
740 ****************************************************/
741
742 static void do_longrange(FILE *log,t_topology *top,t_forcerec *fr,
743 int ngid,t_mdatoms *md,int icg,
744 int jgid,int nlr,
745 atom_id lr[],t_excl bexcl[],int shift,
746 rvec x[],rvec box_size,t_nrnb *nrnb,
747 real lambda,real *dvdlambda,
748 t_groups *grps,bool bCoulOnly,
749 bool bDoForces,bool bHaveLJ[])
750 {
751 int i;
752
753 for(i=0; (i<eNL_NR); i++) {
754 if ((fr->nlist_lr[i].nri > fr->nlist_lr[i].maxnri-32) || bDoForces) {
755 close_neighbor_list(fr,TRUE,i);
756
757 /* Evaluate the forces */
758 do_fnbf(log,fr,x,fr->flr,md,
759 grps->estat.ee[egLJLR],grps->estat.ee[egLR],box_size,
760 nrnb,lambda,dvdlambda,TRUE,i);
761
762 reset_neighbor_list(fr,TRUE,i);
763 }
764 }
765
766 if (!bDoForces) {
767 /* Put the long range particles in a list */
768 put_in_list(bHaveLJ,fr->nWater,
769 ngid,md,icg,jgid,nlr,lr,top->blocks[ebCGS].index,
770 top->blocks[ebCGS].a,bexcl,shift,fr,TRUE,bCoulOnly);
771 }
772 }
773
774 static int ns5_core(FILE *log,t_forcerec *fr,int cg_index[],
775 rvec box_size,int ngid,
776 t_topology *top,t_groups *grps,
777 t_grid *grid,rvec x[],t_excl bexcl[],
778 t_nrnb *nrnb,t_mdatoms *md,
779 real lambda,real *dvdlambda,
780 bool bHaveLJ[])
781 {
782 static atom_id **nl_lr_ljc,**nl_lr_coul,**nl_sr=NULL;
783 static int *nlr_ljc,*nlr_coul,*nsr;
784
785 t_block *cgs=&(top->blocks[ebCGS]);
786 unsigned short *gid=md->cENER;
787 atom_id *i_atoms,*cgsatoms=cgs->a,*cgsindex=cgs->index;
788 int tx,ty,tz,cx,cy,cz,dx,dy,dz,cj;
789 int dx0,dx1,dy0,dy1,dz0,dz1;
790 int Nx,Ny,Nz,delta,shift=-1,j,nrj,nns,nn=-1;
791 int icg=-1,iicg,cgsnr,i0,nri,naaj,min_icg,icg_naaj,jjcg,cgj0,jgid;
792 int *grida,*gridnra,*gridind;
793 rvec xi,*cgcm,*svec;
794 real r2,rs2,rvdw2,rcoul2,XI,YI,ZI;
795
796 cgsnr = cgs->nr;
797 rs2 = sqr(fr->rlist);
798 rvdw2 = sqr(fr->rvdw);
799 rcoul2 = sqr(fr->rcoulomb);
800
801 if (nl_sr == NULL) {
802 /* Short range buffers */
803 snew(nl_sr,ngid);
804 /* Counters */
805 snew(nsr,ngid);
806 snew(nlr_ljc,ngid);
807 snew(nlr_coul,ngid);
808
809 if (rvdw2 > rs2)
810 /* Long range LJC buffers */
811 snew(nl_lr_ljc,ngid);
812
813 if (rcoul2 > rvdw2)
814 /* Long range Coul buffers */
815 snew(nl_lr_coul,ngid);
816
817 for(j=0; (j<ngid); j++) {
818 snew(nl_sr[j],MAX_CG);
819 if (rvdw2 > rs2)
820 snew(nl_lr_ljc[j],MAX_CG);
821 if (rcoul2 > rvdw2)
822 snew(nl_lr_coul[j],MAX_CG);
823 }
824 if (debug)
825 fprintf(debug,"ns5_core: rs2 = %g, rvdw2 = %g, rcoul2 = %g (nm^2)\n",
826 rs2,rvdw2,rcoul2);
827 }
828
829 /* Unpack arrays */
830 cgcm = fr->cg_cm;
831 svec = fr->shift_vec;
832 Nx = grid->nrx;
833 Ny = grid->nry;
834 Nz = grid->nrz;
835 grida = grid->a;
836 gridind = grid->index;
837 gridnra = grid->nra;
838 delta = grid->delta;
839 nns = 0;
840
841 /* Loop over charge groups */
842 for(iicg=fr->cg0; (iicg < fr->hcg); iicg++) {
843 icg = cg_index[iicg];
844 /* Consistency check */
845 if (icg != iicg)
846 fatal_error(0,"icg = %d, iicg = %d, file %s, line %d",icg,iicg,__FILE__,
847 __LINE__);
848 i0 = cgsindex[icg];
849 nri = cgsindex[icg+1]-i0;
850 i_atoms = &(cgsatoms[i0]);
851
852 /* Set the exclusions for the atoms in charge group icg using
853 * a bitmask
854 */
855 setexcl(nri,i_atoms,&top->atoms.excl,TRUE,bexcl);
856
857 /* Compute the number of charge groups that fall within the control
858 * of this one (icg)
859 */
860 naaj = calc_naaj(icg,cgsnr);
861 icg_naaj = icg+naaj;
862 min_icg = icg_naaj-cgsnr;
863
864 /* Changed iicg to icg, DvdS 990115
865 * (but see consistency check above, DvdS 990330)
866 */
867 ci2xyz(grid,icg,&cx,&cy,&cz);
868 #ifdef NS5DB
869 fprintf(log,"icg=%5d, naaj=%5d, cx=%2d, cy=%2d, cz=%2d\n",
870 icg,naaj,cx,cy,cz);
871 #endif
872 /* Loop over shift vectors in three dimensions */
873 for (tx=-1; (tx<=1); tx++) {
874 /* Calculate range of cells in X direction that have the shift tx */
875 if (get_dx(cx,Nx,tx,delta,&dx0,&dx1))
876 continue;
877 for (ty=-1; (ty<=1); ty++) {
878 /* Calculate range of cells in Y direction that have the shift ty */
879 if (get_dx(cy,Ny,ty,delta,&dy0,&dy1))
880 continue;
881 for (tz=-1; (tz<=1); tz++) {
882 /* Calculate range of cells in Z direction that have the shift tz */
883 if (get_dx(cz,Nz,tz,delta,&dz0,&dz1))
884 continue;
885
886 /* Get shift vector */
887 shift=XYZ2IS(tx,ty,tz);
888 #ifdef NS5DB
889 assert(shift >= 0);
890 assert(shift < SHIFTS);
891 #endif
892 XI = cgcm[icg][XX]+svec[shift][XX];
893 YI = cgcm[icg][YY]+svec[shift][YY];
894 ZI = cgcm[icg][ZZ]+svec[shift][ZZ];
895 for(nn=0; (nn<ngid); nn++) {
896 nsr[nn] = 0;
897 nlr_ljc[nn] = 0;
898 nlr_coul[nn] = 0;
899 }
900 #ifdef NS5DB
901 fprintf(log,"shift: %2d, dx0,1: %2d,%2d, dy0,1: %2d,%2d, dz0,1: %2d,%2d\n",
902 shift,dx0,dx1,dy0,dy1,dz0,dz1);
903 fprintf(log,"cgcm: %8.3f %8.3f %8.3f\n",cgcm[icg][XX],
904 cgcm[icg][YY],cgcm[icg][ZZ]);
905 fprintf(log,"xi: %8.3f %8.3f %8.3f\n",XI,YI,ZI);
906 #endif
907 for (dx=dx0; (dx<=dx1); dx++) {
908 for (dy=dy0; (dy<=dy1); dy++) {
909 for (dz=dz0; (dz<=dz1); dz++) {
910 /* Find the grid-cell cj in which possible neighbours are */
911 cj = xyz2ci(Ny,Nz,dx,dy,dz);
912
913 /* Check out how many cgs (nrj) there in this cell */
914 nrj = gridnra[cj];
915
916 /* Find the offset in the cg list */
917 cgj0 = gridind[cj];
918
919 /* Loop over cgs */
920 for (j=0; (j<nrj); j++) {
921 jjcg = grida[cgj0+j];
922
923 /* check whether this guy is in range! */
924 if (((jjcg >= icg) && (jjcg < icg_naaj)) ||
925 ((jjcg < min_icg))) {
926 r2=calc_dx2(XI,YI,ZI,cgcm[jjcg]);
927 if (r2 < rcoul2) {
928 jgid = gid[cgsatoms[cgsindex[jjcg]]];
929 if (r2 < rs2) {
930 if (nsr[jgid] >= MAX_CG) {
931 put_in_list(bHaveLJ,fr->nWater,
932 ngid,md,icg,jgid,nsr[jgid],nl_sr[jgid],
933 cgsindex,cgsatoms,bexcl,
934 shift,fr,FALSE,FALSE);
935 nsr[jgid]=0;
936 }
937 nl_sr[jgid][nsr[jgid]++]=jjcg;
938 }
939 else if (r2 < rvdw2) {
940 if (nlr_ljc[jgid] >= MAX_CG) {
941 do_longrange(log,top,fr,ngid,md,icg,jgid,
942 nlr_ljc[jgid],
943 nl_lr_ljc[jgid],bexcl,shift,x,
944 box_size,nrnb,
945 lambda,dvdlambda,grps,FALSE,FALSE,
946 bHaveLJ);
947 nlr_ljc[jgid]=0;
948 }
949 nl_lr_ljc[jgid][nlr_ljc[jgid]++]=jjcg;
950 }
951 else {
952 if (nlr_coul[jgid] >= MAX_CG) {
953 do_longrange(log,top,fr,ngid,md,icg,jgid,
954 nlr_coul[jgid],
955 nl_lr_coul[jgid],bexcl,shift,x,
956 box_size,nrnb,
957 lambda,dvdlambda,grps,TRUE,FALSE,
958 bHaveLJ);
959 nlr_coul[jgid]=0;
960 }
961 nl_lr_coul[jgid][nlr_coul[jgid]++]=jjcg;
962 }
963 }
964 nns++;
965 }
966 }
967 }
968 }
969 }
970 /* CHECK whether there is anything left in the buffers */
971 for(nn=0; (nn<ngid); nn++) {
972 if (nsr[nn] > 0)
973 put_in_list(bHaveLJ,fr->nWater,
974 ngid,md,icg,nn,nsr[nn],nl_sr[nn],
975 cgsindex,cgsatoms,bexcl,shift,fr,FALSE,FALSE);
976
977 if (nlr_ljc[nn] > 0)
978 do_longrange(log,top,fr,ngid,md,icg,nn,nlr_ljc[nn],
979 nl_lr_ljc[nn],bexcl,shift,x,box_size,nrnb,
980 lambda,dvdlambda,grps,FALSE,FALSE,bHaveLJ);
981
982 if (nlr_coul[nn] > 0)
983 do_longrange(log,top,fr,ngid,md,icg,nn,nlr_coul[nn],
984 nl_lr_coul[nn],bexcl,shift,x,box_size,nrnb,
985 lambda,dvdlambda,grps,TRUE,FALSE,bHaveLJ);
986 }
987 }
988 }
989 }
990 setexcl(nri,i_atoms,&top->atoms.excl,FALSE,bexcl);
991 }
992 /* Perform any left over force calculations */
993 for (nn=0; (nn<ngid); nn++) {
994 if (rvdw2 > rs2)
995 do_longrange(log,top,fr,0,md,icg,nn,nlr_ljc[nn],
996 nl_lr_ljc[nn],bexcl,shift,x,box_size,nrnb,
997 lambda,dvdlambda,grps,FALSE,TRUE,bHaveLJ);
998 if (rcoul2 > rvdw2)
999 do_longrange(log,top,fr,0,md,icg,nn,nlr_coul[nn],
1000 nl_lr_coul[nn],bexcl,shift,x,box_size,nrnb,
1001 lambda,dvdlambda,grps,TRUE,TRUE,bHaveLJ);
1002 }
1003 /* Close off short range neighbourlists */
1004 close_neighbor_list(fr,FALSE,-1);
1005
1006 return nns;
1007 }
1008
1009 static rvec *sptr;
1010 static int sdim;
1011
1012 static int rv_comp(const void *a,const void *b)
1013 {
1014 int ia = *(int *)a;
1015 int ib = *(int *)b;
1016 real diff;
1017
1018 diff = sptr[ia][sdim] - sptr[ib][sdim];
1019 if (diff < 0)
1020 return -1;
1021 else if (diff == 0)
1022 return 0;
1023 else
1024 return 1;
1025 }
1026
1027 static void sort_charge_groups(t_commrec *cr,int cg_index[],int slab_index[],
1028 rvec cg_cm[],int Dimension)
1029 {
1030 int i,nrcg,cgind;
1031
1032 nrcg = slab_index[cr->nprocs];
1033 sptr = cg_cm;
1034 sdim = Dimension;
1035 qsort(cg_index,nrcg,sizeof(cg_index[0]),rv_comp);
1036
1037 if (debug) {
1038 fprintf(debug,"Just sorted the cg_cm array on dimension %d\n",Dimension);
1039 fprintf(debug,"Index: Coordinates of cg_cm\n");
1040 for(i=0; (i<nrcg); i++) {
1041 cgind = cg_index[i];
1042 fprintf(debug,"%8d%10.3f%10.3f%10.3f\n",cgind,
1043 cg_cm[cgind][XX],cg_cm[cgind][YY],cg_cm[cgind][ZZ]);
1044 }
1045 }
1046 sptr = NULL;
1047 sdim = -1;
1048 }
1049
1050 int search_neighbours(FILE *log,t_forcerec *fr,
1051 rvec x[],matrix box,
1052 t_topology *top,t_groups *grps,
1053 t_commrec *cr,t_nsborder *nsb,
1054 t_nrnb *nrnb,t_mdatoms *md,
1055 real lambda,real *dvdlambda)
1056 {
1057 static bool bFirst=TRUE;
1058 static t_grid *grid=NULL;
1059 static t_excl *bexcl;
1060 static bool *bHaveLJ;
1061 static t_ns_buf **ns_buf=NULL;
1062 static int *cg_index=NULL,*slab_index=NULL;
1063 static bool bSwitched=FALSE;
1064
1065 t_block *cgs=&(top->blocks[ebCGS]);
1066 rvec box_size;
1067 int i,j,m,ngid;
1068
1069 int nsearch;
1070 bool bGrid;
1071 char *ptr;
1072
1073 /* Set some local variables */
1074 bGrid=fr->bGrid;
1075 ngid=top->atoms.grps[egcENER].nr;
1076
1077 for(m=0; (m<DIM); m++)
1078 box_size[m]=box[m][m];
1079
1080 /* First time initiation of arrays etc. */
1081 if (bFirst) {
1082 int icg,nr_in_cg,maxcg;
1083
1084 /* Compute largest charge groups size (# atoms) */
1085 nr_in_cg=1;
1086 for (icg=0; (icg < cgs->nr); icg++)
1087 nr_in_cg=max(nr_in_cg,(int)(cgs->index[icg+1]-cgs->index[icg]));
1088
1089 /* Verify whether largest charge group is <= max cg.
1090 * This is determined by the type of the local exclusion type
1091 * Exclusions are stored in bits. (If the type is not large
1092 * enough, enlarge it, unsigned char -> unsigned short -> unsigned long)
1093 */
1094 maxcg=sizeof(t_excl)*8;
1095 if (nr_in_cg > maxcg)
1096 fatal_error(0,"Max #atoms in a charge group: %d > %d\n",
1097 nr_in_cg,maxcg);
1098
1099 snew(bexcl,cgs->nra);
1100 debug_gmx();
1101
1102 if ((ptr=getenv("NLIST")) != NULL) {
1103 sscanf(ptr,"%d",&NLJ_INC);
1104
1105 fprintf(log,"%s: I will increment J-lists by %d\n",
1106 __FILE__,NLJ_INC);
1107 }
1108
1109 /* Check whether we have to do domain decomposition,
1110 * if so set local variables for the charge group index and the
1111 * slab index.
1112 */
1113 if (fr->bDomDecomp) {
1114 snew(slab_index,cr->nprocs+1);
1115 for(i=0; (i<=cr->nprocs); i++)
1116 slab_index[i] = i*((real) cgs->nr/((real) cr->nprocs));
1117 fr->cg0 = slab_index[cr->pid];
1118 fr->hcg = slab_index[cr->pid+1] - fr->cg0;
1119 if (debug)
1120 fprintf(debug,"Will use DOMAIN DECOMPOSITION, "
1121 "from charge group index %d to %d on cpu %d\n",
1122 fr->cg0,fr->cg0+fr->hcg,cr->pid);
1123 }
1124 snew(cg_index,cgs->nr+1);
1125 for(i=0; (i<=cgs->nr); i++)
1126 cg_index[i] = i;
1127
1128 if (bGrid) {
1129 snew(grid,1);
1130 init_grid(log,grid,fr->ndelta,box,fr->rlistlong,cgs->nr);
1131 }
1132
1133 /* Create array that determines whether or not atoms have LJ */
1134 snew(bHaveLJ,fr->ntype);
1135 for(i=0; (i<fr->ntype); i++) {
1136 for(j=0; (j<fr->ntype); j++) {
1137 bHaveLJ[i] = (bHaveLJ[i] ||
1138 (fr->bBHAM ?
1139 ((BHAMA(fr->nbfp,fr->ntype,i,j) != 0) ||
1140 (BHAMB(fr->nbfp,fr->ntype,i,j) != 0) ||
1141 (BHAMC(fr->nbfp,fr->ntype,i,j) != 0)) :
1142 ((C6(fr->nbfp,fr->ntype,i,j) != 0) ||
1143 (C12(fr->nbfp,fr->ntype,i,j) != 0))));
1144 }
1145 }
1146 if (debug)
1147 pr_ivec(debug,0,"bHaveLJ",bHaveLJ,fr->ntype);
1148
1149 bFirst=FALSE;
1150 }
1151 debug_gmx();
1152
1153 /* Reset the neighbourlists */
1154 reset_neighbor_list(fr,FALSE,-1);
1155
1156 if (bGrid) {
1157 grid_first(log,grid,box,fr->rlistlong);
1158 /* Check if box is big enough to do grid searching... */
1159 if ( !( (grid->nrx >= 2*grid->delta+1) &&
1160 (grid->nry >= 2*grid->delta+1) &&
1161 (grid->nrz >= 2*grid->delta+1) ) ) {
1162 if (!bSwitched)
1163 fprintf(log,"WARNING: Box too small for grid-search, "
1164 "switching to simple neighboursearch.\n");
1165 if (fr->bTwinRange)
1166 fatal_error(0,"TWIN-RANGE cut-off with Simple "
1167 "neighboursearching not implemented.\n"
1168 "Use grid neighboursearching, and make (rlong < 0.4 box)");
1169 bGrid=FALSE;
1170 bSwitched=TRUE;
1171 } else {
1172 if (bSwitched)
1173 fprintf(log,"WARNING: Box large enough again for grid-search\n");
1174 bSwitched=FALSE;
1175 }
1176 }
1177 debug_gmx();
1178
1179 if (bGrid) {
1180 /* Don't know why this all is... (DvdS 3/99) */
1181 #ifndef SEGV
1182 int start = 0;
1183 int end = cgs->nr;
1184 #else
1185 int start = fr->cg0;
1186 int end = (cgs->nr+1)/2;
1187 #endif
1188
1189 if (fr->bDomDecomp)
1190 sort_charge_groups(cr,cg_index,slab_index,fr->cg_cm,fr->Dimension);
1191 debug_gmx();
1192
1193 fill_grid(log,fr->bDomDecomp,cg_index,
1194 grid,box,cgs->nr,fr->cg0,fr->hcg,fr->cg_cm);
1195 debug_gmx();
1196
1197 if (PAR(cr))
1198 mv_grid(cr,fr->bDomDecomp,cg_index,grid,nsb->workload);
1199 debug_gmx();
1200
1201 calc_elemnr(log,fr->bDomDecomp,cg_index,grid,start,end,cgs->nr);
1202 calc_ptrs(grid);
1203 grid_last(log,fr->bDomDecomp,cg_index,grid,start,end,cgs->nr);
1204
1205 if (debug) {
1206 check_grid(debug,grid);
1207 print_grid(debug,grid,fr->bDomDecomp,cg_index);
1208 }
1209 }
1210 debug_gmx();
1211
1212 /* Do the core! */
1213 if (bGrid)
1214 nsearch = ns5_core(log,fr,cg_index,box_size,ngid,top,grps,
1215 grid,x,bexcl,nrnb,md,lambda,dvdlambda,bHaveLJ);
1216 else {
1217 /* Only allocate this when necessary, saves 100 kb */
1218 if (ns_buf == NULL) {
1219 snew(ns_buf,ngid);
1220 for(i=0; (i<ngid); i++)
1221 snew(ns_buf[i],SHIFTS);
1222 }
1223 nsearch = ns_simple_core(fr,top,md,box_size,
1224 bexcl,ngid,ns_buf,bHaveLJ);
1225 }
1226 debug_gmx();
1227
1228 #ifdef DEBUG
1229 pr_nsblock(log);
1230 #endif
1231
1232 inc_nrnb(nrnb,eNR_NS,nsearch);
1233 /* inc_nrnb(nrnb,eNR_LR,fr->nlr); */
1234
1235 return nsearch;
1236 }
1237
1238
got merge a long time ago