Python numpy.complex() Examples

The following are 30 code examples of numpy.complex(). 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: eom_kccsd_ghf.py    From pyscf with Apache License 2.0 6 votes vote down vote up
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 
Example #3
Source File: dhf.py    From pyscf with Apache License 2.0 6 votes vote down vote up
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 
Example #4
Source File: uintermediates_slow.py    From pyscf with Apache License 2.0 6 votes vote down vote up
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 
Example #5
Source File: uintermediates_slow.py    From pyscf with Apache License 2.0 6 votes vote down vote up
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 
Example #6
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 #7
Source File: linalg_helper.py    From pyscf with Apache License 2.0 6 votes vote down vote up
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 
Example #8
Source File: test_utils.py    From me-ica with GNU Lesser General Public License v2.1 6 votes vote down vote up
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) 
Example #9
Source File: eom_kccsd_ghf.py    From pyscf with Apache License 2.0 6 votes vote down vote up
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 
Example #10
Source File: dhf.py    From pyscf with Apache License 2.0 6 votes vote down vote up
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 
Example #11
Source File: test_moleintor.py    From pyscf with Apache License 2.0 6 votes vote down vote up
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 
Example #12
Source File: __init__.py    From tenpy with GNU General Public License v3.0 6 votes vote down vote up
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 
Example #13
Source File: dhf.py    From pyscf with Apache License 2.0 5 votes vote down vote up
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 
Example #14
Source File: test_arraywriters.py    From me-ica with GNU Lesser General Public License v2.1 5 votes vote down vote up
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) 
Example #15
Source File: test_ufunclike.py    From recruit with Apache License 2.0 5 votes vote down vote up
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) 
Example #16
Source File: test_indexing.py    From recruit with Apache License 2.0 5 votes vote down vote up
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 
Example #17
Source File: test_ufunclike.py    From recruit with Apache License 2.0 5 votes vote down vote up
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) 
Example #18
Source File: test_common.py    From recruit with Apache License 2.0 5 votes vote down vote up
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])) 
Example #19
Source File: test_arraywriters.py    From me-ica with GNU Lesser General Public License v2.1 5 votes vote down vote up
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) 
Example #20
Source File: gaussian.py    From vampyre with MIT License 5 votes vote down vote up
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) 
Example #21
Source File: dhf.py    From pyscf with Apache License 2.0 5 votes vote down vote up
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 
Example #22
Source File: dhf.py    From pyscf with Apache License 2.0 5 votes vote down vote up
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 
Example #23
Source File: germselection.py    From pyGSTi with Apache License 2.0 5 votes vote down vote up
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 
Example #24
Source File: dhf.py    From pyscf with Apache License 2.0 5 votes vote down vote up
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) 
Example #25
Source File: dhf.py    From pyscf with Apache License 2.0 5 votes vote down vote up
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 
Example #26
Source File: dhf.py    From pyscf with Apache License 2.0 5 votes vote down vote up
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 
Example #27
Source File: dhf.py    From pyscf with Apache License 2.0 5 votes vote down vote up
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 
Example #28
Source File: test_incore.py    From pyscf with Apache License 2.0 5 votes vote down vote up
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)) 
Example #29
Source File: test_gridao.py    From pyscf with Apache License 2.0 5 votes vote down vote up
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 
Example #30
Source File: numpy_helper.py    From pyscf with Apache License 2.0 5 votes vote down vote up
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