Bug Summary

File:gromacs/pulling/pull.c
Location:line 302, column 5
Description:Value stored to 'pgrp1' 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#ifdef HAVE_CONFIG_H1
38#include <config.h>
39#endif
40
41#include <math.h>
42#include <stdio.h>
43#include <stdlib.h>
44#include <string.h>
45
46#include "gromacs/utility/futil.h"
47#include "index.h"
48#include "typedefs.h"
49#include "types/commrec.h"
50#include "network.h"
51#include "pull.h"
52#include "names.h"
53#include "pbc.h"
54#include "mtop_util.h"
55#include "mdrun.h"
56#include "gmx_ga2la.h"
57#include "copyrite.h"
58#include "macros.h"
59
60#include "gromacs/fileio/filenm.h"
61#include "gromacs/fileio/gmxfio.h"
62#include "gromacs/fileio/xvgr.h"
63#include "gromacs/math/vec.h"
64#include "gromacs/utility/smalloc.h"
65
66static void pull_print_group_x(FILE *out, ivec dim, const t_pull_group *pgrp)
67{
68 int m;
69
70 for (m = 0; m < DIM3; m++)
71 {
72 if (dim[m])
73 {
74 fprintf(out, "\t%g", pgrp->x[m]);
75 }
76 }
77}
78
79static void pull_print_coord_dr(FILE *out, ivec dim, const t_pull_coord *pcrd)
80{
81 int m;
82
83 for (m = 0; m < DIM3; m++)
84 {
85 if (dim[m])
86 {
87 fprintf(out, "\t%g", pcrd->dr[m]);
88 }
89 }
90}
91
92static void pull_print_x(FILE *out, t_pull *pull, double t)
93{
94 int c;
95 const t_pull_coord *pcrd;
96
97 fprintf(out, "%.4f", t);
98
99 for (c = 0; c < pull->ncoord; c++)
100 {
101 pcrd = &pull->coord[c];
102
103 if (pull->bPrintRef)
104 {
105 if (PULL_CYL(pull)((pull)->eGeom == epullgCYL))
106 {
107 pull_print_group_x(out, pull->dim, &pull->dyna[c]);
108 }
109 else
110 {
111 pull_print_group_x(out, pull->dim, &pull->group[pcrd->group[0]]);
112 }
113 }
114 pull_print_coord_dr(out, pull->dim, pcrd);
115 }
116 fprintf(out, "\n");
117}
118
119static void pull_print_f(FILE *out, t_pull *pull, double t)
120{
121 int c, d;
122
123 fprintf(out, "%.4f", t);
124
125 for (c = 0; c < pull->ncoord; c++)
126 {
127 fprintf(out, "\t%g", pull->coord[c].f_scal);
128 }
129 fprintf(out, "\n");
130}
131
132void pull_print_output(t_pull *pull, gmx_int64_t step, double time)
133{
134 if ((pull->nstxout != 0) && (step % pull->nstxout == 0))
135 {
136 pull_print_x(pull->out_x, pull, time);
137 }
138
139 if ((pull->nstfout != 0) && (step % pull->nstfout == 0))
140 {
141 pull_print_f(pull->out_f, pull, time);
142 }
143}
144
145static FILE *open_pull_out(const char *fn, t_pull *pull, const output_env_t oenv,
146 gmx_bool bCoord, unsigned long Flags)
147{
148 FILE *fp;
149 int nsets, c, m;
150 char **setname, buf[10];
151
152 if (Flags & MD_APPENDFILES(1<<15))
153 {
154 fp = gmx_fio_fopen(fn, "a+");
155 }
156 else
157 {
158 fp = gmx_fio_fopen(fn, "w+");
159 if (bCoord)
160 {
161 xvgr_header(fp, "Pull COM", "Time (ps)", "Position (nm)",
162 exvggtXNY, oenv);
163 }
164 else
165 {
166 xvgr_header(fp, "Pull force", "Time (ps)", "Force (kJ/mol/nm)",
167 exvggtXNY, oenv);
168 }
169
170 snew(setname, 2*pull->ncoord*DIM)(setname) = save_calloc("setname", "/home/alexxy/Develop/gromacs/src/gromacs/pulling/pull.c"
, 170, (2*pull->ncoord*3), sizeof(*(setname)));
171 nsets = 0;
172 for (c = 0; c < pull->ncoord; c++)
173 {
174 if (bCoord)
175 {
176 if (pull->bPrintRef)
177 {
178 for (m = 0; m < DIM3; m++)
179 {
180 if (pull->dim[m])
181 {
182 sprintf(buf, "%d %s%c", c+1, "c", 'X'+m);
183 setname[nsets] = strdup(buf)(__extension__ (__builtin_constant_p (buf) && ((size_t
)(const void *)((buf) + 1) - (size_t)(const void *)(buf) == 1
) ? (((const char *) (buf))[0] == '\0' ? (char *) calloc ((size_t
) 1, (size_t) 1) : ({ size_t __len = strlen (buf) + 1; char *
__retval = (char *) malloc (__len); if (__retval != ((void*)0
)) __retval = (char *) memcpy (__retval, buf, __len); __retval
; })) : __strdup (buf)));
184 nsets++;
185 }
186 }
187 }
188 for (m = 0; m < DIM3; m++)
189 {
190 if (pull->dim[m])
191 {
192 sprintf(buf, "%d %s%c", c+1, "d", 'X'+m);
193 setname[nsets] = strdup(buf)(__extension__ (__builtin_constant_p (buf) && ((size_t
)(const void *)((buf) + 1) - (size_t)(const void *)(buf) == 1
) ? (((const char *) (buf))[0] == '\0' ? (char *) calloc ((size_t
) 1, (size_t) 1) : ({ size_t __len = strlen (buf) + 1; char *
__retval = (char *) malloc (__len); if (__retval != ((void*)0
)) __retval = (char *) memcpy (__retval, buf, __len); __retval
; })) : __strdup (buf)));
194 nsets++;
195 }
196 }
197 }
198 else
199 {
200 sprintf(buf, "%d", c+1);
201 setname[nsets] = strdup(buf)(__extension__ (__builtin_constant_p (buf) && ((size_t
)(const void *)((buf) + 1) - (size_t)(const void *)(buf) == 1
) ? (((const char *) (buf))[0] == '\0' ? (char *) calloc ((size_t
) 1, (size_t) 1) : ({ size_t __len = strlen (buf) + 1; char *
__retval = (char *) malloc (__len); if (__retval != ((void*)0
)) __retval = (char *) memcpy (__retval, buf, __len); __retval
; })) : __strdup (buf)));
202 nsets++;
203 }
204 }
205 if (nsets > 1)
206 {
207 xvgr_legend(fp, nsets, (const char**)setname, oenv);
208 }
209 for (c = 0; c < nsets; c++)
210 {
211 sfree(setname[c])save_free("setname[c]", "/home/alexxy/Develop/gromacs/src/gromacs/pulling/pull.c"
, 211, (setname[c]));
212 }
213 sfree(setname)save_free("setname", "/home/alexxy/Develop/gromacs/src/gromacs/pulling/pull.c"
, 213, (setname));
214 }
215
216 return fp;
217}
218
219/* Apply forces in a mass weighted fashion */
220static void apply_forces_grp(const t_pull_group *pgrp, const t_mdatoms *md,
221 const dvec f_pull, int sign, rvec *f)
222{
223 int i, ii, m;
224 double wmass, inv_wm;
225
226 inv_wm = pgrp->wscale*pgrp->invtm;
227
228 for (i = 0; i < pgrp->nat_loc; i++)
229 {
230 ii = pgrp->ind_loc[i];
231 wmass = md->massT[ii];
232 if (pgrp->weight_loc)
233 {
234 wmass *= pgrp->weight_loc[i];
235 }
236
237 for (m = 0; m < DIM3; m++)
238 {
239 f[ii][m] += sign * wmass * f_pull[m] * inv_wm;
240 }
241 }
242}
243
244/* Apply forces in a mass weighted fashion */
245static void apply_forces(t_pull * pull, t_mdatoms * md, rvec *f)
246{
247 int c;
248 const t_pull_coord *pcrd;
249
250 for (c = 0; c < pull->ncoord; c++)
251 {
252 pcrd = &pull->coord[c];
253
254 if (PULL_CYL(pull)((pull)->eGeom == epullgCYL))
255 {
256 apply_forces_grp(&pull->dyna[c], md, pcrd->f, -1, f);
257 }
258 else
259 {
260 if (pull->group[pcrd->group[0]].nat > 0)
261 {
262 apply_forces_grp(&pull->group[pcrd->group[0]], md, pcrd->f, -1, f);
263 }
264 }
265 apply_forces_grp(&pull->group[pcrd->group[1]], md, pcrd->f, 1, f);
266 }
267}
268
269static double max_pull_distance2(const t_pull *pull, const t_pbc *pbc)
270{
271 double max_d2;
272 int m;
273
274 max_d2 = GMX_DOUBLE_MAX1.79769312E+308;
275
276 if (pull->eGeom != epullgDIRPBC)
277 {
278 for (m = 0; m < pbc->ndim_ePBC; m++)
279 {
280 if (pull->dim[m] != 0)
281 {
282 max_d2 = min(max_d2, norm2(pbc->box[m]))(((max_d2) < (norm2(pbc->box[m]))) ? (max_d2) : (norm2(
pbc->box[m])) );
283 }
284 }
285 }
286
287 return 0.25*max_d2;
288}
289
290static void low_get_pull_coord_dr(const t_pull *pull,
291 const t_pull_coord *pcrd,
292 const t_pbc *pbc, double t,
293 dvec xg, dvec xref, double max_dist2,
294 dvec dr)
295{
296 const t_pull_group *pgrp0, *pgrp1;
297 int m;
298 dvec xrefr, dref = {0, 0, 0};
299 double dr2;
300
301 pgrp0 = &pull->group[pcrd->group[0]];
302 pgrp1 = &pull->group[pcrd->group[1]];
Value stored to 'pgrp1' is never read
303
304 /* Only the first group can be an absolute reference, in that case nat=0 */
305 if (pgrp0->nat == 0)
306 {
307 for (m = 0; m < DIM3; m++)
308 {
309 xref[m] = pcrd->origin[m];
310 }
311 }
312
313 copy_dvec(xref, xrefr);
314
315 if (pull->eGeom == epullgDIRPBC)
316 {
317 for (m = 0; m < DIM3; m++)
318 {
319 dref[m] = (pcrd->init + pcrd->rate*t)*pcrd->vec[m];
320 }
321 /* Add the reference position, so we use the correct periodic image */
322 dvec_inc(xrefr, dref);
323 }
324
325 pbc_dx_d(pbc, xg, xrefr, dr);
326 dr2 = 0;
327 for (m = 0; m < DIM3; m++)
328 {
329 dr[m] *= pull->dim[m];
330 dr2 += dr[m]*dr[m];
331 }
332 if (max_dist2 >= 0 && dr2 > 0.98*0.98*max_dist2)
333 {
334 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/pulling/pull.c",
334, "Distance between pull groups %d and %d (%f nm) is larger than 0.49 times the box size (%f)",
335 pcrd->group[0], pcrd->group[1], sqrt(dr2), sqrt(max_dist2));
336 }
337
338 if (pull->eGeom == epullgDIRPBC)
339 {
340 dvec_inc(dr, dref);
341 }
342}
343
344static void get_pull_coord_dr(const t_pull *pull,
345 int coord_ind,
346 const t_pbc *pbc, double t,
347 dvec dr)
348{
349 double md2;
350 const t_pull_coord *pcrd;
351
352 if (pull->eGeom == epullgDIRPBC)
353 {
354 md2 = -1;
355 }
356 else
357 {
358 md2 = max_pull_distance2(pull, pbc);
359 }
360
361 pcrd = &pull->coord[coord_ind];
362
363 low_get_pull_coord_dr(pull, pcrd, pbc, t,
364 pull->group[pcrd->group[1]].x,
365 PULL_CYL(pull)((pull)->eGeom == epullgCYL) ? pull->dyna[coord_ind].x : pull->group[pcrd->group[0]].x,
366 md2,
367 dr);
368}
369
370void get_pull_coord_distance(const t_pull *pull,
371 int coord_ind,
372 const t_pbc *pbc, double t,
373 dvec dr, double *dev)
374{
375 static gmx_bool bWarned = FALSE0; /* TODO: this should be fixed for thread-safety,
376 but is fairly benign */
377 const t_pull_coord *pcrd;
378 int m;
379 double ref, drs, inpr;
380
381 pcrd = &pull->coord[coord_ind];
382
383 get_pull_coord_dr(pull, coord_ind, pbc, t, dr);
384
385 ref = pcrd->init + pcrd->rate*t;
386
387 switch (pull->eGeom)
388 {
389 case epullgDIST:
390 /* Pull along the vector between the com's */
391 if (ref < 0 && !bWarned)
392 {
393 fprintf(stderrstderr, "\nPull reference distance for coordinate %d is negative (%f)\n", coord_ind+1, ref);
394 bWarned = TRUE1;
395 }
396 drs = dnorm(dr);
397 if (drs == 0)
398 {
399 /* With no vector we can not determine the direction for the force,
400 * so we set the force to zero.
401 */
402 *dev = 0;
403 }
404 else
405 {
406 /* Determine the deviation */
407 *dev = drs - ref;
408 }
409 break;
410 case epullgDIR:
411 case epullgDIRPBC:
412 case epullgCYL:
413 /* Pull along vec */
414 inpr = 0;
415 for (m = 0; m < DIM3; m++)
416 {
417 inpr += pcrd->vec[m]*dr[m];
418 }
419 *dev = inpr - ref;
420 break;
421 }
422}
423
424void clear_pull_forces(t_pull *pull)
425{
426 int i;
427
428 /* Zeroing the forces is only required for constraint pulling.
429 * It can happen that multiple constraint steps need to be applied
430 * and therefore the constraint forces need to be accumulated.
431 */
432 for (i = 0; i < pull->ncoord; i++)
433 {
434 clear_dvec(pull->coord[i].f);
435 pull->coord[i].f_scal = 0;
436 }
437}
438
439/* Apply constraint using SHAKE */
440static void do_constraint(t_pull *pull, t_pbc *pbc,
441 rvec *x, rvec *v,
442 gmx_bool bMaster, tensor vir,
443 double dt, double t)
444{
445
446 dvec *r_ij; /* x[i] com of i in prev. step. Obeys constr. -> r_ij[i] */
447 dvec unc_ij; /* xp[i] com of i this step, before constr. -> unc_ij */
448 dvec *rnew; /* current 'new' positions of the groups */
449 double *dr_tot; /* the total update of the coords */
450 double ref;
451 dvec vec;
452 double d0, inpr;
453 double lambda, rm, mass, invdt = 0;
454 gmx_bool bConverged_all, bConverged = FALSE0;
455 int niter = 0, g, c, ii, j, m, max_iter = 100;
456 double a;
457 dvec f; /* the pull force */
458 dvec tmp, tmp3;
459 t_pull_group *pdyna, *pgrp0, *pgrp1;
460 t_pull_coord *pcrd;
461
462 snew(r_ij, pull->ncoord)(r_ij) = save_calloc("r_ij", "/home/alexxy/Develop/gromacs/src/gromacs/pulling/pull.c"
, 462, (pull->ncoord), sizeof(*(r_ij)));
463 snew(dr_tot, pull->ncoord)(dr_tot) = save_calloc("dr_tot", "/home/alexxy/Develop/gromacs/src/gromacs/pulling/pull.c"
, 463, (pull->ncoord), sizeof(*(dr_tot)));
464
465 snew(rnew, pull->ngroup)(rnew) = save_calloc("rnew", "/home/alexxy/Develop/gromacs/src/gromacs/pulling/pull.c"
, 465, (pull->ngroup), sizeof(*(rnew)));
466
467 /* copy the current unconstrained positions for use in iterations. We
468 iterate until rinew[i] and rjnew[j] obey the constraints. Then
469 rinew - pull.x_unc[i] is the correction dr to group i */
470 for (g = 0; g < pull->ngroup; g++)
471 {
472 copy_dvec(pull->group[g].xp, rnew[g]);
473 }
474 if (PULL_CYL(pull)((pull)->eGeom == epullgCYL))
475 {
476 /* There is only one pull coordinate and reference group */
477 copy_dvec(pull->dyna[0].xp, rnew[pull->coord[0].group[0]]);
478 }
479
480 /* Determine the constraint directions from the old positions */
481 for (c = 0; c < pull->ncoord; c++)
482 {
483 get_pull_coord_dr(pull, c, pbc, t, r_ij[c]);
484 /* Store the difference vector at time t for printing */
485 copy_dvec(r_ij[c], pull->coord[c].dr);
486 if (debug)
487 {
488 fprintf(debug, "Pull coord %d dr %f %f %f\n",
489 c, r_ij[c][XX0], r_ij[c][YY1], r_ij[c][ZZ2]);
490 }
491
492 if (pull->eGeom == epullgDIR || pull->eGeom == epullgDIRPBC)
493 {
494 /* Select the component along vec */
495 a = 0;
496 for (m = 0; m < DIM3; m++)
497 {
498 a += pull->coord[c].vec[m]*r_ij[c][m];
499 }
500 for (m = 0; m < DIM3; m++)
501 {
502 r_ij[c][m] = a*pull->coord[c].vec[m];
503 }
504 }
505 }
506
507 bConverged_all = FALSE0;
508 while (!bConverged_all && niter < max_iter)
509 {
510 bConverged_all = TRUE1;
511
512 /* loop over all constraints */
513 for (c = 0; c < pull->ncoord; c++)
514 {
515 dvec dr0, dr1;
516
517 pcrd = &pull->coord[c];
518 pgrp0 = &pull->group[pcrd->group[0]];
519 pgrp1 = &pull->group[pcrd->group[1]];
520
521 /* Get the current difference vector */
522 low_get_pull_coord_dr(pull, pcrd, pbc, t,
523 rnew[pcrd->group[1]],
524 rnew[pcrd->group[0]],
525 -1, unc_ij);
526
527 ref = pcrd->init + pcrd->rate*t;
528
529 if (debug)
530 {
531 fprintf(debug, "Pull coord %d, iteration %d\n", c, niter);
532 }
533
534 rm = 1.0/(pgrp0->invtm + pgrp1->invtm);
535
536 switch (pull->eGeom)
537 {
538 case epullgDIST:
539 if (ref <= 0)
540 {
541 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/pulling/pull.c",
541, "The pull constraint reference distance for group %d is <= 0 (%f)", c, ref);
542 }
543
544 {
545 double q, c_a, c_b, c_c;
546
547 c_a = diprod(r_ij[c], r_ij[c]);
548 c_b = diprod(unc_ij, r_ij[c])*2;
549 c_c = diprod(unc_ij, unc_ij) - dsqr(ref);
550
551 if (c_b < 0)
552 {
553 q = -0.5*(c_b - sqrt(c_b*c_b - 4*c_a*c_c));
554 lambda = -q/c_a;
555 }
556 else
557 {
558 q = -0.5*(c_b + sqrt(c_b*c_b - 4*c_a*c_c));
559 lambda = -c_c/q;
560 }
561
562 if (debug)
563 {
564 fprintf(debug,
565 "Pull ax^2+bx+c=0: a=%e b=%e c=%e lambda=%e\n",
566 c_a, c_b, c_c, lambda);
567 }
568 }
569
570 /* The position corrections dr due to the constraints */
571 dsvmul(-lambda*rm*pgrp1->invtm, r_ij[c], dr1);
572 dsvmul( lambda*rm*pgrp0->invtm, r_ij[c], dr0);
573 dr_tot[c] += -lambda*dnorm(r_ij[c]);
574 break;
575 case epullgDIR:
576 case epullgDIRPBC:
577 case epullgCYL:
578 /* A 1-dimensional constraint along a vector */
579 a = 0;
580 for (m = 0; m < DIM3; m++)
581 {
582 vec[m] = pcrd->vec[m];
583 a += unc_ij[m]*vec[m];
584 }
585 /* Select only the component along the vector */
586 dsvmul(a, vec, unc_ij);
587 lambda = a - ref;
588 if (debug)
589 {
590 fprintf(debug, "Pull inpr %e lambda: %e\n", a, lambda);
591 }
592
593 /* The position corrections dr due to the constraints */
594 dsvmul(-lambda*rm*pgrp1->invtm, vec, dr1);
595 dsvmul( lambda*rm*pgrp0->invtm, vec, dr0);
596 dr_tot[c] += -lambda;
597 break;
598 }
599
600 /* DEBUG */
601 if (debug)
602 {
603 int g0, g1;
604
605 g0 = pcrd->group[0];
606 g1 = pcrd->group[1];
607 low_get_pull_coord_dr(pull, pcrd, pbc, t, rnew[g1], rnew[g0], -1, tmp);
608 low_get_pull_coord_dr(pull, pcrd, pbc, t, dr1, dr0, -1, tmp3);
609 fprintf(debug,
610 "Pull cur %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
611 rnew[g0][0], rnew[g0][1], rnew[g0][2],
612 rnew[g1][0], rnew[g1][1], rnew[g1][2], dnorm(tmp));
613 fprintf(debug,
614 "Pull ref %8s %8s %8s %8s %8s %8s d: %8.5f\n",
615 "", "", "", "", "", "", ref);
616 fprintf(debug,
617 "Pull cor %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
618 dr0[0], dr0[1], dr0[2],
619 dr1[0], dr1[1], dr1[2],
620 dnorm(tmp3));
621 } /* END DEBUG */
622
623 /* Update the COMs with dr */
624 dvec_inc(rnew[pcrd->group[1]], dr1);
625 dvec_inc(rnew[pcrd->group[0]], dr0);
626 }
627
628 /* Check if all constraints are fullfilled now */
629 for (c = 0; c < pull->ncoord; c++)
630 {
631 pcrd = &pull->coord[c];
632
633 low_get_pull_coord_dr(pull, pcrd, pbc, t,
634 rnew[pcrd->group[1]],
635 rnew[pcrd->group[0]],
636 -1, unc_ij);
637
638 switch (pull->eGeom)
639 {
640 case epullgDIST:
641 bConverged = fabs(dnorm(unc_ij) - ref) < pull->constr_tol;
642 break;
643 case epullgDIR:
644 case epullgDIRPBC:
645 case epullgCYL:
646 for (m = 0; m < DIM3; m++)
647 {
648 vec[m] = pcrd->vec[m];
649 }
650 inpr = diprod(unc_ij, vec);
651 dsvmul(inpr, vec, unc_ij);
652 bConverged =
653 fabs(diprod(unc_ij, vec) - ref) < pull->constr_tol;
654 break;
655 }
656
657 if (!bConverged)
658 {
659 if (debug)
660 {
661 fprintf(debug, "NOT CONVERGED YET: Group %d:"
662 "d_ref = %f, current d = %f\n",
663 g, ref, dnorm(unc_ij));
664 }
665
666 bConverged_all = FALSE0;
667 }
668 }
669
670 niter++;
671 /* if after all constraints are dealt with and bConverged is still TRUE
672 we're finished, if not we do another iteration */
673 }
674 if (niter > max_iter)
675 {
676 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/pulling/pull.c",
676, "Too many iterations for constraint run: %d", niter);
677 }
678
679 /* DONE ITERATING, NOW UPDATE COORDINATES AND CALC. CONSTRAINT FORCES */
680
681 if (v)
682 {
683 invdt = 1/dt;
684 }
685
686 /* update atoms in the groups */
687 for (g = 0; g < pull->ngroup; g++)
688 {
689 const t_pull_group *pgrp;
690 dvec dr;
691
692 if (PULL_CYL(pull)((pull)->eGeom == epullgCYL) && g == pull->coord[0].group[0])
693 {
694 pgrp = &pull->dyna[0];
695 }
696 else
697 {
698 pgrp = &pull->group[g];
699 }
700
701 /* get the final constraint displacement dr for group g */
702 dvec_sub(rnew[g], pgrp->xp, dr);
703 /* select components of dr */
704 for (m = 0; m < DIM3; m++)
705 {
706 dr[m] *= pull->dim[m];
707 }
708
709 /* update the atom positions */
710 copy_dvec(dr, tmp);
711 for (j = 0; j < pgrp->nat_loc; j++)
712 {
713 ii = pgrp->ind_loc[j];
714 if (pgrp->weight_loc)
715 {
716 dsvmul(pgrp->wscale*pgrp->weight_loc[j], dr, tmp);
717 }
718 for (m = 0; m < DIM3; m++)
719 {
720 x[ii][m] += tmp[m];
721 }
722 if (v)
723 {
724 for (m = 0; m < DIM3; m++)
725 {
726 v[ii][m] += invdt*tmp[m];
727 }
728 }
729 }
730 }
731
732 /* calculate the constraint forces, used for output and virial only */
733 for (c = 0; c < pull->ncoord; c++)
734 {
735 pcrd = &pull->coord[c];
736 pcrd->f_scal = dr_tot[c]/((pull->group[pcrd->group[0]].invtm + pull->group[pcrd->group[1]].invtm)*dt*dt);
737
738 if (vir && bMaster)
739 {
740 double f_invr;
741
742 /* Add the pull contribution to the virial */
743 f_invr = pcrd->f_scal/dnorm(r_ij[c]);
744
745 for (j = 0; j < DIM3; j++)
746 {
747 for (m = 0; m < DIM3; m++)
748 {
749 vir[j][m] -= 0.5*f_invr*r_ij[c][j]*r_ij[c][m];
750 }
751 }
752 }
753 }
754
755 /* finished! I hope. Give back some memory */
756 sfree(r_ij)save_free("r_ij", "/home/alexxy/Develop/gromacs/src/gromacs/pulling/pull.c"
, 756, (r_ij));
757 sfree(dr_tot)save_free("dr_tot", "/home/alexxy/Develop/gromacs/src/gromacs/pulling/pull.c"
, 757, (dr_tot));
758 sfree(rnew)save_free("rnew", "/home/alexxy/Develop/gromacs/src/gromacs/pulling/pull.c"
, 758, (rnew));
759}
760
761/* Pulling with a harmonic umbrella potential or constant force */
762static void do_pull_pot(int ePull,
763 t_pull *pull, t_pbc *pbc, double t, real lambda,
764 real *V, tensor vir, real *dVdl)
765{
766 int c, j, m;
767 double dev, ndr, invdr;
768 real k, dkdl;
769 t_pull_coord *pcrd;
770
771 /* loop over the pull coordinates */
772 *V = 0;
773 *dVdl = 0;
774 for (c = 0; c < pull->ncoord; c++)
775 {
776 pcrd = &pull->coord[c];
777
778 get_pull_coord_distance(pull, c, pbc, t, pcrd->dr, &dev);
779
780 k = (1.0 - lambda)*pcrd->k + lambda*pcrd->kB;
781 dkdl = pcrd->kB - pcrd->k;
782
783 switch (pull->eGeom)
784 {
785 case epullgDIST:
786 ndr = dnorm(pcrd->dr);
787 invdr = 1/ndr;
788 if (ePull == epullUMBRELLA)
789 {
790 pcrd->f_scal = -k*dev;
791 *V += 0.5* k*dsqr(dev);
792 *dVdl += 0.5*dkdl*dsqr(dev);
793 }
794 else
795 {
796 pcrd->f_scal = -k;
797 *V += k*ndr;
798 *dVdl += dkdl*ndr;
799 }
800 for (m = 0; m < DIM3; m++)
801 {
802 pcrd->f[m] = pcrd->f_scal*pcrd->dr[m]*invdr;
803 }
804 break;
805 case epullgDIR:
806 case epullgDIRPBC:
807 case epullgCYL:
808 if (ePull == epullUMBRELLA)
809 {
810 pcrd->f_scal = -k*dev;
811 *V += 0.5* k*dsqr(dev);
812 *dVdl += 0.5*dkdl*dsqr(dev);
813 }
814 else
815 {
816 ndr = 0;
817 for (m = 0; m < DIM3; m++)
818 {
819 ndr += pcrd->vec[m]*pcrd->dr[m];
820 }
821 pcrd->f_scal = -k;
822 *V += k*ndr;
823 *dVdl += dkdl*ndr;
824 }
825 for (m = 0; m < DIM3; m++)
826 {
827 pcrd->f[m] = pcrd->f_scal*pcrd->vec[m];
828 }
829 break;
830 }
831
832 if (vir)
833 {
834 /* Add the pull contribution to the virial */
835 for (j = 0; j < DIM3; j++)
836 {
837 for (m = 0; m < DIM3; m++)
838 {
839 vir[j][m] -= 0.5*pcrd->f[j]*pcrd->dr[m];
840 }
841 }
842 }
843 }
844}
845
846real pull_potential(int ePull, t_pull *pull, t_mdatoms *md, t_pbc *pbc,
847 t_commrec *cr, double t, real lambda,
848 rvec *x, rvec *f, tensor vir, real *dvdlambda)
849{
850 real V, dVdl;
851
852 pull_calc_coms(cr, pull, md, pbc, t, x, NULL((void*)0));
853
854 do_pull_pot(ePull, pull, pbc, t, lambda,
855 &V, pull->bVirial && MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1)) ? vir : NULL((void*)0), &dVdl);
856
857 /* Distribute forces over pulled groups */
858 apply_forces(pull, md, f);
859
860 if (MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1)))
861 {
862 *dvdlambda += dVdl;
863 }
864
865 return (MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1)) ? V : 0.0);
866}
867
868void pull_constraint(t_pull *pull, t_mdatoms *md, t_pbc *pbc,
869 t_commrec *cr, double dt, double t,
870 rvec *x, rvec *xp, rvec *v, tensor vir)
871{
872 pull_calc_coms(cr, pull, md, pbc, t, x, xp);
873
874 do_constraint(pull, pbc, xp, v, pull->bVirial && MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1)), vir, dt, t);
875}
876
877static void make_local_pull_group(gmx_ga2la_t ga2la,
878 t_pull_group *pg, int start, int end)
879{
880 int i, ii;
881
882 pg->nat_loc = 0;
883 for (i = 0; i < pg->nat; i++)
884 {
885 ii = pg->ind[i];
886 if (ga2la)
887 {
888 if (!ga2la_get_home(ga2la, ii, &ii))
889 {
890 ii = -1;
891 }
892 }
893 if (ii >= start && ii < end)
894 {
895 /* This is a home atom, add it to the local pull group */
896 if (pg->nat_loc >= pg->nalloc_loc)
897 {
898 pg->nalloc_loc = over_alloc_dd(pg->nat_loc+1);
899 srenew(pg->ind_loc, pg->nalloc_loc)(pg->ind_loc) = save_realloc("pg->ind_loc", "/home/alexxy/Develop/gromacs/src/gromacs/pulling/pull.c"
, 899, (pg->ind_loc), (pg->nalloc_loc), sizeof(*(pg->
ind_loc)));
900 if (pg->epgrppbc == epgrppbcCOS || pg->weight)
901 {
902 srenew(pg->weight_loc, pg->nalloc_loc)(pg->weight_loc) = save_realloc("pg->weight_loc", "/home/alexxy/Develop/gromacs/src/gromacs/pulling/pull.c"
, 902, (pg->weight_loc), (pg->nalloc_loc), sizeof(*(pg->
weight_loc)));
903 }
904 }
905 pg->ind_loc[pg->nat_loc] = ii;
906 if (pg->weight)
907 {
908 pg->weight_loc[pg->nat_loc] = pg->weight[i];
909 }
910 pg->nat_loc++;
911 }
912 }
913}
914
915void dd_make_local_pull_groups(gmx_domdec_t *dd, t_pull *pull, t_mdatoms *md)
916{
917 gmx_ga2la_t ga2la;
918 int g;
919
920 if (dd)
921 {
922 ga2la = dd->ga2la;
923 }
924 else
925 {
926 ga2la = NULL((void*)0);
927 }
928
929 for (g = 0; g < pull->ngroup; g++)
930 {
931 make_local_pull_group(ga2la, &pull->group[g],
932 0, md->homenr);
933 }
934}
935
936static void init_pull_group_index(FILE *fplog, t_commrec *cr,
937 int g, t_pull_group *pg, ivec pulldims,
938 gmx_mtop_t *mtop, t_inputrec *ir, real lambda)
939{
940 int i, ii, d, nfrozen, ndim;
941 real m, w, mbd;
942 double tmass, wmass, wwmass;
943 gmx_groups_t *groups;
944 gmx_mtop_atomlookup_t alook;
945 t_atom *atom;
946
947 if (EI_ENERGY_MINIMIZATION(ir->eI)((ir->eI) == eiSteep || (ir->eI) == eiCG || (ir->eI)
== eiLBFGS) || ir->eI == eiBD)
948 {
949 /* There are no masses in the integrator.
950 * But we still want to have the correct mass-weighted COMs.
951 * So we store the real masses in the weights.
952 * We do not set nweight, so these weights do not end up in the tpx file.
953 */
954 if (pg->nweight == 0)
955 {
956 snew(pg->weight, pg->nat)(pg->weight) = save_calloc("pg->weight", "/home/alexxy/Develop/gromacs/src/gromacs/pulling/pull.c"
, 956, (pg->nat), sizeof(*(pg->weight)));
957 }
958 }
959
960 if (cr && PAR(cr)((cr)->nnodes > 1))
961 {
962 pg->nat_loc = 0;
963 pg->nalloc_loc = 0;
964 pg->ind_loc = NULL((void*)0);
965 pg->weight_loc = NULL((void*)0);
966 }
967 else
968 {
969 pg->nat_loc = pg->nat;
970 pg->ind_loc = pg->ind;
971 if (pg->epgrppbc == epgrppbcCOS)
972 {
973 snew(pg->weight_loc, pg->nat)(pg->weight_loc) = save_calloc("pg->weight_loc", "/home/alexxy/Develop/gromacs/src/gromacs/pulling/pull.c"
, 973, (pg->nat), sizeof(*(pg->weight_loc)));
974 }
975 else
976 {
977 pg->weight_loc = pg->weight;
978 }
979 }
980
981 groups = &mtop->groups;
982
983 alook = gmx_mtop_atomlookup_init(mtop);
984
985 nfrozen = 0;
986 tmass = 0;
987 wmass = 0;
988 wwmass = 0;
989 for (i = 0; i < pg->nat; i++)
990 {
991 ii = pg->ind[i];
992 gmx_mtop_atomnr_to_atom(alook, ii, &atom);
993 if (ir->opts.nFreeze)
994 {
995 for (d = 0; d < DIM3; d++)
996 {
997 if (pulldims[d] && ir->opts.nFreeze[ggrpnr(groups, egcFREEZE, ii)((groups)->grpnr[egcFREEZE] ? (groups)->grpnr[egcFREEZE
][ii] : 0)][d])
998 {
999 nfrozen++;
1000 }
1001 }
1002 }
1003 if (ir->efep == efepNO)
1004 {
1005 m = atom->m;
1006 }
1007 else
1008 {
1009 m = (1 - lambda)*atom->m + lambda*atom->mB;
1010 }
1011 if (pg->nweight > 0)
1012 {
1013 w = pg->weight[i];
1014 }
1015 else
1016 {
1017 w = 1;
1018 }
1019 if (EI_ENERGY_MINIMIZATION(ir->eI)((ir->eI) == eiSteep || (ir->eI) == eiCG || (ir->eI)
== eiLBFGS))
1020 {
1021 /* Move the mass to the weight */
1022 w *= m;
1023 m = 1;
1024 pg->weight[i] = w;
1025 }
1026 else if (ir->eI == eiBD)
1027 {
1028 if (ir->bd_fric)
1029 {
1030 mbd = ir->bd_fric*ir->delta_t;
1031 }
1032 else
1033 {
1034 if (groups->grpnr[egcTC] == NULL((void*)0))
1035 {
1036 mbd = ir->delta_t/ir->opts.tau_t[0];
1037 }
1038 else
1039 {
1040 mbd = ir->delta_t/ir->opts.tau_t[groups->grpnr[egcTC][ii]];
1041 }
1042 }
1043 w *= m/mbd;
1044 m = mbd;
1045 pg->weight[i] = w;
1046 }
1047 tmass += m;
1048 wmass += m*w;
1049 wwmass += m*w*w;
1050 }
1051
1052 gmx_mtop_atomlookup_destroy(alook);
1053
1054 if (wmass == 0)
1055 {
1056 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/pulling/pull.c",
1056, "The total%s mass of pull group %d is zero",
1057 pg->weight ? " weighted" : "", g);
1058 }
1059 if (fplog)
1060 {
1061 fprintf(fplog,
1062 "Pull group %d: %5d atoms, mass %9.3f", g, pg->nat, tmass);
1063 if (pg->weight || EI_ENERGY_MINIMIZATION(ir->eI)((ir->eI) == eiSteep || (ir->eI) == eiCG || (ir->eI)
== eiLBFGS) || ir->eI == eiBD)
1064 {
1065 fprintf(fplog, ", weighted mass %9.3f", wmass*wmass/wwmass);
1066 }
1067 if (pg->epgrppbc == epgrppbcCOS)
1068 {
1069 fprintf(fplog, ", cosine weighting will be used");
1070 }
1071 fprintf(fplog, "\n");
1072 }
1073
1074 if (nfrozen == 0)
1075 {
1076 /* A value > 0 signals not frozen, it is updated later */
1077 pg->invtm = 1.0;
1078 }
1079 else
1080 {
1081 ndim = 0;
1082 for (d = 0; d < DIM3; d++)
1083 {
1084 ndim += pulldims[d]*pg->nat;
1085 }
1086 if (fplog && nfrozen > 0 && nfrozen < ndim)
1087 {
1088 fprintf(fplog,
1089 "\nWARNING: In pull group %d some, but not all of the degrees of freedom\n"
1090 " that are subject to pulling are frozen.\n"
1091 " For pulling the whole group will be frozen.\n\n",
1092 g);
1093 }
1094 pg->invtm = 0.0;
1095 pg->wscale = 1.0;
1096 }
1097}
1098
1099void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[],
1100 gmx_mtop_t *mtop, t_commrec *cr, const output_env_t oenv, real lambda,
1101 gmx_bool bOutFile, unsigned long Flags)
1102{
1103 t_pull *pull;
1104 t_pull_group *pgrp;
1105 int c, g, start = 0, end = 0, m;
1106
1107 pull = ir->pull;
1108
1109 pull->ePBC = ir->ePBC;
1110 switch (pull->ePBC)
1111 {
1112 case epbcNONE: pull->npbcdim = 0; break;
1113 case epbcXY: pull->npbcdim = 2; break;
1114 default: pull->npbcdim = 3; break;
1115 }
1116
1117 if (fplog)
1118 {
1119 gmx_bool bAbs, bCos;
1120
1121 bAbs = FALSE0;
1122 for (c = 0; c < pull->ncoord; c++)
1123 {
1124 if (pull->group[pull->coord[c].group[0]].nat == 0 ||
1125 pull->group[pull->coord[c].group[1]].nat == 0)
1126 {
1127 bAbs = TRUE1;
1128 }
1129 }
1130
1131 fprintf(fplog, "\nWill apply %s COM pulling in geometry '%s'\n",
1132 EPULLTYPE(ir->ePull)((((ir->ePull) < 0) || ((ir->ePull) >= (epullNR))
) ? "UNDEFINED" : (epull_names)[ir->ePull]), EPULLGEOM(pull->eGeom)((((pull->eGeom) < 0) || ((pull->eGeom) >= (epullgNR
))) ? "UNDEFINED" : (epullg_names)[pull->eGeom]));
1133 fprintf(fplog, "with %d pull coordinate%s and %d group%s\n",
1134 pull->ncoord, pull->ncoord == 1 ? "" : "s",
1135 pull->ngroup, pull->ngroup == 1 ? "" : "s");
1136 if (bAbs)
1137 {
1138 fprintf(fplog, "with an absolute reference\n");
1139 }
1140 bCos = FALSE0;
1141 for (g = 0; g < pull->ngroup; g++)
1142 {
1143 if (pull->group[g].nat > 1 &&
1144 pull->group[g].pbcatom < 0)
1145 {
1146 /* We are using cosine weighting */
1147 fprintf(fplog, "Cosine weighting is used for group %d\n", g);
1148 bCos = TRUE1;
1149 }
1150 }
1151 if (bCos)
1152 {
1153 please_cite(fplog, "Engin2010");
1154 }
1155 }
1156
1157 /* We always add the virial contribution,
1158 * except for geometry = direction_periodic where this is impossible.
1159 */
1160 pull->bVirial = (pull->eGeom != epullgDIRPBC);
1161 if (getenv("GMX_NO_PULLVIR") != NULL((void*)0))
1162 {
1163 if (fplog)
1164 {
1165 fprintf(fplog, "Found env. var., will not add the virial contribution of the COM pull forces\n");
1166 }
1167 pull->bVirial = FALSE0;
1168 }
1169
1170 pull->rbuf = NULL((void*)0);
1171 pull->dbuf = NULL((void*)0);
1172 pull->dbuf_cyl = NULL((void*)0);
1173 pull->bRefAt = FALSE0;
1174 pull->cosdim = -1;
1175 for (g = 0; g < pull->ngroup; g++)
1176 {
1177 pgrp = &pull->group[g];
1178 pgrp->epgrppbc = epgrppbcNONE;
1179 if (pgrp->nat > 0)
1180 {
1181 /* Determine if we need to take PBC into account for calculating
1182 * the COM's of the pull groups.
1183 */
1184 for (m = 0; m < pull->npbcdim; m++)
1185 {
1186 if (pull->dim[m] && pgrp->nat > 1)
1187 {
1188 if (pgrp->pbcatom >= 0)
1189 {
1190 pgrp->epgrppbc = epgrppbcREFAT;
1191 pull->bRefAt = TRUE1;
1192 }
1193 else
1194 {
1195 if (pgrp->weight)
1196 {
1197 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/pulling/pull.c",
1197, "Pull groups can not have relative weights and cosine weighting at same time");
1198 }
1199 pgrp->epgrppbc = epgrppbcCOS;
1200 if (pull->cosdim >= 0 && pull->cosdim != m)
1201 {
1202 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/pulling/pull.c",
1202, "Can only use cosine weighting with pulling in one dimension (use mdp option pull_dim)");
1203 }
1204 pull->cosdim = m;
1205 }
1206 }
1207 }
1208 /* Set the indices */
1209 init_pull_group_index(fplog, cr, g, pgrp, pull->dim, mtop, ir, lambda);
1210 if (PULL_CYL(pull)((pull)->eGeom == epullgCYL) && pgrp->invtm == 0)
1211 {
1212 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/pulling/pull.c",
1212, "Can not have frozen atoms in a cylinder pull group");
1213 }
1214 }
1215 else
1216 {
1217 /* Absolute reference, set the inverse mass to zero */
1218 pgrp->invtm = 0;
1219 pgrp->wscale = 1;
1220 }
1221 }
1222
1223 /* if we use dynamic reference groups, do some initialising for them */
1224 if (PULL_CYL(pull)((pull)->eGeom == epullgCYL))
1225 {
1226 if (ir->ePull == epullCONSTRAINT && pull->ncoord > 1)
1227 {
1228 /* We can't easily update the single reference group with multiple
1229 * constraints. This would require recalculating COMs.
1230 */
1231 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/pulling/pull.c",
1231, "Constraint COM pulling supports only one coordinate with geometry=cylinder, you can use umbrella pulling with multiple coordinates");
1232 }
1233
1234 for (c = 0; c < pull->ncoord; c++)
1235 {
1236 if (pull->group[pull->coord[c].group[0]].nat == 0)
1237 {
1238 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/pulling/pull.c",
1238, "Dynamic reference groups are not supported when using absolute reference!\n");
1239 }
1240 }
1241
1242 snew(pull->dyna, pull->ncoord)(pull->dyna) = save_calloc("pull->dyna", "/home/alexxy/Develop/gromacs/src/gromacs/pulling/pull.c"
, 1242, (pull->ncoord), sizeof(*(pull->dyna)));
1243 }
1244
1245 /* Only do I/O when we are doing dynamics and if we are the MASTER */
1246 pull->out_x = NULL((void*)0);
1247 pull->out_f = NULL((void*)0);
1248 if (bOutFile)
1249 {
1250 if (pull->nstxout > 0)
1251 {
1252 pull->out_x = open_pull_out(opt2fn("-px", nfile, fnm), pull, oenv, TRUE1, Flags);
1253 }
1254 if (pull->nstfout > 0)
1255 {
1256 pull->out_f = open_pull_out(opt2fn("-pf", nfile, fnm), pull, oenv,
1257 FALSE0, Flags);
1258 }
1259 }
1260}
1261
1262void finish_pull(t_pull *pull)
1263{
1264 if (pull->out_x)
1265 {
1266 gmx_fio_fclose(pull->out_x);
1267 }
1268 if (pull->out_f)
1269 {
1270 gmx_fio_fclose(pull->out_f);
1271 }
1272}