diff --git a/.circleci/config.yml b/.circleci/config.yml index 301e7d7e..328109e6 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -46,9 +46,9 @@ jobs: conda install -y -c conda-forge openbabel=3.1.1 # install dependencies conda install -y -c conda-forge xtb=6.7.1 - conda install -y -c conda-forge crest=3.0.2 - pip install ase - conda install pytorch torchvision torchaudio cpuonly -c pytorch + conda install -y -c conda-forge crest=2.12 + pip install ase==3.24.0 + pip install torch==2.5.1 torchvision==0.20.1 torchaudio==2.5.1 --index-url https://download.pytorch.org/whl/cpu pip install torchani==2.2.3 # install the dependencies of AQME with pip pip install . diff --git a/aqme/qdescp.py b/aqme/qdescp.py index e56dfd34..1c88e1ef 100644 --- a/aqme/qdescp.py +++ b/aqme/qdescp.py @@ -863,7 +863,7 @@ def run_sp_xtb(self, file, xyz_file, charge, mult, name, destination): "--uhf", str(mult - 1), "--etemp", - str(self.args.qdescp_temp), + '5000', "-P", "1", ] # for FOD diff --git a/aqme/qdescp_utils.py b/aqme/qdescp_utils.py index 14c6a4d7..7df17ee7 100644 --- a/aqme/qdescp_utils.py +++ b/aqme/qdescp_utils.py @@ -417,7 +417,7 @@ def extract_scc_energy(lines, filename): if "SCC energy" in line: return float(line.split()[3]) - self.args.log.write(f"x WARNING! Could not find SCC energy value in the file: {filename}") + self.args.log.write(f"x WARNING! Could not find SCC energy value in the file: {os.path.basename(filename)}. Some CDFT-based descriptors will be missing.") return None try: @@ -443,16 +443,17 @@ def extract_scc_energy(lines, filename): scc_energy_Nplus1 = extract_scc_energy(data3, file_Nplus1) scc_energy_Nplus2 = extract_scc_energy(data4, file_Nplus2) - if None in [scc_energy, scc_energy_Nminus1, scc_energy_Nminus2, scc_energy_Nplus1, scc_energy_Nplus2]: - self.args.log.write(f"x WARNING! Missing SCC energy in one or more files. Cannot calculate global CDFT descriptors.") - return None - # Convert SCC energies to Hartree - scc_energy *= Hartree - scc_energy_Nminus1 *= Hartree - scc_energy_Nminus2 *= Hartree - scc_energy_Nplus1 *= Hartree - scc_energy_Nplus2 *= Hartree + if scc_energy is not None: + scc_energy *= Hartree + if scc_energy_Nminus1 is not None: + scc_energy_Nminus1 *= Hartree + if scc_energy_Nminus2 is not None: + scc_energy_Nminus2 *= Hartree + if scc_energy_Nplus1 is not None: + scc_energy_Nplus1 *= Hartree + if scc_energy_Nplus2 is not None: + scc_energy_Nplus2 *= Hartree # Initialize variables delta_SCC_IP, delta_SCC_EA, electrophilicity_index = None, None, None @@ -466,60 +467,66 @@ def extract_scc_energy(lines, filename): Electrophilic_descriptor, w_cubic = None, None # Calculate Global CDFT descriptors - delta_SCC_IP = round(((scc_energy_Nminus1 - corr_xtb) - scc_energy),4) - delta_SCC_EA = round((scc_energy - (scc_energy_Nplus1 + corr_xtb)),4) - chemical_hardness = round((delta_SCC_IP - delta_SCC_EA), 4) - chemical_potential = round(-(delta_SCC_IP + delta_SCC_EA) / 2, 4) - electrophilicity_index = (chemical_potential**2)/(2*chemical_hardness) - mulliken_electronegativity = round(-chemical_potential, 4) - electrofugality = round(-delta_SCC_EA + electrophilicity_index, 4) - nucleofugality = round(delta_SCC_IP + electrophilicity_index, 4) - electrodonating_maximum_electron_flow = round((-(chemical_potential/chemical_hardness)),4) - electrodonating_chemical_potential = round(((1/4)*((-3*delta_SCC_IP) - delta_SCC_EA)),4) - lectrodonating_maximum_electron_flow = round((-(electrodonating_chemical_potential/chemical_hardness)),4) - Vertical_second_IP = round((((scc_energy_Nminus2 - scc_energy_Nminus1) - corr_xtb)), 4) - Vertical_second_EA = round((((scc_energy_Nplus1 - scc_energy_Nplus2) + corr_xtb)), 4) - hyper_hardness = round((-((0.5) * (delta_SCC_IP + delta_SCC_EA - Vertical_second_IP - Vertical_second_EA))), 4) - - if chemical_hardness != 0: + if None not in [scc_energy_Nminus1,scc_energy]: + delta_SCC_IP = round(((scc_energy_Nminus1 - corr_xtb) - scc_energy),4) + if None not in [scc_energy_Nplus1,scc_energy]: + delta_SCC_EA = round((scc_energy - (scc_energy_Nplus1 + corr_xtb)),4) + if None not in [delta_SCC_IP,delta_SCC_EA]: + chemical_hardness = round((delta_SCC_IP - delta_SCC_EA), 4) + chemical_potential = round(-(delta_SCC_IP + delta_SCC_EA) / 2, 4) + electrophilicity_index = (chemical_potential**2)/(2*chemical_hardness) + mulliken_electronegativity = round(-chemical_potential, 4) + electrofugality = round(-delta_SCC_EA + electrophilicity_index, 4) + nucleofugality = round(delta_SCC_IP + electrophilicity_index, 4) + electrodonating_maximum_electron_flow = round((-(chemical_potential/chemical_hardness)),4) + electrodonating_chemical_potential = round(((1/4)*((-3*delta_SCC_IP) - delta_SCC_EA)),4) + electrodonating_maximum_electron_flow = round((-(electrodonating_chemical_potential/chemical_hardness)),4) + if None not in [scc_energy_Nminus2,scc_energy_Nminus1]: + Vertical_second_IP = round((((scc_energy_Nminus2 - scc_energy_Nminus1) - corr_xtb)), 4) + if None not in [scc_energy_Nplus1,scc_energy_Nplus2]: + Vertical_second_EA = round((((scc_energy_Nplus1 - scc_energy_Nplus2) + corr_xtb)), 4) + if None not in [delta_SCC_IP,delta_SCC_EA,Vertical_second_IP,Vertical_second_EA]: + hyper_hardness = round((-((0.5) * (delta_SCC_IP + delta_SCC_EA - Vertical_second_IP - Vertical_second_EA))), 4) + + if chemical_hardness is not None and chemical_hardness != 0: chemical_softness = round(1 / chemical_hardness, 4) electrodonating_power_index = round(((delta_SCC_IP + 3 * delta_SCC_EA)**2) / (8 * chemical_hardness), 4) electroaccepting_power_index = round(((3 * delta_SCC_IP + delta_SCC_EA)**2) / (8 * chemical_hardness), 4) intrinsic_reactivity_index = round((delta_SCC_IP + delta_SCC_EA) / chemical_hardness, 4) - Global_hypersoftness = round((hyper_hardness / ((chemical_hardness) ** 3)), 4) + if hyper_hardness is not None: + Global_hypersoftness = round((hyper_hardness / ((chemical_hardness) ** 3)), 4) if electroaccepting_power_index != 0: nucleophilicity_index = round(10 / electroaccepting_power_index, 4) - net_electrophilicity = round((electrodonating_power_index - electroaccepting_power_index), 4) - - # For lectrophilic descriptor calculations - A = ((scc_energy_Nplus1 - scc_energy) + corr_xtb) - c = (Vertical_second_IP - (2 * delta_SCC_IP) + A) / ((2 * Vertical_second_IP) - delta_SCC_IP - A) - a = -((delta_SCC_IP + A) / 2) + (((delta_SCC_IP - A) / 2) * c) - b = ((delta_SCC_IP - A) / 2) - (((delta_SCC_IP + A) / 2) * c) - Gamma = (-3 * c) * (b - (a * c)) - Eta = 2 * (b - (a * c)) - chi = -a - Mu = a - - discriminant = Eta ** 2 - (2 * Gamma * Mu) # Checking the square root - if discriminant < 0: - self.args.log.write(f"x WARNING! Negative discriminant encountered, skipping Electrophilic descriptor calculation.") - else: - inter_phi = math.sqrt(discriminant) - Phi = inter_phi - Eta - Electrophilic_descriptor = round(((chi * (Phi / Gamma)) - (((Phi / Gamma) ** 2) * ((Eta / 2) + (Phi / 6)))), 4) + if None not in [electrodonating_power_index,electroaccepting_power_index]: + net_electrophilicity = round((electrodonating_power_index - electroaccepting_power_index), 4) + + # For electrophilic descriptor calculations + if None not in [scc_energy_Nplus1,scc_energy,Vertical_second_IP,delta_SCC_IP]: + A = ((scc_energy_Nplus1 - scc_energy) + corr_xtb) + c = (Vertical_second_IP - (2 * delta_SCC_IP) + A) / ((2 * Vertical_second_IP) - delta_SCC_IP - A) + a = -((delta_SCC_IP + A) / 2) + (((delta_SCC_IP - A) / 2) * c) + b = ((delta_SCC_IP - A) / 2) - (((delta_SCC_IP + A) / 2) * c) + Gamma = (-3 * c) * (b - (a * c)) + Eta = 2 * (b - (a * c)) + chi = -a + Mu = a + + discriminant = Eta ** 2 - (2 * Gamma * Mu) # Checking the square root + if discriminant >= 0: + inter_phi = math.sqrt(discriminant) + Phi = inter_phi - Eta + Electrophilic_descriptor = round(((chi * (Phi / Gamma)) - (((Phi / Gamma) ** 2) * ((Eta / 2) + (Phi / 6)))), 4) # For cubic electrophilicity index - Gamma_cubic = 2 * delta_SCC_IP - Vertical_second_IP - delta_SCC_EA - Eta_cubic = delta_SCC_IP - delta_SCC_EA + if None not in [delta_SCC_IP,Vertical_second_IP,delta_SCC_EA]: + Gamma_cubic = 2 * delta_SCC_IP - Vertical_second_IP - delta_SCC_EA + Eta_cubic = delta_SCC_IP - delta_SCC_EA if Eta_cubic != 0: Mu_cubic = (1 / 6) * ((-2 * delta_SCC_EA) - (5 * delta_SCC_IP) + Vertical_second_IP) w_cubic = round(((Mu_cubic ** 2) / (2 * Eta_cubic)) * (1 + ((Mu_cubic / (3 * (Eta_cubic) ** 2)) * Gamma_cubic)), 4) - else: - self.args.log.write(f"x WARNING! Eta_cubic is zero, skipping cub. electrophilicity idx calculation.") # Return the calculated descriptors cdft_descriptors = { @@ -545,7 +552,7 @@ def extract_scc_energy(lines, filename): "cub. electrophilicity idx": w_cubic, "max. electron flow": electrodonating_maximum_electron_flow, "Electrodon. Chem. potential": electrodonating_chemical_potential, - "Electrodon. max. electron flow": lectrodonating_maximum_electron_flow + "Electrodon. max. electron flow": electrodonating_maximum_electron_flow } return cdft_descriptors @@ -558,7 +565,13 @@ def calculate_local_CDFT_descriptors(file_fukui, cdft_descriptors,self): with open(file_fukui, "r") as f: data = f.readlines() + # Initialize variables f_pos, f_negs, f_rads = [], [], [] + dual_descriptor,s_pos,s_negs,s_rads = None,None,None,None + Relative_nucleophilicity,Relative_electrophilicity,Grand_canonical_dual_descriptor = None,None,None + w_pos,w_negs,w_rads,Multiphilic_descriptor = None,None,None,None + Nu_pos,Nu_negs,Nu_rads = None,None,None + start, end = None, None for i, line in enumerate(data): @@ -578,50 +591,46 @@ def calculate_local_CDFT_descriptors(file_fukui, cdft_descriptors,self): f_rads.append(f_rad) except ValueError: continue - else: - self.args.log.write("WARNING: Fukui data not found in the file. Please check the '.fukui' file.") - return None - if not f_pos or not f_negs or not f_rads: - self.args.log.write("WARNING: Fukui data lists are empty. Please check the '.fukui' file.") - return None - - if None in [cdft_descriptors]: - self.args.log.write("x WARNING! Missing required CDFT descriptors (Softness, Hypersoftness, Electrophil. idx or Nucleophilicity idx).") - return None + if f_pos == [] or f_negs == [] or f_rads == []: + self.args.log.write("x WARNING! Fukui indices did not generate, please check the '.fukui' file.") chemical_softness = cdft_descriptors.get("Softness") Global_hypersoftness = cdft_descriptors.get("Hypersoftness") electrophilicity_index = cdft_descriptors.get("Electrophil. idx") nucleophilicity_index = cdft_descriptors.get("Nucleophilicity idx") - if None in [chemical_softness, Global_hypersoftness, electrophilicity_index, nucleophilicity_index]: - self.args.log.write("x WARNING! Missing required CDFT descriptors (Softness, Hypersoftness, Electrophil. idx or Nucleophilicity idx).") - return None - - # Calculating local descriptors: - # 1) dual descrip. - dual_descriptor = [round(f_po - f_neg, 4) for f_po, f_neg in zip(f_pos, f_negs)] - # 2) softness+, softness- and softness0 - s_pos = [round(chemical_softness * f_po, 4) for f_po in f_pos] - s_negs = [round(chemical_softness * f_neg, 4) for f_neg in f_negs] - s_rads = [round(chemical_softness * f_rad, 4) for f_rad in f_rads] - # 3) Rel. nucleophilicity - Relative_nucleophilicity = [round(s_neg / s_po, 4) if s_po != 0 else None for s_neg, s_po in zip(s_negs, s_pos)] - # 4) Rel. electrophilicity - Relative_electrophilicity = [round(s_po / s_neg, 4) if s_neg != 0 else None for s_neg, s_po in zip(s_negs, s_pos)] - # 5) GC Dual Descrip. - Grand_canonical_dual_descriptor = [round(Global_hypersoftness * dual, 4) for dual in dual_descriptor] - # 6) softness+, softness- and softness0 - w_pos = [round(electrophilicity_index * f_po, 4) for f_po in f_pos] - w_negs = [round(electrophilicity_index * f_neg, 4) for f_negs in f_negs] - w_rads = [round(electrophilicity_index * f_rad, 4) for f_rad in f_rads] - # 7) softness+, softness- and softness0 - Multiphilic_descriptor = [round(electrophilicity_index * dual, 4) for dual in dual_descriptor] - # 8) softness+, softness- and softness0 - Nu_pos = [round(nucleophilicity_index * f_po, 4) for f_po in f_pos] - Nu_negs = [round(nucleophilicity_index * f_neg, 4) for f_neg in f_negs] - Nu_rads = [round(nucleophilicity_index * f_rad, 4) for f_rad in f_rads] + # Calculating local descriptors + if chemical_softness is not None and f_pos != [] and f_negs != [] and f_rads != []: + # 1) dual descrip. + dual_descriptor = [round(f_po - f_neg, 4) for f_po, f_neg in zip(f_pos, f_negs)] + # 2) softness+, softness- and softness0 + s_pos = [round(chemical_softness * f_po, 4) for f_po in f_pos] + s_negs = [round(chemical_softness * f_neg, 4) for f_neg in f_negs] + s_rads = [round(chemical_softness * f_rad, 4) for f_rad in f_rads] + # 3) Rel. nucleophilicity + Relative_nucleophilicity = [round(s_neg / s_po, 4) if s_po != 0 else None for s_neg, s_po in zip(s_negs, s_pos)] + # 4) Rel. electrophilicity + Relative_electrophilicity = [round(s_po / s_neg, 4) if s_neg != 0 else None for s_neg, s_po in zip(s_negs, s_pos)] + + if Global_hypersoftness is not None and dual_descriptor is not None: + # 5) GC Dual Descrip. + Grand_canonical_dual_descriptor = [round(Global_hypersoftness * dual, 4) for dual in dual_descriptor] + + if electrophilicity_index is not None and f_pos != [] and f_negs != [] and f_rads != []: + # 6) softness+, softness- and softness0 + w_pos = [round(electrophilicity_index * f_po, 4) for f_po in f_pos] + w_negs = [round(electrophilicity_index * f_neg, 4) for f_neg in f_negs] + w_rads = [round(electrophilicity_index * f_rad, 4) for f_rad in f_rads] + if dual_descriptor is not None: + # 7) softness+, softness- and softness0 + Multiphilic_descriptor = [round(electrophilicity_index * dual, 4) for dual in dual_descriptor] + + if nucleophilicity_index is not None and f_pos != [] and f_negs != [] and f_rads != []: + # 8) softness+, softness- and softness0 + Nu_pos = [round(nucleophilicity_index * f_po, 4) for f_po in f_pos] + Nu_negs = [round(nucleophilicity_index * f_neg, 4) for f_neg in f_negs] + Nu_rads = [round(nucleophilicity_index * f_rad, 4) for f_rad in f_rads] localDescriptors = { "fukui+": f_pos, @@ -1598,7 +1607,10 @@ def update_atom_props_json(sorted_indices,match_names,atom_props,json_data,prefi idx_xtb = atom_idx for prop in atom_props: try: - json_data[f'{match_name}_{prop}'] = json_data[prop][idx_xtb] + if json_data[prop] is not None: + json_data[f'{match_name}_{prop}'] = json_data[prop][idx_xtb] + else: + json_data[f'{match_name}_{prop}'] = None if f'{match_name}_' not in prefixes_atom_prop: prefixes_atom_prop.append(f'{match_name}_') except (KeyError,IndexError): # prevents missing values @@ -1610,8 +1622,12 @@ def update_atom_props_json(sorted_indices,match_names,atom_props,json_data,prefi prop_values = [] for prop_name in match_names: prop_values.append(json_data[f'{prop_name}_{prop}']) - json_data[f'{pattern}_max_{prop}'] = max(prop_values) - json_data[f'{pattern}_min_{prop}'] = min(prop_values) + if None not in prop_values: + json_data[f'{pattern}_max_{prop}'] = max(prop_values) + json_data[f'{pattern}_min_{prop}'] = min(prop_values) + else: + json_data[f'{pattern}_max_{prop}'] = None + json_data[f'{pattern}_min_{prop}'] = None if f'{pattern}_max_{prop}' not in prefixes_atom_prop: prefixes_atom_prop.append(f'{pattern}_max_') if f'{pattern}_min_{prop}' not in prefixes_atom_prop: diff --git a/aqme/utils.py b/aqme/utils.py index 18e70874..96bd6e82 100644 --- a/aqme/utils.py +++ b/aqme/utils.py @@ -25,11 +25,11 @@ T = 298.15 obabel_version = "3.1.1" # this MUST match the meta.yaml -aqme_version = "1.7.2" +aqme_version = "1.7.1" time_run = time.strftime("%Y/%m/%d %H:%M:%S", time.localtime()) aqme_ref = f"AQME v {aqme_version}, Alegre-Requena, J. V.; Sowndarya, S.; Perez-Soto, R.; Alturaifi, T.; Paton, R. AQME: Automated Quantum Mechanical Environments for Researchers and Educators. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2023, 13, e1663 (DOI: 10.1002/wcms.1663)." xtb_version = '6.7.1' -crest_version = '3.0.2' +crest_version = '2.12' RDLogger.DisableLog("rdApp.*") @@ -1022,11 +1022,16 @@ def check_version(self, program, version_line, target_version, n_split, install_ lower_version = True if int(version_found.split('.')[0]) == int(target_version.split('.')[0]): - if len(target_version.split('.')) >= 2 and int(version_found.split('.')[1]) == int(target_version.split('.')[1]): - if len(target_version.split('.')) >= 3 \ - and int(version_found.split('.')[2]) >= int(target_version.split('.')[2]) \ - and int(version_found.split('.')[2]) < int(target_version.split('.')[2])+5: # up to four versions ahead - lower_version = False + if len(target_version.split('.')) == 2 \ + and int(version_found.split('.')[1]) >= int(target_version.split('.')[1]) \ + and int(version_found.split('.')[1]) < int(target_version.split('.')[1])+3: # up to two versions ahead + lower_version = False + + if len(target_version.split('.')) == 3 \ + and int(version_found.split('.')[1]) == int(target_version.split('.')[1]) \ + and int(version_found.split('.')[2]) >= int(target_version.split('.')[2]) \ + and int(version_found.split('.')[2]) < int(target_version.split('.')[2])+5: # up to four versions ahead + lower_version = False if version_found == '0.0.0': version_found = 'Unknown' diff --git a/docs/Misc/versions.rst b/docs/Misc/versions.rst index a96c84e3..1e267499 100644 --- a/docs/Misc/versions.rst +++ b/docs/Misc/versions.rst @@ -9,7 +9,7 @@ Version 1.7.1 [`url `__] - CSEARCH uses multithreading to ensure reproducbility acorss runs - Added T1-S0 gap (not only S0-T1) - PTB is now used to calculate HOMO-LUMO gaps, charges, dipole moments, and other descriptors - - Using updated versions of xTB & CREST + - Using xTB v6.7.1 and CREST v2.12 - Fixed bug from MORFEUS descriptors in QDESCP - Improved SDF reader compatible with SDF files from GaussView - Refactoring of the QDESCP module and its utils (faster code execution) diff --git a/docs/README.rst b/docs/README.rst index 16a3f9ff..9483d063 100644 --- a/docs/README.rst +++ b/docs/README.rst @@ -130,7 +130,7 @@ Extra requirements if xTB or CREST are used (compatible with MacOS and Linux onl .. code-block:: shell - conda install -y -c conda-forge crest=3.0.2 + conda install -y -c conda-forge crest=2.12 Extra requirements if `CMIN` is used with ANI models: