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