20.7. DAISpY: A Domain Assignment and Interface Solution in pYthon for Charge-Transfer Analysis - Charge Transfer Analysis Property

The Domain Assignment and Interface Solution in pYthon (DAISpY) module provides domain-based charge-transfer (CT) analysis of electronically excited states within PyBEST. Given an excited-state CI expansion from an EOM-CC or related calculation, DAISpY constructs a domain-to-domain CT matrix that quantifies how electron density is redistributed between user-defined molecular fragments upon excitation.

The module is designed to be method-agnostic: it operates on CI-like amplitudes and molecular orbital coefficients produced by any of the excited-state solvers available in PyBEST (see Quick Guide), including EOM-pCCD+S and EOM-CCSD.

For more details and if you use this module, please cite [szczuczko2026a] and [szczuczko2026b].

20.7.1. Background and Theory 

Understanding the spatial character of electronic excitations—in particular, the distinction between local and charge-transfer contributions—is essential for the analysis of excited-state phenomena in organic electronics, light-harvesting devices, and energy-conversion technologies. DAISpY formalizes this analysis by introducing the possibility to specify chemically motivated atom-based domains and resolving excitation amplitudes with respect to those domains.

Note

The current PyBEST-integrated DAISpY implementation supports spin-free excitation vectors. Spin-resolved excitation analysis is outside the scope of the present implementation.

The overall workflow is organized into four layers:

Domain definition — atoms are grouped into labeled spatial domains (e.g., Donor, Acceptor, Bridge).
Orbital-to-domain assignment — each active-space MO is assigned to a domain based on its AO composition.
CT matrix construction — CI excitation weights are accumulated into a domain-to-domain matrix.
Derived descriptors — scalar quantities (total CT fraction, intra-/inter-domain ratios) are extracted.

20.7.1.1. Physical Input Quantities 

The CT analysis requires the following quantities from a preceding electronic-structure calculation:

Quantity	Symbol	Physical meaning
AO overlap matrix	\(\mathbf{S}\)	Overlap of non-orthogonal atomic orbitals
MO coefficient matrix	\(\mathbf{C}\)	Expansion of molecular orbitals in the AO basis
CI expansion coefficients	\(\mathbf{c}^{(k)}\)	Excited-state wave function amplitudes
Excitation energies (optional)	\(E^{(k)}\)	Electronic energies of excited states
Atomic coordinates	\(\{\mathbf{R}_A\}\)	Cartesian nuclear positions
Shell–atom mapping	—	Assignment of AO shells to atomic centres
Shell type identifiers	—	Angular-momentum character of basis functions
Active-space metadata	—	\(n_\mathrm{core}\), \(n_\mathrm{acto}\), \(n_\mathrm{actv}\)

In a typical PyBEST workflow, all of these quantities are provided automatically by the IOData container returned by the excited-state solver (e.g., REOMpCCDS or REOMCCSD).

20.7.1.2. Domain Definition 

Each atom in the molecular system is assigned to exactly one domain, inducing a partition of the molecule into chemically meaningful regions. Typical choices include molecular fragments, donor/acceptor units, or individual chromophores.

Domain labels are provided in one of the following ways:

Augmented XYZ file — a standard XYZ file where each atom line carries an additional last column containing the domain label:
```
5
comment line
C   0.0  0.0  0.0  Donor
H   0.0  0.0  1.0  Donor
O   1.2  0.0  0.0  Acceptor
H   1.8  0.7  0.0  Acceptor
H   1.8 -0.7  0.0  Acceptor
```
This path is recommended for headless or batch runs.
Python dictionary — a mapping of 0-based atom indices to domain label strings, passed directly via the domains option.

For example:
```
property_options={
    "domains": {
        0: "Donor",
        1: "Donor",
        2: "Acceptor",
        3: "Acceptor",
        4: "Acceptor",
    },
}
```
The atom indices are 0-based and must follow the atom ordering used in the molecular geometry.
GUI picker — if labels are missing or incomplete, the DAISpY graphical domain picker (PySide6 + 3Dmol.js) opens automatically. The picker supports click-to-assign, group selection, right-click unassign, domain rename/recolor, and enforces that all atoms are assigned before confirmation.

