MLWFTransform.C 8.97 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 23 24 25 26 27 28 29 30 31 32 33
// MLWFTransform.C
//
////////////////////////////////////////////////////////////////////////////////

#include <iostream>
#include <iomanip>
#include <complex>
#include <cassert>

#include "MLWFTransform.h"
#include "D3vector.h"
#include "Basis.h"
#include "SlaterDet.h"
#include "UnitCell.h"
#include "jade.h"
#include "blas.h"
using namespace std;

////////////////////////////////////////////////////////////////////////////////
34
MLWFTransform::MLWFTransform(const SlaterDet& sd) : sd_(sd),
Francois Gygi committed
35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61
cell_(sd.basis().cell()), ctxt_(sd.context()),  bm_(BasisMapping(sd.basis()))
{
  a_.resize(6);
  adiag_.resize(6);
  const int n = sd.c().n();
  const int nb = sd.c().nb();
  for ( int k = 0; k < 6; k++ )
  {
    a_[k] = new DoubleMatrix(ctxt_,n,n,nb,nb);
    adiag_[k].resize(n);
  }
  u_ = new DoubleMatrix(ctxt_,n,n,nb,nb);
}

////////////////////////////////////////////////////////////////////////////////
MLWFTransform::~MLWFTransform(void)
{
  for ( int k = 0; k < 6; k++ )
    delete a_[k];
  delete u_;
}

