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