Visualizes halogen bonds.
11823
def show_halogen(self):
        """"""
        grp = self.getPseudoBondGroup("HalogenBonds-%i" % self.tid, associateWith=[self.model])
        grp.lineWidth = 3
        for i in self.plcomplex.halogen_bonds:
            b = grp.newPseudoBond(self.atoms[i[0]], self.atoms[i[1]])
            b.color = self.colorbyname('turquoise')

            self.bs_res_ids.append(i.acc_id)
664785
public void placeAliphaticHeavyChain(IAtomContainer molecule, IAtomContainer chain) throws CDKException {
        //logger.debug("******** Place aliphatic Chain *********");
        int[] first = new int[2];
        int counter = 1;
        int nextAtomNr = 0;
        String id1 = "";
        String id2 = "";
        String id3 = "";
        first = findHeavyAtomsInChain(molecule, chain);
        distances = new double[first[1]];
        firstAtoms = new int[first[1]];
        angles = new double[first[1]];
        secondAtoms = new int[first[1]];
        dihedrals = new double[first[1]];
        thirdAtoms = new int[first[1]];
        firstAtoms[0] = first[0];
        molecule.getAtom(firstAtoms[0]).setFlag(CDKConstants.VISITED, true);
        int hybridisation = 0;
        for (int i = 0; i < chain.getAtomCount(); i++) {
            if (isHeavyAtom(chain.getAtom(i))) {
                if (!chain.getAtom(i).getFlag(CDKConstants.VISITED)) {
                    //logger.debug("Counter:" + counter);
                    nextAtomNr = molecule.indexOf(chain.getAtom(i));
                    id2 = molecule.getAtom(firstAtoms[counter - 1]).getAtomTypeName();
                    id1 = molecule.getAtom(nextAtomNr).getAtomTypeName();

                    if (molecule.getBond(molecule.getAtom(firstAtoms[counter - 1]), molecule.getAtom(nextAtomNr)) == null)
                        throw new CDKException("atoms do not form a chain, please use ModelBuilder3D");

                    distances[counter] = getBondLengthValue(id1, id2);
                    //logger.debug(" Distance:" + distances[counter]);
                    firstAtoms[counter] = nextAtomNr;
                    secondAtoms[counter] = firstAtoms[counter - 1];
                    if (counter > 1) {
                        id3 = molecule.getAtom(firstAtoms[counter - 2]).getAtomTypeName();
                        hybridisation = getHybridisationState(molecule.getAtom(firstAtoms[counter - 1]));
                        angles[counter] = getAngleValue(id1, id2, id3);
                        //Check if sp,sp2
                        if (angles[counter] == -1) {
                            if (hybridisation == 3) {
                                angles[counter] = DEFAULT_SP3_ANGLE;
                            } else if (hybridisation == 2) {
                                angles[counter] = DEFAULT_SP2_ANGLE;
                            } else if (hybridisation == 1) {
                                angles[counter] = DEFAULT_SP_ANGLE;
                            }
                        }
                        thirdAtoms[counter] = firstAtoms[counter - 2];
                        //logger.debug(" Angle:" + angles[counter]);
                    } else {
                        angles[counter] = -1;
                        thirdAtoms[counter] = -1;
                    }
                    if (counter > 2) {
                        //double bond
                        try {
                            if (getDoubleBondConfiguration2D(
                                    molecule.getBond(molecule.getAtom(firstAtoms[counter - 1]),
                                            molecule.getAtom(firstAtoms[counter - 2])),
                                    (molecule.getAtom(firstAtoms[counter])).getPoint2d(),
                                    (molecule.getAtom(firstAtoms[counter - 1])).getPoint2d(),
                                    (molecule.getAtom(firstAtoms[counter - 2])).getPoint2d(),
                                    (molecule.getAtom(firstAtoms[counter - 3])).getPoint2d()) == 5) {
                                dihedrals[counter] = DIHEDRAL_BRANCHED_CHAIN;
                            } else {
                                dihedrals[counter] = DIHEDRAL_EXTENDED_CHAIN;
                            }
                        } catch (CDKException ex1) {
                            dihedrals[counter] = DIHEDRAL_EXTENDED_CHAIN;
                        }
                    } else {
                        dihedrals[counter] = -1;
                    }
                    counter++;
                }
            }
        }
    }
