API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Dec 6th, 2019
143
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 12.72 KB | None | 0 0
raw download clone embed print report
  1. /*****************************************************************************
  2. *
  3. * ALPS MPS DMRG Project
  4. *
  5. * Copyright (C) 2014 Institute for Theoretical Physics, ETH Zurich
  6. *               2011-2011 by Bela Bauer <bauerb@phys.ethz.ch>
  7. *               2011-2013    Michele Dolfi <dolfim@phys.ethz.ch>
  8. *               2014-2014    Sebastian Keller <sebkelle@phys.ethz.ch>
  9. *               2019         Leon Freitag <lefreita@ethz.ch>
  10. *
  11. * This software is part of the ALPS Applications, published under the ALPS
  12. * Application License; you can use, redistribute it and/or modify it under
  13. * the terms of the license, either version 1 or (at your option) any later
  14. * version.
  15. *
  16. * You should have received a copy of the ALPS Application License along with
  17. * the ALPS Applications; see the file LICENSE.txt. If not, the license is also
  18. * available from http://alps.comp-phys.org/.
  19. *
  20. * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  21. * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  22. * FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
  23. * SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
  24. * FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
  25. * ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
  26. * DEALINGS IN THE SOFTWARE.
  27. *
  28. *****************************************************************************/
  29. #ifndef NRDM_XVECTOR_H
  30. #define NRDM_XVECTOR_H
  31.  
  32. // Measurement of transition density matrices for the variation of the MPS tensor
  33. // expressed in nonredundant MPS parameters (XVector)
  34. #include "dmrg/models/measurements/tagged_nrankrdm.h"
  35. #include "dmrg/mp_tensors/xvector.h"
  36. #include "dmrg/models/measurements/measurement_map.h"
  37. #include "dmrg/utils/load_msmps.h"
  38.  
  39.     namespace measurements
  40.     {
  41.         template <class Matrix, class SymmGroup>
  42.         class NRDMXVector : public TaggedNRankRDM<Matrix, SymmGroup>
  43.         {
  44.         typedef TaggedNRankRDM<Matrix, SymmGroup> base;
  45.         typedef typename base::tag_vec tag_vec;
  46.         typedef typename base::scaled_bond_term scaled_bond_term;
  47.         typedef typename base::positions_type positions_type;
  48.         typedef typename base::meas_type meas_type;
  49.  
  50.         using base::ext_labels;
  51.         using base::do_measure_1rdm;
  52.         using base::do_measure_2rdm;
  53.         using base::lattice;
  54.  
  55.         public:
  56.             // Specialization to call the correct constructor of the base class
  57.             // TODO: Check if this compiles with disabled SU2U1/SU2U1PG!
  58.             // symm_traits::HasSU2<SU2U1> and symm_traits::HasSU2<SU2U1PG> yields boost::true_type
  59.             NRDMXVector(boost::true_type, std::string const& chkp, std::string const& aux_filename, std::string const& msmps_names, std::string const& name_,  int canonization_site, const Lattice & lat,
  60.                        boost::shared_ptr<TagHandler<Matrix, SymmGroup> > tag_handler_,
  61.                        typename TermMakerSU2<Matrix, SymmGroup>::OperatorCollection const & op_collection_,
  62.                        positions_type const& positions_ = positions_type(), bool twosite_ = true)
  63.                        : base(name_, lat, tag_handler_, op_collection_, positions_), chkp_(chkp), aux_filename_(aux_filename), msmps_names_(msmps_names), canonization_site_(canonization_site) {};
  64.  
  65.             // 2U1 and other symmetry groups
  66.             NRDMXVector(int lr_site_, boost::false_type, std::string const& chkp, std::string const& aux_filename, std::string const& msmps_names, std::string const& name_, int canonization_site, const Lattice & lat,
  67.                        boost::shared_ptr<TagHandler<Matrix, SymmGroup> > tag_handler_,
  68.                        tag_vec const & identities_, tag_vec const & fillings_, std::vector<scaled_bond_term> const& ops_,
  69.                        bool half_only_, positions_type const& positions_ = positions_type(), bool twosite_ = true)
  70.                        : base(name_, lat, tag_handler_, identities_, fillings_, ops_, half_only_, positions_), chkp_(chkp), aux_filename_(aux_filename), msmps_names_(msmps_names), canonization_site_(canonization_site) {};
  71.  
  72.             virtual void evaluate(MPS<Matrix, SymmGroup> const& ket_mps, boost::optional<reduced_mps<Matrix, SymmGroup> const&> rmps = boost::none)
  73.             {
  74.                 const std::vector<MPS<Matrix, SymmGroup> > & msmps = load_msmps<Matrix, SymmGroup>(msmps_names_);
  75.  
  76.                 lr::XVector<Matrix, SymmGroup> xvec;
  77.                 xvec.load(this->chkp_, canonization_site_);  // Load the block matrix structure
  78.                 xvec.assign_mps(msmps);  // Specify the reference MPS for the xvector
  79.  
  80.                 if (!this->aux_filename_.empty()) // If we have an aux filename specified, load the data from the auxiliary text file (used in the MOLCAS interface)
  81.                     xvec.update_from_textfile(this->aux_filename_);
  82.  
  83.                 if (this->name() == "onerdmxvector")
  84.                     measure_xvector_1rdm(msmps, xvec);
  85.  
  86.                 if (this->name() == "twordmxvector")
  87.                     measure_xvector_2rdm(msmps, xvec);
  88.  
  89.             }
  90.         protected:
  91.             measurement<Matrix, SymmGroup>* do_clone() const
  92.             {
  93.                 return new NRDMXVector(*this);
  94.             }
  95.  
  96.             // Measure transition densities wrt variation of the MPSTensors
  97.             // The transition density wrt the variation of the MPSTensor is the sum of all transition densities
  98.             // for each site, where the corresponding MPSTensor has been replaced by its variation
  99.             void measure_xvector_1rdm(std::vector<MPS<Matrix, SymmGroup> > const & msmps, lr::XVector<Matrix, SymmGroup> const & xvec)
  100.             {
  101.                 onerdm_map<meas_type> rdm_full;
  102.                 MPS<Matrix, SymmGroup> mps_aux = msmps[0];
  103.  
  104.                 for (int site = 0; site < msmps[0].length(); site++)
  105.                 {
  106.                     if (!xvec.has_no_redundant(site))
  107.                     {
  108.                         onerdm_map<meas_type> rdm_partial;
  109.  
  110.                         if (site == canonization_site_)
  111.                         {
  112.                             for (int state = 0; state < msmps.size(); state++)
  113.                             {
  114.                                 assert(msmps[0].length() > 0);
  115.                                 mps_aux[site] = xvec.getB(site,state);
  116.                                 mps_aux[site].make_left_paired(); // required for functioning of expval in parallel
  117.  
  118.                                 bool first_measurement = rdm_full.empty();
  119.                                 onerdm_map<meas_type> & rdm_tomeas = first_measurement ? rdm_full : rdm_partial;
  120.  
  121.                                 rdm_tomeas = do_measure_1rdm(mps_aux, msmps[state]);
  122.  
  123.                                 // divide by SA weight
  124.                                 for (auto&& el: rdm_tomeas)
  125.                                         el.second /= msmps.size();
  126.  
  127.                                 if (!first_measurement)
  128.                                     for (typename onerdm_map<meas_type>::iterator it = rdm_partial.begin(); it != rdm_partial.end(); it++)
  129.                                         rdm_full.at(it->first) += it->second;
  130.  
  131.                                 mps_aux[site] = msmps[state][site];
  132.                             }
  133.                         }
  134.                         else
  135.                         {
  136.                             for (int state = 0; state < msmps.size(); state++)
  137.                                 {
  138.                                 // Measure RDM only once and multiply the result by the number of states
  139.                                 mps_aux[site] = xvec.getB(site,state);
  140.                                 mps_aux[canonization_site_] = msmps[state][canonization_site_];
  141.                                 mps_aux[site].make_left_paired(); // required for functioning of expval in parallel
  142.  
  143.                                 bool first_measurement = rdm_full.empty();
  144.                                 onerdm_map<meas_type> & rdm_tomeas = (first_measurement) ? rdm_full : rdm_partial;
  145.  
  146.                                 rdm_tomeas = do_measure_1rdm(mps_aux, msmps[state]);
  147.  
  148.                                 // divide by SA weight
  149.                                 for (auto&& el: rdm_tomeas)
  150.                                         el.second /= msmps.size();
  151.  
  152.                                 if (!first_measurement)
  153.                                     for (typename onerdm_map<meas_type>::iterator it = rdm_partial.begin(); it != rdm_partial.end(); it++)
  154.                                         rdm_full.at(it->first) += it->second;
  155.  
  156.                                 mps_aux[site] = msmps[state][site];
  157.                             }
  158.                         }
  159.                     }
  160.  
  161.                 }
  162.  
  163.                 // Convert the total map to the old-style results with labels and results
  164.                 std::tie(this->labels, this->vector_results) = measurements_details::convert_meas_map_to_old_results(rdm_full, lattice, ext_labels);
  165.             }
  166.  
  167.             // and the same for 2-rdm (unfortunately copy-pasted from 1-rdm, maybe there's an easy way to make it more generic)
  168.             void measure_xvector_2rdm(std::vector<MPS<Matrix, SymmGroup> > const & msmps, lr::XVector<Matrix, SymmGroup> const & xvec)
  169.             {
  170.                 twordm_map<meas_type> rdm_full;
  171.                 MPS<Matrix, SymmGroup> mps_aux = msmps[0];
  172.  
  173.                 for (int site = 0; site < msmps[0].length(); site++)
  174.                 {
  175.                     if (!xvec.has_no_redundant(site))
  176.                     {
  177.                         twordm_map<meas_type> rdm_partial;
  178.  
  179.                         if (site == canonization_site_)
  180.                         {
  181.                             for (int state = 0; state < msmps.size(); state++)
  182.                             {
  183.                                 assert(msmps[0].length() > 0);
  184.                                 mps_aux[site] = xvec.getB(site,state);
  185.                                 mps_aux[site].make_left_paired(); // required for functioning of expval in parallel
  186.  
  187.                                 bool first_measurement = rdm_full.empty();
  188.                                 twordm_map<meas_type> & rdm_tomeas = (first_measurement) ? rdm_full : rdm_partial;
  189.  
  190.                                 rdm_tomeas = do_measure_2rdm(mps_aux, msmps[state]);
  191.  
  192.                                 for (auto&& el: rdm_tomeas)
  193.                                     el.second /= msmps.size();
  194.  
  195.                                 if (!first_measurement)
  196.                                     for (typename twordm_map<meas_type>::iterator it = rdm_partial.begin(); it != rdm_partial.end(); it++)
  197.                                         rdm_full.at(it->first) += it->second;
  198.  
  199.                                 mps_aux[site] = msmps[state][site];
  200.                             }
  201.                         }
  202.                         else
  203.                         {
  204.                             for (int state = 0; state < msmps.size(); state++)
  205.                             {
  206.                                 // Measure RDM only once and multiply the result by the number of states
  207.                                 mps_aux[site] = xvec.getB(site,state);
  208.                                 mps_aux[canonization_site_] = msmps[state][canonization_site_];
  209.                                 mps_aux[site].make_left_paired(); // required for functioning of expval in parallel
  210.  
  211.                                 bool first_measurement = rdm_full.empty();
  212.                                 twordm_map<meas_type> & rdm_tomeas = (first_measurement) ? rdm_full : rdm_partial;
  213.  
  214.                                 rdm_tomeas = do_measure_2rdm(mps_aux, msmps[state]);
  215.  
  216.                                 // divide by SA weight
  217.                                 for (auto&& el: rdm_tomeas)
  218.                                     el.second /= msmps.size();
  219.  
  220.                                 if (!first_measurement)
  221.                                     for (typename twordm_map<meas_type>::iterator it = rdm_partial.begin(); it != rdm_partial.end(); it++)
  222.                                         rdm_full.at(it->first) += it->second;
  223.  
  224.                                 mps_aux[site] = msmps[state][site];
  225.                             }
  226.                         }
  227.                     }
  228.                 }
  229.  
  230.                 // Convert the total map to the old-style results with labels and results
  231.                 std::tie(this->labels, this->vector_results) = measurements_details::convert_meas_map_to_old_results(rdm_full, lattice, ext_labels);
  232.             }
  233.  
  234.             std::string chkp_;
  235.             std::string aux_filename_;
  236.             std::string msmps_names_;
  237.             int canonization_site_;
  238.         };
  239.     }
  240. #endif
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!