Evaluate the PSF containment angle at a sequence of energies.
325985
def containment_angle(self, energies=None, fraction=0.68, scale_fn=None):
        """"""

        if energies is None:
            energies = self.energies

        vals = self.interp(energies[np.newaxis, :], self.dtheta[:, np.newaxis],
                           scale_fn=scale_fn)
        dtheta = np.radians(self.dtheta[:, np.newaxis] * np.ones(vals.shape))
        return self._calc_containment(dtheta, vals, fraction)
142622
def stSpectralFlux(X, X_prev):
    
    # compute the spectral flux as the sum of square distances:
    sumX = numpy.sum(X + eps)
    sumPrevX = numpy.sum(X_prev + eps)
    F = numpy.sum((X / sumX - X_prev/sumPrevX) ** 2)

    return F
526909
def _JRAxiIntegrand(r,E,L,pot):
    """"""
    return nu.sqrt(2.*(E-potentialAxi(r,pot))-L**2./r**2.)
188894
def phase2t(self, psi):
        
        def _deg(rads, domain_lower_limit):
            # Convert rads to degrees in [0, 360) domain
            degs = degrees(rads % (2*pi))

            # Convert to [domain_lower_limit, domain_lower_limit + 360) domain
            k = domain_lower_limit // 360
            degs += k * 360
            if degs < domain_lower_limit:
                degs += 360
            return degs

        if self.delta > 0:
            degs = _deg(psi, domain_lower_limit=self.theta)
        else:
            degs = _deg(psi, domain_lower_limit=self.theta)
        return (degs - self.theta)/self.delta
1239252
def containment_angle_bin(self, egy_bins, fraction=0.68, scale_fn=None):
        """"""

        vals = self.interp_bin(egy_bins, self.dtheta, scale_fn=scale_fn)
        dtheta = np.radians(self.dtheta[:, np.newaxis] * np.ones(vals.shape))
        return self._calc_containment(dtheta, vals, fraction)
631837
def hugoniot_p(rho, rho0, c0, s):
    
    eta = 1. - (rho0 / rho)
    Ph = rho0 * c0 * c0 * eta / np.power((1. - s * eta), 2.)
    return Ph
500277
private int getNearestBondtoAGivenAtom(IAtomContainer mol, IAtom atom, IBond bond) {
        int nearestBond = 0;
        double[] values;
        double distance = 0;
        IAtom atom0 = bond.getBegin();
        List<IBond> bondsAtLeft = mol.getConnectedBondsList(atom0);
        int partial;
        for (int i = 0; i < bondsAtLeft.size(); i++) {
            IBond curBond = bondsAtLeft.get(i);
            values = calculateDistanceBetweenAtomAndBond(atom, curBond);
            partial = mol.indexOf(curBond);
            if (i == 0) {
                nearestBond = mol.indexOf(curBond);
                distance = values[0];
            } else {
                if (values[0] < distance) {
                    nearestBond = partial;
                }
            }
        }
        return nearestBond;

    }
769600
def hessian_component(self, index1, index2):
        """"""
        result = np.zeros((3, 3), float)
        if index1 == index2:
            for index3 in range(self.numc):
                if self.scaling[index1, index3] > 0:
                    d_1 = 1/self.distances[index1, index3]
                    for (se, ve), (sg, vg), (sh, vh) in zip(
                        self.yield_pair_energies(index1, index3),
                        self.yield_pair_gradients(index1, index3),
                        self.yield_pair_hessians(index1, index3)
                    ):
                        result += (
                            +sh*self.dirouters[index1, index3]*ve
                            +sg*(np.identity(3, float) - self.dirouters[index1, index3])*ve*d_1
                            +sg*np.outer(self.directions[index1, index3],  vg)
                            +sg*np.outer(vg, self.directions[index1, index3])
                            +se*vh
                        )*self.scaling[index1, index3]
        elif self.scaling[index1, index2] > 0:
            d_1 = 1/self.distances[index1, index2]
            for (se, ve), (sg, vg), (sh, vh) in zip(
                self.yield_pair_energies(index1, index2),
                self.yield_pair_gradients(index1, index2),
                self.yield_pair_hessians(index1, index2)
            ):
                result -= (
                    +sh*self.dirouters[index1, index2]*ve
                    +sg*(np.identity(3, float) - self.dirouters[index1, index2])*ve*d_1
                    +sg*np.outer(self.directions[index1, index2],  vg)
                    +sg*np.outer(vg, self.directions[index1, index2])
                    +se*vh
                )*self.scaling[index1, index2]
        return result
