Example with dynamical charge and spin interactions

Below is a short example demonstrating how to set up a computation in the presence of dynamical interactions in the charge and longitudinal and transverse spin channels.

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

gf_struct = [("up",[0]), ("down",[0])]
S = Solver(beta = 20.0,  gf_struct = gf_struct)
#set noninteracting Green's function
mu, eps0, V = 2.0, 0.3, 1.0
S.G0_iw << inverse(iOmega_n + mu - V**2 * inverse(iOmega_n - eps0))

#set dynamical charge and longitudinal spin interactions
w_ch, w_sp = 1.0, 0.8 #screening frequencies
g_ch, g_sp = 0.4, 0.7 #coupling strengths
sigz = lambda s1, s2 : 1 if s1==s2 else -1
import itertools
from pytriqs.gf.descriptors import Function
for b1, b2 in itertools.product(dict(gf_struct).keys(), repeat=2):
 S.D0_iw[b1+"|"+b2] << Function(lambda w : g_ch**2*2*w_ch/(w**2-w_ch**2)\
                                         + g_sp**2*2*w_sp/(w**2-w_sp**2)*sigz(b1,b2) )
#set transverse spin interaction
S.Jperp_iw << Function(lambda w : 4 * g_sp**2 * 2*w_sp/(w**2-w_sp**2))

from pytriqs.operators import *
U = 2.0
S.solve(h_int = U * n('up',0)*n('down',0),
        n_cycles = 10000, length_cycle = 100, n_warmup_cycles = 100)

from pytriqs.archive import *
if mpi.is_master_node():
 with HDFArchive("dyn_interactions.h5",'w') as A:
  A['G_tau'] = S.G_tau

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

Module: OpenFabrics (openib)
  Host: 9a2e883de69c

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 = [2.7725,2.7725]
U = 
[[0,2.905]
 [2.905,0]]
dynamical_U
jperp_interactions with full spin rotational invariance

Warming up ...



Accumulating ...
18:50:58   4% ETA 00:00:02 cycle 468 of 10000
18:51:00  95% ETA 00:00:00 cycle 9592 of 10000


[Rank 0] Collect results: Waiting for all mpi-threads to finish accumulating...
[Rank 0] Timings for all measures:
Measure                    | seconds   
G(tau) measurement         | 0.0496031 
measurement of the MC sign | 0.00073307
Total measure time         | 0.0503361 
[Rank 0] Acceptance rate for all moves:
Move  segment insertion: 0.4986
Move  segment removal: 0.500748
Move  insert composite segment: 0.453249
Move  remove composite segment: 0.451513
Move  swap bosonic lines: 0.105738
[Rank 0] Warmup lasted: 0.0212322 seconds [00:00:00]
[Rank 0] Simulation lasted: 2.21189 seconds [00:00:02]
[Rank 0] Number of measures: 10000
Total number of measures: 10000
Starting run with 1 MPI rank(s) at : 2020-08-31 18:50:58.497778

Let us decompose the various 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 = 20.0,  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}\Delta^\sigma(i\omega_{n}) =\frac{V^{2}}{i\omega_{n}-\varepsilon_{0}}\\\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))

Similarly, we want to define dynamical interactions:

\[\begin{split}\begin{align} \mathcal{D}_{0}^{\sigma\sigma'}(i\Omega) & =\frac{2g_{\mathrm{ch}}^{2}\omega_{\mathrm{ch}}}{(i\Omega)^{2}-\omega_{\mathrm{ch}}^{2}}+\left(-\right)^{\sigma\sigma'}\frac{2g_{\mathrm{sp}}^{2}\omega_{\mathrm{sp}}}{(i\Omega)^{2}-\omega_{\mathrm{sp}}^{2}},\label{eq:D0_example}\\ \mathcal{J}_{\perp}(i\Omega) & =4\cdot\frac{2g_{\mathrm{sp}}^{2}\omega_{\mathrm{sp}}}{(i\Omega)^{2}-\omega_{\mathrm{sp}}^{2}}.\label{eq:Jperp} \end{align}\end{split}\]

This is done through the D0_iw and Jperp_iw accessors:

for b1, b2 in itertools.product(dict(gf_struct).keys(), repeat=2):
 S.D0_iw[b1+"|"+b2] << Function(lambda w : g_ch**2*2*w_ch/(w**2-w_ch**2)\
                                         + g_sp**2*2*w_sp/(w**2-w_sp**2)*sigz(b1,b2) )

S.Jperp_iw << Function(lambda w : 4 * g_sp**2 * 2*w_sp/(w**2-w_sp**2))

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),
        n_cycles = 10000, length_cycle = 100, n_warmup_cycles = 100)

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.

The fact that Jperp_iw is nonzero is automatically detected by the code and the corresponding Monte-Carlo updates are turned on.

The solve() 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: