Examples

The following examples present some common use cases that show how to run qupled and how to post-process the results.

Setup a scheme and analyze the output

This example solves the STLS and shows two ways in which to access the results: Directly from the object used to perform the calculation or from the database used to store the results.

from qupled.postprocess.output import DataBase, OutputType
from qupled.schemes import stls

# Define teh scheme
scheme = stls.Solver()

# Define the input parameters
inputs = stls.Input(10.0, 1.0, mixing=0.5)

# Solve the scheme
scheme.compute(inputs)

# Access the internal energy from the database
results = DataBase.read_results(scheme.run_id, names=["uint"])
print("Internal energy from the output file: ")
print(results["uint"])

# Access the internal energy from the result class
print("Internal energy from the result class: ")
print(scheme.results.uint)

Define an initial guess

Qupled allows to specify a custom initial guess for any scheme. This example shows how to define an initial guess for the STLS scheme, but the same approach can be used for any other scheme.

from qupled.schemes import stls

# Define the object used to solve the scheme
scheme = stls.Solver()

# Define the input parameters
inputs = stls.Input(10.0, 1.0, mixing=0.2)

# Solve scheme
scheme.compute(inputs)

# Create a custom initial guess from the output files of the previous run
inputs.guess = scheme.get_initial_guess(scheme.run_id)

# Solve the scheme again with the new initial guess
scheme.compute(inputs)

Read the results from the database

Qupled stores the results of the calculations in a database that can be easily accessed to retrieve any quantity of interest. This example shows how to read the results from the database in order to plot the static structure factor for two different schemes.

import matplotlib.pyplot as plt
from qupled.postprocess.output import DataBase
from qupled.schemes import esa, rpa

# Solve two schemes: RPA and ESA
scheme = rpa.Solver()
scheme.compute(rpa.Input(10.0, 1.0))
run_id_rpa = scheme.run_id

scheme = esa.Solver()
scheme.compute(esa.Input(10.0, 1.0))
run_id_esa = scheme.run_id

# Retrieve information from the database
rpa_run = DataBase.read_run(run_id_rpa)
esa_run = DataBase.read_run(run_id_esa)


# Plot the static structure factor for the two schemes
plt.plot(
    rpa_run.results["wvg"],
    rpa_run.results["ssf"],
    color="b",
    label=rpa_run.inputs["theory"],
)
plt.plot(
    esa_run.results["wvg"],
    esa_run.results["ssf"],
    color="r",
    label=esa_run.inputs["theory"],
)
plt.legend(loc="lower right")
plt.xlabel("Wave vector")
plt.ylabel("Static structure factor")
plt.title(
    "State point : (coupling = "
    + str(rpa_run.inputs["coupling"])
    + ", degeneracy = "
    + str(rpa_run.inputs["degeneracy"])
    + ")"
)
plt.show()

Post-process the results

Qupled also includes a set of post-processing tools that can be used to further analyze the results of the calculations and extract valuable information about the system.

Compute additional quantities from an existing object

This example shows how the post-processing tools can be used to compute the radial distribution function and the imaginary-time correlation function, two quantities that are not directly computed in the iterative calculations but that can be easily obtained from the results of the calculation.

import matplotlib.pyplot as plt
from matplotlib import colormaps as cm
from qupled.schemes import stls

# Solve the STLS scheme in the ground state for rs=10.0
solver = stls.Solver()
solver.compute(stls.Input(10.0, 0.0, mixing=0.7))

# Check the default values for the quantites to be computed in the post-processing step
print(f"----- Default values for the post-processing quantities -----")
print(f"Radial distribution function: {solver.results.rdf}")
print(f"Imaginary-time correlation function: {solver.results.itcf}")
print("-------------------------------------------------------------")

# Compute the radial distribution function and the imaginary-time correlation function
solver.compute_rdf()
solver.compute_itcf(tau=[0.0, 0.1, 0.2, 0.5, 1.0, 2.0, 5.0])

# Plot the radial distribution function
plt.plot(solver.results.rdf_grid, solver.results.rdf, color="b")
plt.xlabel("Interparticle distance")
plt.ylabel("Radial distribution function")
plt.show()

# Plot the imaginary-time correlation function
colormap = cm["viridis"].reversed()
# nl = min(solver.results.itcf.shape[1], 10)
nl = solver.results.itcf.shape[1]
for i in range(nl):
    color = colormap(1.0 - 1.0 * i / nl)
    plt.plot(
        solver.results.wvg,
        solver.results.itcf[:, i],
        color=color,
        label=f"$\\tau^*={solver.results.tau[i]:.2f}$",
    )
plt.legend(loc="center right")
plt.xlabel("Wave-vector")
plt.ylabel("Imaginary-time correlation function")
plt.show()

Compute additional quantities from the database

The post-processing tools can also be used to compute additional quantities by reading the necessary data from the database. This is particularly useful when revisiting previous calculations that were stored without computing the additional quantities. This example shows how to compute the radial distribution function starting from run results stored in the database.

import matplotlib.pyplot as plt
from qupled.util.dimension import Dimension
from qupled.postprocess.correlation_functions import compute_rdf
from qupled.postprocess.output import DataBase
from qupled.schemes import stls

# Solve an STLS scheme
scheme = stls.Solver()
scheme.compute(stls.Input(2.0, 1.0, dimension=Dimension._2D, mixing=0.7))
run_id = scheme.run_id

# After solving, the RDF is not yet in the database
run = DataBase.read_run(run_id, result_names=["rdf", "rdf_grid"])
print(
    "Radial distribution function available before computation: ", "rdf" in run.results
)