No restrictions are imposed on the number, size, or spatial extent of the domains.

For complete runnable examples showing how domain specifications are used in practice, see Example Usage.

20.7.1.3. Orbital-to-Domain Assignment 

Once domains have been defined, each active-space molecular orbital is analyzed in terms of its atomic-orbital composition. Since atomic orbitals are not orthogonal in general, DAISpY first transforms the MO coefficients to an orthonormal AO representation,

\[C_{\mu p}^{\mathrm{ortho}} = \sum_{\nu} S_{\mu\nu}^{1/2} C_{\nu p},\]

where \(S_{\mu\nu}\) is the AO overlap matrix and \(C_{\nu p}\) is the coefficient of AO \(\nu\) in MO \(p\).

For a given orbital \(p\) and domain \(D\), DAISpY defines a scheme-dependent domain factor,

\[\begin{split}\Omega_D(p) = \begin{cases} 1, & \text{if } D = \arg\max_{D'} \sum_{\mu \in D'} \left| C_{\mu p}^{\mathrm{ortho}} \right|^2, \\[4pt] 0, & \text{otherwise} \end{cases} \quad \text{(hard scheme)}\end{split}\]

\[\Omega_D(p) = \sum_{\mu \in D} \left| C_{\mu p}^{\mathrm{ortho}} \right|^2 \quad \text{(weighted scheme)}.\]

In the hard scheme, \(\Omega_D(p)\) is either 0 or 1 and indicates whether orbital \(p\) is assigned to domain \(D\). In the weighted scheme, \(\Omega_D(p)\) is a continuous value between 0 and 1 and measures how strongly orbital \(p\) is localized on domain \(D\).

Equivalently, DAISpY first computes the continuous MO-domain weights,

\[w_p^D = \sum_{\mu \in D} \left| C_{\mu p}^{\mathrm{ortho}} \right|^2.\]

For ct_matrix_mode="hard", these weights are converted into a dominant-domain assignment,

\[D(p) = \arg\max_D w_p^D.\]

This dominant-domain assignment is controlled by two user-configurable thresholds:

major_thresh (default 0.10) — the minimum weight a domain must contribute to be considered.
tie_break (default 0.50) — when multiple domains exceed major_thresh, the winning domain must reach at least this weight to resolve the tie.

If strict=True, the analysis raises an error for ambiguous assignments; otherwise the maximum-weight domain is chosen and a confidence value is stored.

For ct_matrix_mode="weighted", the continuous MO-domain weights are used directly in the CT matrix accumulation, so the hard-assignment thresholds do not change the weighted CT values.

Note

The major_thresh, tie_break, and strict options affect the hard dominant-domain assignment. They are relevant for ct_matrix_mode="hard" and for reporting the selected MO-domain labels and confidence values. In ct_matrix_mode="weighted", the CT matrix is accumulated from the continuous MO-domain weights directly.

20.7.1.4. CT Matrix Construction 

The central output of the DAISpY analysis is the domain-to-domain CT matrix \(\mathbf{M}\), constructed separately for each excited state. Each element \(M_{D_h D_p}\) quantifies the contribution of excitations moving electron density from hole domain \(D_h\) to particle domain \(D_p\).

Diagonal elements correspond to local (intra-domain) excitations.
Off-diagonal elements correspond to inter-domain charge transfer.

20.7.1.4.1. Hard (discrete) accumulation

20.7.1.4.1.1. Singles

For single excitations, the CT matrix element is accumulated as

\[M_{D_h D_p}^{(S)} = \sum_{i,a} \left| c_i^a \right|^2 \delta_{D(i),D_h}\, \delta_{D(a),D_p},\]

where \(c_i^a\) is the amplitude for the single excitation \(i \to a\), and \(D(i)\), \(D(a)\) are the domains assigned to the occupied and virtual orbitals.

