1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * For more info, check our website at http://www.gromacs.org
34 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
42 #include <sys/types.h>
47 #include "gmx_fatal.h"
51 #include "calc_verletbuf.h"
52 #include "../mdlib/nbnxn_consts.h"
55 /* The include below sets the SIMD instruction type (precision+width)
56 * for all nbnxn SIMD search and non-bonded kernel code.
58 #ifdef GMX_NBNXN_HALF_WIDTH_SIMD
59 #define GMX_USE_HALF_WIDTH_SIMD_HERE
61 #include "gmx_simd_macros.h"
65 /* The code in this file estimates a pairlist buffer length
66 * given a target energy drift per atom per picosecond.
67 * This is done by estimating the drift given a buffer length.
68 * Ideally we would like to have a tight overestimate of the drift,
69 * but that can be difficult to achieve.
71 * Significant approximations used:
73 * Uniform particle density. UNDERESTIMATES the drift by rho_global/rho_local.
75 * Interactions don't affect particle motion. OVERESTIMATES the drift on longer
76 * time scales. This approximation probably introduces the largest errors.
78 * Only take one constraint per particle into account: OVERESTIMATES the drift.
80 * For rotating constraints assume the same functional shape for time scales
81 * where the constraints rotate significantly as the exact expression for
82 * short time scales. OVERESTIMATES the drift on long time scales.
84 * For non-linear virtual sites use the mass of the lightest constructing atom
85 * to determine the displacement. OVER/UNDERESTIMATES the drift, depending on
86 * the geometry and masses of constructing atoms.
88 * Note that the formulas for normal atoms and linear virtual sites are exact,
89 * apart from the first two approximations.
91 * Note that apart from the effect of the above approximations, the actual
92 * drift of the total energy of a system can be order of magnitude smaller
93 * due to cancellation of positive and negative drift for different pairs.
97 /* Struct for unique atom type for calculating the energy drift.
98 * The atom displacement depends on mass and constraints.
99 * The energy jump for given distance depend on LJ type and q.
103 real mass; /* mass */
104 int type; /* type (used for LJ parameters) */
106 gmx_bool bConstr; /* constrained, if TRUE, use #DOF=2 iso 3 */
107 real con_mass; /* mass of heaviest atom connected by constraints */
108 real con_len; /* constraint length to the heaviest atom */
109 } atom_nonbonded_kinetic_prop_t;
111 /* Struct for unique atom type for calculating the energy drift.
112 * The atom displacement depends on mass and constraints.
113 * The energy jump for given distance depend on LJ type and q.
117 atom_nonbonded_kinetic_prop_t prop; /* non-bonded and kinetic atom prop. */
118 int n; /* #atoms of this type in the system */
119 } verletbuf_atomtype_t;
121 void verletbuf_get_list_setup(gmx_bool bGPU,
122 verletbuf_list_setup_t *list_setup)
124 list_setup->cluster_size_i = NBNXN_CPU_CLUSTER_I_SIZE;
128 list_setup->cluster_size_j = NBNXN_GPU_CLUSTER_SIZE;
132 #ifndef GMX_NBNXN_SIMD
133 list_setup->cluster_size_j = NBNXN_CPU_CLUSTER_I_SIZE;
135 list_setup->cluster_size_j = GMX_SIMD_WIDTH_HERE;
136 #ifdef GMX_NBNXN_SIMD_2XNN
137 /* We assume the smallest cluster size to be on the safe side */
138 list_setup->cluster_size_j /= 2;
145 atom_nonbonded_kinetic_prop_equal(const atom_nonbonded_kinetic_prop_t *prop1,
146 const atom_nonbonded_kinetic_prop_t *prop2)
148 return (prop1->mass == prop2->mass &&
149 prop1->type == prop2->type &&
150 prop1->q == prop2->q &&
151 prop1->bConstr == prop2->bConstr &&
152 prop1->con_mass == prop2->con_mass &&
153 prop1->con_len == prop2->con_len);
156 static void add_at(verletbuf_atomtype_t **att_p, int *natt_p,
157 const atom_nonbonded_kinetic_prop_t *prop,
160 verletbuf_atomtype_t *att;
165 /* Ignore massless particles */
173 while (i < natt && !atom_nonbonded_kinetic_prop_equal(prop, &att[i].prop))
185 srenew(*att_p, *natt_p);
186 (*att_p)[i].prop = *prop;
187 (*att_p)[i].n = nmol;
191 static void get_vsite_masses(const gmx_moltype_t *moltype,
192 const gmx_ffparams_t *ffparams,
201 /* Check for virtual sites, determine mass from constructing atoms */
202 for (ft = 0; ft < F_NRE; ft++)
206 il = &moltype->ilist[ft];
208 for (i = 0; i < il->nr; i += 1+NRAL(ft))
211 real cam[5], inv_mass, m_aj;
212 int a1, j, aj, coeff;
214 ip = &ffparams->iparams[il->iatoms[i]];
216 a1 = il->iatoms[i+1];
220 for (j = 1; j < NRAL(ft); j++)
222 cam[j] = moltype->atoms.atom[il->iatoms[i+1+j]].m;
225 cam[j] = vsite_m[il->iatoms[i+1+j]];
229 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.",
231 interaction_function[ft].longname,
232 il->iatoms[i+1+j]+1);
241 vsite_m[a1] = (cam[1]*cam[2])/(cam[2]*sqr(1-ip->vsite.a) + cam[1]*sqr(ip->vsite.a));
245 vsite_m[a1] = (cam[1]*cam[2]*cam[3])/(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));
250 for (j = 0; j < 3*ip->vsiten.n; j += 3)
252 aj = il->iatoms[i+j+2];
253 coeff = ip[il->iatoms[i+j]].vsiten.a;
254 if (moltype->atoms.atom[aj].ptype == eptVSite)
260 m_aj = moltype->atoms.atom[aj].m;
264 gmx_incons("The mass of a vsiten constructing atom is <= 0");
266 inv_mass += coeff*coeff/m_aj;
268 vsite_m[a1] = 1/inv_mass;
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.
284 vsite_m[a1] = cam[1];
285 for (j = 2; j < NRAL(ft); j++)
287 vsite_m[a1] = min(vsite_m[a1], cam[j]);
294 fprintf(debug, "atom %4d %-20s mass %6.3f\n",
295 a1, interaction_function[ft].longname, vsite_m[a1]);
302 static void get_verlet_buffer_atomtypes(const gmx_mtop_t *mtop,
303 verletbuf_atomtype_t **att_p,
307 verletbuf_atomtype_t *att;
309 int mb, nmol, ft, i, a1, a2, a3, a;
310 const t_atoms *atoms;
313 atom_nonbonded_kinetic_prop_t *prop;
315 int n_nonlin_vsite_mol;
320 if (n_nonlin_vsite != NULL)
325 for (mb = 0; mb < mtop->nmolblock; mb++)
327 nmol = mtop->molblock[mb].nmol;
329 atoms = &mtop->moltype[mtop->molblock[mb].type].atoms;
331 /* Check for constraints, as they affect the kinetic energy.
332 * For virtual sites we need the masses and geometry of
333 * the constructing atoms to determine their velocity distribution.
335 snew(prop, atoms->nr);
336 snew(vsite_m, atoms->nr);
338 for (ft = F_CONSTR; ft <= F_CONSTRNC; ft++)
340 il = &mtop->moltype[mtop->molblock[mb].type].ilist[ft];
342 for (i = 0; i < il->nr; i += 1+NRAL(ft))
344 ip = &mtop->ffparams.iparams[il->iatoms[i]];
345 a1 = il->iatoms[i+1];
346 a2 = il->iatoms[i+2];
347 if (atoms->atom[a2].m > prop[a1].con_mass)
349 prop[a1].con_mass = atoms->atom[a2].m;
350 prop[a1].con_len = ip->constr.dA;
352 if (atoms->atom[a1].m > prop[a2].con_mass)
354 prop[a2].con_mass = atoms->atom[a1].m;
355 prop[a2].con_len = ip->constr.dA;
360 il = &mtop->moltype[mtop->molblock[mb].type].ilist[F_SETTLE];
362 for (i = 0; i < il->nr; i += 1+NRAL(F_SETTLE))
364 ip = &mtop->ffparams.iparams[il->iatoms[i]];
365 a1 = il->iatoms[i+1];
366 a2 = il->iatoms[i+2];
367 a3 = il->iatoms[i+3];
368 /* Usually the mass of a1 (usually oxygen) is larger than a2/a3.
369 * If this is not the case, we overestimate the displacement,
370 * which leads to a larger buffer (ok since this is an exotic case).
372 prop[a1].con_mass = atoms->atom[a2].m;
373 prop[a1].con_len = ip->settle.doh;
375 prop[a2].con_mass = atoms->atom[a1].m;
376 prop[a2].con_len = ip->settle.doh;
378 prop[a3].con_mass = atoms->atom[a1].m;
379 prop[a3].con_len = ip->settle.doh;
382 get_vsite_masses(&mtop->moltype[mtop->molblock[mb].type],
385 &n_nonlin_vsite_mol);
386 if (n_nonlin_vsite != NULL)
388 *n_nonlin_vsite += nmol*n_nonlin_vsite_mol;
391 for (a = 0; a < atoms->nr; a++)
393 if (atoms->atom[a].ptype == eptVSite)
395 prop[a].mass = vsite_m[a];
399 prop[a].mass = atoms->atom[a].m;
401 prop[a].type = atoms->atom[a].type;
402 prop[a].q = atoms->atom[a].q;
403 /* We consider an atom constrained, #DOF=2, when it is
404 * connected with constraints to (at least one) atom with
405 * a mass of more than 0.4x its own mass. This is not a critical
406 * parameter, since with roughly equal masses the unconstrained
407 * and constrained displacement will not differ much (and both
408 * overestimate the displacement).
410 prop[a].bConstr = (prop[a].con_mass > 0.4*prop[a].mass);
412 add_at(&att, &natt, &prop[a], nmol);
421 for (a = 0; a < natt; a++)
423 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",
424 a, att[a].prop.mass, att[a].prop.type, att[a].prop.q,
425 att[a].prop.bConstr, att[a].prop.con_mass, att[a].prop.con_len,
434 /* This function computes two components of the estimate of the variance
435 * in the displacement of one atom in a system of two constrained atoms.
436 * Returns in sigma2_2d the variance due to rotation of the constrained
437 * atom around the atom to which it constrained.
438 * Returns in sigma2_3d the variance due to displacement of the COM
439 * of the whole system of the two constrained atoms.
441 * Note that we only take a single constraint (the one to the heaviest atom)
442 * into account. If an atom has multiple constraints, this will result in
443 * an overestimate of the displacement, which gives a larger drift and buffer.
445 static void constrained_atom_sigma2(real kT_fac,
446 const atom_nonbonded_kinetic_prop_t *prop,
455 /* Here we decompose the motion of a constrained atom into two
456 * components: rotation around the COM and translation of the COM.
459 /* Determine the variance for the displacement of the rotational mode */
460 sigma2_rot = kT_fac/(prop->mass*(prop->mass + prop->con_mass)/prop->con_mass);
462 /* The distance from the atom to the COM, i.e. the rotational arm */
463 com_dist = prop->con_len*prop->con_mass/(prop->mass + prop->con_mass);
465 /* The variance relative to the arm */
466 sigma2_rel = sigma2_rot/(com_dist*com_dist);
467 /* At 6 the scaling formula has slope 0,
468 * so we keep sigma2_2d constant after that.
472 /* A constrained atom rotates around the atom it is constrained to.
473 * This results in a smaller linear displacement than for a free atom.
474 * For a perfectly circular displacement, this lowers the displacement
475 * by: 1/arcsin(arc_length)
476 * and arcsin(x) = 1 + x^2/6 + ...
477 * For sigma2_rel<<1 the displacement distribution is erfc
478 * (exact formula is provided below). For larger sigma, it is clear
479 * that the displacement can't be larger than 2*com_dist.
480 * It turns out that the distribution becomes nearly uniform.
481 * For intermediate sigma2_rel, scaling down sigma with the third
482 * order expansion of arcsin with argument sigma_rel turns out
483 * to give a very good approximation of the distribution and variance.
484 * Even for larger values, the variance is only slightly overestimated.
485 * Note that the most relevant displacements are in the long tail.
486 * This rotation approximation always overestimates the tail (which
487 * runs to infinity, whereas it should be <= 2*com_dist).
488 * Thus we always overestimate the drift and the buffer size.
490 scale = 1/(1 + sigma2_rel/6);
491 *sigma2_2d = sigma2_rot*scale*scale;
495 /* sigma_2d is set to the maximum given by the scaling above.
496 * For large sigma2 the real displacement distribution is close
497 * to uniform over -2*con_len to 2*com_dist.
498 * Our erfc with sigma_2d=sqrt(1.5)*com_dist (which means the sigma
499 * of the erfc output distribution is con_dist) overestimates
500 * the variance and additionally has a long tail. This means
501 * we have a (safe) overestimation of the drift.
503 *sigma2_2d = 1.5*com_dist*com_dist;
506 /* The constrained atom also moves (in 3D) with the COM of both atoms */
507 *sigma2_3d = kT_fac/(prop->mass + prop->con_mass);
510 static void get_atom_sigma2(real kT_fac,
511 const atom_nonbonded_kinetic_prop_t *prop,
517 /* Complicated constraint calculation in a separate function */
518 constrained_atom_sigma2(kT_fac, prop, sigma2_2d, sigma2_3d);
522 /* Unconstrained atom: trivial */
524 *sigma2_3d = kT_fac/prop->mass;
528 static void approx_2dof(real s2, real x, real *shift, real *scale)
530 /* A particle with 1 DOF constrained has 2 DOFs instead of 3.
531 * This code is also used for particles with multiple constraints,
532 * in which case we overestimate the displacement.
533 * The 2DOF distribution is sqrt(pi/2)*erfc(r/(sqrt(2)*s))/(2*s).
534 * We approximate this with scale*Gaussian(s,r+shift),
535 * by matching the distribution value and derivative at x.
536 * This is a tight overestimate for all r>=0 at any s and x.
540 ex = exp(-x*x/(2*s2));
541 er = gmx_erfc(x/sqrt(2*s2));
543 *shift = -x + sqrt(2*s2/M_PI)*ex/er;
544 *scale = 0.5*M_PI*exp(ex*ex/(M_PI*er*er))*er;
547 static real ener_drift(const verletbuf_atomtype_t *att, int natt,
548 const gmx_ffparams_t *ffp,
550 real md_ljd, real md_ljr, real md_el, real dd_el,
552 real rlist, real boxvol)
554 double drift_tot, pot1, pot2, pot;
556 real s2i_2d, s2i_3d, s2j_2d, s2j_3d, s2, s;
560 double c_exp, c_erfc;
564 /* Loop over the different atom type pairs */
565 for (i = 0; i < natt; i++)
567 get_atom_sigma2(kT_fac, &att[i].prop, &s2i_2d, &s2i_3d);
568 ti = att[i].prop.type;
570 for (j = i; j < natt; j++)
572 get_atom_sigma2(kT_fac, &att[j].prop, &s2j_2d, &s2j_3d);
573 tj = att[j].prop.type;
575 /* Add up the up to four independent variances */
576 s2 = s2i_2d + s2i_3d + s2j_2d + s2j_3d;
578 /* Note that attractive and repulsive potentials for individual
579 * pairs will partially cancel.
581 /* -dV/dr at the cut-off for LJ + Coulomb */
583 md_ljd*ffp->iparams[ti*ffp->atnr+tj].lj.c6 +
584 md_ljr*ffp->iparams[ti*ffp->atnr+tj].lj.c12 +
585 md_el*att[i].prop.q*att[j].prop.q;
587 /* d2V/dr2 at the cut-off for Coulomb, we neglect LJ */
588 dd = dd_el*att[i].prop.q*att[j].prop.q;
592 /* For constraints: adapt r and scaling for the Gaussian */
593 if (att[i].prop.bConstr)
597 approx_2dof(s2i_2d, r_buffer*s2i_2d/s2, &sh, &sc);
601 if (att[j].prop.bConstr)
605 approx_2dof(s2j_2d, r_buffer*s2j_2d/s2, &sh, &sc);
610 /* Exact contribution of an atom pair with Gaussian displacement
611 * with sigma s to the energy drift for a potential with
612 * derivative -md and second derivative dd at the cut-off.
613 * The only catch is that for potentials that change sign
614 * near the cut-off there could be an unlucky compensation
615 * of positive and negative energy drift.
616 * Such potentials are extremely rare though.
618 * Note that pot has unit energy*length, as the linear
619 * atom density still needs to be put in.
621 c_exp = exp(-rsh*rsh/(2*s2))/sqrt(2*M_PI);
622 c_erfc = 0.5*gmx_erfc(rsh/(sqrt(2*s2)));
626 md/2*((rsh*rsh + s2)*c_erfc - rsh*s*c_exp);
628 dd/6*(s*(rsh*rsh + 2*s2)*c_exp - rsh*(rsh*rsh + 3*s2)*c_erfc);
633 fprintf(debug, "n %d %d d s %.3f %.3f %.3f %.3f con %d md %8.1e dd %8.1e pot1 %8.1e pot2 %8.1e pot %8.1e\n",
635 sqrt(s2i_2d), sqrt(s2i_3d),
636 sqrt(s2j_2d), sqrt(s2j_3d),
637 att[i].prop.bConstr+att[j].prop.bConstr,
638 md, dd, pot1, pot2, pot);
641 /* Multiply by the number of atom pairs */
644 pot *= (double)att[i].n*(att[i].n - 1)/2;
648 pot *= (double)att[i].n*att[j].n;
650 /* We need the line density to get the energy drift of the system.
651 * The effective average r^2 is close to (rlist+sigma)^2.
653 pot *= 4*M_PI*sqr(rlist + s)/boxvol;
655 /* Add the unsigned drift to avoid cancellation of errors */
656 drift_tot += fabs(pot);
663 static real surface_frac(int cluster_size, real particle_distance, real rlist)
667 if (rlist < 0.5*particle_distance)
669 /* We have non overlapping spheres */
673 /* Half the inter-particle distance relative to rlist */
674 d = 0.5*particle_distance/rlist;
676 /* Determine the area of the surface at distance rlist to the closest
677 * particle, relative to surface of a sphere of radius rlist.
678 * The formulas below assume close to cubic cells for the pair search grid,
679 * which the pair search code tries to achieve.
680 * Note that in practice particle distances will not be delta distributed,
681 * but have some spread, often involving shorter distances,
682 * as e.g. O-H bonds in a water molecule. Thus the estimates below will
683 * usually be slightly too high and thus conservative.
685 switch (cluster_size)
688 /* One particle: trivial */
692 /* Two particles: two spheres at fractional distance 2*a */
696 /* We assume a perfect, symmetric tetrahedron geometry.
697 * The surface around a tetrahedron is too complex for a full
698 * analytical solution, so we use a Taylor expansion.
700 area_rel = (1.0 + 1/M_PI*(6*acos(1/sqrt(3))*d +
704 83.0/756.0*d*d*d*d*d*d)));
707 gmx_incons("surface_frac called with unsupported cluster_size");
711 return area_rel/cluster_size;
714 void calc_verlet_buffer_size(const gmx_mtop_t *mtop, real boxvol,
715 const t_inputrec *ir, real drift_target,
716 const verletbuf_list_setup_t *list_setup,
723 real particle_distance;
724 real nb_clust_frac_pairs_not_in_list_at_cutoff;
726 verletbuf_atomtype_t *att = NULL;
729 real md_ljd, md_ljr, md_el, dd_el;
731 real kT_fac, mass_min;
736 /* Resolution of the buffer size */
739 env = getenv("GMX_VERLET_BUFFER_RES");
742 sscanf(env, "%lf", &resolution);
745 /* In an atom wise pair-list there would be no pairs in the list
746 * beyond the pair-list cut-off.
747 * However, we use a pair-list of groups vs groups of atoms.
748 * For groups of 4 atoms, the parallelism of SSE instructions, only
749 * 10% of the atoms pairs are not in the list just beyond the cut-off.
750 * As this percentage increases slowly compared to the decrease of the
751 * Gaussian displacement distribution over this range, we can simply
752 * reduce the drift by this fraction.
753 * For larger groups, e.g. of 8 atoms, this fraction will be lower,
754 * so then buffer size will be on the conservative (large) side.
756 * Note that the formulas used here do not take into account
757 * cancellation of errors which could occur by missing both
758 * attractive and repulsive interactions.
760 * The only major assumption is homogeneous particle distribution.
761 * For an inhomogeneous system, such as a liquid-vapor system,
762 * the buffer will be underestimated. The actual energy drift
763 * will be higher by the factor: local/homogeneous particle density.
765 * The results of this estimate have been checked againt simulations.
766 * In most cases the real drift differs by less than a factor 2.
769 /* Worst case assumption: HCP packing of particles gives largest distance */
770 particle_distance = pow(boxvol*sqrt(2)/mtop->natoms, 1.0/3.0);
772 get_verlet_buffer_atomtypes(mtop, &att, &natt, n_nonlin_vsite);
773 assert(att != NULL && natt >= 0);
777 fprintf(debug, "particle distance assuming HCP packing: %f nm\n",
779 fprintf(debug, "energy drift atom types: %d\n", natt);
782 reppow = mtop->ffparams.reppow;
785 if (ir->vdwtype == evdwCUT)
787 /* -dV/dr of -r^-6 and r^-repporw */
788 md_ljd = -6*pow(ir->rvdw, -7.0);
789 md_ljr = reppow*pow(ir->rvdw, -(reppow+1));
790 /* The contribution of the second derivative is negligible */
794 gmx_fatal(FARGS, "Energy drift calculation is only implemented for plain cut-off Lennard-Jones interactions");
797 elfac = ONE_4PI_EPS0/ir->epsilon_r;
799 /* Determine md=-dV/dr and dd=d^2V/dr^2 */
802 if (ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype))
806 if (ir->coulombtype == eelCUT)
813 eps_rf = ir->epsilon_rf/ir->epsilon_r;
816 k_rf = pow(ir->rcoulomb, -3.0)*(eps_rf - ir->epsilon_r)/(2*eps_rf + ir->epsilon_r);
820 /* epsilon_rf = infinity */
821 k_rf = 0.5*pow(ir->rcoulomb, -3.0);
827 md_el = elfac*(pow(ir->rcoulomb, -2.0) - 2*k_rf*ir->rcoulomb);
829 dd_el = elfac*(2*pow(ir->rcoulomb, -3.0) + 2*k_rf);
831 else if (EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD)
835 b = calc_ewaldcoeff(ir->rcoulomb, ir->ewald_rtol);
838 md_el = elfac*(b*exp(-br*br)*M_2_SQRTPI/rc + gmx_erfc(br)/(rc*rc));
839 dd_el = elfac/(rc*rc)*(2*b*(1 + br*br)*exp(-br*br)*M_2_SQRTPI + 2*gmx_erfc(br)/rc);
843 gmx_fatal(FARGS, "Energy drift calculation is only implemented for Reaction-Field and Ewald electrostatics");
846 /* Determine the variance of the atomic displacement
847 * over nstlist-1 steps: kT_fac
848 * For inertial dynamics (not Brownian dynamics) the mass factor
849 * is not included in kT_fac, it is added later.
853 /* Get the displacement distribution from the random component only.
854 * With accurate integration the systematic (force) displacement
855 * should be negligible (unless nstlist is extremely large, which
856 * you wouldn't do anyhow).
858 kT_fac = 2*BOLTZ*ir->opts.ref_t[0]*(ir->nstlist-1)*ir->delta_t;
861 /* This is directly sigma^2 of the displacement */
862 kT_fac /= ir->bd_fric;
864 /* Set the masses to 1 as kT_fac is the full sigma^2,
865 * but we divide by m in ener_drift().
867 for (i = 0; i < natt; i++)
869 att[i].prop.mass = 1;
876 /* Per group tau_t is not implemented yet, use the maximum */
877 tau_t = ir->opts.tau_t[0];
878 for (i = 1; i < ir->opts.ngtc; i++)
880 tau_t = 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*ir->opts.ref_t[0]*sqr((ir->nstlist-1)*ir->delta_t);
892 mass_min = att[0].prop.mass;
893 for (i = 1; i < natt; i++)
895 mass_min = min(mass_min, att[i].prop.mass);
900 fprintf(debug, "md_ljd %e md_ljr %e\n", md_ljd, md_ljr);
901 fprintf(debug, "md_el %e dd_el %e\n", md_el, dd_el);
902 fprintf(debug, "sqrt(kT_fac) %f\n", sqrt(kT_fac));
903 fprintf(debug, "mass_min %f\n", mass_min);
906 /* Search using bisection */
908 /* The drift will be neglible at 5 times the max sigma */
909 ib1 = (int)(5*2*sqrt(kT_fac/mass_min)/resolution) + 1;
910 while (ib1 - ib0 > 1)
914 rl = max(ir->rvdw, ir->rcoulomb) + rb;
916 /* Calculate the average energy drift at the last step
917 * of the nstlist steps at which the pair-list is used.
919 drift = ener_drift(att, natt, &mtop->ffparams,
921 md_ljd, md_ljr, md_el, dd_el, rb,
924 /* Correct for the fact that we are using a Ni x Nj particle pair list
925 * and not a 1 x 1 particle pair list. This reduces the drift.
927 /* We don't have a formula for 8 (yet), use 4 which is conservative */
928 nb_clust_frac_pairs_not_in_list_at_cutoff =
929 surface_frac(min(list_setup->cluster_size_i, 4),
930 particle_distance, rl)*
931 surface_frac(min(list_setup->cluster_size_j, 4),
932 particle_distance, rl);
933 drift *= nb_clust_frac_pairs_not_in_list_at_cutoff;
935 /* Convert the drift to drift per unit time per atom */
936 drift /= ir->nstlist*ir->delta_t*mtop->natoms;
940 fprintf(debug, "ib %3d %3d %3d rb %.3f %dx%d fac %.3f drift %f\n",
942 list_setup->cluster_size_i, list_setup->cluster_size_j,
943 nb_clust_frac_pairs_not_in_list_at_cutoff,
947 if (fabs(drift) > drift_target)
959 *rlist = max(ir->rvdw, ir->rcoulomb) + ib1*resolution;