1310321
def rotational_invariance(force_constants,
                          supercell,
                          primitive,
                          symprec=1e-5):
    

    print("Check rotational invariance ...")

    fc = force_constants
    p2s = primitive.get_primitive_to_supercell_map()

    smallest_vectors, multiplicity = primitive.get_smallest_vectors()

    abc = "xyz"

    for pi, p in enumerate(p2s):
        for i in range(3):
            mat = np.zeros((3, 3), dtype='double')
            for s in range(supercell.get_number_of_atoms()):
                m = multiplicity[s, pi]
                vecs = smallest_vectors[s, pi, :m]
                v = np.dot(vecs.sum(axis=0) / m, primitive.get_cell())
                for j in range(3):
                    for k in range(3):
                        mat[j, k] += (fc[p, s, i, j] * v[k] -
                                      fc[p, s, i, k] * v[j])

            print("Atom %d %s" % (p + 1, abc[i]))
            for vec in mat:
                print("%10.5f %10.5f %10.5f" % tuple(vec))
1220674
def fswap(p, q):
    

    yield cirq.ISWAP(q, p), cirq.Z(p) ** 1.5
    yield cirq.Z(q) ** 1.5
42526
def _parse_forces(line, lines):
    """"""
    units = line.split()[4].rstrip(":")
    next(lines)
    newline = next(lines)
    total = []; non_local = []; ionic = []; local = []; core_correction = []
    hubbard = []; scf = []; types = []
    while (not "The non-local contrib." in newline) and len(newline.split()) > 0:
        if "=" in newline:
            total.append([float(x) for x in newline.partition("=")[2].split()])
            types.append(int(newline.split()[3]))
        newline = next(lines)

    if len(newline.split()) > 0:
        while not "The ionic contribution" in newline:
            if "=" in newline:
                non_local.append([float(x) for x in newline.partition("=")[2].split()])
            newline = next(lines)
        while not "The local contribution" in newline:
            if "=" in newline:
                ionic.append([float(x) for x in newline.partition("=")[2].split()])
            newline = next(lines)
        while not "The core correction contribution" in newline:
            if "=" in newline:
                local.append([float(x) for x in newline.partition("=")[2].split()])
            newline = next(lines)
        while not "The Hubbard contrib." in newline:
            if "=" in newline:
                core_correction.append([float(x) for x in newline.partition("=")[2].split()])
            newline = next(lines)
        while not "The SCF correction term" in newline:
            if "=" in newline:
                hubbard.append([float(x) for x in newline.partition("=")[2].split()])
            newline = next(lines)
        while len(newline.split()) > 0:
            if "=" in newline:
                scf.append([float(x) for x in newline.partition("=")[2].split()])
            newline = next(lines)
    newline = next(lines)
    total_force = float(newline.split()[3])
    total_scf = float(newline.split()[8])
    return {
        "force units": units,
        "total force": total_force,
        "total SCF correction": total_scf,
        "forces": total,
        "non-local contribution to forces": non_local,
        "ionic contribution to forces": ionic,
        "local contribution to forces": local,
        "core corrections to forces": core_correction,
        "Hubbard contribution to forces": hubbard,
        "SCF correction term to forces": scf,
        "Atomic species index for forces": types
    }
763259
def cal_pst(self, v):
        
        p_h = self._hugoniot_p(v)
        p_th_h = self._hugoniot_pth(v)
        p_st = p_h - p_th_h
        return p_st
785381
def dpsi2_dtheta(self, dL_dpsi2, Z, mu, S, target):
        """"""
        self._psi_computations(Z, mu, S)
        d_var = 2.*self._psi2 / self.variance
        # d_length = 2.*self._psi2[:, :, :, None] * (self._psi2_Zdist_sq * self._psi2_denom + self._psi2_mudist_sq + S[:, None, None, :] / self.lengthscale2) / (self.lengthscale * self._psi2_denom)
        d_length = -2.*self._psi2[:, :, :, None] * (self._psi2_Zdist_sq * self._psi2_denom + self._psi2_mudist_sq + S[:, None, None, :] * self.inv_lengthscale2) / (self.inv_lengthscale * self._psi2_denom)
        target[0] += np.sum(dL_dpsi2 * d_var)
        dpsi2_dlength = d_length * dL_dpsi2[:, :, :, None]
        if not self.ARD:
            target[1] += dpsi2_dlength.sum() # *(-self.lengthscale2)
        else:
            target[1:] += dpsi2_dlength.sum(0).sum(0).sum(0)
