-
- sfree(ir);
-}
-
-
-typedef struct
-{
- int nkx, nky, nkz; /* Fourier grid */
- real fac; /* actual scaling factor */
- real fs; /* Fourierspacing */
-} t_pmegrid;
-
-
-static void copy_grid(t_pmegrid *ingrid, t_pmegrid *outgrid)
-{
- outgrid->nkx = ingrid->nkx;
- outgrid->nky = ingrid->nky;
- outgrid->nkz = ingrid->nkz;
- outgrid->fac = ingrid->fac;
- outgrid->fs = ingrid->fs;
-}
-
-/* Removes entry 'index' from the t_pmegrid list */
-static void remove_from_list(
- t_pmegrid gridlist[],
- int *nlist, /* Length of the list */
- int index) /* Index to remove from the list */
-{
- int j;
-
-
- for (j = index; j < (*nlist - 1); j++)
- {
- copy_grid(&gridlist[j+1], &gridlist[j]);
- }
- *nlist = *nlist - 1;
-}
-
-
-/* Returns the index of the least necessary grid in the list.
- *
- * This is the one where the following conditions hold for the scaling factor:
- * 1. this grid has the smallest distance to its neighboring grid (distance
- * measured by fac) -> this condition is true for two grids at the same time
- * 2. this grid (of the two) has the smaller distance to the other neighboring
- * grid
- *
- * The extreme grids (the ones with the smallest and largest
- * scaling factor) are never thrown away.
- */
-static int get_throwaway_index(
- t_pmegrid gridlist[],
- int nlist) /* Length of the list */
-{
- int i,index = -1;
- real dist,mindist,d_left,d_right;
-
-
- /* Find the two factors with the smallest mutual distance */
- mindist = GMX_FLOAT_MAX;
- for (i = 1; i < nlist; i++)
- {
- dist = fabs(gridlist[i].fac - gridlist[i-1].fac);
- if (dist < mindist)
- {
- index = i;
- mindist = dist;
- }
- }
- /* index and index-1 have the smallest mutual distance */
- if (index == 1)
- {
- /* Never return the first index, i.e. index == 0 */
- ;
- }
- else
- {
- d_left = fabs(gridlist[index-1].fac - gridlist[index-2].fac);
- d_right = fabs(gridlist[index+1].fac - gridlist[index ].fac);
-
- /* Return the left index if its distance to its left neighbor is shorter
- * than the distance of the right index to its right neighbor */
- if (d_left < d_right)
- index--;
- }
- /* Never return the last index */
- if (index == nlist-1)
- index--;
-
- return index;
-}
-
-
-static void make_grid_list(
- real fmin, /* minimum scaling factor (downscaling fac) */
- real fmax, /* maximum scaling factor (upscaling fac) */
- int *ntprs, /* Number of tpr files to test */
- matrix box, /* The box */
- t_pmegrid *griduse[], /* List of grids that have to be tested */
- real fs) /* Requested fourierspacing at a scaling factor
- of unity */
-{
- real req_fac,act_fac=0,act_fs,eps;
- rvec box_size;
- int i,jmin=0,d,ngridall=0;
- int nkx=0,nky=0,nkz=0;
- int nkx_old=0,nky_old=0,nkz_old=0;
- t_pmegrid *gridall;
- int gridalloc,excess;
-
-
- /* Determine length of triclinic box vectors */
- for(d=0; d<DIM; d++)
- {
- box_size[d] = 0;
- for(i=0;i<DIM;i++)
- box_size[d] += box[d][i]*box[d][i];
- box_size[d] = sqrt(box_size[d]);
- }
-
- gridalloc = 25;
- snew(gridall, gridalloc);
-
- fprintf(stdout, "Possible PME grid settings (apart from input file settings):\n");
- fprintf(stdout, " nkx nky nkz max spacing scaling factor\n");
-
- /* eps should provide a fine enough spacing not to miss any grid */
- if (fmax != fmin)
- eps = (fmax-fmin)/(100.0*(*ntprs - 1));
- else
- eps = 1.0/max( (*griduse)[0].nkz, max( (*griduse)[0].nkx, (*griduse)[0].nky ) );
-
- for (req_fac = fmin; act_fac < fmax; req_fac += eps)
- {
- nkx=0;
- nky=0;
- nkz=0;
- calc_grid(NULL,box,fs*req_fac,&nkx,&nky,&nkz);
- act_fs = max(box_size[XX]/nkx,max(box_size[YY]/nky,box_size[ZZ]/nkz));
- act_fac = act_fs/fs;
- if ( ! ( nkx==nkx_old && nky==nky_old && nkz==nkz_old ) /* Exclude if grid is already in list */
- && ! ( nkx==(*griduse)[0].nkx && nky==(*griduse)[0].nky && nkz==(*griduse)[0].nkz ) ) /* Exclude input file grid */
- {
- /* We found a new grid that will do */
- nkx_old = nkx;
- nky_old = nky;
- nkz_old = nkz;
- gridall[ngridall].nkx = nkx;
- gridall[ngridall].nky = nky;
- gridall[ngridall].nkz = nkz;
- gridall[ngridall].fac = act_fac;
- gridall[ngridall].fs = act_fs;
- fprintf(stdout, "%5d%5d%5d %12f %12f\n",nkx,nky,nkz,act_fs,act_fac);
- ngridall++;
- if (ngridall >= gridalloc)
- {
- gridalloc += 25;
- srenew(gridall, gridalloc);
- }
- }
- }
-
- /* excess is the number of grids we have to get rid of */
- excess = ngridall+1 - *ntprs;
-
- /* If we found less grids than tpr files were requested, simply test all
- * the grid there are ... */
- if (excess < 0)
- {
- fprintf(stdout, "NOTE: You requested %d tpr files, but apart from the input file grid only the\n"
- " above listed %d suitable PME grids were found. Will test all suitable settings.\n",
- *ntprs, ngridall);
- *ntprs = ngridall+1;
- }
- else
- {
- if (2 == *ntprs)
- {
- /* We can only keep the input tpr file plus one extra tpr file.
- * We make that choice depending on the values of -upfac and -downfac */
- if (fmin < 1.0)
- {
- /* Keep the one with the -downfac as scaling factor. This is already
- * stored in gridall[0] */
- ;
- }
- else
- {
- /* Keep the one with -upfac as scaling factor */
- copy_grid(&(gridall[ngridall-1]), &(gridall[0]));
- }
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff