src/gmxlib/testlr.c

   1 /*
   2  * $Id$
   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.1
  11  * Copyright (c) 1991-2001, University of Groningen, The Netherlands
  12  * This program is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU General Public License
  14  * as published by the Free Software Foundation; either version 2
  15  * of the License, or (at your option) any later version.
  16  *
  17  * If you want to redistribute modifications, please consider that
  18  * scientific software is very special. Version control is crucial -
  19  * bugs must be traceable. We will be happy to consider code for
  20  * inclusion in the official distribution, but derived work must not
  21  * be called official GROMACS. Details are found in the README & COPYING
  22  * files - if they are missing, get the official version at www.gromacs.org.
  23  *
  24  * To help us fund GROMACS development, we humbly ask that you cite
  25  * the papers on the package - you can find them in the top README file.
  26  *
  27  * For more info, check our website at http://www.gromacs.org
  28  *
  29  * And Hey:
  30  * Gnomes, ROck Monsters And Chili Sauce
  31  */
  32
  33 #include <math.h>
  34 #include <string.h>
  35 #include "typedefs.h"
  36 #include "vec.h"
  37 #include "physics.h"
  38 #include "macros.h"
  39 #include "names.h"
  40 #include "smalloc.h"
  41 #include "tpxio.h"
  42 #include "statutil.h"
  43 #include "writeps.h"
  44 #include "assert.h"
  45 #include "copyrite.h"
  46 #include "xvgr.h"
  47 #include "minvert.h"
  48 #include "pppm.h"
  49 #include "readinp.h"
  50 #include "main.h"
  51 #include "force.h"
  52 #include "nrnb.h"
  53 #include "shift_util.h"
  54 #include "ewald_util.h"
  55 #include "mshift.h"
  56 #include "poisson.h"
  57 #include "mdatoms.h"
  58
  59 void pr_f(char *fn,int natoms,rvec f[])
  60 {
  61   FILE *fp;
  62   int  i;
  63
  64   fp=ffopen(fn,"w");
  65   for(i=0; (i<natoms); i++)
  66     fprintf(fp,"  %12.5e\n  %12.5e\n  %12.5e\n",f[i][XX],f[i][YY],f[i][ZZ]);
  67   fclose(fp);
  68 }
  69
  70 void test_pppm(FILE *log,       bool bVerbose,
  71                bool bGenerGhat, char *ghatfn,
  72                t_atoms *atoms,  t_inputrec *ir,
  73                rvec x[],        rvec f[],
  74                real charge[],   rvec box,
  75                real phi[],      real phi_s[],
  76                int nmol,        t_commrec *cr,
  77                bool bOld,       t_block *cgs)
  78 {
  79   char       buf[256];
  80   real       ener;
  81   int        i;
  82   t_nrnb     nrnb;
  83   t_nsborder nsb;
  84
  85   init_nrnb(&nrnb);
  86   calc_nsb(cgs,1,&nsb,0);
  87
  88   /* First time only setup is done! */
  89   init_pppm(log,cr,&nsb,bVerbose,bOld,box,ghatfn,ir);
  90
  91   ener = do_pppm(log,bVerbose,x,f,charge,box,phi,cr,&nsb,&nrnb);
  92   fprintf(log,"Vpppm = %g\n",ener);
  93
  94   sprintf(buf,"PPPM-%d.pdb",ir->nkx);
  95   write_pqr(buf,atoms,x,phi,0);
  96
  97   pr_f("pppm-force",atoms->nr,f);
  98
  99   calc_ener(log,buf,FALSE,nmol,atoms->nr,phi,charge,&atoms->excl);
 100
 101   for(i=0; (i<atoms->nr); i++)
 102     phi[i]+=phi_s[i];
 103   sprintf(buf,"PPPM-%d+SR",ir->nkx);
 104   calc_ener(log,buf,FALSE,nmol,atoms->nr,phi,charge,&atoms->excl);
 105   strcat(buf,".pdb");
 106   write_pqr(buf,atoms,x,phi,0);
 107 }
 108
 109 void test_poisson(FILE *log,       bool bVerbose,
 110                   t_atoms *atoms,  t_inputrec *ir,
 111                   rvec x[],        rvec f[],
 112                   real charge[],   rvec box,
 113                   real phi[],      real phi_s[],
 114                   int nmol,        t_commrec *cr,
 115                   bool bFour,      rvec f_four[],
 116                   real phi_f[],    bool bOld)
 117 {
 118   char buf[256];
 119   real ener;
 120   rvec beta;
 121   int  i,nit;
 122   t_nrnb nrnb;
 123
 124   init_nrnb(&nrnb);
 125
 126   /* First time only setup is done! */
 127   if (bFour) {
 128     for(i=0; (i<atoms->nr); i++)
 129       phi_f[i] -= phi_s[i];
 130     ener = do_optimize_poisson(log,bVerbose,ir,atoms->nr,x,f,charge,
 131                                box,phi,cr,&nrnb,f_four,phi_f,beta,bOld);
 132     for(i=0; (i<atoms->nr); i++)
 133       phi_f[i] += phi_s[i];
 134     nit = 0;
 135   }
 136   else {
 137     ener = do_poisson(log,bVerbose,ir,atoms->nr,x,f,charge,box,phi,
 138                       cr,&nrnb,&nit,bOld);
 139   }
 140
 141   fprintf(log,"Vpoisson = %g, nit = %d\n",ener,nit);
 142
 143   sprintf(buf,"POISSON-%d.pdb",ir->nkx);
 144   write_pqr(buf,atoms,x,phi,0);
 145
 146   pr_f("poisson-force",atoms->nr,f);
 147
 148   calc_ener(log,buf,FALSE,nmol,atoms->nr,phi,charge,&atoms->excl);
 149
 150   for(i=0; (i<atoms->nr); i++)
 151     phi[i]+=phi_s[i];
 152   sprintf(buf,"POISSON-%d+SR",ir->nkx);
 153   calc_ener(log,buf,FALSE,nmol,atoms->nr,phi,charge,&atoms->excl);
 154   strcat(buf,".pdb");
 155   write_pqr(buf,atoms,x,phi,0);
 156 }
 157
 158 void test_four(FILE *log,int NFILE,t_filenm fnm[],t_atoms *atoms,
 159                t_inputrec *ir,rvec x[],rvec f[],rvec box,real charge[],
 160                real phi_f[],real phi_s[],int nmol,t_commrec *cr,
 161                bool bOld,bool bOldEwald)
 162 {
 163   int  i;
 164   real energy;
 165
 166   if (bOldEwald)
 167     energy = do_ewald(log,ir,atoms->nr,x,f,charge,box,phi_f,cr,bOld);
 168   else
 169     energy = do_ewald_new(log,ir,atoms->nr,x,f,charge,box,phi_f,cr,bOld);
 170
 171   /*symmetrize_phi(log,atoms->nr,phi_f,bVerbose);*/
 172
 173   /* Plot the long range fourier solution in a matrix */
 174   plot_phi(opt2fn("-elf",NFILE,fnm),box,atoms->nr,x,phi_f);
 175   print_phi(opt2fn("-of",NFILE,fnm),atoms->nr,x,phi_f);
 176   calc_ener(log,"Fourier",FALSE,nmol,atoms->nr,phi_f,charge,&atoms->excl);
 177   write_pqr(opt2fn("-pf",NFILE,fnm),atoms,x,phi_f,0.0/*1.5*box[XX]*/);
 178   pr_f("four-force",atoms->nr,f);
 179
 180   /* Calc and plot the total potential */
 181   for(i=0; (i<atoms->nr); i++) {
 182     phi_f[i]+=phi_s[i];
 183     /*clear_rvec(f[i]);*/
 184   }
 185   calc_ener(log,"Fourier+SR",FALSE,nmol,atoms->nr,phi_f,charge,&atoms->excl);
 186 }
 187
 188 static void print_opts(FILE *fp,t_inputrec *ir,bool bFour)
 189 {
 190   fprintf(fp,"Ewald solution: %s\n",bool_names[bFour]);
 191   fprintf(fp,"r1:       %10.3e\n",ir->rcoulomb_switch);
 192   fprintf(fp,"rc:       %10.3e\n",ir->rcoulomb);
 193   if (bFour)
 194     fprintf(fp,"KVectors: %10d  %10d  %10d\n",ir->nkx,ir->nky,ir->nkz);
 195   fprintf(fp,"\n");
 196 }
 197
 198 int main(int argc,char *argv[])
 199 {
 200   static char *desc[] = {
 201     "testlr tests the PPPM and Ewald method for the",
 202     "long range electrostatics problem."
 203   };
 204   static t_filenm  fnm[] = {
 205     { efTPX, NULL,   NULL,       ffREAD },
 206     { efHAT, "-g",   "ghat",     ffOPTRD },
 207     { efOUT, "-o",   "rho",      ffOPTWR },
 208     { efOUT, "-op",  "lr-pb",    ffOPTWR },
 209     { efOUT, "-of",  "lr-four",  ffOPTWR },
 210     { efOUT, "-opt", "tot-pb",   ffOPTWR },
 211     { efOUT, "-oft", "tot-four", ffOPTWR },
 212     { efOUT, "-fin", "lr-four",  ffOPTWR },
 213     { efEPS, "-es",  "sr",       ffOPTWR },
 214     { efEPS, "-elf", "lr-four",  ffOPTWR },
 215     { efEPS, "-etf", "tot-four", ffOPTWR },
 216     { efEPS, "-qr",  "qk-real",  ffOPTWR },
 217     { efEPS, "-qi",  "qk-im",    ffOPTWR },
 218     { efEPS, "-elp", "lr-pb",    ffOPTWR },
 219     { efEPS, "-etp", "tot-pb",   ffOPTWR },
 220     { efEPS, "-rho", "rho",      ffOPTWR },
 221     { efEPS, "-qq",  "charge",   ffOPTWR },
 222     { efXVG, "-gt",  "gk-tab",   ffOPTWR },
 223     { efXVG, "-fcorr","fcorr",   ffWRITE },
 224     { efXVG, "-pcorr","pcorr",   ffWRITE },
 225     { efXVG, "-ftotcorr","ftotcorr",   ffWRITE },
 226     { efXVG, "-ptotcorr","ptotcorr",   ffWRITE },
 227     { efLOG, "-l",   "fptest",   ffWRITE },
 228     { efXVG, "-gr",  "spread",   ffOPTWR },
 229     { efPDB, "-pf",  "pqr-four", ffOPTWR },
 230     { efPDB, "-phitot", "pppm-phitot", ffOPTWR }
 231   };
 232 #define NFILE asize(fnm)
 233   FILE         *log;
 234   t_topology   top;
 235   t_tpxheader  stath;
 236   t_inputrec   ir;
 237   t_block      *excl;
 238   t_forcerec   *fr;
 239   t_commrec    *cr;
 240   t_mdatoms    *mdatoms;
 241   t_graph      *graph;
 242   int          i,step,nre,natoms,nmol;
 243   rvec         *x,*f_sr,*f_excl,*f_four,*f_pppm,*f_pois,box_size,hbox;
 244   matrix       box;
 245   real         t,lambda,vsr,*charge,*phi_f,*phi_pois,*phi_s,*phi_p3m,*rho;
 246
 247   static bool bFour=FALSE,bVerbose=FALSE,bGGhat=FALSE,bPPPM=TRUE,
 248     bPoisson=FALSE,bOld=FALSE,bOldEwald=TRUE;
 249   static int nprocs = 1;
 250   static t_pargs pa[] = {
 251     { "-np",     FALSE, etINT,  &nprocs,  "Do it in parallel" },
 252     { "-ewald",  FALSE, etBOOL, &bFour,   "Do an Ewald solution"},
 253     { "-pppm",   FALSE, etBOOL, &bPPPM,   "Do a PPPM solution" },
 254     { "-poisson",FALSE, etBOOL, &bPoisson,"Do a Poisson solution" },
 255     {    "-v",   FALSE, etBOOL, &bVerbose,"Verbose on"},
 256     { "-ghat",   FALSE, etBOOL, &bGGhat,  "Generate Ghat function"},
 257     { "-old",    FALSE, etBOOL, &bOld,    "Use old function types"},
 258     { "-oldewald",FALSE,etBOOL, &bOldEwald,"Use old Ewald code"}
 259   };
 260
 261   CopyRight(stderr,argv[0]);
 262   parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_CAN_VIEW,TRUE,
 263                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL);
 264
 265   if (nprocs > 1) {
 266     cr = init_par(&argc,argv);
 267     open_log(ftp2fn(efLOG,NFILE,fnm),cr);
 268     log = stdlog;
 269   }
 270   else {
 271     cr     = init_par(&argc,argv);
 272     log    = ftp2FILE(efLOG,NFILE,fnm,"w");
 273     stdlog = log;  }
 274
 275
 276   /* Read topology and coordinates */
 277   read_tpxheader(ftp2fn(efTPX,NFILE,fnm),&stath,FALSE);
 278   snew(x,stath.natoms);
 279   snew(f_sr,stath.natoms);
 280   snew(f_excl,stath.natoms);
 281   snew(f_four,stath.natoms);
 282   snew(f_pppm,stath.natoms);
 283   snew(f_pois,stath.natoms);
 284   read_tpx(ftp2fn(efTPX,NFILE,fnm),&step,&t,&lambda,&ir,
 285            box,&natoms,x,NULL,NULL,&top);
 286   excl=&(top.atoms.excl);
 287   nmol=top.blocks[ebMOLS].nr;
 288
 289   /* Allocate space for potential, charges and rho (charge density) */
 290   snew(charge,stath.natoms);
 291   snew(phi_f,stath.natoms);
 292   snew(phi_p3m,stath.natoms);
 293   snew(phi_pois,stath.natoms);
 294   snew(phi_s,stath.natoms);
 295   snew(rho,stath.natoms);
 296
 297   /* Set the charges */
 298   for(i=0; (i<natoms); i++)
 299     charge[i]=top.atoms.atom[i].q;
 300
 301   /* Make a simple box vector instead of tensor */
 302   for(i=0; (i<DIM); i++)
 303     box_size[i]=box[i][i];
 304
 305   /* Set some constants */
 306   fr      = mk_forcerec();
 307   mdatoms = atoms2md(&(top.atoms),FALSE,FALSE);
 308
 309   set_LRconsts(log,ir.rcoulomb_switch,ir.rcoulomb,box_size,fr);
 310   init_forcerec(log,fr,&ir,&(top.blocks[ebMOLS]),cr,
 311                 &(top.blocks[ebCGS]),&(top.idef),mdatoms,box,FALSE);
 312   calc_shifts(box,box_size,fr->shift_vec,FALSE);
 313
 314   /* Periodicity stuff */
 315   graph = mk_graph(&(top.idef),top.atoms.nr,FALSE,FALSE);
 316   shift_self(graph,fr->shift_vec,x);
 317
 318   calc_LRcorrections(log,0,natoms,ir.rcoulomb_switch,
 319                      ir.rcoulomb,charge,excl,x,f_excl,bOld);
 320   pr_f("f_excl.dat",natoms,f_excl);
 321
 322   /* Compute the short range potential */
 323   put_atoms_in_box(natoms,box,x);
 324   vsr=phi_sr(log,natoms,x,charge,ir.rcoulomb,
 325              ir.rcoulomb_switch,box_size,phi_s,excl,f_sr,bOld);
 326   pr_f("f_sr.dat",natoms,f_sr);
 327
 328   /* Plot the short range potential in a matrix */
 329   calc_ener(log,"Short Range",TRUE,nmol,natoms,phi_s,charge,excl);
 330
 331
 332   if (bFour)
 333     test_four(log,NFILE,fnm,&(top.atoms),&ir,x,f_four,box_size,charge,phi_f,
 334               phi_s,nmol,cr,bOld,bOldEwald);
 335
 336   if (bPPPM)
 337     test_pppm(log,bVerbose,bGGhat,opt2fn("-g",NFILE,fnm),
 338               &(top.atoms),&ir,x,f_pppm,charge,box_size,phi_p3m,phi_s,nmol,
 339               cr,bOld,&(top.blocks[ebCGS]));
 340
 341   if (bPoisson)
 342     test_poisson(log,bVerbose,
 343                  &(top.atoms),&ir,x,f_pois,charge,box_size,phi_pois,
 344                  phi_s,nmol,cr,bFour,f_four,phi_f,bOld);
 345
 346   if (bPPPM && bFour)
 347     analyse_diff(log,"PPPM",
 348                  top.atoms.nr,f_four,f_pppm,phi_f,phi_p3m,phi_s,
 349                  opt2fn("-fcorr",NFILE,fnm),
 350                  opt2fn("-pcorr",NFILE,fnm),
 351                  opt2fn("-ftotcorr",NFILE,fnm),
 352                  opt2fn("-ptotcorr",NFILE,fnm));
 353
 354   if (bPoisson && bFour)
 355     analyse_diff(log,"Poisson",
 356                  top.atoms.nr,f_four,f_pois,phi_f,phi_pois,phi_s,
 357                  opt2fn("-fcorr",NFILE,fnm),
 358                  opt2fn("-pcorr",NFILE,fnm),
 359                  opt2fn("-ftotcorr",NFILE,fnm),
 360                  opt2fn("-ptotcorr",NFILE,fnm));
 361
 362   fclose(log);
 363
 364   thanx(stderr);
 365
 366   return 0;
 367 }
 368