src/gromacs/gmxlib/disre.c

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