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.
42 #include "gromacs/random/random.h"
43 #include "gromacs/utility/smalloc.h"
45 #include "sortwater.h"
47 static rvec *xptr, box_1;
52 void randwater(int astart, int nwater, int nwatom, rvec x[], rvec v[],
55 int i, j, wi, wj, *tab;
59 for (i = 0; (i < nwater); i++)
63 for (j = 0; (j < 23*nwater); j++)
65 wi = (int) (nwater*gmx_rng_uniform_real(rng)) % nwater;
68 wj = (int) (nwater*gmx_rng_uniform_real(rng)) % nwater;
71 wi = astart+wi*nwatom;
72 wj = astart+wj*nwatom;
73 /* Swap coords for wi and wj */
74 for (i = 0; (i < nwatom); i++)
76 copy_rvec(x[wi+i], buf);
77 copy_rvec(x[wj+i], x[wi+i]);
78 copy_rvec(buf, x[wj+i]);
81 copy_rvec(v[wi+i], buf);
82 copy_rvec(v[wj+i], v[wi+i]);
83 copy_rvec(buf, v[wj+i]);
90 static int rvcomp(const void *a, const void *b)
94 aa = nwat*(*((int *)a));
95 bb = nwat*(*((int *)b));
96 if (xptr[aa][XX] < xptr[bb][XX])
100 else if (xptr[aa][XX] > xptr[bb][XX])
110 static int block_index(rvec x, ivec nbox)
115 for (m = 0; (m < DIM); m++)
117 ixyz[m] = ((int)((1+x[m]*box_1[m])*nbox[m])) % nbox[m];
119 return ixyz[XX]*nbox[YY]*nbox[ZZ]+ixyz[YY]*nbox[ZZ]+ixyz[ZZ];
122 static int blockcomp(const void *a, const void *b)
124 int aa, bb, aind, bind;
126 aa = nwat*(*((int *)a));
127 bb = nwat*(*((int *)b));
129 aind = block_index(xptr[aa], NBOX);
130 bind = block_index(xptr[bb], NBOX);
134 if (xptr[aa][XX] < xptr[bb][XX])
138 else if (xptr[aa][XX] > xptr[bb][XX])
153 static void lo_sortwater(int astart, int nwater, int nwatom, rvec x[], rvec v[],
160 /* Sort indices to rvecs */
161 snew(rvindex, nwater);
162 for (i = 0; (i < nwater); i++)
169 qsort(rvindex, nwater, sizeof(rvindex[0]), bBlock ? blockcomp : rvcomp);
172 for (i = 0; (i < nwater); i++)
174 rvi = rvindex[i]*nwatom;
175 fprintf(debug, "rvindex[%5d] = %5d (x = %8.3f %8.3f %8.3f)\n",
176 i, rvi, x[astart+rvi][XX], x[astart+rvi][YY], x[astart+rvi][ZZ]);
179 snew(tmp, nwater*nwatom);
181 for (i = 0; (i < nwater); i++)
183 i0 = astart+nwatom*rvindex[i];
184 for (j = 0; (j < nwatom); j++)
186 copy_rvec(x[i0+j], tmp[nwatom*i+j]);
189 for (i = 0; (i < nwater*nwatom); i++)
191 copy_rvec(tmp[i], x[astart+i]);
194 for (i = 0; (i < nwater); i++)
196 i0 = astart+nwatom*rvindex[i];
197 for (j = 0; (j < nwatom); j++)
199 copy_rvec(v[i0+j], tmp[nwatom*i+j]);
202 for (i = 0; (i < nwater*nwatom); i++)
204 copy_rvec(tmp[i], v[astart+i]);
211 void sortwater(int astart, int nwater, int nwatom, rvec x[], rvec v[])
213 lo_sortwater(astart, nwater, nwatom, x, v, FALSE);
216 static void factorize(int nn, int fac[])
220 for (i = 0; (i <= n); i++)
225 for (i = 2; (i <= n); )
239 fprintf(debug, "Factorizing %d into primes:\n", nn);
240 for (i = 2; (i <= nn); i++)
244 fprintf(debug, "%d ^ %d\n", i, fac[i]);
250 static int ipow(int base, int exp)
254 for (ip = 1, i = 0; (i < exp); i++)
261 static int iv_comp(const void *a, const void *b)
267 if (ia[XX] != ib[XX])
269 return (ia[XX] - ib[XX]);
271 else if (ia[YY] != ib[YY])
273 return (ia[YY] - ib[YY]);
277 return (ia[ZZ] - ib[ZZ]);
281 static int add_bb(ivec BB[], int n, ivec b)
283 #define SWPX(vv, xx, yy) { int tmp; tmp = vv[xx]; vv[xx] = vv[yy]; vv[yy] = tmp; }
284 copy_ivec(b, BB[n++]); /* x y z */
286 copy_ivec(b, BB[n++]); /* y x z */
288 copy_ivec(b, BB[n++]); /* z x y */
290 copy_ivec(b, BB[n++]); /* x z y */
292 copy_ivec(b, BB[n++]); /* y z x */
294 copy_ivec(b, BB[n++]); /* z y x */
295 SWPX(b, XX, ZZ); /* Back to normal */
300 static real box_weight(ivec nbox, matrix box)
305 /* Calculate area of subbox */
306 for (m = 0; (m < DIM); m++)
308 lx[m] = box[m][m]/nbox[m];
310 return 2*(lx[XX]*lx[YY]+lx[XX]*lx[ZZ]+lx[YY]*lx[ZZ]);
313 static int w_comp(const void *a, const void *b)
321 wa = box_weight(ia, BOX);
322 wb = box_weight(ib, BOX);
323 if (fabs(wa - wb) < 1e-4)
325 return (iiprod(ia, ia) - iiprod(ib, ib));
337 static void buildbox(int nnode, ivec nbox, matrix box)
340 int i, j, m, n, n3, ny, *fx, *fy, nbb;
342 n3 = ipow(nnode, 3)*6;
347 factorize(nnode, fx);
348 for (i = 0; (i <= nnode); i++)
350 for (m = 1; (m <= fx[i]); m++)
352 bxyz[XX] = ipow(i, m);
355 for (j = 0; (j <= ny); j++)
357 for (n = 1; (n <= fy[j]); n++)
359 bxyz[YY] = ipow(j, n);
360 bxyz[ZZ] = ny/bxyz[YY];
363 nbb = add_bb(BB, nbb, bxyz);
369 /* Sort boxes and remove doubles */
370 qsort(BB, nbb, sizeof(BB[0]), iv_comp);
372 for (i = 1; (i < nbb); i++)
374 if ((BB[i][XX] != BB[j][XX]) ||
375 (BB[i][YY] != BB[j][YY]) ||
376 (BB[i][ZZ] != BB[j][ZZ]))
379 copy_ivec(BB[i], BB[j]);
383 /* Sort boxes according to weight */
385 qsort(BB, nbb, sizeof(BB[0]), w_comp);
386 for (i = 0; (i < nbb); i++)
388 fprintf(stderr, "nbox = %2d %2d %2d [ prod %3d ] area = %12.5f (nm^2)\n",
389 BB[i][XX], BB[i][YY], BB[i][ZZ],
390 BB[i][XX]*BB[i][YY]*BB[i][ZZ],
391 box_weight(BB[i], box));
393 copy_ivec(BB[0], nbox);
399 void mkcompact(int astart, int nwater, int nwatom, rvec x[], rvec v[],
400 int nnode, matrix box)
402 /* Make a compact configuration for each processor.
403 * Divide the computational box in near cubic boxes and spread them
404 * evenly over processors.
414 buildbox(nnode, NBOX, box);
415 /* copy_ivec(nbox,NBOX); */
416 for (m = 0; (m < DIM); m++)
418 box_1[m] = 1.0/box[m][m];
421 lo_sortwater(astart, nwater, nwatom, x, v, TRUE);