Python numpy.complex() Examples

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 

def get_init_guess(self, kshift, nroots=1, koopmans=False, diag=None, **kwargs):
        """Initial guess vectors of R coefficients"""
        size = self.vector_size(kshift)
        dtype = getattr(diag, 'dtype', np.complex)
        nroots = min(nroots, size)
        guess = []
        # TODO do Koopmans later
        if koopmans:
            raise NotImplementedError
        else:
            idx = diag.argsort()[:nroots]
            for i in idx:
                g = np.zeros(int(size), dtype=dtype)
                g[i] = 1.0
                # TODO do mask_frozen later
                guess.append(g)
        return guess 

def make_s10(mol, gauge_orig=None, mb='RMB'):
    '''First order overlap matrix wrt external magnetic field.
    Note the side effects of set_common_origin.
    '''
    if gauge_orig is not None:
        mol.set_common_origin(gauge_orig)
    n2c = mol.nao_2c()
    n4c = n2c * 2
    c = lib.param.LIGHT_SPEED
    s1 = numpy.zeros((3, n4c, n4c), complex)
    if mb.upper() == 'RMB':
        if gauge_orig is None:
            t1 = mol.intor('int1e_giao_sa10sp_spinor', 3)
        else:
            t1 = mol.intor('int1e_cg_sa10sp_spinor', 3)
        for i in range(3):
            t1cc = t1[i] + t1[i].conj().T
            s1[i,n2c:,n2c:] = t1cc * (.25/c**2)

    if gauge_orig is None:
        sg = mol.intor('int1e_govlp_spinor', 3)
        tg = mol.intor('int1e_spgsp_spinor', 3)
        s1[:,:n2c,:n2c] += sg
        s1[:,n2c:,n2c:] += tg * (.25/c**2)
    return s1 

def cc_Wvvvv(t1,t2,eris):
    tau = make_tau(t2,t1,t1)
    #eris_vovv = np.array(eris.ovvv).transpose(1,0,3,2)
    #tmp = einsum('mb,amef->abef',t1,eris_vovv)
    #Wabef = eris.vvvv - tmp + tmp.transpose(1,0,2,3)
    #Wabef += 0.25*einsum('mnab,mnef->abef',tau,eris.oovv)
    if t1.dtype == np.complex: ds_type = 'c16'
    else: ds_type = 'f8'
    _tmpfile1 = tempfile.NamedTemporaryFile(dir=lib.param.TMPDIR)
    fimd = h5py.File(_tmpfile1.name)
    nocc, nvir = t1.shape
    Wabef = fimd.create_dataset('vvvv', (nvir,nvir,nvir,nvir), ds_type)
    for a in range(nvir):
        Wabef[a] = eris.vvvv[a] 
        Wabef[a] -= einsum('mb,mfe->bef',t1,eris.ovvv[:,a,:,:]) 
        Wabef[a] += einsum('m,mbfe->bef',t1[:,a],eris.ovvv) 
        Wabef[a] += 0.25*einsum('mnb,mnef->bef',tau[:,:,a,:],eris.oovv)
    return Wabef 

def Wvvvv(t1,t2,eris):
    tau = make_tau(t2,t1,t1)
    #Wabef = cc_Wvvvv(t1,t2,eris) + 0.25*einsum('mnab,mnef->abef',tau,eris.oovv)
    if t1.dtype == np.complex: ds_type = 'c16'
    else: ds_type = 'f8'
    _tmpfile1 = tempfile.NamedTemporaryFile(dir=lib.param.TMPDIR)
    fimd = h5py.File(_tmpfile1.name)
    nocc, nvir = t1.shape
    Wabef = fimd.create_dataset('vvvv', (nvir,nvir,nvir,nvir), ds_type)
    #_cc_Wvvvv = cc_Wvvvv(t1,t2,eris)
    for a in range(nvir):
        #Wabef[a] = _cc_Wvvvv[a]
        Wabef[a] = eris.vvvv[a] 
        Wabef[a] -= einsum('mb,mfe->bef',t1,eris.ovvv[:,a,:,:]) 
        Wabef[a] += einsum('m,mbfe->bef',t1[:,a],eris.ovvv) 
        #Wabef[a] += 0.25*einsum('mnb,mnef->bef',tau[:,:,a,:],eris.oovv)

        #Wabef[a] += 0.25*einsum('mnb,mnef->bef',tau[:,:,a,:],eris.oovv) 
        Wabef[a] += 0.5*einsum('mnb,mnef->bef',tau[:,:,a,:],eris.oovv) 
    return Wabef 

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) 

