Extract octahedra as lists of sequence numbers of corner atoms
171056
def get_octahedra(self, atoms, periodicity=3):
        '''
        
        '''
        octahedra = []
        for n, i in enumerate(atoms):
            found = []
            if i.symbol in Perovskite_Structure.B:
                for m, j in enumerate(self.virtual_atoms):
                    if j.symbol in Perovskite_Structure.C and self.virtual_atoms.get_distance(n, m) <= self.OCTAHEDRON_BOND_LENGTH_LIMIT:
                        found.append(m)

            if (periodicity == 3 and len(found) == 6) or (periodicity == 2 and len(found) in [5, 6]):
                octahedra.append([n, found])

        if not len(octahedra): raise ModuleError("Cannot extract valid octahedra: not enough corner atoms found!")
        return octahedra
566874
def wc_lha2dict(lha):
    """"""
    C = OrderedDict()
    # try to read all WCs with 0, 2, or 4 fermions; if not found, set to zero
    for k, (block, i) in WC_dict_0f.items():
        try:
            C[k] = dict(lha['BLOCK'][block]['values'])[i]
        except KeyError:
            C[k] = 0
    for k in definitions.WC_keys_2f:
        try:
            C[k] = lha2matrix(lha['BLOCK']['WC' + k.upper()]['values'], (3,3)).real
        except KeyError:
            C[k] = np.zeros((3,3))
        try: # try to add imaginary part
            C[k] = C[k] + 1j*lha2matrix(lha['BLOCK']['IMWC' + k.upper()]['values'], (3,3))
        except KeyError:
            pass
    for k in definitions.WC_keys_4f:
        try:
            C[k] = lha2matrix(lha['BLOCK']['WC' + k.upper()]['values'], (3,3,3,3))
        except KeyError:
            C[k] = np.zeros((3,3,3,3))
        try: # try to add imaginary part
            C[k] = C[k] + 1j*lha2matrix(lha['BLOCK']['IMWC' + k.upper()]['values'], (3,3,3,3))
        except KeyError:
            pass
    return C
639116
List<IStereoElement> recognise(Set<Projection> projections) {
        
        if (!projections.contains(Projection.Fischer))
            return Collections.emptyList();

        // build atom index and only recognize 2D depictions
        Map<IAtom,Integer> atomToIndex = new HashMap<IAtom, Integer>();
        for (IAtom atom : container.atoms()) {
            if (atom.getPoint2d() == null)
                return Collections.emptyList();
            atomToIndex.put(atom, atomToIndex.size());
        }
        
        RingSearch ringSearch = new RingSearch(container, graph);
        
        final List<IStereoElement> elements = new ArrayList<IStereoElement>(5);

        for (int v = 0; v < container.getAtomCount(); v++) {

            IAtom    focus = container.getAtom(v);
            Elements elem  = Elements.ofNumber(focus.getAtomicNumber());

            if (elem != Carbon)
                continue;
            if (ringSearch.cyclic(v))
                continue;
            if (stereocenters.elementType(v) != Tetracoordinate)
                continue;
            if (!stereocenters.isStereocenter(v))
                continue;

            ITetrahedralChirality element = newTetrahedralCenter(focus,
                                                                 neighbors(v, graph, bonds));

            if (element == null)
                continue;

            // east/west bonds must be to terminal atoms
            IAtom east = element.getLigands()[EAST];
            IAtom west = element.getLigands()[WEST];
            
            if (!east.equals(focus) && !isTerminal(east, atomToIndex))
                continue;
            if (!west.equals(focus) && !isTerminal(west, atomToIndex))
                continue;
            
            elements.add(element);
        }

        return elements;
    }
51727
private static Atom[] getAtoms(Atom[] ca, int[] positions, int length,
			boolean clone) {

		List<Atom> atoms = new ArrayList<Atom>();
		for (int i = 0; i < length; i++) {
			int p = positions[i];
			Atom a;
			if (clone) {
				a = (Atom) ca[p].clone();
				a.setGroup((Group) ca[p].getGroup().clone());
			} else {
				a = ca[p];
			}
			atoms.add(a);
		}
		return atoms.toArray(new Atom[atoms.size()]);
	}
