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