src/gromacs/mdlib/domdec_con.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *
   4  * This file is part of Gromacs        Copyright (c) 1991-2008
   5  * David van der Spoel, Erik Lindahl, Berk Hess, University of Groningen.
   6  *
   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.
  11  *
  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
  14  *
  15  * And Hey:
  16  * Gnomes, ROck Monsters And Chili Sauce
  17  */
  18
  19 #ifdef HAVE_CONFIG_H
  20 #include <config.h>
  21 #endif
  22 #include <assert.h>
  23
  24 #include "smalloc.h"
  25 #include "vec.h"
  26 #include "constr.h"
  27 #include "domdec.h"
  28 #include "domdec_network.h"
  29 #include "mtop_util.h"
  30 #include "gmx_ga2la.h"
  31
  32 typedef struct {
  33     int nsend;
  34     int *a;
  35     int a_nalloc;
  36     int nrecv;
  37 } gmx_specatsend_t;
  38
  39 typedef struct gmx_domdec_specat_comm {
  40     /* The atom indices we need from the surrounding cells */
  41     int  nind_req;
  42     int  *ind_req;
  43     int  ind_req_nalloc;
  44     /* The number of indices to receive during the setup */
  45     int  nreq[DIM][2][2];
  46     /* The atoms to send */
  47     gmx_specatsend_t spas[DIM][2];
  48     gmx_bool *bSendAtom;
  49     int   bSendAtom_nalloc;
  50     /* Send buffers */
  51     int  *ibuf;
  52     int  ibuf_nalloc;
  53     rvec *vbuf;
  54     int  vbuf_nalloc;
  55     rvec *vbuf2;
  56     int  vbuf2_nalloc;
  57     /* The range in the local buffer(s) for received atoms */
  58     int  at_start;
  59     int  at_end;
  60 } gmx_domdec_specat_comm_t;
  61
  62 typedef struct gmx_domdec_constraints {
  63     int  *molb_con_offset;
  64     int  *molb_ncon_mol;
  65     /* The fully local and connected constraints */
  66     int  ncon;
  67     /* The global constraint number, only required for clearing gc_req */
  68     int  *con_gl;
  69     int  *con_nlocat;
  70     int  con_nalloc;
  71     /* Boolean that tells if a global constraint index has been requested */
  72     char *gc_req;
  73     /* Global to local communicated constraint atom only index */
  74     int  *ga2la;
  75 } gmx_domdec_constraints_t;
  76
  77
  78 static void dd_move_f_specat(gmx_domdec_t *dd,gmx_domdec_specat_comm_t *spac,
  79                              rvec *f,rvec *fshift)
  80 {
  81     gmx_specatsend_t *spas;
  82     rvec *vbuf;
  83     int  n,n0,n1,d,dim,dir,i;
  84     ivec vis;
  85     int  is;
  86     gmx_bool bPBC,bScrew;
  87
  88     n = spac->at_end;
  89     for(d=dd->ndim-1; d>=0; d--)
  90     {
  91         dim = dd->dim[d];
  92         if (dd->nc[dim] > 2)
  93         {
  94             /* Pulse the grid forward and backward */
  95             spas = spac->spas[d];
  96             n0 = spas[0].nrecv;
  97             n1 = spas[1].nrecv;
  98             n -= n1 + n0;
  99             vbuf = spac->vbuf;
 100             /* Send and receive the coordinates */
 101             dd_sendrecv2_rvec(dd,d,
 102                               f+n+n1,n0,vbuf              ,spas[0].nsend,
 103                               f+n   ,n1,vbuf+spas[0].nsend,spas[1].nsend);
 104             for(dir=0; dir<2; dir++)
 105             {
 106                 bPBC   = ((dir == 0 && dd->ci[dim] == 0) ||
 107                           (dir == 1 && dd->ci[dim] == dd->nc[dim]-1));
 108                 bScrew = (bPBC && dd->bScrewPBC && dim == XX);
 109
 110                 spas = &spac->spas[d][dir];
 111                 /* Sum the buffer into the required forces */
 112                 if (!bPBC || (!bScrew && fshift == NULL))
 113                 {
 114                     for(i=0; i<spas->nsend; i++)
 115                     {
 116                         rvec_inc(f[spas->a[i]],*vbuf);
 117                         vbuf++;
 118                     }
 119                 }
 120                 else
 121                 {
 122                     clear_ivec(vis);
 123                     vis[dim] = (dir==0 ? 1 : -1);
 124                     is = IVEC2IS(vis);
 125                     if (!bScrew)
 126                     {
 127                         /* Sum and add to shift forces */
 128                         for(i=0; i<spas->nsend; i++)
 129                         {
 130                             rvec_inc(f[spas->a[i]],*vbuf);
 131                             rvec_inc(fshift[is],*vbuf);
 132                             vbuf++;
 133                         }
 134                     }
 135                     else
 136                     {
 137                         /* Rotate the forces */
 138                         for(i=0; i<spas->nsend; i++)
 139                         {
 140                             f[spas->a[i]][XX] += (*vbuf)[XX];
 141                             f[spas->a[i]][YY] -= (*vbuf)[YY];
 142                             f[spas->a[i]][ZZ] -= (*vbuf)[ZZ];
 143                             if (fshift)
 144                             {
 145                                 rvec_inc(fshift[is],*vbuf);
 146                             }
 147                             vbuf++;
 148                         }
 149                     }
 150                 }
 151             }
 152         }
 153         else
 154         {
 155             /* Two cells, so we only need to communicate one way */
 156             spas = &spac->spas[d][0];
 157             n -= spas->nrecv;
 158             /* Send and receive the coordinates */
 159             dd_sendrecv_rvec(dd,d,dddirForward,
 160                              f+n,spas->nrecv,spac->vbuf,spas->nsend);
 161             /* Sum the buffer into the required forces */
 162             if (dd->bScrewPBC && dim == XX &&
 163                 (dd->ci[dim] == 0 ||
 164                  dd->ci[dim] == dd->nc[dim]-1))
 165             {
 166                 for(i=0; i<spas->nsend; i++)
 167                 {
 168                     /* Rotate the force */
 169                     f[spas->a[i]][XX] += spac->vbuf[i][XX];
 170                     f[spas->a[i]][YY] -= spac->vbuf[i][YY];
 171                     f[spas->a[i]][ZZ] -= spac->vbuf[i][ZZ];
 172                 }
 173             }
 174             else
 175             {
 176                 for(i=0; i<spas->nsend; i++)
 177                 {
 178                     rvec_inc(f[spas->a[i]],spac->vbuf[i]);
 179                 }
 180             }
 181         }
 182     }
 183 }
 184
 185 void dd_move_f_vsites(gmx_domdec_t *dd,rvec *f,rvec *fshift)
 186 {
 187     if (dd->vsite_comm)
 188     {
 189         dd_move_f_specat(dd,dd->vsite_comm,f,fshift);
 190     }
 191 }
 192
 193 void dd_clear_f_vsites(gmx_domdec_t *dd,rvec *f)
 194 {
 195     int i;
 196
 197     if (dd->vsite_comm)
 198     {
 199         for(i=dd->vsite_comm->at_start; i<dd->vsite_comm->at_end; i++)
 200         {
 201             clear_rvec(f[i]);
 202         }
 203     }
 204 }
 205
 206 static void dd_move_x_specat(gmx_domdec_t *dd,gmx_domdec_specat_comm_t *spac,
 207                              matrix box,rvec *x0,rvec *x1)
 208 {
 209     gmx_specatsend_t *spas;
 210     rvec *x,*vbuf,*rbuf;
 211     int  nvec,v,n,nn,ns0,ns1,nr0,nr1,nr,d,dim,dir,i;
 212     gmx_bool bPBC,bScrew=FALSE;
 213     rvec shift={0,0,0};
 214
 215     nvec = 1;
 216     if (x1)
 217     {
 218         nvec++;
 219     }
 220
 221     n = spac->at_start;
 222     for(d=0; d<dd->ndim; d++)
 223     {
 224         dim = dd->dim[d];
 225         if (dd->nc[dim] > 2)
 226         {
 227             /* Pulse the grid forward and backward */
 228             vbuf = spac->vbuf;
 229             for(dir=0; dir<2; dir++)
 230             {
 231                 if (dir == 0 && dd->ci[dim] == 0)
 232                 {
 233                     bPBC   = TRUE;
 234                     bScrew = (dd->bScrewPBC && dim == XX);
 235                     copy_rvec(box[dim],shift);
 236                 }
 237                 else if (dir == 1 && dd->ci[dim] == dd->nc[dim]-1)
 238                 {
 239                     bPBC = TRUE;
 240                     bScrew = (dd->bScrewPBC && dim == XX);
 241                     for(i=0; i<DIM; i++)
 242                     {
 243                         shift[i] = -box[dim][i];
 244                     }
 245                 }
 246                 else
 247                 {
 248                     bPBC = FALSE;
 249                     bScrew = FALSE;
 250                 }
 251                 spas = &spac->spas[d][dir];
 252                 for(v=0; v<nvec; v++)
 253                 {
 254                     x = (v == 0 ? x0 : x1);
 255                     /* Copy the required coordinates to the send buffer */
 256                     if (!bPBC)
 257                     {
 258                         /* Only copy */
 259                         for(i=0; i<spas->nsend; i++)
 260                         {
 261                             copy_rvec(x[spas->a[i]],*vbuf);
 262                             vbuf++;
 263                         }
 264                     }
 265                     else if (!bScrew)
 266                     {
 267                         /* Shift coordinates */
 268                         for(i=0; i<spas->nsend; i++)
 269                         {
 270                             rvec_add(x[spas->a[i]],shift,*vbuf);
 271                             vbuf++;
 272                         }
 273                     }
 274                     else
 275                     {
 276                         /* Shift and rotate coordinates */
 277                         for(i=0; i<spas->nsend; i++)
 278                         {
 279                             (*vbuf)[XX] =               x[spas->a[i]][XX] + shift[XX];
 280                             (*vbuf)[YY] = box[YY][YY] - x[spas->a[i]][YY] + shift[YY];
 281                             (*vbuf)[ZZ] = box[ZZ][ZZ] - x[spas->a[i]][ZZ] + shift[ZZ];
 282                             vbuf++;
 283                         }
 284                     }
 285                 }
 286             }
 287             /* Send and receive the coordinates */
 288             spas = spac->spas[d];
 289             ns0 = spas[0].nsend;
 290             nr0 = spas[0].nrecv;
 291             ns1 = spas[1].nsend;
 292             nr1 = spas[1].nrecv;
 293             if (nvec == 1)
 294             {
 295                 dd_sendrecv2_rvec(dd,d,
 296                                   spac->vbuf+ns0,ns1,x0+n    ,nr1,
 297                                   spac->vbuf    ,ns0,x0+n+nr1,nr0);
 298             }
 299             else
 300             {
 301                 /* Communicate both vectors in one buffer */
 302                 rbuf = spac->vbuf2;
 303                 dd_sendrecv2_rvec(dd,d,
 304                                   spac->vbuf+2*ns0,2*ns1,rbuf      ,2*nr1,
 305                                   spac->vbuf      ,2*ns0,rbuf+2*nr1,2*nr0);
 306                 /* Split the buffer into the two vectors */
 307                 nn = n;
 308                 for(dir=1; dir>=0; dir--)
 309                 {
 310                     nr = spas[dir].nrecv;
 311                     for(v=0; v<2; v++)
 312                     {
 313                         x = (v == 0 ? x0 : x1);
 314                         for(i=0; i<nr; i++)
 315                         {
 316                             copy_rvec(*rbuf,x[nn+i]);
 317                             rbuf++;
 318                         }
 319                     }
 320                     nn += nr;
 321                 }
 322             }
 323             n += nr0 + nr1;
 324         }
 325         else
 326         {
 327             spas = &spac->spas[d][0];
 328             /* Copy the required coordinates to the send buffer */
 329             vbuf = spac->vbuf;
 330             for(v=0; v<nvec; v++)
 331             {
 332                 x = (v == 0 ? x0 : x1);
 333                 if (dd->bScrewPBC && dim == XX &&
 334                     (dd->ci[XX] == 0 || dd->ci[XX] == dd->nc[XX]-1))
 335                 {
 336                     /* Here we only perform the rotation, the rest of the pbc
 337                      * is handled in the constraint or viste routines.
 338                      */
 339                     for(i=0; i<spas->nsend; i++)
 340                     {
 341                         (*vbuf)[XX] =               x[spas->a[i]][XX];
 342                         (*vbuf)[YY] = box[YY][YY] - x[spas->a[i]][YY];
 343                         (*vbuf)[ZZ] = box[ZZ][ZZ] - x[spas->a[i]][ZZ];
 344                         vbuf++;
 345                     }
 346                 }
 347                 else
 348                 {
 349                     for(i=0; i<spas->nsend; i++)
 350                     {
 351                         copy_rvec(x[spas->a[i]],*vbuf);
 352                         vbuf++;
 353                     }
 354                 }
 355             }
 356             /* Send and receive the coordinates */
 357             if (nvec == 1)
 358             {
 359                 dd_sendrecv_rvec(dd,d,dddirBackward,
 360                                  spac->vbuf,spas->nsend,x0+n,spas->nrecv);
 361             }
 362             else
 363             {
 364                 /* Communicate both vectors in one buffer */
 365                 rbuf = spac->vbuf2;
 366                 dd_sendrecv_rvec(dd,d,dddirBackward,
 367                                  spac->vbuf,2*spas->nsend,rbuf,2*spas->nrecv);
 368                 /* Split the buffer into the two vectors */
 369                 nr = spas[0].nrecv;
 370                 for(v=0; v<2; v++)
 371                 {
 372                     x = (v == 0 ? x0 : x1);
 373                     for(i=0; i<nr; i++)
 374                     {
 375                         copy_rvec(*rbuf,x[n+i]);
 376                         rbuf++;
 377                     }
 378                 }
 379             }
 380             n += spas->nrecv;
 381         }
 382     }
 383 }
 384
 385 void dd_move_x_constraints(gmx_domdec_t *dd,matrix box,rvec *x0,rvec *x1)
 386 {
 387     if (dd->constraint_comm)
 388     {
 389         dd_move_x_specat(dd,dd->constraint_comm,box,x0,x1);
 390     }
 391 }
 392
 393 void dd_move_x_vsites(gmx_domdec_t *dd,matrix box,rvec *x)
 394 {
 395     if (dd->vsite_comm)
 396     {
 397         dd_move_x_specat(dd,dd->vsite_comm,box,x,NULL);
 398     }
 399 }
 400
 401 int *dd_constraints_nlocalatoms(gmx_domdec_t *dd)
 402 {
 403     if (dd->constraints)
 404     {
 405         return dd->constraints->con_nlocat;
 406     }
 407     else
 408     {
 409         return NULL;
 410     }
 411 }
 412
 413 void dd_clear_local_constraint_indices(gmx_domdec_t *dd)
 414 {
 415     gmx_domdec_constraints_t *dc;
 416     int i;
 417
 418     dc = dd->constraints;
 419
 420     for(i=0; i<dc->ncon; i++)
 421     {
 422         dc->gc_req[dc->con_gl[i]] = 0;
 423     }
 424
 425     if (dd->constraint_comm)
 426     {
 427         for(i=dd->constraint_comm->at_start; i<dd->constraint_comm->at_end; i++)
 428         {
 429             dc->ga2la[dd->gatindex[i]] = -1;
 430         }
 431     }
 432 }
 433
 434 void dd_clear_local_vsite_indices(gmx_domdec_t *dd)
 435 {
 436     int i;
 437
 438     if (dd->vsite_comm)
 439     {
 440         for(i=dd->vsite_comm->at_start; i<dd->vsite_comm->at_end; i++)
 441         {
 442             dd->ga2la_vsite[dd->gatindex[i]] = -1;
 443         }
 444     }
 445 }
 446
 447 static int setup_specat_communication(gmx_domdec_t *dd,
 448                                       gmx_domdec_specat_comm_t *spac,
 449                                       int *ga2la_specat,
 450                                       int at_start,
 451                                       int vbuf_fac,
 452                                       const char *specat_type,
 453                                       const char *add_err)
 454 {
 455     int  nsend[2],nlast,nsend_zero[2]={0,0},*nsend_ptr;
 456     int  d,dim,ndir,dir,nr,ns,i,nrecv_local,n0,start,ireq,ind,buf[2];
 457     int  nat_tot_specat,nat_tot_prev,nalloc_old;
 458     gmx_bool bPBC,bFirst;
 459     gmx_specatsend_t *spas;
 460
 461     if (debug)
 462     {
 463         fprintf(debug,"Begin setup_specat_communication for %s\n",specat_type);
 464     }
 465
 466     /* nsend[0]: the number of atoms requested by this node only,
 467      *           we communicate this for more efficients checks
 468      * nsend[1]: the total number of requested atoms
 469      */
 470     nsend[0] = spac->nind_req;
 471     nsend[1] = nsend[0];
 472     nlast    = nsend[1];
 473     for(d=dd->ndim-1; d>=0; d--)
 474     {
 475         /* Pulse the grid forward and backward */
 476         dim = dd->dim[d];
 477         bPBC = (dim < dd->npbcdim);
 478         if (dd->nc[dim] == 2)
 479         {
 480             /* Only 2 cells, so we only need to communicate once */
 481             ndir = 1;
 482         }
 483         else
 484         {
 485             ndir = 2;
 486         }
 487         for(dir=0; dir<ndir; dir++)
 488         {
 489             if (!bPBC &&
 490                 dd->nc[dim] > 2 &&
 491                 ((dir == 0 && dd->ci[dim] == dd->nc[dim] - 1) ||
 492                  (dir == 1 && dd->ci[dim] == 0)))
 493             {
 494                 /* No pbc: the fist/last cell should not request atoms */
 495                 nsend_ptr = nsend_zero;
 496             }
 497             else
 498             {
 499                 nsend_ptr = nsend;
 500             }
 501             /* Communicate the number of indices */
 502             dd_sendrecv_int(dd,d,dir==0 ? dddirForward : dddirBackward,
 503                             nsend_ptr,2,spac->nreq[d][dir],2);
 504             nr = spac->nreq[d][dir][1];
 505             if (nlast+nr > spac->ind_req_nalloc)
 506             {
 507                 spac->ind_req_nalloc = over_alloc_dd(nlast+nr);
 508                 srenew(spac->ind_req,spac->ind_req_nalloc);
 509             }
 510             /* Communicate the indices */
 511             dd_sendrecv_int(dd,d,dir==0 ? dddirForward : dddirBackward,
 512                             spac->ind_req,nsend_ptr[1],spac->ind_req+nlast,nr);
 513             nlast += nr;
 514         }
 515         nsend[1] = nlast;
 516     }
 517     if (debug)
 518     {
 519         fprintf(debug,"Communicated the counts\n");
 520     }
 521
 522     /* Search for the requested atoms and communicate the indices we have */
 523     nat_tot_specat = at_start;
 524     nrecv_local = 0;
 525     for(d=0; d<dd->ndim; d++)
 526     {
 527         bFirst = (d == 0);
 528         /* Pulse the grid forward and backward */
 529         if (dd->dim[d] >= dd->npbcdim || dd->nc[dd->dim[d]] > 2)
 530         {
 531             ndir = 2;
 532         }
 533         else
 534         {
 535             ndir = 1;
 536         }
 537         nat_tot_prev = nat_tot_specat;
 538         for(dir=ndir-1; dir>=0; dir--)
 539         {
 540             if (nat_tot_specat > spac->bSendAtom_nalloc)
 541             {
 542                 nalloc_old = spac->bSendAtom_nalloc;
 543                 spac->bSendAtom_nalloc = over_alloc_dd(nat_tot_specat);
 544                 srenew(spac->bSendAtom,spac->bSendAtom_nalloc);
 545                 for(i=nalloc_old; i<spac->bSendAtom_nalloc; i++)
 546                 {
 547                     spac->bSendAtom[i] = FALSE;
 548                 }
 549             }
 550             spas = &spac->spas[d][dir];
 551             n0 = spac->nreq[d][dir][0];
 552             nr = spac->nreq[d][dir][1];
 553             if (debug)
 554             {
 555                 fprintf(debug,"dim=%d, dir=%d, searching for %d atoms\n",
 556                         d,dir,nr);
 557             }
 558             start = nlast - nr;
 559             spas->nsend = 0;
 560             nsend[0] = 0;
 561             for(i=0; i<nr; i++)
 562             {
 563                 ireq = spac->ind_req[start+i];
 564                 ind = -1;
 565                 /* Check if this is a home atom and if so ind will be set */
 566                 if (!ga2la_get_home(dd->ga2la,ireq,&ind))
 567                 {
 568                     /* Search in the communicated atoms */
 569                     ind = ga2la_specat[ireq];
 570                 }
 571                 if (ind >= 0)
 572                 {
 573                     if (i < n0 || !spac->bSendAtom[ind])
 574                     {
 575                         if (spas->nsend+1 > spas->a_nalloc)
 576                         {
 577                             spas->a_nalloc = over_alloc_large(spas->nsend+1);
 578                             srenew(spas->a,spas->a_nalloc);
 579                         }
 580                         /* Store the local index so we know which coordinates
 581                          * to send out later.
 582                          */
 583                         spas->a[spas->nsend] = ind;
 584                         spac->bSendAtom[ind] = TRUE;
 585                         if (spas->nsend+1 > spac->ibuf_nalloc)
 586                         {
 587                             spac->ibuf_nalloc = over_alloc_large(spas->nsend+1);
 588                             srenew(spac->ibuf,spac->ibuf_nalloc);
 589                         }
 590                         /* Store the global index so we can send it now */
 591                         spac->ibuf[spas->nsend] = ireq;
 592                         if (i < n0)
 593                         {
 594                             nsend[0]++;
 595                         }
 596                         spas->nsend++;
 597                     }
 598                 }
 599             }
 600             nlast = start;
 601             /* Clear the local flags */
 602             for(i=0; i<spas->nsend; i++)
 603             {
 604                 spac->bSendAtom[spas->a[i]] = FALSE;
 605             }
 606             /* Send and receive the number of indices to communicate */
 607             nsend[1] = spas->nsend;
 608             dd_sendrecv_int(dd,d,dir==0 ? dddirBackward : dddirForward,
 609                             nsend,2,buf,2);
 610             if (debug)
 611             {
 612                 fprintf(debug,"Send to node %d, %d (%d) indices, "
 613                         "receive from node %d, %d (%d) indices\n",
 614                         dd->neighbor[d][1-dir],nsend[1],nsend[0],
 615                         dd->neighbor[d][dir],buf[1],buf[0]);
 616                 if (gmx_debug_at)
 617                 {
 618                     for(i=0; i<spas->nsend; i++)
 619                     {
 620                         fprintf(debug," %d",spac->ibuf[i]+1);
 621                     }
 622                     fprintf(debug,"\n");
 623                 }
 624             }
 625             nrecv_local += buf[0];
 626             spas->nrecv  = buf[1];
 627             if (nat_tot_specat + spas->nrecv > dd->gatindex_nalloc)
 628             {
 629                 dd->gatindex_nalloc =
 630                     over_alloc_dd(nat_tot_specat + spas->nrecv);
 631                 srenew(dd->gatindex,dd->gatindex_nalloc);
 632             }
 633             /* Send and receive the indices */
 634             dd_sendrecv_int(dd,d,dir==0 ? dddirBackward : dddirForward,
 635                             spac->ibuf,spas->nsend,
 636                             dd->gatindex+nat_tot_specat,spas->nrecv);
 637             nat_tot_specat += spas->nrecv;
 638         }
 639
 640         /* Allocate the x/f communication buffers */
 641         ns = spac->spas[d][0].nsend;
 642         nr = spac->spas[d][0].nrecv;
 643         if (ndir == 2)
 644         {
 645             ns += spac->spas[d][1].nsend;
 646             nr += spac->spas[d][1].nrecv;
 647         }
 648         if (vbuf_fac*ns > spac->vbuf_nalloc)
 649         {
 650             spac->vbuf_nalloc = over_alloc_dd(vbuf_fac*ns);
 651             srenew(spac->vbuf,spac->vbuf_nalloc);
 652         }
 653         if (vbuf_fac == 2 && vbuf_fac*nr > spac->vbuf2_nalloc)
 654         {
 655             spac->vbuf2_nalloc = over_alloc_dd(vbuf_fac*nr);
 656             srenew(spac->vbuf2,spac->vbuf2_nalloc);
 657         }
 658
 659         /* Make a global to local index for the communication atoms */
 660         for(i=nat_tot_prev; i<nat_tot_specat; i++)
 661         {
 662             ga2la_specat[dd->gatindex[i]] = i;
 663         }
 664     }
 665
 666     /* Check that in the end we got the number of atoms we asked for */
 667     if (nrecv_local != spac->nind_req)
 668     {
 669         if (debug)
 670         {
 671             fprintf(debug,"Requested %d, received %d (tot recv %d)\n",
 672                     spac->nind_req,nrecv_local,nat_tot_specat-at_start);
 673             if (gmx_debug_at)
 674             {
 675                 for(i=0; i<spac->nind_req; i++)
 676                 {
 677                     fprintf(debug," %s%d",
 678                             ga2la_specat[spac->ind_req[i]]>=0 ? "" : "!",
 679                             spac->ind_req[i]+1);
 680                 }
 681                 fprintf(debug,"\n");
 682             }
 683         }
 684         fprintf(stderr,"\nDD cell %d %d %d: Neighboring cells do not have atoms:",
 685                 dd->ci[XX],dd->ci[YY],dd->ci[ZZ]);
 686         for(i=0; i<spac->nind_req; i++)
 687         {
 688             if (ga2la_specat[spac->ind_req[i]] < 0)
 689             {
 690                 fprintf(stderr," %d",spac->ind_req[i]+1);
 691             }
 692         }
 693         fprintf(stderr,"\n");
 694         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.",
 695                   dd->ci[XX],dd->ci[YY],dd->ci[ZZ],
 696                   nrecv_local,spac->nind_req,specat_type,
 697                   specat_type,add_err,
 698                   dd->bGridJump ? " or use the -rcon option of mdrun" : "");
 699     }
 700
 701     spac->at_start = at_start;
 702     spac->at_end   = nat_tot_specat;
 703
 704     if (debug)
 705     {
 706         fprintf(debug,"Done setup_specat_communication\n");
 707     }
 708
 709     return nat_tot_specat;
 710 }
 711
 712 static void walk_out(int con,int con_offset,int a,int offset,int nrec,
 713                      int ncon1,const t_iatom *ia1,const t_iatom *ia2,
 714                      const t_blocka *at2con,
 715                      const gmx_ga2la_t ga2la,gmx_bool bHomeConnect,
 716                      gmx_domdec_constraints_t *dc,
 717                      gmx_domdec_specat_comm_t *dcc,
 718                      t_ilist *il_local)
 719 {
 720     int a1_gl,a2_gl,a_loc,i,coni,b;
 721     const t_iatom *iap;
 722
 723     if (dc->gc_req[con_offset+con] == 0)
 724     {
 725         /* Add this non-home constraint to the list */
 726         if (dc->ncon+1 > dc->con_nalloc)
 727         {
 728             dc->con_nalloc = over_alloc_large(dc->ncon+1);
 729             srenew(dc->con_gl,dc->con_nalloc);
 730             srenew(dc->con_nlocat,dc->con_nalloc);
 731         }
 732         dc->con_gl[dc->ncon] = con_offset + con;
 733         dc->con_nlocat[dc->ncon] = (bHomeConnect ? 1 : 0);
 734         dc->gc_req[con_offset+con] = 1;
 735         if (il_local->nr + 3 > il_local->nalloc)
 736         {
 737             il_local->nalloc = over_alloc_dd(il_local->nr+3);
 738             srenew(il_local->iatoms,il_local->nalloc);
 739         }
 740         iap = constr_iatomptr(ncon1,ia1,ia2,con);
 741         il_local->iatoms[il_local->nr++] = iap[0];
 742         a1_gl = offset + iap[1];
 743         a2_gl = offset + iap[2];
 744         /* The following indexing code can probably be optizimed */
 745         if (ga2la_get_home(ga2la,a1_gl,&a_loc))
 746         {
 747             il_local->iatoms[il_local->nr++] = a_loc;
 748         }
 749         else
 750         {
 751             /* We set this index later */
 752             il_local->iatoms[il_local->nr++] = -a1_gl - 1;
 753         }
 754         if (ga2la_get_home(ga2la,a2_gl,&a_loc))
 755         {
 756             il_local->iatoms[il_local->nr++] = a_loc;
 757         }
 758         else
 759         {
 760             /* We set this index later */
 761             il_local->iatoms[il_local->nr++] = -a2_gl - 1;
 762         }
 763         dc->ncon++;
 764     }
 765     /* Check to not ask for the same atom more than once */
 766     if (dc->ga2la[offset+a] == -1)
 767     {
 768         assert(dcc);
 769         /* Add this non-home atom to the list */
 770         if (dcc->nind_req+1 > dcc->ind_req_nalloc)
 771         {
 772             dcc->ind_req_nalloc = over_alloc_large(dcc->nind_req+1);
 773             srenew(dcc->ind_req,dcc->ind_req_nalloc);
 774         }
 775         dcc->ind_req[dcc->nind_req++] = offset + a;
 776         /* Temporarily mark with -2, we get the index later */
 777         dc->ga2la[offset+a] = -2;
 778     }
 779
 780     if (nrec > 0)
 781     {
 782         for(i=at2con->index[a]; i<at2con->index[a+1]; i++)
 783         {
 784             coni = at2con->a[i];
 785             if (coni != con)
 786             {
 787                 /* Walk further */
 788                 iap = constr_iatomptr(ncon1,ia1,ia2,coni);
 789                 if (a == iap[1])
 790                 {
 791                     b = iap[2];
 792                 }
 793                 else
 794                 {
 795                     b = iap[1];
 796                 }
 797                 if (!ga2la_get_home(ga2la,offset+b,&a_loc))
 798                 {
 799                     walk_out(coni,con_offset,b,offset,nrec-1,
 800                              ncon1,ia1,ia2,at2con,
 801                              ga2la,FALSE,dc,dcc,il_local);
 802                 }
 803             }
 804         }
 805     }
 806 }
 807
 808 int dd_make_local_constraints(gmx_domdec_t *dd,int at_start,
 809                               gmx_mtop_t *mtop,
 810                               gmx_constr_t constr,int nrec,
 811                               t_ilist *il_local)
 812 {
 813     t_blocka *at2con_mt,*at2con;
 814     gmx_ga2la_t ga2la;
 815     int ncon1,ncon2;
 816     gmx_molblock_t *molb;
 817     t_iatom *ia1,*ia2,*iap;
 818     int nhome,a,a_gl,a_mol,a_loc,b_lo,offset,mb,molnr,b_mol,i,con,con_offset;
 819     gmx_domdec_constraints_t *dc;
 820     int at_end,*ga2la_specat,j;
 821
 822     dc = dd->constraints;
 823
 824     at2con_mt = atom2constraints_moltype(constr);
 825     ga2la  = dd->ga2la;
 826
 827     dc->ncon     = 0;
 828     il_local->nr = 0;
 829     nhome = 0;
 830     if (dd->constraint_comm)
 831     {
 832         dd->constraint_comm->nind_req = 0;
 833     }
 834     for(a=0; a<dd->nat_home; a++)
 835     {
 836         a_gl = dd->gatindex[a];
 837
 838         gmx_mtop_atomnr_to_molblock_ind(mtop,a_gl,&mb,&molnr,&a_mol);
 839         molb = &mtop->molblock[mb];
 840
 841         ncon1 = mtop->moltype[molb->type].ilist[F_CONSTR].nr/3;
 842         ncon2 = mtop->moltype[molb->type].ilist[F_CONSTRNC].nr/3;
 843         if (ncon1 > 0 || ncon2 > 0)
 844         {
 845             ia1 = mtop->moltype[molb->type].ilist[F_CONSTR].iatoms;
 846             ia2 = mtop->moltype[molb->type].ilist[F_CONSTRNC].iatoms;
 847
 848             /* Calculate the global constraint number offset for the molecule.
 849              * This is only required for the global index to make sure
 850              * that we use each constraint only once.
 851              */
 852             con_offset = dc->molb_con_offset[mb] + molnr*dc->molb_ncon_mol[mb];
 853
 854             /* The global atom number offset for this molecule */
 855             offset = a_gl - a_mol;
 856             at2con = &at2con_mt[molb->type];
 857             for(i=at2con->index[a_mol]; i<at2con->index[a_mol+1]; i++)
 858             {
 859                 con = at2con->a[i];
 860                 iap = constr_iatomptr(ncon1,ia1,ia2,con);
 861                 if (a_mol == iap[1])
 862                 {
 863                     b_mol = iap[2];
 864                 }
 865                 else
 866                 {
 867                     b_mol = iap[1];
 868                 }
 869                 if (ga2la_get_home(ga2la,offset+b_mol,&a_loc))
 870                 {
 871                     /* Add this fully home constraint at the first atom */
 872                     if (a_mol < b_mol)
 873                     {
 874                         if (dc->ncon+1 > dc->con_nalloc)
 875                         {
 876                             dc->con_nalloc = over_alloc_large(dc->ncon+1);
 877                             srenew(dc->con_gl,dc->con_nalloc);
 878                             srenew(dc->con_nlocat,dc->con_nalloc);
 879                         }
 880                         dc->con_gl[dc->ncon] = con_offset + con;
 881                         dc->con_nlocat[dc->ncon] = 2;
 882                         if (il_local->nr + 3 > il_local->nalloc)
 883                         {
 884                             il_local->nalloc = over_alloc_dd(il_local->nr + 3);
 885                             srenew(il_local->iatoms,il_local->nalloc);
 886                         }
 887                         b_lo = a_loc;
 888                         il_local->iatoms[il_local->nr++] = iap[0];
 889                         il_local->iatoms[il_local->nr++] = (a_gl == iap[1] ? a    : b_lo);
 890                         il_local->iatoms[il_local->nr++] = (a_gl == iap[1] ? b_lo : a   );
 891                         dc->ncon++;
 892                         nhome++;
 893                     }
 894                 }
 895                 else
 896                 {
 897                     /* We need the nrec constraints coupled to this constraint,
 898                      * so we need to walk out of the home cell by nrec+1 atoms,
 899                      * since already atom bg is not locally present.
 900                      * Therefore we call walk_out with nrec recursions to go
 901                      * after this first call.
 902                      */
 903                     walk_out(con,con_offset,b_mol,offset,nrec,
 904                              ncon1,ia1,ia2,at2con,
 905                              dd->ga2la,TRUE,dc,dd->constraint_comm,il_local);
 906                 }
 907             }
 908         }
 909     }
 910
 911     if (debug)
 912     {
 913         fprintf(debug,
 914                 "Constraints: home %3d border %3d atoms: %3d\n",
 915                 nhome,dc->ncon-nhome,
 916                 dd->constraint_comm ? dd->constraint_comm->nind_req : 0);
 917     }
 918
 919     if (dd->constraint_comm) {
 920         at_end =
 921             setup_specat_communication(dd,dd->constraint_comm,
 922                                        dd->constraints->ga2la,
 923                                        at_start,2,
 924                                        "constraint"," or lincs-order");
 925
 926         /* Fill in the missing indices */
 927         ga2la_specat = dd->constraints->ga2la;
 928         for(i=0; i<il_local->nr; i+=3)
 929         {
 930             iap = il_local->iatoms + i;
 931             for(j=1; j<3; j++)
 932             {
 933                 if (iap[j] < 0)
 934                 {
 935                     iap[j] = ga2la_specat[-iap[j]-1];
 936                 }
 937             }
 938         }
 939     }
 940     else
 941     {
 942         at_end = at_start;
 943     }
 944
 945     return at_end;
 946 }
 947
 948 int dd_make_local_vsites(gmx_domdec_t *dd,int at_start,t_ilist *lil)
 949 {
 950     gmx_domdec_specat_comm_t *spac;
 951     int  *ga2la_specat;
 952     int  ftype,nral,i,j,gat,a;
 953     t_ilist *lilf;
 954     t_iatom *iatoms;
 955     int  at_end;
 956
 957     spac         = dd->vsite_comm;
 958     ga2la_specat = dd->ga2la_vsite;
 959
 960     spac->nind_req = 0;
 961     /* Loop over all the home vsites */
 962     for(ftype=0; ftype<F_NRE; ftype++)
 963     {
 964         if (interaction_function[ftype].flags & IF_VSITE)
 965         {
 966             nral = NRAL(ftype);
 967             lilf = &lil[ftype];
 968             for(i=0; i<lilf->nr; i+=1+nral)
 969             {
 970                 iatoms = lilf->iatoms + i;
 971                 /* Check if we have the other atoms */
 972                 for(j=1; j<1+nral; j++)
 973                 {
 974                     if (iatoms[j] < 0) {
 975                         /* This is not a home atom,
 976                          * we need to ask our neighbors.
 977                          */
 978                         a = -iatoms[j] - 1;
 979                         /* Check to not ask for the same atom more than once */
 980                         if (ga2la_specat[a] == -1)
 981                         {
 982                             /* Add this non-home atom to the list */
 983                             if (spac->nind_req+1 > spac->ind_req_nalloc)
 984                             {
 985                                 spac->ind_req_nalloc =
 986                                     over_alloc_small(spac->nind_req+1);
 987                                 srenew(spac->ind_req,spac->ind_req_nalloc);
 988                             }
 989                             spac->ind_req[spac->nind_req++] = a;
 990                             /* Temporarily mark with -2,
 991                              * we get the index later.
 992                              */
 993                             ga2la_specat[a] = -2;
 994                         }
 995                     }
 996                 }
 997             }
 998         }
 999     }
1000
1001     at_end = setup_specat_communication(dd,dd->vsite_comm,ga2la_specat,
1002                                         at_start,1,"vsite","");
1003
1004     /* Fill in the missing indices */
1005     for(ftype=0; ftype<F_NRE; ftype++)
1006     {
1007         if (interaction_function[ftype].flags & IF_VSITE)
1008         {
1009             nral = NRAL(ftype);
1010             lilf = &lil[ftype];
1011             for(i=0; i<lilf->nr; i+=1+nral)
1012             {
1013                 iatoms = lilf->iatoms + i;
1014                 for(j=1; j<1+nral; j++)
1015                 {
1016                     if (iatoms[j] < 0)
1017                     {
1018                         iatoms[j] = ga2la_specat[-iatoms[j]-1];
1019                     }
1020                 }
1021             }
1022         }
1023     }
1024
1025     return at_end;
1026 }
1027
1028 void init_domdec_constraints(gmx_domdec_t *dd,
1029                              int natoms,gmx_mtop_t *mtop,
1030                              gmx_constr_t constr)
1031 {
1032     gmx_domdec_constraints_t *dc;
1033     gmx_molblock_t *molb;
1034     int mb,ncon,c,a;
1035
1036     if (debug)
1037     {
1038         fprintf(debug,"Begin init_domdec_constraints\n");
1039     }
1040
1041     snew(dd->constraints,1);
1042     dc = dd->constraints;
1043
1044     snew(dc->molb_con_offset,mtop->nmolblock);
1045     snew(dc->molb_ncon_mol,mtop->nmolblock);
1046
1047     ncon = 0;
1048     for(mb=0; mb<mtop->nmolblock; mb++)
1049     {
1050         molb = &mtop->molblock[mb];
1051         dc->molb_con_offset[mb] = ncon;
1052         dc->molb_ncon_mol[mb] =
1053             mtop->moltype[molb->type].ilist[F_CONSTR].nr/3 +
1054             mtop->moltype[molb->type].ilist[F_CONSTRNC].nr/3;
1055         ncon += molb->nmol*dc->molb_ncon_mol[mb];
1056     }
1057
1058     snew(dc->gc_req,ncon);
1059     for(c=0; c<ncon; c++)
1060     {
1061         dc->gc_req[c] = 0;
1062     }
1063
1064     snew(dc->ga2la,natoms);
1065     for(a=0; a<natoms; a++)
1066     {
1067         dc->ga2la[a] = -1;
1068     }
1069
1070     snew(dd->constraint_comm,1);
1071 }
1072
1073 void init_domdec_vsites(gmx_domdec_t *dd,int natoms)
1074 {
1075     int i;
1076     gmx_domdec_constraints_t *dc;
1077
1078     if (debug)
1079     {
1080         fprintf(debug,"Begin init_domdec_vsites\n");
1081     }
1082
1083     snew(dd->ga2la_vsite,natoms);
1084     for(i=0; i<natoms; i++)
1085     {
1086         dd->ga2la_vsite[i] = -1;
1087     }
1088
1089     snew(dd->vsite_comm,1);
1090 }