Difference between revisions of "User:Remig/plico/plicoCommonNT"

From Jmol
Jump to navigation Jump to search
(Handle multiple frames)
(Avoid "axis," a newly reserved word)
 
(5 intermediate revisions by the same user not shown)
Line 3: Line 3:
 
Copy and paste the following into a text editor and save in your scripts folder as plicoNTcommon.spt.
 
Copy and paste the following into a text editor and save in your scripts folder as plicoNTcommon.spt.
 
<pre>#  plicoNTcommon - Jmol script by Ron Mignery
 
<pre>#  plicoNTcommon - Jmol script by Ron Mignery
#  v1.4 beta    7/8/2014 -handle multiple frames
+
#  v1.10 beta    4/12/2016 -axis is now a reserved word
 
#
 
#
 
#  Routines and values common to Plico suite scripts that work with nucleotides
 
#  Routines and values common to Plico suite scripts that work with nucleotides
#  Must be present in the same directory as other Plico scripts that use it
+
#  Must be present in the same directory as other Plico scripts that use them
kNTcommon = 4
+
kNTcommon = 6
 
kC5O5PO3B = -71.0
 
kC5O5PO3B = -71.0
 
kO5PO3C3B = -107.0
 
kO5PO3C3B = -107.0
 
kPO3C3C4B = -161.5
 
kPO3C3C4B = -161.5
kO3C3C4C5B = 125.5
+
kO3C3C4C5B = 140.0
 
kC3C4C5O5B = 55.65
 
kC3C4C5O5B = 55.65
 
kC4C5O5PB = 169.0
 
kC4C5O5PB = 169.0
  
 
kO4C4C3C2B = 15.92
 
kO4C4C3C2B = 15.92
kC4O4C1C2B = -41.7
+
kC4O4C1C2B = -19.9 #-41.7 1bna minimized
kC4O4C1NxB = -159.0
+
kC2O4C1NxB = -122.6 #-159.0 1bna minimized
kC5C4O4C1B = 146.3
+
kC5C4O4C1B = 122.2 #146.3 1bna minimized
 
kC3C1C2O2B = 120.5
 
kC3C1C2O2B = 120.5
  
Line 33: Line 33:
 
kO4C4C3C2A = -35.55
 
kO4C4C3C2A = -35.55
 
kC4O4C1C2A = 3.8
 
kC4O4C1C2A = 3.8
kC4O4C1NxA = -117.4
+
kC2O4C1NxA = -131.0
 
kC5C4O4C1A = 144.85
 
kC5C4O4C1A = 144.85
 
kC3C1C2O2A = 116.3
 
kC3C1C2O2A = 116.3
Line 46: Line 46:
 
function force_p_res(cres, iChain) {
 
function force_p_res(cres, iChain) {
 
     var pres = cres-1
 
     var pres = cres-1
     var aP =  atom_rcn( cres, iChain, "P")
+
     var aP =  get_atom_rcn( cres, iChain, "P")
     var aO5 =  atom_rcn( cres, iChain, "O5\'")
+
     var aO5 =  get_atom_rcn( cres, iChain, "O5\'")
     var aC5 =  atom_rcn( cres, iChain, "C5\'")
+
     var aC5 =  get_atom_rcn( cres, iChain, "C5\'")
     var aC4 =  atom_rcn( cres, iChain, "C4\'")
+
     var aC4 =  get_atom_rcn( cres, iChain, "C4\'")
     var aOP1 = atom_rcn( cres, iChain, "OP1")
+
     var aOP1 = get_atom_rcn( cres, iChain, "OP1")
     var aOP2 = atom_rcn( cres, iChain, "OP2")
+
     var aOP2 = get_atom_rcn( cres, iChain, "OP2")
     var aO3p = atom_rcn( pres, iChain, "O3\'")
+
     var aO3p = get_atom_rcn( pres, iChain, "O3\'")
     var aC3p = atom_rcn( pres, iChain, "C3\'")
+
     var aC3p = get_atom_rcn( pres, iChain, "C3\'")
     if (aO3p.size > 0) {
+
     if (aO3p) {
 
 
        set_distance_atoms(aC5, aO3p, 3.1)
 
  
 +
        var selsave = {selected}
 +
        set_distance_atoms(aP3p, aC5, 3.1)
 
         select aO5
 
         select aO5
 
         var dist = distance(aO3p, aO5)
 
         var dist = distance(aO3p, aO5)
 
         var widen = (dist < 2.85)
 
         var widen = (dist < 2.85)
 
         var dir = (widen ? -1 : 1)
 
         var dir = (widen ? -1 : 1)
         var first = TRUE
+
         var first = true
 
         while (abs(dist-2.85) > kDtolerance) {
 
         while (abs(dist-2.85) > kDtolerance) {
 
             rotateSelected @aC4 @aC5 @dir
 
             rotateSelected @aC4 @aC5 @dir
Line 76: Line 76:
 
             }
 
             }
 
             dist=newdist
 
             dist=newdist
             first = FALSE
+
             first = false
 
         }
 
         }
 
         select aP
 
         select aP
Line 85: Line 85:
 
         aOP2.xyz = get_tet_idx(aO3p.atomIndex, aP.atomIndex, aO5.atomIndex, 1.73)
 
         aOP2.xyz = get_tet_idx(aO3p.atomIndex, aP.atomIndex, aO5.atomIndex, 1.73)
 
         aOP1.xyz = get_tet_idx(aO5.atomIndex, aP.atomIndex, aO3p.atomIndex, 1.73)
 
         aOP1.xyz = get_tet_idx(aO5.atomIndex, aP.atomIndex, aO3p.atomIndex, 1.73)
         minimize select {connected(aP) or aP}
+
         plico_minimize( {(connected(aP) or aP) and not aO3p})
 +
        select selsave
 
     }
 
     }
 
}
 
}
  
 
function fix_p_res(cres, iChain, force) {
 
function fix_p_res(cres, iChain, force) {
    var f = (_frameID/1000000)
+
print format("fix_p_res(cres=%d, ichain=%s, force=%s)", cres, iChain, force)
    var m = (_frameID%1000000)
 
 
     var pres = cres-1
 
     var pres = cres-1
     var aP  = atom_rcn( cres, iChain, "P")
+
     var aP  = get_atom_rcn( cres, iChain, "P")
     var aO5  = atom_rcn( cres, iChain, "O5\'")
+
     var aO5  = get_atom_rcn( cres, iChain, "O5\'")
     var aC5  = atom_rcn( cres, iChain, "C5\'")
+
     var aC5  = get_atom_rcn( cres, iChain, "C5\'")
     var aC4  = atom_rcn( cres, iChain, "C4\'")
+
     var aC4  = get_atom_rcn( cres, iChain, "C4\'")
     var aOP1 = atom_rcn( cres, iChain, "OP1")
+
    var aC1  = get_atom_rcn( cres, iChain, "C1\'")
     var aOP2 = atom_rcn( cres, iChain, "OP2")
+
     var aOP1 = get_atom_rcn( cres, iChain, "OP1")
     var aO3p = atom_rcn( pres, iChain, "O3\'")
+
     var aOP2 = get_atom_rcn( cres, iChain, "OP2")
     var aC3p = atom_rcn( pres, iChain, "C3\'")
+
     var aO3p = get_atom_rcn( pres, iChain, "O3\'")
     if ((aO3p.size > 0)and (aC4.size > 0)) {
+
     var aC3p = get_atom_rcn( pres, iChain, "C3\'")
 +
    var aC4p  = get_atom_rcn( pres, iChain, "C4\'")
 +
     if (aO3p.size and aC4.size) {
 +
        var selsave = {selected}
  
 
         # If collision
 
         # If collision
Line 107: Line 110:
 
             # Push away
 
             # Push away
 
             select {(resno <= @{aC5.resno}) and (chain=iChain)
 
             select {(resno <= @{aC5.resno}) and (chain=iChain)
                 and (file=f) and (model=m)}
+
                 and thisModel}
 
             set_distance_atoms(aC3p, aC5, kCtolerance)
 
             set_distance_atoms(aC3p, aC5, kCtolerance)
 
         }
 
         }
  
         select aO5
+
        # Rotate C4'-C5' until P-O3' is 1.59
         var dist = distance(aO3p, aO5)
+
         select aP
         var widen = (dist < 2.85)
+
        set_distance_atoms(aO5, aP, 1.59)
         var dir = (widen ? -1 : 1)
+
         set_angle_atoms(aC5, aO5, aP, 109)
         var first = TRUE
+
         set_dihedral_atoms(aC4, aC5, aO5, aP, 180)
        while (abs(dist-2.85) > kDtolerance) {
+
        select add aO5
             rotateSelected @aC4 @aC5 @dir
+
         var dist = distance(aP, aO3p)
             var newdist = distance(aO3p, aO5)
+
         if (dist > 1.59) {
            if (widen ? (newdist < dist) : (newdist > dist)) {
+
             var dir = 1.0
                if (first) {
+
             for (var i = 0; i < 180; i++) {
                    dir = -dir
+
                dist = distance(aP, aO3p)
                     rotateSelected @aC4 @aC5 @dir
+
                if ((dist-1.59) < 0.1) {
 +
                     break
 
                 }
 
                 }
                 else {
+
                 rotateSelected @aC4 @aC5 @dir
                     break
+
                var newdist = distance(aP, aO3p)
 +
                if (newDist > dist) {
 +
                    rotateSelected @aC4 @aC5 @{-dir}
 +
                    if (dir > 0) {
 +
                        dir = -dir
 +
                    }
 +
                    else {
 +
                        break
 +
                     }
 
                 }
 
                 }
             }
+
             } #endfor 180
            dist=newdist
 
            first = FALSE
 
 
         }
 
         }
         if (force and (abs(distance(aO3p, aC5)-4.11) < kDtolerance)) {
+
         else {
             # Push away
+
             var dir = -1.0
             select {(resno <= @{aC5.resno}) and (chain=iChain)
+
             for (var i = 0; i < 180; i++) {
                 and (file=f) and (model=m)}
+
                 dist = distance(aP, aO3p)
            set_distance_atoms(aO3p, aC5, 4.11)
+
                if ((1.59-dist) < 0.1) {
        }
+
                     break
        select aP
 
        set_distance_atoms(aO5, aP, 1.73)
 
        set_angle_atoms(aC5, aO5, aP, 110.1)
 
 
 
        dist = distance(aO3p, aP)
 
        widen = (dist < 1.73)
 
        dir = (widen ? -1 : 1)
 
        first = TRUE
 
        while (abs(dist-1.73) > kDtolerance) {
 
            rotateSelected @aC5 @aO5 @dir
 
            var newdist = distance(aO3p, aP)
 
            if (widen ? (newdist < dist) : (newdist > dist)) {
 
                if (first) {
 
                     dir = -dir
 
                    rotateSelected @aC5 @aO5 @dir
 
 
                 }
 
                 }
                 else {
+
                 rotateSelected @aC4 @aC5 @dir
                     break
+
                var newdist = distance(aP, aO3p)
 +
                if (newDist < dist) {
 +
                    rotateSelected @aC4 @aC5 @{-dir}
 +
                    if (dir < 0) {
 +
                        dir = -dir
 +
                    }
 +
                    else {
 +
                        break
 +
                     }
 
                 }
 
                 }
             }
+
             } #endfor 180
            dist=newdist
 
            first = FALSE
 
 
         }
 
         }
 
        #set_dihedral_atoms(aC4, aC5, aO5, aP, 150.3)
 
 
         aOP2.xyz = get_tet_idx(aO3p.atomIndex, aP.atomIndex, aO5.atomIndex, 1.73)
 
         aOP2.xyz = get_tet_idx(aO3p.atomIndex, aP.atomIndex, aO5.atomIndex, 1.73)
 
         aOP1.xyz = get_tet_idx(aO5.atomIndex, aP.atomIndex, aO3p.atomIndex, 1.73)
 
         aOP1.xyz = get_tet_idx(aO5.atomIndex, aP.atomIndex, aO3p.atomIndex, 1.73)
 
         if (force) {
 
         if (force) {
             minimize select {connected(aP) or aP}
+
             select {aP or aOP1 or aOP2}
 +
            plico_minimize( {selected})
 
         }
 
         }
 +
        select selsave
 
     }
 
     }
 
}
 
}
  
 +
function fix_p_res_range(res5, res3, iChain, force) {
 +
    for (var i = res5; i <= res3; i++) {
 +
        fix_p_res(i, iChain, force)
 +
    }
 +
}
 +
 +
function fix_all_nt_collisions( iChain) {
 +
    var selsave = {selected}
 +
    chset = count_collisions(true)
 +
    for (var i = 1; i <= chset.size; i++) {
 +
        c = chset[i]
 +
        cset = (within(kCtolerance, c) and  not c
 +
            and not connected(chset[i]))
 +
        rset = [{resno=@{c.resno}}]
 +
        for (var j = 1; j <= cset.size; j++) {
 +
            rset += {resno=@{cset[j].resno}}
 +
        }
 +
       
 +
        if ({c and cset and not base}) {
 +
            select {@{rset} and base}
 +
            plico_minimize( {selected})
 +
        }
 +
        else if (c.atomname[1][2] == "OP") {
 +
            fix_p_res(chset[i].resno, iChain, true)
 +
        }
 +
        else {
 +
            plico_minimize(rset)
 +
        }
 +
    }
 +
    measure off
 +
    select selsave
 +
}
 +
 +
# The following functions position one nt relative to another:
 +
# Common positioning functions:
 
