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
@@ -128,11 +155,7 @@ void scaleTrajCheck(std::string_view kind,
128155 std::vector<double > basebox (9 );
129156 auto scaled= [&]() {
130157 std::unique_ptr<PLMD ::AtomDistribution> d;
131- if (kind == " sphere-reply212" ) {
132- d = AtomDistribution::getAtomDistribution (" sphere|reply 2 1 2" );
133- } else {
134158 d = AtomDistribution::getAtomDistribution (kind);
135- }
136159 d->frame (baseatoms,basebox,0 ,rng);
137160
138161 auto mod=" scale " + std::to_string (scale) + " " ;
@@ -169,6 +192,81 @@ void scaleTrajCheck(std::string_view kind,
169192 }
170193}
171194
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+ }
172270int main () {
173271 std::ofstream ofs (" output" );
174272 ofs << std::boolalpha;
@@ -177,7 +275,11 @@ int main() {
177275 " cube" ,
178276 " sphere" ,
179277 " globs" ,
180- " sc"
278+ " sc" ,
279+ " fcc" ,
280+ " bcc" ,
281+ " ifcc" ,
282+ " ibcc"
181283 }) {
182284 basecheck (kind,ofs);
183285 ofs << " \n " ;
@@ -196,8 +298,17 @@ int main() {
196298 scaleTrajCheck (" sphere" ,num,ofs);
197299 scaleTrajCheck (" sc" ,num,ofs);
198300 scaleTrajCheck (" globs" ,num,ofs);
199- scaleTrajCheck (" sphere-reply212 " ,num,ofs);
301+ scaleTrajCheck (" sphere|reply 2 1 2 " ,num,ofs);
200302 ofs << " \n " ;
201303 }
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);
202313 return 0 ;
203314}
0 commit comments