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