package require BigDCD

proc gauche {frame} {
    
global acum count id molid nangles gdefect DOT name1 num_steps
	set fid1 [open $name1 a+]
# Delete all existing Dihdral labels so that the one we add has index 0.
	label delete Dihedrals
	$DOT frame $frame
	$DOT update
	set atoms [$DOT get index]
	set atom_list [list]
	for {set j 0} {$j < [llength $atoms]} {incr j 32} {
		set temp_c [lrange $atoms $j [expr $j + 31]]
		lappend atom_list $temp_c
	}

	if {$frame < [expr $num_steps]} {
		if {[llength $atom_list] != 0} {
			for {set j 0} {$j < [llength $atom_list]} {incr j} { 
				for {set k 0} {$k < 29} {incr k} {
					set l [expr $k + 3]
					set Ci [lindex [lindex $atom_list $j] $k]
					set Cj [lindex [lindex $atom_list $j] [expr $k + 1]]
					set Ce [lindex [lindex $atom_list $j] [expr $l - 1]]
					set Cf [lindex [lindex $atom_list $j] $l]
				# Add the label
					label add Dihedrals $molid/$Ci $molid/$Cj $molid/$Ce $molid/$Cf
				# Get its value
					set angle [lindex [lindex [label list Dihedrals] 0] 4]

					label delete Dihedrals

				# determine if angle corresponds to a gauche defect
					if {[expr 50.0 < $angle]} {
						 if {[expr $angle < 70.0]} {
							set gdefect($k) [expr $gdefect($k) + 1.0]
						}
					}
					if {[expr -70.0 < $angle]} {
						 if {[expr $angle < -50.0]} {
							set gdefect($k) [expr $gdefect($k) + 1.0]
						}
					}
				}
			}
		}
	}
#
# compute time average
#

	for {set id 0} {$id < $nangles} {incr id} {
		set acum($id) [expr ($acum($id) * $count + $gdefect($id)) / ($count + 1)]
		puts $fid1 "$id $acum($id)"
		set gdefect($id)    0
	}
	set count [expr $count + 1]
	close $fid1
#	puts "[expr $num_steps - $frame] frames left."
}

mol load psf ../DOT_SiO.psf pdb ../DOT_SiO.pdb
# Use the top molecule
set molid [molinfo top]
set nangles 29
set count 0
set DOT [atomselect top "resname DOT and not hydrogen"]
set num_steps 2000

for {set id 0} {$id < $nangles} {incr id} {
	set acum($id)		0
	set gdefect($id)	0
}

bigdcd gauche 350K_19-20ns.dcd
set name1 "gauche_defects.dat"
set fid2 [open $name1 w]
close $fid2