Here, \(\delta\) denotes the Kronecker delta. It acts as a selector: \(\delta_{D(i),D_h}=1\) only when the domain assigned to orbital \(i\) is the same as the hole-domain index \(D_h\); otherwise it is zero. Similarly, \(\delta_{D(a),D_p}\) selects excitations whose virtual orbital \(a\) belongs to the particle domain \(D_p\).

Therefore, the product of the two delta terms ensures that a given excitation contributes only to the CT matrix element corresponding to its assigned donor–acceptor domain pair.

20.7.1.4.1.2. Doubles

Double excitations \((i,j) \to (a,b)\) do not define a unique donor–acceptor pair. DAISpY decomposes each double into two effective hole–particle channels \(i \to a\) and \(j \to b\), each contributing half the total weight:

\[M_{D_h D_p}^{(D)} = \sum_{i,j,a,b} \frac{1}{2}\left| c_{ij}^{ab} \right|^2 \left[ \delta_{D(i),D_h}\,\delta_{D(a),D_p} + \delta_{D(j),D_h}\,\delta_{D(b),D_p} \right].\]

This preserves normalization and avoids double counting.

20.7.1.4.1.3. Pairs (electron-pair excitations)

For the special case of pair excitations \((i,\bar{i}) \to (a,\bar{a})\), both spin components map to the same spatial orbital domain, so the pair contribution reduces to

\[M_{D_h D_p}^{(P)} = \sum_{i,a} \left| c_{i\bar{i}}^{a\bar{a}} \right|^2 \delta_{D(i),D_h}\,\delta_{D(a),D_p},\]

which is consistent with the spin-free accumulation rule.

20.7.1.4.2. Weighted accumulation

As an alternative, DAISpY supports a weighted accumulation scheme that distributes excitation weights according to the continuous orbital-domain weights rather than discrete assignments. For singles:

\[M_{D_h D_p}^{(w)} = \sum_{i,a} \left| c_i^a \right|^2 \, w_{iD_h} \, w_{aD_p},\]

and for doubles:

\[M_{D_h D_p}^{(w,D)} = \sum_{i,j,a,b} \frac{\left| c_{ij}^{ab} \right|^2}{2} \left( w_{iD_h}\,w_{aD_p} + w_{jD_h}\,w_{bD_p} \right).\]

The weighted scheme reflects the underlying orbital delocalization and can be especially useful for small molecules or many small domains. As domain sizes increase, the two schemes converge.

The accumulation mode is controlled by the ct_matrix_mode option:

"hard" — discrete dominant-domain assignment (Eqs. \(M^{(S)}, M^{(D)}, M^{(P)}\))
"weighted" — continuous weight-based accumulation (Eqs. \(M^{(w)}, M^{(w,D)}\))

The default is "weighted".

20.7.1.4.3. Excitation Indexing

Internally, single excitations span the dense occupied–virtual product space of dimension \(N_S = N_\mathrm{occ} \times N_\mathrm{virt}\). A composite index \(k\) is decoded as

\[a = k \bmod N_\mathrm{virt}, \qquad i = \left\lfloor \frac{k}{N_\mathrm{virt}} \right\rfloor.\]

Double excitations are stored as unordered pairs \((k_1, k_2)\) with \(k_1 \le k_2\), yielding \(N_D = N_S(N_S + 1)/2\) configurations.

The CI vector layout (after optional ground-state removal) follows one of these patterns:

\(N_S\) — singles only (S)
\(2 N_S\) — singles + pairs (S+P)
\(N_D\) — doubles only (D)
\(N_S + N_D\) — singles + doubles (S+D)

DAISpY automatically detects the layout from the CI vector length and the active-space sizes.

20.7.1.5. Derived Descriptors 

From the CT matrix, several scalar quantities can be derived:

Total CT fraction — the sum of off-diagonal elements, indicating inter-domain electron redistribution.
Local excitation fraction — the sum of diagonal elements (1 − total CT).
Dominant donor–acceptor pair — the largest off-diagonal entry.
Intra-/inter-domain ratio — particularly useful when doubles are included.

