mirror of
https://github.com/LCPQ/quantum_package
synced 2024-11-14 01:53:55 +01:00
294 lines
7.3 KiB
Python
Executable File
294 lines
7.3 KiB
Python
Executable File
#!/usr/bin/env python
|
|
|
|
import sys,os
|
|
try:
|
|
QPACKAGE_ROOT = os.environ["QPACKAGE_ROOT"]
|
|
except:
|
|
print "Error: QPACKAGE_ROOT environment variable not found."
|
|
sys.exit(1)
|
|
|
|
sys.path = [ QPACKAGE_ROOT+"/EZFIO/Python", QPACKAGE_ROOT+"/resultsFile" ]+sys.path
|
|
from ezfio import ezfio
|
|
import ezfio as ez
|
|
print "EZFIO: ", os.path.dirname(ez.__file__)
|
|
|
|
try:
|
|
from resultsFile import *
|
|
except:
|
|
print "Error: resultsFile Python library not installed"
|
|
sys.exit(1)
|
|
|
|
|
|
def write_ezfioFile(res,filename):
|
|
res.clean_uncontractions()
|
|
ezfio.set_file(filename)
|
|
|
|
# Electrons
|
|
ezfio.set_electrons_elec_alpha_num(res.num_alpha)
|
|
ezfio.set_electrons_elec_beta_num(res.num_beta)
|
|
|
|
# Nuclei
|
|
ezfio.set_nuclei_nucl_num(len(res.geometry))
|
|
charge = []
|
|
coord = []
|
|
coord_x = []
|
|
coord_y = []
|
|
coord_z = []
|
|
for a in res.geometry:
|
|
charge.append(a.charge)
|
|
if res.units == 'BOHR':
|
|
coord_x.append(a.coord[0])
|
|
coord_y.append(a.coord[1])
|
|
coord_z.append(a.coord[2])
|
|
else:
|
|
coord_x.append(a.coord[0]/a0)
|
|
coord_y.append(a.coord[1]/a0)
|
|
coord_z.append(a.coord[2]/a0)
|
|
ezfio.set_nuclei_nucl_charge(charge)
|
|
label = map(lambda x: x.name, res.geometry)
|
|
ezfio.set_nuclei_nucl_label(label)
|
|
ezfio.set_nuclei_nucl_coord(coord_x+coord_y+coord_z)
|
|
|
|
# Basis
|
|
basis = res.uncontracted_basis
|
|
geom = res.geometry
|
|
|
|
# MO
|
|
MoTag = res.determinants_mo_type
|
|
ezfio.set_mo_basis_mo_label('Orthonormalized')
|
|
MO_type = MoTag
|
|
allMOs = res.uncontracted_mo_sets[MO_type]
|
|
|
|
|
|
closed = [ (allMOs[i].eigenvalue,i) for i in res.closed_mos ]
|
|
active = [ (allMOs[i].eigenvalue,i) for i in res.active_mos ]
|
|
virtual =[ (allMOs[i].eigenvalue,i) for i in res.virtual_mos ]
|
|
|
|
# closed.sort()
|
|
# active.sort()
|
|
# virtual.sort()
|
|
closed = map( lambda x: x[1], closed)
|
|
active = map( lambda x: x[1], active)
|
|
virtual = map( lambda x: x[1], virtual)
|
|
MOindices = closed + active + virtual
|
|
|
|
MOs = []
|
|
for i in MOindices:
|
|
MOs.append(allMOs[i])
|
|
|
|
mo_tot_num = len(MOs)
|
|
while len(MOindices) < mo_tot_num:
|
|
MOindices.append(len(MOindices))
|
|
|
|
MOmap = list(MOindices)
|
|
for i in range(len(MOindices)):
|
|
MOmap[i] = MOindices.index(i)
|
|
|
|
energies = []
|
|
for i in xrange(mo_tot_num):
|
|
energies.append(MOs[i].eigenvalue)
|
|
|
|
if res.occ_num is not None:
|
|
OccNum = []
|
|
for i in MOindices:
|
|
OccNum.append(res.occ_num[MO_type][i])
|
|
|
|
while len(OccNum) < mo_tot_num:
|
|
OccNum.append(0.)
|
|
|
|
MoMatrix = []
|
|
sym0 = [ i.sym for i in res.mo_sets[MO_type] ]
|
|
sym = [ i.sym for i in res.mo_sets[MO_type] ]
|
|
for i in xrange(len(sym)):
|
|
sym[MOmap[i]] = sym0[i]
|
|
|
|
MoMatrix = []
|
|
for i in xrange(len(MOs)):
|
|
m = MOs[i]
|
|
for coef in m.vector:
|
|
MoMatrix.append(coef)
|
|
|
|
while len(MoMatrix) < len(MOs[0].vector)**2:
|
|
MoMatrix.append(0.)
|
|
|
|
ezfio.set_mo_basis_mo_tot_num(mo_tot_num)
|
|
ezfio.set_mo_basis_mo_occ(OccNum)
|
|
|
|
|
|
res.clean_contractions()
|
|
# AO Basis
|
|
import string
|
|
at = []
|
|
num_prim = []
|
|
magnetic_number = []
|
|
angular_number = []
|
|
power_x = []
|
|
power_y = []
|
|
power_z = []
|
|
coefficient = []
|
|
exponent = []
|
|
res.convert_to_cartesian()
|
|
for b in res.basis:
|
|
c = b.center
|
|
for i,atom in enumerate(res.geometry):
|
|
if atom.coord == c:
|
|
at.append(i+1)
|
|
num_prim.append(len(b.prim))
|
|
s = b.sym
|
|
power_x.append( string.count(s,"x") )
|
|
power_y.append( string.count(s,"y") )
|
|
power_z.append( string.count(s,"z") )
|
|
coefficient.append( b.coef )
|
|
exponent.append( [ p.expo for p in b.prim ] )
|
|
ezfio.set_ao_basis_ao_num(len(res.basis))
|
|
ezfio.set_ao_basis_ao_nucl(at)
|
|
ezfio.set_ao_basis_ao_prim_num(num_prim)
|
|
ezfio.set_ao_basis_ao_power(power_x+power_y+power_z)
|
|
prim_num_max = ezfio.get_ao_basis_ao_prim_num_max()
|
|
len_res_basis = len(res.basis)
|
|
for i in range(len(res.basis)):
|
|
coefficient[i] += [ 0. for j in range(len(coefficient[i]),prim_num_max) ]
|
|
exponent[i] += [ 0. for j in range(len(exponent[i]),prim_num_max) ]
|
|
coefficient = reduce(lambda x, y: x+y, coefficient, [])
|
|
exponent = reduce(lambda x, y: x+y, exponent, [])
|
|
coef = []
|
|
expo = []
|
|
for i in range(prim_num_max):
|
|
for j in range(i,len(coefficient),prim_num_max):
|
|
coef.append ( coefficient[j] )
|
|
expo.append ( exponent[j] )
|
|
ezfio.set_ao_basis_ao_coef(coef)
|
|
ezfio.set_ao_basis_ao_expo(expo)
|
|
ezfio.set_ao_basis_ao_basis("Read by resultsFile")
|
|
|
|
# Apply threshold to determinants
|
|
if len(res.determinants) == 1:
|
|
sorted_determinants = [ (-1.,1.,res.determinants[0]) ]
|
|
else:
|
|
sorted_determinants = []
|
|
for i,j in zip(res.det_coefficients[0],res.determinants):
|
|
sorted_determinants.append((-abs(i),i,j))
|
|
sorted_determinants.sort()
|
|
norm = 0.0
|
|
for length, (a,b,c) in enumerate(sorted_determinants):
|
|
if -a < det_threshold:
|
|
length -=1
|
|
break
|
|
norm += a**2
|
|
norm = sqrt(norm)
|
|
length += 1
|
|
for i in xrange(length):
|
|
a = sorted_determinants[i]
|
|
sorted_determinants[i] = (a[0],a[1]/norm,a[2])
|
|
sorted_determinants = sorted_determinants[:length]
|
|
|
|
# MOs
|
|
mo_tot_num = len(res.mo_sets[MoTag])
|
|
closed_mos = res.closed_mos
|
|
active_mos = res.active_mos
|
|
virtual_mos = res.virtual_mos
|
|
to_remove = []
|
|
to_add = []
|
|
for i in active_mos:
|
|
found = False
|
|
for (a,b,c) in sorted_determinants:
|
|
if i in c['alpha']+c['beta']:
|
|
found = True
|
|
break
|
|
if not found:
|
|
to_remove.append(i)
|
|
to_add.append(i)
|
|
virtual_mos = to_add + virtual_mos
|
|
for i in active_mos:
|
|
always = True
|
|
for (a,b,c) in sorted_determinants:
|
|
if not (i in c['alpha'] and i in c['beta']):
|
|
always = False
|
|
break
|
|
if always:
|
|
to_remove.append(i)
|
|
closed_mos.append(i)
|
|
for i in to_remove:
|
|
active_mos.remove(i)
|
|
|
|
|
|
MOindices = closed_mos + active_mos + virtual_mos
|
|
while len(MOindices) < mo_tot_num:
|
|
MOindices.append(len(MOindices))
|
|
MOmap = list(MOindices)
|
|
for i in range(len(MOindices)):
|
|
MOmap[i] = MOindices.index(i)
|
|
|
|
|
|
ezfio.set_mo_basis_mo_tot_num(mo_tot_num)
|
|
|
|
mo = []
|
|
for i in MOindices:
|
|
mo.append(res.mo_sets[MoTag][i])
|
|
|
|
if len(mo) < mo_tot_num:
|
|
newmo = orbital()
|
|
newmo.eigenvalue = 0.
|
|
newmo.vector = [0. for i in range(mo_tot_num)]
|
|
newmo.vector[len(mo)] = 1.
|
|
while len(mo) < mo_tot_num:
|
|
mo.append(newmo)
|
|
Energies = [ m.eigenvalue for m in mo ]
|
|
|
|
if res.occ_num is not None:
|
|
OccNum = []
|
|
for i in MOindices:
|
|
OccNum.append(res.occ_num[MoTag][i])
|
|
|
|
while len(OccNum) < mo_tot_num:
|
|
OccNum.append(0.)
|
|
ezfio.set_mo_basis_mo_occ(OccNum)
|
|
|
|
cls = [ "v" for i in mo ]
|
|
for i in closed_mos:
|
|
cls[MOmap[i]] = 'c'
|
|
for i in active_mos:
|
|
cls[MOmap[i]] = 'a'
|
|
|
|
sym0 = [ i.sym for i in res.mo_sets[MoTag] ]
|
|
sym = [ i.sym for i in res.mo_sets[MoTag] ]
|
|
for i in xrange(len(sym)):
|
|
sym[MOmap[i]] = sym0[i]
|
|
|
|
MoMatrix = []
|
|
for m in mo:
|
|
for coef in m.vector:
|
|
MoMatrix.append(coef)
|
|
while len(MoMatrix) < len(mo[0].vector)**2:
|
|
MoMatrix.append(0.)
|
|
ezfio.set_mo_basis_mo_coef(MoMatrix)
|
|
del MoMatrix
|
|
|
|
|
|
|
|
if __name__ == '__main__':
|
|
# Check command line
|
|
|
|
det_threshold = 0.
|
|
|
|
if len(sys.argv) == 2:
|
|
State=0
|
|
elif len(sys.argv) == 3:
|
|
State=int(sys.argv[2])
|
|
else:
|
|
print "usage: "+sys.argv[0]+" file.out [state]"
|
|
sys.exit(2)
|
|
|
|
firstArg = sys.argv[1]
|
|
|
|
file = getFile(firstArg)
|
|
print firstArg, 'recognized as', str(file).split('.')[-1].split()[0]
|
|
|
|
filename = firstArg+".ezfio"
|
|
write_ezfioFile(file,filename)
|
|
|
|
|
|
|
|
|
|
|