Add support for computing adsorption energies

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:

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.
 

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)}

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 <components/simulate.html#levels>`_.
 
-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)?

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:

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
 

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__)
 

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):

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):

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

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"""

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__)
 

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):

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):

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):

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

examol/store/recipes/__init__.py

@@ -0,0 +1 @@
+"""Tools for computing the properties of molecules from their record"""