11#include " plumed/tools/AtomDistribution.h"
22#include " plumed/tools/Random.h"
33#include " plumed/tools/Vector.h"
4+ #include " plumed/tools/Tools.h"
5+ #include " plumed/tools/Pbc.h"
46#include < array>
57#include < fstream>
68#include < memory>
@@ -15,24 +17,49 @@ void atomsInBoxCheck(
1517 const std::vector<double > &box,
1618 const std::string_view header,
1719 std::ostream & ofs) {
18- Vector lowbound = atoms[0 ];
19- Vector upbound= atoms[0 ];
20- for (unsigned i =1 ; i < atoms.size (); ++i) {
21- for (unsigned j=0 ; j<3 ; ++j) {
22- if (atoms[i][j] < lowbound[j]) {
23- lowbound[j] = atoms[i][j];
20+ Pbc mybox;
21+ mybox.setBox (Tensor{
22+ box[0 ],
23+ box[1 ],
24+ box[2 ],
25+ box[3 ],
26+ box[4 ],
27+ box[5 ],
28+ box[6 ],
29+ box[7 ],
30+ box[8 ]});
31+ ofs << " Atoms are within box dimensions:\n " ;
32+ if (!mybox.isOrthorombic ()) {
33+ bool inbox = true ;
34+ for (const auto & atom : atoms) {
35+ auto scaled = mybox.realToScaled (atom);
36+ inbox &= scaled[0 ]<1.0 ;
37+ inbox &= scaled[1 ]<1.0 ;
38+ inbox &= scaled[2 ]<1.0 ;
39+ if (!inbox) {
40+ break ;
2441 }
25- if (atoms[i][j] > upbound[j]) {
26- upbound[j] = atoms[i][j];
42+ }
43+ ofs <<header <<" all atoms within the non orthorombic box: " << inbox << " \n " ;
44+ } else {
45+ Vector lowbound = atoms[0 ];
46+ Vector upbound= atoms[0 ];
47+ for (unsigned i =1 ; i < atoms.size (); ++i) {
48+ for (unsigned j=0 ; j<3 ; ++j) {
49+ if (atoms[i][j] < lowbound[j]) {
50+ lowbound[j] = atoms[i][j];
51+ }
52+ if (atoms[i][j] > upbound[j]) {
53+ upbound[j] = atoms[i][j];
54+ }
2755 }
2856 }
29- }
30- // box is orhtorombic an starts in 0,0,0:
31- // shifting lowbound and upbound to chek the box:
32- upbound-=lowbound;
33- ofs << " Atoms are within box dimensions:\n " ;
34- for (unsigned j=0 ; j<3 ; ++j) {
35- ofs <<header <<" all atoms in box along " <<xyz[j]<<" : " << (upbound[j]<box[j*3 +j]) << " \n " ;
57+ // box is orhtorombic an starts in 0,0,0:
58+ // shifting lowbound and upbound to chek the box:
59+ upbound-=lowbound;
60+ for (unsigned j=0 ; j<3 ; ++j) {
61+ ofs <<header <<" all atoms in box along " <<xyz[j]<<" : " << (upbound[j]<box[j*3 +j]) << " \n " ;
62+ }
3663 }
3764}
3865
@@ -63,14 +90,15 @@ void replyTrajCheck(std::string_view kind,
6390 unsigned nat = atoms.size ();
6491 const auto oldNat=nat;
6592
66- std::unique_ptr< PLMD ::AtomDistribution> rep= std::make_unique<repliedTrajectory>( [&]() {
93+ auto rep= [&]() {
6794 auto d = AtomDistribution::getAtomDistribution (kind);
6895 d->frame (atoms,basebox,0 ,rng);
69- return d;
96+ auto mod=" reply " + std::to_string (num[0 ]) + " "
97+ + std::to_string (num[1 ]) + " "
98+ + std::to_string (num[2 ]);
99+ return AtomDistribution::decorateAtomDistribution (std::move (d),mod);
70100 }
71- (),
72- num[0 ], num[1 ], num[2 ],
73- nat);
101+ ();
74102
75103 // this must return true
76104 ofs <<header<< " rep->overrideNat(nat)=" <<
@@ -125,20 +153,16 @@ void scaleTrajCheck(std::string_view kind,
125153 rng.setSeed (12345 );
126154 std::vector<Vector> baseatoms (200 );
127155 std::vector<double > basebox (9 );
128- std::unique_ptr< PLMD ::AtomDistribution> scaled= std::make_unique<scaledTrajectory>( [&]() {
156+ auto scaled= [&]() {
129157 std::unique_ptr<PLMD ::AtomDistribution> d;
130- if (kind == " sphere-reply212" ) {
131- d = std::make_unique<repliedTrajectory>(
132- AtomDistribution::getAtomDistribution (" sphere" ),
133- 2 ,1 ,2 ,baseatoms.size ()/(2 *1 *2 ));
134- } else {
135158 d = AtomDistribution::getAtomDistribution (kind);
136- }
137159 d->frame (baseatoms,basebox,0 ,rng);
138- return d;
160+
161+ auto mod=" scale " + std::to_string (scale) + " " ;
162+ return
163+ AtomDistribution::decorateAtomDistribution (std::move (d),mod);
139164 }
140- (),
141- scale);
165+ ();
142166
143167 std::vector<Vector> atoms (200 );
144168 std::vector<double > box (9 );
@@ -151,7 +175,8 @@ void scaleTrajCheck(std::string_view kind,
151175 ofs << header << " The atoms are scaled correctly:\t " ;
152176 bool correct = true ;
153177 for (unsigned i =0 ; i< atoms.size () && correct; ++i) {
154- correct = (abs (scale*baseatoms[i][0 ] - atoms[i][0 ]) < 1000 *PLMD ::epsilon)&&
178+ correct = (abs (scale*baseatoms[i][0 ] - atoms[i][0 ]) < 1000 *PLMD ::epsilon)
179+ &&
155180 (abs (scale*baseatoms[i][1 ] - atoms[i][1 ]) < 1000 *PLMD ::epsilon)&&
156181 (abs (scale*baseatoms[i][2 ] - atoms[i][2 ]) < 1000 *PLMD ::epsilon);
157182
@@ -161,11 +186,87 @@ void scaleTrajCheck(std::string_view kind,
161186 atomsInBoxCheck (atoms,box,header,ofs);
162187 ofs << " New box has the correct dimensions:\n " ;
163188 for (unsigned j=0 ; j<3 ; ++j) {
164- ofs <<header <<" correct box dimension " <<xyz[j]<<" : "
189+ ofs <<header <<" correct box dimension " <<xyz[j]
190+ <<" : "
165191 << ((scale*basebox[j*3 +j] - box[j*3 +j]) < 1000 *PLMD ::epsilon) << " \n " ;
166192 }
167193}
168194
195+ void fixTrajCheck (std::string_view kind,
196+ std::ostream & ofs) {
197+ std::stringstream ss;
198+ ss << " [fixTrajCheck -" << kind << " -]:" ;
199+ auto header = ss.str ();
200+ ofs << header << " \n " ;
201+ // reinitialized each time for stability
202+ Random rng;
203+ rng.setSeed (12345 );
204+ std::vector<Vector> baseatoms (200 );
205+ std::vector<double > basebox (9 );
206+ // sphere generates a new configuration at each step
207+ auto d = AtomDistribution::getAtomDistribution (std::string (kind)+" |fix" );
208+ d->frame (baseatoms,basebox,0 ,rng);
209+ std::vector<Vector> atoms (200 );
210+ std::vector<double > box (9 );
211+ // generating the next frame
212+ d->frame (atoms,box,0 ,rng);
213+
214+ ofs << header << " The atoms are fixed correctly:\t " ;
215+ bool correct = true ;
216+ for (unsigned i =0 ; i< atoms.size () && correct; ++i) {
217+ correct = (abs (baseatoms[i][0 ] - atoms[i][0 ]) < 1000 *PLMD ::epsilon)
218+ &&
219+ (abs (baseatoms[i][1 ] - atoms[i][1 ]) < 1000 *PLMD ::epsilon)&&
220+ (abs (baseatoms[i][2 ] - atoms[i][2 ]) < 1000 *PLMD ::epsilon);
221+
222+ }
223+ ofs << correct << " \n " ;
224+ }
225+
226+ void forceBoxCheck (std::string_view kind, std::ostream & ofs) {
227+ std::string header= " [forceBoxCheck -" + std::string (kind) + " -]:" ;
228+ ofs << header << " \n " ;
229+ {
230+ auto d = AtomDistribution::getAtomDistribution (std::string (kind)+" |box 1 2 3" );
231+ // reinitialized each time for stability
232+ Random rng;
233+ std::vector<Vector> atoms (200 );
234+ std::vector<double > box (9 );
235+ d->frame (atoms,box,0 ,rng);
236+ ofs << header << " The box is changed as asked (ortho): " ;
237+
238+ bool success = (box[0 ] - 1.0 )< 1000 * PLMD ::epsilon &&
239+ box[1 ] < 1000 * PLMD ::epsilon &&
240+ box[2 ] < 1000 * PLMD ::epsilon &&
241+ box[3 ] < 1000 * PLMD ::epsilon &&
242+ (box[4 ] - 2.0 ) < 1000 * PLMD ::epsilon &&
243+ box[5 ] < 1000 * PLMD ::epsilon &&
244+ box[6 ] < 1000 * PLMD ::epsilon &&
245+ box[7 ] < 1000 * PLMD ::epsilon &&
246+ (box[8 ] - 3.0 ) < 1000 * PLMD ::epsilon;
247+ ofs << success << " \n " ;
248+ }
249+ {
250+ auto d = AtomDistribution::getAtomDistribution (std::string (kind)+" |box 1 2 3 4 5 6 7 8 9" );
251+ // reinitialized each time for stability
252+ Random rng;
253+ std::vector<Vector> atoms (200 );
254+ std::vector<double > box (9 );
255+ d->frame (atoms,box,0 ,rng);
256+ ofs << header << " The box is changed as asked (9 elements): " ;
257+
258+ bool success = (box[0 ] - 1.0 )< 1000 * PLMD ::epsilon &&
259+ (box[1 ] - 2.0 ) < 1000 * PLMD ::epsilon &&
260+ (box[2 ] - 3.0 ) < 1000 * PLMD ::epsilon &&
261+ (box[3 ] - 4.0 ) < 1000 * PLMD ::epsilon &&
262+ (box[4 ] - 5.0 ) < 1000 * PLMD ::epsilon &&
263+ (box[5 ] - 6.0 ) < 1000 * PLMD ::epsilon &&
264+ (box[6 ] - 7.0 ) < 1000 * PLMD ::epsilon &&
265+ (box[7 ] - 8.0 ) < 1000 * PLMD ::epsilon &&
266+ (box[8 ] - 9.0 ) < 1000 * PLMD ::epsilon;
267+ ofs << success << " \n " ;
268+ }
269+ }
169270int main () {
170271 std::ofstream ofs (" output" );
171272 ofs << std::boolalpha;
@@ -174,7 +275,11 @@ int main() {
174275 " cube" ,
175276 " sphere" ,
176277 " globs" ,
177- " sc"
278+ " sc" ,
279+ " fcc" ,
280+ " bcc" ,
281+ " ifcc" ,
282+ " ibcc"
178283 }) {
179284 basecheck (kind,ofs);
180285 ofs << " \n " ;
@@ -193,8 +298,17 @@ int main() {
193298 scaleTrajCheck (" sphere" ,num,ofs);
194299 scaleTrajCheck (" sc" ,num,ofs);
195300 scaleTrajCheck (" globs" ,num,ofs);
196- scaleTrajCheck (" sphere-reply212 " ,num,ofs);
301+ scaleTrajCheck (" sphere|reply 2 1 2 " ,num,ofs);
197302 ofs << " \n " ;
198303 }
304+ forceBoxCheck (" sc" ,ofs);
305+ forceBoxCheck (" fcc" ,ofs);
306+ // I actually do not know how to test wiggle in a sensible way
307+ // cube, globs and sphere are generate a randm configuration at each step
308+ fixTrajCheck (" cube" ,ofs);
309+ fixTrajCheck (" globs" ,ofs);
310+ fixTrajCheck (" sphere" ,ofs);
311+ fixTrajCheck (" sphere|reply 1 1 2" ,ofs);
312+ fixTrajCheck (" ifcc|wiggle 0.5" ,ofs);
199313 return 0 ;
200314}
0 commit comments