src/contrib/testlr.c

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