function get_rotors_res(res) {
 
function get_rotors_res(res) {
    var f = (_frameID/1000000)
 
    var m = (_frameID%1000000)
 
 
     var rotors = array()
 
     var rotors = array()
 
     var sRes = res
 
     var sRes = res
 
     var mRes = sRes-1
 
     var mRes = sRes-1
 
     var iChain = {(resno=res) and (atomName="P")
 
     var iChain = {(resno=res) and (atomName="P")
         and (file=f) and (model=m)}.chain
+
         and thisModel}.chain
     var mC4 = atom_rcn( mRes, iChain, "C4\'")
+
     var mC4 = get_atom_rcn( mRes, iChain, "C4\'")
     var mC3 = atom_rcn( mRes, iChain, "C3\'")
+
     var mC3 = get_atom_rcn( mRes, iChain, "C3\'")
     var mO3 = atom_rcn( mRes, iChain, "O3\'")
+
     var mO3 = get_atom_rcn( mRes, iChain, "O3\'")
     var sP =  atom_rcn( sRes, iChain, "P"  )
+
     var sP =  get_atom_rcn( sRes, iChain, "P"  )
     var sO5 = atom_rcn( sRes, iChain, "O5\'")
+
     var sO5 = get_atom_rcn( sRes, iChain, "O5\'")
     var sC5 = atom_rcn( sRes, iChain, "C5\'")
+
     var sC5 = get_atom_rcn( sRes, iChain, "C5\'")
     var sC4 = atom_rcn( sRes, iChain, "C4\'")
+
     var sC4 = get_atom_rcn( sRes, iChain, "C4\'")
     var sC3 = atom_rcn( sRes, iChain, "C3\'")
+
     var sC3 = get_atom_rcn( sRes, iChain, "C3\'")
  
 
     rotors += [mC4.atomIndex, mC3.atomIndex, mO3.atomIndex, sP.atomIndex]
 
     rotors += [mC4.atomIndex, mC3.atomIndex, mO3.atomIndex, sP.atomIndex]
Line 197: Line 231:
 
function get_nt_chi_rotor_res(res, iChain) {
 
function get_nt_chi_rotor_res(res, iChain) {
 
     var rotors = array()
 
     var rotors = array()
     var aO4 = atom_rcn( res, iChain, "O4\'")
+
     var aO4 = get_atom_rcn( res, iChain, "O4\'")
     var aC1 = atom_rcn( res, iChain, "C1\'")
+
     var aC1 = get_atom_rcn( res, iChain, "C1\'")
     var isR  = ((aC1 and {purine}).size > 0)
+
     var isR  = (aC1 and {purine})
 
     var N1or9 = (isR ? "N9" : "N1")
 
     var N1or9 = (isR ? "N9" : "N1")
 
     var C6or8 = (isR ? "C8" : "C6")
 
     var C6or8 = (isR ? "C8" : "C6")
  
     var aN = atom_rcn(res, iChain, N1or9)
+
     var aN = get_atom_rcn(res, iChain, N1or9)
     var aC = atom_rcn(res, iChain, C6or8)
+
     var aC = get_atom_rcn(res, iChain, C6or8)
  
 
     rotors = [aO4.atomIndex, aC1.atomIndex, aN.atomIndex, aC.atomIndex]
 
     rotors = [aO4.atomIndex, aC1.atomIndex, aN.atomIndex, aC.atomIndex]
Line 212: Line 246:
 
function get_nt_ab_rotor_res(res, iChain) {
 
function get_nt_ab_rotor_res(res, iChain) {
 
     var rotors = array()
 
     var rotors = array()
     var aC5 = atom_rcn(res, iChain, "C5\'")
+
     var aC5 = get_atom_rcn(res, iChain, "C5\'")
     var aC4 = atom_rcn(res, iChain, "C4\'")
+
     var aC4 = get_atom_rcn(res, iChain, "C4\'")
     var aC3 = atom_rcn(res, iChain, "C3\'")
+
     var aC3 = get_atom_rcn(res, iChain, "C3\'")
     var aO3 = atom_rcn(res, iChain, "O3\'")
+
     var aO3 = get_atom_rcn(res, iChain, "O3\'")
  
 
     rotors = [aO3.atomIndex, aC3.atomIndex, aC4.atomIndex, aC5.atomIndex]
 
     rotors = [aO3.atomIndex, aC3.atomIndex, aC4.atomIndex, aC5.atomIndex]
Line 229: Line 263:
 
}
 
}
  
 +
# ri=moved rj=fixed
 +
function position_nt_by_vs(ri, rj, ares, vs, iChain, jChain) {
 +
    if (ri > rj) {
 +
        var as = gen_as(ri, rj, iChain, jChain)
 +
        select {(resno < ares) and (resno >= ri)
 +
            and (chain=iChain) and thisModel}
 +
    }
 +
    else {
 +
        var as = gen_as(rj, ri, jChain, iChain)
 +
        select {(resno > ares) and (resno <= ri)
 +
            and (chain=iChain) and thisModel}
 +
    }
 +
    move_it(as, vs)
 +
}
 +
               
 
# Moved object must be selected, fixed object not
 
# Moved object must be selected, fixed object not
# as[6] = fixed[1-3] moved[4-6]
+
# as[6] = fixed[1-3] moved[4-6] chis [7-8]
 
# vs[6] = [distance(as[3-4]), angle(as[2-4]),
 
# vs[6] = [distance(as[3-4]), angle(as[2-4]),
 
#  dihedral(as[1-4]), angle(as[5-3], dihedral(as[6-3],
 
#  dihedral(as[1-4]), angle(as[5-3], dihedral(as[6-3],
Line 247: Line 296:
 
     # Dihedral sets TBD
 
     # Dihedral sets TBD
 
     set_dihedral_atoms(as[2], as[3], as[4], as[5], vs[6])
 
     set_dihedral_atoms(as[2], as[3], as[4], as[5], vs[6])
 +
   
 +
    # If chis
 +
    if (vs.size > 6) {
 +
        var sel = {selected}
 +
        select {(resno=@{as[4].resno}) and base}
 +
        set_dihedral_atoms(as[6], as[5], as[4], as[8], vs[8])
 +
    }
 
}
 
}
  
function gen_as(res5, res3, iChain, jChain) {
+
# ri=moved rj=fixed
 +
function gen_as(ri, rj, iChain, jChain) {
 
     var as = array()
 
     var as = array()
     as[1] = atom_rcn(res3, jChain, "C4\'")
+
     as[1] = get_atom_rcn(rj, jChain, "C4\'")
     as[2] = atom_rcn(res3, jChain, "C1\'")
+
     as[2] = get_atom_rcn(rj, jChain, "C1\'")
 
     as[3] = connected(as[2]) and {element="N"}
 
     as[3] = connected(as[2]) and {element="N"}
     as[5] = atom_rcn(res5, iChain, "C1\'")
+
     as[5] = get_atom_rcn(ri, iChain, "C1\'")
     as[6] = atom_rcn(res5, iChain, "C4\'")
+
     as[6] = get_atom_rcn(ri, iChain, "C4\'")
 
     as[4] = connected(as[5]) and {element="N"}
 
     as[4] = connected(as[5]) and {element="N"}
 
+
   
 +
    as[7] = get_atom_rcn(rj, jChain,
 +
        ((as[1] and {purine}) ? "C8" : "C6"))
 +
    as[8] = get_atom_rcn(ri, iChain,
 +
        ((as[6] and {purine}) ? "C8" : "C6"))
 
     return as
 
     return as
 
}
 
}
  
# Pair res5 on res3 moving res <= res3
+
# Specific positioning functions:
 +
# Pair res5 on res3 moving res5 <= res3
 
function pair_it_res(res5, res3, ares, iChain, jChain) {
 
function pair_it_res(res5, res3, ares, iChain, jChain) {
    var f = (_frameID/1000000)
 
    var m = (_frameID%1000000)
 
 
     var as = gen_as(res5, res3, iChain, jChain)
 
     var as = gen_as(res5, res3, iChain, jChain)
 
     var isA = is_form_a(res5, iChain)
 
     var isA = is_form_a(res5, iChain)
 
     var vs = array()
 
     var vs = array()
     vs[1] = 9.00    # distance iN-jN (1tna[11-24]=9.00)
+
     vs[1] = 8.83                    # distance res5 N9or1 and res3 N9or1     
     vs[2] = 124.# angle iN-jN-jC1 (1tna[11-24]=99.6)//126.0
+
     vs[2] = 126.12                  # angle res5 N9or1 and res3 N9or1 C1     
     vs[3] = (isA ? 160.0 : -140.0)   # dihedral iN-jN-jC1-jC4 (1tna[11-24]=160.0)//-142.6
+
     vs[3] = (isA ? 160.0 : -134.97) # dihedral res5 N9or1 and res3 N9or1 C1 C4
     vs[4] = 124.# angle iC1-iN-jN (1tna[11-24]=124.6)
+
     vs[4] = 125.32                  # angle res5 C1 N9or1 and res3 N9or1     
     vs[5] = (isA ? 160.0 : -140.0)   # dihedral iC4-iN-iC1-jN (1tna[11-24]=160.0)//-138.6
+
     vs[5] = (isA ? 160.0 : -141.46) # dihedral res5 C4 N9or1 C1 and res3 N9or1
     vs[6] = (isA ? -5.0 : -1.6)   # dihedral iN-iC1-jN-jC1 (1tna[11-24]=-5.0)//-1.6
+
    vs[6] = (isA ? -5.0 : -17.87)  # dihedral res5 N9or1 C1 and res3 N9or1 C1
 
+
     vs[7] = (isA ? -20.0 : 38)     # dihedral chi res3
 +
    vs[8] = (isA ? -20.0 : 38)     # dihedral chi res5
 
     if (ares < 0) {
 
     if (ares < 0) {
         select ((resno=res5) and (chain=iChain) and (file=f) and (model=m))
+
         select ((resno=res5) and (chain=iChain) and thisModel)
 
     }
 
     }
     else {
+
     else if (ares > 0) {
         select ((resno <= ares) and (chain=iChain) and (file=f) and (model=m))
+
         select ((resno <= ares) and (chain=iChain) and thisModel)
 
     }
 
     }
 
     move_it(as, vs)
 
     move_it(as, vs)
     fix_p_res(res5, iChain, TRUE)
+
     fix_p_res(res5, iChain, true)
     fix_p_res(res3, jChain, TRUE)
+
     fix_p_res(res5+1, iChain, true)
}
 
 
 
# Unstack res5 on res3 moving just res5
 
function single_unstack_res5_on_res3(res5, res3, iChain, jChain) {
 
    var f = (_frameID/1000000)
 
    var m = (_frameID%1000000)
 
    var as = gen_as(res5, res3, iChain, jChain)
 
 
 
    var vs = array()
 
    vs[1] = 4.31    # distance iN-jN (1tna[38-39]=4.31)
 
    vs[2] = 99.6    # angle iN-jN-jC1 (1tna[38-39]=99.6)
 
    vs[3] = -61.1  # dihedral iN-jN-jC1-jC4 (1tna[38-39]=-61.1)
 
    vs[4] = 128.8  # angle iC1-iN-jN (1tna[38-39]=128.8)
 
    vs[5] = 58.4    # dihedral iC4-iN-iC1-jN (1tna[38-39]=58.4)
 
    vs[6] = 78.4    # dihedral iN-iC1-jN-jC1 (1tna[38-39]=78.4)
 
 
 
    select {(resno=res5)  and (chain=iChain) and (file=f) and (model=m)}
 
    move_it(as, vs)
 
    force_p_res(res3, jChain)
 
    move_it(as, vs)
 
    #fix_p_res(res5, iChain, FALSE)
 
    force_p_res(res3, jChain)
 
 
}
 
}
  
 
# Flatstack res5 on res3 moving just res5
 
# Flatstack res5 on res3 moving just res5
 
function single_flatstack_res5_on_res3(res5, res3, iChain, jChain) {
 
function single_flatstack_res5_on_res3(res5, res3, iChain, jChain) {
    var f = (_frameID/1000000)
 
    var m = (_frameID%1000000)
 
 
     var as = gen_as(res5, res3, iChain, jChain)
 
     var as = gen_as(res5, res3, iChain, jChain)
  
 
     var vs = array()
 
     var vs = array()
     vs[1] = 7.00    # distance iN-jN
+
     vs[1] = 7.00    # distance res5 N9or1 and res3 N9or1     
     vs[2] = 89.1    # angle iN-jN-jC1
+
     vs[2] = 89.1    # angle res5 N9or1 and res3 N9or1 C1     
     vs[3] = -49.9  # dihedral iN-jN-jC1-jC4
+
     vs[3] = -49.9  # dihedral res5 N9or1 and res3 N9or1 C1 C4
     vs[4] = 83.4    # angle iC1-iN-jN
+
     vs[4] = 83.4    # angle res5 C1 N9or1 and res3 N9or1     
     vs[5] = 125.7  # dihedral iC4-iN-iC1-jN
+
     vs[5] = 125.7  # dihedral res5 C4 N9or1 C1 and res3 N9or1
     vs[6] = 5.8    # dihedral iN-iC1-jN-jC1
+
     vs[6] = 5.8    # dihedral res5 N9or1 C1 and res3 N9or1 C1
  
     select {(resno=res5) and (chain=iChain) and (file=f) and (model=m)}
+
     select {(resno=res5) and (chain=iChain) and thisModel}
 
     move_it(as, vs)
 
     move_it(as, vs)
     force_p_res(res3, jChain)
+
     #force_p_res(res3, jChain)
     move_it(as, vs)
+
     #move_it(as, vs)
     #fix_p_res(res3, jChain, TRUE)
+
     #fix_p_res(res3, jChain, true)
     force_p_res(res3, jChain)
+
     #force_p_res(res3, jChain)
 
}
 
}
  
 
# Outstack res5 on res3 moving just res5
 
# Outstack res5 on res3 moving just res5
 
function single_outstack_res5_on_res3(res5, res3, iChain, jChain) {
 
function single_outstack_res5_on_res3(res5, res3, iChain, jChain) {
    var f = (_frameID/1000000)
 
    var m = (_frameID%1000000)
 
 
     var as = gen_as(res5, res3, iChain, jChain)
 
     var as = gen_as(res5, res3, iChain, jChain)
  
 
     var vs = array()
 
     var vs = array()
     vs[1] = 8.23  # distance iN-jN
+
     vs[1] = 8.23  # distance res5 N9or1 and res3 N9or1     
     vs[2] = 32.4  # angle iN-jN-jC1
+
     vs[2] = 32.4  # angle res5 N9or1 and res3 N9or1 C1     
     vs[3] = -26.8  # dihedral iN-jN-jC1-jC4
+
     vs[3] = -26.8  # dihedral res5 N9or1 and res3 N9or1 C1 C4
     vs[4] = 99.6  # angle iC1-iN-jN
+
     vs[4] = 99.6  # angle res5 C1 N9or1 and res3 N9or1     
     vs[5] = 57.4  # dihedral iC4-iN-iC1-jN
+
     vs[5] = 57.4  # dihedral res5 C4 N9or1 C1 and res3 N9or1
     vs[6] = 179.1  # dihedral iN-iC1-jN-jC1
+
     vs[6] = 179.1  # dihedral res5 N9or1 C1 and res3 N9or1 C1
  
     select {(resno=res5) and (chain=iChain) and (file=f) and (model=m)}
+
     select {(resno=res5) and (chain=iChain) and thisModel}
 
     move_it(as, vs)
 
     move_it(as, vs)
 
     force_p_res(res3, jChain)
 
     force_p_res(res3, jChain)
 
     move_it(as, vs)
 
     move_it(as, vs)
     fix_p_res(res3, jChain, TRUE)
+
     ##fix_p_res(res3, jChain, true)
 
}
 
}
  
# Flatstack res3 on res5 moving just res5
+
# Flatstack res5 on res3 moving just res5
function single_flatstack_res3_on_res5(res5, res3, iChain, jChain) {
+
function single_flatstack_res5_on_res3(res3, res5, iChain, jChain) {
    var f = (_frameID/1000000)
 
    var m = (_frameID%1000000)
 
 
     var as = gen_as(res5, res3, iChain, jChain)
 
     var as = gen_as(res5, res3, iChain, jChain)
  
 
     vs = array()
 
     vs = array()
     vs[1] = 6.00   #4.6# distance iN-jN
+
     vs[1] = 6.00     # distance res5 N9or1 and res3 N9or1     
     vs[2] = 90#75.1   # angle iN-jN-jC1
+
     vs[2] = 90#75.1   # angle res5 N9or1 and res3 N9or1 C1     
     vs[3] = 90#135.3   # dihedral iN-jN-jC1-jC4
+
     vs[3] = 90#135.3 # dihedral res5 N9or1 and res3 N9or1 C1 C4
     vs[4] = 90#89.9   # angle iC1-iN-jN
+
     vs[4] = 90#89.9   # angle res5 C1 N9or1 and res3 N9or1     
     vs[5] = -90#-47.3   # dihedral iC4-iN-iC1-jN
+
     vs[5] = -90#-47.3 # dihedral res5 C4 N9or1 C1 and res3 N9or1
     vs[6] = 0#1.7    # dihedral iN-iC1-jN-jC1
+
     vs[6] = 0#1.7    # dihedral res5 N9or1 C1 and res3 N9or1 C1
  
     select {(resno=res5) and (chain=iChain) and (file=f) and (model=m)}
+
     select {(resno=res5) and (chain=iChain) and thisModel}
 
     move_it(as, vs)
 
     move_it(as, vs)
 
     force_p_res(res5, jChain)
 
     force_p_res(res5, jChain)
Line 372: Line 405:
 
# Outstack res3 on res5 moving just res5
 
# Outstack res3 on res5 moving just res5
 
function single_outstack_res3_on_res5(res5, res3, iChain, jChain) {
 
function single_outstack_res3_on_res5(res5, res3, iChain, jChain) {
    var f = (_frameID/1000000)
 
    var m = (_frameID%1000000)
 
 
     var as = gen_as(res5, res3, iChain, jChain)
 
     var as = gen_as(res5, res3, iChain, jChain)
  
 
     var vs = array()
 
     var vs = array()
     vs[1] = 8.9    # distance iN-jN
+
     vs[1] = 8.9    # distance res5 N9or1 and res3 N9or1     
     vs[2] = 65.3  # angle iN-jN-jC1
+
     vs[2] = 65.3  # angle res5 N9or1 and res3 N9or1 C1     
     vs[3] = 55.7  # dihedral iN-jN-jC1-jC4
+
     vs[3] = 55.7  # dihedral res5 N9or1 and res3 N9or1 C1 C4
     vs[4] = 61.2  # angle iC1-iN-jN
+
     vs[4] = 61.2  # angle res5 C1 N9or1 and res3 N9or1     
     vs[5] = -41.2  # dihedral iC4-iN-iC1-jN
+
     vs[5] = -41.2  # dihedral res5 C4 N9or1 C1 and res3 N9or1
     vs[6] = -138.4 # dihedral iN-iC1-jN-jC1
+
     vs[6] = -138.4 # dihedral res5 N9or1 C1 and res3 N9or1 C1
  
     select {(resno=res5) and (chain=iChain) and (file=f) and (model=m)}
+
     select {(resno=res5) and (chain=iChain) and thisModel}
 
     move_it(as, vs)
 
     move_it(as, vs)
 
     force_p_res(res5, jChain)
 
     force_p_res(res5, jChain)
 
     move_it(as, vs)
 
     move_it(as, vs)
     fix_p_res(res5, jChain, TRUE)
+
     ##fix_p_res(res5, jChain, true)
 
}
 
}
  
# Pair A res5 on A res3 Hogsteen (N6-N7)2 moving res5 => res3
+
function make_major_groove_triplex(res5, res3, iChain, jChain) {
function pair_it_h_aa(res5, res3, ares, iChain, jChain) {
+
     var as = gen_as(res5, res3, iChain, jChain)
    var f = (_frameID/1000000)
 
    var m = (_frameID%1000000)
 
     var as = array()
 
 
     var vs = array()
 
     var vs = array()
    as[1] = atom_rcn(res3, jChain, "N6")
 
    as[2] = atom_rcn(res3, jChain, "C1\'")
 
    as[3] = atom_rcn(res3, jChain, "N7")
 
    as[4] = atom_rcn(res5, iChain, "N6")
 
    as[5] = atom_rcn(res5, iChain, "C1\'")
 
    as[6] = atom_rcn(res5, iChain, "N7")
 
    var cp = as[5].xyz
 
    select {(resno = res5) and (chain=iChain) and (file=f) and (model=m)
 
        and base} or @{as[5]}
 
 
    # Set distance of iN6 from jN7 (1tna=2.92)
 
    vs[1] = 2.92
 
 
    # Set angle of iN6 from jN7 and jC1 (1tna=123.1)
 
    vs[2] = 123.1
 
 
    # Set dihedral of iN6 from jN7 and jC1 and jN6 (1tna= 154.9)
 
    vs[3] =  154.9
 
  
     # Set angle of iN7 from iN6 and jC1 (1ana=98.2)
+
     vs[1] = 8.11  # distance res5 N9or1 and res3 N9or1
     vs[4] =  98.2
+
    vs[2] = 166.2  # angle res5 N9or1 and res3 N9or1 C1
 +
    vs[3] = 0.3    # dihedral res5 N9or1 and res3 N9or1 C1 C4
 +
     vs[4] = 162.8  # angle res5 C1 N9or1 and res3 N9or1
 +
    vs[5] = 114.6 # dihedral res5 C4 N9or1 C1 and res3 N9or1
 +
    vs[6] = 138.1  # dihedral res5 N9or1 C1 and res3 N9or1 C1
  
     # Set dihedral of iC4 from iN and iC1 and jN (1tna=18.2)
+
     # Move the nt into final position
    vs[5] = 18.2
+
    select {(resno <= res5) and (chain=iChain) and thisModel}
 
 
    # Set dihedral of iN7 from iN6 and jN7 and jC1 (1tna=177.6)
 
    vs[6] = 177.6
 
 
 
    # Move the base and C1' into final position
 
 
     move_it(as, vs)
 
     move_it(as, vs)
 +
   
 +
    # Rotate ribose to hbond O2' to res5+1 N7
 +
    select {(resno = res5) and not base
 +
        and (chain=iChain) and thisModel}
 +
    var aC1 = get_atom_rcn( res5, iChain, "C1\'")
 +
    var aN9 = get_atom_rcn( res5, iChain, "N9")
 +
    rotate selected @aC1 @aN9 -40.0
  
     # Mark C1' xyz and move it back to its orginal position
+
     # Fix up      
     var pt = as[5].xyz
+
     select {((resno=res5) or (resno=@{res5+1}))
     as[5].xyz = cp
+
        and (chain=iChain) and thisModel}
 
+
    plico_minimize( {selected})
    # Collect available P rotors
+
    fix_p_res(res5+1, iChain, true)
    var rotors = gen_nt_rotors(res5, ares, iChain)
+
}
 
 
    # Until there (4 tries)
 
    var cnt = 0
 
    while (distance(as[5], pt) > kDtolerance) {
 
 
 
        # Rotate on rotor set to move C1' to its new position
 
        move_atom_nt( as[5].atomIndex, pt, 0, rotors)
 
 
 
        # Rotate on anchor chi
 
        select {(resno >= res5) and (resno <= ares and (chain=iChain)
 
            and (file=f) and (model=m)) and not base}
 
        rotate_chi_for_distance_atoms(ares, iChain, as[5], pt, kDtolerance)
 
 
 
        cnt++
 
        if (cnt > 3) {
 
            break
 
        }
 
    }
 
  
    select {(resno=res5) and (chain=iChain) and (file=f) and (model=m)
 
        and not base}
 
    set_distance_atoms(pt, as[5] kDtolerance)
 
  
    force_p_res(ares, iChain)
+
# Pair U res5 on A res3 Hoogsteen N3-N7, N6-O2, O4-O1p moving res5 => res3
}
+
function make_hoogsteen_pair_yr(res5, res3, iChain, jChain) {
 
 
# Pair U res5 on A res3 Hogsteen N3-N7, N6-O2, O4-O1p moving res5 => res3
 
function pair_it_h_ua(res5, res3, ares, iChain, jChain) {
 
    var f = (_frameID/1000000)
 
    var m = (_frameID%1000000)
 
 
     var as = gen_as(res5, res3, iChain, jChain)
 
     var as = gen_as(res5, res3, iChain, jChain)
 
     var vs = array()
 
     var vs = array()
 
     var cp = as[6].xyz
 
     var cp = as[6].xyz
    select {(resno = res5) and (chain=iChain) and (file=f) and (model=m)
 
        and base} or @{as[6]}
 
 
    # Set distance of iN1or3 from jN1or3 (1tna=5.77)
 
    vs[1] = 5.77
 
 
    # Set angle of iN1or3 from jN1or3 and jN9or1 (1tna=115.7)
 
    vs[2] = 115.7
 
  
     # Set dihedral of iN1or3 from jN1or3 and jN9or1 and jC1' (1tna= 19.9)
+
     vs[1] = 7.05  # distance res5 N9or1 and res3 N9or1
     vs[3] = 19.9
+
    vs[2] = 150.7  # angle res5 N9or1 and res3 N9or1 C1
 +
    vs[3] = -33.1  # dihedral res5 N9or1 and res3 N9or1 C1 C4
 +
    vs[4] = 143.0  # angle res5 C1 N9or1 and res3 N9or1
 +
    vs[5] = -173.9 # dihedral res5 C4 N9or1 C1 and res3 N9or1
 +
     vs[6] = -179.8 # dihedral res5 N9or1 C1 and res3 N9or1 C1
  
     # Set angle of iN9or1 from iN1or3 and jN1or3 (1ana=57.3)
+
     # Move the nt into final position
     vs[4] =  57.3
+
     select {(resno <= res5) and (chain=iChain) and thisModel}
 
 
    # Set dihedral of iC1' from iN9or1 and iN1or3 and jN1or3 (1tna=-177.2)
 
    vs[5] = -177.2
 
 
 
    # Set dihedral of iN9or1 from iN1or3 and jN1or3 and jN9or1 (1tna=168.1)
 
    vs[6] = 168.1
 
 
 
    # Move the base and C1' into final position
 
 
     move_it(as, vs)
 
     move_it(as, vs)
  
     # Mark C1' xyz and move it back to its orginal position
+
     # Rotate 5 end out of the way
     var pt = as[6].xyz
+
    select {(resno < res5) and (chain=iChain) and thisModel}
     as[6].xyz = cp
+
    var aC4 = get_atom_rcn( res5, iChain, "C4\'")
 +
     var aC5 = get_atom_rcn( res5, iChain, "C5\'")
 +
     rotate selected @aC4 @aC5 160.0
  
     # Collect available P rotors
+
     # Fix up   
     var rotors = gen_nt_rotors(res5, ares, iChain)
+
     select {((resno=res5) or (resno=@{res5+1}))
 +
    and (chain=iChain) and thisModel}
 +
    plico_minimize( {selected})
 +
    fix_p_res(res5+1, iChain, true)
 +
    fix_p_res(res5, iChain, true)
 +
}
  
     # Until there (4 tries)
+
function level_base(rMove, rFixed, iChain, jChain)  {
     var cnt = 0
+
    var selsave = {selected}
     while (distance(as[5], pt) > kDtolerance) {
+
    var mC1 = get_atom_rcn(rMove, iChain, "C1\'")
 +
    var mIsR = {mC1 and purine}
 +
    var m9or1 = (mIsR ? "N9" : "N1")
 +
     var m6or8 = (mIsR ? "C8" : "C6")
 +
     var m4or2 = (mIsR ? "C4" : "C2")
 +
     var mN  = get_atom_rcn(rMove, iChain, m9or1)
 +
    var mC6or8  = get_atom_rcn(rMove, iChain, m6or8)
 +
    var mC4or2  = get_atom_rcn(rMove, iChain, m4or2)
 +
   
 +
    var fC1 = get_atom_rcn(rFixed, jChain, "C1\'")
 +
    var fIsR = {fC1 and purine}
 +
    var f9or1 = (fIsR ? "N9" : "N1")
 +
    var f6or8 = (fIsR ? "C8" : "C6")
 +
    var f4or2 = (fIsR ? "C4" : "C2")
 +
    var fN  = get_atom_rcn(rFixed, jChain, f9or1)
 +
    var fC6or8  = get_atom_rcn(rFixed, jChain, f6or8)
 +
    var fC4or2  = get_atom_rcn(rFixed, jChain, f4or2)
  
        # Rotate on rotor set to move C1' to its new position
+
    var dist = abs(distance(fc4or2, mc4or2) - distance(fc6or8, mc6or8))
        move_atom_nt( as[5].atomIndex, pt, 0, rotors)
+
    var newdist = dist
 
+
    var dir = 0.1
        # Rotate on anchor chi
+
    select {(resno=rMove) and (chain=iChain) and base and thisModel}
        select {(resno >= res5) and (resno <= ares) and (chain=iChain)
+
    while(newdist > 0.01) {
            and (file=f) and (model=m) and not base}
+
         if (newdist > dist) {
         rotate_chi_for_distance_atoms(ares, iChain, as[6], pt, kDtolerance)
+
            if (dir == 0.1) {
 
+
                dir = -0.1
        cnt++
+
            }
        if (cnt > 3) {
+
             else {
             break
+
                rotateSelected @mC1 @mN @{-dir}
 +
                break
 +
            }
 
         }
 
         }
 +
        dist = newdist
 +
        rotateSelected @mC1 @mN @{dir}
 +
        newdist = abs(distance(fc4or2, mc4or2) - distance(fc6or8, mc6or8))
 
     }
 
     }
 
+
     select selsave   
     select {(resno=res5) and (chain=iChain) and (file=f) and (model=m)
 
        and not base}
 
    set_distance_atoms(pt, as[6] kDtolerance)
 
 
 
    force_p_res(ares, iChain)
 
 
}
 
}
  
 
# Stack res rMove on res rFixed
 
# Stack res rMove on res rFixed
function base_stack_res( rMove, rFixed, iChain, jChain, sep , ang) {
+
function base_stack_res( rMove, rFixed, iChain, jChain, sep , ang, single) {
    var f = (_frameID/1000000)
 
    var m = (_frameID%1000000)
 
 
     var isA = is_form_a(rMove, iChain)
 
     var isA = is_form_a(rMove, iChain)
 +
    var is3on5 = (rMove > rFixed)
 
     var j = rFixed
 
     var j = rFixed
 
     var i = rMove
 
     var i = rMove
 
     var as = array()
 
     var as = array()
 
     var vs = array()
 
     var vs = array()
     as[1] = atom_rcn(j, jChain, "O3\'")
+
     as = gen_as(rMove, rFixed, iChain, jChain)
    as[2] = atom_rcn(j, jChain, "C5\'")
+
     if (single) {
     var Njx = (((as[1] and {purine}).size > 0) ? "N9" : "N1")
+
        select {(resno = i) and (chain=iChain) and thisModel}
    as[3] = atom_rcn(j, jChain, Njx,    m)
+
     }
    as[5] = atom_rcn(i, iChain, "C5\'")
+
     else {
    as[6] = atom_rcn(i, iChain, "O3\'")
+
        if (is3on5) {
     var Nix = (((as[5] and {purine}).size > 0) ? "N9" : "N1")
+
            select {(resno >= i) and (chain=iChain) and thisModel}
     as[4] = atom_rcn(i, iChain, Nix,    m)
+
        }
    select {(resno <= i) and (chain=iChain) and (file=f) and (model=m)}
+
        else {
 
+
            select {(resno <= i) and (chain=iChain) and thisModel}
     # Set distance of iNx from jNx (1tna=4.2)
+
        }
 +
    }
 +
   
 +
     # Set distance of fres N1or9 from mres N1or9 (1tna=4.2)
 
     vs[1] = sep
 
     vs[1] = sep
  
     # Set angle Njx Nix C5i (1ana=85.7)
+
     # Set angle fres C1' N1or9 and mres N1or9 (A=6tna B=1ana)
     vs[2] =  (isA ? 85.7 : 77.8)
+
     vs[2] =  (isA ? (is3on5 ? 92.8 : 113.9) : (is3on5 ? 83.3 : 115.23))#78.3 : 110.23
  
     # Set dihedral Njx Nix C5i O3i (1tna=179.9)
+
     # Set dihedral fres C4' C1' N1or9 and mres N1or9 (A=6tna B=1ana)
     vs[3] = (isA ? 179.9 : 173.3)
+
     vs[3] = (isA ? (is3on5 ? 110.0 : -71.2) : (is3on5 ? 165.92 : -28.31))
  
     # Set angle C5j Njx Nxif (1tna=112.9)
+
     # Set angle fres N1or9 and mres N1or9 C1' (A=6tna B=1ana)
     vs[4] = (isA ? 112.9 : 124.2)
+
     vs[4] = (isA ? (is3on5 ? 113.9 : 92.8) : (is3on5 ? 115.23 : 83.3))
  
     # Set dihedral O3j C5j Njx Nix (1tna= -20)
+
     # Set dihedral fres N1or9 and mres N1or9 C1' C4' (A=6tna B=1ana)
     vs[5] =  (isA ? -20 : -14.2)
+
     vs[5] =  (isA ? (is3on5 ? -71.2 : 110.0) : (is3on5 ? -28.31 : 165.92))
  
     # Set dihedral of C5j Njx Nix C5i (1tna=20)
+
     # Set dihedral of fres C5 N1or9 and mres N1or9 C5 (1tna=20)
 
     vs[6] = ang
 
     vs[6] = ang
 
 
     move_it(as, vs)
 
     move_it(as, vs)
 +
    #select {(resno=rMove) or (resno=rFixed)}
 +
    #minimize {selected}
 
     #force_p_res(i, iChain)
 
     #force_p_res(i, iChain)
     fix_p_res(i+1, iChain, TRUE)
+
     ##fix_p_res(i+1, iChain, true)
 
}
 
}
  
function single_unpair_yy( rMove, rFixed, iChain, jChain) {
 
}
 
function single_unpair_ry( rMove, rFixed, iChain, jChain) {
 
}
 
function single_unpair_yr( rMove, rFixed, iChain, jChain) {
 
}
 
function single_unpair_rr( res5, res3, ares, iChain, jChain) {
 
    var f = (_frameID/1000000)
 
    var m = (_frameID%1000000)
 
 
    # Push aside the r stacked stemward on ares
 
    var isYR5 =  (res3 = (ares-1))
 
    var as = array()
 
    var vs = array()
 
    if (isYR5) {
 
        to_ab_nt_res(res3, res3-1, jChain, FALSE)
 
        force_p_res(res3, jChain)
 
        as = gen_as(res3, ares, jChain, jChain)
 
        select {(resno=res3) and (chain=iChain) and (file=f) and (model=m)}
 
    }
 
    else {
 
        to_ab_nt_res(res5, res5-1, iChain, FALSE)
 
        force_p_res(res5, iChain)
 
        as = gen_as(res5, res5-1, iChain, iChain)
 
        select {(resno=res5) and (chain=iChain) and (file=f) and (model=m)}
 
    }
 
 
    # Set distance of iN from jN (1ana=6.14)
 
    vs[1] = 6.14
 
 
    # Set angle of iN from jN and jC1 (1ana=102.5)
 
    vs[2] = 102.5
 
 
    # Set dihedral of iN from jN and jC1 and jC4 (1ana=-66.2)
 
    vs[3] = -66.2
 
 
    # Set angle of iC1 from iN nad jN (1ana=72.5)
 
    vs[4] = 72.5
 
 
    # Set dihedral of iC4 from iN and iC1 and jN (1ana=160.0)
 
    vs[5] = 176.5
 
 
    # Set dihedral of iN from iC1 and jN and jC1 (1ana=-5.0)
 
    vs[6] = -1.3
 
 
    move_it(as, vs)
 
    if (isYR5) {
 
        force_p_res(res3+1, jChain)
 
        fix_p_res(res3, jChain, FALSE)
 
    }
 
    else {
 
        force_p_res(res5-1, iChain)
 
        fix_p_res(res5, iChain, FALSE)
 
    }
 
 
    #########
 
    # Stack the other r on ares
 
    var as = gen_as(res5, ares, iChain, jChain)
 
    select {(resno = res5) and (chain=iChain) and (file=f) and (model=m)}
 
 
    # Set distance of iN from jN (1ana=9.00)
 
    vs[1] = 9.00
 
 
    # Set angle of iN from jN and jC1 (1ana=124.6)
 
    vs[2] = 124.6
 
 
    # Set dihedral of iN from jN and jC1 and jC4 (1ana=160.0)
 
    vs[3] = 160.0
 
 
    # Set angle of iC1 from iN nad jN (1ana=124.6)
 
    vs[4] = 124.6
 
 
    # Set dihedral of iC4 from iN and iC1 and jN (1ana=160.0)
 
    vs[5] = 160.0
 
 
    # Set dihedral of iN from iC1 and jN and jC1 (1ana=-5.0)
 
    vs[6] = -5.0
 
 
    move_it(as, vs)
 
    fix_p_res(res5, iChain, FALSE)
 
}
 
  
 
# Rotate rotor set to move target atom to its proper place
 
# Rotate rotor set to move target atom to its proper place
function move_atom_nt(targetIdx, targetPt, ares, iRotors) {
+
# ares is the 5ward res limit exclusive of the mobile
     var f = (_frameID/1000000)
+
function move_atom_nt(targetIdx, targetPt, ares, rotors, force) {
    var m = (_frameID%1000000)
+
     var set3
 
     var pt = targetPt
 
     var pt = targetPt
    var rotors = iRotors
 
 
     var targetNo = {atomIndex=targetIdx}.atomno
 
     var targetNo = {atomIndex=targetIdx}.atomno
 
     var targetRes = {atomIndex=targetIdx}.resno
 
     var targetRes = {atomIndex=targetIdx}.resno
 
     var iChain = {atomIndex=targetIdx}.chain
 
     var iChain = {atomIndex=targetIdx}.chain
     gOK = FALSE
+
     gOK = false
 
     var dist = distance(pt, {atomIndex=targetIdx}.xyz)
 
     var dist = distance(pt, {atomIndex=targetIdx}.xyz)
  
Line 687: Line 609:
 
                         v2 = {atomIndex=i3}.xyz - {atomIndex=i2}.xyz
 
                         v2 = {atomIndex=i3}.xyz - {atomIndex=i2}.xyz
  
                         s = sin(abs(angle(v1, {0 0 0}, v2)))
+
                         var s = sin(abs(angle(v1, {0 0 0}, v2)))
 
                         if (s > smax) {
 
                         if (s > smax) {
 
                             smax = s
 
                             smax = s
Line 718: Line 640:
 
                 res3 = {atomIndex=i1}.resno
 
                 res3 = {atomIndex=i1}.resno
 
             }
 
             }
 +
            select remove tbase
 
             #***************************************************
 
             #***************************************************
 
             rotateSelected {atomIndex=i3} {atomIndex=i2} @{-dt}
 
             rotateSelected {atomIndex=i3} {atomIndex=i2} @{-dt}
  
 
             # If collisions
 
             # If collisions
             var res5 = res3-1
+
             set3 = (within(kCtolerance, {selected}) and not {selected}
            var set3 = {(resno=res3) and (atomName!="P") and (atomName!="OP1")
+
            and not connected({selected}))
                and (file=f) and (model=m)}
+
             if ((force == false) and (set3)) {
             var set5 = {(resno=res5) and (atomName!="P") and (atomName!="OP1")
 
                and (file=f) and (model=m)}
 
            if ((set5 and within(kCtolerance, set3)).size > 0) {
 
 
 
 
                 # Binary undo until fixed
 
                 # Binary undo until fixed
 
                 while ((abs(dt) > kDtolerance)
 
                 while ((abs(dt) > kDtolerance)
                     and ((set5 and within(kCtolerance, set3)).size > 0)) {
+
                     and ((set5 and within(kCtolerance, set3)))) {
 
                     dt /= 2.0
 
                     dt /= 2.0
 
                     rotateSelected {atomIndex=i3} {atomIndex=i2} @{dt}
 
                     rotateSelected {atomIndex=i3} {atomIndex=i2} @{dt}
 
                 }
 
                 }
 
                 while ((abs(dt) > kDtolerance)
 
                 while ((abs(dt) > kDtolerance)
                     and ((set5 and within(kCtolerance, set3)).size > 0)) {
+
                     and ((set5 and within(kCtolerance, set3)))) {
 
                     dt /= 2.0
 
                     dt /= 2.0
 
                     rotateSelected {atomIndex=i3} {atomIndex=i2} @{-dt}
 
                     rotateSelected {atomIndex=i3} {atomIndex=i2} @{-dt}
Line 746: Line 665:
 
             dist = distance(pt, {atomIndex=targetIdx})
 
             dist = distance(pt, {atomIndex=targetIdx})
 
             if (dist < kDtolerance) {
 
             if (dist < kDtolerance) {
                 gOK = TRUE
+
                 gOK = true
 
                 gTargetPt = pt
 
                 gTargetPt = pt
 
                 break
 
                 break
Line 760: Line 679:
 
         }
 
         }
 
     }  # endfor 20 passes
 
     }  # endfor 20 passes
 +
    return set3
 +
}
 +
 +
# Counter rotate rotor set to move target atom to its proper place
 +
function move_atom_by_cr_nt(targetIdx, targetPt, ares, iRotors) {
 +
    var pt = targetPt
 +
    var rotors = iRotors
 +
    var targetNo = {atomIndex=targetIdx}.atomno
 +
    var targetRes = {atomIndex=targetIdx}.resno
 +
    var iChain = {atomIndex=targetIdx}.chain
 +
    gOK = false
 +
    var dist = distance(pt, {atomIndex=targetIdx}.xyz)
 +
 +
    # If target is a C1' atom, collect its base
 +
    var tBase = ({})
 +
    var i1 = 0
 +
    var i2 = 0
 +
    var i3 = 0
 +
    var i4 = 0
 +
    if ({atomIndex=targetIdx}.atomName == "C1\'") {
 +
        tBase = {(resno = targetRes) and base}
 +
    }
 +
 +
    # For all C4'-C5' axes
 +
    for (var ri = 1; ri < rotors.size; ri += 4) {
 +
        if ({atomIndex=@{rotors[ri]}}.atomName == "C4\'") {
 +
       
 +
            # While distance lessens
 +
            var dist = distance(pt, {atomIndex=targetIdx})
 +
            var first = true
 +
            var dt = 5.0
 +
            while (dist > kDtolerance) {
 +
           
 +
                # Counter rotate C4'-C5' and O5'-P axes
 +
                var i1 = rotors[ri+8]
 +
                var i2 = rotors[ri+9]
 +
                var i3 = rotors[ri+10]
 +
                var i4 = rotors[ri+11]
 +
                var x2 = rotors[ri+17]
 +
                var x3 = rotors[ri+18]
 +
                var x4 = rotors[ri+19]
 +
                var res3 = 0
 +
               
 +
                # Select and rotate
 +
                if (ares > targetRes) {
 +
                    select_3ward_atom({atomIndex=i3}, ares, iChain)
 +
                    res3 = {atomIndex=i4}.resno
 +
                }
 +
                else {
 +
                    select_5ward_atom({atomIndex=i3}, ares, iChain)
 +
                    res3 = {atomIndex=i1}.resno
 +
                }
 +
                select remove tbase
 +
                rotateSelected {atomIndex=i3} {atomIndex=i2} @{dt}
 +
               
 +
                # Select and counter rotate
 +
                if (ares > targetRes) {
 +
                    select_3ward_atom({atomIndex=x3}, ares, iChain)
 +
                }
 +
                else {
 +
                    select_5ward_atom({atomIndex=x3}, ares, iChain)
 +
                }
 +
                select remove tbase
 +
                rotateSelected {atomIndex=x3} {atomIndex=x2} @{-dt}
 +
               
 +
                # If first and worse, reverse
 +
                var newdist = distance(pt, {atomIndex=targetIdx})
 +
                if (newdist > dist) {
 +
                    if (first) {
 +
                        dt = -dt
 +
                    }
 +
                    else {
 +
                        break
 +
                    }
 +
                }
 +
                first = false
 +
                dist = newdist
 +
 +
                /*# If collisions
 +
                var res5 = res3-1
 +
                var set3 = {(resno=res3) and (atomName!="P") and (atomName!="OP1")
 +
                    and thisModel}
 +
                var set5 = {(resno=res5) and (atomName!="P") and (atomName!="OP1")
 +
                    and thisModel}
 +
                if ((set5 and within(kCtolerance, set3))) {
 +
                    # Binary undo until fixed
 +
                    while ((abs(dt) > kDtolerance)
 +
                        and ((set5 and within(kCtolerance, set3)))) {
 +
                        dt /= 2.0
 +
                        rotateSelected {atomIndex=i3} {atomIndex=i2} @{dt}
 +
                    }
 +
                    while ((abs(dt) > kDtolerance)
 +
                        and ((set5 and within(kCtolerance, set3)))) {
 +
                        dt /= 2.0
 +
                        rotateSelected {atomIndex=i3} {atomIndex=i2} @{-dt}
 +
                    }
 +
                    rotateSelected {atomIndex=i3} {atomIndex=i2} @{dt}
 +
                }*/
 +
   
 +
            }  # endwhile
 +
   
 +
        }
 +
    }  # endfor rotors
 
}
 
}
  
 
# If ares < 0 then adjust iRes only
 
# If ares < 0 then adjust iRes only
 
function to_ab_nt_res(res, ares, iChain, toA) {
 
function to_ab_nt_res(res, ares, iChain, toA) {
     var f = (_frameID/1000000)
+
     var selsave = {selected}
    var m = (_frameID%1000000)
+
     var aO3 =  get_atom_rcn( res, iChain, "O3\'")
     var aO3 =  atom_rcn( res, iChain, "O3\'")
+
     var aC3 =  get_atom_rcn( res, iChain, "C3\'")
     var aC3 =  atom_rcn( res, iChain, "C3\'")
+
     var aC4 =  get_atom_rcn( res, iChain, "C4\'")
     var aC4 =  atom_rcn( res, iChain, "C4\'")
+
     var aC5 =  get_atom_rcn( res, iChain, "C5\'")
     var aC5 =  atom_rcn( res, iChain, "C5\'")
+
     var aC1 =  get_atom_rcn( res, iChain, "C1\'")
     var aC1 =  atom_rcn( res, iChain, "C1\'")
+
     var aC2 =  get_atom_rcn( res, iChain, "C2\'")
     var aC2 =  atom_rcn( res, iChain, "C2\'")
+
     var aO2 =  get_atom_rcn( res, iChain, "O2\'")
     var aO2 =  atom_rcn( res, iChain, "O2\'")
+
     var aO4 =  get_atom_rcn( res, iChain, "O4\'")
     var aO4 =  atom_rcn( res, iChain, "O4\'")
 
  
 
     if (ares < 0) {
 
     if (ares < 0) {
         select ((resno=res) and (chain=iChain) and (file=f) and (model=m)
+
         select ((resno=res) and (chain=iChain) and thisModel
 
             and not aO3 and not aC3 and not aC4)
 
             and not aO3 and not aC3 and not aC4)
 
     }
 
     }
 
     else {
 
     else {
         select ((resno <= ares) and (chain=iChain) and (file=f) and (model=m)
+
         select ((resno >= ares) and (resno <= res) and (chain=iChain)
             and not aO3 and not aC3 and not aC4)
+
             and thisModel and not aO3 and not aC3 and not aC4)
 
     }
 
     }
     set_dihedral_atoms(aO3, aC3, aC4, aC5, (toA ? kO3C3C4C5A : kO3C3C4C5A))
+
     set_dihedral_atoms(aO3, aC3, aC4, aC5, (toA ? kO3C3C4C5A : kO3C3C4C5B))
  
 
     # Set chi
 
     # Set chi
Line 789: Line 810:
 
     var aCx = -1
 
     var aCx = -1
 
     var ang = 0.0
 
     var ang = 0.0
     select {(resno=res) and (chain=iChain) and (file=f) and (model=m) and base}
+
     select {(resno=res) and (chain=iChain) and thisModel and base}
     if ((aC1 and {purine}).size > 0) {
+
     if (aC1 and {purine}) {
         aNx =  atom_rcn( res, iChain, "N9")
+
         aNx =  get_atom_rcn( res, iChain, "N9")
         aCx =  atom_rcn( res, iChain, "C8")
+
         aCx =  get_atom_rcn( res, iChain, "C8")
 
         ang = (toA ? kPuA : kPuB)
 
         ang = (toA ? kPuA : kPuB)
 
     }
 
     }
 
     else {
 
     else {
         aNx =  atom_rcn(res, iChain, "N1")
+
         aNx =  get_atom_rcn(res, iChain, "N1")
         aCx =  atom_rcn(res, iChain, "C6")
+
         aCx =  get_atom_rcn(res, iChain, "C6")
 
         ang = (toA ? kPyA : kPyB)
 
         ang = (toA ? kPyA : kPyB)
 
     }
 
     }
Line 805: Line 826:
 
     var pSet = {aC1 or aC2 or aO2}
 
     var pSet = {aC1 or aC2 or aO2}
 
     select pSet or {(resno=res) and (chain=iChain)
 
     select pSet or {(resno=res) and (chain=iChain)
         and (file=f) and (model=m) and base}
+
         and thisModel and base}
 
     set_dihedral_atoms(aC5, aC4, aO4, aC1, (toA ? kC5C4O4C1A : kC5C4O4C1B))
 
     set_dihedral_atoms(aC5, aC4, aO4, aC1, (toA ? kC5C4O4C1A : kC5C4O4C1B))
    set_dihedral_atoms(aC4, aO4, aC1, aNx, (toA ? kC4O4C1NxA : kC4O4C1NxB))
 
 
     set_dihedral_atoms(aC4, aO4, aC1, aC2, (toA ? kC4O4C1C2A : kC4O4C1C2B))
 
     set_dihedral_atoms(aC4, aO4, aC1, aC2, (toA ? kC4O4C1C2A : kC4O4C1C2B))
     if (aO2.size > 0) {
+
    set_dihedral_atoms(aC2, aO4, aC1, aNx, (toA ? kC2O4C1NxA : kC2O4C1NxB))
 +
     if (aO2) {
 
         select aO2 or aC2
 
         select aO2 or aC2
 
         ang = (toA ? kC3C1C2O2A : kC3C1C2O2B)
 
         ang = (toA ? kC3C1C2O2A : kC3C1C2O2B)
Line 816: Line 837:
 
     set_distance_atoms(aC3, aC2, 1.52)
 
     set_distance_atoms(aC3, aC2, 1.52)
 
     set_distance_atoms(aC1, aC2, 1.52)
 
     set_distance_atoms(aC1, aC2, 1.52)
 +
    select selsave
 
}
 
}
  
 
function adjust_nts(res5, res3, iChain, toab, a, s) {
 
function adjust_nts(res5, res3, iChain, toab, a, s) {
    var savemt = useMinimizationThread
 
    set useMinimizationThread FALSE
 
  
 
     # Collect any pairing
 
     # Collect any pairing
Line 830: Line 850:
 
     # Twist and turn
 
     # Twist and turn
 
     for (var i = res3; i >= res5; i--) {
 
     for (var i = res3; i >= res5; i--) {
         if (toab.size > 0) {
+
        var j = i-res5+1
 +
         if (toab) {
 
             to_ab_nt_res(i, -1, iChain, (toab == "A"))
 
             to_ab_nt_res(i, -1, iChain, (toab == "A"))
             if ((w[i])[1] >= 0) {
+
             if ((w[j])[1] >= 0) {
                 to_ab_nt_res((w[i])[1], -1, (w[i])[2], (toab == "A"))
+
                 to_ab_nt_res((w[j])[1], -1, (w[j])[2], (toab == "A"))
 
             }
 
             }
 
         }
 
         }
        if (i < res3) {
+
    }
            base_stack_res(i, i+1, iChain, iChain, s, a)
+
    for (var i = res5; i < res3; i++) {
        }
+
        base_stack_res(i, i+1, iChain, iChain, s, a, false)
 
     }
 
     }
  
 
     # Restore pairings
 
     # Restore pairings
 
     for (var i = res3; i >= res5; i--) {
 
     for (var i = res3; i >= res5; i--) {
         if ((w[i])[1] >= 0) {
+
        var j = i-res5+1
             pair_it_res((w[i])[1], i, -1, (w[i])[2], iChain)
+
         if ((w[j])[1] >= 0) {
 +
             pair_it_res((w[j])[1], i, -1, (w[j])[2], iChain)
 
         }
 
         }
 
     }
 
     }
Line 850: Line 872:
 
     # Clean up
 
     # Clean up
 
     for (var i = res3; i >= res5; i--) {
 
     for (var i = res3; i >= res5; i--) {
         fix_p_res(i, iChain, TRUE)
+
        var j = i-res5+1
         if ((w[i])[1] >= 0) {
+
         fix_p_res(i, iChain, true)
             fix_p_res((w[i])[1], (w[i])[2], TRUE)
+
         if ((w[j])[1] >= 0) {
 +
             fix_p_res((w[j])[1], (w[j])[2], true)
 
         }
 
         }
 
     }
 
     }
    set useMinimizationThread savemt
 
 
}
 
}
  
 
#########################################################
 
#########################################################
### STAND ALONE GENERAL PURPOSE FUNCTIONS                   ###
+
### STAND ALONE GENERAL PURPOSE FUNCTIONS             ###
 
#########################################################
 
#########################################################
 
function is_form_a( iResno, iChain) {
 
function is_form_a( iResno, iChain) {
     var aO4 = atom_rcn( iResno, iChain, "O4\'")
+
     var aO4 = get_atom_rcn( iResno, iChain, "O4\'")
     var aC1 = atom_rcn( iResno, iChain, "C1\'")
+
     var aC1 = get_atom_rcn( iResno, iChain, "C1\'")
     var aC2 = atom_rcn( iResno, iChain, "C2\'")
+
     var aC2 = get_atom_rcn( iResno, iChain, "C2\'")
     var aC3 = atom_rcn( iResno, iChain, "C3\'")
+
     var aC3 = get_atom_rcn( iResno, iChain, "C3\'")
 
     return (angle(aO4, aC1, aC2, aC3) < 0.0)
 
     return (angle(aO4, aC1, aC2, aC3) < 0.0)
 
}
 
}
  
 
function is_r_res( iResno, iChain) {
 
function is_r_res( iResno, iChain) {
    var f = (_frameID/1000000)
+
     return ({(resno=iResno) and (chain=iChain) and thisModel and purine})
    var m = (_frameID%1000000)
 
     return ({(resno=iResno) and (chain=iChain) and (file=f)
 
        and (model=m) and purine}.size > 0)
 
 
}
 
}
  
 
function repair_p_res(res, iChain) {
 
function repair_p_res(res, iChain) {
     var aP = atom_rcn( res, iChain, "P")
+
     var aP = get_atom_rcn( res, iChain, "P")
     minimize select {connected(aP) or aP}
+
     plico_minimize( {connected(aP) or aP})
 
}
 
}
  
 
function who_pairs(iRes, iChain) {
 
function who_pairs(iRes, iChain) {
     var aC4or6 =  atom_rcn( iRes, iChain, "C4")
+
     var aC4or6 =  get_atom_rcn( iRes, iChain, "C4")
     var aN1or3 =  atom_rcn( iRes, iChain, "N1")
+
     var aN1or3 =  get_atom_rcn( iRes, iChain, "N1")
 
     if ({aN1or3 and purine}.size = 0) {
 
     if ({aN1or3 and purine}.size = 0) {
         aC4or6 =  atom_rcn( iRes, iChain, "C6")
+
         aC4or6 =  get_atom_rcn( iRes, iChain, "C6")
         aN1or3 =  atom_rcn( iRes, iChain, "N3")
+
         aN1or3 =  get_atom_rcn( iRes, iChain, "N3")
 
     }
 
     }
     var near = within(3.1, aN1or3) and {resno!=iRes} and {element="N"}
+
    if (aN1or3) {
    for (var i = 1; i <= near.size; i++) {
+
     var near = within(3.2, aN1or3) and {resno!=iRes} and {element="N"}
        if (angle(near[i], aN1or3, aC4or6) > 150) {
+
        for (var i = 1; i <= near.size; i++) {
            return [near[i].resno, near[i].chain]
+
            var dist = distance(near[i], aN1or3)
 +
            var ang = abs(angle(near[i], aN1or3, aC4or6))
 +
            if (ang > 150) {
 +
                return [near[i].resno, near[i].chain, dist, ang]
 +
            }
 
         }
 
         }
 
     }
 
     }
     return [-1, aC4or6.chain]
+
     return [-1, aC4or6.chain, -1, -1]
 +
}
 +
 
 +
function who_almost_pairs(iRes, iChain) {
 +
    var aC2 =  get_atom_rcn( iRes, iChain, "C2")
 +
    var aC4or6 =  get_atom_rcn( iRes, iChain, "C4")
 +
    var aN1or3 =  get_atom_rcn( iRes, iChain, "N1")
 +
    var pname = "C6"
 +
    if ({aN1or3 and purine}.size = 0) {
 +
        aC4or6 =  get_atom_rcn( iRes, iChain, "C6")
 +
        aN1or3 =  get_atom_rcn( iRes, iChain, "N3")
 +
        pname = "C4"
 +
    }
 +
    if (aN1or3) {
 +
        var near = within(3.4, aN1or3) and {resno!=iRes} and {element="N"}
 +
        for (var i = 1; i <= near.size; i++) {
 +
            var aC6or4 = get_atom_rcn(near[i].resno, near[i].chain, pname)
 +
            var aC2p = get_atom_rcn(near[i].resno, near[i].chain, "C2")
 +
            var dist = distance(near[i], aN1or3)
 +
            var puang = abs(angle(near[i], aN1or3, aC4or6))
 +
            var pyang = abs(angle(aC6or4, near[i], aN1or3))
 +
            var dihedral = abs(angle(aC2p, near[i], aN1or3, aC2))
 +
            if ((puang > 110) and (pyang > 110) and (dihedral < 40)) {
 +
                return [near[i].resno, near[i].chain, dist, puang]
 +
            }
 +
        }
 +
    }
 +
    return [-1, aC4or6.chain, -1, -1]
 
}
 
}
  
Line 916: Line 966:
 
             var a4 = angle(bset[2], bset[3], bset[4])
 
             var a4 = angle(bset[2], bset[3], bset[4])
  
             var isStacked = TRUE
+
             var isStacked = true
  
 
             # Bases are parallel as sin(a1) = sin(a2) and sin(a3) = sin(a4)
 
             # Bases are parallel as sin(a1) = sin(a2) and sin(a3) = sin(a4)
 
             if (abs(sin(a1)-sin(a2)) > 20) {
 
             if (abs(sin(a1)-sin(a2)) > 20) {
                 isStacked = FALSE
+
                 isStacked = false
 
             }
 
             }
 
             if (abs(sin(a3)-sin(a4)) > 20) {
 
             if (abs(sin(a3)-sin(a4)) > 20) {
                 isStacked = FALSE
+
                 isStacked = false
 
             }
 
             }
  
 
             # Bases are stacked as d*sin(a1) < 6.0 and d3 = 0.0
 
             # Bases are stacked as d*sin(a1) < 6.0 and d3 = 0.0
 
             if (d*sin(a1) > 6.2) {
 
             if (d*sin(a1) > 6.2) {
                 isStacked = FALSE
+
                 isStacked = false
 
             }
 
             }
 
             if (abs(dh) > 30) {
 
             if (abs(dh) > 30) {
                 #isStacked = FALSE
+
                 #isStacked = false
 
             }
 
             }
  
Line 945: Line 995:
  
 
function match_nt(mask, nt) {
 
function match_nt(mask, nt) {
     var ret = FALSE
+
     var ret = false
 
     switch (mask) {
 
     switch (mask) {
 
     case "A":
 
     case "A":
Line 954: Line 1,004:
 
         break
 
         break
 
     case "N":
 
     case "N":
         ret = TRUE
+
         ret = true
 +
        break
 +
    case "M":
 +
        ret = ((nt=="A") or (nt=="C"))
 
         break
 
         break
 
     case "Y":
 
     case "Y":
Line 967: Line 1,020:
  
 
# Calls function match_nt above
 
# Calls function match_nt above
function select_seqs(seq, iChain) {
+
function select_seqs(seq, r5, r3, iChain, f, m) {
    var f = (_frameID/1000000)
 
    var m = (_frameID%1000000)
 
 
     select none
 
     select none
     for (var i = {(chain=iChain) and (file=f) and (model=m)}.resno.min;
+
     for (var i = r5; i < r3; i++) {
        i <= {(chain=iChain) and (file=f) and (model=m)}.resno.max; i++) {
+
         var j = 0
         var j = 1
+
         for (; j < seq.size; j++) {
         for (; j <= seq.size; j++) {
+
             if ((gNTres[i+j])[2] >= 0) {
             var nt = {(chain=iChain) and (file=f) and (model=m) and (resno=@{i+j-1})
+
                break
                 and (atomName="C1\'")}.group[1]
+
            }
             if (not match_nt(seq[j], nt)) { # <== external call
+
            if ((i+j) >= r3) {
 +
                 break
 +
            }
 +
             if (not match_nt(seq[j+1], gSeq[i+j])) { # CALL
 
                 break
 
                 break
 
             }
 
             }
 
         }
 
         }
         if (j > seq.size) {
+
         if (j == seq.size) {
             print format("%s at %d (%s-%s-%s)", seq, i,
+
             print format("%s at %d (%s-%s-%s)%s", seq, i,
 
                 gSeq[i-1],
 
                 gSeq[i-1],
                 gSeq[i][i+seq.size-1],
+
                 gSeq[i][i+j-1],
                 gSeq[i+seq.size])
+
                 gSeq[i+j], aster)
             var rset = {(resno=i) and (chain=iChain) and (file=f) and (model=m)}
+
             var rset = {(resno=i) and (chain=iChain) and thisModel}
             rset.selected = TRUE
+
            rset.selected = true
 +
        }
 +
    }
 +
}
 +
 
 +
function find_tetras(seq, r5, r3, based) {
 +
    select none
 +
    for (var i = r5; i < r3; i++) {
 +
        var j = 0
 +
        for (; j < seq.size; j++) {
 +
            if ((i+j) >= r3) {
 +
                break
 +
            }
 +
             if (not match_nt(seq[j+1], gSeq[i+j])) { # CALL
 +
                break
 +
            }
 +
        }
 +
        if (j == seq.size) {
 +
            if (based) {
 +
                var ends = gSeq[i-1] + gSeq[i+j]
 +
                var bads = ["AA","GG","AG","GA"]
 +
                if (bads.count(ends) > 0) {
 +
                    continue
 +
                }
 +
            }
 +
            print format("%s at %d (%s-%s-%s)%s", seq, i,
 +
                gSeq[i-1],
 +
                gSeq[i][i+j-1],
 +
                gSeq[i+j], aster)
 
         }
 
         }
 
     }
 
     }
Line 994: Line 1,076:
 
# Calls is_form_a
 
# Calls is_form_a
 
function select_b_form_nts(iChain) {
 
function select_b_form_nts(iChain) {
    var f = (_frameID/1000000)
 
    var m = (_frameID%1000000)
 
 
     select none
 
     select none
     for (var i = {(chain=iChain) and (file=f) and (model=m)}.resno.min;
+
     for (var i = get_resno_min(iChain); i <= get_resno_max(iChain); i++) {
        i <= {(chain=iChain) and (file=f) and (model=m)}.resno.max; i++) {
 
 
         if (not is_form_a(i, iChain)) { # <== external call
 
         if (not is_form_a(i, iChain)) { # <== external call
 
             print format("Res %d is form B", i)
 
             print format("Res %d is form B", i)
             var rset = {(resno=i) and (chain=iChain) and (file=f) and (model=m)}
+
             var rset = {(resno=i) and (chain=iChain) and thisModel}
             rset.selected = TRUE
+
             rset.selected = true
        }
 
    }
 
}
 
 
 
# From 2LU0
 
function make_uncg_loop(res5, ares, iChain) {
 
    to_ab_nt_res(res5+1, ares, iChain, FALSE)
 
    to_ab_nt_res(res5+2, ares, iChain, FALSE)
 
    var vs = [
 
        [144.8,  28.8, 138.6,  140.9, 149.6]
 
        [-85.1, -148.4,  56.3, -143.6, 135.7]
 
        [-145.9, -25.6, -93.0, -148.1,  49.6]
 
        [-84.8,  145.9,  -7.0, -130.0, -176.2]
 
        [-133.5,  -78.0, -60.2,  132.6,  99.0]]
 
    make_tetra_loop(res5, ares, iChain, m, vs) # <== external call
 
}
 
 
 
# From 2LU0
 
function make_gnra_loop(res5, ares, iChain) {
 
    to_ab_nt_res(res5+1, ares, iChain, FALSE)
 
    var vs = [
 
        [147.0,  23.9,  144.6,  144.5, 151.0]
 
        [-118.2, -90.5,  109.0,  166.6, 100.1]
 
        [-155.6, -34.4, -167.7,  115.4, 131.2]
 
        [-136.5, -69.9,  -86.8, -170.6,  56.0]
 
        [-147.1, -57.5,  -76.8,  147.7,  94.7]]
 
    make_tetra_loop(res5, ares, iChain, m, vs) # <== external call
 
}
 
 
 
function make_tetra_loop(res5, ares, iChain, m, vs) {
 
    var f = (_frameID/1000000)
 
    var m = (_frameID%1000000)
 
    for (var i = (res5+4); i >= ires5; i--) {
 
        var j = i-1
 
        var as = array()
 
        as += atom_rcn( j, iChain, "C4\'")
 
        as += atom_rcn( j, iChain, "C3\'")
 
        as += atom_rcn( j, iChain, "O3\'")
 
        as += atom_rcn( i, iChain, "P")
 
        as += atom_rcn( i, iChain, "O5\'")
 
        as += atom_rcn( i, iChain, "C5\'")
 
        as += atom_rcn( i, iChain, "C4\'")
 
        as += atom_rcn( i, iChain, "C3\'")
 
        as += atom_rcn( i, iChain, "OP1")
 
        as += atom_rcn( i, iChain, "OP2")
 
        for (var k = 5; k > 0; k--) {
 
            pset = ((k>2) ? (connected(@{as[4]}) or @{as[4]}) : ({}))
 
            select {((resno<i) and (resno>ares) and (chain=iChain)
 
                and (file=f) and (model=m))
 
                or pset}
 
            set_dihedral_atoms(as[k+3], as[k+2], as[k+1], as[k],
 
                (vs[i-res5+1])[k])
 
 
         }
 
         }
 
     }
 
     }
Line 1,060: Line 1,087:
  
 
function select_3ward_atom(ar3, ares, iChain) {
 
function select_3ward_atom(ar3, ares, iChain) {
    var f = (_frameID/1000000)
 
    var m = (_frameID%1000000)
 
 
     var i = ar3.resno
 
     var i = ar3.resno
     var aP = atom_rcn( i, iChain, "P")
+
     var aP = get_atom_rcn( i, iChain, "P")
 
     switch(ar3.atomName) {
 
     switch(ar3.atomName) {
 
     case "O3\'" :
 
     case "O3\'" :
 
         select {(resno>i) and (resno<ares) and (chain=iChain)
 
         select {(resno>i) and (resno<ares) and (chain=iChain)
         and (file=f) and (model=m)}
+
         and thisModel}
 
         break
 
         break
 
     case "P" :
 
     case "P" :
 
         select {(resno>=i) and (resno<ares) and (chain=iChain)
 
         select {(resno>=i) and (resno<ares) and (chain=iChain)
             and (file=f) and (model=m)}
+
             and thisModel}
 
         break
 
         break
 
     case "O5\'" :
 
     case "O5\'" :
Line 1,077: Line 1,102:
 
     case "C4\'" :
 
     case "C4\'" :
 
         select {(resno>=i) and (resno<ares) and (chain=iChain)
 
         select {(resno>=i) and (resno<ares) and (chain=iChain)
             and (file=f) and (model=m) and not (connected(aP) or aP)}
+
             and thisModel and not (connected(aP) or aP)}
 
         break
 
         break
 
     case "C3\'" :
 
     case "C3\'" :
         var aO3 = atom_rcn( i, iChain, "O3\'")
+
         var aO3 = get_atom_rcn( i, iChain, "O3\'")
 
         select {((resno>i) and (resno<ares) and (chain=iChain)
 
         select {((resno>i) and (resno<ares) and (chain=iChain)
             and (file=f) and (model=m)) or aO3}
+
             and thisModel) or aO3}
 
         break
 
         break
 
     }
 
     }
Line 1,088: Line 1,113:
  
 
function select_5ward_atom(ar5, ares, iChain) {
 
function select_5ward_atom(ar5, ares, iChain) {
    var f = (_frameID/1000000)
 
    var m = (_frameID%1000000)
 
 
     var i = ar5.resno
 
     var i = ar5.resno
     var aP = atom_rcn( i, iChain, "P")
+
     var aP = get_atom_rcn( i, iChain, "P")
 
     switch(ar5.atomName) {
 
     switch(ar5.atomName) {
 
     case "O3\'" :
 
     case "O3\'" :
 
         select {(resno<=i) and (resno>ares) and (chain=iChain)
 
         select {(resno<=i) and (resno>ares) and (chain=iChain)
             and (file=f) and (model=m)}
+
             and thisModel}
 
         break
 
         break
 
     case "P" :
 
     case "P" :
Line 1,101: Line 1,124:
 
     case "C5\'" :
 
     case "C5\'" :
 
         select {((resno<i) and (resno>ares) and (chain=iChain)
 
         select {((resno<i) and (resno>ares) and (chain=iChain)
             and (file=f) and (model=m)) or (connected(aP) or aP)}
+
             and thisModel) or (connected(aP) or aP)}
 
         break
 
         break
 
     case "C4\'" :
 
     case "C4\'" :
         var aC5 = atom_rcn( i, iChain, "C5\'")
+
         var aC5 = get_atom_rcn( i, iChain, "C5\'")
 
         select {((resno<i) and (resno>ares) and (chain=iChain)
 
         select {((resno<i) and (resno>ares) and (chain=iChain)
             and (file=f) and (model=m)) or (connected(aP) or aP or aC5)}
+
             and thisModel) or (connected(aP) or aP or aC5)}
 
         break
 
         break
 
     }
 
     }
 
}
 
}
 +
 +
function plot_ab_chi( iChain) {
 +
    select none
 +
    for (var i = get_resno_min(iChain); i <= get_resno_max(iChain); i++) {
 +
        var aO4 = get_atom_rcn(i, iChain, "O4\'")
 +
        var aC1 = get_atom_rcn(i, iChain, "C1\'")
 +
        var isR = {aC1 and purine}
 +
        var a1or9 = (isR ? "N9" : "N1")
 +
        var a6or8 = (isR ? "C8" : "C6")
 +
        var aN  = get_atom_rcn(i, iChain, a1or9)
 +
        var aC  = get_atom_rcn(i, iChain, a6or8)
 +
        var aO3 = get_atom_rcn(i, iChain, "O3\'")
 +
        var aC3 = get_atom_rcn(i, iChain, "C3\'")
 +
        var aC4 = get_atom_rcn(i, iChain, "C4\'")
 +
        var aC5 = get_atom_rcn(i, iChain, "C5\'")
 +
       
 +
        var chi = angle(aO4, aC1, aN, aC)
 +
        aN.vx = chi
 +
        var aorb = angle(aO3, aC3, aC4, aC5)
 +
        aN.vy = aorb
 +
        select ADD aN
 +
    }
 +
    plot properties vx vy resno
 +
    set echo top left
 +
    echo "vx = base chi angle        vy = a ==> b form"
 +
}
 +
 +
function move_p_to_close_o3( aP, aO3, ares, force) {
 +
    var cp = aP.xyz
 +
    var aC3 = {connected(aO3) and (atomname="C3\'")}
 +
    select aP
 +
    set_distance_atoms(aO3, aP, 1.74)
 +
    set_angle_atoms(aC3, aO3, aP, 109.0)
 +
    var pt = aP.xyz
 +
    aP.xyz = cp
 +
    var rotors = gen_nt_rotors(aP.resno, ares, aP.chain)
 +
    move_atom_nt(aP.atomIndex, pt, aP.resno-1, rotors, force)
 +
    ##fix_p_res(aP.resno, aP.chain, true)
 +
}
 +
 +
function move_o3_to_close_p( aO3, aP, ares, force) {
 +
    var cp = aO3.xyz
 +
    var aO5 = {connected(aP) and (atomname="O5\'")}
 +
    select aO3
 +
    set_distance_atoms(aP, aO3, 1.74)
 +
    set_angle_atoms(aO5, aP, aO3, 109.0)
 +
    var pt = aO3.xyz
 +
    aO3.xyz = cp
 +
    var rotors = gen_nt_rotors(ares, aO3.resno, aO3.chain)
 +
    move_atom_nt(aO3.atomIndex, pt, aO3.resno+1, rotors, force)
 +
    ##fix_p_res(aP.resno, aP.chain, true)
 +
}
 +
 +
# Select mobile before calling
 +
function pivot_to_close_atoms( aMov, aStat, aPivot, dist) {
 +
    var caxis = cross(aStat.xyz, aMov.xyz) - aPivot.xyz
 +
    var dir = 1
 +
    var d = distance(aMov, aStat)
 +
    if (d > dist) {
 +
        while (d > dist) {
 +
            rotateSelected @aPivot @caxis @dir
 +
            var nd = distance(aMov, aStat)
 +
            if (nd > d) {
 +
                rotateSelected @aPivot @caxis @{-dir}
 +
                if (dir == 1) {
 +
                    rotateSelected @aPivot @caxis @{-dir}
 +
                    dir = -1
 +
                    continue
 +
                }
 +
                else {
 +
                    break
 +
                }
 +
            }
 +
            d = nd
 +
        }
 +
    }
 +
    else {
 +
        while (d < dist) {
 +
            rotateSelected @aPivot @caxis @dir
 +
            var nd = distance(aMov, aStat)
 +
            if (nd < d) {
 +
                rotateSelected @aPivot @caxis @{-dir}
 +
                if (dir == 1) {
 +
                    rotateSelected @aPivot @caxis @{-dir}
 +
                    dir = -1
 +
                    continue
 +
                }
 +
                else {
 +
                    break
 +
                }
 +
            }
 +
            d = nd
 +
        }
 +
    }
 +
}
 +
 +
function print_adjacent_vs( iChain) {
 +
    var tNN = 0
 +
    var taA = 0
 +
    var tbA = 0
 +
    var taD = 0
 +
    var tbD = 0
 +
    var tabD = 0
 +
    var tc = 0.0
 +
    var rmin = get_resno_min(iChain)
 +
    var rmax = get_resno_max(iChain)
 +
    for (var i = rmin; i < rmax; i++) {
 +
        var aC4 = get_atom_rcn(i, iChain, "C4\'")
 +
        var aC1 = get_atom_rcn(i, iChain, "C1\'")
 +
        var isR = {aC1 and purine}
 +
        var a1or9 = (isR ? "N9" : "N1")
 +
        var aN  = get_atom_rcn(i, iChain, a1or9)
 +
       
 +
        var aC4p = get_atom_rcn(i+1, iChain, "C4\'")
 +
        var aC1p = get_atom_rcn(i+1, iChain, "C1\'")
 +
        isR = {aC1p and purine}
 +
        a1or9 = (isR ? "N9" : "N1")
 +
        var aNp  = get_atom_rcn(i+1, iChain, a1or9)
 +
       
 +
        if (aC4 and aC4p) {
 +
           
 +
            var NN = distance(aN, aNp)
 +
            tNN += NN
 +
            var aA = angle(aC1, aN, aNp)
 +
            taA += aA
 +
            var bA = angle(aC1p, aNp, aN)
 +
            tbA += bA
 +
           
 +
            var aD = angle(aC4, aC1, aN, aNp)
 +
            if ((aD < 0) and (taD > 0)) {
 +
                aD += 360
 +
            }
 +
            taD += aD
 +
            var bD = angle(aC4p, aC1p, aNp, aN)
 +
            if ((bD < 0) and (tbD > 0)) {
 +
                bD += 360
 +
            }
 +
            tbD += bD
 +
            var abD = angle(aC1, aN, aNp, aC1p)
 +
            tabD += abD
 +
            tc++
 +
            print format("%6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %d%s %s",
 +
                NN, aA, aD, bA, bD, abD, tc, aC1.group, aC1p.group)
 +
        }
 +
    }
 +
   
 +
    print format("v1=%6.2f v2=%6.2f v3=%6.2f v4=%6.2f v5=%6.2f v6=%6.2f",
 +
        tNN/tc, taA/tc, taD/tc, tbA/tc, tbD/tc, tabD/tc)
 +
}
 +
 +
function measure_adjacent_vs(r5, iChain) {
 +
    var aC4 = get_atom_rcn(r5, iChain, "C4\'")
 +
    var aC1 = get_atom_rcn(r5, iChain, "C1\'")
 +
    var isR = {aC1 and purine}
 +
    var aN1or9  = get_atom_rcn(r5, iChain, (isR ? "N9" : "N1"))
 +
    var aC8or6 =  get_atom_rcn(r5, iChain, (isR ? "C8" : "C6"))
 +
   
 +
    var aC4p = get_atom_rcn(r5+1, iChain, "C4\'")
 +
    var aC1p = get_atom_rcn(r5+1, iChain, "C1\'")
 +
    isR = {aC1p and purine}
 +
    var aN1or9p  = get_atom_rcn(r5+1, iChain, (isR ? "N9" : "N1"))
 +
    var aC8or6p =  get_atom_rcn(r5+1, iChain, (isR ? "C8" : "C6"))
 +
 +
    measure @aN1or9 @aN1or9p
 +
    measure @aC1 @aN1or9 @aN1or9p       
 +
    measure @aC4 @aC1 @aN1or9 @aN1or9p       
 +
    measure @aC1p @aN1or9p @aN1or9       
 +
    measure @aC4p @aC1p @aN1or9p @aN1or9       
 +
    measure @aC1p @aN1or9p @aN1or9 @aC1     
 +
   
 +
    measure @aC4 @aC1 @aN1or9 @aC8or6       
 +
    measure @aC4p @aC1p @aN1or9p @aC8or6p       
 +
}
 +
 +
function print_pair_vs( iChain) {
 +
    var tNN = 0.0
 +
    var taA = 0.0
 +
    var tbA = 0.0
 +
    var taD = 0.0
 +
    var tbD = 0.0
 +
    var tabD = 0.0
 +
    var tc = 0.0
 +
    var taChi = 0.0
 +
    var tbChi = 0.0
 +
    var rmin = get_resno_min(iChain)
 +
    var rmax = get_resno_max(iChain)
 +
    for (var i = rmin; i < rmax; i++) {
 +
        var aC4 = get_atom_rcn(i, iChain, "C4\'")
 +
        if (aC4) {
 +
            var w = who_pairs(i, iChain) # CALL
 +
            var aC1 = get_atom_rcn(i, iChain, "C1\'")
 +
            var isR = {aC1 and purine}
 +
            var a1or9 = (isR ? "N9" : "N1")
 +
            var aN  = get_atom_rcn(i, iChain, a1or9)
 +
            var a6or8 = (isR ? "C8" : "C6")
 +
            var aC  = get_atom_rcn(i, iChain, a6or8)
 +
            var aC4p = get_atom_rcn(w[1], w[2], "C4\'")
 +
            var aC1p = get_atom_rcn(w[1], w[2], "C1\'")
 +
           
 +
            a1or9 = (isR ? "N1" : "N9") # rev
 +
            var aNp  = get_atom_rcn(w[1], w[2], a1or9)
 +
            a6or8 = (isR ? "C6" : "C8") # rev
 +
            var aCp  = get_atom_rcn(w[1], w[2], a6or8)
 +
           
 +
            var NN = distance(aN, aNp)
 +
            tNN += NN
 +
            var aA = angle(aC1, aN, aNp)
 +
            taA += aA
 +
            var bA = angle(aC1p, aNp, aN)
 +
            tbA += bA
 +
           
 +
            var aD = angle(aC4, aC1, aN, aNp)
 +
            if ((aD < 0) and (taD > 0)) {
 +
                aD += 360
 +
            }
 +
            taD += aD
 +
            var bD = angle(aC4p, aC1p, aNp, aN)
 +
            if ((bD < 0) and (tbD > 0)) {
 +
                bD += 360
 +
            }
 +
            tbD += bD
 +
            var abD = angle(aC1, aN, aNp, aC1p)
 +
            tabD += abD
 +
           
 +
            var aChi = angle(aC4, aC1, aN, aC)
 +
            taChi += aChi
 +
            var bChi = angle(aC4p, aC1p, aNp, aCp)
 +
            tbChi += bChi
 +
           
 +
            tc++
 +
            print format(
 +
                "%6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %d%s %d%s",
 +
                NN, aA, aD, bA, bD, abD, aChi, bChi, tc, aC1.group, w[1],
 +
                aC1p.group)
 +
        }
 +
    }
 +
   
 +
    print format(
 +
        "v1=%6.2f v2=%6.2f v3=%6.2f v4=%6.2f v5=%6.2f v6=%6.2f v7=%6.2f v8=%6.2f",
 +
        tNN/tc, taA/tc, taD/tc, tbA/tc, tbD/tc, tabD/tc, taChi/tc, tbChi/tc)
 +
}
 +
 +
function measure_pair_vs(i, j, iChain, jChain) {
 +
    var aC4 = get_atom_rcn(i, iChain, "C4\'")
 +
    var aC4p = get_atom_rcn(j, jChain, "C4\'")
 +
    if (aC4.size and aC4p.size) {
 +
        var aC1 = get_atom_rcn(i, iChain, "C1\'")
 +
        var isR = {aC1 and purine}
 +
        var a1or9 = (isR ? "N9" : "N1")
 +
        var aN  = get_atom_rcn(i, iChain, a1or9)
 +
        var a6or8 = (isR ? "C8" : "C6")
 +
        var aC  = get_atom_rcn(i, iChain, a6or8)
 +
        var aC1p = get_atom_rcn(j, jChain, "C1\'")
 +
        isR = {aC1p and purine}
 +
        a1or9 = (isR ? "N9" : "N1")
 +
        var aNp  = get_atom_rcn(j, jChain, a1or9)
 +
        a6or8 = (isR ? "C8" : "C6")
 +
        var aCp  = get_atom_rcn(j, jChain, a6or8)
 +
   
 +
        measure @aN @aNp
 +
        measure @aC1p @aNp @aN       
 +
        measure @aC4p @aC1p @aNp @aN       
 +
        measure @aC1 @aN @aNp       
 +
        measure @aC4 @aC1 @aN @aNp       
 +
        measure @aC1p @aNp @aN @aC1
 +
        measure @aC4p @aC1p @aNp @aCp
 +
        measure @aC4 @aC1 @aN @aC
 +
    }
 +
    else {
 +
        print "No pair found"
 +
    }     
 +
}
 +
 +
function minimize_for_collision( r, iChain) {
 +
    fix_p_res( r, iChain, true)
 +
 +
    # External
 +
    var cset = (within(kCtolerance, {resno=r}) and not {resno=r} and not
 +
        connected({resno=r}) and {(chain=iChain) and thisModel})
 +
    if (cset) {
 +
        for (var i = 0; i < cset.size; i++) {
 +
            plico_minimize( {(resno=r) or (resno=@{cset[i].resno})})
 +
        }
 +
    }
 +
 +
    # Internal
 +
    cset = (within(kCtolerance, {(resno=r) and (atomName="OP?")}) and not
 +
        {(atomName="OP?") or (atomName="P")}
 +
        and {(chain=iChain) and thisModel})
 +
    if (cset) {
 +
        for (var i = 0; i < cset.size; i++) {
 +
            plico_minimize( {(resno=r) or (resno=@{cset[i].resno})})
 +
        }
 +
    }
 +
}
 +
 +
function eval_pairing( seq, r5, r3, len) {
 +
    var val = 0
 +
    for (var i = 0; i < len; i++) {
 +
        var c5 = seq[r5+i]
 +
        var c3 = seq[r3-i]
 +
        if ((c5 == "A") and (c3 == "U")) {
 +
            val += 2
 +
        }
 +
        else if ((c5 == "U") and (c3 == "A")) {
 +
            val += 2
 +
        }
 +
        else if ((c5 == "G") and (c3 == "C")) {
 +
            val += 3
 +
        }
 +
        else if ((c5 == "C") and (c3 == "G")) {
 +
            val += 3
 +
        }
 +
        else if ((c5 == "G") and (c3 == "U")) {
 +
            val += 2
 +
        }
 +
        else if ((c5 == "U") and (c3 == "G")) {
 +
            val += 1
 +
        }
 +
        else if ((c5 == "A") and (c3 == "A")) {
 +
            val -= 3
 +
        }
 +
        else if ((c5 == "G") and (c3 == "G")) {
 +
            val -= 3
 +
        }
 +
    }
 +
    return val
 +
}
 +
 +
function update_atomnos(iChain) {
 +
    print "Update atomnos..."
 +
    var b = {(chain=iChain) and thisModel}.atomno.min
 +
    for (var r = get_resno_min(iChain); r <= get_resno_max(iChain); r++) {
 +
        var rset = {(resno=r) and (chain=iChain) and thisModel}
 +
        for (var n = rset.atomno.min; n <= rset.atomno.max; n++) {
 +
            var a = {(resno=r) and (chain=iChain) and thisModel and (atomno=n)}
 +
            a.atomno = -b
 +
            b++
 +
        }
 +
    }
 +
    for (var i = -b; i < 0; i++) {
 +
        var a = {(chain=iChain) and thisModel
 +
            and (atomno=i)}
 +
        a.atomno *= -1
 +
    }
 +
    print "Updated"   
 +
}
 +
 +
function replane_bases(iChain) {
 +
    for (var i = get_resno_min(iChain); i <= get_resno_max(iChain); i++) {
 +
        plico_minimize( {resno=i and base})
 +
    }
 +
}
 +
 +
function measure_p_dihedrals(i, iChain) {
 +
        var j = i-1
 +
        var as = array()
 +
        as += get_atom_rcn( j, iChain, "C4\'")
 +
        as += get_atom_rcn( j, iChain, "C3\'")
 +
        as += get_atom_rcn( j, iChain, "O3\'")
 +
        as += get_atom_rcn( i, iChain, "P")
 +
        as += get_atom_rcn( i, iChain, "O5\'")
 +
        as += get_atom_rcn( i, iChain, "C5\'")
 +
        as += get_atom_rcn( i, iChain, "C4\'")
 +
        as += get_atom_rcn( i, iChain, "C3\'")
 +
        for (var k = 1; k <= (as.size-3); k++) {
 +
            measure @{as[k+0]} @{as[k+1]} @{as[k+2]} @{as[k+3]}
 +
        }
 +
}
 +
 +
function pair_stem( r5, r3, ar5, iChain) {
 +
 +
    # Pair entire stem
 +
    var c = r3 - r5 + 1
 +
    for (var k = (c\2)-1 ; k >= 0; k--) {
 +
        pair_it_res(r5+k, r3-k, ar5, iChain, iChain) # CALL
 +
    }
 +
}
 +
 +
# From 2LU0
 +
function make_uncg_loop(res5, ares5, iChain) {
 +
    var va = array()
 +
    va[1] = [4.46, 100.5, -83.4, 114.4, 99.8, 39.2, -22.2, -11.8]
 +
    va[2] = [7.89, 93.6, 50.7, 29.7, 75.3, 163.6, 14.2 -22.2]
 +
    va[3] = [7.57, 90.9, -11.9, 59.5, 75.6, 74.1, 33.5, 14.2]
 +
    va[4] = [6.09, 124.9, -5.8, 46.4, -72.6, -27.1, -138.3, 33.5]
 +
    va[5] = [5.1, 64.0, -102.5, 98.2, 87.0, 58.6, -9.2, -138.3]
 +
    for (var i = 1; i <= 5; i++) {
 +
        var as = gen_as(res5+i-2, res5+i-1, iChain, iChain)
 +
   
 +
        var vs = array()
 +
        vs[1] = (va[i])[1]  # distance res5 N9or1 and res3 N9or1     
 +
        vs[2] = (va[i])[2]  # angle res5 N9or1 and res3 N9or1 C1     
 +
        vs[3] = (va[i])[3]  # dihedral res5 N9or1 and res3 N9or1 C1 C4
 +
        vs[4] = (va[i])[4]  # angle res5 C1 N9or1 and res3 N9or1     
 +
        vs[5] = (va[i])[5]  # dihedral res5 C4 N9or1 C1 and res3 N9or1
 +
        vs[6] = (va[i])[6]  # dihedral res5 N9or1 C1 and res3 N9or1 C1
 +
        vs[7] = (va[i])[7]  # dihedral chi res3
 +
        vs[8] = (va[i])[8]  # dihedral chi res5
 +
   
 +
        select {(resno < @{res5+i-1}) and (resno > ares5)
 +
            and (chain=iChain) and thisModel}
 +
        move_it(as, vs)
 +
        fix_p_res(res5+i-1, iChain, true)
 +
    }
 +
}
 +
 +
# From 2LU0
 +
function make_gnra_loop(res5, ares5, iChain) {
 +
    var va = array()
 +
    va[1] = [4.37, 93.9, -84.9, 117.9, 99.1, 42.0, -19.4, 3.2]
 +
    va[2] = [7.94, 77.5, 47.5, 58.9, 87.1, 179.5, 4.5 -19.4]
 +
    va[3] = [4.9, 102.9, -98.6, 79.2, 112.4, 61.2, -29.6, 4.5]
 +
    va[4] = [4.97, 124.5, -70.0, 99.8, 93.2, 34.0, 2.5, -29.6]
 +
    va[5] = [4.31, 76.0, -98.4, 106.6, 97.2, 45.6, -15.5, 2.5]
 +
    for (var i = 1; i <= 5; i++) {
 +
        var as = gen_as(res5+i-2, res5+i-1, iChain, iChain)
 +
   
 +
        var vs = array()
 +
        vs[1] = (va[i])[1]  # distance res5 N9or1 and res3 N9or1     
 +
        vs[2] = (va[i])[2]  # angle res5 N9or1 and res3 N9or1 C1     
 +
        vs[3] = (va[i])[3]  # dihedral res5 N9or1 and res3 N9or1 C1 C4
 +
        vs[4] = (va[i])[4]  # angle res5 C1 N9or1 and res3 N9or1     
 +
        vs[5] = (va[i])[5]  # dihedral res5 C4 N9or1 C1 and res3 N9or1
 +
        vs[6] = (va[i])[6]  # dihedral res5 N9or1 C1 and res3 N9or1 C1
 +
        vs[7] = (va[i])[7]  # dihedral chi res3
 +
        vs[8] = (va[i])[8]  # dihedral chi res5
 +
   
 +
        select {(resno < @{res5+i-1}) and (resno > ares5)
 +
            and (chain=iChain) and thisModel}
 +
        move_it(as, vs)
 +
        fix_p_res(res5+i-1, iChain, true)       
 +
    }
 +
}
 +
  
 
# end of plicoNTcommon.spt
 
# end of plicoNTcommon.spt
 
</pre>
 
</pre>

Latest revision as of 17:11, 12 April 2016

This script contains routines used by some other scripts of the Plico suite involved with polynucleotide manipulation. It must be located in the same directory as any script that uses these routines.

Copy and paste the following into a text editor and save in your scripts folder as plicoNTcommon.spt.

#   plicoNTcommon - Jmol script by Ron Mignery
#   v1.10 beta    4/12/2016 -axis is now a reserved word
#
#   Routines and values common to Plico suite scripts that work with nucleotides
#   Must be present in the same directory as other Plico scripts that use them
kNTcommon = 6
kC5O5PO3B = -71.0
kO5PO3C3B = -107.0
kPO3C3C4B = -161.5
kO3C3C4C5B = 140.0
kC3C4C5O5B = 55.65
kC4C5O5PB = 169.0

kO4C4C3C2B = 15.92
kC4O4C1C2B = -19.9 #-41.7 1bna minimized
kC2O4C1NxB = -122.6 #-159.0 1bna minimized
kC5C4O4C1B = 122.2 #146.3 1bna minimized
kC3C1C2O2B = 120.5

kPuB = 59.0
kPyB = 61.0

kC5O5PO3A = -59.3
kO5PO3C3A = -63.1
kPO3C3C4A = -157.4
kO3C3C4C5A = 75.5
kC3C4C5O5A = 49.55
kC4C5O5PA = 169.2

kO4C4C3C2A = -35.55
kC4O4C1C2A = 3.8
kC2O4C1NxA = -131.0
kC5C4O4C1A = 144.85
kC3C1C2O2A = 116.3

kPuA = 13.5
kPyA = 16.5

gChain1 = "A"
gChain2 = ""

# Select before calling
function force_p_res(cres, iChain) {
    var pres = cres-1
    var aP =   get_atom_rcn( cres, iChain, "P")
    var aO5 =  get_atom_rcn( cres, iChain, "O5\'")
    var aC5 =  get_atom_rcn( cres, iChain, "C5\'")
    var aC4 =  get_atom_rcn( cres, iChain, "C4\'")
    var aOP1 = get_atom_rcn( cres, iChain, "OP1")
    var aOP2 = get_atom_rcn( cres, iChain, "OP2")
    var aO3p = get_atom_rcn( pres, iChain, "O3\'")
    var aC3p = get_atom_rcn( pres, iChain, "C3\'")
    if (aO3p) {

        var selsave = {selected}
        set_distance_atoms(aP3p, aC5, 3.1)
        select aO5
        var dist = distance(aO3p, aO5)
        var widen = (dist < 2.85)
        var dir = (widen ? -1 : 1)
        var first = true
        while (abs(dist-2.85) > kDtolerance) {
            rotateSelected @aC4 @aC5 @dir
            var newdist = distance(aO3p, aO5)
            if (widen ? (newdist < dist) : (newdist > dist)) {
                if (first) {
                    dir = -dir
                    rotateSelected @aC4 @aC5 @dir
                }
                else {
                    break
                }
            }
            dist=newdist
            first = false
        }
        select aP
        set_distance_atoms(aO5, aP, 1.73)
        set_angle_atoms(aC5, aO5, aP, 110.1)

        #set_dihedral_atoms(aC4, aC5, aO5, aP, 150.3)
        aOP2.xyz = get_tet_idx(aO3p.atomIndex, aP.atomIndex, aO5.atomIndex, 1.73)
        aOP1.xyz = get_tet_idx(aO5.atomIndex, aP.atomIndex, aO3p.atomIndex, 1.73)
        plico_minimize( {(connected(aP) or aP) and not aO3p})
        select selsave
    }
}

function fix_p_res(cres, iChain, force) {
print format("fix_p_res(cres=%d, ichain=%s, force=%s)", cres, iChain, force)
    var pres = cres-1
    var aP   = get_atom_rcn( cres, iChain, "P")
    var aO5  = get_atom_rcn( cres, iChain, "O5\'")
    var aC5  = get_atom_rcn( cres, iChain, "C5\'")
    var aC4  = get_atom_rcn( cres, iChain, "C4\'")
    var aC1  = get_atom_rcn( cres, iChain, "C1\'")
    var aOP1 = get_atom_rcn( cres, iChain, "OP1")
    var aOP2 = get_atom_rcn( cres, iChain, "OP2")
    var aO3p = get_atom_rcn( pres, iChain, "O3\'")
    var aC3p = get_atom_rcn( pres, iChain, "C3\'")
    var aC4p  = get_atom_rcn( pres, iChain, "C4\'")
    if (aO3p.size and aC4.size) {
        var selsave = {selected}

        # If collision
        if (force and distance(aC3p, aC5) <= kCtolerance) {
            # Push away
            select {(resno <= @{aC5.resno}) and (chain=iChain)
                and thisModel}
            set_distance_atoms(aC3p, aC5, kCtolerance)
        }

        # Rotate C4'-C5' until P-O3' is 1.59
        select aP
        set_distance_atoms(aO5, aP, 1.59)
        set_angle_atoms(aC5, aO5, aP, 109)
        set_dihedral_atoms(aC4, aC5, aO5, aP, 180)
        select add aO5
        var dist = distance(aP, aO3p)
        if (dist > 1.59) {
            var dir = 1.0
            for (var i = 0; i < 180; i++) {
                dist = distance(aP, aO3p)
                if ((dist-1.59) < 0.1) {
                    break
                }
                rotateSelected @aC4 @aC5 @dir
                var newdist = distance(aP, aO3p)
                if (newDist > dist) {
                    rotateSelected @aC4 @aC5 @{-dir}
                    if (dir > 0) {
                        dir = -dir
                    }
                    else {
                        break
                    }
                }
            } #endfor 180
        }
        else {
            var dir = -1.0
            for (var i = 0; i < 180; i++) {
                dist = distance(aP, aO3p)
                if ((1.59-dist) < 0.1) {
                    break
                }
                rotateSelected @aC4 @aC5 @dir
                var newdist = distance(aP, aO3p)
                if (newDist < dist) {
                    rotateSelected @aC4 @aC5 @{-dir}
                    if (dir < 0) {
                        dir = -dir
                    }
                    else {
                        break
                    }
                }
            } #endfor 180
        }
        aOP2.xyz = get_tet_idx(aO3p.atomIndex, aP.atomIndex, aO5.atomIndex, 1.73)
        aOP1.xyz = get_tet_idx(aO5.atomIndex, aP.atomIndex, aO3p.atomIndex, 1.73)
        if (force) {
            select {aP or aOP1 or aOP2}
            plico_minimize( {selected})
        }
        select selsave
    }
}

function fix_p_res_range(res5, res3, iChain, force) {
    for (var i = res5; i <= res3; i++) {
        fix_p_res(i, iChain, force)
    }
}

function fix_all_nt_collisions( iChain) {
    var selsave = {selected}
    chset = count_collisions(true)
    for (var i = 1; i <= chset.size; i++) {
        c = chset[i]
        cset = (within(kCtolerance, c) and  not c
            and not connected(chset[i]))
        rset = [{resno=@{c.resno}}]
        for (var j = 1; j <= cset.size; j++) {
            rset += {resno=@{cset[j].resno}}
        }
        
        if ({c and cset and not base}) {
            select {@{rset} and base}
            plico_minimize( {selected})
        }
        else if (c.atomname[1][2] == "OP") {
            fix_p_res(chset[i].resno, iChain, true)
        }
        else {
            plico_minimize(rset)
        }
    }
    measure off
    select selsave
}

# The following functions position one nt relative to another:
# Common positioning functions:
function get_rotors_res(res) {
    var rotors = array()
    var sRes = res
    var mRes = sRes-1
    var iChain = {(resno=res) and (atomName="P")
        and thisModel}.chain
    var mC4 = get_atom_rcn( mRes, iChain, "C4\'")
    var mC3 = get_atom_rcn( mRes, iChain, "C3\'")
    var mO3 = get_atom_rcn( mRes, iChain, "O3\'")
    var sP =  get_atom_rcn( sRes, iChain, "P"   )
    var sO5 = get_atom_rcn( sRes, iChain, "O5\'")
    var sC5 = get_atom_rcn( sRes, iChain, "C5\'")
    var sC4 = get_atom_rcn( sRes, iChain, "C4\'")
    var sC3 = get_atom_rcn( sRes, iChain, "C3\'")

    rotors += [mC4.atomIndex, mC3.atomIndex, mO3.atomIndex, sP.atomIndex]
    rotors += [mC3.atomIndex, mO3.atomIndex, sP.atomIndex, sO5.atomIndex]
    rotors += [mO3.atomIndex, sP.atomIndex, sO5.atomIndex, sC5.atomIndex]
    rotors += [sP.atomIndex, sO5.atomIndex, sC5.atomIndex, sC4.atomIndex]
    rotors += [sO5.atomIndex, sC5.atomIndex, sC4.atomIndex, sC3.atomIndex]
    return rotors
}

function get_nt_chi_rotor_res(res, iChain) {
    var rotors = array()
    var aO4 = get_atom_rcn( res, iChain, "O4\'")
    var aC1 = get_atom_rcn( res, iChain, "C1\'")
    var isR  = (aC1 and {purine})
    var N1or9 = (isR ? "N9" : "N1")
    var C6or8 = (isR ? "C8" : "C6")

    var aN = get_atom_rcn(res, iChain, N1or9)
    var aC = get_atom_rcn(res, iChain, C6or8)

    rotors = [aO4.atomIndex, aC1.atomIndex, aN.atomIndex, aC.atomIndex]
    return rotors
}

function get_nt_ab_rotor_res(res, iChain) {
    var rotors = array()
    var aC5 = get_atom_rcn(res, iChain, "C5\'")
    var aC4 = get_atom_rcn(res, iChain, "C4\'")
    var aC3 = get_atom_rcn(res, iChain, "C3\'")
    var aO3 = get_atom_rcn(res, iChain, "O3\'")

    rotors = [aO3.atomIndex, aC3.atomIndex, aC4.atomIndex, aC5.atomIndex]
    return rotors
}

function gen_nt_rotors(res5, res3, iChain) {
   var rotors = array()
    for (var i = res5+1; i <= res3; i++) {
        rotors += get_rotors_res(i, iChain)
    }
    return rotors
}

# ri=moved rj=fixed
function position_nt_by_vs(ri, rj, ares, vs, iChain, jChain) {
    if (ri > rj) {
        var as = gen_as(ri, rj, iChain, jChain)
        select {(resno < ares) and (resno >= ri)
            and (chain=iChain) and thisModel}
    }
    else {
        var as = gen_as(rj, ri, jChain, iChain)
        select {(resno > ares) and (resno <= ri)
            and (chain=iChain) and thisModel}
    }
    move_it(as, vs)
}
                
# Moved object must be selected, fixed object not
# as[6] = fixed[1-3] moved[4-6] chis [7-8]
# vs[6] = [distance(as[3-4]), angle(as[2-4]),
#  dihedral(as[1-4]), angle(as[5-3], dihedral(as[6-3],
#  dihedral(as[2-5]
function move_it(as, vs) {

    # Distance, angle, dihedral positions atom[4] to a point
    set_distance_atoms(as[3], as[4], vs[1])
    set_angle_atoms(as[2], as[3], as[4], vs[2])
    set_dihedral_atoms(as[1], as[2], as[3], as[4], vs[3])

    # Angle and dihedral orients atom[4]'s object
    set_angle_atoms(as[3], as[4], as[5], vs[4])
    set_dihedral_atoms(as[3], as[4], as[5], as[6], vs[5])

    # Dihedral sets TBD
    set_dihedral_atoms(as[2], as[3], as[4], as[5], vs[6])
    
    # If chis
    if (vs.size > 6) {
        var sel = {selected}
        select {(resno=@{as[4].resno}) and base}
        set_dihedral_atoms(as[6], as[5], as[4], as[8], vs[8])
    }
}

# ri=moved rj=fixed
function gen_as(ri, rj, iChain, jChain) {
    var as = array()
    as[1] = get_atom_rcn(rj, jChain, "C4\'")
    as[2] = get_atom_rcn(rj, jChain, "C1\'")
    as[3] = connected(as[2]) and {element="N"}
    as[5] = get_atom_rcn(ri, iChain, "C1\'")
    as[6] = get_atom_rcn(ri, iChain, "C4\'")
    as[4] = connected(as[5]) and {element="N"}
    
    as[7] = get_atom_rcn(rj, jChain,
        ((as[1] and {purine}) ? "C8" : "C6"))
    as[8] = get_atom_rcn(ri, iChain,
        ((as[6] and {purine}) ? "C8" : "C6"))
    return as
}

# Specific positioning functions:
# Pair res5 on res3 moving res5 <= res3
function pair_it_res(res5, res3, ares, iChain, jChain) {
    var as = gen_as(res5, res3, iChain, jChain)
    var isA = is_form_a(res5, iChain)
    var vs = array()
    vs[1] = 8.83                    # distance res5 N9or1 and res3 N9or1      
    vs[2] = 126.12                  # angle res5 N9or1 and res3 N9or1 C1      
    vs[3] = (isA ? 160.0 : -134.97) # dihedral res5 N9or1 and res3 N9or1 C1 C4
    vs[4] = 125.32                  # angle res5 C1 N9or1 and res3 N9or1      
    vs[5] = (isA ? 160.0 : -141.46) # dihedral res5 C4 N9or1 C1 and res3 N9or1
    vs[6] = (isA ? -5.0 : -17.87)   # dihedral res5 N9or1 C1 and res3 N9or1 C1
    vs[7] = (isA ? -20.0 : 38)      # dihedral chi res3
    vs[8] = (isA ? -20.0 : 38)      # dihedral chi res5
    if (ares < 0) {
        select ((resno=res5) and (chain=iChain) and thisModel)
    }
    else if (ares > 0) {
        select ((resno <= ares) and (chain=iChain) and thisModel)
    }
    move_it(as, vs)
    fix_p_res(res5, iChain, true)
    fix_p_res(res5+1, iChain, true)
}

# Flatstack res5 on res3 moving just res5
function single_flatstack_res5_on_res3(res5, res3, iChain, jChain) {
    var as = gen_as(res5, res3, iChain, jChain)

    var vs = array()
    vs[1] = 7.00    # distance res5 N9or1 and res3 N9or1      
    vs[2] = 89.1    # angle res5 N9or1 and res3 N9or1 C1      
    vs[3] = -49.9   # dihedral res5 N9or1 and res3 N9or1 C1 C4
    vs[4] = 83.4    # angle res5 C1 N9or1 and res3 N9or1      
    vs[5] = 125.7   # dihedral res5 C4 N9or1 C1 and res3 N9or1
    vs[6] = 5.8     # dihedral res5 N9or1 C1 and res3 N9or1 C1

    select {(resno=res5) and (chain=iChain) and thisModel}
    move_it(as, vs)
    #force_p_res(res3, jChain)
    #move_it(as, vs)
    #fix_p_res(res3, jChain, true)
    #force_p_res(res3, jChain)
}

# Outstack res5 on res3 moving just res5
function single_outstack_res5_on_res3(res5, res3, iChain, jChain) {
    var as = gen_as(res5, res3, iChain, jChain)

    var vs = array()
    vs[1] = 8.23   # distance res5 N9or1 and res3 N9or1      
    vs[2] = 32.4   # angle res5 N9or1 and res3 N9or1 C1      
    vs[3] = -26.8  # dihedral res5 N9or1 and res3 N9or1 C1 C4
    vs[4] = 99.6   # angle res5 C1 N9or1 and res3 N9or1      
    vs[5] = 57.4   # dihedral res5 C4 N9or1 C1 and res3 N9or1
    vs[6] = 179.1  # dihedral res5 N9or1 C1 and res3 N9or1 C1

    select {(resno=res5) and (chain=iChain) and thisModel}
    move_it(as, vs)
    force_p_res(res3, jChain)
    move_it(as, vs)
    ##fix_p_res(res3, jChain, true)
}

# Flatstack res5 on res3 moving just res5
function single_flatstack_res5_on_res3(res3, res5, iChain, jChain) {
    var as = gen_as(res5, res3, iChain, jChain)

    vs = array()
    vs[1] = 6.00      # distance res5 N9or1 and res3 N9or1      
    vs[2] = 90#75.1   # angle res5 N9or1 and res3 N9or1 C1      
    vs[3] = 90#135.3  # dihedral res5 N9or1 and res3 N9or1 C1 C4
    vs[4] = 90#89.9   # angle res5 C1 N9or1 and res3 N9or1      
    vs[5] = -90#-47.3 # dihedral res5 C4 N9or1 C1 and res3 N9or1
    vs[6] = 0#1.7     # dihedral res5 N9or1 C1 and res3 N9or1 C1

    select {(resno=res5) and (chain=iChain) and thisModel}
    move_it(as, vs)
    force_p_res(res5, jChain)
}

# Outstack res3 on res5 moving just res5
function single_outstack_res3_on_res5(res5, res3, iChain, jChain) {
    var as = gen_as(res5, res3, iChain, jChain)

    var vs = array()
    vs[1] = 8.9    # distance res5 N9or1 and res3 N9or1      
    vs[2] = 65.3   # angle res5 N9or1 and res3 N9or1 C1      
    vs[3] = 55.7   # dihedral res5 N9or1 and res3 N9or1 C1 C4
    vs[4] = 61.2   # angle res5 C1 N9or1 and res3 N9or1      
    vs[5] = -41.2  # dihedral res5 C4 N9or1 C1 and res3 N9or1
    vs[6] = -138.4 # dihedral res5 N9or1 C1 and res3 N9or1 C1

    select {(resno=res5) and (chain=iChain) and thisModel}
    move_it(as, vs)
    force_p_res(res5, jChain)
    move_it(as, vs)
    ##fix_p_res(res5, jChain, true)
}

function make_major_groove_triplex(res5, res3, iChain, jChain) {
    var as = gen_as(res5, res3, iChain, jChain)
    var vs = array()

    vs[1] = 8.11   # distance res5 N9or1 and res3 N9or1
    vs[2] = 166.2  # angle res5 N9or1 and res3 N9or1 C1
    vs[3] = 0.3    # dihedral res5 N9or1 and res3 N9or1 C1 C4
    vs[4] = 162.8  # angle res5 C1 N9or1 and res3 N9or1
    vs[5] = 114.6  # dihedral res5 C4 N9or1 C1 and res3 N9or1
    vs[6] = 138.1  # dihedral res5 N9or1 C1 and res3 N9or1 C1

    # Move the nt into final position
    select {(resno <= res5) and (chain=iChain) and thisModel}
    move_it(as, vs)
    
    # Rotate ribose to hbond O2' to res5+1 N7
    select {(resno = res5) and not base
        and (chain=iChain) and thisModel}
    var aC1 = get_atom_rcn( res5, iChain, "C1\'")
    var aN9 = get_atom_rcn( res5, iChain, "N9")
    rotate selected @aC1 @aN9 -40.0

    # Fix up     
    select {((resno=res5) or (resno=@{res5+1}))
        and (chain=iChain) and thisModel}
    plico_minimize( {selected})
    fix_p_res(res5+1, iChain, true)
}


# Pair U res5 on A res3 Hoogsteen N3-N7, N6-O2, O4-O1p moving res5 => res3
function make_hoogsteen_pair_yr(res5, res3, iChain, jChain) {
    var as = gen_as(res5, res3, iChain, jChain)
    var vs = array()
    var cp = as[6].xyz

    vs[1] = 7.05   # distance res5 N9or1 and res3 N9or1
    vs[2] = 150.7  # angle res5 N9or1 and res3 N9or1 C1
    vs[3] = -33.1  # dihedral res5 N9or1 and res3 N9or1 C1 C4
    vs[4] = 143.0  # angle res5 C1 N9or1 and res3 N9or1
    vs[5] = -173.9 # dihedral res5 C4 N9or1 C1 and res3 N9or1
    vs[6] = -179.8 # dihedral res5 N9or1 C1 and res3 N9or1 C1

    # Move the nt into final position
    select {(resno <= res5) and (chain=iChain) and thisModel}
    move_it(as, vs)

    # Rotate 5 end out of the way
    select {(resno < res5) and (chain=iChain) and thisModel}
    var aC4 = get_atom_rcn( res5, iChain, "C4\'")
    var aC5 = get_atom_rcn( res5, iChain, "C5\'")
    rotate selected @aC4 @aC5 160.0

    # Fix up     
    select {((resno=res5) or (resno=@{res5+1}))
     and (chain=iChain) and thisModel}
    plico_minimize( {selected})
    fix_p_res(res5+1, iChain, true)
    fix_p_res(res5, iChain, true)
}

function level_base(rMove, rFixed, iChain, jChain)  {
    var selsave = {selected}
    var mC1 = get_atom_rcn(rMove, iChain, "C1\'")
    var mIsR = {mC1 and purine}
    var m9or1 = (mIsR ? "N9" : "N1")
    var m6or8 = (mIsR ? "C8" : "C6")
    var m4or2 = (mIsR ? "C4" : "C2")
    var mN  = get_atom_rcn(rMove, iChain, m9or1)
    var mC6or8  = get_atom_rcn(rMove, iChain, m6or8)
    var mC4or2  = get_atom_rcn(rMove, iChain, m4or2)
    
    var fC1 = get_atom_rcn(rFixed, jChain, "C1\'")
    var fIsR = {fC1 and purine}
    var f9or1 = (fIsR ? "N9" : "N1")
    var f6or8 = (fIsR ? "C8" : "C6")
    var f4or2 = (fIsR ? "C4" : "C2")
    var fN  = get_atom_rcn(rFixed, jChain, f9or1)
    var fC6or8  = get_atom_rcn(rFixed, jChain, f6or8)
    var fC4or2  = get_atom_rcn(rFixed, jChain, f4or2)

    var dist = abs(distance(fc4or2, mc4or2) - distance(fc6or8, mc6or8))
    var newdist = dist
    var dir = 0.1
    select {(resno=rMove) and (chain=iChain) and base and thisModel}
    while(newdist > 0.01) {
        if (newdist > dist) {
            if (dir == 0.1) {
                dir = -0.1
            }
            else {
                rotateSelected @mC1 @mN @{-dir}
                break
            }
        }
        dist = newdist
        rotateSelected @mC1 @mN @{dir}
        newdist = abs(distance(fc4or2, mc4or2) - distance(fc6or8, mc6or8))
    }
    select selsave    
}

# Stack res rMove on res rFixed
function base_stack_res( rMove, rFixed, iChain, jChain, sep , ang, single) {
    var isA = is_form_a(rMove, iChain)
    var is3on5 = (rMove > rFixed)
    var j = rFixed
    var i = rMove
    var as = array()
    var vs = array()
    as = gen_as(rMove, rFixed, iChain, jChain)
    if (single) {
        select {(resno = i) and (chain=iChain) and thisModel}
    }
    else {
        if (is3on5) {
            select {(resno >= i) and (chain=iChain) and thisModel}
        }
        else {
            select {(resno <= i) and (chain=iChain) and thisModel}
        }
    }
    
    # Set distance of fres N1or9 from mres N1or9 (1tna=4.2)
    vs[1] = sep

    # Set angle fres C1' N1or9 and mres N1or9 (A=6tna B=1ana)
    vs[2] =  (isA ? (is3on5 ? 92.8 : 113.9) : (is3on5 ? 83.3 : 115.23))#78.3 : 110.23

    # Set dihedral fres C4' C1' N1or9 and mres N1or9 (A=6tna B=1ana)
    vs[3] = (isA ? (is3on5 ? 110.0 : -71.2) : (is3on5 ? 165.92 : -28.31))

    # Set angle fres N1or9 and mres N1or9 C1' (A=6tna B=1ana)
    vs[4] = (isA ? (is3on5 ? 113.9 : 92.8) : (is3on5 ? 115.23 : 83.3))

    # Set dihedral fres N1or9 and mres N1or9 C1' C4' (A=6tna B=1ana)
    vs[5] =  (isA ? (is3on5 ? -71.2 : 110.0) : (is3on5 ? -28.31 : 165.92))

    # Set dihedral of fres C5 N1or9 and mres N1or9 C5 (1tna=20)
    vs[6] = ang
    move_it(as, vs)
    #select {(resno=rMove) or (resno=rFixed)}
    #minimize {selected} 
    #force_p_res(i, iChain)
    ##fix_p_res(i+1, iChain, true)
}


# Rotate rotor set to move target atom to its proper place
# ares is the 5ward res limit exclusive of the mobile 
function move_atom_nt(targetIdx, targetPt, ares, rotors, force) {
    var set3
    var pt = targetPt
    var targetNo = {atomIndex=targetIdx}.atomno
    var targetRes = {atomIndex=targetIdx}.resno
    var iChain = {atomIndex=targetIdx}.chain
    gOK = false
    var dist = distance(pt, {atomIndex=targetIdx}.xyz)

    # If target is a C1' atom, collect its base
    var tBase = ({})
    var i1 = 0
    var i2 = 0
    var i3 = 0
    var i4 = 0
    if ({atomIndex=targetIdx}.atomName == "C1\'") {
        tBase = {(resno = targetRes) and base}
    }

    # For idx number of passes
    for (var pass1 = 0; pass1 < 20; pass1++) {
        var blocked = ({})
        for (var pass2 = 0; pass2 < (rotors.size/4); pass2++) {

            var v1 = {atomIndex=targetIdx}.xyz - pt

            # Find the most orthgonal unused rotor
            var imax = 0
            var smax = 0.5
            for (var ri = 1; ri < rotors.size; ri += 4) {
                i2 = rotors[ri+1]
                i3 = rotors[ri+2]
                i4 = rotors[ri+3]
                if ((i2 != targetIdx) and (i3 != targetIdx) and (i4 != targetIdx)) {
                    if ({blocked and {atomIndex=i2}}.count == 0) {
                        v2 = {atomIndex=i3}.xyz - {atomIndex=i2}.xyz

                        var s = sin(abs(angle(v1, {0 0 0}, v2)))
                        if (s > smax) {
                            smax = s
                            imax = ri
                        }
                    }
                }
            }

            # If no more rotors, break to next full try
            if (imax == 0) {
               break
            }
            i1 = rotors[imax+0]
            i2 = rotors[imax+1]
            i3 = rotors[imax+2]
            i4 = rotors[imax+3]

            # Get dihedral of rotor with target point
            var dt = (angle({atomIndex=targetIdx}, {atomIndex=i2},
                {atomIndex=i3}, pt)/(rotors.size/20))

            # Select and rotate
            if (ares > targetRes) {
                select_3ward_atom({atomIndex=i3}, ares, iChain)
                res3 = {atomIndex=i4}.resno
            }
            else {
                select_5ward_atom({atomIndex=i3}, ares, iChain)
                res3 = {atomIndex=i1}.resno
            }
            select remove tbase
            #***************************************************
            rotateSelected {atomIndex=i3} {atomIndex=i2} @{-dt}

            # If collisions
            set3 = (within(kCtolerance, {selected}) and not {selected}
             and not connected({selected}))
            if ((force == false) and (set3)) {
                # Binary undo until fixed
                while ((abs(dt) > kDtolerance)
                    and ((set5 and within(kCtolerance, set3)))) {
                    dt /= 2.0
                    rotateSelected {atomIndex=i3} {atomIndex=i2} @{dt}
                }
                while ((abs(dt) > kDtolerance)
                    and ((set5 and within(kCtolerance, set3)))) {
                    dt /= 2.0
                    rotateSelected {atomIndex=i3} {atomIndex=i2} @{-dt}
                }
                rotateSelected {atomIndex=i3} {atomIndex=i2} @{dt}
            }

            # If close enough, stop
            dist = distance(pt, {atomIndex=targetIdx})
            if (dist < kDtolerance) {
                gOK = true
                gTargetPt = pt
                break
            }

            # Block rotor
            blocked |= {atomIndex=i2}

        }   # endfor num rotors passes

        if (gOK) {
            break
        }
    }   # endfor 20 passes
    return set3
}

# Counter rotate rotor set to move target atom to its proper place
function move_atom_by_cr_nt(targetIdx, targetPt, ares, iRotors) {
    var pt = targetPt
    var rotors = iRotors
    var targetNo = {atomIndex=targetIdx}.atomno
    var targetRes = {atomIndex=targetIdx}.resno
    var iChain = {atomIndex=targetIdx}.chain
    gOK = false
    var dist = distance(pt, {atomIndex=targetIdx}.xyz)

    # If target is a C1' atom, collect its base
    var tBase = ({})
    var i1 = 0
    var i2 = 0
    var i3 = 0
    var i4 = 0
    if ({atomIndex=targetIdx}.atomName == "C1\'") {
        tBase = {(resno = targetRes) and base}
    }

    # For all C4'-C5' axes
    for (var ri = 1; ri < rotors.size; ri += 4) {
        if ({atomIndex=@{rotors[ri]}}.atomName == "C4\'") {
        
            # While distance lessens
            var dist = distance(pt, {atomIndex=targetIdx})
            var first = true
            var dt = 5.0
            while (dist > kDtolerance) {
            
                # Counter rotate C4'-C5' and O5'-P axes
                var i1 = rotors[ri+8]
                var i2 = rotors[ri+9]
                var i3 = rotors[ri+10]
                var i4 = rotors[ri+11]
                var x2 = rotors[ri+17]
                var x3 = rotors[ri+18]
                var x4 = rotors[ri+19]
                var res3 = 0
                
                # Select and rotate
                if (ares > targetRes) {
                    select_3ward_atom({atomIndex=i3}, ares, iChain)
                    res3 = {atomIndex=i4}.resno
                }
                else {
                    select_5ward_atom({atomIndex=i3}, ares, iChain)
                    res3 = {atomIndex=i1}.resno
                }
                select remove tbase
                rotateSelected {atomIndex=i3} {atomIndex=i2} @{dt}
                
                # Select and counter rotate
                if (ares > targetRes) {
                    select_3ward_atom({atomIndex=x3}, ares, iChain)
                }
                else {
                    select_5ward_atom({atomIndex=x3}, ares, iChain)
                }
                select remove tbase
                rotateSelected {atomIndex=x3} {atomIndex=x2} @{-dt}
                
                # If first and worse, reverse
                var newdist = distance(pt, {atomIndex=targetIdx})
                if (newdist > dist) {
                    if (first) {
                        dt = -dt
                    }
                    else {
                        break
                    }
                }
                first = false
                dist = newdist

                /*# If collisions
                var res5 = res3-1
                var set3 = {(resno=res3) and (atomName!="P") and (atomName!="OP1")
                    and thisModel}
                var set5 = {(resno=res5) and (atomName!="P") and (atomName!="OP1")
                    and thisModel}
                if ((set5 and within(kCtolerance, set3))) {
                    # Binary undo until fixed
                    while ((abs(dt) > kDtolerance)
                        and ((set5 and within(kCtolerance, set3)))) {
                        dt /= 2.0
                        rotateSelected {atomIndex=i3} {atomIndex=i2} @{dt}
                    }
                    while ((abs(dt) > kDtolerance)
                        and ((set5 and within(kCtolerance, set3)))) {
                        dt /= 2.0
                        rotateSelected {atomIndex=i3} {atomIndex=i2} @{-dt}
                    }
                    rotateSelected {atomIndex=i3} {atomIndex=i2} @{dt}
                }*/
    
            }   # endwhile
    
        }
    }   # endfor rotors
}

# If ares < 0 then adjust iRes only
function to_ab_nt_res(res, ares, iChain, toA) {
    var selsave = {selected}
    var aO3 =  get_atom_rcn( res, iChain, "O3\'")
    var aC3 =  get_atom_rcn( res, iChain, "C3\'")
    var aC4 =  get_atom_rcn( res, iChain, "C4\'")
    var aC5 =  get_atom_rcn( res, iChain, "C5\'")
    var aC1 =  get_atom_rcn( res, iChain, "C1\'")
    var aC2 =  get_atom_rcn( res, iChain, "C2\'")
    var aO2 =  get_atom_rcn( res, iChain, "O2\'")
    var aO4 =  get_atom_rcn( res, iChain, "O4\'")

    if (ares < 0) {
        select ((resno=res) and (chain=iChain) and thisModel
            and not aO3 and not aC3 and not aC4)
    }
    else {
        select ((resno >= ares) and (resno <= res) and (chain=iChain)
            and thisModel and not aO3 and not aC3 and not aC4)
    }
    set_dihedral_atoms(aO3, aC3, aC4, aC5, (toA ? kO3C3C4C5A : kO3C3C4C5B))

    # Set chi
    var aNx = -1
    var aCx = -1
    var ang = 0.0
    select {(resno=res) and (chain=iChain) and thisModel and base}
    if (aC1 and {purine}) {
        aNx =  get_atom_rcn( res, iChain, "N9")
        aCx =  get_atom_rcn( res, iChain, "C8")
        ang = (toA ? kPuA : kPuB)
    }
    else {
        aNx =  get_atom_rcn(res, iChain, "N1")
        aCx =  get_atom_rcn(res, iChain, "C6")
        ang = (toA ? kPyA : kPyB)
    }
    set_dihedral_atoms(aO4, aC1, aNx, aCx, ang)

    # Set pucker 3' endo or 2' endo
    var pSet = {aC1 or aC2 or aO2}
    select pSet or {(resno=res) and (chain=iChain)
        and thisModel and base}
    set_dihedral_atoms(aC5, aC4, aO4, aC1, (toA ? kC5C4O4C1A : kC5C4O4C1B))
    set_dihedral_atoms(aC4, aO4, aC1, aC2, (toA ? kC4O4C1C2A : kC4O4C1C2B))
    set_dihedral_atoms(aC2, aO4, aC1, aNx, (toA ? kC2O4C1NxA : kC2O4C1NxB))
    if (aO2) {
        select aO2 or aC2
        ang = (toA ? kC3C1C2O2A : kC3C1C2O2B)
        set_dihedral_atoms(aC3, aC1, aC2, aO2, (toA ? kC3C1C2O2A : kC3C1C2O2B))
    }
    set_distance_atoms(aC3, aC2, 1.52)
    set_distance_atoms(aC1, aC2, 1.52)
    select selsave
}

function adjust_nts(res5, res3, iChain, toab, a, s) {

    # Collect any pairing
    var w = array()
    for (var i = res5; i <= res3; i++) {
        w = w + [who_pairs(i, iChain)]
    }

    # Twist and turn
    for (var i = res3; i >= res5; i--) {
        var j = i-res5+1
        if (toab) {
            to_ab_nt_res(i, -1, iChain, (toab == "A"))
            if ((w[j])[1] >= 0) {
                to_ab_nt_res((w[j])[1], -1, (w[j])[2], (toab == "A"))
            }
        }
    }
    for (var i = res5; i < res3; i++) {
        base_stack_res(i, i+1, iChain, iChain, s, a, false)
    }

    # Restore pairings
    for (var i = res3; i >= res5; i--) {
        var j = i-res5+1
        if ((w[j])[1] >= 0) {
            pair_it_res((w[j])[1], i, -1, (w[j])[2], iChain)
        }
    }

    # Clean up
    for (var i = res3; i >= res5; i--) {
        var j = i-res5+1
        fix_p_res(i, iChain, true)
        if ((w[j])[1] >= 0) {
            fix_p_res((w[j])[1], (w[j])[2], true)
        }
    }
}

#########################################################
### STAND ALONE GENERAL PURPOSE FUNCTIONS             ###
#########################################################
function is_form_a( iResno, iChain) {
    var aO4 = get_atom_rcn( iResno, iChain, "O4\'")
    var aC1 = get_atom_rcn( iResno, iChain, "C1\'")
    var aC2 = get_atom_rcn( iResno, iChain, "C2\'")
    var aC3 = get_atom_rcn( iResno, iChain, "C3\'")
    return (angle(aO4, aC1, aC2, aC3) < 0.0)
}

function is_r_res( iResno, iChain) {
    return ({(resno=iResno) and (chain=iChain) and thisModel and purine})
}

function repair_p_res(res, iChain) {
    var aP = get_atom_rcn( res, iChain, "P")
    plico_minimize( {connected(aP) or aP})
}

function who_pairs(iRes, iChain) {
    var aC4or6 =  get_atom_rcn( iRes, iChain, "C4")
    var aN1or3 =  get_atom_rcn( iRes, iChain, "N1")
    if ({aN1or3 and purine}.size = 0) {
        aC4or6 =  get_atom_rcn( iRes, iChain, "C6")
        aN1or3 =  get_atom_rcn( iRes, iChain, "N3")
    }
    if (aN1or3) {
    var near = within(3.2, aN1or3) and {resno!=iRes} and {element="N"}
        for (var i = 1; i <= near.size; i++) {
            var dist = distance(near[i], aN1or3)
            var ang = abs(angle(near[i], aN1or3, aC4or6))
            if (ang > 150) {
                return [near[i].resno, near[i].chain, dist, ang]
            }
        }
    }
    return [-1, aC4or6.chain, -1, -1]
}

function who_almost_pairs(iRes, iChain) {
    var aC2 =  get_atom_rcn( iRes, iChain, "C2")
    var aC4or6 =  get_atom_rcn( iRes, iChain, "C4")
    var aN1or3 =  get_atom_rcn( iRes, iChain, "N1")
    var pname = "C6"
    if ({aN1or3 and purine}.size = 0) {
        aC4or6 =  get_atom_rcn( iRes, iChain, "C6")
        aN1or3 =  get_atom_rcn( iRes, iChain, "N3")
        pname = "C4"
    }
    if (aN1or3) {
        var near = within(3.4, aN1or3) and {resno!=iRes} and {element="N"}
        for (var i = 1; i <= near.size; i++) {
            var aC6or4 = get_atom_rcn(near[i].resno, near[i].chain, pname)
            var aC2p = get_atom_rcn(near[i].resno, near[i].chain, "C2")
            var dist = distance(near[i], aN1or3)
            var puang = abs(angle(near[i], aN1or3, aC4or6))
            var pyang = abs(angle(aC6or4, near[i], aN1or3))
            var dihedral = abs(angle(aC2p, near[i], aN1or3, aC2))
            if ((puang > 110) and (pyang > 110) and (dihedral < 40)) {
                return [near[i].resno, near[i].chain, dist, puang]
            }
        }
    }
    return [-1, aC4or6.chain, -1, -1]
}

function who_stacks(iRes, iChain) {
    var ret = array()
    var aNear = ((within(4.0, {(resno=iRes) and base}) )
        and {base} and {not resno=iRes})
    var done = array()
    for (var i = 1; i <= aNear.size; i++) {
        var jRes = aNear[i].resno
        if (not done.find(jRes)) {
            var jChain = aNear[i].chain
            var as = gen_as(iRes, jRes, iChain, jChain)
            var d = distance({(resno=iRes) and base}, {(resno=jRes) and base})
            var a1 = angle(as[2], as[3], as[4])
            var a2 = angle(as[5], as[4], as[3])
            var dh = angle(as[5], as[4], as[3], as[2])
            var bset = ((connected(as[3]) and not as[2])
                or (connected(as[4]) and not as[5]))
            var a3 = angle(bset[1], bset[2], bset[3])
            var a4 = angle(bset[2], bset[3], bset[4])

            var isStacked = true

            # Bases are parallel as sin(a1) = sin(a2) and sin(a3) = sin(a4)
            if (abs(sin(a1)-sin(a2)) > 20) {
                isStacked = false
            }
            if (abs(sin(a3)-sin(a4)) > 20) {
                isStacked = false
            }

            # Bases are stacked as d*sin(a1) < 6.0 and d3 = 0.0
            if (d*sin(a1) > 6.2) {
                isStacked = false
            }
            if (abs(dh) > 30) {
                #isStacked = false
            }

            if (isStacked) {
                ret += aNear[i].resno
            }
            done += jRes
        }
    }

    return ret
}

function match_nt(mask, nt) {
    var ret = false
    switch (mask) {
    case "A":
    case "U":
    case "C":
    case "G":
        ret = (mask = nt)
        break
    case "N":
        ret = true
        break
    case "M":
        ret = ((nt=="A") or (nt=="C"))
        break
    case "Y":
        ret = ((nt=="U") or (nt=="C"))
        break
    case "R":
        ret = ((nt=="A") or (nt=="G"))
        break
    }
    return ret
}

# Calls function match_nt above
function select_seqs(seq, r5, r3, iChain, f, m) {
    select none
    for (var i = r5; i < r3; i++) {
        var j = 0
        for (; j < seq.size; j++) {
            if ((gNTres[i+j])[2] >= 0) {
                break
            }
            if ((i+j) >= r3) {
                break
            }
            if (not match_nt(seq[j+1], gSeq[i+j])) { # CALL
                break
            }
        }
        if (j == seq.size) {
            print format("%s at %d (%s-%s-%s)%s", seq, i,
                gSeq[i-1],
                gSeq[i][i+j-1],
                gSeq[i+j], aster)
            var rset = {(resno=i) and (chain=iChain) and thisModel}
            rset.selected = true
        }
    }
}

function find_tetras(seq, r5, r3, based) {
    select none
    for (var i = r5; i < r3; i++) {
        var j = 0
        for (; j < seq.size; j++) {
            if ((i+j) >= r3) {
                break
            }
            if (not match_nt(seq[j+1], gSeq[i+j])) { # CALL
                break
            }
        }
        if (j == seq.size) {
            if (based) {
                var ends = gSeq[i-1] + gSeq[i+j]
                var bads = ["AA","GG","AG","GA"]
                if (bads.count(ends) > 0) {
                    continue
                }
            }
            print format("%s at %d (%s-%s-%s)%s", seq, i,
                gSeq[i-1],
                gSeq[i][i+j-1],
                gSeq[i+j], aster)
        }
    }
}

# Calls is_form_a
function select_b_form_nts(iChain) {
    select none
    for (var i = get_resno_min(iChain); i <= get_resno_max(iChain); i++) {
        if (not is_form_a(i, iChain)) { # <== external call
            print format("Res %d is form B", i)
            var rset = {(resno=i) and (chain=iChain) and thisModel}
            rset.selected = true
        }
    }
}

function select_3ward_atom(ar3, ares, iChain) {
    var i = ar3.resno
    var aP = get_atom_rcn( i, iChain, "P")
    switch(ar3.atomName) {
    case "O3\'" :
        select {(resno>i) and (resno<ares) and (chain=iChain)
         and thisModel}
        break
    case "P" :
        select {(resno>=i) and (resno<ares) and (chain=iChain)
            and thisModel}
        break
    case "O5\'" :
    case "C5\'" :
    case "C4\'" :
        select {(resno>=i) and (resno<ares) and (chain=iChain)
            and thisModel and not (connected(aP) or aP)}
        break
    case "C3\'" :
        var aO3 = get_atom_rcn( i, iChain, "O3\'")
        select {((resno>i) and (resno<ares) and (chain=iChain)
            and thisModel) or aO3}
        break
    }
}

function select_5ward_atom(ar5, ares, iChain) {
    var i = ar5.resno
    var aP = get_atom_rcn( i, iChain, "P")
    switch(ar5.atomName) {
    case "O3\'" :
        select {(resno<=i) and (resno>ares) and (chain=iChain)
            and thisModel}
        break
    case "P" :
    case "O5\'" :
    case "C5\'" :
        select {((resno<i) and (resno>ares) and (chain=iChain)
            and thisModel) or (connected(aP) or aP)}
        break
    case "C4\'" :
        var aC5 = get_atom_rcn( i, iChain, "C5\'")
        select {((resno<i) and (resno>ares) and (chain=iChain)
            and thisModel) or (connected(aP) or aP or aC5)}
        break
    }
}

function plot_ab_chi( iChain) {
    select none
    for (var i = get_resno_min(iChain); i <= get_resno_max(iChain); i++) {
        var aO4 = get_atom_rcn(i, iChain, "O4\'")
        var aC1 = get_atom_rcn(i, iChain, "C1\'")
        var isR = {aC1 and purine}
        var a1or9 = (isR ? "N9" : "N1")
        var a6or8 = (isR ? "C8" : "C6")
        var aN  = get_atom_rcn(i, iChain, a1or9)
        var aC  = get_atom_rcn(i, iChain, a6or8)
        var aO3 = get_atom_rcn(i, iChain, "O3\'")
        var aC3 = get_atom_rcn(i, iChain, "C3\'")
        var aC4 = get_atom_rcn(i, iChain, "C4\'")
        var aC5 = get_atom_rcn(i, iChain, "C5\'")
        
        var chi = angle(aO4, aC1, aN, aC)
        aN.vx = chi
        var aorb = angle(aO3, aC3, aC4, aC5)
        aN.vy = aorb
        select ADD aN
    }
    plot properties vx vy resno
    set echo top left
    echo "vx = base chi angle        vy = a ==> b form"
}

function move_p_to_close_o3( aP, aO3, ares, force) {
    var cp = aP.xyz
    var aC3 = {connected(aO3) and (atomname="C3\'")}
    select aP
    set_distance_atoms(aO3, aP, 1.74)
    set_angle_atoms(aC3, aO3, aP, 109.0)
    var pt = aP.xyz
    aP.xyz = cp 
    var rotors = gen_nt_rotors(aP.resno, ares, aP.chain)
    move_atom_nt(aP.atomIndex, pt, aP.resno-1, rotors, force)
    ##fix_p_res(aP.resno, aP.chain, true)
}

function move_o3_to_close_p( aO3, aP, ares, force) {
    var cp = aO3.xyz
    var aO5 = {connected(aP) and (atomname="O5\'")}
    select aO3
    set_distance_atoms(aP, aO3, 1.74)
    set_angle_atoms(aO5, aP, aO3, 109.0)
    var pt = aO3.xyz
    aO3.xyz = cp 
    var rotors = gen_nt_rotors(ares, aO3.resno, aO3.chain)
    move_atom_nt(aO3.atomIndex, pt, aO3.resno+1, rotors, force)
    ##fix_p_res(aP.resno, aP.chain, true)
}

# Select mobile before calling
function pivot_to_close_atoms( aMov, aStat, aPivot, dist) {
    var caxis = cross(aStat.xyz, aMov.xyz) - aPivot.xyz
    var dir = 1
    var d = distance(aMov, aStat)
    if (d > dist) {
        while (d > dist) {
            rotateSelected @aPivot @caxis @dir
            var nd = distance(aMov, aStat)
            if (nd > d) {
                rotateSelected @aPivot @caxis @{-dir}
                if (dir == 1) {
                    rotateSelected @aPivot @caxis @{-dir}
                    dir = -1
                    continue
                }
                else {
                    break
                }
            }
            d = nd
        }
    }
    else {
        while (d < dist) {
            rotateSelected @aPivot @caxis @dir
            var nd = distance(aMov, aStat)
            if (nd < d) {
                rotateSelected @aPivot @caxis @{-dir}
                if (dir == 1) {
                    rotateSelected @aPivot @caxis @{-dir}
                    dir = -1
                    continue
                }
                else {
                    break
                }
            }
            d = nd
        }
    }
}

function print_adjacent_vs( iChain) {
    var tNN = 0
    var taA = 0
    var tbA = 0
    var taD = 0
    var tbD = 0
    var tabD = 0
    var tc = 0.0
    var rmin = get_resno_min(iChain)
    var rmax = get_resno_max(iChain)
    for (var i = rmin; i < rmax; i++) {
        var aC4 = get_atom_rcn(i, iChain, "C4\'")
        var aC1 = get_atom_rcn(i, iChain, "C1\'")
        var isR = {aC1 and purine}
        var a1or9 = (isR ? "N9" : "N1")
        var aN  = get_atom_rcn(i, iChain, a1or9)
        
        var aC4p = get_atom_rcn(i+1, iChain, "C4\'")
        var aC1p = get_atom_rcn(i+1, iChain, "C1\'")
        isR = {aC1p and purine}
        a1or9 = (isR ? "N9" : "N1")
        var aNp  = get_atom_rcn(i+1, iChain, a1or9)
        
        if (aC4 and aC4p) {
            
            var NN = distance(aN, aNp)
            tNN += NN
            var aA = angle(aC1, aN, aNp)
            taA += aA
            var bA = angle(aC1p, aNp, aN)
            tbA += bA
            
            var aD = angle(aC4, aC1, aN, aNp)
            if ((aD < 0) and (taD > 0)) {
                aD += 360
            }
            taD += aD
            var bD = angle(aC4p, aC1p, aNp, aN)
            if ((bD < 0) and (tbD > 0)) {
                bD += 360
            }
            tbD += bD
            var abD = angle(aC1, aN, aNp, aC1p)
            tabD += abD
            tc++
            print format("%6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %d%s %s",
                NN, aA, aD, bA, bD, abD, tc, aC1.group, aC1p.group)
        }
    }
    
    print format("v1=%6.2f v2=%6.2f v3=%6.2f v4=%6.2f v5=%6.2f v6=%6.2f",
        tNN/tc, taA/tc, taD/tc, tbA/tc, tbD/tc, tabD/tc)
}

function measure_adjacent_vs(r5, iChain) {
    var aC4 = get_atom_rcn(r5, iChain, "C4\'")
    var aC1 = get_atom_rcn(r5, iChain, "C1\'")
    var isR = {aC1 and purine}
    var aN1or9  = get_atom_rcn(r5, iChain, (isR ? "N9" : "N1"))
    var aC8or6 =  get_atom_rcn(r5, iChain, (isR ? "C8" : "C6"))
    
    var aC4p = get_atom_rcn(r5+1, iChain, "C4\'")
    var aC1p = get_atom_rcn(r5+1, iChain, "C1\'")
    isR = {aC1p and purine}
    var aN1or9p  = get_atom_rcn(r5+1, iChain, (isR ? "N9" : "N1"))
    var aC8or6p =  get_atom_rcn(r5+1, iChain, (isR ? "C8" : "C6"))

    measure @aN1or9 @aN1or9p
    measure @aC1 @aN1or9 @aN1or9p        
    measure @aC4 @aC1 @aN1or9 @aN1or9p        
    measure @aC1p @aN1or9p @aN1or9        
    measure @aC4p @aC1p @aN1or9p @aN1or9        
    measure @aC1p @aN1or9p @aN1or9 @aC1       
    
    measure @aC4 @aC1 @aN1or9 @aC8or6        
    measure @aC4p @aC1p @aN1or9p @aC8or6p        
}

function print_pair_vs( iChain) {
    var tNN = 0.0
    var taA = 0.0
    var tbA = 0.0
    var taD = 0.0
    var tbD = 0.0
    var tabD = 0.0
    var tc = 0.0
    var taChi = 0.0
    var tbChi = 0.0
    var rmin = get_resno_min(iChain)
    var rmax = get_resno_max(iChain)
    for (var i = rmin; i < rmax; i++) {
        var aC4 = get_atom_rcn(i, iChain, "C4\'")
        if (aC4) {
            var w = who_pairs(i, iChain) # CALL
            var aC1 = get_atom_rcn(i, iChain, "C1\'")
            var isR = {aC1 and purine}
            var a1or9 = (isR ? "N9" : "N1")
            var aN  = get_atom_rcn(i, iChain, a1or9)
            var a6or8 = (isR ? "C8" : "C6")
            var aC  = get_atom_rcn(i, iChain, a6or8)
            var aC4p = get_atom_rcn(w[1], w[2], "C4\'")
            var aC1p = get_atom_rcn(w[1], w[2], "C1\'")
            
            a1or9 = (isR ? "N1" : "N9") # rev
            var aNp  = get_atom_rcn(w[1], w[2], a1or9)
            a6or8 = (isR ? "C6" : "C8") # rev
            var aCp  = get_atom_rcn(w[1], w[2], a6or8)
            
            var NN = distance(aN, aNp)
            tNN += NN
            var aA = angle(aC1, aN, aNp)
            taA += aA
            var bA = angle(aC1p, aNp, aN)
            tbA += bA
            
            var aD = angle(aC4, aC1, aN, aNp)
            if ((aD < 0) and (taD > 0)) {
                aD += 360
            }
            taD += aD
            var bD = angle(aC4p, aC1p, aNp, aN)
            if ((bD < 0) and (tbD > 0)) {
                bD += 360
            }
            tbD += bD
            var abD = angle(aC1, aN, aNp, aC1p)
            tabD += abD
            
            var aChi = angle(aC4, aC1, aN, aC)
            taChi += aChi
            var bChi = angle(aC4p, aC1p, aNp, aCp)
            tbChi += bChi
            
            tc++
            print format(
                "%6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %d%s %d%s",
                NN, aA, aD, bA, bD, abD, aChi, bChi, tc, aC1.group, w[1],
                aC1p.group)
        }
    }
    
    print format(
        "v1=%6.2f v2=%6.2f v3=%6.2f v4=%6.2f v5=%6.2f v6=%6.2f v7=%6.2f v8=%6.2f",
        tNN/tc, taA/tc, taD/tc, tbA/tc, tbD/tc, tabD/tc, taChi/tc, tbChi/tc)
}

function measure_pair_vs(i, j, iChain, jChain) {
    var aC4 = get_atom_rcn(i, iChain, "C4\'")
    var aC4p = get_atom_rcn(j, jChain, "C4\'")
    if (aC4.size and aC4p.size) {
        var aC1 = get_atom_rcn(i, iChain, "C1\'")
        var isR = {aC1 and purine}
        var a1or9 = (isR ? "N9" : "N1")
        var aN  = get_atom_rcn(i, iChain, a1or9)
        var a6or8 = (isR ? "C8" : "C6")
        var aC  = get_atom_rcn(i, iChain, a6or8)
        var aC1p = get_atom_rcn(j, jChain, "C1\'")
        isR = {aC1p and purine}
        a1or9 = (isR ? "N9" : "N1")
        var aNp  = get_atom_rcn(j, jChain, a1or9)
        a6or8 = (isR ? "C8" : "C6")
        var aCp  = get_atom_rcn(j, jChain, a6or8)
    
        measure @aN @aNp
        measure @aC1p @aNp @aN        
        measure @aC4p @aC1p @aNp @aN        
        measure @aC1 @aN @aNp        
        measure @aC4 @aC1 @aN @aNp        
        measure @aC1p @aNp @aN @aC1
        measure @aC4p @aC1p @aNp @aCp
        measure @aC4 @aC1 @aN @aC
    }
    else {
        print "No pair found"
    }       
}

function minimize_for_collision( r, iChain) {
    fix_p_res( r, iChain, true)

    # External
    var cset = (within(kCtolerance, {resno=r}) and not {resno=r} and not
        connected({resno=r}) and {(chain=iChain) and thisModel})
    if (cset) {
        for (var i = 0; i < cset.size; i++) {
            plico_minimize( {(resno=r) or (resno=@{cset[i].resno})})
        }
    }

    # Internal
    cset = (within(kCtolerance, {(resno=r) and (atomName="OP?")}) and not
        {(atomName="OP?") or (atomName="P")}
        and {(chain=iChain) and thisModel})
    if (cset) {
        for (var i = 0; i < cset.size; i++) {
            plico_minimize( {(resno=r) or (resno=@{cset[i].resno})})
        }
    }
}

function eval_pairing( seq, r5, r3, len) {
    var val = 0
    for (var i = 0; i < len; i++) {
        var c5 = seq[r5+i]
        var c3 = seq[r3-i]
        if ((c5 == "A") and (c3 == "U")) {
            val += 2
        }
        else if ((c5 == "U") and (c3 == "A")) {
            val += 2
        }
        else if ((c5 == "G") and (c3 == "C")) {
            val += 3
        }
        else if ((c5 == "C") and (c3 == "G")) {
            val += 3
        }
        else if ((c5 == "G") and (c3 == "U")) {
            val += 2
        }
        else if ((c5 == "U") and (c3 == "G")) {
            val += 1
        }
        else if ((c5 == "A") and (c3 == "A")) {
            val -= 3
        }
        else if ((c5 == "G") and (c3 == "G")) {
            val -= 3
        }
    }
    return val
}
 
function update_atomnos(iChain) {
    print "Update atomnos..."
    var b = {(chain=iChain) and thisModel}.atomno.min
    for (var r = get_resno_min(iChain); r <= get_resno_max(iChain); r++) {
        var rset = {(resno=r) and (chain=iChain) and thisModel}
        for (var n = rset.atomno.min; n <= rset.atomno.max; n++) {
            var a = {(resno=r) and (chain=iChain) and thisModel and (atomno=n)}
            a.atomno = -b
            b++
        }
    }
    for (var i = -b; i < 0; i++) {
        var a = {(chain=iChain) and thisModel
            and (atomno=i)}
        a.atomno *= -1
    }
    print "Updated"    
}

function replane_bases(iChain) {
    for (var i = get_resno_min(iChain); i <= get_resno_max(iChain); i++) {
        plico_minimize( {resno=i and base})
    }
}

function measure_p_dihedrals(i, iChain) {
        var j = i-1
        var as = array()
        as += get_atom_rcn( j, iChain, "C4\'")
        as += get_atom_rcn( j, iChain, "C3\'")
        as += get_atom_rcn( j, iChain, "O3\'")
        as += get_atom_rcn( i, iChain, "P")
        as += get_atom_rcn( i, iChain, "O5\'")
        as += get_atom_rcn( i, iChain, "C5\'")
        as += get_atom_rcn( i, iChain, "C4\'")
        as += get_atom_rcn( i, iChain, "C3\'")
        for (var k = 1; k <= (as.size-3); k++) {
            measure @{as[k+0]} @{as[k+1]} @{as[k+2]} @{as[k+3]}
        }
}

function pair_stem( r5, r3, ar5, iChain) {

    # Pair entire stem
    var c = r3 - r5 + 1
    for (var k = (c\2)-1 ; k >= 0; k--) {
        pair_it_res(r5+k, r3-k, ar5, iChain, iChain) # CALL
    }
}

# From 2LU0
function make_uncg_loop(res5, ares5, iChain) {
    var va = array()
    va[1] = [4.46, 100.5, -83.4, 114.4, 99.8, 39.2, -22.2, -11.8]
    va[2] = [7.89, 93.6, 50.7, 29.7, 75.3, 163.6, 14.2 -22.2]
    va[3] = [7.57, 90.9, -11.9, 59.5, 75.6, 74.1, 33.5, 14.2]
    va[4] = [6.09, 124.9, -5.8, 46.4, -72.6, -27.1, -138.3, 33.5]
    va[5] = [5.1, 64.0, -102.5, 98.2, 87.0, 58.6, -9.2, -138.3]
    for (var i = 1; i <= 5; i++) {
        var as = gen_as(res5+i-2, res5+i-1, iChain, iChain)
    
        var vs = array()
        vs[1] = (va[i])[1]   # distance res5 N9or1 and res3 N9or1      
        vs[2] = (va[i])[2]   # angle res5 N9or1 and res3 N9or1 C1      
        vs[3] = (va[i])[3]   # dihedral res5 N9or1 and res3 N9or1 C1 C4
        vs[4] = (va[i])[4]   # angle res5 C1 N9or1 and res3 N9or1      
        vs[5] = (va[i])[5]   # dihedral res5 C4 N9or1 C1 and res3 N9or1
        vs[6] = (va[i])[6]   # dihedral res5 N9or1 C1 and res3 N9or1 C1
        vs[7] = (va[i])[7]   # dihedral chi res3
        vs[8] = (va[i])[8]   # dihedral chi res5
    
        select {(resno < @{res5+i-1}) and (resno > ares5)
            and (chain=iChain) and thisModel}
        move_it(as, vs)
        fix_p_res(res5+i-1, iChain, true)
    }
}

# From 2LU0
function make_gnra_loop(res5, ares5, iChain) {
    var va = array()
    va[1] = [4.37, 93.9, -84.9, 117.9, 99.1, 42.0, -19.4, 3.2]
    va[2] = [7.94, 77.5, 47.5, 58.9, 87.1, 179.5, 4.5 -19.4]
    va[3] = [4.9, 102.9, -98.6, 79.2, 112.4, 61.2, -29.6, 4.5]
    va[4] = [4.97, 124.5, -70.0, 99.8, 93.2, 34.0, 2.5, -29.6]
    va[5] = [4.31, 76.0, -98.4, 106.6, 97.2, 45.6, -15.5, 2.5]
    for (var i = 1; i <= 5; i++) {
        var as = gen_as(res5+i-2, res5+i-1, iChain, iChain)
    
        var vs = array()
        vs[1] = (va[i])[1]   # distance res5 N9or1 and res3 N9or1      
        vs[2] = (va[i])[2]   # angle res5 N9or1 and res3 N9or1 C1      
        vs[3] = (va[i])[3]   # dihedral res5 N9or1 and res3 N9or1 C1 C4
        vs[4] = (va[i])[4]   # angle res5 C1 N9or1 and res3 N9or1      
        vs[5] = (va[i])[5]   # dihedral res5 C4 N9or1 C1 and res3 N9or1
        vs[6] = (va[i])[6]   # dihedral res5 N9or1 C1 and res3 N9or1 C1
        vs[7] = (va[i])[7]   # dihedral chi res3
        vs[8] = (va[i])[8]   # dihedral chi res5
    
        select {(resno < @{res5+i-1}) and (resno > ares5)
            and (chain=iChain) and thisModel}
        move_it(as, vs)
        fix_p_res(res5+i-1, iChain, true)        
    }
}


# end of plicoNTcommon.spt

Contributors

Remig