1051020
@Override
    public DescriptorValue calculate(IAtomContainer container) {
        if (!GeometryUtil.has3DCoordinates(container))
            return getDummyDescriptorValue(new CDKException("Molecule must have 3D coordinates"));

        IsotopeFactory factory = null;
        double mass1;
        double mass2;
        try {
            factory = Isotopes.getInstance();
        } catch (Exception e) {
            logger.debug(e);
        }

        double sum = 0;
        for (int i = 0; i < container.getBondCount(); i++) {
            IBond bond = container.getBond(i);

            if (bond.getAtomCount() != 2) {
                return getDummyDescriptorValue(new CDKException("GravitationalIndex: Only handles 2 center bonds"));
            }

            mass1 = factory.getMajorIsotope(bond.getBegin().getSymbol()).getMassNumber();
            mass2 = factory.getMajorIsotope(bond.getEnd().getSymbol()).getMassNumber();

            Point3d p1 = bond.getBegin().getPoint3d();
            Point3d p2 = bond.getEnd().getPoint3d();

            double x1 = p1.x;
            double y1 = p1.y;
            double z1 = p1.z;
            double x2 = p2.x;
            double y2 = p2.y;
            double z2 = p2.z;

            double dist = (x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2) + (z1 - z2) * (z1 - z2);
            sum += (mass1 * mass2) / dist;
        }

        // heavy atoms only
        double heavysum = 0;
        for (int i = 0; i < container.getBondCount(); i++) {
            IBond b = container.getBond(i);

            if (b.getAtomCount() != 2) {
                return getDummyDescriptorValue(new CDKException("GravitationalIndex: Only handles 2 center bonds"));
            }

            if (b.getBegin().getSymbol().equals("H") || b.getEnd().getSymbol().equals("H")) continue;

            mass1 = factory.getMajorIsotope(b.getBegin().getSymbol()).getMassNumber();
            mass2 = factory.getMajorIsotope(b.getEnd().getSymbol()).getMassNumber();

            Point3d point0 = b.getBegin().getPoint3d();
            Point3d point1 = b.getEnd().getPoint3d();

            double x1 = point0.x;
            double y1 = point0.y;
            double z1 = point0.z;
            double x2 = point1.x;
            double y2 = point1.y;
            double z2 = point1.z;

            double dist = (x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2) + (z1 - z2) * (z1 - z2);
            heavysum += (mass1 * mass2) / dist;
        }

        // all pairs
        ArrayList<Integer> x = new ArrayList<Integer>();
        for (int i = 0; i < container.getAtomCount(); i++) {
            if (!container.getAtom(i).getSymbol().equals("H")) x.add(i);
        }
        int npair = x.size() * (x.size() - 1) / 2;
        pair[] p = new pair[npair];
        for (int i = 0; i < npair; i++)
            p[i] = new pair();
        int pcount = 0;
        for (int i = 0; i < x.size() - 1; i++) {
            for (int j = i + 1; j < x.size(); j++) {
                int present = 0;
                int a = x.get(i);
                int b = x.get(j);
                for (int k = 0; k < pcount; k++) {
                    if ((p[k].x == a && p[k].y == b) || (p[k].y == a && p[k].x == b)) present = 1;
                }
                if (present == 1) continue;
                p[pcount].x = a;
                p[pcount].y = b;
                pcount += 1;
            }
        }
        double allheavysum = 0;
        for (pair aP : p) {
            int atomNumber1 = aP.x;
            int atomNumber2 = aP.y;

            mass1 = factory.getMajorIsotope(container.getAtom(atomNumber1).getSymbol()).getMassNumber();
            mass2 = factory.getMajorIsotope(container.getAtom(atomNumber2).getSymbol()).getMassNumber();

            double x1 = container.getAtom(atomNumber1).getPoint3d().x;
            double y1 = container.getAtom(atomNumber1).getPoint3d().y;
            double z1 = container.getAtom(atomNumber1).getPoint3d().z;
            double x2 = container.getAtom(atomNumber2).getPoint3d().x;
            double y2 = container.getAtom(atomNumber2).getPoint3d().y;
            double z2 = container.getAtom(atomNumber2).getPoint3d().z;

            double dist = (x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2) + (z1 - z2) * (z1 - z2);
            allheavysum += (mass1 * mass2) / dist;
        }

        DoubleArrayResult retval = new DoubleArrayResult(9);
        retval.add(heavysum);
        retval.add(Math.sqrt(heavysum));
        retval.add(Math.pow(heavysum, 1.0 / 3.0));

        retval.add(sum);
        retval.add(Math.sqrt(sum));
        retval.add(Math.pow(sum, 1.0 / 3.0));

        retval.add(allheavysum);
        retval.add(Math.sqrt(allheavysum));
        retval.add(Math.pow(allheavysum, 1.0 / 3.0));

        return new DescriptorValue(getSpecification(), getParameterNames(), getParameters(), retval,
                getDescriptorNames());
    }
