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