EnergyFunctional.C 32.8 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
// EnergyFunctional.C
//
////////////////////////////////////////////////////////////////////////////////

#include "EnergyFunctional.h"
#include "Sample.h"
21
#include "Species.h"
Francois Gygi committed
22
#include "Wavefunction.h"
23
#include "ChargeDensity.h"
Francois Gygi committed
24 25 26 27
#include "SlaterDet.h"
#include "Basis.h"
#include "FourierTransform.h"
#include "StructureFactor.h"
Francois Gygi committed
28
#include "XCOperator.h"
Francois Gygi committed
29
#include "NonLocalPotential.h"
30
#include "ConfinementPotential.h"
Francois Gygi committed
31 32 33 34 35 36 37 38 39 40

#include "Timer.h"
#include "blas.h"

#include <iostream>
#include <iomanip>
#include <algorithm> // fill()
using namespace std;

////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
41
EnergyFunctional::EnergyFunctional( Sample& s, const ChargeDensity& cd)
Francois Gygi committed
42
 : s_(s), cd_(cd)
Francois Gygi committed
43 44 45
{
  const AtomSet& atoms = s_.atoms;
  const Wavefunction& wf = s_.wf;
46

47 48 49 50 51 52 53 54
  sigma_ekin.resize(6);
  sigma_econf.resize(6);
  sigma_eps.resize(6);
  sigma_ehart.resize(6);
  sigma_exc.resize(6);
  sigma_enl.resize(6);
  sigma_esr.resize(6);
  sigma.resize(6);
55

Francois Gygi committed
56
  vbasis_ = cd_.vbasis();
57

Francois Gygi committed
58 59
  // define FT's on vbasis contexts
  vft = cd_.vft();
60 61 62
  int np0v = vft->np0();
  int np1v = vft->np1();
  int np2v = vft->np2();
63

Francois Gygi committed
64 65 66 67 68 69 70 71 72
  v_r.resize(wf.nspin());
  for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
  {
    v_r[ispin].resize(vft->np012loc());
  }
  tmp_r.resize(vft->np012loc());

  if ( s_.ctxt_.onpe0() )
  {
73 74 75 76
    cout << "  EnergyFunctional: np0v,np1v,np2v: " << np0v << " "
         << np1v << " " << np2v << endl;
    cout << "  EnergyFunctional: vft->np012(): "
         << vft->np012() << endl;
Francois Gygi committed
77
  }
78

79
  const int ngloc = vbasis_->localsize();
80

Francois Gygi committed
81
  nsp_ = atoms.nsp();
82

Francois Gygi committed
83 84 85
  vps.resize(nsp_);
  dvps.resize(nsp_);
  rhops.resize(nsp_);
86

Francois Gygi committed
87 88 89 90
  zv_.resize(nsp_);
  rcps_.resize(nsp_);
  na_.resize(nsp_);
  namax_ = 0;
91

Francois Gygi committed
92 93 94 95 96 97 98
  for ( int is = 0; is < nsp_; is++ )
  {
    vps[is].resize(ngloc);
    dvps[is].resize(ngloc);
    rhops[is].resize(ngloc);
    if ( atoms.na(is) > namax_ ) namax_ = atoms.na(is);
  }
99

Francois Gygi committed
100 101 102 103
  xco = new XCOperator(s_, cd_);

  nlp.resize(wf.nspin());
  for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
Francois Gygi committed
104
  {
Francois Gygi committed
105 106 107 108 109
    nlp[ispin].resize(wf.nkp());
    for ( int ikp = 0; ikp < wf.nkp(); ikp++ )
    {
      nlp[ispin][ikp] = new NonLocalPotential(s_.atoms, *wf.sd(ispin,ikp));
    }
Francois Gygi committed
110
  }
111

112 113 114 115 116 117 118
  vion_local_g.resize(ngloc);
  dvion_local_g.resize(ngloc);
  vlocal_g.resize(ngloc);
  vtemp.resize(ngloc);
  rhoelg.resize(ngloc);
  rhogt.resize(ngloc);
  rhopst.resize(ngloc);
119

120 121
  tau0.resize(nsp_);
  fion_esr.resize(nsp_);
Francois Gygi committed
122
  fext.resize(nsp_);
123
  ftmp.resize(3*namax_);
124

125
  eself_ = 0.0;
126

127 128 129
  for ( int is = 0; is < nsp_; is++ )
  {
    Species *s = atoms.species_list[is];
130

131 132 133
    const int na = atoms.na(is);
    tau0[is].resize(3*na);
    fion_esr[is].resize(3*na);
Francois Gygi committed
134
    fext[is].resize(3*na);
135

136 137
    eself_ += na * s->eself();
    na_[is] = na;
138

139 140 141
    zv_[is] = s->zval();
    rcps_[is] = s->rcps();
  }
142

Francois Gygi committed
143 144 145 146 147 148
  // FT for interpolation of wavefunctions on the fine grid
  ft.resize(wf.nkp());
  for ( int ikp = 0; ikp < wf.nkp(); ikp++ )
  {
    ft[ikp] = cd_.ft(ikp);
  }
149

150 151 152 153 154
  // Confinement potentials
  cfp.resize(wf.nkp());
  for ( int ikp = 0; ikp < wf.nkp(); ikp++ )
  {
    cfp[ikp] = 0;
Francois Gygi committed
155 156 157 158 159
    const double facs = 2.0;
    const double sigmas = 0.5;
    cfp[ikp] =
      new ConfinementPotential(s_.ctrl.ecuts,facs,sigmas,
        wf.sd(0,ikp)->basis());
160
  }
161

162
  sf.init(tau0,*vbasis_);
163

164
  cell_moved();
165

166
  atoms_moved();
Francois Gygi committed
167

Francois Gygi committed
168 169 170 171 172
}

////////////////////////////////////////////////////////////////////////////////
EnergyFunctional::~EnergyFunctional(void)
{
Francois Gygi committed
173
  delete xco;
Francois Gygi committed
174 175 176
  for ( int ikp = 0; ikp < s_.wf.nkp(); ikp++ )
  {
    delete cfp[ikp];
Francois Gygi committed
177 178
    for ( int ispin = 0; ispin < nlp.size(); ispin++ )
      delete nlp[ispin][ikp];
Francois Gygi committed
179
  }
180

Francois Gygi committed
181 182 183 184 185 186 187 188 189
  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);
    if ( s_.ctxt_.myproc()==0 )
    {
190 191 192 193 194
      cout << "<timing name=\""
           << setw(15) << (*i).first << "\""
           << " min=\"" << setprecision(3) << setw(9) << tmin << "\""
           << " max=\"" << setprecision(3) << setw(9) << tmax << "\"/>"
           << endl;
Francois Gygi committed
195 196 197 198
    }
  }
}

