Compute the electrostatic potential at each atom due to other atoms
328758
def efield(self):
        """"""
        result = np.zeros((self.numc,3), float)
        for index1 in range(self.numc):
            result[index1] = self.efield_component(index1)
        return result
716800
private static double getRadiusForAmino(AminoAcid amino, Atom atom) {

		if (atom.getElement().equals(Element.H)) return Element.H.getVDWRadius();
		// some unusual entries (e.g. 1tes) contain Deuterium atoms in standard aminoacids
		if (atom.getElement().equals(Element.D)) return Element.D.getVDWRadius();

		String atomCode = atom.getName();
		char aa = amino.getAminoType();

		// here we use the values that Chothia gives in his paper (as NACCESS does)
		if (atom.getElement()==Element.O) {
			return OXIGEN_VDW;
		}
		else if (atom.getElement()==Element.S) {
			return SULFUR_VDW;
		}
		else if (atom.getElement()==Element.N) {
			if (atomCode.equals("NZ")) return TETRAHEDRAL_NITROGEN_VDW; // tetrahedral Nitrogen
			return TRIGONAL_NITROGEN_VDW;								// trigonal Nitrogen
		}
		else if (atom.getElement()==Element.C) { // it must be a carbon
			if (atomCode.equals("C") ||
					atomCode.equals("CE1") || atomCode.equals("CE2") || atomCode.equals("CE3") ||
					atomCode.equals("CH2") ||
					atomCode.equals("CZ") || atomCode.equals("CZ2") || atomCode.equals("CZ3")) {
				return TRIGONAL_CARBON_VDW; 							// trigonal Carbon
			}
			else if (atomCode.equals("CA") || atomCode.equals("CB") ||
					atomCode.equals("CE") ||
					atomCode.equals("CG1") || atomCode.equals("CG2")) {
				return TETRAHEDRAL_CARBON_VDW;							// tetrahedral Carbon
			}
			// the rest of the cases (CD, CD1, CD2, CG) depend on amino acid
			else {
				switch (aa) {
				case 'F':
				case 'W':
				case 'Y':
				case 'H':
				case 'D':
				case 'N':
					return TRIGONAL_CARBON_VDW;

				case 'P':
				case 'K':
				case 'R':
				case 'M':
				case 'I':
				case 'L':
					return TETRAHEDRAL_CARBON_VDW;

				case 'Q':
				case 'E':
					if (atomCode.equals("CD")) return TRIGONAL_CARBON_VDW;
					else if (atomCode.equals("CG")) return TETRAHEDRAL_CARBON_VDW;

				default:
					logger.info("Unexpected carbon atom "+atomCode+" for aminoacid "+aa+", assigning its standard vdw radius");
					return Element.C.getVDWRadius();
				}
			}

			// not any of the expected atoms
		} else {
			// non standard aas, (e.g. MSE, LLP) will always have this problem,
			logger.info("Unexpected atom "+atomCode+" for aminoacid "+aa+ " ("+amino.getPDBName()+"), assigning its standard vdw radius");


			return atom.getElement().getVDWRadius();
		}
	}
352925
def from_ase_atoms(cls, atoms):
        
        return cls(atoms=atoms.get_chemical_symbols(), coords=atoms.positions)
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
228239
def neighbor_atoms(self, atom_symbol=None):
        
        if not atom_symbol:
            return self.neighbors
        else:
            return [atom for atom in self.neighbors if atom['atom_symbol'] == atom_symbol]
882779
def potentials_difference(p1, p2):
    

    max_rdiff = 0.0
    np = len(p1)
    if len(p2) != np:
        print("Different number of potentials")
        return float('inf')

    pots1 = sort_potentials(p1)
    pots2 = sort_potentials(p2)

    for n in range(np):
        pot1 = pots1[n]
        pot2 = pots2[n]

        if pot1['angular_momentum'] != pot2['angular_momentum']:
            print("Different angular momentum for potential {}".format(n))
            return float('inf')

        nprim = len(pot1['gaussian_exponents'])
        if len(pot2['gaussian_exponents']) != nprim:
            print("Different number of primitives for potential {}".format(n))
            return float('inf')

        ngen = len(pot1['coefficients'])
        if len(pot2['coefficients']) != ngen:
            print("Different number of general contractions for potential {}".format(n))
            return float('inf')

        for p in range(nprim):
            e1 = pot1['gaussian_exponents'][p]
            e2 = pot2['gaussian_exponents'][p]
            r = _reldiff(e1, e2)
            if r > 0.0:
                print("   Gaussian Exponent {:3}: {:20} {:20} -> {:16.8e}".format(p, e1, e2, r))
            max_rdiff = max(max_rdiff, r)

            e1 = pot1['r_exponents'][p]
            e2 = pot2['r_exponents'][p]
            r = _reldiff(e1, e2)
            if r > 0.0:
                print("          R Exponent {:3}: {:20} {:20} -> {:16.8e}".format(p, e1, e2, r))
            max_rdiff = max(max_rdiff, r)

            for g in range(ngen):
                c1 = pot1['coefficients'][g][p]
                c2 = pot2['coefficients'][g][p]
                r = _reldiff(c1, c2)
                if r > 0.0:
                    print("         Coefficient {:3}: {:20} {:20} -> {:16.8e}".format(p, c1, c2, r))
                max_rdiff = max(max_rdiff, r)

    print()
    print("Max relative difference for these potentials: {}".format(max_rdiff))
    return max_rdiff
