restrict to Rmax is now a parameter

TCL/polarDensity_for_DTA.tcl

@@ -327,10 +327,18 @@ proc leaflet_detector {atsel_in head tail frame_i leaflet_sorting_algorithm} {
 ;# Calculates the total number of lipids and beads of the given species in each leaflet 
 ;# Assigns the leaflet to user2 
 ;# Returns the following list : [["lower" lower_leaflet_beads lower_leaflet_lipids] ["upper" upper_leaflet_beads upper_leaflet_lipids]] 
-proc frame_leaflet_assignment {species headname tailname lipidbeads_selstr frame_i frame_f} {
+proc frame_leaflet_assignment {species headname tailname lipidbeads_selstr frame_i frame_f {restrict_to_Rmax 0}} {
     global params
-    set outer_r2 [expr $params(Rmax)**2]
-    set sel [ atomselect top "(($species)) and $lipidbeads_selstr and same resid as ((x*x + y*y < $outer_r2))"  frame $frame_i]
+    if {$restrict_to_Rmax == 1} {
+        set outer_r2 [expr $params(Rmax)**2]
+        set sel [ atomselect top "(($species)) and $lipidbeads_selstr and same resid as ((x*x + y*y < $outer_r2))"  frame $frame_i]
+    } elseif {$restrict_to_Rmax == 0} {
+        set sel [ atomselect top "(($species)) and $lipidbeads_selstr"  frame $frame_i]
+    } else {
+        puts "restrict_leaflet_sorter_to_Rmax must be 1 or 0."
+        error
+    }
+
     set sel_num [llength [lsort -unique [$sel get resid] ] ]
     set sel_resid_list [lsort -unique [$sel get resid] ]
     set totals {}
@@ -368,7 +376,7 @@ proc trajectory_leaflet_assignment {species headname tailname lipidbeads_selstr}
     global params
     set num_reassignments 0
     for {set update_frame $params(start_frame)} {$update_frame < $params(end_frame)} {incr update_frame $params(dt)} {
-        frame_leaflet_assignment $species $headname $tailname $lipidbeads_selstr $update_frame [expr $update_frame + $params(dt)] 
+        frame_leaflet_assignment $species $headname $tailname $lipidbeads_selstr $update_frame [expr $update_frame + $params(dt)] $params(restrict_leaflet_sorter_to_Rmax)
         incr num_reassignments
     }
     puts "Checked for leaflet reassignments $num_reassignments times."
@@ -534,7 +542,8 @@ proc set_parameters { config_file_script } {
         Rmax 20.
         Rmin 0.
         dr 5 
-        Ntheta 50        
+        Ntheta 50
+        restrict_leaflet_sorter_to_Rmax 0
     }
     set nframes [molinfo top get numframes] 
     array set params [list end_frame $nframes]
@@ -625,7 +634,7 @@ proc polarDensityBin { config_file_script } {
         set upp_f [open "${outfile}.upp.dat" w]
         set low_f_avg [open "${outfile}.low.avg.dat" w]
         set upp_f_avg [open "${outfile}.upp.avg.dat" w]
-        set totals [frame_leaflet_assignment "resname $species" $headname $tailname $lipidbeads_selstr $params(start_frame) $params(start_frame)]
+        set totals [frame_leaflet_assignment "resname $species" $headname $tailname $lipidbeads_selstr $params(end_frame) $params(end_frame)]        
 
         foreach lu [list $low_f $upp_f] avgfile [list $low_f_avg $upp_f_avg] leaf_total $totals {
             set leaflet_str [lindex $leaf_total 0]
@@ -637,7 +646,7 @@ proc polarDensityBin { config_file_script } {
                 puts $avgfile "#Lipid species $species in $leaflet_str leaflet: ${expected_lipids} molecules, Num beads : ${expected_beads} beads,  Average Area : [format {%0.0f} $area] A^2, Expected Bead Density : [format {%0.5f} [expr $expected_bead_density]]/A^2, Average Chain : [avg_acyl_chain_len "resname $species" $acylchain_selstr] beads, dr*dtheta : [format {%0.5f} [expr $params(dr)*[DtoR $params(dtheta)]]] "
         }
         puts "Processing frames, starting at frame $params(start_frame) and ending at frame $params(end_frame)."
-        trajectory_leaflet_assignment "resname $species" $headname $tailname $lipidbeads_selstr         
+        trajectory_leaflet_assignment "resname $species" $headname $tailname $lipidbeads_selstr   
         ;#the core calculation 
         loop_over_shells $species $headname $tailname $lipidbeads_selstr $low_f $upp_f $low_f_avg $upp_f_avg