LCPQ/resultsFile
10
0
Fork 0
mirror of https://gitlab.com/scemama/resultsFile.git synced 2024-12-22 04:13:41 +01:00
Code Releases Activity

Merge pull request #2 from scemama/master

Browse Source 
Forgot lib
This commit is contained in:
Anthony Scemama 2015-01-19 17:06:39 +01:00
commit b467f187b0
11 changed files with 1421 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
.gitignore vendored
View File

@ -13,14 +13,6 @@ develop-eggs/
dist/
downloads/
eggs/
lib/
lib64/
parts/
sdist/
var/
*.egg-info/
.installed.cfg
*.egg
# PyInstaller
# Usually these files are written by a python script from a template

46
resultsFile/lib/__init__.py Executable file
View File

@ -0,0 +1,46 @@
#!/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
wd = os.path.dirname(__file__)
all = [ i[:-3] for i in os.listdir(wd) if i.endswith(".py") ]
for mod in all:
try:
exec 'from '+mod+' import *'
except:
print "Error importing module", mod
pass

217
resultsFile/lib/atom.py Executable file
View File

@ -0,0 +1,217 @@
#!/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
from library import *
from math import *
class atom(object):
"""Class for an atom."""
def __init__(self):
self._name = None
self._charge = None
self._coord = None
self._basis = None
def __repr__(self):
out = "%8s %10.3f %10.6f %10.6f %10.6f"%tuple(\
[self.name, self.charge]+list(self.coord))
return out
def __repr__debug__(self):
out = ""
out += "Atom:\n"
out += " Name : "+str(self.name)+'\n'
out += " Charge : "+str(self.charge)+'\n'
out += " Coord : "+str(self.coord)+'\n'
out += " Basis set : "+str(self.basis)+'\n'
return out
def __cmp__(self,other):
assert ( isinstance(other,atom) )
if self.charge < other.charge:
return -1
elif self.charge > other.charge:
return 1
elif self.charge == other.charge:
return 0
for i in "name charge coord basis".split():
exec """
def get_%(i)s(self): return self._%(i)s
def set_%(i)s(self,value): self._%(i)s = value
%(i)s = property(fget=get_%(i)s,fset=set_%(i)s) """%locals()
class atomDataElement(object):
"""Atomic database."""
def __init__(self,symbol,charge,covalR,vdwR,valence,color,mass):
self.symbol = symbol
self.charge = int(charge)
self.covalR = float(covalR)
self.vdwR = float(vdwR)
self.mass = float(mass)
self.valence= int(valence)
self.color = int(color)
def __repr__(self):
out = ""
out += "%3d "%(self.charge)
out += "%3s "%(self.symbol)
out += "%10.6f "%(self.covalR)
out += "%10.6f "%(self.vdwR)
out += "%3d "%(self.valence)
out += "%3d "%(self.color)
out += "%10.6f "%(self.mass)
return out
atomDataText = """
0 X 0.0 0.75 0 8 0.0
0 XX 0.0 0.75 0 8 0.0
1 H 0.370000 1.200000 1 0 1.0079
2 He 0.700000 1.700000 0 2 4.00260
3 Li 1.230000 1.700000 1 5 6.941
4 Be 0.890000 1.700000 2 5 9.01218
5 B 0.900000 1.700000 3 5 10.81
6 C 0.850000 1.700000 4 8 12.011
7 N 0.740000 1.550000 3 3 14.0067
8 O 0.740000 1.520000 2 1 15.9994
9 F 0.720000 1.470000 1 7 18.998403
10 Ne 0.700000 1.700000 0 1 20.179
11 Na 1.000000 1.700000 1 6 22.98977
12 Mg 1.360000 1.700000 2 6 24.305
13 Al 1.250000 1.940000 3 6 26.98154
14 Si 1.170000 2.100000 4 5 28.0855
15 P 1.100000 1.800000 3 6 30.97376
16 S 1.100000 1.800000 6 4 32.06
17 Cl 0.990000 1.750000 1 6 35.453
18 Ar 0.700000 1.700000 0 5 39.948
19 K 2.030000 1.700000 1 5 39.0983
20 Ca 1.740000 1.700000 2 5 40.08
21 Sc 1.440000 1.700000 0 5 44.9559
22 Ti 1.320000 1.700000 0 5 47.90
23 V 1.220000 1.980000 0 5 50.9415
24 Cr 0.000000 1.940000 0 5 51.996
25 Mn 1.160000 1.930000 0 5 54.9380
26 Fe 0.000000 1.930000 0 5 55.9332
27 Co 1.150000 1.920000 0 5 58.9332
28 Ni 1.170000 1.700000 0 5 58.70
29 Cu 1.250000 1.700000 0 5 63.546
30 Zn 1.250000 1.700000 0 5 65.38
31 Ga 1.200000 2.020000 0 0 69.72
32 Ge 1.210000 1.700000 0 5 72.59
33 As 1.160000 1.960000 0 1 74.9216
34 Se 0.700000 1.700000 0 5 78.96
35 Br 1.240000 2.100000 0 6 79.904
36 Kr 1.910000 1.700000 0 5 83.80
37 Rb 1.620000 1.700000 0 5 85.4678
38 Sr 1.450000 1.700000 0 5 87.62
39 Y 1.340000 1.700000 0 5 88.9059
40 Zr 1.290000 2.210000 0 5 91.22
41 Nb 1.290000 1.700000 0 5 92.9064
42 Mo 1.240000 2.060000 0 5 95.94
43 Tc 1.250000 1.700000 0 5 98
44 Ru 0.000000 2.010000 0 5 101.07
45 Rh 1.340000 2.010000 0 5 102.9055
46 Pd 1.410000 2.040000 0 5 106.4
47 Ag 1.500000 1.700000 0 5 107.868
48 Cd 1.400000 1.700000 0 5 112.41
49 In 1.410000 1.700000 0 5 114.82
50 Sn 1.370000 1.700000 0 5 118.69
51 Sb 1.330000 1.700000 0 5 121.75
52 Te 0.700000 1.700000 0 5 127.60
53 I 1.330000 2.150000 1 6 126.9045
54 Xe 1.980000 1.700000 0 5 131.30
55 Cs 1.690000 1.700000 0 5 132.9054
56 Ba 1.690000 1.700000 0 5 137.33
57 La 0.000000 0.800000 0 5 138.9055
58 Ce 1.690000 1.700000 0 5 140.12
59 Pr 1.690000 1.700000 0 5 140.9077
60 Nd 1.690000 1.700000 0 5 144.24
61 Pm 1.690000 1.700000 0 5 145
62 Sm 1.690000 1.700000 0 5 150.4
63 Eu 1.690000 1.700000 0 5 151.96
64 Gd 1.690000 1.700000 0 5 157.25
65 Tb 1.690000 1.700000 0 5 158.9254
66 Dy 1.690000 1.700000 0 5 162.50
67 Ho 1.690000 1.700000 0 5 164.9304
68 Er 1.690000 1.700000 0 5 167.26
69 Tm 1.690000 1.700000 0 5 168.9342
70 Yb 1.690000 1.700000 0 5 173.04
71 Lu 1.690000 1.700000 0 5 174.967
72 Hf 1.440000 1.700000 0 5 178.49
73 Ta 1.340000 1.700000 0 5 180.9479
74 W 1.300000 1.700000 0 5 183.85
75 Re 1.280000 1.700000 0 5 186.207
76 Os 1.260000 2.020000 0 5 190.2
77 Ir 1.290000 2.030000 0 5 192.22
78 Pt 1.340000 1.700000 0 5 195.09
79 Au 1.440000 1.700000 0 9 196.9665
80 Hg 1.550000 1.700000 0 5 200.59
81 Tl 1.540000 1.700000 0 5 204.37
82 Pb 1.520000 1.700000 0 5 207.2
83 Bi 1.520000 1.700000 0 5 208.9804
84 Po 1.400000 1.700000 0 5 209
85 At 0.700000 1.700000 0 5 210
86 Rn 2.400000 1.700000 0 5 222
87 Fr 2.000000 1.700000 0 5 223
88 Ra 1.900000 1.700000 0 5 226.0254
89 Ac 1.900000 1.700000 0 5 227.0278
90 Th 1.900000 1.700000 0 5 232.0381
91 Pa 1.900000 1.700000 0 5 231.0359
92 U 1.900000 1.700000 0 5 238.029
93 Np 2.0 2.0 0 5 237.0482
94 Pu 2.0 2.0 0 5 244
95 Am 2.0 2.0 0 5 243
96 Cm 2.0 2.0 0 5 247
97 Bk 2.0 2.0 0 5 247
98 Cf 2.0 2.0 0 5 251
99 Es 2.0 2.0 0 5 254
100 Fm 2.0 2.0 0 5 257
101 Md 2.0 2.0 0 5 258
102 No 2.0 2.0 0 5 259
103 Lr 2.0 2.0 0 5 260
"""
atomData = 110*[None]
for line in atomDataText.splitlines():
try:
charge, symbol, covalR, vdwR, valence, color, mass = line.split()
at = atomDataElement(symbol,charge,covalR,vdwR,valence,color,mass)
atomData[int(at.charge)] = at
except:
pass
if __name__ == '__main__':
for s in atomData:
print s

