Polarizable Embedding “The Frame”

Warning

BETA FEATURE — The FraME model is experimental and may change or be removed in future versions without notice.

Sadly both Polarizable Embedding and Projective-Embedding are often abbreviated as PE, which can lead to confusion. Here, we will instead use Frame for Fragment-based Multiscale Embedding or written-out Polarizable Embedding, while Projective-Embedding remains PE.

The Frame acronym comes is also appropiate as you can think of the Polarizable Embedding as the outer layer - the Frame that holds the solvent/environment response, while Projective-Embedding focuses wavefunction accuracy inside your embedded Quantum Mechanical region.

In this example, we start from a PDB file of acrolein in water, then we demonstrate how to:

  • Generate a polarizable embedding potential (.pot) and DALTON molecule file (.mol) using PyFraME,

  • Convert the DALTON .mol to .xyz,

  • Finally run a Qrunch calculation that combines a polarizable “Frame” (the solvent environment) with our Projection-based Wavefunction-in-DFT embedding (PE) to solve a small active region.

This naturally yields a three-level hierarchy:

  • Frame: polarizable environment (classical, distributed multipoles + polarizabilities)

  • DFT on the full molecule: mean-field description supplying orbitals & density

  • Embedded WF solver: high-level treatment inside the active region

References and background

[Dalton, PyFraME, SRN+18, SRN+19]

Files used in this tutorial

You may preview the PDB interactively with the Dash Molecule 3D Viewer:

  1. Go to the Data tab.

  2. Click the Drag and drop or click to upload file area.

  3. Upload acrolein_water.pdb.

1) Generate the PE Frame potentials with PyFraME

Download the helper

run_pyframe.py

Given the PDB file acrolein_water.pdb in your working directory, run your helper script.

$ python run_pyframe.py

This will create a new folder acrolein_water/ with:

  • acrolein_water.pot — the polarizable embedding potential

  • acrolein_water.mol — a DALTON molecule input file

2) Convert DALTON .mol → .xyz

Qrunch expects an .xyz. So we use the converter:

Download the converter

mol_to_xyz.py

Run it

python mol_to_xyz.py acrolein_water/acrolein_water.mol acrolein_water.xyz

3) Run Qrunch with a PE Frame + Projective-Embedding

We now have the molecule input that define the quantum mechanical region .xyz and the potential .pot file that represents the polarizable environment. We can now run the qrunch script.

Full script

