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