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.
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 void verletbuf_get_list_setup(gmx_bool gmx_unused bSIMD
,
112 verletbuf_list_setup_t
*list_setup
)
114 /* When calling this function we often don't know which kernel type we
115 * are going to use. W choose the kernel type with the smallest possible
116 * i- and j-cluster sizes, so we potentially overestimate, but never
117 * underestimate, the buffer drift.
118 * Note that the current buffer estimation code only handles clusters
119 * of size 1, 2 or 4, so for 4x8 or 8x8 we use the estimate for 4x4.
124 /* The CUDA kernels split the j-clusters in two halves */
125 list_setup
->cluster_size_i
= nbnxn_kernel_to_cluster_i_size(nbnxnk8x8x8_GPU
);
126 list_setup
->cluster_size_j
= nbnxn_kernel_to_cluster_j_size(nbnxnk8x8x8_GPU
)/2;
132 kernel_type
= nbnxnk4x4_PlainC
;
137 #ifdef GMX_NBNXN_SIMD_2XNN
138 /* We use the smallest cluster size to be on the safe side */
139 kernel_type
= nbnxnk4xN_SIMD_2xNN
;
141 kernel_type
= nbnxnk4xN_SIMD_4xN
;
146 list_setup
->cluster_size_i
= nbnxn_kernel_to_cluster_i_size(kernel_type
);
147 list_setup
->cluster_size_j
= nbnxn_kernel_to_cluster_j_size(kernel_type
);
152 atom_nonbonded_kinetic_prop_equal(const atom_nonbonded_kinetic_prop_t
*prop1
,
153 const atom_nonbonded_kinetic_prop_t
*prop2
)
155 return (prop1
->mass
== prop2
->mass
&&
156 prop1
->type
== prop2
->type
&&
157 prop1
->q
== prop2
->q
&&
158 prop1
->bConstr
== prop2
->bConstr
&&
159 prop1
->con_mass
== prop2
->con_mass
&&
160 prop1
->con_len
== prop2
->con_len
);
163 static void add_at(verletbuf_atomtype_t
**att_p
, int *natt_p
,
164 const atom_nonbonded_kinetic_prop_t
*prop
,
167 verletbuf_atomtype_t
*att
;
172 /* Ignore massless particles */
180 while (i
< natt
&& !atom_nonbonded_kinetic_prop_equal(prop
, &att
[i
].prop
))
192 srenew(*att_p
, *natt_p
);
193 (*att_p
)[i
].prop
= *prop
;
194 (*att_p
)[i
].n
= nmol
;
198 static void get_vsite_masses(const gmx_moltype_t
*moltype
,
199 const gmx_ffparams_t
*ffparams
,
208 /* Check for virtual sites, determine mass from constructing atoms */
209 for (ft
= 0; ft
< F_NRE
; ft
++)
213 il
= &moltype
->ilist
[ft
];
215 for (i
= 0; i
< il
->nr
; i
+= 1+NRAL(ft
))
218 real inv_mass
, coeff
, m_aj
;
221 ip
= &ffparams
->iparams
[il
->iatoms
[i
]];
223 a1
= il
->iatoms
[i
+1];
227 /* Only vsiten can have more than four
228 constructing atoms, so NRAL(ft) <= 5 */
231 const int maxj
= NRAL(ft
);
235 for (j
= 1; j
< maxj
; j
++)
237 cam
[j
] = moltype
->atoms
.atom
[il
->iatoms
[i
+1+j
]].m
;
240 cam
[j
] = vsite_m
[il
->iatoms
[i
+1+j
]];
244 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.",
246 interaction_function
[ft
].longname
,
247 il
->iatoms
[i
+1+j
]+1);
255 vsite_m
[a1
] = (cam
[1]*cam
[2])/(cam
[2]*gmx::square(1-ip
->vsite
.a
) + cam
[1]*gmx::square(ip
->vsite
.a
));
259 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
));
262 gmx_incons("Invalid vsite type");
265 /* Use the mass of the lightest constructing atom.
266 * This is an approximation.
267 * If the distance of the virtual site to the
268 * constructing atom is less than all distances
269 * between constructing atoms, this is a safe
270 * over-estimate of the displacement of the vsite.
271 * This condition holds for all H mass replacement
272 * vsite constructions, except for SP2/3 groups.
273 * In SP3 groups one H will have a F_VSITE3
274 * construction, so even there the total drift
275 * estimate shouldn't be far off.
277 vsite_m
[a1
] = cam
[1];
278 for (j
= 2; j
< maxj
; j
++)
280 vsite_m
[a1
] = std::min(vsite_m
[a1
], cam
[j
]);
293 for (j
= 0; j
< 3*ffparams
->iparams
[il
->iatoms
[i
]].vsiten
.n
; j
+= 3)
295 aj
= il
->iatoms
[i
+j
+2];
296 coeff
= ffparams
->iparams
[il
->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 for loop increment of i */
313 i
+= j
- 1 - NRAL(ft
);
317 fprintf(debug
, "atom %4d %-20s mass %6.3f\n",
318 a1
, interaction_function
[ft
].longname
, vsite_m
[a1
]);
325 static void get_verlet_buffer_atomtypes(const gmx_mtop_t
*mtop
,
326 verletbuf_atomtype_t
**att_p
,
330 verletbuf_atomtype_t
*att
;
332 int mb
, nmol
, ft
, i
, a1
, a2
, a3
, a
;
333 const t_atoms
*atoms
;
336 atom_nonbonded_kinetic_prop_t
*prop
;
338 int n_nonlin_vsite_mol
;
343 if (n_nonlin_vsite
!= nullptr)
348 for (mb
= 0; mb
< mtop
->nmolblock
; mb
++)
350 nmol
= mtop
->molblock
[mb
].nmol
;
352 atoms
= &mtop
->moltype
[mtop
->molblock
[mb
].type
].atoms
;
354 /* Check for constraints, as they affect the kinetic energy.
355 * For virtual sites we need the masses and geometry of
356 * the constructing atoms to determine their velocity distribution.
358 snew(prop
, atoms
->nr
);
359 snew(vsite_m
, atoms
->nr
);
361 for (ft
= F_CONSTR
; ft
<= F_CONSTRNC
; ft
++)
363 il
= &mtop
->moltype
[mtop
->molblock
[mb
].type
].ilist
[ft
];
365 for (i
= 0; i
< il
->nr
; i
+= 1+NRAL(ft
))
367 ip
= &mtop
->ffparams
.iparams
[il
->iatoms
[i
]];
368 a1
= il
->iatoms
[i
+1];
369 a2
= il
->iatoms
[i
+2];
370 if (atoms
->atom
[a2
].m
> prop
[a1
].con_mass
)
372 prop
[a1
].con_mass
= atoms
->atom
[a2
].m
;
373 prop
[a1
].con_len
= ip
->constr
.dA
;
375 if (atoms
->atom
[a1
].m
> prop
[a2
].con_mass
)
377 prop
[a2
].con_mass
= atoms
->atom
[a1
].m
;
378 prop
[a2
].con_len
= ip
->constr
.dA
;
383 il
= &mtop
->moltype
[mtop
->molblock
[mb
].type
].ilist
[F_SETTLE
];
385 for (i
= 0; i
< il
->nr
; i
+= 1+NRAL(F_SETTLE
))
387 ip
= &mtop
->ffparams
.iparams
[il
->iatoms
[i
]];
388 a1
= il
->iatoms
[i
+1];
389 a2
= il
->iatoms
[i
+2];
390 a3
= il
->iatoms
[i
+3];
391 /* Usually the mass of a1 (usually oxygen) is larger than a2/a3.
392 * If this is not the case, we overestimate the displacement,
393 * which leads to a larger buffer (ok since this is an exotic case).
395 prop
[a1
].con_mass
= atoms
->atom
[a2
].m
;
396 prop
[a1
].con_len
= ip
->settle
.doh
;
398 prop
[a2
].con_mass
= atoms
->atom
[a1
].m
;
399 prop
[a2
].con_len
= ip
->settle
.doh
;
401 prop
[a3
].con_mass
= atoms
->atom
[a1
].m
;
402 prop
[a3
].con_len
= ip
->settle
.doh
;
405 get_vsite_masses(&mtop
->moltype
[mtop
->molblock
[mb
].type
],
408 &n_nonlin_vsite_mol
);
409 if (n_nonlin_vsite
!= nullptr)
411 *n_nonlin_vsite
+= nmol
*n_nonlin_vsite_mol
;
414 for (a
= 0; a
< atoms
->nr
; a
++)
416 if (atoms
->atom
[a
].ptype
== eptVSite
)
418 prop
[a
].mass
= vsite_m
[a
];
422 prop
[a
].mass
= atoms
->atom
[a
].m
;
424 prop
[a
].type
= atoms
->atom
[a
].type
;
425 prop
[a
].q
= atoms
->atom
[a
].q
;
426 /* We consider an atom constrained, #DOF=2, when it is
427 * connected with constraints to (at least one) atom with
428 * a mass of more than 0.4x its own mass. This is not a critical
429 * parameter, since with roughly equal masses the unconstrained
430 * and constrained displacement will not differ much (and both
431 * overestimate the displacement).
433 prop
[a
].bConstr
= (prop
[a
].con_mass
> 0.4*prop
[a
].mass
);
435 add_at(&att
, &natt
, &prop
[a
], nmol
);
444 for (a
= 0; a
< natt
; a
++)
446 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",
447 a
, att
[a
].prop
.mass
, att
[a
].prop
.type
, att
[a
].prop
.q
,
448 att
[a
].prop
.bConstr
, att
[a
].prop
.con_mass
, att
[a
].prop
.con_len
,
457 /* This function computes two components of the estimate of the variance
458 * in the displacement of one atom in a system of two constrained atoms.
459 * Returns in sigma2_2d the variance due to rotation of the constrained
460 * atom around the atom to which it constrained.
461 * Returns in sigma2_3d the variance due to displacement of the COM
462 * of the whole system of the two constrained atoms.
464 * Note that we only take a single constraint (the one to the heaviest atom)
465 * into account. If an atom has multiple constraints, this will result in
466 * an overestimate of the displacement, which gives a larger drift and buffer.
468 void constrained_atom_sigma2(real kT_fac
,
469 const atom_nonbonded_kinetic_prop_t
*prop
,
473 /* Here we decompose the motion of a constrained atom into two
474 * components: rotation around the COM and translation of the COM.
477 /* Determine the variance of the arc length for the two rotational DOFs */
478 real massFraction
= prop
->con_mass
/(prop
->mass
+ prop
->con_mass
);
479 real sigma2_rot
= kT_fac
*massFraction
/prop
->mass
;
481 /* The distance from the atom to the COM, i.e. the rotational arm */
482 real comDistance
= prop
->con_len
*massFraction
;
484 /* The variance relative to the arm */
485 real sigma2_rel
= sigma2_rot
/gmx::square(comDistance
);
487 /* For sigma2_rel << 1 we don't notice the rotational effect and
488 * we have a normal, Gaussian displacement distribution.
489 * For larger sigma2_rel the displacement is much less, in fact it can
490 * not exceed 2*comDistance. We can calculate MSD/arm^2 as:
491 * integral_x=0-inf distance2(x) x/sigma2_rel exp(-x^2/(2 sigma2_rel)) dx
492 * where x is angular displacement and distance2(x) is the distance^2
493 * between points at angle 0 and x:
494 * distance2(x) = (sin(x) - sin(0))^2 + (cos(x) - cos(0))^2
495 * The limiting value of this MSD is 2, which is also the value for
496 * a uniform rotation distribution that would be reached at long time.
497 * The maximum is 2.5695 at sigma2_rel = 4.5119.
498 * We approximate this integral with a rational polynomial with
499 * coefficients from a Taylor expansion. This approximation is an
500 * overestimate for all values of sigma2_rel. Its maximum value
501 * of 2.6491 is reached at sigma2_rel = sqrt(45/2) = 4.7434.
502 * We keep the approximation constant after that.
503 * We use this approximate MSD as the variance for a Gaussian distribution.
505 * NOTE: For any sensible buffer tolerance this will result in a (large)
506 * overestimate of the buffer size, since the Gaussian has a long tail,
507 * whereas the actual distribution can not reach values larger than 2.
509 /* Coeffients obtained from a Taylor expansion */
510 const real a
= 1.0/3.0;
511 const real b
= 2.0/45.0;
513 /* Our approximation is constant after sigma2_rel = 1/sqrt(b) */
514 sigma2_rel
= std::min(sigma2_rel
, 1/std::sqrt(b
));
516 /* Compute the approximate sigma^2 for 2D motion due to the rotation */
517 *sigma2_2d
= gmx::square(comDistance
)*
518 sigma2_rel
/(1 + a
*sigma2_rel
+ b
*gmx::square(sigma2_rel
));
520 /* The constrained atom also moves (in 3D) with the COM of both atoms */
521 *sigma2_3d
= kT_fac
/(prop
->mass
+ prop
->con_mass
);
524 static void get_atom_sigma2(real kT_fac
,
525 const atom_nonbonded_kinetic_prop_t
*prop
,
531 /* Complicated constraint calculation in a separate function */
532 constrained_atom_sigma2(kT_fac
, prop
, sigma2_2d
, sigma2_3d
);
536 /* Unconstrained atom: trivial */
538 *sigma2_3d
= kT_fac
/prop
->mass
;
542 static void approx_2dof(real s2
, real x
, real
*shift
, real
*scale
)
544 /* A particle with 1 DOF constrained has 2 DOFs instead of 3.
545 * This code is also used for particles with multiple constraints,
546 * in which case we overestimate the displacement.
547 * The 2DOF distribution is sqrt(pi/2)*erfc(r/(sqrt(2)*s))/(2*s).
548 * We approximate this with scale*Gaussian(s,r+shift),
549 * by matching the distribution value and derivative at x.
550 * This is a tight overestimate for all r>=0 at any s and x.
554 ex
= std::exp(-x
*x
/(2*s2
));
555 er
= std::erfc(x
/std::sqrt(2*s2
));
557 *shift
= -x
+ std::sqrt(2*s2
/M_PI
)*ex
/er
;
558 *scale
= 0.5*M_PI
*std::exp(ex
*ex
/(M_PI
*er
*er
))*er
;
561 // Returns an (over)estimate of the energy drift for a single atom pair,
562 // given the kinetic properties, displacement variances and list buffer.
563 static real
energyDriftAtomPair(const atom_nonbonded_kinetic_prop_t
*prop_i
,
564 const atom_nonbonded_kinetic_prop_t
*prop_j
,
565 real s2
, real s2i_2d
, real s2j_2d
,
567 const pot_derivatives_t
*der
)
569 // For relatively small arguments erfc() is so small that if will be 0.0
570 // when stored in a float. We set an argument limit of 8 (Erfc(8)=1e-29),
571 // such that we can divide by erfc and have some space left for arithmetic.
572 const real erfc_arg_max
= 8.0;
579 if (rsh
*rsh
> 2*s2
*erfc_arg_max
*erfc_arg_max
)
581 // Below we calculate c_erfc = 0.5*erfc(rsh/sqrt(2*s2))
582 // When rsh/sqrt(2*s2) increases, this erfc will be the first
583 // result that underflows and becomes 0.0. To avoid this,
584 // we set c_exp=0 and c_erfc=0 for large arguments.
585 // This also avoids NaN in approx_2dof().
586 // In any relevant case this has no effect on the results,
587 // since c_exp < 6e-29, so the displacement is completely
588 // negligible for such atom pairs (and an overestimate).
589 // In nearly all use cases, there will be other atom pairs
590 // that contribute much more to the total, so zeroing
591 // this particular contribution has no effect at all.
597 /* For constraints: adapt r and scaling for the Gaussian */
602 approx_2dof(s2i_2d
, r_buffer
*s2i_2d
/s2
, &sh
, &sc
);
610 approx_2dof(s2j_2d
, r_buffer
*s2j_2d
/s2
, &sh
, &sc
);
615 /* Exact contribution of an atom pair with Gaussian displacement
616 * with sigma s to the energy drift for a potential with
617 * derivative -md and second derivative dd at the cut-off.
618 * The only catch is that for potentials that change sign
619 * near the cut-off there could be an unlucky compensation
620 * of positive and negative energy drift.
621 * Such potentials are extremely rare though.
623 * Note that pot has unit energy*length, as the linear
624 * atom density still needs to be put in.
626 c_exp
= std::exp(-rsh
*rsh
/(2*s2
))/std::sqrt(2*M_PI
);
627 c_erfc
= 0.5*std::erfc(rsh
/(std::sqrt(2*s2
)));
629 real s
= std::sqrt(s2
);
633 der
->md1
/2*((rsh2
+ s2
)*c_erfc
- rsh
*s
*c_exp
);
635 der
->d2
/6*(s
*(rsh2
+ 2*s2
)*c_exp
- rsh
*(rsh2
+ 3*s2
)*c_erfc
);
637 der
->md3
/24*((rsh2
*rsh2
+ 6*rsh2
*s2
+ 3*s2
*s2
)*c_erfc
- rsh
*s
*(rsh2
+ 5*s2
)*c_exp
);
639 return pot1
+ pot2
+ pot3
;
642 static real
energyDrift(const verletbuf_atomtype_t
*att
, int natt
,
643 const gmx_ffparams_t
*ffp
,
645 const pot_derivatives_t
*ljDisp
,
646 const pot_derivatives_t
*ljRep
,
647 const pot_derivatives_t
*elec
,
648 real rlj
, real rcoulomb
,
649 real rlist
, real boxvol
)
651 double drift_tot
= 0;
655 /* No atom displacements: no drift, avoid division by 0 */
659 // Here add up the contribution of all atom pairs in the system to
660 // (estimated) energy drift by looping over all atom type pairs.
661 for (int i
= 0; i
< natt
; i
++)
663 // Get the thermal displacement variance for the i-atom type
664 const atom_nonbonded_kinetic_prop_t
*prop_i
= &att
[i
].prop
;
666 get_atom_sigma2(kT_fac
, prop_i
, &s2i_2d
, &s2i_3d
);
668 for (int j
= i
; j
< natt
; j
++)
670 // Get the thermal displacement variance for the j-atom type
671 const atom_nonbonded_kinetic_prop_t
*prop_j
= &att
[j
].prop
;
673 get_atom_sigma2(kT_fac
, prop_j
, &s2j_2d
, &s2j_3d
);
675 /* Add up the up to four independent variances */
676 real s2
= s2i_2d
+ s2i_3d
+ s2j_2d
+ s2j_3d
;
678 // Set -V', V'' and -V''' at the cut-off for LJ */
679 real c6
= ffp
->iparams
[prop_i
->type
*ffp
->atnr
+ prop_j
->type
].lj
.c6
;
680 real c12
= ffp
->iparams
[prop_i
->type
*ffp
->atnr
+ prop_j
->type
].lj
.c12
;
681 pot_derivatives_t lj
;
682 lj
.md1
= c6
*ljDisp
->md1
+ c12
*ljRep
->md1
;
683 lj
.d2
= c6
*ljDisp
->d2
+ c12
*ljRep
->d2
;
684 lj
.md3
= c6
*ljDisp
->md3
+ c12
*ljRep
->md3
;
686 real pot_lj
= energyDriftAtomPair(prop_i
, prop_j
,
691 // Set -V' and V'' at the cut-off for Coulomb
692 pot_derivatives_t elec_qq
;
693 elec_qq
.md1
= elec
->md1
*prop_i
->q
*prop_j
->q
;
694 elec_qq
.d2
= elec
->d2
*prop_i
->q
*prop_j
->q
;
697 real pot_q
= energyDriftAtomPair(prop_i
, prop_j
,
702 // Note that attractive and repulsive potentials for individual
703 // pairs can partially cancel.
704 real pot
= pot_lj
+ pot_q
;
706 /* Multiply by the number of atom pairs */
709 pot
*= (double)att
[i
].n
*(att
[i
].n
- 1)/2;
713 pot
*= (double)att
[i
].n
*att
[j
].n
;
715 /* We need the line density to get the energy drift of the system.
716 * The effective average r^2 is close to (rlist+sigma)^2.
718 pot
*= 4*M_PI
*gmx::square(rlist
+ std::sqrt(s2
))/boxvol
;
720 /* Add the unsigned drift to avoid cancellation of errors */
721 drift_tot
+= std::abs(pot
);
728 static real
surface_frac(int cluster_size
, real particle_distance
, real rlist
)
732 if (rlist
< 0.5*particle_distance
)
734 /* We have non overlapping spheres */
738 /* Half the inter-particle distance relative to rlist */
739 d
= 0.5*particle_distance
/rlist
;
741 /* Determine the area of the surface at distance rlist to the closest
742 * particle, relative to surface of a sphere of radius rlist.
743 * The formulas below assume close to cubic cells for the pair search grid,
744 * which the pair search code tries to achieve.
745 * Note that in practice particle distances will not be delta distributed,
746 * but have some spread, often involving shorter distances,
747 * as e.g. O-H bonds in a water molecule. Thus the estimates below will
748 * usually be slightly too high and thus conservative.
750 switch (cluster_size
)
753 /* One particle: trivial */
757 /* Two particles: two spheres at fractional distance 2*a */
761 /* We assume a perfect, symmetric tetrahedron geometry.
762 * The surface around a tetrahedron is too complex for a full
763 * analytical solution, so we use a Taylor expansion.
765 area_rel
= (1.0 + 1/M_PI
*(6*std::acos(1/std::sqrt(3))*d
+
766 std::sqrt(3)*d
*d
*(1.0 +
769 83.0/756.0*d
*d
*d
*d
*d
*d
)));
772 gmx_incons("surface_frac called with unsupported cluster_size");
776 return area_rel
/cluster_size
;
779 /* Returns the negative of the third derivative of a potential r^-p
780 * with a force-switch function, evaluated at the cut-off rc.
782 static real
md3_force_switch(real p
, real rswitch
, real rc
)
784 /* The switched force function is:
785 * p*r^-(p+1) + a*(r - rswitch)^2 + b*(r - rswitch)^3
788 real md3_pot
, md3_sw
;
790 a
= -((p
+ 4)*rc
- (p
+ 1)*rswitch
)/(pow(rc
, p
+2)*gmx::square(rc
-rswitch
));
791 b
= ((p
+ 3)*rc
- (p
+ 1)*rswitch
)/(pow(rc
, p
+2)*gmx::power3(rc
-rswitch
));
793 md3_pot
= (p
+ 2)*(p
+ 1)*p
*pow(rc
, p
+3);
794 md3_sw
= 2*a
+ 6*b
*(rc
- rswitch
);
796 return md3_pot
+ md3_sw
;
799 void calc_verlet_buffer_size(const gmx_mtop_t
*mtop
, real boxvol
,
800 const t_inputrec
*ir
,
801 real reference_temperature
,
802 const verletbuf_list_setup_t
*list_setup
,
809 real particle_distance
;
810 real nb_clust_frac_pairs_not_in_list_at_cutoff
;
812 verletbuf_atomtype_t
*att
= nullptr;
815 real kT_fac
, mass_min
;
820 if (!EI_DYNAMICS(ir
->eI
))
822 gmx_incons("Can only determine the Verlet buffer size for integrators that perform dynamics");
824 if (ir
->verletbuf_tol
<= 0)
826 gmx_incons("The Verlet buffer tolerance needs to be larger than zero");
829 if (reference_temperature
< 0)
831 if (EI_MD(ir
->eI
) && ir
->etc
== etcNO
)
833 /* This case should be handled outside calc_verlet_buffer_size */
834 gmx_incons("calc_verlet_buffer_size called with an NVE ensemble and reference_temperature < 0");
837 /* We use the maximum temperature with multiple T-coupl groups.
838 * We could use a per particle temperature, but since particles
839 * interact, this might underestimate the buffer size.
841 reference_temperature
= 0;
842 for (i
= 0; i
< ir
->opts
.ngtc
; i
++)
844 if (ir
->opts
.tau_t
[i
] >= 0)
846 reference_temperature
= std::max(reference_temperature
,
852 /* Resolution of the buffer size */
855 env
= getenv("GMX_VERLET_BUFFER_RES");
858 sscanf(env
, "%lf", &resolution
);
861 /* In an atom wise pair-list there would be no pairs in the list
862 * beyond the pair-list cut-off.
863 * However, we use a pair-list of groups vs groups of atoms.
864 * For groups of 4 atoms, the parallelism of SSE instructions, only
865 * 10% of the atoms pairs are not in the list just beyond the cut-off.
866 * As this percentage increases slowly compared to the decrease of the
867 * Gaussian displacement distribution over this range, we can simply
868 * reduce the drift by this fraction.
869 * For larger groups, e.g. of 8 atoms, this fraction will be lower,
870 * so then buffer size will be on the conservative (large) side.
872 * Note that the formulas used here do not take into account
873 * cancellation of errors which could occur by missing both
874 * attractive and repulsive interactions.
876 * The only major assumption is homogeneous particle distribution.
877 * For an inhomogeneous system, such as a liquid-vapor system,
878 * the buffer will be underestimated. The actual energy drift
879 * will be higher by the factor: local/homogeneous particle density.
881 * The results of this estimate have been checked againt simulations.
882 * In most cases the real drift differs by less than a factor 2.
885 /* Worst case assumption: HCP packing of particles gives largest distance */
886 particle_distance
= std::cbrt(boxvol
*std::sqrt(2)/mtop
->natoms
);
888 get_verlet_buffer_atomtypes(mtop
, &att
, &natt
, n_nonlin_vsite
);
889 assert(att
!= NULL
&& natt
>= 0);
893 fprintf(debug
, "particle distance assuming HCP packing: %f nm\n",
895 fprintf(debug
, "energy drift atom types: %d\n", natt
);
898 pot_derivatives_t ljDisp
= { 0, 0, 0 };
899 pot_derivatives_t ljRep
= { 0, 0, 0 };
900 real repPow
= mtop
->ffparams
.reppow
;
902 if (ir
->vdwtype
== evdwCUT
)
904 real sw_range
, md3_pswf
;
906 switch (ir
->vdw_modifier
)
909 case eintmodPOTSHIFT
:
910 /* -dV/dr of -r^-6 and r^-reppow */
911 ljDisp
.md1
= -6*std::pow(ir
->rvdw
, -7.0);
912 ljRep
.md1
= repPow
*std::pow(ir
->rvdw
, -(repPow
+ 1));
913 /* The contribution of the higher derivatives is negligible */
915 case eintmodFORCESWITCH
:
916 /* At the cut-off: V=V'=V''=0, so we use only V''' */
917 ljDisp
.md3
= -md3_force_switch(6.0, ir
->rvdw_switch
, ir
->rvdw
);
918 ljRep
.md3
= md3_force_switch(repPow
, ir
->rvdw_switch
, ir
->rvdw
);
920 case eintmodPOTSWITCH
:
921 /* At the cut-off: V=V'=V''=0.
922 * V''' is given by the original potential times
923 * the third derivative of the switch function.
925 sw_range
= ir
->rvdw
- ir
->rvdw_switch
;
926 md3_pswf
= 60.0/gmx::power3(sw_range
);
928 ljDisp
.md3
= -std::pow(ir
->rvdw
, -6.0 )*md3_pswf
;
929 ljRep
.md3
= std::pow(ir
->rvdw
, -repPow
)*md3_pswf
;
932 gmx_incons("Unimplemented VdW modifier");
935 else if (EVDW_PME(ir
->vdwtype
))
937 real b
= calc_ewaldcoeff_lj(ir
->rvdw
, ir
->ewald_rtol_lj
);
943 // -dV/dr of g(br)*r^-6 [where g(x) = exp(-x^2)(1+x^2+x^4/2),
944 // see LJ-PME equations in manual] and r^-reppow
945 ljDisp
.md1
= -std::exp(-br2
)*(br6
+ 3.0*br4
+ 6.0*br2
+ 6.0)*std::pow(r
, -7.0);
946 ljRep
.md1
= repPow
*pow(r
, -(repPow
+ 1));
947 // The contribution of the higher derivatives is negligible
951 gmx_fatal(FARGS
, "Energy drift calculation is only implemented for plain cut-off Lennard-Jones interactions");
954 elfac
= ONE_4PI_EPS0
/ir
->epsilon_r
;
956 // Determine the 1st and 2nd derivative for the electostatics
957 pot_derivatives_t elec
= { 0, 0, 0 };
959 if (ir
->coulombtype
== eelCUT
|| EEL_RF(ir
->coulombtype
))
963 if (ir
->coulombtype
== eelCUT
)
970 eps_rf
= ir
->epsilon_rf
/ir
->epsilon_r
;
973 k_rf
= (eps_rf
- ir
->epsilon_r
)/( gmx::power3(ir
->rcoulomb
) * (2*eps_rf
+ ir
->epsilon_r
) );
977 /* epsilon_rf = infinity */
978 k_rf
= 0.5/gmx::power3(ir
->rcoulomb
);
984 elec
.md1
= elfac
*(1.0/gmx::square(ir
->rcoulomb
) - 2*k_rf
*ir
->rcoulomb
);
986 elec
.d2
= elfac
*(2.0/gmx::power3(ir
->rcoulomb
) + 2*k_rf
);
988 else if (EEL_PME(ir
->coulombtype
) || ir
->coulombtype
== eelEWALD
)
992 b
= calc_ewaldcoeff_q(ir
->rcoulomb
, ir
->ewald_rtol
);
995 elec
.md1
= elfac
*(b
*std::exp(-br
*br
)*M_2_SQRTPI
/rc
+ std::erfc(br
)/(rc
*rc
));
996 elec
.d2
= elfac
/(rc
*rc
)*(2*b
*(1 + br
*br
)*std::exp(-br
*br
)*M_2_SQRTPI
+ 2*std::erfc(br
)/rc
);
1000 gmx_fatal(FARGS
, "Energy drift calculation is only implemented for Reaction-Field and Ewald electrostatics");
1003 /* Determine the variance of the atomic displacement
1004 * over nstlist-1 steps: kT_fac
1005 * For inertial dynamics (not Brownian dynamics) the mass factor
1006 * is not included in kT_fac, it is added later.
1010 /* Get the displacement distribution from the random component only.
1011 * With accurate integration the systematic (force) displacement
1012 * should be negligible (unless nstlist is extremely large, which
1013 * you wouldn't do anyhow).
1015 kT_fac
= 2*BOLTZ
*reference_temperature
*(ir
->nstlist
-1)*ir
->delta_t
;
1016 if (ir
->bd_fric
> 0)
1018 /* This is directly sigma^2 of the displacement */
1019 kT_fac
/= ir
->bd_fric
;
1021 /* Set the masses to 1 as kT_fac is the full sigma^2,
1022 * but we divide by m in ener_drift().
1024 for (i
= 0; i
< natt
; i
++)
1026 att
[i
].prop
.mass
= 1;
1033 /* Per group tau_t is not implemented yet, use the maximum */
1034 tau_t
= ir
->opts
.tau_t
[0];
1035 for (i
= 1; i
< ir
->opts
.ngtc
; i
++)
1037 tau_t
= std::max(tau_t
, ir
->opts
.tau_t
[i
]);
1041 /* This kT_fac needs to be divided by the mass to get sigma^2 */
1046 kT_fac
= BOLTZ
*reference_temperature
*gmx::square((ir
->nstlist
-1)*ir
->delta_t
);
1049 mass_min
= att
[0].prop
.mass
;
1050 for (i
= 1; i
< natt
; i
++)
1052 mass_min
= std::min(mass_min
, att
[i
].prop
.mass
);
1057 fprintf(debug
, "Derivatives of non-bonded potentials at the cut-off:\n");
1058 fprintf(debug
, "LJ disp. -V' %9.2e V'' %9.2e -V''' %9.2e\n", ljDisp
.md1
, ljDisp
.d2
, ljDisp
.md3
);
1059 fprintf(debug
, "LJ rep. -V' %9.2e V'' %9.2e -V''' %9.2e\n", ljRep
.md1
, ljRep
.d2
, ljRep
.md3
);
1060 fprintf(debug
, "Electro. -V' %9.2e V'' %9.2e\n", elec
.md1
, elec
.d2
);
1061 fprintf(debug
, "sqrt(kT_fac) %f\n", std::sqrt(kT_fac
));
1062 fprintf(debug
, "mass_min %f\n", mass_min
);
1065 /* Search using bisection */
1067 /* The drift will be neglible at 5 times the max sigma */
1068 ib1
= (int)(5*2*std::sqrt(kT_fac
/mass_min
)/resolution
) + 1;
1069 while (ib1
- ib0
> 1)
1073 rl
= std::max(ir
->rvdw
, ir
->rcoulomb
) + rb
;
1075 /* Calculate the average energy drift at the last step
1076 * of the nstlist steps at which the pair-list is used.
1078 drift
= energyDrift(att
, natt
, &mtop
->ffparams
,
1080 &ljDisp
, &ljRep
, &elec
,
1081 ir
->rvdw
, ir
->rcoulomb
,
1084 /* Correct for the fact that we are using a Ni x Nj particle pair list
1085 * and not a 1 x 1 particle pair list. This reduces the drift.
1087 /* We don't have a formula for 8 (yet), use 4 which is conservative */
1088 nb_clust_frac_pairs_not_in_list_at_cutoff
=
1089 surface_frac(std::min(list_setup
->cluster_size_i
, 4),
1090 particle_distance
, rl
)*
1091 surface_frac(std::min(list_setup
->cluster_size_j
, 4),
1092 particle_distance
, rl
);
1093 drift
*= nb_clust_frac_pairs_not_in_list_at_cutoff
;
1095 /* Convert the drift to drift per unit time per atom */
1096 drift
/= ir
->nstlist
*ir
->delta_t
*mtop
->natoms
;
1100 fprintf(debug
, "ib %3d %3d %3d rb %.3f %dx%d fac %.3f drift %.1e\n",
1102 list_setup
->cluster_size_i
, list_setup
->cluster_size_j
,
1103 nb_clust_frac_pairs_not_in_list_at_cutoff
,
1107 if (std::abs(drift
) > ir
->verletbuf_tol
)
1119 *rlist
= std::max(ir
->rvdw
, ir
->rcoulomb
) + ib1
*resolution
;