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

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

////////////////////////////////////////////////////////////////////////////////
32 33
XCPotential::XCPotential(const ChargeDensity& cd, const string functional_name,
  const Control& ctrl): cd_(cd), vft_(*cd_.vft()), vbasis_(*cd_.vbasis())
Francois Gygi committed
34 35 36 37 38
{
  if ( functional_name == "LDA" )
  {
    xcf_ = new LDAFunctional(cd_.rhor);
  }
Francois Gygi committed
39 40 41 42
  else if ( functional_name == "VWN" )
  {
    xcf_ = new VWNFunctional(cd_.rhor);
  }
Francois Gygi committed
43 44 45 46 47 48 49 50
  else if ( functional_name == "PBE" )
  {
    xcf_ = new PBEFunctional(cd_.rhor);
  }
  else if ( functional_name == "BLYP" )
  {
    xcf_ = new BLYPFunctional(cd_.rhor);
  }
Francois Gygi committed
51 52
  else if ( functional_name == "PBE0" )
  {
53
    const double x_coeff = 1.0 - ctrl.alpha_PBE0;
Francois Gygi committed
54 55 56 57 58 59 60
    const double c_coeff = 1.0;
    xcf_ = new PBEFunctional(cd_.rhor,x_coeff,c_coeff);
  }
  else if ( functional_name == "B3LYP" )
  {
    xcf_ = new B3LYPFunctional(cd_.rhor);
  }
Francois Gygi committed
61 62 63 64 65 66 67
  else
  {
    throw XCPotentialException("unknown functional name");
  }
  nspin_ = cd_.rhor.size();
  ngloc_ = vbasis_.localsize();
  np012loc_ = vft_.np012loc();
68

Francois Gygi committed
69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86
  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
87 88 89 90 91 92
////////////////////////////////////////////////////////////////////////////////
bool XCPotential::isGGA(void)
{
  return xcf_->isGGA();
}

Francois Gygi committed
93
////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
94
void XCPotential::update(vector<vector<double> >& vr)
Francois Gygi committed
95 96
{
  // compute exchange-correlation energy and add vxc potential to vr[ispin][ir]
97

Francois Gygi committed
98 99 100 101 102
  // 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
103

Francois Gygi committed
104 105
  // Output: (through member function xcf())
  //
106
  // exc_, dxc_
Francois Gygi committed
107 108
  //
  // LDA Functional:
109
  //   exc_, dxc_
Francois Gygi committed
110 111 112 113
  //   spin unpolarized: xcf()->exc, xcf()->vxc1
  //   spin polarized:   xcf()->exc, xcf()->vxc1_up, xcf()->vxc1_dn
  //
  // GGA Functional: (through member function xcf())
114
  //   exc_, dxc_
Francois Gygi committed
115
  //   spin unpolarized: xcf()->exc, xcf()->vxc1, xcf()->vxc2
116
  //   spin polarized:   xcf()->exc_up, xcf()->exc_dn,
Francois Gygi committed
117
  //                     xcf()->vxc1_up, xcf()->vxc1_dn
118
  //                     xcf()->vxc2_upup, xcf()->vxc2_dndn,
Francois Gygi committed
119
  //                     xcf()->vxc2_updn, xcf()->vxc2_dnup
120

Francois Gygi committed
121 122 123
  if ( !xcf_->isGGA() )
  {
    // LDA functional
124

Francois Gygi committed
125
    xcf_->setxc();
126

Francois Gygi committed
127
    exc_ = 0.0;
128
    dxc_ = 0.0;
Francois Gygi committed
129
    const double *const e = xcf_->exc;
Francois Gygi committed
130
    const int size = xcf_->np();
131

Francois Gygi committed
132 133 134 135 136 137 138
    if ( nspin_ == 1 )
    {
      // unpolarized
      const double *const rh = xcf_->rho;
      const double *const v = xcf_->vxc1;
      for ( int i = 0; i < size; i++ )
      {
139 140 141 142 143 144
        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
145 146 147 148 149 150 151 152 153
      }
    }
    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;
154 155
      for ( int i = 0; i < size; i++ )
      {
156 157 158
        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
159 160
        vr[0][i] += v_up[i];
        vr[1][i] += v_dn[i];
161
      }
Francois Gygi committed
162
    }
163 164 165 166 167 168
    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
169 170 171 172 173
  }
  else
  {
    // GGA functional
    exc_ = 0.0;
174

Francois Gygi committed
175 176 177 178 179 180 181 182 183 184 185 186 187 188 189
    // 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
190
        dcopy(&np012loc_,(double*)&tmpr[0],&inc2,grj,&inc1);
Francois Gygi committed
191 192 193 194 195 196 197 198 199 200 201 202 203 204 205
      }
    }
    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++;
