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 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
46 #include "sortwater.h"
48 static rvec *xptr,box_1;
53 void randwater(int astart,int nwater,int nwatom,rvec x[],rvec v[],int *seed)
59 for(i=0; (i<nwater); i++)
61 for(j=0; (j<23*nwater); j++) {
62 wi = (int) (nwater*rando(seed)) % nwater;
64 wj = (int) (nwater*rando(seed)) % nwater;
66 wi = astart+wi*nwatom;
67 wj = astart+wj*nwatom;
68 /* Swap coords for wi and wj */
69 for(i=0; (i<nwatom); i++) {
70 copy_rvec(x[wi+i],buf);
71 copy_rvec(x[wj+i],x[wi+i]);
72 copy_rvec(buf,x[wj+i]);
74 copy_rvec(v[wi+i],buf);
75 copy_rvec(v[wj+i],v[wi+i]);
76 copy_rvec(buf,v[wj+i]);
83 static int rvcomp(const void *a,const void *b)
87 aa = nwat*(*((int *)a));
88 bb = nwat*(*((int *)b));
89 if (xptr[aa][XX] < xptr[bb][XX])
91 else if (xptr[aa][XX] > xptr[bb][XX])
97 static int block_index(rvec x,ivec nbox)
102 for(m=0; (m<DIM); m++)
103 ixyz[m] = ((int)((1+x[m]*box_1[m])*nbox[m])) % nbox[m];
104 return ixyz[XX]*nbox[YY]*nbox[ZZ]+ixyz[YY]*nbox[ZZ]+ixyz[ZZ];
107 static int blockcomp(const void *a,const void *b)
111 aa = nwat*(*((int *)a));
112 bb = nwat*(*((int *)b));
114 aind = block_index(xptr[aa],NBOX);
115 bind = block_index(xptr[bb],NBOX);
118 if (xptr[aa][XX] < xptr[bb][XX])
120 else if (xptr[aa][XX] > xptr[bb][XX])
129 static void lo_sortwater(int astart,int nwater,int nwatom,rvec x[],rvec v[],
136 /* Sort indices to rvecs */
137 snew(rvindex,nwater);
138 for(i=0; (i<nwater); i++)
143 qsort(rvindex,nwater,sizeof(rvindex[0]),bBlock ? blockcomp : rvcomp);
145 for(i=0; (i<nwater); i++) {
146 rvi = rvindex[i]*nwatom;
147 fprintf(debug,"rvindex[%5d] = %5d (x = %8.3f %8.3f %8.3f)\n",
148 i,rvi,x[astart+rvi][XX],x[astart+rvi][YY],x[astart+rvi][ZZ]);
150 snew(tmp,nwater*nwatom);
152 for(i=0; (i<nwater); i++) {
153 i0 = astart+nwatom*rvindex[i];
154 for(j=0; (j<nwatom); j++)
155 copy_rvec(x[i0+j],tmp[nwatom*i+j]);
157 for(i=0; (i<nwater*nwatom); i++)
158 copy_rvec(tmp[i],x[astart+i]);
160 for(i=0; (i<nwater); i++) {
161 i0 = astart+nwatom*rvindex[i];
162 for(j=0; (j<nwatom); j++)
163 copy_rvec(v[i0+j],tmp[nwatom*i+j]);
165 for(i=0; (i<nwater*nwatom); i++)
166 copy_rvec(tmp[i],v[astart+i]);
172 void sortwater(int astart,int nwater,int nwatom,rvec x[],rvec v[])
174 lo_sortwater(astart,nwater,nwatom,x,v,FALSE);
177 static void factorize(int nn,int fac[])
181 for(i=0; (i<=n); i++)
193 fprintf(debug,"Factorizing %d into primes:\n",nn);
194 for(i=2; (i<=nn); i++) {
196 fprintf(debug,"%d ^ %d\n",i,fac[i]);
201 static int ipow(int base,int exp)
205 for(ip=1,i=0; (i<exp); i++) {
211 static int iv_comp(const void *a,const void *b)
217 if (ia[XX] != ib[XX])
218 return (ia[XX] - ib[XX]);
219 else if (ia[YY] != ib[YY])
220 return (ia[YY] - ib[YY]);
222 return (ia[ZZ] - ib[ZZ]);
225 static int add_bb(ivec BB[],int n,ivec b)
227 #define SWPX(vv,xx,yy) { int tmp; tmp=vv[xx]; vv[xx] = vv[yy]; vv[yy] = tmp; }
228 copy_ivec(b,BB[n++]); /* x y z */
230 copy_ivec(b,BB[n++]); /* y x z */
232 copy_ivec(b,BB[n++]); /* z x y */
234 copy_ivec(b,BB[n++]); /* x z y */
236 copy_ivec(b,BB[n++]); /* y z x */
238 copy_ivec(b,BB[n++]); /* z y x */
239 SWPX(b,XX,ZZ); /* Back to normal */
244 static real box_weight(ivec nbox,matrix box)
249 /* Calculate area of subbox */
250 for(m=0; (m<DIM); m++)
251 lx[m] = box[m][m]/nbox[m];
252 return 2*(lx[XX]*lx[YY]+lx[XX]*lx[ZZ]+lx[YY]*lx[ZZ]);
255 static int w_comp(const void *a,const void *b)
263 wa = box_weight(ia,BOX);
264 wb = box_weight(ib,BOX);
265 if (fabs(wa - wb) < 1e-4)
266 return (iiprod(ia,ia) - iiprod(ib,ib));
273 static void buildbox(int nnode,ivec nbox,matrix box)
276 int i,j,m,n,n3,ny,*fx,*fy,nbb;
278 n3 = ipow(nnode,3)*6;
284 for(i=0; (i<=nnode); i++) {
285 for(m=1; (m<=fx[i]); m++) {
286 bxyz[XX] = ipow(i,m);
289 for(j=0; (j<=ny); j++) {
290 for(n=1; (n<=fy[j]); n++) {
291 bxyz[YY] = ipow(j,n);
292 bxyz[ZZ] = ny/bxyz[YY];
294 nbb = add_bb(BB,nbb,bxyz);
300 /* Sort boxes and remove doubles */
301 qsort(BB,nbb,sizeof(BB[0]),iv_comp);
303 for(i=1; (i<nbb); i++) {
304 if ((BB[i][XX] != BB[j][XX]) ||
305 (BB[i][YY] != BB[j][YY]) ||
306 (BB[i][ZZ] != BB[j][ZZ])) {
308 copy_ivec(BB[i],BB[j]);
312 /* Sort boxes according to weight */
314 qsort(BB,nbb,sizeof(BB[0]),w_comp);
315 for(i=0; (i<nbb); i++) {
316 fprintf(stderr,"nbox = %2d %2d %2d [ prod %3d ] area = %12.5f (nm^2)\n",
317 BB[i][XX],BB[i][YY],BB[i][ZZ],
318 BB[i][XX]*BB[i][YY]*BB[i][ZZ],
319 box_weight(BB[i],box));
321 copy_ivec(BB[0],nbox);
327 void mkcompact(int astart,int nwater,int nwatom,rvec x[],rvec v[],
328 int nnode,matrix box)
330 /* Make a compact configuration for each processor.
331 * Divide the computational box in near cubic boxes and spread them
332 * evenly over processors.
340 buildbox(nnode,NBOX,box);
341 /* copy_ivec(nbox,NBOX); */
342 for(m=0; (m<DIM); m++)
343 box_1[m] = 1.0/box[m][m];
345 lo_sortwater(astart,nwater,nwatom,x,v,TRUE);