src/gromacs/gmxlib/disre.c

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