Bug Summary

File:gromacs/gmxlib/restcbt.c
Location:line 120, column 5
Description:Value stored to 'sine_phi0' 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) 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#include <math.h>
39#include <assert.h>
40#include "physics.h"
41#include "gromacs/math/vec.h"
42#include "gromacs/math/utilities.h"
43#include "txtdump.h"
44#include "bondf.h"
45#include "gromacs/utility/smalloc.h"
46#include "pbc.h"
47#include "ns.h"
48#include "macros.h"
49#include "names.h"
50#include "mshift.h"
51#include "disre.h"
52#include "orires.h"
53#include "force.h"
54#include "nonbonded.h"
55
56/* This function computes factors needed for restricted angle potential.
57 * For explanations on formula used see file "restcbt.h" */
58
59void compute_factors_restangles(int type, const t_iparams forceparams[],
60 rvec delta_ante, rvec delta_post,
61 real *prefactor, real *ratio_ante, real *ratio_post, real *v)
62{
63 real theta_equil, k_bending;
64 real cosine_theta_equil;
65 real c_ante, c_cros, c_post;
66 real norm;
67 real delta_cosine, cosine_theta;
68 real sine_theta_sq;
69 real term_theta_theta_equil;
70
71 k_bending = forceparams[type].harmonic.krA;
72 theta_equil = forceparams[type].harmonic.rA*DEG2RAD(3.14159265358979323846/180.0);
73 theta_equil = M_PI3.14159265358979323846 - theta_equil;
74 cosine_theta_equil = cos(theta_equil);
75
76 c_ante = iprod(delta_ante, delta_ante);
77 c_cros = iprod(delta_ante, delta_post);
78 c_post = iprod(delta_post, delta_post);
79
80 norm = gmx_invsqrt(c_ante * c_post)gmx_software_invsqrt(c_ante * c_post);
81 cosine_theta = c_cros * norm;
82 sine_theta_sq = 1 - cosine_theta * cosine_theta;
83
84 *ratio_ante = c_cros / c_ante;
85 *ratio_post = c_cros / c_post;
86
87 delta_cosine = cosine_theta - cosine_theta_equil;
88 term_theta_theta_equil = 1 - cosine_theta * cosine_theta_equil;
89 *prefactor = -(k_bending) * delta_cosine * norm * term_theta_theta_equil / (sine_theta_sq * sine_theta_sq);
90
91 *v = k_bending * 0.5 * delta_cosine * delta_cosine / sine_theta_sq;
92
93}
94
95
96/* Compute factors for restricted dihedral potential
97 * For explanations on formula used see file "restcbt.h" */
98void compute_factors_restrdihs(int type, const t_iparams forceparams[],
99 rvec delta_ante, rvec delta_crnt, rvec delta_post,
100 real *factor_phi_ai_ante, real *factor_phi_ai_crnt, real *factor_phi_ai_post,
101 real *factor_phi_aj_ante, real *factor_phi_aj_crnt, real *factor_phi_aj_post,
102 real *factor_phi_ak_ante, real *factor_phi_ak_crnt, real *factor_phi_ak_post,
103 real *factor_phi_al_ante, real *factor_phi_al_crnt, real *factor_phi_al_post,
104 real *prefactor_phi, real *v)
105{
106
107 real phi0, sine_phi0, cosine_phi0;
108 real k_torsion;
109 real c_self_ante, c_self_crnt, c_self_post;
110 real c_cros_ante, c_cros_acrs, c_cros_post;
111 real c_prod, d_post, d_ante;
112 real sine_phi_sq, cosine_phi;
113 real delta_cosine, term_phi_phi0;
114 real ratio_phi_ante, ratio_phi_post;
115 real cos_phi, norm_phi;
116
117 /* Read parameters phi0 and k_torsion */
118 phi0 = forceparams[type].pdihs.phiA * DEG2RAD(3.14159265358979323846/180.0);
119 cosine_phi0 = cos(phi0);
120 sine_phi0 = sin(phi0);
Value stored to 'sine_phi0' is never read
121 k_torsion = forceparams[type].pdihs.cpA;
122
123 /* Computation of the cosine of the dihedral angle. The scalar ("dot") product method
124 * is used. c_*_* cummulate the scalar products of the differences of particles
125 * positions while c_prod, d_ante and d_post are differences of products of scalar
126 * terms that are parts of the derivatives of forces */
127 c_self_ante = iprod(delta_ante, delta_ante);
128 c_self_crnt = iprod(delta_crnt, delta_crnt);
129 c_self_post = iprod(delta_post, delta_post);
130 c_cros_ante = iprod(delta_ante, delta_crnt);
131 c_cros_acrs = iprod(delta_ante, delta_post);
132 c_cros_post = iprod(delta_crnt, delta_post);
133 c_prod = c_cros_ante * c_cros_post - c_self_crnt * c_cros_acrs;
134 d_ante = c_self_ante * c_self_crnt - c_cros_ante * c_cros_ante;
135 d_post = c_self_post * c_self_crnt - c_cros_post * c_cros_post;
136
137 /* When three consecutive beads align, we obtain values close to zero.
138 * Here we avoid small values to prevent round-off errors. */
139 if (d_ante < GMX_REAL_EPS5.96046448E-08)
140 {
141 d_ante = GMX_REAL_EPS5.96046448E-08;
142 }
143 if (d_post < GMX_REAL_EPS5.96046448E-08)
144 {
145 d_post = GMX_REAL_EPS5.96046448E-08;
146 }
147
148 /* Computes the square of the sinus of phi in sine_phi_sq */
149 norm_phi = gmx_invsqrt(d_ante * d_post)gmx_software_invsqrt(d_ante * d_post);
150 cosine_phi = c_prod * norm_phi;
151 sine_phi_sq = 1.0 - cosine_phi * cosine_phi;
152
153 /* It is possible that cosine_phi is slightly bigger than 1.0 due to round-off errors. */
154 if (sine_phi_sq < 0.0)
155 {
156 sine_phi_sq = 0.0;
157 }
158
159 /* Computation of the differences of cosines (delta_cosine) and a term (term_phi_phi0)
160 * that is part of the common prefactor_phi */
161
162 delta_cosine = cosine_phi - cosine_phi0;
163 term_phi_phi0 = 1 - cosine_phi * cosine_phi0;
164
165
166 /* Computation of ratios */
167 ratio_phi_ante = c_prod / d_ante;
168 ratio_phi_post = c_prod / d_post;
169
170 /* Computation of the prefactor - common term for all forces */
171 *prefactor_phi = -(k_torsion) * delta_cosine * norm_phi * term_phi_phi0 / (sine_phi_sq * sine_phi_sq);
172
173 /* Computation of force factors. Factors factor_phi_* are coming from the
174 * derivatives of the torsion angle (phi) with respect to the beads ai, aj, al, ak,
175 * (four) coordinates and they are multiplied in the force computations with the
176 * differences of the particles positions stored in parameters delta_ante,
177 * delta_crnt, delta_post. For formulas see file "restcbt.h" */
178
179 *factor_phi_ai_ante = ratio_phi_ante * c_self_crnt;
180 *factor_phi_ai_crnt = -c_cros_post - ratio_phi_ante * c_cros_ante;
181 *factor_phi_ai_post = c_self_crnt;
182 *factor_phi_aj_ante = -c_cros_post - ratio_phi_ante * (c_self_crnt + c_cros_ante);
183 *factor_phi_aj_crnt = c_cros_post + c_cros_acrs * 2.0 + ratio_phi_ante * (c_self_ante + c_cros_ante) + ratio_phi_post * c_self_post;
184 *factor_phi_aj_post = -(c_cros_ante + c_self_crnt) - ratio_phi_post * c_cros_post;
185 *factor_phi_ak_ante = c_cros_post + c_self_crnt + ratio_phi_ante * c_cros_ante;
186 *factor_phi_ak_crnt = -(c_cros_ante + c_cros_acrs * 2.0)- ratio_phi_ante * c_self_ante - ratio_phi_post * (c_self_post + c_cros_post);
187 *factor_phi_ak_post = c_cros_ante + ratio_phi_post * (c_self_crnt + c_cros_post);
188 *factor_phi_al_ante = -c_self_crnt;
189 *factor_phi_al_crnt = c_cros_ante + ratio_phi_post * c_cros_post;
190 *factor_phi_al_post = -ratio_phi_post * c_self_crnt;
191
192 /* Contribution to energy - see formula in file "restcbt.h"*/
193 *v = k_torsion * 0.5 * delta_cosine * delta_cosine / sine_phi_sq;
194}
195
196
197
198/* Compute factors for CBT potential
199 * For explanations on formula used see file "restcbt.h" */
200
201void compute_factors_cbtdihs(int type, const t_iparams forceparams[],
202 rvec delta_ante, rvec delta_crnt, rvec delta_post,
203 rvec f_phi_ai, rvec f_phi_aj, rvec f_phi_ak, rvec f_phi_al,
204 rvec f_theta_ante_ai, rvec f_theta_ante_aj, rvec f_theta_ante_ak,
205 rvec f_theta_post_aj, rvec f_theta_post_ak, rvec f_theta_post_al,
206 real * v)
207{
208 int j, d;
209 real torsion_coef[NR_CBTDIHS6];
210 real c_self_ante, c_self_crnt, c_self_post;
211 real c_cros_ante, c_cros_acrs, c_cros_post;
212 real c_prod, d_ante, d_post;
213 real norm_phi, norm_theta_ante, norm_theta_post;
214 real cosine_phi, cosine_theta_ante, cosine_theta_post;
215 real sine_theta_ante_sq, sine_theta_post_sq;
216 real sine_theta_ante, sine_theta_post;
217 real prefactor_phi;
218 real ratio_phi_ante, ratio_phi_post;
219 real r1, r2;
220 real factor_phi_ai_ante, factor_phi_ai_crnt, factor_phi_ai_post;
221 real factor_phi_aj_ante, factor_phi_aj_crnt, factor_phi_aj_post;
222 real factor_phi_ak_ante, factor_phi_ak_crnt, factor_phi_ak_post;
223 real factor_phi_al_ante, factor_phi_al_crnt, factor_phi_al_post;
224 real prefactor_theta_ante, ratio_theta_ante_ante, ratio_theta_ante_crnt;
225 real prefactor_theta_post, ratio_theta_post_crnt, ratio_theta_post_post;
226
227 /* The formula for combined bending-torsion potential (see file "restcbt.h") contains
228 * in its expression not only the dihedral angle \f[\phi\f] but also \f[\theta_{i-1}\f]
229 * (theta_ante bellow) and \f[\theta_{i}\f] (theta_post bellow)--- the adjacent bending
230 * angles. The forces for the particles ai, aj, ak, al have components coming from the
231 * derivatives of the potential with respect to all three angles.
232 * This function is organised in 4 parts
233 * PART 1 - Computes force factors common to all the derivatives for the four particles
234 * PART 2 - Computes the force components due to the derivatives of dihedral angle Phi
235 * PART 3 - Computes the force components due to the derivatives of bending angle Theta_Ante
236 * PART 4 - Computes the force components due to the derivatives of bending angle Theta_Post
237 * Bellow we will respct thuis structure */
238
239
240 /* PART 1 - COMPUTES FORCE FACTORS COMMON TO ALL DERIVATIVES FOR THE FOUR PARTICLES */
241
242
243 for (j = 0; (j < NR_CBTDIHS6); j++)
244 {
245 torsion_coef[j] = forceparams[type].cbtdihs.cbtcA[j];
246 }
247
248 /* Computation of the cosine of the dihedral angle. The scalar ("dot") product method
249 * is used. c_*_* cummulate the scalar products of the differences of particles
250 * positions while c_prod, d_ante and d_post are differences of products of scalar
251 * terms that are parts of the derivatives of forces */
252
253 c_self_ante = iprod(delta_ante, delta_ante);
254 c_self_crnt = iprod(delta_crnt, delta_crnt);
255 c_self_post = iprod(delta_post, delta_post);
256 c_cros_ante = iprod(delta_ante, delta_crnt);
257 c_cros_acrs = iprod(delta_ante, delta_post);
258 c_cros_post = iprod(delta_crnt, delta_post);
259 c_prod = c_cros_ante * c_cros_post - c_self_crnt * c_cros_acrs;
260 d_ante = c_self_ante * c_self_crnt - c_cros_ante * c_cros_ante;
261 d_post = c_self_post * c_self_crnt - c_cros_post * c_cros_post;
262
263 /* When three consecutive beads align, we obtain values close to zero.
264 Here we avoid small values to prevent round-off errors. */
265 if (d_ante < GMX_REAL_EPS5.96046448E-08)
266 {
267 d_ante = GMX_REAL_EPS5.96046448E-08;
268 }
269 if (d_post < GMX_REAL_EPS5.96046448E-08)
270 {
271 d_post = GMX_REAL_EPS5.96046448E-08;
272 }
273
274 /* Computations of cosines */
275 norm_phi = gmx_invsqrt(d_ante * d_post)gmx_software_invsqrt(d_ante * d_post);
276 norm_theta_ante = gmx_invsqrt(c_self_ante * c_self_crnt)gmx_software_invsqrt(c_self_ante * c_self_crnt);
277 norm_theta_post = gmx_invsqrt(c_self_crnt * c_self_post)gmx_software_invsqrt(c_self_crnt * c_self_post);
278 cosine_phi = c_prod * norm_phi;
279 cosine_theta_ante = c_cros_ante * norm_theta_ante;
280 cosine_theta_post = c_cros_post * norm_theta_post;
281 sine_theta_ante_sq = 1 - cosine_theta_ante * cosine_theta_ante;
282 sine_theta_post_sq = 1 - cosine_theta_post * cosine_theta_post;
283
284 /* It is possible that cosine_theta is slightly bigger than 1.0 due to round-off errors. */
285 if (sine_theta_ante_sq < 0.0)
286 {
287 sine_theta_ante_sq = 0.0;
288 }
289 if (sine_theta_post_sq < 0.0)
290 {
291 sine_theta_post_sq = 0.0;
292 }
293
294 sine_theta_ante = sqrt(sine_theta_ante_sq);
295 sine_theta_post = sqrt(sine_theta_post_sq);
296
297 /* PART 2 - COMPUTES FORCE COMPONENTS DUE TO DERIVATIVES TO DIHEDRAL ANGLE PHI */
298
299 /* Computation of ratios */
300 ratio_phi_ante = c_prod / d_ante;
301 ratio_phi_post = c_prod / d_post;
302
303 /* Computation of the prefactor */
304 /* Computing 2nd power */
305 r1 = cosine_phi;
306
307 prefactor_phi = -torsion_coef[0] * norm_phi * (torsion_coef[2] + torsion_coef[3] * 2.0 * cosine_phi + torsion_coef[4] * 3.0 * (r1 * r1) + 4*torsion_coef[5]*r1*r1*r1) *
308 sine_theta_ante_sq * sine_theta_ante * sine_theta_post_sq * sine_theta_post;
309
310 /* Computation of factors (important for gaining speed). Factors factor_phi_* are coming from the
311 * derivatives of the torsion angle (phi) with respect to the beads ai, aj, al, ak,
312 * (four) coordinates and they are multiplied in the force computations with the
313 * differences of the particles positions stored in parameters delta_ante,
314 * delta_crnt, delta_post. For formulas see file "restcbt.h" */
315
316 factor_phi_ai_ante = ratio_phi_ante * c_self_crnt;
317 factor_phi_ai_crnt = -c_cros_post - ratio_phi_ante * c_cros_ante;
318 factor_phi_ai_post = c_self_crnt;
319 factor_phi_aj_ante = -c_cros_post - ratio_phi_ante * (c_self_crnt + c_cros_ante);
320 factor_phi_aj_crnt = c_cros_post + c_cros_acrs * 2.0 + ratio_phi_ante * (c_self_ante + c_cros_ante) + ratio_phi_post * c_self_post;
321 factor_phi_aj_post = -(c_cros_ante + c_self_crnt) - ratio_phi_post * c_cros_post;
322 factor_phi_ak_ante = c_cros_post + c_self_crnt + ratio_phi_ante * c_cros_ante;
323 factor_phi_ak_crnt = -(c_cros_ante + c_cros_acrs * 2.0) - ratio_phi_ante * c_self_ante - ratio_phi_post * (c_self_post + c_cros_post);
324 factor_phi_ak_post = c_cros_ante + ratio_phi_post * (c_self_crnt + c_cros_post);
325 factor_phi_al_ante = -c_self_crnt;
326 factor_phi_al_crnt = c_cros_ante + ratio_phi_post * c_cros_post;
327 factor_phi_al_post = -ratio_phi_post * c_self_crnt;
328
329 /* Computation of forces due to the derivatives of dihedral angle phi*/
330 for (d = 0; d < DIM3; d++)
331 {
332 f_phi_ai[d] = prefactor_phi * (factor_phi_ai_ante * delta_ante[d] + factor_phi_ai_crnt * delta_crnt[d] + factor_phi_ai_post * delta_post[d]);
333 f_phi_aj[d] = prefactor_phi * (factor_phi_aj_ante * delta_ante[d] + factor_phi_aj_crnt * delta_crnt[d] + factor_phi_aj_post * delta_post[d]);
334 f_phi_ak[d] = prefactor_phi * (factor_phi_ak_ante * delta_ante[d] + factor_phi_ak_crnt * delta_crnt[d] + factor_phi_ak_post * delta_post[d]);
335 f_phi_al[d] = prefactor_phi * (factor_phi_al_ante * delta_ante[d] + factor_phi_al_crnt * delta_crnt[d] + factor_phi_al_post * delta_post[d]);
336 }
337
338 /* PART 3 - COMPUTES THE FORCE COMPONENTS DUE TO THE DERIVATIVES OF BENDING ANGLE THETHA_ANTHE */
339 /* Computation of ratios */
340 ratio_theta_ante_ante = c_cros_ante / c_self_ante;
341 ratio_theta_ante_crnt = c_cros_ante / c_self_crnt;
342
343 /* Computation of the prefactor */
344 /* Computing 2nd power */
345 r1 = cosine_phi;
346 /* Computing 3rd power */
347 r2 = cosine_phi;
348
349 prefactor_theta_ante = -torsion_coef[0] * norm_theta_ante * ( torsion_coef[1] + torsion_coef[2] * cosine_phi + torsion_coef[3] * (r1 * r1) +
350 torsion_coef[4] * (r2 * (r2 * r2))+ torsion_coef[5] * (r2 * (r2 * (r2 * r2)))) * (-3.0) * cosine_theta_ante * sine_theta_ante * sine_theta_post_sq * sine_theta_post;
351
352
353 /* Computation of forces due to the derivatives of bending angle theta_ante */
354 for (d = 0; d < DIM3; d++)
355 {
356 f_theta_ante_ai[d] = prefactor_theta_ante * (ratio_theta_ante_ante * delta_ante[d] - delta_crnt[d]);
357 f_theta_ante_aj[d] = prefactor_theta_ante * ((ratio_theta_ante_crnt + 1.0) * delta_crnt[d] - (ratio_theta_ante_ante + 1.0) * delta_ante[d]);
358 f_theta_ante_ak[d] = prefactor_theta_ante * (delta_ante[d] - ratio_theta_ante_crnt * delta_crnt[d]);
359 }
360
361 /* PART 4 - COMPUTES THE FORCE COMPONENTS DUE TO THE DERIVATIVES OF THE BENDING ANGLE THETA_POST */
362
363 /* Computation of ratios */
364 ratio_theta_post_crnt = c_cros_post / c_self_crnt;
365 ratio_theta_post_post = c_cros_post / c_self_post;
366
367 /* Computation of the prefactor */
368 /* Computing 2nd power */
369 r1 = cosine_phi;
370 /* Computing 3rd power */
371 r2 = cosine_phi;
372
373 prefactor_theta_post = -torsion_coef[0] * norm_theta_post * (torsion_coef[1] + torsion_coef[2] * cosine_phi + torsion_coef[3] * (r1 * r1) +
374 torsion_coef[4] * (r2 * (r2 * r2)) + torsion_coef[5] * (r2 * (r2 * (r2 * r2)))) * sine_theta_ante_sq * sine_theta_ante * (-3.0) * cosine_theta_post * sine_theta_post;
375
376
377 /* Computation of forces due to the derivatives of bending angle Theta_Post */
378 for (d = 0; d < DIM3; d++)
379 {
380 f_theta_post_aj[d] = prefactor_theta_post * (ratio_theta_post_crnt * delta_crnt[d] - delta_post[d]);
381 f_theta_post_ak[d] = prefactor_theta_post * ((ratio_theta_post_post + 1.0) * delta_post[d] - (ratio_theta_post_crnt + 1.0) * delta_crnt[d]);
382 f_theta_post_al[d] = prefactor_theta_post * (delta_crnt[d] - ratio_theta_post_post * delta_post[d]);
383 }
384 r1 = cosine_phi;
385 r2 = cosine_phi;
386
387 /* Contribution to energy - for formula see file "restcbt.h" */
388 *v = torsion_coef[0] * (torsion_coef[1] + torsion_coef[2] * cosine_phi + torsion_coef[3] * (r1 * r1) +
389 torsion_coef[4] * (r2 * (r2 * r2)) + torsion_coef[5] * (r2 * (r2 * (r2 * r2)))) * sine_theta_ante_sq *
390 sine_theta_ante * sine_theta_post_sq * sine_theta_post;
391
392
393}