Guest User

Untitled

a guest
Dec 6th, 2019
114
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  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
RAW Paste Data