calculate Debye energy using old fortran routine

    :params x: Debye x value
    :return: Debye energy
498564
def debye_E_single(x):
    
    # make the function handles both scalar and array
    if ((x > 0.0) & (x <= 0.1)):
        result = 1. - 0.375 * x + x * x * \
            (0.05 - (5.952380953e-4) * x * x)
    # for 0.1 < x <= 7.25
    if ((x > 0.1) & (x <= 7.25)):
        result = ((((.0946173 * x - 4.432582) * x +
                    85.07724) * x - 800.6087) * x +
                  3953.632) / ((((x + 15.121491) * x +
                                 143.155337) * x + 682.0012) *
                               x + 3953.632)
    # for x > 7.25
    # it appears there might be error for this part, but never been exposed
    # because of rarity of such high x value.
    if (x > 7.25):
        exx = np.exp(-x)
        nn = np.round(25. / x)
        n = nn.astype(np.int64)
        temp = 0.
        if (n > 0):
            temp2 = 1.
            end = n + 1
            for i in range(1, end):
                temps = i * 1.
                temp2 = temp2 * exx
                x3 = temps * x
                temp = temp + temp2 * \
                    (6. + x3 * (6. + x3 * (3. + x3))) / \
                    (temps * temps * temps * temps)
        result = 3.0 * (6.493939402 - temp) / (x * x * x)
    return result
1324465
def energy(self):
        """"""
        e = 0
        for i in range(len(self.state)):
            e += self.distance_matrix[self.state[i-1]][self.state[i]]
        return e
300399
def _calc_ML(sampler, modelidx=0, e_range=None, e_npoints=100):
    

    ML, MLp, MLerr, ML_model = find_ML(sampler, modelidx)

    if e_range is not None:
        # prepare bogus data for calculation
        e_range = validate_array(
            "e_range", u.Quantity(e_range), physical_type="energy"
        )
        e_unit = e_range.unit
        energy = (
            np.logspace(
                np.log10(e_range[0].value),
                np.log10(e_range[1].value),
                e_npoints,
            )
            * e_unit
        )
        data = {
            "energy": energy,
            "flux": np.zeros(energy.shape) * sampler.data["flux"].unit,
        }
        modelout = sampler.modelfn(MLp, data)

        if isinstance(modelout, np.ndarray):
            blob = modelout
        else:
            blob = modelout[modelidx]

        if isinstance(blob, u.Quantity):
            modelx = data["energy"].copy()
            model_ML = blob.copy()
        elif len(blob) == 2:
            modelx = blob[0].copy()
            model_ML = blob[1].copy()
        else:
            raise TypeError("Model {0} has wrong blob format".format(modelidx))

        ML_model = (modelx, model_ML)

    return ML, MLp, MLerr, ML_model
1272364
def energy(mesh, uniform_density=False):
    
    # E = 1/(d+1) sum_i ||x_i||^2 |omega_i| - int_Omega_i ||x||^2
    dim = mesh.cells["nodes"].shape[1] - 1

    star_volume = numpy.zeros(mesh.node_coords.shape[0])
    for i in range(3):
        idx = mesh.cells["nodes"][:, i]
        if uniform_density:
            # rho = 1,
            # int_{star} phi_i * rho = 1/(d+1) sum_{triangles in star} |triangle|
            fastfunc.add.at(star_volume, idx, mesh.cell_volumes)
        else:
            # rho = 1 / tau_j,
            # int_{star} phi_i * rho = 1/(d+1) |num triangles in star|
            fastfunc.add.at(star_volume, idx, numpy.ones(idx.shape, dtype=float))
    x2 = numpy.einsum("ij,ij->i", mesh.node_coords, mesh.node_coords)
    out = 1 / (dim + 1) * numpy.dot(star_volume, x2)

    # could be cached
    assert dim == 2
    x = mesh.node_coords[:, :2]
    triangles = numpy.moveaxis(x[mesh.cells["nodes"]], 0, 1)
    val = quadpy.triangle.integrate(
        lambda x: x[0] ** 2 + x[1] ** 2,
        triangles,
        # Take any scheme with order 2
        quadpy.triangle.Dunavant(2),
    )
    if uniform_density:
        val = numpy.sum(val)
    else:
        rho = 1.0 / mesh.cell_volumes
        val = numpy.dot(val, rho)

    assert out >= val

    return out - val
