EnergyFunctional.C 36.9 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"
Francois Gygi committed
33 34 35 36 37 38

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

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

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

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

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

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

Francois Gygi committed
67 68 69 70 71 72 73 74 75
  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() )
  {
76 77 78 79
    cout << "  EnergyFunctional: np0v,np1v,np2v: " << np0v << " "
         << np1v << " " << np2v << endl;
    cout << "  EnergyFunctional: vft->np012(): "
         << vft->np012() << endl;
Francois Gygi committed
80
  }
81

82
  const int ngloc = vbasis_->localsize();
83

Francois Gygi committed
84
  nsp_ = atoms.nsp();
85

Francois Gygi committed
86 87 88
  vps.resize(nsp_);
  dvps.resize(nsp_);
  rhops.resize(nsp_);
89

Francois Gygi committed
90 91 92
  zv_.resize(nsp_);
  rcps_.resize(nsp_);
  na_.resize(nsp_);
93

Francois Gygi committed
94 95 96 97 98 99
  for ( int is = 0; is < nsp_; is++ )
  {
    vps[is].resize(ngloc);
    dvps[is].resize(ngloc);
    rhops[is].resize(ngloc);
  }
100

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

  nlp.resize(wf.nspin());
  for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
Francois Gygi committed
105
  {
Francois Gygi committed
106 107 108 109 110
    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
111
  }
112

113 114 115 116 117 118 119
  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);
120

121 122
  tau0.resize(nsp_);
  fion_esr.resize(nsp_);
Francois Gygi committed
123
  fext.resize(nsp_);
124

125
  eself_ = 0.0;
126

127 128
  // core_charge_: true if any species has a core charge density
  core_charge_ = false;
129 130 131
  for ( int is = 0; is < nsp_; is++ )
  {
    Species *s = atoms.species_list[is];
132

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

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

141
    zv_[is] = s->ztot();
142
    rcps_[is] = s->rcps();
143 144 145 146 147 148 149 150 151 152 153 154 155 156

    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);
    rhocore_g.resize(ngloc);
    rhocore_r.resize(vft->np012loc());
    // set rhocore_r ptr in ChargeDensity
    cd_.rhocore_r = &rhocore_r[0];
157
  }
158

Francois Gygi committed
159 160 161 162 163 164
  // 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);
  }
165

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

178 179
  // Electric enthalpy
  el_enth_ = 0;
180
  if ( s_.ctrl.polarization != "OFF" )
181 182
    el_enth_ = new ElectricEnthalpy(s_);

183
  sf.init(tau0,*vbasis_);
184

185
  cell_moved();
186

187
  atoms_moved();
Francois Gygi committed
188

Francois Gygi committed
189 190 191 192 193
}

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

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

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

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

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

253 254
  // update XC operator
  // compute xc energy, update self-energy operator and potential
Francois Gygi committed
255 256
  tmap["exc"].start();
  for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
257 258
    memset((void*)&v_r[ispin][0], 0, vft->np012loc()*sizeof(double));

259
  xco->update(v_r, compute_stress);
Francois Gygi committed
260
  exc_ = xco->exc();
261
  dxc_ = xco->dxc();
262 263
  if ( compute_stress )
    xco->compute_stress(sigma_exc);
264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283

  if ( core_charge_ )
  {
    // compute Fourier coefficients of Vxc
    for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
    {
      copy(v_r[ispin].begin(),v_r[ispin].end(),tmp_r.begin());
      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] );
      }
    }
  }

Francois Gygi committed
284
  tmap["exc"].stop();
285 286

  // compute local potential energy:
Francois Gygi committed
287 288
  // integral of el. charge times ionic local pot.
  int len=2*ngloc,inc1=1;
289
  sum[0] = 2.0 * ddot(&len,(double*)&rhoelg[0],&inc1,
Francois Gygi committed
290 291
         (double*)&vion_local_g[0],&inc1);
  // remove double counting for G=0
292
  if ( vbasis_->mype() == 0 )
Francois Gygi committed
293
  {
294
    sum[0] -= real(conj(rhoelg[0])*vion_local_g[0]);
Francois Gygi committed
295
  }
296
  sum[0] *= omega; // sum[0] contains eps
297

298
  // Hartree energy and electron-electron energy (without pseudocharges)
Francois Gygi committed
299
  double ehsum = 0.0;
