src/gromacs/pbcutil/pbc.c

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