Skip to content

Commit

Permalink
Merge branch 'release-3.1.1' into 'master'
Browse files Browse the repository at this point in the history 
Release 3.1.1

See merge request dmrg/maquis-dmrg!25

 - Fixed "pre existing term" bug
 - Fixed overlap bug between different MPSs
  • Loading branch information
@AlbertoBaiardi
Alberto Baiardi authored and AlbertoBaiardi committed Dec 7, 2021
2 parents 3efd0f5 + 5463f1f commit 72fa28a
Show file tree
Hide file tree
Showing 57 changed files with 2,732 additions and 5,122 deletions.
5 changes: 5 additions & 0 deletions CHANGELOG.md
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,10 @@
# Changelog

## Release 3.1.1

- Fixed bug raising the "pre-existing term" message when constructing the MPO
- Enhancement of the overlap calculation in MPS-SI

## Release 3.1.0

- Added support for TD-DMRG with time-independent Hamiltonian
Expand Down
6 changes: 4 additions & 2 deletions dmrg/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -258,7 +258,7 @@ endif()
set(DMRG_YEAR 2021)
set(DMRG_VERSION_MAJOR 3)
set(DMRG_VERSION_MINOR 1)
set(DMRG_VERSION_PATCH 0)
set(DMRG_VERSION_PATCH 1)
set(DMRG_VERSION_BUILD "")

