src/gmxlib/pbc.c

   1 /*
   2  * $Id$
   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 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * Green Red Orange Magenta Azure Cyan Skyblue
  28  */
  29 static char *SRCID_pbc_c = "$Id$";
  30
  31 #include <math.h>
  32 #include "sysstuff.h"
  33 #include "typedefs.h"
  34 #include "vec.h"
  35 #include "maths.h"
  36 #include "main.h"
  37 #include "pbc.h"
  38 #include "smalloc.h"
  39 #include "txtdump.h"
  40 #include "fatal.h"
  41
  42 /*****************************************
  43  *     PERIODIC BOUNDARY CONDITIONS
  44  *****************************************/
  45
  46 /* Global variables */
  47 static bool bInit=FALSE;
  48 static bool bTrunc;
  49 /* static bool bTriclinic; */
  50 static real side,side_2,side3_4;
  51 /* static real diag; */
  52 static rvec gl_fbox,gl_hbox,gl_mhbox;
  53
  54 void init_pbc(matrix box,bool bTruncOct)
  55 {
  56   int i;
  57 /*   int j; */
  58
  59   bInit   = TRUE;
  60   bTrunc  = bTruncOct;
  61   side    = box[0][0];
  62   side_2  = 0.5*side;
  63   side3_4 = 0.75*side;
  64   /*
  65   diag    = side_2*sqrt(3.0);
  66   bTriclinic=FALSE;
  67   for(i=0; (i<DIM); i++)
  68     for(j=0; (j<DIM); j++)
  69       if ((i != j) && (box[i][j] != 0.0))
  70         bTriclinic=TRUE;
  71         */
  72   for(i=0; (i<DIM); i++) {
  73     gl_fbox[i]  =  box[i][i];
  74     gl_hbox[i]  =  gl_fbox[i]*0.5;
  75     gl_mhbox[i] = -gl_hbox[i];
  76   }
  77 }
  78
  79 static void do_trunc(rvec x)
  80 {
  81   int i;
  82
  83   if (fabs(x[XX]-side_2)+fabs(x[YY]-side_2)+fabs(x[ZZ]-side_2) > side3_4)
  84     for(i=0; (i<DIM); i++)
  85       x[i] -= sign(side_2,x[i]);
  86 }
  87
  88 void pbc_dx(rvec x1, rvec x2, rvec dx)
  89 {
  90   int i;
  91
  92   if (!bInit)
  93     fatal_error(0,"pbc_dx called before init_pbc");
  94   rvec_sub(x1,x2,dx);
  95   for(i=0; (i<DIM); i++) {
  96     if (dx[i] > gl_hbox[i])
  97       dx[i] -= gl_fbox[i];
  98     else if (dx[i] <= gl_mhbox[i])
  99       dx[i] += gl_fbox[i];
 100   }
 101   if (bTrunc)
 102     do_trunc(dx);
 103 }
 104
 105 bool image_rect(ivec xi,ivec xj,ivec box_size,real rlong2,int *shift,real *r2)
 106 {
 107   int   m,t;
 108   int   dx,b,b_2;
 109   real  dxr,rij2;
 110
 111   rij2=0.0;
 112   t=0;
 113   for(m=0; (m<DIM); m++) {
 114     dx=xi[m]-xj[m];
 115     t*=DIM;
 116     b=box_size[m];
 117     b_2=b/2;
 118     if (dx < -b_2) {
 119       t+=2;
 120       dx+=b;
 121     }
 122     else if (dx > b_2)
 123       dx-=b;
 124     else
 125       t+=1;
 126     dxr=dx;
 127     rij2+=dxr*dxr;
 128     if (rij2 >= rlong2)
 129       return FALSE;
 130   }
 131
 132   *shift = t;
 133   *r2 = rij2;
 134   return TRUE;
 135 }
 136
 137 bool image_cylindric(ivec xi,ivec xj,ivec box_size,real rlong2,
 138                      int *shift,real *r2)
 139 {
 140   int   m,t;
 141   int   dx,b,b_2;
 142   real  dxr,rij2;
 143
 144   rij2=0.0;
 145   t=0;
 146   for(m=0; (m<DIM); m++) {
 147     dx=xi[m]-xj[m];
 148     t*=DIM;
 149     b=box_size[m];
 150     b_2=b/2;
 151     if (dx < -b_2) {
 152       t+=2;
 153       dx+=b;
 154     }
 155     else if (dx > b_2)
 156       dx-=b;
 157     else
 158       t+=1;
 159
 160     dxr=dx;
 161     rij2+=dxr*dxr;
 162     if (m < ZZ) {
 163       if (rij2 >= rlong2)
 164         return FALSE;
 165     }
 166   }
 167
 168   *shift = t;
 169   *r2 = rij2;
 170   return TRUE;
 171 }
 172
 173 bool image_tri(ivec xi,ivec xj,imatrix box,real rlong2,int *shift,real *r2)
 174 {
 175   fatal_error(0,"image_tri not implemented");
 176
 177   return FALSE; /* To make the compiler happy */
 178 }
 179
 180 void calc_shifts(matrix box,rvec box_size,rvec shift_vec[],bool bTruncOct)
 181 {
 182   int k,l,m,d,n,test;
 183
 184   init_pbc(box,bTruncOct);
 185   for (m=0; (m<DIM); m++)
 186     box_size[m]=box[m][m];
 187
 188   n=0;
 189   for(k = -D_BOX; k <= D_BOX; k++)
 190     for(l = -D_BOX; l <= D_BOX; l++)
 191       for(m = -D_BOX; m <= D_BOX; m++,n++) {
 192         test=XYZ2IS(k,l,m);
 193         if (n != test)
 194           fprintf(stdlog,"n=%d, test=%d\n",n,test);
 195         if (bTrunc && (k!=0) && (l!=0) && (m!=0)) {
 196           /* Truncated edges... */
 197           shift_vec[n][XX]=k*side_2;
 198           shift_vec[n][YY]=l*side_2;
 199           shift_vec[n][ZZ]=m*side_2;
 200         }
 201         else if (!bTrunc || (((k+l+m) % 2)!=0) || ((k==0) && (l==0) && (m==0)))
 202           for(d = 0; d < DIM; d++)
 203             shift_vec[n][d]=k*box[XX][d] + l*box[YY][d] + m*box[ZZ][d];
 204       }
 205 }
 206
 207 void calc_cgcm(FILE *log,int cg0,int cg1,t_block *cgs,
 208                rvec pos[],rvec cg_cm[])
 209 {
 210   int  icg,ai,k,k0,k1,d;
 211   rvec cg;
 212   real nrcg,inv_ncg;
 213   atom_id *cga,*cgindex;
 214
 215 #ifdef DEBUG
 216   fprintf(log,"Calculating centre of geometry for charge groups %d to %d\n",
 217           cg0,cg1);
 218 #endif
 219   cga     = cgs->a;
 220   cgindex = cgs->index;
 221
 222   /* Compute the center of geometry for all charge groups */
 223   for(icg=cg0; (icg<cg1); icg++) {
 224     k0      = cgindex[icg];
 225     k1      = cgindex[icg+1];
 226     nrcg    = k1-k0;
 227     if (nrcg == 1) {
 228       ai = cga[k0];
 229       copy_rvec(pos[ai],cg_cm[icg]);
 230     }
 231     else {
 232       inv_ncg = 1.0/nrcg;
 233
 234       clear_rvec(cg);
 235       for(k=k0; (k<k1); k++)  {
 236         ai     = cga[k];
 237         for(d=0; (d<DIM); d++)
 238           cg[d] += pos[ai][d];
 239       }
 240       for(d=0; (d<DIM); d++)
 241         cg_cm[icg][d] = inv_ncg*cg[d];
 242     }
 243   }
 244 }
 245
 246 void put_charge_groups_in_box(FILE *log,int cg0,int cg1,bool bTruncOct,
 247                               matrix box,rvec box_size,t_block *cgs,
 248                               rvec pos[],rvec shift_vec[],rvec cg_cm[])
 249
 250 {
 251   int  icg,ai,k,k0,k1,d;
 252   rvec cg;
 253   real nrcg,inv_ncg;
 254   atom_id *cga,*cgindex;
 255
 256 #ifdef DEBUG
 257   fprintf(log,"Putting cgs %d to %d in box\n",cg0,cg1);
 258 #endif
 259   cga     = cgs->a;
 260   cgindex = cgs->index;
 261
 262   for(icg=cg0; (icg<cg1); icg++) {
 263     /* First compute the center of geometry for this charge group */
 264     k0      = cgindex[icg];
 265     k1      = cgindex[icg+1];
 266     nrcg    = k1-k0;
 267     inv_ncg = 1.0/nrcg;
 268
 269     clear_rvec(cg);
 270     for(k=k0; (k<k1); k++)  {
 271       ai     = cga[k];
 272       for(d=0; d<DIM; d++)
 273         cg[d] += inv_ncg*pos[ai][d];
 274     }
 275     /* Now check pbc for this cg */
 276     for(d=0; d<DIM; d++) {
 277       while(cg[d] < 0) {
 278         cg[d] += box[d][d];
 279         for(k=k0; (k<k1); k++)
 280           pos[cga[k]][d] += box[d][d];
 281       }
 282       while(cg[d] >= box[d][d]) {
 283         cg[d] -= box[d][d];
 284         for(k=k0; (k<k1); k++)
 285           pos[cga[k]][d] -= box[d][d];
 286       }
 287       cg_cm[icg][d] = cg[d];
 288     }
 289 #ifdef DEBUG_PBC
 290     for(d=0; (d<DIM); d++) {
 291       if ((cg_cm[icg][d] < 0) || (cg_cm[icg][d] >= box[d][d]))
 292         fatal_error(0,"cg_cm[%d] = %15f  %15f  %15f\n"
 293                     "box = %15f  %15f  %15f\n",
 294                     icg,cg_cm[icg][XX],cg_cm[icg][YY],cg_cm[icg][ZZ],
 295                     box[XX][XX],box[YY][YY],box[ZZ][ZZ]);
 296     }
 297 #endif
 298   }
 299 }
 300
 301 void put_atoms_in_box(int natoms,matrix box,rvec x[])
 302 {
 303   int i,m;
 304
 305   for(i=0; (i<natoms); i++)
 306     for(m=0; m < DIM; m++) {
 307       while (x[i][m] < 0)
 308         x[i][m] += box[m][m];
 309       while (x[i][m] >= box[m][m])
 310         x[i][m] -= box[m][m];
 311     }
 312 }