# Script to align all the residues in the selection (suggesting all the residues have the same atoms number), sort by RMSD, and save them as a trajectory (one residue - one time frame. It imitates dynamics of the single residue based on the single multiresidue snapshot. For multiple timesteps it repeats the operation for every timestep and appends them all in one long DCD trajectory.
# Andriy Anishkin (anishkin@icqmail.com) UMCP

set originalMolID [molinfo top]; # Original molecule with multiple residues
# set selection [atomselect $originalMolID "resname C16"];	#Residues of this selection will be converted to trajectory
set selection [atomselect $originalMolID "resname POPC"];	#Residues of this selection will be converted to trajectory
set outputFileNameSuffix "_sel-traj";	# Output filename suffix
set outputStructureFileName [lindex [molinfo $originalMolID get name] 0]${outputFileNameSuffix}.pdb
set outputTrajectoryFileName [lindex [molinfo $originalMolID get name] 0]${outputFileNameSuffix}.dcd

set framesNumber [molinfo $originalMolID get numframes]

$selection frame 0
$selection update

set residuesList [lsort -unique [$selection get residue]]
set residuesNumber [llength $residuesList]
foreach residue $residuesList {; # Create separate selection for every residue
	set residueSelection($residue) [atomselect $originalMolID "residue $residue"]
}

set firstResidueSelection $residueSelection([lindex $residuesList 0]); # Find the first residue

#Move the geometry center of the first residue to the origin
$firstResidueSelection moveby [vecscale -1.0 [measure center $firstResidueSelection]]

animate write pdb $outputStructureFileName beg 0 end 0 waitfor all sel $firstResidueSelection $originalMolID; # Save the structure for the output DCD trajectory
mol new $outputStructureFileName type pdb waitfor all; # Load newly created trajectory structure
set outputMolID [molinfo top]; # Original molecule with multiple residues
set lastResidueSelection [atomselect $outputMolID all]; # This will be the last (current) residue in the output trajectory for the comparison. All other residues will be compared with it


for {set frame 0} {$frame < $framesNumber} {incr frame} {;	# Go through the timesteps, compare molecules and write frames
	puts "Processing frame $frame of [expr {$framesNumber-1}]..."
	set currentResiduesList $residuesList;	# Temporary residues list for comparison with the last reidue in output trajectory. Processed residues (for the original frame) will be marked with -1
	foreach residue $residuesList {; # Update separate selection for every residue
		$residueSelection($residue) frame $frame
	}
	
	for {set r 0} {$r < $residuesNumber} {incr r} {; # Go through the residues list
		set minRMSD {}
		set closestResidue {}
		set closestResidueNumber {}
		for {set cr 0} {$cr < $residuesNumber} {incr cr} {; # Compare all the remaining residues and compare with the last one in the output trajectory
			set residue [lindex $currentResiduesList $cr]
			if {$residue!=-1} {; #Unprocessed residue
				$residueSelection($residue) move [measure fit $residueSelection($residue) $lastResidueSelection]; # do the alignment
				set RMSD [measure rmsd $residueSelection($residue) $lastResidueSelection]; # Measure RMSD
				if {($RMSD<$minRMSD)||($minRMSD=={})} {; # This residue is more similar to the last one (or it is the first residues compared)
					set minRMSD $RMSD
					set closestResidue $residue
					set closestResidueNumber $cr
				}
			}; #Unprocessed residue
		}; # Compare all the remaining residues and compare with the last one in the output trajectory
		
		# Write the most similar residue in a new frame in the DCD
		animate dup frame [expr {[molinfo $outputMolID get numframes]-1}] $outputMolID;	# Duplicate the last frame in the DCD
		$lastResidueSelection frame [molinfo $outputMolID get numframes]; # Update the selection with the last frame in the output
		$residueSelection($closestResidue) move [measure fit $residueSelection($residue) $firstResidueSelection]; # do the alignment relative to the first residue
		$lastResidueSelection set {x y z} [$residueSelection($closestResidue) get {x y z}]; # Insert coordinates of the new closest conformation into the last frame of the output trajectory
# 		puts "residue $r of [expr {$residuesNumber-1}].   Closest residue $closestResidueNumber (resid $closestResidue)"
		set currentResiduesList [lreplace $currentResiduesList $closestResidueNumber $closestResidueNumber -1]
	}; # Go through the residues list
};	# Go through the timesteps, compare molecules and write frames

puts "Writing output trajectory ..."
animate write dcd $outputTrajectoryFileName beg 1 end -1 waitfor all sel $lastResidueSelection $outputMolID; # Save the output DCD trajectory, all except the first (duplicated) residue

foreach residue $residuesList {; # Remove all the selections
	$residueSelection($residue) delete 
}; # Remove all the selections

puts "Finished"