双光子相关性质的计算

遇到了TPA计算的需求，稍微了解了一下发现好像不难算，只是没现成工具。这里记录一下摸索过程

双光子吸收截面

分两种，一种是垂直的，这个比较好做；还有一种是考虑FCHT，这个得计算激发态的优化然后做振动态求和。

垂直

首先要做态求和拿到二光子跃迁张量，Multiwfn可以做这一步（18-5-xx.log-3）。

然后基于导出的SOS.txt，自己编个脚本来做计算。比如：

#!/usr/bin/env python3
"""Vertical two-photon absorption from a Multiwfn SOS.txt file.

This script reads the SOS.txt generated by Multiwfn's excited-state transition dipole
module, computes the state-resolved two-photon transition tensor S_ab(0->f; omega),
writes FCclasses3-compatible derip_form, and optionally estimates a vertical TPA
cross section in GM using a normalized Lorentzian final-state line shape.

SOS formula, atomic units:
  S_ab = sum_k [mu_0k^a mu_kf^b + mu_0k^b mu_kf^a] / (E_k - omega)

Raw TPA invariants follow the FCclasses3 convention:
  D = sum_ab [F*S_aa*S_bb + G*S_ab*S_ab + H*S_ab*S_ba]
  linear:   F=G=H=2
  circular: F=-2, G=H=3

GM estimate:
  delta_TPA(GM) = 8.35158e-4 * omega^2 * g(2*omega, E_f, Gamma) * D
  where omega, E_f, Gamma and g are in atomic units and
  g = Gamma / [pi*((2*omega - E_f)^2 + Gamma^2)].
"""
from __future__ import annotations

import argparse
from pathlib import Path
import sys
import numpy as np

AU2EV = 27.211386245988
EV_NM = 1239.8419843320026
GM_PREFACTOR = 8.35158e-4


def read_sos(path: Path) -> tuple[np.ndarray, np.ndarray]:
    lines: list[str] = []
    for raw in path.read_text(errors="replace").splitlines():
        line = raw.split("//", 1)[0].strip()
        if line:
            lines.append(line)
    if not lines:
        raise ValueError("empty SOS file")
    nstates_raw = int(float(lines[0].split()[0]))
    if nstates_raw < 0:
        raise ValueError("This SOS.txt was generated for alpha-only calculations; excited-state transition dipoles are absent")
    nstates = nstates_raw
    if len(lines) < 1 + nstates:
        raise ValueError("SOS file is truncated before the excitation-energy block")

    energy_au = np.zeros(nstates + 1, dtype=float)
    for n in range(1, nstates + 1):
        toks = lines[n].split()
        istate = int(float(toks[0]))
        if not (1 <= istate <= nstates):
            raise ValueError(f"state index out of range in energy line: {lines[n]!r}")
        energy_au[istate] = float(toks[1]) / AU2EV

    mu = np.full((nstates + 1, nstates + 1, 3), np.nan, dtype=float)
    for line in lines[1 + nstates:]:
        toks = line.split()
        if len(toks) < 5:
            continue
        i = int(float(toks[0])); j = int(float(toks[1]))
        if not (0 <= i <= nstates and 0 <= j <= nstates):
            raise ValueError(f"state index out of range in dipole line: {line!r}")
        vec = np.array([float(toks[2]), float(toks[3]), float(toks[4])], dtype=float)
        mu[i, j, :] = vec
        mu[j, i, :] = vec

    missing = np.argwhere(np.isnan(mu[:, :, 0]))
    if missing.size:
        sample = ", ".join(f"({i},{j})" for i, j in missing[:10])
        raise ValueError(f"missing transition/permanent dipoles for state pairs, e.g. {sample}")
    return energy_au, mu


