Lattice.hpp

Go to the documentation of this file.

 
#pragma once
#include <array>
#include <vector>
#include <tuple>
#include "lib/Geometry.hpp"
#include "lib/Assert.hpp"
 
struct SpinModel;
struct Lattice;
namespace LatticeModelFactory {
    struct LatticeUnitCell;
    struct SpinModelUnitCell;
 
    std::pair<Lattice *, SpinModel *> newLatticeModel(const LatticeModelFactory::LatticeUnitCell &uc, const LatticeModelFactory::SpinModelUnitCell &spinModelDefinition, const int latticeRange, const std::string &ldfPath);
};
 
enum struct SpinComponent : int
{
    X = 0, 
    Y = 1, 
    Z = 2, 
    None = 3 
};
 
struct LatticeOverlap
{
    friend std::pair<Lattice *, SpinModel *> LatticeModelFactory::newLatticeModel(const LatticeModelFactory::LatticeUnitCell &uc, const LatticeModelFactory::SpinModelUnitCell &spinModelDefinition, const int latticeRange, const std::string &ldfPath);
    friend struct Lattice;
 
public:
    LatticeOverlap() : rid1(nullptr), rid2(nullptr), transformedX1(nullptr), transformedY1(nullptr), transformedZ1(nullptr), transformedX2(nullptr), transformedY2(nullptr), transformedZ2(nullptr), size(0) {}
 
    LatticeOverlap(const int size) : size(size)
    {
        rid1 = new int[size];
        rid2 = new int[size];
        transformedX1 = new SpinComponent[size];
        transformedY1 = new SpinComponent[size];
        transformedZ1 = new SpinComponent[size];
        transformedX2 = new SpinComponent[size];
        transformedY2 = new SpinComponent[size];
        transformedZ2 = new SpinComponent[size];
    }
 
    LatticeOverlap(const LatticeOverlap &rhs) : LatticeOverlap(rhs.size)
    {
        memcpy(rid1, rhs.rid1, size * sizeof(int));
        memcpy(rid2, rhs.rid2, size * sizeof(int));
        memcpy(transformedX1, rhs.transformedX1, size * sizeof(SpinComponent));
        memcpy(transformedY1, rhs.transformedY1, size * sizeof(SpinComponent));
        memcpy(transformedZ1, rhs.transformedZ1, size * sizeof(SpinComponent));
        memcpy(transformedX2, rhs.transformedX2, size * sizeof(SpinComponent));
        memcpy(transformedY2, rhs.transformedY2, size * sizeof(SpinComponent));
        memcpy(transformedZ2, rhs.transformedZ2, size * sizeof(SpinComponent));
    }
 
    int *rid1;
    int *rid2;
    SpinComponent *transformedX1;
    SpinComponent *transformedY1;
    SpinComponent *transformedZ1;
    SpinComponent *transformedX2;
    SpinComponent *transformedY2;
    SpinComponent *transformedZ2;
    int size;
 
protected:
    ~LatticeOverlap()
    {
        delete[] rid1;
        delete[] rid2;
        delete[] transformedX1;
        delete[] transformedY1;
        delete[] transformedZ1;
        delete[] transformedX2;
        delete[] transformedY2;
        delete[] transformedZ2;
    }
 
    LatticeOverlap &operator=(const LatticeOverlap &rhs)
    {
        size = rhs.size;
        delete[] rid1;
        rid1 = new int[size];
        memcpy(rid1, rhs.rid1, size * sizeof(int));
        delete[] rid2;
        rid2 = new int[size];
        memcpy(rid2, rhs.rid2, size * sizeof(int));
        delete[] transformedX1;
        transformedX1 = new SpinComponent[size];
        memcpy(transformedX1, rhs.transformedX1, size * sizeof(SpinComponent));
        delete[] transformedY1;
        transformedY1 = new SpinComponent[size];
        memcpy(transformedY1, rhs.transformedY1, size * sizeof(SpinComponent));
        delete[] transformedZ1;
        transformedZ1 = new SpinComponent[size];
        memcpy(transformedZ1, rhs.transformedZ1, size * sizeof(SpinComponent));
        delete[] transformedX2;
        transformedX2 = new SpinComponent[size];
        memcpy(transformedX2, rhs.transformedX2, size * sizeof(SpinComponent));
        delete[] transformedY2;
        transformedY2 = new SpinComponent[size];
        memcpy(transformedY2, rhs.transformedY2, size * sizeof(SpinComponent));
        delete[] transformedZ2;
        transformedZ2 = new SpinComponent[size];
        memcpy(transformedZ2, rhs.transformedZ2, size * sizeof(SpinComponent));
 
        return *this;
    }
};
 
