XCPotential.C 16.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
// XCPotential.C
//
////////////////////////////////////////////////////////////////////////////////

#include "XCPotential.h"
20
#include "LDAFunctional.h"
Francois Gygi committed
21
#include "VWNFunctional.h"
22 23
#include "PBEFunctional.h"
#include "BLYPFunctional.h"
24
#include "HSEFunctional.h"
Francois Gygi committed
25
#include "RSHFunctional.h"
Francois Gygi committed
26
#include "B3LYPFunctional.h"
27
#include "BHandHLYPFunctional.h"
Francois Gygi committed
28 29 30 31 32 33 34
#include "Basis.h"
#include "FourierTransform.h"
#include "blas.h" // daxpy, dcopy
#include <cassert>
using namespace std;

////////////////////////////////////////////////////////////////////////////////
35 36
XCPotential::XCPotential(const ChargeDensity& cd, const string functional_name,
  const Control& ctrl): cd_(cd), vft_(*cd_.vft()), vbasis_(*cd_.vbasis())
Francois Gygi committed
37 38 39 40 41
{
  if ( functional_name == "LDA" )
  {
    xcf_ = new LDAFunctional(cd_.rhor);
  }
Francois Gygi committed
42 43 44 45
  else if ( functional_name == "VWN" )
  {
    xcf_ = new VWNFunctional(cd_.rhor);
  }
Francois Gygi committed
46 47 48 49 50 51 52 53
  else if ( functional_name == "PBE" )
  {
    xcf_ = new PBEFunctional(cd_.rhor);
  }
  else if ( functional_name == "BLYP" )
  {
    xcf_ = new BLYPFunctional(cd_.rhor);
  }
Francois Gygi committed
54 55
  else if ( functional_name == "PBE0" )
  {
56
    const double x_coeff = 1.0 - ctrl.alpha_PBE0;
Francois Gygi committed
57 58 59
    const double c_coeff = 1.0;
    xcf_ = new PBEFunctional(cd_.rhor,x_coeff,c_coeff);
  }
60 61 62 63
  else if ( functional_name == "HSE" )
  {
    xcf_ = new HSEFunctional(cd_.rhor);
  }
Francois Gygi committed
64 65 66 67
  else if ( functional_name == "RSH" )
  {
    xcf_ = new RSHFunctional(cd_.rhor,ctrl.alpha_RSH,ctrl.beta_RSH,ctrl.mu_RSH);
  }
Francois Gygi committed
68 69 70 71
  else if ( functional_name == "B3LYP" )
  {
    xcf_ = new B3LYPFunctional(cd_.rhor);
  }
72 73 74 75
  else if ( functional_name == "BHandHLYP" )
  {
    xcf_ = new BHandHLYPFunctional(cd_.rhor);
  }
Francois Gygi committed
76 77 78 79 80 81 82
  else
  {
    throw XCPotentialException("unknown functional name");
  }
  nspin_ = cd_.rhor.size();
  ngloc_ = vbasis_.localsize();
  np012loc_ = vft_.np012loc();
83

Francois Gygi committed
84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101
  if ( xcf_->isGGA() )
  {
    tmp1.resize(ngloc_);
    if ( nspin_ > 1 )
      tmp2.resize(ngloc_);
    vxctmp.resize(nspin_);
    for ( int ispin = 0; ispin < nspin_; ispin++ )
      vxctmp[ispin].resize(np012loc_);
    tmpr.resize(np012loc_);
  }
}

////////////////////////////////////////////////////////////////////////////////
XCPotential::~XCPotential(void)
{
  delete xcf_;
}

Francois Gygi committed
102 103 104 105 106 107
////////////////////////////////////////////////////////////////////////////////
bool XCPotential::isGGA(void)
{
  return xcf_->isGGA();
}

