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