binaries/src/ViennaRNA/lib/params.c

   1 /*
   2
   3                   c Ivo Hofacker
   4
   5                   Vienna RNA package
   6 */
   7 #include <config.h>
   8 #include <stdio.h>
   9 #include <stdlib.h>
  10 #include <math.h>
  11 #include <string.h>
  12 #include "energy_par.h"
  13 #include "fold_vars.h"
  14 #include "utils.h"
  15 #include "params.h"
  16 /**
  17 *** \file params.c
  18 *** <P>
  19 *** This file provides functions that return temperature scaled energy parameters and
  20 *** Boltzmann weights packed in datastructures
  21 *** </P>
  22 ***/
  23
  24 /*@unused@*/
  25 static char rcsid[] UNUSED = "$Id: params.c,v 1.9 2008/07/04 14:29:14 ivo Exp $";
  26
  27 PRIVATE paramT p;
  28 PRIVATE int id=-1;
  29 /* variables for partition function */
  30 PRIVATE pf_paramT pf;
  31 PRIVATE int pf_id=-1;
  32
  33 #ifdef _OPENMP
  34 #pragma omp threadprivate(id, pf_id)
  35 #endif
  36
  37 PUBLIC paramT *scale_parameters(void){
  38   model_detailsT  md;
  39   set_model_details(&md);
  40   return get_scaled_parameters(temperature, md);
  41 }
  42
  43 PUBLIC paramT *get_scaled_parameters( double temp,
  44                                       model_detailsT md){
  45
  46   unsigned int i,j,k,l;
  47   double tempf;
  48   paramT *params;
  49
  50   params  = (paramT *)space(sizeof(paramT));
  51
  52   /* store the model details */
  53   params->model_details = md;
  54   params->temperature   = temp;
  55   tempf                 = ((params->temperature+K0)/Tmeasure);
  56
  57   for(i = VRNA_GQUAD_MIN_STACK_SIZE; i <= VRNA_GQUAD_MAX_STACK_SIZE; i++)
  58     for(j = 3*VRNA_GQUAD_MIN_LINKER_LENGTH; j <= 3*VRNA_GQUAD_MAX_LINKER_LENGTH; j++){
  59       double GQuadAlpha_T = (double)GQuadAlphadH - (double)(GQuadAlphadH - GQuadAlpha37) * tempf;
  60       double GQuadBeta_T = (double)GQuadBetadH - (double)(GQuadBetadH - GQuadBeta37) * tempf;
  61       params->gquad[i][j] = (int)GQuadAlpha_T*(i-1) + (int)(((double)GQuadBeta_T)*log(j - 2));
  62     }
  63
  64   for (i=0; i<31; i++)
  65     params->hairpin[i]  = hairpindH[i] - (hairpindH[i] - hairpin37[i])*tempf;
  66   for (i=0; i<=MIN2(30,MAXLOOP); i++) {
  67     params->bulge[i]          = bulgedH[i] - (bulgedH[i] - bulge37[i]) * tempf;
  68     params->internal_loop[i]  = internal_loopdH[i] - (internal_loopdH[i] - internal_loop37[i]) * tempf;
  69   }
  70   params->lxc = lxc37*tempf;
  71   for (; i<=MAXLOOP; i++) {
  72     params->bulge[i] = params->bulge[30]+(int)(params->lxc*log((double)(i)/30.));
  73     params->internal_loop[i] = params->internal_loop[30]+(int)(params->lxc*log((double)(i)/30.));
  74   }
  75
  76   params->ninio[2] = niniodH - (niniodH - ninio37) * tempf;
  77
  78   params->TripleC = TripleCdH - (TripleCdH - TripleC37) * tempf;
  79   params->MultipleCA = MultipleCAdH - (MultipleCAdH - MultipleCA37) * tempf;
  80   params->MultipleCB = MultipleCBdH - (MultipleCBdH - MultipleCB37) * tempf;
  81
  82   for (i=0; (i*7)<strlen(Tetraloops); i++)
  83     params->Tetraloop_E[i] = TetraloopdH[i] - (TetraloopdH[i]-Tetraloop37[i])*tempf;
  84   for (i=0; (i*5)<strlen(Triloops); i++)
  85     params->Triloop_E[i] =  TriloopdH[i] - (TriloopdH[i]-Triloop37[i])*tempf;
  86   for (i=0; (i*9)<strlen(Hexaloops); i++)
  87     params->Hexaloop_E[i] =  HexaloopdH[i] - (HexaloopdH[i]-Hexaloop37[i])*tempf;
  88
  89   params->TerminalAU = TerminalAUdH - (TerminalAUdH - TerminalAU37) * tempf;
  90
  91   params->DuplexInit = DuplexInitdH - (DuplexInitdH - DuplexInit37) *tempf;
  92
  93   params->MLbase = ML_BASEdH - (ML_BASEdH - ML_BASE37) * tempf;
  94
  95   for (i=0; i<=NBPAIRS; i++)
  96     params->MLintern[i] = ML_interndH - (ML_interndH - ML_intern37) * tempf;
  97
  98   params->MLclosing = ML_closingdH - (ML_closingdH - ML_closing37) * tempf;
  99
 100
 101   /* stacks    G(T) = H - [H - G(T0)]*T/T0 */
 102   for (i=0; i<=NBPAIRS; i++)
 103     for (j=0; j<=NBPAIRS; j++)
 104       params->stack[i][j] = stackdH[i][j] - (stackdH[i][j] - stack37[i][j])*tempf;
 105
 106   /* mismatches */
 107   for (i=0; i<=NBPAIRS; i++)
 108     for (j=0; j<5; j++)
 109       for (k=0; k<5; k++) {
 110         int mm;
 111         params->mismatchI[i][j][k]    = mismatchIdH[i][j][k] - (mismatchIdH[i][j][k] - mismatchI37[i][j][k])*tempf;
 112         params->mismatchH[i][j][k]    = mismatchHdH[i][j][k] - (mismatchHdH[i][j][k] - mismatchH37[i][j][k])*tempf;
 113         params->mismatch1nI[i][j][k]  = mismatch1nIdH[i][j][k]-(mismatch1nIdH[i][j][k]-mismatch1nI37[i][j][k])*tempf;/* interior nx1 loops */
 114         params->mismatch23I[i][j][k]  = mismatch23IdH[i][j][k]-(mismatch23IdH[i][j][k]-mismatch23I37[i][j][k])*tempf;/* interior 2x3 loops */
 115         if(md.dangles){
 116           mm                      = mismatchMdH[i][j][k] - (mismatchMdH[i][j][k] - mismatchM37[i][j][k])*tempf;
 117           params->mismatchM[i][j][k]    = (mm > 0) ? 0 : mm;
 118           mm                      = mismatchExtdH[i][j][k] - (mismatchExtdH[i][j][k] - mismatchExt37[i][j][k])*tempf;
 119           params->mismatchExt[i][j][k]  = (mm > 0) ? 0 : mm;
 120         }
 121         else{
 122           params->mismatchM[i][j][k] = params->mismatchExt[i][j][k] = 0;
 123         }
 124       }
 125
 126   /* dangles */
 127   for (i=0; i<=NBPAIRS; i++)
 128     for (j=0; j<5; j++) {
 129       int dd;
 130       dd = dangle5_dH[i][j] - (dangle5_dH[i][j] - dangle5_37[i][j])*tempf;
 131       params->dangle5[i][j] = (dd>0) ? 0 : dd;  /* must be <= 0 */
 132       dd = dangle3_dH[i][j] - (dangle3_dH[i][j] - dangle3_37[i][j])*tempf;
 133       params->dangle3[i][j] = (dd>0) ? 0 : dd;  /* must be <= 0 */
 134     }
 135   /* interior 1x1 loops */
 136   for (i=0; i<=NBPAIRS; i++)
 137     for (j=0; j<=NBPAIRS; j++)
 138       for (k=0; k<5; k++)
 139         for (l=0; l<5; l++)
 140           params->int11[i][j][k][l] = int11_dH[i][j][k][l] - (int11_dH[i][j][k][l] - int11_37[i][j][k][l])*tempf;
 141
 142   /* interior 2x1 loops */
 143   for (i=0; i<=NBPAIRS; i++)
 144     for (j=0; j<=NBPAIRS; j++)
 145       for (k=0; k<5; k++)
 146         for (l=0; l<5; l++) {
 147           int m;
 148           for (m=0; m<5; m++)
 149             params->int21[i][j][k][l][m] = int21_dH[i][j][k][l][m] - (int21_dH[i][j][k][l][m] - int21_37[i][j][k][l][m])*tempf;
 150         }
 151   /* interior 2x2 loops */
 152   for (i=0; i<=NBPAIRS; i++)
 153     for (j=0; j<=NBPAIRS; j++)
 154       for (k=0; k<5; k++)
 155         for (l=0; l<5; l++) {
 156           int m,n;
 157           for (m=0; m<5; m++)
 158             for (n=0; n<5; n++)
 159               params->int22[i][j][k][l][m][n] = int22_dH[i][j][k][l][m][n] - (int22_dH[i][j][k][l][m][n]-int22_37[i][j][k][l][m][n])*tempf;
 160         }
 161
 162   strncpy(params->Tetraloops, Tetraloops, 281);
 163   strncpy(params->Triloops, Triloops, 241);
 164   strncpy(params->Hexaloops, Hexaloops, 361);
 165
 166   params->id = ++id;
 167   return params;
 168 }
 169
 170
 171 /*------------------------------------------------------------------------*/
 172 #define SCALE 10
 173 /**
 174 *** dangling ends should never be destabilizing, i.e. expdangle>=1<BR>
 175 *** specific heat needs smooth function (2nd derivative)<BR>
 176 *** we use a*(sin(x+b)+1)^2, with a=2/(3*sqrt(3)), b=Pi/6-sqrt(3)/2,
 177 *** in the interval b<x<sqrt(3)/2
 178 */
 179 #define SMOOTH(X) ((X)/SCALE<-1.2283697)?0:(((X)/SCALE>0.8660254)?(X):\
 180           SCALE*0.38490018*(sin((X)/SCALE-0.34242663)+1)*(sin((X)/SCALE-0.34242663)+1))
 181
 182 /* #define SMOOTH(X) ((X)<0 ? 0 : (X)) */
 183
 184
 185 PUBLIC pf_paramT *get_scaled_pf_parameters(void){
 186   model_detailsT  md;
 187   set_model_details(&md);
 188   return get_boltzmann_factors(temperature, 1.0, md, pf_scale);
 189 }
 190
 191 PUBLIC pf_paramT *get_boltzmann_factors(double temp,
 192                                         double betaScale,
 193                                         model_detailsT md,
 194                                         double pf_scale){
 195
 196   unsigned  int i, j, k, l;
 197   double        kT, TT;
 198   double        GT;
 199   pf_paramT     *pf;
 200
 201   pf                = (pf_paramT *)space(sizeof(pf_paramT));
 202   pf->model_details = md;
 203   pf->temperature   = temp;
 204   pf->alpha         = betaScale;
 205   pf->kT = kT       = betaScale*(temp+K0)*GASCONST;   /* kT in cal/mol  */
 206   pf->pf_scale      = pf_scale;
 207   TT                = (temp+K0)/(Tmeasure);
 208
 209   for(i = VRNA_GQUAD_MIN_STACK_SIZE; i <= VRNA_GQUAD_MAX_STACK_SIZE; i++)
 210     for(j = 3*VRNA_GQUAD_MIN_LINKER_LENGTH; j <= 3*VRNA_GQUAD_MAX_LINKER_LENGTH; j++){
 211       double GQuadAlpha_T = (double)GQuadAlphadH - (double)(GQuadAlphadH - GQuadAlpha37) * TT;
 212       double GQuadBeta_T = (double)GQuadBetadH - (double)(GQuadBetadH - GQuadBeta37) * TT;
 213       GT = ((double)GQuadAlpha_T)*((double)(i-1)) + ((double)GQuadBeta_T)*log(((double)j) - 2.);
 214       pf->expgquad[i][j] = exp( -GT*10./kT);
 215     }
 216
 217   /* loop energies: hairpins, bulges, interior, mulit-loops */
 218   for (i=0; i<31; i++){
 219     GT  = hairpindH[i] - (hairpindH[i] - hairpin37[i])*TT;
 220     pf->exphairpin[i] = exp( -GT*10./kT);
 221   }
 222
 223   for (i=0; i<=MIN2(30, MAXLOOP); i++) {
 224     GT =  bulgedH[i]- (bulgedH[i] - bulge37[i])*TT;
 225     pf->expbulge[i] = exp( -GT*10./kT);
 226     GT =  internal_loopdH[i] - (internal_loopdH[i] - internal_loop37[i])*TT;
 227     pf->expinternal[i] = exp( -GT*10./kT);
 228   }
 229   /* special case of size 2 interior loops (single mismatch) */
 230   if (james_rule) pf->expinternal[2] = exp ( -80*10./kT);
 231
 232   pf->lxc = lxc37*TT;
 233
 234   GT =  DuplexInitdH - (DuplexInitdH - DuplexInit37)*TT;
 235   pf->expDuplexInit = exp( -GT*10./kT);
 236
 237   for (i=31; i<=MAXLOOP; i++) {
 238     GT = bulge37[30]*TT + (pf->lxc*log( i/30.));
 239     pf->expbulge[i] = exp( -GT*10./kT);
 240     GT = internal_loop37[30]*TT + (pf->lxc*log( i/30.));
 241     pf->expinternal[i] = exp( -GT*10./kT);
 242   }
 243
 244   GT = niniodH - (niniodH - ninio37)*TT;
 245   for (j=0; j<=MAXLOOP; j++)
 246       pf->expninio[2][j]=exp(-MIN2(MAX_NINIO,j*GT)*10./kT);
 247
 248   for (i=0; (i*7)<strlen(Tetraloops); i++) {
 249     GT = TetraloopdH[i] - (TetraloopdH[i]-Tetraloop37[i])*TT;
 250     pf->exptetra[i] = exp( -GT*10./kT);
 251   }
 252   for (i=0; (i*5)<strlen(Triloops); i++) {
 253     GT = TriloopdH[i] - (TriloopdH[i]-Triloop37[i])*TT;
 254     pf->exptri[i] = exp( -GT*10./kT);
 255   }
 256   for (i=0; (i*9)<strlen(Hexaloops); i++) {
 257     GT = HexaloopdH[i] - (HexaloopdH[i]-Hexaloop37[i])*TT;
 258     pf->exphex[i] = exp( -GT*10./kT);
 259   }
 260   GT =  ML_closingdH - (ML_closingdH - ML_closing37)*TT;
 261   pf->expMLclosing = exp( -GT*10./kT);
 262
 263   for (i=0; i<=NBPAIRS; i++) {
 264     GT =  ML_interndH - (ML_interndH - ML_intern37)*TT;
 265     /* if (i>2) GT += TerminalAU; */
 266     pf->expMLintern[i] = exp( -GT*10./kT);
 267   }
 268   GT = TerminalAUdH - (TerminalAUdH - TerminalAU37)*TT;
 269   pf->expTermAU = exp(-GT*10./kT);
 270
 271   GT = ML_BASEdH - (ML_BASEdH - ML_BASE37)*TT;
 272
 273   pf->expMLbase=exp(-10.*GT/kT);
 274
 275
 276   /* if dangles==0 just set their energy to 0,
 277      don't let dangle energies become > 0 (at large temps),
 278      but make sure go smoothly to 0                        */
 279   for (i=0; i<=NBPAIRS; i++)
 280     for (j=0; j<=4; j++) {
 281       if (md.dangles) {
 282         GT = dangle5_dH[i][j] - (dangle5_dH[i][j] - dangle5_37[i][j])*TT;
 283         pf->expdangle5[i][j] = exp(SMOOTH(-GT)*10./kT);
 284         GT = dangle3_dH[i][j] - (dangle3_dH[i][j] - dangle3_37[i][j])*TT;
 285         pf->expdangle3[i][j] =  exp(SMOOTH(-GT)*10./kT);
 286       } else
 287         pf->expdangle3[i][j] = pf->expdangle5[i][j] = 1;
 288     }
 289
 290   /* stacking energies */
 291   for (i=0; i<=NBPAIRS; i++)
 292     for (j=0; j<=NBPAIRS; j++) {
 293       GT =  stackdH[i][j] - (stackdH[i][j] - stack37[i][j])*TT;
 294       pf->expstack[i][j] = exp( -GT*10./kT);
 295     }
 296
 297   /* mismatch energies */
 298   for (i=0; i<=NBPAIRS; i++)
 299     for (j=0; j<5; j++)
 300       for (k=0; k<5; k++) {
 301         GT =  mismatchIdH[i][j][k] - ( mismatchIdH[i][j][k] - mismatchI37[i][j][k])*TT;
 302         pf->expmismatchI[i][j][k] = exp(-GT*10.0/kT);
 303         GT = mismatch1nIdH[i][j][k] - (mismatch1nIdH[i][j][k] - mismatch1nI37[i][j][k])*TT;
 304         pf->expmismatch1nI[i][j][k] = exp(-GT*10.0/kT);
 305         GT = mismatchHdH[i][j][k] - (mismatchHdH[i][j][k] - mismatchH37[i][j][k])*TT;
 306         pf->expmismatchH[i][j][k] = exp(-GT*10.0/kT);
 307         if (md.dangles) {
 308           GT = mismatchMdH[i][j][k] - (mismatchMdH[i][j][k] - mismatchM37[i][j][k])*TT;
 309           pf->expmismatchM[i][j][k] = exp(SMOOTH(-GT)*10.0/kT);
 310           GT = mismatchExtdH[i][j][k] - (mismatchExtdH[i][j][k] - mismatchExt37[i][j][k])*TT;
 311           pf->expmismatchExt[i][j][k] = exp(SMOOTH(-GT)*10.0/kT);
 312         }
 313         else{
 314           pf->expmismatchM[i][j][k] = pf->expmismatchExt[i][j][k] = 1.;
 315         }
 316         GT = mismatch23IdH[i][j][k] - (mismatch23IdH[i][j][k] - mismatch23I37[i][j][k])*TT;
 317         pf->expmismatch23I[i][j][k] = exp(-GT*10.0/kT);
 318       }
 319
 320   /* interior lops of length 2 */
 321   for (i=0; i<=NBPAIRS; i++)
 322     for (j=0; j<=NBPAIRS; j++)
 323       for (k=0; k<5; k++)
 324         for (l=0; l<5; l++) {
 325           GT = int11_dH[i][j][k][l] -
 326             (int11_dH[i][j][k][l] - int11_37[i][j][k][l])*TT;
 327           pf->expint11[i][j][k][l] = exp(-GT*10./kT);
 328         }
 329   /* interior 2x1 loops */
 330   for (i=0; i<=NBPAIRS; i++)
 331     for (j=0; j<=NBPAIRS; j++)
 332       for (k=0; k<5; k++)
 333         for (l=0; l<5; l++) {
 334           int m;
 335           for (m=0; m<5; m++) {
 336             GT = int21_dH[i][j][k][l][m] -
 337               (int21_dH[i][j][k][l][m] - int21_37[i][j][k][l][m])*TT;
 338             pf->expint21[i][j][k][l][m] = exp(-GT*10./kT);
 339           }
 340         }
 341
 342   /* interior 2x2 loops */
 343   for (i=0; i<=NBPAIRS; i++)
 344     for (j=0; j<=NBPAIRS; j++)
 345       for (k=0; k<5; k++)
 346         for (l=0; l<5; l++) {
 347           int m,n;
 348           for (m=0; m<5; m++)
 349             for (n=0; n<5; n++) {
 350               GT = int22_dH[i][j][k][l][m][n] -
 351                 (int22_dH[i][j][k][l][m][n]-int22_37[i][j][k][l][m][n])*TT;
 352               pf->expint22[i][j][k][l][m][n] = exp(-GT*10./kT);
 353             }
 354         }
 355
 356   strncpy(pf->Tetraloops, Tetraloops, 281);
 357   strncpy(pf->Triloops, Triloops, 241);
 358   strncpy(pf->Hexaloops, Hexaloops, 361);
 359
 360   return pf;
 361 }
 362
 363 PUBLIC pf_paramT *get_scaled_alipf_parameters(unsigned int n_seq){
 364   model_detailsT  md;
 365   set_model_details(&md);
 366   return get_boltzmann_factors_ali(n_seq, temperature, 1.0, md, pf_scale);
 367 }
 368
 369 PUBLIC pf_paramT *get_boltzmann_factors_ali(unsigned int n_seq,
 370                                             double temperature,
 371                                             double betaScale,
 372                                             model_detailsT md,
 373                                             double pf_scale){
 374
 375   /* scale energy parameters and pre-calculate Boltzmann weights */
 376   unsigned int  i, j, k, l;
 377   double        kTn, TT;
 378   double        GT;
 379   pf_paramT     *pf;
 380
 381   pf                = (pf_paramT *)space(sizeof(pf_paramT));
 382   pf->model_details = md;
 383   pf->alpha         = betaScale;
 384   pf->temperature   = temperature;
 385   pf->pf_scale      = pf_scale;
 386   pf->kT = kTn      = ((double)n_seq)*betaScale*(temperature+K0)*GASCONST;   /* kT in cal/mol  */
 387   TT                = (temperature+K0)/(Tmeasure);
 388
 389
 390    /* loop energies: hairpins, bulges, interior, mulit-loops */
 391   for (i=0; i<31; i++) {
 392     GT =  hairpindH[i] - (hairpindH[i] - hairpin37[i])*TT;
 393     pf->exphairpin[i] = exp( -GT*10./kTn);
 394   }
 395   /*add penalty for too short hairpins*/
 396   for (i=0; i<3; i++) {
 397     GT= 600/*Penalty*/*TT;
 398     pf->exphairpin[i] = exp( -GT*10./kTn);
 399   }
 400
 401   for (i=0; i<=MIN2(30, MAXLOOP); i++) {
 402     GT =  bulgedH[i]- (bulgedH[i] - bulge37[i])*TT;
 403     pf->expbulge[i] = exp( -GT*10./kTn);
 404     GT =  internal_loopdH[i] - (internal_loopdH[i] - internal_loop37[i])*TT;
 405     pf->expinternal[i] = exp( -GT*10./kTn);
 406   }
 407   /* special case of size 2 interior loops (single mismatch) */
 408   if (james_rule) pf->expinternal[2] = exp ( -80*10./kTn);
 409
 410   pf->lxc = lxc37*TT;
 411
 412   GT =  DuplexInitdH - (DuplexInitdH - DuplexInit37)*TT;
 413   pf->expDuplexInit = exp( -GT*10./kTn);
 414
 415   for (i=31; i<=MAXLOOP; i++) {
 416     GT = bulge37[30]*TT + (pf->lxc*log( i/30.));
 417     pf->expbulge[i] = exp( -GT*10./kTn);
 418     GT = internal_loop37[30]*TT + (pf->lxc*log( i/30.));
 419     pf->expinternal[i] = exp( -GT*10./kTn);
 420   }
 421
 422   GT = niniodH - (niniodH - ninio37)*TT;
 423   for (j=0; j<=MAXLOOP; j++)
 424     pf->expninio[2][j]=exp(-MIN2(MAX_NINIO,j*GT)*10./kTn);
 425
 426   for (i=0; (i*7)<strlen(Tetraloops); i++) {
 427     GT = TetraloopdH[i] - (TetraloopdH[i]-Tetraloop37[i])*TT;
 428     pf->exptetra[i] = exp( -GT*10./kTn);
 429   }
 430   for (i=0; (i*5)<strlen(Triloops); i++) {
 431     GT = TriloopdH[i] - (TriloopdH[i]-Triloop37[i])*TT;
 432     pf->exptri[i] = exp( -GT*10./kTn);
 433   }
 434   for (i=0; (i*9)<strlen(Hexaloops); i++) {
 435     GT = HexaloopdH[i] - (HexaloopdH[i]-Hexaloop37[i])*TT;
 436     pf->exphex[i] = exp( -GT*10./kTn);
 437   }
 438   GT =  ML_closingdH - (ML_closingdH - ML_closing37)*TT;
 439   pf->expMLclosing = exp( -GT*10./kTn);
 440
 441   for (i=0; i<=NBPAIRS; i++) { /* includes AU penalty */
 442     GT =  ML_interndH - (ML_interndH - ML_intern37)*TT;
 443     /* if (i>2) GT += TerminalAU; */
 444     pf->expMLintern[i] = exp( -GT*10./kTn);
 445   }
 446   GT = TerminalAUdH - (TerminalAUdH - TerminalAU37)*TT;
 447   pf->expTermAU = exp(-GT*10./kTn);
 448
 449   GT = ML_BASEdH - (ML_BASEdH - ML_BASE37)*TT;
 450   pf->expMLbase=exp(-10.*GT/(kTn/n_seq));
 451
 452
 453   /* if dangle_model==0 just set their energy to 0,
 454      don't let dangle energies become > 0 (at large temps),
 455      but make sure go smoothly to 0                        */
 456   for (i=0; i<=NBPAIRS; i++)
 457     for (j=0; j<=4; j++) {
 458       if (md.dangles) {
 459         GT = dangle5_dH[i][j] - (dangle5_dH[i][j] - dangle5_37[i][j])*TT;
 460         pf->expdangle5[i][j] = exp(SMOOTH(-GT)*10./kTn);
 461         GT = dangle3_dH[i][j] - (dangle3_dH[i][j] - dangle3_37[i][j])*TT;
 462         pf->expdangle3[i][j] =  exp(SMOOTH(-GT)*10./kTn);
 463       } else
 464         pf->expdangle3[i][j] = pf->expdangle5[i][j] = 1;
 465     }
 466
 467   /* stacking energies */
 468   for (i=0; i<=NBPAIRS; i++)
 469     for (j=0; j<=NBPAIRS; j++) {
 470       GT =  stackdH[i][j] - (stackdH[i][j] - stack37[i][j])*TT;
 471       pf->expstack[i][j] = exp( -GT*10./kTn);
 472     }
 473
 474   /* mismatch energies */
 475   for (i=0; i<=NBPAIRS; i++)
 476     for (j=0; j<5; j++)
 477       for (k=0; k<5; k++) {
 478         GT =  mismatchIdH[i][j][k] - ( mismatchIdH[i][j][k] - mismatchI37[i][j][k])*TT;
 479         pf->expmismatchI[i][j][k] = exp(-GT*10.0/kTn);
 480         GT = mismatch1nIdH[i][j][k] - (mismatch1nIdH[i][j][k] - mismatch1nI37[i][j][k])*TT;
 481         pf->expmismatch1nI[i][j][k] = exp(-GT*10.0/kTn);
 482         GT = mismatchHdH[i][j][k] - (mismatchHdH[i][j][k] - mismatchH37[i][j][k])*TT;
 483         pf->expmismatchH[i][j][k] = exp(-GT*10.0/kTn);
 484         if (md.dangles) {
 485           GT = mismatchMdH[i][j][k] - (mismatchMdH[i][j][k] - mismatchM37[i][j][k])*TT;
 486           pf->expmismatchM[i][j][k] = exp(SMOOTH(-GT)*10.0/kTn);
 487           GT = mismatchExtdH[i][j][k] - (mismatchExtdH[i][j][k] - mismatchExt37[i][j][k])*TT;
 488           pf->expmismatchExt[i][j][k] = exp(SMOOTH(-GT)*10.0/kTn);
 489         }
 490         else{
 491           pf->expmismatchM[i][j][k] = pf->expmismatchExt[i][j][k] = 1.;
 492         }
 493         GT = mismatch23IdH[i][j][k] - (mismatch23IdH[i][j][k] - mismatch23I37[i][j][k])*TT;
 494         pf->expmismatch23I[i][j][k] = exp(-GT*10.0/kTn);
 495       }
 496
 497
 498   /* interior lops of length 2 */
 499   for (i=0; i<=NBPAIRS; i++)
 500     for (j=0; j<=NBPAIRS; j++)
 501       for (k=0; k<5; k++)
 502         for (l=0; l<5; l++) {
 503           GT = int11_dH[i][j][k][l] -
 504             (int11_dH[i][j][k][l] - int11_37[i][j][k][l])*TT;
 505           pf->expint11[i][j][k][l] = exp(-GT*10./kTn);
 506         }
 507   /* interior 2x1 loops */
 508   for (i=0; i<=NBPAIRS; i++)
 509     for (j=0; j<=NBPAIRS; j++)
 510       for (k=0; k<5; k++)
 511         for (l=0; l<5; l++) {
 512           int m;
 513           for (m=0; m<5; m++) {
 514             GT = int21_dH[i][j][k][l][m] -
 515               (int21_dH[i][j][k][l][m] - int21_37[i][j][k][l][m])*TT;
 516             pf->expint21[i][j][k][l][m] = exp(-GT*10./kTn);
 517           }
 518         }
 519
 520   /* interior 2x2 loops */
 521   for (i=0; i<=NBPAIRS; i++)
 522     for (j=0; j<=NBPAIRS; j++)
 523       for (k=0; k<5; k++)
 524         for (l=0; l<5; l++) {
 525           int m,n;
 526           for (m=0; m<5; m++)
 527             for (n=0; n<5; n++) {
 528               GT = int22_dH[i][j][k][l][m][n] -
 529                 (int22_dH[i][j][k][l][m][n]-int22_37[i][j][k][l][m][n])*TT;
 530               pf->expint22[i][j][k][l][m][n] = exp(-GT*10./kTn);
 531             }
 532         }
 533
 534   strncpy(pf->Tetraloops, Tetraloops, 281);
 535   strncpy(pf->Triloops, Triloops, 241);
 536   strncpy(pf->Hexaloops, Hexaloops, 361);
 537
 538   return pf;
 539 }
 540
 541 PUBLIC pf_paramT *get_boltzmann_factor_copy(pf_paramT *par){
 542   pf_paramT *copy = NULL;
 543   if(par){
 544     copy = (pf_paramT *) space(sizeof(pf_paramT));
 545     memcpy(copy, par, sizeof(pf_paramT));
 546   }
 547   return copy;
 548 }
 549
 550 PUBLIC paramT *get_parameter_copy(paramT *par){
 551   paramT *copy = NULL;
 552   if(par){
 553     copy = (paramT *) space(sizeof(paramT));
 554     memcpy(copy, par, sizeof(paramT));
 555   }
 556   return copy;
 557 }
 558
 559 /*###########################################*/
 560 /*# deprecated functions below              #*/
 561 /*###########################################*/
 562
 563 PUBLIC paramT *copy_parameters(void){
 564   paramT *copy;
 565   if (p.id != id) scale_parameters();
 566   copy = (paramT *) space(sizeof(paramT));
 567   memcpy(copy, &p, sizeof(paramT));
 568   return copy;
 569 }
 570
 571 PUBLIC paramT *set_parameters(paramT *dest){
 572   memcpy(&p, dest, sizeof(paramT));
 573   return &p;
 574 }
 575
 576 PUBLIC pf_paramT *copy_pf_param(void){
 577   pf_paramT *copy;
 578   if (pf.id != pf_id) scale_pf_parameters();
 579   copy = (pf_paramT *) space(sizeof(pf_paramT));
 580   memcpy(copy, &pf, sizeof(pf_paramT));
 581   return copy;
 582 }
 583
 584 PUBLIC pf_paramT *set_pf_param(paramT *dest){
 585   memcpy(&pf, dest, sizeof(pf_paramT));
 586   return &pf;
 587 }
 588
 589 PUBLIC pf_paramT *scale_pf_parameters(void){
 590   return get_scaled_pf_parameters();
 591 #if 0
 592   /* scale energy parameters and pre-calculate Boltzmann weights */
 593   unsigned int i, j, k, l;
 594   double  kT, TT;
 595   double  GT;
 596
 597   /* scale pf_params() in partfunc.c is only a wrapper, that calls
 598      this functions !! */
 599
 600   pf.temperature = temperature;
 601   kT = (pf.temperature+K0)*GASCONST;   /* kT in cal/mol  */
 602   TT = (pf.temperature+K0)/(Tmeasure);
 603
 604    /* loop energies: hairpins, bulges, interior, mulit-loops */
 605   for (i=0; i<31; i++) {
 606     GT =  hairpin37[i]*TT;
 607     pf.exphairpin[i] = exp( -GT*10./kT);
 608   }
 609   for (i=0; i<=MIN2(30, MAXLOOP); i++) {
 610     GT =  bulge37[i]*TT;
 611     pf.expbulge[i] = exp( -GT*10./kT);
 612     GT =  internal_loop37[i]*TT;
 613     pf.expinternal[i] = exp( -GT*10./kT);
 614   }
 615   /* special case of size 2 interior loops (single mismatch) */
 616   if (james_rule) pf.expinternal[2] = exp ( -80*10./kT);
 617
 618   pf.lxc = lxc37*TT;
 619
 620   GT =  DuplexInitdH - (DuplexInitdH - DuplexInit37)*TT;
 621   pf.expDuplexInit = exp( -GT*10./kT);
 622
 623   for (i=31; i<=MAXLOOP; i++) {
 624     GT = bulge37[30]*TT + (pf.lxc*log( i/30.));
 625     pf.expbulge[i] = exp( -GT*10./kT);
 626     GT = internal_loop37[30]*TT + (pf.lxc*log( i/30.));
 627     pf.expinternal[i] = exp( -GT*10./kT);
 628   }
 629
 630   GT = niniodH - (niniodH - ninio37)*TT;
 631   for (j=0; j<=MAXLOOP; j++)
 632       pf.expninio[2][j]=exp(-MIN2(MAX_NINIO,j*GT)*10./kT);
 633
 634   for (i=0; (i*7)<strlen(Tetraloops); i++) {
 635     GT = TetraloopdH[i] - (TetraloopdH[i]-Tetraloop37[i])*TT;
 636     pf.exptetra[i] = exp( -GT*10./kT);
 637   }
 638   for (i=0; (i*5)<strlen(Triloops); i++) {
 639     GT = TriloopdH[i] - (TriloopdH[i]-Triloop37[i])*TT;
 640     pf.exptri[i] = exp( -GT*10./kT);
 641   }
 642   for (i=0; (i*9)<strlen(Hexaloops); i++) {
 643     GT = HexaloopdH[i] - (HexaloopdH[i]-Hexaloop37[i])*TT;
 644     pf.exphex[i] = exp( -GT*10./kT);
 645   }
 646   GT =  ML_closing37*TT;
 647   pf.expMLclosing = exp( -GT*10./kT);
 648
 649   for (i=0; i<=NBPAIRS; i++) { /* includes AU penalty */
 650     GT =  ML_intern37*TT;
 651     /* if (i>2) GT += TerminalAU; */
 652     pf.expMLintern[i] = exp( -GT*10./kT);
 653   }
 654   GT = TerminalAUdH - (TerminalAUdH - TerminalAU37)*TT;
 655   pf.expTermAU = exp(-GT*10./kT);
 656
 657   GT = ML_BASE37*TT;
 658   pf.expMLbase=exp(-10.*GT/kT);
 659
 660
 661   /* if dangle_model==0 just set their energy to 0,
 662      don't let dangle energies become > 0 (at large temps),
 663      but make sure go smoothly to 0                        */
 664   for (i=0; i<=NBPAIRS; i++)
 665     for (j=0; j<=4; j++) {
 666       if (dangles) {
 667         GT = dangle5_dH[i][j] - (dangle5_dH[i][j] - dangle5_37[i][j])*TT;
 668         pf.expdangle5[i][j] = exp(SMOOTH(-GT)*10./kT);
 669         GT = dangle3_dH[i][j] - (dangle3_dH[i][j] - dangle3_37[i][j])*TT;
 670         pf.expdangle3[i][j] =  exp(SMOOTH(-GT)*10./kT);
 671       } else
 672         pf.expdangle3[i][j] = pf.expdangle5[i][j] = 1;
 673       if (i>2) /* add TermAU penalty into dangle3 */
 674         pf.expdangle3[i][j] *= pf.expTermAU;
 675     }
 676
 677   /* stacking energies */
 678   for (i=0; i<=NBPAIRS; i++)
 679     for (j=0; j<=NBPAIRS; j++) {
 680       GT =  stackdH[i][j] - (stackdH[i][j] - stack37[i][j])*TT;
 681       pf.expstack[i][j] = exp( -GT*10./kT);
 682     }
 683
 684   /* mismatch energies */
 685   for (i=0; i<=NBPAIRS; i++)
 686     for (j=0; j<5; j++)
 687       for (k=0; k<5; k++) {
 688         GT =  mismatchIdH[i][j][k] - ( mismatchIdH[i][j][k] - mismatchI37[i][j][k])*TT;
 689         pf.expmismatchI[i][j][k] = exp(-GT*10./kT);
 690         GT = mismatch1nIdH[i][j][k] - (mismatch1nIdH[i][j][k] - mismatch1nI37[i][j][k])*TT;
 691         pf.expmismatch1nI[i][j][k] = exp(-GT*10./kT);
 692         GT = mismatchHdH[i][j][k] - (mismatchHdH[i][j][k] - mismatchH37[i][j][k])*TT;
 693         pf.expmismatchH[i][j][k] = exp(-GT*10./kT);
 694         GT = mismatch23IdH[i][j][k] - (mismatch23IdH[i][j][k] - mismatch23I37[i][j][k])*TT;
 695         pf.expmismatch23I[i][j][k] = exp(-GT*10./kT);
 696         if (dangles) {
 697           GT = mismatchMdH[i][j][k] - (mismatchMdH[i][j][k] - mismatchM37[i][j][k])*TT;
 698           pf.expmismatchM[i][j][k] = exp(-GT*10./kT);
 699           GT =  mismatchExtdH[i][j][k] - ( mismatchExtdH[i][j][k] - mismatchExt37[i][j][k])*TT;
 700           pf.expmismatchExt[i][j][k] = exp(-GT*10./kT);
 701         }
 702         else{
 703           pf.expmismatchM[i][j][k] = pf.expmismatchExt[i][j][k] = 1.;
 704         }
 705       }
 706
 707
 708   /* interior lops of length 2 */
 709   for (i=0; i<=NBPAIRS; i++)
 710     for (j=0; j<=NBPAIRS; j++)
 711       for (k=0; k<5; k++)
 712         for (l=0; l<5; l++) {
 713           GT = int11_dH[i][j][k][l] -
 714             (int11_dH[i][j][k][l] - int11_37[i][j][k][l])*TT;
 715           pf.expint11[i][j][k][l] = exp(-GT*10./kT);
 716         }
 717   /* interior 2x1 loops */
 718   for (i=0; i<=NBPAIRS; i++)
 719     for (j=0; j<=NBPAIRS; j++)
 720       for (k=0; k<5; k++)
 721         for (l=0; l<5; l++) {
 722           int m;
 723           for (m=0; m<5; m++) {
 724             GT = int21_dH[i][j][k][l][m] -
 725               (int21_dH[i][j][k][l][m] - int21_37[i][j][k][l][m])*TT;
 726             pf.expint21[i][j][k][l][m] = exp(-GT*10./kT);
 727           }
 728         }
 729
 730   /* interior 2x2 loops */
 731   for (i=0; i<=NBPAIRS; i++)
 732     for (j=0; j<=NBPAIRS; j++)
 733       for (k=0; k<5; k++)
 734         for (l=0; l<5; l++) {
 735           int m,n;
 736           for (m=0; m<5; m++)
 737             for (n=0; n<5; n++) {
 738               GT = int22_dH[i][j][k][l][m][n] -
 739                 (int22_dH[i][j][k][l][m][n]-int22_37[i][j][k][l][m][n])*TT;
 740               pf.expint22[i][j][k][l][m][n] = exp(-GT*10./kT);
 741             }
 742         }
 743
 744   strncpy(pf.Tetraloops, Tetraloops, 281);
 745   strncpy(pf.Triloops, Triloops, 241);
 746   strncpy(pf.Hexaloops, Hexaloops, 361);
 747
 748   pf.id = ++pf_id;
 749   return &pf;
 750 #endif
 751 }