How to use the pyscf.gto.M function in pyscf

#!/usr/bin/env python

'''
Mulliken population analysis with NAO
'''

import numpy
from pyscf import gto, scf, lo
from scipy._lib.six import reduce

x = .63
mol = gto.M(atom=[['C', (0, 0, 0)],
                  ['H', (x ,  x,  x)],
                  ['H', (-x, -x,  x)],
                  ['H', (-x,  x, -x)],
                  ['H', ( x, -x, -x)]],
            basis='ccpvtz')
mf = scf.RHF(mol).run()

# C matrix stores the AO to localized orbital coefficients
C = lo.orth_ao(mf, 'nao')

# C is orthogonal wrt to the AO overlap matrix.  C^T S C  is an identity matrix.
print(abs(reduce(numpy.dot, (C.T, mf.get_ovlp(), C)) -
          numpy.eye(mol.nao_nr())).max())  # should be close to 0

# The following linear equation can also be solved using the matrix
# multiplication reduce(numpy.dot (C.T, mf.get_ovlp(), mf.mo_coeff))

'''
CCSD(T) does not have the interface to the geometry optimizer berny_solver.
You need to define a function to compute CCSD(T) total energy and gradients
then use "as_pyscf_method" to pass them to berny_solver.

See also  examples/geomopt/02-as_pyscf_method.py
'''

from pyscf import gto
from pyscf import scf
from pyscf import cc
from pyscf.cc import ccsd_t_lambda_slow as ccsd_t_lambda
from pyscf.grad import ccsd_t as ccsd_t_grad
from pyscf.geomopt import berny_solver

mol = gto.M(
    verbose = 3,
    atom = [
        ['O' , (0. , 0.     , 0.    )],
        ['H' , (0. ,-0.757  ,-0.587)],
        ['H' , (0. , 0.757  ,-0.587)]],
    basis = 'ccpvdz'
)
mf = scf.RHF(mol)
cc_scan = cc.CCSD(mf).as_scanner()

def f(mol):
    # Compute CCSD(T) energy
    mf = scf.RHF(mol).run()
    mycc = cc.CCSD(mf).run()
    et_correction = mycc.ccsd_t()
    e_tot = mycc.e_tot + et_correction

mol = gto.M(atom='O 0 0 0; O 0 0 1.2', spin=2, basis='sto3g',
            verbose=0)
mf = scf.RHF(mol).run()
mci = fci.FCI(mol, mf.mo_coeff)
mci = fci.addons.fix_spin_(mci, ss=0)
e, civec = mci.kernel(nelec=nelec)
print('Singlet E = %.12f  2S+1 = %.7f' %
      (e, mci.spin_square(civec, mf.mo_coeff.shape[1], nelec)[1]))


#
# Be careful with the trick of energy penalty.  Numerical problems may be
# observed when function fix_spin_ was applied.
#

mol = gto.M(atom='O 0 0 0', basis='6-31G')
m = scf.RHF(mol).run()
norb = m.mo_coeff.shape[1]
nelec = mol.nelec
fs = fci.addons.fix_spin_(fci.FCI(mol, m.mo_coeff), .5)
fs.nroots = 15
e, fcivec = fs.kernel(verbose=5)
# The first 5 states should be degenerated. The degeneracy may be broken.
for i, c in enumerate(fcivec):
    print('state = %d, E = %.9f, S^2=%.4f' %
          (i, e[i], fci.spin_op.spin_square(c, norb, nelec)[0]))

#
# Author: Qiming Sun 
#

import numpy
from pyscf import gto, symm, scf

'''
Symmetrize orbital space
'''

mol = gto.M(atom = '''C  0  0  0
                      H  1  1  1
                      H -1 -1  1
                      H  1 -1 -1
                      H -1  1 -1''',
            basis = 'sto3g',
            verbose = 0)
mf = scf.RHF(mol).run()
mo = mf.mo_coeff

#
# call mol.build to construct symmetry adapted basis.
#
# NOTE the molecule orientation.  Simply assigning mol.build(symmetry=True)
# may change the orientation of the molecule.  If the orientation is different
# to the orientation of the orbital space, the orbitals cannot be symmetrized.
# To keep the molecule orientation fixed, we need specify the molecular

mcs = mp.as_scanner()
    mol.set_geom_([
            ["O" , (0. , 0.     , 0.001)],
            [1   , (0. ,-0.757  , 0.587)],
            [1   , (0. , 0.757  , 0.587)]])
    e1 = mcs(mol)
    mol.set_geom_([
            ["O" , (0. , 0.     ,-0.001)],
            [1   , (0. ,-0.757  , 0.587)],
            [1   , (0. , 0.757  , 0.587)]])
    e2 = mcs(mol)
    print(g1[0,2], (e1-e2)/0.002*lib.param.BOHR)

    print('-----------------------------------')
    mol = gto.M(
        atom = [
            ["O" , (0. , 0.     , 0.)],
            [1   , (0. ,-0.757  , 0.587)],
            [1   , (0. , 0.757  , 0.587)]],
        basis = '631g'
    )
    mf = scf.RHF(mol).run()
    mp = mp2.MP2(mf)
    mp.frozen = [0,1,10,11,12]
    mp.max_memory = 1
    mp.kernel()
    g1 = Gradients(mp).kernel()
# O    -0.0000000000    -0.0000000000     0.0037319667
# H    -0.0000000000    -0.0897959298    -0.0018659834
# H     0.0000000000     0.0897959298    -0.0018659834
    print(lib.finger(g1) - 0.12458103614793946)

