ColorCoord 1.0
	VMD Version 1.8.3 or greater
	psfgen plugin


	Script for setting the user field of a residue according to the number
	of hbonds it has with its neighbors.
	Initially intended for ice nucleation analysis - it allows to trace down
	the crystallization nucleus in the water phase

	Script uses VMD measure hbonds command to find all hydrogen bonds for a given
	atom (water oxygen). 
	It implements a support for PBC - firstly 27 replicas of the unit simulation
	cell are created, then hydrogen bonds are analyzed (measure hbonds command
	doesn't support hbonds going over the boundary - this is supposed to change in
	the future).
	The resulting data can be stored/retrieved from a file, averaged, and user
	field can be set according to the obtained values.

	GetMolHbonds distance angle first last
  		calculates the 'coordination numbers' and return an "array" with the data		
	WriteMolHbondsFile n_hbonds_array filename
		writes the data to file
	ReadMolHbondsFile filename
		retrieves the data from file
	AverageMolHbonds n_hbonds_array
		averages the data
	SetMolHbondsUser n_hbonds_array
		sets the user field according to data

	The script defines some global variables that need to be changed to fit your
	purpose. The file is heavily commented so you shouldn't have any problems, however
	please contact me if you encounter any.
	firstly, set all CC_water variables according to your simulation. in my case i have
	6-center water model with OW, HW1, HW2, EP, LP1, LP2; i don't care about the
	extra charges...

	# names for atom selections
	set CC_water_resname "NE6"
	set CC_water_oxygen_name "OW"
	# you can set more, e.g., the lone pairs, ...
	set CC_water_all "OW \"HW.\""
	# total number of particles per residue according to CC_water_all selection
	set CC_water_all_nparticles 3

	# you need to create psf file for your system and namd topology
	# then set the following variables to fit your needs
	# 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 "_ne6_namd.psf"
	set CC_topname ""
	# you can keep this as is - consult the code if you need to know more
	set CC_pbc_threshold 4.0	

	# create the array containing number of hbonds per water molecule in the unit cell
	# this can take really long time!
	# parameters are distance/angle criterion for hbonds, first and last frame (starting
	# from 1)
	array set n_hbonds_array [GetMolBonds 3.0 20 1 10]

	# save to a file if you need
	WriteMolHbondsFile [array get n_hbonds_array] n_hbonds_filename

	# retrieve data if you need
	# it will set all the global variables to values that were used at the time
	# when the data were saved
	array set loaded_n_bonds_array [ReadMolHbondsFile n_hbonds_filename]

	# perform a kind of averaging (pseudo-running average) on the arrays with hbonds
	# if you need - in this case 5-sample averages are done
	# check the code of the proc if you want to know more
	array set n_hbonds_array_averaged [AverageMolHbonds [array get n_hbonds_array] 5]

	# finally set the user field according to the number of hbonds
	SetMolHbondsUser [array get n_hbonds_array]
	# or 
	SetMolHbondsUser [array get n_hbonds_array_averaged]

	# then set a selection 
	# user > 3.5
	# in the selection field for graphical representation and see... :)


	Lubos Vrbka (lubos (dot) vrbka (at) gmail (dot) com)