src/gromacs/gmxlib/orires.c

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