484
resultsFile/lib/basis.py Executable file
View File

@ -0,0 +1,484 @@
#!/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
#import pdb
from library import *
from math import *
import sys
def normalize_basis_name(s):
l = list(s)
is_sphe = False
l.sort()
for i in range(10):
if str(i) in l:
is_sphe = True
if is_sphe:
if l[0] == '0':
l.insert(0,'+')
l = [l[-1]]+l[:-1]
result = ""
for i in l:
result += i
result = result.lower()
return result
class primitive(object):
"""Class for a primitive basis function."""
def __init__(self):
_center = None
_expo = None
_sym = None
def __eq__(self,other):
thr = 1.e-6
result = True
result = result and abs(self.center[0] - other.center[0]) < thr
result = result and abs(self.center[1] - other.center[1]) < thr
result = result and abs(self.center[2] - other.center[2]) < thr
result = result and abs(self.expo - other.expo) < thr
result = result and self.sym == other.sym
return result
def __repr__(self):
out = "%6s %10.6f %10.6f %10.6f %16.6e"%tuple(
[self.sym]+list(self.center)+[self.expo])
return out
def __repr__debug__(self):
out = ""
out += str(self.sym)+'\t'
out += str(self.center)+'\t'
out += str(self.expo)
return out
def __cmp__(self,other):
assert ( isinstance(other,primitive) )
"""Primitive functions are sorted according to the exponents."""
if not isinstance(other,primitive):
raise TypeError
if self.expo < other.expo:
return -1
elif self.expo > other.expo:
return 1
elif self.expo == other.expo:
return 0
def get_norm(self):
result = sqrt(self.overlap(self))
return result
norm = property (get_norm,None,doc="Sqrt( Integral f^2(R) dR ).")
for i in "center expo sym".split():
exec """
def get_%(i)s(self): return self._%(i)s
def set_%(i)s(self,value): self._%(i)s = value
%(i)s = property(fget=get_%(i)s,fset=set_%(i)s) """%locals()
#-----------------------
powersave = { 's':(0,0,0) }
def powers(sym):
if sym in powersave:
return powersave[sym]
result = (sym.count('x'),sym.count('y'),sym.count('z'))
powersave[sym] = result
return result
fact_ = [1.]
def fact(n):
global fact_
nstart = len(fact_)
if n < nstart :
return fact_[n]
else:
for i in range(nstart,n+1):
fact_.append(fact_[i-1]*float(i))
return fact_[n]
def binom(n,m):
return fact(n)/(fact(m)*fact(n-m))
def rintgauss(n):
def ddfact2(n):
if n%2 == 0: print 'error in ddfact2'
res=1.
for i in range(1,n+1,2):
res*=float(i)
return res
res = sqrt(pi)
if n == 0: return res
elif n == 1: return 0.
elif n%2 == 1: return 0.
res /= 2.**(n/2)
res *= ddfact2(n-1)
return res
def Goverlap(fA,fB):
fA = spheToCart(fA)
fB = spheToCart(fB)
if isinstance(fA,contraction):
result = fA.overlap(fB)
elif isinstance(fB,contraction):
result = fB.overlap(fA)
else:
result = GoverlapCart(fA,fB)
return result
def GoverlapCart(fA,fB):
gamA=fA.expo
gamB=fB.expo
gamtot = gamA+gamB
SAB=1.0
A = fA.center
B = fB.center
nA = powers(fA.sym)
nB = powers(fB.sym)
for l in range(3):
Al = A[l]
Bl = B[l]
nAl = nA[l]
nBl = nB[l]
u=gamA/gamtot*Al+gamB/gamtot*Bl
arg=gamtot*u*u-gamA*Al*Al-gamB*Bl*Bl
alpha=exp(arg)/gamtot**((1.+float(nAl)+float(nBl))/2.)
temp = sqrt(gamtot)
wA=temp*(u-Al)
wB=temp*(u-Bl)
accu=0.
for n in range (nAl+1):
wAn = wA**n * binom(nAl,n)
for m in range (nBl+1):
integ=nAl+nBl-n-m
accu+=wAn*wB**m*binom(nBl,m)*rintgauss(integ)
SAB*=accu*alpha
return SAB
def GoverlapCartNorm2(fA,fB):
gamA=fA.expo
gamB=fB.expo
gamtot = gamA+gamB
SAB=1.0
nA = powers(fA.sym)
nB = powers(fB.sym)
for l in range(3):
nAl = nA[l]
nBl = nB[l]
SAB*=rintgauss(nAl+nBl)/(gamA+gamB)**((1.+float(nAl)+float(nBl))/2.)
return SAB
def GoverlapCartNorm(fA):
gamA=fA.expo
SAB=1.0
nA = powers(fA.sym)
for l in range(3):
nAl = nA[l]
SAB*=rintgauss(2*nAl)/(2.*gamA)**(0.5+float(nAl))
return SAB
angular_momentum = {}
for i,l in enumerate("spdfghijklmno"):
angular_momentum[l]=i
def get_lm(sym):
if 'x' in sym or 'y' in sym or 'z' in sym or sym == 's':
return None, None
else:
return angular_momentum[sym[0]], int(sym[1:])
def xyz_from_lm(l,m):
"""Returns the xyz powers and the coefficients of a spherical function
expressed in the cartesian basis"""
power = []
coef = []
absm = abs(m)
nb2 = absm
nb1 = (l-absm)/2
clmt = [ (-0.25)**t * binom(l,t) * binom(l-t,absm+t) for t in xrange(nb1+1) ]
mod_absm_2 = absm % 2
if m>=0:
nb2_start = mod_absm_2
elif m<0:
nb2_start = (absm+1) % 2
if m != 0:
norm = 2.**(-absm)/fact(l) * sqrt(2.*abs(fact(l+m)*fact(l-m)))
else:
norm = 1.
for n1 in range(nb2_start,nb2+1,2):
k = (absm-n1)/2
factor = (-1.)**k * binom(absm,n1) * norm
for t in xrange(nb1+1):
for n2 in xrange(t+1):
coe = clmt[t] * factor * binom(t,n2)
ipw = ( n1+2*n2 , absm-n1+2*(t-n2), l-2*t-absm )
done = False
for i,c in enumerate(coef):
if c != 0.:
if not done and ipw in power:
idx = power.index(ipw)
coef[idx] += coe
done = True
if not done:
power.append(ipw)
coef.append(coe)
return power,coef
def spheToCart(fA):
l,m = get_lm(fA.sym)
if l is None:
return fA
power, coef = xyz_from_lm(l,m)
contr = contraction()
for p,c in zip(power,coef):
gauss = gaussian()
gauss.center = fA.center
gauss.expo = fA.expo
gauss.sym = ''
for l,letter in enumerate('xyz'):
gauss.sym += p[l]*letter
contr.append(c,gauss)
return contr
class gaussian(primitive):
"""Gaussian primitive function."""
def __init__(self):
primitive.__init__(self)
def overlap(self,other):
result = Goverlap(self,other)
return result
def value(self,r):
"""Value at r."""
x, y, z = r
x -= self.center[0]
y -= self.center[1]
z -= self.center[2]
r2 = x**2 + y**2 + z**2
px, py, pz = powers(self.sym)
P = x**px * y**py * z**pz
return P*exp(-self.expo*r2)
def valuer2(self,r2):
"""Value at sqrt(r2)."""
return exp(-self.expo*r2)
#-----------------------
class contraction(object):
"""Contraction of primitive functions."""
def __init__(self):
self._center = None
self._prim = []
self._coef = []
self._sym = None
def __eq__(self,other):
thr = 1.e-6
result = True
result = result and abs(self._center[0] - other._center[0]) < thr
result = result and abs(self._center[1] - other._center[1]) < thr
result = result and abs(self._center[2] - other._center[2]) < thr
result = result and self._sym == other._sym
result = result and len(self._coef) == len(other._coef)
if result:
for i in range(len(self._coef)):
result = result and abs(self._coef[i] - other._coef[i]) < thr
if result:
for i in range(len(self._prim)):
result = result and self._prim[i] == other._prim[i]
return result
def __repr__(self):
out = "%6s %10.6f %10.6f %10.6f\n"%tuple(
[self.sym]+list(self.center))
for i,a in enumerate(self.prim):
out += " %16.6e %16.6e\n"%(a.expo,self.coef[i])
return out
def append(self,c,func):
if self.sym is None:
self.sym = func.sym
if self.center is None:
self.center = func.center
assert ( isinstance(func,primitive) )
assert ( self.center == func.center )
self._prim.append(func)
self._coef.append(c)
def get_norm(self):
sum=0.
for i, ci in enumerate(self.coef):
for j, cj in enumerate(self.coef):
sum += ci*cj*self.prim[i].overlap(self.prim[j])/(self.prim[i].norm*self.prim[j].norm)
result = sqrt(sum)
return result
def overlap(self,other0):
if isinstance(other0,contraction):
other=other0
else:
other=contraction()
other.sym = other0.sym
other.center = other0.center
other.append(1.,other0)
prim = self.prim
oprim = other.prim
sum = 0.
for i, ci in enumerate(self.coef):
for j, cj in enumerate(other.coef):
sum += ci*cj*prim[i].overlap(oprim[j])
return sum
def get_coef(self):
return self._coef
coef = property(fget=get_coef)
def get_center(self):
return self._center
center = property(fget=get_center)
def get_prim(self):
return self._prim
prim = property(fget=get_prim)
def normalize(self):
coef = self.coef
prim = self.prim
n = self.norm
for i, ci in enumerate(coef):
coef[i] /= n
def value(self,r):
"""Value at r."""
x, y, z = r
x -= self.center[0]
y -= self.center[1]
z -= self.center[2]
r2 = x**2 + y**2 + z**2
px, py, pz = powers(self.sym)
P = x**px * y**py * z**pz
result = 0.
for c,chi in zip(self.coef,self.prim):
result += c*chi.valuer2(r2)
return P*result
for i in "center prim sym coef".split():
exec """
def get_%(i)s(self): return self._%(i)s
def set_%(i)s(self,value): self._%(i)s = value
%(i)s = property(fget=get_%(i)s,fset=set_%(i)s) """%locals()
norm = property (get_norm,None,doc="Sqrt( Integral f^2(R) dR ).")
try:
import lib_cython
fact = lib_cython.fact
binom = lib_cython.binom
rintgauss = lib_cython.rintgauss
powers = lib_cython.powers
GoverlapCart = lib_cython.GoverlapCart
GoverlapCartNorm2 = lib_cython.GoverlapCartNorm2
except ImportError:
pass
#--------------
if __name__ == '__main__':
for l in range(6):
print '---'
for m in range(-l,l+1):
print "%3d %3d :"%(l,m), xyz_from_lm(l,m)
sys.exit(0)
if __name__ == '__main__':
a = gaussian()
a.sym = 'f+0'
a.expo = 3.
a.center = [0.,0.,0.]
b = gaussian()
b.sym = 'g+2'
b.expo = 4.
b.center = [0.,0.,0.]
x = spheToCart(a)
y = spheToCart(a)
print x
sys.exit(0)
print GoverlapCart(a,b)
print GoverlapCartNorm2(a,b)
print ''
for i in range(0,10):
print rintgauss(i)
sys.exit(0)
c = gaussian()
c.sym = 's'
c.expo = 1.159
c.center = [0.,0.,0.]
d = contraction()
d.append(0.025494863235, a)
d.append(0.190362765893, b)
d.append(0.852162022245, c)
d.normalize()
print d
sys.exit(0)
a = gaussian()
a.sym = 'x'
a.expo = 1.0
a.center = [0.,0.,0.]
b = gaussian()
b.sym = 'x'
b.expo = 2.0
b.center = [0.,0.,0.]
c = contraction()
c.append(1.,a)
c.append(2.,b)
print c.value( (1.,1.,0.) )
print b.overlap(a)
print a.overlap(b)
print a.norm
print c.norm, c.coef
c.normalize()
print c.norm, c.coef
c.normalize()
print c.norm, c.coef
print Goverlap(a,b)

