From a88a846bea7eab6b48eb1c3304c20fc99b33d7e3 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 19 Jan 2015 17:04:24 +0100 Subject: [PATCH] Forgot lib --- .gitignore | 8 - resultsFile/lib/__init__.py | 46 +++ resultsFile/lib/atom.py | 217 ++++++++++++++ resultsFile/lib/basis.py | 484 +++++++++++++++++++++++++++++++ resultsFile/lib/csf.py | 76 +++++ resultsFile/lib/cython_setup | 51 ++++ resultsFile/lib/fortranBinary.py | 119 ++++++++ resultsFile/lib/integral.py | 74 +++++ resultsFile/lib/lib_cython.pyx | 89 ++++++ resultsFile/lib/library.py | 132 +++++++++ resultsFile/lib/orbital.py | 133 +++++++++ 11 files changed, 1421 insertions(+), 8 deletions(-) create mode 100755 resultsFile/lib/__init__.py create mode 100755 resultsFile/lib/atom.py create mode 100755 resultsFile/lib/basis.py create mode 100755 resultsFile/lib/csf.py create mode 100755 resultsFile/lib/cython_setup create mode 100755 resultsFile/lib/fortranBinary.py create mode 100755 resultsFile/lib/integral.py create mode 100644 resultsFile/lib/lib_cython.pyx create mode 100755 resultsFile/lib/library.py create mode 100755 resultsFile/lib/orbital.py diff --git a/.gitignore b/.gitignore index db4561e..1cd79fd 100644 --- a/.gitignore +++ b/.gitignore @@ -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 diff --git a/resultsFile/lib/__init__.py b/resultsFile/lib/__init__.py new file mode 100755 index 0000000..ec6d1ec --- /dev/null +++ b/resultsFile/lib/__init__.py @@ -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 " +__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 + + diff --git a/resultsFile/lib/atom.py b/resultsFile/lib/atom.py new file mode 100755 index 0000000..d99037f --- /dev/null +++ b/resultsFile/lib/atom.py @@ -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 + diff --git a/resultsFile/lib/basis.py b/resultsFile/lib/basis.py new file mode 100755 index 0000000..20a0c55 --- /dev/null +++ b/resultsFile/lib/basis.py @@ -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) + diff --git a/resultsFile/lib/csf.py b/resultsFile/lib/csf.py new file mode 100755 index 0000000..a0edac9 --- /dev/null +++ b/resultsFile/lib/csf.py @@ -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() + diff --git a/resultsFile/lib/cython_setup b/resultsFile/lib/cython_setup new file mode 100755 index 0000000..4ca692e --- /dev/null +++ b/resultsFile/lib/cython_setup @@ -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 + ) + diff --git a/resultsFile/lib/fortranBinary.py b/resultsFile/lib/fortranBinary.py new file mode 100755 index 0000000..91f0c15 --- /dev/null +++ b/resultsFile/lib/fortranBinary.py @@ -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] diff --git a/resultsFile/lib/integral.py b/resultsFile/lib/integral.py new file mode 100755 index 0000000..81aac93 --- /dev/null +++ b/resultsFile/lib/integral.py @@ -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= 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 + diff --git a/resultsFile/lib/library.py b/resultsFile/lib/library.py new file mode 100755 index 0000000..680a940 --- /dev/null +++ b/resultsFile/lib/library.py @@ -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] + + diff --git a/resultsFile/lib/orbital.py b/resultsFile/lib/orbital.py new file mode 100755 index 0000000..78e7648 --- /dev/null +++ b/resultsFile/lib/orbital.py @@ -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() + + +