From 05db4892934e524633f0c25cbafd326d508d1632 Mon Sep 17 00:00:00 2001 From: Jesse Sandberg Date: Fri, 6 Dec 2024 14:31:08 -0500 Subject: [PATCH 01/10] leaflet sorter now only runs on lipids within rmax --- TCL/polarDensity_for_DTA.tcl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/TCL/polarDensity_for_DTA.tcl b/TCL/polarDensity_for_DTA.tcl index 9018a52..2f8bec4 100644 --- a/TCL/polarDensity_for_DTA.tcl +++ b/TCL/polarDensity_for_DTA.tcl @@ -329,7 +329,8 @@ proc leaflet_detector {atsel_in head tail frame_i leaflet_sorting_algorithm} { ;# 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} { global params - set sel [ atomselect top "(($species)) and $lipidbeads_selstr" frame $frame_i] + set outer_r2 [expr $params(Rmax)**2] + set sel [ atomselect top "(($species)) and $lipidbeads_selstr and ((x*x + y*y < $outer_r2))" frame $frame_i] set sel_num [llength [lsort -unique [$sel get resid] ] ] set sel_resid_list [lsort -unique [$sel get resid] ] set totals {} From 84e34e222f876edaae44e161117a0a3e464dd914 Mon Sep 17 00:00:00 2001 From: Jesse Sandberg Date: Mon, 9 Dec 2024 10:39:42 -0500 Subject: [PATCH 02/10] fixed issue 104 --- TCL/polarDensity_for_DTA.tcl | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/TCL/polarDensity_for_DTA.tcl b/TCL/polarDensity_for_DTA.tcl index 2f8bec4..9d7d5d1 100644 --- a/TCL/polarDensity_for_DTA.tcl +++ b/TCL/polarDensity_for_DTA.tcl @@ -330,7 +330,7 @@ proc leaflet_detector {atsel_in head tail frame_i leaflet_sorting_algorithm} { proc frame_leaflet_assignment {species headname tailname lipidbeads_selstr frame_i frame_f} { global params set outer_r2 [expr $params(Rmax)**2] - set sel [ atomselect top "(($species)) and $lipidbeads_selstr and ((x*x + y*y < $outer_r2))" frame $frame_i] + set sel [ atomselect top "(($species)) and $lipidbeads_selstr and same resid as ((x*x + y*y < $outer_r2))" frame $frame_i] set sel_num [llength [lsort -unique [$sel get resid] ] ] set sel_resid_list [lsort -unique [$sel get resid] ] set totals {} @@ -487,7 +487,7 @@ proc loop_over_frames {shell species headname tailname lipidbeads_selstr start_f proc loop_over_shells {species headname tailname lipidbeads_selstr low_f upp_f low_f_avg upp_f_avg} { global params set delta_frame [expr ($params(end_frame) - $params(start_frame)) / $params(dt)] - for {set ri $params(Rmin)} { $ri<=$params(Rmax)} { set ri [expr $ri + $params(dr)]} { + for {set ri $params(Rmin)} { $ri<$params(Rmax)} { set ri [expr $ri + $params(dr)]} { #loop over shells puts "Now on shell {$ri [expr ${ri}+$params(dr)]}" set rf [expr $ri + $params(dr)] @@ -574,6 +574,12 @@ proc polarDensityBin { config_file_script } { global params set_parameters $config_file_script source $params(utils)/BinTools.tcl + set divisibility_test [expr [expr int([expr $params(Rmax) * 10000])] % [expr int([expr $params(dr) * 10000])]] + if {$divisibility_test != 0} { + puts "Rmax must be evenly divisible by dr." + puts "Exiting polarDensityBin early." + return + } if {$params(use_qwrap) == 1} {load $params(utils)/qwrap.so} set backbone_selstr $params(backbone_selstr) ;#only necessary for backwards compatibility set protein_selstr $params(protein_selstr) ;#only necessary for backwards compatibility From dc3c70cb60ff6e657bc571ca6403add409850654 Mon Sep 17 00:00:00 2001 From: Jesse Sandberg Date: Mon, 9 Dec 2024 10:47:00 -0500 Subject: [PATCH 03/10] added comment --- TCL/polarDensity_for_DTA.tcl | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/TCL/polarDensity_for_DTA.tcl b/TCL/polarDensity_for_DTA.tcl index 9d7d5d1..f13f7e0 100644 --- a/TCL/polarDensity_for_DTA.tcl +++ b/TCL/polarDensity_for_DTA.tcl @@ -574,6 +574,12 @@ proc polarDensityBin { config_file_script } { global params set_parameters $config_file_script source $params(utils)/BinTools.tcl + + ;# check to make sure Rmax is evenly divisible by dr. + ;# Cannot ensure that Rmax and dr are both ints, so divisibility_test + ;# multiplies both by 10,000 before converting to int as a reasonable + ;# precaution. Specifying Rmax and dr to greater than 5 decimal places is + ;# a questionable life choice. set divisibility_test [expr [expr int([expr $params(Rmax) * 10000])] % [expr int([expr $params(dr) * 10000])]] if {$divisibility_test != 0} { puts "Rmax must be evenly divisible by dr." From db3a187c72078adaf772a81857b90b34e525186c Mon Sep 17 00:00:00 2001 From: Jesse Sandberg Date: Mon, 9 Dec 2024 17:53:27 -0500 Subject: [PATCH 04/10] restrict to Rmax is now a parameter --- TCL/polarDensity_for_DTA.tcl | 23 ++++++++++++++++------- 1 file changed, 16 insertions(+), 7 deletions(-) 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 From c4a1177760eff131a3ec9ad1ca7f76faef6fd6e7 Mon Sep 17 00:00:00 2001 From: Jesse Sandberg Date: Tue, 10 Dec 2024 09:04:40 -0500 Subject: [PATCH 05/10] error message now complete --- TCL/polarDensity_for_DTA.tcl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/TCL/polarDensity_for_DTA.tcl b/TCL/polarDensity_for_DTA.tcl index 8a5bd7c..8f7f15a 100644 --- a/TCL/polarDensity_for_DTA.tcl +++ b/TCL/polarDensity_for_DTA.tcl @@ -335,8 +335,7 @@ proc frame_leaflet_assignment {species headname tailname lipidbeads_selstr frame } 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 + error "restrict_leaflet_sorter_to_Rmax must be 1 or 0" } set sel_num [llength [lsort -unique [$sel get resid] ] ] From a534ebb24a6bdba432c0e0e6941bfae641e303e1 Mon Sep 17 00:00:00 2001 From: Jesse Sandberg Date: Tue, 10 Dec 2024 09:05:47 -0500 Subject: [PATCH 06/10] switched puts to error --- TCL/polarDensity_for_DTA.tcl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/TCL/polarDensity_for_DTA.tcl b/TCL/polarDensity_for_DTA.tcl index 8f7f15a..fdabed8 100644 --- a/TCL/polarDensity_for_DTA.tcl +++ b/TCL/polarDensity_for_DTA.tcl @@ -590,10 +590,9 @@ proc polarDensityBin { config_file_script } { ;# a questionable life choice. set divisibility_test [expr [expr int([expr $params(Rmax) * 10000])] % [expr int([expr $params(dr) * 10000])]] if {$divisibility_test != 0} { - puts "Rmax must be evenly divisible by dr." - puts "Exiting polarDensityBin early." - return + error "Rmax must be evenly divisible by dr." } + if {$params(use_qwrap) == 1} {load $params(utils)/qwrap.so} set backbone_selstr $params(backbone_selstr) ;#only necessary for backwards compatibility set protein_selstr $params(protein_selstr) ;#only necessary for backwards compatibility From c95c7e4971fc0f6e10382228589fc11469a7fb4f Mon Sep 17 00:00:00 2001 From: Jesse Sandberg Date: Tue, 7 Jan 2025 14:59:36 -0500 Subject: [PATCH 07/10] fixed 0 0 bug --- TCL/polarDensity_for_DTA.tcl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/TCL/polarDensity_for_DTA.tcl b/TCL/polarDensity_for_DTA.tcl index fdabed8..758a0ca 100644 --- a/TCL/polarDensity_for_DTA.tcl +++ b/TCL/polarDensity_for_DTA.tcl @@ -337,12 +337,11 @@ proc frame_leaflet_assignment {species headname tailname lipidbeads_selstr frame } else { error "restrict_leaflet_sorter_to_Rmax must be 1 or 0" } - set sel_num [llength [lsort -unique [$sel get resid] ] ] set sel_resid_list [lsort -unique [$sel get resid] ] set totals {} if {$sel_num < 1} { - set totals [[list 0 0] [list 0 0]] + set totals [list "lower 0 0" "upper 0 0"] } else { #assign leaflets from $frame_i to user2 field of each bead for this species foreach sel_resid $sel_resid_list { From c7cd234199bd4603056dd19bd9742ff39f0933c6 Mon Sep 17 00:00:00 2001 From: Jesse Sandberg Date: Thu, 23 Jan 2025 13:14:12 -0500 Subject: [PATCH 08/10] fixed divisibility test per discussion with GB --- TCL/polarDensity_for_DTA.tcl | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/TCL/polarDensity_for_DTA.tcl b/TCL/polarDensity_for_DTA.tcl index 758a0ca..f152c98 100644 --- a/TCL/polarDensity_for_DTA.tcl +++ b/TCL/polarDensity_for_DTA.tcl @@ -583,12 +583,9 @@ proc polarDensityBin { config_file_script } { source $params(utils)/BinTools.tcl ;# check to make sure Rmax is evenly divisible by dr. - ;# Cannot ensure that Rmax and dr are both ints, so divisibility_test - ;# multiplies both by 10,000 before converting to int as a reasonable - ;# precaution. Specifying Rmax and dr to greater than 5 decimal places is - ;# a questionable life choice. - set divisibility_test [expr [expr int([expr $params(Rmax) * 10000])] % [expr int([expr $params(dr) * 10000])]] - if {$divisibility_test != 0} { + set divisibility_test [expr [expr double($params(Rmax))] % [expr double($params(dr))]] + set TOLERANCE [expr 10**-12] + if {$divisibility_test >= $TOLERANCE} { error "Rmax must be evenly divisible by dr." } From 8e99ee0adb9895bc56347b4d32c00eaf14b9050b Mon Sep 17 00:00:00 2001 From: Jesse Sandberg Date: Thu, 23 Jan 2025 13:58:31 -0500 Subject: [PATCH 09/10] made divisibility test its own proc --- TCL/polarDensity_for_DTA.tcl | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/TCL/polarDensity_for_DTA.tcl b/TCL/polarDensity_for_DTA.tcl index f152c98..d4cfdff 100644 --- a/TCL/polarDensity_for_DTA.tcl +++ b/TCL/polarDensity_for_DTA.tcl @@ -394,6 +394,19 @@ proc clean_leaflet_assignments {species lipidbeads_selstr} { $sel delete } +#test to see if two floats are evenly divisible. Return 1 if evenly divisible. +#Return 0 if not evenly divisible. +proc test_if_evenly_divisible {dividend divisor} { + set TOLERANCE [expr 10.0**-12] + set float_quotient [expr [expr double($dividend)] / [expr double($divisor)]] + set int_quotient [expr int($float_quotient)] + set diff [expr $float_quotient - $int_quotient] + if {$diff <= $TOLERANCE} { + return 1 + } else { + return 0 + } +} #write radial and theta bin output to file proc output_bins {fl ri rf bins} { @@ -583,9 +596,7 @@ proc polarDensityBin { config_file_script } { source $params(utils)/BinTools.tcl ;# check to make sure Rmax is evenly divisible by dr. - set divisibility_test [expr [expr double($params(Rmax))] % [expr double($params(dr))]] - set TOLERANCE [expr 10**-12] - if {$divisibility_test >= $TOLERANCE} { + if {[test_if_evenly_divisible $params(Rmax) $params(dr)] != 1} { error "Rmax must be evenly divisible by dr." } From 5cefeb3f57888802b11399a5b6367007143d327f Mon Sep 17 00:00:00 2001 From: Jesse Sandberg Date: Fri, 24 Jan 2025 16:24:36 -0500 Subject: [PATCH 10/10] removed superfluous expr and double --- TCL/polarDensity_for_DTA.tcl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/TCL/polarDensity_for_DTA.tcl b/TCL/polarDensity_for_DTA.tcl index d4cfdff..df31e79 100644 --- a/TCL/polarDensity_for_DTA.tcl +++ b/TCL/polarDensity_for_DTA.tcl @@ -398,7 +398,7 @@ proc clean_leaflet_assignments {species lipidbeads_selstr} { #Return 0 if not evenly divisible. proc test_if_evenly_divisible {dividend divisor} { set TOLERANCE [expr 10.0**-12] - set float_quotient [expr [expr double($dividend)] / [expr double($divisor)]] + set float_quotient [expr $dividend / double($divisor)] set int_quotient [expr int($float_quotient)] set diff [expr $float_quotient - $int_quotient] if {$diff <= $TOLERANCE} {