76
resultsFile/lib/csf.py Executable file
View File

@ -0,0 +1,76 @@
#!/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
from library import *
class CSF(object):
"""Class for an configuration state function."""
def __init__(self):
self._determinants = []
self._coefficients = []
def __repr__(self):
out = ""
for d,c in zip(self.determinants,self.coefficients):
out += "%10.8f\n"%(c,)
for spin in d:
for orb in d[spin]:
out += "%4d "%(orb,)
out += "\n"
return out
def __repr__debug__(self):
out = ""
out += "CSF:\n"
out += " Determinants: "+str(self.determinants)+'\n'
out += " Coefficients: "+str(self.coefficients)+'\n'
return out
def __cmp__(self,other):
assert ( isinstance(other,CSF) )
if len(self.coefficients) < len(other.coefficients):
return -1
elif len(self.coefficients) > len(other.coefficients):
return 1
elif len(self.coefficients) == len(other.coefficients):
return 0
def append(self,c,up,dn):
det = {}
det['alpha'] = up
det['beta'] = dn
self._determinants.append( det )
self._coefficients.append( c )
for i in "determinants coefficients".split():
exec """
def get_%(i)s(self): return self._%(i)s
def set_%(i)s(self,value): self._%(i)s = value
%(i)s = property(fget=get_%(i)s,fset=set_%(i)s) """%locals()

