import numpy as np
from esipy.tools import find_dis, find_di, find_di_no, find_lis, find_ns, find_distances, av1245_pairs, wf_type
from esipy import mci as _mci
########## Iring ###########
# Computing the Iring (Restricted and Unrestricted)
[docs]def compute_iring(arr, aom):
"""
Calculation of the Iring aromaticity index.
:param arr: Contains the indices defining the ring connectivity.
:type arr: list of int
:param aom: Specifies the Atomic Overlap Matrices (AOMs) in the MO basis.
:type aom: list of matrices
:returns: The Iring for the given ring connectivity.
:rtype: float
"""
product = np.identity(aom[0].shape[0])
for i in arr:
product = np.dot(product, aom[i - 1])
iring = (2 ** (len(arr) - 1)) * np.trace(product)
return iring
[docs]def compute_iring_no(arr, aom):
"""
Calculation of the Iring aromaticity index for correlated wavefunctions.
:param arr: Contains the indices defining the ring connectivity.
:type arr: list of int
:param aom: Specifies the Atomic Overlap Matrices (AOMs) in the MO basis.
:type aom: list of matrices
:returns: The Iring for the given ring connectivity.
:rtype: float
"""
aom, occ = aom
product = np.identity(aom[0].shape[0])
for i in arr:
product = np.dot(product, np.dot(occ, aom[i - 1]))
return np.trace(product)
########### MCI ###########
[docs]def sequential_mci(arr, aom, partition):
"""
Computes the MCI sequentially by computing the Iring without storing the permutations.
Default option if no number of cores is specified.
"""
# Use esipy.mci module
return _mci.compute_mci(arr, aom, partition=partition, n_cores=1)
[docs]def sequential_mci_no(arr, aom, partition):
"""
Computes the MCI for correlated wavefunctions sequentially by computing the Iring without storing the permutations.
"""
# Use esipy.mci module
return _mci.compute_mci(arr, aom, partition=partition, n_cores=1)
[docs]def multiprocessing_mci(arr, aom, ncores, partition):
"""
Computes the MCI by generating all the permutations for a later distribution along the specified number of cores.
"""
# Use esipy.mci module
return _mci.compute_mci(arr, aom, partition=partition, n_cores=ncores)
[docs]def multiprocessing_mci_no(arr, aom, ncores, partition):
"""
Computes the MCI for correlated wavefunctions by generating all the permutations for a later distribution along the specified number of cores.
"""
# Use esipy.mci module
return _mci.compute_mci(arr, aom, partition=partition, n_cores=ncores)
########### AV1245 ###########
# Calculation of the AV1245 index (Restricted and Unrestricted)
[docs]def compute_av1245(arr, aom, partition):
"""
Computes the AV1245 and AVmin indices. Not available for rings smaller than 6 members.
:param arr: Contains the indices defining the ring connectivity.
:type arr: list of int
:param aom: Specifies the Atomic Overlap Matrices (AOMs) in the MO basis.
:type aom: list of matrices
:param partition: Specifies the atom-in-molecule partition scheme. Options include 'mulliken', 'lowdin', 'meta-lowdin', 'nao', and 'iao'.
:type partition: str
:returns: The AV1245 index, the AVmin index, and each of the AV1245 in a list for the output, respectively.
:rtype: tuple
"""
products = []
for cp in av1245_pairs(arr.copy()):
product = sequential_mci(list(cp), aom, partition)
products.append(1000 * product / 3)
av1245_value = np.mean(products)
avmin_value = min(products, key=abs)
return av1245_value, avmin_value, products
[docs]def compute_av1245_no(arr, aom, partition):
"""
Computes the AV1245 and AVmin indices for correlated wavefunctions. Not available for rings smaller than 6 members.
:param arr: Contains the indices defining the ring connectivity.
:type arr: list of int
:param aom: Specifies the Atomic Overlap Matrices (AOMs) in the MO basis.
:type aom: list of matrices
:param partition: Specifies the atom-in-molecule partition scheme. Options include 'mulliken', 'lowdin', 'meta-lowdin', 'nao', and 'iao'.
:type partition: str
:returns: The AV1245 index, the AVmin index, and each of the AV1245 in a list for the output, respectively.
:rtype: tuple
"""
products = []
for cp in av1245_pairs(arr):
product = sequential_mci_no(list(cp), aom, partition)
products.append(1000 * product / 3)
av1245_value = np.mean(products)
avmin_value = min(products, key=abs)
return av1245_value, avmin_value, products
########### PDI ###########
# Calculation of the PDI (Restricted and Unrestricted)
[docs]def compute_pdi(arr, aom):
"""
Computes the PDI for the given 6-membered ring connectivity. Only computed for rings n=6.
:param arr: Contains the indices defining the ring connectivity.
:type arr: list of int
:param aom: Specifies the Atomic Overlap Matrices (AOMs) in the MO basis.
:type aom: list of matrices
:returns: The PDI value and each of the DIs in para position.
:rtype: tuple
"""
if len(arr) == 6:
pdi_a = 2 * np.trace(np.dot(aom[arr[0] - 1], aom[arr[3] - 1]))
pdi_b = 2 * np.trace(np.dot(aom[arr[1] - 1], aom[arr[4] - 1]))
pdi_c = 2 * np.trace(np.dot(aom[arr[2] - 1], aom[arr[5] - 1]))
pdi_value = (pdi_a + pdi_b + pdi_c) / 3
return pdi_value, [pdi_a, pdi_b, pdi_c]
else:
return None
[docs]def compute_pdi_no(arr, aom):
"""
Computes the PDI for the given 6-membered ring connectivity. Only computed for rings n=6.
:param arr: Contains the indices defining the ring connectivity.
:type arr: list of int
:param aom: Specifies the Atomic Overlap Matrices (AOMs) in the MO basis.
:type aom: list of matrices
:returns: The PDI value and each of the DIs in para position.
:rtype: tuple
"""
aom, occ = aom
if len(arr) == 6:
i1, i2, i3, i4, i5, i6 = arr[0] - 1, arr[1] - 1, arr[2] - 1, arr[3] - 1, arr[4] - 1, arr[5] - 1
pdi_a = 2 * np.trace(np.linalg.multi_dot((occ ** (1 / 2), aom[i1], occ ** (1 / 2), aom[i4])))
pdi_b = 2 * np.trace(np.linalg.multi_dot((occ ** (1 / 2), aom[i2], occ ** (1 / 2), aom[i5])))
pdi_c = 2 * np.trace(np.linalg.multi_dot((occ ** (1 / 2), aom[i3], occ ** (1 / 2), aom[i6])))
pdi_value = (pdi_a + pdi_b + pdi_c) / 3
return pdi_value, [pdi_a, pdi_b, pdi_c]
else:
return None
########### FLU ###########
# Calculation of the FLU (Restricted and Unrestricted)
[docs]def find_flurefs(partition=None):
"""
Sets the reference of the FLU index based on the provided partition.
The available options are "CC" from benzene, "CN" from pyridine,
"BN" from borazine, "NN" from pyridazine and "CS" from thiophene,
all obtained from optimized and single-point calculations at HF/6-31G(d)
level of theory.
:param partition: Specifies the atom-in-molecule partition scheme.
Options include 'mulliken', 'lowdin', 'meta-lowdin', 'nao', and 'iao'.
:type partition: str
:returns: Contains the reference DI for each bond.
:rtype: dict
"""
if partition == "qtaim":
return {"CC": 1.3994, "CN": 1.1957, "BN": 0.3935, "NN": 1.5253, "CS": 1.2366, "CO": 0.8757}
elif partition == "mulliken":
return {"CC": 1.4536, "CN": 1.4051, "BN": 1.1042, "NN": 1.2831, "CS": 1.1013, "CO": 0.9840}
elif partition == "lowdin":
return {"CC": 1.4910, "CN": 1.5864, "BN": 1.3165, "NN": 1.5277, "CS": 1.2213, "CO": 1.2368}
elif partition == "meta-lowdin":
return {"CC": 1.4392, "CN": 1.4521, "BN": 1.2048, "NN": 1.3700, "CS": 1.1455, "CO": 1.0802}
elif partition == "nao":
return {"CC": 1.4347, "CN": 1.4115, "BN": 0.9239, "NN": 1.3708, "CS": 1.1625, "CO": 0.9912}
elif partition == "iao":
return {"CC": 1.4378, "CN": 1.4386, "BN": 1.1631, "NN": 1.3606, "CS": 1.1433, "CO": 1.0708}
else:
return {"CC": 1.4378, "CN": 1.4386, "BN": 1.1631, "NN": 1.3606, "CS": 1.1433, "CO": 1.0708}
[docs]def compute_flu(arr, molinfo, aom, flurefs=None, partition=None):
"""
Computes the FLU index.
:param arr: Contains the indices defining the ring connectivity.
:type arr: list of int
:param molinfo: Contains the molecular information.
:type molinfo: dict
:param aom: Specifies the Atomic Overlap Matrices (AOMs) in the MO basis.
:type aom: list of matrices
:param flurefs: User-provided references for the FLU index.
:type flurefs: dict, optional
:param partition: Specifies the atom-in-molecule partition scheme. Options include 'mulliken', 'lowdin', 'meta-lowdin', 'nao', and 'iao'.
:type partition: str, optional
:returns: The FLU value for the given ring connectivity.
:rtype: float
"""
flu_value, flu_polar = 0, 0
symbols = molinfo["symbols"]
atom_symbols = [str(symbols[int(i) - 1]).strip().capitalize() for i in arr]
bond_types = ["".join(sorted([atom_symbols[i], atom_symbols[(i + 1) % len(arr)]]))
for i in range(len(arr))]
# Setting and update of the reference values
flu_refs = find_flurefs(partition)
if flurefs is not None:
flu_refs.update(flurefs)
dis = find_dis(arr, aom)
lis = find_lis(arr, aom)
ns = find_ns(arr, aom)
for i in range(len(arr)):
if bond_types[i] not in flu_refs:
print(f" | No parameters found for bond type {bond_types[i]}")
return None
flu_deloc = (dis[i] - flu_refs[bond_types[i]]) / flu_refs[bond_types[i]]
a_to_b = dis[i] / 2 * (ns[i] - lis[i])
b_to_a = dis[i] / 2 * (ns[(i + 1) % len(arr)] - lis[(i + 1) % len(arr)])
flu_polar = a_to_b / b_to_a
if flu_polar < 1:
flu_polar = 1 / flu_polar
flu_value += float(flu_deloc * flu_polar) ** 2
return flu_value / len(arr)
########### BOA ###########
# Calculation of the BOA (Restricted and Unrestricted)
[docs]def compute_boa(arr, aom):
"""
Computes the BOA and BOA\_c indices.
:param arr: Contains the indices defining the ring connectivity.
:type arr: list of int
:param aom: Specifies the Atomic Overlap Matrices (AOMs) in the MO basis.
:type aom: list of matrices
:returns: Contains the BOA and the BOA\_c indices, respectively.
:rtype: tuple
"""
n1 = len([i for i in arr if i % 2 != 0])
n2 = len([i for i in arr if i % 2 == 1])
sum_odd = sum(find_di(aom, arr[i - 1], arr[i]) for i in range(0, len(arr), 2))
sum_even = sum(find_di(aom, arr[i + 1], arr[i]) for i in range(0, len(arr) - 1, 2))
boa = abs(sum_odd / n1 - sum_even / n2)
boa_c = 0
for i in range(len(arr)):
diff_di = abs(
find_di(aom, arr[i - 1], arr[i]) - find_di(aom, arr[(i + 1) % len(arr) - 1], arr[(i + 1) % len(arr)]))
boa_c += diff_di / len(arr)
return boa, boa_c
[docs]def compute_boa_no(arr, aom):
"""
Computes the BOA and BOA\_c indices for correlated wavefunctions.
:param arr: Contains the indices defining the ring connectivity.
:type arr: list of int
:param aom: Specifies the Atomic Overlap Matrices (AOMs) in the MO basis.
:type aom: list of matrices
:returns: Contains the BOA and the BOA\_c indices, respectively.
:rtype: tuple
"""
n1 = len([i for i in arr if i % 2 != 0])
n2 = len([i for i in arr if i % 2 == 1])
sum_odd = sum(find_di_no(aom, arr[i - 1], arr[i]) for i in range(0, len(arr), 2))
sum_even = sum(find_di_no(aom, arr[i + 1], arr[i]) for i in range(0, len(arr) - 1, 2))
boa = abs(sum_odd / n1 - sum_even / n2)
boa_c = 0
for i in range(len(arr)):
diff_di = abs(
find_di_no(aom, arr[i - 1], arr[i]) - find_di_no(aom, arr[(i + 1) % len(arr) - 1], arr[(i + 1) % len(arr)]))
boa_c += diff_di / len(arr)
return boa, boa_c
######## GEOMETRIC INDICES ########
# Calculation of the HOMA and/or HOMER indices (Restricted and Unrestricted)
[docs]def compute_homer(arr, molinfo, homerrefs=None):
"""
Computes the HOMER index.
:param arr: Contains the indices defining the ring connectivity.
:type arr: list of int
:param molinfo: Contains the molecular information.
:type molinfo: dict
:param homerrefs: User-provided references for the HOMER index.
:type homerrefs: dict, optional
:returns: HOMER value for the given ring connectivity.
:rtype: float
"""
geom = molinfo["geom"]
refs = {
"CC": {"alpha": 950.74, "r_opt": 1.437},
"CN": {"alpha": 506.43, "r_opt": 1.390},
"NN": {"alpha": 187.36, "r_opt": 1.375},
"CO": {"alpha": 164.96, "r_opt": 1.379}
}
if homerrefs is not None:
refs.update(homerrefs)
symbols = molinfo["symbols"]
atom_symbols = [str(symbols[int(i) - 1]).strip().capitalize() for i in arr]
bond_types = ["".join(sorted([atom_symbols[arr[i] - 1], atom_symbols[arr[(i + 1) % len(arr)] - 1]]))
for i in range(len(arr))]
for i in range(len(arr)):
if bond_types[i] not in refs:
print(f"No parameters found for bond type {bond_types[i]}")
return None
alpha = refs[bond_types[i]]["alpha"]
r_opt = refs[bond_types[i]]["r_opt"]
distances = find_distances(arr, geom)
diff = np.mean([r_opt - distances[i] for i in range(len(arr))])
homer_value = 1 - alpha * diff ** 2
return homer_value
[docs]def compute_homa(arr, molinfo, homarefs=None):
"""
Computes the HOMA index.
:param arr: Contains the indices defining the ring connectivity.
:type arr: list of int
:param molinfo: Contains the molecular information.
:type molinfo: dict
:param homarefs: User-provided references for the HOMA index.
:type homarefs: dict, optional
:returns: HOMA value for the given ring connectivity.
:rtype: float
"""
refs = {
"CC": {"n_opt": 1.590, "c": 0.1702, "r1": 1.467},
"CN": {"n_opt": 1.589, "c": 0.2828, "r1": 1.465},
"CO": {"n_opt": 1.602, "c": 0.2164, "r1": 1.367},
"CP": {"n_opt": 1.587, "c": 0.2510, "r1": 1.814},
"CS": {"n_opt": 1.584, "c": 0.2828, "r1": 1.807},
"NN": {"n_opt": 1.590, "c": 0.2395, "r1": 1.420},
"NO": {"n_opt": 1.586, "c": 0.3621, "r1": 1.415},
"CSe": {"n_opt": 1.590, "c": 0.2970, "r1": 1.959}, # n_opt taken as C-C
"BB": {"n_opt": 1.590, "c": 0.2510, "r1": 1.814}, # n_opt taken as C-C
"BC": {"n_opt": 1.590, "c": 0.1752, "r1": 1.647}, # n_opt taken as C-C
"BN": {"n_opt": 1.590, "c": 0.2900, "r1": 1.564}, # n_opt taken as C-C
"alpha": 257.7, "r_opt": 1.388
}
if homarefs is not None:
refs.update(homarefs)
geom = molinfo["geom"]
if geom is None:
return None
symbols = molinfo["symbols"]
atom_symbols = [str(symbols[int(i) - 1]).strip().capitalize() for i in arr]
bond_types = ["".join(sorted([atom_symbols[i], atom_symbols[(i + 1) % len(arr)]]))
for i in range(len(arr))]
for i in range(len(arr)):
if bond_types[i] not in refs:
print(f" | No parameters found for bond type {bond_types[i]}")
return None
distances = find_distances(arr, geom)
alpha = refs["alpha"]
r_opt = refs["r_opt"]
ravs, bonds = [], []
for i in range(len(arr)):
c = refs[bond_types[i]]["c"]
r1 = refs[bond_types[i]]["r1"]
bond = np.exp((r1 - distances[i]) / c)
rn = 1.467 - 0.1702 * np.log(bond)
ravs.append(rn)
rav = sum(ravs) / len(arr)
if np.mean(rav) > r_opt:
EN = alpha * (r_opt - rav) ** 2
else:
EN = -alpha * (r_opt - rav) ** 2
GEO = 0
for i in range(len(arr)):
GEO += (rav - ravs[i]) ** 2
GEO = GEO * alpha / len(arr)
homa_value = 1 - (EN + GEO)
return homa_value, EN, GEO
# Calculation of the BLA (Restricted and Unrestricted)
[docs]def compute_bla(arr, molinfo):
"""
Computes the BLA and BLA\_c indices.
:param arr: Contains the indices defining the ring connectivity.
:type arr: list of int
:param molinfo: Contains the molecular information.
:type molinfo: dict
:returns: Contains the BLA and the BLA\_c indices, respectively.
:rtype: tuple
"""
if molinfo["geom"] is None:
return None
distances = find_distances(arr, molinfo["geom"])
sum1 = sum(distances[i] for i in range(0, len(arr), 2))
sum2 = sum(distances[i] for i in range(1, len(arr), 2))
bla = abs(sum1 / (len(arr) // 2) - sum2 / (len(arr) - len(arr) // 2))
bla_c = 0
for i in range(len(arr)):
bla_c += abs(distances[i] - distances[(i + 1) % len(distances)]) / len(distances)
return bla, bla_c
