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