def allocateVecs(self):
        self.subH = np.zeros( shape=(self.maxM,self.maxM), dtype=complex )
        self.sol = np.zeros( shape=(self.maxM), dtype=complex )
        self.dgks = np.zeros( shape=(self.maxM), dtype=complex )
        self.nConv = np.zeros( shape=(self.maxM), dtype=int )
        self.eigs = np.zeros( shape=(self.maxM), dtype=complex )
        self.evecs = np.zeros( shape=(self.maxM,self.maxM), dtype=complex )
        self.oldeigs = np.zeros( shape=(self.maxM), dtype=complex )
        self.deigs = np.zeros( shape=(self.maxM), dtype=complex )
        self.outeigs = np.zeros( shape=(self.nEigen), dtype=complex )
        self.outevecs = np.zeros( shape=(self.size,self.nEigen), dtype=complex)
        self.currentSize = 0

        self.Ax = np.zeros( shape=(self.size), dtype=complex )
        self.res = np.zeros( shape=(self.size), dtype=complex )
        self.vlist = np.zeros( shape=(self.maxM,self.size), dtype=complex )
        self.cv = np.zeros( shape = (self.size), dtype = complex )
        self.cAv = np.zeros( shape = (self.size), dtype = complex )
        self.Avlist = np.zeros( shape=(self.maxM,self.size), dtype=complex )
        self.dres = 999.9
        self.resnorm = 999.9
        self.cvEig = 0.1
        self.ciEig = 0
        self.deflated = 0 

def test_better_float():
    # Better float function
    def check_against(f1, f2):
        return f1 if FLOAT_TYPES.index(f1) >= FLOAT_TYPES.index(f2) else f2
    for first in FLOAT_TYPES:
        for other in IUINT_TYPES + np.sctypes['complex']:
            assert_equal(better_float_of(first, other), first)
            assert_equal(better_float_of(other, first), first)
            for other2 in IUINT_TYPES + np.sctypes['complex']:
                assert_equal(better_float_of(other, other2), np.float32)
                assert_equal(better_float_of(other, other2, np.float64),
                             np.float64)
        for second in FLOAT_TYPES:
            assert_equal(better_float_of(first, second),
                         check_against(first, second))
    # Check codes and dtypes work
    assert_equal(better_float_of('f4', 'f8', 'f4'), np.float64)
    assert_equal(better_float_of('i4', 'i8', 'f8'), np.float64) 

def get_init_guess(self, kshift, nroots=1, koopmans=False, diag=None):
        size = self.vector_size()
        dtype = getattr(diag, 'dtype', np.complex)
        nroots = min(nroots, size)
        guess = []
        if koopmans:
            for n in self.nonzero_opadding[kshift][::-1][:nroots]:
                g = np.zeros(int(size), dtype=dtype)
                g[n] = 1.0
                g = self.mask_frozen(g, kshift, const=0.0)
                guess.append(g)
        else:
            idx = diag.argsort()[:nroots]
            for i in idx:
                g = np.zeros(int(size), dtype=dtype)
                g[i] = 1.0
                g = self.mask_frozen(g, kshift, const=0.0)
                guess.append(g)
        return guess 

