Python numpy.conj() Examples

The following are 30 code examples of numpy.conj(). 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: m_c2r.py    From pyscf with Apache License 2.0 6 votes vote down vote up
def __init__(self, j):
    self._j = j
    self._c2r = np.zeros( (2*j+1, 2*j+1), dtype=np.complex128)
    self._c2r[j,j]=1.0
    for m in range(1,j+1):
      self._c2r[m+j, m+j] = sgn[m] * np.sqrt(0.5) 
      self._c2r[m+j,-m+j] = np.sqrt(0.5) 
      self._c2r[-m+j,-m+j]= 1j*np.sqrt(0.5)
      self._c2r[-m+j, m+j]= -sgn[m] * 1j * np.sqrt(0.5)
    
    self._hc_c2r = np.conj(self._c2r).transpose()
    self._conj_c2r = np.conjugate(self._c2r) # what is the difference ? conj and conjugate
    self._tr_c2r = np.transpose(self._c2r)
    
    #print(abs(self._hc_c2r.conj().transpose()-self._c2r).sum())

  #
  #
  # 
Example #2
Source File: optimal_givens_decomposition_test.py    From OpenFermion-Cirq with Apache License 2.0 6 votes vote down vote up
def test_givens_inverse():
    r"""
    The Givens rotation in OpenFermion is defined as

    .. math::

        \begin{pmatrix}
            \cos(\theta) & -e^{i \varphi} \sin(\theta) \\
            \sin(\theta) &     e^{i \varphi} \cos(\theta)
        \end{pmatrix}.

    confirm numerically its hermitian conjugate is it's inverse
    """
    a = numpy.random.random() + 1j * numpy.random.random()
    b = numpy.random.random() + 1j * numpy.random.random()
    ab_rotation = givens_matrix_elements(a, b, which='right')

    assert numpy.allclose(ab_rotation.dot(numpy.conj(ab_rotation).T),
                          numpy.eye(2))
    assert numpy.allclose(numpy.conj(ab_rotation).T.dot(ab_rotation),
                          numpy.eye(2)) 
Example #3
Source File: optimal_givens_decomposition_test.py    From OpenFermion-Cirq with Apache License 2.0 6 votes vote down vote up
def test_circuit_generation_and_accuracy():
    for dim in range(2, 10):
        qubits = cirq.LineQubit.range(dim)
        u_generator = numpy.random.random(
            (dim, dim)) + 1j * numpy.random.random((dim, dim))
        u_generator = u_generator - numpy.conj(u_generator).T
        assert numpy.allclose(-1 * u_generator, numpy.conj(u_generator).T)

        unitary = scipy.linalg.expm(u_generator)
        circuit = cirq.Circuit()
        circuit.append(optimal_givens_decomposition(qubits, unitary))

        fermion_generator = QubitOperator(()) * 0.0
        for i, j in product(range(dim), repeat=2):
            fermion_generator += jordan_wigner(
                FermionOperator(((i, 1), (j, 0)), u_generator[i, j]))

        true_unitary = scipy.linalg.expm(
            get_sparse_operator(fermion_generator).toarray())
        assert numpy.allclose(true_unitary.conj().T.dot(true_unitary),
                              numpy.eye(2 ** dim, dtype=complex))

        test_unitary = cirq.unitary(circuit)
        assert numpy.isclose(
            abs(numpy.trace(true_unitary.conj().T.dot(test_unitary))), 2 ** dim) 
Example #4
Source File: matrix.py    From vampyre with MIT License 6 votes vote down vote up
def svd_dotH(self,s1,q1):
        """
        Performs diagonal matrix multiplication conjugate
        
        Implements :math:`q_0 = \\mathrm{diag}(s_1)^* q_1`.
        
        :param s1: diagonal parameters
        :param q1: input to the diagonal multiplication
        :returns: :code:`q0` diagonal multiplication output
        """
        srep = repeat_axes(np.conj(s1),self.sshape,self.srep_axes,rep=False)
        q0 = srep*q1
        return q0
        
    #def __str__(self):
    #    string = str(self.name) + '\n'\
    #              + 'Input: ' + str(self.shape0) + ',' + str(self.dtype0)
    #    return string 
