Python numpy.real() Examples

The following are 30 code examples of numpy.real(). You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example. You may also want to check out all available functions/classes of the module numpy , or try the search function .
Example #1
Source File: linalg_helper.py    From pyscf with Apache License 2.0 7 votes vote down vote up
def diagonalize_asymm(H):
    """
    Diagonalize a real, *asymmetric* matrix and return sorted results.

    Return the eigenvalues and eigenvectors (column matrix)
    sorted from lowest to highest eigenvalue.
    """
    E,C = np.linalg.eig(H)
    #if np.allclose(E.imag, 0*E.imag):
    #    E = np.real(E)
    #else:
    #    print "WARNING: Eigenvalues are complex, will be returned as such."

    idx = E.real.argsort()
    E = E[idx]
    C = C[:,idx]

    return E,C 
Example #2
Source File: grid.py    From simnibs with GNU General Public License v3.0 7 votes vote down vote up
def quadrature_cc_1D(N):
    """ Computes the Clenshaw Curtis nodes and weights """
    N = np.int(N)        
    if N == 1:
        knots = 0
        weights = 2
    else:
        n = N - 1
        C = np.zeros((N,2))
        k = 2*(1+np.arange(np.floor(n/2)))
        C[::2,0] = 2/np.hstack((1, 1-k*k))
        C[1,1] = -n
        V = np.vstack((C,np.flipud(C[1:n,:])))
        F = np.real(ifft(V, n=None, axis=0))
        knots = F[0:N,1]
        weights = np.hstack((F[0,0],2*F[1:n,0],F[n,0]))
            
    return knots, weights 
Example #3
Source File: test_xrft.py    From xrft with MIT License 6 votes vote down vote up
def test_dft_real_2d(self):
        """
        Test the real discrete Fourier transform function on one-dimensional
        data. Non-trivial because we need to keep only some of the negative
        frequencies.
        """
        Nx, Ny = 16, 32
        da = xr.DataArray(np.random.rand(Nx, Ny), dims=['x', 'y'],
                          coords={'x': range(Nx), 'y': range(Ny)})
        dx = float(da.x[1] - da.x[0])
        dy = float(da.y[1] - da.y[0])

        daft = xrft.dft(da, real='x')
        npt.assert_almost_equal(daft.values,
                               np.fft.rfftn(da.transpose('y','x')).transpose())
        npt.assert_almost_equal(daft.values,
                               xrft.dft(da, dim=['y'], real='x'))

        actual_freq_x = daft.coords['freq_x'].values
        expected_freq_x = np.fft.rfftfreq(Nx, dx)
        npt.assert_almost_equal(actual_freq_x, expected_freq_x)

        actual_freq_y = daft.coords['freq_y'].values
        expected_freq_y = np.fft.fftfreq(Ny, dy)
        npt.assert_almost_equal(actual_freq_y, expected_freq_y) 
Example #4
Source File: interpolation.py    From scarlet with MIT License 6 votes vote down vote up
def fft_convolve(*images):
    """Use FFT's to convove an image with a kernel

    Parameters
    ----------
    images: list of array-like
        A list of images to convolve.

    Returns
    -------
    result: array
        The convolution in pixel space of `img` with `kernel`.
    """
    from autograd.numpy.numpy_boxes import ArrayBox

    Images = [np.fft.fft2(np.fft.ifftshift(img)) for img in images]
    if np.any([isinstance(img, ArrayBox) for img in images]):
        Convolved = Images[0]
        for img in Images[1:]:
            Convolved = Convolved * img
    else:
        Convolved = np.prod(Images, 0)
    convolved = np.fft.ifft2(Convolved)
    return np.fft.fftshift(np.real(convolved)) 
Example #5
Source File: ewald.py    From pyqmc with MIT License 6 votes vote down vote up
def energy(self, configs):
        r"""
        Compute Coulomb energy for a set of configs.  

        .. math:: E_{\rm Coulomb} &= E_{\rm real+reciprocal}^{ee} 
                + E_{\rm self+charged}^{ee} 
                \\&+ E_{\rm real+reciprocal}^{e\text{-ion}} 
                + E_{\rm self+charged}^{e\text{-ion}} 
                \\&+ E_{\rm real+reciprocal}^{\text{ion-ion}} 
                + E_{\rm self+charged}^{\text{ion-ion}}
        
        Inputs:
            configs: pyqmc PeriodicConfigs object of shape (nconf, nelec, ndim)
        Returns: 
            ee: electron-electron part
            ei: electron-ion part
            ii: ion-ion part
        """
        nelec = configs.configs.shape[1]
        ee, ei = self.ewald_electron(configs)
        ee += self.ee_const(nelec)
        ei += self.ei_const(nelec)
        ii = self.ion_ion + self.ii_const
        return ee, ei, ii 
