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