VMD-L Mailing List
From: Norman Geist (norman.geist_at_uni-greifswald.de)
Date: Thu Apr 17 2014 - 01:29:04 CDT
- Next message: Norge Cruz Hernández: "Re: convert AMBER format file to CHARMM format file or NAMD format file"
- Previous message: Jeremiah Babcock: "code precision problem"
- In reply to: Jeremiah Babcock: "code precision problem"
- Next in thread: JC Gumbart: "Re: code precision problem"
- Messages sorted by: [ date ] [ thread ] [ subject ] [ author ] [ attachment ]
You have likely done a failure while coding somewhere as changing precision,
so the number digits behind the comma, can't change the values like that:
0.1433 -> 0.01433
I guess you are using a rounding function somewhere like:
expr {double(round(100*$mynum))/100}
where the number of zeros for the "100" define the number of digits behind
the comma in the result. You might have done an error like:
expr {double(round(10*$mynum))/100}
--------------------------------!-----------------------
This would show you observed behavior and you might want to better do
something like:
proc prec_round {$num $prec} {
return [expr {double(round(1E$prec*$num))/1E$prec}]
}
set result [prec_round 1.3444 2]
will output 1.34
set result [prec_round 1.3444 4]
will output 1.3444
Norman Geist.
Von: owner-vmd-l_at_ks.uiuc.edu [mailto:owner-vmd-l_at_ks.uiuc.edu] Im Auftrag von
Jeremiah Babcock
Gesendet: Mittwoch, 16. April 2014 22:06
An: vmd-l_at_ks.uiuc.edu
Betreff: vmd-l: code precision problem
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
#####################################################
--- Diese E-Mail ist frei von Viren und Malware, denn der avast! Antivirus Schutz ist aktiv. http://www.avast.com
- Next message: Norge Cruz Hernández: "Re: convert AMBER format file to CHARMM format file or NAMD format file"
- Previous message: Jeremiah Babcock: "code precision problem"
- In reply to: Jeremiah Babcock: "code precision problem"
- Next in thread: JC Gumbart: "Re: code precision problem"
- Messages sorted by: [ date ] [ thread ] [ subject ] [ author ] [ attachment ]