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