From: faride badalkhani (farideh.khamseh_at_gmail.com)
Date: Tue Jun 07 2016 - 16:37:13 CDT
Dear NAMD experts,
I have performed a NPT simulation and plotted the Rg per time step. Now, I
need to plot autocorrelation function of the Rg as a function of time. I
prepared the following script, but the result is not reasonable. So, I
think something is wrong in the script. Could you help me, please?
The simulation includes 1019 time step.
set file [open rg.dat r]
while { [gets $file line] != -1 } {
lappend rglist $line
}
close $file
set lc [llength $rglist]
unset rglist
set data " "
set file [open rg.dat r]
while { [gets $file line] != -1 } {
set data "$data $line"
}
close $file
set endlag 1019
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 timestep($k) [lindex $data [expr 2*$k]]
set rg($k) [lindex $data [expr 2*$k+1]]
}
unset data
set avg_rg [avg rg]
for {set k 0} {$k < $lc} {incr k} {
set rg_adj($k) [expr $rg($k) - $avg_rg]
}
set file [open auto-corr.dat w]
for {set lag 0} {$lag <= $endlag} {incr lag} {
for {set k 0} {$k < [expr $lc-$lag]} {incr k} {
set data1($k) $rg_adj($k)
set data2($k) $rg_adj([expr $k+$lag])
set data2sq($k) [expr $data2($k) * $data2($k)]
set dataprod($k) [expr $data1($k) * $data2($k)]
}
set mean1 [avg data1]
set mean2 [avg data2sq]
set meanprod [avg dataprod]
# Calculate the Autocorrelation Function
puts $file "$lag\t [expr ($meanprod - $mean1*$mean1)/($mean2 -
$mean1*$mean1)]"
}
close $file
Regards,
Farideh
This archive was generated by hypermail 2.1.6 : Tue Dec 27 2016 - 23:22:15 CST