1#! /usr/bin/env python3
2
3import sys
4import os
5from pyscf.pbc import scf, gto, tools
6import numpy
7import time
8import h5py
9
10alat0 = 3.6
11nks = 1
12ngs = 12
13
14cell = gto.Cell()
15cell.a = (numpy.ones((3,3))-numpy.eye(3))*alat0/2.0
16cell.atom = (('C',0,0,0),('C',numpy.array([0.25,0.25,0.25])*alat0))
17cell.basis = 'gth-szv'
18cell.pseudo = 'gth-pade'
19cell.gs = [ngs]*3  # 10 grids on postive x direction, => 21^3 grids in total
20cell.verbose = 5
21cell.build()
22
23nk = numpy.array([nks,nks,nks])
24kpt = cell.make_kpts(nk)
25mf = scf.KRHF(cell,kpt)
26mf.chkfile = "scf.dump"
27ehf = mf.kernel()
28
29from qmctools import orthoAO
30hcore = mf.get_hcore(kpts=kpt)                                   # obtain and store core hamiltonian
31fock = (hcore + mf.get_veff(kpts=kpt))                           # store fock matrix (required with orthoAO)
32LINDEP = 1e-8
33X,nmo_per_kpt = orthoAO.getOrthoAORotation(cell,kpt,LINDEP)      # store rotation to orthogonal PAO basis
34with h5py.File(mf.chkfile) as fh5:
35  fh5['scf/hcore'] = hcore
36  fh5['scf/fock'] = fock
37  fh5['scf/orthoAORot'] = X
38  fh5['scf/nmo_per_kpt'] = nmo_per_kpt
39