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