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