ElectricEnthalpy.C 21.3 KB
Newer Older
Francois Gygi committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41
////////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2014 The Regents of the University of California
//
// This file is part of Qbox
//
// Qbox is distributed under the terms of the GNU General Public License
// as published by the Free Software Foundation, either version 2 of
// 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/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// ElectricEnthalpy.C
//
////////////////////////////////////////////////////////////////////////////////

#include <iostream>
#include <iomanip>
#include <complex>
#include <cassert>
#include <cmath>
#include <algorithm> // fill

#include "Timer.h"
#include "Context.h"
#include "Matrix.h"
#include "Sample.h"
#include "D3vector.h"
#include "ElectricEnthalpy.h"
#include "MLWFTransform.h"
#include "Wavefunction.h"
#include "D3vector.h"
#include "Basis.h"
#include "SlaterDet.h"
#include "ChargeDensity.h"
#include "FourierTransform.h"
#include "UnitCell.h"
using namespace std;

42
///////////////////////////////////////////////////////////////////////////////
43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70
double ElectricEnthalpy::vsst(double x) const
{
  // smooth sawtooth periodic potential function
  // x in [-1/2, 1/2]
  // The slope of vsst is 1 at x=0
  //
  // The function vsst approximates the identity function x->x
  // in the interval [-1/2+xcut, 1/2-xcut]
  const double xcut = 0.05;
  const double xcut2 = xcut*xcut;
  // The function vsst is well represented by a
  // discrete Fourier transform of length np = 2*ng
  // Some aliasing error will arise if np < 2*ng
  // The error is determined by the product xcut*ng
  const int ng = 32;
  double v = 0.0;
  for ( int ig = 1; ig < ng; ig++ )
  {
    const double g = 2 * M_PI * ig;
    const double arg = -0.25 * g * g * xcut2;
    // next line: factor sgn to shift origin by 0.5
    const int sgn = 1 - 2*(ig%2);
    const double c = -2.0 * sgn * exp ( arg ) / g;
    v += c * sin(x*g);
  }
  return v;
}

Francois Gygi committed
71
///////////////////////////////////////////////////////////////////////////////
72
ElectricEnthalpy::ElectricEnthalpy(Sample& s): s_(s), wf_(s.wf),
Francois Gygi committed
73
  sd_(*(s.wf.sd(0,0))), ctxt_(s.wf.sd(0,0)->context()),
74
  basis_(s.wf.sd(0,0)->basis())
