2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2012,2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
37 #include "calc_verletbuf.h"
46 #include "gromacs/ewald/ewald-utils.h"
47 #include "gromacs/math/functions.h"
48 #include "gromacs/math/units.h"
49 #include "gromacs/math/vec.h"
50 #include "gromacs/mdlib/nb_verlet.h"
51 #include "gromacs/mdlib/nbnxn_simd.h"
52 #include "gromacs/mdlib/nbnxn_util.h"
53 #include "gromacs/mdtypes/inputrec.h"
54 #include "gromacs/mdtypes/md_enums.h"
55 #include "gromacs/topology/ifunc.h"
56 #include "gromacs/topology/topology.h"
57 #include "gromacs/utility/fatalerror.h"
58 #include "gromacs/utility/smalloc.h"
60 /* The code in this file estimates a pairlist buffer length
61 * given a target energy drift per atom per picosecond.
62 * This is done by estimating the drift given a buffer length.
63 * Ideally we would like to have a tight overestimate of the drift,
64 * but that can be difficult to achieve.
66 * Significant approximations used:
68 * Uniform particle density. UNDERESTIMATES the drift by rho_global/rho_local.
70 * Interactions don't affect particle motion. OVERESTIMATES the drift on longer
71 * time scales. This approximation probably introduces the largest errors.
73 * Only take one constraint per particle into account: OVERESTIMATES the drift.
75 * For rotating constraints assume the same functional shape for time scales
76 * where the constraints rotate significantly as the exact expression for
77 * short time scales. OVERESTIMATES the drift on long time scales.
79 * For non-linear virtual sites use the mass of the lightest constructing atom
80 * to determine the displacement. OVER/UNDERESTIMATES the drift, depending on
81 * the geometry and masses of constructing atoms.
83 * Note that the formulas for normal atoms and linear virtual sites are exact,
84 * apart from the first two approximations.
86 * Note that apart from the effect of the above approximations, the actual
87 * drift of the total energy of a system can be orders of magnitude smaller
88 * due to cancellation of positive and negative drift for different pairs.
92 /* Struct for unique atom type for calculating the energy drift.
93 * The atom displacement depends on mass and constraints.
94 * The energy jump for given distance depend on LJ type and q.
98 atom_nonbonded_kinetic_prop_t prop
; /* non-bonded and kinetic atom prop. */
99 int n
; /* #atoms of this type in the system */
100 } verletbuf_atomtype_t
;
102 // Struct for derivatives of a non-bonded interaction potential
105 real md1
; // -V' at the cutoff
106 real d2
; // V'' at the cutoff
107 real md3
; // -V''' at the cutoff
110 VerletbufListSetup
verletbufGetListSetup(int nbnxnKernelType
)
112 /* Note that the current buffer estimation code only handles clusters
113 * of size 1, 2 or 4, so for 4x8 or 8x8 we use the estimate for 4x4.
115 VerletbufListSetup listSetup
;
117 listSetup
.cluster_size_i
= nbnxn_kernel_to_cluster_i_size(nbnxnKernelType
);
118 listSetup
.cluster_size_j
= nbnxn_kernel_to_cluster_j_size(nbnxnKernelType
);
120 if (nbnxnKernelType
== nbnxnk8x8x8_GPU
||
121 nbnxnKernelType
== nbnxnk8x8x8_PlainC
)
123 /* The GPU kernels (except for OpenCL) split the j-clusters in two halves */
124 listSetup
.cluster_size_j
/= 2;
130 VerletbufListSetup
verletbufGetSafeListSetup(ListSetupType listType
)
132 /* When calling this function we often don't know which kernel type we
133 * are going to use. We choose the kernel type with the smallest possible
134 * i- and j-cluster sizes, so we potentially overestimate, but never
135 * underestimate, the buffer drift.
139 if (listType
== ListSetupType::Gpu
)
141 nbnxnKernelType
= nbnxnk8x8x8_GPU
;
143 else if (GMX_SIMD
&& listType
== ListSetupType::CpuSimdWhenSupported
)
145 #ifdef GMX_NBNXN_SIMD_2XNN
146 /* We use the smallest cluster size to be on the safe side */
147 nbnxnKernelType
= nbnxnk4xN_SIMD_2xNN
;
149 nbnxnKernelType
= nbnxnk4xN_SIMD_4xN
;
154 nbnxnKernelType
= nbnxnk4x4_PlainC
;
157 return verletbufGetListSetup(nbnxnKernelType
);
161 atom_nonbonded_kinetic_prop_equal(const atom_nonbonded_kinetic_prop_t
*prop1
,
162 const atom_nonbonded_kinetic_prop_t
*prop2
)
164 return (prop1
->mass
== prop2
->mass
&&
165 prop1
->type
== prop2
->type
&&
166 prop1
->q
== prop2
->q
&&
167 prop1
->bConstr
== prop2
->bConstr
&&
168 prop1
->con_mass
== prop2
->con_mass
&&
169 prop1
->con_len
== prop2
->con_len
);
172 static void add_at(verletbuf_atomtype_t
**att_p
, int *natt_p
,
173 const atom_nonbonded_kinetic_prop_t
*prop
,
176 verletbuf_atomtype_t
*att
;
181 /* Ignore massless particles */
189 while (i
< natt
&& !atom_nonbonded_kinetic_prop_equal(prop
, &att
[i
].prop
))
201 srenew(*att_p
, *natt_p
);
202 (*att_p
)[i
].prop
= *prop
;
203 (*att_p
)[i
].n
= nmol
;
207 static void get_vsite_masses(const gmx_moltype_t
*moltype
,
208 const gmx_ffparams_t
*ffparams
,
217 /* Check for virtual sites, determine mass from constructing atoms */
218 for (ft
= 0; ft
< F_NRE
; ft
++)
222 il
= &moltype
->ilist
[ft
];
224 for (i
= 0; i
< il
->nr
; i
+= 1+NRAL(ft
))
227 real inv_mass
, coeff
, m_aj
;
230 ip
= &ffparams
->iparams
[il
->iatoms
[i
]];
232 a1
= il
->iatoms
[i
+1];
236 /* Only vsiten can have more than four
237 constructing atoms, so NRAL(ft) <= 5 */
240 const int maxj
= NRAL(ft
);
244 for (j
= 1; j
< maxj
; j
++)
246 cam
[j
] = moltype
->atoms
.atom
[il
->iatoms
[i
+1+j
]].m
;
249 cam
[j
] = vsite_m
[il
->iatoms
[i
+1+j
]];
253 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.",
255 interaction_function
[ft
].longname
,
256 il
->iatoms
[i
+1+j
]+1);
264 vsite_m
[a1
] = (cam
[1]*cam
[2])/(cam
[2]*gmx::square(1-ip
->vsite
.a
) + cam
[1]*gmx::square(ip
->vsite
.a
));
268 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
));
271 gmx_incons("Invalid vsite type");
274 /* Use the mass of the lightest constructing atom.
275 * This is an approximation.
276 * If the distance of the virtual site to the
277 * constructing atom is less than all distances
278 * between constructing atoms, this is a safe
279 * over-estimate of the displacement of the vsite.
280 * This condition holds for all H mass replacement
281 * vsite constructions, except for SP2/3 groups.
282 * In SP3 groups one H will have a F_VSITE3
283 * construction, so even there the total drift
284 * estimate shouldn't be far off.
286 vsite_m
[a1
] = cam
[1];
287 for (j
= 2; j
< maxj
; j
++)
289 vsite_m
[a1
] = std::min(vsite_m
[a1
], cam
[j
]);
302 for (j
= 0; j
< 3*ffparams
->iparams
[il
->iatoms
[i
]].vsiten
.n
; j
+= 3)
304 aj
= il
->iatoms
[i
+j
+2];
305 coeff
= ffparams
->iparams
[il
->iatoms
[i
+j
]].vsiten
.a
;
306 if (moltype
->atoms
.atom
[aj
].ptype
== eptVSite
)
312 m_aj
= moltype
->atoms
.atom
[aj
].m
;
316 gmx_incons("The mass of a vsiten constructing atom is <= 0");
318 inv_mass
+= coeff
*coeff
/m_aj
;
320 vsite_m
[a1
] = 1/inv_mass
;
321 /* Correct for loop increment of i */
322 i
+= j
- 1 - NRAL(ft
);
326 fprintf(debug
, "atom %4d %-20s mass %6.3f\n",
327 a1
, interaction_function
[ft
].longname
, vsite_m
[a1
]);
334 static void get_verlet_buffer_atomtypes(const gmx_mtop_t
*mtop
,
335 verletbuf_atomtype_t
**att_p
,
339 verletbuf_atomtype_t
*att
;
341 int ft
, i
, a1
, a2
, a3
, a
;
344 atom_nonbonded_kinetic_prop_t
*prop
;
346 int n_nonlin_vsite_mol
;
351 if (n_nonlin_vsite
!= nullptr)
356 for (const gmx_molblock_t
&molblock
: mtop
->molblock
)
358 int nmol
= molblock
.nmol
;
359 const gmx_moltype_t
&moltype
= mtop
->moltype
[molblock
.type
];
360 const t_atoms
*atoms
= &moltype
.atoms
;
362 /* Check for constraints, as they affect the kinetic energy.
363 * For virtual sites we need the masses and geometry of
364 * the constructing atoms to determine their velocity distribution.
366 snew(prop
, atoms
->nr
);
367 snew(vsite_m
, atoms
->nr
);
369 for (ft
= F_CONSTR
; ft
<= F_CONSTRNC
; ft
++)
371 il
= &moltype
.ilist
[ft
];
373 for (i
= 0; i
< il
->nr
; i
+= 1+NRAL(ft
))
375 ip
= &mtop
->ffparams
.iparams
[il
->iatoms
[i
]];
376 a1
= il
->iatoms
[i
+1];
377 a2
= il
->iatoms
[i
+2];
378 if (atoms
->atom
[a2
].m
> prop
[a1
].con_mass
)
380 prop
[a1
].con_mass
= atoms
->atom
[a2
].m
;
381 prop
[a1
].con_len
= ip
->constr
.dA
;
383 if (atoms
->atom
[a1
].m
> prop
[a2
].con_mass
)
385 prop
[a2
].con_mass
= atoms
->atom
[a1
].m
;
386 prop
[a2
].con_len
= ip
->constr
.dA
;
391 il
= &moltype
.ilist
[F_SETTLE
];
393 for (i
= 0; i
< il
->nr
; i
+= 1+NRAL(F_SETTLE
))
395 ip
= &mtop
->ffparams
.iparams
[il
->iatoms
[i
]];
396 a1
= il
->iatoms
[i
+1];
397 a2
= il
->iatoms
[i
+2];
398 a3
= il
->iatoms
[i
+3];
399 /* Usually the mass of a1 (usually oxygen) is larger than a2/a3.
400 * If this is not the case, we overestimate the displacement,
401 * which leads to a larger buffer (ok since this is an exotic case).
403 prop
[a1
].con_mass
= atoms
->atom
[a2
].m
;
404 prop
[a1
].con_len
= ip
->settle
.doh
;
406 prop
[a2
].con_mass
= atoms
->atom
[a1
].m
;
407 prop
[a2
].con_len
= ip
->settle
.doh
;
409 prop
[a3
].con_mass
= atoms
->atom
[a1
].m
;
410 prop
[a3
].con_len
= ip
->settle
.doh
;
413 get_vsite_masses(&moltype
,
416 &n_nonlin_vsite_mol
);
417 if (n_nonlin_vsite
!= nullptr)
419 *n_nonlin_vsite
+= nmol
*n_nonlin_vsite_mol
;
422 for (a
= 0; a
< atoms
->nr
; a
++)
424 if (atoms
->atom
[a
].ptype
== eptVSite
)
426 prop
[a
].mass
= vsite_m
[a
];
430 prop
[a
].mass
= atoms
->atom
[a
].m
;
432 prop
[a
].type
= atoms
->atom
[a
].type
;
433 prop
[a
].q
= atoms
->atom
[a
].q
;
434 /* We consider an atom constrained, #DOF=2, when it is
435 * connected with constraints to (at least one) atom with
436 * a mass of more than 0.4x its own mass. This is not a critical
437 * parameter, since with roughly equal masses the unconstrained
438 * and constrained displacement will not differ much (and both
439 * overestimate the displacement).
441 prop
[a
].bConstr
= (prop
[a
].con_mass
> 0.4*prop
[a
].mass
);
443 add_at(&att
, &natt
, &prop
[a
], nmol
);
452 for (a
= 0; a
< natt
; a
++)
454 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",
455 a
, att
[a
].prop
.mass
, att
[a
].prop
.type
, att
[a
].prop
.q
,
456 att
[a
].prop
.bConstr
, att
[a
].prop
.con_mass
, att
[a
].prop
.con_len
,
465 /* This function computes two components of the estimate of the variance
466 * in the displacement of one atom in a system of two constrained atoms.
467 * Returns in sigma2_2d the variance due to rotation of the constrained
468 * atom around the atom to which it constrained.
469 * Returns in sigma2_3d the variance due to displacement of the COM
470 * of the whole system of the two constrained atoms.
472 * Note that we only take a single constraint (the one to the heaviest atom)
473 * into account. If an atom has multiple constraints, this will result in
474 * an overestimate of the displacement, which gives a larger drift and buffer.
476 void constrained_atom_sigma2(real kT_fac
,
477 const atom_nonbonded_kinetic_prop_t
*prop
,
481 /* Here we decompose the motion of a constrained atom into two
482 * components: rotation around the COM and translation of the COM.
485 /* Determine the variance of the arc length for the two rotational DOFs */
486 real massFraction
= prop
->con_mass
/(prop
->mass
+ prop
->con_mass
);
487 real sigma2_rot
= kT_fac
*massFraction
/prop
->mass
;
489 /* The distance from the atom to the COM, i.e. the rotational arm */
490 real comDistance
= prop
->con_len
*massFraction
;
492 /* The variance relative to the arm */
493 real sigma2_rel
= sigma2_rot
/gmx::square(comDistance
);
495 /* For sigma2_rel << 1 we don't notice the rotational effect and
496 * we have a normal, Gaussian displacement distribution.
497 * For larger sigma2_rel the displacement is much less, in fact it can
498 * not exceed 2*comDistance. We can calculate MSD/arm^2 as:
499 * integral_x=0-inf distance2(x) x/sigma2_rel exp(-x^2/(2 sigma2_rel)) dx
500 * where x is angular displacement and distance2(x) is the distance^2
501 * between points at angle 0 and x:
502 * distance2(x) = (sin(x) - sin(0))^2 + (cos(x) - cos(0))^2
503 * The limiting value of this MSD is 2, which is also the value for
504 * a uniform rotation distribution that would be reached at long time.
505 * The maximum is 2.5695 at sigma2_rel = 4.5119.
506 * We approximate this integral with a rational polynomial with
507 * coefficients from a Taylor expansion. This approximation is an
508 * overestimate for all values of sigma2_rel. Its maximum value
509 * of 2.6491 is reached at sigma2_rel = sqrt(45/2) = 4.7434.
510 * We keep the approximation constant after that.
511 * We use this approximate MSD as the variance for a Gaussian distribution.
513 * NOTE: For any sensible buffer tolerance this will result in a (large)
514 * overestimate of the buffer size, since the Gaussian has a long tail,
515 * whereas the actual distribution can not reach values larger than 2.
517 /* Coeffients obtained from a Taylor expansion */
518 const real a
= 1.0/3.0;
519 const real b
= 2.0/45.0;
521 /* Our approximation is constant after sigma2_rel = 1/sqrt(b) */
522 sigma2_rel
= std::min(sigma2_rel
, 1/std::sqrt(b
));
524 /* Compute the approximate sigma^2 for 2D motion due to the rotation */
525 *sigma2_2d
= gmx::square(comDistance
)*
526 sigma2_rel
/(1 + a
*sigma2_rel
+ b
*gmx::square(sigma2_rel
));
528 /* The constrained atom also moves (in 3D) with the COM of both atoms */
529 *sigma2_3d
= kT_fac
/(prop
->mass
+ prop
->con_mass
);
532 static void get_atom_sigma2(real kT_fac
,
533 const atom_nonbonded_kinetic_prop_t
*prop
,
539 /* Complicated constraint calculation in a separate function */
540 constrained_atom_sigma2(kT_fac
, prop
, sigma2_2d
, sigma2_3d
);
544 /* Unconstrained atom: trivial */
546 *sigma2_3d
= kT_fac
/prop
->mass
;
550 static void approx_2dof(real s2
, real x
, real
*shift
, real
*scale
)
552 /* A particle with 1 DOF constrained has 2 DOFs instead of 3.
553 * This code is also used for particles with multiple constraints,
554 * in which case we overestimate the displacement.
555 * The 2DOF distribution is sqrt(pi/2)*erfc(r/(sqrt(2)*s))/(2*s).
556 * We approximate this with scale*Gaussian(s,r+shift),
557 * by matching the distribution value and derivative at x.
558 * This is a tight overestimate for all r>=0 at any s and x.
562 ex
= std::exp(-x
*x
/(2*s2
));
563 er
= std::erfc(x
/std::sqrt(2*s2
));
565 *shift
= -x
+ std::sqrt(2*s2
/M_PI
)*ex
/er
;
566 *scale
= 0.5*M_PI
*std::exp(ex
*ex
/(M_PI
*er
*er
))*er
;
569 // Returns an (over)estimate of the energy drift for a single atom pair,
570 // given the kinetic properties, displacement variances and list buffer.
571 static real
energyDriftAtomPair(const atom_nonbonded_kinetic_prop_t
*prop_i
,
572 const atom_nonbonded_kinetic_prop_t
*prop_j
,
573 real s2
, real s2i_2d
, real s2j_2d
,
575 const pot_derivatives_t
*der
)
577 // For relatively small arguments erfc() is so small that if will be 0.0
578 // when stored in a float. We set an argument limit of 8 (Erfc(8)=1e-29),
579 // such that we can divide by erfc and have some space left for arithmetic.
580 const real erfc_arg_max
= 8.0;
587 if (rsh
*rsh
> 2*s2
*erfc_arg_max
*erfc_arg_max
)
589 // Below we calculate c_erfc = 0.5*erfc(rsh/sqrt(2*s2))
590 // When rsh/sqrt(2*s2) increases, this erfc will be the first
591 // result that underflows and becomes 0.0. To avoid this,
592 // we set c_exp=0 and c_erfc=0 for large arguments.
593 // This also avoids NaN in approx_2dof().
594 // In any relevant case this has no effect on the results,
595 // since c_exp < 6e-29, so the displacement is completely
596 // negligible for such atom pairs (and an overestimate).
597 // In nearly all use cases, there will be other atom pairs
598 // that contribute much more to the total, so zeroing
599 // this particular contribution has no effect at all.
605 /* For constraints: adapt r and scaling for the Gaussian */
610 approx_2dof(s2i_2d
, r_buffer
*s2i_2d
/s2
, &sh
, &sc
);
618 approx_2dof(s2j_2d
, r_buffer
*s2j_2d
/s2
, &sh
, &sc
);
623 /* Exact contribution of an atom pair with Gaussian displacement
624 * with sigma s to the energy drift for a potential with
625 * derivative -md and second derivative dd at the cut-off.
626 * The only catch is that for potentials that change sign
627 * near the cut-off there could be an unlucky compensation
628 * of positive and negative energy drift.
629 * Such potentials are extremely rare though.
631 * Note that pot has unit energy*length, as the linear
632 * atom density still needs to be put in.
634 c_exp
= std::exp(-rsh
*rsh
/(2*s2
))/std::sqrt(2*M_PI
);
635 c_erfc
= 0.5*std::erfc(rsh
/(std::sqrt(2*s2
)));
637 real s
= std::sqrt(s2
);
641 der
->md1
/2*((rsh2
+ s2
)*c_erfc
- rsh
*s
*c_exp
);
643 der
->d2
/6*(s
*(rsh2
+ 2*s2
)*c_exp
- rsh
*(rsh2
+ 3*s2
)*c_erfc
);
645 der
->md3
/24*((rsh2
*rsh2
+ 6*rsh2
*s2
+ 3*s2
*s2
)*c_erfc
- rsh
*s
*(rsh2
+ 5*s2
)*c_exp
);
647 return pot1
+ pot2
+ pot3
;
650 static real
energyDrift(const verletbuf_atomtype_t
*att
, int natt
,
651 const gmx_ffparams_t
*ffp
,
653 const pot_derivatives_t
*ljDisp
,
654 const pot_derivatives_t
*ljRep
,
655 const pot_derivatives_t
*elec
,
656 real rlj
, real rcoulomb
,
657 real rlist
, real boxvol
)
659 double drift_tot
= 0;
663 /* No atom displacements: no drift, avoid division by 0 */
667 // Here add up the contribution of all atom pairs in the system to
668 // (estimated) energy drift by looping over all atom type pairs.
669 for (int i
= 0; i
< natt
; i
++)
671 // Get the thermal displacement variance for the i-atom type
672 const atom_nonbonded_kinetic_prop_t
*prop_i
= &att
[i
].prop
;
674 get_atom_sigma2(kT_fac
, prop_i
, &s2i_2d
, &s2i_3d
);
676 for (int j
= i
; j
< natt
; j
++)
678 // Get the thermal displacement variance for the j-atom type
679 const atom_nonbonded_kinetic_prop_t
*prop_j
= &att
[j
].prop
;
681 get_atom_sigma2(kT_fac
, prop_j
, &s2j_2d
, &s2j_3d
);
683 /* Add up the up to four independent variances */
684 real s2
= s2i_2d
+ s2i_3d
+ s2j_2d
+ s2j_3d
;
686 // Set -V', V'' and -V''' at the cut-off for LJ */
687 real c6
= ffp
->iparams
[prop_i
->type
*ffp
->atnr
+ prop_j
->type
].lj
.c6
;
688 real c12
= ffp
->iparams
[prop_i
->type
*ffp
->atnr
+ prop_j
->type
].lj
.c12
;
689 pot_derivatives_t lj
;
690 lj
.md1
= c6
*ljDisp
->md1
+ c12
*ljRep
->md1
;
691 lj
.d2
= c6
*ljDisp
->d2
+ c12
*ljRep
->d2
;
692 lj
.md3
= c6
*ljDisp
->md3
+ c12
*ljRep
->md3
;
694 real pot_lj
= energyDriftAtomPair(prop_i
, prop_j
,
699 // Set -V' and V'' at the cut-off for Coulomb
700 pot_derivatives_t elec_qq
;
701 elec_qq
.md1
= elec
->md1
*prop_i
->q
*prop_j
->q
;
702 elec_qq
.d2
= elec
->d2
*prop_i
->q
*prop_j
->q
;
705 real pot_q
= energyDriftAtomPair(prop_i
, prop_j
,
710 // Note that attractive and repulsive potentials for individual
711 // pairs can partially cancel.
712 real pot
= pot_lj
+ pot_q
;
714 /* Multiply by the number of atom pairs */
717 pot
*= (double)att
[i
].n
*(att
[i
].n
- 1)/2;
721 pot
*= (double)att
[i
].n
*att
[j
].n
;
723 /* We need the line density to get the energy drift of the system.
724 * The effective average r^2 is close to (rlist+sigma)^2.
726 pot
*= 4*M_PI
*gmx::square(rlist
+ std::sqrt(s2
))/boxvol
;
728 /* Add the unsigned drift to avoid cancellation of errors */
729 drift_tot
+= std::abs(pot
);
736 static real
surface_frac(int cluster_size
, real particle_distance
, real rlist
)
740 if (rlist
< 0.5*particle_distance
)
742 /* We have non overlapping spheres */
746 /* Half the inter-particle distance relative to rlist */
747 d
= 0.5*particle_distance
/rlist
;
749 /* Determine the area of the surface at distance rlist to the closest
750 * particle, relative to surface of a sphere of radius rlist.
751 * The formulas below assume close to cubic cells for the pair search grid,
752 * which the pair search code tries to achieve.
753 * Note that in practice particle distances will not be delta distributed,
754 * but have some spread, often involving shorter distances,
755 * as e.g. O-H bonds in a water molecule. Thus the estimates below will
756 * usually be slightly too high and thus conservative.
758 switch (cluster_size
)
761 /* One particle: trivial */
765 /* Two particles: two spheres at fractional distance 2*a */
769 /* We assume a perfect, symmetric tetrahedron geometry.
770 * The surface around a tetrahedron is too complex for a full
771 * analytical solution, so we use a Taylor expansion.
773 area_rel
= (1.0 + 1/M_PI
*(6*std::acos(1/std::sqrt(3))*d
+
774 std::sqrt(3)*d
*d
*(1.0 +
777 83.0/756.0*d
*d
*d
*d
*d
*d
)));
780 gmx_incons("surface_frac called with unsupported cluster_size");
784 return area_rel
/cluster_size
;
787 /* Returns the negative of the third derivative of a potential r^-p
788 * with a force-switch function, evaluated at the cut-off rc.
790 static real
md3_force_switch(real p
, real rswitch
, real rc
)
792 /* The switched force function is:
793 * p*r^-(p+1) + a*(r - rswitch)^2 + b*(r - rswitch)^3
796 real md3_pot
, md3_sw
;
798 a
= -((p
+ 4)*rc
- (p
+ 1)*rswitch
)/(pow(rc
, p
+2)*gmx::square(rc
-rswitch
));
799 b
= ((p
+ 3)*rc
- (p
+ 1)*rswitch
)/(pow(rc
, p
+2)*gmx::power3(rc
-rswitch
));
801 md3_pot
= (p
+ 2)*(p
+ 1)*p
*pow(rc
, p
+3);
802 md3_sw
= 2*a
+ 6*b
*(rc
- rswitch
);
804 return md3_pot
+ md3_sw
;
807 void calc_verlet_buffer_size(const gmx_mtop_t
*mtop
, real boxvol
,
808 const t_inputrec
*ir
,
811 real reference_temperature
,
812 const VerletbufListSetup
*list_setup
,
819 real particle_distance
;
820 real nb_clust_frac_pairs_not_in_list_at_cutoff
;
822 verletbuf_atomtype_t
*att
= nullptr;
825 real kT_fac
, mass_min
;
830 if (!EI_DYNAMICS(ir
->eI
))
832 gmx_incons("Can only determine the Verlet buffer size for integrators that perform dynamics");
834 if (ir
->verletbuf_tol
<= 0)
836 gmx_incons("The Verlet buffer tolerance needs to be larger than zero");
839 if (reference_temperature
< 0)
841 if (EI_MD(ir
->eI
) && ir
->etc
== etcNO
)
843 /* This case should be handled outside calc_verlet_buffer_size */
844 gmx_incons("calc_verlet_buffer_size called with an NVE ensemble and reference_temperature < 0");
847 /* We use the maximum temperature with multiple T-coupl groups.
848 * We could use a per particle temperature, but since particles
849 * interact, this might underestimate the buffer size.
851 reference_temperature
= 0;
852 for (i
= 0; i
< ir
->opts
.ngtc
; i
++)
854 if (ir
->opts
.tau_t
[i
] >= 0)
856 reference_temperature
= std::max(reference_temperature
,
862 /* Resolution of the buffer size */
865 env
= getenv("GMX_VERLET_BUFFER_RES");
868 sscanf(env
, "%lf", &resolution
);
871 /* In an atom wise pair-list there would be no pairs in the list
872 * beyond the pair-list cut-off.
873 * However, we use a pair-list of groups vs groups of atoms.
874 * For groups of 4 atoms, the parallelism of SSE instructions, only
875 * 10% of the atoms pairs are not in the list just beyond the cut-off.
876 * As this percentage increases slowly compared to the decrease of the
877 * Gaussian displacement distribution over this range, we can simply
878 * reduce the drift by this fraction.
879 * For larger groups, e.g. of 8 atoms, this fraction will be lower,
880 * so then buffer size will be on the conservative (large) side.
882 * Note that the formulas used here do not take into account
883 * cancellation of errors which could occur by missing both
884 * attractive and repulsive interactions.
886 * The only major assumption is homogeneous particle distribution.
887 * For an inhomogeneous system, such as a liquid-vapor system,
888 * the buffer will be underestimated. The actual energy drift
889 * will be higher by the factor: local/homogeneous particle density.
891 * The results of this estimate have been checked againt simulations.
892 * In most cases the real drift differs by less than a factor 2.
895 /* Worst case assumption: HCP packing of particles gives largest distance */
896 particle_distance
= std::cbrt(boxvol
*std::sqrt(2)/mtop
->natoms
);
898 get_verlet_buffer_atomtypes(mtop
, &att
, &natt
, n_nonlin_vsite
);
899 assert(att
!= NULL
&& natt
>= 0);
903 fprintf(debug
, "particle distance assuming HCP packing: %f nm\n",
905 fprintf(debug
, "energy drift atom types: %d\n", natt
);
908 pot_derivatives_t ljDisp
= { 0, 0, 0 };
909 pot_derivatives_t ljRep
= { 0, 0, 0 };
910 real repPow
= mtop
->ffparams
.reppow
;
912 if (ir
->vdwtype
== evdwCUT
)
914 real sw_range
, md3_pswf
;
916 switch (ir
->vdw_modifier
)
919 case eintmodPOTSHIFT
:
920 /* -dV/dr of -r^-6 and r^-reppow */
921 ljDisp
.md1
= -6*std::pow(ir
->rvdw
, -7.0);
922 ljRep
.md1
= repPow
*std::pow(ir
->rvdw
, -(repPow
+ 1));
923 /* The contribution of the higher derivatives is negligible */
925 case eintmodFORCESWITCH
:
926 /* At the cut-off: V=V'=V''=0, so we use only V''' */
927 ljDisp
.md3
= -md3_force_switch(6.0, ir
->rvdw_switch
, ir
->rvdw
);
928 ljRep
.md3
= md3_force_switch(repPow
, ir
->rvdw_switch
, ir
->rvdw
);
930 case eintmodPOTSWITCH
:
931 /* At the cut-off: V=V'=V''=0.
932 * V''' is given by the original potential times
933 * the third derivative of the switch function.
935 sw_range
= ir
->rvdw
- ir
->rvdw_switch
;
936 md3_pswf
= 60.0/gmx::power3(sw_range
);
938 ljDisp
.md3
= -std::pow(ir
->rvdw
, -6.0 )*md3_pswf
;
939 ljRep
.md3
= std::pow(ir
->rvdw
, -repPow
)*md3_pswf
;
942 gmx_incons("Unimplemented VdW modifier");
945 else if (EVDW_PME(ir
->vdwtype
))
947 real b
= calc_ewaldcoeff_lj(ir
->rvdw
, ir
->ewald_rtol_lj
);
953 // -dV/dr of g(br)*r^-6 [where g(x) = exp(-x^2)(1+x^2+x^4/2),
954 // see LJ-PME equations in manual] and r^-reppow
955 ljDisp
.md1
= -std::exp(-br2
)*(br6
+ 3.0*br4
+ 6.0*br2
+ 6.0)*std::pow(r
, -7.0);
956 ljRep
.md1
= repPow
*pow(r
, -(repPow
+ 1));
957 // The contribution of the higher derivatives is negligible
961 gmx_fatal(FARGS
, "Energy drift calculation is only implemented for plain cut-off Lennard-Jones interactions");
964 elfac
= ONE_4PI_EPS0
/ir
->epsilon_r
;
966 // Determine the 1st and 2nd derivative for the electostatics
967 pot_derivatives_t elec
= { 0, 0, 0 };
969 if (ir
->coulombtype
== eelCUT
|| EEL_RF(ir
->coulombtype
))
973 if (ir
->coulombtype
== eelCUT
)
980 eps_rf
= ir
->epsilon_rf
/ir
->epsilon_r
;
983 k_rf
= (eps_rf
- ir
->epsilon_r
)/( gmx::power3(ir
->rcoulomb
) * (2*eps_rf
+ ir
->epsilon_r
) );
987 /* epsilon_rf = infinity */
988 k_rf
= 0.5/gmx::power3(ir
->rcoulomb
);
994 elec
.md1
= elfac
*(1.0/gmx::square(ir
->rcoulomb
) - 2*k_rf
*ir
->rcoulomb
);
996 elec
.d2
= elfac
*(2.0/gmx::power3(ir
->rcoulomb
) + 2*k_rf
);
998 else if (EEL_PME(ir
->coulombtype
) || ir
->coulombtype
== eelEWALD
)
1002 b
= calc_ewaldcoeff_q(ir
->rcoulomb
, ir
->ewald_rtol
);
1005 elec
.md1
= elfac
*(b
*std::exp(-br
*br
)*M_2_SQRTPI
/rc
+ std::erfc(br
)/(rc
*rc
));
1006 elec
.d2
= elfac
/(rc
*rc
)*(2*b
*(1 + br
*br
)*std::exp(-br
*br
)*M_2_SQRTPI
+ 2*std::erfc(br
)/rc
);
1010 gmx_fatal(FARGS
, "Energy drift calculation is only implemented for Reaction-Field and Ewald electrostatics");
1013 /* Determine the variance of the atomic displacement
1014 * over list_lifetime steps: kT_fac
1015 * For inertial dynamics (not Brownian dynamics) the mass factor
1016 * is not included in kT_fac, it is added later.
1020 /* Get the displacement distribution from the random component only.
1021 * With accurate integration the systematic (force) displacement
1022 * should be negligible (unless nstlist is extremely large, which
1023 * you wouldn't do anyhow).
1025 kT_fac
= 2*BOLTZ
*reference_temperature
*list_lifetime
*ir
->delta_t
;
1026 if (ir
->bd_fric
> 0)
1028 /* This is directly sigma^2 of the displacement */
1029 kT_fac
/= ir
->bd_fric
;
1031 /* Set the masses to 1 as kT_fac is the full sigma^2,
1032 * but we divide by m in ener_drift().
1034 for (i
= 0; i
< natt
; i
++)
1036 att
[i
].prop
.mass
= 1;
1043 /* Per group tau_t is not implemented yet, use the maximum */
1044 tau_t
= ir
->opts
.tau_t
[0];
1045 for (i
= 1; i
< ir
->opts
.ngtc
; i
++)
1047 tau_t
= std::max(tau_t
, ir
->opts
.tau_t
[i
]);
1051 /* This kT_fac needs to be divided by the mass to get sigma^2 */
1056 kT_fac
= BOLTZ
*reference_temperature
*gmx::square(list_lifetime
*ir
->delta_t
);
1059 mass_min
= att
[0].prop
.mass
;
1060 for (i
= 1; i
< natt
; i
++)
1062 mass_min
= std::min(mass_min
, att
[i
].prop
.mass
);
1067 fprintf(debug
, "Derivatives of non-bonded potentials at the cut-off:\n");
1068 fprintf(debug
, "LJ disp. -V' %9.2e V'' %9.2e -V''' %9.2e\n", ljDisp
.md1
, ljDisp
.d2
, ljDisp
.md3
);
1069 fprintf(debug
, "LJ rep. -V' %9.2e V'' %9.2e -V''' %9.2e\n", ljRep
.md1
, ljRep
.d2
, ljRep
.md3
);
1070 fprintf(debug
, "Electro. -V' %9.2e V'' %9.2e\n", elec
.md1
, elec
.d2
);
1071 fprintf(debug
, "sqrt(kT_fac) %f\n", std::sqrt(kT_fac
));
1072 fprintf(debug
, "mass_min %f\n", mass_min
);
1075 /* Search using bisection */
1077 /* The drift will be neglible at 5 times the max sigma */
1078 ib1
= (int)(5*2*std::sqrt(kT_fac
/mass_min
)/resolution
) + 1;
1079 while (ib1
- ib0
> 1)
1083 rl
= std::max(ir
->rvdw
, ir
->rcoulomb
) + rb
;
1085 /* Calculate the average energy drift at the last step
1086 * of the nstlist steps at which the pair-list is used.
1088 drift
= energyDrift(att
, natt
, &mtop
->ffparams
,
1090 &ljDisp
, &ljRep
, &elec
,
1091 ir
->rvdw
, ir
->rcoulomb
,
1094 /* Correct for the fact that we are using a Ni x Nj particle pair list
1095 * and not a 1 x 1 particle pair list. This reduces the drift.
1097 /* We don't have a formula for 8 (yet), use 4 which is conservative */
1098 nb_clust_frac_pairs_not_in_list_at_cutoff
=
1099 surface_frac(std::min(list_setup
->cluster_size_i
, 4),
1100 particle_distance
, rl
)*
1101 surface_frac(std::min(list_setup
->cluster_size_j
, 4),
1102 particle_distance
, rl
);
1103 drift
*= nb_clust_frac_pairs_not_in_list_at_cutoff
;
1105 /* Convert the drift to drift per unit time per atom */
1106 drift
/= nstlist
*ir
->delta_t
*mtop
->natoms
;
1110 fprintf(debug
, "ib %3d %3d %3d rb %.3f %dx%d fac %.3f drift %.1e\n",
1112 list_setup
->cluster_size_i
, list_setup
->cluster_size_j
,
1113 nb_clust_frac_pairs_not_in_list_at_cutoff
,
1117 if (std::abs(drift
) > ir
->verletbuf_tol
)
1129 *rlist
= std::max(ir
->rvdw
, ir
->rcoulomb
) + ib1
*resolution
;