The script below shows the full script, which are then explained section by section.

  1"""Demo script for FAST-VQE with Polarizable Embedding and Projection based Wavefunction-in-DFT embedding."""
  2# ruff: noqa # ruff: ignore[noqa-comments]
  3
  4from pathlib import Path
  5
  6from matplotlib import pyplot as plt
  7import numpy as np
  8
  9import qrunch as qc
 10
 11
 12def plot_energies(
 13    result: dict[str, list[float]],
 14    *,
 15    outdir: Path = Path("dist"),
 16    show: bool = False,
 17) -> None:
 18    """Plot the energies from FAST-VQE and CAS-CI."""
 19    energies = np.asarray(result["energies"])
 20    reference_energy = result["casci_energy"][0]
 21
 22    iterations = np.asarray(list(range(len(energies))))
 23    # Create the plot
 24    fig1 = plt.figure()
 25    plt.plot(iterations, energies, color="darkgreen", marker="o", label="FAST-VQE")
 26
 27    # Add a dashed horizontal line at the CAS-CI value.
 28    plt.axhline(y=reference_energy, color="black", linestyle="--", label="CAS-CI")
 29
 30    plt.xlabel("Iteration number")
 31    plt.ylabel("Total energy [Hartree]")
 32    plt.title("FAST-VQE convergence")
 33    plt.legend()
 34
 35    convergence_plot = outdir / "frame_vqe_convergence.png"
 36    fig1.savefig(convergence_plot, dpi=200, bbox_inches="tight")
 37
 38    # -----------------------------
 39    # Figure 2: |Energy - Reference| (log scale)
 40    # -----------------------------
 41    # Use absolute error; add a tiny epsilon to avoid log(0) if we hit the reference exactly.
 42    eps = np.finfo(float).eps
 43    abs_error = np.abs(energies - reference_energy) + eps
 44
 45    fig2 = plt.figure()
 46    plt.semilogy(iterations, abs_error, "-*", color="darkgreen", label="|E - E_ref|")
 47    plt.xlabel("Iteration number")
 48    plt.ylabel("Absolute error [Hartree]")
 49    plt.title("FAST-VQE error vs reference")
 50    plt.legend()
 51    plt.tight_layout()
 52
 53    error_convergence_plot = outdir / "frame_vqe_error_semilogy.png"
 54    fig2.savefig(error_convergence_plot, dpi=200, bbox_inches="tight")
 55
 56    if show:
 57        plt.show()
 58    else:
 59        # Close figures to free memory when running in batch contexts.
 60        plt.close(fig1)
 61        plt.close(fig2)
 62
 63
 64
 65
 66def main(
 67    *,
 68    calculate_casci: bool = False,
 69    max_iterations: int = 200,
 70    path_to_molecule_xyz_file: Path = Path.cwd() / "acrolein_water.xyz",
 71    path_to_molecule_pot_file: Path = Path.cwd() / "acrolein_water.pot",
 72    outdir: Path = Path("dist"),
 73) -> dict[str, list[float]]:
 74    """
 75    Run FAST-VQE on acrolein in water using embedding.
 76
 77    First the Polarizable Embedding is used to model the solvent (water) environment,
 78    and then we apply the Projection based Wavefunction-in-DFT embedding scheme
 79    to treat a specific part of the acrolein molecule using our
 80    quantum computer algorithm.
 81    """
 82    solvent = qc.solvent_creator().frame(potential_file=Path(path_to_molecule_pot_file)).create()
 83
 84    # Create a Molecule configuration (molecule + basis, spin, etc)
 85    # This originate from a large pdb file, so the outer region is
 86    # now defined by the Solvent potential_file, and the
 87    # embedded_atoms define the "active region"
 88    # So here we chose atom index 0 and 1: Oxygen and Carbon
 89    molecular_configuration = qc.build_molecular_configuration(
 90        molecule=Path(path_to_molecule_xyz_file),
 91        basis_set="sto3g",
 92        charge=0,
 93        spin_difference=0,
 94        units="angstrom",
 95        solvent=solvent,
 96        embedded_atoms=[0, 1],
 97    )
 98
 99    pe_problem_builder = (
100        qc
101        .problem_builder_creator()
102        .ground_state()
103        .projective_embedding()
104        .add_problem_modifier()
105        # We further define a Complete Active Space (CAS) here we define CAS(10,10)
106        .active_space(
107            number_of_active_spatial_orbitals=10,
108            number_of_active_alpha_electrons=5,
109        )
110        .choose_orbital_assigner()
111        .total_weight(assignment_tolerance=0.5)
112        .create()
113    )
114
115    # We apply the Projection based Wavefunction-in-DFT embedding scheme on the molecule.
116    pe_problem = pe_problem_builder.build_restricted(molecular_configuration=molecular_configuration)
117
118    # Create an adaptive VQE based ground state calculator.
119    fast_vqe_calculator = (
120        qc
121        .calculator_creator()
122        .vqe()
123        .iterative()
124        .standard()
125        .with_options(
126            options=qc.options.IterativeVqeOptions(
127                max_iterations=max_iterations,
128            )
129        )
130        .choose_stopping_criterion()
131        .patience(patience=5, threshold=1e-4)
132        .create()
133    )
134
135    # Calculate the ground state energy using the FAST-VQE calculator using our simulator
136    energy_result = fast_vqe_calculator.calculate(pe_problem)
137
138    # You could create a standard CAS-CI ground state calculator.
139    if calculate_casci:
140        casci_calculator = (
141            qc
142            .calculator_creator()
143            .configuration_interaction()
144            .standard()  # FCI (full space) or CASCI (if problem has an active space)
145            .create()
146        )
147        casci_result = casci_calculator.calculate(pe_problem)
148        casci_energy = casci_result.total_energy.value
149    else:
150        # Here we cheat and provide the energy directly
151        casci_energy = -187.72474201929262  # CASCI-in-DFT value
152
153    energies = energy_result.total_energy_per_macro_iteration_with_initial_energy_and_final_energy.values
154    output: dict[str, list[float]] = {
155        "energies": energies,
156        "casci_energy": [casci_energy],
157    }
158    plot_energies(output, outdir=outdir)
159
160    return output
161
162
163if __name__ == "__main__":
164    output_dict = main()

Explanation

  1. Import Kvantify Qrunch

from pathlib import Path

from matplotlib import pyplot as plt
import numpy as np

import qrunch as qc

import qrunch as qc loads the public, stable API. For most users this is the only import needed, and the one guaranteed to remain stable across minor versions.

  1. The plotting method

