Polarizable Embedding “The Frame”

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

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

Explanation

  1. Import the QDK library

# pyright: reportUnknownVariableType=false, reportMissingImports=false
# ruff: noqa
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 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.