testUnitCell.C 4.18 KB
Newer Older
Francois Gygi committed
1 2
////////////////////////////////////////////////////////////////////////////////
//
Francois Gygi committed
3 4 5 6
// Copyright (c) 2008 The Regents of the University of California
//
// This file is part of Qbox
//
Francois Gygi committed
7 8
// Qbox is distributed under the terms of the GNU General Public License
// as published by the Free Software Foundation, either version 2 of
Francois Gygi committed
9 10 11 12 13 14
// the License, or (at your option) any later version.
// See the file COPYING in the root directory of this distribution
// or <http://www.gnu.org/licenses/>.
//
////////////////////////////////////////////////////////////////////////////////
//
Francois Gygi committed
15 16 17
// testUnitCell.C
//
////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
18
// $Id: testUnitCell.C,v 1.5 2008-09-08 15:56:20 fgygi Exp $
Francois Gygi committed
19 20 21 22 23 24 25 26 27 28

#include "UnitCell.h"
#include <iostream>
#include <iomanip>
#include <cstdlib>
using namespace std;

int main()
{
  UnitCell cell;
29

Francois Gygi committed
30
  cout << " orthorhombic cell: " << endl;
31
  cell.set(D3vector(4.,0.,0.),D3vector(0.,5.,0.),D3vector(0.,0.,7.));
Francois Gygi committed
32
  cout << cell << endl;
33

Francois Gygi committed
34
  D3vector v;
35

Francois Gygi committed
36 37 38
  const int n = 100000;
  const double d = 10.0;
  int count;
39

Francois Gygi committed
40 41 42 43 44 45 46 47
  count = 0;
  for ( int i = 0; i < n; i++ )
  {
    // generate a random vector in a box of size d x d x d
    double x0 = drand48() - 0.5;
    double x1 = drand48() - 0.5;
    double x2 = drand48() - 0.5;
    v = d * D3vector(x0,x1,x2);
48

Francois Gygi committed
49 50 51 52 53 54 55 56
    if ( cell.in_ws(v) )
      count++;
    cell.fold_in_ws(v);
    if ( !cell.in_ws(v) )
      cout << v << endl;
  }
  cout << " volume ratio: " << cell.volume() / (d*d*d)
       << "  computed: " << count/((double) n) << endl << endl;;
57

Francois Gygi committed
58
  cout << " FCC cell: " << endl;
59
  cell.set(D3vector(4,4,0),D3vector(0,4,4),D3vector(4,0,4));
Francois Gygi committed
60
  cout << cell << endl;
61

Francois Gygi committed
62 63 64 65 66 67 68 69
  count = 0;
  for ( int i = 0; i < n; i++ )
  {
    // generate a random vector in a box of size d x d x d
    double x0 = drand48() - 0.5;
    double x1 = drand48() - 0.5;
    double x2 = drand48() - 0.5;
    v = d * D3vector(x0,x1,x2);
70

Francois Gygi committed
71 72 73 74 75 76 77 78
    if ( cell.in_ws(v) )
      count++;
    cell.fold_in_ws(v);
    if ( !cell.in_ws(v) )
      cout << v << endl;
  }
  cout << " volume ratio: " << cell.volume() / (d*d*d)
       << "  computed: " << count/((double) n) << endl << endl;;
79

Francois Gygi committed
80
  cout << " BCC cell: " << endl;
81
  cell.set(D3vector(4,4,-4),D3vector(-4, 4,4),D3vector(4,-4,4));
Francois Gygi committed
82
  cout << cell << endl;
83

Francois Gygi committed
84 85 86 87 88 89 90 91
  count = 0;
  for ( int i = 0; i < n; i++ )
  {
    // generate a random vector in a box of size d x d x d
    double x0 = drand48() - 0.5;
    double x1 = drand48() - 0.5;
    double x2 = drand48() - 0.5;
    v = d * D3vector(x0,x1,x2);
92

Francois Gygi committed
93 94 95 96 97 98 99 100
    if ( cell.in_ws(v) )
      count++;
    cell.fold_in_ws(v);
    if ( !cell.in_ws(v) )
      cout << v << endl;
  }
  cout << " volume ratio: " << cell.volume() / (d*d*d)
       << "  computed: " << count/((double) n) << endl << endl;;
101

Francois Gygi committed
102
  cout << " Monoclinic cell: " << endl;
103
  cell.set(D3vector(4,0,0.1),D3vector(0.2, 4,0),D3vector(0.1,0.3,4));
Francois Gygi committed
104
  cout << cell << endl;
105

106 107 108 109 110
  // Check matrix multiplication function
  double ztest[9];
  cell.matmult3x3(cell.amat(),cell.amat_inv(),ztest);
  for ( int i = 0; i < 9; i++ )
    cout << " ztest[" << i << "]=" << ztest[i] << endl;
111

Francois Gygi committed
112 113 114 115 116 117 118 119
  count = 0;
  for ( int i = 0; i < n; i++ )
  {
    // generate a random vector in a box of size d x d x d
    double x0 = drand48() - 0.5;
    double x1 = drand48() - 0.5;
    double x2 = drand48() - 0.5;
    v = d * D3vector(x0,x1,x2);
120

Francois Gygi committed
121 122 123 124 125 126 127 128
    if ( cell.in_ws(v) )
      count++;
    cell.fold_in_ws(v);
    if ( !cell.in_ws(v) )
      cout << v << endl;
  }
  cout << " volume ratio: " << cell.volume() / (d*d*d)
       << "  computed: " << count/((double) n) << endl << endl;
129

Francois Gygi committed
130 131
  // test UnitCell::included_in function
  UnitCell rcell;
132
  rcell.set(D3vector(4,4,0),D3vector(0,4,4),D3vector(4,0,4));
Francois Gygi committed
133 134 135 136 137 138
  cell.set(D3vector(4,5,0),D3vector(0,4,4),D3vector(4,0,4));
  cout << " rcell: " << rcell << endl;
  cout << " cell: " << cell << endl;
  cout << " cell is";
  if ( !rcell.encloses(cell) ) cout << " not";
  cout << " included in rcell" << endl;
139 140

  rcell.set(D3vector(4,4,0),D3vector(0,4,4),D3vector(4,0,4));
Francois Gygi committed
141 142 143 144 145
  cout << " rcell: " << rcell << endl;
  cell.set(D3vector(3.8,3.8,0),D3vector(0,3.8,3.8),D3vector(3.8,0,3.8));
  cout << " cell: " << cell << endl;
  cout << " cell is";
  if ( !rcell.encloses(cell) ) cout << " not";
146
  cout << " included in rcell" << endl;
Francois Gygi committed
147
}