VMD-L Mailing List

From: Ryan Woltz (rlwoltz_at_ucdavis.edu)
Date: Tue Sep 13 2022 - 17:25:08 CDT

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