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