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