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