UnitCell.C 10.5 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/>.
//
////////////////////////////////////////////////////////////////////////////////
//
15
// UnitCell.C
Francois Gygi committed
16 17
//
////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
18
// $Id: UnitCell.C,v 1.19 2008-09-15 14:59:58 fgygi Exp $
Francois Gygi committed
19 20 21 22 23 24 25 26 27 28

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

////////////////////////////////////////////////////////////////////////////////
void UnitCell::set(const D3vector& a0, const D3vector& a1, const D3vector& a2)
{
  a_[0] = a0; a_[1] = a1, a_[2] = a2;
29 30 31 32 33 34 35 36 37
  amat_[0] = a0.x;
  amat_[1] = a0.y;
  amat_[2] = a0.z;
  amat_[3] = a1.x;
  amat_[4] = a1.y;
  amat_[5] = a1.z;
  amat_[6] = a2.x;
  amat_[7] = a2.y;
  amat_[8] = a2.z;
38

39
  // volume = det(A)
40
  volume_ = fabs(a0 * ( a1 ^ a2 ));
Francois Gygi committed
41 42
  if ( volume_ > 0.0 )
  {
43 44 45 46 47
    // Compute rows of A-1 (columns of A^-T)
    double fac = 1.0 / volume_;
    D3vector amt0 = fac * a1 ^ a2;
    D3vector amt1 = fac * a2 ^ a0;
    D3vector amt2 = fac * a0 ^ a1;
48

49 50 51 52 53 54 55 56 57
    amat_inv_[0] = amt0.x;
    amat_inv_[1] = amt1.x;
    amat_inv_[2] = amt2.x;
    amat_inv_[3] = amt0.y;
    amat_inv_[4] = amt1.y;
    amat_inv_[5] = amt2.y;
    amat_inv_[6] = amt0.z;
    amat_inv_[7] = amt1.z;
    amat_inv_[8] = amt2.z;
58

59 60 61 62 63 64 65 66 67
    amat_inv_t_[0] = amt0.x;
    amat_inv_t_[1] = amt0.y;
    amat_inv_t_[2] = amt0.z;
    amat_inv_t_[3] = amt1.x;
    amat_inv_t_[4] = amt1.y;
    amat_inv_t_[5] = amt1.z;
    amat_inv_t_[6] = amt2.x;
    amat_inv_t_[7] = amt2.y;
    amat_inv_t_[8] = amt2.z;
68

69 70 71 72 73 74 75 76 77 78 79 80 81 82
    // B = 2 pi A^-T
    b_[0] = 2.0 * M_PI * amt0;
    b_[1] = 2.0 * M_PI * amt1;
    b_[2] = 2.0 * M_PI * amt2;

    bmat_[0] = b_[0].x;
    bmat_[1] = b_[0].y;
    bmat_[2] = b_[0].z;
    bmat_[3] = b_[1].x;
    bmat_[4] = b_[1].y;
    bmat_[5] = b_[1].z;
    bmat_[6] = b_[2].x;
    bmat_[7] = b_[2].y;
    bmat_[8] = b_[2].z;
Francois Gygi committed
83 84 85 86
  }
  else
  {
    b_[0] = b_[1] = b_[2] = D3vector(0.0,0.0,0.0);
87 88
    amat_inv_[0] =  amat_inv_[1] =  amat_inv_[2] =
      amat_inv_[3] =  amat_inv_[4] =  amat_inv_[5] =
89
      amat_inv_[6] =  amat_inv_[7] =  amat_inv_[8] = 0.0;
90 91
    bmat_[0] =  bmat_[1] =  bmat_[2] =
      bmat_[3] =  bmat_[4] =  bmat_[5] =
92
      bmat_[6] =  bmat_[7] =  bmat_[8] = 0.0;
Francois Gygi committed
93
  }
94 95 96



Francois Gygi committed
97 98 99 100 101 102 103 104 105 106 107 108 109
  an_[0]  = a_[0];
  an_[1]  = a_[1];
  an_[2]  = a_[2];
  an_[3]  = a_[0] + a_[1];
  an_[4]  = a_[0] - a_[1];
  an_[5]  = a_[1] + a_[2];
  an_[6]  = a_[1] - a_[2];
  an_[7]  = a_[2] + a_[0];
  an_[8]  = a_[2] - a_[0];
  an_[9]  = a_[0] + a_[1] + a_[2];
  an_[10] = a_[0] - a_[1] - a_[2];
  an_[11] = a_[0] + a_[1] - a_[2];
  an_[12] = a_[0] - a_[1] + a_[2];
110

Francois Gygi committed
111 112 113 114 115 116 117 118 119 120 121 122 123
  bn_[0]  = b_[0];
  bn_[1]  = b_[1];
  bn_[2]  = b_[2];
  bn_[3]  = b_[0] + b_[1];
  bn_[4]  = b_[0] - b_[1];
  bn_[5]  = b_[1] + b_[2];
  bn_[6]  = b_[1] - b_[2];
  bn_[7]  = b_[2] + b_[0];
  bn_[8]  = b_[2] - b_[0];
  bn_[9]  = b_[0] + b_[1] + b_[2];
  bn_[10] = b_[0] - b_[1] - b_[2];
  bn_[11] = b_[0] + b_[1] - b_[2];
  bn_[12] = b_[0] - b_[1] + b_[2];
124

Francois Gygi committed
125 126
  for ( int i = 0; i < 13; i++ )
  {
127 128
    an2h_[i] = 0.5 * norm2(an_[i]);
    bn2h_[i] = 0.5 * norm2(bn_[i]);
Francois Gygi committed
129
  }
Francois Gygi committed
130 131 132 133 134 135 136 137 138 139 140 141

  for ( int i = 0; i < 3; i++ )
    a_norm_[i] = length(a_[i]);
  double sp = normalized(a_[1]) * normalized(a_[2]);
  double c = max(-1.0,min(1.0,sp));
  alpha_ = (180.0/M_PI)*acos(c);
  sp = normalized(a_[0]) * normalized(a_[2]);
  c = max(-1.0,min(1.0,sp));
  beta_ = (180.0/M_PI)*acos(c);
  sp = normalized(a_[0]) * normalized(a_[1]);
  c = max(-1.0,min(1.0,sp));
  gamma_ = (180.0/M_PI)*acos(c);
Francois Gygi committed
142
}
143

