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