2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2012,2013,2014,2015,2016,2017, 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/math/calculate-ewald-splitting-coefficient.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 order 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.
99 int type
; /* type (used for LJ parameters) */
101 gmx_bool bConstr
; /* constrained, if TRUE, use #DOF=2 iso 3 */
102 real con_mass
; /* mass of heaviest atom connected by constraints */
103 real con_len
; /* constraint length to the heaviest atom */
104 } atom_nonbonded_kinetic_prop_t
;
106 /* Struct for unique atom type for calculating the energy drift.
107 * The atom displacement depends on mass and constraints.
108 * The energy jump for given distance depend on LJ type and q.
112 atom_nonbonded_kinetic_prop_t prop
; /* non-bonded and kinetic atom prop. */
113 int n
; /* #atoms of this type in the system */
114 } verletbuf_atomtype_t
;
116 // Struct for derivatives of a non-bonded interaction potential
119 real md1
; // -V' at the cutoff
120 real d2
; // V'' at the cutoff
121 real md3
; // -V''' at the cutoff
124 void verletbuf_get_list_setup(gmx_bool gmx_unused bSIMD
,
126 verletbuf_list_setup_t
*list_setup
)
128 /* When calling this function we often don't know which kernel type we
129 * are going to use. W choose the kernel type with the smallest possible
130 * i- and j-cluster sizes, so we potentially overestimate, but never
131 * underestimate, the buffer drift.
132 * Note that the current buffer estimation code only handles clusters
133 * of size 1, 2 or 4, so for 4x8 or 8x8 we use the estimate for 4x4.
138 /* The CUDA kernels split the j-clusters in two halves */
139 list_setup
->cluster_size_i
= nbnxn_kernel_to_cluster_i_size(nbnxnk8x8x8_GPU
);
140 list_setup
->cluster_size_j
= nbnxn_kernel_to_cluster_j_size(nbnxnk8x8x8_GPU
)/2;
146 kernel_type
= nbnxnk4x4_PlainC
;
151 #ifdef GMX_NBNXN_SIMD_2XNN
152 /* We use the smallest cluster size to be on the safe side */
153 kernel_type
= nbnxnk4xN_SIMD_2xNN
;
155 kernel_type
= nbnxnk4xN_SIMD_4xN
;
160 list_setup
->cluster_size_i
= nbnxn_kernel_to_cluster_i_size(kernel_type
);
161 list_setup
->cluster_size_j
= nbnxn_kernel_to_cluster_j_size(kernel_type
);
166 atom_nonbonded_kinetic_prop_equal(const atom_nonbonded_kinetic_prop_t
*prop1
,
167 const atom_nonbonded_kinetic_prop_t
*prop2
)
169 return (prop1
->mass
== prop2
->mass
&&
170 prop1
->type
== prop2
->type
&&
171 prop1
->q
== prop2
->q
&&
172 prop1
->bConstr
== prop2
->bConstr
&&
173 prop1
->con_mass
== prop2
->con_mass
&&
174 prop1
->con_len
== prop2
->con_len
);
177 static void add_at(verletbuf_atomtype_t
**att_p
, int *natt_p
,
178 const atom_nonbonded_kinetic_prop_t
*prop
,
181 verletbuf_atomtype_t
*att
;
186 /* Ignore massless particles */
194 while (i
< natt
&& !atom_nonbonded_kinetic_prop_equal(prop
, &att
[i
].prop
))
206 srenew(*att_p
, *natt_p
);
207 (*att_p
)[i
].prop
= *prop
;
208 (*att_p
)[i
].n
= nmol
;
212 static void get_vsite_masses(const gmx_moltype_t
*moltype
,
213 const gmx_ffparams_t
*ffparams
,
222 /* Check for virtual sites, determine mass from constructing atoms */
223 for (ft
= 0; ft
< F_NRE
; ft
++)
227 il
= &moltype
->ilist
[ft
];
229 for (i
= 0; i
< il
->nr
; i
+= 1+NRAL(ft
))
232 real inv_mass
, coeff
, m_aj
;
235 ip
= &ffparams
->iparams
[il
->iatoms
[i
]];
237 a1
= il
->iatoms
[i
+1];
241 /* Only vsiten can have more than four
242 constructing atoms, so NRAL(ft) <= 5 */
245 const int maxj
= NRAL(ft
);
249 for (j
= 1; j
< maxj
; j
++)
251 cam
[j
] = moltype
->atoms
.atom
[il
->iatoms
[i
+1+j
]].m
;
254 cam
[j
] = vsite_m
[il
->iatoms
[i
+1+j
]];
258 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.",
260 interaction_function
[ft
].longname
,
261 il
->iatoms
[i
+1+j
]+1);
269 vsite_m
[a1
] = (cam
[1]*cam
[2])/(cam
[2]*gmx::square(1-ip
->vsite
.a
) + cam
[1]*gmx::square(ip
->vsite
.a
));
273 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
));
276 gmx_incons("Invalid vsite type");
279 /* Use the mass of the lightest constructing atom.
280 * This is an approximation.
281 * If the distance of the virtual site to the
282 * constructing atom is less than all distances
283 * between constructing atoms, this is a safe
284 * over-estimate of the displacement of the vsite.
285 * This condition holds for all H mass replacement
286 * vsite constructions, except for SP2/3 groups.
287 * In SP3 groups one H will have a F_VSITE3
288 * construction, so even there the total drift
289 * estimate shouldn't be far off.
291 vsite_m
[a1
] = cam
[1];
292 for (j
= 2; j
< maxj
; j
++)
294 vsite_m
[a1
] = std::min(vsite_m
[a1
], cam
[j
]);
307 for (j
= 0; j
< 3*ffparams
->iparams
[il
->iatoms
[i
]].vsiten
.n
; j
+= 3)
309 aj
= il
->iatoms
[i
+j
+2];
310 coeff
= ffparams
->iparams
[il
->iatoms
[i
+j
]].vsiten
.a
;
311 if (moltype
->atoms
.atom
[aj
].ptype
== eptVSite
)
317 m_aj
= moltype
->atoms
.atom
[aj
].m
;
321 gmx_incons("The mass of a vsiten constructing atom is <= 0");
323 inv_mass
+= coeff
*coeff
/m_aj
;
325 vsite_m
[a1
] = 1/inv_mass
;
326 /* Correct for loop increment of i */
327 i
+= j
- 1 - NRAL(ft
);
331 fprintf(debug
, "atom %4d %-20s mass %6.3f\n",
332 a1
, interaction_function
[ft
].longname
, vsite_m
[a1
]);
339 static void get_verlet_buffer_atomtypes(const gmx_mtop_t
*mtop
,
340 verletbuf_atomtype_t
**att_p
,
344 verletbuf_atomtype_t
*att
;
346 int mb
, nmol
, ft
, i
, a1
, a2
, a3
, a
;
347 const t_atoms
*atoms
;
350 atom_nonbonded_kinetic_prop_t
*prop
;
352 int n_nonlin_vsite_mol
;
357 if (n_nonlin_vsite
!= nullptr)
362 for (mb
= 0; mb
< mtop
->nmolblock
; mb
++)
364 nmol
= mtop
->molblock
[mb
].nmol
;
366 atoms
= &mtop
->moltype
[mtop
->molblock
[mb
].type
].atoms
;
368 /* Check for constraints, as they affect the kinetic energy.
369 * For virtual sites we need the masses and geometry of
370 * the constructing atoms to determine their velocity distribution.
372 snew(prop
, atoms
->nr
);
373 snew(vsite_m
, atoms
->nr
);
375 for (ft
= F_CONSTR
; ft
<= F_CONSTRNC
; ft
++)
377 il
= &mtop
->moltype
[mtop
->molblock
[mb
].type
].ilist
[ft
];
379 for (i
= 0; i
< il
->nr
; i
+= 1+NRAL(ft
))
381 ip
= &mtop
->ffparams
.iparams
[il
->iatoms
[i
]];
382 a1
= il
->iatoms
[i
+1];
383 a2
= il
->iatoms
[i
+2];
384 if (atoms
->atom
[a2
].m
> prop
[a1
].con_mass
)
386 prop
[a1
].con_mass
= atoms
->atom
[a2
].m
;
387 prop
[a1
].con_len
= ip
->constr
.dA
;
389 if (atoms
->atom
[a1
].m
> prop
[a2
].con_mass
)
391 prop
[a2
].con_mass
= atoms
->atom
[a1
].m
;
392 prop
[a2
].con_len
= ip
->constr
.dA
;
397 il
= &mtop
->moltype
[mtop
->molblock
[mb
].type
].ilist
[F_SETTLE
];
399 for (i
= 0; i
< il
->nr
; i
+= 1+NRAL(F_SETTLE
))
401 ip
= &mtop
->ffparams
.iparams
[il
->iatoms
[i
]];
402 a1
= il
->iatoms
[i
+1];
403 a2
= il
->iatoms
[i
+2];
404 a3
= il
->iatoms
[i
+3];
405 /* Usually the mass of a1 (usually oxygen) is larger than a2/a3.
406 * If this is not the case, we overestimate the displacement,
407 * which leads to a larger buffer (ok since this is an exotic case).
409 prop
[a1
].con_mass
= atoms
->atom
[a2
].m
;
410 prop
[a1
].con_len
= ip
->settle
.doh
;
412 prop
[a2
].con_mass
= atoms
->atom
[a1
].m
;
413 prop
[a2
].con_len
= ip
->settle
.doh
;
415 prop
[a3
].con_mass
= atoms
->atom
[a1
].m
;
416 prop
[a3
].con_len
= ip
->settle
.doh
;
419 get_vsite_masses(&mtop
->moltype
[mtop
->molblock
[mb
].type
],
422 &n_nonlin_vsite_mol
);
423 if (n_nonlin_vsite
!= nullptr)
425 *n_nonlin_vsite
+= nmol
*n_nonlin_vsite_mol
;
428 for (a
= 0; a
< atoms
->nr
; a
++)
430 if (atoms
->atom
[a
].ptype
== eptVSite
)
432 prop
[a
].mass
= vsite_m
[a
];
436 prop
[a
].mass
= atoms
->atom
[a
].m
;
438 prop
[a
].type
= atoms
->atom
[a
].type
;
439 prop
[a
].q
= atoms
->atom
[a
].q
;
440 /* We consider an atom constrained, #DOF=2, when it is
441 * connected with constraints to (at least one) atom with
442 * a mass of more than 0.4x its own mass. This is not a critical
443 * parameter, since with roughly equal masses the unconstrained
444 * and constrained displacement will not differ much (and both
445 * overestimate the displacement).
447 prop
[a
].bConstr
= (prop
[a
].con_mass
> 0.4*prop
[a
].mass
);
449 add_at(&att
, &natt
, &prop
[a
], nmol
);
458 for (a
= 0; a
< natt
; a
++)
460 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",
461 a
, att
[a
].prop
.mass
, att
[a
].prop
.type
, att
[a
].prop
.q
,
462 att
[a
].prop
.bConstr
, att
[a
].prop
.con_mass
, att
[a
].prop
.con_len
,
471 /* This function computes two components of the estimate of the variance
472 * in the displacement of one atom in a system of two constrained atoms.
473 * Returns in sigma2_2d the variance due to rotation of the constrained
474 * atom around the atom to which it constrained.
475 * Returns in sigma2_3d the variance due to displacement of the COM
476 * of the whole system of the two constrained atoms.
478 * Note that we only take a single constraint (the one to the heaviest atom)
479 * into account. If an atom has multiple constraints, this will result in
480 * an overestimate of the displacement, which gives a larger drift and buffer.
482 static void constrained_atom_sigma2(real kT_fac
,
483 const atom_nonbonded_kinetic_prop_t
*prop
,
492 /* Here we decompose the motion of a constrained atom into two
493 * components: rotation around the COM and translation of the COM.
496 /* Determine the variance for the displacement of the rotational mode */
497 sigma2_rot
= kT_fac
/(prop
->mass
*(prop
->mass
+ prop
->con_mass
)/prop
->con_mass
);
499 /* The distance from the atom to the COM, i.e. the rotational arm */
500 com_dist
= prop
->con_len
*prop
->con_mass
/(prop
->mass
+ prop
->con_mass
);
502 /* The variance relative to the arm */
503 sigma2_rel
= sigma2_rot
/(com_dist
*com_dist
);
504 /* At 6 the scaling formula has slope 0,
505 * so we keep sigma2_2d constant after that.
509 /* A constrained atom rotates around the atom it is constrained to.
510 * This results in a smaller linear displacement than for a free atom.
511 * For a perfectly circular displacement, this lowers the displacement
512 * by: 1/arcsin(arc_length)
513 * and arcsin(x) = 1 + x^2/6 + ...
514 * For sigma2_rel<<1 the displacement distribution is erfc
515 * (exact formula is provided below). For larger sigma, it is clear
516 * that the displacement can't be larger than 2*com_dist.
517 * It turns out that the distribution becomes nearly uniform.
518 * For intermediate sigma2_rel, scaling down sigma with the third
519 * order expansion of arcsin with argument sigma_rel turns out
520 * to give a very good approximation of the distribution and variance.
521 * Even for larger values, the variance is only slightly overestimated.
522 * Note that the most relevant displacements are in the long tail.
523 * This rotation approximation always overestimates the tail (which
524 * runs to infinity, whereas it should be <= 2*com_dist).
525 * Thus we always overestimate the drift and the buffer size.
527 scale
= 1/(1 + sigma2_rel
/6);
528 *sigma2_2d
= sigma2_rot
*scale
*scale
;
532 /* sigma_2d is set to the maximum given by the scaling above.
533 * For large sigma2 the real displacement distribution is close
534 * to uniform over -2*con_len to 2*com_dist.
535 * Our erfc with sigma_2d=sqrt(1.5)*com_dist (which means the sigma
536 * of the erfc output distribution is con_dist) overestimates
537 * the variance and additionally has a long tail. This means
538 * we have a (safe) overestimation of the drift.
540 *sigma2_2d
= 1.5*com_dist
*com_dist
;
543 /* The constrained atom also moves (in 3D) with the COM of both atoms */
544 *sigma2_3d
= kT_fac
/(prop
->mass
+ prop
->con_mass
);
547 static void get_atom_sigma2(real kT_fac
,
548 const atom_nonbonded_kinetic_prop_t
*prop
,
554 /* Complicated constraint calculation in a separate function */
555 constrained_atom_sigma2(kT_fac
, prop
, sigma2_2d
, sigma2_3d
);
559 /* Unconstrained atom: trivial */
561 *sigma2_3d
= kT_fac
/prop
->mass
;
565 static void approx_2dof(real s2
, real x
, real
*shift
, real
*scale
)
567 /* A particle with 1 DOF constrained has 2 DOFs instead of 3.
568 * This code is also used for particles with multiple constraints,
569 * in which case we overestimate the displacement.
570 * The 2DOF distribution is sqrt(pi/2)*erfc(r/(sqrt(2)*s))/(2*s).
571 * We approximate this with scale*Gaussian(s,r+shift),
572 * by matching the distribution value and derivative at x.
573 * This is a tight overestimate for all r>=0 at any s and x.
577 ex
= std::exp(-x
*x
/(2*s2
));
578 er
= std::erfc(x
/std::sqrt(2*s2
));
580 *shift
= -x
+ std::sqrt(2*s2
/M_PI
)*ex
/er
;
581 *scale
= 0.5*M_PI
*std::exp(ex
*ex
/(M_PI
*er
*er
))*er
;
584 // Returns an (over)estimate of the energy drift for a single atom pair,
585 // given the kinetic properties, displacement variances and list buffer.
586 static real
energyDriftAtomPair(const atom_nonbonded_kinetic_prop_t
*prop_i
,
587 const atom_nonbonded_kinetic_prop_t
*prop_j
,
588 real s2
, real s2i_2d
, real s2j_2d
,
590 const pot_derivatives_t
*der
)
592 // For relatively small arguments erfc() is so small that if will be 0.0
593 // when stored in a float. We set an argument limit of 8 (Erfc(8)=1e-29),
594 // such that we can divide by erfc and have some space left for arithmetic.
595 const real erfc_arg_max
= 8.0;
602 if (rsh
*rsh
> 2*s2
*erfc_arg_max
*erfc_arg_max
)
604 // Below we calculate c_erfc = 0.5*erfc(rsh/sqrt(2*s2))
605 // When rsh/sqrt(2*s2) increases, this erfc will be the first
606 // result that underflows and becomes 0.0. To avoid this,
607 // we set c_exp=0 and c_erfc=0 for large arguments.
608 // This also avoids NaN in approx_2dof().
609 // In any relevant case this has no effect on the results,
610 // since c_exp < 6e-29, so the displacement is completely
611 // negligible for such atom pairs (and an overestimate).
612 // In nearly all use cases, there will be other atom pairs
613 // that contribute much more to the total, so zeroing
614 // this particular contribution has no effect at all.
620 /* For constraints: adapt r and scaling for the Gaussian */
625 approx_2dof(s2i_2d
, r_buffer
*s2i_2d
/s2
, &sh
, &sc
);
633 approx_2dof(s2j_2d
, r_buffer
*s2j_2d
/s2
, &sh
, &sc
);
638 /* Exact contribution of an atom pair with Gaussian displacement
639 * with sigma s to the energy drift for a potential with
640 * derivative -md and second derivative dd at the cut-off.
641 * The only catch is that for potentials that change sign
642 * near the cut-off there could be an unlucky compensation
643 * of positive and negative energy drift.
644 * Such potentials are extremely rare though.
646 * Note that pot has unit energy*length, as the linear
647 * atom density still needs to be put in.
649 c_exp
= std::exp(-rsh
*rsh
/(2*s2
))/std::sqrt(2*M_PI
);
650 c_erfc
= 0.5*std::erfc(rsh
/(std::sqrt(2*s2
)));
652 real s
= std::sqrt(s2
);
656 der
->md1
/2*((rsh2
+ s2
)*c_erfc
- rsh
*s
*c_exp
);
658 der
->d2
/6*(s
*(rsh2
+ 2*s2
)*c_exp
- rsh
*(rsh2
+ 3*s2
)*c_erfc
);
660 der
->md3
/24*((rsh2
*rsh2
+ 6*rsh2
*s2
+ 3*s2
*s2
)*c_erfc
- rsh
*s
*(rsh2
+ 5*s2
)*c_exp
);
662 return pot1
+ pot2
+ pot3
;
665 static real
energyDrift(const verletbuf_atomtype_t
*att
, int natt
,
666 const gmx_ffparams_t
*ffp
,
668 const pot_derivatives_t
*ljDisp
,
669 const pot_derivatives_t
*ljRep
,
670 const pot_derivatives_t
*elec
,
671 real rlj
, real rcoulomb
,
672 real rlist
, real boxvol
)
674 double drift_tot
= 0;
678 /* No atom displacements: no drift, avoid division by 0 */
682 // Here add up the contribution of all atom pairs in the system to
683 // (estimated) energy drift by looping over all atom type pairs.
684 for (int i
= 0; i
< natt
; i
++)
686 // Get the thermal displacement variance for the i-atom type
687 const atom_nonbonded_kinetic_prop_t
*prop_i
= &att
[i
].prop
;
689 get_atom_sigma2(kT_fac
, prop_i
, &s2i_2d
, &s2i_3d
);
691 for (int j
= i
; j
< natt
; j
++)
693 // Get the thermal displacement variance for the j-atom type
694 const atom_nonbonded_kinetic_prop_t
*prop_j
= &att
[j
].prop
;
696 get_atom_sigma2(kT_fac
, prop_j
, &s2j_2d
, &s2j_3d
);
698 /* Add up the up to four independent variances */
699 real s2
= s2i_2d
+ s2i_3d
+ s2j_2d
+ s2j_3d
;
701 // Set -V', V'' and -V''' at the cut-off for LJ */
702 real c6
= ffp
->iparams
[prop_i
->type
*ffp
->atnr
+ prop_j
->type
].lj
.c6
;
703 real c12
= ffp
->iparams
[prop_i
->type
*ffp
->atnr
+ prop_j
->type
].lj
.c12
;
704 pot_derivatives_t lj
;
705 lj
.md1
= c6
*ljDisp
->md1
+ c12
*ljRep
->md1
;
706 lj
.d2
= c6
*ljDisp
->d2
+ c12
*ljRep
->d2
;
707 lj
.md3
= c6
*ljDisp
->md3
+ c12
*ljRep
->md3
;
709 real pot_lj
= energyDriftAtomPair(prop_i
, prop_j
,
714 // Set -V' and V'' at the cut-off for Coulomb
715 pot_derivatives_t elec_qq
;
716 elec_qq
.md1
= elec
->md1
*prop_i
->q
*prop_j
->q
;
717 elec_qq
.d2
= elec
->d2
*prop_i
->q
*prop_j
->q
;
720 real pot_q
= energyDriftAtomPair(prop_i
, prop_j
,
725 // Note that attractive and repulsive potentials for individual
726 // pairs can partially cancel.
727 real pot
= pot_lj
+ pot_q
;
729 /* Multiply by the number of atom pairs */
732 pot
*= (double)att
[i
].n
*(att
[i
].n
- 1)/2;
736 pot
*= (double)att
[i
].n
*att
[j
].n
;
738 /* We need the line density to get the energy drift of the system.
739 * The effective average r^2 is close to (rlist+sigma)^2.
741 pot
*= 4*M_PI
*gmx::square(rlist
+ std::sqrt(s2
))/boxvol
;
743 /* Add the unsigned drift to avoid cancellation of errors */
744 drift_tot
+= std::abs(pot
);
751 static real
surface_frac(int cluster_size
, real particle_distance
, real rlist
)
755 if (rlist
< 0.5*particle_distance
)
757 /* We have non overlapping spheres */
761 /* Half the inter-particle distance relative to rlist */
762 d
= 0.5*particle_distance
/rlist
;
764 /* Determine the area of the surface at distance rlist to the closest
765 * particle, relative to surface of a sphere of radius rlist.
766 * The formulas below assume close to cubic cells for the pair search grid,
767 * which the pair search code tries to achieve.
768 * Note that in practice particle distances will not be delta distributed,
769 * but have some spread, often involving shorter distances,
770 * as e.g. O-H bonds in a water molecule. Thus the estimates below will
771 * usually be slightly too high and thus conservative.
773 switch (cluster_size
)
776 /* One particle: trivial */
780 /* Two particles: two spheres at fractional distance 2*a */
784 /* We assume a perfect, symmetric tetrahedron geometry.
785 * The surface around a tetrahedron is too complex for a full
786 * analytical solution, so we use a Taylor expansion.
788 area_rel
= (1.0 + 1/M_PI
*(6*std::acos(1/std::sqrt(3))*d
+
789 std::sqrt(3)*d
*d
*(1.0 +
792 83.0/756.0*d
*d
*d
*d
*d
*d
)));
795 gmx_incons("surface_frac called with unsupported cluster_size");
799 return area_rel
/cluster_size
;
802 /* Returns the negative of the third derivative of a potential r^-p
803 * with a force-switch function, evaluated at the cut-off rc.
805 static real
md3_force_switch(real p
, real rswitch
, real rc
)
807 /* The switched force function is:
808 * p*r^-(p+1) + a*(r - rswitch)^2 + b*(r - rswitch)^3
811 real md3_pot
, md3_sw
;
813 a
= -((p
+ 4)*rc
- (p
+ 1)*rswitch
)/(pow(rc
, p
+2)*gmx::square(rc
-rswitch
));
814 b
= ((p
+ 3)*rc
- (p
+ 1)*rswitch
)/(pow(rc
, p
+2)*gmx::power3(rc
-rswitch
));
816 md3_pot
= (p
+ 2)*(p
+ 1)*p
*pow(rc
, p
+3);
817 md3_sw
= 2*a
+ 6*b
*(rc
- rswitch
);
819 return md3_pot
+ md3_sw
;
822 void calc_verlet_buffer_size(const gmx_mtop_t
*mtop
, real boxvol
,
823 const t_inputrec
*ir
,
824 real reference_temperature
,
825 const verletbuf_list_setup_t
*list_setup
,
832 real particle_distance
;
833 real nb_clust_frac_pairs_not_in_list_at_cutoff
;
835 verletbuf_atomtype_t
*att
= nullptr;
838 real kT_fac
, mass_min
;
843 if (!EI_DYNAMICS(ir
->eI
))
845 gmx_incons("Can only determine the Verlet buffer size for integrators that perform dynamics");
847 if (ir
->verletbuf_tol
<= 0)
849 gmx_incons("The Verlet buffer tolerance needs to be larger than zero");
852 if (reference_temperature
< 0)
854 if (EI_MD(ir
->eI
) && ir
->etc
== etcNO
)
856 /* This case should be handled outside calc_verlet_buffer_size */
857 gmx_incons("calc_verlet_buffer_size called with an NVE ensemble and reference_temperature < 0");
860 /* We use the maximum temperature with multiple T-coupl groups.
861 * We could use a per particle temperature, but since particles
862 * interact, this might underestimate the buffer size.
864 reference_temperature
= 0;
865 for (i
= 0; i
< ir
->opts
.ngtc
; i
++)
867 if (ir
->opts
.tau_t
[i
] >= 0)
869 reference_temperature
= std::max(reference_temperature
,
875 /* Resolution of the buffer size */
878 env
= getenv("GMX_VERLET_BUFFER_RES");
881 sscanf(env
, "%lf", &resolution
);
884 /* In an atom wise pair-list there would be no pairs in the list
885 * beyond the pair-list cut-off.
886 * However, we use a pair-list of groups vs groups of atoms.
887 * For groups of 4 atoms, the parallelism of SSE instructions, only
888 * 10% of the atoms pairs are not in the list just beyond the cut-off.
889 * As this percentage increases slowly compared to the decrease of the
890 * Gaussian displacement distribution over this range, we can simply
891 * reduce the drift by this fraction.
892 * For larger groups, e.g. of 8 atoms, this fraction will be lower,
893 * so then buffer size will be on the conservative (large) side.
895 * Note that the formulas used here do not take into account
896 * cancellation of errors which could occur by missing both
897 * attractive and repulsive interactions.
899 * The only major assumption is homogeneous particle distribution.
900 * For an inhomogeneous system, such as a liquid-vapor system,
901 * the buffer will be underestimated. The actual energy drift
902 * will be higher by the factor: local/homogeneous particle density.
904 * The results of this estimate have been checked againt simulations.
905 * In most cases the real drift differs by less than a factor 2.
908 /* Worst case assumption: HCP packing of particles gives largest distance */
909 particle_distance
= std::cbrt(boxvol
*std::sqrt(2)/mtop
->natoms
);
911 get_verlet_buffer_atomtypes(mtop
, &att
, &natt
, n_nonlin_vsite
);
912 assert(att
!= NULL
&& natt
>= 0);
916 fprintf(debug
, "particle distance assuming HCP packing: %f nm\n",
918 fprintf(debug
, "energy drift atom types: %d\n", natt
);
921 pot_derivatives_t ljDisp
= { 0, 0, 0 };
922 pot_derivatives_t ljRep
= { 0, 0, 0 };
923 real repPow
= mtop
->ffparams
.reppow
;
925 if (ir
->vdwtype
== evdwCUT
)
927 real sw_range
, md3_pswf
;
929 switch (ir
->vdw_modifier
)
932 case eintmodPOTSHIFT
:
933 /* -dV/dr of -r^-6 and r^-reppow */
934 ljDisp
.md1
= -6*std::pow(ir
->rvdw
, -7.0);
935 ljRep
.md1
= repPow
*std::pow(ir
->rvdw
, -(repPow
+ 1));
936 /* The contribution of the higher derivatives is negligible */
938 case eintmodFORCESWITCH
:
939 /* At the cut-off: V=V'=V''=0, so we use only V''' */
940 ljDisp
.md3
= -md3_force_switch(6.0, ir
->rvdw_switch
, ir
->rvdw
);
941 ljRep
.md3
= md3_force_switch(repPow
, ir
->rvdw_switch
, ir
->rvdw
);
943 case eintmodPOTSWITCH
:
944 /* At the cut-off: V=V'=V''=0.
945 * V''' is given by the original potential times
946 * the third derivative of the switch function.
948 sw_range
= ir
->rvdw
- ir
->rvdw_switch
;
949 md3_pswf
= 60.0/gmx::power3(sw_range
);
951 ljDisp
.md3
= -std::pow(ir
->rvdw
, -6.0 )*md3_pswf
;
952 ljRep
.md3
= std::pow(ir
->rvdw
, -repPow
)*md3_pswf
;
955 gmx_incons("Unimplemented VdW modifier");
958 else if (EVDW_PME(ir
->vdwtype
))
960 real b
= calc_ewaldcoeff_lj(ir
->rvdw
, ir
->ewald_rtol_lj
);
966 // -dV/dr of g(br)*r^-6 [where g(x) = exp(-x^2)(1+x^2+x^4/2),
967 // see LJ-PME equations in manual] and r^-reppow
968 ljDisp
.md1
= -std::exp(-br2
)*(br6
+ 3.0*br4
+ 6.0*br2
+ 6.0)*std::pow(r
, -7.0);
969 ljRep
.md1
= repPow
*pow(r
, -(repPow
+ 1));
970 // The contribution of the higher derivatives is negligible
974 gmx_fatal(FARGS
, "Energy drift calculation is only implemented for plain cut-off Lennard-Jones interactions");
977 elfac
= ONE_4PI_EPS0
/ir
->epsilon_r
;
979 // Determine the 1st and 2nd derivative for the electostatics
980 pot_derivatives_t elec
= { 0, 0, 0 };
982 if (ir
->coulombtype
== eelCUT
|| EEL_RF(ir
->coulombtype
))
986 if (ir
->coulombtype
== eelCUT
)
993 eps_rf
= ir
->epsilon_rf
/ir
->epsilon_r
;
996 k_rf
= (eps_rf
- ir
->epsilon_r
)/( gmx::power3(ir
->rcoulomb
) * (2*eps_rf
+ ir
->epsilon_r
) );
1000 /* epsilon_rf = infinity */
1001 k_rf
= 0.5/gmx::power3(ir
->rcoulomb
);
1007 elec
.md1
= elfac
*(1.0/gmx::square(ir
->rcoulomb
) - 2*k_rf
*ir
->rcoulomb
);
1009 elec
.d2
= elfac
*(2.0/gmx::power3(ir
->rcoulomb
) + 2*k_rf
);
1011 else if (EEL_PME(ir
->coulombtype
) || ir
->coulombtype
== eelEWALD
)
1015 b
= calc_ewaldcoeff_q(ir
->rcoulomb
, ir
->ewald_rtol
);
1018 elec
.md1
= elfac
*(b
*std::exp(-br
*br
)*M_2_SQRTPI
/rc
+ std::erfc(br
)/(rc
*rc
));
1019 elec
.d2
= elfac
/(rc
*rc
)*(2*b
*(1 + br
*br
)*std::exp(-br
*br
)*M_2_SQRTPI
+ 2*std::erfc(br
)/rc
);
1023 gmx_fatal(FARGS
, "Energy drift calculation is only implemented for Reaction-Field and Ewald electrostatics");
1026 /* Determine the variance of the atomic displacement
1027 * over nstlist-1 steps: kT_fac
1028 * For inertial dynamics (not Brownian dynamics) the mass factor
1029 * is not included in kT_fac, it is added later.
1033 /* Get the displacement distribution from the random component only.
1034 * With accurate integration the systematic (force) displacement
1035 * should be negligible (unless nstlist is extremely large, which
1036 * you wouldn't do anyhow).
1038 kT_fac
= 2*BOLTZ
*reference_temperature
*(ir
->nstlist
-1)*ir
->delta_t
;
1039 if (ir
->bd_fric
> 0)
1041 /* This is directly sigma^2 of the displacement */
1042 kT_fac
/= ir
->bd_fric
;
1044 /* Set the masses to 1 as kT_fac is the full sigma^2,
1045 * but we divide by m in ener_drift().
1047 for (i
= 0; i
< natt
; i
++)
1049 att
[i
].prop
.mass
= 1;
1056 /* Per group tau_t is not implemented yet, use the maximum */
1057 tau_t
= ir
->opts
.tau_t
[0];
1058 for (i
= 1; i
< ir
->opts
.ngtc
; i
++)
1060 tau_t
= std::max(tau_t
, ir
->opts
.tau_t
[i
]);
1064 /* This kT_fac needs to be divided by the mass to get sigma^2 */
1069 kT_fac
= BOLTZ
*reference_temperature
*gmx::square((ir
->nstlist
-1)*ir
->delta_t
);
1072 mass_min
= att
[0].prop
.mass
;
1073 for (i
= 1; i
< natt
; i
++)
1075 mass_min
= std::min(mass_min
, att
[i
].prop
.mass
);
1080 fprintf(debug
, "Derivatives of non-bonded potentials at the cut-off:\n");
1081 fprintf(debug
, "LJ disp. -V' %9.2e V'' %9.2e -V''' %9.2e\n", ljDisp
.md1
, ljDisp
.d2
, ljDisp
.md3
);
1082 fprintf(debug
, "LJ rep. -V' %9.2e V'' %9.2e -V''' %9.2e\n", ljRep
.md1
, ljRep
.d2
, ljRep
.md3
);
1083 fprintf(debug
, "Electro. -V' %9.2e V'' %9.2e\n", elec
.md1
, elec
.d2
);
1084 fprintf(debug
, "sqrt(kT_fac) %f\n", std::sqrt(kT_fac
));
1085 fprintf(debug
, "mass_min %f\n", mass_min
);
1088 /* Search using bisection */
1090 /* The drift will be neglible at 5 times the max sigma */
1091 ib1
= (int)(5*2*std::sqrt(kT_fac
/mass_min
)/resolution
) + 1;
1092 while (ib1
- ib0
> 1)
1096 rl
= std::max(ir
->rvdw
, ir
->rcoulomb
) + rb
;
1098 /* Calculate the average energy drift at the last step
1099 * of the nstlist steps at which the pair-list is used.
1101 drift
= energyDrift(att
, natt
, &mtop
->ffparams
,
1103 &ljDisp
, &ljRep
, &elec
,
1104 ir
->rvdw
, ir
->rcoulomb
,
1107 /* Correct for the fact that we are using a Ni x Nj particle pair list
1108 * and not a 1 x 1 particle pair list. This reduces the drift.
1110 /* We don't have a formula for 8 (yet), use 4 which is conservative */
1111 nb_clust_frac_pairs_not_in_list_at_cutoff
=
1112 surface_frac(std::min(list_setup
->cluster_size_i
, 4),
1113 particle_distance
, rl
)*
1114 surface_frac(std::min(list_setup
->cluster_size_j
, 4),
1115 particle_distance
, rl
);
1116 drift
*= nb_clust_frac_pairs_not_in_list_at_cutoff
;
1118 /* Convert the drift to drift per unit time per atom */
1119 drift
/= ir
->nstlist
*ir
->delta_t
*mtop
->natoms
;
1123 fprintf(debug
, "ib %3d %3d %3d rb %.3f %dx%d fac %.3f drift %.1e\n",
1125 list_setup
->cluster_size_i
, list_setup
->cluster_size_j
,
1126 nb_clust_frac_pairs_not_in_list_at_cutoff
,
1130 if (std::abs(drift
) > ir
->verletbuf_tol
)
1142 *rlist
= std::max(ir
->rvdw
, ir
->rcoulomb
) + ib1
*resolution
;