# Compute the RDF and store it back under the same run_id
compute_rdf(run_id)

# The RDF is now available in the database
run = DataBase.read_run(run_id)
print(
    "Radial distribution function available after the computation: ",
    "rdf" in run.results,
)
plt.plot(run.results["rdf_grid"], run.results["rdf"], color="b")
plt.xlabel("Interparticle distance")
plt.ylabel("Radial distribution function")
plt.show()

Finite Size Correction

Qupled also includes a module that can be used to compute the finite size correction to the free energy. This example shows how to use this module to compute the finite size correction for a 3D system at \(r_s=5\) and \(\Theta=1\).

from qupled.postprocess import output, finite_size_correction as fsc
from qupled.schemes import rpa

# Define the scheme
solver = rpa.Solver()
inputs = rpa.Input(coupling=1.0, degeneracy=1.0, cutoff=50.0)

# Compute the finite size correction for different number of particles
run_ids = []
finite_size_correction = fsc.FiniteSizeCorrection()
for number_of_particles in [10, 100, 1000]:
    fsc_inputs = fsc.Input(
        number_of_particles=number_of_particles,
        drs=0.1,
        scheme=inputs,
    )
    finite_size_correction.compute(solver, fsc_inputs)
    run_ids.append(finite_size_correction.run_id)

# Fetch and print the results
col1, col2, col3 = 21, 27, 26
header = f"│ {'Number of Particles':>{col1}} │ {'uint correction':>{col2}} │ {'fxc correction':>{col3}} │"
top = f"┌{'─' * (col1 + 2)}┬{'─' * (col2 + 2)}┬{'─' * (col3 + 2)}┐"
mid = f"├{'─' * (col1 + 2)}┼{'─' * (col2 + 2)}┼{'─' * (col3 + 2)}┤"
bot = f"└{'─' * (col1 + 2)}┴{'─' * (col2 + 2)}┴{'─' * (col3 + 2)}┘"
print("\nFinite Size Correction Results:")
print(top)
print(header)
print(mid)
for run_id in run_ids:
    run_data = output.DataBase.read_run(
        run_id, type=output.OutputType.FINITE_SIZE_CORRECTION
    )
    number_of_particles = run_data.inputs["number_of_particles"]
    uint = run_data.results["uint"]
    fxc = run_data.results["fxc"]
    print(f"│ {number_of_particles:>{col1}} │ {uint:>{col2}.6e} │ {fxc:>{col3}.6e} │")
print(bot)

Computationally intestive schemes

The quantum and the VS schemes are computationally more demanding than the classical schemes. The following examples show how what can be done to speed up the calculations for these schemes.

VS Schemes: Pre-compute the free energy integrand

This example demonstrates how to solve the classical VS-STLS scheme at finite temperature. The calculation is first carried out up to \(r_s=2\) and then resumed up to \(r_s=5\) reusing the pre-computed free energy integrand. VS-type schemes can be numerically demanding, so it is often convenient to be able to restart from a known state point using previously computed quantities.

from qupled.schemes import vsstls

# Define the object used to solve the scheme
scheme = vsstls.Solver()

# Define the input parameters
inputs = vsstls.Input(2.0, 1.0, mixing=0.5, alpha=[-0.2, 0.2])

# Compute
scheme.compute(inputs)

# Load the free energy integrand computed for rs = 2.0
fxci = scheme.get_free_energy_integrand(scheme.run_id)

# Setup a new VSStls simulation for rs = 5.0
inputs.coupling = 5.0
inputs.alpha = [0.5, 0.7]
inputs.free_energy_integrand = fxci

# Compute
scheme.compute(inputs)

Quantum schemes: parallelization and pre-computation

There are two strategies that can be employed to speed up the calculations of quantum schemes:

Parallelization: qupled supports both multithreaded calculations with OpenMP and multiprocessors computations with MPI. OpenMP and MPI can be used concurrently by setting both the number of threads and the number of processes in the input dataclasses. The following example shows how to solve the quantum schemes using OpenMP parallelization.

from qupled.schemes import qstls

# Define a Qstls object to solve the QSTLS scheme
scheme = qstls.Solver()

# Define the input parameters
inputs = qstls.Input(10.0, 1.0, mixing=0.5, matsubara=16, threads=16)

# Solve the QSTLS scheme
scheme.compute(inputs)

Pre-computation: The calculations for the quantum schemes can be made significantly faster if part of the calculation of the auxiliary density response can be skipped. Qupled will look into the database used to store the results to try to find the necessary data to skip the full calculation of the auxiliary density response. Try to run the following example and notice how the second calculation is much faster than the first one.

from qupled.schemes import qstls

# Define the object used to solve the scheme
scheme = qstls.Solver()

# Define the input parameters
inputs = qstls.Input(10.0, 1.0, mixing=0.5, matsubara=16, threads=16)

# Solve the QSTLS scheme and store the internal energy (v1 calculation)
scheme.compute(inputs)
uInt1 = scheme.results.uint

# Repeat the calculation and recompute the internal energy (v2 calculation)
scheme.compute(inputs)
uInt2 = scheme.results.uint

# Compare the internal energies obtained with the two methods
print("Internal energy (v1) = %.8f" % uInt1)
print("Internal energy (v2) = %.8f" % uInt2)

# Change the coupling parameter
inputs.coupling = 20.0

# Compute with the updated coupling parameter
scheme.compute(inputs)

# Change the degeneracy parameter
inputs.degeneracy = 2.0

# Compute with the update degeneracy parameter (this will recompute the fixed component)
scheme.compute(inputs)