# This script counts atoms within selection for all the frames of the trajectory and records results as a separate text files corresponding to the files in the filelist
# Also script extracts 1 parameter value (e.g. Z coordinate) for every atom in selection for every timestep
# Andriy Anishkin (anishkin@icqmail.com) UMCP


#Open file and load part of DCD trajectory
#mol new mscs-cry-chnl-big-oct-cor_wtr_f.psf type psf
mol new mscs-cry-chnl-big-oct-cor_wtr_f_salt.psf type psf
#mol new mscs-cry-chnl-big-oct-cor-l109s_wtr_f_salt.psf type psf
mol off top
set first_frame 0
set last_frame -1
#set filelist {sim_mscs-cry-chnl-big-oct-cor_wtr_f_wrp_0.5.dcd sim_mscs-cry-chnl-big-oct-cor_wtr_f_wrp_1.0.dcd sim_mscs-cry-chnl-big-oct-cor_wtr_f_wrp_1.5.dcd sim_mscs-cry-chnl-big-oct-cor_wtr_f_wrp_2.0.dcd sim_mscs-cry-chnl-big-oct-cor_wtr_f_wrp_2.5.dcd sim_mscs-cry-chnl-big-oct-cor_wtr_f_wrp_3.0.dcd sim_mscs-cry-chnl-big-oct-cor_wtr_f_wrp_3.5.dcd sim_mscs-cry-chnl-big-oct-cor_wtr_f_wrp_4.0.dcd sim_mscs-cry-chnl-big-oct-cor_wtr_f_wrp_4.5.dcd sim_mscs-cry-chnl-big-oct-cor_wtr_f_wrp_5.0.dcd sim_mscs-cry-chnl-big-oct-cor_wtr_f_wrp_5.5.dcd sim_mscs-cry-chnl-big-oct-cor_wtr_f_wrp_6.0.dcd}
set filelist {sim_mscs-cry-chnl-big-oct-cor_wtr_f_salt_wrp_0.5.dcd sim_mscs-cry-chnl-big-oct-cor_wtr_f_salt_wrp_1.0.dcd sim_mscs-cry-chnl-big-oct-cor_wtr_f_salt_wrp_1.5.dcd sim_mscs-cry-chnl-big-oct-cor_wtr_f_salt_wrp_2.0.dcd sim_mscs-cry-chnl-big-oct-cor_wtr_f_salt_wrp_2.5.dcd sim_mscs-cry-chnl-big-oct-cor_wtr_f_salt_wrp_3.0.dcd sim_mscs-cry-chnl-big-oct-cor_wtr_f_salt_wrp_3.5.dcd sim_mscs-cry-chnl-big-oct-cor_wtr_f_salt_wrp_4.0.dcd sim_mscs-cry-chnl-big-oct-cor_wtr_f_salt_wrp_4.5.dcd sim_mscs-cry-chnl-big-oct-cor_wtr_f_salt_wrp_5.0.dcd sim_mscs-cry-chnl-big-oct-cor_wtr_f_salt_wrp_5.5.dcd sim_mscs-cry-chnl-big-oct-cor_wtr_f_salt_wrp_6.0.dcd}
#set filelist {sim_mscs-cry-chnl-big-oct-cor-l109s_wtr_f_salt_wrp_0.5.dcd sim_mscs-cry-chnl-big-oct-cor-l109s_wtr_f_salt_wrp_1.0.dcd sim_mscs-cry-chnl-big-oct-cor-l109s_wtr_f_salt_wrp_1.5.dcd sim_mscs-cry-chnl-big-oct-cor-l109s_wtr_f_salt_wrp_2.0.dcd sim_mscs-cry-chnl-big-oct-cor-l109s_wtr_f_salt_wrp_2.5.dcd sim_mscs-cry-chnl-big-oct-cor-l109s_wtr_f_salt_wrp_3.0.dcd sim_mscs-cry-chnl-big-oct-cor-l109s_wtr_f_salt_wrp_3.5.dcd sim_mscs-cry-chnl-big-oct-cor-l109s_wtr_f_salt_wrp_4.0.dcd sim_mscs-cry-chnl-big-oct-cor-l109s_wtr_f_salt_wrp_4.5.dcd sim_mscs-cry-chnl-big-oct-cor-l109s_wtr_f_salt_wrp_5.0.dcd sim_mscs-cry-chnl-big-oct-cor-l109s_wtr_f_salt_wrp_5.5.dcd sim_mscs-cry-chnl-big-oct-cor-l109s_wtr_f_salt_wrp_6.0.dcd}


#Set selection to look for atoms
set selection [atomselect top "name OH2 and (x^2 + y ^2) < 100 and z > -26 and z < -18"]
set selection_prm [atomselect top "name OH2 and (x^2 + y ^2) < 100 and z > -30 and z < -15"]
#set selection_prm $selection

#Set the parameter to extract
set parameter z


#set filelist {test_wtr.dcd}

foreach crnt_file $filelist {
	#load file
	animate read dcd $crnt_file beg 0 end -1 waitfor all



	#Extract frames from file
	set num_steps [molinfo top get numframes]

	#Open files for writing
	set filename $crnt_file.sel-atm-cnt.txt
	set filename_prm $crnt_file.sel-atm-prm.txt
	set fid [open $filename w]
	set fid_prm [open $filename_prm w]

	for {set frame 0} {$frame < $num_steps} {incr frame} {
		#Update the frame
		$selection frame $frame
		$selection update
		$selection_prm frame $frame
		$selection_prm update
		
		#Count atoms in selection
		set found_atoms [$selection num]
		
		#Extract parameter for atoms
		set extracted_parameter [$selection_prm get $parameter]

		#compose frame number to write to the file
		set frame_write [format "%04d" [expr {$frame + $first_frame}]]
			
			
		#Write results into file
		puts $fid "$frame_write\t$found_atoms"
		puts $fid_prm "$frame_write $extracted_parameter"
		
#		$selection delete
		puts "global frame [expr {$frame + $first_frame}]     file $crnt_file       frame $frame of [expr {($num_steps - 1)}] finished"
	}
	#Close written files
	close $fid
	close $fid_prm

	set first_frame [expr {$first_frame + $num_steps}]
	animate delete all
}


bell
puts "Finished !!!"