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