From: Josh Vermaas (vermaas2_at_illinois.edu)
Date: Wed Aug 26 2015 - 08:04:56 CDT

Hi Mahad,

You are on the right track. Here is how I would do it:

set sel [atomselect top "resid $i and name O4' C1' N9 C4"]
if { [$sel num] != 4 } {
#This is a pyrimidine. Pyrimdines don't have a N9, so their selection
would contain 3 rather than 4 atoms.
$sel delete
set sel [atomselect top "resid $i and name O4' C1' N1 C2"]
}
set idxlist [$sel get index] ; #These will be ordered based on the order
they appear in the pdb. Conviently, this happens to be the order you
want (I think. Please double check this).
set curval [measure dihedral [lindex $idxlist 0] [lindex $idxlist 1]
[lindex $idxlist 2] [lindex $idxlist 3]]
$sel delete

This eliminates many of the atomselections made within a loop, so in
principle it should run faster. It also deletes the atomselections it
makes, so it won't run into memory issues if you run this on a very
large number of frames.

-Josh Vermaas

On 08/26/2015 05:41 AM, Mahad Gatti wrote:
> Hi everyone,
>
> I'm writing a script to calculate a dihedral angle (chi torsional) in
> a series of dcd files, using big_dcd (working with RNA).
> So far so good, I was able to write all the instructions as you will
> see here:
>
> source bigdcd.tcl
> proc myrmsd { frame } {
>
> for {set i 636} {$i<=651} {incr i} { #this range is the one I want
>
> set fout1 [open "dihe_$i.dat" a+]
> set molid 0
>
> set a1 [[atomselect top "resid $i and name O4'"] get index]
> set a2 [[atomselect top "resid $i and name C1'"] get index]
> set a3 [[atomselect top "resid $i and name N1"] get index]
> set a4 [[atomselect top "resid $i and name C2"] get index]
>
> label delete Dihedrals
> set molid [molinfo top]
> set curframe [molinfo top get frame]
> molinfo top set frame $frame
> label add Dihedrals $molid/$a1 $molid/$a2 $molid/$a3 $molid/$a4
> # Get its value
> set curval [lindex [lindex [label list Dihedrals] 0] 4]
> molinfo top set frame $frame
> puts $fout1 [format "%f" $curval]
> close $fout1
> label delete Dihedrals
> }
>
> }
> set mol [mol new mutated.prmtop type parm7 waitfor all]
> # MFL DIST
>
> mol addfile mutated.inpcrd type rst7 waitfor all
> bigdcd myrmsd dcd test-1.dcd ....test-n.dcd
>
> The problem is related to the differences between pyrimidines and
> purines, because the atoms involved (to be more precise the
> nomenclature behind) are different, thus I wanted to add a flag like
> an if cycle to control and choose the proper selection for both cases.
>
> Pyrimidines = O4'-C1'-N1-C2
> Purines = O4'-C1'-N9-C4
>
> Thing is that I don't know how to use an if cycle in tcl that takes
> the returning value of a atomselect function (like [$a1 get resname] )
> and use that information for the right selection.
> My idea was to use the get resname command and check if the string is
> U/A/G/C, but I get errors regarding the syntax probably. I will try to
> write it down in pseudocode.
>
> If {[$a1 get resname] equals (U or C)}
> then atomselect O4'-C1'-N1-C2 (pyrimidines)
> else
> atomselect O4'-C1'-N9-C4 (purines)
> }
>
> Hope that I was clear enough.
>
> Thanks and regards.
>
> Mahad