src/gromacs/listed_forces/disre.cpp

   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,2015,2016,2017 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 /* This file is completely threadsafe - keep it that way! */
  39 #include "gmxpre.h"
  40
  41 #include "disre.h"
  42
  43 #include "config.h"
  44
  45 #include <cmath>
  46 #include <cstdlib>
  47 #include <cstring>
  48
  49 #include <algorithm>
  50
  51 #include "gromacs/gmxlib/network.h"
  52 #include "gromacs/math/functions.h"
  53 #include "gromacs/math/vec.h"
  54 #include "gromacs/mdrunutility/multisim.h"
  55 #include "gromacs/mdtypes/commrec.h"
  56 #include "gromacs/mdtypes/fcdata.h"
  57 #include "gromacs/mdtypes/inputrec.h"
  58 #include "gromacs/mdtypes/mdatom.h"
  59 #include "gromacs/mdtypes/md_enums.h"
  60 #include "gromacs/mdtypes/state.h"
  61 #include "gromacs/pbcutil/ishift.h"
  62 #include "gromacs/pbcutil/pbc.h"
  63 #include "gromacs/topology/mtop_util.h"
  64 #include "gromacs/topology/topology.h"
  65 #include "gromacs/utility/arrayref.h"
  66 #include "gromacs/utility/fatalerror.h"
  67 #include "gromacs/utility/futil.h"
  68 #include "gromacs/utility/gmxassert.h"
  69 #include "gromacs/utility/pleasecite.h"
  70 #include "gromacs/utility/smalloc.h"
  71
  72 void init_disres(FILE*                 fplog,
  73                  const gmx_mtop_t&     mtop,
  74                  t_inputrec*           ir,
  75                  DisResRunMode         disResRunMode,
  76                  DDRole                ddRole,
  77                  NumRanks              numRanks,
  78                  MPI_Comm              communicator,
  79                  const gmx_multisim_t* ms,
  80                  t_disresdata*         dd,
  81                  t_state*              state,
  82                  gmx_bool              bIsREMD)
  83 {
  84     int                  fa, nmol, npair, np;
  85     history_t*           hist;
  86     gmx_mtop_ilistloop_t iloop;
  87     char*                ptr;
  88     int                  type_min, type_max;
  89
  90     if (gmx_mtop_ftype_count(mtop, F_DISRES) == 0)
  91     {
  92         dd->nres = 0;
  93
  94         return;
  95     }
  96
  97     if (fplog)
  98     {
  99         fprintf(fplog, "Initializing the distance restraints\n");
 100     }
 101
 102     dd->dr_weighting = ir->eDisreWeighting;
 103     dd->dr_fc        = ir->dr_fc;
 104     if (EI_DYNAMICS(ir->eI))
 105     {
 106         dd->dr_tau = ir->dr_tau;
 107     }
 108     else
 109     {
 110         dd->dr_tau = 0.0;
 111     }
 112     if (dd->dr_tau == 0.0)
 113     {
 114         dd->dr_bMixed = FALSE;
 115         dd->ETerm     = 0.0;
 116     }
 117     else
 118     {
 119         /* We store the r^-6 time averages in an array that is indexed
 120          * with the local disres iatom index, so this doesn't work with DD.
 121          * Note that DD is not initialized yet here, so we check that we are on multiple ranks,
 122          * but there are probably also issues with e.g. NM MPI parallelization.
 123          */
 124         if ((disResRunMode == DisResRunMode::MDRun) && (numRanks == NumRanks::Multiple))
 125         {
 126             gmx_fatal(FARGS,
 127                       "Time-averaged distance restraints are not supported with MPI "
 128                       "parallelization. You can use OpenMP parallelization on a single node.");
 129         }
 130
 131         dd->dr_bMixed = ir->bDisreMixed;
 132         dd->ETerm     = std::exp(-(ir->delta_t / ir->dr_tau));
 133     }
 134     dd->ETerm1 = 1.0 - dd->ETerm;
 135
 136     dd->nres  = 0;
 137     dd->npair = 0;
 138     type_min  = INT_MAX;
 139     type_max  = 0;
 140     iloop     = gmx_mtop_ilistloop_init(mtop);
 141     while (const InteractionLists* il = gmx_mtop_ilistloop_next(iloop, &nmol))
 142     {
 143         if (nmol > 1 && !(*il)[F_DISRES].empty() && ir->eDisre != DistanceRestraintRefinement::Ensemble)
 144         {
 145             gmx_fatal(FARGS,
 146                       "NMR distance restraints with multiple copies of the same molecule are "
 147                       "currently only supported with ensemble averaging. If you just want to "
 148                       "restrain distances between atom pairs using a flat-bottomed potential, use "
 149                       "a restraint potential (bonds type 10) instead.");
 150         }
 151
 152         np = 0;
 153         for (fa = 0; fa < (*il)[F_DISRES].size(); fa += 3)
 154         {
 155             int type;
 156
 157             type = (*il)[F_DISRES].iatoms[fa];
 158
 159             np++;
 160             npair = mtop.ffparams.iparams[type].disres.npair;
 161             if (np == npair)
 162             {
 163                 dd->nres += (ir->eDisre == DistanceRestraintRefinement::Ensemble ? 1 : nmol);
 164                 dd->npair += nmol * npair;
 165                 np = 0;
 166
 167                 type_min = std::min(type_min, type);
 168                 type_max = std::max(type_max, type);
 169             }
 170         }
 171     }
 172
 173     if ((disResRunMode == DisResRunMode::MDRun) && (numRanks == NumRanks::Multiple) && ir->nstdisreout > 0)
 174     {
 175         /* With DD we currently only have local pair information available */
 176         gmx_fatal(FARGS,
 177                   "With MPI parallelization distance-restraint pair output is not supported. Use "
 178                   "nstdisreout=0 or use OpenMP parallelization on a single node.");
 179     }
 180
 181     /* For communicating and/or reducing (sums of) r^-6 for pairs over threads
 182      * we use multiple arrays in t_disresdata. We need to have unique indices
 183      * for each restraint that work over threads and MPI ranks. To this end
 184      * we use the type index. These should all be in one block and there should
 185      * be dd->nres types, but we check for this here.
 186      * This setup currently does not allow for multiple copies of the same
 187      * molecule without ensemble averaging, this is check for above.
 188      */
 189     GMX_RELEASE_ASSERT(
 190             type_max - type_min + 1 == dd->nres,
 191             "All distance restraint parameter entries in the topology should be consecutive");
 192
 193     dd->type_min = type_min;
 194
 195     snew(dd->rt, dd->npair);
 196
 197     if (dd->dr_tau != 0.0)
 198     {
 199         GMX_RELEASE_ASSERT(state != nullptr,
 200                            "We need a valid state when using time-averaged distance restraints");
 201
 202         hist = &state->hist;
 203         /* Set the "history lack" factor to 1 */
 204         state->flags |= enumValueToBitMask(StateEntry::DisreInitF);
 205         hist->disre_initf = 1.0;
 206         /* Allocate space for the r^-3 time averages */
 207         state->flags |= enumValueToBitMask(StateEntry::DisreRm3Tav);
 208         hist->ndisrepairs = dd->npair;
 209         snew(hist->disre_rm3tav, hist->ndisrepairs);
 210     }
 211     /* Allocate space for a copy of rm3tav,
 212      * so we can call do_force without modifying the state.
 213      */
 214     snew(dd->rm3tav, dd->npair);
 215
 216     /* Allocate Rt_6 and Rtav_6 consecutively in memory so they can be
 217      * averaged over the processors in one call (in calc_disre_R_6)
 218      */
 219     snew(dd->Rt_6, 2 * dd->nres);
 220     dd->Rtav_6 = &(dd->Rt_6[dd->nres]);
 221
 222     ptr = getenv("GMX_DISRE_ENSEMBLE_SIZE");
 223     if ((disResRunMode == DisResRunMode::MDRun) && ms != nullptr && ptr != nullptr && !bIsREMD)
 224     {
 225 #if GMX_MPI
 226         dd->nsystems = 0;
 227         sscanf(ptr, "%d", &dd->nsystems);
 228         if (fplog)
 229         {
 230             fprintf(fplog, "Found GMX_DISRE_ENSEMBLE_SIZE set to %d systems per ensemble\n", dd->nsystems);
 231         }
 232         /* This check is only valid on MASTER(cr), so probably
 233          * ensemble-averaged distance restraints are broken on more
 234          * than one processor per simulation system. */
 235         if (ddRole == DDRole::Master)
 236         {
 237             check_multi_int(fplog, ms, dd->nsystems, "the number of systems per ensemble", FALSE);
 238         }
 239         gmx_bcast(sizeof(int), &dd->nsystems, communicator);
 240
 241         /* We use to allow any value of nsystems which was a divisor
 242          * of ms->numSimulations_. But this required an extra communicator which
 243          * pulled in mpi.h in nearly all C files.
 244          */
 245         if (!(ms->numSimulations_ == 1 || ms->numSimulations_ == dd->nsystems))
 246         {
 247             gmx_fatal(FARGS,
 248                       "GMX_DISRE_ENSEMBLE_SIZE (%d) is not equal to 1 or the number of systems "
 249                       "(option -multidir) %d",
 250                       dd->nsystems,
 251                       ms->numSimulations_);
 252         }
 253         if (fplog)
 254         {
 255             fprintf(fplog, "Our ensemble consists of systems:");
 256             for (int i = 0; i < dd->nsystems; i++)
 257             {
 258                 fprintf(fplog, " %d", (ms->simulationIndex_ / dd->nsystems) * dd->nsystems + i);
 259             }
 260             fprintf(fplog, "\n");
 261         }
 262 #else
 263         GMX_UNUSED_VALUE(communicator);
 264 #endif
 265     }
 266     else
 267     {
 268         dd->nsystems = 1;
 269     }
 270
 271     if (dd->nsystems == 1)
 272     {
 273         dd->Rtl_6 = dd->Rt_6;
 274     }
 275     else
 276     {
 277         snew(dd->Rtl_6, dd->nres);
 278     }
 279
 280     if (dd->npair > 0)
 281     {
 282         if (fplog)
 283         {
 284             fprintf(fplog, "There are %d distance restraints involving %d atom pairs\n", dd->nres, dd->npair);
 285         }
 286         /* Have to avoid g_disre de-referencing cr blindly, mdrun not
 287          * doing consistency checks for ensemble-averaged distance
 288          * restraints when that's not happening, and only doing those
 289          * checks from appropriate processes (since check_multi_int is
 290          * too broken to check whether the communication will
 291          * succeed...) */
 292         if ((disResRunMode == DisResRunMode::MDRun) && ms && dd->nsystems > 1 && (ddRole == DDRole::Master))
 293         {
 294             check_multi_int(fplog, ms, dd->nres, "the number of distance restraints", FALSE);
 295         }
 296         please_cite(fplog, "Tropp80a");
 297         please_cite(fplog, "Torda89a");
 298     }
 299 }
 300
 301 void calc_disres_R_6(const t_commrec*      cr,
 302                      const gmx_multisim_t* ms,
 303                      int                   nfa,
 304                      const t_iatom         forceatoms[],
 305                      const rvec            x[],
 306                      const t_pbc*          pbc,
 307                      t_disresdata*         dd,
 308                      const history_t*      hist)
 309 {
 310     rvec     dx;
 311     real *   rt, *rm3tav, *Rtl_6, *Rt_6, *Rtav_6;
 312     real     ETerm, ETerm1, cf1 = 0, cf2 = 0;
 313     gmx_bool bTav;
 314
 315     bTav   = (dd->dr_tau != 0);
 316     ETerm  = dd->ETerm;
 317     ETerm1 = dd->ETerm1;
 318     rt     = dd->rt;
 319     rm3tav = dd->rm3tav;
 320     Rtl_6  = dd->Rtl_6;
 321     Rt_6   = dd->Rt_6;
 322     Rtav_6 = dd->Rtav_6;
 323
 324     if (bTav)
 325     {
 326         /* scaling factor to smoothly turn on the restraint forces *
 327          * when using time averaging                               */
 328         dd->exp_min_t_tau = hist->disre_initf * ETerm;
 329
 330         cf1 = dd->exp_min_t_tau;
 331         cf2 = 1.0 / (1.0 - dd->exp_min_t_tau);
 332     }
 333
 334     for (int res = 0; res < dd->nres; res++)
 335     {
 336         Rtav_6[res] = 0.0;
 337         Rt_6[res]   = 0.0;
 338     }
 339
 340     /* 'loop' over all atom pairs (pair_nr=fa/3) involved in restraints, *
 341      * the total number of atoms pairs is nfa/3                          */
 342     for (int fa = 0; fa < nfa; fa += 3)
 343     {
 344         int type = forceatoms[fa];
 345         int res  = type - dd->type_min;
 346         int pair = fa / 3;
 347         int ai   = forceatoms[fa + 1];
 348         int aj   = forceatoms[fa + 2];
 349
 350         if (pbc)
 351         {
 352             pbc_dx_aiuc(pbc, x[ai], x[aj], dx);
 353         }
 354         else
 355         {
 356             rvec_sub(x[ai], x[aj], dx);
 357         }
 358         real rt2  = iprod(dx, dx);
 359         real rt_1 = gmx::invsqrt(rt2);
 360         real rt_3 = rt_1 * rt_1 * rt_1;
 361
 362         rt[pair] = rt2 * rt_1;
 363         if (bTav)
 364         {
 365             /* Here we update rm3tav in t_disresdata using the data
 366              * in history_t.
 367              * Thus the results stay correct when this routine
 368              * is called multiple times.
 369              */
 370             rm3tav[pair] = cf2 * ((ETerm - cf1) * hist->disre_rm3tav[pair] + ETerm1 * rt_3);
 371         }
 372         else
 373         {
 374             rm3tav[pair] = rt_3;
 375         }
 376
 377         /* Currently divide_bondeds_over_threads() ensures that pairs within
 378          * the same restraint get assigned to the same thread, so we could
 379          * run this loop thread-parallel.
 380          */
 381         Rt_6[res] += rt_3 * rt_3;
 382         Rtav_6[res] += rm3tav[pair] * rm3tav[pair];
 383     }
 384
 385     /* NOTE: Rt_6 and Rtav_6 are stored consecutively in memory */
 386     if (cr && DOMAINDECOMP(cr))
 387     {
 388         gmx_sum(2 * dd->nres, dd->Rt_6, cr);
 389     }
 390
 391     if (dd->nsystems > 1)
 392     {
 393         real invn = 1.0 / dd->nsystems;
 394
 395         for (int res = 0; res < dd->nres; res++)
 396         {
 397             Rtl_6[res] = Rt_6[res];
 398             Rt_6[res] *= invn;
 399             Rtav_6[res] *= invn;
 400         }
 401
 402         GMX_ASSERT(cr != nullptr && ms != nullptr, "We need multisim with nsystems>1");
 403         gmx_sum_sim(2 * dd->nres, dd->Rt_6, ms);
 404
 405         if (DOMAINDECOMP(cr))
 406         {
 407             gmx_bcast(2 * dd->nres, dd->Rt_6, cr->mpi_comm_mygroup);
 408         }
 409     }
 410
 411     /* Store the base forceatoms pointer, so we can re-calculate the pair
 412      * index in ta_disres() for indexing pair data in t_disresdata when
 413      * using thread parallelization.
 414      */
 415     dd->forceatomsStart = forceatoms;
 416
 417     dd->sumviol = 0;
 418 }
 419
 420 real ta_disres(int              nfa,
 421                const t_iatom*   forceatoms,
 422                const t_iparams* ip,
 423                const rvec*      x,
 424                rvec4*           f,
 425                rvec*            fshift,
 426                const t_pbc*     pbc,
 427                real gmx_unused lambda,
 428                real gmx_unused* dvdlambda,
 429                gmx::ArrayRef<const real> /*charge*/,
 430                t_fcdata gmx_unused* fcd,
 431                t_disresdata*        disresdata,
 432                t_oriresdata gmx_unused* oriresdata,
 433                int gmx_unused* global_atom_index)
 434 {
 435     const real seven_three = 7.0 / 3.0;
 436
 437     rvec     dx;
 438     real     weight_rt_1;
 439     real     smooth_fc, Rt, Rtav, rt2, *Rtl_6, *Rt_6, *Rtav_6;
 440     real     k0, f_scal = 0, fmax_scal, fk_scal, fij;
 441     real     tav_viol, instant_viol, mixed_viol, violtot, vtot;
 442     real     tav_viol_Rtav7, instant_viol_Rtav7;
 443     real     up1, up2, low;
 444     gmx_bool bConservative, bMixed, bViolation;
 445     gmx_bool dr_bMixed;
 446
 447     DistanceRestraintWeighting dr_weighting = disresdata->dr_weighting;
 448     dr_bMixed                               = disresdata->dr_bMixed;
 449     Rtl_6                                   = disresdata->Rtl_6;
 450     Rt_6                                    = disresdata->Rt_6;
 451     Rtav_6                                  = disresdata->Rtav_6;
 452
 453     tav_viol = instant_viol = mixed_viol = tav_viol_Rtav7 = instant_viol_Rtav7 = 0;
 454
 455     smooth_fc = disresdata->dr_fc;
 456     if (disresdata->dr_tau != 0)
 457     {
 458         /* scaling factor to smoothly turn on the restraint forces *
 459          * when using time averaging                               */
 460         smooth_fc *= (1.0 - disresdata->exp_min_t_tau);
 461     }
 462
 463     violtot = 0;
 464     vtot    = 0;
 465
 466     /* 'loop' over all atom pairs (pair_nr=fa/3) involved in restraints, *
 467      * the total number of atoms pairs is nfa/3                          */
 468     int faOffset = static_cast<int>(forceatoms - disresdata->forceatomsStart);
 469     for (int fa = 0; fa < nfa; fa += 3)
 470     {
 471         int type  = forceatoms[fa];
 472         int npair = ip[type].disres.npair;
 473         up1       = ip[type].disres.up1;
 474         up2       = ip[type].disres.up2;
 475         low       = ip[type].disres.low;
 476         k0        = smooth_fc * ip[type].disres.kfac;
 477
 478         int res = type - disresdata->type_min;
 479
 480         /* save some flops when there is only one pair */
 481         if (ip[type].disres.type != 2)
 482         {
 483             bConservative = (dr_weighting == DistanceRestraintWeighting::Conservative) && (npair > 1);
 484             bMixed        = dr_bMixed;
 485             Rt            = gmx::invsixthroot(Rt_6[res]);
 486             Rtav          = gmx::invsixthroot(Rtav_6[res]);
 487         }
 488         else
 489         {
 490             /* When rtype=2 use instantaneous not ensemble averaged distance */
 491             bConservative = (npair > 1);
 492             bMixed        = FALSE;
 493             Rt            = gmx::invsixthroot(Rtl_6[res]);
 494             Rtav          = Rt;
 495         }
 496
 497         if (Rtav > up1)
 498         {
 499             bViolation = TRUE;
 500             tav_viol   = Rtav - up1;
 501         }
 502         else if (Rtav < low)
 503         {
 504             bViolation = TRUE;
 505             tav_viol   = Rtav - low;
 506         }
 507         else
 508         {
 509             bViolation = FALSE;
 510         }
 511
 512         if (bViolation)
 513         {
 514             /* Add 1/npair energy and violation for each of the npair pairs */
 515             real pairFac = 1 / static_cast<real>(npair);
 516
 517             /* NOTE:
 518              * there is no real potential when time averaging is applied
 519              */
 520             vtot += 0.5 * k0 * gmx::square(tav_viol) * pairFac;
 521             if (!bMixed)
 522             {
 523                 f_scal = -k0 * tav_viol;
 524                 violtot += fabs(tav_viol) * pairFac;
 525             }
 526             else
 527             {
 528                 if (Rt > up1)
 529                 {
 530                     if (tav_viol > 0)
 531                     {
 532                         instant_viol = Rt - up1;
 533                     }
 534                     else
 535                     {
 536                         bViolation = FALSE;
 537                     }
 538                 }
 539                 else if (Rt < low)
 540                 {
 541                     if (tav_viol < 0)
 542                     {
 543                         instant_viol = Rt - low;
 544                     }
 545                     else
 546                     {
 547                         bViolation = FALSE;
 548                     }
 549                 }
 550                 else
 551                 {
 552                     bViolation = FALSE;
 553                 }
 554                 if (bViolation)
 555                 {
 556                     mixed_viol = std::sqrt(tav_viol * instant_viol);
 557                     f_scal     = -k0 * mixed_viol;
 558                     violtot += mixed_viol * pairFac;
 559                 }
 560             }
 561         }
 562
 563         if (bViolation)
 564         {
 565             fmax_scal = -k0 * (up2 - up1);
 566             /* Correct the force for the number of restraints */
 567             if (bConservative)
 568             {
 569                 f_scal = std::max(f_scal, fmax_scal);
 570                 if (!bMixed)
 571                 {
 572                     f_scal *= Rtav / Rtav_6[res];
 573                 }
 574                 else
 575                 {
 576                     f_scal /= 2 * mixed_viol;
 577                     tav_viol_Rtav7     = tav_viol * Rtav / Rtav_6[res];
 578                     instant_viol_Rtav7 = instant_viol * Rt / Rt_6[res];
 579                 }
 580             }
 581             else
 582             {
 583                 f_scal /= npair;
 584                 f_scal = std::max(f_scal, fmax_scal);
 585             }
 586
 587             /* Exert the force ... */
 588
 589             int pair = (faOffset + fa) / 3;
 590             int ai   = forceatoms[fa + 1];
 591             int aj   = forceatoms[fa + 2];
 592             int ki   = CENTRAL;
 593             if (pbc)
 594             {
 595                 ki = pbc_dx_aiuc(pbc, x[ai], x[aj], dx);
 596             }
 597             else
 598             {
 599                 rvec_sub(x[ai], x[aj], dx);
 600             }
 601             rt2 = iprod(dx, dx);
 602
 603             weight_rt_1 = gmx::invsqrt(rt2);
 604
 605             if (bConservative)
 606             {
 607                 if (!dr_bMixed)
 608                 {
 609                     weight_rt_1 *= std::pow(disresdata->rm3tav[pair], seven_three);
 610                 }
 611                 else
 612                 {
 613                     weight_rt_1 *=
 614                             tav_viol_Rtav7 * std::pow(disresdata->rm3tav[pair], seven_three)
 615                             + instant_viol_Rtav7
 616                                       / (disresdata->rt[pair] * gmx::power6(disresdata->rt[pair]));
 617                 }
 618             }
 619
 620             fk_scal = f_scal * weight_rt_1;
 621
 622             for (int m = 0; m < DIM; m++)
 623             {
 624                 fij = fk_scal * dx[m];
 625
 626                 f[ai][m] += fij;
 627                 f[aj][m] -= fij;
 628                 if (fshift)
 629                 {
 630                     fshift[ki][m] += fij;
 631                     fshift[CENTRAL][m] -= fij;
 632                 }
 633             }
 634         }
 635     }
 636
 637 #pragma omp atomic
 638     disresdata->sumviol += violtot;
 639
 640     /* Return energy */
 641     return vtot;
 642 }
 643
 644 void update_disres_history(const t_disresdata& dd, history_t* hist)
 645 {
 646     if (dd.dr_tau != 0)
 647     {
 648         /* Copy the new time averages that have been calculated
 649          * in calc_disres_R_6.
 650          */
 651         hist->disre_initf = dd.exp_min_t_tau;
 652         for (int pair = 0; pair < dd.npair; pair++)
 653         {
 654             hist->disre_rm3tav[pair] = dd.rm3tav[pair];
 655         }
 656     }
 657 }