LCPQ
/
resultsFile
mirror of https://gitlab.com/scemama/resultsFile.git
10
0
Fork 0
Code Releases Activity

Initial repository

This commit is contained in:
Anthony Scemama 2015-01-19 15:05:03 +01:00
parent 202f1fd307
commit d8a6b9f82a
16 changed files with 7371 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
Makefile Normal file
View File

@ -0,0 +1,7 @@
VERSION=1.1
default:
- rm -rf build dist
- export VERSION=$(VERSION) ;\
./setup.py --quiet bdist_rpm
- sudo rpm -e resultsFile
- sudo rpm -hiv dist/resultsFile*.noarch.rpm

49
resultsFile/Modules/__init__.py Executable file
View File

@ -0,0 +1,49 @@
#!/usr/bin/env python
# resultsFile is a library which allows to read output files of quantum
# chemistry codes and write input files.
# Copyright (C) 2007 Anthony SCEMAMA
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#
# Anthony Scemama
# LCPQ - IRSAMC
# Universite Paul Sabatier
# 118, route de Narbonne
# 31062 Toulouse Cedex 4
# scemama@irsamc.ups-tlse.fr
"""resultsFile library."""
__author__ = "Anthony SCEMAMA <scemama@irsamc.ups-tlse.fr>"
__date__ = "25 Oct 2007"
import os
wd = os.path.dirname(__file__)
all = [ i[:-3] for i in os.listdir(wd) if i.endswith("File.py") ]
if "q5CostFile" in all:
all.remove("q5CostFile")
all += [ "q5CostFile" ]
for mod in all:
try:
exec 'from '+mod+' import *'
except:
print "Error importing module", mod
pass

1565
resultsFile/Modules/gamessFile.py Executable file
View File

File diff suppressed because it is too large Load Diff

1040
resultsFile/Modules/gaussianFile.py Executable file
View File

File diff suppressed because it is too large Load Diff

39
resultsFile/Modules/include.py Executable file
View File

@ -0,0 +1,39 @@
# resultsFile is a library which allows to read output files of quantum
# chemistry codes and write input files.
# Copyright (C) 2007 Anthony SCEMAMA
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#
# Anthony Scemama
# LCPQ - IRSAMC
# Universite Paul Sabatier
# 118, route de Narbonne
# 31062 Toulouse Cedex 4
# scemama@irsamc.ups-tlse.fr
code = compile("""
try:
from resultsFile import *
except ImportError:
import os
import sys
wd = os.path.dirname(__file__)
wd = '/'.join ( wd.split('/')[:-1] )
sys.path += [ wd ]
from resultsFile import *
""", 'include.py','single')

938
resultsFile/Modules/molproFile.py Executable file
View File

