mol new dna_base.pdb
mol new 5GZB.pdb

set protAll [atomselect 1 all]
set protProt [atomselect 1 "protein"]
set protDnaAtoms [ atomselect 1 "chain B and resid 4 to 10 and name P"]

set dnaAtoms1 [atomselect 0 "chain A and resid 8 to 14 and name P"]
set dnaAtoms2 [atomselect 0 "chain A and resid 17 to 23 and name P"]

set protSe [atomselect 1 "name SE"]
$protSe set name "S"
set protMse [atomselect 1 "resname MSE"]
$protMse set resname "MET"

#Begin moving the protein. 

set fm [measure fit $protDnaAtoms $dnaAtoms1]
$protAll move $fm
$protProt set chain A
$protProt writepdb prot_1.pdb

set fm [measure fit $protDnaAtoms $dnaAtoms2]
$protAll move $fm
$protProt set chain B
$protProt writepdb prot_2.pdb


set dnaAde [atomselect 0 "resname DA"]
$dnaAde set resname ADE
set dnaCyt [atomselect 0 "resname DC"]
$dnaCyt set resname CYT
set dnaGua [atomselect 0 "resname DG"]
$dnaGua set resname GUA
set dnaThy [atomselect 0 "resname DT"]
$dnaThy set resname THY

set dnaA [atomselect 0 "chain A and noh"]
set dnaB [atomselect 0 "chain B and noh"]
$dnaA set chain C
$dnaB set chain D
$dnaA writepdb dnaC.pdb
$dnaB writepdb dnaD.pdb

set patchesA [lsort -unique [$dnaA get {resid resname}]]
set patchesB [lsort -unique [$dnaB get {resid resname}]]




package require psfgen
set paramDir "/n/projects/cm2363/tead4Simulation/params/charmm27"
topology "$paramDir/top_all27_prot_na.rtf"
#topology "$paramDir/top_all27_na.rtf"
alias residue HIS HSE

segment AP {
    first NTER
    last CTER
    pdb prot_1.pdb
}
coordpdb prot_1.pdb AP

segment BP {
    first NTER
    last CTER
    pdb prot_2.pdb
}
coordpdb prot_2.pdb BP

segment CN {
    first 5TER 
    last 3TER
    pdb dnaC.pdb
}
foreach record $patchesA {
    #This is an ugly way to assign variables. Yeesh!
    foreach {resid resname} $record { break }
    if { $resname == "THY" || $resname == "CYT" } {
        patch DEO1 CN:$resid
    }
    if { $resname == "ADE" || $resname == "GUA" } { 
        patch DEO2 CN:$resid
    }
}

coordpdb dnaC.pdb CN

segment DN {
    first 5TER 
    last 3TER
    pdb dnaD.pdb
}
foreach record $patchesB {
    foreach {resid resname} $record { break }
    if { $resname == "THY" || $resname == "CYT" } {
        patch DEO1 DN:$resid
    }
    if { $resname == "ADE" || $resname == "GUA" } { 
        patch DEO2 DN:$resid
    }
}
coordpdb dnaD.pdb DN


guesscoord
writepdb combined.pdb
writepsf combined.psf


#Now center the molecule and measure the max distance.
mol new combined.psf
mol addfile combined.pdb
set sel [atomselect top all]
$sel moveby [vecinvert [measure center $sel]]
$sel writepdb combined.pdb

set nohSel [atomselect top "noh"]
set nohIndexes [$nohSel get index]
set numAtoms [llength $nohIndexes]
set maxDist 0
set maxi 0
set maxj 0
for {set i 0} {$i < $numAtoms - 1} {incr i} {
    puts -nonewline "\r $i / $numAtoms          "
    for {set j [expr {$i + 1}]} {$j < $numAtoms} {incr j} {
        set idxI [lindex $nohIndexes $i]
        set idxJ [lindex $nohIndexes $j]
        set curDist [measure bond [list $idxI $idxJ] ]
        if {$curDist > $maxDist} { 
            puts "$curDist $idxI $idxJ"
            set maxi $idxI
            set maxj $idxJ
            set maxDist $curDist
        }
    }
}
puts "$maxi $maxj $maxDist"

exit