Francois Gygi committed
75 76 77
{
  onpe0_ = ctxt_.onpe0();
  e_field_ = s.ctrl.e_field;
78
  finite_field_ = norm2(e_field_) != 0.0;
Francois Gygi committed
79
  compute_quadrupole_ = false;
80

81 82 83
  if ( s.ctrl.polarization == "OFF" )
    pol_type_ = off;
  else if ( s.ctrl.polarization == "BERRY" )
84
    pol_type_ = berry;
85
  else if ( s.ctrl.polarization == "MLWF" )
86
    pol_type_ = mlwf;
87 88 89 90
  else if ( s.ctrl.polarization == "MLWF_REF" )
    pol_type_ = mlwf_ref;
  else if ( s.ctrl.polarization == "MLWF_REF_Q" )
  {
91
    pol_type_ = mlwf_ref;
92 93
    compute_quadrupole_ = true;
  }
94 95
  else
  {
96
    cerr << "ElectricEnthalpy: invalid polarization option" << endl;
97 98
    ctxt_.abort(1);
  }
Francois Gygi committed
99

100 101 102 103 104 105 106
  // do not allocate further objects if polarization is off
  if ( pol_type_ == off ) return;

  assert(wf_.nkp()==1);
  assert(wf_.nspin()==1);

  dwf_ = new Wavefunction(s.wf);
Francois Gygi committed
107 108 109 110 111
  mlwft_ = new MLWFTransform(sd_);
  mlwft_->set_tol(1.e-10);

  smat_[0] = smat_[1] = smat_[2] = 0;
  rwf_[0] = rwf_[1] = rwf_[2] = 0;
112
  int nst = sd_.nst();
Francois Gygi committed
113

114
  if ( pol_type_ == mlwf_ref || pol_type_ == mlwf_ref_q )
Francois Gygi committed
115
  {
116
    // allocate real space wf arrays for MLWF refinement
117 118
    for ( int idir = 0; idir < 3; idir++ )
      rwf_[idir] = new Wavefunction(wf_);
119
    correction_.resize(nst);
Francois Gygi committed
120
  }
121
  else if ( pol_type_ == berry )
Francois Gygi committed
122 123 124 125
  {
    // allocate complex Berry phase matrix
    int n = sd_.c().n();
    int nb = sd_.c().nb();
126 127
    for ( int idir = 0; idir < 3; idir++ )
      smat_[idir] = new ComplexMatrix(ctxt_,n,n,nb,nb);
Francois Gygi committed
128 129
  }

130
  if ( (pol_type_ != off) && onpe0_ )
Francois Gygi committed
131 132 133 134 135 136 137 138 139 140 141 142 143 144 145
  {
    cout.setf(ios::fixed,ios::floatfield);
    cout.setf(ios::right,ios::adjustfield);
    cout.precision(8);
    cout << "<e_field> " << e_field_ << " </e_field>" << endl;
  }

  mlwfc_.resize(nst);
  mlwfs_.resize(nst);
  quad_.resize(nst);
}

///////////////////////////////////////////////////////////////////////////////
ElectricEnthalpy::~ElectricEnthalpy(void)
{
146 147 148
  if ( pol_type_ == off ) return;

  delete dwf_;
Francois Gygi committed
149
  delete mlwft_;
150
  for ( int idir = 0; idir < 3; idir++ )
Francois Gygi committed
151
  {
152 153
    delete rwf_[idir];
    delete smat_[idir];
Francois Gygi committed
154 155 156 157 158 159 160 161 162
  }

  for ( TimerMap::iterator i = tmap.begin(); i != tmap.end(); i++ )
  {
    double time = (*i).second.real();
    double tmin = time;
    double tmax = time;
    s_.ctxt_.dmin(1,1,&tmin,1);
    s_.ctxt_.dmax(1,1,&tmax,1);
163
    if ( pol_type_ != off && s_.ctxt_.myproc()==0 )
Francois Gygi committed
164
    {
165 166
      string s = "name=\"" + (*i).first + "\"";
      cout << "<timing " << left << setw(22) << s
167 168
           << " min=\"" << setprecision(3) << tmin << "\""
           << " max=\"" << setprecision(3) << tmax << "\"/>"
Francois Gygi committed
169 170 171 172
           << endl;
    }
  }
}
173

