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