def hcore_generator(self, mol):
        aoslices = mol.aoslice_2c_by_atom()
        h1 = self.get_hcore(mol)
        n2c = mol.nao_2c()
        n4c = n2c * 2
        c = lib.param.LIGHT_SPEED
        def hcore_deriv(atm_id):
            shl0, shl1, p0, p1 = aoslices[atm_id]
            with mol.with_rinv_at_nucleus(atm_id):
                z = -mol.atom_charge(atm_id)
                vn = z * mol.intor('int1e_iprinv_spinor', comp=3)
                wn = z * mol.intor('int1e_ipsprinvsp_spinor', comp=3)

            v = numpy.zeros((3,n4c,n4c), numpy.complex)
            v[:,:n2c,:n2c] = vn
            v[:,n2c:,n2c:] = wn * (.25/c**2)
            v[:,p0:p1]         += h1[:,p0:p1]
            v[:,n2c+p0:n2c+p1] += h1[:,n2c+p0:n2c+p1]
            return v + v.conj().transpose(0,2,1)
        return hcore_deriv 

def test_intor_r2e(self):
        mol1 = gto.M(atom=[[1   , (0. , -0.7 , 0.)],
                           [1   , (0. ,  0.7 , 0.)]],
                     basis = '631g')
        nao = mol1.nao_2c()
        eri0 = numpy.empty((3,nao,nao,nao,nao), dtype=numpy.complex)
        ip = 0
        for i in range(mol1.nbas):
            jp = 0
            for j in range(mol1.nbas):
                kp = 0
                for k in range(mol1.nbas):
                    lp = 0
                    for l in range(mol1.nbas):
                        buf = mol1.intor_by_shell('int2e_ip1_spinor', (i,j,k,l), comp=3)
                        di,dj,dk,dl = buf.shape[1:]
                        eri0[:,ip:ip+di,jp:jp+dj,kp:kp+dk,lp:lp+dl] = buf
                        lp += dl
                    kp += dk
                jp += dj
            ip += di 

def _patch_cython():
    # "monkey patch" some objects to avoid cyclic import structure
    from . import _npc_helper
    _npc_helper._charges = charges
    _npc_helper._np_conserved = np_conserved
    assert _npc_helper.QTYPE == charges.QTYPE
    # check types
    warn = False
    import numpy as np
    check_types = [
        (np.float64, np.complex128),
        (np.ones([1]).dtype, (1.j * np.ones([1])).dtype),
        (np.array(1.).dtype, np.array(1.j).dtype),
        (np.array(1., dtype=np.float).dtype, np.array(1., dtype=np.complex).dtype),
    ]
    types_ok = [
        _npc_helper._float_complex_are_64_bit(dt_float, dt_real)
        for dt_float, dt_real in check_types
    ]
    if not np.all(types_ok):
        import warnings
        warnings.warn("(Some of) the default dtypes are not 64-bit. "
                      "Using the compiled cython code (as you do) might make it slower.")
    # done 

def para(mol, mo10, mo_coeff, mo_occ, shielding_nuc=None):
    if shielding_nuc is None:
        shielding_nuc = range(mol.natm)
    n4c = mo_coeff.shape[1]
    n2c = n4c // 2
    msc_para = numpy.zeros((len(shielding_nuc),3,3))
    para_neg = numpy.zeros((len(shielding_nuc),3,3))
    para_occ = numpy.zeros((len(shielding_nuc),3,3))
    h01 = numpy.zeros((3, n4c, n4c), complex)
    orbo = mo_coeff[:,mo_occ>0]
    for n, atm_id in enumerate(shielding_nuc):
        mol.set_rinv_origin(mol.atom_coord(atm_id))
        t01 = mol.intor('int1e_sa01sp_spinor', 3)
        for m in range(3):
            h01[m,:n2c,n2c:] = .5 * t01[m]
            h01[m,n2c:,:n2c] = .5 * t01[m].conj().T
        h01_mo = [reduce(numpy.dot, (mo_coeff.conj().T, x, orbo)) for x in h01]
        for b in range(3):
            for m in range(3):
                # + c.c.
                p = numpy.einsum('ij,ij->i', mo10[b].conj(), h01_mo[m]).real * 2
                msc_para[n,b,m] = p.sum()
                para_neg[n,b,m] = p[:n2c].sum()
                para_occ[n,b,m] = p[mo_occ>0].sum()
    para_pos = msc_para - para_neg - para_occ
    return msc_para, para_pos, para_neg, para_occ 