Francois Gygi committed
174 175 176
///////////////////////////////////////////////////////////////////////////////
void ElectricEnthalpy::update(void)
{
177 178
  if ( pol_type_ == off ) return;

179
  const UnitCell& cell = sd_.basis().cell();
Francois Gygi committed
180 181 182 183 184 185 186 187 188 189 190 191
  // compute cos and sin matrices
  tmap["mlwf_update"].start();
  mlwft_->update();
  tmap["mlwf_update"].stop();
  vector<SlaterDet*> sdcos(3), sdsin(3);
  sdcos[0] = mlwft_->sdcosx();
  sdcos[1] = mlwft_->sdcosy();
  sdcos[2] = mlwft_->sdcosz();
  sdsin[0] = mlwft_->sdsinx();
  sdsin[1] = mlwft_->sdsiny();
  sdsin[2] = mlwft_->sdsinz();

192 193
  dipole_ion_ = s_.atoms.dipole();
  dipole_el_ = D3vector(0,0,0);
194

195
  if ( pol_type_ == mlwf || pol_type_ == mlwf_ref || pol_type_ == mlwf_ref_q )
Francois Gygi committed
196 197 198 199 200 201 202 203 204 205 206
  {
    tmap["mlwf_trans"].start();
    mlwft_->compute_transform();
    mlwft_->apply_transform(sd_);
    tmap["mlwf_trans"].stop();
    for ( int i = 0; i < sd_.nst(); i++ )
    {
      mlwfc_[i] = mlwft_->center(i);
      mlwfs_[i] = mlwft_->spread(i);
    }

207
    if ( pol_type_ == mlwf_ref || pol_type_ == mlwf_ref_q )
208
    {
209 210 211
      tmap["correction"].start();
      compute_correction();
      tmap["correction"].stop();
212
    }
Francois Gygi committed
213

214 215
    for ( int i = 0; i < sd_.nst(); i++ )
    {
216
      dipole_el_ -= 2.0 * mlwfc_[i];
217
      if ( pol_type_ == mlwf_ref || pol_type_ == mlwf_ref_q )
218
        dipole_el_ -= 2.0 * correction_[i];
219 220 221
    }

    // compute gradient
222
    if ( finite_field_ )
Francois Gygi committed
223
    {
224 225
      dwf_->clear();
      for ( int idir = 0; idir < 3; idir++ )
Francois Gygi committed
226
      {
227
        if ( e_field_[idir] != 0.0 )
228
        {
229 230
          // MLWF part
          if ( pol_type_ == mlwf )
231
          {
232 233 234 235 236 237
            const double nst = sd_.nst();
            std::vector<double> adiag_inv_real(nst,0),adiag_inv_imag(nst,0);
            for ( int ist = 0; ist < nst; ist ++ )
            {
              const std::complex<double>
                adiag( mlwft_->adiag(idir*2,ist),mlwft_->adiag(idir*2+1,ist) );
238

239 240
              adiag_inv_real[ist] = real( std::complex<double>(1,0) / adiag );
              adiag_inv_imag[ist] = imag( std::complex<double>(1,0) / adiag );
241

242
            }
243

244 245 246
            DoubleMatrix ccos(sdcos[idir]->c());
            DoubleMatrix csin(sdsin[idir]->c());
            DoubleMatrix cp(dwf_->sd(0,0)->c());
247

248 249 250
            int nloc = cp.nloc();
            int mloc = cp.mloc();
            int ione = 1;
251

252
            const double fac = cell.a_norm(idir)
253
                               * e_field_[idir] / ( 2.0 * M_PI );
254

255 256 257 258 259
            for (int in = 0; in < nloc; in++)
            {
              int ist = cp.jglobal(in);
              double fac1 = adiag_inv_real[ist] * fac;
              double fac2 = adiag_inv_imag[ist] * fac;
260

261 262 263
              double *ptr1 = &cp[in*mloc],
                     *ptrcos = &ccos[in*mloc],
                     *ptrsin = &csin[in*mloc];
264

265 266
              daxpy(&mloc, &fac2, ptrcos, &ione, ptr1, &ione);
              daxpy(&mloc, &fac1, ptrsin, &ione, ptr1, &ione);
267

268
            }
269
          }
270 271 272 273 274 275 276 277 278 279 280 281 282 283
          else if ( pol_type_ == mlwf_ref || pol_type_ == mlwf_ref_q )
          {
            // MLWF_REF part: real-space correction
            DoubleMatrix cc(rwf_[idir]->sd(0,0)->c());
            DoubleMatrix cp(dwf_->sd(0,0)->c());

            int size = cc.size();
            double alpha = e_field_[idir];
            int ione = 1;
            daxpy (&size, &alpha, cc.valptr(), &ione, cp.valptr(), &ione);
          } // if pol_type_
        } // if e_field_[idir]
      } // for idir
    } // if finite_field_
Francois Gygi committed
284
  }
285
  else if ( pol_type_ == berry )
Francois Gygi committed
286
  {
287 288 289 290
    // proxy matrix
    DoubleMatrix gradp(dwf_->sd(0,0)->c());
    if ( finite_field_ )
      dwf_->clear();
Francois Gygi committed
291 292 293

    for ( int idir = 0; idir < 3; idir++ )
    {
294
      const double fac = cell.a_norm(idir)/( 2.0*M_PI );
295
      complex<double>* val = smat_[idir]->valptr();
Francois Gygi committed
296

297 298 299 300
      const double* re = mlwft_->a(idir*2)->cvalptr();
      const double* im = mlwft_->a(idir*2+1)->cvalptr();
      for ( int i = 0; i < smat_[idir]->size(); i++ )
        val[i] = complex<double>(re[i],im[i]);
Francois Gygi committed
301

302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340
      // compute determinant of S
      ComplexMatrix& sm = *smat_[idir];
      const Context& ctxt = sm.context();

      // perform LU decomposition of S
      valarray<int> ipiv;
      sm.lu(ipiv);

      int n = sm.n();
      int nb = sm.nb();
      // note: m == n, mb == nb

      // compute determinant from diagonal values of U  (stored in diag of S)
      // get full diagonal of the matrix in array diag
      valarray<complex<double> > diag(n);
      for ( int ii = 0; ii < n; ii++ )
      {
        int iii = sm.l(ii) * nb + sm.x(ii);
        int jjj = sm.m(ii) * nb + sm.y(ii);
        if ( sm.pr(ii) == ctxt.myrow() &&
             sm.pc(ii) == ctxt.mycol() )
          diag[ii] = val[iii+sm.mloc()*jjj];
      }
      ctxt.dsum(2*n,1,(double*)&diag[0],2*n);

      // sum the argument of diagonal elements
      double sumarg = 0.0;
      for ( int ii = 0; ii < n; ii++ )
        sumarg += arg(diag[ii]);

      // compute determinant
      complex<double> det = 1.0;
      for ( int ii = 0; ii < n; ii++ )
        det *= diag[ii];

      const int sign = sm.signature(ipiv);
      if ( sign == -1 )
        sumarg += M_PI;

341
      // assume occupation number of 2.0
342
      dipole_el_[idir] = - 2.0 * fac * sumarg;
343

344
      if ( finite_field_ )
345
      {
346 347
        // compute inverse of smat
        sm.inverse_from_lu(ipiv);
Francois Gygi committed
348 349 350 351

        // real and img part of S^{-1}
        DoubleMatrix s_real(ctxt_,n,n,nb,nb);
        DoubleMatrix s_img(ctxt_,n,n,nb,nb);
352
        DoubleMatrix sp(sm);
Francois Gygi committed
353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369

        int size = s_real.size();
        int ione = 1, itwo = 2;

        // copy real and imaginary parts of s to s_real and s_img
        dcopy(&size,sp.valptr(),&itwo,s_real.valptr(),&ione);
        dcopy(&size,sp.valptr()+1,&itwo,s_img.valptr(),&ione);

        // proxy Matrix for cosx|psi> and sinx|psi>
        DoubleMatrix cosp(sdcos[idir]->c());
        DoubleMatrix sinp(sdsin[idir]->c());

        // alpha = E_i * L_i / 2pi
        double alpha = fac * e_field_[idir];

        gradp.gemm('n','n',alpha,cosp,s_img,1.0);
        gradp.gemm('n','n',alpha,sinp,s_real,1.0);
370 371
      }
    } // for idir
372
  }
373

374 375 376 377
  dipole_total_ = dipole_ion_ + dipole_el_;
  cell.fold_in_ws(dipole_ion_);
  cell.fold_in_ws(dipole_el_);
  cell.fold_in_ws(dipole_total_);
Francois Gygi committed
378
}
379

