2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
41 #include "gromacs/legacyheaders/constr.h"
42 #include "gromacs/legacyheaders/domdec.h"
43 #include "gromacs/legacyheaders/domdec_network.h"
44 #include "gromacs/legacyheaders/gmx_ga2la.h"
45 #include "gromacs/legacyheaders/gmx_hash.h"
46 #include "gromacs/legacyheaders/gmx_omp_nthreads.h"
47 #include "gromacs/legacyheaders/macros.h"
48 #include "gromacs/legacyheaders/types/commrec.h"
49 #include "gromacs/math/vec.h"
50 #include "gromacs/pbcutil/ishift.h"
51 #include "gromacs/topology/mtop_util.h"
52 #include "gromacs/utility/fatalerror.h"
53 #include "gromacs/utility/gmxassert.h"
54 #include "gromacs/utility/smalloc.h"
69 typedef struct gmx_domdec_specat_comm {
70 /* The number of indices to receive during the setup */
72 /* The atoms to send */
73 gmx_specatsend_t spas[DIM][2];
83 /* The range in the local buffer(s) for received atoms */
87 /* The atom indices we need from the surrounding cells.
88 * We can gather the indices over nthread threads.
92 } gmx_domdec_specat_comm_t;
94 typedef struct gmx_domdec_constraints {
97 /* The fully local and connected constraints */
99 /* The global constraint number, only required for clearing gc_req */
103 /* Boolean that tells if a global constraint index has been requested */
105 /* Global to local communicated constraint atom only index */
108 /* Multi-threading stuff */
111 } gmx_domdec_constraints_t;
114 static void dd_move_f_specat(gmx_domdec_t *dd, gmx_domdec_specat_comm_t *spac,
115 rvec *f, rvec *fshift)
117 gmx_specatsend_t *spas;
119 int n, n0, n1, d, dim, dir, i;
122 gmx_bool bPBC, bScrew;
125 for (d = dd->ndim-1; d >= 0; d--)
130 /* Pulse the grid forward and backward */
131 spas = spac->spas[d];
136 /* Send and receive the coordinates */
137 dd_sendrecv2_rvec(dd, d,
138 f+n+n1, n0, vbuf, spas[0].nsend,
139 f+n, n1, vbuf+spas[0].nsend, spas[1].nsend);
140 for (dir = 0; dir < 2; dir++)
142 bPBC = ((dir == 0 && dd->ci[dim] == 0) ||
143 (dir == 1 && dd->ci[dim] == dd->nc[dim]-1));
144 bScrew = (bPBC && dd->bScrewPBC && dim == XX);
146 spas = &spac->spas[d][dir];
147 /* Sum the buffer into the required forces */
148 if (!bPBC || (!bScrew && fshift == NULL))
150 for (i = 0; i < spas->nsend; i++)
152 rvec_inc(f[spas->a[i]], *vbuf);
159 vis[dim] = (dir == 0 ? 1 : -1);
163 /* Sum and add to shift forces */
164 for (i = 0; i < spas->nsend; i++)
166 rvec_inc(f[spas->a[i]], *vbuf);
167 rvec_inc(fshift[is], *vbuf);
173 /* Rotate the forces */
174 for (i = 0; i < spas->nsend; i++)
176 f[spas->a[i]][XX] += (*vbuf)[XX];
177 f[spas->a[i]][YY] -= (*vbuf)[YY];
178 f[spas->a[i]][ZZ] -= (*vbuf)[ZZ];
181 rvec_inc(fshift[is], *vbuf);
191 /* Two cells, so we only need to communicate one way */
192 spas = &spac->spas[d][0];
194 /* Send and receive the coordinates */
195 dd_sendrecv_rvec(dd, d, dddirForward,
196 f+n, spas->nrecv, spac->vbuf, spas->nsend);
197 /* Sum the buffer into the required forces */
198 if (dd->bScrewPBC && dim == XX &&
200 dd->ci[dim] == dd->nc[dim]-1))
202 for (i = 0; i < spas->nsend; i++)
204 /* Rotate the force */
205 f[spas->a[i]][XX] += spac->vbuf[i][XX];
206 f[spas->a[i]][YY] -= spac->vbuf[i][YY];
207 f[spas->a[i]][ZZ] -= spac->vbuf[i][ZZ];
212 for (i = 0; i < spas->nsend; i++)
214 rvec_inc(f[spas->a[i]], spac->vbuf[i]);
221 void dd_move_f_vsites(gmx_domdec_t *dd, rvec *f, rvec *fshift)
225 dd_move_f_specat(dd, dd->vsite_comm, f, fshift);
229 void dd_clear_f_vsites(gmx_domdec_t *dd, rvec *f)
235 for (i = dd->vsite_comm->at_start; i < dd->vsite_comm->at_end; i++)
242 static void dd_move_x_specat(gmx_domdec_t *dd, gmx_domdec_specat_comm_t *spac,
245 rvec *x1, gmx_bool bX1IsCoord)
247 gmx_specatsend_t *spas;
248 rvec *x, *vbuf, *rbuf;
249 int nvec, v, n, nn, ns0, ns1, nr0, nr1, nr, d, dim, dir, i;
250 gmx_bool bPBC, bScrew = FALSE;
251 rvec shift = {0, 0, 0};
260 for (d = 0; d < dd->ndim; d++)
265 /* Pulse the grid forward and backward */
267 for (dir = 0; dir < 2; dir++)
269 if (dir == 0 && dd->ci[dim] == 0)
272 bScrew = (dd->bScrewPBC && dim == XX);
273 copy_rvec(box[dim], shift);
275 else if (dir == 1 && dd->ci[dim] == dd->nc[dim]-1)
278 bScrew = (dd->bScrewPBC && dim == XX);
279 for (i = 0; i < DIM; i++)
281 shift[i] = -box[dim][i];
289 spas = &spac->spas[d][dir];
290 for (v = 0; v < nvec; v++)
292 x = (v == 0 ? x0 : x1);
293 /* Copy the required coordinates to the send buffer */
294 if (!bPBC || (v == 1 && !bX1IsCoord))
297 for (i = 0; i < spas->nsend; i++)
299 copy_rvec(x[spas->a[i]], *vbuf);
305 /* Shift coordinates */
306 for (i = 0; i < spas->nsend; i++)
308 rvec_add(x[spas->a[i]], shift, *vbuf);
314 /* Shift and rotate coordinates */
315 for (i = 0; i < spas->nsend; i++)
317 (*vbuf)[XX] = x[spas->a[i]][XX] + shift[XX];
318 (*vbuf)[YY] = box[YY][YY] - x[spas->a[i]][YY] + shift[YY];
319 (*vbuf)[ZZ] = box[ZZ][ZZ] - x[spas->a[i]][ZZ] + shift[ZZ];
325 /* Send and receive the coordinates */
326 spas = spac->spas[d];
333 dd_sendrecv2_rvec(dd, d,
334 spac->vbuf+ns0, ns1, x0+n, nr1,
335 spac->vbuf, ns0, x0+n+nr1, nr0);
339 /* Communicate both vectors in one buffer */
341 dd_sendrecv2_rvec(dd, d,
342 spac->vbuf+2*ns0, 2*ns1, rbuf, 2*nr1,
343 spac->vbuf, 2*ns0, rbuf+2*nr1, 2*nr0);
344 /* Split the buffer into the two vectors */
346 for (dir = 1; dir >= 0; dir--)
348 nr = spas[dir].nrecv;
349 for (v = 0; v < 2; v++)
351 x = (v == 0 ? x0 : x1);
352 for (i = 0; i < nr; i++)
354 copy_rvec(*rbuf, x[nn+i]);
365 spas = &spac->spas[d][0];
366 /* Copy the required coordinates to the send buffer */
368 for (v = 0; v < nvec; v++)
370 x = (v == 0 ? x0 : x1);
371 if (dd->bScrewPBC && dim == XX &&
372 (dd->ci[XX] == 0 || dd->ci[XX] == dd->nc[XX]-1))
374 /* Here we only perform the rotation, the rest of the pbc
375 * is handled in the constraint or viste routines.
377 for (i = 0; i < spas->nsend; i++)
379 (*vbuf)[XX] = x[spas->a[i]][XX];
380 (*vbuf)[YY] = box[YY][YY] - x[spas->a[i]][YY];
381 (*vbuf)[ZZ] = box[ZZ][ZZ] - x[spas->a[i]][ZZ];
387 for (i = 0; i < spas->nsend; i++)
389 copy_rvec(x[spas->a[i]], *vbuf);
394 /* Send and receive the coordinates */
397 dd_sendrecv_rvec(dd, d, dddirBackward,
398 spac->vbuf, spas->nsend, x0+n, spas->nrecv);
402 /* Communicate both vectors in one buffer */
404 dd_sendrecv_rvec(dd, d, dddirBackward,
405 spac->vbuf, 2*spas->nsend, rbuf, 2*spas->nrecv);
406 /* Split the buffer into the two vectors */
408 for (v = 0; v < 2; v++)
410 x = (v == 0 ? x0 : x1);
411 for (i = 0; i < nr; i++)
413 copy_rvec(*rbuf, x[n+i]);
423 void dd_move_x_constraints(gmx_domdec_t *dd, matrix box,
424 rvec *x0, rvec *x1, gmx_bool bX1IsCoord)
426 if (dd->constraint_comm)
428 dd_move_x_specat(dd, dd->constraint_comm, box, x0, x1, bX1IsCoord);
432 void dd_move_x_vsites(gmx_domdec_t *dd, matrix box, rvec *x)
436 dd_move_x_specat(dd, dd->vsite_comm, box, x, NULL, FALSE);
440 int *dd_constraints_nlocalatoms(gmx_domdec_t *dd)
444 return dd->constraints->con_nlocat;
452 void dd_clear_local_constraint_indices(gmx_domdec_t *dd)
454 gmx_domdec_constraints_t *dc;
457 dc = dd->constraints;
459 for (i = 0; i < dc->ncon; i++)
461 dc->gc_req[dc->con_gl[i]] = 0;
464 if (dd->constraint_comm)
466 gmx_hash_clear_and_optimize(dc->ga2la);
470 void dd_clear_local_vsite_indices(gmx_domdec_t *dd)
474 gmx_hash_clear_and_optimize(dd->ga2la_vsite);
478 static int setup_specat_communication(gmx_domdec_t *dd,
480 gmx_domdec_specat_comm_t *spac,
481 gmx_hash_t ga2la_specat,
484 const char *specat_type,
487 int nsend[2], nlast, nsend_zero[2] = {0, 0}, *nsend_ptr;
488 int d, dim, ndir, dir, nr, ns, i, nrecv_local, n0, start, indr, ind, buf[2];
489 int nat_tot_specat, nat_tot_prev, nalloc_old;
491 gmx_specatsend_t *spas;
495 fprintf(debug, "Begin setup_specat_communication for %s\n", specat_type);
498 /* nsend[0]: the number of atoms requested by this node only,
499 * we communicate this for more efficients checks
500 * nsend[1]: the total number of requested atoms
505 for (d = dd->ndim-1; d >= 0; d--)
507 /* Pulse the grid forward and backward */
509 bPBC = (dim < dd->npbcdim);
510 if (dd->nc[dim] == 2)
512 /* Only 2 cells, so we only need to communicate once */
519 for (dir = 0; dir < ndir; dir++)
523 ((dir == 0 && dd->ci[dim] == dd->nc[dim] - 1) ||
524 (dir == 1 && dd->ci[dim] == 0)))
526 /* No pbc: the fist/last cell should not request atoms */
527 nsend_ptr = nsend_zero;
533 /* Communicate the number of indices */
534 dd_sendrecv_int(dd, d, dir == 0 ? dddirForward : dddirBackward,
535 nsend_ptr, 2, spac->nreq[d][dir], 2);
536 nr = spac->nreq[d][dir][1];
537 if (nlast+nr > ireq->nalloc)
539 ireq->nalloc = over_alloc_dd(nlast+nr);
540 srenew(ireq->ind, ireq->nalloc);
542 /* Communicate the indices */
543 dd_sendrecv_int(dd, d, dir == 0 ? dddirForward : dddirBackward,
544 ireq->ind, nsend_ptr[1], ireq->ind+nlast, nr);
551 fprintf(debug, "Communicated the counts\n");
554 /* Search for the requested atoms and communicate the indices we have */
555 nat_tot_specat = at_start;
557 for (d = 0; d < dd->ndim; d++)
559 /* Pulse the grid forward and backward */
560 if (dd->dim[d] >= dd->npbcdim || dd->nc[dd->dim[d]] > 2)
568 nat_tot_prev = nat_tot_specat;
569 for (dir = ndir-1; dir >= 0; dir--)
571 if (nat_tot_specat > spac->bSendAtom_nalloc)
573 nalloc_old = spac->bSendAtom_nalloc;
574 spac->bSendAtom_nalloc = over_alloc_dd(nat_tot_specat);
575 srenew(spac->bSendAtom, spac->bSendAtom_nalloc);
576 for (i = nalloc_old; i < spac->bSendAtom_nalloc; i++)
578 spac->bSendAtom[i] = FALSE;
581 spas = &spac->spas[d][dir];
582 n0 = spac->nreq[d][dir][0];
583 nr = spac->nreq[d][dir][1];
586 fprintf(debug, "dim=%d, dir=%d, searching for %d atoms\n",
592 for (i = 0; i < nr; i++)
594 indr = ireq->ind[start+i];
596 /* Check if this is a home atom and if so ind will be set */
597 if (!ga2la_get_home(dd->ga2la, indr, &ind))
599 /* Search in the communicated atoms */
600 ind = gmx_hash_get_minone(ga2la_specat, indr);
604 if (i < n0 || !spac->bSendAtom[ind])
606 if (spas->nsend+1 > spas->a_nalloc)
608 spas->a_nalloc = over_alloc_large(spas->nsend+1);
609 srenew(spas->a, spas->a_nalloc);
611 /* Store the local index so we know which coordinates
614 spas->a[spas->nsend] = ind;
615 spac->bSendAtom[ind] = TRUE;
616 if (spas->nsend+1 > spac->ibuf_nalloc)
618 spac->ibuf_nalloc = over_alloc_large(spas->nsend+1);
619 srenew(spac->ibuf, spac->ibuf_nalloc);
621 /* Store the global index so we can send it now */
622 spac->ibuf[spas->nsend] = indr;
632 /* Clear the local flags */
633 for (i = 0; i < spas->nsend; i++)
635 spac->bSendAtom[spas->a[i]] = FALSE;
637 /* Send and receive the number of indices to communicate */
638 nsend[1] = spas->nsend;
639 dd_sendrecv_int(dd, d, dir == 0 ? dddirBackward : dddirForward,
643 fprintf(debug, "Send to rank %d, %d (%d) indices, "
644 "receive from rank %d, %d (%d) indices\n",
645 dd->neighbor[d][1-dir], nsend[1], nsend[0],
646 dd->neighbor[d][dir], buf[1], buf[0]);
649 for (i = 0; i < spas->nsend; i++)
651 fprintf(debug, " %d", spac->ibuf[i]+1);
653 fprintf(debug, "\n");
656 nrecv_local += buf[0];
657 spas->nrecv = buf[1];
658 if (nat_tot_specat + spas->nrecv > dd->gatindex_nalloc)
660 dd->gatindex_nalloc =
661 over_alloc_dd(nat_tot_specat + spas->nrecv);
662 srenew(dd->gatindex, dd->gatindex_nalloc);
664 /* Send and receive the indices */
665 dd_sendrecv_int(dd, d, dir == 0 ? dddirBackward : dddirForward,
666 spac->ibuf, spas->nsend,
667 dd->gatindex+nat_tot_specat, spas->nrecv);
668 nat_tot_specat += spas->nrecv;
671 /* Allocate the x/f communication buffers */
672 ns = spac->spas[d][0].nsend;
673 nr = spac->spas[d][0].nrecv;
676 ns += spac->spas[d][1].nsend;
677 nr += spac->spas[d][1].nrecv;
679 if (vbuf_fac*ns > spac->vbuf_nalloc)
681 spac->vbuf_nalloc = over_alloc_dd(vbuf_fac*ns);
682 srenew(spac->vbuf, spac->vbuf_nalloc);
684 if (vbuf_fac == 2 && vbuf_fac*nr > spac->vbuf2_nalloc)
686 spac->vbuf2_nalloc = over_alloc_dd(vbuf_fac*nr);
687 srenew(spac->vbuf2, spac->vbuf2_nalloc);
690 /* Make a global to local index for the communication atoms */
691 for (i = nat_tot_prev; i < nat_tot_specat; i++)
693 gmx_hash_change_or_set(ga2la_specat, dd->gatindex[i], i);
697 /* Check that in the end we got the number of atoms we asked for */
698 if (nrecv_local != ireq->n)
702 fprintf(debug, "Requested %d, received %d (tot recv %d)\n",
703 ireq->n, nrecv_local, nat_tot_specat-at_start);
706 for (i = 0; i < ireq->n; i++)
708 ind = gmx_hash_get_minone(ga2la_specat, ireq->ind[i]);
709 fprintf(debug, " %s%d",
710 (ind >= 0) ? "" : "!",
713 fprintf(debug, "\n");
716 fprintf(stderr, "\nDD cell %d %d %d: Neighboring cells do not have atoms:",
717 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
718 for (i = 0; i < ireq->n; i++)
720 if (gmx_hash_get_minone(ga2la_specat, ireq->ind[i]) < 0)
722 fprintf(stderr, " %d", ireq->ind[i]+1);
725 fprintf(stderr, "\n");
726 gmx_fatal(FARGS, "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.",
727 dd->ci[XX], dd->ci[YY], dd->ci[ZZ],
728 nrecv_local, ireq->n, specat_type,
729 specat_type, add_err,
730 dd->bGridJump ? " or use the -rcon option of mdrun" : "");
733 spac->at_start = at_start;
734 spac->at_end = nat_tot_specat;
738 fprintf(debug, "Done setup_specat_communication\n");
741 return nat_tot_specat;
744 static void walk_out(int con, int con_offset, int a, int offset, int nrec,
745 int ncon1, const t_iatom *ia1, const t_iatom *ia2,
746 const t_blocka *at2con,
747 const gmx_ga2la_t ga2la, gmx_bool bHomeConnect,
748 gmx_domdec_constraints_t *dc,
749 gmx_domdec_specat_comm_t *dcc,
753 int a1_gl, a2_gl, a_loc, i, coni, b;
756 if (dc->gc_req[con_offset+con] == 0)
758 /* Add this non-home constraint to the list */
759 if (dc->ncon+1 > dc->con_nalloc)
761 dc->con_nalloc = over_alloc_large(dc->ncon+1);
762 srenew(dc->con_gl, dc->con_nalloc);
763 srenew(dc->con_nlocat, dc->con_nalloc);
765 dc->con_gl[dc->ncon] = con_offset + con;
766 dc->con_nlocat[dc->ncon] = (bHomeConnect ? 1 : 0);
767 dc->gc_req[con_offset+con] = 1;
768 if (il_local->nr + 3 > il_local->nalloc)
770 il_local->nalloc = over_alloc_dd(il_local->nr+3);
771 srenew(il_local->iatoms, il_local->nalloc);
773 iap = constr_iatomptr(ncon1, ia1, ia2, con);
774 il_local->iatoms[il_local->nr++] = iap[0];
775 a1_gl = offset + iap[1];
776 a2_gl = offset + iap[2];
777 /* The following indexing code can probably be optizimed */
778 if (ga2la_get_home(ga2la, a1_gl, &a_loc))
780 il_local->iatoms[il_local->nr++] = a_loc;
784 /* We set this index later */
785 il_local->iatoms[il_local->nr++] = -a1_gl - 1;
787 if (ga2la_get_home(ga2la, a2_gl, &a_loc))
789 il_local->iatoms[il_local->nr++] = a_loc;
793 /* We set this index later */
794 il_local->iatoms[il_local->nr++] = -a2_gl - 1;
798 /* Check to not ask for the same atom more than once */
799 if (gmx_hash_get_minone(dc->ga2la, offset+a) == -1)
802 /* Add this non-home atom to the list */
803 if (ireq->n+1 > ireq->nalloc)
805 ireq->nalloc = over_alloc_large(ireq->n+1);
806 srenew(ireq->ind, ireq->nalloc);
808 ireq->ind[ireq->n++] = offset + a;
809 /* Temporarily mark with -2, we get the index later */
810 gmx_hash_set(dc->ga2la, offset+a, -2);
815 for (i = at2con->index[a]; i < at2con->index[a+1]; i++)
821 iap = constr_iatomptr(ncon1, ia1, ia2, coni);
830 if (!ga2la_get_home(ga2la, offset+b, &a_loc))
832 walk_out(coni, con_offset, b, offset, nrec-1,
833 ncon1, ia1, ia2, at2con,
834 ga2la, FALSE, dc, dcc, il_local, ireq);
841 static void atoms_to_settles(gmx_domdec_t *dd,
842 const gmx_mtop_t *mtop,
844 const int **at2settle_mt,
845 int cg_start, int cg_end,
850 gmx_mtop_atomlookup_t alook;
853 int cg, a, a_gl, a_glsa, a_gls[3], a_locs[3];
854 int mb, molnr, a_mol, offset;
855 const gmx_molblock_t *molb;
863 alook = gmx_mtop_atomlookup_settle_init(mtop);
865 nral = NRAL(F_SETTLE);
867 for (cg = cg_start; cg < cg_end; cg++)
869 if (GET_CGINFO_SETTLE(cginfo[cg]))
871 for (a = dd->cgindex[cg]; a < dd->cgindex[cg+1]; a++)
873 a_gl = dd->gatindex[a];
875 gmx_mtop_atomnr_to_molblock_ind(alook, a_gl, &mb, &molnr, &a_mol);
876 molb = &mtop->molblock[mb];
878 settle = at2settle_mt[molb->type][a_mol];
882 offset = a_gl - a_mol;
884 ia1 = mtop->moltype[molb->type].ilist[F_SETTLE].iatoms;
888 for (sa = 0; sa < nral; sa++)
890 a_glsa = offset + ia1[settle*(1+nral)+1+sa];
892 a_home[sa] = ga2la_get_home(ga2la, a_glsa, &a_locs[sa]);
895 if (nlocal == 0 && a_gl == a_glsa)
905 if (ils_local->nr+1+nral > ils_local->nalloc)
907 ils_local->nalloc = over_alloc_dd(ils_local->nr+1+nral);
908 srenew(ils_local->iatoms, ils_local->nalloc);
911 ils_local->iatoms[ils_local->nr++] = ia1[settle*4];
913 for (sa = 0; sa < nral; sa++)
915 if (ga2la_get_home(ga2la, a_gls[sa], &a_locs[sa]))
917 ils_local->iatoms[ils_local->nr++] = a_locs[sa];
921 ils_local->iatoms[ils_local->nr++] = -a_gls[sa] - 1;
922 /* Add this non-home atom to the list */
923 if (ireq->n+1 > ireq->nalloc)
925 ireq->nalloc = over_alloc_large(ireq->n+1);
926 srenew(ireq->ind, ireq->nalloc);
928 ireq->ind[ireq->n++] = a_gls[sa];
929 /* A check on double atom requests is
930 * not required for settle.
940 gmx_mtop_atomlookup_destroy(alook);
943 static void atoms_to_constraints(gmx_domdec_t *dd,
944 const gmx_mtop_t *mtop,
946 const t_blocka *at2con_mt, int nrec,
950 const t_blocka *at2con;
952 gmx_mtop_atomlookup_t alook;
954 gmx_molblock_t *molb;
955 t_iatom *ia1, *ia2, *iap;
956 int nhome, cg, a, a_gl, a_mol, a_loc, b_lo, offset, mb, molnr, b_mol, i, con, con_offset;
957 gmx_domdec_constraints_t *dc;
958 gmx_domdec_specat_comm_t *dcc;
960 dc = dd->constraints;
961 dcc = dd->constraint_comm;
965 alook = gmx_mtop_atomlookup_init(mtop);
968 for (cg = 0; cg < dd->ncg_home; cg++)
970 if (GET_CGINFO_CONSTR(cginfo[cg]))
972 for (a = dd->cgindex[cg]; a < dd->cgindex[cg+1]; a++)
974 a_gl = dd->gatindex[a];
976 gmx_mtop_atomnr_to_molblock_ind(alook, a_gl, &mb, &molnr, &a_mol);
977 molb = &mtop->molblock[mb];
979 ncon1 = mtop->moltype[molb->type].ilist[F_CONSTR].nr/NRAL(F_SETTLE);
981 ia1 = mtop->moltype[molb->type].ilist[F_CONSTR].iatoms;
982 ia2 = mtop->moltype[molb->type].ilist[F_CONSTRNC].iatoms;
984 /* Calculate the global constraint number offset for the molecule.
985 * This is only required for the global index to make sure
986 * that we use each constraint only once.
989 dc->molb_con_offset[mb] + molnr*dc->molb_ncon_mol[mb];
991 /* The global atom number offset for this molecule */
992 offset = a_gl - a_mol;
993 at2con = &at2con_mt[molb->type];
994 for (i = at2con->index[a_mol]; i < at2con->index[a_mol+1]; i++)
997 iap = constr_iatomptr(ncon1, ia1, ia2, con);
1006 if (ga2la_get_home(ga2la, offset+b_mol, &a_loc))
1008 /* Add this fully home constraint at the first atom */
1011 if (dc->ncon+1 > dc->con_nalloc)
1013 dc->con_nalloc = over_alloc_large(dc->ncon+1);
1014 srenew(dc->con_gl, dc->con_nalloc);
1015 srenew(dc->con_nlocat, dc->con_nalloc);
1017 dc->con_gl[dc->ncon] = con_offset + con;
1018 dc->con_nlocat[dc->ncon] = 2;
1019 if (ilc_local->nr + 3 > ilc_local->nalloc)
1021 ilc_local->nalloc = over_alloc_dd(ilc_local->nr + 3);
1022 srenew(ilc_local->iatoms, ilc_local->nalloc);
1025 ilc_local->iatoms[ilc_local->nr++] = iap[0];
1026 ilc_local->iatoms[ilc_local->nr++] = (a_gl == iap[1] ? a : b_lo);
1027 ilc_local->iatoms[ilc_local->nr++] = (a_gl == iap[1] ? b_lo : a );
1034 /* We need the nrec constraints coupled to this constraint,
1035 * so we need to walk out of the home cell by nrec+1 atoms,
1036 * since already atom bg is not locally present.
1037 * Therefore we call walk_out with nrec recursions to go
1038 * after this first call.
1040 walk_out(con, con_offset, b_mol, offset, nrec,
1041 ncon1, ia1, ia2, at2con,
1042 dd->ga2la, TRUE, dc, dcc, ilc_local, ireq);
1049 gmx_mtop_atomlookup_destroy(alook);
1054 "Constraints: home %3d border %3d atoms: %3d\n",
1055 nhome, dc->ncon-nhome,
1056 dd->constraint_comm ? ireq->n : 0);
1060 int dd_make_local_constraints(gmx_domdec_t *dd, int at_start,
1061 const gmx_mtop_t *mtop,
1063 gmx_constr_t constr, int nrec,
1066 gmx_domdec_constraints_t *dc;
1067 t_ilist *ilc_local, *ils_local;
1069 const t_blocka *at2con_mt;
1070 const int **at2settle_mt;
1071 gmx_hash_t ga2la_specat;
1075 // This code should not be called unless this condition is true,
1076 // because that's the only time init_domdec_constraints is
1078 GMX_RELEASE_ASSERT(dd->bInterCGcons || dd->bInterCGsettles, "dd_make_local_constraints called when there are no local constraints");
1079 // ... and init_domdec_constraints always sets
1080 // dd->constraint_comm...
1081 GMX_RELEASE_ASSERT(dd->constraint_comm, "Invalid use of dd_make_local_constraints before construction of constraint_comm");
1082 // ... which static analysis needs to be reassured about, because
1083 // otherwise, when dd->bInterCGsettles is
1084 // true. dd->constraint_comm is unilaterally dereferenced before
1085 // the call to atoms_to_settles.
1087 dc = dd->constraints;
1089 ilc_local = &il_local[F_CONSTR];
1090 ils_local = &il_local[F_SETTLE];
1094 if (dd->constraint_comm)
1096 at2con_mt = atom2constraints_moltype(constr);
1097 ireq = &dd->constraint_comm->ireq[0];
1102 // Currently unreachable
1107 if (dd->bInterCGsettles)
1109 at2settle_mt = atom2settle_moltype(constr);
1114 /* Settle works inside charge groups, we assigned them already */
1115 at2settle_mt = NULL;
1118 if (at2settle_mt == NULL)
1120 atoms_to_constraints(dd, mtop, cginfo, at2con_mt, nrec,
1128 /* Do the constraints, if present, on the first thread.
1129 * Do the settles on all other threads.
1131 t0_set = ((at2con_mt != NULL && dc->nthread > 1) ? 1 : 0);
1133 #pragma omp parallel for num_threads(dc->nthread) schedule(static)
1134 for (thread = 0; thread < dc->nthread; thread++)
1136 if (at2con_mt && thread == 0)
1138 atoms_to_constraints(dd, mtop, cginfo, at2con_mt, nrec,
1142 if (thread >= t0_set)
1148 /* Distribute the settle check+assignments over
1149 * dc->nthread or dc->nthread-1 threads.
1151 cg0 = (dd->ncg_home*(thread-t0_set ))/(dc->nthread-t0_set);
1152 cg1 = (dd->ncg_home*(thread-t0_set+1))/(dc->nthread-t0_set);
1154 if (thread == t0_set)
1160 ilst = &dc->ils[thread];
1164 ireqt = &dd->constraint_comm->ireq[thread];
1170 atoms_to_settles(dd, mtop, cginfo, at2settle_mt,
1176 /* Combine the generate settles and requested indices */
1177 for (thread = 1; thread < dc->nthread; thread++)
1183 if (thread > t0_set)
1185 ilst = &dc->ils[thread];
1186 if (ils_local->nr + ilst->nr > ils_local->nalloc)
1188 ils_local->nalloc = over_alloc_large(ils_local->nr + ilst->nr);
1189 srenew(ils_local->iatoms, ils_local->nalloc);
1191 for (ia = 0; ia < ilst->nr; ia++)
1193 ils_local->iatoms[ils_local->nr+ia] = ilst->iatoms[ia];
1195 ils_local->nr += ilst->nr;
1198 ireqt = &dd->constraint_comm->ireq[thread];
1199 if (ireq->n+ireqt->n > ireq->nalloc)
1201 ireq->nalloc = over_alloc_large(ireq->n+ireqt->n);
1202 srenew(ireq->ind, ireq->nalloc);
1204 for (ia = 0; ia < ireqt->n; ia++)
1206 ireq->ind[ireq->n+ia] = ireqt->ind[ia];
1208 ireq->n += ireqt->n;
1213 fprintf(debug, "Settles: total %3d\n", ils_local->nr/4);
1217 if (dd->constraint_comm)
1222 setup_specat_communication(dd, ireq, dd->constraint_comm,
1223 dd->constraints->ga2la,
1225 "constraint", " or lincs-order");
1227 /* Fill in the missing indices */
1228 ga2la_specat = dd->constraints->ga2la;
1230 nral1 = 1 + NRAL(F_CONSTR);
1231 for (i = 0; i < ilc_local->nr; i += nral1)
1233 iap = ilc_local->iatoms + i;
1234 for (j = 1; j < nral1; j++)
1238 iap[j] = gmx_hash_get_minone(ga2la_specat, -iap[j]-1);
1243 nral1 = 1 + NRAL(F_SETTLE);
1244 for (i = 0; i < ils_local->nr; i += nral1)
1246 iap = ils_local->iatoms + i;
1247 for (j = 1; j < nral1; j++)
1251 iap[j] = gmx_hash_get_minone(ga2la_specat, -iap[j]-1);
1258 // Currently unreachable
1265 int dd_make_local_vsites(gmx_domdec_t *dd, int at_start, t_ilist *lil)
1267 gmx_domdec_specat_comm_t *spac;
1269 gmx_hash_t ga2la_specat;
1270 int ftype, nral, i, j, a;
1275 spac = dd->vsite_comm;
1276 ireq = &spac->ireq[0];
1277 ga2la_specat = dd->ga2la_vsite;
1280 /* Loop over all the home vsites */
1281 for (ftype = 0; ftype < F_NRE; ftype++)
1283 if (interaction_function[ftype].flags & IF_VSITE)
1287 for (i = 0; i < lilf->nr; i += 1+nral)
1289 iatoms = lilf->iatoms + i;
1290 /* Check if we have the other atoms */
1291 for (j = 1; j < 1+nral; j++)
1295 /* This is not a home atom,
1296 * we need to ask our neighbors.
1299 /* Check to not ask for the same atom more than once */
1300 if (gmx_hash_get_minone(dd->ga2la_vsite, a) == -1)
1302 /* Add this non-home atom to the list */
1303 if (ireq->n+1 > ireq->nalloc)
1305 ireq->nalloc = over_alloc_large(ireq->n+1);
1306 srenew(ireq->ind, ireq->nalloc);
1308 ireq->ind[ireq->n++] = a;
1309 /* Temporarily mark with -2,
1310 * we get the index later.
1312 gmx_hash_set(ga2la_specat, a, -2);
1320 at_end = setup_specat_communication(dd, ireq, dd->vsite_comm, ga2la_specat,
1321 at_start, 1, "vsite", "");
1323 /* Fill in the missing indices */
1324 for (ftype = 0; ftype < F_NRE; ftype++)
1326 if (interaction_function[ftype].flags & IF_VSITE)
1330 for (i = 0; i < lilf->nr; i += 1+nral)
1332 iatoms = lilf->iatoms + i;
1333 for (j = 1; j < 1+nral; j++)
1337 iatoms[j] = gmx_hash_get_minone(ga2la_specat, -iatoms[j]-1);
1347 static gmx_domdec_specat_comm_t *specat_comm_init(int nthread)
1349 gmx_domdec_specat_comm_t *spac;
1352 spac->nthread = nthread;
1353 snew(spac->ireq, spac->nthread);
1358 void init_domdec_constraints(gmx_domdec_t *dd,
1361 gmx_domdec_constraints_t *dc;
1362 gmx_molblock_t *molb;
1367 fprintf(debug, "Begin init_domdec_constraints\n");
1370 snew(dd->constraints, 1);
1371 dc = dd->constraints;
1373 snew(dc->molb_con_offset, mtop->nmolblock);
1374 snew(dc->molb_ncon_mol, mtop->nmolblock);
1377 for (mb = 0; mb < mtop->nmolblock; mb++)
1379 molb = &mtop->molblock[mb];
1380 dc->molb_con_offset[mb] = ncon;
1381 dc->molb_ncon_mol[mb] =
1382 mtop->moltype[molb->type].ilist[F_CONSTR].nr/3 +
1383 mtop->moltype[molb->type].ilist[F_CONSTRNC].nr/3;
1384 ncon += molb->nmol*dc->molb_ncon_mol[mb];
1389 snew(dc->gc_req, ncon);
1390 for (c = 0; c < ncon; c++)
1396 /* Use a hash table for the global to local index.
1397 * The number of keys is a rough estimate, it will be optimized later.
1399 dc->ga2la = gmx_hash_init(std::min(mtop->natoms/20,
1400 mtop->natoms/(2*dd->nnodes)));
1402 dc->nthread = gmx_omp_nthreads_get(emntDomdec);
1403 snew(dc->ils, dc->nthread);
1405 dd->constraint_comm = specat_comm_init(dc->nthread);
1408 void init_domdec_vsites(gmx_domdec_t *dd, int n_intercg_vsite)
1412 fprintf(debug, "Begin init_domdec_vsites\n");
1415 /* Use a hash table for the global to local index.
1416 * The number of keys is a rough estimate, it will be optimized later.
1418 dd->ga2la_vsite = gmx_hash_init(std::min(n_intercg_vsite/20,
1419 n_intercg_vsite/(2*dd->nnodes)));
1421 dd->vsite_comm = specat_comm_init(1);