def test_byte_orders():
    arr = np.arange(10, dtype=np.int32)
    # Test endian read/write of types not requiring scaling
    for tp in (np.uint64, np.float, np.complex):
        dt = np.dtype(tp)
        for code in '<>':
            ndt = dt.newbyteorder(code)
            for klass in (SlopeInterArrayWriter, SlopeArrayWriter,
                          ArrayWriter):
                aw = klass(arr, ndt)
                data_back = round_trip(aw)
                assert_array_almost_equal(arr, data_back) 

def test_isposinf(self):
        a = nx.array([nx.inf, -nx.inf, nx.nan, 0.0, 3.0, -3.0])
        out = nx.zeros(a.shape, bool)
        tgt = nx.array([True, False, False, False, False, False])

        res = ufl.isposinf(a)
        assert_equal(res, tgt)
        res = ufl.isposinf(a, out)
        assert_equal(res, tgt)
        assert_equal(out, tgt)

        a = a.astype(np.complex)
        with assert_raises(TypeError):
            ufl.isposinf(a) 

def test_setitem_ndarray_1d(self):
        # GH5508

        # len of indexer vs length of the 1d ndarray
        df = DataFrame(index=Index(lrange(1, 11)))
        df['foo'] = np.zeros(10, dtype=np.float64)
        df['bar'] = np.zeros(10, dtype=np.complex)

        # invalid
        with pytest.raises(ValueError):
            df.loc[df.index[2:5], 'bar'] = np.array([2.33j, 1.23 + 0.1j,
                                                     2.2, 1.0])

        # valid
        df.loc[df.index[2:6], 'bar'] = np.array([2.33j, 1.23 + 0.1j,
                                                 2.2, 1.0])

        result = df.loc[df.index[2:6], 'bar']
        expected = Series([2.33j, 1.23 + 0.1j, 2.2, 1.0], index=[3, 4, 5, 6],
                          name='bar')
        tm.assert_series_equal(result, expected)

        # dtype getting changed?
        df = DataFrame(index=Index(lrange(1, 11)))
        df['foo'] = np.zeros(10, dtype=np.float64)
        df['bar'] = np.zeros(10, dtype=np.complex)

        with pytest.raises(ValueError):
            df[2:5] = np.arange(1, 4) * 1j 

def test_isneginf(self):
        a = nx.array([nx.inf, -nx.inf, nx.nan, 0.0, 3.0, -3.0])
        out = nx.zeros(a.shape, bool)
        tgt = nx.array([False, True, False, False, False, False])

        res = ufl.isneginf(a)
        assert_equal(res, tgt)
        res = ufl.isneginf(a, out)
        assert_equal(res, tgt)
        assert_equal(out, tgt)

        a = a.astype(np.complex)
        with assert_raises(TypeError):
            ufl.isneginf(a) 

def test_is_complex_dtype():
    assert not com.is_complex_dtype(int)
    assert not com.is_complex_dtype(str)
    assert not com.is_complex_dtype(pd.Series([1, 2]))
    assert not com.is_complex_dtype(np.array(['a', 'b']))

    assert com.is_complex_dtype(np.complex)
    assert com.is_complex_dtype(np.array([1 + 1j, 5])) 

def test_arraywriters():
    # Test initialize
    # Simple cases
    if machine() == 'sparc64' and python_compiler().startswith('GCC'):
        # bus errors on at least np 1.4.1 through 1.6.1 for complex
        test_types = FLOAT_TYPES + IUINT_TYPES
    else:
        test_types = NUMERIC_TYPES
    for klass in (SlopeInterArrayWriter, SlopeArrayWriter, ArrayWriter):
        for type in test_types:
            arr = np.arange(10, dtype=type)
            aw = klass(arr)
            assert_true(aw.array is arr)
            assert_equal(aw.out_dtype, arr.dtype)
            assert_array_equal(arr, round_trip(aw))
            # Byteswapped is OK
            bs_arr = arr.byteswap().newbyteorder('S')
            bs_aw = klass(bs_arr)
            # assert against original array because POWER7 was running into
            # trouble using the byteswapped array (bs_arr)
            assert_array_equal(arr, round_trip(bs_aw))
            bs_aw2 = klass(bs_arr, arr.dtype)
            assert_array_equal(arr, round_trip(bs_aw2))
            # 2D array
            arr2 = np.reshape(arr, (2, 5))
            a2w = klass(arr2)
            # Default out - in order is Fortran
            arr_back = round_trip(a2w)
            assert_array_equal(arr2, arr_back)
            arr_back = round_trip(a2w, 'F')
            assert_array_equal(arr2, arr_back)
            # C order works as well
            arr_back = round_trip(a2w, 'C')
            assert_array_equal(arr2, arr_back)
            assert_true(arr_back.flags.c_contiguous) 