@ -0,0 +1,938 @@
#!/usr/bin/python
# resultsFile is a library which allows to read output files of quantum
# chemistry codes and write input files.
# Copyright (C) 2007 Anthony SCEMAMA
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#
# Anthony Scemama
# LCPQ - IRSAMC
# Universite Paul Sabatier
# 118, route de Narbonne
# 31062 Toulouse Cedex 4
# scemama@irsamc.ups-tlse.fr
import include
eval(include.code)
import struct
import re
molproFile_defined_vars = [ "date", "version", "machine", "memory", "disk",\
"cpu_time", "author", "title", "units", "methods", \
"spin_restrict", "conv_threshs", "energies", \
"one_e_energies", "two_e_energies", \
"kin_energies", "virials", "point_group", "num_elec", \
"charge", "multiplicity","nuclear_energy","dipole","geometry",\
"basis","mo_sets","mo_types",\
"determinants", "num_alpha", "num_beta",\
"closed_mos", "active_mos", "virtual_mos", \
"determinants_mo_type", "det_coefficients", \
"csf_mo_type", "csf_coefficients", "symmetries", "occ_num", \
"csf", "num_states", "num_orb_sym" ]
methods_with_orbitals = ['RHF-SCF', 'MULTI']
class molproFile(resultsFile):
""" Class defining the molpro file.
"""
local_vars = list(local_vars)
defined_vars = list(molproFile_defined_vars)
exec get_data('date',"DATE: ",'-3:-2') in locals()
exec get_data('point_group',"Point group",'2:') in locals()
exec get_data('version'," Version",'1:2') in locals()
exec get_data('nuclear_energy',"NUCLEAR REPULSION ENERGY", '3:4',"float") in locals()
def get_num_elec(self):
if self._num_elec is None:
self._num_elec = self.num_alpha + self.num_beta
return self._num_elec
def get_charge(self):
if self._charge is None:
self._charge = 0.
for a in self.geometry:
self._charge += a.charge
self._charge -= self.num_elec
return self._charge
def get_multiplicity(self):
if self._multiplicity is None:
self._multiplicity = self.num_alpha - self.num_beta + 1
return self._multiplicity
def get_author(self):
if self._author is None:
buffer = ' '.join(self.text[-10:])
regexp = re.compile(r"user=(.+),")
matches = regexp.findall(buffer)
if matches != []:
self._author = matches[0]
return self._author
def get_machine(self):
if self._machine is None:
regexp = re.compile(r"[^/]+/([^\(]+).+DATE:")
for line in self.text:
matches = regexp.findall(line)
if matches != []:
self._machine = matches[0]
return self._machine
return self._machine
def get_disk(self):
if self._disk is None:
buffer = ' '.join(self.text[-100:])
regexp = re.compile(r"DISK USED +\* +([0-9]+\.[0-9]* +[^ \n]+)")
matches = regexp.findall(buffer)
try:
self._disk = matches[-1].title()
except IndexError:
return None
return self._disk
def get_cpu_time(self):
if self._cpu_time is None:
buffer = ' '.join(self.text[-100:])
regexp = re.compile(r"CPU TIMES +\* +([0-9]+\.[0-9]*) +[^ \n]+")
matches = regexp.findall(buffer)
try:
self._cpu_time = matches[-1] + " s"
except IndexError:
return None
return self._cpu_time
def get_memory(self):
if self._memory is None:
buffer = ' '.join(self.text[-100:])
regexp = re.compile(r"SF USED +\* +([0-9]+\.[0-9]* +[^ \n]+)")
matches = regexp.findall(buffer)
try:
self._memory = matches[-1].title()
except IndexError:
return None
return self._memory
def get_title(self):
if self._title is None:
try:
self.find_string("LABEL")
except IndexError:
return None
pos = self._pos
self._title = ' '.join(self.text[pos].split()[2:]).strip()
return self._title
def get_symmetries(self):
if self._symmetries is None:
try:
self.find_string("NUMBER OF CONTRACTIONS")
except IndexError:
return None
pos = self._pos
regexp = re.compile(r" +([0-9]+)([A-Z][^ ]*) +")
self._symmetries = [ [j,int(i)] for i,j in regexp.findall(self.text[pos]) ]
return self._symmetries
def get_units(self):
if self._units is None:
self._units = "BOHR"
return self._units
def get_methods(self):
if self._methods is None:
regexp = re.compile(r"^1PROGRAM *\* *([^ ]+)")
methods = []
for line in self.text:
program = regexp.findall(line)
if program != []:
if program[0] in methods_with_orbitals:
methods += program
self._methods = methods
return self._methods
def get_occ_num(self):
if self._occ_num is None:
occ = {}
program = []
progFound=False
doRead=False
regexp = re.compile(r"^1PROGRAM *\* *([^ ]+)")
regexp2 = re.compile(r"^ +Orb +Occ +")
regexp3 = re.compile(r"^ +[0-9]*\.[0-9]+ +([0-2]*\.[0-9]+|[0-2]|\+|\-) +[^ ]")
for line in self.text:
if not progFound:
program = regexp.findall(line)
if program != []:
if program[0] in methods_with_orbitals:
while program[0] in occ.keys():
program[0] += 'x'
occ[program[0]] = []
progFound=True
else:
buffer = line[:25]
if doRead:
if buffer.strip() != "":
matches = regexp3.findall(buffer)
if matches != []:
if matches[0] == '+':
res = 1.
else:
res = float(matches[0])
occ[program[0]].append(res)
else:
doRead = False
progFound = False
elif regexp2.match(buffer):
doRead = True
if occ != {}:
self._occ_num = occ
return self._occ_num
def get_spin_restrict(self):
if self._spin_restrict is None:
self._spin_restrict = True
for m in self.methods:
if "U" in m:
self._spin_restrict = False
return self._spin_restrict
def get_conv_threshs(self):
if self._conv_threshs is None:
self._conv_threshs = []
regexp = \
re.compile(r" conv.* thre.* +.*([0-9]+\.[0-9]*[^ ]*) *\(energy\)",
re.I)
for line in self.text:
matches = regexp.findall(line)
if matches != []:
self._conv_threshs += [float(matches[0])]
return self._conv_threshs
def get_energies(self):
if self._energies is None:
self._energies = []
self._one_e_energies = []
self._two_e_energies = []
self._ee_pot_energies = []
self._Ne_pot_energies = []
self._pot_energies = []
self._kin_energies = []
self._virials = []
pos = len(self.text)-1
while not self.text[pos].startswith(" PROGRAMS") and pos>0:
pos -=1
while not self.text[pos].startswith(" ****"):
pos += 1
pos +=1
regexp = re.compile(r"^\s+([|\-][0-9]+\.[0-9]+\s*)+")
found = False
while not found:
result = regexp.match(self.text[pos])
if result != None:
found = True
buffer = self.text[pos].split()
for e in buffer:
self._energies.insert(0,float(e))
pos += 1
program = []
progFound=False
regexp = re.compile(r"^1PROGRAM *\* *([^ ]+)")
regexp2 = re.compile(r"^ (Virial|One.elec|Two.elec|Kinetic)[^ ]* +[^ ]+ +(.+)")
progindex = -1
for line in self.text:
if not progFound:
program = regexp.findall(line)
if program != []:
if program[0] in methods_with_orbitals:
progindex += 1
progFound=True
else:
try:
quantity, value = regexp2.findall(line)[0]
if quantity == "Virial":
quantity = self._virials
progFound=False
elif quantity[0:3] == "One":
quantity = self._one_e_energies
elif quantity[0:3] == "Two":
quantity = self._two_e_energies
elif quantity == "Kinetic":
quantity = self._kin_energies
while len(quantity) <= progindex:
quantity.append(None)
quantity[progindex] = float(value)
except:
pass
return self._energies
def get_one_e_energies(self):
if self._one_e_energies is None:
self.get_energies()
return self._one_e_energies
def get_two_e_energies(self):
if self._two_e_energies is None:
self.get_energies()
return self._two_e_energies
def get_kin_energies(self):
if self._kin_energies is None:
self.get_energies()
return self._kin_energies
def get_virials(self):
if self._virials is None:
self.get_energies()
return self._virials
def get_num_alpha(self):
if self._num_alpha is None:
try:
self.find_string(" NUMBER OF ELECTRONS")
pos = self._pos
buffer = self.text[pos].split()
self._num_alpha = int(buffer[3][:-1])
self._num_beta = int(buffer[4][:-1])
except IndexError:
pass
return self._num_alpha
def get_num_beta(self):
if self._num_beta is None:
self.get_num_alpha()
return self._num_beta
def get_geometry(self):
if self._geometry is None:
try:
self.find_string(" ATOMIC COORDINATES")
pos = self._pos
pos += 4
self._geometry = []
buffer = self.text[pos].split()
while len(buffer) > 1:
temp = atom()
temp.name = buffer[1]
temp.charge = float(buffer[2])
temp.coord = (float(buffer[3]), float(buffer[4]), float(buffer[5]))
temp.basis = []
self._geometry.append(temp)
pos += 1
buffer = self.text[pos].split()
except IndexError:
pass
return self._geometry
def get_dipole(self):
if self._dipole is None:
try:
self.find_last_string("DIPOLE MOMENTS:")
pos = self._pos
self._dipole = []
buffer = self.text[pos].split(':')[1]
for d in buffer.split()[:3]:
self._dipole.append( float(d) )
except IndexError:
pass
return self._dipole
def get_num_orb_sym(self):
try:
if self._num_orb_sym is None:
self.get_basis()
except:
self._num_orb_sym = None
self.get_basis()
return self._num_orb_sym
def get_basis(self):
if self._basis is None:
try:
self.find_string("BASIS DATA")
except IndexError:
return None
pos = self._pos
try:
self.find_next_string("NUCLEAR CHARGE")
except IndexError:
return None
end = self._pos-1
num_orb_sym = []
for k in self.symmetries:
num_orb_sym.append([0,0])
pos += 4
basis_read = []
index, sym, nuc, type, expo, coef = "", "", "", "", "", ""
basis = []
currfunc = 0
salcao2mo = []
contrmax=0
buffer = self.text[pos:end+1]
text = []
for line in buffer:
if line[:41].strip() == "":
text[-1] = text[-1][:-1] + line[41:]
else:
text.append(line)
pos=0
end = len(text)
text.append(80*" ")
while pos < end:
begin = pos
line = text[begin]
nfunc = len( line[41:].split() )
nuc = []
try:
idx4 = line.find('.')
idx5 = idx4+1
idx6 = idx4+2
currfunc = int(line[:idx4])
geom = int( line[idx5:idx6] )
except:
line = str(currfunc)+'.'+str(geom)+line[8:]
end2 = pos
while end2<end and int( line[:idx4] ) < currfunc+nfunc and \
int( line[idx5:idx6] ) == geom:
if line[11:15].strip() != '':
nuc.append( (int(line[11:15])-1, line[15:22].strip() ) )
end2 += 1
line = text[end2]
if line[:idx4] == ' ':
line = "%4d.%d "%(currfunc,geom)+line[8:]
line = text[begin]
contr = []
for k in range(nfunc):
for NucType in nuc:
c = contraction()
contr.append(c)
thiscontr=0
if nfunc > contrmax: contrmax = nfunc
for k in range(nfunc):
smo = []
num_orb_sym[geom-1][0] += 1
num_orb_sym[geom-1][1] += len(nuc)
for NucType in nuc:
c = contr[thiscontr]
thiscontr += 1
sign=1.
pos = begin
line = text[pos]
atom = self.geometry[NucType[0]]
type = str(NucType[1])
if type[0] == '-':
type = type[1:]
sign = -1.
if type[0].isdigit():
type = type[1:]
try:
if not type[1].isdigit() :
if type[0].lower() in "pdfghi": type = type[1:]
except:
if type[0].lower() in "pdfghi": type = type[1:]
smo.append(sign/sqrt(float(len(nuc))))
c.center = atom.coord
type = normalize_basis_name(type)
c.sym = type
while pos<end2:
expo = line[24:40].strip()
coef = line[41:].split()
if expo == '':
pos = end2
break
gauss = gaussian()
gauss.center = atom.coord
gauss.expo = float(expo)
gauss.sym = type
if float(coef[k]) != 0.:
c.append(float(coef[k]),gauss)
pos += 1
line = text[pos]
salcao2mo.append( smo )
for c in contr:
# c.normalize()
basis.append(c)
self._num_orb_sym = num_orb_sym
self._basis = basis
self.salcao2mo = salcao2mo
for f in basis:
for at in self.geometry:
if f.center is at.coord:
at.basis.append(f)
if (len(self.symmetries)>1) and contrmax > 1:
file = open('basis.molpro','w')
molpro_write_input(self,file)
file.close()
print >>sys.stderr, """
Warning: there is a bug in the molpro output. Run the 'basis.molpro'
input file with molpro to generate a correct output for the basis set.
"""
return self._basis
def get_mo_types(self):
if self._mo_types is None:
self.get_mo_sets()
return self._mo_types
def get_mo_sets(self):
if self._mo_sets is None:
global methods_with_orbitals
self._mo_sets = {}
self._mo_types = []
doRead=False
regexp = re.compile(r"^1PROGRAM *\* *([^ ]+)")
regexp2 = re.compile(r"^ +Orb +Occ +")
rstring = r"^ {,4}([0-9]+\.[0-9]+) +" # Orb
rstring += r"([0-2| ]\.[0-9]{,5}|[0-2]|\+|\-) *" # Occ
rstring += r"(-*[0-9|-]*\.[0-9]+) +" # Energy
regexp3 = re.compile(rstring)
rstring = r"(([0-9|\-]*\.[0-9]{6} *)+)" # Coefs
regexp4 = re.compile(rstring)
rstring = r"(^ {30,}([0-9] [123456789xyzspdfghijk\+\-]+ *)+)" # Labels
regexp5 = re.compile(rstring)
rstring = r"(^ \*{10,}|orbital dump)" # *********
regexp6 = re.compile(rstring)
rstring = " +[0-9|-]*\.[0-9]+\.[0-9]*"
regexp7 = re.compile(rstring)
text = list(self.text)
curr_sym = self.basis[0].sym
v=None
ReadingLabels = False
PreviousBlank = False
for iline, line in enumerate(text):
buffer = regexp7.findall(line)
while buffer != []: # Insert spaces in "7851.548715699.7413"
if ":" in line:
break
begin = line.find(buffer[0])
end = begin+len(buffer[0])
line = line[0:end-10]+" "+line[end-10:]
buffer = regexp7.findall(line)
line = line[:11]+line[11:32].replace ("**********"," 999999.9999") \
+line[32:].replace('-',' -')
# Find method
program = regexp.findall(line)
if program != []:
if program[0] in methods_with_orbitals:
add_method = False
while program[0] in self._mo_sets:
program[0] += 'x'
add_method = True
if add_method and program[0] not in methods_with_orbitals:
methods_with_orbitals += [program[0]]
index = program[0]
self._mo_types.append(index)
self._mo_sets[index] = []
if doRead:
# Find labels
buffer = regexp6.findall(line)
if buffer != []:
doRead=False
elif line.strip() == "":
if symcount == len(self.symmetries):
if PreviousBlank:
doRead = False
PreviousBlank = False
else:
if not ReadingLabels:
PreviousBlank = True
if v is not None:
for l,s in enumerate(self.num_orb_sym):
if l > symcount-1 :
for k in range(s[1]):
v.vector.append(0.)
self._mo_sets[index].append(v)
v = None
elif line.strip() != "":
buffer = regexp5.findall(line)
PreviousBlank = False
if buffer != []:
if not ReadingLabels:
ReadingLabels = True
else:
ReadingLabels = False
# Find "Orb Occ Energy" values to initiate a new orbital
buffer = regexp3.findall(line)
if buffer != []:
buffer = buffer[0]
symcount = int(buffer[0].split('.')[1])
bf = 0
v = orbital()
v.set = index
v.basis = self.basis
v.sym = self.symmetries[ int(buffer[0].split('.')[1])-1 ][0]
v.eigenvalue = float(buffer[2])
if v.eigenvalue == 999999.9999:
print "Warning line", iline+1, ": '**********' in orbital eigenvalues"
for l,s in enumerate(self.num_orb_sym):
if l < symcount-1 :
bf += s[0]
for k in range(0,s[1]):
v.vector.append(0.)
# Find coefficient values to continue the orbital
buffer = regexp7.findall(line[30:])
while buffer != []: # Insert spaces in "7851.548715699.7413"
begin = line[0:].find(buffer[0])
end = begin+len(buffer[0])
line = line[0:end-10]+" "+line[end-10:]
buffer = regexp7.findall(line)
buffer = regexp4.findall(line[30:])
if buffer != []:
for x in buffer[0][0].split():
for k in self.salcao2mo[bf]:
try:
x2 = k*float(x)
except ValueError:
print "Error line", iline+1, ": orbital coefficients"
sys.exit(1)
v.vector.append(x2)
bf+=1
elif regexp2.match(line):
if self._mo_sets[index] == []:
doRead = True
symcount=0
return self._mo_sets
def get_determinants_mo_type(self):
if self._determinants_mo_type is None:
self._determinants_mo_type = self.mo_types[-1]
return self._determinants_mo_type
def get_csf_mo_type(self):
if self._csf_mo_type is None:
self._csf_mo_type = self.determinants_mo_type
return self._determinants_mo_type
def get_determinants(self):
if self._determinants is None:
determinants = []
if self.csf is not None:
for csf in self.csf:
for new_det in csf.determinants:
determinants.append(new_det)
else:
pass
if determinants != []:
self._determinants_mo_type = self.mo_types[-1]
self._determinants = determinants
return self._determinants
def get_csf(self):
if self._csf is None:
csf = []
csfcoef = [[]]
index = len(self.methods)-1
while self.mo_types[index] not in methods_with_orbitals and index >0:
index -= 1
# Mono-determinant case
if self.mo_types[index].startswith("RHF-SCF"):
new_csf = CSF()
new_spin_det_a = []
new_spin_det_b = []
occ = self.occ_num[self.determinants_mo_type]
for i,v in enumerate(occ):
if v == 1.:
new_spin_det_a.append(i)
elif v == 2.:
new_spin_det_a.append(i)
new_spin_det_b.append(i)
else:
assert v == 0.
new_csf.append(1.,new_spin_det_a,new_spin_det_b)
csf.append(new_csf)
csfcoef[0].append(1.)
self._determinants_mo_type = self.mo_types[index]
# Multi-determinant case
elif self.mo_types[index].startswith("MULTI"):
try:
self.find_last_string('Number of closed-shell orbitals')
pos = self._pos
buffer = self.text[pos].split()
ncore0 = [ int(a) for a in buffer[6:-1] ]
except IndexError:
ncore0 = [ 0 for a in self.symmetries ]
self.find_last_string('Number of active')
pos = self._pos
buffer = self.text[pos].split()
nact0 = [ int(a) for a in buffer[6:-1] ]
self.find_next_string('Number of external')
mo_count = [ 0 for i in self.num_orb_sym ]
j = 0
# curr_sym = self.mo_sets[self.determinants_mo_type][0].sym[0]
curr_sym = self.mo_sets[self.determinants_mo_type][0].sym
for i in self.mo_sets[self.determinants_mo_type]:
if i.sym is not None:
if i.sym != curr_sym:
curr_sym = i.sym
# if i.sym[0] != curr_sym:
# curr_sym = i.sym[0]
j+=1
mo_count[j] += 1
next0 = [ mo_count[i] - nact0[i] - ncore0[i] for i in range(len(self.num_orb_sym)) ]
self.find_next_string('TOTAL ENERGIES')
end = self._pos -2
self.find_prev_string('CI vector')
pos = self._pos +2
while pos < end:
pos += 1
this_csf = CSF()
tempcsf_a = []
tempcsf_b = []
buffer = [ ''.join(self.text[pos].split()[:-1]),self.text[pos].split()[-1]]
mo, coef = buffer
coef = float(coef)
ncoreold=0
nactold=0
mo2=""
for ((ncore,nact),next) in zip(zip(ncore0,nact0),next0):
for i in range(ncore):
tempcsf_a.append(ncoreold+i)
tempcsf_b.append(ncoreold+i)
mo2+="2"
ncoreold += ncore
for k,i in enumerate(mo[nactold:nactold+nact]):
if i in "2a+/":
tempcsf_a.append(ncoreold+k)
if i in "2b-\\":
tempcsf_b.append(ncoreold+k)
mo2 += i
ncoreold += nact+next
nactold += nact
# Correction de la phase due au changement d'ordre des orbitales
p_count=0
for k,i in enumerate(mo2):
if i in "2b-\\":
for j in mo2[k+1:]:
if j in "2a+/":
p_count+=1
this_csf.append((-1.0)**p_count,tempcsf_a,tempcsf_b)
csf.append(this_csf)
csfcoef[0].append(coef)
if csf != []:
self._csf = csf
self._csf_coefficients = csfcoef
return self._csf
def get_csf_coefficients(self):
if self._csf_coefficients is None:
self.get_csf()
return self._csf_coefficients
def get_closed_mos(self):
if self._closed_mos is None:
result = []
occ = filter(lambda x: x>0., self.occ_num[self.determinants_mo_type])
maxmo = len(occ)
for orb in range(maxmo):
present = True
for det in self.determinants:
if present:
for spin_det in det:
if orb not in det[spin_det]: present = False
else:
break
if present: result.append(orb)
self._closed_mos = result
return self._closed_mos
def get_virtual_mos(self):
if self._virtual_mos is None:
result = []
occ = filter(lambda x: x>0., self.occ_num[self.determinants_mo_type])
maxmo = len(self.mo_sets[self.determinants_mo_type])
for orb in range(maxmo):
present = False
for det in self.determinants:
if not present:
for spin_det in det:
if orb in det[spin_det]: present = True
else:
break
if not present: result.append(orb)
self._virtual_mos = result
self._virtual_mos = result
return self._virtual_mos
def get_active_mos(self):
if self._active_mos is None:
cl = self.closed_mos
vi = self.virtual_mos
maxmo = len(self.mo_sets[self.determinants_mo_type])
result = []
for i in range(maxmo):
present = i in cl or i in vi
if not present:
result.append(i)
self._active_mos = result
return self._active_mos
def get_det_coefficients(self):
if self._det_coefficients is None:
index = len(self.methods)-1
while self.mo_types[index] not in methods_with_orbitals and index >0:
index -= 1
# Mono-determinant case
if self.mo_types[index].startswith("RHF-SCF"):
self._det_coefficients = [ [1.] ]
# Multi-determinant case
if self.mo_types[index].startswith("MULTI"):
if self.csf is not None:
self._det_coefficients = []
csf = self.csf
for state_coef in self.csf_coefficients:
vector = []
for i,c in enumerate(state_coef):
for d in csf[i].coefficients:
vector.append(c*d)
self._det_coefficients.append(vector)
return self._det_coefficients
def get_num_states(self):
if self._num_states is None:
self._num_states=1
return self._num_states
# Properties
# ----------
exec build_property("num_orb_sym","Number of SALCAO/sym") in locals()
to_remove = []
for i, j in local_vars:
if i in resultsFile_defined_vars:
to_remove.append( (i,j) )
for i in to_remove:
local_vars.remove(i)
for i, j in local_vars:
if i not in defined_vars:
exec build_get_funcs(i) in locals()
exec build_property(i,j) in locals()
del to_remove, i, j
# Output Routines
# ---------------
def molpro_write_input(res,file):
print >>file, " ***,", res.title
geom = res.geometry
print >>file, """
print,basis;
print,orbitals;
gprint,civector;
gthresh,printci=0.;
geomtyp=xyz
geometry={
""", len(geom)
print >>file, res.title
i=1
for at in geom:
coord = []
for x in at.coord: coord.append(x*a0)
print >>file, "%10s %15.8f %15.8f %15.8f" % tuple( \
[(at.name+str(i)).ljust(10) ]+coord )
i+=1
print >>file, "}\n"
print >>file, "basis={"
lines = []
for idx,at in enumerate(geom):
# Find atom
doPrint = True
for contr in at.basis:
line = ""
# Find Sym
sym = contr.sym
if sym == 's' :
sym = 's'
doPrint = True
elif sym == 'x' :
sym = 'p'
doPrint = True
elif sym == 'd+0' :
sym = 'd'
doPrint = True
elif sym == 'xx' :
sym = 'd'
doPrint = True
elif sym == 'f+0' :
sym = 'f'
doPrint = True
elif sym == 'xxx' :
sym = 'f'
doPrint = True
elif sym == 'g+0' :
sym = 'g'
doPrint = True
elif sym == 'xxxx' :
sym = 'g'
doPrint = True
elif sym == 'h+0' :
sym = 'h'
doPrint = True
elif sym == 'xxxxx' :
sym = 'h'
doPrint = True
elif sym == 'i+0' :
sym = 'i'
doPrint = True
elif sym == 'xxxxxx' :
sym = 'i'
doPrint = True
elif sym == 'xxxxxxx' :
sym = 'j'
doPrint = True
else:
doPrint = False
if doPrint:
# Find exponents
line += "%1s,%s%d"%(sym, at.name, idx+1)
for i, p in enumerate(contr.prim):
if not isinstance(p,gaussian):
raise TypeError("Basis function is not a gaussian")
line += ",%f"%(float(p.expo))
line += "\n"
# Find coefs
line += "c,1."+str(len(contr.prim))
for c in contr.coef:
line += ",%g"%(c)
if line not in lines:
lines.append(line)
for line in lines:
print >>file, line
print >>file, "}\n",
print >>file, "rhf;"
print >>file, "---"
fileTypes.append(molproFile)
if __name__ == '__main__':
main(molproFile)