Example #5
Source File: fourier.py    From ibllib with MIT License 6 votes vote down vote up
def fexpand(x, ns=1, axis=None):
    """
    Reconstructs full spectrum from positive frequencies
    Works on the last dimension (contiguous in c-stored array)

    :param x: numpy.ndarray
    :param axis: axis along which to perform reduction (last axis by default)
    :return: numpy.ndarray
    """
    if axis is None:
        axis = x.ndim - 1
    # dec = int(ns % 2) * 2 - 1
    # xcomp = np.conj(np.flip(x[..., 1:x.shape[-1] + dec], axis=axis))
    ilast = int((ns + (ns % 2)) / 2)
    xcomp = np.conj(np.flip(np.take(x, np.arange(1, ilast), axis=axis), axis=axis))
    return np.concatenate((x, xcomp), axis=axis) 
Example #6
Source File: wishart.py    From SAR-change-detection with MIT License 6 votes vote down vote up
def azimuthal_symmetry(X, Y, n, m):
    p1 = 2
    p2 = 1
    p = np.sqrt(p1**2 + p2**2)

    detX = np.real(X.hvhv*(X.hhhh*X.vvvv - X.hhvv*np.conj(X.hhvv)))
    detY = np.real(Y.hvhv*(Y.hhhh*Y.vvvv - Y.hhvv*np.conj(Y.hhvv)))
    detXY = np.real((X.hvhv+Y.hvhv) * ((X.hhhh+Y.hhhh)*(X.vvvv+Y.vvvv) - (X.hhvv+Y.hhvv)*(np.conj(X.hhvv)+np.conj(Y.hhvv))))

    lnq = (p*(n+m)*np.log(n+m) - p*n*np.log(n) - p*m*np.log(m)
            + n*np.log(detX) + m*np.log(detY) - (n+m)*np.log(detXY))

    rho1 = 1 - (2*p1**2 - 1)/(6*p1) * (1/n + 1/m - 1/(n+m))
    rho2 = 1 - (2*p2**2 - 1)/(6*p2) * (1/n + 1/m - 1/(n+m))
    rho = 1/p**2 * (p1**2 * rho1 + p2**2 * rho2)

    w2 = - p**2/4 * (1-1/rho)**2 + (p1**2*(p1**2-1) + p2**2*(p2**2-1))/24 * (1/n**2 + 1/m**2 - 1/(n+m)**2) * 1/rho**2

    return lnq, rho, w2 
Example #7
Source File: kmp2.py    From pyscf with Apache License 2.0 6 votes vote down vote up
def _gamma1_intermediates(mp, t2=None):
    # Memory optimization should be here
    if t2 is None:
        t2 = mp.t2
    if t2 is None:
        raise NotImplementedError("Run kmp2.kernel with `with_t2=True`")
    nmo = mp.nmo
    nocc = mp.nocc
    nvir = nmo - nocc
    nkpts = mp.nkpts
    dtype = t2.dtype

    dm1occ = np.zeros((nkpts, nocc, nocc), dtype=dtype)
    dm1vir = np.zeros((nkpts, nvir, nvir), dtype=dtype)

    for ki in range(nkpts):
        for kj in range(nkpts):
            for ka in range(nkpts):
                kb = mp.khelper.kconserv[ki, ka, kj]

                dm1vir[kb] += einsum('ijax,ijay->yx', t2[ki][kj][ka].conj(), t2[ki][kj][ka]) * 2 -\
                              einsum('ijax,ijya->yx', t2[ki][kj][ka].conj(), t2[ki][kj][kb])
                dm1occ[kj] += einsum('ixab,iyab->xy', t2[ki][kj][ka].conj(), t2[ki][kj][ka]) * 2 -\
                              einsum('ixab,iyba->xy', t2[ki][kj][ka].conj(), t2[ki][kj][kb])
    return -dm1occ, dm1vir 
