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