300 301 302
  double ehesum = 0.0;
  double ehepsum = 0.0;
  double ehpsum = 0.0;
Francois Gygi committed
303 304
  for ( int ig = 0; ig < ngloc; ig++ )
  {
305 306 307 308 309 310 311
    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];
    rhogt[ig] = r+rp;
Francois Gygi committed
312 313 314
  }
  // factor 1/2 from definition of Ehart cancels with half sum over G
  // Note: rhogt[ig] includes a factor 1/Omega
315
  // Factor omega in next line yields prefactor 4 pi / omega in Ehart
316
  sum[1] = omega * fpi * ehsum;
317 318 319
  sum[2] = omega * fpi * ehesum;
  sum[3] = omega * fpi * ehepsum;
  sum[4] = omega * fpi * ehpsum;
320

321 322 323 324 325 326 327
  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];
328

Francois Gygi committed
329 330 331 332 333 334
  // 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];
  }
335

Francois Gygi committed
336 337
  // FT to tmpr_r
  vft->backward(&vlocal_g[0],&tmp_r[0]);
338

Francois Gygi committed
339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357
  // 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;
    }
  }
358 359 360

  if ( el_enth_ )
    el_enth_->update();
Francois Gygi committed
361 362
}

Francois Gygi committed
363 364 365
////////////////////////////////////////////////////////////////////////////////
double EnergyFunctional::energy(bool compute_hpsi, Wavefunction& dwf,
              bool compute_forces, vector<vector<double> >& fion,
366
              bool compute_stress, valarray<double>& sigma)
Francois Gygi committed
367
{
368
  const bool debug_stress = compute_stress &&
369
    s_.ctrl.debug.find("STRESS") != string::npos;
Francois Gygi committed
370 371 372 373 374 375 376
  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
377
  const double *const g2i = vbasis_->g2i_ptr();
378

379
  const bool use_confinement = s_.ctrl.ecuts > 0.0;
380

381
  for ( int is = 0; is < fion.size(); is++ )
Francois Gygi committed
382 383
    for ( int i = 0; i < fion[is].size(); i++ )
      fion[is][i] = 0.0;
384

Francois Gygi committed
385 386 387 388 389 390
  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();
  }
391

Francois Gygi committed
392 393
  // kinetic energy
  tmap["ekin"].start();
394

395 396 397 398 399
  // 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
400
  ekin_ = 0.0;
401 402 403 404
  econf_ = 0.0;
  sigma_ekin = 0.0;
  sigma_econf = 0.0;
  valarray<double> sum(0.0,14), tsum(0.0,14);
Francois Gygi committed
405 406 407 408
  for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
  {
    for ( int ikp = 0; ikp < wf.nkp(); ikp++ )
    {
409
      const double weight = wf.weight(ikp);
Francois Gygi committed
410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428
      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++ )
429
      {
Francois Gygi committed
430 431
        double tmpsum = 0.0;
        for ( int n = 0; n < nloc; n++ )
432
        {
Francois Gygi committed
433 434
          const double psi2 = norm(p[ig+n*mloc]);
          tmpsum += fac * occ[n] * psi2;
435
        }
Francois Gygi committed
436 437
        psi2sum[ig] = tmpsum;
      }
438

Francois Gygi committed
439 440 441 442 443 444
      // 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;
445

Francois Gygi committed
446 447 448
      for ( int ig = 0; ig < ngwloc; ig++ )
      {
        const double psi2s = psi2sum[ig];
449

Francois Gygi committed
450 451
        // tsum[0]: ekin partial sum
        tsum[0] += psi2s * kpg2[ig];
452

Francois Gygi committed
453 454 455 456 457
        if ( compute_stress )
        {
          const double tkpgx = kpg_x[ig];
          const double tkpgy = kpg_y[ig];
          const double tkpgz = kpg_z[ig];
458

Francois Gygi committed
459
          const double fac_ekin = 2.0 * psi2s;
460

Francois Gygi committed
461 462 463 464 465 466
          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;
467

Francois Gygi committed
468 469 470 471
        }
        // tsum[0-6] contains the contributions to
        // ekin, sigma_ekin, from vector ig
      } // ig
472

Francois Gygi committed
473 474 475 476 477
      if ( use_confinement )
      {
        const valarray<double>& fstress = cfp[ikp]->fstress();
        const valarray<double>& dfstress = cfp[ikp]->dfstress();
        for ( int ig = 0; ig < ngwloc; ig++ )
478
        {
Francois Gygi committed
479 480 481
          const double psi2s = psi2sum[ig];
          // tsum[7]: econf partial sum
          tsum[7] += psi2s * fstress[ig];
482

Francois Gygi committed
483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499
          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
500
      }
Francois Gygi committed
501 502

      sum += weight * tsum;
503 504
    } // ikp
  } // ispin
