src/gromacs/listed-forces/orires.cpp

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