4
$\begingroup$

A simple question about PySCF, but nonetheless I could not find out how to do it: I have a system of atoms for which I already calculated a STO-3G basis in the usual form of MOs: $$\chi = \sum_i c_i \phi_i\tag{1}$$ with $\phi_i$ being Gaussian primitives. How can I then parse this basis into a form so that I can use it for integral calculations in PySCF?

Thank you for your time!

$\endgroup$
1
  • 1
    $\begingroup$ +1 Welcome to our new community and thank you for contributing your question here! We hope to see much more of you in the future !!! $\endgroup$ Jan 17, 2023 at 0:35

1 Answer 1

6
$\begingroup$

Instead of:

    basis = STO-3G

you can get PySCF to read your own "custom" basis set file, for example:

    basis = {
        'H' : parse_gaussian.load('H_av2z.gbs', 'H')
    },

The .gbs file in this case, is written in Gaussian's basis format (make sure you have from pyscf.gto.basis import parse_gaussian at the top of your script).

An example of a .gbs file is this one:

spherical

****
H     0
S    4   1.00
      1.301000D+01           1.968500D-02
      1.962000D+00           1.379770D-01
      4.446000D-01           4.781480D-01
      1.220000D-01           5.012400D-01
S    1   1.00
      1.220000D-01           1.000000D+00
S    1   1.00
      0.0297400              1.0000000
P    1   1.00
      7.270000D-01           1.0000000
P    1   1.00
      0.1410000              1.0000000
****

The full input file that works is:

#!/usr/bin/env python
# Author: Nike Dattani, [email protected]

'''
H in aug-cc-pV2Z
'''

import pyscf
from pyscf.gto.basis import parse_gaussian

mol = pyscf.M(
    atom = 'H 0 0 0',
    unit = 'bohr',
    basis = {
        'H' : parse_gaussian.load('avdz.gbs', 'H')
    },
    spin = 1,
    verbose = 4,
    symmetry = True,
    symmetry_subgroup = 'Dooh'
)
mf = mol.UHF()
mf.init_guess='1e'
mf.kernel()

and the output file will be:

#INFO: **** input file is /project/6004934/nike/atomic/H/av2z/pyscf/H_av2z.py ****
#!/usr/bin/env python
# Author: Nike Dattani, [email protected]

'''
H in aug-cc-pV2Z
'''

import pyscf
from pyscf.gto.basis import parse_gaussian

mol = pyscf.M(
    atom = 'H 0 0 0',
    unit = 'bohr',
    basis = {
        'H' : parse_gaussian.load('avdz.gbs', 'H')
    },
    spin = 1,
    verbose = 4,
    symmetry = True,
    symmetry_subgroup = 'Dooh'
)
mf = mol.UHF()
mf.init_guess='1e'
mf.kernel()

#INFO: ******************** input file end ********************


System: uname_result(system='Linux', node='cedar5.cedar.computecanada.ca', release='3.10.0-1160.53.1.el7.x86_64', version='#1 SMP Fri Jan 14 13:59:45 UTC 2022', machine='x86_64', processor='')  Threads 64
Python 3.8.10 (default, Jun 16 2021, 14:20:20) 
[GCC 9.3.0]
numpy 1.23.0  scipy 1.9.3
Date: Sun Dec  4 15:57:07 2022
PySCF version 2.1.0
PySCF path  /home/nike/.local/lib/python3.8/site-packages/pyscf

[CONFIG] conf_file None
[INPUT] verbose = 4
[INPUT] num. atoms = 1
[INPUT] num. electrons = 1
[INPUT] charge = 0
[INPUT] spin (= nelec alpha-beta = 2S) = 1
[INPUT] symmetry True subgroup Dooh
[INPUT] Mole.unit = bohr
[INPUT] Symbol           X                Y                Z      unit          X                Y                Z       unit  Magmom
[INPUT]  1 H      0.000000000000   0.000000000000   0.000000000000 AA    0.000000000000   0.000000000000   0.000000000000 Bohr   0.0

nuclear repulsion = 0
point group symmetry = SO3, use subgroup Dooh
symmetry origin: [0. 0. 0.]
symmetry axis x: [1. 0. 0.]
symmetry axis y: [0. 1. 0.]
symmetry axis z: [0. 0. 1.]
num. orbitals of irrep A1g = 3
num. orbitals of irrep A1u = 2
num. orbitals of irrep E1uy = 2
num. orbitals of irrep E1ux = 2
number of shells = 5
number of NR pGTOs = 12
number of NR cGTOs = 
basis = {'H': [[0, [13.01, 0.019685], [1.962, 0.137977], [0.4446, 0.478148], [0.122, 0.50124]], [0, [0.122, 1.0]], [0, [0.02974, 1.0]], [1, [0.727, 1.0]], [1, [0.141, 1.0]]]}
ecp = {}
CPU time:         1.66


******** 1 electron system ********
alpha HOMO (A1g) = -0.499334315439586  LUMO (A1g) = -0.121867194807439
converged SCF energy = -0.499334315439586  <S^2> = 0.75  2S+1 = 2
```
$\endgroup$
0

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .