| File: | gromacs/mdlib/domdec_con.c |
| Location: | line 560, column 9 |
| Description: | Value stored to 'bFirst' is never read |
| 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); |
Value stored to 'bFirst' is never read | |
| 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 | } |