Francois Gygi committed
108
////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
109
void XCPotential::update(vector<vector<double> >& vr)
Francois Gygi committed
110 111
{
  // compute exchange-correlation energy and add vxc potential to vr[ispin][ir]
112

Francois Gygi committed
113 114 115 116 117
  // Input: total electronic density in:
  //   vector<vector<double> >           cd_.rhor[ispin][ir] (real space)
  //   vector<vector<complex<double> > > cd_.rhog[ispin][ig] (Fourier coeffs)
  // The array cd_.rhog is only used if xcf->isGGA() == true
  // to compute the density gradients
118

Francois Gygi committed
119 120
  // Output: (through member function xcf())
  //
121
  // exc_, dxc_
Francois Gygi committed
122 123
  //
  // LDA Functional:
124
  //   exc_, dxc_
Francois Gygi committed
125 126 127 128
  //   spin unpolarized: xcf()->exc, xcf()->vxc1
  //   spin polarized:   xcf()->exc, xcf()->vxc1_up, xcf()->vxc1_dn
  //
  // GGA Functional: (through member function xcf())
129
  //   exc_, dxc_
Francois Gygi committed
130
  //   spin unpolarized: xcf()->exc, xcf()->vxc1, xcf()->vxc2
131
  //   spin polarized:   xcf()->exc_up, xcf()->exc_dn,
Francois Gygi committed
132
  //                     xcf()->vxc1_up, xcf()->vxc1_dn
133
  //                     xcf()->vxc2_upup, xcf()->vxc2_dndn,
Francois Gygi committed
134
  //                     xcf()->vxc2_updn, xcf()->vxc2_dnup
135

Francois Gygi committed
136 137 138
  if ( !xcf_->isGGA() )
  {
    // LDA functional
139

Francois Gygi committed
140
    xcf_->setxc();
141

Francois Gygi committed
142
    exc_ = 0.0;
143
    dxc_ = 0.0;
Francois Gygi committed
144
    const double *const e = xcf_->exc;
Francois Gygi committed
145
    const int size = xcf_->np();
146

Francois Gygi committed
147 148 149 150 151 152 153
    if ( nspin_ == 1 )
    {
      // unpolarized
      const double *const rh = xcf_->rho;
      const double *const v = xcf_->vxc1;
      for ( int i = 0; i < size; i++ )
      {
154 155 156 157 158 159
        const double e_i = e[i];
        const double v_i = v[i];
        const double rh_i = rh[i];
        exc_ += rh_i * e_i;
        dxc_ += rh_i * ( e_i - v_i );
        vr[0][i] += v_i;
Francois Gygi committed
160 161 162 163 164 165 166 167 168
      }
    }
    else
    {
      // spin polarized
      const double *const rh_up = xcf_->rho_up;
      const double *const rh_dn = xcf_->rho_dn;
      const double *const v_up = xcf_->vxc1_up;
      const double *const v_dn = xcf_->vxc1_dn;
169 170
      for ( int i = 0; i < size; i++ )
      {
171 172 173
        const double r_i = rh_up[i] + rh_dn[i];
        exc_ += r_i * e[i];
        dxc_ += r_i * e[i] - rh_up[i] * v_up[i] - rh_dn[i] * v_dn[i];
Francois Gygi committed
174 175
        vr[0][i] += v_up[i];
        vr[1][i] += v_dn[i];
176
      }
Francois Gygi committed
177
    }
178 179 180 181 182 183
    double sum[2],tsum[2];
    sum[0] = exc_ * vbasis_.cell().volume() / vft_.np012();
    sum[1] = dxc_ * vbasis_.cell().volume() / vft_.np012();
    MPI_Allreduce(&sum,&tsum,2,MPI_DOUBLE,MPI_SUM,vbasis_.comm());
    exc_ = tsum[0];
    dxc_ = tsum[1];
Francois Gygi committed
184 185 186 187 188
  }
  else
  {
    // GGA functional
    exc_ = 0.0;
189

Francois Gygi committed
190 191 192 193 194 195 196 197 198 199 200 201 202 203 204
    // compute grad_rho
    const double omega_inv = 1.0 / vbasis_.cell().volume();
    if ( nspin_ == 1 )
    {
      for ( int j = 0; j < 3; j++ )
      {
        const double *const gxj = vbasis_.gx_ptr(j);
        for ( int ig = 0; ig < ngloc_; ig++ )
        {
          /* i*G_j*c(G) */
          tmp1[ig] = complex<double>(0.0,omega_inv*gxj[ig]) * cd_.rhog[0][ig];
        }
        vft_.backward(&tmp1[0],&tmpr[0]);
        int inc2=2, inc1=1;
        double *grj = xcf_->grad_rho[j];
Francois Gygi committed
205
        dcopy(&np012loc_,(double*)&tmpr[0],&inc2,grj,&inc1);
Francois Gygi committed
206 207 208 209 210 211 212 213 214 215 216 217 218 219 220
      }
    }
    else
    {
      for ( int j = 0; j < 3; j++ )
      {
        const double *const gxj = vbasis_.gx_ptr(j);
        const complex<double>* rhg0 = &cd_.rhog[0][0];
        const complex<double>* rhg1 = &cd_.rhog[1][0];
        for ( int ig = 0; ig < ngloc_; ig++ )
        {
          /* i*G_j*c(G) */
          const complex<double> igxj(0.0,omega_inv*gxj[ig]);
          const complex<double> c0 = *rhg0++;
          const complex<double> c1 = *rhg1++;
221 222
          tmp1[ig] = igxj * c0;
          tmp2[ig] = igxj * c1;
Francois Gygi committed
223 224 225 226 227 228
        }
        vft_.backward(&tmp1[0],&tmp2[0],&tmpr[0]);
        double *grj_up = xcf_->grad_rho_up[j];
        double *grj_dn = xcf_->grad_rho_dn[j];
        int inc2=2, inc1=1;
        double* p = (double*) &tmpr[0];
Francois Gygi committed
229 230
        dcopy(&np012loc_,p,  &inc2,grj_up,&inc1);
        dcopy(&np012loc_,p+1,&inc2,grj_dn,&inc1);
Francois Gygi committed
231 232
      } // j
    }
233

Francois Gygi committed
234
    xcf_->setxc();
235

Francois Gygi committed
236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264
    // compute xc potential
    // take divergence of grad(rho)*vxc2

    // compute components of grad(rho) * vxc2
    if ( nspin_ == 1 )
    {
      for ( int j = 0; j < 3; j++ )
      {
        const double *const gxj = vbasis_.gx_ptr(j);
        const double *const grj = xcf_->grad_rho[j];
        const double *const v2 = xcf_->vxc2;
        for ( int ir = 0; ir < np012loc_; ir++ )
        {
          tmpr[ir] = grj[ir] * v2[ir];
        }
        // derivative
        vft_.forward(&tmpr[0],&tmp1[0]);
        for ( int ig = 0; ig < ngloc_; ig++ )
        {
          // i*G_j*c(G)
          tmp1[ig] *= complex<double>(0.0,gxj[ig]);
        }
        // back to real space
        vft_.backward(&tmp1[0],&tmpr[0]);
        // accumulate div(vxc2*grad_rho) in vxctmp
        double one = 1.0;
        int inc1 = 1, inc2 = 2;
        if ( j == 0 )
        {
Francois Gygi committed
265
          dcopy(&np012loc_,(double*)&tmpr[0],&inc2,&vxctmp[0][0],&inc1);
Francois Gygi committed
266 267 268
        }
        else
        {
Francois Gygi committed
269
          daxpy(&np012loc_,&one,(double*)&tmpr[0],&inc2,&vxctmp[0][0],&inc1);
Francois Gygi committed
270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305
        }
      }
    }
    else
    {
      double *v2_upup = xcf_->vxc2_upup;
      double *v2_updn = xcf_->vxc2_updn;
      double *v2_dnup = xcf_->vxc2_dnup;
      double *v2_dndn = xcf_->vxc2_dndn;
      for ( int j = 0; j < 3; j++ )
      {
        const double *gxj = vbasis_.gx_ptr(j);
        const double *grj_up = xcf_->grad_rho_up[j];
        const double *grj_dn = xcf_->grad_rho_dn[j];
        for ( int ir = 0; ir < np012loc_; ir++ )
        {
          const double re = v2_upup[ir] * grj_up[ir] + v2_updn[ir] * grj_dn[ir];
          const double im = v2_dnup[ir] * grj_up[ir] + v2_dndn[ir] * grj_dn[ir];
          tmpr[ir] = complex<double>(re,im);
        }
        // derivative
        vft_.forward(&tmpr[0],&tmp1[0],&tmp2[0]);
        for ( int ig = 0; ig < ngloc_; ig++ )
        {
          // i*G_j*c(G)
          const complex<double> igxj(0.0,gxj[ig]);
          tmp1[ig] *= igxj;
          tmp2[ig] *= igxj;
        }
        vft_.backward(&tmp1[0],&tmp2[0],&tmpr[0]);
        // accumulate div(vxc2*grad_rho) in vxctmp
        double one = 1.0;
        int inc1 = 1, inc2 = 2;
        double* p = (double*) &tmpr[0];
        if ( j == 0 )
        {
Francois Gygi committed
306 307
          dcopy(&np012loc_,p  ,&inc2,&vxctmp[0][0],&inc1);
          dcopy(&np012loc_,p+1,&inc2,&vxctmp[1][0],&inc1);
Francois Gygi committed
308 309 310
        }
        else
        {
Francois Gygi committed
311 312
          daxpy(&np012loc_,&one,p  ,&inc2,&vxctmp[0][0],&inc1);
          daxpy(&np012loc_,&one,p+1,&inc2,&vxctmp[1][0],&inc1);
Francois Gygi committed
313 314 315
        }
      } // j
    }
316

Francois Gygi committed
317 318 319
    // add xc potential to local potential in vr[i]
    // div(vxc2*grad_rho) is stored in vxctmp[ispin][ir]

Francois Gygi committed
320
    double esum=0.0;
321
    double dsum=0.0;
Francois Gygi committed
322 323 324 325 326 327 328 329
    if ( nspin_ == 1 )
    {
      const double *const e = xcf_->exc;
      const double *const v1 = xcf_->vxc1;
      const double *const rh = xcf_->rho;
      {
        for ( int ir = 0; ir < np012loc_; ir++ )
        {
330 331 332 333 334 335
          const double e_i = e[ir];
          const double rh_i = rh[ir];
          const double v_i = v1[ir] + vxctmp[0][ir];
          esum += rh_i * e_i;
          dsum += rh_i * ( e_i - v_i );
          vr[0][ir] += v_i;
Francois Gygi committed
336 337 338 339 340 341 342 343 344 345 346
        }
      }
    }
    else
    {
      const double *const v1_up = xcf_->vxc1_up;
      const double *const v1_dn = xcf_->vxc1_dn;
      const double *const eup = xcf_->exc_up;
      const double *const edn = xcf_->exc_dn;
      const double *const rh_up = xcf_->rho_up;
      const double *const rh_dn = xcf_->rho_dn;
347
      for ( int ir = 0; ir < np012loc_; ir++ )
Francois Gygi committed
348
      {
349 350 351 352 353 354 355 356
        const double r_up_i = rh_up[ir];
        const double r_dn_i = rh_dn[ir];
        esum += r_up_i * eup[ir] + r_dn_i * edn[ir];
        const double v_up = v1_up[ir] + vxctmp[0][ir];
        const double v_dn = v1_dn[ir] + vxctmp[1][ir];
        dsum += r_up_i * ( eup[ir] - v_up ) + r_dn_i * ( edn[ir] - v_dn );
        vr[0][ir] += v_up;
        vr[1][ir] += v_dn;
Francois Gygi committed
357 358
      }
    }
359 360 361 362 363 364
    double sum[2], tsum[2];
    sum[0] = esum * vbasis_.cell().volume() / vft_.np012();
    sum[1] = dsum * vbasis_.cell().volume() / vft_.np012();
    MPI_Allreduce(&sum,&tsum,2,MPI_DOUBLE,MPI_SUM,vbasis_.comm());
    exc_ = tsum[0];
    dxc_ = tsum[1];
Francois Gygi committed
365 366 367 368 369 370
  }
}
////////////////////////////////////////////////////////////////////////////////
void XCPotential::compute_stress(valarray<double>& sigma_exc)
{
  // compute exchange-correlation contributions to the stress tensor
371

Francois Gygi committed
372 373 374
  if ( !xcf_->isGGA() )
  {
    // LDA functional
375

376
    double dsum = 0.0;
Francois Gygi committed
377
    const double *const e = xcf_->exc;
Francois Gygi committed
378
    const int size = xcf_->np();
379

Francois Gygi committed
380 381 382 383 384 385 386
    if ( nspin_ == 1 )
    {
      // unpolarized
      const double *const rh = xcf_->rho;
      const double *const v = xcf_->vxc1;
      for ( int i = 0; i < size; i++ )
      {
387
        dsum += rh[i] * (e[i] - v[i]);
Francois Gygi committed
388 389 390 391 392 393 394 395 396
      }
    }
    else
    {
      // spin polarized
      const double *const rh_up = xcf_->rho_up;
      const double *const rh_dn = xcf_->rho_dn;
      const double *const v_up = xcf_->vxc1_up;
      const double *const v_dn = xcf_->vxc1_dn;
397 398 399
      for ( int i = 0; i < size; i++ )
      {
        const double rh = rh_up[i] + rh_dn[i];
400
        dsum += rh * e[i] - rh_up[i] * v_up[i] - rh_dn[i] * v_dn[i];
401
      }
Francois Gygi committed
402 403
    }
    const double fac = 1.0 / vft_.np012();
404
    double sum, tsum;
405
    // Next line: factor omega in volume element cancels 1/omega in
Francois Gygi committed
406
    // definition of sigma_exc
407
    sum = - fac * dsum;
408
    MPI_Allreduce(&sum,&tsum,1,MPI_DOUBLE,MPI_SUM,vbasis_.comm());
409

Francois Gygi committed
410 411 412 413 414 415 416 417 418 419 420
    // Note: contribution to sigma_exc is a multiple of the identity
    sigma_exc[0] = tsum;
    sigma_exc[1] = tsum;
    sigma_exc[2] = tsum;
    sigma_exc[3] = 0.0;
    sigma_exc[4] = 0.0;
    sigma_exc[5] = 0.0;
  }
  else
  {
    // GGA functional
421

Francois Gygi committed
422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466
    double dsum=0.0,sum0=0.0,sum1=0.0,sum2=0.0,
           sum3=0.0,sum4=0.0,sum5=0.0;
    if ( nspin_ == 1 )
    {
      const double *const e = xcf_->exc;
      const double *const v1 = xcf_->vxc1;
      const double *const v2 = xcf_->vxc2;
      const double *const rh = xcf_->rho;
      for ( int ir = 0; ir < np012loc_; ir++ )
      {
        dsum += rh[ir] * ( e[ir] - v1[ir] );
        const double grx = xcf_->grad_rho[0][ir];
        const double gry = xcf_->grad_rho[1][ir];
        const double grz = xcf_->grad_rho[2][ir];
        const double grx2 = grx * grx;
        const double gry2 = gry * gry;
        const double grz2 = grz * grz;
        const double grad2 = grx2 + gry2 + grz2;
        const double v2t = v2[ir];
        sum0 += ( grad2 + grx2 ) * v2t;
        sum1 += ( grad2 + gry2 ) * v2t;
        sum2 += ( grad2 + grz2 ) * v2t;
        sum3 += grx * gry * v2t;
        sum4 += gry * grz * v2t;
        sum5 += grx * grz * v2t;
      }
    }
    else
    {
      const double *const v1_up = xcf_->vxc1_up;
      const double *const v1_dn = xcf_->vxc1_dn;
      const double *const v2_upup = xcf_->vxc2_upup;
      const double *const v2_updn = xcf_->vxc2_updn;
      const double *const v2_dnup = xcf_->vxc2_dnup;
      const double *const v2_dndn = xcf_->vxc2_dndn;
      const double *const eup = xcf_->exc_up;
      const double *const edn = xcf_->exc_dn;
      const double *const rh_up = xcf_->rho_up;
      const double *const rh_dn = xcf_->rho_dn;
      for ( int ir = 0; ir < np012loc_; ir++ )
      {
        const double r_up = rh_up[ir];
        const double r_dn = rh_dn[ir];
        dsum += r_up * ( eup[ir] - v1_up[ir] ) +
                r_dn * ( edn[ir] - v1_dn[ir] );
467

Francois Gygi committed
468 469 470 471 472 473 474
        const double grx_up = xcf_->grad_rho_up[0][ir];
        const double gry_up = xcf_->grad_rho_up[1][ir];
        const double grz_up = xcf_->grad_rho_up[2][ir];
        const double grx2_up = grx_up * grx_up;
        const double gry2_up = gry_up * gry_up;
        const double grz2_up = grz_up * grz_up;
        const double grad2_up = grx2_up + gry2_up + grz2_up;
475

Francois Gygi committed
476 477 478 479 480 481 482
        const double grx_dn = xcf_->grad_rho_dn[0][ir];
        const double gry_dn = xcf_->grad_rho_dn[1][ir];
        const double grz_dn = xcf_->grad_rho_dn[2][ir];
        const double grx2_dn = grx_dn * grx_dn;
        const double gry2_dn = gry_dn * gry_dn;
        const double grz2_dn = grz_dn * grz_dn;
        const double grad2_dn = grx2_dn + gry2_dn + grz2_dn;
483

Francois Gygi committed
484 485 486
        const double grad_up_grad_dn = grx_up * grx_dn +
                                       gry_up * gry_dn +
                                       grz_up * grz_dn;
Francois Gygi committed
487

Francois Gygi committed
488 489 490 491
        const double v2_upup_ir = v2_upup[ir];
        const double v2_updn_ir = v2_updn[ir];
        const double v2_dnup_ir = v2_dnup[ir];
        const double v2_dndn_ir = v2_dndn[ir];
492

Francois Gygi committed
493 494 495 496
        sum0 += v2_upup_ir * ( grad2_up + grx2_up ) +
                v2_updn_ir * ( grad_up_grad_dn + grx_up * grx_dn ) +
                v2_dnup_ir * ( grad_up_grad_dn + grx_dn * grx_up ) +
                v2_dndn_ir * ( grad2_dn + grx2_dn );
497

Francois Gygi committed
498 499 500 501
        sum1 += v2_upup_ir * ( grad2_up + gry2_up ) +
                v2_updn_ir * ( grad_up_grad_dn + gry_up * gry_dn ) +
                v2_dnup_ir * ( grad_up_grad_dn + gry_dn * gry_up ) +
                v2_dndn_ir * ( grad2_dn + gry2_dn );
502

Francois Gygi committed
503 504 505 506
        sum2 += v2_upup_ir * ( grad2_up + grz2_up ) +
                v2_updn_ir * ( grad_up_grad_dn + grz_up * grz_dn ) +
                v2_dnup_ir * ( grad_up_grad_dn + grz_dn * grz_up ) +
                v2_dndn_ir * ( grad2_dn + grz2_dn );
507

Francois Gygi committed
508 509 510 511
        sum3 += v2_upup_ir * grx_up * gry_up +
                v2_updn_ir * grx_up * gry_dn +
                v2_dnup_ir * grx_dn * gry_up +
                v2_dndn_ir * grx_dn * gry_dn;
512

Francois Gygi committed
513 514 515 516
        sum4 += v2_upup_ir * gry_up * grz_up +
                v2_updn_ir * gry_up * grz_dn +
                v2_dnup_ir * gry_dn * grz_up +
                v2_dndn_ir * gry_dn * grz_dn;
517

Francois Gygi committed
518 519 520 521
        sum5 += v2_upup_ir * grx_up * grz_up +
                v2_updn_ir * grx_up * grz_dn +
                v2_dnup_ir * grx_dn * grz_up +
                v2_dndn_ir * grx_dn * grz_dn;
Francois Gygi committed
522 523
      }
    }
524
    double fac = 1.0 / vft_.np012();
525
    double sum[6],tsum[6];
526
    // Next line: factor omega in volume element cancels 1/omega in
527
    // definition of sigma_exc
528 529 530 531 532 533 534
    sum[0] = - fac * ( dsum + sum0 );
    sum[1] = - fac * ( dsum + sum1 );
    sum[2] = - fac * ( dsum + sum2 );
    sum[3] = - fac * sum3;
    sum[4] = - fac * sum4;
    sum[5] = - fac * sum5;
    MPI_Allreduce(sum,tsum,6,MPI_DOUBLE,MPI_SUM,vbasis_.comm());
535

Francois Gygi committed
536 537 538 539 540 541
    sigma_exc[0] = tsum[0];
    sigma_exc[1] = tsum[1];
    sigma_exc[2] = tsum[2];
    sigma_exc[3] = tsum[3];
    sigma_exc[4] = tsum[4];
    sigma_exc[5] = tsum[5];
Francois Gygi committed
542 543
  }
}