src/gmxlib/calcgrid.c

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