src/gmxlib/disre.c

   1 /*
   2  * $Id$
   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 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * Green Red Orange Magenta Azure Cyan Skyblue
  28  */
  29 static char *SRCID_disre_c = "$Id$";
  30
  31 #include <math.h>
  32 #include "typedefs.h"
  33 #include "sysstuff.h"
  34 #include "smalloc.h"
  35 #include "macros.h"
  36 #include "vec.h"
  37 #include "futil.h"
  38 #include "xvgr.h"
  39 #include "fatal.h"
  40 #include "bondf.h"
  41 #include "copyrite.h"
  42 #include "disre.h"
  43
  44 #include "main.h"
  45
  46 /* Some local variables */
  47 static int  dr_weighting; /* Weighting of pairs in one restraint        */
  48 static bool dr_bMixed;    /* Use sqrt of the instantaneous times        *
  49                            * the time averaged violation                */
  50 static real dr_fc;        /* Force constant for disres,                 *
  51                            * which is multiplied by a (possibly)        *
  52                            * different factor for each restraint        */
  53 static real dr_tau;       /* Time constant for disres                   */
  54 static real ETerm,ETerm1; /* Exponential constants                      */
  55
  56 static t_drblock   drblock;
  57
  58 void init_disres(FILE *log,int nfa,t_inputrec *ir)
  59 {
  60   dr_weighting = ir->eDisreWeighting;
  61   dr_fc        = ir->dr_fc;
  62   dr_tau       = ir->dr_tau;
  63   if (dr_tau == 0.0) {
  64     dr_bMixed = FALSE;
  65     ETerm = 0.0;
  66   } else {
  67     dr_bMixed = ir->bDisreMixed;
  68     ETerm = exp(-(ir->delta_t/dr_tau));
  69   }
  70   ETerm1 = 1.0 - ETerm;
  71
  72   drblock.ndr = nfa/(interaction_function[F_DISRES].nratoms+1);
  73   snew(drblock.rav, drblock.ndr);
  74   snew(drblock.rt,  drblock.ndr);
  75
  76   if (drblock.ndr > 0) {
  77     fprintf(log,"Done init_disre, ndr=%d\n",drblock.ndr);
  78     please_cite(log,"Torda89a");
  79   }
  80 }
  81
  82 t_drblock *get_drblock(void)
  83 {
  84   return &drblock;
  85 }
  86
  87 real ta_disres(int nfa,t_iatom forceatoms[],t_iparams ip[],
  88                rvec x[],rvec f[],t_forcerec *fr,t_graph *g,
  89                matrix box,real lambda,real *dvdlambda,
  90                t_mdatoms *md,int ngrp,real egnb[],real egcoul[])
  91 {
  92 #define MAX_DRPAIRS 1000
  93   static real sixth=1.0/6.0;
  94   static real seven_three=7.0/3.0;
  95
  96   atom_id     ai,aj;
  97   int         fa,fa_start,i,pair,ki,kj,m;
  98   int         restraint,pair_nr,pair_nr_start,rindex,npairs,type;
  99   rvec        dx[MAX_DRPAIRS];
 100   real        weight_rt_1[MAX_DRPAIRS];
 101   rvec        *fshift;
 102   real        smooth_fc,rt,Rt,Rav,rav_3,rt_1,rt_3,Rav_6,Rt_6,rt2;
 103   real        k0,f_scal,fmax_scal,fk_scal,fij;
 104   real        tav_viol,instant_viol,mixed_viol,violtot;
 105   real        tav_viol_Rav7,instant_viol_Rav7;
 106   real        up1,up2,low;
 107   bool        bConservative,bMixed,bViolation;
 108   static real exp_min_t_tau=1.0;
 109
 110   tav_viol=instant_viol=mixed_viol=tav_viol_Rav7=instant_viol_Rav7=0;
 111   /* scaling factor to smoothly turn on the restraint forces *
 112    * when using time averaging                               */
 113   exp_min_t_tau *= ETerm;
 114   smooth_fc      = dr_fc * (1.0 - exp_min_t_tau);
 115
 116   fshift  = fr->fshift;
 117   violtot = 0;
 118   pair_nr = 0;
 119
 120   /* 'loop' over all atom pairs (pair_nr=fa/3) involved in restraints, *
 121    * the total number of atoms pairs is nfa/3                          */
 122   fa=0;
 123   while (fa < nfa) {
 124     restraint = forceatoms[fa];
 125     rindex    = ip[restraint].disres.index;
 126
 127     fa_start      = fa;
 128     pair_nr_start = pair_nr;
 129
 130     Rav_6  = 0.0;
 131     Rt_6   = 0.0;
 132
 133     /* Loop over the atom pairs of 'this' restraint */
 134     for(pair=0; (fa<nfa) && (ip[forceatoms[fa]].disres.index == rindex);
 135         pair++,pair_nr++,fa+=3) {
 136       if (pair >= MAX_DRPAIRS)
 137         fatal_error(0,"Too many distance restraint pairs");
 138
 139       ai   = forceatoms[fa+1];
 140       aj   = forceatoms[fa+2];
 141
 142       rvec_sub(x[ai],x[aj],dx[pair]);
 143       rt2  = iprod(dx[pair],dx[pair]);
 144       rt   = sqrt(rt2);
 145       rt_1 = invsqrt(rt2);
 146       rt_3 = rt_1*rt_1*rt_1;
 147       if (drblock.rav[pair_nr] == 0)
 148         rav_3 = rt_3;
 149       else
 150         rav_3 = ETerm*drblock.rav[pair_nr] + ETerm1*rt_3;
 151
 152       weight_rt_1[pair] = rt_1;
 153
 154       drblock.rt[pair_nr]  = rt;
 155       drblock.rav[pair_nr] = rav_3;
 156
 157       Rav_6  += rav_3*rav_3;
 158       Rt_6   += rt_3*rt_3;
 159     }
 160     npairs=pair;
 161
 162     /* Take action depending on restraint, calculate scalar force */
 163     type = ip[restraint].disres.type;
 164     up1  = ip[restraint].disres.up1;
 165     up2  = ip[restraint].disres.up2;
 166     low  = ip[restraint].disres.low;
 167     k0   = smooth_fc*ip[restraint].disres.fac;
 168
 169     /* save some flops when there is only one pair */
 170     if (type != 2) {
 171       bConservative = (dr_weighting == edrwConservative) && (npairs>1);
 172       bMixed        = dr_bMixed;
 173     } else {
 174       bConservative = npairs>1;
 175       bMixed        = FALSE;
 176     }
 177
 178     fmax_scal = -k0*(up2-up1);
 179
 180     Rt   = pow(Rt_6,-sixth);
 181     /* Compute rav from ensemble / restraint averaging */
 182     if (type != 2)
 183       Rav  = pow(Rav_6,-sixth);
 184     else
 185       /* Dirty trick to not use time averaging when type=2 */
 186       Rav  = Rt;
 187
 188     if (Rav > up1) {
 189       bViolation = TRUE;
 190       tav_viol = Rav - up1;
 191     } else if (Rav < low) {
 192       bViolation = TRUE;
 193       tav_viol = Rav - low;
 194     } else
 195       bViolation = FALSE;
 196     if (bViolation) {
 197       if (!bMixed) {
 198         f_scal   = -k0*tav_viol;
 199         violtot += fabs(tav_viol);
 200       } else {
 201         if (Rt > up1) {
 202           if (tav_viol > 0)
 203             instant_viol = Rt - up1;
 204           else
 205             bViolation = FALSE;
 206         } else if (Rt < low) {
 207           if (tav_viol < 0)
 208             instant_viol = Rt - low;
 209           else
 210             bViolation = FALSE;
 211         } else
 212           bViolation = FALSE;
 213         if (bViolation) {
 214           mixed_viol   = sqrt(tav_viol*instant_viol);
 215           f_scal       = -k0*mixed_viol;
 216           violtot     += mixed_viol;
 217         }
 218       }
 219     }
 220     /* save some flops and don't use variables when they are not calculated */
 221     if (bViolation) {
 222       /* Correct the force for the number of restraints */
 223       if (bConservative) {
 224         f_scal  = max(f_scal,fmax_scal);
 225         if (!bMixed)
 226           f_scal *= Rav/Rav_6;
 227         else {
 228           f_scal /= 2*mixed_viol;
 229           tav_viol_Rav7     = tav_viol*Rav/Rav_6;
 230           instant_viol_Rav7 = instant_viol*Rt/Rt_6;
 231         }
 232       } else {
 233         f_scal /= (real)npairs;
 234         f_scal  = max(f_scal,fmax_scal);
 235       }
 236
 237       /* Exert the force ... */
 238       i = fa_start;
 239
 240       for(pair=0; (pair<npairs); pair++) {
 241         restraint= forceatoms[i++];
 242         ai       = forceatoms[i++];
 243         aj       = forceatoms[i++];
 244         ki       = SHIFT_INDEX(g,ai);
 245         kj       = SHIFT_INDEX(g,aj);
 246
 247         if (bConservative) {
 248           rav_3 = drblock.rav[pair_nr_start+pair];
 249           if (!dr_bMixed)
 250             weight_rt_1[pair] *= pow(rav_3,seven_three);
 251           else {
 252             rt = drblock.rt[pair_nr_start+pair];
 253             weight_rt_1[pair] *= tav_viol_Rav7*pow(rav_3,seven_three)+
 254               instant_viol_Rav7*pow(rt,-7);
 255           }
 256         }
 257
 258         fk_scal  = f_scal*weight_rt_1[pair];
 259
 260         for(m=0; (m<DIM); m++) {
 261           fij            = fk_scal*dx[pair][m];
 262
 263           f[ai][m]      += fij;
 264           f[aj][m]      -= fij;
 265           fshift[ki][m] += fij;
 266           fshift[kj][m] -= fij;
 267         }
 268       }
 269     }
 270   }
 271   /* Return violation */
 272   return violtot;
 273 }
 274