src/kernel/do_gct.c

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