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;
  88   double cost_bond,cost_pp,cost_spread,cost_fft,cost_solve,cost_pme;
  89   float fq,fqlj,flj,fljtab,fqljw,fqw,fqspread,ffft,fsolve,fbond;
  90   float ratio;
  91   t_iparams *iparams;
  92   gmx_moltype_t *molt;
  93
  94   bBHAM = (mtop->ffparams.functype[0] == F_BHAM);
  95
  96   bLJcut = ((ir->vdwtype == evdwCUT) && !bBHAM);
  97
  98   /* Computational cost of bonded, non-bonded and PME calculations.
  99    * This will be machine dependent.
 100    * The numbers here are accurate for Intel Core2 and AMD Athlon 64
 101    * in single precision. In double precision PME mesh is slightly cheaper,
 102    * although not so much that the numbers need to be adjusted.
 103    */
 104   fq    = 1.5;
 105   fqlj  = (bLJcut ? 1.5  : 2.0 );
 106   flj   = (bLJcut ? 1.0  : 1.75);
 107   /* Cost of 1 water with one Q/LJ atom */
 108   fqljw = (bLJcut ? 2.0  : 2.25);
 109   /* Cost of 1 water with one Q atom or with 1/3 water (LJ negligible) */
 110   fqw   = 1.75;
 111   /* Cost of q spreading and force interpolation per charge */
 112   fqspread = 35.0;
 113   /* Cost of fft's, will be multiplied with N log(N) */
 114   ffft     =  0.25;
 115   /* Cost of pme_solve, will be multiplied with N */
 116   fsolve   =  1.75;
 117   /* Cost of a bonded interaction divided by the number of (pbc_)dx required */
 118   fbond = 5.0;
 119
 120   iparams = mtop->ffparams.iparams;
 121   atnr = mtop->ffparams.atnr;
 122   nw   = 0;
 123   nqlj = 0;
 124   nq   = 0;
 125   nlj  = 0;
 126   bChargePerturbed = FALSE;
 127   for(mb=0; mb<mtop->nmolblock; mb++) {
 128     molt = &mtop->moltype[mtop->molblock[mb].type];
 129     atom = molt->atoms.atom;
 130     nmol = mtop->molblock[mb].nmol;
 131     a = 0;
 132     for(cg=0; cg<molt->cgs.nr; cg++) {
 133       bWater = !bBHAM;
 134       ncqlj = 0;
 135       ncq   = 0;
 136       nclj  = 0;
 137       a0    = a;
 138       while (a < molt->cgs.index[cg+1]) {
 139         bQ  = (atom[a].q != 0 || atom[a].qB != 0);
 140         bLJ = (iparams[(atnr+1)*atom[a].type].lj.c6  != 0 ||
 141                iparams[(atnr+1)*atom[a].type].lj.c12 != 0);
 142         if (atom[a].q != atom[a].qB) {
 143           bChargePerturbed = TRUE;
 144         }
 145         /* This if this atom fits into water optimization */
 146         if (!((a == a0   &&  bQ &&  bLJ) ||
 147               (a == a0+1 &&  bQ && !bLJ) ||
 148               (a == a0+2 &&  bQ && !bLJ && atom[a].q == atom[a-1].q) ||
 149               (a == a0+3 && !bQ &&  bLJ)))
 150           bWater = FALSE;
 151         if (bQ && bLJ) {
 152           ncqlj++;
 153         } else {
 154           if (bQ)
 155             ncq++;
 156           if (bLJ)
 157             nclj++;
 158         }
 159         a++;
 160       }
 161       if (bWater) {
 162         nw   += nmol;
 163       } else {
 164         nqlj += nmol*ncqlj;
 165         nq   += nmol*ncq;
 166         nlj  += nmol*nclj;
 167       }
 168     }
 169   }
 170   if (debug)
 171     fprintf(debug,"nw %g nqlj %g nq %g nlj %g\n",nw,nqlj,nq,nlj);
 172
 173   cost_bond = fbond*n_bonded_dx(mtop,TRUE);
 174
 175   /* For the PP non-bonded cost it is (unrealistically) assumed
 176    * that all atoms are distributed homogeneously in space.
 177    */
 178   cost_pp = 0.5*(fqljw*nw*nqlj +
 179                  fqw  *nw*(3*nw + nq) +
 180                  fqlj *nqlj*nqlj +
 181                  fq   *nq*(3*nw + nqlj + nq) +
 182                  flj  *nlj*(nw + nqlj + nlj))
 183     *4/3*M_PI*ir->rlist*ir->rlist*ir->rlist/det(box);
 184
 185   cost_spread = fqspread*(3*nw + nqlj + nq);
 186   cost_fft    = ffft*ir->nkx*ir->nky*ir->nkz*log(ir->nkx*ir->nky*ir->nkz);
 187   cost_solve  = fsolve*ir->nkx*ir->nky*ir->nkz;
 188
 189   if (ir->efep != efepNO && bChargePerturbed) {
 190     /* All PME work, except the spline coefficient calculation, doubles */
 191     cost_spread *= 2;
 192     cost_fft    *= 2;
 193     cost_solve  *= 2;
 194   }
 195
 196   cost_pme = cost_spread + cost_fft + cost_solve;
 197
 198   ratio = cost_pme/(cost_bond + cost_pp + cost_pme);
 199
 200   if (debug) {
 201     fprintf(debug,
 202             "cost_bond   %f\n"
 203             "cost_pp     %f\n"
 204             "cost_spread %f\n"
 205             "cost_fft    %f\n"
 206             "cost_solve  %f\n",
 207             cost_bond,cost_pp,cost_spread,cost_fft,cost_solve);
 208
 209     fprintf(debug,"Estimate for relative PME load: %.3f\n",ratio);
 210   }
 211
 212   return ratio;
 213 }