1357
resultsFile/Modules/qmcchemFile.py Executable file
View File

File diff suppressed because it is too large Load Diff

297
resultsFile/Modules/qmcchem_newFile.py Executable file
View File

@ -0,0 +1,297 @@
#!/usr/bin/python
# resultsFile is a library which allows to read output files of quantum
# chemistry codes and write input files.
# Copyright (C) 2007 Anthony SCEMAMA
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#
# Anthony Scemama
# LCPQ - IRSAMC
# Universite Paul Sabatier
# 118, route de Narbonne
# 31062 Toulouse Cedex 4
# scemama@irsamc.ups-tlse.fr
import include
eval(include.code)
import struct
import re
import os
QMCCHEM_PATH = os.getenv("QMCCHEM_PATH",default="")
if QMCCHEM_PATH == "":
print "QmcChem new files are not handled."
class qmcchem_newFile(resultsFile):
pass
else:
sys.path = [ QMCCHEM_PATH+"/scripts" ]+sys.path
from ezfio import ezfio
qmcchem_newFile_defined_vars = [ "date", "version", \
"title", "units", "methods", \
"point_group", "num_elec", \
"charge", "multiplicity","geometry",\
"basis","mo_sets","mo_types",\
"determinants", "num_alpha", "num_beta",\
"closed_mos", "active_mos", "virtual_mos", \
"determinants_mo_type", "det_coefficients", \
"csf_mo_type", "csf_coefficients", "occ_num", \
"csf" ]
class qmcchem_newFile(resultsFile):
""" Class defining the qmcchem_new file.
"""
local_vars = list(local_vars)
defined_vars = list(qmcchem_newFile_defined_vars)
def __init__(self,name):
resultsFile.__init__(self,name)
ezfio.set_filename(self.filename)
def get_version(self):
if self._version is None:
self._version = ezfio.get_version()
return self._version
def get_date(self):
if self._date is None:
self._date = ezfio.get_ezfio_creation()
return self._date
def get_num_elec(self):
if self._num_elec is None:
self._num_elec = self.num_alpha + self.num_beta
return self._num_elec
def get_multiplicity(self):
if self._multiplicity is None:
self._multiplicity = self.num_alpha - self.num_beta + 1
return self._multiplicity
def get_charge(self):
if self._charge is None:
self._charge = sum(ezfio.get_nuclei_nucl_charge())-float(self.num_elec)
return self._charge
def get_title(self):
if self._title is None:
self._title = self.filename
return self._title
def get_units(self):
if self._units is None:
self._units = 'BOHR'
#for a gamess use the units is give by an option on gamess_write_contrl
return self._units
def get_methods(self):
if self._methods is None:
self._methods = ['QMC']
return self._methods
def get_point_group(self):
if self._point_group is None:
self._point_group = "C1"
return self._point_group
def get_geometry(self):
if self._geometry is None:
self.get_geometryX()
self.get_basisX()
return self._geometry
def get_basis(self):
if self._basis is None:
self.get_geometry()
return self._basis
def get_geometryX(self):
self._geometry = []
charge = ezfio.get_nuclei_nucl_charge()
coord = ezfio.get_nuclei_nucl_coord()
num = ezfio.get_nuclei_nucl_num()
for i in range(num):
temp = atom()
temp.charge = charge[i]
temp.coord = (coord[0][i], coord[1][i], coord[2][i])
temp.name = 'X'
temp.basis = []
self._geometry.append(temp)
def get_basisX(self):
coef = ezfio.get_ao_basis_ao_coef()
expo = ezfio.get_ao_basis_ao_expo()
nucl = ezfio.get_ao_basis_ao_nucl()
num = ezfio.get_ao_basis_ao_num()
power= ezfio.get_ao_basis_ao_power()
prim_num= ezfio.get_ao_basis_ao_prim_num()
self._basis = []
for i in range(num):
contr = contraction()
for j in range(prim_num[i]):
gauss = gaussian()
atom = self._geometry[nucl[i]-1]
gauss.center = atom.coord
gauss.expo = expo[j][i]
name = normalize_basis_name('x'*power[0][i]+'y'*power[1][i]+'z'*power[2][i])
if name == '': name = 's'
gauss.sym = name
contr.append(coef[j][i],gauss)
self._geometry[nucl[i]-1].basis.append(contr)
self._basis.append(contr)
def get_mo_types(self):
if self._mo_types is None:
self._mo_types = ['QMC']
return self._mo_types
def get_mo_sets(self):
if self._mo_sets is None:
self._mo_sets = {}
self._mo_sets['QMC'] = []
coef = ezfio.get_mo_basis_mo_coef()
energy = ezfio.get_mo_basis_mo_energy()
num = ezfio.get_mo_basis_mo_tot_num()
for i in range(num):
v = orbital()
v.basis = self.basis
v.set = 'QMC'
v.eigenvalue = energy[i]
v.vector = coef[i]
self.mo_sets['QMC'].append(v)
return self._mo_sets
def get_num_alpha(self):
if self._num_alpha is None:
self._num_alpha = ezfio.get_electrons_elec_alpha_num()
return self._num_alpha
def get_num_beta(self):
if self._num_beta is None:
self._num_beta = ezfio.get_electrons_elec_beta_num()
return self._num_beta
def get_determinants_mo_type(self):
if self._determinants_mo_type is None:
self._determinants_mo_type = 'QMC'
return self._determinants_mo_type
def get_csf_mo_type(self):
if self._csf_mo_type is None:
self._csf_mo_type = 'QMC'
return self._csf_mo_type
def get_determinants(self):
if self._determinants is None:
determinants = []
if self.csf is not None:
for csf in self.csf:
for new_det in csf.determinants:
determinants.append(new_det)
else:
pass
if determinants != []:
self._determinants_mo_type = self.mo_types[-1]
self._determinants = determinants
return self._determinants
def get_csf(self):
method = self.methods[0]
if self._csf is None:
csf = []
ncore = ezfio.get_mo_basis_mo_closed_num()
nact = ezfio.get_mo_basis_mo_active_num()
core_a = []
core_b = []
for i in range(ncore):
core_a.append(i)
core_b.append(i)
num = ezfio.get_determinants_det_num()
occ = ezfio.get_determinants_det_occ()
if occ == []:
occ = [[[0]],[[0]]]
for i in range(num):
this_csf = CSF()
tempcsf_a = core_a + map(lambda x: x-1, occ[0][i])
tempcsf_b = core_b + map(lambda x: x-1, occ[1][i])
this_csf.append(1.,tempcsf_a,tempcsf_b)
csf.append(this_csf)
if csf != []:
self._csf = csf
return self._csf
def get_closed_mos(self):
if self._closed_mos is None:
cls = ezfio.get_mo_basis_mo_classif()
self._closed_mos = []
self._virtual_mos = []
self._active_mos = []
for i in range(len(cls)):
if cls[i] == 'c':
self._closed_mos.append(i)
elif cls[i] == 'a':
self._active_mos.append(i)
elif cls[i] == 'v':
self._virtual_mos.append(i)
return self._closed_mos
def get_virtual_mos(self):
if self._virtual_mos is None:
self.get_closed_mos()
return self._virtual_mos
def get_active_mos(self):
if self._active_mos is None:
self.get_closed_mos()
return self._active_mos
def get_det_coefficients(self):
if self._det_coefficients is None:
self._det_coefficients = []
csf = self.csf
for state_coef in self.csf_coefficients:
vector = []
for i,c in enumerate(state_coef):
for d in csf[i].coefficients:
vector.append(c*d)
self._det_coefficients.append(vector)
return self._det_coefficients
def get_csf_coefficients(self):
if self._csf_coefficients is None:
self._csf_coefficients = [ [] ]
self._csf_coefficients[0] = ezfio.get_determinants_det_coef()
return self._csf_coefficients
def get_occ_num(self):
if self._occ_num is None:
self._occ_num = {}
motype = 'QMC'
self._occ_num[motype] = ezfio.get_mo_basis_mo_occ()
return self._occ_num
fileTypes.append(qmcchem_newFile)
if __name__ == '__main__':
main(qmcchem_newFile)

