EnergyFunctional.C 38 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"
31 32
#include "D3vector.h"
#include "ElectricEnthalpy.h"
33
#include "ExternalPotential.h"
Francois Gygi committed
34 35 36 37 38 39

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

#include <iostream>
#include <iomanip>
40 41
#include <algorithm> // fill(), copy()
#include <functional>
Francois Gygi committed
42 43 44
using namespace std;

////////////////////////////////////////////////////////////////////////////////
45
EnergyFunctional::EnergyFunctional(Sample& s, ChargeDensity& cd)
Francois Gygi committed
46
 : s_(s), cd_(cd)
Francois Gygi committed
47 48 49
{
  const AtomSet& atoms = s_.atoms;
  const Wavefunction& wf = s_.wf;
50

51 52 53 54 55 56 57 58
  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);
59

Francois Gygi committed
60
  vbasis_ = cd_.vbasis();
61

Francois Gygi committed
62 63
  // define FT's on vbasis contexts
  vft = cd_.vft();
64 65 66
  int np0v = vft->np0();
  int np1v = vft->np1();
  int np2v = vft->np2();
67

Francois Gygi committed
68
  v_r.resize(wf.nspin());
69
  vxc_r.resize(wf.nspin());
Francois Gygi committed
70 71 72
  for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
  {
    v_r[ispin].resize(vft->np012loc());
73
    vxc_r[ispin].resize(vft->np012loc());
Francois Gygi committed
74 75 76 77 78
  }
  tmp_r.resize(vft->np012loc());

  if ( s_.ctxt_.onpe0() )
  {
79 80 81 82
    cout << " EnergyFunctional: <np0v> " << np0v << " </np0v>  "
         << "<np1v> " << np1v << " </np1v>  "
         << "<np2v> " << np2v << " </np2v>" << endl;
    cout << " EnergyFunctional: vft->np012(): "
83
         << vft->np012() << endl;
Francois Gygi committed
84
  }
85

86
  // external potential
87
  if ( s_.vext )
88 89
    s_.vext->update(cd_);

90
  const int ngloc = vbasis_->localsize();
91

Francois Gygi committed
92
  nsp_ = atoms.nsp();
93

Francois Gygi committed
94 95 96
  vps.resize(nsp_);
  dvps.resize(nsp_);
  rhops.resize(nsp_);
97

Francois Gygi committed
98 99 100
  zv_.resize(nsp_);
  rcps_.resize(nsp_);
  na_.resize(nsp_);
101

Francois Gygi committed
102 103 104 105 106 107
  for ( int is = 0; is < nsp_; is++ )
  {
    vps[is].resize(ngloc);
    dvps[is].resize(ngloc);
    rhops[is].resize(ngloc);
  }
108

Francois Gygi committed
109 110
  nlp.resize(wf.nspin());
  for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
Francois Gygi committed
111
  {
Francois Gygi committed
112 113 114 115 116
    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
117
  }
118

119 120 121 122 123 124 125
  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);
126

127 128
  tau0.resize(nsp_);
  fion_esr.resize(nsp_);
Francois Gygi committed
129
  fext.resize(nsp_);
130

131
  eself_ = 0.0;
132

133 134
  // core_charge_: true if any species has a core charge density
  core_charge_ = false;
135 136 137
  for ( int is = 0; is < nsp_; is++ )
  {
    Species *s = atoms.species_list[is];
138

139 140 141
    const int na = atoms.na(is);
    tau0[is].resize(3*na);
    fion_esr[is].resize(3*na);
Francois Gygi committed
142
    fext[is].resize(3*na);
143

144 145
    eself_ += na * s->eself();
    na_[is] = na;
146

147
    zv_[is] = s->ztot();
148
    rcps_[is] = s->rcps();
149 150 151 152 153 154 155 156 157 158

    core_charge_ |= s->has_nlcc();
  }

  if ( core_charge_ )
  {
    vxc_g.resize(ngloc);
    rhocore_sp_g.resize(nsp_);
    for ( int is = 0; is < nsp_; is++ )
      rhocore_sp_g[is].resize(ngloc);
159 160
    cd_.rhocore_g.resize(ngloc);
    cd_.rhocore_r.resize(vft->np012loc());
161
  }
162

Francois Gygi committed
163 164 165 166 167 168
  // 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);
  }
169

