EnergyFunctional.C 36.5 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 213 214 215 216
      cout << "<timing name=\""
           << setw(15) << (*i).first << "\""
           << " min=\"" << setprecision(3) << setw(9) << tmin << "\""
           << " max=\"" << setprecision(3) << setw(9) << tmax << "\"/>"
           << 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[2], tsum[2];
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 262
  if ( compute_stress )
    xco->compute_stress(sigma_exc);
263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282

  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
283
  tmap["exc"].stop();
284 285

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

Francois Gygi committed
297 298 299 300 301 302 303 304 305 306 307 308
  // Hartree energy
  ehart_ = 0.0;
  double ehsum = 0.0;
  for ( int ig = 0; ig < ngloc; ig++ )
  {
    const complex<double> tmp = rhoelg[ig] + rhopst[ig];
    ehsum += norm(tmp) * g2i[ig];
    rhogt[ig] = tmp;
  }
  // factor 1/2 from definition of Ehart cancels with half sum over G
  // Note: rhogt[ig] includes a factor 1/Omega
  // Factor omega in next line yields prefactor 4 pi / omega in
309
  sum[1] = omega * fpi * ehsum;
Francois Gygi committed
310
  // tsum[1] contains ehart
311

312
  MPI_Allreduce(sum,tsum,2,MPI_DOUBLE,MPI_SUM,vbasis_->comm());
Francois Gygi committed
313 314
  eps_   = tsum[0];
  ehart_ = tsum[1];
315

Francois Gygi committed
316 317 318 319 320 321
  // 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];
  }
322

Francois Gygi committed
323 324
  // FT to tmpr_r
  vft->backward(&vlocal_g[0],&tmp_r[0]);
325

Francois Gygi committed
326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344
  // 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;
    }
  }
345 346 347

  if ( el_enth_ )
    el_enth_->update();
Francois Gygi committed
348 349
}

Francois Gygi committed
350 351 352
////////////////////////////////////////////////////////////////////////////////
double EnergyFunctional::energy(bool compute_hpsi, Wavefunction& dwf,
              bool compute_forces, vector<vector<double> >& fion,
353
              bool compute_stress, valarray<double>& sigma)
Francois Gygi committed
354
{
355
  const bool debug_stress = compute_stress &&
356
    s_.ctrl.debug.find("STRESS") != string::npos;
Francois Gygi committed
357 358 359 360 361 362 363
  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
364
  const double *const g2i = vbasis_->g2i_ptr();
365

366
  const bool use_confinement = s_.ctrl.ecuts > 0.0;
367

368
  for ( int is = 0; is < fion.size(); is++ )
Francois Gygi committed
369 370
    for ( int i = 0; i < fion[is].size(); i++ )
      fion[is][i] = 0.0;
371

Francois Gygi committed
372 373 374 375 376 377
  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();
  }
378

Francois Gygi committed
379 380
  // kinetic energy
  tmap["ekin"].start();
381

382 383 384 385 386
  // 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
387
  ekin_ = 0.0;
388 389 390 391
  econf_ = 0.0;
  sigma_ekin = 0.0;
  sigma_econf = 0.0;
  valarray<double> sum(0.0,14), tsum(0.0,14);
Francois Gygi committed
392 393 394 395
  for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
  {
    for ( int ikp = 0; ikp < wf.nkp(); ikp++ )
    {
396
      const double weight = wf.weight(ikp);
Francois Gygi committed
397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415
      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++ )
416
      {
Francois Gygi committed
417 418
        double tmpsum = 0.0;
        for ( int n = 0; n < nloc; n++ )
419
        {
Francois Gygi committed
420 421
          const double psi2 = norm(p[ig+n*mloc]);
          tmpsum += fac * occ[n] * psi2;
422
        }
Francois Gygi committed
423 424
        psi2sum[ig] = tmpsum;
      }
425

Francois Gygi committed
426 427 428 429 430 431
      // 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;
432

Francois Gygi committed
433 434 435
      for ( int ig = 0; ig < ngwloc; ig++ )
      {
        const double psi2s = psi2sum[ig];
436

Francois Gygi committed
437 438
        // tsum[0]: ekin partial sum
        tsum[0] += psi2s * kpg2[ig];
439

Francois Gygi committed
440 441 442 443 444
        if ( compute_stress )
        {
          const double tkpgx = kpg_x[ig];
          const double tkpgy = kpg_y[ig];
          const double tkpgz = kpg_z[ig];
445

Francois Gygi committed
446
          const double fac_ekin = 2.0 * psi2s;
447

Francois Gygi committed
448 449 450 451 452 453
          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;
454

Francois Gygi committed
455 456 457 458
        }
        // tsum[0-6] contains the contributions to
        // ekin, sigma_ekin, from vector ig
      } // ig
459

Francois Gygi committed
460 461 462 463 464
      if ( use_confinement )
      {
        const valarray<double>& fstress = cfp[ikp]->fstress();
        const valarray<double>& dfstress = cfp[ikp]->dfstress();
        for ( int ig = 0; ig < ngwloc; ig++ )
465
        {
Francois Gygi committed
466 467 468
          const double psi2s = psi2sum[ig];
          // tsum[7]: econf partial sum
          tsum[7] += psi2s * fstress[ig];
469

Francois Gygi committed
470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486
          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
487
      }
Francois Gygi committed
488 489

      sum += weight * tsum;
490 491
    } // ikp
  } // ispin
