diff --git a/arc/main.py b/arc/main.py index d7fb8131c9..e8d7ae2557 100644 --- a/arc/main.py +++ b/arc/main.py @@ -81,6 +81,7 @@ class ARC(object): instead (e.g., ``opt_level``). composite_method (str, dict, Level, optional): Composite method. conformer_level (str, dict, Level, optional): Level of theory for conformer searches. + conf_generation_level (str, dict, Level, optional): Level of theory for conformer generation. opt_level (str, dict, Level, optional): Level of theory for geometry optimization. freq_level (str, dict, Level, optional): Level of theory for frequency calculations. sp_level (str, dict, Level, optional): Level of theory for single point calculations. @@ -167,6 +168,7 @@ class ARC(object): level_of_theory (str): A shortcut representing either sp//geometry levels or a composite method. composite_method (Level): Composite method. conformer_level (Level): Level of theory for conformer searches. + conf_generation_level (Level): Level of theory for conformer generation. opt_level (Level): Level of theory for geometry optimization. freq_level (Level): Level of theory for frequency calculations. sp_level (Level): Level of theory for single point calculations. @@ -242,6 +244,7 @@ def __init__(self, compute_thermo: bool = True, compute_transport: bool = False, conformer_level: Optional[Union[str, dict, Level]] = None, + conf_generation_level: Optional[Union[str, dict, Level]] = None, dont_gen_confs: List[str] = None, e_confs: float = 5.0, ess_settings: Dict[str, Union[str, List[str]]] = None, @@ -342,6 +345,7 @@ def __init__(self, self.level_of_theory = level_of_theory self.composite_method = composite_method or None self.conformer_level = conformer_level or None + self.conf_generation_level = conf_generation_level or None self.opt_level = opt_level or None self.freq_level = freq_level or None self.sp_level = sp_level or None @@ -485,6 +489,8 @@ def as_dict(self) -> dict: restart_dict['compute_transport'] = self.compute_transport if self.conformer_level is not None: restart_dict['conformer_level'] = self.conformer_level.as_dict() + if self.conf_generation_level is not None: + restart_dict['conf_generation_level'] = self.conf_generation_level.as_dict() if self.dont_gen_confs: restart_dict['dont_gen_confs'] = self.dont_gen_confs if self.ts_adapters: @@ -569,15 +575,21 @@ def execute(self) -> dict: Status dictionary indicating which species converged successfully. """ logger.info('\n') + considered_list = list() for species in self.species: if not isinstance(species, ARCSpecies): raise ValueError(f'All species must be ARCSpecies objects. Got {type(species)}') if species.is_ts: logger.info(f'Considering transition state: {species.label}') else: - logger.info(f'Considering species: {species.label}') - if species.mol is not None: - display(species.mol.copy(deep=True)) + if not species.multi_species: + logger.info(f'Considering species: {species.label}') + if species.mol is not None: + display(species.mol.copy(deep=True)) + elif species.multi_species not in considered_list: + logger.info(f'Considering species: {species.multi_species}') + considered_list.append(species.multi_species) + logger.info('\n') for rxn in self.reactions: if not isinstance(rxn, ARCReaction): @@ -587,6 +599,7 @@ def execute(self) -> dict: rxn_list=self.reactions, composite_method=self.composite_method, conformer_level=self.conformer_level, + conf_generation_level=self.conf_generation_level, opt_level=self.opt_level, freq_level=self.freq_level, sp_level=self.sp_level, @@ -963,6 +976,11 @@ def set_levels_of_theory(self): self.conformer_level = Level(repr=self.conformer_level) logger.info(f'Conformers:{default_flag} {self.conformer_level}') + if self.conf_generation_level is not None: + default_flag = '' + self.conf_generation_level = Level(repr=self.conf_generation_level) + logger.info(f'Conformers_generation:{default_flag} {self.conf_generation_level}') + if self.reactions or any([spc.is_ts for spc in self.species]): if not self.ts_guess_level: self.ts_guess_level = default_levels_of_theory['ts_guesses'] diff --git a/arc/scheduler.py b/arc/scheduler.py index 1d1ec71868..e6240e7eab 100644 --- a/arc/scheduler.py +++ b/arc/scheduler.py @@ -130,6 +130,7 @@ class Scheduler(object): project_directory (str): Folder path for the project: the input file path or ARC/Projects/project-name. composite_method (str, optional): A composite method to use. conformer_level (Union[str, dict], optional): The level of theory to use for conformer comparisons. + conf_generation_level (Union[str, dict], optional): The level of theory to use for conformer generation. opt_level (Union[str, dict], optional): The level of theory to use for geometry optimizations. freq_level (Union[str, dict], optional): The level of theory to use for frequency calculations. sp_level (Union[str, dict], optional): The level of theory to use for single point energy calculations. @@ -201,6 +202,7 @@ class Scheduler(object): Allowed values are He, Ne, Ar, Kr, H2, N2, O2. composite_method (str): A composite method to use. conformer_level (dict): The level of theory to use for conformer comparisons. + conf_generation_level (dict): The level of theory to use for conformer generation. opt_level (dict): The level of theory to use for geometry optimizations. freq_level (dict): The level of theory to use for frequency calculations. sp_level (dict): The level of theory to use for single point energy calculations. @@ -229,6 +231,7 @@ def __init__(self, project_directory: str, composite_method: Optional[Level] = None, conformer_level: Optional[Level] = None, + conf_generation_level: Optional[Level] = None, opt_level: Optional[Level] = None, freq_level: Optional[Level] = None, sp_level: Optional[Level] = None, @@ -308,6 +311,7 @@ def __init__(self, self.servers = list() self.composite_method = composite_method self.conformer_level = conformer_level + self.conf_generation_level = conf_generation_level self.ts_guess_level = ts_guess_level self.opt_level = opt_level self.freq_level = freq_level @@ -424,8 +428,9 @@ def __init__(self, if not self.job_types['conformers'] and len(species.conformers) == 1: # conformers weren't asked for, assign initial_xyz species.initial_xyz = species.conformers[0] - if species.label not in self.running_jobs: - self.running_jobs[species.label if not species.multi_species else species.multi_species] = list() + label = species.label if not species.multi_species else species.multi_species + if label not in self.running_jobs.keys(): + self.running_jobs[label] = list() if self.output[species.label]['convergence']: continue if species.is_monoatomic(): @@ -563,9 +568,9 @@ def schedule_jobs(self): self.determine_most_stable_conformer(label) # also checks isomorphism if self.species_dict[label].initial_xyz is not None: # if initial_xyz is None, then we're probably troubleshooting conformers, don't opt - if not self.composite_method: + if not self.composite_method and (self.job_types['opt'] or self.job_types['fine']): self.run_opt_job(label, fine=self.fine_only) - else: + elif self.composite_method and self.job_types['composite']: self.run_composite_job(label) self.timer = False break @@ -1098,7 +1103,10 @@ def run_conformer_jobs(self, labels: Optional[List[str]] = None): n_confs=n_confs, e_confs=self.e_confs, plot_path=os.path.join(self.project_directory, 'output', 'Species', - label, 'geometry', 'conformers')) + label, 'geometry', 'conformers'), + conf_generation_level=self.conf_generation_level if self.conf_generation_level is not None else None, + conf_path=os.path.join(self.project_directory, 'calcs', 'Species', f'{label}_multi'), + ) self.process_conformers(label) # TSs: elif self.species_dict[label].is_ts \ diff --git a/arc/settings/settings.py b/arc/settings/settings.py index 5d5d86d5fe..20bea287ac 100644 --- a/arc/settings/settings.py +++ b/arc/settings/settings.py @@ -175,6 +175,9 @@ } default_levels_of_theory = {'conformer': 'wb97xd/def2svp', # it's recommended to choose a method with dispersion + 'conf_generation': {'method': 'UFF', + 'solvation_method': 'SMD', + 'solvent': 'water'}, # good default for Gaussian 'ts_guesses': 'wb97xd/def2svp', 'opt': 'wb97xd/def2tzvp', # good default for Gaussian # 'opt': 'wb97m-v/def2tzvp', # good default for QChem diff --git a/arc/species/conformers.py b/arc/species/conformers.py index d427435bf1..8b787aee59 100644 --- a/arc/species/conformers.py +++ b/arc/species/conformers.py @@ -38,6 +38,7 @@ import copy import logging +import os import sys import time from itertools import product @@ -53,16 +54,21 @@ from rmgpy.molecule.molecule import Atom, Bond, Molecule from rmgpy.molecule.element import C as C_ELEMENT, H as H_ELEMENT, F as F_ELEMENT, Cl as Cl_ELEMENT, I as I_ELEMENT -from arc.common import (convert_list_index_0_to_1, +from arc.common import (ARC_PATH, + convert_list_index_0_to_1, determine_top_group_indices, get_single_bond_length, generate_resonance_structures, logger, + save_yaml_file, + read_yaml_file, ) from arc.exceptions import ConformerError, InputError +from arc.job.local import execute_command import arc.plotter +from arc.level import Level from arc.species import converter, vectors - +from arc.parser import parse_xyz_from_file # The number of conformers to generate per range of heavy atoms in the molecule # (will be increased if there are chiral centers) @@ -164,6 +170,8 @@ def generate_conformers(mol_list: Union[List[Molecule], Molecule], return_all_conformers=False, plot_path=None, print_logs=True, + conf_generation_level=None, + conf_path=None, ) -> Optional[Union[list, Tuple[list, list]]]: """ Generate conformers for (non-TS) species starting from a list of RMG Molecules. @@ -195,6 +203,8 @@ def generate_conformers(mol_list: Union[List[Molecule], Molecule], If None, the plot will not be shown (nor saved). print_logs (bool, optional): Whether define a logger so logs are also printed to stdout. Useful when run outside of ARC. True to print. + conf_generation_level (str, optional): The level of conformer generation to perform. + conf_path (str, optional): A folder path in which the conformers will be saved. Raises: ConformerError: If something goes wrong. @@ -267,7 +277,7 @@ def generate_conformers(mol_list: Union[List[Molecule], Molecule], torsions, tops = determine_rotors(mol_list) conformers = generate_force_field_conformers( mol_list=mol_list, label=label, xyzs=xyzs, torsion_num=len(torsions), charge=charge, multiplicity=multiplicity, - num_confs=num_confs_to_generate, force_field=force_field) + num_confs=num_confs_to_generate, force_field=force_field, conf_generation_level=conf_generation_level, conf_path=conf_path) lowest_confs = list() if len(conformers): @@ -630,6 +640,8 @@ def generate_force_field_conformers(label, xyzs=None, num_confs=None, force_field='MMFF94s', + conf_generation_level=None, + conf_path=None, ) -> List[dict]: """ Generate conformers using RDKit and OpenBabel and optimize them using a force field @@ -644,6 +656,8 @@ def generate_force_field_conformers(label, multiplicity (int): The species spin multiplicity. num_confs (int, optional): The number of conformers to generate. force_field (str, optional): The type of force field to use. + conf_generation_level (str, optional): The level of conformer generation to perform. + conf_path (str, optional): A folder path in which the conformers will be saved. Raises: ConformerError: If xyzs is given and it is not a list, or its entries are not strings. @@ -673,6 +687,12 @@ def generate_force_field_conformers(label, except ValueError as e: logger.warning(f'Could not generate conformers for {label}, failed with: {e}') if ff_xyzs: + if conf_generation_level is not None and conf_generation_level.solvent is not None: + ff_xyzs, ff_energies = get_force_field_energies_solvation(label, + ff_xyzs, + conf_generation_level=conf_generation_level, + conf_path=conf_path, + ) for xyz, energy in zip(ff_xyzs, ff_energies): conformers.append({'xyz': xyz, 'index': len(conformers), @@ -1153,6 +1173,75 @@ def get_force_field_energies(label: str, return xyzs, energies +def get_force_field_energies_solvation(label: str, + ff_xyzs: list, + conf_generation_level: Level, + conf_path: str, + ) -> Tuple[list, list]: + """ + Determine force field energies with solvation effect using Gaussian. + + Args: + label (str): The species' label. + ff_xyzs (list): Entries are xyz coordinates in dict format. + conf_generation_level (dict): The level of solvation to use. + conf_path (str): A folder path in which the conformers will be saved. + + Returns: + Tuple[list, list]: + - Entries are xyz coordinates, each in a dict format. + - Entries are the FF energies (in kJ/mol). + """ + ARC_child_path = conf_path + content = dict() + content['project'] = f'conf_gen_multi_{label}' + content['opt_level'] = {'method': conf_generation_level.method, + 'solvation_method': conf_generation_level.solvation_method or 'SMD', + 'solvent': conf_generation_level.solvent, + } + content['report_e_elect'] = True + content['job_types'] = {'opt': True, + 'sp': False, + 'fine': False, + 'freq': False, + 'rotors': False, + 'conformers': False, + } + content['n_confs'] = 1 + content['compute_thermo'] = False + content['species'] = list() + species_per_job = 500 + for i, xyz in enumerate(ff_xyzs): + spc_dict = dict() + spc_dict['label'] = f'{label}_multi_{i}' + spc_dict['xyz'] = xyz + spc_dict['multi_species'] = f'{label}_multi_cluster_{i//species_per_job}' + content['species'].append(spc_dict) + save_yaml_file(path=os.path.join(ARC_child_path, 'input.yml'), content=content) + commands = ['source ~/.bashrc', 'conda activate arc_env', f'python {ARC_PATH}/ARC.py {ARC_child_path}/input.yml'] + command = '; '.join(commands) + stdout, stderr = execute_command(command=command, shell=True, executable='/bin/bash') + with open(os.path.join(ARC_child_path, 'stdout.txt'), 'w') as f: + stdout_str = '\n'.join(str(x) for x in stdout) + f.write(stdout_str) + with open(os.path.join(ARC_child_path, 'stderr.txt'), 'w') as f: + stderr_str = '\n'.join(str(x) for x in stderr) + f.write(stderr_str) + if len(stderr) > 0 or len(stdout) == 0: + logger.warning(f'Got the following error when trying to submit job:\n{stderr}.') + xyzs = list() + energies = list() + content = dict() + energy_dict = read_yaml_file(os.path.join(ARC_child_path, 'output', 'e_elect_summary.yml')) + for i in range(len(ff_xyzs)): + xyzs.append(parse_xyz_from_file(os.path.join(ARC_child_path, 'output', 'Species', f'{label}_multi_{i}', 'geometry', f'{label}_multi_{i}.xyz'))) + energies.append(energy_dict[f'{label}_multi_{i}']) + content[f'{label}_multi_{i}'] = {'xyz': xyzs[-1], 'energy': energies[-1]} + save_yaml_file(path=os.path.join(ARC_child_path, 'output', 'energy_geometry_summary.yml'), content=content) + logger.info(f'{label} conformer with solvation effect are spawned from a subprocess.') + return xyzs, energies + + def openbabel_force_field_on_rdkit_conformers(label, rd_mol, force_field='MMFF94s', optimize=True): """ Optimize RDKit conformers by OpenBabel using a force field (MMFF94 or MMFF94s are recommended). diff --git a/arc/species/species.py b/arc/species/species.py index 2ae0968a10..727aab738f 100644 --- a/arc/species/species.py +++ b/arc/species/species.py @@ -1059,6 +1059,8 @@ def generate_conformers(self, n_confs: int = 10, e_confs: float = 5, plot_path: str = None, + conf_generation_level: dict = None, + conf_path: str = None, ) -> None: """ Generate conformers. @@ -1070,6 +1072,8 @@ def generate_conformers(self, (unique) generated conformers will be stored in the .conformers attribute. plot_path (str, optional): A folder path in which the plot will be saved. If None, the plot will not be shown (nor saved). + conf_generation_level (dict, optional): The level of theory to use for conformer generation. + conf_path (str, optional): A folder path in which the conformers will be saved. """ if self.is_ts: return @@ -1095,6 +1099,8 @@ def generate_conformers(self, return_all_conformers=False, plot_path=plot_path, diastereomers=diastereomers, + conf_generation_level=conf_generation_level, + conf_path=conf_path, ) if len(lowest_confs): self.conformers.extend([conf['xyz'] for conf in lowest_confs])