src/gmxlib/pbc.c

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