lattice: clean old classes, add bz_mesh, gf on bz

pytriqs/lattice/lattice_tools.pyx

@@ -10,19 +10,17 @@ from cython.operator cimport dereference as deref, preincrement as inc #derefere
 
 import numpy
 
-cdef extern from "triqs/lattice/bravais_lattice_and_brillouin_zone.hpp" namespace "triqs::gf" : 
+cdef extern from "triqs/lattice/brillouin_zone.hpp" : 
-    cdef cppclass bravais_lattice_c " triqs::lattice_tools::bravais_lattice" : 
+    cdef cppclass bravais_lattice_c " triqs::lattice::bravais_lattice" : 
-        bravais_lattice_c(matrix[double] & units)
+        bravais_lattice_c(matrix[double] & units, vector[tqa_vector[double]] &,  vector[std_string] &) except +
         int n_orbitals()
         matrix[double] units()
         int dim()
-        void push_back(std_string &, tqa_vector[double] &)
         tqa_vector[double] lattice_to_real_coordinates (tqa_vector_view[double] &)
 
 cdef extern from "triqs/lattice/tight_binding.hpp" namespace "triqs::gf" : 
-    cdef cppclass tight_binding " triqs::lattice_tools::tight_binding" : 
+    cdef cppclass tight_binding " triqs::lattice::tight_binding" : 
-        tight_binding(bravais_lattice_c & units)
+        tight_binding(bravais_lattice_c & units, vector[vector[long]] &, vector[matrix[dcomplex]] &) except +
-        void push_back (vector[long] &, matrix[double] &) 
 
     #Fix the name conflict pv wiht a pxd, cf doc....?
     array_view[complex,THREE] hopping_stack_c "hopping_stack" (tight_binding  & TB, array_view[double,TWO] & k_stack)
@@ -37,10 +35,13 @@ cdef class BravaisLattice :
     cdef bravais_lattice_c * _c
     def __init__(self, units, orbital_positions = None ) :
         """  """
-        self._c = new bravais_lattice_c( matrix[double](units))
+        cdef vector[std_string] names
+        cdef vector[tqa_vector[double]] atom_pos 
         orbital_positions  = orbital_positions if orbital_positions else dict ("", (0,0,0) )
         for name, pos in orbital_positions.items():
-            self._c.push_back( name, tqa_vector[double](pos))
+            names.push_back(name)
+            atom_pos.push_back(tqa_vector[double](pos))
+        self._c = new bravais_lattice_c( matrix[double](units),atom_pos, names)
 
     def __dealloc__(self):
         del self._c
@@ -69,11 +70,16 @@ cdef class TightBinding:
     cdef tight_binding * _c
     def __init__(self, BravaisLattice bravais_lattice, hopping ) : 
         """  """
-        self._c = new tight_binding(deref(bravais_lattice._c))
+        #self._c = new tight_binding(deref(bravais_lattice._c))
-        cdef vector[long] d
+        #cdef vector[long] d
+        cdef vector[matrix[dcomplex]] mats
+        cdef vector[vector[long]] displs
         for displ, mat in hopping.items() :
-            d = displ
+            displs.push_back(displ)
-            self._c.push_back( d,  matrix[double] (numpy.array(mat, dtype =float)))
+            #d = displ
+            mats.push_back(matrix[dcomplex] (numpy.array(mat, dtype =float)))
+            #self._c.push_back( d,  matrix[double] (numpy.array(mat, dtype =float)))
+        self._c = new tight_binding(deref(bravais_lattice._c), displs, mats)
 
     def __dealloc__(self):
         del self._c

test/triqs/gfs/g_k_om.cpp

@@ -0,0 +1,36 @@
+#define TRIQS_ARRAYS_ENFORCE_BOUNDCHECK
+#include <triqs/gfs.hpp>
+#include <triqs/gfs/bz.hpp>
+
+using namespace triqs::gfs;
+using namespace triqs::clef;
+using namespace triqs::arrays;
+using namespace triqs::lattice;
+
+#define TEST(X) std::cout << BOOST_PP_STRINGIZE((X)) << " ---> " << (X) << std::endl << std::endl;
+
+int main() {
+ try {
+  double beta = 1;
+
+  auto bz_ = brillouin_zone{bravais_lattice{make_unit_matrix<double>(2)}};
+
+  auto g_eps = gf<bz>{{bz_, 1000}, {1, 1}};
+
+  auto G = gf<cartesian_product<bz, imfreq>>{{{bz_, 100}, {beta, Fermion, 100}}, {1, 1}};
+
+  // try to assign to expression 
+  placeholder<0> k_;
+  placeholder<1> w_;
+  auto eps = make_expr( [](k_t const& k) { return -2 * (cos(k(0)) + cos(k(1))); }) ;
+
+  G(k_, w_) << 1 / (w_ - eps(k_) - 1 / (w_ + 2));
+
+  // hdf5 
+  H5::H5File file("ess_g_k_om.h5", H5F_ACC_TRUNC );
+  h5_write(file, "g", G);
+
+ 
+ }
+ TRIQS_CATCH_AND_ABORT;
+}

triqs/gfs/bz.hpp

