From 0c696bc03111325619a60d9d1a78e400340c96ff Mon Sep 17 00:00:00 2001 From: Patrick Melix Date: Wed, 19 Jun 2024 13:54:22 +0200 Subject: [PATCH 01/29] first draft of bandfragment and nocv multijob recipes --- recipes/bandfragment.py | 102 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 102 insertions(+) create mode 100644 recipes/bandfragment.py diff --git a/recipes/bandfragment.py b/recipes/bandfragment.py new file mode 100644 index 00000000..2056442f --- /dev/null +++ b/recipes/bandfragment.py @@ -0,0 +1,102 @@ +from scm.plams.core.functions import log +from scm.plams.interfaces.adfsuite.ams import AMSJob +from scm.plams.core.settings import Settings +from scm.plams.recipes.adffragment import ADFFragmentJob, ADFFragmentResults +from scm.plams.tools.units import Units + +from os.path import join as opj +from typing import Dict + +# Check if ASE is present as it is not a default dependency of PLAMS. +__all__ = ['BANDFragmentJob', 'BandFragmentResults'] + + +class BandFragmentResults(ADFFragmentResults): + + def get_energy_decomposition(self, unit='kJ/mol') -> Dict[str, float]: + res=self.job.full.results + res1=self.job.f2.results + res2=self.job.f1.results + + res = {} + + pos = 2 + res['E_int'] = Units.convert(float(res.grep_output('E_int')[-1].split()[pos]), 'au', unit) + res['E_disp'] = Units.convert(float(res.grep_output('E_disp')[-1].split()[pos]), 'au', unit) + res['E_Pauli'] = Units.convert(float(res.grep_output('E_Pauli')[-1].split()[pos]), 'au', unit) + res['E_elstat'] = Units.convert(float(res.grep_output('E_elstat')[-1].split()[pos]), 'au', unit) + res['E_orb'] = Units.convert(float(res.grep_output('E_orb')[-1].split()[pos]), 'au', unit) + + res['E_1'] = res1.get_energy(unit=unit) + res['E_2'] = res2.get_energy(unit=unit) + + return res + + +class BANDFragmentJob(ADFFragmentJob): + _result_type = BandFragmentResults + + def create_mapping_setting(self): + if "fragment" in self.full_settings.input.band: + log("Fragment already present in full_settings. Assuming that the user has already set up the mapping. Skipping the mapping setup.", level=1) + return + set1 = Settings() + set1.filename = opj(self.f1.path,'band.rkf') + # if no mapping is saved, generate a basic mapping, + # first fragment 1 then fragment 2 + set2 = Settings() + set2.filename = opj(self.f2.path,'band.rkf') + if self.fragment_mapping: + set1.atommapping = { i+1: i+1 for i in range(len(self.fragment1)) } + set2.atommapping = { i+1: i+1+len(self.fragment1) for i in range(len(self.fragment2)) } + else: + set1.atommapping = self.fragment_mapping[0].copy() + set1.atommapping = self.fragment_mapping[1].copy() + self.full_settings.input.band.fragment = [set1, set2] + + + def prerun(self): + self.f1 = AMSJob(name='frag1', molecule=self.fragment1, settings=self.settings) + self.f2 = AMSJob(name='frag2', molecule=self.fragment2, settings=self.settings) + self.children=[self.f1,self.f2] + + # create the correct mapping settings for the full job + self.create_mapping_setting() + + self.full = AMSJob(name = 'full', + molecule = self.fragment1 + self.fragment2, + settings = self.settings + self.full_settings) + self.full.depend += [self.f1,self.f2] + + self.children.append(self.full) + + +class NOCVBandFragmentJob(BANDFragmentJob): + _result_type = BandFragmentResults + + def __init__(self, nocv_settings=None, **kwargs): + """Create a BANDFragmentJob with NOCV calculation. + + Args: + nocv_settings (Settings, optional): Settings for the NOCV calculation. Defaults to None. + + """ + super().__init__(**kwargs) + self.nocv_settings = nocv_settings or Settings() + + def prerun(self): + # prerun normal BANDFragmentJob + super().prerun() + # add NOCV run + self.nocv=AMSJob(name="NOCV", molecule= self.fragment1 + self.fragment2, + settings = self.settings + self.full_settings+self.nocv_settings) + # make sure we at least have the restart section to put the file name + if not hasattr(self.nocv.settings, 'input'): + self.nocv.settings.input = Settings() + if not hasattr(self.nocv.settings.input, 'band'): + self.nocv.settings.input.band = Settings() + if not hasattr(self.nocv.settings.input.band, 'restart'): + self.nocv.settings.input.band.restart = Settings() + if not hasattr(self.nocv.settings.input.band.restart, 'file'): + self.nocv.settings.input.band.Restart.File = opj(self.full.path,'band.rkf') + self.children.append(self.nocv) \ No newline at end of file From 1d2240e60ae797d89004698f954a1cd041b74e80 Mon Sep 17 00:00:00 2001 From: Patrick Melix Date: Wed, 31 Jul 2024 09:40:53 +0200 Subject: [PATCH 02/29] fix depreceated sphinx add_stylesheet --- doc/source/conf.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/doc/source/conf.py b/doc/source/conf.py index 04068765..6ce38ce8 100644 --- a/doc/source/conf.py +++ b/doc/source/conf.py @@ -20,7 +20,7 @@ def modify_signature(app, what, name, obj, options, signature, return_annotation def setup(app): if not tags.has("scm_theme"): - app.add_stylesheet("boxes.css") + app.add_css_file("boxes.css") app.add_directive("warning", Danger) app.add_directive("technical", Important) app.connect("autodoc-process-signature", modify_signature) @@ -325,6 +325,8 @@ def setup(app): .. |VibrationsJob| replace:: :class:`~scm.plams.recipes.vibration.VibrationsJob` .. |IRJob| replace:: :class:`~scm.plams.recipes.vibration.IRJob` .. |VibrationsResults| replace:: :class:`~scm.plams.recipes.vibration.VibrationsResults` +.. |ADFFragmentJob| replace:: :class:`~scm.plams.recipes.adffragment.ADFFragmentJob` +.. |BANDFragmentJob| replace:: :class:`~scm.plams.recipes.bandfragment.BANDFragmentJob` .. |RPM| replace:: :ref:`rerun-prevention` .. |cleaning| replace:: :ref:`cleaning` From 0a5ee75801023fcbe4b857136bd750f691eb3b67 Mon Sep 17 00:00:00 2001 From: Patrick Melix Date: Wed, 31 Jul 2024 10:07:29 +0200 Subject: [PATCH 03/29] more work on bandfragment --- recipes/bandfragment.py | 199 +++++++++++++++++++++++++++++----------- 1 file changed, 144 insertions(+), 55 deletions(-) diff --git a/recipes/bandfragment.py b/recipes/bandfragment.py index 2056442f..e73f0b62 100644 --- a/recipes/bandfragment.py +++ b/recipes/bandfragment.py @@ -3,78 +3,159 @@ from scm.plams.core.settings import Settings from scm.plams.recipes.adffragment import ADFFragmentJob, ADFFragmentResults from scm.plams.tools.units import Units +from scm.plams.core.errors import FileError +from scm.plams.mol.molecule import Molecule + from os.path import join as opj -from typing import Dict +from os.path import relpath, abspath, isdir, basename +from typing import Dict, List, Union -# Check if ASE is present as it is not a default dependency of PLAMS. -__all__ = ['BANDFragmentJob', 'BandFragmentResults'] +__all__ = ['BANDFragmentJob', 'BANDFragmentResults'] -class BandFragmentResults(ADFFragmentResults): +class BANDFragmentResults(ADFFragmentResults): + """Subclass of |ADFFragmentResults| for BAND calculations.""" def get_energy_decomposition(self, unit='kJ/mol') -> Dict[str, float]: + """Get the energy decomposition of the fragment calculation. + + Args: + unit (str, optional): The unit of the energy. Defaults to 'kJ/mol'. + + Returns: + Dict[str, float]: The energy decomposition. + """ res=self.job.full.results res1=self.job.f2.results res2=self.job.f1.results - - res = {} - + ret = {} pos = 2 - res['E_int'] = Units.convert(float(res.grep_output('E_int')[-1].split()[pos]), 'au', unit) - res['E_disp'] = Units.convert(float(res.grep_output('E_disp')[-1].split()[pos]), 'au', unit) - res['E_Pauli'] = Units.convert(float(res.grep_output('E_Pauli')[-1].split()[pos]), 'au', unit) - res['E_elstat'] = Units.convert(float(res.grep_output('E_elstat')[-1].split()[pos]), 'au', unit) - res['E_orb'] = Units.convert(float(res.grep_output('E_orb')[-1].split()[pos]), 'au', unit) - - res['E_1'] = res1.get_energy(unit=unit) - res['E_2'] = res2.get_energy(unit=unit) + # E_int appears in a comment below the PEDA Table of the output + ret['E_int'] = Units.convert(float(res.grep_output('E_int')[-2].split()[pos]), 'au', unit) + ret['E_int_disp'] = Units.convert(float(res.grep_output('E_disp')[-1].split()[pos]), 'au', unit) + ret['E_Pauli'] = Units.convert(float(res.grep_output('E_Pauli')[-1].split()[pos]), 'au', unit) + ret['E_elstat'] = Units.convert(float(res.grep_output('E_elstat')[-1].split()[pos]), 'au', unit) + ret['E_orb'] = Units.convert(float(res.grep_output('E_orb')[-1].split()[pos]), 'au', unit) + + ret['E_1'] = res1.get_energy(unit=unit) + ret['E_2'] = res2.get_energy(unit=unit) + + if hasattr(self.job, 'f1_opt'): + ret['E_1_opt'] = self.job.f1_opt.results.get_energy(unit=unit) + ret['E_prep_1'] = ret['E_1'] - ret['E_1_opt'] + if hasattr(self.job, 'f2_opt'): + ret['E_2_opt'] = self.job.f2_opt.results.get_energy(unit=unit) + ret['E_prep_2'] = ret['E_2'] - ret['E_2_opt'] + + if ('E_1_opt' in ret) and ('E_2_opt' in ret): + ret['E_prep'] = ret['E_prep_1'] + ret['E_prep_2'] + ret['E_bond'] = ret['E_int'] + ret['E_prep'] - return res + return ret class BANDFragmentJob(ADFFragmentJob): - _result_type = BandFragmentResults + _result_type = BANDFragmentResults - def create_mapping_setting(self): + def __init__(self, fragment1_opt: Molecule=None, + fragment2_opt: Molecule=None, **kwargs) -> None: + """Create a BANDFragmentJob with the given fragments. + + The optimized fragment structures can be given as arguments. + Single point calculations on these structures will then be performed + to obtain the preparation energy. + + Subclass of |ADFFragmentJob|. + + Args: + fragment1_opt (Molecule, optional): The optimized fragment 1. + fragment2_opt (Molecule, optional): The optimized fragment 2. + **kwargs: Further keyword arguments for |ADFFragmentJob|. + """ + super().__init__(**kwargs) + if isinstance(fragment1_opt, Molecule): + self.fragment1_opt = fragment1_opt.copy() + self.f1_opt = AMSJob(name='frag1_opt', molecule=self.fragment1_opt, + settings=self.settings) + self.children.append(self.f1_opt) + if isinstance(fragment2_opt, Molecule): + self.fragment2_opt = fragment2_opt.copy() + self.f2_opt = AMSJob(name='frag2_opt', molecule=self.fragment2_opt, + settings=self.settings) + self.children.append(self.f2_opt) + + def create_mapping_setting(self) -> None: if "fragment" in self.full_settings.input.band: log("Fragment already present in full_settings. Assuming that the user has already set up the mapping. Skipping the mapping setup.", level=1) return - set1 = Settings() - set1.filename = opj(self.f1.path,'band.rkf') - # if no mapping is saved, generate a basic mapping, # first fragment 1 then fragment 2 + set1 = Settings() + set1.atommapping = { str(i+1): str(i+1) for i in range(len(self.fragment1)) } set2 = Settings() - set2.filename = opj(self.f2.path,'band.rkf') - if self.fragment_mapping: - set1.atommapping = { i+1: i+1 for i in range(len(self.fragment1)) } - set2.atommapping = { i+1: i+1+len(self.fragment1) for i in range(len(self.fragment2)) } - else: - set1.atommapping = self.fragment_mapping[0].copy() - set1.atommapping = self.fragment_mapping[1].copy() + set2.atommapping = { str(i+1): str(i+1+len(self.fragment1)) for i in range(len(self.fragment2))} + # get the correct restart files, working with relative paths + set1.filename = opj("../", self.f1.name, "band.rkf") + set2.filename = opj("../", self.f2.name, "band.rkf") self.full_settings.input.band.fragment = [set1, set2] - def prerun(self): + def new_children(self) -> Union[None, List[AMSJob]]: + """After the first round, add the full job to the children list.""" + if hasattr(self, 'full'): + return None + else: + # create the correct mapping settings for the full job + self.create_mapping_setting() + # create the full job + self.full = AMSJob(name = 'full', + molecule = self.fragment1 + self.fragment2, + settings = self.settings + self.full_settings) + # dependencies are optional, but let's set them up + self.full.depend += [self.f1,self.f2] + return [self.full] + + + def prerun(self) -> None: + """Creates the fragment jobs.""" self.f1 = AMSJob(name='frag1', molecule=self.fragment1, settings=self.settings) self.f2 = AMSJob(name='frag2', molecule=self.fragment2, settings=self.settings) - self.children=[self.f1,self.f2] + self.children += [self.f1,self.f2] + + + @classmethod + def load_external(cls, path: str, jobname: str=None) -> "BANDFragmentJob": + """Load the results of the BANDFragmentJob job from an external path. + + Args: + path (str): The path to the job. It should at least have the + subfolders 'frag1', 'frag2' and 'full'. + jobname (str, optional): The name of the job. Defaults to None. + + Returns: + BANDFragmentJob: The job with the loaded results. + """ + if not isdir(path): + raise FileError("Path {} does not exist, cannot load from it.".format(path)) + path = abspath(path) + jobname = basename(path) if jobname is None else str(jobname) - # create the correct mapping settings for the full job - self.create_mapping_setting() + job = cls(name=jobname) + job.path = path + job.status = "copied" - self.full = AMSJob(name = 'full', - molecule = self.fragment1 + self.fragment2, - settings = self.settings + self.full_settings) - self.full.depend += [self.f1,self.f2] + job.f1 = AMSJob.load_external(opj(path, 'frag1')) + job.f2 = AMSJob.load_external(opj(path, 'frag2')) + job.full = AMSJob.load_external(opj(path, 'full')) + + return job - self.children.append(self.full) class NOCVBandFragmentJob(BANDFragmentJob): - _result_type = BandFragmentResults + _result_type = BANDFragmentResults - def __init__(self, nocv_settings=None, **kwargs): + def __init__(self, nocv_settings=None, **kwargs) -> None: """Create a BANDFragmentJob with NOCV calculation. Args: @@ -84,19 +165,27 @@ def __init__(self, nocv_settings=None, **kwargs): super().__init__(**kwargs) self.nocv_settings = nocv_settings or Settings() - def prerun(self): - # prerun normal BANDFragmentJob - super().prerun() - # add NOCV run - self.nocv=AMSJob(name="NOCV", molecule= self.fragment1 + self.fragment2, - settings = self.settings + self.full_settings+self.nocv_settings) - # make sure we at least have the restart section to put the file name - if not hasattr(self.nocv.settings, 'input'): - self.nocv.settings.input = Settings() - if not hasattr(self.nocv.settings.input, 'band'): - self.nocv.settings.input.band = Settings() - if not hasattr(self.nocv.settings.input.band, 'restart'): - self.nocv.settings.input.band.restart = Settings() - if not hasattr(self.nocv.settings.input.band.restart, 'file'): - self.nocv.settings.input.band.Restart.File = opj(self.full.path,'band.rkf') - self.children.append(self.nocv) \ No newline at end of file + + def new_children(self) -> Union[None, List[AMSJob]]: + """After the first round, add the full job to the children list. + After the second round, add the NOCV job to the children list.""" + # new_children of BANDFragmentJob creates the full job + ret = super().prerun() + if ret is not None: + return ret + if hasattr(self, 'nocv'): + return None + else: + # add NOCV run + self.nocv=AMSJob(name="NOCV", molecule= self.fragment1 + self.fragment2, + settings = self.settings + self.full_settings+self.nocv_settings) + # make sure we at least have the restart section to put the file name + if not hasattr(self.nocv.settings, 'input'): + self.nocv.settings.input = Settings() + if not hasattr(self.nocv.settings.input, 'band'): + self.nocv.settings.input.band = Settings() + if not hasattr(self.nocv.settings.input.band, 'restart'): + self.nocv.settings.input.band.restart = Settings() + if not hasattr(self.nocv.settings.input.band.restart, 'file'): + self.nocv.settings.input.band.Restart.File = opj("../", self.full.name,'band.rkf') + return [self.nocv] \ No newline at end of file From af4d7614995dbbb68857e787481168ab7957fd08 Mon Sep 17 00:00:00 2001 From: Patrick Melix Date: Wed, 31 Jul 2024 10:07:49 +0200 Subject: [PATCH 04/29] documentation --- doc/source/conf.py | 2 ++ doc/source/examples/adffragment.rst | 8 ++++++-- doc/source/examples/bandfragment.rst | 27 +++++++++++++++++++++++++++ doc/source/examples/examples.rst | 1 + recipes/adffragment.py | 15 +++++++++++++++ 5 files changed, 51 insertions(+), 2 deletions(-) create mode 100644 doc/source/examples/bandfragment.rst diff --git a/doc/source/conf.py b/doc/source/conf.py index 6ce38ce8..b6740b7d 100644 --- a/doc/source/conf.py +++ b/doc/source/conf.py @@ -326,7 +326,9 @@ def setup(app): .. |IRJob| replace:: :class:`~scm.plams.recipes.vibration.IRJob` .. |VibrationsResults| replace:: :class:`~scm.plams.recipes.vibration.VibrationsResults` .. |ADFFragmentJob| replace:: :class:`~scm.plams.recipes.adffragment.ADFFragmentJob` +.. |ADFFragmentResults| replace:: :class:`~scm.plams.recipes.adffragment.ADFFragmentResults` .. |BANDFragmentJob| replace:: :class:`~scm.plams.recipes.bandfragment.BANDFragmentJob` +.. |BANDFragmentResults| replace:: :class:`~scm.plams.recipes.bandfragment.BANDFragmentResults` .. |RPM| replace:: :ref:`rerun-prevention` .. |cleaning| replace:: :ref:`cleaning` diff --git a/doc/source/examples/adffragment.rst b/doc/source/examples/adffragment.rst index 8d09dde8..72cca160 100644 --- a/doc/source/examples/adffragment.rst +++ b/doc/source/examples/adffragment.rst @@ -2,6 +2,7 @@ ADF fragment job -------------------- +.. currentmodule:: scm.plams.recipes.adffragment In this module a dedicated job type for ADF fragment analysis is defined. Such an analysis is performed on a molecular system divided into 2 fragments and consists of 3 separate ADF runs: one for each fragment and one for full system. @@ -15,9 +16,12 @@ They are then added to the ``children`` list. The dedicated |Results| subclass for ``ADFFragmentJob`` does not provide too much additional functionality. It simply redirects the usual |AMSResults| methods to the results of the full system calculation. -The source code of the whole module with both abovementioned classes: -.. literalinclude:: ../../../recipes/adffragment.py +API +~~~~~~~~~~~~~~~~~~~~~~~~~ + +.. autoclass:: ADFFragmentJob() +.. autoclass:: ADFFragmentResults() An example usage: (:download:`ethene.xyz <../../../../../../examples/scripting/plams_adffrag/ethene.xyz>`, :download:`butadiene.xyz <../../../../../../examples/scripting/plams_adffrag/butadiene.xyz>`, diff --git a/doc/source/examples/bandfragment.rst b/doc/source/examples/bandfragment.rst new file mode 100644 index 00000000..9e352ef5 --- /dev/null +++ b/doc/source/examples/bandfragment.rst @@ -0,0 +1,27 @@ +.. _band-fragment-recipe: + +BAND fragment job +-------------------- + +.. currentmodule:: scm.plams.recipes.bandfragment + +In this module a dedicated job type for Energy Decomposition Analysis in BAND is defined. +Such an analysis is performed on a periodic system divided into 2 fragments and consists of a minimum of 3 separate Band runs: one for each fragment and one for full system. See + +We define a new job type |BANDFragmentJob|, by subclassing |ADFFragmentJob|, which in turn is a subclass of |MultiJob|. +The constructor (``__init__``) of this new job takes 2 more arguments (``fragment1`` and ``fragment2``) and one optional argument ``full_settings`` for additional input keywords that are used **only** in the full system calculation. Furthermore, you can specify +the optimized fragment geometries using ``fragment1_opt`` and ``fragment2_opt`` for single-point calculations to also obtain the preparation energies. + +In the |prerun| method two fragment jobs and the full system job are created with the proper settings and molecules. +They are then added to the ``children`` list. + +The dedicated |Results| subclass for |BANDFragmentJob| does not provide too much additional functionality. +It simply redirects the usual |AMSResults| methods to the results of the full system calculation. The pEDA results can be obtained using the ``get_energy_decomposition`` method. It will return a dictionary with the available energy decomposition terms. + +The source code of the whole module with both abovementioned classes: + +API +~~~~~~~~~~~~~~~~~~~~~~~~~ + +.. autoclass:: BANDFragmentJob() +.. autoclass:: BANDFragmentResults() \ No newline at end of file diff --git a/doc/source/examples/examples.rst b/doc/source/examples/examples.rst index 7a378444..20ae2123 100644 --- a/doc/source/examples/examples.rst +++ b/doc/source/examples/examples.rst @@ -131,6 +131,7 @@ The source code of ``recipes`` modules is presented here to demonstrate how easy ADFCOSMORSConformers MDJobs adffragment + bandfragment ReorganizationEnergy adfnbo numgrad diff --git a/recipes/adffragment.py b/recipes/adffragment.py index e63b7b85..c17d2bc1 100644 --- a/recipes/adffragment.py +++ b/recipes/adffragment.py @@ -10,21 +10,27 @@ class ADFFragmentResults(Results): def get_properties(self): + """Redirect to |ADFResults| of the full calculation.""" return self.job.full.results.get_properties() def get_main_molecule(self): + """Redirect to |ADFResults| of the full calculation.""" return self.job.full.results.get_main_molecule() def get_input_molecule(self): + """Redirect to |ADFResults| of the full calculation.""" return self.job.full.results.get_input_molecule() def get_energy(self, unit="au"): + """Redirect to |ADFResults| of the full calculation.""" return self.job.full.results.get_energy(unit) def get_dipole_vector(self, unit="au"): + """Redirect to |ADFResults| of the full calculation.""" return self.job.full.results.get_dipole_vector(unit) def get_energy_decomposition(self): + """Get the energy decomposition of the fragment calculation.""" energy_section = self.job.full.results.read_rkf_section("Energy", file="adf") ret = {} for k in ["Electrostatic Energy", "Kinetic Energy", "Elstat Interaction", "XC Energy"]: @@ -33,15 +39,24 @@ def get_energy_decomposition(self): class ADFFragmentJob(MultiJob): + """Subclass of |MultiJob| for ADF fragment calculations.""" _result_type = ADFFragmentResults def __init__(self, fragment1=None, fragment2=None, full_settings=None, **kwargs): + """ + Args: + fragment1 (Molecule): The first fragment. + fragment2 (Molecule): The second fragment. + full_settings (Settings): The settings for the full calculation. + **kwargs: Further keyword arguments for |MultiJob|. + """ MultiJob.__init__(self, **kwargs) self.fragment1 = fragment1.copy() if isinstance(fragment1, Molecule) else fragment1 self.fragment2 = fragment2.copy() if isinstance(fragment2, Molecule) else fragment2 self.full_settings = full_settings or Settings() def prerun(self): # noqa F811 + """Prepare the fragments and the full calculation.""" self.f1 = AMSJob(name="frag1", molecule=self.fragment1, settings=self.settings) self.f2 = AMSJob(name="frag2", molecule=self.fragment2, settings=self.settings) From 9024089aef4adfbbefd37db86bb8ec88f2641a67 Mon Sep 17 00:00:00 2001 From: Patrick Melix Date: Thu, 8 Aug 2024 10:51:05 +0200 Subject: [PATCH 05/29] code example --- doc/source/examples/bandfragment.rst | 56 +++++++++++++++++++++++++++- 1 file changed, 55 insertions(+), 1 deletion(-) diff --git a/doc/source/examples/bandfragment.rst b/doc/source/examples/bandfragment.rst index 9e352ef5..2f40f7ac 100644 --- a/doc/source/examples/bandfragment.rst +++ b/doc/source/examples/bandfragment.rst @@ -24,4 +24,58 @@ API ~~~~~~~~~~~~~~~~~~~~~~~~~ .. autoclass:: BANDFragmentJob() -.. autoclass:: BANDFragmentResults() \ No newline at end of file +.. autoclass:: BANDFragmentResults() + + +Example +~~~~~~~~~~~~~~~~~~~~~~~~~ + +A basic example using `ASE to build a surface slab `_ and perform a BAND fragment calculation: + +.. code-block:: python + + #build the surface + from ase.build import fcc111, add_adsorbate + mol = fcc111('Au', size=(2,2,3)) + add_adsorbate(mol, 'H', 1.5, 'ontop') + mol.center(vacuum=10.0, axis=2) + + # separate the fragments + surface = mol.copy() + symbols = surface.get_chemical_symbols() + del surface[[i for i in range(len(symbols)) if symbols[i] != 'Au']] + adsorbate = mol.copy() + del adsorbate[[i for i in range(len(symbols)) if symbols[i] == 'Au']] + + # optional: load the optimized molecules + from ase import io + surface_opt = io.read('surface_opt.xyz') + adsorbate_opt = io.read('adsorbate_opt.xyz') + assert len(surface_opt) == len(surface) + assert len(adsorbate_opt) == len(adsorbate) + + # settings for job + base_settings = Settings() + base_settings.input.ams.task = 'SinglePoint' + base_settings.input.band.basis.type = 'TZP' + base_settings.input.band.basis.core = 'Medium' + base_settings.input.band.dos.calcdos = 'No' + base_settings.input.band.kspace.regular.numberofpoints = '5 5 1' + base_settings.input.band.beckegrid.quality = 'Good' + base_settings.input.band.zlmfit.quality = 'Good' + base_settings.input.band.usesymmetry = False + base_settings.input.band.xc.gga = 'PBE' + base_settings.input.band.xc.dispersion = 'Grimme4' + + eda_settings = plams.Settings() + eda_settings.input.band.peda = '' + + eda_job = BANDFragmentJob(fragment1=fromASE(surface), fragment2=fromASE(adsorbate), + settings=base_settings, full_settings=eda_settings, + fragment1_opt=fromASE(surface_opt), fragment2_opt=fromASE(adsorbate_opt)) + + results = eda_job.run() + eda_res = job.results.get_energy_decomposition() + print("{:<20} {:>10}".format('Term', 'Energy [kJ/mol]')) + for key, value in eda_res.items(): + print("{:<20} {:>10.4f}".format(key, value)) From 4c9f248b67a5abeafdc4d8ff98cbecf9678340ad Mon Sep 17 00:00:00 2001 From: Patrick Melix Date: Thu, 8 Aug 2024 10:51:46 +0200 Subject: [PATCH 06/29] fix error in results, add loading of optimized sp, add loading of NOCV --- recipes/bandfragment.py | 40 ++++++++++++++++++++++++++++++++++------ 1 file changed, 34 insertions(+), 6 deletions(-) diff --git a/recipes/bandfragment.py b/recipes/bandfragment.py index e73f0b62..fbe26afd 100644 --- a/recipes/bandfragment.py +++ b/recipes/bandfragment.py @@ -8,7 +8,7 @@ from os.path import join as opj -from os.path import relpath, abspath, isdir, basename +from os.path import abspath, isdir, basename from typing import Dict, List, Union __all__ = ['BANDFragmentJob', 'BANDFragmentResults'] @@ -27,8 +27,8 @@ def get_energy_decomposition(self, unit='kJ/mol') -> Dict[str, float]: Dict[str, float]: The energy decomposition. """ res=self.job.full.results - res1=self.job.f2.results - res2=self.job.f1.results + res1=self.job.f1.results + res2=self.job.f2.results ret = {} pos = 2 # E_int appears in a comment below the PEDA Table of the output @@ -147,6 +147,14 @@ def load_external(cls, path: str, jobname: str=None) -> "BANDFragmentJob": job.f1 = AMSJob.load_external(opj(path, 'frag1')) job.f2 = AMSJob.load_external(opj(path, 'frag2')) job.full = AMSJob.load_external(opj(path, 'full')) + job.children = [job.f1, job.f2, job.full] + + if isdir(opj(path, 'frag1_opt')): + job.f1_opt = AMSJob.load_external(opj(path, 'frag1_opt')) + job.children.append(job.f1_opt) + if isdir(opj(path, 'frag2_opt')): + job.f2_opt = AMSJob.load_external(opj(path, 'frag2_opt')) + job.children.append(job.f2_opt) return job @@ -177,8 +185,10 @@ def new_children(self) -> Union[None, List[AMSJob]]: return None else: # add NOCV run - self.nocv=AMSJob(name="NOCV", molecule= self.fragment1 + self.fragment2, - settings = self.settings + self.full_settings+self.nocv_settings) + set = self.settings + self.full_settings + self.nocv_settings + self.nocv = AMSJob(name="NOCV", + molecule= self.fragment1 + self.fragment2, + settings = set) # make sure we at least have the restart section to put the file name if not hasattr(self.nocv.settings, 'input'): self.nocv.settings.input = Settings() @@ -188,4 +198,22 @@ def new_children(self) -> Union[None, List[AMSJob]]: self.nocv.settings.input.band.restart = Settings() if not hasattr(self.nocv.settings.input.band.restart, 'file'): self.nocv.settings.input.band.Restart.File = opj("../", self.full.name,'band.rkf') - return [self.nocv] \ No newline at end of file + return [self.nocv] + + + @classmethod + def load_external(cls, path: str, jobname: str=None) -> "NOCVBandFragmentJob": + """Load the results of the BANDFragmentJob job from an external path. + + Args: + path (str): The path to the job. It should at least have the + subfolders 'frag1', 'frag2', 'full' and 'NOCV'. + jobname (str, optional): The name of the job. Defaults to None. + + Returns: + NOCVBandFragmentJob: The job with the loaded results. + """ + job = super().load_external(path, jobname) + job.nocv = AMSJob.load_external(opj(path, 'NOCV')) + job.children.append(job.nocv) + return job \ No newline at end of file From e05a3652311d5541e9ea09a2ee465c90df397798 Mon Sep 17 00:00:00 2001 From: Patrick Melix Date: Wed, 14 Aug 2024 11:06:14 +0200 Subject: [PATCH 07/29] fix relative restart paths not working by using symlinks --- recipes/bandfragment.py | 25 +++++++++++++++++++++---- 1 file changed, 21 insertions(+), 4 deletions(-) diff --git a/recipes/bandfragment.py b/recipes/bandfragment.py index fbe26afd..d0a188c5 100644 --- a/recipes/bandfragment.py +++ b/recipes/bandfragment.py @@ -5,10 +5,12 @@ from scm.plams.tools.units import Units from scm.plams.core.errors import FileError from scm.plams.mol.molecule import Molecule +from scm.plams import add_to_instance from os.path import join as opj -from os.path import abspath, isdir, basename +from os.path import abspath, isdir, basename, relpath +from os import symlink from typing import Dict, List, Union __all__ = ['BANDFragmentJob', 'BANDFragmentResults'] @@ -94,9 +96,11 @@ def create_mapping_setting(self) -> None: set1.atommapping = { str(i+1): str(i+1) for i in range(len(self.fragment1)) } set2 = Settings() set2.atommapping = { str(i+1): str(i+1+len(self.fragment1)) for i in range(len(self.fragment2))} - # get the correct restart files, working with relative paths - set1.filename = opj("../", self.f1.name, "band.rkf") - set2.filename = opj("../", self.f2.name, "band.rkf") + # get the correct restart files + # working with relative paths does not work for unknown reasons + # using symlinks instead + set1.filename = f"{self.f1.name}.rkf" + set2.filename = f"{self.f2.name}.rkf" self.full_settings.input.band.fragment = [set1, set2] @@ -113,6 +117,19 @@ def new_children(self) -> Union[None, List[AMSJob]]: settings = self.settings + self.full_settings) # dependencies are optional, but let's set them up self.full.depend += [self.f1,self.f2] + # save the fragment paths for the prerun of the full job + self.full.frag_paths = [] + for job in [self.f1, self.f2]: + self.full.frag_paths.append(job.path) + # edit full prerun to create symlinks + @add_to_instance(self.full) + def prerun(self): + """Create symlinks for the restart files.""" + for i, job in enumerate(['frag1', 'frag2']): + rel_path = relpath(self.frag_paths[i], self.path) + symlink(opj(rel_path, 'band.rkf'), opj(self.path, f"{job}.rkf")) + + return [self.full] From b99ff4a8b945aca9df5728066ca9281921abc175 Mon Sep 17 00:00:00 2001 From: Patrick Melix Date: Mon, 19 Aug 2024 15:18:42 +0200 Subject: [PATCH 08/29] fix non-working adffragment --- recipes/adffragment.py | 39 ++++++++++++++++++++++++++++++++------- 1 file changed, 32 insertions(+), 7 deletions(-) diff --git a/recipes/adffragment.py b/recipes/adffragment.py index c17d2bc1..c9b32504 100644 --- a/recipes/adffragment.py +++ b/recipes/adffragment.py @@ -3,6 +3,12 @@ from scm.plams.core.settings import Settings from scm.plams.interfaces.adfsuite.ams import AMSJob from scm.plams.mol.molecule import Molecule +from scm.plams import add_to_instance + +from os.path import join as opj +from os.path import relpath +from os import symlink +from typing import List, Union __all__ = ["ADFFragmentJob", "ADFFragmentResults"] @@ -65,11 +71,30 @@ def prerun(self): # noqa F811 for at in self.fragment2: at.properties.suffix = "adf.f=subsystem2" - self.full = AMSJob( - name="full", molecule=self.fragment1 + self.fragment2, settings=self.settings + self.full_settings - ) - - self.full.settings.input.adf.fragments.subsystem1 = (self.f1, "adf") - self.full.settings.input.adf.fragments.subsystem2 = (self.f2, "adf") + self.children = [self.f1, self.f2] + + def new_children(self) -> Union[None, List[AMSJob]]: + """After the first round, add the full job to the children list.""" + if hasattr(self, 'full'): + return None + else: + settings = self.settings + self.full_settings + settings.input.adf.fragments.subsystem1 = f"{self.f1.name}.rkf" + settings.input.adf.fragments.subsystem2 = f"{self.f2.name}.rkf" + self.full = AMSJob(name="full", molecule=self.fragment1 + self.fragment2, + settings=settings) + self.full.depend += [self.f1,self.f2] + # save the fragment paths for the prerun of the full job + self.full.frag_paths = [] + for job in [self.f1, self.f2]: + self.full.frag_paths.append(job.path) + # edit full prerun to create symlinks + @add_to_instance(self.full) + def prerun(self): + """Create symlinks for the restart files.""" + for i, job in enumerate(['frag1', 'frag2']): + rel_path = relpath(self.frag_paths[i], self.path) + symlink(opj(rel_path, 'adf.rkf'), opj(self.path, f"{job}.rkf")) + + return [self.full] - self.children = [self.f1, self.f2, self.full] From 609102213dc5d17ebee5d42ee25db45b75ddb436 Mon Sep 17 00:00:00 2001 From: Patrick Melix Date: Mon, 26 Aug 2024 10:24:23 +0200 Subject: [PATCH 09/29] Update results collector of ADFFragmentResults, now get the full eda results --- recipes/adffragment.py | 74 +++++++++++++++++++++++++++++++++++++++--- 1 file changed, 69 insertions(+), 5 deletions(-) diff --git a/recipes/adffragment.py b/recipes/adffragment.py index c9b32504..fa08835d 100644 --- a/recipes/adffragment.py +++ b/recipes/adffragment.py @@ -2,13 +2,16 @@ from scm.plams.core.results import Results from scm.plams.core.settings import Settings from scm.plams.interfaces.adfsuite.ams import AMSJob +from scm.plams.tools.units import Units +from scm.plams.core.errors import FileError from scm.plams.mol.molecule import Molecule from scm.plams import add_to_instance + from os.path import join as opj -from os.path import relpath +from os.path import abspath, isdir, basename, relpath from os import symlink -from typing import List, Union +from typing import Dict, List, Union __all__ = ["ADFFragmentJob", "ADFFragmentResults"] @@ -35,12 +38,46 @@ def get_dipole_vector(self, unit="au"): """Redirect to |ADFResults| of the full calculation.""" return self.job.full.results.get_dipole_vector(unit) - def get_energy_decomposition(self): - """Get the energy decomposition of the fragment calculation.""" + def get_energy_decomposition(self, unit='kJ/mol') -> Dict[str, float]: + """Get the energy decomposition of the fragment calculation. + + Args: + unit (str, optional): The unit of the energy. Defaults to 'kJ/mol'. + + Returns: + Dict[str, float]: The energy decomposition. + """ energy_section = self.job.full.results.read_rkf_section("Energy", file="adf") ret = {} for k in ["Electrostatic Energy", "Kinetic Energy", "Elstat Interaction", "XC Energy"]: - ret[k] = energy_section[k] + ret[k] = Units.convert(energy_section[k], 'au', unit) + + # most information not available from the KF file + res=self.job.full.results + res1=self.job.f1.results + res2=self.job.f2.results + pos = -4 # position of the energy in au the output + # E_int appears in a comment below the PEDA Table of the output + ret['E_int'] = Units.convert(float(res.grep_output('Total Bonding Energy:')[-2].split()[pos]), 'au', unit) + ret['E_int_disp'] = Units.convert(float(res.grep_output('Dispersion Energy:')[-1].split()[pos]), 'au', unit) + ret['E_Pauli'] = Units.convert(float(res.grep_output('Pauli Repulsion (Delta')[-1].split()[pos]), 'au', unit) + ret['E_elstat'] = Units.convert(float(res.grep_output('Electrostatic Interaction:')[-1].split()[pos]), 'au', unit) + ret['E_orb'] = Units.convert(float(res.grep_output('Total Orbital Interactions:')[-1].split()[pos]), 'au', unit) + + ret['E_1'] = res1.get_energy(unit=unit) + ret['E_2'] = res2.get_energy(unit=unit) + + if hasattr(self.job, 'f1_opt'): + ret['E_1_opt'] = self.job.f1_opt.results.get_energy(unit=unit) + ret['E_prep_1'] = ret['E_1'] - ret['E_1_opt'] + if hasattr(self.job, 'f2_opt'): + ret['E_2_opt'] = self.job.f2_opt.results.get_energy(unit=unit) + ret['E_prep_2'] = ret['E_2'] - ret['E_2_opt'] + + if ('E_1_opt' in ret) and ('E_2_opt' in ret): + ret['E_prep'] = ret['E_prep_1'] + ret['E_prep_2'] + ret['E_bond'] = ret['E_int'] + ret['E_prep'] + return ret @@ -98,3 +135,30 @@ def prerun(self): return [self.full] + @classmethod + def load_external(cls, path: str, jobname: str=None) -> "ADFFragmentJob": + """Load the results of the ADFFragmentJob job from an external path. + + Args: + path (str): The path to the job. It should at least have the + subfolders 'frag1', 'frag2' and 'full'. + jobname (str, optional): The name of the job. Defaults to None. + + Returns: + ADFFragmentJob: The job with the loaded results. + """ + if not isdir(path): + raise FileError("Path {} does not exist, cannot load from it.".format(path)) + path = abspath(path) + jobname = basename(path) if jobname is None else str(jobname) + + job = cls(name=jobname) + job.path = path + job.status = "copied" + + job.f1 = AMSJob.load_external(opj(path, 'frag1')) + job.f2 = AMSJob.load_external(opj(path, 'frag2')) + job.full = AMSJob.load_external(opj(path, 'full')) + job.children = [job.f1, job.f2, job.full] + + return job \ No newline at end of file From 6f7198d68c0ccbd291991921f09e6036b6bf3bba Mon Sep 17 00:00:00 2001 From: Patrick Melix Date: Tue, 27 Aug 2024 11:25:44 +0200 Subject: [PATCH 10/29] some work on NOCV --- recipes/bandfragment.py | 30 +++++++++++++++++++----------- 1 file changed, 19 insertions(+), 11 deletions(-) diff --git a/recipes/bandfragment.py b/recipes/bandfragment.py index d0a188c5..345fb18d 100644 --- a/recipes/bandfragment.py +++ b/recipes/bandfragment.py @@ -181,7 +181,7 @@ class NOCVBandFragmentJob(BANDFragmentJob): _result_type = BANDFragmentResults def __init__(self, nocv_settings=None, **kwargs) -> None: - """Create a BANDFragmentJob with NOCV calculation. + """Create a BANDFragmentJob with final NOCV calculation. Args: nocv_settings (Settings, optional): Settings for the NOCV calculation. Defaults to None. @@ -195,26 +195,34 @@ def new_children(self) -> Union[None, List[AMSJob]]: """After the first round, add the full job to the children list. After the second round, add the NOCV job to the children list.""" # new_children of BANDFragmentJob creates the full job - ret = super().prerun() + ret = super().new_children() if ret is not None: return ret if hasattr(self, 'nocv'): return None else: # add NOCV run + # settings for the NOCV calculation set = self.settings + self.full_settings + self.nocv_settings + # set the restart file + set.input.band.restart.file = "full.rkf" + # copy the fragment settings + set.input.band.fragment = [ fset.copy() for fset in self.full.settings.input.band.fragment ] + self.nocv = AMSJob(name="NOCV", molecule= self.fragment1 + self.fragment2, settings = set) - # make sure we at least have the restart section to put the file name - if not hasattr(self.nocv.settings, 'input'): - self.nocv.settings.input = Settings() - if not hasattr(self.nocv.settings.input, 'band'): - self.nocv.settings.input.band = Settings() - if not hasattr(self.nocv.settings.input.band, 'restart'): - self.nocv.settings.input.band.restart = Settings() - if not hasattr(self.nocv.settings.input.band.restart, 'file'): - self.nocv.settings.input.band.Restart.File = opj("../", self.full.name,'band.rkf') + + self.nocv.frag_paths = [] + for job in [self.f1, self.f2, self.full]: + self.nocv.frag_paths.append(job.path) + # edit NOCV prerun to create symlinks + @add_to_instance(self.nocv) + def prerun(self): + """Create symlinks for the restart files.""" + for i, job in enumerate(['frag1', 'frag2', 'full']): + rel_path = relpath(self.frag_paths[i], self.path) + symlink(opj(rel_path, 'band.rkf'), opj(self.path, f"{job}.rkf")) return [self.nocv] From a8ce54bd70dd4125e3075433e99a744197546962 Mon Sep 17 00:00:00 2001 From: Patrick Melix Date: Tue, 27 Aug 2024 11:26:31 +0200 Subject: [PATCH 11/29] fix add_to_instance import as suggested --- recipes/adffragment.py | 2 +- recipes/bandfragment.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/recipes/adffragment.py b/recipes/adffragment.py index fa08835d..be7a1431 100644 --- a/recipes/adffragment.py +++ b/recipes/adffragment.py @@ -5,7 +5,7 @@ from scm.plams.tools.units import Units from scm.plams.core.errors import FileError from scm.plams.mol.molecule import Molecule -from scm.plams import add_to_instance +from scm.plams.core.functions import add_to_instance from os.path import join as opj diff --git a/recipes/bandfragment.py b/recipes/bandfragment.py index 345fb18d..a78786f1 100644 --- a/recipes/bandfragment.py +++ b/recipes/bandfragment.py @@ -5,7 +5,7 @@ from scm.plams.tools.units import Units from scm.plams.core.errors import FileError from scm.plams.mol.molecule import Molecule -from scm.plams import add_to_instance +from scm.plams.core.functions import add_to_instance from os.path import join as opj From f5893bea9b0842b60194eec586c0adf780c052e9 Mon Sep 17 00:00:00 2001 From: Patrick Melix Date: Tue, 27 Aug 2024 11:30:18 +0200 Subject: [PATCH 12/29] formatting with black --- recipes/adffragment.py | 83 +++++++++++----------- recipes/bandfragment.py | 150 +++++++++++++++++++--------------------- 2 files changed, 115 insertions(+), 118 deletions(-) diff --git a/recipes/adffragment.py b/recipes/adffragment.py index be7a1431..9644f004 100644 --- a/recipes/adffragment.py +++ b/recipes/adffragment.py @@ -38,51 +38,54 @@ def get_dipole_vector(self, unit="au"): """Redirect to |ADFResults| of the full calculation.""" return self.job.full.results.get_dipole_vector(unit) - def get_energy_decomposition(self, unit='kJ/mol') -> Dict[str, float]: + def get_energy_decomposition(self, unit="kJ/mol") -> Dict[str, float]: """Get the energy decomposition of the fragment calculation. - + Args: unit (str, optional): The unit of the energy. Defaults to 'kJ/mol'. - + Returns: Dict[str, float]: The energy decomposition. """ energy_section = self.job.full.results.read_rkf_section("Energy", file="adf") ret = {} for k in ["Electrostatic Energy", "Kinetic Energy", "Elstat Interaction", "XC Energy"]: - ret[k] = Units.convert(energy_section[k], 'au', unit) - + ret[k] = Units.convert(energy_section[k], "au", unit) + # most information not available from the KF file - res=self.job.full.results - res1=self.job.f1.results - res2=self.job.f2.results + res = self.job.full.results + res1 = self.job.f1.results + res2 = self.job.f2.results pos = -4 # position of the energy in au the output # E_int appears in a comment below the PEDA Table of the output - ret['E_int'] = Units.convert(float(res.grep_output('Total Bonding Energy:')[-2].split()[pos]), 'au', unit) - ret['E_int_disp'] = Units.convert(float(res.grep_output('Dispersion Energy:')[-1].split()[pos]), 'au', unit) - ret['E_Pauli'] = Units.convert(float(res.grep_output('Pauli Repulsion (Delta')[-1].split()[pos]), 'au', unit) - ret['E_elstat'] = Units.convert(float(res.grep_output('Electrostatic Interaction:')[-1].split()[pos]), 'au', unit) - ret['E_orb'] = Units.convert(float(res.grep_output('Total Orbital Interactions:')[-1].split()[pos]), 'au', unit) - - ret['E_1'] = res1.get_energy(unit=unit) - ret['E_2'] = res2.get_energy(unit=unit) - - if hasattr(self.job, 'f1_opt'): - ret['E_1_opt'] = self.job.f1_opt.results.get_energy(unit=unit) - ret['E_prep_1'] = ret['E_1'] - ret['E_1_opt'] - if hasattr(self.job, 'f2_opt'): - ret['E_2_opt'] = self.job.f2_opt.results.get_energy(unit=unit) - ret['E_prep_2'] = ret['E_2'] - ret['E_2_opt'] - - if ('E_1_opt' in ret) and ('E_2_opt' in ret): - ret['E_prep'] = ret['E_prep_1'] + ret['E_prep_2'] - ret['E_bond'] = ret['E_int'] + ret['E_prep'] - + ret["E_int"] = Units.convert(float(res.grep_output("Total Bonding Energy:")[-2].split()[pos]), "au", unit) + ret["E_int_disp"] = Units.convert(float(res.grep_output("Dispersion Energy:")[-1].split()[pos]), "au", unit) + ret["E_Pauli"] = Units.convert(float(res.grep_output("Pauli Repulsion (Delta")[-1].split()[pos]), "au", unit) + ret["E_elstat"] = Units.convert( + float(res.grep_output("Electrostatic Interaction:")[-1].split()[pos]), "au", unit + ) + ret["E_orb"] = Units.convert(float(res.grep_output("Total Orbital Interactions:")[-1].split()[pos]), "au", unit) + + ret["E_1"] = res1.get_energy(unit=unit) + ret["E_2"] = res2.get_energy(unit=unit) + + if hasattr(self.job, "f1_opt"): + ret["E_1_opt"] = self.job.f1_opt.results.get_energy(unit=unit) + ret["E_prep_1"] = ret["E_1"] - ret["E_1_opt"] + if hasattr(self.job, "f2_opt"): + ret["E_2_opt"] = self.job.f2_opt.results.get_energy(unit=unit) + ret["E_prep_2"] = ret["E_2"] - ret["E_2_opt"] + + if ("E_1_opt" in ret) and ("E_2_opt" in ret): + ret["E_prep"] = ret["E_prep_1"] + ret["E_prep_2"] + ret["E_bond"] = ret["E_int"] + ret["E_prep"] + return ret class ADFFragmentJob(MultiJob): """Subclass of |MultiJob| for ADF fragment calculations.""" + _result_type = ADFFragmentResults def __init__(self, fragment1=None, fragment2=None, full_settings=None, **kwargs): @@ -112,33 +115,33 @@ def prerun(self): # noqa F811 def new_children(self) -> Union[None, List[AMSJob]]: """After the first round, add the full job to the children list.""" - if hasattr(self, 'full'): + if hasattr(self, "full"): return None else: settings = self.settings + self.full_settings settings.input.adf.fragments.subsystem1 = f"{self.f1.name}.rkf" settings.input.adf.fragments.subsystem2 = f"{self.f2.name}.rkf" - self.full = AMSJob(name="full", molecule=self.fragment1 + self.fragment2, - settings=settings) - self.full.depend += [self.f1,self.f2] + self.full = AMSJob(name="full", molecule=self.fragment1 + self.fragment2, settings=settings) + self.full.depend += [self.f1, self.f2] # save the fragment paths for the prerun of the full job self.full.frag_paths = [] for job in [self.f1, self.f2]: self.full.frag_paths.append(job.path) + # edit full prerun to create symlinks @add_to_instance(self.full) def prerun(self): """Create symlinks for the restart files.""" - for i, job in enumerate(['frag1', 'frag2']): + for i, job in enumerate(["frag1", "frag2"]): rel_path = relpath(self.frag_paths[i], self.path) - symlink(opj(rel_path, 'adf.rkf'), opj(self.path, f"{job}.rkf")) + symlink(opj(rel_path, "adf.rkf"), opj(self.path, f"{job}.rkf")) return [self.full] @classmethod - def load_external(cls, path: str, jobname: str=None) -> "ADFFragmentJob": + def load_external(cls, path: str, jobname: str = None) -> "ADFFragmentJob": """Load the results of the ADFFragmentJob job from an external path. - + Args: path (str): The path to the job. It should at least have the subfolders 'frag1', 'frag2' and 'full'. @@ -156,9 +159,9 @@ def load_external(cls, path: str, jobname: str=None) -> "ADFFragmentJob": job.path = path job.status = "copied" - job.f1 = AMSJob.load_external(opj(path, 'frag1')) - job.f2 = AMSJob.load_external(opj(path, 'frag2')) - job.full = AMSJob.load_external(opj(path, 'full')) + job.f1 = AMSJob.load_external(opj(path, "frag1")) + job.f2 = AMSJob.load_external(opj(path, "frag2")) + job.full = AMSJob.load_external(opj(path, "full")) job.children = [job.f1, job.f2, job.full] - return job \ No newline at end of file + return job diff --git a/recipes/bandfragment.py b/recipes/bandfragment.py index a78786f1..3aabafca 100644 --- a/recipes/bandfragment.py +++ b/recipes/bandfragment.py @@ -13,57 +13,56 @@ from os import symlink from typing import Dict, List, Union -__all__ = ['BANDFragmentJob', 'BANDFragmentResults'] +__all__ = ["BANDFragmentJob", "BANDFragmentResults"] class BANDFragmentResults(ADFFragmentResults): """Subclass of |ADFFragmentResults| for BAND calculations.""" - def get_energy_decomposition(self, unit='kJ/mol') -> Dict[str, float]: + def get_energy_decomposition(self, unit="kJ/mol") -> Dict[str, float]: """Get the energy decomposition of the fragment calculation. - + Args: unit (str, optional): The unit of the energy. Defaults to 'kJ/mol'. - + Returns: Dict[str, float]: The energy decomposition. """ - res=self.job.full.results - res1=self.job.f1.results - res2=self.job.f2.results + res = self.job.full.results + res1 = self.job.f1.results + res2 = self.job.f2.results ret = {} pos = 2 # E_int appears in a comment below the PEDA Table of the output - ret['E_int'] = Units.convert(float(res.grep_output('E_int')[-2].split()[pos]), 'au', unit) - ret['E_int_disp'] = Units.convert(float(res.grep_output('E_disp')[-1].split()[pos]), 'au', unit) - ret['E_Pauli'] = Units.convert(float(res.grep_output('E_Pauli')[-1].split()[pos]), 'au', unit) - ret['E_elstat'] = Units.convert(float(res.grep_output('E_elstat')[-1].split()[pos]), 'au', unit) - ret['E_orb'] = Units.convert(float(res.grep_output('E_orb')[-1].split()[pos]), 'au', unit) - - ret['E_1'] = res1.get_energy(unit=unit) - ret['E_2'] = res2.get_energy(unit=unit) - - if hasattr(self.job, 'f1_opt'): - ret['E_1_opt'] = self.job.f1_opt.results.get_energy(unit=unit) - ret['E_prep_1'] = ret['E_1'] - ret['E_1_opt'] - if hasattr(self.job, 'f2_opt'): - ret['E_2_opt'] = self.job.f2_opt.results.get_energy(unit=unit) - ret['E_prep_2'] = ret['E_2'] - ret['E_2_opt'] - - if ('E_1_opt' in ret) and ('E_2_opt' in ret): - ret['E_prep'] = ret['E_prep_1'] + ret['E_prep_2'] - ret['E_bond'] = ret['E_int'] + ret['E_prep'] - + ret["E_int"] = Units.convert(float(res.grep_output("E_int")[-2].split()[pos]), "au", unit) + ret["E_int_disp"] = Units.convert(float(res.grep_output("E_disp")[-1].split()[pos]), "au", unit) + ret["E_Pauli"] = Units.convert(float(res.grep_output("E_Pauli")[-1].split()[pos]), "au", unit) + ret["E_elstat"] = Units.convert(float(res.grep_output("E_elstat")[-1].split()[pos]), "au", unit) + ret["E_orb"] = Units.convert(float(res.grep_output("E_orb")[-1].split()[pos]), "au", unit) + + ret["E_1"] = res1.get_energy(unit=unit) + ret["E_2"] = res2.get_energy(unit=unit) + + if hasattr(self.job, "f1_opt"): + ret["E_1_opt"] = self.job.f1_opt.results.get_energy(unit=unit) + ret["E_prep_1"] = ret["E_1"] - ret["E_1_opt"] + if hasattr(self.job, "f2_opt"): + ret["E_2_opt"] = self.job.f2_opt.results.get_energy(unit=unit) + ret["E_prep_2"] = ret["E_2"] - ret["E_2_opt"] + + if ("E_1_opt" in ret) and ("E_2_opt" in ret): + ret["E_prep"] = ret["E_prep_1"] + ret["E_prep_2"] + ret["E_bond"] = ret["E_int"] + ret["E_prep"] + return ret class BANDFragmentJob(ADFFragmentJob): _result_type = BANDFragmentResults - def __init__(self, fragment1_opt: Molecule=None, - fragment2_opt: Molecule=None, **kwargs) -> None: + def __init__(self, fragment1_opt: Molecule = None, fragment2_opt: Molecule = None, **kwargs) -> None: """Create a BANDFragmentJob with the given fragments. - + The optimized fragment structures can be given as arguments. Single point calculations on these structures will then be performed to obtain the preparation energy. @@ -78,24 +77,25 @@ def __init__(self, fragment1_opt: Molecule=None, super().__init__(**kwargs) if isinstance(fragment1_opt, Molecule): self.fragment1_opt = fragment1_opt.copy() - self.f1_opt = AMSJob(name='frag1_opt', molecule=self.fragment1_opt, - settings=self.settings) + self.f1_opt = AMSJob(name="frag1_opt", molecule=self.fragment1_opt, settings=self.settings) self.children.append(self.f1_opt) if isinstance(fragment2_opt, Molecule): self.fragment2_opt = fragment2_opt.copy() - self.f2_opt = AMSJob(name='frag2_opt', molecule=self.fragment2_opt, - settings=self.settings) + self.f2_opt = AMSJob(name="frag2_opt", molecule=self.fragment2_opt, settings=self.settings) self.children.append(self.f2_opt) def create_mapping_setting(self) -> None: if "fragment" in self.full_settings.input.band: - log("Fragment already present in full_settings. Assuming that the user has already set up the mapping. Skipping the mapping setup.", level=1) + log( + "Fragment already present in full_settings. Assuming that the user has already set up the mapping. Skipping the mapping setup.", + level=1, + ) return # first fragment 1 then fragment 2 set1 = Settings() - set1.atommapping = { str(i+1): str(i+1) for i in range(len(self.fragment1)) } + set1.atommapping = {str(i + 1): str(i + 1) for i in range(len(self.fragment1))} set2 = Settings() - set2.atommapping = { str(i+1): str(i+1+len(self.fragment1)) for i in range(len(self.fragment2))} + set2.atommapping = {str(i + 1): str(i + 1 + len(self.fragment1)) for i in range(len(self.fragment2))} # get the correct restart files # working with relative paths does not work for unknown reasons # using symlinks instead @@ -103,47 +103,44 @@ def create_mapping_setting(self) -> None: set2.filename = f"{self.f2.name}.rkf" self.full_settings.input.band.fragment = [set1, set2] - def new_children(self) -> Union[None, List[AMSJob]]: """After the first round, add the full job to the children list.""" - if hasattr(self, 'full'): + if hasattr(self, "full"): return None else: # create the correct mapping settings for the full job self.create_mapping_setting() # create the full job - self.full = AMSJob(name = 'full', - molecule = self.fragment1 + self.fragment2, - settings = self.settings + self.full_settings) + self.full = AMSJob( + name="full", molecule=self.fragment1 + self.fragment2, settings=self.settings + self.full_settings + ) # dependencies are optional, but let's set them up - self.full.depend += [self.f1,self.f2] + self.full.depend += [self.f1, self.f2] # save the fragment paths for the prerun of the full job self.full.frag_paths = [] for job in [self.f1, self.f2]: self.full.frag_paths.append(job.path) + # edit full prerun to create symlinks @add_to_instance(self.full) def prerun(self): """Create symlinks for the restart files.""" - for i, job in enumerate(['frag1', 'frag2']): + for i, job in enumerate(["frag1", "frag2"]): rel_path = relpath(self.frag_paths[i], self.path) - symlink(opj(rel_path, 'band.rkf'), opj(self.path, f"{job}.rkf")) - + symlink(opj(rel_path, "band.rkf"), opj(self.path, f"{job}.rkf")) return [self.full] - def prerun(self) -> None: """Creates the fragment jobs.""" - self.f1 = AMSJob(name='frag1', molecule=self.fragment1, settings=self.settings) - self.f2 = AMSJob(name='frag2', molecule=self.fragment2, settings=self.settings) - self.children += [self.f1,self.f2] - + self.f1 = AMSJob(name="frag1", molecule=self.fragment1, settings=self.settings) + self.f2 = AMSJob(name="frag2", molecule=self.fragment2, settings=self.settings) + self.children += [self.f1, self.f2] @classmethod - def load_external(cls, path: str, jobname: str=None) -> "BANDFragmentJob": + def load_external(cls, path: str, jobname: str = None) -> "BANDFragmentJob": """Load the results of the BANDFragmentJob job from an external path. - + Args: path (str): The path to the job. It should at least have the subfolders 'frag1', 'frag2' and 'full'. @@ -161,36 +158,34 @@ def load_external(cls, path: str, jobname: str=None) -> "BANDFragmentJob": job.path = path job.status = "copied" - job.f1 = AMSJob.load_external(opj(path, 'frag1')) - job.f2 = AMSJob.load_external(opj(path, 'frag2')) - job.full = AMSJob.load_external(opj(path, 'full')) + job.f1 = AMSJob.load_external(opj(path, "frag1")) + job.f2 = AMSJob.load_external(opj(path, "frag2")) + job.full = AMSJob.load_external(opj(path, "full")) job.children = [job.f1, job.f2, job.full] - if isdir(opj(path, 'frag1_opt')): - job.f1_opt = AMSJob.load_external(opj(path, 'frag1_opt')) + if isdir(opj(path, "frag1_opt")): + job.f1_opt = AMSJob.load_external(opj(path, "frag1_opt")) job.children.append(job.f1_opt) - if isdir(opj(path, 'frag2_opt')): - job.f2_opt = AMSJob.load_external(opj(path, 'frag2_opt')) + if isdir(opj(path, "frag2_opt")): + job.f2_opt = AMSJob.load_external(opj(path, "frag2_opt")) job.children.append(job.f2_opt) return job - class NOCVBandFragmentJob(BANDFragmentJob): _result_type = BANDFragmentResults def __init__(self, nocv_settings=None, **kwargs) -> None: """Create a BANDFragmentJob with final NOCV calculation. - + Args: nocv_settings (Settings, optional): Settings for the NOCV calculation. Defaults to None. - + """ super().__init__(**kwargs) self.nocv_settings = nocv_settings or Settings() - def new_children(self) -> Union[None, List[AMSJob]]: """After the first round, add the full job to the children list. After the second round, add the NOCV job to the children list.""" @@ -198,38 +193,37 @@ def new_children(self) -> Union[None, List[AMSJob]]: ret = super().new_children() if ret is not None: return ret - if hasattr(self, 'nocv'): + if hasattr(self, "nocv"): return None else: # add NOCV run - # settings for the NOCV calculation + # settings for the NOCV calculation set = self.settings + self.full_settings + self.nocv_settings # set the restart file set.input.band.restart.file = "full.rkf" # copy the fragment settings - set.input.band.fragment = [ fset.copy() for fset in self.full.settings.input.band.fragment ] + set.input.band.fragment = [fset.copy() for fset in self.full.settings.input.band.fragment] + + self.nocv = AMSJob(name="NOCV", molecule=self.fragment1 + self.fragment2, settings=set) - self.nocv = AMSJob(name="NOCV", - molecule= self.fragment1 + self.fragment2, - settings = set) - self.nocv.frag_paths = [] for job in [self.f1, self.f2, self.full]: self.nocv.frag_paths.append(job.path) + # edit NOCV prerun to create symlinks @add_to_instance(self.nocv) def prerun(self): """Create symlinks for the restart files.""" - for i, job in enumerate(['frag1', 'frag2', 'full']): + for i, job in enumerate(["frag1", "frag2", "full"]): rel_path = relpath(self.frag_paths[i], self.path) - symlink(opj(rel_path, 'band.rkf'), opj(self.path, f"{job}.rkf")) + symlink(opj(rel_path, "band.rkf"), opj(self.path, f"{job}.rkf")) + return [self.nocv] - @classmethod - def load_external(cls, path: str, jobname: str=None) -> "NOCVBandFragmentJob": + def load_external(cls, path: str, jobname: str = None) -> "NOCVBandFragmentJob": """Load the results of the BANDFragmentJob job from an external path. - + Args: path (str): The path to the job. It should at least have the subfolders 'frag1', 'frag2', 'full' and 'NOCV'. @@ -239,6 +233,6 @@ def load_external(cls, path: str, jobname: str=None) -> "NOCVBandFragmentJob": NOCVBandFragmentJob: The job with the loaded results. """ job = super().load_external(path, jobname) - job.nocv = AMSJob.load_external(opj(path, 'NOCV')) + job.nocv = AMSJob.load_external(opj(path, "NOCV")) job.children.append(job.nocv) - return job \ No newline at end of file + return job From 109250316e0af3cbce3dfb9a03363d84f35ba8d7 Mon Sep 17 00:00:00 2001 From: Patrick Melix Date: Tue, 27 Aug 2024 11:38:36 +0200 Subject: [PATCH 13/29] move example to file --- doc/source/examples/bandfragment.rst | 50 ++-------------------------- examples/bandfrag.py | 49 +++++++++++++++++++++++++++ 2 files changed, 51 insertions(+), 48 deletions(-) create mode 100644 examples/bandfrag.py diff --git a/doc/source/examples/bandfragment.rst b/doc/source/examples/bandfragment.rst index 2f40f7ac..7d819377 100644 --- a/doc/source/examples/bandfragment.rst +++ b/doc/source/examples/bandfragment.rst @@ -30,52 +30,6 @@ API Example ~~~~~~~~~~~~~~~~~~~~~~~~~ -A basic example using `ASE to build a surface slab `_ and perform a BAND fragment calculation: +A basic example using `ASE to build a surface slab `_ and perform a BAND fragment calculation: :download:`bandfrag_test.py <../../../../../../examples/scripting/bandfrag_test.py>` -.. code-block:: python - - #build the surface - from ase.build import fcc111, add_adsorbate - mol = fcc111('Au', size=(2,2,3)) - add_adsorbate(mol, 'H', 1.5, 'ontop') - mol.center(vacuum=10.0, axis=2) - - # separate the fragments - surface = mol.copy() - symbols = surface.get_chemical_symbols() - del surface[[i for i in range(len(symbols)) if symbols[i] != 'Au']] - adsorbate = mol.copy() - del adsorbate[[i for i in range(len(symbols)) if symbols[i] == 'Au']] - - # optional: load the optimized molecules - from ase import io - surface_opt = io.read('surface_opt.xyz') - adsorbate_opt = io.read('adsorbate_opt.xyz') - assert len(surface_opt) == len(surface) - assert len(adsorbate_opt) == len(adsorbate) - - # settings for job - base_settings = Settings() - base_settings.input.ams.task = 'SinglePoint' - base_settings.input.band.basis.type = 'TZP' - base_settings.input.band.basis.core = 'Medium' - base_settings.input.band.dos.calcdos = 'No' - base_settings.input.band.kspace.regular.numberofpoints = '5 5 1' - base_settings.input.band.beckegrid.quality = 'Good' - base_settings.input.band.zlmfit.quality = 'Good' - base_settings.input.band.usesymmetry = False - base_settings.input.band.xc.gga = 'PBE' - base_settings.input.band.xc.dispersion = 'Grimme4' - - eda_settings = plams.Settings() - eda_settings.input.band.peda = '' - - eda_job = BANDFragmentJob(fragment1=fromASE(surface), fragment2=fromASE(adsorbate), - settings=base_settings, full_settings=eda_settings, - fragment1_opt=fromASE(surface_opt), fragment2_opt=fromASE(adsorbate_opt)) - - results = eda_job.run() - eda_res = job.results.get_energy_decomposition() - print("{:<20} {:>10}".format('Term', 'Energy [kJ/mol]')) - for key, value in eda_res.items(): - print("{:<20} {:>10.4f}".format(key, value)) +.. literalinclude:: ../../../../../../examples/scripting/bandfrag_test.py \ No newline at end of file diff --git a/examples/bandfrag.py b/examples/bandfrag.py new file mode 100644 index 00000000..b0016564 --- /dev/null +++ b/examples/bandfrag.py @@ -0,0 +1,49 @@ +from scm.plams import Settings +from scm.plams import fromASE +from scm.plams.recipes.bandfragment import BANDFragmentJob + +#build the surface +from ase.build import fcc111, add_adsorbate +mol = fcc111('Au', size=(2,2,3)) +add_adsorbate(mol, 'H', 1.5, 'ontop') +mol.center(vacuum=10.0, axis=2) + +# separate the fragments +surface = mol.copy() +symbols = surface.get_chemical_symbols() +del surface[[i for i in range(len(symbols)) if symbols[i] != 'Au']] +adsorbate = mol.copy() +del adsorbate[[i for i in range(len(symbols)) if symbols[i] == 'Au']] + +# optional: load the optimized molecules +from ase import io +surface_opt = io.read('surface_opt.xyz') +adsorbate_opt = io.read('adsorbate_opt.xyz') +assert len(surface_opt) == len(surface) +assert len(adsorbate_opt) == len(adsorbate) + +# settings for job +base_settings = Settings() +base_settings.input.ams.task = 'SinglePoint' +base_settings.input.band.basis.type = 'TZP' +base_settings.input.band.basis.core = 'Medium' +base_settings.input.band.dos.calcdos = 'No' +base_settings.input.band.kspace.regular.numberofpoints = '5 5 1' +base_settings.input.band.beckegrid.quality = 'Good' +base_settings.input.band.zlmfit.quality = 'Good' +base_settings.input.band.usesymmetry = False +base_settings.input.band.xc.gga = 'PBE' +base_settings.input.band.xc.dispersion = 'Grimme4' + +eda_settings = Settings() +eda_settings.input.band.peda = '' + +eda_job = BANDFragmentJob(fragment1=fromASE(surface), fragment2=fromASE(adsorbate), + settings=base_settings, full_settings=eda_settings, + fragment1_opt=fromASE(surface_opt), fragment2_opt=fromASE(adsorbate_opt)) + +results = eda_job.run() +eda_res = eda_job.results.get_energy_decomposition() +print("{:<20} {:>10}".format('Term', 'Energy [kJ/mol]')) +for key, value in eda_res.items(): + print("{:<20} {:>10.4f}".format(key, value)) \ No newline at end of file From e6aafb6b255ac604731920baaeb7a34b00ded289 Mon Sep 17 00:00:00 2001 From: Patrick Melix Date: Tue, 27 Aug 2024 11:41:41 +0200 Subject: [PATCH 14/29] fix relative path of doc example --- doc/source/examples/bandfragment.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/source/examples/bandfragment.rst b/doc/source/examples/bandfragment.rst index 7d819377..346fa65d 100644 --- a/doc/source/examples/bandfragment.rst +++ b/doc/source/examples/bandfragment.rst @@ -30,6 +30,6 @@ API Example ~~~~~~~~~~~~~~~~~~~~~~~~~ -A basic example using `ASE to build a surface slab `_ and perform a BAND fragment calculation: :download:`bandfrag_test.py <../../../../../../examples/scripting/bandfrag_test.py>` +A basic example using `ASE to build a surface slab `_ and perform a BAND fragment calculation: :download:`bandfrag_test.py <../../../examples/scripting/bandfrag_test.py>` -.. literalinclude:: ../../../../../../examples/scripting/bandfrag_test.py \ No newline at end of file +.. literalinclude:: ../../../examples/scripting/bandfrag_test.py \ No newline at end of file From 6c9af4e15e11de2b6a60fa5712eab4e553dc0cc5 Mon Sep 17 00:00:00 2001 From: Patrick Melix Date: Tue, 27 Aug 2024 11:45:36 +0200 Subject: [PATCH 15/29] add NOCVBandFragmentJob to doc --- doc/source/conf.py | 1 + doc/source/examples/bandfragment.rst | 5 ++++- 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/doc/source/conf.py b/doc/source/conf.py index cc66167c..b3eb7b80 100644 --- a/doc/source/conf.py +++ b/doc/source/conf.py @@ -331,6 +331,7 @@ def setup(app): .. |ADFFragmentResults| replace:: :class:`~scm.plams.recipes.adffragment.ADFFragmentResults` .. |BANDFragmentJob| replace:: :class:`~scm.plams.recipes.bandfragment.BANDFragmentJob` .. |BANDFragmentResults| replace:: :class:`~scm.plams.recipes.bandfragment.BANDFragmentResults` +.. |NOCVBandFragmentJob| replace:: :class:`~scm.plams.recipes.bandfragment.NOCVBandFragmentJob` .. |RPM| replace:: :ref:`rerun-prevention` .. |cleaning| replace:: :ref:`cleaning` diff --git a/doc/source/examples/bandfragment.rst b/doc/source/examples/bandfragment.rst index 346fa65d..aff16dd2 100644 --- a/doc/source/examples/bandfragment.rst +++ b/doc/source/examples/bandfragment.rst @@ -6,7 +6,7 @@ BAND fragment job .. currentmodule:: scm.plams.recipes.bandfragment In this module a dedicated job type for Energy Decomposition Analysis in BAND is defined. -Such an analysis is performed on a periodic system divided into 2 fragments and consists of a minimum of 3 separate Band runs: one for each fragment and one for full system. See +Such an analysis is performed on a periodic system divided into 2 fragments and consists of a minimum of 3 separate BAND runs: one for each fragment and one for full system. See also |ADFFragmentJob|. We define a new job type |BANDFragmentJob|, by subclassing |ADFFragmentJob|, which in turn is a subclass of |MultiJob|. The constructor (``__init__``) of this new job takes 2 more arguments (``fragment1`` and ``fragment2``) and one optional argument ``full_settings`` for additional input keywords that are used **only** in the full system calculation. Furthermore, you can specify @@ -18,6 +18,8 @@ They are then added to the ``children`` list. The dedicated |Results| subclass for |BANDFragmentJob| does not provide too much additional functionality. It simply redirects the usual |AMSResults| methods to the results of the full system calculation. The pEDA results can be obtained using the ``get_energy_decomposition`` method. It will return a dictionary with the available energy decomposition terms. +A derived subclass |NOCVBandFragmentJob| is also provided. It can be usefull for generating NOCV plots after the PEDA-NOCV calculation. + The source code of the whole module with both abovementioned classes: API @@ -25,6 +27,7 @@ API .. autoclass:: BANDFragmentJob() .. autoclass:: BANDFragmentResults() +.. autoclass:: NOCVBandFragmentJob() Example From aa3e06605736ad804c5294aae4fb3bc2b87cea13 Mon Sep 17 00:00:00 2001 From: Patrick Melix Date: Tue, 27 Aug 2024 11:49:07 +0200 Subject: [PATCH 16/29] run black on bandfrag.py --- examples/bandfrag.py | 51 +++++++++++++++++++++++++------------------- 1 file changed, 29 insertions(+), 22 deletions(-) diff --git a/examples/bandfrag.py b/examples/bandfrag.py index b0016564..5bb93e65 100644 --- a/examples/bandfrag.py +++ b/examples/bandfrag.py @@ -2,48 +2,55 @@ from scm.plams import fromASE from scm.plams.recipes.bandfragment import BANDFragmentJob -#build the surface +# build the surface from ase.build import fcc111, add_adsorbate -mol = fcc111('Au', size=(2,2,3)) -add_adsorbate(mol, 'H', 1.5, 'ontop') + +mol = fcc111("Au", size=(2, 2, 3)) +add_adsorbate(mol, "H", 1.5, "ontop") mol.center(vacuum=10.0, axis=2) # separate the fragments surface = mol.copy() symbols = surface.get_chemical_symbols() -del surface[[i for i in range(len(symbols)) if symbols[i] != 'Au']] +del surface[[i for i in range(len(symbols)) if symbols[i] != "Au"]] adsorbate = mol.copy() -del adsorbate[[i for i in range(len(symbols)) if symbols[i] == 'Au']] +del adsorbate[[i for i in range(len(symbols)) if symbols[i] == "Au"]] # optional: load the optimized molecules from ase import io -surface_opt = io.read('surface_opt.xyz') -adsorbate_opt = io.read('adsorbate_opt.xyz') + +surface_opt = io.read("surface_opt.xyz") +adsorbate_opt = io.read("adsorbate_opt.xyz") assert len(surface_opt) == len(surface) assert len(adsorbate_opt) == len(adsorbate) # settings for job base_settings = Settings() -base_settings.input.ams.task = 'SinglePoint' -base_settings.input.band.basis.type = 'TZP' -base_settings.input.band.basis.core = 'Medium' -base_settings.input.band.dos.calcdos = 'No' -base_settings.input.band.kspace.regular.numberofpoints = '5 5 1' -base_settings.input.band.beckegrid.quality = 'Good' -base_settings.input.band.zlmfit.quality = 'Good' +base_settings.input.ams.task = "SinglePoint" +base_settings.input.band.basis.type = "TZP" +base_settings.input.band.basis.core = "Medium" +base_settings.input.band.dos.calcdos = "No" +base_settings.input.band.kspace.regular.numberofpoints = "5 5 1" +base_settings.input.band.beckegrid.quality = "Good" +base_settings.input.band.zlmfit.quality = "Good" base_settings.input.band.usesymmetry = False -base_settings.input.band.xc.gga = 'PBE' -base_settings.input.band.xc.dispersion = 'Grimme4' +base_settings.input.band.xc.gga = "PBE" +base_settings.input.band.xc.dispersion = "Grimme4" eda_settings = Settings() -eda_settings.input.band.peda = '' +eda_settings.input.band.peda = "" -eda_job = BANDFragmentJob(fragment1=fromASE(surface), fragment2=fromASE(adsorbate), - settings=base_settings, full_settings=eda_settings, - fragment1_opt=fromASE(surface_opt), fragment2_opt=fromASE(adsorbate_opt)) +eda_job = BANDFragmentJob( + fragment1=fromASE(surface), + fragment2=fromASE(adsorbate), + settings=base_settings, + full_settings=eda_settings, + fragment1_opt=fromASE(surface_opt), + fragment2_opt=fromASE(adsorbate_opt), +) results = eda_job.run() eda_res = eda_job.results.get_energy_decomposition() -print("{:<20} {:>10}".format('Term', 'Energy [kJ/mol]')) +print("{:<20} {:>10}".format("Term", "Energy [kJ/mol]")) for key, value in eda_res.items(): - print("{:<20} {:>10.4f}".format(key, value)) \ No newline at end of file + print("{:<20} {:>10.4f}".format(key, value)) From beefc24e8da3aac3cf272d79d3cafe333772fef8 Mon Sep 17 00:00:00 2001 From: Patrick Melix Date: Thu, 12 Sep 2024 13:05:27 +0200 Subject: [PATCH 17/29] Revert "fix non-working adffragment" This reverts commit b99ff4a8b945aca9df5728066ca9281921abc175. --- recipes/adffragment.py | 43 +++++++++++++++++------------------------- 1 file changed, 17 insertions(+), 26 deletions(-) diff --git a/recipes/adffragment.py b/recipes/adffragment.py index 9644f004..11ad9648 100644 --- a/recipes/adffragment.py +++ b/recipes/adffragment.py @@ -5,6 +5,7 @@ from scm.plams.tools.units import Units from scm.plams.core.errors import FileError from scm.plams.mol.molecule import Molecule +<<<<<<< HEAD from scm.plams.core.functions import add_to_instance @@ -12,6 +13,15 @@ from os.path import abspath, isdir, basename, relpath from os import symlink from typing import Dict, List, Union +||||||| b99ff4a (fix non-working adffragment) +from scm.plams import add_to_instance + +from os.path import join as opj +from os.path import relpath +from os import symlink +from typing import List, Union +======= +>>>>>>> parent of b99ff4a (fix non-working adffragment) __all__ = ["ADFFragmentJob", "ADFFragmentResults"] @@ -111,32 +121,13 @@ def prerun(self): # noqa F811 for at in self.fragment2: at.properties.suffix = "adf.f=subsystem2" - self.children = [self.f1, self.f2] - - def new_children(self) -> Union[None, List[AMSJob]]: - """After the first round, add the full job to the children list.""" - if hasattr(self, "full"): - return None - else: - settings = self.settings + self.full_settings - settings.input.adf.fragments.subsystem1 = f"{self.f1.name}.rkf" - settings.input.adf.fragments.subsystem2 = f"{self.f2.name}.rkf" - self.full = AMSJob(name="full", molecule=self.fragment1 + self.fragment2, settings=settings) - self.full.depend += [self.f1, self.f2] - # save the fragment paths for the prerun of the full job - self.full.frag_paths = [] - for job in [self.f1, self.f2]: - self.full.frag_paths.append(job.path) - - # edit full prerun to create symlinks - @add_to_instance(self.full) - def prerun(self): - """Create symlinks for the restart files.""" - for i, job in enumerate(["frag1", "frag2"]): - rel_path = relpath(self.frag_paths[i], self.path) - symlink(opj(rel_path, "adf.rkf"), opj(self.path, f"{job}.rkf")) - - return [self.full] + self.full = AMSJob( + name="full", molecule=self.fragment1 + self.fragment2, settings=self.settings + self.full_settings + ) + + self.full.settings.input.adf.fragments.subsystem1 = (self.f1, "adf") + self.full.settings.input.adf.fragments.subsystem2 = (self.f2, "adf") + self.children = [self.f1, self.f2, self.full] @classmethod def load_external(cls, path: str, jobname: str = None) -> "ADFFragmentJob": From 3674ebaa467de78227f6f1aed1e8f97270f6cfec Mon Sep 17 00:00:00 2001 From: Patrick Melix Date: Thu, 12 Sep 2024 13:07:52 +0200 Subject: [PATCH 18/29] fix revert --- recipes/adffragment.py | 16 ++-------------- 1 file changed, 2 insertions(+), 14 deletions(-) diff --git a/recipes/adffragment.py b/recipes/adffragment.py index 11ad9648..445f2988 100644 --- a/recipes/adffragment.py +++ b/recipes/adffragment.py @@ -5,23 +5,11 @@ from scm.plams.tools.units import Units from scm.plams.core.errors import FileError from scm.plams.mol.molecule import Molecule -<<<<<<< HEAD -from scm.plams.core.functions import add_to_instance from os.path import join as opj -from os.path import abspath, isdir, basename, relpath -from os import symlink -from typing import Dict, List, Union -||||||| b99ff4a (fix non-working adffragment) -from scm.plams import add_to_instance - -from os.path import join as opj -from os.path import relpath -from os import symlink -from typing import List, Union -======= ->>>>>>> parent of b99ff4a (fix non-working adffragment) +from os.path import abspath, isdir, basename +from typing import Dict __all__ = ["ADFFragmentJob", "ADFFragmentResults"] From dd3be2d502c3880b2c2f9c9c195fbbd0dcc9272c Mon Sep 17 00:00:00 2001 From: Patrick Melix Date: Wed, 25 Sep 2024 14:37:52 +0200 Subject: [PATCH 19/29] move optimized fragments from BAND into ADF, enable frag specific settings --- recipes/adffragment.py | 55 +++++++++++++++++++++++++++++++++++++---- recipes/bandfragment.py | 24 ------------------ 2 files changed, 50 insertions(+), 29 deletions(-) diff --git a/recipes/adffragment.py b/recipes/adffragment.py index 445f2988..2555b4ec 100644 --- a/recipes/adffragment.py +++ b/recipes/adffragment.py @@ -86,23 +86,51 @@ class ADFFragmentJob(MultiJob): _result_type = ADFFragmentResults - def __init__(self, fragment1=None, fragment2=None, full_settings=None, **kwargs): - """ + def __init__(self, fragment1: Molecule = None, fragment2: Molecule = None, + fragment1_opt: Molecule = None, fragment2_opt: Molecule = None, + full_settings: Settings = None, + frag1_settings: Settings = None, frag2_settings: Settings = None, + **kwargs): + """Create an ADFFragmentJob with the given fragments. + + The optimized fragment structures can be given as arguments. + Single point calculations on these structures will then be performed + to obtain the preparation energy. + + Subclass of |MultiJob|. + Args: fragment1 (Molecule): The first fragment. fragment2 (Molecule): The second fragment. + fragment1_opt (Molecule, optional): The optimized first fragment. + fragment2_opt (Molecule, optional): The optimized second fragment. full_settings (Settings): The settings for the full calculation. + frag1_settings (Settings): Specific settings for the first fragment calculation. + frag2_settings (Settings): Specific settings for the second fragment calculation. **kwargs: Further keyword arguments for |MultiJob|. """ MultiJob.__init__(self, **kwargs) self.fragment1 = fragment1.copy() if isinstance(fragment1, Molecule) else fragment1 + self.frag1_settings = frag1_settings or Settings() self.fragment2 = fragment2.copy() if isinstance(fragment2, Molecule) else fragment2 + self.frag2_settings = frag2_settings or Settings() self.full_settings = full_settings or Settings() + if isinstance(fragment1_opt, Molecule): + self.fragment1_opt = fragment1_opt.copy() + elif fragment1_opt is not None: + raise TypeError("fragment1_opt should be a Molecule object or None") + if isinstance(fragment2_opt, Molecule): + self.fragment2_opt = fragment2_opt.copy() + elif fragment2_opt is not None: + raise TypeError("fragment2_opt should be a Molecule object or None") + def prerun(self): # noqa F811 """Prepare the fragments and the full calculation.""" - self.f1 = AMSJob(name="frag1", molecule=self.fragment1, settings=self.settings) - self.f2 = AMSJob(name="frag2", molecule=self.fragment2, settings=self.settings) + self.f1 = AMSJob(name="frag1", molecule=self.fragment1, + settings=self.settings+self.frag1_settings) + self.f2 = AMSJob(name="frag2", molecule=self.fragment2, + settings=self.settings+self.frag2_settings) for at in self.fragment1: at.properties.suffix = "adf.f=subsystem1" @@ -110,13 +138,23 @@ def prerun(self): # noqa F811 at.properties.suffix = "adf.f=subsystem2" self.full = AMSJob( - name="full", molecule=self.fragment1 + self.fragment2, settings=self.settings + self.full_settings + name="full", molecule=self.fragment1 + self.fragment2, + settings=self.settings + self.full_settings ) self.full.settings.input.adf.fragments.subsystem1 = (self.f1, "adf") self.full.settings.input.adf.fragments.subsystem2 = (self.f2, "adf") self.children = [self.f1, self.f2, self.full] + if hasattr(self, 'fragment1_opt'): + self.f1_opt = AMSJob(name="frag1_opt", molecule=self.fragment1_opt, + settings=self.settings+self.frag1_settings) + self.children.append(self.f1_opt) + if hasattr(self, 'fragment2_opt'): + self.f2_opt = AMSJob(name="frag2_opt", molecule=self.fragment2_opt, + settings=self.settings+self.frag2_settings) + self.children.append(self.f2_opt) + @classmethod def load_external(cls, path: str, jobname: str = None) -> "ADFFragmentJob": """Load the results of the ADFFragmentJob job from an external path. @@ -143,4 +181,11 @@ def load_external(cls, path: str, jobname: str = None) -> "ADFFragmentJob": job.full = AMSJob.load_external(opj(path, "full")) job.children = [job.f1, job.f2, job.full] + if isdir(opj(path, "frag1_opt")): + job.f1_opt = AMSJob.load_external(opj(path, "frag1_opt")) + job.children.append(job.f1_opt) + if isdir(opj(path, "frag2_opt")): + job.f2_opt = AMSJob.load_external(opj(path, "frag2_opt")) + job.children.append(job.f2_opt) + return job diff --git a/recipes/bandfragment.py b/recipes/bandfragment.py index 3aabafca..a0f2ee02 100644 --- a/recipes/bandfragment.py +++ b/recipes/bandfragment.py @@ -60,30 +60,6 @@ def get_energy_decomposition(self, unit="kJ/mol") -> Dict[str, float]: class BANDFragmentJob(ADFFragmentJob): _result_type = BANDFragmentResults - def __init__(self, fragment1_opt: Molecule = None, fragment2_opt: Molecule = None, **kwargs) -> None: - """Create a BANDFragmentJob with the given fragments. - - The optimized fragment structures can be given as arguments. - Single point calculations on these structures will then be performed - to obtain the preparation energy. - - Subclass of |ADFFragmentJob|. - - Args: - fragment1_opt (Molecule, optional): The optimized fragment 1. - fragment2_opt (Molecule, optional): The optimized fragment 2. - **kwargs: Further keyword arguments for |ADFFragmentJob|. - """ - super().__init__(**kwargs) - if isinstance(fragment1_opt, Molecule): - self.fragment1_opt = fragment1_opt.copy() - self.f1_opt = AMSJob(name="frag1_opt", molecule=self.fragment1_opt, settings=self.settings) - self.children.append(self.f1_opt) - if isinstance(fragment2_opt, Molecule): - self.fragment2_opt = fragment2_opt.copy() - self.f2_opt = AMSJob(name="frag2_opt", molecule=self.fragment2_opt, settings=self.settings) - self.children.append(self.f2_opt) - def create_mapping_setting(self) -> None: if "fragment" in self.full_settings.input.band: log( From 3e52686ffd6bf7fac1dd75044874ac869d6c13b4 Mon Sep 17 00:00:00 2001 From: Patrick Melix Date: Tue, 22 Oct 2024 13:28:06 +0200 Subject: [PATCH 20/29] adffragment check routine that does something useful --- recipes/adffragment.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/recipes/adffragment.py b/recipes/adffragment.py index 2555b4ec..8e1f58ce 100644 --- a/recipes/adffragment.py +++ b/recipes/adffragment.py @@ -16,6 +16,17 @@ class ADFFragmentResults(Results): + def check(self): + """Check if the calculation finished successfully. + + Overwriting the method of |MultiJob| because it only checks if the job + finished, not how it finished. + """ + for job in self.job.children: + if not job.check(): + return False + return True + def get_properties(self): """Redirect to |ADFResults| of the full calculation.""" return self.job.full.results.get_properties() From b9db468ca8d8de334d95466bd71fa20a8cf0bd10 Mon Sep 17 00:00:00 2001 From: Patrick Melix Date: Tue, 22 Oct 2024 13:30:20 +0200 Subject: [PATCH 21/29] run black on adffragment --- recipes/adffragment.py | 96 ++++++++++++++++++++++++++++++------------ 1 file changed, 69 insertions(+), 27 deletions(-) diff --git a/recipes/adffragment.py b/recipes/adffragment.py index 8e1f58ce..bff1b4e8 100644 --- a/recipes/adffragment.py +++ b/recipes/adffragment.py @@ -18,7 +18,7 @@ class ADFFragmentResults(Results): def check(self): """Check if the calculation finished successfully. - + Overwriting the method of |MultiJob| because it only checks if the job finished, not how it finished. """ @@ -58,7 +58,12 @@ def get_energy_decomposition(self, unit="kJ/mol") -> Dict[str, float]: """ energy_section = self.job.full.results.read_rkf_section("Energy", file="adf") ret = {} - for k in ["Electrostatic Energy", "Kinetic Energy", "Elstat Interaction", "XC Energy"]: + for k in [ + "Electrostatic Energy", + "Kinetic Energy", + "Elstat Interaction", + "XC Energy", + ]: ret[k] = Units.convert(energy_section[k], "au", unit) # most information not available from the KF file @@ -67,13 +72,27 @@ def get_energy_decomposition(self, unit="kJ/mol") -> Dict[str, float]: res2 = self.job.f2.results pos = -4 # position of the energy in au the output # E_int appears in a comment below the PEDA Table of the output - ret["E_int"] = Units.convert(float(res.grep_output("Total Bonding Energy:")[-2].split()[pos]), "au", unit) - ret["E_int_disp"] = Units.convert(float(res.grep_output("Dispersion Energy:")[-1].split()[pos]), "au", unit) - ret["E_Pauli"] = Units.convert(float(res.grep_output("Pauli Repulsion (Delta")[-1].split()[pos]), "au", unit) + ret["E_int"] = Units.convert( + float(res.grep_output("Total Bonding Energy:")[-2].split()[pos]), "au", unit + ) + ret["E_int_disp"] = Units.convert( + float(res.grep_output("Dispersion Energy:")[-1].split()[pos]), "au", unit + ) + ret["E_Pauli"] = Units.convert( + float(res.grep_output("Pauli Repulsion (Delta")[-1].split()[pos]), + "au", + unit, + ) ret["E_elstat"] = Units.convert( - float(res.grep_output("Electrostatic Interaction:")[-1].split()[pos]), "au", unit + float(res.grep_output("Electrostatic Interaction:")[-1].split()[pos]), + "au", + unit, + ) + ret["E_orb"] = Units.convert( + float(res.grep_output("Total Orbital Interactions:")[-1].split()[pos]), + "au", + unit, ) - ret["E_orb"] = Units.convert(float(res.grep_output("Total Orbital Interactions:")[-1].split()[pos]), "au", unit) ret["E_1"] = res1.get_energy(unit=unit) ret["E_2"] = res2.get_energy(unit=unit) @@ -97,11 +116,17 @@ class ADFFragmentJob(MultiJob): _result_type = ADFFragmentResults - def __init__(self, fragment1: Molecule = None, fragment2: Molecule = None, - fragment1_opt: Molecule = None, fragment2_opt: Molecule = None, - full_settings: Settings = None, - frag1_settings: Settings = None, frag2_settings: Settings = None, - **kwargs): + def __init__( + self, + fragment1: Molecule = None, + fragment2: Molecule = None, + fragment1_opt: Molecule = None, + fragment2_opt: Molecule = None, + full_settings: Settings = None, + frag1_settings: Settings = None, + frag2_settings: Settings = None, + **kwargs + ): """Create an ADFFragmentJob with the given fragments. The optimized fragment structures can be given as arguments. @@ -109,7 +134,7 @@ def __init__(self, fragment1: Molecule = None, fragment2: Molecule = None, to obtain the preparation energy. Subclass of |MultiJob|. - + Args: fragment1 (Molecule): The first fragment. fragment2 (Molecule): The second fragment. @@ -121,9 +146,13 @@ def __init__(self, fragment1: Molecule = None, fragment2: Molecule = None, **kwargs: Further keyword arguments for |MultiJob|. """ MultiJob.__init__(self, **kwargs) - self.fragment1 = fragment1.copy() if isinstance(fragment1, Molecule) else fragment1 + self.fragment1 = ( + fragment1.copy() if isinstance(fragment1, Molecule) else fragment1 + ) self.frag1_settings = frag1_settings or Settings() - self.fragment2 = fragment2.copy() if isinstance(fragment2, Molecule) else fragment2 + self.fragment2 = ( + fragment2.copy() if isinstance(fragment2, Molecule) else fragment2 + ) self.frag2_settings = frag2_settings or Settings() self.full_settings = full_settings or Settings() @@ -138,10 +167,16 @@ def __init__(self, fragment1: Molecule = None, fragment2: Molecule = None, def prerun(self): # noqa F811 """Prepare the fragments and the full calculation.""" - self.f1 = AMSJob(name="frag1", molecule=self.fragment1, - settings=self.settings+self.frag1_settings) - self.f2 = AMSJob(name="frag2", molecule=self.fragment2, - settings=self.settings+self.frag2_settings) + self.f1 = AMSJob( + name="frag1", + molecule=self.fragment1, + settings=self.settings + self.frag1_settings, + ) + self.f2 = AMSJob( + name="frag2", + molecule=self.fragment2, + settings=self.settings + self.frag2_settings, + ) for at in self.fragment1: at.properties.suffix = "adf.f=subsystem1" @@ -149,21 +184,28 @@ def prerun(self): # noqa F811 at.properties.suffix = "adf.f=subsystem2" self.full = AMSJob( - name="full", molecule=self.fragment1 + self.fragment2, - settings=self.settings + self.full_settings + name="full", + molecule=self.fragment1 + self.fragment2, + settings=self.settings + self.full_settings, ) self.full.settings.input.adf.fragments.subsystem1 = (self.f1, "adf") self.full.settings.input.adf.fragments.subsystem2 = (self.f2, "adf") self.children = [self.f1, self.f2, self.full] - if hasattr(self, 'fragment1_opt'): - self.f1_opt = AMSJob(name="frag1_opt", molecule=self.fragment1_opt, - settings=self.settings+self.frag1_settings) + if hasattr(self, "fragment1_opt"): + self.f1_opt = AMSJob( + name="frag1_opt", + molecule=self.fragment1_opt, + settings=self.settings + self.frag1_settings, + ) self.children.append(self.f1_opt) - if hasattr(self, 'fragment2_opt'): - self.f2_opt = AMSJob(name="frag2_opt", molecule=self.fragment2_opt, - settings=self.settings+self.frag2_settings) + if hasattr(self, "fragment2_opt"): + self.f2_opt = AMSJob( + name="frag2_opt", + molecule=self.fragment2_opt, + settings=self.settings + self.frag2_settings, + ) self.children.append(self.f2_opt) @classmethod From f4a60d1492c1eeb350f8e95c79af9a2126623a49 Mon Sep 17 00:00:00 2001 From: Patrick Melix Date: Tue, 22 Oct 2024 13:41:00 +0200 Subject: [PATCH 22/29] now with correct black options --- recipes/adffragment.py | 18 +++++------------- 1 file changed, 5 insertions(+), 13 deletions(-) diff --git a/recipes/adffragment.py b/recipes/adffragment.py index bff1b4e8..f412e5fb 100644 --- a/recipes/adffragment.py +++ b/recipes/adffragment.py @@ -72,12 +72,8 @@ def get_energy_decomposition(self, unit="kJ/mol") -> Dict[str, float]: res2 = self.job.f2.results pos = -4 # position of the energy in au the output # E_int appears in a comment below the PEDA Table of the output - ret["E_int"] = Units.convert( - float(res.grep_output("Total Bonding Energy:")[-2].split()[pos]), "au", unit - ) - ret["E_int_disp"] = Units.convert( - float(res.grep_output("Dispersion Energy:")[-1].split()[pos]), "au", unit - ) + ret["E_int"] = Units.convert(float(res.grep_output("Total Bonding Energy:")[-2].split()[pos]), "au", unit) + ret["E_int_disp"] = Units.convert(float(res.grep_output("Dispersion Energy:")[-1].split()[pos]), "au", unit) ret["E_Pauli"] = Units.convert( float(res.grep_output("Pauli Repulsion (Delta")[-1].split()[pos]), "au", @@ -125,7 +121,7 @@ def __init__( full_settings: Settings = None, frag1_settings: Settings = None, frag2_settings: Settings = None, - **kwargs + **kwargs, ): """Create an ADFFragmentJob with the given fragments. @@ -146,13 +142,9 @@ def __init__( **kwargs: Further keyword arguments for |MultiJob|. """ MultiJob.__init__(self, **kwargs) - self.fragment1 = ( - fragment1.copy() if isinstance(fragment1, Molecule) else fragment1 - ) + self.fragment1 = fragment1.copy() if isinstance(fragment1, Molecule) else fragment1 self.frag1_settings = frag1_settings or Settings() - self.fragment2 = ( - fragment2.copy() if isinstance(fragment2, Molecule) else fragment2 - ) + self.fragment2 = fragment2.copy() if isinstance(fragment2, Molecule) else fragment2 self.frag2_settings = frag2_settings or Settings() self.full_settings = full_settings or Settings() From 483b7924e8fc66ed32300e4fe2f043ecb3020105 Mon Sep 17 00:00:00 2001 From: Patrick Melix Date: Tue, 22 Oct 2024 13:44:13 +0200 Subject: [PATCH 23/29] remove unused import --- recipes/bandfragment.py | 1 - 1 file changed, 1 deletion(-) diff --git a/recipes/bandfragment.py b/recipes/bandfragment.py index a0f2ee02..92ce9844 100644 --- a/recipes/bandfragment.py +++ b/recipes/bandfragment.py @@ -4,7 +4,6 @@ from scm.plams.recipes.adffragment import ADFFragmentJob, ADFFragmentResults from scm.plams.tools.units import Units from scm.plams.core.errors import FileError -from scm.plams.mol.molecule import Molecule from scm.plams.core.functions import add_to_instance From cb5f32cb64cdfd0e85ffcfc124e16f6d65df5949 Mon Sep 17 00:00:00 2001 From: Patrick Melix Date: Fri, 22 Nov 2024 13:37:10 +0100 Subject: [PATCH 24/29] adffragment add eigenvalue and oi extraction --- recipes/adffragment.py | 57 +++++++++++++++++++++++++++++++++++++++++- 1 file changed, 56 insertions(+), 1 deletion(-) diff --git a/recipes/adffragment.py b/recipes/adffragment.py index f412e5fb..d66bed91 100644 --- a/recipes/adffragment.py +++ b/recipes/adffragment.py @@ -9,7 +9,7 @@ from os.path import join as opj from os.path import abspath, isdir, basename -from typing import Dict +from typing import Dict, List, Tuple, Union __all__ = ["ADFFragmentJob", "ADFFragmentResults"] @@ -106,6 +106,61 @@ def get_energy_decomposition(self, unit="kJ/mol") -> Dict[str, float]: return ret + def get_nocv_eigenvalues(self) -> Union[List[float], Tuple[List[float]]]: + """Get the NOCV eigenvalues of the full calculation.""" + keys = self.job.full.results.get_rkf_skeleton(file='adf') + if 'NOCV' not in keys: + raise KeyError("NOCV section not found in the rkf file") + keys = keys['NOCV'] + # restricted calculation + if 'NOCV_eigenvalues_restricted' in keys: + return self.job.full.results.readrkf('NOCV', + 'NOCV_eigenvalues_restricted', file='adf') + # unrestricted calculation + elif 'NOCV_eigenvalues_alpha' in keys: + tpl = (self.job.full.results.readrkf('NOCV', + 'NOCV_eigenvalues_alpha', file='adf'), + self.job.full.results.readrkf('NOCV', + 'NOCV_eigenvalues_beta', file='adf')) + return tpl + else: + raise KeyError("NOCV eigenvalues not found in the rkf file") + + def get_nocv_orbital_interaction(self, unit='kcal/mol' + ) -> Union[List[float], + Tuple[List[float]]]: + """Get the NOCV orbital interactions of the full calculation.""" + def _calc_energies(oi) -> List[float]: + ret = [] + for i in range(int(len(oi)/2)): + a = oi[i] + b = oi[len(oi) - i - 1] + ret.append(a + b) + return ret + + keys = self.job.full.results.get_rkf_skeleton(file='adf') + if 'NOCV' not in keys: + raise KeyError("NOCV section not found in the rkf file") + keys = keys['NOCV'] + # restricted calculation + if 'NOCV_oi_restricted' in keys: + oi = self.job.full.results.readrkf('NOCV', + 'NOCV_oi_restricted', file='adf') + # unrestricted calculation + elif 'NOCV_oi_alpha' in keys: + tpl = (self.job.full.results.readrkf('NOCV', + 'NOCV_oi_alpha', file='adf'), + self.job.full.results.readrkf('NOCV', + 'NOCV_oi_beta', file='adf')) + else: + raise KeyError("NOCV orbital interaction not found in rkf file") + if isinstance(oi, tuple): + energies = (_calc_energies(tpl[0]), _calc_energies(tpl[1])) + else: + energies = _calc_energies(oi) + return Units.convert(energies, 'kcal/mol', unit) + + class ADFFragmentJob(MultiJob): """Subclass of |MultiJob| for ADF fragment calculations.""" From 5e8947ada3105952e662e66881b9146eb7fe82dd Mon Sep 17 00:00:00 2001 From: Patrick Melix Date: Wed, 4 Dec 2024 13:56:34 +0100 Subject: [PATCH 25/29] add check_nocv_eigenvalues --- recipes/adffragment.py | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/recipes/adffragment.py b/recipes/adffragment.py index d66bed91..8ab5ab79 100644 --- a/recipes/adffragment.py +++ b/recipes/adffragment.py @@ -6,6 +6,7 @@ from scm.plams.core.errors import FileError from scm.plams.mol.molecule import Molecule +import numpy as np from os.path import join as opj from os.path import abspath, isdir, basename @@ -126,6 +127,33 @@ def get_nocv_eigenvalues(self) -> Union[List[float], Tuple[List[float]]]: else: raise KeyError("NOCV eigenvalues not found in the rkf file") + def check_nocv_eigenvalues(self) -> bool: + """Basic NOCV check. + + The eigenvalue #i should be equal to -eigenvalue #-i. + + Returns: + bool: True if the check passes, False otherwise. + """ + def _check_list(eigenvalues): + for i in range(int(len(eigenvalues)/2)): + if not np.isclose(eigenvalues[i], -1*eigenvalues[-i-1], atol=1e-4): + return False + return True + + NOCV_eigenvalues = self.job.full.results.get_nocv_eigenvalues() + # spin-restricted calculation + if isinstance(NOCV_eigenvalues, list): + return _check_list(NOCV_eigenvalues) + elif isinstance(NOCV_eigenvalues, tuple): + for eigenvalues in NOCV_eigenvalues: + if not _check_list(eigenvalues): + return False + return True + else: + raise TypeError("Unexpected type for NOCV eigenvalues") + + def get_nocv_orbital_interaction(self, unit='kcal/mol' ) -> Union[List[float], Tuple[List[float]]]: From 8e00e64b63110053193c97249101771eb43a8122 Mon Sep 17 00:00:00 2001 From: Patrick Melix Date: Wed, 4 Dec 2024 14:20:27 +0100 Subject: [PATCH 26/29] fix bad .full. insertion --- recipes/adffragment.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/recipes/adffragment.py b/recipes/adffragment.py index 8ab5ab79..125d32cf 100644 --- a/recipes/adffragment.py +++ b/recipes/adffragment.py @@ -141,7 +141,7 @@ def _check_list(eigenvalues): return False return True - NOCV_eigenvalues = self.job.full.results.get_nocv_eigenvalues() + NOCV_eigenvalues = self.job.results.get_nocv_eigenvalues() # spin-restricted calculation if isinstance(NOCV_eigenvalues, list): return _check_list(NOCV_eigenvalues) From ae452605c09a4a676389cb09de355d6d911e6586 Mon Sep 17 00:00:00 2001 From: Patrick Melix Date: Wed, 4 Dec 2024 14:34:36 +0100 Subject: [PATCH 27/29] run black --- recipes/adffragment.py | 62 ++++++++++++++++++++---------------------- 1 file changed, 29 insertions(+), 33 deletions(-) diff --git a/recipes/adffragment.py b/recipes/adffragment.py index 125d32cf..0cdfef87 100644 --- a/recipes/adffragment.py +++ b/recipes/adffragment.py @@ -109,38 +109,38 @@ def get_energy_decomposition(self, unit="kJ/mol") -> Dict[str, float]: def get_nocv_eigenvalues(self) -> Union[List[float], Tuple[List[float]]]: """Get the NOCV eigenvalues of the full calculation.""" - keys = self.job.full.results.get_rkf_skeleton(file='adf') - if 'NOCV' not in keys: + keys = self.job.full.results.get_rkf_skeleton(file="adf") + if "NOCV" not in keys: raise KeyError("NOCV section not found in the rkf file") - keys = keys['NOCV'] + keys = keys["NOCV"] # restricted calculation - if 'NOCV_eigenvalues_restricted' in keys: - return self.job.full.results.readrkf('NOCV', - 'NOCV_eigenvalues_restricted', file='adf') + if "NOCV_eigenvalues_restricted" in keys: + return self.job.full.results.readrkf("NOCV", "NOCV_eigenvalues_restricted", file="adf") # unrestricted calculation - elif 'NOCV_eigenvalues_alpha' in keys: - tpl = (self.job.full.results.readrkf('NOCV', - 'NOCV_eigenvalues_alpha', file='adf'), - self.job.full.results.readrkf('NOCV', - 'NOCV_eigenvalues_beta', file='adf')) + elif "NOCV_eigenvalues_alpha" in keys: + tpl = ( + self.job.full.results.readrkf("NOCV", "NOCV_eigenvalues_alpha", file="adf"), + self.job.full.results.readrkf("NOCV", "NOCV_eigenvalues_beta", file="adf"), + ) return tpl else: raise KeyError("NOCV eigenvalues not found in the rkf file") def check_nocv_eigenvalues(self) -> bool: """Basic NOCV check. - + The eigenvalue #i should be equal to -eigenvalue #-i. Returns: bool: True if the check passes, False otherwise. """ + def _check_list(eigenvalues): - for i in range(int(len(eigenvalues)/2)): - if not np.isclose(eigenvalues[i], -1*eigenvalues[-i-1], atol=1e-4): + for i in range(int(len(eigenvalues) / 2)): + if not np.isclose(eigenvalues[i], -1 * eigenvalues[-i - 1], atol=1e-4): return False return True - + NOCV_eigenvalues = self.job.results.get_nocv_eigenvalues() # spin-restricted calculation if isinstance(NOCV_eigenvalues, list): @@ -153,41 +153,37 @@ def _check_list(eigenvalues): else: raise TypeError("Unexpected type for NOCV eigenvalues") - - def get_nocv_orbital_interaction(self, unit='kcal/mol' - ) -> Union[List[float], - Tuple[List[float]]]: + def get_nocv_orbital_interaction(self, unit="kcal/mol") -> Union[List[float], Tuple[List[float]]]: """Get the NOCV orbital interactions of the full calculation.""" + def _calc_energies(oi) -> List[float]: ret = [] - for i in range(int(len(oi)/2)): + for i in range(int(len(oi) / 2)): a = oi[i] b = oi[len(oi) - i - 1] ret.append(a + b) return ret - keys = self.job.full.results.get_rkf_skeleton(file='adf') - if 'NOCV' not in keys: + keys = self.job.full.results.get_rkf_skeleton(file="adf") + if "NOCV" not in keys: raise KeyError("NOCV section not found in the rkf file") - keys = keys['NOCV'] + keys = keys["NOCV"] # restricted calculation - if 'NOCV_oi_restricted' in keys: - oi = self.job.full.results.readrkf('NOCV', - 'NOCV_oi_restricted', file='adf') + if "NOCV_oi_restricted" in keys: + oi = self.job.full.results.readrkf("NOCV", "NOCV_oi_restricted", file="adf") # unrestricted calculation - elif 'NOCV_oi_alpha' in keys: - tpl = (self.job.full.results.readrkf('NOCV', - 'NOCV_oi_alpha', file='adf'), - self.job.full.results.readrkf('NOCV', - 'NOCV_oi_beta', file='adf')) + elif "NOCV_oi_alpha" in keys: + tpl = ( + self.job.full.results.readrkf("NOCV", "NOCV_oi_alpha", file="adf"), + self.job.full.results.readrkf("NOCV", "NOCV_oi_beta", file="adf"), + ) else: raise KeyError("NOCV orbital interaction not found in rkf file") if isinstance(oi, tuple): energies = (_calc_energies(tpl[0]), _calc_energies(tpl[1])) else: energies = _calc_energies(oi) - return Units.convert(energies, 'kcal/mol', unit) - + return Units.convert(energies, "kcal/mol", unit) class ADFFragmentJob(MultiJob): From d3af82007cc8b826918ec7bf864d5a6e5ebe2fec Mon Sep 17 00:00:00 2001 From: Patrick Melix Date: Wed, 4 Dec 2024 14:42:08 +0100 Subject: [PATCH 28/29] fixing doc issues --- recipes/adffragment.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/recipes/adffragment.py b/recipes/adffragment.py index 0cdfef87..d61f06ee 100644 --- a/recipes/adffragment.py +++ b/recipes/adffragment.py @@ -29,23 +29,23 @@ def check(self): return True def get_properties(self): - """Redirect to |ADFResults| of the full calculation.""" + """Redirect to |AMSResults| of the full calculation.""" return self.job.full.results.get_properties() def get_main_molecule(self): - """Redirect to |ADFResults| of the full calculation.""" + """Redirect to |AMSResults| of the full calculation.""" return self.job.full.results.get_main_molecule() def get_input_molecule(self): - """Redirect to |ADFResults| of the full calculation.""" + """Redirect to |AMSResults| of the full calculation.""" return self.job.full.results.get_input_molecule() def get_energy(self, unit="au"): - """Redirect to |ADFResults| of the full calculation.""" + """Redirect to |AMSResults| of the full calculation.""" return self.job.full.results.get_energy(unit) def get_dipole_vector(self, unit="au"): - """Redirect to |ADFResults| of the full calculation.""" + """Redirect to |AMSResults| of the full calculation.""" return self.job.full.results.get_dipole_vector(unit) def get_energy_decomposition(self, unit="kJ/mol") -> Dict[str, float]: @@ -211,10 +211,10 @@ def __init__( Subclass of |MultiJob|. Args: - fragment1 (Molecule): The first fragment. - fragment2 (Molecule): The second fragment. - fragment1_opt (Molecule, optional): The optimized first fragment. - fragment2_opt (Molecule, optional): The optimized second fragment. + fragment1 (|Molecule|): The first fragment. + fragment2 (|Molecule|): The second fragment. + fragment1_opt (|Molecule|, optional): The optimized first fragment. + fragment2_opt (|Molecule|, optional): The optimized second fragment. full_settings (Settings): The settings for the full calculation. frag1_settings (Settings): Specific settings for the first fragment calculation. frag2_settings (Settings): Specific settings for the second fragment calculation. From 067fe8ed3f3bcb6408beea61c6eba2469ad8ab9a Mon Sep 17 00:00:00 2001 From: Patrick Melix Date: Wed, 4 Dec 2024 14:48:57 +0100 Subject: [PATCH 29/29] fix doc relative link --- doc/source/examples/bandfragment.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/source/examples/bandfragment.rst b/doc/source/examples/bandfragment.rst index aff16dd2..f3ba207b 100644 --- a/doc/source/examples/bandfragment.rst +++ b/doc/source/examples/bandfragment.rst @@ -33,6 +33,6 @@ API Example ~~~~~~~~~~~~~~~~~~~~~~~~~ -A basic example using `ASE to build a surface slab `_ and perform a BAND fragment calculation: :download:`bandfrag_test.py <../../../examples/scripting/bandfrag_test.py>` +A basic example using `ASE to build a surface slab `_ and perform a BAND fragment calculation: :download:`bandfrag_test.py <../../../examples/bandfrag.py>` -.. literalinclude:: ../../../examples/scripting/bandfrag_test.py \ No newline at end of file +.. literalinclude:: ../../../examples/bandfrag.py \ No newline at end of file