These descriptors are printed automatically at the selected verbosity level and stored in the checkpoint under ct/* keys.

20.7.2. Implemented Features 

The CT analysis is accessible through the CTProperty class, which follows the standard PyBEST PropertyBase lifecycle. It supports:

Singles (S), pairs (P), and doubles (D) excitation blocks, individually or combined
Hard and weighted CT matrix accumulation modes
Single-state and multi-state analysis (single index, comma lists, ranges, or "all")
Augmented XYZ, Python dict, and GUI domain specification paths
Input directory mode for running without a full PyBEST checkpoint
Three verbosity levels (0 = tabular/spreadsheet, 1 = basic report, 2 = extended report)

20.7.3. Quick Guide: CT Analysis with EOM-pCCD+S 

This section provides a concise introduction to running a CT analysis. For a complete description of all options, see the Keyword Arguments section below.

We assume the following standard PyBEST objects are available (see Preliminaries and Quick Guide for details):

lf:

a LinalgFactory instance

occ_model:

an Aufbau occupation model of the AufbauOccModel class

The basic workflow is:

Run the ground-state and excited-state calculations.
Prepare a domain specification (e.g., augmented XYZ file).
Construct a CTProperty and call it with the excited-state output.

from pybest.properties.ct import CTProperty

# Assume eom_output is an IOData container from a previous EOM calculation
# and domains_xyz points to a labeled XYZ file

ct = CTProperty(lf, occ_model)
ct_output = ct(
    eom_output,
    property_options={
        "domains": "domains.xyz",
        "excitations": ("SP",),   # singles + pairs
    },
)

The CTProperty instance handles the entire pipeline automatically:

Reads and validates inputs (CI vectors, overlap, orbital coefficients).
Computes MO-to-domain weights and assigns MOs to domains.
Aggregates the CT matrix and prints results.

The results are returned as an IOData container and saved to the checkpoint.

20.7.4. Domain Specification Options 

DAISpY supports several ways to specify domains.

20.7.4.1. Augmented XYZ file (recommended for batch/headless)

Provide a path to an XYZ file where each atom line has an additional last column with the domain label. This is the recommended approach for non-interactive environments:

property_options={
    "domains": "path/to/domains.xyz",
}

20.7.4.2. Python dictionary 

Pass a mapping of 0-based atom indices to labels:

property_options={
    "domains": {0: "Donor", 1: "Donor", 2: "Acceptor", 3: "Acceptor", 4: "Acceptor"},
}

20.7.4.3. GUI picker (interactive)

If domains is omitted or points to an unlabeled XYZ, the graphical domain picker opens automatically. It supports:

Click-to-assign atoms to the current domain
Right-click to unassign
Rectangle and lasso group selection
Add/delete/rename domains with custom colors
Confirmation only enabled when all atoms are assigned

After confirmation, the domain mapping can optionally be exported to an augmented XYZ file via the export_domains option:

property_options={
    "domains": "geometry.xyz",         # unlabeled → GUI opens
    "export_domains": "labeled.xyz",   # save the result
}

20.7.5. Excitation Block Selection 

The excitations option controls which CI excitation blocks contribute to the CT matrix:

"S" — singles
"P" — pairs (electron-pair excitations)
"D" — doubles (general two-electron excitations)

Multiple blocks can be combined:

# Any of these forms are equivalent:
"excitations": ("S",)          # singles only
"excitations": ("SP",)         # singles + pairs
"excitations": ("S", "D")      # singles + doubles
"excitations": ("S", "P", "D") # all three

The parser normalizes input by uppercasing and splitting combined strings, so "SP", ("S", "P"), and "S,P" are all equivalent.

If a requested block is not present in the CI data, a warning is printed and that block is skipped. If none of the requested blocks are available, the analysis raises an error.

Note

When both pairs (P) and doubles (D) are requested, pairs are implicitly included in the doubles decomposition. DAISpY will hide P to avoid double counting and print a notification.

20.7.6. State Selection 

The state_index option selects which excited states to analyze. Internally, CT uses 1-based indexing for excited states (0 is reserved for the ground state in the CI layout):

# Single state (default)
"state_index": 1

# Comma-separated list
"state_index": "1,3,7"

# Range
"state_index": "2-5"

# All excited states
"state_index": "all"

If some requested states are out of range, they are silently dropped with a warning. If none of the requested states are valid, an error is raised.

20.7.7. Keyword Arguments 

All options are passed as a dictionary to the property_options argument. The full set of supported keys is:

20.7.7.1. Domain and Input Options 

domains:

Domain specification source. Accepts a path to an augmented XYZ file, a Python dict mapping atom indices to labels, a DomainSpec instance, or None (triggers GUI if geometry is available).

export_domains:

(str or None) Path to export the augmented XYZ file with domain labels after GUI selection. Ignored if domains were already fully specified.

iodata:

(IOData) Explicit excited-state data container. When using the standard call signature ct(eom_output, ...), this is set automatically.

input_dir:

(str or None) Path to a user-prepared input directory containing meta.json and array files (see User Input Directory Mode).

20.7.7.2. Analysis Options 

state_index:

(int, str) Excited state(s) to analyze. Accepts a single index, comma-separated list, range (e.g., "2-5"), or "all" (default: 1).

excitations:

(tuple of str) Excitation blocks to include: "S" (singles), "P" (pairs), "D" (doubles). Combined forms like "SP" or "SD" are accepted (default: ("S",)).

ct_matrix_mode:

(str) Accumulation scheme: "hard" for discrete orbital-domain assignment, "weighted" for continuous weight-based accumulation (default: "weighted").

20.7.7.3. MO-Domain Assignment Options 

major_thresh:

(float) Minimum absolute weight required for a domain to be considered as an assignment candidate (default: 0.10).

tie_break:

(float) Minimum absolute weight required for a winning domain when multiple candidates exceed major_thresh (default: 0.50).

strict:

(bool) If True, raise an error for ambiguous MO assignments instead of silently choosing the dominant domain (default: False).

Note

These options affect the hard MO-domain assignment step. They control how the continuous MO-domain weights are converted into a single domain label per active-space MO. They do not define an independent weighted CT accumulation scheme in the current PyBEST implementation.

20.7.7.4. Reporting Options 

verbose:

(int) Verbosity level for printed output:

0 — simplified tabular output (spreadsheet-friendly)

1 — basic report with per-state CT matrix and excitation-block sums (default)

2 — extended report including doubles-only matrix and intra-/inter-domain summary

20.7.8. Output and Results 

The CTProperty stores its results in the checkpoint under ct/* keys. The most important entries are:

ct/matrix:

The domain-to-domain CT matrix for the last analyzed state.

ct/domain_labels:

The ordered list of domain labels defining the matrix axes.

ct/domains:

The full domain mapping (labels + atom-to-domain dict).

ct/mo_weights:

Per-active-MO weights over domains (shape: n_active_MOs × n_domains).

ct/mo_labels:

The chosen domain index per active-space MO.

ct/mo_confidence:

Confidence (weight) of each MO-domain assignment.

ct/used_excitations:

Excitation blocks actually used in the aggregation.

ct/e_ee:

Excitation energies (if available).

ct/matrix_mode:

The accumulation mode used ("hard" or "weighted").

ct/assigner:

Assignment metadata: type, major, tie.

orb_a_ortho:

Orthogonalized orbital coefficients used for the MO-domain weight computation.

A typical printed report at verbose=1 looks like:

Domain-based CT values

     MO→domain assignment (active space):
       A: 28/54 (51.9%)
       B: 26/54 (48.1%)

[INFO] Included excitation blocks: singles (S)

State index: 1
     Excitation energy: 0.260015 au (7.0754 eV)

          A -> A         :     0.00 %
          A -> B         :     0.04 %
          B -> A         :     0.24 %
          B -> B         :    99.72 %

     Sum singles (S):   100.00 % (without GS contribution)

20.7.9. User Input Directory Mode 

For advanced use cases where a full PyBEST checkpoint is not available, DAISpY supports a user-prepared input directory containing the raw numerical data needed for CT analysis.

20.7.9.1. Required files 

meta.json — active-space metadata (see below)
olp.npy or olp.txt — AO overlap matrix, shape (nbasis, nbasis)
orb_a.npy or orb_a.txt — AO-MO coefficient matrix, shape (nbasis, nmo) where nmo = ncore + nacto + nactv
civ_ee.npy / civ_ee.txt or civ_r_ee.npy / civ_r_ee.txt — CI amplitudes

20.7.9.2. Optional files 

e_ee.npy or e_ee.txt — excitation energies (used only for printing)

20.7.9.3. Required `meta.json` keys 

{
    "ncore": 1,
    "nacto": 5,
    "nactv": 19,
    "atom": [8, 1, 1],
    "coordinates": [[0.0, 0.0, 0.0], [0.0, 0.7, 0.5], [0.0, -0.7, 0.5]],
    "shell2atom": [0, 0, 0, 0, 0, 0, 1, 1, 1, 2, 2, 2],
    "shell_types": [0, 0, 0, 1, 1, 2, 0, 0, 1, 0, 0, 1]
}

Where:

“ncore”: number of frozen core orbitals (not included in the excited state calculation),
“nacto”: number of active occupied orbitals,
“nactv”: number of active virtual orbitals,
“atom”: the list of atoms as ordered in the coordinates, labeled by their atomic number/charge,
“coordinates”: the xyz coordinates in Angstrom,
“shell2atom”: a list translating the AOs to atoms; [0, 0, 0, 0, 0, 0, 1, 1, 1, 2, 2, 2] translates to [atom0, atom0, atom0, atom0, atom0, atom0, atom1, atom1, atom1, atom2, atom2, atom2]
“shell_types”: a list containing the shell type (s, p, d, etc.) of each AO listed in “shell2atom”

The shell2atom and shell_types arrays define the mapping from AO shells to atoms and are required for the AO-walk that computes MO-domain weights. Both .npy and whitespace-delimited .txt formats are supported for all array files.

To construct a CTProperty from a user input directory without an existing checkpoint:

ct = CTProperty.from_input_dir("path/to/input_dir")
ct_output = ct(
    property_options={
        "input_dir": "path/to/input_dir",
        "domains": "domains.xyz",
        "excitations": ("S",),
    },
)

20.7.10. Common Pitfalls 

AO counting mismatch — usually caused by inconsistent shell2atom/shell_types in meta.json relative to the olp/orb_a dimensions.
Wrong orb_a column count — the number of columns must equal ncore + nacto + nactv exactly; the column ordering must be [core MOs][active occupied MOs][active virtual MOs].
GUI dependency in headless environments — provide a fully labeled augmented XYZ file to avoid GUI requirements on clusters or CI systems.
Pairs + Doubles double counting — when both P and D are requested, pairs are included implicitly in the doubles decomposition. DAISpY handles this automatically.

20.7.11. Example Usage 

The example files for CT analysis are located in the directory data/examples/properties.

20.7.11.1. Example for CT analysis with EOM-pCCD+S 

This example runs a minimal PyBEST workflow for the water molecule: RHF → ROOpCCD → REOMpCCDS → CTProperty. The CT analysis uses singles and pairs excitation blocks with domain labels provided via a labeled XYZ file.

Listing 20.12 data/examples/properties/water_ct_analysis_eompccds.py

"""
Charge-Transfer (CT) analysis for water (H2O) with cc-pVDZ.

Pipeline:
RHF → ROOpCCD → REOMpCCDS (excited states) → CTProperty

CTProperty is called with the IOData container produced by the excited-state
step and a path to an XYZ file that already contains per-atom domain labels.
If labels are incomplete, the CT GUI picker will open; here we keep it
non-interactive by using a labeled XYZ.
"""

from pybest import context
from pybest.ee_eom import REOMpCCDS
from pybest.gbasis import (
    compute_eri,
    compute_kinetic,
    compute_nuclear,
    compute_nuclear_repulsion,
    compute_overlap,
    get_gobasis,
)
from pybest.geminals import ROOpCCD
from pybest.linalg import DenseLinalgFactory
from pybest.occ_model import AufbauOccModel
from pybest.properties.ct import CTProperty
from pybest.wrappers import RHF

# Molecule and basis
fn_xyz = context.get_fn("test/water.xyz")
basis = get_gobasis("cc-pvdz", fn_xyz)

# Linalg and occupation model
lf = DenseLinalgFactory(basis.nbasis)
occ_model = AufbauOccModel(basis)

# AO quantities
orb_a = lf.create_orbital(basis.nbasis)
olp = compute_overlap(basis)
kin = compute_kinetic(basis)
ne = compute_nuclear(basis)
eri = compute_eri(basis)
external = compute_nuclear_repulsion(basis)

# RHF
rhf = RHF(lf, occ_model)
hf_output = rhf(kin, ne, eri, external, olp, orb_a)

# ROOpCCD (ground state)
pccd = ROOpCCD(lf, occ_model)
pccd_output = pccd(kin, ne, eri, hf_output)

# Excited states (EOM-pCCD+S)
eompccds = REOMpCCDS(lf, occ_model)
eom_output = eompccds(kin, ne, eri, pccd_output)

# CTProperty with minimal inputs: IOData + domains xyz
domains_xyz = context.get_fn("test/water_domains.xyz")
ct = CTProperty(lf, occ_model)
ct_output = ct(
    eom_output,  # IOData-like container with excited-state vectors
    property_options={
        "domains": domains_xyz,  # xyz contains domain labels; GUI opens only if incomplete
        "excitations": (
            "SP",
        ),  # use ("S",), ("P",), ("D",), or a combination of these
    },
)

20.7.11.2. Example for CT analysis with EOM-CCSD 

This example uses an EOM-CCSD excited-state calculation: RHF → ROOpCCD → RCCSD → REOMCCSD → CTProperty. The CT analysis includes both singles and doubles excitation blocks.

Listing 20.13 data/examples/properties/water_ct_analysis_eomccsd.py

Charge-Transfer (CT) analysis for water (H2O) with cc-pVDZ.

Pipeline:
RHF → ROOpCCD → REOMpCCDS (excited states) → CTProperty

CTProperty is called with the IOData container produced by the excited-state
step and a path to an XYZ file that already contains per-atom domain labels.
If labels are incomplete, the CT GUI picker will open; here we keep it
non-interactive by using a labeled XYZ.
"""

from pybest import context
from pybest.cc import RCCSD
from pybest.ee_eom import REOMCCSD
from pybest.gbasis import (
    compute_eri,
    compute_kinetic,
    compute_nuclear,
    compute_nuclear_repulsion,
    compute_overlap,
    get_gobasis,
)
from pybest.geminals import ROOpCCD
from pybest.linalg import DenseLinalgFactory
from pybest.occ_model import AufbauOccModel
from pybest.properties.ct import CTProperty
from pybest.wrappers import RHF

# Molecule and basis
fn_xyz = context.get_fn("test/water.xyz")
basis = get_gobasis("cc-pvdz", fn_xyz)

# Linalg and occupation model
lf = DenseLinalgFactory(basis.nbasis)
occ_model = AufbauOccModel(basis)

# AO quantities
orb_a = lf.create_orbital(basis.nbasis)
olp = compute_overlap(basis)
kin = compute_kinetic(basis)
ne = compute_nuclear(basis)
eri = compute_eri(basis)
external = compute_nuclear_repulsion(basis)

# RHF
rhf = RHF(lf, occ_model)
hf_output = rhf(kin, ne, eri, external, olp, orb_a)

# ROOpCCD (ground state)
pccd = ROOpCCD(lf, occ_model)
pccd_output = pccd(kin, ne, eri, hf_output)

# RCCSD (ground state with ROOpCCD orbitals)
rccsd = RCCSD(lf, occ_model)
rccsd_output = rccsd(kin, ne, eri, pccd_output)

# Excited states (EOM-CCSD)
eomccsd = REOMCCSD(lf, occ_model)
eom_output = eomccsd(kin, ne, eri, rccsd_output, nroot=2)

# CTProperty with minimal inputs: IOData + domains xyz
domains_xyz = context.get_fn("test/water_domains.xyz")
ct = CTProperty(lf, occ_model)
ct_output = ct(
    eom_output,  # IOData-like container with excited-state vectors
    property_options={
        "domains": domains_xyz,  # xyz contains domain labels; GUI opens only if incomplete
        "excitations": (
            "SD",
        ),  # use ("S",), ("P",), ("D",), or a combination of these
    },
)

20.7.11.3. Example with multi-state analysis and all options 

This snippet demonstrates advanced usage with multiple states, custom thresholds, and the extended report:

from pybest.properties.ct import CTProperty

ct = CTProperty(lf, occ_model)
ct_output = ct(
    eom_output,
    property_options={
        "domains": "domains.xyz",
        "state_index": "1-5",         # analyze states 1 through 5
        "excitations": ("S", "D"),     # singles + doubles
        "ct_matrix_mode": "hard",      # discrete dominant-domain assignment
        "major_thresh": 0.10,          # MO assignment threshold
        "tie_break": 0.50,             # tie-break threshold
        "strict": False,               # allow ambiguous assignments
        "verbose": 2,                  # extended report
    },
)

20.7.11.4. CLI Usage Examples 

The CT analysis can also be run from the command line using the pybest-ct.py script. This is useful for post-processing existing PyBEST checkpoints, running headless analyses on clusters, exporting GUI-selected domains, or analyzing a user-prepared input directory. The CLI frontend supports checkpoint-based runs, --input-dir mode, flexible state selection, excitation-block selection, and different reporting verbosity levels.

20.7.11.4.1. Basic checkpoint-based run with labeled domains

Use this mode when a PyBEST checkpoint already contains the excited-state data and the domains are provided in an augmented XYZ file:

python scripts/pybest-ct.py \
    --checkpoint calc.h5 \
    --domains domains.xyz

20.7.11.4.2. Checkpoint run with explicit state and excitation selection

The --state-index option accepts a single state, a comma-separated list, a range, or all. The --excitations option accepts S for singles, P for pairs, and D for doubles:

python scripts/pybest-ct.py \
    --checkpoint calc.h5 \
    --domains domains.xyz \
    --state-index 1-5 \
    --excitations S D \
    --ct-matrix-mode weighted \
    --verbose 2

Equivalent compact excitation forms are also accepted:

python scripts/pybest-ct.py \
    --checkpoint calc.h5 \
    --domains domains.xyz \
    --state-index all \
    --excitations SDP

20.7.11.4.3. GUI-based domain assignment from checkpoint geometry

If --domains is omitted, DAISpY attempts to obtain the molecular geometry from the checkpoint and opens the graphical domain picker:

python scripts/pybest-ct.py \
    --checkpoint calc.h5 \
    --state-index 1 \
    --excitations S

20.7.11.4.4. Exporting domains only

To create a labeled augmented XYZ file without running the CT calculation, provide an unlabeled geometry file and an output path with --export-domains:

python scripts/pybest-ct.py \
    --domains geometry.xyz \
    --export-domains domains.xyz

This opens the GUI picker, lets the user assign domains, and writes the labeled XYZ file after confirmation.

20.7.11.4.5. Input-directory mode without a checkpoint

For workflows where the data are provided as meta.json plus array files instead of a PyBEST checkpoint, use --input-dir:

python scripts/pybest-ct.py \
    --input-dir ct_input_example \
    --domains domains.xyz \
    --state-index 1 \
    --excitations S \
    --verbose 1

20.7.11.4.6. Spreadsheet-friendly output

For compact tabular output that can be copied into a spreadsheet, use --verbose 0:

python scripts/pybest-ct.py \
    --checkpoint calc.h5 \
    --domains domains.xyz \
    --state-index all \
    --excitations S \
    --verbose 0