170 171 172 173 174
  // Confinement potentials
  cfp.resize(wf.nkp());
  for ( int ikp = 0; ikp < wf.nkp(); ikp++ )
  {
    cfp[ikp] = 0;
Francois Gygi committed
175 176 177 178 179
    const double facs = 2.0;
    const double sigmas = 0.5;
    cfp[ikp] =
      new ConfinementPotential(s_.ctrl.ecuts,facs,sigmas,
        wf.sd(0,ikp)->basis());
180
  }
181

182 183
  // Electric enthalpy
  el_enth_ = 0;
184
  if ( s_.ctrl.polarization != "OFF" )
185 186
    el_enth_ = new ElectricEnthalpy(s_);

187
  sf.init(tau0,*vbasis_);
188
  xco = new XCOperator(s_, cd_);
189 190
  cell_moved();
  atoms_moved();
Francois Gygi committed
191 192 193 194 195
}

////////////////////////////////////////////////////////////////////////////////
EnergyFunctional::~EnergyFunctional(void)
{
196
  delete el_enth_;
Francois Gygi committed
197
  delete xco;
Francois Gygi committed
198 199 200
  for ( int ikp = 0; ikp < s_.wf.nkp(); ikp++ )
  {
    delete cfp[ikp];
Francois Gygi committed
201 202
    for ( int ispin = 0; ispin < nlp.size(); ispin++ )
      delete nlp[ispin][ikp];
Francois Gygi committed
203
  }
204

Francois Gygi committed
205 206 207 208 209 210 211 212 213
  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 )
    {
214 215
      string s = "name=\"" + (*i).first + "\"";
      cout << "<timing " << left << setw(22) << s
216 217
           << " min=\"" << setprecision(3) << tmin << "\""
           << " max=\"" << setprecision(3) << tmax << "\"/>"
218
           << endl;
Francois Gygi committed
219 220 221 222
    }
  }
}

Francois Gygi committed
223
////////////////////////////////////////////////////////////////////////////////
224 225
void EnergyFunctional::update_vhxc(bool compute_stress, bool update_vh,
  bool update_vxc)
Francois Gygi committed
226 227 228 229
{
  // called when the charge density has changed
  // update Hartree and xc potentials using the charge density cd_
  // compute Hartree and xc energies
230

Francois Gygi committed
231 232 233 234 235 236 237
  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();
238
  double sum[5], tsum[5];
239

Francois Gygi committed
240 241 242 243 244 245 246 247 248 249 250 251 252 253 254
  // 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] );
    }
  }
255

256 257
  // update XC operator
  // compute xc energy, update self-energy operator and potential
Francois Gygi committed
258
  tmap["exc"].start();
259

260 261
  if ( update_vxc )
    xco->update(vxc_r, compute_stress);
He Ma committed
262

Francois Gygi committed
263
  exc_ = xco->exc();
264
  dxc_ = xco->dxc();
265 266
  if ( compute_stress )
    xco->compute_stress(sigma_exc);
267 268 269 270 271 272

  if ( core_charge_ )
  {
    // compute Fourier coefficients of Vxc
    for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
    {
273
      copy(vxc_r[ispin].begin(),vxc_r[ispin].end(),tmp_r.begin());
274 275 276 277 278 279 280 281 282 283 284
      if ( ispin == 0 )
      {
        vft->forward(&tmp_r[0],&vxc_g[0]);
      }
      else
      {
        vft->forward(&tmp_r[0],&vtemp[0]);
        for ( int ig = 0; ig < ngloc; ig++ )
          vxc_g[ig] = 0.5 * ( vxc_g[ig] + vtemp[ig] );
      }
    }
285 286 287 288 289 290 291 292 293 294 295 296
    // correct dxc_ to include the core charge exchange-correlation term
    // in Harris-Foulkes estimates
    // dxc_correction = sum_ig vxc_g[ig] * rhocore_g[ig]
    double sum = 0.0;
    complex<double> *rh = &cd_.rhocore_g[0];
    int len=2*ngloc,inc1=1;
    sum = 2.0 * ddot(&len,(double*)&vxc_g[0],&inc1,(double*)&rh[0],&inc1);
    // remove double counting for G=0
    if ( vbasis_->mype() == 0 )
      sum -= real(conj(vxc_g[0])*rh[0]);
    double tsum = 0.0;
    MPI_Allreduce(&sum,&tsum,1,MPI_DOUBLE,MPI_SUM,vbasis_->comm());
297 298
    // factor nspin: rhocore_g contains half the core charge if nspin == 2
    dxc_ += wf.nspin() * tsum;
299 300
  }

Francois Gygi committed
301
  tmap["exc"].stop();
