#
# THE MAIN BODY
#

## The molecules: 0 compared to 1 (possibly, different trajectories)

mol new Nonwater.psf type psf waitfor all
mol new Nonwater.psf type psf waitfor all

set dcd_file(0) production1-7_nonwater_Voronoi.dcd
set numframes(0) 55983

set dcd_file(1) production1-7_nonwater_Voronoi.dcd
set numframes(1) 55983

set same_flag 1
set frame_step 10

## RMSD data file to write:

set rmsd_file_P [open "RMSD_map_P.dat" "w"]
set rmsd_file_noterm [open "RMSD_map_noterm.dat" "w"]
set rmsd_file_core_bases [open "RMSD_map_core_bases.dat" "w"]
set rmsd_file_core_bckb [open "RMSD_map_core_bckbdat" "w"]
set rmsd_file_all [open "RMSD_map_all.dat" "w"]

## Atom selections to be used:

foreach i { 0 1 } {
# set all_atoms($i)     [atomselect $i "all"]
set dna_all($i)      [atomselect $i "noh and not name SOD"]
set dna_noterm($i)     [atomselect $i "noh and not name SOD and not resid 1 25"]
set dna_P($i)        [atomselect $i "name P"]
set dna_core_bases($i) [atomselect $i "((segid ADNA and resid 11 to 16) or (segid BDNA and resid 10 to 15)) and noh and not name C2' C3' O3' C4' O4' C5' O5' P O1P O2P"]
set dna_core_bckb($i)  [atomselect $i "((segid ADNA and resid 11 to 16) or (segid BDNA and resid 10 to 15)) and name C1' C2' C3' O3' C4' O4' C5' O5' P"]
}

#
# The main loop.

for { set i 0 } { $i < $numframes(0) } { incr i $frame_step } {

    animate delete all 0
    mol addfile $dcd_file(0) first $i last $i waitfor all 0

    if { $same_flag } then {
      set jfirst [expr $i+$frame_step]
    } else {
      set jfirst 0
    }
for { set j $jfirst } { $j < $numframes(1) } { incr j $frame_step } {

    puts "======= COMPARISON FRAMES $i-$j  =========="

    animate delete all 1
    mol addfile $dcd_file(1) first $j last $j waitfor all 1

# Phosphates alignment
	set transformation_matrix [measure fit $dna_P(1) $dna_P(0)]
	$dna_all(1) move $transformation_matrix
	set rmsd [measure rmsd $dna_P(1) $dna_P(0)]
	puts $rmsd_file_P [format "%5.2f" $rmsd]

# Core 1 (no termini) atoms alignment
	set transformation_matrix [measure fit $dna_noterm(1) $dna_noterm(0)]
	$dna_all(1) move $transformation_matrix
	set rmsd [measure rmsd $dna_noterm(1) $dna_noterm(0)]
	puts $rmsd_file_noterm [format "%5.2f" $rmsd]

# Core 2 atoms alignment (backbone)
        set transformation_matrix [measure fit $dna_core_bckb(1) $dna_core_bckb(0)]
        $dna_all(1) move $transformation_matrix
        set rmsd [measure rmsd $dna_core_bckb(1) $dna_core_bckb(0)]
        puts $rmsd_file_core_bckb [format "%5.2f" $rmsd]

# Core 2 atoms alignment (bases)
        set transformation_matrix [measure fit $dna_core_bases(1) $dna_core_bases(0)]
        $dna_all(1) move $transformation_matrix
        set rmsd [measure rmsd $dna_core_bases(1) $dna_core_bases(0)]
        puts $rmsd_file_core_bases [format "%5.2f" $rmsd]

# All atoms alignment
	set transformation_matrix [measure fit $dna_all(1) $dna_all(0)]
	$dna_all(1) move $transformation_matrix
	set rmsd [measure rmsd $dna_all(1) $dna_all(0)]
	puts $rmsd_file_all [format "%5.2f" $rmsd]
    }
}

close $rmsd_file_P
close $rmsd_file_noterm
close $rmsd_file_core_bckb
close $rmsd_file_core_bases
close $rmsd_file_all