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