ADF fragment job

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.

We define a new job type ADFFragmentJob by extending 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.

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

from scm.plams.core.basejob import MultiJob
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.mol.molecule import Molecule

__all__ = ["ADFFragmentJob", "ADFFragmentResults"]


class ADFFragmentResults(Results):

    def get_properties(self):
        return self.job.full.results.get_properties()

    def get_main_molecule(self):
        return self.job.full.results.get_main_molecule()

    def get_input_molecule(self):
        return self.job.full.results.get_input_molecule()

    def get_energy(self, unit="au"):
        return self.job.full.results.get_energy(unit)

    def get_dipole_vector(self, unit="au"):
        return self.job.full.results.get_dipole_vector(unit)

    def get_energy_decomposition(self):
        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]
        return ret


class ADFFragmentJob(MultiJob):
    _result_type = ADFFragmentResults

    def __init__(self, fragment1=None, fragment2=None, full_settings=None, **kwargs):
        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
        self.f1 = AMSJob(name="frag1", molecule=self.fragment1, settings=self.settings)
        self.f2 = AMSJob(name="frag2", molecule=self.fragment2, settings=self.settings)

        for at in self.fragment1:
            at.properties.suffix = "adf.f=subsystem1"
        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, self.full]

An example usage: (ethene.xyz, butadiene.xyz, adffrag_test.py)

6

C      -0.75086900       1.37782400      -2.43303700
C      -0.05392100       2.51281000      -2.41769100
H      -1.78964800       1.33942600      -2.09651100
H      -0.30849400       0.43896500      -2.76734700
H      -0.49177100       3.45043100      -2.06789100
H       0.98633900       2.54913500      -2.74329400

10

C       0.14667300      -0.21503500       0.40053800
C       1.45297400      -0.07836900       0.12424400
C       2.23119700       1.15868100       0.12912100
C       1.78331500       2.39701500       0.38779700
H      -0.48348000       0.63110600       0.67664100
H      -0.33261900      -1.19332100       0.35411600
H       2.01546300      -0.97840100      -0.14506700
H       3.29046200       1.03872500      -0.12139700
H       2.45728900       3.25301000       0.35150400
H       0.74193400       2.60120700       0.64028800

from scm.plams import Settings, Molecule
from scm.plams.recipes.adffragment import ADFFragmentJob

common = Settings()  # common settings for all 3 jobs
common.input.ams.Task = "SinglePoint"
common.input.adf.basis.type = "DZP"
common.input.adf.xc.gga = "PBE"
common.input.adf.symmetry = "NOSYM"

full = Settings()  # additional settings for full system calculation
full.input.adf.etsnocv  # empty block
full.input.adf.print = "etslowdin"

mol1 = Molecule("ethene.xyz")
mol2 = Molecule("butadiene.xyz")

j = ADFFragmentJob(fragment1=mol1, fragment2=mol2, settings=common, full_settings=full)
r = j.run()

print("Energy decomposition:")
decom = r.get_energy_decomposition()
for k in decom:
    print(k, decom[k])
print(j.full.results.readrkf("NOCV", "NOCV_eigenvalues_restricted", "engine"))