src/gmxlib/sortwater.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 "typedefs.h"
  43 #include "random.h"
  44 #include "smalloc.h"
  45 #include "vec.h"
  46 #include "sortwater.h"
  47
  48 static rvec   *xptr,box_1;
  49 static int    nwat;
  50 static matrix BOX;
  51 static ivec   NBOX;
  52
  53 void randwater(int astart,int nwater,int nwatom,rvec x[],rvec v[],int *seed)
  54 {
  55   int  i,j,wi,wj,*tab;
  56   rvec buf;
  57
  58   snew(tab,nwater);
  59   for(i=0; (i<nwater); i++)
  60     tab[i]=i;
  61   for(j=0; (j<23*nwater); j++) {
  62     wi = (int) (nwater*rando(seed)) % nwater;
  63     do {
  64       wj = (int) (nwater*rando(seed)) % nwater;
  65     } while (wi == wj);
  66     wi = astart+wi*nwatom;
  67     wj = astart+wj*nwatom;
  68     /* Swap coords for wi and wj */
  69     for(i=0; (i<nwatom); i++) {
  70       copy_rvec(x[wi+i],buf);
  71       copy_rvec(x[wj+i],x[wi+i]);
  72       copy_rvec(buf,x[wj+i]);
  73       if (v) {
  74         copy_rvec(v[wi+i],buf);
  75         copy_rvec(v[wj+i],v[wi+i]);
  76         copy_rvec(buf,v[wj+i]);
  77       }
  78     }
  79   }
  80   sfree(tab);
  81 }
  82
  83 static int rvcomp(const void *a,const void *b)
  84 {
  85   int aa,bb;
  86
  87   aa = nwat*(*((int *)a));
  88   bb = nwat*(*((int *)b));
  89   if (xptr[aa][XX] < xptr[bb][XX])
  90     return -1;
  91   else if (xptr[aa][XX] > xptr[bb][XX])
  92     return 1;
  93   else
  94     return 0;
  95 }
  96
  97 static int block_index(rvec x,ivec nbox)
  98 {
  99   ivec ixyz;
 100   int  m;
 101
 102   for(m=0; (m<DIM); m++)
 103     ixyz[m] = ((int)((1+x[m]*box_1[m])*nbox[m])) % nbox[m];
 104   return ixyz[XX]*nbox[YY]*nbox[ZZ]+ixyz[YY]*nbox[ZZ]+ixyz[ZZ];
 105 }
 106
 107 static int blockcomp(const void *a,const void *b)
 108 {
 109   int aa,bb,aind,bind;
 110
 111   aa   = nwat*(*((int *)a));
 112   bb   = nwat*(*((int *)b));
 113
 114   aind = block_index(xptr[aa],NBOX);
 115   bind = block_index(xptr[bb],NBOX);
 116
 117   if (aind == bind) {
 118     if (xptr[aa][XX] < xptr[bb][XX])
 119       return -1;
 120     else if (xptr[aa][XX] > xptr[bb][XX])
 121       return 1;
 122     else
 123       return 0;
 124   }
 125   else
 126     return aind-bind;
 127 }
 128
 129 static void lo_sortwater(int astart,int nwater,int nwatom,rvec x[],rvec v[],
 130                          gmx_bool bBlock)
 131 {
 132   int  i,j,i0,rvi;
 133   int  *rvindex;
 134   rvec *tmp;
 135
 136   /* Sort indices to rvecs */
 137   snew(rvindex,nwater);
 138   for(i=0; (i<nwater); i++)
 139     rvindex[i] = i;
 140   xptr = x+astart;
 141   nwat = nwatom;
 142
 143   qsort(rvindex,nwater,sizeof(rvindex[0]),bBlock ? blockcomp : rvcomp);
 144   if (debug)
 145     for(i=0; (i<nwater); i++) {
 146       rvi = rvindex[i]*nwatom;
 147       fprintf(debug,"rvindex[%5d] = %5d (x = %8.3f  %8.3f  %8.3f)\n",
 148               i,rvi,x[astart+rvi][XX],x[astart+rvi][YY],x[astart+rvi][ZZ]);
 149     }
 150   snew(tmp,nwater*nwatom);
 151
 152   for(i=0; (i<nwater); i++) {
 153     i0 = astart+nwatom*rvindex[i];
 154     for(j=0; (j<nwatom); j++)
 155       copy_rvec(x[i0+j],tmp[nwatom*i+j]);
 156   }
 157   for(i=0; (i<nwater*nwatom); i++)
 158     copy_rvec(tmp[i],x[astart+i]);
 159
 160   for(i=0; (i<nwater); i++) {
 161     i0 = astart+nwatom*rvindex[i];
 162     for(j=0; (j<nwatom); j++)
 163       copy_rvec(v[i0+j],tmp[nwatom*i+j]);
 164   }
 165   for(i=0; (i<nwater*nwatom); i++)
 166     copy_rvec(tmp[i],v[astart+i]);
 167
 168   sfree(tmp);
 169   sfree(rvindex);
 170 }
 171
 172 void sortwater(int astart,int nwater,int nwatom,rvec x[],rvec v[])
 173 {
 174   lo_sortwater(astart,nwater,nwatom,x,v,FALSE);
 175 }
 176
 177 static void factorize(int nn,int fac[])
 178 {
 179   int i,n=nn;
 180
 181   for(i=0; (i<=n); i++)
 182     fac[i] = 0;
 183   fac[1] = 1;
 184   for(i=2; (i<=n); ) {
 185     if ((n % i) == 0) {
 186       fac[i]++;
 187       n = n/i;
 188     }
 189     else
 190       i++;
 191   }
 192   if (debug) {
 193     fprintf(debug,"Factorizing %d into primes:\n",nn);
 194     for(i=2; (i<=nn); i++) {
 195       if (fac[i])
 196         fprintf(debug,"%d ^ %d\n",i,fac[i]);
 197     }
 198   }
 199 }
 200
 201 static int ipow(int base,int exp)
 202 {
 203   int i,ip;
 204
 205   for(ip=1,i=0; (i<exp); i++) {
 206     ip *= base;
 207   }
 208   return ip;
 209 }
 210
 211 static int iv_comp(const void *a,const void *b)
 212 {
 213   int *ia,*ib;
 214
 215   ia = (int *)a;
 216   ib = (int *)b;
 217   if (ia[XX] != ib[XX])
 218     return (ia[XX] - ib[XX]);
 219   else if (ia[YY] != ib[YY])
 220     return (ia[YY] - ib[YY]);
 221   else
 222     return (ia[ZZ] - ib[ZZ]);
 223 }
 224
 225 static int add_bb(ivec BB[],int n,ivec b)
 226 {
 227 #define SWPX(vv,xx,yy) { int tmp; tmp=vv[xx]; vv[xx] = vv[yy]; vv[yy] = tmp; }
 228   copy_ivec(b,BB[n++]); /* x y z */
 229   SWPX(b,XX,YY);
 230   copy_ivec(b,BB[n++]); /* y x z */
 231   SWPX(b,XX,ZZ);
 232   copy_ivec(b,BB[n++]); /* z x y */
 233   SWPX(b,XX,YY);
 234   copy_ivec(b,BB[n++]); /* x z y */
 235   SWPX(b,XX,ZZ);
 236   copy_ivec(b,BB[n++]); /* y z x */
 237   SWPX(b,XX,YY);
 238   copy_ivec(b,BB[n++]); /* z y x */
 239   SWPX(b,XX,ZZ); /* Back to normal */
 240 #undef SWPX
 241   return n;
 242 }
 243
 244 static real box_weight(ivec nbox,matrix box)
 245 {
 246   rvec lx;
 247   int  m;
 248
 249   /* Calculate area of subbox */
 250   for(m=0; (m<DIM); m++)
 251     lx[m] = box[m][m]/nbox[m];
 252   return 2*(lx[XX]*lx[YY]+lx[XX]*lx[ZZ]+lx[YY]*lx[ZZ]);
 253 }
 254
 255 static int w_comp(const void *a,const void *b)
 256 {
 257   int *ia,*ib;
 258   real wa,wb;
 259
 260   ia = (int *)a;
 261   ib = (int *)b;
 262
 263   wa = box_weight(ia,BOX);
 264   wb = box_weight(ib,BOX);
 265   if (fabs(wa - wb) < 1e-4)
 266     return (iiprod(ia,ia) - iiprod(ib,ib));
 267   else if (wa < wb)
 268     return -1;
 269   else
 270     return 1;
 271 }
 272
 273 static void buildbox(int nnode,ivec nbox,matrix box)
 274 {
 275   ivec *BB,bxyz;
 276   int  i,j,m,n,n3,ny,*fx,*fy,nbb;
 277
 278   n3 = ipow(nnode,3)*6;
 279   snew(BB,n3);
 280   nbb=0;
 281   snew(fx,nnode+1);
 282   snew(fy,nnode+1);
 283   factorize(nnode,fx);
 284   for(i=0; (i<=nnode); i++) {
 285     for(m=1; (m<=fx[i]); m++) {
 286       bxyz[XX] = ipow(i,m);
 287       ny = nnode/bxyz[XX];
 288       factorize(ny,fy);
 289       for(j=0; (j<=ny); j++) {
 290         for(n=1; (n<=fy[j]); n++) {
 291           bxyz[YY] = ipow(j,n);
 292           bxyz[ZZ] = ny/bxyz[YY];
 293           if (bxyz[ZZ] > 0) {
 294             nbb = add_bb(BB,nbb,bxyz);
 295           }
 296         }
 297       }
 298     }
 299   }
 300   /* Sort boxes and remove doubles */
 301   qsort(BB,nbb,sizeof(BB[0]),iv_comp);
 302   j = 0;
 303   for(i=1; (i<nbb); i++) {
 304     if ((BB[i][XX] != BB[j][XX]) ||
 305         (BB[i][YY] != BB[j][YY]) ||
 306         (BB[i][ZZ] != BB[j][ZZ])) {
 307       j++;
 308       copy_ivec(BB[i],BB[j]);
 309     }
 310   }
 311   nbb = ++j;
 312   /* Sort boxes according to weight */
 313   copy_mat(box,BOX);
 314   qsort(BB,nbb,sizeof(BB[0]),w_comp);
 315   for(i=0; (i<nbb); i++) {
 316     fprintf(stderr,"nbox = %2d %2d %2d [ prod %3d ] area = %12.5f (nm^2)\n",
 317             BB[i][XX],BB[i][YY],BB[i][ZZ],
 318             BB[i][XX]*BB[i][YY]*BB[i][ZZ],
 319             box_weight(BB[i],box));
 320   }
 321   copy_ivec(BB[0],nbox);
 322   sfree(BB);
 323   sfree(fy);
 324   sfree(fx);
 325 }
 326
 327 void mkcompact(int astart,int nwater,int nwatom,rvec x[],rvec v[],
 328                int nnode,matrix box)
 329 {
 330   /* Make a compact configuration for each processor.
 331    * Divide the computational box in near cubic boxes and spread them
 332    * evenly over processors.
 333    */
 334 /*   ivec nbox; */
 335   int  m;
 336
 337   if (nnode <= 1)
 338     return;
 339
 340   buildbox(nnode,NBOX,box);
 341   /* copy_ivec(nbox,NBOX); */
 342   for(m=0; (m<DIM); m++)
 343     box_1[m] = 1.0/box[m][m];
 344
 345   lo_sortwater(astart,nwater,nwatom,x,v,TRUE);
 346 }