How to use pyscf - 10 common examples

# bar{C}_{pi}^1 ~= C_{pi}^1 - 
    if sscobj.mb.upper() == 'RMB':
        orbo = mo_coeff[:,mo_occ> 0]
        orbv = mo_coeff[:,mo_occ==0]
        n4c = mo_coeff.shape[0]
        n2c = n4c // 2
        c = lib.param.LIGHT_SPEED
        orbvS_T = orbv[n2c:].conj().T
        for ia in atmlst:
            mol.set_rinv_origin(mol.atom_coord(ia))
            a01int = mol.intor('int1e_sa01sp_spinor', 3)
            for k in range(3):
                s1 = orbvS_T.dot(a01int[k].conj().T).dot(orbo[n2c:])
                mo1[ia*3+k] -= s1 * (.25/c**2)

    logger.timer(sscobj, 'solving mo1 eqn', *cput1)
    return mo1, mo_e1

def get_pnucp(mydf, kpts=None):
    cell = mydf.cell
    if kpts is None:
        kpts_lst = numpy.zeros((1,3))
    else:
        kpts_lst = numpy.reshape(kpts, (-1,3))

    log = logger.Logger(mydf.stdout, mydf.verbose)
    t1 = t0 = (time.clock(), time.time())

    nkpts = len(kpts_lst)
    nao = cell.nao_nr()
    nao_pair = nao * (nao+1) // 2

    Gv, Gvbase, kws = cell.get_Gv_weights(mydf.mesh)
    charge = -cell.atom_charges()
    kpt_allow = numpy.zeros(3)
    coulG = tools.get_coulG(cell, kpt_allow, mesh=mydf.mesh, Gv=Gv)
    coulG *= kws
    if mydf.eta == 0:
        wj = numpy.zeros((nkpts,nao_pair), dtype=numpy.complex128)
        SI = cell.get_SI(Gv)
        vG = numpy.einsum('i,ix->x', charge, SI) * coulG
        if cell.dimension == 1 or cell.dimension == 2:

nd = (l + 1) * (l + 2) // 2
            else:
                nd = l * 2 + 1
            p0, p1 = p1, p1 + nd * nc
            if l <= 4:
                idx.append(range(p0, p1))

        idx = numpy.hstack(idx)
        return pmol, mo_coeff[idx]



if __name__ == '__main__':
    from pyscf import scf
    import tempfile
    mol = gto.Mole()
    mol.verbose = 5
    mol.output = None#'out_gho'
    mol.atom = [['C', (0.,0.,0.)],
                ['H', ( 1, 1, 1)],
                ['H', (-1,-1, 1)],
                ['H', ( 1,-1,-1)],
                ['H', (-1, 1,-1)], ]
    mol.basis = {
        'C': 'sto-3g',
        'H': 'sto-3g'}
    mol.build(dump_input=False)
    m = scf.RHF(mol)
    m.scf()
    header(mol, mol.stdout)
    print(order_ao_index(mol))
    orbital_coeff(mol, mol.stdout, m.mo_coeff)

# Author: Qiming Sun 
#

r'''
Applying creation or annihilation operators on FCI wavefunction
        a |0>

Compute density matrices by 
        gamma_{ij} = <0| i^+ j |0>
        Gamma_{ij,kl} = <0| i^+ j^+ l k |0>
'''

import numpy
from pyscf import gto, scf, fci

mol = gto.M(atom='H 0 0 0; Li 0 0 1.1', basis='sto3g')
m = scf.RHF(mol).run()
fs = fci.FCI(mol, m.mo_coeff)
e, c = fs.kernel()

norb = m.mo_energy.size
neleca = nelecb = mol.nelectron // 2

#
# Spin-free 1-particle density matrix
# 
#
dm1 = numpy.zeros((norb,norb))
for i in range(norb):
    for j in range(norb):
        tmp = fci.addons.des_a(c  , norb, (neleca  ,nelecb), j)
        tmp = fci.addons.cre_a(tmp, norb, (neleca-1,nelecb), i)

