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