How to use the pyscf.scf.RHF function in pyscf

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)

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

def test_info_functions_mol():
    from pyscf import gto, scf
    from pyqmc.tbdm import TBDMAccumulator

    mol = gto.Mole()
    mol.atom = """He 0.00 0.00 0.00 """
    mol.basis = "ccpvdz"
    mol.build()

    mf = scf.RHF(mol)
    ehf = mf.kernel()

    wf, to_opt = pyqmc.default_sj(mol, mf)
    accumulators = {
        "pgrad": pyqmc.gradient_generator(mol, wf, to_opt),
        "obdm": OBDMAccumulator(mol, orb_coeff=mf.mo_coeff),
        "tbdm_updown": TBDMAccumulator(mol, np.asarray([mf.mo_coeff] * 2), (0, 1)),
    }
    info_functions(mol, wf, accumulators)

dmrgscf.settings.MPIPREFIX ='mpirun -n 3'

b = 1.4
mol = gto.Mole()
mol.build(
    verbose = 4,
    output = 'fci_nevpt2.out',
    atom = [['H', (0.,0.,i-3.5)] for i in range(8)],
    basis = 'sto-3g',
#
# Note if symmetry is specified, the molecule has to be placed on the proper
# orientation to match the given symmetry.
#
    symmetry = 'd2h',
)
m = scf.RHF(mol)
m.kernel()

#
# FCI-based CASCI + NEVPT2.  Two roots are computed.  mc.ci holds the two CI
# wave functions.  Root ID needs to be specified for the state-specific NEVPT2
# calculation.  By default the lowest root is computed.
#
mc = mcscf.CASCI(m, 4, 4)
mc.fcisolver.nroots = 2
mc.kernel()

ci_nevpt_e1 = mrpt.NEVPT(mc, root=0).kernel()
ci_nevpt_e2 = mrpt.NEVPT(mc, root=1).kernel()

print('CI NEVPT = %.15g %.15g' % (ci_nevpt_e1, ci_nevpt_e2))

Bohr unit, one method is to loop over the internal structure and multiply the
coordinates with 0.52917721092 (pyscf.lib.parameters) to convert Bohr to Angstrom.
We might add another attribute for Mole object to indicate the unit in the future.

Note in this example, the diffuse functions in the basis might cause linear
dependency in the AO space.  It can affect the accuracy of the SCF calculation.
See also example  42-remove_linear_dep.py  to remove the linear dependency.
'''

mol = gto.Mole()
mol.verbose = 3
mol.atom = [['H', (0, 0, i*1.0)] for i in range(10)]
mol.basis = 'aug-ccpvdz'
mol.build()

mf = scf.RHF(mol)
print('condition number', numpy.linalg.cond(mf.get_ovlp()))
mf.scf()

if __name__ == '__main__':
    from pyscf import gto
    from pyscf import scf
    from pyscf.lo import nao
    mol = gto.Mole()
    mol.verbose = 1
    mol.output = 'out_orth'
    mol.atom.extend([
        ['O' , (0. , 0.     , 0.)],
        [1   , (0. , -0.757 , 0.587)],
        [1   , (0. , 0.757  , 0.587)] ])
    mol.basis = {'H': '6-31g',
                 'O': '6-31g',}
    mol.build()

    mf = scf.RHF(mol)
    mf.scf()

    c0 = nao.prenao(mol, mf.make_rdm1())
    c = orth_ao(mol, 'meta_lowdin', c0)

    s = mol.intor_symmetric('cint1e_ovlp_sph')
    p = reduce(numpy.dot, (s, mf.make_rdm1(), s))
    print(reduce(numpy.dot, (c.T, p, c)).diagonal())

from pyscf.dmrgscf import settings
if 'SLURMD_NODENAME' in os.environ:  # slurm system
    settings.MPIPREFIX = 'srun'
elif 'PBS_NODEFILE' in os.environ:   # PBS system
    settings.MPIPREFIX = 'mpirun'
else:  # MPI on single node
    settings.MPIPREFIX = 'mpirun -np 4'

b = 1.2
mol = gto.M(
    verbose = 4,
    atom = 'N 0 0 0; N 0 0 %f'%b,
    basis = 'cc-pvdz',
    symmetry = True,
)
mf = scf.RHF(mol)
mf.kernel()

mc = dmrgscf.DMRGSCF(mf, 8, 8)
emc = mc.kernel()[0]
print(emc)

mc = dmrgscf.DMRGSCF(mf, 8, 8)
mc.state_average_([0.5, 0.5])
mc.kernel()
print(mc.e_tot)

'''
Applying scalar relativistic effects for CASSCF by decorating the SCF or
CASSCF object with .x2c method.

