An script for ACF

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 : Sun Dec 31 2017 - 23:20:32 CST