def to_scf(filename, molpro_orbsym=MOLPRO_ORBSYM, mf=None, **kwargs):
    '''Use the Hamiltonians defined by FCIDUMP to build an SCF object'''
    ctx = read(filename, molpro_orbsym)
    mol = gto.M()
    mol.nelectron = ctx['NELEC']
    mol.spin = ctx['MS2']
    norb = mol.nao = ctx['NORB']
    if 'ECORE' in ctx:
        mol.energy_nuc = lambda *args: ctx['ECORE']
    mol.incore_anyway = True

    if 'ORBSYM' in ctx:
        mol.symmetry = True
        mol.groupname = 'N/A'
        orbsym = numpy.asarray(ctx['ORBSYM'])
        mol.irrep_id = list(set(orbsym))
        mol.irrep_name = [('IR%d' % ir) for ir in mol.irrep_id]
        so = numpy.eye(norb)
        mol.symm_orb = []
        for ir in mol.irrep_id:

from pyscf import gto, scf, mcscf

'''
Tune CASSCF parameters to reduce the FCI solver cost.

The default CASSCF settings are optimized for large system so that the cost of
IO and integral transformation are minimal.  But they are not optimal for
small systems which call expensive FCI solver.  Taking high order expansion for
matrix elements (.ci_update_dep) and more micro iterations (.max_cycle_micro)
can increase the cost of integration but reduce the total time needed FCI solver.

These settings are useful when DMRG-SCF is executed for small systems.
'''

mol = gto.M(atom='''
C   0       0          0
H  .990138, -0.436705  0
H -.990138, -0.436705  0''',
            basis = 'ccpvdz',
            symmetry = 1,
            spin = 2)

mc = mol.apply(scf.RHF).run().apply(mcscf.CASSCF, 14, 6).run(verbose=4, with_dep4=True, max_cycle_micro=10)

mf.kernel()

#
# In the ROHF method, one alpha electron needs to be put in the degenerated
# pi_x*, pi_y* orbitals.
#
mol = gto.M(atom='O 0 0 0; O 0 0 1', charge=2, spin=2)
mf = scf.rohf.ROHF(mol)
mf.verbose = 4
mf = scf.addons.frac_occ(mf)
mf.kernel()

#
# One alpha electron in the degenerated pi_x*, pi_y* orbitals for UHF method.
#
mol = gto.M(atom='O 0 0 0; O 0 0 1', charge=1, spin=1)
mf = scf.UHF(mol)
mf.verbose = 4
mf = scf.addons.frac_occ(mf)
mf.kernel()

detmp = numpy.einsum('ng,nxg->nx', eta_nj, den1)
            de[:,ja] += detmp
            de[:,ia] -= detmp

    psi1_dm1 -= numpy.einsum('nil,aziljm->nazjm', LS0, L1)
    LS1 = numpy.linalg.solve(L.reshape(natm*nlm,-1).T, psi1_dm1.reshape(-1,natm*nlm).T)
    LS1 = LS1.T.reshape(n_dm,natm,3,natm,nlm)
    de += numpy.einsum('nazjx,njx->naz', LS1, phi0_dm2)

    if is_single_dm:
        de = de[0]
    return de


if __name__ == '__main__':
    mol0 = gto.M(atom='H  0.  0.  1.804; F  0.  0.  0.', verbose=0, unit='B')
    mol1 = gto.M(atom='H  0.  0.  1.803; F  0.  0.  0.', verbose=0, unit='B')
    mol2 = gto.M(atom='H  0.  0.  1.805; F  0.  0.  0.', verbose=0, unit='B')

    # TDA with equilibrium_solvation
    mf = mol0.RHF().ddCOSMO().run()
    td = mf.TDA().ddCOSMO().run(equilibrium_solvation=True)
    g1 = td.nuc_grad_method().kernel() # 0  0  -0.5116214042

    mf = mol1.RHF().ddCOSMO().run()
    td1 = mf.TDA().ddCOSMO().run(equilibrium_solvation=True)
    mf = mol2.RHF().ddCOSMO().run()
    td2 = mf.TDA().ddCOSMO().run(equilibrium_solvation=True)
    print((td2.e_tot[0]-td1.e_tot[0])/0.002, g1[0,2])
    print((td2.e_tot[0]-td1.e_tot[0])/0.002 - g1[0,2])

    # TDA without equilibrium_solvation

if kernel is None:
        return vxc_pack(self.sv, dm, xc_code, deriv=2, kernel=kernel, ao_log=self.prod_log, **kvargs)
    else:
        vxc_pack(self.sv, dm, xc_code, deriv=2, kernel=kernel, ao_log=self.prod_log, **kvargs)


#
#
#
if __name__=='__main__':
  from pyscf.nao import prod_basis_c, system_vars_c, comp_overlap_coo
  from pyscf import gto
  import numpy as np
  
  mol = gto.M(atom='O 0 0 0; H 0 0 0.5; H 0 0.5 0', basis='ccpvdz') # coordinates in Angstrom!
  sv = system_vars_c(gto=mol)
  print(sv.atom2s)
  s_ref = comp_overlap_coo(sv).todense()
  pb = prod_basis_c(sv)
  mom0,mom1=pb.comp_moments()
  pab2v = pb.get_vertex_array()
  s_chk = einsum('pab,p->ab', pab2v,mom0)
  print(abs(s_chk-s_ref).sum()/s_chk.size, abs(s_chk-s_ref).max())