Skip to content

Commit

Permalink
Change some names and types
Browse files Browse the repository at this point in the history
  • Loading branch information
j-otsuki committed May 8, 2018
1 parent 19eab2a commit 70721e4
Show file tree
Hide file tree
Showing 10 changed files with 157 additions and 137 deletions.
33 changes: 18 additions & 15 deletions c++/g2_parameters.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@
using triqs::hilbert_space::gf_struct_t;
using indices_t = triqs::hilbert_space::fundamental_operator_set::indices_t;
using four_indices_t = std::tuple<indices_t, indices_t, indices_t, indices_t>;
using three_freqs_t = std::vector< std::tuple<int, int, int> >;
using three_freqs_t = std::tuple<int, int, int>;

namespace pomerol2triqs {

Expand All @@ -21,11 +21,13 @@ namespace pomerol2triqs {
/// Inverse temperature
double beta;

/// Channel in which Matsubara frequency representation is defined.
/// Channel in which Matsubara frequency representation is defined
channel_t channel = PH;

/// Number of bosonic and fermionic Matsubara frequencies.
int n_b, n_f;
/// Number of non-negative bosonic Matsubara frequencies
int n_b;
/// Number of positive fermionic Matsubara frequencies
int n_f;

/// indices of operators in TRIQS convention: (block_name, inner_index)
indices_t index1, index2, index3, index4;
Expand All @@ -42,33 +44,34 @@ namespace pomerol2triqs {
/// Inverse temperature
double beta;

/// Channel in which Matsubara frequency representation is defined.
/// Channel in which Matsubara frequency representation is defined
channel_t channel = PH;

/// Number of bosonic and fermionic Matsubara frequencies.
int n_b, n_f;
/// Number of non-negative bosonic Matsubara frequencies
int n_b;
/// Number of positive fermionic Matsubara frequencies
int n_f;

/// set of indices of four operators in TRIQS convention: (block_name, inner_index)*4
std::vector<four_indices_t> vec_four_indices;
/// four operator indices in TRIQS convention: (block_name, inner_index)*4
std::vector<four_indices_t> four_indices;
};

struct g2_iw_freq_vec_params_t {
struct g2_iw_freq_fix_params_t {

/// Block structure of GF
gf_struct_t gf_struct;

/// Inverse temperature
double beta;

/// Channel in which Matsubara frequency representation is defined.
/// Channel in which Matsubara frequency representation is defined
channel_t channel = PH;

/// three frequencies (wb, wf1, wf2).
three_freqs_t three_freqs;
std::vector<three_freqs_t> three_freqs;

/// set of indices of four operators in TRIQS convention: (block_name, inner_index)*4
// indices_t index1, index2, index3, index4;
std::vector<four_indices_t> vec_four_indices;
/// four operator indices in TRIQS convention: (block_name, inner_index)*4
std::vector<four_indices_t> four_indices;
};

/*
Expand Down
44 changes: 21 additions & 23 deletions c++/pomerol_ed.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -469,7 +469,7 @@ namespace pomerol2triqs {
*/

// core function for computing G2
std::vector<std::complex<double> > pomerol_ed::compute_g2(gf_struct_t const &gf_struct, double beta, channel_t channel, indices_t index1, indices_t index2, indices_t index3, indices_t index4, three_freqs_t const &three_freqs){
triqs::arrays::array<std::complex<double>, 1> pomerol_ed::compute_g2_core(gf_struct_t const &gf_struct, double beta, channel_t channel, indices_t index1, indices_t index2, indices_t index3, indices_t index4, std::vector<three_freqs_t> const &three_freqs){

if (!matrix_h) TRIQS_RUNTIME_ERROR << "G2_iw: no Hamiltonian has been diagonalized";
compute_rho(beta);
Expand All @@ -492,34 +492,35 @@ namespace pomerol2triqs {
pom_g2.compute(false, {}, comm);

// compute g2 value using MPI
g2_iw_freq_vec_t g2( three_freqs.size() );
g2_iw_freq_fix_t g2( three_freqs.size() );
for(int i=comm.rank(); i<three_freqs.size(); i+=comm.size()){
int wb = std::get<0>(three_freqs[i]);
int wf1 = std::get<1>(three_freqs[i]);
int wf2 = std::get<2>(three_freqs[i]);
g2[i] = -pom_g2(wb+wf1, wf2, wf1);
g2(i) = -pom_g2(wb+wf1, wf2, wf1);
}

// broadcast results
for(int i=0; i<three_freqs.size(); i++){
int sender = i % comm.size();
boost::mpi::broadcast(comm, g2[i], sender);
boost::mpi::broadcast(comm, g2(i), sender);
}

// std::cout << "End freq loop: rank" << comm.rank() << std::endl;
return g2;
}

std::vector< std::vector<std::complex<double> > > pomerol_ed::compute_g2_indices_loop(gf_struct_t const &gf_struct, double beta, channel_t channel, std::vector<four_indices_t> const &vec_four_indices, three_freqs_t const &three_freqs){

std::vector<g2_iw_freq_vec_t> vec_g2;
std::vector< triqs::arrays::array<std::complex<double>, 1> > pomerol_ed::compute_g2(gf_struct_t const &gf_struct, double beta, channel_t channel, std::vector<four_indices_t> const &four_indices, std::vector<three_freqs_t> const &three_freqs){

std::vector<g2_iw_freq_fix_t> vec_g2;
// TODO: MPI
for( auto four_indices : vec_four_indices ){
indices_t index1 = std::get<0>(four_indices);
indices_t index2 = std::get<1>(four_indices);
indices_t index3 = std::get<2>(four_indices);
indices_t index4 = std::get<3>(four_indices);
vec_g2.push_back( compute_g2(gf_struct, beta, channel, index1, index2, index3, index4, three_freqs) );
for( auto ind4 : four_indices ){
indices_t index1 = std::get<0>(ind4);
indices_t index2 = std::get<1>(ind4);
indices_t index3 = std::get<2>(ind4);
indices_t index4 = std::get<3>(ind4);
vec_g2.push_back( compute_g2_core(gf_struct, beta, channel, index1, index2, index3, index4, three_freqs) );
}
return vec_g2;
}
Expand All @@ -532,8 +533,8 @@ namespace pomerol2triqs {
p2.channel = p.channel;
p2.n_b = p.n_b;
p2.n_f = p.n_f;
four_indices_t four_indices = std::make_tuple(p.index1, p.index2, p.index3, p.index4);
p2.vec_four_indices.push_back( four_indices );
four_indices_t ind4 = std::make_tuple(p.index1, p.index2, p.index3, p.index4);
p2.four_indices.push_back( ind4 );

std::vector<g2_iw_freq_box_t> vec_g2 = G2_iw_freq_box(p2);
return vec_g2[0];
Expand All @@ -542,7 +543,7 @@ namespace pomerol2triqs {
auto pomerol_ed::G2_iw_freq_box(g2_iw_freq_box_params_t const &p) -> std::vector<g2_iw_freq_box_t> {

// create a list of three frequencies, (wb, wf1, wf2)
three_freqs_t three_freqs;
std::vector<three_freqs_t> three_freqs;
std::vector< std::tuple<int, int, int> > three_indices; // (ib, if1, if2)
{
// indices of fermionic Matsubara frequencies [-n_f:n_f)
Expand All @@ -566,8 +567,8 @@ namespace pomerol2triqs {
}

// compute g2 values
// g2_iw_freq_vec_t g2_three_freqs = compute_g2(p.gf_struct, p.beta, p.channel, p.index1, p.index2, p.index3, p.index4, three_freqs);
std::vector<g2_iw_freq_vec_t> vec_g2_freq_vec = compute_g2_indices_loop(p.gf_struct, p.beta, p.channel, p.vec_four_indices, three_freqs);
// g2_iw_freq_fix_t g2_three_freqs = compute_g2(p.gf_struct, p.beta, p.channel, p.index1, p.index2, p.index3, p.index4, three_freqs);
std::vector<g2_iw_freq_fix_t> vec_g2_freq_vec = compute_g2(p.gf_struct, p.beta, p.channel, p.four_indices, three_freqs);

// reshape G2 (from freq_vec to freq_box)
std::vector<g2_iw_freq_box_t> vec_g2_freq_box;
Expand All @@ -577,7 +578,7 @@ namespace pomerol2triqs {
int ib = std::get<0>(three_indices[i]);
int if1 = std::get<1>(three_indices[i]);
int if2 = std::get<2>(three_indices[i]);
g2(ib, if1, if2) = g2_freq_vec[i];
g2(ib, if1, if2) = g2_freq_vec(i);
}
vec_g2_freq_box.push_back(g2);
}
Expand All @@ -586,11 +587,8 @@ namespace pomerol2triqs {
}


// auto pomerol_ed::G2_iw_three_freqs(g2_three_freqs_params_t const &p) -> g2_iw_freq_vec_t {
// return compute_g2(p.gf_struct, p.beta, p.channel, p.index1, p.index2, p.index3, p.index4, p.three_freqs);
// }
auto pomerol_ed::G2_iw_freqs_vec(g2_iw_freq_vec_params_t const &p) -> std::vector<g2_iw_freq_vec_t> {
return compute_g2_indices_loop(p.gf_struct, p.beta, p.channel, p.vec_four_indices, p.three_freqs);
auto pomerol_ed::G2_iw_freq_fix(g2_iw_freq_fix_params_t const &p) -> std::vector<g2_iw_freq_fix_t> {
return compute_g2(p.gf_struct, p.beta, p.channel, p.four_indices, p.three_freqs);
}
/*
auto pomerol_ed::G2_iw_three_freqs(g2_three_freqs_params_t const &p) -> g2_three_freqs_t {
Expand Down
16 changes: 8 additions & 8 deletions c++/pomerol_ed.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -83,14 +83,14 @@ namespace pomerol2triqs {
// using w_nu_nup_t = cartesian_product<imfreq, imfreq, imfreq>;
// using w_l_lp_t = cartesian_product<imfreq, legendre, legendre>;
using g2_iw_freq_box_t = triqs::arrays::array<std::complex<double>, 3>;
using g2_iw_freq_vec_t = std::vector<std::complex<double> >;
using g2_iw_freq_fix_t = triqs::arrays::array<std::complex<double>, 1>;
// template <typename Mesh, typename Filler>
// block2_gf<Mesh, tensor_valued<4>> compute_g2(gf_struct_t const &gf_struct, gf_mesh<Mesh> const &mesh, block_order_t block_order,
// g2_blocks_t const &g2_blocks, Filler filler) const;
g2_iw_freq_vec_t compute_g2(gf_struct_t const &gf_struct, double beta, channel_t channel, indices_t index1, indices_t index2, indices_t index3, indices_t index4, three_freqs_t const &three_freqs);

// std::vector<g2_iw_freq_vec_t> compute_g2_indices_loop(gf_struct_t const &gf_struct, double beta, channel_t channel, std::vector<four_indices_t> const &vec_four_indices, three_freqs_t const &three_freqs);
std::vector< g2_iw_freq_vec_t > compute_g2_indices_loop(gf_struct_t const &gf_struct, double beta, channel_t channel, std::vector<four_indices_t> const &vec_four_indices, three_freqs_t const &three_freqs);
g2_iw_freq_fix_t compute_g2_core(gf_struct_t const &gf_struct, double beta, channel_t channel, indices_t index1, indices_t index2, indices_t index3, indices_t index4, std::vector<three_freqs_t> const &three_freqs);

std::vector<g2_iw_freq_fix_t> compute_g2(gf_struct_t const &gf_struct, double beta, channel_t channel, std::vector<four_indices_t> const &four_indices, std::vector<three_freqs_t> const &three_freqs);

double density_matrix_cutoff = 1e-15;

Expand Down Expand Up @@ -122,17 +122,17 @@ namespace pomerol2triqs {
/// Retarded Green's function on real energy axis
block_gf<refreq> G_w(gf_struct_t const &gf_struct, double beta, std::pair<double, double> const &energy_window, int n_w, double im_shift = 0);

/// Two-particle Green's function. Specify frequency cutoff, n_b and n_f.
/// Two-particle Green's function, Legacy support
TRIQS_WRAP_ARG_AS_DICT
g2_iw_freq_box_t G2_iw_legacy(g2_iw_legacy_params_t const &p);

/// Two-particle Green's function. Specify frequency cutoff, n_b and n_f.
/// Two-particle Green's function, in a low-frequency box with cutoff n_b and n_f
TRIQS_WRAP_ARG_AS_DICT
std::vector<g2_iw_freq_box_t> G2_iw_freq_box(g2_iw_freq_box_params_t const &p);

/// Two-particle Green's function. Specify three frequencies (wb, wf1, wf2).
/// Two-particle Green's function, for fixed frequencies (wb, wf1, wf2)
TRIQS_WRAP_ARG_AS_DICT
std::vector<g2_iw_freq_vec_t> G2_iw_freqs_vec(g2_iw_freq_vec_params_t const &p);
std::vector<g2_iw_freq_fix_t> G2_iw_freq_fix(g2_iw_freq_fix_params_t const &p);

/// Two-particle Green's function, Matsubara frequencies
// TRIQS_WRAP_ARG_AS_DICT
Expand Down
10 changes: 6 additions & 4 deletions example/2band.atom.py
Original file line number Diff line number Diff line change
Expand Up @@ -96,10 +96,12 @@
###############################

# G2_iw = ed.G2_iw_legacy( index1=('up',0), index2=('dn',0), index3=('dn',1), index4=('up',1), **common_g2_params )
G2_iw = ed.G2_iw_freq_box( vec_four_indices=[(('up',0), ('dn',0), ('dn',1), ('up',1)),], **common_g2_params )[0]
print type(G2_iw)
print G2_iw.shape

four_indices = [(('up',0), ('dn',0), ('dn',1), ('up',1))]
G2_iw = ed.G2_iw_freq_box( four_indices=four_indices, **common_g2_params )

print "G2_iw :", type(G2_iw), "of size", len(G2_iw)
print "G2_iw[0] :", type(G2_iw[0]), "of size", G2_iw[0].shape

# # Compute G^{(2),ph}(i\omega;i\nu,i\nu'), AABB block order
# G2_iw_inu_inup_ph_AABB = ed.G2_iw_inu_inup(channel = "PH",
Expand All @@ -126,4 +128,4 @@
ar['G_iw'] = G_iw
ar['G_tau'] = G_tau
ar['G_w'] = G_w
ar['G2_ph'] = G2_iw
ar['G2_ph'] = G2_iw[0]
31 changes: 24 additions & 7 deletions example/slater.py
Original file line number Diff line number Diff line change
Expand Up @@ -137,19 +137,36 @@

common_g2_params = {'channel' : "PH",
'gf_struct' : gf_struct,
'beta' : beta,
'n_f' : g2_n_wf,
'n_b' : g2_n_wb, }
'beta' : beta,}

###############################
# G^{(2)}(i\omega;i\nu,i\nu') #
###############################

# G2_iw = ed.G2_iw_legacy( index1=('up',0), index2=('dn',0), index3=('dn',1), index4=('up',1), **common_g2_params )
G2_iw = ed.G2_iw_freq_box( vec_four_indices=[(('up',0), ('dn',0), ('dn',1), ('up',1)),], **common_g2_params )[0]

# four-operators indices
four_indices = []
four_indices.append( (('up',0), ('dn',0), ('dn',1), ('up',1)) )
four_indices.append( (('up',0), ('up',0), ('dn',1), ('dn',1)) )

# compute G2 in a low-frequency box
G2_iw = ed.G2_iw_freq_box( four_indices=four_indices, n_f=g2_n_wf, n_b=g2_n_wb, **common_g2_params )

if mpi.is_master_node():
print "G2_iw :", type(G2_iw), "of size", len(G2_iw)
print "G2_iw[0] :", type(G2_iw[0]), "of size", G2_iw[0].shape
with HDFArchive('slater_gf.out.h5', 'a') as ar:
ar['G2_iw_freq_box'] = G2_iw

# compute G2 for given freqs, (wb, wf1, wf2)
three_freqs = []
three_freqs.append( (0, 1, 2) )
three_freqs.append( (0, 1, -2) )
G2_iw = ed.G2_iw_freq_fix( four_indices=four_indices, three_freqs=three_freqs, **common_g2_params )

if mpi.is_master_node():
print type(G2_iw)
print G2_iw.shape
print "G2_iw :", type(G2_iw), "of size", len(G2_iw)
print "G2_iw[0] :", type(G2_iw[0]), "of size", G2_iw[0].shape
with HDFArchive('slater_gf.out.h5', 'a') as ar:
ar['G2_iw'] = G2_iw
ar['G2_iw_freq_fix'] = G2_iw
Loading

0 comments on commit 70721e4

Please sign in to comment.