src/gmxlib/sortwater.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 3.1
  11  * Copyright (c) 1991-2001, University of Groningen, The Netherlands
  12  * This program is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU General Public License
  14  * as published by the Free Software Foundation; either version 2
  15  * of the License, or (at your option) any later version.
  16  *
  17  * If you want to redistribute modifications, please consider that
  18  * scientific software is very special. Version control is crucial -
  19  * bugs must be traceable. We will be happy to consider code for
  20  * inclusion in the official distribution, but derived work must not
  21  * be called official GROMACS. Details are found in the README & COPYING
  22  * files - if they are missing, get the official version at www.gromacs.org.
  23  *
  24  * To help us fund GROMACS development, we humbly ask that you cite
  25  * the papers on the package - you can find them in the top README file.
  26  *
  27  * For more info, check our website at http://www.gromacs.org
  28  *
  29  * And Hey:
  30  * Gnomes, ROck Monsters And Chili Sauce
  31  */
  32
  33 #include "typedefs.h"
  34 #include "random.h"
  35 #include "smalloc.h"
  36 #include "vec.h"
  37 #include "sortwater.h"
  38
  39 static rvec   *xptr,box_1;
  40 static int    nwat;
  41 static matrix BOX;
  42 static ivec   NBOX;
  43
  44 void randwater(int astart,int nwater,int nwatom,rvec x[],rvec v[],int *seed)
  45 {
  46   int  i,j,wi,wj,*tab;
  47   rvec buf;
  48
  49   snew(tab,nwater);
  50   for(i=0; (i<nwater); i++)
  51     tab[i]=i;
  52   for(j=0; (j<23*nwater); j++) {
  53     wi = (int) (nwater*rando(seed)) % nwater;
  54     do {
  55       wj = (int) (nwater*rando(seed)) % nwater;
  56     } while (wi == wj);
  57     wi = astart+wi*nwatom;
  58     wj = astart+wj*nwatom;
  59     /* Swap coords for wi and wj */
  60     for(i=0; (i<nwatom); i++) {
  61       copy_rvec(x[wi+i],buf);
  62       copy_rvec(x[wj+i],x[wi+i]);
  63       copy_rvec(buf,x[wj+i]);
  64       if (v) {
  65         copy_rvec(v[wi+i],buf);
  66         copy_rvec(v[wj+i],v[wi+i]);
  67         copy_rvec(buf,v[wj+i]);
  68       }
  69     }
  70   }
  71   sfree(tab);
  72 }
  73
  74 static int rvcomp(const void *a,const void *b)
  75 {
  76   int aa,bb;
  77
  78   aa = nwat*(*((int *)a));
  79   bb = nwat*(*((int *)b));
  80   if (xptr[aa][XX] < xptr[bb][XX])
  81     return -1;
  82   else if (xptr[aa][XX] > xptr[bb][XX])
  83     return 1;
  84   else
  85     return 0;
  86 }
  87
  88 static int block_index(rvec x,ivec nbox)
  89 {
  90   ivec ixyz;
  91   int  m;
  92
  93   for(m=0; (m<DIM); m++)
  94     ixyz[m] = ((int)((1+x[m]*box_1[m])*nbox[m])) % nbox[m];
  95   return ixyz[XX]*nbox[YY]*nbox[ZZ]+ixyz[YY]*nbox[ZZ]+ixyz[ZZ];
  96 }
  97
  98 static int blockcomp(const void *a,const void *b)
  99 {
 100   int aa,bb,aind,bind;
 101
 102   aa   = nwat*(*((int *)a));
 103   bb   = nwat*(*((int *)b));
 104
 105   aind = block_index(xptr[aa],NBOX);
 106   bind = block_index(xptr[bb],NBOX);
 107
 108   if (aind == bind) {
 109     if (xptr[aa][XX] < xptr[bb][XX])
 110       return -1;
 111     else if (xptr[aa][XX] > xptr[bb][XX])
 112       return 1;
 113     else
 114       return 0;
 115   }
 116   else
 117     return aind-bind;
 118 }
 119
 120 static void lo_sortwater(int astart,int nwater,int nwatom,rvec x[],rvec v[],
 121                          bool bBlock)
 122 {
 123   int  i,j,i0,rvi;
 124   int  *rvindex;
 125   rvec *tmp;
 126
 127   /* Sort indices to rvecs */
 128   snew(rvindex,nwater);
 129   for(i=0; (i<nwater); i++)
 130     rvindex[i] = i;
 131   xptr = x+astart;
 132   nwat = nwatom;
 133
 134   qsort(rvindex,nwater,sizeof(rvindex[0]),bBlock ? blockcomp : rvcomp);
 135   if (debug)
 136     for(i=0; (i<nwater); i++) {
 137       rvi = rvindex[i]*nwatom;
 138       fprintf(debug,"rvindex[%5d] = %5d (x = %8.3f  %8.3f  %8.3f)\n",
 139               i,rvi,x[astart+rvi][XX],x[astart+rvi][YY],x[astart+rvi][ZZ]);
 140     }
 141   snew(tmp,nwater*nwatom);
 142
 143   for(i=0; (i<nwater); i++) {
 144     i0 = astart+nwatom*rvindex[i];
 145     for(j=0; (j<nwatom); j++)
 146       copy_rvec(x[i0+j],tmp[nwatom*i+j]);
 147   }
 148   for(i=0; (i<nwater*nwatom); i++)
 149     copy_rvec(tmp[i],x[astart+i]);
 150
 151   for(i=0; (i<nwater); i++) {
 152     i0 = astart+nwatom*rvindex[i];
 153     for(j=0; (j<nwatom); j++)
 154       copy_rvec(v[i0+j],tmp[nwatom*i+j]);
 155   }
 156   for(i=0; (i<nwater*nwatom); i++)
 157     copy_rvec(tmp[i],v[astart+i]);
 158
 159   sfree(tmp);
 160   sfree(rvindex);
 161 }
 162
 163 void sortwater(int astart,int nwater,int nwatom,rvec x[],rvec v[])
 164 {
 165   lo_sortwater(astart,nwater,nwatom,x,v,FALSE);
 166 }
 167
 168 static void factorize(int nn,int fac[])
 169 {
 170   int i,n=nn;
 171
 172   for(i=0; (i<=n); i++)
 173     fac[i] = 0;
 174   fac[1] = 1;
 175   for(i=2; (i<=n); ) {
 176     if ((n % i) == 0) {
 177       fac[i]++;
 178       n = n/i;
 179     }
 180     else
 181       i++;
 182   }
 183   if (debug) {
 184     fprintf(debug,"Factorizing %d into primes:\n",nn);
 185     for(i=2; (i<=nn); i++) {
 186       if (fac[i])
 187         fprintf(debug,"%d ^ %d\n",i,fac[i]);
 188     }
 189   }
 190 }
 191
 192 static int ipow(int base,int exp)
 193 {
 194   int i,ip;
 195
 196   for(ip=1,i=0; (i<exp); i++) {
 197     ip *= base;
 198   }
 199   return ip;
 200 }
 201
 202 static int iv_comp(const void *a,const void *b)
 203 {
 204   int *ia,*ib;
 205
 206   ia = (int *)a;
 207   ib = (int *)b;
 208   if (ia[XX] != ib[XX])
 209     return (ia[XX] - ib[XX]);
 210   else if (ia[YY] != ib[YY])
 211     return (ia[YY] - ib[YY]);
 212   else
 213     return (ia[ZZ] - ib[ZZ]);
 214 }
 215
 216 static int add_bb(ivec BB[],int n,ivec b)
 217 {
 218 #define SWPX(vv,xx,yy) { int tmp; tmp=vv[xx]; vv[xx] = vv[yy]; vv[yy] = tmp; }
 219   copy_ivec(b,BB[n++]); /* x y z */
 220   SWPX(b,XX,YY);
 221   copy_ivec(b,BB[n++]); /* y x z */
 222   SWPX(b,XX,ZZ);
 223   copy_ivec(b,BB[n++]); /* z x y */
 224   SWPX(b,XX,YY);
 225   copy_ivec(b,BB[n++]); /* x z y */
 226   SWPX(b,XX,ZZ);
 227   copy_ivec(b,BB[n++]); /* y z x */
 228   SWPX(b,XX,YY);
 229   copy_ivec(b,BB[n++]); /* z y x */
 230   SWPX(b,XX,ZZ); /* Back to normal */
 231 #undef SWPX
 232   return n;
 233 }
 234
 235 static real box_weight(ivec nbox,matrix box)
 236 {
 237   rvec lx;
 238   int  m;
 239
 240   /* Calculate area of subbox */
 241   for(m=0; (m<DIM); m++)
 242     lx[m] = box[m][m]/nbox[m];
 243   return 2*(lx[XX]*lx[YY]+lx[XX]*lx[ZZ]+lx[YY]*lx[ZZ]);
 244 }
 245
 246 static int w_comp(const void *a,const void *b)
 247 {
 248   int *ia,*ib;
 249   real wa,wb;
 250
 251   ia = (int *)a;
 252   ib = (int *)b;
 253
 254   wa = box_weight(ia,BOX);
 255   wb = box_weight(ib,BOX);
 256   if (fabs(wa - wb) < 1e-4)
 257     return (iiprod(ia,ia) - iiprod(ib,ib));
 258   else if (wa < wb)
 259     return -1;
 260   else
 261     return 1;
 262 }
 263
 264 static void buildbox(int nnode,ivec nbox,matrix box)
 265 {
 266   ivec *BB,bxyz;
 267   int  i,j,m,n,n3,ny,*fx,*fy,nbb;
 268
 269   n3 = ipow(nnode,3)*6;
 270   snew(BB,n3);
 271   nbb=0;
 272   snew(fx,nnode+1);
 273   snew(fy,nnode+1);
 274   factorize(nnode,fx);
 275   for(i=0; (i<=nnode); i++) {
 276     for(m=1; (m<=fx[i]); m++) {
 277       bxyz[XX] = ipow(i,m);
 278       ny = nnode/bxyz[XX];
 279       factorize(ny,fy);
 280       for(j=0; (j<=ny); j++) {
 281         for(n=1; (n<=fy[j]); n++) {
 282           bxyz[YY] = ipow(j,n);
 283           bxyz[ZZ] = ny/bxyz[YY];
 284           if (bxyz[ZZ] > 0) {
 285             nbb = add_bb(BB,nbb,bxyz);
 286           }
 287         }
 288       }
 289     }
 290   }
 291   /* Sort boxes and remove doubles */
 292   qsort(BB,nbb,sizeof(BB[0]),iv_comp);
 293   j = 0;
 294   for(i=1; (i<nbb); i++) {
 295     if ((BB[i][XX] != BB[j][XX]) ||
 296         (BB[i][YY] != BB[j][YY]) ||
 297         (BB[i][ZZ] != BB[j][ZZ])) {
 298       j++;
 299       copy_ivec(BB[i],BB[j]);
 300     }
 301   }
 302   nbb = ++j;
 303   /* Sort boxes according to weight */
 304   copy_mat(box,BOX);
 305   qsort(BB,nbb,sizeof(BB[0]),w_comp);
 306   for(i=0; (i<nbb); i++) {
 307     fprintf(stderr,"nbox = %2d %2d %2d [ prod %3d ] area = %12.5f (nm^2)\n",
 308             BB[i][XX],BB[i][YY],BB[i][ZZ],
 309             BB[i][XX]*BB[i][YY]*BB[i][ZZ],
 310             box_weight(BB[i],box));
 311   }
 312   copy_ivec(BB[0],nbox);
 313   sfree(BB);
 314   sfree(fy);
 315   sfree(fx);
 316 }
 317
 318 void mkcompact(int astart,int nwater,int nwatom,rvec x[],rvec v[],
 319                int nnode,matrix box)
 320 {
 321   /* Make a compact configuration for each processor.
 322    * Divide the computational box in near cubic boxes and spread them
 323    * evenly over processors.
 324    */
 325 /*   ivec nbox; */
 326   int  m;
 327
 328   if (nnode <= 1)
 329     return;
 330
 331   buildbox(nnode,NBOX,box);
 332   /* copy_ivec(nbox,NBOX); */
 333   for(m=0; (m<DIM); m++)
 334     box_1[m] = 1.0/box[m][m];
 335
 336   lo_sortwater(astart,nwater,nwatom,x,v,TRUE);
 337 }