/home/alexxy/Develop/gromacs/src/gromacs/mdlib/ns.c

Bug Summary

File:	gromacs/mdlib/ns.c
Location:	line 1168, column 5
Description:	Value stored to 'chargeB' is never read

Annotated Source Code

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->nWatMol3; /
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->maxnrjMAX_CGCGSIZE)(nlist->excl) = save_realloc("nlist->excl", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/ns.c" , 527, (nlist->excl), (nlist->maxnrj32), 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) ? (igidngid+jgid) : (jgidngid+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) ? (igidngid+jgid) : (jgidngid+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) ? (igidngid+jgid) : (jgidngid+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) ? (igidngid+jgid) : (jgidngid+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;
	Value stored to 'chargeB' is never read
1169	type = md->typeA;
1170	typeB = md->typeB;
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) ? (igidngid+jgid) : (jgidngid+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) ? (igidngid+jgid) : (jgidngid+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) ? (igidngid+jgid) : (jgidngid+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) ? (igidngid+jgid) : (jgidngid+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 = (dxdx) + (dydy) + (dz*dz);
1642
1643	shift = XYZ2IS(tx, ty, tz)((22 +1)((21 +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 = (ddxddx) + (ddyddy) + (ddz*ddz);
1679
1680	shift = XYZ2IS(tx, ty, tz)((22 +1)((21 +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(((21 +1)(21 +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(((21 +1)(21 +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(((21 +1)(21 +1)(2*2 +1))/2)], nrj, cg_j,
1780	bHaveVdW, ngid, md, icg, jgid, cgs, bexcl, CENTRAL(((21 +1)(21 +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 + ngidmd->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((21 +1)(21 +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]+tybox[YY1][YY1]+tzbox[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]+txbox[XX0][XX0]+tybox[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)((22 +1)((21 +1)((tz)+1)+(ty)+1)+(tx)+2);
2416	#ifdef NS5DB
2417	range_check(shift, 0, SHIFTS)_range_check(shift, 0, ((21 +1)(21 +1)(22 +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, (((21 +1)(21 +1)(22 +1))), sizeof((ns->ns_buf [i])));
2643	}
2644	ncg = ncg_mtop(mtop);
2645	snew(ns->simple_aaj, 2ncg)(ns->simple_aaj) = save_calloc("ns->simple_aaj", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/ns.c" , 2645, (2ncg), 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(scaleir->rlist - rint))((0.5(scaleir->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	}