# -*- coding: iso-8859-1 -*-
'''Module for reading the output files of quantum chemical software.
'''
'''
orbkit
Gunter Hermann, Vincent Pohl, and Axel Schild
Institut fuer Chemie und Biochemie, Freie Universitaet Berlin, 14195 Berlin, Germany
This file is part of orbkit.
orbkit is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as
published by the Free Software Foundation, either version 3 of
the License, or any later version.
orbkit 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 Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with orbkit. If not, see <http://www.gnu.org/licenses/>.
'''
import os
import copy
import numpy
from orbkit.core import l_deg, lquant, orbit, exp, exp_wfn,create_mo_coeff
from orbkit.display import display
from orbkit.qcinfo import QCinfo, get_atom_symbol, load
[docs]def main_read(filename,itype='molden',all_mo=False,spin=None,cclib_parser=None,
**kwargs):
'''Calls the requested read function.
**Parameters:**
filename : str
Specifies the filename for the input file.
itype : str, choices={'molden', 'gamess', 'gaussian.log', 'gaussian.fchk', 'aomix', 'cclib'}
Specifies the type of the input file.
all_mo : bool, optional
If True, all molecular orbitals are returned.
spin : {None, 'alpha', or 'beta'}, optional
If not None, returns exclusively 'alpha' or 'beta' molecular orbitals.
cclib_parser : str
If itype is 'cclib', specifies the cclib.parser.
**Returns:**
qc (class QCinfo) with attributes geo_spec, geo_info, ao_spec, mo_spec, etot :
See :ref:`Central Variables` for details.
**Note:**
All additional keyword arguments are forwarded to the reading functions.
'''
display('Opened \n\t%s\n' % filename)
# What kind of input file has to be read?
reader = {'molden': read_molden,
'gamess': read_gamess,
'gaussian.log': read_gaussian_log,
'gaussian.fchk': read_gaussian_fchk,
'aomix' : read_aomix,
'cclib' : read_with_cclib,
'wfn': read_wfn,
'wfx': read_wfx,
'orbkit.dump': load}
display('Loading %s file...' % itype)
if itype not in reader.keys():
display('Available reader (`itype`) are:\n ' + ', '.join(reader.keys()))
raise NotImplementedError("itype='%s' not implemented!"%itype)
# Return required data
return reader[itype](filename, all_mo=all_mo, spin=spin,
cclib_parser=cclib_parser,**kwargs)
# main_read
[docs]def read_molden(filename, all_mo=False, spin=None, i_md=-1, interactive=True,
**kwargs):
'''Reads all information desired from a molden file.
**Parameters:**
filename : str
Specifies the filename for the input file.
all_mo : bool, optional
If True, all molecular orbitals are returned.
spin : {None, 'alpha', or 'beta'}, optional
If not None, returns exclusively 'alpha' or 'beta' molecular orbitals.
i_md : int, default=-1
Selects the `[Molden Format]` section of the output file.
interactive : bool
If True, the user is asked to select the different sets.
**Returns:**
qc (class QCinfo) with attributes geo_spec, geo_info, ao_spec, mo_spec, etot :
See :ref:`Central Variables` for details.
'''
fid = open(filename,'r') # Open the file
flines = fid.readlines() # Read the WHOLE file into RAM
fid.close() # Close the file
# Is this really a molden file?
ismd = False
for line in flines:
if '[molden format]' in line.lower():
ismd = True
if not ismd:
display('The input file %s is no valid molden file!\n\nIt does' % filename+
' not contain the keyword: [Molden Format]\n')
raise IOError('Not a valid input file')
def check_sel(count,i,interactive=False):
if count == 0:
raise IndexError
elif count == 1:
return 0
message = '\tPlease give an integer from 0 to %d: ' % (count-1)
try:
if interactive:
i = int(raw_input(message))
i = range(count)[i]
except (IndexError,ValueError):
raise IOError(message.replace(':','!'))
else:
display('\tSelecting the %s' %
('last element.' if (i == count-1) else 'element %d.' % i))
return i
has_alpha = []
has_beta = []
restricted = []
cartesian_basis = []
mixed_warning = []
by_orca = []
count = 0
# Go through the file line by line
for il in range(len(flines)):
line = flines[il] # The current line as string
# Check the file for keywords
if '[molden format]' in line.lower():
count += 1
has_alpha.append(False)
has_beta.append(False)
restricted.append(False)
cartesian_basis.append(True)
mixed_warning.append(False)
by_orca.append(False)
if 'orca' in line.lower():
by_orca[-1] = True
if '[5d]' in line.lower() or '[5d7f]' in line.lower():
cartesian_basis[-1] = False
if '[5d10f]' in line.lower():
mixed_warning[-1] = '5D, 10F'
cartesian_basis[-1] = False
if '[7f]' in line.lower():
mixed_warning[-1] = '6D, 7F'
cartesian_basis[-1] = True
if 'Spin' in line and 'alpha' in line.lower():
has_alpha[-1] = True
if 'Spin' in line and 'beta' in line.lower():
has_beta[-1] = True
if 'Occup' in line:
restricted[-1] = restricted[-1] or (float(line.split('=')[1]) > 1.+1e-4)
if count == 0:
display('The input file %s is no valid molden file!\n\nIt does' % filename +
' not contain the keyword: [Molden Format]\n')
raise IOError('Not a valid input file')
else:
if count > 1:
display('\nContent of the molden file:')
display('\tFound %d [Molden Format] keywords, i.e., ' % count +
'this file contains %d molden files.' % count)
i_md = check_sel(count,i_md,interactive=interactive)
if spin is not None:
if restricted[i_md]:
raise IOError('The keyword `spin` is only supported for unrestricted calculations.')
if spin != 'alpha' and spin != 'beta':
raise IOError('`spin=%s` is not a valid option' % spin)
elif spin == 'alpha' and has_alpha[i_md]:
display('Reading only molecular orbitals of spin alpha.')
elif spin == 'beta' and has_beta[i_md]:
display('Reading only molecular orbitals of spin beta.')
elif (not has_alpha[i_md]) and (not has_beta[i_md]):
raise IOError(
'Molecular orbitals in `molden` file do not contain `Spin=` keyword')
elif ((spin == 'alpha' and not has_alpha[i_md]) or
(spin == 'beta' and not has_beta[i_md])):
raise IOError('You requested `%s` orbitals, but None of them are present.'
% spin)
# Set a counter for the AOs
basis_count = 0
sym = {}
# Declare synonyms for molden keywords
synonyms = {'Sym': 'sym',
'Ene': 'energy',
'Occup': 'occ_num',
'Spin': 'spin'
}
MO_keys = synonyms.keys()
count = 0
max_l = 0
start_reading = False
# Go through the file line by line
for il in range(len(flines)):
line = flines[il] # The current line as string
thisline = line.split() # The current line split into segments
# Check the file for keywords
if '[molden format]' in line.lower():
# A new file begins
# Initialize the variables
if i_md == count:
qc = QCinfo()
sec_flag = False # A Flag specifying the current section
start_reading = True # Found the selected section
else:
start_reading = False
count += 1
continue
if start_reading:
if '_ENERGY=' in line:
try:
qc.etot = float(thisline[1])
except IndexError:
pass
elif '[atoms]' in line.lower():
# The section containing information about
# the molecular geometry begins
sec_flag = 'geo_info'
if 'Angs' in line:
# The length are given in Angstroem
# and have to be converted to Bohr radii --
aa_to_au = 1/0.52917720859
else:
# The length are given in Bohr radii
aa_to_au = 1.0
elif '[gto]' in line.lower():
# The section containing information about
# the atomic orbitals begins
sec_flag = 'ao_info'
bNew = True # Indication for start of new AO section
elif '[mo]' in line.lower():
# The section containing information about
# the molecular orbitals begins
sec_flag = 'mo_info'
bNew = True # Indication for start of new MO section
elif '[sto]' in line.lower():
# The orbkit does not support Slater type orbitals
display('orbkit does not work for STOs!\nEXIT\n');
raise IOError('Not a valid input file')
elif '[' in line:
sec_flag = None
else:
# Check if we are in a specific section
if sec_flag == 'geo_info' and thisline != []:
# Geometry section
qc.geo_info.append(thisline[0:3])
qc.geo_spec.append([float(ii)*aa_to_au for ii in thisline[3:]])
if sec_flag == 'ao_info':
# Atomic orbital section
def check_int(i):
try:
int(i)
return True
except ValueError:
return False
if thisline == []:
# There is a blank line after every AO
bNew = True
elif bNew:
# The following AOs are for which atom?
bNew = False
at_num = int(thisline[0]) - 1
ao_num = 0
elif len(thisline) == 3 and check_int(thisline[1]):
# AO information section
# Initialize a new dict for this AO
ao_num = 0 # Initialize number of atomic orbiatls
ao_type = thisline[0].lower() # Which type of atomic orbital do we have
pnum = int(thisline[1]) # Number of primatives
# Calculate the degeneracy of this AO and increase basis_count
for i_ao in ao_type:
# Calculate the degeneracy of this AO and increase basis_count
basis_count += l_deg(lquant[i_ao],cartesian_basis=cartesian_basis[i_md])
max_l = max(max_l,lquant[i_ao])
qc.ao_spec.append({'atom': at_num,
'type': i_ao,
'pnum': -pnum if by_orca[i_md] else pnum,
'coeffs': numpy.zeros((pnum, 2))
})
else:
# Append the AO coefficients
coeffs = numpy.array(line.replace('D','e').split(), dtype=numpy.float64)
for i_ao in range(len(ao_type)):
qc.ao_spec[-len(ao_type)+i_ao]['coeffs'][ao_num,:] = [coeffs[0],
coeffs[1+i_ao]]
ao_num += 1
if sec_flag == 'mo_info':
# Molecular orbital section
if '=' in line:
# MO information section
if bNew:
# Create a numpy array for the MO coefficients and
# for backward compability create a simple counter for 'sym'
qc.mo_spec.append({'coeffs': numpy.zeros(basis_count),
'sym': '%d.1' % (len(qc.mo_spec)+1)})
bNew = False
# Append information to dict of this MO
info = line.replace('\n','').replace(' ','')
info = info.split('=')
if info[0] in MO_keys:
if info[0] == 'Spin':
info[1] = info[1].lower()
elif info[0] != 'Sym':
info[1] = float(info[1])
elif not '.' in info[1]:
from re import search
try:
a = search(r'\d+', info[1]).group()
if a == info[1]:
info[1] = '%s.1' % a
elif info[1].startswith(a):
info[1] = info[1].replace(a, '%s.' % a,1)
else:
raise AttributeError
except AttributeError:
if info[1] not in sym.keys(): sym[info[1]] = 1
else: sym[info[1]] += 1
info[1] = '%d.%s' % (sym[info[1]],info[1])
qc.mo_spec[-1][synonyms[info[0]]] = info[1]
else:
if ('[' or ']') in line:
# start of another section that is not (yet) read
sec_flag = None
else:
# Append the MO coefficients
bNew = True # Reset bNew
index = int(thisline[0])-1
try:
# Try to convert coefficient to float
qc.mo_spec[-1]['coeffs'][index] = float(thisline[1])
except ValueError:
# If it cannot be converted print error message
display('Error in coefficient %d of MO %s!' % (index,
qc.mo_spec[-1]['sym']) + '\nSetting this coefficient to zero...')
# Spherical basis?
if not cartesian_basis[i_md]:
qc.ao_spherical = get_ao_spherical(qc.ao_spec,p=[1,0])
if max_l > 2 and mixed_warning[i_md]:
display('='*80)
display('The input file %s contains ' % filename +
'mixed spherical and Cartesian function (%s).' % mixed_warning[i_md] +
'ORBKIT does not support these basis functions yet. '+
'Pleas contact us, if you need this feature!')
display('='*80)
# Are all MOs requested for the calculation?
if not all_mo:
for i in range(len(qc.mo_spec))[::-1]:
if qc.mo_spec[i]['occ_num'] < 0.0000001:
del qc.mo_spec[i]
# Only molecular orbitals of one spin requested?
if spin is not None:
for i in range(len(qc.mo_spec))[::-1]:
if qc.mo_spec[i]['spin'] != spin:
del qc.mo_spec[i]
if restricted[i_md]:
# Closed shell calculation
for mo in qc.mo_spec:
del mo['spin']
else:
# Rename MOs according to spin
for mo in qc.mo_spec:
mo['sym'] += '_%s' % mo['spin'][0]
# Orca uses for all molecular orbitals the same name
sym = [i['sym'] for i in qc.mo_spec]
if sym[1:] == sym[:-1]:
sym = sym[0].split('.')[-1]
for i in range(len(qc.mo_spec)):
qc.mo_spec[i]['sym'] = '%d.%s' % (i+1,sym)
# Convert geo_info and geo_spec to numpy.ndarrays
qc.format_geo()
# Check the normalization
from orbkit.analytical_integrals import get_ao_overlap,get_lxlylz
norm = numpy.diagonal(get_ao_overlap(qc.geo_spec,qc.geo_spec,qc.ao_spec))
if max(numpy.abs(norm-1.)) > 1e-5:
display('The atomic orbitals are not normalized correctly, renormalizing...\n')
if not by_orca[i_md]:
j = 0
for i in range(len(qc.ao_spec)):
qc.ao_spec[i]['coeffs'][:,1] /= numpy.sqrt(norm[j])
for n in range(l_deg(lquant[qc.ao_spec[i]['type']],cartesian_basis=True)):
j += 1
else:
qc.ao_spec[0]['N'] = 1/numpy.sqrt(norm[:,numpy.newaxis])
if cartesian_basis[i_md]:
from orbkit.cy_overlap import ommited_cca_norm
cca = ommited_cca_norm(get_lxlylz(qc.ao_spec))
for mo in qc.mo_spec:
mo['coeffs'] *= cca
return qc
# read_molden
[docs]def read_gamess(filename, all_mo=False, spin=None, read_properties=False,
**kwargs):
'''Reads all information desired from a Gamess-US output file.
**Parameters:**
filename : str
Specifies the filename for the input file.
all_mo : bool, optional
If True, all molecular orbitals are returned.
**Returns:**
qc (class QCinfo) with attributes geo_spec, geo_info, ao_spec, mo_spec, etot :
See :ref:`Central Variables` for details.
'''
fid = open(filename,'r') # Open the file
flines = fid.readlines() # Read the WHOLE file into RAM
fid.close() # Close the filename
# Initialize the variables
qc = QCinfo()
has_alpha = False # Flag for alpha electron set
has_beta = False # Flag for beta electron set
restricted = True # Flag for restricted calculation
sec_flag = None # A Flag specifying the current section
is_pop_ana = True # Flag for population analysis for ground state
keyword = [' ATOM ATOMIC COORDINATES','\n']
# Keywords for single point calculation and
# geometry optimization
mokey = 'EIGENVECTORS' # Keyword for MOs
unrestopt = False # Flag for unrestricted optimization
bopt = False # Flag for geometry optimization
sym={} # Symmetry of MOs
geo_skip = 1 # Number of lines to skip in geometry section
for il in range(len(flines)):
line = flines[il] # The current line as string
thisline = line.split() # The current line split into segments
# Check the file for keywords
if 'RUNTYP=OPTIMIZE' in line:
keyword = [' COORDINATES OF ALL ATOMS ARE',
'***** EQUILIBRIUM GEOMETRY LOCATED *****']
geo_skip = 2
bopt = True
if 'SCFTYP=UHF' in line:
mokey = ' SET ****'
restricted = False
else:
mokey = 'MOLECULAR ORBITALS'
elif keyword[0] in line and keyword[1] in flines[il-1]:
# The section containing information about
# the molecular geometry begins
sec_flag = 'geo_info'
atom_count = 0 # Counter for Atoms
if '(BOHR)' in line:
# The length are given in Bohr radii
aa_to_au = 1.0
else:
# The length are given in Angstroem
# and have to be converted to Bohr radii
aa_to_au = 1/0.52917720859
elif 'ATOMIC BASIS SET' in line:
# The section containing information about
# the atomic orbitals begins
sec_flag = 'ao_info'
ao_skip = 6 # Number of lines to skip
AO = [] # Atomic orbitals
elif '----- ALPHA SET ' in line:
# The section for alpha electrons
has_alpha = True
has_beta = False
restricted = False
elif '----- BETA SET ' in line:
# The section for alpha electrons
restricted = False
has_alpha = False
has_beta = True
elif mokey in line and len(thisline) < 3:
# The section containing information about
# the molecular orbitals begins
sec_flag = 'mo_info'
mo_skip = 1
len_mo = 0 # Number of MOs
init_mo = False # Initialize new MO section
info_key = None # A Flag specifying the energy and symmetry section
exp_list = []
if 'ALPHA' in line:
has_alpha = True
mo_skip = 0
elif 'BETA' in line:
has_beta = True
has_alpha = False
mo_skip = 0
elif 'NATURAL ORBITALS' in line and len(thisline) <= 3:
display('The natural orbitals are not extracted.')
elif ' NUMBER OF OCCUPIED ORBITALS (ALPHA) =' in line:
occ = [] # occupation number of molecular orbitals
occ.append(int(thisline[-1]))
elif ' NUMBER OF OCCUPIED ORBITALS (BETA ) =' in line:
occ.append(int(thisline[-1]))
# elif 'ECP POTENTIALS' in line:
# sec_flag = 'ecp_info'
# ecp = ''
elif ' NUMBER OF OCCUPIED ORBITALS (ALPHA) KEPT IS =' in line:
occ = [] # occupation number of molecular orbitals
occ.append(int(thisline[-1]))
elif ' NUMBER OF OCCUPIED ORBITALS (BETA ) KEPT IS =' in line:
occ.append(int(thisline[-1]))
elif 'NUMBER OF STATES REQUESTED' in line and read_properties:
# get the number of excited states and initialize variables for
# transition dipole moment and energies
exc_states = int(line.split('=')[1]) # Number of excited states
# Dipole moments matrix: Diagonal elements -> permanent dipole moments
# Off-diagonal elements -> transition dipole moments
qc.dipole_moments = numpy.zeros(((exc_states+1),(exc_states+1),3))
# Multiplicity of ground and excited states
qc.states['multiplicity'] = numpy.zeros(exc_states+1)
# Energies of ground and excited states
qc.states['energy'] = numpy.zeros(exc_states+1)
qc.states['energy'][0] = qc.etot
qc.states['multiplicity'][0] = gs_multi
dm_flag = None # Flag specifying the dipole moments section
elif 'TRANSITION DIPOLE MOMENTS' in line and read_properties:
# Section containing energies of excited states
sec_flag = 'dm_info'
# Energy and Multiplicity for ground state
elif 'SPIN MULTIPLICITY' in line and read_properties:
# odd way to get gound state multiplicity
gs_multi = int(line.split()[3])
elif 'FINAL' in line and read_properties:
# get (last) energy
qc.etot = float(line.split()[4])
elif 'TOTAL MULLIKEN AND LOWDIN ATOMIC POPULATIONS' in line and is_pop_ana == True and read_properties:
# Read Mulliken and Lowdin Atomic Populations
sec_flag = 'pop_info'
pop_skip = 1
is_pop_ana == False
qc.pop_ana['Lowdin'] = []
qc.pop_ana['Mulliken'] = []
else:
# Check if we are in a specific section
if sec_flag == 'geo_info':
if not geo_skip:
if line == '\n':
sec_flag = None
else:
qc.geo_info.append([thisline[0],atom_count+1,thisline[1]])
qc.geo_spec.append([float(ii)*aa_to_au for ii in thisline[2:]])
atom_count += 1
elif geo_skip:
geo_skip -= 1
elif sec_flag == 'ao_info':
if not ao_skip:
if ' TOTAL NUMBER OF BASIS SET SHELLS' in line:
sec_flag = None
else:
if len(thisline) == 1:
# Read atom type
at_type = thisline[0]
AO.append([])
new_ao = False
elif len(thisline) == 0 and new_ao == False:
new_ao = True
else:
coeffs = [float(ii) for ii in thisline[3:]]
if new_ao:
ao_type = thisline[1].lower().replace('l','sp')
for i_ao,t_ao in enumerate(ao_type):
AO[-1].append({'atom_type': at_type,
'type': t_ao,
'pnum': 1,
'coeffs': [[coeffs[0],coeffs[1+i_ao]]]})
new_ao = False
else:
for i_ao in range(len(ao_type)):
AO[-1][-len(ao_type)+i_ao]['coeffs'].append([coeffs[0],
coeffs[1+i_ao]])
AO[-1][-len(ao_type)+i_ao]['pnum'] += 1
elif ao_skip:
ao_skip -= 1
elif sec_flag == 'mo_info':
if not mo_skip:
if 'END OF' in line and 'CALCULATION' in line or '-----------' in line:
sec_flag = None
has_alpha = False
has_beta = False
else:
if thisline == []:
info_key = None
init_mo = True
try:
int(flines[il+1].split()[0])
except ValueError:
sec_flag = None
init_mo = False
elif init_mo:
init_len = len(thisline)
exp_list = []
for ii in range(len(thisline)):
if has_alpha == True or has_beta == True:
qc.mo_spec.append({'coeffs': [],
'energy': 0.0,
'occ_num': 0.0,
'sym': '',
'spin': ''
})
else:
qc.mo_spec.append({'coeffs': [],
'energy': 0.0,
'occ_num': 0.0,
'sym': ''
})
init_mo = False
info_key = 'energy'
elif len(thisline) == init_len and info_key == 'energy':
for ii in range(init_len,0,-1):
qc.mo_spec[-ii]['energy'] = float(thisline[init_len-ii])
info_key = 'symmetry'
elif len(thisline) == init_len and info_key == 'symmetry':
for ii in range(init_len,0,-1):
len_mo += 1
a = thisline[init_len-ii]
if a not in sym.keys(): sym[a] = 1
else: sym[a] = len_mo
if has_alpha:
qc.mo_spec[-ii]['sym'] = '%d.%s_a' % (sym[a], thisline[init_len-ii])
qc.mo_spec[-ii]['spin'] = 'alpha'
elif has_beta:
qc.mo_spec[-ii]['sym'] = '%d.%s_b' % (sym[a], thisline[init_len-ii])
qc.mo_spec[-ii]['spin'] = 'beta'
else:
qc.mo_spec[-ii]['sym'] = '%d.%s' % (sym[a], thisline[init_len-ii])
info_key = 'coeffs'
elif thisline != [] and info_key == 'coeffs':
exp_list.append((line[11:17]))
for ii in range(init_len,0,-1):
qc.mo_spec[-ii]['coeffs'].append(float(line[16:].split()[init_len-ii]))
elif mo_skip:
mo_skip -= 1
elif sec_flag == 'ecp_info':
if 'THE ECP RUN REMOVES' in line:
sec_flag = None
elif 'PARAMETERS FOR' in line:
if line[17:25].split()[0] != ecp:
ecp = line[17:25].split()[0]
zcore = float(line[51:55].split()[0])
ii_geo = int(line[35:41].split()[0])-1
qc.geo_info[ii_geo][2] = str(float(qc.geo_info[ii_geo][2]) - zcore)
else:
ii_geo = int(line[35:41].split()[0])-1
qc.geo_info[ii_geo][2] = str(float(qc.geo_info[ii_geo][2]) - zcore)
elif sec_flag == 'dm_info':
# instead of giving the output in a useful human and machine readable
# way, gamess output syntax differs for transitions involving the
# ground state compared to transitions between excited states...
if 'GROUND STATE (SCF) DIPOLE=' in line:
# ground state dipole is in debye...convert to atomic units
for ii in range(3):
qc.dipole_moments[0][0][ii] = float(thisline[ii+4])*0.393430307
if 'EXPECTATION VALUE DIPOLE MOMENT FOR EXCITED STATE' in line:
state = (int(line.replace('STATE', 'STATE ').split()[7]))
dm_flag = 'state_info'
if 'TRANSITION FROM THE GROUND STATE TO EXCITED STATE' in line:
state = [0,
int(line.replace('STATE', 'STATE ').split()[8])]
dm_flag = 'transition_info'
if 'TRANSITION BETWEEN EXCITED STATES' in line:
state = [int(thisline[4]),
int(line.replace('AND', 'AND ').split()[6])]
dm_flag = 'transition_info'
if 'NATURAL ORBITAL OCCUPATION NUMBERS FOR EXCITED STATE' in line:
sec_flag = None
dm_flag = None
if dm_flag == 'state_info':
if 'STATE MULTIPLICITY' in line:
qc.states['multiplicity'][state] = int(line.split('=')[1])
if 'STATE ENERGY' in line:
qc.states['energy'][state] = float(line.split('=')[1])
if 'STATE DIPOLE' and 'E*BOHR' in line:
for ii in range(3):
qc.dipole_moments[state][state][ii] = float(thisline[ii+3])
elif dm_flag == 'transition_info':
if 'TRANSITION DIPOLE' and 'E*BOHR' in line:
for ii in range(3):
qc.dipole_moments[state[0]][state[1]][ii] = float(thisline[ii+3])
qc.dipole_moments[state[1]][state[0]][ii] = float(thisline[ii+3])
elif sec_flag == 'pop_info':
if not pop_skip:
if line == '\n':
sec_flag = None
else:
qc.pop_ana['Lowdin'].append(float(thisline[5]))
qc.pop_ana['Mulliken'].append(float(thisline[3]))
elif pop_skip:
pop_skip -= 1
# Check usage of same atomic basis sets
basis_set = {}
for ii in range(len(AO)):
if not AO[ii][0]['atom_type'] in basis_set.keys():
basis_set[AO[ii][0]['atom_type']] = AO[ii]
else:
for jj in range(len(AO[ii])):
if AO[ii][jj]['coeffs'] != basis_set[AO[ii][0]['atom_type']][jj]['coeffs']:
raise IOError('Different basis sets for the same atom.')
# Numpy array
for ii in basis_set.keys():
for jj in range(len(basis_set[ii])):
basis_set[ii][jj]['coeffs'] = numpy.array(basis_set[ii][jj]['coeffs'])
for kk in range(len(qc.mo_spec)):
qc.mo_spec[kk]['coeffs'] = numpy.array(qc.mo_spec[kk]['coeffs'])
# Complement atomic basis sets
for kk in range(len(qc.geo_info)):
for ll in range(len(basis_set[qc.geo_info[kk][0]])):
qc.ao_spec.append({'atom': qc.geo_info[kk][1]-1,
'type': basis_set[qc.geo_info[kk][0]][ll]['type'],
'pnum': basis_set[qc.geo_info[kk][0]][ll]['pnum'],
'coeffs': basis_set[qc.geo_info[kk][0]][ll]['coeffs'],
'exp_list': None
})
# Reconstruct exponents list for ao_spec
count = 0
for i,j in enumerate(qc.ao_spec):
l = l_deg(lquant[j['type']])
j['exp_list'] = []
for i in range(l):
j['exp_list'].append((exp_list[count].lower().count('x'),
exp_list[count].lower().count('y'),
exp_list[count].lower().count('z')))
count += 1
j['exp_list'] = numpy.array(j['exp_list'],dtype=numpy.int64)
if restricted:
for ii in range(len(qc.mo_spec)):
if occ[0] and occ[1]:
qc.mo_spec[ii]['occ_num'] += 2.0
occ[0] -= 1
occ[1] -= 1
if not occ[0] and occ[1]:
qc.mo_spec[ii]['occ_num'] += 1.0
occ[1] -= 1
if not occ[1] and occ[0]:
qc.mo_spec[ii]['occ_num'] += 1.0
occ[0] -= 1
if restricted == False:
for ii in range(len(qc.mo_spec)):
if qc.mo_spec[ii]['spin'] == 'alpha' and occ[0] > 0:
qc.mo_spec[ii]['occ_num'] += 1.0
occ[0] -= 1
has_alpha = True
elif qc.mo_spec[ii]['spin'] == 'beta' and occ[1] > 0:
qc.mo_spec[ii]['occ_num'] += 1.0
occ[1] -= 1
has_beta = True
if spin is not None:
if restricted:
raise IOError('The keyword `spin` is only supported for unrestricted calculations.')
if spin != 'alpha' and spin != 'beta':
raise IOError('`spin=%s` is not a valid option' % spin)
elif spin == 'alpha' and has_alpha == True:
display('Reading only molecular orbitals of spin alpha.')
elif spin == 'beta' and has_beta == True:
display('Reading only molecular orbitals of spin beta.')
elif (not has_alpha) and (not has_beta):
raise IOError(
'No spin molecular orbitals available')
elif ((spin == 'alpha' and not has_alpha) or
(spin == 'beta' and not has_beta)):
raise IOError('You requested `%s` orbitals, but None of them are present.'
% spin)
# Are all MOs requested for the calculation?
if not all_mo:
for i in range(len(qc.mo_spec))[::-1]:
if qc.mo_spec[i]['occ_num'] < 0.0000001:
del qc.mo_spec[i]
# Only molecular orbitals of one spin requested?
if spin is not None:
for i in range(len(qc.mo_spec))[::-1]:
if qc.mo_spec[i]['spin'] != spin:
del qc.mo_spec[i]
# Convert geo_info and geo_spec to numpy.ndarrays
qc.format_geo()
return qc
# read_gamess
[docs]def read_gaussian_fchk(filename, all_mo=False, spin=None, **kwargs):
'''Reads all information desired from a Gaussian FChk file.
**Parameters:**
filename : str
Specifies the filename for the input file.
all_mo : bool, optional
If True, all molecular orbitals are returned.
**Returns:**
qc (class QCinfo) with attributes geo_spec, geo_info, ao_spec, mo_spec, etot :
See :ref:`Central Variables` for details.
'''
fid = open(filename,'r') # Open the file
flines = fid.readlines() # Read the WHOLE file into RAM
fid.close() # Close the file
# Is this an unrestricted calculation?
has_beta = False
is_6D = False
is_10F = False
for line in flines:
if 'beta mo coefficients' in line.lower():
has_beta = True
if 'Pure/Cartesian d shells' in line:
is_6D = int(line.split()[-1]) == 1
if 'Pure/Cartesian f shells' in line:
is_10F = int(line.split()[-1]) == 1
cartesian_basis = (is_6D and is_10F)
if ((not is_6D) and is_10F) or (is_6D and (not is_10F)):
raise IOError('Please apply a Spherical Harmonics (5D, 7F) or '+
'a Cartesian Gaussian Basis Set (6D, 10F)!')
if spin is not None:
if spin != 'alpha' and spin != 'beta':
raise IOError('`spin=%s` is not a valid option' % spin)
elif has_beta:
display('Reading only molecular orbitals of spin %s.' % spin)
else:
raise IOError('The keyword `spin` is only supported for unrestricted calculations.')
restricted = (not has_beta)
sec_flag = None
el_num = [0,0]
mo_i0 = {'alpha': 0, 'beta': 0}
what = 'alpha'
index = 0
at_num = 0
ao_num = 0
ao_sp_coeffs = {}
switch = 0
qc = QCinfo()
qc.geo_info = [[],[],[]]
if not cartesian_basis:
qc.ao_spherical = []
# Go through the file line by line
for il in range(len(flines)):
line = flines[il] # The current line as string
thisline = line.split() # The current line split into segments
# Check the file for keywords
if 'Number of alpha electrons' in line:
el_num[0] = int(thisline[5])
elif 'Number of beta electrons' in line:
el_num[1] = int(thisline[5])
elif 'Number of basis functions' in line:
basis_number = int(thisline[5])
elif 'Atomic numbers' in line:
sec_flag = 'geo_info'
index = 0
at_num = int(thisline[-1])
count = 0
qc.geo_info[1] = list(range(1,at_num+1))
elif 'Nuclear charges' in line:
sec_flag = 'geo_info'
index = 2
at_num = int(thisline[-1])
count = 0
elif 'Total Energy' in line:
qc.etot = float(thisline[3])
elif 'Current cartesian coordinates' in line:
at_num = int(thisline[-1])/3
sec_flag = 'geo_pos'
qc.geo_spec = []
count = 0
xyz = []
elif 'Shell types' in line:
sec_flag = 'ao_info'
index = 'type'
ao_num = int(thisline[-1])
count = 0
if qc.ao_spec == []:
for ii in range(ao_num):
qc.ao_spec.append({})
elif 'Number of primitives per shell' in line:
sec_flag = 'ao_info'
index = 'pnum'
ao_num = int(thisline[-1])
count = 0
if qc.ao_spec == []:
for ii in range(ao_num):
qc.ao_spec.append({})
elif 'Shell to atom map' in line:
sec_flag = 'ao_info'
index = 'atom'
ao_num = int(thisline[-1])
count = 0
if qc.ao_spec == []:
for ii in range(ao_num):
qc.ao_spec.append({})
elif 'Primitive exponents' in line:
sec_flag = 'ao_coeffs'
ao_num = int(thisline[-1])
count = 0
switch = 0
index = 0
if qc.ao_spec == []:
raise IOError('Shell types need to be defined before the AO exponents!')
if not 'coeffs' in qc.ao_spec[0].keys():
for ii in range(len(qc.ao_spec)):
pnum = qc.ao_spec[ii]['pnum']
qc.ao_spec[ii]['coeffs'] = numpy.zeros((pnum, 2))
elif 'Contraction coefficients' in line:
if 'P(S=P)' not in line:
sec_flag = 'ao_coeffs'
else:
sec_flag = 'ao_sp_coeffs'
ao_sp_coeffs = {0: []}
ao_num = int(thisline[-1])
count = 0
switch = 1
index = 0
if qc.ao_spec == []:
raise IOError('Shell types need to be defined before the AO exponents!')
if not 'coeffs' in qc.ao_spec[0].keys():
for ii in range(len(qc.ao_spec)):
pnum = qc.ao_spec[ii]['pnum']
qc.ao_spec[ii]['coeffs'] = numpy.zeros((pnum, 2))
elif 'Orbital Energies' in line:
sec_flag = 'mo_eorb'
mo_num = int(thisline[-1])
mo_i0[thisline[0].lower()] = len(qc.mo_spec)
if restricted:
if el_num[0] == el_num[1]:
i = el_num[0]
occ = 2
else:
i = el_num[0 if 'Alpha' in line else 1]
occ = 1
else:
i = el_num[0 if 'Alpha' in line else 1]
occ = 1
for ii in range(mo_num):
qc.mo_spec.append({'coeffs': numpy.zeros(basis_number),
'energy': 0.0,
'occ_num': float(occ if ii < i else 0),
'sym': '%i.1' % (ii+1),
'spin':thisline[0].lower()
})
elif 'MO coefficients' in line:
sec_flag = 'mo_coeffs'
count = 0
index = 0
mo_num = int(thisline[-1])
what = thisline[0].lower()
else:
# Check if we are in a specific section
if sec_flag == 'geo_info':
for ii in thisline:
qc.geo_info[index].append(ii)
count += 1
if count == at_num:
sec_flag = None
elif sec_flag == 'geo_pos':
for ii in thisline:
xyz.append(float(ii))
if len(xyz) == 3:
qc.geo_spec.append(xyz)
xyz = []
count += 1
if count == at_num:
sec_flag = None
elif sec_flag == 'ao_info':
for ii in thisline:
ii = int(ii)
if index is 'type':
ii = orbit[abs(ii)]
l = lquant[ii]
if not cartesian_basis:
for m in (range(0,l+1) if l != 1 else [1,0]):
qc.ao_spherical.append([count,(l,m)])
if m != 0:
qc.ao_spherical.append([count,(l,-m)])
elif index is 'atom':
ii -= 1
qc.ao_spec[count][index] = ii
count += 1
if count == ao_num:
sec_flag = None
elif sec_flag == 'ao_coeffs':
for ii in thisline:
qc.ao_spec[index]['coeffs'][count,switch] = float(ii)
count += 1
ao_num -= 1
if count == qc.ao_spec[index]['pnum']:
index += 1
count = 0
if not ao_num:
sec_flag = None
elif sec_flag == 'ao_sp_coeffs':
for ii in thisline:
ao_sp_coeffs[index].append(float(ii))
count += 1
ao_num -= 1
if count == qc.ao_spec[index]['pnum']:
index += 1
ao_sp_coeffs[index] = []
count = 0
if not ao_num:
sec_flag = None
elif sec_flag == 'mo_eorb':
for ii in thisline:
qc.mo_spec[count]['energy'] = float(ii)
count += 1
if index != 0 and not count % basis_number:
sec_flag = None
elif sec_flag == 'mo_coeffs':
for ii in thisline:
qc.mo_spec[mo_i0[what]+index]['coeffs'][count] = float(ii)
count += 1
if count == basis_number:
count = 0
index += 1
if index != 0 and not index % basis_number:
sec_flag = None
# Look for SP atomic orbitals
if ao_sp_coeffs:
ao_new = []
ao_spherical_new = []
count = 0
for i,ao in enumerate(qc.ao_spec):
if ao['type'] == 'p' and sum(numpy.abs(ao_sp_coeffs[i])) > 0:
ao_new.append(copy.deepcopy(ao))
ao_new[-1]['type'] = 's'
ao_new.append(ao)
ao_new[-1]['type'] = 'p'
ao_new[-1]['coeffs'][:,1] = numpy.array(ao_sp_coeffs[i])
else:
ao_new.append(ao)
if not cartesian_basis:
ao_spherical_new = []
count = 0
for i,ao in enumerate(qc.ao_spec):
if ao['type'] == 'p' and sum(numpy.abs(ao_sp_coeffs[i])) > 0:
ao_spherical_new.append([count,(0,0)])
count += 1
for m in (1, -1, 0):
ao_spherical_new.append([count,(1,m)])
count += 1
else:
l = lquant[ao['type']]
for m in (range(0,l+1) if l != 1 else [1,0]):
ao_spherical_new.append([count,(l,m)])
if m != 0:
ao_spherical_new.append([count,(l,-m)])
count += 1
assert count == len(ao_new)
qc.ao_spherical = ao_spherical_new
qc.ao_spec = ao_new
# Are all MOs requested for the calculation?
if not all_mo:
for i in range(len(qc.mo_spec))[::-1]:
if qc.mo_spec[i]['occ_num'] < 0.0000001:
del qc.mo_spec[i]
# Only molecular orbitals of one spin requested?
if spin is not None:
for i in range(len(qc.mo_spec))[::-1]:
if qc.mo_spec[i]['spin'] != spin:
del qc.mo_spec[i]
if restricted:
# Closed shell calculation
for mo in qc.mo_spec:
del mo['spin']
else:
# Rename MOs according to spin
for mo in qc.mo_spec:
mo['sym'] += '_%s' % mo['spin'][0]
# Check for natural orbital occupations
energy_sum = sum([abs(i['energy']) for i in qc.mo_spec])
if energy_sum < 0.0000001:
display('Attention!\n\tThis FChk file contains natural orbitals. '+
'(There are no energy eigenvalues.)\n\t' +
'In this case, Gaussian does not print the respective natural' +
'occupation numbers!' )
qc.geo_info = numpy.array(qc.geo_info).T
# Convert geo_info and geo_spec to numpy.ndarrays
qc.format_geo()
return qc
# read_gaussian_fchk
[docs]def read_gaussian_log(filename,all_mo=False,spin=None,orientation='standard',
i_link=-1,i_geo=-1,i_ao=-1,i_mo=-1,interactive=True,
**kwargs):
'''Reads all information desired from a Gaussian .log file.
**Parameters:**
filename : str
Specifies the filename for the input file.
all_mo : bool, optional
If True, all molecular orbitals are returned.
spin : {None, 'alpha', or 'beta'}, optional
If not None, returns exclusively 'alpha' or 'beta' molecular orbitals.
orientation : string, choices={'input', 'standard'}, optional
Specifies orientation of the molecule in Gaussian nomenclature. [#first]_
i_link : int, default=-1
Selects the file for linked Gaussian jobs.
i_geo : int, default=-1
Selects the geometry section of the output file.
i_ao : int, default=-1
Selects the atomic orbital section of the output file.
i_mo : int, default=-1
Selects the molecular orbital section of the output file.
interactive : bool
If True, the user is asked to select the different sets.
**Returns:**
qc (class QCinfo) with attributes geo_spec, geo_info, ao_spec, ao_spherical, mo_spec, etot :
See :ref:`Central Variables` for details.
.. [#first] Attention: The MOs in the output are only valid for the standard orientation!
'''
fid = open(filename,'r') # Open the file
flines = fid.readlines() # Read the WHOLE file into RAM
fid.close() # Close the file
# Search the file the specific sections
count = {'link': 0, 'geometry': 0, 'geometry_input': 0, 'atomic orbitals': 0,
'molecular orbitals': [], 'state': []}
def check_sel(count,i,interactive=False,default=-1):
if count == 0:
raise IndexError
elif count == 1:
return 0
message = '\tPlease give an integer from 0 to {0} (default: {0}): '.format(count-1)
try:
if interactive:
i = raw_input(message)
i = default if i == '' else int(i)
i = range(count)[i]
except (IndexError,ValueError):
raise IOError(message.replace(':','!'))
else:
display('\tSelecting the %s' %
('last element.' if (i == count-1) else 'element %d.' % i))
return i
# Go through the file line by line
for il in range(len(flines)):
line = flines[il] # The current line as string
# Check the file for keywords
if ' Entering Link 1' in line:
count['link'] += 1
try:
display('\tFound %d linked GAUSSIAN files.' % count['link'])
i_link = check_sel(count['link'],i_link,interactive=interactive)
except IndexError:
raise IOError('Found no `Entering Link 1` keyword!')
cartesian_basis = True
c_link = 0
# Go through the file line by line
for il in range(len(flines)):
line = flines[il] # The current line as string
thisline = line.split() # The current line split into segments
# Check the file for keywords
if ' Entering Link 1' in line:
c_link += 1
if i_link == (c_link-1):
if ' orientation:' in line:
if '%s orientation:' % orientation in line.lower():
count['geometry'] += 1
if 'input orientation:' in line.lower():
count['geometry_input'] += 1
elif 'Standard basis:' in line or 'General basis read from cards:' in line:
# Check if a cartesian basis has been applied
if '(5D, 7F)' in line:
cartesian_basis = False
elif '(6D, 10F)' not in line:
raise IOError('Please apply a Spherical Harmonics (5D, 7F) or '+
'a Cartesian Gaussian Basis Set (6D, 10F)!')
elif 'AO basis set' in line:
count['atomic orbitals'] += 1
elif 'The electronic state is ' in line:
count['state'].append(thisline[-1][:-1])
elif 'Orbital Coefficients:' in line:
mo_type = thisline[0]
if mo_type != 'Beta':
count['molecular orbitals'].append(mo_type)
else:
count['molecular orbitals'][-1] = 'Alpha&Beta'
display('\nContent of the GAUSSIAN .log file:')
display('\tFound %d geometry section(s). (%s orientation)' %
(count['geometry'], orientation))
try:
i_geo = check_sel(count['geometry'],i_geo,interactive=interactive)
except IndexError:
count['geometry'] = count['geometry_input']
orientation = 'input'
display('\Looking for "Input orientation": \n' +
'\tFound %d geometry section(s). (%s orientation)' %
(count['geometry'], orientation))
try:
i_geo = check_sel(count['geometry'],i_geo,interactive=interactive)
except IndexError:
raise IOError('Found no geometry section!'+
' Are you sure this is a GAUSSIAN .log file?')
try:
display('\tFound %d atomic orbitals section(s) %s.' %
(count['atomic orbitals'],
'(6D, 10F)' if cartesian_basis else '(5D, 7F)'))
i_ao = check_sel(count['atomic orbitals'],i_ao,interactive=interactive)
except IndexError:
raise IOError('Write GFINPUT in your GAUSSIAN route section to print' +
' the basis set information!')
try:
display('\tFound the following %d molecular orbitals section(s):' %
len(count['molecular orbitals']))
except IndexError:
raise IOError('Write IOP(6/7=3) in your GAUSSIAN route section to print\n' +
' all molecular orbitals!')
for i,j in enumerate(count['molecular orbitals']):
string = '\t\tSection %d: %s Orbitals'% (i,j)
try:
string += ' (electronic state: %s)' % count['state'][i]
except IndexError:
pass
display(string)
i_mo = check_sel(len(count['molecular orbitals']),i_mo,interactive=interactive)
if spin is not None:
if spin != 'alpha' and spin != 'beta':
raise IOError('`spin=%s` is not a valid option' % spin)
else:
display('Reading only molecular orbitals of spin %s.' % spin)
aa_to_au = 1/0.52917720859 # conversion factor for Angstroem to bohr radii
# Set a counter for the AOs
basis_count = 0
# Initialize some variables
sec_flag = None
skip = 0
c_link = 0
c_geo = 0
c_ao = 0
c_mo = 0
c_sao = 0
old_ao = -1
orb_sym = []
qc = QCinfo()
index = []
# Go through the file line by line
for il in range(len(flines)):
line = flines[il] # The current line as string
thisline = line.split() # The current line split into segments
# Check the file for keywords
if ' Entering Link 1' in line:
c_link += 1
if i_link == (c_link-1):
if '%s orientation:' % orientation in line.lower():
# The section containing information about
# the molecular geometry begins
if i_geo == c_geo:
qc.geo_info = []
qc.geo_spec = []
sec_flag = 'geo_info'
c_geo += 1
skip = 4
elif 'Standard basis:' in line or 'General basis read from cards:' in line:
# Check if a cartesian basis has been applied
if '(5D, 7F)' in line:
cartesian_basis = False
elif '(6D, 10F)' not in line:
raise IOError('Please apply a Spherical Harmonics (5D, 7F) or '+
'a Cartesian Gaussian Basis Sets (6D, 10F)!')
elif 'AO basis set' in line:
# The section containing information about
# the atomic orbitals begins
if i_ao == c_ao:
qc.ao_spec = []
if not cartesian_basis:
qc.ao_spherical = []
sec_flag = 'ao_info'
c_ao += 1
basis_count = 0
bNew = True # Indication for start of new AO section
elif 'Orbital symmetries:' in line:
sec_flag = 'mo_sym'
add = ''
orb_sym = []
elif 'Orbital Coefficients:' in line:
# The section containing information about
# the molecular orbitals begins
if (i_mo == c_mo):
sec_flag = 'mo_info'
mo_type = count['molecular orbitals'][i_mo]
qc.mo_spec = []
offset = 0
add = ''
orb_spin = []
if orb_sym == []:
if 'Alpha' in mo_type:
add = '_a'
orb_spin = ['alpha'] * basis_count
orb_sym = ['A1'+add] * basis_count
if 'Beta' in mo_type:
add = '_b'
orb_spin += ['beta'] * basis_count
orb_sym += ['A1'+add] * basis_count
for i in range(len(orb_sym)):
# for numpy version < 1.6
c = ((numpy.array(orb_sym[:i+1]) == orb_sym[i]) != 0).sum()
# for numpy version >= 1.6 this could be used:
#c = numpy.count_nonzero(numpy.array(orb_sym[:i+1]) == orb_sym[i])
qc.mo_spec.append({'coeffs': numpy.zeros(basis_count),
'energy': 0.,
'sym': '%d.%s' % (c,orb_sym[i])})
if orb_spin != []:
qc.mo_spec[-1]['spin'] = orb_spin[i]
if mo_type != 'Beta':
c_mo += 1
bNew = True # Indication for start of new MO section
elif 'E(' in line:
qc.etot = float(line.split('=')[1].split()[0])
else:
# Check if we are in a specific section
if sec_flag == 'geo_info':
if not skip:
qc.geo_info.append([thisline[1],thisline[0],thisline[1]])
qc.geo_spec.append([aa_to_au*float(ij) for ij in thisline[3:]])
if '-----------' in flines[il+1]:
sec_flag = None
else:
skip -= 1
if sec_flag == 'ao_info':
# Atomic orbital section
if ' ****' in line:
# There is a line with stars after every AO
bNew = True
# If there is an additional blank line, the AO section is complete
if flines[il+1].split() == []:
sec_flag = None
elif bNew:
# The following AOs are for which atom?
bNew = False
at_num = int(thisline[0]) - 1
ao_num = 0
elif len(thisline) == 4:
# AO information section
# Initialize a new dict for this AO
ao_num = 0 # Initialize number of atomic orbiatls
ao_type = thisline[0].lower() # Type of atomic orbital
pnum = int(thisline[1]) # Number of primatives
for i_ao in ao_type:
# Calculate the degeneracy of this AO and increase basis_count
basis_count += l_deg(lquant[i_ao],cartesian_basis=cartesian_basis)
qc.ao_spec.append({'atom': at_num,
'type': i_ao,
'pnum': pnum,
'coeffs': numpy.zeros((pnum, 2))
})
else:
# Append the AO coefficients
coeffs = numpy.array(line.replace('D','e').split(), dtype=numpy.float64)
for i_ao in range(len(ao_type)):
qc.ao_spec[-len(ao_type)+i_ao]['coeffs'][ao_num,:] = [coeffs[0],
coeffs[1+i_ao]]
ao_num += 1
if sec_flag == 'mo_sym':
if 'electronic state' in line:
sec_flag = None
else:
info = line[18:].replace('(','').replace(')','').split()
if 'Alpha' in line:
add = '_a'
elif 'Beta' in line:
add = '_b'
for i in info:
orb_sym.append(i + add)
if sec_flag == 'mo_info':
# Molecular orbital section
info = line[:21].split()
if info == []:
coeffs = line[21:].split()
if bNew:
index = [offset+i for i in range(len(coeffs))]
bNew = False
else:
for i,j in enumerate(index):
qc.mo_spec[j]['occ_num'] = int('O' in coeffs[i])
if mo_type not in 'Alpha&Beta':
qc.mo_spec[j]['occ_num'] *= 2
elif 'Eigenvalues' in info:
coeffs = line[21:].replace('-',' -').split()
if mo_type == 'Natural':
key = 'occ_num'
else:
key = 'energy'
for i,j in enumerate(index):
qc.mo_spec[j][key] = float(coeffs[i])
else:
coeffs = line[21:].replace('-',' -').split()
if not cartesian_basis and offset == 0:
if old_ao != line[:14].split()[-1] or len(line[:14].split()) == 4:
old_ao = line[:14].split()[-1]
c_sao += 1
i = c_sao-1
l = lquant[line[13].lower()]
m = line[14:21].replace(' ', '').lower()
p = 'yzx'.find(m) if len(m) == 1 else -1
if p != -1:
m = p - 1
elif m == '':
m = 0
else:
m = int(m)
qc.ao_spherical.append([i,(l,m)])
for i,j in enumerate(index):
qc.mo_spec[j]['coeffs'][int(info[0])-1] = float(coeffs[i])
if int(info[0]) == basis_count:
bNew = True
offset = index[-1]+1
if index[-1]+1 == len(orb_sym):
sec_flag = None
orb_sym = []
# Are all MOs requested for the calculation?
if not all_mo:
for i in range(len(qc.mo_spec))[::-1]:
if qc.mo_spec[i]['occ_num'] < 0.0000001:
del qc.mo_spec[i]
if spin is not None:
if orb_spin == []:
raise IOError('You requested `%s` orbitals, but None of them are present.'
% spin)
else:
for i in range(len(qc.mo_spec))[::-1]:
if qc.mo_spec[i]['spin'] != spin:
del qc.mo_spec[i]
# Convert geo_info and geo_spec to numpy.ndarrays
qc.format_geo()
return qc
# read_gaussian_log
[docs]def spin_check(spin,restricted,has_alpha,has_beta):
'''Check if `spin` keyword is valid.
'''
if spin is not None:
if restricted:
raise IOError('The keyword `spin` is only supported for unrestricted calculations.')
if spin != 'alpha' and spin != 'beta':
raise IOError('`spin=%s` is not a valid option' % spin)
elif spin == 'alpha' and has_alpha:
display('Reading only molecular orbitals of spin alpha.')
elif spin == 'beta' and has_beta:
display('Reading only molecular orbitals of spin beta.')
elif (not has_alpha) and (not has_beta):
raise IOError(
'Molecular orbitals in the input file do not contain `Spin=` keyword')
elif ((spin == 'alpha' and not has_alpha) or
(spin == 'beta' and not has_beta)):
raise IOError('You requested `%s` orbitals, but None of them are present.'
% spin)
[docs]def read_aomix(filename, all_mo=False, spin=None, i_md=-1, interactive=True,
created_by_tmol=True, **kwargs):
'''Reads all information desired from a aomix file.
**Parameters:**
filename : str
Specifies the filename for the input file.
all_mo : bool, optional
If True, all molecular orbitals are returned.
spin : {None, 'alpha', or 'beta'}, optional
If not None, returns exclusively 'alpha' or 'beta' molecular orbitals.
i_md : int, default=-1
Selects the `[AOMix Format]` section of the output file.
interactive : bool
If True, the user is asked to select the different sets.
created_by_tmol : bool
If True and if Cartesian basis set is found, the molecular orbital
coefficients will be converted.
**Returns:**
qc (class QCinfo) with attributes geo_spec, geo_info, ao_spec, mo_spec, etot :
See :ref:`Central Variables` for details.
'''
fid = open(filename,'r') # Open the file
flines = fid.readlines() # Read the WHOLE file into RAM
fid.close() # Close the file
# Is this really a aomix file?
if not '[AOMix Format]\n' in flines:
display('The input file %s is no valid aomix file!\n\nIt does' % filename+
' not contain the keyword: [AOMix Format]\n')
raise IOError('Not a valid input file')
def check_sel(count,i,interactive=False):
if count == 0:
raise IndexError
elif count == 1:
return 0
message = '\tPlease give an integer from 0 to %d: ' % (count-1)
try:
if interactive:
i = int(raw_input(message))
i = range(count)[i]
except (IndexError,ValueError):
raise IOError(message.replace(':','!'))
else:
display('\tSelecting the %s' %
('last element.' if (i == count-1) else 'element %d.' % i))
return i
has_alpha = []
has_beta = []
restricted = []
count = 0
# Go through the file line by line
for il in range(len(flines)):
line = flines[il] # The current line as string
# Check the file for keywords
if '[AOMix Format]' in line:
count += 1
has_alpha.append(False)
has_beta.append(False)
restricted.append(False)
if 'Spin' in line and 'alpha' in line.lower():
has_alpha[-1] = True
if 'Spin' in line and 'beta' in line.lower():
has_beta[-1] = True
if 'Occup' in line:
restricted[-1] = restricted[-1] or (float(line.split('=')[1]) > 1.+1e-4)
if count == 0:
display('The input file %s is no valid aomix file!\n\nIt does' % filename +
' not contain the keyword: [AOMix Format]\n')
raise IOError('Not a valid input file')
else:
if count > 1:
display('\nContent of the aomix file:')
display('\tFound %d [AOMix Format] keywords, i.e., ' % count +
'this file contains %d aomix files.' % count)
i_md = check_sel(count,i_md,interactive=interactive)
spin_check(spin,restricted[i_md],has_alpha[i_md],has_beta[i_md])
# Set a counter for the AOs
basis_count = 0
# Declare synonyms for molden keywords
synonyms = {'Sym': 'sym',
'Ene': 'energy',
'Occup': 'occ_num',
'Spin': 'spin'
}
MO_keys = synonyms.keys()
exp_list = []
count = 0
start_reading = False
# Go through the file line by line
for il in range(len(flines)):
line = flines[il] # The current line as string
thisline = line.split() # The current line split into segments
# Check the file for keywords
if '[aomix format]' in line.lower():
# A new file begins
# Initialize the variables
if i_md == count:
qc = QCinfo()
sec_flag = False # A Flag specifying the current section
start_reading = True # Found the selected section
else:
start_reading = False
count += 1
continue
if start_reading:
if '[SCF Energy / Hartree]' in line:
try:
qc.etot = float(flines[il+1].split()[0])
except IndexError:
pass
elif '[atoms]' in line.lower():
# The section containing information about
# the molecular geometry begins
sec_flag = 'geo_info'
if 'Angs' in line:
# The length are given in Angstroem
# and have to be converted to Bohr radii --
aa_to_au = 1/0.52917720859
else:
# The length are given in Bohr radii
aa_to_au = 1.0
elif '[gto]' in line.lower():
# The section containing information about
# the atomic orbitals begins
sec_flag = 'ao_info'
bNew = True # Indication for start of new AO section
elif '[mo]' in line.lower():
# The section containing information about
# the molecular orbitals begins
sec_flag = 'mo_info'
bNew = True # Indication for start of new MO section
elif '[sto]' in line.lower():
# The orbkit does not support Slater type orbitals
display('orbkit does not work for STOs!\nEXIT\n');
raise IOError('Not a valid input file')
else:
# Check if we are in a specific section
if sec_flag == 'geo_info':
# Geometry section
qc.geo_info.append(thisline[0:3])
qc.geo_spec.append([float(ii)*aa_to_au for ii in thisline[3:]])
if sec_flag == 'ao_info':
# Atomic orbital section
def check_int(i):
try:
int(i)
return True
except ValueError:
return False
if thisline == []:
# There is a blank line after every AO
bNew = True
elif bNew:
# The following AOs are for which atom?
bNew = False
at_num = int(thisline[0]) - 1
ao_num = 0
elif len(thisline) == 3 and check_int(thisline[1]):
# AO information section
# Initialize a new dict for this AO
ao_num = 0 # Initialize number of atomic orbiatls
ao_type = thisline[0] # Which type of atomic orbital do we have
pnum = int(thisline[1]) # Number of primatives
# Calculate the degeneracy of this AO and increase basis_count
for i_ao in ao_type:
# Calculate the degeneracy of this AO and increase basis_count
basis_count += l_deg(lquant[i_ao])
qc.ao_spec.append({'atom': at_num,
'type': i_ao,
'pnum': pnum,
'coeffs': numpy.zeros((pnum, 2))
})
else:
# Append the AO coefficients
coeffs = numpy.array(line.replace('D','e').split(), dtype=numpy.float64)
for i_ao in range(len(ao_type)):
qc.ao_spec[-len(ao_type)+i_ao]['coeffs'][ao_num,:] = [coeffs[0],
coeffs[1+i_ao]]
ao_num += 1
if sec_flag == 'mo_info':
# Molecular orbital section
if '=' in line:
# MO information section
if bNew:
# Create a numpy array for the MO coefficients and
# for backward compability create a simple counter for 'sym'
qc.mo_spec.append({'coeffs': numpy.zeros(basis_count),
'sym': '%d.1' % (len(qc.mo_spec)+1)})
bNew = False
# Append information to dict of this MO
info = line.replace('\n','').replace(' ','')
info = info.split('=')
if info[0] in MO_keys:
if info[0] == 'Spin':
info[1] = info[1].lower()
elif info[0] != 'Sym':
info[1] = float(info[1])
elif not '.' in info[1]:
from re import search
a = search(r'\d+', info[1]).group()
if a == info[1]:
info[1] = '%s.1' % a
else:
info[1] = info[1].replace(a, '%s.' % a, 1)
qc.mo_spec[-1][synonyms[info[0]]] = info[1]
else:
if ('[' or ']') in line:
# start of another section that is not (yet) read
sec_flag = None
else:
# Append the MO coefficients
bNew = True # Reset bNew
index = int(thisline[0])-1
try:
# Try to convert coefficient to float
qc.mo_spec[-1]['coeffs'][index] = float(thisline[-1])
if len(qc.mo_spec) == 1:
exp_list.append(thisline[-2])
except ValueError:
# If it cannot be converted print error message
display('Error in coefficient %d of MO %s!' % (index,
qc.mo_spec[-1]['sym']) + '\nSetting this coefficient to zero...')
# Check usage of same atomic basis sets
for ii in range(len(exp_list)):
s = exp_list[ii]
exp = [0,0,0]
c_last = None
for jj in s[1:]:
try:
c = int(jj)
exp[c_last] += (c-1)
except ValueError:
for kk,ll in enumerate('xyz'):
if jj == ll:
exp[kk] += 1
c_last = kk
exp_list[ii] = exp
count = 0
for i,j in enumerate(qc.ao_spec):
l = l_deg(lquant[j['type']])
j['exp_list'] = []
for i in range(l):
j['exp_list'].append((exp_list[count][0],
exp_list[count][1],
exp_list[count][2]))
count += 1
j['exp_list'] = numpy.array(j['exp_list'],dtype=numpy.int64)
# For Cartesian basis sets in Turbomole, the molecular orbital coefficients
# have to be converted.
is_tmol_cart = not (len(qc.mo_spec) % len(qc.mo_spec[0]['coeffs']))
# Are all MOs requested for the calculation?
if not all_mo:
for i in range(len(qc.mo_spec))[::-1]:
if qc.mo_spec[i]['occ_num'] < 0.0000001:
del qc.mo_spec[i]
# Modify qc.mo_spec to support spin
qc.select_spin(restricted[i_md],spin=spin)
# Convert geo_info and geo_spec to numpy.ndarrays
qc.format_geo()
if is_tmol_cart and created_by_tmol:
display('\nFound a Cartesian basis set in the AOMix file.')
display('We assume that this file has been created by Turbomole.')
display('Applying a conversion to the molecular orbital coefficients, ')
display('in order to get normalized orbitals.')
# Convert MO coefficients
from orbkit.analytical_integrals import create_mo_coeff, get_lxlylz
def dfact(n):
if n <= 0:
return 1
else:
return n * dfact(n-2)
mo = create_mo_coeff(qc.mo_spec)
for i,j in enumerate(get_lxlylz(qc.ao_spec)):
norm = (dfact(2*j[0] - 1) * dfact(2*j[1] - 1) * dfact(2*j[2] - 1))
j = sum(j)
if j >1:
mo[:,i] *= numpy.sqrt(norm)
for ii in range(len(qc.mo_spec)):
qc.mo_spec[ii]['coeffs'] = mo[ii]
return qc
# read_aomix
[docs]def read_wfx(filename, all_mo=False, spin=None, **kwargs):
'''Reads all information desired from a wfn file.
**Parameters:**
filename : str
Specifies the filename for the input file.
all_mo : bool, optional
If True, all molecular orbitals are returned.
spin : {None, 'alpha', or 'beta'}, optional
If not None, returns exclusively 'alpha' or 'beta' molecular orbitals.
**Returns:**
qc (class QCinfo) with attributes geo_spec, geo_info, ao_spec, mo_spec, etot :
See :ref:`Central Variables` for details.
'''
# Initialize the variables
qc = QCinfo()
exp_list = []
for j in exp_wfn:
exp_list.extend(j)
exp_list = numpy.array(exp_list,dtype=numpy.int64)
fid = open(filename,'r') # Open the file
flines = fid.readlines() # Read the WHOLE file into RAM
fid.close() # Close the file
is_valid = False
for il in range(len(flines)):
if '<Keywords>' in flines[il] and 'GTO' in flines[il+1]:
is_valid = True
if not is_valid: raise IOError('No valid .wfx file!\nMissing:\n' +
'<Keywords>\n GTO\n</Keywords>')
sec_flag = None # A Flag specifying the current section
at_num = None
mo_num = None
ao_num = None
restricted = True
count = 0
# Go through the file line by line
for il in range(len(flines)):
line = flines[il] # The current line as string
if '<Number of Nuclei>' in line:
at_num = int(flines[il+1])
qc.geo_info = [[None, i+1, None] for i in range(at_num)]
qc.geo_spec = []
elif '<Nuclear Names>' in line:
if not at_num: raise IOError('`<Number of Nuclei>` has to be found ' +
'before `<Nuclear Names>`.')
for i in range(at_num):
qc.geo_info[i][0] = flines[il+i+1].replace(' ','').replace('\n','')
elif '<Atomic Numbers>' in line:
if not at_num: raise IOError('`<Number of Nuclei>` has to be found ' +
'before `<Atomic Numbers>`.')
for i in range(at_num):
qc.geo_info[i][2] = flines[il+i+1].replace(' ','').replace('\n','')
elif '<Nuclear Cartesian Coordinates>' in line:
if not at_num: raise IOError('`<Number of Nuclei>` has to be found ' +
'before `<Nuclear Cartesian Coordinates>`.')
for i in range(at_num):
qc.geo_spec.append(flines[il+i+1].split())
elif '<Number of Primitives>' in line:
ao_num = int(flines[il+1])
qc.ao_spec = [{'atom': None,
'pnum': -1,
'coeffs': None,
'exp_list': None,
} for i in range(ao_num)]
elif '<Primitive Centers>' in line:
sec_flag = 'ao_center'
count = 0
elif '<Primitive Types>' in line:
sec_flag = 'ao_type'
count = 0
elif '<Primitive Exponents>' in line:
sec_flag = 'ao_exp'
count = 0
elif '<Number of Occupied Molecular Orbitals>' in line:
mo_num = int(flines[il+1])
qc.mo_spec = [{'coeffs': numpy.zeros(ao_num),
'energy': None,
'occ_num': None,
'spin': None,
'sym': '%s.1' % (i+1)
} for i in range(mo_num)]
elif '<Molecular Orbital Occupation Numbers>' in line:
for i in range(mo_num):
qc.mo_spec[i]['occ_num'] = float(flines[il+1+i])
elif '<Molecular Orbital Energies>' in line:
for i in range(mo_num):
qc.mo_spec[i]['energy'] = float(flines[il+1+i])
elif '<Molecular Orbital Spin Types>' in line:
for i in range(mo_num):
qc.mo_spec[i]['spin'] = (flines[il+1+i].replace(' ','').replace('\n','')
).replace('and','_').lower()
restricted = restricted and ('_' in qc.mo_spec[i]['spin'])
elif '<MO Number>' in line:
index = int(flines[il+1])-1
for i in range(ao_num):
qc.mo_spec[index]['coeffs'][i] = float(flines[il+3+i])
elif '</' in line:
sec_flag = None
elif sec_flag is not None:
if sec_flag == 'ao_center':
for i in line.split():
qc.ao_spec[count]['atom'] = int(i)-1
count += 1
if sec_flag == 'ao_type':
for i in line.split():
qc.ao_spec[count]['exp_list'] = exp_list[int(i)-1][numpy.newaxis]
count += 1
if sec_flag == 'ao_exp':
for i in line.split():
qc.ao_spec[count]['coeffs'] = numpy.array([[float(i),1.0]])
count += 1
has_alpha = any([i['spin'] == 'alpha' for i in qc.mo_spec])
has_beta = any([i['spin'] == 'beta' for i in qc.mo_spec])
spin_check(spin,restricted,has_alpha,has_beta)
qc.select_spin(restricted,spin=spin)
# Remove numbers from atom names
for i in qc.geo_info:
i[0] = ''.join([k for k in i[0] if not k.isdigit()])
# Convert geo_info and geo_spec to numpy.ndarrays
qc.format_geo()
return qc
[docs]def read_wfn(filename, all_mo=False, spin=None, **kwargs):
'''Reads all information desired from a wfn file.
**Parameters:**
filename : str
Specifies the filename for the input file.
all_mo : bool, optional
If True, all molecular orbitals are returned.
**Returns:**
qc (class QCinfo) with attributes geo_spec, geo_info, ao_spec, mo_spec, etot :
See :ref:`Central Variables` for details.
'''
if spin is not None:
raise IOError('The option `spin` is not supported for the `.wfn` reader.')
# Initialize the variables
qc = QCinfo()
sec_flag = None # A Flag specifying the current section
is_wfn = False # Check type of file
ao_num = 0 # Number of AO
mo_num = 0 # Number of MO
at_num = 0 # Number of atoms
c_type = 0 # Counting variable for AO type
c_exp = 0 # Counting variable for AO exponents
exp_list = []
for j in exp_wfn:
exp_list.extend(j)
exp_list = numpy.array(exp_list,dtype=numpy.int64)
with open(filename) as fileobject:
for line in fileobject:
thisline = line.split() # The current line split into segments
# Check the file for keywords
if 'GAUSSIAN' in line or 'GTO' in line:
if len(thisline) == 8:
mo_num = int(thisline[1])
ao_num = int(thisline[4])
at_num = int(thisline[6])
sec_flag = 'geo_info'
elif 'CENTRE ASSIGNMENTS' in line:
thisline = line[20:].split()
for i in range(len(thisline)):
qc.ao_spec.append({'atom': int(thisline[i])-1,
'pnum': -1,
'coeffs': None,
'exp_list': None,
})
elif 'TYPE ASSIGNMENTS' in line:
thisline = line[18:].split()
for i in range(len(thisline)):
qc.ao_spec[c_type]['exp_list'] = exp_list[int(thisline[i])-1][numpy.newaxis]
c_type += 1
elif 'EXPONENTS' in line:
thisline = line.replace('EXPONENTS','').replace('D','E').split()
for i in thisline:
qc.ao_spec[c_exp]['coeffs'] = numpy.array([[float(i),1.0]])
c_exp += 1
elif 'MO' in line and 'OCC NO =' in line and 'ORB. ENERGY =' in line:
qc.mo_spec.append({'coeffs': numpy.zeros(ao_num),
'energy': float(line[25:].split()[7]),
'occ_num': float(line[25:].split()[3]),
'sym': '%s.1' % thisline[1]
})
sec_flag = 'mo_info'
c_mo = 0 # Counting variable for MOs
else:
if sec_flag == 'geo_info':
if not at_num:
sec_flag = None
elif at_num:
qc.geo_info.append([thisline[0],thisline[-7][:-1],thisline[-1]])
qc.geo_spec.append([float(ii) for ii in thisline[-6:-3]])
at_num -= 1
elif sec_flag == 'mo_info':
for i in thisline:
if (c_mo) < ao_num:
qc.mo_spec[-1]['coeffs'][c_mo] = numpy.array(
float(i.replace('D','E')))
c_mo += 1
if (c_mo) == ao_num:
sec_flag = None
# Remove numbers from atom names
for i in qc.geo_info:
i[0] = ''.join([k for k in i[0] if not k.isdigit()])
# Convert geo_info and geo_spec to numpy.ndarrays
qc.format_geo()
return qc
[docs]def read_with_cclib(filename, cclib_parser=None, all_mo=False, spin=None,
**kwargs):
'''Reads all information desired using cclib.
**Parameters:**
filename : str
Specifies the filename for the input file.
cclib_parser : str
If itype is 'cclib', specifies the cclib.parser.
all_mo : bool, optional
If True, all molecular orbitals are returned.
spin : {None, 'alpha', or 'beta'}, optional
If not None, returns exclusively 'alpha' or 'beta' molecular orbitals.
**Returns:**
qc (class QCinfo) with attributes geo_spec, geo_info, ao_spec, mo_spec, etot :
See :ref:`Central Variables` for details.
'''
if not isinstance(cclib_parser,str):
raise IOError('cclib requires the specification of parser, e.g., ' +
'cclib_parser="Gaussian".')
if cclib_parser == 'Molpro':
display('\nThe Molpro basis set is not properly read by the cclib parser.')
display('Please create a molden file with Molpro, i.e., ' +
'\n\tput,molden,output.molden,NEW;\n')
exec('from cclib.parser import %s as parser' % cclib_parser)
ccData = parser(filename).parse()
return convert_cclib(ccData, all_mo=all_mo, spin=spin)
[docs]def convert_cclib(ccData, all_mo=False, spin=None):
'''Converts a ccData class created by cclib to an instance of
orbkit's QCinfo class.
**Parameters:**
ccData : class
Contains the input data created by cclib.
all_mo : bool, optional
If True, all molecular orbitals are returned.
spin : {None, 'alpha', or 'beta'}, optional
If not None, returns exclusively 'alpha' or 'beta' molecular orbitals.
**Returns:**
qc (class QCinfo) with attributes geo_spec, geo_info, ao_spec, mo_spec, etot :
See :ref:`Central Variables` for details.
'''
aa_to_au = 1/0.52917720859
# Initialize the variables
qc = QCinfo()
# Converting all information concerning atoms and geometry
qc.geo_spec = ccData.atomcoords[0] * aa_to_au
for ii in range(ccData.natom):
symbol = get_atom_symbol(atom=ccData.atomnos[ii])
qc.geo_info.append([symbol,str(ii+1),str(ccData.atomnos[ii])])
# Convert geo_info and geo_spec to numpy.ndarrays
qc.format_geo()
# Converting all information about atomic basis set
for ii in range(ccData.natom):
for jj in range(len(ccData.gbasis[ii])):
pnum = len(ccData.gbasis[ii][jj][1])
qc.ao_spec.append({'atom': ii,
'type': str(ccData.gbasis[ii][jj][0]).lower(),
'pnum': pnum,
'coeffs': numpy.zeros((pnum, 2))
})
for kk in range(pnum):
qc.ao_spec[-1]['coeffs'][kk][0] = ccData.gbasis[ii][jj][1][kk][0]
qc.ao_spec[-1]['coeffs'][kk][1] = ccData.gbasis[ii][jj][1][kk][1]
if hasattr(ccData,'aonames'):
# Reconstruct exponents list for ao_spec
cartesian_basis = True
for i in ccData.aonames:
if '+' in i or '-' in i:
cartesian_basis = False
if not cartesian_basis:
qc.ao_spherical = []
count = 0
for i,ao in enumerate(qc.ao_spec):
l = l_deg(lquant[ao['type']],cartesian_basis=cartesian_basis)
if cartesian_basis:
ao['exp_list'] = []
for ll in range(l):
if cartesian_basis:
ao['exp_list'].append((ccData.aonames[count].lower().count('x'),
ccData.aonames[count].lower().count('y'),
ccData.aonames[count].lower().count('z')))
else:
m = ccData.aonames[count].lower().split('_')[-1]
m = m.replace('+',' +').replace('-',' -').replace('s','s 0').split(' ')
p = 'yzx'.find(m[0][-1])
if p != -1:
m = p - 1
else:
m = int(m[-1])
qc.ao_spherical.append([i,(lquant[ao['type']],m)])
count += 1
# Converting all information about molecular orbitals
ele_num = numpy.sum(ccData.atomnos) - numpy.sum(ccData.coreelectrons) - ccData.charge
ue = (ccData.mult-1)
# Check for natural orbitals and occupation numbers
is_natorb = False
if hasattr(ccData,'nocoeffs'):
if not hasattr(ccData,'nooccnos'):
raise IOError('There are natural orbital coefficients (`nocoeffs`) in the cclib' +
' ccData, but no natural occupation numbers (`nooccnos`)!')
is_natorb = True
restricted = (len(ccData.mosyms) == 1)
if spin is not None:
if spin != 'alpha' and spin != 'beta':
raise IOError('`spin=%s` is not a valid option' % spin)
elif restricted:
raise IOError('The keyword `spin` is only supported for unrestricted calculations.')
else:
display('Converting only molecular orbitals of spin %s.' % spin)
sym = {}
if len(ccData.mosyms) == 1:
add = ['']
orb_sym = [None]
else:
add = ['_a','_b']
orb_sym = ['alpha','beta']
nmo = ccData.nmo if hasattr(ccData,'nmo') else len(ccData.mocoeffs[0])
for ii in range(nmo):
for i,j in enumerate(add):
a = '%s%s' % (ccData.mosyms[i][ii],j)
if a not in sym.keys(): sym[a] = 1
else: sym[a] += 1
if is_natorb:
occ_num = ccData.nooccnos[ii]
elif not restricted:
occ_num = 1.0 if ii <= ccData.homos[i] else 0.0
elif ele_num > ue:
occ_num = 2.0
ele_num -= 2.0
elif ele_num > 0.0 and ele_num <= ue:
occ_num = 1.0
ele_num -= 1.0
ue -= 1.0
else:
occ_num = 0.0
qc.mo_spec.append({'coeffs': (ccData.nocoeffs if is_natorb else ccData.mocoeffs[i])[ii],
'energy': 0.0 if is_natorb else ccData.moenergies[i][ii],
'occ_num': occ_num,
'sym': '%d.%s' %(sym[a],a)
})
if orb_sym[i] is not None:
qc.mo_spec[-1]['spin'] = orb_sym[i]
if spin is not None and spin != orb_sym[i]:
del qc.mo_spec[-1]
# Use default order for atomic basis functions if aonames is not present
if not hasattr(ccData,'aonames'):
display('The attribute `aonames` is not present in the parsed data.')
display('Using the default order of basis functions.')
# Check which basis functions have been used
c_cart = sum([l_deg(l=ao['type'], cartesian_basis=True) for ao in qc.ao_spec])
c_sph = sum([l_deg(l=ao['type'], cartesian_basis=False) for ao in qc.ao_spec])
c = create_mo_coeff(qc.mo_spec,'').shape[-1]
if c != c_cart and c == c_sph: # Spherical basis
qc.ao_spherical = get_ao_spherical(qc.ao_spec,p=[0,1])
elif c != c_cart:
display('Warning: The basis set type does not match with pure spherical ' +
'or pure Cartesian basis!')
display('Please specify qc.mo_spec["exp_list"] and/or qc.ao_spherical by your self.')
# Are all MOs requested for the calculation?
if not all_mo:
for i in range(len(qc.mo_spec))[::-1]:
if qc.mo_spec[i]['occ_num'] < 0.0000001:
del qc.mo_spec[i]
return qc
def get_ao_spherical(ao_spec,p=[1,0]):
ao_spherical = []
for i,ao in enumerate(ao_spec):
ii = ao['type']
l = lquant[ii]
for m in (range(0,l+1) if l != 1 else p):
ao_spherical.append([i,(l,m)])
if m != 0:
ao_spherical.append([i,(l,-m)])
#for m in (range(1,l+1) if l != 1 else p):
#if m != 0:
#ao_spherical.append([i,(l,-m)])
return ao_spherical
[docs]def mo_select(mo_spec, fid_mo_list, strict=False):
'''Selects molecular orbitals from an external file or a list of molecular
orbital labels.
**Parameters:**
mo_spec :
See :ref:`Central Variables` for details.
strict : bool, optional
If True, orbkit will follow strictly the fid_mo_list, i.e., the order of
the molecular orbitals will be kept and multiple occurrences of items
will evoke multiple calculations of the respective molecular orbitals.
fid_mo_list : str, `'all_mo'`, or list
| If fid_mo_list is a str, specifies the filename of the molecular orbitals list.
| If fid_mo_list is 'all_mo', creates a list containing all molecular orbitals.
| If fid_mo_list is a list, provides a list (or a list of lists) of molecular
orbital labels.
**Supported Formats:**
Integer List (Counting from **ONE**!)::
1 2 3
5 4
homo lumo+2:lumo+4
List with Symmetry Labels::
1.1 2.1 1.3
1.1 4.1
4.1 2.3 2.1
**Returns:**
Dictionary with following Members:
:mo: - List of molecular orbital labels.
:mo_ii: - List of molecular orbital indices.
:mo_spec: - Selected elements of mo_spec. See :ref:`Central Variables` for details.
:mo_in_file: - List of molecular orbital labels within the fid_mo_list file.
:sym_select: - If True, symmetry labels have been used.
..attention:
For **unrestricted** calculations, orbkit adds `_a` (alpha) or `_b` (beta) to
the symmetry labels, e.g., `1.1_a`.
If you have specified the option `spin=alpha` or `spin=beta`, only the
alpha or the beta orbitals are taken into account for the counting
within the Integer List.
'''
import re
display('\nProcessing molecular orbital list...')
mo_in_file = []
selected_mo = []
sym_select = False
def assign_selected_mo(selected_mo,mo_spec,strict=False,
what=lambda x,y: y[x]['sym']):
selected_mo_spec = []
selected_mo_ii = []
for i in selected_mo:
is_present = False
for k in range(len(mo_spec)):
if (what(k,mo_spec) == i):
is_present = True
if strict or (i not in selected_mo_ii):
selected_mo_spec.append(mo_spec[k])
selected_mo_ii.append(what(k,mo_spec))
if not is_present:
raise IOError('Cannot find %s in mo_spec' % i)
selected_mo_ii = numpy.array(selected_mo_ii)
return selected_mo_spec,selected_mo_ii
def expr2int(expr):
if isinstance(expr,int):
return expr
x = 0
for i in re.findall(r'\d+|[+-]\d+',expr):
x += int(i)
return x
def get_selection(selected_mo):
mo_occup = numpy.array([i['occ_num'] for i in mo_spec])
homo = (mo_occup>0.).nonzero()[0][-1] + 1 # molden numbering
lumo = (mo_occup>0.).nonzero()[0][-1]+1 + 1 # molden numbering
mo_energy = numpy.array([i['energy'] for i in mo_spec])
last_bound = sum(mo_energy<=0.0) # molden numbering
sel = []
for i in selected_mo:
i = i.lower().replace('homo',str(homo)).replace('lumo',str(lumo))
i = i.replace('last_bound',str(last_bound))
if ':' in i:
k = [1,len(mo_spec)+1,1]
i = i.split(':')
for ik,j in enumerate(i):
if j != '': k[ik] = j
i = list(range(*[expr2int(j) for j in k]))
sel.extend(i)
else:
sel.append(expr2int(i))
return sel
if isinstance(fid_mo_list,str) and fid_mo_list.lower() == 'all_mo':
selected_mo = numpy.array(numpy.arange(len(mo_spec))+1, dtype=numpy.str)
mo_in_file = [selected_mo]
selected_mo_spec = mo_spec
selected_mo_ii = numpy.array([i['sym'] for i in selected_mo_spec])
else:
if isinstance(fid_mo_list,str) and not os.path.exists(fid_mo_list):
if ',' in fid_mo_list:
fid_mo_list = fid_mo_list.split(',')
else:
fid_mo_list = [fid_mo_list]
if isinstance(fid_mo_list, list):
for i in fid_mo_list:
if not isinstance(i, list):
i = i.split(',') if isinstance(i,str) else [i]
selected_mo.extend(list(map(str,i)))
mo_in_file.append(list(map(str,i)))
else:
try:
fid=open(fid_mo_list,'r')
flines = fid.readlines()
fid.close()
for line in flines:
integer = line.replace(',',' ').split()
mo_in_file.append(integer)
selected_mo.extend(integer)
except:
raise IOError('The selected mo-list (%(m)s) is not valid!' %
{'m': fid_mo_list} + '\ne.g.\n\t1\t3\n\t2\t7\t9\n')
# Print some information
for i,j in enumerate(mo_in_file):
display('\tLine %d: %s' % (i+1,', '.join(j)))
# Check if the molecular orbitals are specified by symmetry
# (e.g. 1.1 in MOLPRO nomenclature) or
# by the number in the input file (e.g. 1)
try: # Try to convert selections into integer
for i in selected_mo:
if isinstance(i,int):
continue
i = i.replace('homo','1').replace('lumo','2').replace('last_bound','3')
for r in ['-','+',':']:
i = i.replace(r,'')
int(i)
except ValueError:
sym_select = True
errors = []
for i in range(len(selected_mo)):
if not '.' in selected_mo[i]:
errors.append(i)
if errors:
err = [selected_mo[i] for i in errors]
raise IOError('`%s` are no valid labels according '% ', '.join(err) +
'to the MOLPRO nomenclature, e.g., `5.1` or `5.A1`.' +
'\n\tHint: You cannot mix integer numbering and MOLPRO\'s ' +
'symmetry labels')
if sym_select:
what = lambda x,y: y[x]['sym']
selected_mo_spec,selected_mo_ii = assign_selected_mo(selected_mo,
mo_spec,
strict=strict,
what=what)
else:
selected_mo = get_selection(selected_mo)
if not strict:
selected_mo = list(map(int, selected_mo))
selected_mo.sort()
selected_mo = list(map(str, selected_mo))
what = lambda x,y: str(x+1)
selected_mo_spec,selected_mo_ii = assign_selected_mo(selected_mo,
mo_spec,
strict=strict,
what=what)
selected_mo = selected_mo_ii
for i in range(len(mo_in_file)):
mo_in_file[i] = list(map(str, get_selection(mo_in_file[i])))
# Print some information
display('\nThe following orbitals will be considered...')
for i,j in enumerate(mo_in_file):
display('\tLine %d: %s' % (i+1,', '.join(j)))
display('')
return {'mo': selected_mo, 'mo_ii': selected_mo_ii,
'mo_spec': selected_mo_spec,
'mo_in_file': mo_in_file, 'sym_select': sym_select}