Example #8
Source File: optools.py    From pyGSTi with Apache License 2.0 6 votes vote down vote up
def unitary_to_process_mx(U):
    """
    Compute the super-operator which acts on (row)-vectorized
    density matrices from a unitary operator (matrix) U which
    acts on state vectors.  This super-operator is given by
    the tensor product of U and conjugate(U), i.e. kron(U,U.conj).

    Parameters
    ----------
    U : numpy array
        The unitary matrix which acts on state vectors.

    Returns
    -------
    numpy array
       The super-operator process matrix.
    """
    # U -> kron(U,Uc) since U rho U_dag -> kron(U,Uc)
    #  since AXB --row-vectorize--> kron(A,B.T)*vec(X)
    return _np.kron(U, _np.conjugate(U)) 
Example #9
Source File: test_fft.py    From pyscf with Apache License 2.0 6 votes vote down vote up
def get_j(cell, dm, hermi=1, vhfopt=None, kpt=np.zeros(3), kpts_band=None):
    dm = np.asarray(dm)
    nao = dm.shape[-1]

    coords = gen_grid.gen_uniform_grids(cell)
    if kpts_band is None:
        kpt1 = kpt2 = kpt
        aoR_k1 = aoR_k2 = numint.eval_ao(cell, coords, kpt)
    else:
        kpt1 = kpts_band
        kpt2 = kpt
        aoR_k1 = numint.eval_ao(cell, coords, kpt1)
        aoR_k2 = numint.eval_ao(cell, coords, kpt2)
    ngrids, nao = aoR_k1.shape

    def contract(dm):
        vjR_k2 = get_vjR(cell, dm, aoR_k2)
        vj = (cell.vol/ngrids) * np.dot(aoR_k1.T.conj(), vjR_k2.reshape(-1,1)*aoR_k1)
        return vj

    if dm.ndim == 2:
        vj = contract(dm)
    else:
        vj = lib.asarray([contract(x) for x in dm.reshape(-1,nao,nao)])
    return vj.reshape(dm.shape) 
Example #10
Source File: test_fft.py    From pyscf with Apache License 2.0 6 votes vote down vote up
def get_vkR(mf, cell, aoR_k1, aoR_k2, kpt1, kpt2):
    '''Get the real-space 2-index "exchange" potential V_{i,k1; j,k2}(r)
    where {i,k1} = exp^{i k1 r) |i> , {j,k2} = exp^{-i k2 r) <j|
    '''
    coords = gen_grid.gen_uniform_grids(cell)
    ngrids, nao = aoR_k1.shape

    expmikr = np.exp(-1j*np.dot(kpt1-kpt2,coords.T))
    coulG = tools.get_coulG(cell, kpt1-kpt2, exx=True, mf=mf)
    def prod(ij):
        i, j = divmod(ij, nao)
        rhoR = aoR_k1[:,i] * aoR_k2[:,j].conj()
        rhoG = tools.fftk(rhoR, cell.mesh, expmikr)
        vG = coulG*rhoG
        vR = tools.ifftk(vG, cell.mesh, expmikr.conj())
        return vR

    if aoR_k1.dtype == np.double and aoR_k2.dtype == np.double:
        vR = numpy.asarray([prod(ij).real for ij in range(nao**2)])
    else:
        vR = numpy.asarray([prod(ij) for ij in range(nao**2)])
    return vR.reshape(nao,nao,-1).transpose(2,0,1) 
