# Script for autodetection of the protein atoms on the outside-facing  surface (lipid-facing surface of the membrane channel, atoms will be marked by 1 in occupancy field) and inside-facing surface (pore-facing surface of the channel or protein core, marked by -1). Channel is supposed to be oriented such that the direction normal to membrane plane will be along z axis.
# Andriy Anishkin (anishkin@icqmail.com) UMCP


set filelist {
C:/simulations/mscl_c-o-int/pull-tm1-tm2-ends-g10_cls-opn/Hydration_energy_profiles/temp/min-tgt_mscl_hourglass_c1.0_fr0000.pdb
C:/simulations/mscl_c-o-int/pull-tm1-tm2-ends-g10_cls-opn/Hydration_energy_profiles/temp/min-tgt_mscl_hourglass_c1.0_fr0050.pdb
C:/simulations/mscl_c-o-int/pull-tm1-tm2-ends-g10_cls-opn/Hydration_energy_profiles/temp/min-tgt_mscl_hourglass_c1.0_fr0100.pdb
}



proc mark {} {
	set prot [atomselect top all]
	$prot set occupancy 0
	
	set num_atoms [molinfo top get numatoms]
	#Search radius to estimate atom location relative to the other atoms
	#set search_radius 11
	set search_radius 9
	
	for { set prat 0 } { $prat < $num_atoms} { incr prat } {
	
	        set prats [atomselect top "index $prat"]
	        puts "Finished [expr {$prat*100/$num_atoms}]% of [molinfo top get name]"
	
	        #find local local protein
	        set lprot [atomselect top "within $search_radius of (index $prat)"]
	
	        #if there is some protein
	        if {[$lprot num] > 0} {
	           #find the center of the atom's local local protein
	           lassign [measure center $lprot] wx wy wz
	
	           #find the atom's coordinates
	           set ax [$prats get x]
	           set ay [$prats get y]
	
	           #if the local protein is closer to the center
	           if { [ expr {( ( $wx * $wx ) + ( $wy * $wy ) ) }] <= [ expr {( ( $ax * $ax ) + ( $ay * $ay ) ) }] } {
	              #this is outer atom
	              $prats set occupancy 1
	           } else {
	              #this is inner atom
	              $prats set occupancy -1
	           }
	        }
	        $prats delete
	        $lprot delete
	        unset ax ay wx wy wz
	}
	
	
	
	$prot writepdb  [lindex [molinfo top get name] 0]_O-IO.pdb

}

# molecule new {c:\Andrey\Mechanosensing\Articles\Grant_2004\MscL_s3-prefilter_closed-expanded-open_Hydration_Energy\mscl_s3-prefilter_closed.pdb} filetype pdb
# mark
# molecule top delete

# molecule new {c:\Andrey\Mechanosensing\Articles\Grant_2004\MscL_s3-prefilter_closed-expanded-open_Hydration_Energy\mscl_s3-prefilter_expanded.pdb} filetype pdb
# mark
# molecule top delete

# molecule new {c:\Andrey\Mechanosensing\Articles\Grant_2004\MscL_s3-prefilter_closed-expanded-open_Hydration_Energy\mscl_s3-prefilter_open.pdb} filetype pdb
# mark
# molecule top delete
foreach file $filelist {
	molecule new $file
	mark
	molecule delete [molinfo top]
}

exit