2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2012,2013,2014,2015, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
37 #include "calc_verletbuf.h"
46 #include "gromacs/math/calculate-ewald-splitting-coefficient.h"
47 #include "gromacs/math/functions.h"
48 #include "gromacs/math/units.h"
49 #include "gromacs/math/vec.h"
50 #include "gromacs/mdlib/nb_verlet.h"
51 #include "gromacs/mdlib/nbnxn_simd.h"
52 #include "gromacs/mdlib/nbnxn_util.h"
53 #include "gromacs/mdtypes/inputrec.h"
54 #include "gromacs/mdtypes/md_enums.h"
55 #include "gromacs/topology/ifunc.h"
56 #include "gromacs/topology/topology.h"
57 #include "gromacs/utility/fatalerror.h"
58 #include "gromacs/utility/smalloc.h"
60 /* The code in this file estimates a pairlist buffer length
61 * given a target energy drift per atom per picosecond.
62 * This is done by estimating the drift given a buffer length.
63 * Ideally we would like to have a tight overestimate of the drift,
64 * but that can be difficult to achieve.
66 * Significant approximations used:
68 * Uniform particle density. UNDERESTIMATES the drift by rho_global/rho_local.
70 * Interactions don't affect particle motion. OVERESTIMATES the drift on longer
71 * time scales. This approximation probably introduces the largest errors.
73 * Only take one constraint per particle into account: OVERESTIMATES the drift.
75 * For rotating constraints assume the same functional shape for time scales
76 * where the constraints rotate significantly as the exact expression for
77 * short time scales. OVERESTIMATES the drift on long time scales.
79 * For non-linear virtual sites use the mass of the lightest constructing atom
80 * to determine the displacement. OVER/UNDERESTIMATES the drift, depending on
81 * the geometry and masses of constructing atoms.
83 * Note that the formulas for normal atoms and linear virtual sites are exact,
84 * apart from the first two approximations.
86 * Note that apart from the effect of the above approximations, the actual
87 * drift of the total energy of a system can be order of magnitude smaller
88 * due to cancellation of positive and negative drift for different pairs.
92 /* Struct for unique atom type for calculating the energy drift.
93 * The atom displacement depends on mass and constraints.
94 * The energy jump for given distance depend on LJ type and q.
99 int type; /* type (used for LJ parameters) */
101 gmx_bool bConstr; /* constrained, if TRUE, use #DOF=2 iso 3 */
102 real con_mass; /* mass of heaviest atom connected by constraints */
103 real con_len; /* constraint length to the heaviest atom */
104 } atom_nonbonded_kinetic_prop_t;
106 /* Struct for unique atom type for calculating the energy drift.
107 * The atom displacement depends on mass and constraints.
108 * The energy jump for given distance depend on LJ type and q.
112 atom_nonbonded_kinetic_prop_t prop; /* non-bonded and kinetic atom prop. */
113 int n; /* #atoms of this type in the system */
114 } verletbuf_atomtype_t;
116 void verletbuf_get_list_setup(gmx_bool gmx_unused bSIMD,
118 verletbuf_list_setup_t *list_setup)
120 /* When calling this function we often don't know which kernel type we
121 * are going to use. W choose the kernel type with the smallest possible
122 * i- and j-cluster sizes, so we potentially overestimate, but never
123 * underestimate, the buffer drift.
124 * Note that the current buffer estimation code only handles clusters
125 * of size 1, 2 or 4, so for 4x8 or 8x8 we use the estimate for 4x4.
130 /* The CUDA kernels split the j-clusters in two halves */
131 list_setup->cluster_size_i = nbnxn_kernel_to_cluster_i_size(nbnxnk8x8x8_GPU);
132 list_setup->cluster_size_j = nbnxn_kernel_to_cluster_j_size(nbnxnk8x8x8_GPU)/2;
138 kernel_type = nbnxnk4x4_PlainC;
143 #ifdef GMX_NBNXN_SIMD_2XNN
144 /* We use the smallest cluster size to be on the safe side */
145 kernel_type = nbnxnk4xN_SIMD_2xNN;
147 kernel_type = nbnxnk4xN_SIMD_4xN;
152 list_setup->cluster_size_i = nbnxn_kernel_to_cluster_i_size(kernel_type);
153 list_setup->cluster_size_j = nbnxn_kernel_to_cluster_j_size(kernel_type);
158 atom_nonbonded_kinetic_prop_equal(const atom_nonbonded_kinetic_prop_t *prop1,
159 const atom_nonbonded_kinetic_prop_t *prop2)
161 return (prop1->mass == prop2->mass &&
162 prop1->type == prop2->type &&
163 prop1->q == prop2->q &&
164 prop1->bConstr == prop2->bConstr &&
165 prop1->con_mass == prop2->con_mass &&
166 prop1->con_len == prop2->con_len);
169 static void add_at(verletbuf_atomtype_t **att_p, int *natt_p,
170 const atom_nonbonded_kinetic_prop_t *prop,
173 verletbuf_atomtype_t *att;
178 /* Ignore massless particles */
186 while (i < natt && !atom_nonbonded_kinetic_prop_equal(prop, &att[i].prop))
198 srenew(*att_p, *natt_p);
199 (*att_p)[i].prop = *prop;
200 (*att_p)[i].n = nmol;
204 static void get_vsite_masses(const gmx_moltype_t *moltype,
205 const gmx_ffparams_t *ffparams,
214 /* Check for virtual sites, determine mass from constructing atoms */
215 for (ft = 0; ft < F_NRE; ft++)
219 il = &moltype->ilist[ft];
221 for (i = 0; i < il->nr; i += 1+NRAL(ft))
224 real inv_mass, coeff, m_aj;
227 ip = &ffparams->iparams[il->iatoms[i]];
229 a1 = il->iatoms[i+1];
233 /* Only vsiten can have more than four
234 constructing atoms, so NRAL(ft) <= 5 */
237 const int maxj = NRAL(ft);
241 for (j = 1; j < maxj; j++)
243 cam[j] = moltype->atoms.atom[il->iatoms[i+1+j]].m;
246 cam[j] = vsite_m[il->iatoms[i+1+j]];
250 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.",
252 interaction_function[ft].longname,
253 il->iatoms[i+1+j]+1);
261 vsite_m[a1] = (cam[1]*cam[2])/(cam[2]*gmx::square(1-ip->vsite.a) + cam[1]*gmx::square(ip->vsite.a));
265 vsite_m[a1] = (cam[1]*cam[2]*cam[3])/(cam[2]*cam[3]*gmx::square(1-ip->vsite.a-ip->vsite.b) + cam[1]*cam[3]*gmx::square(ip->vsite.a) + cam[1]*cam[2]*gmx::square(ip->vsite.b));
268 gmx_incons("Invalid vsite type");
271 /* Use the mass of the lightest constructing atom.
272 * This is an approximation.
273 * If the distance of the virtual site to the
274 * constructing atom is less than all distances
275 * between constructing atoms, this is a safe
276 * over-estimate of the displacement of the vsite.
277 * This condition holds for all H mass replacement
278 * vsite constructions, except for SP2/3 groups.
279 * In SP3 groups one H will have a F_VSITE3
280 * construction, so even there the total drift
281 * estimate shouldn't be far off.
283 vsite_m[a1] = cam[1];
284 for (j = 2; j < maxj; j++)
286 vsite_m[a1] = std::min(vsite_m[a1], cam[j]);
299 for (j = 0; j < 3*ffparams->iparams[il->iatoms[i]].vsiten.n; j += 3)
301 aj = il->iatoms[i+j+2];
302 coeff = ffparams->iparams[il->iatoms[i+j]].vsiten.a;
303 if (moltype->atoms.atom[aj].ptype == eptVSite)
309 m_aj = moltype->atoms.atom[aj].m;
313 gmx_incons("The mass of a vsiten constructing atom is <= 0");
315 inv_mass += coeff*coeff/m_aj;
317 vsite_m[a1] = 1/inv_mass;
318 /* Correct for loop increment of i */
319 i += j - 1 - NRAL(ft);
323 fprintf(debug, "atom %4d %-20s mass %6.3f\n",
324 a1, interaction_function[ft].longname, vsite_m[a1]);
331 static void get_verlet_buffer_atomtypes(const gmx_mtop_t *mtop,
332 verletbuf_atomtype_t **att_p,
336 verletbuf_atomtype_t *att;
338 int mb, nmol, ft, i, a1, a2, a3, a;
339 const t_atoms *atoms;
342 atom_nonbonded_kinetic_prop_t *prop;
344 int n_nonlin_vsite_mol;
349 if (n_nonlin_vsite != NULL)
354 for (mb = 0; mb < mtop->nmolblock; mb++)
356 nmol = mtop->molblock[mb].nmol;
358 atoms = &mtop->moltype[mtop->molblock[mb].type].atoms;
360 /* Check for constraints, as they affect the kinetic energy.
361 * For virtual sites we need the masses and geometry of
362 * the constructing atoms to determine their velocity distribution.
364 snew(prop, atoms->nr);
365 snew(vsite_m, atoms->nr);
367 for (ft = F_CONSTR; ft <= F_CONSTRNC; ft++)
369 il = &mtop->moltype[mtop->molblock[mb].type].ilist[ft];
371 for (i = 0; i < il->nr; i += 1+NRAL(ft))
373 ip = &mtop->ffparams.iparams[il->iatoms[i]];
374 a1 = il->iatoms[i+1];
375 a2 = il->iatoms[i+2];
376 if (atoms->atom[a2].m > prop[a1].con_mass)
378 prop[a1].con_mass = atoms->atom[a2].m;
379 prop[a1].con_len = ip->constr.dA;
381 if (atoms->atom[a1].m > prop[a2].con_mass)
383 prop[a2].con_mass = atoms->atom[a1].m;
384 prop[a2].con_len = ip->constr.dA;
389 il = &mtop->moltype[mtop->molblock[mb].type].ilist[F_SETTLE];
391 for (i = 0; i < il->nr; i += 1+NRAL(F_SETTLE))
393 ip = &mtop->ffparams.iparams[il->iatoms[i]];
394 a1 = il->iatoms[i+1];
395 a2 = il->iatoms[i+2];
396 a3 = il->iatoms[i+3];
397 /* Usually the mass of a1 (usually oxygen) is larger than a2/a3.
398 * If this is not the case, we overestimate the displacement,
399 * which leads to a larger buffer (ok since this is an exotic case).
401 prop[a1].con_mass = atoms->atom[a2].m;
402 prop[a1].con_len = ip->settle.doh;
404 prop[a2].con_mass = atoms->atom[a1].m;
405 prop[a2].con_len = ip->settle.doh;
407 prop[a3].con_mass = atoms->atom[a1].m;
408 prop[a3].con_len = ip->settle.doh;
411 get_vsite_masses(&mtop->moltype[mtop->molblock[mb].type],
414 &n_nonlin_vsite_mol);
415 if (n_nonlin_vsite != NULL)
417 *n_nonlin_vsite += nmol*n_nonlin_vsite_mol;
420 for (a = 0; a < atoms->nr; a++)
422 if (atoms->atom[a].ptype == eptVSite)
424 prop[a].mass = vsite_m[a];
428 prop[a].mass = atoms->atom[a].m;
430 prop[a].type = atoms->atom[a].type;
431 prop[a].q = atoms->atom[a].q;
432 /* We consider an atom constrained, #DOF=2, when it is
433 * connected with constraints to (at least one) atom with
434 * a mass of more than 0.4x its own mass. This is not a critical
435 * parameter, since with roughly equal masses the unconstrained
436 * and constrained displacement will not differ much (and both
437 * overestimate the displacement).
439 prop[a].bConstr = (prop[a].con_mass > 0.4*prop[a].mass);
441 add_at(&att, &natt, &prop[a], nmol);
450 for (a = 0; a < natt; a++)
452 fprintf(debug, "type %d: m %5.2f t %d q %6.3f con %d con_m %5.3f con_l %5.3f n %d\n",
453 a, att[a].prop.mass, att[a].prop.type, att[a].prop.q,
454 att[a].prop.bConstr, att[a].prop.con_mass, att[a].prop.con_len,
463 /* This function computes two components of the estimate of the variance
464 * in the displacement of one atom in a system of two constrained atoms.
465 * Returns in sigma2_2d the variance due to rotation of the constrained
466 * atom around the atom to which it constrained.
467 * Returns in sigma2_3d the variance due to displacement of the COM
468 * of the whole system of the two constrained atoms.
470 * Note that we only take a single constraint (the one to the heaviest atom)
471 * into account. If an atom has multiple constraints, this will result in
472 * an overestimate of the displacement, which gives a larger drift and buffer.
474 static void constrained_atom_sigma2(real kT_fac,
475 const atom_nonbonded_kinetic_prop_t *prop,
484 /* Here we decompose the motion of a constrained atom into two
485 * components: rotation around the COM and translation of the COM.
488 /* Determine the variance for the displacement of the rotational mode */
489 sigma2_rot = kT_fac/(prop->mass*(prop->mass + prop->con_mass)/prop->con_mass);
491 /* The distance from the atom to the COM, i.e. the rotational arm */
492 com_dist = prop->con_len*prop->con_mass/(prop->mass + prop->con_mass);
494 /* The variance relative to the arm */
495 sigma2_rel = sigma2_rot/(com_dist*com_dist);
496 /* At 6 the scaling formula has slope 0,
497 * so we keep sigma2_2d constant after that.
501 /* A constrained atom rotates around the atom it is constrained to.
502 * This results in a smaller linear displacement than for a free atom.
503 * For a perfectly circular displacement, this lowers the displacement
504 * by: 1/arcsin(arc_length)
505 * and arcsin(x) = 1 + x^2/6 + ...
506 * For sigma2_rel<<1 the displacement distribution is erfc
507 * (exact formula is provided below). For larger sigma, it is clear
508 * that the displacement can't be larger than 2*com_dist.
509 * It turns out that the distribution becomes nearly uniform.
510 * For intermediate sigma2_rel, scaling down sigma with the third
511 * order expansion of arcsin with argument sigma_rel turns out
512 * to give a very good approximation of the distribution and variance.
513 * Even for larger values, the variance is only slightly overestimated.
514 * Note that the most relevant displacements are in the long tail.
515 * This rotation approximation always overestimates the tail (which
516 * runs to infinity, whereas it should be <= 2*com_dist).
517 * Thus we always overestimate the drift and the buffer size.
519 scale = 1/(1 + sigma2_rel/6);
520 *sigma2_2d = sigma2_rot*scale*scale;
524 /* sigma_2d is set to the maximum given by the scaling above.
525 * For large sigma2 the real displacement distribution is close
526 * to uniform over -2*con_len to 2*com_dist.
527 * Our erfc with sigma_2d=sqrt(1.5)*com_dist (which means the sigma
528 * of the erfc output distribution is con_dist) overestimates
529 * the variance and additionally has a long tail. This means
530 * we have a (safe) overestimation of the drift.
532 *sigma2_2d = 1.5*com_dist*com_dist;
535 /* The constrained atom also moves (in 3D) with the COM of both atoms */
536 *sigma2_3d = kT_fac/(prop->mass + prop->con_mass);
539 static void get_atom_sigma2(real kT_fac,
540 const atom_nonbonded_kinetic_prop_t *prop,
546 /* Complicated constraint calculation in a separate function */
547 constrained_atom_sigma2(kT_fac, prop, sigma2_2d, sigma2_3d);
551 /* Unconstrained atom: trivial */
553 *sigma2_3d = kT_fac/prop->mass;
557 static void approx_2dof(real s2, real x, real *shift, real *scale)
559 /* A particle with 1 DOF constrained has 2 DOFs instead of 3.
560 * This code is also used for particles with multiple constraints,
561 * in which case we overestimate the displacement.
562 * The 2DOF distribution is sqrt(pi/2)*erfc(r/(sqrt(2)*s))/(2*s).
563 * We approximate this with scale*Gaussian(s,r+shift),
564 * by matching the distribution value and derivative at x.
565 * This is a tight overestimate for all r>=0 at any s and x.
569 ex = exp(-x*x/(2*s2));
570 er = std::erfc(x/std::sqrt(2*s2));
572 *shift = -x + std::sqrt(2*s2/M_PI)*ex/er;
573 *scale = 0.5*M_PI*std::exp(ex*ex/(M_PI*er*er))*er;
576 static real ener_drift(const verletbuf_atomtype_t *att, int natt,
577 const gmx_ffparams_t *ffp,
579 real md1_ljd, real d2_ljd, real md3_ljd,
580 real md1_ljr, real d2_ljr, real md3_ljr,
581 real md1_el, real d2_el,
583 real rlist, real boxvol)
585 /* Erfc(8)=1e-29, use this limit so we have some space for arithmetic
586 * on the result when using float precision.
588 const real erfc_arg_max = 8.0;
590 double drift_tot, pot1, pot2, pot3, pot;
592 real s2i_2d, s2i_3d, s2j_2d, s2j_3d, s2, s;
595 real sc_fac, rsh, rsh2;
596 double c_exp, c_erfc;
600 /* Loop over the different atom type pairs */
601 for (i = 0; i < natt; i++)
603 get_atom_sigma2(kT_fac, &att[i].prop, &s2i_2d, &s2i_3d);
604 ti = att[i].prop.type;
606 for (j = i; j < natt; j++)
608 get_atom_sigma2(kT_fac, &att[j].prop, &s2j_2d, &s2j_3d);
609 tj = att[j].prop.type;
611 /* Add up the up to four independent variances */
612 s2 = s2i_2d + s2i_3d + s2j_2d + s2j_3d;
614 /* Note that attractive and repulsive potentials for individual
615 * pairs will partially cancel.
617 /* -dV/dr at the cut-off for LJ + Coulomb */
619 md1_ljd*ffp->iparams[ti*ffp->atnr+tj].lj.c6 +
620 md1_ljr*ffp->iparams[ti*ffp->atnr+tj].lj.c12 +
621 md1_el*att[i].prop.q*att[j].prop.q;
623 /* d2V/dr2 at the cut-off for LJ + Coulomb */
625 d2_ljd*ffp->iparams[ti*ffp->atnr+tj].lj.c6 +
626 d2_ljr*ffp->iparams[ti*ffp->atnr+tj].lj.c12 +
627 d2_el*att[i].prop.q*att[j].prop.q;
629 /* -d3V/dr3 at the cut-off for LJ, we neglect Coulomb */
631 md3_ljd*ffp->iparams[ti*ffp->atnr+tj].lj.c6 +
632 md3_ljr*ffp->iparams[ti*ffp->atnr+tj].lj.c12;
637 if (rsh*rsh > 2*s2*erfc_arg_max*erfc_arg_max)
639 /* Erfc might run out of float and become 0, somewhat before
640 * c_exp becomes 0. To avoid this and to avoid NaN in
641 * approx_2dof, we set both c_expc and c_erfc to zero.
642 * In any relevant case this has no effect on the results,
643 * since c_exp < 6e-29, so the displacement is completely
644 * negligible for such atom pairs (and an overestimate).
645 * In nearly all use cases, there will be other atom
646 * pairs that contribute much more to the total, so zeroing
647 * this particular contribution has no effect at all.
654 /* For constraints: adapt r and scaling for the Gaussian */
655 if (att[i].prop.bConstr)
659 approx_2dof(s2i_2d, r_buffer*s2i_2d/s2, &sh, &sc);
663 if (att[j].prop.bConstr)
667 approx_2dof(s2j_2d, r_buffer*s2j_2d/s2, &sh, &sc);
672 /* Exact contribution of an atom pair with Gaussian displacement
673 * with sigma s to the energy drift for a potential with
674 * derivative -md and second derivative dd at the cut-off.
675 * The only catch is that for potentials that change sign
676 * near the cut-off there could be an unlucky compensation
677 * of positive and negative energy drift.
678 * Such potentials are extremely rare though.
680 * Note that pot has unit energy*length, as the linear
681 * atom density still needs to be put in.
683 c_exp = std::exp(-rsh*rsh/(2*s2))/std::sqrt(2*M_PI);
684 c_erfc = 0.5*std::erfc(rsh/(std::sqrt(2*s2)));
690 md1/2*((rsh2 + s2)*c_erfc - rsh*s*c_exp);
692 d2/6*(s*(rsh2 + 2*s2)*c_exp - rsh*(rsh2 + 3*s2)*c_erfc);
694 md3/24*((rsh2*rsh2 + 6*rsh2*s2 + 3*s2*s2)*c_erfc - rsh*s*(rsh2 + 5*s2)*c_exp);
695 pot = pot1 + pot2 + pot3;
699 fprintf(debug, "n %d %d d s %.3f %.3f %.3f %.3f con %d -d1 %8.1e d2 %8.1e -d3 %8.1e pot1 %8.1e pot2 %8.1e pot3 %8.1e pot %8.1e\n",
701 std::sqrt(s2i_2d), std::sqrt(s2i_3d),
702 std::sqrt(s2j_2d), std::sqrt(s2j_3d),
703 att[i].prop.bConstr+att[j].prop.bConstr,
705 pot1, pot2, pot3, pot);
708 /* Multiply by the number of atom pairs */
711 pot *= (double)att[i].n*(att[i].n - 1)/2;
715 pot *= (double)att[i].n*att[j].n;
717 /* We need the line density to get the energy drift of the system.
718 * The effective average r^2 is close to (rlist+sigma)^2.
720 pot *= 4*M_PI*gmx::square(rlist + s)/boxvol;
722 /* Add the unsigned drift to avoid cancellation of errors */
723 drift_tot += std::abs(pot);
730 static real surface_frac(int cluster_size, real particle_distance, real rlist)
734 if (rlist < 0.5*particle_distance)
736 /* We have non overlapping spheres */
740 /* Half the inter-particle distance relative to rlist */
741 d = 0.5*particle_distance/rlist;
743 /* Determine the area of the surface at distance rlist to the closest
744 * particle, relative to surface of a sphere of radius rlist.
745 * The formulas below assume close to cubic cells for the pair search grid,
746 * which the pair search code tries to achieve.
747 * Note that in practice particle distances will not be delta distributed,
748 * but have some spread, often involving shorter distances,
749 * as e.g. O-H bonds in a water molecule. Thus the estimates below will
750 * usually be slightly too high and thus conservative.
752 switch (cluster_size)
755 /* One particle: trivial */
759 /* Two particles: two spheres at fractional distance 2*a */
763 /* We assume a perfect, symmetric tetrahedron geometry.
764 * The surface around a tetrahedron is too complex for a full
765 * analytical solution, so we use a Taylor expansion.
767 area_rel = (1.0 + 1/M_PI*(6*std::acos(1/std::sqrt(3))*d +
768 std::sqrt(3)*d*d*(1.0 +
771 83.0/756.0*d*d*d*d*d*d)));
774 gmx_incons("surface_frac called with unsupported cluster_size");
778 return area_rel/cluster_size;
781 /* Returns the negative of the third derivative of a potential r^-p
782 * with a force-switch function, evaluated at the cut-off rc.
784 static real md3_force_switch(real p, real rswitch, real rc)
786 /* The switched force function is:
787 * p*r^-(p+1) + a*(r - rswitch)^2 + b*(r - rswitch)^3
790 real md3_pot, md3_sw;
792 a = -((p + 4)*rc - (p + 1)*rswitch)/(pow(rc, p+2)*gmx::square(rc-rswitch));
793 b = ((p + 3)*rc - (p + 1)*rswitch)/(pow(rc, p+2)*gmx::power3(rc-rswitch));
795 md3_pot = (p + 2)*(p + 1)*p*pow(rc, p+3);
796 md3_sw = 2*a + 6*b*(rc - rswitch);
798 return md3_pot + md3_sw;
801 void calc_verlet_buffer_size(const gmx_mtop_t *mtop, real boxvol,
802 const t_inputrec *ir,
803 real reference_temperature,
804 const verletbuf_list_setup_t *list_setup,
811 real particle_distance;
812 real nb_clust_frac_pairs_not_in_list_at_cutoff;
814 verletbuf_atomtype_t *att = NULL;
817 real md1_ljd, d2_ljd, md3_ljd;
818 real md1_ljr, d2_ljr, md3_ljr;
821 real kT_fac, mass_min;
826 if (reference_temperature < 0)
828 if (EI_MD(ir->eI) && ir->etc == etcNO)
830 /* This case should be handled outside calc_verlet_buffer_size */
831 gmx_incons("calc_verlet_buffer_size called with an NVE ensemble and reference_temperature < 0");
834 /* We use the maximum temperature with multiple T-coupl groups.
835 * We could use a per particle temperature, but since particles
836 * interact, this might underestimate the buffer size.
838 reference_temperature = 0;
839 for (i = 0; i < ir->opts.ngtc; i++)
841 if (ir->opts.tau_t[i] >= 0)
843 reference_temperature = std::max(reference_temperature,
849 /* Resolution of the buffer size */
852 env = getenv("GMX_VERLET_BUFFER_RES");
855 sscanf(env, "%lf", &resolution);
858 /* In an atom wise pair-list there would be no pairs in the list
859 * beyond the pair-list cut-off.
860 * However, we use a pair-list of groups vs groups of atoms.
861 * For groups of 4 atoms, the parallelism of SSE instructions, only
862 * 10% of the atoms pairs are not in the list just beyond the cut-off.
863 * As this percentage increases slowly compared to the decrease of the
864 * Gaussian displacement distribution over this range, we can simply
865 * reduce the drift by this fraction.
866 * For larger groups, e.g. of 8 atoms, this fraction will be lower,
867 * so then buffer size will be on the conservative (large) side.
869 * Note that the formulas used here do not take into account
870 * cancellation of errors which could occur by missing both
871 * attractive and repulsive interactions.
873 * The only major assumption is homogeneous particle distribution.
874 * For an inhomogeneous system, such as a liquid-vapor system,
875 * the buffer will be underestimated. The actual energy drift
876 * will be higher by the factor: local/homogeneous particle density.
878 * The results of this estimate have been checked againt simulations.
879 * In most cases the real drift differs by less than a factor 2.
882 /* Worst case assumption: HCP packing of particles gives largest distance */
883 particle_distance = std::cbrt(boxvol*std::sqrt(2)/mtop->natoms);
885 get_verlet_buffer_atomtypes(mtop, &att, &natt, n_nonlin_vsite);
886 assert(att != NULL && natt >= 0);
890 fprintf(debug, "particle distance assuming HCP packing: %f nm\n",
892 fprintf(debug, "energy drift atom types: %d\n", natt);
895 reppow = mtop->ffparams.reppow;
902 if (ir->vdwtype == evdwCUT)
904 real sw_range, md3_pswf;
906 switch (ir->vdw_modifier)
909 case eintmodPOTSHIFT:
910 /* -dV/dr of -r^-6 and r^-reppow */
911 md1_ljd = -6/(ir->rvdw*gmx::power6(ir->rvdw));
912 md1_ljr = reppow*pow(ir->rvdw, -(reppow+1));
913 /* The contribution of the higher derivatives is negligible */
915 case eintmodFORCESWITCH:
916 /* At the cut-off: V=V'=V''=0, so we use only V''' */
917 md3_ljd = -md3_force_switch(6.0, ir->rvdw_switch, ir->rvdw);
918 md3_ljr = md3_force_switch(reppow, ir->rvdw_switch, ir->rvdw);
920 case eintmodPOTSWITCH:
921 /* At the cut-off: V=V'=V''=0.
922 * V''' is given by the original potential times
923 * the third derivative of the switch function.
925 sw_range = ir->rvdw - ir->rvdw_switch;
926 md3_pswf = 60.0/gmx::power3(sw_range);
928 md3_ljd = -1.0/gmx::power6(ir->rvdw)*md3_pswf;
929 md3_ljr = pow(ir->rvdw, -reppow)*md3_pswf;
932 gmx_incons("Unimplemented VdW modifier");
935 else if (EVDW_PME(ir->vdwtype))
937 real b, r, br, br2, br4, br6;
938 b = calc_ewaldcoeff_lj(ir->rvdw, ir->ewald_rtol_lj);
944 /* -dV/dr of g(br)*r^-6 [where g(x) = exp(-x^2)(1+x^2+x^4/2), see LJ-PME equations in manual] and r^-reppow */
945 md1_ljd = -std::exp(-br2)*(br6 + 3.0*br4 + 6.0*br2 + 6.0)/(r*gmx::power6(r));
946 md1_ljr = reppow*std::pow(r, -(reppow+1));
947 /* The contribution of the higher derivatives is negligible */
951 gmx_fatal(FARGS, "Energy drift calculation is only implemented for plain cut-off Lennard-Jones interactions");
954 elfac = ONE_4PI_EPS0/ir->epsilon_r;
956 /* Determine md=-dV/dr and dd=d^2V/dr^2 */
958 if (ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype))
962 if (ir->coulombtype == eelCUT)
969 eps_rf = ir->epsilon_rf/ir->epsilon_r;
972 k_rf = (eps_rf - ir->epsilon_r)/( gmx::power3(ir->rcoulomb) * (2*eps_rf + ir->epsilon_r) );
976 /* epsilon_rf = infinity */
977 k_rf = 0.5/gmx::power3(ir->rcoulomb);
983 md1_el = elfac*(1.0/gmx::square(ir->rcoulomb) - 2*k_rf*ir->rcoulomb);
985 d2_el = elfac*(2.0/gmx::power3(ir->rcoulomb) + 2*k_rf);
987 else if (EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD)
991 b = calc_ewaldcoeff_q(ir->rcoulomb, ir->ewald_rtol);
994 md1_el = elfac*(b*exp(-br*br)*M_2_SQRTPI/rc + std::erfc(br)/(rc*rc));
995 d2_el = elfac/(rc*rc)*(2*b*(1 + br*br)*exp(-br*br)*M_2_SQRTPI + 2*std::erfc(br)/rc);
999 gmx_fatal(FARGS, "Energy drift calculation is only implemented for Reaction-Field and Ewald electrostatics");
1002 /* Determine the variance of the atomic displacement
1003 * over nstlist-1 steps: kT_fac
1004 * For inertial dynamics (not Brownian dynamics) the mass factor
1005 * is not included in kT_fac, it is added later.
1009 /* Get the displacement distribution from the random component only.
1010 * With accurate integration the systematic (force) displacement
1011 * should be negligible (unless nstlist is extremely large, which
1012 * you wouldn't do anyhow).
1014 kT_fac = 2*BOLTZ*reference_temperature*(ir->nstlist-1)*ir->delta_t;
1015 if (ir->bd_fric > 0)
1017 /* This is directly sigma^2 of the displacement */
1018 kT_fac /= ir->bd_fric;
1020 /* Set the masses to 1 as kT_fac is the full sigma^2,
1021 * but we divide by m in ener_drift().
1023 for (i = 0; i < natt; i++)
1025 att[i].prop.mass = 1;
1032 /* Per group tau_t is not implemented yet, use the maximum */
1033 tau_t = ir->opts.tau_t[0];
1034 for (i = 1; i < ir->opts.ngtc; i++)
1036 tau_t = std::max(tau_t, ir->opts.tau_t[i]);
1040 /* This kT_fac needs to be divided by the mass to get sigma^2 */
1045 kT_fac = BOLTZ*reference_temperature*gmx::square((ir->nstlist-1)*ir->delta_t);
1048 mass_min = att[0].prop.mass;
1049 for (i = 1; i < natt; i++)
1051 mass_min = std::min(mass_min, att[i].prop.mass);
1056 fprintf(debug, "md1_ljd %9.2e d2_ljd %9.2e md3_ljd %9.2e\n", md1_ljd, d2_ljd, md3_ljd);
1057 fprintf(debug, "md1_ljr %9.2e d2_ljr %9.2e md3_ljr %9.2e\n", md1_ljr, d2_ljr, md3_ljr);
1058 fprintf(debug, "md1_el %9.2e d2_el %9.2e\n", md1_el, d2_el);
1059 fprintf(debug, "sqrt(kT_fac) %f\n", std::sqrt(kT_fac));
1060 fprintf(debug, "mass_min %f\n", mass_min);
1063 /* Search using bisection */
1065 /* The drift will be neglible at 5 times the max sigma */
1066 ib1 = (int)(5*2*std::sqrt(kT_fac/mass_min)/resolution) + 1;
1067 while (ib1 - ib0 > 1)
1071 rl = std::max(ir->rvdw, ir->rcoulomb) + rb;
1073 /* Calculate the average energy drift at the last step
1074 * of the nstlist steps at which the pair-list is used.
1076 drift = ener_drift(att, natt, &mtop->ffparams,
1078 md1_ljd, d2_ljd, md3_ljd,
1079 md1_ljr, d2_ljr, md3_ljr,
1084 /* Correct for the fact that we are using a Ni x Nj particle pair list
1085 * and not a 1 x 1 particle pair list. This reduces the drift.
1087 /* We don't have a formula for 8 (yet), use 4 which is conservative */
1088 nb_clust_frac_pairs_not_in_list_at_cutoff =
1089 surface_frac(std::min(list_setup->cluster_size_i, 4),
1090 particle_distance, rl)*
1091 surface_frac(std::min(list_setup->cluster_size_j, 4),
1092 particle_distance, rl);
1093 drift *= nb_clust_frac_pairs_not_in_list_at_cutoff;
1095 /* Convert the drift to drift per unit time per atom */
1096 drift /= ir->nstlist*ir->delta_t*mtop->natoms;
1100 fprintf(debug, "ib %3d %3d %3d rb %.3f %dx%d fac %.3f drift %.1e\n",
1102 list_setup->cluster_size_i, list_setup->cluster_size_j,
1103 nb_clust_frac_pairs_not_in_list_at_cutoff,
1107 if (std::abs(drift) > ir->verletbuf_tol)
1119 *rlist = std::max(ir->rvdw, ir->rcoulomb) + ib1*resolution;