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