Francois Gygi committed
199
////////////////////////////////////////////////////////////////////////////////
200
void EnergyFunctional::update_vhxc(bool compute_stress)
Francois Gygi committed
201 202 203 204
{
  // called when the charge density has changed
  // update Hartree and xc potentials using the charge density cd_
  // compute Hartree and xc energies
205

Francois Gygi committed
206 207 208 209 210 211 212 213
  const Wavefunction& wf = s_.wf;
  const UnitCell& cell = wf.cell();
  const double omega = cell.volume();
  const double omega_inv = 1.0 / omega;
  const double *const g2i = vbasis_->g2i_ptr();
  const double fpi = 4.0 * M_PI;
  const int ngloc = vbasis_->localsize();
  double tsum[2];
214

Francois Gygi committed
215 216 217 218 219 220 221 222 223 224 225 226 227 228 229
  // compute total electronic density: rhoelg = rho_up + rho_dn
  if ( wf.nspin() == 1 )
  {
    for ( int ig = 0; ig < ngloc; ig++ )
    {
      rhoelg[ig] = omega_inv * cd_.rhog[0][ig];
    }
  }
  else
  {
    for ( int ig = 0; ig < ngloc; ig++ )
    {
      rhoelg[ig] = omega_inv * ( cd_.rhog[0][ig] + cd_.rhog[1][ig] );
    }
  }
230

231 232
  // update XC operator
  // compute xc energy, update self-energy operator and potential
Francois Gygi committed
233 234 235
  tmap["exc"].start();
  for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
    fill(v_r[ispin].begin(),v_r[ispin].end(),0.0);
236
  xco->update(v_r, compute_stress);
Francois Gygi committed
237
  exc_ = xco->exc();
238 239
  if ( compute_stress )
    xco->compute_stress(sigma_exc);
Francois Gygi committed
240
  tmap["exc"].stop();
241 242

  // compute local potential energy:
Francois Gygi committed
243 244 245 246 247 248 249 250 251 252
  // integral of el. charge times ionic local pot.
  int len=2*ngloc,inc1=1;
  tsum[0] = 2.0 * ddot(&len,(double*)&rhoelg[0],&inc1,
         (double*)&vion_local_g[0],&inc1);
  // remove double counting for G=0
  if ( vbasis_->context().myrow() == 0 )
  {
    tsum[0] -= real(conj(rhoelg[0])*vion_local_g[0]);
  }
  tsum[0] *= omega; // tsum[0] contains eps
253

Francois Gygi committed
254 255 256 257 258 259 260 261 262 263 264 265 266 267
  // Hartree energy
  ehart_ = 0.0;
  double ehsum = 0.0;
  for ( int ig = 0; ig < ngloc; ig++ )
  {
    const complex<double> tmp = rhoelg[ig] + rhopst[ig];
    ehsum += norm(tmp) * g2i[ig];
    rhogt[ig] = tmp;
  }
  // factor 1/2 from definition of Ehart cancels with half sum over G
  // Note: rhogt[ig] includes a factor 1/Omega
  // Factor omega in next line yields prefactor 4 pi / omega in
  tsum[1] = omega * fpi * ehsum;
  // tsum[1] contains ehart
268

Francois Gygi committed
269 270 271
  vbasis_->context().dsum(2,1,&tsum[0],2);
  eps_   = tsum[0];
  ehart_ = tsum[1];
272

Francois Gygi committed
273 274 275 276 277 278
  // compute vlocal_g = vion_local_g + vhart_g
  // where vhart_g = 4 * pi * (rhoelg + rhopst) * g2i
  for ( int ig = 0; ig < ngloc; ig++ )
  {
    vlocal_g[ig] = vion_local_g[ig] + fpi * rhogt[ig] * g2i[ig];
  }
279

Francois Gygi committed
280 281
  // FT to tmpr_r
  vft->backward(&vlocal_g[0],&tmp_r[0]);
282

Francois Gygi committed
283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303
  // add local potential in tmp_r to v_r[ispin][i]
  // v_r contains the xc potential
  const int size = tmp_r.size();
  if ( wf.nspin() == 1 )
  {
    for ( int i = 0; i < size; i++ )
    {
      v_r[0][i] += real(tmp_r[i]);
    }
  }
  else
  {
    for ( int i = 0; i < size; i++ )
    {
      const double vloc = real(tmp_r[i]);
      v_r[0][i] += vloc;
      v_r[1][i] += vloc;
    }
  }
}

Francois Gygi committed
304 305 306
////////////////////////////////////////////////////////////////////////////////
double EnergyFunctional::energy(bool compute_hpsi, Wavefunction& dwf,
              bool compute_forces, vector<vector<double> >& fion,
307
              bool compute_stress, valarray<double>& sigma)
