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