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