Example #11
Source File: tdscf.py    From pyscf with Apache License 2.0 6 votes vote down vote up
def get_k(self,dm):
        """
        Update Exact Exchange Matrix

        Args:
            dm: float or complex
                AO density matrix.
        Returns:
            kmat: float or complex
                Exact Exchange in AO basis
        """
        naux = self.auxmol.nao_nr()
        nao = self.ks.mol.nao_nr()
        kpj = np.einsum("ijp,jk->ikp", self.eri3c, dm)
        pik = np.linalg.solve(self.eri2c, kpj.reshape(-1,naux).T.conj())
        kmat = np.einsum("pik,kjp->ij", pik.reshape(naux,nao,nao), self.eri3c)
        return kmat 
Example #12
Source File: scalarnl.py    From vampyre with MIT License 5 votes vote down vote up
def fnl_aug(self,z,r,rvar):
        """
        Augmented nonlinear function and its gradient.
        
        Given the penalty function :math:`f(z)`, the function returns
        the value and gradient of the augmented function,
        
        :math:`f_{aug}(z,r) = f(z) + (1/\\tau_r)|r-z|^2            
        """
        
        # Evaluate function
        if self.hess_avail:        
            f0, fgrad0, fhess0 = self.fnl(z)
        else:
            f0, fgrad0 = self.fnl(z)
                    
        # Add the augmenting term
        rvar_rep = common.repeat_axes(rvar,self.shape,self.var_axes,rep=False)
        aug = np.abs(z-r)**2/rvar_rep
        aug_grad = 2*np.conj(z-r)/rvar_rep
        aug_hess = 2/rvar_rep
        if not self.is_complex:
            aug /= 2
            aug_grad /= 2
            aug_hess /= 2
        faug = np.sum(f0 + aug)
        faug_grad = fgrad0 + aug_grad
        if self.hess_avail:
            faug_hess = fhess0 + aug_hess                         
            return faug, faug_grad, faug_hess
        else:
            return faug, faug_grad 
Example #13
Source File: test_umath_complex.py    From recruit with Apache License 2.0 5 votes vote down vote up
def test_simple_conjugate(self):
        ref = np.conj(np.sqrt(complex(1, 1)))

        def f(z):
            return np.sqrt(np.conj(z))

        check_complex_value(f, 1, 1, ref.real, ref.imag, False)

    #def test_branch_cut(self):
    #    _check_branch_cut(f, -1, 0, 1, -1) 
Example #14
Source File: test_umath_complex.py    From recruit with Apache License 2.0 5 votes vote down vote up
def test_cabs_inf_nan(self):
        x, y = [], []

        # cabs(+-nan + nani) returns nan
        x.append(np.nan)
        y.append(np.nan)
        check_real_value(np.abs,  np.nan, np.nan, np.nan)

        x.append(np.nan)
        y.append(-np.nan)
        check_real_value(np.abs, -np.nan, np.nan, np.nan)

        # According to C99 standard, if exactly one of the real/part is inf and
        # the other nan, then cabs should return inf
        x.append(np.inf)
        y.append(np.nan)
        check_real_value(np.abs,  np.inf, np.nan, np.inf)

        x.append(-np.inf)
        y.append(np.nan)
        check_real_value(np.abs, -np.inf, np.nan, np.inf)

        # cabs(conj(z)) == conj(cabs(z)) (= cabs(z))
        def f(a):
            return np.abs(np.conj(a))

        def g(a, b):
            return np.abs(complex(a, b))

        xa = np.array(x, dtype=complex)
        for i in range(len(xa)):
            ref = g(x[i], y[i])
            check_real_value(f, x[i], y[i], ref) 
Example #15
Source File: sar_data.py    From SAR-change-detection with MIT License 5 votes vote down vote up
def determinant(self):
        "Determinants of the covariance matrices in a SARData object"
        return np.real((self.hhhh*self.hvhv*self.vvvv
            + self.hhhv*self.hvvv*np.conj(self.hhvv)
            + self.hhvv*np.conj(self.hhhv)*np.conj(self.hvvv)
            - self.hhvv*self.hvhv*np.conj(self.hhvv)
            - self.hhhv*np.conj(self.hhhv)*self.vvvv
            - self.hhhh*self.hvvv*np.conj(self.hvvv))) 
