Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2025 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "AtomDistribution.h"
23 : #include "Exception.h"
24 : #include "View.h"
25 : #include "Tools.h"
26 : #include <cassert>
27 : #include <variant>
28 : #include <optional>
29 : #include <functional>
30 :
31 :
32 : namespace PLMD {
33 :
34 : using unpackedLine=gch::small_vector<std::string_view>;
35 :
36 196 : void AtomDistribution::frame(std::vector<Vector>& posToUpdate,
37 : std::vector<double>& box,
38 : unsigned step,
39 : Random& rng) {
40 : assert(box.size()>=9);
41 : //call the specialized frame
42 196 : frame(make_view(posToUpdate),
43 : View<double,9>(box.data()),
44 : step,
45 : rng);
46 196 : }
47 :
48 : class UniformSphericalVector {
49 : double rCub;
50 :
51 : public:
52 : //assuming rmin=0
53 40 : explicit UniformSphericalVector(const double rmax):
54 40 : rCub (rmax*rmax*rmax/*-rminCub*/) {}
55 7000 : PLMD::Vector operator()(Random& rng) const {
56 7000 : double rho = std::cbrt (/*rminCub + */rng.RandU01()*rCub);
57 7000 : double theta =std::acos (2.0*rng.RandU01() -1.0);
58 7000 : double phi = 2.0 * PLMD::pi * rng.RandU01();
59 : return Vector (
60 7000 : rho * sin (theta) * cos (phi),
61 7000 : rho * sin (theta) * sin (phi),
62 7000 : rho * cos (theta));
63 : }
64 : };
65 :
66 : ///A wiggly line of atoms
67 1 : struct theLine:public AtomDistribution {
68 : static constexpr auto id="line";
69 : static constexpr auto doc="A line of atoms, that wiggles randomly at each step";
70 1 : void frame(View<Vector> posToUpdate,
71 : View<double,9> box,
72 : unsigned /*step*/,
73 : Random& /*rng*/) override {
74 : auto nat = posToUpdate.size();
75 :
76 201 : for (unsigned i=0; i<nat; ++i) {
77 200 : posToUpdate[i] = Vector(i, 0, 0);
78 : }
79 1 : box[0]=nat;
80 1 : box[1]=0.0;
81 1 : box[2]=0.0;
82 1 : box[3]=0.0;
83 1 : box[4]=1.75;
84 1 : box[5]=0.0;
85 1 : box[6]=0.0;
86 1 : box[7]=0.0;
87 1 : box[8]=1.75;
88 1 : }
89 : };
90 :
91 : ///Atom randomly distributed in a sphere
92 17 : struct uniformSphere:public AtomDistribution {
93 : static constexpr auto id="sphere";
94 : static constexpr auto doc="Atoms are displaced uniformly in a sphere, the desity is ccirca 1/unit^3";
95 31 : void frame(View<Vector> posToUpdate,
96 : View<double,9> box,
97 : unsigned /*step*/,
98 : Random& rng) override {
99 : //giving more or less a cubic udm of volume for each atom: V=nat
100 31 : const double rmax= std::cbrt ((3.0/(4.0*PLMD::pi)) * posToUpdate.size());
101 :
102 : UniformSphericalVector usv(rmax);
103 : auto s=posToUpdate.begin();
104 : auto e=posToUpdate.end();
105 : //I am using the iterators:this is slightly faster,
106 : // enough to overcome the cost of the vtable that I added
107 5231 : for (unsigned i=0; s!=e; ++s,++i) {
108 5200 : *s = usv (rng);
109 : }
110 :
111 31 : box[0]=2.0*rmax;
112 31 : box[1]=0.0;
113 31 : box[2]=0.0;
114 31 : box[3]=0.0;
115 31 : box[4]=2.0*rmax;
116 31 : box[5]=0.0;
117 31 : box[6]=0.0;
118 31 : box[7]=0.0;
119 31 : box[8]=2.0*rmax;
120 31 : }
121 : };
122 :
123 : ///Atom randomly distributed between two not overlapping a spheres
124 5 : struct twoGlobs: public AtomDistribution {
125 : static constexpr auto id="globs";
126 : static constexpr auto doc="Two spheres of uniformly distributed atosm";
127 8 : void frame(View<Vector> posToUpdate,
128 : View<double,9> box,
129 : unsigned /*step*/,
130 : Random&rng) override {
131 : //I am using two unigform spheres and 2V=n
132 8 : const double rmax= std::cbrt ((3.0/(8.0*PLMD::pi)) * posToUpdate.size());
133 :
134 : UniformSphericalVector usv(rmax);
135 8 : const std::array<const Vector,2> centers{
136 : PLMD::Vector{0.0,0.0,0.0},
137 : //so they do not overlap
138 : PLMD::Vector{2.0*rmax,2.0*rmax,2.0*rmax}
139 8 : };
140 1608 : std::generate(posToUpdate.begin(),posToUpdate.end(),[&]() {
141 : //RandInt is only declared
142 : // return usv (rng) + centers[rng.RandInt(1)];
143 1600 : return usv (rng) + centers[rng.RandU01()>0.5];
144 : });
145 :
146 8 : box[0]=4.0 *rmax;
147 8 : box[1]=0.0;
148 8 : box[2]=0.0;
149 8 : box[3]=0.0;
150 8 : box[4]=4.0 *rmax;
151 8 : box[5]=0.0;
152 8 : box[6]=0.0;
153 8 : box[7]=0.0;
154 8 : box[8]=4.0 *rmax;
155 8 : }
156 : };
157 :
158 2 : struct uniformCube:public AtomDistribution {
159 : static constexpr auto id="cube";
160 : static constexpr auto doc="Randomly displaced atoms in a cube, the density is circa 1/unit^3";
161 2 : void frame(View<Vector> posToUpdate,
162 : View<double,9> box,
163 : unsigned /*step*/,
164 : Random& rng) override {
165 : //giving more or less a cubic udm of volume for each atom: V = nat
166 2 : const double rmax = std::cbrt(static_cast<double>(posToUpdate.size()));
167 :
168 : auto s=posToUpdate.begin();
169 : auto e=posToUpdate.end();
170 : //I am using the iterators:this is slightly faster,
171 : // enough to overcome the cost of the vtable that I added
172 402 : for (unsigned i=0; s!=e; ++s,++i) {
173 400 : *s = Vector (rng.RandU01()*rmax,rng.RandU01()*rmax,rng.RandU01()*rmax);
174 : }
175 : //+0.05 to avoid overlap
176 2 : box[0]=rmax+0.05;
177 2 : box[1]=0.0;
178 2 : box[2]=0.0;
179 2 : box[3]=0.0;
180 2 : box[4]=rmax+0.05;
181 2 : box[5]=0.0;
182 2 : box[6]=0.0;
183 2 : box[7]=0.0;
184 2 : box[8]=rmax+0.05;
185 2 : }
186 : };
187 :
188 : /// A simple correction to getting (in ceil mode) the correct nearest perfect cube
189 146 : unsigned ceiledPerfectCube(const unsigned N) {
190 : // In some case a perfect cube is not well represented by double
191 : // and will get uncorrectly ceiled up (for example 24389 will result in a 30 and not in 29)
192 :
193 : // This assumes that N!=0
194 146 : const size_t x = std::ceil(std::cbrt(static_cast<double>(N)))-1;
195 146 : if ( x*x*x >= N) {
196 22 : return x;
197 : }
198 124 : return x+1;
199 : }
200 :
201 92 : struct tiledSimpleCubic:public AtomDistribution {
202 : static constexpr auto id="sc";
203 : static constexpr auto doc="Simple Cubic";
204 139 : void frame(View<Vector> posToUpdate,
205 : View<double,9> box,
206 : unsigned /*step*/,
207 : Random& rng) override {
208 : //Tiling the space in this way will not tests 100% the pbc, but
209 : //I do not think that write a spacefilling curve, like Hilbert, Peano or Morton
210 : //could be a good idea, in this case
211 :
212 139 : const unsigned rmax =ceiledPerfectCube(posToUpdate.size());
213 :
214 : auto s=posToUpdate.begin();
215 : auto e=posToUpdate.end();
216 : //I am using the iterators:this is slightly faster,
217 : // enough to overcome the cost of the vtable that I added
218 759 : for (unsigned k=0; k<rmax&&s!=e; ++k) {
219 3464 : for (unsigned j=0; j<rmax&&s!=e; ++j) {
220 16258 : for (unsigned i=0; i<rmax&&s!=e; ++i) {
221 13414 : *s = Vector (i,j,k);
222 13414 : ++s;
223 : }
224 : }
225 : }
226 139 : box[0]=rmax;
227 139 : box[1]=0.0;
228 139 : box[2]=0.0;
229 139 : box[3]=0.0;
230 139 : box[4]=rmax;
231 139 : box[5]=0.0;
232 139 : box[6]=0.0;
233 139 : box[7]=0.0;
234 139 : box[8]=rmax;
235 139 : }
236 : };
237 :
238 : /// This FCC is contained in a orthogonal cubic box
239 : ///
240 : /// Fun facts: full cubes at 4, 32, 108, 256, 500 atoms and so on.
241 : ///
242 : /// At certain sizes above the ones above (36, 504), you get
243 : /// good (100) slabs quite separated on th z axis
244 2 : struct inscribedFaceCenteredCubic:public AtomDistribution {
245 : static constexpr auto id="ifcc";
246 : static constexpr auto doc="FCC, but in a cubic box";
247 2 : void frame(View<Vector> posToUpdate,
248 : View<double,9> box,
249 : unsigned /*step*/,
250 : Random& rng) override {
251 : //Here we are exploiting a little trick:
252 : // if you remove the even or odd atoms from a simple cubic built with our algorithm,
253 : //you get an fcc boxed in a cube
254 2 : const unsigned rmax = [&] {
255 2 : auto x = ceiledPerfectCube(2*posToUpdate.size());
256 : //rmax needs to be even for this trick to work
257 2 : if ( x%2==0) {
258 : return x;
259 : } else {
260 0 : return x+1;
261 : }
262 :
263 2 : }();
264 :
265 : auto s=posToUpdate.begin();
266 : auto e=posToUpdate.end();
267 : //I am using the iterators:this is slightly faster,
268 : // enough to overcome the cost of the vtable that I added
269 16 : for (unsigned k=0; k<rmax&&s!=e; ++k) {
270 114 : for (unsigned j=0; j<rmax&&s!=e; ++j) {
271 : //we choose to show only the atoms with even index: (i+j+k)%2
272 : //Like this we skip steps and lots of divisions
273 500 : for (unsigned i=(j+k)%2; i<rmax&&s!=e; i+=2) {
274 400 : *s = 0.5*Vector (i,j,k);
275 400 : ++s;
276 : }
277 : }
278 : }
279 2 : box[0]=0.5*rmax;
280 2 : box[1]=0.0;
281 2 : box[2]=0.0;
282 2 : box[3]=0.0;
283 2 : box[4]=0.5*rmax;
284 2 : box[5]=0.0;
285 2 : box[6]=0.0;
286 2 : box[7]=0.0;
287 2 : box[8]=0.5*rmax;
288 2 : }
289 : };
290 :
291 : /// This FCC is contained in a "standard" FCC box
292 3 : struct tiledFaceCenteredCubic:public AtomDistribution {
293 : static constexpr auto id="fcc";
294 : static constexpr auto doc="Face Centered Cubic";
295 3 : void frame(View<Vector> posToUpdate,
296 : View<double,9> box,
297 : unsigned /*step*/,
298 : Random& rng) override {
299 :
300 3 : const unsigned rmax =ceiledPerfectCube(posToUpdate.size());
301 :
302 : auto s=posToUpdate.begin();
303 : auto e=posToUpdate.end();
304 : #define X PLMD::Versors::xp<double>
305 : #define Y PLMD::Versors::yp<double>
306 : #define Z PLMD::Versors::zp<double>
307 : const auto a=sqrt(2)*0.5*(X+Y);
308 : const auto b=sqrt(2)*0.5*(X+Z);
309 : const auto c=sqrt(2)*0.5*(Y+Z);
310 : #undef X
311 : #undef Y
312 : #undef Z
313 : //I am using the iterators:this is slightly faster,
314 : // enough to overcome the cost of the vtable that I added
315 21 : for (unsigned k=0; k<rmax&&s!=e; ++k) {
316 120 : for (unsigned j=0; j<rmax&&s!=e; ++j) {
317 702 : for (unsigned i=0; i<rmax&&s!=e; ++i) {
318 : *s = i*a
319 : +j*b
320 600 : +k*c;
321 600 : ++s;
322 : }
323 : }
324 : }
325 :
326 : box.subview<0,3>() = a*rmax;
327 : box.subview<3,3>() = b*rmax;
328 : box.subview<6,3>() = c*rmax;
329 3 : }
330 : };
331 :
332 : /// This BCC is contained in a "standard" FCC box
333 1 : struct tiledBodyCenteredCubic:public AtomDistribution {
334 : static constexpr auto id="bcc";
335 : static constexpr auto doc="Body Centered Cubic";
336 1 : void frame(View<Vector> posToUpdate,
337 : View<double,9> box,
338 : unsigned /*step*/,
339 : Random& rng) override {
340 : //For the base functionality of this see the comment in the tiledSimpleCubic
341 :
342 1 : const unsigned rmax =ceiledPerfectCube(posToUpdate.size());
343 :
344 : auto s=posToUpdate.begin();
345 : auto e=posToUpdate.end();
346 : #define X PLMD::Versors::xp<double>
347 : #define Y PLMD::Versors::yp<double>
348 : #define Z PLMD::Versors::zp<double>
349 : const auto a=0.5*(-X+Y+Z);
350 : const auto b=0.5*( X-Y+Z);
351 : const auto c=0.5*( X+Y-Z);
352 : #undef X
353 : #undef Y
354 : #undef Z
355 : //I am using the iterators:this is slightly faster,
356 : // enough to overcome the cost of the vtable that I added
357 7 : for (unsigned k=0; k<rmax&&s!=e; ++k) {
358 40 : for (unsigned j=0; j<rmax&&s!=e; ++j) {
359 234 : for (unsigned i=0; i<rmax&&s!=e; ++i) {
360 : *s = i*a
361 : +j*b
362 200 : +k*c;
363 200 : ++s;
364 : }
365 : }
366 : }
367 :
368 : box.subview<0,3>() = a*rmax;
369 : box.subview<3,3>() = b*rmax;
370 : box.subview<6,3>() = c*rmax;
371 1 : }
372 : };
373 :
374 : /// This BCC is contained in a orthogonal cubic box
375 1 : struct inscribedBodyCenteredCubic:public AtomDistribution {
376 : static constexpr auto id="ibcc";
377 : static constexpr auto doc="BCC, but in a cubic box";
378 1 : void frame(View<Vector> posToUpdate,
379 : View<double,9> box,
380 : unsigned /*step*/,
381 : Random& rng) override {
382 : //Here we are exploiting a little trick:
383 : // if you remove the even or odd atoms from a simple cubic built with our algorithm,
384 : //you get an bcc boxed in a cube
385 1 : const unsigned rmax = [&] {
386 : //we take one atoms every 4
387 1 : auto x = ceiledPerfectCube(4*posToUpdate.size());
388 : //rmax needs to be even for this trick to work
389 1 : if ( x%2==0) {
390 : return x;
391 : } else {
392 0 : return x+1;
393 : }
394 :
395 1 : }();
396 :
397 : auto s=posToUpdate.begin();
398 : auto e=posToUpdate.end();
399 : //I am using the iterators:this is slightly faster,
400 : // enough to overcome the cost of the vtable that I added
401 9 : for (unsigned k=0; k<rmax&&s!=e; ++k) {
402 : //we alternate on the z axis the choice between atoms with both i and j even or odd
403 : //The +=2 skips an extra check in the inner body
404 48 : for (unsigned j=k%2; j<rmax&&s!=e; j+=2) {
405 240 : for (unsigned i=k%2; i<rmax&&s!=e; i+=2) {
406 200 : *s = 0.5*Vector (i,j,k);
407 200 : ++s;
408 : }
409 : }
410 : }
411 1 : box[0]=0.5*rmax;
412 1 : box[1]=0.0;
413 1 : box[2]=0.0;
414 1 : box[3]=0.0;
415 1 : box[4]=0.5*rmax;
416 1 : box[5]=0.0;
417 1 : box[6]=0.0;
418 1 : box[7]=0.0;
419 1 : box[8]=0.5*rmax;
420 1 : }
421 : };
422 :
423 :
424 : ///a decorator for replicating the atomic distribution
425 : class repliedTrajectory: public AtomDistribution {
426 : std::unique_ptr<AtomDistribution> distribution;
427 : unsigned rX=1;
428 : unsigned rY=1;
429 : unsigned rZ=1;
430 : public:
431 : static constexpr auto id="reply";
432 : static constexpr auto doc=R"=(replicate the box by the given in all the directions:
433 : usage:
434 : reply x y z)=";
435 :
436 16 : static std::unique_ptr<AtomDistribution> decorate(std::unique_ptr<AtomDistribution>&& d,
437 : std::string_view cmd) {
438 : unpackedLine lines;
439 16 : Tools::getWordsSimple(lines,cmd);
440 : unsigned repeatX;
441 : unsigned repeatY;
442 : unsigned repeatZ;
443 16 : plumed_assert(lines.size() ==4) << id << " supports exacly three inputs: x y z";
444 :
445 32 : Tools::convert(std::string(lines[1]),repeatX);
446 32 : Tools::convert(std::string(lines[2]),repeatY);
447 16 : Tools::convert(std::string(lines[3]),repeatZ);
448 16 : return std::make_unique<repliedTrajectory>(std::move(d),
449 : repeatX,
450 : repeatY,
451 16 : repeatZ);
452 : }
453 :
454 : repliedTrajectory(std::unique_ptr<AtomDistribution>&& d,
455 : const unsigned repeatX,
456 : const unsigned repeatY,
457 16 : const unsigned repeatZ) :
458 : distribution(std::move(d)),
459 16 : rX(repeatX),
460 16 : rY(repeatY),
461 16 : rZ(repeatZ)
462 : {}
463 :
464 : ///Generates the inner trajectory of `(nat)/(rX*rY*rZ)`, where `nat` is the dimension of the vector view, then it replicate the atoms in the various directions
465 19 : void frame(View<Vector> posToUpdate,
466 : View<double,9> box,
467 : unsigned step,
468 : Random& rng) override {
469 19 : plumed_assert((posToUpdate.size() % (rX*rY*rZ)==0));
470 19 : const auto nat = posToUpdate.size()/(rX*rY*rZ);
471 19 : auto coordinates=posToUpdate.subview(0,nat);
472 19 : distribution->frame(coordinates,box,step,rng);
473 :
474 : // repetitions
475 19 : auto p = posToUpdate.begin()+nat;
476 :
477 : assert((rX*rY*rZ)*nat == posToUpdate.size());
478 :
479 19 : Vector boxX(box[0],box[1],box[2]);
480 19 : Vector boxY(box[3],box[4],box[5]);
481 19 : Vector boxZ(box[6],box[7],box[8]);
482 53 : for (unsigned x=0; x<rX; ++x) {
483 85 : for (unsigned y=0; y<rY; ++y) {
484 201 : for (unsigned z=0; z<rZ; ++z) {
485 150 : if(x==0&&y==0&&z==0) {
486 19 : continue;
487 : }
488 23531 : for (unsigned i=0; i<nat; ++i) {
489 : *p=coordinates[i]
490 : + x * boxX
491 : + y * boxY
492 23400 : + z * boxZ;
493 23400 : ++p;
494 : }
495 : }
496 : }
497 : }
498 19 : box[0]*=rX;
499 19 : box[1]*=rX;
500 19 : box[2]*=rX;
501 :
502 19 : box[3]*=rY;
503 19 : box[4]*=rY;
504 19 : box[5]*=rY;
505 :
506 19 : box[6]*=rZ;
507 19 : box[7]*=rZ;
508 19 : box[8]*=rZ;
509 19 : }
510 :
511 : ///See the documentation the AtomDistribution
512 12 : bool overrideNat(unsigned& natoms) override {
513 12 : distribution->overrideNat(natoms);
514 12 : natoms *= (rX*rY*rZ);
515 12 : return true;
516 : }
517 : };
518 :
519 : ///a decorator for scaling the atomic positions
520 : class scaledTrajectory: public AtomDistribution {
521 : std::unique_ptr<AtomDistribution> distribution;
522 : double multiplier;
523 : public:
524 : static constexpr auto id="scale";
525 : static constexpr auto doc=R"=(scales the atoms by a certain amount
526 : usage:
527 : scale (optional: mults=2.0))=";
528 :
529 12 : static std::unique_ptr<AtomDistribution> decorate(std::unique_ptr<AtomDistribution>&& d,
530 : std::string_view cmd) {
531 : unpackedLine lines;
532 12 : Tools::getWordsSimple(lines,cmd);
533 12 : plumed_assert(lines.size() <=2) << id << " supports maximum one input";
534 12 : if (lines.size()==2) {
535 12 : double mult=2.0;
536 12 : Tools::convert(std::string(lines[1]), mult);
537 12 : return std::make_unique<scaledTrajectory>(std::move(d),mult);
538 : }
539 :
540 : //the default value is specified in the header
541 0 : return std::make_unique<scaledTrajectory>(std::move(d));
542 : }
543 :
544 : scaledTrajectory(std::unique_ptr<AtomDistribution>&& d,
545 12 : double mult=2.0):
546 : distribution(std::move(d)),
547 12 : multiplier(mult)
548 : {}
549 :
550 12 : void frame(View<Vector> posToUpdate,
551 : View<double,9> box,
552 : unsigned step,
553 : Random& rng) override {
554 12 : distribution->frame(posToUpdate,box,step,rng);
555 :
556 2412 : for (auto& p:posToUpdate) {
557 2400 : p*=multiplier;
558 : }
559 :
560 120 : for (auto& b:box) {
561 108 : b*=multiplier;
562 : }
563 12 : }
564 :
565 0 : bool overrideNat(unsigned& natoms) override {
566 0 : return distribution->overrideNat(natoms);
567 : }
568 : };
569 :
570 :
571 :
572 : /// A decorator for displacing the atoms from the generated underline distribution
573 : class wiggleTrajectory: public AtomDistribution {
574 : std::unique_ptr<AtomDistribution> distribution;
575 : double radius;
576 : public:
577 : static constexpr auto id="wiggle";
578 : static constexpr auto doc=R"=(displaces all the atoms by a certain amount
579 : usage:
580 : wiggle (optional: sphere radius=0.1))=";
581 1 : static std::unique_ptr<AtomDistribution> decorate(
582 : std::unique_ptr<AtomDistribution>&& d,
583 : std::string_view cmd) {
584 : unpackedLine lines;
585 1 : Tools::getWordsSimple(lines,cmd);
586 1 : plumed_assert(lines.size() <=2) << id << " supports maximum one input";
587 1 : if (lines.size()==2) {
588 1 : double displacement=0.1;
589 1 : Tools::convert(std::string(lines[1]), displacement);
590 1 : return std::make_unique<wiggleTrajectory>(std::move(d),displacement);
591 : }
592 : //the default value is specified in the header
593 0 : return std::make_unique<wiggleTrajectory>(std::move(d));
594 : }
595 :
596 : wiggleTrajectory(std::unique_ptr<AtomDistribution>&& d,
597 1 : double amount=0.1) :
598 : distribution(std::move(d)),
599 1 : radius(amount)
600 : {}
601 :
602 1 : void frame(View<Vector> posToUpdate,
603 : View<double,9> box,
604 : unsigned step,
605 : Random& rng) override {
606 1 : distribution->frame(posToUpdate,box,step,rng);
607 :
608 1 : UniformSphericalVector usv(radius);
609 201 : for (auto& v: posToUpdate) {
610 200 : v+= usv(rng);
611 : }
612 1 : }
613 0 : bool overrideNat(unsigned& natoms) override {
614 0 : return distribution->overrideNat(natoms);
615 : }
616 : };
617 :
618 :
619 :
620 : /// A decorator for forcing a box in the trajectory
621 : class forceBoxTrajectory: public AtomDistribution {
622 : std::unique_ptr<AtomDistribution> distribution;
623 : std::array<double,9> fixedbox;
624 : public:
625 : static constexpr auto id="box";
626 : static constexpr auto doc=R"=(Forces a box on the atom distribution
627 : usage:
628 : box xx yy zz
629 : or
630 : box xx xy xy yx yy yz zx zy zz)=";
631 :
632 5 : static std::unique_ptr<AtomDistribution> decorate(
633 : std::unique_ptr<AtomDistribution>&& d,
634 : std::string_view cmd) {
635 : unpackedLine lines;
636 5 : Tools::getWordsSimple(lines,cmd);
637 5 : plumed_assert(lines.size() == 4 || lines.size() == 10) << id << " supports either 3 or 9 numbers as input";
638 5 : std::array<double,9> newBox= {0.0,0.0,0.0,
639 : 0.0,0.0,0.0,
640 : 0.0,0.0,0.0
641 : };
642 5 : if (lines.size()==4) {
643 6 : Tools::convert(std::string(lines[1]), newBox[0]);
644 6 : Tools::convert(std::string(lines[2]), newBox[4]);
645 6 : Tools::convert(std::string(lines[3]), newBox[8]);
646 : } else {// we now is 9
647 4 : Tools::convert(std::string(lines[1]), newBox[0]);
648 4 : Tools::convert(std::string(lines[2]), newBox[1]);
649 4 : Tools::convert(std::string(lines[3]), newBox[2]);
650 4 : Tools::convert(std::string(lines[4]), newBox[3]);
651 4 : Tools::convert(std::string(lines[5]), newBox[4]);
652 4 : Tools::convert(std::string(lines[6]), newBox[5]);
653 4 : Tools::convert(std::string(lines[7]), newBox[6]);
654 4 : Tools::convert(std::string(lines[8]), newBox[7]);
655 4 : Tools::convert(std::string(lines[9]), newBox[8]);
656 : }
657 : //the default value is specified in the header
658 5 : return std::make_unique<forceBoxTrajectory>(std::move(d),newBox);
659 : }
660 : forceBoxTrajectory(std::unique_ptr<AtomDistribution>&& d,
661 5 : const std::array<double,9>& mybox):
662 : distribution(std::move(d)),
663 5 : fixedbox(mybox)
664 : {}
665 :
666 7 : void frame(View<Vector> posToUpdate,
667 : View<double,9> box,
668 : unsigned step,
669 : Random& rng) override {
670 7 : distribution->frame(posToUpdate,box,step,rng);
671 7 : box[0] = fixedbox[0];
672 7 : box[1] = fixedbox[1];
673 7 : box[2] = fixedbox[2];
674 7 : box[3] = fixedbox[3];
675 7 : box[4] = fixedbox[4];
676 7 : box[5] = fixedbox[5];
677 7 : box[6] = fixedbox[6];
678 7 : box[7] = fixedbox[7];
679 7 : box[8] = fixedbox[8];
680 7 : }
681 :
682 1 : bool overrideNat(unsigned& natoms) override {
683 1 : return distribution->overrideNat(natoms);
684 : }
685 : };
686 :
687 : /// A decorator for fixing the state of the trajectory
688 : class fixedTrajectory: public AtomDistribution {
689 : std::unique_ptr<AtomDistribution> distribution;
690 : std::array<double,9> fixedBox= {0.0,0.0,0.0,
691 : 0.0,0.0,0.0,
692 : 0.0,0.0,0.0
693 : };
694 : std::vector<Vector> positions= {};
695 : unsigned stride=0;
696 : bool generate=true;
697 : public:
698 : static constexpr auto id="fix";
699 : static constexpr auto doc=R"=(Fixes the trajectory for stride steps, or forever if stride is 0
700 : usage:
701 : fix (optional stride=0))=";
702 :
703 5 : static std::unique_ptr<AtomDistribution> decorate(
704 : std::unique_ptr<AtomDistribution>&& d,
705 : std::string_view cmd) {
706 : unpackedLine lines;
707 5 : Tools::getWordsSimple(lines,cmd);
708 5 : plumed_assert(lines.size() <= 2) << id << " supports only one optional parameter for the stride";
709 :
710 5 : if (lines.size()==2) {
711 0 : unsigned newStride=0;
712 0 : Tools::convert(std::string(lines[1]), newStride);
713 0 : return std::make_unique<fixedTrajectory>(std::move(d),newStride);
714 : }
715 : //the default value is specified in the header
716 5 : return std::make_unique<fixedTrajectory>(std::move(d));
717 : }
718 : fixedTrajectory(std::unique_ptr<AtomDistribution>&& d,
719 5 : unsigned newStride=0):
720 : distribution(std::move(d)),
721 5 : stride(newStride)
722 : {}
723 :
724 10 : void frame(View<Vector> posToUpdate,
725 : View<double,9> box,
726 : unsigned step,
727 : Random& rng) override {
728 10 : generate = (!generate && stride>0) ?
729 0 : step%stride==0
730 : : generate;
731 10 : if (generate) {
732 5 : positions.resize(posToUpdate.size());
733 5 : distribution->frame(
734 : make_view(positions),
735 : View<double,9>(fixedBox.data()),
736 : step,rng);
737 5 : generate=false;
738 : }
739 10 : plumed_assert(positions.size() == posToUpdate.size()) << "fixed atom distributions: the expected number of atoms changed";
740 : std::copy(positions.begin(),positions.end(),posToUpdate.begin());
741 10 : box[0] = fixedBox[0];
742 10 : box[1] = fixedBox[1];
743 10 : box[2] = fixedBox[2];
744 10 : box[3] = fixedBox[3];
745 10 : box[4] = fixedBox[4];
746 10 : box[5] = fixedBox[5];
747 10 : box[6] = fixedBox[6];
748 10 : box[7] = fixedBox[7];
749 10 : box[8] = fixedBox[8];
750 10 : }
751 :
752 0 : bool overrideNat(unsigned& natoms) override {
753 0 : return distribution->overrideNat(natoms);
754 : }
755 : };
756 : //****************************Helpers and managements***************************
757 :
758 :
759 : namespace ADHelpers {
760 : //How to use this:
761 : //
762 : //The following variants act as "registers" for the distributions and the decorators
763 : //
764 : //To register a new distribution, siply write the struct in the file above
765 : //then add it to the rigth list
766 : //
767 : //Distributions and decorators must have a declared
768 : //`constexpr auto id="name"` and `constexpr auto doc="documentation"` members,
769 : //to be correcly listed and parsed
770 : //
771 : //To work a distribution simply needs to have a frame method, and no constructor.
772 : //
773 : //To add a complex distribution that needs a constructor you need to
774 : //do some more work like the file one, in its own source AtomDistributionFiles.cpp/h.
775 : //It is still compatible with the decorators.
776 : //
777 : //The decorators needs some more work:
778 : //a decorator needs to store the unique_ptr of the distribution that modifies
779 : //a decorator may need to store some extra settings
780 : //a decorator needs a static method "decorate"
781 : // that accepts a rvalue unique_ptr of the distribution that it will modify,
782 : // and a string of data to be parsed,
783 : // this method will return a unique_ptr with the modified distribution.
784 : // Note that the decorator will be the new owner of the original distribution pointer
785 : // and that the command string will contain the declared id of the decorator ("id settings")
786 : //a decorator will needd to override the overrideNat method, ideally only like this:
787 : //```
788 : //bool forceBoxTrajectory::overrideNat(unsigned& natoms) {
789 : // return distribution->overrideNat(natoms);
790 : //}
791 : //```
792 : // To propagate the changes in the number of atoms by all the decorators
793 : // applied to a distribution, it can be more complex (like in "reply")
794 : //a decorator must overload the frame method, and call the one of the
795 : // stored distribution and eventually modify the results of it
796 : //
797 : //Note that as time of writing this, there is no check implemented to verify
798 : // if two distributuin have the same id
799 :
800 : //add a new base distribution here, and the compiler will automagically
801 : //set up the parser and the documentation for you with some template shenanigans
802 : using baseDistributions = std::variant<
803 : theLine,
804 : uniformCube,
805 : uniformSphere,
806 : twoGlobs,
807 : tiledSimpleCubic,
808 : tiledBodyCenteredCubic,
809 : tiledFaceCenteredCubic,
810 : inscribedBodyCenteredCubic,
811 : inscribedFaceCenteredCubic>;
812 :
813 : //add new decorators here:
814 : using decoratorDistribuitions = std::variant<
815 : repliedTrajectory,
816 : scaledTrajectory,
817 : wiggleTrajectory,
818 : forceBoxTrajectory,
819 : fixedTrajectory>;
820 :
821 : template <size_t I=0>
822 589 : std::optional<std::unique_ptr<AtomDistribution>> getAD(std::string_view atomicDistr) {
823 : if constexpr ( I < std::variant_size_v<baseDistributions>) {
824 : using myAD=typename std::variant_alternative_t<I,baseDistributions>;
825 589 : if (atomicDistr==myAD::id) {
826 124 : return std::make_unique<myAD>();
827 : }
828 465 : return getAD<I+1>(atomicDistr);
829 : }
830 0 : return std::nullopt;
831 : }
832 :
833 : using parser=std::function<std::unique_ptr<AtomDistribution>(std::unique_ptr<AtomDistribution>&& d,
834 : std::string_view cmd)>;
835 :
836 : template <size_t I=0>
837 88 : std::optional<parser> getDecorator( std::string_view decoratorDistr) {
838 : if constexpr ( I < std::variant_size_v<decoratorDistribuitions>) {
839 88 : auto name=decoratorDistr.substr(0,decoratorDistr.find(' '));
840 : using myAD=typename std::variant_alternative_t<I,decoratorDistribuitions>;
841 88 : if (name==myAD::id) {
842 : return myAD::decorate;
843 : }
844 49 : return getDecorator<I+1>(decoratorDistr);
845 : }
846 0 : return std::nullopt;
847 : }
848 :
849 : template <size_t I=0>
850 0 : void getDistributionList(std::vector<std::string>& list) {
851 : if constexpr ( I < std::variant_size_v<baseDistributions>) {
852 : using myAD=typename std::variant_alternative_t<I,baseDistributions>;
853 0 : list.push_back(myAD::id);
854 0 : getDistributionList<I+1>(list);
855 : }
856 0 : }
857 :
858 : template <size_t I=0>
859 9 : void getDistributionDocumentation( std::vector<AtomDistribution::documentation>& list) {
860 : if constexpr ( I < std::variant_size_v<baseDistributions>) {
861 : using myAD=typename std::variant_alternative_t<I,baseDistributions>;
862 9 : list.push_back({myAD::id,myAD::doc});
863 8 : getDistributionDocumentation<I+1>(list);
864 : }
865 18 : }
866 :
867 : template <size_t I=0>
868 0 : void getDecoratorList(std::vector<std::string>& list) {
869 : if constexpr ( I < std::variant_size_v<decoratorDistribuitions>) {
870 : using myAD=typename std::variant_alternative_t<I,decoratorDistribuitions>;
871 0 : list.push_back(myAD::id);
872 0 : getDecoratorList<I+1>(list);
873 : }
874 0 : }
875 :
876 : template <size_t I=0>
877 5 : void getDecoratorDocumentation(std::vector<AtomDistribution::documentation>& list) {
878 : if constexpr ( I < std::variant_size_v<decoratorDistribuitions>) {
879 : using myAD=typename std::variant_alternative_t<I,decoratorDistribuitions>;
880 5 : list.push_back({myAD::id,myAD::doc});
881 4 : getDecoratorDocumentation<I+1>(list);
882 : }
883 10 : }
884 :
885 0 : std::string getDistributionList_forError() {
886 0 : auto list = AtomDistribution::getDistributionList();
887 0 : std::string toRet="";
888 0 : std::string pre="\"";
889 0 : for (const auto & distr:list) {
890 0 : toRet += pre+distr+ "\"";
891 : pre=", \"";
892 : }
893 0 : return toRet;
894 0 : }
895 :
896 0 : std::string getDecoratorList_forError() {
897 0 : auto list = AtomDistribution::getDecoratorsList();
898 0 : std::string toRet="";
899 0 : std::string pre="\"";
900 0 : for (const auto& distr:list) {
901 0 : toRet += pre+distr+ "\"";
902 : pre=", \"";
903 : }
904 0 : return toRet;
905 0 : }
906 : } // namespace ADHelpers
907 :
908 124 : std::unique_ptr<AtomDistribution> AtomDistribution::getAtomDistribution(std::string_view atomicDistr) {
909 124 : auto pipePos = atomicDistr.find("|");
910 124 : auto base=atomicDistr.substr(0,pipePos);
911 124 : std::unique_ptr<AtomDistribution> distribution;
912 124 : if(auto d = ADHelpers::getAD(base); d) {
913 : distribution = std::move(*d);
914 : } else {
915 0 : plumed_error() << "The atomic distribution can be only "<<ADHelpers::getDistributionList_forError()
916 : << ", but the input was \""
917 0 : << base <<'"';
918 : }
919 124 : if (pipePos != std::string_view::npos) {
920 : return decorateAtomDistribution(std::move(distribution),
921 12 : atomicDistr.substr(pipePos+1));
922 : }
923 : return distribution;
924 124 : }
925 :
926 37 : std::unique_ptr<AtomDistribution> AtomDistribution::decorateAtomDistribution(
927 : std::unique_ptr<AtomDistribution> && ad,
928 : std::string_view decoratorDistr) {
929 : unpackedLine lines;
930 37 : Tools::getWordsSimple(lines,decoratorDistr,"|");
931 : std::unique_ptr<AtomDistribution> distribution=std::move(ad);
932 76 : for(auto i=0u; i < lines.size(); ++i) {
933 39 : auto decorator=ADHelpers::getDecorator(lines[i]);
934 39 : if (decorator) {
935 : auto tmp=std::move(distribution);
936 78 : distribution = (*decorator)(std::move(tmp),lines[i]);
937 39 : } else {
938 :
939 0 : plumed_error() << "The atomic distribution decorators be only "<<ADHelpers::getDecoratorList_forError()
940 : << ", but the input was \""
941 0 : << lines[i] <<'"';
942 : }
943 :
944 : }
945 37 : return distribution;
946 0 : }
947 :
948 0 : std::vector<std::string> AtomDistribution::getDistributionList() {
949 0 : std::vector<std::string> list;
950 0 : ADHelpers::getDistributionList(list);
951 0 : return list;
952 0 : }
953 1 : std::vector<AtomDistribution::documentation> AtomDistribution::getDistributionDocumentation() {
954 1 : std::vector<documentation> list;
955 1 : ADHelpers::getDistributionDocumentation(list);
956 1 : return list;
957 0 : }
958 :
959 0 : std::vector<std::string> AtomDistribution::getDecoratorsList() {
960 0 : std::vector<std::string> list;
961 0 : ADHelpers::getDecoratorList(list);
962 0 : return list;
963 0 : }
964 :
965 1 : std::vector<AtomDistribution::documentation> AtomDistribution::getDecoratorsDocumentation() {
966 1 : std::vector<documentation> list;
967 1 : ADHelpers::getDecoratorDocumentation(list);
968 1 : return list;
969 0 : }
970 :
971 : } //namespace PLMD
|