def __init__(self, zmean, zvar, shape,name=None,\
                 var_axes = (0,), zmean_axes='all',\
                 is_complex=False, map_est=False, tune_zvar=False,\
                 tune_rvar=False):
                
        if np.isscalar(shape):
            shape = (shape,)
        if is_complex:
            dtype = np.double
        else:
            dtype = np.complex
             
        BaseEst.__init__(self,shape=shape,dtype=dtype,name=name,
                         var_axes=var_axes,type_name='GaussEst', cost_avail=True)
        self.zmean = zmean
        self.zvar = zvar
        self.cost_avail = True  
        self.is_complex = is_complex  
        self.map_est = map_est         
        self.zmean_axes = zmean_axes
        self.tune_zvar = tune_zvar
        self.tune_rvar = tune_rvar
        
        ndim = len(self.shape)
        if self.zmean_axes == 'all':
            self.zmean_axes = tuple(range(ndim))
            
        # If zvar is a scalar, then repeat it to the required shape,
        # which are all the dimensions not being averaged over
        if np.isscalar(self.zvar):
            var_shape = get_var_shape(self.shape, self.var_axes)
            self.zvar = np.tile(self.zvar, var_shape) 

def _call_giao_vhf1(mol, dm):
    c1 = .5 / lib.param.LIGHT_SPEED
    n2c = dm.shape[0] // 2
    dmll = dm[:n2c,:n2c].copy()
    dmls = dm[:n2c,n2c:].copy()
    dmsl = dm[n2c:,:n2c].copy()
    dmss = dm[n2c:,n2c:].copy()
    vj = numpy.zeros((3,n2c*2,n2c*2), dtype=numpy.complex)
    vk = numpy.zeros((3,n2c*2,n2c*2), dtype=numpy.complex)
    vx = _vhf.rdirect_mapdm('int2e_g1_spinor', 'a4ij',
                            ('lk->s2ij', 'jk->s1il'), dmll, 3,
                            mol._atm, mol._bas, mol._env)
    vj[:,:n2c,:n2c] = vx[0]
    vk[:,:n2c,:n2c] = vx[1]
    vx = _vhf.rdirect_mapdm('int2e_spgsp1spsp2_spinor', 'a4ij',
                            ('lk->s2ij', 'jk->s1il'), dmss, 3,
                            mol._atm, mol._bas, mol._env) * c1**4
    vj[:,n2c:,n2c:] = vx[0]
    vk[:,n2c:,n2c:] = vx[1]
    vx = _vhf.rdirect_bindm('int2e_g1spsp2_spinor', 'a4ij',
                            ('lk->s2ij', 'jk->s1il'), (dmss,dmls), 3,
                            mol._atm, mol._bas, mol._env) * c1**2
    vj[:,:n2c,:n2c] += vx[0]
    vk[:,:n2c,n2c:] += vx[1]
    vx = _vhf.rdirect_bindm('int2e_spgsp1_spinor', 'a4ij',
                            ('lk->s2ij', 'jk->s1il'), (dmll,dmsl), 3,
                            mol._atm, mol._bas, mol._env) * c1**2
    vj[:,n2c:,n2c:] += vx[0]
    vk[:,n2c:,:n2c] += vx[1]
    for i in range(3):
        vj[i] = lib.hermi_triu(vj[i], 1)
        vk[i] = vk[i] + vk[i].T.conj()
    return vj, vk 

