3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
29 * For more info, check our website at http://www.gromacs.org
32 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
43 #include "mtop_util.h"
45 int n_bonded_dx(gmx_mtop_t *mtop,gmx_bool bExcl)
47 int mb,nmol,ftype,ndxb,ndx_excl;
51 /* Count the number of pbc_rvec_sub calls required for bonded interactions.
52 * This number is also roughly proportional to the computational cost.
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) {
63 case F_FBPOSRES: ndxb = 1; break;
64 case F_CONNBONDS: ndxb = 0; break;
65 default: ndxb = NRAL(ftype) - 1; break;
67 ndx += nmol*ndxb*molt->ilist[ftype].nr/(1 + NRAL(ftype));
71 ndx_excl += nmol*(molt->excls.nra - molt->atoms.nr)/2;
78 fprintf(debug,"ndx bonded %d exclusions %d\n",ndx,ndx_excl);
85 float pme_load_estimate(gmx_mtop_t *mtop,t_inputrec *ir,matrix box)
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;
97 bBHAM = (mtop->ffparams.functype[0] == F_BHAM);
99 bLJcut = ((ir->vdwtype == evdwCUT) && !bBHAM);
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.
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) */
114 /* Cost of q spreading and force interpolation per charge (mainly memory) */
116 /* Cost of fft's, will be multiplied with N log(N) */
118 /* Cost of pme_solve, will be multiplied with N */
120 /* Cost of a bonded interaction divided by the number of (pbc_)dx nrequired */
123 iparams = mtop->ffparams.iparams;
124 atnr = mtop->ffparams.atnr;
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;
135 for(cg=0; cg<molt->cgs.nr; cg++) {
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;
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)))
174 fprintf(debug,"nw %g nqlj %g nq %g nlj %g\n",nw,nqlj,nq,nlj);
176 cost_bond = fbond*n_bonded_dx(mtop,TRUE);
178 /* For the PP non-bonded cost it is (unrealistically) assumed
179 * that all atoms are distributed homogeneously in space.
181 cost_pp = 0.5*(fqljw*nw*nqlj +
182 fqw *nw*(3*nw + nq) +
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);
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;
192 if (ir->efep != efepNO && bChargePerturbed) {
193 /* All PME work, except the spline coefficient calculation, doubles */
199 cost_pme = cost_spread + cost_fft + cost_solve;
201 ratio = cost_pme/(cost_bond + cost_pp + cost_pme);
210 cost_bond,cost_pp,cost_spread,cost_fft,cost_solve);
212 fprintf(debug,"Estimate for relative PME load: %.3f\n",ratio);