Example #6
Source File: linalg_helper.py    From pyscf with Apache License 2.0 6 votes vote down vote up
def __init__(self,matr_multiply,xStart,inPreCon,nroots=1,tol=1e-10):
        self.matrMultiply = matr_multiply
        self.size = xStart.shape[0]
        self.nEigen = min(nroots, self.size)
        self.maxM = min(30, self.size)
        self.maxOuterLoop = 10
        self.tol = tol

        #
        #  Creating initial guess and preconditioner
        #
        self.x0 = xStart.real.copy()

        self.iteration = 0
        self.totalIter = 0
        self.converged = False
        self.preCon = inPreCon.copy()
        #
        #  Allocating other vectors
        #
        self.allocateVecs() 
Example #7
Source File: lti.py    From python-control with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def damp(self):
        '''Natural frequency, damping ratio of system poles

        Returns
        -------
        wn : array
            Natural frequencies for each system pole
        zeta : array
            Damping ratio for each system pole
        poles : array
            Array of system poles
        '''
        poles = self.pole()

        if isdtime(self, strict=True):
            splane_poles = np.log(poles)/self.dt
        else:
            splane_poles = poles
        wn = absolute(splane_poles)
        Z = -real(splane_poles)/wn
        return wn, Z, poles 
Example #8
Source File: rlocus.py    From python-control with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def _break_points(num, den):
    """Extract break points over real axis and gains given these locations"""
    # type: (np.poly1d, np.poly1d) -> (np.array, np.array)
    dnum = num.deriv(m=1)
    dden = den.deriv(m=1)
    polynom = den * dnum - num * dden
    real_break_pts = polynom.r
    # don't care about infinite break points
    real_break_pts = real_break_pts[num(real_break_pts) != 0]
    k_break = -den(real_break_pts) / num(real_break_pts)
    idx = k_break >= 0   # only positives gains
    k_break = k_break[idx]
    real_break_pts = real_break_pts[idx]
    if len(k_break) == 0:
        k_break = [0]
        real_break_pts = den.roots
    return k_break, real_break_pts 
Example #9
Source File: linalg_helper.py    From pyscf with Apache License 2.0 6 votes vote down vote up
def solveSubspace(self):
        w, v = scipy.linalg.eig(self.subH[:self.currentSize,:self.currentSize])
        idx = w.real.argsort()
        v = v[:,idx]
        w = w[idx].real
#
        imag_norm = np.linalg.norm(w.imag)
        if imag_norm > 1e-12:
            print(" *************************************************** ")
            print(" WARNING  IMAGINARY EIGENVALUE OF NORM %.15g " % (imag_norm))
            print(" *************************************************** ")
        #print "Imaginary norm eigenvectors = ", np.linalg.norm(v.imag)
        #print "Imaginary norm eigenvalue   = ", np.linalg.norm(w.imag)
        #print "eigenvalues = ", w[:min(self.currentSize,7)]
#
        self.sol[:self.currentSize] = v[:,self.ciEig]
        self.evecs[:self.currentSize,:self.currentSize] = v
        self.eigs[:self.currentSize] = w[:self.currentSize]
        self.outeigs[:self.nEigen] = w[:self.nEigen]
        self.cvEig = self.eigs[self.ciEig] 
Example #10
Source File: steering-gainsched.py    From python-control with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def control_output(t, x, u, params):
    # Get the controller parameters
    longpole = params.get('longpole', -2.)
    latpole1 = params.get('latpole1', -1/2 + sqrt(-7)/2)
    latpole2 = params.get('latpole2', -1/2 - sqrt(-7)/2)
    l = params.get('wheelbase', 3)
    
    # Extract the system inputs
    ex, ey, etheta, vd, phid = u

    # Determine the controller gains
    alpha1 = -np.real(latpole1 + latpole2)
    alpha2 = np.real(latpole1 * latpole2)

    # Compute and return the control law
    v = -longpole * ex          # Note: no feedfwd (to make plot interesting)
    if vd != 0:
        phi = phid + (alpha1 * l) / vd * ey + (alpha2 * l) / vd * etheta
    else:
        # We aren't moving, so don't turn the steering wheel
        phi = phid
    
    return  np.array([v, phi])

