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.
41 #include "gromacs/legacyheaders/typedefs.h"
42 #include "gromacs/math/vec.h"
43 #include "sortwater.h"
45 #include "gromacs/random/random.h"
46 #include "gromacs/utility/smalloc.h"
48 static rvec *xptr, box_1;
53 void randwater(int astart, int nwater, int nwatom, rvec x[], rvec v[],
56 int i, j, wi, wj, *tab;
60 for (i = 0; (i < nwater); i++)
64 for (j = 0; (j < 23*nwater); j++)
66 wi = (int) (nwater*gmx_rng_uniform_real(rng)) % nwater;
69 wj = (int) (nwater*gmx_rng_uniform_real(rng)) % nwater;
72 wi = astart+wi*nwatom;
73 wj = astart+wj*nwatom;
74 /* Swap coords for wi and wj */
75 for (i = 0; (i < nwatom); i++)
77 copy_rvec(x[wi+i], buf);
78 copy_rvec(x[wj+i], x[wi+i]);
79 copy_rvec(buf, x[wj+i]);
82 copy_rvec(v[wi+i], buf);
83 copy_rvec(v[wj+i], v[wi+i]);
84 copy_rvec(buf, v[wj+i]);
91 static int rvcomp(const void *a, const void *b)
95 aa = nwat*(*((int *)a));
96 bb = nwat*(*((int *)b));
97 if (xptr[aa][XX] < xptr[bb][XX])
101 else if (xptr[aa][XX] > xptr[bb][XX])
111 static int block_index(rvec x, ivec nbox)
116 for (m = 0; (m < DIM); m++)
118 ixyz[m] = ((int)((1+x[m]*box_1[m])*nbox[m])) % nbox[m];
120 return ixyz[XX]*nbox[YY]*nbox[ZZ]+ixyz[YY]*nbox[ZZ]+ixyz[ZZ];
123 static int blockcomp(const void *a, const void *b)
125 int aa, bb, aind, bind;
127 aa = nwat*(*((int *)a));
128 bb = nwat*(*((int *)b));
130 aind = block_index(xptr[aa], NBOX);
131 bind = block_index(xptr[bb], NBOX);
135 if (xptr[aa][XX] < xptr[bb][XX])
139 else if (xptr[aa][XX] > xptr[bb][XX])
154 static void lo_sortwater(int astart, int nwater, int nwatom, rvec x[], rvec v[],
161 /* Sort indices to rvecs */
162 snew(rvindex, nwater);
163 for (i = 0; (i < nwater); i++)
170 qsort(rvindex, nwater, sizeof(rvindex[0]), bBlock ? blockcomp : rvcomp);
173 for (i = 0; (i < nwater); i++)
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]);
180 snew(tmp, nwater*nwatom);
182 for (i = 0; (i < nwater); i++)
184 i0 = astart+nwatom*rvindex[i];
185 for (j = 0; (j < nwatom); j++)
187 copy_rvec(x[i0+j], tmp[nwatom*i+j]);
190 for (i = 0; (i < nwater*nwatom); i++)
192 copy_rvec(tmp[i], x[astart+i]);
195 for (i = 0; (i < nwater); i++)
197 i0 = astart+nwatom*rvindex[i];
198 for (j = 0; (j < nwatom); j++)
200 copy_rvec(v[i0+j], tmp[nwatom*i+j]);
203 for (i = 0; (i < nwater*nwatom); i++)
205 copy_rvec(tmp[i], v[astart+i]);
212 void sortwater(int astart, int nwater, int nwatom, rvec x[], rvec v[])
214 lo_sortwater(astart, nwater, nwatom, x, v, FALSE);
217 static void factorize(int nn, int fac[])
221 for (i = 0; (i <= n); i++)
226 for (i = 2; (i <= n); )
240 fprintf(debug, "Factorizing %d into primes:\n", nn);
241 for (i = 2; (i <= nn); i++)
245 fprintf(debug, "%d ^ %d\n", i, fac[i]);
251 static int ipow(int base, int exp)
255 for (ip = 1, i = 0; (i < exp); i++)
262 static int iv_comp(const void *a, const void *b)
268 if (ia[XX] != ib[XX])
270 return (ia[XX] - ib[XX]);
272 else if (ia[YY] != ib[YY])
274 return (ia[YY] - ib[YY]);
278 return (ia[ZZ] - ib[ZZ]);
282 static int add_bb(ivec BB[], int n, ivec b)
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 */
287 copy_ivec(b, BB[n++]); /* y x z */
289 copy_ivec(b, BB[n++]); /* z x y */
291 copy_ivec(b, BB[n++]); /* x z y */
293 copy_ivec(b, BB[n++]); /* y z x */
295 copy_ivec(b, BB[n++]); /* z y x */
296 SWPX(b, XX, ZZ); /* Back to normal */
301 static real box_weight(ivec nbox, matrix box)
306 /* Calculate area of subbox */
307 for (m = 0; (m < DIM); m++)
309 lx[m] = box[m][m]/nbox[m];
311 return 2*(lx[XX]*lx[YY]+lx[XX]*lx[ZZ]+lx[YY]*lx[ZZ]);
314 static int w_comp(const void *a, const void *b)
322 wa = box_weight(ia, BOX);
323 wb = box_weight(ib, BOX);
324 if (fabs(wa - wb) < 1e-4)
326 return (iiprod(ia, ia) - iiprod(ib, ib));
338 static void buildbox(int nnode, ivec nbox, matrix box)
341 int i, j, m, n, n3, ny, *fx, *fy, nbb;
343 n3 = ipow(nnode, 3)*6;
348 factorize(nnode, fx);
349 for (i = 0; (i <= nnode); i++)
351 for (m = 1; (m <= fx[i]); m++)
353 bxyz[XX] = ipow(i, m);
356 for (j = 0; (j <= ny); j++)
358 for (n = 1; (n <= fy[j]); n++)
360 bxyz[YY] = ipow(j, n);
361 bxyz[ZZ] = ny/bxyz[YY];
364 nbb = add_bb(BB, nbb, bxyz);
370 /* Sort boxes and remove doubles */
371 qsort(BB, nbb, sizeof(BB[0]), iv_comp);
373 for (i = 1; (i < nbb); i++)
375 if ((BB[i][XX] != BB[j][XX]) ||
376 (BB[i][YY] != BB[j][YY]) ||
377 (BB[i][ZZ] != BB[j][ZZ]))
380 copy_ivec(BB[i], BB[j]);
384 /* Sort boxes according to weight */
386 qsort(BB, nbb, sizeof(BB[0]), w_comp);
387 for (i = 0; (i < nbb); i++)
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));
394 copy_ivec(BB[0], nbox);
400 void mkcompact(int astart, int nwater, int nwatom, rvec x[], rvec v[],
401 int nnode, matrix box)
403 /* Make a compact configuration for each processor.
404 * Divide the computational box in near cubic boxes and spread them
405 * evenly over processors.
415 buildbox(nnode, NBOX, box);
416 /* copy_ivec(nbox,NBOX); */
417 for (m = 0; (m < DIM); m++)
419 box_1[m] = 1.0/box[m][m];
422 lo_sortwater(astart, nwater, nwatom, x, v, TRUE);