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