struct LatticeSiteDescriptor
{
    int rid; 
    SpinComponent spinPermutation[3]; 
};
 
struct LatticeIterator
{
friend struct Lattice;
 
public:
    LatticeIterator() : id(0) {}
 
    LatticeIterator(int id) : id(id) {}
 
    bool operator==(const LatticeIterator &rhs) const
    {
        return id == rhs.id;
    }
 
    bool operator!=(const LatticeIterator &rhs) const
    {
        return id != rhs.id;
    }
 
    int operator-(const LatticeIterator &rhs) const
    {
        return id - rhs.id;
    }
 
    virtual LatticeIterator& operator++()
    {
        ++id;
        return *this;
    }
 
    virtual LatticeIterator& operator--()
    {
        --id;
        return *this;
    }
 
protected:
    int id; 
};
 
struct SublatticeIterator : public LatticeIterator
{
public:
    SublatticeIterator(int *allowedIds) : LatticeIterator(*allowedIds), offset(0), allowedIds(allowedIds) {}
 
    SublatticeIterator &operator++()
    {
        id = allowedIds[++offset];
        return *this;
    }
 
    SublatticeIterator &operator--()
    {
        id = allowedIds[--offset];
        return *this;
    }
 
protected:
    int offset; 
    int *allowedIds; 
};
 
struct Lattice
{
    friend std::pair<Lattice *, SpinModel *> LatticeModelFactory::newLatticeModel(const LatticeModelFactory::LatticeUnitCell &uc, const LatticeModelFactory::SpinModelUnitCell &spinModelDefinition, const int latticeRange, const std::string &ldfPath);
 
protected:
    Lattice() : size(0), _dataSize(0), _symmetryTable(nullptr), _bufferSites(nullptr), _bufferInvertedSites(nullptr), _bufferOverlapMatrices(nullptr), _bufferBasis(nullptr), _bufferLatticeRange(nullptr) {};
 
public:
    ~Lattice()
    {
        delete[] _bufferBasis;
        for (size_t b = 0; b < _basis.size(); ++b) delete[] _bufferLatticeRange[b];
        delete[] _bufferLatticeRange;
        delete[] _symmetryTable;
        delete[] _bufferSites;
        delete[] _bufferInvertedSites;
        delete[] _bufferOverlapMatrices;
    }
 
    LatticeIterator begin() const
    {
        return LatticeIterator(0);
    }
 
    LatticeIterator zero() const
    {
        return LatticeIterator(0);
    }
 
    LatticeIterator end() const
    {
        return LatticeIterator(_dataSize);
    }
 
    LatticeIterator fromParametrization(const int rid) const
    {
        ASSERT(rid < size);
        return LatticeIterator(rid);
    }
 
    std::tuple<int, int, int, int> getSiteParameters(const LatticeIterator &site) const
    {
        ASSERT(site - begin() < _dataSize);
        return _geometryTable[site.id];
    }
 
    geometry::Vec3<double> getSitePosition(const LatticeIterator &site) const
    {
        ASSERT(site - begin() < _dataSize);
        return double(std::get<0>(_geometryTable[site.id])) * _bravaisLattice[0] + double(std::get<1>(_geometryTable[site.id])) * _bravaisLattice[1] + double(std::get<2>(_geometryTable[site.id])) * _bravaisLattice[2] + _basis[std::get<3>(_geometryTable[site.id])];
    }
 
    int symmetryTransform(const LatticeIterator& i1, const LatticeIterator& i2) const
    {
        ASSERT(_symmetryTable[i1.id * _dataSize + i2.id].rid != -1);
        return _symmetryTable[i1.id * _dataSize + i2.id].rid;
    }
 