302 303

  // compute local potential energy:
Francois Gygi committed
304 305
  // integral of el. charge times ionic local pot.
  int len=2*ngloc,inc1=1;
306
  sum[0] = 2.0 * ddot(&len,(double*)&rhoelg[0],&inc1,
Francois Gygi committed
307 308
         (double*)&vion_local_g[0],&inc1);
  // remove double counting for G=0
309
  if ( vbasis_->mype() == 0 )
Francois Gygi committed
310
  {
311
    sum[0] -= real(conj(rhoelg[0])*vion_local_g[0]);
Francois Gygi committed
312
  }
313
  sum[0] *= omega; // sum[0] contains eps
314

315
  // Hartree energy and electron-electron energy (without pseudocharges)
Francois Gygi committed
316
  double ehsum = 0.0;
317 318 319
  double ehesum = 0.0;
  double ehepsum = 0.0;
  double ehpsum = 0.0;
Francois Gygi committed
320 321
  for ( int ig = 0; ig < ngloc; ig++ )
  {
322 323 324 325 326 327
    const complex<double> r = rhoelg[ig];
    const complex<double> rp = rhopst[ig];
    ehsum  += norm(r+rp) * g2i[ig];
    ehesum += norm(r) * g2i[ig];
    ehepsum += 2.0*real(conj(r)*rp * g2i[ig]);
    ehpsum += norm(rp) * g2i[ig];
328
    if ( update_vh ) rhogt[ig] = r+rp;
Francois Gygi committed
329 330 331
  }
  // factor 1/2 from definition of Ehart cancels with half sum over G
  // Note: rhogt[ig] includes a factor 1/Omega
332
  // Factor omega in next line yields prefactor 4 pi / omega in Ehart
333
  sum[1] = omega * fpi * ehsum;
334 335 336
  sum[2] = omega * fpi * ehesum;
  sum[3] = omega * fpi * ehepsum;
  sum[4] = omega * fpi * ehpsum;
337

338 339 340 341 342 343 344
  MPI_Allreduce(sum,tsum,5,MPI_DOUBLE,MPI_SUM,vbasis_->comm());

  eps_      = tsum[0];
  ehart_    = tsum[1];
  ehart_e_  = tsum[2];
  ehart_ep_ = tsum[3];
  ehart_p_  = tsum[4];
345

Francois Gygi committed
346 347 348 349 350 351
  // 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];
  }
352

353
  // FT to tmp_r
Francois Gygi committed
354
  vft->backward(&vlocal_g[0],&tmp_r[0]);
355

356
  // add external potential vext to tmp_r
357
  if ( s_.vext )
358 359 360
  {
    for ( int i = 0; i < tmp_r.size(); i++ )
      tmp_r[i] += s_.vext->v(i);
Francois Gygi committed
361 362
    // update eext_
    eext_ = s_.vext->compute_eext(cd_);
363 364
  }

365 366
  // compute local potential v_r[ispin][i]
  // vxc_r contains the xc potential
Francois Gygi committed
367 368 369 370 371
  const int size = tmp_r.size();
  if ( wf.nspin() == 1 )
  {
    for ( int i = 0; i < size; i++ )
    {
372
      v_r[0][i] = vxc_r[0][i] + real(tmp_r[i]);
Francois Gygi committed
373 374 375 376 377 378 379
    }
  }
  else
  {
    for ( int i = 0; i < size; i++ )
    {
      const double vloc = real(tmp_r[i]);
380 381
      v_r[0][i] = vxc_r[0][i] + vloc;
      v_r[1][i] = vxc_r[1][i] + vloc;
Francois Gygi committed
382 383
    }
  }
384 385 386

  if ( el_enth_ )
    el_enth_->update();
Francois Gygi committed
387 388
}

Francois Gygi committed
389 390 391
////////////////////////////////////////////////////////////////////////////////
double EnergyFunctional::energy(bool compute_hpsi, Wavefunction& dwf,
              bool compute_forces, vector<vector<double> >& fion,
392
              bool compute_stress, valarray<double>& sigma)
Francois Gygi committed
393
{
394
  const bool debug_stress = compute_stress &&
395
    s_.ctrl.debug.find("STRESS") != string::npos;
Francois Gygi committed
396 397 398 399 400 401 402
  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
403
  const double *const g2i = vbasis_->g2i_ptr();
404

405
  const bool use_confinement = s_.ctrl.ecuts > 0.0;
406

407
  for ( int is = 0; is < fion.size(); is++ )
Francois Gygi committed
408 409
    for ( int i = 0; i < fion[is].size(); i++ )
      fion[is][i] = 0.0;
410

Francois Gygi committed
411 412 413 414 415 416
  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();
  }
