Bug Summary

File:gromacs/mdlib/shakef.c
Location:line 326, column 13
Description:Value stored to 'toler' is never read

Annotated Source Code

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#ifdef HAVE_CONFIG_H1
38#include <config.h>
39#endif
40
41#include <math.h>
42#include "typedefs.h"
43#include "gromacs/utility/smalloc.h"
44#include "pbc.h"
45#include "txtdump.h"
46#include "gromacs/math/vec.h"
47#include "nrnb.h"
48#include "constr.h"
49
50typedef struct gmx_shakedata
51{
52 rvec *rij;
53 real *M2;
54 real *tt;
55 real *dist2;
56 int nalloc;
57 /* SOR stuff */
58 real delta;
59 real omega;
60 real gamma;
61} t_gmx_shakedata;
62
63gmx_shakedata_t shake_init()
64{
65 gmx_shakedata_t d;
66
67 snew(d, 1)(d) = save_calloc("d", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/shakef.c"
, 67, (1), sizeof(*(d)));
68
69 d->nalloc = 0;
70 d->rij = NULL((void*)0);
71 d->M2 = NULL((void*)0);
72 d->tt = NULL((void*)0);
73 d->dist2 = NULL((void*)0);
74
75 /* SOR initialization */
76 d->delta = 0.1;
77 d->omega = 1.0;
78 d->gamma = 1000000;
79
80 return d;
81}
82
83static void pv(FILE *log, char *s, rvec x)
84{
85 int m;
86
87 fprintf(log, "%5s:", s);
88 for (m = 0; (m < DIM3); m++)
89 {
90 fprintf(log, " %10.3f", x[m]);
91 }
92 fprintf(log, "\n");
93 fflush(log);
94}
95
96void cshake(atom_id iatom[], int ncon, int *nnit, int maxnit,
97 real dist2[], real xp[], real rij[], real m2[], real omega,
98 real invmass[], real tt[], real lagr[], int *nerror)
99{
100 /*
101 * r.c. van schaik and w.f. van gunsteren
102 * eth zuerich
103 * june 1992
104 * Adapted for use with Gromacs by David van der Spoel november 92 and later.
105 */
106 /* default should be increased! MRS 8/4/2009 */
107 const real mytol = 1e-10;
108
109 int ll, i, j, i3, j3, l3;
110 int ix, iy, iz, jx, jy, jz;
111 real toler, rpij2, rrpr, tx, ty, tz, diff, acor, im, jm;
112 real xh, yh, zh, rijx, rijy, rijz;
113 real tix, tiy, tiz;
114 real tjx, tjy, tjz;
115 int nit, error, nconv;
116 real iconvf;
117
118 error = 0;
119 nconv = 1;
120 for (nit = 0; (nit < maxnit) && (nconv != 0) && (error == 0); nit++)
121 {
122 nconv = 0;
123 for (ll = 0; (ll < ncon) && (error == 0); ll++)
124 {
125 l3 = 3*ll;
126 rijx = rij[l3+XX0];
127 rijy = rij[l3+YY1];
128 rijz = rij[l3+ZZ2];
129 i = iatom[l3+1];
130 j = iatom[l3+2];
131 i3 = 3*i;
132 j3 = 3*j;
133 ix = i3+XX0;
134 iy = i3+YY1;
135 iz = i3+ZZ2;
136 jx = j3+XX0;
137 jy = j3+YY1;
138 jz = j3+ZZ2;
139
140 tx = xp[ix]-xp[jx];
141 ty = xp[iy]-xp[jy];
142 tz = xp[iz]-xp[jz];
143 rpij2 = tx*tx+ty*ty+tz*tz;
144 toler = dist2[ll];
145 diff = toler-rpij2;
146
147 /* iconvf is less than 1 when the error is smaller than a bound */
148 /* But if tt is too big, then it will result in looping in iconv */
149
150 iconvf = fabs(diff)*tt[ll];
151
152 if (iconvf > 1)
153 {
154 nconv = iconvf;
155 rrpr = rijx*tx+rijy*ty+rijz*tz;
156
157 if (rrpr < toler*mytol)
158 {
159 error = ll+1;
160 }
161 else
162 {
163 acor = omega*diff*m2[ll]/rrpr;
164 lagr[ll] += acor;
165 xh = rijx*acor;
166 yh = rijy*acor;
167 zh = rijz*acor;
168 im = invmass[i];
169 jm = invmass[j];
170 xp[ix] += xh*im;
171 xp[iy] += yh*im;
172 xp[iz] += zh*im;
173 xp[jx] -= xh*jm;
174 xp[jy] -= yh*jm;
175 xp[jz] -= zh*jm;
176 }
177 }
178 }
179 }
180 *nnit = nit;
181 *nerror = error;
182}
183
184int vec_shakef(FILE *fplog, gmx_shakedata_t shaked,
185 real invmass[], int ncon,
186 t_iparams ip[], t_iatom *iatom,
187 real tol, rvec x[], rvec prime[], real omega,
188 gmx_bool bFEP, real lambda, real lagr[],
189 real invdt, rvec *v,
190 gmx_bool bCalcVir, tensor vir_r_m_dr, int econq,
191 t_vetavars *vetavar)
192{
193 rvec *rij;
194 real *M2, *tt, *dist2;
195 int maxnit = 1000;
196 int nit = 0, ll, i, j, type;
197 t_iatom *ia;
198 real L1, tol2, toler;
199 real mm = 0., tmp;
200 int error = 0;
201 real g, vscale, rscale, rvscale;
202
203 if (ncon > shaked->nalloc)
204 {
205 shaked->nalloc = over_alloc_dd(ncon);
206 srenew(shaked->rij, shaked->nalloc)(shaked->rij) = save_realloc("shaked->rij", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/shakef.c"
, 206, (shaked->rij), (shaked->nalloc), sizeof(*(shaked
->rij)));
207 srenew(shaked->M2, shaked->nalloc)(shaked->M2) = save_realloc("shaked->M2", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/shakef.c"
, 207, (shaked->M2), (shaked->nalloc), sizeof(*(shaked->
M2)));
208 srenew(shaked->tt, shaked->nalloc)(shaked->tt) = save_realloc("shaked->tt", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/shakef.c"
, 208, (shaked->tt), (shaked->nalloc), sizeof(*(shaked->
tt)));
209 srenew(shaked->dist2, shaked->nalloc)(shaked->dist2) = save_realloc("shaked->dist2", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/shakef.c"
, 209, (shaked->dist2), (shaked->nalloc), sizeof(*(shaked
->dist2)));
210 }
211 rij = shaked->rij;
212 M2 = shaked->M2;
213 tt = shaked->tt;
214 dist2 = shaked->dist2;
215
216 L1 = 1.0-lambda;
217 tol2 = 2.0*tol;
218 ia = iatom;
219 for (ll = 0; (ll < ncon); ll++, ia += 3)
220 {
221 type = ia[0];
222 i = ia[1];
223 j = ia[2];
224
225 mm = 2*(invmass[i]+invmass[j]);
226 rij[ll][XX0] = x[i][XX0]-x[j][XX0];
227 rij[ll][YY1] = x[i][YY1]-x[j][YY1];
228 rij[ll][ZZ2] = x[i][ZZ2]-x[j][ZZ2];
229 M2[ll] = 1.0/mm;
230 if (bFEP)
231 {
232 toler = sqr(L1*ip[type].constr.dA + lambda*ip[type].constr.dB);
233 }
234 else
235 {
236 toler = sqr(ip[type].constr.dA);
237 }
238 dist2[ll] = toler;
239 tt[ll] = 1.0/(toler*tol2);
240 }
241
242 switch (econq)
243 {
244 case econqCoord:
245 cshake(iatom, ncon, &nit, maxnit, dist2, prime[0], rij[0], M2, omega, invmass, tt, lagr, &error);
246 break;
247 case econqVeloc:
248 crattle(iatom, ncon, &nit, maxnit, dist2, prime[0], rij[0], M2, omega, invmass, tt, lagr, &error, invdt, vetavar);
249 break;
250 }
251
252 if (nit >= maxnit)
253 {
254 if (fplog)
255 {
256 fprintf(fplog, "Shake did not converge in %d steps\n", maxnit);
257 }
258 fprintf(stderrstderr, "Shake did not converge in %d steps\n", maxnit);
259 nit = 0;
260 }
261 else if (error != 0)
262 {
263 if (fplog)
264 {
265 fprintf(fplog, "Inner product between old and new vector <= 0.0!\n"
266 "constraint #%d atoms %u and %u\n",
267 error-1, iatom[3*(error-1)+1]+1, iatom[3*(error-1)+2]+1);
268 }
269 fprintf(stderrstderr, "Inner product between old and new vector <= 0.0!\n"
270 "constraint #%d atoms %u and %u\n",
271 error-1, iatom[3*(error-1)+1]+1, iatom[3*(error-1)+2]+1);
272 nit = 0;
273 }
274
275 /* Constraint virial and correct the lagrange multipliers for the length */
276
277 ia = iatom;
278
279 for (ll = 0; (ll < ncon); ll++, ia += 3)
280 {
281
282 if ((econq == econqCoord) && v != NULL((void*)0))
283 {
284 /* Correct the velocities */
285 mm = lagr[ll]*invmass[ia[1]]*invdt/vetavar->rscale;
286 for (i = 0; i < DIM3; i++)
287 {
288 v[ia[1]][i] += mm*rij[ll][i];
289 }
290 mm = lagr[ll]*invmass[ia[2]]*invdt/vetavar->rscale;
291 for (i = 0; i < DIM3; i++)
292 {
293 v[ia[2]][i] -= mm*rij[ll][i];
294 }
295 /* 16 flops */
296 }
297
298 /* constraint virial */
299 if (bCalcVir)
300 {
301 if (econq == econqCoord)
302 {
303 mm = lagr[ll]/vetavar->rvscale;
304 }
305 if (econq == econqVeloc)
306 {
307 mm = lagr[ll]/(vetavar->vscale*vetavar->vscale_nhc[0]);
308 }
309 for (i = 0; i < DIM3; i++)
310 {
311 tmp = mm*rij[ll][i];
312 for (j = 0; j < DIM3; j++)
313 {
314 vir_r_m_dr[i][j] -= tmp*rij[ll][j];
315 }
316 }
317 /* 21 flops */
318 }
319
320 /* Correct the lagrange multipliers for the length */
321 /* (more details would be useful here . . . )*/
322
323 type = ia[0];
324 if (bFEP)
325 {
326 toler = L1*ip[type].constr.dA + lambda*ip[type].constr.dB;
Value stored to 'toler' is never read
327 }
328 else
329 {
330 toler = ip[type].constr.dA;
331 lagr[ll] *= toler;
332 }
333 }
334
335 return nit;
336}
337
338static void check_cons(FILE *log, int nc, rvec x[], rvec prime[], rvec v[],
339 t_iparams ip[], t_iatom *iatom,
340 real invmass[], int econq)
341{
342 t_iatom *ia;
343 int ai, aj;
344 int i;
345 real d, dp;
346 rvec dx, dv;
347
348 fprintf(log,
349 " i mi j mj before after should be\n");
350 ia = iatom;
351 for (i = 0; (i < nc); i++, ia += 3)
352 {
353 ai = ia[1];
354 aj = ia[2];
355 rvec_sub(x[ai], x[aj], dx);
356 d = norm(dx);
357
358 switch (econq)
359 {
360 case econqCoord:
361 rvec_sub(prime[ai], prime[aj], dx);
362 dp = norm(dx);
363 fprintf(log, "%5d %5.2f %5d %5.2f %10.5f %10.5f %10.5f\n",
364 ai+1, 1.0/invmass[ai],
365 aj+1, 1.0/invmass[aj], d, dp, ip[ia[0]].constr.dA);
366 break;
367 case econqVeloc:
368 rvec_sub(v[ai], v[aj], dv);
369 d = iprod(dx, dv);
370 rvec_sub(prime[ai], prime[aj], dv);
371 dp = iprod(dx, dv);
372 fprintf(log, "%5d %5.2f %5d %5.2f %10.5f %10.5f %10.5f\n",
373 ai+1, 1.0/invmass[ai],
374 aj+1, 1.0/invmass[aj], d, dp, 0.);
375 break;
376 }
377 }
378}
379
380gmx_bool bshakef(FILE *log, gmx_shakedata_t shaked,
381 real invmass[], int nblocks, int sblock[],
382 t_idef *idef, t_inputrec *ir, rvec x_s[], rvec prime[],
383 t_nrnb *nrnb, real *lagr, real lambda, real *dvdlambda,
384 real invdt, rvec *v, gmx_bool bCalcVir, tensor vir_r_m_dr,
385 gmx_bool bDumpOnError, int econq, t_vetavars *vetavar)
386{
387 t_iatom *iatoms;
388 real *lam, dt_2, dvdl;
389 int i, n0, ncons, blen, type;
390 int tnit = 0, trij = 0;
391
392#ifdef DEBUG
393 fprintf(log, "nblocks=%d, sblock[0]=%d\n", nblocks, sblock[0]);
394#endif
395
396 ncons = idef->il[F_CONSTR].nr/3;
397
398 for (i = 0; i < ncons; i++)
399 {
400 lagr[i] = 0;
401 }
402
403 iatoms = &(idef->il[F_CONSTR].iatoms[sblock[0]]);
404 lam = lagr;
405 for (i = 0; (i < nblocks); )
406 {
407 blen = (sblock[i+1]-sblock[i]);
408 blen /= 3;
409 n0 = vec_shakef(log, shaked, invmass, blen, idef->iparams,
410 iatoms, ir->shake_tol, x_s, prime, shaked->omega,
411 ir->efep != efepNO, lambda, lam, invdt, v, bCalcVir, vir_r_m_dr,
412 econq, vetavar);
413
414#ifdef DEBUGSHAKE
415 check_cons(log, blen, x_s, prime, v, idef->iparams, iatoms, invmass, econq);
416#endif
417
418 if (n0 == 0)
419 {
420 if (bDumpOnError && log)
421 {
422 {
423 check_cons(log, blen, x_s, prime, v, idef->iparams, iatoms, invmass, econq);
424 }
425 }
426 return FALSE0;
427 }
428 tnit += n0*blen;
429 trij += blen;
430 iatoms += 3*blen; /* Increment pointer! */
431 lam += blen;
432 i++;
433 }
434 /* only for position part? */
435 if (econq == econqCoord)
436 {
437 if (ir->efep != efepNO)
438 {
439 real bondA, bondB;
440 dt_2 = 1/sqr(ir->delta_t);
441 dvdl = 0;
442 for (i = 0; i < ncons; i++)
443 {
444 type = idef->il[F_CONSTR].iatoms[3*i];
445
446 /* dh/dl contribution from constraint force is dh/dr (constraint force) dot dr/dl */
447 /* constraint force is -\sum_i lagr_i* d(constraint)/dr, with constrant = r^2-d^2 */
448 /* constraint force is -\sum_i lagr_i* 2 r */
449 /* so dh/dl = -\sum_i lagr_i* 2 r * dr/dl */
450 /* However, by comparison with lincs and with
451 comparison with a full thermodynamics cycle (see
452 redmine issue #1255), this is off by a factor of
453 two -- the 2r should apparently just be r. Further
454 investigation should be done at some point to
455 understand why and see if there is something deeper
456 we are missing */
457
458 bondA = idef->iparams[type].constr.dA;
459 bondB = idef->iparams[type].constr.dB;
460 dvdl += lagr[i] * dt_2 * ((1.0-lambda)*bondA + lambda*bondB) * (bondB-bondA);
461 }
462 *dvdlambda += dvdl;
463 }
464 }
465#ifdef DEBUG
466 fprintf(log, "tnit: %5d omega: %10.5f\n", tnit, omega);
467#endif
468 if (ir->bShakeSOR)
469 {
470 if (tnit > shaked->gamma)
471 {
472 shaked->delta *= -0.5;
473 }
474 shaked->omega += shaked->delta;
475 shaked->gamma = tnit;
476 }
477 inc_nrnb(nrnb, eNR_SHAKE, tnit)(nrnb)->n[eNR_SHAKE] += tnit;
478 inc_nrnb(nrnb, eNR_SHAKE_RIJ, trij)(nrnb)->n[eNR_SHAKE_RIJ] += trij;
479 if (v)
480 {
481 inc_nrnb(nrnb, eNR_CONSTR_V, trij*2)(nrnb)->n[eNR_CONSTR_V] += trij*2;
482 }
483 if (bCalcVir)
484 {
485 inc_nrnb(nrnb, eNR_CONSTR_VIR, trij)(nrnb)->n[eNR_CONSTR_VIR] += trij;
486 }
487
488 return TRUE1;
489}
490
491void crattle(atom_id iatom[], int ncon, int *nnit, int maxnit,
492 real dist2[], real vp[], real rij[], real m2[], real omega,
493 real invmass[], real tt[], real lagr[], int *nerror, real invdt, t_vetavars *vetavar)
494{
495 /*
496 * r.c. van schaik and w.f. van gunsteren
497 * eth zuerich
498 * june 1992
499 * Adapted for use with Gromacs by David van der Spoel november 92 and later.
500 * rattle added by M.R. Shirts, April 2004, from code written by Jay Ponder in TINKER
501 * second part of rattle algorithm
502 */
503
504 const real mytol = 1e-10;
505
506 int ll, i, j, i3, j3, l3, ii;
507 int ix, iy, iz, jx, jy, jz;
508 real toler, rijd, vpijd, vx, vy, vz, diff, acor, xdotd, fac, im, jm, imdt, jmdt;
509 real xh, yh, zh, rijx, rijy, rijz;
510 real tix, tiy, tiz;
511 real tjx, tjy, tjz;
512 int nit, error, nconv;
513 real veta, vscale_nhc, iconvf;
514
515 veta = vetavar->veta;
516 vscale_nhc = vetavar->vscale_nhc[0]; /* for now, just use the first state */
517
518 error = 0;
519 nconv = 1;
520 for (nit = 0; (nit < maxnit) && (nconv != 0) && (error == 0); nit++)
521 {
522 nconv = 0;
523 for (ll = 0; (ll < ncon) && (error == 0); ll++)
524 {
525 l3 = 3*ll;
526 rijx = rij[l3+XX0];
527 rijy = rij[l3+YY1];
528 rijz = rij[l3+ZZ2];
529 i = iatom[l3+1];
530 j = iatom[l3+2];
531 i3 = 3*i;
532 j3 = 3*j;
533 ix = i3+XX0;
534 iy = i3+YY1;
535 iz = i3+ZZ2;
536 jx = j3+XX0;
537 jy = j3+YY1;
538 jz = j3+ZZ2;
539 vx = vp[ix]-vp[jx];
540 vy = vp[iy]-vp[jy];
541 vz = vp[iz]-vp[jz];
542
543 vpijd = vx*rijx+vy*rijy+vz*rijz;
544 toler = dist2[ll];
545 /* this is r(t+dt) \dotproduct \dot{r}(t+dt) */
546 xdotd = vpijd*vscale_nhc + veta*toler;
547
548 /* iconv is zero when the error is smaller than a bound */
549 iconvf = fabs(xdotd)*(tt[ll]/invdt);
550
551 if (iconvf > 1)
552 {
553 nconv = iconvf;
554 fac = omega*2.0*m2[ll]/toler;
555 acor = -fac*xdotd;
556 lagr[ll] += acor;
557
558 xh = rijx*acor;
559 yh = rijy*acor;
560 zh = rijz*acor;
561
562 im = invmass[i]/vscale_nhc;
563 jm = invmass[j]/vscale_nhc;
564
565 vp[ix] += xh*im;
566 vp[iy] += yh*im;
567 vp[iz] += zh*im;
568 vp[jx] -= xh*jm;
569 vp[jy] -= yh*jm;
570 vp[jz] -= zh*jm;
571 }
572 }
573 }
574 *nnit = nit;
575 *nerror = error;
576}