1 #ifndef LIBVIPSTER_STEP_H
2 #define LIBVIPSTER_STEP_H
3 
4 #include "stepbase.h"
5 
6 #include <vector>
7 #include <memory>
8 
9 namespace Vipster {
10 
11 namespace detail{
12 
13 /*
14  * Basic serial Atom container
15  *
16  * Stores atom in separate vectors
17  */
18 struct AtomList{
19     template<bool isConst>
20     class AtomView;
21     using atom = AtomView<false>;
22     using const_atom = AtomView<true>;
23     template<bool isConst>
24     using AtomIterator = detail::AtomIterator<AtomView, isConst>;
25     using iterator = AtomIterator<false>;
26     using const_iterator = AtomIterator<true>;
27 
28     AtomContext ctxt;
29     // Coordinates
30     std::vector<Vec> coordinates{};
31     // Pointers to PeriodicTable entries
32     std::vector<PeriodicTable::value_type*> elements{};
33     // Properties
34     std::vector<AtomProperties> properties{};
35     // complete Step-interface
getNatAtomList36     size_t getNat() const noexcept {return elements.size();}
37 
AtomListAtomList38     AtomList(AtomFmt fmt)
39         : ctxt{fmt} {}
40     // specialize copy constructor because AtomList is supposed to own its cell
AtomListAtomList41     AtomList(const AtomList &rhs)
42         : ctxt{rhs.ctxt.fmt,
43                rhs.ctxt.pte,
44                std::make_shared<AtomContext::CellData>(*rhs.ctxt.cell)},
45           coordinates{rhs.coordinates},
46           elements{rhs.elements},
47           properties{rhs.properties}
48     {}
49     AtomList(AtomList&&) = default;
50     AtomList& operator=(const AtomList&) = default;
51     AtomList& operator=(AtomList&&) = default;
52     // enable creation from other Atom-sources
53     template<typename T>
AtomListAtomList54     AtomList(const T &rhs)
55         : ctxt{rhs.ctxt.fmt,
56                rhs.ctxt.pte,
57                std::make_shared<AtomContext::CellData>(*rhs.ctxt.cell)},
58           coordinates{rhs.coordinates},
59           elements{rhs.elements},
60           properties{rhs.properties}
61     {}
62 
63     template<bool isConst>
64     class AtomView
65     {
66         /* Wrapper class that wraps distributed SoA data into a singular Atom-interface
67          *
68          * should behave like a reference
69          */
70         Vec *v;
71         PeriodicTable::value_type **elem;
72         AtomProperties *prop;
73 
AtomViewAtomList74         AtomView()
75             : source{nullptr},
76               v{nullptr},
77               elem{nullptr},
78               prop{nullptr}
79         {}
80         template<bool> friend class AtomView;
81         class _Vec{
82         public:
_VecAtomList83             _Vec(AtomView &a): a{a} {}
84             // can only be explicitely constructed to reference an Atom
85             _Vec(const _Vec&) = delete;
86             // assigning changes the origin
87             _Vec& operator=(const _Vec& rhs){
88                 return operator=(static_cast<const Vec&>(rhs));
89             }
90             _Vec& operator=(const Vec& rhs){
91                 *a.v = rhs;
92                 return *this;
93             }
94             // convert to regular Vec-reference
95             operator Vec&() {return *a.v;}
96             operator const Vec&() const {return *a.v;}
97             // formatted Vec access
asFmtAtomList98             Vec asFmt(const AtomContext &ctxt) const
99             {
100                 return makeConverter(a.source->ctxt, ctxt)(static_cast<const Vec&>(*this));
101             }
102             // array access
103             Vec::value_type& operator[](std::size_t i){
104                 return (*a.v)[i];
105             }
106             const Vec::value_type& operator[](std::size_t i) const {
107                 return (*a.v)[i];
108             }
109             // comparison
110             bool operator==(const Vec& rhs) const {return *a.v == rhs;}
111         private:
112             AtomView &a;
113         };
114         class _Name{
115         public:
_NameAtomList116             _Name(AtomView &a): a{a} {}
117             // can only be explicitely constructed to reference an Atom
118             _Name(const _Name&) =delete;
119             // assigning a string makes this point to the
120             // appropriate entry in the periodic table
121             _Name& operator=(const _Name& rhs){
122                 *a.elem = &*a.source->ctxt.pte->find_or_fallback(rhs);
123                 return *this;
124             }
125             _Name& operator=(const std::string& rhs){
126                 *a.elem = &*a.source->ctxt.pte->find_or_fallback(rhs);
127                 return *this;
128             }
129             // convert to regular string-reference
130             operator const std::string&() const {return (*a.elem)->first;}
string_viewAtomList131             operator std::string_view() const {return (*a.elem)->first;}
132             // comparison
133             bool operator==(const std::string& rhs) const {return (*a.elem)->first == rhs;}
134             // const char* access
c_strAtomList135             const char* c_str() const{
136                 return (*a.elem)->first.c_str();
137             }
138             // I/O
139             friend std::ostream& operator<<(std::ostream &s, const _Name &n){
140                 return s << n.c_str();
141             }
142             friend std::istream& operator>>(std::istream &s, _Name &n){
143                 std::string tmp;
144                 s >> tmp;
145                 n = tmp;
146                 return s;
147             }
148         private:
149             AtomView &a;
150         };
151         class _Element{
152         public:
_ElementAtomList153             _Element(AtomView &a): a{a} {}
154             // can only be explicitely constructed to reference an Atom
155             _Element(const _Element& at) =delete;
156             // Allow member-access via pointer-syntax
157             Element* operator->() { return &(*a.elem)->second;}
158             const Element* operator->() const{return &(*a.elem)->second;}
159             // convert to reference
160             operator const Element&() const {return (*a.elem)->second;}
161             // comparison
162             bool operator==(const Element& rhs) const { return (*a.elem)->second == rhs;}
163         private:
164             AtomView &a;
165         };
166         class _Properties{
167         public:
_PropertiesAtomList168             _Properties(AtomView &a): a{a} {}
169             // can only be explicitely constructed to reference an Atom
170             _Properties(const _Properties&) =delete;
171             // like real reference, assigning changes the origin
172             _Properties& operator=(const _Properties& rhs){
173                 *a.prop = static_cast<const AtomProperties&>(rhs);
174                 return *this;
175             }
176             _Properties& operator=(const AtomProperties& rhs){
177                 *a.prop = rhs;
178                 return *this;
179             }
180             // convert to regular AtomProperties-reference
181             operator const AtomProperties&() const {return *a.prop;}
182             // comparison
183             bool operator==(const AtomProperties &rhs) const {return *a.prop == rhs;}
184             // Allow pointer-syntax
185             const AtomProperties* operator->() const {return a.prop;}
186             AtomProperties* operator->() {return a.prop;}
187         private:
188             AtomView &a;
189         };
190     protected:
191         using Source = AtomList;
192         Source *source;
193         AtomView& operator+=(ptrdiff_t i){
194             v += i;
195             elem += i;
196             prop += i;
197             return *this;
198         }
pointToAtomList199         void pointTo(const AtomView& rhs){
200             source = rhs.source;
201             v = rhs.v;
202             elem = rhs.elem;
203             prop = rhs.prop;
204         }
205     public:
206         // "Data", encapsulated in wrapper objects
207         std::conditional_t<isConst, const _Vec, _Vec> coord{*this};
208         std::conditional_t<isConst, const _Name, _Name> name{*this};
209         std::conditional_t<isConst, const _Element, _Element> type{*this};
210         std::conditional_t<isConst, const _Properties, _Properties> properties{*this};
211 
AtomViewAtomList212         AtomView(AtomList &al,
213                  size_t i)
214             : source{&al},
215               v{&al.coordinates[i]},
216               elem{&al.elements[i]},
217               prop{&al.properties[i]}
218         {}
219         // copying is templated to allow conversion to const
220         // copy constructor creates new object pointing to same data
AtomViewAtomList221         AtomView(const AtomView &rhs)
222             : source{rhs.source},
223               v{rhs.v},
224               elem{rhs.elem},
225               prop{rhs.prop}
226         {}
227         template<bool B, bool t=isConst, typename = typename std::enable_if<t>::type>
AtomViewAtomList228         AtomView(const AtomView<B> &rhs)
229             : source{rhs.source},
230               v{rhs.v},
231               elem{rhs.elem},
232               prop{rhs.prop}
233         {}
234         // copy assignment changes data
235         AtomView& operator=(const AtomView &rhs){
236             coord = rhs.coord.asFmt(source->ctxt);
237             name = rhs.name;
238             properties = rhs.properties;
239             return *this;
240         }
241         template<bool B, bool t=isConst, typename = typename std::enable_if<t>::type>
242         AtomView& operator=(const AtomView<B> &rhs){
243             coord = rhs.coord.asFmt(source->ctxt);
244             name = rhs.name;
245             properties = rhs.properties;
246             return *this;
247         }
248         template<bool B, template<bool> typename AV>
249         AtomView& operator=(const AV<B> &rhs){
250             coord = rhs.coord.asFmt(source->ctxt);
251             name = rhs.name;
252             properties = rhs.properties;
253             return *this;
254         }
255         virtual ~AtomView() = default;
256 
257         bool operator==(const AtomView &rhs) const noexcept{
258             return std::tie(coord, name, properties)
259                     ==
260                    std::tie(rhs.coord, rhs.name, rhs.properties);
261         }
262         bool operator!=(const AtomView &rhs) const noexcept{
263             return !operator==(rhs);
264         }
265     };
266 };
267 
268 }
269 
270 /*
271  * Main Step-class
272  *
273  * with AtomList as source, this is the main class to use for atom-storage
274  */
275 
276 class Step: public StepMutable<detail::AtomList>
277 {
278 private:
279     friend class Molecule;
280     // Periodic table
281     void setPTE(std::shared_ptr<PeriodicTable> newPTE);
282 public:
283     Step(AtomFmt at_fmt=AtomFmt::Angstrom,
284          const std::string &comment="");
285     Step(const Step& s);
286     Step(Step&& s);
287     Step& operator=(const Step& s);
288     Step& operator=(Step&& s);
289     template<typename T>
Step(const StepConst<T> & s)290     Step(const StepConst<T>& s)
291         : StepMutable{std::make_shared<detail::AtomList>(s.getFmt()),
292                       std::make_shared<BondList>(),
293                       std::make_shared<std::string>(s.getComment())}
294     {
295         // copy Cell
296         if(s.hasCell()){
297             setCellDim(s.getCellDim(AtomFmt::Bohr), AtomFmt::Bohr);
298             setCellVec(s.getCellVec());
299         }
300         // copy Atoms
301         newAtoms(s);
302         // copy Bonds
303         for(const auto& b: s.getBonds()){
304             addBond(b.at1, b.at2, b.diff, b.type ? b.type->first : "");
305         }
306     }
307 
308     // Format
309     void setFmt(AtomFmt fmt, bool scale=true);
310 
311     // Atoms
312     void    newAtom(std::string name="",
313                     Vec coord=Vec{},
314                     AtomProperties prop=AtomProperties{});
315     template<template<bool> typename T, bool B>
newAtom(const T<B> & at)316     void    newAtom(const T<B>& at){
317         newAtom();
318         back() = at;
319     }
320     void    newAtoms(size_t i);
321     template<typename T>
newAtoms(const StepConst<T> & s)322     void    newAtoms(const StepConst<T>& s)
323     {
324         const size_t oldNat = this->getNat();
325         newAtoms(s.getNat());
326         auto tgt = begin()+oldNat;
327         auto src = s.begin();
328         for(; tgt != end(); ++tgt, ++src){
329             *tgt = *src;
330         }
331     }
332     void    delAtom(size_t i);
333     template<typename T>
delAtoms(StepConst<detail::Selection<T>> & s)334     void delAtoms(StepConst<detail::Selection<T>>& s)
335     {
336         std::set<size_t> indices{};
337         for(const auto [idx, _]: s.getAtoms().indices){
338             indices.insert(idx);
339         }
340         for(auto it = indices.rbegin(); it != indices.rend(); ++it)
341         {
342             if(*it < getNat()){
343                 delAtom(*it);
344             }
345         }
346         s = select({});
347     }
348 
349     // Cell
350     void enableCell(bool val) noexcept;
351     void setCellDim(double cdm, AtomFmt fmt, bool scale=false);
352     void setCellVec(const Mat &vec, bool scale=false);
353 
354     // Modifier functions
355     void modWrap();
356     void modCrop();
357     void modMultiply(size_t x, size_t y, size_t z);
358     void modAlign(uint8_t step_dir, uint8_t target_dir);
359     void modReshape(Mat newMat, double newDim, AtomFmt Fmt);
360 };
361 
362 }
363 #endif // LIBVIPSTER_STEP_H
364