# Define the controller as an input/output system 
Example #11
Source File: steering-gainsched.py    From python-control with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def control_output(t, x, u, params):
    # Get the controller parameters
    longpole = params.get('longpole', -2.)
    latpole1 = params.get('latpole1', -1/2 + sqrt(-7)/2)
    latpole2 = params.get('latpole2', -1/2 - sqrt(-7)/2)
    l = params.get('wheelbase', 3)
    
    # Extract the system inputs
    ex, ey, etheta, vd, phid = u

    # Determine the controller gains
    alpha1 = -np.real(latpole1 + latpole2)
    alpha2 = np.real(latpole1 * latpole2)

    # Compute and return the control law
    v = -longpole * ex          # Note: no feedfwd (to make plot interesting)
    if vd != 0:
        phi = phid + (alpha1 * l) / vd * ey + (alpha2 * l) / vd * etheta
    else:
        # We aren't moving, so don't turn the steering wheel
        phi = phid
    
    return  np.array([v, phi])

# Define the controller as an input/output system 
Example #12
Source File: EasyTL.py    From transferlearning with MIT License 6 votes vote down vote up
def get_cosine_dist(A, B):
    B = np.reshape(B, (1, -1))
    
    if A.shape[1] == 1:
        A = np.hstack((A, np.zeros((A.shape[0], 1))))
        B = np.hstack((B, np.zeros((B.shape[0], 1))))
    
    aa = np.sum(np.multiply(A, A), axis=1).reshape(-1, 1)
    bb = np.sum(np.multiply(B, B), axis=1).reshape(-1, 1)
    ab = A @ B.T
    
    # to avoid NaN for zero norm
    aa[aa==0] = 1
    bb[bb==0] = 1
    
    D = np.real(np.ones((A.shape[0], B.shape[0])) - np.multiply((1/np.sqrt(np.kron(aa, bb.T))), ab))
    
    return D 
Example #13
Source File: polynomial.py    From pyGSTi with Apache License 2.0 6 votes vote down vote up
def compact(self, complex_coeff_tape=True):
        """
        Generate a compact form of this polynomial designed for fast evaluation.

        The resulting "tapes" can be evaluated using
        :function:`opcalc.bulk_eval_compact_polys`.

        Parameters
        ----------
        complex_coeff_tape : bool, optional
            Whether the `ctape` returned array is forced to be of complex type.
            If False, the real part of all coefficients is taken (even if they're
            complex).

        Returns
        -------
        vtape, ctape : numpy.ndarray
            These two 1D arrays specify an efficient means for evaluating this
            polynomial.
        """
        if complex_coeff_tape:
            return self._rep.compact_complex()
        else:
            return self._rep.compact_real() 
Example #14
Source File: __init__.py    From pyGSTi with Apache License 2.0 6 votes vote down vote up
def safe_bulk_eval_compact_polys(vtape, ctape, paramvec, dest_shape):
    """Typechecking wrapper for :function:`bulk_eval_compact_polys`.

    The underlying method has two implementations: one for real-valued
    `ctape`, and one for complex-valued. This wrapper will dynamically
    dispatch to the appropriate implementation method based on the
    type of `ctape`. If the type of `ctape` is known prior to calling,
    it's slightly faster to call the appropriate implementation method
    directly; if not.
    """
    if _np.iscomplexobj(ctape):
        ret = bulk_eval_compact_polys_complex(vtape, ctape, paramvec, dest_shape)
        im_norm = _np.linalg.norm(_np.imag(ret))
        if im_norm > 1e-6:
            print("WARNING: norm(Im part) = {:g}".format(im_norm))
    else:
        ret = bulk_eval_compact_polys(vtape, ctape, paramvec, dest_shape)
    return _np.real(ret) 