487
resultsFile/Modules/wfnFile.py Executable file
View File

@ -0,0 +1,487 @@
#!/usr/bin/python
# resultsFile is a library which allows to read output files of quantum
# chemistry codes and write input files.
# Copyright (C) 2007 Anthony SCEMAMA
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#
# Anthony Scemama
# LCPQ - IRSAMC
# Universite Paul Sabatier
# 118, route de Narbonne
# 31062 Toulouse Cedex 4
# scemama@irsamc.ups-tlse.fr
import include
eval(include.code)
import struct
wfnFile_defined_vars = [
"energies", "title", "units",\
"virials", "num_elec", \
"charge", "geometry",\
"basis","mo_sets","mo_types", \
"num_alpha", "num_beta",\
"closed_mos", "active_mos", "virtual_mos", \
"determinants",
"determinants_mo_type", "det_coefficients", \
"csf_mo_type", "csf_coefficients", "occ_num", \
"csf", "num_states"]
class wfnFile(resultsFile):
""" Class defining the wfn file.
"""
local_vars = list(local_vars)
defined_vars = list(wfnFile_defined_vars)
def get_title(self):
if self._title is None:
self._title = self.text[0].strip()
return self._title
def get_units(self):
if self._units is None:
self._units = "BOHR"
return self._units
def get_num_elec(self):
if self._num_elec is None:
self._num_elec = int(sum(self.occ_num["Wfn"]))
return self._num_elec
def get_occ_num(self):
if self._occ_num is None:
result = {}
occ = []
for line in self.text:
if line.startswith("MO"):
buffer = line.split()
while buffer[0] != "OCC":
buffer.pop(0)
occ.append(float(buffer[3]))
result["Wfn"] = occ
self._occ_num = result
return self._occ_num
def get_energies(self):
if self._energies is None:
buffer = self.text[-1]
self._energies = [float(buffer.split("ENERGY")[1].split()[1])]
self._virials = [float(buffer.split("VIRIAL(-V/T)=")[1].split()[0])]
return self._energies
def get_virials(self):
if self._virials is None:
self.get_energies()
return self._virials
def get_charge(self):
if self._charge is None:
result = 0.
for line in self.text:
if "CHARGE" in line:
result += float(line.split("CHARGE =")[1].strip())
result -= self.num_elec
self._charge = result
return self._charge
def get_geometry(self):
if self._geometry is None:
Natoms = int(self.text[1].split()[6])
pos = 2
begin = pos
self._geometry = []
buffer = self.text[pos].split()
while pos<begin+Natoms:
temp = atom()
temp.name = buffer[0]
temp.charge = float(buffer[9])
temp.coord = (float(buffer[4]), float(buffer[5]), float(buffer[6]))
temp.basis = []
self._geometry.append(temp)
pos += 1
buffer = self.text[pos].split()
b = self.basis
if b is not None:
try:
for f in b:
for at in self._geometry:
if f.center is at.coord:
at.basis.append(f)
except IndexError:
pass
return self._geometry
def get_basis(self):
if self._basis is None:
center = []
type = []
expo = []
conversion = [ "S", "X", "Y", "Z" ]
for line in self.text:
if line.startswith('CENTRE'):
buffer = line[20:].split()
for b in buffer:
center.append(int(b)-1)
elif line.startswith('TYPE'):
buffer = line[20:].split()
for b in buffer:
type.append(int(b))
elif line.startswith('EXPONENTS'):
buffer = line[10:].split()
for b in buffer:
expo.append(float(b.replace('D','e')))
if 10 in type: # Cartesian d orbitals
conversion += [ "XX", "YY", "ZZ", "XY", "XZ", "YZ" ]
else: # Spherical d
conversion += [ "D2-", "D1-", "D0", "D1+", "D2+", "" ]
if 20 in type: # Cartesian f orbitals
conversion += [ "XXX", "YYY", "ZZZ", "YYX", "XXY", "XXZ",
"ZZX", "ZZY", "YYZ", "XYZ" ]
else: # Spherical f
conversion += [ "F3-", "F2-", "F1-", "F0", "F1+", "F2+",
"F3+", "", "", "" ]
for i in range(len(type)):
type[i] = conversion[type[i]-1]
assert ( len(expo) == len(type) )
assert ( len(expo) == len(center) )
self._basis = []
for i in range(len(expo)):
contr = contraction()
gauss = gaussian()
atom = self.geometry[center[i]]
gauss.center = atom.coord
gauss.expo = expo[i]
gauss.sym = type[i]
contr.append(1.,gauss)
self._basis.append(contr)
return self._basis
def get_mo_types(self):
if self._mo_types is None:
self.get_mo_sets()
return self._mo_types
def get_mo_sets(self):
if self._mo_sets is None:
self._mo_sets = {}
index = 'Wfn'
self._mo_types = [index]
posend = {}
k=0
vectors = []
for line in self.text:
if line.startswith('MO'):
if k>0:
vectors.append(v)
v = orbital()
v.set = index
v.basis = self.basis
v.eigenvalue = float(line.split()[-1])
v.sym = None
k+=1
elif line.startswith('END'):
vectors.append(v)
k=0
elif k>0:
for c in line.split():
v.vector.append(float(c.replace('D','e')))
basis = self.basis
for v in vectors:
for i in range(len(basis)):
v.vector[i] *= basis[i].norm
self._mo_sets[index] = vectors
return self._mo_sets
def get_num_alpha(self):
if self._num_alpha is None:
self._num_alpha = self.num_elec-self.num_beta
return self._num_alpha
def get_num_beta(self):
if self._num_beta is None:
self._num_beta = self.num_elec/2
return self._num_beta
def get_determinants_mo_type(self):
if self._determinants_mo_type is None:
self._determinants_mo_type = self.mo_types[-1]
return self._determinants_mo_type
def get_csf_mo_type(self):
if self._csf_mo_type is None:
self._csf_mo_type = self.determinants_mo_type
return self._csf_mo_type
def get_determinants(self):
if self._determinants is None:
determinants = []
if self.csf is not None:
for csf in self.csf:
for new_det in csf.determinants:
determinants.append(new_det)
else:
pass
if determinants != []:
self._determinants_mo_type = self.mo_types[-1]
self._determinants = determinants
return self._determinants
def get_csf(self):
if self._csf is None:
csf = []
new_csf = CSF()
new_spin_det_a = []
new_spin_det_b = []
ntot = 0
for i,n in enumerate(self.occ_num['Wfn']):
if n == 1.0:
if ntot < self.num_alpha:
new_spin_det_a.append(i)
else:
new_spin_det_b.append(i)
ntot += 1
elif n == 2.0:
new_spin_det_a.append(i)
new_spin_det_b.append(i)
new_csf.append(1.,new_spin_det_a,new_spin_det_b)
csf.append(new_csf)
self._determinants_mo_type = self.mo_types[-1]
if csf != []:
self._csf = csf
return self._csf
def get_closed_mos(self):
if self._closed_mos is None:
result = []
maxmo = len(self.mo_sets[self.determinants_mo_type])
for orb in range(maxmo):
present = True
for det in self.determinants:
for spin_det in det:
if orb not in det[spin_det]: present = False
if present: result.append(orb)
self._closed_mos = result
return self._closed_mos
def get_virtual_mos(self):
if self._virtual_mos is None:
result = []
minmo = len(self.closed_mos)
maxmo = len(self.mo_sets[self.determinants_mo_type])
for orb in range(minmo,maxmo):
present = False
for det in self.determinants:
for spin_det in det:
if orb in det[spin_det]: present = True
if not present: result.append(orb)
self._virtual_mos = result
return self._virtual_mos
def get_active_mos(self):
if self._active_mos is None:
cl = self.closed_mos
vi = self.virtual_mos
maxmo = len(self.mo_sets[self.determinants_mo_type])
result = []
for i in range(maxmo):
present = i in cl or i in vi
if not present:
result.append(i)
self._active_mos = result
return self._active_mos
def get_det_coefficients(self):
if self._det_coefficients is None:
self._det_coefficients = [ [1.] ]
return self._det_coefficients
def get_csf_coefficients(self):
if self._csf_coefficients is None:
self._csf_coefficients = [ [1.] ]
return self._csf_coefficients
def get_num_states(self):
if self._num_states is None:
self._num_states=1
return self._num_states
# Properties
# ----------
to_remove = []
for i, j in local_vars:
if i in resultsFile_defined_vars:
to_remove.append( (i,j) )
for i in to_remove:
local_vars.remove(i)
for i, j in local_vars:
if i not in defined_vars:
exec build_get_funcs(i) in locals()
exec build_property(i,j) in locals()
del to_remove, i, j
# Output Routines
# ---------------
import string
basis_correspond = {}
basis_correspond['s' ] = 1
basis_correspond['x' ] = 2
basis_correspond['y' ] = 3
basis_correspond['z' ] = 4
basis_correspond['xx' ] = 5
basis_correspond['yy' ] = 6
basis_correspond['zz' ] = 7
basis_correspond['xy' ] = 8
basis_correspond['xz' ] = 9
basis_correspond['yz' ] = 10
basis_correspond['xxx'] = 11
basis_correspond['yyy'] = 12
basis_correspond['zzz'] = 13
basis_correspond['xyy'] = 14
basis_correspond['xxy'] = 15
basis_correspond['xxz'] = 16
basis_correspond['xzz'] = 17
basis_correspond['yzz'] = 18
basis_correspond['yyz'] = 19
basis_correspond['xyz'] = 20
basis_correspond['d-2'] = 5
basis_correspond['d-1'] = 6
basis_correspond['d+0'] = 7
basis_correspond['d+1'] = 8
basis_correspond['d+2'] = 9
basis_correspond['f-3'] = 11
basis_correspond['f-2'] = 12
basis_correspond['f-1'] = 13
basis_correspond['f+0'] = 14
basis_correspond['f+1'] = 15
basis_correspond['f+2'] = 16
basis_correspond['f+3'] = 17
def wfn_write(res,file,MO_type=None):
print >>file, " "+res.title.strip()
print >>file, "GAUSSIAN %14d MOL ORBITALS %6d PRIMITIVES %8d NUCLEI" \
%(len(res.closed_mos+res.active_mos), len(res.uncontracted_basis), len(res.geometry)),
# write geom
geom = res.geometry
f=1.
if res.units == 'ANGS':
f=1./a0
for i,atom in enumerate(geom):
name = atom.name
for d in string.digits:
name = name.replace(d,'')
print >>file, "\n%3s %4d (CENTRE%3d) %12.8f%12.8f%12.8f CHARGE = %4.1f" \
%(name, i+1, i+1, atom.coord[0]*f, atom.coord[1]*f, atom.coord[2]*f, atom.charge),
# write basis
basis = res.uncontracted_basis
center = []
type = []
expo = []
for prim in basis:
c = prim.center
for i,atom in enumerate(geom):
atom_c = atom.coord
if atom_c == c:
atom_id = i+1
sym = basis_correspond[prim.sym]
center.append(atom_id)
type.append(sym)
expo.append(prim.expo)
k=0
for c in center:
if k%20 == 0:
print "\nCENTRE ASSIGNMENTS ",
print "%2d"%(c),
k+=1
k=0
for t in type:
if k%20 == 0:
print "\nTYPE ASSIGNMENTS ",
print "%2d"%(t),
k+=1
k=0
for e in expo:
if k%5 == 0:
print "\nEXPONENTS ",
print ("%13.7e"%(e)).replace('e','D'),
k+=1
print '\n',
# MOs
if MO_type is None:
MO_type = res.determinants_mo_type
try:
allMOs = res.uncontracted_mo_sets[MO_type]
except KeyError:
print "MO type not present in "+res.filename
return
mos = []
occ = res.occ_num[MO_type]
for i,mo in enumerate(allMOs):
mos.append( (mo.eigenvalue, mo, i) )
mos.sort()
for i,orb2 in enumerate(mos):
if occ[orb2[2]] > 0.:
orb = orb2[1]
print "MO %4d MO 0.0 OCC NO = %9.7f ORB. ENERGY =%12.6f" \
%(i+1, occ[orb2[2]], orb.eigenvalue),
k=0
for c in orb.vector:
if k%5 == 0:
print "\n",
print "%15.8e"%(c),
k+=1
print "\n",
print "END DATA"
idx = res.mo_types.index(MO_type)
try:
print " THE HF ENERGY =%20.12f THE VIRIAL(-V/T)="%(res.energies[idx]),
print "%12.8f"%(res.virials[idx])
except TypeError:
print " THE HF ENERGY =%20.12f THE VIRIAL(-V/T)=%12.8f"%(0.,0.)
def write_wfnFile(file,out=sys.stdout):
if "basis" in file.mo_sets.keys() :
MoType = "basis"
else:
MoType = None
file.clean_uncontractions()
wfn_write(file,out,MoType)
fileTypes.append(wfnFile)
if __name__ == '__main__':
main(wfnFile)

