src/gromacs/pulling/pull.c

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