Example #16
Source File: emd.py    From magpy with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def calc_inst_info(modes,samplerate):
    """
Calculate the instantaneous frequency, amplitude, and phase of
each mode.
    """

    amp=np.zeros(modes.shape,np.float32)
    phase=np.zeros(modes.shape,np.float32)
    f=np.zeros(modes.shape,np.float32)

    print("Mode 1:", len(modes), samplerate)

    for m in range(len(modes)):
        h=scipy.signal.hilbert(modes[m])
        print(len(modes[m]))
        print("Mean Amplitude of mode ", m, np.mean(np.abs(h)))
        print("Mean Phase of mode ", m, np.mean(np.angle(h)))
        phase[m,:]=np.angle(h)
        print("Frequ", np.diff(np.unwrap(phase[:,np.r_[0,0:len(modes[m])]]))/(2*np.pi)*samplerate)
        amp[m,:]=np.abs(h)
        phase[m,:]=np.angle(h)
        f[m,:] = np.r_[np.nan,
                      0.5*(np.angle(-h[2:]*np.conj(h[0:-2]))+np.pi)/(2*np.pi) * samplerate,
                      np.nan]
        print("Mean Frequ of mode ", m, np.mean(np.diff(np.unwrap(phase[:,np.r_[0,0:len(modes[0])]]))/(2*np.pi)*samplerate))

        #f(m,:) = [nan 0.5*(angle(-h(t+1).*conj(h(t-1)))+pi)/(2*pi) * sr nan];

    # calc the freqs (old way)
    #f=np.diff(np.unwrap(phase[:,np.r_[0,0:len(modes[0])]]))/(2*np.pi)*samplerate

    # clip the freqs so they don't go below zero
    #f = f.clip(0,f.max())

    return f,amp,phase 
Example #17
Source File: unitary_svd.py    From vampyre with MIT License 5 votes vote down vote up
def svd_dotH(self,s1,q1):
        """
        Performs diagonal matrix multiplication conjugate

        Implements :math:`q_0 = \\mathrm{diag}(s_1)^* q_1`.

        :param s1: diagonal parameters
        :param q1: input to the diagonal multiplication
        :returns: :code:`q0` diagonal multiplication output
        """
        srep = common.repeat_axes(np.conj(s1),self.shape0,(0,1),rep=False)
        q0 = srep*q1
        return q0 
Example #18
Source File: haldane_C3.py    From tenpy with GNU General Public License v3.0 5 votes vote down vote up
def init_terms(self, model_params):
        t = np.asarray(model_params.get('t', -1.))
        V = np.asarray(model_params.get('V', 0))
        phi_ext = 2 * np.pi * model_params.get('phi_ext', 0.)

        t1 = t
        t2 = 0.39 * t * 1j
        t3 = -0.34 * t

        for u1, u2, dx in self.lat.NN:
            t1_phi = self.coupling_strength_add_ext_flux(t1, dx, [0, phi_ext])
            self.add_coupling(t1_phi, u1, 'CdA', u2, 'CB', dx, 'JW', True)
            self.add_coupling(np.conj(t1_phi), u2, 'CdB', u1, 'CA', -dx, 'JW', True)
            self.add_coupling(V, u1, 'Ntot', u2, 'Ntot', dx)

        for u1, u2, dx in self.lat.nNNA:
            t2_phi = self.coupling_strength_add_ext_flux(t2, dx, [0, phi_ext])
            self.add_coupling(t2_phi, u1, 'CdA', u2, 'CA', dx, 'JW', True)
            self.add_coupling(np.conj(t2_phi), u2, 'CdA', u1, 'CA', -dx, 'JW', True)

        for u1, u2, dx in self.lat.nNNB:
            t2_phi = self.coupling_strength_add_ext_flux(t2, dx, [0, phi_ext])
            self.add_coupling(t2_phi, u1, 'CdB', u2, 'CB', dx, 'JW', True)
            self.add_coupling(np.conj(t2_phi), u2, 'CdB', u1, 'CB', -dx, 'JW', True)

        for u1, u2, dx in self.lat.nnNN:
            t3_phi = self.coupling_strength_add_ext_flux(t3, dx, [0, phi_ext])
            self.add_coupling(t3_phi, u1, 'CdA', u2, 'CB', dx, 'JW', True)
            self.add_coupling(np.conj(t3_phi), u2, 'CdB', u1, 'CA', -dx, 'JW', True) 
