Screening Molecules for Promising Electronic Excitations¶

Identifying systems with certain physical properties in a large molecular database is a common automation task. This example shows how to implement such a screening workflow with a simple PLAMS script.

In this toy study, the goal is to find molecules that absorb light within a target energy window. A straightforward approach would be to run time-dependent DFT with ADF for every molecule, but that quickly becomes expensive for large databases.

A more efficient strategy is to pre-screen the structures with DFTB and only rerun the promising candidates with ADF. The workflow therefore first uses TD-DFTB to optimize all molecules and keep only candidates with excitations in a broad energy window and nonzero oscillator strength, and then refines that subset with ADF TD-DFT to identify the molecules that actually absorb in the desired range.

The focus here is on the scripting pattern rather than on choosing production-quality spectroscopy settings.

Downloads: Notebook | Script ?

Requires: AMS2027 or later

Initial Imports¶

from pathlib import Path
from scm.base import ChemicalSystem, Units
from scm.inputs import ADF, AMS, DFTB
from scm.plams import AMSResults, AMSJob

Helper Functions¶

Set up a couple of useful functions for extracting results.

def get_excitations(results):
    """Returns excitation energies (in eV) and oscillator strengths (in Debye)."""
    if results.job.ok():
        exci_energies_au = results.readrkf("Excitations SS A", "excenergies", file="engine")
        oscillator_str_au = results.readrkf(
            "Excitations SS A", "oscillator strengths", file="engine"
        )
        # The results are stored in atomic units. Convert them to more convenient units:
        au_to_ev = Units.conversion_factor("au", "eV")
        au_to_debye = Units.conversion_factor("au", "Debye")
        exci_energies = [energy * au_to_ev for energy in exci_energies_au]
        oscillator_str = [strength * au_to_debye for strength in oscillator_str_au]
        return exci_energies, oscillator_str
    else:
        return [], []

def has_good_excitations(results, min_energy, max_energy, oscillator_str_threshold=1e-4):
    """Returns True if there is at least one excitation with non-vanishing oscillator strength
    in the energy range [min_energy, max_energy]. Unit for min_energy and max energy: eV."""
    exci_energies, oscillator_str = get_excitations(results)
    for e, o in zip(exci_energies, oscillator_str):
        if min_energy < e < max_energy and o > oscillator_str_threshold:
            return True
    return False

Calculation settings¶

Configure the settings for the various jobs.

# Settings for geometry optimization with the AMS driver:
go_sett = AMS(Task="GeometryOptimization")
go_sett.GeometryOptimization.Convergence.Gradients = 1.0e-4

# Settings for single point calculation with the AMS driver
sp_sett = AMS(Task="SinglePoint")

# Settings for the DFTB engine (including excitations)
dftb_sett = AMS()
dftb_sett.Engine = DFTB(Model="SCC-DFTB", ResourcesDir="QUASINANO2015")
dftb_sett.Engine.Properties.Excitations.TDDFTB.Calc = "singlet"
dftb_sett.Engine.Properties.Excitations.TDDFTB.Lowest = 10
dftb_sett.Engine.Occupation.Temperature = 5.0

# Settings for the geometry optimization with the ADF engine
adf_sett = AMS()
adf_sett.Engine = ADF()
adf_sett.Engine.Basis.Type = "DZP"
adf_sett.Engine.NumericalQuality = "Basic"

# Settings for the excitation calculation using the ADF engine
adf_exci_sett = AMS()
adf_exci_sett.Engine = ADF()
adf_exci_sett.Engine.Basis.Type = "TZP"
adf_exci_sett.Engine.XC.GGA = "PBE"
adf_exci_sett.Engine.NumericalQuality = "Basic"
adf_exci_sett.Engine.Symmetry = "NoSym"
adf_exci_sett.Engine.Excitations.Lowest = [10]
adf_exci_sett.Engine.Excitations.OnlySing = True

Load Molecules¶

Import all xyz files in the folder ‘molecules’.

from pathlib import Path

molecules_dir = Path("molecules")
molecules_dir.mkdir(exist_ok=True)

