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