Bug Summary

File:gromacs/mdlib/adress.c
Location:line 549, column 21
Description:Value stored to 'fscal' is never read

Annotated Source Code

1/*
2 * This file is part of the GROMACS molecular simulation package.
3 *
4 * Copyright (c) 2009 Christoph Junghans, Brad Lambeth.
5 * Copyright (c) 2011,2012,2013,2014, by the GROMACS development team, led by
6 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
7 * and including many others, as listed in the AUTHORS file in the
8 * top-level source directory and at http://www.gromacs.org.
9 *
10 * GROMACS is free software; you can redistribute it and/or
11 * modify it under the terms of the GNU Lesser General Public License
12 * as published by the Free Software Foundation; either version 2.1
13 * of the License, or (at your option) any later version.
14 *
15 * GROMACS is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 * Lesser General Public License for more details.
19 *
20 * You should have received a copy of the GNU Lesser General Public
21 * License along with GROMACS; if not, see
22 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
23 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 *
25 * If you want to redistribute modifications to GROMACS, please
26 * consider that scientific software is very special. Version
27 * control is crucial - bugs must be traceable. We will be happy to
28 * consider code for inclusion in the official distribution, but
29 * derived work must not be called official GROMACS. Details are found
30 * in the README & COPYING files - if they are missing, get the
31 * official version at http://www.gromacs.org.
32 *
33 * To help us fund GROMACS development, we humbly ask that you cite
34 * the research papers on the package. Check out http://www.gromacs.org.
35 */
36
37#include "adress.h"
38#include "gromacs/math/utilities.h"
39#include "pbc.h"
40#include "types/simple.h"
41#include "typedefs.h"
42#include "gromacs/math/vec.h"
43
44#include "gromacs/utility/fatalerror.h"
45
46real
47adress_weight(rvec x,
48 int adresstype,
49 real adressr,
50 real adressw,
51 rvec * ref,
52 t_pbc * pbc,
53 t_forcerec * fr )
54{
55 int i;
56 real l2 = adressr+adressw;
57 real sqr_dl, dl;
58 real tmp;
59 rvec dx;
60
61 sqr_dl = 0.0;
62
63 if (pbc)
64 {
65 pbc_dx(pbc, (*ref), x, dx);
66 }
67 else
68 {
69 rvec_sub((*ref), x, dx);
70 }
71
72 switch (adresstype)
73 {
74 case eAdressOff:
75 /* default to explicit simulation */
76 return 1;
77 case eAdressConst:
78 /* constant value for weighting function = adressw */
79 return fr->adress_const_wf;
80 case eAdressXSplit:
81 /* plane through center of ref, varies in x direction */
82 sqr_dl = dx[0]*dx[0];
83 break;
84 case eAdressSphere:
85 /* point at center of ref, assuming cubic geometry */
86 for (i = 0; i < 3; i++)
87 {
88 sqr_dl += dx[i]*dx[i];
89 }
90 break;
91 default:
92 /* default to explicit simulation */
93 return 1;
94 }
95
96 dl = sqrt(sqr_dl);
97
98 /* molecule is coarse grained */
99 if (dl > l2)
100 {
101 return 0;
102 }
103 /* molecule is explicit */
104 else if (dl < adressr)
105 {
106 return 1;
107 }
108 /* hybrid region */
109 else
110 {
111 tmp = cos((dl-adressr)*M_PI3.14159265358979323846/2/adressw);
112 return tmp*tmp;
113 }
114}
115
116void
117update_adress_weights_com(FILE gmx_unused__attribute__ ((unused)) * fplog,
118 int cg0,
119 int cg1,
120 t_block * cgs,
121 rvec x[],
122 t_forcerec * fr,
123 t_mdatoms * mdatoms,
124 t_pbc * pbc)
125{
126 int icg, k, k0, k1, d;
127 real nrcg, inv_ncg, mtot, inv_mtot;
128 atom_id * cgindex;
129 rvec ix;
130 int adresstype;
131 real adressr, adressw;
132 rvec * ref;
133 real * massT;
134 real * wf;
135
136
137 int n_hyb, n_ex, n_cg;
138
139 n_hyb = 0;
140 n_cg = 0;
141 n_ex = 0;
142
143 adresstype = fr->adress_type;
144 adressr = fr->adress_ex_width;
145 adressw = fr->adress_hy_width;
146 massT = mdatoms->massT;
147 wf = mdatoms->wf;
148 ref = &(fr->adress_refs);
149
150
151 /* Since this is center of mass AdResS, the vsite is not guaranteed
152 * to be on the same node as the constructing atoms. Therefore we
153 * loop over the charge groups, calculate their center of mass,
154 * then use this to calculate wf for each atom. This wastes vsite
155 * construction, but it's the only way to assure that the explicit
156 * atoms have the same wf as their vsite. */
157
158#ifdef DEBUG
159 fprintf(fplog, "Calculating center of mass for charge groups %d to %d\n",
160 cg0, cg1);
161#endif
162 cgindex = cgs->index;
163
164 /* Compute the center of mass for all charge groups */
165 for (icg = cg0; (icg < cg1); icg++)
166 {
167 k0 = cgindex[icg];
168 k1 = cgindex[icg+1];
169 nrcg = k1-k0;
170 if (nrcg == 1)
171 {
172 wf[k0] = adress_weight(x[k0], adresstype, adressr, adressw, ref, pbc, fr);
173 if (wf[k0] == 0)
174 {
175 n_cg++;
176 }
177 else if (wf[k0] == 1)
178 {
179 n_ex++;
180 }
181 else
182 {
183 n_hyb++;
184 }
185 }
186 else
187 {
188 mtot = 0.0;
189 for (k = k0; (k < k1); k++)
190 {
191 mtot += massT[k];
192 }
193 if (mtot > 0.0)
194 {
195 inv_mtot = 1.0/mtot;
196
197 clear_rvec(ix);
198 for (k = k0; (k < k1); k++)
199 {
200 for (d = 0; (d < DIM3); d++)
201 {
202 ix[d] += x[k][d]*massT[k];
203 }
204 }
205 for (d = 0; (d < DIM3); d++)
206 {
207 ix[d] *= inv_mtot;
208 }
209 }
210 /* Calculate the center of gravity if the charge group mtot=0 (only vsites) */
211 else
212 {
213 inv_ncg = 1.0/nrcg;
214
215 clear_rvec(ix);
216 for (k = k0; (k < k1); k++)
217 {
218 for (d = 0; (d < DIM3); d++)
219 {
220 ix[d] += x[k][d];
221 }
222 }
223 for (d = 0; (d < DIM3); d++)
224 {
225 ix[d] *= inv_ncg;
226 }
227 }
228
229 /* Set wf of all atoms in charge group equal to wf of com */
230 wf[k0] = adress_weight(ix, adresstype, adressr, adressw, ref, pbc, fr);
231
232 if (wf[k0] == 0)
233 {
234 n_cg++;
235 }
236 else if (wf[k0] == 1)
237 {
238 n_ex++;
239 }
240 else
241 {
242 n_hyb++;
243 }
244
245 for (k = (k0+1); (k < k1); k++)
246 {
247 wf[k] = wf[k0];
248 }
249 }
250 }
251}
252
253void update_adress_weights_atom_per_atom(
254 int cg0,
255 int cg1,
256 t_block * cgs,
257 rvec x[],
258 t_forcerec * fr,
259 t_mdatoms * mdatoms,
260 t_pbc * pbc)
261{
262 int icg, k, k0, k1, d;
263 real nrcg, inv_ncg, mtot, inv_mtot;
264 atom_id * cgindex;
265 rvec ix;
266 int adresstype;
267 real adressr, adressw;
268 rvec * ref;
269 real * massT;
270 real * wf;
271
272
273 int n_hyb, n_ex, n_cg;
274
275 n_hyb = 0;
276 n_cg = 0;
277 n_ex = 0;
278
279 adresstype = fr->adress_type;
280 adressr = fr->adress_ex_width;
281 adressw = fr->adress_hy_width;
282 massT = mdatoms->massT;
283 wf = mdatoms->wf;
284 ref = &(fr->adress_refs);
285
286 cgindex = cgs->index;
287
288 /* Weighting function is determined for each atom individually.
289 * This is an approximation
290 * as in the theory requires an interpolation based on the center of masses.
291 * Should be used with caution */
292
293 for (icg = cg0; (icg < cg1); icg++)
294 {
295 k0 = cgindex[icg];
296 k1 = cgindex[icg + 1];
297 nrcg = k1 - k0;
298
299 for (k = (k0); (k < k1); k++)
300 {
301 wf[k] = adress_weight(x[k], adresstype, adressr, adressw, ref, pbc, fr);
302 if (wf[k] == 0)
303 {
304 n_cg++;
305 }
306 else if (wf[k] == 1)
307 {
308 n_ex++;
309 }
310 else
311 {
312 n_hyb++;
313 }
314 }
315
316 }
317}
318
319void
320update_adress_weights_cog(t_iparams ip[],
321 t_ilist ilist[],
322 rvec x[],
323 t_forcerec * fr,
324 t_mdatoms * mdatoms,
325 t_pbc * pbc)
326{
327 int i, j, k, nr, nra, inc;
328 int ftype, adresstype;
329 t_iatom avsite, ai, aj, ak, al;
330 t_iatom * ia;
331 real adressr, adressw;
332 rvec * ref;
333 real * wf;
334 int n_hyb, n_ex, n_cg;
335
336 adresstype = fr->adress_type;
337 adressr = fr->adress_ex_width;
338 adressw = fr->adress_hy_width;
339 wf = mdatoms->wf;
340 ref = &(fr->adress_refs);
341
342
343 n_hyb = 0;
344 n_cg = 0;
345 n_ex = 0;
346
347
348 /* Since this is center of geometry AdResS, we know the vsite
349 * is in the same charge group node as the constructing atoms.
350 * Loop over vsite types, calculate the weight of the vsite,
351 * then assign that weight to the constructing atoms. */
352
353 for (ftype = 0; (ftype < F_NRE); ftype++)
354 {
355 if (interaction_function[ftype].flags & IF_VSITE1<<1)
356 {
357 nra = interaction_function[ftype].nratoms;
358 nr = ilist[ftype].nr;
359 ia = ilist[ftype].iatoms;
360
361 for (i = 0; (i < nr); )
362 {
363 /* The vsite and first constructing atom */
364 avsite = ia[1];
365 ai = ia[2];
366 wf[avsite] = adress_weight(x[avsite], adresstype, adressr, adressw, ref, pbc, fr);
367 wf[ai] = wf[avsite];
368
369 if (wf[ai] == 0)
370 {
371 n_cg++;
372 }
373 else if (wf[ai] == 1)
374 {
375 n_ex++;
376 }
377 else
378 {
379 n_hyb++;
380 }
381
382 /* Assign the vsite wf to rest of constructing atoms depending on type */
383 inc = nra+1;
384 switch (ftype)
385 {
386 case F_VSITE2:
387 aj = ia[3];
388 wf[aj] = wf[avsite];
389 break;
390 case F_VSITE3:
391 aj = ia[3];
392 wf[aj] = wf[avsite];
393 ak = ia[4];
394 wf[ak] = wf[avsite];
395 break;
396 case F_VSITE3FD:
397 aj = ia[3];
398 wf[aj] = wf[avsite];
399 ak = ia[4];
400 wf[ak] = wf[avsite];
401 break;
402 case F_VSITE3FAD:
403 aj = ia[3];
404 wf[aj] = wf[avsite];
405 ak = ia[4];
406 wf[ak] = wf[avsite];
407 break;
408 case F_VSITE3OUT:
409 aj = ia[3];
410 wf[aj] = wf[avsite];
411 ak = ia[4];
412 wf[ak] = wf[avsite];
413 break;
414 case F_VSITE4FD:
415 aj = ia[3];
416 wf[aj] = wf[avsite];
417 ak = ia[4];
418 wf[ak] = wf[avsite];
419 al = ia[5];
420 wf[al] = wf[avsite];
421 break;
422 case F_VSITE4FDN:
423 aj = ia[3];
424 wf[aj] = wf[avsite];
425 ak = ia[4];
426 wf[ak] = wf[avsite];
427 al = ia[5];
428 wf[al] = wf[avsite];
429 break;
430 case F_VSITEN:
431 inc = 3*ip[ia[0]].vsiten.n;
432 for (j = 3; j < inc; j += 3)
433 {
434 ai = ia[j+2];
435 wf[ai] = wf[avsite];
436 }
437 break;
438 default:
439 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/mdlib/adress.c",
439, "No such vsite type %d in %s, line %d",
440 ftype, __FILE__"/home/alexxy/Develop/gromacs/src/gromacs/mdlib/adress.c", __LINE__440);
441 }
442
443 /* Increment loop variables */
444 i += inc;
445 ia += inc;
446 }
447 }
448 }
449}
450
451void
452update_adress_weights_atom(int cg0,
453 int cg1,
454 t_block * cgs,
455 rvec x[],
456 t_forcerec * fr,
457 t_mdatoms * mdatoms,
458 t_pbc * pbc)
459{
460 int icg, k, k0, k1;
461 atom_id * cgindex;
462 int adresstype;
463 real adressr, adressw;
464 rvec * ref;
465 real * massT;
466 real * wf;
467
468 adresstype = fr->adress_type;
469 adressr = fr->adress_ex_width;
470 adressw = fr->adress_hy_width;
471 massT = mdatoms->massT;
472 wf = mdatoms->wf;
473 ref = &(fr->adress_refs);
474 cgindex = cgs->index;
475
476 /* Only use first atom in charge group.
477 * We still can't be sure that the vsite and constructing
478 * atoms are on the same processor, so we must calculate
479 * in the same way as com adress. */
480
481 for (icg = cg0; (icg < cg1); icg++)
482 {
483 k0 = cgindex[icg];
484 k1 = cgindex[icg+1];
485 wf[k0] = adress_weight(x[k0], adresstype, adressr, adressw, ref, pbc, fr);
486
487 /* Set wf of all atoms in charge group equal to wf of first atom in charge group*/
488 for (k = (k0+1); (k < k1); k++)
489 {
490 wf[k] = wf[k0];
491 }
492 }
493}
494
495void
496adress_thermo_force(int start,
497 int homenr,
498 t_block * cgs,
499 rvec x[],
500 rvec f[],
501 t_forcerec * fr,
502 t_mdatoms * mdatoms,
503 t_pbc * pbc)
504{
505 int iatom, n0, nnn, nrcg, i;
506 int adresstype;
507 real adressw, adressr;
508 atom_id * cgindex;
509 unsigned short * ptype;
510 rvec * ref;
511 real * wf;
512 real tabscale;
513 real * ATFtab;
514 rvec dr;
515 real w, wsq, wmin1, wmin1sq, wp, wt, rinv, sqr_dl, dl;
516 real eps, eps2, F, Geps, Heps2, Fp, dmu_dwp, dwp_dr, fscal;
517
518 adresstype = fr->adress_type;
519 adressw = fr->adress_hy_width;
520 adressr = fr->adress_ex_width;
521 cgindex = cgs->index;
522 ptype = mdatoms->ptype;
523 ref = &(fr->adress_refs);
524 wf = mdatoms->wf;
525
526 for (iatom = start; (iatom < start+homenr); iatom++)
527 {
528 if (egp_coarsegrained(fr, mdatoms->cENER[iatom]))
529 {
530 if (ptype[iatom] == eptVSite)
531 {
532 w = wf[iatom];
533 /* is it hybrid or apply the thermodynamics force everywhere?*/
534 if (mdatoms->tf_table_index[iatom] != NO_TF_TABLE255)
535 {
536 if (fr->n_adress_tf_grps > 0)
537 {
538 /* multi component tf is on, select the right table */
539 ATFtab = fr->atf_tabs[mdatoms->tf_table_index[iatom]].data;
540 tabscale = fr->atf_tabs[mdatoms->tf_table_index[iatom]].scale;
541 }
542 else
543 {
544 /* just on component*/
545 ATFtab = fr->atf_tabs[DEFAULT_TF_TABLE0].data;
546 tabscale = fr->atf_tabs[DEFAULT_TF_TABLE0].scale;
547 }
548
549 fscal = 0;
Value stored to 'fscal' is never read
550 if (pbc)
551 {
552 pbc_dx(pbc, (*ref), x[iatom], dr);
553 }
554 else
555 {
556 rvec_sub((*ref), x[iatom], dr);
557 }
558
559
560
561
562 /* calculate distace to adress center again */
563 sqr_dl = 0.0;
564 switch (adresstype)
565 {
566 case eAdressXSplit:
567 /* plane through center of ref, varies in x direction */
568 sqr_dl = dr[0]*dr[0];
569 rinv = gmx_invsqrt(dr[0]*dr[0])gmx_software_invsqrt(dr[0]*dr[0]);
570 break;
571 case eAdressSphere:
572 /* point at center of ref, assuming cubic geometry */
573 for (i = 0; i < 3; i++)
574 {
575 sqr_dl += dr[i]*dr[i];
576 }
577 rinv = gmx_invsqrt(sqr_dl)gmx_software_invsqrt(sqr_dl);
578 break;
579 default:
580 /* This case should not happen */
581 rinv = 0.0;
582 }
583
584 dl = sqrt(sqr_dl);
585 /* table origin is adress center */
586 wt = dl*tabscale;
587 n0 = wt;
588 eps = wt-n0;
589 eps2 = eps*eps;
590 nnn = 4*n0;
591 F = ATFtab[nnn+1];
592 Geps = eps*ATFtab[nnn+2];
593 Heps2 = eps2*ATFtab[nnn+3];
594 Fp = F+Geps+Heps2;
595 F = (Fp+Geps+2.0*Heps2)*tabscale;
596
597 fscal = F*rinv;
598
599 f[iatom][0] += fscal*dr[0];
600 if (adresstype != eAdressXSplit)
601 {
602 f[iatom][1] += fscal*dr[1];
603 f[iatom][2] += fscal*dr[2];
604 }
605 }
606 }
607 }
608 }
609}
610
611gmx_bool egp_explicit(t_forcerec * fr, int egp_nr)
612{
613 return fr->adress_group_explicit[egp_nr];
614}
615
616gmx_bool egp_coarsegrained(t_forcerec * fr, int egp_nr)
617{
618 return !fr->adress_group_explicit[egp_nr];
619}