alh3_xyz = molecules_dir / "AlH3.xyz"
cscl2_xyz = molecules_dir / "CSCl2.xyz"
h20_xyz = molecules_dir / "H2O.xyz"
nh3_xyz = molecules_dir / "NH3.xyz"
s2cl2_xyz = molecules_dir / "S2Cl2.xyz"

alh3_xyz.write_text("""\
4

Al  0.0  0.0  0.0
H   1.0  0.0  0.0
H   0.5  1.0  0.0
H   0.5 -1.0  0.0
""")

cscl2_xyz.write_text("""\
4

C   0.0  0.0  0.0
S   0.0  0.0 -1.3
Cl  0.0  1.0  1.0
Cl  0.0 -1.0  1.0
""")

h20_xyz.write_text("""\
3

O 0.0 0.0 0.0
H 1.0 0.0 0.0
H 0.0 1.0 0.0
""")

nh3_xyz.write_text("""\
4

N 0.0 0.0 0.1
H 1.0 0.0 0.0
H 0.5 1.0 0.0
H 0.5 -1.0 0.0
""")

s2cl2_xyz.write_text("""\
4

S   0.0  0.0  1.3
S   0.0  0.0 -1.3
Cl  0.1  0.6  2.3
Cl  0.1 -0.6 -2.3
""")

molecules_dir = Path("molecules")
molecules = {path.stem: ChemicalSystem.from_xyz(str(path)) for path in sorted(molecules_dir.glob("*.xyz"))}

DFTB Prescreen¶

Perform an initial prescreen of all molecules with DFTB.

promising_molecules = {}

for name, mol in molecules.items():
    dftb_job = AMSJob(name="DFTB_" + name, molecule=mol, settings=go_sett + dftb_sett)
    dftb_results = dftb_job.run()

    if has_good_excitations(dftb_results, 1, 6):
        promising_molecules[name] = dftb_results.get_main_system()

[19.10|12:19:12] JOB DFTB_AlH3 STARTED
[19.10|12:19:12] JOB DFTB_AlH3 RUNNING
[19.10|12:19:13] JOB DFTB_AlH3 FINISHED
[19.10|12:19:13] JOB DFTB_AlH3 SUCCESSFUL
[19.10|12:19:13] JOB DFTB_CSCl2 STARTED
[19.10|12:19:13] JOB DFTB_CSCl2 RUNNING
[19.10|12:19:14] JOB DFTB_CSCl2 FINISHED
[19.10|12:19:14] JOB DFTB_CSCl2 SUCCESSFUL
[19.10|12:19:14] JOB DFTB_H2O STARTED
[19.10|12:19:14] JOB DFTB_H2O RUNNING
[19.10|12:19:15] JOB DFTB_H2O FINISHED
[19.10|12:19:15] JOB DFTB_H2O SUCCESSFUL
[19.10|12:19:15] JOB DFTB_NH3 STARTED
[19.10|12:19:15] JOB DFTB_NH3 RUNNING
[19.10|12:19:17] JOB DFTB_NH3 FINISHED
[19.10|12:19:17] JOB DFTB_NH3 SUCCESSFUL
[19.10|12:19:17] JOB DFTB_S2Cl2 STARTED
[19.10|12:19:17] JOB DFTB_S2Cl2 RUNNING
[19.10|12:19:18] JOB DFTB_S2Cl2 FINISHED
[19.10|12:19:18] JOB DFTB_S2Cl2 SUCCESSFUL

print(f"Found {len(promising_molecules)} promising molecules with DFTB")

Found 2 promising molecules with DFTB

Optimization and excitations calculation with ADF¶

For each of the molecules identified in the prescreen, run a further calculation with ADF.

for name, mol in promising_molecules.items():
    adf_go_job = AMSJob(name="ADF_GO_" + name, molecule=mol, settings=go_sett + adf_sett)
    adf_go_job.run()

    optimized_mol = adf_go_job.results.get_main_system()

    adf_exci_job = AMSJob(
        name="ADF_exci_" + name, molecule=optimized_mol, settings=sp_sett + adf_exci_sett
    )
    adf_exci_results = adf_exci_job.run()

    if has_good_excitations(adf_exci_results, 2, 4):
        print(f"Molecule {name} has excitation(s) satisfying our criteria!")
        print(optimized_mol)
        exci_energies, oscillator_str = get_excitations(adf_exci_results)
        print("Excitation energy [eV], oscillator strength:")
        for e, o in zip(exci_energies, oscillator_str):
            print(f"{e:8.4f}, {o:8.4f}")

