2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
45 #include "sortwater.h"
47 #include "gromacs/random/random.h"
48 #include "gromacs/utility/smalloc.h"
50 static rvec *xptr, box_1;
55 void randwater(int astart, int nwater, int nwatom, rvec x[], rvec v[],
58 int i, j, wi, wj, *tab;
62 for (i = 0; (i < nwater); i++)
66 for (j = 0; (j < 23*nwater); j++)
68 wi = (int) (nwater*gmx_rng_uniform_real(rng)) % nwater;
71 wj = (int) (nwater*gmx_rng_uniform_real(rng)) % nwater;
74 wi = astart+wi*nwatom;
75 wj = astart+wj*nwatom;
76 /* Swap coords for wi and wj */
77 for (i = 0; (i < nwatom); i++)
79 copy_rvec(x[wi+i], buf);
80 copy_rvec(x[wj+i], x[wi+i]);
81 copy_rvec(buf, x[wj+i]);
84 copy_rvec(v[wi+i], buf);
85 copy_rvec(v[wj+i], v[wi+i]);
86 copy_rvec(buf, v[wj+i]);
93 static int rvcomp(const void *a, const void *b)
97 aa = nwat*(*((int *)a));
98 bb = nwat*(*((int *)b));
99 if (xptr[aa][XX] < xptr[bb][XX])
103 else if (xptr[aa][XX] > xptr[bb][XX])
113 static int block_index(rvec x, ivec nbox)
118 for (m = 0; (m < DIM); m++)
120 ixyz[m] = ((int)((1+x[m]*box_1[m])*nbox[m])) % nbox[m];
122 return ixyz[XX]*nbox[YY]*nbox[ZZ]+ixyz[YY]*nbox[ZZ]+ixyz[ZZ];
125 static int blockcomp(const void *a, const void *b)
127 int aa, bb, aind, bind;
129 aa = nwat*(*((int *)a));
130 bb = nwat*(*((int *)b));
132 aind = block_index(xptr[aa], NBOX);
133 bind = block_index(xptr[bb], NBOX);
137 if (xptr[aa][XX] < xptr[bb][XX])
141 else if (xptr[aa][XX] > xptr[bb][XX])
156 static void lo_sortwater(int astart, int nwater, int nwatom, rvec x[], rvec v[],
163 /* Sort indices to rvecs */
164 snew(rvindex, nwater);
165 for (i = 0; (i < nwater); i++)
172 qsort(rvindex, nwater, sizeof(rvindex[0]), bBlock ? blockcomp : rvcomp);
175 for (i = 0; (i < nwater); i++)
177 rvi = rvindex[i]*nwatom;
178 fprintf(debug, "rvindex[%5d] = %5d (x = %8.3f %8.3f %8.3f)\n",
179 i, rvi, x[astart+rvi][XX], x[astart+rvi][YY], x[astart+rvi][ZZ]);
182 snew(tmp, nwater*nwatom);
184 for (i = 0; (i < nwater); i++)
186 i0 = astart+nwatom*rvindex[i];
187 for (j = 0; (j < nwatom); j++)
189 copy_rvec(x[i0+j], tmp[nwatom*i+j]);
192 for (i = 0; (i < nwater*nwatom); i++)
194 copy_rvec(tmp[i], x[astart+i]);
197 for (i = 0; (i < nwater); i++)
199 i0 = astart+nwatom*rvindex[i];
200 for (j = 0; (j < nwatom); j++)
202 copy_rvec(v[i0+j], tmp[nwatom*i+j]);
205 for (i = 0; (i < nwater*nwatom); i++)
207 copy_rvec(tmp[i], v[astart+i]);
214 void sortwater(int astart, int nwater, int nwatom, rvec x[], rvec v[])
216 lo_sortwater(astart, nwater, nwatom, x, v, FALSE);
219 static void factorize(int nn, int fac[])
223 for (i = 0; (i <= n); i++)
228 for (i = 2; (i <= n); )
242 fprintf(debug, "Factorizing %d into primes:\n", nn);
243 for (i = 2; (i <= nn); i++)
247 fprintf(debug, "%d ^ %d\n", i, fac[i]);
253 static int ipow(int base, int exp)
257 for (ip = 1, i = 0; (i < exp); i++)
264 static int iv_comp(const void *a, const void *b)
270 if (ia[XX] != ib[XX])
272 return (ia[XX] - ib[XX]);
274 else if (ia[YY] != ib[YY])
276 return (ia[YY] - ib[YY]);
280 return (ia[ZZ] - ib[ZZ]);
284 static int add_bb(ivec BB[], int n, ivec b)
286 #define SWPX(vv, xx, yy) { int tmp; tmp = vv[xx]; vv[xx] = vv[yy]; vv[yy] = tmp; }
287 copy_ivec(b, BB[n++]); /* x y z */
289 copy_ivec(b, BB[n++]); /* y x z */
291 copy_ivec(b, BB[n++]); /* z x y */
293 copy_ivec(b, BB[n++]); /* x z y */
295 copy_ivec(b, BB[n++]); /* y z x */
297 copy_ivec(b, BB[n++]); /* z y x */
298 SWPX(b, XX, ZZ); /* Back to normal */
303 static real box_weight(ivec nbox, matrix box)
308 /* Calculate area of subbox */
309 for (m = 0; (m < DIM); m++)
311 lx[m] = box[m][m]/nbox[m];
313 return 2*(lx[XX]*lx[YY]+lx[XX]*lx[ZZ]+lx[YY]*lx[ZZ]);
316 static int w_comp(const void *a, const void *b)
324 wa = box_weight(ia, BOX);
325 wb = box_weight(ib, BOX);
326 if (fabs(wa - wb) < 1e-4)
328 return (iiprod(ia, ia) - iiprod(ib, ib));
340 static void buildbox(int nnode, ivec nbox, matrix box)
343 int i, j, m, n, n3, ny, *fx, *fy, nbb;
345 n3 = ipow(nnode, 3)*6;
350 factorize(nnode, fx);
351 for (i = 0; (i <= nnode); i++)
353 for (m = 1; (m <= fx[i]); m++)
355 bxyz[XX] = ipow(i, m);
358 for (j = 0; (j <= ny); j++)
360 for (n = 1; (n <= fy[j]); n++)
362 bxyz[YY] = ipow(j, n);
363 bxyz[ZZ] = ny/bxyz[YY];
366 nbb = add_bb(BB, nbb, bxyz);
372 /* Sort boxes and remove doubles */
373 qsort(BB, nbb, sizeof(BB[0]), iv_comp);
375 for (i = 1; (i < nbb); i++)
377 if ((BB[i][XX] != BB[j][XX]) ||
378 (BB[i][YY] != BB[j][YY]) ||
379 (BB[i][ZZ] != BB[j][ZZ]))
382 copy_ivec(BB[i], BB[j]);
386 /* Sort boxes according to weight */
388 qsort(BB, nbb, sizeof(BB[0]), w_comp);
389 for (i = 0; (i < nbb); i++)
391 fprintf(stderr, "nbox = %2d %2d %2d [ prod %3d ] area = %12.5f (nm^2)\n",
392 BB[i][XX], BB[i][YY], BB[i][ZZ],
393 BB[i][XX]*BB[i][YY]*BB[i][ZZ],
394 box_weight(BB[i], box));
396 copy_ivec(BB[0], nbox);
402 void mkcompact(int astart, int nwater, int nwatom, rvec x[], rvec v[],
403 int nnode, matrix box)
405 /* Make a compact configuration for each processor.
406 * Divide the computational box in near cubic boxes and spread them
407 * evenly over processors.
417 buildbox(nnode, NBOX, box);
418 /* copy_ivec(nbox,NBOX); */
419 for (m = 0; (m < DIM); m++)
421 box_1[m] = 1.0/box[m][m];
424 lo_sortwater(astart, nwater, nwatom, x, v, TRUE);