664256
def getbfs(coords, gbasis):
    """"""

    sym2powerlist = {
        'S' : [(0,0,0)],
        'P' : [(1,0,0),(0,1,0),(0,0,1)],
        'D' : [(2,0,0),(0,2,0),(0,0,2),(1,1,0),(0,1,1),(1,0,1)],
        'F' : [(3,0,0),(2,1,0),(2,0,1),(1,2,0),(1,1,1),(1,0,2),
               (0,3,0),(0,2,1),(0,1,2), (0,0,3)]
        }

    bfs = []
    for i, at_coords in enumerate(coords):
        bs = gbasis[i]
        for sym,prims in bs:
            for power in sym2powerlist[sym]:
                bf = cgbf(at_coords,power)
                for expnt,coef in prims:
                    bf.add_pgbf(expnt,coef)
                bf.normalize()
                bfs.append(bf)

    return bfs
1325181
private void calcInterfacesCrystal(StructureInterfaceList set, double cutoff) {


		// initialising debugging vars
		long start = -1;
		long end = -1;
		int trialCount = 0;
		int skippedRedundant = 0;
		int skippedAUsNoOverlap = 0;
		int skippedChainsNoOverlap = 0;
		int skippedSelfEquivalent = 0;

		// The bounding boxes of all AUs of the unit cell
		UnitCellBoundingBox bbGrid = new UnitCellBoundingBox(numOperatorsSg, numPolyChainsAu);;
		// we calculate all the bounds of each of the asym units, those will then be reused and translated
		bbGrid.setBbs(structure, ops, INCLUDE_HETATOMS);


		// if not crystallographic there's no search to do in other cells, only chains within "AU" will be checked
		if (!searchBeyondAU) numCells = 0;

		boolean verbose = logger.isDebugEnabled();

		if (verbose) {
			trialCount = 0;
			start= System.currentTimeMillis();
			int neighbors = (2*numCells+1)*(2*numCells+1)*(2*numCells+1)-1;
			int auTrials = (numPolyChainsAu*(numPolyChainsAu-1))/2;
			int trials = numPolyChainsAu*numOperatorsSg*numPolyChainsAu*neighbors;
			logger.debug("Chain clash trials within original AU: "+auTrials);
			logger.debug(
					"Chain clash trials between the original AU and the neighbouring "+neighbors+
					" whole unit cells ("+numCells+" neighbours)" +
					"(2x"+numPolyChainsAu+"chains x "+numOperatorsSg+"AUs x "+neighbors+"cells) : "+trials);
			logger.debug("Total trials: "+(auTrials+trials));
		}

		List<Chain> polyChains = structure.getPolyChains();

		for (int a=-numCells;a<=numCells;a++) {
			for (int b=-numCells;b<=numCells;b++) {
				for (int c=-numCells;c<=numCells;c++) {

					Point3i trans = new Point3i(a,b,c);
					Vector3d transOrth = new Vector3d(a,b,c);
					if (a!=0 || b!=0 || c!=0) {
						// we avoid doing the transformation for 0,0,0 (in case it's not crystallographic)
						this.crystallographicInfo.getCrystalCell().transfToOrthonormal(transOrth);
					}

					UnitCellBoundingBox bbGridTrans = bbGrid.getTranslatedBbs(transOrth);

					for (int n=0;n<numOperatorsSg;n++) {

						// short-cut strategies
						// 1) we skip first of all if the bounding boxes of the AUs don't overlap
						if (!bbGrid.getAuBoundingBox(0).overlaps(bbGridTrans.getAuBoundingBox(n), cutoff)) {
							skippedAUsNoOverlap++;
							continue;
						}

						// 2) we check if we didn't already see its equivalent symmetry operator partner
						CrystalTransform tt = new CrystalTransform(this.crystallographicInfo.getSpaceGroup(), n);
						tt.translate(trans);
						if (isRedundantTransform(tt)) {
							skippedRedundant++;
							continue;
						}
						addVisitedTransform(tt);


						boolean selfEquivalent = false;

						// 3) an operator can be "self redundant" if it is the inverse of itself (involutory, e.g. all pure 2-folds with no translation)
						if (tt.isEquivalent(tt)) {
							logger.debug("Transform {} is equivalent to itself, will skip half of i-chains to j-chains comparisons", tt.toString());
							// in this case we can't skip the operator, but we can skip half of the matrix comparisons e.g. j>i
							// we set a flag and do that within the loop below
							selfEquivalent = true;
						}

						StringBuilder builder = null;
						if (verbose) builder = new StringBuilder(String.valueOf(tt)).append(" ");

						// Now that we know that boxes overlap and operator is not redundant, we have to go to the details
						int contactsFound = 0;

						for (int j=0;j<numPolyChainsAu;j++) {

							for (int i=0;i<numPolyChainsAu;i++) { // we only have to compare the original asymmetric unit to every full cell around

								if(selfEquivalent && (j>i)) {
									// in case of self equivalency of the operator we can safely skip half of the matrix
									skippedSelfEquivalent++;
									continue;
								}
								// special case of original AU, we don't compare a chain to itself
								if (n==0 && a==0 && b==0 && c==0 && i==j) continue;

								// before calculating the AtomContactSet we check for overlap, then we save putting atoms into the grid
								if (!bbGrid.getChainBoundingBox(0,i).overlaps(bbGridTrans.getChainBoundingBox(n,j),cutoff)) {
									skippedChainsNoOverlap++;
									if (verbose) {
										builder.append(".");
									}
									continue;
								}

								trialCount++;

								// finally we've gone through all short-cuts and the 2 chains seem to be close enough:
								// we do the calculation of contacts
								Chain chaini = polyChains.get(i);
								Chain chainj = polyChains.get(j);

								if (n!=0 || a!=0 || b!=0 || c!=0) {
									Matrix4d mJCryst = new Matrix4d(ops[n]);
									translate(mJCryst, transOrth);
									chainj = (Chain)chainj.clone();
									Calc.transform(chainj,mJCryst);
								}

								StructureInterface interf = calcContacts(chaini, chainj, cutoff, tt, builder);
								if (interf == null) {
									continue;
								}

								contactsFound++;
								if(this.hasNcsOps()) {
									StructureInterface interfNcsRef = findNcsRef(interf);
									set.addNcsEquivalent(interf,interfNcsRef);
								} else {
									set.add(interf);
								}
							}
						}

						if( verbose ) {
							if (a==0 && b==0 && c==0 && n==0)
								builder.append(" "+contactsFound+"("+(numPolyChainsAu*(numPolyChainsAu-1))/2+")");
							else if (selfEquivalent)
								builder.append(" "+contactsFound+"("+(numPolyChainsAu*(numPolyChainsAu+1))/2+")");
							else
								builder.append(" "+contactsFound+"("+numPolyChainsAu*numPolyChainsAu+")");

							logger.debug(builder.toString());
						}
					}
				}
			}
		}

		end = System.currentTimeMillis();
		logger.debug("\n"+trialCount+" chain-chain clash trials done. Time "+(end-start)/1000+"s");
		logger.debug("  skipped (not overlapping AUs)       : "+skippedAUsNoOverlap);
		logger.debug("  skipped (not overlapping chains)    : "+skippedChainsNoOverlap);
		logger.debug("  skipped (sym redundant op pairs)    : "+skippedRedundant);
		logger.debug("  skipped (sym redundant self op)     : "+skippedSelfEquivalent);
		logger.debug("Found "+set.size()+" interfaces.");
	}
705375
def spectral_rolloff(data, coeff):
    """"""
    data = np.mean(data, axis=1)
    nFFT = len(data) // 2
    X = FFT(data, nFFT)

    totalEnergy = np.sum(X ** 2)
    fftLength = len(X)
    Thres = coeff * totalEnergy

    # Find the spectral rolloff as the frequency position where the
    # respective spectral energy is equal to c*totalEnergy
    CumSum = np.cumsum(X ** 2) + EPSILON
    [a, ] = np.nonzero(CumSum > Thres)
    if len(a) > 0:
        mC = np.float64(a[0]) / (float(fftLength))
    else:
        mC = 0.0

    return mC