////////////////////////////////////////////////////////////////////////////////
void MLWFTransform::compute_transform(void)
{
  const int maxsweep = 50;
  const double tol = 1.e-8;
62 63

  SlaterDet sdcosx(sd_), sdsinx(sd_),
Francois Gygi committed
64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90
            sdcosy(sd_), sdsiny(sd_),
            sdcosz(sd_), sdsinz(sd_);
  const ComplexMatrix& c = sd_.c();
  ComplexMatrix& ccosx = sdcosx.c();
  ComplexMatrix& csinx = sdsinx.c();
  ComplexMatrix& ccosy = sdcosy.c();
  ComplexMatrix& csiny = sdsiny.c();
  ComplexMatrix& ccosz = sdcosz.c();
  ComplexMatrix& csinz = sdsinz.c();
  // proxy real matrices cr, cc, cs
  DoubleMatrix cr(c);
  DoubleMatrix ccx(ccosx);
  DoubleMatrix csx(csinx);
  DoubleMatrix ccy(ccosy);
  DoubleMatrix csy(csiny);
  DoubleMatrix ccz(ccosz);
  DoubleMatrix csz(csinz);
  vector<complex<double> > zvec(bm_.zvec_size()),
    zvec_cos(bm_.zvec_size()), zvec_sin(bm_.zvec_size()),
    ct(bm_.np012loc()), ct_cos(bm_.np012loc()), ct_sin(bm_.np012loc());

  for ( int i = 0; i < 6; i++ )
  {
    a_[i]->resize(c.n(), c.n(), c.nb(), c.nb());
    adiag_[i].resize(c.n());
  }
  u_->resize(c.n(), c.n(), c.nb(), c.nb());
91

Francois Gygi committed
92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107
  // loop over all local states
  const int np0 = bm_.np0();
  const int np1 = bm_.np1();
  const int np2 = bm_.np2();
  const int np01 = np0 * np1;
  const int np2loc = bm_.np2loc();
  const int nvec = bm_.nvec();
  for ( int n = 0; n < c.nloc(); n++ )
  {
    const complex<double>* f = c.cvalptr(n*c.mloc());
    complex<double>* fcx = ccosx.valptr(n*c.mloc());
    complex<double>* fsx = csinx.valptr(n*c.mloc());
    complex<double>* fcy = ccosy.valptr(n*c.mloc());
    complex<double>* fsy = csiny.valptr(n*c.mloc());
    complex<double>* fcz = ccosz.valptr(n*c.mloc());
    complex<double>* fsz = csinz.valptr(n*c.mloc());
108

Francois Gygi committed
109 110 111
    // direction z
    // map state to array zvec_
    bm_.vector_to_zvec(&f[0],&zvec[0]);
112

Francois Gygi committed
113 114 115 116 117 118
    for ( int ivec = 0; ivec < nvec; ivec++ )
    {
      const int ibase = ivec * np2;
      compute_sincos(np2,&zvec[ibase],&zvec_cos[ibase],&zvec_sin[ibase]);
    }
    // map back zvec_cos to sdcos and zvec_sin to sdsin
119
    bm_.zvec_to_vector(&zvec_cos[0],&fcz[0]);
Francois Gygi committed
120
    bm_.zvec_to_vector(&zvec_sin[0],&fsz[0]);
121 122 123

    // x direction
    // map zvec to ct
Francois Gygi committed
124
    bm_.transpose_fwd(&zvec[0],&ct[0]);
125

Francois Gygi committed
126 127 128 129 130 131 132 133 134
    for ( int iz = 0; iz < np2loc; iz++ )
    {
      for ( int iy = 0; iy < np1; iy++ )
      {
        const int ibase = iz * np01 + iy * np0;
        compute_sincos(np0,&ct[ibase],&ct_cos[ibase],&ct_sin[ibase]);
      }
    }
    // transpose back ct_cos to zvec_cos
135
    bm_.transpose_bwd(&ct_cos[0],&zvec_cos[0]);
Francois Gygi committed
136
    // transpose back ct_sin to zvec_sin
137 138
    bm_.transpose_bwd(&ct_sin[0],&zvec_sin[0]);

Francois Gygi committed
139
    // map back zvec_cos to sdcos and zvec_sin to sdsin
140
    bm_.zvec_to_vector(&zvec_cos[0],&fcx[0]);
Francois Gygi committed
141
    bm_.zvec_to_vector(&zvec_sin[0],&fsx[0]);
142

Francois Gygi committed
143 144 145 146 147 148 149 150 151 152
    // y direction
    vector<complex<double> > c_tmp(np1),ccos_tmp(np1),csin_tmp(np1);
    int one = 1;
    int len = np1;
    int stride = np0;
    for ( int iz = 0; iz < np2loc; iz++ )
    {
      for ( int ix = 0; ix < np0; ix++ )
      {
        const int ibase = iz * np01 + ix;
153
        zcopy(&len,&ct[ibase],&stride,&c_tmp[0],&one);
Francois Gygi committed
154
        compute_sincos(np1,&c_tmp[0],&ccos_tmp[0],&csin_tmp[0]);
155 156
        zcopy(&len,&ccos_tmp[0],&one,&ct_cos[ibase],&stride);
        zcopy(&len,&csin_tmp[0],&one,&ct_sin[ibase],&stride);
Francois Gygi committed
157 158 159
      }
    }
    // transpose back ct_cos to zvec_cos
160
    bm_.transpose_bwd(&ct_cos[0],&zvec_cos[0]);
Francois Gygi committed
161
    // transpose back ct_sin to zvec_sin
162 163
    bm_.transpose_bwd(&ct_sin[0],&zvec_sin[0]);

Francois Gygi committed
164
    // map back zvec_cos and zvec_sin
165
    bm_.zvec_to_vector(&zvec_cos[0],&fcy[0]);
Francois Gygi committed
166
    bm_.zvec_to_vector(&zvec_sin[0],&fsy[0]);
167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185
  }

  // dot products a_[0] = <cos x>, a_[1] = <sin x>
  a_[0]->gemm('t','n',2.0,cr,ccx,0.0);
  a_[0]->ger(-1.0,cr,0,ccx,0);
  a_[1]->gemm('t','n',2.0,cr,csx,0.0);
  a_[1]->ger(-1.0,cr,0,csx,0);

  // dot products a_[2] = <cos y>, a_[3] = <sin y>
  a_[2]->gemm('t','n',2.0,cr,ccy,0.0);
  a_[2]->ger(-1.0,cr,0,ccy,0);
  a_[3]->gemm('t','n',2.0,cr,csy,0.0);
  a_[3]->ger(-1.0,cr,0,csy,0);

  // dot products a_[4] = <cos z>, a_[5] = <sin z>
  a_[4]->gemm('t','n',2.0,cr,ccz,0.0);
  a_[4]->ger(-1.0,cr,0,ccz,0);
  a_[5]->gemm('t','n',2.0,cr,csz,0.0);
  a_[5]->ger(-1.0,cr,0,csz,0);
Francois Gygi committed
186 187 188 189 190

  int nsweep = jade(maxsweep,tol,a_,*u_,adiag_);
}