print(abs(l2new-l2new.transpose(1,0,3,2)).sum())
    print(numpy.dot(l2new.flatten(), numpy.arange(35**2)) - 48427109.5409886)
    print(numpy.dot(l2new.flatten(), numpy.sin(numpy.arange(35**2)))-137.758016736487)
    print(numpy.dot(numpy.sin(l2new.flatten()), numpy.arange(35**2))-507.656936701192)


    mol = gto.Mole()
    mol.verbose = 0
    mol.atom = [
        [8 , (0. , 0.     , 0.)],
        [1 , (0. , -0.757 , 0.587)],
        [1 , (0. , 0.757  , 0.587)]]

    mol.basis = 'cc-pvdz'
    mol.build()
    rhf = scf.RHF(mol)
    rhf.conv_tol = 1e-16
    rhf.scf()

    mcc = ccsd.CCSD(rhf)
    mcc.conv_tol = 1e-12
    ecc, t1, t2 = mcc.kernel()

    nmo = rhf.mo_energy.size
    fock0 = numpy.diag(rhf.mo_energy)
    nocc = mol.nelectron // 2
    nvir = nmo - nocc

    eris = mcc.ao2mo()
    conv, l1, l2 = kernel(mcc, eris, t1, t2, tol=1e-8)
    print(numpy.linalg.norm(l1)-0.0132626841292)
    print(numpy.linalg.norm(l2)-0.212575609057)

def amplitudes_to_vector_singlet(r1, r2, kconserv):
    '''Transform 3- and 7-dimensional arrays, r1 and r2, to a 1-dimensional
    array with unique indices.

    For example:
        r1: t_{i k_i}^{a k_a}
        r2: t_{i k_i, j k_j}^{a k_a, b k_b}
        return: a vector with all r1 elements, and r2 elements whose indices
    satisfy (i k_i a k_a) >= (j k_j b k_b)
    '''
    cput0 = (time.clock(), time.time())
    log = logger.Logger(sys.stdout, logger.DEBUG)
    # r1 indices: k_i, i, a
    nkpts, nocc, nvir = np.asarray(r1.shape)[[0, 1, 2]]
    nov = nocc * nvir

    # r2 indices (old): k_i, k_J, k_a, i, J, a, B
    # r2 indices (new): (k_i, k_a), (k_J), (i, a), (J, B)
    r2 = r2.transpose(0,2,1,3,5,4,6).reshape(nkpts**2, nkpts, nov, nov)

    idx, idy = np.tril_indices(nov)
    nov2_tril = nov * (nov + 1) // 2
    nov2 = nov * nov

    vector = np.empty(r2.size, dtype=r2.dtype)
    offset = 0
    for ki, ka, kj in kpts_helper.loop_kkk(nkpts):
        kb = kconserv[ki, ka, kj]

if isinstance(nelec, (int, numpy.integer)):
        if shciobj.spin is None:
            nelecb = nelec // 2
        else:
            nelecb = (nelec - shciobj.spin) // 2
        neleca = nelec - nelecb
    else :
        neleca, nelecb = nelec

    if shciobj.groupname is not None and shciobj.orbsym is not []:
# First removing the symmetry forbidden integrals. This has been done using
# the pyscf internal irrep-IDs (stored in shciobj.orbsym)
        orbsym = numpy.asarray(shciobj.orbsym) % 10
        pair_irrep = (orbsym.reshape(-1,1) ^ orbsym)[numpy.tril_indices(ncas)]
        sym_forbid = pair_irrep.reshape(-1,1) != pair_irrep.ravel()
        eri_cas = ao2mo.restore(4, eri_cas, ncas)
        eri_cas[sym_forbid] = 0
        eri_cas = ao2mo.restore(8, eri_cas, ncas)
        # Convert the pyscf internal irrep-ID to molpro irrep-ID
        orbsym = numpy.asarray(symmetry.convert_orbsym(shciobj.groupname, orbsym))
    else:
        orbsym = []
        eri_cas = ao2mo.restore(8, eri_cas, ncas)

    if not os.path.exists(shciobj.runtimedir):
        os.makedirs(shciobj.runtimedir)

    # The name of the FCIDUMP file, default is "FCIDUMP".
    integralFile = os.path.join(shciobj.runtimedir, shciobj.integralfile)
    tools.fcidump.from_integrals(integralFile, h1eff, eri_cas, ncas,
                                 neleca+nelecb, ecore, ms=abs(neleca-nelecb),
                                 orbsym=orbsym)