def make_h10rkb(mol, dm0, gauge_orig=None, with_gaunt=False,
                verbose=logger.WARN):
    '''Note the side effects of set_common_origin'''

    if gauge_orig is not None:
        mol.set_common_origin(gauge_orig)
    if isinstance(verbose, logger.Logger):
        log = verbose
    else:
        log = logger.Logger(mol.stdout, verbose)
    log.debug('first order Fock matrix / RKB')
    if gauge_orig is None:
        t1 = mol.intor('int1e_giao_sa10sp_spinor', 3)
    else:
        t1 = mol.intor('int1e_cg_sa10sp_spinor', 3)
    if with_gaunt:
        sys.stderr('NMR gaunt part not implemented\n')
    n2c = t1.shape[2]
    n4c = n2c * 2
    h1 = numpy.zeros((3, n4c, n4c), complex)
    for i in range(3):
        h1[i,:n2c,n2c:] += .5 * t1[i]
        h1[i,n2c:,:n2c] += .5 * t1[i].conj().T
    return h1

#TODO the uncouupled force 

def calc_bulk_twirled_DDD(model, germsList, eps=1e-6, check=False,
                          germLengths=None, comm=None):
    """Calculate the positive squares of the germ Jacobians.
    twirledDerivDaggerDeriv == array J.H*J contributions from each germ
    (J=Jacobian) indexed by (iGerm, iModelParam1, iModelParam2)
    size (nGerms, vec_model_dim, vec_model_dim)
    """
    if germLengths is None:
        germLengths = _np.array([len(germ) for germ in germsList])

    twirledDeriv = bulk_twirled_deriv(model, germsList, eps, check, comm) / germLengths[:, None, None]

    #OLD: slow, I think because conjugate *copies* a large tensor, causing a memory bottleneck
    #twirledDerivDaggerDeriv = _np.einsum('ijk,ijl->ikl',
    #                                     _np.conjugate(twirledDeriv),
    #                                     twirledDeriv)

    #NEW: faster, one-germ-at-a-time computation requires less memory.
    nGerms, _, vec_model_dim = twirledDeriv.shape
    twirledDerivDaggerDeriv = _np.empty((nGerms, vec_model_dim, vec_model_dim),
                                        dtype=_np.complex)
    for i in range(nGerms):
        twirledDerivDaggerDeriv[i, :, :] = _np.dot(
            twirledDeriv[i, :, :].conjugate().T, twirledDeriv[i, :, :])

    return twirledDerivDaggerDeriv 

def dia(mol, dm0, gauge_orig=None, shielding_nuc=None, mb='RMB'):
    '''Note the side effects of set_common_origin'''

    if shielding_nuc is None:
        shielding_nuc = range(mol.natm)
    if gauge_orig is not None:
        mol.set_common_origin(gauge_orig)

    n4c = dm0.shape[0]
    n2c = n4c // 2
    msc_dia = []
    for n, atm_id in enumerate(shielding_nuc):
        mol.set_rinv_origin(mol.atom_coord(atm_id))
        if mb.upper() == 'RMB':
            if gauge_orig is None:
                t11 = mol.intor('int1e_giao_sa10sa01_spinor', 9)
                t11 += mol.intor('int1e_spgsa01_spinor', 9)
            else:
                t11 = mol.intor('int1e_cg_sa10sa01_spinor', 9)
        elif gauge_orig is None:
            t11 = mol.intor('int1e_spgsa01_spinor', 9)
        else:
            t11 = numpy.zeros(9)
        h11 = numpy.zeros((9, n4c, n4c), complex)
        for i in range(9):
            h11[i,n2c:,:n2c] = t11[i] * .5
            h11[i,:n2c,n2c:] = t11[i].conj().T * .5
        a11 = numpy.real(numpy.einsum('xij,ji->x', h11, dm0))
        #    XX, XY, XZ, YX, YY, YZ, ZX, ZY, ZZ = 1..9
        #  => [[XX, XY, XZ], [YX, YY, YZ], [ZX, ZY, ZZ]]
        msc_dia.append(a11)
    return numpy.array(msc_dia).reshape(-1, 3, 3) 