505

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

509
  ekin_  = sum[0];
510 511 512 513 514 515
  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];
516

517
  econf_ = sum[7];
518 519 520 521 522 523
  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];
524

525 526
  sigma_ekin *= omega_inv;
  sigma_econf *= omega_inv;
527

Francois Gygi committed
528
  tmap["ekin"].stop();
529

530 531 532 533
  // Stress from Eps
  sigma_eps = 0.0;
  if ( compute_stress )
  {
534
    sum = 0.0;
535 536 537 538 539 540 541 542
    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)
543
      const double fac = 2.0 * gi[ig] *
544
        real( conj(rhoelg[ig]) * dvion_local_g[ig] );
545

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

550 551 552 553 554 555
      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;
556
    }
557
    MPI_Allreduce(&sum[0],&tsum[0],6,MPI_DOUBLE,MPI_SUM,vbasis_->comm());
Francois Gygi committed
558

Francois Gygi committed
559 560 561 562 563 564
    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
565
  }
566

567 568 569
  // Stress from Hartree energy
  if ( compute_stress )
  {
570
    sum = 0.0;
571 572 573
    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);
574

575 576 577 578 579 580 581
    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];

582 583 584 585 586 587
      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;
588 589 590 591 592 593 594 595 596 597 598 599
    }

    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];
600 601 602
        Species *s = s_.atoms.species_list[is];
        const double g = vbasis_->g_ptr()[ig];
        const double gi = vbasis_->gi_ptr()[ig];
603
        // next line: keep only real part
604 605 606 607 608 609 610 611 612 613 614 615
        // 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
          temp -= (vxc_g[ig].real() * sg.real() + vxc_g[ig].imag() * sg.imag())
               * drhoc_g * gi * omega_inv / fpi;
        }
616

617 618 619 620
        const double tgx = g_x[ig];
        const double tgy = g_y[ig];
        const double tgz = g_z[ig];

621 622 623 624 625 626
        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;
627 628
      }
    }
629
    MPI_Allreduce(&sum[0],&tsum[0],6,MPI_DOUBLE,MPI_SUM,vbasis_->comm());
Francois Gygi committed
630 631 632 633 634 635 636 637 638 639
    // 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];
640
  } // compute_stress
641

642
  // Non local energy and forces
Francois Gygi committed
643
  tmap["nonlocal"].start();
Francois Gygi committed
644 645
  // modify next loop to span only local ikp
  enl_ = 0.0;
646 647 648 649 650 651
  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
652
  for ( int ikp = 0; ikp < wf.nkp(); ikp++ )
653
  {
Francois Gygi committed
654
    for ( int ispin = 0; ispin < nlp.size(); ispin++ )
655
    {
Francois Gygi committed
656 657 658
      enl_ += wf.weight(ikp) * nlp[ispin][ikp]->energy(compute_hpsi,
              *dwf.sd(ispin,ikp), compute_forces, fion_enl, compute_stress,
              sigma_enl_kp);
659

660 661 662 663
      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];
664

665 666 667
      if ( compute_stress )
        sigma_enl += wf.weight(ikp) * sigma_enl_kp;
    }
668
  }
Francois Gygi committed
669 670 671
  tmap["nonlocal"].stop();

  ecoul_ = ehart_ + esr_ - eself_;
672 673 674 675 676 677 678
  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;
  }
679
  etotal_ = ekin_ + econf_ + eps_ + enl_ + ecoul_ + exc_ + ets_ + eexf_;
680
  enthalpy_ = etotal_;
681 682 683 684 685 686

  // Electric enthalpy
  eefield_ = 0.0;
  if ( el_enth_ )
  {
    tmap["el_enth_energy"].start();
687
    eefield_ = el_enth_->enthalpy(dwf,compute_hpsi);
688
    tmap["el_enth_energy"].stop();
689
    enthalpy_ += eefield_;
690 691 692 693 694 695 696 697 698 699 700 701 702

    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;
        }
    }
  }