417

Francois Gygi committed
418 419
  // kinetic energy
  tmap["ekin"].start();
420

421 422 423 424 425
  // 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
426
  ekin_ = 0.0;
427 428 429 430
  econf_ = 0.0;
  sigma_ekin = 0.0;
  sigma_econf = 0.0;
  valarray<double> sum(0.0,14), tsum(0.0,14);
Francois Gygi committed
431 432 433 434
  for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
  {
    for ( int ikp = 0; ikp < wf.nkp(); ikp++ )
    {
435
      const double weight = wf.weight(ikp);
Francois Gygi committed
436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453
      const SlaterDet& sd = *(wf.sd(ispin,ikp));
      const Basis& wfbasis = sd.basis();
      // 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++ )
454
      {
Francois Gygi committed
455 456
        double tmpsum = 0.0;
        for ( int n = 0; n < nloc; n++ )
457
        {
Francois Gygi committed
458 459
          const double psi2 = norm(p[ig+n*mloc]);
          tmpsum += fac * occ[n] * psi2;
460
        }
Francois Gygi committed
461 462
        psi2sum[ig] = tmpsum;
      }
463

Francois Gygi committed
464 465 466 467 468 469
      // 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;
470

Francois Gygi committed
471 472 473
      for ( int ig = 0; ig < ngwloc; ig++ )
      {
        const double psi2s = psi2sum[ig];
474

Francois Gygi committed
475 476
        // tsum[0]: ekin partial sum
        tsum[0] += psi2s * kpg2[ig];
477

Francois Gygi committed
478 479 480 481 482
        if ( compute_stress )
        {
          const double tkpgx = kpg_x[ig];
          const double tkpgy = kpg_y[ig];
          const double tkpgz = kpg_z[ig];
483

Francois Gygi committed
484
          const double fac_ekin = 2.0 * psi2s;
485

Francois Gygi committed
486 487 488 489 490 491
          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;
492

Francois Gygi committed
493 494 495 496
        }
        // tsum[0-6] contains the contributions to
        // ekin, sigma_ekin, from vector ig
      } // ig
497

Francois Gygi committed
498 499 500 501 502
      if ( use_confinement )
      {
        const valarray<double>& fstress = cfp[ikp]->fstress();
        const valarray<double>& dfstress = cfp[ikp]->dfstress();
        for ( int ig = 0; ig < ngwloc; ig++ )
503
        {
Francois Gygi committed
504 505 506
          const double psi2s = psi2sum[ig];
          // tsum[7]: econf partial sum
          tsum[7] += psi2s * fstress[ig];
507

Francois Gygi committed
508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524
          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
525
      }
Francois Gygi committed
526 527

      sum += weight * tsum;
528 529
    } // ikp
  } // ispin
530

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

534
  ekin_  = sum[0];
535 536 537 538 539 540
  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];
541

542
  econf_ = sum[7];
543 544 545 546 547 548
  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];
549

550 551
  sigma_ekin *= omega_inv;
  sigma_econf *= omega_inv;
552

Francois Gygi committed
553
  tmap["ekin"].stop();
554

555 556 557 558
  // Stress from Eps
  sigma_eps = 0.0;
  if ( compute_stress )
  {
559
    sum = 0.0;
560 561 562 563 564 565 566 567
    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)
568
      const double fac = 2.0 * gi[ig] *
569
        real( conj(rhoelg[ig]) * dvion_local_g[ig] );
570

571 572 573
      const double tgx = g_x[ig];
      const double tgy = g_y[ig];
      const double tgz = g_z[ig];
574

575 576 577 578 579 580
      sum[0] += fac * tgx * tgx;
      sum[1] += fac * tgy * tgy;
      sum[2] += fac * tgz * tgz;
      sum[3] += fac * tgx * tgy;
      sum[4] += fac * tgy * tgz;
      sum[5] += fac * tgx * tgz;
581
    }
582
    MPI_Allreduce(&sum[0],&tsum[0],6,MPI_DOUBLE,MPI_SUM,vbasis_->comm());
Francois Gygi committed
583

Francois Gygi committed
584 585 586 587 588 589
    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
590
  }
591

592 593 594
  // Stress from Hartree energy
  if ( compute_stress )
  {
595
    sum = 0.0;
596 597 598
    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);
599

