VMD-L Mailing List

From: Nuno Loureiro Ferreira (nunolf_at_ci.uc.pt)
Date: Tue Jul 17 2007 - 02:57:40 CDT

Hi Sally

 From off of my head, this script should do the job.
Check it, it should work the way you want.

Notice that contact pairs are appended to pairsList var. If the number
of pairs is huge, probably you need more memory.
Or, instead of doing the -unique sorting after all contacts have been
measured (all traj frames) ,
you could update the uList var after each frame computation, but it will
slowww down the proc.
I'm still not very proficient on memory usage requirements, so give it a
try.

Cheers,
Nuno

#######################################
proc num_contact {mol sel1 sel2 dist fname} {
    puts "started at : [exec date]"
    set nf [molinfo top get numframes]
    set f [open $fname w]
    set ftmp [open ${fname}.tmp w]
    puts "measuring contacts ..."
    for {set i 0} {$i < $nf} {incr i} {
        $sel1 frame $i
        $sel1 update
        $sel2 frame $i
        $sel2 update
        set cList [measure contacts $dist $sel1 $sel2]
        foreach j [lindex $cList 0] k [lindex $cList 1] {
            lappend pairsList "$j $k" ;# list of contact pairs, frame
$i; adds up
            puts $ftmp "$j $k"
        }
    }
    close $ftmp
    puts "sorting pairs of contacts ..."
    set uList [lsort -unique $pairsList] ;# unique contact pairs
    # counter for each unique pair of contacts
    foreach item $uList {
        set elem1 [lindex $item 0]
        set elem2 [lindex $item 1]
        puts $f "contact $item : [lindex [exec grep "$elem1 $elem2"
${fname}.tmp | wc] 0] \toccurrence(s)"
    }
    close $f
    file delete ${fname}.tmp
    puts "ended at : [exec date]"
}

sally cii wrote:
> Dear all,
>
> Thank you for all your responses.
>
> I have $sel1 and $sel2 imported before I loaded this script. By
> correcting with set list1 [lsort -increasing $lists 0], I can get the
> list sorted and saved into the file $f.
>
> This script is small part of what I really need to do, what I really
> want is a script that does:
>
> For all the atoms in $sel1 there is a counter, every time an atom in
> $sel2 get within certain distance ($dist) to an atom in $sel1 during
> $nf frames, the corresponding counter to that atom adds 1. And the
> final output I want is the values of counter for every atom in $sel1.
> As I don't know how to accomplish that in one script, therefore I
> tried the previous script to sort the list of $sel1 and planned to
> count the index afterwards. I really appreciate if anyone can shed
> some light on how to do the counting in the script for me. Thank you.
>
> Best regards,
> Sally Cii
>
>
> On 7/14/07, *Nuno Loureiro Ferreira* <nunolf_at_ci.uc.pt
> <mailto:nunolf_at_ci.uc.pt>> wrote:
>
> Has Axel said, you should import those sel1 sel2 vars.
> But you say that you get the result on the console, so I deduce
> that you
> are feeding those vars some way to the proc.
>
> About the no output ... because your script reports an error before
> reaching the "puts $f".
> Specifically:
>
> set list1 [lsort -increasing $lists 0]
>
> What's the "0" doing there? ;-)
> And by the way, why are you sorting a list of 2 elements? Be sure you
> know what you are doing, since the results
> from the measure contacts command, are indices, you cannot sort these
> lists, since you the correct correspondence between the elements
> of both
> element lists.
>
> Probably you want to sort the pairs of contacts, using -increasing and
> -index 0?
> Check the following:
>
> vmd > set c [measure contacts 3 $sel1 $sel2]
> {25 25 21 22 21 22 34 22 22 34 34 35 35 23} {36 39 36 36 39 39 39
> 38 56
> 38 37 36 37 39}
> vmd > set a [lindex $c 0]
> 25 25 21 22 21 22 34 22 22 34 34 35 35 23
> vmd > set b [lindex $c 1]
> 36 39 36 36 39 39 39 38 56 38 37 36 37 39
> vmd > foreach i $a j $b {lappend nlist "$i $j"}
> vmd > puts $nlist
> {25 36} {25 39} {21 36} {22 36} {21 39} {22 39} {34 39} {22 38}
> {22 56}
> {34 38} {34 37} {35 36} {35 37} {23 39}
> vmd > lsort -increasing -index 0 $nlist
> {21 36} {21 39} {22 36} {22 39} {22 38} {22 56} {23 39} {25 36}
> {25 39}
> {34 39} {34 38} {34 37} {35 36} {35 37}
>
> Is this what you want?
>
> Cheers,
> Nuno
>
>
>
> sally cii wrote:
> > Dear all,
> > I want to write a tcl script that can return the list of atom
> index in
> > $sel1 that within a certain distance contact with atoms in $sel2
> > during a trajectory. Here is the trial script, after I load this
> > script, I got the atom index printed out in tcl console together
> with
> > complain about 'bad option', but not saved in the file I wanted.
> > Can anyone tell me how to return a list in a file? thank you in
> advance!
> >
> > #####################################################
> > proc num_contact {mol dist fname} {
> >
> > set nf [molinfo top get numframes]
> >
> > set f [open $fname "w"]
> >
> > for {set i 0} {$i < $nf} {incr i} {
> >
> > $sel1 frame $i
> >
> > $sel1 update
> >
> > $sel2 frame $i
> >
> > $sel2 update
> >
> > set lists [measure contacts $dist $sel1 $sel2]
> >
> > set list1 [lsort -increasing $lists 0]
> >
> > puts $f "$list1"}
> >
> > close $f
> >
> > $sel1 delete
> >
> > $sel2 delete
> >
> > }
> > ####################################################
> >
> > Sincerely,
> > Sally Cii
> >
> ------------------------------------------------------------------------
> >
> > No virus found in this incoming message.
> > Checked by AVG Free Edition.
> > Version: 7.5.476 / Virus Database: 269.10.4/898 - Release Date:
> 12-07-2007 16:08
> >
>
>
> ------------------------------------------------------------------------
>
> No virus found in this incoming message.
> Checked by AVG Free Edition.
> Version: 7.5.476 / Virus Database: 269.10.6/902 - Release Date: 15-07-2007 14:21
>