492

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

496
  ekin_  = sum[0];
497 498 499 500 501 502
  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];
503

504
  econf_ = sum[7];
505 506 507 508 509 510
  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];
511

512 513
  sigma_ekin *= omega_inv;
  sigma_econf *= omega_inv;
514

Francois Gygi committed
515
  tmap["ekin"].stop();
516

517 518 519 520
  // Stress from Eps
  sigma_eps = 0.0;
  if ( compute_stress )
  {
521
    sum = 0.0;
522 523 524 525 526 527 528 529
    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)
530
      const double fac = 2.0 * gi[ig] *
531
        real( conj(rhoelg[ig]) * dvion_local_g[ig] );
532

533 534 535
      const double tgx = g_x[ig];
      const double tgy = g_y[ig];
      const double tgz = g_z[ig];
536

537 538 539 540 541 542
      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;
543
    }
544
    MPI_Allreduce(&sum[0],&tsum[0],6,MPI_DOUBLE,MPI_SUM,vbasis_->comm());
Francois Gygi committed
545

Francois Gygi committed
546 547 548 549 550 551
    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
552
  }
553

554 555 556
  // Stress from Hartree energy
  if ( compute_stress )
  {
557
    sum = 0.0;
558 559 560
    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);
561

562 563 564 565 566 567 568
    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];

569 570 571 572 573 574
      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;
575 576 577 578 579 580 581 582 583 584 585 586
    }

    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];
587 588 589
        Species *s = s_.atoms.species_list[is];
        const double g = vbasis_->g_ptr()[ig];
        const double gi = vbasis_->gi_ptr()[ig];
590
        // next line: keep only real part
591 592 593 594 595 596 597 598 599 600 601 602
        // 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;
        }
603

604 605 606 607
        const double tgx = g_x[ig];
        const double tgy = g_y[ig];
        const double tgz = g_z[ig];

608 609 610 611 612 613
        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;
614 615
      }
    }
616
    MPI_Allreduce(&sum[0],&tsum[0],6,MPI_DOUBLE,MPI_SUM,vbasis_->comm());
Francois Gygi committed
617 618 619 620 621 622 623 624 625 626
    // 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];
627
  } // compute_stress
628

629
  // Non local energy and forces
Francois Gygi committed
630
  tmap["nonlocal"].start();
Francois Gygi committed
631 632
  // modify next loop to span only local ikp
  enl_ = 0.0;
633 634 635 636 637 638
  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
639
  for ( int ikp = 0; ikp < wf.nkp(); ikp++ )
640
  {
Francois Gygi committed
641
    for ( int ispin = 0; ispin < nlp.size(); ispin++ )
642
    {
Francois Gygi committed
643 644 645
      enl_ += wf.weight(ikp) * nlp[ispin][ikp]->energy(compute_hpsi,
              *dwf.sd(ispin,ikp), compute_forces, fion_enl, compute_stress,
              sigma_enl_kp);
646

647 648 649 650
      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];
651

652 653 654
      if ( compute_stress )
        sigma_enl += wf.weight(ikp) * sigma_enl_kp;
    }
655
  }
Francois Gygi committed
656 657 658
  tmap["nonlocal"].stop();

  ecoul_ = ehart_ + esr_ - eself_;
659 660 661 662 663 664 665
  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;
  }
666
  etotal_ = ekin_ + econf_ + eps_ + enl_ + ecoul_ + exc_ + ets_ + eexf_;
667
  enthalpy_ = etotal_;
