1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This file is part of Gromacs Copyright (c) 1991-2008
5 * David van der Spoel, Erik Lindahl, Berk Hess, University of Groningen.
7 * This program is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU General Public License
9 * as published by the Free Software Foundation; either version 2
10 * of the License, or (at your option) any later version.
12 * To help us fund GROMACS development, we humbly ask that you cite
13 * the research papers on the package. Check out http://www.gromacs.org
16 * Gnomes, ROck Monsters And Chili Sauce
28 #include "domdec_network.h"
29 #include "mtop_util.h"
30 #include "gmx_ga2la.h"
32 #include "gmx_omp_nthreads.h"
48 typedef struct gmx_domdec_specat_comm {
49 /* The number of indices to receive during the setup */
51 /* The atoms to send */
52 gmx_specatsend_t spas[DIM][2];
62 /* The range in the local buffer(s) for received atoms */
66 /* The atom indices we need from the surrounding cells.
67 * We can gather the indices over nthread threads.
71 } gmx_domdec_specat_comm_t;
73 typedef struct gmx_domdec_constraints {
76 /* The fully local and connected constraints */
78 /* The global constraint number, only required for clearing gc_req */
82 /* Boolean that tells if a global constraint index has been requested */
84 /* Global to local communicated constraint atom only index */
87 /* Multi-threading stuff */
90 } gmx_domdec_constraints_t;
93 static void dd_move_f_specat(gmx_domdec_t *dd, gmx_domdec_specat_comm_t *spac,
94 rvec *f, rvec *fshift)
96 gmx_specatsend_t *spas;
98 int n, n0, n1, d, dim, dir, i;
101 gmx_bool bPBC, bScrew;
104 for (d = dd->ndim-1; d >= 0; d--)
109 /* Pulse the grid forward and backward */
110 spas = spac->spas[d];
115 /* Send and receive the coordinates */
116 dd_sendrecv2_rvec(dd, d,
117 f+n+n1, n0, vbuf, spas[0].nsend,
118 f+n, n1, vbuf+spas[0].nsend, spas[1].nsend);
119 for (dir = 0; dir < 2; dir++)
121 bPBC = ((dir == 0 && dd->ci[dim] == 0) ||
122 (dir == 1 && dd->ci[dim] == dd->nc[dim]-1));
123 bScrew = (bPBC && dd->bScrewPBC && dim == XX);
125 spas = &spac->spas[d][dir];
126 /* Sum the buffer into the required forces */
127 if (!bPBC || (!bScrew && fshift == NULL))
129 for (i = 0; i < spas->nsend; i++)
131 rvec_inc(f[spas->a[i]], *vbuf);
138 vis[dim] = (dir == 0 ? 1 : -1);
142 /* Sum and add to shift forces */
143 for (i = 0; i < spas->nsend; i++)
145 rvec_inc(f[spas->a[i]], *vbuf);
146 rvec_inc(fshift[is], *vbuf);
152 /* Rotate the forces */
153 for (i = 0; i < spas->nsend; i++)
155 f[spas->a[i]][XX] += (*vbuf)[XX];
156 f[spas->a[i]][YY] -= (*vbuf)[YY];
157 f[spas->a[i]][ZZ] -= (*vbuf)[ZZ];
160 rvec_inc(fshift[is], *vbuf);
170 /* Two cells, so we only need to communicate one way */
171 spas = &spac->spas[d][0];
173 /* Send and receive the coordinates */
174 dd_sendrecv_rvec(dd, d, dddirForward,
175 f+n, spas->nrecv, spac->vbuf, spas->nsend);
176 /* Sum the buffer into the required forces */
177 if (dd->bScrewPBC && dim == XX &&
179 dd->ci[dim] == dd->nc[dim]-1))
181 for (i = 0; i < spas->nsend; i++)
183 /* Rotate the force */
184 f[spas->a[i]][XX] += spac->vbuf[i][XX];
185 f[spas->a[i]][YY] -= spac->vbuf[i][YY];
186 f[spas->a[i]][ZZ] -= spac->vbuf[i][ZZ];
191 for (i = 0; i < spas->nsend; i++)
193 rvec_inc(f[spas->a[i]], spac->vbuf[i]);
200 void dd_move_f_vsites(gmx_domdec_t *dd, rvec *f, rvec *fshift)
204 dd_move_f_specat(dd, dd->vsite_comm, f, fshift);
208 void dd_clear_f_vsites(gmx_domdec_t *dd, rvec *f)
214 for (i = dd->vsite_comm->at_start; i < dd->vsite_comm->at_end; i++)
221 static void dd_move_x_specat(gmx_domdec_t *dd, gmx_domdec_specat_comm_t *spac,
222 matrix box, rvec *x0, rvec *x1)
224 gmx_specatsend_t *spas;
225 rvec *x, *vbuf, *rbuf;
226 int nvec, v, n, nn, ns0, ns1, nr0, nr1, nr, d, dim, dir, i;
227 gmx_bool bPBC, bScrew = FALSE;
228 rvec shift = {0, 0, 0};
237 for (d = 0; d < dd->ndim; d++)
242 /* Pulse the grid forward and backward */
244 for (dir = 0; dir < 2; dir++)
246 if (dir == 0 && dd->ci[dim] == 0)
249 bScrew = (dd->bScrewPBC && dim == XX);
250 copy_rvec(box[dim], shift);
252 else if (dir == 1 && dd->ci[dim] == dd->nc[dim]-1)
255 bScrew = (dd->bScrewPBC && dim == XX);
256 for (i = 0; i < DIM; i++)
258 shift[i] = -box[dim][i];
266 spas = &spac->spas[d][dir];
267 for (v = 0; v < nvec; v++)
269 x = (v == 0 ? x0 : x1);
270 /* Copy the required coordinates to the send buffer */
274 for (i = 0; i < spas->nsend; i++)
276 copy_rvec(x[spas->a[i]], *vbuf);
282 /* Shift coordinates */
283 for (i = 0; i < spas->nsend; i++)
285 rvec_add(x[spas->a[i]], shift, *vbuf);
291 /* Shift and rotate coordinates */
292 for (i = 0; i < spas->nsend; i++)
294 (*vbuf)[XX] = x[spas->a[i]][XX] + shift[XX];
295 (*vbuf)[YY] = box[YY][YY] - x[spas->a[i]][YY] + shift[YY];
296 (*vbuf)[ZZ] = box[ZZ][ZZ] - x[spas->a[i]][ZZ] + shift[ZZ];
302 /* Send and receive the coordinates */
303 spas = spac->spas[d];
310 dd_sendrecv2_rvec(dd, d,
311 spac->vbuf+ns0, ns1, x0+n, nr1,
312 spac->vbuf, ns0, x0+n+nr1, nr0);
316 /* Communicate both vectors in one buffer */
318 dd_sendrecv2_rvec(dd, d,
319 spac->vbuf+2*ns0, 2*ns1, rbuf, 2*nr1,
320 spac->vbuf, 2*ns0, rbuf+2*nr1, 2*nr0);
321 /* Split the buffer into the two vectors */
323 for (dir = 1; dir >= 0; dir--)
325 nr = spas[dir].nrecv;
326 for (v = 0; v < 2; v++)
328 x = (v == 0 ? x0 : x1);
329 for (i = 0; i < nr; i++)
331 copy_rvec(*rbuf, x[nn+i]);
342 spas = &spac->spas[d][0];
343 /* Copy the required coordinates to the send buffer */
345 for (v = 0; v < nvec; v++)
347 x = (v == 0 ? x0 : x1);
348 if (dd->bScrewPBC && dim == XX &&
349 (dd->ci[XX] == 0 || dd->ci[XX] == dd->nc[XX]-1))
351 /* Here we only perform the rotation, the rest of the pbc
352 * is handled in the constraint or viste routines.
354 for (i = 0; i < spas->nsend; i++)
356 (*vbuf)[XX] = x[spas->a[i]][XX];
357 (*vbuf)[YY] = box[YY][YY] - x[spas->a[i]][YY];
358 (*vbuf)[ZZ] = box[ZZ][ZZ] - x[spas->a[i]][ZZ];
364 for (i = 0; i < spas->nsend; i++)
366 copy_rvec(x[spas->a[i]], *vbuf);
371 /* Send and receive the coordinates */
374 dd_sendrecv_rvec(dd, d, dddirBackward,
375 spac->vbuf, spas->nsend, x0+n, spas->nrecv);
379 /* Communicate both vectors in one buffer */
381 dd_sendrecv_rvec(dd, d, dddirBackward,
382 spac->vbuf, 2*spas->nsend, rbuf, 2*spas->nrecv);
383 /* Split the buffer into the two vectors */
385 for (v = 0; v < 2; v++)
387 x = (v == 0 ? x0 : x1);
388 for (i = 0; i < nr; i++)
390 copy_rvec(*rbuf, x[n+i]);
400 void dd_move_x_constraints(gmx_domdec_t *dd, matrix box, rvec *x0, rvec *x1)
402 if (dd->constraint_comm)
404 dd_move_x_specat(dd, dd->constraint_comm, box, x0, x1);
408 void dd_move_x_vsites(gmx_domdec_t *dd, matrix box, rvec *x)
412 dd_move_x_specat(dd, dd->vsite_comm, box, x, NULL);
416 int *dd_constraints_nlocalatoms(gmx_domdec_t *dd)
420 return dd->constraints->con_nlocat;
428 void dd_clear_local_constraint_indices(gmx_domdec_t *dd)
430 gmx_domdec_constraints_t *dc;
433 dc = dd->constraints;
435 for (i = 0; i < dc->ncon; i++)
437 dc->gc_req[dc->con_gl[i]] = 0;
440 if (dd->constraint_comm)
442 gmx_hash_clear_and_optimize(dc->ga2la);
446 void dd_clear_local_vsite_indices(gmx_domdec_t *dd)
452 gmx_hash_clear_and_optimize(dd->ga2la_vsite);
456 static int setup_specat_communication(gmx_domdec_t *dd,
458 gmx_domdec_specat_comm_t *spac,
459 gmx_hash_t ga2la_specat,
462 const char *specat_type,
465 int nsend[2], nlast, nsend_zero[2] = {0, 0}, *nsend_ptr;
466 int d, dim, ndir, dir, nr, ns, i, nrecv_local, n0, start, indr, ind, buf[2];
467 int nat_tot_specat, nat_tot_prev, nalloc_old;
468 gmx_bool bPBC, bFirst;
469 gmx_specatsend_t *spas;
473 fprintf(debug, "Begin setup_specat_communication for %s\n", specat_type);
476 /* nsend[0]: the number of atoms requested by this node only,
477 * we communicate this for more efficients checks
478 * nsend[1]: the total number of requested atoms
483 for (d = dd->ndim-1; d >= 0; d--)
485 /* Pulse the grid forward and backward */
487 bPBC = (dim < dd->npbcdim);
488 if (dd->nc[dim] == 2)
490 /* Only 2 cells, so we only need to communicate once */
497 for (dir = 0; dir < ndir; dir++)
501 ((dir == 0 && dd->ci[dim] == dd->nc[dim] - 1) ||
502 (dir == 1 && dd->ci[dim] == 0)))
504 /* No pbc: the fist/last cell should not request atoms */
505 nsend_ptr = nsend_zero;
511 /* Communicate the number of indices */
512 dd_sendrecv_int(dd, d, dir == 0 ? dddirForward : dddirBackward,
513 nsend_ptr, 2, spac->nreq[d][dir], 2);
514 nr = spac->nreq[d][dir][1];
515 if (nlast+nr > ireq->nalloc)
517 ireq->nalloc = over_alloc_dd(nlast+nr);
518 srenew(ireq->ind, ireq->nalloc);
520 /* Communicate the indices */
521 dd_sendrecv_int(dd, d, dir == 0 ? dddirForward : dddirBackward,
522 ireq->ind, nsend_ptr[1], ireq->ind+nlast, nr);
529 fprintf(debug, "Communicated the counts\n");
532 /* Search for the requested atoms and communicate the indices we have */
533 nat_tot_specat = at_start;
535 for (d = 0; d < dd->ndim; d++)
538 /* Pulse the grid forward and backward */
539 if (dd->dim[d] >= dd->npbcdim || dd->nc[dd->dim[d]] > 2)
547 nat_tot_prev = nat_tot_specat;
548 for (dir = ndir-1; dir >= 0; dir--)
550 if (nat_tot_specat > spac->bSendAtom_nalloc)
552 nalloc_old = spac->bSendAtom_nalloc;
553 spac->bSendAtom_nalloc = over_alloc_dd(nat_tot_specat);
554 srenew(spac->bSendAtom, spac->bSendAtom_nalloc);
555 for (i = nalloc_old; i < spac->bSendAtom_nalloc; i++)
557 spac->bSendAtom[i] = FALSE;
560 spas = &spac->spas[d][dir];
561 n0 = spac->nreq[d][dir][0];
562 nr = spac->nreq[d][dir][1];
565 fprintf(debug, "dim=%d, dir=%d, searching for %d atoms\n",
571 for (i = 0; i < nr; i++)
573 indr = ireq->ind[start+i];
575 /* Check if this is a home atom and if so ind will be set */
576 if (!ga2la_get_home(dd->ga2la, indr, &ind))
578 /* Search in the communicated atoms */
579 ind = gmx_hash_get_minone(ga2la_specat, indr);
583 if (i < n0 || !spac->bSendAtom[ind])
585 if (spas->nsend+1 > spas->a_nalloc)
587 spas->a_nalloc = over_alloc_large(spas->nsend+1);
588 srenew(spas->a, spas->a_nalloc);
590 /* Store the local index so we know which coordinates
593 spas->a[spas->nsend] = ind;
594 spac->bSendAtom[ind] = TRUE;
595 if (spas->nsend+1 > spac->ibuf_nalloc)
597 spac->ibuf_nalloc = over_alloc_large(spas->nsend+1);
598 srenew(spac->ibuf, spac->ibuf_nalloc);
600 /* Store the global index so we can send it now */
601 spac->ibuf[spas->nsend] = indr;
611 /* Clear the local flags */
612 for (i = 0; i < spas->nsend; i++)
614 spac->bSendAtom[spas->a[i]] = FALSE;
616 /* Send and receive the number of indices to communicate */
617 nsend[1] = spas->nsend;
618 dd_sendrecv_int(dd, d, dir == 0 ? dddirBackward : dddirForward,
622 fprintf(debug, "Send to node %d, %d (%d) indices, "
623 "receive from node %d, %d (%d) indices\n",
624 dd->neighbor[d][1-dir], nsend[1], nsend[0],
625 dd->neighbor[d][dir], buf[1], buf[0]);
628 for (i = 0; i < spas->nsend; i++)
630 fprintf(debug, " %d", spac->ibuf[i]+1);
632 fprintf(debug, "\n");
635 nrecv_local += buf[0];
636 spas->nrecv = buf[1];
637 if (nat_tot_specat + spas->nrecv > dd->gatindex_nalloc)
639 dd->gatindex_nalloc =
640 over_alloc_dd(nat_tot_specat + spas->nrecv);
641 srenew(dd->gatindex, dd->gatindex_nalloc);
643 /* Send and receive the indices */
644 dd_sendrecv_int(dd, d, dir == 0 ? dddirBackward : dddirForward,
645 spac->ibuf, spas->nsend,
646 dd->gatindex+nat_tot_specat, spas->nrecv);
647 nat_tot_specat += spas->nrecv;
650 /* Allocate the x/f communication buffers */
651 ns = spac->spas[d][0].nsend;
652 nr = spac->spas[d][0].nrecv;
655 ns += spac->spas[d][1].nsend;
656 nr += spac->spas[d][1].nrecv;
658 if (vbuf_fac*ns > spac->vbuf_nalloc)
660 spac->vbuf_nalloc = over_alloc_dd(vbuf_fac*ns);
661 srenew(spac->vbuf, spac->vbuf_nalloc);
663 if (vbuf_fac == 2 && vbuf_fac*nr > spac->vbuf2_nalloc)
665 spac->vbuf2_nalloc = over_alloc_dd(vbuf_fac*nr);
666 srenew(spac->vbuf2, spac->vbuf2_nalloc);
669 /* Make a global to local index for the communication atoms */
670 for (i = nat_tot_prev; i < nat_tot_specat; i++)
672 gmx_hash_change_or_set(ga2la_specat, dd->gatindex[i], i);
676 /* Check that in the end we got the number of atoms we asked for */
677 if (nrecv_local != ireq->n)
681 fprintf(debug, "Requested %d, received %d (tot recv %d)\n",
682 ireq->n, nrecv_local, nat_tot_specat-at_start);
685 for (i = 0; i < ireq->n; i++)
687 ind = gmx_hash_get_minone(ga2la_specat, ireq->ind[i]);
688 fprintf(debug, " %s%d",
689 (ind >= 0) ? "" : "!",
692 fprintf(debug, "\n");
695 fprintf(stderr, "\nDD cell %d %d %d: Neighboring cells do not have atoms:",
696 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
697 for (i = 0; i < ireq->n; i++)
699 if (gmx_hash_get_minone(ga2la_specat, ireq->ind[i]) < 0)
701 fprintf(stderr, " %d", ireq->ind[i]+1);
704 fprintf(stderr, "\n");
705 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.",
706 dd->ci[XX], dd->ci[YY], dd->ci[ZZ],
707 nrecv_local, ireq->n, specat_type,
708 specat_type, add_err,
709 dd->bGridJump ? " or use the -rcon option of mdrun" : "");
712 spac->at_start = at_start;
713 spac->at_end = nat_tot_specat;
717 fprintf(debug, "Done setup_specat_communication\n");
720 return nat_tot_specat;
723 static void walk_out(int con, int con_offset, int a, int offset, int nrec,
724 int ncon1, const t_iatom *ia1, const t_iatom *ia2,
725 const t_blocka *at2con,
726 const gmx_ga2la_t ga2la, gmx_bool bHomeConnect,
727 gmx_domdec_constraints_t *dc,
728 gmx_domdec_specat_comm_t *dcc,
732 int a1_gl, a2_gl, a_loc, i, coni, b;
735 if (dc->gc_req[con_offset+con] == 0)
737 /* Add this non-home constraint to the list */
738 if (dc->ncon+1 > dc->con_nalloc)
740 dc->con_nalloc = over_alloc_large(dc->ncon+1);
741 srenew(dc->con_gl, dc->con_nalloc);
742 srenew(dc->con_nlocat, dc->con_nalloc);
744 dc->con_gl[dc->ncon] = con_offset + con;
745 dc->con_nlocat[dc->ncon] = (bHomeConnect ? 1 : 0);
746 dc->gc_req[con_offset+con] = 1;
747 if (il_local->nr + 3 > il_local->nalloc)
749 il_local->nalloc = over_alloc_dd(il_local->nr+3);
750 srenew(il_local->iatoms, il_local->nalloc);
752 iap = constr_iatomptr(ncon1, ia1, ia2, con);
753 il_local->iatoms[il_local->nr++] = iap[0];
754 a1_gl = offset + iap[1];
755 a2_gl = offset + iap[2];
756 /* The following indexing code can probably be optizimed */
757 if (ga2la_get_home(ga2la, a1_gl, &a_loc))
759 il_local->iatoms[il_local->nr++] = a_loc;
763 /* We set this index later */
764 il_local->iatoms[il_local->nr++] = -a1_gl - 1;
766 if (ga2la_get_home(ga2la, a2_gl, &a_loc))
768 il_local->iatoms[il_local->nr++] = a_loc;
772 /* We set this index later */
773 il_local->iatoms[il_local->nr++] = -a2_gl - 1;
777 /* Check to not ask for the same atom more than once */
778 if (gmx_hash_get_minone(dc->ga2la, offset+a) == -1)
781 /* Add this non-home atom to the list */
782 if (ireq->n+1 > ireq->nalloc)
784 ireq->nalloc = over_alloc_large(ireq->n+1);
785 srenew(ireq->ind, ireq->nalloc);
787 ireq->ind[ireq->n++] = offset + a;
788 /* Temporarily mark with -2, we get the index later */
789 gmx_hash_set(dc->ga2la, offset+a, -2);
794 for (i = at2con->index[a]; i < at2con->index[a+1]; i++)
800 iap = constr_iatomptr(ncon1, ia1, ia2, coni);
809 if (!ga2la_get_home(ga2la, offset+b, &a_loc))
811 walk_out(coni, con_offset, b, offset, nrec-1,
812 ncon1, ia1, ia2, at2con,
813 ga2la, FALSE, dc, dcc, il_local, ireq);
820 static void atoms_to_settles(gmx_domdec_t *dd,
821 const gmx_mtop_t *mtop,
823 const int **at2settle_mt,
824 int cg_start, int cg_end,
829 gmx_mtop_atomlookup_t alook;
832 int cg, a, a_gl, a_glsa, a_gls[3], a_locs[3];
833 int mb, molnr, a_mol, offset;
834 const gmx_molblock_t *molb;
842 alook = gmx_mtop_atomlookup_settle_init(mtop);
844 nral = NRAL(F_SETTLE);
846 for (cg = cg_start; cg < cg_end; cg++)
848 if (GET_CGINFO_SETTLE(cginfo[cg]))
850 for (a = dd->cgindex[cg]; a < dd->cgindex[cg+1]; a++)
852 a_gl = dd->gatindex[a];
854 gmx_mtop_atomnr_to_molblock_ind(alook, a_gl, &mb, &molnr, &a_mol);
855 molb = &mtop->molblock[mb];
857 settle = at2settle_mt[molb->type][a_mol];
861 offset = a_gl - a_mol;
863 ia1 = mtop->moltype[molb->type].ilist[F_SETTLE].iatoms;
867 for (sa = 0; sa < nral; sa++)
869 a_glsa = offset + ia1[settle*(1+nral)+1+sa];
871 a_home[sa] = ga2la_get_home(ga2la, a_glsa, &a_locs[sa]);
874 if (nlocal == 0 && a_gl == a_glsa)
884 if (ils_local->nr+1+nral > ils_local->nalloc)
886 ils_local->nalloc = over_alloc_dd(ils_local->nr+1+nral);
887 srenew(ils_local->iatoms, ils_local->nalloc);
890 ils_local->iatoms[ils_local->nr++] = ia1[settle*4];
892 for (sa = 0; sa < nral; sa++)
894 if (ga2la_get_home(ga2la, a_gls[sa], &a_locs[sa]))
896 ils_local->iatoms[ils_local->nr++] = a_locs[sa];
900 ils_local->iatoms[ils_local->nr++] = -a_gls[sa] - 1;
901 /* Add this non-home atom to the list */
902 if (ireq->n+1 > ireq->nalloc)
904 ireq->nalloc = over_alloc_large(ireq->n+1);
905 srenew(ireq->ind, ireq->nalloc);
907 ireq->ind[ireq->n++] = a_gls[sa];
908 /* A check on double atom requests is
909 * not required for settle.
919 gmx_mtop_atomlookup_destroy(alook);
922 static void atoms_to_constraints(gmx_domdec_t *dd,
923 const gmx_mtop_t *mtop,
925 const t_blocka *at2con_mt, int nrec,
929 const t_blocka *at2con;
931 gmx_mtop_atomlookup_t alook;
933 gmx_molblock_t *molb;
934 t_iatom *ia1, *ia2, *iap;
935 int nhome, cg, a, a_gl, a_mol, a_loc, b_lo, offset, mb, molnr, b_mol, i, con, con_offset;
936 gmx_domdec_constraints_t *dc;
937 gmx_domdec_specat_comm_t *dcc;
939 dc = dd->constraints;
940 dcc = dd->constraint_comm;
944 alook = gmx_mtop_atomlookup_init(mtop);
947 for (cg = 0; cg < dd->ncg_home; cg++)
949 if (GET_CGINFO_CONSTR(cginfo[cg]))
951 for (a = dd->cgindex[cg]; a < dd->cgindex[cg+1]; a++)
953 a_gl = dd->gatindex[a];
955 gmx_mtop_atomnr_to_molblock_ind(alook, a_gl, &mb, &molnr, &a_mol);
956 molb = &mtop->molblock[mb];
958 ncon1 = mtop->moltype[molb->type].ilist[F_CONSTR].nr/NRAL(F_SETTLE);
960 ia1 = mtop->moltype[molb->type].ilist[F_CONSTR].iatoms;
961 ia2 = mtop->moltype[molb->type].ilist[F_CONSTRNC].iatoms;
963 /* Calculate the global constraint number offset for the molecule.
964 * This is only required for the global index to make sure
965 * that we use each constraint only once.
968 dc->molb_con_offset[mb] + molnr*dc->molb_ncon_mol[mb];
970 /* The global atom number offset for this molecule */
971 offset = a_gl - a_mol;
972 at2con = &at2con_mt[molb->type];
973 for (i = at2con->index[a_mol]; i < at2con->index[a_mol+1]; i++)
976 iap = constr_iatomptr(ncon1, ia1, ia2, con);
985 if (ga2la_get_home(ga2la, offset+b_mol, &a_loc))
987 /* Add this fully home constraint at the first atom */
990 if (dc->ncon+1 > dc->con_nalloc)
992 dc->con_nalloc = over_alloc_large(dc->ncon+1);
993 srenew(dc->con_gl, dc->con_nalloc);
994 srenew(dc->con_nlocat, dc->con_nalloc);
996 dc->con_gl[dc->ncon] = con_offset + con;
997 dc->con_nlocat[dc->ncon] = 2;
998 if (ilc_local->nr + 3 > ilc_local->nalloc)
1000 ilc_local->nalloc = over_alloc_dd(ilc_local->nr + 3);
1001 srenew(ilc_local->iatoms, ilc_local->nalloc);
1004 ilc_local->iatoms[ilc_local->nr++] = iap[0];
1005 ilc_local->iatoms[ilc_local->nr++] = (a_gl == iap[1] ? a : b_lo);
1006 ilc_local->iatoms[ilc_local->nr++] = (a_gl == iap[1] ? b_lo : a );
1013 /* We need the nrec constraints coupled to this constraint,
1014 * so we need to walk out of the home cell by nrec+1 atoms,
1015 * since already atom bg is not locally present.
1016 * Therefore we call walk_out with nrec recursions to go
1017 * after this first call.
1019 walk_out(con, con_offset, b_mol, offset, nrec,
1020 ncon1, ia1, ia2, at2con,
1021 dd->ga2la, TRUE, dc, dcc, ilc_local, ireq);
1028 gmx_mtop_atomlookup_destroy(alook);
1033 "Constraints: home %3d border %3d atoms: %3d\n",
1034 nhome, dc->ncon-nhome,
1035 dd->constraint_comm ? ireq->n : 0);
1039 int dd_make_local_constraints(gmx_domdec_t *dd, int at_start,
1040 const gmx_mtop_t *mtop,
1042 gmx_constr_t constr, int nrec,
1045 gmx_domdec_constraints_t *dc;
1046 t_ilist *ilc_local, *ils_local;
1048 const t_blocka *at2con_mt;
1049 const int **at2settle_mt;
1050 gmx_hash_t ga2la_specat;
1054 dc = dd->constraints;
1056 ilc_local = &il_local[F_CONSTR];
1057 ils_local = &il_local[F_SETTLE];
1061 if (dd->constraint_comm)
1063 at2con_mt = atom2constraints_moltype(constr);
1064 ireq = &dd->constraint_comm->ireq[0];
1073 if (dd->bInterCGsettles)
1075 at2settle_mt = atom2settle_moltype(constr);
1080 /* Settle works inside charge groups, we assigned them already */
1081 at2settle_mt = NULL;
1084 if (at2settle_mt == NULL)
1086 atoms_to_constraints(dd, mtop, cginfo, at2con_mt, nrec,
1094 /* Do the constraints, if present, on the first thread.
1095 * Do the settles on all other threads.
1097 t0_set = ((at2con_mt != NULL && dc->nthread > 1) ? 1 : 0);
1099 #pragma omp parallel for num_threads(dc->nthread) schedule(static)
1100 for (thread = 0; thread < dc->nthread; thread++)
1102 if (at2con_mt && thread == 0)
1104 atoms_to_constraints(dd, mtop, cginfo, at2con_mt, nrec,
1108 if (thread >= t0_set)
1114 /* Distribute the settle check+assignments over
1115 * dc->nthread or dc->nthread-1 threads.
1117 cg0 = (dd->ncg_home*(thread-t0_set ))/(dc->nthread-t0_set);
1118 cg1 = (dd->ncg_home*(thread-t0_set+1))/(dc->nthread-t0_set);
1120 if (thread == t0_set)
1126 ilst = &dc->ils[thread];
1130 ireqt = &dd->constraint_comm->ireq[thread];
1136 atoms_to_settles(dd, mtop, cginfo, at2settle_mt,
1142 /* Combine the generate settles and requested indices */
1143 for (thread = 1; thread < dc->nthread; thread++)
1149 if (thread > t0_set)
1151 ilst = &dc->ils[thread];
1152 if (ils_local->nr + ilst->nr > ils_local->nalloc)
1154 ils_local->nalloc = over_alloc_large(ils_local->nr + ilst->nr);
1155 srenew(ils_local->iatoms, ils_local->nalloc);
1157 for (ia = 0; ia < ilst->nr; ia++)
1159 ils_local->iatoms[ils_local->nr+ia] = ilst->iatoms[ia];
1161 ils_local->nr += ilst->nr;
1164 ireqt = &dd->constraint_comm->ireq[thread];
1165 if (ireq->n+ireqt->n > ireq->nalloc)
1167 ireq->nalloc = over_alloc_large(ireq->n+ireqt->n);
1168 srenew(ireq->ind, ireq->nalloc);
1170 for (ia = 0; ia < ireqt->n; ia++)
1172 ireq->ind[ireq->n+ia] = ireqt->ind[ia];
1174 ireq->n += ireqt->n;
1179 fprintf(debug, "Settles: total %3d\n", ils_local->nr/4);
1183 if (dd->constraint_comm)
1188 setup_specat_communication(dd, ireq, dd->constraint_comm,
1189 dd->constraints->ga2la,
1191 "constraint", " or lincs-order");
1193 /* Fill in the missing indices */
1194 ga2la_specat = dd->constraints->ga2la;
1196 nral1 = 1 + NRAL(F_CONSTR);
1197 for (i = 0; i < ilc_local->nr; i += nral1)
1199 iap = ilc_local->iatoms + i;
1200 for (j = 1; j < nral1; j++)
1204 iap[j] = gmx_hash_get_minone(ga2la_specat, -iap[j]-1);
1209 nral1 = 1 + NRAL(F_SETTLE);
1210 for (i = 0; i < ils_local->nr; i += nral1)
1212 iap = ils_local->iatoms + i;
1213 for (j = 1; j < nral1; j++)
1217 iap[j] = gmx_hash_get_minone(ga2la_specat, -iap[j]-1);
1230 int dd_make_local_vsites(gmx_domdec_t *dd, int at_start, t_ilist *lil)
1232 gmx_domdec_specat_comm_t *spac;
1234 gmx_hash_t ga2la_specat;
1235 int ftype, nral, i, j, gat, a;
1240 spac = dd->vsite_comm;
1241 ireq = &spac->ireq[0];
1242 ga2la_specat = dd->ga2la_vsite;
1245 /* Loop over all the home vsites */
1246 for (ftype = 0; ftype < F_NRE; ftype++)
1248 if (interaction_function[ftype].flags & IF_VSITE)
1252 for (i = 0; i < lilf->nr; i += 1+nral)
1254 iatoms = lilf->iatoms + i;
1255 /* Check if we have the other atoms */
1256 for (j = 1; j < 1+nral; j++)
1260 /* This is not a home atom,
1261 * we need to ask our neighbors.
1264 /* Check to not ask for the same atom more than once */
1265 if (gmx_hash_get_minone(dd->ga2la_vsite, a) == -1)
1267 /* Add this non-home atom to the list */
1268 if (ireq->n+1 > ireq->nalloc)
1270 ireq->nalloc = over_alloc_large(ireq->n+1);
1271 srenew(ireq->ind, ireq->nalloc);
1273 ireq->ind[ireq->n++] = a;
1274 /* Temporarily mark with -2,
1275 * we get the index later.
1277 gmx_hash_set(ga2la_specat, a, -2);
1285 at_end = setup_specat_communication(dd, ireq, dd->vsite_comm, ga2la_specat,
1286 at_start, 1, "vsite", "");
1288 /* Fill in the missing indices */
1289 for (ftype = 0; ftype < F_NRE; ftype++)
1291 if (interaction_function[ftype].flags & IF_VSITE)
1295 for (i = 0; i < lilf->nr; i += 1+nral)
1297 iatoms = lilf->iatoms + i;
1298 for (j = 1; j < 1+nral; j++)
1302 iatoms[j] = gmx_hash_get_minone(ga2la_specat, -iatoms[j]-1);
1312 static gmx_domdec_specat_comm_t *specat_comm_init(int nthread)
1314 gmx_domdec_specat_comm_t *spac;
1317 spac->nthread = nthread;
1318 snew(spac->ireq, spac->nthread);
1323 void init_domdec_constraints(gmx_domdec_t *dd,
1325 gmx_constr_t constr)
1327 gmx_domdec_constraints_t *dc;
1328 gmx_molblock_t *molb;
1333 fprintf(debug, "Begin init_domdec_constraints\n");
1336 snew(dd->constraints, 1);
1337 dc = dd->constraints;
1339 snew(dc->molb_con_offset, mtop->nmolblock);
1340 snew(dc->molb_ncon_mol, mtop->nmolblock);
1343 for (mb = 0; mb < mtop->nmolblock; mb++)
1345 molb = &mtop->molblock[mb];
1346 dc->molb_con_offset[mb] = ncon;
1347 dc->molb_ncon_mol[mb] =
1348 mtop->moltype[molb->type].ilist[F_CONSTR].nr/3 +
1349 mtop->moltype[molb->type].ilist[F_CONSTRNC].nr/3;
1350 ncon += molb->nmol*dc->molb_ncon_mol[mb];
1355 snew(dc->gc_req, ncon);
1356 for (c = 0; c < ncon; c++)
1362 /* Use a hash table for the global to local index.
1363 * The number of keys is a rough estimate, it will be optimized later.
1365 dc->ga2la = gmx_hash_init(min(mtop->natoms/20,
1366 mtop->natoms/(2*dd->nnodes)));
1368 dc->nthread = gmx_omp_nthreads_get(emntDomdec);
1369 snew(dc->ils, dc->nthread);
1371 dd->constraint_comm = specat_comm_init(dc->nthread);
1374 void init_domdec_vsites(gmx_domdec_t *dd, int n_intercg_vsite)
1377 gmx_domdec_constraints_t *dc;
1381 fprintf(debug, "Begin init_domdec_vsites\n");
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.
1387 dd->ga2la_vsite = gmx_hash_init(min(n_intercg_vsite/20,
1388 n_intercg_vsite/(2*dd->nnodes)));
1390 dd->vsite_comm = specat_comm_init(1);