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