////////////////////////////////////////////////////////////////////////////////
191
void MLWFTransform::compute_sincos(const int n, const complex<double>* f,
Francois Gygi committed
192 193 194 195 196 197
  complex<double>* fc, complex<double>* fs)
{
  // fc[i] =     0.5 * ( f[i-1] + f[i+1] )
  // fs[i] = (0.5/i) * ( f[i-1] - f[i+1] )

  // i = 0
198 199 200 201 202
  complex<double> zp = f[n-1];
  complex<double> zm = f[1];
  fc[0] = 0.5 * ( zp + zm );
  complex<double> zdiff = zp - zm;
  fs[0] = 0.5 * complex<double>(imag(zdiff),-real(zdiff));
Francois Gygi committed
203 204 205 206 207
  for ( int i = 1; i < n-1; i++ )
  {
    const complex<double> zzp = f[i-1];
    const complex<double> zzm = f[i+1];
    fc[i] = 0.5 * ( zzp + zzm );
208
    const complex<double> zzdiff = zzp - zzm;
Francois Gygi committed
209 210 211
    fs[i] = 0.5 * complex<double>(imag(zzdiff),-real(zzdiff));
  }
  // i = n-1
212 213
  zp = f[n-2];
  zm = f[0];
Francois Gygi committed
214
  fc[n-1] = 0.5 * ( zp + zm );
215
  zdiff = zp - zm;
Francois Gygi committed
216 217
  fs[n-1] = 0.5 * complex<double>(imag(zdiff),-real(zdiff));
}
218

Francois Gygi committed
219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236
////////////////////////////////////////////////////////////////////////////////
D3vector MLWFTransform::center(int i)
{
  assert(i>=0 && i<sd_.nst());
  const double cx = adiag_[0][i];
  const double sx = adiag_[1][i];
  const double cy = adiag_[2][i];
  const double sy = adiag_[3][i];
  const double cz = adiag_[4][i];
  const double sz = adiag_[5][i];
  // Next lines: M_1_PI = 1.0/pi
  const double itwopi = 1.0 / ( 2.0 * M_PI );
  const double t0 = itwopi * atan2(sx,cx);
  const double t1 = itwopi * atan2(sy,cy);
  const double t2 = itwopi * atan2(sz,cz);
  const double x = t0*cell_.a(0).x + t1*cell_.a(1).x + t2*cell_.a(2).x;
  const double y = t0*cell_.a(0).y + t1*cell_.a(1).y + t2*cell_.a(2).y;
  const double z = t0*cell_.a(0).z + t1*cell_.a(1).z + t2*cell_.a(2).z;
237

Francois Gygi committed
238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281
  return D3vector(x,y,z);
}

////////////////////////////////////////////////////////////////////////////////
double MLWFTransform::spread2(int i, int j)
{
  assert(i>=0 && i<sd_.nst());
  assert(j>=0 && j<3);
  const double c = adiag_[2*j][i];
  const double s = adiag_[2*j+1][i];
  // Next line: M_1_PI = 1.0/pi
  const double fac = 1.0 / length(cell_.b(j));
  return fac*fac * ( 1.0 - c*c - s*s );
}

////////////////////////////////////////////////////////////////////////////////
double MLWFTransform::spread2(int i)
{
  assert(i>=0 & i<sd_.nst());
  return spread2(i,0) + spread2(i,1) + spread2(i,2);
}

////////////////////////////////////////////////////////////////////////////////
double MLWFTransform::spread(int i)
{
  return sqrt(spread2(i));
}

////////////////////////////////////////////////////////////////////////////////
double MLWFTransform::spread2(void)
{
  double sum = 0.0;
  for ( int i = 0; i < sd_.nst(); i++ )
    sum += spread2(i);
  return sum;
}

////////////////////////////////////////////////////////////////////////////////
double MLWFTransform::spread(void)
{
  return sqrt(spread2());
}

////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
282
D3vector MLWFTransform::dipole(void)
Francois Gygi committed
283
{
Francois Gygi committed
284
  // total electronic dipole
Francois Gygi committed
285 286
  D3vector sum(0.0,0.0,0.0);
  for ( int i = 0; i < sd_.nst(); i++ )
Francois Gygi committed
287
    sum -= sd_.occ(i) * center(i);
Francois Gygi committed
288 289 290 291 292 293 294 295 296 297 298
  return sum;
}

////////////////////////////////////////////////////////////////////////////////
void MLWFTransform::apply_transform(SlaterDet& sd)
{
  // proxy double matrix c
  DoubleMatrix c(sd.c());
  DoubleMatrix cp(c);
  c.gemm('n','n',1.0,cp,*u_,0.0);
}