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, 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.
44 #include "gromacs/utility/smalloc.h"
46 #include "gromacs/math/utilities.h"
48 #include "types/commrec.h"
52 #include "nonbonded.h"
56 #include "gmx_fatal.h"
59 #include "mtop_util.h"
66 * E X C L U S I O N H A N D L I N G
70 static void SETEXCL_(t_excl e
[], atom_id i
, atom_id j
)
74 static void RMEXCL_(t_excl e
[], atom_id i
, atom_id j
)
76 e
[j
] = e
[j
] & ~(1<<i
);
78 static gmx_bool
ISEXCL_(t_excl e
[], atom_id i
, atom_id j
)
80 return (gmx_bool
)(e
[j
] & (1<<i
));
82 static gmx_bool
NOTEXCL_(t_excl e
[], atom_id i
, atom_id j
)
84 return !(ISEXCL(e
, i
, j
));
87 #define SETEXCL(e, i, j) (e)[((atom_id) (j))] |= (1<<((atom_id) (i)))
88 #define RMEXCL(e, i, j) (e)[((atom_id) (j))] &= (~(1<<((atom_id) (i))))
89 #define ISEXCL(e, i, j) (gmx_bool) ((e)[((atom_id) (j))] & (1<<((atom_id) (i))))
90 #define NOTEXCL(e, i, j) !(ISEXCL(e, i, j))
94 round_up_to_simd_width(int length
, int simd_width
)
96 int offset
, newlength
;
98 offset
= (simd_width
> 0) ? length
% simd_width
: 0;
100 return (offset
== 0) ? length
: length
-offset
+simd_width
;
102 /************************************************
104 * U T I L I T I E S F O R N S
106 ************************************************/
108 void reallocate_nblist(t_nblist
*nl
)
112 fprintf(debug
, "reallocating neigborlist (ielec=%d, ivdw=%d, igeometry=%d, type=%d), maxnri=%d\n",
113 nl
->ielec
, nl
->ivdw
, nl
->igeometry
, nl
->type
, nl
->maxnri
);
115 srenew(nl
->iinr
, nl
->maxnri
);
116 if (nl
->igeometry
== GMX_NBLIST_GEOMETRY_CG_CG
)
118 srenew(nl
->iinr_end
, nl
->maxnri
);
120 srenew(nl
->gid
, nl
->maxnri
);
121 srenew(nl
->shift
, nl
->maxnri
);
122 srenew(nl
->jindex
, nl
->maxnri
+1);
126 static void init_nblist(FILE *log
, t_nblist
*nl_sr
, t_nblist
*nl_lr
,
127 int maxsr
, int maxlr
,
128 int ivdw
, int ivdwmod
,
129 int ielec
, int ielecmod
,
130 int igeometry
, int type
,
131 gmx_bool bElecAndVdwSwitchDiffers
)
137 for (i
= 0; (i
< 2); i
++)
139 nl
= (i
== 0) ? nl_sr
: nl_lr
;
140 homenr
= (i
== 0) ? maxsr
: maxlr
;
148 /* Set coul/vdw in neighborlist, and for the normal loops we determine
149 * an index of which one to call.
152 nl
->ivdwmod
= ivdwmod
;
154 nl
->ielecmod
= ielecmod
;
156 nl
->igeometry
= igeometry
;
158 if (nl
->type
== GMX_NBLIST_INTERACTION_FREE_ENERGY
)
160 nl
->igeometry
= GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE
;
163 /* This will also set the simd_padding_width field */
164 gmx_nonbonded_set_kernel_pointers( (i
== 0) ? log
: NULL
, nl
, bElecAndVdwSwitchDiffers
);
166 /* maxnri is influenced by the number of shifts (maximum is 8)
167 * and the number of energy groups.
168 * If it is not enough, nl memory will be reallocated during the run.
169 * 4 seems to be a reasonable factor, which only causes reallocation
170 * during runs with tiny and many energygroups.
172 nl
->maxnri
= homenr
*4;
182 reallocate_nblist(nl
);
187 fprintf(debug
, "Initiating neighbourlist (ielec=%d, ivdw=%d, type=%d) for %s interactions,\nwith %d SR, %d LR atoms.\n",
188 nl
->ielec
, nl
->ivdw
, nl
->type
, gmx_nblist_geometry_names
[nl
->igeometry
], maxsr
, maxlr
);
193 void init_neighbor_list(FILE *log
, t_forcerec
*fr
, int homenr
)
195 /* Make maxlr tunable! (does not seem to be a big difference though)
196 * This parameter determines the number of i particles in a long range
197 * neighbourlist. Too few means many function calls, too many means
200 int maxsr
, maxsr_wat
, maxlr
, maxlr_wat
;
201 int ielec
, ivdw
, ielecmod
, ivdwmod
, type
;
203 int igeometry_def
, igeometry_w
, igeometry_ww
;
205 gmx_bool bElecAndVdwSwitchDiffers
;
208 /* maxsr = homenr-fr->nWatMol*3; */
213 gmx_fatal(FARGS
, "%s, %d: Negative number of short range atoms.\n"
214 "Call your Gromacs dealer for assistance.", __FILE__
, __LINE__
);
216 /* This is just for initial allocation, so we do not reallocate
217 * all the nlist arrays many times in a row.
218 * The numbers seem very accurate, but they are uncritical.
220 maxsr_wat
= min(fr
->nWatMol
, (homenr
+2)/3);
224 maxlr_wat
= min(maxsr_wat
, maxlr
);
228 maxlr
= maxlr_wat
= 0;
231 /* Determine the values for ielec/ivdw. */
232 ielec
= fr
->nbkernel_elec_interaction
;
233 ivdw
= fr
->nbkernel_vdw_interaction
;
234 ielecmod
= fr
->nbkernel_elec_modifier
;
235 ivdwmod
= fr
->nbkernel_vdw_modifier
;
236 type
= GMX_NBLIST_INTERACTION_STANDARD
;
237 bElecAndVdwSwitchDiffers
= ( (fr
->rcoulomb_switch
!= fr
->rvdw_switch
) || (fr
->rcoulomb
!= fr
->rvdw
));
239 fr
->ns
.bCGlist
= (getenv("GMX_NBLISTCG") != 0);
242 igeometry_def
= GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE
;
246 igeometry_def
= GMX_NBLIST_GEOMETRY_CG_CG
;
249 fprintf(log
, "\nUsing charge-group - charge-group neighbor lists and kernels\n\n");
253 if (fr
->solvent_opt
== esolTIP4P
)
255 igeometry_w
= GMX_NBLIST_GEOMETRY_WATER4_PARTICLE
;
256 igeometry_ww
= GMX_NBLIST_GEOMETRY_WATER4_WATER4
;
260 igeometry_w
= GMX_NBLIST_GEOMETRY_WATER3_PARTICLE
;
261 igeometry_ww
= GMX_NBLIST_GEOMETRY_WATER3_WATER3
;
264 for (i
= 0; i
< fr
->nnblists
; i
++)
266 nbl
= &(fr
->nblists
[i
]);
268 if ((fr
->adress_type
!= eAdressOff
) && (i
>= fr
->nnblists
/2))
270 type
= GMX_NBLIST_INTERACTION_ADRESS
;
272 init_nblist(log
, &nbl
->nlist_sr
[eNL_VDWQQ
], &nbl
->nlist_lr
[eNL_VDWQQ
],
273 maxsr
, maxlr
, ivdw
, ivdwmod
, ielec
, ielecmod
, igeometry_def
, type
, bElecAndVdwSwitchDiffers
);
274 init_nblist(log
, &nbl
->nlist_sr
[eNL_VDW
], &nbl
->nlist_lr
[eNL_VDW
],
275 maxsr
, maxlr
, ivdw
, ivdwmod
, GMX_NBKERNEL_ELEC_NONE
, eintmodNONE
, igeometry_def
, type
, bElecAndVdwSwitchDiffers
);
276 init_nblist(log
, &nbl
->nlist_sr
[eNL_QQ
], &nbl
->nlist_lr
[eNL_QQ
],
277 maxsr
, maxlr
, GMX_NBKERNEL_VDW_NONE
, eintmodNONE
, ielec
, ielecmod
, igeometry_def
, type
, bElecAndVdwSwitchDiffers
);
278 init_nblist(log
, &nbl
->nlist_sr
[eNL_VDWQQ_WATER
], &nbl
->nlist_lr
[eNL_VDWQQ_WATER
],
279 maxsr_wat
, maxlr_wat
, ivdw
, ivdwmod
, ielec
, ielecmod
, igeometry_w
, type
, bElecAndVdwSwitchDiffers
);
280 init_nblist(log
, &nbl
->nlist_sr
[eNL_QQ_WATER
], &nbl
->nlist_lr
[eNL_QQ_WATER
],
281 maxsr_wat
, maxlr_wat
, GMX_NBKERNEL_VDW_NONE
, eintmodNONE
, ielec
, ielecmod
, igeometry_w
, type
, bElecAndVdwSwitchDiffers
);
282 init_nblist(log
, &nbl
->nlist_sr
[eNL_VDWQQ_WATERWATER
], &nbl
->nlist_lr
[eNL_VDWQQ_WATERWATER
],
283 maxsr_wat
, maxlr_wat
, ivdw
, ivdwmod
, ielec
, ielecmod
, igeometry_ww
, type
, bElecAndVdwSwitchDiffers
);
284 init_nblist(log
, &nbl
->nlist_sr
[eNL_QQ_WATERWATER
], &nbl
->nlist_lr
[eNL_QQ_WATERWATER
],
285 maxsr_wat
, maxlr_wat
, GMX_NBKERNEL_VDW_NONE
, eintmodNONE
, ielec
, ielecmod
, igeometry_ww
, type
, bElecAndVdwSwitchDiffers
);
287 /* Did we get the solvent loops so we can use optimized water kernels? */
288 if (nbl
->nlist_sr
[eNL_VDWQQ_WATER
].kernelptr_vf
== NULL
289 || nbl
->nlist_sr
[eNL_QQ_WATER
].kernelptr_vf
== NULL
290 #ifndef DISABLE_WATERWATER_NLIST
291 || nbl
->nlist_sr
[eNL_VDWQQ_WATERWATER
].kernelptr_vf
== NULL
292 || nbl
->nlist_sr
[eNL_QQ_WATERWATER
].kernelptr_vf
== NULL
296 fr
->solvent_opt
= esolNO
;
299 fprintf(log
, "Note: The available nonbonded kernels do not support water optimization - disabling.\n");
303 if (fr
->efep
!= efepNO
)
305 init_nblist(log
, &nbl
->nlist_sr
[eNL_VDWQQ_FREE
], &nbl
->nlist_lr
[eNL_VDWQQ_FREE
],
306 maxsr
, maxlr
, ivdw
, ivdwmod
, ielec
, ielecmod
, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE
, GMX_NBLIST_INTERACTION_FREE_ENERGY
, bElecAndVdwSwitchDiffers
);
307 init_nblist(log
, &nbl
->nlist_sr
[eNL_VDW_FREE
], &nbl
->nlist_lr
[eNL_VDW_FREE
],
308 maxsr
, maxlr
, ivdw
, ivdwmod
, GMX_NBKERNEL_ELEC_NONE
, eintmodNONE
, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE
, GMX_NBLIST_INTERACTION_FREE_ENERGY
, bElecAndVdwSwitchDiffers
);
309 init_nblist(log
, &nbl
->nlist_sr
[eNL_QQ_FREE
], &nbl
->nlist_lr
[eNL_QQ_FREE
],
310 maxsr
, maxlr
, GMX_NBKERNEL_VDW_NONE
, eintmodNONE
, ielec
, ielecmod
, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE
, GMX_NBLIST_INTERACTION_FREE_ENERGY
, bElecAndVdwSwitchDiffers
);
314 if (fr
->bQMMM
&& fr
->qr
->QMMMscheme
!= eQMMMschemeoniom
)
316 init_nblist(log
, &fr
->QMMMlist
, NULL
,
317 maxsr
, maxlr
, 0, 0, ielec
, ielecmod
, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE
, GMX_NBLIST_INTERACTION_STANDARD
, bElecAndVdwSwitchDiffers
);
325 fr
->ns
.nblist_initialized
= TRUE
;
328 static void reset_nblist(t_nblist
*nl
)
338 static void reset_neighbor_lists(t_forcerec
*fr
, gmx_bool bResetSR
, gmx_bool bResetLR
)
344 /* only reset the short-range nblist */
345 reset_nblist(&(fr
->QMMMlist
));
348 for (n
= 0; n
< fr
->nnblists
; n
++)
350 for (i
= 0; i
< eNL_NR
; i
++)
354 reset_nblist( &(fr
->nblists
[n
].nlist_sr
[i
]) );
358 reset_nblist( &(fr
->nblists
[n
].nlist_lr
[i
]) );
367 static gmx_inline
void new_i_nblist(t_nblist
*nlist
, atom_id i_atom
, int shift
, int gid
)
369 int i
, k
, nri
, nshift
;
373 /* Check whether we have to increase the i counter */
375 (nlist
->iinr
[nri
] != i_atom
) ||
376 (nlist
->shift
[nri
] != shift
) ||
377 (nlist
->gid
[nri
] != gid
))
379 /* This is something else. Now see if any entries have
380 * been added in the list of the previous atom.
383 ((nlist
->jindex
[nri
+1] > nlist
->jindex
[nri
]) &&
384 (nlist
->gid
[nri
] != -1)))
386 /* If so increase the counter */
389 if (nlist
->nri
>= nlist
->maxnri
)
391 nlist
->maxnri
+= over_alloc_large(nlist
->nri
);
392 reallocate_nblist(nlist
);
395 /* Set the number of neighbours and the atom number */
396 nlist
->jindex
[nri
+1] = nlist
->jindex
[nri
];
397 nlist
->iinr
[nri
] = i_atom
;
398 nlist
->gid
[nri
] = gid
;
399 nlist
->shift
[nri
] = shift
;
403 /* Adding to previous list. First remove possible previous padding */
404 if (nlist
->simd_padding_width
> 1)
406 while (nlist
->nrj
> 0 && nlist
->jjnr
[nlist
->nrj
-1] < 0)
414 static gmx_inline
void close_i_nblist(t_nblist
*nlist
)
416 int nri
= nlist
->nri
;
421 /* Add elements up to padding. Since we allocate memory in units
422 * of the simd_padding width, we do not have to check for possible
423 * list reallocation here.
425 while ((nlist
->nrj
% nlist
->simd_padding_width
) != 0)
427 /* Use -4 here, so we can write forces for 4 atoms before real data */
428 nlist
->jjnr
[nlist
->nrj
++] = -4;
430 nlist
->jindex
[nri
+1] = nlist
->nrj
;
432 len
= nlist
->nrj
- nlist
->jindex
[nri
];
433 /* If there are no j-particles we have to reduce the
434 * number of i-particles again, to prevent errors in the
437 if ((len
== 0) && (nlist
->nri
> 0))
444 static gmx_inline
void close_nblist(t_nblist
*nlist
)
446 /* Only close this nblist when it has been initialized.
447 * Avoid the creation of i-lists with no j-particles.
451 /* Some assembly kernels do not support empty lists,
452 * make sure here that we don't generate any empty lists.
453 * With the current ns code this branch is taken in two cases:
454 * No i-particles at all: nri=-1 here
455 * There are i-particles, but no j-particles; nri=0 here
461 /* Close list number nri by incrementing the count */
466 static gmx_inline
void close_neighbor_lists(t_forcerec
*fr
, gmx_bool bMakeQMMMnblist
)
472 close_nblist(&(fr
->QMMMlist
));
475 for (n
= 0; n
< fr
->nnblists
; n
++)
477 for (i
= 0; (i
< eNL_NR
); i
++)
479 close_nblist(&(fr
->nblists
[n
].nlist_sr
[i
]));
480 close_nblist(&(fr
->nblists
[n
].nlist_lr
[i
]));
486 static gmx_inline
void add_j_to_nblist(t_nblist
*nlist
, atom_id j_atom
, gmx_bool bLR
)
488 int nrj
= nlist
->nrj
;
490 if (nlist
->nrj
>= nlist
->maxnrj
)
492 nlist
->maxnrj
= round_up_to_simd_width(over_alloc_small(nlist
->nrj
+ 1), nlist
->simd_padding_width
);
496 fprintf(debug
, "Increasing %s nblist (ielec=%d,ivdw=%d,type=%d,igeometry=%d) j size to %d\n",
497 bLR
? "LR" : "SR", nlist
->ielec
, nlist
->ivdw
, nlist
->type
, nlist
->igeometry
, nlist
->maxnrj
);
500 srenew(nlist
->jjnr
, nlist
->maxnrj
);
503 nlist
->jjnr
[nrj
] = j_atom
;
507 static gmx_inline
void add_j_to_nblist_cg(t_nblist
*nlist
,
508 atom_id j_start
, int j_end
,
509 t_excl
*bexcl
, gmx_bool i_is_j
,
512 int nrj
= nlist
->nrj
;
515 if (nlist
->nrj
>= nlist
->maxnrj
)
517 nlist
->maxnrj
= over_alloc_small(nlist
->nrj
+ 1);
520 fprintf(debug
, "Increasing %s nblist (ielec=%d,ivdw=%d,type=%d,igeometry=%d) j size to %d\n",
521 bLR
? "LR" : "SR", nlist
->ielec
, nlist
->ivdw
, nlist
->type
, nlist
->igeometry
, nlist
->maxnrj
);
524 srenew(nlist
->jjnr
, nlist
->maxnrj
);
525 srenew(nlist
->jjnr_end
, nlist
->maxnrj
);
526 srenew(nlist
->excl
, nlist
->maxnrj
*MAX_CGCGSIZE
);
529 nlist
->jjnr
[nrj
] = j_start
;
530 nlist
->jjnr_end
[nrj
] = j_end
;
532 if (j_end
- j_start
> MAX_CGCGSIZE
)
534 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
);
537 /* Set the exclusions */
538 for (j
= j_start
; j
< j_end
; j
++)
540 nlist
->excl
[nrj
*MAX_CGCGSIZE
+ j
- j_start
] = bexcl
[j
];
544 /* Avoid double counting of intra-cg interactions */
545 for (j
= 1; j
< j_end
-j_start
; j
++)
547 nlist
->excl
[nrj
*MAX_CGCGSIZE
+ j
] |= (1<<j
) - 1;
555 put_in_list_t (gmx_bool bHaveVdW
[],
572 put_in_list_at(gmx_bool bHaveVdW
[],
588 /* The a[] index has been removed,
589 * to put it back in i_atom should be a[i0] and jj should be a[jj].
594 t_nblist
* vdwc_free
= NULL
;
595 t_nblist
* vdw_free
= NULL
;
596 t_nblist
* coul_free
= NULL
;
597 t_nblist
* vdwc_ww
= NULL
;
598 t_nblist
* coul_ww
= NULL
;
600 int i
, j
, jcg
, igid
, gid
, nbl_ind
, ind_ij
;
601 atom_id jj
, jj0
, jj1
, i_atom
;
606 real
*charge
, *chargeB
;
607 real qi
, qiB
, qq
, rlj
;
608 gmx_bool bFreeEnergy
, bFree
, bFreeJ
, bNotEx
, *bPert
;
609 gmx_bool bDoVdW_i
, bDoCoul_i
, bDoCoul_i_sol
;
613 /* Copy some pointers */
615 charge
= md
->chargeA
;
616 chargeB
= md
->chargeB
;
619 bPert
= md
->bPerturbed
;
623 nicg
= index
[icg
+1]-i0
;
625 /* Get the i charge group info */
626 igid
= GET_CGINFO_GID(cginfo
[icg
]);
628 iwater
= (solvent_opt
!= esolNO
) ? GET_CGINFO_SOLOPT(cginfo
[icg
]) : esolNO
;
633 /* Check if any of the particles involved are perturbed.
634 * If not we can do the cheaper normal put_in_list
635 * and use more solvent optimization.
637 for (i
= 0; i
< nicg
; i
++)
639 bFreeEnergy
|= bPert
[i0
+i
];
641 /* Loop over the j charge groups */
642 for (j
= 0; (j
< nj
&& !bFreeEnergy
); j
++)
647 /* Finally loop over the atoms in the j-charge group */
648 for (jj
= jj0
; jj
< jj1
; jj
++)
650 bFreeEnergy
|= bPert
[jj
];
655 /* Unpack pointers to neighbourlist structs */
656 if (fr
->nnblists
== 1)
662 nbl_ind
= fr
->gid2nblists
[GID(igid
, jgid
, ngid
)];
666 nlist
= fr
->nblists
[nbl_ind
].nlist_lr
;
670 nlist
= fr
->nblists
[nbl_ind
].nlist_sr
;
673 if (iwater
!= esolNO
)
675 vdwc
= &nlist
[eNL_VDWQQ_WATER
];
676 vdw
= &nlist
[eNL_VDW
];
677 coul
= &nlist
[eNL_QQ_WATER
];
678 #ifndef DISABLE_WATERWATER_NLIST
679 vdwc_ww
= &nlist
[eNL_VDWQQ_WATERWATER
];
680 coul_ww
= &nlist
[eNL_QQ_WATERWATER
];
685 vdwc
= &nlist
[eNL_VDWQQ
];
686 vdw
= &nlist
[eNL_VDW
];
687 coul
= &nlist
[eNL_QQ
];
692 if (iwater
!= esolNO
)
694 /* Loop over the atoms in the i charge group */
696 gid
= GID(igid
, jgid
, ngid
);
697 /* Create new i_atom for each energy group */
698 if (bDoCoul
&& bDoVdW
)
700 new_i_nblist(vdwc
, i_atom
, shift
, gid
);
701 #ifndef DISABLE_WATERWATER_NLIST
702 new_i_nblist(vdwc_ww
, i_atom
, shift
, gid
);
707 new_i_nblist(vdw
, i_atom
, shift
, gid
);
711 new_i_nblist(coul
, i_atom
, shift
, gid
);
712 #ifndef DISABLE_WATERWATER_NLIST
713 new_i_nblist(coul_ww
, i_atom
, shift
, gid
);
716 /* Loop over the j charge groups */
717 for (j
= 0; (j
< nj
); j
++)
727 jwater
= GET_CGINFO_SOLOPT(cginfo
[jcg
]);
729 if (iwater
== esolSPC
&& jwater
== esolSPC
)
731 /* Interaction between two SPC molecules */
734 /* VdW only - only first atoms in each water interact */
735 add_j_to_nblist(vdw
, jj0
, bLR
);
739 #ifdef DISABLE_WATERWATER_NLIST
740 /* Add entries for the three atoms - only do VdW if we need to */
743 add_j_to_nblist(coul
, jj0
, bLR
);
747 add_j_to_nblist(vdwc
, jj0
, bLR
);
749 add_j_to_nblist(coul
, jj0
+1, bLR
);
750 add_j_to_nblist(coul
, jj0
+2, bLR
);
752 /* One entry for the entire water-water interaction */
755 add_j_to_nblist(coul_ww
, jj0
, bLR
);
759 add_j_to_nblist(vdwc_ww
, jj0
, bLR
);
764 else if (iwater
== esolTIP4P
&& jwater
== esolTIP4P
)
766 /* Interaction between two TIP4p molecules */
769 /* VdW only - only first atoms in each water interact */
770 add_j_to_nblist(vdw
, jj0
, bLR
);
774 #ifdef DISABLE_WATERWATER_NLIST
775 /* Add entries for the four atoms - only do VdW if we need to */
778 add_j_to_nblist(vdw
, jj0
, bLR
);
780 add_j_to_nblist(coul
, jj0
+1, bLR
);
781 add_j_to_nblist(coul
, jj0
+2, bLR
);
782 add_j_to_nblist(coul
, jj0
+3, bLR
);
784 /* One entry for the entire water-water interaction */
787 add_j_to_nblist(coul_ww
, jj0
, bLR
);
791 add_j_to_nblist(vdwc_ww
, jj0
, bLR
);
798 /* j charge group is not water, but i is.
799 * Add entries to the water-other_atom lists; the geometry of the water
800 * molecule doesn't matter - that is taken care of in the nonbonded kernel,
801 * so we don't care if it is SPC or TIP4P...
808 for (jj
= jj0
; (jj
< jj1
); jj
++)
812 add_j_to_nblist(coul
, jj
, bLR
);
818 for (jj
= jj0
; (jj
< jj1
); jj
++)
820 if (bHaveVdW
[type
[jj
]])
822 add_j_to_nblist(vdw
, jj
, bLR
);
828 /* _charge_ _groups_ interact with both coulomb and LJ */
829 /* Check which atoms we should add to the lists! */
830 for (jj
= jj0
; (jj
< jj1
); jj
++)
832 if (bHaveVdW
[type
[jj
]])
836 add_j_to_nblist(vdwc
, jj
, bLR
);
840 add_j_to_nblist(vdw
, jj
, bLR
);
843 else if (charge
[jj
] != 0)
845 add_j_to_nblist(coul
, jj
, bLR
);
852 close_i_nblist(coul
);
853 close_i_nblist(vdwc
);
854 #ifndef DISABLE_WATERWATER_NLIST
855 close_i_nblist(coul_ww
);
856 close_i_nblist(vdwc_ww
);
861 /* no solvent as i charge group */
862 /* Loop over the atoms in the i charge group */
863 for (i
= 0; i
< nicg
; i
++)
866 gid
= GID(igid
, jgid
, ngid
);
869 /* Create new i_atom for each energy group */
870 if (bDoVdW
&& bDoCoul
)
872 new_i_nblist(vdwc
, i_atom
, shift
, gid
);
876 new_i_nblist(vdw
, i_atom
, shift
, gid
);
880 new_i_nblist(coul
, i_atom
, shift
, gid
);
882 bDoVdW_i
= (bDoVdW
&& bHaveVdW
[type
[i_atom
]]);
883 bDoCoul_i
= (bDoCoul
&& qi
!= 0);
885 if (bDoVdW_i
|| bDoCoul_i
)
887 /* Loop over the j charge groups */
888 for (j
= 0; (j
< nj
); j
++)
892 /* Check for large charge groups */
903 /* Finally loop over the atoms in the j-charge group */
904 for (jj
= jj0
; jj
< jj1
; jj
++)
906 bNotEx
= NOTEXCL(bExcl
, i
, jj
);
914 add_j_to_nblist(coul
, jj
, bLR
);
919 if (bHaveVdW
[type
[jj
]])
921 add_j_to_nblist(vdw
, jj
, bLR
);
926 if (bHaveVdW
[type
[jj
]])
930 add_j_to_nblist(vdwc
, jj
, bLR
);
934 add_j_to_nblist(vdw
, jj
, bLR
);
937 else if (charge
[jj
] != 0)
939 add_j_to_nblist(coul
, jj
, bLR
);
947 close_i_nblist(coul
);
948 close_i_nblist(vdwc
);
954 /* we are doing free energy */
955 vdwc_free
= &nlist
[eNL_VDWQQ_FREE
];
956 vdw_free
= &nlist
[eNL_VDW_FREE
];
957 coul_free
= &nlist
[eNL_QQ_FREE
];
958 /* Loop over the atoms in the i charge group */
959 for (i
= 0; i
< nicg
; i
++)
962 gid
= GID(igid
, jgid
, ngid
);
964 qiB
= chargeB
[i_atom
];
966 /* Create new i_atom for each energy group */
967 if (bDoVdW
&& bDoCoul
)
969 new_i_nblist(vdwc
, i_atom
, shift
, gid
);
973 new_i_nblist(vdw
, i_atom
, shift
, gid
);
977 new_i_nblist(coul
, i_atom
, shift
, gid
);
980 new_i_nblist(vdw_free
, i_atom
, shift
, gid
);
981 new_i_nblist(coul_free
, i_atom
, shift
, gid
);
982 new_i_nblist(vdwc_free
, i_atom
, shift
, gid
);
984 bDoVdW_i
= (bDoVdW
&&
985 (bHaveVdW
[type
[i_atom
]] || bHaveVdW
[typeB
[i_atom
]]));
986 bDoCoul_i
= (bDoCoul
&& (qi
!= 0 || qiB
!= 0));
987 /* For TIP4P the first atom does not have a charge,
988 * but the last three do. So we should still put an atom
989 * without LJ but with charge in the water-atom neighborlist
990 * for a TIP4p i charge group.
991 * For SPC type water the first atom has LJ and charge,
992 * so there is no such problem.
994 if (iwater
== esolNO
)
996 bDoCoul_i_sol
= bDoCoul_i
;
1000 bDoCoul_i_sol
= bDoCoul
;
1003 if (bDoVdW_i
|| bDoCoul_i_sol
)
1005 /* Loop over the j charge groups */
1006 for (j
= 0; (j
< nj
); j
++)
1010 /* Check for large charge groups */
1021 /* Finally loop over the atoms in the j-charge group */
1022 bFree
= bPert
[i_atom
];
1023 for (jj
= jj0
; (jj
< jj1
); jj
++)
1025 bFreeJ
= bFree
|| bPert
[jj
];
1026 /* Complicated if, because the water H's should also
1027 * see perturbed j-particles
1029 if (iwater
== esolNO
|| i
== 0 || bFreeJ
)
1031 bNotEx
= NOTEXCL(bExcl
, i
, jj
);
1039 if (charge
[jj
] != 0 || chargeB
[jj
] != 0)
1041 add_j_to_nblist(coul_free
, jj
, bLR
);
1044 else if (!bDoCoul_i
)
1046 if (bHaveVdW
[type
[jj
]] || bHaveVdW
[typeB
[jj
]])
1048 add_j_to_nblist(vdw_free
, jj
, bLR
);
1053 if (bHaveVdW
[type
[jj
]] || bHaveVdW
[typeB
[jj
]])
1055 if (charge
[jj
] != 0 || chargeB
[jj
] != 0)
1057 add_j_to_nblist(vdwc_free
, jj
, bLR
);
1061 add_j_to_nblist(vdw_free
, jj
, bLR
);
1064 else if (charge
[jj
] != 0 || chargeB
[jj
] != 0)
1066 add_j_to_nblist(coul_free
, jj
, bLR
);
1072 /* This is done whether or not bWater is set */
1073 if (charge
[jj
] != 0)
1075 add_j_to_nblist(coul
, jj
, bLR
);
1078 else if (!bDoCoul_i_sol
)
1080 if (bHaveVdW
[type
[jj
]])
1082 add_j_to_nblist(vdw
, jj
, bLR
);
1087 if (bHaveVdW
[type
[jj
]])
1089 if (charge
[jj
] != 0)
1091 add_j_to_nblist(vdwc
, jj
, bLR
);
1095 add_j_to_nblist(vdw
, jj
, bLR
);
1098 else if (charge
[jj
] != 0)
1100 add_j_to_nblist(coul
, jj
, bLR
);
1108 close_i_nblist(vdw
);
1109 close_i_nblist(coul
);
1110 close_i_nblist(vdwc
);
1111 close_i_nblist(vdw_free
);
1112 close_i_nblist(coul_free
);
1113 close_i_nblist(vdwc_free
);
1119 put_in_list_adress(gmx_bool bHaveVdW
[],
1135 /* The a[] index has been removed,
1136 * to put it back in i_atom should be a[i0] and jj should be a[jj].
1141 t_nblist
* vdwc_adress
= NULL
;
1142 t_nblist
* vdw_adress
= NULL
;
1143 t_nblist
* coul_adress
= NULL
;
1144 t_nblist
* vdwc_ww
= NULL
;
1145 t_nblist
* coul_ww
= NULL
;
1147 int i
, j
, jcg
, igid
, gid
, nbl_ind
, nbl_ind_adress
;
1148 atom_id jj
, jj0
, jj1
, i_atom
;
1153 real
*charge
, *chargeB
;
1155 real qi
, qiB
, qq
, rlj
;
1156 gmx_bool bFreeEnergy
, bFree
, bFreeJ
, bNotEx
, *bPert
;
1157 gmx_bool bDoVdW_i
, bDoCoul_i
, bDoCoul_i_sol
;
1159 gmx_bool j_all_atom
;
1161 t_nblist
*nlist
, *nlist_adress
;
1162 gmx_bool bEnergyGroupCG
;
1164 /* Copy some pointers */
1165 cginfo
= fr
->cginfo
;
1166 charge
= md
->chargeA
;
1167 chargeB
= md
->chargeB
;
1170 bPert
= md
->bPerturbed
;
1173 /* Get atom range */
1175 nicg
= index
[icg
+1]-i0
;
1177 /* Get the i charge group info */
1178 igid
= GET_CGINFO_GID(cginfo
[icg
]);
1180 iwater
= (solvent_opt
!= esolNO
) ? GET_CGINFO_SOLOPT(cginfo
[icg
]) : esolNO
;
1184 gmx_fatal(FARGS
, "AdResS does not support free energy pertubation\n");
1187 /* Unpack pointers to neighbourlist structs */
1188 if (fr
->nnblists
== 2)
1195 nbl_ind
= fr
->gid2nblists
[GID(igid
, jgid
, ngid
)];
1196 nbl_ind_adress
= nbl_ind
+fr
->nnblists
/2;
1200 nlist
= fr
->nblists
[nbl_ind
].nlist_lr
;
1201 nlist_adress
= fr
->nblists
[nbl_ind_adress
].nlist_lr
;
1205 nlist
= fr
->nblists
[nbl_ind
].nlist_sr
;
1206 nlist_adress
= fr
->nblists
[nbl_ind_adress
].nlist_sr
;
1210 vdwc
= &nlist
[eNL_VDWQQ
];
1211 vdw
= &nlist
[eNL_VDW
];
1212 coul
= &nlist
[eNL_QQ
];
1214 vdwc_adress
= &nlist_adress
[eNL_VDWQQ
];
1215 vdw_adress
= &nlist_adress
[eNL_VDW
];
1216 coul_adress
= &nlist_adress
[eNL_QQ
];
1218 /* We do not support solvent optimization with AdResS for now.
1219 For this we would need hybrid solvent-other kernels */
1221 /* no solvent as i charge group */
1222 /* Loop over the atoms in the i charge group */
1223 for (i
= 0; i
< nicg
; i
++)
1226 gid
= GID(igid
, jgid
, ngid
);
1227 qi
= charge
[i_atom
];
1229 /* Create new i_atom for each energy group */
1230 if (bDoVdW
&& bDoCoul
)
1232 new_i_nblist(vdwc
, i_atom
, shift
, gid
);
1233 new_i_nblist(vdwc_adress
, i_atom
, shift
, gid
);
1238 new_i_nblist(vdw
, i_atom
, shift
, gid
);
1239 new_i_nblist(vdw_adress
, i_atom
, shift
, gid
);
1244 new_i_nblist(coul
, i_atom
, shift
, gid
);
1245 new_i_nblist(coul_adress
, i_atom
, shift
, gid
);
1247 bDoVdW_i
= (bDoVdW
&& bHaveVdW
[type
[i_atom
]]);
1248 bDoCoul_i
= (bDoCoul
&& qi
!= 0);
1250 /* Here we find out whether the energy groups interaction belong to a
1251 * coarse-grained (vsite) or atomistic interaction. Note that, beacuse
1252 * interactions between coarse-grained and other (atomistic) energygroups
1253 * are excluded automatically by grompp, it is sufficient to check for
1254 * the group id of atom i (igid) */
1255 bEnergyGroupCG
= !egp_explicit(fr
, igid
);
1257 if (bDoVdW_i
|| bDoCoul_i
)
1259 /* Loop over the j charge groups */
1260 for (j
= 0; (j
< nj
); j
++)
1264 /* Check for large charge groups */
1275 /* Finally loop over the atoms in the j-charge group */
1276 for (jj
= jj0
; jj
< jj1
; jj
++)
1278 bNotEx
= NOTEXCL(bExcl
, i
, jj
);
1280 /* Now we have to exclude interactions which will be zero
1281 * anyway due to the AdResS weights (in previous implementations
1282 * this was done in the force kernel). This is necessary as
1283 * pure interactions (those with b_hybrid=false, i.e. w_i*w_j==1 or 0)
1284 * are put into neighbour lists which will be passed to the
1285 * standard (optimized) kernels for speed. The interactions with
1286 * b_hybrid=true are placed into the _adress neighbour lists and
1287 * processed by the generic AdResS kernel.
1289 if ( (bEnergyGroupCG
&&
1290 wf
[i_atom
] >= 1-GMX_REAL_EPS
&& wf
[jj
] >= 1-GMX_REAL_EPS
) ||
1291 ( !bEnergyGroupCG
&& wf
[jj
] <= GMX_REAL_EPS
) )
1296 b_hybrid
= !((wf
[i_atom
] >= 1-GMX_REAL_EPS
&& wf
[jj
] >= 1-GMX_REAL_EPS
) ||
1297 (wf
[i_atom
] <= GMX_REAL_EPS
&& wf
[jj
] <= GMX_REAL_EPS
));
1303 if (charge
[jj
] != 0)
1307 add_j_to_nblist(coul
, jj
, bLR
);
1311 add_j_to_nblist(coul_adress
, jj
, bLR
);
1315 else if (!bDoCoul_i
)
1317 if (bHaveVdW
[type
[jj
]])
1321 add_j_to_nblist(vdw
, jj
, bLR
);
1325 add_j_to_nblist(vdw_adress
, jj
, bLR
);
1331 if (bHaveVdW
[type
[jj
]])
1333 if (charge
[jj
] != 0)
1337 add_j_to_nblist(vdwc
, jj
, bLR
);
1341 add_j_to_nblist(vdwc_adress
, jj
, bLR
);
1348 add_j_to_nblist(vdw
, jj
, bLR
);
1352 add_j_to_nblist(vdw_adress
, jj
, bLR
);
1357 else if (charge
[jj
] != 0)
1361 add_j_to_nblist(coul
, jj
, bLR
);
1365 add_j_to_nblist(coul_adress
, jj
, bLR
);
1374 close_i_nblist(vdw
);
1375 close_i_nblist(coul
);
1376 close_i_nblist(vdwc
);
1377 close_i_nblist(vdw_adress
);
1378 close_i_nblist(coul_adress
);
1379 close_i_nblist(vdwc_adress
);
1385 put_in_list_qmmm(gmx_bool gmx_unused bHaveVdW
[],
1387 t_mdatoms gmx_unused
* md
,
1397 gmx_bool gmx_unused bDoVdW
,
1398 gmx_bool gmx_unused bDoCoul
,
1399 int gmx_unused solvent_opt
)
1402 int i
, j
, jcg
, igid
, gid
;
1403 atom_id jj
, jj0
, jj1
, i_atom
;
1407 /* Get atom range */
1409 nicg
= index
[icg
+1]-i0
;
1411 /* Get the i charge group info */
1412 igid
= GET_CGINFO_GID(fr
->cginfo
[icg
]);
1414 coul
= &fr
->QMMMlist
;
1416 /* Loop over atoms in the ith charge group */
1417 for (i
= 0; i
< nicg
; i
++)
1420 gid
= GID(igid
, jgid
, ngid
);
1421 /* Create new i_atom for each energy group */
1422 new_i_nblist(coul
, i_atom
, shift
, gid
);
1424 /* Loop over the j charge groups */
1425 for (j
= 0; j
< nj
; j
++)
1429 /* Charge groups cannot have QM and MM atoms simultaneously */
1434 /* Finally loop over the atoms in the j-charge group */
1435 for (jj
= jj0
; jj
< jj1
; jj
++)
1437 bNotEx
= NOTEXCL(bExcl
, i
, jj
);
1440 add_j_to_nblist(coul
, jj
, bLR
);
1445 close_i_nblist(coul
);
1450 put_in_list_cg(gmx_bool gmx_unused bHaveVdW
[],
1452 t_mdatoms gmx_unused
* md
,
1462 gmx_bool gmx_unused bDoVdW
,
1463 gmx_bool gmx_unused bDoCoul
,
1464 int gmx_unused solvent_opt
)
1467 int igid
, gid
, nbl_ind
;
1471 cginfo
= fr
->cginfo
[icg
];
1473 igid
= GET_CGINFO_GID(cginfo
);
1474 gid
= GID(igid
, jgid
, ngid
);
1476 /* Unpack pointers to neighbourlist structs */
1477 if (fr
->nnblists
== 1)
1483 nbl_ind
= fr
->gid2nblists
[gid
];
1487 vdwc
= &fr
->nblists
[nbl_ind
].nlist_lr
[eNL_VDWQQ
];
1491 vdwc
= &fr
->nblists
[nbl_ind
].nlist_sr
[eNL_VDWQQ
];
1494 /* Make a new neighbor list for charge group icg.
1495 * Currently simply one neighbor list is made with LJ and Coulomb.
1496 * If required, zero interactions could be removed here
1497 * or in the force loop.
1499 new_i_nblist(vdwc
, index
[icg
], shift
, gid
);
1500 vdwc
->iinr_end
[vdwc
->nri
] = index
[icg
+1];
1502 for (j
= 0; (j
< nj
); j
++)
1505 /* Skip the icg-icg pairs if all self interactions are excluded */
1506 if (!(jcg
== icg
&& GET_CGINFO_EXCL_INTRA(cginfo
)))
1508 /* Here we add the j charge group jcg to the list,
1509 * exclusions are also added to the list.
1511 add_j_to_nblist_cg(vdwc
, index
[jcg
], index
[jcg
+1], bExcl
, icg
== jcg
, bLR
);
1515 close_i_nblist(vdwc
);
1518 static void setexcl(atom_id start
, atom_id end
, t_blocka
*excl
, gmx_bool b
,
1525 for (i
= start
; i
< end
; i
++)
1527 for (k
= excl
->index
[i
]; k
< excl
->index
[i
+1]; k
++)
1529 SETEXCL(bexcl
, i
-start
, excl
->a
[k
]);
1535 for (i
= start
; i
< end
; i
++)
1537 for (k
= excl
->index
[i
]; k
< excl
->index
[i
+1]; k
++)
1539 RMEXCL(bexcl
, i
-start
, excl
->a
[k
]);
1545 int calc_naaj(int icg
, int cgtot
)
1549 if ((cgtot
% 2) == 1)
1551 /* Odd number of charge groups, easy */
1552 naaj
= 1 + (cgtot
/2);
1554 else if ((cgtot
% 4) == 0)
1556 /* Multiple of four is hard */
1593 fprintf(log
, "naaj=%d\n", naaj
);
1599 /************************************************
1601 * S I M P L E C O R E S T U F F
1603 ************************************************/
1605 static real
calc_image_tric(rvec xi
, rvec xj
, matrix box
,
1606 rvec b_inv
, int *shift
)
1608 /* This code assumes that the cut-off is smaller than
1609 * a half times the smallest diagonal element of the box.
1611 const real h25
= 2.5;
1616 /* Compute diff vector */
1617 dz
= xj
[ZZ
] - xi
[ZZ
];
1618 dy
= xj
[YY
] - xi
[YY
];
1619 dx
= xj
[XX
] - xi
[XX
];
1621 /* Perform NINT operation, using trunc operation, therefore
1622 * we first add 2.5 then subtract 2 again
1624 tz
= dz
*b_inv
[ZZ
] + h25
;
1626 dz
-= tz
*box
[ZZ
][ZZ
];
1627 dy
-= tz
*box
[ZZ
][YY
];
1628 dx
-= tz
*box
[ZZ
][XX
];
1630 ty
= dy
*b_inv
[YY
] + h25
;
1632 dy
-= ty
*box
[YY
][YY
];
1633 dx
-= ty
*box
[YY
][XX
];
1635 tx
= dx
*b_inv
[XX
]+h25
;
1637 dx
-= tx
*box
[XX
][XX
];
1639 /* Distance squared */
1640 r2
= (dx
*dx
) + (dy
*dy
) + (dz
*dz
);
1642 *shift
= XYZ2IS(tx
, ty
, tz
);
1647 static real
calc_image_rect(rvec xi
, rvec xj
, rvec box_size
,
1648 rvec b_inv
, int *shift
)
1650 const real h15
= 1.5;
1656 /* Compute diff vector */
1657 dx
= xj
[XX
] - xi
[XX
];
1658 dy
= xj
[YY
] - xi
[YY
];
1659 dz
= xj
[ZZ
] - xi
[ZZ
];
1661 /* Perform NINT operation, using trunc operation, therefore
1662 * we first add 1.5 then subtract 1 again
1664 tx
= dx
*b_inv
[XX
] + h15
;
1665 ty
= dy
*b_inv
[YY
] + h15
;
1666 tz
= dz
*b_inv
[ZZ
] + h15
;
1671 /* Correct diff vector for translation */
1672 ddx
= tx
*box_size
[XX
] - dx
;
1673 ddy
= ty
*box_size
[YY
] - dy
;
1674 ddz
= tz
*box_size
[ZZ
] - dz
;
1676 /* Distance squared */
1677 r2
= (ddx
*ddx
) + (ddy
*ddy
) + (ddz
*ddz
);
1679 *shift
= XYZ2IS(tx
, ty
, tz
);
1684 static void add_simple(t_ns_buf
*nsbuf
, int nrj
, atom_id cg_j
,
1685 gmx_bool bHaveVdW
[], int ngid
, t_mdatoms
*md
,
1686 int icg
, int jgid
, t_block
*cgs
, t_excl bexcl
[],
1687 int shift
, t_forcerec
*fr
, put_in_list_t
*put_in_list
)
1689 if (nsbuf
->nj
+ nrj
> MAX_CG
)
1691 put_in_list(bHaveVdW
, ngid
, md
, icg
, jgid
, nsbuf
->ncg
, nsbuf
->jcg
,
1692 cgs
->index
, bexcl
, shift
, fr
, FALSE
, TRUE
, TRUE
, fr
->solvent_opt
);
1693 /* Reset buffer contents */
1694 nsbuf
->ncg
= nsbuf
->nj
= 0;
1696 nsbuf
->jcg
[nsbuf
->ncg
++] = cg_j
;
1700 static void ns_inner_tric(rvec x
[], int icg
, int *i_egp_flags
,
1701 int njcg
, atom_id jcg
[],
1702 matrix box
, rvec b_inv
, real rcut2
,
1703 t_block
*cgs
, t_ns_buf
**ns_buf
,
1704 gmx_bool bHaveVdW
[], int ngid
, t_mdatoms
*md
,
1705 t_excl bexcl
[], t_forcerec
*fr
,
1706 put_in_list_t
*put_in_list
)
1710 int *cginfo
= fr
->cginfo
;
1711 atom_id cg_j
, *cgindex
;
1714 cgindex
= cgs
->index
;
1716 for (j
= 0; (j
< njcg
); j
++)
1719 nrj
= cgindex
[cg_j
+1]-cgindex
[cg_j
];
1720 if (calc_image_tric(x
[icg
], x
[cg_j
], box
, b_inv
, &shift
) < rcut2
)
1722 jgid
= GET_CGINFO_GID(cginfo
[cg_j
]);
1723 if (!(i_egp_flags
[jgid
] & EGP_EXCL
))
1725 add_simple(&ns_buf
[jgid
][shift
], nrj
, cg_j
,
1726 bHaveVdW
, ngid
, md
, icg
, jgid
, cgs
, bexcl
, shift
, fr
,
1733 static void ns_inner_rect(rvec x
[], int icg
, int *i_egp_flags
,
1734 int njcg
, atom_id jcg
[],
1735 gmx_bool bBox
, rvec box_size
, rvec b_inv
, real rcut2
,
1736 t_block
*cgs
, t_ns_buf
**ns_buf
,
1737 gmx_bool bHaveVdW
[], int ngid
, t_mdatoms
*md
,
1738 t_excl bexcl
[], t_forcerec
*fr
,
1739 put_in_list_t
*put_in_list
)
1743 int *cginfo
= fr
->cginfo
;
1744 atom_id cg_j
, *cgindex
;
1747 cgindex
= cgs
->index
;
1751 for (j
= 0; (j
< njcg
); j
++)
1754 nrj
= cgindex
[cg_j
+1]-cgindex
[cg_j
];
1755 if (calc_image_rect(x
[icg
], x
[cg_j
], box_size
, b_inv
, &shift
) < rcut2
)
1757 jgid
= GET_CGINFO_GID(cginfo
[cg_j
]);
1758 if (!(i_egp_flags
[jgid
] & EGP_EXCL
))
1760 add_simple(&ns_buf
[jgid
][shift
], nrj
, cg_j
,
1761 bHaveVdW
, ngid
, md
, icg
, jgid
, cgs
, bexcl
, shift
, fr
,
1769 for (j
= 0; (j
< njcg
); j
++)
1772 nrj
= cgindex
[cg_j
+1]-cgindex
[cg_j
];
1773 if ((rcut2
== 0) || (distance2(x
[icg
], x
[cg_j
]) < rcut2
))
1775 jgid
= GET_CGINFO_GID(cginfo
[cg_j
]);
1776 if (!(i_egp_flags
[jgid
] & EGP_EXCL
))
1778 add_simple(&ns_buf
[jgid
][CENTRAL
], nrj
, cg_j
,
1779 bHaveVdW
, ngid
, md
, icg
, jgid
, cgs
, bexcl
, CENTRAL
, fr
,
1787 /* ns_simple_core needs to be adapted for QMMM still 2005 */
1789 static int ns_simple_core(t_forcerec
*fr
,
1790 gmx_localtop_t
*top
,
1792 matrix box
, rvec box_size
,
1793 t_excl bexcl
[], atom_id
*aaj
,
1794 int ngid
, t_ns_buf
**ns_buf
,
1795 put_in_list_t
*put_in_list
, gmx_bool bHaveVdW
[])
1799 int nsearch
, icg
, jcg
, igid
, i0
, nri
, nn
;
1802 /* atom_id *i_atoms; */
1803 t_block
*cgs
= &(top
->cgs
);
1804 t_blocka
*excl
= &(top
->excls
);
1807 gmx_bool bBox
, bTriclinic
;
1810 rlist2
= sqr(fr
->rlist
);
1812 bBox
= (fr
->ePBC
!= epbcNONE
);
1815 for (m
= 0; (m
< DIM
); m
++)
1817 b_inv
[m
] = divide_err(1.0, box_size
[m
]);
1819 bTriclinic
= TRICLINIC(box
);
1826 cginfo
= fr
->cginfo
;
1829 for (icg
= fr
->cg0
; (icg
< fr
->hcg
); icg
++)
1832 i0 = cgs->index[icg];
1833 nri = cgs->index[icg+1]-i0;
1834 i_atoms = &(cgs->a[i0]);
1835 i_eg_excl = fr->eg_excl + ngid*md->cENER[*i_atoms];
1836 setexcl(nri,i_atoms,excl,TRUE,bexcl);
1838 igid
= GET_CGINFO_GID(cginfo
[icg
]);
1839 i_egp_flags
= fr
->egp_flags
+ ngid
*igid
;
1840 setexcl(cgs
->index
[icg
], cgs
->index
[icg
+1], excl
, TRUE
, bexcl
);
1842 naaj
= calc_naaj(icg
, cgs
->nr
);
1845 ns_inner_tric(fr
->cg_cm
, icg
, i_egp_flags
, naaj
, &(aaj
[icg
]),
1846 box
, b_inv
, rlist2
, cgs
, ns_buf
,
1847 bHaveVdW
, ngid
, md
, bexcl
, fr
, put_in_list
);
1851 ns_inner_rect(fr
->cg_cm
, icg
, i_egp_flags
, naaj
, &(aaj
[icg
]),
1852 bBox
, box_size
, b_inv
, rlist2
, cgs
, ns_buf
,
1853 bHaveVdW
, ngid
, md
, bexcl
, fr
, put_in_list
);
1857 for (nn
= 0; (nn
< ngid
); nn
++)
1859 for (k
= 0; (k
< SHIFTS
); k
++)
1861 nsbuf
= &(ns_buf
[nn
][k
]);
1864 put_in_list(bHaveVdW
, ngid
, md
, icg
, nn
, nsbuf
->ncg
, nsbuf
->jcg
,
1865 cgs
->index
, bexcl
, k
, fr
, FALSE
, TRUE
, TRUE
, fr
->solvent_opt
);
1866 nsbuf
->ncg
= nsbuf
->nj
= 0;
1870 /* setexcl(nri,i_atoms,excl,FALSE,bexcl); */
1871 setexcl(cgs
->index
[icg
], cgs
->index
[icg
+1], excl
, FALSE
, bexcl
);
1873 close_neighbor_lists(fr
, FALSE
);
1878 /************************************************
1880 * N S 5 G R I D S T U F F
1882 ************************************************/
1884 static gmx_inline
void get_dx(int Nx
, real gridx
, real rc2
, int xgi
, real x
,
1885 int *dx0
, int *dx1
, real
*dcx2
)
1913 for (i
= xgi0
; i
>= 0; i
--)
1915 dcx
= (i
+1)*gridx
-x
;
1924 for (i
= xgi1
; i
< Nx
; i
++)
1937 static gmx_inline
void get_dx_dd(int Nx
, real gridx
, real rc2
, int xgi
, real x
,
1938 int ncpddc
, int shift_min
, int shift_max
,
1939 int *g0
, int *g1
, real
*dcx2
)
1942 int g_min
, g_max
, shift_home
;
1975 g_min
= (shift_min
== shift_home
? 0 : ncpddc
);
1976 g_max
= (shift_max
== shift_home
? ncpddc
- 1 : Nx
- 1);
1983 else if (shift_max
< 0)
1998 /* Check one grid cell down */
1999 dcx
= ((*g0
- 1) + 1)*gridx
- x
;
2011 /* Check one grid cell up */
2012 dcx
= (*g1
+ 1)*gridx
- x
;
2024 #define sqr(x) ((x)*(x))
2025 #define calc_dx2(XI, YI, ZI, y) (sqr(XI-y[XX]) + sqr(YI-y[YY]) + sqr(ZI-y[ZZ]))
2026 #define calc_cyl_dx2(XI, YI, y) (sqr(XI-y[XX]) + sqr(YI-y[YY]))
2027 /****************************************************
2029 * F A S T N E I G H B O R S E A R C H I N G
2031 * Optimized neighboursearching routine using grid
2032 * at least 1x1x1, see GROMACS manual
2034 ****************************************************/
2037 static void get_cutoff2(t_forcerec
*fr
, gmx_bool bDoLongRange
,
2038 real
*rvdw2
, real
*rcoul2
,
2039 real
*rs2
, real
*rm2
, real
*rl2
)
2041 *rs2
= sqr(fr
->rlist
);
2043 if (bDoLongRange
&& fr
->bTwinRange
)
2045 /* With plain cut-off or RF we need to make the list exactly
2046 * up to the cut-off and the cut-off's can be different,
2047 * so we can not simply set them to rlistlong.
2048 * To keep this code compatible with (exotic) old cases,
2049 * we also create lists up to rvdw/rcoulomb for PME and Ewald.
2050 * The interaction check should correspond to:
2051 * !ir_vdw/coulomb_might_be_zero_at_cutoff from inputrec.c.
2053 if (((fr
->vdwtype
== evdwCUT
|| fr
->vdwtype
== evdwPME
) &&
2054 fr
->vdw_modifier
== eintmodNONE
) ||
2055 fr
->rvdw
<= fr
->rlist
)
2057 *rvdw2
= sqr(fr
->rvdw
);
2061 *rvdw2
= sqr(fr
->rlistlong
);
2063 if (((fr
->eeltype
== eelCUT
||
2064 (EEL_RF(fr
->eeltype
) && fr
->eeltype
!= eelRF_ZERO
) ||
2065 fr
->eeltype
== eelPME
||
2066 fr
->eeltype
== eelEWALD
) &&
2067 fr
->coulomb_modifier
== eintmodNONE
) ||
2068 fr
->rcoulomb
<= fr
->rlist
)
2070 *rcoul2
= sqr(fr
->rcoulomb
);
2074 *rcoul2
= sqr(fr
->rlistlong
);
2079 /* Workaround for a gcc -O3 or -ffast-math problem */
2083 *rm2
= min(*rvdw2
, *rcoul2
);
2084 *rl2
= max(*rvdw2
, *rcoul2
);
2087 static void init_nsgrid_lists(t_forcerec
*fr
, int ngid
, gmx_ns_t
*ns
)
2089 real rvdw2
, rcoul2
, rs2
, rm2
, rl2
;
2092 get_cutoff2(fr
, TRUE
, &rvdw2
, &rcoul2
, &rs2
, &rm2
, &rl2
);
2094 /* Short range buffers */
2095 snew(ns
->nl_sr
, ngid
);
2097 snew(ns
->nsr
, ngid
);
2098 snew(ns
->nlr_ljc
, ngid
);
2099 snew(ns
->nlr_one
, ngid
);
2101 /* Always allocate both list types, since rcoulomb might now change with PME load balancing */
2102 /* Long range VdW and Coul buffers */
2103 snew(ns
->nl_lr_ljc
, ngid
);
2104 /* Long range VdW or Coul only buffers */
2105 snew(ns
->nl_lr_one
, ngid
);
2107 for (j
= 0; (j
< ngid
); j
++)
2109 snew(ns
->nl_sr
[j
], MAX_CG
);
2110 snew(ns
->nl_lr_ljc
[j
], MAX_CG
);
2111 snew(ns
->nl_lr_one
[j
], MAX_CG
);
2116 "ns5_core: rs2 = %g, rm2 = %g, rl2 = %g (nm^2)\n",
2121 static int nsgrid_core(t_commrec
*cr
, t_forcerec
*fr
,
2122 matrix box
, int ngid
,
2123 gmx_localtop_t
*top
,
2125 t_excl bexcl
[], gmx_bool
*bExcludeAlleg
,
2127 put_in_list_t
*put_in_list
,
2128 gmx_bool bHaveVdW
[],
2129 gmx_bool bDoLongRange
, gmx_bool bMakeQMMMnblist
)
2132 atom_id
**nl_lr_ljc
, **nl_lr_one
, **nl_sr
;
2133 int *nlr_ljc
, *nlr_one
, *nsr
;
2134 gmx_domdec_t
*dd
= NULL
;
2135 t_block
*cgs
= &(top
->cgs
);
2136 int *cginfo
= fr
->cginfo
;
2137 /* atom_id *i_atoms,*cgsindex=cgs->index; */
2139 int cell_x
, cell_y
, cell_z
;
2140 int d
, tx
, ty
, tz
, dx
, dy
, dz
, cj
;
2141 #ifdef ALLOW_OFFDIAG_LT_HALFDIAG
2142 int zsh_ty
, zsh_tx
, ysh_tx
;
2144 int dx0
, dx1
, dy0
, dy1
, dz0
, dz1
;
2145 int Nx
, Ny
, Nz
, shift
= -1, j
, nrj
, nns
, nn
= -1;
2146 real gridx
, gridy
, gridz
, grid_x
, grid_y
, grid_z
;
2147 real
*dcx2
, *dcy2
, *dcz2
;
2149 int cg0
, cg1
, icg
= -1, cgsnr
, i0
, igid
, nri
, naaj
, max_jcg
;
2150 int jcg0
, jcg1
, jjcg
, cgj0
, jgid
;
2151 int *grida
, *gridnra
, *gridind
;
2152 gmx_bool rvdw_lt_rcoul
, rcoul_lt_rvdw
;
2153 rvec xi
, *cgcm
, grid_offset
;
2154 real r2
, rs2
, rvdw2
, rcoul2
, rm2
, rl2
, XI
, YI
, ZI
, dcx
, dcy
, dcz
, tmp1
, tmp2
;
2156 gmx_bool bDomDec
, bTriclinicX
, bTriclinicY
;
2161 bDomDec
= DOMAINDECOMP(cr
);
2167 bTriclinicX
= ((YY
< grid
->npbcdim
&&
2168 (!bDomDec
|| dd
->nc
[YY
] == 1) && box
[YY
][XX
] != 0) ||
2169 (ZZ
< grid
->npbcdim
&&
2170 (!bDomDec
|| dd
->nc
[ZZ
] == 1) && box
[ZZ
][XX
] != 0));
2171 bTriclinicY
= (ZZ
< grid
->npbcdim
&&
2172 (!bDomDec
|| dd
->nc
[ZZ
] == 1) && box
[ZZ
][YY
] != 0);
2176 get_cutoff2(fr
, bDoLongRange
, &rvdw2
, &rcoul2
, &rs2
, &rm2
, &rl2
);
2178 rvdw_lt_rcoul
= (rvdw2
>= rcoul2
);
2179 rcoul_lt_rvdw
= (rcoul2
>= rvdw2
);
2181 if (bMakeQMMMnblist
)
2189 nl_lr_ljc
= ns
->nl_lr_ljc
;
2190 nl_lr_one
= ns
->nl_lr_one
;
2191 nlr_ljc
= ns
->nlr_ljc
;
2192 nlr_one
= ns
->nlr_one
;
2200 gridind
= grid
->index
;
2201 gridnra
= grid
->nra
;
2204 gridx
= grid
->cell_size
[XX
];
2205 gridy
= grid
->cell_size
[YY
];
2206 gridz
= grid
->cell_size
[ZZ
];
2210 copy_rvec(grid
->cell_offset
, grid_offset
);
2211 copy_ivec(grid
->ncpddc
, ncpddc
);
2216 #ifdef ALLOW_OFFDIAG_LT_HALFDIAG
2217 zsh_ty
= floor(-box
[ZZ
][YY
]/box
[YY
][YY
]+0.5);
2218 zsh_tx
= floor(-box
[ZZ
][XX
]/box
[XX
][XX
]+0.5);
2219 ysh_tx
= floor(-box
[YY
][XX
]/box
[XX
][XX
]+0.5);
2220 if (zsh_tx
!= 0 && ysh_tx
!= 0)
2222 /* This could happen due to rounding, when both ratios are 0.5 */
2231 /* We only want a list for the test particle */
2240 /* Set the shift range */
2241 for (d
= 0; d
< DIM
; d
++)
2245 /* Check if we need periodicity shifts.
2246 * Without PBC or with domain decomposition we don't need them.
2248 if (d
>= ePBC2npbcdim(fr
->ePBC
) || (bDomDec
&& dd
->nc
[d
] > 1))
2255 box
[XX
][XX
] - fabs(box
[YY
][XX
]) - fabs(box
[ZZ
][XX
]) < sqrt(rl2
))
2266 /* Loop over charge groups */
2267 for (icg
= cg0
; (icg
< cg1
); icg
++)
2269 igid
= GET_CGINFO_GID(cginfo
[icg
]);
2270 /* Skip this charge group if all energy groups are excluded! */
2271 if (bExcludeAlleg
[igid
])
2276 i0
= cgs
->index
[icg
];
2278 if (bMakeQMMMnblist
)
2280 /* Skip this charge group if it is not a QM atom while making a
2281 * QM/MM neighbourlist
2283 if (md
->bQM
[i0
] == FALSE
)
2285 continue; /* MM particle, go to next particle */
2288 /* Compute the number of charge groups that fall within the control
2291 naaj
= calc_naaj(icg
, cgsnr
);
2298 /* make a normal neighbourlist */
2302 /* Get the j charge-group and dd cell shift ranges */
2303 dd_get_ns_ranges(cr
->dd
, icg
, &jcg0
, &jcg1
, sh0
, sh1
);
2308 /* Compute the number of charge groups that fall within the control
2311 naaj
= calc_naaj(icg
, cgsnr
);
2317 /* The i-particle is awlways the test particle,
2318 * so we want all j-particles
2320 max_jcg
= cgsnr
- 1;
2324 max_jcg
= jcg1
- cgsnr
;
2329 i_egp_flags
= fr
->egp_flags
+ igid
*ngid
;
2331 /* Set the exclusions for the atoms in charge group icg using a bitmask */
2332 setexcl(i0
, cgs
->index
[icg
+1], &top
->excls
, TRUE
, bexcl
);
2334 ci2xyz(grid
, icg
, &cell_x
, &cell_y
, &cell_z
);
2336 /* Changed iicg to icg, DvdS 990115
2337 * (but see consistency check above, DvdS 990330)
2340 fprintf(log
, "icg=%5d, naaj=%5d, cell %d %d %d\n",
2341 icg
, naaj
, cell_x
, cell_y
, cell_z
);
2343 /* Loop over shift vectors in three dimensions */
2344 for (tz
= -shp
[ZZ
]; tz
<= shp
[ZZ
]; tz
++)
2346 ZI
= cgcm
[icg
][ZZ
]+tz
*box
[ZZ
][ZZ
];
2347 /* Calculate range of cells in Z direction that have the shift tz */
2348 zgi
= cell_z
+ tz
*Nz
;
2351 get_dx(Nz
, gridz
, rl2
, zgi
, ZI
, &dz0
, &dz1
, dcz2
);
2353 get_dx_dd(Nz
, gridz
, rl2
, zgi
, ZI
-grid_offset
[ZZ
],
2354 ncpddc
[ZZ
], sh0
[ZZ
], sh1
[ZZ
], &dz0
, &dz1
, dcz2
);
2360 for (ty
= -shp
[YY
]; ty
<= shp
[YY
]; ty
++)
2362 YI
= cgcm
[icg
][YY
]+ty
*box
[YY
][YY
]+tz
*box
[ZZ
][YY
];
2363 /* Calculate range of cells in Y direction that have the shift ty */
2366 ygi
= (int)(Ny
+ (YI
- grid_offset
[YY
])*grid_y
) - Ny
;
2370 ygi
= cell_y
+ ty
*Ny
;
2373 get_dx(Ny
, gridy
, rl2
, ygi
, YI
, &dy0
, &dy1
, dcy2
);
2375 get_dx_dd(Ny
, gridy
, rl2
, ygi
, YI
-grid_offset
[YY
],
2376 ncpddc
[YY
], sh0
[YY
], sh1
[YY
], &dy0
, &dy1
, dcy2
);
2382 for (tx
= -shp
[XX
]; tx
<= shp
[XX
]; tx
++)
2384 XI
= cgcm
[icg
][XX
]+tx
*box
[XX
][XX
]+ty
*box
[YY
][XX
]+tz
*box
[ZZ
][XX
];
2385 /* Calculate range of cells in X direction that have the shift tx */
2388 xgi
= (int)(Nx
+ (XI
- grid_offset
[XX
])*grid_x
) - Nx
;
2392 xgi
= cell_x
+ tx
*Nx
;
2395 get_dx(Nx
, gridx
, rl2
, xgi
*Nx
, XI
, &dx0
, &dx1
, dcx2
);
2397 get_dx_dd(Nx
, gridx
, rl2
, xgi
, XI
-grid_offset
[XX
],
2398 ncpddc
[XX
], sh0
[XX
], sh1
[XX
], &dx0
, &dx1
, dcx2
);
2404 /* Adress: an explicit cg that has a weigthing function of 0 is excluded
2405 * from the neigbour list as it will not interact */
2406 if (fr
->adress_type
!= eAdressOff
)
2408 if (md
->wf
[cgs
->index
[icg
]] <= GMX_REAL_EPS
&& egp_explicit(fr
, igid
))
2413 /* Get shift vector */
2414 shift
= XYZ2IS(tx
, ty
, tz
);
2416 range_check(shift
, 0, SHIFTS
);
2418 for (nn
= 0; (nn
< ngid
); nn
++)
2425 fprintf(log
, "shift: %2d, dx0,1: %2d,%2d, dy0,1: %2d,%2d, dz0,1: %2d,%2d\n",
2426 shift
, dx0
, dx1
, dy0
, dy1
, dz0
, dz1
);
2427 fprintf(log
, "cgcm: %8.3f %8.3f %8.3f\n", cgcm
[icg
][XX
],
2428 cgcm
[icg
][YY
], cgcm
[icg
][ZZ
]);
2429 fprintf(log
, "xi: %8.3f %8.3f %8.3f\n", XI
, YI
, ZI
);
2431 for (dx
= dx0
; (dx
<= dx1
); dx
++)
2433 tmp1
= rl2
- dcx2
[dx
];
2434 for (dy
= dy0
; (dy
<= dy1
); dy
++)
2436 tmp2
= tmp1
- dcy2
[dy
];
2439 for (dz
= dz0
; (dz
<= dz1
); dz
++)
2441 if (tmp2
> dcz2
[dz
])
2443 /* Find grid-cell cj in which possible neighbours are */
2444 cj
= xyz2ci(Ny
, Nz
, dx
, dy
, dz
);
2446 /* Check out how many cgs (nrj) there in this cell */
2449 /* Find the offset in the cg list */
2452 /* Check if all j's are out of range so we
2453 * can skip the whole cell.
2454 * Should save some time, especially with DD.
2457 (grida
[cgj0
] >= max_jcg
&&
2458 (grida
[cgj0
] >= jcg1
|| grida
[cgj0
+nrj
-1] < jcg0
)))
2464 for (j
= 0; (j
< nrj
); j
++)
2466 jjcg
= grida
[cgj0
+j
];
2468 /* check whether this guy is in range! */
2469 if ((jjcg
>= jcg0
&& jjcg
< jcg1
) ||
2472 r2
= calc_dx2(XI
, YI
, ZI
, cgcm
[jjcg
]);
2475 /* jgid = gid[cgsatoms[cgsindex[jjcg]]]; */
2476 jgid
= GET_CGINFO_GID(cginfo
[jjcg
]);
2477 /* check energy group exclusions */
2478 if (!(i_egp_flags
[jgid
] & EGP_EXCL
))
2482 if (nsr
[jgid
] >= MAX_CG
)
2484 /* Add to short-range list */
2485 put_in_list(bHaveVdW
, ngid
, md
, icg
, jgid
,
2486 nsr
[jgid
], nl_sr
[jgid
],
2487 cgs
->index
, /* cgsatoms, */ bexcl
,
2488 shift
, fr
, FALSE
, TRUE
, TRUE
, fr
->solvent_opt
);
2491 nl_sr
[jgid
][nsr
[jgid
]++] = jjcg
;
2495 if (nlr_ljc
[jgid
] >= MAX_CG
)
2497 /* Add to LJ+coulomb long-range list */
2498 put_in_list(bHaveVdW
, ngid
, md
, icg
, jgid
,
2499 nlr_ljc
[jgid
], nl_lr_ljc
[jgid
], top
->cgs
.index
,
2500 bexcl
, shift
, fr
, TRUE
, TRUE
, TRUE
, fr
->solvent_opt
);
2503 nl_lr_ljc
[jgid
][nlr_ljc
[jgid
]++] = jjcg
;
2507 if (nlr_one
[jgid
] >= MAX_CG
)
2509 /* Add to long-range list with only coul, or only LJ */
2510 put_in_list(bHaveVdW
, ngid
, md
, icg
, jgid
,
2511 nlr_one
[jgid
], nl_lr_one
[jgid
], top
->cgs
.index
,
2512 bexcl
, shift
, fr
, TRUE
, rvdw_lt_rcoul
, rcoul_lt_rvdw
, fr
->solvent_opt
);
2515 nl_lr_one
[jgid
][nlr_one
[jgid
]++] = jjcg
;
2527 /* CHECK whether there is anything left in the buffers */
2528 for (nn
= 0; (nn
< ngid
); nn
++)
2532 put_in_list(bHaveVdW
, ngid
, md
, icg
, nn
, nsr
[nn
], nl_sr
[nn
],
2533 cgs
->index
, /* cgsatoms, */ bexcl
,
2534 shift
, fr
, FALSE
, TRUE
, TRUE
, fr
->solvent_opt
);
2537 if (nlr_ljc
[nn
] > 0)
2539 put_in_list(bHaveVdW
, ngid
, md
, icg
, nn
, nlr_ljc
[nn
],
2540 nl_lr_ljc
[nn
], top
->cgs
.index
,
2541 bexcl
, shift
, fr
, TRUE
, TRUE
, TRUE
, fr
->solvent_opt
);
2544 if (nlr_one
[nn
] > 0)
2546 put_in_list(bHaveVdW
, ngid
, md
, icg
, nn
, nlr_one
[nn
],
2547 nl_lr_one
[nn
], top
->cgs
.index
,
2548 bexcl
, shift
, fr
, TRUE
, rvdw_lt_rcoul
, rcoul_lt_rvdw
, fr
->solvent_opt
);
2554 /* setexcl(nri,i_atoms,&top->atoms.excl,FALSE,bexcl); */
2555 setexcl(cgs
->index
[icg
], cgs
->index
[icg
+1], &top
->excls
, FALSE
, bexcl
);
2557 /* No need to perform any left-over force calculations anymore (as we used to do here)
2558 * since we now save the proper long-range lists for later evaluation.
2563 /* Close neighbourlists */
2564 close_neighbor_lists(fr
, bMakeQMMMnblist
);
2569 void ns_realloc_natoms(gmx_ns_t
*ns
, int natoms
)
2573 if (natoms
> ns
->nra_alloc
)
2575 ns
->nra_alloc
= over_alloc_dd(natoms
);
2576 srenew(ns
->bexcl
, ns
->nra_alloc
);
2577 for (i
= 0; i
< ns
->nra_alloc
; i
++)
2584 void init_ns(FILE *fplog
, const t_commrec
*cr
,
2585 gmx_ns_t
*ns
, t_forcerec
*fr
,
2586 const gmx_mtop_t
*mtop
)
2588 int mt
, icg
, nr_in_cg
, maxcg
, i
, j
, jcg
, ngid
, ncg
;
2592 /* Compute largest charge groups size (# atoms) */
2594 for (mt
= 0; mt
< mtop
->nmoltype
; mt
++)
2596 cgs
= &mtop
->moltype
[mt
].cgs
;
2597 for (icg
= 0; (icg
< cgs
->nr
); icg
++)
2599 nr_in_cg
= max(nr_in_cg
, (int)(cgs
->index
[icg
+1]-cgs
->index
[icg
]));
2603 /* Verify whether largest charge group is <= max cg.
2604 * This is determined by the type of the local exclusion type
2605 * Exclusions are stored in bits. (If the type is not large
2606 * enough, enlarge it, unsigned char -> unsigned short -> unsigned long)
2608 maxcg
= sizeof(t_excl
)*8;
2609 if (nr_in_cg
> maxcg
)
2611 gmx_fatal(FARGS
, "Max #atoms in a charge group: %d > %d\n",
2615 ngid
= mtop
->groups
.grps
[egcENER
].nr
;
2616 snew(ns
->bExcludeAlleg
, ngid
);
2617 for (i
= 0; i
< ngid
; i
++)
2619 ns
->bExcludeAlleg
[i
] = TRUE
;
2620 for (j
= 0; j
< ngid
; j
++)
2622 if (!(fr
->egp_flags
[i
*ngid
+j
] & EGP_EXCL
))
2624 ns
->bExcludeAlleg
[i
] = FALSE
;
2632 ns
->grid
= init_grid(fplog
, fr
);
2633 init_nsgrid_lists(fr
, ngid
, ns
);
2638 snew(ns
->ns_buf
, ngid
);
2639 for (i
= 0; (i
< ngid
); i
++)
2641 snew(ns
->ns_buf
[i
], SHIFTS
);
2643 ncg
= ncg_mtop(mtop
);
2644 snew(ns
->simple_aaj
, 2*ncg
);
2645 for (jcg
= 0; (jcg
< ncg
); jcg
++)
2647 ns
->simple_aaj
[jcg
] = jcg
;
2648 ns
->simple_aaj
[jcg
+ncg
] = jcg
;
2652 /* Create array that determines whether or not atoms have VdW */
2653 snew(ns
->bHaveVdW
, fr
->ntype
);
2654 for (i
= 0; (i
< fr
->ntype
); i
++)
2656 for (j
= 0; (j
< fr
->ntype
); j
++)
2658 ns
->bHaveVdW
[i
] = (ns
->bHaveVdW
[i
] ||
2660 ((BHAMA(fr
->nbfp
, fr
->ntype
, i
, j
) != 0) ||
2661 (BHAMB(fr
->nbfp
, fr
->ntype
, i
, j
) != 0) ||
2662 (BHAMC(fr
->nbfp
, fr
->ntype
, i
, j
) != 0)) :
2663 ((C6(fr
->nbfp
, fr
->ntype
, i
, j
) != 0) ||
2664 (C12(fr
->nbfp
, fr
->ntype
, i
, j
) != 0))));
2669 pr_bvec(debug
, 0, "bHaveVdW", ns
->bHaveVdW
, fr
->ntype
, TRUE
);
2674 if (!DOMAINDECOMP(cr
))
2676 ns_realloc_natoms(ns
, mtop
->natoms
);
2679 ns
->nblist_initialized
= FALSE
;
2681 /* nbr list debug dump */
2683 char *ptr
= getenv("GMX_DUMP_NL");
2686 ns
->dump_nl
= strtol(ptr
, NULL
, 10);
2689 fprintf(fplog
, "GMX_DUMP_NL = %d", ns
->dump_nl
);
2700 int search_neighbours(FILE *log
, t_forcerec
*fr
,
2702 gmx_localtop_t
*top
,
2703 gmx_groups_t
*groups
,
2705 t_nrnb
*nrnb
, t_mdatoms
*md
,
2707 gmx_bool bDoLongRangeNS
)
2709 t_block
*cgs
= &(top
->cgs
);
2710 rvec box_size
, grid_x0
, grid_x1
;
2712 real min_size
, grid_dens
;
2716 gmx_bool
*i_egp_flags
;
2717 int cg_start
, cg_end
, start
, end
;
2720 gmx_domdec_zones_t
*dd_zones
;
2721 put_in_list_t
*put_in_list
;
2725 /* Set some local variables */
2727 ngid
= groups
->grps
[egcENER
].nr
;
2729 for (m
= 0; (m
< DIM
); m
++)
2731 box_size
[m
] = box
[m
][m
];
2734 if (fr
->ePBC
!= epbcNONE
)
2736 if (sqr(fr
->rlistlong
) >= max_cutoff2(fr
->ePBC
, box
))
2738 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.");
2742 min_size
= min(box_size
[XX
], min(box_size
[YY
], box_size
[ZZ
]));
2743 if (2*fr
->rlistlong
>= min_size
)
2745 gmx_fatal(FARGS
, "One of the box diagonal elements has become smaller than twice the cut-off length.");
2750 if (DOMAINDECOMP(cr
))
2752 ns_realloc_natoms(ns
, cgs
->index
[cgs
->nr
]);
2756 /* Reset the neighbourlists */
2757 reset_neighbor_lists(fr
, TRUE
, TRUE
);
2759 if (bGrid
&& bFillGrid
)
2763 if (DOMAINDECOMP(cr
))
2765 dd_zones
= domdec_zones(cr
->dd
);
2771 get_nsgrid_boundaries(grid
->nboundeddim
, box
, NULL
, NULL
, NULL
, NULL
,
2772 cgs
->nr
, fr
->cg_cm
, grid_x0
, grid_x1
, &grid_dens
);
2774 grid_first(log
, grid
, NULL
, NULL
, box
, grid_x0
, grid_x1
,
2775 fr
->rlistlong
, grid_dens
);
2782 if (DOMAINDECOMP(cr
))
2785 fill_grid(dd_zones
, grid
, end
, -1, end
, fr
->cg_cm
);
2787 grid
->icg1
= dd_zones
->izone
[dd_zones
->nizone
-1].cg1
;
2791 fill_grid(NULL
, grid
, cgs
->nr
, fr
->cg0
, fr
->hcg
, fr
->cg_cm
);
2792 grid
->icg0
= fr
->cg0
;
2793 grid
->icg1
= fr
->hcg
;
2797 calc_elemnr(grid
, start
, end
, cgs
->nr
);
2799 grid_last(grid
, start
, end
, cgs
->nr
);
2804 print_grid(debug
, grid
);
2809 /* Set the grid cell index for the test particle only.
2810 * The cell to cg index is not corrected, but that does not matter.
2812 fill_grid(NULL
, ns
->grid
, fr
->hcg
, fr
->hcg
-1, fr
->hcg
, fr
->cg_cm
);
2816 if (fr
->adress_type
== eAdressOff
)
2818 if (!fr
->ns
.bCGlist
)
2820 put_in_list
= put_in_list_at
;
2824 put_in_list
= put_in_list_cg
;
2829 put_in_list
= put_in_list_adress
;
2836 nsearch
= nsgrid_core(cr
, fr
, box
, ngid
, top
,
2837 grid
, ns
->bexcl
, ns
->bExcludeAlleg
,
2838 md
, put_in_list
, ns
->bHaveVdW
,
2839 bDoLongRangeNS
, FALSE
);
2841 /* neighbour searching withouth QMMM! QM atoms have zero charge in
2842 * the classical calculation. The charge-charge interaction
2843 * between QM and MM atoms is handled in the QMMM core calculation
2844 * (see QMMM.c). The VDW however, we'd like to compute classically
2845 * and the QM MM atom pairs have just been put in the
2846 * corresponding neighbourlists. in case of QMMM we still need to
2847 * fill a special QMMM neighbourlist that contains all neighbours
2848 * of the QM atoms. If bQMMM is true, this list will now be made:
2850 if (fr
->bQMMM
&& fr
->qr
->QMMMscheme
!= eQMMMschemeoniom
)
2852 nsearch
+= nsgrid_core(cr
, fr
, box
, ngid
, top
,
2853 grid
, ns
->bexcl
, ns
->bExcludeAlleg
,
2854 md
, put_in_list_qmmm
, ns
->bHaveVdW
,
2855 bDoLongRangeNS
, TRUE
);
2860 nsearch
= ns_simple_core(fr
, top
, md
, box
, box_size
,
2861 ns
->bexcl
, ns
->simple_aaj
,
2862 ngid
, ns
->ns_buf
, put_in_list
, ns
->bHaveVdW
);
2870 inc_nrnb(nrnb
, eNR_NS
, nsearch
);
2871 /* inc_nrnb(nrnb,eNR_LR,fr->nlr); */
2876 int natoms_beyond_ns_buffer(t_inputrec
*ir
, t_forcerec
*fr
, t_block
*cgs
,
2877 matrix scale_tot
, rvec
*x
)
2879 int cg0
, cg1
, cg
, a0
, a1
, a
, i
, j
;
2880 real rint
, hbuf2
, scale
;
2882 gmx_bool bIsotropic
;
2887 rint
= max(ir
->rcoulomb
, ir
->rvdw
);
2888 if (ir
->rlist
< rint
)
2890 gmx_fatal(FARGS
, "The neighbor search buffer has negative size: %f nm",
2898 if (!EI_DYNAMICS(ir
->eI
) || !DYNAMIC_BOX(*ir
))
2900 hbuf2
= sqr(0.5*(ir
->rlist
- rint
));
2901 for (cg
= cg0
; cg
< cg1
; cg
++)
2903 a0
= cgs
->index
[cg
];
2904 a1
= cgs
->index
[cg
+1];
2905 for (a
= a0
; a
< a1
; a
++)
2907 if (distance2(cg_cm
[cg
], x
[a
]) > hbuf2
)
2917 scale
= scale_tot
[0][0];
2918 for (i
= 1; i
< DIM
; i
++)
2920 /* With anisotropic scaling, the original spherical ns volumes become
2921 * ellipsoids. To avoid costly transformations we use the minimum
2922 * eigenvalue of the scaling matrix for determining the buffer size.
2923 * Since the lower half is 0, the eigenvalues are the diagonal elements.
2925 scale
= min(scale
, scale_tot
[i
][i
]);
2926 if (scale_tot
[i
][i
] != scale_tot
[i
-1][i
-1])
2930 for (j
= 0; j
< i
; j
++)
2932 if (scale_tot
[i
][j
] != 0)
2938 hbuf2
= sqr(0.5*(scale
*ir
->rlist
- rint
));
2941 for (cg
= cg0
; cg
< cg1
; cg
++)
2943 svmul(scale
, cg_cm
[cg
], cgsc
);
2944 a0
= cgs
->index
[cg
];
2945 a1
= cgs
->index
[cg
+1];
2946 for (a
= a0
; a
< a1
; a
++)
2948 if (distance2(cgsc
, x
[a
]) > hbuf2
)
2957 /* Anistropic scaling */
2958 for (cg
= cg0
; cg
< cg1
; cg
++)
2960 /* Since scale_tot contains the transpose of the scaling matrix,
2961 * we need to multiply with the transpose.
2963 tmvmul_ur0(scale_tot
, cg_cm
[cg
], cgsc
);
2964 a0
= cgs
->index
[cg
];
2965 a1
= cgs
->index
[cg
+1];
2966 for (a
= a0
; a
< a1
; a
++)
2968 if (distance2(cgsc
, x
[a
]) > hbuf2
)