1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (c) 2008 The Regents of the University of California
4 //
5 // This file is part of Qbox
6 //
7 // Qbox is distributed under the terms of the GNU General Public License
8 // as published by the Free Software Foundation, either version 2 of
9 // the License, or (at your option) any later version.
10 // See the file COPYING in the root directory of this distribution
11 // or <http://www.gnu.org/licenses/>.
12 //
13 ////////////////////////////////////////////////////////////////////////////////
14 //
15 // testUnitCell.cpp
16 //
17 ////////////////////////////////////////////////////////////////////////////////
18
19 #include "UnitCell.h"
20 #include <iostream>
21 #include <iomanip>
22 #include <cstdlib>
23 using namespace std;
24
main()25 int main()
26 {
27 UnitCell cell;
28
29 cout << " orthorhombic cell: " << endl;
30 cell.set(D3vector(4.,0.,0.),D3vector(0.,5.,0.),D3vector(0.,0.,7.));
31 cout << cell << endl;
32
33 D3vector v;
34
35 const int n = 100000;
36 const double d = 10.0;
37 int count;
38
39 count = 0;
40 for ( int i = 0; i < n; i++ )
41 {
42 // generate a random vector in a box of size d x d x d
43 double x0 = drand48() - 0.5;
44 double x1 = drand48() - 0.5;
45 double x2 = drand48() - 0.5;
46 v = d * D3vector(x0,x1,x2);
47
48 if ( cell.in_ws(v) )
49 count++;
50 cell.fold_in_ws(v);
51 if ( !cell.in_ws(v) )
52 cout << v << endl;
53 }
54 cout << " volume ratio: " << cell.volume() / (d*d*d)
55 << " computed: " << count/((double) n) << endl << endl;;
56
57 cout << " FCC cell: " << endl;
58 cell.set(D3vector(4,4,0),D3vector(0,4,4),D3vector(4,0,4));
59 cout << cell << endl;
60
61 count = 0;
62 for ( int i = 0; i < n; i++ )
63 {
64 // generate a random vector in a box of size d x d x d
65 double x0 = drand48() - 0.5;
66 double x1 = drand48() - 0.5;
67 double x2 = drand48() - 0.5;
68 v = d * D3vector(x0,x1,x2);
69
70 if ( cell.in_ws(v) )
71 count++;
72 cell.fold_in_ws(v);
73 if ( !cell.in_ws(v) )
74 cout << v << endl;
75 }
76 cout << " volume ratio: " << cell.volume() / (d*d*d)
77 << " computed: " << count/((double) n) << endl << endl;;
78
79 cout << " BCC cell: " << endl;
80 cell.set(D3vector(4,4,-4),D3vector(-4, 4,4),D3vector(4,-4,4));
81 cout << cell << endl;
82
83 count = 0;
84 for ( int i = 0; i < n; i++ )
85 {
86 // generate a random vector in a box of size d x d x d
87 double x0 = drand48() - 0.5;
88 double x1 = drand48() - 0.5;
89 double x2 = drand48() - 0.5;
90 v = d * D3vector(x0,x1,x2);
91
92 if ( cell.in_ws(v) )
93 count++;
94 cell.fold_in_ws(v);
95 if ( !cell.in_ws(v) )
96 cout << v << endl;
97 }
98 cout << " volume ratio: " << cell.volume() / (d*d*d)
99 << " computed: " << count/((double) n) << endl << endl;;
100
101 cout << " Monoclinic cell: " << endl;
102 cell.set(D3vector(4,0,0.1),D3vector(0.2, 4,0),D3vector(0.1,0.3,4));
103 cout << cell << endl;
104
105 // Check matrix multiplication function
106 double ztest[9];
107 cell.matmult3x3(cell.amat(),cell.amat_inv(),ztest);
108 for ( int i = 0; i < 9; i++ )
109 cout << " ztest[" << i << "]=" << ztest[i] << endl;
110
111 count = 0;
112 for ( int i = 0; i < n; i++ )
113 {
114 // generate a random vector in a box of size d x d x d
115 double x0 = drand48() - 0.5;
116 double x1 = drand48() - 0.5;
117 double x2 = drand48() - 0.5;
118 v = d * D3vector(x0,x1,x2);
119
120 if ( cell.in_ws(v) )
121 count++;
122 cell.fold_in_ws(v);
123 if ( !cell.in_ws(v) )
124 cout << v << endl;
125 }
126 cout << " volume ratio: " << cell.volume() / (d*d*d)
127 << " computed: " << count/((double) n) << endl << endl;
128
129 // test UnitCell::included_in function
130 UnitCell rcell;
131 rcell.set(D3vector(4,4,0),D3vector(0,4,4),D3vector(4,0,4));
132 cell.set(D3vector(4,5,0),D3vector(0,4,4),D3vector(4,0,4));
133 cout << " rcell: " << rcell << endl;
134 cout << " cell: " << cell << endl;
135 cout << " cell is";
136 if ( !rcell.encloses(cell) ) cout << " not";
137 cout << " included in rcell" << endl;
138
139 rcell.set(D3vector(4,4,0),D3vector(0,4,4),D3vector(4,0,4));
140 cout << " rcell: " << rcell << endl;
141 cell.set(D3vector(3.8,3.8,0),D3vector(0,3.8,3.8),D3vector(3.8,0,3.8));
142 cout << " cell: " << cell << endl;
143 cout << " cell is";
144 if ( !rcell.encloses(cell) ) cout << " not";
145 cout << " included in rcell" << endl;
146 }
147