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                         FALSE);
 229         check_multi_int(fplog,ms,od->nref,
 230                         "the number of fit atoms for orientation restraining",
 231                         FALSE);
 232         check_multi_int(fplog,ms,ir->nsteps,"nsteps",FALSE);
 233         /* Copy the reference coordinates from the master to the other nodes */
 234         gmx_sum_sim(DIM*od->nref,od->xref[0],ms);
 235     }
 236
 237     please_cite(fplog,"Hess2003");
 238 }
 239
 240 void diagonalize_orires_tensors(t_oriresdata *od)
 241 {
 242     int           ex,i,j,nrot,ord[DIM],t;
 243     matrix        S,TMP;
 244
 245     if (od->M == NULL)
 246     {
 247         snew(od->M,DIM);
 248         for(i=0; i<DIM; i++)
 249         {
 250             snew(od->M[i],DIM);
 251         }
 252         snew(od->eig_diag,DIM);
 253         snew(od->v,DIM);
 254         for(i=0; i<DIM; i++)
 255         {
 256             snew(od->v[i],DIM);
 257         }
 258     }
 259
 260     for(ex=0; ex<od->nex; ex++)
 261     {
 262         /* Rotate the S tensor back to the reference frame */
 263         mmul(od->R,od->S[ex],TMP);
 264         mtmul(TMP,od->R,S);
 265         for(i=0; i<DIM; i++)
 266         {
 267             for(j=0; j<DIM; j++)
 268             {
 269                 od->M[i][j] = S[i][j];
 270             }
 271         }
 272
 273         jacobi(od->M,DIM,od->eig_diag,od->v,&nrot);
 274
 275         for(i=0; i<DIM; i++)
 276         {
 277             ord[i] = i;
 278         }
 279         for(i=0; i<DIM; i++)
 280         {
 281             for(j=i+1; j<DIM; j++)
 282             {
 283                 if (sqr(od->eig_diag[ord[j]]) > sqr(od->eig_diag[ord[i]]))
 284                 {
 285                     t = ord[i];
 286                     ord[i] = ord[j];
 287                     ord[j] = t;
 288                 }
 289             }
 290         }
 291
 292         for(i=0; i<DIM; i++)
 293         {
 294             od->eig[ex*12 + i] = od->eig_diag[ord[i]];
 295         }
 296         for(i=0; i<DIM; i++)
 297         {
 298             for(j=0; j<DIM; j++)
 299             {
 300                 od->eig[ex*12 + 3 + 3*i + j] = od->v[j][ord[i]];
 301             }
 302         }
 303     }
 304 }
 305
 306 void print_orires_log(FILE *log,t_oriresdata *od)
 307 {
 308     int  ex,i;
 309     real *eig;
 310
 311     diagonalize_orires_tensors(od);
 312
 313     for(ex=0; ex<od->nex; ex++)
 314     {
 315         eig = od->eig + ex*12;
 316         fprintf(log,"  Orientation experiment %d:\n",ex+1);
 317         fprintf(log,"    order parameter: %g\n",eig[0]);
 318         for(i=0; i<DIM; i++)
 319         {
 320             fprintf(log,"    eig: %6.3f   %6.3f %6.3f %6.3f\n",
 321                     (eig[0] != 0) ? eig[i]/eig[0] : eig[i],
 322                     eig[DIM+i*DIM+XX],
 323                     eig[DIM+i*DIM+YY],
 324                     eig[DIM+i*DIM+ZZ]);
 325         }
 326         fprintf(log,"\n");
 327     }
 328 }
 329
 330 real calc_orires_dev(const gmx_multisim_t *ms,
 331                      int nfa,const t_iatom forceatoms[],const t_iparams ip[],
 332                      const t_mdatoms *md,const rvec x[],const t_pbc *pbc,
 333                      t_fcdata *fcd,history_t *hist)
 334 {
 335     int          fa,d,i,j,type,ex,nref;
 336     real         edt,edt_1,invn,pfac,r2,invr,corrfac,weight,wsv2,sw,dev;
 337     tensor       *S,R,TMP;
 338     rvec5        *Dinsl,*Dins,*Dtav,*rhs;
 339     real         *mref,***T;
 340     double       mtot;
 341     rvec         *xref,*xtmp,com,r_unrot,r;
 342     t_oriresdata *od;
 343     gmx_bool         bTAV;
 344     const real   two_thr=2.0/3.0;
 345
 346     od = &(fcd->orires);
 347
 348     if (od->nr == 0)
 349     {
 350         /* This means that this is not the master node */
 351         gmx_fatal(FARGS,"Orientation restraints are only supported on the master node, use less processors");
 352     }
 353
 354     bTAV = (od->edt != 0);
 355     edt  = od->edt;
 356     edt_1= od->edt_1;
 357     S    = od->S;
 358     Dinsl= od->Dinsl;
 359     Dins = od->Dins;
 360     Dtav = od->Dtav;
 361     T    = od->TMP;
 362     rhs  = od->tmp;
 363     nref = od->nref;
 364     mref = od->mref;
 365     xref = od->xref;
 366     xtmp = od->xtmp;
 367
 368     if (bTAV)
 369     {
 370         od->exp_min_t_tau = hist->orire_initf*edt;
 371
 372         /* Correction factor to correct for the lack of history
 373          * at short times.
 374          */
 375         corrfac = 1.0/(1.0 - od->exp_min_t_tau);
 376     }
 377     else
 378     {
 379         corrfac = 1.0;
 380     }
 381
 382     if (ms)
 383     {
 384         invn = 1.0/ms->nsim;
 385     }
 386     else
 387     {
 388         invn = 1.0;
 389     }
 390
 391     clear_rvec(com);
 392     mtot = 0;
 393     j=0;
 394     for(i=0; i<md->nr; i++)
 395     {
 396         if (md->cORF[i] == 0)
 397         {
 398             copy_rvec(x[i],xtmp[j]);
 399             mref[j] = md->massT[i];
 400             for(d=0; d<DIM; d++)
 401             {
 402                 com[d] += mref[j]*xref[j][d];
 403             }
 404             mtot += mref[j];
 405             j++;
 406         }
 407     }
 408     svmul(1.0/mtot,com,com);
 409     for(j=0; j<nref; j++)
 410     {
 411         rvec_dec(xtmp[j],com);
 412     }
 413     /* Calculate the rotation matrix to rotate x to the reference orientation */
 414     calc_fit_R(DIM,nref,mref,xref,xtmp,R);
 415     copy_mat(R,od->R);
 416
 417     d = 0;
 418     for(fa=0; fa<nfa; fa+=3)
 419     {
 420         type = forceatoms[fa];
 421         if (pbc)
 422         {
 423             pbc_dx_aiuc(pbc,x[forceatoms[fa+1]],x[forceatoms[fa+2]],r_unrot);
 424         }
 425         else
 426         {
 427             rvec_sub(x[forceatoms[fa+1]],x[forceatoms[fa+2]],r_unrot);
 428         }
 429         mvmul(R,r_unrot,r);
 430         r2   = norm2(r);
 431         invr = gmx_invsqrt(r2);
 432         /* Calculate the prefactor for the D tensor, this includes the factor 3! */
 433         pfac = ip[type].orires.c*invr*invr*3;
 434         for(i=0; i<ip[type].orires.power; i++)
 435         {
 436             pfac *= invr;
 437         }
 438         Dinsl[d][0] = pfac*(2*r[0]*r[0] + r[1]*r[1] - r2);
 439         Dinsl[d][1] = pfac*(2*r[0]*r[1]);
 440         Dinsl[d][2] = pfac*(2*r[0]*r[2]);
 441         Dinsl[d][3] = pfac*(2*r[1]*r[1] + r[0]*r[0] - r2);
 442         Dinsl[d][4] = pfac*(2*r[1]*r[2]);
 443
 444         if (ms)
 445         {
 446             for(i=0; i<5; i++)
 447             {
 448                 Dins[d][i] = Dinsl[d][i]*invn;
 449             }
 450         }
 451
 452         d++;
 453     }
 454
 455     if (ms)
 456     {
 457         gmx_sum_sim(5*od->nr,Dins[0],ms);
 458     }
 459
 460     /* Calculate the order tensor S for each experiment via optimization */
 461     for(ex=0; ex<od->nex; ex++)
 462     {
 463         for(i=0; i<5; i++)
 464         {
 465             rhs[ex][i] = 0;
 466             for(j=0; j<=i; j++)
 467             {
 468                 T[ex][i][j] = 0;
 469             }
 470         }
 471     }
 472     d = 0;
 473     for(fa=0; fa<nfa; fa+=3)
 474     {
 475         if (bTAV)
 476         {
 477             /* Here we update Dtav in t_fcdata using the data in history_t.
 478              * Thus the results stay correct when this routine
 479              * is called multiple times.
 480              */
 481             for(i=0; i<5; i++)
 482             {
 483                 Dtav[d][i] = edt*hist->orire_Dtav[d*5+i] + edt_1*Dins[d][i];
 484             }
 485         }
 486
 487         type   = forceatoms[fa];
 488         ex     = ip[type].orires.ex;
 489         weight = ip[type].orires.kfac;
 490         /* Calculate the vector rhs and half the matrix T for the 5 equations */
 491         for(i=0; i<5; i++) {
 492             rhs[ex][i] += Dtav[d][i]*ip[type].orires.obs*weight;
 493             for(j=0; j<=i; j++)
 494             {
 495                 T[ex][i][j] += Dtav[d][i]*Dtav[d][j]*weight;
 496             }
 497         }
 498         d++;
 499     }
 500     /* Now we have all the data we can calculate S */
 501     for(ex=0; ex<od->nex; ex++)
 502     {
 503         /* Correct corrfac and copy one half of T to the other half */
 504         for(i=0; i<5; i++)
 505         {
 506             rhs[ex][i]  *= corrfac;
 507             T[ex][i][i] *= sqr(corrfac);
 508             for(j=0; j<i; j++)
 509             {
 510                 T[ex][i][j] *= sqr(corrfac);
 511                 T[ex][j][i]  = T[ex][i][j];
 512             }
 513         }
 514         m_inv_gen(T[ex],5,T[ex]);
 515         /* Calculate the orientation tensor S for this experiment */
 516         S[ex][0][0] = 0;
 517         S[ex][0][1] = 0;
 518         S[ex][0][2] = 0;
 519         S[ex][1][1] = 0;
 520         S[ex][1][2] = 0;
 521         for(i=0; i<5; i++)
 522         {
 523             S[ex][0][0] += 1.5*T[ex][0][i]*rhs[ex][i];
 524             S[ex][0][1] += 1.5*T[ex][1][i]*rhs[ex][i];
 525             S[ex][0][2] += 1.5*T[ex][2][i]*rhs[ex][i];
 526             S[ex][1][1] += 1.5*T[ex][3][i]*rhs[ex][i];
 527             S[ex][1][2] += 1.5*T[ex][4][i]*rhs[ex][i];
 528         }
 529         S[ex][1][0] = S[ex][0][1];
 530         S[ex][2][0] = S[ex][0][2];
 531         S[ex][2][1] = S[ex][1][2];
 532         S[ex][2][2] = -S[ex][0][0] - S[ex][1][1];
 533     }
 534
 535     wsv2 = 0;
 536     sw   = 0;
 537
 538     d = 0;
 539     for(fa=0; fa<nfa; fa+=3)
 540     {
 541         type = forceatoms[fa];
 542         ex = ip[type].orires.ex;
 543
 544         od->otav[d] = two_thr*
 545             corrfac*(S[ex][0][0]*Dtav[d][0] + S[ex][0][1]*Dtav[d][1] +
 546                      S[ex][0][2]*Dtav[d][2] + S[ex][1][1]*Dtav[d][3] +
 547                      S[ex][1][2]*Dtav[d][4]);
 548         if (bTAV)
 549         {
 550             od->oins[d] = two_thr*(S[ex][0][0]*Dins[d][0] + S[ex][0][1]*Dins[d][1] +
 551                                    S[ex][0][2]*Dins[d][2] + S[ex][1][1]*Dins[d][3] +
 552                                    S[ex][1][2]*Dins[d][4]);
 553         }
 554         if (ms)
 555         {
 556             /* When ensemble averaging is used recalculate the local orientation
 557              * for output to the energy file.
 558              */
 559             od->oinsl[d] = two_thr*
 560                 (S[ex][0][0]*Dinsl[d][0] + S[ex][0][1]*Dinsl[d][1] +
 561                  S[ex][0][2]*Dinsl[d][2] + S[ex][1][1]*Dinsl[d][3] +
 562                  S[ex][1][2]*Dinsl[d][4]);
 563         }
 564
 565         dev = od->otav[d] - ip[type].orires.obs;
 566
 567         wsv2 += ip[type].orires.kfac*sqr(dev);
 568         sw   += ip[type].orires.kfac;
 569
 570         d++;
 571     }
 572     od->rmsdev = sqrt(wsv2/sw);
 573
 574     /* Rotate the S matrices back, so we get the correct grad(tr(S D)) */
 575     for(ex=0; ex<od->nex; ex++)
 576     {
 577         tmmul(R,S[ex],TMP);
 578         mmul(TMP,R,S[ex]);
 579     }
 580
 581     return od->rmsdev;
 582
 583     /* Approx. 120*nfa/3 flops */
 584 }
 585
 586 real orires(int nfa,const t_iatom forceatoms[],const t_iparams ip[],
 587             const rvec x[],rvec f[],rvec fshift[],
 588             const t_pbc *pbc,const t_graph *g,
 589             real lambda,real *dvdlambda,
 590             const t_mdatoms *md,t_fcdata *fcd,
 591             int *global_atom_index)
 592 {
 593     atom_id      ai,aj;
 594     int          fa,d,i,type,ex,power,ki=CENTRAL;
 595     ivec         dt;
 596     real         r2,invr,invr2,fc,smooth_fc,dev,devins,pfac;
 597     rvec         r,Sr,fij;
 598     real         vtot;
 599     const t_oriresdata *od;
 600     gmx_bool         bTAV;
 601
 602     vtot = 0;
 603     od = &(fcd->orires);
 604
 605     if (od->fc != 0)
 606     {
 607         bTAV = (od->edt != 0);
 608
 609         smooth_fc = od->fc;
 610         if (bTAV)
 611         {
 612             /* Smoothly switch on the restraining when time averaging is used */
 613             smooth_fc *= (1.0 - od->exp_min_t_tau);
 614         }
 615
 616         d = 0;
 617         for(fa=0; fa<nfa; fa+=3)
 618         {
 619             type  = forceatoms[fa];
 620             ai    = forceatoms[fa+1];
 621             aj    = forceatoms[fa+2];
 622             if (pbc)
 623             {
 624                 ki = pbc_dx_aiuc(pbc,x[ai],x[aj],r);
 625             }
 626             else
 627             {
 628                 rvec_sub(x[ai],x[aj],r);
 629             }
 630             r2    = norm2(r);
 631             invr  = gmx_invsqrt(r2);
 632             invr2 = invr*invr;
 633             ex    = ip[type].orires.ex;
 634             power = ip[type].orires.power;
 635             fc    = smooth_fc*ip[type].orires.kfac;
 636             dev   = od->otav[d] - ip[type].orires.obs;
 637
 638             /* NOTE:
 639              * there is no real potential when time averaging is applied
 640              */
 641             vtot += 0.5*fc*sqr(dev);
 642
 643             if (bTAV)
 644             {
 645                 /* Calculate the force as the sqrt of tav times instantaneous */
 646                 devins = od->oins[d] - ip[type].orires.obs;
 647                 if (dev*devins <= 0)
 648                 {
 649                     dev = 0;
 650                 }
 651                 else
 652                 {
 653                     dev = sqrt(dev*devins);
 654                     if (devins < 0)
 655                     {
 656                         dev = -dev;
 657                     }
 658                 }
 659             }
 660
 661             pfac  = fc*ip[type].orires.c*invr2;
 662             for(i=0; i<power; i++)
 663             {
 664                 pfac *= invr;
 665             }
 666             mvmul(od->S[ex],r,Sr);
 667             for(i=0; i<DIM; i++)
 668             {
 669                 fij[i] =
 670                     -pfac*dev*(4*Sr[i] - 2*(2+power)*invr2*iprod(Sr,r)*r[i]);
 671             }
 672
 673             if (g)
 674             {
 675                 ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
 676                 ki=IVEC2IS(dt);
 677             }
 678
 679             for(i=0; i<DIM; i++)
 680             {
 681                 f[ai][i]           += fij[i];
 682                 f[aj][i]           -= fij[i];
 683                 fshift[ki][i]      += fij[i];
 684                 fshift[CENTRAL][i] -= fij[i];
 685             }
 686             d++;
 687         }
 688     }
 689
 690     return vtot;
 691
 692     /* Approx. 80*nfa/3 flops */
 693 }
 694
 695 void update_orires_history(t_fcdata *fcd,history_t *hist)
 696 {
 697     t_oriresdata *od;
 698     int pair,i;
 699
 700     od = &(fcd->orires);
 701     if (od->edt != 0)
 702     {
 703         /* Copy the new time averages that have been calculated
 704          *  in calc_orires_dev.
 705          */
 706         hist->orire_initf = od->exp_min_t_tau;
 707         for(pair=0; pair<od->nr; pair++)
 708         {
 709             for(i=0; i<5; i++)
 710             {
 711                 hist->orire_Dtav[pair*5+i] = od->Dtav[pair][i];
 712             }
 713         }
 714     }
 715 }