src/gmxlib/sortwater.c

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