src/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 "perf_est.h"
  39 #include "physics.h"
  40 #include "vec.h"
  41 #include "mtop_util.h"
  42
  43 int n_bonded_dx(gmx_mtop_t *mtop,bool bExcl)
  44 {
  45   int mb,nmol,ftype,ndxb,ndx_excl;
  46   int ndx;
  47   gmx_moltype_t *molt;
  48
  49   /* Count the number of pbc_rvec_sub calls required for bonded interactions.
  50    * This number is also roughly proportional to the computational cost.
  51    */
  52   ndx = 0;
  53   ndx_excl = 0;
  54   for(mb=0; mb<mtop->nmolblock; mb++) {
  55     molt = &mtop->moltype[mtop->molblock[mb].type];
  56     nmol = mtop->molblock[mb].nmol;
  57     for(ftype=0; ftype<F_NRE; ftype++) {
  58       if (interaction_function[ftype].flags & IF_BOND) {
  59         switch (ftype) {
  60                 case F_POSRES:    ndxb = 1; break;
  61                 case F_CONNBONDS: ndxb = 0; break;
  62                 default:     ndxb = NRAL(ftype) - 1; break;
  63         }
  64         ndx += nmol*ndxb*molt->ilist[ftype].nr/(1 + NRAL(ftype));
  65       }
  66     }
  67     if (bExcl) {
  68       ndx_excl += nmol*(molt->excls.nra - molt->atoms.nr)/2;
  69     } else {
  70       ndx_excl = 0;
  71     }
  72   }
  73
  74   if (debug)
  75     fprintf(debug,"ndx bonded %d exclusions %d\n",ndx,ndx_excl);
  76
  77   ndx += ndx_excl;
  78
  79   return ndx;
  80 }
  81
  82 float pme_load_estimate(gmx_mtop_t *mtop,t_inputrec *ir,matrix box)
  83 {
  84   t_atom *atom;
  85   int  mb,nmol,atnr,cg,a,a0,ncqlj,ncq,nclj;
  86   bool bBHAM,bLJcut,bChargePerturbed,bWater,bQ,bLJ;
  87   double nw,nqlj,nq,nlj,cost_bond,cost_pp,cost_spread,cost_fft;
  88   float fq,fqlj,flj,fljtab,fqljw,fqw,fqspread,ffft,fbond;
  89   float ratio;
  90   t_iparams *iparams;
  91   gmx_moltype_t *molt;
  92
  93   bBHAM = (mtop->ffparams.functype[0] == F_BHAM);
  94
  95   bLJcut = ((ir->vdwtype == evdwCUT) && !bBHAM);
  96
  97   /* Computational cost relative to a tabulated q-q interaction.
  98    * This will be machine dependent.
  99    * The numbers here are accurate for Intel Core2 and AMD Athlon 64
 100    * in single precision. In double precision PME mesh is slightly cheaper,
 101    * although not so much that the numbers need to be adjusted.
 102    */
 103   fq    = 1.0;
 104   fqlj  = (bLJcut ? 1.5  : 2.0 );
 105   flj   = (bLJcut ? 0.5  : 1.5 );
 106   /* Cost of 1 water with one Q/LJ atom */
 107   fqljw = (bLJcut ? 1.75 : 2.25);
 108   /* Cost of 1 water with one Q atom or with 1/3 water (LJ negligible) */
 109   fqw   = 1.5;
 110   /* Cost of q spreading and force interpolation per charge */
 111   fqspread = 25.0;
 112   /* Cost of fft's + pme_solve, will be multiplied with N log(N) */
 113   ffft     =  0.4;
 114   /* Cost of a bonded interaction divided by the number of (pbc_)dx required */
 115   fbond = 5.0;
 116
 117   iparams = mtop->ffparams.iparams;
 118   atnr = mtop->ffparams.atnr;
 119   nw   = 0;
 120   nqlj = 0;
 121   nq   = 0;
 122   nlj  = 0;
 123   bChargePerturbed = FALSE;
 124   for(mb=0; mb<mtop->nmolblock; mb++) {
 125     molt = &mtop->moltype[mtop->molblock[mb].type];
 126     atom = molt->atoms.atom;
 127     nmol = mtop->molblock[mb].nmol;
 128     a = 0;
 129     for(cg=0; cg<molt->cgs.nr; cg++) {
 130       bWater = !bBHAM;
 131       ncqlj = 0;
 132       ncq   = 0;
 133       nclj  = 0;
 134       a0    = a;
 135       while (a < molt->cgs.index[cg+1]) {
 136         bQ  = (atom[a].q != 0 || atom[a].qB != 0);
 137         bLJ = (iparams[(atnr+1)*atom[a].type].lj.c6  != 0 ||
 138                iparams[(atnr+1)*atom[a].type].lj.c12 != 0);
 139         if (atom[a].q != atom[a].qB) {
 140           bChargePerturbed = TRUE;
 141         }
 142         /* This if this atom fits into water optimization */
 143         if (!((a == a0   &&  bQ &&  bLJ) ||
 144               (a == a0+1 &&  bQ && !bLJ) ||
 145               (a == a0+2 &&  bQ && !bLJ && atom[a].q == atom[a-1].q) ||
 146               (a == a0+3 && !bQ &&  bLJ)))
 147           bWater = FALSE;
 148         if (bQ && bLJ) {
 149           ncqlj++;
 150         } else {
 151           if (bQ)
 152             ncq++;
 153           if (bLJ)
 154             nclj++;
 155         }
 156         a++;
 157       }
 158       if (bWater) {
 159         nw   += nmol;
 160       } else {
 161         nqlj += nmol*ncqlj;
 162         nq   += nmol*ncq;
 163         nlj  += nmol*nclj;
 164       }
 165     }
 166   }
 167   if (debug)
 168     fprintf(debug,"nw %g nqlj %g nq %g nlj %g\n",nw,nqlj,nq,nlj);
 169
 170   cost_bond = fbond*n_bonded_dx(mtop,TRUE);
 171
 172   /* For the PP non-bonded cost it is (unrealistically) assumed
 173    * that all atoms are distributed homogeneously in space.
 174    */
 175   cost_pp = 0.5*(fqljw*nw*nqlj +
 176                  fqw  *nw*(3*nw + nq) +
 177                  fqlj *nqlj*nqlj +
 178                  fq   *nq*(3*nw + nqlj + nq) +
 179                  flj  *nlj*(nw + nqlj + nlj))
 180     *4/3*M_PI*ir->rlist*ir->rlist*ir->rlist/det(box);
 181
 182   cost_spread = fqspread*(3*nw + nqlj + nq);
 183   cost_fft    = ffft*ir->nkx*ir->nky*ir->nkz*log(ir->nkx*ir->nky*ir->nkz);
 184
 185   if (ir->efep != efepNO && bChargePerturbed) {
 186     /* All PME work, except the spline coefficient calculation, doubles */
 187     cost_spread *= 2;
 188     cost_fft    *= 2;
 189   }
 190
 191   ratio =
 192     (cost_spread + cost_fft)/(cost_bond + cost_pp + cost_spread + cost_fft);
 193
 194   if (debug) {
 195     fprintf(debug,
 196             "cost_bond   %f\n"
 197             "cost_pp     %f\n"
 198             "cost_spread %f\n"
 199             "cost_fft    %f\n",
 200             cost_bond,cost_pp,cost_spread,cost_fft);
 201
 202     fprintf(debug,"Estimate for relative PME load: %.3f\n",ratio);
 203   }
 204
 205   return ratio;
 206 }