VMD-L Mailing List

From: Jeremiah Babcock (zhc605_at_my.utsa.edu)
Date: Mon Apr 21 2014 - 12:32:30 CDT

My apologies on the late reply. I had to do some extensive testing of the
script and incorporate suggested amendments. Maxim, your advice was well
founded, as was Norman's and JC's. The problem was indeed the array to
function procedure. Thank you all for your help.

Regards,
Jeremiah B.

On Thu, Apr 17, 2014 at 1:36 PM, Maxim Belkin <mbelkin_at_ks.uiuc.edu> wrote:

> Hi,
> The problem might be with passing array to a function in tcl.
>
> - delete 'global M_PI' - I think you can use it without this line.
> - delete '$first update' and '$second update' from the loop, your
> selections do not change in time.
> - change 'set time [expr $j]' to 'set time $j'
> - change 'global $dataarray' in avg proc to 'upvar $dataarray localarray'
> and then in the code of avg procedure use 'localarray'.
> - again, in avg proc you can do
>
> foreach name [array names localarray] { set total [expr {$total +
> $localarray($name)}] }
>
> - later in the code:
> for {set k 0} {$k <= $lc} {incr k} {
> set angle($k) [lindex $data $k]
> }
>
> - use decimal points, e.g. change 3 to 3.0: 'set P2($k) [expr
> 0.5*(3.0*$muDot($k)*$muDot($k) - 1.0)]'
>
>
> Maxim
>
> On Apr 16, 2014, at 3:06 PM, Jeremiah Babcock <zhc605_at_my.utsa.edu> wrote:
>
> > Hello,
> >
> > I wrote a script for calculating the second order Legendre polynomial
> time correlation function for selected atomic vectors for the form <
> P2[u(t)*u(t+T)]>. My problem is that depending on the precision I select,
> my output is not dependable. For precision <=4, I get odd outputs for each
> frame:
> >
> > Precision 4:
> > 0.1433
> > 0.1433
> > 0.1433
> > 0.1433
> > ...
> >
> >
> > Precision 3:
> > 0.01433
> > 0.01433
> > 0.01433
> > 0.01433
> > ...
> >
> >
> > and so forth. For >4, I get good numbers, but regardless of which atoms
> I select, it always returns the same output values. This can not be
> correct, because some residues diffuse faster than others.
> >
> > I'm not a coding expert, but I've put a lot of time into researching
> this. I admit that I can't figure it out, so I'm moving on to Ambertools
> ptraj. However, I've posted this in the hope that maybe it can be saved
> and it can be used to compare outputs from ptraj. Thank you for your time.
> >
> > Jeremiah Babcock
> >
> >
> > #####################################################
> >
> > # Write a file where the coordinates of selected atoms are saved
> > set file [open ExcitationDipole.dat w]
> > set nfram [molinfo top get numframes]
> > set first [atomselect top "index 1"]
> > set second [atomselect top "index 2"]
> >
> > # Calculation of normalized vectors between #
> > # the selected atomic coordinates #
> >
> > set tcl_precision 5
> > global M_PI
> >
> > for {set j 0} {$j < $nfram } { incr j } {
> > $first frame $j
> > $second frame $j
> > $first update
> > $second update
> >
> > puts "Frame ............................ $j"
> > puts " "
> >
> > set time [expr $j]
> > set r1 [lindex [$first get {x y z}] 0]
> > set r2 [lindex [$second get {x y z}] 0]
> > set diff [vecsub $r2 $r1]
> > set mu [vecnorm $diff]
> > set type [list $mu]
> > puts $file "$type"
> > puts "$time $mu"
> > puts " "
> > }
> >
> > close $file
> >
> > # Reading coordinates
> >
> > set file [open ExcitationDipole.dat r]
> > while { [gets $file line] != -1 } {
> > lappend templist $line
> > }
> > close $file
> >
> > set lc [llength $templist]
> > unset templist
> >
> > set data " "
> > set file [open ExcitationDipole.dat r]
> > while { [gets $file line] != -1 } {
> > set data "$data $line"
> > }
> > close $file
> >
> > # $endlag represents the final coordinate frame step
> > # taken from the list length $lc
> > set endlag $lc
> > set lc [expr $lc -1]
> >
> > # procedure to calculate average of list
> >
> > proc avg {dataarray} {
> > global $dataarray
> > set numdata [array size $dataarray]
> > set total 0
> > for {set k 0} {$k < $numdata} {incr k} {
> > set total [expr "$[list $dataarray]($k)" + $total]
> > }
> > return [expr $total/$numdata]
> > }
> >
> > for {set k 0} {$k <= $lc} {incr k} {
> > set angle($k) [lindex $data [expr $k]]
> > }
> > unset data
> >
> > # Write a file with output values
> > set file1 [open Anisotropy.dat w]
> >
> > # Second order Legendre time correlation function
> > # calculation <P2[u(t)*u(t+T)]>
> >
> > for {set lag 0} {$lag <= $endlag} {incr lag} {
> >
> > for {set k 0} {$k < [expr $lc-$lag]} {incr k} {
> > set data1($k) $angle($k)
> > set data2($k) $angle([expr $k+$lag])
> > set muDot($k) [vecdot $data1($k) $data2($k)]
> > set P2($k) [expr (3*pow($muDot($k),2)-1)/2]
> > }
> > set meanP2 [avg P2]
> > puts "$lag $meanP2"
> > puts $file1 "$lag $meanP2"
> > }
> >
> > close $file1
> >
> > #####################################################
>
>