Bug Summary

File:gromacs/gmxlib/disre.c
Location:line 112, column 5
Description:Value stored to 'ip' 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/* 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 "typedefs.h"
46#include "types/commrec.h"
47#include "gromacs/utility/smalloc.h"
48#include "macros.h"
49#include "gromacs/math/vec.h"
50#include "gromacs/utility/futil.h"
51#include "gromacs/utility/fatalerror.h"
52#include "bondf.h"
53#include "copyrite.h"
54#include "disre.h"
55#include "main.h"
56#include "mtop_util.h"
57
58void init_disres(FILE *fplog, const gmx_mtop_t *mtop,
59 t_inputrec *ir, const t_commrec *cr,
60 t_fcdata *fcd, t_state *state, gmx_bool bIsREMD)
61{
62 int fa, nmol, i, npair, np;
63 t_iparams *ip;
64 t_disresdata *dd;
65 history_t *hist;
66 gmx_mtop_ilistloop_t iloop;
67 t_ilist *il;
68 char *ptr;
69
70 dd = &(fcd->disres);
71
72 if (gmx_mtop_ftype_count(mtop, F_DISRES) == 0)
73 {
74 dd->nres = 0;
75
76 return;
77 }
78
79 if (fplog)
80 {
81 fprintf(fplog, "Initializing the distance restraints\n");
82 }
83
84
85 if (ir->eDisre == edrEnsemble)
86 {
87 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/disre.c",
87, "Sorry, distance restraints with ensemble averaging over multiple molecules in one system are not functional in this version of GROMACS");
88 }
89
90 dd->dr_weighting = ir->eDisreWeighting;
91 dd->dr_fc = ir->dr_fc;
92 if (EI_DYNAMICS(ir->eI)(((ir->eI) == eiMD || ((ir->eI) == eiVV || (ir->eI) ==
eiVVAK)) || ((ir->eI) == eiSD1 || (ir->eI) == eiSD2) ||
(ir->eI) == eiBD))
93 {
94 dd->dr_tau = ir->dr_tau;
95 }
96 else
97 {
98 dd->dr_tau = 0.0;
99 }
100 if (dd->dr_tau == 0.0)
101 {
102 dd->dr_bMixed = FALSE0;
103 dd->ETerm = 0.0;
104 }
105 else
106 {
107 dd->dr_bMixed = ir->bDisreMixed;
108 dd->ETerm = exp(-(ir->delta_t/ir->dr_tau));
109 }
110 dd->ETerm1 = 1.0 - dd->ETerm;
111
112 ip = mtop->ffparams.iparams;
Value stored to 'ip' is never read
113
114 dd->nres = 0;
115 dd->npair = 0;
116 iloop = gmx_mtop_ilistloop_init(mtop);
117 while (gmx_mtop_ilistloop_next(iloop, &il, &nmol))
118 {
119 np = 0;
120 for (fa = 0; fa < il[F_DISRES].nr; fa += 3)
121 {
122 np++;
123 npair = mtop->ffparams.iparams[il[F_DISRES].iatoms[fa]].disres.npair;
124 if (np == npair)
125 {
126 dd->nres += (ir->eDisre == edrEnsemble ? 1 : nmol)*npair;
127 dd->npair += nmol*npair;
128 np = 0;
129 }
130 }
131 }
132
133 if (cr && PAR(cr)((cr)->nnodes > 1))
134 {
135 /* Temporary check, will be removed when disre is implemented with DD */
136 const char *notestr = "NOTE: atoms involved in distance restraints should be within the same domain. If this is not the case mdrun generates a fatal error. If you encounter this, use a single MPI rank (Verlet+OpenMP+GPUs work fine).";
137
138 if (MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1)))
139 {
140 fprintf(stderrstderr, "\n%s\n\n", notestr);
141 }
142 if (fplog)
143 {
144 fprintf(fplog, "%s\n", notestr);
145 }
146
147 if (dd->dr_tau != 0 || ir->eDisre == edrEnsemble || cr->ms != NULL((void*)0) ||
148 dd->nres != dd->npair)
149 {
150 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/disre.c",
150, "Time or ensemble averaged or multiple pair distance restraints do not work (yet) with domain decomposition, use a single MPI rank%s", cr->ms ? " per simulation" : "");
151 }
152 if (ir->nstdisreout != 0)
153 {
154 if (fplog)
155 {
156 fprintf(fplog, "\nWARNING: Can not write distance restraint data to energy file with domain decomposition\n\n");
157 }
158 if (MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1)))
159 {
160 fprintf(stderrstderr, "\nWARNING: Can not write distance restraint data to energy file with domain decomposition\n");
161 }
162 ir->nstdisreout = 0;
163 }
164 }
165
166 snew(dd->rt, dd->npair)(dd->rt) = save_calloc("dd->rt", "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/disre.c"
, 166, (dd->npair), sizeof(*(dd->rt)));
167
168 if (dd->dr_tau != 0.0)
169 {
170 hist = &state->hist;
171 /* Set the "history lack" factor to 1 */
172 state->flags |= (1<<estDISRE_INITF);
173 hist->disre_initf = 1.0;
174 /* Allocate space for the r^-3 time averages */
175 state->flags |= (1<<estDISRE_RM3TAV);
176 hist->ndisrepairs = dd->npair;
177 snew(hist->disre_rm3tav, hist->ndisrepairs)(hist->disre_rm3tav) = save_calloc("hist->disre_rm3tav"
, "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/disre.c", 177
, (hist->ndisrepairs), sizeof(*(hist->disre_rm3tav)));
178 }
179 /* Allocate space for a copy of rm3tav,
180 * so we can call do_force without modifying the state.
181 */
182 snew(dd->rm3tav, dd->npair)(dd->rm3tav) = save_calloc("dd->rm3tav", "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/disre.c"
, 182, (dd->npair), sizeof(*(dd->rm3tav)));
183
184 /* Allocate Rt_6 and Rtav_6 consecutively in memory so they can be
185 * averaged over the processors in one call (in calc_disre_R_6)
186 */
187 snew(dd->Rt_6, 2*dd->nres)(dd->Rt_6) = save_calloc("dd->Rt_6", "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/disre.c"
, 187, (2*dd->nres), sizeof(*(dd->Rt_6)));
188 dd->Rtav_6 = &(dd->Rt_6[dd->nres]);
189
190 ptr = getenv("GMX_DISRE_ENSEMBLE_SIZE");
191 if (cr && cr->ms != NULL((void*)0) && ptr != NULL((void*)0) && !bIsREMD)
192 {
193#ifdef GMX_MPI
194 dd->nsystems = 0;
195 sscanf(ptr, "%d", &dd->nsystems);
196 if (fplog)
197 {
198 fprintf(fplog, "Found GMX_DISRE_ENSEMBLE_SIZE set to %d systems per ensemble\n", dd->nsystems);
199 }
200 /* This check is only valid on MASTER(cr), so probably
201 * ensemble-averaged distance restraints are broken on more
202 * than one processor per simulation system. */
203 if (MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1)))
204 {
205 check_multi_int(fplog, cr->ms, dd->nsystems,
206 "the number of systems per ensemble",
207 FALSE0);
208 }
209 gmx_bcast_sim(sizeof(int), &dd->nsystems, cr);
210
211 /* We use to allow any value of nsystems which was a divisor
212 * of ms->nsim. But this required an extra communicator which
213 * was stored in t_fcdata. This pulled in mpi.h in nearly all C files.
214 */
215 if (!(cr->ms->nsim == 1 || cr->ms->nsim == dd->nsystems))
216 {
217 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/disre.c",
217, "GMX_DISRE_ENSEMBLE_SIZE (%d) is not equal to 1 or the number of systems (option -multi) %d", dd->nsystems, cr->ms->nsim);
218 }
219 if (fplog)
220 {
221 fprintf(fplog, "Our ensemble consists of systems:");
222 for (i = 0; i < dd->nsystems; i++)
223 {
224 fprintf(fplog, " %d",
225 (cr->ms->sim/dd->nsystems)*dd->nsystems+i);
226 }
227 fprintf(fplog, "\n");
228 }
229 snew(dd->Rtl_6, dd->nres)(dd->Rtl_6) = save_calloc("dd->Rtl_6", "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/disre.c"
, 229, (dd->nres), sizeof(*(dd->Rtl_6)));
230#endif
231 }
232 else
233 {
234 dd->nsystems = 1;
235 dd->Rtl_6 = dd->Rt_6;
236 }
237
238 if (dd->npair > 0)
239 {
240 if (fplog)
241 {
242 fprintf(fplog, "There are %d distance restraints involving %d atom pairs\n", dd->nres, dd->npair);
243 }
244 /* Have to avoid g_disre de-referencing cr blindly, mdrun not
245 * doing consistency checks for ensemble-averaged distance
246 * restraints when that's not happening, and only doing those
247 * checks from appropriate processes (since check_multi_int is
248 * too broken to check whether the communication will
249 * succeed...) */
250 if (cr && cr->ms && dd->nsystems > 1 && MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1)))
251 {
252 check_multi_int(fplog, cr->ms, fcd->disres.nres,
253 "the number of distance restraints",
254 FALSE0);
255 }
256 please_cite(fplog, "Tropp80a");
257 please_cite(fplog, "Torda89a");
258 }
259}
260
261void calc_disres_R_6(int nfa, const t_iatom forceatoms[], const t_iparams ip[],
262 const rvec x[], const t_pbc *pbc,
263 t_fcdata *fcd, history_t *hist)
264{
265 atom_id ai, aj;
266 int fa, res, i, pair, ki, kj, m;
267 int type, npair, np;
268 rvec dx;
269 real *rt, *rm3tav, *Rtl_6, *Rt_6, *Rtav_6;
270 real rt_1, rt_3, rt2;
271 ivec it, jt, dt;
272 t_disresdata *dd;
273 real ETerm, ETerm1, cf1 = 0, cf2 = 0, invn = 0;
274 gmx_bool bTav;
275
276 dd = &(fcd->disres);
277 bTav = (dd->dr_tau != 0);
278 ETerm = dd->ETerm;
279 ETerm1 = dd->ETerm1;
280 rt = dd->rt;
281 rm3tav = dd->rm3tav;
282 Rtl_6 = dd->Rtl_6;
283 Rt_6 = dd->Rt_6;
284 Rtav_6 = dd->Rtav_6;
285
286 if (bTav)
287 {
288 /* scaling factor to smoothly turn on the restraint forces *
289 * when using time averaging */
290 dd->exp_min_t_tau = hist->disre_initf*ETerm;
291
292 cf1 = dd->exp_min_t_tau;
293 cf2 = 1.0/(1.0 - dd->exp_min_t_tau);
294 }
295
296 if (dd->nsystems > 1)
297 {
298 invn = 1.0/dd->nsystems;
299 }
300
301 /* 'loop' over all atom pairs (pair_nr=fa/3) involved in restraints, *
302 * the total number of atoms pairs is nfa/3 */
303 res = 0;
304 fa = 0;
305 while (fa < nfa)
306 {
307 type = forceatoms[fa];
308 npair = ip[type].disres.npair;
309
310 Rtav_6[res] = 0.0;
311 Rt_6[res] = 0.0;
312
313 /* Loop over the atom pairs of 'this' restraint */
314 np = 0;
315 while (fa < nfa && np < npair)
316 {
317 pair = fa/3;
318 ai = forceatoms[fa+1];
319 aj = forceatoms[fa+2];
320
321 if (pbc)
322 {
323 pbc_dx_aiuc(pbc, x[ai], x[aj], dx);
324 }
325 else
326 {
327 rvec_sub(x[ai], x[aj], dx);
328 }
329 rt2 = iprod(dx, dx);
330 rt_1 = gmx_invsqrt(rt2)gmx_software_invsqrt(rt2);
331 rt_3 = rt_1*rt_1*rt_1;
332
333 rt[pair] = sqrt(rt2);
334 if (bTav)
335 {
336 /* Here we update rm3tav in t_fcdata using the data
337 * in history_t.
338 * Thus the results stay correct when this routine
339 * is called multiple times.
340 */
341 rm3tav[pair] = cf2*((ETerm - cf1)*hist->disre_rm3tav[pair] +
342 ETerm1*rt_3);
343 }
344 else
345 {
346 rm3tav[pair] = rt_3;
347 }
348
349 Rt_6[res] += rt_3*rt_3;
350 Rtav_6[res] += rm3tav[pair]*rm3tav[pair];
351
352 fa += 3;
353 np++;
354 }
355 if (dd->nsystems > 1)
356 {
357 Rtl_6[res] = Rt_6[res];
358 Rt_6[res] *= invn;
359 Rtav_6[res] *= invn;
360 }
361
362 res++;
363 }
364}
365
366real ta_disres(int nfa, const t_iatom forceatoms[], const t_iparams ip[],
367 const rvec x[], rvec f[], rvec fshift[],
368 const t_pbc *pbc, const t_graph *g,
369 real gmx_unused__attribute__ ((unused)) lambda, real gmx_unused__attribute__ ((unused)) *dvdlambda,
370 const t_mdatoms gmx_unused__attribute__ ((unused)) *md, t_fcdata *fcd,
371 int gmx_unused__attribute__ ((unused)) *global_atom_index)
372{
373 const real sixth = 1.0/6.0;
374 const real seven_three = 7.0/3.0;
375
376 atom_id ai, aj;
377 int fa, res, npair, p, pair, ki = CENTRAL(((2*1 +1)*(2*1 +1)*(2*2 +1))/2), m;
378 int type;
379 rvec dx;
380 real weight_rt_1;
381 real smooth_fc, Rt, Rtav, rt2, *Rtl_6, *Rt_6, *Rtav_6;
382 real k0, f_scal = 0, fmax_scal, fk_scal, fij;
383 real tav_viol, instant_viol, mixed_viol, violtot, vtot;
384 real tav_viol_Rtav7, instant_viol_Rtav7;
385 real up1, up2, low;
386 gmx_bool bConservative, bMixed, bViolation;
387 ivec it, jt, dt;
388 t_disresdata *dd;
389 int dr_weighting;
390 gmx_bool dr_bMixed;
391
392 dd = &(fcd->disres);
393 dr_weighting = dd->dr_weighting;
394 dr_bMixed = dd->dr_bMixed;
395 Rtl_6 = dd->Rtl_6;
396 Rt_6 = dd->Rt_6;
397 Rtav_6 = dd->Rtav_6;
398
399 tav_viol = instant_viol = mixed_viol = tav_viol_Rtav7 = instant_viol_Rtav7 = 0;
400
401 smooth_fc = dd->dr_fc;
402 if (dd->dr_tau != 0)
403 {
404 /* scaling factor to smoothly turn on the restraint forces *
405 * when using time averaging */
406 smooth_fc *= (1.0 - dd->exp_min_t_tau);
407 }
408
409 violtot = 0;
410 vtot = 0;
411
412 /* 'loop' over all atom pairs (pair_nr=fa/3) involved in restraints, *
413 * the total number of atoms pairs is nfa/3 */
414 res = 0;
415 fa = 0;
416 while (fa < nfa)
417 {
418 type = forceatoms[fa];
419 /* Take action depending on restraint, calculate scalar force */
420 npair = ip[type].disres.npair;
421 up1 = ip[type].disres.up1;
422 up2 = ip[type].disres.up2;
423 low = ip[type].disres.low;
424 k0 = smooth_fc*ip[type].disres.kfac;
425
426 /* save some flops when there is only one pair */
427 if (ip[type].disres.type != 2)
428 {
429 bConservative = (dr_weighting == edrwConservative) && (npair > 1);
430 bMixed = dr_bMixed;
431 Rt = pow(Rt_6[res], -sixth);
432 Rtav = pow(Rtav_6[res], -sixth);
433 }
434 else
435 {
436 /* When rtype=2 use instantaneous not ensemble avereged distance */
437 bConservative = (npair > 1);
438 bMixed = FALSE0;
439 Rt = pow(Rtl_6[res], -sixth);
440 Rtav = Rt;
441 }
442
443 if (Rtav > up1)
444 {
445 bViolation = TRUE1;
446 tav_viol = Rtav - up1;
447 }
448 else if (Rtav < low)
449 {
450 bViolation = TRUE1;
451 tav_viol = Rtav - low;
452 }
453 else
454 {
455 bViolation = FALSE0;
456 }
457
458 if (bViolation)
459 {
460 /* NOTE:
461 * there is no real potential when time averaging is applied
462 */
463 vtot += 0.5*k0*sqr(tav_viol);
464 if (1/vtot == 0)
465 {
466 printf("vtot is inf: %f\n", vtot);
467 }
468 if (!bMixed)
469 {
470 f_scal = -k0*tav_viol;
471 violtot += fabs(tav_viol);
472 }
473 else
474 {
475 if (Rt > up1)
476 {
477 if (tav_viol > 0)
478 {
479 instant_viol = Rt - up1;
480 }
481 else
482 {
483 bViolation = FALSE0;
484 }
485 }
486 else if (Rt < low)
487 {
488 if (tav_viol < 0)
489 {
490 instant_viol = Rt - low;
491 }
492 else
493 {
494 bViolation = FALSE0;
495 }
496 }
497 else
498 {
499 bViolation = FALSE0;
500 }
501 if (bViolation)
502 {
503 mixed_viol = sqrt(tav_viol*instant_viol);
504 f_scal = -k0*mixed_viol;
505 violtot += mixed_viol;
506 }
507 }
508 }
509
510 if (bViolation)
511 {
512 fmax_scal = -k0*(up2-up1);
513 /* Correct the force for the number of restraints */
514 if (bConservative)
515 {
516 f_scal = max(f_scal, fmax_scal)(((f_scal) > (fmax_scal)) ? (f_scal) : (fmax_scal) );
517 if (!bMixed)
518 {
519 f_scal *= Rtav/Rtav_6[res];
520 }
521 else
522 {
523 f_scal /= 2*mixed_viol;
524 tav_viol_Rtav7 = tav_viol*Rtav/Rtav_6[res];
525 instant_viol_Rtav7 = instant_viol*Rt/Rt_6[res];
526 }
527 }
528 else
529 {
530 f_scal /= (real)npair;
531 f_scal = max(f_scal, fmax_scal)(((f_scal) > (fmax_scal)) ? (f_scal) : (fmax_scal) );
532 }
533
534 /* Exert the force ... */
535
536 /* Loop over the atom pairs of 'this' restraint */
537 for (p = 0; p < npair; p++)
538 {
539 pair = fa/3;
540 ai = forceatoms[fa+1];
541 aj = forceatoms[fa+2];
542
543 if (pbc)
544 {
545 ki = pbc_dx_aiuc(pbc, x[ai], x[aj], dx);
546 }
547 else
548 {
549 rvec_sub(x[ai], x[aj], dx);
550 }
551 rt2 = iprod(dx, dx);
552
553 weight_rt_1 = gmx_invsqrt(rt2)gmx_software_invsqrt(rt2);
554
555 if (bConservative)
556 {
557 if (!dr_bMixed)
558 {
559 weight_rt_1 *= pow(dd->rm3tav[pair], seven_three);
560 }
561 else
562 {
563 weight_rt_1 *= tav_viol_Rtav7*pow(dd->rm3tav[pair], seven_three)+
564 instant_viol_Rtav7*pow(dd->rt[pair], -7);
565 }
566 }
567
568 fk_scal = f_scal*weight_rt_1;
569
570 if (g)
571 {
572 ivec_sub(SHIFT_IVEC(g, ai)((g)->ishift[ai]), SHIFT_IVEC(g, aj)((g)->ishift[aj]), dt);
573 ki = IVEC2IS(dt)(((2*2 +1)*((2*1 +1)*(((dt)[2])+1)+((dt)[1])+1)+((dt)[0])+2));
574 }
575
576 for (m = 0; m < DIM3; m++)
577 {
578 fij = fk_scal*dx[m];
579
580 f[ai][m] += fij;
581 f[aj][m] -= fij;
582 fshift[ki][m] += fij;
583 fshift[CENTRAL(((2*1 +1)*(2*1 +1)*(2*2 +1))/2)][m] -= fij;
584 }
585 fa += 3;
586 }
587 }
588 else
589 {
590 /* No violation so force and potential contributions */
591 fa += 3*npair;
592 }
593 res++;
594 }
595
596 dd->sumviol = violtot;
597
598 /* Return energy */
599 return vtot;
600}
601
602void update_disres_history(t_fcdata *fcd, history_t *hist)
603{
604 t_disresdata *dd;
605 int pair;
606
607 dd = &(fcd->disres);
608 if (dd->dr_tau != 0)
609 {
610 /* Copy the new time averages that have been calculated
611 * in calc_disres_R_6.
612 */
613 hist->disre_initf = dd->exp_min_t_tau;
614 for (pair = 0; pair < dd->npair; pair++)
615 {
616 hist->disre_rm3tav[pair] = dd->rm3tav[pair];
617 }
618 }
619}