puts "tools for the analysis of the number of h-bonded neighbors"
puts "lubos vrbka, 2006"
puts ""
puts "defined procedures/functions are:"
puts "GetMolHbonds {distance angle first last} (returns array)"
puts "WriteMolHbondsFile {n_hbonds_array filename}"
puts "ReadMolHbondsFile {filename} (returns array)"
puts "AverageMolHbonds {n_hbonds_array average} (returns array)"
puts "SetMolHbondsUser {n_hbonds_array}"
puts ""
puts "check the source/docs for usage / more info"

#to check if array exists from the previous runs if yes it has to be deleted
set outfile [open test.txt w]
proc existence {variable} {
    upvar $variable testVar
    if { [info exists testVar] } {
	puts "$variable Exists"
	#array unset $variable *
	array unset $variable
	set $variable {}
    } else {
	puts "$variable Does Not Exist"
    }
}
existence inbulkArray

# the following global variables control the analysis
# identification of frames in the trajectory
# these are overwritten from the procedure GetMolHbonds/ReadMolHbondsFile
set CC_first_frame 0
set CC_last_frame 1
set CC_total_frames 1
set CC_mol_id 0

# you need to change the following to change the names, ...
# names for atom selections
set CC_water_resname "TIP3"
set CC_osm_resname "GLYC"


# variables for PBC handling and analysis
# name of the temporary file - you need to have a write access to its location!
set CC_tmpname "_temporary.pdb"
set CC_psfname "extended_withGlyc2M.psf"
set CC_topname "top_all27_prot_lipid_ret_csff_adam_spc.rtf"



# *********************************************************************************
# *********************************************************************************


  # note that first frame has number 0
  global total_osm
  global CC_first_frame
  global CC_last_frame
  global CC_total_frames
  global CC_mol_id
  global CC_water_resname
  global CC_tmpname
  global CC_psfname
  global CC_topname
  

  # get ID of the currently active molecule
  # if you need another molecule -> set it here
  set mol_ID [molinfo top]
  set CC_mol_id $mol_ID
  set total_osm [atomselect top "resname GLYC and name C1 "]
  # ************************************************************************

  # handle frame selection
  # set global variables to zero
  set CC_first_frame 0
  set CC_last_frame 0
  set CC_total_frames 0
  
  set totframes [molinfo top get numframes]
  set startframe 0
  set endframe [expr $totframes - 1]	
  # set the global variables
  set CC_first_frame $startframe
  set CC_last_frame $endframe
  set CC_total_frames $totframes

 

  # ************************************************************************

  # set some selections
  set wat [atomselect $mol_ID "resname $CC_water_resname"]
  set osm [atomselect $mol_ID "resname $CC_osm_resname"] 

  
  # repeat analysis for all required frames
  for {set i $startframe} {$i <= $endframe} {incr i} {
    # update selections for the original geometry and get some residue mapping
    $wat frame $i
    $wat update
    $osm frame $i
    $osm update
    set wat_res [$wat get resid]
    set wat_name [$wat get name]
    set osm_res [$osm get resid]
    set osm_name [$osm get name]
     
    set box [molinfo $mol_ID get {a b c}]

    # ************************************************************************
    # handle the PBC - create the images of the central unit cell
    # requires psfgen to work
    package require psfgen

    # read the psf with topology for our system
    # has to be created beforehand
    mol load psf $CC_psfname
    # read the topology (was used to construct the psf file)
    topology $CC_topname

    set psf_ID [molinfo top]

    set n 0
    set seglist {}

    # force segid original structure - if empty then it doesn't work
    # $wat set segid R$n

    # for all 27 different combinations of central cell shift
    foreach x [list 0.0 -1.0 1.0] {
      foreach y [list 0.0 -1.0 1.0] {
        foreach z [list 0.0 -1.0 1.0] {
          set vec [list [expr {$x * [lindex $box 0]} ] [expr {$y * [lindex $box 1]} ] [expr {$z * [lindex $box 2]}]]
          $wat moveby $vec
          $osm moveby $vec
          puts "shifting by $vec"
          segment R$n {
            first NONE
            last NONE
           
          }
	  
          lappend seglist R$n
	  puts "i'm here\n"
          foreach resid $wat_res name $wat_name pos [$wat get {x y z}] {
            coord R$n $resid $name $pos
          }
	  foreach resid $osm_res name $osm_name pos [$osm get {x y z}] {
            coord R$n $resid $name $pos
          }
          incr n
          $wat moveby [vecinvert $vec]
	  $osm moveby [vecinvert $vec]
	  # clean up
	  unset vec
        }
      }
    }

    # just to be sure and clean up
    unset seglist
 
    # write the structure to a temporary file
    writepdb $CC_tmpname
    # load the temporary file
    mol load pdb $CC_tmpname

    # ************************************************************************
    # do the analysis
	set j 1
	while {$j<30} {
   		set incutoff [atomselect top "resname GLYC and name C1 and same residue as (within $j of protein)"]
    		set inbulk [ expr [$total_osm num] - [$incutoff num] ]
		if { ![info exists inbulkArray($j)] } {
	    		set inbulkArray($j) $inbulk
	
		} else {
	    		set inbulkArray($j) [expr $inbulkArray($j)+$inbulk]
		} 

		incr j
	}



    # clean up
    resetpsf
    unset wat_res
    unset wat_name
    unset osm_res
    unset osm_name
    # remove the molecule with PBC and psf
    mol delete top
    mol delete $psf_ID

  # repeat for another frame
  }
  #to go over the inbulk array and delete the values to the number of frames
  foreach {dist count} [array get inbulkArray] {
   puts $outfile "dist: $dist Count: $count"
}

close $outfile
  # clean up
  $wat delete
  $osm delete
  
  
  puts "done"
 



# *********************************************************************************
# *********************************************************************************