599145
def get_sequences_from_cluster(c1, c2, data):
    """"""
    seqs1 = data[c1]['seqs']
    seqs2 = data[c2]['seqs']
    seqs = list(set(seqs1 + seqs2))
    names = []
    for s in seqs:
        if s in seqs1 and s in seqs2:
            names.append("both")
        elif s in seqs1:
            names.append(c1)
        else:
            names.append(c2)
    return seqs, names
454110
def process_parsed_coords(coords):
    
    geometry = np.zeros(shape=(len(coords), 3), dtype=float)
    for ii, entry in enumerate(coords):
        for jj in range(3):
            geometry[ii, jj] = float(entry[jj])
    return geometry
95162
def shapesides(inputtocheck, inputtype='shape'):
    

    # Define the array of sides to a shape
    shapestosides = {
        'triangle': 3,
        'square': 4,
        'pentagon': 5,
        'hexagon': 6,
        'heptagon': 7,
        'octagon': 8,
        'nonagon': 9,
        'decagon': 10,
        'hendecagon': 11,
        'dodecagon': 12,
        'triskaidecagon': 13,
        'tetrakaidecagon': 14,
        'pentadecagon': 15,
        'hexakaidecagon': 16,
        'heptadecagon': 17,
        'octakaidecagon': 18,
        'enneadecagon': 19,
        'icosagon': 20,
        'triacontagon': 30,
        'tetracontagon': 40,
        'pentacontagon': 50,
        'hexacontagon': 60,
        'heptacontagon': 70,
        'octacontagon': 80,
        'enneacontagon': 90,
        'hectagon': 100,
        'chiliagon': 1000,
        'myriagon': 10000,
        'megagon': 1000000,
        'googolgon': pow(10, 100),
        'ngon': 'n'
    }

    # Define an array with the flipped version of the sides to a shape
    sidestoshapes = dictflip(shapestosides)

    # If the lowercase version of the input type is 'shape'
    if inputtype.lower() == 'shape':
        # If the lowercase version of the shape is in the array
        if inputtocheck.lower() in shapestosides:
            # Return the corresponding sides
            return shapestosides[inputtocheck.lower()]

        # Return 'n'
        return shapestosides['n']

    if inputtype.lower() == 'sides':
        # If the lowercase version of the shape is in the array
        if inputtocheck.lower() in sidestoshapes:
            # Return the corresponding sides
            return sidestoshapes[inputtocheck.lower()]

        # Return 'ngon'
        return sidestoshapes['ngon']

    # Raise a warning
    raise ValueError("Invalid input type.")
616386
def from_string(contents):
        
        lines = contents.split("\n")
        num_sites = int(lines[0])
        coords = []
        sp = []
        prop = []
        coord_patt = re.compile(
            r"(\w+)\s+([0-9\-\.]+)\s+([0-9\-\.]+)\s+([0-9\-\.]+)\s+" +
            r"([0-9\-\.]+)"
        )
        for i in range(2, 2 + num_sites):
            m = coord_patt.search(lines[i])
            if m:
                sp.append(m.group(1))  # this is 1-indexed
                # coords.append(map(float, m.groups()[1:4]))  # this is 0-indexed
                coords.append([float(j)
                               for j in [m.group(i) for i in [3, 4, 2]]])
                prop.append(float(m.group(5)))
        return ZeoVoronoiXYZ(
            Molecule(sp, coords, site_properties={'voronoi_radius': prop})
        )