Francois Gygi committed
308
{
309
  const bool debug_stress = compute_stress &&
310
    s_.ctrl.debug.find("STRESS") != string::npos;
Francois Gygi committed
311 312 313 314 315 316 317
  const double fpi = 4.0 * M_PI;

  const Wavefunction& wf = s_.wf;
  const UnitCell& cell(wf.cell());
  const double omega = cell.volume();
  const double omega_inv = 1.0 / omega;
  const int ngloc = vbasis_->localsize();
Francois Gygi committed
318
  const double *const g2i = vbasis_->g2i_ptr();
319

320
  const bool use_confinement = s_.ctrl.ecuts > 0.0;
321

322
  for ( int is = 0; is < fion.size(); is++ )
Francois Gygi committed
323 324
    for ( int i = 0; i < fion[is].size(); i++ )
      fion[is][i] = 0.0;
325

Francois Gygi committed
326 327 328 329 330 331
  if ( compute_hpsi )
  {
    for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
      for ( int ikp = 0; ikp < wf.nkp(); ikp++ )
        dwf.sd(ispin,ikp)->c().clear();
  }
332

Francois Gygi committed
333 334
  // kinetic energy
  tmap["ekin"].start();
335

336 337 338 339 340
  // compute ekin, confinement energy, stress from ekin and econf
  // ekin = sum_G |c_G|^2  G^2
  // econf = sum_G |c_G|^2 fstress[G]
  // stress_ekin_ij = (1/Omega) sum_G |c_G|^2 * 2 * G_i * G_j
  // stress_econf_ij = (1/Omega) sum_G |c_G|^2 * dfstress[G] * G_i * G_j
Francois Gygi committed
341
  ekin_ = 0.0;
342 343 344 345
  econf_ = 0.0;
  sigma_ekin = 0.0;
  sigma_econf = 0.0;
  valarray<double> sum(0.0,14), tsum(0.0,14);
Francois Gygi committed
346 347 348 349
  for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
  {
    for ( int ikp = 0; ikp < wf.nkp(); ikp++ )
    {
350
      const double weight = wf.weight(ikp);
Francois Gygi committed
351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369
      const SlaterDet& sd = *(wf.sd(ispin,ikp));
      const Basis& wfbasis = sd.basis();
      const D3vector kp = wfbasis.kpoint();
      // factor fac in next lines: 2.0 for G and -G (if basis is real) and
      // 0.5 from 1/(2m)
      const double fac = wfbasis.real() ? 1.0 : 0.5;
      const ComplexMatrix& c = sd.c();
      const Context& sdctxt = sd.context();

      // compute psi2sum(G) = fac * sum_G occ(n) psi2(n,G)
      const int ngwloc = wfbasis.localsize();
      valarray<double> psi2sum(ngwloc);
      const complex<double>* p = c.cvalptr();
      const int mloc = c.mloc();
      const int nloc = c.nloc();
      // nn = global n index
      const int nnbase = sdctxt.mycol() * c.nb();
      const double * const occ = sd.occ_ptr(nnbase);
      for ( int ig = 0; ig < ngwloc; ig++ )
370
      {
Francois Gygi committed
371 372
        double tmpsum = 0.0;
        for ( int n = 0; n < nloc; n++ )
373
        {
Francois Gygi committed
374 375
          const double psi2 = norm(p[ig+n*mloc]);
          tmpsum += fac * occ[n] * psi2;
376
        }
Francois Gygi committed
377 378
        psi2sum[ig] = tmpsum;
      }
379

Francois Gygi committed
380 381 382 383 384 385
      // accumulate contributions to ekin,econf,sigma_ekin,sigma_econf in tsum
      const double *const kpg2  = wfbasis.kpg2_ptr();
      const double *const kpg_x = wfbasis.kpgx_ptr(0);
      const double *const kpg_y = wfbasis.kpgx_ptr(1);
      const double *const kpg_z = wfbasis.kpgx_ptr(2);
      tsum = 0.0;
386

Francois Gygi committed
387 388 389
      for ( int ig = 0; ig < ngwloc; ig++ )
      {
        const double psi2s = psi2sum[ig];
390

Francois Gygi committed
391 392
        // tsum[0]: ekin partial sum
        tsum[0] += psi2s * kpg2[ig];
393

Francois Gygi committed
394 395 396 397 398
        if ( compute_stress )
        {
          const double tkpgx = kpg_x[ig];
          const double tkpgy = kpg_y[ig];
          const double tkpgz = kpg_z[ig];
399

Francois Gygi committed
400
          const double fac_ekin = 2.0 * psi2s;
401

Francois Gygi committed
402 403 404 405 406 407
          tsum[1]  += fac_ekin * tkpgx * tkpgx;
          tsum[2]  += fac_ekin * tkpgy * tkpgy;
          tsum[3]  += fac_ekin * tkpgz * tkpgz;
          tsum[4]  += fac_ekin * tkpgx * tkpgy;
          tsum[5]  += fac_ekin * tkpgy * tkpgz;
          tsum[6]  += fac_ekin * tkpgx * tkpgz;
408

Francois Gygi committed
409 410 411 412
        }
        // tsum[0-6] contains the contributions to
        // ekin, sigma_ekin, from vector ig
      } // ig
413

Francois Gygi committed
414 415 416 417 418
      if ( use_confinement )
      {
        const valarray<double>& fstress = cfp[ikp]->fstress();
        const valarray<double>& dfstress = cfp[ikp]->dfstress();
        for ( int ig = 0; ig < ngwloc; ig++ )
419
        {
Francois Gygi committed
420 421 422
          const double psi2s = psi2sum[ig];
          // tsum[7]: econf partial sum
          tsum[7] += psi2s * fstress[ig];
423

Francois Gygi committed
424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440
          if ( compute_stress )
          {
            const double tkpgx = kpg_x[ig];
            const double tkpgy = kpg_y[ig];
            const double tkpgz = kpg_z[ig];

            const double fac_econf = psi2s * dfstress[ig];
            tsum[8]  += fac_econf * tkpgx * tkpgx;
            tsum[9]  += fac_econf * tkpgy * tkpgy;
            tsum[10] += fac_econf * tkpgz * tkpgz;
            tsum[11] += fac_econf * tkpgx * tkpgy;
            tsum[12] += fac_econf * tkpgy * tkpgz;
            tsum[13] += fac_econf * tkpgx * tkpgz;
          }
          // tsum[7-13] contains the contributions to
          // econf,sigma_econf from vector ig
        } // ig
441
      }
Francois Gygi committed
442 443

      sum += weight * tsum;
444 445
    } // ikp
  } // ispin
446

447
  // sum contains the contributions to ekin, etc.. from this task
Francois Gygi committed
448
  wf.context().dsum(14,1,&sum[0],14);
449

450
  ekin_  = sum[0];
451 452 453 454 455 456
  sigma_ekin[0] = sum[1];
  sigma_ekin[1] = sum[2];
  sigma_ekin[2] = sum[3];
  sigma_ekin[3] = sum[4];
  sigma_ekin[4] = sum[5];
  sigma_ekin[5] = sum[6];
457

458
  econf_ = sum[7];
459 460 461 462 463 464
  sigma_econf[0] = sum[8];
  sigma_econf[1] = sum[9];
  sigma_econf[2] = sum[10];
  sigma_econf[3] = sum[11];
  sigma_econf[4] = sum[12];
  sigma_econf[5] = sum[13];
465

466 467
  sigma_ekin *= omega_inv;
  sigma_econf *= omega_inv;
468

Francois Gygi committed
469
  tmap["ekin"].stop();
470

471 472 473 474
  // Stress from Eps
  sigma_eps = 0.0;
  if ( compute_stress )
  {
Francois Gygi committed
475
    tsum = 0.0;
476 477 478 479 480 481 482 483
    const double *const gi  = vbasis_->gi_ptr();
    const double *const g_x = vbasis_->gx_ptr(0);
    const double *const g_y = vbasis_->gx_ptr(1);
    const double *const g_z = vbasis_->gx_ptr(2);
    for ( int ig = 0; ig < ngloc; ig++ )
    {
      // factor of 2 in next line: G and -G
      // note: gi[0] == 0.0 in next line (no division by zero)
484
      const double fac = 2.0 * gi[ig] *
485
        real( conj(rhoelg[ig]) * dvion_local_g[ig] );
486

487 488 489
      const double tgx = g_x[ig];
      const double tgy = g_y[ig];
      const double tgz = g_z[ig];
490

Francois Gygi committed
491 492 493 494 495 496
      tsum[0] += fac * tgx * tgx;
      tsum[1] += fac * tgy * tgy;
      tsum[2] += fac * tgz * tgz;
      tsum[3] += fac * tgx * tgy;
      tsum[4] += fac * tgy * tgz;
      tsum[5] += fac * tgx * tgz;
497
    }
Francois Gygi committed
498
    vbasis_->context().dsum(6,1,&tsum[0],6);
Francois Gygi committed
499

Francois Gygi committed
500 501 502 503 504 505
    sigma_eps[0] = eps_ * omega_inv + tsum[0];
    sigma_eps[1] = eps_ * omega_inv + tsum[1];
    sigma_eps[2] = eps_ * omega_inv + tsum[2];
    sigma_eps[3] = tsum[3];
    sigma_eps[4] = tsum[4];
    sigma_eps[5] = tsum[5];
Francois Gygi committed
506
  }
507

508 509 510
  // Stress from Hartree energy
  if ( compute_stress )
  {
Francois Gygi committed
511
    tsum = 0.0;
512 513 514
    const double *const g_x = vbasis_->gx_ptr(0);
    const double *const g_y = vbasis_->gx_ptr(1);
    const double *const g_z = vbasis_->gx_ptr(2);
515

516 517 518 519 520 521 522
    for ( int ig = 0; ig < ngloc; ig++ )
    {
      const double temp = norm(rhogt[ig]) * g2i[ig] * g2i[ig];
      const double tgx = g_x[ig];
      const double tgy = g_y[ig];
      const double tgz = g_z[ig];

Francois Gygi committed
523 524 525 526 527 528
      tsum[0] += temp * tgx * tgx;
      tsum[1] += temp * tgy * tgy;
      tsum[2] += temp * tgz * tgz;
      tsum[3] += temp * tgx * tgy;
      tsum[4] += temp * tgy * tgz;
      tsum[5] += temp * tgx * tgz;
529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545
    }

    for ( int is = 0; is < nsp_; is++ )
    {
      // Factor 2 in next line: 2*real(x) = ( x + h.c. )
      // Factor 0.5 in next line: definition of sigma
      const double fac2 = 2.0 * 0.25 * rcps_[is]*rcps_[is];
      complex<double> *s = &sf.sfac[is][0];
      for ( int ig = 0; ig < ngloc; ig++ )
      {
        const complex<double> sg = s[ig];
        const complex<double> rg = rhogt[ig];
        // next line: keep only real part
        const double temp = fac2 *
        ( rg.real() * sg.real() +
          rg.imag() * sg.imag() )
        * rhops[is][ig] * g2i[ig];
546

547 548 549 550
        const double tgx = g_x[ig];
        const double tgy = g_y[ig];
        const double tgz = g_z[ig];

Francois Gygi committed
551 552 553 554 555 556
        tsum[0] += temp * tgx * tgx;
        tsum[1] += temp * tgy * tgy;
        tsum[2] += temp * tgz * tgz;
        tsum[3] += temp * tgx * tgy;
        tsum[4] += temp * tgy * tgz;
        tsum[5] += temp * tgx * tgz;
557 558
      }
    }
Francois Gygi committed
559 560 561 562 563 564 565 566 567 568 569
    vbasis_->context().dsum(6,1,&tsum[0],6);
    // Factor in next line:
    //  factor 2 from G and -G
    //  factor fpi from definition of sigma
    //  no factor 1/Omega^2 (already included in rhogt[ig] above)
    sigma_ehart[0] = ehart_ * omega_inv - 2.0 * fpi * tsum[0];
    sigma_ehart[1] = ehart_ * omega_inv - 2.0 * fpi * tsum[1];
    sigma_ehart[2] = ehart_ * omega_inv - 2.0 * fpi * tsum[2];
    sigma_ehart[3] = - 2.0 * fpi * tsum[3];
    sigma_ehart[4] = - 2.0 * fpi * tsum[4];
    sigma_ehart[5] = - 2.0 * fpi * tsum[5];
570
  } // compute_stress
571

572
  // Non local energy and forces
Francois Gygi committed
573
  tmap["nonlocal"].start();
Francois Gygi committed
574 575
  // modify next loop to span only local ikp
  enl_ = 0.0;
576 577 578 579 580 581
  vector<vector<double> > fion_enl;
  fion_enl.resize(nsp_);
  for ( int is = 0; is < nsp_; is++ )
    fion_enl[is].resize(3*na_[is]);
  valarray<double> sigma_enl_kp(6);
  sigma_enl = 0.0;
Francois Gygi committed
582
  for ( int ikp = 0; ikp < wf.nkp(); ikp++ )
583
  {
Francois Gygi committed
584
    for ( int ispin = 0; ispin < nlp.size(); ispin++ )
585
    {
Francois Gygi committed
586 587 588
      enl_ += wf.weight(ikp) * nlp[ispin][ikp]->energy(compute_hpsi,
              *dwf.sd(ispin,ikp), compute_forces, fion_enl, compute_stress,
              sigma_enl_kp);
589

590 591 592 593
      if ( compute_forces )
        for ( int is = 0; is < nsp_; is++ )
          for ( int i = 0; i < 3*na_[is]; i++ )
            fion[is][i] += wf.weight(ikp) * fion_enl[is][i];
594

595 596 597
      if ( compute_stress )
        sigma_enl += wf.weight(ikp) * sigma_enl_kp;
    }
598
  }
Francois Gygi committed
599 600 601
  tmap["nonlocal"].stop();

  ecoul_ = ehart_ + esr_ - eself_;
602 603 604 605 606 607 608
  ets_ = 0.0;
  if ( s_.ctrl.fermi_temp > 0.0 )
  {
    const double wf_entropy = wf.entropy();
    const double boltz = 1.0 / ( 11605.0 * 2.0 * 13.6058 );
    ets_ = - wf_entropy * s_.ctrl.fermi_temp * boltz;
  }
609
  etotal_ = ekin_ + econf_ + eps_ + enl_ + ecoul_ + exc_ + ets_ + eexf_;
610

Francois Gygi committed
611 612 613 614 615 616 617
  if ( compute_hpsi )
  {
    tmap["hpsi"].start();
    for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
    {
      for ( int ikp = 0; ikp < wf.nkp(); ikp++ )
      {
Francois Gygi committed
618 619 620 621 622 623 624 625 626 627 628
        const SlaterDet& sd = *(wf.sd(ispin,ikp));
        SlaterDet& sdp = *(dwf.sd(ispin,ikp));
        const ComplexMatrix& c = sd.c();
        const Basis& wfbasis = sd.basis();
        ComplexMatrix& cp = dwf.sd(ispin,ikp)->c();
        const int mloc = cp.mloc();
        const double* kpg2 = wfbasis.kpg2_ptr();
        const int ngwloc = wfbasis.localsize();

        // Laplacian
        if ( use_confinement )
Francois Gygi committed
629
        {
Francois Gygi committed
630
          for ( int n = 0; n < sd.nstloc(); n++ )
Francois Gygi committed
631
          {
Francois Gygi committed
632 633
            const valarray<double>& fstress = cfp[ikp]->fstress();
            for ( int ig = 0; ig < ngwloc; ig++ )
634
            {
Francois Gygi committed
635 636
              cp[ig+mloc*n] += 0.5 * ( kpg2[ig] + fstress[ig] ) *
                               c[ig+mloc*n];
637
            }
638
          }
Francois Gygi committed
639 640 641 642
        }
        else
        {
          for ( int n = 0; n < sd.nstloc(); n++ )
643
          {
Francois Gygi committed
644
            for ( int ig = 0; ig < ngwloc; ig++ )
Francois Gygi committed
645
            {
Francois Gygi committed
646
              cp[ig+mloc*n] += 0.5 * kpg2[ig] * c[ig+mloc*n];
Francois Gygi committed
647 648
            }
          }
Francois Gygi committed
649
        }
Francois Gygi committed
650 651

        // local potential
Francois Gygi committed
652
        sd.rs_mul_add(*ft[ikp], &v_r[ispin][0], sdp);
Francois Gygi committed
653 654 655 656
      } //ikp
    } //ispin

    // apply self-energy operator
657
    xco->apply_self_energy(dwf);
Francois Gygi committed
658

Francois Gygi committed
659 660
    tmap["hpsi"].stop();
  } // if compute_hpsi
661

Francois Gygi committed
662 663 664 665 666 667
  if ( compute_forces )
  {
    const int* idx = vbasis_->idx_ptr();
    const double* gx0 = vbasis_->gx_ptr(0);
    const double* gx1 = vbasis_->gx_ptr(1);
    const double* gx2 = vbasis_->gx_ptr(2);
668

Francois Gygi committed
669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689
    for ( int is = 0; is < nsp_; is++ )
    {
      for ( int ig = 0; ig < ngloc; ig++ )
      {
        double tmp = fpi * rhops[is][ig] * g2i[ig];
        vtemp[ig] =  tmp * conj(rhogt[ig]) + vps[is][ig] * conj(rhoelg[ig]);
      }
      memset((void*)&ftmp[0],0,3*namax_*sizeof(double));
      // loop over atoms of species is
      for ( int ia = 0; ia < na_[is]; ia++ )
      {
        double sum0=0.0,sum1=0.0,sum2=0.0;
        double *c0 = sf.cos0_ptr(is,ia);
        double *c1 = sf.cos1_ptr(is,ia);
        double *c2 = sf.cos2_ptr(is,ia);
        double *s0 = sf.sin0_ptr(is,ia);
        double *s1 = sf.sin1_ptr(is,ia);
        double *s2 = sf.sin2_ptr(is,ia);
        for ( int ig = 0; ig < ngloc; ig++ )
        {
          // compute Exp[igr] in 3D as a product of 1D contributions
Francois Gygi committed
690 691 692
          // complex<double> teigr = ei0[kv[3*ig+0]] *
          //                         ei1[kv[3*ig+1]] *
          //                         ei2[kv[3*ig+2]];
Francois Gygi committed
693 694 695 696
          const int iii = ig+ig+ig;
          const int kx = idx[iii];
          const int ky = idx[iii+1];
          const int kz = idx[iii+2];
697

Francois Gygi committed
698 699 700
          const double cos_a = c0[kx];
          const double cos_b = c1[ky];
          const double cos_c = c2[kz];
701

Francois Gygi committed
702 703 704
          const double sin_a = s0[kx];
          const double sin_b = s1[ky];
          const double sin_c = s2[kz];
705

Francois Gygi committed
706 707
          // Next line: exp(-i*gr) =
          // (cos_a - I sin_a)*(cos_b - I sin_b)*(cos_c - I sin_c)
708
          double teigr_re =
Francois Gygi committed
709 710
            cos_a*cos_b*cos_c - sin_a*sin_b*cos_c -
            sin_a*cos_b*sin_c - cos_a*sin_b*sin_c;
711
          double teigr_im =
Francois Gygi committed
712 713
            sin_a*sin_b*sin_c - sin_a*cos_b*cos_c -
            cos_a*sin_b*cos_c - cos_a*cos_b*sin_c;
714

Francois Gygi committed
715
          /* fion is real */
716
          double tmp = teigr_re * vtemp[ig].imag() +
Francois Gygi committed
717 718 719 720 721 722 723 724 725 726 727
                       teigr_im * vtemp[ig].real();

          sum0 += tmp * gx0[ig];
          sum1 += tmp * gx1[ig];
          sum2 += tmp * gx2[ig];
        }
        ftmp[3*ia]   = sum0;
        ftmp[3*ia+1] = sum1;
        ftmp[3*ia+2] = sum2;

      }
728

Francois Gygi committed
729 730 731 732 733 734 735 736
      int len = 3*na_[is];
      vbasis_->context().dsum(len,1,&ftmp[0],len);

      double fac = -2.0 * omega;
      for ( int i = 0; i < 3*na_[is]; i++ )
      {
        fion[is][i] += fion_esr[is][i] + fac * ftmp[i];
      }
Francois Gygi committed
737 738 739 740 741 742 743 744 745

      // add external forces
      if ( s_.extforces.size() > 0 )
      {
        assert(fion.size()==fext.size());
        assert(fion[is].size()==fext[is].size());
        for ( int i = 0; i < fion[is].size(); i++ )
          fion[is][i] += fext[is][i];
      }
Francois Gygi committed
746 747
    }
  }
748

749 750 751 752
  if ( compute_stress )
    sigma = sigma_ekin + sigma_econf + sigma_eps + sigma_enl +
            sigma_ehart + sigma_exc + sigma_esr;

753
  if ( debug_stress && s_.ctxt_.onpe0() )
754
  {
Francois Gygi committed
755
    //const double gpa = 29421.5;
756 757 758
    cout.setf(ios::fixed,ios::floatfield);
    cout.setf(ios::right,ios::adjustfield);
    cout << setprecision(8);
759
    cout << " <stress_tensor unit=\"atomic_units\">\n"
760 761 762 763 764 765 766 767 768 769 770 771
         << "   <sigma_ekin_xx> " << setw(12)
         << sigma_ekin[0] << " </sigma_ekin_xx>\n"
         << "   <sigma_ekin_yy> " << setw(12)
         << sigma_ekin[1] << " </sigma_ekin_yy>\n"
         << "   <sigma_ekin_zz> " << setw(12)
         << sigma_ekin[2] << " </sigma_ekin_zz>\n"
         << "   <sigma_ekin_xy> " << setw(12)
         << sigma_ekin[3] << " </sigma_ekin_xy>\n"
         << "   <sigma_ekin_yz> " << setw(12)
         << sigma_ekin[4] << " </sigma_ekin_yz>\n"
         << "   <sigma_ekin_xz> " << setw(12)
         << sigma_ekin[5] << " </sigma_ekin_xz>\n"
772
         << endl
773 774 775 776 777 778 779 780 781 782 783 784
         << "   <sigma_econf_xx> " << setw(12)
         << sigma_econf[0] << " </sigma_econf_xx>\n"
         << "   <sigma_econf_yy> " << setw(12)
         << sigma_econf[1] << " </sigma_econf_yy>\n"
         << "   <sigma_econf_zz> " << setw(12)
         << sigma_econf[2] << " </sigma_econf_zz>\n"
         << "   <sigma_econf_xy> " << setw(12)
         << sigma_econf[3] << " </sigma_econf_xy>\n"
         << "   <sigma_econf_yz> " << setw(12)
         << sigma_econf[4] << " </sigma_econf_yz>\n"
         << "   <sigma_econf_xz> " << setw(12)
         << sigma_econf[5] << " </sigma_econf_xz>\n"
785
         << endl
786 787 788 789 790 791 792 793 794 795 796 797
         << "   <sigma_eps_xx> " << setw(12)
         << sigma_eps[0] << " </sigma_eps_xx>\n"
         << "   <sigma_eps_yy> " << setw(12)
         << sigma_eps[1] << " </sigma_eps_yy>\n"
         << "   <sigma_eps_zz> " << setw(12)
         << sigma_eps[2] << " </sigma_eps_zz>\n"
         << "   <sigma_eps_xy> " << setw(12)
         << sigma_eps[3] << " </sigma_eps_xy>\n"
         << "   <sigma_eps_yz> " << setw(12)
         << sigma_eps[4] << " </sigma_eps_yz>\n"
         << "   <sigma_eps_xz> " << setw(12)
         << sigma_eps[5] << " </sigma_eps_xz>\n"
798
         << endl
799 800 801 802 803 804 805 806 807 808 809 810
         << "   <sigma_enl_xx> " << setw(12)
         << sigma_enl[0] << " </sigma_enl_xx>\n"
         << "   <sigma_enl_yy> " << setw(12)
         << sigma_enl[1] << " </sigma_enl_yy>\n"
         << "   <sigma_enl_zz> " << setw(12)
         << sigma_enl[2] << " </sigma_enl_zz>\n"
         << "   <sigma_enl_xy> " << setw(12)
         << sigma_enl[3] << " </sigma_enl_xy>\n"
         << "   <sigma_enl_yz> " << setw(12)
         << sigma_enl[4] << " </sigma_enl_yz>\n"
         << "   <sigma_enl_xz> " << setw(12)
         << sigma_enl[5] << " </sigma_enl_xz>\n"
811
         << endl
812 813 814 815 816 817 818 819 820 821 822 823
         << "   <sigma_ehart_xx> " << setw(12)
         << sigma_ehart[0] << " </sigma_ehart_xx>\n"
         << "   <sigma_ehart_yy> " << setw(12)
         << sigma_ehart[1] << " </sigma_ehart_yy>\n"
         << "   <sigma_ehart_zz> " << setw(12)
         << sigma_ehart[2] << " </sigma_ehart_zz>\n"
         << "   <sigma_ehart_xy> " << setw(12)
         << sigma_ehart[3] << " </sigma_ehart_xy>\n"
         << "   <sigma_ehart_yz> " << setw(12)
         << sigma_ehart[4] << " </sigma_ehart_yz>\n"
         << "   <sigma_ehart_xz> " << setw(12)
         << sigma_ehart[5] << " </sigma_ehart_xz>\n"
824
         << endl
825 826 827 828 829 830 831 832 833 834 835 836
         << "   <sigma_exc_xx> " << setw(12)
         << sigma_exc[0] << " </sigma_exc_xx>\n"
         << "   <sigma_exc_yy> " << setw(12)
         << sigma_exc[1] << " </sigma_exc_yy>\n"
         << "   <sigma_exc_zz> " << setw(12)
         << sigma_exc[2] << " </sigma_exc_zz>\n"
         << "   <sigma_exc_xy> " << setw(12)
         << sigma_exc[3] << " </sigma_exc_xy>\n"
         << "   <sigma_exc_yz> " << setw(12)
         << sigma_exc[4] << " </sigma_exc_yz>\n"
         << "   <sigma_exc_xz> " << setw(12)
         << sigma_exc[5] << " </sigma_exc_xz>\n"
837
         << endl
838 839 840 841 842 843 844 845 846 847 848 849
         << "   <sigma_esr_xx> " << setw(12)
         << sigma_esr[0] << " </sigma_esr_xx>\n"
         << "   <sigma_esr_yy> " << setw(12)
         << sigma_esr[1] << " </sigma_esr_yy>\n"
         << "   <sigma_esr_zz> " << setw(12)
         << sigma_esr[2] << " </sigma_esr_zz>\n"
         << "   <sigma_esr_xy> " << setw(12)
         << sigma_esr[3] << " </sigma_esr_xy>\n"
         << "   <sigma_esr_yz> " << setw(12)
         << sigma_esr[4] << " </sigma_esr_yz>\n"
         << "   <sigma_esr_xz> " << setw(12)
         << sigma_esr[5] << " </sigma_esr_xz>\n"
850
         << endl
851 852 853 854 855 856
         << "   <sigma_eks_xx> " << setw(12) << sigma[0] << " </sigma_eks_xx>\n"
         << "   <sigma_eks_yy> " << setw(12) << sigma[1] << " </sigma_eks_yy>\n"
         << "   <sigma_eks_zz> " << setw(12) << sigma[2] << " </sigma_eks_zz>\n"
         << "   <sigma_eks_xy> " << setw(12) << sigma[3] << " </sigma_eks_xy>\n"
         << "   <sigma_eks_yz> " << setw(12) << sigma[4] << " </sigma_eks_yz>\n"
         << "   <sigma_eks_xz> " << setw(12) << sigma[5] << " </sigma_eks_xz>\n"
857 858
         << " </stress_tensor>" << endl;
  }
Francois Gygi committed
859 860 861 862 863 864 865 866 867
  return etotal_;
}

