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