if(NOT DMRG_VERSION_BUILD AND EXISTS ${PROJECT_SOURCE_DIR}/.svn)
Expand Down Expand Up @@ -393,7 +393,9 @@ if(QCMAQUIS_TESTS)
add_test(NAME Test_Relativistic COMMAND test_rel)
add_test(NAME Test_Integral_Map COMMAND test_integral_map)
add_test(NAME Test_HiRDM COMMAND test_hirdm)
add_test(NAME Test_MPS_MPO_Ops COMMAND test_mps_mpo_ops)
add_test(NAME Test_MPS_MPO_Ops_TwoU1 COMMAND test_mps_mpo_ops_TwoU1)
add_test(NAME Test_MPS_MPO_Ops_Electronic COMMAND test_mps_mpo_ops_electronic)
add_test(NAME Test_MPS_Overlap_Electronic COMMAND test_mps_overlap_electronic)
add_test(NAME Test_SiteProblem COMMAND test_siteproblem)
add_test(NAME Test_MPS_Join COMMAND test_mpsjoin)
add_test(NAME Test_Wigner COMMAND test_wigner)
Expand Down
4 changes: 0 additions & 4 deletions dmrg/framework/dmrg/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,10 +28,6 @@ foreach(SYMM ${BUILD_SYMMETRIES})
if(EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/${TRIAL_INCLUDE_FACTORY}")
add_line_to(MAQUIS_INCLUDE_FACTORIES_${GROUP_NAME} "#include \"dmrg/${TRIAL_INCLUDE_FACTORY}\"")
endif()
set(TRIAL_INCLUDE_FACTORY "models/continuum/factory_${SYMM_SUFFIX}.ipp")
if(EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/${TRIAL_INCLUDE_FACTORY}")
add_line_to(MAQUIS_INCLUDE_FACTORIES_${GROUP_NAME} "#include \"dmrg/${TRIAL_INCLUDE_FACTORY}\"")
endif()
endforeach(SYMM)

configure_symm_file("models/model_factory_symm/model_factory_tpl.cpp.in" "${CMAKE_CURRENT_BINARY_DIR}/models/model_factory_symm/model_factory_{SYMM}.cpp" CMAKE_SYMM_GROUP)
Expand Down
86 changes: 45 additions & 41 deletions dmrg/framework/dmrg/block_matrix/block_matrix_algorithms.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,52 +45,46 @@

#include "dmrg/utils/parallel.hpp"

/** @brief Struct storing the results of an MPS truncation */
struct truncation_results {
std::size_t bond_dimension; // new bond dimension
double truncated_weight; // sum of discarded eigenvalues (square of singuler values)
double truncated_fraction; // sum of discarded singular values
double smallest_ev; // smallest eigenvalue kept

/** @brief Empty constructor */
truncation_results() { }

/** @brief Constructor taking input data */
truncation_results(std::size_t m, double tw, double tf, double se)
: bond_dimension(m)
, truncated_weight(tw)
, truncated_fraction(tf)
, smallest_ev(se)
: bond_dimension(m), truncated_weight(tw), truncated_fraction(tf), smallest_ev(se)
{ }
};

//template<class Matrix1, class Matrix2, class Matrix3, class SymmGroup, class Scheduler = parallel::scheduler_nop>
/** @brief Multiplication between block matrices */
template<class Matrix1, class Matrix2, class Matrix3, class SymmGroup, class Scheduler>
void gemm(block_matrix<Matrix1, SymmGroup> const & A,
block_matrix<Matrix2, SymmGroup> const & B,
block_matrix<Matrix3, SymmGroup> & C,
const Scheduler& scheduler = Scheduler())
void gemm(block_matrix<Matrix1, SymmGroup> const & A, block_matrix<Matrix2, SymmGroup> const & B,
block_matrix<Matrix3, SymmGroup> & C, const Scheduler& scheduler = Scheduler())
{
// Types definition
using charge = typename SymmGroup::charge;
using const_iterator = typename DualIndex<SymmGroup>::const_iterator;
//
C.clear();
assert(B.basis().is_sorted());

typedef typename SymmGroup::charge charge;
typedef typename DualIndex<SymmGroup>::const_iterator const_iterator;
const_iterator B_begin = B.basis().begin();
const_iterator B_end = B.basis().end();
for (std::size_t k = 0; k < A.n_blocks(); ++k) {

for (int k = 0; k < A.n_blocks(); ++k) {
charge ar = A.basis().right_charge(k);
const_iterator it = B.basis().left_lower_bound(ar);

for ( ; it != B_end && it->lc == ar; ++it)
{
std::size_t matched_block = std::distance(B_begin, it);
for ( ; it != B_end && it->lc == ar; ++it) {
auto matched_block = std::distance(B_begin, it);
Matrix3 tmp(num_rows(A[k]), it->rs);

parallel::guard proc(scheduler(k));
gemm(A[k], B[matched_block], tmp);
C.match_and_add_block(tmp, A.basis().left_charge(k), it->rc);
}
}

if(scheduler.propagate()){
Index<SymmGroup> B_left_basis = B.left_basis();
C.size_index.resize(C.n_blocks()); // propagating A size_index onto C - otherwise might C.index_sizes();
Expand All @@ -102,6 +96,7 @@ void gemm(block_matrix<Matrix1, SymmGroup> const & A,
}
}

/** @brief Wrapper around gemm that avoids setting the scheduler */
template<class Matrix1, class Matrix2, class Matrix3, class SymmGroup>
void gemm(block_matrix<Matrix1, SymmGroup> const & A,
block_matrix<Matrix2, SymmGroup> const & B,
Expand All @@ -110,8 +105,18 @@ void gemm(block_matrix<Matrix1, SymmGroup> const & A,
gemm(A, B, C, parallel::scheduler_nop());
}

// ref_left_basis: basis of B equivalent in the bra for contractions where bra != ket
// for bra == ket it is equal to the left basis of B
/**
* @brief Matrix-matrix multiplication with left basis trimming.
*
* In addition to performing the matrix-matrix multiplication, this routine
* also "trims" the row index by keeping only the blocks that are present
* in [ref_left_basis].
*
* @param A First input matrix
* @param B Second input matrix
* @param C Output matrix
* @param ref_left_basis Index used as a reference for trimming
*/
template<class Matrix1, class Matrix2, class Matrix3, class SymmGroup>
void gemm_trim_left(block_matrix<Matrix1, SymmGroup> const & A,
block_matrix<Matrix2, SymmGroup> const & B,
Expand All @@ -120,29 +125,33 @@ void gemm_trim_left(block_matrix<Matrix1, SymmGroup> const & A,
{
parallel::scheduler_size_indexed scheduler(A);
C.clear();

typedef typename SymmGroup::charge charge;
Index<SymmGroup> B_left_basis = B.left_basis();
for (std::size_t k = 0; k < A.n_blocks(); ++k) {
std::size_t matched_block = B_left_basis.position(A.basis().right_charge(k));

auto matched_block = B_left_basis.position(A.basis().right_charge(k));
// Match right basis of A with left basis of B
if ( matched_block == B.n_blocks() )
continue;

if ( !ref_left_basis.has(A.basis().left_charge(k)) )
continue;

std::size_t new_block = C.insert_block(new Matrix3(num_rows(A[k]), num_cols(B[matched_block])),
A.basis().left_charge(k), B.basis().right_charge(matched_block));

auto new_block = C.insert_block(new Matrix3(num_rows(A[k]), num_cols(B[matched_block])),
A.basis().left_charge(k), B.basis().right_charge(matched_block));
parallel::guard proc(scheduler(k));
gemm(A[k], B[matched_block], C[new_block]);
}
}

// ref_right_basis: basis of A equivalent in the bra for contractions where bra != ket
// for bra == ket it is equal to the right basis of A
/**
* @brief Matrix-matrix multiplication with right basis trimming.
*
* In addition to performing the matrix-matrix multiplication, this routine
* also "trims" the column index by keeping only the blocks that are present
* in [ref_right_basis].
*
* @param A First input matrix
* @param B Second input matrix
* @param C Output matrix
* @param ref_left_basis Index used as a reference for trimming
*/
template<class Matrix1, class Matrix2, class Matrix3, class SymmGroup>
void gemm_trim_right(block_matrix<Matrix1, SymmGroup> const & A,
block_matrix<Matrix2, SymmGroup> const & B,
Expand All @@ -154,21 +163,16 @@ void gemm_trim_right(block_matrix<Matrix1, SymmGroup> const & A,

typedef typename SymmGroup::charge charge;
Index<SymmGroup> A_right_basis = A.right_basis();
for (std::size_t k = 0; k < B.n_blocks(); ++k) {
std::size_t matched_block = A_right_basis.position(B.basis().left_charge(k));

for (int k = 0; k < B.n_blocks(); ++k) {
auto matched_block = A_right_basis.position(B.basis().left_charge(k));
// Match right basis of A with left basis of B
if ( matched_block == A.n_blocks() )
continue;

// Also match right_basis of bra with B.right_basis()

if ( !ref_right_basis.has(B.basis().right_charge(k)) )
continue;

std::size_t new_block = C.insert_block(new Matrix3(num_rows(A[matched_block]), num_cols(B[k])),
A.basis().left_charge(matched_block), B.basis().right_charge(k));

A.basis().left_charge(matched_block), B.basis().right_charge(k));
parallel::guard proc(scheduler(k));
gemm(A[matched_block], B[k], C[new_block]);
}
Expand Down
4 changes: 2 additions & 2 deletions dmrg/framework/dmrg/models/chem/2u1/chem_helper.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,8 +101,8 @@ namespace detail {
term_descriptor term;
term.is_fermionic = false;
term.coeff = scale * ptag1.second * ptag2.second;
term.push_back(boost::make_tuple(p1, ptag1.first));
term.push_back(boost::make_tuple(p2, ptag2.first));
term.push_back(std::make_pair(p1, ptag1.first));
term.push_back(std::make_pair(p2, ptag2.first));

IndexTuple id(p1, p2, ptag1.first, ptag2.first);

Expand Down
8 changes: 4 additions & 4 deletions dmrg/framework/dmrg/models/chem/2u1/model.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -211,7 +211,7 @@ void qc_model<Matrix, SymmGroup>::create_terms()

term_descriptor term;
term.coeff = matrix_elements[m];
term.push_back( boost::make_tuple(0, ident[lat.get_prop<typename SymmGroup::subcharge>("type", 0)]));
term.push_back( std::make_pair(0, ident[lat.get_prop<typename SymmGroup::subcharge>("type", 0)]));
this->terms_.push_back(term);

used_elements[m] += 1;
Expand All @@ -222,13 +222,13 @@ void qc_model<Matrix, SymmGroup>::create_terms()
{
term_descriptor term;
term.coeff = matrix_elements[m];
term.push_back( boost::make_tuple(i, count_up[lat.get_prop<typename SymmGroup::subcharge>("type", i)]));
term.push_back( std::make_pair(i, count_up[lat.get_prop<typename SymmGroup::subcharge>("type", i)]));
this->terms_.push_back(term);
}
{
term_descriptor term;
term.coeff = matrix_elements[m];
term.push_back( boost::make_tuple(i, count_down[lat.get_prop<typename SymmGroup::subcharge>("type", i)]));
term.push_back( std::make_pair(i, count_down[lat.get_prop<typename SymmGroup::subcharge>("type", i)]));
this->terms_.push_back(term);
}

Expand Down Expand Up @@ -260,7 +260,7 @@ void qc_model<Matrix, SymmGroup>::create_terms()

term_descriptor term;
term.coeff = matrix_elements[m];
term.push_back(boost::make_tuple(i, docc[lat.get_prop<typename SymmGroup::subcharge>("type", 0)]));
term.push_back(std::make_pair(i, docc[lat.get_prop<typename SymmGroup::subcharge>("type", 0)]));
this->terms_.push_back(term);

used_elements[m] += 1;
Expand Down
48 changes: 23 additions & 25 deletions dmrg/framework/dmrg/models/chem/2u1/term_maker.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,6 @@
* Copyright (C) 2015 Laboratory for Physical Chemistry, ETH Zurich
* 2012-2015 by Sebastian Keller <[email protected]>
*
*
* This software is part of the ALPS Applications, published under the ALPS
* Application License; you can use, redistribute it and/or modify it under
* the terms of the license, either version 1 or (at your option) any later
Expand Down Expand Up @@ -39,10 +38,9 @@ struct TermMaker {
typedef typename term_descriptor::value_type pos_op_t;
typedef typename S::subcharge sc_t;

static bool compare_tag(pos_op_t p1,
pos_op_t p2)
static bool compare_tag(const pos_op_t& p1, const pos_op_t& p2)
{
return boost::tuples::get<0>(p1) < boost::tuples::get<0>(p2);
return p1.first < p2.first;
}

static term_descriptor two_term(bool sign, std::vector<tag_type> const & fill_op, value_type scale, pos_t i, pos_t j,
Expand All @@ -53,8 +51,8 @@ struct TermMaker {
term_descriptor term;
term.is_fermionic = sign;
term.coeff = scale;
term.push_back(boost::make_tuple(i, op1[lat.get_prop<sc_t>("type", i)]));
term.push_back(boost::make_tuple(j, op2[lat.get_prop<sc_t>("type", j)]));
term.push_back(std::make_pair(i, op1[lat.get_prop<sc_t>("type", i)]));
term.push_back(std::make_pair(j, op2[lat.get_prop<sc_t>("type", j)]));
return term;
}

Expand All @@ -70,14 +68,14 @@ struct TermMaker {
std::pair<tag_type, value_type> ptag;
if (i < j) {
ptag = op_table->get_product_tag(fill_op[lat.get_prop<sc_t>("type", i)], op1[lat.get_prop<sc_t>("type", i)]);
term.push_back(boost::make_tuple(i, ptag.first));
term.push_back(boost::make_tuple(j, op2[lat.get_prop<sc_t>("type", j)]));
term.push_back(std::make_pair(i, ptag.first));
term.push_back(std::make_pair(j, op2[lat.get_prop<sc_t>("type", j)]));
term.coeff *= ptag.second;
}
else {
ptag = op_table->get_product_tag(fill_op[lat.get_prop<sc_t>("type", j)], op2[lat.get_prop<sc_t>("type", j)]);
term.push_back(boost::make_tuple(i, op1[lat.get_prop<sc_t>("type", i)]));
term.push_back(boost::make_tuple(j, ptag.first));
term.push_back(std::make_pair(i, op1[lat.get_prop<sc_t>("type", i)]));
term.push_back(std::make_pair(j, ptag.first));
term.coeff *= -ptag.second;
}
return term;
Expand All @@ -99,14 +97,14 @@ struct TermMaker {
std::pair<tag_type, value_type> ptag;
if (i < j) {
ptag = op_table->get_product_tag(fill_op[lat.get_prop<sc_t>("type", i)], pre_ptag.first);
term.push_back(boost::make_tuple(i, ptag.first));
term.push_back(boost::make_tuple(j, op3[lat.get_prop<sc_t>("type", j)]));
term.push_back(std::make_pair(i, ptag.first));
term.push_back(std::make_pair(j, op3[lat.get_prop<sc_t>("type", j)]));
term.coeff *= ptag.second * pre_ptag.second;
}
else {
ptag = op_table->get_product_tag(fill_op[lat.get_prop<sc_t>("type", j)], op3[lat.get_prop<sc_t>("type", j)]);
term.push_back(boost::make_tuple(i, pre_ptag.first));
term.push_back(boost::make_tuple(j, ptag.first));
term.push_back(std::make_pair(i, pre_ptag.first));
term.push_back(std::make_pair(j, ptag.first));
term.coeff *= -ptag.second * pre_ptag.second;
}
return term;
Expand Down Expand Up @@ -155,9 +153,9 @@ struct TermMaker {
}

std::vector<pos_op_t> sterm;
sterm.push_back( boost::make_tuple(pb, boson_op) );
sterm.push_back( boost::make_tuple(p1, op1) );
sterm.push_back( boost::make_tuple(p2, op2) );
sterm.push_back( std::make_pair(pb, boson_op) );
sterm.push_back( std::make_pair(p1, op1) );
sterm.push_back( std::make_pair(p2, op2) );
std::sort(sterm.begin(), sterm.end(), compare_tag);

term.push_back(sterm[0]);
Expand Down Expand Up @@ -186,18 +184,18 @@ struct TermMaker {
if(idx[c1] > idx[c2]) inv_count++;

std::vector<pos_op_t> sterm;
sterm.push_back(boost::make_tuple(i, op_i[lat.get_prop<sc_t>("type", i)]));
sterm.push_back(boost::make_tuple(j, op_j[lat.get_prop<sc_t>("type", j)]));
sterm.push_back(boost::make_tuple(k, op_k[lat.get_prop<sc_t>("type", k)]));
sterm.push_back(boost::make_tuple(l, op_l[lat.get_prop<sc_t>("type", l)]));
sterm.push_back(std::make_pair(i, op_i[lat.get_prop<sc_t>("type", i)]));
sterm.push_back(std::make_pair(j, op_j[lat.get_prop<sc_t>("type", j)]));
sterm.push_back(std::make_pair(k, op_k[lat.get_prop<sc_t>("type", k)]));
sterm.push_back(std::make_pair(l, op_l[lat.get_prop<sc_t>("type", l)]));
std::sort(sterm.begin(), sterm.end(), compare_tag);

std::pair<tag_type, value_type> ptag;
ptag = op_table->get_product_tag(fill_op[lat.get_prop<sc_t>("type", boost::tuples::get<0>(sterm[0]))], boost::tuples::get<1>(sterm[0]));
boost::tuples::get<1>(sterm[0]) = ptag.first;
ptag = op_table->get_product_tag(fill_op[lat.get_prop<sc_t>("type", sterm[0].first)], sterm[0].second);
sterm[0].second = ptag.first;
term.coeff *= ptag.second;
ptag = op_table->get_product_tag(fill_op[lat.get_prop<sc_t>("type", boost::tuples::get<0>(sterm[2]))], boost::tuples::get<1>(sterm[2]));
boost::tuples::get<1>(sterm[2]) = ptag.first;
ptag = op_table->get_product_tag(fill_op[lat.get_prop<sc_t>("type", sterm[2].first)], sterm[2].second);
sterm[2].second = ptag.first;
term.coeff *= ptag.second;

if (inv_count % 2)
Expand Down
4 changes: 2 additions & 2 deletions dmrg/framework/dmrg/models/chem/rel/model.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -114,7 +114,7 @@ void rel_qc_model<Matrix, SymmGroup>::create_terms()

term_descriptor term;
term.coeff = matrix_elements[m];
term.push_back( boost::make_tuple(0, ident[lat.get_prop<typename SymmGroup::subcharge>("type", 0)]) );
term.push_back( std::make_pair(0, ident[lat.get_prop<typename SymmGroup::subcharge>("type", 0)]) );
this->terms_.push_back(term);

used_elements[m] += 1;
Expand All @@ -127,7 +127,7 @@ void rel_qc_model<Matrix, SymmGroup>::create_terms()
{
term_descriptor term;
term.coeff = matrix_elements[m];
term.push_back( boost::make_tuple(i, count[lat.get_prop<typename SymmGroup::subcharge>("type", i)]));
term.push_back( std::make_pair(i, count[lat.get_prop<typename SymmGroup::subcharge>("type", i)]));
this->terms_.push_back(term);
}
used_elements[m] += 1;
Expand Down
Loading

0 comments on commit 72fa28a

Please sign in to comment.