600 601 602 603 604 605 606
    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];

607 608 609 610 611 612
      sum[0] += temp * tgx * tgx;
      sum[1] += temp * tgy * tgy;
      sum[2] += temp * tgz * tgz;
      sum[3] += temp * tgx * tgy;
      sum[4] += temp * tgy * tgz;
      sum[5] += temp * tgx * tgz;
613 614 615 616 617 618 619 620 621 622 623 624
    }

    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];
625 626 627
        Species *s = s_.atoms.species_list[is];
        const double g = vbasis_->g_ptr()[ig];
        const double gi = vbasis_->gi_ptr()[ig];
628
        // next line: keep only real part
629 630 631 632 633 634 635 636 637
        // contribution of pseudocharge of ion
        double temp = fac2 * (rg.real() * sg.real() + rg.imag() * sg.imag())
          * rhops[is][ig] * g2i[ig];
        if ( core_charge_ )
        {
          double rhoc_g, drhoc_g;
          s->drhocoreg(g,rhoc_g,drhoc_g);
          // next line: keep only real part
          // contribution of core correction
638
          temp -= (vxc_g[ig].real()*sg.real() + vxc_g[ig].imag()*sg.imag())
639 640
               * drhoc_g * gi * omega_inv / fpi;
        }
641

642 643 644 645
        const double tgx = g_x[ig];
        const double tgy = g_y[ig];
        const double tgz = g_z[ig];

646 647 648 649 650 651
        sum[0] += temp * tgx * tgx;
        sum[1] += temp * tgy * tgy;
        sum[2] += temp * tgz * tgz;
        sum[3] += temp * tgx * tgy;
        sum[4] += temp * tgy * tgz;
        sum[5] += temp * tgx * tgz;
652 653
      }
    }
654
    MPI_Allreduce(&sum[0],&tsum[0],6,MPI_DOUBLE,MPI_SUM,vbasis_->comm());
Francois Gygi committed
655 656 657 658 659 660 661 662 663 664
    // 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];
665
  } // compute_stress
666

667
  // Non local energy and forces
Francois Gygi committed
668
  tmap["nonlocal"].start();
Francois Gygi committed
669 670
  // modify next loop to span only local ikp
  enl_ = 0.0;
671 672 673 674 675 676
  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
677
  for ( int ikp = 0; ikp < wf.nkp(); ikp++ )
678
  {
Francois Gygi committed
679
    for ( int ispin = 0; ispin < nlp.size(); ispin++ )
680
    {
Francois Gygi committed
681 682 683
      enl_ += wf.weight(ikp) * nlp[ispin][ikp]->energy(compute_hpsi,
              *dwf.sd(ispin,ikp), compute_forces, fion_enl, compute_stress,
              sigma_enl_kp);
684

685 686 687 688
      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];
689

690 691 692
      if ( compute_stress )
        sigma_enl += wf.weight(ikp) * sigma_enl_kp;
    }
693
  }
Francois Gygi committed
694 695 696
  tmap["nonlocal"].stop();

  ecoul_ = ehart_ + esr_ - eself_;
697 698 699 700 701 702 703
  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;
  }
704
  etotal_ = ekin_ + econf_ + eps_ + enl_ + ecoul_ + exc_ + ets_ + eexf_;
Francois Gygi committed
705 706
  if ( s_.vext )
    etotal_ += eext_;
707
  enthalpy_ = etotal_;
708 709 710 711 712 713

  // Electric enthalpy
  eefield_ = 0.0;
  if ( el_enth_ )
  {
    tmap["el_enth_energy"].start();
714
    eefield_ = el_enth_->enthalpy(dwf,compute_hpsi);
715
    tmap["el_enth_energy"].stop();
716
    enthalpy_ += eefield_;
717 718 719 720 721 722 723 724 725 726 727 728 729

    if ( compute_forces )
    {
      for ( int is = 0; is < nsp_; is++ )
        for ( int ia = 0; ia < na_[is]; ia++ )
        {
          D3vector f = zv_[is] * s_.ctrl.e_field;
          fion[is][3*ia]   += f.x;
          fion[is][3*ia+1] += f.y;
          fion[is][3*ia+2] += f.z;
        }
    }
  }
730

731
  epv_ = 0.0;
732 733 734
  if ( compute_stress )
  {
    valarray<double> sigma_ext(s_.ctrl.ext_stress,6);
Francois Gygi committed
735 736
    const double gpa = 29421.5;
    const double pext = (sigma_ext[0]+sigma_ext[1]+sigma_ext[2])/(3.0*gpa);
737 738 739
    epv_ = pext * omega;
    enthalpy_ += epv_;
  }