Example #15
Source File: reportableqty.py    From pyGSTi with Apache License 2.0 6 votes vote down vote up
def absdiff(self, constant_value, separate_re_im=False):
        """
        Returns a ReportableQty that is the (element-wise in the vector case)
        difference between `constant_value` and this one given by:

        `abs(self - constant_value)`.
        """
        if separate_re_im:
            re_v = _np.fabs(_np.real(self.value) - _np.real(constant_value))
            im_v = _np.fabs(_np.imag(self.value) - _np.imag(constant_value))
            if self.has_eb():
                return (ReportableQty(re_v, _np.fabs(_np.real(self.errbar)), self.nonMarkovianEBs),
                        ReportableQty(im_v, _np.fabs(_np.imag(self.errbar)), self.nonMarkovianEBs))
            else:
                return ReportableQty(re_v), ReportableQty(im_v)

        else:
            v = _np.absolute(self.value - constant_value)
            if self.has_eb():
                return ReportableQty(v, _np.absolute(self.errbar), self.nonMarkovianEBs)
            else:
                return ReportableQty(v) 
Example #16
Source File: reportableqty.py    From pyGSTi with Apache License 2.0 6 votes vote down vote up
def infidelity_diff(self, constant_value):
        """
        Returns a ReportableQty that is the (element-wise in the vector case)
        difference between `constant_value` and this one given by:

        `1.0 - Re(conjugate(constant_value) * self )`
        """
        # let diff(x) = 1.0 - Re(const.C * x) = 1.0 - (const.re * x.re + const.im * x.im)
        # so d(diff)/dx.re = -const.re, d(diff)/dx.im = -const.im
        # diff(x + dx) = diff(x) + d(diff)/dx * dx
        # diff(x + dx) - diff(x) =  - (const.re * dx.re + const.im * dx.im)
        v = 1.0 - _np.real(_np.conjugate(constant_value) * self.value)
        if self.has_eb():
            eb = abs(_np.real(constant_value) * _np.real(self.errbar)
                     + _np.imag(constant_value) * _np.real(self.errbar))
            return ReportableQty(v, eb, self.nonMarkovianEBs)
        else:
            return ReportableQty(v) 
Example #17
Source File: circuit.py    From pyGSTi with Apache License 2.0 6 votes vote down vote up
def _num_to_rqc_str(num):
    """Convert float to string to be included in RQC quil DEFGATE block
    (as written by _np_to_quil_def_str)."""
    num = _np.complex(_np.real_if_close(num))
    if _np.imag(num) == 0:
        output = str(_np.real(num))
        return output
    else:
        real_part = _np.real(num)
        imag_part = _np.imag(num)
        if imag_part < 0:
            sgn = '-'
            imag_part = imag_part * -1
        elif imag_part > 0:
            sgn = '+'
        else:
            assert False
        return '{}{}{}i'.format(real_part, sgn, imag_part) 
Example #18
Source File: point_cloud.py    From FRIDA with MIT License 6 votes vote down vote up
def classical_mds(self, D):
        ''' 
        Classical multidimensional scaling

        Parameters
        ----------
        D : square 2D ndarray
            Euclidean Distance Matrix (matrix containing squared distances between points
        '''

        # Apply MDS algorithm for denoising
        n = D.shape[0]
        J = np.eye(n) - np.ones((n,n))/float(n)
        G = -0.5*np.dot(J, np.dot(D, J))

        s, U = np.linalg.eig(G)

        # we need to sort the eigenvalues in decreasing order
        s = np.real(s)
        o = np.argsort(s)
        s = s[o[::-1]]
        U = U[:,o[::-1]]

        S = np.diag(s)[0:self.dim,:]
        self.X = np.dot(np.sqrt(S),U.T) 
Example #19
Source File: slowreplib.py    From pyGSTi with Apache License 2.0 5 votes vote down vote up
def probability(self, state):
        scratch = _np.empty(self.dim, 'd')
        Edense = self.todense(scratch)
        return _np.dot(Edense, state.base)  # not vdot b/c data is *real* 
Example #20
Source File: ewald.py    From pyqmc with MIT License 5 votes vote down vote up
def set_lattice_displacements(self, nlatvec):
        """
        Generates list of lattice-vector displacements to add together for real-space sum, going from `-nlatvec` to `nlatvec` in each lattice direction.
        """
        XYZ = np.meshgrid(*[np.arange(-nlatvec, nlatvec + 1)] * 3, indexing="ij")
        xyz = np.stack(XYZ, axis=-1).reshape((-1, 3))
        self.lattice_displacements = np.dot(xyz, self.latvec) 
Example #21
Source File: slowreplib.py    From pyGSTi with Apache License 2.0 5 votes vote down vote up
def probability(self, state):
        # can assume state is a DMStateRep
        return _np.dot(self.base, state.base)  # not vdot b/c *real* data 
