2 * This file is part of the GROMACS molecular simulation package.
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,2017, 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.
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.
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.
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.
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.
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.
49 #include "gromacs/domdec/domdec.h"
50 #include "gromacs/domdec/domdec_struct.h"
51 #include "gromacs/gmxlib/network.h"
52 #include "gromacs/gmxlib/nrnb.h"
53 #include "gromacs/gmxlib/nonbonded/nonbonded.h"
54 #include "gromacs/math/functions.h"
55 #include "gromacs/math/utilities.h"
56 #include "gromacs/math/vec.h"
57 #include "gromacs/math/vecdump.h"
58 #include "gromacs/mdlib/force.h"
59 #include "gromacs/mdlib/nsgrid.h"
60 #include "gromacs/mdlib/qmmm.h"
61 #include "gromacs/mdtypes/commrec.h"
62 #include "gromacs/mdtypes/group.h"
63 #include "gromacs/mdtypes/inputrec.h"
64 #include "gromacs/mdtypes/md_enums.h"
65 #include "gromacs/pbcutil/ishift.h"
66 #include "gromacs/pbcutil/pbc.h"
67 #include "gromacs/topology/mtop_util.h"
68 #include "gromacs/utility/fatalerror.h"
69 #include "gromacs/utility/smalloc.h"
72 * E X C L U S I O N H A N D L I N G
76 static void SETEXCL_(t_excl e
[], int i
, int j
)
80 static void RMEXCL_(t_excl e
[], int i
, int j
)
82 e
[j
] = e
[j
] & ~(1<<i
);
84 static gmx_bool
ISEXCL_(t_excl e
[], int i
, int j
)
86 return (gmx_bool
)(e
[j
] & (1<<i
));
88 static gmx_bool
NOTEXCL_(t_excl e
[], int i
, int j
)
90 return !(ISEXCL(e
, i
, j
));
93 #define SETEXCL(e, i, j) (e)[((int) (j))] |= (1<<((int) (i)))
94 #define RMEXCL(e, i, j) (e)[((int) (j))] &= (~(1<<((int) (i))))
95 #define ISEXCL(e, i, j) (gmx_bool) ((e)[((int) (j))] & (1<<((int) (i))))
96 #define NOTEXCL(e, i, j) !(ISEXCL(e, i, j))
100 round_up_to_simd_width(int length
, int simd_width
)
104 offset
= (simd_width
> 0) ? length
% simd_width
: 0;
106 return (offset
== 0) ? length
: length
-offset
+simd_width
;
108 /************************************************
110 * U T I L I T I E S F O R N S
112 ************************************************/
114 void reallocate_nblist(t_nblist
*nl
)
118 fprintf(debug
, "reallocating neigborlist (ielec=%d, ivdw=%d, igeometry=%d, type=%d), maxnri=%d\n",
119 nl
->ielec
, nl
->ivdw
, nl
->igeometry
, nl
->type
, nl
->maxnri
);
121 srenew(nl
->iinr
, nl
->maxnri
);
122 if (nl
->igeometry
== GMX_NBLIST_GEOMETRY_CG_CG
)
124 srenew(nl
->iinr_end
, nl
->maxnri
);
126 srenew(nl
->gid
, nl
->maxnri
);
127 srenew(nl
->shift
, nl
->maxnri
);
128 srenew(nl
->jindex
, nl
->maxnri
+1);
132 static void init_nblist(FILE *log
, t_nblist
*nl_sr
,
134 int ivdw
, int ivdwmod
,
135 int ielec
, int ielecmod
,
136 int igeometry
, int type
,
137 gmx_bool bElecAndVdwSwitchDiffers
)
152 /* Set coul/vdw in neighborlist, and for the normal loops we determine
153 * an index of which one to call.
156 nl
->ivdwmod
= ivdwmod
;
158 nl
->ielecmod
= ielecmod
;
160 nl
->igeometry
= igeometry
;
162 if (nl
->type
== GMX_NBLIST_INTERACTION_FREE_ENERGY
)
164 nl
->igeometry
= GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE
;
167 /* This will also set the simd_padding_width field */
168 gmx_nonbonded_set_kernel_pointers(log
, nl
, bElecAndVdwSwitchDiffers
);
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.
176 nl
->maxnri
= homenr
*4;
183 nl
->jindex
= nullptr;
185 nl
->excl_fep
= nullptr;
186 reallocate_nblist(nl
);
191 fprintf(debug
, "Initiating neighbourlist (ielec=%d, ivdw=%d, type=%d) for %s interactions,\nwith %d SR atoms.\n",
192 nl
->ielec
, nl
->ivdw
, nl
->type
, gmx_nblist_geometry_names
[nl
->igeometry
], maxsr
);
197 void init_neighbor_list(FILE *log
, t_forcerec
*fr
, int homenr
)
199 int maxsr
, maxsr_wat
;
200 int ielec
, ivdw
, ielecmod
, ivdwmod
, type
;
201 int igeometry_def
, igeometry_w
, igeometry_ww
;
203 gmx_bool bElecAndVdwSwitchDiffers
;
206 /* maxsr = homenr-fr->nWatMol*3; */
211 gmx_fatal(FARGS
, "%s, %d: Negative number of short range atoms.\n"
212 "Call your GROMACS dealer for assistance.", __FILE__
, __LINE__
);
214 /* This is just for initial allocation, so we do not reallocate
215 * all the nlist arrays many times in a row.
216 * The numbers seem very accurate, but they are uncritical.
218 maxsr_wat
= std::min(fr
->nWatMol
, (homenr
+2)/3);
220 /* Determine the values for ielec/ivdw. */
221 ielec
= fr
->nbkernel_elec_interaction
;
222 ivdw
= fr
->nbkernel_vdw_interaction
;
223 ielecmod
= fr
->nbkernel_elec_modifier
;
224 ivdwmod
= fr
->nbkernel_vdw_modifier
;
225 type
= GMX_NBLIST_INTERACTION_STANDARD
;
226 bElecAndVdwSwitchDiffers
= ( (fr
->rcoulomb_switch
!= fr
->rvdw_switch
) || (fr
->rcoulomb
!= fr
->rvdw
));
228 fr
->ns
->bCGlist
= (getenv("GMX_NBLISTCG") != nullptr);
229 if (!fr
->ns
->bCGlist
)
231 igeometry_def
= GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE
;
235 igeometry_def
= GMX_NBLIST_GEOMETRY_CG_CG
;
238 fprintf(log
, "\nUsing charge-group - charge-group neighbor lists and kernels\n\n");
242 if (fr
->solvent_opt
== esolTIP4P
)
244 igeometry_w
= GMX_NBLIST_GEOMETRY_WATER4_PARTICLE
;
245 igeometry_ww
= GMX_NBLIST_GEOMETRY_WATER4_WATER4
;
249 igeometry_w
= GMX_NBLIST_GEOMETRY_WATER3_PARTICLE
;
250 igeometry_ww
= GMX_NBLIST_GEOMETRY_WATER3_WATER3
;
253 for (i
= 0; i
< fr
->nnblists
; i
++)
255 nbl
= &(fr
->nblists
[i
]);
257 init_nblist(log
, &nbl
->nlist_sr
[eNL_VDWQQ
],
258 maxsr
, ivdw
, ivdwmod
, ielec
, ielecmod
, igeometry_def
, type
, bElecAndVdwSwitchDiffers
);
259 init_nblist(log
, &nbl
->nlist_sr
[eNL_VDW
],
260 maxsr
, ivdw
, ivdwmod
, GMX_NBKERNEL_ELEC_NONE
, eintmodNONE
, igeometry_def
, type
, bElecAndVdwSwitchDiffers
);
261 init_nblist(log
, &nbl
->nlist_sr
[eNL_QQ
],
262 maxsr
, GMX_NBKERNEL_VDW_NONE
, eintmodNONE
, ielec
, ielecmod
, igeometry_def
, type
, bElecAndVdwSwitchDiffers
);
263 init_nblist(log
, &nbl
->nlist_sr
[eNL_VDWQQ_WATER
],
264 maxsr_wat
, ivdw
, ivdwmod
, ielec
, ielecmod
, igeometry_w
, type
, bElecAndVdwSwitchDiffers
);
265 init_nblist(log
, &nbl
->nlist_sr
[eNL_QQ_WATER
],
266 maxsr_wat
, GMX_NBKERNEL_VDW_NONE
, eintmodNONE
, ielec
, ielecmod
, igeometry_w
, type
, bElecAndVdwSwitchDiffers
);
267 init_nblist(log
, &nbl
->nlist_sr
[eNL_VDWQQ_WATERWATER
],
268 maxsr_wat
, ivdw
, ivdwmod
, ielec
, ielecmod
, igeometry_ww
, type
, bElecAndVdwSwitchDiffers
);
269 init_nblist(log
, &nbl
->nlist_sr
[eNL_QQ_WATERWATER
],
270 maxsr_wat
, GMX_NBKERNEL_VDW_NONE
, eintmodNONE
, ielec
, ielecmod
, igeometry_ww
, type
, bElecAndVdwSwitchDiffers
);
272 /* Did we get the solvent loops so we can use optimized water kernels? */
273 if (nbl
->nlist_sr
[eNL_VDWQQ_WATER
].kernelptr_vf
== nullptr
274 || nbl
->nlist_sr
[eNL_QQ_WATER
].kernelptr_vf
== nullptr
275 || nbl
->nlist_sr
[eNL_VDWQQ_WATERWATER
].kernelptr_vf
== nullptr
276 || nbl
->nlist_sr
[eNL_QQ_WATERWATER
].kernelptr_vf
== nullptr)
278 fr
->solvent_opt
= esolNO
;
281 fprintf(log
, "Note: The available nonbonded kernels do not support water optimization - disabling.\n");
285 if (fr
->efep
!= efepNO
)
287 init_nblist(log
, &nbl
->nlist_sr
[eNL_VDWQQ_FREE
],
288 maxsr
, ivdw
, ivdwmod
, ielec
, ielecmod
, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE
, GMX_NBLIST_INTERACTION_FREE_ENERGY
, bElecAndVdwSwitchDiffers
);
289 init_nblist(log
, &nbl
->nlist_sr
[eNL_VDW_FREE
],
290 maxsr
, ivdw
, ivdwmod
, GMX_NBKERNEL_ELEC_NONE
, eintmodNONE
, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE
, GMX_NBLIST_INTERACTION_FREE_ENERGY
, bElecAndVdwSwitchDiffers
);
291 init_nblist(log
, &nbl
->nlist_sr
[eNL_QQ_FREE
],
292 maxsr
, GMX_NBKERNEL_VDW_NONE
, eintmodNONE
, ielec
, ielecmod
, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE
, GMX_NBLIST_INTERACTION_FREE_ENERGY
, bElecAndVdwSwitchDiffers
);
296 if (fr
->bQMMM
&& fr
->qr
->QMMMscheme
!= eQMMMschemeoniom
)
298 if (nullptr == fr
->QMMMlist
)
300 snew(fr
->QMMMlist
, 1);
302 init_nblist(log
, fr
->QMMMlist
,
303 maxsr
, 0, 0, ielec
, ielecmod
, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE
, GMX_NBLIST_INTERACTION_STANDARD
, bElecAndVdwSwitchDiffers
);
311 fr
->ns
->nblist_initialized
= TRUE
;
314 static void reset_nblist(t_nblist
*nl
)
324 static void reset_neighbor_lists(t_forcerec
*fr
)
330 /* only reset the short-range nblist */
331 reset_nblist(fr
->QMMMlist
);
334 for (n
= 0; n
< fr
->nnblists
; n
++)
336 for (i
= 0; i
< eNL_NR
; i
++)
338 reset_nblist( &(fr
->nblists
[n
].nlist_sr
[i
]) );
346 static gmx_inline
void new_i_nblist(t_nblist
*nlist
, int i_atom
, int shift
, int gid
)
348 int nri
= nlist
->nri
;
350 /* Check whether we have to increase the i counter */
352 (nlist
->iinr
[nri
] != i_atom
) ||
353 (nlist
->shift
[nri
] != shift
) ||
354 (nlist
->gid
[nri
] != gid
))
356 /* This is something else. Now see if any entries have
357 * been added in the list of the previous atom.
360 ((nlist
->jindex
[nri
+1] > nlist
->jindex
[nri
]) &&
361 (nlist
->gid
[nri
] != -1)))
363 /* If so increase the counter */
366 if (nlist
->nri
>= nlist
->maxnri
)
368 nlist
->maxnri
+= over_alloc_large(nlist
->nri
);
369 reallocate_nblist(nlist
);
372 /* Set the number of neighbours and the atom number */
373 nlist
->jindex
[nri
+1] = nlist
->jindex
[nri
];
374 nlist
->iinr
[nri
] = i_atom
;
375 nlist
->gid
[nri
] = gid
;
376 nlist
->shift
[nri
] = shift
;
380 /* Adding to previous list. First remove possible previous padding */
381 if (nlist
->simd_padding_width
> 1)
383 while (nlist
->nrj
> 0 && nlist
->jjnr
[nlist
->nrj
-1] < 0)
391 static gmx_inline
void close_i_nblist(t_nblist
*nlist
)
393 int nri
= nlist
->nri
;
398 /* Add elements up to padding. Since we allocate memory in units
399 * of the simd_padding width, we do not have to check for possible
400 * list reallocation here.
402 while ((nlist
->nrj
% nlist
->simd_padding_width
) != 0)
404 /* Use -4 here, so we can write forces for 4 atoms before real data */
405 nlist
->jjnr
[nlist
->nrj
++] = -4;
407 nlist
->jindex
[nri
+1] = nlist
->nrj
;
409 len
= nlist
->nrj
- nlist
->jindex
[nri
];
410 /* If there are no j-particles we have to reduce the
411 * number of i-particles again, to prevent errors in the
414 if ((len
== 0) && (nlist
->nri
> 0))
421 static gmx_inline
void close_nblist(t_nblist
*nlist
)
423 /* Only close this nblist when it has been initialized.
424 * Avoid the creation of i-lists with no j-particles.
428 /* Some assembly kernels do not support empty lists,
429 * make sure here that we don't generate any empty lists.
430 * With the current ns code this branch is taken in two cases:
431 * No i-particles at all: nri=-1 here
432 * There are i-particles, but no j-particles; nri=0 here
438 /* Close list number nri by incrementing the count */
443 static gmx_inline
void close_neighbor_lists(t_forcerec
*fr
, gmx_bool bMakeQMMMnblist
)
449 close_nblist(fr
->QMMMlist
);
452 for (n
= 0; n
< fr
->nnblists
; n
++)
454 for (i
= 0; (i
< eNL_NR
); i
++)
456 close_nblist(&(fr
->nblists
[n
].nlist_sr
[i
]));
462 static gmx_inline
void add_j_to_nblist(t_nblist
*nlist
, int j_atom
)
464 int nrj
= nlist
->nrj
;
466 if (nlist
->nrj
>= nlist
->maxnrj
)
468 nlist
->maxnrj
= round_up_to_simd_width(over_alloc_small(nlist
->nrj
+ 1), nlist
->simd_padding_width
);
472 fprintf(debug
, "Increasing SR nblist (ielec=%d,ivdw=%d,type=%d,igeometry=%d) j size to %d\n",
473 nlist
->ielec
, nlist
->ivdw
, nlist
->type
, nlist
->igeometry
, nlist
->maxnrj
);
476 srenew(nlist
->jjnr
, nlist
->maxnrj
);
479 nlist
->jjnr
[nrj
] = j_atom
;
483 static gmx_inline
void add_j_to_nblist_cg(t_nblist
*nlist
,
484 int j_start
, int j_end
,
485 t_excl
*bexcl
, gmx_bool i_is_j
)
487 int nrj
= nlist
->nrj
;
490 if (nlist
->nrj
>= nlist
->maxnrj
)
492 nlist
->maxnrj
= over_alloc_small(nlist
->nrj
+ 1);
495 fprintf(debug
, "Increasing SR nblist (ielec=%d,ivdw=%d,type=%d,igeometry=%d) j size to %d\n",
496 nlist
->ielec
, nlist
->ivdw
, nlist
->type
, nlist
->igeometry
, nlist
->maxnrj
);
499 srenew(nlist
->jjnr
, nlist
->maxnrj
);
500 srenew(nlist
->jjnr_end
, nlist
->maxnrj
);
501 srenew(nlist
->excl
, nlist
->maxnrj
*MAX_CGCGSIZE
);
504 nlist
->jjnr
[nrj
] = j_start
;
505 nlist
->jjnr_end
[nrj
] = j_end
;
507 if (j_end
- j_start
> MAX_CGCGSIZE
)
509 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
);
512 /* Set the exclusions */
513 for (j
= j_start
; j
< j_end
; j
++)
515 nlist
->excl
[nrj
*MAX_CGCGSIZE
+ j
- j_start
] = bexcl
[j
];
519 /* Avoid double counting of intra-cg interactions */
520 for (j
= 1; j
< j_end
-j_start
; j
++)
522 nlist
->excl
[nrj
*MAX_CGCGSIZE
+ j
] |= (1<<j
) - 1;
530 put_in_list_t (gmx_bool bHaveVdW
[],
546 put_in_list_at(gmx_bool bHaveVdW
[],
561 /* The a[] index has been removed,
562 * to put it back in i_atom should be a[i0] and jj should be a[jj].
567 t_nblist
* vdwc_free
= nullptr;
568 t_nblist
* vdw_free
= nullptr;
569 t_nblist
* coul_free
= nullptr;
570 t_nblist
* vdwc_ww
= nullptr;
571 t_nblist
* coul_ww
= nullptr;
573 int i
, j
, jcg
, igid
, gid
, nbl_ind
;
574 int jj
, jj0
, jj1
, i_atom
;
579 real
*charge
, *chargeB
;
581 gmx_bool bFreeEnergy
, bFree
, bFreeJ
, bNotEx
, *bPert
;
582 gmx_bool bDoVdW_i
, bDoCoul_i
, bDoCoul_i_sol
;
586 /* Copy some pointers */
588 charge
= md
->chargeA
;
589 chargeB
= md
->chargeB
;
592 bPert
= md
->bPerturbed
;
596 nicg
= index
[icg
+1]-i0
;
598 /* Get the i charge group info */
599 igid
= GET_CGINFO_GID(cginfo
[icg
]);
601 iwater
= (solvent_opt
!= esolNO
) ? GET_CGINFO_SOLOPT(cginfo
[icg
]) : esolNO
;
606 /* Check if any of the particles involved are perturbed.
607 * If not we can do the cheaper normal put_in_list
608 * and use more solvent optimization.
610 for (i
= 0; i
< nicg
; i
++)
612 bFreeEnergy
|= bPert
[i0
+i
];
614 /* Loop over the j charge groups */
615 for (j
= 0; (j
< nj
&& !bFreeEnergy
); j
++)
620 /* Finally loop over the atoms in the j-charge group */
621 for (jj
= jj0
; jj
< jj1
; jj
++)
623 bFreeEnergy
|= bPert
[jj
];
628 /* Unpack pointers to neighbourlist structs */
629 if (fr
->nnblists
== 1)
635 nbl_ind
= fr
->gid2nblists
[GID(igid
, jgid
, ngid
)];
637 nlist
= fr
->nblists
[nbl_ind
].nlist_sr
;
639 if (iwater
!= esolNO
)
641 vdwc
= &nlist
[eNL_VDWQQ_WATER
];
642 vdw
= &nlist
[eNL_VDW
];
643 coul
= &nlist
[eNL_QQ_WATER
];
644 vdwc_ww
= &nlist
[eNL_VDWQQ_WATERWATER
];
645 coul_ww
= &nlist
[eNL_QQ_WATERWATER
];
649 vdwc
= &nlist
[eNL_VDWQQ
];
650 vdw
= &nlist
[eNL_VDW
];
651 coul
= &nlist
[eNL_QQ
];
656 if (iwater
!= esolNO
)
658 /* Loop over the atoms in the i charge group */
660 gid
= GID(igid
, jgid
, ngid
);
661 /* Create new i_atom for each energy group */
662 if (bDoCoul
&& bDoVdW
)
664 new_i_nblist(vdwc
, i_atom
, shift
, gid
);
665 new_i_nblist(vdwc_ww
, i_atom
, shift
, gid
);
669 new_i_nblist(vdw
, i_atom
, shift
, gid
);
673 new_i_nblist(coul
, i_atom
, shift
, gid
);
674 new_i_nblist(coul_ww
, i_atom
, shift
, gid
);
676 /* Loop over the j charge groups */
677 for (j
= 0; (j
< nj
); j
++)
687 jwater
= GET_CGINFO_SOLOPT(cginfo
[jcg
]);
689 if (iwater
== esolSPC
&& jwater
== esolSPC
)
691 /* Interaction between two SPC molecules */
694 /* VdW only - only first atoms in each water interact */
695 add_j_to_nblist(vdw
, jj0
);
699 /* One entry for the entire water-water interaction */
702 add_j_to_nblist(coul_ww
, jj0
);
706 add_j_to_nblist(vdwc_ww
, jj0
);
710 else if (iwater
== esolTIP4P
&& jwater
== esolTIP4P
)
712 /* Interaction between two TIP4p molecules */
715 /* VdW only - only first atoms in each water interact */
716 add_j_to_nblist(vdw
, jj0
);
720 /* One entry for the entire water-water interaction */
723 add_j_to_nblist(coul_ww
, jj0
);
727 add_j_to_nblist(vdwc_ww
, jj0
);
733 /* j charge group is not water, but i is.
734 * Add entries to the water-other_atom lists; the geometry of the water
735 * molecule doesn't matter - that is taken care of in the nonbonded kernel,
736 * so we don't care if it is SPC or TIP4P...
743 for (jj
= jj0
; (jj
< jj1
); jj
++)
747 add_j_to_nblist(coul
, jj
);
753 for (jj
= jj0
; (jj
< jj1
); jj
++)
755 if (bHaveVdW
[type
[jj
]])
757 add_j_to_nblist(vdw
, jj
);
763 /* _charge_ _groups_ interact with both coulomb and LJ */
764 /* Check which atoms we should add to the lists! */
765 for (jj
= jj0
; (jj
< jj1
); jj
++)
767 if (bHaveVdW
[type
[jj
]])
771 add_j_to_nblist(vdwc
, jj
);
775 add_j_to_nblist(vdw
, jj
);
778 else if (charge
[jj
] != 0)
780 add_j_to_nblist(coul
, jj
);
787 close_i_nblist(coul
);
788 close_i_nblist(vdwc
);
789 close_i_nblist(coul_ww
);
790 close_i_nblist(vdwc_ww
);
794 /* no solvent as i charge group */
795 /* Loop over the atoms in the i charge group */
796 for (i
= 0; i
< nicg
; i
++)
799 gid
= GID(igid
, jgid
, ngid
);
802 /* Create new i_atom for each energy group */
803 if (bDoVdW
&& bDoCoul
)
805 new_i_nblist(vdwc
, i_atom
, shift
, gid
);
809 new_i_nblist(vdw
, i_atom
, shift
, gid
);
813 new_i_nblist(coul
, i_atom
, shift
, gid
);
815 bDoVdW_i
= (bDoVdW
&& bHaveVdW
[type
[i_atom
]]);
816 bDoCoul_i
= (bDoCoul
&& qi
!= 0);
818 if (bDoVdW_i
|| bDoCoul_i
)
820 /* Loop over the j charge groups */
821 for (j
= 0; (j
< nj
); j
++)
825 /* Check for large charge groups */
836 /* Finally loop over the atoms in the j-charge group */
837 for (jj
= jj0
; jj
< jj1
; jj
++)
839 bNotEx
= NOTEXCL(bExcl
, i
, jj
);
847 add_j_to_nblist(coul
, jj
);
852 if (bHaveVdW
[type
[jj
]])
854 add_j_to_nblist(vdw
, jj
);
859 if (bHaveVdW
[type
[jj
]])
863 add_j_to_nblist(vdwc
, jj
);
867 add_j_to_nblist(vdw
, jj
);
870 else if (charge
[jj
] != 0)
872 add_j_to_nblist(coul
, jj
);
880 close_i_nblist(coul
);
881 close_i_nblist(vdwc
);
887 /* we are doing free energy */
888 vdwc_free
= &nlist
[eNL_VDWQQ_FREE
];
889 vdw_free
= &nlist
[eNL_VDW_FREE
];
890 coul_free
= &nlist
[eNL_QQ_FREE
];
891 /* Loop over the atoms in the i charge group */
892 for (i
= 0; i
< nicg
; i
++)
895 gid
= GID(igid
, jgid
, ngid
);
897 qiB
= chargeB
[i_atom
];
899 /* Create new i_atom for each energy group */
900 if (bDoVdW
&& bDoCoul
)
902 new_i_nblist(vdwc
, i_atom
, shift
, gid
);
906 new_i_nblist(vdw
, i_atom
, shift
, gid
);
910 new_i_nblist(coul
, i_atom
, shift
, gid
);
913 new_i_nblist(vdw_free
, i_atom
, shift
, gid
);
914 new_i_nblist(coul_free
, i_atom
, shift
, gid
);
915 new_i_nblist(vdwc_free
, i_atom
, shift
, gid
);
917 bDoVdW_i
= (bDoVdW
&&
918 (bHaveVdW
[type
[i_atom
]] || bHaveVdW
[typeB
[i_atom
]]));
919 bDoCoul_i
= (bDoCoul
&& (qi
!= 0 || qiB
!= 0));
920 /* For TIP4P the first atom does not have a charge,
921 * but the last three do. So we should still put an atom
922 * without LJ but with charge in the water-atom neighborlist
923 * for a TIP4p i charge group.
924 * For SPC type water the first atom has LJ and charge,
925 * so there is no such problem.
927 if (iwater
== esolNO
)
929 bDoCoul_i_sol
= bDoCoul_i
;
933 bDoCoul_i_sol
= bDoCoul
;
936 if (bDoVdW_i
|| bDoCoul_i_sol
)
938 /* Loop over the j charge groups */
939 for (j
= 0; (j
< nj
); j
++)
943 /* Check for large charge groups */
954 /* Finally loop over the atoms in the j-charge group */
955 bFree
= bPert
[i_atom
];
956 for (jj
= jj0
; (jj
< jj1
); jj
++)
958 bFreeJ
= bFree
|| bPert
[jj
];
959 /* Complicated if, because the water H's should also
960 * see perturbed j-particles
962 if (iwater
== esolNO
|| i
== 0 || bFreeJ
)
964 bNotEx
= NOTEXCL(bExcl
, i
, jj
);
972 if (charge
[jj
] != 0 || chargeB
[jj
] != 0)
974 add_j_to_nblist(coul_free
, jj
);
979 if (bHaveVdW
[type
[jj
]] || bHaveVdW
[typeB
[jj
]])
981 add_j_to_nblist(vdw_free
, jj
);
986 if (bHaveVdW
[type
[jj
]] || bHaveVdW
[typeB
[jj
]])
988 if (charge
[jj
] != 0 || chargeB
[jj
] != 0)
990 add_j_to_nblist(vdwc_free
, jj
);
994 add_j_to_nblist(vdw_free
, jj
);
997 else if (charge
[jj
] != 0 || chargeB
[jj
] != 0)
999 add_j_to_nblist(coul_free
, jj
);
1005 /* This is done whether or not bWater is set */
1006 if (charge
[jj
] != 0)
1008 add_j_to_nblist(coul
, jj
);
1011 else if (!bDoCoul_i_sol
)
1013 if (bHaveVdW
[type
[jj
]])
1015 add_j_to_nblist(vdw
, jj
);
1020 if (bHaveVdW
[type
[jj
]])
1022 if (charge
[jj
] != 0)
1024 add_j_to_nblist(vdwc
, jj
);
1028 add_j_to_nblist(vdw
, jj
);
1031 else if (charge
[jj
] != 0)
1033 add_j_to_nblist(coul
, jj
);
1041 close_i_nblist(vdw
);
1042 close_i_nblist(coul
);
1043 close_i_nblist(vdwc
);
1044 close_i_nblist(vdw_free
);
1045 close_i_nblist(coul_free
);
1046 close_i_nblist(vdwc_free
);
1052 put_in_list_qmmm(gmx_bool gmx_unused bHaveVdW
[],
1054 t_mdatoms gmx_unused
* md
,
1063 gmx_bool gmx_unused bDoVdW
,
1064 gmx_bool gmx_unused bDoCoul
,
1065 int gmx_unused solvent_opt
)
1068 int i
, j
, jcg
, igid
, gid
;
1069 int jj
, jj0
, jj1
, i_atom
;
1073 /* Get atom range */
1075 nicg
= index
[icg
+1]-i0
;
1077 /* Get the i charge group info */
1078 igid
= GET_CGINFO_GID(fr
->cginfo
[icg
]);
1080 coul
= fr
->QMMMlist
;
1082 /* Loop over atoms in the ith charge group */
1083 for (i
= 0; i
< nicg
; i
++)
1086 gid
= GID(igid
, jgid
, ngid
);
1087 /* Create new i_atom for each energy group */
1088 new_i_nblist(coul
, i_atom
, shift
, gid
);
1090 /* Loop over the j charge groups */
1091 for (j
= 0; j
< nj
; j
++)
1095 /* Charge groups cannot have QM and MM atoms simultaneously */
1100 /* Finally loop over the atoms in the j-charge group */
1101 for (jj
= jj0
; jj
< jj1
; jj
++)
1103 bNotEx
= NOTEXCL(bExcl
, i
, jj
);
1106 add_j_to_nblist(coul
, jj
);
1111 close_i_nblist(coul
);
1116 put_in_list_cg(gmx_bool gmx_unused bHaveVdW
[],
1118 t_mdatoms gmx_unused
* md
,
1127 gmx_bool gmx_unused bDoVdW
,
1128 gmx_bool gmx_unused bDoCoul
,
1129 int gmx_unused solvent_opt
)
1132 int igid
, gid
, nbl_ind
;
1136 cginfo
= fr
->cginfo
[icg
];
1138 igid
= GET_CGINFO_GID(cginfo
);
1139 gid
= GID(igid
, jgid
, ngid
);
1141 /* Unpack pointers to neighbourlist structs */
1142 if (fr
->nnblists
== 1)
1148 nbl_ind
= fr
->gid2nblists
[gid
];
1150 vdwc
= &fr
->nblists
[nbl_ind
].nlist_sr
[eNL_VDWQQ
];
1152 /* Make a new neighbor list for charge group icg.
1153 * Currently simply one neighbor list is made with LJ and Coulomb.
1154 * If required, zero interactions could be removed here
1155 * or in the force loop.
1157 new_i_nblist(vdwc
, index
[icg
], shift
, gid
);
1158 vdwc
->iinr_end
[vdwc
->nri
] = index
[icg
+1];
1160 for (j
= 0; (j
< nj
); j
++)
1163 /* Skip the icg-icg pairs if all self interactions are excluded */
1164 if (!(jcg
== icg
&& GET_CGINFO_EXCL_INTRA(cginfo
)))
1166 /* Here we add the j charge group jcg to the list,
1167 * exclusions are also added to the list.
1169 add_j_to_nblist_cg(vdwc
, index
[jcg
], index
[jcg
+1], bExcl
, icg
== jcg
);
1173 close_i_nblist(vdwc
);
1176 static void setexcl(int start
, int end
, t_blocka
*excl
, gmx_bool b
,
1183 for (i
= start
; i
< end
; i
++)
1185 for (k
= excl
->index
[i
]; k
< excl
->index
[i
+1]; k
++)
1187 SETEXCL(bexcl
, i
-start
, excl
->a
[k
]);
1193 for (i
= start
; i
< end
; i
++)
1195 for (k
= excl
->index
[i
]; k
< excl
->index
[i
+1]; k
++)
1197 RMEXCL(bexcl
, i
-start
, excl
->a
[k
]);
1203 int calc_naaj(int icg
, int cgtot
)
1207 if ((cgtot
% 2) == 1)
1209 /* Odd number of charge groups, easy */
1210 naaj
= 1 + (cgtot
/2);
1212 else if ((cgtot
% 4) == 0)
1214 /* Multiple of four is hard */
1251 fprintf(log
, "naaj=%d\n", naaj
);
1257 /************************************************
1259 * S I M P L E C O R E S T U F F
1261 ************************************************/
1263 static real
calc_image_tric(rvec xi
, rvec xj
, matrix box
,
1264 rvec b_inv
, int *shift
)
1266 /* This code assumes that the cut-off is smaller than
1267 * a half times the smallest diagonal element of the box.
1269 const real h25
= 2.5;
1274 /* Compute diff vector */
1275 dz
= xj
[ZZ
] - xi
[ZZ
];
1276 dy
= xj
[YY
] - xi
[YY
];
1277 dx
= xj
[XX
] - xi
[XX
];
1279 /* Perform NINT operation, using trunc operation, therefore
1280 * we first add 2.5 then subtract 2 again
1282 tz
= static_cast<int>(dz
*b_inv
[ZZ
] + h25
);
1284 dz
-= tz
*box
[ZZ
][ZZ
];
1285 dy
-= tz
*box
[ZZ
][YY
];
1286 dx
-= tz
*box
[ZZ
][XX
];
1288 ty
= static_cast<int>(dy
*b_inv
[YY
] + h25
);
1290 dy
-= ty
*box
[YY
][YY
];
1291 dx
-= ty
*box
[YY
][XX
];
1293 tx
= static_cast<int>(dx
*b_inv
[XX
]+h25
);
1295 dx
-= tx
*box
[XX
][XX
];
1297 /* Distance squared */
1298 r2
= (dx
*dx
) + (dy
*dy
) + (dz
*dz
);
1300 *shift
= XYZ2IS(tx
, ty
, tz
);
1305 static real
calc_image_rect(rvec xi
, rvec xj
, rvec box_size
,
1306 rvec b_inv
, int *shift
)
1308 const real h15
= 1.5;
1314 /* Compute diff vector */
1315 dx
= xj
[XX
] - xi
[XX
];
1316 dy
= xj
[YY
] - xi
[YY
];
1317 dz
= xj
[ZZ
] - xi
[ZZ
];
1319 /* Perform NINT operation, using trunc operation, therefore
1320 * we first add 1.5 then subtract 1 again
1322 tx
= static_cast<int>(dx
*b_inv
[XX
] + h15
);
1323 ty
= static_cast<int>(dy
*b_inv
[YY
] + h15
);
1324 tz
= static_cast<int>(dz
*b_inv
[ZZ
] + h15
);
1329 /* Correct diff vector for translation */
1330 ddx
= tx
*box_size
[XX
] - dx
;
1331 ddy
= ty
*box_size
[YY
] - dy
;
1332 ddz
= tz
*box_size
[ZZ
] - dz
;
1334 /* Distance squared */
1335 r2
= (ddx
*ddx
) + (ddy
*ddy
) + (ddz
*ddz
);
1337 *shift
= XYZ2IS(tx
, ty
, tz
);
1342 static void add_simple(t_ns_buf
* nsbuf
, int nrj
, int cg_j
,
1343 gmx_bool bHaveVdW
[], int ngid
, t_mdatoms
*md
,
1344 int icg
, int jgid
, t_block
*cgs
, t_excl bexcl
[],
1345 int shift
, t_forcerec
*fr
, put_in_list_t
*put_in_list
)
1347 if (nsbuf
->nj
+ nrj
> MAX_CG
)
1349 put_in_list(bHaveVdW
, ngid
, md
, icg
, jgid
, nsbuf
->ncg
, nsbuf
->jcg
,
1350 cgs
->index
, bexcl
, shift
, fr
, TRUE
, TRUE
, fr
->solvent_opt
);
1351 /* Reset buffer contents */
1352 nsbuf
->ncg
= nsbuf
->nj
= 0;
1354 nsbuf
->jcg
[nsbuf
->ncg
++] = cg_j
;
1358 static void ns_inner_tric(rvec x
[], int icg
, int *i_egp_flags
,
1359 int njcg
, int jcg
[],
1360 matrix box
, rvec b_inv
, real rcut2
,
1361 t_block
*cgs
, t_ns_buf
**ns_buf
,
1362 gmx_bool bHaveVdW
[], int ngid
, t_mdatoms
*md
,
1363 t_excl bexcl
[], t_forcerec
*fr
,
1364 put_in_list_t
*put_in_list
)
1368 int *cginfo
= fr
->cginfo
;
1371 cgindex
= cgs
->index
;
1373 for (j
= 0; (j
< njcg
); j
++)
1376 nrj
= cgindex
[cg_j
+1]-cgindex
[cg_j
];
1377 if (calc_image_tric(x
[icg
], x
[cg_j
], box
, b_inv
, &shift
) < rcut2
)
1379 jgid
= GET_CGINFO_GID(cginfo
[cg_j
]);
1380 if (!(i_egp_flags
[jgid
] & EGP_EXCL
))
1382 add_simple(&ns_buf
[jgid
][shift
], nrj
, cg_j
,
1383 bHaveVdW
, ngid
, md
, icg
, jgid
, cgs
, bexcl
, shift
, fr
,
1390 static void ns_inner_rect(rvec x
[], int icg
, int *i_egp_flags
,
1391 int njcg
, int jcg
[],
1392 gmx_bool bBox
, rvec box_size
, rvec b_inv
, real rcut2
,
1393 t_block
*cgs
, t_ns_buf
**ns_buf
,
1394 gmx_bool bHaveVdW
[], int ngid
, t_mdatoms
*md
,
1395 t_excl bexcl
[], t_forcerec
*fr
,
1396 put_in_list_t
*put_in_list
)
1400 int *cginfo
= fr
->cginfo
;
1403 cgindex
= cgs
->index
;
1407 for (j
= 0; (j
< njcg
); j
++)
1410 nrj
= cgindex
[cg_j
+1]-cgindex
[cg_j
];
1411 if (calc_image_rect(x
[icg
], x
[cg_j
], box_size
, b_inv
, &shift
) < rcut2
)
1413 jgid
= GET_CGINFO_GID(cginfo
[cg_j
]);
1414 if (!(i_egp_flags
[jgid
] & EGP_EXCL
))
1416 add_simple(&ns_buf
[jgid
][shift
], nrj
, cg_j
,
1417 bHaveVdW
, ngid
, md
, icg
, jgid
, cgs
, bexcl
, shift
, fr
,
1425 for (j
= 0; (j
< njcg
); j
++)
1428 nrj
= cgindex
[cg_j
+1]-cgindex
[cg_j
];
1429 if ((rcut2
== 0) || (distance2(x
[icg
], x
[cg_j
]) < rcut2
))
1431 jgid
= GET_CGINFO_GID(cginfo
[cg_j
]);
1432 if (!(i_egp_flags
[jgid
] & EGP_EXCL
))
1434 add_simple(&ns_buf
[jgid
][CENTRAL
], nrj
, cg_j
,
1435 bHaveVdW
, ngid
, md
, icg
, jgid
, cgs
, bexcl
, CENTRAL
, fr
,
1443 /* ns_simple_core needs to be adapted for QMMM still 2005 */
1445 static int ns_simple_core(t_forcerec
*fr
,
1446 gmx_localtop_t
*top
,
1448 matrix box
, rvec box_size
,
1449 t_excl bexcl
[], int *aaj
,
1450 int ngid
, t_ns_buf
**ns_buf
,
1451 put_in_list_t
*put_in_list
, gmx_bool bHaveVdW
[])
1455 int nsearch
, icg
, igid
, nn
;
1459 t_block
*cgs
= &(top
->cgs
);
1460 t_blocka
*excl
= &(top
->excls
);
1463 gmx_bool bBox
, bTriclinic
;
1466 rlist2
= gmx::square(fr
->rlist
);
1468 bBox
= (fr
->ePBC
!= epbcNONE
);
1471 for (m
= 0; (m
< DIM
); m
++)
1473 if (gmx_numzero(box_size
[m
]))
1475 gmx_fatal(FARGS
, "Dividing by zero box size!");
1477 b_inv
[m
] = 1.0/box_size
[m
];
1479 bTriclinic
= TRICLINIC(box
);
1486 cginfo
= fr
->cginfo
;
1489 for (icg
= fr
->cg0
; (icg
< fr
->hcg
); icg
++)
1492 i0 = cgs->index[icg];
1493 nri = cgs->index[icg+1]-i0;
1494 i_atoms = &(cgs->a[i0]);
1495 i_eg_excl = fr->eg_excl + ngid*md->cENER[*i_atoms];
1496 setexcl(nri,i_atoms,excl,TRUE,bexcl);
1498 igid
= GET_CGINFO_GID(cginfo
[icg
]);
1499 i_egp_flags
= fr
->egp_flags
+ ngid
*igid
;
1500 setexcl(cgs
->index
[icg
], cgs
->index
[icg
+1], excl
, TRUE
, bexcl
);
1502 naaj
= calc_naaj(icg
, cgs
->nr
);
1505 ns_inner_tric(fr
->cg_cm
, icg
, i_egp_flags
, naaj
, &(aaj
[icg
]),
1506 box
, b_inv
, rlist2
, cgs
, ns_buf
,
1507 bHaveVdW
, ngid
, md
, bexcl
, fr
, put_in_list
);
1511 ns_inner_rect(fr
->cg_cm
, icg
, i_egp_flags
, naaj
, &(aaj
[icg
]),
1512 bBox
, box_size
, b_inv
, rlist2
, cgs
, ns_buf
,
1513 bHaveVdW
, ngid
, md
, bexcl
, fr
, put_in_list
);
1517 for (nn
= 0; (nn
< ngid
); nn
++)
1519 for (k
= 0; (k
< SHIFTS
); k
++)
1521 nsbuf
= &(ns_buf
[nn
][k
]);
1524 put_in_list(bHaveVdW
, ngid
, md
, icg
, nn
, nsbuf
->ncg
, nsbuf
->jcg
,
1525 cgs
->index
, bexcl
, k
, fr
, TRUE
, TRUE
, fr
->solvent_opt
);
1526 nsbuf
->ncg
= nsbuf
->nj
= 0;
1530 /* setexcl(nri,i_atoms,excl,FALSE,bexcl); */
1531 setexcl(cgs
->index
[icg
], cgs
->index
[icg
+1], excl
, FALSE
, bexcl
);
1533 close_neighbor_lists(fr
, FALSE
);
1538 /************************************************
1540 * N S 5 G R I D S T U F F
1542 ************************************************/
1544 static gmx_inline
void get_dx_dd(int Nx
, real gridx
, real rc2
, int xgi
, real x
,
1545 int ncpddc
, int shift_min
, int shift_max
,
1546 int *g0
, int *g1
, real
*dcx2
)
1549 int g_min
, g_max
, shift_home
;
1582 g_min
= (shift_min
== shift_home
? 0 : ncpddc
);
1583 g_max
= (shift_max
== shift_home
? ncpddc
- 1 : Nx
- 1);
1590 else if (shift_max
< 0)
1605 /* Check one grid cell down */
1606 dcx
= ((*g0
- 1) + 1)*gridx
- x
;
1618 /* Check one grid cell up */
1619 dcx
= (*g1
+ 1)*gridx
- x
;
1631 #define calc_dx2(XI, YI, ZI, y) (gmx::square(XI-y[XX]) + gmx::square(YI-y[YY]) + gmx::square(ZI-y[ZZ]))
1632 #define calc_cyl_dx2(XI, YI, y) (gmx::square(XI-y[XX]) + gmx::square(YI-y[YY]))
1633 /****************************************************
1635 * F A S T N E I G H B O R S E A R C H I N G
1637 * Optimized neighboursearching routine using grid
1638 * at least 1x1x1, see GROMACS manual
1640 ****************************************************/
1643 static void get_cutoff2(t_forcerec
*fr
, real
*rs2
)
1645 *rs2
= gmx::square(fr
->rlist
);
1648 static void init_nsgrid_lists(t_forcerec
*fr
, int ngid
, gmx_ns_t
*ns
)
1653 get_cutoff2(fr
, &rs2
);
1655 /* Short range buffers */
1656 snew(ns
->nl_sr
, ngid
);
1658 snew(ns
->nsr
, ngid
);
1660 for (j
= 0; (j
< ngid
); j
++)
1662 snew(ns
->nl_sr
[j
], MAX_CG
);
1667 "ns5_core: rs2 = %g (nm^2)\n",
1672 static int nsgrid_core(t_commrec
*cr
, t_forcerec
*fr
,
1673 matrix box
, int ngid
,
1674 gmx_localtop_t
*top
,
1676 t_excl bexcl
[], gmx_bool
*bExcludeAlleg
,
1678 put_in_list_t
*put_in_list
,
1679 gmx_bool bHaveVdW
[],
1680 gmx_bool bMakeQMMMnblist
)
1686 t_block
*cgs
= &(top
->cgs
);
1687 int *cginfo
= fr
->cginfo
;
1688 /* int *i_atoms,*cgsindex=cgs->index; */
1690 int cell_x
, cell_y
, cell_z
;
1691 int d
, tx
, ty
, tz
, dx
, dy
, dz
, cj
;
1692 #ifdef ALLOW_OFFDIAG_LT_HALFDIAG
1693 int zsh_ty
, zsh_tx
, ysh_tx
;
1695 int dx0
, dx1
, dy0
, dy1
, dz0
, dz1
;
1696 int Nx
, Ny
, Nz
, shift
= -1, j
, nrj
, nns
, nn
= -1;
1697 real gridx
, gridy
, gridz
, grid_x
, grid_y
;
1698 real
*dcx2
, *dcy2
, *dcz2
;
1700 int cg0
, cg1
, icg
= -1, cgsnr
, i0
, igid
, naaj
, max_jcg
;
1701 int jcg0
, jcg1
, jjcg
, cgj0
, jgid
;
1702 int *grida
, *gridnra
, *gridind
;
1703 rvec
*cgcm
, grid_offset
;
1704 real r2
, rs2
, XI
, YI
, ZI
, tmp1
, tmp2
;
1706 gmx_bool bDomDec
, bTriclinicX
, bTriclinicY
;
1711 bDomDec
= DOMAINDECOMP(cr
);
1714 bTriclinicX
= ((YY
< grid
->npbcdim
&&
1715 (!bDomDec
|| dd
->nc
[YY
] == 1) && box
[YY
][XX
] != 0) ||
1716 (ZZ
< grid
->npbcdim
&&
1717 (!bDomDec
|| dd
->nc
[ZZ
] == 1) && box
[ZZ
][XX
] != 0));
1718 bTriclinicY
= (ZZ
< grid
->npbcdim
&&
1719 (!bDomDec
|| dd
->nc
[ZZ
] == 1) && box
[ZZ
][YY
] != 0);
1723 get_cutoff2(fr
, &rs2
);
1734 gridind
= grid
->index
;
1735 gridnra
= grid
->nra
;
1738 gridx
= grid
->cell_size
[XX
];
1739 gridy
= grid
->cell_size
[YY
];
1740 gridz
= grid
->cell_size
[ZZ
];
1743 copy_rvec(grid
->cell_offset
, grid_offset
);
1744 copy_ivec(grid
->ncpddc
, ncpddc
);
1749 #ifdef ALLOW_OFFDIAG_LT_HALFDIAG
1750 zsh_ty
= floor(-box
[ZZ
][YY
]/box
[YY
][YY
]+0.5);
1751 zsh_tx
= floor(-box
[ZZ
][XX
]/box
[XX
][XX
]+0.5);
1752 ysh_tx
= floor(-box
[YY
][XX
]/box
[XX
][XX
]+0.5);
1753 if (zsh_tx
!= 0 && ysh_tx
!= 0)
1755 /* This could happen due to rounding, when both ratios are 0.5 */
1762 /* We only want a list for the test particle */
1771 /* Set the shift range */
1772 for (d
= 0; d
< DIM
; d
++)
1776 /* Check if we need periodicity shifts.
1777 * Without PBC or with domain decomposition we don't need them.
1779 if (d
>= ePBC2npbcdim(fr
->ePBC
) || (bDomDec
&& dd
->nc
[d
] > 1))
1786 box
[XX
][XX
] - std::abs(box
[YY
][XX
]) - std::abs(box
[ZZ
][XX
]) < std::sqrt(rs2
))
1797 /* Loop over charge groups */
1798 for (icg
= cg0
; (icg
< cg1
); icg
++)
1800 igid
= GET_CGINFO_GID(cginfo
[icg
]);
1801 /* Skip this charge group if all energy groups are excluded! */
1802 if (bExcludeAlleg
[igid
])
1807 i0
= cgs
->index
[icg
];
1809 if (bMakeQMMMnblist
)
1811 /* Skip this charge group if it is not a QM atom while making a
1812 * QM/MM neighbourlist
1814 if (md
->bQM
[i0
] == FALSE
)
1816 continue; /* MM particle, go to next particle */
1819 /* Compute the number of charge groups that fall within the control
1822 naaj
= calc_naaj(icg
, cgsnr
);
1829 /* make a normal neighbourlist */
1833 /* Get the j charge-group and dd cell shift ranges */
1834 dd_get_ns_ranges(cr
->dd
, icg
, &jcg0
, &jcg1
, sh0
, sh1
);
1839 /* Compute the number of charge groups that fall within the control
1842 naaj
= calc_naaj(icg
, cgsnr
);
1848 /* The i-particle is awlways the test particle,
1849 * so we want all j-particles
1851 max_jcg
= cgsnr
- 1;
1855 max_jcg
= jcg1
- cgsnr
;
1860 i_egp_flags
= fr
->egp_flags
+ igid
*ngid
;
1862 /* Set the exclusions for the atoms in charge group icg using a bitmask */
1863 setexcl(i0
, cgs
->index
[icg
+1], &top
->excls
, TRUE
, bexcl
);
1865 ci2xyz(grid
, icg
, &cell_x
, &cell_y
, &cell_z
);
1867 /* Changed iicg to icg, DvdS 990115
1868 * (but see consistency check above, DvdS 990330)
1871 fprintf(log
, "icg=%5d, naaj=%5d, cell %d %d %d\n",
1872 icg
, naaj
, cell_x
, cell_y
, cell_z
);
1874 /* Loop over shift vectors in three dimensions */
1875 for (tz
= -shp
[ZZ
]; tz
<= shp
[ZZ
]; tz
++)
1877 ZI
= cgcm
[icg
][ZZ
]+tz
*box
[ZZ
][ZZ
];
1878 /* Calculate range of cells in Z direction that have the shift tz */
1879 zgi
= cell_z
+ tz
*Nz
;
1880 get_dx_dd(Nz
, gridz
, rs2
, zgi
, ZI
-grid_offset
[ZZ
],
1881 ncpddc
[ZZ
], sh0
[ZZ
], sh1
[ZZ
], &dz0
, &dz1
, dcz2
);
1886 for (ty
= -shp
[YY
]; ty
<= shp
[YY
]; ty
++)
1888 YI
= cgcm
[icg
][YY
]+ty
*box
[YY
][YY
]+tz
*box
[ZZ
][YY
];
1889 /* Calculate range of cells in Y direction that have the shift ty */
1892 ygi
= (int)(Ny
+ (YI
- grid_offset
[YY
])*grid_y
) - Ny
;
1896 ygi
= cell_y
+ ty
*Ny
;
1898 get_dx_dd(Ny
, gridy
, rs2
, ygi
, YI
-grid_offset
[YY
],
1899 ncpddc
[YY
], sh0
[YY
], sh1
[YY
], &dy0
, &dy1
, dcy2
);
1904 for (tx
= -shp
[XX
]; tx
<= shp
[XX
]; tx
++)
1906 XI
= cgcm
[icg
][XX
]+tx
*box
[XX
][XX
]+ty
*box
[YY
][XX
]+tz
*box
[ZZ
][XX
];
1907 /* Calculate range of cells in X direction that have the shift tx */
1910 xgi
= (int)(Nx
+ (XI
- grid_offset
[XX
])*grid_x
) - Nx
;
1914 xgi
= cell_x
+ tx
*Nx
;
1916 get_dx_dd(Nx
, gridx
, rs2
, xgi
, XI
-grid_offset
[XX
],
1917 ncpddc
[XX
], sh0
[XX
], sh1
[XX
], &dx0
, &dx1
, dcx2
);
1922 /* Get shift vector */
1923 shift
= XYZ2IS(tx
, ty
, tz
);
1925 range_check(shift
, 0, SHIFTS
);
1927 for (nn
= 0; (nn
< ngid
); nn
++)
1932 fprintf(log
, "shift: %2d, dx0,1: %2d,%2d, dy0,1: %2d,%2d, dz0,1: %2d,%2d\n",
1933 shift
, dx0
, dx1
, dy0
, dy1
, dz0
, dz1
);
1934 fprintf(log
, "cgcm: %8.3f %8.3f %8.3f\n", cgcm
[icg
][XX
],
1935 cgcm
[icg
][YY
], cgcm
[icg
][ZZ
]);
1936 fprintf(log
, "xi: %8.3f %8.3f %8.3f\n", XI
, YI
, ZI
);
1938 for (dx
= dx0
; (dx
<= dx1
); dx
++)
1940 tmp1
= rs2
- dcx2
[dx
];
1941 for (dy
= dy0
; (dy
<= dy1
); dy
++)
1943 tmp2
= tmp1
- dcy2
[dy
];
1946 for (dz
= dz0
; (dz
<= dz1
); dz
++)
1948 if (tmp2
> dcz2
[dz
])
1950 /* Find grid-cell cj in which possible neighbours are */
1951 cj
= xyz2ci(Ny
, Nz
, dx
, dy
, dz
);
1953 /* Check out how many cgs (nrj) there in this cell */
1956 /* Find the offset in the cg list */
1959 /* Check if all j's are out of range so we
1960 * can skip the whole cell.
1961 * Should save some time, especially with DD.
1964 (grida
[cgj0
] >= max_jcg
&&
1965 (grida
[cgj0
] >= jcg1
|| grida
[cgj0
+nrj
-1] < jcg0
)))
1971 for (j
= 0; (j
< nrj
); j
++)
1973 jjcg
= grida
[cgj0
+j
];
1975 /* check whether this guy is in range! */
1976 if ((jjcg
>= jcg0
&& jjcg
< jcg1
) ||
1979 r2
= calc_dx2(XI
, YI
, ZI
, cgcm
[jjcg
]);
1982 /* jgid = gid[cgsatoms[cgsindex[jjcg]]]; */
1983 jgid
= GET_CGINFO_GID(cginfo
[jjcg
]);
1984 /* check energy group exclusions */
1985 if (!(i_egp_flags
[jgid
] & EGP_EXCL
))
1987 if (nsr
[jgid
] >= MAX_CG
)
1989 /* Add to short-range list */
1990 put_in_list(bHaveVdW
, ngid
, md
, icg
, jgid
,
1991 nsr
[jgid
], nl_sr
[jgid
],
1992 cgs
->index
, /* cgsatoms, */ bexcl
,
1993 shift
, fr
, TRUE
, TRUE
, fr
->solvent_opt
);
1996 nl_sr
[jgid
][nsr
[jgid
]++] = jjcg
;
2007 /* CHECK whether there is anything left in the buffers */
2008 for (nn
= 0; (nn
< ngid
); nn
++)
2012 put_in_list(bHaveVdW
, ngid
, md
, icg
, nn
, nsr
[nn
], nl_sr
[nn
],
2013 cgs
->index
, /* cgsatoms, */ bexcl
,
2014 shift
, fr
, TRUE
, TRUE
, fr
->solvent_opt
);
2020 setexcl(cgs
->index
[icg
], cgs
->index
[icg
+1], &top
->excls
, FALSE
, bexcl
);
2022 /* No need to perform any left-over force calculations anymore (as we used to do here)
2023 * since we now save the proper long-range lists for later evaluation.
2026 /* Close neighbourlists */
2027 close_neighbor_lists(fr
, bMakeQMMMnblist
);
2032 void ns_realloc_natoms(gmx_ns_t
*ns
, int natoms
)
2036 if (natoms
> ns
->nra_alloc
)
2038 ns
->nra_alloc
= over_alloc_dd(natoms
);
2039 srenew(ns
->bexcl
, ns
->nra_alloc
);
2040 for (i
= 0; i
< ns
->nra_alloc
; i
++)
2047 void init_ns(FILE *fplog
, const t_commrec
*cr
,
2048 gmx_ns_t
*ns
, t_forcerec
*fr
,
2049 const gmx_mtop_t
*mtop
)
2051 int mt
, icg
, nr_in_cg
, maxcg
, i
, j
, jcg
, ngid
, ncg
;
2054 /* Compute largest charge groups size (# atoms) */
2056 for (mt
= 0; mt
< mtop
->nmoltype
; mt
++)
2058 cgs
= &mtop
->moltype
[mt
].cgs
;
2059 for (icg
= 0; (icg
< cgs
->nr
); icg
++)
2061 nr_in_cg
= std::max(nr_in_cg
, (int)(cgs
->index
[icg
+1]-cgs
->index
[icg
]));
2065 /* Verify whether largest charge group is <= max cg.
2066 * This is determined by the type of the local exclusion type
2067 * Exclusions are stored in bits. (If the type is not large
2068 * enough, enlarge it, unsigned char -> unsigned short -> unsigned long)
2070 maxcg
= sizeof(t_excl
)*8;
2071 if (nr_in_cg
> maxcg
)
2073 gmx_fatal(FARGS
, "Max #atoms in a charge group: %d > %d\n",
2077 ngid
= mtop
->groups
.grps
[egcENER
].nr
;
2078 snew(ns
->bExcludeAlleg
, ngid
);
2079 for (i
= 0; i
< ngid
; i
++)
2081 ns
->bExcludeAlleg
[i
] = TRUE
;
2082 for (j
= 0; j
< ngid
; j
++)
2084 if (!(fr
->egp_flags
[i
*ngid
+j
] & EGP_EXCL
))
2086 ns
->bExcludeAlleg
[i
] = FALSE
;
2094 ns
->grid
= init_grid(fplog
, fr
);
2095 init_nsgrid_lists(fr
, ngid
, ns
);
2100 snew(ns
->ns_buf
, ngid
);
2101 for (i
= 0; (i
< ngid
); i
++)
2103 snew(ns
->ns_buf
[i
], SHIFTS
);
2105 ncg
= ncg_mtop(mtop
);
2106 snew(ns
->simple_aaj
, 2*ncg
);
2107 for (jcg
= 0; (jcg
< ncg
); jcg
++)
2109 ns
->simple_aaj
[jcg
] = jcg
;
2110 ns
->simple_aaj
[jcg
+ncg
] = jcg
;
2114 /* Create array that determines whether or not atoms have VdW */
2115 snew(ns
->bHaveVdW
, fr
->ntype
);
2116 for (i
= 0; (i
< fr
->ntype
); i
++)
2118 for (j
= 0; (j
< fr
->ntype
); j
++)
2120 ns
->bHaveVdW
[i
] = (ns
->bHaveVdW
[i
] ||
2122 ((BHAMA(fr
->nbfp
, fr
->ntype
, i
, j
) != 0) ||
2123 (BHAMB(fr
->nbfp
, fr
->ntype
, i
, j
) != 0) ||
2124 (BHAMC(fr
->nbfp
, fr
->ntype
, i
, j
) != 0)) :
2125 ((C6(fr
->nbfp
, fr
->ntype
, i
, j
) != 0) ||
2126 (C12(fr
->nbfp
, fr
->ntype
, i
, j
) != 0))));
2131 pr_bvec(debug
, 0, "bHaveVdW", ns
->bHaveVdW
, fr
->ntype
, TRUE
);
2135 ns
->bexcl
= nullptr;
2136 if (!DOMAINDECOMP(cr
))
2138 ns_realloc_natoms(ns
, mtop
->natoms
);
2141 ns
->nblist_initialized
= FALSE
;
2143 /* nbr list debug dump */
2145 char *ptr
= getenv("GMX_DUMP_NL");
2148 ns
->dump_nl
= strtol(ptr
, nullptr, 10);
2151 fprintf(fplog
, "GMX_DUMP_NL = %d", ns
->dump_nl
);
2162 int search_neighbours(FILE *log
, t_forcerec
*fr
,
2164 gmx_localtop_t
*top
,
2165 gmx_groups_t
*groups
,
2167 t_nrnb
*nrnb
, t_mdatoms
*md
,
2170 t_block
*cgs
= &(top
->cgs
);
2171 rvec box_size
, grid_x0
, grid_x1
;
2173 real min_size
, grid_dens
;
2179 gmx_domdec_zones_t
*dd_zones
;
2180 put_in_list_t
*put_in_list
;
2184 /* Set some local variables */
2186 ngid
= groups
->grps
[egcENER
].nr
;
2188 for (m
= 0; (m
< DIM
); m
++)
2190 box_size
[m
] = box
[m
][m
];
2193 if (fr
->ePBC
!= epbcNONE
)
2195 if (gmx::square(fr
->rlist
) >= max_cutoff2(fr
->ePBC
, box
))
2197 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.");
2201 min_size
= std::min(box_size
[XX
], std::min(box_size
[YY
], box_size
[ZZ
]));
2202 if (2*fr
->rlist
>= min_size
)
2204 gmx_fatal(FARGS
, "One of the box diagonal elements has become smaller than twice the cut-off length.");
2209 if (DOMAINDECOMP(cr
))
2211 ns_realloc_natoms(ns
, cgs
->index
[cgs
->nr
]);
2214 /* Reset the neighbourlists */
2215 reset_neighbor_lists(fr
);
2217 if (bGrid
&& bFillGrid
)
2221 if (DOMAINDECOMP(cr
))
2223 dd_zones
= domdec_zones(cr
->dd
);
2229 get_nsgrid_boundaries(grid
->nboundeddim
, box
, nullptr, nullptr, nullptr, nullptr,
2230 cgs
->nr
, fr
->cg_cm
, grid_x0
, grid_x1
, &grid_dens
);
2232 grid_first(log
, grid
, nullptr, nullptr, box
, grid_x0
, grid_x1
,
2233 fr
->rlist
, grid_dens
);
2239 if (DOMAINDECOMP(cr
))
2242 fill_grid(dd_zones
, grid
, end
, -1, end
, fr
->cg_cm
);
2244 grid
->icg1
= dd_zones
->izone
[dd_zones
->nizone
-1].cg1
;
2248 fill_grid(nullptr, grid
, cgs
->nr
, fr
->cg0
, fr
->hcg
, fr
->cg_cm
);
2249 grid
->icg0
= fr
->cg0
;
2250 grid
->icg1
= fr
->hcg
;
2253 calc_elemnr(grid
, start
, end
, cgs
->nr
);
2255 grid_last(grid
, start
, end
, cgs
->nr
);
2260 print_grid(debug
, grid
);
2265 /* Set the grid cell index for the test particle only.
2266 * The cell to cg index is not corrected, but that does not matter.
2268 fill_grid(nullptr, ns
->grid
, fr
->hcg
, fr
->hcg
-1, fr
->hcg
, fr
->cg_cm
);
2271 if (!fr
->ns
->bCGlist
)
2273 put_in_list
= put_in_list_at
;
2277 put_in_list
= put_in_list_cg
;
2284 nsearch
= nsgrid_core(cr
, fr
, box
, ngid
, top
,
2285 grid
, ns
->bexcl
, ns
->bExcludeAlleg
,
2286 md
, put_in_list
, ns
->bHaveVdW
,
2289 /* neighbour searching withouth QMMM! QM atoms have zero charge in
2290 * the classical calculation. The charge-charge interaction
2291 * between QM and MM atoms is handled in the QMMM core calculation
2292 * (see QMMM.c). The VDW however, we'd like to compute classically
2293 * and the QM MM atom pairs have just been put in the
2294 * corresponding neighbourlists. in case of QMMM we still need to
2295 * fill a special QMMM neighbourlist that contains all neighbours
2296 * of the QM atoms. If bQMMM is true, this list will now be made:
2298 if (fr
->bQMMM
&& fr
->qr
->QMMMscheme
!= eQMMMschemeoniom
)
2300 nsearch
+= nsgrid_core(cr
, fr
, box
, ngid
, top
,
2301 grid
, ns
->bexcl
, ns
->bExcludeAlleg
,
2302 md
, put_in_list_qmmm
, ns
->bHaveVdW
,
2308 nsearch
= ns_simple_core(fr
, top
, md
, box
, box_size
,
2309 ns
->bexcl
, ns
->simple_aaj
,
2310 ngid
, ns
->ns_buf
, put_in_list
, ns
->bHaveVdW
);
2313 inc_nrnb(nrnb
, eNR_NS
, nsearch
);