User:Remig/plico/plicoCommonNT

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.4 beta    7/8/2014 -handle multiple frames
#
#   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
kNTcommon = 4
kC5O5PO3B = -71.0
kO5PO3C3B = -107.0
kPO3C3C4B = -161.5
kO3C3C4C5B = 125.5
kC3C4C5O5B = 55.65
kC4C5O5PB = 169.0

kO4C4C3C2B = 15.92
kC4O4C1C2B = -41.7
kC4O4C1NxB = -159.0
kC5C4O4C1B = 146.3
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
kC4O4C1NxA = -117.4
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 =   atom_rcn( cres, iChain, "P")
    var aO5 =  atom_rcn( cres, iChain, "O5\'")
    var aC5 =  atom_rcn( cres, iChain, "C5\'")
    var aC4 =  atom_rcn( cres, iChain, "C4\'")
    var aOP1 = atom_rcn( cres, iChain, "OP1")
    var aOP2 = atom_rcn( cres, iChain, "OP2")
    var aO3p = atom_rcn( pres, iChain, "O3\'")
    var aC3p = atom_rcn( pres, iChain, "C3\'")
    if (aO3p.size > 0) {

        set_distance_atoms(aC5, aO3p, 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)
        minimize select {connected(aP) or aP}
    }
}

function fix_p_res(cres, iChain, force) {
    var f = (_frameID/1000000)
    var m = (_frameID%1000000)
    var pres = cres-1
    var aP   = atom_rcn( cres, iChain, "P")
    var aO5  = atom_rcn( cres, iChain, "O5\'")
    var aC5  = atom_rcn( cres, iChain, "C5\'")
    var aC4  = atom_rcn( cres, iChain, "C4\'")
    var aOP1 = atom_rcn( cres, iChain, "OP1")
    var aOP2 = atom_rcn( cres, iChain, "OP2")
    var aO3p = atom_rcn( pres, iChain, "O3\'")
    var aC3p = atom_rcn( pres, iChain, "C3\'")
    if ((aO3p.size > 0)and (aC4.size > 0)) {

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

        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
        }
        if (force and (abs(distance(aO3p, aC5)-4.11) < kDtolerance)) {
            # Push away
            select {(resno <= @{aC5.resno}) and (chain=iChain)
                and (file=f) and (model=m)}
            set_distance_atoms(aO3p, aC5, 4.11)
        }
        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 {
                    break
                }
            }
            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)
        aOP1.xyz = get_tet_idx(aO5.atomIndex, aP.atomIndex, aO3p.atomIndex, 1.73)
        if (force) {
            minimize select {connected(aP) or aP}
        }
    }
}

function get_rotors_res(res) {
    var f = (_frameID/1000000)
    var m = (_frameID%1000000)
    var rotors = array()
    var sRes = res
    var mRes = sRes-1
    var iChain = {(resno=res) and (atomName="P")
        and (file=f) and (model=m)}.chain
    var mC4 = atom_rcn( mRes, iChain, "C4\'")
    var mC3 = atom_rcn( mRes, iChain, "C3\'")
    var mO3 = atom_rcn( mRes, iChain, "O3\'")
    var sP =  atom_rcn( sRes, iChain, "P"   )
    var sO5 = atom_rcn( sRes, iChain, "O5\'")
    var sC5 = atom_rcn( sRes, iChain, "C5\'")
    var sC4 = atom_rcn( sRes, iChain, "C4\'")
    var sC3 = 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 = atom_rcn( res, iChain, "O4\'")
    var aC1 = atom_rcn( res, iChain, "C1\'")
    var isR  = ((aC1 and {purine}).size > 0)
    var N1or9 = (isR ? "N9" : "N1")
    var C6or8 = (isR ? "C8" : "C6")

    var aN = atom_rcn(res, iChain, N1or9)
    var aC = 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 = atom_rcn(res, iChain, "C5\'")
    var aC4 = atom_rcn(res, iChain, "C4\'")
    var aC3 = atom_rcn(res, iChain, "C3\'")
    var aO3 = 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
}

# Moved object must be selected, fixed object not
# as[6] = fixed[1-3] moved[4-6]
# 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])
}

function gen_as(res5, res3, iChain, jChain) {
    var as = array()
    as[1] = atom_rcn(res3, jChain, "C4\'")
    as[2] = atom_rcn(res3, jChain, "C1\'")
    as[3] = connected(as[2]) and {element="N"}
    as[5] = atom_rcn(res5, iChain, "C1\'")
    as[6] = atom_rcn(res5, iChain, "C4\'")
    as[4] = connected(as[5]) and {element="N"}

    return as
}

