src/programs/mdrun/do_gct.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include "typedefs.h"
  40 #include "xmdrun.h"
  41 #include "futil.h"
  42 #include "xvgr.h"
  43 #include "macros.h"
  44 #include "physics.h"
  45 #include "network.h"
  46 #include "smalloc.h"
  47 #include "mdrun.h"
  48 #include "string2.h"
  49 #include "readinp.h"
  50 #include "filenm.h"
  51 #include "update.h"
  52 #include "vec.h"
  53 #include "main.h"
  54 #include "txtdump.h"
  55
  56 /*#define DEBUGGCT*/
  57 t_coupl_rec *init_coupling(FILE *log, int nfile, const t_filenm fnm[],
  58                            t_commrec *cr, t_forcerec *fr,
  59                            t_mdatoms *md, t_idef *idef)
  60 {
  61     int          i, nc, index, j;
  62     int          ati, atj;
  63     t_coupl_rec *tcr;
  64
  65     snew(tcr, 1);
  66     if (MASTER(cr))
  67     {
  68         read_gct (opt2fn("-j", nfile, fnm), tcr);
  69         write_gct(opt2fn("-jo", nfile, fnm), tcr, idef);
  70     }
  71     /* Update all processors with coupling info */
  72     if (PAR(cr))
  73     {
  74         comm_tcr(log, cr, &tcr);
  75     }
  76
  77     /* Copy information from the coupling to the force field stuff */
  78     copy_ff(tcr, fr, md, idef);
  79
  80     return tcr;
  81 }
  82
  83 static real Ecouple(t_coupl_rec *tcr, real ener[])
  84 {
  85     if (tcr->bInter)
  86     {
  87         return ener[F_COUL_SR]+ener[F_LJ]+ener[F_COUL_LR]+ener[F_LJ_LR]+
  88                ener[F_RF_EXCL]+ener[F_COUL_RECIP];
  89     }
  90     else
  91     {
  92         return ener[F_EPOT];
  93     }
  94 }
  95
  96 static char *mk_gct_nm(const char *fn, int ftp, int ati, int atj)
  97 {
  98     static char buf[256];
  99
 100     strcpy(buf, fn);
 101     if (atj == -1)
 102     {
 103         sprintf(buf+strlen(fn)-4, "%d.%s", ati, ftp2ext(ftp));
 104     }
 105     else
 106     {
 107         sprintf(buf+strlen(fn)-4, "%d_%d.%s", ati, atj, ftp2ext(ftp));
 108     }
 109
 110     return buf;
 111 }
 112
 113 static void pr_ff(t_coupl_rec *tcr, real time, t_idef *idef,
 114                   t_commrec *cr, int nfile, const t_filenm fnm[],
 115                   const output_env_t oenv)
 116 {
 117     static FILE  *prop;
 118     static FILE **out = NULL;
 119     static FILE **qq  = NULL;
 120     static FILE **ip  = NULL;
 121     t_coupl_LJ   *tclj;
 122     t_coupl_BU   *tcbu;
 123     char          buf[256];
 124     const char   *leg[]  =  { "C12", "C6" };
 125     const char   *eleg[] =  { "Epsilon", "Sigma" };
 126     const char   *bleg[] = { "A", "B", "C" };
 127     char        **raleg;
 128     int           i, j, index;
 129
 130     if ((prop == NULL) && (out == NULL) && (qq == NULL) && (ip == NULL))
 131     {
 132         prop = xvgropen(opt2fn("-runav", nfile, fnm),
 133                         "Properties and Running Averages", "Time (ps)", "", oenv);
 134         snew(raleg, 2*eoObsNR);
 135         for (i = j = 0; (i < eoObsNR); i++)
 136         {
 137             if (tcr->bObsUsed[i])
 138             {
 139                 raleg[j++] = strdup(eoNames[i]);
 140                 sprintf(buf, "RA-%s", eoNames[i]);
 141                 raleg[j++] = strdup(buf);
 142             }
 143         }
 144         xvgr_legend(prop, j, (const char**)raleg, oenv);
 145         for (i = 0; (i < j); i++)
 146         {
 147             sfree(raleg[i]);
 148         }
 149         sfree(raleg);
 150
 151         if (tcr->nLJ)
 152         {
 153             snew(out, tcr->nLJ);
 154             for (i = 0; (i < tcr->nLJ); i++)
 155             {
 156                 if (tcr->tcLJ[i].bPrint)
 157                 {
 158                     tclj   = &(tcr->tcLJ[i]);
 159                     out[i] =
 160                         xvgropen(mk_gct_nm(opt2fn("-ffout", nfile, fnm),
 161                                            efXVG, tclj->at_i, tclj->at_j),
 162                                  "General Coupling Lennard Jones", "Time (ps)",
 163                                  "Force constant (units)", oenv);
 164                     fprintf(out[i], "@ subtitle \"Interaction between types %d and %d\"\n",
 165                             tclj->at_i, tclj->at_j);
 166                     if (tcr->combrule == 1)
 167                     {
 168                         xvgr_legend(out[i], asize(leg), leg, oenv);
 169                     }
 170                     else
 171                     {
 172                         xvgr_legend(out[i], asize(eleg), eleg, oenv);
 173                     }
 174                     fflush(out[i]);
 175                 }
 176             }
 177         }
 178         else if (tcr->nBU)
 179         {
 180             snew(out, tcr->nBU);
 181             for (i = 0; (i < tcr->nBU); i++)
 182             {
 183                 if (tcr->tcBU[i].bPrint)
 184                 {
 185                     tcbu   = &(tcr->tcBU[i]);
 186                     out[i] =
 187                         xvgropen(mk_gct_nm(opt2fn("-ffout", nfile, fnm), efXVG,
 188                                            tcbu->at_i, tcbu->at_j),
 189                                  "General Coupling Buckingham", "Time (ps)",
 190                                  "Force constant (units)", oenv);
 191                     fprintf(out[i], "@ subtitle \"Interaction between types %d and %d\"\n",
 192                             tcbu->at_i, tcbu->at_j);
 193                     xvgr_legend(out[i], asize(bleg), bleg, oenv);
 194                     fflush(out[i]);
 195                 }
 196             }
 197         }
 198         snew(qq, tcr->nQ);
 199         for (i = 0; (i < tcr->nQ); i++)
 200         {
 201             if (tcr->tcQ[i].bPrint)
 202             {
 203                 qq[i] = xvgropen(mk_gct_nm(opt2fn("-ffout", nfile, fnm), efXVG,
 204                                            tcr->tcQ[i].at_i, -1),
 205                                  "General Coupling Charge", "Time (ps)", "Charge (e)",
 206                                  oenv);
 207                 fprintf(qq[i], "@ subtitle \"Type %d\"\n", tcr->tcQ[i].at_i);
 208                 fflush(qq[i]);
 209             }
 210         }
 211         snew(ip, tcr->nIP);
 212         for (i = 0; (i < tcr->nIP); i++)
 213         {
 214             sprintf(buf, "gctIP%d", tcr->tIP[i].type);
 215             ip[i] = xvgropen(mk_gct_nm(opt2fn("-ffout", nfile, fnm), efXVG, 0, -1),
 216                              "General Coupling iparams", "Time (ps)", "ip ()", oenv);
 217             index = tcr->tIP[i].type;
 218             fprintf(ip[i], "@ subtitle \"Coupling to %s\"\n",
 219                     interaction_function[idef->functype[index]].longname);
 220             fflush(ip[i]);
 221         }
 222     }
 223     /* Write properties to file */
 224     fprintf(prop, "%10.3f", time);
 225     for (i = 0; (i < eoObsNR); i++)
 226     {
 227         if (tcr->bObsUsed[i])
 228         {
 229             fprintf(prop, "  %10.3e  %10.3e", tcr->act_value[i], tcr->av_value[i]);
 230         }
 231     }
 232     fprintf(prop, "\n");
 233     fflush(prop);
 234
 235     for (i = 0; (i < tcr->nLJ); i++)
 236     {
 237         tclj = &(tcr->tcLJ[i]);
 238         if (tclj->bPrint)
 239         {
 240             if (tcr->combrule == 1)
 241             {
 242                 fprintf(out[i], "%14.7e  %14.7e  %14.7e\n",
 243                         time, tclj->c12, tclj->c6);
 244             }
 245             else
 246             {
 247                 double sigma   = pow(tclj->c12/tclj->c6, 1.0/6.0);
 248                 double epsilon = 0.25*sqr(tclj->c6)/tclj->c12;
 249                 fprintf(out[i], "%14.7e  %14.7e  %14.7e\n",
 250                         time, epsilon, sigma);
 251             }
 252             fflush(out[i]);
 253         }
 254     }
 255     for (i = 0; (i < tcr->nBU); i++)
 256     {
 257         tcbu = &(tcr->tcBU[i]);
 258         if (tcbu->bPrint)
 259         {
 260             fprintf(out[i], "%14.7e  %14.7e  %14.7e  %14.7e\n",
 261                     time, tcbu->a, tcbu->b, tcbu->c);
 262             fflush(out[i]);
 263         }
 264     }
 265     for (i = 0; (i < tcr->nQ); i++)
 266     {
 267         if (tcr->tcQ[i].bPrint)
 268         {
 269             fprintf(qq[i], "%14.7e  %14.7e\n", time, tcr->tcQ[i].Q);
 270             fflush(qq[i]);
 271         }
 272     }
 273     for (i = 0; (i < tcr->nIP); i++)
 274     {
 275         fprintf(ip[i], "%10g  ", time);
 276         index = tcr->tIP[i].type;
 277         switch (idef->functype[index])
 278         {
 279             case F_BONDS:
 280                 fprintf(ip[i], "%10g  %10g\n", tcr->tIP[i].iprint.harmonic.krA,
 281                         tcr->tIP[i].iprint.harmonic.rA);
 282                 break;
 283             default:
 284                 break;
 285         }
 286         fflush(ip[i]);
 287     }
 288 }
 289
 290 static void pr_dev(t_coupl_rec *tcr,
 291                    real t, real dev[eoObsNR], t_commrec *cr, int nfile,
 292                    const t_filenm fnm[], const output_env_t oenv)
 293 {
 294     static FILE *fp = NULL;
 295     char       **ptr;
 296     int          i, j;
 297
 298     if (!fp)
 299     {
 300         fp = xvgropen(opt2fn("-devout", nfile, fnm),
 301                       "Deviations from target value", "Time (ps)", "", oenv);
 302         snew(ptr, eoObsNR);
 303         for (i = j = 0; (i < eoObsNR); i++)
 304         {
 305             if (tcr->bObsUsed[i])
 306             {
 307                 ptr[j++] = strdup(eoNames[i]);
 308             }
 309         }
 310         xvgr_legend(fp, j, (const char**)ptr, oenv);
 311         for (i = 0; i < j; i++)
 312         {
 313             sfree(ptr[i]);
 314         }
 315         sfree(ptr);
 316     }
 317     fprintf(fp, "%10.3f", t);
 318     for (i = 0; (i < eoObsNR); i++)
 319     {
 320         if (tcr->bObsUsed[i])
 321         {
 322             fprintf(fp, "  %10.3e", dev[i]);
 323         }
 324     }
 325     fprintf(fp, "\n");
 326     fflush(fp);
 327 }
 328
 329 static void upd_nbfplj(FILE *log, real *nbfp, int atnr, real f6[], real f12[],
 330                        int combrule)
 331 {
 332     double *sigma, *epsilon, c6, c12, eps, sig, sig6;
 333     int     n, m, k;
 334
 335     /* Update the nonbonded force parameters */
 336     switch (combrule)
 337     {
 338         case 1:
 339             for (k = n = 0; (n < atnr); n++)
 340             {
 341                 for (m = 0; (m < atnr); m++, k++)
 342                 {
 343                     C6 (nbfp, atnr, n, m) *= f6[k];
 344                     C12(nbfp, atnr, n, m) *= f12[k];
 345                 }
 346             }
 347             break;
 348         case 2:
 349         case 3:
 350             /* Convert to sigma and epsilon */
 351             snew(sigma, atnr);
 352             snew(epsilon, atnr);
 353             for (n = 0; (n < atnr); n++)
 354             {
 355                 k   = n*(atnr+1);
 356                 c6  = C6 (nbfp, atnr, n, n) * f6[k];
 357                 c12 = C12(nbfp, atnr, n, n) * f12[k];
 358                 if ((c6 == 0) || (c12 == 0))
 359                 {
 360                     gmx_fatal(FARGS, "You can not use combination rule %d with zero C6 (%f) or C12 (%f)", combrule, c6, c12);
 361                 }
 362                 sigma[n]   = pow(c12/c6, 1.0/6.0);
 363                 epsilon[n] = 0.25*(c6*c6/c12);
 364             }
 365             for (k = n = 0; (n < atnr); n++)
 366             {
 367                 for (m = 0; (m < atnr); m++, k++)
 368                 {
 369                     eps  = sqrt(epsilon[n]*epsilon[m]);
 370                     if (combrule == 2)
 371                     {
 372                         sig  = 0.5*(sigma[n]+sigma[m]);
 373                     }
 374                     else
 375                     {
 376                         sig  = sqrt(sigma[n]*sigma[m]);
 377                     }
 378                     sig6 = pow(sig, 6.0);
 379                     /* nbfp now includes the 6.0/12.0 derivative prefactors */
 380                     C6 (nbfp, atnr, n, m) = 4*eps*sig6/6.0;
 381                     C12(nbfp, atnr, n, m) = 4*eps*sig6*sig6/12.0;
 382                 }
 383             }
 384             sfree(sigma);
 385             sfree(epsilon);
 386             break;
 387         default:
 388             gmx_fatal(FARGS, "Combination rule should be 1,2 or 3 instead of %d",
 389                       combrule);
 390     }
 391 }
 392
 393 static void upd_nbfpbu(FILE *log, real *nbfp, int atnr,
 394                        real fa[], real fb[], real fc[])
 395 {
 396     int n, m, k;
 397
 398     /* Update the nonbonded force parameters */
 399     for (k = n = 0; (n < atnr); n++)
 400     {
 401         for (m = 0; (m < atnr); m++, k++)
 402         {
 403             BHAMA(nbfp, atnr, n, m) *= fa[k];
 404             BHAMB(nbfp, atnr, n, m) *= fb[k];
 405             BHAMC(nbfp, atnr, n, m) *= fc[k];
 406         }
 407     }
 408 }
 409
 410 void gprod(t_commrec *cr, int n, real f[])
 411 {
 412     /* Compute the global product of all elements in an array
 413      * such that after gprod f[i] = PROD_j=1,nprocs f[i][j]
 414      */
 415     static int   nbuf = 0;
 416     static real *buf  = NULL;
 417     int          i;
 418
 419     if (n > nbuf)
 420     {
 421         nbuf = n;
 422         srenew(buf, nbuf);
 423     }
 424
 425 #ifdef GMX_MPI
 426 #ifdef GMX_DOUBLE
 427     MPI_Allreduce(f, buf, n, MPI_DOUBLE, MPI_PROD, cr->mpi_comm_mygroup);
 428 #else
 429     MPI_Allreduce(f, buf, n, MPI_FLOAT, MPI_PROD, cr->mpi_comm_mygroup);
 430 #endif
 431     for (i = 0; (i < n); i++)
 432     {
 433         f[i] = buf[i];
 434     }
 435 #endif
 436 }
 437
 438 static void set_factor_matrix(int ntypes, real f[], real fmult, int ati, int atj)
 439 {
 440 #define FMIN 0.95
 441 #define FMAX 1.05
 442     int i;
 443
 444     fmult = min(FMAX, max(FMIN, fmult));
 445     if (atj != -1)
 446     {
 447         f[ntypes*ati+atj] *= fmult;
 448         f[ntypes*atj+ati] *= fmult;
 449     }
 450     else
 451     {
 452         for (i = 0; (i < ntypes); i++)
 453         {
 454             f[ntypes*ati+i] *= fmult;
 455             f[ntypes*i+ati] *= fmult;
 456         }
 457     }
 458 #undef FMIN
 459 #undef FMAX
 460 }
 461
 462 static real calc_deviation(real xav, real xt, real x0)
 463 {
 464     /* This may prevent overshooting in GCT coupling... */
 465
 466     /* real dev;
 467
 468        if (xav > x0) {
 469        if (xt > x0)
 470         dev = min(xav-x0,xt-x0);
 471        else
 472         dev = 0;
 473        }
 474        else {
 475        if (xt < x0)
 476         dev = max(xav-x0,xt-x0);
 477        else
 478         dev = 0;
 479        }
 480      */
 481     return x0-xav;
 482 }
 483
 484 static real calc_dist(FILE *log, rvec x[])
 485 {
 486     static gmx_bool bFirst = TRUE;
 487     static gmx_bool bDist;
 488     static int      i1, i2;
 489     char           *buf;
 490     rvec            dx;
 491
 492     if (bFirst)
 493     {
 494         if ((buf = getenv("DISTGCT")) == NULL)
 495         {
 496             bDist = FALSE;
 497         }
 498         else
 499         {
 500             bDist  = (sscanf(buf, "%d%d", &i1, &i2) == 2);
 501             if (bDist)
 502             {
 503                 fprintf(log, "GCT: Will couple to distance between %d and %d\n", i1, i2);
 504             }
 505             else
 506             {
 507                 fprintf(log, "GCT: Will not couple to distances\n");
 508             }
 509         }
 510         bFirst = FALSE;
 511     }
 512     if (bDist)
 513     {
 514         rvec_sub(x[i1], x[i2], dx);
 515         return norm(dx);
 516     }
 517     else
 518     {
 519         return 0.0;
 520     }
 521 }
 522
 523 real run_aver(real old, real cur, int step, int nmem)
 524 {
 525     nmem   = max(1, nmem);
 526
 527     return ((nmem-1)*old+cur)/nmem;
 528 }
 529
 530 static void set_act_value(t_coupl_rec *tcr, int index, real val, int step)
 531 {
 532     tcr->act_value[index] = val;
 533     tcr->av_value[index]  = run_aver(tcr->av_value[index], val, step, tcr->nmemory);
 534 }
 535
 536 static void upd_f_value(FILE *log, int atnr, real xi, real dt, real factor,
 537                         real ff[], int ati, int atj)
 538 {
 539     real fff;
 540
 541     if (xi != 0)
 542     {
 543         fff = 1 + (dt/xi)  * factor;
 544         if (fff > 0)
 545         {
 546             set_factor_matrix(atnr, ff, sqrt(fff), ati, atj);
 547         }
 548     }
 549 }
 550
 551 static void dump_fm(FILE *fp, int n, real f[], char *s)
 552 {
 553     int i, j;
 554
 555     fprintf(fp, "Factor matrix (all numbers -1) %s\n", s);
 556     for (i = 0; (i < n); i++)
 557     {
 558         for (j = 0; (j < n); j++)
 559         {
 560             fprintf(fp, "  %10.3e", f[n*i+j]-1.0);
 561         }
 562         fprintf(fp, "\n");
 563     }
 564 }
 565
 566 void do_coupling(FILE *log, const output_env_t oenv, int nfile,
 567                  const t_filenm fnm[], t_coupl_rec *tcr, real t,
 568                  int step, real ener[], t_forcerec *fr, t_inputrec *ir,
 569                  gmx_bool bMaster, t_mdatoms *md, t_idef *idef, real mu_aver, int nmols,
 570                  t_commrec *cr, matrix box, tensor virial,
 571                  tensor pres, rvec mu_tot,
 572                  rvec x[], rvec f[], gmx_bool bDoIt)
 573 {
 574 #define enm2Debye 48.0321
 575 #define d2e(x) (x)/enm2Debye
 576 #define enm2kjmol(x) (x)*0.0143952 /* = 2.0*4.0*M_PI*EPSILON0 */
 577
 578     static real     *f6, *f12, *fa, *fb, *fc, *fq;
 579     static gmx_bool  bFirst = TRUE;
 580
 581     int              i, j, ati, atj, atnr2, type, ftype;
 582     real             deviation[eoObsNR], prdev[eoObsNR], epot0, dist, rmsf;
 583     real             ff6, ff12, ffa, ffb, ffc, ffq, factor, dt, mu_ind;
 584     real             Epol, Eintern, Virial, muabs, xiH = -1, xiS = -1, xi6, xi12;
 585     rvec             fmol[2];
 586     gmx_bool         bTest, bPrint;
 587     t_coupl_LJ      *tclj;
 588     t_coupl_BU      *tcbu;
 589     t_coupl_Q       *tcq;
 590     t_coupl_iparams *tip;
 591
 592     atnr2 = idef->atnr * idef->atnr;
 593     if (bFirst)
 594     {
 595         if (PAR(cr))
 596         {
 597             fprintf(log, "GCT: this is parallel\n");
 598         }
 599         else
 600         {
 601             fprintf(log, "GCT: this is not parallel\n");
 602         }
 603         fflush(log);
 604         snew(f6, atnr2);
 605         snew(f12, atnr2);
 606         snew(fa, atnr2);
 607         snew(fb, atnr2);
 608         snew(fc, atnr2);
 609         snew(fq, idef->atnr);
 610
 611         if (tcr->bVirial)
 612         {
 613             int  nrdf = 0;
 614             real TTT  = 0;
 615             real Vol  = det(box);
 616
 617             for (i = 0; (i < ir->opts.ngtc); i++)
 618             {
 619                 nrdf += ir->opts.nrdf[i];
 620                 TTT  += ir->opts.nrdf[i]*ir->opts.ref_t[i];
 621             }
 622             TTT /= nrdf;
 623
 624             /* Calculate reference virial from reference temperature and pressure */
 625             tcr->ref_value[eoVir] = 0.5*BOLTZ*nrdf*TTT - (3.0/2.0)*
 626                 Vol*tcr->ref_value[eoPres];
 627
 628             fprintf(log, "GCT: TTT = %g, nrdf = %d, vir0 = %g,  Vol = %g\n",
 629                     TTT, nrdf, tcr->ref_value[eoVir], Vol);
 630             fflush(log);
 631         }
 632         bFirst = FALSE;
 633     }
 634
 635     bPrint = MASTER(cr) && do_per_step(step, ir->nstlog);
 636     dt     = ir->delta_t;
 637
 638     /* Initiate coupling to the reference pressure and temperature to start
 639      * coupling slowly.
 640      */
 641     if (step == 0)
 642     {
 643         for (i = 0; (i < eoObsNR); i++)
 644         {
 645             tcr->av_value[i] = tcr->ref_value[i];
 646         }
 647         if ((tcr->ref_value[eoDipole]) != 0.0)
 648         {
 649             mu_ind                 = mu_aver - d2e(tcr->ref_value[eoDipole]); /* in e nm */
 650             Epol                   = mu_ind*mu_ind/(enm2kjmol(tcr->ref_value[eoPolarizability]));
 651             tcr->av_value[eoEpot] -= Epol;
 652             fprintf(log, "GCT: mu_aver = %g(D), mu_ind = %g(D), Epol = %g (kJ/mol)\n",
 653                     mu_aver*enm2Debye, mu_ind*enm2Debye, Epol);
 654         }
 655     }
 656
 657     /* We want to optimize the LJ params, usually to the Vaporization energy
 658      * therefore we only count intermolecular degrees of freedom.
 659      * Note that this is now optional. switch UseEinter to yes in your gct file
 660      * if you want this.
 661      */
 662     dist      = calc_dist(log, x);
 663     muabs     = norm(mu_tot);
 664     Eintern   = Ecouple(tcr, ener)/nmols;
 665     Virial    = virial[XX][XX]+virial[YY][YY]+virial[ZZ][ZZ];
 666
 667     /*calc_force(md->nr,f,fmol);*/
 668     clear_rvec(fmol[0]);
 669
 670     /* Use a memory of tcr->nmemory steps, so we actually couple to the
 671      * average observable over the last tcr->nmemory steps. This may help
 672      * in avoiding local minima in parameter space.
 673      */
 674     set_act_value(tcr, eoPres, ener[F_PRES], step);
 675     set_act_value(tcr, eoEpot, Eintern,     step);
 676     set_act_value(tcr, eoVir,  Virial,      step);
 677     set_act_value(tcr, eoDist, dist,        step);
 678     set_act_value(tcr, eoMu,   muabs,       step);
 679     set_act_value(tcr, eoFx,   fmol[0][XX], step);
 680     set_act_value(tcr, eoFy,   fmol[0][YY], step);
 681     set_act_value(tcr, eoFz,   fmol[0][ZZ], step);
 682     set_act_value(tcr, eoPx,   pres[XX][XX], step);
 683     set_act_value(tcr, eoPy,   pres[YY][YY], step);
 684     set_act_value(tcr, eoPz,   pres[ZZ][ZZ], step);
 685
 686     epot0 = tcr->ref_value[eoEpot];
 687     /* If dipole != 0.0 assume we want to use polarization corrected coupling */
 688     if ((tcr->ref_value[eoDipole]) != 0.0)
 689     {
 690         mu_ind = mu_aver - d2e(tcr->ref_value[eoDipole]); /* in e nm */
 691
 692         Epol = mu_ind*mu_ind/(enm2kjmol(tcr->ref_value[eoPolarizability]));
 693
 694         epot0 -= Epol;
 695         if (debug)
 696         {
 697             fprintf(debug, "mu_ind: %g (%g D) mu_aver: %g (%g D)\n",
 698                     mu_ind, mu_ind*enm2Debye, mu_aver, mu_aver*enm2Debye);
 699             fprintf(debug, "Eref %g Epol %g Erunav %g Eact %g\n",
 700                     tcr->ref_value[eoEpot], Epol, tcr->av_value[eoEpot],
 701                     tcr->act_value[eoEpot]);
 702         }
 703     }
 704
 705     if (bPrint)
 706     {
 707         pr_ff(tcr, t, idef, cr, nfile, fnm, oenv);
 708     }
 709     /* Calculate the deviation of average value from the target value */
 710     for (i = 0; (i < eoObsNR); i++)
 711     {
 712         deviation[i] = calc_deviation(tcr->av_value[i], tcr->act_value[i],
 713                                       tcr->ref_value[i]);
 714         prdev[i]     = tcr->ref_value[i] - tcr->act_value[i];
 715     }
 716     deviation[eoEpot] = calc_deviation(tcr->av_value[eoEpot], tcr->act_value[eoEpot],
 717                                        epot0);
 718     prdev[eoEpot]     = epot0 - tcr->act_value[eoEpot];
 719
 720     if (bPrint)
 721     {
 722         pr_dev(tcr, t, prdev, cr, nfile, fnm, oenv);
 723     }
 724
 725     /* First set all factors to 1 */
 726     for (i = 0; (i < atnr2); i++)
 727     {
 728         f6[i] = f12[i] = fa[i] = fb[i] = fc[i] = 1.0;
 729     }
 730     for (i = 0; (i < idef->atnr); i++)
 731     {
 732         fq[i] = 1.0;
 733     }
 734
 735     /* Now compute the actual coupling compononents */
 736     if (!fr->bBHAM)
 737     {
 738         if (bDoIt)
 739         {
 740             for (i = 0; (i < tcr->nLJ); i++)
 741             {
 742                 tclj = &(tcr->tcLJ[i]);
 743
 744                 ati = tclj->at_i;
 745                 atj = tclj->at_j;
 746
 747                 ff6 = ff12 = 1.0;
 748
 749                 if (tclj->eObs == eoForce)
 750                 {
 751                     gmx_fatal(FARGS, "Hack code for this to work again ");
 752                     if (debug)
 753                     {
 754                         fprintf(debug, "Have computed derivatives: xiH = %g, xiS = %g\n", xiH, xiS);
 755                     }
 756                     if (ati == 1)
 757                     {
 758                         /* Hydrogen */
 759                         ff12 += xiH;
 760                     }
 761                     else if (ati == 2)
 762                     {
 763                         /* Shell */
 764                         ff12 += xiS;
 765                     }
 766                     else
 767                     {
 768                         gmx_fatal(FARGS, "No H, no Shell, edit code at %s, line %d\n",
 769                                   __FILE__, __LINE__);
 770                     }
 771                     if (ff6 > 0)
 772                     {
 773                         set_factor_matrix(idef->atnr, f6, sqrt(ff6), ati, atj);
 774                     }
 775                     if (ff12 > 0)
 776                     {
 777                         set_factor_matrix(idef->atnr, f12, sqrt(ff12), ati, atj);
 778                     }
 779                 }
 780                 else
 781                 {
 782                     if (debug)
 783                     {
 784                         fprintf(debug, "tcr->tcLJ[%d].xi_6 = %g, xi_12 = %g deviation = %g\n", i,
 785                                 tclj->xi_6, tclj->xi_12, deviation[tclj->eObs]);
 786                     }
 787                     factor = deviation[tclj->eObs];
 788
 789                     upd_f_value(log, idef->atnr, tclj->xi_6, dt, factor, f6, ati, atj);
 790                     upd_f_value(log, idef->atnr, tclj->xi_12, dt, factor, f12, ati, atj);
 791                 }
 792             }
 793         }
 794         if (PAR(cr))
 795         {
 796             gprod(cr, atnr2, f6);
 797             gprod(cr, atnr2, f12);
 798 #ifdef DEBUGGCT
 799             dump_fm(log, idef->atnr, f6, "f6");
 800             dump_fm(log, idef->atnr, f12, "f12");
 801 #endif
 802         }
 803         upd_nbfplj(log, fr->nbfp, idef->atnr, f6, f12, tcr->combrule);
 804
 805         /* Copy for printing */
 806         for (i = 0; (i < tcr->nLJ); i++)
 807         {
 808             tclj = &(tcr->tcLJ[i]);
 809             ati  = tclj->at_i;
 810             atj  = tclj->at_j;
 811             if (atj == -1)
 812             {
 813                 atj = ati;
 814             }
 815             /* nbfp now includes the 6.0/12.0 derivative prefactors */
 816             tclj->c6  =  C6(fr->nbfp, fr->ntype, ati, atj)/6.0;
 817             tclj->c12 = C12(fr->nbfp, fr->ntype, ati, atj)/12.0;
 818         }
 819     }
 820     else
 821     {
 822         if (bDoIt)
 823         {
 824             for (i = 0; (i < tcr->nBU); i++)
 825             {
 826                 tcbu   = &(tcr->tcBU[i]);
 827                 factor = deviation[tcbu->eObs];
 828                 ati    = tcbu->at_i;
 829                 atj    = tcbu->at_j;
 830
 831                 upd_f_value(log, idef->atnr, tcbu->xi_a, dt, factor, fa, ati, atj);
 832                 upd_f_value(log, idef->atnr, tcbu->xi_b, dt, factor, fb, ati, atj);
 833                 upd_f_value(log, idef->atnr, tcbu->xi_c, dt, factor, fc, ati, atj);
 834             }
 835         }
 836         if (PAR(cr))
 837         {
 838             gprod(cr, atnr2, fa);
 839             gprod(cr, atnr2, fb);
 840             gprod(cr, atnr2, fc);
 841         }
 842         upd_nbfpbu(log, fr->nbfp, idef->atnr, fa, fb, fc);
 843         /* Copy for printing */
 844         for (i = 0; (i < tcr->nBU); i++)
 845         {
 846             tcbu = &(tcr->tcBU[i]);
 847             ati  = tcbu->at_i;
 848             atj  = tcbu->at_j;
 849             if (atj == -1)
 850             {
 851                 atj = ati;
 852             }
 853             /* nbfp now includes the 6.0 derivative prefactors */
 854             tcbu->a = BHAMA(fr->nbfp, fr->ntype, ati, atj);
 855             tcbu->b = BHAMB(fr->nbfp, fr->ntype, ati, atj);
 856             tcbu->c = BHAMC(fr->nbfp, fr->ntype, ati, atj)/6.0;
 857             if (debug)
 858             {
 859                 fprintf(debug, "buck (type=%d) = %e, %e, %e\n",
 860                         tcbu->at_i, tcbu->a, tcbu->b, tcbu->c);
 861             }
 862         }
 863     }
 864     if (bDoIt)
 865     {
 866         for (i = 0; (i < tcr->nQ); i++)
 867         {
 868             tcq = &(tcr->tcQ[i]);
 869             if (tcq->xi_Q)
 870             {
 871                 ffq = 1.0 + (dt/tcq->xi_Q) * deviation[tcq->eObs];
 872             }
 873             else
 874             {
 875                 ffq = 1.0;
 876             }
 877             fq[tcq->at_i] *= ffq;
 878         }
 879     }
 880     if (PAR(cr))
 881     {
 882         gprod(cr, idef->atnr, fq);
 883     }
 884
 885     for (j = 0; (j < md->nr); j++)
 886     {
 887         md->chargeA[j] *= fq[md->typeA[j]];
 888     }
 889     for (i = 0; (i < tcr->nQ); i++)
 890     {
 891         tcq = &(tcr->tcQ[i]);
 892         for (j = 0; (j < md->nr); j++)
 893         {
 894             if (md->typeA[j] == tcq->at_i)
 895             {
 896                 tcq->Q = md->chargeA[j];
 897                 break;
 898             }
 899         }
 900         if (j == md->nr)
 901         {
 902             gmx_fatal(FARGS, "Coupling type %d not found", tcq->at_i);
 903         }
 904     }
 905     for (i = 0; (i < tcr->nIP); i++)
 906     {
 907         tip    = &(tcr->tIP[i]);
 908         type   = tip->type;
 909         ftype  = idef->functype[type];
 910         factor = dt*deviation[tip->eObs];
 911
 912         switch (ftype)
 913         {
 914             case F_BONDS:
 915                 if (tip->xi.harmonic.krA)
 916                 {
 917                     idef->iparams[type].harmonic.krA *= (1+factor/tip->xi.harmonic.krA);
 918                 }
 919                 if (tip->xi.harmonic.rA)
 920                 {
 921                     idef->iparams[type].harmonic.rA *= (1+factor/tip->xi.harmonic.rA);
 922                 }
 923                 break;
 924             default:
 925                 break;
 926         }
 927         tip->iprint = idef->iparams[type];
 928     }
 929 }