732
resultsFile/Modules/xmvbFile.py Executable file
View File

@ -0,0 +1,732 @@
#!/usr/bin/python
# resultsFile is a library which allows to read output files of quantum
# chemistry codes and write input files.
# Copyright (C) 2007 Anthony SCEMAMA
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#
# Anthony Scemama
# LCPQ - IRSAMC
# Universite Paul Sabatier
# 118, route de Narbonne
# 31062 Toulouse Cedex 4
# scemama@irsamc.ups-tlse.fr
import include
eval(include.code)
import struct
import re
xmvbFile_defined_vars = [ "date", "version", "machine", "memory", "disk",\
"cpu_time", "author", "title", "units", "methods", "options", \
"spin_restrict", "conv_threshs", "energies", \
"ee_pot_energies", \
"Ne_pot_energies", "pot_energies", \
"kin_energies", "point_group", "num_elec", \
"charge", "multiplicity","nuclear_energy","dipole","geometry",\
"basis","mo_sets","mo_types","mulliken_mo","mulliken_ao",\
"mulliken_atom","lowdin_ao", "mulliken_atom","lowdin_atom",\
"two_e_int_ao", "determinants", "num_alpha", "num_beta",\
"closed_mos", "active_mos", "virtual_mos", \
"determinants_mo_type", "det_coefficients", \
"csf_mo_type", "csf_coefficients", "symmetries", "occ_num", \
"csf", "num_states"]
class xmvbFile(resultsFile):
""" Class defining the xmvb file.
"""
local_vars = list(local_vars)
defined_vars = list(xmvbFile_defined_vars)
def get_options(self):
if self._machine is None:
self.find_string("\\")
pos = self._pos
self.find_next_string("\\@")
end = self._pos
buffer = ""
self._options = []
for line in self.text[pos:end+1]:
buffer += line[1:].replace('\n','')
buffer = buffer.split('\\\\')
for l in buffer:
self._options.append(l.split('\\'))
self._options.pop()
return self._options
def get_machine(self):
if self._machine is None:
self._machine = self.options[0][2][5:].lower()
return self._machine
def get_version(self):
if self._version is None:
self._version = self.options[4][0].split('=')[1]
return self._version
def get_author(self):
if self._author is None:
self._author = self.options[0][7].lower()
return self._author
def get_charge(self):
if self._charge is None:
self._charge = float(self.options[3][0].split(',')[0])
return self._charge
def get_disk(self):
if self._disk is None:
try:
self.find_string("File lengths")
except IndexError:
return None
pos = self._pos
line = self.text[pos].split()
disk = 0
for i in line[4::2]:
disk += float(i)
disk = disk/1000.
if disk > 1.:
self._disk = str(disk)+" Gb"
else:
disk *= 1000.
if disk > 1.:
self._disk = str(disk)+" Mb"
else:
disk *= 1000.
self._disk = str(disk)+" kb"
return self._disk
def get_memory(self):
if self._memory is None:
try:
self.find_string("Leave Link")
except IndexError:
return None
pos = self._pos
line = self.text[pos].split()
memory = float(line[10])*8. / 1000000000.
if memory > 1.:
self._memory = str(memory)+" Gb"
else:
memory *= 1000.
if memory > 1.:
self._memory = str(memory)+" Mb"
else:
memory *= 1000.
self._memory = str(memory)+" kb"
return self._memory
def get_symmetries(self):
if self._symmetries is None:
try:
self.find_string("There are")
except IndexError:
return None
pos = self._pos
begin = pos
try:
self.find_next_string("Integral")
except IndexError:
return None
end = self._pos-1
sym = []
for k in range(begin,end):
buffer = self.text[k].split()
sym.append([buffer[8],int(buffer[2])])
self._symmetries = sym
return self._symmetries
def get_units(self):
if self._units is None:
try:
self.find_string("Coordinates")
except IndexError:
return None
pos = self._pos
units = self.text[pos].split()[4][1:-1]
if units != 'Angstroms':
self._units = 'BOHR'
else:
self._units = 'ANGS'
return self._units
def get_methods(self):
if self._methods is None:
methods = []
methods.append(self.options[0][4])
self._methods = methods
return self._methods
def get_spin_restrict(self):
if self._spin_restrict is None:
method = self.methods[0]
self._spin_restrict = True
if method == 'UHF': self._spin_restrict = False
return self._spin_restrict
def get_conv_threshs(self):
if self._conv_threshs is None:
self._conv_threshs = []
for m in self.methods:
if m == 'RHF' or m == 'UHF' or m == 'ROHF':
self.find_string("SCF Done")
pos = self._pos + 1
self._conv_threshs.append(float(self.text[pos].split()[2]))
if m == 'CASSCF':
self.find_string("Enter MCSCF program")
self.find_next_string("USED ACCURACY IN CHECKING CONVEGERGENCE")
pos = self._pos
self._conv_threshs.append(float(self.text[pos].split('=')[1]))
if self._conv_threshs == []:
self._conv_threshs = None
return self._conv_threshs
def get_ee_pot_energies(self):
if self._ee_pot_energies is None:
self._ee_pot_energies = []
for i,e in enumerate(self.kin_energies):
self._ee_pot_energies.append(self.energies[i]\
-self.nuclear_energy\
-self.kin_energies[i]\
-self.Ne_pot_energies[i])
return self._ee_pot_energies
def get_Ne_pot_energies(self):
if self._Ne_pot_energies is None:
self.find_string("N-N")
pos = self._pos
self._Ne_pot_energies = [float(self.text[pos].replace('=','= ').split()[3])]
return self._Ne_pot_energies
def get_point_group(self):
if self._point_group is None:
self._point_group = self.options[4][-1].split()[0].split('=')[1]
return self._point_group
def get_geometry(self):
if self._geometry is None:
self._geometry = []
self._pos = 0
pos=0
try:
while True:
pos = self._pos
self.find_next_string("Number Number Type")
self._pos += 1
except IndexError:
pass
pos +=1
self._pos=pos
self.find_next_string("-----")
end = self._pos
while pos<end:
temp = atom()
buffer = self.text[pos].split()
temp.charge = float(buffer[1])
temp.coord = (float(buffer[3]), float(buffer[4]), float(buffer[5]))
self._geometry.append(temp)
pos += 1
for i,line in enumerate(self.options[3][1:]):
buffer = line.split(',')
self._geometry[i].name = buffer[0]
# try:
# b = self.basis
# for f in b:
# for at in self._geometry:
# if f.center is at.coord:
# at.basis.append(f)
# except IndexError:
# pass
return self._geometry
def get_basis(self):
if self._basis is None:
gfprint=False
gfinput=False
buffer = self.options[1][0].split()
if 'GFPRINT' in buffer: gfprint=True
elif 'GFINPUT' in buffer: gfinput=True
if gfprint:
Polar=False
try:
self.find_string("Standard basis")
except:
self.find_string("General basis")
pos = self._pos
if "5D" in self.text[pos]: Polar=True
try:
self.find_next_string("AO basis set:")
pos = self._pos
except IndexError:
return None
self.find_string("Integral buffers")
end = self._pos
doLoop=True
while (doLoop):
try:
self.find_prev_string("There are")
end = self._pos
except:
doLoop = False
pos += 1
basis_read = []
line = self.text[pos].split()
iatom=0
atom = line[1]
while pos < end:
if line[0] == 'Atom':
index = int(line[3])
sym = line[4]
nfunc = int(line[5])
if atom != line[1]:
iatom += 1
atom = line[1]
bf = []
pos+=1
line = self.text[pos].split()
for k in range(nfunc):
expo = float(line[0].replace('D','E'))
coef = float(line[1].replace('D','E'))
if sym == "SP":
coef2 = float(line[2])
bf.append( [expo,coef,coef2] )
else:
bf.append( [expo,coef] )
pos += 1
line = self.text[pos].split()
if len(bf) > 0:
basis_read.append( [index,sym,bf,iatom] )
else:
print "GFPRINT should be present in the gaussian keywords."
return None
Nmax = basis_read[len(basis_read)-1][0]
basis = [None for i in range(Nmax)]
for b in basis_read:
basis[b[0]-1] = [b[1],b[2],b[3]]
NotNone = 0
ReadNone = False
for i in range(len(basis)-1,-1,-1):
if basis[i] == None:
ReadNone = True
basis[i] = list(basis[i+NotNone])
else:
if ReadNone:
NotNone = 0
ReadNone = False
NotNone += 1
k=0
while k<len(basis):
if basis[k][0] == "S":
mylist = []
elif basis[k][0] == "P":
mylist = [ "X", "Y", "Z" ]
elif basis[k][0] == "SP":
mylist = [ "S", "X", "Y", "Z" ]
elif basis[k][0] == "D":
if not Polar:
mylist = [ "XX", "YY", "ZZ", "XY", "XZ", "YZ" ]
else:
mylist = [ "D0", "D1+", "D1-", "D2+", "D2-" ]
elif basis[k][0] == "F":
if not Polar:
mylist = [ "XXX", "YYY", "ZZZ", "YYX", "XXY", "XXZ", "ZZX", "ZZY",
"YYZ", "XYZ" ]
else:
mylist = [ "F0", "F1+", "F1-", "F2+", "F2-", "F3+", "F3-" ]
elif basis[k][0] == "G":
if not Polar:
mylist = [ "ZZZZ", "ZZZY", "YYZZ", "YYYZ", "YYYY", "ZZZX",
"ZZXY", "YYXZ","YYYX", "XXZZ", "XXYZ",
"XXYY", "XXXZ","XXXY", "XXXX" ]
else:
mylist = [ "G0", "G1+", "G1-", "G2+", "G2+", "G3+", "G3-",
"G4+", "G4-" ]
#elif basis[k][0] == "H":
# mylist = [ "XXXXX", "YYYYY", "ZZZZZ", "XXXXY", "XXXXZ", "YYYYX", \
# "YYYYZ", "ZZZZX", "ZZZZY", "XXXYY", "XXXZZ", \
# "YYYXX", "YYYZZ", "ZZZXX", "ZZZYY", "XXXYZ", \
# "YYYXZ", "ZZZXY", "XXYYZ", "XXZZY", "YYZZX" ]
#elif basis[k][0] == "I":
# mylist = [ "XXXXXX", "YYYYYY", "ZZZZZZ", "XXXXXY", "XXXXXZ", \
# "YYYYYX", "YYYYYZ", "ZZZZZX", "ZZZZZY", "XXXXYY", \
# "XXXXZZ", "YYYYXX", "YYYYZZ", "ZZZZXX", "ZZZZYY", \
# "XXXXYZ", "YYYYXZ", "ZZZZXY", "XXXYYY", "XXXZZZ", \
# "YYYZZZ", "XXXYYZ", "XXXZZY", "YYYXXZ", "YYYZZX", \
# "ZZZXXY", "ZZZYYX" ]
mylist = map(normalize_basis_name,mylist)
for i in mylist[:-1]:
basis[k][0] = i
basis.insert(k,list(basis[k]))
k+=1
for i in mylist[-1:]:
basis[k][0] = i
k+=1
self._basis = []
for buffer in basis:
contr = contraction()
for c in buffer[1]:
gauss = gaussian()
atom = self.geometry[int(buffer[2])]
gauss.center = atom.coord
gauss.expo = c[0]
gauss.sym = buffer[0]
if len(c) == 3:
if gauss.sym == "S":
contr.append(c[1],gauss)
else:
contr.append(c[2],gauss)
else:
contr.append(c[1],gauss)
self._basis.append(contr)
return self._basis
def get_mo_types(self):
if self._mo_types is None:
self.get_mo_sets()
return self._mo_types
def get_mulliken_mo(self):
if self._mulliken_mo is None:
pass
return self._mulliken_mo
def get_mulliken_ao(self):
if self._mulliken_ao is None:
pass
return self._mulliken_ao
def get_lowdin_ao(self):
if self._lowdin_ao is None:
pass
return self._lowdin_ao
def get_mulliken_atom(self):
if self._mulliken_atom is None:
try:
self.find_string("Mulliken atomic charges")
except IndexError:
return None
self._pos += 2
pos = self._pos
self.find_next_string("Sum of Mulliken")
end = self._pos
line = self.text[pos].split()
vm = orbital()
vm.set = "Mulliken_Atom"
vm.eigenvalue = 0.
while pos < end:
value = float(line[2])
vm.vector.append(value)
vm.eigenvalue += value
pos += 1
line = self.text[pos].split()
self._mulliken_atom = vm
return self._mulliken_atom
def get_lowdin_atom(self):
if self._lowdin_atom is None:
pass
return self._lowdin_atom
def get_two_e_int_ao(self):
if self._two_e_int_ao is None:
pass
return self._two_e_int_ao
def get_closed_mos(self):
if self._closed_mos is None:
result = []
maxmo = len(self.mo_sets[self.determinants_mo_type])
for orb in range(maxmo):
present = True
for det in self.determinants:
for spin_det in det:
if orb not in det[spin_det]: present = False
if present: result.append(orb)
self._closed_mos = result
return self._closed_mos
def get_virtual_mos(self):
if self._virtual_mos is None:
result = []
minmo = len(self.closed_mos)
maxmo = len(self.mo_sets[self.determinants_mo_type])
for orb in range(minmo,maxmo):
present = False
for det in self.determinants:
for spin_det in det:
if orb in det[spin_det]: present = True
if not present: result.append(orb)
self._virtual_mos = result
return self._virtual_mos
def get_active_mos(self):
if self._active_mos is None:
cl = self.closed_mos
vi = self.virtual_mos
maxmo = len(self.mo_sets[self.determinants_mo_type])
result = []
for i in range(maxmo):
present = i in cl or i in vi
if not present:
result.append(i)
self._active_mos = result
return self._active_mos
def get_occ_num(self):
if self._occ_num is None:
if self.mulliken_mo is not None:
occ = {}
for motype in self.mo_types:
occ[motype] = [ mo.eigenvalue for mo in self.mulliken_mo ]
if occ != {}:
self._occ_num = occ
return self._occ_num
def get_date(self):
if self._date is None:
self.find_string("Job started at")
pos = self._pos
self._date = self.text[pos][16:]
return self._date
def get_multiplicity(self):
if self._multiplicity is None:
self.find_string("nmul=")
pos = self._pos
self._multiplicity = int(self.text[pos].split('=')[-1])
return self._multiplicity
def get_num_elec(self):
if self._num_elec is None:
self.find_string("nelectron=")
pos = self._pos
self._num_elec= int(self.text[pos].split('=')[1].split()[0])
return self._num_elec
def get_nuclear_energy(self):
if self._nuclear_energy is None:
self.find_string("Nuclear Repulsion Energy:")
pos = self._pos
self._nuclear_energy = float(self.text[pos].split(':')[1])
return self._nuclear_energy
def get_title(self):
if self._title is None:
self.find_string("End of Input")
pos = self._pos+3
self._title = self.text[pos].strip()
return self._title
def get_cpu_time(self):
if self._cpu_time is None:
try:
self.find_string("Cpu for the Job:")
except IndexError:
return None
pos = self._pos
self._cpu_time = self.text[pos].split(':')[1].split()[0]+' s'
return self._cpu_time
def get_dipole(self):
if self._dipole is None:
self.find_string("Dipole moment")
pos = self._pos+2
line = self.text[pos].split()
self._dipole = []
self._dipole.append(float(line[1]))
self._dipole.append(float(line[3]))
self._dipole.append(float(line[5]))
return self._dipole
def get_energies(self):
if self._energies is None:
self._energies = []
self.find_string("Total Energy:")
pos = self._pos
self._energies.append(float(self.text[pos].split(':')[1]))
return self._energies
def get_pot_energies(self):
if self._pot_energies is None:
self._pot_energies = []
self.find_string("Potential energy:")
pos = self._pos
self._pot_energies.append(float(self.text[pos].split(':')[1]))
return self._pot_energies
def get_kin_energies(self):
if self._kin_energies is None:
self._kin_energies = []
self.find_string("Kinetic energy:")
pos = self._pos
self._kin_energies.append(float(self.text[pos].split(':')[1]))
return self._kin_energies
def get_mo_sets(self):
if self._mo_sets is None:
self._mo_types = ['VB']
posend = {}
self.find_string("norb=")
norb = int(self.text[self._pos].split("norb=")[1].split()[0])
self.find_string("OPTIMIZED ORBITALS")
pos = self._pos+3
iorb=0
vectors = []
while(iorb<norb):
lorb = len(self.text[pos].split())
pos += 1
begin = pos
for l in range(lorb):
iorb+=1
pos = begin
line = self.text[pos].split()
v = orbital()
v.set = self.mo_types[0]
v.basis = None
v.eigenvalue = float(iorb)
while len(line) > 0:
v.vector.append(float(line[l+1]))
pos += 1
line = self.text[pos].split()
vectors.append(v)
self._mo_sets = {}
self._mo_sets['VB'] = vectors
return self._mo_sets
def get_num_alpha(self):
if self._num_alpha is None:
self._num_alpha = self.num_elec/2 + (self.multiplicity-1)/2
return self._num_alpha
def get_num_beta(self):
if self._num_beta is None:
self._num_beta = self.num_elec/2 - (self.multiplicity-1)/2
return self._num_beta
def get_determinants_mo_type(self):
if self._determinants_mo_type is None:
self._determinants_mo_type = self.mo_types[-1]
return self._determinants_mo_type
def get_csf_mo_type(self):
if self._csf_mo_type is None:
self._csf_mo_type = self.determinants_mo_type
return self._csf_mo_type
def get_determinants(self):
if self._determinants is None:
determinants = []
if self.csf is not None:
for csf in self.csf:
for new_det in csf.determinants:
determinants.append(new_det)
else:
pass
if determinants != []:
self._determinants_mo_type = self.mo_types[-1]
self._determinants = determinants
return self._determinants
def get_csf(self):
if self._csf is None:
csf = []
csf_coefficients = []
self.find_string('COEFFICIENTS OF DETERMINANTS')
self.find_next_string('1')
pos = self._pos
buffer = self.text[pos]
while buffer.strip() != "":
tempcsf_a = []
tempcsf_b = []
buffer = self.text[pos]
coef = float(buffer.split()[1])
ii=0
while ii < self.num_elec:
buffer = buffer[24:].split()
for i in buffer:
ii+=1
if ii <= self.num_alpha:
tempcsf_a.append(int(i)-1)
else:
tempcsf_b.append(int(i)-1)
pos += 1
buffer = self.text[pos]
this_csf = CSF()
this_csf.append(1.,tempcsf_a,tempcsf_b)
csf.append(this_csf)
csf_coefficients.append([coef])
if csf != []:
self._csf = csf
self._csf_coefficients = csf_coefficients
return self._csf
def get_det_coefficients(self):
if self._det_coefficients is None:
if self.csf is not None:
self._det_coefficients = []
csf = self.csf
vector = []
for state_coef in self.csf_coefficients:
for i,c in enumerate(state_coef):
for d in csf[i].coefficients:
vector.append(c*d)
self._det_coefficients.append(vector)
return self._det_coefficients
def get_csf_coefficients(self):
if self._csf_coefficients is None:
self.get_csf()
return self._csf_coefficients
def get_num_states(self):
if self._num_states is None:
self._num_states=1
return self._num_states
to_remove = []
for i, j in local_vars:
if i in resultsFile_defined_vars:
to_remove.append( (i,j) )
for i in to_remove:
local_vars.remove(i)
for i, j in local_vars:
if i not in defined_vars:
exec build_get_funcs(i) in locals()
exec build_property(i,j) in locals()
del to_remove, i, j
fileTypes.insert(0,xmvbFile)
if __name__ == '__main__':
main(xmvbFile)
###### END #####

