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.1 beta    5/16/2014 -lc all functions
#
#   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 = 2
gChain1 = "A"
gChain2 = ""

function set_p_angle_res(cres, ang, iChain) {
print format("setPangleaRes(%d, %5.2f, %s)", cres, ang, iChain)
    var pres = cres-1
    var stator = {(resno=cres) and (atomName="C1\'") and (chain=iChain)}
    var pivot = {(resno=cres) and (atomName="P") and (chain=iChain)}
    var rotor = {(resno=pres) and (atomName="C1\'") and (chain=iChain)}
    select {(resno < cres) and (chain=iChain)}
    set_angle_atoms(stator, pivot, rotor, ang)
}
                
function adjust_p_angle_res(cres, iChain) {
    var pres = cres-1
    var aP = {(resno=cres) and (atomName="P") and (chain=iChain)}
    var aO5 = {(resno=cres) and (atomName="O5\'") and (chain=iChain)}
    var aC5 = {(resno=cres) and (atomName="C5\'") and (chain=iChain)}
    var aC4 = {(resno=cres) and (atomName="C4\'") and (chain=iChain)}
    var aOP1 = {(resno=cres) and (atomName="OP1") and (chain=iChain)}
    var aOP2 = {(resno=cres) and (atomName="OP2") and (chain=iChain)}
    var aO3p = {(resno=pres) and (atomName="O3\'") and (chain=iChain)}
    
    select aO5
    var dist = distance(aO3p, aO5)
    var widen = (dist < 2.84)
    var dir = (widen ? -1 : 1)
    var first = TRUE
    while (abs(dist-2.84) > kDtolerance) {
        rotateSelected @aC4 @aC5 @dir
        var newdist = distance(aO3p, aO5)
        if (widen ? (newdist < dist) : (newdist > dist)) {
            if (first) {
                dir = -dir
                rotateSelected @aC5 @aO5 @dir
            }
            else {
                break
            }
        }
        dist=newdist
        first = FALSE
    }
    aOP2.xyz = get_tet_idx(aO3p.atomIndex, aP.atomIndex, aO5.atomIndex, 1.5)
    aOP1.xyz = get_tet_idx(aO5.atomIndex, aP.atomIndex, aO3p.atomIndex, 1.5)
    minimize select {connected(aP) or aP}
}
                