Example #22
Source File: reportableqty.py    From pyGSTi with Apache License 2.0 5 votes vote down vote up
def hermitian_to_real(self):
        """
        Returns a ReportableQty that holds the real matrix
        whose upper/lower triangle contains the real/imaginary parts
        of the corresponding off-diagonal matrix elements of the
        *Hermitian* matrix stored in this ReportableQty.

        This is used for display purposes.  If this object doesn't
        contain a Hermitian matrix, `ValueError` is raised.
        """
        if _np.linalg.norm(self.value - _np.conjugate(self.value).T) > 1e-8:
            raise ValueError("Contained value must be Hermitian!")

        def _convert(A):
            ret = _np.empty(A.shape, 'd')
            for i in range(A.shape[0]):
                ret[i, i] = A[i, i].real
                for j in range(i + 1, A.shape[1]):
                    ret[i, j] = A[i, j].real
                    ret[j, i] = A[i, j].imag
            return ret

        v = _convert(self.value)
        if self.has_eb():
            eb = _convert(self.errbar)
            return ReportableQty(v, eb, self.nonMarkovianEBs)
        else:
            return ReportableQty(v) 
Example #23
Source File: ecg_simulate.py    From NeuroKit with MIT License 5 votes vote down vote up
def _ecg_simulate_rrprocess(
    flo=0.1, fhi=0.25, flostd=0.01, fhistd=0.01, lfhfratio=0.5, hrmean=60, hrstd=1, sfrr=1, n=256
):
    w1 = 2 * np.pi * flo
    w2 = 2 * np.pi * fhi
    c1 = 2 * np.pi * flostd
    c2 = 2 * np.pi * fhistd
    sig2 = 1
    sig1 = lfhfratio
    rrmean = 60 / hrmean
    rrstd = 60 * hrstd / (hrmean * hrmean)

    df = sfrr / n
    w = np.arange(n) * 2 * np.pi * df
    dw1 = w - w1
    dw2 = w - w2

    Hw1 = sig1 * np.exp(-0.5 * (dw1 / c1) ** 2) / np.sqrt(2 * np.pi * c1 ** 2)
    Hw2 = sig2 * np.exp(-0.5 * (dw2 / c2) ** 2) / np.sqrt(2 * np.pi * c2 ** 2)
    Hw = Hw1 + Hw2
    Hw0 = np.concatenate((Hw[0 : int(n / 2)], Hw[int(n / 2) - 1 :: -1]))
    Sw = (sfrr / 2) * np.sqrt(Hw0)

    ph0 = 2 * np.pi * np.random.uniform(size=int(n / 2 - 1))
    ph = np.concatenate([[0], ph0, [0], -np.flipud(ph0)])
    SwC = Sw * np.exp(1j * ph)
    x = (1 / n) * np.real(np.fft.ifft(SwC))

    xstd = np.std(x)
    ratio = rrstd / xstd
    return rrmean + x * ratio  # Return RR 
Example #24
Source File: slowreplib.py    From pyGSTi with Apache License 2.0 5 votes vote down vote up
def compact_real(self):
        """
        Returns a real representation of this polynomial as a
        `(variable_tape, coefficient_tape)` 2-tuple of 1D nupy arrays.
        The coefficient tape is *always* a complex array, even if
        none of the polynomial's coefficients are complex.

        Such compact representations are useful for storage and later
        evaluation, but not suited to polynomial manipulation.

        Returns
        -------
        vtape : numpy.ndarray
            A 1D array of integers (variable indices).
        ctape : numpy.ndarray
            A 1D array of *real* coefficients.
        """
        nTerms = len(self)
        vinds = {i: self._int_to_vinds(i) for i in self.keys()}
        nVarIndices = sum(map(len, vinds.values()))
        vtape = _np.empty(1 + nTerms + nVarIndices, _np.int64)  # "variable" tape
        ctape = _np.empty(nTerms, complex)  # "coefficient tape"

        i = 0
        vtape[i] = nTerms; i += 1
        for iTerm, k in enumerate(sorted(self.keys())):
            v = vinds[k]  # so don't need to compute self._int_to_vinds(k)
            l = len(v)
            ctape[iTerm] = self[k]
            vtape[i] = l; i += 1
            vtape[i:i + l] = v; i += l
        assert(i == len(vtape)), "Logic Error!"
        return vtape, ctape 
