src/gromacs/pbcutil/pbc.cpp

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