VMD-L Mailing List

From: dipak sanap (dipak94sanap_at_gmail.com)
Date: Thu Oct 22 2015 - 01:06:39 CDT

Dear All,

I am running simulations which has 96 residue protein. I want to calculate
number of water molecules in hydration shell for carbonyl oxygen of each
residue over whole trajectory. I have written a script, however, it is
giving odd output. Please, any help would be appreciated.
Here is my script:

#!/bin/env wish

set carb [atomselect top "name O"]; #select all carbonyl oxygens for all
residues.

proc hydshell {sel} {

set wat [atomselect top "resname SOL and name OW"] ;#select all water
oxygens
set n [molinfo top get numframes]; #number of frames

for {set i 1} {$i <= $n} {incr i} {

$sel frame $i
$sel update
$wat frame $i
$wat update

foreach co $sel {

set nbrs [measure contacts 4 $wat $co]; #get neighbour water oxygens
around carbonyl oxygen of each residue
set nbrs_list [lindex $nbrs 0] ; #get list of neighbour oxygens
set nbrs_count [llength $nbrs_list]; #count number of neighbours
lappend l1 [list $nbrs_list,] ;#append residue name and neighbours
}
}
return $l1
}
set cn_file [open "cnonly_1-5.dat" w];
puts $cn_file "[hydshell $carb]";
close $cn_file