703

704
  epv_ = 0.0;
705 706 707 708 709 710 711
  if ( compute_stress )
  {
    valarray<double> sigma_ext(s_.ctrl.ext_stress,6);
    const double pext = (sigma_ext[0]+sigma_ext[1]+sigma_ext[2])/3.0;
    epv_ = pext * omega;
    enthalpy_ += epv_;
  }
712

Francois Gygi committed
713 714 715 716 717 718 719
  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
720 721 722 723 724 725 726 727 728 729 730
        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
731
        {
Francois Gygi committed
732
          for ( int n = 0; n < sd.nstloc(); n++ )
Francois Gygi committed
733
          {
Francois Gygi committed
734 735
            const valarray<double>& fstress = cfp[ikp]->fstress();
            for ( int ig = 0; ig < ngwloc; ig++ )
736
            {
Francois Gygi committed
737 738
              cp[ig+mloc*n] += 0.5 * ( kpg2[ig] + fstress[ig] ) *
                               c[ig+mloc*n];
739
            }
740
          }
Francois Gygi committed
741 742 743 744
        }
        else
        {
          for ( int n = 0; n < sd.nstloc(); n++ )
745
          {
Francois Gygi committed
746
            for ( int ig = 0; ig < ngwloc; ig++ )
Francois Gygi committed
747
            {
Francois Gygi committed
748
              cp[ig+mloc*n] += 0.5 * kpg2[ig] * c[ig+mloc*n];
Francois Gygi committed
749 750
            }
          }
Francois Gygi committed
751
        }
Francois Gygi committed
752 753

        // local potential
Francois Gygi committed
754
        sd.rs_mul_add(*ft[ikp], &v_r[ispin][0], sdp);
Francois Gygi committed
755 756 757 758
      } //ikp
    } //ispin

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

Francois Gygi committed
761 762
    tmap["hpsi"].stop();
  } // if compute_hpsi
763

Francois Gygi committed
764 765 766 767 768 769
  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);
770
    valarray<double> fion_nl, fion_nl_tmp;
771

Francois Gygi committed
772 773 774 775 776 777
    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]);
778 779
        if ( core_charge_ )
          vtemp[ig] += conj(vxc_g[ig]) * rhocore_sp_g[is][ig];
Francois Gygi committed
780
      }
781 782 783
      fion_nl.resize(3*na_[is]);
      fion_nl = 0.0;
      fion_nl_tmp.resize(3*na_[is]);
Francois Gygi committed
784 785 786 787 788 789 790 791 792 793 794 795 796
      // 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
797 798 799
          // complex<double> teigr = ei0[kv[3*ig+0]] *
          //                         ei1[kv[3*ig+1]] *
          //                         ei2[kv[3*ig+2]];
Francois Gygi committed
800 801 802 803
          const int iii = ig+ig+ig;
          const int kx = idx[iii];
          const int ky = idx[iii+1];
          const int kz = idx[iii+2];
804

Francois Gygi committed
805 806 807
          const double cos_a = c0[kx];
          const double cos_b = c1[ky];
          const double cos_c = c2[kz];
808

Francois Gygi committed
809 810 811
          const double sin_a = s0[kx];
          const double sin_b = s1[ky];
          const double sin_c = s2[kz];
812

Francois Gygi committed
813 814
          // Next line: exp(-i*gr) =
          // (cos_a - I sin_a)*(cos_b - I sin_b)*(cos_c - I sin_c)
815
          double teigr_re =
Francois Gygi committed
816 817
            cos_a*cos_b*cos_c - sin_a*sin_b*cos_c -
            sin_a*cos_b*sin_c - cos_a*sin_b*sin_c;
818
          double teigr_im =
Francois Gygi committed
819 820
            sin_a*sin_b*sin_c - sin_a*cos_b*cos_c -
            cos_a*sin_b*cos_c - cos_a*cos_b*sin_c;
821

Francois Gygi committed
822
          /* fion is real */
823
          double tmp = teigr_re * vtemp[ig].imag() +
Francois Gygi committed
824 825 826 827 828 829
                       teigr_im * vtemp[ig].real();

          sum0 += tmp * gx0[ig];
          sum1 += tmp * gx1[ig];
          sum2 += tmp * gx2[ig];
        }
