1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * For more info, check our website at http://www.gromacs.org
34 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
42 #include <sys/types.h>
47 #include "gmx_fatal.h"
51 #include "calc_verletbuf.h"
52 #include "../mdlib/nbnxn_consts.h"
54 /* Struct for unique atom type for calculating the energy drift.
55 * The atom displacement depends on mass and constraints.
56 * The energy jump for given distance depend on LJ type and q.
61 int type
; /* type (used for LJ parameters) */
63 int con
; /* constrained: 0, else 1, if 1, use #DOF=2 iso 3 */
64 int n
; /* total #atoms of this type in the system */
65 } verletbuf_atomtype_t
;
68 void verletbuf_get_list_setup(gmx_bool bGPU
,
69 verletbuf_list_setup_t
*list_setup
)
71 list_setup
->cluster_size_i
= NBNXN_CPU_CLUSTER_I_SIZE
;
75 list_setup
->cluster_size_j
= NBNXN_GPU_CLUSTER_SIZE
;
80 list_setup
->cluster_size_j
= NBNXN_CPU_CLUSTER_I_SIZE
;
84 #ifdef GMX_X86_AVX_256
89 list_setup
->cluster_size_j
= simd_width
/(sizeof(real
)*8);
94 static void add_at(verletbuf_atomtype_t
**att_p
,int *natt_p
,
95 real mass
,int type
,real q
,int con
,int nmol
)
97 verletbuf_atomtype_t
*att
;
102 /* Ignore massless particles */
111 !(mass
== att
[i
].mass
&&
112 type
== att
[i
].type
&&
126 srenew(*att_p
,*natt_p
);
127 (*att_p
)[i
].mass
= mass
;
128 (*att_p
)[i
].type
= type
;
130 (*att_p
)[i
].con
= con
;
131 (*att_p
)[i
].n
= nmol
;
135 static void get_verlet_buffer_atomtypes(const gmx_mtop_t
*mtop
,
136 verletbuf_atomtype_t
**att_p
,
140 verletbuf_atomtype_t
*att
;
142 int mb
,nmol
,ft
,i
,j
,a1
,a2
,a3
,a
;
143 const t_atoms
*atoms
;
147 real
*con_m
,*vsite_m
,cam
[5];
152 if (n_nonlin_vsite
!= NULL
)
157 for(mb
=0; mb
<mtop
->nmolblock
; mb
++)
159 nmol
= mtop
->molblock
[mb
].nmol
;
161 atoms
= &mtop
->moltype
[mtop
->molblock
[mb
].type
].atoms
;
163 /* Check for constraints, as they affect the kinetic energy */
164 snew(con_m
,atoms
->nr
);
165 snew(vsite_m
,atoms
->nr
);
167 for(ft
=F_CONSTR
; ft
<=F_CONSTRNC
; ft
++)
169 il
= &mtop
->moltype
[mtop
->molblock
[mb
].type
].ilist
[ft
];
171 for(i
=0; i
<il
->nr
; i
+=1+NRAL(ft
))
173 a1
= il
->iatoms
[i
+1];
174 a2
= il
->iatoms
[i
+2];
175 con_m
[a1
] += atoms
->atom
[a2
].m
;
176 con_m
[a2
] += atoms
->atom
[a1
].m
;
180 il
= &mtop
->moltype
[mtop
->molblock
[mb
].type
].ilist
[F_SETTLE
];
182 for(i
=0; i
<il
->nr
; i
+=1+NRAL(F_SETTLE
))
184 a1
= il
->iatoms
[i
+1];
185 a2
= il
->iatoms
[i
+2];
186 a3
= il
->iatoms
[i
+3];
187 con_m
[a1
] += atoms
->atom
[a2
].m
+ atoms
->atom
[a3
].m
;
188 con_m
[a2
] += atoms
->atom
[a1
].m
+ atoms
->atom
[a3
].m
;
189 con_m
[a3
] += atoms
->atom
[a1
].m
+ atoms
->atom
[a2
].m
;
192 /* Check for virtual sites, determine mass from constructing atoms */
193 for(ft
=0; ft
<F_NRE
; ft
++)
197 il
= &mtop
->moltype
[mtop
->molblock
[mb
].type
].ilist
[ft
];
199 for(i
=0; i
<il
->nr
; i
+=1+NRAL(ft
))
201 ip
= &mtop
->ffparams
.iparams
[il
->iatoms
[i
]];
203 a1
= il
->iatoms
[i
+1];
205 for(j
=1; j
<NRAL(ft
); j
++)
207 cam
[j
] = atoms
->atom
[il
->iatoms
[i
+1+j
]].m
;
210 cam
[j
] = vsite_m
[il
->iatoms
[i
+1+j
]];
214 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.",
215 *mtop
->moltype
[mtop
->molblock
[mb
].type
].name
,
216 interaction_function
[ft
].longname
,
217 il
->iatoms
[i
+1+j
]+1);
224 /* Exact except for ignoring constraints */
225 vsite_m
[a1
] = (cam
[2]*sqr(1-ip
->vsite
.a
) + cam
[1]*sqr(ip
->vsite
.a
))/(cam
[1]*cam
[2]);
228 /* Exact except for ignoring constraints */
229 vsite_m
[a1
] = (cam
[2]*cam
[3]*sqr(1-ip
->vsite
.a
-ip
->vsite
.b
) + cam
[1]*cam
[3]*sqr(ip
->vsite
.a
) + cam
[1]*cam
[2]*sqr(ip
->vsite
.b
))/(cam
[1]*cam
[2]*cam
[3]);
232 /* Use the mass of the lightest constructing atom.
233 * This is an approximation.
234 * If the distance of the virtual site to the
235 * constructing atom is less than all distances
236 * between constructing atoms, this is a safe
237 * over-estimate of the displacement of the vsite.
238 * This condition holds for all H mass replacement
239 * replacement vsite constructions, except for SP2/3
240 * groups. In SP3 groups one H will have a F_VSITE3
241 * construction, so even there the total drift
242 * estimation shouldn't be far off.
245 vsite_m
[a1
] = cam
[1];
246 for(j
=2; j
<NRAL(ft
); j
++)
248 vsite_m
[a1
] = min(vsite_m
[a1
],cam
[j
]);
250 if (n_nonlin_vsite
!= NULL
)
252 *n_nonlin_vsite
+= nmol
;
260 for(a
=0; a
<atoms
->nr
; a
++)
262 at
= &atoms
->atom
[a
];
263 /* We consider an atom constrained, #DOF=2, when it is
264 * connected with constraints to one or more atoms with
265 * total mass larger than 1.5 that of the atom itself.
268 at
->m
,at
->type
,at
->q
,con_m
[a
] > 1.5*at
->m
,nmol
);
277 for(a
=0; a
<natt
; a
++)
279 fprintf(debug
,"type %d: m %5.2f t %d q %6.3f con %d n %d\n",
280 a
,att
[a
].mass
,att
[a
].type
,att
[a
].q
,att
[a
].con
,att
[a
].n
);
288 static void approx_2dof(real s2
,real x
,
289 real
*shift
,real
*scale
)
291 /* A particle with 1 DOF constrained has 2 DOFs instead of 3.
292 * This code is also used for particles with multiple constraints,
293 * in which case we overestimate the displacement.
294 * The 2DOF distribution is sqrt(pi/2)*erfc(r/(sqrt(2)*s))/(2*s).
295 * We approximate this with scale*Gaussian(s,r+shift),
296 * by matching the distribution value and derivative at x.
297 * This is a tight overestimate for all r>=0 at any s and x.
301 ex
= exp(-x
*x
/(2*s2
));
302 er
= gmx_erfc(x
/sqrt(2*s2
));
304 *shift
= -x
+ sqrt(2*s2
/M_PI
)*ex
/er
;
305 *scale
= 0.5*M_PI
*exp(ex
*ex
/(M_PI
*er
*er
))*er
;
308 static real
ener_drift(const verletbuf_atomtype_t
*att
,int natt
,
309 const gmx_ffparams_t
*ffp
,
311 real md_ljd
,real md_ljr
,real md_el
,real dd_el
,
313 real rlist
,real boxvol
)
315 double drift_tot
,pot1
,pot2
,pot
;
325 /* Loop over the different atom type pairs */
326 for(i
=0; i
<natt
; i
++)
328 s2i
= kT_fac
/att
[i
].mass
;
331 for(j
=i
; j
<natt
; j
++)
333 s2j
= kT_fac
/att
[j
].mass
;
336 /* Note that attractive and repulsive potentials for individual
337 * pairs will partially cancel.
339 /* -dV/dr at the cut-off for LJ + Coulomb */
341 md_ljd
*ffp
->iparams
[ti
*ffp
->atnr
+tj
].lj
.c6
+
342 md_ljr
*ffp
->iparams
[ti
*ffp
->atnr
+tj
].lj
.c12
+
343 md_el
*att
[i
].q
*att
[j
].q
;
345 /* d2V/dr2 at the cut-off for Coulomb, we neglect LJ */
346 dd
= dd_el
*att
[i
].q
*att
[j
].q
;
352 /* For constraints: adapt r and scaling for the Gaussian */
356 approx_2dof(s2i
,r_buffer
*s2i
/s2
,&sh
,&sc
);
363 approx_2dof(s2j
,r_buffer
*s2j
/s2
,&sh
,&sc
);
368 /* Exact contribution of an atom pair with Gaussian displacement
369 * with sigma s to the energy drift for a potential with
370 * derivative -md and second derivative dd at the cut-off.
371 * The only catch is that for potentials that change sign
372 * near the cut-off there could be an unlucky compensation
373 * of positive and negative energy drift.
374 * Such potentials are extremely rare though.
376 * Note that pot has unit energy*length, as the linear
377 * atom density still needs to be put in.
379 c_exp
= exp(-rsh
*rsh
/(2*s2
))/sqrt(2*M_PI
);
380 c_erfc
= 0.5*gmx_erfc(rsh
/(sqrt(2*s2
)));
384 md
/2*((rsh
*rsh
+ s2
)*c_erfc
- rsh
*s
*c_exp
);
386 dd
/6*(s
*(rsh
*rsh
+ 2*s2
)*c_exp
- rsh
*(rsh
*rsh
+ 3*s2
)*c_erfc
);
391 fprintf(debug
,"n %d %d d s %.3f %.3f con %d md %8.1e dd %8.1e pot1 %8.1e pot2 %8.1e pot %8.1e\n",
392 att
[i
].n
,att
[j
].n
,sqrt(s2i
),sqrt(s2j
),
393 att
[i
].con
+att
[j
].con
,
394 md
,dd
,pot1
,pot2
,pot
);
397 /* Multiply by the number of atom pairs */
400 pot
*= (double)att
[i
].n
*(att
[i
].n
- 1)/2;
404 pot
*= (double)att
[i
].n
*att
[j
].n
;
406 /* We need the line density to get the energy drift of the system.
407 * The effective average r^2 is close to (rlist+sigma)^2.
409 pot
*= 4*M_PI
*sqr(rlist
+ s
)/boxvol
;
411 /* Add the unsigned drift to avoid cancellation of errors */
412 drift_tot
+= fabs(pot
);
419 static real
surface_frac(int cluster_size
,real particle_distance
,real rlist
)
423 if (rlist
< 0.5*particle_distance
)
425 /* We have non overlapping spheres */
429 /* Half the inter-particle distance relative to rlist */
430 d
= 0.5*particle_distance
/rlist
;
432 /* Determine the area of the surface at distance rlist to the closest
433 * particle, relative to surface of a sphere of radius rlist.
434 * The formulas below assume close to cubic cells for the pair search grid,
435 * which the pair search code tries to achieve.
436 * Note that in practice particle distances will not be delta distributed,
437 * but have some spread, often involving shorter distances,
438 * as e.g. O-H bonds in a water molecule. Thus the estimates below will
439 * usually be slightly too high and thus conservative.
441 switch (cluster_size
)
444 /* One particle: trivial */
448 /* Two particles: two spheres at fractional distance 2*a */
452 /* We assume a perfect, symmetric tetrahedron geometry.
453 * The surface around a tetrahedron is too complex for a full
454 * analytical solution, so we use a Taylor expansion.
456 area_rel
= (1.0 + 1/M_PI
*(6*acos(1/sqrt(3))*d
+
460 83.0/756.0*d
*d
*d
*d
*d
*d
)));
463 gmx_incons("surface_frac called with unsupported cluster_size");
467 return area_rel
/cluster_size
;
470 void calc_verlet_buffer_size(const gmx_mtop_t
*mtop
,real boxvol
,
471 const t_inputrec
*ir
,real drift_target
,
472 const verletbuf_list_setup_t
*list_setup
,
479 real particle_distance
;
480 real nb_clust_frac_pairs_not_in_list_at_cutoff
;
482 verletbuf_atomtype_t
*att
=NULL
;
485 real md_ljd
,md_ljr
,md_el
,dd_el
;
487 real kT_fac
,mass_min
;
492 /* Resolution of the buffer size */
495 env
= getenv("GMX_VERLET_BUFFER_RES");
498 sscanf(env
,"%lf",&resolution
);
501 /* In an atom wise pair-list there would be no pairs in the list
502 * beyond the pair-list cut-off.
503 * However, we use a pair-list of groups vs groups of atoms.
504 * For groups of 4 atoms, the parallelism of SSE instructions, only
505 * 10% of the atoms pairs are not in the list just beyond the cut-off.
506 * As this percentage increases slowly compared to the decrease of the
507 * Gaussian displacement distribution over this range, we can simply
508 * reduce the drift by this fraction.
509 * For larger groups, e.g. of 8 atoms, this fraction will be lower,
510 * so then buffer size will be on the conservative (large) side.
512 * Note that the formulas used here do not take into account
513 * cancellation of errors which could occur by missing both
514 * attractive and repulsive interactions.
516 * The only major assumption is homogeneous particle distribution.
517 * For an inhomogeneous system, such as a liquid-vapor system,
518 * the buffer will be underestimated. The actual energy drift
519 * will be higher by the factor: local/homogeneous particle density.
521 * The results of this estimate have been checked againt simulations.
522 * In most cases the real drift differs by less than a factor 2.
525 /* Worst case assumption: HCP packing of particles gives largest distance */
526 particle_distance
= pow(boxvol
*sqrt(2)/mtop
->natoms
,1.0/3.0);
528 get_verlet_buffer_atomtypes(mtop
,&att
,&natt
,n_nonlin_vsite
);
529 assert(att
!= NULL
&& natt
>= 0);
533 fprintf(debug
,"particle distance assuming HCP packing: %f nm\n",
535 fprintf(debug
,"energy drift atom types: %d\n",natt
);
538 reppow
= mtop
->ffparams
.reppow
;
541 if (ir
->vdwtype
== evdwCUT
)
543 /* -dV/dr of -r^-6 and r^-repporw */
544 md_ljd
= -6*pow(ir
->rvdw
,-7.0);
545 md_ljr
= reppow
*pow(ir
->rvdw
,-(reppow
+1));
546 /* The contribution of the second derivative is negligible */
550 gmx_fatal(FARGS
,"Energy drift calculation is only implemented for plain cut-off Lennard-Jones interactions");
553 elfac
= ONE_4PI_EPS0
/ir
->epsilon_r
;
555 /* Determine md=-dV/dr and dd=d^2V/dr^2 */
558 if (ir
->coulombtype
== eelCUT
|| EEL_RF(ir
->coulombtype
))
562 if (ir
->coulombtype
== eelCUT
)
569 eps_rf
= ir
->epsilon_rf
/ir
->epsilon_r
;
572 k_rf
= pow(ir
->rcoulomb
,-3.0)*(eps_rf
- ir
->epsilon_r
)/(2*eps_rf
+ ir
->epsilon_r
);
576 /* epsilon_rf = infinity */
577 k_rf
= 0.5*pow(ir
->rcoulomb
,-3.0);
583 md_el
= elfac
*(pow(ir
->rcoulomb
,-2.0) - 2*k_rf
*ir
->rcoulomb
);
585 dd_el
= elfac
*(2*pow(ir
->rcoulomb
,-3.0) + 2*k_rf
);
587 else if (EEL_PME(ir
->coulombtype
) || ir
->coulombtype
== eelEWALD
)
591 b
= calc_ewaldcoeff(ir
->rcoulomb
,ir
->ewald_rtol
);
594 md_el
= elfac
*(2*b
*exp(-br
*br
)/(sqrt(M_PI
)*rc
) + gmx_erfc(br
)/(rc
*rc
));
595 dd_el
= elfac
/(rc
*rc
)*(4*b
*(1 + br
*br
)*exp(-br
*br
)/sqrt(M_PI
) + 2*gmx_erfc(br
)/rc
);
599 gmx_fatal(FARGS
,"Energy drift calculation is only implemented for Reaction-Field and Ewald electrostatics");
602 /* Determine the variance of the atomic displacement
603 * over nstlist-1 steps: kT_fac
604 * For inertial dynamics (not Brownian dynamics) the mass factor
605 * is not included in kT_fac, it is added later.
609 /* Get the displacement distribution from the random component only.
610 * With accurate integration the systematic (force) displacement
611 * should be negligible (unless nstlist is extremely large, which
612 * you wouldn't do anyhow).
614 kT_fac
= 2*BOLTZ
*ir
->opts
.ref_t
[0]*(ir
->nstlist
-1)*ir
->delta_t
;
617 /* This is directly sigma^2 of the displacement */
618 kT_fac
/= ir
->bd_fric
;
620 /* Set the masses to 1 as kT_fac is the full sigma^2,
621 * but we divide by m in ener_drift().
623 for(i
=0; i
<natt
; i
++)
632 /* Per group tau_t is not implemented yet, use the maximum */
633 tau_t
= ir
->opts
.tau_t
[0];
634 for(i
=1; i
<ir
->opts
.ngtc
; i
++)
636 tau_t
= max(tau_t
,ir
->opts
.tau_t
[i
]);
640 /* This kT_fac needs to be divided by the mass to get sigma^2 */
645 kT_fac
= BOLTZ
*ir
->opts
.ref_t
[0]*sqr((ir
->nstlist
-1)*ir
->delta_t
);
648 mass_min
= att
[0].mass
;
649 for(i
=1; i
<natt
; i
++)
651 mass_min
= min(mass_min
,att
[i
].mass
);
656 fprintf(debug
,"md_ljd %e md_ljr %e\n",md_ljd
,md_ljr
);
657 fprintf(debug
,"md_el %e dd_el %e\n",md_el
,dd_el
);
658 fprintf(debug
,"sqrt(kT_fac) %f\n",sqrt(kT_fac
));
659 fprintf(debug
,"mass_min %f\n",mass_min
);
662 /* Search using bisection */
664 /* The drift will be neglible at 5 times the max sigma */
665 ib1
= (int)(5*2*sqrt(kT_fac
/mass_min
)/resolution
) + 1;
666 while (ib1
- ib0
> 1)
670 rl
= max(ir
->rvdw
,ir
->rcoulomb
) + rb
;
672 /* Calculate the average energy drift at the last step
673 * of the nstlist steps at which the pair-list is used.
675 drift
= ener_drift(att
,natt
,&mtop
->ffparams
,
677 md_ljd
,md_ljr
,md_el
,dd_el
,rb
,
680 /* Correct for the fact that we are using a Ni x Nj particle pair list
681 * and not a 1 x 1 particle pair list. This reduces the drift.
683 /* We don't have a formula for 8 (yet), use 4 which is conservative */
684 nb_clust_frac_pairs_not_in_list_at_cutoff
=
685 surface_frac(min(list_setup
->cluster_size_i
,4),
686 particle_distance
,rl
)*
687 surface_frac(min(list_setup
->cluster_size_j
,4),
688 particle_distance
,rl
);
689 drift
*= nb_clust_frac_pairs_not_in_list_at_cutoff
;
691 /* Convert the drift to drift per unit time per atom */
692 drift
/= ir
->nstlist
*ir
->delta_t
*mtop
->natoms
;
696 fprintf(debug
,"ib %3d %3d %3d rb %.3f %dx%d fac %.3f drift %f\n",
698 list_setup
->cluster_size_i
,list_setup
->cluster_size_j
,
699 nb_clust_frac_pairs_not_in_list_at_cutoff
,
703 if (fabs(drift
) > drift_target
)
715 *rlist
= max(ir
->rvdw
,ir
->rcoulomb
) + ib1
*resolution
;