Francois Gygi committed
144 145 146 147 148 149 150 151 152 153 154 155
////////////////////////////////////////////////////////////////////////////////
bool UnitCell::in_ws(const D3vector& v) const
{
  bool in = true;
  int i = 0;
  while ( i < 13 && in )
  {
    in = ( abs(v*an_[i]) <= an2h_[i] ) ;
    i++;
  }
  return in;
}
156

Francois Gygi committed
157 158 159
////////////////////////////////////////////////////////////////////////////////
void UnitCell::fold_in_ws(D3vector& v) const
{
160
  const double epsilon = 1.e-10;
Francois Gygi committed
161
  bool done = false;
162 163 164
  const int maxiter = 10;
  int iter = 0;
  while ( !done && iter < maxiter )
Francois Gygi committed
165 166 167 168 169
  {
    done = true;
    for ( int i = 0; (i < 13) && done; i++ )
    {
      const double sp = v*an_[i];
170
      if ( sp > an2h_[i] + epsilon )
Francois Gygi committed
171 172 173 174
      {
        done = false;
        do
          v -= an_[i];
175
        while ( v*an_[i] > an2h_[i] + epsilon );
Francois Gygi committed
176
      }
177
      else if ( sp < -an2h_[i] - epsilon )
Francois Gygi committed
178 179 180 181
      {
        done = false;
        do
          v += an_[i];
182
        while ( v*an_[i] < -an2h_[i] - epsilon );
Francois Gygi committed
183 184
      }
    }
185
    iter++;
Francois Gygi committed
186
  }
187
  assert(iter < maxiter);
Francois Gygi committed
188
}
189

Francois Gygi committed
190 191 192 193 194 195 196 197 198 199 200 201
////////////////////////////////////////////////////////////////////////////////
bool UnitCell::in_bz(const D3vector& k) const
{
  bool in = true;
  int i = 0;
  while ( i < 13 && in )
  {
    in = ( abs(k*bn_[i]) <= bn2h_[i] ) ;
    i++;
  }
  return in;
}
202

