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