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