mirror of
https://github.com/LCPQ/quantum_package
synced 2024-12-25 13:53:49 +01:00
139 lines
4.2 KiB
Python
139 lines
4.2 KiB
Python
#!/usr/bin/env python2
|
|
|
|
import numpy,re,sys
|
|
|
|
def pyscf2QP(cell,mf, kpts=[], int_threshold = 1E-15):
|
|
# The integral will be not printed in they are bellow that
|
|
|
|
|
|
PBC=False
|
|
ComputeMode= re.split('[. ]', str(mf))
|
|
print 'ComputeMode=',ComputeMode
|
|
|
|
for n in ComputeMode:
|
|
if n in ("UHF","KUHF","UKS"):
|
|
sys.exit('Unrestricted calculation unsupported in Quantum Package')
|
|
if n == "pbc":
|
|
PBC=True
|
|
|
|
if PBC and len(kpts) == 0:
|
|
sys.exit("ERROR (read!): You need to specify explicit the list of K-point (including gamma)")
|
|
|
|
print 'Performing PBC?:',PBC
|
|
if PBC:
|
|
from pyscf.pbc import ao2mo
|
|
from pyscf.pbc import tools
|
|
else:
|
|
from pyscf import ao2mo
|
|
|
|
natom = len(cell.atom_coords())
|
|
print 'n_atom', natom
|
|
print 'num_elec', cell.nelectron
|
|
print 'nucl_num', len(cell.atom_coords())
|
|
|
|
|
|
print ''
|
|
mo_coeff = mf.mo_coeff # List of mo_coeff for each k-point
|
|
if not PBC:
|
|
nmo = mo_coeff.shape[1]
|
|
else:
|
|
nmo = mo_coeff[0].shape[1]
|
|
|
|
|
|
# Wrote all the parameter need to creat a dummy EZFIO folder who will containt the integral after.
|
|
# More an implentation detail than a real thing
|
|
with open('param','w') as f:
|
|
f.write(' '.join(map(str,(cell.nelectron, nmo, natom))))
|
|
# _
|
|
# |\ | _ | _ _. ._ |_) _ ._ | _ o _ ._
|
|
# | \| |_| (_ | (/_ (_| | | \ (/_ |_) |_| | _> | (_) | |
|
|
# |
|
|
|
|
print 'mf, cell', mf.energy_nuc(), cell.energy_nuc()
|
|
shift = tools.pbc.madelung(cell, numpy.zeros(3))*cell.nelectron * -.5 if PBC else 0
|
|
e_nuc = cell.energy_nuc() + shift
|
|
|
|
print 'nucl_repul', e_nuc
|
|
with open('e_nuc','w') as f:
|
|
f.write(str(e_nuc))
|
|
|
|
|
|
from itertools import product
|
|
|
|
# ___
|
|
# | ._ _|_ _ _ ._ _. | _ |\/| _ ._ _
|
|
# _|_ | | |_ (/_ (_| | (_| | _> | | (_) | | (_)
|
|
# _|
|
|
|
|
if PBC:
|
|
h_ao = ('kinetic', mf.get_hcore(kpts=kpts) ) # Give only one k point ?
|
|
dummy_ao = ('nuclear', numpy.zeros( (len(kpts),nmo,nmo), dtype=numpy.float ))
|
|
else:
|
|
h_ao = ('kinetic', mf.get_hcore() )
|
|
dummy_ao = ('nuclear', numpy.zeros( (nmo,nmo), dtype=numpy.float ))
|
|
|
|
def gen_mono_MO(mo_coeff,l_int,shift=0):
|
|
# 2Id transfortion Transformation. For now we handle only one or zero K point.
|
|
print 'l_int.shape=',l_int.shape
|
|
|
|
l_int_mo = reduce(numpy.dot, (mo_coeff.T, l_int, mo_coeff)) #This formula is only right for one kpt.
|
|
|
|
print 'l_int_mo=',l_int_mo
|
|
|
|
for i,j in product(range(nmo), repeat=2):
|
|
int_ = l_int_mo[i,j]
|
|
yield (i+1+shift,j+1+shift, int_)
|
|
|
|
# Print
|
|
for name, ao in (h_ao,dummy_ao):
|
|
with open('%s_mo' % name,'w') as f:
|
|
print '%s_mo' % name
|
|
if not PBC:
|
|
for mono in gen_mono_MO(mo_coeff,ao):
|
|
f.write('%s %s %s\n'% mono)
|
|
else:
|
|
for i,(m,a) in enumerate(zip(mo_coeff,ao)):
|
|
for mono in gen_mono_MO(m,a,i):
|
|
f.write('%s %s %s\n'% mono)
|
|
|
|
# ___ _
|
|
# | ._ _|_ _ _ ._ _. | _ |_) o
|
|
# _|_ | | |_ (/_ (_| | (_| | _> |_) |
|
|
# _|
|
|
#
|
|
|
|
def ao2mo_amazing(mo_coeff):
|
|
if PBC:
|
|
eri_4d= mf.with_df.ao2mo(mo_coeff,compact=False)
|
|
else:
|
|
eri_4d= ao2mo.kernel(cell,mo_coeff,compact=False)
|
|
|
|
return eri_4d.reshape((nmo,)*4)
|
|
|
|
|
|
def write_amazing(eri_4d, shift=0):
|
|
|
|
# HANDLE 8 FOLD by Scemama way. Maybe we can use compact=True
|
|
for l in range(nmo):
|
|
for k in range(nmo):
|
|
for j in range(l,nmo):
|
|
for i in range(max(j,k),nmo):
|
|
v = eri_4d[i,k,j,l]
|
|
if abs(v) > int_threshold:
|
|
f.write('%s %s %s %s %s\n' % (i+1+shift,j+1+shift,k+1+shift,l+1+shift,v))
|
|
|
|
|
|
if PBC:
|
|
eri_4d= mf.with_df.ao2mo(mo_coeff[0],compact=False)
|
|
else: #Molecular
|
|
eri_4d= ao2mo.kernel(cell,mo_coeff,compact=False)
|
|
|
|
eri_4d = eri_4d.reshape((nmo,)*4)
|
|
|
|
f = open('bielec_mo','w')
|
|
for i,mc in enumerate(mo_coeff):
|
|
eri = ao2mo_amazing(mc)
|
|
write_amazing(eri, nmo*i)
|
|
|
|
|