627934
def ParallelCalculate(cls,syslst,properties=['energy'],system_changes=all_changes):
        '''
        
        '''
        print('Launching:',end=' ')
        sys.stdout.flush()
        for n,s in enumerate(syslst):
            try :
                s.calc.block=False
                s.calc.calculate(atoms=s,properties=properties,system_changes=system_changes)
            except CalcNotReadyError:
                s.calc.block=True
            print(n+1, end=' ')
            sys.stdout.flush()
        print()
        print('     Done:', end=' ')
        sys.stdout.flush()
        for n,s in enumerate(syslst):
            s.calc.read_results()
            print( n+1, end=' ')
            sys.stdout.flush()
        print()
        return syslst
1173910
def Iaxn(mt, x, n, *args):
    
    return (Sx(mt, x + 1) - Sx(mt, x + n + 1) - n * Nx(mt, x + n + 1)) / Dx(mt, x)
284464
def get_potential_energy(self, a=None):
        """"""
        if a is None:
            a = self.a
        if self.mask is None:
            return -np.sum(np.dot(a.get_positions()[self.mask], self.force))
        else:
            return -np.sum(np.dot(a.get_positions(), self.force))
699384
double getBetheFreeEnergy() {
        // 
        // G_{Bethe} = \sum_a \sum_{x_a} - b(x_a) ln \chi(x_a)
        //              + \sum_a \sum_{x_a} b(x_a) ln b(x_a)
        //              + \sum_i (n_i - 1) \sum_{x_i} b(x_i) ln b(x_i)
        //           = \sum_a \sum_{x_a} b(x_a) ln (b(x_a) / \chi(x_a))
        //              + \sum_i (n_i - 1) \sum_{x_i} b(x_i) ln b(x_i)
        //
        //     where n_i is the number of neighbors of the variable x_i,
        //     b(x_a) and b(x_i) are normalized distributions and x_a is 
        //     the set of variables participating in factor a. 
        //
                
        double bethe = 0.0;
        for (int a=0; a<fg.getFactors().size(); a++) {
            Factor f = fg.getFactors().get(a);
            if (!(f instanceof GlobalFactor)) {
                int numConfigs = f.getVars().calcNumConfigs();
                VarTensor beliefs = getFactorBeliefs(a);
                for (int c=0; c<numConfigs; c++) {                
                    // Since we want multiplication by 0 to always give 0 (not the case for Double.POSITIVE_INFINITY or Double.NaN.
                    double b_c = beliefs.getValue(c);
                    if (b_c != s.zero()) {
                        double r_b_c = s.toReal(b_c);
                        double log_b_c = s.toLogProb(b_c);
                        double log_chi_c = f.getLogUnormalizedScore(c);
                        bethe += r_b_c * (log_b_c - log_chi_c);
                    }
                }
            } else {
                VarTensor[] inMsgs = getMsgs(f.getId(), msgs, IN_MSG);
                bethe += ((GlobalFactor) f).getExpectedLogBelief(inMsgs);
            }
        }
        for (int v=0; v<fg.getVars().size(); v++) {
            Var var = fg.getVars().get(v);
            int numNeighbors = bg.numNbsT1(v);
            VarTensor beliefs = getVarBeliefs(v);
            double sum = 0.0;
            for (int c=0; c<var.getNumStates(); c++) {
                double b_c = beliefs.getValue(c);
                if (b_c != s.zero()) {
                    double r_b_c = s.toReal(b_c);
                    double log_b_c = s.toLogProb(b_c);
                    sum += r_b_c * log_b_c;
                }
            }
            bethe -= (numNeighbors - 1) * sum;
        }
        
        assert !Double.isNaN(bethe);        
        return bethe;
    }
