src/contrib/do_gct.c

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