# Created by Eduard Schreiner (eschrein@ks.uiuc.edu)

proc ww_nearest { macromol what cutoff howmany newtrajname {mol top} } {
    set newmol [mol new]
    set macromol_sel [atomselect $mol "${macromol}"]
    set macromol_inds [$macromol_sel get index]
    set nr_steps [molinfo $mol get numframes ]
    
    for {set t 0} {$t < $nr_steps} {incr t } {
	$macromol_sel frame $t; $macromol_sel update
	#determime what has to be written
	set a [atomselect $mol "noh and (same residue as ((${what}) and within $cutoff of ${macromol}))" frame $t]
	set inds [$a get index]
	$a delete
	#check if we got enough "water" molecules within the specified cutoff
	if { [llength $inds] >= $howmany } {
	    #this list of lists keeps the "water" index and its closest distance to the "protein"
	    set all_dists [list]
	    foreach i $inds {
		set a [atomselect $mol "same residue as index $i" frame $t]
		set cont_w_m [measure contacts $cutoff $a $macromol_sel]
		$a delete
		#determine the closest contact
		#use a temporary list to avoid double counts later
		set tmp_dists [list]
		foreach w [lindex $cont_w_m 0] p [lindex $cont_w_m 1] {
		    set tmp_var [measure bond [list $w $p] molid $mol frame $t]
		    lappend tmp_dists [list $tmp_var $i]
		}
		#get the smallest distance to the protein
		lappend all_dists [lindex [lsort -increasing $tmp_dists] 0]
	    }
	    #sort the resulting distance/index list according to the distance
	    #i.e. the "water" molecules are sorted according to their distance from the "protein"
	    set all_dists_sorted [lsort -increasing $all_dists]
	    #get the first $howmany elements
	    set tmp_inds [list]
	    for {set i 0} {$i < $howmany} {incr i} {
		lappend tmp_inds [lindex $all_dists_sorted $i]
	    }
	    set water_inds [join $tmp_inds " "]
	    set writesel [atomselect $mol "($macromol) or (same residue as index $water_inds)" frame $t]
	    $writesel writepdb tmp-molwat.pdb
	    if {$t == 0} {
		$writesel writepsf ${newtrajname}.psf
	    }
	    $writesel delete
	    mol addfile tmp-molwat.pdb $newmol
	} else { 
	    puts "there are less than $howmany water molecules within the specified cutoff of $cutoff"
	}
    }
    $macromol_sel delete
    
    #write the trimmed trajectory
    animate write dcd  ${newtrajname}.dcd beg 0 end [expr $nr_steps -1] waitfor all $newmol
    mol delete $newmol
    exec rm -f tmp-molwat.pdb
}