File: | gromacs/mdlib/domdec_con.c |
Location: | line 1153, column 26 |
Description: | Access to field 'ireq' results in a dereference of a null pointer (loaded from field 'constraint_comm') |
1 | /* | |||
2 | * This file is part of the GROMACS molecular simulation package. | |||
3 | * | |||
4 | * Copyright (c) 2006,2007,2008,2009,2010,2012,2013,2014, by the GROMACS development team, led by | |||
5 | * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, | |||
6 | * and including many others, as listed in the AUTHORS file in the | |||
7 | * top-level source directory and at http://www.gromacs.org. | |||
8 | * | |||
9 | * GROMACS is free software; you can redistribute it and/or | |||
10 | * modify it under the terms of the GNU Lesser General Public License | |||
11 | * as published by the Free Software Foundation; either version 2.1 | |||
12 | * of the License, or (at your option) any later version. | |||
13 | * | |||
14 | * GROMACS is distributed in the hope that it will be useful, | |||
15 | * but WITHOUT ANY WARRANTY; without even the implied warranty of | |||
16 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |||
17 | * Lesser General Public License for more details. | |||
18 | * | |||
19 | * You should have received a copy of the GNU Lesser General Public | |||
20 | * License along with GROMACS; if not, see | |||
21 | * http://www.gnu.org/licenses, or write to the Free Software Foundation, | |||
22 | * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. | |||
23 | * | |||
24 | * If you want to redistribute modifications to GROMACS, please | |||
25 | * consider that scientific software is very special. Version | |||
26 | * control is crucial - bugs must be traceable. We will be happy to | |||
27 | * consider code for inclusion in the official distribution, but | |||
28 | * derived work must not be called official GROMACS. Details are found | |||
29 | * in the README & COPYING files - if they are missing, get the | |||
30 | * official version at http://www.gromacs.org. | |||
31 | * | |||
32 | * To help us fund GROMACS development, we humbly ask that you cite | |||
33 | * the research papers on the package. Check out http://www.gromacs.org. | |||
34 | */ | |||
35 | ||||
36 | #ifdef HAVE_CONFIG_H1 | |||
37 | #include <config.h> | |||
38 | #endif | |||
39 | #include <assert.h> | |||
40 | ||||
41 | #include "gromacs/utility/smalloc.h" | |||
42 | #include "gromacs/math/vec.h" | |||
43 | #include "constr.h" | |||
44 | #include "types/commrec.h" | |||
45 | #include "domdec.h" | |||
46 | #include "domdec_network.h" | |||
47 | #include "mtop_util.h" | |||
48 | #include "gmx_ga2la.h" | |||
49 | #include "gmx_hash.h" | |||
50 | #include "gmx_omp_nthreads.h" | |||
51 | #include "macros.h" | |||
52 | ||||
53 | #include "gromacs/utility/fatalerror.h" | |||
54 | ||||
55 | typedef struct { | |||
56 | int nsend; | |||
57 | int *a; | |||
58 | int a_nalloc; | |||
59 | int nrecv; | |||
60 | } gmx_specatsend_t; | |||
61 | ||||
62 | typedef struct { | |||
63 | int *ind; | |||
64 | int nalloc; | |||
65 | int n; | |||
66 | } ind_req_t; | |||
67 | ||||
68 | typedef struct gmx_domdec_specat_comm { | |||
69 | /* The number of indices to receive during the setup */ | |||
70 | int nreq[DIM3][2][2]; | |||
71 | /* The atoms to send */ | |||
72 | gmx_specatsend_t spas[DIM3][2]; | |||
73 | gmx_bool *bSendAtom; | |||
74 | int bSendAtom_nalloc; | |||
75 | /* Send buffers */ | |||
76 | int *ibuf; | |||
77 | int ibuf_nalloc; | |||
78 | rvec *vbuf; | |||
79 | int vbuf_nalloc; | |||
80 | rvec *vbuf2; | |||
81 | int vbuf2_nalloc; | |||
82 | /* The range in the local buffer(s) for received atoms */ | |||
83 | int at_start; | |||
84 | int at_end; | |||
85 | ||||
86 | /* The atom indices we need from the surrounding cells. | |||
87 | * We can gather the indices over nthread threads. | |||
88 | */ | |||
89 | int nthread; | |||
90 | ind_req_t *ireq; | |||
91 | } gmx_domdec_specat_comm_t; | |||
92 | ||||
93 | typedef struct gmx_domdec_constraints { | |||
94 | int *molb_con_offset; | |||
95 | int *molb_ncon_mol; | |||
96 | /* The fully local and connected constraints */ | |||
97 | int ncon; | |||
98 | /* The global constraint number, only required for clearing gc_req */ | |||
99 | int *con_gl; | |||
100 | int *con_nlocat; | |||
101 | int con_nalloc; | |||
102 | /* Boolean that tells if a global constraint index has been requested */ | |||
103 | char *gc_req; | |||
104 | /* Global to local communicated constraint atom only index */ | |||
105 | gmx_hash_t ga2la; | |||
106 | ||||
107 | /* Multi-threading stuff */ | |||
108 | int nthread; | |||
109 | t_ilist *ils; | |||
110 | } gmx_domdec_constraints_t; | |||
111 | ||||
112 | ||||
113 | static void dd_move_f_specat(gmx_domdec_t *dd, gmx_domdec_specat_comm_t *spac, | |||
114 | rvec *f, rvec *fshift) | |||
115 | { | |||
116 | gmx_specatsend_t *spas; | |||
117 | rvec *vbuf; | |||
118 | int n, n0, n1, d, dim, dir, i; | |||
119 | ivec vis; | |||
120 | int is; | |||
121 | gmx_bool bPBC, bScrew; | |||
122 | ||||
123 | n = spac->at_end; | |||
124 | for (d = dd->ndim-1; d >= 0; d--) | |||
125 | { | |||
126 | dim = dd->dim[d]; | |||
127 | if (dd->nc[dim] > 2) | |||
128 | { | |||
129 | /* Pulse the grid forward and backward */ | |||
130 | spas = spac->spas[d]; | |||
131 | n0 = spas[0].nrecv; | |||
132 | n1 = spas[1].nrecv; | |||
133 | n -= n1 + n0; | |||
134 | vbuf = spac->vbuf; | |||
135 | /* Send and receive the coordinates */ | |||
136 | dd_sendrecv2_rvec(dd, d, | |||
137 | f+n+n1, n0, vbuf, spas[0].nsend, | |||
138 | f+n, n1, vbuf+spas[0].nsend, spas[1].nsend); | |||
139 | for (dir = 0; dir < 2; dir++) | |||
140 | { | |||
141 | bPBC = ((dir == 0 && dd->ci[dim] == 0) || | |||
142 | (dir == 1 && dd->ci[dim] == dd->nc[dim]-1)); | |||
143 | bScrew = (bPBC && dd->bScrewPBC && dim == XX0); | |||
144 | ||||
145 | spas = &spac->spas[d][dir]; | |||
146 | /* Sum the buffer into the required forces */ | |||
147 | if (!bPBC || (!bScrew && fshift == NULL((void*)0))) | |||
148 | { | |||
149 | for (i = 0; i < spas->nsend; i++) | |||
150 | { | |||
151 | rvec_inc(f[spas->a[i]], *vbuf); | |||
152 | vbuf++; | |||
153 | } | |||
154 | } | |||
155 | else | |||
156 | { | |||
157 | clear_ivec(vis); | |||
158 | vis[dim] = (dir == 0 ? 1 : -1); | |||
159 | is = IVEC2IS(vis)(((2*2 +1)*((2*1 +1)*(((vis)[2])+1)+((vis)[1])+1)+((vis)[0])+ 2)); | |||
160 | if (!bScrew) | |||
161 | { | |||
162 | /* Sum and add to shift forces */ | |||
163 | for (i = 0; i < spas->nsend; i++) | |||
164 | { | |||
165 | rvec_inc(f[spas->a[i]], *vbuf); | |||
166 | rvec_inc(fshift[is], *vbuf); | |||
167 | vbuf++; | |||
168 | } | |||
169 | } | |||
170 | else | |||
171 | { | |||
172 | /* Rotate the forces */ | |||
173 | for (i = 0; i < spas->nsend; i++) | |||
174 | { | |||
175 | f[spas->a[i]][XX0] += (*vbuf)[XX0]; | |||
176 | f[spas->a[i]][YY1] -= (*vbuf)[YY1]; | |||
177 | f[spas->a[i]][ZZ2] -= (*vbuf)[ZZ2]; | |||
178 | if (fshift) | |||
179 | { | |||
180 | rvec_inc(fshift[is], *vbuf); | |||
181 | } | |||
182 | vbuf++; | |||
183 | } | |||
184 | } | |||
185 | } | |||
186 | } | |||
187 | } | |||
188 | else | |||
189 | { | |||
190 | /* Two cells, so we only need to communicate one way */ | |||
191 | spas = &spac->spas[d][0]; | |||
192 | n -= spas->nrecv; | |||
193 | /* Send and receive the coordinates */ | |||
194 | dd_sendrecv_rvec(dd, d, dddirForward, | |||
195 | f+n, spas->nrecv, spac->vbuf, spas->nsend); | |||
196 | /* Sum the buffer into the required forces */ | |||
197 | if (dd->bScrewPBC && dim == XX0 && | |||
198 | (dd->ci[dim] == 0 || | |||
199 | dd->ci[dim] == dd->nc[dim]-1)) | |||
200 | { | |||
201 | for (i = 0; i < spas->nsend; i++) | |||
202 | { | |||
203 | /* Rotate the force */ | |||
204 | f[spas->a[i]][XX0] += spac->vbuf[i][XX0]; | |||
205 | f[spas->a[i]][YY1] -= spac->vbuf[i][YY1]; | |||
206 | f[spas->a[i]][ZZ2] -= spac->vbuf[i][ZZ2]; | |||
207 | } | |||
208 | } | |||
209 | else | |||
210 | { | |||
211 | for (i = 0; i < spas->nsend; i++) | |||
212 | { | |||
213 | rvec_inc(f[spas->a[i]], spac->vbuf[i]); | |||
214 | } | |||
215 | } | |||
216 | } | |||
217 | } | |||
218 | } | |||
219 | ||||
220 | void dd_move_f_vsites(gmx_domdec_t *dd, rvec *f, rvec *fshift) | |||
221 | { | |||
222 | if (dd->vsite_comm) | |||
223 | { | |||
224 | dd_move_f_specat(dd, dd->vsite_comm, f, fshift); | |||
225 | } | |||
226 | } | |||
227 | ||||
228 | void dd_clear_f_vsites(gmx_domdec_t *dd, rvec *f) | |||
229 | { | |||
230 | int i; | |||
231 | ||||
232 | if (dd->vsite_comm) | |||
233 | { | |||
234 | for (i = dd->vsite_comm->at_start; i < dd->vsite_comm->at_end; i++) | |||
235 | { | |||
236 | clear_rvec(f[i]); | |||
237 | } | |||
238 | } | |||
239 | } | |||
240 | ||||
241 | static void dd_move_x_specat(gmx_domdec_t *dd, gmx_domdec_specat_comm_t *spac, | |||
242 | matrix box, | |||
243 | rvec *x0, | |||
244 | rvec *x1, gmx_bool bX1IsCoord) | |||
245 | { | |||
246 | gmx_specatsend_t *spas; | |||
247 | rvec *x, *vbuf, *rbuf; | |||
248 | int nvec, v, n, nn, ns0, ns1, nr0, nr1, nr, d, dim, dir, i; | |||
249 | gmx_bool bPBC, bScrew = FALSE0; | |||
250 | rvec shift = {0, 0, 0}; | |||
251 | ||||
252 | nvec = 1; | |||
253 | if (x1 != NULL((void*)0)) | |||
254 | { | |||
255 | nvec++; | |||
256 | } | |||
257 | ||||
258 | n = spac->at_start; | |||
259 | for (d = 0; d < dd->ndim; d++) | |||
260 | { | |||
261 | dim = dd->dim[d]; | |||
262 | if (dd->nc[dim] > 2) | |||
263 | { | |||
264 | /* Pulse the grid forward and backward */ | |||
265 | vbuf = spac->vbuf; | |||
266 | for (dir = 0; dir < 2; dir++) | |||
267 | { | |||
268 | if (dir == 0 && dd->ci[dim] == 0) | |||
269 | { | |||
270 | bPBC = TRUE1; | |||
271 | bScrew = (dd->bScrewPBC && dim == XX0); | |||
272 | copy_rvec(box[dim], shift); | |||
273 | } | |||
274 | else if (dir == 1 && dd->ci[dim] == dd->nc[dim]-1) | |||
275 | { | |||
276 | bPBC = TRUE1; | |||
277 | bScrew = (dd->bScrewPBC && dim == XX0); | |||
278 | for (i = 0; i < DIM3; i++) | |||
279 | { | |||
280 | shift[i] = -box[dim][i]; | |||
281 | } | |||
282 | } | |||
283 | else | |||
284 | { | |||
285 | bPBC = FALSE0; | |||
286 | bScrew = FALSE0; | |||
287 | } | |||
288 | spas = &spac->spas[d][dir]; | |||
289 | for (v = 0; v < nvec; v++) | |||
290 | { | |||
291 | x = (v == 0 ? x0 : x1); | |||
292 | /* Copy the required coordinates to the send buffer */ | |||
293 | if (!bPBC || (v == 1 && !bX1IsCoord)) | |||
294 | { | |||
295 | /* Only copy */ | |||
296 | for (i = 0; i < spas->nsend; i++) | |||
297 | { | |||
298 | copy_rvec(x[spas->a[i]], *vbuf); | |||
299 | vbuf++; | |||
300 | } | |||
301 | } | |||
302 | else if (!bScrew) | |||
303 | { | |||
304 | /* Shift coordinates */ | |||
305 | for (i = 0; i < spas->nsend; i++) | |||
306 | { | |||
307 | rvec_add(x[spas->a[i]], shift, *vbuf); | |||
308 | vbuf++; | |||
309 | } | |||
310 | } | |||
311 | else | |||
312 | { | |||
313 | /* Shift and rotate coordinates */ | |||
314 | for (i = 0; i < spas->nsend; i++) | |||
315 | { | |||
316 | (*vbuf)[XX0] = x[spas->a[i]][XX0] + shift[XX0]; | |||
317 | (*vbuf)[YY1] = box[YY1][YY1] - x[spas->a[i]][YY1] + shift[YY1]; | |||
318 | (*vbuf)[ZZ2] = box[ZZ2][ZZ2] - x[spas->a[i]][ZZ2] + shift[ZZ2]; | |||
319 | vbuf++; | |||
320 | } | |||
321 | } | |||
322 | } | |||
323 | } | |||
324 | /* Send and receive the coordinates */ | |||
325 | spas = spac->spas[d]; | |||
326 | ns0 = spas[0].nsend; | |||
327 | nr0 = spas[0].nrecv; | |||
328 | ns1 = spas[1].nsend; | |||
329 | nr1 = spas[1].nrecv; | |||
330 | if (nvec == 1) | |||
331 | { | |||
332 | dd_sendrecv2_rvec(dd, d, | |||
333 | spac->vbuf+ns0, ns1, x0+n, nr1, | |||
334 | spac->vbuf, ns0, x0+n+nr1, nr0); | |||
335 | } | |||
336 | else | |||
337 | { | |||
338 | /* Communicate both vectors in one buffer */ | |||
339 | rbuf = spac->vbuf2; | |||
340 | dd_sendrecv2_rvec(dd, d, | |||
341 | spac->vbuf+2*ns0, 2*ns1, rbuf, 2*nr1, | |||
342 | spac->vbuf, 2*ns0, rbuf+2*nr1, 2*nr0); | |||
343 | /* Split the buffer into the two vectors */ | |||
344 | nn = n; | |||
345 | for (dir = 1; dir >= 0; dir--) | |||
346 | { | |||
347 | nr = spas[dir].nrecv; | |||
348 | for (v = 0; v < 2; v++) | |||
349 | { | |||
350 | x = (v == 0 ? x0 : x1); | |||
351 | for (i = 0; i < nr; i++) | |||
352 | { | |||
353 | copy_rvec(*rbuf, x[nn+i]); | |||
354 | rbuf++; | |||
355 | } | |||
356 | } | |||
357 | nn += nr; | |||
358 | } | |||
359 | } | |||
360 | n += nr0 + nr1; | |||
361 | } | |||
362 | else | |||
363 | { | |||
364 | spas = &spac->spas[d][0]; | |||
365 | /* Copy the required coordinates to the send buffer */ | |||
366 | vbuf = spac->vbuf; | |||
367 | for (v = 0; v < nvec; v++) | |||
368 | { | |||
369 | x = (v == 0 ? x0 : x1); | |||
370 | if (dd->bScrewPBC && dim == XX0 && | |||
371 | (dd->ci[XX0] == 0 || dd->ci[XX0] == dd->nc[XX0]-1)) | |||
372 | { | |||
373 | /* Here we only perform the rotation, the rest of the pbc | |||
374 | * is handled in the constraint or viste routines. | |||
375 | */ | |||
376 | for (i = 0; i < spas->nsend; i++) | |||
377 | { | |||
378 | (*vbuf)[XX0] = x[spas->a[i]][XX0]; | |||
379 | (*vbuf)[YY1] = box[YY1][YY1] - x[spas->a[i]][YY1]; | |||
380 | (*vbuf)[ZZ2] = box[ZZ2][ZZ2] - x[spas->a[i]][ZZ2]; | |||
381 | vbuf++; | |||
382 | } | |||
383 | } | |||
384 | else | |||
385 | { | |||
386 | for (i = 0; i < spas->nsend; i++) | |||
387 | { | |||
388 | copy_rvec(x[spas->a[i]], *vbuf); | |||
389 | vbuf++; | |||
390 | } | |||
391 | } | |||
392 | } | |||
393 | /* Send and receive the coordinates */ | |||
394 | if (nvec == 1) | |||
395 | { | |||
396 | dd_sendrecv_rvec(dd, d, dddirBackward, | |||
397 | spac->vbuf, spas->nsend, x0+n, spas->nrecv); | |||
398 | } | |||
399 | else | |||
400 | { | |||
401 | /* Communicate both vectors in one buffer */ | |||
402 | rbuf = spac->vbuf2; | |||
403 | dd_sendrecv_rvec(dd, d, dddirBackward, | |||
404 | spac->vbuf, 2*spas->nsend, rbuf, 2*spas->nrecv); | |||
405 | /* Split the buffer into the two vectors */ | |||
406 | nr = spas[0].nrecv; | |||
407 | for (v = 0; v < 2; v++) | |||
408 | { | |||
409 | x = (v == 0 ? x0 : x1); | |||
410 | for (i = 0; i < nr; i++) | |||
411 | { | |||
412 | copy_rvec(*rbuf, x[n+i]); | |||
413 | rbuf++; | |||
414 | } | |||
415 | } | |||
416 | } | |||
417 | n += spas->nrecv; | |||
418 | } | |||
419 | } | |||
420 | } | |||
421 | ||||
422 | void dd_move_x_constraints(gmx_domdec_t *dd, matrix box, | |||
423 | rvec *x0, rvec *x1, gmx_bool bX1IsCoord) | |||
424 | { | |||
425 | if (dd->constraint_comm) | |||
426 | { | |||
427 | dd_move_x_specat(dd, dd->constraint_comm, box, x0, x1, bX1IsCoord); | |||
428 | } | |||
429 | } | |||
430 | ||||
431 | void dd_move_x_vsites(gmx_domdec_t *dd, matrix box, rvec *x) | |||
432 | { | |||
433 | if (dd->vsite_comm) | |||
434 | { | |||
435 | dd_move_x_specat(dd, dd->vsite_comm, box, x, NULL((void*)0), FALSE0); | |||
436 | } | |||
437 | } | |||
438 | ||||
439 | int *dd_constraints_nlocalatoms(gmx_domdec_t *dd) | |||
440 | { | |||
441 | if (dd->constraints) | |||
442 | { | |||
443 | return dd->constraints->con_nlocat; | |||
444 | } | |||
445 | else | |||
446 | { | |||
447 | return NULL((void*)0); | |||
448 | } | |||
449 | } | |||
450 | ||||
451 | void dd_clear_local_constraint_indices(gmx_domdec_t *dd) | |||
452 | { | |||
453 | gmx_domdec_constraints_t *dc; | |||
454 | int i; | |||
455 | ||||
456 | dc = dd->constraints; | |||
457 | ||||
458 | for (i = 0; i < dc->ncon; i++) | |||
459 | { | |||
460 | dc->gc_req[dc->con_gl[i]] = 0; | |||
461 | } | |||
462 | ||||
463 | if (dd->constraint_comm) | |||
464 | { | |||
465 | gmx_hash_clear_and_optimize(dc->ga2la); | |||
466 | } | |||
467 | } | |||
468 | ||||
469 | void dd_clear_local_vsite_indices(gmx_domdec_t *dd) | |||
470 | { | |||
471 | int i; | |||
472 | ||||
473 | if (dd->vsite_comm) | |||
474 | { | |||
475 | gmx_hash_clear_and_optimize(dd->ga2la_vsite); | |||
476 | } | |||
477 | } | |||
478 | ||||
479 | static int setup_specat_communication(gmx_domdec_t *dd, | |||
480 | ind_req_t *ireq, | |||
481 | gmx_domdec_specat_comm_t *spac, | |||
482 | gmx_hash_t ga2la_specat, | |||
483 | int at_start, | |||
484 | int vbuf_fac, | |||
485 | const char *specat_type, | |||
486 | const char *add_err) | |||
487 | { | |||
488 | int nsend[2], nlast, nsend_zero[2] = {0, 0}, *nsend_ptr; | |||
489 | int d, dim, ndir, dir, nr, ns, i, nrecv_local, n0, start, indr, ind, buf[2]; | |||
490 | int nat_tot_specat, nat_tot_prev, nalloc_old; | |||
491 | gmx_bool bPBC, bFirst; | |||
492 | gmx_specatsend_t *spas; | |||
493 | ||||
494 | if (debug) | |||
495 | { | |||
496 | fprintf(debug, "Begin setup_specat_communication for %s\n", specat_type); | |||
497 | } | |||
498 | ||||
499 | /* nsend[0]: the number of atoms requested by this node only, | |||
500 | * we communicate this for more efficients checks | |||
501 | * nsend[1]: the total number of requested atoms | |||
502 | */ | |||
503 | nsend[0] = ireq->n; | |||
504 | nsend[1] = nsend[0]; | |||
505 | nlast = nsend[1]; | |||
506 | for (d = dd->ndim-1; d >= 0; d--) | |||
507 | { | |||
508 | /* Pulse the grid forward and backward */ | |||
509 | dim = dd->dim[d]; | |||
510 | bPBC = (dim < dd->npbcdim); | |||
511 | if (dd->nc[dim] == 2) | |||
512 | { | |||
513 | /* Only 2 cells, so we only need to communicate once */ | |||
514 | ndir = 1; | |||
515 | } | |||
516 | else | |||
517 | { | |||
518 | ndir = 2; | |||
519 | } | |||
520 | for (dir = 0; dir < ndir; dir++) | |||
521 | { | |||
522 | if (!bPBC && | |||
523 | dd->nc[dim] > 2 && | |||
524 | ((dir == 0 && dd->ci[dim] == dd->nc[dim] - 1) || | |||
525 | (dir == 1 && dd->ci[dim] == 0))) | |||
526 | { | |||
527 | /* No pbc: the fist/last cell should not request atoms */ | |||
528 | nsend_ptr = nsend_zero; | |||
529 | } | |||
530 | else | |||
531 | { | |||
532 | nsend_ptr = nsend; | |||
533 | } | |||
534 | /* Communicate the number of indices */ | |||
535 | dd_sendrecv_int(dd, d, dir == 0 ? dddirForward : dddirBackward, | |||
536 | nsend_ptr, 2, spac->nreq[d][dir], 2); | |||
537 | nr = spac->nreq[d][dir][1]; | |||
538 | if (nlast+nr > ireq->nalloc) | |||
539 | { | |||
540 | ireq->nalloc = over_alloc_dd(nlast+nr); | |||
541 | srenew(ireq->ind, ireq->nalloc)(ireq->ind) = save_realloc("ireq->ind", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/domdec_con.c" , 541, (ireq->ind), (ireq->nalloc), sizeof(*(ireq->ind ))); | |||
542 | } | |||
543 | /* Communicate the indices */ | |||
544 | dd_sendrecv_int(dd, d, dir == 0 ? dddirForward : dddirBackward, | |||
545 | ireq->ind, nsend_ptr[1], ireq->ind+nlast, nr); | |||
546 | nlast += nr; | |||
547 | } | |||
548 | nsend[1] = nlast; | |||
549 | } | |||
550 | if (debug) | |||
551 | { | |||
552 | fprintf(debug, "Communicated the counts\n"); | |||
553 | } | |||
554 | ||||
555 | /* Search for the requested atoms and communicate the indices we have */ | |||
556 | nat_tot_specat = at_start; | |||
557 | nrecv_local = 0; | |||
558 | for (d = 0; d < dd->ndim; d++) | |||
559 | { | |||
560 | bFirst = (d == 0); | |||
561 | /* Pulse the grid forward and backward */ | |||
562 | if (dd->dim[d] >= dd->npbcdim || dd->nc[dd->dim[d]] > 2) | |||
563 | { | |||
564 | ndir = 2; | |||
565 | } | |||
566 | else | |||
567 | { | |||
568 | ndir = 1; | |||
569 | } | |||
570 | nat_tot_prev = nat_tot_specat; | |||
571 | for (dir = ndir-1; dir >= 0; dir--) | |||
572 | { | |||
573 | if (nat_tot_specat > spac->bSendAtom_nalloc) | |||
574 | { | |||
575 | nalloc_old = spac->bSendAtom_nalloc; | |||
576 | spac->bSendAtom_nalloc = over_alloc_dd(nat_tot_specat); | |||
577 | srenew(spac->bSendAtom, spac->bSendAtom_nalloc)(spac->bSendAtom) = save_realloc("spac->bSendAtom", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/domdec_con.c" , 577, (spac->bSendAtom), (spac->bSendAtom_nalloc), sizeof (*(spac->bSendAtom))); | |||
578 | for (i = nalloc_old; i < spac->bSendAtom_nalloc; i++) | |||
579 | { | |||
580 | spac->bSendAtom[i] = FALSE0; | |||
581 | } | |||
582 | } | |||
583 | spas = &spac->spas[d][dir]; | |||
584 | n0 = spac->nreq[d][dir][0]; | |||
585 | nr = spac->nreq[d][dir][1]; | |||
586 | if (debug) | |||
587 | { | |||
588 | fprintf(debug, "dim=%d, dir=%d, searching for %d atoms\n", | |||
589 | d, dir, nr); | |||
590 | } | |||
591 | start = nlast - nr; | |||
592 | spas->nsend = 0; | |||
593 | nsend[0] = 0; | |||
594 | for (i = 0; i < nr; i++) | |||
595 | { | |||
596 | indr = ireq->ind[start+i]; | |||
597 | ind = -1; | |||
598 | /* Check if this is a home atom and if so ind will be set */ | |||
599 | if (!ga2la_get_home(dd->ga2la, indr, &ind)) | |||
600 | { | |||
601 | /* Search in the communicated atoms */ | |||
602 | ind = gmx_hash_get_minone(ga2la_specat, indr); | |||
603 | } | |||
604 | if (ind >= 0) | |||
605 | { | |||
606 | if (i < n0 || !spac->bSendAtom[ind]) | |||
607 | { | |||
608 | if (spas->nsend+1 > spas->a_nalloc) | |||
609 | { | |||
610 | spas->a_nalloc = over_alloc_large(spas->nsend+1)(int)(1.19*(spas->nsend+1) + 1000); | |||
611 | srenew(spas->a, spas->a_nalloc)(spas->a) = save_realloc("spas->a", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/domdec_con.c" , 611, (spas->a), (spas->a_nalloc), sizeof(*(spas->a ))); | |||
612 | } | |||
613 | /* Store the local index so we know which coordinates | |||
614 | * to send out later. | |||
615 | */ | |||
616 | spas->a[spas->nsend] = ind; | |||
617 | spac->bSendAtom[ind] = TRUE1; | |||
618 | if (spas->nsend+1 > spac->ibuf_nalloc) | |||
619 | { | |||
620 | spac->ibuf_nalloc = over_alloc_large(spas->nsend+1)(int)(1.19*(spas->nsend+1) + 1000); | |||
621 | srenew(spac->ibuf, spac->ibuf_nalloc)(spac->ibuf) = save_realloc("spac->ibuf", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/domdec_con.c" , 621, (spac->ibuf), (spac->ibuf_nalloc), sizeof(*(spac ->ibuf))); | |||
622 | } | |||
623 | /* Store the global index so we can send it now */ | |||
624 | spac->ibuf[spas->nsend] = indr; | |||
625 | if (i < n0) | |||
626 | { | |||
627 | nsend[0]++; | |||
628 | } | |||
629 | spas->nsend++; | |||
630 | } | |||
631 | } | |||
632 | } | |||
633 | nlast = start; | |||
634 | /* Clear the local flags */ | |||
635 | for (i = 0; i < spas->nsend; i++) | |||
636 | { | |||
637 | spac->bSendAtom[spas->a[i]] = FALSE0; | |||
638 | } | |||
639 | /* Send and receive the number of indices to communicate */ | |||
640 | nsend[1] = spas->nsend; | |||
641 | dd_sendrecv_int(dd, d, dir == 0 ? dddirBackward : dddirForward, | |||
642 | nsend, 2, buf, 2); | |||
643 | if (debug) | |||
644 | { | |||
645 | fprintf(debug, "Send to node %d, %d (%d) indices, " | |||
646 | "receive from node %d, %d (%d) indices\n", | |||
647 | dd->neighbor[d][1-dir], nsend[1], nsend[0], | |||
648 | dd->neighbor[d][dir], buf[1], buf[0]); | |||
649 | if (gmx_debug_at) | |||
650 | { | |||
651 | for (i = 0; i < spas->nsend; i++) | |||
652 | { | |||
653 | fprintf(debug, " %d", spac->ibuf[i]+1); | |||
654 | } | |||
655 | fprintf(debug, "\n"); | |||
656 | } | |||
657 | } | |||
658 | nrecv_local += buf[0]; | |||
659 | spas->nrecv = buf[1]; | |||
660 | if (nat_tot_specat + spas->nrecv > dd->gatindex_nalloc) | |||
661 | { | |||
662 | dd->gatindex_nalloc = | |||
663 | over_alloc_dd(nat_tot_specat + spas->nrecv); | |||
664 | srenew(dd->gatindex, dd->gatindex_nalloc)(dd->gatindex) = save_realloc("dd->gatindex", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/domdec_con.c" , 664, (dd->gatindex), (dd->gatindex_nalloc), sizeof(*( dd->gatindex))); | |||
665 | } | |||
666 | /* Send and receive the indices */ | |||
667 | dd_sendrecv_int(dd, d, dir == 0 ? dddirBackward : dddirForward, | |||
668 | spac->ibuf, spas->nsend, | |||
669 | dd->gatindex+nat_tot_specat, spas->nrecv); | |||
670 | nat_tot_specat += spas->nrecv; | |||
671 | } | |||
672 | ||||
673 | /* Allocate the x/f communication buffers */ | |||
674 | ns = spac->spas[d][0].nsend; | |||
675 | nr = spac->spas[d][0].nrecv; | |||
676 | if (ndir == 2) | |||
677 | { | |||
678 | ns += spac->spas[d][1].nsend; | |||
679 | nr += spac->spas[d][1].nrecv; | |||
680 | } | |||
681 | if (vbuf_fac*ns > spac->vbuf_nalloc) | |||
682 | { | |||
683 | spac->vbuf_nalloc = over_alloc_dd(vbuf_fac*ns); | |||
684 | srenew(spac->vbuf, spac->vbuf_nalloc)(spac->vbuf) = save_realloc("spac->vbuf", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/domdec_con.c" , 684, (spac->vbuf), (spac->vbuf_nalloc), sizeof(*(spac ->vbuf))); | |||
685 | } | |||
686 | if (vbuf_fac == 2 && vbuf_fac*nr > spac->vbuf2_nalloc) | |||
687 | { | |||
688 | spac->vbuf2_nalloc = over_alloc_dd(vbuf_fac*nr); | |||
689 | srenew(spac->vbuf2, spac->vbuf2_nalloc)(spac->vbuf2) = save_realloc("spac->vbuf2", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/domdec_con.c" , 689, (spac->vbuf2), (spac->vbuf2_nalloc), sizeof(*(spac ->vbuf2))); | |||
690 | } | |||
691 | ||||
692 | /* Make a global to local index for the communication atoms */ | |||
693 | for (i = nat_tot_prev; i < nat_tot_specat; i++) | |||
694 | { | |||
695 | gmx_hash_change_or_set(ga2la_specat, dd->gatindex[i], i); | |||
696 | } | |||
697 | } | |||
698 | ||||
699 | /* Check that in the end we got the number of atoms we asked for */ | |||
700 | if (nrecv_local != ireq->n) | |||
701 | { | |||
702 | if (debug) | |||
703 | { | |||
704 | fprintf(debug, "Requested %d, received %d (tot recv %d)\n", | |||
705 | ireq->n, nrecv_local, nat_tot_specat-at_start); | |||
706 | if (gmx_debug_at) | |||
707 | { | |||
708 | for (i = 0; i < ireq->n; i++) | |||
709 | { | |||
710 | ind = gmx_hash_get_minone(ga2la_specat, ireq->ind[i]); | |||
711 | fprintf(debug, " %s%d", | |||
712 | (ind >= 0) ? "" : "!", | |||
713 | ireq->ind[i]+1); | |||
714 | } | |||
715 | fprintf(debug, "\n"); | |||
716 | } | |||
717 | } | |||
718 | fprintf(stderrstderr, "\nDD cell %d %d %d: Neighboring cells do not have atoms:", | |||
719 | dd->ci[XX0], dd->ci[YY1], dd->ci[ZZ2]); | |||
720 | for (i = 0; i < ireq->n; i++) | |||
721 | { | |||
722 | if (gmx_hash_get_minone(ga2la_specat, ireq->ind[i]) < 0) | |||
723 | { | |||
724 | fprintf(stderrstderr, " %d", ireq->ind[i]+1); | |||
725 | } | |||
726 | } | |||
727 | fprintf(stderrstderr, "\n"); | |||
728 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/domdec_con.c" , 728, "DD cell %d %d %d could only obtain %d of the %d atoms that are connected via %ss from the neighboring cells. This probably means your %s lengths are too long compared to the domain decomposition cell size. Decrease the number of domain decomposition grid cells%s%s.", | |||
729 | dd->ci[XX0], dd->ci[YY1], dd->ci[ZZ2], | |||
730 | nrecv_local, ireq->n, specat_type, | |||
731 | specat_type, add_err, | |||
732 | dd->bGridJump ? " or use the -rcon option of mdrun" : ""); | |||
733 | } | |||
734 | ||||
735 | spac->at_start = at_start; | |||
736 | spac->at_end = nat_tot_specat; | |||
737 | ||||
738 | if (debug) | |||
739 | { | |||
740 | fprintf(debug, "Done setup_specat_communication\n"); | |||
741 | } | |||
742 | ||||
743 | return nat_tot_specat; | |||
744 | } | |||
745 | ||||
746 | static void walk_out(int con, int con_offset, int a, int offset, int nrec, | |||
747 | int ncon1, const t_iatom *ia1, const t_iatom *ia2, | |||
748 | const t_blocka *at2con, | |||
749 | const gmx_ga2la_t ga2la, gmx_bool bHomeConnect, | |||
750 | gmx_domdec_constraints_t *dc, | |||
751 | gmx_domdec_specat_comm_t *dcc, | |||
752 | t_ilist *il_local, | |||
753 | ind_req_t *ireq) | |||
754 | { | |||
755 | int a1_gl, a2_gl, a_loc, i, coni, b; | |||
756 | const t_iatom *iap; | |||
757 | ||||
758 | if (dc->gc_req[con_offset+con] == 0) | |||
759 | { | |||
760 | /* Add this non-home constraint to the list */ | |||
761 | if (dc->ncon+1 > dc->con_nalloc) | |||
762 | { | |||
763 | dc->con_nalloc = over_alloc_large(dc->ncon+1)(int)(1.19*(dc->ncon+1) + 1000); | |||
764 | srenew(dc->con_gl, dc->con_nalloc)(dc->con_gl) = save_realloc("dc->con_gl", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/domdec_con.c" , 764, (dc->con_gl), (dc->con_nalloc), sizeof(*(dc-> con_gl))); | |||
765 | srenew(dc->con_nlocat, dc->con_nalloc)(dc->con_nlocat) = save_realloc("dc->con_nlocat", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/domdec_con.c" , 765, (dc->con_nlocat), (dc->con_nalloc), sizeof(*(dc-> con_nlocat))); | |||
766 | } | |||
767 | dc->con_gl[dc->ncon] = con_offset + con; | |||
768 | dc->con_nlocat[dc->ncon] = (bHomeConnect ? 1 : 0); | |||
769 | dc->gc_req[con_offset+con] = 1; | |||
770 | if (il_local->nr + 3 > il_local->nalloc) | |||
771 | { | |||
772 | il_local->nalloc = over_alloc_dd(il_local->nr+3); | |||
773 | srenew(il_local->iatoms, il_local->nalloc)(il_local->iatoms) = save_realloc("il_local->iatoms", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/domdec_con.c" , 773, (il_local->iatoms), (il_local->nalloc), sizeof(* (il_local->iatoms))); | |||
774 | } | |||
775 | iap = constr_iatomptr(ncon1, ia1, ia2, con)((con) < (ncon1) ? (ia1)+(con)*3 : (ia2)+(con-ncon1)*3); | |||
776 | il_local->iatoms[il_local->nr++] = iap[0]; | |||
777 | a1_gl = offset + iap[1]; | |||
778 | a2_gl = offset + iap[2]; | |||
779 | /* The following indexing code can probably be optizimed */ | |||
780 | if (ga2la_get_home(ga2la, a1_gl, &a_loc)) | |||
781 | { | |||
782 | il_local->iatoms[il_local->nr++] = a_loc; | |||
783 | } | |||
784 | else | |||
785 | { | |||
786 | /* We set this index later */ | |||
787 | il_local->iatoms[il_local->nr++] = -a1_gl - 1; | |||
788 | } | |||
789 | if (ga2la_get_home(ga2la, a2_gl, &a_loc)) | |||
790 | { | |||
791 | il_local->iatoms[il_local->nr++] = a_loc; | |||
792 | } | |||
793 | else | |||
794 | { | |||
795 | /* We set this index later */ | |||
796 | il_local->iatoms[il_local->nr++] = -a2_gl - 1; | |||
797 | } | |||
798 | dc->ncon++; | |||
799 | } | |||
800 | /* Check to not ask for the same atom more than once */ | |||
801 | if (gmx_hash_get_minone(dc->ga2la, offset+a) == -1) | |||
802 | { | |||
803 | assert(dcc)((void) (0)); | |||
804 | /* Add this non-home atom to the list */ | |||
805 | if (ireq->n+1 > ireq->nalloc) | |||
806 | { | |||
807 | ireq->nalloc = over_alloc_large(ireq->n+1)(int)(1.19*(ireq->n+1) + 1000); | |||
808 | srenew(ireq->ind, ireq->nalloc)(ireq->ind) = save_realloc("ireq->ind", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/domdec_con.c" , 808, (ireq->ind), (ireq->nalloc), sizeof(*(ireq->ind ))); | |||
809 | } | |||
810 | ireq->ind[ireq->n++] = offset + a; | |||
811 | /* Temporarily mark with -2, we get the index later */ | |||
812 | gmx_hash_set(dc->ga2la, offset+a, -2); | |||
813 | } | |||
814 | ||||
815 | if (nrec > 0) | |||
816 | { | |||
817 | for (i = at2con->index[a]; i < at2con->index[a+1]; i++) | |||
818 | { | |||
819 | coni = at2con->a[i]; | |||
820 | if (coni != con) | |||
821 | { | |||
822 | /* Walk further */ | |||
823 | iap = constr_iatomptr(ncon1, ia1, ia2, coni)((coni) < (ncon1) ? (ia1)+(coni)*3 : (ia2)+(coni-ncon1)*3); | |||
824 | if (a == iap[1]) | |||
825 | { | |||
826 | b = iap[2]; | |||
827 | } | |||
828 | else | |||
829 | { | |||
830 | b = iap[1]; | |||
831 | } | |||
832 | if (!ga2la_get_home(ga2la, offset+b, &a_loc)) | |||
833 | { | |||
834 | walk_out(coni, con_offset, b, offset, nrec-1, | |||
835 | ncon1, ia1, ia2, at2con, | |||
836 | ga2la, FALSE0, dc, dcc, il_local, ireq); | |||
837 | } | |||
838 | } | |||
839 | } | |||
840 | } | |||
841 | } | |||
842 | ||||
843 | static void atoms_to_settles(gmx_domdec_t *dd, | |||
844 | const gmx_mtop_t *mtop, | |||
845 | const int *cginfo, | |||
846 | const int **at2settle_mt, | |||
847 | int cg_start, int cg_end, | |||
848 | t_ilist *ils_local, | |||
849 | ind_req_t *ireq) | |||
850 | { | |||
851 | gmx_ga2la_t ga2la; | |||
852 | gmx_mtop_atomlookup_t alook; | |||
853 | int settle; | |||
854 | int nral, sa; | |||
855 | int cg, a, a_gl, a_glsa, a_gls[3], a_locs[3]; | |||
856 | int mb, molnr, a_mol, offset; | |||
857 | const gmx_molblock_t *molb; | |||
858 | const t_iatom *ia1; | |||
859 | gmx_bool a_home[3]; | |||
860 | int nlocal; | |||
861 | gmx_bool bAssign; | |||
862 | ||||
863 | ga2la = dd->ga2la; | |||
864 | ||||
865 | alook = gmx_mtop_atomlookup_settle_init(mtop); | |||
866 | ||||
867 | nral = NRAL(F_SETTLE)(interaction_function[(F_SETTLE)].nratoms); | |||
868 | ||||
869 | for (cg = cg_start; cg < cg_end; cg++) | |||
870 | { | |||
871 | if (GET_CGINFO_SETTLE(cginfo[cg])( (cginfo[cg]) & (1<<21))) | |||
872 | { | |||
873 | for (a = dd->cgindex[cg]; a < dd->cgindex[cg+1]; a++) | |||
874 | { | |||
875 | a_gl = dd->gatindex[a]; | |||
876 | ||||
877 | gmx_mtop_atomnr_to_molblock_ind(alook, a_gl, &mb, &molnr, &a_mol); | |||
878 | molb = &mtop->molblock[mb]; | |||
879 | ||||
880 | settle = at2settle_mt[molb->type][a_mol]; | |||
881 | ||||
882 | if (settle >= 0) | |||
883 | { | |||
884 | offset = a_gl - a_mol; | |||
885 | ||||
886 | ia1 = mtop->moltype[molb->type].ilist[F_SETTLE].iatoms; | |||
887 | ||||
888 | bAssign = FALSE0; | |||
889 | nlocal = 0; | |||
890 | for (sa = 0; sa < nral; sa++) | |||
891 | { | |||
892 | a_glsa = offset + ia1[settle*(1+nral)+1+sa]; | |||
893 | a_gls[sa] = a_glsa; | |||
894 | a_home[sa] = ga2la_get_home(ga2la, a_glsa, &a_locs[sa]); | |||
895 | if (a_home[sa]) | |||
896 | { | |||
897 | if (nlocal == 0 && a_gl == a_glsa) | |||
898 | { | |||
899 | bAssign = TRUE1; | |||
900 | } | |||
901 | nlocal++; | |||
902 | } | |||
903 | } | |||
904 | ||||
905 | if (bAssign) | |||
906 | { | |||
907 | if (ils_local->nr+1+nral > ils_local->nalloc) | |||
908 | { | |||
909 | ils_local->nalloc = over_alloc_dd(ils_local->nr+1+nral); | |||
910 | srenew(ils_local->iatoms, ils_local->nalloc)(ils_local->iatoms) = save_realloc("ils_local->iatoms", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/domdec_con.c" , 910, (ils_local->iatoms), (ils_local->nalloc), sizeof (*(ils_local->iatoms))); | |||
911 | } | |||
912 | ||||
913 | ils_local->iatoms[ils_local->nr++] = ia1[settle*4]; | |||
914 | ||||
915 | for (sa = 0; sa < nral; sa++) | |||
916 | { | |||
917 | if (ga2la_get_home(ga2la, a_gls[sa], &a_locs[sa])) | |||
918 | { | |||
919 | ils_local->iatoms[ils_local->nr++] = a_locs[sa]; | |||
920 | } | |||
921 | else | |||
922 | { | |||
923 | ils_local->iatoms[ils_local->nr++] = -a_gls[sa] - 1; | |||
924 | /* Add this non-home atom to the list */ | |||
925 | if (ireq->n+1 > ireq->nalloc) | |||
926 | { | |||
927 | ireq->nalloc = over_alloc_large(ireq->n+1)(int)(1.19*(ireq->n+1) + 1000); | |||
928 | srenew(ireq->ind, ireq->nalloc)(ireq->ind) = save_realloc("ireq->ind", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/domdec_con.c" , 928, (ireq->ind), (ireq->nalloc), sizeof(*(ireq->ind ))); | |||
929 | } | |||
930 | ireq->ind[ireq->n++] = a_gls[sa]; | |||
931 | /* A check on double atom requests is | |||
932 | * not required for settle. | |||
933 | */ | |||
934 | } | |||
935 | } | |||
936 | } | |||
937 | } | |||
938 | } | |||
939 | } | |||
940 | } | |||
941 | ||||
942 | gmx_mtop_atomlookup_destroy(alook); | |||
943 | } | |||
944 | ||||
945 | static void atoms_to_constraints(gmx_domdec_t *dd, | |||
946 | const gmx_mtop_t *mtop, | |||
947 | const int *cginfo, | |||
948 | const t_blocka *at2con_mt, int nrec, | |||
949 | t_ilist *ilc_local, | |||
950 | ind_req_t *ireq) | |||
951 | { | |||
952 | const t_blocka *at2con; | |||
953 | gmx_ga2la_t ga2la; | |||
954 | gmx_mtop_atomlookup_t alook; | |||
955 | int ncon1; | |||
956 | gmx_molblock_t *molb; | |||
957 | t_iatom *ia1, *ia2, *iap; | |||
958 | int nhome, cg, a, a_gl, a_mol, a_loc, b_lo, offset, mb, molnr, b_mol, i, con, con_offset; | |||
959 | gmx_domdec_constraints_t *dc; | |||
960 | gmx_domdec_specat_comm_t *dcc; | |||
961 | ||||
962 | dc = dd->constraints; | |||
963 | dcc = dd->constraint_comm; | |||
964 | ||||
965 | ga2la = dd->ga2la; | |||
966 | ||||
967 | alook = gmx_mtop_atomlookup_init(mtop); | |||
968 | ||||
969 | nhome = 0; | |||
970 | for (cg = 0; cg < dd->ncg_home; cg++) | |||
971 | { | |||
972 | if (GET_CGINFO_CONSTR(cginfo[cg])( (cginfo[cg]) & (1<<20))) | |||
973 | { | |||
974 | for (a = dd->cgindex[cg]; a < dd->cgindex[cg+1]; a++) | |||
975 | { | |||
976 | a_gl = dd->gatindex[a]; | |||
977 | ||||
978 | gmx_mtop_atomnr_to_molblock_ind(alook, a_gl, &mb, &molnr, &a_mol); | |||
979 | molb = &mtop->molblock[mb]; | |||
980 | ||||
981 | ncon1 = mtop->moltype[molb->type].ilist[F_CONSTR].nr/NRAL(F_SETTLE)(interaction_function[(F_SETTLE)].nratoms); | |||
982 | ||||
983 | ia1 = mtop->moltype[molb->type].ilist[F_CONSTR].iatoms; | |||
984 | ia2 = mtop->moltype[molb->type].ilist[F_CONSTRNC].iatoms; | |||
985 | ||||
986 | /* Calculate the global constraint number offset for the molecule. | |||
987 | * This is only required for the global index to make sure | |||
988 | * that we use each constraint only once. | |||
989 | */ | |||
990 | con_offset = | |||
991 | dc->molb_con_offset[mb] + molnr*dc->molb_ncon_mol[mb]; | |||
992 | ||||
993 | /* The global atom number offset for this molecule */ | |||
994 | offset = a_gl - a_mol; | |||
995 | at2con = &at2con_mt[molb->type]; | |||
996 | for (i = at2con->index[a_mol]; i < at2con->index[a_mol+1]; i++) | |||
997 | { | |||
998 | con = at2con->a[i]; | |||
999 | iap = constr_iatomptr(ncon1, ia1, ia2, con)((con) < (ncon1) ? (ia1)+(con)*3 : (ia2)+(con-ncon1)*3); | |||
1000 | if (a_mol == iap[1]) | |||
1001 | { | |||
1002 | b_mol = iap[2]; | |||
1003 | } | |||
1004 | else | |||
1005 | { | |||
1006 | b_mol = iap[1]; | |||
1007 | } | |||
1008 | if (ga2la_get_home(ga2la, offset+b_mol, &a_loc)) | |||
1009 | { | |||
1010 | /* Add this fully home constraint at the first atom */ | |||
1011 | if (a_mol < b_mol) | |||
1012 | { | |||
1013 | if (dc->ncon+1 > dc->con_nalloc) | |||
1014 | { | |||
1015 | dc->con_nalloc = over_alloc_large(dc->ncon+1)(int)(1.19*(dc->ncon+1) + 1000); | |||
1016 | srenew(dc->con_gl, dc->con_nalloc)(dc->con_gl) = save_realloc("dc->con_gl", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/domdec_con.c" , 1016, (dc->con_gl), (dc->con_nalloc), sizeof(*(dc-> con_gl))); | |||
1017 | srenew(dc->con_nlocat, dc->con_nalloc)(dc->con_nlocat) = save_realloc("dc->con_nlocat", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/domdec_con.c" , 1017, (dc->con_nlocat), (dc->con_nalloc), sizeof(*(dc ->con_nlocat))); | |||
1018 | } | |||
1019 | dc->con_gl[dc->ncon] = con_offset + con; | |||
1020 | dc->con_nlocat[dc->ncon] = 2; | |||
1021 | if (ilc_local->nr + 3 > ilc_local->nalloc) | |||
1022 | { | |||
1023 | ilc_local->nalloc = over_alloc_dd(ilc_local->nr + 3); | |||
1024 | srenew(ilc_local->iatoms, ilc_local->nalloc)(ilc_local->iatoms) = save_realloc("ilc_local->iatoms", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/domdec_con.c" , 1024, (ilc_local->iatoms), (ilc_local->nalloc), sizeof (*(ilc_local->iatoms))); | |||
1025 | } | |||
1026 | b_lo = a_loc; | |||
1027 | ilc_local->iatoms[ilc_local->nr++] = iap[0]; | |||
1028 | ilc_local->iatoms[ilc_local->nr++] = (a_gl == iap[1] ? a : b_lo); | |||
1029 | ilc_local->iatoms[ilc_local->nr++] = (a_gl == iap[1] ? b_lo : a ); | |||
1030 | dc->ncon++; | |||
1031 | nhome++; | |||
1032 | } | |||
1033 | } | |||
1034 | else | |||
1035 | { | |||
1036 | /* We need the nrec constraints coupled to this constraint, | |||
1037 | * so we need to walk out of the home cell by nrec+1 atoms, | |||
1038 | * since already atom bg is not locally present. | |||
1039 | * Therefore we call walk_out with nrec recursions to go | |||
1040 | * after this first call. | |||
1041 | */ | |||
1042 | walk_out(con, con_offset, b_mol, offset, nrec, | |||
1043 | ncon1, ia1, ia2, at2con, | |||
1044 | dd->ga2la, TRUE1, dc, dcc, ilc_local, ireq); | |||
1045 | } | |||
1046 | } | |||
1047 | } | |||
1048 | } | |||
1049 | } | |||
1050 | ||||
1051 | gmx_mtop_atomlookup_destroy(alook); | |||
1052 | ||||
1053 | if (debug) | |||
1054 | { | |||
1055 | fprintf(debug, | |||
1056 | "Constraints: home %3d border %3d atoms: %3d\n", | |||
1057 | nhome, dc->ncon-nhome, | |||
1058 | dd->constraint_comm ? ireq->n : 0); | |||
1059 | } | |||
1060 | } | |||
1061 | ||||
1062 | int dd_make_local_constraints(gmx_domdec_t *dd, int at_start, | |||
1063 | const gmx_mtop_t *mtop, | |||
1064 | const int *cginfo, | |||
1065 | gmx_constr_t constr, int nrec, | |||
1066 | t_ilist *il_local) | |||
1067 | { | |||
1068 | gmx_domdec_constraints_t *dc; | |||
1069 | t_ilist *ilc_local, *ils_local; | |||
1070 | ind_req_t *ireq; | |||
1071 | const t_blocka *at2con_mt; | |||
1072 | const int **at2settle_mt; | |||
1073 | gmx_hash_t ga2la_specat; | |||
1074 | int at_end, i, j; | |||
1075 | t_iatom *iap; | |||
1076 | ||||
1077 | dc = dd->constraints; | |||
1078 | ||||
1079 | ilc_local = &il_local[F_CONSTR]; | |||
1080 | ils_local = &il_local[F_SETTLE]; | |||
1081 | ||||
1082 | dc->ncon = 0; | |||
1083 | ilc_local->nr = 0; | |||
1084 | if (dd->constraint_comm) | |||
| ||||
1085 | { | |||
1086 | at2con_mt = atom2constraints_moltype(constr); | |||
1087 | ireq = &dd->constraint_comm->ireq[0]; | |||
1088 | ireq->n = 0; | |||
1089 | } | |||
1090 | else | |||
1091 | { | |||
1092 | at2con_mt = NULL((void*)0); | |||
1093 | ireq = NULL((void*)0); | |||
1094 | } | |||
1095 | ||||
1096 | if (dd->bInterCGsettles) | |||
1097 | { | |||
1098 | at2settle_mt = atom2settle_moltype(constr); | |||
1099 | ils_local->nr = 0; | |||
1100 | } | |||
1101 | else | |||
1102 | { | |||
1103 | /* Settle works inside charge groups, we assigned them already */ | |||
1104 | at2settle_mt = NULL((void*)0); | |||
1105 | } | |||
1106 | ||||
1107 | if (at2settle_mt == NULL((void*)0)) | |||
1108 | { | |||
1109 | atoms_to_constraints(dd, mtop, cginfo, at2con_mt, nrec, | |||
1110 | ilc_local, ireq); | |||
1111 | } | |||
1112 | else | |||
1113 | { | |||
1114 | int t0_set; | |||
1115 | int thread; | |||
1116 | ||||
1117 | /* Do the constraints, if present, on the first thread. | |||
1118 | * Do the settles on all other threads. | |||
1119 | */ | |||
1120 | t0_set = ((at2con_mt != NULL((void*)0) && dc->nthread > 1) ? 1 : 0); | |||
1121 | ||||
1122 | #pragma omp parallel for num_threads(dc->nthread) schedule(static) | |||
1123 | for (thread = 0; thread < dc->nthread; thread++) | |||
1124 | { | |||
1125 | if (at2con_mt && thread == 0) | |||
1126 | { | |||
1127 | atoms_to_constraints(dd, mtop, cginfo, at2con_mt, nrec, | |||
1128 | ilc_local, ireq); | |||
1129 | } | |||
1130 | ||||
1131 | if (thread >= t0_set) | |||
1132 | { | |||
1133 | int cg0, cg1; | |||
1134 | t_ilist *ilst; | |||
1135 | ind_req_t *ireqt; | |||
1136 | ||||
1137 | /* Distribute the settle check+assignments over | |||
1138 | * dc->nthread or dc->nthread-1 threads. | |||
1139 | */ | |||
1140 | cg0 = (dd->ncg_home*(thread-t0_set ))/(dc->nthread-t0_set); | |||
1141 | cg1 = (dd->ncg_home*(thread-t0_set+1))/(dc->nthread-t0_set); | |||
1142 | ||||
1143 | if (thread == t0_set) | |||
1144 | { | |||
1145 | ilst = ils_local; | |||
1146 | } | |||
1147 | else | |||
1148 | { | |||
1149 | ilst = &dc->ils[thread]; | |||
1150 | } | |||
1151 | ilst->nr = 0; | |||
1152 | ||||
1153 | ireqt = &dd->constraint_comm->ireq[thread]; | |||
| ||||
1154 | if (thread > 0) | |||
1155 | { | |||
1156 | ireqt->n = 0; | |||
1157 | } | |||
1158 | ||||
1159 | atoms_to_settles(dd, mtop, cginfo, at2settle_mt, | |||
1160 | cg0, cg1, | |||
1161 | ilst, ireqt); | |||
1162 | } | |||
1163 | } | |||
1164 | ||||
1165 | /* Combine the generate settles and requested indices */ | |||
1166 | for (thread = 1; thread < dc->nthread; thread++) | |||
1167 | { | |||
1168 | t_ilist *ilst; | |||
1169 | ind_req_t *ireqt; | |||
1170 | int ia; | |||
1171 | ||||
1172 | if (thread > t0_set) | |||
1173 | { | |||
1174 | ilst = &dc->ils[thread]; | |||
1175 | if (ils_local->nr + ilst->nr > ils_local->nalloc) | |||
1176 | { | |||
1177 | ils_local->nalloc = over_alloc_large(ils_local->nr + ilst->nr)(int)(1.19*(ils_local->nr + ilst->nr) + 1000); | |||
1178 | srenew(ils_local->iatoms, ils_local->nalloc)(ils_local->iatoms) = save_realloc("ils_local->iatoms", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/domdec_con.c" , 1178, (ils_local->iatoms), (ils_local->nalloc), sizeof (*(ils_local->iatoms))); | |||
1179 | } | |||
1180 | for (ia = 0; ia < ilst->nr; ia++) | |||
1181 | { | |||
1182 | ils_local->iatoms[ils_local->nr+ia] = ilst->iatoms[ia]; | |||
1183 | } | |||
1184 | ils_local->nr += ilst->nr; | |||
1185 | } | |||
1186 | ||||
1187 | ireqt = &dd->constraint_comm->ireq[thread]; | |||
1188 | if (ireq->n+ireqt->n > ireq->nalloc) | |||
1189 | { | |||
1190 | ireq->nalloc = over_alloc_large(ireq->n+ireqt->n)(int)(1.19*(ireq->n+ireqt->n) + 1000); | |||
1191 | srenew(ireq->ind, ireq->nalloc)(ireq->ind) = save_realloc("ireq->ind", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/domdec_con.c" , 1191, (ireq->ind), (ireq->nalloc), sizeof(*(ireq-> ind))); | |||
1192 | } | |||
1193 | for (ia = 0; ia < ireqt->n; ia++) | |||
1194 | { | |||
1195 | ireq->ind[ireq->n+ia] = ireqt->ind[ia]; | |||
1196 | } | |||
1197 | ireq->n += ireqt->n; | |||
1198 | } | |||
1199 | ||||
1200 | if (debug) | |||
1201 | { | |||
1202 | fprintf(debug, "Settles: total %3d\n", ils_local->nr/4); | |||
1203 | } | |||
1204 | } | |||
1205 | ||||
1206 | if (dd->constraint_comm) | |||
1207 | { | |||
1208 | int nral1; | |||
1209 | ||||
1210 | at_end = | |||
1211 | setup_specat_communication(dd, ireq, dd->constraint_comm, | |||
1212 | dd->constraints->ga2la, | |||
1213 | at_start, 2, | |||
1214 | "constraint", " or lincs-order"); | |||
1215 | ||||
1216 | /* Fill in the missing indices */ | |||
1217 | ga2la_specat = dd->constraints->ga2la; | |||
1218 | ||||
1219 | nral1 = 1 + NRAL(F_CONSTR)(interaction_function[(F_CONSTR)].nratoms); | |||
1220 | for (i = 0; i < ilc_local->nr; i += nral1) | |||
1221 | { | |||
1222 | iap = ilc_local->iatoms + i; | |||
1223 | for (j = 1; j < nral1; j++) | |||
1224 | { | |||
1225 | if (iap[j] < 0) | |||
1226 | { | |||
1227 | iap[j] = gmx_hash_get_minone(ga2la_specat, -iap[j]-1); | |||
1228 | } | |||
1229 | } | |||
1230 | } | |||
1231 | ||||
1232 | nral1 = 1 + NRAL(F_SETTLE)(interaction_function[(F_SETTLE)].nratoms); | |||
1233 | for (i = 0; i < ils_local->nr; i += nral1) | |||
1234 | { | |||
1235 | iap = ils_local->iatoms + i; | |||
1236 | for (j = 1; j < nral1; j++) | |||
1237 | { | |||
1238 | if (iap[j] < 0) | |||
1239 | { | |||
1240 | iap[j] = gmx_hash_get_minone(ga2la_specat, -iap[j]-1); | |||
1241 | } | |||
1242 | } | |||
1243 | } | |||
1244 | } | |||
1245 | else | |||
1246 | { | |||
1247 | at_end = at_start; | |||
1248 | } | |||
1249 | ||||
1250 | return at_end; | |||
1251 | } | |||
1252 | ||||
1253 | int dd_make_local_vsites(gmx_domdec_t *dd, int at_start, t_ilist *lil) | |||
1254 | { | |||
1255 | gmx_domdec_specat_comm_t *spac; | |||
1256 | ind_req_t *ireq; | |||
1257 | gmx_hash_t ga2la_specat; | |||
1258 | int ftype, nral, i, j, gat, a; | |||
1259 | t_ilist *lilf; | |||
1260 | t_iatom *iatoms; | |||
1261 | int at_end; | |||
1262 | ||||
1263 | spac = dd->vsite_comm; | |||
1264 | ireq = &spac->ireq[0]; | |||
1265 | ga2la_specat = dd->ga2la_vsite; | |||
1266 | ||||
1267 | ireq->n = 0; | |||
1268 | /* Loop over all the home vsites */ | |||
1269 | for (ftype = 0; ftype < F_NRE; ftype++) | |||
1270 | { | |||
1271 | if (interaction_function[ftype].flags & IF_VSITE1<<1) | |||
1272 | { | |||
1273 | nral = NRAL(ftype)(interaction_function[(ftype)].nratoms); | |||
1274 | lilf = &lil[ftype]; | |||
1275 | for (i = 0; i < lilf->nr; i += 1+nral) | |||
1276 | { | |||
1277 | iatoms = lilf->iatoms + i; | |||
1278 | /* Check if we have the other atoms */ | |||
1279 | for (j = 1; j < 1+nral; j++) | |||
1280 | { | |||
1281 | if (iatoms[j] < 0) | |||
1282 | { | |||
1283 | /* This is not a home atom, | |||
1284 | * we need to ask our neighbors. | |||
1285 | */ | |||
1286 | a = -iatoms[j] - 1; | |||
1287 | /* Check to not ask for the same atom more than once */ | |||
1288 | if (gmx_hash_get_minone(dd->ga2la_vsite, a) == -1) | |||
1289 | { | |||
1290 | /* Add this non-home atom to the list */ | |||
1291 | if (ireq->n+1 > ireq->nalloc) | |||
1292 | { | |||
1293 | ireq->nalloc = over_alloc_large(ireq->n+1)(int)(1.19*(ireq->n+1) + 1000); | |||
1294 | srenew(ireq->ind, ireq->nalloc)(ireq->ind) = save_realloc("ireq->ind", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/domdec_con.c" , 1294, (ireq->ind), (ireq->nalloc), sizeof(*(ireq-> ind))); | |||
1295 | } | |||
1296 | ireq->ind[ireq->n++] = a; | |||
1297 | /* Temporarily mark with -2, | |||
1298 | * we get the index later. | |||
1299 | */ | |||
1300 | gmx_hash_set(ga2la_specat, a, -2); | |||
1301 | } | |||
1302 | } | |||
1303 | } | |||
1304 | } | |||
1305 | } | |||
1306 | } | |||
1307 | ||||
1308 | at_end = setup_specat_communication(dd, ireq, dd->vsite_comm, ga2la_specat, | |||
1309 | at_start, 1, "vsite", ""); | |||
1310 | ||||
1311 | /* Fill in the missing indices */ | |||
1312 | for (ftype = 0; ftype < F_NRE; ftype++) | |||
1313 | { | |||
1314 | if (interaction_function[ftype].flags & IF_VSITE1<<1) | |||
1315 | { | |||
1316 | nral = NRAL(ftype)(interaction_function[(ftype)].nratoms); | |||
1317 | lilf = &lil[ftype]; | |||
1318 | for (i = 0; i < lilf->nr; i += 1+nral) | |||
1319 | { | |||
1320 | iatoms = lilf->iatoms + i; | |||
1321 | for (j = 1; j < 1+nral; j++) | |||
1322 | { | |||
1323 | if (iatoms[j] < 0) | |||
1324 | { | |||
1325 | iatoms[j] = gmx_hash_get_minone(ga2la_specat, -iatoms[j]-1); | |||
1326 | } | |||
1327 | } | |||
1328 | } | |||
1329 | } | |||
1330 | } | |||
1331 | ||||
1332 | return at_end; | |||
1333 | } | |||
1334 | ||||
1335 | static gmx_domdec_specat_comm_t *specat_comm_init(int nthread) | |||
1336 | { | |||
1337 | gmx_domdec_specat_comm_t *spac; | |||
1338 | ||||
1339 | snew(spac, 1)(spac) = save_calloc("spac", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/domdec_con.c" , 1339, (1), sizeof(*(spac))); | |||
1340 | spac->nthread = nthread; | |||
1341 | snew(spac->ireq, spac->nthread)(spac->ireq) = save_calloc("spac->ireq", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/domdec_con.c" , 1341, (spac->nthread), sizeof(*(spac->ireq))); | |||
1342 | ||||
1343 | return spac; | |||
1344 | } | |||
1345 | ||||
1346 | void init_domdec_constraints(gmx_domdec_t *dd, | |||
1347 | gmx_mtop_t *mtop) | |||
1348 | { | |||
1349 | gmx_domdec_constraints_t *dc; | |||
1350 | gmx_molblock_t *molb; | |||
1351 | int mb, ncon, c, a; | |||
1352 | ||||
1353 | if (debug) | |||
1354 | { | |||
1355 | fprintf(debug, "Begin init_domdec_constraints\n"); | |||
1356 | } | |||
1357 | ||||
1358 | snew(dd->constraints, 1)(dd->constraints) = save_calloc("dd->constraints", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/domdec_con.c" , 1358, (1), sizeof(*(dd->constraints))); | |||
1359 | dc = dd->constraints; | |||
1360 | ||||
1361 | snew(dc->molb_con_offset, mtop->nmolblock)(dc->molb_con_offset) = save_calloc("dc->molb_con_offset" , "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/domdec_con.c" , 1361, (mtop->nmolblock), sizeof(*(dc->molb_con_offset ))); | |||
1362 | snew(dc->molb_ncon_mol, mtop->nmolblock)(dc->molb_ncon_mol) = save_calloc("dc->molb_ncon_mol", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/domdec_con.c" , 1362, (mtop->nmolblock), sizeof(*(dc->molb_ncon_mol)) ); | |||
1363 | ||||
1364 | ncon = 0; | |||
1365 | for (mb = 0; mb < mtop->nmolblock; mb++) | |||
1366 | { | |||
1367 | molb = &mtop->molblock[mb]; | |||
1368 | dc->molb_con_offset[mb] = ncon; | |||
1369 | dc->molb_ncon_mol[mb] = | |||
1370 | mtop->moltype[molb->type].ilist[F_CONSTR].nr/3 + | |||
1371 | mtop->moltype[molb->type].ilist[F_CONSTRNC].nr/3; | |||
1372 | ncon += molb->nmol*dc->molb_ncon_mol[mb]; | |||
1373 | } | |||
1374 | ||||
1375 | if (ncon > 0) | |||
1376 | { | |||
1377 | snew(dc->gc_req, ncon)(dc->gc_req) = save_calloc("dc->gc_req", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/domdec_con.c" , 1377, (ncon), sizeof(*(dc->gc_req))); | |||
1378 | for (c = 0; c < ncon; c++) | |||
1379 | { | |||
1380 | dc->gc_req[c] = 0; | |||
1381 | } | |||
1382 | } | |||
1383 | ||||
1384 | /* Use a hash table for the global to local index. | |||
1385 | * The number of keys is a rough estimate, it will be optimized later. | |||
1386 | */ | |||
1387 | dc->ga2la = gmx_hash_init(min(mtop->natoms/20,(((mtop->natoms/20) < (mtop->natoms/(2*dd->nnodes ))) ? (mtop->natoms/20) : (mtop->natoms/(2*dd->nnodes )) ) | |||
1388 | mtop->natoms/(2*dd->nnodes))(((mtop->natoms/20) < (mtop->natoms/(2*dd->nnodes ))) ? (mtop->natoms/20) : (mtop->natoms/(2*dd->nnodes )) )); | |||
1389 | ||||
1390 | dc->nthread = gmx_omp_nthreads_get(emntDomdec); | |||
1391 | snew(dc->ils, dc->nthread)(dc->ils) = save_calloc("dc->ils", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/domdec_con.c" , 1391, (dc->nthread), sizeof(*(dc->ils))); | |||
1392 | ||||
1393 | dd->constraint_comm = specat_comm_init(dc->nthread); | |||
1394 | } | |||
1395 | ||||
1396 | void init_domdec_vsites(gmx_domdec_t *dd, int n_intercg_vsite) | |||
1397 | { | |||
1398 | int i; | |||
1399 | gmx_domdec_constraints_t *dc; | |||
1400 | ||||
1401 | if (debug) | |||
1402 | { | |||
1403 | fprintf(debug, "Begin init_domdec_vsites\n"); | |||
1404 | } | |||
1405 | ||||
1406 | /* Use a hash table for the global to local index. | |||
1407 | * The number of keys is a rough estimate, it will be optimized later. | |||
1408 | */ | |||
1409 | dd->ga2la_vsite = gmx_hash_init(min(n_intercg_vsite/20,(((n_intercg_vsite/20) < (n_intercg_vsite/(2*dd->nnodes ))) ? (n_intercg_vsite/20) : (n_intercg_vsite/(2*dd->nnodes )) ) | |||
1410 | n_intercg_vsite/(2*dd->nnodes))(((n_intercg_vsite/20) < (n_intercg_vsite/(2*dd->nnodes ))) ? (n_intercg_vsite/20) : (n_intercg_vsite/(2*dd->nnodes )) )); | |||
1411 | ||||
1412 | dd->vsite_comm = specat_comm_init(1); | |||
1413 | } |