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