# -*- coding: iso-8859-1 -*-
'''Module for all additional features of ORBKIT.'''
'''
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 general modules
import os
import numpy
from scipy import integrate
# Import orbkit modules
from orbkit import core,grid,output,options
from orbkit.display import display
from orbkit.read import mo_select
[docs]def calc_mo(qc, fid_mo_list, drv=None, otype=None, ofid=None,
numproc=None, slice_length=None):
'''Calculates and saves the selected molecular orbitals or the derivatives thereof.
**Parameters:**
qc.geo_spec, qc.geo_info, qc.ao_spec, qc.mo_spec :
See :ref:`Central Variables` for details.
fid_mo_list : str
Specifies the filename of the molecular orbitals list or list of molecular
orbital labels (cf. :mod:`orbkit.read.mo_select` for details).
If fid_mo_list is 'all_mo', creates a list containing all molecular orbitals.
drv : string or list of strings {None,'x','y', 'z', 'xx', 'xy', ...}, optional
If not None, a derivative calculation of is requested.
otype : str or list of str, optional
Specifies output file type. See :data:`otypes` for details.
ofid : str, optional
Specifies output file name. If None, the filename will be based on
:mod:`orbkit.options.outputname`.
numproc : int, optional
Specifies number of subprocesses for multiprocessing.
If None, uses the value from :mod:`options.numproc`.
slice_length : int, optional
Specifies the number of points per subprocess.
If None, uses the value from :mod:`options.slice_length`.
**Returns:**
mo_list : numpy.ndarray, shape=((NMO,) + N)
Contains the NMO=len(qc.mo_spec) molecular orbitals on a grid.
mo_info : dict
Contains information of the selected molecular orbitals and has 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.
'''
mo_info = mo_select(qc.mo_spec, fid_mo_list, strict=True)
qc_select = qc.todict()
qc_select['mo_spec'] = mo_info['mo_spec']
slice_length = options.slice_length if slice_length is None else slice_length
numproc = options.numproc if numproc is None else numproc
# Calculate the AOs and MOs
mo_list = core.rho_compute(qc_select,
calc_mo=True,
drv=drv,
slice_length=slice_length,
numproc=numproc)
if otype is None:
return mo_list, mo_info
if ofid is None:
ofid = '%s_MO' % (options.outputname)
if not options.no_output:
if 'h5' in otype:
output.main_output(mo_list,qc.geo_info,qc.geo_spec,data_id='MO',
outputname=ofid,
mo_spec=qc_select['mo_spec'],drv=drv,is_mo_output=True)
# Create Output
cube_files = []
for i,j in enumerate(qc_select['mo_spec']):
outputname = '%s_%s' % (ofid,mo_info['mo'][i])
comments = ('%s,Occ=%.1f,E=%+.4f' % (mo_info['mo'][i],
j['occ_num'],
j['energy']))
index = numpy.index_exp[:,i] if drv is not None else i
output_written = output.main_output(mo_list[index],
qc.geo_info,qc.geo_spec,
outputname=outputname,
comments=comments,
otype=otype,omit=['h5','vmd','mayavi'],
drv=drv)
for i in output_written:
if i.endswith('.cb'):
cube_files.append(i)
if 'vmd' in otype and cube_files != []:
display('\nCreating VMD network file...' +
'\n\t%(o)s.vmd' % {'o': ofid})
output.vmd_network_creator(ofid,cube_files=cube_files)
if 'mayavi' in otype:
datalabels = ['MO %(sym)s, Occ=%(occ_num).2f, E=%(energy)+.4f E_h' %
i for i in qc_select['mo_spec']]
if drv is not None:
tmp = []
for i in drv:
for j in datalabels:
tmp.append('d/d%s of %s' % (i,j))
datalabels = tmp
data = mo_list.reshape((-1,) + grid.get_shape())
output.main_output(data,qc.geo_info,qc.geo_spec,
otype='mayavi',datalabels=datalabels)
return mo_list, mo_info
[docs]def mo_set(qc, fid_mo_list, drv=None, laplacian=None,
otype=None, ofid=None, return_all=True,
numproc=None, slice_length=None):
'''Calculates and saves the density or the derivative thereof
using selected molecular orbitals.
**Parameters:**
qc.geo_spec, qc.geo_info, qc.ao_spec, qc.mo_spec :
See :ref:`Central Variables` for details.
fid_mo_list : str
Specifies the filename of the molecular orbitals list or list of molecular
orbital labels (cf. :mod:`orbkit.read.mo_select` for details).
otype : str or list of str, optional
Specifies output file type. See :data:`otypes` for details.
ofid : str, optional
Specifies output file name. If None, the filename will be based on
:mod:`orbkit.options.outputname`.
drv : string or list of strings {None,'x','y', 'z', 'xx', 'xy', ...}, optional
If not None, a derivative calculation is requested.
return_all : bool
If False, no data will be returned.
numproc : int, optional
Specifies number of subprocesses for multiprocessing.
If None, uses the value from :mod:`options.numproc`.
slice_length : int, optional
Specifies the number of points per subprocess.
If None, uses the value from :mod:`options.slice_length`.
**Returns:**
datasets : numpy.ndarray, shape=((NSET,) + N)
Contains the NSET molecular orbital sets on a grid.
delta_datasets : numpy.ndarray, shape=((NSET,NDRV) + N)
Contains the NDRV NSET molecular orbital set on a grid. This is only
present if derivatives are requested.
mo_info : dict
Contains information of the selected molecular orbitals and has 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.
'''
mo_info = mo_select(qc.mo_spec, fid_mo_list, strict=False)
qc_select = qc.todict()
drv = options.drv if drv is None else drv
laplacian = options.laplacian if laplacian is None else laplacian
slice_length = options.slice_length if slice_length is None else slice_length
numproc = options.numproc if numproc is None else numproc
if ofid is None:
ofid = options.outputname
if 'h5' in otype and os.path.exists('%s.h5' % ofid):
raise IOError('%s.h5 already exists!' % ofid)
datasets = []
delta_datasets = []
cube_files = []
for i_file,j_file in enumerate(mo_info['mo_in_file']):
display('Starting with the %d. element of the molecular orbital list (%s)...\n\t' %
(i_file+1,fid_mo_list) + str(j_file) +
'\n\t(Only regarding existing and occupied mos.)\n')
qc_select['mo_spec'] = []
for i_mo,j_mo in enumerate(mo_info['mo']):
if j_mo in j_file:
if mo_info['sym_select']:
ii_mo = numpy.argwhere(mo_info['mo_ii'] == j_mo)
else:
ii_mo = i_mo
qc_select['mo_spec'].append(mo_info['mo_spec'][int(ii_mo)])
data = core.rho_compute(qc_select,
drv=drv,
laplacian=laplacian,
slice_length=slice_length,
numproc=numproc)
datasets.append(data)
if drv is None:
rho = data
elif laplacian:
rho, delta_rho, laplacian_rho = data
delta_datasets.extend(delta_rho)
delta_datasets.append(laplacian_rho)
else:
rho, delta_rho = data
delta_datasets.append(delta_rho)
if options.z_reduced_density:
if grid.is_vector:
display('\nSo far, reducing the density is not supported for vector grids.\n')
elif drv is not None:
display('\nSo far, reducing the density is not supported for the derivative of the density.\n')
else:
rho = integrate.simps(rho, grid.x, dx=grid.delta_[0], axis=0, even='avg')
rho = integrate.simps(rho, grid.y, dx=grid.delta_[1], axis=0, even='avg')
if not options.no_output:
if 'h5' in otype:
display('Saving to Hierarchical Data Format file (HDF5)...')
group = '/mo_set:%03d' % (i_file+1)
display('\n\t%s.h5 in the group "%s"' % (ofid,group))
output.HDF5_creator(rho,ofid,qc.geo_info,qc.geo_spec,data_id='rho',
mode='w',group=group,mo_spec=qc_select['mo_spec'])
if drv is not None:
for i,j in enumerate(drv):
data_id = 'rho_d%s' % j
output.HDF5_creator(delta_rho[i],ofid,qc.geo_info,qc.geo_spec,
data_id=data_id,data_only=True,mode='a',
group=group,mo_spec=qc_select['mo_spec'])
if laplacian:
data_id = 'rho_laplacian'
output.HDF5_creator(laplacian_rho,ofid,qc.geo_info,qc.geo_spec,
data_id=data_id,data_only=True,mode='a',
group=group,mo_spec=qc_select['mo_spec'])
fid = '%s_%03d' % (ofid, i_file+1)
cube_files.append('%s.cb' % fid)
comments = ('mo_set:'+','.join(j_file))
output.main_output(rho,qc.geo_info,qc.geo_spec,outputname=fid,
otype=otype,omit=['h5','vmd','mayavi'],
comments=comments)
if drv is not None:
for i,j in enumerate(drv):
fid = '%s_%03d_d%s' % (ofid, i_file+1, j)
cube_files.append('%s.cb' % fid)
comments = ('d%s_of_mo_set:' % j + ','.join(j_file))
output.main_output(delta_rho[i],qc.geo_info,qc.geo_spec,outputname=fid,
otype=otype,omit=['h5','vmd','mayavi'],
comments=comments)
if laplacian:
fid = '%s_%03d_laplacian' % (ofid, i_file+1)
cube_files.append('%s.cb' % fid)
comments = ('laplacian_of_mo_set:' + ','.join(j_file))
output.main_output(laplacian_rho,qc.geo_info,qc.geo_spec,outputname=fid,
otype=otype,omit=['h5','vmd','mayavi'],
comments=comments)
if 'vmd' in otype and cube_files != []:
display('\nCreating VMD network file...' +
'\n\t%(o)s.vmd' % {'o': ofid})
output.vmd_network_creator(ofid,cube_files=cube_files)
datasets = numpy.array(datasets)
if drv is None:
if 'mayavi' in otype:
output.main_output(datasets,qc.geo_info,qc.geo_spec,
otype='mayavi',datalabels=mo_info['mo_in_file'])
return datasets, mo_info
else:
delta_datasets = numpy.array(delta_datasets)
if 'mayavi' in otype:
datalabels = []
for i in mo_info['mo_in_file']:
datalabels.extend(['d/d%s of %s' % (j,i) for j in drv])
if laplacian:
datalabels.append('laplacian of %s' % i)
output.main_output(delta_datasets.reshape((-1,) + grid.get_shape()),
qc.geo_info,qc.geo_spec,otype='mayavi',datalabels=datalabels)
return datasets, delta_datasets, mo_info
# mo_set
[docs]def calc_ao(qc, drv=None, otype=None, ofid=None,
numproc=None, slice_length=None):
'''Computes and saves all atomic orbital or a derivative thereof.
**Parameters:**
qc.geo_spec, qc.geo_info, qc.ao_spec, qc.mo_spec :
See :ref:`Central Variables` for details.
drv : string or list of strings {None,'x','y', 'z', 'xx', 'xy', ...}, optional
If not None, a derivative calculation of the atomic orbitals
is requested.
otype : str or list of str, optional
Specifies output file type. See :data:`otypes` for details.
ofid : str, optional
Specifies output file name. If None, the filename will be based on
:mod:`orbkit.options.outputname`.
numproc : int, optional
Specifies number of subprocesses for multiprocessing.
If None, uses the value from :mod:`options.numproc`.
slice_length : int, optional
Specifies the number of points per subprocess.
If None, uses the value from :mod:`options.slice_length`.
**Returns:**
ao_list : numpy.ndarray, shape=((NAO,) + N)
Contains the computed NAO atomic orbitals on a grid.
'''
slice_length = options.slice_length if slice_length is None else slice_length
numproc = options.numproc if numproc is None else numproc
if qc.ao_spherical is None:
lxlylz,assign = core.get_lxlylz(qc.ao_spec,get_assign=True)
datalabels = ['lxlylz=%s,atom=%d' % (lxlylz[i],qc.ao_spec[assign[i]]['atom']) for i in range(len(lxlylz))]
else:
datalabels = ['l,m=%s,atom=%d' % (j,qc.ao_spec[i]['atom']) for i,j in qc.ao_spherical]
ao_list = core.rho_compute(qc,
calc_ao=True,
drv=drv,
slice_length=options.slice_length,
numproc=options.numproc)
if otype is None:
return ao_list
if ofid is None:
ofid = '%s_AO' % (options.outputname)
if not options.no_output:
if 'h5' in otype:
output.main_output(mo_list,qc.geo_info,qc.geo_spec,data_id='AO',
outputname=ofid,drv=drv,is_mo_output=False)
# Create Output
cube_files = []
for i in range(len(datalabels)):
outputname = '%s_%03d' % (ofid,i)
index = numpy.index_exp[:,i] if drv is not None else i
output_written = output.main_output(ao_list[index],
qc.geo_info,qc.geo_spec,
outputname=outputname,
comments=datalabels[i].replace('[','').replace(']',''),
otype=otype,omit=['h5','vmd','mayavi'],
drv=drv)
for i in output_written:
if i.endswith('.cb'):
cube_files.append(i)
if 'vmd' in otype and cube_files != []:
display('\nCreating VMD network file...' +
'\n\t%(o)s.vmd' % {'o': ofid})
output.vmd_network_creator(ofid,cube_files=cube_files)
if 'mayavi' in otype:
if drv is not None:
tmp = []
for i in drv:
for j in datalabels:
tmp.append('d/d%s of %s' % (i,j))
datalabels = tmp
data = ao_list.reshape((-1,) + grid.get_shape())
output.main_output(data,qc.geo_info,qc.geo_spec,
otype='mayavi',datalabels=datalabels)
return ao_list
[docs]def atom2index(atom,geo_info=None):
'''Converts a list of atom numbers to indices of :data:`geo_info`.
**Parameters:**
atom : int or list of int
Specifies the indices of the atoms (counting from one).
geo_info : optional
See :ref:`Central Variables` for details.
**Returns:**
atom : list of int
Contains the indices of the atoms (counting from one).
index : numpy.ndarray, dtype=int
Contains the indices of the atoms as occuring in qc.geo_info
(counting from zero).
'''
if not (isinstance(atom,list) or isinstance(atom,numpy.ndarray)):
atom = [atom]
index = []
if geo_info is not None:
geo_info = numpy.array(geo_info)[:,1]
for a in atom:
i = numpy.argwhere(geo_info.astype(int) == a)
if len(i) != 1:
raise ValueError('No or multiple occurence ' +
'of the atom number %d in geo_info!' % a)
index.append(int(i))
index = numpy.array(index, dtype=int)
else:
try:
index = numpy.array(atom,dtype=int)
except ValueError:
raise ValueError('Cannot convert atom to integer array!')
return atom, index
[docs]def gross_atomic_density(atom,qc,
bReturnmo=False,ao_list=None,mo_list=None,drv=None):
r'''Computes the gross atomic density with respect to the selected atoms.
.. math::
\rho^a = \sum_i^{N_{\rm MO}} {\rm occ}_i \cdot \varphi_i^a \cdot \varphi_i
**Parameters:**
atom : 'all' or int or list of int
Specifies the atoms (counting from one) for which the gross atomic density
will be computed.
If (atom == 'all') or (atom == -1), computes the gross atomic density
for all atoms.
qc.geo_spec, qc.geo_info, qc.ao_spec, qc.mo_spec :
See :ref:`Central Variables` for details.
bReturnmo : bool, optional
If True, the gross atomic molecular orbitals are additionally returned.
**Returns:**
rho_atom : list of numpy.ndarrays, shape=(len(atoms,) + N)
Contains the atom gross atomic density on a grid.
mo_atom : list of numpy.ndarrays, shape=(len(atoms,NMO) + N)
Contains the NMO=len(mo_spec) gross atomic molecular orbitals on a grid.
'''
if (atom == 'all' or atom == -1):
atom = range(1, len(qc.geo_info) + 1)
atom, index = atom2index(atom,geo_info=qc.geo_info)
display('Computing the gross atomic density with respect to '+
'the atom(s) (internal numbering)')
outp = '\t['
for i,a in enumerate(atom):
if not i % 10 and i != 0:
outp += '\n\t'
outp += '%d,\t' % a
display('%s]\n' % outp[:-2])
display('\tCalculating ao_list & mo_list')
if ao_list is None:
ao_list = core.ao_creator(qc.geo_spec,qc.ao_spec,ao_spherical=qc.ao_spherical,drv=drv)
if mo_list is None:
mo_list = core.mo_creator(ao_list,qc.mo_spec)
display('\tCalculating the gross atomic density')
N = mo_list.shape[1:]
if bReturnmo: mo_atom = [[] for a in index]
rho_atom = [numpy.zeros(N) for a in index]
for i,a in enumerate(index):
display('\t\tFinished %d of %d' % (i+1,len(index)))
ao_index = []
ao = []
ll = 0
for ii in qc.ao_spec:
for l in range(core.l_deg(l=ii['type'])):
if ii['atom'] == a:
ao_index.append(ll)
ao.append(ao_list[ll])
ll += 1
for ii_mo,spec in enumerate(qc.mo_spec):
mo_info = numpy.zeros(N)
for jj in range(len(ao_index)):
mo_info += spec['coeffs'][ao_index[jj]] * ao[jj]
rho_atom[i] += spec['occ_num'] * mo_list[ii_mo] * mo_info
if bReturnmo: mo_atom[i].append(mo_info)
string = 'Returning the gross atomic density'
if bReturnmo:
display(string + ' and\n\tthe gross atomic molecular orbitals')
return rho_atom, mo_atom
else:
display(string)
return rho_atom
# gross_atomic_density
[docs]def numerical_mulliken_charges(atom,qc,
ao_list=None,mo_list=None,rho_atom=None):
r'''Compute the Mulliken charges and gross atomic populations of the selected
atoms *numerically* using the respective gross atomic densities.
.. math::
P^a = \int \rho^a(x,y,z) {\rm d}^3r
**Parameters:**
atom : 'all' or int or list of int
Specifies the atom (numbering starting from one) to which the numerical
Mulliken charges and gross populations will be computed.
If (atom == 'all') or (atom == -1), computes the charges and populations
for all atoms.
**Returns:**
rho_atom : list of numpy.ndarrays, shape=(len(atoms,) + N)
Contains the atom gross atomic density on a grid.
mulliken_num : dict, (present if not is_vector)
Contains information of Mulliken charge analysis and has following members:
:population: - Mulliken population for each atom.
:charges: - Mulliken charges for each atom.
'''
if (atom == 'all' or atom == -1 or atom == [-1]):
atom = range(1, len(qc.geo_info) + 1)
atom, index = atom2index(atom,geo_info=qc.geo_info)
rho_atom = gross_atomic_density(atom,qc,
ao_list=ao_list,mo_list=mo_list)
if grid.is_vector:
display('Warning: You have applied a vector grid.\n' +
'The integration is not implemented for such a grid!\n' +
'\nReturning only the gross atomic densities...')
return rho_atom
GP_A = numpy.array([core.integration(i) for i in rho_atom])
mulliken_num = {'population': GP_A,
'charge': numpy.array(qc.geo_info[:,2],dtype=float)-GP_A}
display('\nNumerical Mulliken Charges Q and Gross Atomic Populations GP:')
for i,a in enumerate(index):
a = int(a)
display('\tAtom %s (%s):\tQ = %+0.4f ( GP = %0.4f )' %
(qc.geo_info[a][1],qc.geo_info[a][0],
mulliken_num['charge'][i],mulliken_num['population'][i]))
return rho_atom, mulliken_num
[docs]def mo_transition_flux_density(i,j,qc,drv='x',
ao_list=None,mo_list=None,
delta_ao_list=None,delta_mo_list=None):
'''Calculate one component (e.g. drv='x') of the
transition electoronic flux density between the
molecular orbitals i and j.
.. math::
moTEFD_{i,j}(r) = <mo_i|\delta(r-r')\\nabla_x|mo_j>
**Parameters:**
i: int
index of the first mo (Bra)
j: int
index of the second mo (Ket)
drv: {0,1,2,'x','y','z'}
The desired component of the vector field of the
transition electoronic flux density
**Returns:**
mo_tefd : numpy.ndarray
'''
if mo_list is None:
if ao_list is None:
display('\tComputing ao_list and ' +
'mo #%d, since it is not given.' % i)
ao_list = core.ao_creator(qc.geo_spec,qc.ao_spec)
else:
display('\tComputing mo #%d, since it is not given.' % i)
mo = core.mo_creator(ao_list,[qc.mo_spec[i]])[0]
else:
mo = mo_list[i]
if delta_mo_list is None:
if delta_ao_list is None:
display('\tComputing delta_ao_list and the derivative of ' +
'mo #%d, since it is not given.' % j)
delta_ao_list = core.ao_creator(qc.geo_spec,qc.ao_spec,drv=drv)
else:
display('\tComputing mo #%d, since it is not given.' % j)
delta_mo = core.mo_creator(delta_ao_list,[qc.mo_spec[j]])[0]
else:
delta_mo = delta_mo_list[i]
return mo*delta_mo
# mo_transition_flux_density