@@ -0,0 +1,75 @@
+/*******************************************************************************
+ *
+ * TRIQS: a Toolbox for Research in Interacting Quantum Systems
+ *
+ * Copyright (C) 2012-2013 by O. Parcollet
+ *
+ * TRIQS is free software: you can redistribute it and/or modify it under the
+ * terms of the GNU General Public License as published by the Free Software
+ * Foundation, either version 3 of the License, or (at your option) any later
+ * version.
+ *
+ * TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY
+ * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+ * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+ * details.
+ *
+ * You should have received a copy of the GNU General Public License along with
+ * TRIQS. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#pragma once
+#include "./tools.hpp"
+#include "./gf.hpp"
+#include "./local/tail.hpp"
+#include "./domains/R.hpp"
+#include "../lattice/bz_mesh.hpp"
+#include "./evaluators.hpp"
+
+namespace triqs {
+namespace gfs {
+
+ struct bz {};
+
+ template <typename Opt> struct gf_mesh<bz, Opt> : lattice::bz_mesh {
+  template <typename... T> gf_mesh(T &&... x) : lattice::bz_mesh(std::forward<T>(x)...) {}
+ };
+
+ namespace gfs_implementation {
+
+  // h5 name
+  template <typename Opt> struct h5_name<bz, matrix_valued, Opt> {
+   static std::string invoke() { return "BZ"; }
+  };
+
+  /// ---------------------------  data access  ---------------------------------
+  template <typename Opt> struct data_proxy<bz, matrix_valued, Opt> : data_proxy_array<std::complex<double>, 3> {};
+  template <typename Opt> struct data_proxy<bz, scalar_valued, Opt> : data_proxy_array<std::complex<double>, 1> {};
+
+  /// ---------------------------  evaluator ---------------------------------
+
+  // simple evaluation : take the point on the grid...
+  template <> struct evaluator_fnt_on_mesh<bz> {
+   lattice::k_t k;
+   evaluator_fnt_on_mesh() = default;
+   template <typename MeshType> evaluator_fnt_on_mesh(MeshType const &m, lattice::k_t const &k_) { k = k_; }
+   template <typename F> auto operator()(F const &f) const DECL_AND_RETURN(f(k));
+  };
+
+  // ------------- evaluator  -------------------
+  // handle the case where the matsu. freq is out of grid...
+  template <typename Target, typename Opt> struct evaluator<bz, Target, Opt> {
+   static constexpr int arity = 1;
+
+   template <typename G>
+   std::c14::conditional_t<std::is_same<Target, matrix_valued>::value, arrays::matrix_const_view<std::complex<double>>,
+                          std::complex<double>>
+   operator()(G const *g, lattice::k_t const &k) const {
+    auto n = g->mesh().locate_neighbours(k); // TO BE IMPROVED
+    return (*g)[n];
+   }
+  };
+ }
+}
+}
+

triqs/lattice/bravais_lattice.hpp

@@ -0,0 +1,82 @@
+/*******************************************************************************
+ *
+ * TRIQS: a Toolbox for Research in Interacting Quantum Systems
+ *
+ * Copyright (C) 2011-2014 by M. Ferrero, O. Parcollet
+ *
+ * TRIQS is free software: you can redistribute it and/or modify it under the
+ * terms of the GNU General Public License as published by the Free Software
+ * Foundation, either version 3 of the License, or (at your option) any later
+ * version.
+ *
+ * TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY
+ * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+ * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+ * details.
+ *
+ * You should have received a copy of the GNU General Public License along with
+ * TRIQS. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#pragma once
+#include <triqs/arrays.hpp>
+#include <string>
+#include <vector>
+
+namespace triqs {
+namespace lattice {
+
+ using r_t = arrays::vector<double>;
+ using k_t = arrays::vector<double>;
+ using dcomplex = std::complex<double>;
+ using arrays::matrix;
+ using arrays::array;
+ using arrays::range;
+
+ // --------------------------------------------------------
+
+ class bravais_lattice {
+
+  public:
+  using point_t = std::vector<int>; // domain concept. PUT on STACK
+
+  bravais_lattice(matrix<double> const& units, std::vector<r_t> atom_orb_pos, std::vector<std::string> atom_orb_name);
+  bravais_lattice(matrix<double> const& units, std::vector<r_t> atom_orb_pos)
+     : bravais_lattice(units, atom_orb_pos, std::vector<std::string>(atom_orb_pos.size(), "")) {}
+  bravais_lattice(matrix<double> const& units) : bravais_lattice(units, std::vector<r_t>{{0, 0, 0}}) {}
+  bravais_lattice() : bravais_lattice(arrays::make_unit_matrix<double>(3)) {}
+
+  int n_orbitals() const { return atom_orb_name.size(); }
+  arrays::matrix_const_view<double> units() const { return units_; }
+  int dim() const { return dim_; }
+
+  /// Transform into real coordinates.
+  template <typename R> r_t lattice_to_real_coordinates(R const& x) const {
+   r_t res(3);
+   res() = 0;
+   for (int i = 0; i < dim_; i++) res += x(i) * units_(i, range{});
+   return res;
+  }
+
+  /// Write into HDF5
+  friend void h5_write(h5::group fg, std::string subgroup_name, bravais_lattice const& bl);
+
+  /// Read from HDF5
+  friend void h5_read(h5::group fg, std::string subgroup_name, bravais_lattice& bl);
+
+  //  BOOST Serialization
+  friend class boost::serialization::access;
+  template <class Archive> void serialize(Archive& ar, const unsigned int version) {
+   ar& boost::serialization::make_nvp("units", units_);
+   ar& boost::serialization::make_nvp("atom_orb_pos", atom_orb_pos);
+   ar& boost::serialization::make_nvp("atom_orb_name", atom_orb_name);
+  }
+
+  private:
+  matrix<double> units_;
+  std::vector<r_t> atom_orb_pos;          // atom_orb_pos[i] = position of ith atoms/orbitals in the unit cell
+  std::vector<std::string> atom_orb_name; // names of these atoms/orbitals.
+  int dim_;
+ };
+}
+}

triqs/lattice/bravais_lattice_and_brillouin_zone.cpp

@@ -18,79 +18,106 @@
  * TRIQS. If not, see <http://www.gnu.org/licenses/>.
  *
  ******************************************************************************/
-#include "bravais_lattice_and_brillouin_zone.hpp"
+#include "./bravais_lattice.hpp"
+#include "./brillouin_zone.hpp"
 #include <triqs/arrays/blas_lapack/dot.hpp>
 #include <triqs/arrays/linalg/det_and_inverse.hpp>
 #include <triqs/arrays/linalg/cross_product.hpp>