[19.10|12:19:18] JOB ADF_GO_CSCl2 STARTED
[19.10|12:19:18] JOB ADF_GO_CSCl2 RUNNING
[19.10|12:19:24] JOB ADF_GO_CSCl2 FINISHED
[19.10|12:19:24] JOB ADF_GO_CSCl2 SUCCESSFUL
[19.10|12:19:24] JOB ADF_exci_CSCl2 STARTED
[19.10|12:19:24] JOB ADF_exci_CSCl2 RUNNING
[19.10|12:19:29] JOB ADF_exci_CSCl2 FINISHED
[19.10|12:19:29] JOB ADF_exci_CSCl2 SUCCESSFUL
[19.10|12:19:29] JOB ADF_GO_S2Cl2 STARTED
[19.10|12:19:29] JOB ADF_GO_S2Cl2 RUNNING
... output trimmed ....
  3.4107,   0.0114
  3.5386,   0.0160
  3.5400,   0.0011
  3.9864,   0.1105
  4.3225,   0.0049
  4.3513,   0.2551
  4.7544,   0.0011
  4.9414,   0.0105
  5.3188,   0.0036
  5.3272,   0.0721

Python Script¶

#!/usr/bin/env python
# coding: utf-8

# ## Initial Imports

from pathlib import Path
from scm.base import ChemicalSystem, Units
from scm.inputs import ADF, AMS, DFTB
from scm.plams import AMSResults, AMSJob


# ## Helper Functions
# Set up a couple of useful functions for extracting results.

def get_excitations(results):
    """Returns excitation energies (in eV) and oscillator strengths (in Debye)."""
    if results.job.ok():
        exci_energies_au = results.readrkf("Excitations SS A", "excenergies", file="engine")
        oscillator_str_au = results.readrkf(
            "Excitations SS A", "oscillator strengths", file="engine"
        )
        # The results are stored in atomic units. Convert them to more convenient units:
        au_to_ev = Units.conversion_factor("au", "eV")
        au_to_debye = Units.conversion_factor("au", "Debye")
        exci_energies = [energy * au_to_ev for energy in exci_energies_au]
        oscillator_str = [strength * au_to_debye for strength in oscillator_str_au]
        return exci_energies, oscillator_str
    else:
        return [], []


def has_good_excitations(results, min_energy, max_energy, oscillator_str_threshold=1e-4):
    """Returns True if there is at least one excitation with non-vanishing oscillator strength
    in the energy range [min_energy, max_energy]. Unit for min_energy and max energy: eV."""
    exci_energies, oscillator_str = get_excitations(results)
    for e, o in zip(exci_energies, oscillator_str):
        if min_energy < e < max_energy and o > oscillator_str_threshold:
            return True
    return False


# ## Calculation settings
# 
# Configure the settings for the various jobs.

# Settings for geometry optimization with the AMS driver:
go_sett = AMS(Task="GeometryOptimization")
go_sett.GeometryOptimization.Convergence.Gradients = 1.0e-4


# Settings for single point calculation with the AMS driver
sp_sett = AMS(Task="SinglePoint")


# Settings for the DFTB engine (including excitations)
dftb_sett = AMS()
dftb_sett.Engine = DFTB(Model="SCC-DFTB", ResourcesDir="QUASINANO2015")
dftb_sett.Engine.Properties.Excitations.TDDFTB.Calc = "singlet"
dftb_sett.Engine.Properties.Excitations.TDDFTB.Lowest = 10
dftb_sett.Engine.Occupation.Temperature = 5.0


# Settings for the geometry optimization with the ADF engine
adf_sett = AMS()
adf_sett.Engine = ADF()
adf_sett.Engine.Basis.Type = "DZP"
adf_sett.Engine.NumericalQuality = "Basic"


