Python numpy.conj() Examples

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())

  #
  #
  # 

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)) 

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) 

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 

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) 

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 

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 

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)) 

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) 

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) 

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 

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 

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) 

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) 

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))) 

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 

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 

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) 

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) 

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() 

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]] 

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 

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 

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 

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 

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 

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 

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 

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' 

def dotH(self,z1):
        """
        Compute conjugate transpose multiplication:math:`A^*(z1)`
        """
        return self.A.conj().T.dot(z1)