////////////////////////////////////////////////////////////////////////////////
void EnergyFunctional::atoms_moved(void)
{
  const AtomSet& atoms = s_.atoms;
  int ngloc = vbasis_->localsize();

868
  // fill tau0 with values in atom_list
869

Francois Gygi committed
870 871
  atoms.get_positions(tau0);
  sf.update(tau0,*vbasis_);
872

Francois Gygi committed
873 874
  // compute Fourier coefficients of the local potential
  memset( (void*)&vion_local_g[0], 0, 2*ngloc*sizeof(double) );
875
  memset( (void*)&dvion_local_g[0], 0, 2*ngloc*sizeof(double) );
Francois Gygi committed
876 877 878 879 880 881
  memset( (void*)&rhopst[0], 0, 2*ngloc*sizeof(double) );
  for ( int is = 0; is < atoms.nsp(); is++ )
  {
    complex<double> *s = &sf.sfac[is][0];
    for ( int ig = 0; ig < ngloc; ig++ )
    {
882 883 884 885
      const complex<double> sg = s[ig];
      rhopst[ig] += sg * rhops[is][ig];
      vion_local_g[ig] += sg * vps[is][ig];
      dvion_local_g[ig] += sg * dvps[is][ig];
Francois Gygi committed
886 887
    }
  }
888

Francois Gygi committed
889 890
  // compute esr: pseudocharge repulsion energy
  const UnitCell& cell = s_.wf.cell();
891
  const double omega_inv = 1.0 / cell.volume();
Francois Gygi committed
892 893 894 895 896 897 898 899 900 901 902 903
  // distance between opposite planes of the cell
  const double d0 = 2.0 * M_PI / length(cell.b(0));
  const double d1 = 2.0 * M_PI / length(cell.b(1));
  const double d2 = 2.0 * M_PI / length(cell.b(2));
  assert (d0!=0.0);
  assert (d1!=0.0);
  assert (d2!=0.0);
#ifdef DEBUG
  cout << " EnergyFunctional: d0 = " << d0 << endl;
  cout << " EnergyFunctional: d1 = " << d1 << endl;
  cout << " EnergyFunctional: d2 = " << d2 << endl;
#endif
904

Francois Gygi committed
905
  esr_  = 0.0;
906
  sigma_esr = 0.0;
Francois Gygi committed
907 908 909 910
  for ( int is = 0; is < nsp_; is++ )
    for ( int i = 0; i < fion_esr[is].size(); i++ )
      fion_esr[is][i] = 0.0;

Francois Gygi committed
911
  for ( int is1 = 0; is1 < nsp_; is1++ )
Francois Gygi committed
912
  {
Francois Gygi committed
913
    for ( int is2 = is1; is2 < nsp_; is2++ )
Francois Gygi committed
914
    {
Francois Gygi committed
915 916 917
      double rcps12 = sqrt ( rcps_[is1]*rcps_[is1]+rcps_[is2]*rcps_[is2] );
      // convergence criterion for lattice sums:
      // fac * rcps12 < ncell * d
918
      const double fac = 8.0;
Francois Gygi committed
919 920 921 922 923 924 925 926 927 928
      const int ncell0 = (int) (fac * rcps12 / d0);
      const int ncell1 = (int) (fac * rcps12 / d1);
      const int ncell2 = (int) (fac * rcps12 / d2);
#ifdef DEBUG
      const double mindist = min(min(d0,d1),d2);
      cout << " EnergyFunctional: mindist/rcps12: " << mindist/rcps12 << endl;
      cout << " EnergyFunctional: ncell[012]: "
           << ncell0 << " " << ncell1 << " " << ncell2 << endl;
#endif
      for ( int ia1 = 0; ia1 < na_[is1]; ia1++ )
Francois Gygi committed
929
      {
Francois Gygi committed
930
        int ia2min = 0;
931
        if ( is1 == is2 ) ia2min = ia1;
Francois Gygi committed
932
        for ( int ia2 = ia2min; ia2 < na_[is2]; ia2++ )
Francois Gygi committed
933
        {
934
          const bool same_atom = ( is1==is2 && ia1==ia2 );
Francois Gygi committed
935 936 937 938 939 940 941 942 943 944 945
          double x12 = tau0[is1][3*ia1+0] - tau0[is2][3*ia2+0];
          double y12 = tau0[is1][3*ia1+1] - tau0[is2][3*ia2+1];
          double z12 = tau0[is1][3*ia1+2] - tau0[is2][3*ia2+2];
          D3vector v12(x12,y12,z12);
          cell.fold_in_ws(v12);

          // loop over neighboring cells
          for ( int ic0 = -ncell0; ic0 <= ncell0; ic0++ )
            for ( int ic1 = -ncell1; ic1 <= ncell1; ic1++ )
              for ( int ic2 = -ncell2; ic2 <= ncell2; ic2++ )
              {
946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979
                if ( !same_atom || ic0!=0 || ic1!=0 || ic2!=0 )
                {
                  double fac = 1.0;
                  if ( same_atom )
                    fac = 0.5;
                  D3vector s = ic0*cell.a(0) + ic1*cell.a(1) + ic2*cell.a(2);
                  x12 = v12.x + s.x;
                  y12 = v12.y + s.y;
                  z12 = v12.z + s.z;

                  double r12 = sqrt(x12*x12 + y12*y12 + z12*z12);
                  double arg = r12 / rcps12;
                  double esr_term = zv_[is1] * zv_[is2] * erfc(arg) / r12;
                  esr_ += fac * esr_term;

                  double desr_erfc = 2.0 * zv_[is1]*zv_[is2] *
                         exp(-arg*arg)/(rcps12*sqrt(M_PI));

                  // desrdr = (1/r) d Esr / dr
                  double desrdr = - fac * (esr_term+desr_erfc) / ( r12*r12 );
                  fion_esr[is1][3*ia1+0] -= desrdr * x12;
                  fion_esr[is2][3*ia2+0] += desrdr * x12;
                  fion_esr[is1][3*ia1+1] -= desrdr * y12;
                  fion_esr[is2][3*ia2+1] += desrdr * y12;
                  fion_esr[is1][3*ia1+2] -= desrdr * z12;
                  fion_esr[is2][3*ia2+2] += desrdr * z12;

                  sigma_esr[0] += desrdr * x12 * x12;
                  sigma_esr[1] += desrdr * y12 * y12;
                  sigma_esr[2] += desrdr * z12 * z12;
                  sigma_esr[3] += desrdr * x12 * y12;
                  sigma_esr[4] += desrdr * y12 * z12;
                  sigma_esr[5] += desrdr * x12 * z12;
                }
Francois Gygi committed
980
              }
Francois Gygi committed
981 982 983 984
        }
      }
    }
  }
985
  sigma_esr *= - omega_inv;
Francois Gygi committed
986 987 988 989

  // get external forces in fext
  eexf_ = s_.extforces.energy(tau0,fext);

Francois Gygi committed
990 991 992 993 994
}

