2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2012,2013,2014,2015,2016,2017,2018, 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/ewald/ewald-utils.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 orders 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.
98 atom_nonbonded_kinetic_prop_t prop; /* non-bonded and kinetic atom prop. */
99 int n; /* #atoms of this type in the system */
100 } verletbuf_atomtype_t;
102 // Struct for derivatives of a non-bonded interaction potential
105 real md1; // -V' at the cutoff
106 real d2; // V'' at the cutoff
107 real md3; // -V''' at the cutoff
110 VerletbufListSetup verletbufGetListSetup(int nbnxnKernelType)
112 /* Note that the current buffer estimation code only handles clusters
113 * of size 1, 2 or 4, so for 4x8 or 8x8 we use the estimate for 4x4.
115 VerletbufListSetup listSetup;
117 listSetup.cluster_size_i = nbnxn_kernel_to_cluster_i_size(nbnxnKernelType);
118 listSetup.cluster_size_j = nbnxn_kernel_to_cluster_j_size(nbnxnKernelType);
120 if (nbnxnKernelType == nbnxnk8x8x8_GPU ||
121 nbnxnKernelType == nbnxnk8x8x8_PlainC)
123 /* The GPU kernels (except for OpenCL) split the j-clusters in two halves */
124 listSetup.cluster_size_j /= 2;
130 VerletbufListSetup verletbufGetSafeListSetup(ListSetupType listType)
132 /* When calling this function we often don't know which kernel type we
133 * are going to use. We choose the kernel type with the smallest possible
134 * i- and j-cluster sizes, so we potentially overestimate, but never
135 * underestimate, the buffer drift.
139 if (listType == ListSetupType::Gpu)
141 nbnxnKernelType = nbnxnk8x8x8_GPU;
143 else if (GMX_SIMD && listType == ListSetupType::CpuSimdWhenSupported)
145 #ifdef GMX_NBNXN_SIMD_2XNN
146 /* We use the smallest cluster size to be on the safe side */
147 nbnxnKernelType = nbnxnk4xN_SIMD_2xNN;
149 nbnxnKernelType = nbnxnk4xN_SIMD_4xN;
154 nbnxnKernelType = nbnxnk4x4_PlainC;
157 return verletbufGetListSetup(nbnxnKernelType);
161 atom_nonbonded_kinetic_prop_equal(const atom_nonbonded_kinetic_prop_t *prop1,
162 const atom_nonbonded_kinetic_prop_t *prop2)
164 return (prop1->mass == prop2->mass &&
165 prop1->type == prop2->type &&
166 prop1->q == prop2->q &&
167 prop1->bConstr == prop2->bConstr &&
168 prop1->con_mass == prop2->con_mass &&
169 prop1->con_len == prop2->con_len);
172 static void add_at(verletbuf_atomtype_t **att_p, int *natt_p,
173 const atom_nonbonded_kinetic_prop_t *prop,
176 verletbuf_atomtype_t *att;
181 /* Ignore massless particles */
189 while (i < natt && !atom_nonbonded_kinetic_prop_equal(prop, &att[i].prop))
201 srenew(*att_p, *natt_p);
202 (*att_p)[i].prop = *prop;
203 (*att_p)[i].n = nmol;
207 /* Returns the mass of atom atomIndex or 1 when setMassesToOne=true */
208 static real getMass(const t_atoms &atoms,
214 return atoms.atom[atomIndex].m;
222 static void get_vsite_masses(const gmx_moltype_t *moltype,
223 const gmx_ffparams_t *ffparams,
233 /* Check for virtual sites, determine mass from constructing atoms */
234 for (ft = 0; ft < F_NRE; ft++)
238 il = &moltype->ilist[ft];
240 for (i = 0; i < il->nr; i += 1+NRAL(ft))
243 real inv_mass, coeff, m_aj;
246 ip = &ffparams->iparams[il->iatoms[i]];
248 a1 = il->iatoms[i+1];
252 /* Only vsiten can have more than four
253 constructing atoms, so NRAL(ft) <= 5 */
256 const int maxj = NRAL(ft);
260 for (j = 1; j < maxj; j++)
262 int aj = il->iatoms[i + 1 + j];
263 cam[j] = getMass(moltype->atoms, aj, setMassesToOne);
266 cam[j] = vsite_m[aj];
270 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.",
272 interaction_function[ft].longname,
281 vsite_m[a1] = (cam[1]*cam[2])/(cam[2]*gmx::square(1-ip->vsite.a) + cam[1]*gmx::square(ip->vsite.a));
285 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));
288 gmx_incons("Invalid vsite type");
291 /* Use the mass of the lightest constructing atom.
292 * This is an approximation.
293 * If the distance of the virtual site to the
294 * constructing atom is less than all distances
295 * between constructing atoms, this is a safe
296 * over-estimate of the displacement of the vsite.
297 * This condition holds for all H mass replacement
298 * vsite constructions, except for SP2/3 groups.
299 * In SP3 groups one H will have a F_VSITE3
300 * construction, so even there the total drift
301 * estimate shouldn't be far off.
303 vsite_m[a1] = cam[1];
304 for (j = 2; j < maxj; j++)
306 vsite_m[a1] = std::min(vsite_m[a1], cam[j]);
319 for (j = 0; j < 3*ffparams->iparams[il->iatoms[i]].vsiten.n; j += 3)
321 int aj = il->iatoms[i + j + 2];
322 coeff = ffparams->iparams[il->iatoms[i+j]].vsiten.a;
323 if (moltype->atoms.atom[aj].ptype == eptVSite)
329 m_aj = moltype->atoms.atom[aj].m;
333 gmx_incons("The mass of a vsiten constructing atom is <= 0");
335 inv_mass += coeff*coeff/m_aj;
337 vsite_m[a1] = 1/inv_mass;
338 /* Correct for loop increment of i */
339 i += j - 1 - NRAL(ft);
343 fprintf(debug, "atom %4d %-20s mass %6.3f\n",
344 a1, interaction_function[ft].longname, vsite_m[a1]);
351 static void get_verlet_buffer_atomtypes(const gmx_mtop_t *mtop,
353 verletbuf_atomtype_t **att_p,
357 verletbuf_atomtype_t *att;
359 int mb, nmol, ft, i, a1, a2, a3, a;
360 const t_atoms *atoms;
363 atom_nonbonded_kinetic_prop_t *prop;
365 int n_nonlin_vsite_mol;
370 if (n_nonlin_vsite != nullptr)
375 for (mb = 0; mb < mtop->nmolblock; mb++)
377 nmol = mtop->molblock[mb].nmol;
379 atoms = &mtop->moltype[mtop->molblock[mb].type].atoms;
381 /* Check for constraints, as they affect the kinetic energy.
382 * For virtual sites we need the masses and geometry of
383 * the constructing atoms to determine their velocity distribution.
385 snew(prop, atoms->nr);
386 snew(vsite_m, atoms->nr);
388 for (ft = F_CONSTR; ft <= F_CONSTRNC; ft++)
390 il = &mtop->moltype[mtop->molblock[mb].type].ilist[ft];
392 for (i = 0; i < il->nr; i += 1+NRAL(ft))
394 ip = &mtop->ffparams.iparams[il->iatoms[i]];
395 a1 = il->iatoms[i+1];
396 a2 = il->iatoms[i+2];
397 real mass1 = getMass(*atoms, a1, setMassesToOne);
398 real mass2 = getMass(*atoms, a2, setMassesToOne);
399 if (mass2 > prop[a1].con_mass)
401 prop[a1].con_mass = mass2;
402 prop[a1].con_len = ip->constr.dA;
404 if (mass1 > prop[a2].con_mass)
406 prop[a2].con_mass = mass1;
407 prop[a2].con_len = ip->constr.dA;
412 il = &mtop->moltype[mtop->molblock[mb].type].ilist[F_SETTLE];
414 for (i = 0; i < il->nr; i += 1+NRAL(F_SETTLE))
416 ip = &mtop->ffparams.iparams[il->iatoms[i]];
417 a1 = il->iatoms[i+1];
418 a2 = il->iatoms[i+2];
419 a3 = il->iatoms[i+3];
420 /* Usually the mass of a1 (usually oxygen) is larger than a2/a3.
421 * If this is not the case, we overestimate the displacement,
422 * which leads to a larger buffer (ok since this is an exotic case).
424 prop[a1].con_mass = getMass(*atoms, a2, setMassesToOne);
425 prop[a1].con_len = ip->settle.doh;
427 prop[a2].con_mass = getMass(*atoms, a1, setMassesToOne);
428 prop[a2].con_len = ip->settle.doh;
430 prop[a3].con_mass = getMass(*atoms, a1, setMassesToOne);
431 prop[a3].con_len = ip->settle.doh;
434 get_vsite_masses(&mtop->moltype[mtop->molblock[mb].type],
438 &n_nonlin_vsite_mol);
439 if (n_nonlin_vsite != nullptr)
441 *n_nonlin_vsite += nmol*n_nonlin_vsite_mol;
444 for (a = 0; a < atoms->nr; a++)
446 if (atoms->atom[a].ptype == eptVSite)
448 prop[a].mass = vsite_m[a];
452 prop[a].mass = getMass(*atoms, a, setMassesToOne);
454 prop[a].type = atoms->atom[a].type;
455 prop[a].q = atoms->atom[a].q;
456 /* We consider an atom constrained, #DOF=2, when it is
457 * connected with constraints to (at least one) atom with
458 * a mass of more than 0.4x its own mass. This is not a critical
459 * parameter, since with roughly equal masses the unconstrained
460 * and constrained displacement will not differ much (and both
461 * overestimate the displacement).
463 prop[a].bConstr = (prop[a].con_mass > 0.4*prop[a].mass);
465 add_at(&att, &natt, &prop[a], nmol);
474 for (a = 0; a < natt; a++)
476 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",
477 a, att[a].prop.mass, att[a].prop.type, att[a].prop.q,
478 att[a].prop.bConstr, att[a].prop.con_mass, att[a].prop.con_len,
487 /* This function computes two components of the estimate of the variance
488 * in the displacement of one atom in a system of two constrained atoms.
489 * Returns in sigma2_2d the variance due to rotation of the constrained
490 * atom around the atom to which it constrained.
491 * Returns in sigma2_3d the variance due to displacement of the COM
492 * of the whole system of the two constrained atoms.
494 * Note that we only take a single constraint (the one to the heaviest atom)
495 * into account. If an atom has multiple constraints, this will result in
496 * an overestimate of the displacement, which gives a larger drift and buffer.
498 void constrained_atom_sigma2(real kT_fac,
499 const atom_nonbonded_kinetic_prop_t *prop,
503 /* Here we decompose the motion of a constrained atom into two
504 * components: rotation around the COM and translation of the COM.
507 /* Determine the variance of the arc length for the two rotational DOFs */
508 real massFraction = prop->con_mass/(prop->mass + prop->con_mass);
509 real sigma2_rot = kT_fac*massFraction/prop->mass;
511 /* The distance from the atom to the COM, i.e. the rotational arm */
512 real comDistance = prop->con_len*massFraction;
514 /* The variance relative to the arm */
515 real sigma2_rel = sigma2_rot/gmx::square(comDistance);
517 /* For sigma2_rel << 1 we don't notice the rotational effect and
518 * we have a normal, Gaussian displacement distribution.
519 * For larger sigma2_rel the displacement is much less, in fact it can
520 * not exceed 2*comDistance. We can calculate MSD/arm^2 as:
521 * integral_x=0-inf distance2(x) x/sigma2_rel exp(-x^2/(2 sigma2_rel)) dx
522 * where x is angular displacement and distance2(x) is the distance^2
523 * between points at angle 0 and x:
524 * distance2(x) = (sin(x) - sin(0))^2 + (cos(x) - cos(0))^2
525 * The limiting value of this MSD is 2, which is also the value for
526 * a uniform rotation distribution that would be reached at long time.
527 * The maximum is 2.5695 at sigma2_rel = 4.5119.
528 * We approximate this integral with a rational polynomial with
529 * coefficients from a Taylor expansion. This approximation is an
530 * overestimate for all values of sigma2_rel. Its maximum value
531 * of 2.6491 is reached at sigma2_rel = sqrt(45/2) = 4.7434.
532 * We keep the approximation constant after that.
533 * We use this approximate MSD as the variance for a Gaussian distribution.
535 * NOTE: For any sensible buffer tolerance this will result in a (large)
536 * overestimate of the buffer size, since the Gaussian has a long tail,
537 * whereas the actual distribution can not reach values larger than 2.
539 /* Coeffients obtained from a Taylor expansion */
540 const real a = 1.0/3.0;
541 const real b = 2.0/45.0;
543 /* Our approximation is constant after sigma2_rel = 1/sqrt(b) */
544 sigma2_rel = std::min(sigma2_rel, 1/std::sqrt(b));
546 /* Compute the approximate sigma^2 for 2D motion due to the rotation */
547 *sigma2_2d = gmx::square(comDistance)*
548 sigma2_rel/(1 + a*sigma2_rel + b*gmx::square(sigma2_rel));
550 /* The constrained atom also moves (in 3D) with the COM of both atoms */
551 *sigma2_3d = kT_fac/(prop->mass + prop->con_mass);
554 static void get_atom_sigma2(real kT_fac,
555 const atom_nonbonded_kinetic_prop_t *prop,
561 /* Complicated constraint calculation in a separate function */
562 constrained_atom_sigma2(kT_fac, prop, sigma2_2d, sigma2_3d);
566 /* Unconstrained atom: trivial */
568 *sigma2_3d = kT_fac/prop->mass;
572 static void approx_2dof(real s2, real x, real *shift, real *scale)
574 /* A particle with 1 DOF constrained has 2 DOFs instead of 3.
575 * This code is also used for particles with multiple constraints,
576 * in which case we overestimate the displacement.
577 * The 2DOF distribution is sqrt(pi/2)*erfc(r/(sqrt(2)*s))/(2*s).
578 * We approximate this with scale*Gaussian(s,r+shift),
579 * by matching the distribution value and derivative at x.
580 * This is a tight overestimate for all r>=0 at any s and x.
584 ex = std::exp(-x*x/(2*s2));
585 er = std::erfc(x/std::sqrt(2*s2));
587 *shift = -x + std::sqrt(2*s2/M_PI)*ex/er;
588 *scale = 0.5*M_PI*std::exp(ex*ex/(M_PI*er*er))*er;
591 // Returns an (over)estimate of the energy drift for a single atom pair,
592 // given the kinetic properties, displacement variances and list buffer.
593 static real energyDriftAtomPair(bool isConstrained_i,
594 bool isConstrained_j,
595 real s2, real s2i_2d, real s2j_2d,
597 const pot_derivatives_t *der)
599 // For relatively small arguments erfc() is so small that if will be 0.0
600 // when stored in a float. We set an argument limit of 8 (Erfc(8)=1e-29),
601 // such that we can divide by erfc and have some space left for arithmetic.
602 const real erfc_arg_max = 8.0;
609 if (rsh*rsh > 2*s2*erfc_arg_max*erfc_arg_max)
611 // Below we calculate c_erfc = 0.5*erfc(rsh/sqrt(2*s2))
612 // When rsh/sqrt(2*s2) increases, this erfc will be the first
613 // result that underflows and becomes 0.0. To avoid this,
614 // we set c_exp=0 and c_erfc=0 for large arguments.
615 // This also avoids NaN in approx_2dof().
616 // In any relevant case this has no effect on the results,
617 // since c_exp < 6e-29, so the displacement is completely
618 // negligible for such atom pairs (and an overestimate).
619 // In nearly all use cases, there will be other atom pairs
620 // that contribute much more to the total, so zeroing
621 // this particular contribution has no effect at all.
627 /* For constraints: adapt r and scaling for the Gaussian */
632 approx_2dof(s2i_2d, r_buffer*s2i_2d/s2, &sh, &sc);
640 approx_2dof(s2j_2d, r_buffer*s2j_2d/s2, &sh, &sc);
645 /* Exact contribution of an atom pair with Gaussian displacement
646 * with sigma s to the energy drift for a potential with
647 * derivative -md and second derivative dd at the cut-off.
648 * The only catch is that for potentials that change sign
649 * near the cut-off there could be an unlucky compensation
650 * of positive and negative energy drift.
651 * Such potentials are extremely rare though.
653 * Note that pot has unit energy*length, as the linear
654 * atom density still needs to be put in.
656 c_exp = std::exp(-rsh*rsh/(2*s2))/std::sqrt(2*M_PI);
657 c_erfc = 0.5*std::erfc(rsh/(std::sqrt(2*s2)));
659 real s = std::sqrt(s2);
663 der->md1/2*((rsh2 + s2)*c_erfc - rsh*s*c_exp);
665 der->d2/6*(s*(rsh2 + 2*s2)*c_exp - rsh*(rsh2 + 3*s2)*c_erfc);
667 der->md3/24*((rsh2*rsh2 + 6*rsh2*s2 + 3*s2*s2)*c_erfc - rsh*s*(rsh2 + 5*s2)*c_exp);
669 return pot1 + pot2 + pot3;
672 static real energyDrift(const verletbuf_atomtype_t *att, int natt,
673 const gmx_ffparams_t *ffp,
675 const pot_derivatives_t *ljDisp,
676 const pot_derivatives_t *ljRep,
677 const pot_derivatives_t *elec,
678 real rlj, real rcoulomb,
679 real rlist, real boxvol)
681 double drift_tot = 0;
685 /* No atom displacements: no drift, avoid division by 0 */
689 // Here add up the contribution of all atom pairs in the system to
690 // (estimated) energy drift by looping over all atom type pairs.
691 for (int i = 0; i < natt; i++)
693 // Get the thermal displacement variance for the i-atom type
694 const atom_nonbonded_kinetic_prop_t *prop_i = &att[i].prop;
696 get_atom_sigma2(kT_fac, prop_i, &s2i_2d, &s2i_3d);
698 for (int j = i; j < natt; j++)
700 // Get the thermal displacement variance for the j-atom type
701 const atom_nonbonded_kinetic_prop_t *prop_j = &att[j].prop;
703 get_atom_sigma2(kT_fac, prop_j, &s2j_2d, &s2j_3d);
705 /* Add up the up to four independent variances */
706 real s2 = s2i_2d + s2i_3d + s2j_2d + s2j_3d;
708 // Set -V', V'' and -V''' at the cut-off for LJ */
709 real c6 = ffp->iparams[prop_i->type*ffp->atnr + prop_j->type].lj.c6;
710 real c12 = ffp->iparams[prop_i->type*ffp->atnr + prop_j->type].lj.c12;
711 pot_derivatives_t lj;
712 lj.md1 = c6*ljDisp->md1 + c12*ljRep->md1;
713 lj.d2 = c6*ljDisp->d2 + c12*ljRep->d2;
714 lj.md3 = c6*ljDisp->md3 + c12*ljRep->md3;
716 real pot_lj = energyDriftAtomPair(prop_i->bConstr, prop_j->bConstr,
721 // Set -V' and V'' at the cut-off for Coulomb
722 pot_derivatives_t elec_qq;
723 elec_qq.md1 = elec->md1*prop_i->q*prop_j->q;
724 elec_qq.d2 = elec->d2 *prop_i->q*prop_j->q;
727 real pot_q = energyDriftAtomPair(prop_i->bConstr, prop_j->bConstr,
732 // Note that attractive and repulsive potentials for individual
733 // pairs can partially cancel.
734 real pot = pot_lj + pot_q;
736 /* Multiply by the number of atom pairs */
739 pot *= (double)att[i].n*(att[i].n - 1)/2;
743 pot *= (double)att[i].n*att[j].n;
745 /* We need the line density to get the energy drift of the system.
746 * The effective average r^2 is close to (rlist+sigma)^2.
748 pot *= 4*M_PI*gmx::square(rlist + std::sqrt(s2))/boxvol;
750 /* Add the unsigned drift to avoid cancellation of errors */
751 drift_tot += std::abs(pot);
758 static real surface_frac(int cluster_size, real particle_distance, real rlist)
762 if (rlist < 0.5*particle_distance)
764 /* We have non overlapping spheres */
768 /* Half the inter-particle distance relative to rlist */
769 d = 0.5*particle_distance/rlist;
771 /* Determine the area of the surface at distance rlist to the closest
772 * particle, relative to surface of a sphere of radius rlist.
773 * The formulas below assume close to cubic cells for the pair search grid,
774 * which the pair search code tries to achieve.
775 * Note that in practice particle distances will not be delta distributed,
776 * but have some spread, often involving shorter distances,
777 * as e.g. O-H bonds in a water molecule. Thus the estimates below will
778 * usually be slightly too high and thus conservative.
780 switch (cluster_size)
783 /* One particle: trivial */
787 /* Two particles: two spheres at fractional distance 2*a */
791 /* We assume a perfect, symmetric tetrahedron geometry.
792 * The surface around a tetrahedron is too complex for a full
793 * analytical solution, so we use a Taylor expansion.
795 area_rel = (1.0 + 1/M_PI*(6*std::acos(1/std::sqrt(3))*d +
796 std::sqrt(3)*d*d*(1.0 +
799 83.0/756.0*d*d*d*d*d*d)));
802 gmx_incons("surface_frac called with unsupported cluster_size");
806 return area_rel/cluster_size;
809 /* Returns the negative of the third derivative of a potential r^-p
810 * with a force-switch function, evaluated at the cut-off rc.
812 static real md3_force_switch(real p, real rswitch, real rc)
814 /* The switched force function is:
815 * p*r^-(p+1) + a*(r - rswitch)^2 + b*(r - rswitch)^3
818 real md3_pot, md3_sw;
820 a = -((p + 4)*rc - (p + 1)*rswitch)/(pow(rc, p+2)*gmx::square(rc-rswitch));
821 b = ((p + 3)*rc - (p + 1)*rswitch)/(pow(rc, p+2)*gmx::power3(rc-rswitch));
823 md3_pot = (p + 2)*(p + 1)*p*pow(rc, p+3);
824 md3_sw = 2*a + 6*b*(rc - rswitch);
826 return md3_pot + md3_sw;
829 /* Returns the maximum reference temperature over all coupled groups */
830 static real maxReferenceTemperature(const t_inputrec &ir)
832 if (EI_MD(ir.eI) && ir.etc == etcNO)
834 /* This case should be handled outside calc_verlet_buffer_size */
835 gmx_incons("calc_verlet_buffer_size called with an NVE ensemble and reference_temperature < 0");
838 real maxTemperature = 0;
839 for (int i = 0; i < ir.opts.ngtc; i++)
841 if (ir.opts.tau_t[i] >= 0)
843 maxTemperature = std::max(maxTemperature, ir.opts.ref_t[i]);
847 return maxTemperature;
850 /* Returns the variance of the atomic displacement over timePeriod.
852 * Note: When not using BD with a non-mass dependendent friction coefficient,
853 * the return value still needs to be divided by the particle mass.
855 static real displacementVariance(const t_inputrec &ir,
863 /* Get the displacement distribution from the random component only.
864 * With accurate integration the systematic (force) displacement
865 * should be negligible (unless nstlist is extremely large, which
866 * you wouldn't do anyhow).
868 kT_fac = 2*BOLTZ*temperature*timePeriod;
871 /* This is directly sigma^2 of the displacement */
872 kT_fac /= ir.bd_fric;
876 /* Per group tau_t is not implemented yet, use the maximum */
877 real tau_t = ir.opts.tau_t[0];
878 for (int i = 1; i < ir.opts.ngtc; i++)
880 tau_t = std::max(tau_t, ir.opts.tau_t[i]);
884 /* This kT_fac needs to be divided by the mass to get sigma^2 */
889 kT_fac = BOLTZ*temperature*gmx::square(timePeriod);
895 /* Returns the largest sigma of the Gaussian displacement over all particle
896 * types. This ignores constraints, so is an overestimate.
898 static real maxSigma(real kT_fac,
900 const verletbuf_atomtype_t *att)
903 real smallestMass = att[0].prop.mass;
904 for (int i = 1; i < natt; i++)
906 smallestMass = std::min(smallestMass, att[i].prop.mass);
909 return 2*std::sqrt(kT_fac/smallestMass);
912 void calc_verlet_buffer_size(const gmx_mtop_t *mtop, real boxvol,
913 const t_inputrec *ir,
916 real reference_temperature,
917 const VerletbufListSetup *list_setup,
924 real particle_distance;
925 real nb_clust_frac_pairs_not_in_list_at_cutoff;
927 verletbuf_atomtype_t *att = nullptr;
934 if (!EI_DYNAMICS(ir->eI))
936 gmx_incons("Can only determine the Verlet buffer size for integrators that perform dynamics");
938 if (ir->verletbuf_tol <= 0)
940 gmx_incons("The Verlet buffer tolerance needs to be larger than zero");
943 if (reference_temperature < 0)
945 /* We use the maximum temperature with multiple T-coupl groups.
946 * We could use a per particle temperature, but since particles
947 * interact, this might underestimate the buffer size.
949 reference_temperature = maxReferenceTemperature(*ir);
952 /* Resolution of the buffer size */
955 env = getenv("GMX_VERLET_BUFFER_RES");
958 sscanf(env, "%lf", &resolution);
961 /* In an atom wise pair-list there would be no pairs in the list
962 * beyond the pair-list cut-off.
963 * However, we use a pair-list of groups vs groups of atoms.
964 * For groups of 4 atoms, the parallelism of SSE instructions, only
965 * 10% of the atoms pairs are not in the list just beyond the cut-off.
966 * As this percentage increases slowly compared to the decrease of the
967 * Gaussian displacement distribution over this range, we can simply
968 * reduce the drift by this fraction.
969 * For larger groups, e.g. of 8 atoms, this fraction will be lower,
970 * so then buffer size will be on the conservative (large) side.
972 * Note that the formulas used here do not take into account
973 * cancellation of errors which could occur by missing both
974 * attractive and repulsive interactions.
976 * The only major assumption is homogeneous particle distribution.
977 * For an inhomogeneous system, such as a liquid-vapor system,
978 * the buffer will be underestimated. The actual energy drift
979 * will be higher by the factor: local/homogeneous particle density.
981 * The results of this estimate have been checked againt simulations.
982 * In most cases the real drift differs by less than a factor 2.
985 /* Worst case assumption: HCP packing of particles gives largest distance */
986 particle_distance = std::cbrt(boxvol*std::sqrt(2)/mtop->natoms);
988 /* TODO: Obtain masses through (future) integrator functionality
989 * to avoid scattering the code with (or forgetting) checks.
991 const bool setMassesToOne = (ir->eI == eiBD && ir->bd_fric > 0);
992 get_verlet_buffer_atomtypes(mtop, setMassesToOne, &att, &natt, n_nonlin_vsite);
993 assert(att != NULL && natt >= 0);
997 fprintf(debug, "particle distance assuming HCP packing: %f nm\n",
999 fprintf(debug, "energy drift atom types: %d\n", natt);
1002 pot_derivatives_t ljDisp = { 0, 0, 0 };
1003 pot_derivatives_t ljRep = { 0, 0, 0 };
1004 real repPow = mtop->ffparams.reppow;
1006 if (ir->vdwtype == evdwCUT)
1008 real sw_range, md3_pswf;
1010 switch (ir->vdw_modifier)
1013 case eintmodPOTSHIFT:
1014 /* -dV/dr of -r^-6 and r^-reppow */
1015 ljDisp.md1 = -6*std::pow(ir->rvdw, -7.0);
1016 ljRep.md1 = repPow*std::pow(ir->rvdw, -(repPow + 1));
1017 /* The contribution of the higher derivatives is negligible */
1019 case eintmodFORCESWITCH:
1020 /* At the cut-off: V=V'=V''=0, so we use only V''' */
1021 ljDisp.md3 = -md3_force_switch(6.0, ir->rvdw_switch, ir->rvdw);
1022 ljRep.md3 = md3_force_switch(repPow, ir->rvdw_switch, ir->rvdw);
1024 case eintmodPOTSWITCH:
1025 /* At the cut-off: V=V'=V''=0.
1026 * V''' is given by the original potential times
1027 * the third derivative of the switch function.
1029 sw_range = ir->rvdw - ir->rvdw_switch;
1030 md3_pswf = 60.0/gmx::power3(sw_range);
1032 ljDisp.md3 = -std::pow(ir->rvdw, -6.0 )*md3_pswf;
1033 ljRep.md3 = std::pow(ir->rvdw, -repPow)*md3_pswf;
1036 gmx_incons("Unimplemented VdW modifier");
1039 else if (EVDW_PME(ir->vdwtype))
1041 real b = calc_ewaldcoeff_lj(ir->rvdw, ir->ewald_rtol_lj);
1047 // -dV/dr of g(br)*r^-6 [where g(x) = exp(-x^2)(1+x^2+x^4/2),
1048 // see LJ-PME equations in manual] and r^-reppow
1049 ljDisp.md1 = -std::exp(-br2)*(br6 + 3.0*br4 + 6.0*br2 + 6.0)*std::pow(r, -7.0);
1050 ljRep.md1 = repPow*pow(r, -(repPow + 1));
1051 // The contribution of the higher derivatives is negligible
1055 gmx_fatal(FARGS, "Energy drift calculation is only implemented for plain cut-off Lennard-Jones interactions");
1058 elfac = ONE_4PI_EPS0/ir->epsilon_r;
1060 // Determine the 1st and 2nd derivative for the electostatics
1061 pot_derivatives_t elec = { 0, 0, 0 };
1063 if (ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype))
1067 if (ir->coulombtype == eelCUT)
1074 eps_rf = ir->epsilon_rf/ir->epsilon_r;
1077 k_rf = (eps_rf - ir->epsilon_r)/( gmx::power3(ir->rcoulomb) * (2*eps_rf + ir->epsilon_r) );
1081 /* epsilon_rf = infinity */
1082 k_rf = 0.5/gmx::power3(ir->rcoulomb);
1088 elec.md1 = elfac*(1.0/gmx::square(ir->rcoulomb) - 2*k_rf*ir->rcoulomb);
1090 elec.d2 = elfac*(2.0/gmx::power3(ir->rcoulomb) + 2*k_rf);
1092 else if (EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD)
1096 b = calc_ewaldcoeff_q(ir->rcoulomb, ir->ewald_rtol);
1099 elec.md1 = elfac*(b*std::exp(-br*br)*M_2_SQRTPI/rc + std::erfc(br)/(rc*rc));
1100 elec.d2 = elfac/(rc*rc)*(2*b*(1 + br*br)*std::exp(-br*br)*M_2_SQRTPI + 2*std::erfc(br)/rc);
1104 gmx_fatal(FARGS, "Energy drift calculation is only implemented for Reaction-Field and Ewald electrostatics");
1107 /* Determine the variance of the atomic displacement
1108 * over list_lifetime steps: kT_fac
1109 * For inertial dynamics (not Brownian dynamics) the mass factor
1110 * is not included in kT_fac, it is added later.
1112 const real kT_fac = displacementVariance(*ir, reference_temperature,
1113 list_lifetime*ir->delta_t);
1117 fprintf(debug, "Derivatives of non-bonded potentials at the cut-off:\n");
1118 fprintf(debug, "LJ disp. -V' %9.2e V'' %9.2e -V''' %9.2e\n", ljDisp.md1, ljDisp.d2, ljDisp.md3);
1119 fprintf(debug, "LJ rep. -V' %9.2e V'' %9.2e -V''' %9.2e\n", ljRep.md1, ljRep.d2, ljRep.md3);
1120 fprintf(debug, "Electro. -V' %9.2e V'' %9.2e\n", elec.md1, elec.d2);
1121 fprintf(debug, "sqrt(kT_fac) %f\n", std::sqrt(kT_fac));
1124 /* Search using bisection */
1126 /* The drift will be neglible at 5 times the max sigma */
1127 ib1 = (int)(5*maxSigma(kT_fac, natt, att)/resolution) + 1;
1128 while (ib1 - ib0 > 1)
1132 rl = std::max(ir->rvdw, ir->rcoulomb) + rb;
1134 /* Calculate the average energy drift at the last step
1135 * of the nstlist steps at which the pair-list is used.
1137 drift = energyDrift(att, natt, &mtop->ffparams,
1139 &ljDisp, &ljRep, &elec,
1140 ir->rvdw, ir->rcoulomb,
1143 /* Correct for the fact that we are using a Ni x Nj particle pair list
1144 * and not a 1 x 1 particle pair list. This reduces the drift.
1146 /* We don't have a formula for 8 (yet), use 4 which is conservative */
1147 nb_clust_frac_pairs_not_in_list_at_cutoff =
1148 surface_frac(std::min(list_setup->cluster_size_i, 4),
1149 particle_distance, rl)*
1150 surface_frac(std::min(list_setup->cluster_size_j, 4),
1151 particle_distance, rl);
1152 drift *= nb_clust_frac_pairs_not_in_list_at_cutoff;
1154 /* Convert the drift to drift per unit time per atom */
1155 drift /= nstlist*ir->delta_t*mtop->natoms;
1159 fprintf(debug, "ib %3d %3d %3d rb %.3f %dx%d fac %.3f drift %.1e\n",
1161 list_setup->cluster_size_i, list_setup->cluster_size_j,
1162 nb_clust_frac_pairs_not_in_list_at_cutoff,
1166 if (std::abs(drift) > ir->verletbuf_tol)
1178 *rlist = std::max(ir->rvdw, ir->rcoulomb) + ib1*resolution;
1181 /* Returns the pairlist buffer size for use as a minimum buffer size
1183 * Note that this is a rather crude estimate. It is ok for a buffer
1184 * set for good energy conservation or RF electrostatics. But it is
1185 * too small with PME and the buffer set with the default tolerance.
1187 static real minCellSizeFromPairlistBuffer(const t_inputrec &ir)
1189 return ir.rlist - std::max(ir.rvdw, ir.rcoulomb);
1192 real minCellSizeForAtomDisplacement(const gmx_mtop_t &mtop,
1193 const t_inputrec &ir,
1194 real chanceRequested)
1196 if (!EI_DYNAMICS(ir.eI) || (EI_MD(ir.eI) && ir.etc == etcNO))
1198 return minCellSizeFromPairlistBuffer(ir);
1201 /* We use the maximum temperature with multiple T-coupl groups.
1202 * We could use a per particle temperature, but since particles
1203 * interact, this might underestimate the displacements.
1205 const real temperature = maxReferenceTemperature(ir);
1207 const bool setMassesToOne = (ir.eI == eiBD && ir.bd_fric > 0);
1209 verletbuf_atomtype_t *att = nullptr;
1211 get_verlet_buffer_atomtypes(&mtop, setMassesToOne, &att, &natt, nullptr);
1213 const real kT_fac = displacementVariance(ir, temperature,
1214 ir.nstlist*ir.delta_t);
1216 /* Resolution of the cell size */
1217 real resolution = 0.001;
1219 /* Search using bisection, avoid 0 and start at 1 */
1221 /* The chance will be neglible at 10 times the max sigma */
1222 int ib1 = (int)(10*maxSigma(kT_fac, natt, att)/resolution) + 1;
1224 while (ib1 - ib0 > 1)
1226 int ib = (ib0 + ib1)/2;
1227 cellSize = ib*resolution;
1229 /* We assumes atom are distributed uniformly over the cell width.
1230 * Once an atom has moved by more than the cellSize (as passed
1231 * as the buffer argument to energyDriftAtomPair() below),
1232 * the chance of crossing the boundary of the neighbor cell
1233 * thus increases as 1/cellSize with the additional displacement
1234 * on to of cellSize. We thus create a linear interaction with
1235 * derivative = -1/cellSize. Using this in the energyDriftAtomPair
1236 * function will return the chance of crossing the next boundary.
1238 const pot_derivatives_t boundaryInteraction = { 1/cellSize, 0, 0 };
1241 for (int i = 0; i < natt; i++)
1243 const atom_nonbonded_kinetic_prop_t &propAtom = att[i].prop;
1246 get_atom_sigma2(kT_fac, &propAtom, &s2_2d, &s2_3d);
1248 real chancePerAtom = energyDriftAtomPair(propAtom.bConstr, false,
1249 s2_2d + s2_3d, s2_2d, 0,
1251 &boundaryInteraction);
1253 if (propAtom.bConstr)
1255 /* energyDriftAtomPair() uses an unlimited Gaussian displacement
1256 * distribution for constrained atoms, whereas they can
1257 * actually not move more than the COM of the two constrained
1258 * atoms plus twice the distance from the COM.
1259 * Use this maximum, limited displacement when this results in
1260 * a smaller chance (note that this is still an overestimate).
1262 real massFraction = propAtom.con_mass/(propAtom.mass + propAtom.con_mass);
1263 real comDistance = propAtom.con_len*massFraction;
1265 real chanceWithMaxDistance =
1266 energyDriftAtomPair(false, false,
1268 cellSize - 2*comDistance,
1269 &boundaryInteraction);
1270 chancePerAtom = std::min(chancePerAtom, chanceWithMaxDistance);
1273 /* Take into account the line density of the boundary */
1274 chancePerAtom /= cellSize;
1276 chance += att[i].n*chancePerAtom;
1279 /* Note: chance is for every nstlist steps */
1280 if (chance > chanceRequested*ir.nstlist)