Francois Gygi committed
380
///////////////////////////////////////////////////////////////////////////////
381
double ElectricEnthalpy::enthalpy(Wavefunction& dwf, bool compute_hpsi)
Francois Gygi committed
382
{
383 384 385 386
  // return zero if polarization is off, or field is zero
  if ( pol_type_ == off || !finite_field_ )
    return 0.0;

387
  enthalpy_ = - e_field_ * dipole_total_;
Francois Gygi committed
388 389 390 391
  if ( compute_hpsi )
  {
    // assert gamma-only and no spin
    assert(dwf.nkp()==1 && dwf.nspin()==1);
392
    dwf.sd(0,0)->c() += dwf_->sd(0,0)->c();
Francois Gygi committed
393
  }
394
  return enthalpy_;
Francois Gygi committed
395 396 397 398 399
}

///////////////////////////////////////////////////////////////////////////////
// Correction scheme by M. Stengel and N. Spaldin,
// Phys. Rev. B 73, 075121 (2006)
400
// Calculate correction in real space and derivatives of the correction
Francois Gygi committed
401
///////////////////////////////////////////////////////////////////////////////
402
void ElectricEnthalpy::compute_correction(void)
Francois Gygi committed
403
{
404 405 406
  int np0v = basis_.np(0);
  int np1v = basis_.np(1);
  int np2v = basis_.np(2);
Francois Gygi committed
407 408 409
  const ComplexMatrix& c = sd_.c();
  DoubleMatrix cp(c);

410
  FourierTransform ft(basis_,np0v,np1v,np2v);
Francois Gygi committed
411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436

  int np012v = ft.np012();
  int np012loc = ft.np012loc();
  int nst = sd_.nst();
  int nloc = c.nloc();
  int mloc = c.mloc();

  // store (x-x0)|psi> in rwfs
  rwf_[0]->clear();
  rwf_[1]->clear();
  rwf_[2]->clear();

  ComplexMatrix& cx = rwf_[0]->sd(0,0)->c();
  ComplexMatrix& cy = rwf_[1]->sd(0,0)->c();
  ComplexMatrix& cz = rwf_[2]->sd(0,0)->c();

  DoubleMatrix cpx(cx);
  DoubleMatrix cpy(cy);
  DoubleMatrix cpz(cz);

  // calculate refinements
  // ref is scaled by np012v
  vector<double> ref(nst*3);
  if ( compute_quadrupole_ ) ref.resize(nst*9);

  // cell size;
437 438 439 440
  const UnitCell& cell = sd_.basis().cell();
  const double ax = cell.amat(0);
  const double ay = cell.amat(4);
  const double az = cell.amat(8);
Francois Gygi committed
441 442 443 444 445 446 447 448 449 450 451 452

  // half cell size;
  const double ax2 = ax / 2.0;
  const double ay2 = ay / 2.0;
  const double az2 = az / 2.0;

  // grid size;
  const double dx = ax / np0v;
  const double dy = ay / np1v;
  const double dz = az / np2v;

  for ( int i = 0; i < nst; i++ )
453
    correction_[i] = D3vector(0,0,0);
Francois Gygi committed
454

455 456 457 458 459 460 461 462
  int niter = 1;
  // check if override from the debug variable
  // use: set debug MLWF_REF_NITER <value>
  if ( s_.ctrl.debug.find("MLWF_REF_NITER") != string::npos )
  {
    istringstream is(s_.ctrl.debug);
    string s;
    is >> s >> niter;
463
    assert(s=="MLWF_REF_NITER");
464
    if ( onpe0_ )
Francois Gygi committed
465
      cout << " ElectricEnthalpy: override niter value: niter= "
466 467 468 469
           << niter << endl;
    assert(niter > 0);
  }
  for ( int iter = 0; iter < niter; iter++ )
Francois Gygi committed
470 471 472
  {
    fill(ref.begin(),ref.end(),0.0);

473 474 475 476 477 478
    // wavefunctions in real space
    vector<complex<double> > wftmp(np012loc);
    vector<complex<double> > xwftmp(np012loc);
    vector<complex<double> > ywftmp(np012loc);
    vector<complex<double> > zwftmp(np012loc);

Francois Gygi committed
479 480 481 482 483 484 485 486 487
    for ( int in = 0; in < nloc; in++ )
    {
      int n = c.jglobal(in);
      double* pref;
      if ( compute_quadrupole_ )
        pref = &ref[9*n];
      else
        pref = &ref[3*n];

488
      // real space wavefunction in wftmp
Francois Gygi committed
489 490 491 492 493
      tmap["ft"].start();
      ft.backward(c.cvalptr(mloc*in),&wftmp[0]);
      tmap["ft"].stop();

      tmap["real"].start();
494 495 496
      double x0 = mlwfc_[n][0] + correction_[n][0];
      double y0 = mlwfc_[n][1] + correction_[n][1];
      double z0 = mlwfc_[n][2] + correction_[n][2];
497

498 499 500
      // compute shifted sawtooth potentials vx, vy, vz
      vector<double> vx(ft.np0()), vy(ft.np1()), vz(ft.np2());
      for ( int i = 0; i < vx.size(); i++ )
Francois Gygi committed
501
      {
502 503
        double x = i * dx - x0;
        if ( x >  ax2 ) x -= ax;
Francois Gygi committed
504
        if ( x < -ax2 ) x += ax;
505
#ifdef NO_VSST
506
        vx[i] = x;
507 508 509
#else
        vx[i] = ax * vsst(x/ax);
#endif
510 511 512 513 514
      }
      for ( int j = 0; j < vy.size(); j++ )
      {
        double y = j * dy - y0;
        if ( y >  ay2 ) y -= ay;
Francois Gygi committed
515
        if ( y < -ay2 ) y += ay;
516
#ifdef NO_VSST
517
        vy[j] = y;
518 519 520
#else
        vy[j] = ay * vsst(y/ay);
#endif
521 522 523 524 525 526
      }
      for ( int k = 0; k < vz.size(); k++ )
      {
        double z = k * dz - z0;
        if ( z >  az2 ) z -= az;
        if ( z < -az2 ) z += az;
527
#ifdef NO_VSST
528
        vz[k] = z;
529 530 531
#else
        vz[k] = az * vsst(z/az);
#endif
532 533 534 535
      }

      for ( int i = 0; i < np012loc; i++ )
      {
536 537 538
        int ix = ft.i(i);
        int iy = ft.j(i);
        int iz = ft.k(i);
Francois Gygi committed
539

540 541 542 543
        const double wft = real(wftmp[i]);
        const double xwft = wft * vx[ix];
        const double ywft = wft * vy[iy];
        const double zwft = wft * vz[iz];
Francois Gygi committed
544

545 546 547
        pref[0] += wft * xwft;
        pref[1] += wft * ywft;
        pref[2] += wft * zwft;
Francois Gygi committed
548

549
        xwftmp[i] = xwft;
550 551
        ywftmp[i] = ywft;
        zwftmp[i] = zwft;
Francois Gygi committed
552 553 554

        if ( compute_quadrupole_ )
        {
555 556 557 558 559 560
          pref[3] += xwft * xwft;
          pref[4] += ywft * ywft;
          pref[5] += zwft * zwft;
          pref[6] += xwft * ywft;
          pref[7] += ywft * zwft;
          pref[8] += zwft * xwft;
Francois Gygi committed
561
        }
562
      } // for i
Francois Gygi committed
563 564
      tmap["real"].stop();

565
      // ft to get xwf in reciprocal space at the last iteration
566
      if ( iter == niter - 1 )
Francois Gygi committed
567 568
      {
        tmap["ft"].start();
569 570 571 572 573 574
        if ( e_field_[0] != 0.0 )
          ft.forward(&xwftmp[0],cx.valptr(mloc*in));
        if ( e_field_[1] != 0.0 )
          ft.forward(&ywftmp[0],cy.valptr(mloc*in));
        if ( e_field_[2] != 0.0 )
          ft.forward(&zwftmp[0],cz.valptr(mloc*in));
Francois Gygi committed
575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591
        tmap["ft"].stop();
      } // if
    } //for in

    ctxt_.barrier();
    tmap["dsum"].start();
    if ( compute_quadrupole_ )
      ctxt_.dsum(9*nst,1,&ref[0],9*nst);
    else
      ctxt_.dsum(3*nst,1,&ref[0],3*nst);
    tmap["dsum"].stop();

    tmap["real"].start();
    if ( compute_quadrupole_ )
    {
      for ( int ist = 0; ist < nst; ist++ )
      {
592
        D3vector& pcor = correction_[ist];
Francois Gygi committed
593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609
        D3tensor& pquad = quad_[ist];

        pcor[0] += ref[ist*9]/np012v;
        pcor[1] += ref[ist*9+1]/np012v;
        pcor[2] += ref[ist*9+2]/np012v;
        pquad.setdiag ( 0, ref[ist*9+3]/np012v - pcor[0] * pcor[0] );
        pquad.setdiag ( 1, ref[ist*9+4]/np012v - pcor[1] * pcor[1] );
        pquad.setdiag ( 2, ref[ist*9+5]/np012v - pcor[2] * pcor[2] );
        pquad.setoffdiag ( 0, ref[ist*9+6]/np012v - pcor[0] * pcor[1] );
        pquad.setoffdiag ( 1, ref[ist*9+7]/np012v - pcor[1] * pcor[2] );
        pquad.setoffdiag ( 2, ref[ist*9+8]/np012v - pcor[2] * pcor[0] );
      }
    }
    else
    {
      for ( int ist = 0; ist < nst; ist++ )
      {
610
        D3vector& pcor = correction_[ist];
Francois Gygi committed
611 612 613 614 615 616
        pcor[0] += ref[ist*3]/np012v;
        pcor[1] += ref[ist*3+1]/np012v;
        pcor[2] += ref[ist*3+2]/np012v;
      }
    }
    tmap["real"].stop();
617
  } // for iter
Francois Gygi committed
618
}
619 620 621 622

