#! /usr/bin/perl -w
# 
require $ENV{'CYCLIZE_HOME'}."lib/cyclize_lib.pl";

# Build the cyclize parameters hash
%cyclize_parameters = (
	"whole_chains"     => 1e8,
	"nrad_stats"       => 1e7,
	"icalcs"           => 100,
	"radial_cutoff"    => 60,
	"axial_cutoff"     => 40,
	"torsional_cutoff" => 36,
	"nkeepers"         => 10,
);

# Set the default values for generating the DNA
$Temp       = 300;

# Build the parameter hashes defined below
&build_parameter_hashes; 

# Read in the sequence file
$file = "12a156.seq";
%seq = &read_sequence( $file );

# Assign B-form parameters to the whole thing
%seq = &fill_params( \%seq, \%BDNA );

# Find and set the A-tract parameters
%seq = &fill_Atract_params( \%seq, \%Atract, \%Atract_3, \%Atract_5 );

# Find the flex site
$flex_sequence = "ACGCCTATAAACGCCTATAAACGCCTTGCTC";
@flex_positions = &find_positions( \%seq, $flex_sequence );

# Set the flex site parameters
%seq = &fill_params( \%seq, \%FLEX, @flex_positions );

# &print_params( \%seq, "sequence.in" );
# exit;

for $repeat ( 1 ) {
	&Cyclize( \%seq, \%cyclize_parameters, "flex_example_$repeat" );
}

exit;

############################################################
# build_parameters_hashes
# 
# USAGE: 
# &build_parameter_hashes;
#
# -This subroutine initialize the defaults for the "parameter hashes" used
#  to construct the "sequence hash".  
# -A single parameter hash contains six keys: tilt, tilt_flex, roll,
#  roll_flex, twist and twist_flex.
############################################################
sub build_parameter_hashes {

	#####################################################
	# BDNA site parameter definitions
	$BDNA_tormod         = 2.4e-19; # torsional modulus (erg*cm)
	$BDNA_rpb            = 3.4e-8;  # rise/bp (cm)
	$BDNA_helical_repeat = 10.45;   # helical repeat (bp/360 degrees)

	%BDNA = (
		tilt       => 0,
		roll       => 0,
		twist      => 360/$BDNA_helical_repeat,
		tilt_flex  => 4.842,
		roll_flex  => 4.842,
		twist_flex => sqrt( $BDNA_rpb*$k*$Temp/$BDNA_tormod )*180/$PI,
		dz         => 3.4,
	);

	#####################################################
	# FLEX site parameter definitions
	$FLEX_tormod         = 2.4e-19; # torsional modulus (erg*cm)
	$FLEX_rpb            = 3.4e-8;  # rise/bp (cm)
	$FLEX_helical_repeat = 10.45;   # helical repeat (bp/360 degrees)

	%FLEX = (
		tilt       => 0,
		roll       => 0,
		twist      => 360/$FLEX_helical_repeat,
		tilt_flex  => 9.2,
		roll_flex  => 9.2,
		twist_flex => sqrt( $FLEX_rpb*$k*$Temp/$FLEX_tormod )*180/$PI,
		dz         => 3.4,
	);
	
	#####################################################
	# Atract site parameter definitions
	$Atract_tormod         = 2.4e-19; # torsional modulus (erg*cm)
	$Atract_rpb            = 3.4e-8;  # rise/bp (cm)
	$Atract_helical_repeat = 10.33;   # helical repeat (bp/360 degrees)

	# Default Atract values for middle of Atract region
	%Atract = (
		tilt       => 0,        
		roll       => 0,        
		twist      => 360/$Atract_helical_repeat,
		tilt_flex  => 4.842,
		roll_flex  => 4.842,
		twist_flex => sqrt( $Atract_rpb*$k*$Temp/$Atract_tormod )*180/$PI,
		dz         => 3.4,
	);
	$Atract{min} = 4;       # minimum number of A's in Atract
	$Atract{max} = 6;       # maximum number of A's in Atract

	# Default Atract values for 5' junction of Atract
	%Atract_5 = (
		tilt       => -7.7,
		roll       => 4.6,   
		twist      => 360/$Atract_helical_repeat,
		tilt_flex  => 4.842,
		roll_flex  => 4.842,
		twist_flex => sqrt( $Atract_rpb*$k*$Temp/$Atract_tormod )*180/$PI,
		dz         => 3.4,
	);

	# Default Atract values for 3' junction
	%Atract_3 = (
		tilt       => 9.7,
		roll       => 0,
		twist      => 360/$Atract_helical_repeat,
		tilt_flex  => 4.842,
		roll_flex  => 4.842,
		twist_flex => sqrt( $Atract_rpb*$k*$Temp/$Atract_tormod )*180/$PI,
		dz         => 3.4,
	);
}