src/gromacs/gmxlib/calcgrid.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.2.0
  11  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  12  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  13  * Copyright (c) 2001-2004, The GROMACS development team,
  14  * check out http://www.gromacs.org for more information.
  15
  16  * This program is free software; you can redistribute it and/or
  17  * modify it under the terms of the GNU General Public License
  18  * as published by the Free Software Foundation; either version 2
  19  * of the License, or (at your option) any later version.
  20  *
  21  * If you want to redistribute modifications, please consider that
  22  * scientific software is very special. Version control is crucial -
  23  * bugs must be traceable. We will be happy to consider code for
  24  * inclusion in the official distribution, but derived work must not
  25  * be called official GROMACS. Details are found in the README & COPYING
  26  * files - if they are missing, get the official version at www.gromacs.org.
  27  *
  28  * To help us fund GROMACS development, we humbly ask that you cite
  29  * the papers on the package - you can find them in the top README file.
  30  *
  31  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * GROningen Mixture of Alchemy and Childrens' Stories
  35  */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include <math.h>
  41
  42 #include "typedefs.h"
  43 #include "smalloc.h"
  44 #include "gmx_fatal.h"
  45 #include "calcgrid.h"
  46
  47 #define facNR 4
  48 const int factor[facNR] = {2,3,5,7};
  49
  50 static void make_list(int start_fac,int *ng,int ng_max,int *n_list,int **list)
  51 {
  52     int i,fac;
  53
  54     if (*ng < ng_max)
  55     {
  56         if (*n_list % 100 == 0)
  57         {
  58             srenew(*list,*n_list+100);
  59         }
  60         (*list)[*n_list] = *ng;
  61         (*n_list)++;
  62
  63         for(i=start_fac; i<facNR; i++)
  64         {
  65             fac = factor[i];
  66             /* The choice of grid size is based on benchmarks of fftw
  67              * and the need for a lot of factors for nice DD decomposition.
  68              * The base criterion is that a grid size is not included
  69              * when there is a larger grid size that produces a faster 3D FFT.
  70              * Allow any power for 2, two for 3 and 5, but only one for 7.
  71              * Three for 3 are ok when there is also a factor of 2.
  72              * Two factors of 5 are not allowed with a factor of 3 or 7.
  73              * A factor of 7 does not go with a factor of 5, 7 or 9.
  74              */
  75             if ((fac == 2) ||
  76                 (fac == 3 && (*ng % 9 != 0 ||
  77                               (*ng % 2 == 0 && *ng % 27 != 0))) ||
  78                 (fac == 5 && *ng % 15 != 0 && *ng % 25 != 0) ||
  79                 (fac == 7 && *ng % 5 != 0 && *ng % 7 != 0 && *ng % 9 != 0))
  80             {
  81                 *ng *= fac;
  82                 make_list(i,ng,ng_max,n_list,list);
  83                 *ng /= fac;
  84             }
  85         }
  86     }
  87 }
  88
  89 static int list_comp(const void *a,const void *b)
  90 {
  91   return (*((int *)a) - *((int *)b));
  92 }
  93
  94 real calc_grid(FILE *fp,matrix box,real gr_sp,
  95                int *nx,int *ny,int *nz)
  96 {
  97     int  d,n[DIM];
  98     int  i,j,nmin[DIM];
  99     rvec box_size,spacing;
 100     real max_spacing;
 101     int  ng_max,ng;
 102     int  n_list,*list;
 103
 104     if (gr_sp <= 0)
 105     {
 106         gmx_fatal(FARGS,"invalid fourier grid spacing: %g",gr_sp);
 107     }
 108
 109     /* New grid calculation setup:
 110      *
 111      * To maintain similar accuracy for triclinic PME grids as for rectangular
 112      * ones, the max grid spacing should set along the box vectors rather than
 113      * cartesian X/Y/Z directions. This will lead to slightly larger grids, but
 114      * it is much better than having to go to pme_order=6.
 115      *
 116      * Thus, instead of just extracting the diagonal elements to box_size[d], we
 117      * now calculate the cartesian length of the vectors.
 118      *
 119      * /Erik Lindahl, 20060402.
 120      */
 121     for(d=0; d<DIM; d++)
 122     {
 123         box_size[d] = 0;
 124         for(i=0;i<DIM;i++)
 125         {
 126             box_size[d] += box[d][i]*box[d][i];
 127         }
 128         box_size[d] = sqrt(box_size[d]);
 129     }
 130
 131     n[XX] = *nx;
 132     n[YY] = *ny;
 133     n[ZZ] = *nz;
 134
 135     ng = 1;
 136     ng_max = 1;
 137     for(d=0; d<DIM; d++)
 138     {
 139         nmin[d] = (int)(box_size[d]/gr_sp + 0.999);
 140         if (2*nmin[d] > ng_max)
 141         {
 142             ng_max = 2*nmin[d];
 143         }
 144     }
 145     n_list=0;
 146     list=NULL;
 147     make_list(0,&ng,ng_max,&n_list,&list);
 148
 149     if ((*nx<=0) || (*ny<=0) || (*nz<=0))
 150     {
 151         if (NULL != fp)
 152         {
 153             fprintf(fp,"Calculating fourier grid dimensions for%s%s%s\n",
 154                     *nx > 0 ? "":" X",*ny > 0 ? "":" Y",*nz > 0 ? "":" Z");
 155         }
 156     }
 157
 158     qsort(list,n_list,sizeof(list[0]),list_comp);
 159     if (debug)
 160     {
 161         for(i=0; i<n_list; i++)
 162             fprintf(debug,"grid: %d\n",list[i]);
 163     }
 164
 165     for(d=0; d<DIM; d++)
 166     {
 167         for(i=0; (i<n_list) && (n[d]<=0); i++)
 168         {
 169             if (list[i] >= nmin[d])
 170             {
 171                 n[d] = list[i];
 172             }
 173         }
 174     }
 175
 176     sfree(list);
 177
 178     max_spacing = 0;
 179     for(d=0; d<DIM; d++)
 180     {
 181         spacing[d] = box_size[d]/n[d];
 182         if (spacing[d] > max_spacing)
 183         {
 184             max_spacing = spacing[d];
 185         }
 186     }
 187     *nx = n[XX];
 188     *ny = n[YY];
 189     *nz = n[ZZ];
 190     if (NULL != fp)
 191     {
 192         fprintf(fp,"Using a fourier grid of %dx%dx%d, spacing %.3f %.3f %.3f\n",
 193                 *nx,*ny,*nz,spacing[XX],spacing[YY],spacing[ZZ]);
 194     }
 195
 196     return max_spacing;
 197 }
 198
 199