206 207
          tmp1[ig] = igxj * c0;
          tmp2[ig] = igxj * c1;
Francois Gygi committed
208 209 210 211 212 213
        }
        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
214 215
        dcopy(&np012loc_,p,  &inc2,grj_up,&inc1);
        dcopy(&np012loc_,p+1,&inc2,grj_dn,&inc1);
Francois Gygi committed
216 217
      } // j
    }
218

Francois Gygi committed
219
    xcf_->setxc();
220

Francois Gygi committed
221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249
    // 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
250
          dcopy(&np012loc_,(double*)&tmpr[0],&inc2,&vxctmp[0][0],&inc1);
Francois Gygi committed
251 252 253
        }
        else
        {
Francois Gygi committed
254
          daxpy(&np012loc_,&one,(double*)&tmpr[0],&inc2,&vxctmp[0][0],&inc1);
Francois Gygi committed
255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290
        }
      }
    }
    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
291 292
          dcopy(&np012loc_,p  ,&inc2,&vxctmp[0][0],&inc1);
          dcopy(&np012loc_,p+1,&inc2,&vxctmp[1][0],&inc1);
Francois Gygi committed
293 294 295
        }
        else
        {
Francois Gygi committed
296 297
          daxpy(&np012loc_,&one,p  ,&inc2,&vxctmp[0][0],&inc1);
          daxpy(&np012loc_,&one,p+1,&inc2,&vxctmp[1][0],&inc1);
Francois Gygi committed
298 299 300
        }
      } // j
    }
301

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

Francois Gygi committed
305
    double esum=0.0;
306
    double dsum=0.0;
Francois Gygi committed
307 308 309 310 311 312 313 314
    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++ )
        {
315 316 317 318 319 320
          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
321 322 323 324 325 326 327 328 329 330 331
        }
      }
    }
    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;
332
      for ( int ir = 0; ir < np012loc_; ir++ )
Francois Gygi committed
333
      {
334 335 336 337 338 339 340 341
        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
342 343
      }
    }
344 345 346 347 348 349
    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
350 351 352 353 354 355
  }
}
////////////////////////////////////////////////////////////////////////////////
void XCPotential::compute_stress(valarray<double>& sigma_exc)
{
  // compute exchange-correlation contributions to the stress tensor
356

Francois Gygi committed
357 358 359
  if ( !xcf_->isGGA() )
  {
    // LDA functional
360

361
    double dsum = 0.0;
Francois Gygi committed
362
    const double *const e = xcf_->exc;
Francois Gygi committed
363
    const int size = xcf_->np();
364

Francois Gygi committed
365 366 367 368 369 370 371
    if ( nspin_ == 1 )
    {
      // unpolarized
      const double *const rh = xcf_->rho;
      const double *const v = xcf_->vxc1;
      for ( int i = 0; i < size; i++ )
      {
372
        dsum += rh[i] * (e[i] - v[i]);
Francois Gygi committed
373 374 375 376 377 378 379 380 381
      }
    }
    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;
382 383 384
      for ( int i = 0; i < size; i++ )
      {
        const double rh = rh_up[i] + rh_dn[i];
385
        dsum += rh * e[i] - rh_up[i] * v_up[i] - rh_dn[i] * v_dn[i];
386
      }
Francois Gygi committed
387 388
    }
    const double fac = 1.0 / vft_.np012();
389
    double sum, tsum;
390
    // Next line: factor omega in volume element cancels 1/omega in
Francois Gygi committed
391
    // definition of sigma_exc
392
    sum = - fac * dsum;
393
    MPI_Allreduce(&sum,&tsum,1,MPI_DOUBLE,MPI_SUM,vbasis_.comm());
394

Francois Gygi committed
395 396 397 398 399 400 401 402 403 404 405
    // 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
406

Francois Gygi committed
407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 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
    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] );
452

Francois Gygi committed
453 454 455 456 457 458 459
        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;
460

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

Francois Gygi committed
469 470 471
        const double grad_up_grad_dn = grx_up * grx_dn +
                                       gry_up * gry_dn +
                                       grz_up * grz_dn;
Francois Gygi committed
472

Francois Gygi committed
473 474 475 476
        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];
477

Francois Gygi committed
478 479 480 481
        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 );
482

Francois Gygi committed
483 484 485 486
        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 );
487

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

Francois Gygi committed
493 494 495 496
        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;
497

Francois Gygi committed
498 499 500 501
        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;
502

Francois Gygi committed
503 504 505 506
        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
507 508
      }
    }
509
    double fac = 1.0 / vft_.np012();
510
    double sum[6],tsum[6];
511
    // Next line: factor omega in volume element cancels 1/omega in
512
    // definition of sigma_exc
513 514 515 516 517 518 519
    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());
520

Francois Gygi committed
521 522 523 524 525 526
    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
527 528
  }
}