Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Leaflet sorter speed up and shell selection bug fix #102

Merged
Changes from 6 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
34 changes: 27 additions & 7 deletions TCL/polarDensity_for_DTA.tcl
Original file line number Diff line number Diff line change
Expand Up @@ -327,9 +327,17 @@ 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} {
JesseSandberg marked this conversation as resolved.
Show resolved Hide resolved
proc frame_leaflet_assignment {species headname tailname lipidbeads_selstr frame_i frame_f {restrict_to_Rmax 0}} {
global params
set sel [ atomselect top "(($species)) and $lipidbeads_selstr" frame $frame_i]
JesseSandberg marked this conversation as resolved.
Show resolved Hide resolved
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 {
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 {}
Expand Down Expand Up @@ -367,7 +375,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)]
JesseSandberg marked this conversation as resolved.
Show resolved Hide resolved
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."
Expand Down Expand Up @@ -486,7 +494,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)]} {
JesseSandberg marked this conversation as resolved.
Show resolved Hide resolved
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)]
Expand Down Expand Up @@ -533,7 +541,8 @@ proc set_parameters { config_file_script } {
Rmax 20.
Rmin 0.
dr 5
Ntheta 50
JesseSandberg marked this conversation as resolved.
Show resolved Hide resolved
Ntheta 50
restrict_leaflet_sorter_to_Rmax 0
}
set nframes [molinfo top get numframes]
array set params [list end_frame $nframes]
Expand Down Expand Up @@ -573,6 +582,17 @@ 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.
JesseSandberg marked this conversation as resolved.
Show resolved Hide resolved
set divisibility_test [expr [expr int([expr $params(Rmax) * 10000])] % [expr int([expr $params(dr) * 10000])]]
if {$divisibility_test != 0} {
error "Rmax must be evenly divisible by dr."
}
Copy link
Member

@gbrannigan gbrannigan Jan 23, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JesseSandberg

This code is more complicated than it needs to be. Do a simple (i.e. double($Rmax) % $dr) modulo to get the remainder, and then compare the remainder against a small tolerance (not 0!).

Do not make the tolerance a magic number.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Code has been modified per discussion on slack. Thank you!


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
Expand Down Expand Up @@ -612,7 +632,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)]
JesseSandberg marked this conversation as resolved.
Show resolved Hide resolved
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]
Expand All @@ -624,7 +644,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
JesseSandberg marked this conversation as resolved.
Show resolved Hide resolved
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

Expand Down