Calculate mean residence time of water in MD simulation with AMBERTOOLS. Part3

I have published several blogs for the calculation of residence time of water. I thought the distribution of residence time should be ok for the analysis. But the advisor asked me to calculate the mean residence time. This requires the integration of the mean residence time of all chosen water molecules.

This is introduced at http://en.wikipedia.org/wiki/Residence_time.

 MRT = \frac{1}{N}\sum_{i=1}^m t_i n_i

where:

  •  is the total number of groups
  •  is the average time in the body
  •  is the number of molecules in the ith group
  •  is the total number of molecules introduced into the system

For the example, there are 1000 frames. Each frame is 1ps. I chose 100ps for one round calculation. Since there are 1000 frames. I did a set of loops from 1 to 100… to 901 to 1000. The output is the mean residence time in 100ps of one water molecule at each loop. And then I can use excel to get the mean and standard deviation.

name=`echo $1 | cut -f1 -d"."`
rm -f ${name}_total_rt.dat
sed -i /Frame/d $1
for p in {1..901};do
        ((q=$p+99))
        m=`cat $1 | cut -c -10 | grep -n " ${p} " | head -1 | cut -f1 -d ":"`
        n=`cat $1 | cut -c -10 | grep -n " ${q} " | tail -1 | cut -f1 -d ":"`
        a=(`sed -n ${m},${n}p $1 | cut -c 36- | sort -n | uniq -c | sort -nr -k1 | cut -c -7`)
        ((sum=0))
        for ((i=0;i<${#a[*]};i++));do
                ((c=${a[$i]}))
                ((sum=$sum + $c))
        done
        sum=`echo "$sum/${#a[*]}" | bc -l | cut -c -6`
        echo $sum >> ${name}_total_rt.dat 
        unset a 
done
Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s