src/gromacs/gmxlib/pbc.c

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