src/gromacs/mdlib/pull.c

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