475273
def _get_ATOM_sequences(self):
        ''''''

        # Get a list of all residues with ATOM or HETATM records
        atom_sequences = {}
        structural_residue_IDs_set = set() # use a set for a quicker lookup
        ignore_HETATMs = True # todo: fix this if we need to deal with HETATMs

        residue_lines_by_chain = []
        structural_residue_IDs_set = []

        present_chain_ids = {}
        for l in self.structure_lines:
            if len(l) > 21 and l[:3] != 'TER':
                present_chain_ids[l[21]] = present_chain_ids.get(l[21], set())
                present_chain_ids[l[21]].add(l[:6])

        model_index = 0
        residue_lines_by_chain.append([])
        structural_residue_IDs_set.append(set())
        full_code_map = {}
        hetatm_map = {}
        full_atom_map = {}
        for l in self.structure_lines:
            chain_id = None
            if l.startswith("TER   "):
                model_index += 1
                residue_lines_by_chain.append([])
                structural_residue_IDs_set.append(set())
            else:
                residue_id = l[21:27]
                if residue_id not in structural_residue_IDs_set[model_index]:
                    residue_lines_by_chain[model_index].append(l)
                    structural_residue_IDs_set[model_index].add(residue_id)
                if l.startswith('ATOM'):
                    chain_id = l[21]
                    # Only use ATOM records to build the code map as chains can have ligands HETATMs
                    full_code_map[chain_id] = full_code_map.get(chain_id, set())
                    full_code_map[chain_id].add(l[17:20].strip())

                    # Only use ATOM records to build the atom map as CA-only chains can have ligands described in full as HETATMs
                    full_atom_map[chain_id] = full_atom_map.get(chain_id, set())
                    full_atom_map[chain_id].add(l[12:16].strip())
                elif l.startswith('HETATM'):
                    chain_id = l[21]
                    hetatm_map[chain_id] = hetatm_map.get(chain_id, set())
                    hetatm_map[chain_id].add(l[17:20].strip())

        # Get the residues used by the residue lines. These can be used to determine the chain type if the header is missing.
        for chain_id in self.atom_chain_order:
            if full_code_map.get(chain_id):
                # The chains may contain other molecules e.g. MG or HOH so before we decide their type based on residue types alone,
                # we subtract out those non-canonicals
                canonical_molecules = full_code_map[chain_id].intersection(dna_nucleotides.union(rna_nucleotides).union(residue_types_3))

                determined_chain_type = None
                if canonical_molecules.union(dna_nucleotides) == dna_nucleotides:
                    determined_chain_type = 'DNA'
                elif canonical_molecules.union(rna_nucleotides) == rna_nucleotides:
                    determined_chain_type = 'RNA'
                elif len(full_code_map[chain_id]) == 1 and 'UNK' in full_code_map[chain_id]:
                    determined_chain_type = 'Unknown'
                elif canonical_molecules:
                    if len(full_atom_map[chain_id]) == 1 and 'CA' in full_atom_map[chain_id]:
                        determined_chain_type = 'Protein skeleton'
                    else:
                        determined_chain_type = 'Protein'
                else:
                    determined_chain_type = PDB._determine_heterogen_chain_type(canonical_molecules)

                if self.chain_types.get(chain_id):
                    assert(self.chain_types[chain_id] == determined_chain_type)
                else:
                    self.chain_types[chain_id] = determined_chain_type

        line_types_by_chain = []
        chain_ids = []
        for model_index in range(len(residue_lines_by_chain)):
            line_types = set()
            if residue_lines_by_chain[model_index]:
                if missing_chain_ids.get(self.pdb_id):
                    chain_ids.append(missing_chain_ids[self.pdb_id])
                else:
                    chain_ids.append(residue_lines_by_chain[model_index][0][21])
            for l in residue_lines_by_chain[model_index]:
                line_types.add(l[0:6])
            if line_types == set(['ATOM']):
                line_types_by_chain.append('ATOM')
            elif line_types == set(['HETATM']):
                line_types_by_chain.append('HETATM')
            else:
                line_types_by_chain.append('Mixed')

        for x in range(0, len(residue_lines_by_chain)):
            residue_lines = residue_lines_by_chain[x]
            line_types = line_types_by_chain[x]
            if ignore_HETATMs and line_types == 'HETATM':
                continue

            for y in range(len(residue_lines)):
                l = residue_lines[y]
                residue_type = l[17:20].strip()
                if l.startswith("HETATM"):
                    if self.modified_residue_mapping_3.get(residue_type):
                        residue_type = self.modified_residue_mapping_3[residue_type]
                    elif y == (len(residue_lines) - 1):
                        # last residue in the chain
                        if residue_type == 'NH2':
                            residue_type = 'UNK' # fixes a few cases e.g. 1MBG, 1K9Q, 1KA6
                        elif ignore_HETATMs:
                            continue

                    elif ignore_HETATMs:
                        continue

                residue_id = l[21:27]
                chain_id = l[21]
                if missing_chain_ids.get(self.pdb_id):
                    chain_id = missing_chain_ids[self.pdb_id]

                if chain_id in self.chain_types:
                    # This means the pdb had SEQRES and we constructed atom_sequences
                    chain_type = self.chain_types[chain_id]
                else:
                    # Otherwise assume this is protein
                    chain_type = 'Protein'

                atom_sequences[chain_id] = atom_sequences.get(chain_id, Sequence(chain_type))

                residue_type = self.modified_residue_mapping_3.get(residue_type, residue_type)

                short_residue_type = None
                if residue_type == 'UNK':
                    short_residue_type = 'X'
                elif chain_type == 'Unknown':
                    assert(False) # we should not reach here - Unknown chains should only contain UNK records
                elif chain_type == 'Protein' or chain_type == 'Protein skeleton':
                    short_residue_type = residue_type_3to1_map.get(residue_type) or protonated_residue_type_3to1_map.get(residue_type) or non_canonical_amino_acids.get(residue_type)
                elif chain_type == 'DNA':
                    short_residue_type = dna_nucleotides_2to1_map.get(residue_type) or non_canonical_dna.get(residue_type)
                elif chain_type == 'RNA':
                    short_residue_type = non_canonical_rna.get(residue_type) or residue_type

                if not short_residue_type:
                    if l.startswith("ATOM") and l[12:16] == ' OH2' and l[17:20] == 'TIP':
                        continue
                    elif not self.strict:
                        short_residue_type = 'X'
                    else:
                        raise NonCanonicalResidueException("Unrecognized residue type %s in PDB file '%s', residue ID '%s'." % (residue_type, str(self.pdb_id), str(residue_id)))

                #structural_residue_IDs.append((residue_id, short_residue_type))
                # KAB - way to allow for multiresidue noncanonical AA's
                if len(short_residue_type) == 1:
                    atom_sequences[chain_id].add(PDBResidue(residue_id[0], residue_id[1:], short_residue_type, chain_type))
                else:
                    for char in short_residue_type:
                        atom_sequences[chain_id].add(PDBResidue(residue_id[0], residue_id[1:], char, chain_type))

        # Assign 'Ligand' or 'Heterogen' to all HETATM-only chains
        for chain_id in present_chain_ids.keys():
            if chain_id not in self.chain_types:
                assert('ATOM  ' not in present_chain_ids[chain_id])
                self.chain_types[chain_id] = PDB._determine_heterogen_chain_type(hetatm_map.get(chain_id, set()))
        self.atom_sequences = atom_sequences