def _call_vhf1(mol, dm):
    c1 = .5 / lib.param.LIGHT_SPEED
    n2c = dm.shape[0] // 2
    dmll = dm[:n2c,:n2c].copy()
    dmls = dm[:n2c,n2c:].copy()
    dmsl = dm[n2c:,:n2c].copy()
    dmss = dm[n2c:,n2c:].copy()
    vj = numpy.zeros((3,n2c*2,n2c*2), dtype=numpy.complex)
    vk = numpy.zeros((3,n2c*2,n2c*2), dtype=numpy.complex)
    vj[:,:n2c,:n2c], vk[:,:n2c,:n2c] = \
            _vhf.rdirect_mapdm('int2e_ip1_spinor', 's2kl',
                               ('lk->s1ij', 'jk->s1il'), dmll, 3,
                               mol._atm, mol._bas, mol._env)
    vj[:,n2c:,n2c:], vk[:,n2c:,n2c:] = \
            _vhf.rdirect_mapdm('int2e_ipspsp1spsp2_spinor', 's2kl',
                               ('lk->s1ij', 'jk->s1il'), dmss, 3,
                               mol._atm, mol._bas, mol._env) * c1**4
    vx = _vhf.rdirect_bindm('int2e_ipspsp1_spinor', 's2kl',
                            ('lk->s1ij', 'jk->s1il'), (dmll, dmsl), 3,
                            mol._atm, mol._bas, mol._env) * c1**2
    vj[:,n2c:,n2c:] += vx[0]
    vk[:,n2c:,:n2c] += vx[1]
    vx = _vhf.rdirect_bindm('int2e_ip1spsp2_spinor', 's2kl',
                            ('lk->s1ij', 'jk->s1il'), (dmss, dmls), 3,
                            mol._atm, mol._bas, mol._env) * c1**2
    vj[:,:n2c,:n2c] += vx[0]
    vk[:,:n2c,n2c:] += vx[1]
    return vj, vk 

def get_ovlp(mol):
    n2c = mol.nao_2c()
    n4c = n2c * 2
    c = lib.param.LIGHT_SPEED

    s  = mol.intor('int1e_ipovlp_spinor', comp=3)
    t  = mol.intor('int1e_ipkin_spinor', comp=3)
    s1e = numpy.zeros((3,n4c,n4c), numpy.complex)
    s1e[:,:n2c,:n2c] = s
    s1e[:,n2c:,n2c:] = t * (.5/c**2)
    return -s1e 

def get_hcore(mol):
    n2c = mol.nao_2c()
    n4c = n2c * 2
    c = lib.param.LIGHT_SPEED

    t  = mol.intor('int1e_ipkin_spinor', comp=3)
    vn = mol.intor('int1e_ipnuc_spinor', comp=3)
    wn = mol.intor('int1e_ipspnucsp_spinor', comp=3)
    h1e = numpy.zeros((3,n4c,n4c), numpy.complex)
    h1e[:,:n2c,:n2c] = vn
    h1e[:,n2c:,:n2c] = t
    h1e[:,:n2c,n2c:] = t
    h1e[:,n2c:,n2c:] = wn * (.25/c**2) - t
    return -h1e 

def test_r_incore(self):
        j3c = df.r_incore.aux_e2(mol, auxmol, intor='int3c2e_spinor', aosym='s1')
        nao = mol.nao_2c()
        naoaux = auxmol.nao_nr()
        j3c = j3c.reshape(nao,nao,naoaux)

        eri0 = numpy.empty((nao,nao,naoaux), dtype=numpy.complex)
        pi = 0
        for i in range(mol.nbas):
            pj = 0
            for j in range(mol.nbas):
                pk = 0
                for k in range(mol.nbas, mol.nbas+auxmol.nbas):
                    shls = (i, j, k)
                    buf = gto.moleintor.getints_by_shell('int3c2e_spinor',
                                                         shls, atm, bas, env)
                    di, dj, dk = buf.shape
                    eri0[pi:pi+di,pj:pj+dj,pk:pk+dk] = buf
                    pk += dk
                pj += dj
            pi += di
        self.assertTrue(numpy.allclose(eri0, j3c))
        eri1 = df.r_incore.aux_e2(mol, auxmol, intor='int3c2e_spinor',
                                  aosym='s2ij')
        for i in range(naoaux):
            j3c[:,:,i] = lib.unpack_tril(eri1[:,i])
        self.assertTrue(numpy.allclose(eri0, j3c)) 