44
resultsFile/__init__.py Executable file
View File

@ -0,0 +1,44 @@
#!/usr/bin/env python
# resultsFile is a library which allows to read output files of quantum
# chemistry codes and write input files.
# Copyright (C) 2007 Anthony SCEMAMA
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#
# Anthony Scemama
# LCPQ - IRSAMC
# Universite Paul Sabatier
# 118, route de Narbonne
# 31062 Toulouse Cedex 4
# scemama@irsamc.ups-tlse.fr
"""resultsFile library."""
__author__ = "Anthony SCEMAMA <scemama@irsamc.ups-tlse.fr>"
__date__ = "20 Nov 2007"
import os
all = [ "resultsFile", "getFile", "lib", "Modules" ]
for mod in all:
try:
exec 'from '+mod+' import *'
except:
print "Error importing module", mod
pass

51
resultsFile/cython_setup Executable file
View File

@ -0,0 +1,51 @@
#!/usr/bin/python
# resultsFile is a library which allows to read output files of quantum
# chemistry codes and write input files.
# Copyright (C) 2007 Anthony SCEMAMA
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#
# Anthony Scemama
# LCPQ - IRSAMC
# Universite Paul Sabatier
# 118, route de Narbonne
# 31062 Toulouse Cedex 4
# scemama@irsamc.ups-tlse.fr
from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext
"""
./cython_setup build_ext --inplace
"""
ext_modules = [Extension("resultsFile_cython", ["resultsFile_cython.pyx"])]
import os
setup(name="resultsFile",
version=os.getenv("VERSION","1.0"),
author="Anthony Scemama",
author_email="scemama@irsamc.ups-tlse.fr",
license="gpl-license",
description="Module for I/O on Quantum Chemistry files.",
packages=["resultsFile","resultsFile.lib","resultsFile.Modules"],
cmdclass = {'build_ext': build_ext},
ext_modules = ext_modules
)