Francois Gygi committed
203 204 205
////////////////////////////////////////////////////////////////////////////////
void UnitCell::fold_in_bz(D3vector& k) const
{
206
  const double epsilon = 1.e-10;
Francois Gygi committed
207
  bool done = false;
208 209 210
  const int maxiter = 10;
  int iter = 0;
  while ( !done && iter < maxiter )
Francois Gygi committed
211 212 213 214 215
  {
    done = true;
    for ( int i = 0; (i < 13) && done; i++ )
    {
      double sp = k*bn_[i];
216
      if ( sp > bn2h_[i] + epsilon )
Francois Gygi committed
217 218 219 220
      {
        done = false;
        do
          k -= bn_[i];
221
        while ( k*bn_[i] > bn2h_[i] + epsilon );
Francois Gygi committed
222
      }
223
      else if ( sp < -bn2h_[i] - epsilon )
Francois Gygi committed
224 225 226 227
      {
        done = false;
        do
          k += bn_[i];
228
        while ( k*bn_[i] < -bn2h_[i] - epsilon );
Francois Gygi committed
229 230
      }
    }
231
    iter++;
Francois Gygi committed
232
  }
233
  assert(iter < maxiter);
Francois Gygi committed
234
}
235

Francois Gygi committed
236 237 238 239 240 241 242
////////////////////////////////////////////////////////////////////////////////
bool UnitCell::encloses(const UnitCell& c) const
{
  bool in = true;
  int i = 0;
  while ( i < 13 && in )
  {
243 244 245 246
    in = ( contains(c.an_[i]) );
    if ( !in )
      cout << "UnitCell::encloses: " << c.an_[i] << " not in cell "
      << c << endl;
Francois Gygi committed
247 248 249 250 251 252 253 254 255 256 257 258
    i++;
  }
  return in;
}

////////////////////////////////////////////////////////////////////////////////
bool UnitCell::contains(D3vector v) const
{
  const double fac = 0.5 / ( 2.0 * M_PI );
  const double p0 = fac * v * b_[0];
  const double p1 = fac * v * b_[1];
  const double p2 = fac * v * b_[2];
259 260 261
  return ( (p0 > 0.0) && (p0 <= 1.0) &&
           (p1 > 0.0) && (p1 <= 1.0) &&
           (p2 > 0.0) && (p2 <= 1.0) );
Francois Gygi committed
262
}
263

Francois Gygi committed
264 265 266 267
////////////////////////////////////////////////////////////////////////////////
void UnitCell::print(ostream& os) const
{
  os.setf(ios::fixed,ios::floatfield);
Francois Gygi committed
268
  os << setprecision(8);
269 270
  os << "<unit_cell " << endl;
  os << "    a=\"" << setw(12) << a_[0].x << " "
271
                   << setw(12) << a_[0].y << " "
272 273
                   << setw(12) << a_[0].z << "\"" << endl;
  os << "    b=\"" << setw(12) << a_[1].x << " "
274
                   << setw(12) << a_[1].y << " "
275 276
                   << setw(12) << a_[1].z << "\"" << endl;
  os << "    c=\"" << setw(12) << a_[2].x << " "
277
                   << setw(12) << a_[2].y << " "
278
                   << setw(12) << a_[2].z << "\"" << " />" << endl;
Francois Gygi committed
279
}
280

Francois Gygi committed
281 282 283 284 285
////////////////////////////////////////////////////////////////////////////////
bool UnitCell::operator==(const UnitCell& c) const
{
  return ( a_[0]==c.a_[0] && a_[1]==c.a_[1] && a_[2]==c.a_[2] );
}
286

287 288 289 290 291 292 293 294 295 296 297 298
////////////////////////////////////////////////////////////////////////////////
bool UnitCell::operator!=(const UnitCell& c) const
{
  return ! ( *this == c );
}

////////////////////////////////////////////////////////////////////////////////
void UnitCell::vecmult3x3(const double* x, const double* y, double *z) const
{
  //  | z0 |     | x0 x3 x6 |     | y0 |
  //  | z1 |  =  | x1 x4 x7 |  *  | y1 |
  //  | z2 |     | x2 x5 x8 |     | y2 |
299

300 301 302 303 304 305 306
  const double z0 = x[0]*y[0]+x[3]*y[1]+x[6]*y[2];
  const double z1 = x[1]*y[0]+x[4]*y[1]+x[7]*y[2];
  const double z2 = x[2]*y[0]+x[5]*y[1]+x[8]*y[2];
  z[0] = z0;
  z[1] = z1;
  z[2] = z2;
}
307

308 309 310 311 312