830 831 832
        fion_nl[3*ia]   = sum0;
        fion_nl[3*ia+1] = sum1;
        fion_nl[3*ia+2] = sum2;
Francois Gygi committed
833 834

      }
835

836 837
      MPI_Allreduce(&fion_nl[0],&fion_nl_tmp[0],3*na_[is],
                    MPI_DOUBLE,MPI_SUM,vbasis_->comm());
Francois Gygi committed
838 839 840 841

      double fac = -2.0 * omega;
      for ( int i = 0; i < 3*na_[is]; i++ )
      {
842
        fion[is][i] += fion_esr[is][i] + fac * fion_nl_tmp[i];
Francois Gygi committed
843
      }
Francois Gygi committed
844 845 846 847 848 849 850 851 852

      // 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
853 854
    }
  }
855

856
  if ( compute_stress )
857
  {
858 859
    sigma = sigma_ekin + sigma_econf + sigma_eps + sigma_enl +
            sigma_ehart + sigma_exc + sigma_esr;
860
  }
861

862
  if ( debug_stress && s_.ctxt_.onpe0() )
863
  {
Francois Gygi committed
864
    //const double gpa = 29421.5;
865 866 867
    cout.setf(ios::fixed,ios::floatfield);
    cout.setf(ios::right,ios::adjustfield);
    cout << setprecision(8);
868
    cout << " <stress_tensor unit=\"atomic_units\">\n"
869 870 871 872 873 874 875 876 877 878 879 880
         << "   <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"
881
         << endl
882 883 884 885 886 887 888 889 890 891 892 893
         << "   <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"
894
         << endl
895 896 897 898 899 900 901 902 903 904 905 906
         << "   <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"
907
         << endl
908 909 910 911 912 913 914 915 916 917 918 919
         << "   <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"
920
         << endl
921 922 923 924 925 926 927 928 929 930 931 932
         << "   <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"
933
         << endl
934 935 936 937 938 939 940 941 942 943 944 945
         << "   <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"
946
         << endl
947 948 949 950 951 952 953 954 955 956 957 958
         << "   <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"
959
         << endl
960 961 962 963 964 965
         << "   <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"
966 967
         << " </stress_tensor>" << endl;
  }
Francois Gygi committed
968 969 970 971 972 973 974
  return etotal_;
}

////////////////////////////////////////////////////////////////////////////////
void EnergyFunctional::atoms_moved(void)
{
  const AtomSet& atoms = s_.atoms;
975
  const Wavefunction& wf = s_.wf;
Francois Gygi committed
976 977
  int ngloc = vbasis_->localsize();

978
  // fill tau0 with values in atom_list
979

Francois Gygi committed
980 981
  atoms.get_positions(tau0);
  sf.update(tau0,*vbasis_);
982

Francois Gygi committed
983 984
  // compute Fourier coefficients of the local potential
  memset( (void*)&vion_local_g[0], 0, 2*ngloc*sizeof(double) );
985
  memset( (void*)&dvion_local_g[0], 0, 2*ngloc*sizeof(double) );
Francois Gygi committed
986
  memset( (void*)&rhopst[0], 0, 2*ngloc*sizeof(double) );
987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009

  if ( core_charge_ )
  {
    // recalculate Fourier coefficients of the core charge
    assert(rhocore_g.size()==ngloc);
    memset( (void*)&rhocore_g[0], 0, 2*ngloc*sizeof(double) );
    const double spin_fac = wf.cell().volume() / wf.nspin();
    for ( int is = 0; is < atoms.nsp(); is++ )
    {
      complex<double> *s = &sf.sfac[is][0];
      for ( int ig = 0; ig < ngloc; ig++ )
      {
        const complex<double> sg = s[ig];
        // divide core charge by two if two spins
        rhocore_g[ig] += spin_fac * sg * rhocore_sp_g[is][ig];
      }
    }
    // recalculate real-space core charge
    vft->backward(&rhocore_g[0],&tmp_r[0]);
    for ( int i = 0; i < tmp_r.size(); i++ )
      rhocore_r[i] = tmp_r[i].real();
  }

Francois Gygi committed
1010 1011 1012 1013 1014
  for ( int is = 0; is < atoms.nsp(); is++ )
  {
    complex<double> *s = &sf.sfac[is][0];
    for ( int ig = 0; ig < ngloc; ig++ )
    {
1015 1016 1017 1018
      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
1019 1020
    }
  }
1021

Francois Gygi committed
1022 1023
  // compute esr: pseudocharge repulsion energy
  const UnitCell& cell = s_.wf.cell();
1024
  const double omega_inv = 1.0 / cell.volume();
Francois Gygi committed
1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036
  // 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
1037

Francois Gygi committed
1038
  esr_  = 0.0;
1039
  sigma_esr = 0.0;
Francois Gygi committed
1040 1041 1042 1043
  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
1044
  for ( int is1 = 0; is1 < nsp_; is1++ )
Francois Gygi committed
1045
  {
Francois Gygi committed
1046
    for ( int is2 = is1; is2 < nsp_; is2++ )
Francois Gygi committed
1047
    {
Francois Gygi committed
1048 1049 1050
      double rcps12 = sqrt ( rcps_[is1]*rcps_[is1]+rcps_[is2]*rcps_[is2] );
      // convergence criterion for lattice sums:
      // fac * rcps12 < ncell * d
1051
      const double fac = 8.0;
Francois Gygi committed
1052 1053 1054 1055 1056 1057 1058 1059 1060 1061
      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
1062
      {
Francois Gygi committed
1063
        int ia2min = 0;
1064
        if ( is1 == is2 ) ia2min = ia1;
Francois Gygi committed
1065
        for ( int ia2 = ia2min; ia2 < na_[is2]; ia2++ )
Francois Gygi committed
1066
        {
1067
          const bool same_atom = ( is1==is2 && ia1==ia2 );
Francois Gygi committed
1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078
          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++ )
              {
1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112
                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
1113
              }
Francois Gygi committed
1114 1115 1116 1117
        }
      }
    }
  }