54
resultsFile/getFile.py Executable file
View File

@ -0,0 +1,54 @@
#!/usr/bin/python
# resultsFile is a library which allows to read output files of quantum
# chemistry codes and write input files.
# Copyright (C) 2007 Anthony SCEMAMA
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#
# Anthony Scemama
# LCPQ - IRSAMC
# Universite Paul Sabatier
# 118, route de Narbonne
# 31062 Toulouse Cedex 4
# scemama@irsamc.ups-tlse.fr
from resultsFile import *
import sys
# Find fileType
def getFile(name):
MyList = fileTypes+[None]
for fileType in MyList:
try:
file = fileType(name)
n = file.num_elec
if n is None or n <=0 or not isinstance(n,int):
raise TypeError
break
except:
pass
if fileType is None:
print "File type not recognized."
sys.exit(1)
return file
if __name__=='__main__':
getFile(sys.argv[1])

416
resultsFile/resultsFile.py Executable file
View File

@ -0,0 +1,416 @@
#!/usr/bin/python
# resultsFile is a library which allows to read output files of quantum
# chemistry codes and write input files.
# Copyright (C) 2007 Anthony SCEMAMA
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#
# Anthony Scemama
# LCPQ - IRSAMC
# Universite Paul Sabatier
# 118, route de Narbonne
# 31062 Toulouse Cedex 4
# scemama@irsamc.ups-tlse.fr
from lib import *
import lib.basis as Basis
import sys
import copy
fileTypes = []
local_vars = [ \
# File properties
( 'filename' , "Name of the results file."),
( 'text' , "Text of the results file."),
( 'pos' , "Position in the results file."),
( 'date' , "When the calculation was performed."),
( 'version' , "Version of the code generating the file."),
( 'author' , "Who ran the calculation."),
( 'machine' , "Machine where the calculation was run."),
( 'memory' , "Requested memory for the calculation."),
( 'disk' , "Requested disk space for the calculation."),
( 'num_proc' , "Number of processors used."),
( 'cpu_time' , "CPU time."),
# General properties
( 'title' , "Title of the run."),
( 'options' , "Options given in the input file."),
( 'units' , "Units for the geometry (au or angstroms)."),
( 'methods' , "List of calculation methods."),
( 'spin_restrict' , "Open-shell or closed-shell calculations."),
( 'conv_threshs' , "List of convergence thresholds."),
( 'energies' , "List of energies."),
( 'one_e_energies', "List of one electron energies."),
( 'two_e_energies', "List of two electron energies."),
( 'ee_pot_energies',"List of electron-electron potential energies."),
( 'Ne_pot_energies',"List of nucleus-electron potential energies."),
( 'pot_energies' , "List of potential energies."),
( 'kin_energies' , "List of kinetic energies."),
( 'virials' , "Virial ratios."),
( 'mulliken_mo' , "Mulliken atomic population in each MO."),
( 'mulliken_ao' , "Mulliken atomic population in each AO."),
( 'lowdin_ao' , "Lowdin atomic population in each AO."),
( 'mulliken_atom' , "Mulliken atomic population."),
( 'lowdin_atom' , "Lowdin atomic population."),
( 'dipole' , "Dipole moment"),
( 'quadrupole' , "Quadrupole moment"),
( 'num_states' , "Number of electronic states"),
# Geometry properties
( 'point_group' , "Symmetry used."),
( 'geometry' , "Atom types and coordinates."),
( 'symmetries' , "Irreducible representations"),
( 'num_elec' , "Number of electrons."),
( 'num_alpha' , "Number of Alpha electrons."),
( 'num_beta' , "Number of Beta electrons."),
( 'charge' , "Charge of the system."),
( 'multiplicity' , "Spin multiplicity of the system."),
( 'nuclear_energy', "Repulsion of the nuclei."),
( 'gradient_energy', "Gradient of the Energy wrt nucl coord."),
# Basis set
( 'basis' , "Basis set definition"),
( 'uncontracted_basis', "Uncontracted Basis set"),
# Orbitals
( 'mo_sets' , "List of molecular orbitals"),
( 'mo_types' , "Types of molecular orbitals (canonical, natural,...)"),
( 'occ_num' , "Occupation numbers"),
( 'uncontracted_mo_sets', "List of molecular orbitals in the uncontracted basis set."),
# Determinants
( 'closed_mos' , "Closed shell molecular orbitals"),
( 'active_mos' , "Active molecular orbitals"),
( 'virtual_mos' , "Virtual molecular orbitals"),
( 'csf' , "List of Configuration State Functions"),
( 'determinants' , "List of Determinants"),
( 'csf_mo_type' , "MO type of the determinants"),
( 'determinants_mo_type' , "MO type of the determinants"),
( 'csf_coefficients', "Coefficients of the CSFs"),
( 'det_coefficients', "Coefficients of the determinants"),
# Integrals
( 'one_e_int_ao' , "One electron integrals in AO basis"),
( 'two_e_int_ao' , "Two electron integrals in AO basis"),
( 'one_e_int_mo' , "One electron integrals in MO basis"),
( 'two_e_int_mo' , "Two electron integrals in MO basis"),
]
resultsFile_defined_vars = ["text","uncontracted_basis", "uncontracted_mo_sets"]
class resultsFile(object):
""" Class containing the definition of files.
"""
local_vars = list(local_vars)
defined_vars = list(resultsFile_defined_vars)
def __init__(self,name):
"""All local variables are set to None, except filename.
"""
for var, doc in local_vars:
setattr(self,'_'+var,None)
self._filename = name
def get_text(self):
"""Reads the whole file.
"""
if self._text is None:
try:
file = open(self.filename,"r")
self._text = file.readlines()
except IOError:
print "Unable to open "+self.filename
sys.exit(1)
file.close()
return self._text
def get_uncontracted_basis(self):
if self._uncontracted_basis is None:
try:
from resultsFile_cython import get_uncontracted_basis as f
has_cython = True
except ImportError:
has_cython = False
basis = self.basis
if has_cython:
self._uncontracted_basis = f(basis)
else:
uncontr = []
for contr in basis:
for b in contr.prim:
uncontr.append(b)
self._uncontracted_basis = uncontr
return self._uncontracted_basis
def get_uncontracted_mo_sets(self):
if self._uncontracted_mo_sets is None:
try:
from resultsFile_cython import get_uncontracted_mo_sets as f
has_cython = True
except ImportError:
has_cython = False
if has_cython:
self._uncontracted_mo_sets = f(self.basis,
self.uncontracted_basis,
self.mo_sets,
self.mo_types)
else:
uncontr = {}
basis = self.basis
for motype in self.mo_types:
uncontr[motype] = []
for mo in self.mo_sets[motype]:
lenmovector = len(mo.vector)
monew = orbital()
monew.basis = self.uncontracted_basis
monew.eigenvalue = mo.eigenvalue
monew.set = motype
v = []
for i, contr in enumerate(basis):
if i<lenmovector:
ci = mo.vector[i]
if ci == 0.:
for p in range(len(contr.prim)):
v.append(0.)
else:
for p, c in zip(contr.prim,contr.coef):
v.append(c*ci/p.norm)
monew.vector = v
uncontr[motype].append(monew)
self._uncontracted_mo_sets = uncontr
#endif
return self._uncontracted_mo_sets
def clean_contractions(self):
basis = self.basis
newbasis = []
idx = range(len(basis))
for k,b1 in enumerate(basis):
addBasis=True
for l, b2 in enumerate(basis[:k]):
if b2 == b1:
idx[k] = l
addBasis=False
break
if addBasis:
newbasis.append(b1)
self._basis = newbasis
mo_sets = self.mo_sets
for motype in self.mo_types:
for mo in mo_sets[motype]:
lenmovector = len(mo.vector)
newvec = [None for i in idx]
for i in idx:
newvec[i] = 0.
for k,l in enumerate(idx):
if k < lenmovector:
newvec[l] += mo.vector[k]
mo.vector = []
for c in newvec:
if c is not None:
mo.vector.append(c)
def clean_uncontractions(self):
basis = self.uncontracted_basis
newbasis = []
idx = range(len(basis))
for k,b1 in enumerate(basis):
addBasis=True
for l, b2 in enumerate(basis[:k]):
if b2 == b1:
idx[k] = l
addBasis=False
break
if addBasis:
newbasis.append(b1)
self._uncontracted_basis = newbasis
mo_sets = self.uncontracted_mo_sets
for motype in self.mo_types:
for mo in mo_sets[motype]:
lenmovector = len(mo.vector)
newvec = [None for i in idx]
for i in idx:
newvec[i] = 0.
for k,l in enumerate(idx):
if k < lenmovector:
newvec[l] += mo.vector[k]
mo.vector = []
for c in newvec:
if c is not None:
mo.vector.append(c)
def convert_to_cartesian(self):
basis = self.basis
newbasis = []
idx = range(len(basis))
map = []
weight = []
for i,b in enumerate(basis):
l, m = Basis.get_lm(b.sym)
if l is None:
newbasis.append(b)
map.append(i)
weight.append(1.)
else:
powers, coefs = xyz_from_lm(l,m)
for j,prim in enumerate(b.prim):
b.coef[j] /= prim.norm
for c, p in zip(coefs, powers):
contr = copy.deepcopy(b)
sym = ''
for l,letter in enumerate('xyz'):
sym += p[l]*letter
contr.sym = sym
for j,prim in enumerate(contr.prim):
prim.sym = sym
contr.coef[j] *= prim.norm
newbasis.append(contr)
map.append(i)
weight.append(c)
mo_sets = self.mo_sets
for motype in self.mo_types:
for mo in mo_sets[motype]:
newvec = []
vec = mo.vector
for i,w in zip(map,weight):
newvec.append(vec[i]*w)
mo.vector = newvec
same_as = {}
for i,b1 in enumerate(newbasis):
for j,b2 in enumerate(newbasis[:i]):
if b1 == b2:
same_as[i] = j
weight[j] += weight[i]
break
to_remove = same_as.keys()
to_remove.sort()
to_remove.reverse()
for i in to_remove:
newbasis.pop(i)
weight.pop(i)
map.pop(i)
for motype in self.mo_types:
for mo in mo_sets[motype]:
for i in to_remove:
index = same_as[i]
value = mo.vector.pop(i)
mo.vector[index] += value
self._basis = newbasis
self._mo_sets = mo_sets
def find_string(self,chars):
"""Finds the 1st occurence of chars.
"""
self._pos = 0
self.find_next_string(chars)
def find_last_string(self,chars):
"""Finds the 1st occurence of chars.
"""
self._pos = len(self.text)-1
self.find_prev_string(chars)
def find_next_string(self,chars):
"""Finds the next occurence of chars.
"""
pos = self._pos
text = self.text
found = False
while not found and pos < len(text):
if chars in text[pos]:
found = True
else:
pos += 1
if not found:
raise IndexError
self._pos = pos
def find_prev_string(self,chars):
"""Finds the next occurence of chars.
"""
pos = self._pos
text = self.text
found = False
while not found and pos < len(text):
if chars in text[pos]:
found = True
else:
pos -= 1
if not found:
raise IndexError
self._pos = pos
for i, j in local_vars:
if i not in defined_vars:
exec build_get_funcs(i) in locals()
exec build_property(i,j) in locals()
del i,j
def main(fileType):
import getopt
print """
resultsFile version 1.0, Copyright (C) 2007 Anthony SCEMAMA
resultsFile comes with ABSOLUTELY NO WARRANTY; for details see the
gpl-license file.
This is free software, and you are welcome to redistribute it
under certain conditions; for details see the gpl-license file."""
full_list = fileType.defined_vars + resultsFile.defined_vars
try:
opts, args = getopt.gnu_getopt(sys.argv[1:],'',full_list)
except getopt.GetoptError:
args = []
if len(args) == 0:
usage(fileType)
sys.exit(2)
f = fileType(args[0])
for o,a in opts:
print o[2:]
print ''.join(['-' for k in o[2:]])
PrettyPrint = prettyPrint
exec 'PrettyPrint(f.'+o[2:]+')' in locals()
print ""
sys.exit(0)
def usage(fileType):
print ""
print "Usage:"
print "------"
print ""
print sys.argv[0], '[options] file'
print ""
print "Options:"
print "--------"
print ""
for o in fileType.defined_vars + resultsFile.defined_vars:
line = (" --"+o).ljust(30)+': '
for l in fileType.local_vars:
if l[0] == o:
line += l[1]
break
print line
print ""
if __name__ == '__main__':
main(resultsFile)

