Guided tour

The ‘cthyb-segment’ solver is a continuous-time Quantum Monte Carlo solver aimed at solving an Anderson impurity model with static and dynamical density-density interactions.

Below is a very short example demonstrating the basic usage of the solver.

from pytriqs.gf import *
from pytriqs.archive import *
from pytriqs.operators import *
import pytriqs.utility.mpi as mpi
from triqs_ctseg import Solver

#Green's function: two 1x1 blocks
gf_struct = [("up",[0]), ("down",[0])]
#initialize the solver
S = Solver(
     beta = 10.0, #inverse temperature
     gf_struct = gf_struct )

#set up the static interaction
U = 2.0 #static interaction

#impurity non-interacting Green's function
mu = 1.0 #chemical potential
eps0 = 0.0 #fermionic bath level
V = 0.5 #coupling to fermionic level
S.G0_iw << \
   inverse(iOmega_n + mu \
           -V**2 * inverse(iOmega_n - eps0)\
           -V**2 * inverse(iOmega_n + eps0))

#dynamical charge interactions
w0 = 1.0 #bosonic bath level
lambd = 0.4 #coupling to bosonic level
from itertools import product
for b1,b2 in product(dict(gf_struct).keys(), repeat=2):
 S.D0_iw[b1+"|"+b2] << \
   lambd**2*(inverse(iOmega_n - w0) \
                 - inverse(iOmega_n + w0))

#run the Monte-Carlo simulation
S.solve(
 h_int = U * n('up',0)*n('down',0), #local Hamiltonian
 n_cycles  = 10000, #Monte-Carlo cycles
 length_cycle = 100, #length of a cycle
 n_warmup_cycles = 1000,#thermalization cycles
 n_w_b_vertex = 2, #bosonic freqs. in vertex
 n_w_f_vertex = 20,#fermionic freqs. in vertex
 measure_nn = True, #density-density corr.
 measure_nnw = True, #density-density corr. fct
 measure_gw = True, #one-particle Green fct
 measure_fw = True, #one-particle improved est.
 measure_g2w = True,#3-leg corr.function
 measure_g3w = True,#4-leg corr. function
)

#store results in a HDF5 archive
if mpi.is_master_node():
 with HDFArchive("ct_seg.h5",'a') as A:
  A['G_tau'] = S.G_tau
  A['nn'] = S.nn
  A['nn_iw'] = S.nn_iw
  A['Sigma_iw']  = S.Sigma_iw
  A["G_2w"] = S.G_2w
  A["G_3w"] = S.G_3w

 ---Output:---
--------------------------------------------------------------------------
[[54840,1],0]: A high-performance Open MPI point-to-point messaging module
was unable to find any relevant network interfaces:

Module: OpenFabrics (openib)
  Host: 91f1d395389c

Another transport will be used instead, although this may result in
lower performance.

NOTE: You can disable this warning by setting the MCA parameter
btl_base_warn_component_unused to 0.
--------------------------------------------------------------------------

╔╦╗╦═╗╦╔═╗ ╔═╗  ┌─┐┌┬┐┌─┐┌─┐┌─┐
 ║ ╠╦╝║║═╬╗╚═╗  │   │ └─┐├┤ │ ┬
 ╩ ╩╚═╩╚═╝╚╚═╝  └─┘ ┴ └─┘└─┘└─┘
mu = [1.16,1.16]
U = 
[[0,1.68]
 [1.68,0]]
dynamical_U


Warming up ...
18:53:27  67% ETA 00:00:00 cycle 670 of 1000



Accumulating ...
18:53:27   0% ETA 00:00:18 cycle 52 of 10000
18:53:29  11% ETA 00:00:17 cycle 1100 of 10000
18:53:31  24% ETA 00:00:14 cycle 2433 of 10000
18:53:34  41% ETA 00:00:11 cycle 4109 of 10000
18:53:38  61% ETA 00:00:07 cycle 6190 of 10000
18:53:43  87% ETA 00:00:02 cycle 8761 of 10000


[Rank 0] Collect results: Waiting for all mpi-threads to finish accumulating...
rank 	 n,m 	 average +/- error 	 autocorrelation_time
0	 0,0 	 0.595634 +/- 0.0106237	0.0562727
0	 10,0 	 0.0962339 +/- 0.0215287	0.0210831
0	 -10,0 	 0.103848 +/- 0.0214855	0.0894778
[Rank 0] Timings for all measures:
Measure                     | seconds   
G(omega) measurement        | 4.71295   
G(tau) measurement          | 0.0203364 
density measurement         | 0.00136901
measurement of the MC sign  | 0.000896263
nn(omega) measurement       | 0.0238941 
three-frequency measurement | 3.49349   
two-frequency measurement   | 0.0284891 
Total measure time          | 8.28143   
[Rank 0] Acceptance rate for all moves:
Move  segment insertion: 0.524754
Move  segment removal: 0.524461
[Rank 0] Warmup lasted: 0.148543 seconds [00:00:00]
[Rank 0] Simulation lasted: 19.1076 seconds [00:00:19]
[Rank 0] Number of measures: 10000
Total number of measures: 10000
Starting run with 1 MPI rank(s) at : 2020-08-31 18:53:26.671155

To use the solver, one has to go through the following steps:

Import the Solver class:

from triqs_ctseg import Solver

Declare and construct a Solver instance:

gf_struct = [("up",[0]), ("down",[0])]
S = Solver(
     beta = 10.0, #inverse temperature
     gf_struct = gf_struct )

Here, gf_struct is a Python dictionary that specifies the block structure of the Green’s function, here one \(1\times 1\) ([0] is of length 1) for each spin (up, down).

Some construction parameters are mandatory, others are optional. See here for a complete list.

Initialize the inputs of the solver. In this script, we want to define the noninteracting Green’s function as:

\[\begin{equation} \left[\mathcal{G}_{0}\right]^{\sigma}(i\omega_{n})=\left[i\omega_{n}-\varepsilon_{a\sigma}-\Delta_{ab}^{\sigma}(i\omega_{n})\right]^{-1},\label{eq:G0_def} \end{equation}\]

with

\[\begin{split}\begin{align*} \Delta^\sigma(i\omega_{n}) & =\frac{V^{2}}{i\omega_{n}-\varepsilon_{0}} +\frac{V^{2}}{i\omega_{n}+\varepsilon_{0}},\\ \end{align*}\end{split}\]

In the script, this is done via an accessor called G0_iw. It is set by the user in the following way:

S.G0_iw << \
   inverse(iOmega_n + mu \
           -V**2 * inverse(iOmega_n - eps0)\
           -V**2 * inverse(iOmega_n + eps0))

Similarly, we want to define dynamical interactions:

\[\begin{align*} \mathcal{D}_{0}^{\sigma,\sigma'}(i\Omega_{m}) & =\frac{2\lambda^{2}\omega_0}{(i\Omega_{m})^{2}-\omega_{0}^2}. \end{align*}\]

This is done through the D0_iw accessor:

for b1,b2 in product(dict(gf_struct).keys(), repeat=2):
 S.D0_iw[b1+"|"+b2] << \
   lambd**2*(inverse(iOmega_n - w0) \
                 - inverse(iOmega_n + w0))

Other accessors are used to read out the observables after running the solver, like at the end of script:

A["G_tau"] = S.G_tau
A['nn'] = S.nn
A['nn_iw'] = S.nn_iw
A['Sigma_iw']  = S.Sigma_iw
A["G_2w"] = S.G_2w
A["G_3w"] = S.G_3w

Here, G_tau is the one-particle Green’s function \(G^\sigma(\tau)\), nn the density-density correlator \(\chi^{\sigma\sigma'}_{ab}\), nn_iw the density-density correlation function \(\chi^{\sigma\sigma'}_{ab}(i\Omega)\), Sigma_iw the self-energy \(\Sigma^\sigma(i\omega)\), G_2w the three-point correlation function \(\chi^{\sigma\sigma'}_{abc}(i\omega,i\Omega)\) and G_3w the four-point correlation function \(G^{2,\sigma\sigma'}_{abcd}(i\omega,i\omega',i\Omega)\). All these observables are defined here or there.

For a complete list of the accessors (both input and output), see here.

Call the solve() method:

S.solve(
 h_int = U * n('up',0)*n('down',0), #local Hamiltonian
 n_cycles  = 100000, #Monte-Carlo cycles
 length_cycle = 100, #length of a cycle
 n_warmup_cycles = 1000,#thermalization cycles
 n_w_b_vertex = 2, #bosonic freqs. in vertex
 n_w_f_vertex = 20,#fermionic freqs. in vertex
 measure_nn = True, #density-density corr.
 measure_nnw = True, #density-density corr. fct
 measure_gw = True, #one-particle Green fct
 measure_fw = True, #one-particle improved est.
 measure_g2w = True,#3-leg corr.function
 measure_g3w = True,#4-leg corr. function
)

h_int is the local many-body Hamiltonian, it is an operator expression. n_cycles, length_cycle and n_warmup_cycles respectively set the number of configuration measurements, the number of Monte-Carlo updates between two measurements and the number of thermalization cycles. n_w_b_vertex and n_w_f_vertex set the number of bosonic and fermionic Matsubara frequencies in the multi-variable correlation functions, and the last six options switch on the measurement of the corresponding observables.

This method comes with mandatory and optional parameters, whose complete list can be found here.

To know more about the solver, visit one of the following pages: