diff --git a/TCL/polarDensity_for_DTA.tcl b/TCL/polarDensity_for_DTA.tcl index f13f7e0..8a5bd7c 100644 --- a/TCL/polarDensity_for_DTA.tcl +++ b/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