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