?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
removed a wrong message about cg-cg kernels in ns.c
[gromacs/adressmacs.git] / src / mdlib / ns.c
blob61bdfb93e01cf905aa81f03b13d6f90ca1bfeaac
1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
2 *
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 3.2.0
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
15
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
20 *
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
27 *
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
30 *
31 * For more info, check our website at http://www.gromacs.org
32 *
33 * And Hey:
34 * GROwing Monsters And Cloning Shrimps
35 */
36 #ifdef HAVE_CONFIG_H
37 #include <config.h>
38 #endif
39
40 #ifdef GMX_THREAD_SHM_FDECOMP
41 #include <pthread.h>
42 #endif
43
44 #include <math.h>
45 #include <string.h>
46 #include "sysstuff.h"
47 #include "smalloc.h"
48 #include "macros.h"
49 #include "maths.h"
50 #include "vec.h"
51 #include "network.h"
52 #include "nsgrid.h"
53 #include "force.h"
54 #include "nonbonded.h"
55 #include "ns.h"
56 #include "pbc.h"
57 #include "names.h"
58 #include "gmx_fatal.h"
59 #include "nrnb.h"
60 #include "txtdump.h"
61 #include "mtop_util.h"
62
63 #include "domdec.h"
64 #include "adress.h"
65
66
67 /* Uncomment the define below to use the generic charge group - charge group
68 * inner loops (in src/gmxlib/nonbonded/nb_generic_cg.c).
69 * All intra charge-group interaction should be excluded.
70 * There should be no inter charge group exclusions.
71 * There can be no perturbed LJ types or charges.
72 * mdrun currently does NOT check for this.
73 */
74 /* #define GMX_CG_INNERLOOP */
75
76
77 /*
78 * E X C L U S I O N H A N D L I N G
79 */
80
81 #ifdef DEBUG
82 static void SETEXCL_(t_excl e[],atom_id i,atom_id j)
83 { e[j] = e[j] | (1<<i); }
84 static void RMEXCL_(t_excl e[],atom_id i,atom_id j)
85 { e[j]=e[j] & ~(1<<i); }
86 static bool ISEXCL_(t_excl e[],atom_id i,atom_id j)
87 { return (bool)(e[j] & (1<<i)); }
88 static bool NOTEXCL_(t_excl e[],atom_id i,atom_id j)
89 { return !(ISEXCL(e,i,j)); }
90 #else
91 #define SETEXCL(e,i,j) (e)[((atom_id) (j))] |= (1<<((atom_id) (i)))
92 #define RMEXCL(e,i,j) (e)[((atom_id) (j))] &= (~(1<<((atom_id) (i))))
93 #define ISEXCL(e,i,j) (bool) ((e)[((atom_id) (j))] & (1<<((atom_id) (i))))
94 #define NOTEXCL(e,i,j) !(ISEXCL(e,i,j))
95 #endif
96
97 /************************************************
98 *
99 * U T I L I T I E S F O R N S
100 *
101 ************************************************/
102
103 static void reallocate_nblist(t_nblist *nl)
104 {
105 if (gmx_debug_at)
106 {
107 fprintf(debug,"reallocating neigborlist il_code=%d, maxnri=%d\n",
108 nl->il_code,nl->maxnri);
109 }
110 srenew(nl->iinr, nl->maxnri);
111 if (nl->enlist == enlistCG_CG)
112 {
113 srenew(nl->iinr_end,nl->maxnri);
114 }
115 srenew(nl->gid, nl->maxnri);
116 srenew(nl->shift, nl->maxnri);
117 srenew(nl->jindex, nl->maxnri+1);
118 }
119
120 /* ivdw/icoul are used to determine the type of interaction, so we
121 * can set an innerloop index here. The obvious choice for this would have
122 * been the vdwtype/coultype values in the forcerecord, but unfortunately
123 * those types are braindead - for instance both Buckingham and normal
124 * Lennard-Jones use the same value (evdwCUT), and a separate boolean variable
125 * to determine which interaction is used. There is further no special value
126 * for 'no interaction'. For backward compatibility with old TPR files we won't
127 * change this in the 3.x series, so when calling this routine you should use:
128 *
129 * icoul=0 no coulomb interaction
130 * icoul=1 cutoff standard coulomb
131 * icoul=2 reaction-field coulomb
132 * icoul=3 tabulated coulomb
133 *
134 * ivdw=0 no vdw interaction
135 * ivdw=1 standard L-J interaction
136 * ivdw=2 Buckingham
137 * ivdw=3 tabulated vdw.
138 *
139 * Kind of ugly, but it works.
140 */
141 static void init_nblist(t_nblist *nl_sr,t_nblist *nl_lr,
142 int maxsr,int maxlr,
143 int ivdw, int icoul,
144 bool bfree, int enlist)
145 {
146 t_nblist *nl;
147 int homenr;
148 int i,nn;
149
150 int inloop[20] =
151 {
152 eNR_NBKERNEL_NONE,
153 eNR_NBKERNEL010,
154 eNR_NBKERNEL020,
155 eNR_NBKERNEL030,
156 eNR_NBKERNEL100,
157 eNR_NBKERNEL110,
158 eNR_NBKERNEL120,
159 eNR_NBKERNEL130,
160 eNR_NBKERNEL200,
161 eNR_NBKERNEL210,
162 eNR_NBKERNEL220,
163 eNR_NBKERNEL230,
164 eNR_NBKERNEL300,
165 eNR_NBKERNEL310,
166 eNR_NBKERNEL320,
167 eNR_NBKERNEL330,
168 eNR_NBKERNEL400,
169 eNR_NBKERNEL410,
170 eNR_NBKERNEL_NONE,
171 eNR_NBKERNEL430
172 };
173
174 for(i=0; (i<2); i++)
175 {
176 nl = (i == 0) ? nl_sr : nl_lr;
177 homenr = (i == 0) ? maxsr : maxlr;
178
179 if (nl == NULL)
180 {
181 continue;
182 }
183
184 /* Set coul/vdw in neighborlist, and for the normal loops we determine
185 * an index of which one to call.
186 */
187 nl->ivdw = ivdw;
188 nl->icoul = icoul;
189 nl->free_energy = bfree;
190
191 if (bfree)
192 {
193 nl->enlist = enlistATOM_ATOM;
194 nl->il_code = eNR_NBKERNEL_FREE_ENERGY;
195 }
196 else
197 {
198 nl->enlist = enlist;
199
200 nn = inloop[4*icoul + ivdw];
201
202 /* solvent loops follow directly after the corresponding
203 * ordinary loops, in the order:
204 *
205 * SPC, SPC-SPC, TIP4p, TIP4p-TIP4p
206 *
207 */
208 switch (enlist) {
209 case enlistATOM_ATOM:
210 case enlistCG_CG:
211 break;
212 case enlistSPC_ATOM: nn += 1; break;
213 case enlistSPC_SPC: nn += 2; break;
214 case enlistTIP4P_ATOM: nn += 3; break;
215 case enlistTIP4P_TIP4P: nn += 4; break;
216 }
217
218 nl->il_code = nn;
219 }
220
221 if (debug)
222 fprintf(debug,"Initiating neighbourlist type %d for %s interactions,\nwith %d SR, %d LR atoms.\n",
223 nl->il_code,ENLISTTYPE(enlist),maxsr,maxlr);
224
225 /* maxnri is influenced by the number of shifts (maximum is 8)
226 * and the number of energy groups.
227 * If it is not enough, nl memory will be reallocated during the run.
228 * 4 seems to be a reasonable factor, which only causes reallocation
229 * during runs with tiny and many energygroups.
230 */
231 nl->maxnri = homenr*4;
232 nl->maxnrj = 0;
233 nl->maxlen = 0;
234 nl->nri = -1;
235 nl->nrj = 0;
236 nl->iinr = NULL;
237 nl->gid = NULL;
238 nl->shift = NULL;
239 nl->jindex = NULL;
240 reallocate_nblist(nl);
241 nl->jindex[0] = 0;
242 #ifdef GMX_THREAD_SHM_FDECOMP
243 nl->counter = 0;
244 snew(nl->mtx,1);
245 pthread_mutex_init(nl->mtx,NULL);
246 #endif
247 }
248 }
249
250 void init_neighbor_list(FILE *log,t_forcerec *fr,int homenr)
251 {
252 /* Make maxlr tunable! (does not seem to be a big difference though)
253 * This parameter determines the number of i particles in a long range
254 * neighbourlist. Too few means many function calls, too many means
255 * cache trashing.
256 */
257 int maxsr,maxsr_wat,maxlr,maxlr_wat;
258 int icoul,icoulf,ivdw;
259 int solvent;
260 int enlist_def,enlist_w,enlist_ww;
261 int i;
262 t_nblists *nbl;
263
264 /* maxsr = homenr-fr->nWatMol*3; */
265 maxsr = homenr;
266
267 if (maxsr < 0)
268 {
269 gmx_fatal(FARGS,"%s, %d: Negative number of short range atoms.\n"
270 "Call your Gromacs dealer for assistance.",__FILE__,__LINE__);
271 }
272 /* This is just for initial allocation, so we do not reallocate
273 * all the nlist arrays many times in a row.
274 * The numbers seem very accurate, but they are uncritical.
275 */
276 maxsr_wat = min(fr->nWatMol,(homenr+2)/3);
277 if (fr->bTwinRange)
278 {
279 maxlr = 50;
280 maxlr_wat = min(maxsr_wat,maxlr);
281 }
282 else
283 {
284 maxlr = maxlr_wat = 0;
285 }
286
287 /* Determine the values for icoul/ivdw. */
288 /* Start with GB */
289 if(fr->bGB)
290 {
291 icoul=4;
292 }
293 else if (fr->bcoultab)
294 {
295 icoul = 3;
296 }
297 else if (EEL_RF(fr->eeltype))
298 {
299 icoul = 2;
300 }
301 else
302 {
303 icoul = 1;
304 }
305
306 if (fr->bvdwtab)
307 {
308 ivdw = 3;
309 }
310 else if (fr->bBHAM)
311 {
312 ivdw = 2;
313 }
314 else
315 {
316 ivdw = 1;
317 }
318
319 fr->ns.bCGlist = (getenv("GMX_NBLISTCG") != 0);
320 if (!fr->ns.bCGlist)
321 {
322 enlist_def = enlistATOM_ATOM;
323 }
324 else
325 {
326 enlist_def = enlistCG_CG;
327 if (log != NULL)
328 {
329 fprintf(log,"\nUsing charge-group - charge-group neighbor lists and kernels\n\n");
330 }
331 }
332
333 if (fr->solvent_opt == esolTIP4P) {
334 enlist_w = enlistTIP4P_ATOM;
335 enlist_ww = enlistTIP4P_TIP4P;
336 } else {
337 enlist_w = enlistSPC_ATOM;
338 enlist_ww = enlistSPC_SPC;
339 }
340
341 for(i=0; i<fr->nnblists; i++)
342 {
343 nbl = &(fr->nblists[i]);
344 init_nblist(&nbl->nlist_sr[eNL_VDWQQ],&nbl->nlist_lr[eNL_VDWQQ],
345 maxsr,maxlr,ivdw,icoul,FALSE,enlist_def);
346 init_nblist(&nbl->nlist_sr[eNL_VDW],&nbl->nlist_lr[eNL_VDW],
347 maxsr,maxlr,ivdw,0,FALSE,enlist_def);
348 init_nblist(&nbl->nlist_sr[eNL_QQ],&nbl->nlist_lr[eNL_QQ],
349 maxsr,maxlr,0,icoul,FALSE,enlist_def);
350 init_nblist(&nbl->nlist_sr[eNL_VDWQQ_WATER],&nbl->nlist_lr[eNL_VDWQQ_WATER],
351 maxsr_wat,maxlr_wat,ivdw,icoul, FALSE,enlist_w);
352 init_nblist(&nbl->nlist_sr[eNL_QQ_WATER],&nbl->nlist_lr[eNL_QQ_WATER],
353 maxsr_wat,maxlr_wat,0,icoul, FALSE,enlist_w);
354 init_nblist(&nbl->nlist_sr[eNL_VDWQQ_WATERWATER],&nbl->nlist_lr[eNL_VDWQQ_WATERWATER],
355 maxsr_wat,maxlr_wat,ivdw,icoul, FALSE,enlist_ww);
356 init_nblist(&nbl->nlist_sr[eNL_QQ_WATERWATER],&nbl->nlist_lr[eNL_QQ_WATERWATER],
357 maxsr_wat,maxlr_wat,0,icoul, FALSE,enlist_ww);
358
359 if (fr->efep != efepNO)
360 {
361 if (fr->bEwald)
362 {
363 icoulf = 5;
364 }
365 else
366 {
367 icoulf = icoul;
368 }
369
370 init_nblist(&nbl->nlist_sr[eNL_VDWQQ_FREE],&nbl->nlist_lr[eNL_VDWQQ_FREE],
371 maxsr,maxlr,ivdw,icoulf,TRUE,enlistATOM_ATOM);
372 init_nblist(&nbl->nlist_sr[eNL_VDW_FREE],&nbl->nlist_lr[eNL_VDW_FREE],
373 maxsr,maxlr,ivdw,0,TRUE,enlistATOM_ATOM);
374 init_nblist(&nbl->nlist_sr[eNL_QQ_FREE],&nbl->nlist_lr[eNL_QQ_FREE],
375 maxsr,maxlr,0,icoulf,TRUE,enlistATOM_ATOM);
376 }
377 }
378 /* QMMM MM list */
379 if (fr->bQMMM && fr->qr->QMMMscheme != eQMMMschemeoniom)
380 {
381 init_nblist(&fr->QMMMlist,NULL,
382 maxsr,maxlr,0,icoul,FALSE,enlistATOM_ATOM);
383 }
384
385 fr->ns.nblist_initialized=TRUE;
386 }
387
388 static void reset_nblist(t_nblist *nl)
389 {
390 nl->nri = -1;
391 nl->nrj = 0;
392 nl->maxlen = 0;
393 if (nl->jindex)
394 {
395 nl->jindex[0] = 0;
396 }
397 }
398
399 static void reset_neighbor_list(t_forcerec *fr,bool bLR,int nls,int eNL)
400 {
401 int n,i;
402
403 if (bLR)
404 {
405 reset_nblist(&(fr->nblists[nls].nlist_lr[eNL]));
406 }
407 else
408 {
409 for(n=0; n<fr->nnblists; n++)
410 {
411 for(i=0; i<eNL_NR; i++)
412 {
413 reset_nblist(&(fr->nblists[n].nlist_sr[i]));
414 }
415 }
416 if (fr->bQMMM)
417 {
418 /* only reset the short-range nblist */
419 reset_nblist(&(fr->QMMMlist));
420 }
421 }
422 }
423
424
425
426
427 static inline void new_i_nblist(t_nblist *nlist,
428 bool bLR,atom_id i_atom,int shift,int gid)
429 {
430 int i,k,nri,nshift;
431
432 nri = nlist->nri;
433
434 /* Check whether we have to increase the i counter */
435 if ((nri == -1) ||
436 (nlist->iinr[nri] != i_atom) ||
437 (nlist->shift[nri] != shift) ||
438 (nlist->gid[nri] != gid))
439 {
440 /* This is something else. Now see if any entries have
441 * been added in the list of the previous atom.
442 */
443 if ((nri == -1) ||
444 ((nlist->jindex[nri+1] > nlist->jindex[nri]) &&
445 (nlist->gid[nri] != -1)))
446 {
447 /* If so increase the counter */
448 nlist->nri++;
449 nri++;
450 if (nlist->nri >= nlist->maxnri)
451 {
452 nlist->maxnri += over_alloc_large(nlist->nri);
453 reallocate_nblist(nlist);
454 }
455 }
456 /* Set the number of neighbours and the atom number */
457 nlist->jindex[nri+1] = nlist->jindex[nri];
458 nlist->iinr[nri] = i_atom;
459 nlist->gid[nri] = gid;
460 nlist->shift[nri] = shift;
461 }
462 }
463
464 static inline void close_i_nblist(t_nblist *nlist)
465 {
466 int nri = nlist->nri;
467 int len;
468
469 if (nri >= 0)
470 {
471 nlist->jindex[nri+1] = nlist->nrj;
472
473 len=nlist->nrj - nlist->jindex[nri];
474
475 /* nlist length for water i molecules is treated statically
476 * in the innerloops
477 */
478 if (len > nlist->maxlen)
479 {
480 nlist->maxlen = len;
481 }
482 }
483 }
484
485 static inline void close_nblist(t_nblist *nlist)
486 {
487 /* Only close this nblist when it has been initialized */
488 if (nlist->jindex)
489 {
490 nlist->nri++;
491 }
492 }
493
494 static inline void close_neighbor_list(t_forcerec *fr,bool bLR,int nls,int eNL,
495 bool bMakeQMMMnblist)
496 {
497 int n,i;
498
499 if (bMakeQMMMnblist) {
500 if (!bLR)
501 {
502 close_nblist(&(fr->QMMMlist));
503 }
504 }
505 else
506 {
507 if (bLR)
508 {
509 close_nblist(&(fr->nblists[nls].nlist_lr[eNL]));
510 }
511 else
512 {
513 for(n=0; n<fr->nnblists; n++)
514 {
515 for(i=0; (i<eNL_NR); i++)
516 {
517 close_nblist(&(fr->nblists[n].nlist_sr[i]));
518 }
519 }
520 }
521 }
522 }
523
524 static inline void add_j_to_nblist(t_nblist *nlist,atom_id j_atom,bool bLR)
525 {
526 int nrj=nlist->nrj;
527
528 if (nlist->nrj >= nlist->maxnrj)
529 {
530 nlist->maxnrj = over_alloc_small(nlist->nrj + 1);
531 if (gmx_debug_at)
532 fprintf(debug,"Increasing %s nblist %s j size to %d\n",
533 bLR ? "LR" : "SR",nrnb_str(nlist->il_code),nlist->maxnrj);
534
535 srenew(nlist->jjnr,nlist->maxnrj);
536 }
537
538 nlist->jjnr[nrj] = j_atom;
539 nlist->nrj ++;
540 }
541
542 static inline void add_j_to_nblist_cg(t_nblist *nlist,
543 atom_id j_start,int j_end,
544 t_excl *bexcl,bool bLR)
545 {
546 int nrj=nlist->nrj;
547 int j;
548
549 if (nlist->nrj >= nlist->maxnrj)
550 {
551 nlist->maxnrj = over_alloc_small(nlist->nrj + 1);
552 if (gmx_debug_at)
553 fprintf(debug,"Increasing %s nblist %s j size to %d\n",
554 bLR ? "LR" : "SR",nrnb_str(nlist->il_code),nlist->maxnrj);
555
556 srenew(nlist->jjnr ,nlist->maxnrj);
557 srenew(nlist->jjnr_end,nlist->maxnrj);
558 srenew(nlist->excl ,nlist->maxnrj*MAX_CGCGSIZE);
559 }
560
561 nlist->jjnr[nrj] = j_start;
562 nlist->jjnr_end[nrj] = j_end;
563
564 if (j_end - j_start > MAX_CGCGSIZE)
565 {
566 gmx_fatal(FARGS,"The charge-group - charge-group neighborlist do not support charge groups larger than %d, found a charge group of size %d",MAX_CGCGSIZE,j_end-j_start);
567 }
568
569 /* Set the exclusions */
570 for(j=j_start; j<j_end; j++)
571 {
572 nlist->excl[nrj*MAX_CGCGSIZE + j - j_start] = bexcl[j];
573 }
574
575 nlist->nrj ++;
576 }
577
578 typedef void
579 put_in_list_t(bool bHaveVdW[],
580 int ngid,
581 t_mdatoms * md,
582 int icg,
583 int jgid,
584 int nj,
585 atom_id jjcg[],
586 atom_id index[],
587 t_excl bExcl[],
588 int shift,
589 t_forcerec * fr,
590 bool bLR,
591 bool bDoVdW,
592 bool bDoCoul);
593
594 static void
595 put_in_list_at(bool bHaveVdW[],
596 int ngid,
597 t_mdatoms * md,
598 int icg,
599 int jgid,
600 int nj,
601 atom_id jjcg[],
602 atom_id index[],
603 t_excl bExcl[],
604 int shift,
605 t_forcerec * fr,
606 bool bLR,
607 bool bDoVdW,
608 bool bDoCoul)
609 {
610 /* The a[] index has been removed,
611 * to put it back in i_atom should be a[i0] and jj should be a[jj].
612 */
613 t_nblist * vdwc;
614 t_nblist * vdw;
615 t_nblist * coul;
616 t_nblist * vdwc_free = NULL;
617 t_nblist * vdw_free = NULL;
618 t_nblist * coul_free = NULL;
619 t_nblist * vdwc_ww = NULL;
620 t_nblist * coul_ww = NULL;
621
622 int i,j,jcg,igid,gid,nbl_ind,ind_ij;
623 atom_id jj,jj0,jj1,i_atom;
624 int i0,nicg,len;
625
626 int *cginfo;
627 int *type,*typeB;
628 real *charge,*chargeB;
629 real qi,qiB,qq,rlj;
630 bool bFreeEnergy,bFree,bFreeJ,bNotEx,*bPert;
631 bool bDoVdW_i,bDoCoul_i,bDoCoul_i_sol;
632 int iwater,jwater;
633 t_nblist *nlist;
634
635 /* Copy some pointers */
636 cginfo = fr->cginfo;
637 charge = md->chargeA;
638 chargeB = md->chargeB;
639 type = md->typeA;
640 typeB = md->typeB;
641 bPert = md->bPerturbed;
642
643 /* Get atom range */
644 i0 = index[icg];
645 nicg = index[icg+1]-i0;
646
647 /* Get the i charge group info */
648 igid = GET_CGINFO_GID(cginfo[icg]);
649 iwater = GET_CGINFO_SOLOPT(cginfo[icg]);
650
651 bFreeEnergy = FALSE;
652 if (md->nPerturbed)
653 {
654 /* Check if any of the particles involved are perturbed.
655 * If not we can do the cheaper normal put_in_list
656 * and use more solvent optimization.
657 */
658 for(i=0; i<nicg; i++)
659 {
660 bFreeEnergy |= bPert[i0+i];
661 }
662 /* Loop over the j charge groups */
663 for(j=0; (j<nj && !bFreeEnergy); j++)
664 {
665 jcg = jjcg[j];
666 jj0 = index[jcg];
667 jj1 = index[jcg+1];
668 /* Finally loop over the atoms in the j-charge group */
669 for(jj=jj0; jj<jj1; jj++)
670 {
671 bFreeEnergy |= bPert[jj];
672 }
673 }
674 }
675
676 /* Unpack pointers to neighbourlist structs */
677 if (fr->nnblists == 1)
678 {
679 nbl_ind = 0;
680 }
681 else
682 {
683 nbl_ind = fr->gid2nblists[GID(igid,jgid,ngid)];
684 }
685 if (bLR)
686 {
687 nlist = fr->nblists[nbl_ind].nlist_lr;
688 }
689 else
690 {
691 nlist = fr->nblists[nbl_ind].nlist_sr;
692 }
693
694 if (iwater != esolNO)
695 {
696 vdwc = &nlist[eNL_VDWQQ_WATER];
697 vdw = &nlist[eNL_VDW];
698 coul = &nlist[eNL_QQ_WATER];
699 #ifndef DISABLE_WATERWATER_NLIST
700 vdwc_ww = &nlist[eNL_VDWQQ_WATERWATER];
701 coul_ww = &nlist[eNL_QQ_WATERWATER];
702 #endif
703 }
704 else
705 {
706 vdwc = &nlist[eNL_VDWQQ];
707 vdw = &nlist[eNL_VDW];
708 coul = &nlist[eNL_QQ];
709 }
710
711 if (!bFreeEnergy)
712 {
713 if (iwater != esolNO)
714 {
715 /* Loop over the atoms in the i charge group */
716 i_atom = i0;
717 gid = GID(igid,jgid,ngid);
718 /* Create new i_atom for each energy group */
719 if (bDoCoul && bDoVdW)
720 {
721 new_i_nblist(vdwc,bLR,i_atom,shift,gid);
722 #ifndef DISABLE_WATERWATER_NLIST
723 new_i_nblist(vdwc_ww,bLR,i_atom,shift,gid);
724 #endif
725 }
726 if (bDoVdW)
727 {
728 new_i_nblist(vdw,bLR,i_atom,shift,gid);
729 }
730 if (bDoCoul)
731 {
732 new_i_nblist(coul,bLR,i_atom,shift,gid);
733 #ifndef DISABLE_WATERWATER_NLIST
734 new_i_nblist(coul_ww,bLR,i_atom,shift,gid);
735 #endif
736 }
737 /* Loop over the j charge groups */
738 for(j=0; (j<nj); j++)
739 {
740 jcg=jjcg[j];
741
742 if (jcg == icg)
743 {
744 continue;
745 }
746
747 jj0 = index[jcg];
748 jwater = GET_CGINFO_SOLOPT(cginfo[jcg]);
749
750 if (iwater == esolSPC && jwater == esolSPC)
751 {
752 /* Interaction between two SPC molecules */
753 if (!bDoCoul)
754 {
755 /* VdW only - only first atoms in each water interact */
756 add_j_to_nblist(vdw,jj0,bLR);
757 }
758 else
759 {
760 #ifdef DISABLE_WATERWATER_NLIST
761 /* Add entries for the three atoms - only do VdW if we need to */
762 if (!bDoVdW)
763 {
764 add_j_to_nblist(coul,jj0,bLR);
765 }
766 else
767 {
768 add_j_to_nblist(vdwc,jj0,bLR);
769 }
770 add_j_to_nblist(coul,jj0+1,bLR);
771 add_j_to_nblist(coul,jj0+2,bLR);
772 #else
773 /* One entry for the entire water-water interaction */
774 if (!bDoVdW)
775 {
776 add_j_to_nblist(coul_ww,jj0,bLR);
777 }
778 else
779 {
780 add_j_to_nblist(vdwc_ww,jj0,bLR);
781 }
782 #endif
783 }
784 }
785 else if (iwater == esolTIP4P && jwater == esolTIP4P)
786 {
787 /* Interaction between two TIP4p molecules */
788 if (!bDoCoul)
789 {
790 /* VdW only - only first atoms in each water interact */
791 add_j_to_nblist(vdw,jj0,bLR);
792 }
793 else
794 {
795 #ifdef DISABLE_WATERWATER_NLIST
796 /* Add entries for the four atoms - only do VdW if we need to */
797 if (bDoVdW)
798 {
799 add_j_to_nblist(vdw,jj0,bLR);
800 }
801 add_j_to_nblist(coul,jj0+1,bLR);
802 add_j_to_nblist(coul,jj0+2,bLR);
803 add_j_to_nblist(coul,jj0+3,bLR);
804 #else
805 /* One entry for the entire water-water interaction */
806 if (!bDoVdW)
807 {
808 add_j_to_nblist(coul_ww,jj0,bLR);
809 }
810 else
811 {
812 add_j_to_nblist(vdwc_ww,jj0,bLR);
813 }
814 #endif
815 }
816 }
817 else
818 {
819 /* j charge group is not water, but i is.
820 * Add entries to the water-other_atom lists; the geometry of the water
821 * molecule doesn't matter - that is taken care of in the nonbonded kernel,
822 * so we don't care if it is SPC or TIP4P...
823 */
824
825 jj1 = index[jcg+1];
826
827 if (!bDoVdW)
828 {
829 for(jj=jj0; (jj<jj1); jj++)
830 {
831 if (charge[jj] != 0)
832 {
833 add_j_to_nblist(coul,jj,bLR);
834 }
835 }
836 }
837 else if (!bDoCoul)
838 {
839 for(jj=jj0; (jj<jj1); jj++)
840 {
841 if (bHaveVdW[type[jj]])
842 {
843 add_j_to_nblist(vdw,jj,bLR);
844 }
845 }
846 }
847 else
848 {
849 /* _charge_ _groups_ interact with both coulomb and LJ */
850 /* Check which atoms we should add to the lists! */
851 for(jj=jj0; (jj<jj1); jj++)
852 {
853 if (bHaveVdW[type[jj]])
854 {
855 if (charge[jj] != 0)
856 {
857 add_j_to_nblist(vdwc,jj,bLR);
858 }
859 else
860 {
861 add_j_to_nblist(vdw,jj,bLR);
862 }
863 }
864 else if (charge[jj] != 0)
865 {
866 add_j_to_nblist(coul,jj,bLR);
867 }
868 }
869 }
870 }
871 }
872 close_i_nblist(vdw);
873 close_i_nblist(coul);
874 close_i_nblist(vdwc);
875 #ifndef DISABLE_WATERWATER_NLIST
876 close_i_nblist(coul_ww);
877 close_i_nblist(vdwc_ww);
878 #endif
879 }
880 else
881 {
882 /* no solvent as i charge group */
883 /* Loop over the atoms in the i charge group */
884 for(i=0; i<nicg; i++)
885 {
886 i_atom = i0+i;
887 gid = GID(igid,jgid,ngid);
888 qi = charge[i_atom];
889
890 /* Create new i_atom for each energy group */
891 if (bDoVdW && bDoCoul)
892 {
893 new_i_nblist(vdwc,bLR,i_atom,shift,gid);
894 }
895 if (bDoVdW)
896 {
897 new_i_nblist(vdw,bLR,i_atom,shift,gid);
898 }
899 if (bDoCoul)
900 {
901 new_i_nblist(coul,bLR,i_atom,shift,gid);
902 }
903 bDoVdW_i = (bDoVdW && bHaveVdW[type[i_atom]]);
904 bDoCoul_i = (bDoCoul && qi!=0);
905
906 if (bDoVdW_i || bDoCoul_i)
907 {
908 /* Loop over the j charge groups */
909 for(j=0; (j<nj); j++)
910 {
911 jcg=jjcg[j];
912
913 /* Check for large charge groups */
914 if (jcg == icg)
915 {
916 jj0 = i0 + i + 1;
917 }
918 else
919 {
920 jj0 = index[jcg];
921 }
922
923 jj1=index[jcg+1];
924 /* Finally loop over the atoms in the j-charge group */
925 for(jj=jj0; jj<jj1; jj++)
926 {
927 bNotEx = NOTEXCL(bExcl,i,jj);
928
929 if (bNotEx)
930 {
931 if (!bDoVdW_i)
932 {
933 if (charge[jj] != 0)
934 {
935 add_j_to_nblist(coul,jj,bLR);
936 }
937 }
938 else if (!bDoCoul_i)
939 {
940 if (bHaveVdW[type[jj]])
941 {
942 add_j_to_nblist(vdw,jj,bLR);
943 }
944 }
945 else
946 {
947 if (bHaveVdW[type[jj]])
948 {
949 if (charge[jj] != 0)
950 {
951 add_j_to_nblist(vdwc,jj,bLR);
952 }
953 else
954 {
955 add_j_to_nblist(vdw,jj,bLR);
956 }
957 }
958 else if (charge[jj] != 0)
959 {
960 add_j_to_nblist(coul,jj,bLR);
961 }
962 }
963 }
964 }
965 }
966 }
967 close_i_nblist(vdw);
968 close_i_nblist(coul);
969 close_i_nblist(vdwc);
970 }
971 }
972 }
973 else
974 {
975 /* we are doing free energy */
976 vdwc_free = &nlist[eNL_VDWQQ_FREE];
977 vdw_free = &nlist[eNL_VDW_FREE];
978 coul_free = &nlist[eNL_QQ_FREE];
979 /* Loop over the atoms in the i charge group */
980 for(i=0; i<nicg; i++)
981 {
982 i_atom = i0+i;
983 gid = GID(igid,jgid,ngid);
984 qi = charge[i_atom];
985 qiB = chargeB[i_atom];
986
987 /* Create new i_atom for each energy group */
988 if (bDoVdW && bDoCoul)
989 new_i_nblist(vdwc,bLR,i_atom,shift,gid);
990 if (bDoVdW)
991 new_i_nblist(vdw,bLR,i_atom,shift,gid);
992 if (bDoCoul)
993 new_i_nblist(coul,bLR,i_atom,shift,gid);
994
995 new_i_nblist(vdw_free,bLR,i_atom,shift,gid);
996 new_i_nblist(coul_free,bLR,i_atom,shift,gid);
997 new_i_nblist(vdwc_free,bLR,i_atom,shift,gid);
998
999 bDoVdW_i = (bDoVdW &&
1000 (bHaveVdW[type[i_atom]] || bHaveVdW[typeB[i_atom]]));
1001 bDoCoul_i = (bDoCoul && (qi!=0 || qiB!=0));
1002 /* For TIP4P the first atom does not have a charge,
1003 * but the last three do. So we should still put an atom
1004 * without LJ but with charge in the water-atom neighborlist
1005 * for a TIP4p i charge group.
1006 * For SPC type water the first atom has LJ and charge,
1007 * so there is no such problem.
1008 */
1009 if (iwater == esolNO)
1010 {
1011 bDoCoul_i_sol = bDoCoul_i;
1012 }
1013 else
1014 {
1015 bDoCoul_i_sol = bDoCoul;
1016 }
1017
1018 if (bDoVdW_i || bDoCoul_i_sol)
1019 {
1020 /* Loop over the j charge groups */
1021 for(j=0; (j<nj); j++)
1022 {
1023 jcg=jjcg[j];
1024
1025 /* Check for large charge groups */
1026 if (jcg == icg)
1027 {
1028 jj0 = i0 + i + 1;
1029 }
1030 else
1031 {
1032 jj0 = index[jcg];
1033 }
1034
1035 jj1=index[jcg+1];
1036 /* Finally loop over the atoms in the j-charge group */
1037 bFree = bPert[i_atom];
1038 for(jj=jj0; (jj<jj1); jj++)
1039 {
1040 bFreeJ = bFree || bPert[jj];
1041 /* Complicated if, because the water H's should also
1042 * see perturbed j-particles
1043 */
1044 if (iwater==esolNO || i==0 || bFreeJ)
1045 {
1046 bNotEx = NOTEXCL(bExcl,i,jj);
1047
1048 if (bNotEx)
1049 {
1050 if (bFreeJ)
1051 {
1052 if (!bDoVdW_i)
1053 {
1054 if (charge[jj]!=0 || chargeB[jj]!=0)
1055 {
1056 add_j_to_nblist(coul_free,jj,bLR);
1057 }
1058 }
1059 else if (!bDoCoul_i)
1060 {
1061 if (bHaveVdW[type[jj]] || bHaveVdW[typeB[jj]])
1062 {
1063 add_j_to_nblist(vdw_free,jj,bLR);
1064 }
1065 }
1066 else
1067 {
1068 if (bHaveVdW[type[jj]] || bHaveVdW[typeB[jj]])
1069 {
1070 if (charge[jj]!=0 || chargeB[jj]!=0)
1071 {
1072 add_j_to_nblist(vdwc_free,jj,bLR);
1073 }
1074 else
1075 {
1076 add_j_to_nblist(vdw_free,jj,bLR);
1077 }
1078 }
1079 else if (charge[jj]!=0 || chargeB[jj]!=0)
1080 add_j_to_nblist(coul_free,jj,bLR);
1081 }
1082 }
1083 else if (!bDoVdW_i)
1084 {
1085 /* This is done whether or not bWater is set */
1086 if (charge[jj] != 0)
1087 {
1088 add_j_to_nblist(coul,jj,bLR);
1089 }
1090 }
1091 else if (!bDoCoul_i_sol)
1092 {
1093 if (bHaveVdW[type[jj]])
1094 {
1095 add_j_to_nblist(vdw,jj,bLR);
1096 }
1097 }
1098 else
1099 {
1100 if (bHaveVdW[type[jj]])
1101 {
1102 if (charge[jj] != 0)
1103 {
1104 add_j_to_nblist(vdwc,jj,bLR);
1105 }
1106 else
1107 {
1108 add_j_to_nblist(vdw,jj,bLR);
1109 }
1110 }
1111 else if (charge[jj] != 0)
1112 {
1113 add_j_to_nblist(coul,jj,bLR);
1114 }
1115 }
1116 }
1117 }
1118 }
1119 }
1120 }
1121 close_i_nblist(vdw);
1122 close_i_nblist(coul);
1123 close_i_nblist(vdwc);
1124 close_i_nblist(vdw_free);
1125 close_i_nblist(coul_free);
1126 close_i_nblist(vdwc_free);
1127 }
1128 }
1129 }
1130
1131 static void
1132 put_in_list_qmmm(bool bHaveVdW[],
1133 int ngid,
1134 t_mdatoms * md,
1135 int icg,
1136 int jgid,
1137 int nj,
1138 atom_id jjcg[],
1139 atom_id index[],
1140 t_excl bExcl[],
1141 int shift,
1142 t_forcerec * fr,
1143 bool bLR,
1144 bool bDoVdW,
1145 bool bDoCoul)
1146 {
1147 t_nblist * coul;
1148 int i,j,jcg,igid,gid;
1149 atom_id jj,jj0,jj1,i_atom;
1150 int i0,nicg;
1151 bool bNotEx;
1152
1153 /* Get atom range */
1154 i0 = index[icg];
1155 nicg = index[icg+1]-i0;
1156
1157 /* Get the i charge group info */
1158 igid = GET_CGINFO_GID(fr->cginfo[icg]);
1159
1160 coul = &fr->QMMMlist;
1161
1162 /* Loop over atoms in the ith charge group */
1163 for (i=0;i<nicg;i++)
1164 {
1165 i_atom = i0+i;
1166 gid = GID(igid,jgid,ngid);
1167 /* Create new i_atom for each energy group */
1168 new_i_nblist(coul,bLR,i_atom,shift,gid);
1169
1170 /* Loop over the j charge groups */
1171 for (j=0;j<nj;j++)
1172 {
1173 jcg=jjcg[j];
1174
1175 /* Charge groups cannot have QM and MM atoms simultaneously */
1176 if (jcg!=icg)
1177 {
1178 jj0 = index[jcg];
1179 jj1 = index[jcg+1];
1180 /* Finally loop over the atoms in the j-charge group */
1181 for(jj=jj0; jj<jj1; jj++)
1182 {
1183 bNotEx = NOTEXCL(bExcl,i,jj);
1184 if(bNotEx)
1185 add_j_to_nblist(coul,jj,bLR);
1186 }
1187 }
1188 }
1189 close_i_nblist(coul);
1190 }
1191 }
1192
1193 static void
1194 put_in_list_cg(bool bHaveVdW[],
1195 int ngid,
1196 t_mdatoms * md,
1197 int icg,
1198 int jgid,
1199 int nj,
1200 atom_id jjcg[],
1201 atom_id index[],
1202 t_excl bExcl[],
1203 int shift,
1204 t_forcerec * fr,
1205 bool bLR,
1206 bool bDoVdW,
1207 bool bDoCoul)
1208 {
1209 int cginfo;
1210 int igid,gid,nbl_ind;
1211 t_nblist * vdwc;
1212 int j,jcg;
1213
1214 cginfo = fr->cginfo[icg];
1215
1216 igid = GET_CGINFO_GID(cginfo);
1217 gid = GID(igid,jgid,ngid);
1218
1219 /* Unpack pointers to neighbourlist structs */
1220 if (fr->nnblists == 1)
1221 {
1222 nbl_ind = 0;
1223 }
1224 else
1225 {
1226 nbl_ind = fr->gid2nblists[gid];
1227 }
1228 if (bLR)
1229 {
1230 vdwc = &fr->nblists[nbl_ind].nlist_lr[eNL_VDWQQ];
1231 }
1232 else
1233 {
1234 vdwc = &fr->nblists[nbl_ind].nlist_sr[eNL_VDWQQ];
1235 }
1236
1237 /* Make a new neighbor list for charge group icg.
1238 * Currently simply one neighbor list is made with LJ and Coulomb.
1239 * If required, zero interactions could be removed here
1240 * or in the force loop.
1241 */
1242 new_i_nblist(vdwc,bLR,index[icg],shift,gid);
1243 vdwc->iinr_end[vdwc->nri] = index[icg+1];
1244
1245 for(j=0; (j<nj); j++)
1246 {
1247 jcg = jjcg[j];
1248 /* Skip the icg-icg pairs if all self interactions are excluded */
1249 if (!(jcg == icg && GET_CGINFO_EXCL_INTRA(cginfo)))
1250 {
1251 /* Here we add the j charge group jcg to the list,
1252 * exclusions are also added to the list.
1253 */
1254 add_j_to_nblist_cg(vdwc,index[jcg],index[jcg+1],bExcl,bLR);
1255 }
1256 }
1257
1258 close_i_nblist(vdwc);
1259 }
1260
1261 static void setexcl(atom_id start,atom_id end,t_blocka *excl,bool b,
1262 t_excl bexcl[])
1263 {
1264 atom_id i,k;
1265
1266 if (b)
1267 {
1268 for(i=start; i<end; i++)
1269 {
1270 for(k=excl->index[i]; k<excl->index[i+1]; k++)
1271 {
1272 SETEXCL(bexcl,i-start,excl->a[k]);
1273 }
1274 }
1275 }
1276 else
1277 {
1278 for(i=start; i<end; i++)
1279 {
1280 for(k=excl->index[i]; k<excl->index[i+1]; k++)
1281 {
1282 RMEXCL(bexcl,i-start,excl->a[k]);
1283 }
1284 }
1285 }
1286 }
1287
1288 int calc_naaj(int icg,int cgtot)
1289 {
1290 int naaj;
1291
1292 if ((cgtot % 2) == 1)
1293 {
1294 /* Odd number of charge groups, easy */
1295 naaj = 1 + (cgtot/2);
1296 }
1297 else if ((cgtot % 4) == 0)
1298 {
1299 /* Multiple of four is hard */
1300 if (icg < cgtot/2)
1301 {
1302 if ((icg % 2) == 0)
1303 {
1304 naaj=1+(cgtot/2);
1305 }
1306 else
1307 {
1308 naaj=cgtot/2;
1309 }
1310 }
1311 else
1312 {
1313 if ((icg % 2) == 1)
1314 {
1315 naaj=1+(cgtot/2);
1316 }
1317 else
1318 {
1319 naaj=cgtot/2;
1320 }
1321 }
1322 }
1323 else
1324 {
1325 /* cgtot/2 = odd */
1326 if ((icg % 2) == 0)
1327 {
1328 naaj=1+(cgtot/2);
1329 }
1330 else
1331 {
1332 naaj=cgtot/2;
1333 }
1334 }
1335 #ifdef DEBUG
1336 fprintf(log,"naaj=%d\n",naaj);
1337 #endif
1338
1339 return naaj;
1340 }
1341
1342 /************************************************
1343 *
1344 * S I M P L E C O R E S T U F F
1345 *
1346 ************************************************/
1347
1348 static real calc_image_tric(rvec xi,rvec xj,matrix box,
1349 rvec b_inv,int *shift)
1350 {
1351 /* This code assumes that the cut-off is smaller than
1352 * a half times the smallest diagonal element of the box.
1353 */
1354 const real h25=2.5;
1355 real dx,dy,dz;
1356 real r2;
1357 int tx,ty,tz;
1358
1359 /* Compute diff vector */
1360 dz = xj[ZZ] - xi[ZZ];
1361 dy = xj[YY] - xi[YY];
1362 dx = xj[XX] - xi[XX];
1363
1364 /* Perform NINT operation, using trunc operation, therefore
1365 * we first add 2.5 then subtract 2 again
1366 */
1367 tz = dz*b_inv[ZZ] + h25;
1368 tz -= 2;
1369 dz -= tz*box[ZZ][ZZ];
1370 dy -= tz*box[ZZ][YY];
1371 dx -= tz*box[ZZ][XX];
1372
1373 ty = dy*b_inv[YY] + h25;
1374 ty -= 2;
1375 dy -= ty*box[YY][YY];
1376 dx -= ty*box[YY][XX];
1377
1378 tx = dx*b_inv[XX]+h25;
1379 tx -= 2;
1380 dx -= tx*box[XX][XX];
1381
1382 /* Distance squared */
1383 r2 = (dx*dx) + (dy*dy) + (dz*dz);
1384
1385 *shift = XYZ2IS(tx,ty,tz);
1386
1387 return r2;
1388 }
1389
1390 static real calc_image_rect(rvec xi,rvec xj,rvec box_size,
1391 rvec b_inv,int *shift)
1392 {
1393 const real h15=1.5;
1394 real ddx,ddy,ddz;
1395 real dx,dy,dz;
1396 real r2;
1397 int tx,ty,tz;
1398
1399 /* Compute diff vector */
1400 dx = xj[XX] - xi[XX];
1401 dy = xj[YY] - xi[YY];
1402 dz = xj[ZZ] - xi[ZZ];
1403
1404 /* Perform NINT operation, using trunc operation, therefore
1405 * we first add 1.5 then subtract 1 again
1406 */
1407 tx = dx*b_inv[XX] + h15;
1408 ty = dy*b_inv[YY] + h15;
1409 tz = dz*b_inv[ZZ] + h15;
1410 tx--;
1411 ty--;
1412 tz--;
1413
1414 /* Correct diff vector for translation */
1415 ddx = tx*box_size[XX] - dx;
1416 ddy = ty*box_size[YY] - dy;
1417 ddz = tz*box_size[ZZ] - dz;
1418
1419 /* Distance squared */
1420 r2 = (ddx*ddx) + (ddy*ddy) + (ddz*ddz);
1421
1422 *shift = XYZ2IS(tx,ty,tz);
1423
1424 return r2;
1425 }
1426
1427 static void add_simple(t_ns_buf *nsbuf,int nrj,atom_id cg_j,
1428 bool bHaveVdW[],int ngid,t_mdatoms *md,
1429 int icg,int jgid,t_block *cgs,t_excl bexcl[],
1430 int shift,t_forcerec *fr,put_in_list_t *put_in_list)
1431 {
1432 if (nsbuf->nj + nrj > MAX_CG)
1433 {
1434 put_in_list(bHaveVdW,ngid,md,icg,jgid,nsbuf->ncg,nsbuf->jcg,
1435 cgs->index,bexcl,shift,fr,FALSE,TRUE,TRUE);
1436 /* Reset buffer contents */
1437 nsbuf->ncg = nsbuf->nj = 0;
1438 }
1439 nsbuf->jcg[nsbuf->ncg++] = cg_j;
1440 nsbuf->nj += nrj;
1441 }
1442
1443 static void ns_inner_tric(rvec x[],int icg,int *i_egp_flags,
1444 int njcg,atom_id jcg[],
1445 matrix box,rvec b_inv,real rcut2,
1446 t_block *cgs,t_ns_buf **ns_buf,
1447 bool bHaveVdW[],int ngid,t_mdatoms *md,
1448 t_excl bexcl[],t_forcerec *fr,
1449 put_in_list_t *put_in_list)
1450 {
1451 int shift;
1452 int j,nrj,jgid;
1453 int *cginfo=fr->cginfo;
1454 atom_id cg_j,*cgindex;
1455 t_ns_buf *nsbuf;
1456
1457 cgindex = cgs->index;
1458 shift = CENTRAL;
1459 for(j=0; (j<njcg); j++)
1460 {
1461 cg_j = jcg[j];
1462 nrj = cgindex[cg_j+1]-cgindex[cg_j];
1463 if (calc_image_tric(x[icg],x[cg_j],box,b_inv,&shift) < rcut2)
1464 {
1465 jgid = GET_CGINFO_GID(cginfo[cg_j]);
1466 if (!(i_egp_flags[jgid] & EGP_EXCL))
1467 {
1468 add_simple(&ns_buf[jgid][shift],nrj,cg_j,
1469 bHaveVdW,ngid,md,icg,jgid,cgs,bexcl,shift,fr,
1470 put_in_list);
1471 }
1472 }
1473 }
1474 }
1475
1476 static void ns_inner_rect(rvec x[],int icg,int *i_egp_flags,
1477 int njcg,atom_id jcg[],
1478 bool bBox,rvec box_size,rvec b_inv,real rcut2,
1479 t_block *cgs,t_ns_buf **ns_buf,
1480 bool bHaveVdW[],int ngid,t_mdatoms *md,
1481 t_excl bexcl[],t_forcerec *fr,
1482 put_in_list_t *put_in_list)
1483 {
1484 int shift;
1485 int j,nrj,jgid;
1486 int *cginfo=fr->cginfo;
1487 atom_id cg_j,*cgindex;
1488 t_ns_buf *nsbuf;
1489
1490 cgindex = cgs->index;
1491 if (bBox)
1492 {
1493 shift = CENTRAL;
1494 for(j=0; (j<njcg); j++)
1495 {
1496 cg_j = jcg[j];
1497 nrj = cgindex[cg_j+1]-cgindex[cg_j];
1498 if (calc_image_rect(x[icg],x[cg_j],box_size,b_inv,&shift) < rcut2)
1499 {
1500 jgid = GET_CGINFO_GID(cginfo[cg_j]);
1501 if (!(i_egp_flags[jgid] & EGP_EXCL))
1502 {
1503 add_simple(&ns_buf[jgid][shift],nrj,cg_j,
1504 bHaveVdW,ngid,md,icg,jgid,cgs,bexcl,shift,fr,
1505 put_in_list);
1506 }
1507 }
1508 }
1509 }
1510 else
1511 {
1512 for(j=0; (j<njcg); j++)
1513 {
1514 cg_j = jcg[j];
1515 nrj = cgindex[cg_j+1]-cgindex[cg_j];
1516 if ((rcut2 == 0) || (distance2(x[icg],x[cg_j]) < rcut2)) {
1517 jgid = GET_CGINFO_GID(cginfo[cg_j]);
1518 if (!(i_egp_flags[jgid] & EGP_EXCL))
1519 {
1520 add_simple(&ns_buf[jgid][CENTRAL],nrj,cg_j,
1521 bHaveVdW,ngid,md,icg,jgid,cgs,bexcl,CENTRAL,fr,
1522 put_in_list);
1523 }
1524 }
1525 }
1526 }
1527 }
1528
1529 /* ns_simple_core needs to be adapted for QMMM still 2005 */
1530
1531 static int ns_simple_core(t_forcerec *fr,
1532 gmx_localtop_t *top,
1533 t_mdatoms *md,
1534 matrix box,rvec box_size,
1535 t_excl bexcl[],atom_id *aaj,
1536 int ngid,t_ns_buf **ns_buf,
1537 put_in_list_t *put_in_list,bool bHaveVdW[])
1538 {
1539 int naaj,k;
1540 real rlist2;
1541 int nsearch,icg,jcg,igid,i0,nri,nn;
1542 int *cginfo;
1543 t_ns_buf *nsbuf;
1544 /* atom_id *i_atoms; */
1545 t_block *cgs=&(top->cgs);
1546 t_blocka *excl=&(top->excls);
1547 rvec b_inv;
1548 int m;
1549 bool bBox,bTriclinic;
1550 int *i_egp_flags;
1551
1552 rlist2 = sqr(fr->rlist);
1553
1554 bBox = (fr->ePBC != epbcNONE);
1555 if (bBox)
1556 {
1557 for(m=0; (m<DIM); m++)
1558 {
1559 b_inv[m] = divide(1.0,box_size[m]);
1560 }
1561 bTriclinic = TRICLINIC(box);
1562 }
1563 else
1564 {
1565 bTriclinic = FALSE;
1566 }
1567
1568 cginfo = fr->cginfo;
1569
1570 nsearch=0;
1571 for (icg=fr->cg0; (icg<fr->hcg); icg++)
1572 {
1573 /*
1574 i0 = cgs->index[icg];
1575 nri = cgs->index[icg+1]-i0;
1576 i_atoms = &(cgs->a[i0]);
1577 i_eg_excl = fr->eg_excl + ngid*md->cENER[*i_atoms];
1578 setexcl(nri,i_atoms,excl,TRUE,bexcl);
1579 */
1580 igid = GET_CGINFO_GID(cginfo[icg]);
1581 i_egp_flags = fr->egp_flags + ngid*igid;
1582 setexcl(cgs->index[icg],cgs->index[icg+1],excl,TRUE,bexcl);
1583
1584 naaj=calc_naaj(icg,cgs->nr);
1585 if (bTriclinic)
1586 {
1587 ns_inner_tric(fr->cg_cm,icg,i_egp_flags,naaj,&(aaj[icg]),
1588 box,b_inv,rlist2,cgs,ns_buf,
1589 bHaveVdW,ngid,md,bexcl,fr,put_in_list);
1590 }
1591 else
1592 {
1593 ns_inner_rect(fr->cg_cm,icg,i_egp_flags,naaj,&(aaj[icg]),
1594 bBox,box_size,b_inv,rlist2,cgs,ns_buf,
1595 bHaveVdW,ngid,md,bexcl,fr,put_in_list);
1596 }
1597 nsearch += naaj;
1598
1599 for(nn=0; (nn<ngid); nn++)
1600 {
1601 for(k=0; (k<SHIFTS); k++)
1602 {
1603 nsbuf = &(ns_buf[nn][k]);
1604 if (nsbuf->ncg > 0)
1605 {
1606 put_in_list(bHaveVdW,ngid,md,icg,nn,nsbuf->ncg,nsbuf->jcg,
1607 cgs->index,bexcl,k,fr,FALSE,TRUE,TRUE);
1608 nsbuf->ncg=nsbuf->nj=0;
1609 }
1610 }
1611 }
1612 /* setexcl(nri,i_atoms,excl,FALSE,bexcl); */
1613 setexcl(cgs->index[icg],cgs->index[icg+1],excl,FALSE,bexcl);
1614 }
1615 close_neighbor_list(fr,FALSE,-1,-1,FALSE);
1616
1617 return nsearch;
1618 }
1619
1620 /************************************************
1621 *
1622 * N S 5 G R I D S T U F F
1623 *
1624 ************************************************/
1625
1626 static inline void get_dx(int Nx,real gridx,real rc2,int xgi,real x,
1627 int *dx0,int *dx1,real *dcx2)
1628 {
1629 real dcx,tmp;
1630 int xgi0,xgi1,i;
1631
1632 if (xgi < 0)
1633 {
1634 *dx0 = 0;
1635 xgi0 = -1;
1636 *dx1 = -1;
1637 xgi1 = 0;
1638 }
1639 else if (xgi >= Nx)
1640 {
1641 *dx0 = Nx;
1642 xgi0 = Nx-1;
1643 *dx1 = Nx-1;
1644 xgi1 = Nx;
1645 }
1646 else
1647 {
1648 dcx2[xgi] = 0;
1649 *dx0 = xgi;
1650 xgi0 = xgi-1;
1651 *dx1 = xgi;
1652 xgi1 = xgi+1;
1653 }
1654
1655 for(i=xgi0; i>=0; i--)
1656 {
1657 dcx = (i+1)*gridx-x;
1658 tmp = dcx*dcx;
1659 if (tmp >= rc2)
1660 break;
1661 *dx0 = i;
1662 dcx2[i] = tmp;
1663 }
1664 for(i=xgi1; i<Nx; i++)
1665 {
1666 dcx = i*gridx-x;
1667 tmp = dcx*dcx;
1668 if (tmp >= rc2)
1669 {
1670 break;
1671 }
1672 *dx1 = i;
1673 dcx2[i] = tmp;
1674 }
1675 }
1676
1677 static inline void get_dx_dd(int Nx,real gridx,real rc2,int xgi,real x,
1678 int ncpddc,int shift_min,int shift_max,
1679 int *g0,int *g1,real *dcx2)
1680 {
1681 real dcx,tmp;
1682 int g_min,g_max,shift_home;
1683
1684 if (xgi < 0)
1685 {
1686 g_min = 0;
1687 g_max = Nx - 1;
1688 *g0 = 0;
1689 *g1 = -1;
1690 }
1691 else if (xgi >= Nx)
1692 {
1693 g_min = 0;
1694 g_max = Nx - 1;
1695 *g0 = Nx;
1696 *g1 = Nx - 1;
1697 }
1698 else
1699 {
1700 if (ncpddc == 0)
1701 {
1702 g_min = 0;
1703 g_max = Nx - 1;
1704 }
1705 else
1706 {
1707 if (xgi < ncpddc)
1708 {
1709 shift_home = 0;
1710 }
1711 else
1712 {
1713 shift_home = -1;
1714 }
1715 g_min = (shift_min == shift_home ? 0 : ncpddc);
1716 g_max = (shift_max == shift_home ? ncpddc - 1 : Nx - 1);
1717 }
1718 if (shift_min > 0)
1719 {
1720 *g0 = g_min;
1721 *g1 = g_min - 1;
1722 }
1723 else if (shift_max < 0)
1724 {
1725 *g0 = g_max + 1;
1726 *g1 = g_max;
1727 }
1728 else
1729 {
1730 *g0 = xgi;
1731 *g1 = xgi;
1732 dcx2[xgi] = 0;
1733 }
1734 }
1735
1736 while (*g0 > g_min)
1737 {
1738 /* Check one grid cell down */
1739 dcx = ((*g0 - 1) + 1)*gridx - x;
1740 tmp = dcx*dcx;
1741 if (tmp >= rc2)
1742 {
1743 break;
1744 }
1745 (*g0)--;
1746 dcx2[*g0] = tmp;
1747 }
1748
1749 while (*g1 < g_max)
1750 {
1751 /* Check one grid cell up */
1752 dcx = (*g1 + 1)*gridx - x;
1753 tmp = dcx*dcx;
1754 if (tmp >= rc2)
1755 {
1756 break;
1757 }
1758 (*g1)++;
1759 dcx2[*g1] = tmp;
1760 }
1761 }
1762
1763
1764 #define sqr(x) ((x)*(x))
1765 #define calc_dx2(XI,YI,ZI,y) (sqr(XI-y[XX]) + sqr(YI-y[YY]) + sqr(ZI-y[ZZ]))
1766 #define calc_cyl_dx2(XI,YI,y) (sqr(XI-y[XX]) + sqr(YI-y[YY]))
1767 /****************************************************
1768 *
1769 * F A S T N E I G H B O R S E A R C H I N G
1770 *
1771 * Optimized neighboursearching routine using grid
1772 * at least 1x1x1, see GROMACS manual
1773 *
1774 ****************************************************/
1775
1776 static void do_longrange(t_commrec *cr,gmx_localtop_t *top,t_forcerec *fr,
1777 int ngid,t_mdatoms *md,int icg,
1778 int jgid,int nlr,
1779 atom_id lr[],t_excl bexcl[],int shift,
1780 rvec x[],rvec box_size,t_nrnb *nrnb,
1781 real lambda,real *dvdlambda,
1782 gmx_grppairener_t *grppener,
1783 bool bDoVdW,bool bDoCoul,
1784 bool bEvaluateNow,put_in_list_t *put_in_list,
1785 bool bHaveVdW[],
1786 bool bDoForces,rvec *f)
1787 {
1788 int n,i;
1789 t_nblist *nl;
1790
1791 for(n=0; n<fr->nnblists; n++)
1792 {
1793 for(i=0; (i<eNL_NR); i++)
1794 {
1795 nl = &fr->nblists[n].nlist_lr[i];
1796 if ((nl->nri > nl->maxnri-32) || bEvaluateNow)
1797 {
1798 close_neighbor_list(fr,TRUE,n,i,FALSE);
1799 /* Evaluate the energies and forces */
1800 do_nonbonded(cr,fr,x,f,md,NULL,
1801 grppener->ener[fr->bBHAM ? egBHAMLR : egLJLR],
1802 grppener->ener[egCOULLR],
1803 grppener->ener[egGB],box_size,
1804 nrnb,lambda,dvdlambda,n,i,
1805 GMX_DONB_LR | GMX_DONB_FORCES);
1806
1807 reset_neighbor_list(fr,TRUE,n,i);
1808 }
1809 }
1810 }
1811
1812 if (!bEvaluateNow)
1813 {
1814 /* Put the long range particles in a list */
1815 /* do_longrange is never called for QMMM */
1816 put_in_list(bHaveVdW,ngid,md,icg,jgid,nlr,lr,top->cgs.index,
1817 bexcl,shift,fr,TRUE,bDoVdW,bDoCoul);
1818 }
1819 }
1820
1821 static void get_cutoff2(t_forcerec *fr,bool bDoLongRange,
1822 real *rvdw2,real *rcoul2,
1823 real *rs2,real *rm2,real *rl2)
1824 {
1825 *rs2 = sqr(fr->rlist);
1826 if (bDoLongRange && fr->bTwinRange)
1827 {
1828 /* The VdW and elec. LR cut-off's could be different,
1829 * so we can not simply set them to rlistlong.
1830 */
1831 if (EVDW_MIGHT_BE_ZERO_AT_CUTOFF(fr->vdwtype) &&
1832 fr->rvdw > fr->rlist)
1833 {
1834 *rvdw2 = sqr(fr->rlistlong);
1835 }
1836 else
1837 {
1838 *rvdw2 = sqr(fr->rvdw);
1839 }
1840 if (EEL_MIGHT_BE_ZERO_AT_CUTOFF(fr->eeltype) &&
1841 fr->rcoulomb > fr->rlist)
1842 {
1843 *rcoul2 = sqr(fr->rlistlong);
1844 }
1845 else
1846 {
1847 *rcoul2 = sqr(fr->rcoulomb);
1848 }
1849 }
1850 else
1851 {
1852 /* Workaround for a gcc -O3 or -ffast-math problem */
1853 *rvdw2 = *rs2;
1854 *rcoul2 = *rs2;
1855 }
1856 *rm2 = min(*rvdw2,*rcoul2);
1857 *rl2 = max(*rvdw2,*rcoul2);
1858 }
1859
1860 static void init_nsgrid_lists(t_forcerec *fr,int ngid,gmx_ns_t *ns)
1861 {
1862 real rvdw2,rcoul2,rs2,rm2,rl2;
1863 int j;
1864
1865 get_cutoff2(fr,TRUE,&rvdw2,&rcoul2,&rs2,&rm2,&rl2);
1866
1867 /* Short range buffers */
1868 snew(ns->nl_sr,ngid);
1869 /* Counters */
1870 snew(ns->nsr,ngid);
1871 snew(ns->nlr_ljc,ngid);
1872 snew(ns->nlr_one,ngid);
1873
1874 if (rm2 > rs2)
1875 {
1876 /* Long range VdW and Coul buffers */
1877 snew(ns->nl_lr_ljc,ngid);
1878 }
1879 if (rl2 > rm2)
1880 {
1881 /* Long range VdW or Coul only buffers */
1882 snew(ns->nl_lr_one,ngid);
1883 }
1884 for(j=0; (j<ngid); j++) {
1885 snew(ns->nl_sr[j],MAX_CG);
1886 if (rm2 > rs2)
1887 {
1888 snew(ns->nl_lr_ljc[j],MAX_CG);
1889 }
1890 if (rl2 > rm2)
1891 {
1892 snew(ns->nl_lr_one[j],MAX_CG);
1893 }
1894 }
1895 if (debug)
1896 {
1897 fprintf(debug,
1898 "ns5_core: rs2 = %g, rm2 = %g, rl2 = %g (nm^2)\n",
1899 rs2,rm2,rl2);
1900 }
1901 }
1902
1903 static int nsgrid_core(FILE *log,t_commrec *cr,t_forcerec *fr,
1904 matrix box,rvec box_size,int ngid,
1905 gmx_localtop_t *top,
1906 t_grid *grid,rvec x[],
1907 t_excl bexcl[],bool *bExcludeAlleg,
1908 t_nrnb *nrnb,t_mdatoms *md,
1909 real lambda,real *dvdlambda,
1910 gmx_grppairener_t *grppener,
1911 put_in_list_t *put_in_list,
1912 bool bHaveVdW[],
1913 bool bDoLongRange,bool bDoForces,rvec *f,
1914 bool bMakeQMMMnblist)
1915 {
1916 gmx_ns_t *ns;
1917 atom_id **nl_lr_ljc,**nl_lr_one,**nl_sr;
1918 int *nlr_ljc,*nlr_one,*nsr;
1919 gmx_domdec_t *dd=NULL;
1920 t_block *cgs=&(top->cgs);
1921 int *cginfo=fr->cginfo;
1922 /* atom_id *i_atoms,*cgsindex=cgs->index; */
1923 ivec sh0,sh1,shp;
1924 int cell_x,cell_y,cell_z;
1925 int d,tx,ty,tz,dx,dy,dz,cj;
1926 #ifdef ALLOW_OFFDIAG_LT_HALFDIAG
1927 int zsh_ty,zsh_tx,ysh_tx;
1928 #endif
1929 int dx0,dx1,dy0,dy1,dz0,dz1;
1930 int Nx,Ny,Nz,shift=-1,j,nrj,nns,nn=-1;
1931 real gridx,gridy,gridz,grid_x,grid_y,grid_z;
1932 real *dcx2,*dcy2,*dcz2;
1933 int zgi,ygi,xgi;
1934 int cg0,cg1,icg=-1,cgsnr,i0,igid,nri,naaj,max_jcg;
1935 int jcg0,jcg1,jjcg,cgj0,jgid;
1936 int *grida,*gridnra,*gridind;
1937 bool rvdw_lt_rcoul,rcoul_lt_rvdw;
1938 rvec xi,*cgcm,grid_offset;
1939 real r2,rs2,rvdw2,rcoul2,rm2,rl2,XI,YI,ZI,dcx,dcy,dcz,tmp1,tmp2;
1940 int *i_egp_flags;
1941 bool bDomDec,bTriclinicX,bTriclinicY;
1942 ivec ncpddc;
1943
1944 ns = &fr->ns;
1945
1946 bDomDec = DOMAINDECOMP(cr);
1947 if (bDomDec)
1948 {
1949 dd = cr->dd;
1950 }
1951
1952 bTriclinicX = ((YY < grid->npbcdim &&
1953 (!bDomDec || dd->nc[YY]==1) && box[YY][XX] != 0) ||
1954 (ZZ < grid->npbcdim &&
1955 (!bDomDec || dd->nc[ZZ]==1) && box[ZZ][XX] != 0));
1956 bTriclinicY = (ZZ < grid->npbcdim &&
1957 (!bDomDec || dd->nc[ZZ]==1) && box[ZZ][YY] != 0);
1958
1959 cgsnr = cgs->nr;
1960
1961 get_cutoff2(fr,bDoLongRange,&rvdw2,&rcoul2,&rs2,&rm2,&rl2);
1962
1963 rvdw_lt_rcoul = (rvdw2 >= rcoul2);
1964 rcoul_lt_rvdw = (rcoul2 >= rvdw2);
1965
1966 if (bMakeQMMMnblist)
1967 {
1968 rm2 = rl2;
1969 rs2 = rl2;
1970 }
1971
1972 nl_sr = ns->nl_sr;
1973 nsr = ns->nsr;
1974 nl_lr_ljc = ns->nl_lr_ljc;
1975 nl_lr_one = ns->nl_lr_one;
1976 nlr_ljc = ns->nlr_ljc;
1977 nlr_one = ns->nlr_one;
1978
1979 /* Unpack arrays */
1980 cgcm = fr->cg_cm;
1981 Nx = grid->n[XX];
1982 Ny = grid->n[YY];
1983 Nz = grid->n[ZZ];
1984 grida = grid->a;
1985 gridind = grid->index;
1986 gridnra = grid->nra;
1987 nns = 0;
1988
1989 gridx = grid->cell_size[XX];
1990 gridy = grid->cell_size[YY];
1991 gridz = grid->cell_size[ZZ];
1992 grid_x = 1/gridx;
1993 grid_y = 1/gridy;
1994 grid_z = 1/gridz;
1995 copy_rvec(grid->cell_offset,grid_offset);
1996 copy_ivec(grid->ncpddc,ncpddc);
1997 dcx2 = grid->dcx2;
1998 dcy2 = grid->dcy2;
1999 dcz2 = grid->dcz2;
2000
2001 #ifdef ALLOW_OFFDIAG_LT_HALFDIAG
2002 zsh_ty = floor(-box[ZZ][YY]/box[YY][YY]+0.5);
2003 zsh_tx = floor(-box[ZZ][XX]/box[XX][XX]+0.5);
2004 ysh_tx = floor(-box[YY][XX]/box[XX][XX]+0.5);
2005 if (zsh_tx!=0 && ysh_tx!=0)
2006 {
2007 /* This could happen due to rounding, when both ratios are 0.5 */
2008 ysh_tx = 0;
2009 }
2010 #endif
2011
2012 debug_gmx();
2013
2014 if (fr->n_tpi)
2015 {
2016 /* We only want a list for the test particle */
2017 cg0 = cgsnr - 1;
2018 }
2019 else
2020 {
2021 cg0 = grid->icg0;
2022 }
2023 cg1 = grid->icg1;
2024
2025 /* Set the shift range */
2026 for(d=0; d<DIM; d++)
2027 {
2028 sh0[d] = -1;
2029 sh1[d] = 1;
2030 /* Check if we need periodicity shifts.
2031 * Without PBC or with domain decomposition we don't need them.
2032 */
2033 if (d >= ePBC2npbcdim(fr->ePBC) || (bDomDec && dd->nc[d] > 1))
2034 {
2035 shp[d] = 0;
2036 }
2037 else
2038 {
2039 if (d == XX &&
2040 box[XX][XX] - fabs(box[YY][XX]) - fabs(box[ZZ][XX]) < sqrt(rl2))
2041 {
2042 shp[d] = 2;
2043 }
2044 else
2045 {
2046 shp[d] = 1;
2047 }
2048 }
2049 }
2050
2051 /* Loop over charge groups */
2052 for(icg=cg0; (icg < cg1); icg++)
2053 {
2054 igid = GET_CGINFO_GID(cginfo[icg]);
2055 /* Skip this charge group if all energy groups are excluded! */
2056 if (bExcludeAlleg[igid])
2057 {
2058 continue;
2059 }
2060
2061 i0 = cgs->index[icg];
2062
2063 if (bMakeQMMMnblist)
2064 {
2065 /* Skip this charge group if it is not a QM atom while making a
2066 * QM/MM neighbourlist
2067 */
2068 if (md->bQM[i0]==FALSE)
2069 {
2070 continue; /* MM particle, go to next particle */
2071 }
2072
2073 /* Compute the number of charge groups that fall within the control
2074 * of this one (icg)
2075 */
2076 naaj = calc_naaj(icg,cgsnr);
2077 jcg0 = icg;
2078 jcg1 = icg + naaj;
2079 max_jcg = cgsnr;
2080 }
2081 else
2082 {
2083 /* make a normal neighbourlist */
2084
2085 if (bDomDec)
2086 {
2087 /* Get the j charge-group and dd cell shift ranges */
2088 dd_get_ns_ranges(cr->dd,icg,&jcg0,&jcg1,sh0,sh1);
2089 max_jcg = 0;
2090 }
2091 else
2092 {
2093 /* Compute the number of charge groups that fall within the control
2094 * of this one (icg)
2095 */
2096 naaj = calc_naaj(icg,cgsnr);
2097 jcg0 = icg;
2098 jcg1 = icg + naaj;
2099
2100 if (fr->n_tpi)
2101 {
2102 /* The i-particle is awlways the test particle,
2103 * so we want all j-particles
2104 */
2105 max_jcg = cgsnr - 1;
2106 }
2107 else
2108 {
2109 max_jcg = jcg1 - cgsnr;
2110 }
2111 }
2112 }
2113
2114 i_egp_flags = fr->egp_flags + igid*ngid;
2115
2116 /* Set the exclusions for the atoms in charge group icg using a bitmask */
2117 setexcl(i0,cgs->index[icg+1],&top->excls,TRUE,bexcl);
2118
2119 ci2xyz(grid,icg,&cell_x,&cell_y,&cell_z);
2120
2121 /* Changed iicg to icg, DvdS 990115
2122 * (but see consistency check above, DvdS 990330)
2123 */
2124 #ifdef NS5DB
2125 fprintf(log,"icg=%5d, naaj=%5d, cell %d %d %d\n",
2126 icg,naaj,cell_x,cell_y,cell_z);
2127 #endif
2128 /* Loop over shift vectors in three dimensions */
2129 for (tz=-shp[ZZ]; tz<=shp[ZZ]; tz++)
2130 {
2131 ZI = cgcm[icg][ZZ]+tz*box[ZZ][ZZ];
2132 /* Calculate range of cells in Z direction that have the shift tz */
2133 zgi = cell_z + tz*Nz;
2134 #define FAST_DD_NS
2135 #ifndef FAST_DD_NS
2136 get_dx(Nz,gridz,rl2,zgi,ZI,&dz0,&dz1,dcz2);
2137 #else
2138 get_dx_dd(Nz,gridz,rl2,zgi,ZI-grid_offset[ZZ],
2139 ncpddc[ZZ],sh0[ZZ],sh1[ZZ],&dz0,&dz1,dcz2);
2140 #endif
2141 if (dz0 > dz1)
2142 {
2143 continue;
2144 }
2145 for (ty=-shp[YY]; ty<=shp[YY]; ty++)
2146 {
2147 YI = cgcm[icg][YY]+ty*box[YY][YY]+tz*box[ZZ][YY];
2148 /* Calculate range of cells in Y direction that have the shift ty */
2149 if (bTriclinicY)
2150 {
2151 ygi = (int)(Ny + (YI - grid_offset[YY])*grid_y) - Ny;
2152 }
2153 else
2154 {
2155 ygi = cell_y + ty*Ny;
2156 }
2157 #ifndef FAST_DD_NS
2158 get_dx(Ny,gridy,rl2,ygi,YI,&dy0,&dy1,dcy2);
2159 #else
2160 get_dx_dd(Ny,gridy,rl2,ygi,YI-grid_offset[YY],
2161 ncpddc[YY],sh0[YY],sh1[YY],&dy0,&dy1,dcy2);
2162 #endif
2163 if (dy0 > dy1)
2164 {
2165 continue;
2166 }
2167 for (tx=-shp[XX]; tx<=shp[XX]; tx++)
2168 {
2169 XI = cgcm[icg][XX]+tx*box[XX][XX]+ty*box[YY][XX]+tz*box[ZZ][XX];
2170 /* Calculate range of cells in X direction that have the shift tx */
2171 if (bTriclinicX)
2172 {
2173 xgi = (int)(Nx + (XI - grid_offset[XX])*grid_x) - Nx;
2174 }
2175 else
2176 {
2177 xgi = cell_x + tx*Nx;
2178 }
2179 #ifndef FAST_DD_NS
2180 get_dx(Nx,gridx,rl2,xgi*Nx,XI,&dx0,&dx1,dcx2);
2181 #else
2182 get_dx_dd(Nx,gridx,rl2,xgi,XI-grid_offset[XX],
2183 ncpddc[XX],sh0[XX],sh1[XX],&dx0,&dx1,dcx2);
2184 #endif
2185 if (dx0 > dx1)
2186 {
2187 continue;
2188 }
2189 /* Adress: an explicit cg that has a weigthing function of 0 is excluded
2190 * from the neigbour list as it will not interact */
2191 if (fr->adress_type != eAdressOff){
2192 if (md->wf[cgs->index[icg]]==0 && egp_explicit(fr, igid)){
2193 continue;
2194 }
2195 }
2196 /* Get shift vector */
2197 shift=XYZ2IS(tx,ty,tz);
2198 #ifdef NS5DB
2199 range_check(shift,0,SHIFTS);
2200 #endif
2201 for(nn=0; (nn<ngid); nn++)
2202 {
2203 nsr[nn] = 0;
2204 nlr_ljc[nn] = 0;
2205 nlr_one[nn] = 0;
2206 }
2207 #ifdef NS5DB
2208 fprintf(log,"shift: %2d, dx0,1: %2d,%2d, dy0,1: %2d,%2d, dz0,1: %2d,%2d\n",
2209 shift,dx0,dx1,dy0,dy1,dz0,dz1);
2210 fprintf(log,"cgcm: %8.3f %8.3f %8.3f\n",cgcm[icg][XX],
2211 cgcm[icg][YY],cgcm[icg][ZZ]);
2212 fprintf(log,"xi: %8.3f %8.3f %8.3f\n",XI,YI,ZI);
2213 #endif
2214 for (dx=dx0; (dx<=dx1); dx++)
2215 {
2216 tmp1 = rl2 - dcx2[dx];
2217 for (dy=dy0; (dy<=dy1); dy++)
2218 {
2219 tmp2 = tmp1 - dcy2[dy];
2220 if (tmp2 > 0)
2221 {
2222 for (dz=dz0; (dz<=dz1); dz++) {
2223 if (tmp2 > dcz2[dz]) {
2224 /* Find grid-cell cj in which possible neighbours are */
2225 cj = xyz2ci(Ny,Nz,dx,dy,dz);
2226
2227 /* Check out how many cgs (nrj) there in this cell */
2228 nrj = gridnra[cj];
2229
2230 /* Find the offset in the cg list */
2231 cgj0 = gridind[cj];
2232
2233 /* Check if all j's are out of range so we
2234 * can skip the whole cell.
2235 * Should save some time, especially with DD.
2236 */
2237 if (nrj == 0 ||
2238 (grida[cgj0] >= max_jcg &&
2239 (grida[cgj0] >= jcg1 || grida[cgj0+nrj-1] < jcg0)))
2240 {
2241 continue;
2242 }
2243
2244 /* Loop over cgs */
2245 for (j=0; (j<nrj); j++)
2246 {
2247 jjcg = grida[cgj0+j];
2248
2249 /* check whether this guy is in range! */
2250 if ((jjcg >= jcg0 && jjcg < jcg1) ||
2251 (jjcg < max_jcg))
2252 {
2253 r2=calc_dx2(XI,YI,ZI,cgcm[jjcg]);
2254 if (r2 < rl2) {
2255 /* jgid = gid[cgsatoms[cgsindex[jjcg]]]; */
2256 jgid = GET_CGINFO_GID(cginfo[jjcg]);
2257 /* check energy group exclusions */
2258 if (!(i_egp_flags[jgid] & EGP_EXCL))
2259 {
2260 if (r2 < rs2)
2261 {
2262 if (nsr[jgid] >= MAX_CG)
2263 {
2264 put_in_list(bHaveVdW,ngid,md,icg,jgid,
2265 nsr[jgid],nl_sr[jgid],
2266 cgs->index,/* cgsatoms, */ bexcl,
2267 shift,fr,FALSE,TRUE,TRUE);
2268 nsr[jgid]=0;
2269 }
2270 nl_sr[jgid][nsr[jgid]++]=jjcg;
2271 }
2272 else if (r2 < rm2)
2273 {
2274 if (nlr_ljc[jgid] >= MAX_CG)
2275 {
2276 do_longrange(cr,top,fr,ngid,md,icg,jgid,
2277 nlr_ljc[jgid],
2278 nl_lr_ljc[jgid],bexcl,shift,x,
2279 box_size,nrnb,
2280 lambda,dvdlambda,
2281 grppener,
2282 TRUE,TRUE,FALSE,
2283 put_in_list,
2284 bHaveVdW,
2285 bDoForces,f);
2286 nlr_ljc[jgid]=0;
2287 }
2288 nl_lr_ljc[jgid][nlr_ljc[jgid]++]=jjcg;
2289 }
2290 else
2291 {
2292 if (nlr_one[jgid] >= MAX_CG) {
2293 do_longrange(cr,top,fr,ngid,md,icg,jgid,
2294 nlr_one[jgid],
2295 nl_lr_one[jgid],bexcl,shift,x,
2296 box_size,nrnb,
2297 lambda,dvdlambda,
2298 grppener,
2299 rvdw_lt_rcoul,rcoul_lt_rvdw,FALSE,
2300 put_in_list,
2301 bHaveVdW,
2302 bDoForces,f);
2303 nlr_one[jgid]=0;
2304 }
2305 nl_lr_one[jgid][nlr_one[jgid]++]=jjcg;
2306 }
2307 }
2308 }
2309 nns++;
2310 }
2311 }
2312 }
2313 }
2314 }
2315 }
2316 }
2317 /* CHECK whether there is anything left in the buffers */
2318 for(nn=0; (nn<ngid); nn++)
2319 {
2320 if (nsr[nn] > 0)
2321 {
2322 put_in_list(bHaveVdW,ngid,md,icg,nn,nsr[nn],nl_sr[nn],
2323 cgs->index, /* cgsatoms, */ bexcl,
2324 shift,fr,FALSE,TRUE,TRUE);
2325 }
2326
2327 if (nlr_ljc[nn] > 0)
2328 {
2329 do_longrange(cr,top,fr,ngid,md,icg,nn,nlr_ljc[nn],
2330 nl_lr_ljc[nn],bexcl,shift,x,box_size,nrnb,
2331 lambda,dvdlambda,grppener,TRUE,TRUE,FALSE,
2332 put_in_list,bHaveVdW,bDoForces,f);
2333 }
2334
2335 if (nlr_one[nn] > 0)
2336 {
2337 do_longrange(cr,top,fr,ngid,md,icg,nn,nlr_one[nn],
2338 nl_lr_one[nn],bexcl,shift,x,box_size,nrnb,
2339 lambda,dvdlambda,grppener,
2340 rvdw_lt_rcoul,rcoul_lt_rvdw,FALSE,
2341 put_in_list,bHaveVdW,bDoForces,f);
2342 }
2343 }
2344 }
2345 }
2346 }
2347 /* setexcl(nri,i_atoms,&top->atoms.excl,FALSE,bexcl); */
2348 setexcl(cgs->index[icg],cgs->index[icg+1],&top->excls,FALSE,bexcl);
2349 }
2350 /* Perform any left over force calculations */
2351 for (nn=0; (nn<ngid); nn++)
2352 {
2353 if (rm2 > rs2)
2354 {
2355 do_longrange(cr,top,fr,0,md,icg,nn,nlr_ljc[nn],
2356 nl_lr_ljc[nn],bexcl,shift,x,box_size,nrnb,
2357 lambda,dvdlambda,grppener,
2358 TRUE,TRUE,TRUE,put_in_list,bHaveVdW,bDoForces,f);
2359 }
2360 if (rl2 > rm2) {
2361 do_longrange(cr,top,fr,0,md,icg,nn,nlr_one[nn],
2362 nl_lr_one[nn],bexcl,shift,x,box_size,nrnb,
2363 lambda,dvdlambda,grppener,
2364 rvdw_lt_rcoul,rcoul_lt_rvdw,
2365 TRUE,put_in_list,bHaveVdW,bDoForces,f);
2366 }
2367 }
2368 debug_gmx();
2369
2370 /* Close off short range neighbourlists */
2371 close_neighbor_list(fr,FALSE,-1,-1,bMakeQMMMnblist);
2372
2373 return nns;
2374 }
2375
2376 void ns_realloc_natoms(gmx_ns_t *ns,int natoms)
2377 {
2378 int i;
2379
2380 if (natoms > ns->nra_alloc)
2381 {
2382 ns->nra_alloc = over_alloc_dd(natoms);
2383 srenew(ns->bexcl,ns->nra_alloc);
2384 for(i=0; i<ns->nra_alloc; i++)
2385 {
2386 ns->bexcl[i] = 0;
2387 }
2388 }
2389 }
2390
2391 void init_ns(FILE *fplog,const t_commrec *cr,
2392 gmx_ns_t *ns,t_forcerec *fr,
2393 const gmx_mtop_t *mtop,
2394 matrix box)
2395 {
2396 int mt,icg,nr_in_cg,maxcg,i,j,jcg,ngid,ncg;
2397 t_block *cgs;
2398 char *ptr;
2399
2400 /* Compute largest charge groups size (# atoms) */
2401 nr_in_cg=1;
2402 for(mt=0; mt<mtop->nmoltype; mt++) {
2403 cgs = &mtop->moltype[mt].cgs;
2404 for (icg=0; (icg < cgs->nr); icg++)
2405 {
2406 nr_in_cg=max(nr_in_cg,(int)(cgs->index[icg+1]-cgs->index[icg]));
2407 }
2408 }
2409
2410 /* Verify whether largest charge group is <= max cg.
2411 * This is determined by the type of the local exclusion type
2412 * Exclusions are stored in bits. (If the type is not large
2413 * enough, enlarge it, unsigned char -> unsigned short -> unsigned long)
2414 */
2415 maxcg = sizeof(t_excl)*8;
2416 if (nr_in_cg > maxcg)
2417 {
2418 gmx_fatal(FARGS,"Max #atoms in a charge group: %d > %d\n",
2419 nr_in_cg,maxcg);
2420 }
2421
2422 ngid = mtop->groups.grps[egcENER].nr;
2423 snew(ns->bExcludeAlleg,ngid);
2424 for(i=0; i<ngid; i++) {
2425 ns->bExcludeAlleg[i] = TRUE;
2426 for(j=0; j<ngid; j++)
2427 {
2428 if (!(fr->egp_flags[i*ngid+j] & EGP_EXCL))
2429 {
2430 ns->bExcludeAlleg[i] = FALSE;
2431 }
2432 }
2433 }
2434
2435 if (fr->bGrid) {
2436 /* Grid search */
2437 ns->grid = init_grid(fplog,fr);
2438 init_nsgrid_lists(fr,ngid,ns);
2439 }
2440 else
2441 {
2442 /* Simple search */
2443 snew(ns->ns_buf,ngid);
2444 for(i=0; (i<ngid); i++)
2445 {
2446 snew(ns->ns_buf[i],SHIFTS);
2447 }
2448 ncg = ncg_mtop(mtop);
2449 snew(ns->simple_aaj,2*ncg);
2450 for(jcg=0; (jcg<ncg); jcg++)
2451 {
2452 ns->simple_aaj[jcg] = jcg;
2453 ns->simple_aaj[jcg+ncg] = jcg;
2454 }
2455 }
2456
2457 /* Create array that determines whether or not atoms have VdW */
2458 snew(ns->bHaveVdW,fr->ntype);
2459 for(i=0; (i<fr->ntype); i++)
2460 {
2461 for(j=0; (j<fr->ntype); j++)
2462 {
2463 ns->bHaveVdW[i] = (ns->bHaveVdW[i] ||
2464 (fr->bBHAM ?
2465 ((BHAMA(fr->nbfp,fr->ntype,i,j) != 0) ||
2466 (BHAMB(fr->nbfp,fr->ntype,i,j) != 0) ||
2467 (BHAMC(fr->nbfp,fr->ntype,i,j) != 0)) :
2468 ((C6(fr->nbfp,fr->ntype,i,j) != 0) ||
2469 (C12(fr->nbfp,fr->ntype,i,j) != 0))));
2470 }
2471 }
2472 if (debug)
2473 pr_bvec(debug,0,"bHaveVdW",ns->bHaveVdW,fr->ntype,TRUE);
2474
2475 ns->nra_alloc = 0;
2476 ns->bexcl = NULL;
2477 if (!DOMAINDECOMP(cr))
2478 {
2479 /* This could be reduced with particle decomposition */
2480 ns_realloc_natoms(ns,mtop->natoms);
2481 }
2482
2483 ns->nblist_initialized=FALSE;
2484
2485 /* nbr list debug dump */
2486 {
2487 char *ptr=getenv("GMX_DUMP_NL");
2488 if (ptr)
2489 {
2490 ns->dump_nl=strtol(ptr,NULL,10);
2491 if (fplog)
2492 {
2493 fprintf(fplog, "GMX_DUMP_NL = %d", ns->dump_nl);
2494 }
2495 }
2496 else
2497 {
2498 ns->dump_nl=0;
2499 }
2500 }
2501 }
2502
2503
2504 int search_neighbours(FILE *log,t_forcerec *fr,
2505 rvec x[],matrix box,
2506 gmx_localtop_t *top,
2507 gmx_groups_t *groups,
2508 t_commrec *cr,
2509 t_nrnb *nrnb,t_mdatoms *md,
2510 real lambda,real *dvdlambda,
2511 gmx_grppairener_t *grppener,
2512 bool bFillGrid,
2513 bool bDoLongRange,
2514 bool bDoForces,rvec *f)
2515 {
2516 t_block *cgs=&(top->cgs);
2517 rvec box_size,grid_x0,grid_x1;
2518 int i,j,m,ngid;
2519 real min_size,grid_dens;
2520 int nsearch;
2521 bool bGrid;
2522 char *ptr;
2523 bool *i_egp_flags;
2524 int cg_start,cg_end,start,end;
2525 gmx_ns_t *ns;
2526 t_grid *grid;
2527 gmx_domdec_zones_t *dd_zones;
2528 put_in_list_t *put_in_list;
2529
2530 ns = &fr->ns;
2531
2532 /* Set some local variables */
2533 bGrid = fr->bGrid;
2534 ngid = groups->grps[egcENER].nr;
2535
2536 for(m=0; (m<DIM); m++)
2537 {
2538 box_size[m] = box[m][m];
2539 }
2540
2541 if (fr->ePBC != epbcNONE)
2542 {
2543 if (sqr(fr->rlistlong) >= max_cutoff2(fr->ePBC,box))
2544 {
2545 gmx_fatal(FARGS,"One of the box vectors has become shorter than twice the cut-off length or box_yy-|box_zy| or box_zz has become smaller than the cut-off.");
2546 }
2547 if (!bGrid)
2548 {
2549 min_size = min(box_size[XX],min(box_size[YY],box_size[ZZ]));
2550 if (2*fr->rlistlong >= min_size)
2551 gmx_fatal(FARGS,"One of the box diagonal elements has become smaller than twice the cut-off length.");
2552 }
2553 }
2554
2555 if (DOMAINDECOMP(cr))
2556 {
2557 ns_realloc_natoms(ns,cgs->index[cgs->nr]);
2558 }
2559 debug_gmx();
2560
2561 /* Reset the neighbourlists */
2562 reset_neighbor_list(fr,FALSE,-1,-1);
2563
2564 if (bGrid && bFillGrid)
2565 {
2566
2567 grid = ns->grid;
2568 if (DOMAINDECOMP(cr))
2569 {
2570 dd_zones = domdec_zones(cr->dd);
2571 }
2572 else
2573 {
2574 dd_zones = NULL;
2575
2576 get_nsgrid_boundaries(grid,NULL,box,NULL,NULL,NULL,
2577 cgs->nr,fr->cg_cm,grid_x0,grid_x1,&grid_dens);
2578
2579 grid_first(log,grid,NULL,NULL,fr->ePBC,box,grid_x0,grid_x1,
2580 fr->rlistlong,grid_dens);
2581 }
2582 debug_gmx();
2583
2584 /* Don't know why this all is... (DvdS 3/99) */
2585 #ifndef SEGV
2586 start = 0;
2587 end = cgs->nr;
2588 #else
2589 start = fr->cg0;
2590 end = (cgs->nr+1)/2;
2591 #endif
2592
2593 if (DOMAINDECOMP(cr))
2594 {
2595 end = cgs->nr;
2596 fill_grid(log,dd_zones,grid,end,-1,end,fr->cg_cm);
2597 grid->icg0 = 0;
2598 grid->icg1 = dd_zones->izone[dd_zones->nizone-1].cg1;
2599 }
2600 else
2601 {
2602 fill_grid(log,NULL,grid,cgs->nr,fr->cg0,fr->hcg,fr->cg_cm);
2603 grid->icg0 = fr->cg0;
2604 grid->icg1 = fr->hcg;
2605 debug_gmx();
2606
2607 if (PARTDECOMP(cr))
2608 mv_grid(cr,grid);
2609 debug_gmx();
2610 }
2611
2612 calc_elemnr(log,grid,start,end,cgs->nr);
2613 calc_ptrs(grid);
2614 grid_last(log,grid,start,end,cgs->nr);
2615
2616 if (gmx_debug_at)
2617 {
2618 check_grid(debug,grid);
2619 print_grid(debug,grid);
2620 }
2621 }
2622 else if (fr->n_tpi)
2623 {
2624 /* Set the grid cell index for the test particle only.
2625 * The cell to cg index is not corrected, but that does not matter.
2626 */
2627 fill_grid(log,NULL,ns->grid,fr->hcg,fr->hcg-1,fr->hcg,fr->cg_cm);
2628 }
2629 debug_gmx();
2630
2631 if (!fr->ns.bCGlist)
2632 {
2633 put_in_list = put_in_list_at;
2634 }
2635 else
2636 {
2637 put_in_list = put_in_list_cg;
2638 }
2639
2640 /* Do the core! */
2641 if (bGrid)
2642 {
2643 grid = ns->grid;
2644 nsearch = nsgrid_core(log,cr,fr,box,box_size,ngid,top,
2645 grid,x,ns->bexcl,ns->bExcludeAlleg,
2646 nrnb,md,lambda,dvdlambda,grppener,
2647 put_in_list,ns->bHaveVdW,
2648 bDoLongRange,bDoForces,f,
2649 FALSE);
2650
2651 /* neighbour searching withouth QMMM! QM atoms have zero charge in
2652 * the classical calculation. The charge-charge interaction
2653 * between QM and MM atoms is handled in the QMMM core calculation
2654 * (see QMMM.c). The VDW however, we'd like to compute classically
2655 * and the QM MM atom pairs have just been put in the
2656 * corresponding neighbourlists. in case of QMMM we still need to
2657 * fill a special QMMM neighbourlist that contains all neighbours
2658 * of the QM atoms. If bQMMM is true, this list will now be made:
2659 */
2660 if (fr->bQMMM && fr->qr->QMMMscheme!=eQMMMschemeoniom)
2661 {
2662 nsearch += nsgrid_core(log,cr,fr,box,box_size,ngid,top,
2663 grid,x,ns->bexcl,ns->bExcludeAlleg,
2664 nrnb,md,lambda,dvdlambda,grppener,
2665 put_in_list_qmmm,ns->bHaveVdW,
2666 bDoLongRange,bDoForces,f,
2667 TRUE);
2668 }
2669 }
2670 else
2671 {
2672 nsearch = ns_simple_core(fr,top,md,box,box_size,
2673 ns->bexcl,ns->simple_aaj,
2674 ngid,ns->ns_buf,put_in_list,ns->bHaveVdW);
2675 }
2676 debug_gmx();
2677
2678 #ifdef DEBUG
2679 pr_nsblock(log);
2680 #endif
2681
2682 inc_nrnb(nrnb,eNR_NS,nsearch);
2683 /* inc_nrnb(nrnb,eNR_LR,fr->nlr); */
2684
2685 return nsearch;
2686 }
2687
2688 int natoms_beyond_ns_buffer(t_inputrec *ir,t_forcerec *fr,t_block *cgs,
2689 matrix scale_tot,rvec *x)
2690 {
2691 int cg0,cg1,cg,a0,a1,a,i,j;
2692 real rint,hbuf2,scale;
2693 rvec *cg_cm,cgsc;
2694 bool bIsotropic;
2695 int nBeyond;
2696
2697 nBeyond = 0;
2698
2699 rint = max(ir->rcoulomb,ir->rvdw);
2700 if (ir->rlist < rint)
2701 {
2702 gmx_fatal(FARGS,"The neighbor search buffer has negative size: %f nm",
2703 ir->rlist - rint);
2704 }
2705 cg_cm = fr->cg_cm;
2706
2707 cg0 = fr->cg0;
2708 cg1 = fr->hcg;
2709
2710 if (!EI_DYNAMICS(ir->eI) || !DYNAMIC_BOX(*ir))
2711 {
2712 hbuf2 = sqr(0.5*(ir->rlist - rint));
2713 for(cg=cg0; cg<cg1; cg++)
2714 {
2715 a0 = cgs->index[cg];
2716 a1 = cgs->index[cg+1];
2717 for(a=a0; a<a1; a++)
2718 {
2719 if (distance2(cg_cm[cg],x[a]) > hbuf2)
2720 {
2721 nBeyond++;
2722 }
2723 }
2724 }
2725 }
2726 else
2727 {
2728 bIsotropic = TRUE;
2729 scale = scale_tot[0][0];
2730 for(i=1; i<DIM; i++)
2731 {
2732 /* With anisotropic scaling, the original spherical ns volumes become
2733 * ellipsoids. To avoid costly transformations we use the minimum
2734 * eigenvalue of the scaling matrix for determining the buffer size.
2735 * Since the lower half is 0, the eigenvalues are the diagonal elements.
2736 */
2737 scale = min(scale,scale_tot[i][i]);
2738 if (scale_tot[i][i] != scale_tot[i-1][i-1])
2739 {
2740 bIsotropic = FALSE;
2741 }
2742 for(j=0; j<i; j++)
2743 {
2744 if (scale_tot[i][j] != 0)
2745 {
2746 bIsotropic = FALSE;
2747 }
2748 }
2749 }
2750 hbuf2 = sqr(0.5*(scale*ir->rlist - rint));
2751 if (bIsotropic)
2752 {
2753 for(cg=cg0; cg<cg1; cg++)
2754 {
2755 svmul(scale,cg_cm[cg],cgsc);
2756 a0 = cgs->index[cg];
2757 a1 = cgs->index[cg+1];
2758 for(a=a0; a<a1; a++)
2759 {
2760 if (distance2(cgsc,x[a]) > hbuf2)
2761 {
2762 nBeyond++;
2763 }
2764 }
2765 }
2766 }
2767 else
2768 {
2769 /* Anistropic scaling */
2770 for(cg=cg0; cg<cg1; cg++)
2771 {
2772 /* Since scale_tot contains the transpose of the scaling matrix,
2773 * we need to multiply with the transpose.
2774 */
2775 tmvmul_ur0(scale_tot,cg_cm[cg],cgsc);
2776 a0 = cgs->index[cg];
2777 a1 = cgs->index[cg+1];
2778 for(a=a0; a<a1; a++)
2779 {
2780 if (distance2(cgsc,x[a]) > hbuf2)
2781 {
2782 nBeyond++;
2783 }
2784 }
2785 }
2786 }
2787 }
2788
2789 return nBeyond;
2790 }
got merge a long time ago