src/contrib/testlr.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 #include <math.h>
  36 #include <string.h>
  37 #include "typedefs.h"
  38 #include "gromacs/math/vec.h"
  39 #include "gromacs/math/units.h"
  40 #include "macros.h"
  41 #include "names.h"
  42 #include "gromacs/utility/smalloc.h"
  43 #include "gromacs/fileio/tpxio.h"
  44 #include "gromacs/commandline/pargs.h"
  45 #include "gromacs/fileio/writeps.h"
  46 #include "copyrite.h"
  47 #include "pppm.h"
  48 #include "readinp.h"
  49 #include "force.h"
  50 #include "nrnb.h"
  51 #include "coulomb.h"
  52 #include "gromacs/pbcutil/mshift.h"
  53 #include "poisson.h"
  54 #include "mdatoms.h"
  55
  56 static real phi_sr(FILE *log,int nj,rvec x[],real charge[],real rc,real r1,
  57                    rvec box, real phi[],t_block *excl,rvec f_sr[],gmx_bool bOld)
  58 {
  59   int  i,j,k,m,ni,i1,i2;
  60   real pp,r2,R,R_1,R_2,rc2;
  61   real qi,qj,vsr,eps,fscal;
  62   rvec dx;
  63
  64   vsr = 0.0;
  65   eps = ONE_4PI_EPS0;
  66   rc2 = rc*rc;
  67   ni=0;
  68   for(i=0; (i<nj-1); i++) {
  69     qi=charge[i];
  70     for(j=i+1; (j<nj); j++) {
  71       i1=excl->index[i];
  72       i2=excl->index[i+1];
  73       for(k=i1; (k<i2); k++)
  74         if (excl->a[k] == j)
  75           break;
  76       if (k == i2) {
  77         r2=calc_dx2dx(x[i],x[j],box,dx);
  78         if (r2 < rc2) {
  79           qj    = charge[j];
  80           R_1   = invsqrt(r2);
  81           R_2   = R_1*R_1;
  82           R     = invsqrt(R_2);
  83           if (bOld) {
  84             fscal = old_f(R,rc,r1)*R_2;
  85             pp    = old_phi(R,rc,r1);
  86           }
  87           else {
  88             fscal = new_f(R,rc)*R_2;
  89             pp    = new_phi(R,rc);
  90           }
  91           phi[i] += eps*qj*pp;
  92           phi[j] += eps*qi*pp;
  93           vsr    += eps*qj*qi*pp;
  94           for(m=0; (m<DIM); m++) {
  95             f_sr[i][m] += dx[m]*fscal;
  96             f_sr[j][m] -= dx[m]*fscal;
  97           }
  98           ni++;
  99         }
 100       }
 101     }
 102   }
 103   fprintf(log,"There were %d short range interactions, vsr=%g\n",ni,vsr);
 104
 105   return vsr;
 106 }
 107
 108 void calc_ener(FILE *fp,char *title,gmx_bool bHeader,int nmol,
 109                int natoms,real phi[],real charge[],t_block *excl)
 110 {
 111   int  i,i1,i2,j,k;
 112   real qq,qi,vv,V,Vex,Vc,Vt;
 113
 114   qq   = 0;
 115   V    = 0;
 116   for(i=0; (i<natoms); i++) {
 117     vv   = charge[i]*phi[i];
 118     V   += vv;
 119     qq  += charge[i]*charge[i];
 120   }
 121   V  = 0.5*V;
 122   Vc = 0.5*C*ONE_4PI_EPS0*qq;
 123
 124   qq = 0;
 125   for(i=0; (i<excl->nr); i++) {
 126     i1 = excl->index[i];
 127     i2 = excl->index[i+1];
 128     qi = charge[i];
 129     for(j=i1; (j<i2); j++) {
 130       k = excl->a[j];
 131       if (k != i)
 132         qq+=qi*charge[k];
 133     }
 134   }
 135   Vex = qq*0.5*C*ONE_4PI_EPS0;
 136
 137   Vt = V - Vc - Vex;
 138
 139   if (bHeader) {
 140     fprintf(fp,"%12s  %12s  %12s  %12s  %12s  %12s\n",
 141             "","Vphi","Vself","Vexcl","Vtot","1Mol");
 142
 143   }
 144   fprintf(fp,"%12s  %12.5e  %12.5e  %12.5e  %12.5e  %12.5e\n",
 145           title,V,Vc,Vex,Vt,Vt/nmol);
 146 }
 147
 148 void write_pqr(char *fn,t_atoms *atoms,rvec x[],real phi[],real dx)
 149 {
 150   FILE *fp;
 151   int  i,rnr;
 152
 153   fp=gmx_fio_fopen(fn,"w");
 154   for(i=0; (i<atoms->nr); i++) {
 155     rnr=atoms->atom[i].resnr;
 156     fprintf(fp,"%-6s%5d  %-4.4s%3.3s %c%4d    %8.3f%8.3f%8.3f%6.2f%6.2f\n",
 157             "ATOM",(i+1),*atoms->atomname[i],*atoms->resname[rnr],' ',
 158             (rnr+1) % 10000,
 159             10*(dx+x[i][XX]),10*x[i][YY],10*(x[i][ZZ]),0.0,phi[i]);
 160   }
 161   gmx_fio_fclose(fp);
 162 }
 163
 164 void write_grid_pqr(char *fn,int nx,int ny,int nz,real ***phi)
 165 {
 166   FILE *fp;
 167   int  i,j,k,rnr=0;
 168   real fac=4.0;
 169
 170   fp=gmx_fio_fopen(fn,"w");
 171   for(i=0; (i<nx); i++)
 172     for(j=0; (j<ny); j++)
 173       for(k=0; (k<nz); k++,rnr++)
 174         fprintf(fp,"%-6s%5d  %-4.4s%3.3s %c%4d    %8.3f%8.3f%8.3f%6.2f%6.2f\n",
 175                 "ATOM",(i+1),"C","C",' ',
 176                 1+((rnr+1) % 10000),fac*i,fac*j,fac*k,0.0,phi[i][j][k]);
 177   gmx_fio_fclose(fp);
 178 }
 179 void plot_phi(char *fn,rvec box,int natoms,rvec x[],real phi[])
 180 {
 181   t_psdata eps;
 182   real phi_max,rr,gg,bb,fac,dx,x0,y0;
 183   real offset;
 184   int  i;
 185
 186   phi_max=phi[0];
 187   rr=gg=bb=0.0;
 188   for(i=0; (i<natoms); i++)
 189     phi_max=max(phi_max,fabs(phi[i]));
 190
 191   if (phi_max==0.0) {
 192     fprintf(stderr,"All values zero, see .out file\n");
 193     return;
 194   }
 195   offset=20.0;
 196   fac=15.0;
 197 #ifdef DEBUG
 198   fprintf(stderr,"Scaling box by %g\n",fac);
 199 #endif
 200   eps=ps_open(fn,0,0,
 201               (real)(fac*box[XX]+2*offset),(real)(fac*box[YY]+2*offset));
 202   ps_translate(eps,offset,offset);
 203   ps_color(eps,0,0,0);
 204   ps_box(eps,1,1,(real)(fac*box[XX]-1),(real)(fac*box[YY]-1));
 205   dx=0.15*fac;
 206   for(i=0; (i<natoms); i++) {
 207     rr=gg=bb=1.0;
 208     if (phi[i] < 0)
 209       gg=bb=(1.0+(phi[i]/phi_max));
 210     else
 211       rr=gg=(1.0-(phi[i]/phi_max));
 212     rr=rgbset(rr);
 213     gg=rgbset(gg);
 214     bb=rgbset(bb);
 215     ps_color(eps,rr,gg,bb);
 216     x0=fac*x[i][XX];
 217     y0=fac*x[i][YY];
 218     ps_fillbox(eps,(real)(x0-dx),(real)(y0-dx),(real)(x0+dx),(real)(y0+dx));
 219   }
 220   ps_close(eps);
 221 }
 222
 223 void plot_qtab(char *fn,int nx,int ny,int nz,real ***qtab)
 224 {
 225   rvec box;
 226   rvec *xx;
 227   real *phi;
 228   int  i,npt,ix,iy,iz;
 229
 230   box[XX]=nx;
 231   box[YY]=ny;
 232   box[ZZ]=nz;
 233
 234   npt=(box[XX]*box[YY]*box[ZZ]);
 235   snew(xx,npt);
 236   snew(phi,npt);
 237   nx/=2;
 238   ny/=2;
 239   nz/=2;
 240   i=0;
 241   for(ix=-nx; (ix<nx); ix++)
 242     for(iy=-ny; (iy<ny); iy++)
 243       for(iz=-nz; (iz<nz); iz++,i++) {
 244         xx[i][XX]=ix+nx+0.5;
 245         xx[i][YY]=iy+ny+0.5;
 246         xx[i][ZZ]=iz+nz+0.5; /* onzin */
 247         phi[i]=qtab[ix+nx][iy+ny][iz+nz];
 248       }
 249
 250   plot_phi(fn,box,npt,xx,phi);
 251
 252   sfree(xx);
 253   sfree(phi);
 254 }
 255
 256 void print_phi(char *fn,int natoms,rvec x[],real phi[])
 257 {
 258   FILE *fp;
 259   int  i;
 260
 261   fp=gmx_fio_fopen(fn,"w");
 262   for(i=0; (i<natoms); i++)
 263     fprintf(fp,"%10d  %12.5e\n",i,phi[i]);
 264   gmx_fio_fclose(fp);
 265 }
 266
 267 void pr_f(char *fn,int natoms,rvec f[])
 268 {
 269   FILE *fp;
 270   int  i;
 271
 272   fp=gmx_fio_fopen(fn,"w");
 273   for(i=0; (i<natoms); i++)
 274     fprintf(fp,"  %12.5e\n  %12.5e\n  %12.5e\n",f[i][XX],f[i][YY],f[i][ZZ]);
 275   gmx_fio_fclose(fp);
 276 }
 277
 278 void test_pppm(FILE *log,       gmx_bool bVerbose,
 279                gmx_bool bGenerGhat, char *ghatfn,
 280                t_atoms *atoms,  t_inputrec *ir,
 281                rvec x[],        rvec f[],
 282                real charge[],   rvec box,
 283                real phi[],      real phi_s[],
 284                int nmol,        t_commrec *cr,
 285                gmx_bool bOld,       t_block *cgs)
 286 {
 287   char       buf[256];
 288   real       ener;
 289   int        i;
 290   t_nrnb     nrnb;
 291   t_nsborder nsb;
 292
 293   init_nrnb(&nrnb);
 294   calc_nsb(cgs,1,&nsb,0);
 295
 296   /* First time only setup is done! */
 297   init_pppm(log,cr,&nsb,bVerbose,bOld,box,ghatfn,ir);
 298
 299   ener = do_pppm(log,bVerbose,x,f,charge,box,phi,cr,&nsb,&nrnb);
 300   fprintf(log,"Vpppm = %g\n",ener);
 301
 302   sprintf(buf,"PPPM-%d.pdb",ir->nkx);
 303   write_pqr(buf,atoms,x,phi,0);
 304
 305   pr_f("pppm-force",atoms->nr,f);
 306
 307   calc_ener(log,buf,FALSE,nmol,atoms->nr,phi,charge,&atoms->excl);
 308
 309   for(i=0; (i<atoms->nr); i++)
 310     phi[i]+=phi_s[i];
 311   sprintf(buf,"PPPM-%d+SR",ir->nkx);
 312   calc_ener(log,buf,FALSE,nmol,atoms->nr,phi,charge,&atoms->excl);
 313   strcat(buf,".pdb");
 314   write_pqr(buf,atoms,x,phi,0);
 315 }
 316
 317 void test_poisson(FILE *log,       gmx_bool bVerbose,
 318                   t_atoms *atoms,  t_inputrec *ir,
 319                   rvec x[],        rvec f[],
 320                   real charge[],   rvec box,
 321                   real phi[],      real phi_s[],
 322                   int nmol,        t_commrec *cr,
 323                   gmx_bool bFour,      rvec f_four[],
 324                   real phi_f[],    gmx_bool bOld)
 325 {
 326   char buf[256];
 327   real ener;
 328   rvec beta;
 329   int  i,nit;
 330   t_nrnb nrnb;
 331
 332   init_nrnb(&nrnb);
 333
 334   /* First time only setup is done! */
 335   if (bFour) {
 336     for(i=0; (i<atoms->nr); i++)
 337       phi_f[i] -= phi_s[i];
 338     ener = do_optimize_poisson(log,bVerbose,ir,atoms->nr,x,f,charge,
 339                                box,phi,cr,&nrnb,f_four,phi_f,beta,bOld);
 340     for(i=0; (i<atoms->nr); i++)
 341       phi_f[i] += phi_s[i];
 342     nit = 0;
 343   }
 344   else {
 345     ener = do_poisson(log,bVerbose,ir,atoms->nr,x,f,charge,box,phi,
 346                       cr,&nrnb,&nit,bOld);
 347   }
 348
 349   fprintf(log,"Vpoisson = %g, nit = %d\n",ener,nit);
 350
 351   sprintf(buf,"POISSON-%d.pdb",ir->nkx);
 352   write_pqr(buf,atoms,x,phi,0);
 353
 354   pr_f("poisson-force",atoms->nr,f);
 355
 356   calc_ener(log,buf,FALSE,nmol,atoms->nr,phi,charge,&atoms->excl);
 357
 358   for(i=0; (i<atoms->nr); i++)
 359     phi[i]+=phi_s[i];
 360   sprintf(buf,"POISSON-%d+SR",ir->nkx);
 361   calc_ener(log,buf,FALSE,nmol,atoms->nr,phi,charge,&atoms->excl);
 362   strcat(buf,".pdb");
 363   write_pqr(buf,atoms,x,phi,0);
 364 }
 365
 366 void test_four(FILE *log,int NFILE,t_filenm fnm[],t_atoms *atoms,
 367                t_inputrec *ir,rvec x[],rvec f[],rvec box,real charge[],
 368                real phi_f[],real phi_s[],int nmol,t_commrec *cr,
 369                gmx_bool bOld,gmx_bool bOldEwald)
 370 {
 371   int  i;
 372   real energy;
 373   ewald_tab_t et;
 374
 375   init_ewald_tab(&et, NULL, ir, log);
 376
 377   if (bOldEwald)
 378     energy = do_ewald(log,ir,atoms->nr,x,f,charge,box,phi_f,cr,bOld,et);
 379   else
 380     energy = do_ewald_new(log,ir,atoms->nr,x,f,charge,box,phi_f,cr,bOld,et);
 381
 382   /*symmetrize_phi(log,atoms->nr,phi_f,bVerbose);*/
 383
 384   /* Plot the long range fourier solution in a matrix */
 385   plot_phi(opt2fn("-elf",NFILE,fnm),box,atoms->nr,x,phi_f);
 386   print_phi(opt2fn("-of",NFILE,fnm),atoms->nr,x,phi_f);
 387   calc_ener(log,"Fourier",FALSE,nmol,atoms->nr,phi_f,charge,&atoms->excl);
 388   write_pqr(opt2fn("-pf",NFILE,fnm),atoms,x,phi_f,0.0/*1.5*box[XX]*/);
 389   pr_f("four-force",atoms->nr,f);
 390
 391   /* Calc and plot the total potential */
 392   for(i=0; (i<atoms->nr); i++) {
 393     phi_f[i]+=phi_s[i];
 394     /*clear_rvec(f[i]);*/
 395   }
 396   calc_ener(log,"Fourier+SR",FALSE,nmol,atoms->nr,phi_f,charge,&atoms->excl);
 397 }
 398
 399 static void print_opts(FILE *fp,t_inputrec *ir,gmx_bool bFour)
 400 {
 401   fprintf(fp,"Ewald solution: %s\n",gmx_bool_names[bFour]);
 402   fprintf(fp,"r1:       %10.3e\n",ir->rcoulomb_switch);
 403   fprintf(fp,"rc:       %10.3e\n",ir->rcoulomb);
 404   if (bFour)
 405     fprintf(fp,"KVectors: %10d  %10d  %10d\n",ir->nkx,ir->nky,ir->nkz);
 406   fprintf(fp,"\n");
 407 }
 408
 409 int main(int argc,char *argv[])
 410 {
 411   static char *desc[] = {
 412     "[TT]testlr[tt] tests the PPPM and Ewald method for the",
 413     "long range electrostatics problem."
 414   };
 415   static t_filenm  fnm[] = {
 416     { efTPX, NULL,   NULL,       ffREAD },
 417     { efHAT, "-g",   "ghat",     ffOPTRD },
 418     { efOUT, "-o",   "rho",      ffOPTWR },
 419     { efOUT, "-op",  "lr-pb",    ffOPTWR },
 420     { efOUT, "-of",  "lr-four",  ffOPTWR },
 421     { efOUT, "-opt", "tot-pb",   ffOPTWR },
 422     { efOUT, "-oft", "tot-four", ffOPTWR },
 423     { efOUT, "-fin", "lr-four",  ffOPTWR },
 424     { efEPS, "-es",  "sr",       ffOPTWR },
 425     { efEPS, "-elf", "lr-four",  ffOPTWR },
 426     { efEPS, "-etf", "tot-four", ffOPTWR },
 427     { efEPS, "-qr",  "qk-real",  ffOPTWR },
 428     { efEPS, "-qi",  "qk-im",    ffOPTWR },
 429     { efEPS, "-elp", "lr-pb",    ffOPTWR },
 430     { efEPS, "-etp", "tot-pb",   ffOPTWR },
 431     { efEPS, "-rho", "rho",      ffOPTWR },
 432     { efEPS, "-qq",  "charge",   ffOPTWR },
 433     { efXVG, "-gt",  "gk-tab",   ffOPTWR },
 434     { efXVG, "-fcorr","fcorr",   ffWRITE },
 435     { efXVG, "-pcorr","pcorr",   ffWRITE },
 436     { efXVG, "-ftotcorr","ftotcorr",   ffWRITE },
 437     { efXVG, "-ptotcorr","ptotcorr",   ffWRITE },
 438     { efLOG, "-l",   "fptest",   ffWRITE },
 439     { efXVG, "-gr",  "spread",   ffOPTWR },
 440     { efPDB, "-pf",  "pqr-four", ffOPTWR },
 441     { efPDB, "-phitot", "pppm-phitot", ffOPTWR }
 442   };
 443 #define NFILE asize(fnm)
 444   FILE         *log;
 445   t_topology   top;
 446   t_tpxheader  stath;
 447   t_inputrec   ir;
 448   t_block      *excl;
 449   t_forcerec   *fr;
 450   t_commrec    *cr;
 451   t_mdatoms    *mdatoms;
 452   t_graph      *graph;
 453   int          i,step,nre,natoms,nmol;
 454   rvec         *x,*f_sr,*f_excl,*f_four,*f_pppm,*f_pois,box_size,hbox;
 455   matrix       box;
 456   real         t,lambda,vsr,*charge,*phi_f,*phi_pois,*phi_s,*phi_p3m,*rho;
 457   output_env_t oenv;
 458
 459   static gmx_bool bFour=FALSE,bVerbose=FALSE,bGGhat=FALSE,bPPPM=TRUE,
 460     bPoisson=FALSE,bOld=FALSE,bOldEwald=TRUE;
 461   static int nprocs = 1;
 462   static t_pargs pa[] = {
 463     { "-np",     FALSE, etINT,  &nprocs,  "Do it in parallel" },
 464     { "-ewald",  FALSE, etBOOL, &bFour,   "Do an Ewald solution"},
 465     { "-pppm",   FALSE, etBOOL, &bPPPM,   "Do a PPPM solution" },
 466     { "-poisson",FALSE, etBOOL, &bPoisson,"Do a Poisson solution" },
 467     {    "-v",   FALSE, etBOOL, &bVerbose,"Verbose on"},
 468     { "-ghat",   FALSE, etBOOL, &bGGhat,  "Generate Ghat function"},
 469     { "-old",    FALSE, etBOOL, &bOld,    "Use old function types"},
 470     { "-oldewald",FALSE,etBOOL, &bOldEwald,"Use old Ewald code"}
 471   };
 472
 473   CopyRight(stderr,argv[0]);
 474   parse_common_args_r(&argc,argv,PCA_CAN_TIME | PCA_CAN_VIEW,
 475                       NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
 476
 477   if (nprocs > 1) {
 478     cr = init_par(&argc,argv);
 479     open_log(ftp2fn(efLOG,NFILE,fnm),cr);
 480     log = stdlog;
 481   }
 482   else {
 483     cr     = init_par(&argc,argv);
 484     log    = ftp2FILE(efLOG,NFILE,fnm,"w");
 485     stdlog = log;  }
 486
 487
 488   /* Read topology and coordinates */
 489   read_tpxheader(ftp2fn(efTPX,NFILE,fnm),&stath,FALSE);
 490   snew(x,stath.natoms);
 491   snew(f_sr,stath.natoms);
 492   snew(f_excl,stath.natoms);
 493   snew(f_four,stath.natoms);
 494   snew(f_pppm,stath.natoms);
 495   snew(f_pois,stath.natoms);
 496   read_tpx(ftp2fn(efTPX,NFILE,fnm),&step,&t,&lambda,&ir,
 497            box,&natoms,x,NULL,NULL,&top);
 498   excl=&(top.atoms.excl);
 499   nmol=top.blocks[ebMOLS].nr;
 500
 501   /* Allocate space for potential, charges and rho (charge density) */
 502   snew(charge,stath.natoms);
 503   snew(phi_f,stath.natoms);
 504   snew(phi_p3m,stath.natoms);
 505   snew(phi_pois,stath.natoms);
 506   snew(phi_s,stath.natoms);
 507   snew(rho,stath.natoms);
 508
 509   /* Set the charges */
 510   for(i=0; (i<natoms); i++)
 511     charge[i]=top.atoms.atom[i].q;
 512
 513   /* Make a simple box vector instead of tensor */
 514   for(i=0; (i<DIM); i++)
 515     box_size[i]=box[i][i];
 516
 517   /* Set some constants */
 518   fr      = mk_forcerec();
 519   mdatoms = atoms2md(&(top.atoms),FALSE,FALSE);
 520
 521   set_LRconsts(log,ir.rcoulomb_switch,ir.rcoulomb,box_size,fr);
 522   init_forcerec(log,fr,&ir,&(top.blocks[ebMOLS]),cr,
 523                 &(top.blocks[ebCGS]),&(top.idef),mdatoms,box,FALSE);
 524   calc_shifts(box,box_size,fr->shift_vec,FALSE);
 525
 526   /* Periodicity stuff */
 527   graph = mk_graph(&(top.idef),top.atoms.nr,FALSE,FALSE);
 528   shift_self(graph,fr->shift_vec,x);
 529
 530   calc_LRcorrections(log,0,natoms,ir.rcoulomb_switch,
 531                      ir.rcoulomb,charge,excl,x,f_excl,bOld);
 532   pr_f("f_excl.dat",natoms,f_excl);
 533
 534   /* Compute the short range potential */
 535   put_atoms_in_box(natoms,box,x);
 536   vsr=phi_sr(log,natoms,x,charge,ir.rcoulomb,
 537              ir.rcoulomb_switch,box_size,phi_s,excl,f_sr,bOld);
 538   pr_f("f_sr.dat",natoms,f_sr);
 539
 540   /* Plot the short range potential in a matrix */
 541   calc_ener(log,"Short Range",TRUE,nmol,natoms,phi_s,charge,excl);
 542
 543
 544   if (bFour)
 545     test_four(log,NFILE,fnm,&(top.atoms),&ir,x,f_four,box_size,charge,phi_f,
 546               phi_s,nmol,cr,bOld,bOldEwald);
 547
 548   if (bPPPM)
 549     test_pppm(log,bVerbose,bGGhat,opt2fn("-g",NFILE,fnm),
 550               &(top.atoms),&ir,x,f_pppm,charge,box_size,phi_p3m,phi_s,nmol,
 551               cr,bOld,&(top.blocks[ebCGS]));
 552
 553   if (bPoisson)
 554     test_poisson(log,bVerbose,
 555                  &(top.atoms),&ir,x,f_pois,charge,box_size,phi_pois,
 556                  phi_s,nmol,cr,bFour,f_four,phi_f,bOld);
 557
 558   if (bPPPM && bFour)
 559     analyse_diff(log,"PPPM",oenv,
 560                  top.atoms.nr,f_four,f_pppm,phi_f,phi_p3m,phi_s,
 561                  opt2fn("-fcorr",NFILE,fnm),
 562                  opt2fn("-pcorr",NFILE,fnm),
 563                  opt2fn("-ftotcorr",NFILE,fnm),
 564                  opt2fn("-ptotcorr",NFILE,fnm));
 565
 566   if (bPoisson && bFour)
 567     analyse_diff(log,"Poisson",oenv,
 568                  top.atoms.nr,f_four,f_pois,phi_f,phi_pois,phi_s,
 569                  opt2fn("-fcorr",NFILE,fnm),
 570                  opt2fn("-pcorr",NFILE,fnm),
 571                  opt2fn("-ftotcorr",NFILE,fnm),
 572                  opt2fn("-ptotcorr",NFILE,fnm));
 573
 574   gmx_fio_fclose(log);
 575
 576   gmx_thanx(stderr);
 577
 578   return 0;
 579 }
 580