102025
def ρ(self, e):
        
        d = self._dichalcogenide
        at, λ = d.at, d.λ
        return abs(2 * e - λ) * (2 * at**2)**(-1)
813426
def write_cfour(basis):
    '''
    '''

    # March 2019
    # Format determined from http://slater.chemie.uni-mainz.de/cfour/index.php?n=Main.NewFormatOfAnEntryInTheGENBASFile

    # Uncontract all, then make general
    basis = manip.uncontract_spdf(basis, 0, True)
    basis = manip.make_general(basis, False)
    basis = sort.sort_basis(basis, False)

    # Elements for which we have electron basis
    electron_elements = [k for k, v in basis['elements'].items() if 'electron_shells' in v]

    # Elements for which we have ECP
    ecp_elements = [k for k, v in basis['elements'].items() if 'ecp_potentials' in v]

    s = '\n'

    if len(electron_elements) > 0:
        # Electron Basis
        for z in electron_elements:
            data = basis['elements'][z]
            sym = lut.element_sym_from_Z(z).upper()
            nshell = len(data['electron_shells'])

            s += '{}:{}\n'.format(sym, basis['name'])
            s += basis['description'] + '\n'
            s += '\n'
            s += '{:>3}\n'.format(nshell)

            s_am = ''
            s_ngen = ''
            s_nprim = ''
            for sh in data['electron_shells']:
                s_am += '{:>5}'.format(sh['angular_momentum'][0])
                s_ngen += '{:>5}'.format(len(sh['coefficients']))
                s_nprim += '{:>5}'.format(len(sh['exponents']))

            s += s_am + '\n'
            s += s_ngen + '\n'
            s += s_nprim + '\n'
            s += '\n'

            for shell in data['electron_shells']:
                exponents = shell['exponents']
                coefficients = shell['coefficients']
                coefficients = list(map(list, zip(*coefficients)))

                s += '  '.join(exponents).replace("E", "D") + '\n\n'
                for c in coefficients:
                    s += '  '.join(c).replace("E", "D") + '\n'
                s += '\n'

    # Write out ECP
    if len(ecp_elements) > 0:
        s += '\n\n! Effective core Potentials\n'

        for z in ecp_elements:
            data = basis['elements'][z]
            sym = lut.element_sym_from_Z(z).upper()
            max_ecp_am = max([x['angular_momentum'][0] for x in data['ecp_potentials']])
            max_ecp_amchar = lut.amint_to_char([max_ecp_am]).lower()

            # Sort lowest->highest, then put the highest at the beginning
            ecp_list = sorted(data['ecp_potentials'], key=lambda x: x['angular_momentum'])
            ecp_list.insert(0, ecp_list.pop())

            s += '*\n'
            s += '{}:{}\n'.format(sym, basis['name'])
            s += '# ' + basis['description'] + '\n'
            s += '*\n'
            s += '    NCORE = {}    LMAX = {}\n'.format(data['ecp_electrons'], max_ecp_am)

            for pot in ecp_list:
                rexponents = pot['r_exponents']
                gexponents = pot['gaussian_exponents']
                coefficients = pot['coefficients']

                am = pot['angular_momentum']
                amchar = lut.amint_to_char(am).lower()

                if am[0] == max_ecp_am:
                    s += '{}\n'.format(amchar)
                else:
                    s += '{}-{}\n'.format(amchar, max_ecp_amchar)

                point_places = [6, 18, 25]
                s += printing.write_matrix([*coefficients, rexponents, gexponents], point_places)
                #for p in range(len(rexponents)):
                #    s += '{}  {}  {};\n'.format(gexponents[p], rexponents[p], coefficients[0][p])
            s += '*\n'
    return s
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
87029
private void calculate(List<String> paths, IAtomContainer mol) {
        AllPairsShortestPaths apsp     = new AllPairsShortestPaths(mol);
        int                   numAtoms = mol.getAtomCount();
        for (int i = 0; i < numAtoms; i++) {
            if (!include(mol.getAtom(i)))
                continue;
            for (int j = i + 1; j < numAtoms; j++) {
                if (!include(mol.getAtom(j))) continue;
                final int dist = apsp.from(i).distanceTo(j);
                if (dist > MAX_DISTANCE)
                    continue;
                final IAtom beg = mol.getAtom(i);
                final IAtom end = mol.getAtom(j);
                paths.add(encodePath(dist, beg, end));
                paths.add(encodePath(dist, end, beg));
                if (isHalogen(mol.getAtom(i)) || isHalogen(mol.getAtom(j))) {
                    paths.add(encodeHalPath(dist, beg, end));
                    paths.add(encodeHalPath(dist, end, beg));
                }
            }
        }
    }