483379
def _extract_blocks(x, block_h, block_w):
  
  (_, height, width, depth) = common_layers.shape_list(x)
  assert height % block_h == 0
  assert width % block_w == 0
  x = tf.reshape(x, [-1, height//block_h, block_h,
                     width//block_w, block_w, depth])
  return tf.transpose(x, [0, 1, 3, 2, 4, 5])
288926
private static String calculateMolecularFormula(List<AbstractMolecule> molecules) throws BuilderMoleculeException, CTKException, ChemistryException {
    Map<String, Integer> atomNumberMap = new TreeMap<String, Integer>();
    for (AbstractMolecule molecule : molecules) {
      LOG.info(molecule.getMolecule().toString());
      atomNumberMap = generateAtomNumberMap(molecule, atomNumberMap);
    }
    LOG.info("GET map");
    StringBuilder sb = new StringBuilder();
    Set<String> atoms = atomNumberMap.keySet();
    for (Iterator<String> i = atoms.iterator(); i.hasNext();) {
      String atom = i.next();
      String num = atomNumberMap.get(atom).toString();
      if (num.equals("1")) {
        num = "";
      }
      sb.append(atom);
      sb.append(num.toString());
    }
    return sb.toString();
  }
70614
def invariants(mol):
    """"""
    atom_ids = {}
    for a in mol.atoms:
        components = []
        components.append(a.number)
        components.append(len(a.oatoms))
        components.append(a.hcount)
        components.append(a.charge)
        components.append(a.mass)
        if len(a.rings) > 0:
            components.append(1)

        atom_ids[a.index] = gen_hash(components)

    return atom_ids
320084
def getpexts(lsp):
    '''
    
    '''
    lines=lsp.split('\n');
    #unfortunately regex doesn't work here
    lns,planens = zip(
        *[ (i,int(re.search('^ *extract *([0-9]+)',line).group(1)))
           for i,line in enumerate(lines)
           if re.search('^ *extract *[0-9]+', line)]);
    if len(lns) == 0: return [];
    end = lns[-1];
    
    for i,line in enumerate(lines[end+1:]):
        if re.match(' *\[',line): break;
    end += i;
    lineranges = zip(lns,(lns+(end,))[1:]);
    planes=dict()
    for (i,end),plane in zip(lineranges,planens):
        d=dict();
        labels = [
            'species',
            'direction',
            'position',];
        datarx = [
            '^ *species *([0-9]+)',
            '^ *direction *([xXyYzZ])',
            '^ *at *(.*)',];
        convs  = [
            lambda s: int(s),
            lambda i: i,
            lambda s: np.array(
                map(float,s.split(' '))),
        ];
        for line in lines[i:end]:
            for label,rx,conv in zip(labels,datarx,convs):
                if re.match(rx,line):
                    d[label]=conv(re.match(rx,line).group(1));
                pass
            pass
        planes[plane] = d;
    return planes;
1062529
def _ncores(meshes, twod):
    """"""
    nnpb = len(meshes)  # number of nodes per block
    nns = [1, 1, 1]  # number of nodes in x, y, z directions
    if twod is None or 'X' in twod:
        while (nnpb > 1 and
               meshes[nns[0]]['X'][0, 0, 0] ==
               meshes[nns[0] - 1]['X'][-1, 0, 0]):
            nns[0] += 1
            nnpb -= 1
    cpu = lambda icy: icy * nns[0]
    if twod is None or 'Y' in twod:
        while (nnpb > 1 and
               meshes[cpu(nns[1])]['Y'][0, 0, 0] ==
               meshes[cpu(nns[1] - 1)]['Y'][0, -1, 0]):
            nns[1] += 1
            nnpb -= nns[0]
    cpu = lambda icz: icz * nns[0] * nns[1]
    while (nnpb > 1 and
           meshes[cpu(nns[2])]['Z'][0, 0, 0] ==
           meshes[cpu(nns[2] - 1)]['Z'][0, 0, -1]):
        nns[2] += 1
        nnpb -= nns[0] * nns[1]
    return np.array(nns)
809650
def sequence(arcs):
    """"""
    succ = {}
    for (i,j) in arcs:
        succ[i] = j
    curr = 1    # first node being visited
    sol = [curr]
    for i in range(len(arcs)-2):
        curr = succ[curr]
        sol.append(curr)
    return sol
162832
def neighbors2(self, distance, chain_residue, atom = None, resid_list = None):

        #atom = " CA "
        ''''''

        if atom == None:     # consider all atoms
            lines = [line for line in self.atomlines(resid_list) if line[17:20] in allowed_PDB_residues_types]
        else:                # consider only given atoms
            lines = [line for line in self.atomlines(resid_list) if line[17:20] in allowed_PDB_residues_types and line[12:16] == atom]

        shash = spatialhash.SpatialHash(distance)

        for line in lines:
            pos = (float(line[30:38]), float(line[38:46]), float(line[46:54]))
            shash.insert(pos, line[21:26])

        neighbor_list = []
        for line in lines:
            resid = line[21:26]
            if resid == chain_residue:
                pos = (float(line[30:38]), float(line[38:46]), float(line[46:54]))

                for data in shash.nearby(pos, distance):
                    if data[1] not in neighbor_list:
                        neighbor_list.append(data[1])
                neighbor_list.sort()
        return neighbor_list
344895
def get_structure(atoms, cls=None):
        
        symbols = atoms.get_chemical_symbols()
        positions = atoms.get_positions()
        lattice = atoms.get_cell()

        cls = Structure if cls is None else cls
        return cls(lattice, symbols, positions,
                   coords_are_cartesian=True)