////////////////////////////////////////////////////////////////////////////////
void UnitCell::vecsmult3x3(const double* xs, const double* y, double *z) const
{
  // multiply a vector by a symmetric 3x3 matrix
Francois Gygi committed
313

314 315 316
  //  | z0 |     | xs0 xs3 xs5 |     | y0 |
  //  | z1 |  =  | xs3 xs1 xs4 |  *  | y1 |
  //  | z2 |     | xs5 xs4 xs2 |     | y2 |
Francois Gygi committed
317

318 319 320 321 322
  z[0] = xs[0]*y[0]+xs[3]*y[1]+xs[5]*y[2];
  z[1] = xs[3]*y[0]+xs[1]*y[1]+xs[4]*y[2];
  z[2] = xs[5]*y[0]+xs[4]*y[1]+xs[2]*y[2];

}
323 324 325 326 327 328
////////////////////////////////////////////////////////////////////////////////
void UnitCell::matmult3x3(const double* x, const double* y, double *z) const
{
  //  | z0 z3 z6 |     | x0 x3 x6 |     | y0 y3 y6 |
  //  | z1 z4 z7 |  =  | x1 x4 x7 |  *  | y1 y4 y7 |
  //  | z2 z5 z8 |     | x2 x5 x8 |     | y2 y5 y8 |
329

330 331 332
  const double z00 = x[0]*y[0]+x[3]*y[1]+x[6]*y[2];
  const double z10 = x[1]*y[0]+x[4]*y[1]+x[7]*y[2];
  const double z20 = x[2]*y[0]+x[5]*y[1]+x[8]*y[2];
333

334 335 336
  const double z01 = x[0]*y[3]+x[3]*y[4]+x[6]*y[5];
  const double z11 = x[1]*y[3]+x[4]*y[4]+x[7]*y[5];
  const double z21 = x[2]*y[3]+x[5]*y[4]+x[8]*y[5];
337

338 339 340
  const double z02 = x[0]*y[6]+x[3]*y[7]+x[6]*y[8];
  const double z12 = x[1]*y[6]+x[4]*y[7]+x[7]*y[8];
  const double z22 = x[2]*y[6]+x[5]*y[7]+x[8]*y[8];
341

342 343 344
  z[0] = z00;
  z[1] = z10;
  z[2] = z20;
345

346 347 348
  z[3] = z01;
  z[4] = z11;
  z[5] = z21;
349

350 351
  z[6] = z02;
  z[7] = z12;
352
  z[8] = z22;
353
}
354

355 356 357 358 359 360
////////////////////////////////////////////////////////////////////////////////
void UnitCell::smatmult3x3(const double* xs, const double* y, double *z) const
{
  //  | z0 z3 z6 |     | xs0 xs3 xs5 |     | y0 y3 y6 |
  //  | z1 z4 z7 |  =  | xs3 xs1 xs4 |  *  | y1 y4 y7 |
  //  | z2 z5 z8 |     | xs5 xs4 xs2 |     | y2 y5 y8 |
361

362 363 364
  const double z00 = xs[0]*y[0]+xs[3]*y[1]+xs[5]*y[2];
  const double z10 = xs[3]*y[0]+xs[1]*y[1]+xs[4]*y[2];
  const double z20 = xs[5]*y[0]+xs[4]*y[1]+xs[2]*y[2];
365

366 367 368
  const double z01 = xs[0]*y[3]+xs[3]*y[4]+xs[5]*y[5];
  const double z11 = xs[3]*y[3]+xs[1]*y[4]+xs[4]*y[5];
  const double z21 = xs[5]*y[3]+xs[4]*y[4]+xs[2]*y[5];
369

370 371 372
  const double z02 = xs[0]*y[6]+xs[3]*y[7]+xs[5]*y[8];
  const double z12 = xs[3]*y[6]+xs[1]*y[7]+xs[4]*y[8];
  const double z22 = xs[5]*y[6]+xs[4]*y[7]+xs[2]*y[8];
373

374 375 376
  z[0] = z00;
  z[1] = z10;
  z[2] = z20;
377

378 379 380
  z[3] = z01;
  z[4] = z11;
  z[5] = z21;
381

382 383
  z[6] = z02;
  z[7] = z12;
384
  z[8] = z22;
385
}
Francois Gygi committed
386 387
////////////////////////////////////////////////////////////////////////////////
ostream& operator<< ( ostream& os, const UnitCell& cell )
388 389
{
  cell.print(os);
Francois Gygi committed
390 391
  return os;
}