-namespace triqs { namespace lattice_tools { 
+namespace triqs {
+namespace lattice {
 
- using namespace tqa;
+ const double almost_zero = 1e-10;
- using namespace std;
- //using triqs::arrays::blas::dot;
- const double almost_zero(1E-10);
 
- bravais_lattice::bravais_lattice( units_type const & units__) : units_(3,3), dim_(units__.shape(0)) { 
+ bravais_lattice::bravais_lattice(matrix<double> const& units__, std::vector<r_t> atom_orb_pos_,
-  units_(range(0,dim_),range()) = units__();
+                                  std::vector<std::string> atom_orb_name_)
-  units_(range(dim_,3),range()) = 0;
+    : units_(3, 3), atom_orb_pos(atom_orb_pos_), atom_orb_name(atom_orb_name_) {
+  dim_ = first_dim(units__);
+  if ((dim_ < 1) || (dim_ > 3)) TRIQS_RUNTIME_ERROR << " units matrix must be square matrix of size 1, 2 or 3";
+  using arrays::range;
+  auto r = range(0, dim_);
+  units_() = 0;
+  units_(r, r) = units__(r, r);
   // First complete the basis. Add some tests for safety
-  tqa::vector<double> ux(3),uy(3),uz(3); 
+  arrays::vector<double> ux(3), uy(3), uz(3);
-  assert(dim_==2);
+  double delta;
+  auto _ = range{};
   switch (dim_) {
    case 1:
-    ux = units_(0,range());
+    ux = units_(0, _);
-    uz() = 0; uz(1) = 1 ;
+    uz() = 0;
+    uz(1) = 1;
     uz = uz - dot(uz, ux) * ux;
     // no luck, ux was parallel to z, another one must work
     if (sqrt(dot(uz, uz)) < almost_zero) {
-     uz() = 0; uz(2) = 1; // 0,0,1;
+     uz() = 0;
+     uz(2) = 1; // 0,0,1;
      uz = uz - dot(uz, ux) * ux;
     }
     uz /= sqrt(dot(uz, uz));
     uy = cross_product(uz, ux);
     uy = uy / sqrt(dot(uy, uy)); // uy can not be 0
-    units_(1,range()) = uz;
+    units_(1, _) = uz;
-    units_(2,range()) = uy;
+    units_(2, _) = uy;
     break;
    case 2:
-    uy() = 0; uy(2) = 1 ;
+    uy() = 0;
-    uy = cross_product(units_(0,range()),units_(1,range()));
+    uy(2) = 1;
-    double delta = sqrt(dot(uy,uy));
+    uy = cross_product(units_(0, _), units_(1, _));
-    if (abs(delta)<almost_zero) TRIQS_RUNTIME_ERROR<<"Tiling : the 2 vectors of unit are not independent";
+    delta = sqrt(dot(uy, uy));
-    units_(2,range()) = uy /delta;
+    if (abs(delta) < almost_zero) TRIQS_RUNTIME_ERROR << "Bravais Lattice : the 2 vectors of unit are not independent : " << units__;
+    units_(2, _) = uy / delta;
+    break;
+   case 3:
+    TRIQS_RUNTIME_ERROR << " 3d bravais lattice not implemented";
+    break;
   }
- //cerr<<" Units = "<< units_<<endl; 
+ }
+ //------------------------------------------------------------------------------------
+
+ /// Write into HDF5
+ void h5_write(h5::group fg, std::string subgroup_name, bravais_lattice const& bl) {
+  h5::group gr = fg.create_group(subgroup_name);
+  h5_write(gr, "units", bl.units_); // NOT COMPLETE
+ }
+
+ /// Read from HDF5
+ void h5_read(h5::group fg, std::string subgroup_name, bravais_lattice& bl) {
+  h5::group gr = fg.open_group(subgroup_name);
+  matrix<double> u;
+  h5_read(gr, "units", u);
+  bl = bravais_lattice{u}; // NOT COMPLETE
  }
 
+ //------------------------------------------------------------------------------------
  //------------------------------------------------------------------------------------
 
  brillouin_zone::brillouin_zone(bravais_lattice const& bl_) : lattice_(bl_), K_reciprocal(3, 3) {
-  bravais_lattice::units_type Units(lattice().units());
+  using arrays::range;
-    std::cout << Units << std::endl;
+  auto _ = range{};
-  double delta = dot(Units(0,range()), cross_product(Units(1,range()),Units(2,range())));
+  auto Units = lattice().units();
-    std::cout << dot(Units(0,range()), cross_product(Units(1,range()),Units(2,range())))<<std::endl;
+  double delta = dot(Units(0, _), cross_product(Units(1, _), Units(2, _)));
-    std::cout << cross_product(Units(1,range()),Units(2,range()))<<std::endl;
+  if (abs(delta) < almost_zero) TRIQS_RUNTIME_ERROR << "Brillouin Zone : the 3 vectors of Units are not independant" << Units;
-  if (abs(delta)<almost_zero) TRIQS_RUNTIME_ERROR<<"Tiling : the 3 vectors of Units are not independant";
+  K_reciprocal(0, _) = cross_product(Units(1, _), Units(2, _)) / delta;
-  K_reciprocal(0,range()) = cross_product(Units(1,range()),Units(2,range())) / delta;
+  K_reciprocal(1, _) = cross_product(Units(2, _), Units(0, _)) / delta;
-  K_reciprocal(1,range()) = cross_product(Units(2,range()),Units(0,range())) / delta;
+  K_reciprocal(2, _) = cross_product(Units(0, _), Units(1, _)) / delta;
-  K_reciprocal(2,range()) = cross_product(Units(0,range()),Units(1,range())) / delta;
+  K_reciprocal = K_reciprocal * 2 * M_PI;
-  //for (size_t i =0; i< lattice().dim();i++) std::cerr << " K_reciprocal(" << i << ")/(2pi) =  " << K_reciprocal(i,range())<< std::endl;
-  const double pi = acos(-1.0);
-  K_reciprocal = K_reciprocal*2*pi;
   K_reciprocal_inv = inverse(K_reciprocal);
  }
 
- K_view_type brillouin_zone::lattice_to_real_coordinates (K_view_type const & k) const { 
+ //------------------------------------------------------------------------------------
-  if (k.size()!=lattice().dim()) TRIQS_RUNTIME_ERROR<<"latt_to_real_k : dimension of k must be "<<lattice().dim();
+ /// Write into HDF5
-  K_type res(3); res()=0; int dim = lattice().dim();
+ void h5_write(h5::group fg, std::string subgroup_name, brillouin_zone const& bz) {
-  for (int i =0; i< dim;i++) res += k (i) * K_reciprocal(i,range());
+  h5::group gr = fg.create_group(subgroup_name);
-  return(res);
+  h5_write(gr, "bravais_lattice", bz.lattice_);
  }
 
- K_view_type brillouin_zone::real_to_lattice_coordinates (K_view_type const & k) const {
+ /// Read from HDF5
-  if (k.size()!=lattice().dim()) TRIQS_RUNTIME_ERROR<<"latt_to_real_k : dimension of k must be "<<lattice().dim();
+ void h5_read(h5::group fg, std::string subgroup_name, brillouin_zone& bz) {
-  K_type res(3);res()=0; int dim = lattice().dim();
+  h5::group gr = fg.open_group(subgroup_name);
-  for (int i =0; i< dim;i++) res += k (i) * K_reciprocal_inv(i,range());
+  bravais_lattice bl;
-  return(res);
+  h5_read(gr, "bravais_lattice", bl);
+  bz = brillouin_zone{bl};
  }
-}}//namespaces
+}
+} // namespaces
 

triqs/lattice/bravais_lattice_and_brillouin_zone.hpp

@@ -1,99 +0,0 @@
-/*******************************************************************************
- *
- * TRIQS: a Toolbox for Research in Interacting Quantum Systems
- *
- * Copyright (C) 2011 by M. Ferrero, O. Parcollet
- *
- * TRIQS is free software: you can redistribute it and/or modify it under the
- * terms of the GNU General Public License as published by the Free Software
- * Foundation, either version 3 of the License, or (at your option) any later
- * version.
- *
- * TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY
- * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
- * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
- * details.
- *
- * You should have received a copy of the GNU General Public License along with
- * TRIQS. If not, see <http://www.gnu.org/licenses/>.
- *
- ******************************************************************************/
-#ifndef TRIQS_LATTICE_BRAVAIS_LATTICE_H
-#define TRIQS_LATTICE_BRAVAIS_LATTICE_H
-#include <triqs/arrays/array.hpp>
-#include <triqs/arrays/matrix.hpp>
-#include <triqs/arrays/vector.hpp>
-#include <triqs/utility/exceptions.hpp>
-#include <unordered_map>
-#include <map>
-#include <string>
-#include <unordered_map>
-
-namespace triqs { namespace lattice_tools { 
-
- namespace tqa = triqs::arrays; 
- using tqa::array; using tqa::array_view; using tqa::matrix_view;using tqa::matrix;using tqa::range;
-
- typedef tqa::vector <double> R_type;
- typedef tqa::vector <double> K_type;
- typedef tqa::vector_view<double> R_view_type;
- typedef tqa::vector_view<double> K_view_type;
- typedef std::complex<double> dcomplex;
-
- class bravais_lattice { 
-  public : 
-   typedef matrix<double> units_type;
-   typedef std::unordered_map<std::string, R_type> orbital_type; // name -> position in space
-  
-   bravais_lattice(units_type const & units); 
-   
-   size_t n_orbitals() const {return orbitals_.size();}
-   units_type const & units() const { return units_;}
-   size_t dim() const { return dim_; }   
-
-  /**
-    * Push_back mechanism of a pair orbital_name -> position in R^3
-    * VectorType is anything from which a tqa::vector <double> can be constructed
-    */
-   template<typename VectorType>
-    void push_back(std::string const & s, VectorType && v) { orbitals_.insert(std::make_pair(s,std::forward<VectorType>(v)));}
-    
-   /***
-    * Transform into real coordinates.
-    * @param[in] x : coordinates in the basis :unit_i arrays::vector or vector_view
-    * @return  Coordinates in R^3 ("real" coordinates) 
-    */
-   template<typename R> 
-    R_view_type lattice_to_real_coordinates(R const & x) const { 
-     R_type res(3); res() =0;
-     for (size_t i =0; i< dim();i++) res += x (i) * units_(i,range());
-     return(res);
-    }
-
-  protected :
-   units_type units_;
-   size_t dim_;
-   orbital_type orbitals_;
- };
-
- class brillouin_zone { 
-  bravais_lattice lattice_;
-  tqa::matrix<double> K_reciprocal, K_reciprocal_inv;
-  public : 
-
-  brillouin_zone( bravais_lattice const & bl_);
-  bravais_lattice lattice() const { return lattice_;}
-
-  /***
-   * Transform into real coordinates.
-   * @param[in] k : coordinates in the basis (K_reciprocal_i)
-   * @return  Coordinates in R^3 ("real" coordinates) 
-   */
-  K_view_type lattice_to_real_coordinates (K_view_type const & k) const; 
-
-  /// Inverse of latt_to_real_k
-  K_view_type real_to_lattice_coordinates (K_view_type const & k) const;
- };
-
-}}
-#endif

triqs/lattice/brillouin_zone.hpp

@@ -0,0 +1,75 @@
+/*******************************************************************************
+ *
+ * TRIQS: a Toolbox for Research in Interacting Quantum Systems
+ *
+ * Copyright (C) 2011-2014 by M. Ferrero, O. Parcollet
+ *
+ * TRIQS is free software: you can redistribute it and/or modify it under the
+ * terms of the GNU General Public License as published by the Free Software
+ * Foundation, either version 3 of the License, or (at your option) any later
+ * version.
+ *
+ * TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY
+ * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+ * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+ * details.
+ *
+ * You should have received a copy of the GNU General Public License along with
+ * TRIQS. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#pragma once
+#include "./bravais_lattice.hpp"
+
+namespace triqs {
+namespace lattice {
+
+ class brillouin_zone {
+
+  public:
+  using point_t = k_t; // domain concept
+
+  brillouin_zone() { // default construction, 3d cubic lattice
+   K_reciprocal = arrays::make_unit_matrix<double>(3);
+   K_reciprocal_inv() = K_reciprocal;
+  }
+
+  /// Construct from a bravais_lattice
+  brillouin_zone(bravais_lattice const& bl_);
+
+  /// Access to the underlying bravais lattice
+  bravais_lattice lattice() const { return lattice_; }
+
+  /// Transform from lattice to real coordinates
+  template <typename K> k_t lattice_to_real_coordinates(K const& k) const { return _transfo_impl(k, K_reciprocal); }
+
+  /// Transform from real to lattice coordinates
+  template <typename K> k_t real_to_lattice_coordinates(K const& k) const { return _transfo_impl(k, K_reciprocal_inv); }
+
+  /// Write into HDF5
+  friend void h5_write(h5::group fg, std::string subgroup_name, brillouin_zone const& bz);
+
+  /// Read from HDF5
+  friend void h5_read(h5::group fg, std::string subgroup_name, brillouin_zone& bz);
+
+  //  BOOST Serialization
+  friend class boost::serialization::access;
+  template <class Archive> void serialize(Archive& ar, const unsigned int version) {
+   ar& boost::serialization::make_nvp("bravais_lattice", lattice_);
+  }
+
+  private:
+  bravais_lattice lattice_;
+  arrays::matrix<double> K_reciprocal, K_reciprocal_inv;
+
+  template <typename K> k_t _transfo_impl(K const& k, arrays::matrix<double> const& K_base) const {
+   if (first_dim(k) != lattice().dim()) TRIQS_RUNTIME_ERROR << "latt_to_real_k : dimension of k must be " << lattice().dim();
+   k_t res(3);
+   res() = 0;
+   int dim = lattice().dim();
+   for (int i = 0; i < dim; i++) res += k(i) * K_reciprocal_inv(i, range{});
+   return res;
+  }
+ };
+}
+}

triqs/lattice/bz_mesh.cpp

@@ -0,0 +1,75 @@
+/*******************************************************************************
+ *
+ * TRIQS: a Toolbox for Research in Interacting Quantum Systems
+ *
+ * Copyright (C) 2011 by M. Ferrero, O. Parcollet
+ *
+ * TRIQS is free software: you can redistribute it and/or modify it under the
+ * terms of the GNU General Public License as published by the Free Software
+ * Foundation, either version 3 of the License, or (at your option) any later
+ * version.
+ *
+ * TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY
+ * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+ * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+ * details.
+ *
+ * You should have received a copy of the GNU General Public License along with
+ * TRIQS. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#include "./bz_mesh.hpp"
+#include "./grid_generator.hpp"
+namespace triqs {
+namespace lattice {
+
+ bz_mesh::bz_mesh(brillouin_zone const &bz, int n_l) : bz(bz) {
+  // compute the k points
+  for (auto grid = grid_generator{bz.lattice().dim(), n_l}; grid; ++grid) k_pt_stack.push_back(*grid);
+ }
+
+ // locate the closest point : VERY PRIMITIVE : TO BE IMPROVED (kd-tree ?)
+ long bz_mesh::locate_neighbours(k_t k) const {
+  // reduction to the BZ ?
+  // very slow : improve with a kd-tree ?
+  double d_min = 100; // > pi
+  auto sqr = [](double x) { return x * x; };
+  int pos = 0, pos_min = 0;
+  for (auto const &p : k_pt_stack) {
+   int dim = k.size();
+   double d = 0;
+   for (int i = 0; i < dim; ++i) d += sqr(k(i) - p(i));
+   if (d < d_min) {
+    d_min = d;
+    pos_min = pos;
+   }
+   pos++;
+  }
+  return pos;
+ }
+
+ /// Write into HDF5
+ void h5_write(h5::group fg, std::string subgroup_name, bz_mesh const &m) {
+  h5::group gr = fg.create_group(subgroup_name);
+  h5_write(gr, "domain", m.domain());
+  auto kp = matrix<double>(m.k_pt_stack.size(), 3);
+  kp() = 0;
+  for (int i = 0; i < m.k_pt_stack.size(); ++i) kp(i, range{}) = m.k_pt_stack[i];
+  h5_write(gr, "k_pt_stack", kp);
+ }
+
+ /// Read from HDF5
+ void h5_read(h5::group fg, std::string subgroup_name, bz_mesh &m) {
+  h5::group gr = fg.open_group(subgroup_name);
+  brillouin_zone dom;
+  h5_read(gr, "domain", dom);
+  matrix<double> kp;
+  h5_read(gr, "k_pt_stack", kp);
+  std::vector<k_t> k_pt_stack;
+  int s = first_dim(kp);
+  for (int i = 0; i < s; ++i) k_pt_stack[i] = kp(i, range{});
+  m = bz_mesh(std::move(dom), k_pt_stack);
+ }
+}
+} // namespaces
+

triqs/lattice/bz_mesh.hpp

@@ -0,0 +1,109 @@
+/*******************************************************************************
+ *
+ * TRIQS: a Toolbox for Research in Interacting Quantum Systems
+ *
+ * Copyright (C) 2012 by M. Ferrero, O. Parcollet
+ *
+ * TRIQS is free software: you can redistribute it and/or modify it under the
+ * terms of the GNU General Public License as published by the Free Software
+ * Foundation, either version 3 of the License, or (at your option) any later
+ * version.
+ *
+ * TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY
+ * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+ * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+ * details.
+ *
+ * You should have received a copy of the GNU General Public License along with
+ * TRIQS. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#pragma once
+#include "./brillouin_zone.hpp"
+#include "../gfs/tools.hpp"
+#include "../gfs/meshes/mesh_tools.hpp"
+
+namespace triqs {
+namespace lattice {
+
+ struct bz_mesh {
+
+  using domain_t = brillouin_zone;
+  using index_t = long;
+  using domain_pt_t = typename domain_t::point_t;
+
+  bz_mesh() = default;
+  bz_mesh(brillouin_zone const &bz, int n_l);
+  bz_mesh(brillouin_zone const &bz, std::vector<k_t> k_pt_stack) : bz(bz), k_pt_stack(std::move(k_pt_stack)) {}
+
+  domain_t const &domain() const { return bz; }
+  long size() const { return k_pt_stack.size(); }
+
+  /// Conversions point <-> index <-> linear_index
+  domain_pt_t const &index_to_point(index_t i) const { return k_pt_stack[i]; }
+
+  long index_to_linear(index_t ind) const { return ind; }
+
+  // f (k) -> void where k is a k_t, a point in the BZ
+  template <typename F> friend void foreach(bz_mesh const &m, F f) {
+   for (auto const &k : m.k_pt_stack) f(k);
+  }
+
+  // locate the closest point : VERY PRIMITIVE : TO BE IMPROVED (kd-tree ?) 
+  long locate_neighbours(k_t k) const;
+
+  /// The wrapper for the mesh point
+  class mesh_point_t : gfs::tag::mesh_point, public utility::arithmetic_ops_by_cast<mesh_point_t, domain_pt_t> {
+   bz_mesh const *m;
+   index_t _index;
+
+   public:
+   mesh_point_t() : m(nullptr) {}
+   mesh_point_t(bz_mesh const &mesh, index_t const &index_) : m(&mesh), _index(index_) {}
+   mesh_point_t(bz_mesh const &mesh) : mesh_point_t(mesh, 0) {}
+   void advance() { ++_index; }
+   using cast_t = domain_pt_t;
+   operator cast_t() const { return m->index_to_point(_index); }
+   long linear_index() const { return _index; }
+   long index() const { return _index; }
+   bool at_end() const { return (_index == m->size()); }
+   void reset() { _index = 0; }
+  };
+
+  /// Accessing a point of the mesh
+  mesh_point_t operator[](index_t i) const { return mesh_point_t(*this, i); }
+
+  /// Iterating on all the points...
+  using const_iterator = gfs::mesh_pt_generator<bz_mesh>;
+  const_iterator begin() const { return const_iterator(this); }
+  const_iterator end() const { return const_iterator(this, true); }
+  const_iterator cbegin() const { return const_iterator(this); }
+  const_iterator cend() const { return const_iterator(this, true); }
+
+  /// Mesh comparison
+  bool operator==(bz_mesh const &M) const { return ( (size() == M.size())); }
+  //bool operator==(bz_mesh const &M) const { return ((bz == M.bz) && (size() == M.size())); }
+  bool operator!=(bz_mesh const &M) const { return !(operator==(M)); }
+
+  /// Write into HDF5
+  friend void h5_write(h5::group fg, std::string subgroup_name, bz_mesh const &m);
+
+  /// Read from HDF5
+  friend void h5_read(h5::group fg, std::string subgroup_name, bz_mesh &m);
+
+  //  BOOST Serialization
+  friend class boost::serialization::access;
+  template <class Archive> void serialize(Archive &ar, const unsigned int version) {
+   ar &boost::serialization::make_nvp("domain", bz);
+   ar &boost::serialization::make_nvp("k_pt_stack", k_pt_stack);
+  }
+
+  friend std::ostream &operator<<(std::ostream &sout, bz_mesh const &m) { return sout << "Mesh over BZ "; }
+
+  private:
+  domain_t bz;
+  std::vector<k_t> k_pt_stack;
+ };
+}
+}
+

triqs/lattice/concepts.rst

@@ -1,28 +0,0 @@
-
-
-Concept Function 
-return_type
-arg_type
-return_type operator()(arg_type const &) const  : the evaluation of the function
-
-
-Concept FunctionOnBravaisLattice
-
-refines Function
-bravais_lattice const & lattice() const 
-
-
-Concept ShortRangeFunctionOnBravaisLattice
-
-* refines FunctionOnBravaisLattice
-* derived from Tag::ShortRangeFunctionOnBravaisLattice  ??
-* map_type
-* begin, end: iterator -> pair<arg_type, return_type>
-
-
-Concept FunctionOnBrillouinZone
-
-* refines Function
-brillouin_zone const & bz() const
-
-

triqs/lattice/functors.hpp

@@ -1,90 +0,0 @@
-
-/*******************************************************************************
- *
- * TRIQS: a Toolbox for Research in Interacting Quantum Systems
- *
- * Copyright (C) 2011 by M. Ferrero, O. Parcollet
- *
- * TRIQS is free software: you can redistribute it and/or modify it under the
- * terms of the GNU General Public License as published by the Free Software
- * Foundation, either version 3 of the License, or (at your option) any later
- * version.
- *
- * TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY
- * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
- * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
- * details.
- *
- * You should have received a copy of the GNU General Public License along with
- * TRIQS. If not, see <http://www.gnu.org/licenses/>.
- *
- ******************************************************************************/
-#ifndef TRIQS_LATTICE_FUNCTORS_H
-#define TRIQS_LATTICE_FUNCTORS_H
-
-namespace triqs { namespace lattice_tools { 
-
- const double pi = acos(-1.0);
- const std::complex<double> I(0,1);
-
- template<typename F> struct minus_chech_impl {
-  typedef typename F::arg_type arg_type;
-  typedef typename F::return_type return_type;
-  F f;
-  minus_chech_impl(F const & f_):f(f_){}
-  brillouin_zone const & bz() const {return f.bz();}
-  return_type operator()(arg_type const & x) const { return_type res(f(x)); res *=-1; return res;} 
- };
-
- namespace result_of { template<typename F> struct minus_chech{ typedef minus_chech_impl<F> type;}; }
-
- /**
-  *  Given f of type F which models FunctionOnBravaisLattice, minus_check(f) : 
-  *   - returns -f(-args)
-  *   - its type models Function
-  *
-  */
- template<typename F> minus_chech_impl<F> minus_chech(F const & f) { return minus_chech_impl<F> (f);}
-
- template<typename F >//, typename Enabler = boost::enable_if< Tag::check<Tag::ShortRangeFunctionOnBravaisLattice, F> > >
-  class fourier_impl {
-   F f;
-   brillouin_zone bz_;
-  // deduce the return type from decltype(begin()->second)
-   public: 
-   typedef typename regular_type_if_exists_else_type< decltype(f.begin()->second)>::type return_construct_type;
-   typedef typename view_type_if_exists_else_type<return_construct_type>::type return_type;
-   typedef K_view_type arg_type;
-
-   fourier_impl (F f_):f(f_), bz_(f_.lattice()), res(f.n_bands(),f.n_bands()) {}
-
-   //brillouin_zone const & bz() const { return bz_; }
-
-   return_type operator()(K_view_type const & k) {
-    res()=0;
-    for (auto const & pdm : f) { res += pdm.second * exp( 2*pi*I* this->dot_product(k,pdm.first)); }
-    return res;
-   }
-
-   protected:
-   inline double dot_product(K_view_type const & a, typename F::arg_type const & b) const {
-    assert(b.size()>= this->bz_.lattice().dim()); 
-    double r=0; for (size_t i=0; i< this->bz_.lattice().dim();++i) r += a(i) * b[i];
-    return r;
-   }
-   return_construct_type res;
-   //typename F::return_construct_type res;
-  };
-
- /**
-  *  Given f of type F which models ShortRangeFunctionOnBravaisLattice, Fourier(f) returns 
-  *   - a type that models FunctionOnBravaisLattice
-  *   - and returns the Fourier transform f(k)
-  */
- template<typename F> fourier_impl<F> Fourier(F f) { return fourier_impl<F> (f);}
- //namespace result_of { template<typename F> struct Fourier { typedef fourier_impl<F> type;}; }
-
-}}
-
-#endif
-

triqs/lattice/grid_generator.hpp

@@ -1,4 +1,3 @@
-
 /*******************************************************************************
  *
  * TRIQS: a Toolbox for Research in Interacting Quantum Systems
@@ -19,55 +18,77 @@
  * TRIQS. If not, see <http://www.gnu.org/licenses/>.
  *
  ******************************************************************************/
+#pragma once
+#include "./brillouin_zone.hpp"
 
-#ifndef TRIQS_GRID_GENERATOR_H
+namespace triqs {
-#define TRIQS_GRID_GENERATOR_H
+namespace lattice {
-
-#include <triqs/arrays/array.hpp>
-
-namespace triqs { namespace lattice_tools { 
 
  /**
   * Generate the point in a cuboid as an array<double,1> const &
   */
- class grid_generator : 
+ class grid_generator : public boost::iterator_facade<grid_generator, k_t const &, boost::forward_traversal_tag, k_t const &> {
-  public boost::iterator_facade< grid_generator, K_type const &, boost::forward_traversal_tag, K_type const & > {
   friend class boost::iterator_core_access;
-   size_t dim, nkpts,nx,ny,nz,N_X,N_Y,N_Z,index_;
+  int dim, nkpts, nx, ny, nz, N_X, N_Y, N_Z, index_;
   double step_x, step_y, step_z;
   bool at_end;
-   K_type pt;
+  k_t pt;
   void init() {
-    nx=0; ny=0; nz=0;at_end = false;index_ =0;
+   nx = 0;
+   ny = 0;
+   nz = 0;
+   at_end = false;
+   index_ = 0;
    N_X = nkpts;
    N_Y = (dim > 1 ? nkpts : 1);
    N_Z = (dim > 2 ? nkpts : 1);
-    step_x =  1.0/double(N_X); step_y =  1.0/double(N_Y); step_z =  1.0/double(N_Z);
+   step_x = 1.0 / double(N_X);
-    pt(0)=step_x/2; pt(1)=step_y/2; pt(2)=step_z/2;
+   step_y = 1.0 / double(N_Y);
+   step_z = 1.0 / double(N_Z);
+   pt(0) = step_x / 2;
+   pt(1) = step_y / 2;
+   pt(2) = step_z / 2;
   }
 
   void increment() {
-    if (nx<N_X-1) { ++nx; pt(0) += step_x; ++index_; return;}
+   if (nx < N_X - 1) {
-    pt(0) = step_x/2; nx = 0;
+    ++nx;
-    if (ny<N_Y-1) { ++ny; pt(1) += step_y; ++index_; return;}
+    pt(0) += step_x;
-    pt(1) = step_y/2; ny = 0;
+    ++index_;
-    if (nz<N_Z-1) { ++nz; pt(2) += step_z; ++index_; return;}
+    return;
+   }
+   pt(0) = step_x / 2;
+   nx = 0;
+   if (ny < N_Y - 1) {
+    ++ny;
+    pt(1) += step_y;
+    ++index_;
+    return;
+   }
+   pt(1) = step_y / 2;
+   ny = 0;
+   if (nz < N_Z - 1) {
+    ++nz;
+    pt(2) += step_z;
+    ++index_;
+    return;
+   }
    at_end = true;
   }
 
   value_type dereference() const { return pt; }
-   bool equal(grid_generator const & other) const { return ((other.dim == dim) && (other.index_==index_) && (other.nkpts==nkpts));}
+  bool equal(grid_generator const &other) const {
+   return ((other.dim == dim) && (other.index_ == index_) && (other.nkpts == nkpts));
+  }
 
   public:
   /// dim : dimension, nkpts : number of k point in each dimension
-   grid_generator(size_t dim_, size_t nkpts_): dim(dim_),nkpts(nkpts_), pt(3) {init();}
+  grid_generator(int dim_, int nkpts_) : dim(dim_), nkpts(nkpts_), pt(3) { init(); }
   grid_generator() : dim(3), nkpts(0), pt(3) { init(); }
-   size_t size () const { return (N_X * N_Y * N_Z);}
+  int size() const { return (N_X * N_Y * N_Z); }
-   size_t index() const { return index_;}
+  int index() const { return index_; }
   operator bool() const { return !(at_end); }
  };
-
+}
-}}
+}
-#endif
-
 

triqs/lattice/tight_binding.cpp

@@ -22,59 +22,72 @@
 #include <triqs/arrays/algorithms.hpp>
 #include <triqs/arrays/linalg/eigenelements.hpp>
 #include "grid_generator.hpp"
-#include "functors.hpp"
+namespace triqs {
-namespace triqs { namespace lattice_tools { 
+namespace lattice {
-using namespace std;
+
- using namespace tqa;
+ using namespace arrays;
+
+ tight_binding::tight_binding(bravais_lattice const& bl, std::vector<std::vector<long>> all_disp,
+                              std::vector<matrix<dcomplex>> all_matrices)
+    : bl_(bl), all_disp(std::move(all_disp)), all_matrices(std::move(all_matrices)) {
+
+  // checking inputs
+  if (all_disp.size() != all_matrices.size()) TRIQS_RUNTIME_ERROR << " Number of displacements != Number of matrices";
+  for (int i = 0; i < all_disp.size(); ++i) {
+   if (all_disp[i].size() != bl_.dim())
+    TRIQS_RUNTIME_ERROR << "displacement of incorrect size : got " << all_disp[i].size() << "instead of " << bl_.dim();
+   if (first_dim(all_matrices[i]) != n_bands())
+    TRIQS_RUNTIME_ERROR << "the first dim matrix is of size " << first_dim(all_matrices[i]) << " instead of " << n_bands();
+   if (second_dim(all_matrices[i]) != n_bands())
+    TRIQS_RUNTIME_ERROR << "the second dim matrix is of size " << second_dim(all_matrices[i]) << " instead of " << n_bands();
+  }
+ }
 
  //------------------------------------------------------
-
+ array<dcomplex, 3> hopping_stack(tight_binding const& TB, arrays::array_const_view<double, 2> k_stack) {
- array_view <dcomplex,3> hopping_stack (tight_binding const & TB, array_view<double,2> const & k_stack) {
+  auto TK = fourier(TB);
-  auto TK = Fourier(TB);  
   array<dcomplex, 3> res(TB.n_bands(), TB.n_bands(), k_stack.shape(1));
-  for(size_t i = 0; i<k_stack.shape(1); ++i) res(range(), range(), i) = TK(k_stack(range(),i));
+  for (int i = 0; i < k_stack.shape(1); ++i) res(range(), range(), i) = TK(k_stack(range(), i));
   return res;
  }
 
  //------------------------------------------------------
-
+ array<double, 2> energies_on_bz_path(tight_binding const& TB, k_t const& K1, k_t const& K2, int n_pts) {
- array_view<double,2> energies_on_bz_path(tight_binding const & TB, K_view_type K1, K_view_type K2, size_t n_pts) {
+  auto TK = fourier(TB);
-  auto TK = Fourier(TB);
+  int norb = TB.lattice().n_orbitals();
-  const size_t norb=TB.lattice().n_orbitals();
+  int ndim = TB.lattice().dim();
-  const size_t ndim=TB.lattice().dim();
   array<double, 2> eval(norb, n_pts);
-  K_type dk = (K2 - K1)/double(n_pts), k = K1;
+  k_t dk = (K2 - K1) / double(n_pts), k = K1;
-  for (size_t i =0; i<n_pts; ++i, k += dk) {
+  for (int i = 0; i < n_pts; ++i, k += dk) {
-   eval(range(),i) = linalg::eigenvalues( TK( k (range(0,ndim))), false);
+   eval(range(), i) = linalg::eigenvalues(TK(k(range(0, ndim)))(), false);
   }
   return eval;
  }
 
  //------------------------------------------------------
- array_view<double,2> energies_on_bz_grid(tight_binding const & TB, size_t n_pts) {
+ array<double, 2> energies_on_bz_grid(tight_binding const& TB, int n_pts) {
 
-  auto TK = Fourier(TB);
+  auto TK = fourier(TB);
-  const size_t norb=TB.lattice().n_orbitals();
+  int norb = TB.lattice().n_orbitals();
-  const size_t ndim=TB.lattice().dim();
+  int ndim = TB.lattice().dim();
   grid_generator grid(ndim, n_pts);
   array<double, 2> eval(norb, grid.size());
   for (; grid; ++grid) {
-   eval(range(),grid.index()) = linalg::eigenvalues( TK( (*grid) (range(0,ndim))), false);
+   eval(range(), grid.index()) = linalg::eigenvalues(TK((*grid)(range(0, ndim)))(), false);
   }
   return eval;
  }
 
  //------------------------------------------------------
 
- std::pair<array<double,1>, array<double,2> > dos(tight_binding const & TB, size_t nkpts, size_t neps) {
+ std::pair<array<double, 1>, array<double, 2>> dos(tight_binding const& TB, int nkpts, int neps) {
 
-  // The Fourier transform of TK
+  // The fourier transform of TK
-  // auto TK = Fourier(TB); // C++0x .... 
+  auto TK = fourier(TB);
-  auto TK = Fourier(TB);
 
   // loop on the BZ
-  const size_t ndim=TB.lattice().dim();
+  int ndim = TB.lattice().dim();
-  const size_t norb=TB.lattice().n_orbitals();
+  int norb = TB.lattice().n_orbitals();
   grid_generator grid(ndim, nkpts);
   array<double, 1> tempeval(norb);
   array<dcomplex, 3> evec(norb, norb, grid.size());
@@ -96,19 +109,20 @@ using namespace std;
 
   // define the epsilon mesh, etc.
   array<double, 1> epsilon(neps);
-  double epsmax = tqa::max_element(eval);
+  double epsmax = max_element(eval);
-  double epsmin = tqa::min_element(eval);
+  double epsmin = min_element(eval);
   double deps = (epsmax - epsmin) / neps;
-  //for (size_t i =0; i< neps; ++i) epsilon(i)= epsmin+i/(neps-1.0)*(epsmax-epsmin);
+  // for (int i =0; i< neps; ++i) epsilon(i)= epsmin+i/(neps-1.0)*(epsmax-epsmin);
-  for (size_t i =0; i< neps; ++i) epsilon(i)=epsmin+(i+0.5)*deps;
+  for (int i = 0; i < neps; ++i) epsilon(i) = epsmin + (i + 0.5) * deps;
 
   // bin the eigenvalues according to their energy
   // NOTE: a is defined as an integer. it is the index for the DOS.
   // REPORT <<"Starting Binning ...."<<endl;
-  array<double,2> rho (neps,norb);rho()=0;
+  array<double, 2> rho(neps, norb);
-  for(size_t l=0;l<norb;l++){
+  rho() = 0;
-   for (size_t j=0;j<grid.size();j++){
+  for (int l = 0; l < norb; l++) {
-    for (size_t k=0;k<norb;k++){
+   for (int j = 0; j < grid.size(); j++) {
+    for (int k = 0; k < norb; k++) {
      int a = int((eval(k, j) - epsmin) / deps);
      if (a == int(neps)) a = a - 1;
      rho(a, l) += real(conj(evec(l, k, j)) * evec(l, k, j));
@@ -123,26 +137,28 @@ using namespace std;
 
  //----------------------------------------------------------------------------------
 
- std::pair<array<double,1>, array<double,1> > dos_patch(tight_binding const & TB, const array<double,2> & triangles, size_t neps, size_t ndiv) {
+ std::pair<array<double, 1>, array<double, 1>> dos_patch(tight_binding const& TB, const array<double, 2>& triangles, int neps,
+                                                         int ndiv) {
   // WARNING: This version only works for a single band Hamiltonian in 2 dimensions!!!!
   // triangles is an array of points defining the triangles of the patch
   // neps in the number of bins in energy
   // ndiv in the number of divisions used to divide the triangles
 
-  //const size_t ndim=TB.lattice().dim();
+  // int ndim=TB.lattice().dim();
-  //const size_t norb=TB.lattice().n_orbitals();
+  // int norb=TB.lattice().n_orbitals();
   int ntri = triangles.shape(0) / 3;
   array<double, 1> dos(neps);
 
   // Check consistency
-  const size_t ndim=TB.lattice().dim();
+  int ndim = TB.lattice().dim();
-  //const size_t norb=TB.lattice().n_orbitals();
+  // int norb=TB.lattice().n_orbitals();
   if (ndim != 2) TRIQS_RUNTIME_ERROR << "dos_patch : dimension 2 only !";
-  if (triangles.shape(1) != ndim) TRIQS_RUNTIME_ERROR<<"dos_patch : the second dimension of the 'triangle' array in not "<<ndim;
+  if (triangles.shape(1) != ndim)
+   TRIQS_RUNTIME_ERROR << "dos_patch : the second dimension of the 'triangle' array in not " << ndim;
 
   // Every triangle has ndiv*ndiv k points
-  size_t nk = ntri*ndiv*ndiv;
+  int nk = ntri * ndiv * ndiv;
-  size_t k_index = 0;
+  int k_index = 0;
   double epsmax = -100000, epsmin = 100000;
   array<dcomplex, 2> thop(1, 1);
   array<double, 1> energ(nk), weight(nk);
@@ -153,8 +169,8 @@ using namespace std;
   int pt = 0;
   double s, t;
 
-  // The Fourier transform of TK
+  // The fourier transform of TK
-  auto TK = Fourier(TB); 
+  auto TK = fourier(TB);
 
   // loop over the triangles
   for (int tri = 0; tri < ntri; tri++) {
@@ -170,11 +186,11 @@ using namespace std;
    // so I use it to weight the sum in the dos
    double area = abs(0.5 * ((b(0) - a(0)) * (c(1) - a(1)) - (b(1) - a(1)) * (c(0) - a(0))) / (ndiv * ndiv));
 
-   for (size_t i = 0; i<ndiv; i++) {
+   for (int i = 0; i < ndiv; i++) {
     s = i / double(ndiv);
-    for (size_t j = 0; j<ndiv-i; j++) {
+    for (int j = 0; j < ndiv - i; j++) {
      t = j / double(ndiv);
-     for (size_t k = 0; k<2; k++) {
+     for (int k = 0; k < 2; k++) {
 
       rv = a + s * (b - a) + t * (c - a) + (k + 1.0) * g;
 
@@ -198,25 +214,21 @@ using namespace std;
   // define the epsilon mesh, etc.
   array<double, 1> epsilon(neps);
   double deps = (epsmax - epsmin) / neps;
-  for (size_t i =0; i< neps; ++i) epsilon(i)= epsmin+i/(neps-1.0)*(epsmax-epsmin);
+  for (int i = 0; i < neps; ++i) epsilon(i) = epsmin + i / (neps - 1.0) * (epsmax - epsmin);
 
   // bin the eigenvalues according to their energy
   int ind;
   double totalweight(0.0);
   dos() = 0.0;
-  for (size_t j = 0; j < nk; j++) {
+  for (int j = 0; j < nk; j++) {
    ind = int((energ(j) - epsmin) / deps);
    if (ind == int(neps)) ind--;
    dos(ind) += weight(j);
    totalweight += weight(j);
   }
   dos /= deps; // Normalize the DOS
-  return std::make_pair(epsilon, dos);
+  return {std::move(epsilon), std::move(dos)};
+ }
+}
 }
 
-
-
-
-
- }}
-

triqs/lattice/tight_binding.hpp

@@ -18,74 +18,82 @@
  * TRIQS. If not, see <http://www.gnu.org/licenses/>.
  *
  ******************************************************************************/
-#ifndef TRIQS_LATTICE_TIGHTBINDINGS_H
+#pragma once
-#define TRIQS_LATTICE_TIGHTBINDINGS_H
+#include "brillouin_zone.hpp"
-#include "bravais_lattice_and_brillouin_zone.hpp"
+#include <triqs/utility/complex_ops.hpp>
 
-namespace triqs { namespace lattice_tools { 
+namespace triqs {
+namespace lattice {
 
  /**
    For tightbinding Hamiltonian with fully localised orbitals
-   Model the ShortRangeFunctionOnBravaisLattice concept.
    Overlap between orbital is taken as unit matrix.
   */
  class tight_binding {
-   typedef std::vector<std::pair<std::vector<long>, matrix<dcomplex>>> displ_value_stack_t;
-   displ_value_stack_t displ_value_stack;
-   bravais_lattice bl_;
-  public : 
-   typedef std::vector<long>                                      arg_type;
 
+  bravais_lattice bl_;
+  std::vector<std::vector<long>> all_disp;
+  std::vector<matrix<dcomplex>> all_matrices;
+
+  public:
   ///
-   tight_binding (bravais_lattice const & bl) : bl_(bl) {};
+  tight_binding(bravais_lattice const& bl, std::vector<std::vector<long>> all_disp, std::vector<matrix<dcomplex>> all_matrices);
 
   /// Underlying lattice
   bravais_lattice const& lattice() const { return bl_; }
 
   /// Number of bands, i.e. size of the matrix t(k)
-   size_t n_bands() const { return bl_.n_orbitals();}
+  int n_bands() const { return bl_.n_orbitals(); }
 
-   /**
+  // calls F(R, t(R)) for all R
-    * Push_back mechanism of a pair displacement -> matrix
+  template <typename F> friend void foreach(tight_binding const& tb, F f) {
-    * VectorIntType is anything from which a std::vector<long> can be constructed
+   int n = tb.all_disp.size();
-    * MatrixDComplexType is anything from which a matrix<dcomplex> can be constructed
+   for (int i = 0; i < n; ++i) f(tb.all_disp[i], tb.all_matrices[i]);
-    */
-   template<typename VectorIntType, typename MatrixDComplexType>
-    void push_back(VectorIntType && v, MatrixDComplexType && m) {
-     std::vector<long> V(std::forward<VectorIntType>(v)); 
-     if (v.size() != bl_.dim()) TRIQS_RUNTIME_ERROR<<"tight_binding : displacement of incorrect size : got "<< v.size() << "instead of "<< bl_.dim();
-     matrix<dcomplex> M(std::forward<MatrixDComplexType>(m));
-     if (first_dim(M) != n_bands()) TRIQS_RUNTIME_ERROR<<"tight_binding : the first dim matrix is of size "<<  first_dim(M) <<" instead of "<< n_bands();
-     if (second_dim(M) != n_bands()) TRIQS_RUNTIME_ERROR<<"tight_binding : the first dim matrix is of size "<<  second_dim(M) <<" instead of "<< n_bands();
-     displ_value_stack.push_back(std::make_pair(std::move(V), std::move(M)));
   }
 
-   void reserve(size_t n) { displ_value_stack.reserve(n);}
+  // a simple function on the domain brillouin_zone
+  struct fourier_impl {
+   tight_binding const& tb;
+   const int nb;
 
-   // iterators
+   using domain_t = brillouin_zone;
-   typedef displ_value_stack_t::const_iterator const_iterator; 
+   brillouin_zone domain() const {
-   typedef displ_value_stack_t::iterator iterator; 
+    return {tb.lattice()};
-   const_iterator begin() const { return displ_value_stack.begin();} 
+   }
-   const_iterator end() const   { return displ_value_stack.end();} 
+
-   iterator begin() { return displ_value_stack.begin();} 
+   template <typename K> matrix<dcomplex> operator()(K const& k) {
-   iterator end()   { return displ_value_stack.end();} 
+    matrix<dcomplex> res(nb, nb);
+    res() = 0;
+    foreach(tb, [&](std::vector<long> const& displ, matrix<dcomplex> const& m) {
+     double dot_prod = 0;
+     int imax = displ.size();
+     for (int i = 0; i < imax; ++i) dot_prod += k(i) * displ[i];
+     //double dot_prod = k(0) * displ[0] + k(1) * displ[1];
+     res += m * exp(2_j * M_PI * dot_prod);
+    });
+    return res;
+   }
   };
 
+  fourier_impl friend fourier(tight_binding const& tb) {
+   return {tb, tb.n_bands()};
+  }
+
+ }; // tight_binding
+
  /**
    Factorized version of hopping (for speed)
    k_in[:,n] is the nth vector
    In the result, R[:,:,n] is the corresponding hopping t(k)
    */
- array_view <dcomplex,3> hopping_stack (tight_binding const & TB, array_view<double,2> const & k_stack);
+ array<dcomplex, 3> hopping_stack(tight_binding const& TB, arrays::array_const_view<double, 2> k_stack);
  // not optimal ordering here
 
- std::pair<array<double,1>, array<double,2> > dos(tight_binding const & TB, size_t nkpts, size_t neps); 
+ std::pair<array<double, 1>, array<double, 2>> dos(tight_binding const& TB, int nkpts, int neps);
- std::pair<array<double,1>, array<double,1> > dos_patch(tight_binding const & TB, const array<double,2> & triangles, size_t neps, size_t ndiv);
+ std::pair<array<double, 1>, array<double, 1>> dos_patch(tight_binding const& TB, const array<double, 2>& triangles, int neps,
- array_view<double,2> energies_on_bz_path(tight_binding const & TB, K_view_type K1, K_view_type K2, size_t n_pts);
+                                                         int ndiv);
- array_view<double,2> energies_on_bz_grid(tight_binding const & TB, size_t n_pts);
+ array<double, 2> energies_on_bz_path(tight_binding const& TB, k_t const& K1, k_t const& K2, int n_pts);
-
+ array<double, 2> energies_on_bz_grid(tight_binding const& TB, int n_pts);
-}}
+}
-
+}
-#endif
-