# Pair res5 on res3 moving res <= res3
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 isA = is_form_a(res5, iChain)
    var vs = array()
    vs[1] = 9.00    # distance iN-jN (1tna[11-24]=9.00)
    vs[2] = 124.6   # angle iN-jN-jC1 (1tna[11-24]=99.6)//126.0
    vs[3] = (isA ? 160.0 : -140.0)   # dihedral iN-jN-jC1-jC4 (1tna[11-24]=160.0)//-142.6
    vs[4] = 124.6   # angle iC1-iN-jN (1tna[11-24]=124.6)
    vs[5] = (isA ? 160.0 : -140.0)   # dihedral iC4-iN-iC1-jN (1tna[11-24]=160.0)//-138.6
    vs[6] = (isA ? -5.0 : -1.6)    # dihedral iN-iC1-jN-jC1 (1tna[11-24]=-5.0)//-1.6

    if (ares < 0) {
        select ((resno=res5) and (chain=iChain) and (file=f) and (model=m))
    }
    else {
        select ((resno <= ares) and (chain=iChain) and (file=f) and (model=m))
    }
    move_it(as, vs)
    fix_p_res(res5, iChain, TRUE)
    fix_p_res(res3, jChain, 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
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 vs = array()
    vs[1] = 7.00    # distance iN-jN
    vs[2] = 89.1    # angle iN-jN-jC1
    vs[3] = -49.9   # dihedral iN-jN-jC1-jC4
    vs[4] = 83.4    # angle iC1-iN-jN
    vs[5] = 125.7   # dihedral iC4-iN-iC1-jN
    vs[6] = 5.8     # dihedral iN-iC1-jN-jC1

    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(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 f = (_frameID/1000000)
    var m = (_frameID%1000000)
    var as = gen_as(res5, res3, iChain, jChain)

    var vs = array()
    vs[1] = 8.23   # distance iN-jN
    vs[2] = 32.4   # angle iN-jN-jC1
    vs[3] = -26.8  # dihedral iN-jN-jC1-jC4
    vs[4] = 99.6   # angle iC1-iN-jN
    vs[5] = 57.4   # dihedral iC4-iN-iC1-jN
    vs[6] = 179.1  # dihedral iN-iC1-jN-jC1

    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(res3, jChain, TRUE)
}

# Flatstack res3 on res5 moving just res5
function single_flatstack_res3_on_res5(res5, res3, iChain, jChain) {
    var f = (_frameID/1000000)
    var m = (_frameID%1000000)
    var as = gen_as(res5, res3, iChain, jChain)

    vs = array()
    vs[1] = 6.00    #4.6# distance iN-jN
    vs[2] = 90#75.1    # angle iN-jN-jC1
    vs[3] = 90#135.3   # dihedral iN-jN-jC1-jC4
    vs[4] = 90#89.9    # angle iC1-iN-jN
    vs[5] = -90#-47.3   # dihedral iC4-iN-iC1-jN
    vs[6] = 0#1.7     # dihedral iN-iC1-jN-jC1

    select {(resno=res5) and (chain=iChain) and (file=f) and (model=m)}
    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 f = (_frameID/1000000)
    var m = (_frameID%1000000)
    var as = gen_as(res5, res3, iChain, jChain)

    var vs = array()
    vs[1] = 8.9    # distance iN-jN
    vs[2] = 65.3   # angle iN-jN-jC1
    vs[3] = 55.7   # dihedral iN-jN-jC1-jC4
    vs[4] = 61.2   # angle iC1-iN-jN
    vs[5] = -41.2  # dihedral iC4-iN-iC1-jN
    vs[6] = -138.4 # dihedral iN-iC1-jN-jC1

    select {(resno=res5) and (chain=iChain) and (file=f) and (model=m)}
    move_it(as, vs)
    force_p_res(res5, jChain)
    move_it(as, vs)
    fix_p_res(res5, jChain, TRUE)
}

# Pair A res5 on A res3 Hogsteen (N6-N7)2 moving res5 => res3
function pair_it_h_aa(res5, res3, ares, iChain, jChain) {
    var f = (_frameID/1000000)
    var m = (_frameID%1000000)
    var as = 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[4] =  98.2

    # Set dihedral of iC4 from iN and iC1 and jN (1tna=18.2)
    vs[5] = 18.2

    # 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)

    # Mark C1' xyz and move it back to its orginal position
    var pt = as[5].xyz
    as[5].xyz = cp

    # Collect available P rotors
    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 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 vs = array()
    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[3] =  19.9

    # Set angle of iN9or1 from iN1or3 and jN1or3 (1ana=57.3)
    vs[4] =  57.3

    # 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)

    # Mark C1' xyz and move it back to its orginal position
    var pt = as[6].xyz
    as[6].xyz = cp

    # Collect available P rotors
    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[6], 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[6] kDtolerance)

    force_p_res(ares, iChain)
}