668 669 670 671 672 673

  // Electric enthalpy
  eefield_ = 0.0;
  if ( el_enth_ )
  {
    tmap["el_enth_energy"].start();
674
    eefield_ = el_enth_->enthalpy(dwf,compute_hpsi);
675
    tmap["el_enth_energy"].stop();
676
    enthalpy_ += eefield_;
677 678 679 680 681 682 683 684 685 686 687 688 689

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

691
  epv_ = 0.0;
692 693 694 695 696 697 698
  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_;
  }
699

Francois Gygi committed
700 701 702 703 704 705 706
  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
707 708 709 710 711 712 713 714 715 716 717
        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
718
        {
Francois Gygi committed
719
          for ( int n = 0; n < sd.nstloc(); n++ )
Francois Gygi committed
720
          {
Francois Gygi committed
721 722
            const valarray<double>& fstress = cfp[ikp]->fstress();
            for ( int ig = 0; ig < ngwloc; ig++ )
723
            {
Francois Gygi committed
724 725
              cp[ig+mloc*n] += 0.5 * ( kpg2[ig] + fstress[ig] ) *
                               c[ig+mloc*n];
726
            }
727
          }
Francois Gygi committed
728 729 730 731
        }
        else
        {
          for ( int n = 0; n < sd.nstloc(); n++ )
732
          {
Francois Gygi committed
733
            for ( int ig = 0; ig < ngwloc; ig++ )
Francois Gygi committed
734
            {
Francois Gygi committed
735
              cp[ig+mloc*n] += 0.5 * kpg2[ig] * c[ig+mloc*n];
Francois Gygi committed
736 737
            }
          }
Francois Gygi committed
738
        }
Francois Gygi committed
739 740

        // local potential
Francois Gygi committed
741
        sd.rs_mul_add(*ft[ikp], &v_r[ispin][0], sdp);
Francois Gygi committed
742 743 744 745
      } //ikp
    } //ispin

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

Francois Gygi committed
748 749
    tmap["hpsi"].stop();
  } // if compute_hpsi
750

Francois Gygi committed
751 752 753 754 755 756
  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);
757
    valarray<double> fion_nl, fion_nl_tmp;
758

Francois Gygi committed
759 760 761 762 763 764
    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]);
765 766
        if ( core_charge_ )
          vtemp[ig] += conj(vxc_g[ig]) * rhocore_sp_g[is][ig];
Francois Gygi committed
767
      }
768 769 770
      fion_nl.resize(3*na_[is]);
      fion_nl = 0.0;
      fion_nl_tmp.resize(3*na_[is]);
Francois Gygi committed
771 772 773 774 775 776 777 778 779 780 781 782 783
      // 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
784 785 786
          // complex<double> teigr = ei0[kv[3*ig+0]] *
          //                         ei1[kv[3*ig+1]] *
          //                         ei2[kv[3*ig+2]];
Francois Gygi committed
787 788 789 790
          const int iii = ig+ig+ig;
          const int kx = idx[iii];
          const int ky = idx[iii+1];
          const int kz = idx[iii+2];
791

Francois Gygi committed
792 793 794
          const double cos_a = c0[kx];
          const double cos_b = c1[ky];
          const double cos_c = c2[kz];
795

Francois Gygi committed
796 797 798
          const double sin_a = s0[kx];
          const double sin_b = s1[ky];
          const double sin_c = s2[kz];
799

Francois Gygi committed
800 801
          // Next line: exp(-i*gr) =
          // (cos_a - I sin_a)*(cos_b - I sin_b)*(cos_c - I sin_c)
802
          double teigr_re =
Francois Gygi committed
803 804
            cos_a*cos_b*cos_c - sin_a*sin_b*cos_c -
            sin_a*cos_b*sin_c - cos_a*sin_b*sin_c;
805
          double teigr_im =
Francois Gygi committed
806 807
            sin_a*sin_b*sin_c - sin_a*cos_b*cos_c -
            cos_a*sin_b*cos_c - cos_a*cos_b*sin_c;
808

Francois Gygi committed
809
          /* fion is real */
810
          double tmp = teigr_re * vtemp[ig].imag() +
Francois Gygi committed
811 812 813 814 815 816
                       teigr_im * vtemp[ig].real();

          sum0 += tmp * gx0[ig];
          sum1 += tmp * gx1[ig];
          sum2 += tmp * gx2[ig];
        }
817 818 819
        fion_nl[3*ia]   = sum0;
        fion_nl[3*ia+1] = sum1;
        fion_nl[3*ia+2] = sum2;