def test_sp_spinor(self):
        ao = eval_gto(mol, 'GTOval_sp_spinor', coords)
        self.assertAlmostEqual(finger(ao), 26.212937567473656-68.970076521029782j, 9)

        numpy.random.seed(1)
        rs = numpy.random.random((213,3))
        rs = (rs-.5)**2 * 30
        ao1 = eval_gto(mol1, 'GTOval_spinor', rs, shls_slice=(0,mol1.nbas//2))
        ao2 = eval_gto(mol1, 'GTOval_sp_spinor', rs, shls_slice=(mol1.nbas//2,mol1.nbas))
        self.assertAlmostEqual(abs(ao1-ao2*1j).sum(), 0, 9)
#
#        t = gto.cart2j_kappa(0, 2)
#        aonr = eval_gto(mol1, 'GTOval_cart', rs, shls_slice=(1,2))
#        aa = numpy.zeros((2,12))
#        aa[:1,:6] = aonr[:,:6]
#        aa[1:,6:] = aonr[:,:6]
#        print aa.dot(t[:,:4])
#
#        t = gto.cart2j_kappa(0, 1)/0.488602511902919921
#        aonr = eval_gto(mol1, 'GTOval_ip_cart', rs, comp=3, shls_slice=(mol1.nbas//2+1,mol1.nbas//2+2))
#        print 'x', aonr[0,0,:3]
#        print 'y', aonr[1,0,:3]
#        print 'z', aonr[2,0,:3]
#        aa = numpy.zeros((3,2,6),dtype=numpy.complex)
#        aa[0,:1,3:] = aonr[0,0,:3]
#        aa[0,1:,:3] = aonr[0,0,:3]
#        aa[1,:1,3:] =-aonr[1,0,:3]*1j
#        aa[1,1:,:3] = aonr[1,0,:3]*1j
#        aa[2,:1,:3] = aonr[2,0,:3]
#        aa[2,1:,3:] =-aonr[2,0,:3]
#        aa = (aa[0]*-1j + aa[1]*-1j + aa[2]*-1j)
#        print 'p',aa.dot(t[:,2:])*1j 

def takebak_2d(out, a, idx, idy):
    '''Reverse operation of take_2d.  Equivalent to out[idx[:,None],idy] += a
    for a 2D array.

    Examples:

    >>> out = numpy.zeros((3,3))
    >>> takebak_2d(out, numpy.ones((2,2)), [0,2], [0,2])
    [[ 1.  0.  1.]
     [ 0.  0.  0.]
     [ 1.  0.  1.]]
    '''
    assert(out.flags.c_contiguous)
    a = numpy.asarray(a, order='C')
    idx = numpy.asarray(idx, dtype=numpy.int32)
    idy = numpy.asarray(idy, dtype=numpy.int32)
    if a.dtype == numpy.double:
        fn = _np_helper.NPdtakebak_2d
    elif a.dtype == numpy.complex:
        fn = _np_helper.NPztakebak_2d
    else:
        out[idx[:,None], idy] += a
        return out
    fn(out.ctypes.data_as(ctypes.c_void_p),
       a.ctypes.data_as(ctypes.c_void_p),
       idx.ctypes.data_as(ctypes.c_void_p),
       idy.ctypes.data_as(ctypes.c_void_p),
       ctypes.c_int(out.shape[1]), ctypes.c_int(a.shape[1]),
       ctypes.c_int(idx.size), ctypes.c_int(idy.size))
    return out