VMD-L Mailing List

From: Axel Kohlmeyer (akohlmey_at_gmail.com)
Date: Sun May 22 2022 - 12:24:04 CDT

Your problem seems to be the usual one that atom selections created with
the atomselect command persist until they go out of scope, but that will
not happen within loops.

So *everywhere* you have an [atomselect ....] you will need to do:
# at the beginning of the loop iteration
set sel1 [atomselect ...]
set sel2 [atomselect ...]
..
# now use $sel1 $sel2
..
# at the end of the loop iteration
$sel1 delete
$sel2 delete
..

If you are looping over a long trajectory and are already short on RAM, not
doing this can put you over the edge.

On Sun, May 22, 2022 at 12:55 PM Ackad, Edward <eackad_at_siue.edu> wrote:

> Hi all,
> I'm trying to write a script to renormalize the SASA measurement from GaMD
> runs (the scripts provided only do the reweighting for other measures). So
> I am trying to calculate the SASA for a small set of residues (using
> -restrict, which is why I can't use mdtraj's version). I can fit the
> stride=1 trajectory into memory without issue but once I start to calculate
> the SASA the memory blows up. I use a set for each frame and output to a
> file so I don't know why the memory requirements keep escalating. Is there
> a free or wipe memory option or some issue with my code (shown below)? Or
> any other ideas to make this not-so-memory intensive?
> Thanks!
>
>
> proc run_sasa {oname last_10ps_frame dcdstride } {
> # create the dictionary of the values for each reside
> set hydrophob [dict create ILE -9.73 ALA -9.57 PHE -9.23 LEU -8.66 MET
> -7.65 PRO -7.55 VAL -8.26 TRP -8.62 CYS -5.63 GLY -4.57 THR -4.06 SER -3.56
> TYR -0.05 GLN 2.2 ASN 2.76 HIS 1.41 >
>
> set rdis 1.4
> set fname "$oname.dat"
> set outfile [open $fname w]
> set molid top
> set nes "protein and resid 40 to 53"
> set nesphobic "protein and resid 40 to 53 and hydrophobic"
> set bucket "protein and resid 48 51"
>
> for {set frame 0} {$frame < [molinfo $molid get numframes]} {incr frame}
> {
> set printline [measure sasa $rdis [atomselect $molid "protein" frame
> $frame]]
> set printline2 [measure sasa $rdis [atomselect $molid "protein" frame
> $frame] -restrict [atomselect $molid $bucket frame $frame ] ]
> set printline3 [measure sasa $rdis [atomselect $molid "protein" frame
> $frame] -restrict [atomselect $molid $nes frame $frame ] ]
> set printline4 [measure sasa $rdis [atomselect $molid "protein" frame
> $frame] -restrict [atomselect $molid $nesphobic frame $frame ] ]
>
> # calculate the sasa's for each hydrophobic res
> set totalhydrophobic 0.0
> set totalsasa 0.0
> # loop over all the NES residues
> for { set i 40 } {$i < 54 } {incr i} {
> #use just one atom (alpha == name CA) so we only get 1 name back
> instead of a list
> set ressel [atomselect top "resid $i and alpha"]
> #determine what the coefficent is using the name of the residue as
> the index for the dictionary
> set coef [dict get $hydrophob [$ressel get resname]]
> set thisres [measure sasa $rdis [atomselect $molid "protein" frame
> $frame] -restrict [atomselect $molid "protein and resid $i" frame $frame ] ]
> # sum each residues value scaled by the coef
> set totalsasa [expr $totalhydrophobic + $thisres*$coef]
> # if the residue is hydrophoic (coef <0 then sub it so we can get a
> positive number)
> if {$coef <0 } {
> set totalhydrophobic [expr $totalhydrophobic - $thisres*$coef]
> }
> }
>
> #set the actual frame of the frame to match with the output in the log
> file
> if { $frame <= $last_10ps_frame} {
> set outframe [expr $dcdstride*$frame*10000]
> } else {
> set outframe [expr $dcdstride*$frame*100000]
> }
>
> puts $outfile "$outframe $printline $printline2 $printline3
> $printline4 $totalsasa $totalhydrophobic"
> }
> close $outfile
>
> }
>
> 

-- 
Dr. Axel Kohlmeyer  akohlmey_at_gmail.com  https://urldefense.com/v3/__http://goo.gl/1wk0__;!!DZ3fjg!7fc3SyhOcXB3vlxz4zHBTan1lZ7pvDKdrFIxnGeiiYFgtlwTaU0m42YCJNDxab2tcBIH0OVpVozFQzhvAg$ 
College of Science & Technology, Temple University, Philadelphia PA, USA
International Centre for Theoretical Physics, Trieste. Italy.