# Stack res rMove on res rFixed
function base_stack_res( rMove, rFixed, iChain, jChain, sep , ang) {
    var f = (_frameID/1000000)
    var m = (_frameID%1000000)
    var isA = is_form_a(rMove, iChain)
    var j = rFixed
    var i = rMove
    var as = array()
    var vs = array()
    as[1] = atom_rcn(j, jChain, "O3\'")
    as[2] = atom_rcn(j, jChain, "C5\'")
    var Njx = (((as[1] and {purine}).size > 0) ? "N9" : "N1")
    as[3] = atom_rcn(j, jChain, Njx,    m)
    as[5] = atom_rcn(i, iChain, "C5\'")
    as[6] = atom_rcn(i, iChain, "O3\'")
    var Nix = (((as[5] and {purine}).size > 0) ? "N9" : "N1")
    as[4] = atom_rcn(i, iChain, Nix,    m)
    select {(resno <= i) and (chain=iChain) and (file=f) and (model=m)}

    # Set distance of iNx from jNx (1tna=4.2)
    vs[1] = sep

    # Set angle Njx Nix C5i (1ana=85.7)
    vs[2] =  (isA ? 85.7 : 77.8)

    # Set dihedral Njx Nix C5i O3i (1tna=179.9)
    vs[3] = (isA ? 179.9 : 173.3)

    # Set angle C5j Njx Nxif (1tna=112.9)
    vs[4] = (isA ? 112.9 : 124.2)

    # Set dihedral O3j C5j Njx Nix (1tna= -20)
    vs[5] =  (isA ? -20 : -14.2)

    # Set dihedral of C5j Njx Nix C5i (1tna=20)
    vs[6] = ang

    move_it(as, vs)
    #force_p_res(i, iChain)
    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
function move_atom_nt(targetIdx, targetPt, ares, iRotors) {
    var f = (_frameID/1000000)
    var m = (_frameID%1000000)
    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 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

                        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
            }
            #***************************************************
            rotateSelected {atomIndex=i3} {atomIndex=i2} @{-dt}

            # If collisions
            var res5 = res3-1
            var set3 = {(resno=res3) and (atomName!="P") and (atomName!="OP1")
                and (file=f) and (model=m)}
            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
                while ((abs(dt) > kDtolerance)
                    and ((set5 and within(kCtolerance, set3)).size > 0)) {
                    dt /= 2.0
                    rotateSelected {atomIndex=i3} {atomIndex=i2} @{dt}
                }
                while ((abs(dt) > kDtolerance)
                    and ((set5 and within(kCtolerance, set3)).size > 0)) {
                    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
}

# If ares < 0 then adjust iRes only
function to_ab_nt_res(res, ares, iChain, toA) {
    var f = (_frameID/1000000)
    var m = (_frameID%1000000)
    var aO3 =  atom_rcn( res, iChain, "O3\'")
    var aC3 =  atom_rcn( res, iChain, "C3\'")
    var aC4 =  atom_rcn( res, iChain, "C4\'")
    var aC5 =  atom_rcn( res, iChain, "C5\'")
    var aC1 =  atom_rcn( res, iChain, "C1\'")
    var aC2 =  atom_rcn( res, iChain, "C2\'")
    var aO2 =  atom_rcn( res, iChain, "O2\'")
    var aO4 =  atom_rcn( res, iChain, "O4\'")

    if (ares < 0) {
        select ((resno=res) and (chain=iChain) and (file=f) and (model=m)
            and not aO3 and not aC3 and not aC4)
    }
    else {
        select ((resno <= ares) and (chain=iChain) and (file=f) and (model=m)
            and not aO3 and not aC3 and not aC4)
    }
    set_dihedral_atoms(aO3, aC3, aC4, aC5, (toA ? kO3C3C4C5A : kO3C3C4C5A))

    # Set chi
    var aNx = -1
    var aCx = -1
    var ang = 0.0
    select {(resno=res) and (chain=iChain) and (file=f) and (model=m) and base}
    if ((aC1 and {purine}).size > 0) {
        aNx =  atom_rcn( res, iChain, "N9")
        aCx =  atom_rcn( res, iChain, "C8")
        ang = (toA ? kPuA : kPuB)
    }
    else {
        aNx =  atom_rcn(res, iChain, "N1")
        aCx =  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 (file=f) and (model=m) and base}
    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))
    if (aO2.size > 0) {
        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)
}