    int symmetryTransform(const LatticeIterator &i1, const LatticeIterator &i2, SpinComponent &spinComponent) const
    {
        ASSERT(_symmetryTable[i1.id * _dataSize + i2.id].rid != -1);
        if (spinComponent != SpinComponent::None) spinComponent = _symmetryTable[i1.id * _dataSize + i2.id].spinPermutation[static_cast<int>(spinComponent)];
        return _symmetryTable[i1.id * _dataSize + i2.id].rid;
    }
 
    int symmetryTransform(const LatticeIterator &i1, const LatticeIterator &i2, SpinComponent &spinComponent1, SpinComponent &spinComponent2) const
    {
        ASSERT(_symmetryTable[i1.id * _dataSize + i2.id].rid != -1);
        ASSERT(&spinComponent1 != &spinComponent2);
 
        if (spinComponent1 != SpinComponent::None) spinComponent1 = _symmetryTable[i1.id * _dataSize + i2.id].spinPermutation[static_cast<int>(spinComponent1)];
        if (spinComponent2 != SpinComponent::None) spinComponent2 = _symmetryTable[i1.id * _dataSize + i2.id].spinPermutation[static_cast<int>(spinComponent2)];
        return _symmetryTable[i1.id * _dataSize + i2.id].rid;
    }
 
    int symmetryTransform(const LatticeIterator &i1, const LatticeIterator &i2, SpinComponent &spinComponent1, SpinComponent &spinComponent2, SpinComponent &spinComponent3) const
    {
        ASSERT(_symmetryTable[i1.id * _dataSize + i2.id].rid != -1);
        ASSERT(&spinComponent1 != &spinComponent2);
        ASSERT(&spinComponent2 != &spinComponent3);
        ASSERT(&spinComponent1 != &spinComponent3);
 
        if (spinComponent1 != SpinComponent::None) spinComponent1 = _symmetryTable[i1.id * _dataSize + i2.id].spinPermutation[static_cast<int>(spinComponent1)];
        if (spinComponent2 != SpinComponent::None) spinComponent2 = _symmetryTable[i1.id * _dataSize + i2.id].spinPermutation[static_cast<int>(spinComponent2)];
        if (spinComponent3 != SpinComponent::None) spinComponent3 = _symmetryTable[i1.id * _dataSize + i2.id].spinPermutation[static_cast<int>(spinComponent3)];
        return _symmetryTable[i1.id * _dataSize + i2.id].rid;
    }
 
    const LatticeOverlap &getOverlap(const int rid) const
    {
        ASSERT(rid < size);
        return _bufferOverlapMatrices[rid];
    }
 
    const LatticeSiteDescriptor *getInvertedSites() const
    {
        return _bufferInvertedSites;
    }
 
    const LatticeSiteDescriptor *getSites() const
    {
        return _bufferSites;
    }
 
    SublatticeIterator getRange(const LatticeIterator &b) const
    {
        ASSERT(std::get<0>(getSiteParameters(b)) == 0);
        ASSERT(std::get<1>(getSiteParameters(b)) == 0);
        ASSERT(std::get<2>(getSiteParameters(b)) == 0);
        ASSERT(std::get<3>(getSiteParameters(b)) < _basis.size());
        
        return SublatticeIterator(_bufferLatticeRange[std::get<3>(getSiteParameters(b))]);
    }
 
    SublatticeIterator getRange(const int b) const
    {
        ASSERT(b < _basis.size());
        return SublatticeIterator(_bufferLatticeRange[b]);
    }
 
    SublatticeIterator getBasis() const
    {
        return SublatticeIterator(_bufferBasis);
    }
 
    std::vector<geometry::Vec3<double> > _bravaisLattice; 
    std::vector<geometry::Vec3<double> > _basis; 
 
    int size; 
 
protected: 
    std::vector<std::tuple<int, int, int, int>> _geometryTable; 
 
    int _dataSize; 
    
    LatticeSiteDescriptor *_symmetryTable; 
    LatticeSiteDescriptor *_bufferSites; 
    LatticeSiteDescriptor *_bufferInvertedSites; 
    LatticeOverlap *_bufferOverlapMatrices; 
 
    int *_bufferBasis; 
    int **_bufferLatticeRange; 
};