src/gromacs/gmxlib/orires.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   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.2.0
  11  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  12  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  13  * Copyright (c) 2001-2004, The GROMACS development team,
  14  * check out http://www.gromacs.org for more information.
  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  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * GROningen Mixture of Alchemy and Childrens' Stories
  35  */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include "typedefs.h"
  41 #include "smalloc.h"
  42 #include "vec.h"
  43 #include "nrjac.h"
  44 #include "network.h"
  45 #include "orires.h"
  46 #include "do_fit.h"
  47 #include "main.h"
  48 #include "copyrite.h"
  49 #include "pbc.h"
  50 #include "mtop_util.h"
  51
  52 void init_orires(FILE *fplog,const gmx_mtop_t *mtop,
  53                  rvec xref[],
  54                  const t_inputrec *ir,
  55                  const gmx_multisim_t *ms,t_oriresdata *od,
  56                  t_state *state)
  57 {
  58     int    i,j,d,ex,nmol,nr,*nr_ex;
  59     double mtot;
  60     rvec   com;
  61     gmx_mtop_ilistloop_t iloop;
  62     t_ilist *il;
  63     gmx_mtop_atomloop_all_t aloop;
  64     t_atom *atom;
  65
  66     od->fc  = ir->orires_fc;
  67     od->nex = 0;
  68     od->S   = NULL;
  69
  70     od->M=NULL;
  71     od->eig=NULL;
  72     od->v=NULL;
  73
  74     od->nr = gmx_mtop_ftype_count(mtop,F_ORIRES);
  75     if (od->nr == 0)
  76     {
  77         return;
  78     }
  79
  80     nr_ex = NULL;
  81
  82     iloop = gmx_mtop_ilistloop_init(mtop);
  83     while (gmx_mtop_ilistloop_next(iloop,&il,&nmol))
  84     {
  85         for(i=0; i<il[F_ORIRES].nr; i+=3)
  86         {
  87             ex = mtop->ffparams.iparams[il[F_ORIRES].iatoms[i]].orires.ex;
  88             if (ex >= od->nex)
  89             {
  90                 srenew(nr_ex,ex+1);
  91                 for(j=od->nex; j<ex+1; j++)
  92                 {
  93                     nr_ex[j] = 0;
  94             }
  95                 od->nex = ex+1;
  96             }
  97             nr_ex[ex]++;
  98         }
  99     }
 100     snew(od->S,od->nex);
 101     /* When not doing time averaging, the instaneous and time averaged data
 102      * are indentical and the pointers can point to the same memory.
 103      */
 104     snew(od->Dinsl,od->nr);
 105     if (ms)
 106     {
 107         snew(od->Dins,od->nr);
 108     }
 109     else
 110     {
 111         od->Dins = od->Dinsl;
 112     }
 113
 114     if (ir->orires_tau == 0)
 115     {
 116         od->Dtav = od->Dins;
 117         od->edt  = 0.0;
 118         od->edt_1= 1.0;
 119     }
 120     else
 121     {
 122         snew(od->Dtav,od->nr);
 123         od->edt  = exp(-ir->delta_t/ir->orires_tau);
 124         od->edt_1= 1.0 - od->edt;
 125
 126         /* Extend the state with the orires history */
 127         state->flags |= (1<<estORIRE_INITF);
 128         state->hist.orire_initf = 1;
 129         state->flags |= (1<<estORIRE_DTAV);
 130         state->hist.norire_Dtav = od->nr*5;
 131         snew(state->hist.orire_Dtav,state->hist.norire_Dtav);
 132     }
 133
 134     snew(od->oinsl,od->nr);
 135     if (ms)
 136     {
 137         snew(od->oins,od->nr);
 138     }
 139     else
 140     {
 141         od->oins = od->oinsl;
 142     }
 143     if (ir->orires_tau == 0) {
 144         od->otav = od->oins;
 145     }
 146     else
 147     {
 148         snew(od->otav,od->nr);
 149     }
 150     snew(od->tmp,od->nex);
 151     snew(od->TMP,od->nex);
 152     for(ex=0; ex<od->nex; ex++)
 153     {
 154         snew(od->TMP[ex],5);
 155         for(i=0; i<5; i++)
 156         {
 157             snew(od->TMP[ex][i],5);
 158         }
 159     }
 160
 161     od->nref = 0;
 162     for(i=0; i<mtop->natoms; i++)
 163     {
 164         if (ggrpnr(&mtop->groups,egcORFIT,i) == 0)
 165         {
 166             od->nref++;
 167         }
 168     }
 169     snew(od->mref,od->nref);
 170     snew(od->xref,od->nref);
 171     snew(od->xtmp,od->nref);
 172
 173     snew(od->eig,od->nex*12);
 174
 175     /* Determine the reference structure on the master node.
 176      * Copy it to the other nodes after checking multi compatibility,
 177      * so we are sure the subsystems match before copying.
 178      */
 179     clear_rvec(com);
 180     mtot = 0.0;
 181     j = 0;
 182     aloop = gmx_mtop_atomloop_all_init(mtop);
 183     while(gmx_mtop_atomloop_all_next(aloop,&i,&atom))
 184     {
 185         if (mtop->groups.grpnr[egcORFIT] == NULL ||
 186             mtop->groups.grpnr[egcORFIT][i] == 0)
 187         {
 188             /* Not correct for free-energy with changing masses */
 189             od->mref[j] = atom->m;
 190             if (ms==NULL || MASTERSIM(ms))
 191             {
 192                 copy_rvec(xref[i],od->xref[j]);
 193                 for(d=0; d<DIM; d++)
 194                 {
 195                     com[d] += od->mref[j]*xref[i][d];
 196                 }
 197             }
 198             mtot += od->mref[j];
 199             j++;
 200         }
 201     }
 202     svmul(1.0/mtot,com,com);
 203     if (ms==NULL || MASTERSIM(ms))
 204     {
 205         for(j=0; j<od->nref; j++)
 206         {
 207             rvec_dec(od->xref[j],com);
 208         }
 209     }
 210
 211     fprintf(fplog,"Found %d orientation experiments\n",od->nex);
 212     for(i=0; i<od->nex; i++)
 213     {
 214         fprintf(fplog,"  experiment %d has %d restraints\n",i+1,nr_ex[i]);
 215     }
 216
 217     sfree(nr_ex);
 218
 219     fprintf(fplog,"  the fit group consists of %d atoms and has total mass %g\n",
 220             od->nref,mtot);
 221
 222     if (ms)
 223     {
 224         fprintf(fplog,"  the orientation restraints are ensemble averaged over %d systems\n",ms->nsim);
 225
 226         check_multi_int(fplog,ms,od->nr,
 227                         "the number of orientation restraints");
 228         check_multi_int(fplog,ms,od->nref,
 229                         "the number of fit atoms for orientation restraining");
 230         check_multi_int(fplog,ms,ir->nsteps,"nsteps");
 231         /* Copy the reference coordinates from the master to the other nodes */
 232         gmx_sum_sim(DIM*od->nref,od->xref[0],ms);
 233     }
 234
 235     please_cite(fplog,"Hess2003");
 236 }
 237
 238 void diagonalize_orires_tensors(t_oriresdata *od)
 239 {
 240     int           ex,i,j,nrot,ord[DIM],t;
 241     matrix        S,TMP;
 242
 243     if (od->M == NULL)
 244     {
 245         snew(od->M,DIM);
 246         for(i=0; i<DIM; i++)
 247         {
 248             snew(od->M[i],DIM);
 249         }
 250         snew(od->eig_diag,DIM);
 251         snew(od->v,DIM);
 252         for(i=0; i<DIM; i++)
 253         {
 254             snew(od->v[i],DIM);
 255         }
 256     }
 257
 258     for(ex=0; ex<od->nex; ex++)
 259     {
 260         /* Rotate the S tensor back to the reference frame */
 261         mmul(od->R,od->S[ex],TMP);
 262         mtmul(TMP,od->R,S);
 263         for(i=0; i<DIM; i++)
 264         {
 265             for(j=0; j<DIM; j++)
 266             {
 267                 od->M[i][j] = S[i][j];
 268             }
 269         }
 270
 271         jacobi(od->M,DIM,od->eig_diag,od->v,&nrot);
 272
 273         for(i=0; i<DIM; i++)
 274         {
 275             ord[i] = i;
 276         }
 277         for(i=0; i<DIM; i++)
 278         {
 279             for(j=i+1; j<DIM; j++)
 280             {
 281                 if (sqr(od->eig_diag[ord[j]]) > sqr(od->eig_diag[ord[i]]))
 282                 {
 283                     t = ord[i];
 284                     ord[i] = ord[j];
 285                     ord[j] = t;
 286                 }
 287             }
 288         }
 289
 290         for(i=0; i<DIM; i++)
 291         {
 292             od->eig[ex*12 + i] = od->eig_diag[ord[i]];
 293         }
 294         for(i=0; i<DIM; i++)
 295         {
 296             for(j=0; j<DIM; j++)
 297             {
 298                 od->eig[ex*12 + 3 + 3*i + j] = od->v[j][ord[i]];
 299             }
 300         }
 301     }
 302 }
 303
 304 void print_orires_log(FILE *log,t_oriresdata *od)
 305 {
 306     int  ex,i;
 307     real *eig;
 308
 309     diagonalize_orires_tensors(od);
 310
 311     for(ex=0; ex<od->nex; ex++)
 312     {
 313         eig = od->eig + ex*12;
 314         fprintf(log,"  Orientation experiment %d:\n",ex+1);
 315         fprintf(log,"    order parameter: %g\n",eig[0]);
 316         for(i=0; i<DIM; i++)
 317         {
 318             fprintf(log,"    eig: %6.3f   %6.3f %6.3f %6.3f\n",
 319                     (eig[0] != 0) ? eig[i]/eig[0] : eig[i],
 320                     eig[DIM+i*DIM+XX],
 321                     eig[DIM+i*DIM+YY],
 322                     eig[DIM+i*DIM+ZZ]);
 323         }
 324         fprintf(log,"\n");
 325     }
 326 }
 327
 328 real calc_orires_dev(const gmx_multisim_t *ms,
 329                      int nfa,const t_iatom forceatoms[],const t_iparams ip[],
 330                      const t_mdatoms *md,const rvec x[],const t_pbc *pbc,
 331                      t_fcdata *fcd,history_t *hist)
 332 {
 333     int          fa,d,i,j,type,ex,nref;
 334     real         edt,edt_1,invn,pfac,r2,invr,corrfac,weight,wsv2,sw,dev;
 335     tensor       *S,R,TMP;
 336     rvec5        *Dinsl,*Dins,*Dtav,*rhs;
 337     real         *mref,***T;
 338     double       mtot;
 339     rvec         *xref,*xtmp,com,r_unrot,r;
 340     t_oriresdata *od;
 341     gmx_bool         bTAV;
 342     const real   two_thr=2.0/3.0;
 343
 344     od = &(fcd->orires);
 345
 346     if (od->nr == 0)
 347     {
 348         /* This means that this is not the master node */
 349         gmx_fatal(FARGS,"Orientation restraints are only supported on the master node, use less processors");
 350     }
 351
 352     bTAV = (od->edt != 0);
 353     edt  = od->edt;
 354     edt_1= od->edt_1;
 355     S    = od->S;
 356     Dinsl= od->Dinsl;
 357     Dins = od->Dins;
 358     Dtav = od->Dtav;
 359     T    = od->TMP;
 360     rhs  = od->tmp;
 361     nref = od->nref;
 362     mref = od->mref;
 363     xref = od->xref;
 364     xtmp = od->xtmp;
 365
 366     if (bTAV)
 367     {
 368         od->exp_min_t_tau = hist->orire_initf*edt;
 369
 370         /* Correction factor to correct for the lack of history
 371          * at short times.
 372          */
 373         corrfac = 1.0/(1.0 - od->exp_min_t_tau);
 374     }
 375     else
 376     {
 377         corrfac = 1.0;
 378     }
 379
 380     if (ms)
 381     {
 382         invn = 1.0/ms->nsim;
 383     }
 384     else
 385     {
 386         invn = 1.0;
 387     }
 388
 389     clear_rvec(com);
 390     mtot = 0;
 391     j=0;
 392     for(i=0; i<md->nr; i++)
 393     {
 394         if (md->cORF[i] == 0)
 395         {
 396             copy_rvec(x[i],xtmp[j]);
 397             mref[j] = md->massT[i];
 398             for(d=0; d<DIM; d++)
 399             {
 400                 com[d] += mref[j]*xref[j][d];
 401             }
 402             mtot += mref[j];
 403             j++;
 404         }
 405     }
 406     svmul(1.0/mtot,com,com);
 407     for(j=0; j<nref; j++)
 408     {
 409         rvec_dec(xtmp[j],com);
 410     }
 411     /* Calculate the rotation matrix to rotate x to the reference orientation */
 412     calc_fit_R(DIM,nref,mref,xref,xtmp,R);
 413     copy_mat(R,od->R);
 414
 415     d = 0;
 416     for(fa=0; fa<nfa; fa+=3)
 417     {
 418         type = forceatoms[fa];
 419         if (pbc)
 420         {
 421             pbc_dx_aiuc(pbc,x[forceatoms[fa+1]],x[forceatoms[fa+2]],r_unrot);
 422         }
 423         else
 424         {
 425             rvec_sub(x[forceatoms[fa+1]],x[forceatoms[fa+2]],r_unrot);
 426         }
 427         mvmul(R,r_unrot,r);
 428         r2   = norm2(r);
 429         invr = gmx_invsqrt(r2);
 430         /* Calculate the prefactor for the D tensor, this includes the factor 3! */
 431         pfac = ip[type].orires.c*invr*invr*3;
 432         for(i=0; i<ip[type].orires.power; i++)
 433         {
 434             pfac *= invr;
 435         }
 436         Dinsl[d][0] = pfac*(2*r[0]*r[0] + r[1]*r[1] - r2);
 437         Dinsl[d][1] = pfac*(2*r[0]*r[1]);
 438         Dinsl[d][2] = pfac*(2*r[0]*r[2]);
 439         Dinsl[d][3] = pfac*(2*r[1]*r[1] + r[0]*r[0] - r2);
 440         Dinsl[d][4] = pfac*(2*r[1]*r[2]);
 441
 442         if (ms)
 443         {
 444             for(i=0; i<5; i++)
 445             {
 446                 Dins[d][i] = Dinsl[d][i]*invn;
 447             }
 448         }
 449
 450         d++;
 451     }
 452
 453     if (ms)
 454     {
 455         gmx_sum_sim(5*od->nr,Dins[0],ms);
 456     }
 457
 458     /* Calculate the order tensor S for each experiment via optimization */
 459     for(ex=0; ex<od->nex; ex++)
 460     {
 461         for(i=0; i<5; i++)
 462         {
 463             rhs[ex][i] = 0;
 464             for(j=0; j<=i; j++)
 465             {
 466                 T[ex][i][j] = 0;
 467             }
 468         }
 469     }
 470     d = 0;
 471     for(fa=0; fa<nfa; fa+=3)
 472     {
 473         if (bTAV)
 474         {
 475             /* Here we update Dtav in t_fcdata using the data in history_t.
 476              * Thus the results stay correct when this routine
 477              * is called multiple times.
 478              */
 479             for(i=0; i<5; i++)
 480             {
 481                 Dtav[d][i] = edt*hist->orire_Dtav[d*5+i] + edt_1*Dins[d][i];
 482             }
 483         }
 484
 485         type   = forceatoms[fa];
 486         ex     = ip[type].orires.ex;
 487         weight = ip[type].orires.kfac;
 488         /* Calculate the vector rhs and half the matrix T for the 5 equations */
 489         for(i=0; i<5; i++) {
 490             rhs[ex][i] += Dtav[d][i]*ip[type].orires.obs*weight;
 491             for(j=0; j<=i; j++)
 492             {
 493                 T[ex][i][j] += Dtav[d][i]*Dtav[d][j]*weight;
 494             }
 495         }
 496         d++;
 497     }
 498     /* Now we have all the data we can calculate S */
 499     for(ex=0; ex<od->nex; ex++)
 500     {
 501         /* Correct corrfac and copy one half of T to the other half */
 502         for(i=0; i<5; i++)
 503         {
 504             rhs[ex][i]  *= corrfac;
 505             T[ex][i][i] *= sqr(corrfac);
 506             for(j=0; j<i; j++)
 507             {
 508                 T[ex][i][j] *= sqr(corrfac);
 509                 T[ex][j][i]  = T[ex][i][j];
 510             }
 511         }
 512         m_inv_gen(T[ex],5,T[ex]);
 513         /* Calculate the orientation tensor S for this experiment */
 514         S[ex][0][0] = 0;
 515         S[ex][0][1] = 0;
 516         S[ex][0][2] = 0;
 517         S[ex][1][1] = 0;
 518         S[ex][1][2] = 0;
 519         for(i=0; i<5; i++)
 520         {
 521             S[ex][0][0] += 1.5*T[ex][0][i]*rhs[ex][i];
 522             S[ex][0][1] += 1.5*T[ex][1][i]*rhs[ex][i];
 523             S[ex][0][2] += 1.5*T[ex][2][i]*rhs[ex][i];
 524             S[ex][1][1] += 1.5*T[ex][3][i]*rhs[ex][i];
 525             S[ex][1][2] += 1.5*T[ex][4][i]*rhs[ex][i];
 526         }
 527         S[ex][1][0] = S[ex][0][1];
 528         S[ex][2][0] = S[ex][0][2];
 529         S[ex][2][1] = S[ex][1][2];
 530         S[ex][2][2] = -S[ex][0][0] - S[ex][1][1];
 531     }
 532
 533     wsv2 = 0;
 534     sw   = 0;
 535
 536     d = 0;
 537     for(fa=0; fa<nfa; fa+=3)
 538     {
 539         type = forceatoms[fa];
 540         ex = ip[type].orires.ex;
 541
 542         od->otav[d] = two_thr*
 543             corrfac*(S[ex][0][0]*Dtav[d][0] + S[ex][0][1]*Dtav[d][1] +
 544                      S[ex][0][2]*Dtav[d][2] + S[ex][1][1]*Dtav[d][3] +
 545                      S[ex][1][2]*Dtav[d][4]);
 546         if (bTAV)
 547         {
 548             od->oins[d] = two_thr*(S[ex][0][0]*Dins[d][0] + S[ex][0][1]*Dins[d][1] +
 549                                    S[ex][0][2]*Dins[d][2] + S[ex][1][1]*Dins[d][3] +
 550                                    S[ex][1][2]*Dins[d][4]);
 551         }
 552         if (ms)
 553         {
 554             /* When ensemble averaging is used recalculate the local orientation
 555              * for output to the energy file.
 556              */
 557             od->oinsl[d] = two_thr*
 558                 (S[ex][0][0]*Dinsl[d][0] + S[ex][0][1]*Dinsl[d][1] +
 559                  S[ex][0][2]*Dinsl[d][2] + S[ex][1][1]*Dinsl[d][3] +
 560                  S[ex][1][2]*Dinsl[d][4]);
 561         }
 562
 563         dev = od->otav[d] - ip[type].orires.obs;
 564
 565         wsv2 += ip[type].orires.kfac*sqr(dev);
 566         sw   += ip[type].orires.kfac;
 567
 568         d++;
 569     }
 570     od->rmsdev = sqrt(wsv2/sw);
 571
 572     /* Rotate the S matrices back, so we get the correct grad(tr(S D)) */
 573     for(ex=0; ex<od->nex; ex++)
 574     {
 575         tmmul(R,S[ex],TMP);
 576         mmul(TMP,R,S[ex]);
 577     }
 578
 579     return od->rmsdev;
 580
 581     /* Approx. 120*nfa/3 flops */
 582 }
 583
 584 real orires(int nfa,const t_iatom forceatoms[],const t_iparams ip[],
 585             const rvec x[],rvec f[],rvec fshift[],
 586             const t_pbc *pbc,const t_graph *g,
 587             real lambda,real *dvdlambda,
 588             const t_mdatoms *md,t_fcdata *fcd,
 589             int *global_atom_index)
 590 {
 591     atom_id      ai,aj;
 592     int          fa,d,i,type,ex,power,ki=CENTRAL;
 593     ivec         dt;
 594     real         r2,invr,invr2,fc,smooth_fc,dev,devins,pfac;
 595     rvec         r,Sr,fij;
 596     real         vtot;
 597     const t_oriresdata *od;
 598     gmx_bool         bTAV;
 599
 600     vtot = 0;
 601     od = &(fcd->orires);
 602
 603     if (od->fc != 0)
 604     {
 605         bTAV = (od->edt != 0);
 606
 607         smooth_fc = od->fc;
 608         if (bTAV)
 609         {
 610             /* Smoothly switch on the restraining when time averaging is used */
 611             smooth_fc *= (1.0 - od->exp_min_t_tau);
 612         }
 613
 614         d = 0;
 615         for(fa=0; fa<nfa; fa+=3)
 616         {
 617             type  = forceatoms[fa];
 618             ai    = forceatoms[fa+1];
 619             aj    = forceatoms[fa+2];
 620             if (pbc)
 621             {
 622                 ki = pbc_dx_aiuc(pbc,x[ai],x[aj],r);
 623             }
 624             else
 625             {
 626                 rvec_sub(x[ai],x[aj],r);
 627             }
 628             r2    = norm2(r);
 629             invr  = gmx_invsqrt(r2);
 630             invr2 = invr*invr;
 631             ex    = ip[type].orires.ex;
 632             power = ip[type].orires.power;
 633             fc    = smooth_fc*ip[type].orires.kfac;
 634             dev   = od->otav[d] - ip[type].orires.obs;
 635
 636             /* NOTE:
 637              * there is no real potential when time averaging is applied
 638              */
 639             vtot += 0.5*fc*sqr(dev);
 640
 641             if (bTAV)
 642             {
 643                 /* Calculate the force as the sqrt of tav times instantaneous */
 644                 devins = od->oins[d] - ip[type].orires.obs;
 645                 if (dev*devins <= 0)
 646                 {
 647                     dev = 0;
 648                 }
 649                 else
 650                 {
 651                     dev = sqrt(dev*devins);
 652                     if (devins < 0)
 653                     {
 654                         dev = -dev;
 655                     }
 656                 }
 657             }
 658
 659             pfac  = fc*ip[type].orires.c*invr2;
 660             for(i=0; i<power; i++)
 661             {
 662                 pfac *= invr;
 663             }
 664             mvmul(od->S[ex],r,Sr);
 665             for(i=0; i<DIM; i++)
 666             {
 667                 fij[i] =
 668                     -pfac*dev*(4*Sr[i] - 2*(2+power)*invr2*iprod(Sr,r)*r[i]);
 669             }
 670
 671             if (g)
 672             {
 673                 ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
 674                 ki=IVEC2IS(dt);
 675             }
 676
 677             for(i=0; i<DIM; i++)
 678             {
 679                 f[ai][i]           += fij[i];
 680                 f[aj][i]           -= fij[i];
 681                 fshift[ki][i]      += fij[i];
 682                 fshift[CENTRAL][i] -= fij[i];
 683             }
 684             d++;
 685         }
 686     }
 687
 688     return vtot;
 689
 690     /* Approx. 80*nfa/3 flops */
 691 }
 692
 693 void update_orires_history(t_fcdata *fcd,history_t *hist)
 694 {
 695     t_oriresdata *od;
 696     int pair,i;
 697
 698     od = &(fcd->orires);
 699     if (od->edt != 0)
 700     {
 701         /* Copy the new time averages that have been calculated
 702          *  in calc_orires_dev.
 703          */
 704         hist->orire_initf = od->exp_min_t_tau;
 705         for(pair=0; pair<od->nr; pair++)
 706         {
 707             for(i=0; i<5; i++)
 708             {
 709                 hist->orire_Dtav[pair*5+i] = od->Dtav[pair][i];
 710             }
 711         }
 712     }
 713 }