src/gromacs/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  * Copyright (c) 2010, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #ifdef HAVE_CONFIG_H
  38 #include <config.h>
  39 #endif
  40
  41 #include "typedefs.h"
  42 #include "random.h"
  43 #include "smalloc.h"
  44 #include "vec.h"
  45 #include "sortwater.h"
  46
  47 static rvec   *xptr, box_1;
  48 static int     nwat;
  49 static matrix  BOX;
  50 static ivec    NBOX;
  51
  52 void randwater(int astart, int nwater, int nwatom, rvec x[], rvec v[], int *seed)
  53 {
  54     int  i, j, wi, wj, *tab;
  55     rvec buf;
  56
  57     snew(tab, nwater);
  58     for (i = 0; (i < nwater); i++)
  59     {
  60         tab[i] = i;
  61     }
  62     for (j = 0; (j < 23*nwater); j++)
  63     {
  64         wi = (int) (nwater*rando(seed)) % nwater;
  65         do
  66         {
  67             wj = (int) (nwater*rando(seed)) % nwater;
  68         }
  69         while (wi == wj);
  70         wi = astart+wi*nwatom;
  71         wj = astart+wj*nwatom;
  72         /* Swap coords for wi and wj */
  73         for (i = 0; (i < nwatom); i++)
  74         {
  75             copy_rvec(x[wi+i], buf);
  76             copy_rvec(x[wj+i], x[wi+i]);
  77             copy_rvec(buf, x[wj+i]);
  78             if (v)
  79             {
  80                 copy_rvec(v[wi+i], buf);
  81                 copy_rvec(v[wj+i], v[wi+i]);
  82                 copy_rvec(buf, v[wj+i]);
  83             }
  84         }
  85     }
  86     sfree(tab);
  87 }
  88
  89 static int rvcomp(const void *a, const void *b)
  90 {
  91     int aa, bb;
  92
  93     aa = nwat*(*((int *)a));
  94     bb = nwat*(*((int *)b));
  95     if (xptr[aa][XX] < xptr[bb][XX])
  96     {
  97         return -1;
  98     }
  99     else if (xptr[aa][XX] > xptr[bb][XX])
 100     {
 101         return 1;
 102     }
 103     else
 104     {
 105         return 0;
 106     }
 107 }
 108
 109 static int block_index(rvec x, ivec nbox)
 110 {
 111     ivec ixyz;
 112     int  m;
 113
 114     for (m = 0; (m < DIM); m++)
 115     {
 116         ixyz[m] = ((int)((1+x[m]*box_1[m])*nbox[m])) % nbox[m];
 117     }
 118     return ixyz[XX]*nbox[YY]*nbox[ZZ]+ixyz[YY]*nbox[ZZ]+ixyz[ZZ];
 119 }
 120
 121 static int blockcomp(const void *a, const void *b)
 122 {
 123     int aa, bb, aind, bind;
 124
 125     aa   = nwat*(*((int *)a));
 126     bb   = nwat*(*((int *)b));
 127
 128     aind = block_index(xptr[aa], NBOX);
 129     bind = block_index(xptr[bb], NBOX);
 130
 131     if (aind == bind)
 132     {
 133         if (xptr[aa][XX] < xptr[bb][XX])
 134         {
 135             return -1;
 136         }
 137         else if (xptr[aa][XX] > xptr[bb][XX])
 138         {
 139             return 1;
 140         }
 141         else
 142         {
 143             return 0;
 144         }
 145     }
 146     else
 147     {
 148         return aind-bind;
 149     }
 150 }
 151
 152 static void lo_sortwater(int astart, int nwater, int nwatom, rvec x[], rvec v[],
 153                          gmx_bool bBlock)
 154 {
 155     int   i, j, i0, rvi;
 156     int  *rvindex;
 157     rvec *tmp;
 158
 159     /* Sort indices to rvecs */
 160     snew(rvindex, nwater);
 161     for (i = 0; (i < nwater); i++)
 162     {
 163         rvindex[i] = i;
 164     }
 165     xptr = x+astart;
 166     nwat = nwatom;
 167
 168     qsort(rvindex, nwater, sizeof(rvindex[0]), bBlock ? blockcomp : rvcomp);
 169     if (debug)
 170     {
 171         for (i = 0; (i < nwater); i++)
 172         {
 173             rvi = rvindex[i]*nwatom;
 174             fprintf(debug, "rvindex[%5d] = %5d (x = %8.3f  %8.3f  %8.3f)\n",
 175                     i, rvi, x[astart+rvi][XX], x[astart+rvi][YY], x[astart+rvi][ZZ]);
 176         }
 177     }
 178     snew(tmp, nwater*nwatom);
 179
 180     for (i = 0; (i < nwater); i++)
 181     {
 182         i0 = astart+nwatom*rvindex[i];
 183         for (j = 0; (j < nwatom); j++)
 184         {
 185             copy_rvec(x[i0+j], tmp[nwatom*i+j]);
 186         }
 187     }
 188     for (i = 0; (i < nwater*nwatom); i++)
 189     {
 190         copy_rvec(tmp[i], x[astart+i]);
 191     }
 192
 193     for (i = 0; (i < nwater); i++)
 194     {
 195         i0 = astart+nwatom*rvindex[i];
 196         for (j = 0; (j < nwatom); j++)
 197         {
 198             copy_rvec(v[i0+j], tmp[nwatom*i+j]);
 199         }
 200     }
 201     for (i = 0; (i < nwater*nwatom); i++)
 202     {
 203         copy_rvec(tmp[i], v[astart+i]);
 204     }
 205
 206     sfree(tmp);
 207     sfree(rvindex);
 208 }
 209
 210 void sortwater(int astart, int nwater, int nwatom, rvec x[], rvec v[])
 211 {
 212     lo_sortwater(astart, nwater, nwatom, x, v, FALSE);
 213 }
 214
 215 static void factorize(int nn, int fac[])
 216 {
 217     int i, n = nn;
 218
 219     for (i = 0; (i <= n); i++)
 220     {
 221         fac[i] = 0;
 222     }
 223     fac[1] = 1;
 224     for (i = 2; (i <= n); )
 225     {
 226         if ((n % i) == 0)
 227         {
 228             fac[i]++;
 229             n = n/i;
 230         }
 231         else
 232         {
 233             i++;
 234         }
 235     }
 236     if (debug)
 237     {
 238         fprintf(debug, "Factorizing %d into primes:\n", nn);
 239         for (i = 2; (i <= nn); i++)
 240         {
 241             if (fac[i])
 242             {
 243                 fprintf(debug, "%d ^ %d\n", i, fac[i]);
 244             }
 245         }
 246     }
 247 }
 248
 249 static int ipow(int base, int exp)
 250 {
 251     int i, ip;
 252
 253     for (ip = 1, i = 0; (i < exp); i++)
 254     {
 255         ip *= base;
 256     }
 257     return ip;
 258 }
 259
 260 static int iv_comp(const void *a, const void *b)
 261 {
 262     int *ia, *ib;
 263
 264     ia = (int *)a;
 265     ib = (int *)b;
 266     if (ia[XX] != ib[XX])
 267     {
 268         return (ia[XX] - ib[XX]);
 269     }
 270     else if (ia[YY] != ib[YY])
 271     {
 272         return (ia[YY] - ib[YY]);
 273     }
 274     else
 275     {
 276         return (ia[ZZ] - ib[ZZ]);
 277     }
 278 }
 279
 280 static int add_bb(ivec BB[], int n, ivec b)
 281 {
 282 #define SWPX(vv, xx, yy) { int tmp; tmp = vv[xx]; vv[xx] = vv[yy]; vv[yy] = tmp; }
 283     copy_ivec(b, BB[n++]); /* x y z */
 284     SWPX(b, XX, YY);
 285     copy_ivec(b, BB[n++]); /* y x z */
 286     SWPX(b, XX, ZZ);
 287     copy_ivec(b, BB[n++]); /* z x y */
 288     SWPX(b, XX, YY);
 289     copy_ivec(b, BB[n++]); /* x z y */
 290     SWPX(b, XX, ZZ);
 291     copy_ivec(b, BB[n++]); /* y z x */
 292     SWPX(b, XX, YY);
 293     copy_ivec(b, BB[n++]); /* z y x */
 294     SWPX(b, XX, ZZ);       /* Back to normal */
 295 #undef SWPX
 296     return n;
 297 }
 298
 299 static real box_weight(ivec nbox, matrix box)
 300 {
 301     rvec lx;
 302     int  m;
 303
 304     /* Calculate area of subbox */
 305     for (m = 0; (m < DIM); m++)
 306     {
 307         lx[m] = box[m][m]/nbox[m];
 308     }
 309     return 2*(lx[XX]*lx[YY]+lx[XX]*lx[ZZ]+lx[YY]*lx[ZZ]);
 310 }
 311
 312 static int w_comp(const void *a, const void *b)
 313 {
 314     int *ia, *ib;
 315     real wa, wb;
 316
 317     ia = (int *)a;
 318     ib = (int *)b;
 319
 320     wa = box_weight(ia, BOX);
 321     wb = box_weight(ib, BOX);
 322     if (fabs(wa - wb) < 1e-4)
 323     {
 324         return (iiprod(ia, ia) - iiprod(ib, ib));
 325     }
 326     else if (wa < wb)
 327     {
 328         return -1;
 329     }
 330     else
 331     {
 332         return 1;
 333     }
 334 }
 335
 336 static void buildbox(int nnode, ivec nbox, matrix box)
 337 {
 338     ivec *BB, bxyz;
 339     int   i, j, m, n, n3, ny, *fx, *fy, nbb;
 340
 341     n3 = ipow(nnode, 3)*6;
 342     snew(BB, n3);
 343     nbb = 0;
 344     snew(fx, nnode+1);
 345     snew(fy, nnode+1);
 346     factorize(nnode, fx);
 347     for (i = 0; (i <= nnode); i++)
 348     {
 349         for (m = 1; (m <= fx[i]); m++)
 350         {
 351             bxyz[XX] = ipow(i, m);
 352             ny       = nnode/bxyz[XX];
 353             factorize(ny, fy);
 354             for (j = 0; (j <= ny); j++)
 355             {
 356                 for (n = 1; (n <= fy[j]); n++)
 357                 {
 358                     bxyz[YY] = ipow(j, n);
 359                     bxyz[ZZ] = ny/bxyz[YY];
 360                     if (bxyz[ZZ] > 0)
 361                     {
 362                         nbb = add_bb(BB, nbb, bxyz);
 363                     }
 364                 }
 365             }
 366         }
 367     }
 368     /* Sort boxes and remove doubles */
 369     qsort(BB, nbb, sizeof(BB[0]), iv_comp);
 370     j = 0;
 371     for (i = 1; (i < nbb); i++)
 372     {
 373         if ((BB[i][XX] != BB[j][XX]) ||
 374             (BB[i][YY] != BB[j][YY]) ||
 375             (BB[i][ZZ] != BB[j][ZZ]))
 376         {
 377             j++;
 378             copy_ivec(BB[i], BB[j]);
 379         }
 380     }
 381     nbb = ++j;
 382     /* Sort boxes according to weight */
 383     copy_mat(box, BOX);
 384     qsort(BB, nbb, sizeof(BB[0]), w_comp);
 385     for (i = 0; (i < nbb); i++)
 386     {
 387         fprintf(stderr, "nbox = %2d %2d %2d [ prod %3d ] area = %12.5f (nm^2)\n",
 388                 BB[i][XX], BB[i][YY], BB[i][ZZ],
 389                 BB[i][XX]*BB[i][YY]*BB[i][ZZ],
 390                 box_weight(BB[i], box));
 391     }
 392     copy_ivec(BB[0], nbox);
 393     sfree(BB);
 394     sfree(fy);
 395     sfree(fx);
 396 }
 397
 398 void mkcompact(int astart, int nwater, int nwatom, rvec x[], rvec v[],
 399                int nnode, matrix box)
 400 {
 401     /* Make a compact configuration for each processor.
 402      * Divide the computational box in near cubic boxes and spread them
 403      * evenly over processors.
 404      */
 405 /*   ivec nbox; */
 406     int  m;
 407
 408     if (nnode <= 1)
 409     {
 410         return;
 411     }
 412
 413     buildbox(nnode, NBOX, box);
 414     /* copy_ivec(nbox,NBOX); */
 415     for (m = 0; (m < DIM); m++)
 416     {
 417         box_1[m] = 1.0/box[m][m];
 418     }
 419
 420     lo_sortwater(astart, nwater, nwatom, x, v, TRUE);
 421 }