def plot_energies(
    result: dict[str, list[float]],
    *,
    outdir: Path = Path("dist"),
    show: bool = False,
) -> None:
    """Plot the energies from FAST-VQE and CAS-CI."""
    energies = np.asarray(result["energies"])
    reference_energy = result["casci_energy"][0]

    iterations = np.asarray(list(range(len(energies))))
    # Create the plot
    fig1 = plt.figure()
    plt.plot(iterations, energies, color="darkgreen", marker="o", label="FAST-VQE")

    # Add a dashed horizontal line at the CAS-CI value.
    plt.axhline(y=reference_energy, color="black", linestyle="--", label="CAS-CI")

    plt.xlabel("Iteration number")
    plt.ylabel("Total energy [Hartree]")
    plt.title("FAST-VQE convergence")
    plt.legend()

    convergence_plot = outdir / "frame_vqe_convergence.png"
    fig1.savefig(convergence_plot, dpi=200, bbox_inches="tight")

    # -----------------------------
    # Figure 2: |Energy - Reference| (log scale)
    # -----------------------------
    # Use absolute error; add a tiny epsilon to avoid log(0) if we hit the reference exactly.
    eps = np.finfo(float).eps
    abs_error = np.abs(energies - reference_energy) + eps

    fig2 = plt.figure()
    plt.semilogy(iterations, abs_error, "-*", color="darkgreen", label="|E - E_ref|")
    plt.xlabel("Iteration number")
    plt.ylabel("Absolute error [Hartree]")
    plt.title("FAST-VQE error vs reference")
    plt.legend()
    plt.tight_layout()

    error_convergence_plot = outdir / "frame_vqe_error_semilogy.png"
    fig2.savefig(error_convergence_plot, dpi=200, bbox_inches="tight")

    if show:
        plt.show()
    else:
        # Close figures to free memory when running in batch contexts.
        plt.close(fig1)
        plt.close(fig2)

This code plot the convergence of the VQE energy and the error relative to the reference. Complete Active Space Configuration Interaction (CASCI) energy.

The code generates a PNG file under dist/:

  • frame_vqe_convergence.png

  • frame_vqe_error_semilogy.png

  1. Define the molecular configuration

    solvent = qc.solvent_creator().frame(potential_file=Path(path_to_molecule_pot_file)).create()

    # Create a Molecule configuration (molecule + basis, spin, etc)
    # This originate from a large pdb file, so the outer region is
    # now defined by the Solvent potential_file, and the
    # embedded_atoms define the "active region"
    # So here we chose atom index 0 and 1: Oxygen and Carbon
    molecular_configuration = qc.build_molecular_configuration(
        molecule=Path(path_to_molecule_xyz_file),
        basis_set="sto3g",
        charge=0,
        spin_difference=0,
        units="angstrom",
        solvent=solvent,
        embedded_atoms=[0, 1],
    )

We specify the acrolein + water system from the .xyz file, set the STO-3G basis, neutral charge, and zero spin difference. The key addition is the solvent (The SolventFrame) that include the polarizable embedding potential.

Furthermore, the embedded atoms are defined as the first two atoms (O and C of acrolein).

  1. Build the ground-state problem

    pe_problem_builder = (
        qc
        .problem_builder_creator()
        .ground_state()
        .projective_embedding()
        .add_problem_modifier()
        # We further define a Complete Active Space (CAS) here we define CAS(10,10)
        .active_space(
            number_of_active_spatial_orbitals=10,
            number_of_active_alpha_electrons=5,
        )
        .choose_orbital_assigner()
        .total_weight(assignment_tolerance=0.5)
        .create()
    )

    # We apply the Projection based Wavefunction-in-DFT embedding scheme on the molecule.
    pe_problem = pe_problem_builder.build_restricted(molecular_configuration=molecular_configuration)

We create a standard ground-state problem using the problem builder.

  1. Create and run the FAST-VQE calculator

    # Create an adaptive VQE based ground state calculator.
    fast_vqe_calculator = (
        qc
        .calculator_creator()
        .vqe()
        .iterative()
        .standard()
        .with_options(
            options=qc.options.IterativeVqeOptions(
                max_iterations=max_iterations,
            )
        )
        .choose_stopping_criterion()
        .patience(patience=5, threshold=1e-4)
        .create()
    )

    # Calculate the ground state energy using the FAST-VQE calculator using our simulator
    energy_result = fast_vqe_calculator.calculate(pe_problem)

We specify a maximum number of iterations, build and run the FAST-VQE calculator.

  1. Compute a high-accuracy FCI reference

    # You could create a standard CAS-CI ground state calculator.
    if calculate_casci:
        casci_calculator = (
            qc
            .calculator_creator()
            .configuration_interaction()
            .standard()  # FCI (full space) or CASCI (if problem has an active space)
            .create()
        )
        casci_result = casci_calculator.calculate(pe_problem)
        casci_energy = casci_result.total_energy.value
    else:
        # Here we cheat and provide the energy directly
        casci_energy = -187.72474201929262  # CASCI-in-DFT value

We build a CASCI calculator and calculate the CASCI energy for the same problem.

  1. Make the plots

    energies = energy_result.total_energy_per_macro_iteration_with_initial_energy_and_final_energy.values
    output: dict[str, list[float]] = {
        "energies": energies,
        "casci_energy": [casci_energy],
    }
    plot_energies(output, outdir=outdir)

We plot the VQE convergence and the error relative to the CASCI reference, by calling the plotting function defined above.

Running the Example

After saving the script as run_acrolein.py, you can run it directly from the command line:

$ python run_acrolein.py

You should see 2 plots:

FAST-VQE convergence plot FAST-VQE convergence plot

Naturally you can increase the convergence criteria, and max iterations to get a more accurate result.