diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..010a68a --- /dev/null +++ b/Makefile @@ -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 diff --git a/resultsFile/Modules/__init__.py b/resultsFile/Modules/__init__.py new file mode 100755 index 0000000..4baa72e --- /dev/null +++ b/resultsFile/Modules/__init__.py @@ -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 " +__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 + + diff --git a/resultsFile/Modules/gamessFile.py b/resultsFile/Modules/gamessFile.py new file mode 100755 index 0000000..7734f67 --- /dev/null +++ b/resultsFile/Modules/gamessFile.py @@ -0,0 +1,1565 @@ +#!/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 + +gamessFile_defined_vars = [ "date", "version", "machine", "memory", "disk",\ + "cpu_time", "author", "title", "units", "methods", "options", \ + "spin_restrict", "conv_threshs", "energies", \ + "one_e_energies", "two_e_energies", "ee_pot_energies", \ + "Ne_pot_energies", "pot_energies", \ + "kin_energies", "virials", "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", "two_e_int_ao_filename", + "one_e_int_ao_filename", "atom_to_ao_range", "gradient_energy" ] + +class gamessFile(resultsFile): + """ Class defining the gamess file. + """ + + local_vars = list(local_vars) + defined_vars = list(gamessFile_defined_vars) + + exec get_data('date',"EXECUTION OF GAMESS BEGUN",'4:') in locals() + exec get_data('machine',"Files used on the master node",'6:7') in locals() + exec get_data('version',"GAMESS VERSION",'4:8') in locals() + exec get_data('num_proc',"MEMDDI DISTRIBUTED OVER",'3:4',"int") in locals() + exec get_data('num_elec',"NUMBER OF ELECTRONS", '4:5',"int") in locals() + exec get_data('charge',"CHARGE OF MOLECULE", '4:5',"float") in locals() + exec get_data('nuclear_energy',"THE NUCLEAR REPULSION ENERGY", '5:6',"float") in locals() + + def get_multiplicity(self): + if self._multiplicity is None: + try: + self.find_string("SPIN MULTIPLICITY") + except IndexError: + try: + self.find_string("STATE MULTIPLICITY") + except IndexError: + return None + pos = self._pos + if pos is not None: + line = self.text[pos].split() + self._multiplicity = int(' '.join(line[3:4])) + return self._multiplicity + + def get_two_e_int_mo_filename(self): + try: + getattr(self,'_two_e_int_mo_filename') + except AttributeError: + self._two_e_int_mo_filename = None + if self._two_e_int_mo_filename is None: + filename = self.filename.split('.') + filename.pop() + filename.append("F08") + filename = '.'.join(filename) + self._two_e_int_mo_filename = filename + return self._two_e_int_mo_filename + + def set_two_e_int_mo_filename(self,value): + self._two_e_int_mo_filename = value + + def get_two_e_int_ao_filename(self): + try: + getattr(self,'_two_e_int_ao_filename') + except AttributeError: + self._two_e_int_ao_filename = None + if self._two_e_int_ao_filename is None: + filename = self.filename.split('.') + filename.pop() + filename.append("F08") + filename = '.'.join(filename) + self._two_e_int_ao_filename = filename + return self._two_e_int_ao_filename + + def set_two_e_int_ao_filename(self,value): + self._two_e_int_ao_filename = value + + def get_disk(self): + if self._disk is None: + try: + self.find_string("Files used on the master node") + except IndexError: + return None + pos = self._pos + pos += 1 + line = self.text[pos].split() + disk = 0 + while len(line) == 8: + disk += float(line[4]) + pos += 1 + line = self.text[pos].split() + disk = disk/1000000000. + 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_cpu_time(self): + if self._cpu_time is None: + try: + self.find_string("EXECUTION OF GAMESS TERMINATED NORMALLY") + except IndexError: + return None + pos = self._pos + pos -= 3 + line = self.text[pos].split() + self._cpu_time = line[-4]+" s" + return self._cpu_time + + def get_memory(self): + if self._memory is None: + try: + self.find_string("EXECUTION OF GAMESS TERMINATED NORMALLY") + except IndexError: + return None + pos = self._pos + pos -= 1 + line = self.text[pos].split() + memory = float(line[0])*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_author(self): + if self._author is None: + try: + self.find_string("Files used on the master node") + except IndexError: + return None + pos = self._pos + pos += 1 + line = self.text[pos].split() + self._author = line[2] + return self._author + + def get_title(self): + if self._title is None: + try: + self.find_string("RUN TITLE") + except IndexError: + return None + pos = self._pos + pos += 2 + self._title = self.text[pos].strip() + return self._title + + def get_symmetries(self): + if self._symmetries is None: + try: + self.find_string("DIMENSIONS OF THE SYMMETRY SUBSPACES") + except IndexError: + return None + pos = self._pos + begin = pos + 1 + try: + self.find_next_string("DONE") + except IndexError: + return None + end = self._pos-1 + sym = [] + for k in range(begin,end): + buffer = self.text[k].replace('=',' ').split() + for i,j in zip(buffer[0::2],buffer[1::2]): + sym.append([i,int(j)]) + self._symmetries = sym + return self._symmetries + + def get_units(self): + if self._units is None: + try: + self.find_string("ATOM ATOMIC COORDINATES") + except IndexError: + return None + pos = self._pos + line = list(self.text[pos].split().pop()) + line.pop() + line.pop(0) + self._units = ''.join(line) + return self._units + + def get_methods(self): + if self._methods is None: + methods = [] + options = self.options + for m in ['SCFTYP', 'CITYP']: + if options[m] != 'NONE': + methods.append(options[m]) + self._methods = methods + return self._methods + + def get_occ_num(self): + if self._occ_num is None: + occ = {} + for motype in self.mo_types: + occ[motype] = [ 0. for mo in self.mo_sets[motype] ] + if motype == "Natural": + for i,mo in enumerate(self.mo_sets[motype]): + occ[motype][i] = mo.eigenvalue + else: + if self.mulliken_mo is not None: + for i,mo in enumerate(self.mulliken_mo): + occ[motype][i] = mo.eigenvalue + if occ != {}: + self._occ_num = occ + return self._occ_num + + def get_options(self): + if self._options is None: + options = {} + list_of_keys = ["$CONTRL OPTIONS","INTEGRAL TRANSFORMATION OPTIONS", \ + "INTEGRAL INPUT OPTIONS","GUESS OPTIONS"] + for key in list_of_keys: + try: + self.find_string(key) + pos = self._pos + 2 + line = self.text[pos].replace('=',' ') + tokens = line.split()[0::2] + values = line.split()[1::2] + while len(tokens) > 0: + for i,t in enumerate(tokens): + options[t] = values[i] + pos += 1 + line = self.text[pos].replace('FRIEND= ','FRIEND=NONE') + line = line.replace('=',' ') + tokens = line.split()[0::2] + values = line.split()[1::2] + if len(tokens) != len(values): + print "error: ",line + sys.exit(0) + except IndexError: + pass + try: + self.find_string("$SYSTEM OPTIONS") + pos = self._pos + 2 + line = self.text[pos].split() + options['MEMORY'] = line[2] + pos += 1 + line = self.text[pos].split() + options['MEMDDI'] = line[2] + pos += 3 + line = self.text[pos].replace('=',' ').split() + options[line[0]] = line[1] + pos += 1 + line = self.text[pos].replace('=',' ') + tokens = line.split()[0::2] + values = line.split()[1::2] + for i,t in enumerate(tokens): + options[t] = values[i] + except IndexError: + pass + try: + self.find_string("PROPERTIES INPUT") + pos = self._pos + 4 + line = self.text[pos].replace('=',' ') + tokens = line.split()[0::2] + values = line.split()[1::2] + while tokens[0] != "EXTRAPOLATION": + for i,t in enumerate(tokens): + options[t] = values[i] + pos += 1 + line = self.text[pos].replace('=',' ') + tokens = line.split()[0::2] + values = line.split()[1::2] + except IndexError: + pass + try: + self.find_string(" SCF CALCULATION") + pos = self._pos + 4 + line = self.text[pos].replace('=',' ') + tokens = line.split()[0::2] + values = line.split()[1::2] + while tokens[0] != "DENSITY": + for i,t in enumerate(tokens): + options[t] = values[i] + pos += 1 + line = self.text[pos].replace('=',' ') + tokens = line.split()[0::2] + values = line.split()[1::2] + line = self.text[pos].replace('=',' ').split() + options[line[2]] = line[3] + pos += 1 + line = self.text[pos].replace('=',' ').split() + options[line[6]] = line[7] + except IndexError: + pass + try: + self.find_string("GUGA DISTINCT ROW TABLE") + pos = self._pos + 3 + line = self.text[pos].replace('=',' ') + tokens = line.split()[0::2] + values = line.split()[1::2] + while tokens[0] != "-CORE-": + for i,t in enumerate(tokens): + options[t] = values[i] + pos += 1 + line = self.text[pos].replace('=',' ') + tokens = line.split()[0::2] + values = line.split()[1::2] + except IndexError: + pass + try: + self.find_string("GUGA DISTINCT ROW TABLE") + self.find_next_string("-CORE-") + pos = self._pos + 1 + line = self.text[pos].replace('=',' ') + tokens = line.split()[0::2] + values = line.split()[1::2] + while tokens[0] != "THE": + for i,t in enumerate(tokens): + options[t] = values[i] + pos += 1 + line = self.text[pos].replace('=',' ') + tokens = line.split()[0::2] + values = line.split()[1::2] + except IndexError: + pass + if options != {}: + self._options = options + return self._options + + def get_spin_restrict(self): + if self._spin_restrict is None: + options = self.options + if options['SCFTYP'] == 'RHF': + self._spin_restrict = True + elif options['SCFTYP'] == 'ROHF': + self._spin_restrict = True + elif options['SCFTYP'] == 'UHF': + self._spin_restrict = False + elif options['CITYP'] != 'NONE' or options['SCFTYP'] == 'MCSCF': + self._spin_restrict = True + return self._spin_restrict + + + def get_conv_threshs(self): + if self._conv_threshs is None: + options = self.options + self._conv_threshs = [] + for m in self.methods: + if m == 'RHF' or m == 'UHF' or m == 'ROHF': + self._conv_threshs.append(float(options['CONV'])) + if m == 'GUGA': + try: + self.find_string("DAVIDSON METHOD") + self.find_next_string("CONVERGENCE CRITERION") + res = float(self.text[self._pos].split('=')[1]) + self._conv_threshs.append(res) + except IndexError: + pass + return self._conv_threshs + + def get_gradient_energy(self): + if self._gradient_energy is None: + try: + self.find_string("GRADIENT OF THE ENERGY") + self.find_next_string("E'X") + pos = self._pos+1 + self._gradient_energy = [] + for i in self.geometry: + buffer = self.text[pos].split() + label = buffer[1] + assert label == i.name + f = map(float,buffer[2:]) + self._gradient_energy.append(f) + pos +=1 + except IndexError: + self._gradient_energy = [] + return self._gradient_energy + + + 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 = [] + poslist = [] + for m in self.methods: + if m == 'RHF' or m == 'ROHF': + try: + self.find_string("FINAL R") + poslist.append(self._pos) + except IndexError: + pass + elif m == 'UHF': + try: + self.find_string("FINAL U") + poslist.append(self._pos) + except IndexError: + pass + elif m == 'GUGA': + try: + self.find_string("END OF DENSITY MATRIX CALCULATION") + poslist.append(self._pos) + except IndexError: + pass + elif m == 'MCSCF': + try: + self.find_string("FINAL MCSCF ENERGY IS") + poslist.append(self._pos) + except IndexError: + pass + else: + try: + self.find_string("FINAL ENERGY IS") + poslist.append(self._pos) + except IndexError: + pass + for pos in poslist: + try: + self._pos = pos + self.find_next_string("ENERGY COMPONENTS") + self.find_next_string("ONE ELECTRON ENERGY") + pos = self._pos + line = self.text[pos].split() + self._one_e_energies.append(float(line[4])) + self.find_next_string("TWO ELECTRON ENERGY") + pos = self._pos + line = self.text[pos].split() + self._two_e_energies.append(float(line[4])) + self.find_next_string("TOTAL ENERGY") + pos = self._pos + line = self.text[pos].split() + self._energies.append(float(line[3])) + self.find_next_string("ELECTRON-ELECTRON POT") + pos = self._pos + line = self.text[pos].split() + self._ee_pot_energies.append(float(line[4])) + self.find_next_string("NUCLEUS-ELECTRON POT") + pos = self._pos + line = self.text[pos].split() + self._Ne_pot_energies.append(float(line[4])) + self.find_next_string("NUCLEUS-NUCLEUS POT") + pos = self._pos + line = self.text[pos].split() + self._pot_energies.append(float(line[4])) + self.find_next_string("TOTAL KINETIC ENE") + pos = self._pos + line = self.text[pos].split() + self._kin_energies.append(float(line[4])) + self.find_next_string("VIRIAL RATIO") + pos = self._pos + line = self.text[pos].split() + self._virials.append(float(line[4])) + except IndexError: + 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_ee_pot_energies(self): + if self._ee_pot_energies is None: + self.get_energies() + return self._ee_pot_energies + + def get_Ne_pot_energies(self): + if self._Ne_pot_energies is None: + self.get_energies() + return self._Ne_pot_energies + + def get_pot_energies(self): + if self._pot_energies is None: + self.get_energies() + return self._pot_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_point_group(self): + if self._point_group is None: + try: + self.find_string("THE POINT GROUP IS") + pos = self._pos + line = self.text[pos].replace(',',' ').split() + group = line[4] + self._point_group = group.replace('N',line[6]) + except IndexError: + pass + return self._point_group + + def get_geometry(self): + if self._geometry is None: + try: + self.find_string(" ATOM ATOMIC COORDINATES") + pos = self._pos + pos += 2 + self._geometry = [] + buffer = self.text[pos].split() + while len(buffer) > 1: + temp = atom() + temp.name = buffer[0] + temp.charge = float(buffer[1]) + temp.coord = (float(buffer[2]), float(buffer[3]), float(buffer[4])) + temp.basis = [] + self._geometry.append(temp) + pos += 1 + buffer = self.text[pos].split() + except IndexError: + pass + 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_dipole(self): + if self._dipole is None: + try: + self.find_string("ELECTROSTATIC MOMENTS") + pos = self._pos + pos += 6 + self._dipole = [] + for d in self.text[pos].split(): + self._dipole.append( float(d) ) + except IndexError: + pass + return self._dipole + + def get_basis(self): + if self._basis is None: + try: + self.find_string("SHELL TYPE") + except IndexError: + return None + pos = self._pos + try: + self.find_next_string("TOTAL NUMBER OF") + except IndexError: + return None + end = self._pos + pos += 4 + basis_read = [] + while pos < end: + line = self.text[pos].split() + bf = [] + while len(line) > 1: + index = int(line[0]) + sym = line[1] + expo = float(line[3]) + coef = float(line[4]) + if sym == "L": + coef2 = float(line[5]) + 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] ) + pos += 1 + + 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]] + 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 0 and pos < end: + try: + il2 = int(line[2])+shift + except ValueError: + line.insert(2,str(curatom)) + self.text[pos] = ' '.join(line) + il2 = int(line[2])+shift + if il2 == 1: + is_zero = False + elif not is_zero and il2 == 0: + il2 -= shift + shift += 100 + il2 += shift + is_zero = True + if curatom < il2: + self._atom_to_ao_range.append(list(ao_range)) + ao_range[1] += 1 + curatom = il2 + ao_range[0] = int(line[0]) + else: + ao_range[1] = int(line[0]) + pos += 1 + line = self.text[pos].split() + self._atom_to_ao_range.append(ao_range) + self._atom_to_ao_range = list(self._atom_to_ao_range[:len(self.geometry)]) + return self._atom_to_ao_range + + + def get_mo_sets(self): + if self._mo_sets is None: + self._mo_sets = [] + self._mo_types = [] + posend = {} + # UHF + if not self.spin_restrict: + try: + index = self.options['SCFTYP']+'a' + self.find_string(" CONVERGED") + self.find_next_string("EIGENVECTORS") + self.find_next_string("1 2") + posend[index] = [] + posend[index].append(self._pos) + self.find_next_string("-- BETA SET --") + posend[index].append(int(self._pos)-1) + self._mo_types.append(index) + index = self.options['SCFTYP']+'b' + self.find_next_string("EIGENVECTORS") + self.find_next_string("1 2") + posend[index] = [] + posend[index].append(self._pos) + self.find_next_string("END OF "+index[:3]+" CALCULATION") + posend[index].append(int(self._pos)) + self._mo_types.append(index) + except IndexError: + pass + # RHF/ROHF + else: + try: + index = self.options['SCFTYP'] + self.find_string(" CONVERGED") + self.find_next_string("EIGENVECTORS") + self.find_next_string("1 2") + posend[index] = [] + posend[index].append(self._pos) + self.find_next_string("END OF ") + posend[index].append(int(self._pos)) + self._mo_types.append(index) + except IndexError: + pass + try: + index = self.options['LOCAL'] + self.find_string("THE BOYS LOCALIZED ORBITALS ARE") + self.find_next_string("1 2") + posend[index] = [] + posend[index].append(self._pos) + self.find_next_string("END OF ") + posend[index].append(int(self._pos)) + self._mo_types.append(index) + except IndexError: + pass + # GUGA + try: + if self.options['CITYP'] == 'NONE' and self.options['SCFTYP'] != 'MCSCF': + raise IndexError + index = self.options['CITYP'] + self.find_string("INITIAL GUESS ORBITALS") + self.find_next_string("1 2") + posend[index] = [] + posend[index].append(self._pos) + self.find_next_string("END OF INITIAL ORBITAL SELECTION") + posend[index].append(self._pos) + self._mo_types.append(index) + except IndexError: + pass + try: + index = 'Natural' + self.find_string("NATURAL ORBITALS IN ATOMIC ORBITAL BASIS") + self.find_next_string("1 2") + posend[index] = [] + posend[index].append(self._pos) + self.find_next_string("END OF DENSITY MATRIX CALCULATION") + posend[index].append(self._pos) + self._mo_types.append(index) + except IndexError: + pass + # MCSCF + try: + if self.options['SCFTYP'] != 'MCSCF': + raise IndexError + index = self.options['SCFTYP'] + self.find_string("MCSCF OPTIMIZED ORBITALS") + self.find_next_string("1 2") + posend[index] = [] + posend[index].append(self._pos) + self.find_next_string("DONE WITH MCSCF ITERATIONS") + posend[index].append(self._pos) + self._mo_types.append(index) + except IndexError: + pass + try: + index = 'Natural' + self.find_string("MCSCF NATURAL ORBITALS") + self.find_next_string("1 2") + posend[index] = [] + posend[index].append(self._pos) + self.find_next_string("LZ VALUE ANALYSIS FOR MOLECULAR ORBITALS") + posend[index].append(self._pos) + self._mo_types.append(index) + except IndexError: + pass + self._mo_sets = {} + for index in self._mo_types: + pos,end = posend[index] + pos -= 1 + curatom = 1 + vectors = [] + while pos < end: + pos += 1 + orbnum = self.text[pos].split() + pos += 1 + line = self.text[pos].split() + if len(line) == 0: + pos += 1 + line = self.text[pos].split() + try: + eigval = [ float(k) for k in line ] + pos += 1 + syms = self.text[pos].split() + pos += 1 + except: + eigval = [] + syms = [] + pass + line = self.text[pos].split() + begin = pos + # Build vectors + pos = begin + lmax = len(orbnum)+4 + for l in range(4,lmax): + pos = begin + line = self.text[pos] + line = line.replace('-',' -').split() + v = orbital() + v.set = index + v.basis = self.basis + if len(eigval)>0: v.eigenvalue = eigval[l-4] + if len(syms) >0: v.sym = syms[l-4] + while len(line) > 0 and pos < end: + try: + bid = int(line[2]) + except ValueError: + line.insert(2,str(curatom)) + self.text[pos] = ' '.join(line) + v.vector.append(float(line[l])) + pos += 1 + line = self.text[pos] + line = line.replace('-',' -').split() + vectors.append(v) + self._mo_sets[index] = vectors + if 'BOYS' in self._mo_sets: + self._mo_sets['BOYS'] += \ + self._mo_sets['RHF'][len(self._mo_sets['BOYS']):] + + return self._mo_sets + + def get_mulliken_mo(self): + if self._mulliken_mo is None: + posend = [] + vectors = [] + if self.spin_restrict: + try: + self.find_string("MULLIKEN ATOMIC POPULATION IN EACH MOLECULAR ORBITAL") + except IndexError: + return None + pos = self._pos + try: + self.find_string("ATOMIC SPIN POPULATION") + except IndexError: + try: + self.find_next_string("POPULATIONS IN EACH AO") + except IndexError: + return None + end = self._pos + pos += 2 + posend.append([pos,end]) + else: + try: + self.find_string("MULLIKEN ATOMIC POPULATION IN EACH MOLECULAR ORBITAL") + self.find_next_string("ALPHA ORBITALS") + pos = self._pos + self.find_next_string("MULLIKEN ATOMIC POPULATION") + end = self._pos + pos += 2 + posend.append([pos,end]) + self.find_next_string("MULLIKEN ATOMIC POPULATION") + self.find_next_string("BETA ORBITALS") + pos = self._pos + self.find_next_string("ATOMIC SPIN POP") + end = self._pos + pos += 2 + posend.append([pos,end]) + except IndexError: + return None + for pos, end in posend: + while pos < end: + try: + pos += 2 + line = self.text[pos].split() + eigval = [ float(k) for k in line ] + pos += 2 + line = self.text[pos].split() + begin = pos + lmax = len(eigval)+1 + for l in range(1,lmax): + pos = begin + line = self.text[pos].split() + v = orbital() + v.set = "Mulliken_MO" + v.basis = self.basis + v.eigenvalue = eigval[l-1] + while len(line) > 1 and pos < end: + v.vector.append(float(line[l])) + pos += 1 + line = self.text[pos].split() + vectors.append(v) + pos += 1 + except: + return None + self._mulliken_mo = vectors + return self._mulliken_mo + + def get_mulliken_ao(self): + if self._mulliken_ao is None: + try: + self.find_string(" POPULATIONS IN EACH AO") + except IndexError: + return None + self._pos += 2 + pos = self._pos + self.find_next_string("----") + end = self._pos - 1 + line = self.text[pos][26:].split() + vm = orbital() + vl = orbital() + vm.set = "Mulliken_AO" + vl.set = "Lowdin_AO" +# vm.basis = self.basis +# vl.basis = self.basis + vm.eigenvalue = 0. + vl.eigenvalue = 0. + while pos < end: + value = float(line[0]) + vm.vector.append(value) + vm.eigenvalue += value + value = float(line[1]) + vl.vector.append(value) + vl.eigenvalue += value + pos += 1 + line = self.text[pos][26:].split() + self._mulliken_ao = vm + self._lowdin_ao = vl + return self._mulliken_ao + + def get_lowdin_ao(self): + if self._lowdin_ao is None: + self.get_mulliken_ao() + return self._lowdin_ao + + + def get_mulliken_atom(self): + if self._mulliken_atom is None: + try: + self.find_string(" TOTAL MULLIKEN AND LOWDIN ATOMIC POPULATIONS") + except IndexError: + return None + self._pos += 2 + pos = self._pos + self.find_next_string("-----------") + end = self._pos - 1 + line = self.text[pos].split() + vm = orbital() + vl = orbital() + vm.set = "Mulliken_Atom" + vl.set = "Lowdin_Atom" +# vm.basis = self.geometry +# vl.basis = self.geometry + vm.eigenvalue = 0. + vl.eigenvalue = 0. + while pos < end: + value = float(line[2]) + vm.vector.append(value) + vm.eigenvalue += value + value = float(line[4]) + vl.vector.append(value) + vl.eigenvalue += value + pos += 1 + line = self.text[pos].split() + self._mulliken_atom = vm + self._lowdin_atom = vl + return self._mulliken_atom + + def get_lowdin_atom(self): + if self._lowdin_atom is None: + self.get_mulliken_atom() + return self._lowdin_atom + + def get_two_e_int_ao(self): + if self._two_e_int_ao is None: + try: + self.find_string("2 ELECTRON INTEGRALS") + except IndexError: + return None + self.find_next_string("BYTES/INTEGRAL") + pos = self._pos + line = self.text[pos].split() + Nint = int( line[1] ) + Nbytes = int( line[6] ) + assert ( Nbytes == 12 or Nbytes == 16 ) + try: + file = gamessFile.two_e_int_array(self.two_e_int_filename) + except IOError: + return None + if Nbytes == 12: + tc = 1 + tci = 'b' + elif Nbytes == 16: + tc = 2 + tci = 'h' + file.form = 'l'+4*Nint*tci+Nint*'d' + self._two_e_int_ao = file + #TODO + return self._two_e_int_ao + + def get_num_alpha(self): + if self._num_alpha is None: + self._num_alpha = (self.num_elec + 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 - self.multiplicity+1)/2 + return self._num_beta + + def get_determinants_mo_type(self): + if self._determinants_mo_type is None: + # Mono-determinant case + if self.options['SCFTYP'] == 'RHF' or\ + self.options['SCFTYP'] == 'UHF' or\ + self.options['SCFTYP'] == 'ROHF': + self._determinants_mo_type = self.mo_types[-1] + # Multi-determinant case TODO + elif self.options['SCFTYP'] == 'MCSCF': + self._determinants_mo_type = 'MCSCF' + # Multi-determinant case TODO + elif self.options['CITYP'] != '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.determinant_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 = determinants + return self._determinants + + def get_csf(self): + if self._csf is None: + + # Mono-determinant case + csf = [] + if self.options['SCFTYP'] == 'RHF' or\ + self.options['SCFTYP'] == 'UHF' or\ + self.options['SCFTYP'] == 'ROHF': + new_csf = CSF() + new_spin_a = [] + new_spin_b = [] + for i in range(self.num_alpha): + new_spin_a.append(i) + for i in range(self.num_beta): + new_spin_b.append(i) + new_csf.append(1.,new_spin_a,new_spin_b) + csf.append(new_csf) + self._determinants_mo_type = self.mo_types[-1] + + # Multi-determinant case + if self.options['CITYP'] == 'GUGA' or self.options['SCFTYP'] == 'MCSCF': + try: + self.find_string('SYMMETRIES FOR THE') + pos = self._pos + buffer = self.text[pos].split() + ncore = int (buffer[3]) + nact = int (buffer[5]) + self.find_string('TOTAL NUMBER OF INTEGRALS') + end = self._pos + self.find_string('DETERMINANT CONTRIBUTION TO CSF') + self.find_next_string('CASE VECTOR') + pos = self.pos + while pos < end: + pos += 4 + if len (self.text[pos].split()) > 0: + this_csf = CSF() + while not self.text[pos].startswith(' CASE VECTOR') and len(self.text[pos])>1: + buffer = self.text[pos].split('=')[1].split(':') + mostr = buffer[1].replace("-"," -").split() + mo = [] + for i in mostr: + mo.append (int(i)) + p_count=0 + for j in range ( len(mo) ): + if mo[j] < 0: + for l in range ( j+1, len(mo) ): + if mo[l] > 0: + p_count+=1 + coef = float (buffer[0])*(-1.0)**p_count + tempcsf_a = [] + tempcsf_b = [] + for i in range(ncore): + tempcsf_a.append(i) + tempcsf_b.append(i) + for i in range(len(mo)): + orb = int (mo [i]) + if orb > 0: + tempcsf_a.append(orb-1) + else: + tempcsf_b.append(-orb-1) + this_csf.append(coef,tempcsf_a,tempcsf_b) + pos+=1 + csf.append(this_csf) + except IndexError: + pass + + 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: + # Mono-determinant case + if self.options['SCFTYP'] == 'RHF' or\ + self.options['SCFTYP'] == 'UHF' or\ + self.options['SCFTYP'] == 'ROHF': + self._det_coefficients = [ [1.] ] + # Multi-determinant case + 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_csf_coefficients(self): + if self._csf_coefficients is None: + # Mono-determinant case + if self.options['SCFTYP'] == 'RHF' or\ + self.options['SCFTYP'] == 'UHF' or\ + self.options['SCFTYP'] == 'ROHF': + self._csf_coefficients = [ [1.] ] + # Multi-determinant case + elif self.options['CITYP'] == 'GUGA' or self.options['SCFTYP'] == 'MCSCF': + self._csf_coefficients = [ [] for i in range(self.num_states) ] + self._pos = 0 + for state in range(self.num_states): + cur_csf=0 + self.find_next_string("STATE #") + self._pos += 4 + pos = self._pos + try: + while True: + line = self.text[pos].split() + csf=int(line[0]) + for c in range(cur_csf,csf-1): + self._csf_coefficients[state].append(0.) + cur_csf=csf + self._csf_coefficients[state].append(float(line[1])) + pos += 1 + except: + pass + + return self._csf_coefficients + + def get_num_states(self): + if self._num_states is None: + self._num_states=1 + try: + self.find_string("NUMBER OF STATES REQUESTED") + pos = self._pos + buffer = self.text[pos].split('=') + self._num_states = int(buffer[1]) + except IndexError: + pass + 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 + + atom_to_ao_range = property(fget=get_atom_to_ao_range) + + for i in "two_e_int_mo_filename two_e_int_ao_filename".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 two_e_int_array(object): + + def __init__(self,filename): + self._file = fortranBinary(filename,"rb") + + def get_form(self): + return self._file.form + + def set_form(self,form): + self._file.form = form + + def __iter__(self): + for rec in self._file: + Nint = rec.pop(0) + rec_i = rec[:4*Nint] + rec_f = rec[4*Nint:] + del rec + for m in range(Nint): + n = 4*m + i = rec_i[n] + j = rec_i[n+1] + k = rec_i[n+2] + l = rec_i[n+3] + v = rec_f[m] + int = integral() + int.indices = (i, j, k, l) + int.value = v + yield int + + def __getitem__(self,index): + self._file.seek(0) + rec = self._file[0] + Nint = rec.pop(0) + while (index >= Nint): + index -= Nint + rec = self._file.next() + Nint = rec.pop(0) + + rec_i = rec[:4*Nint] + rec_f = rec[4*Nint:] + del rec + m = index + n = 4*m + i = rec_i[n] + j = rec_i[n+1] + k = rec_i[n+2] + l = rec_i[n+3] + v = rec_f[m] + int = integral() + int.indices = (i, j, k, l) + int.value = v + return int + + form = property ( get_form, set_form ) + +# Output Routines +# --------------- + +def gamess_write_contrl(res,file, runtyp="ENERGY",guess="HUCKEL", + exetyp="RUN",scftyp=None,moread=False,NoSym=False,units="BOHR"): + print >>file, " $CONTRL" + print >>file, " EXETYP=",exetyp + print >>file, " COORD= UNIQUE", " UNITS=",units + print >>file, " RUNTYP=", runtyp + if NoSym: + print >>file, " NOSYM=1" + if scftyp is None: + if res.num_alpha == res.num_beta: + scftyp="RHF" + else: + scftyp="ROHF" + print >>file, " SCFTYP=", scftyp + print >>file, " MULT=", res.multiplicity + print >>file, " ICHARG=",str(int(res.charge)) + print >>file, " $END" + print >>file, "" + print >>file, " $GUESS" + if moread: + print >>file, " GUESS=MOREAD" + print >>file, " NORB=",len(res.mo_sets[res.mo_types[-1]]) + else: + print >>file, " GUESS=HUCKEL" + print >>file, " $END" + + +def gamess_write_vec(res, file): +# if motype is None: +# motypelist = res.mo_types +# else: +# motypelist = [motype] +# for motype in motypelist: + motype = res.mo_types[-1] + print >>file, motype, "MOs" + print >>file, " $VEC" + moset = res.mo_sets[motype] + moindex = res.closed_mos + res.active_mos + \ + res.virtual_mos + for idx,moi in enumerate(moindex): + mo = moset[moi] + v = list(mo.vector) + line = 1 + while (len(v) > 0): + fields = [] + monum = "%4d"%(idx+1 ) + monum = monum[-2:] + fields.append ( monum ) + linum = "%4d"%(line) + linum = linum[-3:] + fields.append ( linum ) + for n in range(5): + try: + fields.append ( "%15.8E"%v.pop(0) ) + except IndexError: + pass + print >>file, ''.join(fields) + line += 1 + print >>file, " $END" + +def gamess_write_data(res, file): + print >>file, " $DATA" + print >>file, res.title + try: + point_group = res.point_group + N = point_group[1] + if N != '1': + point_group = point_group[0]+'N'+point_group[2:].capitalize() + point_group += ' '+N+'\n' + except TypeError: + point_group = 'C1' + N = '1' + print >>file, point_group + for at in res.geometry: + print >>file, "%10s %3.1f %15.8f %15.8f %15.8f" % tuple( \ + [at.name.ljust(10), at.charge]+list(at.coord) ) + for contr in at.basis: + sym = contr.sym + if 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 == 'xxxxx' : + sym = 'H' + doPrint = True + elif sym == 'xxxxxx' : + sym = 'I' + doPrint = True + else: + doPrint = False + if doPrint: + print >>file, "%4s%8d"%(sym, len(contr.prim)) + for i, p in enumerate(contr.prim): + if not isinstance(p,gaussian): + raise TypeError("Basis function is not a gaussian") + print >>file, "%6d%26.10f%12.8f"%(i+1,p.expo,contr.coef[i]) + print >>file, "" + print >>file, " $END" + +def gamess_write_integrals(res,file): + filename = file.name.replace('.inp','.moints') + +def write_gamess(file,stdout): + gamess_write_contrl(file,stdout) + x = file.basis + x = file.geometry + gamess_write_data(file,stdout) + gamess_write_vec(file,stdout) + + +fileTypes.append(gamessFile) + +if __name__ == '__main__': + main(gamessFile) diff --git a/resultsFile/Modules/gaussianFile.py b/resultsFile/Modules/gaussianFile.py new file mode 100755 index 0000000..b81411f --- /dev/null +++ b/resultsFile/Modules/gaussianFile.py @@ -0,0 +1,1040 @@ +#!/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 + +gaussianFile_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 gaussianFile(resultsFile): + """ Class defining the gaussian file. + """ + + local_vars = list(local_vars) + defined_vars = list(gaussianFile_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_date(self): + if self._date is None: + self._date = self.options[0][8] + return self._date + + 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_dipole(self): + if self._dipole is None: + self._dipole = [] + dip = self.text[pos] + for x in dip: + self._dipole.append(float(x)) + return self._dipole + + def get_multiplicity(self): + if self._multiplicity is None: + self._multiplicity = int(self.options[3][0].split(',')[1]) + return self._multiplicity + + 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_nuclear_energy(self): + if self._nuclear_energy is None: + self.find_string("N-N") + pos = self._pos + self._nuclear_energy = self.text[pos].replace('=','= ').split()[1] + self._nuclear_energy = float(self._nuclear_energy.replace('D','E')) + return self._nuclear_energy + + def get_title(self): + if self._title is None: + self._title = self.options[2][0].strip() + return self._title + + 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_cpu_time(self): + if self._cpu_time is None: + try: + self.find_string("Job cpu time") + except IndexError: + return None + pos = self._pos + line = self.text[pos].split() + self._cpu_time = float(line[3])*60*60*24 + self._cpu_time += float(line[5])*60*60 + self._cpu_time += float(line[7])*60 + self._cpu_time += float(line[9]) + self._cpu_time = str(self._cpu_time)+' s' + return self._cpu_time + + 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 + 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.startswith("U"): 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_energies(self): + if self._energies is None: + self._energies = [] + for m in self.methods: + if m == 'CASSCF' \ + or m.startswith('R')\ + or m.startswith('U'): + #if self.point_group == "C01": + self._energies.append(float(self.options[4][2].split('=')[1])) + #else: + # self._energies.append(float(self.options[4][3].split('=')[1])) + return self._energies + + 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_pot_energies(self): + if self._pot_energies is None: + self._pot_energies = [] + for i,e in enumerate(self.Ne_pot_energies): + self._pot_energies.append (self.ee_pot_energies[i]+\ + e+self.nuclear_energy) + return self._pot_energies + + def get_kin_energies(self): + if self._kin_energies is None: + self.find_string("N-N") + pos = self._pos + self._kin_energies = [float(self.text[pos].replace('=','= ').split()[5])] + return self._kin_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("X Y Z") + self._pos += 1 + except IndexError: + pass + pos +=1 + self._pos=pos + self.find_next_string("-----") + end = self._pos + while pos 0: + basis_read.append( [index,sym,bf,iatom] ) + if line[0].startswith('==============='): + pos = end + 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 k0: v.sym = syms[l] + while len(line) > 0 and pos < end: + line = self.text[pos][21:].strip().replace('-',' -') + line = line.split() + try: + v.vector.append(float(line[l])) + except: + break + pos += 1 + line = self.text[pos][1:21].strip().split() + if "DENSITY" in line: + line = [] + vectors.append(v) + self._mo_sets[index] = vectors + if self._mo_sets == {}: + self._mo_sets = None + return self._mo_sets + + 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_num_alpha(self): + if self._num_alpha is None: + self.find_string("alpha electrons") + pos = self._pos + self._num_alpha = int(self.text[pos].split()[0]) + return self._num_alpha + + def get_num_beta(self): + if self._num_beta is None: + self.find_string("beta electrons") + pos = self._pos + self._num_beta = int(self.text[pos].split()[3]) + 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): + method = self.methods[0] + if self._csf is None: + + csf = [] + # Mono-determinant case + if method.startswith("R") and not method.startswith("RO"): + new_csf = CSF() + new_spin_det = [] + for i in range(self.num_alpha): + new_spin_det.append(i) + new_csf.append(1.,new_spin_det,new_spin_det) + csf.append(new_csf) + self._determinants_mo_type = self.mo_types[-1] + + elif method.startswith("RO") or method.startswith("U"): + new_csf = CSF() + new_spin_deta = [] + new_spin_detb = [] + for i in range(self.num_alpha): + new_spin_deta.append(i) + for i in range(self.num_beta): + new_spin_detb.append(i) + new_csf.append(1.,new_spin_deta,new_spin_detb) + csf.append(new_csf) + self._determinants_mo_type = self.mo_types[-1] + + # Multi-determinant case + if method == 'CASSCF': + try: + self.find_string('Enter MCSCF program') + self.find_next_string('CORE-ORBITALS') + pos = self._pos + buffer = self.text[pos].split('=') + ncore = int(buffer[-1]) + self.find_string('no. active ELECTRONS') + pos = self._pos + buffer = self.text[pos].split() + nact = int (buffer[4]) + self.find_next_string('PRIMARY BASIS FUNCTION') + pos = self._pos + mostr = self.text[pos].split('=')[1].split() + tempcsf_a = [] + tempcsf_b = [] + for i in range(ncore): + tempcsf_a.append(i) + tempcsf_b.append(i) + to_add = tempcsf_a + old=0 + for i in mostr: + i = int(i) + if i <= old: to_add=tempcsf_b + to_add.append(i+ncore-1) + old=i + this_csf = CSF() + this_csf.append(1.,tempcsf_a,tempcsf_b) + csf.append(this_csf) + pos = self._pos+1 + self.find_next_string("NO OF BASIS FUNCTIONS") + end = self._pos + while pos < end: + this_csf = CSF() + tempcsf_a = [] + tempcsf_b = [] + for tempcsf in [tempcsf_a,tempcsf_b]: + pos += 1 + buffer = self.text[pos].split() + for i in range(ncore): + tempcsf.append(i) + for i in buffer: + orb = int (i)+ncore-1 + tempcsf.append(orb) + this_csf.append(1.,tempcsf_a,tempcsf_b) + pos += 1 + csf.append(this_csf) + except IndexError: + pass + + 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: + # Mono-determinant case + if self.options[0][4].startswith('R') or\ + self.options[0][4].startswith('U'): + self._det_coefficients = [ [1.] ] + # Multi-determinant case TODO + 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_csf_coefficients(self): + if self._csf_coefficients is None: + # Mono-determinant case + if self.options[0][4].startswith('R') or\ + self.options[0][4].startswith('U'): + self._csf_coefficients=[ [1.] ] + # Multi-determinant case TODO + elif self.options[0][4] == 'CASSCF': + self._csf_coefficients = [ [] for i in range(self.num_states) ] + self._pos = 0 + regexp = re.compile("\( *([0-9]+)\)([ |-][0-9]*\.[0-9]*) *") + for state in range(self.num_states): + for i in self.csf: + self._csf_coefficients[state].append(0.) + cur_csf=0 + self.find_string("2ND ORD PT ENERGY") + self._pos += 1 + self.find_next_string("EIGENVALUE ") + pos = self._pos + 1 + matches = [1] + while matches != []: + matches = regexp.findall(self.text[pos]) + for m in matches: + self._csf_coefficients[state][int(m[0])-1] = float(m[1]) + pos += 1 + return self._csf_coefficients + + def get_occ_num(self): + if self._occ_num is None: + occ = {} + method = self.methods[0] + for motype in self.mo_types: + occ[motype] = [ 0. for mo in self.mo_sets[motype] ] + if method.startswith("R") and not method.startswith("RO"): + for i in range(self.num_alpha): + occ[motype][i] = 2. + elif method.startswith("RO"): + for i in range(self.num_beta): + occ[motype][i] = 2. + for i in range(self.num_beta,self.num_alpha): + occ[motype][i] = 1. + elif self.mulliken_mo is not None: + for motype in self.mo_types: + occ[motype] = [ mo.eigenvalue for mo in self.mulliken_mo ] + self._occ_num = occ + return self._occ_num + + def get_num_states(self): + if self._num_states is None: + self._num_states=1 + try: + self.find_string("NUMBER OF STATES REQUESTED") + pos = self._pos + buffer = self.text[pos].split('=') + self._num_states = int(buffer[1]) + except IndexError: + pass + 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 + + + class two_e_int_array(object): + + def __init__(self,filename): + self._file = fortranBinary(filename,"rb") + + def get_form(self): + return self._file.form + + def set_form(self,form): + self._file.form = form + + def __iter__(self): + for rec in self._file: + Nint = rec.pop(0) + rec_i = rec[:4*Nint] + rec_f = rec[4*Nint:] + del rec + for m in range(Nint): + n = 4*m + i = rec_i[n] + j = rec_i[n+1] + k = rec_i[n+2] + l = rec_i[n+3] + v = rec_f[m] + int = integral() + int.indices = (i, j, k, l) + int.value = v + yield int + + def __getitem__(self,index): + self._file.seek(0) + rec = self._file[0] + Nint = rec.pop(0) + while (index >= Nint): + index -= Nint + rec = self._file.next() + Nint = rec.pop(0) + + rec_i = rec[:4*Nint] + rec_f = rec[4*Nint:] + del rec + m = index + n = 4*m + i = rec_i[n] + j = rec_i[n+1] + k = rec_i[n+2] + l = rec_i[n+3] + v = rec_f[m] + int = integral() + int.indices = (i, j, k, l) + int.value = v + return int + + form = property ( get_form, set_form ) + +# Output Routines +# --------------- + +def write_vec(res, file): + for motype in res.mo_types: + print >>file, " $VEC" + moset = res.mo_sets[motype] + for idx,mo in enumerate(moset): + v = list(mo.vector) + line = 1 + while (len(v) > 0): + fields = [] + monum = "%4d"%(idx ) + monum = monum[-2:] + fields.append ( monum ) + linum = "%4d"%(line) + linum = linum[-3:] + fields.append ( linum ) + for n in range(5): + fields.append ( "%15.8E"%v.pop(0) ) + print >>file, ''.join(fields) + line += 1 + print >>file, " $END" + +def write_data(res, file): + print >>file, " $DATA" + print >>file, res.title + point_group = res.point_group + N = point_group[1] + if N != '1': + point_group = point_group[0]+'N'+point_group[2:] + point_group += ' '+N + print >>file, point_group + if N != '1': + print >>file, "" + for at in res.geometry: + print >>file, "%10s %3.1f %15.8f %15.8f %15.8f" % tuple( \ + [at.name.ljust(10), at.charge]+list(at.coord) ) + for contr in at.basis: + sym = contr.sym + if sym == 's' : + doPrint = True + elif sym == 'x' : + sym = 'P' + doPrint = True + elif sym == 'xx' : + sym = 'D' + doPrint = True + elif sym == 'xxx' : + sym = 'F' + doPrint = True + elif sym == 'xxxx' : + sym = 'G' + doPrint = True + elif sym == 'xxxxx' : + sym = 'H' + doPrint = True + else: + doPrint = False + if doPrint: + print >>file, "%4s%8d"%(sym, len(contr.prim)) + for i, p in enumerate(contr.prim): + if not isinstance(p,gaussian): + raise TypeError("Basis function is not a gaussian") + print >>file, "%6d%26.10f%12.8f"%(i+1,p.expo,contr.coef[i]) + print >>file, "" + print >>file, " $END" + + + +fileTypes.append(gaussianFile) + +if __name__ == '__main__': + main(gaussianFile) diff --git a/resultsFile/Modules/include.py b/resultsFile/Modules/include.py new file mode 100755 index 0000000..4e1d786 --- /dev/null +++ b/resultsFile/Modules/include.py @@ -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') + diff --git a/resultsFile/Modules/molproFile.py b/resultsFile/Modules/molproFile.py new file mode 100755 index 0000000..25cb67f --- /dev/null +++ b/resultsFile/Modules/molproFile.py @@ -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 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 pos1) 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) + diff --git a/resultsFile/Modules/qmcchemFile.py b/resultsFile/Modules/qmcchemFile.py new file mode 100755 index 0000000..5c1f552 --- /dev/null +++ b/resultsFile/Modules/qmcchemFile.py @@ -0,0 +1,1357 @@ +#!/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 sys + +qmcchemFile_defined_vars = [] + +class qmcchemFile(resultsFile): + """ Class defining the QMC=Chem file. + """ + + local_vars = list(local_vars) + defined_vars = list(qmcchemFile_defined_vars) + + 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 + + + +JastLink='Dets' +def qmcchem_write_qmcin(resultsFile,method='VMC',irestart=0,version=2007.1,\ + JastType='None',NumWalkers=100,NumSteps=100,NumCalcs=50,\ + TimeStep=0.001,PointsRead=False,PointsWrite=True,DriftInitMax=0.,\ + NumIter=1000, Minimize='Variance', EtaMin=0.9,\ + Theta=0, NumBlocks=100, LenBlocks=100, srmcAlgo='Sorella',\ + NumProj=100,NumStepCorr=2,Jump=1,NumTime=10,state=0,\ + Schema='Langevin',JastLink=JastLink,Adiag=0.001): + + try: + EnerRef = resultsFile.energies[-1] + except TypeError: + EnerRef = 0. + if DriftInitMax == 0.: + DriftInitMax = max(10.,abs(EnerRef)/4.) + + file = open("qmcin","w") + print >>file, "&general" + print >>file, " Irestart=",irestart + print >>file, " Version =%5.1f"%(version) + print >>file, " Method ='"+method+"'" + print >>file, "/\n" + print >>file, "=geometry_data=" + geometry = resultsFile.geometry + for at in geometry: + if resultsFile.units == 'BOHR': + print >>file, "0 %5d %16f %16f %16f 0 4"%tuple([int(at.charge)]+list(at.coord)) + else: + print >>file, "0 %5d %16f %16f %16f 0 4"%tuple([int(at.charge)]+ \ + [i/a0 for i in at.coord]) + print >>file, "=end=\n" + print >>file, "&wfunction" + print >>file, " NumAlpha =", resultsFile.num_alpha + print >>file, " NumBeta =", resultsFile.num_beta + print >>file, " JastType ='"+JastType+"'" +# print >>file, " JastLink ='"+JastLink+"'" +# print >>file, " NormMOs =", 0 + print >>file, "/\n" + print >>file, "&vmc" + print >>file, " Schema = '"+Schema+"'" + print >>file, " NumWalkers =", NumWalkers + print >>file, " LenBlocks =", NumSteps + print >>file, " NumCalcs =", NumCalcs + if Schema == 'Langevin': + TimeStep2 = 0.2 + else: + TimeStep2 = TimeStep + print >>file, " TimeStep =", TimeStep2 + print >>file, " PointsRead = ."+str(PointsRead)+"." + print >>file, " PointsWrite = ."+str(PointsWrite)+"." + print >>file, " EnerRef =", EnerRef + print >>file, " EPLF = .False." + print >>file, " DriftInitMax =", DriftInitMax + print >>file, "/\n" + print >>file, "&optimlinear" + print >>file, " Schema = '"+Schema+"'" + print >>file, " NumWalkers =", NumWalkers + print >>file, " LenBlocks =", LenBlocks + print >>file, " NumCalcs =", NumCalcs + print >>file, " TimeStep =", TimeStep2 + print >>file, " Adiag =", Adiag + print >>file, "/\n" + print >>file, "&dmc" + print >>file, " NumWalkers =", NumWalkers + print >>file, " LenBlocks =", LenBlocks + print >>file, " NumCalcs =", NumCalcs + print >>file, " TimeStep =", TimeStep + print >>file, " PointsRead = ."+str(PointsRead)+"." + print >>file, " PointsWrite = ."+str(PointsWrite)+"." + print >>file, " TimeStep =", TimeStep + print >>file, " EnerRef =", EnerRef + print >>file, " EPLF = .False." + print >>file, "/\n" +# print >>file, "&srmc" +# print >>file, " Algorithm ='"+srmcAlgo+"'" +# print >>file, " NumBlocks =", NumBlocks +# print >>file, " NumWalkers =", NumWalkers +# print >>file, " LenBlocks =", LenBlocks +# print >>file, " NumProj =", NumProj +# print >>file, " TimeStep =", TimeStep +# print >>file, " PointsRead = ."+str(PointsRead)+"." +# print >>file, " PointsWrite = ."+str(PointsWrite)+"." +# print >>file, " EnerRef =", resultsFile.energies[-1] +# print >>file, " DriftInitMax =", DriftInitMax +# print >>file, "/\n" +# print >>file, "&pdmc" +# print >>file, " NumWalkers =", NumWalkers +# print >>file, " LenBlocks =", NumSteps +# print >>file, " TimeStep =", TimeStep +# print >>file, " PointsRead = ."+str(PointsRead)+"." +# print >>file, " PointsWrite = ."+str(PointsWrite)+"." +# print >>file, " EnerRef =", resultsFile.energies[-1] +# print >>file, " DriftInitMax =", DriftInitMax +# print >>file, " NumStepCorr =", NumStepCorr +# print >>file, " Jump =", Jump +# print >>file, "/\n" +# print >>file, "&gfmc" +# print >>file, " NumBlocks =", NumBlocks +# print >>file, " NumWalkers =", NumWalkers +# print >>file, " LenBlocks =", LenBlocks +# print >>file, " NumCalcs =", NumCalcs +# print >>file, " TimeStep =", TimeStep +# print >>file, " PointsRead = ."+str(PointsRead)+"." +# print >>file, " PointsWrite = ."+str(PointsWrite)+"." +# print >>file, " EnerRef =", resultsFile.energies[-1] +# print >>file, " DriftInitMax =", DriftInitMax +# print >>file, " Theta =", Theta +# print >>file, " NumTime =", NumTime +# print >>file, "/\n" +# print >>file, "&plotorb" +# print >>file, " Idim = 1" +# print >>file, " Norb = 1" +# print >>file, " MO =", resultsFile.num_elec/2 +# print >>file, " Npoints = 20" +# print >>file, " x0=0., y0=0., z0=-5." +# print >>file, " x1=0., y1=0., z1= 5." +# print >>file, "/\n" +# print >>file, "&test" +# print >>file, " NumWalkers =", NumWalkers +# print >>file, " TimeStep =", TimeStep +# print >>file, " PointsRead = ."+str(PointsRead)+"." +# print >>file, " PointsWrite = ."+str(PointsWrite)+"." +# print >>file, " Theta =", Theta +# print >>file, "/\n" + +def qmcchem_write_qmc_det(resultsFile,MO_type=None,state=0,det_thr=0., \ + Slater=False,FitCusp=True): + if MO_type is None: + MO_type = resultsFile.determinants_mo_type + try: + MO_set = resultsFile.mo_sets[MO_type] + except KeyError: + print "MO type "+MO_type+" not present in "+resultsFile.filename + return + file = open('qmc.det','w') + print >>file, "=jastrow_per_determinant=" + print >>file, " (ALL,1)" + print >>file, "=end=\n" + print >>file, "=mos_det_data=" + MOindices = resultsFile.closed_mos + resultsFile.active_mos + NumMOs = len(MOindices) + MOs = [] + for i in MOindices: + MOs.append(MO_set[i]) + if Slater: + nslater = 0 + eigthr = [ MOs[i].eigenvalue for i in range(NumMOs) ] + eigthr.sort() + eigthr = eigthr[len(resultsFile.geometry)-1] + + for i in range(NumMOs): + if Slater and MOs[i].eigenvalue <= eigthr: + print >>file, i+1, 'Slater '+str(nslater+1)+' '+ \ + resultsFile.geometry[nslater].name[0]+'1s' + nslater += 1 + else: + if FitCusp: + print >>file, i+1, 'FitCusp 0.2' + else: + print >>file, i+1, 'Gaussian' + print >>file, "=end=\n" + print >>file, "=psi_data=" + NumAlpha = resultsFile.num_alpha + NumBeta = resultsFile.num_beta + for orb0 in resultsFile.closed_mos: +# print >>file, '('+str(MOindices.index(orb0)+1)+',0) ', + print >>file, str(MOindices.index(orb0)+1)+' ', + print >>file, '\n', + MyList = [] + for i,j in zip(resultsFile.det_coefficients[state],\ + resultsFile.determinants): + MyList.append((-abs(i),i,j)) + MyList.sort() + if det_thr >= 1.: + det_thr = MyList[0][1] +# if det_thr > abs(MyList[0][1]): +# det_thr = abs(MyList[0][1]) + for bid, detcoef, det in MyList: + if abs(detcoef) >= det_thr: + print >>file, "0", detcoef + for orb0 in det['alpha']: + if orb0 in resultsFile.active_mos: +# print >>file, '('+str(MOindices.index(orb0)+1)+',0) ', + print >>file, str(MOindices.index(orb0)+1)+' ', + print >>file, '\n', + for orb0 in det['beta']: + if orb0 in resultsFile.active_mos: +# print >>file, '('+str(MOindices.index(orb0)+1)+',0) ', + print >>file, str(MOindices.index(orb0)+1)+' ', + print >>file, '\n', + print >>file, "=end=\n" + file.close() + +def qmcchem_write_qmc_mos(resultsFile,MO_type=None,thresh=0.): + if MO_type is None: + MO_type = resultsFile.determinants_mo_type + try: + allMOs = resultsFile.uncontracted_mo_sets[MO_type] + except KeyError: + print "MO type not present in "+resultsFile.filename + return + file = open('qmc.mos', 'w') + print >>file, "=mos_data=" + MOindices = resultsFile.closed_mos + resultsFile.active_mos+ resultsFile.virtual_mos + NumMOs = len(MOindices) + MOs = [] + for i in MOindices: + MOs.append(allMOs[i]) + for i in range(NumMOs): + to_print = [] + for idx,j in enumerate(MOs[i].vector): + if abs(j) > thresh: + buffer = "0 %d %20.10f\n"%(idx+1,j) + to_print.append(buffer) + buffer = " %d %d ! %f\n"%(i+1,len(to_print),MOs[i].eigenvalue) + to_print.insert(0,buffer) + file.writelines(to_print) + print >>file, "=end=" + file.close() + +def sort_xyz(s): + l = list(s) + l.sort() + result = "" + for i in l: + result += i + return result + +def qmcchem_write_qmc_basis(resultsFile): + basis = resultsFile.uncontracted_basis + geom = resultsFile.geometry + file = open('qmc.basis', 'w') + print >>file, "=ao_basis=" + 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 = prim.sym + print >>file, "%3d %s %f"%(atom_id, sym, prim.expo) + print >>file, "=end=" + file.close() + +def qmcchem_write_qmc_aos(): + file = open('qmc.aos', 'w') + print >>file, "=ao_slater=" + print >>file, "=end=" + file.close() + +def qmcchem_write_qmc_jast_det(resultsFile): + + file = open('qmc.jast_det','w') + print >>file, "=jast_det_data=" + print >>file, 1 + nparm = (4+len(resultsFile.geometry)) + for i in xrange (nparm): + print >>file, 0, 0. + print >>file, "=end=" + file.close() + + file = open('qmc.jast_det.general','w') + print >>file, "=jast_det_data=" + print >>file, 1 + nparm = (19+9*len(resultsFile.geometry)) + for i in xrange (nparm): + print >>file, 0, 0. + print >>file, "=end=" + file.close() + +def qmcchem_write_qmc_seed(): + file = open("qmc.seed","w") + print >>file, """=random_seed= +1 + 1 123456789 123223945 + 2 744725 23082440 + 3 96136469 236392769 + 4 141196565 165692149 + 5 135925013 158284772 + 6 80321813 164064526 + 7 242822421 3950450 + 8 86555925 175193106 + 9 148393237 240761836 + 10 159898901 171522049 + 11 121072917 77799761 + 12 31915285 209381596 + 13 160861461 150208417 + 14 239475989 160552609 + 15 267758869 103276252 + 16 245710101 249137490 + 17 173329685 203048707 + 18 50617621 246253807 + 19 146009365 262586389 + 20 191069461 6905206 + 21 185797909 178835986 + 22 130194709 6852329 + 23 24259861 201065723 + 24 136428821 435527 + 25 198266133 125558990 + 26 209771797 94316688 + 27 170945813 100920717 + 28 81788181 210608325 + 29 210734357 91206455 + 30 20913429 86873573 + 31 49196309 144357837 + 32 27147541 81432560 + 33 223202581 123767117 + 34 100490517 99620582 + 35 195882261 245148089 + 36 240942357 130659015 + 37 235670805 2794256 + 38 180067605 247655316 + 39 74132757 48627027 + 40 186301717 70732108 + 41 248139029 44712064 + 42 259644693 109204463 + 43 220818709 5436569 + 44 131661077 150967166 + 45 260607253 29073821 + 46 70786325 67801080 + 47 99069205 29347725 + 48 77020437 83808605 + 49 4640021 3868264 + 50 150363413 238542765 + 51 245755157 34131501 + 52 22379797 118571753 + 53 17108245 17084127 + 54 229940501 199656207 + 55 124005653 201994363 + 56 236174613 236136481 + 57 29576469 116710147 + 58 41082133 66239007 + 59 2256149 178271862 + 60 181533973 148947207 + 61 42044693 82299604 + 62 120659221 221824219 + 63 148942101 145170461 + 64 126893333 106319258 + 65 54512917 230276690 + 66 200236309 244638532 + 67 27192597 16461170 + 68 72252693 89132506 + 69 66981141 71759229 + 70 11377941 249779547 + 71 173878549 242785907 + 72 17612053 78266823 + 73 79449365 191606869 + 74 90955029 83909406 + 75 52129045 201044770 + 76 231406869 54602080 + 77 91917589 100937434 + 78 170532117 130561166 + 79 198814997 73444221 + 80 176766229 267453607 + 81 104385813 116175116 + 82 250109205 236396971 + 83 77065493 42190726 + 84 122125589 160830363 + 85 116854037 16873195 + 86 61250837 248078966 + 87 223751445 21055291 + 88 67484949 252483132 + 89 129322261 119455863 + 90 140827925 12269293 + 91 102001941 192244382 + 92 12844309 254856330 + 93 141790485 203476400 + 94 220405013 180936465 + 95 248687893 201093549 + 96 226639125 148829828 + 97 154258709 48488086 + 98 31546645 63871715 + 99 126938389 229809258 + 100 171998485 183718956 + 101 166726933 239350569 + 102 111123733 44608097 + 103 5188885 192162516 + 104 117357845 71968129 + 105 179195157 18746217 + 106 190700821 238243212 + 107 151874837 1924330 + 108 62717205 62892675 + 109 191663381 239970134 + 110 1842453 223003749 + 111 30125333 109736622 + 112 8076565 137372466 + 113 204131605 145704688 + 114 81419541 113987306 + 115 176811285 160934942 + 116 221871381 7851917 + 117 216599829 52374071 + 118 160996629 26291484 + 119 55061781 69290300 + 120 167230741 192081814 + 121 229068053 7967019 + 122 240573717 75013883 + 123 201747733 17009158 + 124 112590101 134071116 + 125 241536277 60472268 + 126 51715349 106816648 + 127 79998229 186297982 + 128 57949461 83135151 + 129 254004501 257878554 + 130 131292437 236797377 + 131 226684181 222492322 + 132 3308821 20153791 + 133 266472725 111303701 + 134 210869525 43182759 + 135 104934677 39363188 + 136 217103637 194442363 + 137 10505493 205607358 + 138 22011157 177941307 + 139 251620629 87552498 + 140 162462997 50009829 + 141 22973717 51907347 + 142 101588245 219299707 + 143 129871125 12395806 + 144 107822357 104606972 + 145 35441941 235063317 + 146 181165333 13920104 + 147 8121621 264535031 + 148 53181717 70678208 + 149 47910165 266193092 + 150 260742421 213771010 + 151 154807573 220870268 + 152 266976533 197538864 + 153 60378389 193285408 + 154 71884053 128643658 + 155 33058069 63607983 + 156 212335893 197633358 + 157 72846613 64329001 + 158 151461141 142071102 + 159 179744021 243390094 + 160 157695253 51841561 + 161 85314837 195748063 + 162 231038229 100715487 + 163 57994517 137116699 + 164 103054613 9478801 + 165 97783061 98660418 + 166 42179861 119674414 + 167 204680469 195429716 + 168 48413973 51424950 + 169 110251285 89490258 + 170 121756949 45610025 + 171 82930965 63664699 + 172 262208789 158559879 + 173 122719509 216226319 + 174 201334037 262055377 + 175 229616917 192463566 + 176 207568149 43328006 + 177 135187733 258421881 + 178 12475669 78801703 + 179 107867413 227161871 + 180 152927509 223480114 + 181 147655957 264065680 + 182 92052757 147817513 + 183 254553365 81530620 + 184 98286869 143025163 + 185 160124181 12710996 + 186 171629845 47329496 + 187 132803861 206211735 + 188 43646229 51278481 + 189 172592405 89217477 + 190 251206933 160870708 + 191 11054357 246540767 + 192 257441045 197555395 + 193 185060629 4702947 + 194 62348565 66667838 + 195 157740309 116288723 + 196 202800405 25864867 + 197 197528853 75591598 + 198 141925653 148253940 + 199 35990805 266097525 + 200 148159765 53957680 + 201 209997077 81436710 + 202 221502741 252291159 + 203 182676757 72867267 + 204 93519125 262713706 + 205 222465301 70227019 + 206 32644373 225441640 + 207 60927253 255675327 + 208 38878485 96141905 + 209 234933525 89951261 + 210 112221461 182802981 + 211 207613205 191421799 + 212 252673301 71993060 + 213 247401749 188598172 + 214 191798549 239472783 + 215 85863701 62313149 + 216 198032661 171147045 + 217 259869973 145721032 + 218 2940181 242113191 + 219 232549653 50555839 + 220 143392021 106048275 + 221 3902741 9308578 + 222 82517269 37386347 + 223 110800149 69920879 + 224 88751381 126012078 + 225 16370965 95785000 + 226 162094357 8825308 + 227 257486101 34179275 + 228 34110741 211918326 + 229 28839189 184703579 + 230 241671445 3092218 + 231 135736597 125537493 + 232 247905557 76211434 + 233 41307413 55617595 + 234 52813077 135284678 + 235 13987093 257766540 + 236 193264917 236642188 + 237 53775637 24951240 + 238 132390165 252064830 + 239 160673045 76201967 + 240 138624277 137219547 + 241 66243861 140693250 + 242 211967253 200094819 + 243 38923541 31485696 + 244 83983637 27258839 + 245 78712085 182396905 + 246 23108885 94472246 + 247 185609493 37388733 + 248 29342997 156075392 + 249 91180309 198050941 + 250 102685973 50294709 + 251 63859989 7682088 + 252 243137813 236113621 + 253 103648533 235644094 + 254 182263061 182659809 + 255 210545941 124572223 + 256 188497173 248253400 + 257 116116757 74729644 + 258 261840149 69794234 + 259 88796437 33394692 + 260 133856533 173374600 + 261 128584981 31731783 + 262 72981781 95231041 + 263 235482389 184791413 + 264 79215893 260792549 + 265 141053205 154639247 + 266 152558869 105632372 + 267 113732885 224097940 + 268 24575253 222951663 + 269 153521429 223005316 + 270 232135957 216095828 + 271 260418837 65085279 + 272 238370069 40731813 + 273 165989653 16383270 + 274 43277589 4848098 + 275 138669333 158395352 + 276 183729429 231883785 + 277 178457877 119632757 + 278 122854677 123857692 + 279 16919829 149363710 + 280 129088789 240416538 + 281 190926101 43871601 + 282 202431765 151351299 + 283 163605781 220196816 + 284 74448149 47209944 + 285 203394325 105523994 + 286 13573397 202426520 + 287 41856277 16230224 + 288 19807509 170014787 + 289 215862549 84143216 + 290 93150485 123745497 + 291 188542229 256541308 + 292 233602325 52840026 + 293 228330773 27718003 + 294 172727573 30405831 + 295 66792725 49594710 + 296 178961685 213436447 + 297 240798997 252672547 + 298 252304661 37505122 + 299 213478677 114467804 + 300 124321045 95813009 + 301 253267221 1689216 + 302 63446293 260140971 + 303 91729173 96496144 + 304 69680405 217720496 + 305 265735445 128063114 + 306 143023381 8104608 + 307 238415125 177886192 + 308 15039765 23167868 + 309 9768213 142912066 + 310 222600469 201800002 + 311 116665621 3973502 + 312 228834581 29905908 + 313 22236437 94224806 + 314 33742101 151018386 + 315 263351573 25399992 + 316 174193941 218814490 + 317 34704661 29990071 + 318 113319189 239292814 + 319 141602069 155936672 + 320 119553301 33902573 + 321 47172885 266632053 + 322 192896277 44849975 + 323 19852565 40919093 + 324 64912661 261356397 + 325 59641109 46833120 + 326 4037909 219658382 + 327 166538517 130989174 + 328 10272021 76749466 + 329 72109333 223888376 + 330 83614997 73509265 + 331 44789013 71482469 + 332 224066837 266268019 + 333 84577557 40480189 + 334 163192085 258371137 + 335 191474965 44605440 + 336 169426197 5485562 + 337 97045781 81468207 + 338 242769173 84035230 + 339 69725461 232564553 + 340 114785557 80588334 + 341 109514005 126405710 + 342 53910805 202470057 + 343 216411413 12259902 + 344 60144917 204020751 + 345 121982229 223281434 + 346 133487893 191902304 + 347 94661909 102768865 + 348 5504277 88227229 + 349 134450453 151648659 + 350 213064981 167429572 + 351 241347861 149426992 + 352 219299093 250958551 + 353 146918677 227931577 + 354 24206613 244149462 + 355 119598357 66005293 + 356 164658453 136223679 + 357 159386901 231683468 + 358 103783701 288660 + 359 266284309 34710230 + 360 110017813 261773396 + 361 171855125 210893068 + 362 183360789 87815679 + 363 144534805 237748269 + 364 55377173 71616662 + 365 184323349 213549113 + 366 262937877 84957207 + 367 22785301 52019505 + 368 736533 83504261 + 369 196791573 19204883 + 370 74079509 106810845 + 371 169471253 196601313 + 372 214531349 9880608 + 373 209259797 212720026 + 374 153656597 38735 + 375 47721749 48393791 + 376 159890709 100061033 + 377 221728021 36776910 + 378 233233685 148173934 + 379 194407701 58038857 + 380 105250069 66489951 + 381 234196245 76235183 + 382 44375317 129443131 + 383 72658197 139307521 + 384 50609429 158482690 + 385 246664469 110648125 + 386 123952405 58943924 + 387 219344149 205970789 + 388 264404245 88483665 + 389 259132693 188004472 + 390 203529493 51773914 + 391 97594645 171799671 + 392 209763605 105808206 + 393 3165461 87857505 + 394 14671125 223030702 + 395 244280597 219000629 + 396 155122965 191336184 + 397 15633685 126631414 + 398 94248213 150940974 + 399 122531093 261344673 + 400 100482325 57512015 + 401 28101909 83879480 + 402 173825301 219037787 + 403 781589 212602810 + 404 45841685 222086483 + 405 40570133 7590439 + 406 253402389 5547829 + 407 147467541 254981503 + 408 259636501 129068547 + 409 53038357 214188483 + 410 64544021 162439613 + 411 25718037 33816306 + 412 204995861 27773537 + 413 65506581 214791436 + 414 144121109 267939825 + 415 172403989 268184593 + 416 150355221 167516780 + 417 77974805 57388034 + 418 223698197 168710610 + 419 50654485 66551006 + 420 95714581 260742692 + 421 90443029 58402469 + 422 34839829 248285025 + 423 197340437 147992919 + 424 41073941 19895689 + 425 102911253 265823477 + 426 114416917 84889756 + 427 75590933 157845886 + 428 254868757 231162010 + 429 115379477 190768882 + 430 193994005 62057860 + 431 222276885 9880913 + 432 200228117 70115161 + 433 127847701 149662876 + 434 5135637 26451482 + 435 100527381 154739922 + 436 145587477 54505925 + 437 140315925 190494195 + 438 84712725 93168220 + 439 247213333 237758463 + 440 90946837 165214174 + 441 152784149 92816119 + 442 164289813 108870219 + 443 125463829 172707546 + 444 36306197 114684324 + 445 165252373 173052840 + 446 243866901 188655079 + 447 3714325 141793634 + 448 250101013 152231702 + 449 177720597 210757638 + 450 55008533 179184945 + 451 150400277 58787734 + 452 195460373 258736182 + 453 190188821 253919249 + 454 134585621 195557415 + 455 28650773 105896312 + 456 140819733 146642179 + 457 202657045 82090953 + 458 214162709 84434634 + 459 175336725 196890374 + 460 86179093 65265021 + 461 215125269 11696942 + 462 25304341 229349659 + 463 53587221 245540930 + 464 31538453 263920036 + 465 227593493 90725952 + 466 104881429 208529176 + 467 200273173 165618986 + 468 245333269 186616183 + 469 240061717 98731263 + 470 184458517 137070786 + 471 78523669 139331008 + 472 190692629 82668792 + 473 252529941 83701611 + 474 264035605 130072089 + 475 225209621 80448002 + 476 136051989 201393190 + 477 264998165 93625732 + 478 75177237 34195230 + 479 103460117 171176434 + 480 81411349 255233793 + 481 9030933 176492363 + 482 154754325 232973263 + 483 250146069 56851854 + 484 26770709 225070473 + 485 21499157 111854782 + 486 234331413 36197421 + 487 128396565 188116184 + 488 240565525 91783101 + 489 33967381 216137182 + 490 45473045 95836217 + 491 6647061 210304975 + 492 185924885 104687007 + 493 46435605 457387 + 494 125050133 258551793 + 495 153333013 37189234 + 496 131284245 244662062 + 497 58903829 49675045 + 498 204627221 102570838 + 499 31583509 119410883 + 500 76643605 224152682 + 501 71372053 143343436 + 502 15768853 11426409 + 503 178269461 102305472 + 504 22002965 24038739 + 505 83840277 61015840 + 506 95345941 100216104 + 507 56519957 168079467 + 508 235797781 162071016 + 509 96308501 119116449 + 510 174923029 215602068 + 511 203205909 230503874 + 512 181157141 82258475 + 513 108776725 97198543 + 514 254500117 204246445 + 515 81456405 203349703 + 516 126516501 33916443 + 517 121244949 43250858 + 518 65641749 181246836 + 519 228142357 387960 + 520 71875861 266360248 + 521 133713173 5262130 + 522 145218837 261700839 + 523 106392853 72260567 + 524 17235221 223598850 + 525 146181397 31221095 + 526 224795925 23835143 + 527 253078805 64303074 + 528 231030037 154947576 + 529 158649621 169116489 + 530 35937557 119618261 + 531 131329301 158721947 + 532 176389397 41286300 + 533 171117845 198501592 + 534 115514645 127276879 + 535 9579797 852737 + 536 121748757 131930349 + 537 183586069 167365140 + 538 195091733 161908598 + 539 156265749 41337363 + 540 67108117 139324139 + 541 196054293 123695869 + 542 6233365 70175563 + 543 34516245 193946835 + 544 12467477 44347542 + 545 208522517 115482515 + 546 85810453 235610828 + 547 181202197 104016703 + 548 226262293 96315885 + 549 220990741 190713814 + 550 165387541 236441082 + 551 59452693 222188889 + 552 171621653 7673586 + 553 233458965 128943046 + 554 244964629 187763925 + 555 206138645 193798943 + 556 116981013 27735972 + 557 245927189 246594403 + 558 56106261 204676958 + 559 84389141 201053331 + 560 62340373 137382915 + 561 258395413 54785709 + 562 135683349 133842323 + 563 231075093 157723059 + 564 7699733 49058831 + 565 2428181 138376613 + 566 215260437 90357621 + 567 109325589 246014593 + 568 221494549 12079047 + 569 14896405 8484937 + 570 26402069 189320453 + 571 256011541 111263483 + 572 166853909 7323437 + 573 27364629 249970330 + 574 105979157 8957505 + 575 134262037 204111651 + 576 112213269 15671872 + 577 39832853 105515160 + 578 185556245 201237290 + 579 12512533 169894648 + 580 57572629 18004224 + 581 52301077 159979075 + 582 265133333 75951040 + 583 159198485 190818937 + 584 2931989 263635821 + 585 64769301 192915355 + 586 76274965 16631812 + 587 37448981 180655528 + 588 216726805 196575622 + 589 77237525 252312736 + 590 155852053 138377204 + 591 184134933 53175427 + 592 162086165 66138957 + 593 89705749 117724498 + 594 235429141 19413905 + 595 62385429 259020556 + 596 107445525 121641153 + 597 102173973 105574833 + 598 46570773 43274972 + 599 209071381 175091009 + 600 52804885 75526626 + 601 114642197 263852477 + 602 126147861 56622547 + 603 87321877 252028708 + 604 266599701 177110703 + 605 127110421 103675254 + 606 205724949 174554231 + 607 234007829 135169203 + 608 211959061 138837802 + 609 139578645 209902812 + 610 16866581 243732169 + 611 112258325 6718960 + 612 157318421 210023250 + 613 152046869 93652975 + 614 96443669 110818503 + 615 258944277 48884441 + 616 102677781 103111463 + 617 164515093 71349935 + 618 176020757 159346290 + 619 137194773 175436656 + 620 48037141 67417769 + 621 176983317 190982428 + 622 255597845 235977674 + 623 15445269 31711156 + 624 261831957 83822039 + 625 189451541 232103734 + 626 66739477 187374800 + 627 162131221 68349860 + 628 207191317 133204147 + 629 201919765 242702589 + 630 146316565 128635266 + 631 40381717 199123778 + 632 152550677 196443964 + 633 214387989 2332273 + 634 225893653 174856673 + 635 187067669 69368460 + 636 97910037 254421362 + 637 226856213 95852434 + 638 37035285 172701166 + 639 65318165 129725828 + 640 43269397 19580757 + 641 239324437 34380896 + 642 116612373 237266343 + 643 212004117 25531432 + 644 257064213 9672932 + 645 251792661 134341851 + 646 196189461 215214349 + 647 90254613 207427194 + 648 202423573 205577761 + 649 264260885 175288579 + 650 7331093 221642785 + 651 236940565 52313208 + 652 147782933 51304203 + 653 8293653 205209817 + 654 86908181 103213793 + 655 115191061 10831396 + 656 93142293 64603042 + 657 20761877 3658843 + 658 166485269 243460430 + 659 261877013 265188220 + 660 38501653 226354150 + 661 33230101 155495306 + 662 246062357 220609384 + 663 140127509 192283778 + 664 252296469 249001942 + 665 45698325 171837030 + 666 57203989 149758256 + 667 18378005 242759989 + 668 197655829 113426292 + 669 58166549 100672751 + 670 136781077 146004644 + 671 165063957 61952404 + 672 143015189 68942527 + 673 70634773 258426661 + 674 216358165 56010693 + 675 43314453 100502945 + 676 88374549 96430519 + 677 83102997 156216584 + 678 27499797 263309460 + 679 190000405 3747162 + 680 33733909 176770140 + 681 95571221 110466712 + 682 107076885 77692175 + 683 68250901 222326977 + 684 247528725 22405805 + 685 108039445 169165781 + 686 186653973 151127351 + 687 214936853 133142484 + 688 192888085 151088300 + 689 120507669 111867071 + 690 266231061 61841676 + 691 93187349 186835605 + 692 138247445 6826584 + 693 132975893 254994774 + 694 77372693 193368207 + 695 239873301 28741890 + 696 83606805 107371441 + 697 145444117 109666714 + 698 156949781 123933630 + 699 118123797 109503261 + 700 28966165 165167287 + 701 157912341 260742539 + 702 236526869 237071002 + 703 264809749 74455268 + 704 242760981 161093993 + 705 170380565 219340073 + 706 47668501 111007012 + 707 143060245 105804377 + 708 188120341 76031433 + 709 182848789 33448052 + 710 127245589 129274714 + 711 21310741 117321595 + 712 133479701 159294934 + 713 195317013 19490668 + 714 206822677 138536253 + 715 167996693 22777929 + 716 78839061 123328912 + 717 207785237 225456657 + 718 17964309 253889230 + 719 46247189 4379845 + 720 24198421 217448695 + 721 220253461 162463843 + 722 97541397 53560331 + 723 192933141 244333805 + 724 237993237 154098698 + 725 232721685 146936418 + 726 177118485 189518069 + 727 71183637 119539907 + 728 183352597 182594251 + 729 245189909 226863118 + 730 256695573 239989132 + 731 217869589 80640069 + 732 128711957 15379769 + 733 257658133 181797223 + 734 67837205 51635665 + 735 96120085 41405301 + 736 74071317 170206036 + 737 1690901 59727470 + 738 147414293 7990722 + 739 242806037 184042065 + 740 19430677 91082012 + 741 14159125 177078049 + 742 226991381 224151904 + 743 121056533 153885915 + 744 233225493 27323024 + 745 26627349 44966785 + 746 38133013 9910444 + 747 267742485 133143313 + 748 178584853 228244402 + 749 39095573 248253326 + 750 117710101 17234852 + 751 145992981 35585269 + 752 123944213 137855105 + 753 51563797 29620040 + 754 197287189 92787273 + 755 24243477 43418246 + 756 69303573 5470461 + 757 64032021 242362031 + 758 8428821 83229852 + 759 170929429 70413251 + 760 14662933 80405798 + 761 76500245 129161667 + 762 88005909 103660187 + 763 49179925 30341294 + 764 228457749 75105531 + 765 88968469 6443140 + 766 167582997 740423 + 767 195865877 105408837 + 768 173817109 238884990 + 769 101436693 190630642 + 770 247160085 158003616 + 771 74116373 209386890 + 772 119176469 15753134 + 773 113904917 192841997 + 774 58301717 153676455 + 775 220802325 256046459 + 776 64535829 191896203 + 777 126373141 61065941 + 778 137878805 102856538 + 779 99052821 159158554 + 780 9895189 211323157 + 781 138841365 111726666 + 782 217455893 120641466 + 783 245738773 100929637 + 784 223690005 54913867 + 785 151309589 124377452 + 786 28597525 53693384 + 787 123989269 263566174 + 788 169049365 240419119 + 789 163777813 147007035 + 790 108174613 17109890 + 791 2239765 23968260 + 792 114408725 211847872 + 793 176246037 227604151 + 794 187751701 125988585 + 795 148925717 101213270 + 796 59768085 218515454 + 797 188714261 145722080 + 798 267328789 226991613 + 799 27176213 140636758 + 800 5127445 241301737 + 801 201182485 217785014 + 802 78470421 166781119 + 803 173862165 56009730 + 804 218922261 261086592 + 805 213650709 223346233 + 806 158047509 60454701 + 807 52112661 29538652 + 808 164281621 258749893 + 809 226118933 210394473 + 810 237624597 23109960 + 811 198798613 243429986 + 812 109640981 215171511 + 813 238587157 226918470 + 814 48766229 169844497 + 815 77049109 74583830 + 816 55000341 111231318 + 817 251055381 52471504 + 818 128343317 78884998 + 819 223735061 242077558 + 820 359701 196244642 + 821 263523605 3477767 + 822 207920405 133764520 + 823 101985557 122811268 + 824 214154517 182655898 + 825 7556373 127925996 + 826 19062037 181145208 + 827 248671509 167426878 + 828 159513877 51344960 + 829 20024597 205369469 + 830 98639125 67689204 + 831 126922005 21259942 + 832 104873237 51627155 + 833 32492821 15361467 + 834 178216213 176929565 + 835 5172501 134952379 + 836 50232597 164382355 + 837 44961045 142761638 + 838 257793301 87092979 + 839 151858453 153839740 + 840 264027413 102054975 + 841 57429269 98687806 + 842 68934933 181712503 + 843 30108949 260128491 + 844 209386773 113960345 + 845 69897493 199564163 + 846 148512021 39014823 + 847 176794901 99154182 + 848 154746133 180978336 + 849 82365717 224943989 + 850 228089109 42532996 + 851 55045397 121558735 + 852 100105493 15553364 + 853 94833941 222816020 + 854 39230741 38929167 + 855 201731349 241113156 + 856 45464853 135436213 + 857 107302165 241168992 + 858 118807829 143300934 + 859 79981845 103152999 + 860 259259669 253071298 + 861 119770389 59556185 + 862 198384917 202310442 + 863 226667797 158320182 + 864 204619029 80903037 + 865 132238613 262837247 + 866 9526549 62619836 + 867 104918293 51950259 + 868 149978389 136682213 + 869 144706837 93694546 + 870 89103637 107762170 + 871 251604245 234685148 + 872 95337749 132853242 + 873 157175061 136987730 + 874 168680725 184399589 + 875 129854741 83424947 + 876 40697109 50295996 + 877 169643285 172270079 + 878 248257813 139194237 + 879 8105237 48811575 + 880 254491925 138325802 + 881 182111509 247530329 + 882 59399445 87243715 + 883 154791189 44616039 + 884 199851285 109387078 + 885 194579733 142321760 + 886 138976533 143645621 + 887 33041685 253044805 + 888 145210645 212795151 + 889 207047957 173068564 + 890 218553621 155062100 + 891 179727637 50997967 + 892 90570005 160994437 + 893 219516181 119324021 + 894 29695253 236590753 + 895 57978133 157552903 + 896 35929365 203300264 + 897 231984405 29076867 + 898 109272341 234893722 + 899 204664085 218045163 + 900 249724181 52157047 + 901 244452629 218751294 + 902 188849429 265068608 + 903 82914581 146245757 + 904 195083541 225315572 + 905 256920853 199465126 + 906 268426517 173777555 + 907 229600533 124361147 + 908 140442901 166784798 + 909 953621 19207100 + 910 79568149 76118164 + 911 107851029 66162343 + 912 85802261 125880053 + 913 13421845 262836862 + 914 159145237 87188033 + 915 254536981 153855807 + 916 31161621 83481209 + 917 25890069 173036781 + 918 238722325 53649307 + 919 132787477 32777093 + 920 244956437 20468137 + 921 38358293 66231049 + 922 49863957 90599587 + 923 11037973 153568120 + 924 190315797 186156167 + 925 50826517 258843858 + 926 129441045 44701015 + 927 157723925 161564439 + 928 135675157 24554258 + 929 63294741 261993032 + 930 209018133 31051192 + 931 35974421 238972516 + 932 81034517 53413194 + 933 75762965 123667307 + 934 20159765 164747719 + 935 182660373 31127901 + 936 26393877 253612847 + 937 88231189 160290875 + 938 99736853 24017282 + 939 60910869 257107972 + 940 240188693 69162176 + 941 100699413 151417016 + 942 179313941 260828394 + 943 207596821 25377367 + 944 185548053 17811967 + 945 113167637 145034466 + 946 258891029 184972287 + 947 85847317 55013464 + 948 130907413 80442091 + 949 125635861 189131961 + 950 70032661 179982018 + 951 232533269 259787269 + 952 76266773 237932420 + 953 138104085 63262781 + 954 149609749 92519729 + 955 110783765 16598880 + 956 21626133 202727370 + 957 150572309 83851118 + 958 229186837 37683021 + 959 257469717 44525671 + 960 235420949 224142268 + 961 163040533 30450252 + 962 40328469 130569495 + 963 135720213 257338652 + 964 180780309 14621532 + 965 175508757 219484375 + 966 119905557 217841293 + 967 13970709 31937918 + 968 126139669 91915945 + 969 187976981 162071311 + 970 199482645 146160560 + 971 160656661 87400844 + 972 71499029 168469923 + 973 200445205 174635252 + 974 10624277 30624897 + 975 38907157 69062984 + 976 16858389 225163338 + 977 212913429 36729478 + 978 90201365 254767358 + 979 185593109 159130800 + 980 230653205 242876061 + 981 225381653 64778181 + 982 169778453 128379176 + 983 63843605 2939846 + 984 176012565 202487966 + 985 237849877 38334641 + 986 249355541 34993407 + 987 210529557 51132040 + 988 121371925 84878924 + 989 250318101 5387594 + 990 60497173 89707652 + 991 88780053 217478392 + 992 66731285 139364263 + 993 262786325 13925776 + 994 140074261 139184053 + 995 235466005 147314452 + 996 12090645 78388399 + 997 6819093 111937924 + 998 219651349 30084755 + 999 113716501 22846686 + 1000 225885461 151266659 +=end= +""" + file.close() + +fileTypes.append(qmcchemFile) + +def write_qmcchemFile(f): + assert (isinstance(f,resultsFile)) + if "MCSCF" in f.mo_sets.keys() : + MoType = "MCSCF" + else: + MoType = None + f.clean_uncontractions() + print "Writing qmcin" + qmcchem_write_qmcin(f,method='VMC', JastType='None') + print "Writing qmc.seed" + qmcchem_write_qmc_seed() + print "Writing qmc.basis" + qmcchem_write_qmc_basis(f) + print "Writing qmc.mos" + qmcchem_write_qmc_mos(f,MoType) + print "Writing qmc.jast_det" + qmcchem_write_qmc_jast_det(f) + print "Writing qmc.det" + qmcchem_write_qmc_det(f,MoType) + diff --git a/resultsFile/Modules/qmcchem_newFile.py b/resultsFile/Modules/qmcchem_newFile.py new file mode 100755 index 0000000..70994d5 --- /dev/null +++ b/resultsFile/Modules/qmcchem_newFile.py @@ -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) + diff --git a/resultsFile/Modules/wfnFile.py b/resultsFile/Modules/wfnFile.py new file mode 100755 index 0000000..934876f --- /dev/null +++ b/resultsFile/Modules/wfnFile.py @@ -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 pos0: + 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) + diff --git a/resultsFile/Modules/xmvbFile.py b/resultsFile/Modules/xmvbFile.py new file mode 100755 index 0000000..569231b --- /dev/null +++ b/resultsFile/Modules/xmvbFile.py @@ -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 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 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 ##### + diff --git a/resultsFile/__init__.py b/resultsFile/__init__.py new file mode 100755 index 0000000..815e7a8 --- /dev/null +++ b/resultsFile/__init__.py @@ -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 " +__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 + + diff --git a/resultsFile/cython_setup b/resultsFile/cython_setup new file mode 100755 index 0000000..d2575e3 --- /dev/null +++ b/resultsFile/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("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 + ) + diff --git a/resultsFile/getFile.py b/resultsFile/getFile.py new file mode 100755 index 0000000..a21384b --- /dev/null +++ b/resultsFile/getFile.py @@ -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]) + diff --git a/resultsFile/resultsFile.py b/resultsFile/resultsFile.py new file mode 100755 index 0000000..06e06d2 --- /dev/null +++ b/resultsFile/resultsFile.py @@ -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