////////////////////////////////////////////////////////////////////////////////
void ElectricEnthalpy::print(ostream& os) const
{
623 624
  if ( pol_type_ == off ) return;

625
  os << fixed << right << setprecision(8);
626 627
  // print MLWF centers if pol_type_ == MLWF or MLWF_REF or MLWF_REF_Q
  if ( pol_type_ == mlwf || pol_type_ == mlwf_ref || pol_type_ == mlwf_ref_q )
628
  {
629 630 631 632 633
    int nst = sd_.nst();
    os << " <mlwf_set size=\"" << nst << "\">" << endl;
    for ( int i = 0; i < nst; i++ )
    {
      os << " <mlwf center=\"" << setprecision(8)
634 635 636 637
         << setw(12) << mlwfc_[i].x << " "
         << setw(12) << mlwfc_[i].y << " "
         << setw(12) << mlwfc_[i].z << " \"\n"
         << "       spread=\" " << mlwfs_[i] << " \"/>" << endl;
638
      if ( pol_type_ == mlwf_ref )
639 640
      {
        os << " <mlwf_ref center=\"" << setprecision(8)
641 642
           << setw(12) << mlwfc_[i].x + correction_[i].x << " "
           << setw(12) << mlwfc_[i].y + correction_[i].y << " "
643
           << setw(12) << mlwfc_[i].z + correction_[i].z << " \"";
644 645 646
        if ( compute_quadrupole_ )
        {
          // add spread attribute
647
          os << " \n     spread=\" " << sqrt(quad_[i].trace()) << " \"";
648 649 650 651 652
        }
        os << "/>" << endl;

        if ( compute_quadrupole_ )
          os << "    <quad>"
653 654 655 656 657
             << setw(12) << quad_[i][0] << " "
             << setw(12) << quad_[i][4] << " "
             << setw(12) << quad_[i][8] << " "
             << setw(12) << quad_[i][1] << " "
             << setw(12) << quad_[i][2] << " "
658 659 660 661 662
             << setw(12) << quad_[i][5]
             << " </quad>" << endl;
      }
    }
    os << " </mlwf_set>" << endl;
663
  }
664

665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686
  // print dipole
  os << setprecision(10) << fixed << right;
  os << " <dipole>\n";
  os << "   <dipole_ion>   "
     << setw(14) << dipole_ion_.x << " "
     << setw(14) << dipole_ion_.y << " "
     << setw(14) << dipole_ion_.z << " </dipole_ion>\n";
  os << "   <dipole_el>    "
     << setw(14) << dipole_el_.x << " "
     << setw(14) << dipole_el_.y << " "
     << setw(14) << dipole_el_.z << " </dipole_el>\n";
  os << "   <dipole_total> "
     << setw(14) << dipole_total_.x << " "
     << setw(14) << dipole_total_.y << " "
     << setw(14) << dipole_total_.z << " </dipole_total>\n";
  os << " </dipole>\n";

  if ( compute_quadrupole_ )
  {
    D3tensor q_mlwfc;
    D3tensor q_mlwfs;
    for ( int ist = 0; ist < sd_.nst(); ist++ )
687
    {
688 689
      D3vector ctr = mlwfc_[ist];
      for (int j=0; j<3; j++)
690
      {
691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728
        for (int k = 0; k < 3; k++)
          q_mlwfc[j*3+k] -= 2.0 * ctr[j] * ctr[k];
      }
      q_mlwfs -= quad_[ist] * 2.0;
    } // for ist

    D3tensor q_ion = s_.atoms.quadrupole();
    D3tensor q_mlwf = q_mlwfc + q_mlwfs;
    //total primitive quadrupoles
    D3tensor q_total = q_ion + q_mlwf;
    //traceless quadrupole
    D3tensor q_traceless = q_total;
    q_traceless.traceless();

    os << " <quadrupole> " << endl;
    os << "   <quadrupole_ion> " << endl
       << q_ion
       << "   </quadrupole_ion>" << endl;
    os << "   <quadrupole_el> " << endl
       << q_mlwf
       << "   </quadrupole_el>" << endl;
    os << "   <quadrupole_total> " << endl
       << q_total
       << "   </quadrupole_total>" << endl;
    os << "   <traceless_quadrupole> " << endl
       << q_traceless
       << "   </traceless_quadrupole>" << endl;
    char uplo = 'u';
    D3vector eigval;
    D3tensor eigvec;
    q_traceless.syev(uplo, eigval, eigvec);
    os << "   <traceless_quadrupole_eigval> " << endl
       << "    " << eigval << endl
       << "   </traceless_quadrupole_eigval>" << endl;
    os << "   <traceless_quadrupole_eigvec> " << endl
       << eigvec
       << "   </traceless_quadrupole_eigvec>" << endl;
    os << " </quadrupole> " << endl;
729 730
  }

731 732
}

733 734 735 736 737 738 739
////////////////////////////////////////////////////////////////////////////////
void ElectricEnthalpy::set_e_field(D3vector e_field_val)
{
  e_field_ = e_field_val;
  finite_field_ = norm2(e_field_) != 0.0;
}

740 741 742 743 744 745
////////////////////////////////////////////////////////////////////////////////
ostream& operator<< ( ostream& os, const ElectricEnthalpy& e )
{
  e.print(os);
  return os;
}