1327580
private List<Integer> getBondNosForRingGroup(IAtomContainer molecule, List<Integer> ringGroup, IRingSet ringSet) {
        List<Integer> btc = new ArrayList<Integer>();
        for (Integer i : ringGroup) {
            for (IBond bond : ringSet.getAtomContainer(i).bonds()) {
                if (btc.size() > 0) {
                    if (!btc.contains(molecule.indexOf(bond))) {
                        btc.add(molecule.indexOf(bond));
                    }
                } else {
                    btc.add(molecule.indexOf(bond));
                }
            }
        }
        return btc;
    }
601479
def display_upstream_structure(structure_dict):
    
    graph = _create_graph(structure_dict)
    plt = Image(graph.create_png())
    display(plt)
192047
def zoom_to_ligand(self):
        """"""
        cmd.center(self.ligname)
        cmd.orient(self.ligname)
        cmd.turn('x', 110)  # If the ligand is aligned with the longest axis, aromatic rings are hidden
        if 'AllBSRes' in cmd.get_names("selections"):
            cmd.zoom('%s or AllBSRes' % self.ligname, 3)
        else:
            if self.object_exists(self.ligname):
                cmd.zoom(self.ligname, 3)
        cmd.origin(self.ligname)
811388
def show_sbridges(self):
        """"""
        # Salt Bridges
        grp = self.getPseudoBondGroup("Salt Bridges-%i" % self.tid, associateWith=[self.model])
        grp.lineWidth = 3
        grp.lineType = self.chimera.Dash
        for i, sbridge in enumerate(self.plcomplex.saltbridges):

            m = self.model
            r = m.newResidue("pseudoatoms", " ", 1, " ")
            chargecenter1 = m.newAtom("CHARGE", self.chimera.Element("CHARGE"))
            x, y, z = sbridge.positive_center
            chargecenter1.setCoord(self.chimera.Coord(x, y, z))
            r.addAtom(chargecenter1)

            chargecenter2 = m.newAtom("CHARGE", self.chimera.Element("CHARGE"))
            x, y, z = sbridge.negative_center
            chargecenter2.setCoord(self.chimera.Coord(x, y, z))
            r.addAtom(chargecenter2)

            b = grp.newPseudoBond(chargecenter1, chargecenter2)
            b.color = self.colorbyname('yellow')

            if sbridge.protispos:
                self.bs_res_ids += sbridge.positive_atoms
            else:
                self.bs_res_ids += sbridge.negative_atoms
80235
def find_unpaired_ligand(self):
        
        unpaired_hba, unpaired_hbd, unpaired_hal = [], [], []
        # Unpaired hydrogen bond acceptors/donors in ligand (not used for hydrogen bonds/water, salt bridges/mcomplex)
        involved_atoms = [hbond.a.idx for hbond in self.hbonds_pdon] + [hbond.d.idx for hbond in self.hbonds_ldon]
        [[involved_atoms.append(atom.idx) for atom in sb.negative.atoms] for sb in self.saltbridge_lneg]
        [[involved_atoms.append(atom.idx) for atom in sb.positive.atoms] for sb in self.saltbridge_pneg]
        [involved_atoms.append(wb.a.idx) for wb in self.water_bridges if wb.protisdon]
        [involved_atoms.append(wb.d.idx) for wb in self.water_bridges if not wb.protisdon]
        [involved_atoms.append(mcomplex.target.atom.idx) for mcomplex in self.metal_complexes
         if mcomplex.location == 'ligand']
        for atom in [hba.a for hba in self.ligand.get_hba()]:
            if atom.idx not in involved_atoms:
                unpaired_hba.append(atom)
        for atom in [hbd.d for hbd in self.ligand.get_hbd()]:
            if atom.idx not in involved_atoms:
                unpaired_hbd.append(atom)

        # unpaired halogen bond donors in ligand (not used for the previous + halogen bonds)
        [involved_atoms.append(atom.don.x.idx) for atom in self.halogen_bonds]
        for atom in [haldon.x for haldon in self.ligand.halogenbond_don]:
            if atom.idx not in involved_atoms:
                unpaired_hal.append(atom)

        return unpaired_hba, unpaired_hbd, unpaired_hal