def omega_from_args(args: argparse.Namespace, ef_au: float) -> float:
    supplied = [args.omega_au is not None, args.omega_ev is not None, args.omega_nm is not None]
    if sum(supplied) > 1:
        raise ValueError("use only one of --omega-au, --omega-ev, or --omega-nm")
    if args.omega_au is not None:
        return float(args.omega_au)
    if args.omega_ev is not None:
        return float(args.omega_ev) / AU2EV
    if args.omega_nm is not None:
        return (EV_NM / float(args.omega_nm)) / AU2EV
    return ef_au / 2.0


def compute_tensor(energy_au: np.ndarray, mu: np.ndarray, final_state: int, omega_au: float, ninter: int) -> tuple[np.ndarray, list[tuple[int,float,np.ndarray]]]:
    nstates = len(energy_au) - 1
    if not (1 <= final_state <= nstates):
        raise ValueError(f"final state must be between 1 and {nstates}")
    if ninter <= 0 or ninter > nstates:
        ninter = nstates
    if ninter < final_state:
        ninter = final_state

    S = np.zeros((3, 3), dtype=float)
    contribs: list[tuple[int,float,np.ndarray]] = []
    for k in range(ninter + 1):
        denom = energy_au[k] - omega_au
        if abs(denom) < 1.0e-12:
            raise ZeroDivisionError(
                f"intermediate state k={k} has near-zero denominator E_k-omega={denom:.3e} a.u.; "
                "choose a different omega or implement a complex damping convention"
            )
        C = np.zeros((3, 3), dtype=float)
        for a in range(3):
            for b in range(3):
                C[a, b] = (mu[0, k, a] * mu[k, final_state, b] + mu[0, k, b] * mu[k, final_state, a]) / denom
        S += C
        contribs.append((k, denom, C))
    return S, contribs


def tpa_invariants(S: np.ndarray) -> tuple[float, float]:
    trterm = 0.0
    normterm = 0.0
    exchterm = 0.0
    for a in range(3):
        for b in range(3):
            trterm += S[a, a] * S[b, b]
            normterm += S[a, b] * S[a, b]
            exchterm += S[a, b] * S[b, a]
    return 2.0 * trterm + 2.0 * normterm + 2.0 * exchterm, -2.0 * trterm + 3.0 * normterm + 3.0 * exchterm


