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