| File: | gromacs/mdlib/clincs.c |
| Location: | line 1246, column 9 |
| Description: | Value stored to 'bLocal' is never read |
| 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-2004, 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 | /* This file is completely threadsafe - keep it that way! */ |
| 38 | #ifdef HAVE_CONFIG_H1 |
| 39 | #include <config.h> |
| 40 | #endif |
| 41 | |
| 42 | #include <math.h> |
| 43 | #include <stdlib.h> |
| 44 | |
| 45 | #include "types/commrec.h" |
| 46 | #include "constr.h" |
| 47 | #include "copyrite.h" |
| 48 | #include "physics.h" |
| 49 | #include "gromacs/math/vec.h" |
| 50 | #include "pbc.h" |
| 51 | #include "mdrun.h" |
| 52 | #include "nrnb.h" |
| 53 | #include "domdec.h" |
| 54 | #include "mtop_util.h" |
| 55 | #include "gmx_omp_nthreads.h" |
| 56 | |
| 57 | #include "gromacs/fileio/gmxfio.h" |
| 58 | #include "gromacs/utility/fatalerror.h" |
| 59 | #include "gromacs/utility/gmxomp.h" |
| 60 | #include "gromacs/utility/smalloc.h" |
| 61 | |
| 62 | typedef struct { |
| 63 | int b0; /* first constraint for this thread */ |
| 64 | int b1; /* b1-1 is the last constraint for this thread */ |
| 65 | int nind; /* number of indices */ |
| 66 | int *ind; /* constraint index for updating atom data */ |
| 67 | int nind_r; /* number of indices */ |
| 68 | int *ind_r; /* constraint index for updating atom data */ |
| 69 | int ind_nalloc; /* allocation size of ind and ind_r */ |
| 70 | tensor vir_r_m_dr; /* temporary variable for virial calculation */ |
| 71 | } lincs_thread_t; |
| 72 | |
| 73 | typedef struct gmx_lincsdata { |
| 74 | int ncg; /* the global number of constraints */ |
| 75 | int ncg_flex; /* the global number of flexible constraints */ |
| 76 | int ncg_triangle; /* the global number of constraints in triangles */ |
| 77 | int nIter; /* the number of iterations */ |
| 78 | int nOrder; /* the order of the matrix expansion */ |
| 79 | int nc; /* the number of constraints */ |
| 80 | int nc_alloc; /* the number we allocated memory for */ |
| 81 | int ncc; /* the number of constraint connections */ |
| 82 | int ncc_alloc; /* the number we allocated memory for */ |
| 83 | real matlam; /* the FE lambda value used for filling blc and blmf */ |
| 84 | real *bllen0; /* the reference distance in topology A */ |
| 85 | real *ddist; /* the reference distance in top B - the r.d. in top A */ |
| 86 | int *bla; /* the atom pairs involved in the constraints */ |
| 87 | real *blc; /* 1/sqrt(invmass1 + invmass2) */ |
| 88 | real *blc1; /* as blc, but with all masses 1 */ |
| 89 | int *blnr; /* index into blbnb and blmf */ |
| 90 | int *blbnb; /* list of constraint connections */ |
| 91 | int ntriangle; /* the local number of constraints in triangles */ |
| 92 | int *triangle; /* the list of triangle constraints */ |
| 93 | int *tri_bits; /* the bits tell if the matrix element should be used */ |
| 94 | int ncc_triangle; /* the number of constraint connections in triangles */ |
| 95 | gmx_bool bCommIter; /* communicate before each LINCS interation */ |
| 96 | real *blmf; /* matrix of mass factors for constraint connections */ |
| 97 | real *blmf1; /* as blmf, but with all masses 1 */ |
| 98 | real *bllen; /* the reference bond length */ |
| 99 | int nth; /* The number of threads doing LINCS */ |
| 100 | lincs_thread_t *th; /* LINCS thread division */ |
| 101 | unsigned *atf; /* atom flags for thread parallelization */ |
| 102 | int atf_nalloc; /* allocation size of atf */ |
| 103 | /* arrays for temporary storage in the LINCS algorithm */ |
| 104 | rvec *tmpv; |
| 105 | real *tmpncc; |
| 106 | real *tmp1; |
| 107 | real *tmp2; |
| 108 | real *tmp3; |
| 109 | real *tmp4; |
| 110 | real *mlambda; /* the Lagrange multipliers * -1 */ |
| 111 | /* storage for the constraint RMS relative deviation output */ |
| 112 | real rmsd_data[3]; |
| 113 | } t_gmx_lincsdata; |
| 114 | |
| 115 | real *lincs_rmsd_data(struct gmx_lincsdata *lincsd) |
| 116 | { |
| 117 | return lincsd->rmsd_data; |
| 118 | } |
| 119 | |
| 120 | real lincs_rmsd(struct gmx_lincsdata *lincsd, gmx_bool bSD2) |
| 121 | { |
| 122 | if (lincsd->rmsd_data[0] > 0) |
| 123 | { |
| 124 | return sqrt(lincsd->rmsd_data[bSD2 ? 2 : 1]/lincsd->rmsd_data[0]); |
| 125 | } |
| 126 | else |
| 127 | { |
| 128 | return 0; |
| 129 | } |
| 130 | } |
| 131 | |
| 132 | /* Do a set of nrec LINCS matrix multiplications. |
| 133 | * This function will return with up to date thread-local |
| 134 | * constraint data, without an OpenMP barrier. |
| 135 | */ |
| 136 | static void lincs_matrix_expand(const struct gmx_lincsdata *lincsd, |
| 137 | int b0, int b1, |
| 138 | const real *blcc, |
| 139 | real *rhs1, real *rhs2, real *sol) |
| 140 | { |
| 141 | int nrec, rec, b, j, n, nr0, nr1; |
| 142 | real mvb, *swap; |
| 143 | int ntriangle, tb, bits; |
| 144 | const int *blnr = lincsd->blnr, *blbnb = lincsd->blbnb; |
| 145 | const int *triangle = lincsd->triangle, *tri_bits = lincsd->tri_bits; |
| 146 | |
| 147 | ntriangle = lincsd->ntriangle; |
| 148 | nrec = lincsd->nOrder; |
| 149 | |
| 150 | for (rec = 0; rec < nrec; rec++) |
| 151 | { |
| 152 | #pragma omp barrier |
| 153 | for (b = b0; b < b1; b++) |
| 154 | { |
| 155 | mvb = 0; |
| 156 | for (n = blnr[b]; n < blnr[b+1]; n++) |
| 157 | { |
| 158 | j = blbnb[n]; |
| 159 | mvb = mvb + blcc[n]*rhs1[j]; |
| 160 | } |
| 161 | rhs2[b] = mvb; |
| 162 | sol[b] = sol[b] + mvb; |
| 163 | } |
| 164 | swap = rhs1; |
| 165 | rhs1 = rhs2; |
| 166 | rhs2 = swap; |
| 167 | } /* nrec*(ncons+2*nrtot) flops */ |
| 168 | |
| 169 | if (ntriangle > 0) |
| 170 | { |
| 171 | /* Perform an extra nrec recursions for only the constraints |
| 172 | * involved in rigid triangles. |
| 173 | * In this way their accuracy should come close to those of the other |
| 174 | * constraints, since traingles of constraints can produce eigenvalues |
| 175 | * around 0.7, while the effective eigenvalue for bond constraints |
| 176 | * is around 0.4 (and 0.7*0.7=0.5). |
| 177 | */ |
| 178 | /* We need to copy the temporary array, since only the elements |
| 179 | * for constraints involved in triangles are updated and then |
| 180 | * the pointers are swapped. This saving copying the whole arrary. |
| 181 | * We need barrier as other threads might still be reading from rhs2. |
| 182 | */ |
| 183 | #pragma omp barrier |
| 184 | for (b = b0; b < b1; b++) |
| 185 | { |
| 186 | rhs2[b] = rhs1[b]; |
| 187 | } |
| 188 | #pragma omp barrier |
| 189 | #pragma omp master |
| 190 | { |
| 191 | for (rec = 0; rec < nrec; rec++) |
| 192 | { |
| 193 | for (tb = 0; tb < ntriangle; tb++) |
| 194 | { |
| 195 | b = triangle[tb]; |
| 196 | bits = tri_bits[tb]; |
| 197 | mvb = 0; |
| 198 | nr0 = blnr[b]; |
| 199 | nr1 = blnr[b+1]; |
| 200 | for (n = nr0; n < nr1; n++) |
| 201 | { |
| 202 | if (bits & (1<<(n-nr0))) |
| 203 | { |
| 204 | j = blbnb[n]; |
| 205 | mvb = mvb + blcc[n]*rhs1[j]; |
| 206 | } |
| 207 | } |
| 208 | rhs2[b] = mvb; |
| 209 | sol[b] = sol[b] + mvb; |
| 210 | } |
| 211 | swap = rhs1; |
| 212 | rhs1 = rhs2; |
| 213 | rhs2 = swap; |
| 214 | } |
| 215 | } /* flops count is missing here */ |
| 216 | |
| 217 | /* We need a barrier here as the calling routine will continue |
| 218 | * to operate on the thread-local constraints without barrier. |
| 219 | */ |
| 220 | #pragma omp barrier |
| 221 | } |
| 222 | } |
| 223 | |
| 224 | static void lincs_update_atoms_noind(int ncons, const int *bla, |
| 225 | real prefac, |
| 226 | const real *fac, rvec *r, |
| 227 | const real *invmass, |
| 228 | rvec *x) |
| 229 | { |
| 230 | int b, i, j; |
| 231 | real mvb, im1, im2, tmp0, tmp1, tmp2; |
| 232 | |
| 233 | if (invmass != NULL((void*)0)) |
| 234 | { |
| 235 | for (b = 0; b < ncons; b++) |
| 236 | { |
| 237 | i = bla[2*b]; |
| 238 | j = bla[2*b+1]; |
| 239 | mvb = prefac*fac[b]; |
| 240 | im1 = invmass[i]; |
| 241 | im2 = invmass[j]; |
| 242 | tmp0 = r[b][0]*mvb; |
| 243 | tmp1 = r[b][1]*mvb; |
| 244 | tmp2 = r[b][2]*mvb; |
| 245 | x[i][0] -= tmp0*im1; |
| 246 | x[i][1] -= tmp1*im1; |
| 247 | x[i][2] -= tmp2*im1; |
| 248 | x[j][0] += tmp0*im2; |
| 249 | x[j][1] += tmp1*im2; |
| 250 | x[j][2] += tmp2*im2; |
| 251 | } /* 16 ncons flops */ |
| 252 | } |
| 253 | else |
| 254 | { |
| 255 | for (b = 0; b < ncons; b++) |
| 256 | { |
| 257 | i = bla[2*b]; |
| 258 | j = bla[2*b+1]; |
| 259 | mvb = prefac*fac[b]; |
| 260 | tmp0 = r[b][0]*mvb; |
| 261 | tmp1 = r[b][1]*mvb; |
| 262 | tmp2 = r[b][2]*mvb; |
| 263 | x[i][0] -= tmp0; |
| 264 | x[i][1] -= tmp1; |
| 265 | x[i][2] -= tmp2; |
| 266 | x[j][0] += tmp0; |
| 267 | x[j][1] += tmp1; |
| 268 | x[j][2] += tmp2; |
| 269 | } |
| 270 | } |
| 271 | } |
| 272 | |
| 273 | static void lincs_update_atoms_ind(int ncons, const int *ind, const int *bla, |
| 274 | real prefac, |
| 275 | const real *fac, rvec *r, |
| 276 | const real *invmass, |
| 277 | rvec *x) |
| 278 | { |
| 279 | int bi, b, i, j; |
| 280 | real mvb, im1, im2, tmp0, tmp1, tmp2; |
| 281 | |
| 282 | if (invmass != NULL((void*)0)) |
| 283 | { |
| 284 | for (bi = 0; bi < ncons; bi++) |
| 285 | { |
| 286 | b = ind[bi]; |
| 287 | i = bla[2*b]; |
| 288 | j = bla[2*b+1]; |
| 289 | mvb = prefac*fac[b]; |
| 290 | im1 = invmass[i]; |
| 291 | im2 = invmass[j]; |
| 292 | tmp0 = r[b][0]*mvb; |
| 293 | tmp1 = r[b][1]*mvb; |
| 294 | tmp2 = r[b][2]*mvb; |
| 295 | x[i][0] -= tmp0*im1; |
| 296 | x[i][1] -= tmp1*im1; |
| 297 | x[i][2] -= tmp2*im1; |
| 298 | x[j][0] += tmp0*im2; |
| 299 | x[j][1] += tmp1*im2; |
| 300 | x[j][2] += tmp2*im2; |
| 301 | } /* 16 ncons flops */ |
| 302 | } |
| 303 | else |
| 304 | { |
| 305 | for (bi = 0; bi < ncons; bi++) |
| 306 | { |
| 307 | b = ind[bi]; |
| 308 | i = bla[2*b]; |
| 309 | j = bla[2*b+1]; |
| 310 | mvb = prefac*fac[b]; |
| 311 | tmp0 = r[b][0]*mvb; |
| 312 | tmp1 = r[b][1]*mvb; |
| 313 | tmp2 = r[b][2]*mvb; |
| 314 | x[i][0] -= tmp0; |
| 315 | x[i][1] -= tmp1; |
| 316 | x[i][2] -= tmp2; |
| 317 | x[j][0] += tmp0; |
| 318 | x[j][1] += tmp1; |
| 319 | x[j][2] += tmp2; |
| 320 | } /* 16 ncons flops */ |
| 321 | } |
| 322 | } |
| 323 | |
| 324 | static void lincs_update_atoms(struct gmx_lincsdata *li, int th, |
| 325 | real prefac, |
| 326 | const real *fac, rvec *r, |
| 327 | const real *invmass, |
| 328 | rvec *x) |
| 329 | { |
| 330 | if (li->nth == 1) |
| 331 | { |
| 332 | /* Single thread, we simply update for all constraints */ |
| 333 | lincs_update_atoms_noind(li->nc, li->bla, prefac, fac, r, invmass, x); |
| 334 | } |
| 335 | else |
| 336 | { |
| 337 | /* Update the atom vector components for our thread local |
| 338 | * constraints that only access our local atom range. |
| 339 | * This can be done without a barrier. |
| 340 | */ |
| 341 | lincs_update_atoms_ind(li->th[th].nind, li->th[th].ind, |
| 342 | li->bla, prefac, fac, r, invmass, x); |
| 343 | |
| 344 | if (li->th[li->nth].nind > 0) |
| 345 | { |
| 346 | /* Update the constraints that operate on atoms |
| 347 | * in multiple thread atom blocks on the master thread. |
| 348 | */ |
| 349 | #pragma omp barrier |
| 350 | #pragma omp master |
| 351 | { |
| 352 | lincs_update_atoms_ind(li->th[li->nth].nind, |
| 353 | li->th[li->nth].ind, |
| 354 | li->bla, prefac, fac, r, invmass, x); |
| 355 | } |
| 356 | } |
| 357 | } |
| 358 | } |
| 359 | |
| 360 | /* LINCS projection, works on derivatives of the coordinates */ |
| 361 | static void do_lincsp(rvec *x, rvec *f, rvec *fp, t_pbc *pbc, |
| 362 | struct gmx_lincsdata *lincsd, int th, |
| 363 | real *invmass, |
| 364 | int econq, real *dvdlambda, |
| 365 | gmx_bool bCalcVir, tensor rmdf) |
| 366 | { |
| 367 | int b0, b1, b, i, j, k, n; |
| 368 | real tmp0, tmp1, tmp2, im1, im2, mvb, rlen, len, wfac, lam; |
| 369 | rvec dx; |
| 370 | int *bla, *blnr, *blbnb; |
| 371 | rvec *r; |
| 372 | real *blc, *blmf, *blcc, *rhs1, *rhs2, *sol; |
| 373 | |
| 374 | b0 = lincsd->th[th].b0; |
| 375 | b1 = lincsd->th[th].b1; |
| 376 | |
| 377 | bla = lincsd->bla; |
| 378 | r = lincsd->tmpv; |
| 379 | blnr = lincsd->blnr; |
| 380 | blbnb = lincsd->blbnb; |
| 381 | if (econq != econqForce) |
| 382 | { |
| 383 | /* Use mass-weighted parameters */ |
| 384 | blc = lincsd->blc; |
| 385 | blmf = lincsd->blmf; |
| 386 | } |
| 387 | else |
| 388 | { |
| 389 | /* Use non mass-weighted parameters */ |
| 390 | blc = lincsd->blc1; |
| 391 | blmf = lincsd->blmf1; |
| 392 | } |
| 393 | blcc = lincsd->tmpncc; |
| 394 | rhs1 = lincsd->tmp1; |
| 395 | rhs2 = lincsd->tmp2; |
| 396 | sol = lincsd->tmp3; |
| 397 | |
| 398 | /* Compute normalized i-j vectors */ |
| 399 | if (pbc) |
| 400 | { |
| 401 | for (b = b0; b < b1; b++) |
| 402 | { |
| 403 | pbc_dx_aiuc(pbc, x[bla[2*b]], x[bla[2*b+1]], dx); |
| 404 | unitv(dx, r[b]); |
| 405 | } |
| 406 | } |
| 407 | else |
| 408 | { |
| 409 | for (b = b0; b < b1; b++) |
| 410 | { |
| 411 | rvec_sub(x[bla[2*b]], x[bla[2*b+1]], dx); |
| 412 | unitv(dx, r[b]); |
| 413 | } /* 16 ncons flops */ |
| 414 | } |
| 415 | |
| 416 | #pragma omp barrier |
| 417 | for (b = b0; b < b1; b++) |
| 418 | { |
| 419 | tmp0 = r[b][0]; |
| 420 | tmp1 = r[b][1]; |
| 421 | tmp2 = r[b][2]; |
| 422 | i = bla[2*b]; |
| 423 | j = bla[2*b+1]; |
| 424 | for (n = blnr[b]; n < blnr[b+1]; n++) |
| 425 | { |
| 426 | k = blbnb[n]; |
| 427 | blcc[n] = blmf[n]*(tmp0*r[k][0] + tmp1*r[k][1] + tmp2*r[k][2]); |
| 428 | } /* 6 nr flops */ |
| 429 | mvb = blc[b]*(tmp0*(f[i][0] - f[j][0]) + |
| 430 | tmp1*(f[i][1] - f[j][1]) + |
| 431 | tmp2*(f[i][2] - f[j][2])); |
| 432 | rhs1[b] = mvb; |
| 433 | sol[b] = mvb; |
| 434 | /* 7 flops */ |
| 435 | } |
| 436 | /* Together: 23*ncons + 6*nrtot flops */ |
| 437 | |
| 438 | lincs_matrix_expand(lincsd, b0, b1, blcc, rhs1, rhs2, sol); |
| 439 | /* nrec*(ncons+2*nrtot) flops */ |
| 440 | |
| 441 | if (econq == econqDeriv_FlexCon) |
| 442 | { |
| 443 | /* We only want to constraint the flexible constraints, |
| 444 | * so we mask out the normal ones by setting sol to 0. |
| 445 | */ |
| 446 | for (b = b0; b < b1; b++) |
| 447 | { |
| 448 | if (!(lincsd->bllen0[b] == 0 && lincsd->ddist[b] == 0)) |
| 449 | { |
| 450 | sol[b] = 0; |
| 451 | } |
| 452 | } |
| 453 | } |
| 454 | |
| 455 | /* We multiply sol by blc, so we can use lincs_update_atoms for OpenMP */ |
| 456 | for (b = b0; b < b1; b++) |
| 457 | { |
| 458 | sol[b] *= blc[b]; |
| 459 | } |
| 460 | |
| 461 | /* When constraining forces, we should not use mass weighting, |
| 462 | * so we pass invmass=NULL, which results in the use of 1 for all atoms. |
| 463 | */ |
| 464 | lincs_update_atoms(lincsd, th, 1.0, sol, r, |
| 465 | (econq != econqForce) ? invmass : NULL((void*)0), fp); |
| 466 | |
| 467 | if (dvdlambda != NULL((void*)0)) |
| 468 | { |
| 469 | #pragma omp barrier |
| 470 | for (b = b0; b < b1; b++) |
| 471 | { |
| 472 | *dvdlambda -= sol[b]*lincsd->ddist[b]; |
| 473 | } |
| 474 | /* 10 ncons flops */ |
| 475 | } |
| 476 | |
| 477 | if (bCalcVir) |
| 478 | { |
| 479 | /* Constraint virial, |
| 480 | * determines sum r_bond x delta f, |
| 481 | * where delta f is the constraint correction |
| 482 | * of the quantity that is being constrained. |
| 483 | */ |
| 484 | for (b = b0; b < b1; b++) |
| 485 | { |
| 486 | mvb = lincsd->bllen[b]*sol[b]; |
| 487 | for (i = 0; i < DIM3; i++) |
| 488 | { |
| 489 | tmp1 = mvb*r[b][i]; |
| 490 | for (j = 0; j < DIM3; j++) |
| 491 | { |
| 492 | rmdf[i][j] += tmp1*r[b][j]; |
| 493 | } |
| 494 | } |
| 495 | } /* 23 ncons flops */ |
| 496 | } |
| 497 | } |
| 498 | |
| 499 | static void do_lincs(rvec *x, rvec *xp, matrix box, t_pbc *pbc, |
| 500 | struct gmx_lincsdata *lincsd, int th, |
| 501 | real *invmass, |
| 502 | t_commrec *cr, |
| 503 | gmx_bool bCalcLambda, |
| 504 | real wangle, int *warn, |
| 505 | real invdt, rvec *v, |
| 506 | gmx_bool bCalcVir, tensor vir_r_m_dr) |
| 507 | { |
| 508 | int b0, b1, b, i, j, k, n, iter; |
| 509 | real tmp0, tmp1, tmp2, im1, im2, mvb, rlen, len, len2, dlen2, wfac; |
| 510 | rvec dx; |
| 511 | int *bla, *blnr, *blbnb; |
| 512 | rvec *r; |
| 513 | real *blc, *blmf, *bllen, *blcc, *rhs1, *rhs2, *sol, *blc_sol, *mlambda; |
| 514 | int *nlocat; |
| 515 | |
| 516 | b0 = lincsd->th[th].b0; |
| 517 | b1 = lincsd->th[th].b1; |
| 518 | |
| 519 | bla = lincsd->bla; |
| 520 | r = lincsd->tmpv; |
| 521 | blnr = lincsd->blnr; |
| 522 | blbnb = lincsd->blbnb; |
| 523 | blc = lincsd->blc; |
| 524 | blmf = lincsd->blmf; |
| 525 | bllen = lincsd->bllen; |
| 526 | blcc = lincsd->tmpncc; |
| 527 | rhs1 = lincsd->tmp1; |
| 528 | rhs2 = lincsd->tmp2; |
| 529 | sol = lincsd->tmp3; |
| 530 | blc_sol = lincsd->tmp4; |
| 531 | mlambda = lincsd->mlambda; |
| 532 | |
| 533 | if (DOMAINDECOMP(cr)(((cr)->dd != ((void*)0)) && ((cr)->nnodes > 1)) && cr->dd->constraints) |
| 534 | { |
| 535 | nlocat = dd_constraints_nlocalatoms(cr->dd); |
| 536 | } |
| 537 | else |
| 538 | { |
| 539 | nlocat = NULL((void*)0); |
| 540 | } |
| 541 | |
| 542 | if (pbc) |
| 543 | { |
| 544 | /* Compute normalized i-j vectors */ |
| 545 | for (b = b0; b < b1; b++) |
| 546 | { |
| 547 | pbc_dx_aiuc(pbc, x[bla[2*b]], x[bla[2*b+1]], dx); |
| 548 | unitv(dx, r[b]); |
| 549 | } |
| 550 | #pragma omp barrier |
| 551 | for (b = b0; b < b1; b++) |
| 552 | { |
| 553 | for (n = blnr[b]; n < blnr[b+1]; n++) |
| 554 | { |
| 555 | blcc[n] = blmf[n]*iprod(r[b], r[blbnb[n]]); |
| 556 | } |
| 557 | pbc_dx_aiuc(pbc, xp[bla[2*b]], xp[bla[2*b+1]], dx); |
| 558 | mvb = blc[b]*(iprod(r[b], dx) - bllen[b]); |
| 559 | rhs1[b] = mvb; |
| 560 | sol[b] = mvb; |
| 561 | } |
| 562 | } |
| 563 | else |
| 564 | { |
| 565 | /* Compute normalized i-j vectors */ |
| 566 | for (b = b0; b < b1; b++) |
| 567 | { |
| 568 | i = bla[2*b]; |
| 569 | j = bla[2*b+1]; |
| 570 | tmp0 = x[i][0] - x[j][0]; |
| 571 | tmp1 = x[i][1] - x[j][1]; |
| 572 | tmp2 = x[i][2] - x[j][2]; |
| 573 | rlen = gmx_invsqrt(tmp0*tmp0+tmp1*tmp1+tmp2*tmp2)gmx_software_invsqrt(tmp0*tmp0+tmp1*tmp1+tmp2*tmp2); |
| 574 | r[b][0] = rlen*tmp0; |
| 575 | r[b][1] = rlen*tmp1; |
| 576 | r[b][2] = rlen*tmp2; |
| 577 | } /* 16 ncons flops */ |
| 578 | |
| 579 | #pragma omp barrier |
| 580 | for (b = b0; b < b1; b++) |
| 581 | { |
| 582 | tmp0 = r[b][0]; |
| 583 | tmp1 = r[b][1]; |
| 584 | tmp2 = r[b][2]; |
| 585 | len = bllen[b]; |
| 586 | i = bla[2*b]; |
| 587 | j = bla[2*b+1]; |
| 588 | for (n = blnr[b]; n < blnr[b+1]; n++) |
| 589 | { |
| 590 | k = blbnb[n]; |
| 591 | blcc[n] = blmf[n]*(tmp0*r[k][0] + tmp1*r[k][1] + tmp2*r[k][2]); |
| 592 | } /* 6 nr flops */ |
| 593 | mvb = blc[b]*(tmp0*(xp[i][0] - xp[j][0]) + |
| 594 | tmp1*(xp[i][1] - xp[j][1]) + |
| 595 | tmp2*(xp[i][2] - xp[j][2]) - len); |
| 596 | rhs1[b] = mvb; |
| 597 | sol[b] = mvb; |
| 598 | /* 10 flops */ |
| 599 | } |
| 600 | /* Together: 26*ncons + 6*nrtot flops */ |
| 601 | } |
| 602 | |
| 603 | lincs_matrix_expand(lincsd, b0, b1, blcc, rhs1, rhs2, sol); |
| 604 | /* nrec*(ncons+2*nrtot) flops */ |
| 605 | |
| 606 | for (b = b0; b < b1; b++) |
| 607 | { |
| 608 | mlambda[b] = blc[b]*sol[b]; |
| 609 | } |
| 610 | |
| 611 | /* Update the coordinates */ |
| 612 | lincs_update_atoms(lincsd, th, 1.0, mlambda, r, invmass, xp); |
| 613 | |
| 614 | /* |
| 615 | ******** Correction for centripetal effects ******** |
| 616 | */ |
| 617 | |
| 618 | wfac = cos(DEG2RAD(3.14159265358979323846/180.0)*wangle); |
| 619 | wfac = wfac*wfac; |
| 620 | |
| 621 | for (iter = 0; iter < lincsd->nIter; iter++) |
| 622 | { |
| 623 | if ((lincsd->bCommIter && DOMAINDECOMP(cr)(((cr)->dd != ((void*)0)) && ((cr)->nnodes > 1)) && cr->dd->constraints)) |
| 624 | { |
| 625 | #pragma omp barrier |
| 626 | #pragma omp master |
| 627 | { |
| 628 | /* Communicate the corrected non-local coordinates */ |
| 629 | if (DOMAINDECOMP(cr)(((cr)->dd != ((void*)0)) && ((cr)->nnodes > 1))) |
| 630 | { |
| 631 | dd_move_x_constraints(cr->dd, box, xp, NULL((void*)0), FALSE0); |
| 632 | } |
| 633 | } |
| 634 | } |
| 635 | |
| 636 | #pragma omp barrier |
| 637 | for (b = b0; b < b1; b++) |
| 638 | { |
| 639 | len = bllen[b]; |
| 640 | if (pbc) |
| 641 | { |
| 642 | pbc_dx_aiuc(pbc, xp[bla[2*b]], xp[bla[2*b+1]], dx); |
| 643 | } |
| 644 | else |
| 645 | { |
| 646 | rvec_sub(xp[bla[2*b]], xp[bla[2*b+1]], dx); |
| 647 | } |
| 648 | len2 = len*len; |
| 649 | dlen2 = 2*len2 - norm2(dx); |
| 650 | if (dlen2 < wfac*len2 && (nlocat == NULL((void*)0) || nlocat[b])) |
| 651 | { |
| 652 | *warn = b; |
| 653 | } |
| 654 | if (dlen2 > 0) |
| 655 | { |
| 656 | mvb = blc[b]*(len - dlen2*gmx_invsqrt(dlen2)gmx_software_invsqrt(dlen2)); |
| 657 | } |
| 658 | else |
| 659 | { |
| 660 | mvb = blc[b]*len; |
| 661 | } |
| 662 | rhs1[b] = mvb; |
| 663 | sol[b] = mvb; |
| 664 | } /* 20*ncons flops */ |
| 665 | |
| 666 | lincs_matrix_expand(lincsd, b0, b1, blcc, rhs1, rhs2, sol); |
| 667 | /* nrec*(ncons+2*nrtot) flops */ |
| 668 | |
| 669 | for (b = b0; b < b1; b++) |
| 670 | { |
| 671 | mvb = blc[b]*sol[b]; |
| 672 | blc_sol[b] = mvb; |
| 673 | mlambda[b] += mvb; |
| 674 | } |
| 675 | |
| 676 | /* Update the coordinates */ |
| 677 | lincs_update_atoms(lincsd, th, 1.0, blc_sol, r, invmass, xp); |
| 678 | } |
| 679 | /* nit*ncons*(37+9*nrec) flops */ |
| 680 | |
| 681 | if (v != NULL((void*)0)) |
| 682 | { |
| 683 | /* Update the velocities */ |
| 684 | lincs_update_atoms(lincsd, th, invdt, mlambda, r, invmass, v); |
| 685 | /* 16 ncons flops */ |
| 686 | } |
| 687 | |
| 688 | if (nlocat != NULL((void*)0) && bCalcLambda) |
| 689 | { |
| 690 | /* In lincs_update_atoms thread might cross-read mlambda */ |
| 691 | #pragma omp barrier |
| 692 | |
| 693 | /* Only account for local atoms */ |
| 694 | for (b = b0; b < b1; b++) |
| 695 | { |
| 696 | mlambda[b] *= 0.5*nlocat[b]; |
| 697 | } |
| 698 | } |
| 699 | |
| 700 | if (bCalcVir) |
| 701 | { |
| 702 | /* Constraint virial */ |
| 703 | for (b = b0; b < b1; b++) |
| 704 | { |
| 705 | tmp0 = -bllen[b]*mlambda[b]; |
| 706 | for (i = 0; i < DIM3; i++) |
| 707 | { |
| 708 | tmp1 = tmp0*r[b][i]; |
| 709 | for (j = 0; j < DIM3; j++) |
| 710 | { |
| 711 | vir_r_m_dr[i][j] -= tmp1*r[b][j]; |
| 712 | } |
| 713 | } |
| 714 | } /* 22 ncons flops */ |
| 715 | } |
| 716 | |
| 717 | /* Total: |
| 718 | * 26*ncons + 6*nrtot + nrec*(ncons+2*nrtot) |
| 719 | * + nit * (20*ncons + nrec*(ncons+2*nrtot) + 17 ncons) |
| 720 | * |
| 721 | * (26+nrec)*ncons + (6+2*nrec)*nrtot |
| 722 | * + nit * ((37+nrec)*ncons + 2*nrec*nrtot) |
| 723 | * if nit=1 |
| 724 | * (63+nrec)*ncons + (6+4*nrec)*nrtot |
| 725 | */ |
| 726 | } |
| 727 | |
| 728 | void set_lincs_matrix(struct gmx_lincsdata *li, real *invmass, real lambda) |
| 729 | { |
| 730 | int i, a1, a2, n, k, sign, center; |
| 731 | int end, nk, kk; |
| 732 | const real invsqrt2 = 0.7071067811865475244; |
| 733 | |
| 734 | for (i = 0; (i < li->nc); i++) |
| 735 | { |
| 736 | a1 = li->bla[2*i]; |
| 737 | a2 = li->bla[2*i+1]; |
| 738 | li->blc[i] = gmx_invsqrt(invmass[a1] + invmass[a2])gmx_software_invsqrt(invmass[a1] + invmass[a2]); |
| 739 | li->blc1[i] = invsqrt2; |
| 740 | } |
| 741 | |
| 742 | /* Construct the coupling coefficient matrix blmf */ |
| 743 | li->ntriangle = 0; |
| 744 | li->ncc_triangle = 0; |
| 745 | for (i = 0; (i < li->nc); i++) |
| 746 | { |
| 747 | a1 = li->bla[2*i]; |
| 748 | a2 = li->bla[2*i+1]; |
| 749 | for (n = li->blnr[i]; (n < li->blnr[i+1]); n++) |
| 750 | { |
| 751 | k = li->blbnb[n]; |
| 752 | if (a1 == li->bla[2*k] || a2 == li->bla[2*k+1]) |
| 753 | { |
| 754 | sign = -1; |
| 755 | } |
| 756 | else |
| 757 | { |
| 758 | sign = 1; |
| 759 | } |
| 760 | if (a1 == li->bla[2*k] || a1 == li->bla[2*k+1]) |
| 761 | { |
| 762 | center = a1; |
| 763 | end = a2; |
| 764 | } |
| 765 | else |
| 766 | { |
| 767 | center = a2; |
| 768 | end = a1; |
| 769 | } |
| 770 | li->blmf[n] = sign*invmass[center]*li->blc[i]*li->blc[k]; |
| 771 | li->blmf1[n] = sign*0.5; |
| 772 | if (li->ncg_triangle > 0) |
| 773 | { |
| 774 | /* Look for constraint triangles */ |
| 775 | for (nk = li->blnr[k]; (nk < li->blnr[k+1]); nk++) |
| 776 | { |
| 777 | kk = li->blbnb[nk]; |
| 778 | if (kk != i && kk != k && |
| 779 | (li->bla[2*kk] == end || li->bla[2*kk+1] == end)) |
| 780 | { |
| 781 | if (li->ntriangle == 0 || |
| 782 | li->triangle[li->ntriangle-1] < i) |
| 783 | { |
| 784 | /* Add this constraint to the triangle list */ |
| 785 | li->triangle[li->ntriangle] = i; |
| 786 | li->tri_bits[li->ntriangle] = 0; |
| 787 | li->ntriangle++; |
| 788 | if (li->blnr[i+1] - li->blnr[i] > sizeof(li->tri_bits[0])*8 - 1) |
| 789 | { |
| 790 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/clincs.c", 790, "A constraint is connected to %d constraints, this is more than the %d allowed for constraints participating in triangles", |
| 791 | li->blnr[i+1] - li->blnr[i], |
| 792 | sizeof(li->tri_bits[0])*8-1); |
| 793 | } |
| 794 | } |
| 795 | li->tri_bits[li->ntriangle-1] |= (1<<(n-li->blnr[i])); |
| 796 | li->ncc_triangle++; |
| 797 | } |
| 798 | } |
| 799 | } |
| 800 | } |
| 801 | } |
| 802 | |
| 803 | if (debug) |
| 804 | { |
| 805 | fprintf(debug, "Of the %d constraints %d participate in triangles\n", |
| 806 | li->nc, li->ntriangle); |
| 807 | fprintf(debug, "There are %d couplings of which %d in triangles\n", |
| 808 | li->ncc, li->ncc_triangle); |
| 809 | } |
| 810 | |
| 811 | /* Set matlam, |
| 812 | * so we know with which lambda value the masses have been set. |
| 813 | */ |
| 814 | li->matlam = lambda; |
| 815 | } |
| 816 | |
| 817 | static int count_triangle_constraints(t_ilist *ilist, t_blocka *at2con) |
| 818 | { |
| 819 | int ncon1, ncon_tot; |
| 820 | int c0, a00, a01, n1, c1, a10, a11, ac1, n2, c2, a20, a21; |
| 821 | int ncon_triangle; |
| 822 | gmx_bool bTriangle; |
| 823 | t_iatom *ia1, *ia2, *iap; |
| 824 | |
| 825 | ncon1 = ilist[F_CONSTR].nr/3; |
| 826 | ncon_tot = ncon1 + ilist[F_CONSTRNC].nr/3; |
| 827 | |
| 828 | ia1 = ilist[F_CONSTR].iatoms; |
| 829 | ia2 = ilist[F_CONSTRNC].iatoms; |
| 830 | |
| 831 | ncon_triangle = 0; |
| 832 | for (c0 = 0; c0 < ncon_tot; c0++) |
| 833 | { |
| 834 | bTriangle = FALSE0; |
| 835 | iap = constr_iatomptr(ncon1, ia1, ia2, c0)((c0) < (ncon1) ? (ia1)+(c0)*3 : (ia2)+(c0-ncon1)*3); |
| 836 | a00 = iap[1]; |
| 837 | a01 = iap[2]; |
| 838 | for (n1 = at2con->index[a01]; n1 < at2con->index[a01+1]; n1++) |
| 839 | { |
| 840 | c1 = at2con->a[n1]; |
| 841 | if (c1 != c0) |
| 842 | { |
| 843 | iap = constr_iatomptr(ncon1, ia1, ia2, c1)((c1) < (ncon1) ? (ia1)+(c1)*3 : (ia2)+(c1-ncon1)*3); |
| 844 | a10 = iap[1]; |
| 845 | a11 = iap[2]; |
| 846 | if (a10 == a01) |
| 847 | { |
| 848 | ac1 = a11; |
| 849 | } |
| 850 | else |
| 851 | { |
| 852 | ac1 = a10; |
| 853 | } |
| 854 | for (n2 = at2con->index[ac1]; n2 < at2con->index[ac1+1]; n2++) |
| 855 | { |
| 856 | c2 = at2con->a[n2]; |
| 857 | if (c2 != c0 && c2 != c1) |
| 858 | { |
| 859 | iap = constr_iatomptr(ncon1, ia1, ia2, c2)((c2) < (ncon1) ? (ia1)+(c2)*3 : (ia2)+(c2-ncon1)*3); |
| 860 | a20 = iap[1]; |
| 861 | a21 = iap[2]; |
| 862 | if (a20 == a00 || a21 == a00) |
| 863 | { |
| 864 | bTriangle = TRUE1; |
| 865 | } |
| 866 | } |
| 867 | } |
| 868 | } |
| 869 | } |
| 870 | if (bTriangle) |
| 871 | { |
| 872 | ncon_triangle++; |
| 873 | } |
| 874 | } |
| 875 | |
| 876 | return ncon_triangle; |
| 877 | } |
| 878 | |
| 879 | static gmx_bool more_than_two_sequential_constraints(const t_ilist *ilist, |
| 880 | const t_blocka *at2con) |
| 881 | { |
| 882 | t_iatom *ia1, *ia2, *iap; |
| 883 | int ncon1, ncon_tot, c; |
| 884 | int a1, a2; |
| 885 | gmx_bool bMoreThanTwoSequentialConstraints; |
| 886 | |
| 887 | ncon1 = ilist[F_CONSTR].nr/3; |
| 888 | ncon_tot = ncon1 + ilist[F_CONSTRNC].nr/3; |
| 889 | |
| 890 | ia1 = ilist[F_CONSTR].iatoms; |
| 891 | ia2 = ilist[F_CONSTRNC].iatoms; |
| 892 | |
| 893 | bMoreThanTwoSequentialConstraints = FALSE0; |
| 894 | for (c = 0; c < ncon_tot && !bMoreThanTwoSequentialConstraints; c++) |
| 895 | { |
| 896 | iap = constr_iatomptr(ncon1, ia1, ia2, c)((c) < (ncon1) ? (ia1)+(c)*3 : (ia2)+(c-ncon1)*3); |
| 897 | a1 = iap[1]; |
| 898 | a2 = iap[2]; |
| 899 | /* Check if this constraint has constraints connected at both atoms */ |
| 900 | if (at2con->index[a1+1] - at2con->index[a1] > 1 && |
| 901 | at2con->index[a2+1] - at2con->index[a2] > 1) |
| 902 | { |
| 903 | bMoreThanTwoSequentialConstraints = TRUE1; |
| 904 | } |
| 905 | } |
| 906 | |
| 907 | return bMoreThanTwoSequentialConstraints; |
| 908 | } |
| 909 | |
| 910 | static int int_comp(const void *a, const void *b) |
| 911 | { |
| 912 | return (*(int *)a) - (*(int *)b); |
| 913 | } |
| 914 | |
| 915 | gmx_lincsdata_t init_lincs(FILE *fplog, gmx_mtop_t *mtop, |
| 916 | int nflexcon_global, t_blocka *at2con, |
| 917 | gmx_bool bPLINCS, int nIter, int nProjOrder) |
| 918 | { |
| 919 | struct gmx_lincsdata *li; |
| 920 | int mb; |
| 921 | gmx_moltype_t *molt; |
| 922 | |
| 923 | if (fplog) |
| 924 | { |
| 925 | fprintf(fplog, "\nInitializing%s LINear Constraint Solver\n", |
| 926 | bPLINCS ? " Parallel" : ""); |
| 927 | } |
| 928 | |
| 929 | snew(li, 1)(li) = save_calloc("li", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/clincs.c" , 929, (1), sizeof(*(li))); |
| 930 | |
| 931 | li->ncg = |
| 932 | gmx_mtop_ftype_count(mtop, F_CONSTR) + |
| 933 | gmx_mtop_ftype_count(mtop, F_CONSTRNC); |
| 934 | li->ncg_flex = nflexcon_global; |
| 935 | |
| 936 | li->nIter = nIter; |
| 937 | li->nOrder = nProjOrder; |
| 938 | |
| 939 | li->ncg_triangle = 0; |
| 940 | li->bCommIter = FALSE0; |
| 941 | for (mb = 0; mb < mtop->nmolblock; mb++) |
| 942 | { |
| 943 | molt = &mtop->moltype[mtop->molblock[mb].type]; |
| 944 | li->ncg_triangle += |
| 945 | mtop->molblock[mb].nmol* |
| 946 | count_triangle_constraints(molt->ilist, |
| 947 | &at2con[mtop->molblock[mb].type]); |
| 948 | if (bPLINCS && li->bCommIter == FALSE0) |
| 949 | { |
| 950 | /* Check if we need to communicate not only before LINCS, |
| 951 | * but also before each iteration. |
| 952 | * The check for only two sequential constraints is only |
| 953 | * useful for the common case of H-bond only constraints. |
| 954 | * With more effort we could also make it useful for small |
| 955 | * molecules with nr. sequential constraints <= nOrder-1. |
| 956 | */ |
| 957 | li->bCommIter = (li->nOrder < 1 || more_than_two_sequential_constraints(molt->ilist, &at2con[mtop->molblock[mb].type])); |
| 958 | } |
| 959 | } |
| 960 | if (debug && bPLINCS) |
| 961 | { |
| 962 | fprintf(debug, "PLINCS communication before each iteration: %d\n", |
| 963 | li->bCommIter); |
| 964 | } |
| 965 | |
| 966 | /* LINCS can run on any number of threads. |
| 967 | * Currently the number is fixed for the whole simulation, |
| 968 | * but it could be set in set_lincs(). |
| 969 | */ |
| 970 | li->nth = gmx_omp_nthreads_get(emntLINCS); |
| 971 | if (li->nth == 1) |
| 972 | { |
| 973 | snew(li->th, 1)(li->th) = save_calloc("li->th", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/clincs.c" , 973, (1), sizeof(*(li->th))); |
| 974 | } |
| 975 | else |
| 976 | { |
| 977 | /* Allocate an extra elements for "thread-overlap" constraints */ |
| 978 | snew(li->th, li->nth+1)(li->th) = save_calloc("li->th", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/clincs.c" , 978, (li->nth+1), sizeof(*(li->th))); |
| 979 | } |
| 980 | if (debug) |
| 981 | { |
| 982 | fprintf(debug, "LINCS: using %d threads\n", li->nth); |
| 983 | } |
| 984 | |
| 985 | if (bPLINCS || li->ncg_triangle > 0) |
| 986 | { |
| 987 | please_cite(fplog, "Hess2008a"); |
| 988 | } |
| 989 | else |
| 990 | { |
| 991 | please_cite(fplog, "Hess97a"); |
| 992 | } |
| 993 | |
| 994 | if (fplog) |
| 995 | { |
| 996 | fprintf(fplog, "The number of constraints is %d\n", li->ncg); |
| 997 | if (bPLINCS) |
| 998 | { |
| 999 | fprintf(fplog, "There are inter charge-group constraints,\n" |
| 1000 | "will communicate selected coordinates each lincs iteration\n"); |
| 1001 | } |
| 1002 | if (li->ncg_triangle > 0) |
| 1003 | { |
| 1004 | fprintf(fplog, |
| 1005 | "%d constraints are involved in constraint triangles,\n" |
| 1006 | "will apply an additional matrix expansion of order %d for couplings\n" |
| 1007 | "between constraints inside triangles\n", |
| 1008 | li->ncg_triangle, li->nOrder); |
| 1009 | } |
| 1010 | } |
| 1011 | |
| 1012 | return li; |
| 1013 | } |
| 1014 | |
| 1015 | /* Sets up the work division over the threads */ |
| 1016 | static void lincs_thread_setup(struct gmx_lincsdata *li, int natoms) |
| 1017 | { |
| 1018 | lincs_thread_t *li_m; |
| 1019 | int th; |
| 1020 | unsigned *atf; |
| 1021 | int a; |
| 1022 | |
| 1023 | if (natoms > li->atf_nalloc) |
| 1024 | { |
| 1025 | li->atf_nalloc = over_alloc_large(natoms)(int)(1.19*(natoms) + 1000); |
| 1026 | srenew(li->atf, li->atf_nalloc)(li->atf) = save_realloc("li->atf", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/clincs.c" , 1026, (li->atf), (li->atf_nalloc), sizeof(*(li->atf ))); |
| 1027 | } |
| 1028 | |
| 1029 | atf = li->atf; |
| 1030 | /* Clear the atom flags */ |
| 1031 | for (a = 0; a < natoms; a++) |
| 1032 | { |
| 1033 | atf[a] = 0; |
| 1034 | } |
| 1035 | |
| 1036 | for (th = 0; th < li->nth; th++) |
| 1037 | { |
| 1038 | lincs_thread_t *li_th; |
| 1039 | int b; |
| 1040 | |
| 1041 | li_th = &li->th[th]; |
| 1042 | |
| 1043 | /* The constraints are divided equally over the threads */ |
| 1044 | li_th->b0 = (li->nc* th )/li->nth; |
| 1045 | li_th->b1 = (li->nc*(th+1))/li->nth; |
| 1046 | |
| 1047 | if (th < sizeof(*atf)*8) |
| 1048 | { |
| 1049 | /* For each atom set a flag for constraints from each */ |
| 1050 | for (b = li_th->b0; b < li_th->b1; b++) |
| 1051 | { |
| 1052 | atf[li->bla[b*2] ] |= (1U<<th); |
| 1053 | atf[li->bla[b*2+1]] |= (1U<<th); |
| 1054 | } |
| 1055 | } |
| 1056 | } |
| 1057 | |
| 1058 | #pragma omp parallel for num_threads(li->nth) schedule(static) |
| 1059 | for (th = 0; th < li->nth; th++) |
| 1060 | { |
| 1061 | lincs_thread_t *li_th; |
| 1062 | unsigned mask; |
| 1063 | int b; |
| 1064 | |
| 1065 | li_th = &li->th[th]; |
| 1066 | |
| 1067 | if (li_th->b1 - li_th->b0 > li_th->ind_nalloc) |
| 1068 | { |
| 1069 | li_th->ind_nalloc = over_alloc_large(li_th->b1-li_th->b0)(int)(1.19*(li_th->b1-li_th->b0) + 1000); |
| 1070 | srenew(li_th->ind, li_th->ind_nalloc)(li_th->ind) = save_realloc("li_th->ind", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/clincs.c" , 1070, (li_th->ind), (li_th->ind_nalloc), sizeof(*(li_th ->ind))); |
| 1071 | srenew(li_th->ind_r, li_th->ind_nalloc)(li_th->ind_r) = save_realloc("li_th->ind_r", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/clincs.c" , 1071, (li_th->ind_r), (li_th->ind_nalloc), sizeof(*(li_th ->ind_r))); |
| 1072 | } |
| 1073 | |
| 1074 | if (th < sizeof(*atf)*8) |
| 1075 | { |
| 1076 | mask = (1U<<th) - 1U; |
| 1077 | |
| 1078 | li_th->nind = 0; |
| 1079 | li_th->nind_r = 0; |
| 1080 | for (b = li_th->b0; b < li_th->b1; b++) |
| 1081 | { |
| 1082 | /* We let the constraint with the lowest thread index |
| 1083 | * operate on atoms with constraints from multiple threads. |
| 1084 | */ |
| 1085 | if (((atf[li->bla[b*2]] & mask) == 0) && |
| 1086 | ((atf[li->bla[b*2+1]] & mask) == 0)) |
| 1087 | { |
| 1088 | /* Add the constraint to the local atom update index */ |
| 1089 | li_th->ind[li_th->nind++] = b; |
| 1090 | } |
| 1091 | else |
| 1092 | { |
| 1093 | /* Add the constraint to the rest block */ |
| 1094 | li_th->ind_r[li_th->nind_r++] = b; |
| 1095 | } |
| 1096 | } |
| 1097 | } |
| 1098 | else |
| 1099 | { |
| 1100 | /* We are out of bits, assign all constraints to rest */ |
| 1101 | for (b = li_th->b0; b < li_th->b1; b++) |
| 1102 | { |
| 1103 | li_th->ind_r[li_th->nind_r++] = b; |
| 1104 | } |
| 1105 | } |
| 1106 | } |
| 1107 | |
| 1108 | /* We need to copy all constraints which have not be assigned |
| 1109 | * to a thread to a separate list which will be handled by one thread. |
| 1110 | */ |
| 1111 | li_m = &li->th[li->nth]; |
| 1112 | |
| 1113 | li_m->nind = 0; |
| 1114 | for (th = 0; th < li->nth; th++) |
| 1115 | { |
| 1116 | lincs_thread_t *li_th; |
| 1117 | int b; |
| 1118 | |
| 1119 | li_th = &li->th[th]; |
| 1120 | |
| 1121 | if (li_m->nind + li_th->nind_r > li_m->ind_nalloc) |
| 1122 | { |
| 1123 | li_m->ind_nalloc = over_alloc_large(li_m->nind+li_th->nind_r)(int)(1.19*(li_m->nind+li_th->nind_r) + 1000); |
| 1124 | srenew(li_m->ind, li_m->ind_nalloc)(li_m->ind) = save_realloc("li_m->ind", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/clincs.c" , 1124, (li_m->ind), (li_m->ind_nalloc), sizeof(*(li_m-> ind))); |
| 1125 | } |
| 1126 | |
| 1127 | for (b = 0; b < li_th->nind_r; b++) |
| 1128 | { |
| 1129 | li_m->ind[li_m->nind++] = li_th->ind_r[b]; |
| 1130 | } |
| 1131 | |
| 1132 | if (debug) |
| 1133 | { |
| 1134 | fprintf(debug, "LINCS thread %d: %d constraints\n", |
| 1135 | th, li_th->nind); |
| 1136 | } |
| 1137 | } |
| 1138 | |
| 1139 | if (debug) |
| 1140 | { |
| 1141 | fprintf(debug, "LINCS thread r: %d constraints\n", |
| 1142 | li_m->nind); |
| 1143 | } |
| 1144 | } |
| 1145 | |
| 1146 | |
| 1147 | void set_lincs(t_idef *idef, t_mdatoms *md, |
| 1148 | gmx_bool bDynamics, t_commrec *cr, |
| 1149 | struct gmx_lincsdata *li) |
| 1150 | { |
| 1151 | int start, natoms, nflexcon; |
| 1152 | t_blocka at2con; |
| 1153 | t_iatom *iatom; |
| 1154 | int i, k, ncc_alloc, ni, con, nconnect, concon; |
| 1155 | int type, a1, a2; |
| 1156 | real lenA = 0, lenB; |
| 1157 | gmx_bool bLocal; |
| 1158 | |
| 1159 | li->nc = 0; |
| 1160 | li->ncc = 0; |
| 1161 | /* Zero the thread index ranges. |
| 1162 | * Otherwise without local constraints we could return with old ranges. |
| 1163 | */ |
| 1164 | for (i = 0; i < li->nth; i++) |
| 1165 | { |
| 1166 | li->th[i].b0 = 0; |
| 1167 | li->th[i].b1 = 0; |
| 1168 | li->th[i].nind = 0; |
| 1169 | } |
| 1170 | if (li->nth > 1) |
| 1171 | { |
| 1172 | li->th[li->nth].nind = 0; |
| 1173 | } |
| 1174 | |
| 1175 | /* This is the local topology, so there are only F_CONSTR constraints */ |
| 1176 | if (idef->il[F_CONSTR].nr == 0) |
| 1177 | { |
| 1178 | /* There are no constraints, |
| 1179 | * we do not need to fill any data structures. |
| 1180 | */ |
| 1181 | return; |
| 1182 | } |
| 1183 | |
| 1184 | if (debug) |
| 1185 | { |
| 1186 | fprintf(debug, "Building the LINCS connectivity\n"); |
| 1187 | } |
| 1188 | |
| 1189 | if (DOMAINDECOMP(cr)(((cr)->dd != ((void*)0)) && ((cr)->nnodes > 1))) |
| 1190 | { |
| 1191 | if (cr->dd->constraints) |
| 1192 | { |
| 1193 | dd_get_constraint_range(cr->dd, &start, &natoms); |
| 1194 | } |
| 1195 | else |
| 1196 | { |
| 1197 | natoms = cr->dd->nat_home; |
| 1198 | } |
| 1199 | start = 0; |
| 1200 | } |
| 1201 | else |
| 1202 | { |
| 1203 | start = 0; |
| 1204 | natoms = md->homenr; |
| 1205 | } |
| 1206 | at2con = make_at2con(start, natoms, idef->il, idef->iparams, bDynamics, |
| 1207 | &nflexcon); |
| 1208 | |
| 1209 | |
| 1210 | if (idef->il[F_CONSTR].nr/3 > li->nc_alloc || li->nc_alloc == 0) |
| 1211 | { |
| 1212 | li->nc_alloc = over_alloc_dd(idef->il[F_CONSTR].nr/3); |
| 1213 | srenew(li->bllen0, li->nc_alloc)(li->bllen0) = save_realloc("li->bllen0", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/clincs.c" , 1213, (li->bllen0), (li->nc_alloc), sizeof(*(li->bllen0 ))); |
| 1214 | srenew(li->ddist, li->nc_alloc)(li->ddist) = save_realloc("li->ddist", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/clincs.c" , 1214, (li->ddist), (li->nc_alloc), sizeof(*(li->ddist ))); |
| 1215 | srenew(li->bla, 2*li->nc_alloc)(li->bla) = save_realloc("li->bla", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/clincs.c" , 1215, (li->bla), (2*li->nc_alloc), sizeof(*(li->bla ))); |
| 1216 | srenew(li->blc, li->nc_alloc)(li->blc) = save_realloc("li->blc", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/clincs.c" , 1216, (li->blc), (li->nc_alloc), sizeof(*(li->blc) )); |
| 1217 | srenew(li->blc1, li->nc_alloc)(li->blc1) = save_realloc("li->blc1", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/clincs.c" , 1217, (li->blc1), (li->nc_alloc), sizeof(*(li->blc1 ))); |
| 1218 | srenew(li->blnr, li->nc_alloc+1)(li->blnr) = save_realloc("li->blnr", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/clincs.c" , 1218, (li->blnr), (li->nc_alloc+1), sizeof(*(li->blnr ))); |
| 1219 | srenew(li->bllen, li->nc_alloc)(li->bllen) = save_realloc("li->bllen", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/clincs.c" , 1219, (li->bllen), (li->nc_alloc), sizeof(*(li->bllen ))); |
| 1220 | srenew(li->tmpv, li->nc_alloc)(li->tmpv) = save_realloc("li->tmpv", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/clincs.c" , 1220, (li->tmpv), (li->nc_alloc), sizeof(*(li->tmpv ))); |
| 1221 | srenew(li->tmp1, li->nc_alloc)(li->tmp1) = save_realloc("li->tmp1", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/clincs.c" , 1221, (li->tmp1), (li->nc_alloc), sizeof(*(li->tmp1 ))); |
| 1222 | srenew(li->tmp2, li->nc_alloc)(li->tmp2) = save_realloc("li->tmp2", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/clincs.c" , 1222, (li->tmp2), (li->nc_alloc), sizeof(*(li->tmp2 ))); |
| 1223 | srenew(li->tmp3, li->nc_alloc)(li->tmp3) = save_realloc("li->tmp3", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/clincs.c" , 1223, (li->tmp3), (li->nc_alloc), sizeof(*(li->tmp3 ))); |
| 1224 | srenew(li->tmp4, li->nc_alloc)(li->tmp4) = save_realloc("li->tmp4", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/clincs.c" , 1224, (li->tmp4), (li->nc_alloc), sizeof(*(li->tmp4 ))); |
| 1225 | srenew(li->mlambda, li->nc_alloc)(li->mlambda) = save_realloc("li->mlambda", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/clincs.c" , 1225, (li->mlambda), (li->nc_alloc), sizeof(*(li-> mlambda))); |
| 1226 | if (li->ncg_triangle > 0) |
| 1227 | { |
| 1228 | /* This is allocating too much, but it is difficult to improve */ |
| 1229 | srenew(li->triangle, li->nc_alloc)(li->triangle) = save_realloc("li->triangle", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/clincs.c" , 1229, (li->triangle), (li->nc_alloc), sizeof(*(li-> triangle))); |
| 1230 | srenew(li->tri_bits, li->nc_alloc)(li->tri_bits) = save_realloc("li->tri_bits", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/clincs.c" , 1230, (li->tri_bits), (li->nc_alloc), sizeof(*(li-> tri_bits))); |
| 1231 | } |
| 1232 | } |
| 1233 | |
| 1234 | iatom = idef->il[F_CONSTR].iatoms; |
| 1235 | |
| 1236 | ncc_alloc = li->ncc_alloc; |
| 1237 | li->blnr[0] = 0; |
| 1238 | |
| 1239 | ni = idef->il[F_CONSTR].nr/3; |
| 1240 | |
| 1241 | con = 0; |
| 1242 | nconnect = 0; |
| 1243 | li->blnr[con] = nconnect; |
| 1244 | for (i = 0; i < ni; i++) |
| 1245 | { |
| 1246 | bLocal = TRUE1; |
Value stored to 'bLocal' is never read | |
| 1247 | type = iatom[3*i]; |
| 1248 | a1 = iatom[3*i+1]; |
| 1249 | a2 = iatom[3*i+2]; |
| 1250 | lenA = idef->iparams[type].constr.dA; |
| 1251 | lenB = idef->iparams[type].constr.dB; |
| 1252 | /* Skip the flexible constraints when not doing dynamics */ |
| 1253 | if (bDynamics || lenA != 0 || lenB != 0) |
| 1254 | { |
| 1255 | li->bllen0[con] = lenA; |
| 1256 | li->ddist[con] = lenB - lenA; |
| 1257 | /* Set the length to the topology A length */ |
| 1258 | li->bllen[con] = li->bllen0[con]; |
| 1259 | li->bla[2*con] = a1; |
| 1260 | li->bla[2*con+1] = a2; |
| 1261 | /* Construct the constraint connection matrix blbnb */ |
| 1262 | for (k = at2con.index[a1-start]; k < at2con.index[a1-start+1]; k++) |
| 1263 | { |
| 1264 | concon = at2con.a[k]; |
| 1265 | if (concon != i) |
| 1266 | { |
| 1267 | if (nconnect >= ncc_alloc) |
| 1268 | { |
| 1269 | ncc_alloc = over_alloc_small(nconnect+1)(int)(1.19*(nconnect+1) + 8000); |
| 1270 | srenew(li->blbnb, ncc_alloc)(li->blbnb) = save_realloc("li->blbnb", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/clincs.c" , 1270, (li->blbnb), (ncc_alloc), sizeof(*(li->blbnb))); |
| 1271 | } |
| 1272 | li->blbnb[nconnect++] = concon; |
| 1273 | } |
| 1274 | } |
| 1275 | for (k = at2con.index[a2-start]; k < at2con.index[a2-start+1]; k++) |
| 1276 | { |
| 1277 | concon = at2con.a[k]; |
| 1278 | if (concon != i) |
| 1279 | { |
| 1280 | if (nconnect+1 > ncc_alloc) |
| 1281 | { |
| 1282 | ncc_alloc = over_alloc_small(nconnect+1)(int)(1.19*(nconnect+1) + 8000); |
| 1283 | srenew(li->blbnb, ncc_alloc)(li->blbnb) = save_realloc("li->blbnb", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/clincs.c" , 1283, (li->blbnb), (ncc_alloc), sizeof(*(li->blbnb))); |
| 1284 | } |
| 1285 | li->blbnb[nconnect++] = concon; |
| 1286 | } |
| 1287 | } |
| 1288 | li->blnr[con+1] = nconnect; |
| 1289 | |
| 1290 | if (cr->dd == NULL((void*)0)) |
| 1291 | { |
| 1292 | /* Order the blbnb matrix to optimize memory access */ |
| 1293 | qsort(&(li->blbnb[li->blnr[con]]), li->blnr[con+1]-li->blnr[con], |
| 1294 | sizeof(li->blbnb[0]), int_comp); |
| 1295 | } |
| 1296 | /* Increase the constraint count */ |
| 1297 | con++; |
| 1298 | } |
| 1299 | } |
| 1300 | |
| 1301 | done_blocka(&at2con); |
| 1302 | |
| 1303 | /* This is the real number of constraints, |
| 1304 | * without dynamics the flexible constraints are not present. |
| 1305 | */ |
| 1306 | li->nc = con; |
| 1307 | |
| 1308 | li->ncc = li->blnr[con]; |
| 1309 | if (cr->dd == NULL((void*)0)) |
| 1310 | { |
| 1311 | /* Since the matrix is static, we can free some memory */ |
| 1312 | ncc_alloc = li->ncc; |
| 1313 | srenew(li->blbnb, ncc_alloc)(li->blbnb) = save_realloc("li->blbnb", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/clincs.c" , 1313, (li->blbnb), (ncc_alloc), sizeof(*(li->blbnb))); |
| 1314 | } |
| 1315 | |
| 1316 | if (ncc_alloc > li->ncc_alloc) |
| 1317 | { |
| 1318 | li->ncc_alloc = ncc_alloc; |
| 1319 | srenew(li->blmf, li->ncc_alloc)(li->blmf) = save_realloc("li->blmf", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/clincs.c" , 1319, (li->blmf), (li->ncc_alloc), sizeof(*(li->blmf ))); |
| 1320 | srenew(li->blmf1, li->ncc_alloc)(li->blmf1) = save_realloc("li->blmf1", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/clincs.c" , 1320, (li->blmf1), (li->ncc_alloc), sizeof(*(li->blmf1 ))); |
| 1321 | srenew(li->tmpncc, li->ncc_alloc)(li->tmpncc) = save_realloc("li->tmpncc", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/clincs.c" , 1321, (li->tmpncc), (li->ncc_alloc), sizeof(*(li-> tmpncc))); |
| 1322 | } |
| 1323 | |
| 1324 | if (debug) |
| 1325 | { |
| 1326 | fprintf(debug, "Number of constraints is %d, couplings %d\n", |
| 1327 | li->nc, li->ncc); |
| 1328 | } |
| 1329 | |
| 1330 | if (li->nth == 1) |
| 1331 | { |
| 1332 | li->th[0].b0 = 0; |
| 1333 | li->th[0].b1 = li->nc; |
| 1334 | } |
| 1335 | else |
| 1336 | { |
| 1337 | lincs_thread_setup(li, md->nr); |
| 1338 | } |
| 1339 | |
| 1340 | set_lincs_matrix(li, md->invmass, md->lambda); |
| 1341 | } |
| 1342 | |
| 1343 | static void lincs_warning(FILE *fplog, |
| 1344 | gmx_domdec_t *dd, rvec *x, rvec *xprime, t_pbc *pbc, |
| 1345 | int ncons, int *bla, real *bllen, real wangle, |
| 1346 | int maxwarn, int *warncount) |
| 1347 | { |
| 1348 | int b, i, j; |
| 1349 | rvec v0, v1; |
| 1350 | real wfac, d0, d1, cosine; |
| 1351 | char buf[STRLEN4096]; |
| 1352 | |
| 1353 | wfac = cos(DEG2RAD(3.14159265358979323846/180.0)*wangle); |
| 1354 | |
| 1355 | sprintf(buf, "bonds that rotated more than %g degrees:\n" |
| 1356 | " atom 1 atom 2 angle previous, current, constraint length\n", |
| 1357 | wangle); |
| 1358 | fprintf(stderrstderr, "%s", buf); |
| 1359 | if (fplog) |
| 1360 | { |
| 1361 | fprintf(fplog, "%s", buf); |
| 1362 | } |
| 1363 | |
| 1364 | for (b = 0; b < ncons; b++) |
| 1365 | { |
| 1366 | i = bla[2*b]; |
| 1367 | j = bla[2*b+1]; |
| 1368 | if (pbc) |
| 1369 | { |
| 1370 | pbc_dx_aiuc(pbc, x[i], x[j], v0); |
| 1371 | pbc_dx_aiuc(pbc, xprime[i], xprime[j], v1); |
| 1372 | } |
| 1373 | else |
| 1374 | { |
| 1375 | rvec_sub(x[i], x[j], v0); |
| 1376 | rvec_sub(xprime[i], xprime[j], v1); |
| 1377 | } |
| 1378 | d0 = norm(v0); |
| 1379 | d1 = norm(v1); |
| 1380 | cosine = iprod(v0, v1)/(d0*d1); |
| 1381 | if (cosine < wfac) |
| 1382 | { |
| 1383 | sprintf(buf, " %6d %6d %5.1f %8.4f %8.4f %8.4f\n", |
| 1384 | ddglatnr(dd, i), ddglatnr(dd, j), |
| 1385 | RAD2DEG(180.0/3.14159265358979323846)*acos(cosine), d0, d1, bllen[b]); |
| 1386 | fprintf(stderrstderr, "%s", buf); |
| 1387 | if (fplog) |
| 1388 | { |
| 1389 | fprintf(fplog, "%s", buf); |
| 1390 | } |
| 1391 | if (!gmx_isfinite(d1)) |
| 1392 | { |
| 1393 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/clincs.c", 1393, "Bond length not finite."); |
| 1394 | } |
| 1395 | |
| 1396 | (*warncount)++; |
| 1397 | } |
| 1398 | } |
| 1399 | if (*warncount > maxwarn) |
| 1400 | { |
| 1401 | too_many_constraint_warnings(econtLINCS, *warncount); |
| 1402 | } |
| 1403 | } |
| 1404 | |
| 1405 | static void cconerr(gmx_domdec_t *dd, |
| 1406 | int ncons, int *bla, real *bllen, rvec *x, t_pbc *pbc, |
| 1407 | real *ncons_loc, real *ssd, real *max, int *imax) |
| 1408 | { |
| 1409 | real len, d, ma, ssd2, r2; |
| 1410 | int *nlocat, count, b, im; |
| 1411 | rvec dx; |
| 1412 | |
| 1413 | if (dd && dd->constraints) |
| 1414 | { |
| 1415 | nlocat = dd_constraints_nlocalatoms(dd); |
| 1416 | } |
| 1417 | else |
| 1418 | { |
| 1419 | nlocat = 0; |
| 1420 | } |
| 1421 | |
| 1422 | ma = 0; |
| 1423 | ssd2 = 0; |
| 1424 | im = 0; |
| 1425 | count = 0; |
| 1426 | for (b = 0; b < ncons; b++) |
| 1427 | { |
| 1428 | if (pbc) |
| 1429 | { |
| 1430 | pbc_dx_aiuc(pbc, x[bla[2*b]], x[bla[2*b+1]], dx); |
| 1431 | } |
| 1432 | else |
| 1433 | { |
| 1434 | rvec_sub(x[bla[2*b]], x[bla[2*b+1]], dx); |
| 1435 | } |
| 1436 | r2 = norm2(dx); |
| 1437 | len = r2*gmx_invsqrt(r2)gmx_software_invsqrt(r2); |
| 1438 | d = fabs(len/bllen[b]-1); |
| 1439 | if (d > ma && (nlocat == NULL((void*)0) || nlocat[b])) |
| 1440 | { |
| 1441 | ma = d; |
| 1442 | im = b; |
| 1443 | } |
| 1444 | if (nlocat == NULL((void*)0)) |
| 1445 | { |
| 1446 | ssd2 += d*d; |
| 1447 | count++; |
| 1448 | } |
| 1449 | else |
| 1450 | { |
| 1451 | ssd2 += nlocat[b]*d*d; |
| 1452 | count += nlocat[b]; |
| 1453 | } |
| 1454 | } |
| 1455 | |
| 1456 | *ncons_loc = (nlocat ? 0.5 : 1)*count; |
| 1457 | *ssd = (nlocat ? 0.5 : 1)*ssd2; |
| 1458 | *max = ma; |
| 1459 | *imax = im; |
| 1460 | } |
| 1461 | |
| 1462 | gmx_bool constrain_lincs(FILE *fplog, gmx_bool bLog, gmx_bool bEner, |
| 1463 | t_inputrec *ir, |
| 1464 | gmx_int64_t step, |
| 1465 | struct gmx_lincsdata *lincsd, t_mdatoms *md, |
| 1466 | t_commrec *cr, |
| 1467 | rvec *x, rvec *xprime, rvec *min_proj, |
| 1468 | matrix box, t_pbc *pbc, |
| 1469 | real lambda, real *dvdlambda, |
| 1470 | real invdt, rvec *v, |
| 1471 | gmx_bool bCalcVir, tensor vir_r_m_dr, |
| 1472 | int econq, |
| 1473 | t_nrnb *nrnb, |
| 1474 | int maxwarn, int *warncount) |
| 1475 | { |
| 1476 | char buf[STRLEN4096], buf2[22], buf3[STRLEN4096]; |
| 1477 | int i, warn, p_imax, error; |
| 1478 | real ncons_loc, p_ssd, p_max = 0; |
| 1479 | rvec dx; |
| 1480 | gmx_bool bOK; |
| 1481 | |
| 1482 | bOK = TRUE1; |
| 1483 | |
| 1484 | if (lincsd->nc == 0 && cr->dd == NULL((void*)0)) |
| 1485 | { |
| 1486 | if (bLog || bEner) |
| 1487 | { |
| 1488 | lincsd->rmsd_data[0] = 0; |
| 1489 | if (ir->eI == eiSD2 && v == NULL((void*)0)) |
| 1490 | { |
| 1491 | i = 2; |
| 1492 | } |
| 1493 | else |
| 1494 | { |
| 1495 | i = 1; |
| 1496 | } |
| 1497 | lincsd->rmsd_data[i] = 0; |
| 1498 | } |
| 1499 | |
| 1500 | return bOK; |
| 1501 | } |
| 1502 | |
| 1503 | if (econq == econqCoord) |
| 1504 | { |
| 1505 | if (ir->efep != efepNO) |
| 1506 | { |
| 1507 | if (md->nMassPerturbed && lincsd->matlam != md->lambda) |
| 1508 | { |
| 1509 | set_lincs_matrix(lincsd, md->invmass, md->lambda); |
| 1510 | } |
| 1511 | |
| 1512 | for (i = 0; i < lincsd->nc; i++) |
| 1513 | { |
| 1514 | lincsd->bllen[i] = lincsd->bllen0[i] + lambda*lincsd->ddist[i]; |
| 1515 | } |
| 1516 | } |
| 1517 | |
| 1518 | if (lincsd->ncg_flex) |
| 1519 | { |
| 1520 | /* Set the flexible constraint lengths to the old lengths */ |
| 1521 | if (pbc != NULL((void*)0)) |
| 1522 | { |
| 1523 | for (i = 0; i < lincsd->nc; i++) |
| 1524 | { |
| 1525 | if (lincsd->bllen[i] == 0) |
| 1526 | { |
| 1527 | pbc_dx_aiuc(pbc, x[lincsd->bla[2*i]], x[lincsd->bla[2*i+1]], dx); |
| 1528 | lincsd->bllen[i] = norm(dx); |
| 1529 | } |
| 1530 | } |
| 1531 | } |
| 1532 | else |
| 1533 | { |
| 1534 | for (i = 0; i < lincsd->nc; i++) |
| 1535 | { |
| 1536 | if (lincsd->bllen[i] == 0) |
| 1537 | { |
| 1538 | lincsd->bllen[i] = |
| 1539 | sqrt(distance2(x[lincsd->bla[2*i]], |
| 1540 | x[lincsd->bla[2*i+1]])); |
| 1541 | } |
| 1542 | } |
| 1543 | } |
| 1544 | } |
| 1545 | |
| 1546 | if (bLog && fplog) |
| 1547 | { |
| 1548 | cconerr(cr->dd, lincsd->nc, lincsd->bla, lincsd->bllen, xprime, pbc, |
| 1549 | &ncons_loc, &p_ssd, &p_max, &p_imax); |
| 1550 | } |
| 1551 | |
| 1552 | /* This warn var can be updated by multiple threads |
| 1553 | * at the same time. But as we only need to detect |
| 1554 | * if a warning occured or not, this is not an issue. |
| 1555 | */ |
| 1556 | warn = -1; |
| 1557 | |
| 1558 | /* The OpenMP parallel region of constrain_lincs for coords */ |
| 1559 | #pragma omp parallel num_threads(lincsd->nth) |
| 1560 | { |
| 1561 | int th = gmx_omp_get_thread_num(); |
| 1562 | |
| 1563 | clear_mat(lincsd->th[th].vir_r_m_dr); |
| 1564 | |
| 1565 | do_lincs(x, xprime, box, pbc, lincsd, th, |
| 1566 | md->invmass, cr, |
| 1567 | bCalcVir || (ir->efep != efepNO), |
| 1568 | ir->LincsWarnAngle, &warn, |
| 1569 | invdt, v, bCalcVir, |
| 1570 | th == 0 ? vir_r_m_dr : lincsd->th[th].vir_r_m_dr); |
| 1571 | } |
| 1572 | |
| 1573 | if (ir->efep != efepNO) |
| 1574 | { |
| 1575 | real dt_2, dvdl = 0; |
| 1576 | |
| 1577 | dt_2 = 1.0/(ir->delta_t*ir->delta_t); |
| 1578 | for (i = 0; (i < lincsd->nc); i++) |
| 1579 | { |
| 1580 | dvdl -= lincsd->mlambda[i]*dt_2*lincsd->ddist[i]; |
| 1581 | } |
| 1582 | *dvdlambda += dvdl; |
| 1583 | } |
| 1584 | |
| 1585 | if (bLog && fplog && lincsd->nc > 0) |
| 1586 | { |
| 1587 | fprintf(fplog, " Rel. Constraint Deviation: RMS MAX between atoms\n"); |
| 1588 | fprintf(fplog, " Before LINCS %.6f %.6f %6d %6d\n", |
| 1589 | sqrt(p_ssd/ncons_loc), p_max, |
| 1590 | ddglatnr(cr->dd, lincsd->bla[2*p_imax]), |
| 1591 | ddglatnr(cr->dd, lincsd->bla[2*p_imax+1])); |
| 1592 | } |
| 1593 | if (bLog || bEner) |
| 1594 | { |
| 1595 | cconerr(cr->dd, lincsd->nc, lincsd->bla, lincsd->bllen, xprime, pbc, |
| 1596 | &ncons_loc, &p_ssd, &p_max, &p_imax); |
| 1597 | /* Check if we are doing the second part of SD */ |
| 1598 | if (ir->eI == eiSD2 && v == NULL((void*)0)) |
| 1599 | { |
| 1600 | i = 2; |
| 1601 | } |
| 1602 | else |
| 1603 | { |
| 1604 | i = 1; |
| 1605 | } |
| 1606 | lincsd->rmsd_data[0] = ncons_loc; |
| 1607 | lincsd->rmsd_data[i] = p_ssd; |
| 1608 | } |
| 1609 | else |
| 1610 | { |
| 1611 | lincsd->rmsd_data[0] = 0; |
| 1612 | lincsd->rmsd_data[1] = 0; |
| 1613 | lincsd->rmsd_data[2] = 0; |
| 1614 | } |
| 1615 | if (bLog && fplog && lincsd->nc > 0) |
| 1616 | { |
| 1617 | fprintf(fplog, |
| 1618 | " After LINCS %.6f %.6f %6d %6d\n\n", |
| 1619 | sqrt(p_ssd/ncons_loc), p_max, |
| 1620 | ddglatnr(cr->dd, lincsd->bla[2*p_imax]), |
| 1621 | ddglatnr(cr->dd, lincsd->bla[2*p_imax+1])); |
| 1622 | } |
| 1623 | |
| 1624 | if (warn >= 0) |
| 1625 | { |
| 1626 | if (maxwarn >= 0) |
| 1627 | { |
| 1628 | cconerr(cr->dd, lincsd->nc, lincsd->bla, lincsd->bllen, xprime, pbc, |
| 1629 | &ncons_loc, &p_ssd, &p_max, &p_imax); |
| 1630 | if (MULTISIM(cr)((cr)->ms)) |
| 1631 | { |
| 1632 | sprintf(buf3, " in simulation %d", cr->ms->sim); |
| 1633 | } |
| 1634 | else |
| 1635 | { |
| 1636 | buf3[0] = 0; |
| 1637 | } |
| 1638 | sprintf(buf, "\nStep %s, time %g (ps) LINCS WARNING%s\n" |
| 1639 | "relative constraint deviation after LINCS:\n" |
| 1640 | "rms %.6f, max %.6f (between atoms %d and %d)\n", |
| 1641 | gmx_step_str(step, buf2), ir->init_t+step*ir->delta_t, |
| 1642 | buf3, |
| 1643 | sqrt(p_ssd/ncons_loc), p_max, |
| 1644 | ddglatnr(cr->dd, lincsd->bla[2*p_imax]), |
| 1645 | ddglatnr(cr->dd, lincsd->bla[2*p_imax+1])); |
| 1646 | if (fplog) |
| 1647 | { |
| 1648 | fprintf(fplog, "%s", buf); |
| 1649 | } |
| 1650 | fprintf(stderrstderr, "%s", buf); |
| 1651 | lincs_warning(fplog, cr->dd, x, xprime, pbc, |
| 1652 | lincsd->nc, lincsd->bla, lincsd->bllen, |
| 1653 | ir->LincsWarnAngle, maxwarn, warncount); |
| 1654 | } |
| 1655 | bOK = (p_max < 0.5); |
| 1656 | } |
| 1657 | |
| 1658 | if (lincsd->ncg_flex) |
| 1659 | { |
| 1660 | for (i = 0; (i < lincsd->nc); i++) |
| 1661 | { |
| 1662 | if (lincsd->bllen0[i] == 0 && lincsd->ddist[i] == 0) |
| 1663 | { |
| 1664 | lincsd->bllen[i] = 0; |
| 1665 | } |
| 1666 | } |
| 1667 | } |
| 1668 | } |
| 1669 | else |
| 1670 | { |
| 1671 | /* The OpenMP parallel region of constrain_lincs for derivatives */ |
| 1672 | #pragma omp parallel num_threads(lincsd->nth) |
| 1673 | { |
| 1674 | int th = gmx_omp_get_thread_num(); |
| 1675 | |
| 1676 | do_lincsp(x, xprime, min_proj, pbc, lincsd, th, |
| 1677 | md->invmass, econq, ir->efep != efepNO ? dvdlambda : NULL((void*)0), |
| 1678 | bCalcVir, th == 0 ? vir_r_m_dr : lincsd->th[th].vir_r_m_dr); |
| 1679 | } |
| 1680 | } |
| 1681 | |
| 1682 | if (bCalcVir && lincsd->nth > 1) |
| 1683 | { |
| 1684 | for (i = 1; i < lincsd->nth; i++) |
| 1685 | { |
| 1686 | m_add(vir_r_m_dr, lincsd->th[i].vir_r_m_dr, vir_r_m_dr); |
| 1687 | } |
| 1688 | } |
| 1689 | |
| 1690 | /* count assuming nit=1 */ |
| 1691 | inc_nrnb(nrnb, eNR_LINCS, lincsd->nc)(nrnb)->n[eNR_LINCS] += lincsd->nc; |
| 1692 | inc_nrnb(nrnb, eNR_LINCSMAT, (2+lincsd->nOrder)*lincsd->ncc)(nrnb)->n[eNR_LINCSMAT] += (2+lincsd->nOrder)*lincsd-> ncc; |
| 1693 | if (lincsd->ntriangle > 0) |
| 1694 | { |
| 1695 | inc_nrnb(nrnb, eNR_LINCSMAT, lincsd->nOrder*lincsd->ncc_triangle)(nrnb)->n[eNR_LINCSMAT] += lincsd->nOrder*lincsd->ncc_triangle; |
| 1696 | } |
| 1697 | if (v) |
| 1698 | { |
| 1699 | inc_nrnb(nrnb, eNR_CONSTR_V, lincsd->nc*2)(nrnb)->n[eNR_CONSTR_V] += lincsd->nc*2; |
| 1700 | } |
| 1701 | if (bCalcVir) |
| 1702 | { |
| 1703 | inc_nrnb(nrnb, eNR_CONSTR_VIR, lincsd->nc)(nrnb)->n[eNR_CONSTR_VIR] += lincsd->nc; |
| 1704 | } |
| 1705 | |
| 1706 | return bOK; |
| 1707 | } |