Francois Gygi committed
820 821

      }
822

823 824
      MPI_Allreduce(&fion_nl[0],&fion_nl_tmp[0],3*na_[is],
                    MPI_DOUBLE,MPI_SUM,vbasis_->comm());
Francois Gygi committed
825 826 827 828

      double fac = -2.0 * omega;
      for ( int i = 0; i < 3*na_[is]; i++ )
      {
829
        fion[is][i] += fion_esr[is][i] + fac * fion_nl_tmp[i];
Francois Gygi committed
830
      }
Francois Gygi committed
831 832 833 834 835 836 837 838 839

      // 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
840 841
    }
  }
842

843
  if ( compute_stress )
844
  {
845 846
    sigma = sigma_ekin + sigma_econf + sigma_eps + sigma_enl +
            sigma_ehart + sigma_exc + sigma_esr;
847
  }
848

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

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

965
  // fill tau0 with values in atom_list
966

Francois Gygi committed
967 968
  atoms.get_positions(tau0);
  sf.update(tau0,*vbasis_);
969

Francois Gygi committed
970 971
  // compute Fourier coefficients of the local potential
  memset( (void*)&vion_local_g[0], 0, 2*ngloc*sizeof(double) );
972
  memset( (void*)&dvion_local_g[0], 0, 2*ngloc*sizeof(double) );
Francois Gygi committed
973
  memset( (void*)&rhopst[0], 0, 2*ngloc*sizeof(double) );
974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996

  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
997 998 999 1000 1001
  for ( int is = 0; is < atoms.nsp(); is++ )
  {
    complex<double> *s = &sf.sfac[is][0];
    for ( int ig = 0; ig < ngloc; ig++ )
    {
1002 1003 1004 1005
      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
1006 1007
    }
  }
1008

Francois Gygi committed
1009 1010
  // compute esr: pseudocharge repulsion energy
  const UnitCell& cell = s_.wf.cell();
1011
  const double omega_inv = 1.0 / cell.volume();
Francois Gygi committed
1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023
  // 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
1024

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

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

Francois Gygi committed
1110 1111 1112 1113 1114
}

////////////////////////////////////////////////////////////////////////////////
void EnergyFunctional::cell_moved(void)
{
Francois Gygi committed
1115 1116
  const Wavefunction& wf = s_.wf;
  const UnitCell& cell = wf.cell();
1117
  // resize vbasis_
Francois Gygi committed
1118
  vbasis_->resize(cell,s_.wf.refcell(),4.0*s_.wf.ecut());
1119

1120 1121 1122
  const int ngloc = vbasis_->localsize();
  const double omega = cell.volume();
  assert(omega != 0.0);
1123 1124
  const double omega_inv = 1.0 / omega;

1125 1126
  const AtomSet& atoms = s_.atoms;
  for ( int is = 0; is < nsp_; is++ )
1127
  {
1128 1129
    Species *s = atoms.species_list[is];
    const double * const g = vbasis_->g_ptr();
1130
    double v,dv,rhoc;
1131
    for ( int ig = 0; ig < ngloc; ig++ )
1132
    {
1133 1134 1135 1136
      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;
1137 1138 1139 1140 1141
      if ( core_charge_ )
      {
        s->rhocoreg(g[ig],rhoc);
        rhocore_sp_g[is][ig] = rhoc * omega_inv;
      }
1142
    }
1143
  }
1144

Francois Gygi committed
1145
  // Update confinement potentials and non-local potentials
1146 1147
  for ( int ikp = 0; ikp < wf.nkp(); ikp++ )
  {
Francois Gygi committed
1148
    cfp[ikp]->update();
Francois Gygi committed
1149 1150
    for ( int ispin = 0; ispin < nlp.size(); ispin++ )
      nlp[ispin][ikp]->update_twnl();
1151
  }
1152 1153 1154

  // Update exchange-correlation operator
  xco->cell_moved();
Francois Gygi committed
1155
}
1156 1157 1158 1159 1160 1161 1162

////////////////////////////////////////////////////////////////////////////////
void EnergyFunctional::print(ostream& os) const
{
  os.setf(ios::fixed,ios::floatfield);
  os.setf(ios::right,ios::adjustfield);
  os << setprecision(8);
1163 1164 1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176
  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;
1177
}
1178

1179 1180
////////////////////////////////////////////////////////////////////////////////
ostream& operator<< ( ostream& os, const EnergyFunctional& e )
1181 1182
{
  e.print(os);
1183 1184
  return os;
}