src/mdlib/pppm.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  * GROwing Monsters And Cloning Shrimps
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <stdio.h>
  40 #include <math.h>
  41 #include "physics.h"
  42 #include "typedefs.h"
  43 #include "smalloc.h"
  44 #include "vec.h"
  45 #include "xvgr.h"
  46 #include "gmx_fatal.h"
  47 #include "txtdump.h"
  48 #include "network.h"
  49 #include "nrnb.h"
  50 #include "pppm.h"
  51 #include "coulomb.h"
  52 #include "mdrun.h"
  53 #include "gmx_fft.h"
  54 #include "pme.h"
  55
  56 #ifdef GMX_MPI
  57 #include "gmx_parallel_3dfft.h"
  58 #endif
  59
  60 #define llim  (-1)
  61 #define ulim   (1)
  62 #define llim2 (-3)
  63 #define ulim2  (3)
  64
  65
  66
  67 /* PPPM temporarily disabled while working on 2D PME */
  68 #define DISABLE_PPPM
  69
  70
  71
  72 #ifndef DISABLE_PPPM
  73
  74 /* TODO: fix thread-safety */
  75
  76 static void calc_invh(rvec box,int nx,int ny,int nz,rvec invh)
  77 {
  78   invh[XX] = nx/box[XX];
  79   invh[YY] = ny/box[YY];
  80   invh[ZZ] = nz/box[ZZ];
  81 }
  82
  83 static void calc_weights(int iatom,int nx,int ny,int nz,
  84                          rvec x,rvec box,rvec invh,ivec ixyz,real WXYZ[])
  85 {
  86   const  real half=0.5;
  87   tensor wxyz;
  88   real   abc,ttt,fact;
  89 #ifdef DEBUG
  90   real   wtot;
  91 #endif
  92   ivec   nxyz;
  93   int    it,j,k,m,nm;
  94   real   Wx,Wy,Wzx,Wzy,Wzz;
  95
  96   fact = 0.125;
  97
  98   nxyz[XX] = nx;
  99   nxyz[YY] = ny;
 100   nxyz[ZZ] = nz;
 101   for(m=0; (m<DIM); m++) {
 102     /* Put particle in the box... */
 103     ttt = x[m]*invh[m];
 104
 105     /* Round to nearest grid point. Do the math in integer! */
 106     it  = (ttt+0.5);
 107     nm  = nxyz[m];
 108     if (it < 0) {
 109       ttt+=nm;
 110       it +=nm;
 111     }
 112     else if (it >= nm) {
 113       ttt-=nm;
 114       it -=nm;
 115     }
 116     if ((it < 0) || (it >= nxyz[m]))
 117       gmx_fatal(FARGS,"iatom = %d, it = %d, x=%f, ttt=%f",iatom,it,x[m],ttt);
 118     ixyz[m]    = it;
 119
 120     /* Fraction (offset) from grid point */
 121     abc        = ttt - (real)ixyz[m];
 122
 123     wxyz[m][0] = sqr((real)(half  - abc));
 124     wxyz[m][1] = 1.5 - 2.0*sqr(abc);
 125     wxyz[m][2] = sqr((real)(half  + abc));
 126   }
 127   Wzx=wxyz[ZZ][XX];
 128   Wzy=wxyz[ZZ][YY];
 129   Wzz=wxyz[ZZ][ZZ];
 130   for(j=m=0; (j<DIM); j++) {
 131     Wx = wxyz[XX][j]*fact;
 132     for(k=0; (k<DIM); k++,m+=3) {
 133       Wy        = Wx*wxyz[YY][k];
 134       WXYZ[m]   = Wy*Wzx;
 135       WXYZ[m+1] = Wy*Wzy;
 136       WXYZ[m+2] = Wy*Wzz;
 137     }
 138   }
 139 #ifdef DEBUG
 140   wtot = 0;
 141   for(j=0; (j<27); j++)
 142     wtot+=WXYZ[j];
 143   fprintf(stderr,"wtot = %g\n",wtot);
 144 #endif
 145 #ifdef HACK
 146   for(j=0; (j<27); j++)
 147     WXYZ[j] = 0;
 148   WXYZ[13] = 1.0;
 149 #endif
 150 }
 151
 152 static void calc_nxyz(int nx,int ny,int nz,
 153                       int **nnx,int **nny,int **nnz)
 154 {
 155   int i;
 156
 157   snew(*nnx,3*nx);
 158   snew(*nny,3*ny);
 159   snew(*nnz,3*nz);
 160   for(i=0; (i<3*nx); i++)
 161     (*nnx)[i] = i % nx;
 162   for(i=0; (i<3*ny); i++)
 163     (*nny)[i] = i % ny;
 164   for(i=0; (i<3*nz); i++)
 165     (*nnz)[i] = i % nz;
 166 }
 167
 168 static void spread_q(FILE *log,bool bVerbose,
 169                      int start,int nr,
 170                      rvec x[],real charge[],rvec box,
 171                      t_fftgrid *grid,t_nrnb *nrnb)
 172 {
 173   static bool bFirst = TRUE;
 174   static int  *nnx,*nny,*nnz;
 175   rvec   invh;
 176   real   qi,qwt;
 177 #ifdef DEBUG
 178   real   qt;
 179 #endif
 180   real   WXYZ[27];
 181   ivec   ixyz;
 182   int    i,iX,iY,iZ,index;
 183   int    jx,jy,jz,jcx,jcy,jcz;
 184   int    nxyz;
 185   int    nx,ny,nz,nx2,ny2,nz2,la2,la12;
 186   real *ptr;
 187
 188   unpack_fftgrid(grid,&nx,&ny,&nz,&nx2,&ny2,&nz2,&la2,&la12,TRUE,&ptr);
 189
 190   calc_invh(box,nx,ny,nz,invh);
 191
 192   if (bFirst) {
 193     if (log) {
 194       fprintf(log,
 195               "Spreading Charges using Triangle Shaped on %dx%dx%d grid\n",
 196               nx,ny,nz);
 197       fprintf(log,"invh = %10g,%10g,%10g\n",invh[XX],invh[YY],invh[ZZ]);
 198     }
 199
 200     calc_nxyz(nx,ny,nz,&nnx,&nny,&nnz);
 201
 202     bFirst = FALSE;
 203   }
 204
 205   for(i=start; (i<start+nr); i++) {
 206     qi=charge[i];
 207
 208     /* Each charge is spread over the nearest 27 grid cells,
 209      * thus we loop over -1..1 in X,Y and Z direction
 210      * We apply the TSC (triangle shaped charge)
 211      * see Luty et. al, JCP 103 (1995) 3014
 212      */
 213
 214     if (qi != 0.0) {
 215       calc_weights(i,nx,ny,nz,x[i],box,invh,ixyz,WXYZ);
 216       iX  = ixyz[XX] + nx;
 217       iY  = ixyz[YY] + ny;
 218       iZ  = ixyz[ZZ] + nz;
 219
 220 #ifdef DEBUG
 221       qt=0;
 222 #endif
 223       nxyz = 0;
 224       for(jx=-1; (jx<=1); jx++) {
 225         jcx = nnx[iX + jx];
 226         for(jy=-1; (jy<=1); jy++) {
 227           jcy = nny[iY + jy];
 228           for(jz=-1; (jz<=1); jz++,nxyz++) {
 229             jcz   = nnz[iZ + jz];
 230             index = INDEX(jcx,jcy,jcz);
 231             qwt   = qi*WXYZ[nxyz];
 232             grid->ptr[index]+=qwt;
 233 #ifdef DEBUG
 234             qt   += qwt;
 235             if (bVerbose && log)
 236               fprintf(log,"spread %4d %2d %2d %2d  %10.3e, weight=%10.3e\n",
 237                       index,jcx,jcy,jcz,grid->ptr[index],WXYZ[nxyz]);
 238 #endif
 239           }
 240         }
 241       }
 242 #ifdef DEBUG
 243       if (log)
 244         fprintf(log,"q[%3d] = %6.3f, qwt = %6.3f\n",i,qi,qt);
 245 #endif
 246     }
 247   }
 248   inc_nrnb(nrnb,eNR_SPREADQ,9*nr);
 249   inc_nrnb(nrnb,eNR_WEIGHTS,3*nr);
 250 }
 251
 252 static real gather_inner(int JCXYZ[],real WXYZ[],int ixw[],int iyw[],int izw[],
 253                          int la2,int la12,
 254                   real c1x,real c1y,real c1z,real c2x,real c2y,real c2z,
 255                   real qi,rvec f,real ptr[])
 256 {
 257   real pi,fX,fY,fZ,weight;
 258   int  jxyz,m,jcx,jcy,jcz;
 259   int  jcx0,jcy0,jcz0;
 260
 261   pi = 0.0;
 262   fX = 0.0;
 263   fY = 0.0;
 264   fZ = 0.0;
 265
 266   /* Now loop over 27 surrounding vectors */
 267   for(jxyz=m=0; (jxyz < 27); jxyz++,m+=3) {
 268     jcx    = JCXYZ[m];
 269     jcy    = JCXYZ[m+1];
 270     jcz    = JCXYZ[m+2];
 271     weight = WXYZ[jxyz];
 272
 273     jcx0   = ixw[jcx];
 274     jcy0   = iyw[jcy];
 275     jcz0   = izw[jcz];
 276
 277     /* Electrostatic Potential ! */
 278     pi += weight * ptr[INDEX(jcx0,jcy0,jcz0)];
 279
 280     /* Forces */
 281     fX += weight * ((c1x*(ptr[INDEX(ixw[jcx-1],jcy0,jcz0)] -
 282                           ptr[INDEX(ixw[jcx+1],jcy0,jcz0)] )) +
 283                     (c2x*(ptr[INDEX(ixw[jcx-2],jcy0,jcz0)] -
 284                           ptr[INDEX(ixw[jcx+2],jcy0,jcz0)] )));
 285     fY += weight * ((c1y*(ptr[INDEX(jcx0,iyw[jcy-1],jcz0)] -
 286                           ptr[INDEX(jcx0,iyw[jcy+1],jcz0)] ))  +
 287                     (c2y*(ptr[INDEX(jcx0,iyw[jcy-2],jcz0)] -
 288                           ptr[INDEX(jcx0,iyw[jcy+2],jcz0)] )));
 289     fZ += weight * ((c1z*(ptr[INDEX(jcx0,jcy0,izw[jcz-1])] -
 290                           ptr[INDEX(jcx0,jcy0,izw[jcz+1])] ))  +
 291                     (c2z*(ptr[INDEX(jcx0,jcy0,izw[jcz-2])] -
 292                           ptr[INDEX(jcx0,jcy0,izw[jcz+2])] )));
 293   }
 294   f[XX] += qi*fX;
 295   f[YY] += qi*fY;
 296   f[ZZ] += qi*fZ;
 297
 298   return pi;
 299 }
 300
 301 static real gather_f(FILE *log,bool bVerbose,
 302                      int start,int nr,rvec x[],rvec f[],real charge[],rvec box,
 303                      real pot[],t_fftgrid *grid,rvec beta,t_nrnb *nrnb)
 304 {
 305   static bool bFirst=TRUE;
 306   static int  *nnx,*nny,*nnz;
 307   static int  JCXYZ[81];
 308   int    i,m;
 309   real   energy;
 310   real   qi,pi;
 311   ivec   ixyz;
 312   rvec   invh,c1,c2;
 313   real   WXYZ[27];
 314   real   c1x,c1y,c1z,c2x,c2y,c2z;
 315   int    ixw[7],iyw[7],izw[7];
 316   int    ll;
 317   int    nx,ny,nz,nx2,ny2,nz2,la2,la12;
 318   real *ptr;
 319
 320   unpack_fftgrid(grid,&nx,&ny,&nz,&nx2,&ny2,&nz2,&la2,&la12,TRUE,&ptr);
 321
 322   calc_invh(box,nx,ny,nz,invh);
 323
 324   for(m=0; (m<DIM); m++) {
 325     c1[m] = (beta[m]/2.0)*invh[m];
 326     c2[m] = ((1.0-beta[m])/4.0)*invh[m];
 327   }
 328   c1x = c1[XX];
 329   c1y = c1[YY];
 330   c1z = c1[ZZ];
 331   c2x = c2[XX];
 332   c2y = c2[YY];
 333   c2z = c2[ZZ];
 334
 335   if (bFirst) {
 336     if (log) {
 337       fprintf(log,"Gathering Forces using Triangle Shaped on %dx%dx%d grid\n",
 338               nx,ny,nz);
 339       fprintf(log,"beta = %10g,%10g,%10g\n",beta[XX],beta[YY],beta[ZZ]);
 340       fprintf(log,"c1   = %10g,%10g,%10g\n",c1[XX],c1[YY],c1[ZZ]);
 341       fprintf(log,"c2   = %10g,%10g,%10g\n",c2[XX],c2[YY],c2[ZZ]);
 342       fprintf(log,"invh = %10g,%10g,%10g\n",invh[XX],invh[YY],invh[ZZ]);
 343     }
 344     calc_nxyz(nx,ny,nz,&nnx,&nny,&nnz);
 345
 346     for(i=0; (i<27); i++) {
 347       JCXYZ[3*i]   = 2 + (i/9);
 348       JCXYZ[3*i+1] = 2 + (i/3) % 3;
 349       JCXYZ[3*i+2] = 2 + (i % 3);
 350     }
 351
 352     bFirst = FALSE;
 353   }
 354
 355   energy=0.0;
 356   for(i=start; (i<start+nr); i++) {
 357     /* Each charge is spread over the nearest 27 grid cells,
 358      * thus we loop over -1..1 in X,Y and Z direction
 359      * We apply the TSC (triangle shaped charge)
 360      * see Luty et. al, JCP 103 (1995) 3014
 361      */
 362
 363     calc_weights(i,nx,ny,nz,x[i],box,invh,ixyz,WXYZ);
 364
 365     for(ll=llim2; (ll<=ulim2); ll++) {
 366       ixw[ll-llim2] = nnx[ixyz[XX]+ll+nx];
 367       iyw[ll-llim2] = nny[ixyz[YY]+ll+ny];
 368       izw[ll-llim2] = nnz[ixyz[ZZ]+ll+nz];
 369     }
 370
 371     qi      = charge[i];
 372     pi      = gather_inner(JCXYZ,WXYZ,ixw,iyw,izw,la2,la12,
 373                            c1x,c1y,c1z,c2x,c2y,c2z,
 374                            qi,f[i],ptr);
 375
 376     energy += pi*qi;
 377     pot[i]  = pi;
 378   }
 379
 380   inc_nrnb(nrnb,eNR_GATHERF,27*nr);
 381   inc_nrnb(nrnb,eNR_WEIGHTS,3*nr);
 382
 383   return energy*0.5;
 384 }
 385
 386 static void convolution(FILE *fp,bool bVerbose,t_fftgrid *grid,real ***ghat,
 387                         t_commrec *cr)
 388 {
 389   int      i,j,k,index;
 390   real     gk;
 391   int      nx,ny,nz,nx2,ny2,nz2,la2,la12;
 392   t_complex  *ptr;
 393   int      *nTest;
 394   int jstart,jend;
 395
 396   unpack_fftgrid(grid,&nx,&ny,&nz,&nx2,&ny2,&nz2,
 397                  &la2,&la12,FALSE,(real **)&ptr);
 398   snew(nTest,grid->nptr);
 399
 400   if(PAR(cr)) {
 401 #if (defined GMX_MPI && !defined GMX_WITHOUT_FFTW)
 402     jstart=grid->pfft.local_y_start_after_transpose;
 403     jend=jstart+grid->pfft.local_ny_after_transpose;
 404
 405     for(j=jstart; (j<jend); j++) { /* local cells */
 406         for(i=0; (i<nx); i++) {
 407             for(k=0;k<(nz/2+1); k++) {
 408                 gk    = ghat[i][j][k];
 409                 index = INDEX(j,i,k);
 410                 ptr[index].re *= gk;
 411                 ptr[index].im *= gk;
 412                 nTest[index]++;
 413             }
 414         }
 415     }
 416 #ifdef DEBUG
 417     for(j=jstart; (j<jend); j++) {
 418         for(i=0; (i<nx); i++) {
 419             for(k=0; k<(nz/2+1); k++) {
 420                 index = INDEX(j,i,k);
 421                 if (nTest[index] != 1)
 422                     fprintf(fp,"Index %d sucks, set %d times\n",index,nTest[index]);
 423             }
 424         }
 425     }
 426 #endif /* DEBUG */
 427 #endif /* GMX_MPI */
 428   } else { /* if not running in parallel */
 429       for(i=0; (i<nx); i++) {
 430           for(j=0; (j<ny); j++) {
 431               for(k=0;k<(nz/2+1); k++) {
 432                   gk    = ghat[i][j][k];
 433                   index = INDEX(i,j,k);
 434                   ptr[index].re *= gk;
 435                   ptr[index].im *= gk;
 436                   nTest[index]++;
 437               }
 438           }
 439       }
 440 #ifdef DEBUG
 441       for(i=0; (i<nx); i++) {
 442           for(j=0; (j<ny); j++) {
 443               for(k=0; k<(nz/2+1); k++) {
 444                   index = INDEX(i,j,k);
 445                   if (nTest[index] != 1)
 446                       fprintf(fp,"Index %d sucks, set %d times\n",index,nTest[index]);
 447               }
 448           }
 449       }
 450 #endif
 451   }
 452   sfree(nTest);
 453 }
 454
 455
 456 void solve_pppm(FILE *fp,t_commrec *cr,
 457                 t_fftgrid *grid,real ***ghat,rvec box,
 458                 bool bVerbose,t_nrnb *nrnb)
 459 {
 460   int  ntot,npppm;
 461
 462   /*  if (bVerbose)
 463       print_fftgrid(fp,"Q-Real",grid,grid->nxyz,"qreal.pdb",box,TRUE);*/
 464
 465   gmxfft3D(grid,GMX_FFT_REAL_TO_COMPLEX,cr);
 466
 467   /*  if (bVerbose) {
 468       print_fftgrid(fp,"Q-k",grid,1.0,"qk-re.pdb",box,TRUE);
 469       print_fftgrid(fp,"Q-k",grid,1.0,"qk-im.pdb",box,FALSE);
 470       fprintf(stderr,"Doing convolution\n");
 471       }*/
 472
 473   convolution(fp,bVerbose,grid,ghat,cr);
 474
 475   /*  if (bVerbose)
 476       print_fftgrid(fp,"Convolution",grid,1.0,
 477       "convolute.pdb",box,TRUE);*/
 478
 479   gmxfft3D(grid,GMX_FFT_COMPLEX_TO_REAL,cr);
 480
 481   /*  if (bVerbose)
 482       print_fftgrid(fp,"Potential",grid,1.0,"pot.pdb",box,TRUE);*/
 483
 484   ntot  = grid->nxyz;
 485   npppm = ntot*log((real)ntot)/log(2.0);
 486   inc_nrnb(nrnb,eNR_FFT,2*npppm);
 487   inc_nrnb(nrnb,eNR_CONV,ntot);
 488 }
 489
 490
 491 static rvec      beta;
 492 static real      ***ghat=NULL;
 493 static t_fftgrid *grid=NULL;
 494
 495 #endif
 496
 497
 498 int gmx_pppm_init(FILE *log,      t_commrec *cr,
 499                   const output_env_t oenv, bool bVerbose,
 500                   bool bOld,      matrix box,
 501                   char *ghatfn,   t_inputrec *ir,
 502                   bool bReproducible)
 503 {
 504   int   nx,ny,nz,m,porder;
 505   ivec  grids;
 506   real  r1,rc;
 507   const real tol = 1e-5;
 508   rvec  box_diag,spacing;
 509
 510 #ifdef DISABLE_PPPM
 511     gmx_fatal(FARGS,"PPPM temporarily disabled while working on 2DPME\n");
 512     return -1;
 513 #else
 514
 515 #ifdef GMX_WITHOUT_FFTW
 516   gmx_fatal(FARGS,"PPPM used, but GROMACS was compiled without FFTW support!\n");
 517 #endif
 518
 519   if (log) {
 520     if (cr != NULL) {
 521       if (cr->nnodes > 1)
 522         fprintf(log,"Initializing parallel PPPM.\n");
 523     }
 524     fprintf(log,
 525             "Will use the PPPM algorithm for long-range electrostatics\n");
 526   }
 527
 528   for(m=0; m<DIM; m++)
 529     box_diag[m] = box[m][m];
 530
 531   if (!gmx_fexist(ghatfn)) {
 532     beta[XX]=beta[YY]=beta[ZZ]= 1.85;
 533     nx     = ir->nkx;
 534     ny     = ir->nky;
 535     nz     = ir->nkz;
 536
 537     if (log) {
 538       fprintf(log,"Generating Ghat function\n");
 539       fprintf(log,"Grid size is %d x %d x %d\n",nx,ny,nz);
 540     }
 541
 542     if ((nx < 4) || (ny < 4) || (nz < 4))
 543       gmx_fatal(FARGS,"Grid must be at least 4 points in all directions");
 544
 545     ghat   = mk_rgrid(nx,ny,nz);
 546     mk_ghat(NULL,nx,ny,nz,ghat,box_diag,
 547             ir->rcoulomb_switch,ir->rcoulomb,TRUE,bOld);
 548
 549     if (bVerbose)
 550       pr_scalar_gk("generghat.xvg",oenv,nx,ny,nz,box_diag,ghat);
 551   }
 552   else {
 553     fprintf(stderr,"Reading Ghat function from %s\n",ghatfn);
 554     ghat = rd_ghat(log,oenv,ghatfn,grids,spacing,beta,&porder,&r1,&rc);
 555
 556     /* Check whether cut-offs correspond */
 557     if ((fabs(r1-ir->rcoulomb_switch)>tol) || (fabs(rc-ir->rcoulomb)>tol)) {
 558       if (log) {
 559         fprintf(log,"rcoulomb_switch = %10.3e  rcoulomb = %10.3e"
 560                 "  r1 = %10.3e  rc = %10.3e\n",
 561                 ir->rcoulomb_switch,ir->rcoulomb,r1,rc);
 562         fflush(log);
 563       }
 564       gmx_fatal(FARGS,"Cut-off lengths in tpb file and Ghat file %s "
 565                   "do not match\nCheck your log file!",ghatfn);
 566     }
 567
 568     /* Check whether boxes correspond */
 569     for(m=0; (m<DIM); m++)
 570       if (fabs(box_diag[m]-grids[m]*spacing[m]) > tol) {
 571         if (log) {
 572           pr_rvec(log,0,"box",box_diag,DIM,TRUE);
 573           pr_rvec(log,0,"grid-spacing",spacing,DIM,TRUE);
 574           pr_ivec(log,0,"grid size",grids,DIM,TRUE);
 575           fflush(log);
 576         }
 577         gmx_fatal(FARGS,"Box sizes in tpb file and Ghat file %s do not match\n"
 578                     "Check your log file!",ghatfn);
 579       }
 580
 581     if (porder != 2)
 582       gmx_fatal(FARGS,"porder = %d, should be 2 in %s",porder,ghatfn);
 583
 584     nx = grids[XX];
 585     ny = grids[YY];
 586     nz = grids[ZZ];
 587
 588     if (bVerbose)
 589       pr_scalar_gk("optimghat.xvg",oenv,nx,ny,nz,box_diag,ghat);
 590   }
 591   /* Now setup the FFT things */
 592 #ifdef GMX_MPI
 593   cr->mpi_comm_mygroup=cr->mpi_comm_mysim;
 594 #endif
 595   grid = mk_fftgrid(nx,ny,nz,NULL,NULL,cr,bReproducible);
 596
 597   return 0;
 598 #endif
 599 }
 600
 601 int gmx_pppm_do(FILE *log,       gmx_pme_t pme,
 602                 bool bVerbose,
 603                 rvec x[],        rvec f[],
 604                 real charge[],   rvec box,
 605                 real phi[],      t_commrec *cr,
 606                 int start,       int nr,
 607                 t_nrnb *nrnb,
 608                 int pme_order,   real *energy)
 609 {
 610 #ifdef DISABLE_PPPM
 611     gmx_fatal(FARGS,"PPPM temporarily disabled while working on 2DPME\n");
 612     return -1;
 613 #else
 614
 615     /* Make the grid empty */
 616   clear_fftgrid(grid);
 617
 618   /* First step: spreading the charges over the grid. */
 619   spread_q(log,bVerbose,start,nr,x,charge,box,grid,nrnb);
 620
 621   /* In the parallel code we have to sum the grids from neighbouring nodes */
 622   if (PAR(cr))
 623     gmx_sum_qgrid(pme,cr,grid,GMX_SUM_QGRID_FORWARD);
 624
 625   /* Second step: solving the poisson equation in Fourier space */
 626   solve_pppm(log,cr,grid,ghat,box,bVerbose,nrnb);
 627
 628   /* In the parallel code we have to sum once again... */
 629   if (PAR(cr))
 630     gmx_sum_qgrid(pme,cr,grid,GMX_SUM_QGRID_BACKWARD);
 631
 632   /* Third and last step: gather the forces, energies and potential
 633    * from the grid.
 634    */
 635   *energy = gather_f(log,bVerbose,start,nr,x,f,charge,box,
 636                      phi,grid,beta,nrnb);
 637
 638   return 0;
 639 #endif
 640 }
 641
 642 #ifndef DISABLE_PPPM
 643 static int gmx_pppm_opt_do(FILE *log,       gmx_pme_t pme,
 644                            t_inputrec *ir,  bool bVerbose,
 645                            int natoms,
 646                            rvec x[],        rvec f[],
 647                            real charge[],   rvec box,
 648                            real phi[],      t_commrec *cr,
 649                            t_nrnb *nrnb,    rvec beta,
 650                            t_fftgrid *grid, bool bOld,
 651                            real *energy)
 652 {
 653   real      ***ghat;
 654   int       nx,ny,nz;
 655
 656   if (log)
 657     fprintf(log,"Generating Ghat function\n");
 658   nx     = ir->nkx;
 659   ny     = ir->nky;
 660   nz     = ir->nkz;
 661   ghat   = mk_rgrid(nx,ny,nz);
 662   mk_ghat(NULL,nx,ny,nz,ghat,box,ir->rcoulomb_switch,ir->rcoulomb,TRUE,bOld);
 663
 664   /* pr_scalar_gk("generghat.xvg",nx,ny,nz,box,ghat); */
 665
 666   /* Now start the actual PPPM procedure.
 667    * First step: spreading the charges over the grid.
 668    */
 669   /* Make the grid empty */
 670   clear_fftgrid(grid);
 671
 672   spread_q(log,bVerbose,0,natoms,x,charge,box,grid,nrnb);
 673
 674   /* Second step: solving the poisson equation in Fourier space */
 675   solve_pppm(log,cr,grid,ghat,box,bVerbose,nrnb);
 676
 677   /* Third and last step: gather the forces, energies and potential
 678    * from the grid.
 679    */
 680   *energy = gather_f(log,bVerbose,0,natoms,x,f,charge,box,phi,grid,beta,nrnb);
 681
 682   free_rgrid(ghat,nx,ny);
 683
 684   return 0;
 685 }
 686
 687 #endif