function adjust_nts(res5, res3, iChain, toab, a, s) {
    var savemt = useMinimizationThread
    set useMinimizationThread FALSE

    # 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--) {
        if (toab.size > 0) {
            to_ab_nt_res(i, -1, iChain, (toab == "A"))
            if ((w[i])[1] >= 0) {
                to_ab_nt_res((w[i])[1], -1, (w[i])[2], (toab == "A"))
            }
        }
        if (i < res3) {
            base_stack_res(i, i+1, iChain, iChain, s, a)
        }
    }

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

    # Clean up
    for (var i = res3; i >= res5; i--) {
        fix_p_res(i, iChain, TRUE)
        if ((w[i])[1] >= 0) {
            fix_p_res((w[i])[1], (w[i])[2], TRUE)
        }
    }
    set useMinimizationThread savemt
}

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

function is_r_res( iResno, iChain) {
    var f = (_frameID/1000000)
    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) {
    var aP = atom_rcn( res, iChain, "P")
    minimize select {connected(aP) or aP}
}

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

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 "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, iChain) {
    var f = (_frameID/1000000)
    var m = (_frameID%1000000)
    select none
    for (var i = {(chain=iChain) and (file=f) and (model=m)}.resno.min;
        i <= {(chain=iChain) and (file=f) and (model=m)}.resno.max; i++) {
        var j = 1
        for (; j <= seq.size; j++) {
            var nt = {(chain=iChain) and (file=f) and (model=m) and (resno=@{i+j-1})
                and (atomName="C1\'")}.group[1]
            if (not match_nt(seq[j], nt)) { # <== external call
                break
            }
        }
        if (j > seq.size) {
            print format("%s at %d (%s-%s-%s)", seq, i,
                gSeq[i-1],
                gSeq[i][i+seq.size-1],
                gSeq[i+seq.size])
            var rset = {(resno=i) and (chain=iChain) and (file=f) and (model=m)}
            rset.selected = TRUE
        }
    }
}

# Calls is_form_a
function select_b_form_nts(iChain) {
    var f = (_frameID/1000000)
    var m = (_frameID%1000000)
    select none
    for (var i = {(chain=iChain) and (file=f) and (model=m)}.resno.min;
        i <= {(chain=iChain) and (file=f) and (model=m)}.resno.max; 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 (file=f) and (model=m)}
            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])
        }
    }
}

function select_3ward_atom(ar3, ares, iChain) {
    var f = (_frameID/1000000)
    var m = (_frameID%1000000)
    var i = ar3.resno
    var aP = atom_rcn( i, iChain, "P")
    switch(ar3.atomName) {
    case "O3\'" :
        select {(resno>i) and (resno<ares) and (chain=iChain)
         and (file=f) and (model=m)}
        break
    case "P" :
        select {(resno>=i) and (resno<ares) and (chain=iChain)
            and (file=f) and (model=m)}
        break
    case "O5\'" :
    case "C5\'" :
    case "C4\'" :
        select {(resno>=i) and (resno<ares) and (chain=iChain)
            and (file=f) and (model=m) and not (connected(aP) or aP)}
        break
    case "C3\'" :
        var aO3 = atom_rcn( i, iChain, "O3\'")
        select {((resno>i) and (resno<ares) and (chain=iChain)
            and (file=f) and (model=m)) or aO3}
        break
    }
}

function select_5ward_atom(ar5, ares, iChain) {
    var f = (_frameID/1000000)
    var m = (_frameID%1000000)
    var i = ar5.resno
    var aP = atom_rcn( i, iChain, "P")
    switch(ar5.atomName) {
    case "O3\'" :
        select {(resno<=i) and (resno>ares) and (chain=iChain)
            and (file=f) and (model=m)}
        break
    case "P" :
    case "O5\'" :
    case "C5\'" :
        select {((resno<i) and (resno>ares) and (chain=iChain)
            and (file=f) and (model=m)) or (connected(aP) or aP)}
        break
    case "C4\'" :
        var aC5 = atom_rcn( i, iChain, "C5\'")
        select {((resno<i) and (resno>ares) and (chain=iChain)
            and (file=f) and (model=m)) or (connected(aP) or aP or aC5)}
        break
    }
}

# end of plicoNTcommon.spt

Contributors

Remig