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/mdtypes/inputrec.h"
49 #include "gromacs/mdtypes/md_enums.h"
50 #include "gromacs/nbnxm/nbnxm.h"
51 #include "gromacs/nbnxm/nbnxm_geometry.h"
52 #include "gromacs/nbnxm/nbnxm_simd.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(Nbnxm::KernelType 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
= Nbnxm::IClusterSizePerKernelType
[nbnxnKernelType
];
118 listSetup
.cluster_size_j
= Nbnxm::JClusterSizePerKernelType
[nbnxnKernelType
];
120 if (!Nbnxm::kernelTypeUsesSimplePairlist(nbnxnKernelType
))
122 /* The GPU kernels (except for OpenCL) split the j-clusters in two halves */
123 listSetup
.cluster_size_j
/= 2;
129 VerletbufListSetup
verletbufGetSafeListSetup(ListSetupType listType
)
131 /* When calling this function we often don't know which kernel type we
132 * are going to use. We choose the kernel type with the smallest possible
133 * i- and j-cluster sizes, so we potentially overestimate, but never
134 * underestimate, the buffer drift.
136 Nbnxm::KernelType nbnxnKernelType
;
138 if (listType
== ListSetupType::Gpu
)
140 nbnxnKernelType
= Nbnxm::KernelType::Gpu8x8x8
;
142 else if (GMX_SIMD
&& listType
== ListSetupType::CpuSimdWhenSupported
)
144 #ifdef GMX_NBNXN_SIMD_2XNN
145 /* We use the smallest cluster size to be on the safe side */
146 nbnxnKernelType
= Nbnxm::KernelType::Cpu4xN_Simd_2xNN
;
148 nbnxnKernelType
= Nbnxm::KernelType::Cpu4xN_Simd_4xN
;
153 nbnxnKernelType
= Nbnxm::KernelType::Cpu4x4_PlainC
;
156 return verletbufGetListSetup(nbnxnKernelType
);
159 // Returns whether prop1 and prop2 are identical
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 addAtomtype(std::vector
<VerletbufAtomtype
> *att
,
173 const atom_nonbonded_kinetic_prop_t
&prop
,
178 /* Ignore massless particles */
183 while (i
< att
->size() &&
184 !atom_nonbonded_kinetic_prop_equal(prop
, (*att
)[i
].prop
))
195 att
->push_back({ prop
, nmol
});
199 /* Returns the mass of atom atomIndex or 1 when setMassesToOne=true */
200 static real
getMass(const t_atoms
&atoms
,
206 return atoms
.atom
[atomIndex
].m
;
214 // Set the masses of a vsites in vsite_m and the non-linear vsite count in n_nonlin_vsite
215 static void get_vsite_masses(const gmx_moltype_t
&moltype
,
216 const gmx_ffparams_t
&ffparams
,
218 gmx::ArrayRef
<real
> vsite_m
)
220 int numNonlinearVsites
= 0;
222 /* Check for virtual sites, determine mass from constructing atoms */
223 for (const auto &ilist
: extractILists(moltype
.ilist
, IF_VSITE
))
225 for (size_t i
= 0; i
< ilist
.iatoms
.size(); i
+= ilistStride(ilist
))
227 const t_iparams
&ip
= ffparams
.iparams
[ilist
.iatoms
[i
]];
228 const int a1
= ilist
.iatoms
[i
+ 1];
230 if (ilist
.functionType
!= F_VSITEN
)
232 /* Only vsiten can have more than four
233 constructing atoms, so NRAL(ft) <= 5 */
234 const int maxj
= NRAL(ilist
.functionType
);
235 std::vector
<real
> cam(maxj
, 0);
236 GMX_ASSERT(maxj
<= 5, "This code expect at most 5 atoms in a vsite");
237 for (int j
= 1; j
< maxj
; j
++)
239 const int aj
= ilist
.iatoms
[i
+ 1 + j
];
240 cam
[j
] = getMass(moltype
.atoms
, aj
, setMassesToOne
);
243 cam
[j
] = vsite_m
[aj
];
245 /* A vsite should be constructed from normal atoms or
246 * vsites of lower complexity, which we have processed
247 * in a previous iteration.
249 GMX_ASSERT(cam
[j
] != 0, "We should have a non-zero mass");
252 switch (ilist
.functionType
)
256 vsite_m
[a1
] = (cam
[1]*cam
[2])/(cam
[2]*gmx::square(1 - ip
.vsite
.a
) + cam
[1]*gmx::square(ip
.vsite
.a
));
260 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
));
263 GMX_RELEASE_ASSERT(false, "VsiteN should not end up in this code path");
266 /* Use the mass of the lightest constructing atom.
267 * This is an approximation.
268 * If the distance of the virtual site to the
269 * constructing atom is less than all distances
270 * between constructing atoms, this is a safe
271 * over-estimate of the displacement of the vsite.
272 * This condition holds for all H mass replacement
273 * vsite constructions, except for SP2/3 groups.
274 * In SP3 groups one H will have a F_VSITE3
275 * construction, so even there the total drift
276 * estimate shouldn't be far off.
278 vsite_m
[a1
] = cam
[1];
279 for (int j
= 2; j
< maxj
; j
++)
281 vsite_m
[a1
] = std::min(vsite_m
[a1
], cam
[j
]);
283 numNonlinearVsites
++;
291 int numConstructingAtoms
= ffparams
.iparams
[ilist
.iatoms
[i
]].vsiten
.n
;
292 for (int j
= 0; j
< 3*numConstructingAtoms
; j
+= 3)
294 int aj
= ilist
.iatoms
[i
+ j
+ 2];
295 real coeff
= ffparams
.iparams
[ilist
.iatoms
[i
+ j
]].vsiten
.a
;
297 if (moltype
.atoms
.atom
[aj
].ptype
== eptVSite
)
303 m_aj
= moltype
.atoms
.atom
[aj
].m
;
307 gmx_incons("The mass of a vsiten constructing atom is <= 0");
309 inv_mass
+= coeff
*coeff
/m_aj
;
311 vsite_m
[a1
] = 1/inv_mass
;
312 /* Correct the loop increment of i for processes more than 1 entry */
313 i
+= (numConstructingAtoms
- 1)*ilistStride(ilist
);
317 fprintf(debug
, "atom %4d %-20s mass %6.3f\n",
318 a1
, interaction_function
[ilist
.functionType
].longname
, vsite_m
[a1
]);
323 if (debug
&& numNonlinearVsites
> 0)
325 fprintf(debug
, "The molecule type has %d non-linear virtual constructions\n",
330 static std::vector
<VerletbufAtomtype
>
331 getVerletBufferAtomtypes(const gmx_mtop_t
&mtop
,
332 const bool setMassesToOne
)
334 std::vector
<VerletbufAtomtype
> att
;
335 int ft
, i
, a1
, a2
, a3
, a
;
338 for (const gmx_molblock_t
&molblock
: mtop
.molblock
)
340 int nmol
= molblock
.nmol
;
341 const gmx_moltype_t
&moltype
= mtop
.moltype
[molblock
.type
];
342 const t_atoms
*atoms
= &moltype
.atoms
;
344 /* Check for constraints, as they affect the kinetic energy.
345 * For virtual sites we need the masses and geometry of
346 * the constructing atoms to determine their velocity distribution.
347 * Thus we need a list of properties for all atoms which
348 * we partially fill when looping over constraints.
350 std::vector
<atom_nonbonded_kinetic_prop_t
> prop(atoms
->nr
);
352 for (ft
= F_CONSTR
; ft
<= F_CONSTRNC
; ft
++)
354 const InteractionList
&il
= moltype
.ilist
[ft
];
356 for (i
= 0; i
< il
.size(); i
+= 1+NRAL(ft
))
358 ip
= &mtop
.ffparams
.iparams
[il
.iatoms
[i
]];
361 real mass1
= getMass(*atoms
, a1
, setMassesToOne
);
362 real mass2
= getMass(*atoms
, a2
, setMassesToOne
);
363 if (mass2
> prop
[a1
].con_mass
)
365 prop
[a1
].con_mass
= mass2
;
366 prop
[a1
].con_len
= ip
->constr
.dA
;
368 if (mass1
> prop
[a2
].con_mass
)
370 prop
[a2
].con_mass
= mass1
;
371 prop
[a2
].con_len
= ip
->constr
.dA
;
376 const InteractionList
&il
= moltype
.ilist
[F_SETTLE
];
378 for (i
= 0; i
< il
.size(); i
+= 1+NRAL(F_SETTLE
))
380 ip
= &mtop
.ffparams
.iparams
[il
.iatoms
[i
]];
384 /* Usually the mass of a1 (usually oxygen) is larger than a2/a3.
385 * If this is not the case, we overestimate the displacement,
386 * which leads to a larger buffer (ok since this is an exotic case).
388 prop
[a1
].con_mass
= getMass(*atoms
, a2
, setMassesToOne
);
389 prop
[a1
].con_len
= ip
->settle
.doh
;
391 prop
[a2
].con_mass
= getMass(*atoms
, a1
, setMassesToOne
);
392 prop
[a2
].con_len
= ip
->settle
.doh
;
394 prop
[a3
].con_mass
= getMass(*atoms
, a1
, setMassesToOne
);
395 prop
[a3
].con_len
= ip
->settle
.doh
;
398 std::vector
<real
> vsite_m(atoms
->nr
);
399 get_vsite_masses(moltype
,
404 for (a
= 0; a
< atoms
->nr
; a
++)
406 if (atoms
->atom
[a
].ptype
== eptVSite
)
408 prop
[a
].mass
= vsite_m
[a
];
412 prop
[a
].mass
= getMass(*atoms
, a
, setMassesToOne
);
414 prop
[a
].type
= atoms
->atom
[a
].type
;
415 prop
[a
].q
= atoms
->atom
[a
].q
;
416 /* We consider an atom constrained, #DOF=2, when it is
417 * connected with constraints to (at least one) atom with
418 * a mass of more than 0.4x its own mass. This is not a critical
419 * parameter, since with roughly equal masses the unconstrained
420 * and constrained displacement will not differ much (and both
421 * overestimate the displacement).
423 prop
[a
].bConstr
= (prop
[a
].con_mass
> 0.4*prop
[a
].mass
);
425 addAtomtype(&att
, prop
[a
], nmol
);
431 for (size_t a
= 0; a
< att
.size(); a
++)
433 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",
434 a
, att
[a
].prop
.mass
, att
[a
].prop
.type
, att
[a
].prop
.q
,
435 gmx::boolToString(att
[a
].prop
.bConstr
), att
[a
].prop
.con_mass
, att
[a
].prop
.con_len
,
443 /* This function computes two components of the estimate of the variance
444 * in the displacement of one atom in a system of two constrained atoms.
445 * Returns in sigma2_2d the variance due to rotation of the constrained
446 * atom around the atom to which it constrained.
447 * Returns in sigma2_3d the variance due to displacement of the COM
448 * of the whole system of the two constrained atoms.
450 * Note that we only take a single constraint (the one to the heaviest atom)
451 * into account. If an atom has multiple constraints, this will result in
452 * an overestimate of the displacement, which gives a larger drift and buffer.
454 void constrained_atom_sigma2(real kT_fac
,
455 const atom_nonbonded_kinetic_prop_t
*prop
,
459 /* Here we decompose the motion of a constrained atom into two
460 * components: rotation around the COM and translation of the COM.
463 /* Determine the variance of the arc length for the two rotational DOFs */
464 real massFraction
= prop
->con_mass
/(prop
->mass
+ prop
->con_mass
);
465 real sigma2_rot
= kT_fac
*massFraction
/prop
->mass
;
467 /* The distance from the atom to the COM, i.e. the rotational arm */
468 real comDistance
= prop
->con_len
*massFraction
;
470 /* The variance relative to the arm */
471 real sigma2_rel
= sigma2_rot
/gmx::square(comDistance
);
473 /* For sigma2_rel << 1 we don't notice the rotational effect and
474 * we have a normal, Gaussian displacement distribution.
475 * For larger sigma2_rel the displacement is much less, in fact it can
476 * not exceed 2*comDistance. We can calculate MSD/arm^2 as:
477 * integral_x=0-inf distance2(x) x/sigma2_rel exp(-x^2/(2 sigma2_rel)) dx
478 * where x is angular displacement and distance2(x) is the distance^2
479 * between points at angle 0 and x:
480 * distance2(x) = (sin(x) - sin(0))^2 + (cos(x) - cos(0))^2
481 * The limiting value of this MSD is 2, which is also the value for
482 * a uniform rotation distribution that would be reached at long time.
483 * The maximum is 2.5695 at sigma2_rel = 4.5119.
484 * We approximate this integral with a rational polynomial with
485 * coefficients from a Taylor expansion. This approximation is an
486 * overestimate for all values of sigma2_rel. Its maximum value
487 * of 2.6491 is reached at sigma2_rel = sqrt(45/2) = 4.7434.
488 * We keep the approximation constant after that.
489 * We use this approximate MSD as the variance for a Gaussian distribution.
491 * NOTE: For any sensible buffer tolerance this will result in a (large)
492 * overestimate of the buffer size, since the Gaussian has a long tail,
493 * whereas the actual distribution can not reach values larger than 2.
495 /* Coeffients obtained from a Taylor expansion */
496 const real a
= 1.0/3.0;
497 const real b
= 2.0/45.0;
499 /* Our approximation is constant after sigma2_rel = 1/sqrt(b) */
500 sigma2_rel
= std::min(sigma2_rel
, 1/std::sqrt(b
));
502 /* Compute the approximate sigma^2 for 2D motion due to the rotation */
503 *sigma2_2d
= gmx::square(comDistance
)*
504 sigma2_rel
/(1 + a
*sigma2_rel
+ b
*gmx::square(sigma2_rel
));
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
= std::exp(-x
*x
/(2*s2
));
541 er
= std::erfc(x
/std::sqrt(2*s2
));
543 *shift
= -x
+ std::sqrt(2*s2
/M_PI
)*ex
/er
;
544 *scale
= 0.5*M_PI
*std::exp(ex
*ex
/(M_PI
*er
*er
))*er
;
547 // Returns an (over)estimate of the energy drift for a single atom pair,
548 // given the kinetic properties, displacement variances and list buffer.
549 static real
energyDriftAtomPair(bool isConstrained_i
,
550 bool isConstrained_j
,
551 real s2
, real s2i_2d
, real s2j_2d
,
553 const pot_derivatives_t
*der
)
555 // For relatively small arguments erfc() is so small that if will be 0.0
556 // when stored in a float. We set an argument limit of 8 (Erfc(8)=1e-29),
557 // such that we can divide by erfc and have some space left for arithmetic.
558 const real erfc_arg_max
= 8.0;
565 if (rsh
*rsh
> 2*s2
*erfc_arg_max
*erfc_arg_max
)
567 // Below we calculate c_erfc = 0.5*erfc(rsh/sqrt(2*s2))
568 // When rsh/sqrt(2*s2) increases, this erfc will be the first
569 // result that underflows and becomes 0.0. To avoid this,
570 // we set c_exp=0 and c_erfc=0 for large arguments.
571 // This also avoids NaN in approx_2dof().
572 // In any relevant case this has no effect on the results,
573 // since c_exp < 6e-29, so the displacement is completely
574 // negligible for such atom pairs (and an overestimate).
575 // In nearly all use cases, there will be other atom pairs
576 // that contribute much more to the total, so zeroing
577 // this particular contribution has no effect at all.
583 /* For constraints: adapt r and scaling for the Gaussian */
588 approx_2dof(s2i_2d
, r_buffer
*s2i_2d
/s2
, &sh
, &sc
);
596 approx_2dof(s2j_2d
, r_buffer
*s2j_2d
/s2
, &sh
, &sc
);
601 /* Exact contribution of an atom pair with Gaussian displacement
602 * with sigma s to the energy drift for a potential with
603 * derivative -md and second derivative dd at the cut-off.
604 * The only catch is that for potentials that change sign
605 * near the cut-off there could be an unlucky compensation
606 * of positive and negative energy drift.
607 * Such potentials are extremely rare though.
609 * Note that pot has unit energy*length, as the linear
610 * atom density still needs to be put in.
612 c_exp
= std::exp(-rsh
*rsh
/(2*s2
))/std::sqrt(2*M_PI
);
613 c_erfc
= 0.5*std::erfc(rsh
/(std::sqrt(2*s2
)));
615 real s
= std::sqrt(s2
);
619 der
->md1
/2*((rsh2
+ s2
)*c_erfc
- rsh
*s
*c_exp
);
621 der
->d2
/6*(s
*(rsh2
+ 2*s2
)*c_exp
- rsh
*(rsh2
+ 3*s2
)*c_erfc
);
623 der
->md3
/24*((rsh2
*rsh2
+ 6*rsh2
*s2
+ 3*s2
*s2
)*c_erfc
- rsh
*s
*(rsh2
+ 5*s2
)*c_exp
);
625 return pot1
+ pot2
+ pot3
;
628 // Computes and returns an estimate of the energy drift for the whole system
629 static real
energyDrift(gmx::ArrayRef
<const VerletbufAtomtype
> att
,
630 const gmx_ffparams_t
*ffp
,
632 const pot_derivatives_t
*ljDisp
,
633 const pot_derivatives_t
*ljRep
,
634 const pot_derivatives_t
*elec
,
635 real rlj
, real rcoulomb
,
636 real rlist
, real boxvol
)
638 double drift_tot
= 0;
642 /* No atom displacements: no drift, avoid division by 0 */
646 // Here add up the contribution of all atom pairs in the system to
647 // (estimated) energy drift by looping over all atom type pairs.
648 for (gmx::index i
= 0; i
< att
.ssize(); i
++)
650 // Get the thermal displacement variance for the i-atom type
651 const atom_nonbonded_kinetic_prop_t
*prop_i
= &att
[i
].prop
;
653 get_atom_sigma2(kT_fac
, prop_i
, &s2i_2d
, &s2i_3d
);
655 for (gmx::index j
= i
; j
< att
.ssize(); j
++)
657 // Get the thermal displacement variance for the j-atom type
658 const atom_nonbonded_kinetic_prop_t
*prop_j
= &att
[j
].prop
;
660 get_atom_sigma2(kT_fac
, prop_j
, &s2j_2d
, &s2j_3d
);
662 /* Add up the up to four independent variances */
663 real s2
= s2i_2d
+ s2i_3d
+ s2j_2d
+ s2j_3d
;
665 // Set -V', V'' and -V''' at the cut-off for LJ */
666 real c6
= ffp
->iparams
[prop_i
->type
*ffp
->atnr
+ prop_j
->type
].lj
.c6
;
667 real c12
= ffp
->iparams
[prop_i
->type
*ffp
->atnr
+ prop_j
->type
].lj
.c12
;
668 pot_derivatives_t lj
;
669 lj
.md1
= c6
*ljDisp
->md1
+ c12
*ljRep
->md1
;
670 lj
.d2
= c6
*ljDisp
->d2
+ c12
*ljRep
->d2
;
671 lj
.md3
= c6
*ljDisp
->md3
+ c12
*ljRep
->md3
;
673 real pot_lj
= energyDriftAtomPair(prop_i
->bConstr
, prop_j
->bConstr
,
678 // Set -V' and V'' at the cut-off for Coulomb
679 pot_derivatives_t elec_qq
;
680 elec_qq
.md1
= elec
->md1
*prop_i
->q
*prop_j
->q
;
681 elec_qq
.d2
= elec
->d2
*prop_i
->q
*prop_j
->q
;
684 real pot_q
= energyDriftAtomPair(prop_i
->bConstr
, prop_j
->bConstr
,
689 // Note that attractive and repulsive potentials for individual
690 // pairs can partially cancel.
691 real pot
= pot_lj
+ pot_q
;
693 /* Multiply by the number of atom pairs */
696 pot
*= static_cast<double>(att
[i
].n
)*(att
[i
].n
- 1)/2;
700 pot
*= static_cast<double>(att
[i
].n
)*att
[j
].n
;
702 /* We need the line density to get the energy drift of the system.
703 * The effective average r^2 is close to (rlist+sigma)^2.
705 pot
*= 4*M_PI
*gmx::square(rlist
+ std::sqrt(s2
))/boxvol
;
707 /* Add the unsigned drift to avoid cancellation of errors */
708 drift_tot
+= std::abs(pot
);
715 // Returns the chance that a particle in a cluster is at distance rlist
716 // when the cluster is at distance rlist
717 static real
surface_frac(int cluster_size
, real particle_distance
, real rlist
)
721 if (rlist
< 0.5*particle_distance
)
723 /* We have non overlapping spheres */
727 /* Half the inter-particle distance relative to rlist */
728 d
= 0.5*particle_distance
/rlist
;
730 /* Determine the area of the surface at distance rlist to the closest
731 * particle, relative to surface of a sphere of radius rlist.
732 * The formulas below assume close to cubic cells for the pair search grid,
733 * which the pair search code tries to achieve.
734 * Note that in practice particle distances will not be delta distributed,
735 * but have some spread, often involving shorter distances,
736 * as e.g. O-H bonds in a water molecule. Thus the estimates below will
737 * usually be slightly too high and thus conservative.
739 switch (cluster_size
)
742 /* One particle: trivial */
746 /* Two particles: two spheres at fractional distance 2*a */
750 /* We assume a perfect, symmetric tetrahedron geometry.
751 * The surface around a tetrahedron is too complex for a full
752 * analytical solution, so we use a Taylor expansion.
754 area_rel
= (1.0 + 1/M_PI
*(6*std::acos(1/std::sqrt(3))*d
+
755 std::sqrt(3)*d
*d
*(1.0 +
758 83.0/756.0*d
*d
*d
*d
*d
*d
)));
761 gmx_incons("surface_frac called with unsupported cluster_size");
764 return area_rel
/cluster_size
;
767 /* Returns the negative of the third derivative of a potential r^-p
768 * with a force-switch function, evaluated at the cut-off rc.
770 static real
md3_force_switch(real p
, real rswitch
, real rc
)
772 /* The switched force function is:
773 * p*r^-(p+1) + a*(r - rswitch)^2 + b*(r - rswitch)^3
776 real md3_pot
, md3_sw
;
778 a
= -((p
+ 4)*rc
- (p
+ 1)*rswitch
)/(pow(rc
, p
+2)*gmx::square(rc
-rswitch
));
779 b
= ((p
+ 3)*rc
- (p
+ 1)*rswitch
)/(pow(rc
, p
+2)*gmx::power3(rc
-rswitch
));
781 md3_pot
= (p
+ 2)*(p
+ 1)*p
*pow(rc
, p
+3);
782 md3_sw
= 2*a
+ 6*b
*(rc
- rswitch
);
784 return md3_pot
+ md3_sw
;
787 /* Returns the variance of the atomic displacement over timePeriod.
789 * Note: When not using BD with a non-mass dependendent friction coefficient,
790 * the return value still needs to be divided by the particle mass.
792 static real
displacementVariance(const t_inputrec
&ir
,
800 /* Get the displacement distribution from the random component only.
801 * With accurate integration the systematic (force) displacement
802 * should be negligible (unless nstlist is extremely large, which
803 * you wouldn't do anyhow).
805 kT_fac
= 2*BOLTZ
*temperature
*timePeriod
;
808 /* This is directly sigma^2 of the displacement */
809 kT_fac
/= ir
.bd_fric
;
813 /* Per group tau_t is not implemented yet, use the maximum */
814 real tau_t
= ir
.opts
.tau_t
[0];
815 for (int i
= 1; i
< ir
.opts
.ngtc
; i
++)
817 tau_t
= std::max(tau_t
, ir
.opts
.tau_t
[i
]);
821 /* This kT_fac needs to be divided by the mass to get sigma^2 */
826 kT_fac
= BOLTZ
*temperature
*gmx::square(timePeriod
);
832 /* Returns the largest sigma of the Gaussian displacement over all particle
833 * types. This ignores constraints, so is an overestimate.
835 static real
maxSigma(real kT_fac
,
836 gmx::ArrayRef
<const VerletbufAtomtype
> att
)
838 GMX_ASSERT(!att
.empty(), "We should have at least one type");
839 real smallestMass
= att
[0].prop
.mass
;
840 for (const auto &atomType
: att
)
842 smallestMass
= std::min(smallestMass
, atomType
.prop
.mass
);
846 return 2*std::sqrt(kT_fac
/smallestMass
);
850 calcVerletBufferSize(const gmx_mtop_t
&mtop
,
851 const real boxVolume
,
852 const t_inputrec
&ir
,
854 const int listLifetime
,
855 real referenceTemperature
,
856 const VerletbufListSetup
&listSetup
)
861 real particle_distance
;
862 real nb_clust_frac_pairs_not_in_list_at_cutoff
;
869 if (!EI_DYNAMICS(ir
.eI
))
871 gmx_incons("Can only determine the Verlet buffer size for integrators that perform dynamics");
873 if (ir
.verletbuf_tol
<= 0)
875 gmx_incons("The Verlet buffer tolerance needs to be larger than zero");
878 if (referenceTemperature
< 0)
880 /* We use the maximum temperature with multiple T-coupl groups.
881 * We could use a per particle temperature, but since particles
882 * interact, this might underestimate the buffer size.
884 referenceTemperature
= maxReferenceTemperature(ir
);
886 GMX_RELEASE_ASSERT(referenceTemperature
>= 0, "Without T-coupling we should not end up here");
889 /* Resolution of the buffer size */
892 env
= getenv("GMX_VERLET_BUFFER_RES");
895 sscanf(env
, "%lf", &resolution
);
898 /* In an atom wise pair-list there would be no pairs in the list
899 * beyond the pair-list cut-off.
900 * However, we use a pair-list of groups vs groups of atoms.
901 * For groups of 4 atoms, the parallelism of SSE instructions, only
902 * 10% of the atoms pairs are not in the list just beyond the cut-off.
903 * As this percentage increases slowly compared to the decrease of the
904 * Gaussian displacement distribution over this range, we can simply
905 * reduce the drift by this fraction.
906 * For larger groups, e.g. of 8 atoms, this fraction will be lower,
907 * so then buffer size will be on the conservative (large) side.
909 * Note that the formulas used here do not take into account
910 * cancellation of errors which could occur by missing both
911 * attractive and repulsive interactions.
913 * The only major assumption is homogeneous particle distribution.
914 * For an inhomogeneous system, such as a liquid-vapor system,
915 * the buffer will be underestimated. The actual energy drift
916 * will be higher by the factor: local/homogeneous particle density.
918 * The results of this estimate have been checked againt simulations.
919 * In most cases the real drift differs by less than a factor 2.
922 /* Worst case assumption: HCP packing of particles gives largest distance */
923 particle_distance
= std::cbrt(boxVolume
*std::sqrt(2)/mtop
.natoms
);
925 /* TODO: Obtain masses through (future) integrator functionality
926 * to avoid scattering the code with (or forgetting) checks.
928 const bool setMassesToOne
= (ir
.eI
== eiBD
&& ir
.bd_fric
> 0);
930 getVerletBufferAtomtypes(mtop
, setMassesToOne
);
931 GMX_ASSERT(!att
.empty(), "We expect at least one type");
935 fprintf(debug
, "particle distance assuming HCP packing: %f nm\n",
937 fprintf(debug
, "energy drift atom types: %zu\n", att
.size());
940 pot_derivatives_t ljDisp
= { 0, 0, 0 };
941 pot_derivatives_t ljRep
= { 0, 0, 0 };
942 real repPow
= mtop
.ffparams
.reppow
;
944 if (ir
.vdwtype
== evdwCUT
)
946 real sw_range
, md3_pswf
;
948 switch (ir
.vdw_modifier
)
951 case eintmodPOTSHIFT
:
952 /* -dV/dr of -r^-6 and r^-reppow */
953 ljDisp
.md1
= -6*std::pow(ir
.rvdw
, -7.0);
954 ljRep
.md1
= repPow
*std::pow(ir
.rvdw
, -(repPow
+ 1));
955 /* The contribution of the higher derivatives is negligible */
957 case eintmodFORCESWITCH
:
958 /* At the cut-off: V=V'=V''=0, so we use only V''' */
959 ljDisp
.md3
= -md3_force_switch(6.0, ir
.rvdw_switch
, ir
.rvdw
);
960 ljRep
.md3
= md3_force_switch(repPow
, ir
.rvdw_switch
, ir
.rvdw
);
962 case eintmodPOTSWITCH
:
963 /* At the cut-off: V=V'=V''=0.
964 * V''' is given by the original potential times
965 * the third derivative of the switch function.
967 sw_range
= ir
.rvdw
- ir
.rvdw_switch
;
968 md3_pswf
= 60.0/gmx::power3(sw_range
);
970 ljDisp
.md3
= -std::pow(ir
.rvdw
, -6.0 )*md3_pswf
;
971 ljRep
.md3
= std::pow(ir
.rvdw
, -repPow
)*md3_pswf
;
974 gmx_incons("Unimplemented VdW modifier");
977 else if (EVDW_PME(ir
.vdwtype
))
979 real b
= calc_ewaldcoeff_lj(ir
.rvdw
, ir
.ewald_rtol_lj
);
985 // -dV/dr of g(br)*r^-6 [where g(x) = exp(-x^2)(1+x^2+x^4/2),
986 // see LJ-PME equations in manual] and r^-reppow
987 ljDisp
.md1
= -std::exp(-br2
)*(br6
+ 3.0*br4
+ 6.0*br2
+ 6.0)*std::pow(r
, -7.0);
988 ljRep
.md1
= repPow
*pow(r
, -(repPow
+ 1));
989 // The contribution of the higher derivatives is negligible
993 gmx_fatal(FARGS
, "Energy drift calculation is only implemented for plain cut-off Lennard-Jones interactions");
996 elfac
= ONE_4PI_EPS0
/ir
.epsilon_r
;
998 // Determine the 1st and 2nd derivative for the electostatics
999 pot_derivatives_t elec
= { 0, 0, 0 };
1001 if (ir
.coulombtype
== eelCUT
|| EEL_RF(ir
.coulombtype
))
1005 if (ir
.coulombtype
== eelCUT
)
1012 eps_rf
= ir
.epsilon_rf
/ir
.epsilon_r
;
1015 k_rf
= (eps_rf
- ir
.epsilon_r
)/( gmx::power3(ir
.rcoulomb
) * (2*eps_rf
+ ir
.epsilon_r
) );
1019 /* epsilon_rf = infinity */
1020 k_rf
= 0.5/gmx::power3(ir
.rcoulomb
);
1026 elec
.md1
= elfac
*(1.0/gmx::square(ir
.rcoulomb
) - 2*k_rf
*ir
.rcoulomb
);
1028 elec
.d2
= elfac
*(2.0/gmx::power3(ir
.rcoulomb
) + 2*k_rf
);
1030 else if (EEL_PME(ir
.coulombtype
) || ir
.coulombtype
== eelEWALD
)
1034 b
= calc_ewaldcoeff_q(ir
.rcoulomb
, ir
.ewald_rtol
);
1037 elec
.md1
= elfac
*(b
*std::exp(-br
*br
)*M_2_SQRTPI
/rc
+ std::erfc(br
)/(rc
*rc
));
1038 elec
.d2
= elfac
/(rc
*rc
)*(2*b
*(1 + br
*br
)*std::exp(-br
*br
)*M_2_SQRTPI
+ 2*std::erfc(br
)/rc
);
1042 gmx_fatal(FARGS
, "Energy drift calculation is only implemented for Reaction-Field and Ewald electrostatics");
1045 /* Determine the variance of the atomic displacement
1046 * over list_lifetime steps: kT_fac
1047 * For inertial dynamics (not Brownian dynamics) the mass factor
1048 * is not included in kT_fac, it is added later.
1050 const real kT_fac
= displacementVariance(ir
, referenceTemperature
,
1051 listLifetime
*ir
.delta_t
);
1055 fprintf(debug
, "Derivatives of non-bonded potentials at the cut-off:\n");
1056 fprintf(debug
, "LJ disp. -V' %9.2e V'' %9.2e -V''' %9.2e\n", ljDisp
.md1
, ljDisp
.d2
, ljDisp
.md3
);
1057 fprintf(debug
, "LJ rep. -V' %9.2e V'' %9.2e -V''' %9.2e\n", ljRep
.md1
, ljRep
.d2
, ljRep
.md3
);
1058 fprintf(debug
, "Electro. -V' %9.2e V'' %9.2e\n", elec
.md1
, elec
.d2
);
1059 fprintf(debug
, "sqrt(kT_fac) %f\n", std::sqrt(kT_fac
));
1062 /* Search using bisection */
1064 /* The drift will be neglible at 5 times the max sigma */
1065 ib1
= static_cast<int>(5*maxSigma(kT_fac
, att
)/resolution
) + 1;
1066 while (ib1
- ib0
> 1)
1070 rl
= std::max(ir
.rvdw
, ir
.rcoulomb
) + rb
;
1072 /* Calculate the average energy drift at the last step
1073 * of the nstlist steps at which the pair-list is used.
1075 drift
= energyDrift(att
, &mtop
.ffparams
,
1077 &ljDisp
, &ljRep
, &elec
,
1078 ir
.rvdw
, ir
.rcoulomb
,
1081 /* Correct for the fact that we are using a Ni x Nj particle pair list
1082 * and not a 1 x 1 particle pair list. This reduces the drift.
1084 /* We don't have a formula for 8 (yet), use 4 which is conservative */
1085 nb_clust_frac_pairs_not_in_list_at_cutoff
=
1086 surface_frac(std::min(listSetup
.cluster_size_i
, 4),
1087 particle_distance
, rl
)*
1088 surface_frac(std::min(listSetup
.cluster_size_j
, 4),
1089 particle_distance
, rl
);
1090 drift
*= nb_clust_frac_pairs_not_in_list_at_cutoff
;
1092 /* Convert the drift to drift per unit time per atom */
1093 drift
/= nstlist
*ir
.delta_t
*mtop
.natoms
;
1097 fprintf(debug
, "ib %3d %3d %3d rb %.3f %dx%d fac %.3f drift %.1e\n",
1099 listSetup
.cluster_size_i
, listSetup
.cluster_size_j
,
1100 nb_clust_frac_pairs_not_in_list_at_cutoff
,
1104 if (std::abs(drift
) > ir
.verletbuf_tol
)
1114 return std::max(ir
.rvdw
, ir
.rcoulomb
) + ib1
*resolution
;
1117 /* Returns the pairlist buffer size for use as a minimum buffer size
1119 * Note that this is a rather crude estimate. It is ok for a buffer
1120 * set for good energy conservation or RF electrostatics. But it is
1121 * too small with PME and the buffer set with the default tolerance.
1123 static real
minCellSizeFromPairlistBuffer(const t_inputrec
&ir
)
1125 return ir
.rlist
- std::max(ir
.rvdw
, ir
.rcoulomb
);
1129 chanceOfAtomCrossingCell(gmx::ArrayRef
<const VerletbufAtomtype
> atomtypes
,
1133 /* We assume atoms are distributed uniformly over the cell width.
1134 * Once an atom has moved by more than the cellSize (as passed
1135 * as the buffer argument to energyDriftAtomPair() below),
1136 * the chance of crossing the boundary of the neighbor cell
1137 * thus increases as 1/cellSize with the additional displacement
1138 * on top of cellSize. We thus create a linear interaction with
1139 * derivative = -1/cellSize. Using this in the energyDriftAtomPair
1140 * function will return the chance of crossing the next boundary.
1142 const pot_derivatives_t boundaryInteraction
= { 1/cellSize
, 0, 0 };
1145 for (const VerletbufAtomtype
&att
: atomtypes
)
1147 const atom_nonbonded_kinetic_prop_t
&propAtom
= att
.prop
;
1150 get_atom_sigma2(kT_fac
, &propAtom
, &s2_2d
, &s2_3d
);
1152 real chancePerAtom
= energyDriftAtomPair(propAtom
.bConstr
, false,
1153 s2_2d
+ s2_3d
, s2_2d
, 0,
1155 &boundaryInteraction
);
1157 if (propAtom
.bConstr
)
1159 /* energyDriftAtomPair() uses an unlimited Gaussian displacement
1160 * distribution for constrained atoms, whereas they can
1161 * actually not move more than the COM of the two constrained
1162 * atoms plus twice the distance from the COM.
1163 * Use this maximum, limited displacement when this results in
1164 * a smaller chance (note that this is still an overestimate).
1166 real massFraction
= propAtom
.con_mass
/(propAtom
.mass
+ propAtom
.con_mass
);
1167 real comDistance
= propAtom
.con_len
*massFraction
;
1169 real chanceWithMaxDistance
=
1170 energyDriftAtomPair(false, false,
1172 cellSize
- 2*comDistance
,
1173 &boundaryInteraction
);
1174 chancePerAtom
= std::min(chancePerAtom
, chanceWithMaxDistance
);
1177 /* Take into account the line density of the boundary */
1178 chancePerAtom
/= cellSize
;
1180 chance
+= att
.n
*chancePerAtom
;
1186 /* Struct for storing constraint properties of atoms */
1187 struct AtomConstraintProps
1189 void addConstraint(real length
)
1191 numConstraints
+= 1;
1192 sumLengths
+= length
;
1195 int numConstraints
= 0; /* The number of constraints of an atom */
1196 real sumLengths
= 0; /* The sum of constraint length over the constraints */
1199 /* Constructs and returns a list of constraint properties per atom */
1200 static std::vector
<AtomConstraintProps
>
1201 getAtomConstraintProps(const gmx_moltype_t
&moltype
,
1202 const gmx_ffparams_t
&ffparams
)
1204 const t_atoms
&atoms
= moltype
.atoms
;
1205 std::vector
<AtomConstraintProps
> props(atoms
.nr
);
1207 for (const auto &ilist
: extractILists(moltype
.ilist
, IF_CONSTRAINT
))
1209 // Settles are handled separately
1210 if (ilist
.functionType
== F_SETTLE
)
1215 for (size_t i
= 0; i
< ilist
.iatoms
.size(); i
+= ilistStride(ilist
))
1217 int type
= ilist
.iatoms
[i
];
1218 int a1
= ilist
.iatoms
[i
+ 1];
1219 int a2
= ilist
.iatoms
[i
+ 2];
1220 real length
= ffparams
.iparams
[type
].constr
.dA
;
1221 props
[a1
].addConstraint(length
);
1222 props
[a2
].addConstraint(length
);
1229 /* Return the chance of at least one update group in a molecule crossing a cell of size cellSize */
1231 chanceOfUpdateGroupCrossingCell(const gmx_moltype_t
&moltype
,
1232 const gmx_ffparams_t
&ffparams
,
1233 const gmx::RangePartitioning
&updateGrouping
,
1237 const t_atoms
&atoms
= moltype
.atoms
;
1238 GMX_ASSERT(updateGrouping
.fullRange().end() == atoms
.nr
, "The update groups should match the molecule type");
1240 const pot_derivatives_t boundaryInteraction
= { 1/cellSize
, 0, 0 };
1242 const auto atomConstraintProps
= getAtomConstraintProps(moltype
, ffparams
);
1245 for (int group
= 0; group
< updateGrouping
.numBlocks(); group
++)
1247 const auto &block
= updateGrouping
.block(group
);
1248 /* Determine the number of atoms with constraints and the mass of the COG */
1249 int numAtomsWithConstraints
= 0;
1251 for (const int atom
: block
)
1253 if (atomConstraintProps
[atom
].numConstraints
> 0)
1255 numAtomsWithConstraints
++;
1257 massSum
+= moltype
.atoms
.atom
[atom
].m
;
1259 /* Determine the maximum possible distance between the center of mass
1260 * and the center of geometry of the update group
1262 real maxComCogDistance
= 0;
1263 if (numAtomsWithConstraints
== 2)
1265 for (const int atom
: block
)
1267 if (atomConstraintProps
[atom
].numConstraints
> 0)
1269 GMX_ASSERT(atomConstraintProps
[atom
].numConstraints
== 1,
1270 "Two atoms should be connected by one constraint");
1271 maxComCogDistance
= std::abs(atoms
.atom
[atom
].m
/massSum
- 0.5)*atomConstraintProps
[atom
].sumLengths
;
1276 else if (numAtomsWithConstraints
> 2)
1278 for (const int atom
: block
)
1280 if (atomConstraintProps
[atom
].numConstraints
== numAtomsWithConstraints
- 1)
1282 real comCogDistance
= atomConstraintProps
[atom
].sumLengths
/numAtomsWithConstraints
;
1283 maxComCogDistance
= std::max(maxComCogDistance
, comCogDistance
);
1287 else if (block
.size() > 1)
1289 // All normal atoms must be connected by SETTLE
1290 for (const int atom
: block
)
1292 const auto &ilist
= moltype
.ilist
[F_SETTLE
];
1293 GMX_RELEASE_ASSERT(ilist
.size() > 0, "There should be at least one settle in this moltype");
1294 for (int i
= 0; i
< ilist
.size(); i
+= 1 + NRAL(F_SETTLE
))
1296 if (atom
== ilist
.iatoms
[i
+ 1])
1298 const t_iparams
&iparams
= ffparams
.iparams
[ilist
.iatoms
[i
]];
1299 real dOH
= iparams
.settle
.doh
;
1300 real dHH
= iparams
.settle
.dhh
;
1301 real dOMidH
= std::sqrt(dOH
*dOH
- 0.25_real
*dHH
*dHH
);
1302 maxComCogDistance
= std::abs(atoms
.atom
[atom
].m
/massSum
- 1.0_real
/3.0_real
)*dOMidH
;
1307 real s2_3d
= kT_fac
/massSum
;
1308 chance
+= energyDriftAtomPair(false, false,
1310 cellSize
- 2*maxComCogDistance
,
1311 &boundaryInteraction
);
1317 /* Return the chance of at least one update group in the system crossing a cell of size cellSize */
1319 chanceOfUpdateGroupCrossingCell(const gmx_mtop_t
&mtop
,
1320 PartitioningPerMoltype updateGrouping
,
1324 GMX_RELEASE_ASSERT(updateGrouping
.size() == mtop
.moltype
.size(),
1325 "The update groups should match the topology");
1328 for (const gmx_molblock_t
&molblock
: mtop
.molblock
)
1330 const gmx_moltype_t
&moltype
= mtop
.moltype
[molblock
.type
];
1332 molblock
.nmol
*chanceOfUpdateGroupCrossingCell(moltype
, mtop
.ffparams
, updateGrouping
[molblock
.type
],
1340 minCellSizeForAtomDisplacement(const gmx_mtop_t
&mtop
,
1341 const t_inputrec
&ir
,
1342 PartitioningPerMoltype updateGrouping
,
1343 real chanceRequested
)
1345 if (!EI_DYNAMICS(ir
.eI
) || (EI_MD(ir
.eI
) && ir
.etc
== etcNO
))
1347 return minCellSizeFromPairlistBuffer(ir
);
1350 /* We use the maximum temperature with multiple T-coupl groups.
1351 * We could use a per particle temperature, but since particles
1352 * interact, this might underestimate the displacements.
1354 const real temperature
= maxReferenceTemperature(ir
);
1356 const bool setMassesToOne
= (ir
.eI
== eiBD
&& ir
.bd_fric
> 0);
1358 const auto atomtypes
= getVerletBufferAtomtypes(mtop
, setMassesToOne
);
1360 const real kT_fac
= displacementVariance(ir
, temperature
,
1361 ir
.nstlist
*ir
.delta_t
);
1363 /* Resolution of the cell size */
1364 real resolution
= 0.001;
1366 /* Search using bisection, avoid 0 and start at 1 */
1368 /* The chance will be neglible at 10 times the max sigma */
1369 int ib1
= int(10*maxSigma(kT_fac
, atomtypes
)/resolution
) + 1;
1371 while (ib1
- ib0
> 1)
1373 int ib
= (ib0
+ ib1
)/2;
1374 cellSize
= ib
*resolution
;
1377 if (updateGrouping
.empty())
1379 chance
= chanceOfAtomCrossingCell(atomtypes
, kT_fac
, cellSize
);
1383 chance
= chanceOfUpdateGroupCrossingCell(mtop
, updateGrouping
, kT_fac
, cellSize
);
1386 /* Note: chance is for every nstlist steps */
1387 if (chance
> chanceRequested
*ir
.nstlist
)