diff --git a/docs/api/examol.store.rst b/docs/api/examol.store.rst index 74d343da5..45488e745 100644 --- a/docs/api/examol.store.rst +++ b/docs/api/examol.store.rst @@ -41,3 +41,33 @@ examol.store.recipes .. automodule:: examol.store.recipes :members: :show-inheritance: + +examol.store.recipes.base ++++++++++++++++++++++++++ + +.. automodule:: examol.store.recipes.base + :members: + :show-inheritance: + + +examol.store.recipes.basic +++++++++++++++++++++++++++ + +.. automodule:: examol.store.recipes.basic + :members: + :show-inheritance: + + +examol.store.recipes.redox +++++++++++++++++++++++++++ + +.. automodule:: examol.store.recipes.redox + :members: + :show-inheritance: + +examol.store.recipes.surface +++++++++++++++++++++++++++++ + +.. automodule:: examol.store.recipes.surface + :members: + :show-inheritance: diff --git a/docs/components/store.rst b/docs/components/store.rst index c19ba7d31..3fb64a8ce 100644 --- a/docs/components/store.rst +++ b/docs/components/store.rst @@ -78,7 +78,7 @@ Recipes ------- Recipes define how to compute property of a molecule from multiple energy computations. -All are based on the :class:`~examol.store.recipes.PropertyRecipe` object, and provide a +All are based on the :class:`~examol.store.recipes.base.PropertyRecipe` object, and provide a function to compute the property from a molecule data record and second to generate the list of computations required to complete a computation. diff --git a/docs/conf.py b/docs/conf.py index 8a26f1626..6be4b9322 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -22,6 +22,7 @@ intersphinx_mapping = {'python': ('https://docs.python.org/3', None), 'mongoengine': ('https://docs.mongoengine.org/', None), + 'ase': ('https://wiki.fysik.dtu.dk/ase/', None), 'parsl': ('https://parsl.readthedocs.io/en/stable/', None), 'colmena': ('https://colmena.readthedocs.io/en/stable/', None), 'proxystore': ('https://docs.proxystore.dev/main/', None)} diff --git a/docs/quickstart.rst b/docs/quickstart.rst index 243a592d8..f5bc0f405 100644 --- a/docs/quickstart.rst +++ b/docs/quickstart.rst @@ -105,7 +105,7 @@ Both recipes and simulator are designed to ensure all calculations in a set are ExaMol defines a set quantum chemistry methods, which are accessible via the Simulator and enumerated in `the Simulate documentation `_. -Recipes are based on the :class:`~examol.store.recipes.PropertyRecipe` class +Recipes are based on the :class:`~examol.store.recipes.base.PropertyRecipe` class, and implement methods to compute a certain property and determine which computations are needed. Your specification will contain the details of what you wish to compute (e.g., which solvent for a solvation energy) and the level of accuracy to compute it (e.g., which XC functional)? diff --git a/examol/score/base.py b/examol/score/base.py index a2744dbfe..6c24fc1d7 100644 --- a/examol/score/base.py +++ b/examol/score/base.py @@ -5,7 +5,7 @@ import numpy as np from examol.store.models import MoleculeRecord -from examol.store.recipes import PropertyRecipe +from examol.store.recipes.base import PropertyRecipe def collect_outputs(records: list[MoleculeRecord], recipes: list[PropertyRecipe]) -> np.ndarray: diff --git a/examol/score/nfp.py b/examol/score/nfp.py index 9f10546f4..797d1e6b3 100644 --- a/examol/score/nfp.py +++ b/examol/score/nfp.py @@ -12,7 +12,7 @@ from examol.store.models import MoleculeRecord from examol.utils.conversions import convert_string_to_nx -from examol.store.recipes import PropertyRecipe +from examol.store.recipes.base import PropertyRecipe from .base import Scorer, collect_outputs from .utils.tf import LRLogger, TimeLimitCallback, EpochTimeLogger diff --git a/examol/select/base.py b/examol/select/base.py index f66b891fa..c0056ddef 100644 --- a/examol/select/base.py +++ b/examol/select/base.py @@ -7,7 +7,7 @@ import numpy as np from examol.store.models import MoleculeRecord -from examol.store.recipes import PropertyRecipe +from examol.store.recipes.base import PropertyRecipe logger = logging.getLogger(__name__) diff --git a/examol/select/bayes.py b/examol/select/bayes.py index 725d21001..df9d202b1 100644 --- a/examol/select/bayes.py +++ b/examol/select/bayes.py @@ -6,7 +6,7 @@ from examol.select.base import RankingSelector, _extract_observations from examol.store.models import MoleculeRecord -from examol.store.recipes import PropertyRecipe +from examol.store.recipes.base import PropertyRecipe class ExpectedImprovement(RankingSelector): diff --git a/examol/select/botorch.py b/examol/select/botorch.py index d8abd0b73..823a59b78 100644 --- a/examol/select/botorch.py +++ b/examol/select/botorch.py @@ -19,7 +19,7 @@ from examol.select.base import RankingSelector, _extract_observations from examol.store.models import MoleculeRecord -from examol.store.recipes import PropertyRecipe +from examol.store.recipes.base import PropertyRecipe class _EnsembleCovarianceModel(Model): diff --git a/examol/simulate/ase/__init__.py b/examol/simulate/ase/__init__.py index 5deaa8a17..862d4ea0b 100644 --- a/examol/simulate/ase/__init__.py +++ b/examol/simulate/ase/__init__.py @@ -28,6 +28,7 @@ # See methods in: https://github.com/exalearn/quantum-chemistry-on-polaris/blob/main/cp2k/ # We increase the cutoff slightly to be on the safe side _cutoff_lookup: dict[tuple[str, str], float] = { + ('PBE', 'DZVP-MOLOPT-SR-GTH'): 700., ('BLYP', 'SZV-MOLOPT-GTH'): 700., ('BLYP', 'DZVP-MOLOPT-GTH'): 700., ('B3LYP', 'def2-SVP'): 500., @@ -36,7 +37,7 @@ } # Base input file -_cp2k_inp = """&FORCE_EVAL +_cp2k_nopbc_inp = """&FORCE_EVAL &DFT &XC &XC_FUNCTIONAL $XC$ @@ -69,6 +70,36 @@ &END &END FORCE_EVAL """ +_cp2k_pbc_inp = """ +&FORCE_EVAL +&DFT + &SCF + ADDED_MOS -1 ! Add as many as possible + &SMEAR ON + METHOD FERMI_DIRAC + ELECTRONIC_TEMPERATURE [K] 300 + &END SMEAR + &MIXING + METHOD BROYDEN_MIXING + ALPHA 0.2 + BETA 1.5 + NMIXING 8 + &END MIXING + &END SCF + &XC + &XC_FUNCTIONAL $XC$ + &END XC_FUNCTIONAL + &END XC + &MGRID + NGRIDS 5 + REL_CUTOFF 60 + &END MGRID + &POISSON + PERIODIC XYZ + PSOLVER PERIODIC + &END POISSON +&END DFT +&END FORCE_EVAL""" # Solvent data (solvent name -> (gamma, e0) for CP2K, solvent name - xTB/G name for xTB/G) _solv_data = { @@ -178,6 +209,8 @@ def create_configuration(self, name: str, xyz: str, charge: int, solvent: str | xc_name = xc_name.upper() # Determine the proper basis set, pseudopotential, and method + input_template = _cp2k_nopbc_inp + stresses = True if xc_name in ['B3LYP', 'WB97X-D3']: xc_name = xc_name.replace("-", "_") # Underscores used in LibXC xc_section = f'\n&HYB_GGA_XC_{xc_name}\n&END HYB_GGA_XC_{xc_name}' @@ -199,11 +232,25 @@ def create_configuration(self, name: str, xyz: str, charge: int, solvent: str | pp_file_name = None # Use the default method = 'GPW' + elif xc_name == 'PBE': + # Used only for fully-periodic computations + input_template = _cp2k_pbc_inp + + basis_set_name = f'{basis_set_id}-MOLOPT-SR-GTH'.upper() + basis_set_file = 'BASIS_MOLOPT' + + xc_section = xc_name + + potential = 'GTH-PBE' + pp_file_name = None + + method = 'GPW' + stresses = True else: # pragma: no-coverage raise ValueError(f'XC functional "{xc_name}" not yet supported') # Inject the proper XC functional - inp = _cp2k_inp + inp = input_template inp = inp.replace('$XC$', xc_section) inp = inp.replace('$METHOD$', method) @@ -238,13 +285,13 @@ def create_configuration(self, name: str, xyz: str, charge: int, solvent: str | uks=charge != 0, inp=inp, cutoff=cutoff * units.Ry, - max_scf=32, + max_scf=256 if xc_name == 'PBE' else 32, # More steps for PBE calculations, which lack an OT outer loops basis_set_file=str(basis_set_file), basis_set=basis_set_name, pseudo_potential=potential, potential_file=pp_file_name, poisson_solver=None, - stress_tensor=False, + stress_tensor=stresses, command=self.cp2k_command) } else: # pragma: no-cover @@ -259,7 +306,7 @@ def optimize_structure(self, mol_key: str, xyz: str, config_name: str, charge: i calc_cfg = self.create_configuration(config_name, xyz, charge, solvent) # Parse the XYZ file into atoms - atoms = examol.utils.conversions.read_from_string(xyz, 'xyz') + atoms = examol.utils.conversions.read_from_string(xyz, 'extxyz') # Make the run directory based on a hash of the input configuration run_path = self._make_run_directory('opt', mol_key, xyz, charge, config_name, solvent) @@ -277,7 +324,7 @@ def optimize_structure(self, mol_key: str, xyz: str, config_name: str, charge: i # Overwrite our atoms with th last in the trajectory with Trajectory(traj_path, mode='r') as traj: atoms = traj[-1] - except InvalidULMFileError: + except (InvalidULMFileError, IndexError): traj_path.unlink() pass @@ -337,11 +384,11 @@ def optimize_structure(self, mol_key: str, xyz: str, config_name: str, charge: i # Convert to the output format out_traj = [] - out_strc = examol.utils.conversions.write_to_string(atoms, 'xyz') + out_strc = examol.utils.conversions.write_to_string(atoms, 'extxyz', columns=['symbols', 'positions', 'move_mask']) out_result = SimResult(config_name=config_name, charge=charge, solvent=solvent, xyz=out_strc, energy=atoms.get_potential_energy(), forces=atoms.get_forces()) for atoms in traj_lst: - traj_xyz = examol.utils.conversions.write_to_string(atoms, 'xyz') + traj_xyz = examol.utils.conversions.write_to_string(atoms, 'extxyz', columns=['symbols', 'positions', 'move_mask']) traj_res = SimResult(config_name=config_name, charge=charge, solvent=solvent, xyz=traj_xyz, energy=atoms.get_potential_energy(), forces=atoms.get_forces()) out_traj.append(traj_res) @@ -373,7 +420,10 @@ def _prepare_atoms(self, atoms: ase.Atoms, charge: int, config: dict): config: Configuration detail """ if 'cp2k' in config['name']: - add_vacuum_buffer(atoms, buffer_size=config['buffer_size'], cubic=re.match(r'PSOLVER\s+MT', config['kwargs']['inp'].upper()) is None) + if any(atoms.pbc): + atoms.pbc = True # All or none, never partial PBC + else: + add_vacuum_buffer(atoms, buffer_size=config['buffer_size'], cubic=re.match(r'PSOLVER\s+MT', config['kwargs']['inp'].upper()) is None) elif 'xtb' in config['name'] or 'mopac' in config['name']: utils.initialize_charges(atoms, charge) @@ -389,7 +439,7 @@ def compute_energy(self, mol_key: str, xyz: str, config_name: str, charge: int = run_path = self._make_run_directory('single', mol_key, xyz, charge, config_name, solvent) # Parse the XYZ file into atoms - atoms = examol.utils.conversions.read_from_string(xyz, 'xyz') + atoms = examol.utils.conversions.read_from_string(xyz, 'extxyz') # Run inside a temporary directory old_path = Path.cwd() @@ -410,7 +460,7 @@ def compute_energy(self, mol_key: str, xyz: str, config_name: str, charge: int = # Report the results if self.ase_db_path is not None: self.update_database([atoms], config_name, charge, solvent) - out_strc = examol.utils.conversions.write_to_string(atoms, 'xyz') + out_strc = examol.utils.conversions.write_to_string(atoms, 'extxyz', columns=['symbols', 'positions', 'move_mask']) out_result = SimResult(config_name=config_name, charge=charge, solvent=solvent, xyz=out_strc, energy=energy, forces=forces) succeeded = True # So tht we know whether to delete output directory diff --git a/examol/simulate/base.py b/examol/simulate/base.py index 16cc8817a..afffeeb9d 100644 --- a/examol/simulate/base.py +++ b/examol/simulate/base.py @@ -24,22 +24,26 @@ class SimResult: # Outputs xyz: str = field(repr=False) - """XYZ-format structure, adjusted such that the center of mass is at the origin""" + """XYZ-format structure. Adjusted such that the center of mass is at the origin if there are no PBCs""" energy: float | None = None """Energy of the molecule (units: eV)""" forces: np.ndarray | None = None """Forces acting on each atom (units: eV/Ang)""" def __post_init__(self): - # Ensure the XYZ is centered about zero - atoms = read_from_string(self.xyz, 'xyz') - atoms.center() - self.xyz = write_to_string(atoms, 'xyz') + # Ensure the XYZ is centered about zero if we are not using PBCs + atoms = read_from_string(self.xyz, 'extxyz') + if not any(atoms.pbc): + atoms.center() + self.xyz = write_to_string(atoms, 'xyz') + else: + atoms.calc = None + self.xyz = write_to_string(atoms, 'extxyz', columns=['symbols', 'positions', 'move_mask']) @property def atoms(self) -> ase.Atoms: """ASE Atoms object representation of the structure""" - return read_from_string(self.xyz, 'xyz') + return read_from_string(self.xyz, 'extxyz') def json(self, **kwargs) -> str: """Write the record to JSON format""" diff --git a/examol/specify/__init__.py b/examol/specify/__init__.py index 753cfa0e7..0fd4b8228 100644 --- a/examol/specify/__init__.py +++ b/examol/specify/__init__.py @@ -17,7 +17,7 @@ from examol.steer.base import MoleculeThinker from examol.store.db.base import MoleculeStore from examol.store.db.memory import InMemoryStore -from examol.store.recipes import PropertyRecipe +from examol.store.recipes.base import PropertyRecipe logger = logging.getLogger(__name__) diff --git a/examol/steer/base.py b/examol/steer/base.py index 5ebeb1ee4..b3297657f 100644 --- a/examol/steer/base.py +++ b/examol/steer/base.py @@ -17,7 +17,7 @@ from examol.solution import SolutionSpecification from examol.store.db.base import MoleculeStore from examol.store.models import MoleculeRecord -from examol.store.recipes import PropertyRecipe, SimulationRequest +from examol.store.recipes.base import PropertyRecipe, SimulationRequest class MoleculeThinker(BaseThinker): diff --git a/examol/steer/baseline.py b/examol/steer/baseline.py index 3dd164773..c654a134f 100644 --- a/examol/steer/baseline.py +++ b/examol/steer/baseline.py @@ -9,7 +9,7 @@ from examol.specify import SolutionSpecification from examol.steer.base import MoleculeThinker from examol.store.db.base import MoleculeStore -from examol.store.recipes import PropertyRecipe +from examol.store.recipes.base import PropertyRecipe class BruteForceThinker(MoleculeThinker): diff --git a/examol/steer/single.py b/examol/steer/single.py index 11a6d5247..f3806d59e 100644 --- a/examol/steer/single.py +++ b/examol/steer/single.py @@ -26,7 +26,7 @@ from examol.solution import SingleFidelityActiveLearning from ..store.db.base import MoleculeStore from ..store.models import MoleculeRecord -from ..store.recipes import PropertyRecipe +from ..store.recipes.base import PropertyRecipe class SingleStepThinker(MoleculeThinker): diff --git a/examol/store/models.py b/examol/store/models.py index c1d9b2c72..9951b9ab0 100644 --- a/examol/store/models.py +++ b/examol/store/models.py @@ -10,6 +10,7 @@ from rdkit import Chem from examol.simulate.base import SimResult +from examol.simulate.initialize import generate_inchi_and_xyz from examol.utils.chemistry import parse_from_molecule_string from examol.utils.conversions import read_from_string, write_to_string @@ -267,7 +268,7 @@ def find_lowest_conformer(self, config_name: str, charge: int, solvent: str | No - Lowest-energy conformer - Energy of the structure (eV) Raises: - NoSuchConformer: If we lack a conformer with these settings + MissingData: If we lack a conformer with these settings """ # Output results @@ -287,3 +288,30 @@ def find_lowest_conformer(self, config_name: str, charge: int, solvent: str | No raise MissingData(config_name, charge, solvent) return stable_conformer, lowest_energy + + def find_closest_xyz(self, config_name: str, charge: int) -> tuple[Conformer | None, str]: + """Find the most similar conformer to a certain request + + Prioritizes first by whether a conformer was optimized with the same configuration, + then those with the closest charge, + and then by those created most recently. + + Args: + config_name: Desired computation level + charge: Desired charge + Returns: + - Conformer, if one was matched + - The XYZ closest to the target calculation. Will be generated if no conformers available + """ + # Raise error if there are no conformers + if len(self.conformers) == 0: + return None, generate_inchi_and_xyz(self.identifier.smiles)[1] + + best_conf = None + best_score = (True, float('inf'), float('inf')) + for conf in self.conformers: + my_score = (conf.config_name != config_name, abs(conf.charge - charge), -conf.date_created.timestamp()) + if my_score < best_score: + best_score = my_score + best_conf = conf + return best_conf, best_conf.xyz diff --git a/examol/store/recipes/__init__.py b/examol/store/recipes/__init__.py new file mode 100644 index 000000000..9d7f062f2 --- /dev/null +++ b/examol/store/recipes/__init__.py @@ -0,0 +1 @@ +"""Tools for computing the properties of molecules from their record""" diff --git a/examol/store/recipes.py b/examol/store/recipes/base.py similarity index 50% rename from examol/store/recipes.py rename to examol/store/recipes/base.py index a1196dcd2..ac4b43030 100644 --- a/examol/store/recipes.py +++ b/examol/store/recipes/base.py @@ -1,9 +1,6 @@ -"""Tools for computing the properties of molecules from their record""" +"""Base classes other utility classes for recipes""" from dataclasses import dataclass, field -from ase import units -import numpy as np - from examol.simulate.initialize import generate_inchi_and_xyz from examol.store.models import MoleculeRecord @@ -158,23 +155,7 @@ def suggest_computations(self, record: MoleculeRecord) -> list[SimulationRequest if len(matching_conformers) > 0: _, conformer = min(matching_conformers) else: - # If there are no conformers, make an XYZ to start with - if len(record.conformers) == 0: - _, xyz = generate_inchi_and_xyz(record.identifier.smiles) - output.append( - SimulationRequest(xyz=xyz, config_name=geometry.config_name, charge=geometry.charge, optimize=True, solvent=None) - ) - continue - - # Preference conformers created using the same method followed by same charge, followed by creation date - best_score = (True, float('inf'), float('inf')) - best_xyz = None - for conf in record.conformers: - my_score = (conf.config_name != geometry.config_name, abs(conf.charge - geometry.charge), -conf.date_created.timestamp()) - if my_score < best_score: - best_score = my_score - best_xyz = conf.xyz - + _, best_xyz = record.find_closest_xyz(geometry.config_name, geometry.charge) output.append( SimulationRequest(xyz=best_xyz, config_name=geometry.config_name, charge=geometry.charge, optimize=True, solvent=None) ) @@ -188,125 +169,3 @@ def suggest_computations(self, record: MoleculeRecord) -> list[SimulationRequest ) return output - - -class SolvationEnergy(PropertyRecipe): - """Compute the solvation energy in kcal/mol - - Args: - config_name: Name of the configuration used to compute energy - solvent: Target solvent - """ - - def __init__(self, config_name: str, solvent: str): - super().__init__('solvation_energy', level=f'{config_name}-{solvent}') - self.solvent = solvent - self.config_name = config_name - - @classmethod - def from_name(cls, name: str, level: str) -> 'SolvationEnergy': - config_name, solvent = level.split("-") - return cls(config_name, solvent) - - @property - def recipe(self) -> dict[RequiredGeometry, list[RequiredEnergy]]: - return { - RequiredGeometry(config_name=self.config_name, charge=0): [RequiredEnergy(config_name=self.config_name, charge=0, solvent=self.solvent)] - } - - def compute_property(self, record: MoleculeRecord) -> float: - # Get the lowest-energy conformer with both solvent and solvation energy - vacuum_energy: float = np.inf - output: float | None = None - for conf in record.conformers: - vac_ind = conf.get_energy_index(self.config_name, 0, None) - sol_ind = conf.get_energy_index(self.config_name, 0, self.solvent) - if (vac_ind is not None) and (sol_ind is not None) and conf.energies[vac_ind].energy < vacuum_energy: - output = conf.energies[sol_ind].energy - conf.energies[vac_ind].energy - vacuum_energy = conf.energies[sol_ind].energy - - if output is None: - raise ValueError(f'Missing data for config="{self.config_name}", solvent={self.solvent}') - return output * units.mol / units.kcal - - -class RedoxEnergy(PropertyRecipe): - """Compute the redox energy for a molecule - - The level is named by the configuration used to compute the energy, - whether a solvent was included, and whether we are computing the vertical or adiabatic energy. - - Args: - charge: Amount the charge of the molecule should change by - energy_config: Configuration used to compute the energy - solvent: Solvent in which molecule is dissolved, if any - """ - - def __init__(self, charge: int, energy_config: str, vertical: bool = False, solvent: str | None = None): - # Make the name of the property based on the charge state - assert abs(charge) > 0 - prefix = { - 1: '', - 2: 'double_' - }[abs(charge)] - name = prefix + ('reduction_potential' if charge < 0 else 'oxidation_potential') - - # Name the level based on the energy level and solvent - level = energy_config - if solvent is not None: - level += "-" + solvent - level += ("-vertical" if vertical else "-adiabatic") - super().__init__(name, level) - - # Save the settings - self.energy_config = energy_config - self.charge = charge - self.vertical = vertical - self.solvent = solvent - - @classmethod - def from_name(cls, name: str, level: str) -> 'RedoxEnergy': - # Determine the charge state - charge = -1 if 'reduction' in name else 1 - if 'double' in name: - charge *= 2 - - # Determine the level - if level.count('-') == 1: - solvent = None - config_name, approx = level.split("-") - else: - config_name, solvent, approx = level.split("-") - return cls(charge=charge, energy_config=config_name, vertical=approx == 'vertical', solvent=solvent) - - @property - def recipe(self) -> dict[RequiredGeometry, list[RequiredEnergy]]: - if self.vertical: - return { - RequiredGeometry(config_name=self.energy_config, charge=0): [ - RequiredEnergy(config_name=self.energy_config, charge=0, solvent=self.solvent), - RequiredEnergy(config_name=self.energy_config, charge=self.charge, solvent=self.solvent) - ] - } - else: - return dict( - (RequiredGeometry(config_name=self.energy_config, charge=c), - [RequiredEnergy(config_name=self.energy_config, charge=c, solvent=self.solvent)]) - for c in [0, self.charge] - ) - - def compute_property(self, record: MoleculeRecord) -> float: - # Get the lowest energy conformer of the neutral molecule - neutral_conf, neutral_energy = record.find_lowest_conformer(self.energy_config, 0, self.solvent) - - # Get the charged energy - if self.vertical: - # Get the energy for this conformer - charged_energy = neutral_conf.get_energy(self.energy_config, self.charge, self.solvent) - else: - # Ge the lowest-energy conformer for the charged molecule - charged_conf, charged_energy = record.find_lowest_conformer(self.energy_config, self.charge, self.solvent) - if charged_conf.xyz_hash == neutral_conf.xyz_hash: - raise ValueError('We do not have a relaxed charged molecule') - - return charged_energy - neutral_energy diff --git a/examol/store/recipes/basic.py b/examol/store/recipes/basic.py new file mode 100644 index 000000000..fe57bde83 --- /dev/null +++ b/examol/store/recipes/basic.py @@ -0,0 +1,47 @@ +"""Properties computed for many types of applications""" + +import numpy as np +from ase import units + +from examol.store.models import MoleculeRecord +from examol.store.recipes.base import PropertyRecipe, RequiredGeometry, RequiredEnergy + + +class SolvationEnergy(PropertyRecipe): + """Compute the solvation energy in kcal/mol + + Args: + config_name: Name of the configuration used to compute energy + solvent: Target solvent + """ + + def __init__(self, config_name: str, solvent: str): + super().__init__('solvation_energy', level=f'{config_name}-{solvent}') + self.solvent = solvent + self.config_name = config_name + + @classmethod + def from_name(cls, name: str, level: str) -> 'SolvationEnergy': + config_name, solvent = level.split("-") + return cls(config_name, solvent) + + @property + def recipe(self) -> dict[RequiredGeometry, list[RequiredEnergy]]: + return { + RequiredGeometry(config_name=self.config_name, charge=0): [RequiredEnergy(config_name=self.config_name, charge=0, solvent=self.solvent)] + } + + def compute_property(self, record: MoleculeRecord) -> float: + # Get the lowest-energy conformer with both solvent and solvation energy + vacuum_energy: float = np.inf + output: float | None = None + for conf in record.conformers: + vac_ind = conf.get_energy_index(self.config_name, 0, None) + sol_ind = conf.get_energy_index(self.config_name, 0, self.solvent) + if (vac_ind is not None) and (sol_ind is not None) and conf.energies[vac_ind].energy < vacuum_energy: + output = conf.energies[sol_ind].energy - conf.energies[vac_ind].energy + vacuum_energy = conf.energies[sol_ind].energy + + if output is None: + raise ValueError(f'Missing data for config="{self.config_name}", solvent={self.solvent}') + return output * units.mol / units.kcal diff --git a/examol/store/recipes/redox.py b/examol/store/recipes/redox.py new file mode 100644 index 000000000..7c712739c --- /dev/null +++ b/examol/store/recipes/redox.py @@ -0,0 +1,85 @@ +"""Properties specific to electrochemical applications""" +from examol.store.models import MoleculeRecord +from examol.store.recipes.base import PropertyRecipe, RequiredGeometry, RequiredEnergy + + +class RedoxEnergy(PropertyRecipe): + """Compute the redox energy for a molecule + + The level is named by the configuration used to compute the energy, + whether a solvent was included, and whether we are computing the vertical or adiabatic energy. + + Args: + charge: Amount the charge of the molecule should change by + energy_config: Configuration used to compute the energy + solvent: Solvent in which molecule is dissolved, if any + """ + + def __init__(self, charge: int, energy_config: str, vertical: bool = False, solvent: str | None = None): + # Make the name of the property based on the charge state + assert abs(charge) > 0 + prefix = { + 1: '', + 2: 'double_' + }[abs(charge)] + name = prefix + ('reduction_potential' if charge < 0 else 'oxidation_potential') + + # Name the level based on the energy level and solvent + level = energy_config + if solvent is not None: + level += "-" + solvent + level += ("-vertical" if vertical else "-adiabatic") + super().__init__(name, level) + + # Save the settings + self.energy_config = energy_config + self.charge = charge + self.vertical = vertical + self.solvent = solvent + + @classmethod + def from_name(cls, name: str, level: str) -> 'RedoxEnergy': + # Determine the charge state + charge = -1 if 'reduction' in name else 1 + if 'double' in name: + charge *= 2 + + # Determine the level + if level.count('-') == 1: + solvent = None + config_name, approx = level.split("-") + else: + config_name, solvent, approx = level.split("-") + return cls(charge=charge, energy_config=config_name, vertical=approx == 'vertical', solvent=solvent) + + @property + def recipe(self) -> dict[RequiredGeometry, list[RequiredEnergy]]: + if self.vertical: + return { + RequiredGeometry(config_name=self.energy_config, charge=0): [ + RequiredEnergy(config_name=self.energy_config, charge=0, solvent=self.solvent), + RequiredEnergy(config_name=self.energy_config, charge=self.charge, solvent=self.solvent) + ] + } + else: + return dict( + (RequiredGeometry(config_name=self.energy_config, charge=c), + [RequiredEnergy(config_name=self.energy_config, charge=c, solvent=self.solvent)]) + for c in [0, self.charge] + ) + + def compute_property(self, record: MoleculeRecord) -> float: + # Get the lowest energy conformer of the neutral molecule + neutral_conf, neutral_energy = record.find_lowest_conformer(self.energy_config, 0, self.solvent) + + # Get the charged energy + if self.vertical: + # Get the energy for this conformer + charged_energy = neutral_conf.get_energy(self.energy_config, self.charge, self.solvent) + else: + # Ge the lowest-energy conformer for the charged molecule + charged_conf, charged_energy = record.find_lowest_conformer(self.energy_config, self.charge, self.solvent) + if charged_conf.xyz_hash == neutral_conf.xyz_hash: + raise ValueError('We do not have a relaxed charged molecule') + + return charged_energy - neutral_energy diff --git a/examol/store/recipes/surface.py b/examol/store/recipes/surface.py new file mode 100644 index 000000000..9916dd730 --- /dev/null +++ b/examol/store/recipes/surface.py @@ -0,0 +1,73 @@ +"""Properties related to surface chemistry""" +from dataclasses import dataclass, field + +import ase +from pydantic import BaseModel, Field + +from examol.store.recipes.base import PropertyRecipe, SimulationRequest + + +class SurfaceSimulationRequest(SimulationRequest): + """Request for a simulation involving both a molecule and surface + + The :attr:`xyz` will include both the molecule and the surface. + New fields describe the starting position + """ + + # Metadata describing how the surface was initialized + surface_name: str = ... + """Name of the surface, which should map to the geometry of a slab in the slab library""" + conformer_hash: str = ... + """Hash of the conformer used as the starting geometry for the molecule""" + starting_atom: int = ... + """Index of the atom which starts nearest to the surface""" + orientation: tuple[float, float, float] + """Rotation angles of the molecule around the center of mass starting from the original orientation""" + + +class SurfaceSite(BaseModel): + """Adsorption site for a certain molecule + + These inputs serve as inputs to :meth:`ase.build.add_adsorbate`. + """ + + name: str | None = None + """Optional name of the surface site (e.g., bridge)""" + height: float = 2 + """Starting height of the adsorbate above the slab in Angstrom""" + coordinate: tuple[float, float, float] = ... + """Position of the adsorbate within the unit cell""" + vector: tuple[float, float, float] = ... + """Vector normal to the surface for this site""" + + +class SurfaceRecord(BaseModel): + """Describe the surface on which molecule are adsorbed""" + + # Definition of the cell + name: str = ... + """Name of the surface (e.g., Cl-terminated 111 NaCl)""" + slab: str = Field(repr=False) + """3D geometry of the surface in extXYZ format""" + surface_sites: list[SurfaceSite] = Field(default_factory=list) + """Possible starting positions for the adsorbate""" + + # Properties + energy: float = Field(repr=False) + """Energy of the cell (units: eV)""" + + +class AdsorptionEnergy(PropertyRecipe): + """Compute the adsorption energy of a molecule on a surface + + The slab name corresponds to the name of a JSON file in a directory of "slab information." + Each file contains the relaxed geometry of the slab, and information about how to place adsorbates on it. + """ + + # TODO (wardlt): Be able to support >1 surfaces for different catalysis + # TODO (wardlt): Consider how to prioritize the different required calculations across surfaces + # TODO (wardlt): Determine how to select molecule orientation + def __init__(self, surface_name: str, config_name: str): + super().__init__("adsorption_energy", f"{surface_name}-config_name") + self.surface_name = surface_name + self.config_name = config_name diff --git a/examples/redoxmers-bebop/spec.py b/examples/redoxmers-bebop/spec.py index b579012ca..4e499c6f5 100644 --- a/examples/redoxmers-bebop/spec.py +++ b/examples/redoxmers-bebop/spec.py @@ -12,7 +12,7 @@ from examol.simulate.ase import ASESimulator from examol.start.fast import RandomStarter from examol.steer.single import SingleStepThinker -from examol.store.recipes import RedoxEnergy +from examol.store.recipes.redox import RedoxEnergy from examol.select.bayes import ExpectedImprovement from examol.specify import ExaMolSpecification diff --git a/examples/redoxmers/generate_database.py b/examples/redoxmers/generate_database.py index 0bc78a5c0..42afa87c0 100644 --- a/examples/redoxmers/generate_database.py +++ b/examples/redoxmers/generate_database.py @@ -3,7 +3,7 @@ from examol.simulate.ase import ASESimulator from examol.simulate.initialize import generate_inchi_and_xyz from examol.store.models import MoleculeRecord -from examol.store.recipes import RedoxEnergy +from examol.store.recipes.redox import RedoxEnergy from rdkit import RDLogger from foundry import Foundry from pathlib import Path diff --git a/examples/redoxmers/run/report.md b/examples/redoxmers/run/report.md index 7e567afbd..490e81fc5 100644 --- a/examples/redoxmers/run/report.md +++ b/examples/redoxmers/run/report.md @@ -1,12 +1,12 @@ # Run Report -Report time: 2023-10-18 17:38:06.910593 +Report time: 2023-10-24 06:49:06.396262 ## Task Summary Measures how many tasks have run as part of the application | Task Type | Count | Node Hours | Failures | |-------------|---------|--------------|------------| -| simulation | 2 | 0.00041 | 0 (0.0%) | +| simulation | 2 | 0.00062 | 0 (0.0%) | ## Outcomes over Time The property of the molecules over time. diff --git a/examples/redoxmers/spec.py b/examples/redoxmers/spec.py index 5ddf8abb6..5cc912bcc 100644 --- a/examples/redoxmers/spec.py +++ b/examples/redoxmers/spec.py @@ -13,7 +13,7 @@ from examol.solution import SingleFidelityActiveLearning from examol.start.fast import RandomStarter from examol.steer.single import SingleStepThinker -from examol.store.recipes import RedoxEnergy +from examol.store.recipes.redox import RedoxEnergy from examol.select.baseline import GreedySelector from examol.specify import ExaMolSpecification diff --git a/scripts/lohc/0_make-surfaces/0_converge-cp2k/0_converge-cutoffs.ipynb b/scripts/lohc/0_make-surfaces/0_converge-cp2k/0_converge-cutoffs.ipynb new file mode 100644 index 000000000..072209452 --- /dev/null +++ b/scripts/lohc/0_make-surfaces/0_converge-cp2k/0_converge-cutoffs.ipynb @@ -0,0 +1,240 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "4e3b9beb-1014-4848-837e-51651cc56fa3", + "metadata": {}, + "source": [ + "# Set up Convergence\n", + "Make a certain surface, place a adsorbate and then run a force computation." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "7a4931a3-2a12-4bfe-97bf-587bca9e23d8", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "from examol.utils.conversions import write_to_string\n", + "from ase.build import surface, molecule, add_adsorbate\n", + "from ase.calculators.cp2k import CP2K\n", + "from ase.db import connect\n", + "from ase import units" + ] + }, + { + "cell_type": "markdown", + "id": "11466f8f-aee8-4f1c-8cc0-a250b4d4ba6a", + "metadata": {}, + "source": [ + "Configuration" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "9590bdf8-82f1-4ad1-b670-4aa23fb97b62", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "xc: str = 'PBE'\n", + "basis_set: str = 'DZVP-MOLOPT-SR-GTH'\n", + "cutoff: float = 900 # In Ry\n", + "rel_cutoff: float = 50" + ] + }, + { + "cell_type": "markdown", + "id": "5bbf79bb-5b68-4adc-ba93-36ad48a0f303", + "metadata": {}, + "source": [ + "## Make the surface geometry\n", + "Have a simple, 4-layered Pt (111) slab with a fair amount of vacuum layers and a methane docked up top. For now, let's explore a fully-periodic." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "0fe8a89a-182a-42fd-a798-7524e893a040", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "s1 = surface('Pt', (1, 1, 1), 4)\n", + "s1.center(vacuum=12, axis=2)\n", + "s1 *= [2, 2, 1]\n", + "s1.pbc = True" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "0d986f1f-f672-47bc-915e-d2a31884c043", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "mol = molecule('CH4')" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "895ad5fa-66c2-4851-bf6d-77029c9a8c60", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "add_adsorbate(s1, mol, height=4)" + ] + }, + { + "cell_type": "markdown", + "id": "5f66c02e-6f41-460d-8bfb-ee53e6f7e53d", + "metadata": {}, + "source": [ + "## Make the Calculator\n", + "CP2K with fully periodic boundary conditions and smearing to account for metal ions." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "99c97cd9-6a3a-48a3-95c0-bbe961dbc19d", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "cp2k_opts = dict(\n", + " xc=None,\n", + " inp=f\"\"\"\n", + "&FORCE_EVAL\n", + "&DFT\n", + " &SCF\n", + " ADDED_MOS -1 ! Add as many as possible\n", + " &SMEAR ON \n", + " METHOD FERMI_DIRAC \n", + " ELECTRONIC_TEMPERATURE [K] 300\n", + " &END SMEAR \n", + " &MIXING\n", + " METHOD BROYDEN_MIXING\n", + " &END MIXING\n", + " &END SCF\n", + " &XC\n", + " &XC_FUNCTIONAL {xc}\n", + " &END XC_FUNCTIONAL\n", + " &END XC\n", + " &MGRID\n", + " NGRIDS 5\n", + " REL_CUTOFF {rel_cutoff}\n", + " &END MGRID\n", + " &POISSON\n", + " PERIODIC XYZ\n", + " PSOLVER PERIODIC\n", + " &END POISSON\n", + "&END DFT\n", + "&END FORCE_EVAL\"\"\",\n", + " basis_set_file='BASIS_MOLOPT',\n", + " basis_set=basis_set,\n", + " pseudo_potential=f'GTH-{xc}',\n", + " poisson_solver=None,\n", + ")\n", + "calc = CP2K(directory='conv',\n", + " command='/home/lward/Software/cp2k-2023.2/exe/local_cuda/cp2k_shell.ssmp',\n", + " stress_tensor=True,\n", + " max_scf=250,\n", + " cutoff=cutoff * units.Ry, # Converts to eV\n", + " **cp2k_opts)" + ] + }, + { + "cell_type": "markdown", + "id": "73f71e7f-0ee1-45d0-aa43-9cb6c85e5ab4", + "metadata": {}, + "source": [ + "Run to convergence" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "04b81cf2-4a31-4989-8c91-e25a68cf32d8", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "s1.calc = calc" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "bc352472-1bf7-498e-ab1b-076c49eed5b9", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "forces = s1.get_forces()" + ] + }, + { + "cell_type": "markdown", + "id": "38d7a455-d01e-46e7-ab7c-2b04ee105fbd", + "metadata": {}, + "source": [ + "Save the forces and such" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "24a6b2f3-1c1d-49c9-ba1c-62a37afa3b14", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "with connect('convergence.db') as db:\n", + " db.write(s1, xc=xc, cutoff=cutoff, rel_cutoff=rel_cutoff, basis_set=basis_set)" + ] + }, + { + "cell_type": "markdown", + "id": "8f4392b7-7c14-4c26-9ad2-113b3f0c342b", + "metadata": {}, + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/scripts/lohc/0_make-surfaces/0_converge-cp2k/1_plot-convergence-results.ipynb b/scripts/lohc/0_make-surfaces/0_converge-cp2k/1_plot-convergence-results.ipynb new file mode 100644 index 000000000..16981b31f --- /dev/null +++ b/scripts/lohc/0_make-surfaces/0_converge-cp2k/1_plot-convergence-results.ipynb @@ -0,0 +1,176 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "e669526f-aeaa-44c8-96b0-7c8858462feb", + "metadata": {}, + "source": [ + "# Check Convergence Settings\n", + "Plot the convergence wrt cutoff choices." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "9a1402e5-6002-43b2-a504-ebc1cc9fcf18", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "from matplotlib import pyplot as plt\n", + "from ase.db import connect\n", + "import pandas as pd\n", + "import numpy as np" + ] + }, + { + "cell_type": "markdown", + "id": "32fc289f-7d3c-40e2-b2e2-e4d1aa5c8535", + "metadata": {}, + "source": [ + "## Load in the Data\n", + "Group by structure and different settings" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "f7578057-e386-44a1-8d9d-519643af3c44", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "data = []\n", + "with connect('convergence.db') as db:\n", + " for row in db.select(''):\n", + " atoms = row.toatoms()\n", + " data.append({\n", + " 'composition': atoms.get_chemical_formula(),\n", + " 'energy': atoms.get_potential_energy(),\n", + " 'forces': atoms.get_forces(),\n", + " 'n_atoms': len(atoms),\n", + " **row.key_value_pairs\n", + " })\n", + "data = pd.DataFrame(data)" + ] + }, + { + "cell_type": "markdown", + "id": "2af38fe4-1235-4750-9938-6a5f803947ef", + "metadata": { + "tags": [] + }, + "source": [ + "## Plot Energy and Force Diff wrt Cutoff\n", + "Hold the relative cutoff constant and run" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "b32513b3-ffb5-441c-bbc5-021627c20add", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "rel_cutoff_fixed: float = 50" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "5e92ee96-39f7-47b2-b8dc-d939e0122d16", + "metadata": { + "tags": [] + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoEAAAC+CAYAAABK4NXpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABknElEQVR4nO3deVhUZfsH8O8wsiogCghuiIgLiorivqWWW2lmmnsulJXmmqVlv0wrNcvK8nV5SbFcSs3eNCuUFHcLnQFll32RfWfYmbl/fxw5OALKPjNwf65rLvGcZ87cPJx55p5znkVCRATGGGOMMdas6Gk6AMYYY4wx1vg4CWSMMcYYa4Y4CWSMMcYYa4Y4CWSMMcYYa4Y4CWSMMcYYa4Y4CWSMMcYYa4Y4CWSMMcYYa4Y4CWSMMcYYa4ZaaDoATVCpVEhISICpqSkkEommw2GMNRIiQm5uLtq3bw89vab5HZjbN8aar5q2cc0yCUxISECnTp00HQZjTEPi4uLQsWNHTYfRILh9Y4xVt41rlkmgqakpAKGSzMzMNBxNA1Epgdh/gbwUoKU10HkIoCfVdFSMaVROTg46deoktgFNUbNo3wBu4xirRE3buGaZBJbdIjEzM2uajWTQWcBzA5CTUL7NrD0w6XPAaZrm4mJMSzTl26RNvn0DuI1j7Cmq28Y1zU4xzVnQWeDkq+qNIwDkJArbg85qJi7GGKsP3MYxVm84CWxKVErh2zGokp0Pt3luFMoxxpiu4TaOsXrVLG8HN1kxNyt+O1ZDQM4DoZz9qEYLi1VNqVSipKRE02E0Sfr6+pBKuY9Yk8JtHGP1ipPApkSRXL1y6RHcQGoBhUKB+Ph4EFV2VYPVlUQiQceOHdGqVStNh8LqS3XbuOqWY6yZ4ySwKWnVrnrl/lwPxP0LjNsEmDfNaTK0nVKpRHx8PExMTGBlZdWkBypoAhEhNTUV8fHxcHR05CuCTUV127jqlmOsmeMksCmxGy6MkMtJROV9ZgDo6QOqEsD/FPDc1kYNj5UrKSkBEcHKygrGxsaaDqdJsrKyQnR0NEpKSjgJbCqe2sZJhP12wxs7MsZ0Eg8MaUr0pMIUCQCAx68sSYTHzIPAaxeByTuAVlblu3+eD/y1EUi930jBMqB2U5UoVYRbEek44/cAtyLSoVTx7eTK8NXVJuiJbdxDk3bwfIGMVRMngU2N0zTglR8BM1v17Wbthe1OLwIdXYFBr5XvS70PhJwD/t0H/GcQ8MNUIOgMoOQBC9rGMyARIz+/hLnu/2D1z36Y6/4PRn5+CZ4BifVy/NLSUmzduhU9e/ZE79690bNnTyxbtgx+fn6wtLSsUF4ikUChUKht++GHHyCRSHDu3Dlx2+LFi9GxY0f0798fvXr1whtvvCEOiPn4449RXFwslv3jjz/g6uoKQ0NDrF+/vsJrXrlyBYMGDRLju3XrVr387kxHVNXGQQLM+C/PE8hYDfDt4KbIaRrQ83lhhJwiWegfYze86m/HbbsBC04Dtw8C9z2BqKvCo5UNMHCx8KjQ4LLG5hmQiLeOyivcBEvKLsRbR+XYt2AAJvWp29/Jzc0NGRkZuHXrFiwsLKBSqXD69OlqTzocHx+PAwcOYOjQoRX2bdy4EW+//TYKCwsxduxY7N+/HytXrsSWLVuwfv16GBgYAAAcHR1x8OBBnDp1CoWFhWrHSEhIwKJFi/DXX3+hV69eKCwsrFCGNQNO06DsPgUh/55HQUYcnIO+gmFBMpCbpOnIGNMpdUoCS0pKkJSUhPz8fFhZWaFNmzb1FRerKz1p9UcA6+kB3Z4VHllxgOwwIP8BUCQBV3YAbR2Avq80aLhMkF9cWul2pYrw8dnAKmdHkwD4+GwQRnSzhFRPAj2JBEb6NbslFh4ejlOnTiE2NhYWFhYAAD09PcyaNQvR0dHVOsayZcvw9ddfY8OGDVWWMTIywqhRoxAaGoo333wTADB8+HDo6enhwoUL6N69OwDgf//7X4Xn7t27FwsWLECvXr3EYxkZGdXk12RNgGdAIrb8HoTEbADohFnSF7FO/39IyTZCP00Hx5gOqXESqFAocOzYMfz000/w8fFBUVGRuK9jx46YMGECli1bhkGDBtVroKyRtO4EjP8/YMwGIPgsEHBauIVc5u7PQEEm0G8uYNxaY2E2VU4fna/V8whAUk4hnD++AAAYYt8GJ94YVqNjyOVyODo6VnrbFwCysrLQv3//Kp+/b98+9O7dG0OGDHni62RmZuL8+fNYs2YNlixZggMHDuDmzZvVmsolKCgI9vb2ePbZZ5GWloZRo0bh888/h4mJyVOfy5qGyq6I/6ochd+UI1F6tQX2dU6s8xVxxpqLGiWBX3/9NT777DN06dIF06ZNw8aNG9GhQwcYGxsjIyMDAQEBuHbtGp577jkMHToU3333HRwdHRsqdtaQWhgAzjOFRxmVCri8A8iMAv7eIuwb9BrQvr/GwmSNp3Xr1vDz81PbVjb4IioqCu7u7rhx40aVz9+xYwcOHjwIiUSCl19+GYsXL65xDCUlJbh8+TL+/vtvmJqaYunSpfj444+xc+fOGh+L6R6lirDl96AKV8SVkEIJ4Yr4lt+D8JyTDaR6PDCIsaepURJ48+ZNeHt7w9nZudL9gwcPxtKlS7F//34cPHgQV65c4SSwKSElMHyl0HcwJRDwPSI8OrgCg9yA3i8B+jzdSV0EbZ1Y6XafqAws9rj91OcfXjIIg+3bQK8WI2MHDBiAsLAwpKeno23btjV67q1bt5CQkCDepk1KSoKbmxs+/fRTvP766wDK+wTWhZ2dHVxcXMTb1XPmzOEEsBk5eD0SidlV9wGVohQjFOcRdjUNPZ+Z3YiRMaabajQ6+NSpU1UmgI8yNDTE8uXL8dprrz21LNMhUn0h2XvrBrDEE3CeJcw7+OAO8NtbwiTUrE5MDFpU+hjlaAVbc6OqJsWABICtuRFGOVrBxKBFjfsDAkC3bt3w8ssvw83NDVlZWQCESZd//PFHKJVPXot13rx5SEpKQnR0NKKjozF06FAcPHhQTACfxNTUFNnZ2dWKcd68efD29ha7oXh6eqJfP+4F1hQRESJSFWpTIF0OTX3icxZI/8aX+gfQ4c52Xj+YsWrgKWJYzUkkgN0w4OXvgXXBwPiPAPPOQP/55WUyIoGQPwBl5QMdWM1I9STYPNUJQOUzQALA5qlOdb4FdujQIfTr1w9DhgxB79690bt3b9y8eRO5ubl1Ou6TvPPOOxg3bhz69++PlJQUXL58GR07dsRXX32FAwcOoGPHjjh79iwAYQDJ1KlT0b9/fzg7OyM1NRVbt2p20vOSkhLExcUhNDQUGRkZGo1F1xWWKHE5NAUfnw3EM19exvhdV+AXlynuH9/T+onP/0U5GtlkAlNFFBD8e0OHy5jOk1AdFi4tLCzEvXv3kJKSApVKpbZv2jTtnaspJycH5ubmyM7OrvbUF+wpVEpAoickiADw53uAzwHArKMwxcyAVwFTXsqpTGFhIaKiomBvb1+j0a3loyLLb4nZmhth81Qn7gz/mMrquL7e+9o8QE7X2rd0RRH+CkiCd0gKbkSkobCk/LNEXyrBp9P7YPagzgCEPoEjP7+EpOzCqtYLwYctf4Ob8iRg0xd442p5m8RYM1DT93+tp4jx9PTEq6++irS0tAr7JBLJU28fsSbm8TkIW1kBJm2BnHjA+1Nhqple04TbyXYjuGGupUl9bPGckw18ojKQklsIa1MjDLZvw53gGxEPkKubEqUKeUWlaG0izAsZnZ6PD38LEPfbmBlhbE8rjO1hjRHdLNHSsPxjquyK+FtH5ZBAfeG4sndAlylrAc8/gKR7QJgX0H1Cw/9SjOmoWl8J7NatGyZOnIiPPvoI7drp1hUeXfumrLNKCoWVR+4cBOL+Ld/eZRSw+Fzlz1Epqz/JtQ6r7ZVAVn0NdSVw1qxZ+Oijj57aP7qoqAgHDx6EgYFBo/aP1sb2LTW3CFfup8I7JAVXw1LxQl9bbJ/RF4Bwdc/th9sY1KUNxvawRi9b06cu+ffUK+LnNwG39gCdhgBLz/OXTtZsNNqVwJSUFKxbt07nEkDWiPSNgH6zhUfiPSEZvHcS6DCwvIxKBaQGA+16A0FnAc8NQE5C+X6z9sJaobwUFNMSp06dqla5sgFyzRER4V58NrxDU+AdkoK78eoDf/ziyv8v1ZPg8JLBNTp+xSvihlAqCabG+kKB4SsBn/8KXz5jbgBdRtb5d2KsKap1Ejhz5kxcvnwZDg4O9RkPa6ps+wJTdwPPbVUftRd5CTj6MtDWEUgPq/i8nETg5KsP1z3mRJAxbVVYohRHpUskEqw94YfItDxxf58OZhjXwxrP9LRGv46t6/x6Uj0JhjkIUxntuxyBzz1DMMrREkfchgCmNoDLAiA9HDBoWefXYqypqnUSuGfPHsyaNQvXrl2Ds7Mz9PX11favWrWqzsGxJsjIXP3/yUGARFp5AghAXBTNc6OwHnITvDXMdNe6desq3S6RSGBkZIRu3brhxRdfbJJLahIR7icr4B2agkshKQhJzMHtD5+FYQvhPfpCX1vcT1ZgXE9rPNPDCtZmDdft4YW+tvjifAiuhaXhfnIuurczBSbvFKa1YoxVqdZJ4PHjx3H+/HkYGxvj8uXLan04JBIJJ4GsekasAizshKt9VSIg54HQV7C66yEz1gh8fX0hl8uhVCrRo0cPEBHCwsIglUrRs2dP7N27F++88w6uX78OJycnTYdbgVJFNRpklF9cilsR6bgUkoLLoal4kFWgtv9uXDYG2wsJ77oJPRo09kd1amOCCU428AxMwqHrUdjxcl9OABmrhlongR9++CG2bt2KjRs3Qk+PpxtkdaAsqV45RXLDxqErmsngGV1QdpXPw8ND7ISdk5MDNzc3jBw5Eq+//jrmzZuHtWvX4vz52q0L3VCqO90QEYlf8j1uROOL86HiPsMWehjm0BbjelpjbA9rdGqjuTWc3UbZwzMwCb/6PsC7E3ugbStDYYciVRgk0nc20E77EnHGNKnW2VtxcTFmz57NCSCru1bVHFxU3XJNWdBZ4Js+wA8vAKfdhH+/6SNsrwelpaXYunUrevbsid69e6Nnz55YtmwZ/Pz8YGlpWaG8RCKBQqFQ2/bDDz9AIpHg3LnyEeCLFy9Gx44d0b9/f/Tq1QtvvPEGSkqE5P/jjz9GcXGx2jFOnz4NZ2dn9O7dG05OToiOjlbbn5qainbt2mHmzJnQpC+++AKffPKJ2ig8MzMzcT1jExMTfPTRR5DJZBqMsiLPgES8dVReYQm2pOxCvHVUjm+87uPTc0EYv+syzt4tH6g1toc1OrQ2xsKhdji02BV+H03A4SWD8eqwLhpNAAHA1c4CfTuao7hUhWP/xpbv8NwI3PgGuP6VxmJjTFvVOoNbtGgRTpw4UZ+xsObKbrgwCvhJi6KZtRfKNWdBZ4Xb5o+OngbKB8/UQyLo5uaG27dv49atWwgMDERQUBCee+65ak81Eh8fjwMHDmDo0KEV9m3cuBF+fn7w9fXFvXv3sH//fgDAli1b1JJAX19ffPjhhzh//jwCAwPxzz//wNpafaWI5cuXY8qUKXX4TetHdnY2UlJSKmxPTU1FTk4OAKB169YVklxNUqoIW34PqnSyZXr4+OZiGL6/HoWI1Dy1pdp62Zri+oax+GR6H4zr2Q7GBtpzBVoikcBtpD0A4MdbMSgqfTgAbcTDrkkBp4H0CA1Fx5h2qvXtYKVSiZ07d+L8+fPo27dvhYEhX33F37pYNelJhWlgTr4KVDoFLAHFeUByoDDKuCkrzqt8u0oJ/PUeUOVHt0SYXqfrM0J9SvQAfeMavXR4eDhOnTqF2NhYWFhYAAD09PQwa9asClfiqrJs2TJ8/fXX2LBhQ5VljIyMMGrUKISGhuLNN98EICwHp6enhwsXLmDXrl1455130L59ewCokIAeO3YM7dq1g6urq9rVRk148cUXsXTpUuzatQuDBg2CRCKBj48P1q9fj+nTpwMAfHx80L17d43G+SifqIwKVwArM7q7JeYM6oyRjuVXgJ82f5+mTXG2xfY/Q2BiKMWDzAJ0tWoF2PYDHCcAYReEK4LTvtN0mIxpjVongf7+/nBxcQEABAQEqO3T9oaCaSGnacI0MJXNE2hoJswleOQlYMlfgJX2fKDWu23ta/lEEuptRyfhv3YjgSV/1OgIcrkcjo6Old72BYCsrCz079+/yufv27cPvXv3xpAhQ574OpmZmTh//jzWrFmDJUuW4MCBA7h58yZatWoFAAgKCkLXrl0xZswY5OTk4IUXXsDHH38MqVSKhIQEfPXVV7hy5Qp++eWXGv1+DeHAgQNYu3Yt5syZg9LSUhAR9PX1sWjRIvGLcM+ePfH9999rONJyKblPTwAB4OUBHTHFWbeWItSX6uHUm8PQobUx9B4d4DJqvZAE+v0EjNkAmHfUXJCMaZFaJ4He3t71GQdjQiLY8/mKgx6KFcAP04BEP+DHF4GlfwEWXTQdbbPTunVr+Pn5qW0r+8IXFRUFd3d33Lhxo8rn79ixAwcPHoREIsHLL7+MxYsXV1qupKQEMpkMnp6eICJMmzYNBw4cwPLly/H6669j586dYsKoaa1atYK7uzu+/vprREZGgojg4OCgFt+TEmdNsDat3lQt1S2nbSrtm9h5iLBSUfQ14OZ3wOTPGz8wxrRQrZNAQLgycPDgQQQHB0MikcDJyQlLly6Fubn5059cB1evXsUXX3wBmUyGxMRE/O9//xNvvTAdpyetOA2MkTmw4Ffg8BQgNURIBJd4Ama6dZWiWj5IqHx7zE3gWDUGQcz/RUicJTXv7jtgwACEhYUhPT0dbdu2rdFzb926hYSEBPTq1QsAkJSUBDc3N3z66ad4/fXXAQh9At9+++2nHsvOzg4zZsyAsbFwO3vGjBnw8fHB8uXLcevWLbi5uQEAFAoFCgoKMHHiRI2OvL148SIuXryIlJQUqFQqtX2HDh3SUFRVG2zfBrbmRkjKLqy0c4EEgI25kTjVi64qLFHin8h0PNPjYX/SUeuEJFD2g3BlsJWVZgNkTAvUemDInTt34ODggK+//hoZGRlIS0vDV199BQcHB8jl8vqMsYK8vDz069cPe/bsadDXYVqkZVtg4W/CFcDMaODIdCAvXbMxNQSDlpU/HMZVY/BMB6GcQcsa9wcEhPXAX375Zbi5uSErKwuAMD3Ijz/+CKVS+cTnzps3D0lJSYiOjkZ0dDSGDh2KgwcPigngk5iamiI7u3wZsXnz5uHChQtQqVRQKpXw8vJCv379AAAZGRnia3z55ZeYPHmyRhPALVu2YMKECbh48SLS0tKQmZmp9tBGUj0JNk8Vpkp5/Gwq+//mqU5PnC9Q2+UUlmDEjktYcvg2ospWLek6FrAbAbjMR+V9axlrfmp9JXDt2rWYNm0a3N3d0aKFcJjS0lK89tprWLNmDa5evVpvQT5u8uTJmDx5coMdn2kpM1vg1bPAoUnCFcF/9wPjNmk6qsbx1MEzACbtqPN8gYcOHcKnn36KIUOGoEWLFiAijB49Gn37NtyAnHfeeQfjxo2DsbExLly4gDlz5uDOnTvo3bs3pFIpRo8eXa0riJqwf/9+HD58GAsXLtR0KDUyqY8t9i0YUGGeQJtK5gnURWZG+ujXqTUuhaTg8I0obHmxDyCRAIvOATytGWPlqJaMjIwoODi4wvbAwEAyNjau7WFrDAD973//e2KZwsJCys7OFh9xcXEEgG7evCmWCQwMpNjYWCIiKigoIJlMRjk5OURElJSURH5+fmLZkJAQio6OJiKi4uJikslklJWVRUREKSkpJJfLxbL379+nyMhIIiIqLS0lmUxGGRkZRESUlpZGMpmMVCoVERGFh4dTeHg4ERGpVCqSyWSUlpZGREQZGRkkk8motLSUiIgiIyPp/v374uvI5XJKSUkhIqKsrCySyWRUXFxMRETR0dEUEhIilvXz86OkpCQiIsrJySGZTEYFBQVERBQbG0uBgYFi2Xv37lFCQgIRESkUCpLJZJSfn09ERPHx8RQQECCWDQgIoLi4OCIiys/PJ5lMRrm5uURElJCQQHfv3hXLBgUFUUxMjPj3kclklJ2dLda3r6+vWn1HRUWV1/ffv1LmL+uIlKWUkpJCMplMrb4jIiLU6js9PZ2IiNLT00kmk5FSqSQiooiICAoLCxOfK5PJKDU1Va2+S0pKxPoODQ0Vy/r6+lJycjIREWVnZ5NMJqOioiIiIoqJiVF7b9y9e5cSExOJiCg3N1c8X4KCgig7O1usz7J6KztOaWkpKRQK8W9eVFREhb6niHb1JNpsVv7Y1YuUAb+plS0uLqa8vDy14xYWFhIRkVKpfGLZgoICsaxKpSKFQiHWQ3FxMSkUimqVLSkpIYVCIZ7fhYWF4nlGJJxPZedoWdmyv83jZfPy8sSyZfXyaNlH6zAvL4+KioqooKCAAgIC6M6dO+L+oKAgAiCea3XVpk0b8T2rLbKzs6v9O5YqVXQzPI1+842nm+FpVKpUNUKEjeN6WCrZbThHvf7vL8rKL9Z0OIw1ipq8/4mIap0EWltb0/nz5yts9/T0JGtr69oetsaqkwRu3ry5bPortYetra1YxtnZmVauXElERGFhYQSAvL29iYho586dZGFhIZYdOnQoubm5EZGQ3ACgc+fOERHRnj17yMDAQCw7fvx4mjNnDhGV/3FOnjxJREQeHh4EQPzQnDp1Kk2dOpWIhA9FAOTh4UFERCdPnlT7w86ZM4fGjx8vvo6BgQHt2bOHiIjOnTtHAMTkzc3NjYYOHSqWtbCwoJ07dxIRkbe3NwEQk6GVK1eSs7OzWLZDhw60efNmIiLy8fEhAGIyt3HjRnJwcBDLdu/endavX09EQkL4aKL9ySefkI2NjVjWxcWFli9fTkREUVFRBIC8vLyIiGjXrl1kamoqlh0xYgQtWrSIiIQkGwCdOXOGiIj2799PUqmUqERInCZMmEAzZ84kIiHJAEDHjx8nIqIjR44QADFpeemll2jKlCni6wAgd3d3IiI6ffo0ARATyPnz59OYMWPEsiYmJrR7924iEs55AGICvGzZMnJ1dRXLWlpa0rZt24iI6Nq1awSAgoKCKCgoiCIjI8nf318se/fuXYqPjxfjv337tph0xcfHC3WvLCWKvEqxf3xFyf/+QqQspYKCArp9+7b4xSUhIUEtkQ4MDBS/uBQVFdHt27fFLy5JSUlqiXRISIiYSJeUlNDt27fFLy4pKSl0+/Ztsez9+/fFc0epVNLt27fFLy5paWl0+/ZtMdkMDw9XS6Rv374tfnHJzMyk27dvi4leZGSkWiItk8nERDo7O5tu374t/h2jo6PVvrj4+fnRgwcPqKCggHx9fcnOzk78srJy5cp6TQLfe+892rp1a70cq77U9EOgqVKpVDTx6ytkt+Ec7b/8WKL+QE50ehlRYa5mgmOsgTRaErhy5Urq2LEj/fzzzxQbG0txcXH0008/UceOHWn16tW1PWyN8ZVAQbO6EiiTUWZmplDfiQkk+2IG0U/ziEpLmvyVwKrKPu3qHl8JbJgrgatWraLWrVvT6NGj6e2336a1a9eqPTSBk8ByJ3xiyW7DORq27W8qKRXOF1IqiXb3F66i3/hOswEyVs9q+v6XEFGtesgWFxfj3Xffxf79+1FaWgoA0NfXx1tvvYUdO3bA0NCwNoetMYlEUuPRwTk5OTA3N0d2dna1V0JgWirxLvD9s4CyGOg7B5i+Tyf6/BQWFiIqKgr29vYwMtLNqTi0XWV1XN/v/bFjx1a5TyKR4NKlS3V+jZri9q1cYYkSIz+/hDRFMb6b64Kp/R7Owyn7Afh9FdDKBlh9F9Dn9yBrGmr6/q/1wBADAwPs3r0b27dvR0REBIgI3bp1g4mJZtePZM2MbT9g1g/AiQXAvZ+FkbHP7xI6geuAWn4HY9XQGHXL86VqNyN9KeYPscPui2EIeJBdngT2mwtc+RzIeQDcPQ64LtVsoIxpSK2TwNjYWHTq1AkmJiZwdnausK9z5851Dq4qCoUC4eHh4v+joqLg5+eHNm3aNOjrMi3Vcwow47/A6deAOwcBw1bAs1u0OhHU19eHRCJBamoqrKyseJWdekZESE1NhUQiqbCkZWPx8/PTuomim6NFw7tgaj9bdLM2Ld/YwgAYvkpYoej6N4DLq4C0TtPmMqaTan07WCqVIjExscLC7unp6bC2tn7qvGJ1cfny5UpvwyxatAiHDx9+6vP5dkkTJTsM/L5a+Hnc/wGj12s0nKdRKBSIj4/nq4ENRCKRoGPHjmqrdzT0ez87OxvHjh3D999/j7t37zZoO1gVbt+qqTgf+MYZyE8DXjoA9Juj6YgYq7NGux1MRJVevVAoFA3ex+mZZ57hD05W0cDFQHEecP4D4VaP80ytXl6uVatWcHR0RElJiaZDaZL09fUhldZt3sTqunTpEg4dOoRff/0VdnZ2ePnll3Hw4MFGeW1WfQlZBTBooQfLVoaAgQkwbDlwcStw7SvA+RWd6E/MWH2qcRK4bt06AMK37P/7v/9T6wOoVCrx77//8i0QpjnDVgClhUB7F61OAMtIpdJGS1RY/YqPj8fhw4dx6NAh5OXl4ZVXXkFJSQlOnz4NJycnTYfHHrPvcgS+vBAKt5H2+GCKsLwhBr0G+P8CDHgVICXqsIgWYzqpxkmgr68vAOFKoL+/PwwMDMR9BgYG6NevH9av1+7bcKyJG/WO+v9Li4U+QIzVkylTpuD69et44YUX8N1332HSpEmQSqXYv3+/pkNjVejerhWUKsJPPrFYPd4RLQ1bCOuSv3VTq/sPM9aQapwElo2GW7JkCXbv3s19Tph2SwsDjs0EJn8BdJ+g6WhYE3HhwgWsWrUKb731FhwdHTUdDquGsT2s0dWyJSLT8nDqThwWj7AXdnACyJqxWl/79vDwgJmZGYKCguDp6YmzZ8+qPRjTCv8eADKjgZMLgahrmo6GNRHXrl1Dbm4uXF1dMWTIEOzZswepqamaDos9gZ6eBEtGdAEAeNyMhlL1SL9yZSlw7yRwZoVmgmNMQ2o9OjgqKgrTp0+Hv78/JBKJOFCjbLCIJkbFVRePnmtGlCXAiYXA/b8Ag1bAq2eBjgM1HRXTkPp+7+fn5+Pnn3/GoUOH4OPjA6VSia+++gpLly6Fqanp0w/QALh9q1p+cSmGbruInMJSuL/qiuec2gk7cpOEkcLKYmDxn0CXEZoNlLFaqun7v9ZXAletWgV7e3skJyfDxMQEgYGBuHr1KlxdXXH58uXaHpax+iXVB2YdBuzHAMUK4OgMIDlQ01GxJsLExARLly7F9evX4e/vj3feeQc7duyAtbU1pk2bpunw2GNMDFpg7hBhLtmD1yPLd5jaAC4LhJ+vfamByBjTjFongbdu3cLWrVthZWUFPT096OnpYeTIkdi+fTtWrVpVnzEyVjf6RsCc40DHwUBhFvDjdCAt/GnPYqxGevTogZ07dyI+Ph4//fSTpsNhVVg0rAukehL4x2cjJaewfMeI1YBECkRcAh7INRcgY42o1kmgUqkUJ2G1tLREQkICAMDOzg6hoaH1Ex1j9cWwFTD/FGDjDOSlABc+1HRETId98MEH8PHxqXSfVCrF9OnTuW+0lmrf2hgHFgzEzffHw9rskTltLboAzrOEn6/t0khsjDW2WieBffr0wb179wAAQ4YMwc6dO3Hjxg1s3boVXbt2rbcAGas3xq2BBf8T1g19iafyYLWXmJiIF154Aba2tli2bBn++OMPFBUVaTosVk3POrWDuXElywmOWgdAAoScA1KCGz0uxhpbrZPADz/8ECqVCgDw6aefIiYmBqNGjcKff/6Jb7/9tt4CZKxetbISEkDj1uXblLxiB6sZDw8PJCcn4+TJk2jdujXeeecdWFpaYsaMGTh8+DDS0tI0HSKrBiJCUvYjt4StegC9pgo/X/9aM0Ex1ohqPTq4MhkZGbCwsKh0OTltwqPnmOiffcKKAa/+BhhqZjQnazwN+d4PDg7G77//jjNnzuDOnTsYMmQIpk2bhrlz56JDhw71+lpPwu1b9SRmF+DNIzLEZRbg5sZxMNJ/uHJP4l3grw3C2uPdntVskIzVUIOPDn5SX5g2bdpofQLImCgvHbiyE3hwBzg+Bygp0HRETEf4+flV2NarVy+89957uHHjBuLj47Fo0SJcu3atxoNEtm/fjkGDBsHU1BTW1taYPn0697NuAFatDJGmKEZGXjF+831QvsO2H7DUkxNA1izUOAnkvjCsyWjZFlj4K2BgCsRcB06+Kiwxx9hTDBgwAAMHDsS+ffuQnZ1dYb+VlRXc3Nxw5syZGi+jeeXKFaxYsQL//PMPvLy8UFpaigkTJiAvL6++wmcAWkj1sHh4FwDAwetRqMebYozpjBongdwXhjUp7V2A+SeBFsZA2AXg19eF1QMYe4IbN25gwIAB2LhxI2xtbbFgwQJxSc268vT0xOLFi9G7d2/069cPHh4eiI2NhUwmq5fjs3KzB3dCSwMpwlIUuBr22GdXfgbgvU14MNZE1WpgiEQiwahRo7Bz506EhITAx8cHQ4cOhbu7Ozp06IDRo0fjyy+/xIMHD55+MMY0zW44MOcYIDUAgn4Dfl8FPBz0xFhlhg0bBnd3dyQlJWHfvn2Ij4/Hs88+CwcHB3z22WeIj4+vt9cqu9LYpk2bejsmE5gZ6WOWaycAwtVANQm+wJXPgRvfAnl8cYM1TTVOAhuyLwxjGtNtPDDzkDBZrN8xIPxvTUfEdICxsTEWLVqEy5cv4/79+5g7dy4OHDgAe3t7TJkypc7HJyKsW7cOI0eORJ8+fSotU1RUhJycHLUHq74lI7pAIgGu3k9FWHJu+Q6HccKdgtIC4J+9mguQsQZU4ySwIfvCMKZRvaYC0/cBz30CdJ+g6WiYjnFwcMDGjRuxadMmmJmZ4fz583U+5ttvv4179+498Qv19u3bYW5uLj46depU59dtTuzatsRzvYQ1hH+RP3IFVyIBRr0j/OzjDhRkNX5wjDWwGieBDdkXhjGN6zcbGPHIsofcP5BVw5UrV7Bo0SLY2Njgvffew4wZM3Djxo06HXPlypU4e/YsvL290bFjxyrLvf/++8jOzhYfcXFxdXrd5mjVeEfsXzAQ703sqb6jx/OAVS+gKAe47a6Z4BhrQDVOAhuzLwxjGlWYDfzwgnAVgLHHxMXF4ZNPPoGDgwPGjh2LiIgIfPfdd0hISIC7uzuGDh1aq+MSEd5++238+uuvuHTpEuzt7Z9Y3tDQEGZmZmoPVjN9OphjUh8bSPUem+JMT+/hKiIQ5hQt5hHarGmp9YohDd0XhjGNu3cSiL0F/Lke8OP+razcc889B3t7e+zduxczZ85EcHAwrl+/jiVLlqBly5Z1OvaKFStw9OhRHD9+HKampkhKSkJSUhIKCngey8ZQXKpCifKRgWG9ZwjrCuenA7IfNBYXYw2hRX0cpKwvTKdOnfDBBx/US18YxjRu0GtAejjw737gzHLAoCXgNE3TUTEtYGxsjNOnT+OFF16AVCqt12Pv27cPAPDMM8+obffw8MDixYvr9bWYuiO3ovHtpXC8P7knZgx4eAte2gIY/S4Q+w/QfaJmA2SsntU5Cbxy5QoOHTqE06dPQyqV4pVXXoGbm1t9xMaYZkkkwMTtQJEC8DsK/LIUmPczryTAcPbs2QY7Nk9arDk5haVIzS3CwetReMmlQ/kKWC4LhAdjTUytbgc3VF8YxrSOnh4w7VvAaTqgKgF+XgDE3NJ0VEyLXLt2DQsWLMCwYcPEuVGPHDmC69evazgyVlPzBneGkb4eAhNy4BOVoelwGGtwNU4CG7IvDGNaSU8KzHAHHCcIc4b99iagLNF0VEwLnD59GhMnToSxsTF8fX3FJTRzc3OxbRuvNKFrLFoaiLeBK0weDQBJAcCpxcC9U40bGGMNpMZJYFlfmPj4eHz++efo0aNHQ8TFmHZpYQC88iPQZyYw9wQg1dd0REwLfPrpp9i/fz/c3d2hr19+TgwfPhxyuVyDkbHaWjqiCwDAKzgZMemPjQa+7wkE/g+49iWvKsSahBongWfPnsWLL75Y752hGdN6+sbAzIOA9SNziamUmouHaVxoaChGjx5dYbuZmRmysrIaPyBWZ92sTTGmuxWIAI8b0eo7B78OGJoBqSFA6B8aiY+x+lTrKWIA7gvD6p9SRbgVkY4zfg9wKyIdSpUWd5KPvAzsHwlk8xrZGqNSAlHXAP9fhH8bOSm3tbVFeHh4he3Xr19H165dGzUWVn/cRgpzM566E4f84kcmjDcyFxJBALj6JcCDeJiOq3USqOm+MHv37oW9vT2MjIwwcOBAXLt2rcFfkzUsz4BEjPz8Eua6/4PVP/thrvs/GPn5JXgGJGo6tIpUSsDzAyAlCPjxRUCRqvGEpNkJOgt800eY0Pu0m/DvN32E7Y3kjTfewOrVq/Hvv/9CIpEgISEBx44dw/r167F8+fJGi4PVr1GOllj+jANOvTkcJgaPTaIxdDmgbwIk+gERFzUSH2P1RUK1nI/AxcUFa9euxauvvgpTU1PcvXsXXbt2hZ+fHyZNmoSkpKT6jlV04sQJLFy4EHv37sWIESNw4MABfP/99wgKCkLnzp2f+vycnByYm5sjOzubZ9fXEp4BiXjrqByPn4xl8/fvWzAAk/rYNnZYT5YVCxyaDOTEA+adhdHDuY8krGbtgUmf89yCDSHoLHDyVaCqM+aVHyut94Z472/atAlff/01CgsLAQgreKxfvx6ffPJJvRy/prh9awSe7wP/7AU6DweW/qXpaBgT1fT9X+sk0MTEBEFBQejSpYtaEhgZGQknJyexQWwIQ4YMwYABA8RJVQGgV69emD59OrZv3/7U53MjqV2UKsLIzy8hMbvyc0YCwMbcCNc3jKu4rJOmpUcA/31GWFu0gicnJFpDpQRibgKKZKBVO8BuuDAiWluplMIVv5yEKgpIhAR8jX+F36Oh3vv5+fkICgqCSqWCk5MTWrVqVW/Hrilu3+qfSkXQe7TtyUkAdvcDlMXAkr+E9wxjWqCm7/9aTxZd1hemS5cuatsbui9McXExZDIZNm7cqLZ9woQJuHnzZqXPKSoqEm9XA0IlMe3hE5VRZQIICNd6ErML4ROVgWEObRsvsOqw6AK0MASKKttJACSA50ag5/PamVgFnQU8N6gnVNpwBZNIWKe1KEdYw7nw4b9FOUAL4yckgABAQM4DIbG1H9XgoRYWFiIgIAApKSlQqVRqd0GmTdPi5J89VUJWAb7yuo+ErAIcf/2R+W/N2gPDVwFGZoCNs+YCZKyOap0ElvWFOXTokNgX5tatW1i/fj0++uij+oxRTVpaGpRKJdq1a6e2vV27dlXegt6+fTu2bNnSYDGxuknJrd5V4+qWa1QxN4G81CcUeJiQHBgtfHAYtHz4aCX8O2K10NkcEOYgy02sWMagJaDfUpi4uj5VdUs1J1HYXpcrmMoSQCItjzk9AkgOUE/mHv35ua1AWweh7M3vAK/NAFXRp3L0e9WLQZFcu9hrwNPTEwsXLkR6enqFfRKJBEol9wvVZS30JDjj9wAlSoJfXBb6d2pdvnP8/2ksLsbqS62TwPfeew/Z2dkYO3YsCgsLMXr0aLEvzNtvv12fMVZKXM7nISKqsK3M+++/j3Xr1on/z8nJQadOnRo0PlZ9cRn51SpnbWrUwJHUQnUTjeQA4fG4oSvKf75zCLhzsOpjrPID2gijFvHPPiDg10cSxpbqyeOg14CWlkLZtHCh3+KjSWULY+CvDajYpw4Qr2D++S5gYQ8UK8qvyHWfJFz9AICgM0Dgb5VfrSvJB968Adj0EcoG/gpc+rTq323oW+VJoNSwPAHUayFMyWFkLryukTlgalP1cR7Vqt3Ty9TR22+/jVdeeQUfffRRhS+mTPdZmxlhar/2+FX+AAevR+G7uS6VFyQSlplkTMfUae3gzz77DJs2bWrUvjCWlpaQSqUVrvqlpKRU2QgbGhrC0NCwQeNitde/Y+sn7i/rEzjYvk2jxFMj1U00xmwAWncWbnEWKx7+mwcYPvJ+MWsP2PYv31ecBxTnAvRwUlqDR8qmRwDxPlW/nvOs8iTQ7xhw/asa/VoAAYok4MBI9c2PJnap94XkriqF2eU/W9gDnYYISdyjSV3Zz20e6ULSbzbQa6qwX9+k4oerSilM1puTiMqT2Id9Ahuhn1ZKSgrWrVvHCWAT5jbSHr/KH+BP/0S8P7kn2rc2Vi8QdBa4tgt46YD6HKKM6YA6JYGa6AtjYGCAgQMHwsvLCy+99JK43cvLCy+++GKDvCarP3lFpTj2bwwAYNlo4crPCEdLvDbKHgevCcs0PfqxXvbxv3mqk/YNCgGERMOs/dMTkjEbnt4ncPR64fEoIqC0SEgIjS3Ktw9yA7o+UzGpLFYIV+FMHkmYTdoA1k7q5UqreWvdoBXQyro8WdN7pMlwGCcksY8ndo/+v4zzTOFRHUbm6s99nJ5U6LN48lUIZ0glZ8ykHY3SB3PmzJm4fPkyHBwcGvy1mGb0bm+OoV3b4J/IDPxwKxrvT+6lXuDeCWG6mOtfATP+q5EYGautWo8O1mRfmLIpYvbv349hw4bhv//9L9zd3REYGAg7O7unPp9HzzW+3MIS/HgrBgevRyEjrxitDFvg+oaxaG1iIJbxDEjElt+D1AaJ2JobYfNUJ+2bHuZRYt86oNKERBtHB0dcBo5U40vTonONMriiViod1NJBSACrqO/6fu/n5+dj1qxZsLKygrOzs9rScQCwatWqOr9GTXH7Vv+8gpLx+o93YGbUArfeH4+Who98GUrwFWYIkEiBlbLyLhtMe+jaDAh10GijgzXZF2b27NlIT0/H1q1bkZiYiD59+uDPP/+sVgLIGld2QQl+uBmNg9ejkF1QAgDo0tYEK8Z2U29IAUzqY4vnnGzgE5WBlNxCWJsKt4C18grgo5ymCYlepaNsq05INMp+VPWuYGrz1BdO04RR1xps3I8fP47z58/D2NgYly9fVuuXLJFINJIEsvo3vqc1urQ1QXR6Pk7L4/HqsC7lO9u7AA7jhYmjb3wDTN2tqTBZZbR1BgQtUesrgWZmZvD19dXJ2yD8Tblx/B2UjLUn/JBbJCy75GDVEivHOeKFvrZoIa3nka7aQNe+beriFcw6qu/3vo2NDVatWoWNGzdCr75Hb9cSt28N41d5PB5kFmD+UDu0aWmgvjPmJuAxGZAaAKvvCkkG07xaTiqvyxrtSiD3hWFP08PGFAUlSnRv1worxzliirOt9l/Vqws9qfbeOq2MLl7B1DLFxcWYPXu21iSArOHMGNCx6p12w4XVQ2JvAjf3AJMafulU9hQqpdC2PWkGBG2ew7WR1PpKoDb2haku/qZc/1JyC+F+NRLpimJ8Nbu/uN0/Phu925upz7bPtIuuXcGsg/p+769duxZWVlb44IMP6iG6+sHtm4aE/Q0ce1kY0b4mAGipZRPbNzehfwI/zX16OW3u91wLjXYlkPvCMABIzinE/isROP5vLIpKhalMlo/thm7WwnQmzh2fMMqTaQddu4KpRZRKJXbu3Inz58+jb9++Fb4Mf/VVTafmYdrOKygZ7lcj8fG03nBq/8iHbLfxgMtCoMcU9ZH8rHEoS4GCTKCVlfD/zNjqPa8RJpXXZrVOAj/88ENs3bpVq/rCsMbzIKsA+y9H4MSdOBQ/TP76d2qN1c86wsGqpYajY6xx+Pv7w8VFmEA4IEB9MvCqJq9nuu2M3wP4RGfg0I0ofDmrX/kOiQR4cY/mAmuO8jOA8L+B++eFfzsNAeafFPa16129YzTCpPLarNZJIPeFab6uh6VhyWEflCiFngSDulhg1XhHjOxmyR98rFnx9vbWdAiskS0daY9z9xJx1i8B703qUfVKRryKSMNIDgRC/wLCLgDxt8sn0weAJH+he4uetPpzuGrzDAiNoNYZ3KJFi3DixIn6jIVpsbKrfQAw0M4C5sb6GNa1LX56fShOvjEMoxytOAFkjDV5AzpbwKVzaxQrVTj6TyW3HEsKgCtfAP8ZAhRXb0lM9gSlRer//2sDcOkTIO5fIQG07g2MXAssPQ+s8S/vz1w2qTyA8mUHoP7/STuA2FtC8thM1fpKIPeFaR7CUxT4j3c4QpNycW7lSOjpSWBsIMWfq0dp51q+jDWw2NhYdO7cudrlHzx4gA4dOjRgRKyxuY20x9vHfXHsnxgsf8YBRvqPDKTS0wf8jgKZ0YD8B2FdbFYzmdHA/QtA2Hkg+oaQ3JX19XN6UVj/3HGC8GjdqerjPG0GBOtewPfjhf6EMw8CPSY36K+ljWqdBHJfmKbtfnIuvrsUjnP3ElA2fvxOTKa4fi8ngKy5GjRoEKZNm4bXX38dgwcPrrRMdnY2Tp48id27d+ONN97AypUrGzlK1pAm9bZBh9bGeJBVgDN+DzB70CNfCqQtgBFrgHNrgBvfAq5LgRa8dv0TKUuEK3v3zwuPtFD1/VFXypedHPy68KiuJ00qX5AprNcedUUYSTzhU2DYimZ1G7/Wt4O9vb2rfFy6dKk+Y2SNKCghB28dlWHC11fx+10hAXzOqR1+f3ukmAAy1pwFBwfD3NwckyZNQrt27fD888/j9ddfx8qVK7FgwQIMGDAA1tbWOHz4ML744guNJYCBgYHiz0FBQYiLiwMgrPkul8uRm5sLAEhOTsbdu3fFsqGhoYiJEdb3LikpgVwuR3Z2NgAgNTUVvr6+YtmwsDBERQlrfiuVSsjlcmRmZgIA0tPTIZfLUTYLWUREBCIiIgAARAS5XC4uO5qZmQm5XC4uNxoVFYWwsDDxdXx9fZGamgpASLDlcjlKSoQViGJiYhAaWp403L17F8nJwojP3NxcyOVyFBYKS1HGxcUhKChILOvv74/ExEQAQF5eHuRyOQoKCgAIV3AfrcPAwEDEx8cL9VJchLGWuVAVF+DQ9WgkJCTg3r17YtlgwwGILbUEchNQdPsI5HI5cnJyxPr28/NTq+/o6Gi1+s7KyhLrWy6Xq9V3ZGSkWn1nZGQAADIyMiCXy6FSCV13IiMjER4eLj5XLpcjLS1Nrb5LS0vF+r5//75Y1s/PDykpKQCEKUfkcjmKi4sBCFfCQ0JCxLL37t1DUlISAEChUKjVd3x8vFp9BwQEICFBuCKX/7C+8/PzAf9fkPDdZAT89rWQAEqkCGrRF/F91wLL/0Fht+chl8uhUCgAAElJSWr1HRISgthY4dZ8cXGxWn2npKTA756/MAOC80zcL2mHqBihbKm+KeROHyKzxxwAhLTf3od89zwhKQUQHh4u1rdKpaq0vsvO2cjISLVzVi6Xi+dsVlaW2jkbHR2tds76+fmJ52xZfRcVFYn1HRwcrFbfZedsWX2XnbPx8fFqZauFmqHs7GwCQNnZ2ZoORavcjcskuw3nxMdbR+9Q4AOuI9Z01Od7v6CggE6fPk1r1qyh6dOn08SJE2n+/Pn05Zdfkr+/fz1EWztlv6Otra24zdnZmVauXElERGFhYQSAvL29iYho586dZGFhIZYdOnQoubm5ERFRQkICAaBz584REdGePXvIwMBALDt+/HiaM2eO2uuePHmSiIg8PDwIAJWUlBAR0dSpU2nq1KlERFRSUkIAyMPDg4iITp48qfZ3mTNnDo0fP158HQMDA9qzZw8REZ07d44AUEJCAhERubm50dChQ8WyFhYWtHPnTiIi8vb2JgAUFhZGREQrV64kZ2dnsWyHDh1o8+bNRETk4+NDAOju3btERLRx40ZycHAQy3bv3p3Wr19PREQBAQEEgMa/d4B+842nrVu3ko2NjVjWxcWFls8YRbTZjKL+rxcBIC8vLyIi2rVrF5mamoplR4wYQYsWLSIiopSUFAJAZ86cISKi/fv3k1QqFctOmDCBZs6cSURECoWCANDx48eJiOjIkSMEgAoLC4mI6KWXXqIpU6aIzwVA7u7uRER0+vRpAkDp6elERDR//nwaM2aMWNbExIR2795NRESenp4EgOLi4oiIaNmyZeTq6iqWtbS0pG3bthER0bVr1wgAhYSEEBHR2rVrycnJSSioUpFdx/a0acE4IvfxJPPYQABIJpMR5abQpnGtyc6qFZH/L0T5GeTk5ERr164lIqKQkBACQNeuXSMiom3btpGlpaUYg6urKy1btoyIiOLi4ggAeXp6EhHR7t27ycTERCw7ZswYmj9/PhERpaenEwA6/csvRDe+I/epxgSA6PBUovwMmjJlCr300ktERFRYWEgA6MiRI0REdPz4cQJACoWCiIhmzpxJEyZMEF9HKpXS/v37iYjozJkzBIBSUlKIiGjRokU0YsQIsaypqSnt2rWLiIi8vLwIAEVFRRER0fLly8nFxUUsa2NjQ5988gkREd28eZMAUEBAABERrV+/nhwcHGrUxtUoCYyJialJcYqPj69R+cbCSWC5lJxC8WeVSkUv7rlOK4/LKTQpR4NRMdYwmsN7v+x3vHnzprgtMDCQYmNjiUhIXmUyGeXkCO/xpKQk8vPzE8uGhIRQdHQ0EREVFxeTTCajrKwsIhKSFLlcLpa9f/8+RUZGEhFRaWkpyWQyysjIICKitLQ0kslkpFKpiIgoPDycwsPDiUhoa2QyGaWlpRERUUZGBslkMiotLSUiosjISLp//774OnK5XPwAzcrKIplMRsXFxUREFB0dLSYdRER+fn6UlJREREQ5OTkkk8mooKCAiIhiY2MpMDBQLHvv3j0xmVQoFCSTySg/P5+IhM+vsg9XIiHxK0uE8vPzSSaTUW5uLhEJyXJZ8khEFBQURDFhwUSf21PhJlOS/bxDPOeSkpLI19dXrb7LPvDL6jszM1Osb5lMplbfERERavVdlsilp6eTTCYjpVJJREQRERFi8ktEJJPJKDU1Va2+yxL0yMhICg0NFcv6+vpScnIyEQnnk0wmo6KCfKLIqxRzfi8Fex0hUgp/q7t371JiYiIREeXm5qrVd1x4CAX+vpfozEqiL3uS/1st6cG6VkSbzShv/ySSyWSUl5dHREQP4uLUvjwFBgaK9V12zpbVd2Jiolp9BwcHi/lJUVERyWQysb6Tk5PV6js0NFQ8Z0tKStTO2dRbP5FsuSXRZjMir48pLCxMrG+lUllpfZedsxEREWrnrEwmE8/ZzMxMtXM2KipK7Zz19fUVz9my+i5L5mNiYigoKEgse/fuXfGcLavvsnM2Li6O/vnnnxq1cTVaMaRdu3ZNoi9Mc5hRX6ki+ERlICW3ENamRhhs30Ztybbb0Rn49mIY7sVn4/qGsTA1Egb2FJeqYNCCp/1hTVNzeO83h99RZ1z9Arj0qTCC9a0butvXLOhsFYMrPq96eUllCfCFA1CYXb5N3wTo+kz5oA5zLRwwleQPXP8amL5PJ/tyNuiKIcHBwdi2bRsmTZoEfX19uLq6on379jAyMkJmZiaCgoIQGBgIV1dXfPHFF5g8ufmNtNEGngGJ2PJ7EBKzC8VttuZG+OgFJ5ib6OPbi2H4J1Lo19BCTwKfqAyM7yVMmMkJIGOMVV92QQlO3I5FfrESa57trr5z0OtAViww7G3dTgBPvooKc+3lJArbZ3oALS2FkbwZUcCcY8J+qb4weXNqKNB9IuA4EegyEtDX8kGFNs7AzEPl/1cphYEjDuM0F1MDqtXawYWFhfjzzz9x7do1REdHo6CgAJaWlnBxccHEiRPRp0+fhoi13jTlb8qeAYl466i80qkxH6UvlWCWaye8NcYBndqYNEpsjGlaU37vl2kOv6M2uRWRjrnu/8BIXw+3No6HRUsDTYdUf1RK4Js+6lcAK5BALUFcE1A+bUtRLmDQSncTYAC48CFw8ztg5Dpg3P8BWr5ARqOsHWxkZIQZM2ZgxowZtXk6ayBKFWHL70FPTQAXDO2M5c90Q/vWxo0SF2NNjUKhQKtWrTQdBtMCQ7u2gZOtGYISc3DcJxYrxnarunDZaha6IvrGUxJAACDA0EyYhsVxAmDyyCwShqYNGl6DIwKkD28JX/8KSA8HXjoAGDSdCye1nieQaR+fqAy1W8BVed65PSeAjNWBhYUFEhMTYWlpqelQmIZJJBIsHWmP9afu4sdb0Vg2uiv0pY9dLcqOBy5+AuQmAovOaibQJykpBDIigLT7QOp94d+0+0BqyNOfCwDP7wL6vtKwMWqCRAKM/z+gbTfg7Eog+CyQHQfM/RkwtdF0dPVCu69rshpJyX16AliTcoyxyimVSnE+NgAYMWKEOM8Xa36m9rOFZStDJOcU4U//xEpKSICA00LfsphbjR6fKC9NmDRZdhhIjyjf7n8K2DccOLUYuLwNCPgFSLoHKIurd1xT24aIVnv0nysk78ZtgARfwH0ckHjv6c/TAXwlsImISFXg2D8x1SrLq30wVr/u3buHvLw8TYfBNMSwhRSvDrPDV173cfB6FKb1a6++cpZ5B6D/PGEZuWu7ALtfGj6orFgg8LeHV/XChH8LMsr3v/A10NZB+NmqB2BkDlj2ACy7A5aOwr9tHIAj04UrmJV2NJIIo4Tthjf876NpdsOB1y8Cx2cLdXlkOrD6HmCo291COAnUcYnZBdj9dxhOyeKhVD25N6AEgI25Ea/8wRhj9Wz+kM7Y4x2Oe/HZuBOTiUFdHmtnR6wGfI8A4V7AHQ+hv9yjS5jVVJECSA9Tv33bdzbQ6wVhf1Yc4PV/FZ/XurOQ4LW0Kt/WcRCwIabyARyTP384OvixASB4WHbSDt3q51gXbboCbl7AL0uAAa/qfAII1CIJXLhwIQ4cOAATk6bTMVJXlShVeHHPDaTkCsvLPNurHYY5tMWn54Rleip5u2LzVCe1+QIZY7Vz/PhxjB49Gs7OzgB4zfTmrm0rQ7zi2hE5BaWwMNGvpICDMGVK7C1hXeEyT5pvj0i4JVs2X11aOPDnO8KVvZwHFctbdClPAq16Ar1nCAmfVffyK3uVDWp40rnrNA145ccq5gncUfU8gU2VcWtgwa/qdZYRCZh3EqbF0TE1niJGKpUiMTER1tbWAIA33ngDO3bsgIWFhVimpKQE+vraWxm6PIVCYYkShi30xA+cvZfDcTkkFRsm98BAO+GbZ1XzBG6e6oRJfZp43w3GnqC+3vujR4/G3bt3kZubC319fZSWlmL27NkYNWoUBgwYgH79+sHISDPdLnS5fdN1RFT1l4Ggs8DJhZXseOSKmll79du3aWHAkGXA+I+EMjkJwFe9yp/a0urh7duHD7thQHuXev2dRCql0J9QkVy3K5hNTU6C0EfQ0lFIlo0tnv6chgynhu//GieBenp6SEpKEpNAMzMz+Pn5oWvXrgCExbHt7OzEBaS1kS42kiVKFU7cjsPui2HY+XJfjO0p1L9SRdCTVLwK8bQVQxhrjur7vR8WFgaZTAa5XA6ZTAZfX19kZWWhRYsW6Nmzp9oi941FF9u3Jq9a8+1VoddUYPZR4WciwO+4kHC07aY+HQvTjKirwE9zgWKF8DeZd7K8r6UGNMo8gY+qLIcsLq7miCL2VCoV4Q//ROy6EIro9HwAwLF/Y8QksKrETqonwTCHto0WJ2PNkaOjIxwdHTFnzhxxW1RUFGQyGTw8PDQYGdOk8BQFjtyKxvqJPYQlOWNuVi8BbOsIdBhYPjDDqgdgYV++XyIBXOY3XOCs5uxHA0vPCwNG0sOB78cDs48BXUZoOrJqaZCBIdw3pu6ICFfD0rDTMwSBCTkAgLYtDbByXDfMG2Kn4egYY5XJzs7GX3/9he+//x53797VdDhMA4gIy4/JcD9ZgU5tTPDaqK7CLdTqeGYj4DyzYQNk9c+mD/D6JeDnucADGfDji8DU3TqRsNdqnsDjx49DLpejpKQEACd9DWHTbwFYdMgHgQk5aGXYAuue644r743F4hH2vL4vY1rm0qVLWLBgAWxsbLBlyxbY29s//UmsSZJIJFg8XPj7H74ZLcza0Kpd9Z5c3XJM+5i2Axb/AThNB1QlwJnlwN2fNR3VU9U4mxg5ciQ2b94MV1dXtGrVCvn5+fjggw/wn//8B7du3UJubm5DxNnsjOluBQOpHtxG2uPKu89g1XhHtDLkGX0Y0xbx8fH49NNP4eDggGnTpoGIcPr0aSQkJGDLli2aDo9p0IwBHWBhoo/4zAJcCEwSBlGYtUf5PA2PkwBmHZrHfHtNmb4xMNMDGP0u0K6PsJSelqvxwJAyVXWIBoRvQkqlsj7jrFfa1nG6bK6/HjamWDJC+AZJREjOKYKNOU/szFh9qa/3/pQpU+Dt7Y1x48Zh3rx5mD59Olq2bCnuDwwMRN++fTXSDmpb+9ZcfXk+FHu8w+FqZ4Ff3hr+cHTwqw/3VjKB1ys/Nr/pVpqykgIhKQSEAT0FmY0ykKfRBoZoqkP0Z599hj/++AN+fn4wMDAQE09dlJVfjL2XI3D4ZjSKS1WwMNHHnEGdYWwghUQi4QSQMS3l6emJefPmYc2aNXB1ddV0OEwLLRxmhwNXI3AnJhN347LQj+fba17KEkAAuLEb+PcAMO9nwLaf5mKqRL11LivrEL1t2zZ4enrW12ErKC4uxqxZs/DWW2812Gs0tPziUvzHOxyjdnrjv1cjUVyqwuAubfD9IlcYG/C8S4xpuxs3bsDY2Bjjxo1Djx49sHXrVoSHh2s6LKZF2pkZ4YW+7QEAh25ECRudpgFrAoBF54CXDwr/rvHnBLApKy0C7p0AchOAQ5OAkD80HZGaOieBjd0hesuWLVi7dq04S7+uuRyagjFfXMYX50ORW1iKnjam8Fg8CCfeGCpO9swY027Dhg2Du7s7kpKSsGHDBly4cAE9evTA0KFD8d133yE5uZqjQVmT5jbSHhYm+uhq+cjyYnpSwH6UMArYfhRPuNzUtTAElnoCDuOAknzg5/nAjW+FW8RaoFZJIHeIrj1bc2OkKYrQqY0xvpndH3+uGoWxPa15hDVjOsjExARLly7F9evXERQUhNGjR2Pbtm149tlnNR0a0wJ9Opjjnw/GY/WzjpoOhWmSkTkw7xTg6gaAhDWdz64ESjU/p3KNk8ApU6bA0dERt27dwtatW5GcnIxjx45hypQpkEqlWpnMFBUVIScnR+0BCJ23ywQFBSEuLg4AUFhYCLlcLo50Tk5OVpvzKzQ0FDExMQCEJfLkcjmys7MBAKmpqfD19QUgDO44fuFffPbzZQCAUqlEXkIY/vNyd1xc9wxGdTaCn5+vOOF2REQEIiIixOfK5XKkp6cDADIzMyGXy8WO5lFRUQgLCxNj8vX1RWpqKgDh1vyjU/jExMQgNDRULHv37l3xSkVubi7kcrm4wktcXByCgoLEsv7+/khMTAQA5OXlQS6Xo6CgAADw4MEDtToMDAxEfHw8AKCgoAByuRwKhQIAkJiYqLZ6QnBwMGJjY8W/j1wuF/8uycnJ8PPzU6vv6Ohotfou6wuampoKuVwulg0LC0NkZKRY33K5HBkZGQCAjIwMyOVyqFQqAEBkZKTaLTy5XI60tDS1+i4tLRXr+/79+2JZPz8/pKSkABA64srlcnGS9NjYWISEhIhl7927h6SkJACAQqFQq+/4+Hi1+g4ICEBCgtBfKD8/H3K5HPn5wiThCQkJCAgIEMsGBQWJ9V12zpbVd1JSklp9h4SEiPVdXFysVt8pKSlq9X3//n1ERQm3r0pLSyGXy5GZmQkASEtLU6vv8PBwsb5VKlWl9V12zkZGRqqds3K5XDxns7Ky1M7Z6OhotXPWz89PPGfL6ruoqEis7+DgYLX6Ljtny+r70XO2ofTo0QM7d+5EfHw8fv31Vzz/vPaPDGQNz7AFX+ljAKQtgOd3CetES/QA3yPA8VeAh59HGkM1JJFIaP78+XT79u1K9wcEBJCenl6Njrl582aCMFyqysfjr+fh4UHm5uZ1Or6tra1YxtnZmVauXElERGFhYQSAvL29iYho586dZGFhIZYdOnQoubm5ERFRQkICAaBz584REdGePXvIwMCA5DEZNOfALTKy60cte42miJRcys7OJgB08uRJ8XcAQCUlJURENHXqVJo6dSoREZWUlBAA8vDwICKikydPEgDKzs4mIqI5c+bQ+PHjxZgMDAxoz549RER07tw5AkAJCQlEROTm5kZDhw4Vy1pYWNDOnTuJiMjb25sAUFhYGBERrVy5kpydncWyHTp0oM2bNxMRkY+PDwGgu3fvEhHRxo0bycHBQSzbvXt3Wr9+PREJ5wEAunnzJhERffLJJ2RjYyOWdXFxoeXLlxMRUVRUFAEgLy8vIiLatWsXmZqaimVHjBhBixYtIiKilJQUAkBnzpwhIqL9+/eTVCoVy06YMIFmzpxJREQKhYIA0PHjx4mI6MiRIwSACgsLiYjopZdeoilTpojPBUDu7u5ERHT69GkCQOnp6URENH/+fBozZoxY1sTEhHbv3k1ERJ6engSA4uLiiIho2bJl5OrqKpa1tLSkbdu2ERHRtWvXCACFhIQQEdHatWvJyclJLGtnZ0ebNm0iIiKZTEYASCaTERHRpk2byM7OTizr5OREa9euJSKikJAQAkDXrl0jIqJt27aRpaWlWNbV1ZWWLVtGRERxcXEEgDw9PYmIaPfu3WRiYiKWHTNmDM2fP5+IiNLT0wkAnT59moiI3N3d6dFmY8qUKfTSSy8REVFhYSEBoCNHjhAR0fHjxwkAKRQKIiKaOXMmTZgwQXyuVCql/fv3ExHRmTNnCAClpKQQEdGiRYtoxIgRYllTU1PatWsXERF5eXkRAIqKiiIiouXLl5OLi4tY1sbGhj755BMiIrp58yYBoICAACISzu9H30dNUVk705R/R12jVKroUnAyXQxO0nQoTBuEehJ91p7I5/t6P3RN3/81TgJv3rxJr732GpmamlL37t1py5YtYgJBVLskMDU1lYKDg5/4KCgoUHtOTZLAwsJCys7OFh9lH4JlCQoRUWBgIMXGxhIRUUFBAclkMsrJySEioqSkJPLz8xPLhoSEUHR0NBERFRcXk0wmo6ysLCIi+icwkmZsPUp2G84Jjzfdaa37ecpQFFFpaSnJZDLKyMggIqK0tDSSyWSkUqmIiCg8PJzCw8OJiEilUpFMJqO0tDQiIsrIyCCZTEalpaVERBQZGUn3798XY5LL5eIHaFZWFslkMiouLiYioujoaDHpICLy8/OjpCShMcrJySGZTCbWb2xsLAUGBopl7927JyaTCoWCZDIZ5efnExFRfHy8+OFKJPztyxKh/Px8kslklJubS0RCslyWPBIRBQUFUUxMjPj3kclk4kmblJREvr6+avVd9oFfVt+ZmZlEJCSFZUkSEdH9+/cpIiKCiEis77JELj09nWQyGSmVSiIiioiIUDt3ZTIZpaamqtV3WYIeGRlJoaGhYllfX19KTk4mIuFNJ5PJqKioiIiIYmJiKDg4WCx79+5dSkxMJCKi3NxctfqOi4tTq29/f3968OABERHl5eWRTCajvLw8IiJ68OAB+fv7i2UDAwPF+i47Z8vqOzExUa2+g4ODxfouKipSq+/k5GS1+g4NDaXIyEgiEr6MPHrOpqamqtV3WFiYWN9KpbLS+i47ZyMiItTOWZlMJp6zmZmZaudsVFSU2jnr6+srnrNl9V2WzMfExFBQUJBafZeds2X1XXbOBgUFNfkEiZNA7XPCJ5bsNpyjsV96k1Kp0nQ4TBvkJKr/v7S4Xg5b0/d/recJzM/Px88//4xDhw7h1q1bGDRoEObPn4/evXvjueeea/D5sQ4fPow1a9bUaoqYhphHK01RhJ2eIfhFFg8VAXoSYMaAjlj7XHd0aG389AMwxhpcc5hDrzn8jrpGUVSKYdsuIreoFB6LB4lrvzMGAMjPADymAMNWAAMW1ulQNX3/13p0sKY6RMfGxsLPzw+xsbFQKpXw8/ODn5+f2BdKU6QSCf7yT4KKgOec2sFzzWh8OasfJ4CMMdbMtTJsgdmDOgEADl6P0nA0TOvIPIDUYODs24DXR43aT7Be5glszA7RH330EVxcXLB582YoFAq4uLjAxcUFd+7cqffXUqoItyLSccbvAW5FpAtrQD6UX1yKX2Tx4qAOi5YG+PSlPjj91nC4v+qK7u1M6z0exljzsXfvXtjb28PIyAgDBw7EtWvXNB0Sq4NFw7tATwJcD0/DiduxlX6usGZq5DpgzAbh5xu7gZMLgeI8QKUEoq4B/r8I/6rq/w5rrW8H67LqXC71DEjElt+DkJhdKG6zNTfCpud7ITO/BN9eDENqbhEOLxmEZ3rwpX3GdIGu3Co9ceIEFi5ciL1792LEiBE4cOAAvv/+ewQFBaFz585PfK6u/I7N0Yy9NyCPzVLbZmtuhM1TnTCpj61mgmLa495J4MwKQFkMtLYTJppWJJXvN2svjC5+wuTiNX3/cxJYSSV5BiTiraNyPK1iOrcxwdYXe3MSyJiO0JUEaciQIRgwYAD27dsnbuvVqxemT5+O7du3P/G5uvI7NjeeAYl486i8wvaySdX2LRjAiSADYv8Fjr4MFOdWsvPp60w3Wp/ApkqpImz5PeiJCaCeBPh4mhP+XjeGE0DGWL0qLi6GTCbDhAkT1LZPmDABN2/e1FBUrC7KPlcqU/ZZs+X3IL41zICOroCBSRU7H54fnhvr7dYwJ4GP8YnKULsFXBkVAT3amcGgBVcfY6x+paWlQalUol27dmrb27VrJ046/qiqJsNn2uNpnysEIDG7ED5RGY0XFNNOMTcBxZOWnSQg54FQrh5wFvOYlNwnJ4A1LccYY7Xx+OpLRFTpikzbt2+Hubm5+OjUqVNjhciqiT9XWLU9MQGsRbmn4CTwMdamRvVajjHGasLS0hJSqbTCVb+UlJQKVwcB4P3330d2drb4KFv+kmkP/lxh1daq4nu8TuWegpPAxwy2bwNbcyNUtQKyBMJorsH2bRozLMZYM2FgYICBAwfCy8tLbbuXlxeGDx9eobyhoSHMzMzUHky78OcKqza74cIo4CedLWYdhHL1gJPAx0j1JNg81QlAxT9B2f83T3WCVK+qPxBjjNXNunXr8P333+PQoUMIDg7G2rVrERsbizfffFPTobFa4M8VVm16UmEaGABVni2Tdgjl6uPl6uUoTcykPrbYt2AAbMzVL83bmBvxMH7GWIObPXs2vvnmG2zduhX9+/fH1atX8eeff8LOzk7TobFa4s8VVm1O04RpYMweOyfM2j9xepjaaJbzBGZnZ6N169aIi4t74q0TpYogi85EqqIQVq2MMLCLBX9TY0yH5eTkoFOnTsjKyoK5ubmmw2kQ1W3fmGbw5wqrNpVSmDcwLwVoaQ10HvLUK4A1beNa1FesuiQ3V5iEkUfRMdY85ebmNtkkkNs3xlh127hmeSVQpVIhISEBpqamlU658KiyrFoXv1Xrauy6Gjegu7HratxAzWInIuTm5qJ9+/bQ02uavWFq0r4Bzedvr0047sanq7HXNO6atnHN8kqgnp4eOnbsWKPn6PKoO12NXVfjBnQ3dl2NG6h+7E31CmCZ2rRvQPP422sbjrvx6WrsNYm7Jm1c0/wqzBhjjDHGnoiTQMYYY4yxZoiTwKcwNDTE5s2bYWhoqOlQakxXY9fVuAHdjV1X4wZ0O3ZtoMv1p6uxc9yNT1djb+i4m+XAEMYYY4yx5o6vBDLGGGOMNUOcBDLGGGOMNUOcBDLGGGOMNUOcBALYvn07JBIJ1qxZI24jInz88cdo3749jI2N8cwzzyAwMFDteUVFRVi5ciUsLS3RsmVLTJs2DfHx8Q0a68cffwyJRKL2sLGx0fq4AeDBgwdYsGAB2rZtCxMTE/Tv3x8ymUzrY+/SpUuFOpdIJFixYoVWx11aWooPP/wQ9vb2MDY2RteuXbF161aoVCqxjLbGnpubizVr1sDOzg7GxsYYPnw4bt++rfVxayNu3xrv766LbZyutm8At3H1Ejc1cz4+PtSlSxfq27cvrV69Wty+Y8cOMjU1pdOnT5O/vz/Nnj2bbG1tKScnRyzz5ptvUocOHcjLy4vkcjmNHTuW+vXrR6WlpQ0W7+bNm6l3796UmJgoPlJSUrQ+7oyMDLKzs6PFixfTv//+S1FRUfT3339TeHi41seekpKiVt9eXl4EgLy9vbU67k8//ZTatm1L586do6ioKDp16hS1atWKvvnmG7GMtsb+yiuvkJOTE125coXCwsJo8+bNZGZmRvHx8Vodt7bh9q3x/u662sbpavtGxG1cfcTdrJPA3NxccnR0JC8vLxozZozYSKpUKrKxsaEdO3aIZQsLC8nc3Jz2799PRERZWVmkr69PP//8s1jmwYMHpKenR56eng0W8+bNm6lfv36V7tPmuDds2EAjR46scr82x/641atXk4ODA6lUKq2O+/nnn6elS5eqbZsxYwYtWLCAiLS3zvPz80kqldK5c+fUtvfr1482bdqktXFrG27fGi9uoqbTxulK+0bEbVx9xN2sbwevWLECzz//PJ599lm17VFRUUhKSsKECRPEbYaGhhgzZgxu3rwJAJDJZCgpKVEr0759e/Tp00cs01DCwsLQvn172NvbY86cOYiMjNT6uM+ePQtXV1fMmjUL1tbWcHFxgbu7u7hfm2N/VHFxMY4ePYqlS5dCIpFoddwjR47ExYsXcf/+fQDA3bt3cf36dUyZMgWA9tZ5aWkplEoljIyM1LYbGxvj+vXrWhu3tuH2rXHjbgptnC61bwC3cfURd7NNAn/++WfI5XJs3769wr6kpCQAQLt27dS2t2vXTtyXlJQEAwMDWFhYVFmmIQwZMgQ//vgjzp8/D3d3dyQlJWH48OFIT0/X6rgjIyOxb98+ODo64vz583jzzTexatUq/Pjjj2Jc2hr7o3777TdkZWVh8eLFYkxlMVQVk6bi3rBhA+bOnYuePXtCX18fLi4uWLNmDebOnavVsZuammLYsGH45JNPkJCQAKVSiaNHj+Lff/9FYmKi1satTbh9a/y/e1No43SpfQO4jauPuFvU8XfRSXFxcVi9ejUuXLhQIRN/lEQiUfs/EVXY9rjqlKmLyZMniz87Oztj2LBhcHBwwA8//IChQ4cC0M64VSoVXF1dsW3bNgCAi4sLAgMDsW/fPrz66qtiOW2M/VEHDx7E5MmT0b59e7Xt2hj3iRMncPToURw/fhy9e/eGn58f1qxZg/bt22PRokViOW2M/ciRI1i6dCk6dOgAqVSKAQMGYN68eZDL5WIZbYxbG3D7VlFj/N2bQhunS+0bwG1cZWoad7O8EiiTyZCSkoKBAweiRYsWaNGiBa5cuYJvv/0WLVq0ELPvx7PplJQUcZ+NjQ2Ki4uRmZlZZZnG0LJlSzg7OyMsLEwcRaeNcdva2sLJyUltW69evRAbGyvGBWhn7GViYmLw999/47XXXhO3aXPc7777LjZu3Ig5c+bA2dkZCxcuxNq1a8WrQ9ocu4ODA65cuQKFQoG4uDj4+PigpKQE9vb2Wh23NuD2TTNx63obp2vtG8BtXH3E3SyTwPHjx8Pf3x9+fn7iw9XVFfPnz4efnx+6du0KGxsbeHl5ic8pLi7GlStXMHz4cADAwIEDoa+vr1YmMTERAQEBYpnGUFRUhODgYNja2oonjzbGPWLECISGhqptu3//Puzs7ABAq2Mv4+HhAWtrazz//PPiNm2OOz8/H3p66m9xqVQqTp+gzbGXadmyJWxtbZGZmYnz58/jxRdf1Im4NYnbN83ErettnK61bwC3cfUSd7WHkDRxj46eIxKGZ5ubm9Ovv/5K/v7+NHfu3EqHZ3fs2JH+/vtvksvlNG7cuAYfVv7OO+/Q5cuXKTIykv755x964YUXyNTUlKKjo7U6bh8fH2rRogV99tlnFBYWRseOHSMTExM6evSoWEZbYyciUiqV1LlzZ9qwYUOFfdoa96JFi6hDhw7i9Am//vorWVpa0nvvvaf1sXt6etJff/1FkZGRdOHCBerXrx8NHjyYiouLtTpubcXtW8P/3XW5jdPF9o2I27j6iJuTwIcebyRVKhVt3ryZbGxsyNDQkEaPHk3+/v5qzykoKKC3336b2rRpQ8bGxvTCCy9QbGxsg8ZZNleQvr4+tW/fnmbMmEGBgYFaHzcR0e+//059+vQhQ0ND6tmzJ/33v/9V26/NsZ8/f54AUGhoaIV92hp3Tk4OrV69mjp37kxGRkbUtWtX2rRpExUVFWl97CdOnKCuXbuSgYEB2djY0IoVKygrK0vr49ZW3L41zt9dV9s4XWzfiLiNq4+4JURE9XRVkzHGGGOM6Yhm2SeQMcYYY6y54ySQMcYYY6wZ4iSQMcYYY6wZ4iSQMcYYY6wZ4iSQMcYYY6wZ4iSQMcYYY6wZ4iSQMcYYY6wZ4iSQMcYYY6wZ4iSQNQs3btyAs7Mz9PX1MX369Cq3VSY0NBQ2NjbIzc2tt3hSUlJgZWWFBw8e1NsxGWPNF7dxrDY4CWSNIikpCStXrkTXrl1haGiITp06YerUqbh48WK1j3H58mVIJBJkZWXV+PXXrVuH/v37IyoqCocPH65yW2U2bdqEFStWwNTUVC2Oskfbtm0xbtw43Lhxo9rxWFtbY+HChdi8eXONfxfGmPbhNk4dt3G6gZNA1uCio6MxcOBAXLp0CTt37oS/vz88PT0xduxYrFixolFiiIiIwLhx49CxY0e0bt26ym2Pi4+Px9mzZ7FkyZIK+0JDQ5GYmIjLly/DysoKzz//PFJSUqod05IlS3Ds2DFkZmbW5ldijGkJbuMqx22cDqiHdZAZe6LJkydThw4dSKFQVNiXmZlJRERRUVEEgHx9fdX2ASBvb29x/6OPRYsWERFRYWEhrVy5kqysrMjQ0JBGjBhBPj4+asd99OHh4VHptsrs2rWLXF1d1bZ5e3sTADF2IqJ79+4RADp79iwpFAoyNTWlU6dOqT3v7NmzZGJiQjk5OeK2Ll260MGDB6tZk4wxbcRtnIDbON3DVwJZg8rIyICnpydWrFiBli1bVthf1bfTx3Xq1AmnT58GUP7tdPfu3QCA9957D6dPn8YPP/wAuVyObt26YeLEicjIyECnTp2QmJgIMzMzfPPNN0hMTMSsWbMqbJs9e3alr3v16lW4uro+Mbb8/Hx4eHgAAPT19dGyZUvMmTNH3FbGw8MDM2fOFG+5AMDgwYNx7dq1atUBY0z7cBtXjts43dNC0wGwpi08PBxEhJ49e9bpOFKpFG3atAEg9DUpa1jz8vKwb98+HD58GJMnTwYAuLu7w8vLCwcPHsS7774LGxsbSCQSmJubw8bGBgDQsmXLCtsqU3abpzIdO3YEIDSQRISBAwdi/PjxAIDXXnsNw4cPR0JCAtq3b4+0tDScO3cOXl5easfo0KEDfH19a18xjDGN4jaO2zhdxlcCWYMiIgCARCJpkONHRESgpKQEI0aMELfp6+tj8ODBCA4OrvPxCwoKYGRkVOm+a9euQS6X46effoKdnR0OHz4MfX19AMK33969e+PHH38EABw5cgSdO3fG6NGj1Y5hbGyM/Pz8OsfJGNMMbuO4jdNlnASyBuXo6AiJRPLUxkpPTzgVyxpUACgpKXnq8atqgImoXhplS0vLKjs129vbo3v37pg9eza2bNmCl156CUVFReL+1157Tbxd4uHhgSVLllSIKSMjA1ZWVnWOkzGmGdzGcRunyzgJZA2qTZs2mDhxIv7zn/8gLy+vwv6yqRDKGonExERxn5+fn1pZAwMDAIBSqRS3devWDQYGBrh+/bq4raSkBHfu3EGvXr3qHL+LiwuCgoKeWm7hwoVQqVTYu3evuG3BggWIjY3Ft99+i8DAQCxatKjC8wICAuDi4lLnOBljmsFtHLdxuoyTQNbg9u7dC6VSicGDB+P06dMICwtDcHAwvv32WwwbNgyAcMtg6NCh2LFjB4KCgnD16lV8+OGHasexs7ODRCLBuXPnkJqaCoVCgZYtW+Ktt97Cu+++C09PTwQFBeH1119Hfn4+3Nzc6hz7xIkTcevWLbVGuTJ6enpYs2YNduzYId76sLCwwIwZM/Duu+9iwoQJYv+aMvn5+ZDJZJgwYUKd42SMaQ63cdzG6SzNDEpmzU1CQgKtWLGC7OzsyMDAgDp06EDTpk0jb29vsUxQUBANHTqUjI2NqX///nThwgVx+oQyW7duJRsbG5JIJOL0CQUFBbRy5UqytLSsMH1CGXNz8wpTJFS27XGlpaXUoUMH8vT0FLdVNn0CEZFCoSALCwv6/PPPxW0XL14kAHTy5MkKxz5+/Dj16NHjia/PGNMN3MZxG6eLJESPdFBgjFWwd+9enDlzBufPn6/xc48dO4bVq1cjISFBvNVTZvDgwVizZg3mzZtXX6EyxliNcRvXfPEUMYw9xbJly5CZmYnc3Fy1+a+eJD8/H1FRUdi+fTveeOONCo1jSkoKZs6ciblz5zZEyIwxVm3cxjVffCWQsQbw8ccf47PPPsPo0aNx5swZtGrVStMhMcZYveE2rmngJJAxxhhjrBni0cGMMcYYY80QJ4GMMcYYY80QJ4GMMcYYY80QJ4GMMcYYY80QJ4GMMcYYY80QJ4GMMcYYY80QJ4GMMcYYY80QJ4GMMcYYY80QJ4GMMcYYY83Q/wPDfr2gZX4xFAAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fig, axs = plt.subplots(1, 2, figsize=(6.5, 2.), sharex=True)\n", + "for (comp, xc, basis), group in data.query(f'rel_cutoff=={rel_cutoff_fixed}').groupby(['composition', 'xc', 'basis_set']):\n", + " group = group.sort_values('cutoff')\n", + " n_atoms = group.iloc[0]['n_atoms']\n", + " label = f'{comp}'\n", + " \n", + " # Plot the energy and force diffs\n", + " cutoffs = group['cutoff']\n", + " energy_diffs = group['energy'] - group['energy'].iloc[-1]\n", + " force_diffs = group['forces'].apply(lambda x: np.subtract(x, group['forces'].iloc[-1])).apply(np.abs).apply(np.max)\n", + " \n", + " axs[0].plot(cutoffs, energy_diffs * 1000, '--o', label=label)\n", + " axs[1].plot(cutoffs, force_diffs * 1000, '--o', label=label)\n", + " \n", + "# Place the convergence levels\n", + "for ax in axs:\n", + " ax.set_xlim(ax.get_xlim())\n", + "eng_conv = 1 # meV/Ang\n", + "for i in [-1, 1]:\n", + " axs[0].plot(axs[0].get_xlim(), [eng_conv * i] * 2, 'k:', lw=1)\n", + " \n", + "force_conv = 1 # meV/Ang\n", + "axs[1].plot(axs[1].get_xlim(), [force_conv] * 2, 'k:', lw=1)\n", + "\n", + "\n", + "# Make the labels\n", + "axs[0].legend(fontsize=8)\n", + "axs[0].set_ylabel('$\\Delta E$ (meV/atom)')\n", + "axs[1].set_ylabel('$\\Delta F$ (meV/Ang)')\n", + "for ax in axs:\n", + " ax.set_xlabel('Cutoff (Ry)')\n", + "\n", + "fig.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "id": "f2038a36-1614-4684-b101-56d9b9bebf4a", + "metadata": {}, + "source": [ + "Seems like a cutoff of 7 should be sufficient based on force and energy." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "37c1c1af-fed4-4c95-91a5-b6590d49b85b", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/scripts/lohc/0_make-surfaces/0_converge-cp2k/README.md b/scripts/lohc/0_make-surfaces/0_converge-cp2k/README.md new file mode 100644 index 000000000..79925348a --- /dev/null +++ b/scripts/lohc/0_make-surfaces/0_converge-cp2k/README.md @@ -0,0 +1,8 @@ +# Verify Convergence + +There are two notebooks used here. + +The first runs an example surface and allows you to change the cutoff energies. +Use it by making a different example surface and running a single-point force computation. + +The second plots the results so that you can determine appropriate cutoff energies. diff --git a/scripts/lohc/0_make-surfaces/1_relax-surfaces/0_create-bulk.ipynb b/scripts/lohc/0_make-surfaces/1_relax-surfaces/0_create-bulk.ipynb new file mode 100644 index 000000000..5fd28f3e1 --- /dev/null +++ b/scripts/lohc/0_make-surfaces/1_relax-surfaces/0_create-bulk.ipynb @@ -0,0 +1,205 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "a6a0aada-1b9d-4357-beb3-9e512e530e3d", + "metadata": {}, + "source": [ + "# Relax Bulk with CP2K\n", + "Use CP2K PBE to relax a certain structure" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "897e19e7-ff11-45dc-b2aa-0a6d51225fc4", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "from examol.simulate.ase.utils import make_ephemeral_calculator\n", + "from examol.simulate.ase import ASESimulator\n", + "from examol.utils.conversions import write_to_string\n", + "from pathlib import Path\n", + "from ase.constraints import ExpCellFilter\n", + "from ase.optimize import BFGS\n", + "from ase.build import bulk" + ] + }, + { + "cell_type": "markdown", + "id": "00202102-0409-4e7d-9fba-0136a11c03dd", + "metadata": {}, + "source": [ + "Configuration" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "5fec8734-89c1-41b8-a37e-b8da0bb8e255", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "atoms = bulk('Pd')\n", + "name = 'fcc-Pd'\n", + "out_dir = Path('relaxed-bulk')" + ] + }, + { + "cell_type": "markdown", + "id": "a86a6581-75e8-4323-8870-3d7420be0a02", + "metadata": {}, + "source": [ + "## Make the Bulk Cell\n", + "Do whatever you need to start with a new cell" + ] + }, + { + "cell_type": "markdown", + "id": "c1b73933-908e-4563-b063-8b21e1f6e292", + "metadata": {}, + "source": [ + "## Make the Simulator\n", + "The tool we'll use to make an ASE calculator" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "bae06c56-23c7-4ad4-a829-b8d5bd60706c", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "sim = ASESimulator(cp2k_command='/home/lward/Software/cp2k-2023.2/exe/local_cuda/cp2k_shell.ssmp')" + ] + }, + { + "cell_type": "markdown", + "id": "476ff1bc-f4f7-483a-b768-d2460b97970d", + "metadata": {}, + "source": [ + "## Run the optimization\n", + "We'll use ExaMol to make ASE calculators with our target settings." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "d2ad4510-24cd-40bf-ba24-dba1f46cf43c", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "xyz = write_to_string(atoms, 'extxyz')" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "a0773854-b367-4ed8-9b11-7261535ee991", + "metadata": { + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Initial volume: 14.715967249999999. Stress: 0.05\n", + " Step Time Energy fmax\n", + "BFGS: 0 08:52:40 -3469.093567 0.781361\n", + "BFGS: 1 08:53:02 -3469.099891 0.423853\n", + "BFGS: 2 08:53:21 -3469.102291 0.013849\n", + "BFGS: 3 08:53:36 -3469.102293 0.000236\n", + "Final volume: 14.394562659381776. Stress: 0.00\n" + ] + } + ], + "source": [ + "with make_ephemeral_calculator(sim.create_configuration('cp2k_pbe_dzvp', xyz, 0, None)) as calc:\n", + "\n", + " # Start by computing the stresses acting on the cell\n", + " calc.directory = 'run'\n", + " atoms.calc = calc\n", + " stresses = atoms.get_stress()\n", + " print(f'Initial volume: {atoms.get_volume()}. Stress: {stresses[:3].sum() / 3:.2f}')\n", + " \n", + " # Run the optimization\n", + " opt_atoms = ExpCellFilter(atoms) # Allow ASE to optimize the lattice parameter\n", + " opt = BFGS(opt_atoms)\n", + " opt.run(fmax=0.001)\n", + " stresses = atoms.get_stress()\n", + " print(f'Final volume: {atoms.get_volume()}. Stress: {stresses[:3].sum() / 3:.2f}')" + ] + }, + { + "cell_type": "markdown", + "id": "2ddeff5a-2d11-4fc0-8e23-e5d9903063d1", + "metadata": {}, + "source": [ + "Save the relaxed structure" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "a0532f5b-0dcb-4f70-8f82-9d552cbfa284", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "out_dir.mkdir(exist_ok=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "76b8485d-fcdb-426d-b9f8-d975f6107ccb", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "with out_dir.joinpath(f'{name}.extxyz').open('w') as fp:\n", + " print(write_to_string(atoms, 'extxyz'), file=fp)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a64c5385-22fc-4007-b504-303a15acd7e2", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/scripts/lohc/0_make-surfaces/1_relax-surfaces/1_generate-surfaces.ipynb b/scripts/lohc/0_make-surfaces/1_relax-surfaces/1_generate-surfaces.ipynb new file mode 100644 index 000000000..e77158966 --- /dev/null +++ b/scripts/lohc/0_make-surfaces/1_relax-surfaces/1_generate-surfaces.ipynb @@ -0,0 +1,240 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "fba20aac-3b71-49c7-a1ec-9d441d25e572", + "metadata": {}, + "source": [ + "# Generate Unrelaxed Surfaces\n", + "Use [CatKit](https://catkit.readthedocs.io/en/latest/) to produce initial slab structures and enumerate the adsorbate sites." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "07e0242f-77a4-475b-b560-8c94f511afb9", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "from catkit.gen.adsorption import AdsorptionSites\n", + "from catkit.gen.surface import SlabGenerator\n", + "from pathlib import Path\n", + "from ase.io import read\n", + "from examol.utils.conversions import read_from_string, write_to_string\n", + "import json" + ] + }, + { + "cell_type": "markdown", + "id": "0313ffbe-612c-4902-ad75-7020ec6feed3", + "metadata": {}, + "source": [ + "Configuration" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "63abefda-dd14-4575-80ec-702831ff7543", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "bulk_path: str = './relaxed-bulk/fcc-Pd.extxyz' # Path to a relaxed bulk structure\n", + "facet: tuple[int, int, int] = (1, 1, 1) # Which surface to use\n", + "supercell: tuple[int, int] = (3, 3) # Size of the supercell\n", + "layers: int = 6 # Number of layers for the slab" + ] + }, + { + "cell_type": "markdown", + "id": "a2a1afcf-88ab-4a62-8203-ce310d374070", + "metadata": {}, + "source": [ + "Derived" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "57b83281-d189-466a-84f7-5d5ec290e093", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "bulk_path = Path(bulk_path)\n", + "bulk_name = bulk_path.name.split(\".\")[0]\n", + "output_name = f'{bulk_name}_{\"\".join(map(str, facet))}_{\"x\".join(map(str, supercell))}'\n", + "out_dir = Path('surfaces') / output_name\n", + "out_dir.mkdir(parents=True)" + ] + }, + { + "cell_type": "markdown", + "id": "bba219d4-a66b-4a3f-ab36-5ac40359d449", + "metadata": {}, + "source": [ + "## Load in the Structure\n", + "Into an ase.Atoms object" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "8567b17e-7154-43de-a70e-b12897af635a", + "metadata": { + "tags": [] + }, + "outputs": [ + { + "data": { + "text/plain": [ + "Atoms(symbols='Pd', pbc=True, cell=[[1.0167062034150769e-07, 1.9307360191576008, 1.9307356679659875], [1.930735554740657, 1.0167063715463548e-07, 1.9307356679659706], [1.9307355547405736, 1.9307360191575014, 1.0167072120474846e-07]], forces=..., calculator=SinglePointCalculator(...))" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "atoms = read(bulk_path)\n", + "atoms" + ] + }, + { + "cell_type": "markdown", + "id": "d613bae2-50da-4ab1-a84f-a5b3bb708917", + "metadata": {}, + "source": [ + "## Make the Surface Generator\n", + "CatKit uses a [SlabGenerator](https://catkit-jboes.readthedocs.io/en/latest/_static/frontmatter/catgen.html) object to produce surfaces and detect adsorption sites. \n", + "It takes the desired surface geometry as inputs" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "f9a910ee-b206-4517-afc9-9bee246d1548", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "gen = SlabGenerator(\n", + " atoms,\n", + " miller_index=facet,\n", + " layers=layers,\n", + " fixed=layers // 2, # Fix the bottom half\n", + " layer_type='trim', # Get the desired number of layers\n", + " standardize_bulk=True,\n", + " vacuum=10, # Will expand the vacuum layer\n", + " tol=1e-20,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "ff092911-3a66-4c86-a2f9-d8176ae140c2", + "metadata": {}, + "source": [ + "## Save each slab\n", + "Save the unrelaxed slab and coordinates for the absorption sites" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "b820c114-4dce-4386-bfab-cb020f62df87", + "metadata": { + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found 1 unique terminations\n" + ] + } + ], + "source": [ + "terminations = gen.get_unique_terminations()\n", + "print(f'Found {len(terminations)} unique terminations')" + ] + }, + { + "cell_type": "markdown", + "id": "cdd541b6-c750-4c13-a122-c9a1bc6c0ac2", + "metadata": {}, + "source": [ + "Iterate over the terminations" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "4da6c9a9-f333-464f-9123-65b1fdf59724", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "for term in range(len(terminations)):\n", + " # Create the slab then detect the adsorption sites\n", + " slab = gen.get_slab(size=supercell, iterm=term)\n", + " sites = AdsorptionSites(slab)\n", + " \n", + " # Get both the coordinates and vectors for each sites\n", + " coords = sites.get_coordinates()\n", + " vectors = sites.get_adsorption_vectors()\n", + " conn = sites.get_connectivity()\n", + " \n", + " # Save the slab and site information\n", + " slab_dir = out_dir / f'term={term}'\n", + " slab_dir.mkdir()\n", + " \n", + " slab.write(slab_dir / 'unrelaxed.extxyz', columns=['symbols', 'positions', 'move_mask'])\n", + " with slab_dir.joinpath('site-info.json').open('w') as fp:\n", + " json.dump({\n", + " 'coords': coords.tolist(),\n", + " 'vectors': vectors.tolist(),\n", + " 'connectivity': conn.tolist(),\n", + " }, fp)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0ef8dbbe-d2cb-4e30-9629-bed38a0d76b9", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/scripts/lohc/0_make-surfaces/1_relax-surfaces/2_relax-surfaces.py b/scripts/lohc/0_make-surfaces/1_relax-surfaces/2_relax-surfaces.py new file mode 100644 index 000000000..d57d98cd1 --- /dev/null +++ b/scripts/lohc/0_make-surfaces/1_relax-surfaces/2_relax-surfaces.py @@ -0,0 +1,57 @@ +"""Relax the surfaces in parallel on HPC""" +from argparse import ArgumentParser +from pathlib import Path +import logging +import sys + +from ase.io import read + +from examol.simulate.ase import ASESimulator +from examol.utils.conversions import write_to_string + +if __name__ == "__main__": + # Make the argument parser + parser = ArgumentParser() + parser.add_argument('--config', default='cp2k_pbe_dzvp', help='Configuration to use for computations') + parser.add_argument('initial', nargs='+', help='Directory holding initial structures to relax', type=Path) + args = parser.parse_args() + + # Make the logger + logger = logging.getLogger('main') + handler = logging.StreamHandler(sys.stdout) + handler.setFormatter(logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s')) + logger.addHandler(handler) + logger.setLevel(logging.INFO) + + # Find the structures to relax + all_surfaces = sum([list(Path(i).rglob('unrelaxed.extxyz')) for i in args.initial], []) + logger.info(f'Found {len(all_surfaces)} surfaces') + unrelaxed_surfaces = set( + surface for surface in all_surfaces + if not surface.parent.joinpath('relaxed.extxyz').is_file() + ) + logger.info(f'Found {len(unrelaxed_surfaces)} we have not relaxed yet') + + # Create the simulator + sim = ASESimulator( + scratch_dir='./scratch', + cp2k_command='/home/lward/Software/cp2k-2023.2/exe/local_cuda/cp2k_shell.ssmp', + clean_after_run=False + ) + + # Loop over the surfaces to relax + # TODO (wardlt): Use Parsl to run loop on HPC + for surface in unrelaxed_surfaces: + logger.info(f'Starting to relax {surface.parent}') + # Load the file and write as an extended XYZ file + atoms = read(surface) + xyz = write_to_string(atoms, 'extxyz', columns=['symbols', 'positions', 'move_mask']) + + # Relax the surfaces + run_name = f'{surface.parent.parent.name}-{surface.parent.name}' + result, _, _ = sim.optimize_structure(run_name, xyz, config_name=args.config, charge=0, solvent=None) + + # Save the relaxed structure + out_path = surface.parent.joinpath('relaxed.extxyz') + result.atoms.write(out_path, columns=['symbols', 'positions', 'move_mask']) + logger.info(f'Done with {surface.parent}. Saved result to {out_path}') diff --git a/scripts/lohc/0_make-surfaces/1_relax-surfaces/3_assemble-surface-records.ipynb b/scripts/lohc/0_make-surfaces/1_relax-surfaces/3_assemble-surface-records.ipynb new file mode 100644 index 000000000..af47b3a43 --- /dev/null +++ b/scripts/lohc/0_make-surfaces/1_relax-surfaces/3_assemble-surface-records.ipynb @@ -0,0 +1,187 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "25f2a483-88a2-4638-8965-2d0812137763", + "metadata": {}, + "source": [ + "# Assemble Surface Records\n", + "ExaMol specifies the structure of a surface in a format which combines the surface geometry and information about where to place adsorbates.\n", + "Such data are currently stored in a directory." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "15a03474-b523-48a7-8264-7259eaddb7a9", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "from examol.store.recipes.surface import SurfaceRecord, SurfaceSite\n", + "from pathlib import Path\n", + "from ase.io import read\n", + "import json" + ] + }, + { + "cell_type": "markdown", + "id": "2c1b26d6-c64c-4348-ab3f-a489460c48f8", + "metadata": {}, + "source": [ + "## Make a Conversion Function\n", + "The directory structure holds the name of the surface, the sites are in a \"site-info.json\" file, and the surface slab is in a extXYZ file." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "d8dc425e-5294-42c3-8536-560a3004d64c", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "def convert_directory(path: Path) -> SurfaceRecord:\n", + " \"\"\"Convert a directory into a SurfaceRecord\n", + " \n", + " Args:\n", + " path: Path to directory being converted\n", + " Returns:\n", + " Record describing the surface\n", + " \"\"\"\n", + " \n", + " # Start by processing the sites\n", + " with open(path / 'site-info.json') as fp:\n", + " site_info = json.load(fp)\n", + " sites = [\n", + " SurfaceSite(\n", + " coordinate=coords,\n", + " vector=site_info['vectors'][i]\n", + " ) for i, coords in enumerate(site_info['coords'])\n", + " ]\n", + " \n", + " # Load in the slab\n", + " slab_path = path / 'relaxed.extxyz'\n", + " slab = read(slab_path)\n", + " slab_xyz = slab_path.read_text()\n", + " \n", + " # Make the record\n", + " name = f'{path.parent.name}_{path.name}'\n", + " return SurfaceRecord(\n", + " name=name,\n", + " slab=slab_xyz,\n", + " surface_sites=sites,\n", + " energy=slab.get_potential_energy()\n", + " )" + ] + }, + { + "cell_type": "markdown", + "id": "ca527918-9cfa-4ac2-81a5-bc10bb9affec", + "metadata": {}, + "source": [ + "## Convert all relaxed geometries\n", + "Find all relaxed geometries and store them" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "9f3d46d5-b741-4816-baa6-0cf780f61f84", + "metadata": { + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found 2 completed surfaces\n" + ] + } + ], + "source": [ + "record_paths = list(Path('surfaces/').rglob('relaxed.extxyz'))\n", + "print(f'Found {len(record_paths)} completed surfaces')" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "0e86da13-c72e-4d74-83bd-8985b12f6cb8", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "records = [convert_directory(p.parent) for p in record_paths]" + ] + }, + { + "cell_type": "markdown", + "id": "d1534063-7d5e-4d2a-9264-aee536228a9e", + "metadata": {}, + "source": [ + "Save the records" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "6403e5e0-d728-47af-9ea5-50b96e02354c", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "record_dir = Path('records')\n", + "record_dir.mkdir(exist_ok=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "95460110-b018-49ff-89ba-85157d1f4042", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "for record in records:\n", + " with open(record_dir / f'{record.name}.json', 'w') as fp:\n", + " print(record.json(indent=2), file=fp)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "66928624-c7f4-48a9-8d80-e385274079f4", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/scripts/lohc/0_make-surfaces/1_relax-surfaces/README.md b/scripts/lohc/0_make-surfaces/1_relax-surfaces/README.md new file mode 100644 index 000000000..8a011cf26 --- /dev/null +++ b/scripts/lohc/0_make-surfaces/1_relax-surfaces/README.md @@ -0,0 +1,11 @@ +# Pipeline for Creating the Starting Surfaces + +The result of these notebooks is a starting point for adsorption calculations. +The starting input is a bulk crystal, which we then: + +1. Relax using CP2K with PBE configured to ExaMol's standard settings +2. Produce slab structures and enumerate adsorbate sites with CatKit +3. Relax the slab structures with CP2K +4. Assemble a ExaMol-format description of the structures + +Each notebook uses CatKit to produce surface slabs and a find surface sites. \ No newline at end of file diff --git a/scripts/lohc/0_make-surfaces/1_relax-surfaces/records/fcc-Pd_100_3x3_term=0.json b/scripts/lohc/0_make-surfaces/1_relax-surfaces/records/fcc-Pd_100_3x3_term=0.json new file mode 100644 index 000000000..3087e7aff --- /dev/null +++ b/scripts/lohc/0_make-surfaces/1_relax-surfaces/records/fcc-Pd_100_3x3_term=0.json @@ -0,0 +1,49 @@ +{ + "name": "fcc-Pd_100_3x3_term=0", + "slab": "54\nLattice=\"8.191417821843899 0.0 0.0 0.0 8.191417821843899 0.0 0.0 0.0 29.6536784822636\" Properties=species:S:1:pos:R:3:move_mask:L:1 fixed=\"0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26\" energy=-186963.64655061 stress=\"-0.045863123937479 1.5590786751008e-10 3.3470345317445e-10 1.5590786751008e-10 -0.045863123943449 3.4244654075271e-10 3.3470345317445e-10 3.4244654075271e-10 -0.00020308841365442\" free_energy=-186963.64655061 pbc=\"T T T\"\nPd 1.36523630 1.36523630 10.00000000 F\nPd 4.09570891 1.36523630 10.00000000 F\nPd 6.82618152 1.36523630 10.00000000 F\nPd 1.36523630 4.09570891 10.00000000 F\nPd 4.09570891 4.09570891 10.00000000 F\nPd 6.82618152 4.09570891 10.00000000 F\nPd 1.36523630 6.82618152 10.00000000 F\nPd 4.09570891 6.82618152 10.00000000 F\nPd 6.82618152 6.82618152 10.00000000 F\nPd 0.00000000 0.00000000 11.93073570 F\nPd 2.73047261 0.00000000 11.93073570 F\nPd 5.46094521 0.00000000 11.93073570 F\nPd 0.00000000 2.73047261 11.93073570 F\nPd 2.73047261 2.73047261 11.93073570 F\nPd 5.46094521 2.73047261 11.93073570 F\nPd 0.00000000 5.46094521 11.93073570 F\nPd 2.73047261 5.46094521 11.93073570 F\nPd 5.46094521 5.46094521 11.93073570 F\nPd 1.36523630 1.36523630 13.86147139 F\nPd 4.09570891 1.36523630 13.86147139 F\nPd 1.36523630 4.09570891 13.86147139 F\nPd 4.09570891 4.09570891 13.86147139 F\nPd 1.36523630 6.82618152 13.86147139 F\nPd 4.09570891 6.82618152 13.86147139 F\nPd 6.82618152 1.36523630 13.86147139 F\nPd 6.82618152 4.09570891 13.86147139 F\nPd 6.82618152 6.82618152 13.86147139 F\nPd 0.00000016 0.00000016 15.76359357 T\nPd 2.73047114 0.00000016 15.76359658 T\nPd 5.46094655 0.00000015 15.76359660 T\nPd 0.00000016 2.73047114 15.76359659 T\nPd 2.73047115 2.73047115 15.76359961 T\nPd 5.46094655 2.73047114 15.76359961 T\nPd 0.00000016 5.46094655 15.76359659 T\nPd 2.73047115 5.46094655 15.76359961 T\nPd 5.46094655 5.46094655 15.76359962 T\nPd 1.36523217 1.36523217 17.69367434 T\nPd 4.09570891 1.36523217 17.69367876 T\nPd 6.82618555 1.36523217 17.69367436 T\nPd 1.36523218 4.09570891 17.69367876 T\nPd 4.09570891 4.09570891 17.69368319 T\nPd 6.82618555 4.09570891 17.69367878 T\nPd 1.36523218 6.82618555 17.69367435 T\nPd 4.09570891 6.82618555 17.69367877 T\nPd 6.82618555 6.82618555 17.69367438 T\nPd 0.00000004 0.00000004 19.62929249 T\nPd 2.73047065 0.00000004 19.62927906 T\nPd 5.46094715 0.00000004 19.62927900 T\nPd 0.00000004 2.73047065 19.62927906 T\nPd 2.73047065 2.73047065 19.62926564 T\nPd 5.46094716 2.73047065 19.62926557 T\nPd 0.00000004 5.46094715 19.62927900 T\nPd 2.73047065 5.46094716 19.62926557 T\nPd 5.46094716 5.46094716 19.62926552 T\n", + "surface_sites": [ + { + "name": null, + "height": 2, + "coordinate": [ + 8.510800900557517e-17, + -1.2886862146324043e-33, + 19.6536784822636 + ], + "vector": [ + 5.835271171915818e-16, + -0.0, + 1.0 + ] + }, + { + "name": null, + "height": 2, + "coordinate": [ + 2.7304726072812997, + 1.3652363036406499, + 19.6536784822636 + ], + "vector": [ + -2.2120226672260502e-15, + 3.552713678800501e-15, + 1.0 + ] + }, + { + "name": null, + "height": 2, + "coordinate": [ + 1.3652363036406499, + 1.3652363036406499, + 19.6536784822636 + ], + "vector": [ + -4.207674210941454e-15, + -3.552713678800501e-15, + 1.0 + ] + } + ], + "energy": -186963.64655061 +} diff --git a/scripts/lohc/0_make-surfaces/1_relax-surfaces/records/fcc-Pd_111_3x3_term=0.json b/scripts/lohc/0_make-surfaces/1_relax-surfaces/records/fcc-Pd_111_3x3_term=0.json new file mode 100644 index 000000000..83363be21 --- /dev/null +++ b/scripts/lohc/0_make-surfaces/1_relax-surfaces/records/fcc-Pd_111_3x3_term=0.json @@ -0,0 +1,63 @@ +{ + "name": "fcc-Pd_111_3x3_term=0", + "slab": "54\nLattice=\"8.191417821843899 0.0 0.0 4.0957089109219496 7.09397592672941 0.0 0.0 0.0 31.14710774080997\" Properties=species:S:1:pos:R:3:move_mask:L:1 energy=-186955.46912365 stress=\"-0.076773741561668 -1.8655965556549e-08 -7.391151849393e-09 -1.8655965556549e-08 -0.076773719134593 -4.3429999429157e-09 -7.391151849393e-09 -4.3429999429157e-09 0.0072783409049177\" free_energy=-186955.46912365 pbc=\"T T T\"\nPd 8.19141782 1.57643909 10.00000000 F\nPd 2.73047261 1.57643909 10.00000000 F\nPd 5.46094521 1.57643909 10.00000000 F\nPd 4.09570891 3.94109774 10.00000000 F\nPd 6.82618152 3.94109774 10.00000000 F\nPd 9.55665413 3.94109774 10.00000000 F\nPd 5.46094521 6.30575638 10.00000000 F\nPd 8.19141782 6.30575638 10.00000000 F\nPd 10.92189043 6.30575638 10.00000000 F\nPd 1.36523630 0.78821955 12.22942155 F\nPd 4.09570891 0.78821955 12.22942155 F\nPd 6.82618152 0.78821955 12.22942155 F\nPd 2.73047261 3.15287819 12.22942155 F\nPd 5.46094521 3.15287819 12.22942155 F\nPd 8.19141782 3.15287819 12.22942155 F\nPd 4.09570891 5.51753683 12.22942155 F\nPd 6.82618152 5.51753683 12.22942155 F\nPd 9.55665413 5.51753683 12.22942155 F\nPd 5.46094521 -0.00000000 14.45884310 F\nPd 2.73047261 -0.00000000 14.45884310 F\nPd 0.00000000 -0.00000000 14.45884310 F\nPd 1.36523630 2.36465864 14.45884310 F\nPd 4.09570891 2.36465864 14.45884310 F\nPd 2.73047261 4.72931728 14.45884310 F\nPd 6.82618152 2.36465864 14.45884310 F\nPd 5.46094521 4.72931728 14.45884310 F\nPd 8.19141782 4.72931728 14.45884310 F\nPd 2.73047294 1.57643931 16.70083769 T\nPd 5.46094553 1.57643935 16.70083788 T\nPd 8.19141808 1.57643918 16.70083820 T\nPd 4.09570926 3.94109785 16.70083791 T\nPd 6.82618179 3.94109792 16.70083816 T\nPd 9.55665440 3.94109782 16.70083837 T\nPd 5.46094548 6.30575655 16.70083817 T\nPd 8.19141802 6.30575653 16.70083838 T\nPd 10.92189052 6.30575644 16.70083858 T\nPd 1.36523624 0.78821950 18.87857086 T\nPd 4.09570875 0.78821945 18.87857116 T\nPd 6.82618139 0.78821942 18.87857102 T\nPd 2.73047240 3.15287812 18.87857138 T\nPd 5.46094493 3.15287802 18.87857191 T\nPd 8.19141757 3.15287805 18.87857142 T\nPd 4.09570874 5.51753681 18.87857083 T\nPd 6.82618125 5.51753671 18.87857161 T\nPd 9.55665393 5.51753669 18.87857143 T\nPd 5.46094486 -0.00000006 21.04388067 T\nPd -0.00000029 -0.00000017 21.04388102 T\nPd 2.73047233 -0.00000011 21.04388102 T\nPd 1.36523578 2.36465829 21.04388099 T\nPd 4.09570829 2.36465829 21.04388107 T\nPd 6.82618139 2.36465832 21.04388072 T\nPd 2.73047265 4.72931718 21.04388068 T\nPd 5.46094461 4.72931722 21.04388071 T\nPd 8.19141763 4.72931715 21.04388032 T\n", + "surface_sites": [ + { + "name": null, + "height": 2, + "coordinate": [ + 5.460945214562599, + -1.1102227400285158e-16, + 21.14710774080997 + ], + "vector": [ + -1.6416554233197544e-15, + -0.0, + 1.0 + ] + }, + { + "name": null, + "height": 2, + "coordinate": [ + 0.6826181518203249, + 1.1823293211215682, + 21.14710774080997 + ], + "vector": [ + -8.6051042102145e-15, + -3.552713678800501e-15, + 1.0 + ] + }, + { + "name": null, + "height": 2, + "coordinate": [ + 1.3652363036406499, + 0.7882195474143788, + 21.14710774080997 + ], + "vector": [ + 3.914107470955258e-15, + 3.552713678800501e-15, + 1.0 + ] + }, + { + "name": null, + "height": 2, + "coordinate": [ + 2.7304726072812997, + 1.5764390948287577, + 21.14710774080997 + ], + "vector": [ + -7.182243703660577e-16, + -0.0, + 1.0 + ] + } + ], + "energy": -186955.46912365 +} diff --git a/scripts/lohc/0_make-surfaces/README.md b/scripts/lohc/0_make-surfaces/README.md new file mode 100644 index 000000000..08a9a5f84 --- /dev/null +++ b/scripts/lohc/0_make-surfaces/README.md @@ -0,0 +1,8 @@ +# Seed the Information Needed for Absorption Computations + +The following directories are used to prepare for surface energy computations. + +1. _Optional_: Verify that CP2K is using converged parameters for our target surfaces. + - _TBD_: Update the documentation on using CP2K with PBE. +2. Generate then relax the surfaces. +3. Compute surface adsorption for certain species. \ No newline at end of file diff --git a/scripts/redoxmers/0_check-chemistry-settings/compile-dataset.ipynb b/scripts/redoxmers/0_check-chemistry-settings/compile-dataset.ipynb index 857ff8476..c133ed718 100644 --- a/scripts/redoxmers/0_check-chemistry-settings/compile-dataset.ipynb +++ b/scripts/redoxmers/0_check-chemistry-settings/compile-dataset.ipynb @@ -24,7 +24,7 @@ "outputs": [], "source": [ "from examol.store.models import MoleculeRecord\n", - "from examol.store.recipes import RedoxEnergy, SolvationEnergy\n", + "from examol.store.recipes.redox import RedoxEnergy, SolvationEnergy\n", "from base64 import b64decode\n", "from pathlib import Path\n", "from tqdm import tqdm\n", diff --git a/scripts/redoxmers/2_initial-data/0_gather-initial-dataset.py b/scripts/redoxmers/2_initial-data/0_gather-initial-dataset.py index bece720ce..a0de9aa6e 100644 --- a/scripts/redoxmers/2_initial-data/0_gather-initial-dataset.py +++ b/scripts/redoxmers/2_initial-data/0_gather-initial-dataset.py @@ -21,7 +21,7 @@ from colmena.models import Result from examol.store.models import MoleculeRecord -from examol.store.recipes import RedoxEnergy, SolvationEnergy +from examol.store.recipes.redox import RedoxEnergy, SolvationEnergy import configs diff --git a/scripts/redoxmers/2_initial-data/configs.py b/scripts/redoxmers/2_initial-data/configs.py index a581123f2..d21e11497 100644 --- a/scripts/redoxmers/2_initial-data/configs.py +++ b/scripts/redoxmers/2_initial-data/configs.py @@ -141,7 +141,11 @@ def make_polaris_config() -> tuple[Config, ASESimulator, int, list[str]]: ) sim = ASESimulator( scratch_dir='cp2k-files', +<<<<<<< HEAD + optimization_steps=100, +======= optimization_steps=150, +>>>>>>> main cp2k_command=f'mpiexec -n {nodes_per_cp2k * 4} --ppn 4 --cpu-bind depth --depth 8 -env OMP_NUM_THREADS=8 ' f'--hostfile /tmp/hostfiles/local_hostfile.`printf %02d $PARSL_WORKER_RANK` ' '/lus/grand/projects/CSC249ADCD08/cp2k/set_affinity_gpu_polaris.sh ' diff --git a/scripts/redoxmers/3_initial-models/0_pull-test-datasets.ipynb b/scripts/redoxmers/3_initial-models/0_pull-test-datasets.ipynb index 30e482444..592eacad6 100644 --- a/scripts/redoxmers/3_initial-models/0_pull-test-datasets.ipynb +++ b/scripts/redoxmers/3_initial-models/0_pull-test-datasets.ipynb @@ -42,7 +42,7 @@ "outputs": [], "source": [ "dataset_name: str = 'mdf-mos'\n", - "target_prop: str = 'reduction_potential'\n", + "target_prop: str = 'oxidation_potential'\n", "level: str = 'mopac_pm7-acn-adiabatic'" ] }, @@ -88,14 +88,14 @@ "name": "stderr", "output_type": "stream", "text": [ - "1115110it [00:31, 35706.54it/s]" + "1115110it [00:30, 36128.20it/s]" ] }, { "name": "stdout", "output_type": "stream", "text": [ - "Loaded 150585 matching records. Hash: b29db032ead2c10d02d23947a5f452c8\n" + "Loaded 104122 matching records. Hash: 919289d96001c62a986bd1fedae18a5c\n" ] }, { @@ -171,7 +171,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "Split off 112938 for a training set\n" + "Split off 78091 for a training set\n" ] } ], @@ -190,7 +190,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": null, "id": "dac90dfd-9feb-4d1d-b4ae-a9fa5c6b5064", "metadata": { "tags": [] @@ -200,8 +200,8 @@ "name": "stderr", "output_type": "stream", "text": [ - "100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 112938/112938 [06:07<00:00, 306.95it/s]\n", - "100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 37647/37647 [02:01<00:00, 310.27it/s]\n" + "100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 78091/78091 [04:13<00:00, 308.61it/s]\n", + " 85%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▏ | 22224/26031 [01:12<00:12, 311.20it/s]" ] } ], diff --git a/scripts/redoxmers/3_initial-models/nfp-mpnn/0_run-test.py b/scripts/redoxmers/3_initial-models/nfp-mpnn/0_run-test.py index 30a23e146..07ecb6cf1 100644 --- a/scripts/redoxmers/3_initial-models/nfp-mpnn/0_run-test.py +++ b/scripts/redoxmers/3_initial-models/nfp-mpnn/0_run-test.py @@ -1,5 +1,5 @@ """Train a model and save results into a directory structure""" -from examol.store.recipes import RedoxEnergy, PropertyRecipe +from examol.store.recipes.redox import RedoxEnergy, PropertyRecipe from examol.score.nfp import NFPScorer, make_simple_network from examol.store.models import MoleculeRecord from sklearn import metrics diff --git a/scripts/redoxmers/3_initial-models/nfp-multifi-mpnn/0_run-test.py b/scripts/redoxmers/3_initial-models/nfp-multifi-mpnn/0_run-test.py index f8e12739e..0e56fb1df 100644 --- a/scripts/redoxmers/3_initial-models/nfp-multifi-mpnn/0_run-test.py +++ b/scripts/redoxmers/3_initial-models/nfp-multifi-mpnn/0_run-test.py @@ -1,5 +1,5 @@ """Train a model and save results into a directory structure""" -from examol.store.recipes import RedoxEnergy +from examol.store.recipes.redox import RedoxEnergy from examol.score.nfp import NFPScorer, make_simple_network from examol.store.models import MoleculeRecord from sklearn import metrics diff --git a/scripts/redoxmers/4_run-examol/spec.py b/scripts/redoxmers/4_run-examol/spec.py index 0f822aae5..7921f2d17 100644 --- a/scripts/redoxmers/4_run-examol/spec.py +++ b/scripts/redoxmers/4_run-examol/spec.py @@ -16,8 +16,8 @@ from examol.select.bayes import ExpectedImprovement from examol.simulate.ase import ASESimulator from examol.steer.single import SingleStepThinker -from examol.store.recipes import RedoxEnergy -from examol.specify.spec import ExaMolSpecification +from examol.store.recipes.redox import RedoxEnergy +from examol.specify import ExaMolSpecification logger = logging.getLogger('spec') diff --git a/tests/report/test_reporting.py b/tests/report/test_reporting.py index d90cd3cc9..2e105265b 100644 --- a/tests/report/test_reporting.py +++ b/tests/report/test_reporting.py @@ -8,7 +8,7 @@ from examol.reporting.markdown import MarkdownReporter from examol.steer.single import SingleStepThinker -from examol.store.recipes import RedoxEnergy +from examol.store.recipes.redox import RedoxEnergy example_dir = Path(__file__).parent / 'example' diff --git a/tests/score/conftest.py b/tests/score/conftest.py index 65c7c3c86..0d6b95bc6 100644 --- a/tests/score/conftest.py +++ b/tests/score/conftest.py @@ -1,7 +1,7 @@ from pytest import fixture from examol.store.models import MoleculeRecord -from examol.store.recipes import PropertyRecipe +from examol.store.recipes.base import PropertyRecipe class FakeRecipe(PropertyRecipe): diff --git a/tests/select/conftest.py b/tests/select/conftest.py index b2098d803..86ad1f68c 100644 --- a/tests/select/conftest.py +++ b/tests/select/conftest.py @@ -2,7 +2,7 @@ import numpy as np from examol.store.models import MoleculeRecord -from examol.store.recipes import PropertyRecipe +from examol.store.recipes.base import PropertyRecipe class TestRecipe(PropertyRecipe): diff --git a/tests/simulation/test_ase.py b/tests/simulation/test_ase.py index 42e74c3e2..4295d49f4 100644 --- a/tests/simulation/test_ase.py +++ b/tests/simulation/test_ase.py @@ -6,11 +6,11 @@ import numpy as np from ase import units -from ase.calculators.gaussian import Gaussian from ase.db import connect -from pytest import mark, fixture, raises, param -from ase.build import molecule +from ase.build import molecule, bulk from ase.calculators.lj import LennardJones +from ase.calculators.gaussian import Gaussian +from pytest import mark, fixture, raises, param from examol.simulate.ase import ASESimulator from examol.simulate.ase.utils import make_ephemeral_calculator @@ -19,6 +19,7 @@ try: import xtb # noqa: F401 + has_xtb = True except ImportError: has_xtb = False @@ -67,6 +68,12 @@ def test_cp2k_configs(tmpdir, strc): assert 'GAPW' in config['kwargs']['inp'] assert Path(config['kwargs']['basis_set_file']).is_file() + # With PBE + config = sim.create_configuration('cp2k_pbe_dzvp', strc, charge=0, solvent=None) + assert config['kwargs']['cutoff'] == 700 * units.Ry + assert 'PBE' in config['kwargs']['inp'] + assert 'PERIODIC XYZ' in config['kwargs']['inp'] + # With an undefined basis set with raises(AssertionError): sim.create_configuration('cp2k_blyp_notreal', strc, charge=1, solvent=None) @@ -81,7 +88,8 @@ def test_cp2k_configs(tmpdir, strc): + [(xc, 0, 'acn') for xc in cp2k_configs_to_test] # With a solvent + [(cp2k_configs_to_test[0], -1, 'acn')] # Open shell and a solvent ) -def test_ase_singlepoint(tmpdir, strc, config, charge, solvent): +def test_cp2k_0d_singlepoint(tmpdir, strc, config, charge, solvent): + """Test the computations w/o periodic boundary conditions.""" sim = ASESimulator(scratch_dir=tmpdir) sim.compute_energy('test', strc, config_name=config, charge=charge, solvent=solvent) @@ -246,3 +254,18 @@ def test_opt_failure(tmpdir): opt, steps, _ = sim.optimize_structure('test', strc, 'mopac_pm7', charge=0, solvent=None) assert np.max(opt.forces) <= 0.02, 'Optimization did not finish successfully' assert len(steps) < 100 + + +@mark.skipif(not has_cpk2, reason='Requires CP2K') +def test_3d_structures(tmpdir): + """Test relaxing a 3D structure""" + + # Make a primitive cell of Aluminum and convert it to extxyz + atoms = bulk('Al') + xyz = write_to_string(atoms, 'extxyz') + + # Run a relaxation + sim = ASESimulator(scratch_dir='./3d-test') + result, _, _ = sim.optimize_structure('Al', xyz, config_name='cp2k_pbe_dzvp', charge=0, solvent=None) + assert np.max(result.forces) < 0.01 + assert all(result.atoms.pbc) diff --git a/tests/simulation/test_problem_recipes.py b/tests/simulation/test_problem_recipes.py index 952b45962..230db2ae6 100644 --- a/tests/simulation/test_problem_recipes.py +++ b/tests/simulation/test_problem_recipes.py @@ -6,7 +6,7 @@ from examol.simulate.ase import ASESimulator from examol.store.models import MoleculeRecord -from examol.store.recipes import RedoxEnergy +from examol.store.recipes.redox import RedoxEnergy is_mac = platform.startswith("darwin") diff --git a/tests/specify/test_spec.py b/tests/specify/test_spec.py index c00a04ef9..d44db5119 100644 --- a/tests/specify/test_spec.py +++ b/tests/specify/test_spec.py @@ -18,7 +18,7 @@ from examol.steer.single import SingleStepThinker from examol.store.db.memory import InMemoryStore from examol.store.models import MoleculeRecord -from examol.store.recipes import RedoxEnergy +from examol.store.recipes.redox import RedoxEnergy @fixture() diff --git a/tests/steer/conftest.py b/tests/steer/conftest.py index 685b2efb0..34f31187f 100644 --- a/tests/steer/conftest.py +++ b/tests/steer/conftest.py @@ -14,7 +14,7 @@ from examol.store.db.base import MoleculeStore from examol.store.db.memory import InMemoryStore from examol.store.models import MoleculeRecord -from examol.store.recipes import RedoxEnergy +from examol.store.recipes.redox import RedoxEnergy @fixture() diff --git a/tests/store/test_models.py b/tests/store/test_models.py index f893b468a..a639db57b 100644 --- a/tests/store/test_models.py +++ b/tests/store/test_models.py @@ -140,3 +140,32 @@ def test_match_opt_only(record, sim_result): # Getting the conformer w/o filtering should be the original energy, filtering should yield the new higher energy assert record.find_lowest_conformer(sim_result.config_name, 1, None, optimized_only=False)[1] == original_energy assert record.find_lowest_conformer(sim_result.config_name, 1, None, optimized_only=True)[1] == sim_result.energy + + +def test_best_xyz(record, sim_result): + # Test generating an XYZ + record.conformers.clear() # Start with nothing + conf, new_xyz = record.find_closest_xyz('test', 0) + assert conf is None + + # Add some conformers + record.conformers = [ + Conformer.from_xyz(new_xyz, config_name='test', charge=0), + Conformer.from_xyz(new_xyz, config_name='test', charge=1), + Conformer.from_xyz(new_xyz, config_name='not_test', charge=-1), + Conformer.from_xyz(new_xyz, config_name='not_test', charge=0) + ] + + # We should match to the config_name over charge, and get the closest charge + matched, _ = record.find_closest_xyz(config_name='test', charge=-1) + assert matched.config_name == 'test' + assert matched.charge == 0 + + matched, _ = record.find_closest_xyz(config_name='test', charge=2) + assert matched.config_name == 'test' + assert matched.charge == 1 + + # We should match to the newest if several have the same difference in charge and none match config_name + matched, _ = record.find_closest_xyz(config_name='another', charge=0) + assert matched.config_name == 'not_test' + assert matched.charge == 0 diff --git a/tests/store/test_recipe.py b/tests/store/test_recipe.py index 73c6e3e70..496396c3d 100644 --- a/tests/store/test_recipe.py +++ b/tests/store/test_recipe.py @@ -6,7 +6,8 @@ from ase import units from examol.store.models import MissingData -from examol.store.recipes import SolvationEnergy, RedoxEnergy +from examol.store.recipes.basic import SolvationEnergy +from examol.store.recipes.redox import RedoxEnergy def test_solvation(record, sim_result):