Generating new templates for custom models
Here we discuss how to add new models to the template library.
As an example, we consider the case of two-orbital impurity model
(model 2orb-UJ) that comes predefined in the library of default templates.
The files are located in the directory templates/2orb-UJ/QS.
The model is defined in 2orb-UJ.m:
def2ch[nrimp=2];
impuritybasis = {d[1],d[2]};
Heps1 = eps1 number[d[1]];
Hint1 = U1 hubbard[d[1]];
Hhyb1 = gammaPolCh[1] hop[f[0], d[1]]; (* f[0] is the 0-th site of the Wilson chain for channel 1 *)
Heps2 = eps2 number[d[2]];
Hint2 = U2 hubbard[d[2]];
Hhyb2 = gammaPolCh[2] hop[f[1], d[2]]; (* f[1] is the 0-th site of the Wilson chain for channel 2 *)
H12 = U12 nc[number[d[1]], number[d[2]]]+ J12 spinspin[d[1], d[2]];
Heps = Heps1 + Heps2;
Hint = Hint1 + Hint2 + H12;
Hhyb = Hhyb1 + Hhyb2;
H = H0 + Heps + Hint + Hhyb;
selfopd1 = ( Chop @ Expand @ komutator[Hint, d[#1, 1, #2]] )&;
selfopd2 = ( Chop @ Expand @ komutator[Hint, d[#1, 2, #2]] )&;
The first line declares that this is a two-channel problem with two impurity orbitals.
The second line defines how we label the impurity orbitals. Here we choose consecutively
numbered d[1] and d[2]. The main part of the file then defines various parts
of the Hamiltonian, making use of sneg function
calls such as number, Hubbard, spinspin, etc. Finally, in the last two lines
we define the auxiliary operators which are used in the “self-energy trick” to obtain
a high-quality estimate for the self-energy. They are obtained by taking the commutator
between the particle creation/annihilation operators and the interacting part of the
Hamiltonian.
The input file to NRG template generator is param:
[extra]
U1=
U2=
eps1=
eps2=
U12=
J12=
[param]
options=GENERATE_TEMPLATE
model=2orb-UJ.m
symtype=QS
ops=n_d1 n_d1^2 n_d2 n_d2^2 S_d1S_d2 n_d1n_d2 A_d1 A_d2 self_d1 self_d2 sigma_d1 sigma_d2
data_has_rescaled_energies=false
Here we declare all parameters (although the values need not be given), instruct
nrginit to generate a template file data.in rather than a specific instantiation
data, define the symmetry type (here QS), and list all operators of interest.
The last line is mandatory for generating NRGLjubljana_interface template; it instructs
the code to generate Hamiltonians which are not scaled in terms of the initial energy
scale of the problem, as is customarily done in NRG calculations.
The operators of interest are defined in modeloperators.m:
tt={};
tt = Join[tt, mtSingletOp["n_d1", number[d[1]] ]];
tt = Join[tt, mtSingletOp["n_d2", number[d[2]] ]];
tt = Join[tt, mtSingletOp["n_d1^2", pow[number[d[1]],2] ]];
tt = Join[tt, mtSingletOp["n_d2^2", pow[number[d[2]],2] ]];
tt = Join[tt, mtSingletOp["S_d1S_d2", spinspin[d[1], d[2]] ]];
tt = Join[tt, mtSingletOp["n_d1n_d2", nc[number[d[1]], number[d[2]]] ]];
tt = Join[tt, mtDoubletOp["A_d1", d[1] ]];
tt = Join[tt, mtDoubletOp["A_d2", d[2] ]];
tt = Join[tt, mtDoubletOp["self_d1", selfopd1 ]];
tt = Join[tt, mtDoubletOp["self_d2", selfopd2 ]];
tt = Join[tt, mtTripletOp["sigma_d1", d[1] ]];
tt = Join[tt, mtTripletOp["sigma_d2", d[2] ]];
tt
Here we again make use of sneg functions, as well as of high-level functions defined
in initial.m which is part of the nrginit Mathematica code.
The template files are then obtained by running nrginit. This will generate
data.in, ham_* and op.* template files. The log file with some details about
the generation process is saved as mmalog.
The main template data.in is lengthy and we show only a small chunk:
# Input file for NRG Ljubljana
# symtype QS
# Using sneg version 1.251
#!9
# Number of channels, chain sites, subspaces:
2 0 15
# SCALE 1
# Energies (GS energy subtracted, multiplied by 1/SCALE):
-4 1
1
DIAG ham_-4.1
-3 2
4
DIAG ham_-3.2
-2 1
10
DIAG ham_-2.1
-2 3
6
DIAG ham_-2.3
-1 2
20
DIAG ham_-1.2
-1 4
4
DIAG ham_-1.4
0 1
20
DIAG ham_0.1
0 3
15
DIAG ham_0.3
0 5
1
DIAG ham_0.5
1 2
20
DIAG ham_1.2
1 4
4
DIAG ham_1.4
2 1
10
DIAG ham_2.1
2 3
6
DIAG ham_2.3
3 2
4
DIAG ham_3.2
4 1
1
DIAG ham_4.1
# Irreducible matrix elements for Wilson chains:
f 0 0
20
4 1 3 2
0. 0. 1.4142135623730951 0.
3 2 2 3
0. 1.224744871391589 0. 0. 0. 0.
0. 0. 1.224744871391589 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. -1.224744871391589
3 2 2 1
0. 0. 0. -0.7071067811865475 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. -0.7071067811865475 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 1. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. -0.7071067811865475 0.
The ham_* contain the symbolic expressions for the Hamiltonian matrices in various invariant subspaces.
For instance, in the half-filled \(Q=0\) spin triplet \((2S+1)=3\) subspace,
we have ham_0.3:
{{eps2 + coefzeta[2, 0], 0, -gammaPolCh[1], 0, 0, gammaPolCh[2], 0, 0, 0, 0,
0, 0, 0, 0, 0}, {0, eps1 + coefzeta[2, 0], 0, 0, 0, 0,
gammaPolCh[2]/Sqrt[2], 0, -(gammaPolCh[2]/Sqrt[6]), 0, 0,
(-2*gammaPolCh[2])/Sqrt[3], 0, 0, 0}, {-gammaPolCh[1], 0,
eps1 + eps2 + J12/4 + U12 - coefzeta[1, 0] + coefzeta[2, 0], 0, 0, 0, 0, 0,
0, -gammaPolCh[2], 0, 0, 0, 0, 0}, {0, 0, 0, eps2 + coefzeta[1, 0], 0, 0,
-(gammaPolCh[1]/Sqrt[2]), 0, -(Sqrt[3/2]*gammaPolCh[1]), 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, eps1 + coefzeta[1, 0], 0, 0, -gammaPolCh[1], 0, 0, 0, 0,
gammaPolCh[2], 0, 0}, {gammaPolCh[2], 0, 0, 0, 0, 2*eps2 + U2, 0, 0, 0,
gammaPolCh[1], 0, 0, 0, 0, 0}, {0, gammaPolCh[2]/Sqrt[2], 0,
-(gammaPolCh[1]/Sqrt[2]), 0, 0, eps1 + eps2 - (3*J12)/4 + U12, 0, 0, 0,
-(gammaPolCh[1]/Sqrt[2]), 0, 0, gammaPolCh[2]/Sqrt[2], 0},
{0, 0, 0, 0, -gammaPolCh[1], 0, 0, 2*eps1 + U1, 0, 0, 0, 0, 0, 0,
-gammaPolCh[2]}, {0, -(gammaPolCh[2]/Sqrt[6]), 0,
-(Sqrt[3/2]*gammaPolCh[1]), 0, 0, 0, 0, eps1 + eps2 + J12/4 + U12, 0,
-(Sqrt[3/2]*gammaPolCh[1]), 0, 0, -(gammaPolCh[2]/Sqrt[6]), 0},
{0, 0, -gammaPolCh[2], 0, 0, gammaPolCh[1], 0, 0, 0,
eps1 + 2*eps2 + 2*U12 + U2 - coefzeta[1, 0], 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0, -(gammaPolCh[1]/Sqrt[2]), 0, -(Sqrt[3/2]*gammaPolCh[1]),
0, 2*eps1 + eps2 + U1 + 2*U12 - coefzeta[1, 0], 0, 0, 0, 0},
{0, (-2*gammaPolCh[2])/Sqrt[3], 0, 0, 0, 0, 0, 0, 0, 0, 0,
eps1 + eps2 + J12/4 + U12, 0, (-2*gammaPolCh[2])/Sqrt[3], 0},
{0, 0, 0, 0, gammaPolCh[2], 0, 0, 0, 0, 0, 0, 0,
eps1 + eps2 + J12/4 + U12 + coefzeta[1, 0] - coefzeta[2, 0], 0,
gammaPolCh[1]}, {0, 0, 0, 0, 0, 0, gammaPolCh[2]/Sqrt[2], 0,
-(gammaPolCh[2]/Sqrt[6]), 0, 0, (-2*gammaPolCh[2])/Sqrt[3], 0,
eps1 + 2*eps2 + 2*U12 + U2 - coefzeta[2, 0], 0},
{0, 0, 0, 0, 0, 0, 0, -gammaPolCh[2], 0, 0, 0, 0, gammaPolCh[1], 0,
2*eps1 + eps2 + U1 + 2*U12 - coefzeta[2, 0]}}
In a similar fashion, op.* files contain the symbolic forms of the operators. For instance
one part of the auxiliary operator self_d2 takes the following form:
{{0., 0., 0., 0.}, {0., 0., 0., 0.}, {0., 0., 0., 0.}, {0., 0., 0., 0.},
{0., 0., 0., 0.}, {1.1547005383792517*U2, 0., 0., 0.},
{0., -0.6123724356957946*J12 + 0.816496580927726*U12, 0., 0.},
{0., 0., 0., 0.}, {0., -0.11785113019775795*J12 - 0.4714045207910318*U12,
0., 0.}, {0., 0., 0.2886751345948129*J12 - 1.1547005383792517*U12 -
1.1547005383792517*U2, 0.}, {0., 0., 0., 0.},
{0., -0.08333333333333336*J12 - 0.3333333333333334*U12, 0., 0.},
{0., 0., 0., 0.}, {0., 0., 0., 0.2886751345948129*J12 -
1.1547005383792517*U12 - 1.1547005383792517*U2}, {0., 0., 0., 0.}}
We then need to define the structure of Green’s functions, separately for
Green’s functions and for dynamical susceptibility functions, in gf_stuct
and chi_struct, respectively. Here they are the same:
imp 2
imp 2
The quantities that need to be computed are listed in file info:
ops:n_d1 n_d1^2 n_d2 n_d2^2 S_d1S_d2 n_d1n_d2 A_d1 A_d2 self_d1 self_d2 sigma_d1 sigma_d2
specs:n_d1-n_d1 n_d1-n_d2 n_d2-n_d1 n_d2-n_d2
specd:A_d1-A_d1 A_d2-A_d2 self_d1-A_d1 self_d2-A_d2
spect:sigma_d1-sigma_d1 sigma_d1-sigma_d2 sigma_d2-sigma_d1 sigma_d2-sigma_d2
specq:
specot:
params:U1 U2 eps1 eps2 U12 J12
polarized:
Finally, NRGLjubljana_interface requires a number of initialization, wrapper and processing scripts.
The first, prepare, simply copies the required files to the temporary working directory:
#!/bin/bash
DIR="$( cd "$( dirname "${BASH_SOURCE[0]}" )" >/dev/null 2>&1 && pwd )"
cp $DIR/data.in $DIR/ham* $DIR/op.* $DIR/discretize $DIR/instantiate $DIR/wilson $DIR/process $DIR/*_struct $DIR/info .
The second performs the discretization:
#!/bin/sh
mv param param.orig
cp param.orig param
echo dos=Gamma_imp_00.dat >>param
adapt P
adapt N
mv FSOL.dat FSOL_imp_00.dat
mv FSOLNEG.dat FSOLNEG_imp_00.dat
cp param.orig param
echo dos=Gamma_imp_11.dat >>param
adapt P
adapt N
mv FSOL.dat FSOL_imp_11.dat
mv FSOLNEG.dat FSOLNEG_imp_11.dat
Then comes the most complex one which instantiates the input file for NRG calculation, numerical-valued data,
from the symbolic-expression template data.in:
#!/usr/bin/env perl
# Process a 'data.in' template file: fill in the numerical values,
# diagonalize the Hamiltonian matrices, perform the transformations.
# Rok Zitko, 2009-2019
use strict;
use warnings;
use File::Copy;
use Math::Trig;
# Paths to tools used
my $matrix = "matrix";
my $diag = "diag -B";
my $unitary = "unitary -B";
my $verbose = 1;
# Split param file into sections
my $fnparam = "param";
open (F, "<$fnparam") or die "Can't open $fnparam.";
my $section = "";
while (<F>) {
if (/^\[(.*)\]$/) {
close O if $section ne "";
$section = $1;
my $fnout = "$fnparam.$section";
open (O, ">$fnout") or die "Cant open $fnout.";
}
print O if $section ne "";
}
close O if $section ne "";
close F;
system "./wilson";
my $nr = getparam("Nmax");
### Translate data.in -> data
my $fnin = "data.in";
open(F, "<$fnin") or die "Can't open $fnin.";
my $fnout = "data";
open(OUT, ">$fnout") or die "Can't open $fnout.";
my $mode = 0; # Current parser mode
my $factor = -1; # Rescale factor for energies
my $smallest = 1e9; # Lowest rescaled energy (Egs=smallest/factor)
my %dim; # Dimensions of invariant subspaces
my $subspaces = -1; # Number of different subspaces
my $Egs = -999;
my $perch = 1;
if (getparam("polarized") eq "true") {
$perch = 2;
}
my $channels;
my $coefchannels;
my $fctr = 0;
while (<F>) {
print;
if (/^#/) {
if (/^# SCALE\s+(\S+)$/) {
# Extract the rescale factor
$factor = 1.0/$1;
}
# Copy comments verbatim
print OUT;
next;
}
if ($mode == 0) {
if (!/(\d+)\s+(\d+)\s+(\d+)/) { die "Parse error 1: $_"; }
$channels = $1;
$coefchannels = $channels * $perch;
my $impurities = $2;
$subspaces = $3;
$impurities = $nr; # override!
print "ch=$channels coefch=$coefchannels im=$impurities sub=$subspaces\n";
print OUT "$channels $impurities $subspaces\n";
$mode = 1;
next;
}
if ($mode == 1) {
my @buffer;
my $i;
for ($i = 1; $i <= $subspaces; $i++) {
if ($i != 1) { $_ = <F>; }
if (!/(\S+)\s+(\S+)/) { die "Parse error 2: $_"; }
my $qnspaces = "$1 $2";
my $qn = "$1.$2";
my $size;
chomp($size = <F>);
$dim{$qn} = $size; # store size!
print "[$i] $qn dim=$size\n";
my $line1;
chomp($line1 = <F>);
$line1 =~ /DIAG\s+(\S+)$/;
my $fn = $1;
system("$matrix -c$coefchannels param.extra $fn >ham");
if ($? != 0) { die "matrix call error"; }
copy("ham", "ham.$qn");
($factor > 0) or die "no factor?";
system("$diag -q -s $factor -o val -O vec ham");
my $eig = `cat val`;
push(@buffer, [ "$qnspaces\n$size\n", $eig ]);
if ($eig =~ /^([-0-9.]+)(\s|\n)/) {
my $min = $1;
if ($smallest > $min) { $smallest = $min; }
}
copy("vec", "vec.$qn");
}
$Egs=$smallest/$factor;
print "E_gs=$Egs\n";
@buffer = map { $_->[0].subtract($_->[1]) } @buffer;
print OUT @buffer;
$mode = 2;
next;
}
if ($mode == 2) {
if (!/f\s+(\d+)\s+(\d+)/) { die "Parse error 3: $_"; }
print OUT;
unitary();
$fctr++;
if ($fctr == $channels) {
$mode = 3;
}
next;
}
if ($mode == 3) {
if (!/^e/) { die "Parse error 5: $_"; }
my $egs = <F>;
print OUT "e\n$Egs\n";
$mode = 4;
next;
}
if ($mode == 4) {
if (/^[zZT]/) { # Discretization tables
$mode = 5;
next;
}
if (/^[spdvt]/) {
print OUT;
unitary();
next;
}
die "Parse error 6: $_";
}
}
my $adddisc = 1;
if ($adddisc) {
print OUT "z\n";
for (my $ch = 1; $ch <= $perch*$channels; $ch++) {
print OUT "$nr\n";
my $xi = `cat xi${ch}.dat`;
print OUT "$xi";
}
for (my $ch = 1; $ch <= $perch*$channels; $ch++) {
print OUT "$nr\n";
my $zeta = `cat zeta${ch}.dat`;
print OUT "$zeta";
}
}
sub subtract
{
my $vals = shift;
my @l = split(' ', $vals);
@l = map { $_-$smallest } @l;
return "@l\n";
}
sub unitary
{
chomp(my $nr = <F>);
print "nr=$nr\n";
print OUT "$nr\n";
for (my $i = 1; $i <= $nr; $i++) {
print "unitary i=$i\n";
chomp(my $line = <F>);
if (!($line =~ /(\S+)\s+(\S+)\s+(\S+)\s+(\S+)/)) { die "Parse error 4: $_"; }
my $qn1spaces = "$1 $2";
my $qn1 = "$1.$2";
my $dim1 = $dim{$qn1};
my $qn2spaces = "$3 $4";
my $qn2 = "$3.$4";
my $dim2 = $dim{$qn2};
print "$qn1 $dim1 $qn2 $dim2\n";
$line = <F>;
my $c = substr($line, 0, 1);
if ($c ne "o") {
open(O, ">mat") or die "Can't open mat for writing.";
for (my $n = 1; $n <= $dim1; $n++) {
print O $line;
if ($n != $dim1) { $line = <F>; }
} # for $n
close(O);
} else {
chomp(my $fn = $line);
-e $fn or die "File not found: $fn";
system("$matrix -c$coefchannels param.extra $fn >mat");
if ($? != 0) { die "matrix call error for $fn"; }
}
my $cmd = "$unitary -l -q -o mat.res vec.$qn1 mat vec.$qn2";
system($cmd);
if ($? != 0) { die "unitary call error: $cmd"; }
my $res = `cat mat.res`;
print OUT "$qn1spaces $qn2spaces\n$res";
} # for $i
}
my $target_dir = shift;
if (defined($target_dir) && -d $target_dir) {
system "cp param data $target_dir";
}
sub getparam
{
my $keyword = shift;
my $fn = "param";
open (G, "<", $fn) or die "Can't open $fn for reading: $!\n";
while (<G>) {
if (/^$keyword=(.*)/) {
return $1;
}
}
close G;
die "Parsing failure: $keyword not found.";
}
This script, in turn, calls another one that calculates the Wilson chain coefficients:
#!/usr/bin/env perl
system "mv param param.orig";
prochyb("imp_00",1);
prochyb("imp_11",2);
sub prochyb
{
my $x = shift;
my $n = shift;
my $hybfn = "Gamma_${x}.dat";
-e $hybfn or die "Hybridisation function $hybfn not found.";
-e "FSOL_${x}.dat" or die "FSOL_${x}.dat not found.";
-e "FSOLNEG_${x}.dat" or die "FSOLNEG_{x}.dat not found.";
system "cp param.orig param";
system "echo dos=Gamma_${x}.dat >>param";
system "cp FSOL_${x}.dat FSOL.dat";
system "cp FSOLNEG_${x}.dat FSOLNEG.dat";
system "nrgchain";
-e "xi.dat" && -e "zeta.dat" && -e "theta.dat" or die "Wilson chain coefficients not found.";
system "mv xi.dat xi${n}.dat";
system "mv zeta.dat zeta${n}.dat";
system "mv theta.dat theta${n}.dat";
}
The last script, process, calculates averages over the twist parameter and performs the broadening
of spectral functions:
#!/usr/bin/env perl
use strict;
use warnings;
use Math::Trig;
my $pi = pi();
my $Nz = getparam("Nz");
my $min = getparam("broaden_min");
my $max = getparam("broaden_max");
my $ratio = getparam("broaden_ratio");
my $T = getparam("T");
my $alpha = getparam("alpha");
my $gamma = getparam("gamma");
print "Nz=$Nz min=$min max=$max ratio=$ratio T=$T alpha=$alpha gamma=$gamma\n";
my $pr = "FDM_dens"; # prefix
avg("spec_${pr}_A_d1-A_d1.bin", "A_imp_00");
avg("spec_${pr}_A_d2-A_d2.bin", "A_imp_11");
avg("spec_${pr}_self_d1-A_d1.bin", "B_imp_00");
avg("spec_${pr}_self_d2-A_d2.bin", "B_imp_11");
scaley("A_imp_00", -$pi, "imG_imp_00");
scaley("A_imp_11", -$pi, "imG_imp_11");
scaley("B_imp_00", -$pi, "imF_imp_00");
scaley("B_imp_11", -$pi, "imF_imp_11");
kk("imG_imp_00");
kk("imG_imp_11");
kk("imF_imp_00");
kk("imF_imp_11");
avg("corr_${pr}_n_d1-n_d1.bin", "imNN_imp_00");
avg("corr_${pr}_n_d2-n_d2.bin", "imNN_imp_11");
avg("corr_${pr}_n_d1-n_d2.bin", "imNN_imp_01");
avg("corr_${pr}_n_d2-n_d1.bin", "imNN_imp_10");
avg("spin_${pr}_sigma_d1-sigma_d1.bin", "imSS_imp_00");
avg("spin_${pr}_sigma_d2-sigma_d2.bin", "imSS_imp_11");
avg("spin_${pr}_sigma_d1-sigma_d2.bin", "imSS_imp_01");
avg("spin_${pr}_sigma_d2-sigma_d1.bin", "imSS_imp_10");
kk("imNN_imp_00");
kk("imNN_imp_11");
kk("imNN_imp_01");
kk("imNN_imp_10");
kk("imSS_imp_00");
kk("imSS_imp_11");
kk("imSS_imp_01");
kk("imSS_imp_10");
sub kk {
my $in = shift;
$in .= ".dat";
my $out;
($out = $in) =~ s/^im/re/;
-e $in or die "Input file not found: $in. Stopped";
system("kk $in $out");
-e $out or die "Output file not found:$ out. Stopped";
}
sub avg
{
my $in = shift;
my $out = shift;
$out .= ".dat";
unlink "$out";
system "broaden -x $gamma -m $min -M $max -r $ratio $in $Nz $alpha $T 1e-9";
if ($?) {
system "echo Broadening of $in failed: $? >>>ERROR";
die "broaden failed: $?\n";
}
my $specfn = "spec.dat";
if (!-e $specfn) {
system "echo Failed to produce broadened spectrum $specfn -> $out. >>>ERROR";
die "broaden failed\n";
}
system "mv $specfn $out";
}
sub getparam
{
my $keyword = shift;
my $fn = "param";
open (F, "<", $fn) or die "Can't open $fn for reading: $!\n";
while (<F>) {
if (/^$keyword=(.*)/) {
return $1;
}
}
close F;
die "Parsing failure: $keyword not found.";
}
sub scaley
{
my $fnin = shift;
my $factor = eval(shift);
my $fnout = shift;
$fnin .= ".dat";
$fnout .= ".dat";
open (F, "<", $fnin) or die "Can't open $fnin for reading: $!\n";
open (G, ">", $fnout) or die "Can't open $fnout for writing: $!\n";
while (<F>) {
if (!/^#/) {
chomp;
my @d = split;
$d[1] = $d[1] * $factor;
print G "@d\n";
} else {
print G;
}
}
close G;
close F;
}