1098914
def pication(rings, pos_charged, protcharged):
    
    data = namedtuple(
        'pication', 'ring charge distance offset type restype resnr reschain restype_l resnr_l reschain_l protcharged')
    pairings = []
    if len(rings) == 0 or len(pos_charged) == 0:
        return pairings
    for ring in rings:
        c = ring.center
        for p in pos_charged:
            d = euclidean3d(c, p.center)
            # Project the center of charge into the ring and measure distance to ring center
            proj = projection(ring.normal, ring.center, p.center)
            offset = euclidean3d(proj, ring.center)
            if not config.MIN_DIST < d < config.PICATION_DIST_MAX or not offset < config.PISTACK_OFFSET_MAX:
                continue
            if type(p).__name__ == 'lcharge' and p.fgroup == 'tertamine':
                # Special case here if the ligand has a tertiary amine, check an additional angle
                # Otherwise, we might have have a pi-cation interaction 'through' the ligand
                n_atoms = [a_neighbor for a_neighbor in OBAtomAtomIter(p.atoms[0].OBAtom)]
                n_atoms_coords = [(a.x(), a.y(), a.z()) for a in n_atoms]
                amine_normal = np.cross(vector(n_atoms_coords[0], n_atoms_coords[1]),
                                        vector(n_atoms_coords[2], n_atoms_coords[0]))
                b = vecangle(ring.normal, amine_normal)
                # Smallest of two angles, depending on direction of normal
                a = min(b, 180 - b if not 180 - b < 0 else b)
                if not a > 30.0:
                    resnr, restype = whichresnumber(ring.atoms[0]), whichrestype(ring.atoms[0])
                    reschain = whichchain(ring.atoms[0])
                    resnr_l, restype_l = whichresnumber(p.orig_atoms[0]), whichrestype(p.orig_atoms[0])
                    reschain_l = whichchain(p.orig_atoms[0])
                    contact = data(ring=ring, charge=p, distance=d, offset=offset, type='regular',
                                   restype=restype, resnr=resnr, reschain=reschain,
                                   restype_l=restype_l, resnr_l=resnr_l, reschain_l=reschain_l,
                                   protcharged=protcharged)
                    pairings.append(contact)
                break
            resnr = whichresnumber(p.atoms[0]) if protcharged else whichresnumber(ring.atoms[0])
            resnr_l = whichresnumber(ring.orig_atoms[0]) if protcharged else whichresnumber(p.orig_atoms[0])
            restype = whichrestype(p.atoms[0]) if protcharged else whichrestype(ring.atoms[0])
            restype_l = whichrestype(ring.orig_atoms[0]) if protcharged else whichrestype(p.orig_atoms[0])
            reschain = whichchain(p.atoms[0]) if protcharged else whichchain(ring.atoms[0])
            reschain_l = whichchain(ring.orig_atoms[0]) if protcharged else whichchain(p.orig_atoms[0])
            contact = data(ring=ring, charge=p, distance=d, offset=offset, type='regular', restype=restype,
                           resnr=resnr, reschain=reschain, restype_l=restype_l, resnr_l=resnr_l,
                           reschain_l=reschain_l, protcharged=protcharged)
            pairings.append(contact)
    return filter_contacts(pairings)
513361
def show_sbridges(self):
        """"""
        for i, saltb in enumerate(self.plcomplex.saltbridges):
            if saltb.protispos:
                for patom in saltb.positive_atoms:
                    cmd.select('PosCharge-P', 'PosCharge-P or (id %i & %s)' % (patom, self.protname))
                for latom in saltb.negative_atoms:
                    cmd.select('NegCharge-L', 'NegCharge-L or (id %i & %s)' % (latom, self.ligname))
                for sbgroup in [['ps-sbl-1-%i' % i, 'Chargecenter-P', saltb.positive_center],
                                ['ps-sbl-2-%i' % i, 'Chargecenter-L', saltb.negative_center]]:
                    cmd.pseudoatom(sbgroup[0], pos=sbgroup[2])
                    cmd.pseudoatom(sbgroup[1], pos=sbgroup[2])
                cmd.distance('Saltbridges', 'ps-sbl-1-%i' % i, 'ps-sbl-2-%i' % i)
            else:
                for patom in saltb.negative_atoms:
                    cmd.select('NegCharge-P', 'NegCharge-P or (id %i & %s)' % (patom, self.protname))
                for latom in saltb.positive_atoms:
                    cmd.select('PosCharge-L', 'PosCharge-L or (id %i & %s)' % (latom, self.ligname))
                for sbgroup in [['ps-sbp-1-%i' % i, 'Chargecenter-P', saltb.negative_center],
                                ['ps-sbp-2-%i' % i, 'Chargecenter-L', saltb.positive_center]]:
                    cmd.pseudoatom(sbgroup[0], pos=sbgroup[2])
                    cmd.pseudoatom(sbgroup[1], pos=sbgroup[2])
                cmd.distance('Saltbridges', 'ps-sbp-1-%i' % i, 'ps-sbp-2-%i' % i)

        if self.object_exists('Saltbridges'):
            cmd.set('dash_color', 'yellow', 'Saltbridges')
            cmd.set('dash_gap', 0.5, 'Saltbridges')
1271946
private static Shape createBondHighlight(IBond bond, double radius) {

        double x1 = bond.getBegin().getPoint2d().x;
        double x2 = bond.getEnd().getPoint2d().x;
        double y1 = bond.getBegin().getPoint2d().y;
        double y2 = bond.getEnd().getPoint2d().y;

        double dx = x2 - x1;
        double dy = y2 - y1;

        double mag = Math.sqrt((dx * dx) + (dy * dy));

        dx /= mag;
        dy /= mag;

        double r2 = radius / 2;

        Shape s = new RoundRectangle2D.Double(x1 - r2, y1 - r2, mag + radius, radius, radius, radius);

        double theta = Math.atan2(dy, dx);

        return AffineTransform.getRotateInstance(theta, x1, y1).createTransformedShape(s);
    }
31140
private void identifyBonds() {
        IAtomContainer spt = getSpanningTree();
        IRing ring;
        int nBasicRings = 0;
        for (int i = 0; i < totalEdgeCount; i++) {
            if (!bondsInTree[i]) {
                ring = getRing(spt, molecule.getBond(i));
                for (int b = 0; b < ring.getBondCount(); b++) {
                    int m = molecule.indexOf(ring.getBond(b));
                    cb[nBasicRings][m] = 1;
                }
                nBasicRings++;
            }
        }
        bondsAcyclicCount = 0;
        bondsCyclicCount = 0;
        for (int i = 0; i < totalEdgeCount; i++) {
            int s = 0;
            for (int j = 0; j < nBasicRings; j++) {
                s += cb[j][i];
            }
            switch (s) {
                case (0): {
                    bondsAcyclicCount++;
                    break;
                }
                case (1): {
                    bondsCyclicCount++;
                    break;
                }
                default: {
                    bondsCyclicCount++;
                }
            }
        }
        identifiedBonds = true;
    }
640456
def assign_valence(mol):
    """"""
    for u, v, bond in mol.bonds_iter():
        if bond.order == 2:
            mol.atom(u).pi = 1
            mol.atom(v).pi = 1
            if mol.atom(u).symbol == "O" and not mol.atom(u).charge:
                mol.atom(v).carbonyl_C = 1
            if mol.atom(v).symbol == "O" and not mol.atom(v).charge:
                mol.atom(u).carbonyl_C = 1
        elif bond.order == 3:
            mol.atom(u).pi = mol.atom(v).pi = 2
    max_nbr = {"C": 4, "Si": 4, "N": 3, "P": 3, "As": 3,
               "O": 2, "S": 2, "Se": 2, "F": 1, "Cl": 1, "Br": 1, "I": 1}
    for i, nbrs in mol.neighbors_iter():
        atom = mol.atom(i)
        if len(nbrs) == 2 and all(bond.order == 2 for bond in nbrs.values()):
            atom.pi = 2  # sp (allene, ketene)
        if atom.symbol in max_nbr:
            h_cnt = max_nbr[atom.symbol] - len(nbrs) - atom.pi + atom.charge
            if h_cnt > 0:
                mol.atom(i).add_hydrogen(h_cnt)
    mol.descriptors.add("Valence")
950758
private  void trackHBondEnergy(int i, int j, double energy) {

		if (groups[i].getPDBName().equals("PRO")) {
			logger.debug("Ignore: PRO {}",groups[i].getResidueNumber());
			return;
		}

		SecStrucState stateOne = getSecStrucState(i);
		SecStrucState stateTwo = getSecStrucState(j);

		double acc1e = stateOne.getAccept1().getEnergy();
		double acc2e = stateOne.getAccept2().getEnergy();

		double don1e = stateTwo.getDonor1().getEnergy();
		double don2e = stateTwo.getDonor2().getEnergy();

		//Acceptor: N-H-->O
		if (energy < acc1e) {
			logger.debug("{} < {}",energy,acc1e);
			stateOne.setAccept2(stateOne.getAccept1());

			HBond bond = new HBond();
			bond.setEnergy(energy);
			bond.setPartner(j);

			stateOne.setAccept1(bond);

		} else if ( energy < acc2e ) {
			logger.debug("{} < {}",energy,acc2e);

			HBond bond = new HBond();
			bond.setEnergy(energy);
			bond.setPartner(j);

			stateOne.setAccept2(bond);
		}


		//The other side of the bond: donor O-->N-H
		if (energy <  don1e) {
			logger.debug("{} < {}",energy,don1e);
			stateTwo.setDonor2(stateTwo.getDonor1());

			HBond bond = new HBond();
			bond.setEnergy(energy);
			bond.setPartner(i);

			stateTwo.setDonor1(bond);

		} else if ( energy < don2e ) {
			logger.debug("{} < {}",energy,don2e);

			HBond bond = new HBond();
			bond.setEnergy(energy);
			bond.setPartner(i);

			stateTwo.setDonor2(bond);
		}
	}
518343
def display_molecule(mol, autozoom=True):
    '''

    '''
    s = System([mol])
    display_system(s, autozoom=True)
