diff --git a/pytriqs/lattice/lattice_tools.pyx b/pytriqs/lattice/lattice_tools.pyx index 7998faf3..530cd298 100644 --- a/pytriqs/lattice/lattice_tools.pyx +++ b/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 cppclass bravais_lattice_c " triqs::lattice_tools::bravais_lattice" : - bravais_lattice_c(matrix[double] & units) +cdef extern from "triqs/lattice/brillouin_zone.hpp" : + cdef cppclass bravais_lattice_c " triqs::lattice::bravais_lattice" : + 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" : - tight_binding(bravais_lattice_c & units) - void push_back (vector[long] &, matrix[double] &) + cdef cppclass tight_binding " triqs::lattice::tight_binding" : + tight_binding(bravais_lattice_c & units, vector[vector[long]] &, vector[matrix[dcomplex]] &) except + #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)) - cdef vector[long] d + #self._c = new tight_binding(deref(bravais_lattice._c)) + #cdef vector[long] d + cdef vector[matrix[dcomplex]] mats + cdef vector[vector[long]] displs for displ, mat in hopping.items() : - d = displ - self._c.push_back( d, matrix[double] (numpy.array(mat, dtype =float))) + displs.push_back(displ) + #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 diff --git a/test/triqs/gfs/g_k_om.cpp b/test/triqs/gfs/g_k_om.cpp new file mode 100644 index 00000000..146ea6d0 --- /dev/null +++ b/test/triqs/gfs/g_k_om.cpp @@ -0,0 +1,36 @@ +#define TRIQS_ARRAYS_ENFORCE_BOUNDCHECK +#include +#include + +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(2)}}; + + auto g_eps = gf{{bz_, 1000}, {1, 1}}; + + auto G = gf>{{{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; +} diff --git a/triqs/gfs/bz.hpp b/triqs/gfs/bz.hpp new file mode 100644 index 00000000..bb397348 --- /dev/null +++ b/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 . + * + ******************************************************************************/ +#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 struct gf_mesh : lattice::bz_mesh { + template gf_mesh(T &&... x) : lattice::bz_mesh(std::forward(x)...) {} + }; + + namespace gfs_implementation { + + // h5 name + template struct h5_name { + static std::string invoke() { return "BZ"; } + }; + + /// --------------------------- data access --------------------------------- + template struct data_proxy : data_proxy_array, 3> {}; + template struct data_proxy : data_proxy_array, 1> {}; + + /// --------------------------- evaluator --------------------------------- + + // simple evaluation : take the point on the grid... + template <> struct evaluator_fnt_on_mesh { + lattice::k_t k; + evaluator_fnt_on_mesh() = default; + template evaluator_fnt_on_mesh(MeshType const &m, lattice::k_t const &k_) { k = k_; } + template auto operator()(F const &f) const DECL_AND_RETURN(f(k)); + }; + + // ------------- evaluator ------------------- + // handle the case where the matsu. freq is out of grid... + template struct evaluator { + static constexpr int arity = 1; + + template + std::c14::conditional_t::value, arrays::matrix_const_view>, + std::complex> + operator()(G const *g, lattice::k_t const &k) const { + auto n = g->mesh().locate_neighbours(k); // TO BE IMPROVED + return (*g)[n]; + } + }; + } +} +} + diff --git a/triqs/lattice/bravais_lattice.hpp b/triqs/lattice/bravais_lattice.hpp new file mode 100644 index 00000000..dbed5809 --- /dev/null +++ b/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 . + * + ******************************************************************************/ +#pragma once +#include +#include +#include + +namespace triqs { +namespace lattice { + + using r_t = arrays::vector; + using k_t = arrays::vector; + using dcomplex = std::complex; + using arrays::matrix; + using arrays::array; + using arrays::range; + + // -------------------------------------------------------- + + class bravais_lattice { + + public: + using point_t = std::vector; // domain concept. PUT on STACK + + bravais_lattice(matrix const& units, std::vector atom_orb_pos, std::vector atom_orb_name); + bravais_lattice(matrix const& units, std::vector atom_orb_pos) + : bravais_lattice(units, atom_orb_pos, std::vector(atom_orb_pos.size(), "")) {} + bravais_lattice(matrix const& units) : bravais_lattice(units, std::vector{{0, 0, 0}}) {} + bravais_lattice() : bravais_lattice(arrays::make_unit_matrix(3)) {} + + int n_orbitals() const { return atom_orb_name.size(); } + arrays::matrix_const_view units() const { return units_; } + int dim() const { return dim_; } + + /// Transform into real coordinates. + template 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 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 units_; + std::vector atom_orb_pos; // atom_orb_pos[i] = position of ith atoms/orbitals in the unit cell + std::vector atom_orb_name; // names of these atoms/orbitals. + int dim_; + }; +} +} diff --git a/triqs/lattice/bravais_lattice_and_brillouin_zone.cpp b/triqs/lattice/bravais_lattice_and_brillouin_zone.cpp index dfc00636..fb9a4942 100644 --- a/triqs/lattice/bravais_lattice_and_brillouin_zone.cpp +++ b/triqs/lattice/bravais_lattice_and_brillouin_zone.cpp @@ -18,79 +18,106 @@ * TRIQS. If not, see . * ******************************************************************************/ -#include "bravais_lattice_and_brillouin_zone.hpp" +#include "./bravais_lattice.hpp" +#include "./brillouin_zone.hpp" #include #include #include -namespace triqs { namespace lattice_tools { +namespace triqs { +namespace lattice { - using namespace tqa; - using namespace std; - //using triqs::arrays::blas::dot; - const double almost_zero(1E-10); + const double almost_zero = 1e-10; - bravais_lattice::bravais_lattice( units_type const & units__) : units_(3,3), dim_(units__.shape(0)) { - units_(range(0,dim_),range()) = units__(); - units_(range(dim_,3),range()) = 0; + bravais_lattice::bravais_lattice(matrix const& units__, std::vector atom_orb_pos_, + std::vector atom_orb_name_) + : 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 ux(3),uy(3),uz(3); - assert(dim_==2); + arrays::vector ux(3), uy(3), uz(3); + double delta; + auto _ = range{}; switch (dim_) { - case 1: - ux = units_(0,range()); - uz() = 0; uz(1) = 1 ; - uz = uz - dot(uz,ux)* ux; + case 1: + ux = units_(0, _); + 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)) 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()); - std::cout << Units << std::endl; - double delta = dot(Units(0,range()), cross_product(Units(1,range()),Units(2,range()))); - std::cout << dot(Units(0,range()), cross_product(Units(1,range()),Units(2,range())))<. - * - ******************************************************************************/ -#ifndef TRIQS_LATTICE_BRAVAIS_LATTICE_H -#define TRIQS_LATTICE_BRAVAIS_LATTICE_H -#include -#include -#include -#include -#include -#include -#include -#include - -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 R_type; - typedef tqa::vector K_type; - typedef tqa::vector_view R_view_type; - typedef tqa::vector_view K_view_type; - typedef std::complex dcomplex; - - class bravais_lattice { - public : - typedef matrix units_type; - typedef std::unordered_map 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 can be constructed - */ - template - void push_back(std::string const & s, VectorType && v) { orbitals_.insert(std::make_pair(s,std::forward(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 - 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 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 diff --git a/triqs/lattice/brillouin_zone.hpp b/triqs/lattice/brillouin_zone.hpp new file mode 100644 index 00000000..b9a4887d --- /dev/null +++ b/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 . + * + ******************************************************************************/ +#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(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 k_t lattice_to_real_coordinates(K const& k) const { return _transfo_impl(k, K_reciprocal); } + + /// Transform from real to lattice coordinates + template 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 void serialize(Archive& ar, const unsigned int version) { + ar& boost::serialization::make_nvp("bravais_lattice", lattice_); + } + + private: + bravais_lattice lattice_; + arrays::matrix K_reciprocal, K_reciprocal_inv; + + template k_t _transfo_impl(K const& k, arrays::matrix 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; + } + }; +} +} diff --git a/triqs/lattice/bz_mesh.cpp b/triqs/lattice/bz_mesh.cpp new file mode 100644 index 00000000..4996ef07 --- /dev/null +++ b/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 . + * + ******************************************************************************/ +#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(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 kp; + h5_read(gr, "k_pt_stack", kp); + std::vector 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 + diff --git a/triqs/lattice/bz_mesh.hpp b/triqs/lattice/bz_mesh.hpp new file mode 100644 index 00000000..3d870581 --- /dev/null +++ b/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 . + * + ******************************************************************************/ +#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_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 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 { + 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; + 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 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_pt_stack; + }; +} +} + diff --git a/triqs/lattice/concepts.rst b/triqs/lattice/concepts.rst deleted file mode 100644 index 6aa9ce58..00000000 --- a/triqs/lattice/concepts.rst +++ /dev/null @@ -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 - - -Concept FunctionOnBrillouinZone - -* refines Function -brillouin_zone const & bz() const - - diff --git a/triqs/lattice/functors.hpp b/triqs/lattice/functors.hpp deleted file mode 100644 index 53423770..00000000 --- a/triqs/lattice/functors.hpp +++ /dev/null @@ -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 . - * - ******************************************************************************/ -#ifndef TRIQS_LATTICE_FUNCTORS_H -#define TRIQS_LATTICE_FUNCTORS_H - -namespace triqs { namespace lattice_tools { - - const double pi = acos(-1.0); - const std::complex I(0,1); - - template 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 struct minus_chech{ typedef minus_chech_impl type;}; } - - /** - * Given f of type F which models FunctionOnBravaisLattice, minus_check(f) : - * - returns -f(-args) - * - its type models Function - * - */ - template minus_chech_impl minus_chech(F const & f) { return minus_chech_impl (f);} - - template//, typename Enabler = boost::enable_if< Tag::check > > - 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::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 fourier_impl Fourier(F f) { return fourier_impl (f);} - //namespace result_of { template struct Fourier { typedef fourier_impl type;}; } - -}} - -#endif - diff --git a/triqs/lattice/grid_generator.hpp b/triqs/lattice/grid_generator.hpp index f76e5138..85e20144 100644 --- a/triqs/lattice/grid_generator.hpp +++ b/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 . * ******************************************************************************/ +#pragma once +#include "./brillouin_zone.hpp" -#ifndef TRIQS_GRID_GENERATOR_H -#define TRIQS_GRID_GENERATOR_H +namespace triqs { +namespace lattice { -#include - -namespace triqs { namespace lattice_tools { - - /** + /** * Generate the point in a cuboid as an array const & */ - class grid_generator : - 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_; - double step_x, step_y, step_z; - bool at_end; - K_type pt; - void init() { - 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); - pt(0)=step_x/2; pt(1)=step_y/2; pt(2)=step_z/2; + class grid_generator : public boost::iterator_facade { + friend class boost::iterator_core_access; + int dim, nkpts, nx, ny, nz, N_X, N_Y, N_Z, index_; + double step_x, step_y, step_z; + bool at_end; + k_t pt; + void init() { + 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); + 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; } - - void increment() { - if (nx #include #include "grid_generator.hpp" -#include "functors.hpp" -namespace triqs { namespace lattice_tools { -using namespace std; - using namespace tqa; +namespace triqs { +namespace lattice { + + using namespace arrays; + + tight_binding::tight_binding(bravais_lattice const& bl, std::vector> all_disp, + std::vector> 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_view hopping_stack (tight_binding const & TB, array_view const & k_stack) { - auto TK = Fourier(TB); - array res(TB.n_bands(), TB.n_bands(), k_stack.shape(1)); - for(size_t i = 0; i hopping_stack(tight_binding const& TB, arrays::array_const_view k_stack) { + auto TK = fourier(TB); + array res(TB.n_bands(), TB.n_bands(), k_stack.shape(1)); + for (int i = 0; i < k_stack.shape(1); ++i) res(range(), range(), i) = TK(k_stack(range(), i)); return res; } //------------------------------------------------------ - - array_view energies_on_bz_path(tight_binding const & TB, K_view_type K1, K_view_type K2, size_t n_pts) { - auto TK = Fourier(TB); - const size_t norb=TB.lattice().n_orbitals(); - const size_t ndim=TB.lattice().dim(); - array eval(norb,n_pts); - K_type dk = (K2 - K1)/double(n_pts), k = K1; - for (size_t i =0; i energies_on_bz_path(tight_binding const& TB, k_t const& K1, k_t const& K2, int n_pts) { + auto TK = fourier(TB); + int norb = TB.lattice().n_orbitals(); + int ndim = TB.lattice().dim(); + array eval(norb, n_pts); + k_t dk = (K2 - K1) / double(n_pts), k = K1; + for (int i = 0; i < n_pts; ++i, k += dk) { + eval(range(), i) = linalg::eigenvalues(TK(k(range(0, ndim)))(), false); } return eval; } //------------------------------------------------------ - array_view energies_on_bz_grid(tight_binding const & TB, size_t n_pts) { + array energies_on_bz_grid(tight_binding const& TB, int n_pts) { - auto TK = Fourier(TB); - const size_t norb=TB.lattice().n_orbitals(); - const size_t ndim=TB.lattice().dim(); - grid_generator grid(ndim,n_pts); - array eval(norb,grid.size()); - for (; grid ; ++grid) { - eval(range(),grid.index()) = linalg::eigenvalues( TK( (*grid) (range(0,ndim))), false); + auto TK = fourier(TB); + int norb = TB.lattice().n_orbitals(); + int ndim = TB.lattice().dim(); + grid_generator grid(ndim, n_pts); + array eval(norb, grid.size()); + for (; grid; ++grid) { + eval(range(), grid.index()) = linalg::eigenvalues(TK((*grid)(range(0, ndim)))(), false); } return eval; } //------------------------------------------------------ - std::pair, array > dos(tight_binding const & TB, size_t nkpts, size_t neps) { + std::pair, array> dos(tight_binding const& TB, int nkpts, int neps) { - // The Fourier transform of TK - // auto TK = Fourier(TB); // C++0x .... - auto TK = Fourier(TB); + // The fourier transform of TK + auto TK = fourier(TB); // loop on the BZ - const size_t ndim=TB.lattice().dim(); - const size_t norb=TB.lattice().n_orbitals(); - grid_generator grid(ndim,nkpts); - array tempeval(norb); - array evec(norb,norb,grid.size()); - array eval(norb,grid.size()); - if (norb ==1) - for (; grid ; ++grid) { - double ee = real(TK( (*grid) (range(0,ndim)))(0,0)); - eval(0,grid.index()) =ee; - evec(0,0,grid.index()) =1; + int ndim = TB.lattice().dim(); + int norb = TB.lattice().n_orbitals(); + grid_generator grid(ndim, nkpts); + array tempeval(norb); + array evec(norb, norb, grid.size()); + array eval(norb, grid.size()); + if (norb == 1) + for (; grid; ++grid) { + double ee = real(TK((*grid)(range(0, ndim)))(0, 0)); + eval(0, grid.index()) = ee; + evec(0, 0, grid.index()) = 1; } - else - for (; grid ; ++grid) { - //cerr<<" index = "< eval_sl = eval(range(),grid.index()); - array_view evec_sl = evec(range(),range(),grid.index()); - std::tie (eval_sl,evec_sl) = linalg::eigenelements( TK( (*grid) (range(0,ndim)))); //, true); - //cerr<< " point "<< *grid << " value "<< eval_sl<< endl; //" "<< (*grid) (range(0,ndim)) << endl; + else + for (; grid; ++grid) { + // cerr<<" index = "< eval_sl = eval(range(), grid.index()); + array_view evec_sl = evec(range(), range(), grid.index()); + std::tie(eval_sl, evec_sl) = linalg::eigenelements(TK((*grid)(range(0, ndim)))); //, true); + // cerr<< " point "<< *grid << " value "<< eval_sl<< endl; //" "<< (*grid) (range(0,ndim)) << endl; } // define the epsilon mesh, etc. - array epsilon(neps); - double epsmax = tqa::max_element(eval); - double epsmin = tqa::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 (size_t i =0; i< neps; ++i) epsilon(i)=epsmin+(i+0.5)*deps; + array epsilon(neps); + double epsmax = max_element(eval); + double epsmin = min_element(eval); + double deps = (epsmax - epsmin) / neps; + // for (int 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 + 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 ...."< rho (neps,norb);rho()=0; - for(size_t l=0;l rho(neps, norb); + rho() = 0; + for (int l = 0; l < norb; l++) { + 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)); + // dos(a) += real(conj(evec(l,k,j))*evec(l,k,j)); } } } - //rho = rho / double(grid.size()*deps); - rho /= grid.size()*deps; - return std::make_pair( epsilon, rho); + // rho = rho / double(grid.size()*deps); + rho /= grid.size() * deps; + return std::make_pair(epsilon, rho); } //---------------------------------------------------------------------------------- - std::pair, array > dos_patch(tight_binding const & TB, const array & triangles, size_t neps, size_t ndiv) { + std::pair, array> dos_patch(tight_binding const& TB, const array& 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(); - //const size_t norb=TB.lattice().n_orbitals(); - int ntri = triangles.shape(0)/3; - array dos(neps); + // int ndim=TB.lattice().dim(); + // int norb=TB.lattice().n_orbitals(); + int ntri = triangles.shape(0) / 3; + array dos(neps); // Check consistency - const size_t ndim=TB.lattice().dim(); - //const size_t 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 "< thop(1,1); - array energ(nk), weight(nk); + int nk = ntri * ndiv * ndiv; + int k_index = 0; + double epsmax = -100000, epsmin = 100000; + array thop(1, 1); + array energ(nk), weight(nk); // a, b, c are the corners of the triangle // g the center of gravity taken from a - array a(ndim), b(ndim), c(ndim), g(ndim), rv(ndim); + array a(ndim), b(ndim), c(ndim), g(ndim), rv(ndim); int pt = 0; double s, t; - // The Fourier transform of TK - auto TK = Fourier(TB); + // The fourier transform of TK + auto TK = fourier(TB); // loop over the triangles for (int tri = 0; tri < ntri; tri++) { - a = triangles(pt,range()); + a = triangles(pt, range()); pt++; - b = triangles(pt,range()); + b = triangles(pt, range()); pt++; - c = triangles(pt,range()); + c = triangles(pt, range()); pt++; - g = ((a+b+c)/3.0-a)/double(ndiv); + g = ((a + b + c) / 3.0 - a) / double(ndiv); // the area around a k point might be different from one triangle to the other // 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)); + 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; iepsmax) epsmax=energ(k_index); - if (energ(k_index) epsmax) epsmax = energ(k_index); + if (energ(k_index) < epsmin) epsmin = energ(k_index); k_index++; } } @@ -196,27 +212,23 @@ using namespace std; assert(k_index == nk); // define the epsilon mesh, etc. - array epsilon(neps); - double deps=(epsmax-epsmin)/neps; - for (size_t i =0; i< neps; ++i) epsilon(i)= epsmin+i/(neps-1.0)*(epsmax-epsmin); + array epsilon(neps); + double deps = (epsmax - epsmin) / neps; + 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++) { - ind=int((energ(j)-epsmin)/deps); + 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); + dos /= deps; // Normalize the DOS + return {std::move(epsilon), std::move(dos)}; } - - - - - - }} +} +} diff --git a/triqs/lattice/tight_binding.hpp b/triqs/lattice/tight_binding.hpp index 24180e28..2f618bff 100644 --- a/triqs/lattice/tight_binding.hpp +++ b/triqs/lattice/tight_binding.hpp @@ -18,74 +18,82 @@ * TRIQS. If not, see . * ******************************************************************************/ -#ifndef TRIQS_LATTICE_TIGHTBINDINGS_H -#define TRIQS_LATTICE_TIGHTBINDINGS_H -#include "bravais_lattice_and_brillouin_zone.hpp" +#pragma once +#include "brillouin_zone.hpp" +#include -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. + Overlap between orbital is taken as unit matrix. */ - class tight_binding { - typedef std::vector, matrix>> displ_value_stack_t; - displ_value_stack_t displ_value_stack; - bravais_lattice bl_; - public : - typedef std::vector arg_type; + class tight_binding { - /// - tight_binding (bravais_lattice const & bl) : bl_(bl) {}; - - /// Underlying lattice - bravais_lattice const & lattice() const { return bl_;} + bravais_lattice bl_; + std::vector> all_disp; + std::vector> all_matrices; - /// Number of bands, i.e. size of the matrix t(k) - size_t n_bands() const { return bl_.n_orbitals();} - - /** - * Push_back mechanism of a pair displacement -> matrix - * VectorIntType is anything from which a std::vector can be constructed - * MatrixDComplexType is anything from which a matrix can be constructed - */ - template - void push_back(VectorIntType && v, MatrixDComplexType && m) { - std::vector V(std::forward(v)); - if (v.size() != bl_.dim()) TRIQS_RUNTIME_ERROR<<"tight_binding : displacement of incorrect size : got "<< v.size() << "instead of "<< bl_.dim(); - matrix M(std::forward(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))); - } + public: + /// + tight_binding(bravais_lattice const& bl, std::vector> all_disp, std::vector> all_matrices); - void reserve(size_t n) { displ_value_stack.reserve(n);} + /// Underlying lattice + bravais_lattice const& lattice() const { return bl_; } - // iterators - typedef displ_value_stack_t::const_iterator const_iterator; - typedef displ_value_stack_t::iterator iterator; - 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();} - iterator end() { return displ_value_stack.end();} - }; + /// Number of bands, i.e. size of the matrix t(k) + int n_bands() const { return bl_.n_orbitals(); } + + // calls F(R, t(R)) for all R + template friend void foreach(tight_binding const& tb, F f) { + int n = tb.all_disp.size(); + for (int i = 0; i < n; ++i) f(tb.all_disp[i], tb.all_matrices[i]); + } + + // a simple function on the domain brillouin_zone + struct fourier_impl { + tight_binding const& tb; + const int nb; + + using domain_t = brillouin_zone; + brillouin_zone domain() const { + return {tb.lattice()}; + } + + template matrix operator()(K const& k) { + matrix res(nb, nb); + res() = 0; + foreach(tb, [&](std::vector const& displ, matrix 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 hopping_stack (tight_binding const & TB, array_view const & k_stack); + array hopping_stack(tight_binding const& TB, arrays::array_const_view k_stack); // not optimal ordering here - std::pair, array > dos(tight_binding const & TB, size_t nkpts, size_t neps); - std::pair, array > dos_patch(tight_binding const & TB, const array & triangles, size_t neps, size_t ndiv); - array_view energies_on_bz_path(tight_binding const & TB, K_view_type K1, K_view_type K2, size_t n_pts); - array_view energies_on_bz_grid(tight_binding const & TB, size_t n_pts); - -}} - -#endif - + std::pair, array> dos(tight_binding const& TB, int nkpts, int neps); + std::pair, array> dos_patch(tight_binding const& TB, const array& triangles, int neps, + int ndiv); + array energies_on_bz_path(tight_binding const& TB, k_t const& K1, k_t const& K2, int n_pts); + array energies_on_bz_grid(tight_binding const& TB, int n_pts); +} +}