src/mdlib/fftgrid.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  * GROwing Monsters And Cloning Shrimps
  35  */
  36 /* This file is completely threadsafe - keep it that way! */
  37 #ifdef HAVE_CONFIG_H
  38 #include <config.h>
  39 #endif
  40
  41 #include <stdlib.h>
  42 #include <stdio.h>
  43
  44 #include "typedefs.h"
  45 #include "futil.h"
  46 #include "smalloc.h"
  47 #include "futil.h"
  48 #include "macros.h"
  49 #include "network.h"
  50 #include "fftgrid.h"
  51 #include "gmx_fft.h"
  52 #include "gmx_parallel_3dfft.h"
  53 #include "gmxfio.h"
  54
  55
  56 #ifdef GMX_MPI
  57 static void print_parfft(FILE *fp,char *title,t_parfft *pfft)
  58 {
  59   fprintf(fp,"PARALLEL FFT DATA:\n"
  60           "   local_nx:                 %3d  local_x_start:                 %3d\n"
  61           "   local_ny_after_transpose: %3d  local_y_start_after_transpose  %3d\n",
  62           pfft->local_nx,pfft->local_x_start,pfft->local_ny_after_transpose,
  63           pfft->local_y_start_after_transpose);
  64 }
  65 #endif
  66
  67
  68 void *
  69 gmx_alloc_aligned(size_t size)
  70 {
  71     void *p0,*p;
  72
  73     p0 = malloc(size+32);
  74
  75     if(p0 == NULL)
  76     {
  77         gmx_fatal(FARGS,"Failed to allocated %u bytes of aligned memory.",size+32);
  78     }
  79
  80     p = (void *) (((size_t) p0 + 32) & (~((size_t) 31)));
  81
  82     /* Yeah, yeah, we cannot free this pointer, but who cares... */
  83     return p;
  84 }
  85
  86 t_fftgrid *mk_fftgrid(int          nx,
  87                       int          ny,
  88                       int          nz,
  89                       int          *node2slab,
  90                       int          *slab2grid_x,
  91                       t_commrec *  cr,
  92                       bool         bReproducible)
  93 {
  94 /* parallel runs with non-parallel ffts haven't been tested yet */
  95     int           nnodes;
  96     int           x1,y1,maxlocalsize;
  97     t_fftgrid *   grid;
  98     int           flags;
  99
 100     nnodes = 1;
 101 #ifdef GMX_MPI
 102     if (cr && cr->nnodes > 1) {
 103         MPI_Comm_size(cr->mpi_comm_mygroup,&nnodes);
 104     }
 105 #endif
 106
 107     snew(grid,1);
 108     grid->nx   = nx;
 109     grid->ny   = ny;
 110     grid->nz   = nz;
 111     grid->nxyz = nx*ny*nz;
 112     grid->bParallel = (nnodes > 1);
 113
 114     if (grid->bParallel)
 115     {
 116         grid->la2r = (nz/2+1)*2;
 117     }
 118     else
 119     {
 120         grid->la2r = nz;
 121     }
 122
 123     grid->la2c = (nz/2+1);
 124
 125     grid->la12r = ny*grid->la2r;
 126
 127     if (grid->bParallel)
 128     {
 129         grid->la12c = nx*grid->la2c;
 130     }
 131     else
 132     {
 133         grid->la12c = ny*grid->la2c;
 134     }
 135
 136     /* This code assumes that the when the grid is not divisble by nnodes,
 137      * the maximum difference in local grid sizes is 1.
 138      */
 139     x1 = (nx % nnodes == 0 ? 0 : 1);
 140     y1 = (ny % nnodes == 0 ? 0 : 1);
 141
 142     grid->nptr = (nx + x1)*(ny + y1)*grid->la2c*2;
 143
 144     if (grid->bParallel)
 145     {
 146 #ifdef GMX_MPI
 147         gmx_parallel_3dfft_init(&grid->mpi_fft_setup,nx,ny,nz,
 148                                 node2slab,slab2grid_x,cr->mpi_comm_mygroup,
 149                                 bReproducible);
 150
 151         gmx_parallel_3dfft_limits(grid->mpi_fft_setup,
 152                                   &(grid->pfft.local_x_start),
 153                                   &(grid->pfft.local_nx),
 154                                   &(grid->pfft.local_y_start_after_transpose),
 155                                   &(grid->pfft.local_ny_after_transpose));
 156 #else
 157         gmx_fatal(FARGS,"Parallel FFT supported with MPI only!");
 158 #endif
 159     }
 160     else
 161     {
 162         gmx_fft_init_3d_real(&grid->fft_setup,nx,ny,nz,bReproducible ? GMX_FFT_FLAG_CONSERVATIVE : GMX_FFT_FLAG_NONE);
 163     }
 164     grid->ptr = (real *)gmx_alloc_aligned(grid->nptr*sizeof(*(grid->ptr)));
 165
 166 #ifdef GMX_MPI
 167     if (grid->bParallel && debug)
 168     {
 169         print_parfft(debug,"Plan", &grid->pfft);
 170     }
 171     if (grid->bParallel)
 172     {
 173         maxlocalsize = max((nx/nnodes + x1)*ny*grid->la2c*2,
 174                            (ny/nnodes + y1)*nx*grid->la2c*2);
 175         grid->workspace = (real *)
 176             gmx_alloc_aligned(maxlocalsize*sizeof(*(grid->workspace)));
 177     }
 178     else
 179     {
 180         grid->workspace =
 181             (real*)gmx_alloc_aligned(grid->nptr*sizeof(*(grid->workspace)));
 182     }
 183 #else /* no MPI */
 184     grid->workspace = (real *)gmx_alloc_aligned(grid->nptr*sizeof(*(grid->workspace)));
 185 #endif
 186
 187     return grid;
 188 }
 189
 190 void
 191 pr_fftgrid(FILE *fp,char *title,t_fftgrid *grid)
 192 {
 193   int     i,j,k,index_x,index_xy,ntot=0;
 194   int     nx,ny,nz,nx2,ny2,nz2,la12,la2;
 195   real *  ptr;
 196
 197   /* Unpack structure */
 198   unpack_fftgrid(grid,&nx,&ny,&nz,&nx2,&ny2,&nz2,&la2,&la12,TRUE,&ptr);
 199   for(i=0; (i<nx); i++) {
 200     index_x = la12*i;
 201     for(j=0; (j<ny); j++) {
 202       index_xy = index_x + la2*j;
 203       for(k=0; (k<nz); k++) {
 204         if (ptr[index_xy+k] != 0) {
 205           fprintf(fp,"%-12s  %5d  %5d  %5d  %12.5e\n",
 206                   title,i,j,k,ptr[index_xy+k]);
 207           ntot++;
 208         }
 209       }
 210     }
 211   }
 212   fprintf(fp,"%d non zero elements in %s\n",ntot,title);
 213 }
 214
 215 void done_fftgrid(t_fftgrid *grid)
 216 {
 217   /* memory can't be freed because it is allocated by gmx_alloc_aligned */
 218   if (grid->ptr) {
 219     /* sfree(grid->ptr); */
 220     grid->ptr = NULL;
 221   }
 222
 223   if (grid->workspace) {
 224     /* sfree(grid->workspace); */
 225     grid->workspace=NULL;
 226   }
 227 }
 228
 229
 230 void gmxfft3D(t_fftgrid *grid,enum gmx_fft_direction dir,t_commrec *cr)
 231 {
 232   real *tmp;
 233
 234   if (grid->bParallel)
 235   {
 236 #ifdef GMX_MPI
 237     if( dir == GMX_FFT_REAL_TO_COMPLEX || dir == GMX_FFT_COMPLEX_TO_REAL )
 238     {
 239       gmx_parallel_3dfft(grid->mpi_fft_setup,dir,grid->ptr,grid->ptr);
 240     }
 241     else
 242     {
 243         gmx_fatal(FARGS,"Invalid direction for FFT: %d",dir);
 244     }
 245 #else
 246     gmx_fatal(FARGS,"Parallel FFT supported with MPI only!");
 247 #endif
 248   }
 249   else
 250   {
 251     if( dir == GMX_FFT_REAL_TO_COMPLEX || dir == GMX_FFT_COMPLEX_TO_REAL)
 252     {
 253         gmx_fft_3d_real(grid->fft_setup,dir,grid->ptr,grid->workspace);
 254         tmp = grid->ptr;
 255         grid->ptr = grid->workspace;
 256         grid->workspace = tmp;
 257     }
 258     else
 259     {
 260       gmx_fatal(FARGS,"Invalid direction for FFT: %d",dir);
 261     }
 262   }
 263 }
 264
 265 void clear_fftgrid(t_fftgrid *grid)
 266 {
 267     /* clears the whole grid */
 268   int      i,ngrid;
 269   real *   ptr;
 270
 271   ngrid = grid->nptr;
 272   ptr   = grid->ptr;
 273
 274   for (i=0; (i<ngrid); i++) {
 275     ptr[i] = 0;
 276   }
 277 }
 278
 279 void unpack_fftgrid(t_fftgrid *grid,int *nx,int *ny,int *nz,
 280                     int *nx2,int *ny2,int *nz2,
 281                     int *la2,int *la12,bool bReal,real **ptr)
 282 {
 283   *nx  = grid->nx;
 284   *ny  = grid->ny;
 285   *nz  = grid->nz;
 286   *nx2 = 2*grid->nx;
 287   *ny2 = 2*grid->ny;
 288   *nz2 = 2*grid->nz;
 289   if(bReal) {
 290     *la2 = grid->la2r;
 291     *la12= grid->la12r;
 292   }
 293   else {
 294     *la2 = grid->la2c;
 295     *la12= grid->la12c;
 296   }
 297   *ptr = grid->ptr;
 298 }
 299
 300
 301
 302 /*****************************************************************
 303  *
 304  * For backward compatibility (for testing the ewald code vs. PPPM etc)
 305  * some old grid routines are retained here.
 306  *
 307  ************************************************************************/
 308
 309 real ***mk_rgrid(int nx,int ny,int nz)
 310 {
 311   real *ptr1;
 312   real **ptr2;
 313   real ***ptr3;
 314   int  i,j,n2,n3;
 315
 316   snew(ptr1,nx*ny*nz);
 317   snew(ptr2,nx*ny);
 318   snew(ptr3,nx);
 319
 320   n2=n3=0;
 321   for(i=0; (i<nx); i++) {
 322     ptr3[i]=&(ptr2[n2]);
 323     for(j=0; (j<ny); j++,n2++) {
 324       ptr2[n2] = &(ptr1[n3]);
 325       n3 += nz;
 326     }
 327   }
 328   return ptr3;
 329 }
 330
 331 void free_rgrid(real ***grid,int nx,int ny)
 332 {
 333   int i;
 334
 335   sfree(grid[0][0]);
 336   for(i=0; (i<nx); i++) {
 337     sfree(grid[i]);
 338   }
 339   sfree(grid);
 340 }
 341
 342 real print_rgrid(FILE *fp,char *title,int nx,int ny,int nz,real ***grid)
 343 {
 344   int  ix,iy,iz;
 345   real g,gtot;
 346
 347   gtot=0;
 348   if (fp)
 349     fprintf(fp,"Printing all non-zero real elements of %s\n",title);
 350   for(ix=0; (ix<nx); ix++)
 351     for(iy=0; (iy<ny); iy++)
 352       for(iz=0; (iz<nz); iz++) {
 353         g=grid[ix][iy][iz];
 354         if (fp && (g != 0))
 355           fprintf(fp,"%s[%2d][%2d][%2d] = %12.5e\n",title,ix,iy,iz,g);
 356         gtot+=g;
 357       }
 358   return gtot;
 359 }
 360
 361 void print_rgrid_pdb(char *fn,int nx,int ny,int nz,real ***grid)
 362 {
 363   FILE *fp;
 364   int  ix,iy,iz,n,ig;
 365   real x,y,z,g;
 366
 367   n=1;
 368   fp=gmx_fio_fopen(fn,"w");
 369   for(ix=0; (ix<nx); ix++) {
 370     for(iy=0; (iy<ny); iy++) {
 371       for(iz=0; (iz<nz); iz++) {
 372         g=grid[ix][iy][iz];
 373         ig=g;
 374         if ((ig != 0) || (1)) {
 375           x = 4*ix;
 376           y = 4*iy;
 377           z = 4*iz;
 378           fprintf(fp,"ATOM  %5d  Na   Na     1    %8.3f%8.3f%8.3f%6.2f%6.2f\n",
 379                   n++,x,y,z,0.0,g);
 380         }
 381       }
 382     }
 383   }
 384   gmx_fio_fclose(fp);
 385 }
 386
 387 void clear_rgrid(int nx,int ny,int nz,real ***grid)
 388 {
 389   int i,j,k;
 390
 391   for(i=0; (i<nx); i++)
 392     for(j=0; (j<ny); j++)
 393       for(k=0; (k<nz); k++)
 394         grid[i][j][k] = 0;
 395 }
 396
 397 void clear_cgrid(int nx,int ny,int nz,t_complex ***grid)
 398 {
 399   int i,j,k;
 400
 401   for(i=0; (i<nx); i++)
 402     for(j=0; (j<ny); j++)
 403       for(k=0; (k<nz); k++)
 404         grid[i][j][k] = cnul;
 405 }
 406
 407 t_complex ***mk_cgrid(int nx,int ny,int nz)
 408 {
 409   t_complex *ptr1;
 410   t_complex **ptr2;
 411   t_complex ***ptr3;
 412   int  i,j,n2,n3;
 413
 414   snew(ptr1,nx*ny*nz);
 415   snew(ptr2,nx*ny);
 416   snew(ptr3,nx);
 417
 418   n2=n3=0;
 419   for(i=0; (i<nx); i++) {
 420     ptr3[i]=&(ptr2[n2]);
 421     for(j=0; (j<ny); j++,n2++) {
 422       ptr2[n2] = &(ptr1[n3]);
 423       n3 += nz;
 424     }
 425   }
 426   return ptr3;
 427 }
 428
 429 void free_cgrid(t_complex ***grid,int nx,int ny)
 430 {
 431   int i;
 432
 433   sfree(grid[0][0]);
 434   for(i=0; (i<nx); i++)
 435     sfree(grid[i]);
 436   sfree(grid);
 437 }
 438
 439 t_complex print_cgrid(FILE *fp,char *title,int nx,int ny,int nz,
 440                       t_complex ***grid)
 441 {
 442   int     ix,iy,iz;
 443   t_complex g,gtot;
 444
 445   gtot=cnul;
 446   if (fp)
 447     fprintf(fp,"Printing all non-zero complex elements of %s\n",title);
 448   for(ix=0; (ix<nx); ix++)
 449     for(iy=0; (iy<ny); iy++)
 450       for(iz=0; (iz<nz); iz++) {
 451         g=grid[ix][iy][iz];
 452         if (fp  && ((g.re != 0) || (g.im != 0)))
 453           fprintf(fp,"%s[%2d][%2d][%2d] = %12.5e + i %12.5e\n",
 454                   title,ix,iy,iz,g.re,g.im);
 455         gtot = cadd(gtot,g);
 456       }
 457   return gtot;
 458 }
 459
 460 void print_cgrid_pdb(char *fn,int nx,int ny,int nz,t_complex ***grid)
 461 {
 462   FILE *fp;
 463   int  ix,iy,iz,n;
 464   real x,y,z,g;
 465
 466   n=1;
 467   fp=gmx_fio_fopen(fn,"w");
 468   for(ix=0; (ix<nx); ix++) {
 469     for(iy=0; (iy<ny); iy++) {
 470       for(iz=0; (iz<nz); iz++) {
 471         g=grid[ix][iy][iz].re;
 472         if (g != 0) {
 473           x = 4*ix;
 474           y = 4*iy;
 475           z = 4*iz;
 476           fprintf(fp,"ATOM  %5d  Na   Na     1    %8.3f%8.3f%8.3f%6.2f%6.2f\n",
 477                   n++,x,y,z,0.0,g);
 478         }
 479       }
 480     }
 481   }
 482   gmx_fio_fclose(fp);
 483 }
 484
 485
 486
 487
 488
 489
 490
 491