# Work around - Functions returning values must be in lower case for 14.0.13
function get_interbase_rotors(aP) {
    var rotors = array()
    var sRes = aP.resno
    var mRes = sRes-1
    var iChain = aP.chain
    var mC4 = {(resno=mRes) and (chain=iChain) and (atomName="C4\'")}
    var mC3 = {(resno=mRes) and (chain=iChain) and (atomName="C3\'")}
    var mO3 = {(resno=mRes) and (chain=iChain) and (atomName="O3\'")}
    var sP = {(resno=sRes) and (chain=iChain) and (atomName="P")}
    var sO5 = {(resno=sRes) and (chain=iChain) and (atomName="O5\'")}
    var sC5 = {(resno=sRes) and (chain=iChain) and (atomName="C5\'")}
    var sC4 = {(resno=sRes) and (chain=iChain) and (atomName="C4\'")}
    var sC3 = {(resno=sRes) and (chain=iChain) and (atomName="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
}

# Work around - Functions returning values must be in lower case for 14.0.13
function gen_nt_rotos(res5, res3, iChain) {
    var rotors = array()
    for (var i = res5; i < res3; i++) {
        var aP = {(resno=i) and (chain=iChain) and (atomName="P")}
        rotors += get_interbase_rotors(aP)        
    }
    return rotors
}

function set_res_distance(stator, mover, dist, rotors) {
    var selsave = {selected}
    var cp = mover.xyz
    select mover
    set_distance_atoms(stator, mover, dist)
    var pt = mover.xyz
    mover.xyz = cp
    select selsave
    toab_track_idx(mover.atomIndex, pt, rotors)
    toab_track_idx(mover.atomIndex, pt, rotors)
    toab_track_idx(mover.atomIndex, pt, rotors)
}

# Pair res i on res j
function pair_it(i, j, iChain, jChain) {
    # Calc point offset from res j's C1'
    var jC1 = {(resno=j) and (atomName="C1\'") and (chain=jChain)}
    var jN = {connected(jC1) and (element="N")}
    var jjC1 = {(resno=@{j+1}) and (atomName="C1\'") and (chain=jChain)}
    var jC4 = {(resno=j) and (atomName="C4\'") and (chain=jChain)}
    var iC1 = {(resno=i) and (atomName="C1\'") and (chain=iChain)}
    var iN = {connected(iC1) and (element="N")}
    var iC4 = {(resno=i) and (atomName="C4\'") and (chain=iChain)}
    select {(resno <= i)}

    # First point of a 3-point landing iN
    # Set distance of iN from jN (1ana=9.00)
    set_distance_atoms(jN, iN, 9.00)

    # Set angle of iN from jN and jC1 (1ana=124.6)
    set_angle_atoms(jC1, jN, iN, 124.6)

    # Set dihedral of iN from jN and jC1 and jC4 (1ana=160.0)
    set_dihedral_atoms(jC4, jC1, jN, iN, 160.0)

    # Second point of a 3-point landing iC1
    # Set angle of iC1 from iN nad jN (1ana=124.6)
    set_angle_atoms(jN, iN, iC1, 124.6)

    # Set dihedral of iN from iC1 and jN and jC1 (1ana=-5.0)
    set_dihedral_atoms(jC1, jN, iN, iC1, -5.0)

    # Third point of a 3-point landing iC1
    # Set dihedral of iC4 from iN and iC1 and jN (1ana=160.0)
    set_dihedral_atoms(jN, iN, iC1, iC4, 160.0)

}

function rotate_selected_cd_atoms(a1, a2, ang) {
print format("rotate_selected_cd_atoms(a1=%s, a2=%s, ang=%5.2f", a1, a2, ang)
    #rotateSelected a1 a2 @ang#TBD
    rang = abs(ang)
    iang = 1.0
    dir = ((ang < 0) ? -iang : iang)
    while (rang > 0) {
        rotateSelected @a1 @a2 @dir
        if (is_collision_in_select()) {
        #ca = count_collision_in_select(TRUE)
        #if (ca.size > 0) {
            rotateSelected @a1 @a2 @{-dir}
            break
        }
        rang -= iang
    }
}

# Rotate rotor set to move target atom to its proper place
function toab_track_idx(targetIdx, targetPt, iRotors) {
    var pt = targetPt
    var rotors = iRotors
    var targetNo = {atomIndex=targetIdx}.atomno
    var iChain = {atomIndex=targetIdx}.chain
    gOK = FALSE
    var dist = distance(pt, {atomIndex=targetIdx}.xyz)

    # 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) {
                var i2 = rotors[ri+1]
                var i3 = rotors[ri+2]
                var i4 = rotors[ri+3]
                if ((i2 != targetIdx) and (i3 != targetIdx) and (i4 != targetIdx)) {
                    if ({blocked and {atomIndex=i2}}.count == 0) {
                        var 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
            }
            var i1 = rotors[imax+0]
            var i2 = rotors[imax+1]
            var i3 = rotors[imax+2]
            var i4 = rotors[imax+3]

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

            # Rotate to minimize vector ====================
            select (((atomno >= targetNo) and (chain = iChain) and
                (atomno <= @{{atomIndex=i3}.atomno}))
                or connected({atomIndex=i2}))
            #rotateSelected {atomIndex=i2} {atomIndex=i3} @dt
            rotate_selected_cd_atoms({atomIndex=i2}, {atomIndex=i3}, dt)

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

            # Block rotor
            blocked |= {atomIndex=i2}

        }   # endfor num rotors passes

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

# Work around - Functions returning values must be in lower case for 14.0.13
function is_form_a( iResno, iChain) {
    var aO4 = {(resno=iResno) and (chain=iChain) and (atomName="O4\'")}
    var aC1 = {(resno=iResno) and (chain=iChain) and (atomName="C1\'")}
    var aC2 = {(resno=iResno) and (chain=iChain) and (atomName="C2\'")}
    var aC3 = {(resno=iResno) and (chain=iChain) and (atomName="C3\'")}
    return (angle(aO4, aC1, aC2, aC3) < 0.0)
}

function repair_p_res(res, iChain) {
    var aP = {((resno=res) and (chain=iChain) and (atomName="P"))}
    minimize select {connected(aP) or aP}
}

function pivot_180(res5, res3, iChain) {
    var aO5 = {(resno=res5) and (atomName="O5\'") and (chain=iChain)}
    var bO3 = {(resno=@{res3-1}) and (atomName="O3\'") and (chain=iChain)}
    var aP = {((resno=res5) and (chain=iChain) and (atomName="P"))}
    select {(resno>=res5) and (resno<res3) and not (connected(aP) or aP)}
    fix not selected
    rotateSelected @bO3 @aO5 180.0
    fix none
    repair_p_res(res5, iChain)
}
# end of plicoNTcommon.spt

Contributors

Remig