src/gromacs/gmxlib/pbc.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #ifdef HAVE_CONFIG_H
  38 #include <config.h>
  39 #endif
  40
  41 #include <math.h>
  42 #include "sysstuff.h"
  43 #include "typedefs.h"
  44 #include "vec.h"
  45 #include "maths.h"
  46 #include "main.h"
  47 #include "pbc.h"
  48 #include "smalloc.h"
  49 #include "txtdump.h"
  50 #include "gmx_fatal.h"
  51 #include "names.h"
  52 #include "macros.h"
  53 #include "gmx_omp_nthreads.h"
  54
  55 /* Skip 0 so we have more chance of detecting if we forgot to call set_pbc. */
  56 enum {
  57     epbcdxRECTANGULAR = 1, epbcdxTRICLINIC,
  58     epbcdx2D_RECT,       epbcdx2D_TRIC,
  59     epbcdx1D_RECT,       epbcdx1D_TRIC,
  60     epbcdxSCREW_RECT,    epbcdxSCREW_TRIC,
  61     epbcdxNOPBC,         epbcdxUNSUPPORTED
  62 };
  63
  64 /* Margin factor for error message and correction if the box is too skewed */
  65 #define BOX_MARGIN         1.0010
  66 #define BOX_MARGIN_CORRECT 1.0005
  67
  68 int ePBC2npbcdim(int ePBC)
  69 {
  70     int npbcdim = 0;
  71
  72     switch (ePBC)
  73     {
  74         case epbcXYZ:   npbcdim = 3; break;
  75         case epbcXY:    npbcdim = 2; break;
  76         case epbcSCREW: npbcdim = 3; break;
  77         case epbcNONE:  npbcdim = 0; break;
  78         default: gmx_fatal(FARGS, "Unknown ePBC=%d in ePBC2npbcdim", ePBC);
  79     }
  80
  81     return npbcdim;
  82 }
  83
  84 int inputrec2nboundeddim(t_inputrec *ir)
  85 {
  86     if (ir->nwall == 2 && ir->ePBC == epbcXY)
  87     {
  88         return 3;
  89     }
  90     else
  91     {
  92         return ePBC2npbcdim(ir->ePBC);
  93     }
  94 }
  95
  96 void dump_pbc(FILE *fp, t_pbc *pbc)
  97 {
  98     rvec sum_box;
  99
 100     fprintf(fp, "ePBCDX = %d\n", pbc->ePBCDX);
 101     pr_rvecs(fp, 0, "box", pbc->box, DIM);
 102     pr_rvecs(fp, 0, "fbox_diag", &pbc->fbox_diag, 1);
 103     pr_rvecs(fp, 0, "hbox_diag", &pbc->hbox_diag, 1);
 104     pr_rvecs(fp, 0, "mhbox_diag", &pbc->mhbox_diag, 1);
 105     rvec_add(pbc->hbox_diag, pbc->mhbox_diag, sum_box);
 106     pr_rvecs(fp, 0, "sum of the above two", &sum_box, 1);
 107     fprintf(fp, "max_cutoff2 = %g\n", pbc->max_cutoff2);
 108     fprintf(fp, "bLimitDistance = %s\n", EBOOL(pbc->bLimitDistance));
 109     fprintf(fp, "limit_distance2 = %g\n", pbc->limit_distance2);
 110     fprintf(fp, "ntric_vec = %d\n", pbc->ntric_vec);
 111     if (pbc->ntric_vec > 0)
 112     {
 113         pr_ivecs(fp, 0, "tric_shift", pbc->tric_shift, pbc->ntric_vec, FALSE);
 114         pr_rvecs(fp, 0, "tric_vec", pbc->tric_vec, pbc->ntric_vec);
 115     }
 116 }
 117
 118 const char *check_box(int ePBC, matrix box)
 119 {
 120     const char *ptr;
 121
 122     if (ePBC == -1)
 123     {
 124         ePBC = guess_ePBC(box);
 125     }
 126
 127     if (ePBC == epbcNONE)
 128     {
 129         return NULL;
 130     }
 131
 132     if ((box[XX][YY] != 0) || (box[XX][ZZ] != 0) || (box[YY][ZZ] != 0))
 133     {
 134         ptr = "Only triclinic boxes with the first vector parallel to the x-axis and the second vector in the xy-plane are supported.";
 135     }
 136     else if (ePBC == epbcSCREW && (box[YY][XX] != 0 || box[ZZ][XX] != 0))
 137     {
 138         ptr = "The unit cell can not have off-diagonal x-components with screw pbc";
 139     }
 140     else if (fabs(box[YY][XX]) > BOX_MARGIN*0.5*box[XX][XX] ||
 141              (ePBC != epbcXY &&
 142               (fabs(box[ZZ][XX]) > BOX_MARGIN*0.5*box[XX][XX] ||
 143                fabs(box[ZZ][YY]) > BOX_MARGIN*0.5*box[YY][YY])))
 144     {
 145         ptr = "Triclinic box is too skewed.";
 146     }
 147     else
 148     {
 149         ptr = NULL;
 150     }
 151
 152     return ptr;
 153 }
 154
 155 real max_cutoff2(int ePBC, matrix box)
 156 {
 157     real min_hv2, min_ss;
 158
 159     /* Physical limitation of the cut-off
 160      * by half the length of the shortest box vector.
 161      */
 162     min_hv2 = min(0.25*norm2(box[XX]), 0.25*norm2(box[YY]));
 163     if (ePBC != epbcXY)
 164     {
 165         min_hv2 = min(min_hv2, 0.25*norm2(box[ZZ]));
 166     }
 167
 168     /* Limitation to the smallest diagonal element due to optimizations:
 169      * checking only linear combinations of single box-vectors (2 in x)
 170      * in the grid search and pbc_dx is a lot faster
 171      * than checking all possible combinations.
 172      */
 173     if (ePBC == epbcXY)
 174     {
 175         min_ss = min(box[XX][XX], box[YY][YY]);
 176     }
 177     else
 178     {
 179         min_ss = min(box[XX][XX], min(box[YY][YY]-fabs(box[ZZ][YY]), box[ZZ][ZZ]));
 180     }
 181
 182     return min(min_hv2, min_ss*min_ss);
 183 }
 184
 185 /* this one is mostly harmless... */
 186 static gmx_bool bWarnedGuess = FALSE;
 187
 188 int guess_ePBC(matrix box)
 189 {
 190     int ePBC;
 191
 192     if (box[XX][XX] > 0 && box[YY][YY] > 0 && box[ZZ][ZZ] > 0)
 193     {
 194         ePBC = epbcXYZ;
 195     }
 196     else if (box[XX][XX] > 0 && box[YY][YY] > 0 && box[ZZ][ZZ] == 0)
 197     {
 198         ePBC = epbcXY;
 199     }
 200     else if (box[XX][XX] == 0 && box[YY][YY] == 0 && box[ZZ][ZZ] == 0)
 201     {
 202         ePBC = epbcNONE;
 203     }
 204     else
 205     {
 206         if (!bWarnedGuess)
 207         {
 208             fprintf(stderr, "WARNING: Unsupported box diagonal %f %f %f, "
 209                     "will not use periodic boundary conditions\n\n",
 210                     box[XX][XX], box[YY][YY], box[ZZ][ZZ]);
 211             bWarnedGuess = TRUE;
 212         }
 213         ePBC = epbcNONE;
 214     }
 215
 216     if (debug)
 217     {
 218         fprintf(debug, "Guessed pbc = %s from the box matrix\n", epbc_names[ePBC]);
 219     }
 220
 221     return ePBC;
 222 }
 223
 224 static int correct_box_elem(FILE *fplog, int step, tensor box, int v, int d)
 225 {
 226     int shift, maxshift = 10;
 227
 228     shift = 0;
 229
 230     /* correct elem d of vector v with vector d */
 231     while (box[v][d] > BOX_MARGIN_CORRECT*0.5*box[d][d])
 232     {
 233         if (fplog)
 234         {
 235             fprintf(fplog, "Step %d: correcting invalid box:\n", step);
 236             pr_rvecs(fplog, 0, "old box", box, DIM);
 237         }
 238         rvec_dec(box[v], box[d]);
 239         shift--;
 240         if (fplog)
 241         {
 242             pr_rvecs(fplog, 0, "new box", box, DIM);
 243         }
 244         if (shift <= -maxshift)
 245         {
 246             gmx_fatal(FARGS,
 247                       "Box was shifted at least %d times. Please see log-file.",
 248                       maxshift);
 249         }
 250     }
 251     while (box[v][d] < -BOX_MARGIN_CORRECT*0.5*box[d][d])
 252     {
 253         if (fplog)
 254         {
 255             fprintf(fplog, "Step %d: correcting invalid box:\n", step);
 256             pr_rvecs(fplog, 0, "old box", box, DIM);
 257         }
 258         rvec_inc(box[v], box[d]);
 259         shift++;
 260         if (fplog)
 261         {
 262             pr_rvecs(fplog, 0, "new box", box, DIM);
 263         }
 264         if (shift >= maxshift)
 265         {
 266             gmx_fatal(FARGS,
 267                       "Box was shifted at least %d times. Please see log-file.",
 268                       maxshift);
 269         }
 270     }
 271
 272     return shift;
 273 }
 274
 275 gmx_bool correct_box(FILE *fplog, int step, tensor box, t_graph *graph)
 276 {
 277     int      zy, zx, yx, i;
 278     gmx_bool bCorrected;
 279
 280     /* check if the box still obeys the restrictions, if not, correct it */
 281     zy = correct_box_elem(fplog, step, box, ZZ, YY);
 282     zx = correct_box_elem(fplog, step, box, ZZ, XX);
 283     yx = correct_box_elem(fplog, step, box, YY, XX);
 284
 285     bCorrected = (zy || zx || yx);
 286
 287     if (bCorrected && graph)
 288     {
 289         /* correct the graph */
 290         for (i = graph->at_start; i < graph->at_end; i++)
 291         {
 292             graph->ishift[i][YY] -= graph->ishift[i][ZZ]*zy;
 293             graph->ishift[i][XX] -= graph->ishift[i][ZZ]*zx;
 294             graph->ishift[i][XX] -= graph->ishift[i][YY]*yx;
 295         }
 296     }
 297
 298     return bCorrected;
 299 }
 300
 301 int ndof_com(t_inputrec *ir)
 302 {
 303     int n = 0;
 304
 305     switch (ir->ePBC)
 306     {
 307         case epbcXYZ:
 308         case epbcNONE:
 309             n = 3;
 310             break;
 311         case epbcXY:
 312             n = (ir->nwall == 0 ? 3 : 2);
 313             break;
 314         case epbcSCREW:
 315             n = 1;
 316             break;
 317         default:
 318             gmx_incons("Unknown pbc in calc_nrdf");
 319     }
 320
 321     return n;
 322 }
 323
 324 static void low_set_pbc(t_pbc *pbc, int ePBC, ivec *dd_nc, matrix box)
 325 {
 326     int         order[5] = {0, -1, 1, -2, 2};
 327     int         ii, jj, kk, i, j, k, d, dd, jc, kc, npbcdim, shift;
 328     ivec        bPBC;
 329     real        d2old, d2new, d2new_c;
 330     rvec        trial, pos;
 331     gmx_bool    bXY, bUse;
 332     const char *ptr;
 333
 334     pbc->ePBC      = ePBC;
 335     pbc->ndim_ePBC = ePBC2npbcdim(ePBC);
 336
 337     copy_mat(box, pbc->box);
 338     pbc->bLimitDistance = FALSE;
 339     pbc->max_cutoff2    = 0;
 340     pbc->dim            = -1;
 341
 342     for (i = 0; (i < DIM); i++)
 343     {
 344         pbc->fbox_diag[i]  =  box[i][i];
 345         pbc->hbox_diag[i]  =  pbc->fbox_diag[i]*0.5;
 346         pbc->mhbox_diag[i] = -pbc->hbox_diag[i];
 347     }
 348
 349     ptr = check_box(ePBC, box);
 350     if (ePBC == epbcNONE)
 351     {
 352         pbc->ePBCDX = epbcdxNOPBC;
 353     }
 354     else if (ptr)
 355     {
 356         fprintf(stderr,   "Warning: %s\n", ptr);
 357         pr_rvecs(stderr, 0, "         Box", box, DIM);
 358         fprintf(stderr,   "         Can not fix pbc.\n");
 359         pbc->ePBCDX          = epbcdxUNSUPPORTED;
 360         pbc->bLimitDistance  = TRUE;
 361         pbc->limit_distance2 = 0;
 362     }
 363     else
 364     {
 365         if (ePBC == epbcSCREW && dd_nc)
 366         {
 367             /* This combinated should never appear here */
 368             gmx_incons("low_set_pbc called with screw pbc and dd_nc != NULL");
 369         }
 370
 371         npbcdim = 0;
 372         for (i = 0; i < DIM; i++)
 373         {
 374             if ((dd_nc && (*dd_nc)[i] > 1) || (ePBC == epbcXY && i == ZZ))
 375             {
 376                 bPBC[i] = 0;
 377             }
 378             else
 379             {
 380                 bPBC[i] = 1;
 381                 npbcdim++;
 382             }
 383         }
 384         switch (npbcdim)
 385         {
 386             case 1:
 387                 /* 1D pbc is not an mdp option and it is therefore only used
 388                  * with single shifts.
 389                  */
 390                 pbc->ePBCDX = epbcdx1D_RECT;
 391                 for (i = 0; i < DIM; i++)
 392                 {
 393                     if (bPBC[i])
 394                     {
 395                         pbc->dim = i;
 396                     }
 397                 }
 398                 for (i = 0; i < pbc->dim; i++)
 399                 {
 400                     if (pbc->box[pbc->dim][i] != 0)
 401                     {
 402                         pbc->ePBCDX = epbcdx1D_TRIC;
 403                     }
 404                 }
 405                 break;
 406             case 2:
 407                 pbc->ePBCDX = epbcdx2D_RECT;
 408                 for (i = 0; i < DIM; i++)
 409                 {
 410                     if (!bPBC[i])
 411                     {
 412                         pbc->dim = i;
 413                     }
 414                 }
 415                 for (i = 0; i < DIM; i++)
 416                 {
 417                     if (bPBC[i])
 418                     {
 419                         for (j = 0; j < i; j++)
 420                         {
 421                             if (pbc->box[i][j] != 0)
 422                             {
 423                                 pbc->ePBCDX = epbcdx2D_TRIC;
 424                             }
 425                         }
 426                     }
 427                 }
 428                 break;
 429             case 3:
 430                 if (ePBC != epbcSCREW)
 431                 {
 432                     if (TRICLINIC(box))
 433                     {
 434                         pbc->ePBCDX = epbcdxTRICLINIC;
 435                     }
 436                     else
 437                     {
 438                         pbc->ePBCDX = epbcdxRECTANGULAR;
 439                     }
 440                 }
 441                 else
 442                 {
 443                     pbc->ePBCDX = (box[ZZ][YY] == 0 ? epbcdxSCREW_RECT : epbcdxSCREW_TRIC);
 444                     if (pbc->ePBCDX == epbcdxSCREW_TRIC)
 445                     {
 446                         fprintf(stderr,
 447                                 "Screw pbc is not yet implemented for triclinic boxes.\n"
 448                                 "Can not fix pbc.\n");
 449                         pbc->ePBCDX = epbcdxUNSUPPORTED;
 450                     }
 451                 }
 452                 break;
 453             default:
 454                 gmx_fatal(FARGS, "Incorrect number of pbc dimensions with DD: %d",
 455                           npbcdim);
 456         }
 457         pbc->max_cutoff2 = max_cutoff2(ePBC, box);
 458
 459         if (pbc->ePBCDX == epbcdxTRICLINIC ||
 460             pbc->ePBCDX == epbcdx2D_TRIC ||
 461             pbc->ePBCDX == epbcdxSCREW_TRIC)
 462         {
 463             if (debug)
 464             {
 465                 pr_rvecs(debug, 0, "Box", box, DIM);
 466                 fprintf(debug, "max cutoff %.3f\n", sqrt(pbc->max_cutoff2));
 467             }
 468             pbc->ntric_vec = 0;
 469             /* We will only use single shifts, but we will check a few
 470              * more shifts to see if there is a limiting distance
 471              * above which we can not be sure of the correct distance.
 472              */
 473             for (kk = 0; kk < 5; kk++)
 474             {
 475                 k = order[kk];
 476                 if (!bPBC[ZZ] && k != 0)
 477                 {
 478                     continue;
 479                 }
 480                 for (jj = 0; jj < 5; jj++)
 481                 {
 482                     j = order[jj];
 483                     if (!bPBC[YY] && j != 0)
 484                     {
 485                         continue;
 486                     }
 487                     for (ii = 0; ii < 3; ii++)
 488                     {
 489                         i = order[ii];
 490                         if (!bPBC[XX] && i != 0)
 491                         {
 492                             continue;
 493                         }
 494                         /* A shift is only useful when it is trilinic */
 495                         if (j != 0 || k != 0)
 496                         {
 497                             d2old = 0;
 498                             d2new = 0;
 499                             for (d = 0; d < DIM; d++)
 500                             {
 501                                 trial[d] = i*box[XX][d] + j*box[YY][d] + k*box[ZZ][d];
 502                                 /* Choose the vector within the brick around 0,0,0 that
 503                                  * will become the shortest due to shift try.
 504                                  */
 505                                 if (d == pbc->dim)
 506                                 {
 507                                     trial[d] = 0;
 508                                     pos[d]   = 0;
 509                                 }
 510                                 else
 511                                 {
 512                                     if (trial[d] < 0)
 513                                     {
 514                                         pos[d] = min( pbc->hbox_diag[d], -trial[d]);
 515                                     }
 516                                     else
 517                                     {
 518                                         pos[d] = max(-pbc->hbox_diag[d], -trial[d]);
 519                                     }
 520                                 }
 521                                 d2old += sqr(pos[d]);
 522                                 d2new += sqr(pos[d] + trial[d]);
 523                             }
 524                             if (BOX_MARGIN*d2new < d2old)
 525                             {
 526                                 if (j < -1 || j > 1 || k < -1 || k > 1)
 527                                 {
 528                                     /* Check if there is a single shift vector
 529                                      * that decreases this distance even more.
 530                                      */
 531                                     jc = 0;
 532                                     kc = 0;
 533                                     if (j < -1 || j > 1)
 534                                     {
 535                                         jc = j/2;
 536                                     }
 537                                     if (k < -1 || k > 1)
 538                                     {
 539                                         kc = k/2;
 540                                     }
 541                                     d2new_c = 0;
 542                                     for (d = 0; d < DIM; d++)
 543                                     {
 544                                         d2new_c += sqr(pos[d] + trial[d]
 545                                                        - jc*box[YY][d] - kc*box[ZZ][d]);
 546                                     }
 547                                     if (d2new_c > BOX_MARGIN*d2new)
 548                                     {
 549                                         /* Reject this shift vector, as there is no a priori limit
 550                                          * to the number of shifts that decrease distances.
 551                                          */
 552                                         if (!pbc->bLimitDistance || d2new <  pbc->limit_distance2)
 553                                         {
 554                                             pbc->limit_distance2 = d2new;
 555                                         }
 556                                         pbc->bLimitDistance = TRUE;
 557                                     }
 558                                 }
 559                                 else
 560                                 {
 561                                     /* Check if shifts with one box vector less do better */
 562                                     bUse = TRUE;
 563                                     for (dd = 0; dd < DIM; dd++)
 564                                     {
 565                                         shift = (dd == 0 ? i : (dd == 1 ? j : k));
 566                                         if (shift)
 567                                         {
 568                                             d2new_c = 0;
 569                                             for (d = 0; d < DIM; d++)
 570                                             {
 571                                                 d2new_c += sqr(pos[d] + trial[d] - shift*box[dd][d]);
 572                                             }
 573                                             if (d2new_c <= BOX_MARGIN*d2new)
 574                                             {
 575                                                 bUse = FALSE;
 576                                             }
 577                                         }
 578                                     }
 579                                     if (bUse)
 580                                     {
 581                                         /* Accept this shift vector. */
 582                                         if (pbc->ntric_vec >= MAX_NTRICVEC)
 583                                         {
 584                                             fprintf(stderr, "\nWARNING: Found more than %d triclinic correction vectors, ignoring some.\n"
 585                                                     "  There is probably something wrong with your box.\n", MAX_NTRICVEC);
 586                                             pr_rvecs(stderr, 0, "         Box", box, DIM);
 587                                         }
 588                                         else
 589                                         {
 590                                             copy_rvec(trial, pbc->tric_vec[pbc->ntric_vec]);
 591                                             pbc->tric_shift[pbc->ntric_vec][XX] = i;
 592                                             pbc->tric_shift[pbc->ntric_vec][YY] = j;
 593                                             pbc->tric_shift[pbc->ntric_vec][ZZ] = k;
 594                                             pbc->ntric_vec++;
 595                                         }
 596                                     }
 597                                 }
 598                                 if (debug)
 599                                 {
 600                                     fprintf(debug, "  tricvec %2d = %2d %2d %2d  %5.2f %5.2f  %5.2f %5.2f %5.2f  %5.2f %5.2f %5.2f\n",
 601                                             pbc->ntric_vec, i, j, k,
 602                                             sqrt(d2old), sqrt(d2new),
 603                                             trial[XX], trial[YY], trial[ZZ],
 604                                             pos[XX], pos[YY], pos[ZZ]);
 605                                 }
 606                             }
 607                         }
 608                     }
 609                 }
 610             }
 611         }
 612     }
 613 }
 614
 615 void set_pbc(t_pbc *pbc, int ePBC, matrix box)
 616 {
 617     if (ePBC == -1)
 618     {
 619         ePBC = guess_ePBC(box);
 620     }
 621
 622     low_set_pbc(pbc, ePBC, NULL, box);
 623 }
 624
 625 t_pbc *set_pbc_dd(t_pbc *pbc, int ePBC,
 626                   gmx_domdec_t *dd, gmx_bool bSingleDir, matrix box)
 627 {
 628     ivec nc2;
 629     int  npbcdim, i;
 630
 631     if (dd == NULL)
 632     {
 633         npbcdim = DIM;
 634     }
 635     else
 636     {
 637         if (ePBC == epbcSCREW && dd->nc[XX] > 1)
 638         {
 639             /* The rotation has been taken care of during coordinate communication */
 640             ePBC = epbcXYZ;
 641         }
 642         npbcdim = 0;
 643         for (i = 0; i < DIM; i++)
 644         {
 645             if (dd->nc[i] <= (bSingleDir ? 1 : 2))
 646             {
 647                 nc2[i] = 1;
 648                 if (!(ePBC == epbcXY && i == ZZ))
 649                 {
 650                     npbcdim++;
 651                 }
 652             }
 653             else
 654             {
 655                 nc2[i] = dd->nc[i];
 656             }
 657         }
 658     }
 659
 660     if (npbcdim > 0)
 661     {
 662         low_set_pbc(pbc, ePBC, npbcdim < DIM ? &nc2 : NULL, box);
 663     }
 664
 665     return (npbcdim > 0 ? pbc : NULL);
 666 }
 667
 668 void pbc_dx(const t_pbc *pbc, const rvec x1, const rvec x2, rvec dx)
 669 {
 670     int      i, j;
 671     rvec     dx_start, trial;
 672     real     d2min, d2trial;
 673     gmx_bool bRot;
 674
 675     rvec_sub(x1, x2, dx);
 676
 677     switch (pbc->ePBCDX)
 678     {
 679         case epbcdxRECTANGULAR:
 680             for (i = 0; i < DIM; i++)
 681             {
 682                 while (dx[i] > pbc->hbox_diag[i])
 683                 {
 684                     dx[i] -= pbc->fbox_diag[i];
 685                 }
 686                 while (dx[i] <= pbc->mhbox_diag[i])
 687                 {
 688                     dx[i] += pbc->fbox_diag[i];
 689                 }
 690             }
 691             break;
 692         case epbcdxTRICLINIC:
 693             for (i = DIM-1; i >= 0; i--)
 694             {
 695                 while (dx[i] > pbc->hbox_diag[i])
 696                 {
 697                     for (j = i; j >= 0; j--)
 698                     {
 699                         dx[j] -= pbc->box[i][j];
 700                     }
 701                 }
 702                 while (dx[i] <= pbc->mhbox_diag[i])
 703                 {
 704                     for (j = i; j >= 0; j--)
 705                     {
 706                         dx[j] += pbc->box[i][j];
 707                     }
 708                 }
 709             }
 710             /* dx is the distance in a rectangular box */
 711             d2min = norm2(dx);
 712             if (d2min > pbc->max_cutoff2)
 713             {
 714                 copy_rvec(dx, dx_start);
 715                 d2min = norm2(dx);
 716                 /* Now try all possible shifts, when the distance is within max_cutoff
 717                  * it must be the shortest possible distance.
 718                  */
 719                 i = 0;
 720                 while ((d2min > pbc->max_cutoff2) && (i < pbc->ntric_vec))
 721                 {
 722                     rvec_add(dx_start, pbc->tric_vec[i], trial);
 723                     d2trial = norm2(trial);
 724                     if (d2trial < d2min)
 725                     {
 726                         copy_rvec(trial, dx);
 727                         d2min = d2trial;
 728                     }
 729                     i++;
 730                 }
 731             }
 732             break;
 733         case epbcdx2D_RECT:
 734             for (i = 0; i < DIM; i++)
 735             {
 736                 if (i != pbc->dim)
 737                 {
 738                     while (dx[i] > pbc->hbox_diag[i])
 739                     {
 740                         dx[i] -= pbc->fbox_diag[i];
 741                     }
 742                     while (dx[i] <= pbc->mhbox_diag[i])
 743                     {
 744                         dx[i] += pbc->fbox_diag[i];
 745                     }
 746                 }
 747             }
 748             break;
 749         case epbcdx2D_TRIC:
 750             d2min = 0;
 751             for (i = DIM-1; i >= 0; i--)
 752             {
 753                 if (i != pbc->dim)
 754                 {
 755                     while (dx[i] > pbc->hbox_diag[i])
 756                     {
 757                         for (j = i; j >= 0; j--)
 758                         {
 759                             dx[j] -= pbc->box[i][j];
 760                         }
 761                     }
 762                     while (dx[i] <= pbc->mhbox_diag[i])
 763                     {
 764                         for (j = i; j >= 0; j--)
 765                         {
 766                             dx[j] += pbc->box[i][j];
 767                         }
 768                     }
 769                     d2min += dx[i]*dx[i];
 770                 }
 771             }
 772             if (d2min > pbc->max_cutoff2)
 773             {
 774                 copy_rvec(dx, dx_start);
 775                 d2min = norm2(dx);
 776                 /* Now try all possible shifts, when the distance is within max_cutoff
 777                  * it must be the shortest possible distance.
 778                  */
 779                 i = 0;
 780                 while ((d2min > pbc->max_cutoff2) && (i < pbc->ntric_vec))
 781                 {
 782                     rvec_add(dx_start, pbc->tric_vec[i], trial);
 783                     d2trial = 0;
 784                     for (j = 0; j < DIM; j++)
 785                     {
 786                         if (j != pbc->dim)
 787                         {
 788                             d2trial += trial[j]*trial[j];
 789                         }
 790                     }
 791                     if (d2trial < d2min)
 792                     {
 793                         copy_rvec(trial, dx);
 794                         d2min = d2trial;
 795                     }
 796                     i++;
 797                 }
 798             }
 799             break;
 800         case epbcdxSCREW_RECT:
 801             /* The shift definition requires x first */
 802             bRot = FALSE;
 803             while (dx[XX] > pbc->hbox_diag[XX])
 804             {
 805                 dx[XX] -= pbc->fbox_diag[XX];
 806                 bRot    = !bRot;
 807             }
 808             while (dx[XX] <= pbc->mhbox_diag[XX])
 809             {
 810                 dx[XX] += pbc->fbox_diag[YY];
 811                 bRot    = !bRot;
 812             }
 813             if (bRot)
 814             {
 815                 /* Rotate around the x-axis in the middle of the box */
 816                 dx[YY] = pbc->box[YY][YY] - x1[YY] - x2[YY];
 817                 dx[ZZ] = pbc->box[ZZ][ZZ] - x1[ZZ] - x2[ZZ];
 818             }
 819             /* Normal pbc for y and z */
 820             for (i = YY; i <= ZZ; i++)
 821             {
 822                 while (dx[i] > pbc->hbox_diag[i])
 823                 {
 824                     dx[i] -= pbc->fbox_diag[i];
 825                 }
 826                 while (dx[i] <= pbc->mhbox_diag[i])
 827                 {
 828                     dx[i] += pbc->fbox_diag[i];
 829                 }
 830             }
 831             break;
 832         case epbcdxNOPBC:
 833         case epbcdxUNSUPPORTED:
 834             break;
 835         default:
 836             gmx_fatal(FARGS, "Internal error in pbc_dx, set_pbc has not been called");
 837             break;
 838     }
 839 }
 840
 841 int pbc_dx_aiuc(const t_pbc *pbc, const rvec x1, const rvec x2, rvec dx)
 842 {
 843     int  i, j, is;
 844     rvec dx_start, trial;
 845     real d2min, d2trial;
 846     ivec ishift, ishift_start;
 847
 848     rvec_sub(x1, x2, dx);
 849     clear_ivec(ishift);
 850
 851     switch (pbc->ePBCDX)
 852     {
 853         case epbcdxRECTANGULAR:
 854             for (i = 0; i < DIM; i++)
 855             {
 856                 if (dx[i] > pbc->hbox_diag[i])
 857                 {
 858                     dx[i] -=  pbc->fbox_diag[i];
 859                     ishift[i]--;
 860                 }
 861                 else if (dx[i] <= pbc->mhbox_diag[i])
 862                 {
 863                     dx[i] +=  pbc->fbox_diag[i];
 864                     ishift[i]++;
 865                 }
 866             }
 867             break;
 868         case epbcdxTRICLINIC:
 869             /* For triclinic boxes the performance difference between
 870              * if/else and two while loops is negligible.
 871              * However, the while version can cause extreme delays
 872              * before a simulation crashes due to large forces which
 873              * can cause unlimited displacements.
 874              * Also allowing multiple shifts would index fshift beyond bounds.
 875              */
 876             for (i = DIM-1; i >= 1; i--)
 877             {
 878                 if (dx[i] > pbc->hbox_diag[i])
 879                 {
 880                     for (j = i; j >= 0; j--)
 881                     {
 882                         dx[j] -= pbc->box[i][j];
 883                     }
 884                     ishift[i]--;
 885                 }
 886                 else if (dx[i] <= pbc->mhbox_diag[i])
 887                 {
 888                     for (j = i; j >= 0; j--)
 889                     {
 890                         dx[j] += pbc->box[i][j];
 891                     }
 892                     ishift[i]++;
 893                 }
 894             }
 895             /* Allow 2 shifts in x */
 896             if (dx[XX] > pbc->hbox_diag[XX])
 897             {
 898                 dx[XX] -= pbc->fbox_diag[XX];
 899                 ishift[XX]--;
 900                 if (dx[XX] > pbc->hbox_diag[XX])
 901                 {
 902                     dx[XX] -= pbc->fbox_diag[XX];
 903                     ishift[XX]--;
 904                 }
 905             }
 906             else if (dx[XX] <= pbc->mhbox_diag[XX])
 907             {
 908                 dx[XX] += pbc->fbox_diag[XX];
 909                 ishift[XX]++;
 910                 if (dx[XX] <= pbc->mhbox_diag[XX])
 911                 {
 912                     dx[XX] += pbc->fbox_diag[XX];
 913                     ishift[XX]++;
 914                 }
 915             }
 916             /* dx is the distance in a rectangular box */
 917             d2min = norm2(dx);
 918             if (d2min > pbc->max_cutoff2)
 919             {
 920                 copy_rvec(dx, dx_start);
 921                 copy_ivec(ishift, ishift_start);
 922                 d2min = norm2(dx);
 923                 /* Now try all possible shifts, when the distance is within max_cutoff
 924                  * it must be the shortest possible distance.
 925                  */
 926                 i = 0;
 927                 while ((d2min > pbc->max_cutoff2) && (i < pbc->ntric_vec))
 928                 {
 929                     rvec_add(dx_start, pbc->tric_vec[i], trial);
 930                     d2trial = norm2(trial);
 931                     if (d2trial < d2min)
 932                     {
 933                         copy_rvec(trial, dx);
 934                         ivec_add(ishift_start, pbc->tric_shift[i], ishift);
 935                         d2min = d2trial;
 936                     }
 937                     i++;
 938                 }
 939             }
 940             break;
 941         case epbcdx2D_RECT:
 942             for (i = 0; i < DIM; i++)
 943             {
 944                 if (i != pbc->dim)
 945                 {
 946                     if (dx[i] > pbc->hbox_diag[i])
 947                     {
 948                         dx[i] -= pbc->fbox_diag[i];
 949                         ishift[i]--;
 950                     }
 951                     else if (dx[i] <= pbc->mhbox_diag[i])
 952                     {
 953                         dx[i] += pbc->fbox_diag[i];
 954                         ishift[i]++;
 955                     }
 956                 }
 957             }
 958             break;
 959         case epbcdx2D_TRIC:
 960             d2min = 0;
 961             for (i = DIM-1; i >= 1; i--)
 962             {
 963                 if (i != pbc->dim)
 964                 {
 965                     if (dx[i] > pbc->hbox_diag[i])
 966                     {
 967                         for (j = i; j >= 0; j--)
 968                         {
 969                             dx[j] -= pbc->box[i][j];
 970                         }
 971                         ishift[i]--;
 972                     }
 973                     else if (dx[i] <= pbc->mhbox_diag[i])
 974                     {
 975                         for (j = i; j >= 0; j--)
 976                         {
 977                             dx[j] += pbc->box[i][j];
 978                         }
 979                         ishift[i]++;
 980                     }
 981                     d2min += dx[i]*dx[i];
 982                 }
 983             }
 984             if (pbc->dim != XX)
 985             {
 986                 /* Allow 2 shifts in x */
 987                 if (dx[XX] > pbc->hbox_diag[XX])
 988                 {
 989                     dx[XX] -= pbc->fbox_diag[XX];
 990                     ishift[XX]--;
 991                     if (dx[XX] > pbc->hbox_diag[XX])
 992                     {
 993                         dx[XX] -= pbc->fbox_diag[XX];
 994                         ishift[XX]--;
 995                     }
 996                 }
 997                 else if (dx[XX] <= pbc->mhbox_diag[XX])
 998                 {
 999                     dx[XX] += pbc->fbox_diag[XX];
1000                     ishift[XX]++;
1001                     if (dx[XX] <= pbc->mhbox_diag[XX])
1002                     {
1003                         dx[XX] += pbc->fbox_diag[XX];
1004                         ishift[XX]++;
1005                     }
1006                 }
1007                 d2min += dx[XX]*dx[XX];
1008             }
1009             if (d2min > pbc->max_cutoff2)
1010             {
1011                 copy_rvec(dx, dx_start);
1012                 copy_ivec(ishift, ishift_start);
1013                 /* Now try all possible shifts, when the distance is within max_cutoff
1014                  * it must be the shortest possible distance.
1015                  */
1016                 i = 0;
1017                 while ((d2min > pbc->max_cutoff2) && (i < pbc->ntric_vec))
1018                 {
1019                     rvec_add(dx_start, pbc->tric_vec[i], trial);
1020                     d2trial = 0;
1021                     for (j = 0; j < DIM; j++)
1022                     {
1023                         if (j != pbc->dim)
1024                         {
1025                             d2trial += trial[j]*trial[j];
1026                         }
1027                     }
1028                     if (d2trial < d2min)
1029                     {
1030                         copy_rvec(trial, dx);
1031                         ivec_add(ishift_start, pbc->tric_shift[i], ishift);
1032                         d2min = d2trial;
1033                     }
1034                     i++;
1035                 }
1036             }
1037             break;
1038         case epbcdx1D_RECT:
1039             i = pbc->dim;
1040             if (dx[i] > pbc->hbox_diag[i])
1041             {
1042                 dx[i] -= pbc->fbox_diag[i];
1043                 ishift[i]--;
1044             }
1045             else if (dx[i] <= pbc->mhbox_diag[i])
1046             {
1047                 dx[i] += pbc->fbox_diag[i];
1048                 ishift[i]++;
1049             }
1050             break;
1051         case epbcdx1D_TRIC:
1052             i = pbc->dim;
1053             if (dx[i] > pbc->hbox_diag[i])
1054             {
1055                 rvec_dec(dx, pbc->box[i]);
1056                 ishift[i]--;
1057             }
1058             else if (dx[i] <= pbc->mhbox_diag[i])
1059             {
1060                 rvec_inc(dx, pbc->box[i]);
1061                 ishift[i]++;
1062             }
1063             break;
1064         case epbcdxSCREW_RECT:
1065             /* The shift definition requires x first */
1066             if (dx[XX] > pbc->hbox_diag[XX])
1067             {
1068                 dx[XX] -= pbc->fbox_diag[XX];
1069                 ishift[XX]--;
1070             }
1071             else if (dx[XX] <= pbc->mhbox_diag[XX])
1072             {
1073                 dx[XX] += pbc->fbox_diag[XX];
1074                 ishift[XX]++;
1075             }
1076             if (ishift[XX] == 1 || ishift[XX] == -1)
1077             {
1078                 /* Rotate around the x-axis in the middle of the box */
1079                 dx[YY] = pbc->box[YY][YY] - x1[YY] - x2[YY];
1080                 dx[ZZ] = pbc->box[ZZ][ZZ] - x1[ZZ] - x2[ZZ];
1081             }
1082             /* Normal pbc for y and z */
1083             for (i = YY; i <= ZZ; i++)
1084             {
1085                 if (dx[i] > pbc->hbox_diag[i])
1086                 {
1087                     dx[i] -= pbc->fbox_diag[i];
1088                     ishift[i]--;
1089                 }
1090                 else if (dx[i] <= pbc->mhbox_diag[i])
1091                 {
1092                     dx[i] += pbc->fbox_diag[i];
1093                     ishift[i]++;
1094                 }
1095             }
1096             break;
1097         case epbcdxNOPBC:
1098         case epbcdxUNSUPPORTED:
1099             break;
1100         default:
1101             gmx_fatal(FARGS, "Internal error in pbc_dx_aiuc, set_pbc_dd or set_pbc has not been called");
1102             break;
1103     }
1104
1105     is = IVEC2IS(ishift);
1106     if (debug)
1107     {
1108         range_check_mesg(is, 0, SHIFTS, "PBC shift vector index range check.");
1109     }
1110
1111     return is;
1112 }
1113
1114 void pbc_dx_d(const t_pbc *pbc, const dvec x1, const dvec x2, dvec dx)
1115 {
1116     int      i, j;
1117     dvec     dx_start, trial;
1118     double   d2min, d2trial;
1119     gmx_bool bRot;
1120
1121     dvec_sub(x1, x2, dx);
1122
1123     switch (pbc->ePBCDX)
1124     {
1125         case epbcdxRECTANGULAR:
1126         case epbcdx2D_RECT:
1127             for (i = 0; i < DIM; i++)
1128             {
1129                 if (i != pbc->dim)
1130                 {
1131                     while (dx[i] > pbc->hbox_diag[i])
1132                     {
1133                         dx[i] -= pbc->fbox_diag[i];
1134                     }
1135                     while (dx[i] <= pbc->mhbox_diag[i])
1136                     {
1137                         dx[i] += pbc->fbox_diag[i];
1138                     }
1139                 }
1140             }
1141             break;
1142         case epbcdxTRICLINIC:
1143         case epbcdx2D_TRIC:
1144             d2min = 0;
1145             for (i = DIM-1; i >= 0; i--)
1146             {
1147                 if (i != pbc->dim)
1148                 {
1149                     while (dx[i] > pbc->hbox_diag[i])
1150                     {
1151                         for (j = i; j >= 0; j--)
1152                         {
1153                             dx[j] -= pbc->box[i][j];
1154                         }
1155                     }
1156                     while (dx[i] <= pbc->mhbox_diag[i])
1157                     {
1158                         for (j = i; j >= 0; j--)
1159                         {
1160                             dx[j] += pbc->box[i][j];
1161                         }
1162                     }
1163                     d2min += dx[i]*dx[i];
1164                 }
1165             }
1166             if (d2min > pbc->max_cutoff2)
1167             {
1168                 copy_dvec(dx, dx_start);
1169                 /* Now try all possible shifts, when the distance is within max_cutoff
1170                  * it must be the shortest possible distance.
1171                  */
1172                 i = 0;
1173                 while ((d2min > pbc->max_cutoff2) && (i < pbc->ntric_vec))
1174                 {
1175                     for (j = 0; j < DIM; j++)
1176                     {
1177                         trial[j] = dx_start[j] + pbc->tric_vec[i][j];
1178                     }
1179                     d2trial = 0;
1180                     for (j = 0; j < DIM; j++)
1181                     {
1182                         if (j != pbc->dim)
1183                         {
1184                             d2trial += trial[j]*trial[j];
1185                         }
1186                     }
1187                     if (d2trial < d2min)
1188                     {
1189                         copy_dvec(trial, dx);
1190                         d2min = d2trial;
1191                     }
1192                     i++;
1193                 }
1194             }
1195             break;
1196         case epbcdxSCREW_RECT:
1197             /* The shift definition requires x first */
1198             bRot = FALSE;
1199             while (dx[XX] > pbc->hbox_diag[XX])
1200             {
1201                 dx[XX] -= pbc->fbox_diag[XX];
1202                 bRot    = !bRot;
1203             }
1204             while (dx[XX] <= pbc->mhbox_diag[XX])
1205             {
1206                 dx[XX] += pbc->fbox_diag[YY];
1207                 bRot    = !bRot;
1208             }
1209             if (bRot)
1210             {
1211                 /* Rotate around the x-axis in the middle of the box */
1212                 dx[YY] = pbc->box[YY][YY] - x1[YY] - x2[YY];
1213                 dx[ZZ] = pbc->box[ZZ][ZZ] - x1[ZZ] - x2[ZZ];
1214             }
1215             /* Normal pbc for y and z */
1216             for (i = YY; i <= ZZ; i++)
1217             {
1218                 while (dx[i] > pbc->hbox_diag[i])
1219                 {
1220                     dx[i] -= pbc->fbox_diag[i];
1221                 }
1222                 while (dx[i] <= pbc->mhbox_diag[i])
1223                 {
1224                     dx[i] += pbc->fbox_diag[i];
1225                 }
1226             }
1227             break;
1228         case epbcdxNOPBC:
1229         case epbcdxUNSUPPORTED:
1230             break;
1231         default:
1232             gmx_fatal(FARGS, "Internal error in pbc_dx, set_pbc has not been called");
1233             break;
1234     }
1235 }
1236
1237 gmx_bool image_rect(ivec xi, ivec xj, ivec box_size, real rlong2, int *shift, real *r2)
1238 {
1239     int     m, t;
1240     int     dx, b, b_2;
1241     real    dxr, rij2;
1242
1243     rij2 = 0.0;
1244     t    = 0;
1245     for (m = 0; (m < DIM); m++)
1246     {
1247         dx  = xi[m]-xj[m];
1248         t  *= DIM;
1249         b   = box_size[m];
1250         b_2 = b/2;
1251         if (dx < -b_2)
1252         {
1253             t  += 2;
1254             dx += b;
1255         }
1256         else if (dx > b_2)
1257         {
1258             dx -= b;
1259         }
1260         else
1261         {
1262             t += 1;
1263         }
1264         dxr   = dx;
1265         rij2 += dxr*dxr;
1266         if (rij2 >= rlong2)
1267         {
1268             return FALSE;
1269         }
1270     }
1271
1272     *shift = t;
1273     *r2    = rij2;
1274     return TRUE;
1275 }
1276
1277 gmx_bool image_cylindric(ivec xi, ivec xj, ivec box_size, real rlong2,
1278                          int *shift, real *r2)
1279 {
1280     int     m, t;
1281     int     dx, b, b_2;
1282     real    dxr, rij2;
1283
1284     rij2 = 0.0;
1285     t    = 0;
1286     for (m = 0; (m < DIM); m++)
1287     {
1288         dx  = xi[m]-xj[m];
1289         t  *= DIM;
1290         b   = box_size[m];
1291         b_2 = b/2;
1292         if (dx < -b_2)
1293         {
1294             t  += 2;
1295             dx += b;
1296         }
1297         else if (dx > b_2)
1298         {
1299             dx -= b;
1300         }
1301         else
1302         {
1303             t += 1;
1304         }
1305
1306         dxr   = dx;
1307         rij2 += dxr*dxr;
1308         if (m < ZZ)
1309         {
1310             if (rij2 >= rlong2)
1311             {
1312                 return FALSE;
1313             }
1314         }
1315     }
1316
1317     *shift = t;
1318     *r2    = rij2;
1319     return TRUE;
1320 }
1321
1322 void calc_shifts(matrix box, rvec shift_vec[])
1323 {
1324     int k, l, m, d, n, test;
1325
1326     n = 0;
1327     for (m = -D_BOX_Z; m <= D_BOX_Z; m++)
1328     {
1329         for (l = -D_BOX_Y; l <= D_BOX_Y; l++)
1330         {
1331             for (k = -D_BOX_X; k <= D_BOX_X; k++, n++)
1332             {
1333                 test = XYZ2IS(k, l, m);
1334                 if (n != test)
1335                 {
1336                     gmx_incons("inconsistent shift numbering");
1337                 }
1338                 for (d = 0; d < DIM; d++)
1339                 {
1340                     shift_vec[n][d] = k*box[XX][d] + l*box[YY][d] + m*box[ZZ][d];
1341                 }
1342             }
1343         }
1344     }
1345 }
1346
1347 void calc_box_center(int ecenter, matrix box, rvec box_center)
1348 {
1349     int d, m;
1350
1351     clear_rvec(box_center);
1352     switch (ecenter)
1353     {
1354         case ecenterTRIC:
1355             for (m = 0; (m < DIM); m++)
1356             {
1357                 for (d = 0; d < DIM; d++)
1358                 {
1359                     box_center[d] += 0.5*box[m][d];
1360                 }
1361             }
1362             break;
1363         case ecenterRECT:
1364             for (d = 0; d < DIM; d++)
1365             {
1366                 box_center[d] = 0.5*box[d][d];
1367             }
1368             break;
1369         case ecenterZERO:
1370             break;
1371         default:
1372             gmx_fatal(FARGS, "Unsupported value %d for ecenter", ecenter);
1373     }
1374 }
1375
1376 void calc_triclinic_images(matrix box, rvec img[])
1377 {
1378     int i;
1379
1380     /* Calculate 3 adjacent images in the xy-plane */
1381     copy_rvec(box[0], img[0]);
1382     copy_rvec(box[1], img[1]);
1383     if (img[1][XX] < 0)
1384     {
1385         svmul(-1, img[1], img[1]);
1386     }
1387     rvec_sub(img[1], img[0], img[2]);
1388
1389     /* Get the next 3 in the xy-plane as mirror images */
1390     for (i = 0; i < 3; i++)
1391     {
1392         svmul(-1, img[i], img[3+i]);
1393     }
1394
1395     /* Calculate the first 4 out of xy-plane images */
1396     copy_rvec(box[2], img[6]);
1397     if (img[6][XX] < 0)
1398     {
1399         svmul(-1, img[6], img[6]);
1400     }
1401     for (i = 0; i < 3; i++)
1402     {
1403         rvec_add(img[6], img[i+1], img[7+i]);
1404     }
1405
1406     /* Mirror the last 4 from the previous in opposite rotation */
1407     for (i = 0; i < 4; i++)
1408     {
1409         svmul(-1, img[6 + (2+i) % 4], img[10+i]);
1410     }
1411 }
1412
1413 void calc_compact_unitcell_vertices(int ecenter, matrix box, rvec vert[])
1414 {
1415     rvec img[NTRICIMG], box_center;
1416     int  n, i, j, tmp[4], d;
1417
1418     calc_triclinic_images(box, img);
1419
1420     n = 0;
1421     for (i = 2; i <= 5; i += 3)
1422     {
1423         tmp[0] = i-1;
1424         if (i == 2)
1425         {
1426             tmp[1] = 8;
1427         }
1428         else
1429         {
1430             tmp[1] = 6;
1431         }
1432         tmp[2] = (i+1) % 6;
1433         tmp[3] = tmp[1]+4;
1434         for (j = 0; j < 4; j++)
1435         {
1436             for (d = 0; d < DIM; d++)
1437             {
1438                 vert[n][d] = img[i][d]+img[tmp[j]][d]+img[tmp[(j+1)%4]][d];
1439             }
1440             n++;
1441         }
1442     }
1443     for (i = 7; i <= 13; i += 6)
1444     {
1445         tmp[0] = (i-7)/2;
1446         tmp[1] = tmp[0]+1;
1447         if (i == 7)
1448         {
1449             tmp[2] = 8;
1450         }
1451         else
1452         {
1453             tmp[2] = 10;
1454         }
1455         tmp[3] = i-1;
1456         for (j = 0; j < 4; j++)
1457         {
1458             for (d = 0; d < DIM; d++)
1459             {
1460                 vert[n][d] = img[i][d]+img[tmp[j]][d]+img[tmp[(j+1)%4]][d];
1461             }
1462             n++;
1463         }
1464     }
1465     for (i = 9; i <= 11; i += 2)
1466     {
1467         if (i == 9)
1468         {
1469             tmp[0] = 3;
1470         }
1471         else
1472         {
1473             tmp[0] = 0;
1474         }
1475         tmp[1] = tmp[0]+1;
1476         if (i == 9)
1477         {
1478             tmp[2] = 6;
1479         }
1480         else
1481         {
1482             tmp[2] = 12;
1483         }
1484         tmp[3] = i-1;
1485         for (j = 0; j < 4; j++)
1486         {
1487             for (d = 0; d < DIM; d++)
1488             {
1489                 vert[n][d] = img[i][d]+img[tmp[j]][d]+img[tmp[(j+1)%4]][d];
1490             }
1491             n++;
1492         }
1493     }
1494
1495     calc_box_center(ecenter, box, box_center);
1496     for (i = 0; i < NCUCVERT; i++)
1497     {
1498         for (d = 0; d < DIM; d++)
1499         {
1500             vert[i][d] = vert[i][d]*0.25+box_center[d];
1501         }
1502     }
1503 }
1504
1505 int *compact_unitcell_edges()
1506 {
1507     /* this is an index in vert[] (see calc_box_vertices) */
1508     /*static int edge[NCUCEDGE*2];*/
1509     int             *edge;
1510     static const int hexcon[24] = {
1511         0, 9, 1, 19, 2, 15, 3, 21,
1512         4, 17, 5, 11, 6, 23, 7, 13,
1513         8, 20, 10, 18, 12, 16, 14, 22
1514     };
1515     int              e, i, j;
1516     gmx_bool         bFirst = TRUE;
1517
1518     snew(edge, NCUCEDGE*2);
1519
1520     if (bFirst)
1521     {
1522         e = 0;
1523         for (i = 0; i < 6; i++)
1524         {
1525             for (j = 0; j < 4; j++)
1526             {
1527                 edge[e++] = 4*i + j;
1528                 edge[e++] = 4*i + (j+1) % 4;
1529             }
1530         }
1531         for (i = 0; i < 12*2; i++)
1532         {
1533             edge[e++] = hexcon[i];
1534         }
1535
1536         bFirst = FALSE;
1537     }
1538
1539     return edge;
1540 }
1541
1542 void put_atoms_in_box_omp(int ePBC, matrix box, int natoms, rvec x[])
1543 {
1544     int t, nth;
1545     nth = gmx_omp_nthreads_get(emntDefault);
1546
1547 #pragma omp parallel for num_threads(nth) schedule(static)
1548     for (t = 0; t < nth; t++)
1549     {
1550         int offset, len;
1551
1552         offset = (natoms*t    )/nth;
1553         len    = (natoms*(t + 1))/nth - offset;
1554         put_atoms_in_box(ePBC, box, len, x + offset);
1555     }
1556 }
1557
1558 void put_atoms_in_box(int ePBC, matrix box, int natoms, rvec x[])
1559 {
1560     int npbcdim, i, m, d;
1561
1562     if (ePBC == epbcSCREW)
1563     {
1564         gmx_fatal(FARGS, "Sorry, %s pbc is not yet supported", epbc_names[ePBC]);
1565     }
1566
1567     if (ePBC == epbcXY)
1568     {
1569         npbcdim = 2;
1570     }
1571     else
1572     {
1573         npbcdim = 3;
1574     }
1575
1576     if (TRICLINIC(box))
1577     {
1578         for (i = 0; (i < natoms); i++)
1579         {
1580             for (m = npbcdim-1; m >= 0; m--)
1581             {
1582                 while (x[i][m] < 0)
1583                 {
1584                     for (d = 0; d <= m; d++)
1585                     {
1586                         x[i][d] += box[m][d];
1587                     }
1588                 }
1589                 while (x[i][m] >= box[m][m])
1590                 {
1591                     for (d = 0; d <= m; d++)
1592                     {
1593                         x[i][d] -= box[m][d];
1594                     }
1595                 }
1596             }
1597         }
1598     }
1599     else
1600     {
1601         for (i = 0; i < natoms; i++)
1602         {
1603             for (d = 0; d < npbcdim; d++)
1604             {
1605                 while (x[i][d] < 0)
1606                 {
1607                     x[i][d] += box[d][d];
1608                 }
1609                 while (x[i][d] >= box[d][d])
1610                 {
1611                     x[i][d] -= box[d][d];
1612                 }
1613             }
1614         }
1615     }
1616 }
1617
1618 void put_atoms_in_triclinic_unitcell(int ecenter, matrix box,
1619                                      int natoms, rvec x[])
1620 {
1621     rvec   box_center, shift_center;
1622     real   shm01, shm02, shm12, shift;
1623     int    i, m, d;
1624
1625     calc_box_center(ecenter, box, box_center);
1626
1627     /* The product of matrix shm with a coordinate gives the shift vector
1628        which is required determine the periodic cell position */
1629     shm01 = box[1][0]/box[1][1];
1630     shm02 = (box[1][1]*box[2][0] - box[2][1]*box[1][0])/(box[1][1]*box[2][2]);
1631     shm12 = box[2][1]/box[2][2];
1632
1633     clear_rvec(shift_center);
1634     for (d = 0; d < DIM; d++)
1635     {
1636         rvec_inc(shift_center, box[d]);
1637     }
1638     svmul(0.5, shift_center, shift_center);
1639     rvec_sub(box_center, shift_center, shift_center);
1640
1641     shift_center[0] = shm01*shift_center[1] + shm02*shift_center[2];
1642     shift_center[1] = shm12*shift_center[2];
1643     shift_center[2] = 0;
1644
1645     for (i = 0; (i < natoms); i++)
1646     {
1647         for (m = DIM-1; m >= 0; m--)
1648         {
1649             shift = shift_center[m];
1650             if (m == 0)
1651             {
1652                 shift += shm01*x[i][1] + shm02*x[i][2];
1653             }
1654             else if (m == 1)
1655             {
1656                 shift += shm12*x[i][2];
1657             }
1658             while (x[i][m]-shift < 0)
1659             {
1660                 for (d = 0; d <= m; d++)
1661                 {
1662                     x[i][d] += box[m][d];
1663                 }
1664             }
1665             while (x[i][m]-shift >= box[m][m])
1666             {
1667                 for (d = 0; d <= m; d++)
1668                 {
1669                     x[i][d] -= box[m][d];
1670                 }
1671             }
1672         }
1673     }
1674 }
1675
1676 const char *
1677 put_atoms_in_compact_unitcell(int ePBC, int ecenter, matrix box,
1678                               int natoms, rvec x[])
1679 {
1680     t_pbc pbc;
1681     rvec  box_center, dx;
1682     int   i;
1683
1684     set_pbc(&pbc, ePBC, box);
1685     calc_box_center(ecenter, box, box_center);
1686     for (i = 0; i < natoms; i++)
1687     {
1688         pbc_dx(&pbc, x[i], box_center, dx);
1689         rvec_add(box_center, dx, x[i]);
1690     }
1691
1692     return pbc.bLimitDistance ?
1693            "WARNING: Could not put all atoms in the compact unitcell\n"
1694            : NULL;
1695 }