src/gmxlib/pbc.c

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