Example #25
Source File: signal_psd.py    From NeuroKit with MIT License 5 votes vote down vote up
def _signal_psd_burg(signal, sampling_rate=1000, order=15, criteria="KIC", corrected=True, side="one-sided", norm=True, nperseg=None):

    nfft = int(nperseg * 2)
    ar, rho, ref = _signal_arma_burg(signal, order=order, criteria=criteria, corrected=corrected, side=side, norm=norm)
    psd = _signal_psd_from_arma(ar=ar, rho=rho, sampling_rate=sampling_rate, nfft=nfft, side=side, norm=norm)

    # signal is real, not complex
    if nfft % 2 == 0:
        power = psd[0:int(nfft / 2 + 1)] * 2
    else:
        power = psd[0:int((nfft + 1) / 2)] * 2

    # angular frequencies, w
    # for one-sided psd, w spans [0, pi]
    # for two-sdied psd, w spans [0, 2pi)
    # for dc-centered psd, w spans (-pi, pi] for even nfft, (-pi, pi) for add nfft
    if side == "one-sided":
        w = np.pi * np.linspace(0, 1, len(power))
#    elif side == "two-sided":
#        w = np.pi * np.linspace(0, 2, len(power), endpoint=False)  #exclude last point
#    elif side == "centerdc":
#        if nfft % 2 == 0:
#            w = np.pi * np.linspace(-1, 1, len(power))
#        else:
#            w = np.pi * np.linspace(-1, 1, len(power) + 1, endpoint=False)  # exclude last point
#            w = w[1:]  # exclude first point (extra)

    frequency = (w * sampling_rate) / (2 * np.pi)

    return frequency, power 
Example #26
Source File: cond.py    From simnibs with GNU General Public License v3.0 5 votes vote down vote up
def _get_sorted_eigenv(tensors):
    eig_val, eig_vec = np.linalg.eig(tensors)
    eig_val = np.real(eig_val)
    eig_vec = np.real(eig_vec)
    sort = eig_val.argsort(axis=1)[:, ::-1]
    eig_val = np.sort(eig_val, axis=1)[:, ::-1]
    s_eig_vec = np.zeros_like(tensors)
    for i in range(3):
        s_eig_vec[:, :, i] = eig_vec[np.arange(len(eig_vec)), :, sort[:, i]]
    eig_vec = s_eig_vec
    return eig_val, eig_vec 
Example #27
Source File: energy.py    From pyqmc with MIT License 5 votes vote down vote up
def kinetic(configs, wf):
    nconf, nelec, ndim = configs.configs.shape
    ke = np.zeros(nconf)
    for e in range(nelec):
        ke += -0.5 * np.real(wf.laplacian(e, configs.electron(e)))
    return ke 
Example #28
Source File: ewald.py    From pyqmc with MIT License 5 votes vote down vote up
def ewald_ion(self):
        r"""
        Compute ion contribution to Ewald sums.  Since the ions don't move in our calculations, the ion-ion term only needs to be computed once.

        Note: We ignore the constant term :math:`\frac{1}{2} \sum_{I} Z_I^2 C_{\rm self\ image}` in the real-space ion-ion sum corresponding to the interaction of an ion with its own image in other cells.

        The real-space part:

        .. math:: E_{\rm real\ space}^{\text{ion-ion}} = \sum_{\vec{n}} \sum_{I<J}^{N_{ion}} Z_I Z_J \frac{{\rm erfc}(\alpha |\vec{x}_{IJ}+\vec{n}|)}{|\vec{x}_{IJ}+\vec{n}|} 

        The reciprocal-space part:

        .. math:: E_{\rm reciprocal\ space}^{\text{ion-ion}} = \sum_{\vec{G} > 0 } W_G \left| \sum_{I=1}^{N_{ion}} Z_I e^{-i\vec{G}\cdot\vec{x}_I} \right|^2

        Returns:
            ion_ion: float, ion-ion component of Ewald sum
        """
        # Real space part
        if len(self.atom_charges) == 1:
            ion_ion_real = 0
        else:
            dist = pyqmc.distance.MinimalImageDistance(self.latvec)
            ion_distances, ion_inds = dist.dist_matrix(self.atom_coords[np.newaxis])
            rvec = ion_distances[:, :, np.newaxis, :] + self.lattice_displacements
            r = np.linalg.norm(rvec, axis=-1)
            charge_ij = np.prod(self.atom_charges[np.asarray(ion_inds)], axis=1)
            ion_ion_real = np.einsum("j,ijk->", charge_ij, erfc(self.alpha * r) / r)

        # Reciprocal space part
        GdotR = np.dot(self.gpoints, self.atom_coords.T)
        self.ion_exp = np.dot(np.exp(1j * GdotR), self.atom_charges)
        ion_ion_rec = np.dot(self.gweight, np.abs(self.ion_exp) ** 2)

        ion_ion = ion_ion_real + ion_ion_rec
        return ion_ion 