# Settings for the excitation calculation using the ADF engine
adf_exci_sett = AMS()
adf_exci_sett.Engine = ADF()
adf_exci_sett.Engine.Basis.Type = "TZP"
adf_exci_sett.Engine.XC.GGA = "PBE"
adf_exci_sett.Engine.NumericalQuality = "Basic"
adf_exci_sett.Engine.Symmetry = "NoSym"
adf_exci_sett.Engine.Excitations.Lowest = [10]
adf_exci_sett.Engine.Excitations.OnlySing = True


# ## Load Molecules
# Import all xyz files in the folder 'molecules'.

from pathlib import Path

molecules_dir = Path("molecules")
molecules_dir.mkdir(exist_ok=True)

alh3_xyz = molecules_dir / "AlH3.xyz"
cscl2_xyz = molecules_dir / "CSCl2.xyz"
h20_xyz = molecules_dir / "H2O.xyz"
nh3_xyz = molecules_dir / "NH3.xyz"
s2cl2_xyz = molecules_dir / "S2Cl2.xyz"

alh3_xyz.write_text("""\
4

Al  0.0  0.0  0.0
H   1.0  0.0  0.0
H   0.5  1.0  0.0
H   0.5 -1.0  0.0
""")

cscl2_xyz.write_text("""\
4

C   0.0  0.0  0.0
S   0.0  0.0 -1.3
Cl  0.0  1.0  1.0
Cl  0.0 -1.0  1.0
""")

h20_xyz.write_text("""\
3

O 0.0 0.0 0.0
H 1.0 0.0 0.0
H 0.0 1.0 0.0
""")

nh3_xyz.write_text("""\
4

N 0.0 0.0 0.1
H 1.0 0.0 0.0
H 0.5 1.0 0.0
H 0.5 -1.0 0.0
""")

s2cl2_xyz.write_text("""\
4

S   0.0  0.0  1.3
S   0.0  0.0 -1.3
Cl  0.1  0.6  2.3
Cl  0.1 -0.6 -2.3
""")

molecules_dir = Path("molecules")
molecules = {path.stem: ChemicalSystem.from_xyz(str(path)) for path in sorted(molecules_dir.glob("*.xyz"))}


# ## DFTB Prescreen
# Perform an initial prescreen of all molecules with DFTB.

promising_molecules = {}


for name, mol in molecules.items():
    dftb_job = AMSJob(name="DFTB_" + name, molecule=mol, settings=go_sett + dftb_sett)
    dftb_results = dftb_job.run()

    if has_good_excitations(dftb_results, 1, 6):
        promising_molecules[name] = dftb_results.get_main_system()


print(f"Found {len(promising_molecules)} promising molecules with DFTB")


# ## Optimization and excitations calculation with ADF
# For each of the molecules identified in the prescreen, run a further calculation with ADF.

for name, mol in promising_molecules.items():
    adf_go_job = AMSJob(name="ADF_GO_" + name, molecule=mol, settings=go_sett + adf_sett)
    adf_go_job.run()

    optimized_mol = adf_go_job.results.get_main_system()

    adf_exci_job = AMSJob(
        name="ADF_exci_" + name, molecule=optimized_mol, settings=sp_sett + adf_exci_sett
    )
    adf_exci_results = adf_exci_job.run()

    if has_good_excitations(adf_exci_results, 2, 4):
        print(f"Molecule {name} has excitation(s) satisfying our criteria!")
        print(optimized_mol)
        exci_energies, oscillator_str = get_excitations(adf_exci_results)
        print("Excitation energy [eV], oscillator strength:")
        for e, o in zip(exci_energies, oscillator_str):
            print(f"{e:8.4f}, {o:8.4f}")

Screening Molecules for Promising Electronic Excitations¶

Initial Imports¶

Helper Functions¶

Calculation settings¶

Load Molecules¶

DFTB Prescreen¶

Optimization and excitations calculation with ADF¶

See also¶

Related examples¶

Related tutorials¶

Related documentation¶

Python Script¶