VMD-L Mailing List

From: Josh Vermaas (vermaasj_at_msu.edu)
Date: Wed Sep 14 2022 - 08:05:29 CDT

Hi Ryan,

I think there are going to be a few problems here with what is written
so far. I think measure bond takes two indices as arguments, not
atomselections, and does so by default for all frames. This is what I
might do for the bondmeasure proc.

proc bondmeasure { name idx1 idx2 } {
     set fout [open "bond_$name.dat" w]
     set bonddistancelist [measure bond $idx1 $idx2]
     for {set i 0 } { $i < [molinfo top get numframe] } { incr i } {
         puts $fout "$i [lindex $bonddistancelist $i]"
     }
     $fout close
}

This will open a file, collect all the bond distances, and then write
them to a file in the format you specified. Now how do you get those
indices? Something like this would get you two lists, keeping in mind
that arrays in tcl are 0-indexed.

set sel [atomselect top "sequence SIGYGD and name O"]
set idx1list [lrange [$sel get index] 0 11]
set idx2list [lrange [$sel get index] 12 end]

This will get you two lists, idx1list, and idx2list, which should have
the elements you want in the order you wanted. You'd then have a small
foreach loop to go through all the pairs you were interested in.

set counter 0
foreach idx1 $idx1list idx2 $idx2list {
     bondmeasure "pair$counter" $idx1 $idx2
     incr counter
}

-Josh

On 9/13/22 6:25 PM, Ryan Woltz wrote:
> Dear community,
>
>         I apologize if this is a little beyond the capabilities of tcl
> but from the pages I've read it should be possible to do, I just need
> help with the final step of extracting atom indicies automatically. I
> have a series of MD simulations and I would like to write a script to
> calculate bond length between 6 consecutive backbone oxygens between
> two chain. There are four chains total so I need to do this twice. I'm
> using a script from online called myrmsd which has been great for that
> job so I'm using it as a starting point. This is the function I've
> written using that as a template.
>
> proc bondmeasure { frame name } {
> global ref sel all
> $all move [measure fit $sel $ref]
> set fout [open "bond_$name.dat" a+]
> puts "$frame: [measure bond $sel $ref]"
> puts $fout "$frame      [format "%f" [measure bond $sel $ref]]"
> close $fout
> }
>
> set ref [atomselect top "name CA and segname PROC and sequence SIGYGD"
> frame 0]
> set all [atomselect top all]
> set sel [atomselect top "name CA and segname PROG and sequence SIGYGD"]
> set nframes [molinfo top get numframes]
> for {set i 0} {$i <$nframes} {incr i} {
> animate goto $i
> bondmeasure $i "G_SF"
> }
>
>       If this works like the original it prints a very nice 2 column
> table that I can then import into a graphing script. However, the
> atomselect as you might notice isn't specific enough for a single
> atom. I think the easiest way is to extract the indicies of each atom
> which I can do with this. NOTE: I have four segid PROA, PROC, PROE,
> PROG so I could include segid to get a shorter list.
>
> atomselect top "sequence SIGYGD and name O" (this printed
> atomselect166729, which is unique identifier so might need to save
> this to a variable correct?)
> atomselect166729 list
> 3904 3923 3930 3951 3958 3970 11097 11116 11123 11144 11151 11163
> 18290 18309 18316 18337 18344 18356 25483 25502 25509 25530 25537 25549
>
>        My question is, is there a way to take this list and extract
> the indicies from this and plug into the atomselect of the $sel $ref
> variables above? From the above example I need to compare:
> 3904 18290
> 3923 18309
>
>       or elements 1 and 13, 2 and 14...12 and 24. If this is not
> possible I can alternatively include resids but it'd be a little
> bothersome as each protein being simulated is slightly different in
> length and I'd have to change this numver for each one and hope it is
> correct. Or if I get a student a script like this might be used
> incorrectly if the person is not familiar with my simulations like I
> am and know these resid changes.
>
>      I'm very much open to any other suggestions of writing bond
> lengths for every frame of 12 specific pairs of atoms to 12 separate
> text files that are in a two columns table format and a unique file
> name for each file.
>
> Thank you for any help you can give,
>
> Ryan 

-- 
Josh Vermaas
vermaasj_at_msu.edu
Assistant Professor, Plant Research Laboratory and Biochemistry and Molecular Biology
Michigan State University
vermaaslab.github.io