See pyscf/examples/scf/21-x2c.py
'''

mol = gto.M(
    verbose = 0,
    atom = '''8  0  0.     0
              1  0  -0.757 0.587
              1  0  0.757  0.587''',
    basis = 'ccpvdz',
)

mf = scf.RHF(mol).x2c()
mf.kernel()

mc = mcscf.CASSCF(mf, 6, 8)
mc.kernel()
print('E = %.12f, ref = -76.128478294795' % mc.e_tot)

#
# Decorating CASSCF with .x2c method has the same effects as decorating SCF object
#
mf = scf.RHF(mol)
mf.kernel()

mc = mcscf.CASSCF(mf, 6, 8).x2c()
mc.kernel()
print('E = %.12f, ref = -76.128478294795' % mc.e_tot)

It's important to compare multi-reference calculations based on the comparable
reference states.  Here we compute the SVD eig and the determinant value of
the CAS space overlap to measure how close two CASSCF results are.  If two
CASSCF are close, the SVD eig should be close 1
'''

mol1 = gto.M(atom='O 0 0 0; O 0 0 1', basis='ccpvtz', spin=2, symmetry=1)
mf = scf.RHF(mol1)
mf.kernel()
mc = mcscf.CASSCF(mf, 7, 4)
mc.kernel()
mo1core = mc.mo_coeff[:,:mc.ncore]
mo1cas = mc.mo_coeff[:,mc.ncore:mc.ncore+mc.ncas]

mol2 = gto.M(atom='O 0 0 0; O 0 0 1', basis='ccpvdz', spin=2, symmetry=1)
mf = scf.RHF(mol2)
mf.kernel()
mc = mcscf.CASSCF(mf, 7, (2,2))
mc.kernel()
mo2core = mc.mo_coeff[:,:mc.ncore]
mo2cas = mc.mo_coeff[:,mc.ncore:mc.ncore+mc.ncas]

s = gto.intor_cross('cint1e_ovlp_sph', mol1, mol2)
score = reduce(numpy.dot, (mo1core.T, s, mo2core))
scas = reduce(numpy.dot, (mo1cas.T, s, mo2cas))
numpy.set_printoptions(4)
print(' SVD eig = %s' % numpy.linalg.svd(score)[1])
print('det()    = %s' % numpy.linalg.det(score))
print(' SVD eig = %s' % numpy.linalg.svd(scas)[1])
print('det()    = %s' % numpy.linalg.det(scas))

#

'''
TDDFT analytical nuclear gradients.
'''

from pyscf import gto, scf, dft, tddft

mol = gto.M(
    atom = [
        ['O' , 0. , 0.     , 0],
        ['H' , 0. , -0.757 , 0.587],
        ['H' , 0. ,  0.757 , 0.587]],
    basis = 'ccpvdz')

mf = scf.RHF(mol).run()
postmf = tddft.TDHF(mf).run()
g = postmf.nuc_grad_method()
g.kernel(state=1)

mf = dft.RKS(mol).x2c().set(xc='pbe0').run()
# Switch to xcfun because 3rd order GGA functional derivative is not
# available in libxc
mf._numint.libxc = dft.xcfun
postmf = tddft.TDDFT(mf).run()
# PySCF-1.6.1 and newer supports the .Gradients method to create a grad
# object after grad module was imported. It is equivalent to call the
# .nuc_grad_method method.
from pyscf import grad
g = postmf.Gradients()
g.kernel(state=1)