# -*- coding: utf-8 -*-
#The MIT License (MIT)
#
#Copyright (c) 2014 Lukasz Mentel
#
#Permission is hereby granted, free of charge, to any person obtaining a copy
#of this software and associated documentation files (the "Software"), to deal
#in the Software without restriction, including without limitation the rights
#to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
#copies of the Software, and to permit persons to whom the Software is
#furnished to do so, subject to the following conditions:
#
#The above copyright notice and this permission notice shall be included in all
#copies or substantial portions of the Software.
#
#THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
#IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
#FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
#AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
#LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
#OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
#SOFTWARE.
from __future__ import division, print_function
import json
import os
import pickle
import re
from collections import OrderedDict
from copy import copy, deepcopy
from itertools import chain
import numpy as np
from scipy.linalg import sqrtm, inv
from scipy.special import factorial, factorial2, binom
from chemtools.basisparse import (parse_basis, merge_exponents, CFDTYPE, get_l,
NumpyEncoder)
[docs]class BasisSet(object):
'''
Basis set module supporting basic operation on basis sets.
Args:
name : str
Name of the basis set
element : str
Symbol of the element
Kwargs:
kind : str
Classification of the basis set functions, *diffuse*, *tight*
family : str
basis set family
functions : dict
Dict of functions with *s*, *p*, *d*, *f*, ... as keys
params : list of dicts
Parameters for generating the functions according to the model
'''
def __init__(self, name, element, family=None, kind=None,
functions=None, info=None):
self.name = name
self.element = element
self.family = family
self.kind = kind
self.functions = functions
self.info = info
def __add__(self, other):
'''Add functions from another BasisSet object
Args:
other : BasisSet
BasisSet object whose functions will be added to the existing ones
Returns:
BasisSet instance with functions from `self` and `other` merged
'''
newf = deepcopy(self.functions)
for oshell, ofs in other.functions.items():
if oshell.lower() in self.functions.keys():
i = self.functions[oshell]['e'].size
newf[oshell]['e'] = np.concatenate((self.functions[oshell]['e'],
ofs['e']))
newcf = [np.copy(cf) for cf in ofs['cf']]
for cf in newcf:
cf['idx'] = cf['idx'] + i
newf[oshell]['cf'].extend(newcf)
else:
newf[oshell] = deepcopy(ofs)
return BasisSet(name=self.name, element=self.element, functions=newf)
[docs] @staticmethod
def from_pickle(fname, **kwargs):
'''Read a pickled BasisSet object from a pickle file
Args:
fname : str
File name containing the BasisSet
kwargs : dict
Extra arguments for the `pickle.load` method
Raises:
UnicodeDecodeError
When you try to read a python2 pickle using python3, to fix
that use `encoding='latin1'` option
'''
with open(fname, 'rb') as fil:
return pickle.load(fil, **kwargs)
[docs] @classmethod
def from_optpars(cls, x0, funs=None, name=None, element=None,
explogs=False):
'''
Return a basis set object generated from a sequence based on the
specified arguments.
Args:
x0 : list or numpy.array
Parameters to generate the basis set as a continuous list or array
funs : list of tuples
A list of tuple specifying the shell type, number of functions and
parameters, e.g. `[('s', 'et', 4, (0.5, 2.0)), ('p', 'et', 3,
(1.0, 3.0))]`
name : str
Name of the basis set
element : str
Chemical symbol of the element
Returns:
out : BasisSet
'''
functions = dict()
ni = 0
nt = 0
for shell, seq, nf, params in funs:
functions[shell] = dict()
if seq in ["exp", "exponents"]:
nt += nf
if explogs:
functions[shell]['e'] = generate_exponents(seq, nf, np.exp(x0[ni: nt]))
else:
functions[shell]['e'] = generate_exponents(seq, nf, x0[ni: nt])
ni += nf
else:
# TODO:
# params shouldn't be here since it can have a dummy,
# real values are taken from x0
nt += len(params)
functions[shell]['e'] = generate_exponents(seq, nf, x0[ni:nt])
ni += len(params)
functions = OrderedDict(sorted(functions.items(),
key=lambda x: get_l(x[0])))
bs = cls(name=name, element=element, functions=functions)
bs.uncontract()
return bs
[docs] @classmethod
def from_sequence(cls, funs=None, name=None, element=None):
'''
Return a basis set object generated from a sequence based on the
specified arguments.
Args:
funs : list of tuples
A list of tuple specifying the shell type, number of functions and
parameters, e.g. `[('s', 'et', 4, (0.5, 2.0)), ('p', 'et', 3,
(1.0, 3.0))]`
name : str
Name of the basis set
element : str
Chemical symbol of the element
Returns:
out : BasisSet
'''
functions = dict()
for shell, seq, nf, params in funs:
functions[shell] = dict()
functions[shell]['e'] = generate_exponents(seq, nf, params)
functions = OrderedDict(sorted(functions.items(),
key=lambda x: get_l(x[0])))
bs = cls(name=name, element=element, functions=functions)
bs.uncontract()
bs.sort()
return bs
[docs] @classmethod
def from_file(cls, fname=None, fmt=None, name=None):
'''Read and parse a basis set from file and return a BasisSet object
Args:
fname : str
File name
fmt : str
Format of the basis set in the file (*molpro*, *gamessus*)
name : str
Name of the basis set
Returns:
out : BasisSet or dict
Basisset object parsed from file or dictionary of BasisSet objects
'''
if name is None:
name = os.path.splitext(fname)[0]
with open(fname, 'r') as fobj:
basstr = fobj.read()
return BasisSet.from_str(basstr, fmt=fmt, name=name)
[docs] @classmethod
def from_str(cls, string, fmt=None, name=None):
'''
Parse a basis set from string
Args:
string : str
A string with the basis set
fmt : str
Format of the basis set in the file: *molpro*, *gamessus*
name : str
Name of the basis set
Returns:
out : BasisSet or dict
Basisset object parsed from string or dictionary of BasisSet
objects
'''
res = parse_basis(string, fmt=fmt)
# return a BasisSet object if only one basis parsed or dict of BasisSet
# objects with atomic symbol as key
out = dict()
if len(res) == 1:
atom, fs = list(res.items())[0]
return cls(name=name, element=list(res.keys())[0],
functions=OrderedDict(sorted(fs.items(),
key=lambda x: get_l(x[0]))))
else:
for atom, fs in res.items():
out[atom] = cls(name=name, element=atom,
functions=OrderedDict(sorted(fs.items(),
key=lambda x: get_l(x[0]))))
return out
[docs] @classmethod
def from_json(cls, jsonstring, **kwargs):
'''
Instantiate the `BasisSet` object from a JSON string
Args:
jsonstring: str
A JSON serialized string with the basis set
'''
def json_numpy_obj_hook(dct):
'''
Decodes a previously encoded numpy ndarray with proper shape
and dtype.
Args:
dct: (dict) json encoded ndarray
Returns:
(ndarray) if input was an encoded ndarray
'''
if isinstance(dct, dict) and 'dtype' in dct:
if dct['dtype'] == 'CFDTYPE':
return np.array([tuple(r) for r in dct['data']],
dtype=CFDTYPE)
else:
return np.array(dct['data'], dtype=dct['dtype'])
return dct
loaded = json.loads(jsonstring, object_hook=json_numpy_obj_hook,
**kwargs)
return cls(**loaded)
[docs] def append(self, other):
'''Append functions from another BasisSet object
Args:
other : BasisSet
BasisSet object whose functions will be added to the existing ones
'''
for oshell, ofs in other.functions.items():
if oshell.lower() in self.functions.keys():
i = self.functions[oshell]['e'].size
self.functions[oshell]['e'] = np.concatenate((self.functions[oshell]['e'], ofs['e']))
newcf = [np.copy(cf) for cf in ofs['cf']]
for cf in newcf:
cf['idx'] = cf['idx'] + i
self.functions[oshell]['cf'].extend(newcf)
else:
self.functions[oshell] = deepcopy(ofs)
[docs] def completeness_profile(self, zetas):
'''
Calculate the completeness profile of each shell of the basis set
Args:
zetas : numpy.array
Scaning exponents
Returns:
out : numpy.array
numpy array with values of the profile (shells, zetas)
'''
out = np.zeros((zetas.size, len(self.functions.keys())))
for i, (shell, fs) in enumerate(self.functions.items()):
cc = self.contraction_matrix(shell)
# calculate the overlap matrix
S = self.shell_overlap(shell)
# calculate the inverse square root of S
X = inv(sqrtm(S))
SO = primitive_overlap(get_l(shell), fs['e'], zetas)
J = np.dot(np.dot(cc, X).T, SO)
Y = np.sum(np.square(J), axis=0)
out[:, i] = Y
return out
[docs] def contraction_matrix(self, shell):
'''
Return the contraction coefficients for a given shell in a
matrix form with size `ne * nc`, where `ne` is the number of
exponents and `nc` is the number of contracted functions
Args:
shell : str
shell label, *s*, *p*, *d*, ...
Returns:
out : 2D numpy.array
2D array with contraction coefficients
'''
if shell in self.functions.keys():
fs = self.functions[shell]
out = np.zeros((fs['e'].size, len(fs['cf'])))
for i, cf in enumerate(fs['cf']):
out[cf['idx'], i] = cf['cc']
return out
else:
raise ValueError("shell '{}' is not present in the BasisSet".format(shell))
[docs] def contractions_per_shell(self):
'''
Calculate how many contracted functions are in each shell.
Returns:
out : list of ints
'''
return [len(f['cf']) for s, f in self.functions.items()]
[docs] def contraction_scheme(self):
'''
Return a string describing the contraction scheme.
'''
cs, ec = [], []
for shell, fs in self.functions.items():
cs.append((shell, len(fs['e']), len(fs['cf'])))
ec.append([len(cfs) for cfs in fs['cf']])
return "({p:s}) -> [{c:s}] : {{{ec}}}".format(
p="".join(["{0:d}{1:s}".format(c[1], c[0]) for c in cs]),
c="".join(["{0:d}{1:s}".format(c[2], c[0]) for c in cs]),
ec="/".join([" ".join(["{0:d}".format(c) for c in x]) for x in ec]))
[docs] def contraction_type(self):
'''
Try to determine the contraction type: segmented, general,
uncontracted, unknown.
'''
pps = self.primitives_per_shell()
ppc = self.primitives_per_contraction()
if any(x > 1 for x in pps):
if all(all(x == 1 for x in shell) for shell in ppc):
return "uncontracted"
elif all(all(pinc == np for pinc in shell)
for np, shell in zip(pps, ppc)):
return "general"
else:
return "unknown"
# one function per shell case
if all(all(x == 1 for x in shell) for shell in ppc):
return "uncontracted 1fps"
[docs] def get_exponents(self, asdict=True):
'''
Return the exponents of a given shell or if the shell isn't
specified return all of the available exponents
Args:
asdict (bool): if `True` a dict with exps per shell is
returned
'''
if asdict:
return {shell: f['e'] for shell, f in self.functions.items()}
else:
exps = [self.functions[k]['e'] for k in self.functions.keys()]
array = np.fromiter(chain.from_iterable(exps), dtype=np.float)
return np.sort(array)[::-1]
def get_filename(self, ext=None):
name = self.name.strip().replace(' ', '_') + '-' + self.element
if ext:
return name + '.' + ext
return name
[docs] def nf(self, spherical=True):
'''
Calculate the number of basis functions
Args:
spherical : bool
flag indicating if spherical or cartesian functions should be used,
default: True
Returns:
res : int
number of basis functions
'''
if spherical:
return sum(nspherical(get_l(shell)) * len(fs['cf'])
for shell, fs in self.functions.items())
else:
return sum(ncartesian(get_l(shell)) * len(fs['cf'])
for shell, fs in self.functions.items())
[docs] def normalize(self):
'''
Normalize contraction coefficients for each contracted functions based
on the primitive overlaps so that the norm is equal to 1.
'''
for shell, fs in self.functions.items():
cc = self.contraction_matrix(shell)
po = primitive_overlap(get_l(shell), fs['e'], fs['e'])
for col in range(cc.shape[1]):
norm2 = np.dot(cc[:, col], np.dot(po, cc[:, col]))
fs['cf'][col]['cc'] = cc[fs['cf'][col]['idx'], col] / np.sqrt(norm2)
[docs] def normalization(self):
'''
For each function (contracted) calculate the norm and return a list of
tuples containing the shell, function index and the norm respectively.
'''
out = list()
for shell, fs in self.functions.items():
cc = self.contraction_matrix(shell)
po = primitive_overlap(get_l(shell), fs['e'], fs['e'])
for col in range(cc.shape[1]):
norm2 = np.dot(cc[:, col], np.dot(po, cc[:, col]))
out.append(tuple([shell, col, norm2]))
return out
[docs] def nprimitive(self, spherical=True):
'''
Return the number of primitive functions assuming sphrical or
cartesian gaussians.
Args:
spherical : bool
A flag to select either spherical or cartesian gaussians
Returns:
out : int
Number of primitive function in the basis set
'''
if spherical:
# calculate the number of spherical components per shell
ncomp = [nspherical(get_l(shell))
for shell in self.functions.keys()]
else:
# calculate the number of cartesian components per shell
ncomp = [ncartesian(get_l(shell))
for shell in self.functions.keys()]
return sum([prim * nc for prim, nc in zip(self.primitives_per_shell(),
ncomp)])
[docs] def partial_wave_expand(self):
'''
From a given basis set with shells spdf... return a list of basis sets
that are subsets of the entered basis set with increasing angular
momentum functions included [s, sp, spd, spdf, ...]
'''
res = list()
shells = self.functions.keys()
for i in range(1, len(shells) + 1):
bscopy = copy(self)
bscopy.functions = {k: v for k, v in self.functions.items()
if k in shells[:i]}
res.append(bscopy)
return res
[docs] def primitives_per_shell(self):
'''
Calculate how many primitive functions are in each shell.
Returns:
out : list of ints
'''
return [len(f['e']) for s, f in self.functions.items()]
[docs] def primitives_per_contraction(self):
'''
Calculate how many primities are used in each contracted function.
Returns:
out : list of ints
'''
return [[len(cc) for cc in f['cf']] for s, f in self.functions.items()]
[docs] def uncontract(self, copy=False):
'''
Uncontract the basis set. This replaces the contraction coefficients in
the current object.
Args:
copy : bool
If `True` return an uncontracted copy of the basis set rather than
uncontracting in place, default is `False`.
'''
if copy:
res = deepcopy(self)
for shell, fs in res.functions.items():
fs['cf'] = [np.array([tuple([i, 1.0])], dtype=CFDTYPE)
for i, _ in enumerate(fs['e'])]
return res
else:
for shell, fs in self.functions.items():
fs['cf'] = [np.array([tuple([i, 1.0])], dtype=CFDTYPE)
for i, _ in enumerate(fs['e'])]
[docs] def shell_overlap(self, shell):
'''
Calculate the overlap integrals for a given shell
Args:
shell : str
Shell
Returns:
out : numpy.array
Overlap integral matrix
'''
exps = self.functions[shell]['e']
S = primitive_overlap(get_l(shell), exps, exps)
cc = self.contraction_matrix(shell)
return np.dot(cc.T, np.dot(S, cc))
[docs] def sort(self, reverse=False):
'''
Sort shells in the order of increasing angular momentum and for each
shell sort the exponents.
Args:
reverse : bool
If `False` sort the exponents in each shell in the descending order
(default), else sort exponents in ascending order
'''
self.functions = OrderedDict(sorted(self.functions.items(),
key=lambda x: get_l(x[0])))
for shell, fs in self.functions.items():
if reverse:
idx = np.argsort(fs['e'])
else:
idx = np.argsort(fs['e'])[::-1]
# actually sort the exponents and coefficients
fs['e'] = fs['e'][idx]
for cf in fs['cf']:
cf['idx'] = np.asarray([np.nonzero(idx == y)[0][0]
for y in cf['idx']])
cf.sort(order='idx')
[docs] def to_cfour(self, comment="", efmt="15.8f", cfmt="15.8f"):
'''
Return a string with the basis set in (new) CFOUR format.
Args:
comment : str
comment string
efmt : str
string describing output format for the exponents,
default: "20.10f"
cfmt : str
string describing output format for the contraction
coefficients, default: "15.8f"
Returns:
res : str
basis set string in Cfour format
'''
am, ne, cf = [], [], []
for shell, shellfs in sorted(self.functions.items(),
key=lambda x: get_l(x[0])):
am.append(get_l(shell))
ne.append(len(shellfs["e"]))
cf.append(len(shellfs["cf"]))
res = "\n{e}:{s}\n{c}\n\n".format(e=self.element, s=self.name,
c=comment)
res += "{0:3d}\n".format(len(self.functions.keys()))
res += "".join(["{0:5d}".format(x) for x in am]) + "\n"
res += "".join(["{0:5d}".format(x) for x in cf]) + "\n"
res += "".join(["{0:5d}".format(x) for x in ne]) + "\n"
res += "\n"
for shell, fs in self.functions.items():
for lst in splitlist(fs['e'], 5):
res += "".join(["{0:>{efmt}}".format(e, efmt=efmt)
for e in lst]) + "\n"
res += "\n"
# create an array with all the contraction coefficients for a given shell
cc = self.contraction_matrix(shell)
for row in cc:
for splitrow in splitlist(row, 5):
res += "{c}".format(c="".join(["{0:{cfmt}}".format(c, cfmt=cfmt) for c in splitrow])) + "\n"
res += "\n"
return res
[docs] def to_dalton(self, fmt='prec'):
'''
Return a string with the basis set in DALTON format.
Args:
fmt : str
string describing output format for the exponents and coefficents
- `prec` "20.10f"
- `default` "10.4f"
- or python format string e.g. "15.8f"
Returns:
res : str
basis set string in Dalton format
'''
formats = {'prec': '20.10f', 'default': '10.4f'}
ffmt = formats.get(fmt, fmt)
if fmt == 'prec':
ffmt = '20.10f'
nitems = 3
fmtlabel = 'H'
elif fmt == 'default':
ffmt = '10.4f'
nitems = 7
fmtlabel = ' '
else:
ffmt = fmt
cwidth = int(ffmt.split('.')[0])
nitems = (80 - cwidth) // cwidth
fmtlabel = "{0:d}F{1}.{2}".format(nitems, ffmt.split('.')[0], re.sub('[A-Za-z]+', '', ffmt.split('.')[1]))
res = "! {s}\n".format(s=self.name)
for shell, fs in self.functions.items():
res += "! {s} functions\n".format(s=shell)
res += "{f:1s}{p:>4d}{c:>4d}\n".format(f=fmtlabel, p=len(fs['e']),
c=len(fs['cf']))
# create an array with all the contraction coefficients for a given shell
cc = self.contraction_matrix(shell)
for expt, row in zip(fs['e'], cc):
it = splitlist(row, nitems)
res += "{e:>{efmt}}{c}".format(e=expt, efmt=ffmt, c="".join(["{0:{cfmt}}".format(c, cfmt=ffmt) for c in next(it)])) + "\n"
for row in it:
res += ' ' * int(ffmt.split('.')[0]) + \
"".join(["{0:{cfmt}}".format(c, cfmt=ffmt) for c in row]) + "\n"
return res
[docs] def to_gamessus(self, efmt="20.10f", cfmt="15.8f"):
'''
Return a string with the basis set in GAMESS(US) format.
Args:
efmt : str
string describing output format for the exponents,
default: "20.10f"
cfmt : str
string describing output format for the contraction
coefficients, default: "15.8f"
Returns:
res : str
basis set string in Gamess(US) format
'''
res = ''
for shell, fs in self.functions.items():
for contraction in fs["cf"]:
res += "{s:<1s}{n:>3d}\n".format(s=shell.upper(),
n=len(contraction))
for i, (idx, coeff) in enumerate(contraction, start=1):
res += "{i:3d}{e:>{efmt}}{c:>{cfmt}}".format(i=i,
e=fs['e'][idx],
efmt=efmt,
c=coeff,
cfmt=cfmt) + "\n"
return res + "\n"
[docs] def to_gaussian(self, efmt="20.10f", cfmt="15.8f"):
'''
Return a string with the basis set in Gaussian format.
Args:
efmt : str
string describing output format for the exponents,
default: "20.10f"
cfmt : str
string describing output format for the contraction
coefficients, default: "15.8f"
Returns:
res : str
basis set string in Gaussian format
'''
res = "****\n{e:<7s}0\n".format(e=self.element)
for shell, fs in self.functions.items():
for contraction in fs["cf"]:
res += "{s:<1s}{n:>4d}{i:>7.2f}\n".format(s=shell.upper(),
n=len(contraction),
i=1.0)
for idx, coeff in contraction:
res += "{e:>{efmt}}{c:>{cfmt}}".format(e=fs['e'][idx],
efmt=efmt, c=coeff,
cfmt=cfmt) + "\n"
return res + "****\n"
[docs] def to_json(self, fname=None, **kwargs):
'''
Serizalize the basis set object to JSON format
'''
if fname:
with open(fname, 'w') as fjson:
json.dump(self.__dict__, fjson, cls=NumpyEncoder, **kwargs)
return json.dumps(self.__dict__, cls=NumpyEncoder, **kwargs)
[docs] def to_latex(self, efmt="20.10f", cfmt="15.8f"):
'''
Return a string with the basis set as LaTeX table/
Args:
efmt : str
Output format for the exponents, default: "20.10f"
cfmt : str
Output format for the contraction coefficients,
default: "15.8f"
Returns:
res : str
basis set string in LaTeX format
'''
# get the number of contracted functions per shell
ncf = [sum([len(cfs) > 1 for cfs in fs['cf']])
for sh, fs in self.functions.items()]
idx = np.argmax(ncf)
nccols = ncf[idx]
out = '\\begin{{tabular}}{{{}}}\n'.format('r' * (nccols + 2))
out += 'No. & \multicolumn{1}{c}{Exponent} & ' +\
'\multicolumn{{{0:d}}}{{c}}{{Coefficients }} \\\ \n'.format(nccols)
for shell, fs in self.functions.items():
out += '\hline \n'
out += '\multicolumn{{{0:d}}}{{c}}{{ {1:s} shell }} \\\ \hline \n'.format(nccols + 2, shell)
cc = self.contraction_matrix(shell)
count = 0
# select columns with more than 1 non zero coefficients
nonzerocolmask = np.array([np.count_nonzero(col) > 1
for col in cc.T])
nonzerorowmask = np.array([np.count_nonzero(row) > 0
for row in cc[:, nonzerocolmask]])
if np.any(nonzerocolmask):
for expt, cfc in zip(fs['e'][nonzerorowmask],
cc[np.ix_(nonzerorowmask,
nonzerocolmask)]):
count += 1
coeffrow = " & ".join(["{0:{cfmt}}".format(c, cfmt=cfmt)
for c in cfc])
out += "{i:5d} & {e:>{efmt}} & {c}".format(i=count, e=expt,
efmt=efmt,
c=coeffrow) + " \\\ \n"
if np.any(np.logical_not(nonzerocolmask)):
for colidx in np.where(np.logical_not(nonzerocolmask))[0]:
count += 1
nonzerorowmask = np.array([np.count_nonzero(row) > 0
for row in cc[:, colidx]])
e = fs['e'][nonzerorowmask][0]
c = cc[nonzerorowmask, :][0][colidx]
out += "{i:5d} & {e:>{efmt}} & \\\ \n".format(i=count, e=e,
efmt=efmt)
out += '\end{tabular}'
return out
[docs] def to_molpro(self, withpars=False, efmt="20.10f", cfmt="15.8f"):
'''
Return a string with the basis set in MOLPRO format.
Args:
withpars : bool
A flag to indicate whether to wrap the basis with
`basis={ }` string
efmt : str
Output format for the exponents, default: "20.10f"
cfmt : str
Output format for the contraction coefficients,
default: "15.8f"
Returns:
res : str
basis set string
'''
# reorder indices of cf if necessary without modiyfing the instance
funs = deepcopy(self.functions)
for shell, fs in funs.items():
if not has_consecutive_indices(fs):
funs[shell] = reorder_shell_to_consecutive(fs)
res = ""
for shell, fs in funs.items():
exps = ", ".join(["{0:>{efmt}}".format(e, efmt=efmt).lstrip()
for e in fs['e']]) + '\n'
res += "{s:>s}, {e:>s}, ".format(s=shell, e=self.element) + exps
for icol, cf in enumerate(fs['cf']):
if cf.size == 1:
coeffs = ", ".join(["{0:>{cfmt}}".format(cc,
cfmt=cfmt).lstrip() for cc in cf['cc']])
res += "c, {0:d}.{0:d}, ".format(cf['idx'][0] + 1) + \
coeffs + '\n'
else:
# indices should be consecutive
coeffs = ", ".join(["{0:>{cfmt}}".format(cc,
cfmt=cfmt).lstrip() for cc in cf['cc']])
res += "c, {0:d}.{1:d}, ".format(cf['idx'].min() + 1,
cf['idx'].max() + 1) +\
coeffs + '\n'
if withpars:
res = 'basis={\n' + res + '}'
return res
[docs] def to_nwchem(self, efmt="20.10f", cfmt="15.8f"):
'''
Return a string with the basis set in NWChem format.
Args:
efmt : str
string describing output format for the exponents,
default: "20.10f"
cfmt : str
string describing output format for the contraction
coefficients, default: "15.8f"
Returns:
res : str
basis set string in NwChem format
'''
res = 'BASIS "ao basis" PRINT\n'
for shell, fs in self.functions.items():
# create an array of contraction coefficients for a given shell
cc = self.contraction_matrix(shell)
# select columns with more than 1 non zero coefficients
nonzerocolmask = np.array([np.count_nonzero(col) > 1
for col in cc.T])
nonzerorowmask = np.array([np.count_nonzero(row) > 0
for row in cc[:, nonzerocolmask]])
if np.any(nonzerocolmask):
res += "{e} {s}\n".format(e=self.element, s=shell)
for expt, cfc in zip(fs['e'][nonzerorowmask],
cc[np.ix_(nonzerorowmask, nonzerocolmask)]):
res += "{e:>{efmt}}{c}".format(e=expt, efmt=efmt, c="".join(["{0:{cfmt}}".format(c, cfmt=cfmt) for c in cfc])) + "\n"
if np.any(np.logical_not(nonzerocolmask)):
for colidx in np.where(np.logical_not(nonzerocolmask))[0]:
res += "{e} {s}\n".format(e=self.element, s=shell)
nonzerorowmask = np.array([np.count_nonzero(row) > 0
for row in cc[:, colidx]])
e = fs['e'][nonzerorowmask][0]
c = cc[nonzerorowmask, :][0][colidx]
res += "{e:>{efmt}}{c:>{cfmt}}".format(e=e, efmt=efmt, c=c,
cfmt=cfmt) + "\n"
return res + "END\n"
[docs] def to_pickle(self, fname=None):
'''Save the basis set in pickle format under the filename `fname`
Args:
fname : str
File name
'''
if fname is None:
fname = self.get_filename('pkl')
with open(fname, 'wb') as fbas:
pickle.dump(self, fbas)
[docs] def print_functions(self, efmt="20.10f", cfmt="15.8f"):
'''
Return a string with the basis set.
Args:
efmt : str
string describing output format for the exponents,
default: "20.10f"
cfmt : str
string describing output format for the contraction
coefficients, default: "15.8f"
Returns:
res : str
basis set string
'''
res = ''
for shell, fs in self.functions.items():
# create an array with all the contraction coefficients
# for a given shell
res += "\n" + "{s} shell".format(s=shell).center(40, '=') + "\n"
cc = self.contraction_matrix(shell)
count = 0
# select columns with more than 1 non zero coefficients
nonzerocolmask = np.array([np.count_nonzero(col) > 1
for col in cc.T])
nonzerorowmask = np.array([np.count_nonzero(row) > 0
for row in cc[:, nonzerocolmask]])
if np.any(nonzerocolmask):
res += 'Contracted:\n'
for expt, cfc in zip(fs['e'][nonzerorowmask],
cc[np.ix_(nonzerorowmask,
nonzerocolmask)]):
count += 1
coeffrow = "".join(["{0:{cfmt}}".format(c, cfmt=cfmt)
for c in cfc])
res += "{i:5d}{e:>{efmt}}{c}".format(i=count, e=expt,
efmt=efmt,
c=coeffrow) + "\n"
if np.any(np.logical_not(nonzerocolmask)):
res += 'Uncontracted:\n'
for colidx in np.where(np.logical_not(nonzerocolmask))[0]:
count += 1
nonzerorowmask = np.array([np.count_nonzero(row) > 0
for row in cc[:, colidx]])
e = fs['e'][nonzerorowmask][0]
c = cc[nonzerorowmask, :][0][colidx]
res += "{i:5d}{e:>{efmt}}{c:>{cfmt}}".format(i=count, e=e,
efmt=efmt,
c=c,
cfmt=cfmt) + "\n"
return res
def __repr__(self):
keys = ['name', 'element', 'family', 'kind']
res = "<BasisSet(\n"
for key in keys:
res += "\t{k:<20s} = {v}\n".format(k=key, v=getattr(self, key))
res += self.print_functions()
res += ")>"
return res
def __str__(self):
keys = ['name', 'element', 'family', 'kind']
res = ''
for key in keys:
res += "{k:<20s} = {v}\n".format(k=key.capitalize(),
v=getattr(self, key))
res += 'Functions:\n'
res += self.print_functions()
return res
[docs]def has_consecutive_indices(shell):
'''
Check if all the contracted functions have consecutive indices
Args:
shell : dict
Basis functions for a given shell asa dict with structure
``{'e' : np.array(), 'cf': [np.array(), np.array(), ...]}``
'''
for i, cf in enumerate(shell['cf']):
if len(cf) > 1:
if not np.all(np.sort(cf['idx'])[1:] -
np.sort(cf['idx'])[:-1] == 1):
return False
else:
return True
[docs]def reorder_shell_to_consecutive(shell):
'''
Reorder the exponents so that the indices of the contracted functions
have consecutive inidices.
Args:
shell : dict
Basis functions for a given shell asa dict with structure
``{'e' : np.array(), 'cf': [np.array(), np.array(), ...]}``
Returns:
shell : dict
Same shell as on input but with reordered exponents and relabelled
contracted functions
'''
# get the index of the largest contracted function
cfidx = np.argmax([len(x) for x in shell['cf']])
# get the new indices where the ones corresponding to the contracted
# function are consecutive
oldidx = shell['cf'][cfidx]['idx']
conidx = np.arange(shell['e'].size, dtype=int)
restidx = np.setdiff1d(conidx, oldidx)
newidx = np.concatenate((oldidx, restidx))
funs = {'e': shell['e'][newidx], 'cf': []}
newidxlist = newidx.tolist()
for cf in shell['cf']:
funs['cf'].append(np.array([(newidxlist.index(i), c) for i, c in cf],
dtype=CFDTYPE))
# check if all the contracted functions have consecutive indices
if not has_consecutive_indices(funs):
raise ValueError('cannot reorder functions, check #cf: {}'.format(i))
return funs
[docs]def merge(first, other):
'''
Merge functions from two BasisSet objects
Args:
first : BasisSet
First `BasisSet` object to merge
other : BasisSet
Second `BasisSet` object to merge
Returns:
our : BasisSet
BasisSet instance with functions from `first` and `other` merged
'''
newf = deepcopy(first.functions)
for oshell, ofs in other.functions.items():
if oshell.lower() in first.functions.keys():
exps, idxs, idxo = merge_exponents(first.functions[oshell]['e'], ofs['e'])
newf[oshell]['e'] = exps
for cf in newf[oshell]['cf']:
cf['idx'] = idxs[cf['idx']]
newcf = [np.copy(cf) for cf in ofs['cf']]
for cf in newcf:
cf['idx'] = idxo[cf['idx']]
newf[oshell]['cf'].extend(newcf)
else:
newf[oshell] = deepcopy(ofs)
return BasisSet(name=first.name, element=first.element, functions=newf)
[docs]def primitive_overlap(l, a, b):
'''
Calculate the overlap integrals for a given shell `l` and two sets of
exponents
.. math::
S(\zeta_i, \zeta_j) = \\frac{2^{l + \\frac{3}{2}} \left(\zeta_{1} \zeta_{2}\\right)^{\\frac{l}{2} + \\frac{3}{4}}}{\left(\zeta_{1} + \zeta_{2}\\right)^{l + \\frac{3}{2}}}
Args:
l : int
Angular momentum quantum number of the shell
a (M,) : numpy.array
First vector of exponents
b (N,) : numpy.array
Second vector of exponents
Returns:
out (M, N) : numpy.array
Overlap integrals
'''
return np.power(2.0 * np.sqrt(np.outer(a, b)) / np.add.outer(a, b), l + 1.5)
[docs]def nspherical(l):
'''
Calculate the number of spherical components of a function with a given
angular momentum value `l`.
'''
return 2 * l + 1
[docs]def ncartesian(l):
'''
Calculate the number of cartesian components of a function with a given
angular momentum value *l*.
'''
return int((l + 1) * (l + 2) / 2)
[docs]def zetas2legendre(zetas, kmax):
'''
From a set of exponents (`zetas`), using least square fit calculate the
expansion coefficients into the legendre polynomials of the order `kmax`.
Args:
kmax : int
length of the legendre expansion
zetas : list
list of exponents (floats) to be fitted
Returns:
coeff : numpy.array
numpy array of legendre expansion coeffcients of length kmax
'''
# exponents should be sorted in the ascending order
zetas = sorted(zetas)
c = np.asarray([1.0] * kmax, dtype=float)
leg = np.polynomial.Legendre(c)
a = np.zeros((len(zetas), kmax))
for j in range(len(zetas)):
for k in range(kmax):
arg = (2.0 * (j + 1.0) - 2.0) / (len(zetas) - 1.0) - 1.0
a[j, k] = leg.basis(k)(arg)
return np.linalg.lstsq(a, np.log(zetas))[0]
[docs]def zetas2eventemp(zetas):
'''
From a set of exponents (`zetas`), using least square fit calculate the
approximate $\\alpha$ and $\\beta$ parameters from the even tempered
expansion.
Args:
zetas : list
list of exponents (floats) to be fitted
Returns:
coeff : numpy.array
numpy array of legendre expansion coeffcients of length kmax
'''
return min(zetas), np.power(max(zetas) / min(zetas), 1.0 / (len(zetas) - 1))
def generate_exponents(formula, nf, params):
'''
Generate a sequence of exponents from a specified formula
Args:
formula : str
name of the sequence from which the exponents are generated, one of:
- `et`, `even`, `eventemp`, `even tempered`
- `wt`, `well`, `welltemp`, `well tempered`
- `le`, `legendre`
- `exp`, `exponents`
nf : int
number of exponents
params : tuple
parameters for the formula
Returns:
numpy.array with exponents
'''
if formula.lower() in ["et", "even", "eventemp", "even tempered"]:
return eventemp(nf, params)
elif formula.lower() in ["wt", "well", "welltemp", "well tempered"]:
return welltemp(nf, params)
elif formula.lower() in ["le", "legendre"]:
return legendre(nf, params)
elif formula.lower() in ["exp", "exponents"]:
return np.sort(np.array(params))[::-1]
else:
raise ValueError('unknown sequence name: {}'.format(formula))
def get_num_params(ftuple):
'''
Get the number of parameter for a given 4-tuple desciribing functions
>>> get_num_params(('s', 'et', 20, (0.4, 2.0)))
2
'''
formula = ftuple[1]
if formula.lower() in ["et", "even", "eventemp", "even tempered"]:
return 2
elif formula.lower() in ["wt", "well", "welltemp", "well tempered"]:
return 4
elif formula.lower() in ["le", "legendre"]:
return ftuple[2]
elif formula.lower() in ["exp", "exponents"]:
return ftuple[2]
else:
raise ValueError('unknown sequence name: {}'.format(formula))
[docs]def eventemp(numexp, params):
'''
Generate a sequence of nf even tempered exponents accodring to
the even tempered formula
.. math::
\\zeta_i = \\alpha \cdot \\beta^{i-1}
Args:
numexp : int
Number fo exponents to generate
params : tuple of floats
Alpha and beta parameters
Returns:
res : numpy array
Array of generated exponents (floats)
'''
if not isinstance(numexp, int):
raise TypeError('"nf" variable should be of "int" type, got: {}'.format(type(numexp)))
if len(params) != 2:
raise ValueError('"params" tuple should have exactly 2 entries, got {}'.format(len(params)))
alpha, beta = params
zetas = alpha * np.power(beta, np.arange(numexp))
return zetas[::-1]
[docs]def welltemp(numexp, params):
'''
Generate a sequence of nf well tempered exponents accodring to
the well tempered fromula
.. math::
\\zeta_i = \\alpha \cdot \\beta^{i-1} \cdot \\left[1 + \\gamma \cdot \\left(\\frac{i}{N}\\right)^{\delta}\\right]
Args:
numexp : int
Number fo exponents to generate
params : tuple of floats
Alpha, beta, gamma and delta parameters
Returns:
res : numpy.array
Array of generated exponents (floats)
'''
if not isinstance(numexp, int):
raise TypeError('"nf" variable should be of "int" type, got: {}'.format(type(numexp)))
if len(params) != 4:
raise ValueError('"params" tuple should have exactly 4 entries, got {}'.format(len(params)))
alpha, beta, gamma, delta = params
zetas = alpha * np.power(beta, np.arange(numexp)) * \
(1 + gamma * np.power(np.arange(1, numexp + 1) / numexp, delta))
zetas.sort()
return zetas[::-1]
[docs]def legendre(numexp, coeffs):
'''
Generate a sequence of nf exponents from expansion in the orthonormal
legendre polynomials as described in: Peterson, G. A. et.al J. Chem. Phys.,
Vol. 118, No. 3 (2003), eq. (7).
.. math::
\ln \\zeta_i = \\sum^{k_{\max}}_{k=0} A_k P_k \\left(\\frac{2j-2}{N-1}-1\\right)
Args:
numexp : int
Number fo exponents to generate
params : tuple of floats
Polynomial coefficients (expansion parameters)
Returns:
res : numpy array
Array of generated exponents (floats)
'''
if not isinstance(numexp, int):
raise TypeError('"nf" variable should be of "int" type, got: {}'.format(type(numexp)))
if len(coeffs) < 1:
raise ValueError('"coeffs" tuple should have at least 1 entry, got {}'.format(len(coeffs)))
# special case for one function
if len(coeffs) == 1:
return [np.exp(coeffs[0])]
poly = np.polynomial.legendre.Legendre(coeffs)
zetas = [poly(((2.0 * (i + 1.0) - 2.0) / (numexp - 1.0)) - 1.0) for i in range(numexp)]
return np.exp(zetas[::-1])
def splitlist(l, n):
'''
Split a list into sublists of size `n`
'''
if len(l) % n == 0:
splits = len(l) // n
elif len(l) % n != 0 and len(l) > n:
splits = len(l) // n + 1
else:
splits = 1
for i in range(splits):
yield l[n * i:n * i + n]
def sliceinto(toslice, chunk_sizes):
'''
Slice a list ``toslice`` into chunks of sizes defined in ``chunk_sizes``
>>> sliceinto(list(range(10)), (6, 4))
[[0, 1, 2, 3, 4, 5], [6, 7, 8, 9]]
'''
if len(toslice) != sum(chunk_sizes):
raise ValueError('cannot slice list, size mismatch {} != {}'.format(len(toslice),
sum(chunk_sizes)))
sidx = [sum(chunk_sizes[:i]) for i in range(0, len(chunk_sizes))]
return [toslice[i:i + s] for i, s in zip(sidx, chunk_sizes)]
[docs]def xyzlist(l):
'''
Generate an array of :math:`l_x`, :math:`l_y`, :math:`l_z` components of
cartesian gaussian with a given angular momentum value in canonical order.
For exampe:
- ``l = 0`` generates the result ``array([[0, 0, 0]])``
- ``l = 1`` generates ``array([[1, 0, 0], [0, 1, 0], [0, 0, 1])``
- etc.
The functions are coded by triples of powers of *x*, *y*, *z*, namely
``[1, 2, 3]`` corresponds to :math:`xy^{2}z^{3}`.
Args:
l : int
Angular momentum value
Returns:
out ((l+1)*(l+2)/2, 3) : numpy.array
Array of monomial powers, where columns correspond to *x*, *y*, and *z*
respectively and rows correspond to functions.
'''
ncart = (l + 1) * (l + 2) // 2
out = np.zeros((ncart, 3), dtype=np.int32)
index = 0
for i in range(l + 1):
for j in range(i + 1):
out[index, :] = [l - i, i - j, j]
index += 1
return out
[docs]def zlmtoxyz(l):
'''
Generates the expansion coefficients of the real spherical harmonics in
terms of products of cartesian components. Method based on [Schelegel1995]_
.. [Schelegel1995] Schlegel, H. B., & Frisch, M. J. (1995). "Transformation between
Cartesian and pure spherical harmonic Gaussians". International Journal
of Quantum Chemistry, 54(2), 83–87.
`doi:10.1002/qua.560540202 <http:www.dx.doi.org/10.1002/qua.560540202>`_
Args:
l : int
Angular momentum value
Returns:
out :math:`((l+1)*(l+2)/2, 2*l + 1)` : numpy.array
Expansion coefficients of real spherical harmonics in terms of
cartesian gaussians
'''
ncart = (l + 1) * (l + 2) // 2
nspher = 2 * l + 1
out = np.zeros((ncart, nspher))
cartc = xyzlist(l)
m = 0
for icart in range(ncart):
kx = cartc[icart, 0]
ky = cartc[icart, 1]
kz = cartc[icart, 2]
if kx + ky % 2 == 1:
continue
j = (kx + ky) / 2
tmpc = factorial2(2 * kx - 1) * factorial2(2 * ky - 1) * factorial2(2 * kz - 1) / factorial2(2 * l - 1)
tmpc = np.sqrt(tmpc) / factorial(j)
for i in range(l // 2 + 1):
if j > i:
continue
if kx % 2 == 1:
continue
k = kx // 2
if k > j:
continue
tmpi = tmpc * factorial2(2 * l - 2 * i - 1) / factorial(l - 2 * i) / factorial(i - j) / 2.0**i
if i % 2 == 1:
tmpi = -tmpi
out[icart, l] += binom(j, k) * tmpi
for m in range(1, l + 1):
for icart in range(ncart):
kx = cartc[icart, 0]
ky = cartc[icart, 1]
kz = cartc[icart, 2]
jj = kx + ky - m
if jj % 2 == 1:
continue
if jj < 0:
continue
j = jj // 2
tmpc = factorial2(2 * kx - 1) * factorial2(2 * ky - 1) * factorial2(2 * kz - 1) / factorial2(2 * l - 1)
tmpc = np.sqrt(2.0 * tmpc * factorial(l - m) / factorial(l + m)) / factorial(j)
for i in range((l - m) // 2 + 1):
if j > i:
continue
tmpi = tmpc * factorial2(2 * l - 2 * i - 1) / factorial(l - m - 2 * i) / factorial(i - j) / 2.0**i
if i % 2 == 1:
tmpi = -tmpi
for k in range(j + 1):
kk = kx - 2 * k
if kk < 0 or kk > m:
continue
tmpk = tmpi * binom(j, k) * binom(m, kk)
kkk = m - kk
if kkk % 4 == 0:
out[icart, m + l] += tmpk
elif kkk % 4 == 1:
out[icart, -m + l] += tmpk
elif kkk % 4 == 2:
out[icart, m + l] -= tmpk
else:
out[icart, -m + l] -= tmpk
return out