2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team,
6 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012,2013, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
44 #include <sys/types.h>
49 #include "gmx_fatal.h"
53 #include "calc_verletbuf.h"
54 #include "../mdlib/nbnxn_consts.h"
57 /* The include below sets the SIMD instruction type (precision+width)
58 * for all nbnxn SIMD search and non-bonded kernel code.
60 #ifdef GMX_NBNXN_HALF_WIDTH_SIMD
61 #define GMX_USE_HALF_WIDTH_SIMD_HERE
63 #include "gmx_simd_macros.h"
66 /* Struct for unique atom type for calculating the energy drift.
67 * The atom displacement depends on mass and constraints.
68 * The energy jump for given distance depend on LJ type and q.
73 int type; /* type (used for LJ parameters) */
75 int con; /* constrained: 0, else 1, if 1, use #DOF=2 iso 3 */
76 int n; /* total #atoms of this type in the system */
77 } verletbuf_atomtype_t;
80 void verletbuf_get_list_setup(gmx_bool bGPU,
81 verletbuf_list_setup_t *list_setup)
83 list_setup->cluster_size_i = NBNXN_CPU_CLUSTER_I_SIZE;
87 list_setup->cluster_size_j = NBNXN_GPU_CLUSTER_SIZE;
91 #ifndef GMX_NBNXN_SIMD
92 list_setup->cluster_size_j = NBNXN_CPU_CLUSTER_I_SIZE;
94 list_setup->cluster_size_j = GMX_SIMD_WIDTH_HERE;
95 #ifdef GMX_NBNXN_SIMD_2XNN
96 /* We assume the smallest cluster size to be on the safe side */
97 list_setup->cluster_size_j /= 2;
103 static void add_at(verletbuf_atomtype_t **att_p, int *natt_p,
104 real mass, int type, real q, int con, int nmol)
106 verletbuf_atomtype_t *att;
111 /* Ignore massless particles */
120 !(mass == att[i].mass &&
121 type == att[i].type &&
135 srenew(*att_p, *natt_p);
136 (*att_p)[i].mass = mass;
137 (*att_p)[i].type = type;
139 (*att_p)[i].con = con;
140 (*att_p)[i].n = nmol;
144 static void get_verlet_buffer_atomtypes(const gmx_mtop_t *mtop,
145 verletbuf_atomtype_t **att_p,
149 verletbuf_atomtype_t *att;
151 int mb, nmol, ft, i, j, a1, a2, a3, a;
152 const t_atoms *atoms;
156 real *con_m, *vsite_m, cam[5];
161 if (n_nonlin_vsite != NULL)
166 for (mb = 0; mb < mtop->nmolblock; mb++)
168 nmol = mtop->molblock[mb].nmol;
170 atoms = &mtop->moltype[mtop->molblock[mb].type].atoms;
172 /* Check for constraints, as they affect the kinetic energy */
173 snew(con_m, atoms->nr);
174 snew(vsite_m, atoms->nr);
176 for (ft = F_CONSTR; ft <= F_CONSTRNC; ft++)
178 il = &mtop->moltype[mtop->molblock[mb].type].ilist[ft];
180 for (i = 0; i < il->nr; i += 1+NRAL(ft))
182 a1 = il->iatoms[i+1];
183 a2 = il->iatoms[i+2];
184 con_m[a1] += atoms->atom[a2].m;
185 con_m[a2] += atoms->atom[a1].m;
189 il = &mtop->moltype[mtop->molblock[mb].type].ilist[F_SETTLE];
191 for (i = 0; i < il->nr; i += 1+NRAL(F_SETTLE))
193 a1 = il->iatoms[i+1];
194 a2 = il->iatoms[i+2];
195 a3 = il->iatoms[i+3];
196 con_m[a1] += atoms->atom[a2].m + atoms->atom[a3].m;
197 con_m[a2] += atoms->atom[a1].m + atoms->atom[a3].m;
198 con_m[a3] += atoms->atom[a1].m + atoms->atom[a2].m;
201 /* Check for virtual sites, determine mass from constructing atoms */
202 for (ft = 0; ft < F_NRE; ft++)
206 il = &mtop->moltype[mtop->molblock[mb].type].ilist[ft];
208 for (i = 0; i < il->nr; i += 1+NRAL(ft))
210 ip = &mtop->ffparams.iparams[il->iatoms[i]];
212 a1 = il->iatoms[i+1];
214 for (j = 1; j < NRAL(ft); j++)
216 cam[j] = atoms->atom[il->iatoms[i+1+j]].m;
219 cam[j] = vsite_m[il->iatoms[i+1+j]];
223 gmx_fatal(FARGS, "In molecule type '%s' %s construction involves atom %d, which is a virtual site of equal or high complexity. This is not supported.",
224 *mtop->moltype[mtop->molblock[mb].type].name,
225 interaction_function[ft].longname,
226 il->iatoms[i+1+j]+1);
233 /* Exact except for ignoring constraints */
234 vsite_m[a1] = (cam[2]*sqr(1-ip->vsite.a) + cam[1]*sqr(ip->vsite.a))/(cam[1]*cam[2]);
237 /* Exact except for ignoring constraints */
238 vsite_m[a1] = (cam[2]*cam[3]*sqr(1-ip->vsite.a-ip->vsite.b) + cam[1]*cam[3]*sqr(ip->vsite.a) + cam[1]*cam[2]*sqr(ip->vsite.b))/(cam[1]*cam[2]*cam[3]);
241 /* Use the mass of the lightest constructing atom.
242 * This is an approximation.
243 * If the distance of the virtual site to the
244 * constructing atom is less than all distances
245 * between constructing atoms, this is a safe
246 * over-estimate of the displacement of the vsite.
247 * This condition holds for all H mass replacement
248 * replacement vsite constructions, except for SP2/3
249 * groups. In SP3 groups one H will have a F_VSITE3
250 * construction, so even there the total drift
251 * estimation shouldn't be far off.
254 vsite_m[a1] = cam[1];
255 for (j = 2; j < NRAL(ft); j++)
257 vsite_m[a1] = min(vsite_m[a1], cam[j]);
259 if (n_nonlin_vsite != NULL)
261 *n_nonlin_vsite += nmol;
269 for (a = 0; a < atoms->nr; a++)
271 at = &atoms->atom[a];
272 /* We consider an atom constrained, #DOF=2, when it is
273 * connected with constraints to one or more atoms with
274 * total mass larger than 1.5 that of the atom itself.
277 at->m, at->type, at->q, con_m[a] > 1.5*at->m, nmol);
286 for (a = 0; a < natt; a++)
288 fprintf(debug, "type %d: m %5.2f t %d q %6.3f con %d n %d\n",
289 a, att[a].mass, att[a].type, att[a].q, att[a].con, att[a].n);
297 static void approx_2dof(real s2, real x,
298 real *shift, real *scale)
300 /* A particle with 1 DOF constrained has 2 DOFs instead of 3.
301 * This code is also used for particles with multiple constraints,
302 * in which case we overestimate the displacement.
303 * The 2DOF distribution is sqrt(pi/2)*erfc(r/(sqrt(2)*s))/(2*s).
304 * We approximate this with scale*Gaussian(s,r+shift),
305 * by matching the distribution value and derivative at x.
306 * This is a tight overestimate for all r>=0 at any s and x.
310 ex = exp(-x*x/(2*s2));
311 er = gmx_erfc(x/sqrt(2*s2));
313 *shift = -x + sqrt(2*s2/M_PI)*ex/er;
314 *scale = 0.5*M_PI*exp(ex*ex/(M_PI*er*er))*er;
317 static real ener_drift(const verletbuf_atomtype_t *att, int natt,
318 const gmx_ffparams_t *ffp,
320 real md_ljd, real md_ljr, real md_el, real dd_el,
322 real rlist, real boxvol)
324 double drift_tot, pot1, pot2, pot;
326 real s2i, s2j, s2, s;
330 double c_exp, c_erfc;
334 /* Loop over the different atom type pairs */
335 for (i = 0; i < natt; i++)
337 s2i = kT_fac/att[i].mass;
340 for (j = i; j < natt; j++)
342 s2j = kT_fac/att[j].mass;
345 /* Note that attractive and repulsive potentials for individual
346 * pairs will partially cancel.
348 /* -dV/dr at the cut-off for LJ + Coulomb */
350 md_ljd*ffp->iparams[ti*ffp->atnr+tj].lj.c6 +
351 md_ljr*ffp->iparams[ti*ffp->atnr+tj].lj.c12 +
352 md_el*att[i].q*att[j].q;
354 /* d2V/dr2 at the cut-off for Coulomb, we neglect LJ */
355 dd = dd_el*att[i].q*att[j].q;
361 /* For constraints: adapt r and scaling for the Gaussian */
365 approx_2dof(s2i, r_buffer*s2i/s2, &sh, &sc);
372 approx_2dof(s2j, r_buffer*s2j/s2, &sh, &sc);
377 /* Exact contribution of an atom pair with Gaussian displacement
378 * with sigma s to the energy drift for a potential with
379 * derivative -md and second derivative dd at the cut-off.
380 * The only catch is that for potentials that change sign
381 * near the cut-off there could be an unlucky compensation
382 * of positive and negative energy drift.
383 * Such potentials are extremely rare though.
385 * Note that pot has unit energy*length, as the linear
386 * atom density still needs to be put in.
388 c_exp = exp(-rsh*rsh/(2*s2))/sqrt(2*M_PI);
389 c_erfc = 0.5*gmx_erfc(rsh/(sqrt(2*s2)));
393 md/2*((rsh*rsh + s2)*c_erfc - rsh*s*c_exp);
395 dd/6*(s*(rsh*rsh + 2*s2)*c_exp - rsh*(rsh*rsh + 3*s2)*c_erfc);
400 fprintf(debug, "n %d %d d s %.3f %.3f con %d md %8.1e dd %8.1e pot1 %8.1e pot2 %8.1e pot %8.1e\n",
401 att[i].n, att[j].n, sqrt(s2i), sqrt(s2j),
402 att[i].con+att[j].con,
403 md, dd, pot1, pot2, pot);
406 /* Multiply by the number of atom pairs */
409 pot *= (double)att[i].n*(att[i].n - 1)/2;
413 pot *= (double)att[i].n*att[j].n;
415 /* We need the line density to get the energy drift of the system.
416 * The effective average r^2 is close to (rlist+sigma)^2.
418 pot *= 4*M_PI*sqr(rlist + s)/boxvol;
420 /* Add the unsigned drift to avoid cancellation of errors */
421 drift_tot += fabs(pot);
428 static real surface_frac(int cluster_size, real particle_distance, real rlist)
432 if (rlist < 0.5*particle_distance)
434 /* We have non overlapping spheres */
438 /* Half the inter-particle distance relative to rlist */
439 d = 0.5*particle_distance/rlist;
441 /* Determine the area of the surface at distance rlist to the closest
442 * particle, relative to surface of a sphere of radius rlist.
443 * The formulas below assume close to cubic cells for the pair search grid,
444 * which the pair search code tries to achieve.
445 * Note that in practice particle distances will not be delta distributed,
446 * but have some spread, often involving shorter distances,
447 * as e.g. O-H bonds in a water molecule. Thus the estimates below will
448 * usually be slightly too high and thus conservative.
450 switch (cluster_size)
453 /* One particle: trivial */
457 /* Two particles: two spheres at fractional distance 2*a */
461 /* We assume a perfect, symmetric tetrahedron geometry.
462 * The surface around a tetrahedron is too complex for a full
463 * analytical solution, so we use a Taylor expansion.
465 area_rel = (1.0 + 1/M_PI*(6*acos(1/sqrt(3))*d +
469 83.0/756.0*d*d*d*d*d*d)));
472 gmx_incons("surface_frac called with unsupported cluster_size");
476 return area_rel/cluster_size;
479 void calc_verlet_buffer_size(const gmx_mtop_t *mtop, real boxvol,
480 const t_inputrec *ir, real drift_target,
481 const verletbuf_list_setup_t *list_setup,
488 real particle_distance;
489 real nb_clust_frac_pairs_not_in_list_at_cutoff;
491 verletbuf_atomtype_t *att = NULL;
494 real md_ljd, md_ljr, md_el, dd_el;
496 real kT_fac, mass_min;
501 /* Resolution of the buffer size */
504 env = getenv("GMX_VERLET_BUFFER_RES");
507 sscanf(env, "%lf", &resolution);
510 /* In an atom wise pair-list there would be no pairs in the list
511 * beyond the pair-list cut-off.
512 * However, we use a pair-list of groups vs groups of atoms.
513 * For groups of 4 atoms, the parallelism of SSE instructions, only
514 * 10% of the atoms pairs are not in the list just beyond the cut-off.
515 * As this percentage increases slowly compared to the decrease of the
516 * Gaussian displacement distribution over this range, we can simply
517 * reduce the drift by this fraction.
518 * For larger groups, e.g. of 8 atoms, this fraction will be lower,
519 * so then buffer size will be on the conservative (large) side.
521 * Note that the formulas used here do not take into account
522 * cancellation of errors which could occur by missing both
523 * attractive and repulsive interactions.
525 * The only major assumption is homogeneous particle distribution.
526 * For an inhomogeneous system, such as a liquid-vapor system,
527 * the buffer will be underestimated. The actual energy drift
528 * will be higher by the factor: local/homogeneous particle density.
530 * The results of this estimate have been checked againt simulations.
531 * In most cases the real drift differs by less than a factor 2.
534 /* Worst case assumption: HCP packing of particles gives largest distance */
535 particle_distance = pow(boxvol*sqrt(2)/mtop->natoms, 1.0/3.0);
537 get_verlet_buffer_atomtypes(mtop, &att, &natt, n_nonlin_vsite);
538 assert(att != NULL && natt >= 0);
542 fprintf(debug, "particle distance assuming HCP packing: %f nm\n",
544 fprintf(debug, "energy drift atom types: %d\n", natt);
547 reppow = mtop->ffparams.reppow;
550 if (ir->vdwtype == evdwCUT)
552 /* -dV/dr of -r^-6 and r^-repporw */
553 md_ljd = -6*pow(ir->rvdw, -7.0);
554 md_ljr = reppow*pow(ir->rvdw, -(reppow+1));
555 /* The contribution of the second derivative is negligible */
559 gmx_fatal(FARGS, "Energy drift calculation is only implemented for plain cut-off Lennard-Jones interactions");
562 elfac = ONE_4PI_EPS0/ir->epsilon_r;
564 /* Determine md=-dV/dr and dd=d^2V/dr^2 */
567 if (ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype))
571 if (ir->coulombtype == eelCUT)
578 eps_rf = ir->epsilon_rf/ir->epsilon_r;
581 k_rf = pow(ir->rcoulomb, -3.0)*(eps_rf - ir->epsilon_r)/(2*eps_rf + ir->epsilon_r);
585 /* epsilon_rf = infinity */
586 k_rf = 0.5*pow(ir->rcoulomb, -3.0);
592 md_el = elfac*(pow(ir->rcoulomb, -2.0) - 2*k_rf*ir->rcoulomb);
594 dd_el = elfac*(2*pow(ir->rcoulomb, -3.0) + 2*k_rf);
596 else if (EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD)
600 b = calc_ewaldcoeff(ir->rcoulomb, ir->ewald_rtol);
603 md_el = elfac*(b*exp(-br*br)*M_2_SQRTPI/rc + gmx_erfc(br)/(rc*rc));
604 dd_el = elfac/(rc*rc)*(2*b*(1 + br*br)*exp(-br*br)*M_2_SQRTPI + 2*gmx_erfc(br)/rc);
608 gmx_fatal(FARGS, "Energy drift calculation is only implemented for Reaction-Field and Ewald electrostatics");
611 /* Determine the variance of the atomic displacement
612 * over nstlist-1 steps: kT_fac
613 * For inertial dynamics (not Brownian dynamics) the mass factor
614 * is not included in kT_fac, it is added later.
618 /* Get the displacement distribution from the random component only.
619 * With accurate integration the systematic (force) displacement
620 * should be negligible (unless nstlist is extremely large, which
621 * you wouldn't do anyhow).
623 kT_fac = 2*BOLTZ*ir->opts.ref_t[0]*(ir->nstlist-1)*ir->delta_t;
626 /* This is directly sigma^2 of the displacement */
627 kT_fac /= ir->bd_fric;
629 /* Set the masses to 1 as kT_fac is the full sigma^2,
630 * but we divide by m in ener_drift().
632 for (i = 0; i < natt; i++)
641 /* Per group tau_t is not implemented yet, use the maximum */
642 tau_t = ir->opts.tau_t[0];
643 for (i = 1; i < ir->opts.ngtc; i++)
645 tau_t = max(tau_t, ir->opts.tau_t[i]);
649 /* This kT_fac needs to be divided by the mass to get sigma^2 */
654 kT_fac = BOLTZ*ir->opts.ref_t[0]*sqr((ir->nstlist-1)*ir->delta_t);
657 mass_min = att[0].mass;
658 for (i = 1; i < natt; i++)
660 mass_min = min(mass_min, att[i].mass);
665 fprintf(debug, "md_ljd %e md_ljr %e\n", md_ljd, md_ljr);
666 fprintf(debug, "md_el %e dd_el %e\n", md_el, dd_el);
667 fprintf(debug, "sqrt(kT_fac) %f\n", sqrt(kT_fac));
668 fprintf(debug, "mass_min %f\n", mass_min);
671 /* Search using bisection */
673 /* The drift will be neglible at 5 times the max sigma */
674 ib1 = (int)(5*2*sqrt(kT_fac/mass_min)/resolution) + 1;
675 while (ib1 - ib0 > 1)
679 rl = max(ir->rvdw, ir->rcoulomb) + rb;
681 /* Calculate the average energy drift at the last step
682 * of the nstlist steps at which the pair-list is used.
684 drift = ener_drift(att, natt, &mtop->ffparams,
686 md_ljd, md_ljr, md_el, dd_el, rb,
689 /* Correct for the fact that we are using a Ni x Nj particle pair list
690 * and not a 1 x 1 particle pair list. This reduces the drift.
692 /* We don't have a formula for 8 (yet), use 4 which is conservative */
693 nb_clust_frac_pairs_not_in_list_at_cutoff =
694 surface_frac(min(list_setup->cluster_size_i, 4),
695 particle_distance, rl)*
696 surface_frac(min(list_setup->cluster_size_j, 4),
697 particle_distance, rl);
698 drift *= nb_clust_frac_pairs_not_in_list_at_cutoff;
700 /* Convert the drift to drift per unit time per atom */
701 drift /= ir->nstlist*ir->delta_t*mtop->natoms;
705 fprintf(debug, "ib %3d %3d %3d rb %.3f %dx%d fac %.3f drift %f\n",
707 list_setup->cluster_size_i, list_setup->cluster_size_j,
708 nb_clust_frac_pairs_not_in_list_at_cutoff,
712 if (fabs(drift) > drift_target)
724 *rlist = max(ir->rvdw, ir->rcoulomb) + ib1*resolution;