1118
  sigma_esr *= - omega_inv;
Francois Gygi committed
1119 1120 1121 1122

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

Francois Gygi committed
1123 1124 1125 1126 1127
}

////////////////////////////////////////////////////////////////////////////////
void EnergyFunctional::cell_moved(void)
{
Francois Gygi committed
1128 1129
  const Wavefunction& wf = s_.wf;
  const UnitCell& cell = wf.cell();
1130
  // resize vbasis_
Francois Gygi committed
1131
  vbasis_->resize(cell,s_.wf.refcell(),4.0*s_.wf.ecut());
1132

1133 1134 1135
  const int ngloc = vbasis_->localsize();
  const double omega = cell.volume();
  assert(omega != 0.0);
1136 1137
  const double omega_inv = 1.0 / omega;

1138 1139
  const AtomSet& atoms = s_.atoms;
  for ( int is = 0; is < nsp_; is++ )
1140
  {
1141 1142
    Species *s = atoms.species_list[is];
    const double * const g = vbasis_->g_ptr();
1143
    double v,dv,rhoc;
1144
    for ( int ig = 0; ig < ngloc; ig++ )
1145
    {
1146 1147 1148 1149
      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;
1150 1151 1152 1153 1154
      if ( core_charge_ )
      {
        s->rhocoreg(g[ig],rhoc);
        rhocore_sp_g[is][ig] = rhoc * omega_inv;
      }
1155
    }
1156
  }
1157

Francois Gygi committed
1158
  // Update confinement potentials and non-local potentials
1159 1160
  for ( int ikp = 0; ikp < wf.nkp(); ikp++ )
  {
Francois Gygi committed
1161
    cfp[ikp]->update();
Francois Gygi committed
1162 1163
    for ( int ispin = 0; ispin < nlp.size(); ispin++ )
      nlp[ispin][ikp]->update_twnl();
1164
  }
1165 1166 1167

  // Update exchange-correlation operator
  xco->cell_moved();
Francois Gygi committed
1168
}
1169 1170 1171 1172 1173 1174 1175

////////////////////////////////////////////////////////////////////////////////
void EnergyFunctional::print(ostream& os) const
{
  os.setf(ios::fixed,ios::floatfield);
  os.setf(ios::right,ios::adjustfield);
  os << setprecision(8);
1176 1177 1178 1179 1180 1181 1182 1183 1184 1185 1186 1187 1188 1189
  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;
1190
}
1191

1192 1193
////////////////////////////////////////////////////////////////////////////////
ostream& operator<< ( ostream& os, const EnergyFunctional& e )
1194 1195
{
  e.print(os);
1196 1197
  return os;
}