51
resultsFile/lib/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("lib_cython", ["lib_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
)

119
resultsFile/lib/fortranBinary.py Executable file
View File

@ -0,0 +1,119 @@
#!/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 struct
from library import *
import sys
class fortranBinary(object):
"""
Class for a fortran binary file.
Behaves like an array of records.
"""
def __init__(self,name,mode):
self._name = name
self._mode = mode
self._file = None
self._recl = None
self._form = None
def get_file(self):
if self._file is None:
self._file = open(self._name,self._mode)
return self._file
def seek(self,index):
self.file.seek(index)
def next(self):
if self.form is None:
raise TypeError
data = self.file.read(4) # Rec length at the beginning of record
if not data: raise IndexError
self.recl = struct.unpack('l',data)[0]
data = self.file.read(self._recl)
data = struct.unpack(self.form,data)
self.file.read(4)
return list(data)
def __iter__(self):
if self.form is None:
raise TypeError
self.file.seek(0)
while True:
data = self.file.read(4) # Rec length at the beginning of record
if not data: raise StopIteration
self.recl = struct.unpack('l',data)[0]
data = self.file.read(self._recl)
data = struct.unpack(self.form,data)
self.file.read(4)
yield list(data)
def __getitem__(self,index):
if self.form is None:
raise TypeError
self.file.seek(0)
for i in range(index-1):
data = self.file.read(4) # Rec length at the beginning of record
if not data: raise IndexError
recl = struct.unpack('l',data)[0]
data = self.file.read(recl+4)
data = self.file.read(4) # Rec length at the beginning of record
if not data: raise IndexError
self.recl = struct.unpack('l',data)[0]
data = self.file.read(self._recl)
data = struct.unpack(self.form,data)
self.file.read(4)
return list(data)
def close(self):
self._file.close()
# recl = xproperty ( '_recl', doc="Record length")
# mode = xproperty ( '_mode', doc="R/W mode")
# name = xproperty ( '_name', doc="R/W mode")
# file = property ( get_file, doc="Binary file")
# form = xproperty ( '_form', doc="Format to read the binary record")
for i in "recl mode name file form".split():
exec """
def get_%(i)s(self): return self._%(i)s
def set_%(i)s(self,value): self._%(i)s = value
%(i)s = property(fget=get_%(i)s,fset=set_%(i)s) """%locals()
if __name__ == '__main__':
f = fortranBinary(sys.argv[1],"rb")
f.form = 'l'+4*15000*'h'+15000*'d'
print f.name
print f.mode
print f[0]
print f[10]
print f[0]

74
resultsFile/lib/integral.py Executable file
View File

@ -0,0 +1,74 @@
#!/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 library import *
from math import *
class integral(object):
"""Class for an integral."""
def __init__(self):
self._indices = None
self._value = None
def __repr__(self):
out = ""
i = self.indices
if len(i) == 4:
out += "( %4d %4d | %4d %4d )"% tuple( [int(i[k]) for k in range(4)] )
elif len(i) == 2:
out += "( %4d | h | %4d )"% tuple( [int(i[k]) for k in range(2)] )
out += " : %24.16e" % ( self.value )
return out
def __cmp__(self,other):
assert ( isinstance(other,integral) )
if self._value < other._value:
return -1
elif self._value > other._value:
return 1
elif self._value == other._value:
return 0
for i in "indices value".split():
exec """
def get_%(i)s(self): return self._%(i)s
def set_%(i)s(self,value): self._%(i)s = value
%(i)s = property(fget=get_%(i)s,fset=set_%(i)s) """%locals()
if __name__ == '__main__':
i = integral()
i.indices = [1,4]
i.value = 1.5
print i
j = integral()
j.indices = [4,3,2,1]
j.value = 2.5
print j
print i<j

89
resultsFile/lib/lib_cython.pyx Normal file
View File

@ -0,0 +1,89 @@
#!/usr/bin/python
from math import *
powersave = { 's':(0,0,0) }
def powers(sym):
if sym in powersave:
return powersave[sym]
result = (sym.count('x'),sym.count('y'),sym.count('z'))
powersave[sym] = result
return result
fact_ = [1.]
cpdef fact(int n):
global fact_
cdef int nstart = len(fact_)
cdef int i
if n >= nstart :
for i in range(nstart,n+1):
fact_.append(float(fact_[i-1]*i))
return fact_[n]
def binom(int n,int m):
return fact(n)/(fact(m)*fact(n-m))
cdef ddfact2(int n):
if n%2 == 0: print 'error in ddfact2'
cdef double res=1.
cdef int i
for i in range(1,n+1,2):
res*=float(i)
return res
cdef double sqpi = sqrt(pi)
cpdef rintgauss(int n):
res = sqpi
if n == 0: return res
elif n == 1: return 0.
elif n%2 == 1: return 0.
res /= 2.**(n/2)
res *= ddfact2(n-1)
return res
cpdef GoverlapCart(fA,fB):
cdef double gamA=fA.expo
cdef double gamB=fB.expo
cdef double gamtot = gamA+gamB
cdef double SAB=1.0
cdef int l, n, m
cdef double u, arg, alpha, temp, wA, wB, accu
cdef int integ
A = fA.center
B = fB.center
nA = powers(fA.sym)
nB = powers(fB.sym)
for l in range(3):
Al = A[l]
Bl = B[l]
nAl = nA[l]
nBl = nB[l]
u=gamA/gamtot*Al+gamB/gamtot*Bl
arg=gamtot*u*u-gamA*Al*Al-gamB*Bl*Bl
alpha=exp(arg)/gamtot**((1.+float(nAl)+float(nBl))*0.5)
temp = sqrt(gamtot)
wA=temp*(u-Al)
wB=temp*(u-Bl)
accu=0.
for n in range (nAl+1):
for m in range (nBl+1):
integ=nAl+nBl-n-m
accu+=wA**n*wB**m*binom(nAl,n)*binom(nBl,m)*rintgauss(integ)
SAB*=accu*alpha
return SAB
cpdef GoverlapCartNorm2(fA,fB):
cdef double gamA=fA.expo
cdef double gamB=fB.expo
cdef double gamtot = gamA+gamB
cdef double SAB=1.0
cdef int l
nA = powers(fA.sym)
nB = powers(fB.sym)
for l in range(3):
nAl = nA[l]
nBl = nB[l]
SAB*=rintgauss(nAl+nBl)/(gamA+gamB)**((1.+float(nAl)+float(nBl))/2.)
return SAB

132
resultsFile/lib/library.py Executable file
View File

@ -0,0 +1,132 @@
#!/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 math import *
a0 = 0.529177249
def prettyPrint(obj,shift=" "):
if isinstance(obj,list):
for k in obj:
prettyPrint(k)
elif isinstance(obj,dict):
for k in obj:
print k
prettyPrint(obj[k])
else:
print str(obj)
def xproperty(fget, fset=None, fdel=None, doc=None):
"""Automatically defines the properties if the functions don't exist."""
if isinstance(fget,str):
attr_name = fget
if fset is None:
fset = fget
fget = lambda obj: getattr(obj, attr_name)
if isinstance(fset,str):
attr_name = fset
fset = lambda obj, val: setattr(obj, attr_name, val)
return property(fget, fset, fdel, doc)
def build_get_funcs(var):
line = ""
line += """def get_"""+var+"""(self):
try:
getattr(self,'_"""+var+"""')
except AttributeError:
self._"""+var+""" = None
return getattr(self, '_"""+var+"""')\n"""
return line
def build_property(var,txt):
line = "def _get_"+var+"(self):\n"
line += " return self.get_"+var+"()\n"
line += var+" = property(_get_"+var+", fset=None, fdel=None, doc='"+txt+"')\n"
return line
def get_data(var,txt,pos,type=""):
line = """def get_%(var)s(self):
if self._%(var)s is None:
try:
self.find_string("%(txt)s")
except IndexError:
return None
pos = self._pos
if pos is not None:
line = self.text[pos].split()"""
if type is "":
line += """
self._%(var)s = ' '.join(line[%(pos)s])"""
else:
line += """
self._%(var)s = """+type+"""(' '.join(line[%(pos)s]))"""
line += """
return self._%(var)s"""
line = line % {'var': var, 'txt': txt, 'pos': pos}
return line
def init_variable(var,value):
assert ( value is not None )
line = """def get_"""+var+"""(self):
try:
getattr(self,'_"""+var+"""')
except AttributeError:
self._"""+var+""" = None
if self._"""+var+""" is None:
self._"""+var+""" = """+str(value)+"""
return self._"""+var
return line
def cartesianToSpherical(cart):
x, y, z = cart
r = sqrt(x**2+y**2+z**2)
sign = 1.
if x != 0.:
sign = abs(x)/x
theta = atan(y/x)
else:
sign = 1.
theta = atan(pi/2.)
if r != 0.:
phi = acos(z/r)*sign
else:
phi = 0.
return [r,theta,phi]
def sphericalToCartesian(sphe):
r, theta, phi = sphe
x = r*cos(theta)*sin(phi)
y = r*sin(theta)*sin(phi)
z = r*cos(phi)
return [x,y,z]

133
resultsFile/lib/orbital.py Executable file
View File

@ -0,0 +1,133 @@
#!/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
from library import *
from math import *
class orbital(object):
"""Class for an orbital."""
def __init__(self):
self._set = None
self._eigenvalue = None
self._vector = []
self._basis = None
self._sym = None
def __repr__(self):
out = "%10s %6s %10.6f\n"%(\
self.set, self.sym, self.eigenvalue)
count=0
for c in self.vector:
out += "%17.8e "%(c,)
count += 1
if count == 5:
count = 0
out += '\n'
return out
def __repr__debug__(self):
out = ""
out += "Orbital:\n"
out += " Set : "+str(self.set)+'\n'
out += " Sym : "+str(self.sym)+'\n'
out += " Basis : "+str(self.basis)+'\n'
out += " Eigenvalue : "+str(self.eigenvalue)+'\n'
out += " Vector : "+str(self.vector)+'\n'
return out
def __cmp__(self,other):
assert ( isinstance(other,orbital) )
if self.eigenvalue < other.eigenvalue:
return -1
elif self.eigenvalue > other.eigenvalue:
return 1
elif self.eigenvalue == other.eigenvalue:
return 0
def dot(self,other):
assert ( isinstance(other,orbital) )
result = 0.
for i,j in zip(self.vector,other.vector):
result += i*j
return sqrt(result)
def is_ortho(self,other):
assert ( isinstance(other,orbital) )
if self.dot(other) == 0. :
return True
else:
return False
def norm(self):
return self.dot(self)
def value(self,r):
"""Value at r."""
x, y, z = r
result = 0.
for c,chi in zip(self.vector,self.basis):
result += c*chi.value(r)
return result
for i in "eigenvalue vector basis set sym".split():
exec """
def get_%(i)s(self): return self._%(i)s
def set_%(i)s(self,value): self._%(i)s = value
%(i)s = property(fget=get_%(i)s,fset=set_%(i)s) """%locals()
if __name__ == '__main__':
l = []
l.append(orbital())
l.append(orbital())
l.append(orbital())
l[0].eigenvalue = 2.
l[1].eigenvalue = 1.
l[2].eigenvalue = 3.
for i in l:
print i.eigenvalue
l.sort()
print ""
for i in l:
print i.eigenvalue
l[0].vector = [1., 0., 0.]
l[1].vector = [0., 1., 0.]
l[2].vector = [0., 1., 1.]
print l[0].is_ortho(l[1])
print l[1].is_ortho(l[2])
print l[1].dot(l[2])
print l[1].dot(l[0])
print l[1].norm()
print l[2].norm()