////////////////////////////////////////////////////////////////////////////////
void EnergyFunctional::cell_moved(void)
{
Francois Gygi committed
995 996
  const Wavefunction& wf = s_.wf;
  const UnitCell& cell = wf.cell();
997
  // resize vbasis_
Francois Gygi committed
998
  vbasis_->resize(cell,s_.wf.refcell(),4.0*s_.wf.ecut());
999

1000 1001 1002
  const int ngloc = vbasis_->localsize();
  const double omega = cell.volume();
  assert(omega != 0.0);
1003 1004
  const double omega_inv = 1.0 / omega;

1005 1006
  const AtomSet& atoms = s_.atoms;
  for ( int is = 0; is < nsp_; is++ )
1007
  {
1008 1009
    Species *s = atoms.species_list[is];
    const double * const g = vbasis_->g_ptr();
1010
    double v,dv;
1011
    for ( int ig = 0; ig < ngloc; ig++ )
1012
    {
1013 1014 1015 1016
      rhops[is][ig] = s->rhopsg(g[ig]) * omega_inv;
      s->dvlocg(g[ig],v,dv);
      vps[is][ig] =  v * omega_inv;
      dvps[is][ig] =  dv * omega_inv;
1017
    }
1018
  }
1019

Francois Gygi committed
1020
  // Update confinement potentials and non-local potentials
1021 1022
  for ( int ikp = 0; ikp < wf.nkp(); ikp++ )
  {
Francois Gygi committed
1023
    cfp[ikp]->update();
Francois Gygi committed
1024 1025
    for ( int ispin = 0; ispin < nlp.size(); ispin++ )
      nlp[ispin][ikp]->update_twnl();
1026
  }
Francois Gygi committed
1027
}
1028 1029 1030 1031 1032 1033 1034