Example #19
Source File: qutrit.py    From pyGSTi with Apache License 2.0 5 votes vote down vote up
def _random_rot(scale, arrType=_np.array, seed=None):
    rndm = _np.random.RandomState(seed)
    randH = scale * (rndm.randn(3, 3) + 1j * rndm.randn(3, 3))
    randH = _np.dot(_np.conj(randH.T), randH)
    randU = _linalg.expm(-1j * randH)
    return arrType(randU) 
Example #20
Source File: tfvis.py    From python-control with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def mouse_move(self, event):
        """ Handle mouse movement, redraw pzmap while drag/dropping """
        if (self.move_zero != None and
            event.xdata != None and 
            event.ydata != None):

            if (self.index1 == self.index2):
                # Real pole/zero
                new = event.xdata
                if (self.move_zero == True):
                    self.zeros[self.index1] = new
                elif (self.move_zero == False):
                    self.poles[self.index1] = new
            else:
                # Complex poles/zeros
                new = event.xdata + 1.0j*event.ydata
                if (self.move_zero == True):
                    self.zeros[self.index1] = new
                    self.zeros[self.index2] = conj(new)
                elif (self.move_zero == False):
                    self.poles[self.index1] = new
                    self.poles[self.index2] = conj(new)
            tfcn = None
            if (self.move_zero == True):
                self.tfi.set_zeros(self.zeros)
                tfcn = self.tfi.get_tf()
            elif (self.move_zero == False):
                self.tfi.set_poles(self.poles)
                tfcn = self.tfi.get_tf()
            if (tfcn != None):
                self.draw_pz(tfcn)
                self.canvas_pzmap.draw() 
Example #21
Source File: mps.py    From tenpy with GNU General Public License v3.0 5 votes vote down vote up
def entanglement_spectrum(self, by_charge=False):
        r"""return entanglement energy spectrum.

        Parameters
        ----------
        by_charge : bool
            Wheter we should sort the spectrum on each bond by the possible charges.

        Returns
        -------
        ent_spectrum : list
            For each (non-trivial) bond the entanglement spectrum.
            If `by_charge` is ``False``, return (for each bond) a sorted 1D ndarray
            with the convention :math:`S_i^2 = e^{-\xi_i}`, where :math:`S_i` labels a Schmidt
            value and :math:`\xi_i` labels the entanglement 'energy' in the returned spectrum.
            If `by_charge` is True, return a a list of tuples ``(charge, sub_spectrum)``
            for each possible charge on that bond.
        """
        if by_charge:
            res = []
            for i in range(self.L + 1)[self.nontrivial_bonds]:
                ss = -2. * np.log(self._S[i])
                if i < self.L:
                    leg = self._B[i].get_leg('vL')
                else:  # i == L: segment b.c.
                    leg = self._B[i - 1].get_leg('vR').conj()
                spectrum = [(leg.get_charge(qi), np.sort(ss[leg.get_slice(qi)]))
                            for qi in range(leg.block_number)]
                res.append(spectrum)
            return res
        else:
            return [np.sort(-2. * np.log(ss)) for ss in self._S[self.nontrivial_bonds]] 