740

Francois Gygi committed
741 742 743 744 745 746 747
  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
748 749 750 751 752 753 754 755 756 757 758
        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
759
        {
Francois Gygi committed
760
          for ( int n = 0; n < sd.nstloc(); n++ )
Francois Gygi committed
761
          {
Francois Gygi committed
762 763
            const valarray<double>& fstress = cfp[ikp]->fstress();
            for ( int ig = 0; ig < ngwloc; ig++ )
764
            {
Francois Gygi committed
765 766
              cp[ig+mloc*n] += 0.5 * ( kpg2[ig] + fstress[ig] ) *
                               c[ig+mloc*n];
767
            }
768
          }
Francois Gygi committed
769 770 771 772
        }
        else
        {
          for ( int n = 0; n < sd.nstloc(); n++ )
773
          {
Francois Gygi committed
774
            for ( int ig = 0; ig < ngwloc; ig++ )
Francois Gygi committed
775
            {
Francois Gygi committed
776
              cp[ig+mloc*n] += 0.5 * kpg2[ig] * c[ig+mloc*n];
Francois Gygi committed
777 778
            }
          }
Francois Gygi committed
779
        }
Francois Gygi committed
780 781

        // local potential
Francois Gygi committed
782
        sd.rs_mul_add(*ft[ikp], &v_r[ispin][0], sdp);
Francois Gygi committed
783 784 785 786
      } //ikp
    } //ispin

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

Francois Gygi committed
789 790
    tmap["hpsi"].stop();
  } // if compute_hpsi
791

Francois Gygi committed
792 793 794 795 796 797
  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);
798
    valarray<double> fion_nl, fion_nl_tmp;
799

Francois Gygi committed
800 801 802 803 804 805
    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]);
806 807
        if ( core_charge_ )
          vtemp[ig] += conj(vxc_g[ig]) * rhocore_sp_g[is][ig];
Francois Gygi committed
808
      }
809 810 811
      fion_nl.resize(3*na_[is]);
      fion_nl = 0.0;
      fion_nl_tmp.resize(3*na_[is]);
Francois Gygi committed
812 813 814 815 816 817 818 819 820 821 822 823 824
      // 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
825 826 827
          // complex<double> teigr = ei0[kv[3*ig+0]] *
          //                         ei1[kv[3*ig+1]] *
          //                         ei2[kv[3*ig+2]];
Francois Gygi committed
828 829 830 831
          const int iii = ig+ig+ig;
          const int kx = idx[iii];
          const int ky = idx[iii+1];
          const int kz = idx[iii+2];
832

Francois Gygi committed
833 834 835
          const double cos_a = c0[kx];
          const double cos_b = c1[ky];
          const double cos_c = c2[kz];
836

Francois Gygi committed
837 838 839
          const double sin_a = s0[kx];
          const double sin_b = s1[ky];
          const double sin_c = s2[kz];
840

Francois Gygi committed
841 842
          // Next line: exp(-i*gr) =
          // (cos_a - I sin_a)*(cos_b - I sin_b)*(cos_c - I sin_c)
843
          double teigr_re =
Francois Gygi committed
844 845
            cos_a*cos_b*cos_c - sin_a*sin_b*cos_c -
            sin_a*cos_b*sin_c - cos_a*sin_b*sin_c;
846
          double teigr_im =
Francois Gygi committed
847 848
            sin_a*sin_b*sin_c - sin_a*cos_b*cos_c -
            cos_a*sin_b*cos_c - cos_a*cos_b*sin_c;
849

Francois Gygi committed
850
          /* fion is real */
851
          double tmp = teigr_re * vtemp[ig].imag() +
Francois Gygi committed
852 853 854 855 856 857
                       teigr_im * vtemp[ig].real();

          sum0 += tmp * gx0[ig];
          sum1 += tmp * gx1[ig];
          sum2 += tmp * gx2[ig];
        }
858 859 860
        fion_nl[3*ia]   = sum0;
        fion_nl[3*ia+1] = sum1;
        fion_nl[3*ia+2] = sum2;
Francois Gygi committed
861 862

      }
863

864 865
      MPI_Allreduce(&fion_nl[0],&fion_nl_tmp[0],3*na_[is],
                    MPI_DOUBLE,MPI_SUM,vbasis_->comm());