////////////////////////////////////////////////////////////////////////////////
void EnergyFunctional::print(ostream& os) const
{
  os.setf(ios::fixed,ios::floatfield);
  os.setf(ios::right,ios::adjustfield);
  os << setprecision(8);
Francois Gygi committed
1035 1036 1037 1038 1039 1040 1041 1042 1043 1044
  os << "  <ekin>   " << setw(15) << ekin()   << " </ekin>\n"
     << "  <econf>  " << setw(15) << econf()  << " </econf>\n"
     << "  <eps>    " << setw(15) << eps()    << " </eps>\n"
     << "  <enl>    " << setw(15) << enl()    << " </enl>\n"
     << "  <ecoul>  " << setw(15) << ecoul()  << " </ecoul>\n"
     << "  <exc>    " << setw(15) << exc()    << " </exc>\n"
     << "  <esr>    " << setw(15) << esr()    << " </esr>\n"
     << "  <eself>  " << setw(15) << eself()  << " </eself>\n"
     << "  <ets>    " << setw(15) << ets()    << " </ets>\n"
     << "  <eexf>   " << setw(15) << eexf()   << " </eexf>\n"
1045 1046 1047
     << "  <etotal> " << setw(15) << etotal() << " </etotal>\n"
     << flush;
}
1048

1049 1050
////////////////////////////////////////////////////////////////////////////////
ostream& operator<< ( ostream& os, const EnergyFunctional& e )
1051 1052
{
  e.print(os);
1053 1054
  return os;
}