| File: | gromacs/gmxpreprocess/calc_verletbuf.c |
| Location: | line 998, column 16 |
| Description: | Dereference of null pointer |
| 1 | /* | |||
| 2 | * This file is part of the GROMACS molecular simulation package. | |||
| 3 | * | |||
| 4 | * Copyright (c) 2012,2013,2014, 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. | |||
| 8 | * | |||
| 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. | |||
| 13 | * | |||
| 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. | |||
| 18 | * | |||
| 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. | |||
| 23 | * | |||
| 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. | |||
| 31 | * | |||
| 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. | |||
| 34 | */ | |||
| 35 | #ifdef HAVE_CONFIG_H1 | |||
| 36 | #include <config.h> | |||
| 37 | #endif | |||
| 38 | ||||
| 39 | #include <assert.h> | |||
| 40 | #include <math.h> | |||
| 41 | #include <stdlib.h> | |||
| 42 | ||||
| 43 | #include <sys/types.h> | |||
| 44 | ||||
| 45 | #include "typedefs.h" | |||
| 46 | #include "physics.h" | |||
| 47 | #include "macros.h" | |||
| 48 | #include "gromacs/math/vec.h" | |||
| 49 | #include "coulomb.h" | |||
| 50 | #include "calc_verletbuf.h" | |||
| 51 | #include "../mdlib/nbnxn_consts.h" | |||
| 52 | ||||
| 53 | #include "gromacs/utility/fatalerror.h" | |||
| 54 | #include "gromacs/utility/smalloc.h" | |||
| 55 | ||||
| 56 | #ifdef GMX_NBNXN_SIMD | |||
| 57 | /* The include below sets the SIMD instruction type (precision+width) | |||
| 58 | * for all nbnxn SIMD search and non-bonded kernel code. | |||
| 59 | */ | |||
| 60 | #ifdef GMX_NBNXN_HALF_WIDTH_SIMD | |||
| 61 | #define GMX_USE_HALF_WIDTH_SIMD_HERE | |||
| 62 | #endif | |||
| 63 | #include "gromacs/simd/simd.h" | |||
| 64 | #endif | |||
| 65 | ||||
| 66 | ||||
| 67 | /* The code in this file estimates a pairlist buffer length | |||
| 68 | * given a target energy drift per atom per picosecond. | |||
| 69 | * This is done by estimating the drift given a buffer length. | |||
| 70 | * Ideally we would like to have a tight overestimate of the drift, | |||
| 71 | * but that can be difficult to achieve. | |||
| 72 | * | |||
| 73 | * Significant approximations used: | |||
| 74 | * | |||
| 75 | * Uniform particle density. UNDERESTIMATES the drift by rho_global/rho_local. | |||
| 76 | * | |||
| 77 | * Interactions don't affect particle motion. OVERESTIMATES the drift on longer | |||
| 78 | * time scales. This approximation probably introduces the largest errors. | |||
| 79 | * | |||
| 80 | * Only take one constraint per particle into account: OVERESTIMATES the drift. | |||
| 81 | * | |||
| 82 | * For rotating constraints assume the same functional shape for time scales | |||
| 83 | * where the constraints rotate significantly as the exact expression for | |||
| 84 | * short time scales. OVERESTIMATES the drift on long time scales. | |||
| 85 | * | |||
| 86 | * For non-linear virtual sites use the mass of the lightest constructing atom | |||
| 87 | * to determine the displacement. OVER/UNDERESTIMATES the drift, depending on | |||
| 88 | * the geometry and masses of constructing atoms. | |||
| 89 | * | |||
| 90 | * Note that the formulas for normal atoms and linear virtual sites are exact, | |||
| 91 | * apart from the first two approximations. | |||
| 92 | * | |||
| 93 | * Note that apart from the effect of the above approximations, the actual | |||
| 94 | * drift of the total energy of a system can be order of magnitude smaller | |||
| 95 | * due to cancellation of positive and negative drift for different pairs. | |||
| 96 | */ | |||
| 97 | ||||
| 98 | ||||
| 99 | /* Struct for unique atom type for calculating the energy drift. | |||
| 100 | * The atom displacement depends on mass and constraints. | |||
| 101 | * The energy jump for given distance depend on LJ type and q. | |||
| 102 | */ | |||
| 103 | typedef struct | |||
| 104 | { | |||
| 105 | real mass; /* mass */ | |||
| 106 | int type; /* type (used for LJ parameters) */ | |||
| 107 | real q; /* charge */ | |||
| 108 | gmx_bool bConstr; /* constrained, if TRUE, use #DOF=2 iso 3 */ | |||
| 109 | real con_mass; /* mass of heaviest atom connected by constraints */ | |||
| 110 | real con_len; /* constraint length to the heaviest atom */ | |||
| 111 | } atom_nonbonded_kinetic_prop_t; | |||
| 112 | ||||
| 113 | /* Struct for unique atom type for calculating the energy drift. | |||
| 114 | * The atom displacement depends on mass and constraints. | |||
| 115 | * The energy jump for given distance depend on LJ type and q. | |||
| 116 | */ | |||
| 117 | typedef struct | |||
| 118 | { | |||
| 119 | atom_nonbonded_kinetic_prop_t prop; /* non-bonded and kinetic atom prop. */ | |||
| 120 | int n; /* #atoms of this type in the system */ | |||
| 121 | } verletbuf_atomtype_t; | |||
| 122 | ||||
| 123 | void verletbuf_get_list_setup(gmx_bool bGPU, | |||
| 124 | verletbuf_list_setup_t *list_setup) | |||
| 125 | { | |||
| 126 | list_setup->cluster_size_i = NBNXN_CPU_CLUSTER_I_SIZE4; | |||
| 127 | ||||
| 128 | if (bGPU) | |||
| 129 | { | |||
| 130 | list_setup->cluster_size_j = NBNXN_GPU_CLUSTER_SIZE8; | |||
| 131 | } | |||
| 132 | else | |||
| 133 | { | |||
| 134 | #ifndef GMX_NBNXN_SIMD | |||
| 135 | list_setup->cluster_size_j = NBNXN_CPU_CLUSTER_I_SIZE4; | |||
| 136 | #else | |||
| 137 | list_setup->cluster_size_j = GMX_SIMD_REAL_WIDTH; | |||
| 138 | #ifdef GMX_NBNXN_SIMD_2XNN | |||
| 139 | /* We assume the smallest cluster size to be on the safe side */ | |||
| 140 | list_setup->cluster_size_j /= 2; | |||
| 141 | #endif | |||
| 142 | #endif | |||
| 143 | } | |||
| 144 | } | |||
| 145 | ||||
| 146 | static gmx_bool | |||
| 147 | atom_nonbonded_kinetic_prop_equal(const atom_nonbonded_kinetic_prop_t *prop1, | |||
| 148 | const atom_nonbonded_kinetic_prop_t *prop2) | |||
| 149 | { | |||
| 150 | return (prop1->mass == prop2->mass && | |||
| 151 | prop1->type == prop2->type && | |||
| 152 | prop1->q == prop2->q && | |||
| 153 | prop1->bConstr == prop2->bConstr && | |||
| 154 | prop1->con_mass == prop2->con_mass && | |||
| 155 | prop1->con_len == prop2->con_len); | |||
| 156 | } | |||
| 157 | ||||
| 158 | static void add_at(verletbuf_atomtype_t **att_p, int *natt_p, | |||
| 159 | const atom_nonbonded_kinetic_prop_t *prop, | |||
| 160 | int nmol) | |||
| 161 | { | |||
| 162 | verletbuf_atomtype_t *att; | |||
| 163 | int natt, i; | |||
| 164 | ||||
| 165 | if (prop->mass == 0) | |||
| 166 | { | |||
| 167 | /* Ignore massless particles */ | |||
| 168 | return; | |||
| 169 | } | |||
| 170 | ||||
| 171 | att = *att_p; | |||
| 172 | natt = *natt_p; | |||
| 173 | ||||
| 174 | i = 0; | |||
| 175 | while (i < natt && !atom_nonbonded_kinetic_prop_equal(prop, &att[i].prop)) | |||
| 176 | { | |||
| 177 | i++; | |||
| 178 | } | |||
| 179 | ||||
| 180 | if (i < natt) | |||
| 181 | { | |||
| 182 | att[i].n += nmol; | |||
| 183 | } | |||
| 184 | else | |||
| 185 | { | |||
| 186 | (*natt_p)++; | |||
| 187 | srenew(*att_p, *natt_p)(*att_p) = save_realloc("*att_p", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/calc_verletbuf.c" , 187, (*att_p), (*natt_p), sizeof(*(*att_p))); | |||
| 188 | (*att_p)[i].prop = *prop; | |||
| 189 | (*att_p)[i].n = nmol; | |||
| 190 | } | |||
| 191 | } | |||
| 192 | ||||
| 193 | static void get_vsite_masses(const gmx_moltype_t *moltype, | |||
| 194 | const gmx_ffparams_t *ffparams, | |||
| 195 | real *vsite_m, | |||
| 196 | int *n_nonlin_vsite) | |||
| 197 | { | |||
| 198 | int ft, i; | |||
| 199 | const t_ilist *il; | |||
| 200 | ||||
| 201 | *n_nonlin_vsite = 0; | |||
| 202 | ||||
| 203 | /* Check for virtual sites, determine mass from constructing atoms */ | |||
| 204 | for (ft = 0; ft < F_NRE; ft++) | |||
| 205 | { | |||
| 206 | if (IS_VSITE(ft)(interaction_function[(ft)].flags & 1<<1)) | |||
| 207 | { | |||
| 208 | il = &moltype->ilist[ft]; | |||
| 209 | ||||
| 210 | for (i = 0; i < il->nr; i += 1+NRAL(ft)(interaction_function[(ft)].nratoms)) | |||
| 211 | { | |||
| 212 | const t_iparams *ip; | |||
| 213 | real cam[5], inv_mass, m_aj; | |||
| 214 | int a1, j, aj, coeff; | |||
| 215 | ||||
| 216 | ip = &ffparams->iparams[il->iatoms[i]]; | |||
| 217 | ||||
| 218 | a1 = il->iatoms[i+1]; | |||
| 219 | ||||
| 220 | if (ft != F_VSITEN) | |||
| 221 | { | |||
| 222 | for (j = 1; j < NRAL(ft)(interaction_function[(ft)].nratoms); j++) | |||
| 223 | { | |||
| 224 | cam[j] = moltype->atoms.atom[il->iatoms[i+1+j]].m; | |||
| 225 | if (cam[j] == 0) | |||
| 226 | { | |||
| 227 | cam[j] = vsite_m[il->iatoms[i+1+j]]; | |||
| 228 | } | |||
| 229 | if (cam[j] == 0) | |||
| 230 | { | |||
| 231 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/calc_verletbuf.c" , 231, "In molecule type '%s' %s construction involves atom %d, which is a virtual site of equal or high complexity. This is not supported.", | |||
| 232 | *moltype->name, | |||
| 233 | interaction_function[ft].longname, | |||
| 234 | il->iatoms[i+1+j]+1); | |||
| 235 | } | |||
| 236 | } | |||
| 237 | } | |||
| 238 | ||||
| 239 | switch (ft) | |||
| 240 | { | |||
| 241 | case F_VSITE2: | |||
| 242 | /* Exact */ | |||
| 243 | vsite_m[a1] = (cam[1]*cam[2])/(cam[2]*sqr(1-ip->vsite.a) + cam[1]*sqr(ip->vsite.a)); | |||
| 244 | break; | |||
| 245 | case F_VSITE3: | |||
| 246 | /* Exact */ | |||
| 247 | vsite_m[a1] = (cam[1]*cam[2]*cam[3])/(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)); | |||
| 248 | break; | |||
| 249 | case F_VSITEN: | |||
| 250 | /* Exact */ | |||
| 251 | inv_mass = 0; | |||
| 252 | for (j = 0; j < 3*ip->vsiten.n; j += 3) | |||
| 253 | { | |||
| 254 | aj = il->iatoms[i+j+2]; | |||
| 255 | coeff = ip[il->iatoms[i+j]].vsiten.a; | |||
| 256 | if (moltype->atoms.atom[aj].ptype == eptVSite) | |||
| 257 | { | |||
| 258 | m_aj = vsite_m[aj]; | |||
| 259 | } | |||
| 260 | else | |||
| 261 | { | |||
| 262 | m_aj = moltype->atoms.atom[aj].m; | |||
| 263 | } | |||
| 264 | if (m_aj <= 0) | |||
| 265 | { | |||
| 266 | gmx_incons("The mass of a vsiten constructing atom is <= 0")_gmx_error("incons", "The mass of a vsiten constructing atom is <= 0" , "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/calc_verletbuf.c" , 266); | |||
| 267 | } | |||
| 268 | inv_mass += coeff*coeff/m_aj; | |||
| 269 | } | |||
| 270 | vsite_m[a1] = 1/inv_mass; | |||
| 271 | break; | |||
| 272 | default: | |||
| 273 | /* Use the mass of the lightest constructing atom. | |||
| 274 | * This is an approximation. | |||
| 275 | * If the distance of the virtual site to the | |||
| 276 | * constructing atom is less than all distances | |||
| 277 | * between constructing atoms, this is a safe | |||
| 278 | * over-estimate of the displacement of the vsite. | |||
| 279 | * This condition holds for all H mass replacement | |||
| 280 | * vsite constructions, except for SP2/3 groups. | |||
| 281 | * In SP3 groups one H will have a F_VSITE3 | |||
| 282 | * construction, so even there the total drift | |||
| 283 | * estimate shouldn't be far off. | |||
| 284 | */ | |||
| 285 | assert(j >= 1)((void) (0)); | |||
| 286 | vsite_m[a1] = cam[1]; | |||
| 287 | for (j = 2; j < NRAL(ft)(interaction_function[(ft)].nratoms); j++) | |||
| 288 | { | |||
| 289 | vsite_m[a1] = min(vsite_m[a1], cam[j])(((vsite_m[a1]) < (cam[j])) ? (vsite_m[a1]) : (cam[j]) ); | |||
| 290 | } | |||
| 291 | (*n_nonlin_vsite)++; | |||
| 292 | break; | |||
| 293 | } | |||
| 294 | if (gmx_debug_at) | |||
| 295 | { | |||
| 296 | fprintf(debug, "atom %4d %-20s mass %6.3f\n", | |||
| 297 | a1, interaction_function[ft].longname, vsite_m[a1]); | |||
| 298 | } | |||
| 299 | } | |||
| 300 | } | |||
| 301 | } | |||
| 302 | } | |||
| 303 | ||||
| 304 | static void get_verlet_buffer_atomtypes(const gmx_mtop_t *mtop, | |||
| 305 | verletbuf_atomtype_t **att_p, | |||
| 306 | int *natt_p, | |||
| 307 | int *n_nonlin_vsite) | |||
| 308 | { | |||
| 309 | verletbuf_atomtype_t *att; | |||
| 310 | int natt; | |||
| 311 | int mb, nmol, ft, i, a1, a2, a3, a; | |||
| 312 | const t_atoms *atoms; | |||
| 313 | const t_ilist *il; | |||
| 314 | const t_iparams *ip; | |||
| 315 | atom_nonbonded_kinetic_prop_t *prop; | |||
| 316 | real *vsite_m; | |||
| 317 | int n_nonlin_vsite_mol; | |||
| 318 | ||||
| 319 | att = NULL((void*)0); | |||
| 320 | natt = 0; | |||
| 321 | ||||
| 322 | if (n_nonlin_vsite != NULL((void*)0)) | |||
| 323 | { | |||
| 324 | *n_nonlin_vsite = 0; | |||
| 325 | } | |||
| 326 | ||||
| 327 | for (mb = 0; mb < mtop->nmolblock; mb++) | |||
| 328 | { | |||
| 329 | nmol = mtop->molblock[mb].nmol; | |||
| 330 | ||||
| 331 | atoms = &mtop->moltype[mtop->molblock[mb].type].atoms; | |||
| 332 | ||||
| 333 | /* Check for constraints, as they affect the kinetic energy. | |||
| 334 | * For virtual sites we need the masses and geometry of | |||
| 335 | * the constructing atoms to determine their velocity distribution. | |||
| 336 | */ | |||
| 337 | snew(prop, atoms->nr)(prop) = save_calloc("prop", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/calc_verletbuf.c" , 337, (atoms->nr), sizeof(*(prop))); | |||
| 338 | snew(vsite_m, atoms->nr)(vsite_m) = save_calloc("vsite_m", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/calc_verletbuf.c" , 338, (atoms->nr), sizeof(*(vsite_m))); | |||
| 339 | ||||
| 340 | for (ft = F_CONSTR; ft <= F_CONSTRNC; ft++) | |||
| 341 | { | |||
| 342 | il = &mtop->moltype[mtop->molblock[mb].type].ilist[ft]; | |||
| 343 | ||||
| 344 | for (i = 0; i < il->nr; i += 1+NRAL(ft)(interaction_function[(ft)].nratoms)) | |||
| 345 | { | |||
| 346 | ip = &mtop->ffparams.iparams[il->iatoms[i]]; | |||
| 347 | a1 = il->iatoms[i+1]; | |||
| 348 | a2 = il->iatoms[i+2]; | |||
| 349 | if (atoms->atom[a2].m > prop[a1].con_mass) | |||
| 350 | { | |||
| 351 | prop[a1].con_mass = atoms->atom[a2].m; | |||
| 352 | prop[a1].con_len = ip->constr.dA; | |||
| 353 | } | |||
| 354 | if (atoms->atom[a1].m > prop[a2].con_mass) | |||
| 355 | { | |||
| 356 | prop[a2].con_mass = atoms->atom[a1].m; | |||
| 357 | prop[a2].con_len = ip->constr.dA; | |||
| 358 | } | |||
| 359 | } | |||
| 360 | } | |||
| 361 | ||||
| 362 | il = &mtop->moltype[mtop->molblock[mb].type].ilist[F_SETTLE]; | |||
| 363 | ||||
| 364 | for (i = 0; i < il->nr; i += 1+NRAL(F_SETTLE)(interaction_function[(F_SETTLE)].nratoms)) | |||
| 365 | { | |||
| 366 | ip = &mtop->ffparams.iparams[il->iatoms[i]]; | |||
| 367 | a1 = il->iatoms[i+1]; | |||
| 368 | a2 = il->iatoms[i+2]; | |||
| 369 | a3 = il->iatoms[i+3]; | |||
| 370 | /* Usually the mass of a1 (usually oxygen) is larger than a2/a3. | |||
| 371 | * If this is not the case, we overestimate the displacement, | |||
| 372 | * which leads to a larger buffer (ok since this is an exotic case). | |||
| 373 | */ | |||
| 374 | prop[a1].con_mass = atoms->atom[a2].m; | |||
| 375 | prop[a1].con_len = ip->settle.doh; | |||
| 376 | ||||
| 377 | prop[a2].con_mass = atoms->atom[a1].m; | |||
| 378 | prop[a2].con_len = ip->settle.doh; | |||
| 379 | ||||
| 380 | prop[a3].con_mass = atoms->atom[a1].m; | |||
| 381 | prop[a3].con_len = ip->settle.doh; | |||
| 382 | } | |||
| 383 | ||||
| 384 | get_vsite_masses(&mtop->moltype[mtop->molblock[mb].type], | |||
| 385 | &mtop->ffparams, | |||
| 386 | vsite_m, | |||
| 387 | &n_nonlin_vsite_mol); | |||
| 388 | if (n_nonlin_vsite != NULL((void*)0)) | |||
| 389 | { | |||
| 390 | *n_nonlin_vsite += nmol*n_nonlin_vsite_mol; | |||
| 391 | } | |||
| 392 | ||||
| 393 | for (a = 0; a < atoms->nr; a++) | |||
| 394 | { | |||
| 395 | if (atoms->atom[a].ptype == eptVSite) | |||
| 396 | { | |||
| 397 | prop[a].mass = vsite_m[a]; | |||
| 398 | } | |||
| 399 | else | |||
| 400 | { | |||
| 401 | prop[a].mass = atoms->atom[a].m; | |||
| 402 | } | |||
| 403 | prop[a].type = atoms->atom[a].type; | |||
| 404 | prop[a].q = atoms->atom[a].q; | |||
| 405 | /* We consider an atom constrained, #DOF=2, when it is | |||
| 406 | * connected with constraints to (at least one) atom with | |||
| 407 | * a mass of more than 0.4x its own mass. This is not a critical | |||
| 408 | * parameter, since with roughly equal masses the unconstrained | |||
| 409 | * and constrained displacement will not differ much (and both | |||
| 410 | * overestimate the displacement). | |||
| 411 | */ | |||
| 412 | prop[a].bConstr = (prop[a].con_mass > 0.4*prop[a].mass); | |||
| 413 | ||||
| 414 | add_at(&att, &natt, &prop[a], nmol); | |||
| 415 | } | |||
| 416 | ||||
| 417 | sfree(vsite_m)save_free("vsite_m", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/calc_verletbuf.c" , 417, (vsite_m)); | |||
| 418 | sfree(prop)save_free("prop", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/calc_verletbuf.c" , 418, (prop)); | |||
| 419 | } | |||
| 420 | ||||
| 421 | if (gmx_debug_at) | |||
| 422 | { | |||
| 423 | for (a = 0; a < natt; a++) | |||
| 424 | { | |||
| 425 | 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", | |||
| 426 | a, att[a].prop.mass, att[a].prop.type, att[a].prop.q, | |||
| 427 | att[a].prop.bConstr, att[a].prop.con_mass, att[a].prop.con_len, | |||
| 428 | att[a].n); | |||
| 429 | } | |||
| 430 | } | |||
| 431 | ||||
| 432 | *att_p = att; | |||
| 433 | *natt_p = natt; | |||
| 434 | } | |||
| 435 | ||||
| 436 | /* This function computes two components of the estimate of the variance | |||
| 437 | * in the displacement of one atom in a system of two constrained atoms. | |||
| 438 | * Returns in sigma2_2d the variance due to rotation of the constrained | |||
| 439 | * atom around the atom to which it constrained. | |||
| 440 | * Returns in sigma2_3d the variance due to displacement of the COM | |||
| 441 | * of the whole system of the two constrained atoms. | |||
| 442 | * | |||
| 443 | * Note that we only take a single constraint (the one to the heaviest atom) | |||
| 444 | * into account. If an atom has multiple constraints, this will result in | |||
| 445 | * an overestimate of the displacement, which gives a larger drift and buffer. | |||
| 446 | */ | |||
| 447 | static void constrained_atom_sigma2(real kT_fac, | |||
| 448 | const atom_nonbonded_kinetic_prop_t *prop, | |||
| 449 | real *sigma2_2d, | |||
| 450 | real *sigma2_3d) | |||
| 451 | { | |||
| 452 | real sigma2_rot; | |||
| 453 | real com_dist; | |||
| 454 | real sigma2_rel; | |||
| 455 | real scale; | |||
| 456 | ||||
| 457 | /* Here we decompose the motion of a constrained atom into two | |||
| 458 | * components: rotation around the COM and translation of the COM. | |||
| 459 | */ | |||
| 460 | ||||
| 461 | /* Determine the variance for the displacement of the rotational mode */ | |||
| 462 | sigma2_rot = kT_fac/(prop->mass*(prop->mass + prop->con_mass)/prop->con_mass); | |||
| 463 | ||||
| 464 | /* The distance from the atom to the COM, i.e. the rotational arm */ | |||
| 465 | com_dist = prop->con_len*prop->con_mass/(prop->mass + prop->con_mass); | |||
| 466 | ||||
| 467 | /* The variance relative to the arm */ | |||
| 468 | sigma2_rel = sigma2_rot/(com_dist*com_dist); | |||
| 469 | /* At 6 the scaling formula has slope 0, | |||
| 470 | * so we keep sigma2_2d constant after that. | |||
| 471 | */ | |||
| 472 | if (sigma2_rel < 6) | |||
| 473 | { | |||
| 474 | /* A constrained atom rotates around the atom it is constrained to. | |||
| 475 | * This results in a smaller linear displacement than for a free atom. | |||
| 476 | * For a perfectly circular displacement, this lowers the displacement | |||
| 477 | * by: 1/arcsin(arc_length) | |||
| 478 | * and arcsin(x) = 1 + x^2/6 + ... | |||
| 479 | * For sigma2_rel<<1 the displacement distribution is erfc | |||
| 480 | * (exact formula is provided below). For larger sigma, it is clear | |||
| 481 | * that the displacement can't be larger than 2*com_dist. | |||
| 482 | * It turns out that the distribution becomes nearly uniform. | |||
| 483 | * For intermediate sigma2_rel, scaling down sigma with the third | |||
| 484 | * order expansion of arcsin with argument sigma_rel turns out | |||
| 485 | * to give a very good approximation of the distribution and variance. | |||
| 486 | * Even for larger values, the variance is only slightly overestimated. | |||
| 487 | * Note that the most relevant displacements are in the long tail. | |||
| 488 | * This rotation approximation always overestimates the tail (which | |||
| 489 | * runs to infinity, whereas it should be <= 2*com_dist). | |||
| 490 | * Thus we always overestimate the drift and the buffer size. | |||
| 491 | */ | |||
| 492 | scale = 1/(1 + sigma2_rel/6); | |||
| 493 | *sigma2_2d = sigma2_rot*scale*scale; | |||
| 494 | } | |||
| 495 | else | |||
| 496 | { | |||
| 497 | /* sigma_2d is set to the maximum given by the scaling above. | |||
| 498 | * For large sigma2 the real displacement distribution is close | |||
| 499 | * to uniform over -2*con_len to 2*com_dist. | |||
| 500 | * Our erfc with sigma_2d=sqrt(1.5)*com_dist (which means the sigma | |||
| 501 | * of the erfc output distribution is con_dist) overestimates | |||
| 502 | * the variance and additionally has a long tail. This means | |||
| 503 | * we have a (safe) overestimation of the drift. | |||
| 504 | */ | |||
| 505 | *sigma2_2d = 1.5*com_dist*com_dist; | |||
| 506 | } | |||
| 507 | ||||
| 508 | /* The constrained atom also moves (in 3D) with the COM of both atoms */ | |||
| 509 | *sigma2_3d = kT_fac/(prop->mass + prop->con_mass); | |||
| 510 | } | |||
| 511 | ||||
| 512 | static void get_atom_sigma2(real kT_fac, | |||
| 513 | const atom_nonbonded_kinetic_prop_t *prop, | |||
| 514 | real *sigma2_2d, | |||
| 515 | real *sigma2_3d) | |||
| 516 | { | |||
| 517 | if (prop->bConstr) | |||
| 518 | { | |||
| 519 | /* Complicated constraint calculation in a separate function */ | |||
| 520 | constrained_atom_sigma2(kT_fac, prop, sigma2_2d, sigma2_3d); | |||
| 521 | } | |||
| 522 | else | |||
| 523 | { | |||
| 524 | /* Unconstrained atom: trivial */ | |||
| 525 | *sigma2_2d = 0; | |||
| 526 | *sigma2_3d = kT_fac/prop->mass; | |||
| 527 | } | |||
| 528 | } | |||
| 529 | ||||
| 530 | static void approx_2dof(real s2, real x, real *shift, real *scale) | |||
| 531 | { | |||
| 532 | /* A particle with 1 DOF constrained has 2 DOFs instead of 3. | |||
| 533 | * This code is also used for particles with multiple constraints, | |||
| 534 | * in which case we overestimate the displacement. | |||
| 535 | * The 2DOF distribution is sqrt(pi/2)*erfc(r/(sqrt(2)*s))/(2*s). | |||
| 536 | * We approximate this with scale*Gaussian(s,r+shift), | |||
| 537 | * by matching the distribution value and derivative at x. | |||
| 538 | * This is a tight overestimate for all r>=0 at any s and x. | |||
| 539 | */ | |||
| 540 | real ex, er; | |||
| 541 | ||||
| 542 | ex = exp(-x*x/(2*s2)); | |||
| 543 | er = gmx_erfc(x/sqrt(2*s2))gmx_erfcf(x/sqrt(2*s2)); | |||
| 544 | ||||
| 545 | *shift = -x + sqrt(2*s2/M_PI3.14159265358979323846)*ex/er; | |||
| 546 | *scale = 0.5*M_PI3.14159265358979323846*exp(ex*ex/(M_PI3.14159265358979323846*er*er))*er; | |||
| 547 | } | |||
| 548 | ||||
| 549 | static real ener_drift(const verletbuf_atomtype_t *att, int natt, | |||
| 550 | const gmx_ffparams_t *ffp, | |||
| 551 | real kT_fac, | |||
| 552 | real md1_ljd, real d2_ljd, real md3_ljd, | |||
| 553 | real md1_ljr, real d2_ljr, real md3_ljr, | |||
| 554 | real md1_el, real d2_el, | |||
| 555 | real r_buffer, | |||
| 556 | real rlist, real boxvol) | |||
| 557 | { | |||
| 558 | double drift_tot, pot1, pot2, pot3, pot; | |||
| 559 | int i, j; | |||
| 560 | real s2i_2d, s2i_3d, s2j_2d, s2j_3d, s2, s; | |||
| 561 | int ti, tj; | |||
| 562 | real md1, d2, md3; | |||
| 563 | real sc_fac, rsh, rsh2; | |||
| 564 | double c_exp, c_erfc; | |||
| 565 | ||||
| 566 | drift_tot = 0; | |||
| 567 | ||||
| 568 | /* Loop over the different atom type pairs */ | |||
| 569 | for (i = 0; i < natt; i++) | |||
| 570 | { | |||
| 571 | get_atom_sigma2(kT_fac, &att[i].prop, &s2i_2d, &s2i_3d); | |||
| 572 | ti = att[i].prop.type; | |||
| 573 | ||||
| 574 | for (j = i; j < natt; j++) | |||
| 575 | { | |||
| 576 | get_atom_sigma2(kT_fac, &att[j].prop, &s2j_2d, &s2j_3d); | |||
| 577 | tj = att[j].prop.type; | |||
| 578 | ||||
| 579 | /* Add up the up to four independent variances */ | |||
| 580 | s2 = s2i_2d + s2i_3d + s2j_2d + s2j_3d; | |||
| 581 | ||||
| 582 | /* Note that attractive and repulsive potentials for individual | |||
| 583 | * pairs will partially cancel. | |||
| 584 | */ | |||
| 585 | /* -dV/dr at the cut-off for LJ + Coulomb */ | |||
| 586 | md1 = | |||
| 587 | md1_ljd*ffp->iparams[ti*ffp->atnr+tj].lj.c6 + | |||
| 588 | md1_ljr*ffp->iparams[ti*ffp->atnr+tj].lj.c12 + | |||
| 589 | md1_el*att[i].prop.q*att[j].prop.q; | |||
| 590 | ||||
| 591 | /* d2V/dr2 at the cut-off for LJ + Coulomb */ | |||
| 592 | d2 = | |||
| 593 | d2_ljd*ffp->iparams[ti*ffp->atnr+tj].lj.c6 + | |||
| 594 | d2_ljr*ffp->iparams[ti*ffp->atnr+tj].lj.c12 + | |||
| 595 | d2_el*att[i].prop.q*att[j].prop.q; | |||
| 596 | ||||
| 597 | /* -d3V/dr3 at the cut-off for LJ, we neglect Coulomb */ | |||
| 598 | md3 = | |||
| 599 | md3_ljd*ffp->iparams[ti*ffp->atnr+tj].lj.c6 + | |||
| 600 | md3_ljr*ffp->iparams[ti*ffp->atnr+tj].lj.c12; | |||
| 601 | ||||
| 602 | rsh = r_buffer; | |||
| 603 | sc_fac = 1.0; | |||
| 604 | /* For constraints: adapt r and scaling for the Gaussian */ | |||
| 605 | if (att[i].prop.bConstr) | |||
| 606 | { | |||
| 607 | real sh, sc; | |||
| 608 | ||||
| 609 | approx_2dof(s2i_2d, r_buffer*s2i_2d/s2, &sh, &sc); | |||
| 610 | rsh += sh; | |||
| 611 | sc_fac *= sc; | |||
| 612 | } | |||
| 613 | if (att[j].prop.bConstr) | |||
| 614 | { | |||
| 615 | real sh, sc; | |||
| 616 | ||||
| 617 | approx_2dof(s2j_2d, r_buffer*s2j_2d/s2, &sh, &sc); | |||
| 618 | rsh += sh; | |||
| 619 | sc_fac *= sc; | |||
| 620 | } | |||
| 621 | ||||
| 622 | /* Exact contribution of an atom pair with Gaussian displacement | |||
| 623 | * with sigma s to the energy drift for a potential with | |||
| 624 | * derivative -md and second derivative dd at the cut-off. | |||
| 625 | * The only catch is that for potentials that change sign | |||
| 626 | * near the cut-off there could be an unlucky compensation | |||
| 627 | * of positive and negative energy drift. | |||
| 628 | * Such potentials are extremely rare though. | |||
| 629 | * | |||
| 630 | * Note that pot has unit energy*length, as the linear | |||
| 631 | * atom density still needs to be put in. | |||
| 632 | */ | |||
| 633 | c_exp = exp(-rsh*rsh/(2*s2))/sqrt(2*M_PI3.14159265358979323846); | |||
| 634 | c_erfc = 0.5*gmx_erfc(rsh/(sqrt(2*s2)))gmx_erfcf(rsh/(sqrt(2*s2))); | |||
| 635 | s = sqrt(s2); | |||
| 636 | rsh2 = rsh*rsh; | |||
| 637 | ||||
| 638 | pot1 = sc_fac* | |||
| 639 | md1/2*((rsh2 + s2)*c_erfc - rsh*s*c_exp); | |||
| 640 | pot2 = sc_fac* | |||
| 641 | d2/6*(s*(rsh2 + 2*s2)*c_exp - rsh*(rsh2 + 3*s2)*c_erfc); | |||
| 642 | pot3 = | |||
| 643 | md3/24*((rsh2*rsh2 + 6*rsh2*s2 + 3*s2*s2)*c_erfc - rsh*s*(rsh2 + 5*s2)*c_exp); | |||
| 644 | pot = pot1 + pot2 + pot3; | |||
| 645 | ||||
| 646 | if (gmx_debug_at) | |||
| 647 | { | |||
| 648 | fprintf(debug, "n %d %d d s %.3f %.3f %.3f %.3f con %d -d1 %8.1e d2 %8.1e -d3 %8.1e pot1 %8.1e pot2 %8.1e pot3 %8.1e pot %8.1e\n", | |||
| 649 | att[i].n, att[j].n, | |||
| 650 | sqrt(s2i_2d), sqrt(s2i_3d), | |||
| 651 | sqrt(s2j_2d), sqrt(s2j_3d), | |||
| 652 | att[i].prop.bConstr+att[j].prop.bConstr, | |||
| 653 | md1, d2, md3, | |||
| 654 | pot1, pot2, pot3, pot); | |||
| 655 | } | |||
| 656 | ||||
| 657 | /* Multiply by the number of atom pairs */ | |||
| 658 | if (j == i) | |||
| 659 | { | |||
| 660 | pot *= (double)att[i].n*(att[i].n - 1)/2; | |||
| 661 | } | |||
| 662 | else | |||
| 663 | { | |||
| 664 | pot *= (double)att[i].n*att[j].n; | |||
| 665 | } | |||
| 666 | /* We need the line density to get the energy drift of the system. | |||
| 667 | * The effective average r^2 is close to (rlist+sigma)^2. | |||
| 668 | */ | |||
| 669 | pot *= 4*M_PI3.14159265358979323846*sqr(rlist + s)/boxvol; | |||
| 670 | ||||
| 671 | /* Add the unsigned drift to avoid cancellation of errors */ | |||
| 672 | drift_tot += fabs(pot); | |||
| 673 | } | |||
| 674 | } | |||
| 675 | ||||
| 676 | return drift_tot; | |||
| 677 | } | |||
| 678 | ||||
| 679 | static real surface_frac(int cluster_size, real particle_distance, real rlist) | |||
| 680 | { | |||
| 681 | real d, area_rel; | |||
| 682 | ||||
| 683 | if (rlist < 0.5*particle_distance) | |||
| 684 | { | |||
| 685 | /* We have non overlapping spheres */ | |||
| 686 | return 1.0; | |||
| 687 | } | |||
| 688 | ||||
| 689 | /* Half the inter-particle distance relative to rlist */ | |||
| 690 | d = 0.5*particle_distance/rlist; | |||
| 691 | ||||
| 692 | /* Determine the area of the surface at distance rlist to the closest | |||
| 693 | * particle, relative to surface of a sphere of radius rlist. | |||
| 694 | * The formulas below assume close to cubic cells for the pair search grid, | |||
| 695 | * which the pair search code tries to achieve. | |||
| 696 | * Note that in practice particle distances will not be delta distributed, | |||
| 697 | * but have some spread, often involving shorter distances, | |||
| 698 | * as e.g. O-H bonds in a water molecule. Thus the estimates below will | |||
| 699 | * usually be slightly too high and thus conservative. | |||
| 700 | */ | |||
| 701 | switch (cluster_size) | |||
| 702 | { | |||
| 703 | case 1: | |||
| 704 | /* One particle: trivial */ | |||
| 705 | area_rel = 1.0; | |||
| 706 | break; | |||
| 707 | case 2: | |||
| 708 | /* Two particles: two spheres at fractional distance 2*a */ | |||
| 709 | area_rel = 1.0 + d; | |||
| 710 | break; | |||
| 711 | case 4: | |||
| 712 | /* We assume a perfect, symmetric tetrahedron geometry. | |||
| 713 | * The surface around a tetrahedron is too complex for a full | |||
| 714 | * analytical solution, so we use a Taylor expansion. | |||
| 715 | */ | |||
| 716 | area_rel = (1.0 + 1/M_PI3.14159265358979323846*(6*acos(1/sqrt(3))*d + | |||
| 717 | sqrt(3)*d*d*(1.0 + | |||
| 718 | 5.0/18.0*d*d + | |||
| 719 | 7.0/45.0*d*d*d*d + | |||
| 720 | 83.0/756.0*d*d*d*d*d*d))); | |||
| 721 | break; | |||
| 722 | default: | |||
| 723 | gmx_incons("surface_frac called with unsupported cluster_size")_gmx_error("incons", "surface_frac called with unsupported cluster_size" , "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/calc_verletbuf.c" , 723); | |||
| 724 | area_rel = 1.0; | |||
| 725 | } | |||
| 726 | ||||
| 727 | return area_rel/cluster_size; | |||
| 728 | } | |||
| 729 | ||||
| 730 | /* Returns the negative of the third derivative of a potential r^-p | |||
| 731 | * with a force-switch function, evaluated at the cut-off rc. | |||
| 732 | */ | |||
| 733 | static real md3_force_switch(real p, real rswitch, real rc) | |||
| 734 | { | |||
| 735 | /* The switched force function is: | |||
| 736 | * p*r^-(p+1) + a*(r - rswitch)^2 + b*(r - rswitch)^3 | |||
| 737 | */ | |||
| 738 | real a, b; | |||
| 739 | real md3_pot, md3_sw; | |||
| 740 | ||||
| 741 | a = -((p + 4)*rc - (p + 1)*rswitch)/(pow(rc, p+2)*pow(rc-rswitch, 2)); | |||
| 742 | b = ((p + 3)*rc - (p + 1)*rswitch)/(pow(rc, p+2)*pow(rc-rswitch, 3)); | |||
| 743 | ||||
| 744 | md3_pot = (p + 2)*(p + 1)*p*pow(rc, p+3); | |||
| 745 | md3_sw = 2*a + 6*b*(rc - rswitch); | |||
| 746 | ||||
| 747 | return md3_pot + md3_sw; | |||
| 748 | } | |||
| 749 | ||||
| 750 | void calc_verlet_buffer_size(const gmx_mtop_t *mtop, real boxvol, | |||
| 751 | const t_inputrec *ir, | |||
| 752 | real reference_temperature, | |||
| 753 | const verletbuf_list_setup_t *list_setup, | |||
| 754 | int *n_nonlin_vsite, | |||
| 755 | real *rlist) | |||
| 756 | { | |||
| 757 | double resolution; | |||
| 758 | char *env; | |||
| 759 | ||||
| 760 | real particle_distance; | |||
| 761 | real nb_clust_frac_pairs_not_in_list_at_cutoff; | |||
| 762 | ||||
| 763 | verletbuf_atomtype_t *att = NULL((void*)0); | |||
| 764 | int natt = -1, i; | |||
| 765 | double reppow; | |||
| 766 | real md1_ljd, d2_ljd, md3_ljd; | |||
| 767 | real md1_ljr, d2_ljr, md3_ljr; | |||
| 768 | real md1_el, d2_el; | |||
| 769 | real elfac; | |||
| 770 | real kT_fac, mass_min; | |||
| 771 | int ib0, ib1, ib; | |||
| 772 | real rb, rl; | |||
| 773 | real drift; | |||
| 774 | ||||
| 775 | if (reference_temperature < 0) | |||
| ||||
| 776 | { | |||
| 777 | if (EI_MD(ir->eI)((ir->eI) == eiMD || ((ir->eI) == eiVV || (ir->eI) == eiVVAK)) && ir->etc == etcNO) | |||
| 778 | { | |||
| 779 | /* This case should be handled outside calc_verlet_buffer_size */ | |||
| 780 | gmx_incons("calc_verlet_buffer_size called with an NVE ensemble and reference_temperature < 0")_gmx_error("incons", "calc_verlet_buffer_size called with an NVE ensemble and reference_temperature < 0" , "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/calc_verletbuf.c" , 780); | |||
| 781 | } | |||
| 782 | ||||
| 783 | /* We use the maximum temperature with multiple T-coupl groups. | |||
| 784 | * We could use a per particle temperature, but since particles | |||
| 785 | * interact, this might underestimate the buffer size. | |||
| 786 | */ | |||
| 787 | reference_temperature = 0; | |||
| 788 | for (i = 0; i < ir->opts.ngtc; i++) | |||
| 789 | { | |||
| 790 | if (ir->opts.tau_t[i] >= 0) | |||
| 791 | { | |||
| 792 | reference_temperature = max(reference_temperature,(((reference_temperature) > (ir->opts.ref_t[i])) ? (reference_temperature ) : (ir->opts.ref_t[i]) ) | |||
| 793 | ir->opts.ref_t[i])(((reference_temperature) > (ir->opts.ref_t[i])) ? (reference_temperature ) : (ir->opts.ref_t[i]) ); | |||
| 794 | } | |||
| 795 | } | |||
| 796 | } | |||
| 797 | ||||
| 798 | /* Resolution of the buffer size */ | |||
| 799 | resolution = 0.001; | |||
| 800 | ||||
| 801 | env = getenv("GMX_VERLET_BUFFER_RES"); | |||
| 802 | if (env != NULL((void*)0)) | |||
| 803 | { | |||
| 804 | sscanf(env, "%lf", &resolution); | |||
| 805 | } | |||
| 806 | ||||
| 807 | /* In an atom wise pair-list there would be no pairs in the list | |||
| 808 | * beyond the pair-list cut-off. | |||
| 809 | * However, we use a pair-list of groups vs groups of atoms. | |||
| 810 | * For groups of 4 atoms, the parallelism of SSE instructions, only | |||
| 811 | * 10% of the atoms pairs are not in the list just beyond the cut-off. | |||
| 812 | * As this percentage increases slowly compared to the decrease of the | |||
| 813 | * Gaussian displacement distribution over this range, we can simply | |||
| 814 | * reduce the drift by this fraction. | |||
| 815 | * For larger groups, e.g. of 8 atoms, this fraction will be lower, | |||
| 816 | * so then buffer size will be on the conservative (large) side. | |||
| 817 | * | |||
| 818 | * Note that the formulas used here do not take into account | |||
| 819 | * cancellation of errors which could occur by missing both | |||
| 820 | * attractive and repulsive interactions. | |||
| 821 | * | |||
| 822 | * The only major assumption is homogeneous particle distribution. | |||
| 823 | * For an inhomogeneous system, such as a liquid-vapor system, | |||
| 824 | * the buffer will be underestimated. The actual energy drift | |||
| 825 | * will be higher by the factor: local/homogeneous particle density. | |||
| 826 | * | |||
| 827 | * The results of this estimate have been checked againt simulations. | |||
| 828 | * In most cases the real drift differs by less than a factor 2. | |||
| 829 | */ | |||
| 830 | ||||
| 831 | /* Worst case assumption: HCP packing of particles gives largest distance */ | |||
| 832 | particle_distance = pow(boxvol*sqrt(2)/mtop->natoms, 1.0/3.0); | |||
| 833 | ||||
| 834 | get_verlet_buffer_atomtypes(mtop, &att, &natt, n_nonlin_vsite); | |||
| 835 | assert(att != NULL && natt >= 0)((void) (0)); | |||
| 836 | ||||
| 837 | if (debug) | |||
| 838 | { | |||
| 839 | fprintf(debug, "particle distance assuming HCP packing: %f nm\n", | |||
| 840 | particle_distance); | |||
| 841 | fprintf(debug, "energy drift atom types: %d\n", natt); | |||
| 842 | } | |||
| 843 | ||||
| 844 | reppow = mtop->ffparams.reppow; | |||
| 845 | md1_ljd = 0; | |||
| 846 | d2_ljd = 0; | |||
| 847 | md3_ljd = 0; | |||
| 848 | md1_ljr = 0; | |||
| 849 | d2_ljr = 0; | |||
| 850 | md3_ljr = 0; | |||
| 851 | if (ir->vdwtype == evdwCUT) | |||
| 852 | { | |||
| 853 | real sw_range, md3_pswf; | |||
| 854 | ||||
| 855 | switch (ir->vdw_modifier) | |||
| 856 | { | |||
| 857 | case eintmodNONE: | |||
| 858 | case eintmodPOTSHIFT: | |||
| 859 | /* -dV/dr of -r^-6 and r^-reppow */ | |||
| 860 | md1_ljd = -6*pow(ir->rvdw, -7.0); | |||
| 861 | md1_ljr = reppow*pow(ir->rvdw, -(reppow+1)); | |||
| 862 | /* The contribution of the higher derivatives is negligible */ | |||
| 863 | break; | |||
| 864 | case eintmodFORCESWITCH: | |||
| 865 | /* At the cut-off: V=V'=V''=0, so we use only V''' */ | |||
| 866 | md3_ljd = -md3_force_switch(6.0, ir->rvdw_switch, ir->rvdw); | |||
| 867 | md3_ljr = md3_force_switch(reppow, ir->rvdw_switch, ir->rvdw); | |||
| 868 | break; | |||
| 869 | case eintmodPOTSWITCH: | |||
| 870 | /* At the cut-off: V=V'=V''=0. | |||
| 871 | * V''' is given by the original potential times | |||
| 872 | * the third derivative of the switch function. | |||
| 873 | */ | |||
| 874 | sw_range = ir->rvdw - ir->rvdw_switch; | |||
| 875 | md3_pswf = 60.0*pow(sw_range, -3.0); | |||
| 876 | ||||
| 877 | md3_ljd = -pow(ir->rvdw, -6.0 )*md3_pswf; | |||
| 878 | md3_ljr = pow(ir->rvdw, -reppow)*md3_pswf; | |||
| 879 | break; | |||
| 880 | default: | |||
| 881 | gmx_incons("Unimplemented VdW modifier")_gmx_error("incons", "Unimplemented VdW modifier", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/calc_verletbuf.c" , 881); | |||
| 882 | } | |||
| 883 | } | |||
| 884 | else if (EVDW_PME(ir->vdwtype)((ir->vdwtype) == evdwPME)) | |||
| 885 | { | |||
| 886 | real b, r, br, br2, br4, br6; | |||
| 887 | b = calc_ewaldcoeff_lj(ir->rvdw, ir->ewald_rtol_lj); | |||
| 888 | r = ir->rvdw; | |||
| 889 | br = b*r; | |||
| 890 | br2 = br*br; | |||
| 891 | br4 = br2*br2; | |||
| 892 | br6 = br4*br2; | |||
| 893 | /* -dV/dr of g(br)*r^-6 [where g(x) = exp(-x^2)(1+x^2+x^4/2), see LJ-PME equations in manual] and r^-reppow */ | |||
| 894 | md1_ljd = -exp(-br2)*(br6 + 3.0*br4 + 6.0*br2 + 6.0)*pow(r, -7.0); | |||
| 895 | md1_ljr = reppow*pow(r, -(reppow+1)); | |||
| 896 | /* The contribution of the higher derivatives is negligible */ | |||
| 897 | } | |||
| 898 | else | |||
| 899 | { | |||
| 900 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/calc_verletbuf.c" , 900, "Energy drift calculation is only implemented for plain cut-off Lennard-Jones interactions"); | |||
| 901 | } | |||
| 902 | ||||
| 903 | elfac = ONE_4PI_EPS0((332.0636930*(4.184))*0.1)/ir->epsilon_r; | |||
| 904 | ||||
| 905 | /* Determine md=-dV/dr and dd=d^2V/dr^2 */ | |||
| 906 | md1_el = 0; | |||
| 907 | d2_el = 0; | |||
| 908 | if (ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype)((ir->coulombtype) == eelRF || (ir->coulombtype) == eelGRF || (ir->coulombtype) == eelRF_NEC || (ir->coulombtype) == eelRF_ZERO )) | |||
| 909 | { | |||
| 910 | real eps_rf, k_rf; | |||
| 911 | ||||
| 912 | if (ir->coulombtype == eelCUT) | |||
| 913 | { | |||
| 914 | eps_rf = 1; | |||
| 915 | k_rf = 0; | |||
| 916 | } | |||
| 917 | else | |||
| 918 | { | |||
| 919 | eps_rf = ir->epsilon_rf/ir->epsilon_r; | |||
| 920 | if (eps_rf != 0) | |||
| 921 | { | |||
| 922 | k_rf = pow(ir->rcoulomb, -3.0)*(eps_rf - ir->epsilon_r)/(2*eps_rf + ir->epsilon_r); | |||
| 923 | } | |||
| 924 | else | |||
| 925 | { | |||
| 926 | /* epsilon_rf = infinity */ | |||
| 927 | k_rf = 0.5*pow(ir->rcoulomb, -3.0); | |||
| 928 | } | |||
| 929 | } | |||
| 930 | ||||
| 931 | if (eps_rf > 0) | |||
| 932 | { | |||
| 933 | md1_el = elfac*(pow(ir->rcoulomb, -2.0) - 2*k_rf*ir->rcoulomb); | |||
| 934 | } | |||
| 935 | d2_el = elfac*(2*pow(ir->rcoulomb, -3.0) + 2*k_rf); | |||
| 936 | } | |||
| 937 | else if (EEL_PME(ir->coulombtype)((ir->coulombtype) == eelPME || (ir->coulombtype) == eelPMESWITCH || (ir->coulombtype) == eelPMEUSER || (ir->coulombtype ) == eelPMEUSERSWITCH || (ir->coulombtype) == eelP3M_AD) || ir->coulombtype == eelEWALD) | |||
| 938 | { | |||
| 939 | real b, rc, br; | |||
| 940 | ||||
| 941 | b = calc_ewaldcoeff_q(ir->rcoulomb, ir->ewald_rtol); | |||
| 942 | rc = ir->rcoulomb; | |||
| 943 | br = b*rc; | |||
| 944 | md1_el = elfac*(b*exp(-br*br)*M_2_SQRTPI1.12837916709551257390/rc + gmx_erfc(br)gmx_erfcf(br)/(rc*rc)); | |||
| 945 | d2_el = elfac/(rc*rc)*(2*b*(1 + br*br)*exp(-br*br)*M_2_SQRTPI1.12837916709551257390 + 2*gmx_erfc(br)gmx_erfcf(br)/rc); | |||
| 946 | } | |||
| 947 | else | |||
| 948 | { | |||
| 949 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/calc_verletbuf.c" , 949, "Energy drift calculation is only implemented for Reaction-Field and Ewald electrostatics"); | |||
| 950 | } | |||
| 951 | ||||
| 952 | /* Determine the variance of the atomic displacement | |||
| 953 | * over nstlist-1 steps: kT_fac | |||
| 954 | * For inertial dynamics (not Brownian dynamics) the mass factor | |||
| 955 | * is not included in kT_fac, it is added later. | |||
| 956 | */ | |||
| 957 | if (ir->eI == eiBD) | |||
| 958 | { | |||
| 959 | /* Get the displacement distribution from the random component only. | |||
| 960 | * With accurate integration the systematic (force) displacement | |||
| 961 | * should be negligible (unless nstlist is extremely large, which | |||
| 962 | * you wouldn't do anyhow). | |||
| 963 | */ | |||
| 964 | kT_fac = 2*BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))*reference_temperature*(ir->nstlist-1)*ir->delta_t; | |||
| 965 | if (ir->bd_fric > 0) | |||
| 966 | { | |||
| 967 | /* This is directly sigma^2 of the displacement */ | |||
| 968 | kT_fac /= ir->bd_fric; | |||
| 969 | ||||
| 970 | /* Set the masses to 1 as kT_fac is the full sigma^2, | |||
| 971 | * but we divide by m in ener_drift(). | |||
| 972 | */ | |||
| 973 | for (i = 0; i < natt; i++) | |||
| 974 | { | |||
| 975 | att[i].prop.mass = 1; | |||
| 976 | } | |||
| 977 | } | |||
| 978 | else | |||
| 979 | { | |||
| 980 | real tau_t; | |||
| 981 | ||||
| 982 | /* Per group tau_t is not implemented yet, use the maximum */ | |||
| 983 | tau_t = ir->opts.tau_t[0]; | |||
| 984 | for (i = 1; i < ir->opts.ngtc; i++) | |||
| 985 | { | |||
| 986 | tau_t = max(tau_t, ir->opts.tau_t[i])(((tau_t) > (ir->opts.tau_t[i])) ? (tau_t) : (ir->opts .tau_t[i]) ); | |||
| 987 | } | |||
| 988 | ||||
| 989 | kT_fac *= tau_t; | |||
| 990 | /* This kT_fac needs to be divided by the mass to get sigma^2 */ | |||
| 991 | } | |||
| 992 | } | |||
| 993 | else | |||
| 994 | { | |||
| 995 | kT_fac = BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))*reference_temperature*sqr((ir->nstlist-1)*ir->delta_t); | |||
| 996 | } | |||
| 997 | ||||
| 998 | mass_min = att[0].prop.mass; | |||
| ||||
| 999 | for (i = 1; i < natt; i++) | |||
| 1000 | { | |||
| 1001 | mass_min = min(mass_min, att[i].prop.mass)(((mass_min) < (att[i].prop.mass)) ? (mass_min) : (att[i]. prop.mass) ); | |||
| 1002 | } | |||
| 1003 | ||||
| 1004 | if (debug) | |||
| 1005 | { | |||
| 1006 | fprintf(debug, "md1_ljd %9.2e d2_ljd %9.2e md3_ljd %9.2e\n", md1_ljd, d2_ljd, md3_ljd); | |||
| 1007 | fprintf(debug, "md1_ljr %9.2e d2_ljr %9.2e md3_ljr %9.2e\n", md1_ljr, d2_ljr, md3_ljr); | |||
| 1008 | fprintf(debug, "md1_el %9.2e d2_el %9.2e\n", md1_el, d2_el); | |||
| 1009 | fprintf(debug, "sqrt(kT_fac) %f\n", sqrt(kT_fac)); | |||
| 1010 | fprintf(debug, "mass_min %f\n", mass_min); | |||
| 1011 | } | |||
| 1012 | ||||
| 1013 | /* Search using bisection */ | |||
| 1014 | ib0 = -1; | |||
| 1015 | /* The drift will be neglible at 5 times the max sigma */ | |||
| 1016 | ib1 = (int)(5*2*sqrt(kT_fac/mass_min)/resolution) + 1; | |||
| 1017 | while (ib1 - ib0 > 1) | |||
| 1018 | { | |||
| 1019 | ib = (ib0 + ib1)/2; | |||
| 1020 | rb = ib*resolution; | |||
| 1021 | rl = max(ir->rvdw, ir->rcoulomb)(((ir->rvdw) > (ir->rcoulomb)) ? (ir->rvdw) : (ir ->rcoulomb) ) + rb; | |||
| 1022 | ||||
| 1023 | /* Calculate the average energy drift at the last step | |||
| 1024 | * of the nstlist steps at which the pair-list is used. | |||
| 1025 | */ | |||
| 1026 | drift = ener_drift(att, natt, &mtop->ffparams, | |||
| 1027 | kT_fac, | |||
| 1028 | md1_ljd, d2_ljd, md3_ljd, | |||
| 1029 | md1_ljr, d2_ljr, md3_ljr, | |||
| 1030 | md1_el, d2_el, | |||
| 1031 | rb, | |||
| 1032 | rl, boxvol); | |||
| 1033 | ||||
| 1034 | /* Correct for the fact that we are using a Ni x Nj particle pair list | |||
| 1035 | * and not a 1 x 1 particle pair list. This reduces the drift. | |||
| 1036 | */ | |||
| 1037 | /* We don't have a formula for 8 (yet), use 4 which is conservative */ | |||
| 1038 | nb_clust_frac_pairs_not_in_list_at_cutoff = | |||
| 1039 | surface_frac(min(list_setup->cluster_size_i, 4)(((list_setup->cluster_size_i) < (4)) ? (list_setup-> cluster_size_i) : (4) ), | |||
| 1040 | particle_distance, rl)* | |||
| 1041 | surface_frac(min(list_setup->cluster_size_j, 4)(((list_setup->cluster_size_j) < (4)) ? (list_setup-> cluster_size_j) : (4) ), | |||
| 1042 | particle_distance, rl); | |||
| 1043 | drift *= nb_clust_frac_pairs_not_in_list_at_cutoff; | |||
| 1044 | ||||
| 1045 | /* Convert the drift to drift per unit time per atom */ | |||
| 1046 | drift /= ir->nstlist*ir->delta_t*mtop->natoms; | |||
| 1047 | ||||
| 1048 | if (debug) | |||
| 1049 | { | |||
| 1050 | fprintf(debug, "ib %3d %3d %3d rb %.3f %dx%d fac %.3f drift %f\n", | |||
| 1051 | ib0, ib, ib1, rb, | |||
| 1052 | list_setup->cluster_size_i, list_setup->cluster_size_j, | |||
| 1053 | nb_clust_frac_pairs_not_in_list_at_cutoff, | |||
| 1054 | drift); | |||
| 1055 | } | |||
| 1056 | ||||
| 1057 | if (fabs(drift) > ir->verletbuf_tol) | |||
| 1058 | { | |||
| 1059 | ib0 = ib; | |||
| 1060 | } | |||
| 1061 | else | |||
| 1062 | { | |||
| 1063 | ib1 = ib; | |||
| 1064 | } | |||
| 1065 | } | |||
| 1066 | ||||
| 1067 | sfree(att)save_free("att", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/calc_verletbuf.c" , 1067, (att)); | |||
| 1068 | ||||
| 1069 | *rlist = max(ir->rvdw, ir->rcoulomb)(((ir->rvdw) > (ir->rcoulomb)) ? (ir->rvdw) : (ir ->rcoulomb) ) + ib1*resolution; | |||
| 1070 | } |