Example #22
Source File: tfvis.py    From python-control with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def button_press(self, event):
        """ Handle button presses, detect if we are going to move
        any poles/zeros"""
        # find closest pole/zero
        if event.xdata != None and event.ydata != None:

            new = event.xdata + 1.0j*event.ydata

            tzeros = list(abs(self.zeros-new))
            tpoles = list(abs(self.poles-new))

            if (size(tzeros) > 0):
                minz = min(tzeros)
            else:
                minz = float('inf')
            if (size(tpoles) > 0):
                minp = min(tpoles)
            else:
                minp = float('inf')

            if (minz < 2 or minp < 2):
                if (minz < minp):
                    # Moving zero(s)
                    self.index1 = tzeros.index(minz)
                    self.index2 = list(self.zeros).index(
                        conj(self.zeros[self.index1]))
                    self.move_zero = True
                else:
                    # Moving pole(s)
                    self.index1 = tpoles.index(minp)
                    self.index2 = list(self.poles).index(
                        conj(self.poles[self.index1]))
                    self.move_zero = False 
Example #23
Source File: mps.py    From tenpy with GNU General Public License v3.0 5 votes vote down vote up
def probability_per_charge(self, bond=0):
        """Return probabilites of charge value on the left of a given bond.

        For example for particle number conservation, define
        :math:`N_b = sum_{i<b} n_i` for a given bond `b`.
        This function returns the possible values of `N_b` as rows of `charge_values`,
        and for each row the probabilty that this combination occurs in the given state.

        Parameters
        ----------
        bond : int
            The bond to be considered. The returned charges are summed on the left of this bond.

        Returns
        -------
        charge_values : 2D array
            Columns correspond to the different charges in `self.chinfo`.
            Rows are the different charge fluctuations at this bond
        probabilities : 1D array
            For each row of `charge_values` the probablity for these values of charge fluctuations.
        """
        if self.bc == 'segment' and bond == self.L:
            S = self.get_SR(self.L - 1)**2
            leg = self.get_B(self.L - 1, form=None).get_leg('vR').conj()
        else:  # usually the case
            S = self.get_SL(bond)**2
            leg = self.get_B(bond, form=None).get_leg('vL')
        assert leg.qconj == +1
        if not leg.is_blocked():
            raise ValueError("leg not blocked: can have duplicate entries in charge values")
        ps = []
        for qi in range(leg.block_number):
            sl = leg.get_slice(qi)
            ps.append(np.sum(S[sl]))
        ps = np.array(ps)
        if abs(np.sum(ps) - 1.) > 1.e-10:
            warnings.warn("Probability_per_charge: Sum of probabilites not 1. Canonical form?",
                          stacklevel=2)
        return leg.charges.copy(), ps 
Example #24
Source File: mps.py    From tenpy with GNU General Public License v3.0 5 votes vote down vote up
def full_contraction(self, i0):
        """Calculate the overlap by a full contraction of the network.

        The full contraction of the environments gives the overlap ``<bra|ket>``,
        taking into account :attr:`MPS.norm` of both `bra` and `ket`.
        For this purpose, this function contracts
        ``get_LP(i0+1, store=False)`` and ``get_RP(i0, store=False)`` with appropriate singular
        values in between.

        Parameters
        ----------
        i0 : int
            Site index.
        """
        if self.ket.finite and i0 + 1 == self.L:
            # special case to handle `_to_valid_index` correctly:
            # get_LP(L) is not valid for finite b.c, so we use need to calculate it explicitly.
            LP = self.get_LP(i0, store=False)
            LP = self._contract_LP(i0, LP)
        else:
            LP = self.get_LP(i0 + 1, store=False)
        # multiply with `S`: a bit of a hack: use 'private' MPS._scale_axis_B
        S_bra = self.bra.get_SR(i0).conj()
        LP = self.bra._scale_axis_B(LP, S_bra, form_diff=1., axis_B='vR*', cutoff=0.)
        # cutoff is not used for form_diff = 1
        S_ket = self.ket.get_SR(i0)
        LP = self.bra._scale_axis_B(LP, S_ket, form_diff=1., axis_B='vR', cutoff=0.)
        RP = self.get_RP(i0, store=False)
        contr = npc.inner(LP, RP, axes=[['vR*', 'vR'], ['vL*', 'vL']], do_conj=False)
        return contr * self.bra.norm * self.ket.norm 
