2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2012,2013,2014,2015,2016,2017,2018,2019, 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"
44 #include "gromacs/ewald/ewald_utils.h"
45 #include "gromacs/math/functions.h"
46 #include "gromacs/math/units.h"
47 #include "gromacs/math/vec.h"
48 #include "gromacs/mdlib/nb_verlet.h"
49 #include "gromacs/mdlib/nbnxn_simd.h"
50 #include "gromacs/mdlib/nbnxn_util.h"
51 #include "gromacs/mdtypes/inputrec.h"
52 #include "gromacs/mdtypes/md_enums.h"
53 #include "gromacs/topology/block.h"
54 #include "gromacs/topology/ifunc.h"
55 #include "gromacs/topology/topology.h"
56 #include "gromacs/utility/fatalerror.h"
57 #include "gromacs/utility/real.h"
58 #include "gromacs/utility/strconvert.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.
96 struct VerletbufAtomtype
98 atom_nonbonded_kinetic_prop_t prop
; /* non-bonded and kinetic atom prop. */
99 int n
; /* #atoms of this type in the system */
102 // Struct for derivatives of a non-bonded interaction potential
103 struct pot_derivatives_t
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
);
160 // Returns whether prop1 and prop2 are identical
162 atom_nonbonded_kinetic_prop_equal(const atom_nonbonded_kinetic_prop_t
&prop1
,
163 const atom_nonbonded_kinetic_prop_t
&prop2
)
165 return (prop1
.mass
== prop2
.mass
&&
166 prop1
.type
== prop2
.type
&&
167 prop1
.q
== prop2
.q
&&
168 prop1
.bConstr
== prop2
.bConstr
&&
169 prop1
.con_mass
== prop2
.con_mass
&&
170 prop1
.con_len
== prop2
.con_len
);
173 static void addAtomtype(std::vector
<VerletbufAtomtype
> *att
,
174 const atom_nonbonded_kinetic_prop_t
&prop
,
179 /* Ignore massless particles */
184 while (i
< att
->size() &&
185 !atom_nonbonded_kinetic_prop_equal(prop
, (*att
)[i
].prop
))
196 att
->push_back({ prop
, nmol
});
200 /* Returns the mass of atom atomIndex or 1 when setMassesToOne=true */
201 static real
getMass(const t_atoms
&atoms
,
207 return atoms
.atom
[atomIndex
].m
;
215 // Set the masses of a vsites in vsite_m and the non-linear vsite count in n_nonlin_vsite
216 static void get_vsite_masses(const gmx_moltype_t
&moltype
,
217 const gmx_ffparams_t
&ffparams
,
219 gmx::ArrayRef
<real
> vsite_m
,
222 GMX_RELEASE_ASSERT(n_nonlin_vsite
, "Expect a valid pointer");
226 /* Check for virtual sites, determine mass from constructing atoms */
227 for (const auto &ilist
: extractILists(moltype
.ilist
, IF_VSITE
))
229 for (size_t i
= 0; i
< ilist
.iatoms
.size(); i
+= ilistStride(ilist
))
231 const t_iparams
&ip
= ffparams
.iparams
[ilist
.iatoms
[i
]];
232 const int a1
= ilist
.iatoms
[i
+ 1];
234 if (ilist
.functionType
!= F_VSITEN
)
236 /* Only vsiten can have more than four
237 constructing atoms, so NRAL(ft) <= 5 */
238 const int maxj
= NRAL(ilist
.functionType
);
239 std::vector
<real
> cam(maxj
, 0);
240 GMX_ASSERT(maxj
<= 5, "This code expect at most 5 atoms in a vsite");
241 for (int j
= 1; j
< maxj
; j
++)
243 const int aj
= ilist
.iatoms
[i
+ 1 + j
];
244 cam
[j
] = getMass(moltype
.atoms
, aj
, setMassesToOne
);
247 cam
[j
] = vsite_m
[aj
];
249 /* A vsite should be constructed from normal atoms or
250 * vsites of lower complexity, which we have processed
251 * in a previous iteration.
253 GMX_ASSERT(cam
[j
] != 0, "We should have a non-zero mass");
256 switch (ilist
.functionType
)
260 vsite_m
[a1
] = (cam
[1]*cam
[2])/(cam
[2]*gmx::square(1 - ip
.vsite
.a
) + cam
[1]*gmx::square(ip
.vsite
.a
));
264 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
));
267 GMX_RELEASE_ASSERT(false, "VsiteN should not end up in this code path");
270 /* Use the mass of the lightest constructing atom.
271 * This is an approximation.
272 * If the distance of the virtual site to the
273 * constructing atom is less than all distances
274 * between constructing atoms, this is a safe
275 * over-estimate of the displacement of the vsite.
276 * This condition holds for all H mass replacement
277 * vsite constructions, except for SP2/3 groups.
278 * In SP3 groups one H will have a F_VSITE3
279 * construction, so even there the total drift
280 * estimate shouldn't be far off.
282 vsite_m
[a1
] = cam
[1];
283 for (int j
= 2; j
< maxj
; j
++)
285 vsite_m
[a1
] = std::min(vsite_m
[a1
], cam
[j
]);
295 int numConstructingAtoms
= ffparams
.iparams
[ilist
.iatoms
[i
]].vsiten
.n
;
296 for (int j
= 0; j
< 3*numConstructingAtoms
; j
+= 3)
298 int aj
= ilist
.iatoms
[i
+ j
+ 2];
299 real coeff
= ffparams
.iparams
[ilist
.iatoms
[i
+ j
]].vsiten
.a
;
301 if (moltype
.atoms
.atom
[aj
].ptype
== eptVSite
)
307 m_aj
= moltype
.atoms
.atom
[aj
].m
;
311 gmx_incons("The mass of a vsiten constructing atom is <= 0");
313 inv_mass
+= coeff
*coeff
/m_aj
;
315 vsite_m
[a1
] = 1/inv_mass
;
316 /* Correct the loop increment of i for processes more than 1 entry */
317 i
+= (numConstructingAtoms
- 1)*ilistStride(ilist
);
321 fprintf(debug
, "atom %4d %-20s mass %6.3f\n",
322 a1
, interaction_function
[ilist
.functionType
].longname
, vsite_m
[a1
]);
328 static std::vector
<VerletbufAtomtype
>
329 get_verlet_buffer_atomtypes(const gmx_mtop_t
*mtop
,
333 std::vector
<VerletbufAtomtype
> att
;
334 int ft
, i
, a1
, a2
, a3
, a
;
336 int n_nonlin_vsite_mol
;
338 if (n_nonlin_vsite
!= nullptr)
343 for (const gmx_molblock_t
&molblock
: mtop
->molblock
)
345 int nmol
= molblock
.nmol
;
346 const gmx_moltype_t
&moltype
= mtop
->moltype
[molblock
.type
];
347 const t_atoms
*atoms
= &moltype
.atoms
;
349 /* Check for constraints, as they affect the kinetic energy.
350 * For virtual sites we need the masses and geometry of
351 * the constructing atoms to determine their velocity distribution.
352 * Thus we need a list of properties for all atoms which
353 * we partially fill when looping over constraints.
355 std::vector
<atom_nonbonded_kinetic_prop_t
> prop(atoms
->nr
);
357 for (ft
= F_CONSTR
; ft
<= F_CONSTRNC
; ft
++)
359 const InteractionList
&il
= moltype
.ilist
[ft
];
361 for (i
= 0; i
< il
.size(); i
+= 1+NRAL(ft
))
363 ip
= &mtop
->ffparams
.iparams
[il
.iatoms
[i
]];
366 real mass1
= getMass(*atoms
, a1
, setMassesToOne
);
367 real mass2
= getMass(*atoms
, a2
, setMassesToOne
);
368 if (mass2
> prop
[a1
].con_mass
)
370 prop
[a1
].con_mass
= mass2
;
371 prop
[a1
].con_len
= ip
->constr
.dA
;
373 if (mass1
> prop
[a2
].con_mass
)
375 prop
[a2
].con_mass
= mass1
;
376 prop
[a2
].con_len
= ip
->constr
.dA
;
381 const InteractionList
&il
= moltype
.ilist
[F_SETTLE
];
383 for (i
= 0; i
< il
.size(); i
+= 1+NRAL(F_SETTLE
))
385 ip
= &mtop
->ffparams
.iparams
[il
.iatoms
[i
]];
389 /* Usually the mass of a1 (usually oxygen) is larger than a2/a3.
390 * If this is not the case, we overestimate the displacement,
391 * which leads to a larger buffer (ok since this is an exotic case).
393 prop
[a1
].con_mass
= getMass(*atoms
, a2
, setMassesToOne
);
394 prop
[a1
].con_len
= ip
->settle
.doh
;
396 prop
[a2
].con_mass
= getMass(*atoms
, a1
, setMassesToOne
);
397 prop
[a2
].con_len
= ip
->settle
.doh
;
399 prop
[a3
].con_mass
= getMass(*atoms
, a1
, setMassesToOne
);
400 prop
[a3
].con_len
= ip
->settle
.doh
;
403 std::vector
<real
> vsite_m(atoms
->nr
);
404 get_vsite_masses(moltype
,
408 &n_nonlin_vsite_mol
);
409 if (n_nonlin_vsite
!= nullptr)
411 *n_nonlin_vsite
+= nmol
*n_nonlin_vsite_mol
;
414 for (a
= 0; a
< atoms
->nr
; a
++)
416 if (atoms
->atom
[a
].ptype
== eptVSite
)
418 prop
[a
].mass
= vsite_m
[a
];
422 prop
[a
].mass
= getMass(*atoms
, a
, setMassesToOne
);
424 prop
[a
].type
= atoms
->atom
[a
].type
;
425 prop
[a
].q
= atoms
->atom
[a
].q
;
426 /* We consider an atom constrained, #DOF=2, when it is
427 * connected with constraints to (at least one) atom with
428 * a mass of more than 0.4x its own mass. This is not a critical
429 * parameter, since with roughly equal masses the unconstrained
430 * and constrained displacement will not differ much (and both
431 * overestimate the displacement).
433 prop
[a
].bConstr
= (prop
[a
].con_mass
> 0.4*prop
[a
].mass
);
435 addAtomtype(&att
, prop
[a
], nmol
);
441 for (size_t a
= 0; a
< att
.size(); a
++)
443 fprintf(debug
, "type %zu: m %5.2f t %d q %6.3f con %s con_m %5.3f con_l %5.3f n %d\n",
444 a
, att
[a
].prop
.mass
, att
[a
].prop
.type
, att
[a
].prop
.q
,
445 gmx::boolToString(att
[a
].prop
.bConstr
), att
[a
].prop
.con_mass
, att
[a
].prop
.con_len
,
453 /* This function computes two components of the estimate of the variance
454 * in the displacement of one atom in a system of two constrained atoms.
455 * Returns in sigma2_2d the variance due to rotation of the constrained
456 * atom around the atom to which it constrained.
457 * Returns in sigma2_3d the variance due to displacement of the COM
458 * of the whole system of the two constrained atoms.
460 * Note that we only take a single constraint (the one to the heaviest atom)
461 * into account. If an atom has multiple constraints, this will result in
462 * an overestimate of the displacement, which gives a larger drift and buffer.
464 void constrained_atom_sigma2(real kT_fac
,
465 const atom_nonbonded_kinetic_prop_t
*prop
,
469 /* Here we decompose the motion of a constrained atom into two
470 * components: rotation around the COM and translation of the COM.
473 /* Determine the variance of the arc length for the two rotational DOFs */
474 real massFraction
= prop
->con_mass
/(prop
->mass
+ prop
->con_mass
);
475 real sigma2_rot
= kT_fac
*massFraction
/prop
->mass
;
477 /* The distance from the atom to the COM, i.e. the rotational arm */
478 real comDistance
= prop
->con_len
*massFraction
;
480 /* The variance relative to the arm */
481 real sigma2_rel
= sigma2_rot
/gmx::square(comDistance
);
483 /* For sigma2_rel << 1 we don't notice the rotational effect and
484 * we have a normal, Gaussian displacement distribution.
485 * For larger sigma2_rel the displacement is much less, in fact it can
486 * not exceed 2*comDistance. We can calculate MSD/arm^2 as:
487 * integral_x=0-inf distance2(x) x/sigma2_rel exp(-x^2/(2 sigma2_rel)) dx
488 * where x is angular displacement and distance2(x) is the distance^2
489 * between points at angle 0 and x:
490 * distance2(x) = (sin(x) - sin(0))^2 + (cos(x) - cos(0))^2
491 * The limiting value of this MSD is 2, which is also the value for
492 * a uniform rotation distribution that would be reached at long time.
493 * The maximum is 2.5695 at sigma2_rel = 4.5119.
494 * We approximate this integral with a rational polynomial with
495 * coefficients from a Taylor expansion. This approximation is an
496 * overestimate for all values of sigma2_rel. Its maximum value
497 * of 2.6491 is reached at sigma2_rel = sqrt(45/2) = 4.7434.
498 * We keep the approximation constant after that.
499 * We use this approximate MSD as the variance for a Gaussian distribution.
501 * NOTE: For any sensible buffer tolerance this will result in a (large)
502 * overestimate of the buffer size, since the Gaussian has a long tail,
503 * whereas the actual distribution can not reach values larger than 2.
505 /* Coeffients obtained from a Taylor expansion */
506 const real a
= 1.0/3.0;
507 const real b
= 2.0/45.0;
509 /* Our approximation is constant after sigma2_rel = 1/sqrt(b) */
510 sigma2_rel
= std::min(sigma2_rel
, 1/std::sqrt(b
));
512 /* Compute the approximate sigma^2 for 2D motion due to the rotation */
513 *sigma2_2d
= gmx::square(comDistance
)*
514 sigma2_rel
/(1 + a
*sigma2_rel
+ b
*gmx::square(sigma2_rel
));
516 /* The constrained atom also moves (in 3D) with the COM of both atoms */
517 *sigma2_3d
= kT_fac
/(prop
->mass
+ prop
->con_mass
);
520 static void get_atom_sigma2(real kT_fac
,
521 const atom_nonbonded_kinetic_prop_t
*prop
,
527 /* Complicated constraint calculation in a separate function */
528 constrained_atom_sigma2(kT_fac
, prop
, sigma2_2d
, sigma2_3d
);
532 /* Unconstrained atom: trivial */
534 *sigma2_3d
= kT_fac
/prop
->mass
;
538 static void approx_2dof(real s2
, real x
, real
*shift
, real
*scale
)
540 /* A particle with 1 DOF constrained has 2 DOFs instead of 3.
541 * This code is also used for particles with multiple constraints,
542 * in which case we overestimate the displacement.
543 * The 2DOF distribution is sqrt(pi/2)*erfc(r/(sqrt(2)*s))/(2*s).
544 * We approximate this with scale*Gaussian(s,r+shift),
545 * by matching the distribution value and derivative at x.
546 * This is a tight overestimate for all r>=0 at any s and x.
550 ex
= std::exp(-x
*x
/(2*s2
));
551 er
= std::erfc(x
/std::sqrt(2*s2
));
553 *shift
= -x
+ std::sqrt(2*s2
/M_PI
)*ex
/er
;
554 *scale
= 0.5*M_PI
*std::exp(ex
*ex
/(M_PI
*er
*er
))*er
;
557 // Returns an (over)estimate of the energy drift for a single atom pair,
558 // given the kinetic properties, displacement variances and list buffer.
559 static real
energyDriftAtomPair(bool isConstrained_i
,
560 bool isConstrained_j
,
561 real s2
, real s2i_2d
, real s2j_2d
,
563 const pot_derivatives_t
*der
)
565 // For relatively small arguments erfc() is so small that if will be 0.0
566 // when stored in a float. We set an argument limit of 8 (Erfc(8)=1e-29),
567 // such that we can divide by erfc and have some space left for arithmetic.
568 const real erfc_arg_max
= 8.0;
575 if (rsh
*rsh
> 2*s2
*erfc_arg_max
*erfc_arg_max
)
577 // Below we calculate c_erfc = 0.5*erfc(rsh/sqrt(2*s2))
578 // When rsh/sqrt(2*s2) increases, this erfc will be the first
579 // result that underflows and becomes 0.0. To avoid this,
580 // we set c_exp=0 and c_erfc=0 for large arguments.
581 // This also avoids NaN in approx_2dof().
582 // In any relevant case this has no effect on the results,
583 // since c_exp < 6e-29, so the displacement is completely
584 // negligible for such atom pairs (and an overestimate).
585 // In nearly all use cases, there will be other atom pairs
586 // that contribute much more to the total, so zeroing
587 // this particular contribution has no effect at all.
593 /* For constraints: adapt r and scaling for the Gaussian */
598 approx_2dof(s2i_2d
, r_buffer
*s2i_2d
/s2
, &sh
, &sc
);
606 approx_2dof(s2j_2d
, r_buffer
*s2j_2d
/s2
, &sh
, &sc
);
611 /* Exact contribution of an atom pair with Gaussian displacement
612 * with sigma s to the energy drift for a potential with
613 * derivative -md and second derivative dd at the cut-off.
614 * The only catch is that for potentials that change sign
615 * near the cut-off there could be an unlucky compensation
616 * of positive and negative energy drift.
617 * Such potentials are extremely rare though.
619 * Note that pot has unit energy*length, as the linear
620 * atom density still needs to be put in.
622 c_exp
= std::exp(-rsh
*rsh
/(2*s2
))/std::sqrt(2*M_PI
);
623 c_erfc
= 0.5*std::erfc(rsh
/(std::sqrt(2*s2
)));
625 real s
= std::sqrt(s2
);
629 der
->md1
/2*((rsh2
+ s2
)*c_erfc
- rsh
*s
*c_exp
);
631 der
->d2
/6*(s
*(rsh2
+ 2*s2
)*c_exp
- rsh
*(rsh2
+ 3*s2
)*c_erfc
);
633 der
->md3
/24*((rsh2
*rsh2
+ 6*rsh2
*s2
+ 3*s2
*s2
)*c_erfc
- rsh
*s
*(rsh2
+ 5*s2
)*c_exp
);
635 return pot1
+ pot2
+ pot3
;
638 // Computes and returns an estimate of the energy drift for the whole system
639 static real
energyDrift(gmx::ArrayRef
<const VerletbufAtomtype
> att
,
640 const gmx_ffparams_t
*ffp
,
642 const pot_derivatives_t
*ljDisp
,
643 const pot_derivatives_t
*ljRep
,
644 const pot_derivatives_t
*elec
,
645 real rlj
, real rcoulomb
,
646 real rlist
, real boxvol
)
648 double drift_tot
= 0;
652 /* No atom displacements: no drift, avoid division by 0 */
656 // Here add up the contribution of all atom pairs in the system to
657 // (estimated) energy drift by looping over all atom type pairs.
658 for (int i
= 0; i
< att
.size(); i
++)
660 // Get the thermal displacement variance for the i-atom type
661 const atom_nonbonded_kinetic_prop_t
*prop_i
= &att
[i
].prop
;
663 get_atom_sigma2(kT_fac
, prop_i
, &s2i_2d
, &s2i_3d
);
665 for (int j
= i
; j
< att
.size(); j
++)
667 // Get the thermal displacement variance for the j-atom type
668 const atom_nonbonded_kinetic_prop_t
*prop_j
= &att
[j
].prop
;
670 get_atom_sigma2(kT_fac
, prop_j
, &s2j_2d
, &s2j_3d
);
672 /* Add up the up to four independent variances */
673 real s2
= s2i_2d
+ s2i_3d
+ s2j_2d
+ s2j_3d
;
675 // Set -V', V'' and -V''' at the cut-off for LJ */
676 real c6
= ffp
->iparams
[prop_i
->type
*ffp
->atnr
+ prop_j
->type
].lj
.c6
;
677 real c12
= ffp
->iparams
[prop_i
->type
*ffp
->atnr
+ prop_j
->type
].lj
.c12
;
678 pot_derivatives_t lj
;
679 lj
.md1
= c6
*ljDisp
->md1
+ c12
*ljRep
->md1
;
680 lj
.d2
= c6
*ljDisp
->d2
+ c12
*ljRep
->d2
;
681 lj
.md3
= c6
*ljDisp
->md3
+ c12
*ljRep
->md3
;
683 real pot_lj
= energyDriftAtomPair(prop_i
->bConstr
, prop_j
->bConstr
,
688 // Set -V' and V'' at the cut-off for Coulomb
689 pot_derivatives_t elec_qq
;
690 elec_qq
.md1
= elec
->md1
*prop_i
->q
*prop_j
->q
;
691 elec_qq
.d2
= elec
->d2
*prop_i
->q
*prop_j
->q
;
694 real pot_q
= energyDriftAtomPair(prop_i
->bConstr
, prop_j
->bConstr
,
699 // Note that attractive and repulsive potentials for individual
700 // pairs can partially cancel.
701 real pot
= pot_lj
+ pot_q
;
703 /* Multiply by the number of atom pairs */
706 pot
*= static_cast<double>(att
[i
].n
)*(att
[i
].n
- 1)/2;
710 pot
*= static_cast<double>(att
[i
].n
)*att
[j
].n
;
712 /* We need the line density to get the energy drift of the system.
713 * The effective average r^2 is close to (rlist+sigma)^2.
715 pot
*= 4*M_PI
*gmx::square(rlist
+ std::sqrt(s2
))/boxvol
;
717 /* Add the unsigned drift to avoid cancellation of errors */
718 drift_tot
+= std::abs(pot
);
725 // Returns the chance that a particle in a cluster is at distance rlist
726 // when the cluster is at distance rlist
727 static real
surface_frac(int cluster_size
, real particle_distance
, real rlist
)
731 if (rlist
< 0.5*particle_distance
)
733 /* We have non overlapping spheres */
737 /* Half the inter-particle distance relative to rlist */
738 d
= 0.5*particle_distance
/rlist
;
740 /* Determine the area of the surface at distance rlist to the closest
741 * particle, relative to surface of a sphere of radius rlist.
742 * The formulas below assume close to cubic cells for the pair search grid,
743 * which the pair search code tries to achieve.
744 * Note that in practice particle distances will not be delta distributed,
745 * but have some spread, often involving shorter distances,
746 * as e.g. O-H bonds in a water molecule. Thus the estimates below will
747 * usually be slightly too high and thus conservative.
749 switch (cluster_size
)
752 /* One particle: trivial */
756 /* Two particles: two spheres at fractional distance 2*a */
760 /* We assume a perfect, symmetric tetrahedron geometry.
761 * The surface around a tetrahedron is too complex for a full
762 * analytical solution, so we use a Taylor expansion.
764 area_rel
= (1.0 + 1/M_PI
*(6*std::acos(1/std::sqrt(3))*d
+
765 std::sqrt(3)*d
*d
*(1.0 +
768 83.0/756.0*d
*d
*d
*d
*d
*d
)));
771 gmx_incons("surface_frac called with unsupported cluster_size");
774 return area_rel
/cluster_size
;
777 /* Returns the negative of the third derivative of a potential r^-p
778 * with a force-switch function, evaluated at the cut-off rc.
780 static real
md3_force_switch(real p
, real rswitch
, real rc
)
782 /* The switched force function is:
783 * p*r^-(p+1) + a*(r - rswitch)^2 + b*(r - rswitch)^3
786 real md3_pot
, md3_sw
;
788 a
= -((p
+ 4)*rc
- (p
+ 1)*rswitch
)/(pow(rc
, p
+2)*gmx::square(rc
-rswitch
));
789 b
= ((p
+ 3)*rc
- (p
+ 1)*rswitch
)/(pow(rc
, p
+2)*gmx::power3(rc
-rswitch
));
791 md3_pot
= (p
+ 2)*(p
+ 1)*p
*pow(rc
, p
+3);
792 md3_sw
= 2*a
+ 6*b
*(rc
- rswitch
);
794 return md3_pot
+ md3_sw
;
797 /* Returns the variance of the atomic displacement over timePeriod.
799 * Note: When not using BD with a non-mass dependendent friction coefficient,
800 * the return value still needs to be divided by the particle mass.
802 static real
displacementVariance(const t_inputrec
&ir
,
810 /* Get the displacement distribution from the random component only.
811 * With accurate integration the systematic (force) displacement
812 * should be negligible (unless nstlist is extremely large, which
813 * you wouldn't do anyhow).
815 kT_fac
= 2*BOLTZ
*temperature
*timePeriod
;
818 /* This is directly sigma^2 of the displacement */
819 kT_fac
/= ir
.bd_fric
;
823 /* Per group tau_t is not implemented yet, use the maximum */
824 real tau_t
= ir
.opts
.tau_t
[0];
825 for (int i
= 1; i
< ir
.opts
.ngtc
; i
++)
827 tau_t
= std::max(tau_t
, ir
.opts
.tau_t
[i
]);
831 /* This kT_fac needs to be divided by the mass to get sigma^2 */
836 kT_fac
= BOLTZ
*temperature
*gmx::square(timePeriod
);
842 /* Returns the largest sigma of the Gaussian displacement over all particle
843 * types. This ignores constraints, so is an overestimate.
845 static real
maxSigma(real kT_fac
,
846 gmx::ArrayRef
<const VerletbufAtomtype
> att
)
848 GMX_ASSERT(!att
.empty(), "We should have at least one type");
849 real smallestMass
= att
[0].prop
.mass
;
850 for (int i
= 1; i
< att
.size(); i
++)
852 smallestMass
= std::min(smallestMass
, att
[i
].prop
.mass
);
855 return 2*std::sqrt(kT_fac
/smallestMass
);
858 void calc_verlet_buffer_size(const gmx_mtop_t
*mtop
, real boxvol
,
859 const t_inputrec
*ir
,
862 real reference_temperature
,
863 const VerletbufListSetup
*list_setup
,
870 real particle_distance
;
871 real nb_clust_frac_pairs_not_in_list_at_cutoff
;
878 if (!EI_DYNAMICS(ir
->eI
))
880 gmx_incons("Can only determine the Verlet buffer size for integrators that perform dynamics");
882 if (ir
->verletbuf_tol
<= 0)
884 gmx_incons("The Verlet buffer tolerance needs to be larger than zero");
887 if (reference_temperature
< 0)
889 /* We use the maximum temperature with multiple T-coupl groups.
890 * We could use a per particle temperature, but since particles
891 * interact, this might underestimate the buffer size.
893 reference_temperature
= maxReferenceTemperature(*ir
);
895 GMX_RELEASE_ASSERT(reference_temperature
>= 0, "Without T-coupling we should not end up here");
898 /* Resolution of the buffer size */
901 env
= getenv("GMX_VERLET_BUFFER_RES");
904 sscanf(env
, "%lf", &resolution
);
907 /* In an atom wise pair-list there would be no pairs in the list
908 * beyond the pair-list cut-off.
909 * However, we use a pair-list of groups vs groups of atoms.
910 * For groups of 4 atoms, the parallelism of SSE instructions, only
911 * 10% of the atoms pairs are not in the list just beyond the cut-off.
912 * As this percentage increases slowly compared to the decrease of the
913 * Gaussian displacement distribution over this range, we can simply
914 * reduce the drift by this fraction.
915 * For larger groups, e.g. of 8 atoms, this fraction will be lower,
916 * so then buffer size will be on the conservative (large) side.
918 * Note that the formulas used here do not take into account
919 * cancellation of errors which could occur by missing both
920 * attractive and repulsive interactions.
922 * The only major assumption is homogeneous particle distribution.
923 * For an inhomogeneous system, such as a liquid-vapor system,
924 * the buffer will be underestimated. The actual energy drift
925 * will be higher by the factor: local/homogeneous particle density.
927 * The results of this estimate have been checked againt simulations.
928 * In most cases the real drift differs by less than a factor 2.
931 /* Worst case assumption: HCP packing of particles gives largest distance */
932 particle_distance
= std::cbrt(boxvol
*std::sqrt(2)/mtop
->natoms
);
934 /* TODO: Obtain masses through (future) integrator functionality
935 * to avoid scattering the code with (or forgetting) checks.
937 const bool setMassesToOne
= (ir
->eI
== eiBD
&& ir
->bd_fric
> 0);
939 get_verlet_buffer_atomtypes(mtop
, setMassesToOne
, n_nonlin_vsite
);
940 GMX_ASSERT(!att
.empty(), "We expect at least one type");
944 fprintf(debug
, "particle distance assuming HCP packing: %f nm\n",
946 fprintf(debug
, "energy drift atom types: %zu\n", att
.size());
949 pot_derivatives_t ljDisp
= { 0, 0, 0 };
950 pot_derivatives_t ljRep
= { 0, 0, 0 };
951 real repPow
= mtop
->ffparams
.reppow
;
953 if (ir
->vdwtype
== evdwCUT
)
955 real sw_range
, md3_pswf
;
957 switch (ir
->vdw_modifier
)
960 case eintmodPOTSHIFT
:
961 /* -dV/dr of -r^-6 and r^-reppow */
962 ljDisp
.md1
= -6*std::pow(ir
->rvdw
, -7.0);
963 ljRep
.md1
= repPow
*std::pow(ir
->rvdw
, -(repPow
+ 1));
964 /* The contribution of the higher derivatives is negligible */
966 case eintmodFORCESWITCH
:
967 /* At the cut-off: V=V'=V''=0, so we use only V''' */
968 ljDisp
.md3
= -md3_force_switch(6.0, ir
->rvdw_switch
, ir
->rvdw
);
969 ljRep
.md3
= md3_force_switch(repPow
, ir
->rvdw_switch
, ir
->rvdw
);
971 case eintmodPOTSWITCH
:
972 /* At the cut-off: V=V'=V''=0.
973 * V''' is given by the original potential times
974 * the third derivative of the switch function.
976 sw_range
= ir
->rvdw
- ir
->rvdw_switch
;
977 md3_pswf
= 60.0/gmx::power3(sw_range
);
979 ljDisp
.md3
= -std::pow(ir
->rvdw
, -6.0 )*md3_pswf
;
980 ljRep
.md3
= std::pow(ir
->rvdw
, -repPow
)*md3_pswf
;
983 gmx_incons("Unimplemented VdW modifier");
986 else if (EVDW_PME(ir
->vdwtype
))
988 real b
= calc_ewaldcoeff_lj(ir
->rvdw
, ir
->ewald_rtol_lj
);
994 // -dV/dr of g(br)*r^-6 [where g(x) = exp(-x^2)(1+x^2+x^4/2),
995 // see LJ-PME equations in manual] and r^-reppow
996 ljDisp
.md1
= -std::exp(-br2
)*(br6
+ 3.0*br4
+ 6.0*br2
+ 6.0)*std::pow(r
, -7.0);
997 ljRep
.md1
= repPow
*pow(r
, -(repPow
+ 1));
998 // The contribution of the higher derivatives is negligible
1002 gmx_fatal(FARGS
, "Energy drift calculation is only implemented for plain cut-off Lennard-Jones interactions");
1005 elfac
= ONE_4PI_EPS0
/ir
->epsilon_r
;
1007 // Determine the 1st and 2nd derivative for the electostatics
1008 pot_derivatives_t elec
= { 0, 0, 0 };
1010 if (ir
->coulombtype
== eelCUT
|| EEL_RF(ir
->coulombtype
))
1014 if (ir
->coulombtype
== eelCUT
)
1021 eps_rf
= ir
->epsilon_rf
/ir
->epsilon_r
;
1024 k_rf
= (eps_rf
- ir
->epsilon_r
)/( gmx::power3(ir
->rcoulomb
) * (2*eps_rf
+ ir
->epsilon_r
) );
1028 /* epsilon_rf = infinity */
1029 k_rf
= 0.5/gmx::power3(ir
->rcoulomb
);
1035 elec
.md1
= elfac
*(1.0/gmx::square(ir
->rcoulomb
) - 2*k_rf
*ir
->rcoulomb
);
1037 elec
.d2
= elfac
*(2.0/gmx::power3(ir
->rcoulomb
) + 2*k_rf
);
1039 else if (EEL_PME(ir
->coulombtype
) || ir
->coulombtype
== eelEWALD
)
1043 b
= calc_ewaldcoeff_q(ir
->rcoulomb
, ir
->ewald_rtol
);
1046 elec
.md1
= elfac
*(b
*std::exp(-br
*br
)*M_2_SQRTPI
/rc
+ std::erfc(br
)/(rc
*rc
));
1047 elec
.d2
= elfac
/(rc
*rc
)*(2*b
*(1 + br
*br
)*std::exp(-br
*br
)*M_2_SQRTPI
+ 2*std::erfc(br
)/rc
);
1051 gmx_fatal(FARGS
, "Energy drift calculation is only implemented for Reaction-Field and Ewald electrostatics");
1054 /* Determine the variance of the atomic displacement
1055 * over list_lifetime steps: kT_fac
1056 * For inertial dynamics (not Brownian dynamics) the mass factor
1057 * is not included in kT_fac, it is added later.
1059 const real kT_fac
= displacementVariance(*ir
, reference_temperature
,
1060 list_lifetime
*ir
->delta_t
);
1064 fprintf(debug
, "Derivatives of non-bonded potentials at the cut-off:\n");
1065 fprintf(debug
, "LJ disp. -V' %9.2e V'' %9.2e -V''' %9.2e\n", ljDisp
.md1
, ljDisp
.d2
, ljDisp
.md3
);
1066 fprintf(debug
, "LJ rep. -V' %9.2e V'' %9.2e -V''' %9.2e\n", ljRep
.md1
, ljRep
.d2
, ljRep
.md3
);
1067 fprintf(debug
, "Electro. -V' %9.2e V'' %9.2e\n", elec
.md1
, elec
.d2
);
1068 fprintf(debug
, "sqrt(kT_fac) %f\n", std::sqrt(kT_fac
));
1071 /* Search using bisection */
1073 /* The drift will be neglible at 5 times the max sigma */
1074 ib1
= static_cast<int>(5*maxSigma(kT_fac
, att
)/resolution
) + 1;
1075 while (ib1
- ib0
> 1)
1079 rl
= std::max(ir
->rvdw
, ir
->rcoulomb
) + rb
;
1081 /* Calculate the average energy drift at the last step
1082 * of the nstlist steps at which the pair-list is used.
1084 drift
= energyDrift(att
, &mtop
->ffparams
,
1086 &ljDisp
, &ljRep
, &elec
,
1087 ir
->rvdw
, ir
->rcoulomb
,
1090 /* Correct for the fact that we are using a Ni x Nj particle pair list
1091 * and not a 1 x 1 particle pair list. This reduces the drift.
1093 /* We don't have a formula for 8 (yet), use 4 which is conservative */
1094 nb_clust_frac_pairs_not_in_list_at_cutoff
=
1095 surface_frac(std::min(list_setup
->cluster_size_i
, 4),
1096 particle_distance
, rl
)*
1097 surface_frac(std::min(list_setup
->cluster_size_j
, 4),
1098 particle_distance
, rl
);
1099 drift
*= nb_clust_frac_pairs_not_in_list_at_cutoff
;
1101 /* Convert the drift to drift per unit time per atom */
1102 drift
/= nstlist
*ir
->delta_t
*mtop
->natoms
;
1106 fprintf(debug
, "ib %3d %3d %3d rb %.3f %dx%d fac %.3f drift %.1e\n",
1108 list_setup
->cluster_size_i
, list_setup
->cluster_size_j
,
1109 nb_clust_frac_pairs_not_in_list_at_cutoff
,
1113 if (std::abs(drift
) > ir
->verletbuf_tol
)
1123 *rlist
= std::max(ir
->rvdw
, ir
->rcoulomb
) + ib1
*resolution
;
1126 /* Returns the pairlist buffer size for use as a minimum buffer size
1128 * Note that this is a rather crude estimate. It is ok for a buffer
1129 * set for good energy conservation or RF electrostatics. But it is
1130 * too small with PME and the buffer set with the default tolerance.
1132 static real
minCellSizeFromPairlistBuffer(const t_inputrec
&ir
)
1134 return ir
.rlist
- std::max(ir
.rvdw
, ir
.rcoulomb
);
1138 chanceOfAtomCrossingCell(gmx::ArrayRef
<const VerletbufAtomtype
> atomtypes
,
1142 /* We assume atoms are distributed uniformly over the cell width.
1143 * Once an atom has moved by more than the cellSize (as passed
1144 * as the buffer argument to energyDriftAtomPair() below),
1145 * the chance of crossing the boundary of the neighbor cell
1146 * thus increases as 1/cellSize with the additional displacement
1147 * on top of cellSize. We thus create a linear interaction with
1148 * derivative = -1/cellSize. Using this in the energyDriftAtomPair
1149 * function will return the chance of crossing the next boundary.
1151 const pot_derivatives_t boundaryInteraction
= { 1/cellSize
, 0, 0 };
1154 for (const VerletbufAtomtype
&att
: atomtypes
)
1156 const atom_nonbonded_kinetic_prop_t
&propAtom
= att
.prop
;
1159 get_atom_sigma2(kT_fac
, &propAtom
, &s2_2d
, &s2_3d
);
1161 real chancePerAtom
= energyDriftAtomPair(propAtom
.bConstr
, false,
1162 s2_2d
+ s2_3d
, s2_2d
, 0,
1164 &boundaryInteraction
);
1166 if (propAtom
.bConstr
)
1168 /* energyDriftAtomPair() uses an unlimited Gaussian displacement
1169 * distribution for constrained atoms, whereas they can
1170 * actually not move more than the COM of the two constrained
1171 * atoms plus twice the distance from the COM.
1172 * Use this maximum, limited displacement when this results in
1173 * a smaller chance (note that this is still an overestimate).
1175 real massFraction
= propAtom
.con_mass
/(propAtom
.mass
+ propAtom
.con_mass
);
1176 real comDistance
= propAtom
.con_len
*massFraction
;
1178 real chanceWithMaxDistance
=
1179 energyDriftAtomPair(false, false,
1181 cellSize
- 2*comDistance
,
1182 &boundaryInteraction
);
1183 chancePerAtom
= std::min(chancePerAtom
, chanceWithMaxDistance
);
1186 /* Take into account the line density of the boundary */
1187 chancePerAtom
/= cellSize
;
1189 chance
+= att
.n
*chancePerAtom
;
1195 /* Struct for storing constraint properties of atoms */
1196 struct AtomConstraintProps
1198 void addConstraint(real length
)
1200 numConstraints
+= 1;
1201 sumLengths
+= length
;
1204 int numConstraints
= 0; /* The number of constraints of an atom */
1205 real sumLengths
= 0; /* The sum of constraint length over the constraints */
1208 /* Constructs and returns a list of constraint properties per atom */
1209 static std::vector
<AtomConstraintProps
>
1210 getAtomConstraintProps(const gmx_moltype_t
&moltype
,
1211 const gmx_ffparams_t
&ffparams
)
1213 const t_atoms
&atoms
= moltype
.atoms
;
1214 std::vector
<AtomConstraintProps
> props(atoms
.nr
);
1216 for (const auto &ilist
: extractILists(moltype
.ilist
, IF_CONSTRAINT
))
1218 // Settles are handled separately
1219 if (ilist
.functionType
== F_SETTLE
)
1224 for (size_t i
= 0; i
< ilist
.iatoms
.size(); i
+= ilistStride(ilist
))
1226 int type
= ilist
.iatoms
[i
];
1227 int a1
= ilist
.iatoms
[i
+ 1];
1228 int a2
= ilist
.iatoms
[i
+ 2];
1229 real length
= ffparams
.iparams
[type
].constr
.dA
;
1230 props
[a1
].addConstraint(length
);
1231 props
[a2
].addConstraint(length
);
1238 /* Return the chance of at least one update group in a molecule crossing a cell of size cellSize */
1240 chanceOfUpdateGroupCrossingCell(const gmx_moltype_t
&moltype
,
1241 const gmx_ffparams_t
&ffparams
,
1242 const gmx::RangePartitioning
&updateGrouping
,
1246 const t_atoms
&atoms
= moltype
.atoms
;
1247 GMX_ASSERT(updateGrouping
.fullRange().end() == atoms
.nr
, "The update groups should match the molecule type");
1249 const pot_derivatives_t boundaryInteraction
= { 1/cellSize
, 0, 0 };
1251 const auto atomConstraintProps
= getAtomConstraintProps(moltype
, ffparams
);
1254 for (int group
= 0; group
< updateGrouping
.numBlocks(); group
++)
1256 const auto &block
= updateGrouping
.block(group
);
1257 /* Determine the number of atoms with constraints and the mass of the COG */
1258 int numAtomsWithConstraints
= 0;
1260 for (const int atom
: block
)
1262 if (atomConstraintProps
[atom
].numConstraints
> 0)
1264 numAtomsWithConstraints
++;
1266 massSum
+= moltype
.atoms
.atom
[atom
].m
;
1268 /* Determine the maximum possible distance between the center of mass
1269 * and the center of geometry of the update group
1271 real maxComCogDistance
= 0;
1272 if (numAtomsWithConstraints
== 2)
1274 for (const int atom
: block
)
1276 if (atomConstraintProps
[atom
].numConstraints
> 0)
1278 GMX_ASSERT(atomConstraintProps
[atom
].numConstraints
== 1,
1279 "Two atoms should be connected by one constraint");
1280 maxComCogDistance
= std::abs(atoms
.atom
[atom
].m
/massSum
- 0.5)*atomConstraintProps
[atom
].sumLengths
;
1285 else if (numAtomsWithConstraints
> 2)
1287 for (const int atom
: block
)
1289 if (atomConstraintProps
[atom
].numConstraints
== numAtomsWithConstraints
- 1)
1291 real comCogDistance
= atomConstraintProps
[atom
].sumLengths
/numAtomsWithConstraints
;
1292 maxComCogDistance
= std::max(maxComCogDistance
, comCogDistance
);
1296 else if (block
.size() > 1)
1298 // All normal atoms must be connected by SETTLE
1299 for (const int atom
: block
)
1301 const auto &ilist
= moltype
.ilist
[F_SETTLE
];
1302 GMX_RELEASE_ASSERT(ilist
.size() > 0, "There should be at least one settle in this moltype");
1303 for (int i
= 0; i
< ilist
.size(); i
+= 1 + NRAL(F_SETTLE
))
1305 if (atom
== ilist
.iatoms
[i
+ 1])
1307 const t_iparams
&iparams
= ffparams
.iparams
[ilist
.iatoms
[i
]];
1308 real dOH
= iparams
.settle
.doh
;
1309 real dHH
= iparams
.settle
.dhh
;
1310 real dOMidH
= std::sqrt(dOH
*dOH
- 0.25_real
*dHH
*dHH
);
1311 maxComCogDistance
= std::abs(atoms
.atom
[atom
].m
/massSum
- 1.0_real
/3.0_real
)*dOMidH
;
1316 real s2_3d
= kT_fac
/massSum
;
1317 chance
+= energyDriftAtomPair(false, false,
1319 cellSize
- 2*maxComCogDistance
,
1320 &boundaryInteraction
);
1326 /* Return the chance of at least one update group in the system crossing a cell of size cellSize */
1328 chanceOfUpdateGroupCrossingCell(const gmx_mtop_t
&mtop
,
1329 PartitioningPerMoltype updateGrouping
,
1333 GMX_RELEASE_ASSERT(static_cast<size_t>(updateGrouping
.size()) == mtop
.moltype
.size(),
1334 "The update groups should match the topology");
1337 for (const gmx_molblock_t
&molblock
: mtop
.molblock
)
1339 const gmx_moltype_t
&moltype
= mtop
.moltype
[molblock
.type
];
1341 molblock
.nmol
*chanceOfUpdateGroupCrossingCell(moltype
, mtop
.ffparams
, updateGrouping
[molblock
.type
],
1349 minCellSizeForAtomDisplacement(const gmx_mtop_t
&mtop
,
1350 const t_inputrec
&ir
,
1351 PartitioningPerMoltype updateGrouping
,
1352 real chanceRequested
)
1354 if (!EI_DYNAMICS(ir
.eI
) || (EI_MD(ir
.eI
) && ir
.etc
== etcNO
))
1356 return minCellSizeFromPairlistBuffer(ir
);
1359 /* We use the maximum temperature with multiple T-coupl groups.
1360 * We could use a per particle temperature, but since particles
1361 * interact, this might underestimate the displacements.
1363 const real temperature
= maxReferenceTemperature(ir
);
1365 const bool setMassesToOne
= (ir
.eI
== eiBD
&& ir
.bd_fric
> 0);
1367 const auto atomtypes
= get_verlet_buffer_atomtypes(&mtop
, setMassesToOne
, nullptr);
1369 const real kT_fac
= displacementVariance(ir
, temperature
,
1370 ir
.nstlist
*ir
.delta_t
);
1372 /* Resolution of the cell size */
1373 real resolution
= 0.001;
1375 /* Search using bisection, avoid 0 and start at 1 */
1377 /* The chance will be neglible at 10 times the max sigma */
1378 int ib1
= int(10*maxSigma(kT_fac
, atomtypes
)/resolution
) + 1;
1380 while (ib1
- ib0
> 1)
1382 int ib
= (ib0
+ ib1
)/2;
1383 cellSize
= ib
*resolution
;
1386 if (updateGrouping
.empty())
1388 chance
= chanceOfAtomCrossingCell(atomtypes
, kT_fac
, cellSize
);
1392 chance
= chanceOfUpdateGroupCrossingCell(mtop
, updateGrouping
, kT_fac
, cellSize
);
1395 /* Note: chance is for every nstlist steps */
1396 if (chance
> chanceRequested
*ir
.nstlist
)