src/gromacs/mdlib/pull.c

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