Example #29
Source File: ewald.py    From pyqmc with MIT License 5 votes vote down vote up
def set_up_reciprocal_ewald_sum(self, ewald_gmax):
        r"""
        Determine parameters for Ewald sums. 

        :math:`\alpha` determines the partitioning of the real and reciprocal-space parts.

        We define a weight `gweight` for the part of the reciprocal-space sums that doesn't depend on the coordinates:
        
        .. math:: W_G = \frac{4\pi}{V |\vec{G}|^2} e^{- \frac{|\vec{G}|^2}{ 4\alpha^2}}

        Inputs:
            latvec: (3, 3) array of lattice vectors; latvec[0] is the first
            ewald_gmax: int, max number of reciprocal lattice vectors to check away from 0
        """
        cellvolume = np.linalg.det(self.latvec)
        recvec = np.linalg.inv(self.latvec)
        crossproduct = recvec.T * cellvolume

        # Determine alpha
        tmpheight_i = np.einsum("ij,ij->i", crossproduct, self.latvec)
        length_i = np.linalg.norm(crossproduct, axis=1)
        smallestheight = np.amin(np.abs(tmpheight_i) / length_i)
        self.alpha = 5.0 / smallestheight
        print("Setting Ewald alpha to ", self.alpha)

        # Determine G points to include in reciprocal Ewald sum
        XYZ = np.meshgrid(*[np.arange(-ewald_gmax, ewald_gmax + 1)] * 3, indexing="ij")
        X, Y, Z = [x.ravel() for x in XYZ]
        positive_octants = X + 1e-6 * Y + 1e-12 * Z > 0  # assume ewald_gmax < 1e5
        gpoints = np.stack((X, Y, Z), axis=-1)[positive_octants]
        gpoints = np.dot(gpoints, recvec) * 2 * np.pi
        gsquared = np.sum(gpoints ** 2, axis=1)
        gweight = 4 * np.pi * np.exp(-gsquared / (4 * self.alpha ** 2))
        gweight /= cellvolume * gsquared
        bigweight = gweight > 1e-10
        self.gpoints = gpoints[bigweight]
        self.gweight = gweight[bigweight]

        self.set_ewald_constants(cellvolume) 
Example #30
Source File: optimize_ortho.py    From pyqmc with MIT License 5 votes vote down vote up
def evaluate(wfs, coords, pgrad, sampler, sample_options, warmup):
    return_data, coords = sampler(wfs, coords, pgrad, **sample_options)
    """ 
    For wave functions wfs and coordinate set coords, evaluate the overlap and energy of the last wave function. 

    Returns a dictionary with relevant information.
    """
    avg_data = {}
    for k, it in return_data.items():
        avg_data[k] = np.average(it[warmup:, ...], axis=0)
    N = avg_data["overlap"].diagonal()
    # Derivatives are only for the optimized wave function, so they miss
    # an index
    N_derivative = 2 * np.real(avg_data["overlap_gradient"][-1])
    Nij = np.outer(N, N)
    S = avg_data["overlap"] / np.sqrt(Nij)
    S_derivative = avg_data["overlap_gradient"] / Nij[-1, :, np.newaxis] - np.einsum(
        "j,m->jm", avg_data["overlap"][-1, :] / Nij[-1, :], N_derivative / N[-1]
    )
    energy_derivative = 2.0 * (avg_data["dpH"] - avg_data["total"] * avg_data["dppsi"])
    dp = avg_data["dppsi"]
    condition = np.real(avg_data["dpidpj"] - np.einsum("i,j->ij", dp, dp))

    return {
        "N": N,
        "S": S,
        "S_derivative": S_derivative,
        "energy_derivative": energy_derivative,
        "N_derivative": N_derivative,
        "condition": condition,
        "total": avg_data["total"],
    }