File: | gromacs/mdlib/adress.c |
Location: | line 471, column 5 |
Description: | Value stored to 'massT' is never read |
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 | |
46 | real |
47 | adress_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 | |
116 | void |
117 | update_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 | |
253 | void 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 | |
319 | void |
320 | update_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 | |
451 | void |
452 | update_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; |
Value stored to 'massT' is never read | |
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 | |
495 | void |
496 | adress_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; |
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 | |
611 | gmx_bool egp_explicit(t_forcerec * fr, int egp_nr) |
612 | { |
613 | return fr->adress_group_explicit[egp_nr]; |
614 | } |
615 | |
616 | gmx_bool egp_coarsegrained(t_forcerec * fr, int egp_nr) |
617 | { |
618 | return !fr->adress_group_explicit[egp_nr]; |
619 | } |