src/kernel/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     read_gct (opt2fn("-j",nfile,fnm), tcr);
  68     write_gct(opt2fn("-jo",nfile,fnm),tcr,idef);
  69   }
  70   /* Update all processors with coupling info */
  71   if (PAR(cr))
  72     comm_tcr(log,cr,&tcr);
  73
  74   /* Copy information from the coupling to the force field stuff */
  75   copy_ff(tcr,fr,md,idef);
  76
  77   return tcr;
  78 }
  79
  80 static real Ecouple(t_coupl_rec *tcr,real ener[])
  81 {
  82   if (tcr->bInter)
  83     return ener[F_COUL_SR]+ener[F_LJ]+ener[F_COUL_LR]+ener[F_LJ_LR]+
  84       ener[F_RF_EXCL]+ener[F_COUL_RECIP];
  85   else
  86     return ener[F_EPOT];
  87 }
  88
  89 static char *mk_gct_nm(const char *fn,int ftp,int ati,int atj)
  90 {
  91   static char buf[256];
  92
  93   strcpy(buf,fn);
  94   if (atj == -1)
  95     sprintf(buf+strlen(fn)-4,"%d.%s",ati,ftp2ext(ftp));
  96   else
  97     sprintf(buf+strlen(fn)-4,"%d_%d.%s",ati,atj,ftp2ext(ftp));
  98
  99   return buf;
 100 }
 101
 102 static void pr_ff(t_coupl_rec *tcr,real time,t_idef *idef,
 103                   t_commrec *cr,int nfile,const t_filenm fnm[],
 104                   const output_env_t oenv)
 105 {
 106   static FILE *prop;
 107   static FILE **out=NULL;
 108   static FILE **qq=NULL;
 109   static FILE **ip=NULL;
 110   t_coupl_LJ  *tclj;
 111   t_coupl_BU  *tcbu;
 112   char        buf[256];
 113   const char *leg[]  =  { "C12", "C6" };
 114   const char *eleg[] =  { "Epsilon", "Sigma" };
 115   const char *bleg[] = { "A", "B", "C" };
 116   char        **raleg;
 117   int         i,j,index;
 118
 119   if ((prop == NULL) && (out == NULL) && (qq == NULL) && (ip == NULL)) {
 120     prop=xvgropen(opt2fn("-runav",nfile,fnm),
 121                   "Properties and Running Averages","Time (ps)","",oenv);
 122     snew(raleg,2*eoObsNR);
 123     for(i=j=0; (i<eoObsNR); i++) {
 124       if (tcr->bObsUsed[i]) {
 125         raleg[j++] = strdup(eoNames[i]);
 126         sprintf(buf,"RA-%s",eoNames[i]);
 127         raleg[j++] = strdup(buf);
 128       }
 129     }
 130     xvgr_legend(prop,j,(const char**)raleg,oenv);
 131     for(i=0; (i<j); i++)
 132       sfree(raleg[i]);
 133     sfree(raleg);
 134
 135     if (tcr->nLJ) {
 136       snew(out,tcr->nLJ);
 137       for(i=0; (i<tcr->nLJ); i++) {
 138         if (tcr->tcLJ[i].bPrint) {
 139           tclj   = &(tcr->tcLJ[i]);
 140           out[i] =
 141             xvgropen(mk_gct_nm(opt2fn("-ffout",nfile,fnm),
 142                                efXVG,tclj->at_i,tclj->at_j),
 143                      "General Coupling Lennard Jones","Time (ps)",
 144                      "Force constant (units)",oenv);
 145           fprintf(out[i],"@ subtitle \"Interaction between types %d and %d\"\n",
 146                   tclj->at_i,tclj->at_j);
 147           if (tcr->combrule == 1)
 148             xvgr_legend(out[i],asize(leg),leg,oenv);
 149           else
 150             xvgr_legend(out[i],asize(eleg),eleg,oenv);
 151           fflush(out[i]);
 152         }
 153       }
 154     }
 155     else if (tcr->nBU) {
 156       snew(out,tcr->nBU);
 157       for(i=0; (i<tcr->nBU); i++) {
 158         if (tcr->tcBU[i].bPrint) {
 159           tcbu=&(tcr->tcBU[i]);
 160           out[i] =
 161             xvgropen(mk_gct_nm(opt2fn("-ffout",nfile,fnm),efXVG,
 162                                tcbu->at_i,tcbu->at_j),
 163                      "General Coupling Buckingham","Time (ps)",
 164                      "Force constant (units)",oenv);
 165           fprintf(out[i],"@ subtitle \"Interaction between types %d and %d\"\n",
 166                   tcbu->at_i,tcbu->at_j);
 167           xvgr_legend(out[i],asize(bleg),bleg,oenv);
 168           fflush(out[i]);
 169         }
 170       }
 171     }
 172     snew(qq,tcr->nQ);
 173     for(i=0; (i<tcr->nQ); i++) {
 174       if (tcr->tcQ[i].bPrint) {
 175         qq[i] = xvgropen(mk_gct_nm(opt2fn("-ffout",nfile,fnm),efXVG,
 176                                    tcr->tcQ[i].at_i,-1),
 177                          "General Coupling Charge","Time (ps)","Charge (e)",
 178                          oenv);
 179         fprintf(qq[i],"@ subtitle \"Type %d\"\n",tcr->tcQ[i].at_i);
 180         fflush(qq[i]);
 181       }
 182     }
 183     snew(ip,tcr->nIP);
 184     for(i=0; (i<tcr->nIP); i++) {
 185       sprintf(buf,"gctIP%d",tcr->tIP[i].type);
 186       ip[i]=xvgropen(mk_gct_nm(opt2fn("-ffout",nfile,fnm),efXVG,0,-1),
 187                      "General Coupling iparams","Time (ps)","ip ()",oenv);
 188       index=tcr->tIP[i].type;
 189       fprintf(ip[i],"@ subtitle \"Coupling to %s\"\n",
 190               interaction_function[idef->functype[index]].longname);
 191       fflush(ip[i]);
 192     }
 193   }
 194   /* Write properties to file */
 195   fprintf(prop,"%10.3f",time);
 196   for(i=0; (i<eoObsNR); i++)
 197     if (tcr->bObsUsed[i])
 198       fprintf(prop,"  %10.3e  %10.3e",tcr->act_value[i],tcr->av_value[i]);
 199   fprintf(prop,"\n");
 200   fflush(prop);
 201
 202   for(i=0; (i<tcr->nLJ); i++) {
 203     tclj=&(tcr->tcLJ[i]);
 204     if (tclj->bPrint) {
 205       if (tcr->combrule == 1)
 206         fprintf(out[i],"%14.7e  %14.7e  %14.7e\n",
 207                 time,tclj->c12,tclj->c6);
 208       else {
 209         double sigma   = pow(tclj->c12/tclj->c6,1.0/6.0);
 210         double epsilon = 0.25*sqr(tclj->c6)/tclj->c12;
 211         fprintf(out[i],"%14.7e  %14.7e  %14.7e\n",
 212                 time,epsilon,sigma);
 213       }
 214       fflush(out[i]);
 215     }
 216   }
 217   for(i=0; (i<tcr->nBU); i++) {
 218     tcbu=&(tcr->tcBU[i]);
 219     if (tcbu->bPrint) {
 220       fprintf(out[i],"%14.7e  %14.7e  %14.7e  %14.7e\n",
 221               time,tcbu->a,tcbu->b,tcbu->c);
 222       fflush(out[i]);
 223     }
 224   }
 225   for(i=0; (i<tcr->nQ); i++) {
 226     if (tcr->tcQ[i].bPrint) {
 227       fprintf(qq[i],"%14.7e  %14.7e\n",time,tcr->tcQ[i].Q);
 228       fflush(qq[i]);
 229     }
 230   }
 231   for(i=0; (i<tcr->nIP); i++) {
 232     fprintf(ip[i],"%10g  ",time);
 233     index=tcr->tIP[i].type;
 234     switch(idef->functype[index]) {
 235     case F_BONDS:
 236       fprintf(ip[i],"%10g  %10g\n",tcr->tIP[i].iprint.harmonic.krA,
 237               tcr->tIP[i].iprint.harmonic.rA);
 238       break;
 239     default:
 240       break;
 241     }
 242     fflush(ip[i]);
 243   }
 244 }
 245
 246 static void pr_dev(t_coupl_rec *tcr,
 247                    real t,real dev[eoObsNR],t_commrec *cr,int nfile,
 248                    const t_filenm fnm[],const output_env_t oenv)
 249 {
 250   static FILE *fp=NULL;
 251   char   **ptr;
 252   int    i,j;
 253
 254   if (!fp) {
 255     fp=xvgropen(opt2fn("-devout",nfile,fnm),
 256                 "Deviations from target value","Time (ps)","",oenv);
 257     snew(ptr,eoObsNR);
 258     for(i=j=0; (i<eoObsNR); i++)
 259       if (tcr->bObsUsed[i])
 260                   ptr[j++] = strdup(eoNames[i]);
 261           xvgr_legend(fp,j,(const char**)ptr,oenv);
 262           for(i=0;i<j;i++)
 263                   sfree(ptr[i]);
 264     sfree(ptr);
 265   }
 266   fprintf(fp,"%10.3f",t);
 267   for(i=0; (i<eoObsNR); i++)
 268     if (tcr->bObsUsed[i])
 269       fprintf(fp,"  %10.3e",dev[i]);
 270   fprintf(fp,"\n");
 271   fflush(fp);
 272 }
 273
 274 static void upd_nbfplj(FILE *log,real *nbfp,int atnr,real f6[],real f12[],
 275                        int combrule)
 276 {
 277   double *sigma,*epsilon,c6,c12,eps,sig,sig6;
 278   int n,m,k;
 279
 280   /* Update the nonbonded force parameters */
 281   switch (combrule) {
 282   case 1:
 283     for(k=n=0; (n<atnr); n++) {
 284       for(m=0; (m<atnr); m++,k++) {
 285         C6 (nbfp,atnr,n,m) *= f6[k];
 286         C12(nbfp,atnr,n,m) *= f12[k];
 287       }
 288     }
 289     break;
 290   case 2:
 291   case 3:
 292     /* Convert to sigma and epsilon */
 293     snew(sigma,atnr);
 294     snew(epsilon,atnr);
 295     for(n=0; (n<atnr); n++) {
 296       k = n*(atnr+1);
 297       c6  = C6 (nbfp,atnr,n,n) * f6[k];
 298       c12 = C12(nbfp,atnr,n,n) * f12[k];
 299       if ((c6 == 0) || (c12 == 0))
 300         gmx_fatal(FARGS,"You can not use combination rule %d with zero C6 (%f) or C12 (%f)",combrule,c6,c12);
 301       sigma[n]   = pow(c12/c6,1.0/6.0);
 302       epsilon[n] = 0.25*(c6*c6/c12);
 303     }
 304     for(k=n=0; (n<atnr); n++) {
 305       for(m=0; (m<atnr); m++,k++) {
 306         eps  = sqrt(epsilon[n]*epsilon[m]);
 307         if (combrule == 2)
 308           sig  = 0.5*(sigma[n]+sigma[m]);
 309         else
 310           sig  = sqrt(sigma[n]*sigma[m]);
 311         sig6 = pow(sig,6.0);
 312     /* nbfp now includes the 6.0/12.0 derivative prefactors */
 313         C6 (nbfp,atnr,n,m) = 4*eps*sig6/6.0;
 314         C12(nbfp,atnr,n,m) = 4*eps*sig6*sig6/12.0;
 315       }
 316     }
 317     sfree(sigma);
 318     sfree(epsilon);
 319     break;
 320   default:
 321     gmx_fatal(FARGS,"Combination rule should be 1,2 or 3 instead of %d",
 322               combrule);
 323   }
 324 }
 325
 326 static void upd_nbfpbu(FILE *log,real *nbfp,int atnr,
 327                        real fa[],real fb[],real fc[])
 328 {
 329   int n,m,k;
 330
 331   /* Update the nonbonded force parameters */
 332   for(k=n=0; (n<atnr); n++) {
 333     for(m=0; (m<atnr); m++,k++) {
 334       BHAMA(nbfp,atnr,n,m) *= fa[k];
 335       BHAMB(nbfp,atnr,n,m) *= fb[k];
 336       BHAMC(nbfp,atnr,n,m) *= fc[k];
 337     }
 338   }
 339 }
 340
 341 void gprod(t_commrec *cr,int n,real f[])
 342 {
 343   /* Compute the global product of all elements in an array
 344    * such that after gprod f[i] = PROD_j=1,nprocs f[i][j]
 345    */
 346   static int  nbuf=0;
 347   static real *buf=NULL;
 348   int  i;
 349
 350   if (n > nbuf) {
 351     nbuf = n;
 352     srenew(buf, nbuf);
 353   }
 354
 355 #ifdef GMX_MPI
 356 #ifdef GMX_DOUBLE
 357   MPI_Allreduce(f,buf,n,MPI_DOUBLE,MPI_PROD,cr->mpi_comm_mygroup);
 358 #else
 359   MPI_Allreduce(f,buf,n,MPI_FLOAT, MPI_PROD,cr->mpi_comm_mygroup);
 360 #endif
 361   for(i=0; (i<n); i++)
 362     f[i] = buf[i];
 363 #endif
 364 }
 365
 366 static void set_factor_matrix(int ntypes,real f[],real fmult,int ati,int atj)
 367 {
 368 #define FMIN 0.95
 369 #define FMAX 1.05
 370   int i;
 371
 372   fmult = min(FMAX,max(FMIN,fmult));
 373   if (atj != -1) {
 374     f[ntypes*ati+atj] *= fmult;
 375     f[ntypes*atj+ati] *= fmult;
 376   }
 377   else {
 378     for(i=0; (i<ntypes); i++) {
 379       f[ntypes*ati+i] *= fmult;
 380       f[ntypes*i+ati] *= fmult;
 381     }
 382   }
 383 #undef FMIN
 384 #undef FMAX
 385 }
 386
 387 static real calc_deviation(real xav,real xt,real x0)
 388 {
 389   /* This may prevent overshooting in GCT coupling... */
 390
 391   /* real dev;
 392
 393   if (xav > x0) {
 394     if (xt > x0)
 395       dev = min(xav-x0,xt-x0);
 396     else
 397       dev = 0;
 398   }
 399   else {
 400     if (xt < x0)
 401       dev = max(xav-x0,xt-x0);
 402     else
 403       dev = 0;
 404   }
 405 */
 406   return x0-xav;
 407 }
 408
 409 static real calc_dist(FILE *log,rvec x[])
 410 {
 411   static gmx_bool bFirst=TRUE;
 412   static gmx_bool bDist;
 413   static int i1,i2;
 414   char *buf;
 415   rvec   dx;
 416
 417   if (bFirst) {
 418     if ((buf = getenv("DISTGCT")) == NULL)
 419       bDist = FALSE;
 420     else {
 421       bDist  = (sscanf(buf,"%d%d",&i1,&i2) == 2);
 422       if (bDist)
 423         fprintf(log,"GCT: Will couple to distance between %d and %d\n",i1,i2);
 424       else
 425         fprintf(log,"GCT: Will not couple to distances\n");
 426     }
 427     bFirst = FALSE;
 428   }
 429   if (bDist) {
 430     rvec_sub(x[i1],x[i2],dx);
 431     return norm(dx);
 432   }
 433   else
 434     return 0.0;
 435 }
 436
 437 real run_aver(real old,real cur,int step,int nmem)
 438 {
 439   nmem   = max(1,nmem);
 440
 441   return ((nmem-1)*old+cur)/nmem;
 442 }
 443
 444 static void set_act_value(t_coupl_rec *tcr,int index,real val,int step)
 445 {
 446   tcr->act_value[index] = val;
 447   tcr->av_value[index]  = run_aver(tcr->av_value[index],val,step,tcr->nmemory);
 448 }
 449
 450 static void upd_f_value(FILE *log,int atnr,real xi,real dt,real factor,
 451                         real ff[],int ati,int atj)
 452 {
 453   real fff;
 454
 455   if (xi != 0) {
 456     fff = 1 + (dt/xi)  * factor;
 457     if (fff > 0)
 458       set_factor_matrix(atnr,ff,sqrt(fff),ati,atj);
 459   }
 460 }
 461
 462 static void dump_fm(FILE *fp,int n,real f[],char *s)
 463 {
 464   int i,j;
 465
 466   fprintf(fp,"Factor matrix (all numbers -1) %s\n",s);
 467   for(i=0; (i<n); i++) {
 468     for(j=0; (j<n); j++)
 469       fprintf(fp,"  %10.3e",f[n*i+j]-1.0);
 470     fprintf(fp,"\n");
 471   }
 472 }
 473
 474 void do_coupling(FILE *log,const output_env_t oenv,int nfile,
 475                  const t_filenm fnm[], t_coupl_rec *tcr,real t,
 476                  int step,real ener[], t_forcerec *fr,t_inputrec *ir,
 477                  gmx_bool bMaster, t_mdatoms *md,t_idef *idef,real mu_aver,int nmols,
 478                  t_commrec *cr,matrix box,tensor virial,
 479                  tensor pres,rvec mu_tot,
 480                  rvec x[],rvec f[],gmx_bool bDoIt)
 481 {
 482 #define enm2Debye 48.0321
 483 #define d2e(x) (x)/enm2Debye
 484 #define enm2kjmol(x) (x)*0.0143952 /* = 2.0*4.0*M_PI*EPSILON0 */
 485
 486   static real *f6,*f12,*fa,*fb,*fc,*fq;
 487   static gmx_bool bFirst = TRUE;
 488
 489   int         i,j,ati,atj,atnr2,type,ftype;
 490   real        deviation[eoObsNR],prdev[eoObsNR],epot0,dist,rmsf;
 491   real        ff6,ff12,ffa,ffb,ffc,ffq,factor,dt,mu_ind;
 492   real        Epol,Eintern,Virial,muabs,xiH=-1,xiS=-1,xi6,xi12;
 493   rvec        fmol[2];
 494   gmx_bool        bTest,bPrint;
 495   t_coupl_LJ  *tclj;
 496   t_coupl_BU  *tcbu;
 497   t_coupl_Q   *tcq;
 498   t_coupl_iparams *tip;
 499
 500   atnr2 = idef->atnr * idef->atnr;
 501   if (bFirst) {
 502     if (PAR(cr))
 503       fprintf(log,"GCT: this is parallel\n");
 504     else
 505       fprintf(log,"GCT: this is not parallel\n");
 506     fflush(log);
 507     snew(f6, atnr2);
 508     snew(f12,atnr2);
 509     snew(fa, atnr2);
 510     snew(fb, atnr2);
 511     snew(fc, atnr2);
 512     snew(fq, idef->atnr);
 513
 514     if (tcr->bVirial) {
 515       int  nrdf = 0;
 516       real TTT  = 0;
 517       real Vol  = det(box);
 518
 519       for(i=0; (i<ir->opts.ngtc); i++) {
 520         nrdf += ir->opts.nrdf[i];
 521         TTT  += ir->opts.nrdf[i]*ir->opts.ref_t[i];
 522       }
 523       TTT /= nrdf;
 524
 525       /* Calculate reference virial from reference temperature and pressure */
 526       tcr->ref_value[eoVir] = 0.5*BOLTZ*nrdf*TTT - (3.0/2.0)*
 527         Vol*tcr->ref_value[eoPres];
 528
 529       fprintf(log,"GCT: TTT = %g, nrdf = %d, vir0 = %g,  Vol = %g\n",
 530               TTT,nrdf,tcr->ref_value[eoVir],Vol);
 531       fflush(log);
 532     }
 533     bFirst = FALSE;
 534   }
 535
 536   bPrint = MASTER(cr) && do_per_step(step,ir->nstlog);
 537   dt     = ir->delta_t;
 538
 539   /* Initiate coupling to the reference pressure and temperature to start
 540    * coupling slowly.
 541    */
 542   if (step == 0) {
 543     for(i=0; (i<eoObsNR); i++)
 544       tcr->av_value[i] = tcr->ref_value[i];
 545     if ((tcr->ref_value[eoDipole]) != 0.0) {
 546       mu_ind = mu_aver - d2e(tcr->ref_value[eoDipole]); /* in e nm */
 547       Epol   = mu_ind*mu_ind/(enm2kjmol(tcr->ref_value[eoPolarizability]));
 548       tcr->av_value[eoEpot] -= Epol;
 549       fprintf(log,"GCT: mu_aver = %g(D), mu_ind = %g(D), Epol = %g (kJ/mol)\n",
 550               mu_aver*enm2Debye,mu_ind*enm2Debye,Epol);
 551     }
 552   }
 553
 554   /* We want to optimize the LJ params, usually to the Vaporization energy
 555    * therefore we only count intermolecular degrees of freedom.
 556    * Note that this is now optional. switch UseEinter to yes in your gct file
 557    * if you want this.
 558    */
 559   dist      = calc_dist(log,x);
 560   muabs     = norm(mu_tot);
 561   Eintern   = Ecouple(tcr,ener)/nmols;
 562   Virial    = virial[XX][XX]+virial[YY][YY]+virial[ZZ][ZZ];
 563
 564   /*calc_force(md->nr,f,fmol);*/
 565   clear_rvec(fmol[0]);
 566
 567   /* Use a memory of tcr->nmemory steps, so we actually couple to the
 568    * average observable over the last tcr->nmemory steps. This may help
 569    * in avoiding local minima in parameter space.
 570    */
 571   set_act_value(tcr,eoPres, ener[F_PRES],step);
 572   set_act_value(tcr,eoEpot, Eintern,     step);
 573   set_act_value(tcr,eoVir,  Virial,      step);
 574   set_act_value(tcr,eoDist, dist,        step);
 575   set_act_value(tcr,eoMu,   muabs,       step);
 576   set_act_value(tcr,eoFx,   fmol[0][XX], step);
 577   set_act_value(tcr,eoFy,   fmol[0][YY], step);
 578   set_act_value(tcr,eoFz,   fmol[0][ZZ], step);
 579   set_act_value(tcr,eoPx,   pres[XX][XX],step);
 580   set_act_value(tcr,eoPy,   pres[YY][YY],step);
 581   set_act_value(tcr,eoPz,   pres[ZZ][ZZ],step);
 582
 583   epot0 = tcr->ref_value[eoEpot];
 584   /* If dipole != 0.0 assume we want to use polarization corrected coupling */
 585   if ((tcr->ref_value[eoDipole]) != 0.0) {
 586     mu_ind = mu_aver - d2e(tcr->ref_value[eoDipole]); /* in e nm */
 587
 588     Epol = mu_ind*mu_ind/(enm2kjmol(tcr->ref_value[eoPolarizability]));
 589
 590     epot0 -= Epol;
 591     if (debug) {
 592       fprintf(debug,"mu_ind: %g (%g D) mu_aver: %g (%g D)\n",
 593               mu_ind,mu_ind*enm2Debye,mu_aver,mu_aver*enm2Debye);
 594       fprintf(debug,"Eref %g Epol %g Erunav %g Eact %g\n",
 595               tcr->ref_value[eoEpot],Epol,tcr->av_value[eoEpot],
 596               tcr->act_value[eoEpot]);
 597     }
 598   }
 599
 600   if (bPrint) {
 601     pr_ff(tcr,t,idef,cr,nfile,fnm,oenv);
 602   }
 603   /* Calculate the deviation of average value from the target value */
 604   for(i=0; (i<eoObsNR); i++) {
 605     deviation[i] = calc_deviation(tcr->av_value[i],tcr->act_value[i],
 606                                   tcr->ref_value[i]);
 607     prdev[i]     = tcr->ref_value[i] - tcr->act_value[i];
 608   }
 609   deviation[eoEpot] = calc_deviation(tcr->av_value[eoEpot],tcr->act_value[eoEpot],
 610                                      epot0);
 611   prdev[eoEpot]     = epot0 - tcr->act_value[eoEpot];
 612
 613   if (bPrint)
 614     pr_dev(tcr,t,prdev,cr,nfile,fnm,oenv);
 615
 616   /* First set all factors to 1 */
 617   for(i=0; (i<atnr2); i++) {
 618     f6[i] = f12[i] = fa[i] = fb[i] = fc[i] = 1.0;
 619   }
 620   for(i=0; (i<idef->atnr); i++)
 621     fq[i] = 1.0;
 622
 623   /* Now compute the actual coupling compononents */
 624   if (!fr->bBHAM) {
 625     if (bDoIt) {
 626       for(i=0; (i<tcr->nLJ); i++) {
 627         tclj=&(tcr->tcLJ[i]);
 628
 629         ati=tclj->at_i;
 630         atj=tclj->at_j;
 631
 632         ff6 = ff12 = 1.0;
 633
 634         if (tclj->eObs == eoForce) {
 635           gmx_fatal(FARGS,"Hack code for this to work again ");
 636           if (debug)
 637             fprintf(debug,"Have computed derivatives: xiH = %g, xiS = %g\n",xiH,xiS);
 638           if (ati == 1) {
 639             /* Hydrogen */
 640             ff12 += xiH;
 641           }
 642           else if (ati == 2) {
 643             /* Shell */
 644             ff12 += xiS;
 645           }
 646           else
 647             gmx_fatal(FARGS,"No H, no Shell, edit code at %s, line %d\n",
 648                         __FILE__,__LINE__);
 649           if (ff6 > 0)
 650             set_factor_matrix(idef->atnr,f6, sqrt(ff6), ati,atj);
 651           if (ff12 > 0)
 652             set_factor_matrix(idef->atnr,f12,sqrt(ff12),ati,atj);
 653         }
 654         else {
 655           if (debug)
 656             fprintf(debug,"tcr->tcLJ[%d].xi_6 = %g, xi_12 = %g deviation = %g\n",i,
 657                     tclj->xi_6,tclj->xi_12,deviation[tclj->eObs]);
 658           factor=deviation[tclj->eObs];
 659
 660           upd_f_value(log,idef->atnr,tclj->xi_6, dt,factor,f6, ati,atj);
 661           upd_f_value(log,idef->atnr,tclj->xi_12,dt,factor,f12,ati,atj);
 662         }
 663       }
 664     }
 665     if (PAR(cr)) {
 666       gprod(cr,atnr2,f6);
 667       gprod(cr,atnr2,f12);
 668 #ifdef DEBUGGCT
 669       dump_fm(log,idef->atnr,f6,"f6");
 670       dump_fm(log,idef->atnr,f12,"f12");
 671 #endif
 672     }
 673     upd_nbfplj(log,fr->nbfp,idef->atnr,f6,f12,tcr->combrule);
 674
 675     /* Copy for printing */
 676     for(i=0; (i<tcr->nLJ); i++) {
 677       tclj=&(tcr->tcLJ[i]);
 678       ati = tclj->at_i;
 679       atj = tclj->at_j;
 680       if (atj == -1)
 681       {
 682           atj = ati;
 683       }
 684       /* nbfp now includes the 6.0/12.0 derivative prefactors */
 685       tclj->c6  =  C6(fr->nbfp,fr->ntype,ati,atj)/6.0;
 686       tclj->c12 = C12(fr->nbfp,fr->ntype,ati,atj)/12.0;
 687     }
 688   }
 689   else {
 690     if (bDoIt) {
 691       for(i=0; (i<tcr->nBU); i++) {
 692         tcbu   = &(tcr->tcBU[i]);
 693         factor = deviation[tcbu->eObs];
 694         ati    = tcbu->at_i;
 695         atj    = tcbu->at_j;
 696
 697         upd_f_value(log,idef->atnr,tcbu->xi_a,dt,factor,fa,ati,atj);
 698         upd_f_value(log,idef->atnr,tcbu->xi_b,dt,factor,fb,ati,atj);
 699         upd_f_value(log,idef->atnr,tcbu->xi_c,dt,factor,fc,ati,atj);
 700       }
 701     }
 702     if (PAR(cr)) {
 703       gprod(cr,atnr2,fa);
 704       gprod(cr,atnr2,fb);
 705       gprod(cr,atnr2,fc);
 706     }
 707     upd_nbfpbu(log,fr->nbfp,idef->atnr,fa,fb,fc);
 708     /* Copy for printing */
 709     for(i=0; (i<tcr->nBU); i++) {
 710       tcbu=&(tcr->tcBU[i]);
 711       ati = tcbu->at_i;
 712       atj = tcbu->at_j;
 713       if (atj == -1)
 714       {
 715           atj = ati;
 716       }
 717       /* nbfp now includes the 6.0 derivative prefactors */
 718       tcbu->a = BHAMA(fr->nbfp,fr->ntype,ati,atj);
 719       tcbu->b = BHAMB(fr->nbfp,fr->ntype,ati,atj);
 720       tcbu->c = BHAMC(fr->nbfp,fr->ntype,ati,atj)/6.0;
 721       if (debug)
 722         fprintf(debug,"buck (type=%d) = %e, %e, %e\n",
 723                 tcbu->at_i,tcbu->a,tcbu->b,tcbu->c);
 724     }
 725   }
 726   if (bDoIt) {
 727     for(i=0; (i<tcr->nQ); i++) {
 728       tcq=&(tcr->tcQ[i]);
 729       if (tcq->xi_Q)
 730         ffq = 1.0 + (dt/tcq->xi_Q) * deviation[tcq->eObs];
 731       else
 732         ffq = 1.0;
 733       fq[tcq->at_i] *= ffq;
 734     }
 735   }
 736   if (PAR(cr))
 737     gprod(cr,idef->atnr,fq);
 738
 739   for(j=0; (j<md->nr); j++) {
 740     md->chargeA[j] *= fq[md->typeA[j]];
 741   }
 742   for(i=0; (i<tcr->nQ); i++) {
 743     tcq=&(tcr->tcQ[i]);
 744     for(j=0; (j<md->nr); j++) {
 745       if (md->typeA[j] == tcq->at_i) {
 746         tcq->Q = md->chargeA[j];
 747         break;
 748       }
 749     }
 750     if (j == md->nr)
 751       gmx_fatal(FARGS,"Coupling type %d not found",tcq->at_i);
 752   }
 753   for(i=0; (i<tcr->nIP); i++) {
 754     tip    = &(tcr->tIP[i]);
 755     type   = tip->type;
 756     ftype  = idef->functype[type];
 757     factor = dt*deviation[tip->eObs];
 758
 759     switch(ftype) {
 760     case F_BONDS:
 761       if (tip->xi.harmonic.krA) idef->iparams[type].harmonic.krA *= (1+factor/tip->xi.harmonic.krA);
 762       if (tip->xi.harmonic.rA) idef->iparams[type].harmonic.rA *= (1+factor/tip->xi.harmonic.rA);
 763         break;
 764     default:
 765       break;
 766     }
 767     tip->iprint=idef->iparams[type];
 768   }
 769 }
 770