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