src/gromacs/gmxlib/pbc.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.2.0
  11  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  12  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  13  * Copyright (c) 2001-2004, The GROMACS development team,
  14  * check out http://www.gromacs.org for more information.
  15
  16  * This program is free software; you can redistribute it and/or
  17  * modify it under the terms of the GNU General Public License
  18  * as published by the Free Software Foundation; either version 2
  19  * of the License, or (at your option) any later version.
  20  *
  21  * If you want to redistribute modifications, please consider that
  22  * scientific software is very special. Version control is crucial -
  23  * bugs must be traceable. We will be happy to consider code for
  24  * inclusion in the official distribution, but derived work must not
  25  * be called official GROMACS. Details are found in the README & COPYING
  26  * files - if they are missing, get the official version at www.gromacs.org.
  27  *
  28  * To help us fund GROMACS development, we humbly ask that you cite
  29  * the papers on the package - you can find them in the top README file.
  30  *
  31  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * GROningen Mixture of Alchemy and Childrens' Stories
  35  */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include <math.h>
  41 #include "sysstuff.h"
  42 #include "typedefs.h"
  43 #include "vec.h"
  44 #include "maths.h"
  45 #include "main.h"
  46 #include "pbc.h"
  47 #include "smalloc.h"
  48 #include "txtdump.h"
  49 #include "gmx_fatal.h"
  50 #include "names.h"
  51 #include "macros.h"
  52 #include "gmx_omp_nthreads.h"
  53
  54 /* Skip 0 so we have more chance of detecting if we forgot to call set_pbc. */
  55 enum {
  56     epbcdxRECTANGULAR = 1, epbcdxTRICLINIC,
  57     epbcdx2D_RECT,       epbcdx2D_TRIC,
  58     epbcdx1D_RECT,       epbcdx1D_TRIC,
  59     epbcdxSCREW_RECT,    epbcdxSCREW_TRIC,
  60     epbcdxNOPBC,         epbcdxUNSUPPORTED
  61 };
  62
  63 /* Margin factor for error message and correction if the box is too skewed */
  64 #define BOX_MARGIN         1.0010
  65 #define BOX_MARGIN_CORRECT 1.0005
  66
  67 int ePBC2npbcdim(int ePBC)
  68 {
  69     int npbcdim = 0;
  70
  71     switch (ePBC)
  72     {
  73         case epbcXYZ:   npbcdim = 3; break;
  74         case epbcXY:    npbcdim = 2; break;
  75         case epbcSCREW: npbcdim = 3; break;
  76         case epbcNONE:  npbcdim = 0; break;
  77         default: gmx_fatal(FARGS, "Unknown ePBC=%d in ePBC2npbcdim", ePBC);
  78     }
  79
  80     return npbcdim;
  81 }
  82
  83 int inputrec2nboundeddim(t_inputrec *ir)
  84 {
  85     if (ir->nwall == 2 && ir->ePBC == epbcXY)
  86     {
  87         return 3;
  88     }
  89     else
  90     {
  91         return ePBC2npbcdim(ir->ePBC);
  92     }
  93 }
  94
  95 void dump_pbc(FILE *fp, t_pbc *pbc)
  96 {
  97     rvec sum_box;
  98
  99     fprintf(fp, "ePBCDX = %d\n", pbc->ePBCDX);
 100     pr_rvecs(fp, 0, "box", pbc->box, DIM);
 101     pr_rvecs(fp, 0, "fbox_diag", &pbc->fbox_diag, 1);
 102     pr_rvecs(fp, 0, "hbox_diag", &pbc->hbox_diag, 1);
 103     pr_rvecs(fp, 0, "mhbox_diag", &pbc->mhbox_diag, 1);
 104     rvec_add(pbc->hbox_diag, pbc->mhbox_diag, sum_box);
 105     pr_rvecs(fp, 0, "sum of the above two", &sum_box, 1);
 106     fprintf(fp, "max_cutoff2 = %g\n", pbc->max_cutoff2);
 107     fprintf(fp, "bLimitDistance = %s\n", EBOOL(pbc->bLimitDistance));
 108     fprintf(fp, "limit_distance2 = %g\n", pbc->limit_distance2);
 109     fprintf(fp, "ntric_vec = %d\n", pbc->ntric_vec);
 110     if (pbc->ntric_vec > 0)
 111     {
 112         pr_ivecs(fp, 0, "tric_shift", pbc->tric_shift, pbc->ntric_vec, FALSE);
 113         pr_rvecs(fp, 0, "tric_vec", pbc->tric_vec, pbc->ntric_vec);
 114     }
 115 }
 116
 117 const char *check_box(int ePBC, matrix box)
 118 {
 119     const char *ptr;
 120
 121     if (ePBC == -1)
 122     {
 123         ePBC = guess_ePBC(box);
 124     }
 125
 126     if (ePBC == epbcNONE)
 127     {
 128         return NULL;
 129     }
 130
 131     if ((box[XX][YY] != 0) || (box[XX][ZZ] != 0) || (box[YY][ZZ] != 0))
 132     {
 133         ptr = "Only triclinic boxes with the first vector parallel to the x-axis and the second vector in the xy-plane are supported.";
 134     }
 135     else if (ePBC == epbcSCREW && (box[YY][XX] != 0 || box[ZZ][XX] != 0))
 136     {
 137         ptr = "The unit cell can not have off-diagonal x-components with screw pbc";
 138     }
 139     else if (fabs(box[YY][XX]) > BOX_MARGIN*0.5*box[XX][XX] ||
 140              (ePBC != epbcXY &&
 141               (fabs(box[ZZ][XX]) > BOX_MARGIN*0.5*box[XX][XX] ||
 142                fabs(box[ZZ][YY]) > BOX_MARGIN*0.5*box[YY][YY])))
 143     {
 144         ptr = "Triclinic box is too skewed.";
 145     }
 146     else
 147     {
 148         ptr = NULL;
 149     }
 150
 151     return ptr;
 152 }
 153
 154 real max_cutoff2(int ePBC, matrix box)
 155 {
 156     real min_hv2, min_ss;
 157
 158     /* Physical limitation of the cut-off
 159      * by half the length of the shortest box vector.
 160      */
 161     min_hv2 = min(0.25*norm2(box[XX]), 0.25*norm2(box[YY]));
 162     if (ePBC != epbcXY)
 163     {
 164         min_hv2 = min(min_hv2, 0.25*norm2(box[ZZ]));
 165     }
 166
 167     /* Limitation to the smallest diagonal element due to optimizations:
 168      * checking only linear combinations of single box-vectors (2 in x)
 169      * in the grid search and pbc_dx is a lot faster
 170      * than checking all possible combinations.
 171      */
 172     if (ePBC == epbcXY)
 173     {
 174         min_ss = min(box[XX][XX], box[YY][YY]);
 175     }
 176     else
 177     {
 178         min_ss = min(box[XX][XX], min(box[YY][YY]-fabs(box[ZZ][YY]), box[ZZ][ZZ]));
 179     }
 180
 181     return min(min_hv2, min_ss*min_ss);
 182 }
 183
 184 /* this one is mostly harmless... */
 185 static gmx_bool bWarnedGuess = FALSE;
 186
 187 int guess_ePBC(matrix box)
 188 {
 189     int ePBC;
 190
 191     if (box[XX][XX] > 0 && box[YY][YY] > 0 && box[ZZ][ZZ] > 0)
 192     {
 193         ePBC = epbcXYZ;
 194     }
 195     else if (box[XX][XX] > 0 && box[YY][YY] > 0 && box[ZZ][ZZ] == 0)
 196     {
 197         ePBC = epbcXY;
 198     }
 199     else if (box[XX][XX] == 0 && box[YY][YY] == 0 && box[ZZ][ZZ] == 0)
 200     {
 201         ePBC = epbcNONE;
 202     }
 203     else
 204     {
 205         if (!bWarnedGuess)
 206         {
 207             fprintf(stderr, "WARNING: Unsupported box diagonal %f %f %f, "
 208                     "will not use periodic boundary conditions\n\n",
 209                     box[XX][XX], box[YY][YY], box[ZZ][ZZ]);
 210             bWarnedGuess = TRUE;
 211         }
 212         ePBC = epbcNONE;
 213     }
 214
 215     if (debug)
 216     {
 217         fprintf(debug, "Guessed pbc = %s from the box matrix\n", epbc_names[ePBC]);
 218     }
 219
 220     return ePBC;
 221 }
 222
 223 static int correct_box_elem(FILE *fplog, int step, tensor box, int v, int d)
 224 {
 225     int shift, maxshift = 10;
 226
 227     shift = 0;
 228
 229     /* correct elem d of vector v with vector d */
 230     while (box[v][d] > BOX_MARGIN_CORRECT*0.5*box[d][d])
 231     {
 232         if (fplog)
 233         {
 234             fprintf(fplog, "Step %d: correcting invalid box:\n", step);
 235             pr_rvecs(fplog, 0, "old box", box, DIM);
 236         }
 237         rvec_dec(box[v], box[d]);
 238         shift--;
 239         if (fplog)
 240         {
 241             pr_rvecs(fplog, 0, "new box", box, DIM);
 242         }
 243         if (shift <= -maxshift)
 244         {
 245             gmx_fatal(FARGS,
 246                       "Box was shifted at least %d times. Please see log-file.",
 247                       maxshift);
 248         }
 249     }
 250     while (box[v][d] < -BOX_MARGIN_CORRECT*0.5*box[d][d])
 251     {
 252         if (fplog)
 253         {
 254             fprintf(fplog, "Step %d: correcting invalid box:\n", step);
 255             pr_rvecs(fplog, 0, "old box", box, DIM);
 256         }
 257         rvec_inc(box[v], box[d]);
 258         shift++;
 259         if (fplog)
 260         {
 261             pr_rvecs(fplog, 0, "new box", box, DIM);
 262         }
 263         if (shift >= maxshift)
 264         {
 265             gmx_fatal(FARGS,
 266                       "Box was shifted at least %d times. Please see log-file.",
 267                       maxshift);
 268         }
 269     }
 270
 271     return shift;
 272 }
 273
 274 gmx_bool correct_box(FILE *fplog, int step, tensor box, t_graph *graph)
 275 {
 276     int      zy, zx, yx, i;
 277     gmx_bool bCorrected;
 278
 279     /* check if the box still obeys the restrictions, if not, correct it */
 280     zy = correct_box_elem(fplog, step, box, ZZ, YY);
 281     zx = correct_box_elem(fplog, step, box, ZZ, XX);
 282     yx = correct_box_elem(fplog, step, box, YY, XX);
 283
 284     bCorrected = (zy || zx || yx);
 285
 286     if (bCorrected && graph)
 287     {
 288         /* correct the graph */
 289         for (i = graph->at_start; i < graph->at_end; i++)
 290         {
 291             graph->ishift[i][YY] -= graph->ishift[i][ZZ]*zy;
 292             graph->ishift[i][XX] -= graph->ishift[i][ZZ]*zx;
 293             graph->ishift[i][XX] -= graph->ishift[i][YY]*yx;
 294         }
 295     }
 296
 297     return bCorrected;
 298 }
 299
 300 int ndof_com(t_inputrec *ir)
 301 {
 302     int n = 0;
 303
 304     switch (ir->ePBC)
 305     {
 306         case epbcXYZ:
 307         case epbcNONE:
 308             n = 3;
 309             break;
 310         case epbcXY:
 311             n = (ir->nwall == 0 ? 3 : 2);
 312             break;
 313         case epbcSCREW:
 314             n = 1;
 315             break;
 316         default:
 317             gmx_incons("Unknown pbc in calc_nrdf");
 318     }
 319
 320     return n;
 321 }
 322
 323 static void low_set_pbc(t_pbc *pbc, int ePBC, ivec *dd_nc, matrix box)
 324 {
 325     int         order[5] = {0, -1, 1, -2, 2};
 326     int         ii, jj, kk, i, j, k, d, dd, jc, kc, npbcdim, shift;
 327     ivec        bPBC;
 328     real        d2old, d2new, d2new_c;
 329     rvec        trial, pos;
 330     gmx_bool    bXY, bUse;
 331     const char *ptr;
 332
 333     pbc->ePBC      = ePBC;
 334     pbc->ndim_ePBC = ePBC2npbcdim(ePBC);
 335
 336     copy_mat(box, pbc->box);
 337     pbc->bLimitDistance = FALSE;
 338     pbc->max_cutoff2    = 0;
 339     pbc->dim            = -1;
 340
 341     for (i = 0; (i < DIM); i++)
 342     {
 343         pbc->fbox_diag[i]  =  box[i][i];
 344         pbc->hbox_diag[i]  =  pbc->fbox_diag[i]*0.5;
 345         pbc->mhbox_diag[i] = -pbc->hbox_diag[i];
 346     }
 347
 348     ptr = check_box(ePBC, box);
 349     if (ePBC == epbcNONE)
 350     {
 351         pbc->ePBCDX = epbcdxNOPBC;
 352     }
 353     else if (ptr)
 354     {
 355         fprintf(stderr,   "Warning: %s\n", ptr);
 356         pr_rvecs(stderr, 0, "         Box", box, DIM);
 357         fprintf(stderr,   "         Can not fix pbc.\n");
 358         pbc->ePBCDX          = epbcdxUNSUPPORTED;
 359         pbc->bLimitDistance  = TRUE;
 360         pbc->limit_distance2 = 0;
 361     }
 362     else
 363     {
 364         if (ePBC == epbcSCREW && dd_nc)
 365         {
 366             /* This combinated should never appear here */
 367             gmx_incons("low_set_pbc called with screw pbc and dd_nc != NULL");
 368         }
 369
 370         npbcdim = 0;
 371         for (i = 0; i < DIM; i++)
 372         {
 373             if ((dd_nc && (*dd_nc)[i] > 1) || (ePBC == epbcXY && i == ZZ))
 374             {
 375                 bPBC[i] = 0;
 376             }
 377             else
 378             {
 379                 bPBC[i] = 1;
 380                 npbcdim++;
 381             }
 382         }
 383         switch (npbcdim)
 384         {
 385             case 1:
 386                 /* 1D pbc is not an mdp option and it is therefore only used
 387                  * with single shifts.
 388                  */
 389                 pbc->ePBCDX = epbcdx1D_RECT;
 390                 for (i = 0; i < DIM; i++)
 391                 {
 392                     if (bPBC[i])
 393                     {
 394                         pbc->dim = i;
 395                     }
 396                 }
 397                 for (i = 0; i < pbc->dim; i++)
 398                 {
 399                     if (pbc->box[pbc->dim][i] != 0)
 400                     {
 401                         pbc->ePBCDX = epbcdx1D_TRIC;
 402                     }
 403                 }
 404                 break;
 405             case 2:
 406                 pbc->ePBCDX = epbcdx2D_RECT;
 407                 for (i = 0; i < DIM; i++)
 408                 {
 409                     if (!bPBC[i])
 410                     {
 411                         pbc->dim = i;
 412                     }
 413                 }
 414                 for (i = 0; i < DIM; i++)
 415                 {
 416                     if (bPBC[i])
 417                     {
 418                         for (j = 0; j < i; j++)
 419                         {
 420                             if (pbc->box[i][j] != 0)
 421                             {
 422                                 pbc->ePBCDX = epbcdx2D_TRIC;
 423                             }
 424                         }
 425                     }
 426                 }
 427                 break;
 428             case 3:
 429                 if (ePBC != epbcSCREW)
 430                 {
 431                     if (TRICLINIC(box))
 432                     {
 433                         pbc->ePBCDX = epbcdxTRICLINIC;
 434                     }
 435                     else
 436                     {
 437                         pbc->ePBCDX = epbcdxRECTANGULAR;
 438                     }
 439                 }
 440                 else
 441                 {
 442                     pbc->ePBCDX = (box[ZZ][YY] == 0 ? epbcdxSCREW_RECT : epbcdxSCREW_TRIC);
 443                     if (pbc->ePBCDX == epbcdxSCREW_TRIC)
 444                     {
 445                         fprintf(stderr,
 446                                 "Screw pbc is not yet implemented for triclinic boxes.\n"
 447                                 "Can not fix pbc.\n");
 448                         pbc->ePBCDX = epbcdxUNSUPPORTED;
 449                     }
 450                 }
 451                 break;
 452             default:
 453                 gmx_fatal(FARGS, "Incorrect number of pbc dimensions with DD: %d",
 454                           npbcdim);
 455         }
 456         pbc->max_cutoff2 = max_cutoff2(ePBC, box);
 457
 458         if (pbc->ePBCDX == epbcdxTRICLINIC ||
 459             pbc->ePBCDX == epbcdx2D_TRIC ||
 460             pbc->ePBCDX == epbcdxSCREW_TRIC)
 461         {
 462             if (debug)
 463             {
 464                 pr_rvecs(debug, 0, "Box", box, DIM);
 465                 fprintf(debug, "max cutoff %.3f\n", sqrt(pbc->max_cutoff2));
 466             }
 467             pbc->ntric_vec = 0;
 468             /* We will only use single shifts, but we will check a few
 469              * more shifts to see if there is a limiting distance
 470              * above which we can not be sure of the correct distance.
 471              */
 472             for (kk = 0; kk < 5; kk++)
 473             {
 474                 k = order[kk];
 475                 if (!bPBC[ZZ] && k != 0)
 476                 {
 477                     continue;
 478                 }
 479                 for (jj = 0; jj < 5; jj++)
 480                 {
 481                     j = order[jj];
 482                     if (!bPBC[YY] && j != 0)
 483                     {
 484                         continue;
 485                     }
 486                     for (ii = 0; ii < 3; ii++)
 487                     {
 488                         i = order[ii];
 489                         if (!bPBC[XX] && i != 0)
 490                         {
 491                             continue;
 492                         }
 493                         /* A shift is only useful when it is trilinic */
 494                         if (j != 0 || k != 0)
 495                         {
 496                             d2old = 0;
 497                             d2new = 0;
 498                             for (d = 0; d < DIM; d++)
 499                             {
 500                                 trial[d] = i*box[XX][d] + j*box[YY][d] + k*box[ZZ][d];
 501                                 /* Choose the vector within the brick around 0,0,0 that
 502                                  * will become the shortest due to shift try.
 503                                  */
 504                                 if (d == pbc->dim)
 505                                 {
 506                                     trial[d] = 0;
 507                                     pos[d]   = 0;
 508                                 }
 509                                 else
 510                                 {
 511                                     if (trial[d] < 0)
 512                                     {
 513                                         pos[d] = min( pbc->hbox_diag[d], -trial[d]);
 514                                     }
 515                                     else
 516                                     {
 517                                         pos[d] = max(-pbc->hbox_diag[d], -trial[d]);
 518                                     }
 519                                 }
 520                                 d2old += sqr(pos[d]);
 521                                 d2new += sqr(pos[d] + trial[d]);
 522                             }
 523                             if (BOX_MARGIN*d2new < d2old)
 524                             {
 525                                 if (j < -1 || j > 1 || k < -1 || k > 1)
 526                                 {
 527                                     /* Check if there is a single shift vector
 528                                      * that decreases this distance even more.
 529                                      */
 530                                     jc = 0;
 531                                     kc = 0;
 532                                     if (j < -1 || j > 1)
 533                                     {
 534                                         jc = j/2;
 535                                     }
 536                                     if (k < -1 || k > 1)
 537                                     {
 538                                         kc = k/2;
 539                                     }
 540                                     d2new_c = 0;
 541                                     for (d = 0; d < DIM; d++)
 542                                     {
 543                                         d2new_c += sqr(pos[d] + trial[d]
 544                                                        - jc*box[YY][d] - kc*box[ZZ][d]);
 545                                     }
 546                                     if (d2new_c > BOX_MARGIN*d2new)
 547                                     {
 548                                         /* Reject this shift vector, as there is no a priori limit
 549                                          * to the number of shifts that decrease distances.
 550                                          */
 551                                         if (!pbc->bLimitDistance || d2new <  pbc->limit_distance2)
 552                                         {
 553                                             pbc->limit_distance2 = d2new;
 554                                         }
 555                                         pbc->bLimitDistance = TRUE;
 556                                     }
 557                                 }
 558                                 else
 559                                 {
 560                                     /* Check if shifts with one box vector less do better */
 561                                     bUse = TRUE;
 562                                     for (dd = 0; dd < DIM; dd++)
 563                                     {
 564                                         shift = (dd == 0 ? i : (dd == 1 ? j : k));
 565                                         if (shift)
 566                                         {
 567                                             d2new_c = 0;
 568                                             for (d = 0; d < DIM; d++)
 569                                             {
 570                                                 d2new_c += sqr(pos[d] + trial[d] - shift*box[dd][d]);
 571                                             }
 572                                             if (d2new_c <= BOX_MARGIN*d2new)
 573                                             {
 574                                                 bUse = FALSE;
 575                                             }
 576                                         }
 577                                     }
 578                                     if (bUse)
 579                                     {
 580                                         /* Accept this shift vector. */
 581                                         if (pbc->ntric_vec >= MAX_NTRICVEC)
 582                                         {
 583                                             fprintf(stderr, "\nWARNING: Found more than %d triclinic correction vectors, ignoring some.\n"
 584                                                     "  There is probably something wrong with your box.\n", MAX_NTRICVEC);
 585                                             pr_rvecs(stderr, 0, "         Box", box, DIM);
 586                                         }
 587                                         else
 588                                         {
 589                                             copy_rvec(trial, pbc->tric_vec[pbc->ntric_vec]);
 590                                             pbc->tric_shift[pbc->ntric_vec][XX] = i;
 591                                             pbc->tric_shift[pbc->ntric_vec][YY] = j;
 592                                             pbc->tric_shift[pbc->ntric_vec][ZZ] = k;
 593                                             pbc->ntric_vec++;
 594                                         }
 595                                     }
 596                                 }
 597                                 if (debug)
 598                                 {
 599                                     fprintf(debug, "  tricvec %2d = %2d %2d %2d  %5.2f %5.2f  %5.2f %5.2f %5.2f  %5.2f %5.2f %5.2f\n",
 600                                             pbc->ntric_vec, i, j, k,
 601                                             sqrt(d2old), sqrt(d2new),
 602                                             trial[XX], trial[YY], trial[ZZ],
 603                                             pos[XX], pos[YY], pos[ZZ]);
 604                                 }
 605                             }
 606                         }
 607                     }
 608                 }
 609             }
 610         }
 611     }
 612 }
 613
 614 void set_pbc(t_pbc *pbc, int ePBC, matrix box)
 615 {
 616     if (ePBC == -1)
 617     {
 618         ePBC = guess_ePBC(box);
 619     }
 620
 621     low_set_pbc(pbc, ePBC, NULL, box);
 622 }
 623
 624 t_pbc *set_pbc_dd(t_pbc *pbc, int ePBC,
 625                   gmx_domdec_t *dd, gmx_bool bSingleDir, matrix box)
 626 {
 627     ivec nc2;
 628     int  npbcdim, i;
 629
 630     if (dd == NULL)
 631     {
 632         npbcdim = DIM;
 633     }
 634     else
 635     {
 636         if (ePBC == epbcSCREW && dd->nc[XX] > 1)
 637         {
 638             /* The rotation has been taken care of during coordinate communication */
 639             ePBC = epbcXYZ;
 640         }
 641         npbcdim = 0;
 642         for (i = 0; i < DIM; i++)
 643         {
 644             if (dd->nc[i] <= (bSingleDir ? 1 : 2))
 645             {
 646                 nc2[i] = 1;
 647                 if (!(ePBC == epbcXY && i == ZZ))
 648                 {
 649                     npbcdim++;
 650                 }
 651             }
 652             else
 653             {
 654                 nc2[i] = dd->nc[i];
 655             }
 656         }
 657     }
 658
 659     if (npbcdim > 0)
 660     {
 661         low_set_pbc(pbc, ePBC, npbcdim < DIM ? &nc2 : NULL, box);
 662     }
 663
 664     return (npbcdim > 0 ? 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 void pbc_dx_d(const t_pbc *pbc, const dvec x1, const dvec x2, dvec dx)
1114 {
1115     int      i, j;
1116     dvec     dx_start, trial;
1117     double   d2min, d2trial;
1118     gmx_bool bRot;
1119
1120     dvec_sub(x1, x2, dx);
1121
1122     switch (pbc->ePBCDX)
1123     {
1124         case epbcdxRECTANGULAR:
1125         case epbcdx2D_RECT:
1126             for (i = 0; i < DIM; i++)
1127             {
1128                 if (i != pbc->dim)
1129                 {
1130                     while (dx[i] > pbc->hbox_diag[i])
1131                     {
1132                         dx[i] -= pbc->fbox_diag[i];
1133                     }
1134                     while (dx[i] <= pbc->mhbox_diag[i])
1135                     {
1136                         dx[i] += pbc->fbox_diag[i];
1137                     }
1138                 }
1139             }
1140             break;
1141         case epbcdxTRICLINIC:
1142         case epbcdx2D_TRIC:
1143             d2min = 0;
1144             for (i = DIM-1; i >= 0; i--)
1145             {
1146                 if (i != pbc->dim)
1147                 {
1148                     while (dx[i] > pbc->hbox_diag[i])
1149                     {
1150                         for (j = i; j >= 0; j--)
1151                         {
1152                             dx[j] -= pbc->box[i][j];
1153                         }
1154                     }
1155                     while (dx[i] <= pbc->mhbox_diag[i])
1156                     {
1157                         for (j = i; j >= 0; j--)
1158                         {
1159                             dx[j] += pbc->box[i][j];
1160                         }
1161                     }
1162                     d2min += dx[i]*dx[i];
1163                 }
1164             }
1165             if (d2min > pbc->max_cutoff2)
1166             {
1167                 copy_dvec(dx, dx_start);
1168                 /* Now try all possible shifts, when the distance is within max_cutoff
1169                  * it must be the shortest possible distance.
1170                  */
1171                 i = 0;
1172                 while ((d2min > pbc->max_cutoff2) && (i < pbc->ntric_vec))
1173                 {
1174                     for (j = 0; j < DIM; j++)
1175                     {
1176                         trial[j] = dx_start[j] + pbc->tric_vec[i][j];
1177                     }
1178                     d2trial = 0;
1179                     for (j = 0; j < DIM; j++)
1180                     {
1181                         if (j != pbc->dim)
1182                         {
1183                             d2trial += trial[j]*trial[j];
1184                         }
1185                     }
1186                     if (d2trial < d2min)
1187                     {
1188                         copy_dvec(trial, dx);
1189                         d2min = d2trial;
1190                     }
1191                     i++;
1192                 }
1193             }
1194             break;
1195         case epbcdxSCREW_RECT:
1196             /* The shift definition requires x first */
1197             bRot = FALSE;
1198             while (dx[XX] > pbc->hbox_diag[XX])
1199             {
1200                 dx[XX] -= pbc->fbox_diag[XX];
1201                 bRot    = !bRot;
1202             }
1203             while (dx[XX] <= pbc->mhbox_diag[XX])
1204             {
1205                 dx[XX] += pbc->fbox_diag[YY];
1206                 bRot    = !bRot;
1207             }
1208             if (bRot)
1209             {
1210                 /* Rotate around the x-axis in the middle of the box */
1211                 dx[YY] = pbc->box[YY][YY] - x1[YY] - x2[YY];
1212                 dx[ZZ] = pbc->box[ZZ][ZZ] - x1[ZZ] - x2[ZZ];
1213             }
1214             /* Normal pbc for y and z */
1215             for (i = YY; i <= ZZ; i++)
1216             {
1217                 while (dx[i] > pbc->hbox_diag[i])
1218                 {
1219                     dx[i] -= pbc->fbox_diag[i];
1220                 }
1221                 while (dx[i] <= pbc->mhbox_diag[i])
1222                 {
1223                     dx[i] += pbc->fbox_diag[i];
1224                 }
1225             }
1226             break;
1227         case epbcdxNOPBC:
1228         case epbcdxUNSUPPORTED:
1229             break;
1230         default:
1231             gmx_fatal(FARGS, "Internal error in pbc_dx, set_pbc has not been called");
1232             break;
1233     }
1234 }
1235
1236 gmx_bool image_rect(ivec xi, ivec xj, ivec box_size, real rlong2, int *shift, real *r2)
1237 {
1238     int     m, t;
1239     int     dx, b, b_2;
1240     real    dxr, rij2;
1241
1242     rij2 = 0.0;
1243     t    = 0;
1244     for (m = 0; (m < DIM); m++)
1245     {
1246         dx  = xi[m]-xj[m];
1247         t  *= DIM;
1248         b   = box_size[m];
1249         b_2 = b/2;
1250         if (dx < -b_2)
1251         {
1252             t  += 2;
1253             dx += b;
1254         }
1255         else if (dx > b_2)
1256         {
1257             dx -= b;
1258         }
1259         else
1260         {
1261             t += 1;
1262         }
1263         dxr   = dx;
1264         rij2 += dxr*dxr;
1265         if (rij2 >= rlong2)
1266         {
1267             return FALSE;
1268         }
1269     }
1270
1271     *shift = t;
1272     *r2    = rij2;
1273     return TRUE;
1274 }
1275
1276 gmx_bool image_cylindric(ivec xi, ivec xj, ivec box_size, real rlong2,
1277                          int *shift, real *r2)
1278 {
1279     int     m, t;
1280     int     dx, b, b_2;
1281     real    dxr, rij2;
1282
1283     rij2 = 0.0;
1284     t    = 0;
1285     for (m = 0; (m < DIM); m++)
1286     {
1287         dx  = xi[m]-xj[m];
1288         t  *= DIM;
1289         b   = box_size[m];
1290         b_2 = b/2;
1291         if (dx < -b_2)
1292         {
1293             t  += 2;
1294             dx += b;
1295         }
1296         else if (dx > b_2)
1297         {
1298             dx -= b;
1299         }
1300         else
1301         {
1302             t += 1;
1303         }
1304
1305         dxr   = dx;
1306         rij2 += dxr*dxr;
1307         if (m < ZZ)
1308         {
1309             if (rij2 >= rlong2)
1310             {
1311                 return FALSE;
1312             }
1313         }
1314     }
1315
1316     *shift = t;
1317     *r2    = rij2;
1318     return TRUE;
1319 }
1320
1321 void calc_shifts(matrix box, rvec shift_vec[])
1322 {
1323     int k, l, m, d, n, test;
1324
1325     n = 0;
1326     for (m = -D_BOX_Z; m <= D_BOX_Z; m++)
1327     {
1328         for (l = -D_BOX_Y; l <= D_BOX_Y; l++)
1329         {
1330             for (k = -D_BOX_X; k <= D_BOX_X; k++, n++)
1331             {
1332                 test = XYZ2IS(k, l, m);
1333                 if (n != test)
1334                 {
1335                     gmx_incons("inconsistent shift numbering");
1336                 }
1337                 for (d = 0; d < DIM; d++)
1338                 {
1339                     shift_vec[n][d] = k*box[XX][d] + l*box[YY][d] + m*box[ZZ][d];
1340                 }
1341             }
1342         }
1343     }
1344 }
1345
1346 void calc_box_center(int ecenter, matrix box, rvec box_center)
1347 {
1348     int d, m;
1349
1350     clear_rvec(box_center);
1351     switch (ecenter)
1352     {
1353         case ecenterTRIC:
1354             for (m = 0; (m < DIM); m++)
1355             {
1356                 for (d = 0; d < DIM; d++)
1357                 {
1358                     box_center[d] += 0.5*box[m][d];
1359                 }
1360             }
1361             break;
1362         case ecenterRECT:
1363             for (d = 0; d < DIM; d++)
1364             {
1365                 box_center[d] = 0.5*box[d][d];
1366             }
1367             break;
1368         case ecenterZERO:
1369             break;
1370         default:
1371             gmx_fatal(FARGS, "Unsupported value %d for ecenter", ecenter);
1372     }
1373 }
1374
1375 void calc_triclinic_images(matrix box, rvec img[])
1376 {
1377     int i;
1378
1379     /* Calculate 3 adjacent images in the xy-plane */
1380     copy_rvec(box[0], img[0]);
1381     copy_rvec(box[1], img[1]);
1382     if (img[1][XX] < 0)
1383     {
1384         svmul(-1, img[1], img[1]);
1385     }
1386     rvec_sub(img[1], img[0], img[2]);
1387
1388     /* Get the next 3 in the xy-plane as mirror images */
1389     for (i = 0; i < 3; i++)
1390     {
1391         svmul(-1, img[i], img[3+i]);
1392     }
1393
1394     /* Calculate the first 4 out of xy-plane images */
1395     copy_rvec(box[2], img[6]);
1396     if (img[6][XX] < 0)
1397     {
1398         svmul(-1, img[6], img[6]);
1399     }
1400     for (i = 0; i < 3; i++)
1401     {
1402         rvec_add(img[6], img[i+1], img[7+i]);
1403     }
1404
1405     /* Mirror the last 4 from the previous in opposite rotation */
1406     for (i = 0; i < 4; i++)
1407     {
1408         svmul(-1, img[6 + (2+i) % 4], img[10+i]);
1409     }
1410 }
1411
1412 void calc_compact_unitcell_vertices(int ecenter, matrix box, rvec vert[])
1413 {
1414     rvec img[NTRICIMG], box_center;
1415     int  n, i, j, tmp[4], d;
1416
1417     calc_triclinic_images(box, img);
1418
1419     n = 0;
1420     for (i = 2; i <= 5; i += 3)
1421     {
1422         tmp[0] = i-1;
1423         if (i == 2)
1424         {
1425             tmp[1] = 8;
1426         }
1427         else
1428         {
1429             tmp[1] = 6;
1430         }
1431         tmp[2] = (i+1) % 6;
1432         tmp[3] = tmp[1]+4;
1433         for (j = 0; j < 4; j++)
1434         {
1435             for (d = 0; d < DIM; d++)
1436             {
1437                 vert[n][d] = img[i][d]+img[tmp[j]][d]+img[tmp[(j+1)%4]][d];
1438             }
1439             n++;
1440         }
1441     }
1442     for (i = 7; i <= 13; i += 6)
1443     {
1444         tmp[0] = (i-7)/2;
1445         tmp[1] = tmp[0]+1;
1446         if (i == 7)
1447         {
1448             tmp[2] = 8;
1449         }
1450         else
1451         {
1452             tmp[2] = 10;
1453         }
1454         tmp[3] = i-1;
1455         for (j = 0; j < 4; j++)
1456         {
1457             for (d = 0; d < DIM; d++)
1458             {
1459                 vert[n][d] = img[i][d]+img[tmp[j]][d]+img[tmp[(j+1)%4]][d];
1460             }
1461             n++;
1462         }
1463     }
1464     for (i = 9; i <= 11; i += 2)
1465     {
1466         if (i == 9)
1467         {
1468             tmp[0] = 3;
1469         }
1470         else
1471         {
1472             tmp[0] = 0;
1473         }
1474         tmp[1] = tmp[0]+1;
1475         if (i == 9)
1476         {
1477             tmp[2] = 6;
1478         }
1479         else
1480         {
1481             tmp[2] = 12;
1482         }
1483         tmp[3] = i-1;
1484         for (j = 0; j < 4; j++)
1485         {
1486             for (d = 0; d < DIM; d++)
1487             {
1488                 vert[n][d] = img[i][d]+img[tmp[j]][d]+img[tmp[(j+1)%4]][d];
1489             }
1490             n++;
1491         }
1492     }
1493
1494     calc_box_center(ecenter, box, box_center);
1495     for (i = 0; i < NCUCVERT; i++)
1496     {
1497         for (d = 0; d < DIM; d++)
1498         {
1499             vert[i][d] = vert[i][d]*0.25+box_center[d];
1500         }
1501     }
1502 }
1503
1504 int *compact_unitcell_edges()
1505 {
1506     /* this is an index in vert[] (see calc_box_vertices) */
1507     /*static int edge[NCUCEDGE*2];*/
1508     int             *edge;
1509     static const int hexcon[24] = {
1510         0, 9, 1, 19, 2, 15, 3, 21,
1511         4, 17, 5, 11, 6, 23, 7, 13,
1512         8, 20, 10, 18, 12, 16, 14, 22
1513     };
1514     int              e, i, j;
1515     gmx_bool         bFirst = TRUE;
1516
1517     snew(edge, NCUCEDGE*2);
1518
1519     if (bFirst)
1520     {
1521         e = 0;
1522         for (i = 0; i < 6; i++)
1523         {
1524             for (j = 0; j < 4; j++)
1525             {
1526                 edge[e++] = 4*i + j;
1527                 edge[e++] = 4*i + (j+1) % 4;
1528             }
1529         }
1530         for (i = 0; i < 12*2; i++)
1531         {
1532             edge[e++] = hexcon[i];
1533         }
1534
1535         bFirst = FALSE;
1536     }
1537
1538     return edge;
1539 }
1540
1541 void put_atoms_in_box_omp(int ePBC, matrix box, int natoms, rvec x[])
1542 {
1543     int t, nth;
1544     nth = gmx_omp_nthreads_get(emntDefault);
1545
1546 #pragma omp parallel for num_threads(nth) schedule(static)
1547     for (t = 0; t < nth; t++)
1548     {
1549         int offset, len;
1550
1551         offset = (natoms*t    )/nth;
1552         len    = (natoms*(t + 1))/nth - offset;
1553         put_atoms_in_box(ePBC, box, len, x + offset);
1554     }
1555 }
1556
1557 void put_atoms_in_box(int ePBC, matrix box, int natoms, rvec x[])
1558 {
1559     int npbcdim, i, m, d;
1560
1561     if (ePBC == epbcSCREW)
1562     {
1563         gmx_fatal(FARGS, "Sorry, %s pbc is not yet supported", epbc_names[ePBC]);
1564     }
1565
1566     if (ePBC == epbcXY)
1567     {
1568         npbcdim = 2;
1569     }
1570     else
1571     {
1572         npbcdim = 3;
1573     }
1574
1575     if (TRICLINIC(box))
1576     {
1577         for (i = 0; (i < natoms); i++)
1578         {
1579             for (m = npbcdim-1; m >= 0; m--)
1580             {
1581                 while (x[i][m] < 0)
1582                 {
1583                     for (d = 0; d <= m; d++)
1584                     {
1585                         x[i][d] += box[m][d];
1586                     }
1587                 }
1588                 while (x[i][m] >= box[m][m])
1589                 {
1590                     for (d = 0; d <= m; d++)
1591                     {
1592                         x[i][d] -= box[m][d];
1593                     }
1594                 }
1595             }
1596         }
1597     }
1598     else
1599     {
1600         for (i = 0; i < natoms; i++)
1601         {
1602             for (d = 0; d < npbcdim; d++)
1603             {
1604                 while (x[i][d] < 0)
1605                 {
1606                     x[i][d] += box[d][d];
1607                 }
1608                 while (x[i][d] >= box[d][d])
1609                 {
1610                     x[i][d] -= box[d][d];
1611                 }
1612             }
1613         }
1614     }
1615 }
1616
1617 void put_atoms_in_triclinic_unitcell(int ecenter, matrix box,
1618                                      int natoms, rvec x[])
1619 {
1620     rvec   box_center, shift_center;
1621     real   shm01, shm02, shm12, shift;
1622     int    i, m, d;
1623
1624     calc_box_center(ecenter, box, box_center);
1625
1626     /* The product of matrix shm with a coordinate gives the shift vector
1627        which is required determine the periodic cell position */
1628     shm01 = box[1][0]/box[1][1];
1629     shm02 = (box[1][1]*box[2][0] - box[2][1]*box[1][0])/(box[1][1]*box[2][2]);
1630     shm12 = box[2][1]/box[2][2];
1631
1632     clear_rvec(shift_center);
1633     for (d = 0; d < DIM; d++)
1634     {
1635         rvec_inc(shift_center, box[d]);
1636     }
1637     svmul(0.5, shift_center, shift_center);
1638     rvec_sub(box_center, shift_center, shift_center);
1639
1640     shift_center[0] = shm01*shift_center[1] + shm02*shift_center[2];
1641     shift_center[1] = shm12*shift_center[2];
1642     shift_center[2] = 0;
1643
1644     for (i = 0; (i < natoms); i++)
1645     {
1646         for (m = DIM-1; m >= 0; m--)
1647         {
1648             shift = shift_center[m];
1649             if (m == 0)
1650             {
1651                 shift += shm01*x[i][1] + shm02*x[i][2];
1652             }
1653             else if (m == 1)
1654             {
1655                 shift += shm12*x[i][2];
1656             }
1657             while (x[i][m]-shift < 0)
1658             {
1659                 for (d = 0; d <= m; d++)
1660                 {
1661                     x[i][d] += box[m][d];
1662                 }
1663             }
1664             while (x[i][m]-shift >= box[m][m])
1665             {
1666                 for (d = 0; d <= m; d++)
1667                 {
1668                     x[i][d] -= box[m][d];
1669                 }
1670             }
1671         }
1672     }
1673 }
1674
1675 const char *
1676 put_atoms_in_compact_unitcell(int ePBC, int ecenter, matrix box,
1677                               int natoms, rvec x[])
1678 {
1679     t_pbc pbc;
1680     rvec  box_center, dx;
1681     int   i;
1682
1683     set_pbc(&pbc, ePBC, box);
1684     calc_box_center(ecenter, box, box_center);
1685     for (i = 0; i < natoms; i++)
1686     {
1687         pbc_dx(&pbc, x[i], box_center, dx);
1688         rvec_add(box_center, dx, x[i]);
1689     }
1690
1691     return pbc.bLimitDistance ?
1692            "WARNING: Could not put all atoms in the compact unitcell\n"
1693            : NULL;
1694 }