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