125542
private void label(final ITetrahedralChirality element) {

        final IAtom focus = element.getChiralAtom();
        final IAtom[] atoms = element.getLigands();
        final IBond[] bonds = new IBond[4];

        int p = parity(element.getStereo());
        int n = 0;

        // unspecified centre, no need to assign labels
        if (p == 0) return;

        for (int i = 0; i < 4; i++) {
            if (atoms[i].equals(focus)) {
                p *= indexParity(i); // implicit H, adjust parity
            } else {
                bonds[n] = container.getBond(focus, atoms[i]);
                if (bonds[n] == null)
                    throw new IllegalArgumentException("Inconsistent stereo,"
                                                       + " tetrahedral centre"
                                                       + " contained atom not"
                                                       + " stored in molecule");
                atoms[n] = atoms[i];
                n++;
            }
        }

        // sort coordinates and adjust parity (rank gives us the sorted order)
        int[] rank = new int[n];
        for (int i = 0; i < n; i++)
            rank[i] = i;
        p *= sortClockwise(rank, focus, atoms, n);

        // special case when there are three neighbors are acute and an implicit
        // hydrogen is opposite all three neighbors. The central label needs to
        // be inverted, atoms could be laid out like this automatically, consider
        // CC1C[C@H]2CC[C@@H]1C2
        int invert = -1;
        if (n == 3) {
            // find a triangle of non-sequential neighbors (sorted clockwise)
            // which has anti-clockwise winding
            for (int i = 0; i < n; i++) {
                Point2d a = atoms[rank[i]].getPoint2d();
                Point2d b = focus.getPoint2d();
                Point2d c = atoms[rank[(i + 2) % n]].getPoint2d();
                double det = (a.x - c.x) * (b.y - c.y) - (a.y - c.y) * (b.x - c.x);
                if (det > 0) {
                    invert = rank[(i + 1) % n];
                    break;
                }
            }
        }

        // assign all up/down labels to an auxiliary array
        IBond.Stereo[] labels = new IBond.Stereo[n];
        for (int i = 0; i < n; i++) {
            int v = rank[i];

            // 4 neighbors (invert every other one)
            if (n == 4) p *= -1;

            labels[v] = invert == v ? p > 0 ? DOWN : UP : p > 0 ? UP : DOWN;
        }

        // set the label for the highest priority and available bond
        IBond.Stereo firstlabel      = null;
        boolean      assignTwoLabels = assignTwoLabels(bonds, labels);
        for (int v : priority(atomToIndex.get(focus), atoms, n)) {
            IBond bond = bonds[v];
            if (bond.getStereo() != NONE || bond.getOrder() != SINGLE)
                continue;
            // first label
            if (firstlabel == null) {
                bond.setAtoms(new IAtom[]{focus, atoms[v]}); // avoids UP_INVERTED/DOWN_INVERTED
                bond.setStereo(labels[v]);
                firstlabel = labels[v];
                // don't assign a second label when there are only three ligands
                if (!assignTwoLabels)
                    break;
            }
            // second label
            else if (labels[v] != firstlabel) {
                // don't add if it's possibly a stereo-centre
                if (isSp3Carbon(atoms[v], graph[container.indexOf(atoms[v])].length))
                    break;
                bond.setAtoms(new IAtom[]{focus, atoms[v]}); // avoids UP_INVERTED/DOWN_INVERTED
                bond.setStereo(labels[v]);
                break;
            }
        }

        // it should be possible to always assign labels somewhere -> unchecked exception
        if (firstlabel == null)
            throw new IllegalArgumentException("could not assign non-planar (up/down) labels");
    }
1245191
def draw_figure(self,data_for_color=None, data_for_size=None, data_for_clouds=None, rot_bonds=None, color_for_clouds="Blues", color_type_color="viridis"):
        
        self.molecule = Molecule(self.topol_data)

        self.draw = Draw(self.topol_data,self.molecule,self.hbonds,self.pistacking,self.salt_bridges,self.lig_descr)
        self.draw.draw_molecule(data_for_color, data_for_size, data_for_clouds, rot_bonds, color_for_clouds, color_type_color)

        self.figure = Figure(self.molecule,self.topol_data,self.draw)
        self.figure.add_bigger_box()
        self.figure.manage_the_plots()
        self.figure.draw_white_circles()
        self.figure.put_everything_together()
        self.figure.write_final_draw_file(self.output_name)
1155344
def plot_degbandshalffill():
    """"""
    ulim = [3.45, 5.15, 6.85, 8.55]
    bands = range(1, 5)
    for band, u_int in zip(bands, ulim):
        name = 'Z_half_'+str(band)+'band'
        dop = [0.5]
        data = ssplt.calc_z(band, dop, np.arange(0, u_int, 0.1),0., name)
        plt.plot(data['u_int'], data['zeta'][0, :, 0], label='$N={}$'.format(str(band)))

    ssplt.label_saves('Z_half_multiorb.png')