if unpack:
                buf = numpy.empty((min(blksize, naux), nao * (nao + 1) // 2))
            for b0, b1 in lib.prange(0, naux, blksize):
                LpqR, LpqI = load(j3c, b0, b1, LpqR, LpqI)
                yield LpqR, LpqI, 1

        if cell.dimension == 2 and cell.low_dim_ft_type != 'inf_vacuum':
            # Truncated Coulomb operator is not postive definite. Load the
            # CDERI tensor of negative part.
            LpqR = LpqI = None
            with _load3c(self._cderi, 'j3c-', kpti_kptj, 'j3c-kptij',
                         ignore_key_error=True) as j3c:
                naux = j3c.shape[0]
                if unpack:
                    buf = numpy.empty((min(blksize, naux), nao * (nao + 1) // 2))
                for b0, b1 in lib.prange(0, naux, blksize):
                    LpqR, LpqI = load(j3c, b0, b1, LpqR, LpqI)
                    yield LpqR, LpqI, -1

def test():

    mol = gto.M(
        atom="Li 0. 0. 0.; Li 0. 0. 1.5", basis="sto-3g", unit="bohr", verbose=0
    )
    mf = scf.RHF(mol).run()

    # Lowdin orthogonalized AO basis.
    lowdin = lo.orth_ao(mol, "lowdin")

    # MOs in the Lowdin basis.
    mo = solve(lowdin, mf.mo_coeff)

    # make AO to localized orbital coefficients.
    mfobdm = mf.make_rdm1(mo, mf.mo_occ)

    ### Test OBDM calculation.
    nconf = 500
    nsteps = 400
    warmup = 15
    wf = PySCFSlater(mol, mf)
    configs = initial_guess(mol, nconf)

C    0.000000000000    -1.398696930758     0.000000000000
     H    2.157597486829     1.245660462400     0.000000000000
     C    1.211265339156     0.699329968382     0.000000000000
     H    2.157597486829    -1.245660462400     0.000000000000
     C    1.211265339156    -0.699329968382     0.000000000000
     H   -2.157597486829     1.245660462400     0.000000000000
     C   -1.211265339156     0.699329968382     0.000000000000
     H   -2.157597486829    -1.245660462400     0.000000000000
     C   -1.211265339156    -0.699329968382     0.000000000000
  '''
mol.basis = '6-31g'
mol.symmetry = 0
mol.charge = 0
mol.spin = 0
mol.build()
mf = scf.RHF( mol )
mf.verbose = 0
mf.scf()

filename_mo       = 'benzene-631g-mo.molden'
filename_edmiston = 'benzene-631g-edmiston.molden'

with open( filename_mo, 'w' ) as thefile:
    molden.header( mol, thefile )
    molden.orbital_coeff( mol, thefile, mf.mo_coeff )
print("Molecular orbitals saved in", filename_mo)

# Localize the pi-type orbitals. Counting starts from 0! 12 orbitals as 6-31G is DZ.
tolocalize = np.array([17, 20, 21, 22, 23, 30, 36, 41, 42, 47, 48, 49]) - 1
loc  = localizer.localizer( mol, mf.mo_coeff[:,tolocalize], 'edmiston' )
loc.verbose = param.VERBOSE_DEBUG
new_coeff = loc.optimize()