id	summary	reporter	owner	description	type	status	priority	milestone	component	version	resolution	keywords	cc	blockedby	blocking	notify_on_close	platform	project
2334	Residues added to N-terminus of a chain do not appear in Model.polymers	Tristan Croll	Eric Pettersen	"If I use the (new today, still a bit messy) code below (`add_amino_acid_residue`) to do a quick-and-dirty addition of a new residue to the N-terminus of an existing chain, then *most* things work - bonds and cartoons draw correctly, and simulations run - but the new residue doesn't appear in `Model.polymers(missing_structure_treatment=Model.PMS_NEVER_CONNECTS)`.

Code below. The call is something like `add_amino_acid_residue(model, 'PHE', next_res=selected_residues(session)[0])`.

While playing around with attempts to debug this, I also noticed that `Model.reorder_residues()` segfaults. Will raise a separate ticket.

{{{
def set_new_atom_style(session, atoms):
    from chimerax.atomic import selected_atoms, Atom
    from chimerax.core.commands import run
    atoms.draw_modes = Atom.STICK_STYLE
    current_sel = selected_atoms(session)
    session.selection.clear()
    atoms.selected = True
    run(session, ""color sel bychain; color sel byhetero"")
    session.selection.clear()
    current_sel.selected = True

def fix_amino_acid_protonation_state(residue):
    if residue.name not in ('GLU', 'ASP'):
        raise TypeError('This method is only applicable to GLU and ASP residues!')
    if residue.name == 'GLU':
        a = residue.find_atom('HE2')
        if a is not None:
            a.delete()
    elif residue.name == 'ASP':
        a = residue.find_atom('HD2')
        if a is not None:
            a.delete()

def new_residue_from_template(model, template, chain_id, center, residue_number=None, insert_code = ' ', b_factor = 50):
    '''
    Create a new residue based on a template, and add it to the model.
    '''
    if residue_number is None:
        if chain_id in m.residues.chain_ids:
            residue_number = suggest_new_residue_number_for_ligand(model, chain_id)
        else:
            residue_number = 0
    import numpy
    from chimerax.atomic import Atom
    t_coords = numpy.array([a.coord for a in template.atoms])
    t_center = t_coords.mean(axis=0)
    t_coords += numpy.array(center) - t_center
    tatom_to_atom = {}
    r = model.new_residue(template.name, chain_id, residue_number, insert=insert_code)
    for i, ta in enumerate(template.atoms):
        a = tatom_to_atom[ta] = model.new_atom(ta.name, ta.element)
        a.coord = t_coords[i]
        a.bfactor = b_factor
        r.add_atom(a)
        for tn in ta.neighbors:
            n = tatom_to_atom.get(tn, None)
            if n is not None:
                model.new_bond(a, n)
    set_new_atom_style(model.session, r.atoms)
    return r

def add_amino_acid_residue(model, resname, prev_res=None, next_res=None,
        chain_id=None, number=None, center=None, insertion_code=' ', b_factor=50,
        occupancy=1, phi=-135, psi=135):
    session = model.session
    if (not chain_id or not number or not center) and (not prev_res and not next_res):
        raise TypeError('If no anchor residues are specified, chain ID, '
            'number and center must be provided!')
    if prev_res and next_res:
        raise TypeError('Cannot specify both previous and next residues!')
    if prev_res:
        for n in prev_res.find_atom('C').neighbors:
            if n.residue != prev_res:
                raise TypeError('This residue already has another bonded to its '
                    'C terminus!')
        chain_id = prev_res.chain_id
        oxt = prev_res.find_atom('OXT')
        if oxt is not None:
            oxt.delete()
    elif next_res:
        for n in next_res.find_atom('N').neighbors:
            if n.residue != next_res:
                raise TypeError('This residue already has another bonded to its '
                    'N terminus!')
        chain_id = next_res.chain_id
        for hname in ('H2', 'H3'):
            h = next_res.find_atom(hname)
            if h is not None:
                h.delete()
            h = next_res.find_atom('H1')
            if h is not None:
                h.name='H'
            if next_res.name == 'PRO':
                h = next_res.find_atom('H')
                if h:
                    h.delete()
    if number is None:
        if prev_res:
            number = prev_res.number + 1
        elif next_res:
            number = next_res.number - 1

    from chimerax.atomic import mmcif
    tmpl = mmcif.find_template_residue(session, resname)
    from .place_ligand import new_residue_from_template
    import numpy
    # delete extraneous atoms
    r = new_residue_from_template(model, tmpl, chain_id, [0,0,0], number,
            insert_code=insertion_code, b_factor=b_factor)
    r.atoms[numpy.in1d(r.atoms.names, ['OXT', 'HXT', 'H2'])].delete()

    # Translate and rotate residue to (roughly) match the desired position
    if not next_res and not prev_res:
        r.atoms.coords += numpy.array(center) - r.atoms.coords.mean(axis=0)
    else:
        from chimerax.core.geometry import align_points
        if prev_res:
            model.new_bond(r.find_atom('N'), prev_res.find_atom('C'))
            n_pos = _find_next_N_position(prev_res)
            ca_pos = _find_next_CA_position(n_pos, prev_res)
            c_pos = _find_next_C_position(ca_pos, n_pos, prev_res, phi)
            target_coords = numpy.array([n_pos, ca_pos, c_pos])
            align_coords = numpy.array([r.find_atom(a).coord for a in ['N', 'CA', 'C']])
        elif next_res:
            model.new_bond(r.find_atom('C'), next_res.find_atom('N'))
            c_pos = _find_prev_C_position(next_res, psi)
            ca_pos = _find_prev_CA_position(c_pos, next_res)
            o_pos = _find_prev_O_position(c_pos, next_res)
            target_coords = numpy.array([c_pos, ca_pos, o_pos])
            align_coords = numpy.array([r.find_atom(a).coord for a in ['C', 'CA', 'O']])

        tf = align_points(align_coords, target_coords)[0]
        r.atoms.coords = tf*r.atoms.coords
    if r.name in ('GLU', 'ASP'):
        fix_amino_acid_protonation_state(r)
    if r.name == 'PRO':
        r.atoms[r.atoms.names=='H'].delete()
    r.atoms.bfactors = b_factor
    r.atoms.occupancies = occupancy

    set_new_atom_style(session, r.atoms)
    return r



def _find_next_N_position(prev_res):
    from chimerax.atomic.struct_edit import find_pt
    bond_length = 1.34
    angle = 120
    dihedral = 180
    a1 = prev_res.find_atom('C')
    a2 = prev_res.find_atom('CA')
    a3 = prev_res.find_atom('O')
    return find_pt(*[a.coord for a in [a1, a2, a3]], bond_length, angle, dihedral)

def _find_next_CA_position(n_pos, prev_res):
    from chimerax.atomic.struct_edit import find_pt
    bond_length = 1.48
    angle = 124
    omega = 180
    c = prev_res.find_atom('C')
    ca = prev_res.find_atom('CA')
    return find_pt(n_pos, *[a.coord for a in [c, ca]], bond_length, angle, omega)

def _find_next_C_position(ca_pos, n_pos, prev_res, phi):
    from chimerax.atomic.struct_edit import find_pt
    bond_length = 1.53
    angle = 120
    c = prev_res.find_atom('C')
    return find_pt(ca_pos, n_pos, c.coord, bond_length, angle, phi)


def _find_prev_C_position(next_res, psi):
    from chimerax.atomic.struct_edit import find_pt
    bond_length = 1.34
    angle = 120
    a1 = next_res.find_atom('N')
    a2 = next_res.find_atom('CA')
    a3 = next_res.find_atom('C')
    return find_pt(*[a.coord for a in [a1, a2, a3]], bond_length, angle, psi)

def _find_prev_CA_position(c_pos, next_res):
    from chimerax.atomic.struct_edit import find_pt
    bond_length = 1.53
    angle = 120
    omega = 180
    n = next_res.find_atom('N')
    ca = next_res.find_atom('CA')
    return find_pt(c_pos, *[a.coord for a in [n, ca]], bond_length, angle, omega)

def _find_prev_O_position(c_pos, next_res):
    from chimerax.atomic.struct_edit import find_pt
    bond_length = 1.22
    angle = 120
    dihedral = 0
    n = next_res.find_atom('N')
    ca = next_res.find_atom('CA')
    return find_pt(c_pos, *[a.coord for a in (n, ca)], bond_length, angle, dihedral)
}}}"	defect	closed	major		Structure Editing		fixed						all	ChimeraX