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