Bug Summary

File:gromacs/mdlib/wall.c
Location:line 317, column 43
Description:The left operand of '*' is a garbage value

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-2008, The GROMACS development team.
6 * Copyright (c) 2013,2014, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
10 *
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
15 *
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
20 *
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
25 *
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
33 *
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
36 */
37
38#ifdef HAVE_CONFIG_H1
39#include <config.h>
40#endif
41
42#include <math.h>
43#include <string.h>
44
45#include "typedefs.h"
46#include "macros.h"
47#include "force.h"
48#include "nrnb.h"
49#include "gromacs/math/vec.h"
50
51#include "gromacs/fileio/filenm.h"
52#include "gromacs/math/utilities.h"
53#include "gromacs/utility/cstringutil.h"
54#include "gromacs/utility/smalloc.h"
55
56void make_wall_tables(FILE *fplog, const output_env_t oenv,
57 const t_inputrec *ir, const char *tabfn,
58 const gmx_groups_t *groups,
59 t_forcerec *fr)
60{
61 int w, negp_pp, egp, i, j;
62 int *nm_ind;
63 char buf[STRLEN4096];
64 t_forcetable *tab;
65
66 negp_pp = ir->opts.ngener - ir->nwall;
67 nm_ind = groups->grps[egcENER].nm_ind;
68
69 if (fplog)
70 {
71 fprintf(fplog, "Reading user tables for %d energy groups with %d walls\n",
72 negp_pp, ir->nwall);
73 }
74
75 snew(fr->wall_tab, ir->nwall)(fr->wall_tab) = save_calloc("fr->wall_tab", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/wall.c"
, 75, (ir->nwall), sizeof(*(fr->wall_tab)));
76 for (w = 0; w < ir->nwall; w++)
77 {
78 snew(fr->wall_tab[w], negp_pp)(fr->wall_tab[w]) = save_calloc("fr->wall_tab[w]", "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/wall.c"
, 78, (negp_pp), sizeof(*(fr->wall_tab[w])));
79 for (egp = 0; egp < negp_pp; egp++)
80 {
81 /* If the energy group pair is excluded, we don't need a table */
82 if (!(fr->egp_flags[egp*ir->opts.ngener+negp_pp+w] & EGP_EXCL(1<<0)))
83 {
84 tab = &fr->wall_tab[w][egp];
85 sprintf(buf, "%s", tabfn);
86 sprintf(buf + strlen(tabfn) - strlen(ftp2ext(efXVG)) - 1, "_%s_%s.%s",
87 *groups->grpname[nm_ind[egp]],
88 *groups->grpname[nm_ind[negp_pp+w]],
89 ftp2ext(efXVG));
90 *tab = make_tables(fplog, oenv, fr, FALSE0, buf, 0, GMX_MAKETABLES_FORCEUSER(1<<0));
91 /* Since wall have no charge, we can compress the table */
92 for (i = 0; i <= tab->n; i++)
93 {
94 for (j = 0; j < 8; j++)
95 {
96 tab->data[8*i+j] = tab->data[12*i+4+j];
97 }
98 }
99 }
100 }
101 }
102}
103
104static void wall_error(int a, rvec *x, real r)
105{
106 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/wall.c", 106,
107 "An atom is beyond the wall: coordinates %f %f %f, distance %f\n"
108 "You might want to use the mdp option wall_r_linpot",
109 x[a][XX0], x[a][YY1], x[a][ZZ2], r);
110}
111
112real do_walls(t_inputrec *ir, t_forcerec *fr, matrix box, t_mdatoms *md,
113 rvec x[], rvec f[], real lambda, real Vlj[], t_nrnb *nrnb)
114{
115 int nwall, w, lam, i;
116 int ntw[2], at, ntype, ngid, ggid, *egp_flags, *type;
117 real *nbfp, lamfac, fac_d[2], fac_r[2], Cd, Cr, Vtot, Fwall[2];
118 real wall_z[2], r, mr, r1, r2, r4, Vd, Vr, V = 0, Fd, Fr, F = 0, dvdlambda;
119 dvec xf_z;
120 int n0, nnn;
121 real tabscale, *VFtab, rt, eps, eps2, Yt, Ft, Geps, Heps, Heps2, Fp, VV, FF;
122 unsigned short *gid = md->cENER;
123 t_forcetable *tab;
124
125 nwall = ir->nwall;
126 ngid = ir->opts.ngener;
127 ntype = fr->ntype;
128 nbfp = fr->nbfp;
129 egp_flags = fr->egp_flags;
130
131 for (w = 0; w < nwall; w++)
1
Assuming 'w' is < 'nwall'
→
2
←
Loop condition is true. Entering loop body
→
5
←
Assuming 'w' is < 'nwall'
→
6
←
Loop condition is true. Entering loop body
→
9
←
Assuming 'w' is < 'nwall'
→
10
←
Loop condition is true. Entering loop body
→
13
←
Assuming 'w' is >= 'nwall'
→
14
←
Loop condition is false. Execution continues on line 149
→
132 {
133 ntw[w] = 2*ntype*ir->wall_atomtype[w];
134 switch (ir->wall_type)
3
←
Control jumps to the 'default' case at line 144
→
7
←
Control jumps to the 'default' case at line 144
→
11
←
Control jumps to the 'default' case at line 144
→
135 {
136 case ewt93:
137 fac_d[w] = ir->wall_density[w]*M_PI3.14159265358979323846/6;
138 fac_r[w] = ir->wall_density[w]*M_PI3.14159265358979323846/45;
139 break;
140 case ewt104:
141 fac_d[w] = ir->wall_density[w]*M_PI3.14159265358979323846/2;
142 fac_r[w] = ir->wall_density[w]*M_PI3.14159265358979323846/5;
143 break;
144 default:
145 break;
4
←
 Execution continues on line 147
→
8
←
 Execution continues on line 147
→
12
←
 Execution continues on line 147
→
146 }
147 Fwall[w] = 0;
148 }
149 wall_z[0] = 0;
150 wall_z[1] = box[ZZ2][ZZ2];
151
152 Vtot = 0;
153 dvdlambda = 0;
154 clear_dvec(xf_z);
155 for (lam = 0; lam < (md->nPerturbed ? 2 : 1); lam++)
15
←
'?' condition is false
→
16
←
Loop condition is true. Entering loop body
→
156 {
157 if (md->nPerturbed)
17
←
Taking false branch
→
158 {
159 if (lam == 0)
160 {
161 lamfac = 1 - lambda;
162 type = md->typeA;
163 }
164 else
165 {
166 lamfac = 0;
167 type = md->typeB;
168 }
169 }
170 else
171 {
172 lamfac = 1;
173 type = md->typeA;
174 }
175 for (i = 0; i < md->homenr; i++)
18
←
Loop condition is true. Entering loop body
→
176 {
177 for (w = 0; w < nwall; w++)
19
←
Loop condition is true. Entering loop body
→
21
←
Loop condition is true. Entering loop body
→
23
←
Loop condition is true. Entering loop body
→
178 {
179 /* The wall energy groups are always at the end of the list */
180 ggid = gid[i]*ngid + ngid - nwall + w;
181 at = type[i];
182 /* nbfp now includes the 6.0/12.0 derivative prefactors */
183 Cd = nbfp[ntw[w]+2*at]/6.0;
184 Cr = nbfp[ntw[w]+2*at+1]/12.0;
185 if (!((Cd == 0 && Cr == 0) || (egp_flags[ggid] & EGP_EXCL(1<<0))))
20
←
Taking false branch
→
22
←
Taking false branch
→
24
←
Taking true branch
→
186 {
187 if (w == 0)
25
←
Taking false branch
→
188 {
189 r = x[i][ZZ2];
190 }
191 else
192 {
193 r = wall_z[1] - x[i][ZZ2];
194 }
195 if (r < ir->wall_r_linpot)
26
←
Taking false branch
→
196 {
197 mr = ir->wall_r_linpot - r;
198 r = ir->wall_r_linpot;
199 }
200 else
201 {
202 mr = 0;
203 }
204 switch (ir->wall_type)
27
←
Control jumps to the 'default' case at line 292
→
205 {
206 case ewtTABLE:
207 if (r < 0)
208 {
209 wall_error(i, x, r);
210 }
211 tab = &(fr->wall_tab[w][gid[i]]);
212 tabscale = tab->scale;
213 VFtab = tab->data;
214
215 rt = r*tabscale;
216 n0 = rt;
217 if (n0 >= tab->n)
218 {
219 /* Beyond the table range, set V and F to zero */
220 V = 0;
221 F = 0;
222 }
223 else
224 {
225 eps = rt - n0;
226 eps2 = eps*eps;
227 /* Dispersion */
228 nnn = 8*n0;
229 Yt = VFtab[nnn];
230 Ft = VFtab[nnn+1];
231 Geps = VFtab[nnn+2]*eps;
232 Heps2 = VFtab[nnn+3]*eps2;
233 Fp = Ft + Geps + Heps2;
234 VV = Yt + Fp*eps;
235 FF = Fp + Geps + 2.0*Heps2;
236 Vd = Cd*VV;
237 Fd = Cd*FF;
238 /* Repulsion */
239 nnn = nnn + 4;
240 Yt = VFtab[nnn];
241 Ft = VFtab[nnn+1];
242 Geps = VFtab[nnn+2]*eps;
243 Heps2 = VFtab[nnn+3]*eps2;
244 Fp = Ft + Geps + Heps2;
245 VV = Yt + Fp*eps;
246 FF = Fp + Geps + 2.0*Heps2;
247 Vr = Cr*VV;
248 Fr = Cr*FF;
249 V = Vd + Vr;
250 F = -lamfac*(Fd + Fr)*tabscale;
251 }
252 break;
253 case ewt93:
254 if (r <= 0)
255 {
256 wall_error(i, x, r);
257 }
258 r1 = 1/r;
259 r2 = r1*r1;
260 r4 = r2*r2;
261 Vd = fac_d[w]*Cd*r2*r1;
262 Vr = fac_r[w]*Cr*r4*r4*r1;
263 V = Vr - Vd;
264 F = lamfac*(9*Vr - 3*Vd)*r1;
265 break;
266 case ewt104:
267 if (r <= 0)
268 {
269 wall_error(i, x, r);
270 }
271 r1 = 1/r;
272 r2 = r1*r1;
273 r4 = r2*r2;
274 Vd = fac_d[w]*Cd*r4;
275 Vr = fac_r[w]*Cr*r4*r4*r2;
276 V = Vr - Vd;
277 F = lamfac*(10*Vr - 4*Vd)*r1;
278 break;
279 case ewt126:
280 if (r <= 0)
281 {
282 wall_error(i, x, r);
283 }
284 r1 = 1/r;
285 r2 = r1*r1;
286 r4 = r2*r2;
287 Vd = Cd*r4*r2;
288 Vr = Cr*r4*r4*r4;
289 V = Vr - Vd;
290 F = lamfac*(12*Vr - 6*Vd)*r1;
291 break;
292 default:
293 break;
28
←
 Execution continues on line 295
→
294 }
295 if (mr > 0)
29
←
Taking false branch
→
296 {
297 V += mr*F;
298 }
299 if (w == 1)
30
←
Taking false branch
→
300 {
301 F = -F;
302 }
303 Vlj[ggid] += lamfac*V;
304 Vtot += V;
305 f[i][ZZ2] += F;
306 /* Because of the single sum virial calculation we need
307 * to add the full virial contribution of the walls.
308 * Since the force only has a z-component, there is only
309 * a contribution to the z component of the virial tensor.
310 * We could also determine the virial contribution directly,
311 * which would be cheaper here, but that would require extra
312 * communication for f_novirsum for with virtual sites
313 * in parallel.
314 */
315 xf_z[XX0] -= x[i][XX0]*F;
316 xf_z[YY1] -= x[i][YY1]*F;
317 xf_z[ZZ2] -= wall_z[w]*F;
31
←
The left operand of '*' is a garbage value
318 }
319 }
320 }
321 if (md->nPerturbed)
322 {
323 dvdlambda += (lam == 0 ? -1 : 1)*Vtot;
324 }
325
326 inc_nrnb(nrnb, eNR_WALLS, md->homenr)(nrnb)->n[eNR_WALLS] += md->homenr;
327 }
328
329 for (i = 0; i < DIM3; i++)
330 {
331 fr->vir_wall_z[i] = -0.5*xf_z[i];
332 }
333
334 return dvdlambda;
335}