973571
def energy(data):
    """"""
    data = np.mean(data, axis=1)
    return np.sum(data ** 2) / np.float64(len(data))
138519
def energy(self):
        """"""
        result = 0.0
        for index1 in range(self.numc):
            for index2 in range(index1):
                if self.scaling[index1, index2] > 0:
                    for se, ve in self.yield_pair_energies(index1, index2):
                        result += se*ve*self.scaling[index1, index2]
        return result
1365777
def energy(self, hamiltonian=None):
        
        """

        if self.hamiltonian is None and hamiltonian is None:
            raise ValueError("To compute the total energy, a hamiltonian"
                             " object must be provided!")

        from ..potential import PotentialBase
        if isinstance(hamiltonian, PotentialBase):
            from ..potential import Hamiltonian

            warnings.warn("This function now expects a `Hamiltonian` instance "
                          "instead of a `PotentialBase` subclass instance. If "
                          "you are using a static reference frame, you just "
                          "need to pass your potential object in to the "
                          "Hamiltonian constructor to use, e.g., Hamiltonian"
                          "(potential).", DeprecationWarning)

            hamiltonian = Hamiltonian(hamiltonian)

        if hamiltonian is None:
            hamiltonian = self.hamiltonian

        return hamiltonian(self)
1078245
def energy_ES(q, v):
    """"""
    # Body 0 is the sun, Body 1 is the earth
    m0 = mass[0]
    m1 = mass[1]

    # Positions of sun and earth
    q0: np.ndarray = q[:, slices[0]]
    q1: np.ndarray = q[:, slices[1]]

    # Velocities of sun and earth
    v0: np.ndarray = v[:, slices[0]]
    v1: np.ndarray = v[:, slices[1]]

    # Kinetic energy is 1/2 mv^2
    T0: np.ndarray = 0.5 * m0 * np.sum(v0 * v0, axis=1)
    T1: np.ndarray = 0.5 * m1 * np.sum(v1 * v1, axis=1)
    T: np.ndarray = T0 + T1
    
    # Potential energy is -G m1 m2  / r
    dv_01 = q1 - q0
    r_01 = np.linalg.norm(dv_01, axis=1)
    U_01: np.ndarray = -G * m0 * m1 / r_01
    U: np.ndarray = U_01
    
    # Total energy H = T + U
    H = T + U
    
    return H, T, U
981774
def rms(self, x, params=()):
        
        internal_x, internal_params = self.pre_process(np.asarray(x),
                                                       np.asarray(params))
        if internal_params.ndim > 1:
            raise NotImplementedError("Parameters should be constant.")
        result = np.empty(internal_x.size//self.nx)
        for idx in range(internal_x.shape[0]):
            result[idx] = np.sqrt(np.mean(np.square(self.f_cb(
                internal_x[idx, :], internal_params))))
        return result
1148568
def get_XIJ(nat, i_now, dr_now, dr_old):
    
    # Do an element-wise outer product
    dr_dr = dr_now.reshape(-1,3,1)*dr_old.reshape(-1,1,3)

    xij = np.zeros([nat,3,3])
    for i in range(3):
        for j in range(3):
            # For each atom, sum over all neighbors
            xij[:,i,j] = np.bincount(i_now, weights=dr_dr[:,i,j])

    return xij
1298051
def energy_minimum(self, x):
    
    E      = self.energy
    y_min  = self.find_energy_minimum(E, x)
    return y_min
402427
def mse(x, xhat):
    
    buf_ = x - xhat
    np.square(buf_, out=buf_)  # square in-place
    sum_ = np.sum(buf_)
    sum_ /= x.size  # divide in-place
    return sum_
589794
def energy(self):
        
        return self.uncorrected_energy + np.sum(list(self.corrections.values()))