| File: | gromacs/mdlib/shellfc.c |
| Location: | line 188, column 26 |
| Description: | Dereference of undefined pointer value |
| 1 | /* | |||
| 2 | * This file is part of the GROMACS molecular simulation package. | |||
| 3 | * | |||
| 4 | * Copyright (c) 1991-2000, University of Groningen, The Netherlands. | |||
| 5 | * Copyright (c) 2001-2008, The GROMACS development team. | |||
| 6 | * Copyright (c) 2013,2014, by the GROMACS development team, led by | |||
| 7 | * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, | |||
| 8 | * and including many others, as listed in the AUTHORS file in the | |||
| 9 | * top-level source directory and at http://www.gromacs.org. | |||
| 10 | * | |||
| 11 | * GROMACS is free software; you can redistribute it and/or | |||
| 12 | * modify it under the terms of the GNU Lesser General Public License | |||
| 13 | * as published by the Free Software Foundation; either version 2.1 | |||
| 14 | * of the License, or (at your option) any later version. | |||
| 15 | * | |||
| 16 | * GROMACS is distributed in the hope that it will be useful, | |||
| 17 | * but WITHOUT ANY WARRANTY; without even the implied warranty of | |||
| 18 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |||
| 19 | * Lesser General Public License for more details. | |||
| 20 | * | |||
| 21 | * You should have received a copy of the GNU Lesser General Public | |||
| 22 | * License along with GROMACS; if not, see | |||
| 23 | * http://www.gnu.org/licenses, or write to the Free Software Foundation, | |||
| 24 | * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. | |||
| 25 | * | |||
| 26 | * If you want to redistribute modifications to GROMACS, please | |||
| 27 | * consider that scientific software is very special. Version | |||
| 28 | * control is crucial - bugs must be traceable. We will be happy to | |||
| 29 | * consider code for inclusion in the official distribution, but | |||
| 30 | * derived work must not be called official GROMACS. Details are found | |||
| 31 | * in the README & COPYING files - if they are missing, get the | |||
| 32 | * official version at http://www.gromacs.org. | |||
| 33 | * | |||
| 34 | * To help us fund GROMACS development, we humbly ask that you cite | |||
| 35 | * the research papers on the package. Check out http://www.gromacs.org. | |||
| 36 | */ | |||
| 37 | #ifdef HAVE_CONFIG_H1 | |||
| 38 | #include <config.h> | |||
| 39 | #endif | |||
| 40 | ||||
| 41 | #include <stdlib.h> | |||
| 42 | #include <string.h> | |||
| 43 | ||||
| 44 | #include "typedefs.h" | |||
| 45 | #include "types/commrec.h" | |||
| 46 | #include "gromacs/utility/smalloc.h" | |||
| 47 | #include "gromacs/utility/fatalerror.h" | |||
| 48 | #include "gromacs/math/vec.h" | |||
| 49 | #include "txtdump.h" | |||
| 50 | #include "force.h" | |||
| 51 | #include "mdrun.h" | |||
| 52 | #include "mdatoms.h" | |||
| 53 | #include "vsite.h" | |||
| 54 | #include "network.h" | |||
| 55 | #include "names.h" | |||
| 56 | #include "constr.h" | |||
| 57 | #include "domdec.h" | |||
| 58 | #include "physics.h" | |||
| 59 | #include "shellfc.h" | |||
| 60 | #include "mtop_util.h" | |||
| 61 | #include "chargegroup.h" | |||
| 62 | #include "macros.h" | |||
| 63 | ||||
| 64 | ||||
| 65 | typedef struct { | |||
| 66 | int nnucl; | |||
| 67 | atom_id shell; /* The shell id */ | |||
| 68 | atom_id nucl1, nucl2, nucl3; /* The nuclei connected to the shell */ | |||
| 69 | /* gmx_bool bInterCG; */ /* Coupled to nuclei outside cg? */ | |||
| 70 | real k; /* force constant */ | |||
| 71 | real k_1; /* 1 over force constant */ | |||
| 72 | rvec xold; | |||
| 73 | rvec fold; | |||
| 74 | rvec step; | |||
| 75 | } t_shell; | |||
| 76 | ||||
| 77 | typedef struct gmx_shellfc { | |||
| 78 | int nshell_gl; /* The number of shells in the system */ | |||
| 79 | t_shell *shell_gl; /* All the shells (for DD only) */ | |||
| 80 | int *shell_index_gl; /* Global shell index (for DD only) */ | |||
| 81 | gmx_bool bInterCG; /* Are there inter charge-group shells? */ | |||
| 82 | int nshell; /* The number of local shells */ | |||
| 83 | t_shell *shell; /* The local shells */ | |||
| 84 | int shell_nalloc; /* The allocation size of shell */ | |||
| 85 | gmx_bool bPredict; /* Predict shell positions */ | |||
| 86 | gmx_bool bRequireInit; /* Require initialization of shell positions */ | |||
| 87 | int nflexcon; /* The number of flexible constraints */ | |||
| 88 | rvec *x[2]; /* Array for iterative minimization */ | |||
| 89 | rvec *f[2]; /* Array for iterative minimization */ | |||
| 90 | int x_nalloc; /* The allocation size of x and f */ | |||
| 91 | rvec *acc_dir; /* Acceleration direction for flexcon */ | |||
| 92 | rvec *x_old; /* Old coordinates for flexcon */ | |||
| 93 | int flex_nalloc; /* The allocation size of acc_dir and x_old */ | |||
| 94 | rvec *adir_xnold; /* Work space for init_adir */ | |||
| 95 | rvec *adir_xnew; /* Work space for init_adir */ | |||
| 96 | int adir_nalloc; /* Work space for init_adir */ | |||
| 97 | } t_gmx_shellfc; | |||
| 98 | ||||
| 99 | ||||
| 100 | static void pr_shell(FILE *fplog, int ns, t_shell s[]) | |||
| 101 | { | |||
| 102 | int i; | |||
| 103 | ||||
| 104 | fprintf(fplog, "SHELL DATA\n"); | |||
| 105 | fprintf(fplog, "%5s %8s %5s %5s %5s\n", | |||
| 106 | "Shell", "Force k", "Nucl1", "Nucl2", "Nucl3"); | |||
| 107 | for (i = 0; (i < ns); i++) | |||
| 108 | { | |||
| 109 | fprintf(fplog, "%5d %8.3f %5d", s[i].shell, 1.0/s[i].k_1, s[i].nucl1); | |||
| 110 | if (s[i].nnucl == 2) | |||
| 111 | { | |||
| 112 | fprintf(fplog, " %5d\n", s[i].nucl2); | |||
| 113 | } | |||
| 114 | else if (s[i].nnucl == 3) | |||
| 115 | { | |||
| 116 | fprintf(fplog, " %5d %5d\n", s[i].nucl2, s[i].nucl3); | |||
| 117 | } | |||
| 118 | else | |||
| 119 | { | |||
| 120 | fprintf(fplog, "\n"); | |||
| 121 | } | |||
| 122 | } | |||
| 123 | } | |||
| 124 | ||||
| 125 | static void predict_shells(FILE *fplog, rvec x[], rvec v[], real dt, | |||
| 126 | int ns, t_shell s[], | |||
| 127 | real mass[], gmx_mtop_t *mtop, gmx_bool bInit) | |||
| 128 | { | |||
| 129 | int i, m, s1, n1, n2, n3; | |||
| 130 | real dt_1, dt_2, dt_3, fudge, tm, m1, m2, m3; | |||
| 131 | rvec *ptr; | |||
| 132 | gmx_mtop_atomlookup_t alook = NULL((void*)0); | |||
| 133 | t_atom *atom; | |||
| ||||
| 134 | ||||
| 135 | if (mass == NULL((void*)0)) | |||
| 136 | { | |||
| 137 | alook = gmx_mtop_atomlookup_init(mtop); | |||
| 138 | } | |||
| 139 | ||||
| 140 | /* We introduce a fudge factor for performance reasons: with this choice | |||
| 141 | * the initial force on the shells is about a factor of two lower than | |||
| 142 | * without | |||
| 143 | */ | |||
| 144 | fudge = 1.0; | |||
| 145 | ||||
| 146 | if (bInit) | |||
| 147 | { | |||
| 148 | if (fplog) | |||
| 149 | { | |||
| 150 | fprintf(fplog, "RELAX: Using prediction for initial shell placement\n"); | |||
| 151 | } | |||
| 152 | ptr = x; | |||
| 153 | dt_1 = 1; | |||
| 154 | } | |||
| 155 | else | |||
| 156 | { | |||
| 157 | ptr = v; | |||
| 158 | dt_1 = fudge*dt; | |||
| 159 | } | |||
| 160 | ||||
| 161 | for (i = 0; (i < ns); i++) | |||
| 162 | { | |||
| 163 | s1 = s[i].shell; | |||
| 164 | if (bInit) | |||
| 165 | { | |||
| 166 | clear_rvec(x[s1]); | |||
| 167 | } | |||
| 168 | switch (s[i].nnucl) | |||
| 169 | { | |||
| 170 | case 1: | |||
| 171 | n1 = s[i].nucl1; | |||
| 172 | for (m = 0; (m < DIM3); m++) | |||
| 173 | { | |||
| 174 | x[s1][m] += ptr[n1][m]*dt_1; | |||
| 175 | } | |||
| 176 | break; | |||
| 177 | case 2: | |||
| 178 | n1 = s[i].nucl1; | |||
| 179 | n2 = s[i].nucl2; | |||
| 180 | if (mass) | |||
| 181 | { | |||
| 182 | m1 = mass[n1]; | |||
| 183 | m2 = mass[n2]; | |||
| 184 | } | |||
| 185 | else | |||
| 186 | { | |||
| 187 | /* Not the correct masses with FE, but it is just a prediction... */ | |||
| 188 | m1 = atom[n1].m; | |||
| ||||
| 189 | m2 = atom[n2].m; | |||
| 190 | } | |||
| 191 | tm = dt_1/(m1+m2); | |||
| 192 | for (m = 0; (m < DIM3); m++) | |||
| 193 | { | |||
| 194 | x[s1][m] += (m1*ptr[n1][m]+m2*ptr[n2][m])*tm; | |||
| 195 | } | |||
| 196 | break; | |||
| 197 | case 3: | |||
| 198 | n1 = s[i].nucl1; | |||
| 199 | n2 = s[i].nucl2; | |||
| 200 | n3 = s[i].nucl3; | |||
| 201 | if (mass) | |||
| 202 | { | |||
| 203 | m1 = mass[n1]; | |||
| 204 | m2 = mass[n2]; | |||
| 205 | m3 = mass[n3]; | |||
| 206 | } | |||
| 207 | else | |||
| 208 | { | |||
| 209 | /* Not the correct masses with FE, but it is just a prediction... */ | |||
| 210 | gmx_mtop_atomnr_to_atom(alook, n1, &atom); | |||
| 211 | m1 = atom->m; | |||
| 212 | gmx_mtop_atomnr_to_atom(alook, n2, &atom); | |||
| 213 | m2 = atom->m; | |||
| 214 | gmx_mtop_atomnr_to_atom(alook, n3, &atom); | |||
| 215 | m3 = atom->m; | |||
| 216 | } | |||
| 217 | tm = dt_1/(m1+m2+m3); | |||
| 218 | for (m = 0; (m < DIM3); m++) | |||
| 219 | { | |||
| 220 | x[s1][m] += (m1*ptr[n1][m]+m2*ptr[n2][m]+m3*ptr[n3][m])*tm; | |||
| 221 | } | |||
| 222 | break; | |||
| 223 | default: | |||
| 224 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/shellfc.c" , 224, "Shell %d has %d nuclei!", i, s[i].nnucl); | |||
| 225 | } | |||
| 226 | } | |||
| 227 | ||||
| 228 | if (mass == NULL((void*)0)) | |||
| 229 | { | |||
| 230 | gmx_mtop_atomlookup_destroy(alook); | |||
| 231 | } | |||
| 232 | } | |||
| 233 | ||||
| 234 | gmx_shellfc_t init_shell_flexcon(FILE *fplog, | |||
| 235 | gmx_bool bCutoffSchemeIsVerlet, | |||
| 236 | gmx_mtop_t *mtop, int nflexcon, | |||
| 237 | rvec *x) | |||
| 238 | { | |||
| 239 | struct gmx_shellfc *shfc; | |||
| 240 | t_shell *shell; | |||
| 241 | int *shell_index = NULL((void*)0), *at2cg; | |||
| 242 | t_atom *atom; | |||
| 243 | int n[eptNR], ns, nshell, nsi; | |||
| 244 | int i, j, nmol, type, mb, mt, a_offset, cg, mol, ftype, nra; | |||
| 245 | real qS, alpha; | |||
| 246 | int aS, aN = 0; /* Shell and nucleus */ | |||
| 247 | int bondtypes[] = { F_BONDS, F_HARMONIC, F_CUBICBONDS, F_POLARIZATION, F_ANHARM_POL, F_WATER_POL }; | |||
| 248 | #define NBT((int)(sizeof(bondtypes)/sizeof((bondtypes)[0]))) asize(bondtypes)((int)(sizeof(bondtypes)/sizeof((bondtypes)[0]))) | |||
| 249 | t_iatom *ia; | |||
| 250 | gmx_mtop_atomloop_block_t aloopb; | |||
| 251 | gmx_mtop_atomloop_all_t aloop; | |||
| 252 | gmx_ffparams_t *ffparams; | |||
| 253 | gmx_molblock_t *molb; | |||
| 254 | gmx_moltype_t *molt; | |||
| 255 | t_block *cgs; | |||
| 256 | ||||
| 257 | /* Count number of shells, and find their indices */ | |||
| 258 | for (i = 0; (i < eptNR); i++) | |||
| 259 | { | |||
| 260 | n[i] = 0; | |||
| 261 | } | |||
| 262 | ||||
| 263 | aloopb = gmx_mtop_atomloop_block_init(mtop); | |||
| 264 | while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol)) | |||
| 265 | { | |||
| 266 | n[atom->ptype] += nmol; | |||
| 267 | } | |||
| 268 | ||||
| 269 | if (fplog) | |||
| 270 | { | |||
| 271 | /* Print the number of each particle type */ | |||
| 272 | for (i = 0; (i < eptNR); i++) | |||
| 273 | { | |||
| 274 | if (n[i] != 0) | |||
| 275 | { | |||
| 276 | fprintf(fplog, "There are: %d %ss\n", n[i], ptype_str[i]); | |||
| 277 | } | |||
| 278 | } | |||
| 279 | } | |||
| 280 | ||||
| 281 | nshell = n[eptShell]; | |||
| 282 | ||||
| 283 | if (nshell == 0 && nflexcon == 0) | |||
| 284 | { | |||
| 285 | /* We're not doing shells or flexible constraints */ | |||
| 286 | return NULL((void*)0); | |||
| 287 | } | |||
| 288 | ||||
| 289 | if (bCutoffSchemeIsVerlet) | |||
| 290 | { | |||
| 291 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/shellfc.c" , 291, "The shell code does not work with the Verlet cut-off scheme.\n"); | |||
| 292 | } | |||
| 293 | ||||
| 294 | snew(shfc, 1)(shfc) = save_calloc("shfc", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/shellfc.c" , 294, (1), sizeof(*(shfc))); | |||
| 295 | shfc->nflexcon = nflexcon; | |||
| 296 | ||||
| 297 | if (nshell == 0) | |||
| 298 | { | |||
| 299 | return shfc; | |||
| 300 | } | |||
| 301 | ||||
| 302 | /* We have shells: fill the shell data structure */ | |||
| 303 | ||||
| 304 | /* Global system sized array, this should be avoided */ | |||
| 305 | snew(shell_index, mtop->natoms)(shell_index) = save_calloc("shell_index", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/shellfc.c" , 305, (mtop->natoms), sizeof(*(shell_index))); | |||
| 306 | ||||
| 307 | aloop = gmx_mtop_atomloop_all_init(mtop); | |||
| 308 | nshell = 0; | |||
| 309 | while (gmx_mtop_atomloop_all_next(aloop, &i, &atom)) | |||
| 310 | { | |||
| 311 | if (atom->ptype == eptShell) | |||
| 312 | { | |||
| 313 | shell_index[i] = nshell++; | |||
| 314 | } | |||
| 315 | } | |||
| 316 | ||||
| 317 | snew(shell, nshell)(shell) = save_calloc("shell", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/shellfc.c" , 317, (nshell), sizeof(*(shell))); | |||
| 318 | ||||
| 319 | /* Initiate the shell structures */ | |||
| 320 | for (i = 0; (i < nshell); i++) | |||
| 321 | { | |||
| 322 | shell[i].shell = NO_ATID(atom_id)(~0); | |||
| 323 | shell[i].nnucl = 0; | |||
| 324 | shell[i].nucl1 = NO_ATID(atom_id)(~0); | |||
| 325 | shell[i].nucl2 = NO_ATID(atom_id)(~0); | |||
| 326 | shell[i].nucl3 = NO_ATID(atom_id)(~0); | |||
| 327 | /* shell[i].bInterCG=FALSE; */ | |||
| 328 | shell[i].k_1 = 0; | |||
| 329 | shell[i].k = 0; | |||
| 330 | } | |||
| 331 | ||||
| 332 | ffparams = &mtop->ffparams; | |||
| 333 | ||||
| 334 | /* Now fill the structures */ | |||
| 335 | shfc->bInterCG = FALSE0; | |||
| 336 | ns = 0; | |||
| 337 | a_offset = 0; | |||
| 338 | for (mb = 0; mb < mtop->nmolblock; mb++) | |||
| 339 | { | |||
| 340 | molb = &mtop->molblock[mb]; | |||
| 341 | molt = &mtop->moltype[molb->type]; | |||
| 342 | ||||
| 343 | cgs = &molt->cgs; | |||
| 344 | snew(at2cg, molt->atoms.nr)(at2cg) = save_calloc("at2cg", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/shellfc.c" , 344, (molt->atoms.nr), sizeof(*(at2cg))); | |||
| 345 | for (cg = 0; cg < cgs->nr; cg++) | |||
| 346 | { | |||
| 347 | for (i = cgs->index[cg]; i < cgs->index[cg+1]; i++) | |||
| 348 | { | |||
| 349 | at2cg[i] = cg; | |||
| 350 | } | |||
| 351 | } | |||
| 352 | ||||
| 353 | atom = molt->atoms.atom; | |||
| 354 | for (mol = 0; mol < molb->nmol; mol++) | |||
| 355 | { | |||
| 356 | for (j = 0; (j < NBT((int)(sizeof(bondtypes)/sizeof((bondtypes)[0])))); j++) | |||
| 357 | { | |||
| 358 | ia = molt->ilist[bondtypes[j]].iatoms; | |||
| 359 | for (i = 0; (i < molt->ilist[bondtypes[j]].nr); ) | |||
| 360 | { | |||
| 361 | type = ia[0]; | |||
| 362 | ftype = ffparams->functype[type]; | |||
| 363 | nra = interaction_function[ftype].nratoms; | |||
| 364 | ||||
| 365 | /* Check whether we have a bond with a shell */ | |||
| 366 | aS = NO_ATID(atom_id)(~0); | |||
| 367 | ||||
| 368 | switch (bondtypes[j]) | |||
| 369 | { | |||
| 370 | case F_BONDS: | |||
| 371 | case F_HARMONIC: | |||
| 372 | case F_CUBICBONDS: | |||
| 373 | case F_POLARIZATION: | |||
| 374 | case F_ANHARM_POL: | |||
| 375 | if (atom[ia[1]].ptype == eptShell) | |||
| 376 | { | |||
| 377 | aS = ia[1]; | |||
| 378 | aN = ia[2]; | |||
| 379 | } | |||
| 380 | else if (atom[ia[2]].ptype == eptShell) | |||
| 381 | { | |||
| 382 | aS = ia[2]; | |||
| 383 | aN = ia[1]; | |||
| 384 | } | |||
| 385 | break; | |||
| 386 | case F_WATER_POL: | |||
| 387 | aN = ia[4]; /* Dummy */ | |||
| 388 | aS = ia[5]; /* Shell */ | |||
| 389 | break; | |||
| 390 | default: | |||
| 391 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/shellfc.c" , 391, "Death Horror: %s, %d", __FILE__"/home/alexxy/Develop/gromacs/src/gromacs/mdlib/shellfc.c", __LINE__391); | |||
| 392 | } | |||
| 393 | ||||
| 394 | if (aS != NO_ATID(atom_id)(~0)) | |||
| 395 | { | |||
| 396 | qS = atom[aS].q; | |||
| 397 | ||||
| 398 | /* Check whether one of the particles is a shell... */ | |||
| 399 | nsi = shell_index[a_offset+aS]; | |||
| 400 | if ((nsi < 0) || (nsi >= nshell)) | |||
| 401 | { | |||
| 402 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/shellfc.c" , 402, "nsi is %d should be within 0 - %d. aS = %d", | |||
| 403 | nsi, nshell, aS); | |||
| 404 | } | |||
| 405 | if (shell[nsi].shell == NO_ATID(atom_id)(~0)) | |||
| 406 | { | |||
| 407 | shell[nsi].shell = a_offset + aS; | |||
| 408 | ns++; | |||
| 409 | } | |||
| 410 | else if (shell[nsi].shell != a_offset+aS) | |||
| 411 | { | |||
| 412 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/shellfc.c" , 412, "Weird stuff in %s, %d", __FILE__"/home/alexxy/Develop/gromacs/src/gromacs/mdlib/shellfc.c", __LINE__412); | |||
| 413 | } | |||
| 414 | ||||
| 415 | if (shell[nsi].nucl1 == NO_ATID(atom_id)(~0)) | |||
| 416 | { | |||
| 417 | shell[nsi].nucl1 = a_offset + aN; | |||
| 418 | } | |||
| 419 | else if (shell[nsi].nucl2 == NO_ATID(atom_id)(~0)) | |||
| 420 | { | |||
| 421 | shell[nsi].nucl2 = a_offset + aN; | |||
| 422 | } | |||
| 423 | else if (shell[nsi].nucl3 == NO_ATID(atom_id)(~0)) | |||
| 424 | { | |||
| 425 | shell[nsi].nucl3 = a_offset + aN; | |||
| 426 | } | |||
| 427 | else | |||
| 428 | { | |||
| 429 | if (fplog) | |||
| 430 | { | |||
| 431 | pr_shell(fplog, ns, shell); | |||
| 432 | } | |||
| 433 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/shellfc.c" , 433, "Can not handle more than three bonds per shell\n"); | |||
| 434 | } | |||
| 435 | if (at2cg[aS] != at2cg[aN]) | |||
| 436 | { | |||
| 437 | /* shell[nsi].bInterCG = TRUE; */ | |||
| 438 | shfc->bInterCG = TRUE1; | |||
| 439 | } | |||
| 440 | ||||
| 441 | switch (bondtypes[j]) | |||
| 442 | { | |||
| 443 | case F_BONDS: | |||
| 444 | case F_HARMONIC: | |||
| 445 | shell[nsi].k += ffparams->iparams[type].harmonic.krA; | |||
| 446 | break; | |||
| 447 | case F_CUBICBONDS: | |||
| 448 | shell[nsi].k += ffparams->iparams[type].cubic.kb; | |||
| 449 | break; | |||
| 450 | case F_POLARIZATION: | |||
| 451 | case F_ANHARM_POL: | |||
| 452 | if (!gmx_within_tol(qS, atom[aS].qB, GMX_REAL_EPS5.96046448E-08*10)) | |||
| 453 | { | |||
| 454 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/shellfc.c" , 454, "polarize can not be used with qA(%e) != qB(%e) for atom %d of molecule block %d", qS, atom[aS].qB, aS+1, mb+1); | |||
| 455 | } | |||
| 456 | shell[nsi].k += sqr(qS)*ONE_4PI_EPS0((332.0636930*(4.184))*0.1)/ | |||
| 457 | ffparams->iparams[type].polarize.alpha; | |||
| 458 | break; | |||
| 459 | case F_WATER_POL: | |||
| 460 | if (!gmx_within_tol(qS, atom[aS].qB, GMX_REAL_EPS5.96046448E-08*10)) | |||
| 461 | { | |||
| 462 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/shellfc.c" , 462, "water_pol can not be used with qA(%e) != qB(%e) for atom %d of molecule block %d", qS, atom[aS].qB, aS+1, mb+1); | |||
| 463 | } | |||
| 464 | alpha = (ffparams->iparams[type].wpol.al_x+ | |||
| 465 | ffparams->iparams[type].wpol.al_y+ | |||
| 466 | ffparams->iparams[type].wpol.al_z)/3.0; | |||
| 467 | shell[nsi].k += sqr(qS)*ONE_4PI_EPS0((332.0636930*(4.184))*0.1)/alpha; | |||
| 468 | break; | |||
| 469 | default: | |||
| 470 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/shellfc.c" , 470, "Death Horror: %s, %d", __FILE__"/home/alexxy/Develop/gromacs/src/gromacs/mdlib/shellfc.c", __LINE__470); | |||
| 471 | } | |||
| 472 | shell[nsi].nnucl++; | |||
| 473 | } | |||
| 474 | ia += nra+1; | |||
| 475 | i += nra+1; | |||
| 476 | } | |||
| 477 | } | |||
| 478 | a_offset += molt->atoms.nr; | |||
| 479 | } | |||
| 480 | /* Done with this molecule type */ | |||
| 481 | sfree(at2cg)save_free("at2cg", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/shellfc.c" , 481, (at2cg)); | |||
| 482 | } | |||
| 483 | ||||
| 484 | /* Verify whether it's all correct */ | |||
| 485 | if (ns != nshell) | |||
| 486 | { | |||
| 487 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/shellfc.c" , 487, "Something weird with shells. They may not be bonded to something"); | |||
| 488 | } | |||
| 489 | ||||
| 490 | for (i = 0; (i < ns); i++) | |||
| 491 | { | |||
| 492 | shell[i].k_1 = 1.0/shell[i].k; | |||
| 493 | } | |||
| 494 | ||||
| 495 | if (debug) | |||
| 496 | { | |||
| 497 | pr_shell(debug, ns, shell); | |||
| 498 | } | |||
| 499 | ||||
| 500 | ||||
| 501 | shfc->nshell_gl = ns; | |||
| 502 | shfc->shell_gl = shell; | |||
| 503 | shfc->shell_index_gl = shell_index; | |||
| 504 | ||||
| 505 | shfc->bPredict = (getenv("GMX_NOPREDICT") == NULL((void*)0)); | |||
| 506 | shfc->bRequireInit = FALSE0; | |||
| 507 | if (!shfc->bPredict) | |||
| 508 | { | |||
| 509 | if (fplog) | |||
| 510 | { | |||
| 511 | fprintf(fplog, "\nWill never predict shell positions\n"); | |||
| 512 | } | |||
| 513 | } | |||
| 514 | else | |||
| 515 | { | |||
| 516 | shfc->bRequireInit = (getenv("GMX_REQUIRE_SHELL_INIT") != NULL((void*)0)); | |||
| 517 | if (shfc->bRequireInit && fplog) | |||
| 518 | { | |||
| 519 | fprintf(fplog, "\nWill always initiate shell positions\n"); | |||
| 520 | } | |||
| 521 | } | |||
| 522 | ||||
| 523 | if (shfc->bPredict) | |||
| 524 | { | |||
| 525 | if (x) | |||
| 526 | { | |||
| 527 | predict_shells(fplog, x, NULL((void*)0), 0, shfc->nshell_gl, shfc->shell_gl, | |||
| 528 | NULL((void*)0), mtop, TRUE1); | |||
| 529 | } | |||
| 530 | ||||
| 531 | if (shfc->bInterCG) | |||
| 532 | { | |||
| 533 | if (fplog) | |||
| 534 | { | |||
| 535 | fprintf(fplog, "\nNOTE: there all shells that are connected to particles outside thier own charge group, will not predict shells positions during the run\n\n"); | |||
| 536 | } | |||
| 537 | shfc->bPredict = FALSE0; | |||
| 538 | } | |||
| 539 | } | |||
| 540 | ||||
| 541 | return shfc; | |||
| 542 | } | |||
| 543 | ||||
| 544 | void make_local_shells(t_commrec *cr, t_mdatoms *md, | |||
| 545 | struct gmx_shellfc *shfc) | |||
| 546 | { | |||
| 547 | t_shell *shell; | |||
| 548 | int a0, a1, *ind, nshell, i; | |||
| 549 | gmx_domdec_t *dd = NULL((void*)0); | |||
| 550 | ||||
| 551 | if (DOMAINDECOMP(cr)(((cr)->dd != ((void*)0)) && ((cr)->nnodes > 1))) | |||
| 552 | { | |||
| 553 | dd = cr->dd; | |||
| 554 | a0 = 0; | |||
| 555 | a1 = dd->nat_home; | |||
| 556 | } | |||
| 557 | else | |||
| 558 | { | |||
| 559 | /* Single node: we need all shells, just copy the pointer */ | |||
| 560 | shfc->nshell = shfc->nshell_gl; | |||
| 561 | shfc->shell = shfc->shell_gl; | |||
| 562 | ||||
| 563 | return; | |||
| 564 | } | |||
| 565 | ||||
| 566 | ind = shfc->shell_index_gl; | |||
| 567 | ||||
| 568 | nshell = 0; | |||
| 569 | shell = shfc->shell; | |||
| 570 | for (i = a0; i < a1; i++) | |||
| 571 | { | |||
| 572 | if (md->ptype[i] == eptShell) | |||
| 573 | { | |||
| 574 | if (nshell+1 > shfc->shell_nalloc) | |||
| 575 | { | |||
| 576 | shfc->shell_nalloc = over_alloc_dd(nshell+1); | |||
| 577 | srenew(shell, shfc->shell_nalloc)(shell) = save_realloc("shell", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/shellfc.c" , 577, (shell), (shfc->shell_nalloc), sizeof(*(shell))); | |||
| 578 | } | |||
| 579 | if (dd) | |||
| 580 | { | |||
| 581 | shell[nshell] = shfc->shell_gl[ind[dd->gatindex[i]]]; | |||
| 582 | } | |||
| 583 | else | |||
| 584 | { | |||
| 585 | shell[nshell] = shfc->shell_gl[ind[i]]; | |||
| 586 | } | |||
| 587 | ||||
| 588 | /* With inter-cg shells we can no do shell prediction, | |||
| 589 | * so we do not need the nuclei numbers. | |||
| 590 | */ | |||
| 591 | if (!shfc->bInterCG) | |||
| 592 | { | |||
| 593 | shell[nshell].nucl1 = i + shell[nshell].nucl1 - shell[nshell].shell; | |||
| 594 | if (shell[nshell].nnucl > 1) | |||
| 595 | { | |||
| 596 | shell[nshell].nucl2 = i + shell[nshell].nucl2 - shell[nshell].shell; | |||
| 597 | } | |||
| 598 | if (shell[nshell].nnucl > 2) | |||
| 599 | { | |||
| 600 | shell[nshell].nucl3 = i + shell[nshell].nucl3 - shell[nshell].shell; | |||
| 601 | } | |||
| 602 | } | |||
| 603 | shell[nshell].shell = i; | |||
| 604 | nshell++; | |||
| 605 | } | |||
| 606 | } | |||
| 607 | ||||
| 608 | shfc->nshell = nshell; | |||
| 609 | shfc->shell = shell; | |||
| 610 | } | |||
| 611 | ||||
| 612 | static void do_1pos(rvec xnew, rvec xold, rvec f, real step) | |||
| 613 | { | |||
| 614 | real xo, yo, zo; | |||
| 615 | real dx, dy, dz; | |||
| 616 | ||||
| 617 | xo = xold[XX0]; | |||
| 618 | yo = xold[YY1]; | |||
| 619 | zo = xold[ZZ2]; | |||
| 620 | ||||
| 621 | dx = f[XX0]*step; | |||
| 622 | dy = f[YY1]*step; | |||
| 623 | dz = f[ZZ2]*step; | |||
| 624 | ||||
| 625 | xnew[XX0] = xo+dx; | |||
| 626 | xnew[YY1] = yo+dy; | |||
| 627 | xnew[ZZ2] = zo+dz; | |||
| 628 | } | |||
| 629 | ||||
| 630 | static void do_1pos3(rvec xnew, rvec xold, rvec f, rvec step) | |||
| 631 | { | |||
| 632 | real xo, yo, zo; | |||
| 633 | real dx, dy, dz; | |||
| 634 | ||||
| 635 | xo = xold[XX0]; | |||
| 636 | yo = xold[YY1]; | |||
| 637 | zo = xold[ZZ2]; | |||
| 638 | ||||
| 639 | dx = f[XX0]*step[XX0]; | |||
| 640 | dy = f[YY1]*step[YY1]; | |||
| 641 | dz = f[ZZ2]*step[ZZ2]; | |||
| 642 | ||||
| 643 | xnew[XX0] = xo+dx; | |||
| 644 | xnew[YY1] = yo+dy; | |||
| 645 | xnew[ZZ2] = zo+dz; | |||
| 646 | } | |||
| 647 | ||||
| 648 | static void directional_sd(rvec xold[], rvec xnew[], rvec acc_dir[], | |||
| 649 | int start, int homenr, real step) | |||
| 650 | { | |||
| 651 | int i; | |||
| 652 | ||||
| 653 | for (i = start; i < homenr; i++) | |||
| 654 | { | |||
| 655 | do_1pos(xnew[i], xold[i], acc_dir[i], step); | |||
| 656 | } | |||
| 657 | } | |||
| 658 | ||||
| 659 | static void shell_pos_sd(rvec xcur[], rvec xnew[], rvec f[], | |||
| 660 | int ns, t_shell s[], int count) | |||
| 661 | { | |||
| 662 | const real step_scale_min = 0.8, | |||
| 663 | step_scale_increment = 0.2, | |||
| 664 | step_scale_max = 1.2, | |||
| 665 | step_scale_multiple = (step_scale_max - step_scale_min) / step_scale_increment; | |||
| 666 | int i, shell, d; | |||
| 667 | real dx, df, k_est; | |||
| 668 | #ifdef PRINT_STEP | |||
| 669 | real step_min, step_max; | |||
| 670 | ||||
| 671 | step_min = 1e30; | |||
| 672 | step_max = 0; | |||
| 673 | #endif | |||
| 674 | for (i = 0; (i < ns); i++) | |||
| 675 | { | |||
| 676 | shell = s[i].shell; | |||
| 677 | if (count == 1) | |||
| 678 | { | |||
| 679 | for (d = 0; d < DIM3; d++) | |||
| 680 | { | |||
| 681 | s[i].step[d] = s[i].k_1; | |||
| 682 | #ifdef PRINT_STEP | |||
| 683 | step_min = min(step_min, s[i].step[d])(((step_min) < (s[i].step[d])) ? (step_min) : (s[i].step[d ]) ); | |||
| 684 | step_max = max(step_max, s[i].step[d])(((step_max) > (s[i].step[d])) ? (step_max) : (s[i].step[d ]) ); | |||
| 685 | #endif | |||
| 686 | } | |||
| 687 | } | |||
| 688 | else | |||
| 689 | { | |||
| 690 | for (d = 0; d < DIM3; d++) | |||
| 691 | { | |||
| 692 | dx = xcur[shell][d] - s[i].xold[d]; | |||
| 693 | df = f[shell][d] - s[i].fold[d]; | |||
| 694 | /* -dx/df gets used to generate an interpolated value, but would | |||
| 695 | * cause a NaN if df were binary-equal to zero. Values close to | |||
| 696 | * zero won't cause problems (because of the min() and max()), so | |||
| 697 | * just testing for binary inequality is OK. */ | |||
| 698 | if (0.0 != df) | |||
| 699 | { | |||
| 700 | k_est = -dx/df; | |||
| 701 | /* Scale the step size by a factor interpolated from | |||
| 702 | * step_scale_min to step_scale_max, as k_est goes from 0 to | |||
| 703 | * step_scale_multiple * s[i].step[d] */ | |||
| 704 | s[i].step[d] = | |||
| 705 | step_scale_min * s[i].step[d] + | |||
| 706 | step_scale_increment * min(step_scale_multiple * s[i].step[d], max(k_est, 0))(((step_scale_multiple * s[i].step[d]) < ((((k_est) > ( 0)) ? (k_est) : (0) ))) ? (step_scale_multiple * s[i].step[d] ) : ((((k_est) > (0)) ? (k_est) : (0) )) ); | |||
| 707 | } | |||
| 708 | else | |||
| 709 | { | |||
| 710 | /* Here 0 == df */ | |||
| 711 | if (gmx_numzero(dx)) /* 0 == dx */ | |||
| 712 | { | |||
| 713 | /* Likely this will never happen, but if it does just | |||
| 714 | * don't scale the step. */ | |||
| 715 | } | |||
| 716 | else /* 0 != dx */ | |||
| 717 | { | |||
| 718 | s[i].step[d] *= step_scale_max; | |||
| 719 | } | |||
| 720 | } | |||
| 721 | #ifdef PRINT_STEP | |||
| 722 | step_min = min(step_min, s[i].step[d])(((step_min) < (s[i].step[d])) ? (step_min) : (s[i].step[d ]) ); | |||
| 723 | step_max = max(step_max, s[i].step[d])(((step_max) > (s[i].step[d])) ? (step_max) : (s[i].step[d ]) ); | |||
| 724 | #endif | |||
| 725 | } | |||
| 726 | } | |||
| 727 | copy_rvec(xcur[shell], s[i].xold); | |||
| 728 | copy_rvec(f[shell], s[i].fold); | |||
| 729 | ||||
| 730 | do_1pos3(xnew[shell], xcur[shell], f[shell], s[i].step); | |||
| 731 | ||||
| 732 | if (gmx_debug_at) | |||
| 733 | { | |||
| 734 | fprintf(debug, "shell[%d] = %d\n", i, shell); | |||
| 735 | pr_rvec(debug, 0, "fshell", f[shell], DIM3, TRUE1); | |||
| 736 | pr_rvec(debug, 0, "xold", xcur[shell], DIM3, TRUE1); | |||
| 737 | pr_rvec(debug, 0, "step", s[i].step, DIM3, TRUE1); | |||
| 738 | pr_rvec(debug, 0, "xnew", xnew[shell], DIM3, TRUE1); | |||
| 739 | } | |||
| 740 | } | |||
| 741 | #ifdef PRINT_STEP | |||
| 742 | printf("step %.3e %.3e\n", step_min, step_max); | |||
| 743 | #endif | |||
| 744 | } | |||
| 745 | ||||
| 746 | static void decrease_step_size(int nshell, t_shell s[]) | |||
| 747 | { | |||
| 748 | int i; | |||
| 749 | ||||
| 750 | for (i = 0; i < nshell; i++) | |||
| 751 | { | |||
| 752 | svmul(0.8, s[i].step, s[i].step); | |||
| 753 | } | |||
| 754 | } | |||
| 755 | ||||
| 756 | static void print_epot(FILE *fp, gmx_int64_t mdstep, int count, real epot, real df, | |||
| 757 | int ndir, real sf_dir) | |||
| 758 | { | |||
| 759 | char buf[22]; | |||
| 760 | ||||
| 761 | fprintf(fp, "MDStep=%5s/%2d EPot: %12.8e, rmsF: %6.2e", | |||
| 762 | gmx_step_str(mdstep, buf), count, epot, df); | |||
| 763 | if (ndir) | |||
| 764 | { | |||
| 765 | fprintf(fp, ", dir. rmsF: %6.2e\n", sqrt(sf_dir/ndir)); | |||
| 766 | } | |||
| 767 | else | |||
| 768 | { | |||
| 769 | fprintf(fp, "\n"); | |||
| 770 | } | |||
| 771 | } | |||
| 772 | ||||
| 773 | ||||
| 774 | static real rms_force(t_commrec *cr, rvec f[], int ns, t_shell s[], | |||
| 775 | int ndir, real *sf_dir, real *Epot) | |||
| 776 | { | |||
| 777 | int i, shell, ntot; | |||
| 778 | double buf[4]; | |||
| 779 | ||||
| 780 | buf[0] = *sf_dir; | |||
| 781 | for (i = 0; i < ns; i++) | |||
| 782 | { | |||
| 783 | shell = s[i].shell; | |||
| 784 | buf[0] += norm2(f[shell]); | |||
| 785 | } | |||
| 786 | ntot = ns; | |||
| 787 | ||||
| 788 | if (PAR(cr)((cr)->nnodes > 1)) | |||
| 789 | { | |||
| 790 | buf[1] = ntot; | |||
| 791 | buf[2] = *sf_dir; | |||
| 792 | buf[3] = *Epot; | |||
| 793 | gmx_sumd(4, buf, cr); | |||
| 794 | ntot = (int)(buf[1] + 0.5); | |||
| 795 | *sf_dir = buf[2]; | |||
| 796 | *Epot = buf[3]; | |||
| 797 | } | |||
| 798 | ntot += ndir; | |||
| 799 | ||||
| 800 | return (ntot ? sqrt(buf[0]/ntot) : 0); | |||
| 801 | } | |||
| 802 | ||||
| 803 | static void check_pbc(FILE *fp, rvec x[], int shell) | |||
| 804 | { | |||
| 805 | int m, now; | |||
| 806 | ||||
| 807 | now = shell-4; | |||
| 808 | for (m = 0; (m < DIM3); m++) | |||
| 809 | { | |||
| 810 | if (fabs(x[shell][m]-x[now][m]) > 0.3) | |||
| 811 | { | |||
| 812 | pr_rvecs(fp, 0, "SHELL-X", x+now, 5); | |||
| 813 | break; | |||
| 814 | } | |||
| 815 | } | |||
| 816 | } | |||
| 817 | ||||
| 818 | static void dump_shells(FILE *fp, rvec x[], rvec f[], real ftol, int ns, t_shell s[]) | |||
| 819 | { | |||
| 820 | int i, shell; | |||
| 821 | real ft2, ff2; | |||
| 822 | ||||
| 823 | ft2 = sqr(ftol); | |||
| 824 | ||||
| 825 | for (i = 0; (i < ns); i++) | |||
| 826 | { | |||
| 827 | shell = s[i].shell; | |||
| 828 | ff2 = iprod(f[shell], f[shell]); | |||
| 829 | if (ff2 > ft2) | |||
| 830 | { | |||
| 831 | fprintf(fp, "SHELL %5d, force %10.5f %10.5f %10.5f, |f| %10.5f\n", | |||
| 832 | shell, f[shell][XX0], f[shell][YY1], f[shell][ZZ2], sqrt(ff2)); | |||
| 833 | } | |||
| 834 | check_pbc(fp, x, shell); | |||
| 835 | } | |||
| 836 | } | |||
| 837 | ||||
| 838 | static void init_adir(FILE *log, gmx_shellfc_t shfc, | |||
| 839 | gmx_constr_t constr, t_idef *idef, t_inputrec *ir, | |||
| 840 | t_commrec *cr, int dd_ac1, | |||
| 841 | gmx_int64_t step, t_mdatoms *md, int start, int end, | |||
| 842 | rvec *x_old, rvec *x_init, rvec *x, | |||
| 843 | rvec *f, rvec *acc_dir, | |||
| 844 | gmx_bool bMolPBC, matrix box, | |||
| 845 | real *lambda, real *dvdlambda, t_nrnb *nrnb) | |||
| 846 | { | |||
| 847 | rvec *xnold, *xnew; | |||
| 848 | double w_dt; | |||
| 849 | int gf, ga, gt; | |||
| 850 | real dt, scale; | |||
| 851 | int n, d; | |||
| 852 | unsigned short *ptype; | |||
| 853 | rvec p, dx; | |||
| 854 | ||||
| 855 | if (DOMAINDECOMP(cr)(((cr)->dd != ((void*)0)) && ((cr)->nnodes > 1))) | |||
| 856 | { | |||
| 857 | n = dd_ac1; | |||
| 858 | } | |||
| 859 | else | |||
| 860 | { | |||
| 861 | n = end - start; | |||
| 862 | } | |||
| 863 | if (n > shfc->adir_nalloc) | |||
| 864 | { | |||
| 865 | shfc->adir_nalloc = over_alloc_dd(n); | |||
| 866 | srenew(shfc->adir_xnold, shfc->adir_nalloc)(shfc->adir_xnold) = save_realloc("shfc->adir_xnold", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/shellfc.c" , 866, (shfc->adir_xnold), (shfc->adir_nalloc), sizeof( *(shfc->adir_xnold))); | |||
| 867 | srenew(shfc->adir_xnew, shfc->adir_nalloc)(shfc->adir_xnew) = save_realloc("shfc->adir_xnew", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/shellfc.c" , 867, (shfc->adir_xnew), (shfc->adir_nalloc), sizeof(* (shfc->adir_xnew))); | |||
| 868 | } | |||
| 869 | xnold = shfc->adir_xnold; | |||
| 870 | xnew = shfc->adir_xnew; | |||
| 871 | ||||
| 872 | ptype = md->ptype; | |||
| 873 | ||||
| 874 | dt = ir->delta_t; | |||
| 875 | ||||
| 876 | /* Does NOT work with freeze or acceleration groups (yet) */ | |||
| 877 | for (n = start; n < end; n++) | |||
| 878 | { | |||
| 879 | w_dt = md->invmass[n]*dt; | |||
| 880 | ||||
| 881 | for (d = 0; d < DIM3; d++) | |||
| 882 | { | |||
| 883 | if ((ptype[n] != eptVSite) && (ptype[n] != eptShell)) | |||
| 884 | { | |||
| 885 | xnold[n-start][d] = x[n][d] - (x_init[n][d] - x_old[n][d]); | |||
| 886 | xnew[n-start][d] = 2*x[n][d] - x_old[n][d] + f[n][d]*w_dt*dt; | |||
| 887 | } | |||
| 888 | else | |||
| 889 | { | |||
| 890 | xnold[n-start][d] = x[n][d]; | |||
| 891 | xnew[n-start][d] = x[n][d]; | |||
| 892 | } | |||
| 893 | } | |||
| 894 | } | |||
| 895 | constrain(log, FALSE0, FALSE0, constr, idef, ir, NULL((void*)0), cr, step, 0, md, | |||
| 896 | x, xnold-start, NULL((void*)0), bMolPBC, box, | |||
| 897 | lambda[efptBONDED], &(dvdlambda[efptBONDED]), | |||
| 898 | NULL((void*)0), NULL((void*)0), nrnb, econqCoord, FALSE0, 0, 0); | |||
| 899 | constrain(log, FALSE0, FALSE0, constr, idef, ir, NULL((void*)0), cr, step, 0, md, | |||
| 900 | x, xnew-start, NULL((void*)0), bMolPBC, box, | |||
| 901 | lambda[efptBONDED], &(dvdlambda[efptBONDED]), | |||
| 902 | NULL((void*)0), NULL((void*)0), nrnb, econqCoord, FALSE0, 0, 0); | |||
| 903 | ||||
| 904 | for (n = start; n < end; n++) | |||
| 905 | { | |||
| 906 | for (d = 0; d < DIM3; d++) | |||
| 907 | { | |||
| 908 | xnew[n-start][d] = | |||
| 909 | -(2*x[n][d]-xnold[n-start][d]-xnew[n-start][d])/sqr(dt) | |||
| 910 | - f[n][d]*md->invmass[n]; | |||
| 911 | } | |||
| 912 | clear_rvec(acc_dir[n]); | |||
| 913 | } | |||
| 914 | ||||
| 915 | /* Project the acceleration on the old bond directions */ | |||
| 916 | constrain(log, FALSE0, FALSE0, constr, idef, ir, NULL((void*)0), cr, step, 0, md, | |||
| 917 | x_old, xnew-start, acc_dir, bMolPBC, box, | |||
| 918 | lambda[efptBONDED], &(dvdlambda[efptBONDED]), | |||
| 919 | NULL((void*)0), NULL((void*)0), nrnb, econqDeriv_FlexCon, FALSE0, 0, 0); | |||
| 920 | } | |||
| 921 | ||||
| 922 | int relax_shell_flexcon(FILE *fplog, t_commrec *cr, gmx_bool bVerbose, | |||
| 923 | gmx_int64_t mdstep, t_inputrec *inputrec, | |||
| 924 | gmx_bool bDoNS, int force_flags, | |||
| 925 | gmx_localtop_t *top, | |||
| 926 | gmx_constr_t constr, | |||
| 927 | gmx_enerdata_t *enerd, t_fcdata *fcd, | |||
| 928 | t_state *state, rvec f[], | |||
| 929 | tensor force_vir, | |||
| 930 | t_mdatoms *md, | |||
| 931 | t_nrnb *nrnb, gmx_wallcycle_t wcycle, | |||
| 932 | t_graph *graph, | |||
| 933 | gmx_groups_t *groups, | |||
| 934 | struct gmx_shellfc *shfc, | |||
| 935 | t_forcerec *fr, | |||
| 936 | gmx_bool bBornRadii, | |||
| 937 | double t, rvec mu_tot, | |||
| 938 | gmx_bool *bConverged, | |||
| 939 | gmx_vsite_t *vsite, | |||
| 940 | FILE *fp_field) | |||
| 941 | { | |||
| 942 | int nshell; | |||
| 943 | t_shell *shell; | |||
| 944 | t_idef *idef; | |||
| 945 | rvec *pos[2], *force[2], *acc_dir = NULL((void*)0), *x_old = NULL((void*)0); | |||
| 946 | real Epot[2], df[2]; | |||
| 947 | rvec dx; | |||
| 948 | real sf_dir, invdt; | |||
| 949 | real ftol, xiH, xiS, dum = 0; | |||
| 950 | char sbuf[22]; | |||
| 951 | gmx_bool bCont, bInit; | |||
| 952 | int nat, dd_ac0, dd_ac1 = 0, i; | |||
| 953 | int start = 0, homenr = md->homenr, end = start+homenr, cg0, cg1; | |||
| 954 | int nflexcon, g, number_steps, d, Min = 0, count = 0; | |||
| 955 | #define Try(1-Min) (1-Min) /* At start Try = 1 */ | |||
| 956 | ||||
| 957 | bCont = (mdstep == inputrec->init_step) && inputrec->bContinuation; | |||
| 958 | bInit = (mdstep == inputrec->init_step) || shfc->bRequireInit; | |||
| 959 | ftol = inputrec->em_tol; | |||
| 960 | number_steps = inputrec->niter; | |||
| 961 | nshell = shfc->nshell; | |||
| 962 | shell = shfc->shell; | |||
| 963 | nflexcon = shfc->nflexcon; | |||
| 964 | ||||
| 965 | idef = &top->idef; | |||
| 966 | ||||
| 967 | if (DOMAINDECOMP(cr)(((cr)->dd != ((void*)0)) && ((cr)->nnodes > 1))) | |||
| 968 | { | |||
| 969 | nat = dd_natoms_vsite(cr->dd); | |||
| 970 | if (nflexcon > 0) | |||
| 971 | { | |||
| 972 | dd_get_constraint_range(cr->dd, &dd_ac0, &dd_ac1); | |||
| 973 | nat = max(nat, dd_ac1)(((nat) > (dd_ac1)) ? (nat) : (dd_ac1) ); | |||
| 974 | } | |||
| 975 | } | |||
| 976 | else | |||
| 977 | { | |||
| 978 | nat = state->natoms; | |||
| 979 | } | |||
| 980 | ||||
| 981 | if (nat > shfc->x_nalloc) | |||
| 982 | { | |||
| 983 | /* Allocate local arrays */ | |||
| 984 | shfc->x_nalloc = over_alloc_dd(nat); | |||
| 985 | for (i = 0; (i < 2); i++) | |||
| 986 | { | |||
| 987 | srenew(shfc->x[i], shfc->x_nalloc)(shfc->x[i]) = save_realloc("shfc->x[i]", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/shellfc.c" , 987, (shfc->x[i]), (shfc->x_nalloc), sizeof(*(shfc-> x[i]))); | |||
| 988 | srenew(shfc->f[i], shfc->x_nalloc)(shfc->f[i]) = save_realloc("shfc->f[i]", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/shellfc.c" , 988, (shfc->f[i]), (shfc->x_nalloc), sizeof(*(shfc-> f[i]))); | |||
| 989 | } | |||
| 990 | } | |||
| 991 | for (i = 0; (i < 2); i++) | |||
| 992 | { | |||
| 993 | pos[i] = shfc->x[i]; | |||
| 994 | force[i] = shfc->f[i]; | |||
| 995 | } | |||
| 996 | ||||
| 997 | /* When we had particle decomposition, this code only worked with | |||
| 998 | * PD when all particles involved with each shell were in the same | |||
| 999 | * charge group. Not sure if this is still relevant. */ | |||
| 1000 | if (bDoNS && inputrec->ePBC != epbcNONE && !DOMAINDECOMP(cr)(((cr)->dd != ((void*)0)) && ((cr)->nnodes > 1))) | |||
| 1001 | { | |||
| 1002 | /* This is the only time where the coordinates are used | |||
| 1003 | * before do_force is called, which normally puts all | |||
| 1004 | * charge groups in the box. | |||
| 1005 | */ | |||
| 1006 | cg0 = 0; | |||
| 1007 | cg1 = top->cgs.nr; | |||
| 1008 | put_charge_groups_in_box(fplog, cg0, cg1, fr->ePBC, state->box, | |||
| 1009 | &(top->cgs), state->x, fr->cg_cm); | |||
| 1010 | if (graph) | |||
| 1011 | { | |||
| 1012 | mk_mshift(fplog, graph, fr->ePBC, state->box, state->x); | |||
| 1013 | } | |||
| 1014 | } | |||
| 1015 | ||||
| 1016 | /* After this all coordinate arrays will contain whole molecules */ | |||
| 1017 | if (graph) | |||
| 1018 | { | |||
| 1019 | shift_self(graph, state->box, state->x); | |||
| 1020 | } | |||
| 1021 | ||||
| 1022 | if (nflexcon) | |||
| 1023 | { | |||
| 1024 | if (nat > shfc->flex_nalloc) | |||
| 1025 | { | |||
| 1026 | shfc->flex_nalloc = over_alloc_dd(nat); | |||
| 1027 | srenew(shfc->acc_dir, shfc->flex_nalloc)(shfc->acc_dir) = save_realloc("shfc->acc_dir", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/shellfc.c" , 1027, (shfc->acc_dir), (shfc->flex_nalloc), sizeof(*( shfc->acc_dir))); | |||
| 1028 | srenew(shfc->x_old, shfc->flex_nalloc)(shfc->x_old) = save_realloc("shfc->x_old", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/shellfc.c" , 1028, (shfc->x_old), (shfc->flex_nalloc), sizeof(*(shfc ->x_old))); | |||
| 1029 | } | |||
| 1030 | acc_dir = shfc->acc_dir; | |||
| 1031 | x_old = shfc->x_old; | |||
| 1032 | for (i = 0; i < homenr; i++) | |||
| 1033 | { | |||
| 1034 | for (d = 0; d < DIM3; d++) | |||
| 1035 | { | |||
| 1036 | shfc->x_old[i][d] = | |||
| 1037 | state->x[start+i][d] - state->v[start+i][d]*inputrec->delta_t; | |||
| 1038 | } | |||
| 1039 | } | |||
| 1040 | } | |||
| 1041 | ||||
| 1042 | /* Do a prediction of the shell positions */ | |||
| 1043 | if (shfc->bPredict && !bCont) | |||
| 1044 | { | |||
| 1045 | predict_shells(fplog, state->x, state->v, inputrec->delta_t, nshell, shell, | |||
| 1046 | md->massT, NULL((void*)0), bInit); | |||
| 1047 | } | |||
| 1048 | ||||
| 1049 | /* do_force expected the charge groups to be in the box */ | |||
| 1050 | if (graph) | |||
| 1051 | { | |||
| 1052 | unshift_self(graph, state->box, state->x); | |||
| 1053 | } | |||
| 1054 | ||||
| 1055 | /* Calculate the forces first time around */ | |||
| 1056 | if (gmx_debug_at) | |||
| 1057 | { | |||
| 1058 | pr_rvecs(debug, 0, "x b4 do_force", state->x + start, homenr); | |||
| 1059 | } | |||
| 1060 | do_force(fplog, cr, inputrec, mdstep, nrnb, wcycle, top, groups, | |||
| 1061 | state->box, state->x, &state->hist, | |||
| 1062 | force[Min], force_vir, md, enerd, fcd, | |||
| 1063 | state->lambda, graph, | |||
| 1064 | fr, vsite, mu_tot, t, fp_field, NULL((void*)0), bBornRadii, | |||
| 1065 | (bDoNS ? GMX_FORCE_NS(1<<2) : 0) | force_flags); | |||
| 1066 | ||||
| 1067 | sf_dir = 0; | |||
| 1068 | if (nflexcon) | |||
| 1069 | { | |||
| 1070 | init_adir(fplog, shfc, | |||
| 1071 | constr, idef, inputrec, cr, dd_ac1, mdstep, md, start, end, | |||
| 1072 | shfc->x_old-start, state->x, state->x, force[Min], | |||
| 1073 | shfc->acc_dir-start, | |||
| 1074 | fr->bMolPBC, state->box, state->lambda, &dum, nrnb); | |||
| 1075 | ||||
| 1076 | for (i = start; i < end; i++) | |||
| 1077 | { | |||
| 1078 | sf_dir += md->massT[i]*norm2(shfc->acc_dir[i-start]); | |||
| 1079 | } | |||
| 1080 | } | |||
| 1081 | ||||
| 1082 | Epot[Min] = enerd->term[F_EPOT]; | |||
| 1083 | ||||
| 1084 | df[Min] = rms_force(cr, shfc->f[Min], nshell, shell, nflexcon, &sf_dir, &Epot[Min]); | |||
| 1085 | df[Try(1-Min)] = 0; | |||
| 1086 | if (debug) | |||
| 1087 | { | |||
| 1088 | fprintf(debug, "df = %g %g\n", df[Min], df[Try(1-Min)]); | |||
| 1089 | } | |||
| 1090 | ||||
| 1091 | if (gmx_debug_at) | |||
| 1092 | { | |||
| 1093 | pr_rvecs(debug, 0, "force0", force[Min], md->nr); | |||
| 1094 | } | |||
| 1095 | ||||
| 1096 | if (nshell+nflexcon > 0) | |||
| 1097 | { | |||
| 1098 | /* Copy x to pos[Min] & pos[Try]: during minimization only the | |||
| 1099 | * shell positions are updated, therefore the other particles must | |||
| 1100 | * be set here. | |||
| 1101 | */ | |||
| 1102 | memcpy(pos[Min], state->x, nat*sizeof(state->x[0])); | |||
| 1103 | memcpy(pos[Try(1-Min)], state->x, nat*sizeof(state->x[0])); | |||
| 1104 | } | |||
| 1105 | ||||
| 1106 | if (bVerbose && MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) | |||
| 1107 | { | |||
| 1108 | print_epot(stdoutstdout, mdstep, 0, Epot[Min], df[Min], nflexcon, sf_dir); | |||
| 1109 | } | |||
| 1110 | ||||
| 1111 | if (debug) | |||
| 1112 | { | |||
| 1113 | fprintf(debug, "%17s: %14.10e\n", | |||
| 1114 | interaction_function[F_EKIN].longname, enerd->term[F_EKIN]); | |||
| 1115 | fprintf(debug, "%17s: %14.10e\n", | |||
| 1116 | interaction_function[F_EPOT].longname, enerd->term[F_EPOT]); | |||
| 1117 | fprintf(debug, "%17s: %14.10e\n", | |||
| 1118 | interaction_function[F_ETOT].longname, enerd->term[F_ETOT]); | |||
| 1119 | fprintf(debug, "SHELLSTEP %s\n", gmx_step_str(mdstep, sbuf)); | |||
| 1120 | } | |||
| 1121 | ||||
| 1122 | /* First check whether we should do shells, or whether the force is | |||
| 1123 | * low enough even without minimization. | |||
| 1124 | */ | |||
| 1125 | *bConverged = (df[Min] < ftol); | |||
| 1126 | ||||
| 1127 | for (count = 1; (!(*bConverged) && (count < number_steps)); count++) | |||
| 1128 | { | |||
| 1129 | if (vsite) | |||
| 1130 | { | |||
| 1131 | construct_vsites(vsite, pos[Min], inputrec->delta_t, state->v, | |||
| 1132 | idef->iparams, idef->il, | |||
| 1133 | fr->ePBC, fr->bMolPBC, cr, state->box); | |||
| 1134 | } | |||
| 1135 | ||||
| 1136 | if (nflexcon) | |||
| 1137 | { | |||
| 1138 | init_adir(fplog, shfc, | |||
| 1139 | constr, idef, inputrec, cr, dd_ac1, mdstep, md, start, end, | |||
| 1140 | x_old-start, state->x, pos[Min], force[Min], acc_dir-start, | |||
| 1141 | fr->bMolPBC, state->box, state->lambda, &dum, nrnb); | |||
| 1142 | ||||
| 1143 | directional_sd(pos[Min], pos[Try(1-Min)], acc_dir-start, start, end, | |||
| 1144 | fr->fc_stepsize); | |||
| 1145 | } | |||
| 1146 | ||||
| 1147 | /* New positions, Steepest descent */ | |||
| 1148 | shell_pos_sd(pos[Min], pos[Try(1-Min)], force[Min], nshell, shell, count); | |||
| 1149 | ||||
| 1150 | /* do_force expected the charge groups to be in the box */ | |||
| 1151 | if (graph) | |||
| 1152 | { | |||
| 1153 | unshift_self(graph, state->box, pos[Try(1-Min)]); | |||
| 1154 | } | |||
| 1155 | ||||
| 1156 | if (gmx_debug_at) | |||
| 1157 | { | |||
| 1158 | pr_rvecs(debug, 0, "RELAX: pos[Min] ", pos[Min] + start, homenr); | |||
| 1159 | pr_rvecs(debug, 0, "RELAX: pos[Try] ", pos[Try(1-Min)] + start, homenr); | |||
| 1160 | } | |||
| 1161 | /* Try the new positions */ | |||
| 1162 | do_force(fplog, cr, inputrec, 1, nrnb, wcycle, | |||
| 1163 | top, groups, state->box, pos[Try(1-Min)], &state->hist, | |||
| 1164 | force[Try(1-Min)], force_vir, | |||
| 1165 | md, enerd, fcd, state->lambda, graph, | |||
| 1166 | fr, vsite, mu_tot, t, fp_field, NULL((void*)0), bBornRadii, | |||
| 1167 | force_flags); | |||
| 1168 | ||||
| 1169 | if (gmx_debug_at) | |||
| 1170 | { | |||
| 1171 | pr_rvecs(debug, 0, "RELAX: force[Min]", force[Min] + start, homenr); | |||
| 1172 | pr_rvecs(debug, 0, "RELAX: force[Try]", force[Try(1-Min)] + start, homenr); | |||
| 1173 | } | |||
| 1174 | sf_dir = 0; | |||
| 1175 | if (nflexcon) | |||
| 1176 | { | |||
| 1177 | init_adir(fplog, shfc, | |||
| 1178 | constr, idef, inputrec, cr, dd_ac1, mdstep, md, start, end, | |||
| 1179 | x_old-start, state->x, pos[Try(1-Min)], force[Try(1-Min)], acc_dir-start, | |||
| 1180 | fr->bMolPBC, state->box, state->lambda, &dum, nrnb); | |||
| 1181 | ||||
| 1182 | for (i = start; i < end; i++) | |||
| 1183 | { | |||
| 1184 | sf_dir += md->massT[i]*norm2(acc_dir[i-start]); | |||
| 1185 | } | |||
| 1186 | } | |||
| 1187 | ||||
| 1188 | Epot[Try(1-Min)] = enerd->term[F_EPOT]; | |||
| 1189 | ||||
| 1190 | df[Try(1-Min)] = rms_force(cr, force[Try(1-Min)], nshell, shell, nflexcon, &sf_dir, &Epot[Try(1-Min)]); | |||
| 1191 | ||||
| 1192 | if (debug) | |||
| 1193 | { | |||
| 1194 | fprintf(debug, "df = %g %g\n", df[Min], df[Try(1-Min)]); | |||
| 1195 | } | |||
| 1196 | ||||
| 1197 | if (debug) | |||
| 1198 | { | |||
| 1199 | if (gmx_debug_at) | |||
| 1200 | { | |||
| 1201 | pr_rvecs(debug, 0, "F na do_force", force[Try(1-Min)] + start, homenr); | |||
| 1202 | } | |||
| 1203 | if (gmx_debug_at) | |||
| 1204 | { | |||
| 1205 | fprintf(debug, "SHELL ITER %d\n", count); | |||
| 1206 | dump_shells(debug, pos[Try(1-Min)], force[Try(1-Min)], ftol, nshell, shell); | |||
| 1207 | } | |||
| 1208 | } | |||
| 1209 | ||||
| 1210 | if (bVerbose && MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) | |||
| 1211 | { | |||
| 1212 | print_epot(stdoutstdout, mdstep, count, Epot[Try(1-Min)], df[Try(1-Min)], nflexcon, sf_dir); | |||
| 1213 | } | |||
| 1214 | ||||
| 1215 | *bConverged = (df[Try(1-Min)] < ftol); | |||
| 1216 | ||||
| 1217 | if ((df[Try(1-Min)] < df[Min])) | |||
| 1218 | { | |||
| 1219 | if (debug) | |||
| 1220 | { | |||
| 1221 | fprintf(debug, "Swapping Min and Try\n"); | |||
| 1222 | } | |||
| 1223 | if (nflexcon) | |||
| 1224 | { | |||
| 1225 | /* Correct the velocities for the flexible constraints */ | |||
| 1226 | invdt = 1/inputrec->delta_t; | |||
| 1227 | for (i = start; i < end; i++) | |||
| 1228 | { | |||
| 1229 | for (d = 0; d < DIM3; d++) | |||
| 1230 | { | |||
| 1231 | state->v[i][d] += (pos[Try(1-Min)][i][d] - pos[Min][i][d])*invdt; | |||
| 1232 | } | |||
| 1233 | } | |||
| 1234 | } | |||
| 1235 | Min = Try(1-Min); | |||
| 1236 | } | |||
| 1237 | else | |||
| 1238 | { | |||
| 1239 | decrease_step_size(nshell, shell); | |||
| 1240 | } | |||
| 1241 | } | |||
| 1242 | if (MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1)) && !(*bConverged)) | |||
| 1243 | { | |||
| 1244 | /* Note that the energies and virial are incorrect when not converged */ | |||
| 1245 | if (fplog) | |||
| 1246 | { | |||
| 1247 | fprintf(fplog, | |||
| 1248 | "step %s: EM did not converge in %d iterations, RMS force %.3f\n", | |||
| 1249 | gmx_step_str(mdstep, sbuf), number_steps, df[Min]); | |||
| 1250 | } | |||
| 1251 | fprintf(stderrstderr, | |||
| 1252 | "step %s: EM did not converge in %d iterations, RMS force %.3f\n", | |||
| 1253 | gmx_step_str(mdstep, sbuf), number_steps, df[Min]); | |||
| 1254 | } | |||
| 1255 | ||||
| 1256 | /* Copy back the coordinates and the forces */ | |||
| 1257 | memcpy(state->x, pos[Min], nat*sizeof(state->x[0])); | |||
| 1258 | memcpy(f, force[Min], nat*sizeof(f[0])); | |||
| 1259 | ||||
| 1260 | return count; | |||
| 1261 | } |