proc HBONDING {CutOff} {
	set C {-0.110000 2.000000}
	set CA {-0.070000 1.992400}
	set CC {-0.070000 2.000000}
	set CD {-0.070000 2.000000}
	set CE1 {-0.068000 2.090000}
	set CE2 {-0.064000 2.080000}
	set CM {-0.110000 2.100000}
	set CP1 {-0.020000 2.275000}
	set CP2 {-0.055000 2.175000}
	set CP3 {-0.055000 2.175000}
	set CPA {-0.090000 1.800000}
	set CPB {-0.090000 1.800000}
	set CPH1 {-0.050000 1.800000}
	set CPH2 {-0.050000 1.800000}
	set CPM {-0.090000 1.800000}
	set CPT {-0.090000 1.800000}
	set CS {-0.110000 2.200000}
	set CT1 {-0.020000 2.275000}
	set CT2 {-0.055000 2.175000}
	set CT3 {-0.080000 2.060000}
	set CY {-0.070000 1.992400}
	set H {-0.046000 0.224500}
	set HA {-0.022000 1.320000}
	set HA1 {-0.031000 1.250000}
	set HA2 {-0.026000 1.260000}
	set HB {-0.022000 1.320000}
	set HC {-0.046000 0.224500}
	set HP {-0.030000 1.358200}
	set HR1 {-0.046000 0.900000}
	set HR2 {-0.046000 0.700000}
	set HR3 {-0.007800 1.468000}
	set HS {-0.100000 0.450000}
	set HT {-0.046000 0.224500}
	set N  {-0.200000 1.850000}
	set NC2 {-0.200000 1.850000}
	set NH1 {-0.200000 1.850000}
	set NH2 {-0.200000 1.850000}
	set NH3 {-0.200000 1.850000}
	set NP {-0.200000 1.850000}
	set NPH {-0.200000 1.850000}
	set NR1 {-0.200000 1.850000}
	set NR2 {-0.200000 1.850000}
	set NR3 {-0.200000 1.850000}
	set NY {-0.200000 1.850000}
	set O {-0.120000 1.700000}
	set OB {-0.120000 1.700000}
	set OC {-0.120000 1.700000}
	set OH1 {-0.152100 1.770000}
	set OM {-0.120000 1.700000}
	set OS {-0.152100 1.770000}
	set OT {-0.152100 1.768200}
	set CAL {-0.120000 1.710000}
	set FE { 0.000000 0.650000}
	set S {-0.450000 2.000000}
	set SM {-0.380000 1.975000}
	set SS {-0.470000 2.200000}
	set ZN {-0.250000 1.090000}
	set DUM {-0.000000 0.000000}
	set HE {-0.021270 1.4800}
	set NE {-0.086000 1.5300}
			  
	set ElectT 0.0
	set vdWT 0.0
	set HBET 0.0
	set NHBE 0
	set NI 0
	set ChainA [atomselect top "fragment 0"]
	set ChainB [atomselect top "fragment 1"]
	set all [atomselect top "all"]
	set types [$all get type]
	set masses [$all get mass]
	set charges [$all get charge]
	set positions [$all get {x y z}]
	set contacts [measure contacts $CutOff $ChainA $ChainB] 
	set atomsA [lindex $contacts 0] 
	set atomsB [lindex $contacts 1]
	set outfile [open energy.txt w];
	foreach A $atomsA B $atomsB {
		set NI [expr $NI + 1]
		set typeA [lindex $types $A]
		set typeB [lindex $types $B]
		set massA [lindex $masses $A]
		set massB [lindex $masses $B]
		set posA [lindex $positions $A]
		set posB [lindex $positions $B]
		set charA [lindex $charges $A]
		set charB [lindex $charges $B]
		set eA [lindex [expr $$typeA] 0]
		set eB [lindex [expr $$typeB] 0]
		set rA [lindex [expr $$typeA] 1]
		set rB [lindex [expr $$typeB] 1]
		set epsilon [expr sqrt($eA * $eB)]
		set rmin [expr $rA + $rB]
		set dist [vecdist $posA $posB]
		set rAB [expr ($rmin/$dist)*($rmin/$dist)*($rmin/$dist)*($rmin/$dist)*($rmin/$dist)*($rmin/$dist)]
		set Elect [expr $charA*$charB*55.1407113089098/$dist]
		set ElectT [expr $ElectT + $Elect]
		set vdW [expr $epsilon*($rAB*$rAB - 2*$rAB)]
		set vdWT [expr $vdWT + $vdW]
#Here I need the atom bonded to H in order to restrict HBond count with angle
		set HBE 0
		if {$rA == 0.224500 && $massB > 14.0 || $rB == 0.224500 && $massA > 14.0} {
			set NHBE [expr $NHBE + 1]
			set HBE [expr $charA*$charB*55.1407113089098/$dist + $epsilon*($rAB*$rAB - 2*$rAB)]
			set HBET [expr $HBET + $HBE]
		}
		puts $outfile "$NI $typeA - $typeB ($A - $B) $Elect $vdW $HBE"
	}
	close $outfile
	return [list $NI [expr $Elect*55.1407113089098] $vdW $NHBE $HBET]
}