Example #25
Source File: tdscf.py    From pyscf with Apache License 2.0 5 votes vote down vote up
def fockbuild(self,dm_lao,it = -1):
        """
        Updates Fock matrix

        Args:
            dm_lao: float or complex
                Lowdin AO density matrix.
            it: int
                iterator for SCF DIIS

        Returns:
            fmat: float or complex
                Fock matrix in Lowdin AO basis
            jmat: float or complex
                Coulomb matrix in AO basis
            kmat: float or complex
                Exact Exchange in AO basis
        """
        if self.params["Model"] == "TDHF":
            Pt = 2.0*transmat(dm_lao,self.x,-1)
            jmat,kmat = self.get_jk(Pt)
            veff = 0.5*(jmat+jmat.T.conj()) - 0.5*(0.5*(kmat + kmat.T.conj()))
            if self.adiis and it > 0:
                return transmat(self.adiis.update(self.s,Pt,self.h + veff),\
                self.x), jmat, kmat
            else:
                return  transmat(self.h + veff,self.x), jmat, kmat
        elif self.params["Model"] == "TDDFT":
            Pt = 2 * transmat(dm_lao,self.x,-1)
            jmat = self.J = self.get_j(Pt)
            Veff = self.J.astype(complex)
            Vxc, excsum, kmat = self.get_vxc(Pt)
            Veff += Vxc
            if self.adiis and it > 0:
                return transmat(self.adiis.update(self.s,Pt,self.h + Veff),\
                self.x), jmat, kmat
            else:
                return transmat(self.h + Veff,self.x), jmat, kmat 
Example #26
Source File: tdscf.py    From pyscf with Apache License 2.0 5 votes vote down vote up
def transmat(M,U,inv = 1):
    if inv == 1:
        # U.t() * M * U
        Mtilde = np.dot(np.dot(U.T.conj(),M),U)
    elif inv == -1:
        # U * M * U.t()
        Mtilde = np.dot(np.dot(U,M),U.T.conj())
    return Mtilde 
Example #27
Source File: m_lorentzian.py    From pyscf with Apache License 2.0 5 votes vote down vote up
def llc_imag(x, w1, w2, e1, e2):
  """ Product of two Lorentzians one of which is conjugated """
  res = (lorentzian(x, w1, e1)*np.conj(lorentzian(x, w2, e2))).imag
  return res 
Example #28
Source File: xrft.py    From xrft with MIT License 5 votes vote down vote up
def _power_spectrum(daft, dim, N, density):

    ps = (daft * np.conj(daft)).real

    if density:
        ps /= (np.asarray(N).prod()) ** 2
        for i in dim:
            ps /= daft['freq_' + i + '_spacing']

    return ps 
Example #29
Source File: mps.py    From tenpy with GNU General Public License v3.0 5 votes vote down vote up
def _contract_LP(self, i, LP):
        """Contract LP with the tensors on site `i` to form ``self._LP[i+1]``"""
        LP = npc.tensordot(LP, self.ket.get_B(i, form='A'), axes=('vR', 'vL'))
        axes = (self.ket._get_p_label('*') + ['vL*'], self.ket._p_label + ['vR*'])
        # for a ususal MPS, axes = (['p*', 'vL*'], ['p', 'vR*'])
        LP = npc.tensordot(self.bra.get_B(i, form='A').conj(), LP, axes=axes)
        return LP  # labels 'vR*', 'vR' 
Example #30
Source File: matrix.py    From vampyre with MIT License 5 votes vote down vote up
def dotH(self,z1):
        """
        Compute conjugate transpose multiplication:math:`A^*(z1)`
        """
        return self.A.conj().T.dot(z1)