src/gromacs/mdlib/perf_est.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  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  10  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  11  * Copyright (c) 2001-2008, The GROMACS development team,
  12  * check out http://www.gromacs.org for more information.
  13
  14  * This program is free software; you can redistribute it and/or
  15  * modify it under the terms of the GNU General Public License
  16  * as published by the Free Software Foundation; either version 2
  17  * of the License, or (at your option) any later version.
  18  *
  19  * If you want to redistribute modifications, please consider that
  20  * scientific software is very special. Version control is crucial -
  21  * bugs must be traceable. We will be happy to consider code for
  22  * inclusion in the official distribution, but derived work must not
  23  * be called official GROMACS. Details are found in the README & COPYING
  24  * files - if they are missing, get the official version at www.gromacs.org.
  25  *
  26  * To help us fund GROMACS development, we humbly ask that you cite
  27  * the papers on the package - you can find them in the top README file.
  28  *
  29  * For more info, check our website at http://www.gromacs.org
  30  *
  31  * And Hey:
  32  * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
  33  */
  34 #ifdef HAVE_CONFIG_H
  35 #include <config.h>
  36 #endif
  37
  38 #include <math.h>
  39
  40 #include "perf_est.h"
  41 #include "physics.h"
  42 #include "vec.h"
  43 #include "mtop_util.h"
  44
  45 int n_bonded_dx(gmx_mtop_t *mtop,gmx_bool bExcl)
  46 {
  47   int mb,nmol,ftype,ndxb,ndx_excl;
  48   int ndx;
  49   gmx_moltype_t *molt;
  50
  51   /* Count the number of pbc_rvec_sub calls required for bonded interactions.
  52    * This number is also roughly proportional to the computational cost.
  53    */
  54   ndx = 0;
  55   ndx_excl = 0;
  56   for(mb=0; mb<mtop->nmolblock; mb++) {
  57     molt = &mtop->moltype[mtop->molblock[mb].type];
  58     nmol = mtop->molblock[mb].nmol;
  59     for(ftype=0; ftype<F_NRE; ftype++) {
  60       if (interaction_function[ftype].flags & IF_BOND) {
  61         switch (ftype) {
  62                 case F_POSRES:    ndxb = 1; break;
  63                 case F_CONNBONDS: ndxb = 0; break;
  64                 default:     ndxb = NRAL(ftype) - 1; break;
  65         }
  66         ndx += nmol*ndxb*molt->ilist[ftype].nr/(1 + NRAL(ftype));
  67       }
  68     }
  69     if (bExcl) {
  70       ndx_excl += nmol*(molt->excls.nra - molt->atoms.nr)/2;
  71     } else {
  72       ndx_excl = 0;
  73     }
  74   }
  75
  76   if (debug)
  77     fprintf(debug,"ndx bonded %d exclusions %d\n",ndx,ndx_excl);
  78
  79   ndx += ndx_excl;
  80
  81   return ndx;
  82 }
  83
  84 float pme_load_estimate(gmx_mtop_t *mtop,t_inputrec *ir,matrix box)
  85 {
  86   t_atom *atom;
  87   int  mb,nmol,atnr,cg,a,a0,ncqlj,ncq,nclj;
  88   gmx_bool bBHAM,bLJcut,bChargePerturbed,bWater,bQ,bLJ;
  89   double nw,nqlj,nq,nlj;
  90   double cost_bond,cost_pp,cost_spread,cost_fft,cost_solve,cost_pme;
  91   float fq,fqlj,flj,fljtab,fqljw,fqw,fqspread,ffft,fsolve,fbond;
  92   float ratio;
  93   t_iparams *iparams;
  94   gmx_moltype_t *molt;
  95
  96   bBHAM = (mtop->ffparams.functype[0] == F_BHAM);
  97
  98   bLJcut = ((ir->vdwtype == evdwCUT) && !bBHAM);
  99
 100   /* Computational cost of bonded, non-bonded and PME calculations.
 101    * This will be machine dependent.
 102    * The numbers here are accurate for Intel Core2 and AMD Athlon 64
 103    * in single precision. In double precision PME mesh is slightly cheaper,
 104    * although not so much that the numbers need to be adjusted.
 105    */
 106   fq    = 1.5;
 107   fqlj  = (bLJcut ? 1.5  : 2.0 );
 108   flj   = (bLJcut ? 1.0  : 1.75);
 109   /* Cost of 1 water with one Q/LJ atom */
 110   fqljw = (bLJcut ? 2.0  : 2.25);
 111   /* Cost of 1 water with one Q atom or with 1/3 water (LJ negligible) */
 112   fqw   = 1.75;
 113   /* Cost of q spreading and force interpolation per charge (mainly memory) */
 114   fqspread = 0.55;
 115   /* Cost of fft's, will be multiplied with N log(N) */
 116   ffft     = 0.20;
 117   /* Cost of pme_solve, will be multiplied with N */
 118   fsolve   = 0.80;
 119   /* Cost of a bonded interaction divided by the number of (pbc_)dx nrequired */
 120   fbond = 5.0;
 121
 122   iparams = mtop->ffparams.iparams;
 123   atnr = mtop->ffparams.atnr;
 124   nw   = 0;
 125   nqlj = 0;
 126   nq   = 0;
 127   nlj  = 0;
 128   bChargePerturbed = FALSE;
 129   for(mb=0; mb<mtop->nmolblock; mb++) {
 130     molt = &mtop->moltype[mtop->molblock[mb].type];
 131     atom = molt->atoms.atom;
 132     nmol = mtop->molblock[mb].nmol;
 133     a = 0;
 134     for(cg=0; cg<molt->cgs.nr; cg++) {
 135       bWater = !bBHAM;
 136       ncqlj = 0;
 137       ncq   = 0;
 138       nclj  = 0;
 139       a0    = a;
 140       while (a < molt->cgs.index[cg+1]) {
 141         bQ  = (atom[a].q != 0 || atom[a].qB != 0);
 142         bLJ = (iparams[(atnr+1)*atom[a].type].lj.c6  != 0 ||
 143                iparams[(atnr+1)*atom[a].type].lj.c12 != 0);
 144         if (atom[a].q != atom[a].qB) {
 145           bChargePerturbed = TRUE;
 146         }
 147         /* This if this atom fits into water optimization */
 148         if (!((a == a0   &&  bQ &&  bLJ) ||
 149               (a == a0+1 &&  bQ && !bLJ) ||
 150               (a == a0+2 &&  bQ && !bLJ && atom[a].q == atom[a-1].q) ||
 151               (a == a0+3 && !bQ &&  bLJ)))
 152           bWater = FALSE;
 153         if (bQ && bLJ) {
 154           ncqlj++;
 155         } else {
 156           if (bQ)
 157             ncq++;
 158           if (bLJ)
 159             nclj++;
 160         }
 161         a++;
 162       }
 163       if (bWater) {
 164         nw   += nmol;
 165       } else {
 166         nqlj += nmol*ncqlj;
 167         nq   += nmol*ncq;
 168         nlj  += nmol*nclj;
 169       }
 170     }
 171   }
 172   if (debug)
 173     fprintf(debug,"nw %g nqlj %g nq %g nlj %g\n",nw,nqlj,nq,nlj);
 174
 175   cost_bond = fbond*n_bonded_dx(mtop,TRUE);
 176
 177   /* For the PP non-bonded cost it is (unrealistically) assumed
 178    * that all atoms are distributed homogeneously in space.
 179    */
 180   cost_pp = 0.5*(fqljw*nw*nqlj +
 181                  fqw  *nw*(3*nw + nq) +
 182                  fqlj *nqlj*nqlj +
 183                  fq   *nq*(3*nw + nqlj + nq) +
 184                  flj  *nlj*(nw + nqlj + nlj))
 185     *4/3*M_PI*ir->rlist*ir->rlist*ir->rlist/det(box);
 186
 187   cost_spread = fqspread*(3*nw + nqlj + nq)*pow(ir->pme_order,3);
 188   cost_fft    = ffft*ir->nkx*ir->nky*ir->nkz*log(ir->nkx*ir->nky*ir->nkz);
 189   cost_solve  = fsolve*ir->nkx*ir->nky*ir->nkz;
 190
 191   if (ir->efep != efepNO && bChargePerturbed) {
 192     /* All PME work, except the spline coefficient calculation, doubles */
 193     cost_spread *= 2;
 194     cost_fft    *= 2;
 195     cost_solve  *= 2;
 196   }
 197
 198   cost_pme = cost_spread + cost_fft + cost_solve;
 199
 200   ratio = cost_pme/(cost_bond + cost_pp + cost_pme);
 201
 202   if (debug) {
 203     fprintf(debug,
 204             "cost_bond   %f\n"
 205             "cost_pp     %f\n"
 206             "cost_spread %f\n"
 207             "cost_fft    %f\n"
 208             "cost_solve  %f\n",
 209             cost_bond,cost_pp,cost_spread,cost_fft,cost_solve);
 210
 211     fprintf(debug,"Estimate for relative PME load: %.3f\n",ratio);
 212   }
 213
 214   return ratio;
 215 }