732860
def max_electronegativity(self):
        '''
        
        '''
        maximum = 0
        for e1, e2 in combinations(self.elements, 2):
            if abs(Element(e1).X - Element(e2).X) > maximum:
                maximum = abs(Element(e1).X - Element(e2).X)
        return maximum
13618
private boolean getPresenceOfNitro(IAtomContainer ac, IAtom atom) {
        List neighbours = ac.getConnectedAtomsList(atom);
        List second = null;
        IBond bond = null;
        //int counter = 0;
        for (int i = 0; i < neighbours.size(); i++) {
            IAtom neighbour = (IAtom) neighbours.get(i);
            if (neighbour.getSymbol().equals("N")) {
                second = ac.getConnectedAtomsList(neighbour);
                for (int b = 0; b < second.size(); b++) {
                    IAtom conAtom = (IAtom) second.get(b);
                    if (conAtom.getSymbol().equals("O")) {
                        bond = ac.getBond(neighbour, conAtom);
                        if (bond.getOrder() == IBond.Order.DOUBLE) {
                            return true;
                        }
                    }
                }
            }
        }
        return false;
    }
969439
def dot(self, other):
        """"""
        dot_product = 0
        a = self.elements
        b = other.elements
        a_len = len(a)
        b_len = len(b)
        i = j = 0

        while i < a_len and j < b_len:
            a_val = a[i]
            b_val = b[j]
            if a_val < b_val:
                i += 2
            elif a_val > b_val:
                j += 2
            else:
                dot_product += a[i + 1] * b[j + 1]
                i += 2
                j += 2

        return dot_product
1003711
public void addAtoms(Atom[] atoms, BoundingBox bounds) {
		this.iAtoms = Calc.atomsToPoints(atoms);
		this.iAtomObjects = atoms;

		if (bounds!=null) {
			this.ibounds = bounds;
		} else {
			this.ibounds = new BoundingBox(iAtoms);
		}

		this.jAtoms = null;
		this.jAtomObjects = null;
		this.jbounds = null;

		fillGrid();
	}
1200092
def getPyroleLikeAtoms(cycle):
    """"""    
    result = {}
    # the outgoing bonds might need to be single or aromatic
    for atom in cycle.atoms:
        lookup = (atom.symbol, atom.charge, atom.hcount, len(atom.bonds))
        if PyroleTable.get(lookup, 0):
            result[atom.handle] = atom

    return result
973040
def get_bond_symmetry(site_symmetry,
                      lattice,
                      positions,
                      atom_center,
                      atom_disp,
                      symprec=1e-5):
    
    bond_sym = []
    pos = positions
    for rot in site_symmetry:
        rot_pos = (np.dot(pos[atom_disp] - pos[atom_center], rot.T) +
                   pos[atom_center])
        diff = pos[atom_disp] - rot_pos
        diff -= np.rint(diff)
        dist = np.linalg.norm(np.dot(lattice, diff))
        if dist < symprec:
            bond_sym.append(rot)

    return np.array(bond_sym)
737546
private List<Integer[]> getAtomNoPairsForRingGroup(IAtomContainer molecule, List<Integer> bondsToCheck) {
        List<Integer[]> aptc = new ArrayList<Integer[]>();
        for (Integer i : bondsToCheck) {
            Integer[] aps = new Integer[2];
            aps[0] = molecule.indexOf(molecule.getBond(i).getBegin());
            aps[1] = molecule.indexOf(molecule.getBond(i).getEnd());
            aptc.add(aps);
        }
        return aptc;
    }