StructureFactor.C 4.36 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 18 19 20 21 22
// StructureFactor.C
//
////////////////////////////////////////////////////////////////////////////////

#include "StructureFactor.h"
#include "Basis.h"
#include "UnitCell.h"
#include <cassert>
23
#include <cstring>
24
using namespace std;
Francois Gygi committed
25 26

////////////////////////////////////////////////////////////////////////////////
27
void StructureFactor::init(const vector<vector<double> >& tau,
Francois Gygi committed
28 29 30 31 32 33 34 35 36 37 38
  const Basis& basis)
{
  _k0min = basis.idxmin(0);
  _k1min = basis.idxmin(1);
  _k2min = basis.idxmin(2);
  _k0max = basis.idxmax(0);
  _k1max = basis.idxmax(1);
  _k2max = basis.idxmax(2);
  _k0range = _k0max - _k0min + 1;
  _k1range = _k1max - _k1min + 1;
  _k2range = _k2max - _k2min + 1;
39

Francois Gygi committed
40
  // get dimensions of tau[nsp][3*na[is]]
41

Francois Gygi committed
42 43
  _nsp = tau.size();
  _na.resize(_nsp);
44

Francois Gygi committed
45 46 47 48 49 50 51
  cos0.resize(_nsp);
  cos1.resize(_nsp);
  cos2.resize(_nsp);
  sin0.resize(_nsp);
  sin1.resize(_nsp);
  sin2.resize(_nsp);
  sfac.resize(_nsp);
52

Francois Gygi committed
53
  _ng = basis.localsize();
54

Francois Gygi committed
55 56 57 58 59 60 61 62 63 64 65 66 67 68 69
  for ( int is = 0; is < _nsp; is++ )
  {
    assert( tau[is].size() % 3 == 0 );
    _na[is] = tau[is].size() / 3; // tau[nsp][3*na[is]]
    cos0[is].resize(_na[is]*_k0range);
    cos1[is].resize(_na[is]*_k1range);
    cos2[is].resize(_na[is]*_k2range);
    sin0[is].resize(_na[is]*_k0range);
    sin1[is].resize(_na[is]*_k1range);
    sin2[is].resize(_na[is]*_k2range);
    sfac[is].resize(_ng);
  }
}

////////////////////////////////////////////////////////////////////////////////
70
void StructureFactor::update(const vector<vector<double> >& tau,
Francois Gygi committed
71 72
  const Basis& basis)
{
73
  // it is assumed that the dimensions of tau and the basis have
Francois Gygi committed
74
  // not changed since the last call to StructureFactor::init
75 76 77

  // check that number of species has not changed
  assert(tau.size() == _nsp);
Francois Gygi committed
78
  assert(basis.localsize() == _ng);
79

Francois Gygi committed
80
  const int * const idx = basis.idx_ptr();
81

Francois Gygi committed
82 83 84 85 86 87 88 89 90
  const UnitCell& cell = basis.cell();
  const D3vector b0 = cell.b(0);
  const D3vector b1 = cell.b(1);
  const D3vector b2 = cell.b(2);

  for ( int is = 0; is < _nsp; is++ )
  {
    assert( 3 * _na[is] == tau[is].size() );
    memset( (void*)&sfac[is][0], 0, 2*_ng*sizeof(double) );
91

Francois Gygi committed
92 93 94 95 96 97 98 99
    for ( int ia = 0; ia < _na[is]; ia++ )
    {
      double *c0 = cos0_ptr(is,ia);
      double *c1 = cos1_ptr(is,ia);
      double *c2 = cos2_ptr(is,ia);
      double *s0 = sin0_ptr(is,ia);
      double *s1 = sin1_ptr(is,ia);
      double *s2 = sin2_ptr(is,ia);
100

Francois Gygi committed
101
      const double * const tauptr = &tau[is][3*ia];
102

Francois Gygi committed
103
      const D3vector t(tauptr[0],tauptr[1],tauptr[2]);
104

Francois Gygi committed
105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133
      /* x direction */
      const double fac0 = b0 * t;

      for ( int i = _k0min; i < _k0max+1; i++ )
      {
        const double arg = i * fac0;
        c0[i] = cos(arg);
        s0[i] = sin(arg);
      }

      /* y direction */
      const double fac1 = b1 * t;

      for ( int i = _k1min; i < _k1max+1; i++ )
      {
        const double arg = i * fac1;
        c1[i] = cos(arg);
        s1[i] = sin(arg);
      }

      /* z direction */
      const double fac2 = b2 * t;

      for ( int i = _k2min; i < _k2max+1; i++ )
      {
        const double arg = i * fac2;
        c2[i] = cos(arg);
        s2[i] = sin(arg);
      }
134 135 136

      // compute sfac[is][i]

Francois Gygi committed
137 138 139 140 141 142
      for ( int i = 0; i < _ng; i++ )
      {
        const int iii = i+i+i;
        const int kx = idx[iii];
        const int ky = idx[iii+1];
        const int kz = idx[iii+2];
143

Francois Gygi committed
144 145 146
        const double cos_a = c0[kx];
        const double cos_b = c1[ky];
        const double cos_c = c2[kz];
147

Francois Gygi committed
148 149 150
        const double sin_a = s0[kx];
        const double sin_b = s1[ky];
        const double sin_c = s2[kz];
151 152

        // Next line: exp(-i*gr) =
Francois Gygi committed
153 154
        // (cos_a - I sin_a)*(cos_b - I sin_b)*(cos_c - I sin_c)
        sfac[is][i] += complex<double>(
155
        cos_a*cos_b*cos_c - sin_a*sin_b*cos_c -
Francois Gygi committed
156
        sin_a*cos_b*sin_c - cos_a*sin_b*sin_c,
157 158
        sin_a*sin_b*sin_c - sin_a*cos_b*cos_c -
        cos_a*sin_b*cos_c - cos_a*cos_b*sin_c );
Francois Gygi committed
159 160 161 162
      }
    }
  }
}