241
resultsFile/resultsFile_cython.pyx Executable file
View File

@ -0,0 +1,241 @@
#!/usr/bin/python
# resultsFile is a library which allows to read output files of quantum
# chemistry codes and write input files.
# Copyright (C) 2007 Anthony SCEMAMA
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#
# Anthony Scemama
# LCPQ - IRSAMC
# Universite Paul Sabatier
# 118, route de Narbonne
# 31062 Toulouse Cedex 4
# scemama@irsamc.ups-tlse.fr
from lib import *
import lib.basis as Basis
def get_uncontracted_basis(basis):
uncontr = []
for contr in basis:
for b in contr.prim:
uncontr.append(b)
return uncontr
def get_uncontracted_mo_sets(basis,uncontracted_basis,mo_sets,mo_types):
cdef dict uncontr = {}
cdef int lenmovector
cdef int lenbasis = len(basis)
cdef double ci
cdef int i, j
for motype in mo_types:
uncontr[motype] = []
for mo in mo_sets[motype]:
lenmovector = len(mo.vector)
monew = orbital()
monew.basis = uncontracted_basis
monew.eigenvalue = mo.eigenvalue
monew.set = motype
v = []
for i in range(lenbasis):
contr = basis[i]
if i<lenmovector:
ci = mo.vector[i]
if ci == 0.:
for j in range(len(contr.prim)):
v.append(0.)
else:
for p, c in zip(contr.prim,contr.coef):
v.append(c*ci/p.norm)
monew.vector = v
uncontr[motype].append(monew)
return uncontr
def clean_contractions(basis):
newbasis = []
cdef int i, k, l, lenmovector
idx = range(len(basis))
for k,b1 in enumerate(basis):
addBasis=True
for l, b2 in enumerate(basis[:k]):
if b2 == b1:
idx[k] = l
addBasis=False
break
if addBasis:
newbasis.append(b1)
self._basis = newbasis
mo_sets = self.mo_sets
for motype in self.mo_types:
for mo in mo_sets[motype]:
lenmovector = len(mo.vector)
newvec = [None for i in idx]
for i in idx:
newvec[i] = 0.
for k,l in enumerate(idx):
if k < lenmovector:
newvec[l] += mo.vector[k]
mo.vector = []
for c in newvec:
if c is not None:
mo.vector.append(c)
# def clean_uncontractions(self):
# basis = self.uncontracted_basis
# newbasis = []
# idx = range(len(basis))
# for k,b1 in enumerate(basis):
# addBasis=True
# for l, b2 in enumerate(basis[:k]):
# if b2 == b1:
# idx[k] = l
# addBasis=False
# break
# if addBasis:
# newbasis.append(b1)
# self._uncontracted_basis = newbasis
# mo_sets = self.uncontracted_mo_sets
# for motype in self.mo_types:
# for mo in mo_sets[motype]:
# lenmovector = len(mo.vector)
# newvec = [None for i in idx]
# for i in idx:
# newvec[i] = 0.
# for k,l in enumerate(idx):
# if k < lenmovector:
# newvec[l] += mo.vector[k]
# mo.vector = []
# for c in newvec:
# if c is not None:
# mo.vector.append(c)
#
# def convert_to_cartesian(self):
# basis = self.basis
# newbasis = []
# idx = range(len(basis))
# map = []
# weight = []
# for i,b in enumerate(basis):
# l, m = Basis.get_lm(b.sym)
# if l is None:
# newbasis.append(b)
# map.append(i)
# weight.append(1.)
# else:
# powers, coefs = xyz_from_lm(l,m)
# for j,prim in enumerate(b.prim):
# b.coef[j] /= prim.norm
# for c, p in zip(coefs, powers):
# contr = copy.deepcopy(b)
# sym = ''
# for l,letter in enumerate('xyz'):
# sym += p[l]*letter
# contr.sym = sym
# for j,prim in enumerate(contr.prim):
# prim.sym = sym
# contr.coef[j] *= prim.norm
# newbasis.append(contr)
# map.append(i)
# weight.append(c)
# mo_sets = self.mo_sets
# for motype in self.mo_types:
# for mo in mo_sets[motype]:
# newvec = []
# vec = mo.vector
# for i,w in zip(map,weight):
# newvec.append(vec[i]*w)
# mo.vector = newvec
# same_as = {}
# for i,b1 in enumerate(newbasis):
# for j,b2 in enumerate(newbasis[:i]):
# if b1 == b2:
# same_as[i] = j
# weight[j] += weight[i]
# break
# to_remove = same_as.keys()
# to_remove.sort()
# to_remove.reverse()
# for i in to_remove:
# newbasis.pop(i)
# weight.pop(i)
# map.pop(i)
# for motype in self.mo_types:
# for mo in mo_sets[motype]:
# for i in to_remove:
# index = same_as[i]
# value = mo.vector.pop(i)
# mo.vector[index] += value
# self._basis = newbasis
# self._mo_sets = mo_sets
# def find_string(self,chars):
# """Finds the 1st occurence of chars.
# """
# self._pos = 0
# self.find_next_string(chars)
# def find_last_string(self,chars):
# """Finds the 1st occurence of chars.
# """
# self._pos = len(self.text)-1
# self.find_prev_string(chars)
# def find_next_string(self,chars):
# """Finds the next occurence of chars.
# """
# pos = self._pos
# text = self.text
# found = False
# while not found and pos < len(text):
# if chars in text[pos]:
# found = True
# else:
# pos += 1
# if not found:
# raise IndexError
# self._pos = pos
# def find_prev_string(self,chars):
# """Finds the next occurence of chars.
# """
# pos = self._pos
# text = self.text
# found = False
# while not found and pos < len(text):
# if chars in text[pos]:
# found = True
# else:
# pos -= 1
# if not found:
# raise IndexError
# self._pos = pos
# for i, j in local_vars:
# if i not in defined_vars:
# exec build_get_funcs(i) in locals()
# exec build_property(i,j) in locals()
# del i,j

54
setup.py Executable file
View File

@ -0,0 +1,54 @@
#!/usr/bin/python
# resultsFile is a library which allows to read output files of quantum
# chemistry codes and write input files.
# Copyright (C) 2007 Anthony SCEMAMA
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#
# Anthony Scemama
# LCPQ - IRSAMC
# Universite Paul Sabatier
# 118, route de Narbonne
# 31062 Toulouse Cedex 4
# scemama@irsamc.ups-tlse.fr
from distutils.core import setup
"""
- To create a source distribution, run:
python setup.py sdist
- To install the module, run:
python setup.py install
- To create an rpm distribution
python setup.py bdist_rpm
"""
import os
setup(name="resultsFile",
version=os.getenv("VERSION","1.0"),
author="Anthony Scemama",
author_email="scemama@irsamc.ups-tlse.fr",
license="gpl-license",
description="Module for I/O on Quantum Chemistry files.",
packages=["resultsFile","resultsFile.lib","resultsFile.Modules"]
)