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,2018, 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.
48 #include "gromacs/domdec/domdec.h"
49 #include "gromacs/domdec/domdec_struct.h"
50 #include "gromacs/gmxlib/network.h"
51 #include "gromacs/gmxlib/nrnb.h"
52 #include "gromacs/gmxlib/nonbonded/nonbonded.h"
53 #include "gromacs/math/functions.h"
54 #include "gromacs/math/utilities.h"
55 #include "gromacs/math/vec.h"
56 #include "gromacs/math/vecdump.h"
57 #include "gromacs/mdlib/force.h"
58 #include "gromacs/mdlib/nsgrid.h"
59 #include "gromacs/mdlib/qmmm.h"
60 #include "gromacs/mdtypes/commrec.h"
61 #include "gromacs/mdtypes/group.h"
62 #include "gromacs/mdtypes/inputrec.h"
63 #include "gromacs/mdtypes/md_enums.h"
64 #include "gromacs/pbcutil/ishift.h"
65 #include "gromacs/pbcutil/pbc.h"
66 #include "gromacs/topology/mtop_util.h"
67 #include "gromacs/utility/fatalerror.h"
68 #include "gromacs/utility/smalloc.h"
71 * E X C L U S I O N H A N D L I N G
75 static void SETEXCL_(t_excl e
[], int i
, int j
)
79 static void RMEXCL_(t_excl e
[], int i
, int j
)
81 e
[j
] = e
[j
] & ~(1<<i
);
83 static gmx_bool
ISEXCL_(t_excl e
[], int i
, int j
)
85 return (gmx_bool
)(e
[j
] & (1<<i
));
87 static gmx_bool
NOTEXCL_(t_excl e
[], int i
, int j
)
89 return !(ISEXCL(e
, i
, j
));
92 #define SETEXCL(e, i, j) (e)[((int) (j))] |= (1<<((int) (i)))
93 #define RMEXCL(e, i, j) (e)[((int) (j))] &= (~(1<<((int) (i))))
94 #define ISEXCL(e, i, j) (gmx_bool) ((e)[((int) (j))] & (1<<((int) (i))))
95 #define NOTEXCL(e, i, j) !(ISEXCL(e, i, j))
99 round_up_to_simd_width(int length
, int simd_width
)
103 offset
= (simd_width
> 0) ? length
% simd_width
: 0;
105 return (offset
== 0) ? length
: length
-offset
+simd_width
;
107 /************************************************
109 * U T I L I T I E S F O R N S
111 ************************************************/
113 void reallocate_nblist(t_nblist
*nl
)
117 fprintf(debug
, "reallocating neigborlist (ielec=%d, ivdw=%d, igeometry=%d, type=%d), maxnri=%d\n",
118 nl
->ielec
, nl
->ivdw
, nl
->igeometry
, nl
->type
, nl
->maxnri
);
120 srenew(nl
->iinr
, nl
->maxnri
);
121 if (nl
->igeometry
== GMX_NBLIST_GEOMETRY_CG_CG
)
123 srenew(nl
->iinr_end
, nl
->maxnri
);
125 srenew(nl
->gid
, nl
->maxnri
);
126 srenew(nl
->shift
, nl
->maxnri
);
127 srenew(nl
->jindex
, nl
->maxnri
+1);
131 static void init_nblist(FILE *log
, t_nblist
*nl_sr
,
133 int ivdw
, int ivdwmod
,
134 int ielec
, int ielecmod
,
135 int igeometry
, int type
,
136 gmx_bool bElecAndVdwSwitchDiffers
)
151 /* Set coul/vdw in neighborlist, and for the normal loops we determine
152 * an index of which one to call.
155 nl
->ivdwmod
= ivdwmod
;
157 nl
->ielecmod
= ielecmod
;
159 nl
->igeometry
= igeometry
;
161 if (nl
->type
== GMX_NBLIST_INTERACTION_FREE_ENERGY
)
163 nl
->igeometry
= GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE
;
166 /* This will also set the simd_padding_width field */
167 gmx_nonbonded_set_kernel_pointers(log
, nl
, bElecAndVdwSwitchDiffers
);
169 /* maxnri is influenced by the number of shifts (maximum is 8)
170 * and the number of energy groups.
171 * If it is not enough, nl memory will be reallocated during the run.
172 * 4 seems to be a reasonable factor, which only causes reallocation
173 * during runs with tiny and many energygroups.
175 nl
->maxnri
= homenr
*4;
182 nl
->jindex
= nullptr;
184 nl
->excl_fep
= nullptr;
185 reallocate_nblist(nl
);
190 fprintf(debug
, "Initiating neighbourlist (ielec=%d, ivdw=%d, type=%d) for %s interactions,\nwith %d SR atoms.\n",
191 nl
->ielec
, nl
->ivdw
, nl
->type
, gmx_nblist_geometry_names
[nl
->igeometry
], maxsr
);
196 void init_neighbor_list(FILE *log
, t_forcerec
*fr
, int homenr
)
198 int maxsr
, maxsr_wat
;
199 int ielec
, ivdw
, ielecmod
, ivdwmod
, type
;
200 int igeometry_def
, igeometry_w
, igeometry_ww
;
202 gmx_bool bElecAndVdwSwitchDiffers
;
205 /* maxsr = homenr-fr->nWatMol*3; */
210 gmx_fatal(FARGS
, "%s, %d: Negative number of short range atoms.\n"
211 "Call your GROMACS dealer for assistance.", __FILE__
, __LINE__
);
213 /* This is just for initial allocation, so we do not reallocate
214 * all the nlist arrays many times in a row.
215 * The numbers seem very accurate, but they are uncritical.
217 maxsr_wat
= std::min(fr
->nWatMol
, (homenr
+2)/3);
219 /* Determine the values for ielec/ivdw. */
220 ielec
= fr
->nbkernel_elec_interaction
;
221 ivdw
= fr
->nbkernel_vdw_interaction
;
222 ielecmod
= fr
->nbkernel_elec_modifier
;
223 ivdwmod
= fr
->nbkernel_vdw_modifier
;
224 type
= GMX_NBLIST_INTERACTION_STANDARD
;
225 bElecAndVdwSwitchDiffers
= ( (fr
->ic
->rcoulomb_switch
!= fr
->ic
->rvdw_switch
) || (fr
->ic
->rcoulomb
!= fr
->ic
->rvdw
));
227 fr
->ns
->bCGlist
= (getenv("GMX_NBLISTCG") != nullptr);
228 if (!fr
->ns
->bCGlist
)
230 igeometry_def
= GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE
;
234 igeometry_def
= GMX_NBLIST_GEOMETRY_CG_CG
;
237 fprintf(log
, "\nUsing charge-group - charge-group neighbor lists and kernels\n\n");
241 if (fr
->solvent_opt
== esolTIP4P
)
243 igeometry_w
= GMX_NBLIST_GEOMETRY_WATER4_PARTICLE
;
244 igeometry_ww
= GMX_NBLIST_GEOMETRY_WATER4_WATER4
;
248 igeometry_w
= GMX_NBLIST_GEOMETRY_WATER3_PARTICLE
;
249 igeometry_ww
= GMX_NBLIST_GEOMETRY_WATER3_WATER3
;
252 for (i
= 0; i
< fr
->nnblists
; i
++)
254 nbl
= &(fr
->nblists
[i
]);
256 init_nblist(log
, &nbl
->nlist_sr
[eNL_VDWQQ
],
257 maxsr
, ivdw
, ivdwmod
, ielec
, ielecmod
, igeometry_def
, type
, bElecAndVdwSwitchDiffers
);
258 init_nblist(log
, &nbl
->nlist_sr
[eNL_VDW
],
259 maxsr
, ivdw
, ivdwmod
, GMX_NBKERNEL_ELEC_NONE
, eintmodNONE
, igeometry_def
, type
, bElecAndVdwSwitchDiffers
);
260 init_nblist(log
, &nbl
->nlist_sr
[eNL_QQ
],
261 maxsr
, GMX_NBKERNEL_VDW_NONE
, eintmodNONE
, ielec
, ielecmod
, igeometry_def
, type
, bElecAndVdwSwitchDiffers
);
262 init_nblist(log
, &nbl
->nlist_sr
[eNL_VDWQQ_WATER
],
263 maxsr_wat
, ivdw
, ivdwmod
, ielec
, ielecmod
, igeometry_w
, type
, bElecAndVdwSwitchDiffers
);
264 init_nblist(log
, &nbl
->nlist_sr
[eNL_QQ_WATER
],
265 maxsr_wat
, GMX_NBKERNEL_VDW_NONE
, eintmodNONE
, ielec
, ielecmod
, igeometry_w
, type
, bElecAndVdwSwitchDiffers
);
266 init_nblist(log
, &nbl
->nlist_sr
[eNL_VDWQQ_WATERWATER
],
267 maxsr_wat
, ivdw
, ivdwmod
, ielec
, ielecmod
, igeometry_ww
, type
, bElecAndVdwSwitchDiffers
);
268 init_nblist(log
, &nbl
->nlist_sr
[eNL_QQ_WATERWATER
],
269 maxsr_wat
, GMX_NBKERNEL_VDW_NONE
, eintmodNONE
, ielec
, ielecmod
, igeometry_ww
, type
, bElecAndVdwSwitchDiffers
);
271 /* Did we get the solvent loops so we can use optimized water kernels? */
272 if (nbl
->nlist_sr
[eNL_VDWQQ_WATER
].kernelptr_vf
== nullptr
273 || nbl
->nlist_sr
[eNL_QQ_WATER
].kernelptr_vf
== nullptr
274 || nbl
->nlist_sr
[eNL_VDWQQ_WATERWATER
].kernelptr_vf
== nullptr
275 || nbl
->nlist_sr
[eNL_QQ_WATERWATER
].kernelptr_vf
== nullptr)
277 fr
->solvent_opt
= esolNO
;
280 fprintf(log
, "Note: The available nonbonded kernels do not support water optimization - disabling.\n");
284 if (fr
->efep
!= efepNO
)
286 init_nblist(log
, &nbl
->nlist_sr
[eNL_VDWQQ_FREE
],
287 maxsr
, ivdw
, ivdwmod
, ielec
, ielecmod
, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE
, GMX_NBLIST_INTERACTION_FREE_ENERGY
, bElecAndVdwSwitchDiffers
);
288 init_nblist(log
, &nbl
->nlist_sr
[eNL_VDW_FREE
],
289 maxsr
, ivdw
, ivdwmod
, GMX_NBKERNEL_ELEC_NONE
, eintmodNONE
, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE
, GMX_NBLIST_INTERACTION_FREE_ENERGY
, bElecAndVdwSwitchDiffers
);
290 init_nblist(log
, &nbl
->nlist_sr
[eNL_QQ_FREE
],
291 maxsr
, GMX_NBKERNEL_VDW_NONE
, eintmodNONE
, ielec
, ielecmod
, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE
, GMX_NBLIST_INTERACTION_FREE_ENERGY
, bElecAndVdwSwitchDiffers
);
295 if (fr
->bQMMM
&& fr
->qr
->QMMMscheme
!= eQMMMschemeoniom
)
297 if (nullptr == fr
->QMMMlist
)
299 snew(fr
->QMMMlist
, 1);
301 init_nblist(log
, fr
->QMMMlist
,
302 maxsr
, 0, 0, ielec
, ielecmod
, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE
, GMX_NBLIST_INTERACTION_STANDARD
, bElecAndVdwSwitchDiffers
);
310 fr
->ns
->nblist_initialized
= TRUE
;
313 static void reset_nblist(t_nblist
*nl
)
323 static void reset_neighbor_lists(t_forcerec
*fr
)
329 /* only reset the short-range nblist */
330 reset_nblist(fr
->QMMMlist
);
333 for (n
= 0; n
< fr
->nnblists
; n
++)
335 for (i
= 0; i
< eNL_NR
; i
++)
337 reset_nblist( &(fr
->nblists
[n
].nlist_sr
[i
]) );
345 static inline void new_i_nblist(t_nblist
*nlist
, int i_atom
, int shift
, int gid
)
347 int nri
= nlist
->nri
;
349 /* Check whether we have to increase the i counter */
351 (nlist
->iinr
[nri
] != i_atom
) ||
352 (nlist
->shift
[nri
] != shift
) ||
353 (nlist
->gid
[nri
] != gid
))
355 /* This is something else. Now see if any entries have
356 * been added in the list of the previous atom.
359 ((nlist
->jindex
[nri
+1] > nlist
->jindex
[nri
]) &&
360 (nlist
->gid
[nri
] != -1)))
362 /* If so increase the counter */
365 if (nlist
->nri
>= nlist
->maxnri
)
367 nlist
->maxnri
+= over_alloc_large(nlist
->nri
);
368 reallocate_nblist(nlist
);
371 /* Set the number of neighbours and the atom number */
372 nlist
->jindex
[nri
+1] = nlist
->jindex
[nri
];
373 nlist
->iinr
[nri
] = i_atom
;
374 nlist
->gid
[nri
] = gid
;
375 nlist
->shift
[nri
] = shift
;
379 /* Adding to previous list. First remove possible previous padding */
380 if (nlist
->simd_padding_width
> 1)
382 while (nlist
->nrj
> 0 && nlist
->jjnr
[nlist
->nrj
-1] < 0)
390 static inline void close_i_nblist(t_nblist
*nlist
)
392 int nri
= nlist
->nri
;
397 /* Add elements up to padding. Since we allocate memory in units
398 * of the simd_padding width, we do not have to check for possible
399 * list reallocation here.
401 while ((nlist
->nrj
% nlist
->simd_padding_width
) != 0)
403 /* Use -4 here, so we can write forces for 4 atoms before real data */
404 nlist
->jjnr
[nlist
->nrj
++] = -4;
406 nlist
->jindex
[nri
+1] = nlist
->nrj
;
408 len
= nlist
->nrj
- nlist
->jindex
[nri
];
409 /* If there are no j-particles we have to reduce the
410 * number of i-particles again, to prevent errors in the
413 if ((len
== 0) && (nlist
->nri
> 0))
420 static inline void close_nblist(t_nblist
*nlist
)
422 /* Only close this nblist when it has been initialized.
423 * Avoid the creation of i-lists with no j-particles.
427 /* Some assembly kernels do not support empty lists,
428 * make sure here that we don't generate any empty lists.
429 * With the current ns code this branch is taken in two cases:
430 * No i-particles at all: nri=-1 here
431 * There are i-particles, but no j-particles; nri=0 here
437 /* Close list number nri by incrementing the count */
442 static inline void close_neighbor_lists(t_forcerec
*fr
, gmx_bool bMakeQMMMnblist
)
448 close_nblist(fr
->QMMMlist
);
451 for (n
= 0; n
< fr
->nnblists
; n
++)
453 for (i
= 0; (i
< eNL_NR
); i
++)
455 close_nblist(&(fr
->nblists
[n
].nlist_sr
[i
]));
461 static inline void add_j_to_nblist(t_nblist
*nlist
, int j_atom
)
463 int nrj
= nlist
->nrj
;
465 if (nlist
->nrj
>= nlist
->maxnrj
)
467 nlist
->maxnrj
= round_up_to_simd_width(over_alloc_small(nlist
->nrj
+ 1), nlist
->simd_padding_width
);
471 fprintf(debug
, "Increasing SR nblist (ielec=%d,ivdw=%d,type=%d,igeometry=%d) j size to %d\n",
472 nlist
->ielec
, nlist
->ivdw
, nlist
->type
, nlist
->igeometry
, nlist
->maxnrj
);
475 srenew(nlist
->jjnr
, nlist
->maxnrj
);
478 nlist
->jjnr
[nrj
] = j_atom
;
482 static inline void add_j_to_nblist_cg(t_nblist
*nlist
,
483 int j_start
, int j_end
,
484 t_excl
*bexcl
, gmx_bool i_is_j
)
486 int nrj
= nlist
->nrj
;
489 if (nlist
->nrj
>= nlist
->maxnrj
)
491 nlist
->maxnrj
= over_alloc_small(nlist
->nrj
+ 1);
494 fprintf(debug
, "Increasing SR nblist (ielec=%d,ivdw=%d,type=%d,igeometry=%d) j size to %d\n",
495 nlist
->ielec
, nlist
->ivdw
, nlist
->type
, nlist
->igeometry
, nlist
->maxnrj
);
498 srenew(nlist
->jjnr
, nlist
->maxnrj
);
499 srenew(nlist
->jjnr_end
, nlist
->maxnrj
);
500 srenew(nlist
->excl
, nlist
->maxnrj
*MAX_CGCGSIZE
);
503 nlist
->jjnr
[nrj
] = j_start
;
504 nlist
->jjnr_end
[nrj
] = j_end
;
506 if (j_end
- j_start
> MAX_CGCGSIZE
)
508 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
);
511 /* Set the exclusions */
512 for (j
= j_start
; j
< j_end
; j
++)
514 nlist
->excl
[nrj
*MAX_CGCGSIZE
+ j
- j_start
] = bexcl
[j
];
518 /* Avoid double counting of intra-cg interactions */
519 for (j
= 1; j
< j_end
-j_start
; j
++)
521 nlist
->excl
[nrj
*MAX_CGCGSIZE
+ j
] |= (1<<j
) - 1;
529 put_in_list_t (gmx_bool bHaveVdW
[],
545 put_in_list_at(gmx_bool bHaveVdW
[],
560 /* The a[] index has been removed,
561 * to put it back in i_atom should be a[i0] and jj should be a[jj].
566 t_nblist
* vdwc_free
= nullptr;
567 t_nblist
* vdw_free
= nullptr;
568 t_nblist
* coul_free
= nullptr;
569 t_nblist
* vdwc_ww
= nullptr;
570 t_nblist
* coul_ww
= nullptr;
572 int i
, j
, jcg
, igid
, gid
, nbl_ind
;
573 int jj
, jj0
, jj1
, i_atom
;
578 real
*charge
, *chargeB
;
580 gmx_bool bFreeEnergy
, bFree
, bFreeJ
, bNotEx
, *bPert
;
581 gmx_bool bDoVdW_i
, bDoCoul_i
, bDoCoul_i_sol
;
585 /* Copy some pointers */
587 charge
= md
->chargeA
;
588 chargeB
= md
->chargeB
;
591 bPert
= md
->bPerturbed
;
595 nicg
= index
[icg
+1]-i0
;
597 /* Get the i charge group info */
598 igid
= GET_CGINFO_GID(cginfo
[icg
]);
600 iwater
= (solvent_opt
!= esolNO
) ? GET_CGINFO_SOLOPT(cginfo
[icg
]) : esolNO
;
605 /* Check if any of the particles involved are perturbed.
606 * If not we can do the cheaper normal put_in_list
607 * and use more solvent optimization.
609 for (i
= 0; i
< nicg
; i
++)
611 bFreeEnergy
|= bPert
[i0
+i
];
613 /* Loop over the j charge groups */
614 for (j
= 0; (j
< nj
&& !bFreeEnergy
); j
++)
619 /* Finally loop over the atoms in the j-charge group */
620 for (jj
= jj0
; jj
< jj1
; jj
++)
622 bFreeEnergy
|= bPert
[jj
];
627 /* Unpack pointers to neighbourlist structs */
628 if (fr
->nnblists
== 1)
634 nbl_ind
= fr
->gid2nblists
[GID(igid
, jgid
, ngid
)];
636 nlist
= fr
->nblists
[nbl_ind
].nlist_sr
;
638 if (iwater
!= esolNO
)
640 vdwc
= &nlist
[eNL_VDWQQ_WATER
];
641 vdw
= &nlist
[eNL_VDW
];
642 coul
= &nlist
[eNL_QQ_WATER
];
643 vdwc_ww
= &nlist
[eNL_VDWQQ_WATERWATER
];
644 coul_ww
= &nlist
[eNL_QQ_WATERWATER
];
648 vdwc
= &nlist
[eNL_VDWQQ
];
649 vdw
= &nlist
[eNL_VDW
];
650 coul
= &nlist
[eNL_QQ
];
655 if (iwater
!= esolNO
)
657 /* Loop over the atoms in the i charge group */
659 gid
= GID(igid
, jgid
, ngid
);
660 /* Create new i_atom for each energy group */
661 if (bDoCoul
&& bDoVdW
)
663 new_i_nblist(vdwc
, i_atom
, shift
, gid
);
664 new_i_nblist(vdwc_ww
, i_atom
, shift
, gid
);
668 new_i_nblist(vdw
, i_atom
, shift
, gid
);
672 new_i_nblist(coul
, i_atom
, shift
, gid
);
673 new_i_nblist(coul_ww
, i_atom
, shift
, gid
);
675 /* Loop over the j charge groups */
676 for (j
= 0; (j
< nj
); j
++)
686 jwater
= GET_CGINFO_SOLOPT(cginfo
[jcg
]);
688 if (iwater
== esolSPC
&& jwater
== esolSPC
)
690 /* Interaction between two SPC molecules */
693 /* VdW only - only first atoms in each water interact */
694 add_j_to_nblist(vdw
, jj0
);
698 /* One entry for the entire water-water interaction */
701 add_j_to_nblist(coul_ww
, jj0
);
705 add_j_to_nblist(vdwc_ww
, jj0
);
709 else if (iwater
== esolTIP4P
&& jwater
== esolTIP4P
)
711 /* Interaction between two TIP4p molecules */
714 /* VdW only - only first atoms in each water interact */
715 add_j_to_nblist(vdw
, jj0
);
719 /* One entry for the entire water-water interaction */
722 add_j_to_nblist(coul_ww
, jj0
);
726 add_j_to_nblist(vdwc_ww
, jj0
);
732 /* j charge group is not water, but i is.
733 * Add entries to the water-other_atom lists; the geometry of the water
734 * molecule doesn't matter - that is taken care of in the nonbonded kernel,
735 * so we don't care if it is SPC or TIP4P...
742 for (jj
= jj0
; (jj
< jj1
); jj
++)
746 add_j_to_nblist(coul
, jj
);
752 for (jj
= jj0
; (jj
< jj1
); jj
++)
754 if (bHaveVdW
[type
[jj
]])
756 add_j_to_nblist(vdw
, jj
);
762 /* _charge_ _groups_ interact with both coulomb and LJ */
763 /* Check which atoms we should add to the lists! */
764 for (jj
= jj0
; (jj
< jj1
); jj
++)
766 if (bHaveVdW
[type
[jj
]])
770 add_j_to_nblist(vdwc
, jj
);
774 add_j_to_nblist(vdw
, jj
);
777 else if (charge
[jj
] != 0)
779 add_j_to_nblist(coul
, jj
);
786 close_i_nblist(coul
);
787 close_i_nblist(vdwc
);
788 close_i_nblist(coul_ww
);
789 close_i_nblist(vdwc_ww
);
793 /* no solvent as i charge group */
794 /* Loop over the atoms in the i charge group */
795 for (i
= 0; i
< nicg
; i
++)
798 gid
= GID(igid
, jgid
, ngid
);
801 /* Create new i_atom for each energy group */
802 if (bDoVdW
&& bDoCoul
)
804 new_i_nblist(vdwc
, i_atom
, shift
, gid
);
808 new_i_nblist(vdw
, i_atom
, shift
, gid
);
812 new_i_nblist(coul
, i_atom
, shift
, gid
);
814 bDoVdW_i
= (bDoVdW
&& bHaveVdW
[type
[i_atom
]]);
815 bDoCoul_i
= (bDoCoul
&& qi
!= 0);
817 if (bDoVdW_i
|| bDoCoul_i
)
819 /* Loop over the j charge groups */
820 for (j
= 0; (j
< nj
); j
++)
824 /* Check for large charge groups */
835 /* Finally loop over the atoms in the j-charge group */
836 for (jj
= jj0
; jj
< jj1
; jj
++)
838 bNotEx
= NOTEXCL(bExcl
, i
, jj
);
846 add_j_to_nblist(coul
, jj
);
851 if (bHaveVdW
[type
[jj
]])
853 add_j_to_nblist(vdw
, jj
);
858 if (bHaveVdW
[type
[jj
]])
862 add_j_to_nblist(vdwc
, jj
);
866 add_j_to_nblist(vdw
, jj
);
869 else if (charge
[jj
] != 0)
871 add_j_to_nblist(coul
, jj
);
879 close_i_nblist(coul
);
880 close_i_nblist(vdwc
);
886 /* we are doing free energy */
887 vdwc_free
= &nlist
[eNL_VDWQQ_FREE
];
888 vdw_free
= &nlist
[eNL_VDW_FREE
];
889 coul_free
= &nlist
[eNL_QQ_FREE
];
890 /* Loop over the atoms in the i charge group */
891 for (i
= 0; i
< nicg
; i
++)
894 gid
= GID(igid
, jgid
, ngid
);
896 qiB
= chargeB
[i_atom
];
898 /* Create new i_atom for each energy group */
899 if (bDoVdW
&& bDoCoul
)
901 new_i_nblist(vdwc
, i_atom
, shift
, gid
);
905 new_i_nblist(vdw
, i_atom
, shift
, gid
);
909 new_i_nblist(coul
, i_atom
, shift
, gid
);
912 new_i_nblist(vdw_free
, i_atom
, shift
, gid
);
913 new_i_nblist(coul_free
, i_atom
, shift
, gid
);
914 new_i_nblist(vdwc_free
, i_atom
, shift
, gid
);
916 bDoVdW_i
= (bDoVdW
&&
917 (bHaveVdW
[type
[i_atom
]] || bHaveVdW
[typeB
[i_atom
]]));
918 bDoCoul_i
= (bDoCoul
&& (qi
!= 0 || qiB
!= 0));
919 /* For TIP4P the first atom does not have a charge,
920 * but the last three do. So we should still put an atom
921 * without LJ but with charge in the water-atom neighborlist
922 * for a TIP4p i charge group.
923 * For SPC type water the first atom has LJ and charge,
924 * so there is no such problem.
926 if (iwater
== esolNO
)
928 bDoCoul_i_sol
= bDoCoul_i
;
932 bDoCoul_i_sol
= bDoCoul
;
935 if (bDoVdW_i
|| bDoCoul_i_sol
)
937 /* Loop over the j charge groups */
938 for (j
= 0; (j
< nj
); j
++)
942 /* Check for large charge groups */
953 /* Finally loop over the atoms in the j-charge group */
954 bFree
= bPert
[i_atom
];
955 for (jj
= jj0
; (jj
< jj1
); jj
++)
957 bFreeJ
= bFree
|| bPert
[jj
];
958 /* Complicated if, because the water H's should also
959 * see perturbed j-particles
961 if (iwater
== esolNO
|| i
== 0 || bFreeJ
)
963 bNotEx
= NOTEXCL(bExcl
, i
, jj
);
971 if (charge
[jj
] != 0 || chargeB
[jj
] != 0)
973 add_j_to_nblist(coul_free
, jj
);
978 if (bHaveVdW
[type
[jj
]] || bHaveVdW
[typeB
[jj
]])
980 add_j_to_nblist(vdw_free
, jj
);
985 if (bHaveVdW
[type
[jj
]] || bHaveVdW
[typeB
[jj
]])
987 if (charge
[jj
] != 0 || chargeB
[jj
] != 0)
989 add_j_to_nblist(vdwc_free
, jj
);
993 add_j_to_nblist(vdw_free
, jj
);
996 else if (charge
[jj
] != 0 || chargeB
[jj
] != 0)
998 add_j_to_nblist(coul_free
, jj
);
1004 /* This is done whether or not bWater is set */
1005 if (charge
[jj
] != 0)
1007 add_j_to_nblist(coul
, jj
);
1010 else if (!bDoCoul_i_sol
)
1012 if (bHaveVdW
[type
[jj
]])
1014 add_j_to_nblist(vdw
, jj
);
1019 if (bHaveVdW
[type
[jj
]])
1021 if (charge
[jj
] != 0)
1023 add_j_to_nblist(vdwc
, jj
);
1027 add_j_to_nblist(vdw
, jj
);
1030 else if (charge
[jj
] != 0)
1032 add_j_to_nblist(coul
, jj
);
1040 close_i_nblist(vdw
);
1041 close_i_nblist(coul
);
1042 close_i_nblist(vdwc
);
1043 close_i_nblist(vdw_free
);
1044 close_i_nblist(coul_free
);
1045 close_i_nblist(vdwc_free
);
1051 put_in_list_qmmm(gmx_bool gmx_unused bHaveVdW
[],
1053 t_mdatoms gmx_unused
* md
,
1062 gmx_bool gmx_unused bDoVdW
,
1063 gmx_bool gmx_unused bDoCoul
,
1064 int gmx_unused solvent_opt
)
1067 int i
, j
, jcg
, igid
, gid
;
1068 int jj
, jj0
, jj1
, i_atom
;
1072 /* Get atom range */
1074 nicg
= index
[icg
+1]-i0
;
1076 /* Get the i charge group info */
1077 igid
= GET_CGINFO_GID(fr
->cginfo
[icg
]);
1079 coul
= fr
->QMMMlist
;
1081 /* Loop over atoms in the ith charge group */
1082 for (i
= 0; i
< nicg
; i
++)
1085 gid
= GID(igid
, jgid
, ngid
);
1086 /* Create new i_atom for each energy group */
1087 new_i_nblist(coul
, i_atom
, shift
, gid
);
1089 /* Loop over the j charge groups */
1090 for (j
= 0; j
< nj
; j
++)
1094 /* Charge groups cannot have QM and MM atoms simultaneously */
1099 /* Finally loop over the atoms in the j-charge group */
1100 for (jj
= jj0
; jj
< jj1
; jj
++)
1102 bNotEx
= NOTEXCL(bExcl
, i
, jj
);
1105 add_j_to_nblist(coul
, jj
);
1110 close_i_nblist(coul
);
1115 put_in_list_cg(gmx_bool gmx_unused bHaveVdW
[],
1117 t_mdatoms gmx_unused
* md
,
1126 gmx_bool gmx_unused bDoVdW
,
1127 gmx_bool gmx_unused bDoCoul
,
1128 int gmx_unused solvent_opt
)
1131 int igid
, gid
, nbl_ind
;
1135 cginfo
= fr
->cginfo
[icg
];
1137 igid
= GET_CGINFO_GID(cginfo
);
1138 gid
= GID(igid
, jgid
, ngid
);
1140 /* Unpack pointers to neighbourlist structs */
1141 if (fr
->nnblists
== 1)
1147 nbl_ind
= fr
->gid2nblists
[gid
];
1149 vdwc
= &fr
->nblists
[nbl_ind
].nlist_sr
[eNL_VDWQQ
];
1151 /* Make a new neighbor list for charge group icg.
1152 * Currently simply one neighbor list is made with LJ and Coulomb.
1153 * If required, zero interactions could be removed here
1154 * or in the force loop.
1156 new_i_nblist(vdwc
, index
[icg
], shift
, gid
);
1157 vdwc
->iinr_end
[vdwc
->nri
] = index
[icg
+1];
1159 for (j
= 0; (j
< nj
); j
++)
1162 /* Skip the icg-icg pairs if all self interactions are excluded */
1163 if (!(jcg
== icg
&& GET_CGINFO_EXCL_INTRA(cginfo
)))
1165 /* Here we add the j charge group jcg to the list,
1166 * exclusions are also added to the list.
1168 add_j_to_nblist_cg(vdwc
, index
[jcg
], index
[jcg
+1], bExcl
, icg
== jcg
);
1172 close_i_nblist(vdwc
);
1175 static void setexcl(int start
, int end
, t_blocka
*excl
, gmx_bool b
,
1182 for (i
= start
; i
< end
; i
++)
1184 for (k
= excl
->index
[i
]; k
< excl
->index
[i
+1]; k
++)
1186 SETEXCL(bexcl
, i
-start
, excl
->a
[k
]);
1192 for (i
= start
; i
< end
; i
++)
1194 for (k
= excl
->index
[i
]; k
< excl
->index
[i
+1]; k
++)
1196 RMEXCL(bexcl
, i
-start
, excl
->a
[k
]);
1202 int calc_naaj(int icg
, int cgtot
)
1206 if ((cgtot
% 2) == 1)
1208 /* Odd number of charge groups, easy */
1209 naaj
= 1 + (cgtot
/2);
1211 else if ((cgtot
% 4) == 0)
1213 /* Multiple of four is hard */
1250 fprintf(log
, "naaj=%d\n", naaj
);
1256 /************************************************
1258 * S I M P L E C O R E S T U F F
1260 ************************************************/
1262 static real
calc_image_tric(rvec xi
, rvec xj
, matrix box
,
1263 rvec b_inv
, int *shift
)
1265 /* This code assumes that the cut-off is smaller than
1266 * a half times the smallest diagonal element of the box.
1268 const real h25
= 2.5;
1273 /* Compute diff vector */
1274 dz
= xj
[ZZ
] - xi
[ZZ
];
1275 dy
= xj
[YY
] - xi
[YY
];
1276 dx
= xj
[XX
] - xi
[XX
];
1278 /* Perform NINT operation, using trunc operation, therefore
1279 * we first add 2.5 then subtract 2 again
1281 tz
= static_cast<int>(dz
*b_inv
[ZZ
] + h25
);
1283 dz
-= tz
*box
[ZZ
][ZZ
];
1284 dy
-= tz
*box
[ZZ
][YY
];
1285 dx
-= tz
*box
[ZZ
][XX
];
1287 ty
= static_cast<int>(dy
*b_inv
[YY
] + h25
);
1289 dy
-= ty
*box
[YY
][YY
];
1290 dx
-= ty
*box
[YY
][XX
];
1292 tx
= static_cast<int>(dx
*b_inv
[XX
]+h25
);
1294 dx
-= tx
*box
[XX
][XX
];
1296 /* Distance squared */
1297 r2
= (dx
*dx
) + (dy
*dy
) + (dz
*dz
);
1299 *shift
= XYZ2IS(tx
, ty
, tz
);
1304 static real
calc_image_rect(rvec xi
, rvec xj
, rvec box_size
,
1305 rvec b_inv
, int *shift
)
1307 const real h15
= 1.5;
1313 /* Compute diff vector */
1314 dx
= xj
[XX
] - xi
[XX
];
1315 dy
= xj
[YY
] - xi
[YY
];
1316 dz
= xj
[ZZ
] - xi
[ZZ
];
1318 /* Perform NINT operation, using trunc operation, therefore
1319 * we first add 1.5 then subtract 1 again
1321 tx
= static_cast<int>(dx
*b_inv
[XX
] + h15
);
1322 ty
= static_cast<int>(dy
*b_inv
[YY
] + h15
);
1323 tz
= static_cast<int>(dz
*b_inv
[ZZ
] + h15
);
1328 /* Correct diff vector for translation */
1329 ddx
= tx
*box_size
[XX
] - dx
;
1330 ddy
= ty
*box_size
[YY
] - dy
;
1331 ddz
= tz
*box_size
[ZZ
] - dz
;
1333 /* Distance squared */
1334 r2
= (ddx
*ddx
) + (ddy
*ddy
) + (ddz
*ddz
);
1336 *shift
= XYZ2IS(tx
, ty
, tz
);
1341 static void add_simple(t_ns_buf
* nsbuf
, int nrj
, int cg_j
,
1342 gmx_bool bHaveVdW
[], int ngid
, t_mdatoms
*md
,
1343 int icg
, int jgid
, t_block
*cgs
, t_excl bexcl
[],
1344 int shift
, t_forcerec
*fr
, put_in_list_t
*put_in_list
)
1346 if (nsbuf
->nj
+ nrj
> MAX_CG
)
1348 put_in_list(bHaveVdW
, ngid
, md
, icg
, jgid
, nsbuf
->ncg
, nsbuf
->jcg
,
1349 cgs
->index
, bexcl
, shift
, fr
, TRUE
, TRUE
, fr
->solvent_opt
);
1350 /* Reset buffer contents */
1351 nsbuf
->ncg
= nsbuf
->nj
= 0;
1353 nsbuf
->jcg
[nsbuf
->ncg
++] = cg_j
;
1357 static void ns_inner_tric(rvec x
[], int icg
, int *i_egp_flags
,
1358 int njcg
, int jcg
[],
1359 matrix box
, rvec b_inv
, real rcut2
,
1360 t_block
*cgs
, t_ns_buf
**ns_buf
,
1361 gmx_bool bHaveVdW
[], int ngid
, t_mdatoms
*md
,
1362 t_excl bexcl
[], t_forcerec
*fr
,
1363 put_in_list_t
*put_in_list
)
1367 int *cginfo
= fr
->cginfo
;
1370 cgindex
= cgs
->index
;
1372 for (j
= 0; (j
< njcg
); j
++)
1375 nrj
= cgindex
[cg_j
+1]-cgindex
[cg_j
];
1376 if (calc_image_tric(x
[icg
], x
[cg_j
], box
, b_inv
, &shift
) < rcut2
)
1378 jgid
= GET_CGINFO_GID(cginfo
[cg_j
]);
1379 if (!(i_egp_flags
[jgid
] & EGP_EXCL
))
1381 add_simple(&ns_buf
[jgid
][shift
], nrj
, cg_j
,
1382 bHaveVdW
, ngid
, md
, icg
, jgid
, cgs
, bexcl
, shift
, fr
,
1389 static void ns_inner_rect(rvec x
[], int icg
, int *i_egp_flags
,
1390 int njcg
, int jcg
[],
1391 gmx_bool bBox
, rvec box_size
, rvec b_inv
, real rcut2
,
1392 t_block
*cgs
, t_ns_buf
**ns_buf
,
1393 gmx_bool bHaveVdW
[], int ngid
, t_mdatoms
*md
,
1394 t_excl bexcl
[], t_forcerec
*fr
,
1395 put_in_list_t
*put_in_list
)
1399 int *cginfo
= fr
->cginfo
;
1402 cgindex
= cgs
->index
;
1406 for (j
= 0; (j
< njcg
); j
++)
1409 nrj
= cgindex
[cg_j
+1]-cgindex
[cg_j
];
1410 if (calc_image_rect(x
[icg
], x
[cg_j
], box_size
, b_inv
, &shift
) < rcut2
)
1412 jgid
= GET_CGINFO_GID(cginfo
[cg_j
]);
1413 if (!(i_egp_flags
[jgid
] & EGP_EXCL
))
1415 add_simple(&ns_buf
[jgid
][shift
], nrj
, cg_j
,
1416 bHaveVdW
, ngid
, md
, icg
, jgid
, cgs
, bexcl
, shift
, fr
,
1424 for (j
= 0; (j
< njcg
); j
++)
1427 nrj
= cgindex
[cg_j
+1]-cgindex
[cg_j
];
1428 if ((rcut2
== 0) || (distance2(x
[icg
], x
[cg_j
]) < rcut2
))
1430 jgid
= GET_CGINFO_GID(cginfo
[cg_j
]);
1431 if (!(i_egp_flags
[jgid
] & EGP_EXCL
))
1433 add_simple(&ns_buf
[jgid
][CENTRAL
], nrj
, cg_j
,
1434 bHaveVdW
, ngid
, md
, icg
, jgid
, cgs
, bexcl
, CENTRAL
, fr
,
1442 /* ns_simple_core needs to be adapted for QMMM still 2005 */
1444 static int ns_simple_core(t_forcerec
*fr
,
1445 gmx_localtop_t
*top
,
1447 matrix box
, rvec box_size
,
1448 t_excl bexcl
[], int *aaj
,
1449 int ngid
, t_ns_buf
**ns_buf
,
1450 put_in_list_t
*put_in_list
, gmx_bool bHaveVdW
[])
1454 int nsearch
, icg
, igid
, nn
;
1458 t_block
*cgs
= &(top
->cgs
);
1459 t_blocka
*excl
= &(top
->excls
);
1462 gmx_bool bBox
, bTriclinic
;
1465 rlist2
= gmx::square(fr
->rlist
);
1467 bBox
= (fr
->ePBC
!= epbcNONE
);
1470 for (m
= 0; (m
< DIM
); m
++)
1472 if (gmx_numzero(box_size
[m
]))
1474 gmx_fatal(FARGS
, "Dividing by zero box size!");
1476 b_inv
[m
] = 1.0/box_size
[m
];
1478 bTriclinic
= TRICLINIC(box
);
1485 cginfo
= fr
->cginfo
;
1488 for (icg
= fr
->cg0
; (icg
< fr
->hcg
); icg
++)
1491 i0 = cgs->index[icg];
1492 nri = cgs->index[icg+1]-i0;
1493 i_atoms = &(cgs->a[i0]);
1494 i_eg_excl = fr->eg_excl + ngid*md->cENER[*i_atoms];
1495 setexcl(nri,i_atoms,excl,TRUE,bexcl);
1497 igid
= GET_CGINFO_GID(cginfo
[icg
]);
1498 i_egp_flags
= fr
->egp_flags
+ ngid
*igid
;
1499 setexcl(cgs
->index
[icg
], cgs
->index
[icg
+1], excl
, TRUE
, bexcl
);
1501 naaj
= calc_naaj(icg
, cgs
->nr
);
1504 ns_inner_tric(fr
->cg_cm
, icg
, i_egp_flags
, naaj
, &(aaj
[icg
]),
1505 box
, b_inv
, rlist2
, cgs
, ns_buf
,
1506 bHaveVdW
, ngid
, md
, bexcl
, fr
, put_in_list
);
1510 ns_inner_rect(fr
->cg_cm
, icg
, i_egp_flags
, naaj
, &(aaj
[icg
]),
1511 bBox
, box_size
, b_inv
, rlist2
, cgs
, ns_buf
,
1512 bHaveVdW
, ngid
, md
, bexcl
, fr
, put_in_list
);
1516 for (nn
= 0; (nn
< ngid
); nn
++)
1518 for (k
= 0; (k
< SHIFTS
); k
++)
1520 nsbuf
= &(ns_buf
[nn
][k
]);
1523 put_in_list(bHaveVdW
, ngid
, md
, icg
, nn
, nsbuf
->ncg
, nsbuf
->jcg
,
1524 cgs
->index
, bexcl
, k
, fr
, TRUE
, TRUE
, fr
->solvent_opt
);
1525 nsbuf
->ncg
= nsbuf
->nj
= 0;
1529 /* setexcl(nri,i_atoms,excl,FALSE,bexcl); */
1530 setexcl(cgs
->index
[icg
], cgs
->index
[icg
+1], excl
, FALSE
, bexcl
);
1532 close_neighbor_lists(fr
, FALSE
);
1537 /************************************************
1539 * N S 5 G R I D S T U F F
1541 ************************************************/
1543 static inline void get_dx_dd(int Nx
, real gridx
, real rc2
, int xgi
, real x
,
1544 int ncpddc
, int shift_min
, int shift_max
,
1545 int *g0
, int *g1
, real
*dcx2
)
1548 int g_min
, g_max
, shift_home
;
1581 g_min
= (shift_min
== shift_home
? 0 : ncpddc
);
1582 g_max
= (shift_max
== shift_home
? ncpddc
- 1 : Nx
- 1);
1589 else if (shift_max
< 0)
1604 /* Check one grid cell down */
1605 dcx
= ((*g0
- 1) + 1)*gridx
- x
;
1617 /* Check one grid cell up */
1618 dcx
= (*g1
+ 1)*gridx
- x
;
1630 #define calc_dx2(XI, YI, ZI, y) (gmx::square(XI-y[XX]) + gmx::square(YI-y[YY]) + gmx::square(ZI-y[ZZ]))
1631 #define calc_cyl_dx2(XI, YI, y) (gmx::square(XI-y[XX]) + gmx::square(YI-y[YY]))
1632 /****************************************************
1634 * F A S T N E I G H B O R S E A R C H I N G
1636 * Optimized neighboursearching routine using grid
1637 * at least 1x1x1, see GROMACS manual
1639 ****************************************************/
1642 static void get_cutoff2(t_forcerec
*fr
, real
*rs2
)
1644 *rs2
= gmx::square(fr
->rlist
);
1647 static void init_nsgrid_lists(t_forcerec
*fr
, int ngid
, gmx_ns_t
*ns
)
1652 get_cutoff2(fr
, &rs2
);
1654 /* Short range buffers */
1655 snew(ns
->nl_sr
, ngid
);
1657 snew(ns
->nsr
, ngid
);
1659 for (j
= 0; (j
< ngid
); j
++)
1661 snew(ns
->nl_sr
[j
], MAX_CG
);
1666 "ns5_core: rs2 = %g (nm^2)\n",
1671 static int nsgrid_core(const t_commrec
*cr
, t_forcerec
*fr
,
1672 matrix box
, int ngid
,
1673 gmx_localtop_t
*top
,
1675 t_excl bexcl
[], gmx_bool
*bExcludeAlleg
,
1677 put_in_list_t
*put_in_list
,
1678 gmx_bool bHaveVdW
[],
1679 gmx_bool bMakeQMMMnblist
)
1685 t_block
*cgs
= &(top
->cgs
);
1686 int *cginfo
= fr
->cginfo
;
1687 /* int *i_atoms,*cgsindex=cgs->index; */
1689 int cell_x
, cell_y
, cell_z
;
1690 int d
, tx
, ty
, tz
, dx
, dy
, dz
, cj
;
1691 #ifdef ALLOW_OFFDIAG_LT_HALFDIAG
1692 int zsh_ty
, zsh_tx
, ysh_tx
;
1694 int dx0
, dx1
, dy0
, dy1
, dz0
, dz1
;
1695 int Nx
, Ny
, Nz
, shift
= -1, j
, nrj
, nns
, nn
= -1;
1696 real gridx
, gridy
, gridz
, grid_x
, grid_y
;
1697 real
*dcx2
, *dcy2
, *dcz2
;
1699 int cg0
, cg1
, icg
= -1, cgsnr
, i0
, igid
, naaj
, max_jcg
;
1700 int jcg0
, jcg1
, jjcg
, cgj0
, jgid
;
1701 int *grida
, *gridnra
, *gridind
;
1702 rvec
*cgcm
, grid_offset
;
1703 real r2
, rs2
, XI
, YI
, ZI
, tmp1
, tmp2
;
1705 gmx_bool bDomDec
, bTriclinicX
, bTriclinicY
;
1710 bDomDec
= DOMAINDECOMP(cr
);
1713 bTriclinicX
= ((YY
< grid
->npbcdim
&&
1714 (!bDomDec
|| dd
->nc
[YY
] == 1) && box
[YY
][XX
] != 0) ||
1715 (ZZ
< grid
->npbcdim
&&
1716 (!bDomDec
|| dd
->nc
[ZZ
] == 1) && box
[ZZ
][XX
] != 0));
1717 bTriclinicY
= (ZZ
< grid
->npbcdim
&&
1718 (!bDomDec
|| dd
->nc
[ZZ
] == 1) && box
[ZZ
][YY
] != 0);
1722 get_cutoff2(fr
, &rs2
);
1733 gridind
= grid
->index
;
1734 gridnra
= grid
->nra
;
1737 gridx
= grid
->cell_size
[XX
];
1738 gridy
= grid
->cell_size
[YY
];
1739 gridz
= grid
->cell_size
[ZZ
];
1742 copy_rvec(grid
->cell_offset
, grid_offset
);
1743 copy_ivec(grid
->ncpddc
, ncpddc
);
1748 #ifdef ALLOW_OFFDIAG_LT_HALFDIAG
1749 zsh_ty
= floor(-box
[ZZ
][YY
]/box
[YY
][YY
]+0.5);
1750 zsh_tx
= floor(-box
[ZZ
][XX
]/box
[XX
][XX
]+0.5);
1751 ysh_tx
= floor(-box
[YY
][XX
]/box
[XX
][XX
]+0.5);
1752 if (zsh_tx
!= 0 && ysh_tx
!= 0)
1754 /* This could happen due to rounding, when both ratios are 0.5 */
1761 /* We only want a list for the test particle */
1770 /* Set the shift range */
1771 for (d
= 0; d
< DIM
; d
++)
1775 /* Check if we need periodicity shifts.
1776 * Without PBC or with domain decomposition we don't need them.
1778 if (d
>= ePBC2npbcdim(fr
->ePBC
) || (bDomDec
&& dd
->nc
[d
] > 1))
1785 box
[XX
][XX
] - std::abs(box
[YY
][XX
]) - std::abs(box
[ZZ
][XX
]) < std::sqrt(rs2
))
1796 /* Loop over charge groups */
1797 for (icg
= cg0
; (icg
< cg1
); icg
++)
1799 igid
= GET_CGINFO_GID(cginfo
[icg
]);
1800 /* Skip this charge group if all energy groups are excluded! */
1801 if (bExcludeAlleg
[igid
])
1806 i0
= cgs
->index
[icg
];
1808 if (bMakeQMMMnblist
)
1810 /* Skip this charge group if it is not a QM atom while making a
1811 * QM/MM neighbourlist
1813 if (md
->bQM
[i0
] == FALSE
)
1815 continue; /* MM particle, go to next particle */
1818 /* Compute the number of charge groups that fall within the control
1821 naaj
= calc_naaj(icg
, cgsnr
);
1828 /* make a normal neighbourlist */
1832 /* Get the j charge-group and dd cell shift ranges */
1833 dd_get_ns_ranges(cr
->dd
, icg
, &jcg0
, &jcg1
, sh0
, sh1
);
1838 /* Compute the number of charge groups that fall within the control
1841 naaj
= calc_naaj(icg
, cgsnr
);
1847 /* The i-particle is awlways the test particle,
1848 * so we want all j-particles
1850 max_jcg
= cgsnr
- 1;
1854 max_jcg
= jcg1
- cgsnr
;
1859 i_egp_flags
= fr
->egp_flags
+ igid
*ngid
;
1861 /* Set the exclusions for the atoms in charge group icg using a bitmask */
1862 setexcl(i0
, cgs
->index
[icg
+1], &top
->excls
, TRUE
, bexcl
);
1864 ci2xyz(grid
, icg
, &cell_x
, &cell_y
, &cell_z
);
1866 /* Changed iicg to icg, DvdS 990115
1867 * (but see consistency check above, DvdS 990330)
1870 fprintf(log
, "icg=%5d, naaj=%5d, cell %d %d %d\n",
1871 icg
, naaj
, cell_x
, cell_y
, cell_z
);
1873 /* Loop over shift vectors in three dimensions */
1874 for (tz
= -shp
[ZZ
]; tz
<= shp
[ZZ
]; tz
++)
1876 ZI
= cgcm
[icg
][ZZ
]+tz
*box
[ZZ
][ZZ
];
1877 /* Calculate range of cells in Z direction that have the shift tz */
1878 zgi
= cell_z
+ tz
*Nz
;
1879 get_dx_dd(Nz
, gridz
, rs2
, zgi
, ZI
-grid_offset
[ZZ
],
1880 ncpddc
[ZZ
], sh0
[ZZ
], sh1
[ZZ
], &dz0
, &dz1
, dcz2
);
1885 for (ty
= -shp
[YY
]; ty
<= shp
[YY
]; ty
++)
1887 YI
= cgcm
[icg
][YY
]+ty
*box
[YY
][YY
]+tz
*box
[ZZ
][YY
];
1888 /* Calculate range of cells in Y direction that have the shift ty */
1891 ygi
= (int)(Ny
+ (YI
- grid_offset
[YY
])*grid_y
) - Ny
;
1895 ygi
= cell_y
+ ty
*Ny
;
1897 get_dx_dd(Ny
, gridy
, rs2
, ygi
, YI
-grid_offset
[YY
],
1898 ncpddc
[YY
], sh0
[YY
], sh1
[YY
], &dy0
, &dy1
, dcy2
);
1903 for (tx
= -shp
[XX
]; tx
<= shp
[XX
]; tx
++)
1905 XI
= cgcm
[icg
][XX
]+tx
*box
[XX
][XX
]+ty
*box
[YY
][XX
]+tz
*box
[ZZ
][XX
];
1906 /* Calculate range of cells in X direction that have the shift tx */
1909 xgi
= (int)(Nx
+ (XI
- grid_offset
[XX
])*grid_x
) - Nx
;
1913 xgi
= cell_x
+ tx
*Nx
;
1915 get_dx_dd(Nx
, gridx
, rs2
, xgi
, XI
-grid_offset
[XX
],
1916 ncpddc
[XX
], sh0
[XX
], sh1
[XX
], &dx0
, &dx1
, dcx2
);
1921 /* Get shift vector */
1922 shift
= XYZ2IS(tx
, ty
, tz
);
1924 range_check(shift
, 0, SHIFTS
);
1926 for (nn
= 0; (nn
< ngid
); nn
++)
1931 fprintf(log
, "shift: %2d, dx0,1: %2d,%2d, dy0,1: %2d,%2d, dz0,1: %2d,%2d\n",
1932 shift
, dx0
, dx1
, dy0
, dy1
, dz0
, dz1
);
1933 fprintf(log
, "cgcm: %8.3f %8.3f %8.3f\n", cgcm
[icg
][XX
],
1934 cgcm
[icg
][YY
], cgcm
[icg
][ZZ
]);
1935 fprintf(log
, "xi: %8.3f %8.3f %8.3f\n", XI
, YI
, ZI
);
1937 for (dx
= dx0
; (dx
<= dx1
); dx
++)
1939 tmp1
= rs2
- dcx2
[dx
];
1940 for (dy
= dy0
; (dy
<= dy1
); dy
++)
1942 tmp2
= tmp1
- dcy2
[dy
];
1945 for (dz
= dz0
; (dz
<= dz1
); dz
++)
1947 if (tmp2
> dcz2
[dz
])
1949 /* Find grid-cell cj in which possible neighbours are */
1950 cj
= xyz2ci(Ny
, Nz
, dx
, dy
, dz
);
1952 /* Check out how many cgs (nrj) there in this cell */
1955 /* Find the offset in the cg list */
1958 /* Check if all j's are out of range so we
1959 * can skip the whole cell.
1960 * Should save some time, especially with DD.
1963 (grida
[cgj0
] >= max_jcg
&&
1964 (grida
[cgj0
] >= jcg1
|| grida
[cgj0
+nrj
-1] < jcg0
)))
1970 for (j
= 0; (j
< nrj
); j
++)
1972 jjcg
= grida
[cgj0
+j
];
1974 /* check whether this guy is in range! */
1975 if ((jjcg
>= jcg0
&& jjcg
< jcg1
) ||
1978 r2
= calc_dx2(XI
, YI
, ZI
, cgcm
[jjcg
]);
1981 /* jgid = gid[cgsatoms[cgsindex[jjcg]]]; */
1982 jgid
= GET_CGINFO_GID(cginfo
[jjcg
]);
1983 /* check energy group exclusions */
1984 if (!(i_egp_flags
[jgid
] & EGP_EXCL
))
1986 if (nsr
[jgid
] >= MAX_CG
)
1988 /* Add to short-range list */
1989 put_in_list(bHaveVdW
, ngid
, md
, icg
, jgid
,
1990 nsr
[jgid
], nl_sr
[jgid
],
1991 cgs
->index
, /* cgsatoms, */ bexcl
,
1992 shift
, fr
, TRUE
, TRUE
, fr
->solvent_opt
);
1995 nl_sr
[jgid
][nsr
[jgid
]++] = jjcg
;
2006 /* CHECK whether there is anything left in the buffers */
2007 for (nn
= 0; (nn
< ngid
); nn
++)
2011 put_in_list(bHaveVdW
, ngid
, md
, icg
, nn
, nsr
[nn
], nl_sr
[nn
],
2012 cgs
->index
, /* cgsatoms, */ bexcl
,
2013 shift
, fr
, TRUE
, TRUE
, fr
->solvent_opt
);
2019 setexcl(cgs
->index
[icg
], cgs
->index
[icg
+1], &top
->excls
, FALSE
, bexcl
);
2021 /* No need to perform any left-over force calculations anymore (as we used to do here)
2022 * since we now save the proper long-range lists for later evaluation.
2025 /* Close neighbourlists */
2026 close_neighbor_lists(fr
, bMakeQMMMnblist
);
2031 static void ns_realloc_natoms(gmx_ns_t
*ns
, int natoms
)
2035 if (natoms
> ns
->nra_alloc
)
2037 ns
->nra_alloc
= over_alloc_dd(natoms
);
2038 srenew(ns
->bexcl
, ns
->nra_alloc
);
2039 for (i
= 0; i
< ns
->nra_alloc
; i
++)
2046 void init_ns(FILE *fplog
, const t_commrec
*cr
,
2047 gmx_ns_t
*ns
, t_forcerec
*fr
,
2048 const gmx_mtop_t
*mtop
)
2050 int icg
, nr_in_cg
, maxcg
, i
, j
, jcg
, ngid
, ncg
;
2052 /* Compute largest charge groups size (# atoms) */
2054 for (const gmx_moltype_t
&molt
: mtop
->moltype
)
2056 const t_block
*cgs
= &molt
.cgs
;
2057 for (icg
= 0; (icg
< cgs
->nr
); icg
++)
2059 nr_in_cg
= std::max(nr_in_cg
, (int)(cgs
->index
[icg
+1]-cgs
->index
[icg
]));
2063 /* Verify whether largest charge group is <= max cg.
2064 * This is determined by the type of the local exclusion type
2065 * Exclusions are stored in bits. (If the type is not large
2066 * enough, enlarge it, unsigned char -> unsigned short -> unsigned long)
2068 maxcg
= sizeof(t_excl
)*8;
2069 if (nr_in_cg
> maxcg
)
2071 gmx_fatal(FARGS
, "Max #atoms in a charge group: %d > %d\n",
2075 ngid
= mtop
->groups
.grps
[egcENER
].nr
;
2076 snew(ns
->bExcludeAlleg
, ngid
);
2077 for (i
= 0; i
< ngid
; i
++)
2079 ns
->bExcludeAlleg
[i
] = TRUE
;
2080 for (j
= 0; j
< ngid
; j
++)
2082 if (!(fr
->egp_flags
[i
*ngid
+j
] & EGP_EXCL
))
2084 ns
->bExcludeAlleg
[i
] = FALSE
;
2092 ns
->grid
= init_grid(fplog
, fr
);
2093 init_nsgrid_lists(fr
, ngid
, ns
);
2098 snew(ns
->ns_buf
, ngid
);
2099 for (i
= 0; (i
< ngid
); i
++)
2101 snew(ns
->ns_buf
[i
], SHIFTS
);
2103 ncg
= ncg_mtop(mtop
);
2104 snew(ns
->simple_aaj
, 2*ncg
);
2105 for (jcg
= 0; (jcg
< ncg
); jcg
++)
2107 ns
->simple_aaj
[jcg
] = jcg
;
2108 ns
->simple_aaj
[jcg
+ncg
] = jcg
;
2112 /* Create array that determines whether or not atoms have VdW */
2113 snew(ns
->bHaveVdW
, fr
->ntype
);
2114 for (i
= 0; (i
< fr
->ntype
); i
++)
2116 for (j
= 0; (j
< fr
->ntype
); j
++)
2118 ns
->bHaveVdW
[i
] = (ns
->bHaveVdW
[i
] ||
2120 ((BHAMA(fr
->nbfp
, fr
->ntype
, i
, j
) != 0) ||
2121 (BHAMB(fr
->nbfp
, fr
->ntype
, i
, j
) != 0) ||
2122 (BHAMC(fr
->nbfp
, fr
->ntype
, i
, j
) != 0)) :
2123 ((C6(fr
->nbfp
, fr
->ntype
, i
, j
) != 0) ||
2124 (C12(fr
->nbfp
, fr
->ntype
, i
, j
) != 0))));
2129 pr_bvec(debug
, 0, "bHaveVdW", ns
->bHaveVdW
, fr
->ntype
, TRUE
);
2133 ns
->bexcl
= nullptr;
2134 if (!DOMAINDECOMP(cr
))
2136 ns_realloc_natoms(ns
, mtop
->natoms
);
2139 ns
->nblist_initialized
= FALSE
;
2141 /* nbr list debug dump */
2143 char *ptr
= getenv("GMX_DUMP_NL");
2146 ns
->dump_nl
= strtol(ptr
, nullptr, 10);
2149 fprintf(fplog
, "GMX_DUMP_NL = %d", ns
->dump_nl
);
2159 void done_ns(gmx_ns_t
*ns
, int numEnergyGroups
)
2161 sfree(ns
->bExcludeAlleg
);
2164 for (int i
= 0; i
< numEnergyGroups
; i
++)
2166 sfree(ns
->ns_buf
[i
]);
2170 sfree(ns
->simple_aaj
);
2171 sfree(ns
->bHaveVdW
);
2172 done_grid(ns
->grid
);
2176 int search_neighbours(FILE *log
, t_forcerec
*fr
,
2178 gmx_localtop_t
*top
,
2179 gmx_groups_t
*groups
,
2180 const t_commrec
*cr
,
2181 t_nrnb
*nrnb
, t_mdatoms
*md
,
2184 t_block
*cgs
= &(top
->cgs
);
2185 rvec box_size
, grid_x0
, grid_x1
;
2187 real min_size
, grid_dens
;
2193 gmx_domdec_zones_t
*dd_zones
;
2194 put_in_list_t
*put_in_list
;
2198 /* Set some local variables */
2200 ngid
= groups
->grps
[egcENER
].nr
;
2202 for (m
= 0; (m
< DIM
); m
++)
2204 box_size
[m
] = box
[m
][m
];
2207 if (fr
->ePBC
!= epbcNONE
)
2209 if (gmx::square(fr
->rlist
) >= max_cutoff2(fr
->ePBC
, box
))
2211 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.");
2215 min_size
= std::min(box_size
[XX
], std::min(box_size
[YY
], box_size
[ZZ
]));
2216 if (2*fr
->rlist
>= min_size
)
2218 gmx_fatal(FARGS
, "One of the box diagonal elements has become smaller than twice the cut-off length.");
2223 if (DOMAINDECOMP(cr
))
2225 ns_realloc_natoms(ns
, cgs
->index
[cgs
->nr
]);
2228 /* Reset the neighbourlists */
2229 reset_neighbor_lists(fr
);
2231 if (bGrid
&& bFillGrid
)
2235 if (DOMAINDECOMP(cr
))
2237 dd_zones
= domdec_zones(cr
->dd
);
2243 get_nsgrid_boundaries(grid
->nboundeddim
, box
, nullptr, nullptr, nullptr, nullptr,
2244 cgs
->nr
, fr
->cg_cm
, grid_x0
, grid_x1
, &grid_dens
);
2246 grid_first(log
, grid
, nullptr, nullptr, box
, grid_x0
, grid_x1
,
2247 fr
->rlist
, grid_dens
);
2253 if (DOMAINDECOMP(cr
))
2256 fill_grid(dd_zones
, grid
, end
, -1, end
, fr
->cg_cm
);
2258 grid
->icg1
= dd_zones
->izone
[dd_zones
->nizone
-1].cg1
;
2262 fill_grid(nullptr, grid
, cgs
->nr
, fr
->cg0
, fr
->hcg
, fr
->cg_cm
);
2263 grid
->icg0
= fr
->cg0
;
2264 grid
->icg1
= fr
->hcg
;
2267 calc_elemnr(grid
, start
, end
, cgs
->nr
);
2269 grid_last(grid
, start
, end
, cgs
->nr
);
2274 print_grid(debug
, grid
);
2279 /* Set the grid cell index for the test particle only.
2280 * The cell to cg index is not corrected, but that does not matter.
2282 fill_grid(nullptr, ns
->grid
, fr
->hcg
, fr
->hcg
-1, fr
->hcg
, fr
->cg_cm
);
2285 if (!fr
->ns
->bCGlist
)
2287 put_in_list
= put_in_list_at
;
2291 put_in_list
= put_in_list_cg
;
2298 nsearch
= nsgrid_core(cr
, fr
, box
, ngid
, top
,
2299 grid
, ns
->bexcl
, ns
->bExcludeAlleg
,
2300 md
, put_in_list
, ns
->bHaveVdW
,
2303 /* neighbour searching withouth QMMM! QM atoms have zero charge in
2304 * the classical calculation. The charge-charge interaction
2305 * between QM and MM atoms is handled in the QMMM core calculation
2306 * (see QMMM.c). The VDW however, we'd like to compute classically
2307 * and the QM MM atom pairs have just been put in the
2308 * corresponding neighbourlists. in case of QMMM we still need to
2309 * fill a special QMMM neighbourlist that contains all neighbours
2310 * of the QM atoms. If bQMMM is true, this list will now be made:
2312 if (fr
->bQMMM
&& fr
->qr
->QMMMscheme
!= eQMMMschemeoniom
)
2314 nsearch
+= nsgrid_core(cr
, fr
, box
, ngid
, top
,
2315 grid
, ns
->bexcl
, ns
->bExcludeAlleg
,
2316 md
, put_in_list_qmmm
, ns
->bHaveVdW
,
2322 nsearch
= ns_simple_core(fr
, top
, md
, box
, box_size
,
2323 ns
->bexcl
, ns
->simple_aaj
,
2324 ngid
, ns
->ns_buf
, put_in_list
, ns
->bHaveVdW
);
2327 inc_nrnb(nrnb
, eNR_NS
, nsearch
);