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 "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         np = 0;
 117         for(fa=0; fa<il[F_DISRES].nr; fa+=3)
 118         {
 119             np++;
 120             npair = mtop->ffparams.iparams[il[F_DISRES].iatoms[fa]].disres.npair;
 121             if (np == npair)
 122             {
 123                 dd->nres  += (ir->eDisre==edrEnsemble ? 1 : nmol)*npair;
 124                 dd->npair += nmol*npair;
 125                 np = 0;
 126             }
 127         }
 128     }
 129
 130     if (cr && PAR(cr) && !bPartDecomp)
 131     {
 132         /* Temporary check, will be removed when disre is implemented with DD */
 133         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)";
 134
 135         if (MASTER(cr))
 136             fprintf(stderr,"\n%s\n\n",notestr);
 137         if (fplog)
 138             fprintf(fplog,"%s\n",notestr);
 139
 140         if (dd->dr_tau != 0 || ir->eDisre == edrEnsemble || cr->ms != NULL ||
 141             dd->nres != dd->npair)
 142         {
 143             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)");
 144         }
 145         if (ir->nstdisreout != 0)
 146         {
 147             if (fplog)
 148             {
 149                 fprintf(fplog,"\nWARNING: Can not write distance restraint data to energy file with domain decomposition\n\n");
 150             }
 151             if (MASTER(cr))
 152             {
 153                 fprintf(stderr,"\nWARNING: Can not write distance restraint data to energy file with domain decomposition\n");
 154             }
 155             ir->nstdisreout = 0;
 156         }
 157     }
 158
 159     snew(dd->rt,dd->npair);
 160
 161     if (dd->dr_tau != 0.0)
 162     {
 163         hist = &state->hist;
 164         /* Set the "history lack" factor to 1 */
 165         state->flags |= (1<<estDISRE_INITF);
 166         hist->disre_initf = 1.0;
 167         /* Allocate space for the r^-3 time averages */
 168         state->flags |= (1<<estDISRE_RM3TAV);
 169         hist->ndisrepairs = dd->npair;
 170         snew(hist->disre_rm3tav,hist->ndisrepairs);
 171     }
 172     /* Allocate space for a copy of rm3tav,
 173      * so we can call do_force without modifying the state.
 174      */
 175     snew(dd->rm3tav,dd->npair);
 176
 177     /* Allocate Rt_6 and Rtav_6 consecutively in memory so they can be
 178      * averaged over the processors in one call (in calc_disre_R_6)
 179      */
 180     snew(dd->Rt_6,2*dd->nres);
 181     dd->Rtav_6 = &(dd->Rt_6[dd->nres]);
 182
 183     ptr = getenv("GMX_DISRE_ENSEMBLE_SIZE");
 184     if (cr && cr->ms != NULL && ptr != NULL && !bIsREMD)
 185     {
 186 #ifdef GMX_MPI
 187         dd->nsystems = 0;
 188         sscanf(ptr,"%d",&dd->nsystems);
 189         if (fplog)
 190         {
 191             fprintf(fplog,"Found GMX_DISRE_ENSEMBLE_SIZE set to %d systems per ensemble\n",dd->nsystems);
 192         }
 193         /* This check is only valid on MASTER(cr), so probably
 194          * ensemble-averaged distance restraints are broken on more
 195          * than one processor per simulation system. */
 196         if (MASTER(cr))
 197         {
 198             check_multi_int(fplog,cr->ms,dd->nsystems,
 199                             "the number of systems per ensemble",
 200                             FALSE);
 201         }
 202         gmx_bcast_sim(sizeof(int), &dd->nsystems, cr);
 203
 204         /* We use to allow any value of nsystems which was a divisor
 205          * of ms->nsim. But this required an extra communicator which
 206          * was stored in t_fcdata. This pulled in mpi.h in nearly all C files.
 207          */
 208         if (!(cr->ms->nsim == 1 || cr->ms->nsim == dd->nsystems))
 209         {
 210             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);
 211         }
 212         if (fplog)
 213         {
 214             fprintf(fplog,"Our ensemble consists of systems:");
 215             for(i=0; i<dd->nsystems; i++)
 216             {
 217                 fprintf(fplog," %d",
 218                         (cr->ms->sim/dd->nsystems)*dd->nsystems+i);
 219             }
 220             fprintf(fplog,"\n");
 221             }
 222         snew(dd->Rtl_6,dd->nres);
 223 #endif
 224     }
 225     else
 226     {
 227         dd->nsystems = 1;
 228         dd->Rtl_6 = dd->Rt_6;
 229     }
 230
 231     if (dd->npair > 0)
 232     {
 233         if (fplog) {
 234             fprintf(fplog,"There are %d distance restraints involving %d atom pairs\n",dd->nres,dd->npair);
 235         }
 236         /* Have to avoid g_disre de-referencing cr blindly, mdrun not
 237          * doing consistency checks for ensemble-averaged distance
 238          * restraints when that's not happening, and only doing those
 239          * checks from appropriate processes (since check_multi_int is
 240          * too broken to check whether the communication will
 241          * succeed...) */
 242         if (cr && cr->ms && dd->nsystems > 1 && MASTER(cr))
 243         {
 244             check_multi_int(fplog,cr->ms,fcd->disres.nres,
 245                             "the number of distance restraints",
 246                             FALSE);
 247         }
 248         please_cite(fplog,"Tropp80a");
 249         please_cite(fplog,"Torda89a");
 250     }
 251 }
 252
 253 void calc_disres_R_6(const gmx_multisim_t *ms,
 254                      int nfa,const t_iatom forceatoms[],const t_iparams ip[],
 255                      const rvec x[],const t_pbc *pbc,
 256                      t_fcdata *fcd,history_t *hist)
 257 {
 258     atom_id     ai,aj;
 259     int         fa,res,i,pair,ki,kj,m;
 260     int         type,npair,np;
 261     rvec        dx;
 262     real        *rt,*rm3tav,*Rtl_6,*Rt_6,*Rtav_6;
 263     real        rt_1,rt_3,rt2;
 264     ivec        it,jt,dt;
 265     t_disresdata *dd;
 266     real        ETerm,ETerm1,cf1=0,cf2=0,invn=0;
 267     gmx_bool        bTav;
 268
 269     dd = &(fcd->disres);
 270     bTav         = (dd->dr_tau != 0);
 271     ETerm        = dd->ETerm;
 272     ETerm1       = dd->ETerm1;
 273     rt           = dd->rt;
 274     rm3tav       = dd->rm3tav;
 275     Rtl_6        = dd->Rtl_6;
 276     Rt_6         = dd->Rt_6;
 277     Rtav_6       = dd->Rtav_6;
 278
 279     if (bTav)
 280     {
 281         /* scaling factor to smoothly turn on the restraint forces *
 282          * when using time averaging                               */
 283         dd->exp_min_t_tau = hist->disre_initf*ETerm;
 284
 285         cf1 = dd->exp_min_t_tau;
 286         cf2 = 1.0/(1.0 - dd->exp_min_t_tau);
 287     }
 288
 289     if (dd->nsystems > 1)
 290     {
 291         invn = 1.0/dd->nsystems;
 292     }
 293
 294     /* 'loop' over all atom pairs (pair_nr=fa/3) involved in restraints, *
 295      * the total number of atoms pairs is nfa/3                          */
 296     res = 0;
 297     fa  = 0;
 298     while (fa < nfa)
 299     {
 300         type  = forceatoms[fa];
 301         npair = ip[type].disres.npair;
 302
 303         Rtav_6[res] = 0.0;
 304         Rt_6[res]   = 0.0;
 305
 306         /* Loop over the atom pairs of 'this' restraint */
 307         np = 0;
 308         while (fa < nfa && np < npair)
 309         {
 310             pair = fa/3;
 311             ai   = forceatoms[fa+1];
 312             aj   = forceatoms[fa+2];
 313
 314             if (pbc)
 315             {
 316                 pbc_dx_aiuc(pbc,x[ai],x[aj],dx);
 317             }
 318             else
 319             {
 320                 rvec_sub(x[ai],x[aj],dx);
 321             }
 322             rt2  = iprod(dx,dx);
 323             rt_1 = gmx_invsqrt(rt2);
 324             rt_3 = rt_1*rt_1*rt_1;
 325
 326             rt[pair]         = sqrt(rt2);
 327             if (bTav)
 328             {
 329                 /* Here we update rm3tav in t_fcdata using the data
 330                  * in history_t.
 331                  * Thus the results stay correct when this routine
 332                  * is called multiple times.
 333                  */
 334                 rm3tav[pair] = cf2*((ETerm - cf1)*hist->disre_rm3tav[pair] +
 335                                     ETerm1*rt_3);
 336             }
 337             else
 338             {
 339                 rm3tav[pair] = rt_3;
 340             }
 341
 342             Rt_6[res]       += rt_3*rt_3;
 343             Rtav_6[res]     += rm3tav[pair]*rm3tav[pair];
 344
 345             fa += 3;
 346             np++;
 347         }
 348         if (dd->nsystems > 1)
 349         {
 350             Rtl_6[res]   = Rt_6[res];
 351             Rt_6[res]   *= invn;
 352             Rtav_6[res] *= invn;
 353         }
 354
 355         res++;
 356     }
 357
 358 #ifdef GMX_MPI
 359     if (dd->nsystems > 1)
 360     {
 361         gmx_sum_sim(2*dd->nres,Rt_6,ms);
 362     }
 363 #endif
 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 lambda,real *dvdlambda,
 370                const t_mdatoms *md,t_fcdata *fcd,
 371                int *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 }