Francois Gygi committed
866 867 868 869

      double fac = -2.0 * omega;
      for ( int i = 0; i < 3*na_[is]; i++ )
      {
870
        fion[is][i] += fion_esr[is][i] + fac * fion_nl_tmp[i];
Francois Gygi committed
871
      }
Francois Gygi committed
872 873 874 875 876 877 878 879 880

      // 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
881 882
    }
  }
883

884
  if ( compute_stress )
885
  {
886 887
    sigma = sigma_ekin + sigma_econf + sigma_eps + sigma_enl +
            sigma_ehart + sigma_exc + sigma_esr;
888
  }
889

890
  if ( debug_stress && s_.ctxt_.onpe0() )
891 892 893 894
  {
    cout.setf(ios::fixed,ios::floatfield);
    cout.setf(ios::right,ios::adjustfield);
    cout << setprecision(8);
895
    cout << " <stress_tensor unit=\"atomic_units\">\n"
896 897 898 899 900 901 902 903 904 905 906 907
         << "   <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"
908
         << endl
909 910 911 912 913 914 915 916 917 918 919 920
         << "   <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"
921
         << endl
922 923 924 925 926 927 928 929 930 931 932 933
         << "   <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"
934
         << endl
935 936 937 938 939 940 941 942 943 944 945 946
         << "   <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"
947
         << endl
948 949 950 951 952 953 954 955 956 957 958 959
         << "   <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"
960
         << endl
961 962 963 964 965 966 967 968 969 970 971 972
         << "   <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"
973
         << endl
974 975 976 977 978 979 980 981 982 983 984 985
         << "   <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"
986
         << endl
987 988 989 990 991 992
         << "   <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"
993 994
         << " </stress_tensor>" << endl;
  }
Francois Gygi committed
995 996 997 998 999 1000 1001
  return etotal_;
}

////////////////////////////////////////////////////////////////////////////////
void EnergyFunctional::atoms_moved(void)
{
  const AtomSet& atoms = s_.atoms;
1002
  const Wavefunction& wf = s_.wf;
1003
  const UnitCell& cell = s_.wf.cell();
Francois Gygi committed
1004 1005
  int ngloc = vbasis_->localsize();

1006
  // fill tau0 with values in atom_list
1007

Francois Gygi committed
1008 1009
  atoms.get_positions(tau0);
  sf.update(tau0,*vbasis_);
1010

Francois Gygi committed
1011 1012
  // compute Fourier coefficients of the local potential
  memset( (void*)&vion_local_g[0], 0, 2*ngloc*sizeof(double) );
1013
  memset( (void*)&dvion_local_g[0], 0, 2*ngloc*sizeof(double) );
Francois Gygi committed
1014
  memset( (void*)&rhopst[0], 0, 2*ngloc*sizeof(double) );
1015 1016 1017 1018

  if ( core_charge_ )
  {
    // recalculate Fourier coefficients of the core charge
1019 1020 1021 1022
    assert(cd_.rhocore_g.size()==ngloc);
    memset( (void*)&cd_.rhocore_g[0], 0, 2*ngloc*sizeof(double) );
    // divide core charge by two if two spins
    const double spin_fac = cell.volume() / wf.nspin();
1023 1024 1025
    for ( int is = 0; is < atoms.nsp(); is++ )
    {
      complex<double> *s = &sf.sfac[is][0];
1026
      double *r = &rhocore_sp_g[is][0];
1027
      for ( int ig = 0; ig < ngloc; ig++ )
1028
        cd_.rhocore_g[ig] += spin_fac * s[ig] * r[ig];
1029 1030
    }
    // recalculate real-space core charge
1031 1032
    vft->backward(&cd_.rhocore_g[0],&tmp_r[0]);
    const double omega_inv = 1.0 / cell.volume();
1033
    for ( int i = 0; i < tmp_r.size(); i++ )
1034
      cd_.rhocore_r[i] = tmp_r[i].real() * omega_inv;
1035 1036
  }

Francois Gygi committed
1037 1038 1039 1040 1041
  for ( int is = 0; is < atoms.nsp(); is++ )
  {
    complex<double> *s = &sf.sfac[is][0];
    for ( int ig = 0; ig < ngloc; ig++ )
    {
1042 1043 1044 1045
      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
1046 1047
    }
  }
1048

Francois Gygi committed
1049
  // compute esr: pseudocharge repulsion energy
1050
  const double omega_inv = 1.0 / cell.volume();
Francois Gygi committed
1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062
  // 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
1063

Francois Gygi committed
1064
  esr_  = 0.0;
1065
  sigma_esr = 0.0;
Francois Gygi committed
1066 1067 1068 1069
  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
1070
  for ( int is1 = 0; is1 < nsp_; is1++ )
Francois Gygi committed
1071
  {
Francois Gygi committed
1072
    for ( int is2 = is1; is2 < nsp_; is2++ )
Francois Gygi committed
1073
    {
Francois Gygi committed
1074 1075 1076
      double rcps12 = sqrt ( rcps_[is1]*rcps_[is1]+rcps_[is2]*rcps_[is2] );
      // convergence criterion for lattice sums:
      // fac * rcps12 < ncell * d
1077
      const double fac = 8.0;
Francois Gygi committed
1078 1079 1080 1081 1082 1083 1084 1085 1086 1087
      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
1088
      {
Francois Gygi committed
1089
        int ia2min = 0;
1090
        if ( is1 == is2 ) ia2min = ia1;
Francois Gygi committed
1091
        for ( int ia2 = ia2min; ia2 < na_[is2]; ia2++ )
Francois Gygi committed
1092
        {
1093
          const bool same_atom = ( is1==is2 && ia1==ia2 );
Francois Gygi committed
1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104
          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++ )
              {
1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137 1138
                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
1139
              }
Francois Gygi committed
1140 1141 1142 1143
        }
      }
    }
  }
