src/contrib/pmetest.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-2006, 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 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include "typedefs.h"
  43 #include "macros.h"
  44 #include "smalloc.h"
  45 #include "copyrite.h"
  46 #include "main.h"
  47 #include "nrnb.h"
  48 #include "txtdump.h"
  49 #include "tpxio.h"
  50 #include "statutil.h"
  51 #include "futil.h"
  52 #include "gmx_fatal.h"
  53 #include "vec.h"
  54 #include "mdatoms.h"
  55 #include "coulomb.h"
  56 #include "nsb.h"
  57 #include "rmpbc.h"
  58 #include "pme.h"
  59 #include "force.h"
  60 #include "xvgr.h"
  61 #include "pbc.h"
  62
  63 #ifdef GMX_LIB_MPI
  64 #include <mpi.h>
  65 #endif
  66 #ifdef GMX_THREAD_MPI
  67 #include "tmpi.h"
  68 #endif
  69
  70 #include "block_tx.h"
  71
  72 rvec *xptr=NULL;
  73
  74 static int comp_xptr(const void *a,const void *b)
  75 {
  76   int  va,vb;
  77   real dx;
  78
  79   va = *(int *)a;
  80   vb = *(int *)b;
  81
  82   if ((dx = (xptr[va][XX] - xptr[vb][XX])) < 0)
  83     return -1;
  84   else if (dx > 0)
  85     return 1;
  86   else
  87     return 0;
  88 }
  89
  90 static void do_my_pme(FILE *fp,real tm,gmx_bool bVerbose,t_inputrec *ir,
  91                       rvec x[],rvec xbuf[],rvec f[],
  92                       real charge[],real qbuf[],real qqbuf[],
  93                       matrix box,gmx_bool bSort,
  94                       t_commrec *cr,t_nsborder *nsb,t_nrnb *nrnb,
  95                       t_block *excl,real qtot,
  96                       t_forcerec *fr,int index[],FILE *fp_xvg,
  97                       int ngroups,unsigned short cENER[])
  98 {
  99   real   ener,vcorr,q,xx,dvdl=0,vdip,vcharge;
 100   tensor vir,vir_corr,vir_tot;
 101   rvec   mu_tot[2];
 102   int    i,m,ii,ig,jg;
 103   real   **epme,*qptr;
 104
 105   /* Initiate local variables */
 106   fr->f_el_recip = f;
 107   clear_mat(vir);
 108   clear_mat(vir_corr);
 109
 110   if (ngroups > 1) {
 111     fprintf(fp,"There are %d energy groups\n",ngroups);
 112     snew(epme,ngroups);
 113     for(i=0; (i<ngroups); i++)
 114       snew(epme[i],ngroups);
 115   }
 116
 117   /* Put x is in the box, this part needs to be parallellized properly */
 118   /*put_atoms_in_box(box,nsb->natoms,x);*/
 119   /* Here sorting of X (and q) is done.
 120    * Alternatively, one could just put the atoms in one of the
 121    * cr->nnodes slabs. That is much cheaper than sorting.
 122    */
 123   for(i=0; (i<nsb->natoms); i++)
 124     index[i] = i;
 125   if (bSort) {
 126     xptr = x;
 127     qsort(index,nsb->natoms,sizeof(index[0]),comp_xptr);
 128     xptr = NULL; /* To trap unintentional use of the ptr */
 129   }
 130   /* After sorting we only need the part that is to be computed on
 131    * this processor. We also compute the mu_tot here (system dipole)
 132    */
 133   clear_rvec(mu_tot[0]);
 134   for(i=START(nsb); (i<START(nsb)+HOMENR(nsb)); i++) {
 135     ii      = index[i];
 136     q       = charge[ii];
 137     qbuf[i] = q;
 138     for(m=0; (m<DIM); m++) {
 139       xx         = x[ii][m];
 140       xbuf[i][m] = xx;
 141       mu_tot[0][m] += q*xx;
 142     }
 143     clear_rvec(f[ii]);
 144   }
 145   copy_rvec(mu_tot[0],mu_tot[1]);
 146   if (debug) {
 147     pr_rvec(debug,0,"qbuf",qbuf,nsb->natoms,TRUE);
 148     pr_rvecs(debug,0,"xbuf",xbuf,nsb->natoms);
 149     pr_rvecs(debug,0,"box",box,DIM);
 150   }
 151   for(ig=0; (ig<ngroups); ig++) {
 152     for(jg=ig; (jg<ngroups); jg++) {
 153       if (ngroups > 1) {
 154         for(i=START(nsb); (i<START(nsb)+HOMENR(nsb)); i++) {
 155           if ((cENER[i] == ig) || (cENER[i] == jg))
 156             qqbuf[i] = qbuf[i];
 157           else
 158             qqbuf[i] = 0;
 159         }
 160         qptr = qqbuf;
 161       }
 162       else
 163         qptr = qbuf;
 164       ener  = do_pme(fp,bVerbose,ir,xbuf,f,qptr,qptr,box,cr,
 165                      nsb,nrnb,vir,fr->ewaldcoeff,FALSE,0,&dvdl,FALSE);
 166       vcorr = ewald_LRcorrection(fp,nsb,cr,fr,qptr,qptr,excl,xbuf,box,mu_tot,
 167                                  ir->ewald_geometry,ir->epsilon_surface,
 168                                  0,&dvdl,&vdip,&vcharge);
 169       gmx_sum(1,&ener,cr);
 170       gmx_sum(1,&vcorr,cr);
 171       if (ngroups > 1)
 172         epme[ig][jg] = ener+vcorr;
 173     }
 174   }
 175   if (ngroups > 1) {
 176     if (fp_xvg)
 177       fprintf(fp_xvg,"%10.3f",tm);
 178     for(ig=0; (ig<ngroups); ig++) {
 179       for(jg=ig; (jg<ngroups); jg++) {
 180         if (ig != jg)
 181           epme[ig][jg] -= epme[ig][ig]+epme[jg][jg];
 182         if (fp_xvg)
 183           fprintf(fp_xvg,"  %12.5e",epme[ig][jg]);
 184       }
 185     }
 186     if (fp_xvg)
 187       fprintf(fp_xvg,"\n");
 188   }
 189   else {
 190     fprintf(fp,"Time: %10.3f Energy: %12.5e  Correction: %12.5e  Total: %12.5e\n",
 191             tm,ener,vcorr,ener+vcorr);
 192     if (fp_xvg)
 193       fprintf(fp_xvg,"%10.3f %12.5e %12.5e %12.5e\n",tm,ener+vcorr,vdip,vcharge);
 194     if (bVerbose) {
 195       m_add(vir,vir_corr,vir_tot);
 196       gmx_sum(9,vir_tot[0],cr);
 197       pr_rvecs(fp,0,"virial",vir_tot,DIM);
 198     }
 199     fflush(fp);
 200   }
 201 }
 202
 203 int main(int argc,char *argv[])
 204 {
 205   static char *desc[] = {
 206     "The [TT]pmetest[tt] program tests the scaling of the PME code. When only given",
 207     "a [TT].tpr[tt] file it will compute PME for one frame. When given a trajectory",
 208     "it will do so for all the frames in the trajectory. Before the PME",
 209     "routine is called the coordinates are sorted along the X-axis.[PAR]",
 210     "As an extra service to the public the program can also compute",
 211     "long-range Coulomb energies for components of the system. When the",
 212     "[TT]-groups[tt] flag is given to the program the energy groups",
 213     "from the [TT].tpr[tt] file will be read, and half an energy matrix computed."
 214   };
 215   t_commrec    *cr,*mcr;
 216   static t_filenm fnm[] = {
 217     { efTPX, NULL,      NULL,       ffREAD  },
 218     { efTRN, "-o",      NULL,       ffWRITE },
 219     { efLOG, "-g",      "pme",      ffWRITE },
 220     { efTRX, "-f",      NULL,       ffOPTRD },
 221     { efXVG, "-x",      "ener-pme", ffWRITE }
 222   };
 223 #define NFILE asize(fnm)
 224
 225   /* Command line options ! */
 226   static gmx_bool bVerbose=FALSE;
 227   static gmx_bool bOptFFT=FALSE;
 228   static gmx_bool bSort=FALSE;
 229   static int  ewald_geometry=eewg3D;
 230   static int  nnodes=1;
 231   static int  pme_order=0;
 232   static rvec grid = { -1, -1, -1 };
 233   static real rc   = 0.0;
 234   static real dtol = 0.0;
 235   static gmx_bool bGroups = FALSE;
 236   static t_pargs pa[] = {
 237     { "-np",      FALSE, etINT, {&nnodes},
 238       "Number of nodes, must be the same as used for [TT]grompp[tt]" },
 239     { "-v",       FALSE, etBOOL,{&bVerbose},
 240       "Be loud and noisy" },
 241     { "-sort",    FALSE, etBOOL,{&bSort},
 242       "Sort coordinates. Crucial for domain decomposition." },
 243     { "-grid",    FALSE, etRVEC,{&grid},
 244       "Number of grid cells in X, Y, Z dimension (if -1 use from [TT].tpr[tt])" },
 245     { "-order",   FALSE, etINT, {&pme_order},
 246       "Order of the PME spreading algorithm" },
 247     { "-groups",  FALSE, etBOOL, {&bGroups},
 248       "Compute half an energy matrix based on the energy groups in your [TT].tpr[tt] file" },
 249     { "-rc",      FALSE, etREAL, {&rc},
 250       "Rcoulomb for Ewald summation" },
 251     { "-tol",     FALSE, etREAL, {&dtol},
 252       "Tolerance for Ewald summation" }
 253   };
 254   FILE        *fp;
 255   t_inputrec  *ir;
 256   t_topology  top;
 257   t_tpxheader tpx;
 258   t_nrnb      nrnb;
 259   t_nsborder  *nsb;
 260   t_forcerec  *fr;
 261   t_mdatoms   *mdatoms;
 262   char        title[STRLEN];
 263   int         natoms,step,status,i,ncg,root;
 264   real        t,lambda,ewaldcoeff,qtot;
 265   rvec        *x,*f,*xbuf;
 266   int         *index;
 267   gmx_bool        bCont;
 268   real        *charge,*qbuf,*qqbuf;
 269   matrix      box;
 270
 271   /* Start the actual parallel code if necessary */
 272   cr   = init_par(&argc,&argv);
 273   root = 0;
 274
 275   if (MASTER(cr))
 276     CopyRight(stderr,argv[0]);
 277
 278   /* Parse command line on all processors, arguments are passed on in
 279    * init_par (see above)
 280    */
 281   parse_common_args(&argc,argv,
 282                     PCA_KEEP_ARGS | PCA_NOEXIT_ON_ARGS | PCA_BE_NICE |
 283                     PCA_CAN_SET_DEFFNM | (MASTER(cr) ? 0 : PCA_QUIET),
 284                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL);
 285
 286 #ifndef GMX_MPI
 287   if (nnodes > 1)
 288     gmx_fatal(FARGS,"GROMACS compiled without MPI support - can't do parallel runs");
 289 #endif
 290
 291   /* Open log files on all processors */
 292   open_log(ftp2fn(efLOG,NFILE,fnm),cr);
 293   snew(ir,1);
 294
 295   if (MASTER(cr)) {
 296     /* Read tpr file etc. */
 297     read_tpxheader(ftp2fn(efTPX,NFILE,fnm),&tpx,FALSE,NULL,NULL);
 298     snew(x,tpx.natoms);
 299     read_tpx(ftp2fn(efTPX,NFILE,fnm),&step,&t,&lambda,ir,
 300              box,&natoms,x,NULL,NULL,&top);
 301     /* Charges */
 302     qtot = 0;
 303     snew(charge,natoms);
 304     for(i=0; (i<natoms); i++) {
 305       charge[i] = top.atoms.atom[i].q;
 306       qtot += charge[i];
 307     }
 308
 309     /* Grid stuff */
 310     if (opt2parg_bSet("-grid",asize(pa),pa)) {
 311       ir->nkx = grid[XX];
 312       ir->nky = grid[YY];
 313       ir->nkz = grid[ZZ];
 314     }
 315     /* Check command line parameters for consistency */
 316     if ((ir->nkx <= 0) || (ir->nky <= 0) || (ir->nkz <= 0))
 317       gmx_fatal(FARGS,"PME grid = %d %d %d",ir->nkx,ir->nky,ir->nkz);
 318     if (opt2parg_bSet("-rc",asize(pa),pa))
 319       ir->rcoulomb = rc;
 320     if (ir->rcoulomb <= 0)
 321       gmx_fatal(FARGS,"rcoulomb should be > 0 (not %f)",ir->rcoulomb);
 322     if (opt2parg_bSet("-order",asize(pa),pa))
 323       ir->pme_order = pme_order;
 324     if (ir->pme_order <= 0)
 325       gmx_fatal(FARGS,"pme_order should be > 0 (not %d)",ir->pme_order);
 326     if (opt2parg_bSet("-tol",asize(pa),pa))
 327       ir->ewald_rtol = dtol;
 328     if (ir->ewald_rtol <= 0)
 329       gmx_fatal(FARGS,"ewald_tol should be > 0 (not %f)",ir->ewald_rtol);
 330   }
 331   else {
 332     init_top(&top);
 333   }
 334
 335   /* Add parallellization code here */
 336   snew(nsb,1);
 337   if (MASTER(cr)) {
 338     ncg = top.blocks[ebCGS].multinr[0];
 339     for(i=0; (i<cr->nnodes-1); i++)
 340       top.blocks[ebCGS].multinr[i] = min(ncg,(ncg*(i+1))/cr->nnodes);
 341     for( ; (i<MAXNODES); i++)
 342       top.blocks[ebCGS].multinr[i] = ncg;
 343   }
 344   if (PAR(cr)) {
 345     /* Set some variables to zero to avoid core dumps */
 346     ir->opts.ngtc = ir->opts.ngacc = ir->opts.ngfrz = ir->opts.ngener = 0;
 347 #ifdef GMX_MPI
 348     /* Distribute the data over processors */
 349     MPI_Bcast(&natoms,1,MPI_INT,root,MPI_COMM_WORLD);
 350     MPI_Bcast(ir,sizeof(*ir),MPI_BYTE,root,MPI_COMM_WORLD);
 351     MPI_Bcast(&qtot,1,GMX_MPI_REAL,root,MPI_COMM_WORLD);
 352 #endif
 353
 354     /* Call some dedicated communication routines, master sends n-1 times */
 355     if (MASTER(cr)) {
 356       for(i=1; (i<cr->nnodes); i++) {
 357         mv_block(i,&(top.blocks[ebCGS]));
 358         mv_block(i,&(top.atoms.excl));
 359       }
 360     }
 361     else {
 362       ld_block(root,&(top.blocks[ebCGS]));
 363       ld_block(root,&(top.atoms.excl));
 364     }
 365     if (!MASTER(cr)) {
 366       snew(charge,natoms);
 367       snew(x,natoms);
 368     }
 369 #ifdef GMX_MPI
 370     MPI_Bcast(charge,natoms,GMX_MPI_REAL,root,MPI_COMM_WORLD);
 371 #endif
 372   }
 373   ewaldcoeff = calc_ewaldcoeff(ir->rcoulomb,ir->ewald_rtol);
 374
 375
 376   if (bVerbose)
 377     pr_inputrec(stdlog,0,"Inputrec",ir);
 378
 379   /* Allocate memory for temp arrays etc. */
 380   snew(xbuf,natoms);
 381   snew(f,natoms);
 382   snew(qbuf,natoms);
 383   snew(qqbuf,natoms);
 384   snew(index,natoms);
 385
 386   /* Initialize the PME code */
 387   init_pme(stdlog,cr,ir->nkx,ir->nky,ir->nkz,ir->pme_order,
 388            natoms,FALSE,bOptFFT,ewald_geometry);
 389
 390   /* MFlops accounting */
 391   init_nrnb(&nrnb);
 392
 393   /* Initialize the work division */
 394   calc_nsb(stdlog,&(top.blocks[ebCGS]),cr->nnodes,nsb,0);
 395   nsb->nodeid = cr->nodeid;
 396   print_nsb(stdlog,"pmetest",nsb);
 397
 398   /* Initiate forcerec */
 399   mdatoms = atoms2md(stdlog,&top.atoms,ir->opts.nFreeze,ir->eI,
 400                      ir->delta_t,0,ir->opts.tau_t,FALSE,FALSE);
 401   snew(fr,1);
 402   init_forcerec(stdlog,fr,ir,&top,cr,mdatoms,nsb,box,FALSE,NULL,NULL,FALSE);
 403
 404   /* First do PME based on coordinates in tpr file, send them to
 405    * other processors if needed.
 406    */
 407   if (MASTER(cr))
 408     fprintf(stdlog,"-----\n"
 409             "Results based on tpr file %s\n",ftp2fn(efTPX,NFILE,fnm));
 410 #ifdef GMX_MPI
 411   if (PAR(cr)) {
 412     MPI_Bcast(x[0],natoms*DIM,GMX_MPI_REAL,root,MPI_COMM_WORLD);
 413     MPI_Bcast(box[0],DIM*DIM,GMX_MPI_REAL,root,MPI_COMM_WORLD);
 414     MPI_Bcast(&t,1,GMX_MPI_REAL,root,MPI_COMM_WORLD);
 415   }
 416 #endif
 417   do_my_pme(stdlog,0,bVerbose,ir,x,xbuf,f,charge,qbuf,qqbuf,box,bSort,
 418             cr,nsb,&nrnb,&(top.atoms.excl),qtot,fr,index,NULL,
 419             bGroups ? ir->opts.ngener : 1,mdatoms->cENER);
 420
 421   /* If we have a trajectry file, we will read the frames in it and compute
 422    * the PME energy.
 423    */
 424   if (ftp2bSet(efTRX,NFILE,fnm)) {
 425     fprintf(stdlog,"-----\n"
 426             "Results based on trx file %s\n",ftp2fn(efTRX,NFILE,fnm));
 427     if (MASTER(cr)) {
 428       sfree(x);
 429       natoms = read_first_x(&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
 430       if (natoms != top.atoms.nr)
 431         gmx_fatal(FARGS,"natoms in trx = %d, in tpr = %d",natoms,top.atoms.nr);
 432       fp = xvgropen(ftp2fn(efXVG,NFILE,fnm),"PME Energy","Time (ps)","E (kJ/mol)");
 433     }
 434     else
 435       fp = NULL;
 436     do {
 437       /* Send coordinates, box and time to the other nodes */
 438 #ifdef GMX_MPI
 439       if (PAR(cr)) {
 440         MPI_Bcast(x[0],natoms*DIM,GMX_MPI_REAL,root,MPI_COMM_WORLD);
 441         MPI_Bcast(box[0],DIM*DIM,GMX_MPI_REAL,root,MPI_COMM_WORLD);
 442         MPI_Bcast(&t,1,GMX_MPI_REAL,root,MPI_COMM_WORLD);
 443       }
 444 #endif
 445       rm_pbc(&top.idef,nsb->natoms,box,x,x);
 446       /* Call the PME wrapper function */
 447       do_my_pme(stdlog,t,bVerbose,ir,x,xbuf,f,charge,qbuf,qqbuf,box,bSort,cr,
 448                 nsb,&nrnb,&(top.atoms.excl),qtot,fr,index,fp,
 449                 bGroups ? ir->opts.ngener : 1,mdatoms->cENER);
 450       /* Only the master processor reads more data */
 451       if (MASTER(cr))
 452           bCont = read_next_x(status,&t,natoms,x,box);
 453       /* Check whether we need to continue */
 454 #ifdef GMX_MPI
 455       if (PAR(cr))
 456           MPI_Bcast(&bCont,1,MPI_INT,root,MPI_COMM_WORLD);
 457 #endif
 458
 459     } while (bCont);
 460
 461     /* Finish I/O, close files */
 462     if (MASTER(cr)) {
 463       close_trx(status);
 464       ffclose(fp);
 465     }
 466   }
 467
 468   if (bVerbose) {
 469     /* Do some final I/O about performance, might be useful in debugging */
 470     fprintf(stdlog,"-----\n");
 471     print_nrnb(stdlog,&nrnb);
 472   }
 473
 474   /* Finish the parallel stuff */
 475   if (gmx_parallel_env_initialized())
 476     gmx_finalize(cr);
 477
 478   /* Thank the audience, as usual */
 479   if (MASTER(cr))
 480     thanx(stderr);
 481
 482   return 0;
 483 }
 484