--- /dev/null
+/*
+ Heat Capacity of RNA molecule
+
+ c Ivo Hofacker and Peter Stadler
+ Vienna RNA package
+
+
+ calculates specific heat using C = - T d^2/dT^2 G(T)
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <ctype.h>
+#include <math.h>
+#include <unistd.h>
+#include "utils.h"
+#include "fold_vars.h"
+#include "fold.h"
+#include "params.h"
+#include "part_func.h"
+#include "read_epars.h"
+#include "RNAheat_cmdl.h"
+
+/*@unused@*/
+static char rcsid[] = "$Id: RNAheat.c,v 1.12 2000/09/28 11:23:14 ivo Rel $";
+
+#define MAXWIDTH 201
+#define GASCONST 1.98717 /* in [cal/K] */
+#define K0 273.15
+
+PRIVATE float F[MAXWIDTH];
+PRIVATE float ddiff(float f[], float h, int m);
+PRIVATE void heat_capacity(char *string, float T_min, float T_max, float h, int m);
+
+int main(int argc, char *argv[]){
+ struct RNAheat_args_info args_info;
+ char *string, *input_string, *ns_bases, *c, *ParamFile, *rec_sequence, *rec_id, **rec_rest, *orig_sequence;
+ int i, length, l, sym;
+ float T_min, T_max, h;
+ int mpoints, istty, noconv = 0;
+ unsigned int input_type;
+ unsigned int rec_type, read_opt;
+
+ string = ParamFile = ns_bases = NULL;
+ T_min = 0.;
+ T_max = 100.;
+ h = 1;
+ mpoints = 2;
+ dangles = 2; /* dangles can be 0 (no dangles) or 2, default is 2 */
+ rec_type = read_opt = 0;
+ rec_id = rec_sequence = orig_sequence = NULL;
+ rec_rest = NULL;
+
+ /*
+ #############################################
+ # check the command line parameters
+ #############################################
+ */
+ if(RNAheat_cmdline_parser(argc, argv, &args_info) != 0) exit(1);
+
+ /* do not take special tetra loop energies into account */
+ if(args_info.noTetra_given) tetra_loop=0;
+ /* set dangle model */
+ if(args_info.dangles_given){
+ if((args_info.dangles_arg != 0) && (args_info.dangles_arg != 2))
+ warn_user("required dangle model not implemented, falling back to default dangles=2");
+ else
+ dangles = args_info.dangles_arg;
+ }
+ /* do not allow weak pairs */
+ if(args_info.noLP_given) noLonelyPairs = 1;
+ /* do not allow wobble pairs (GU) */
+ if(args_info.noGU_given) noGU = 1;
+ /* do not allow weak closing pairs (AU,GU) */
+ if(args_info.noClosingGU_given) no_closingGU = 1;
+ /* do not convert DNA nucleotide "T" to appropriate RNA "U" */
+ if(args_info.noconv_given) noconv = 1;
+ /* set energy model */
+ if(args_info.energyModel_given) energy_set = args_info.energyModel_arg;
+ /* take another energy parameter set */
+ if(args_info.paramFile_given) ParamFile = strdup(args_info.paramFile_arg);
+ /* Allow other pairs in addition to the usual AU,GC,and GU pairs */
+ if(args_info.nsp_given) ns_bases = strdup(args_info.nsp_arg);
+ /* Tmin */
+ if(args_info.Tmin_given) T_min = args_info.Tmin_arg;
+ /* Tmax */
+ if(args_info.Tmax_given) T_max = args_info.Tmax_arg;
+ /* step size */
+ if(args_info.stepsize_given) h = args_info.stepsize_arg;
+ /* ipoints */
+ if(args_info.ipoints_given){
+ mpoints = args_info.ipoints_arg;
+ if (mpoints < 1) mpoints = 1;
+ if (mpoints > 100) mpoints = 100;
+ }
+
+ /* free allocated memory of command line data structure */
+ RNAheat_cmdline_parser_free (&args_info);
+
+ /*
+ #############################################
+ # begin initializing
+ #############################################
+ */
+ if (ParamFile!=NULL) read_parameter_file(ParamFile);
+
+ if (ns_bases != NULL) {
+ nonstandards = space(33);
+ c=ns_bases;
+ i=sym=0;
+ if (*c=='-') {
+ sym=1; c++;
+ }
+ while (*c!='\0') {
+ if (*c!=',') {
+ nonstandards[i++]=*c++;
+ nonstandards[i++]=*c;
+ if ((sym)&&(*c!=*(c-1))) {
+ nonstandards[i++]=*c;
+ nonstandards[i++]=*(c-1);
+ }
+ }
+ c++;
+ }
+ }
+
+ istty = isatty(fileno(stdout))&&isatty(fileno(stdin));
+
+ read_opt |= VRNA_INPUT_NO_REST;
+ if(istty){
+ print_tty_input_seq();
+ read_opt |= VRNA_INPUT_NOSKIP_BLANK_LINES;
+ }
+
+ /*
+ #############################################
+ # main loop: continue until end of file
+ #############################################
+ */
+ while(
+ !((rec_type = read_record(&rec_id, &rec_sequence, &rec_rest, read_opt))
+ & (VRNA_INPUT_ERROR | VRNA_INPUT_QUIT))){
+ /*
+ ########################################################
+ # init everything according to the data we've read
+ ########################################################
+ */
+ if(rec_id && !istty) printf("%s\n", rec_id);
+
+ length = (int)strlen(rec_sequence);
+
+ /* convert DNA alphabet to RNA if not explicitely switched off */
+ if(!noconv) str_DNA2RNA(rec_sequence);
+ /* store case-unmodified sequence */
+ orig_sequence = strdup(rec_sequence);
+ /* convert sequence to uppercase letters only */
+ str_uppercase(rec_sequence);
+
+ if(istty) printf("length = %d\n", length);
+ /*
+ ########################################################
+ # done with 'stdin' handling
+ ########################################################
+ */
+
+ heat_capacity(rec_sequence, T_min, T_max, h, mpoints);
+ (void) fflush(stdout);
+
+ /* clean up */
+ if(rec_id) free(rec_id);
+ free(rec_sequence);
+ free(orig_sequence);
+ rec_id = rec_sequence = orig_sequence = NULL;
+ rec_rest = NULL;
+ /* print user help for the next round if we get input from tty */
+
+ if(istty) print_tty_input_seq();
+ }
+ return EXIT_SUCCESS;
+}
+
+PRIVATE void heat_capacity(char *string, float T_min, float T_max,
+ float h, int m)
+{
+ int length, i;
+ char *structure;
+ float hc, kT, min_en;
+
+ length = (int) strlen(string);
+
+ do_backtrack = 0;
+
+ temperature = T_min -m*h;
+ /* initialize_fold(length); <- obsolete */
+ structure = (char *) space((unsigned) length+1);
+ min_en = fold(string, structure);
+ free(structure); free_arrays();
+ kT = (temperature+K0)*GASCONST/1000; /* in kcal */
+ pf_scale = exp(-(1.07*min_en)/kT/length );
+ /* init_pf_fold(length); <- obsolete */
+ pf_paramT *pf_parameters = NULL;
+ model_detailsT md;
+ set_model_details(&md);
+
+ for (i=0; i<2*m+1; i++) {
+ if(pf_parameters) free(pf_parameters);
+ pf_parameters = get_boltzmann_factors(temperature, 1.0, md, pf_scale);
+ F[i] = pf_fold_par(string, NULL, pf_parameters, 0, 0, 0); /* T_min -2h */
+ temperature += h;
+ kT = (temperature+K0)*GASCONST/1000;
+ pf_scale=exp(-(F[i]/length +h*0.00727)/kT); /* try to extrapolate F */
+ update_pf_params(length);
+ }
+ while (temperature <= (T_max+m*h+h)) {
+
+ hc = - ddiff(F,h,m)* (temperature +K0 - m*h -h);
+ printf("%g %g\n", (temperature-m*h-h), hc);
+
+ for (i=0; i<2*m; i++)
+ F[i] = F[i+1];
+
+ F[2*m] = pf_fold_par(string, NULL, pf_parameters, 0, 0, 0);
+/* printf("%g\n", F[2*m]);*/
+ temperature += h;
+ kT = (temperature+K0)*GASCONST/1000;
+ pf_scale=exp(-(F[i]/length +h*0.00727)/kT);
+ if(pf_parameters) free(pf_parameters);
+ pf_parameters = get_boltzmann_factors(temperature, 1.0, md, pf_scale);
+ update_pf_params(length);
+ }
+ free_pf_arrays();
+}
+
+/* ------------------------------------------------------------------------- */
+
+PRIVATE float ddiff(float f[], float h, int m)
+{
+ float fp;
+ int i;
+ float A, B;
+ A = (float)(m*(m+1)*(2*m+1)/3 ); /* 2*sum(x^2) */
+ B = (float)(m*(m+1)*(2*m+1) ) * (float)(3*m*m+3*m-1) /15.; /* 2*sum(x^4) */
+
+ fp=0.;
+ for (i=0; i<2*m+1; i++)
+ fp += f[i]*( A - (float) ( (2*m+1)*(i-m)*(i-m)) );
+
+ fp /= ( ( A*A - B*( (float)(2*m+1) ) )*h*h/2. );
+ return (float)fp;
+
+}
+
git://source.jalview.org / jabaws.git / blobdiff