1144
  sigma_esr *= - omega_inv;
Francois Gygi committed
1145 1146 1147 1148

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

Francois Gygi committed
1149 1150 1151 1152 1153
}

////////////////////////////////////////////////////////////////////////////////
void EnergyFunctional::cell_moved(void)
{
Francois Gygi committed
1154 1155
  const Wavefunction& wf = s_.wf;
  const UnitCell& cell = wf.cell();
1156
  // resize vbasis_
Francois Gygi committed
1157
  vbasis_->resize(cell,s_.wf.refcell(),4.0*s_.wf.ecut());
1158

1159 1160 1161
  const int ngloc = vbasis_->localsize();
  const double omega = cell.volume();
  assert(omega != 0.0);
1162 1163
  const double omega_inv = 1.0 / omega;

1164 1165
  const AtomSet& atoms = s_.atoms;
  for ( int is = 0; is < nsp_; is++ )
1166
  {
1167 1168
    Species *s = atoms.species_list[is];
    const double * const g = vbasis_->g_ptr();
1169
    double v,dv,rhoc;
1170
    for ( int ig = 0; ig < ngloc; ig++ )
1171
    {
1172 1173 1174 1175
      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;
1176 1177 1178 1179 1180
      if ( core_charge_ )
      {
        s->rhocoreg(g[ig],rhoc);
        rhocore_sp_g[is][ig] = rhoc * omega_inv;
      }
1181
    }
1182
  }
1183

Francois Gygi committed
1184
  // Update confinement potentials and non-local potentials
1185 1186
  for ( int ikp = 0; ikp < wf.nkp(); ikp++ )
  {
Francois Gygi committed
1187
    cfp[ikp]->update();
Francois Gygi committed
1188 1189
    for ( int ispin = 0; ispin < nlp.size(); ispin++ )
      nlp[ispin][ikp]->update_twnl();
1190
  }
1191 1192 1193

  // Update exchange-correlation operator
  xco->cell_moved();
Francois Gygi committed
1194
}
1195 1196 1197 1198 1199 1200 1201

////////////////////////////////////////////////////////////////////////////////
void EnergyFunctional::print(ostream& os) const
{
  os.setf(ios::fixed,ios::floatfield);
  os.setf(ios::right,ios::adjustfield);
  os << setprecision(8);
1202 1203 1204 1205 1206 1207 1208 1209 1210 1211 1212 1213 1214 1215
  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"
     << "  <etotal>  " << setw(15) << etotal() << " </etotal>\n"
     << "  <epv>     " << setw(15) << epv() << " </epv>\n"
     << "  <eefield> " << setw(15) << eefield() << " </eefield>\n"
     << "  <enthalpy>" << setw(15) << enthalpy() << " </enthalpy>" << endl;
Francois Gygi committed
1216 1217
  if ( s_.vext )
    os << "  <eext>    " << setw(15) << eext() << " </eext>" << endl;
1218
}
1219

1220 1221
////////////////////////////////////////////////////////////////////////////////
ostream& operator<< ( ostream& os, const EnergyFunctional& e )
1222 1223
{
  e.print(os);
1224 1225
  return os;
}