def main() -> int:
    ap = argparse.ArgumentParser(description="Compute vertical TPA tensor and optional GM estimate from Multiwfn SOS.txt")
    ap.add_argument("sos", type=Path, help="Multiwfn SOS.txt generated with excited-state transition dipoles")
    ap.add_argument("--state", type=int, required=True, help="target final excited state index f")
    ap.add_argument("--ninter", type=int, default=0, help="number of intermediate excited states; 0 = all")
    ap.add_argument("--omega-au", type=float, help="incident photon frequency in a.u.; default = E_f/2")
    ap.add_argument("--omega-ev", type=float, help="incident photon energy in eV")
    ap.add_argument("--omega-nm", type=float, help="incident wavelength in nm")
    ap.add_argument("--hwhm-ev", type=float, default=0.0, help="Lorentzian HWHM Gamma in eV for GM estimate; 0 = skip GM")
    ap.add_argument("--out", type=Path, default=Path("derip_form"), help="output FCclasses3 derip_form")
    ap.add_argument("--info", type=Path, default=Path("tpa_vertical_result.txt"), help="diagnostic result file")
    ap.add_argument("--contrib", type=Path, default=Path("tpa_S_contrib.txt"), help="per-intermediate-state contribution file")
    args = ap.parse_args()

    energy_au, mu = read_sos(args.sos)
    nstates = len(energy_au) - 1
    if not (1 <= args.state <= nstates):
        raise ValueError(f"--state must be between 1 and {nstates}")
    omega_au = omega_from_args(args, energy_au[args.state])
    ninter = args.ninter if args.ninter > 0 else nstates
    if ninter < args.state:
        ninter = args.state
    if ninter > nstates:
        ninter = nstates

    S, contribs = compute_tensor(energy_au, mu, args.state, omega_au, ninter)
    dlin, dcir = tpa_invariants(S)

    gtp = None
    gm_lin = None
    gm_cir = None
    if args.hwhm_ev and args.hwhm_ev > 0:
        gamma_au = args.hwhm_ev / AU2EV
        gtp = gamma_au / (np.pi * ((2.0 * omega_au - energy_au[args.state]) ** 2 + gamma_au ** 2))
        gm_lin = GM_PREFACTOR * omega_au**2 * gtp * dlin
        gm_cir = GM_PREFACTOR * omega_au**2 * gtp * dcir

    np.savetxt(args.out, S, fmt="%20.10E")
    with args.contrib.open("w") as fh:
        fh.write("# k E_k(a.u.) E_k(eV) denom(a.u.) dSxx dSxy dSxz dSyx dSyy dSyz dSzx dSzy dSzz\n")
        for k, denom, C in contribs:
            vals = [k, energy_au[k], energy_au[k] * AU2EV, denom, *C.reshape(-1)]
            fh.write("{:6d} ".format(vals[0]) + " ".join(f"{x:18.8E}" for x in vals[1:]) + "\n")

    with args.info.open("w") as fh:
        fh.write(f"final_state {args.state}\n")
        fh.write(f"n_intermediate_excited_states {ninter}\n")
        fh.write(f"excitation_energy_au {energy_au[args.state]:.12E}\n")
        fh.write(f"excitation_energy_eV {energy_au[args.state] * AU2EV:.12E}\n")
        fh.write(f"omega_au {omega_au:.12E}\n")
        fh.write(f"omega_eV {omega_au * AU2EV:.12E}\n")
        fh.write(f"omega_nm {EV_NM / (omega_au * AU2EV):.12E}\n")
        fh.write("S_tensor_au\n")
        np.savetxt(fh, S, fmt="%20.10E")
        fh.write(f"D_TPA_linear {dlin:.12E}\n")
        fh.write(f"D_TPA_circular {dcir:.12E}\n")
        if gtp is not None:
            fh.write(f"hwhm_eV {args.hwhm_ev:.12E}\n")
            fh.write(f"hwhm_au {args.hwhm_ev / AU2EV:.12E}\n")
            fh.write(f"lorentzian_g_au_inv {gtp:.12E}\n")
            fh.write(f"delta_TPA_linear_GM {gm_lin:.12E}\n")
            fh.write(f"delta_TPA_circular_GM {gm_cir:.12E}\n")

    print(f"Wrote {args.out}, {args.info}, and {args.contrib}")
    print(f"state={args.state}, E_f={energy_au[args.state] * AU2EV:.8f} eV, omega={omega_au * AU2EV:.8f} eV, ninter={ninter}")
    print(f"D_TPA linear={dlin:.8E}, circular={dcir:.8E}")
    if gm_lin is not None:
        print(f"delta_TPA linear={gm_lin:.8E} GM, circular={gm_cir:.8E} GM")
    return 0


if __name__ == "__main__":
    try:
        raise SystemExit(main())
    except Exception as exc:
        print(f"ERROR: {exc}", file=sys.stderr)
        raise SystemExit(1)

Ref:

• 使用Multiwfn基于完全态求和(SOS)方法计算极化率和超极化率
• 求助：Gaussiain计算双光子吸收截面的问题

考虑FC近似

这里要做振动态求和还是用专业软件比较好。以FCclasses为例，需要额外准备一个二光子跃迁张量文件derip_form，其格式为:

Sxx  Sxy  Sxz
Syx  Syy  Syz
Szx  Szy  Szz

可以从Multiwfn导出的SOS.txt里整理得到这些数据。随后利用FCclasses计算TPA：

$$$
PROPERTY = TPA
MODEL = AH
DIPOLE = FC
TEMP = 298.15
BROADFUN = GAU
HWHM = 0.10
METHOD = TD

NORMALMODES = COMPUTE
COORDS = CARTESIAN

STATE1_FILE = S0.fcc
STATE2_FILE = S1.fcc