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