Commit 26449546 by Francois Gygi

merge --reintegrate of oncv branch into trunk

git-svn-id: http://qboxcode.org/svn/qb/trunk@1673 cba15fb0-1239-40c8-b417-11db7ca47a34
parent 40e559d5
......@@ -24,12 +24,12 @@
#include <iomanip>
#include <algorithm> // fill
#include <functional>
using namespace std;
////////////////////////////////////////////////////////////////////////////////
ChargeDensity::ChargeDensity(const Wavefunction& wf) : ctxt_(wf.context()),
wf_(wf), vcomm_(wf.sd(0,0)->basis().comm())
{
vbasis_ = new Basis(vcomm_, D3vector(0,0,0));
vbasis_->resize(wf.cell(),wf.refcell(),4.0*wf.ecut());
......@@ -89,6 +89,8 @@ wf_(wf), vcomm_(wf.sd(0,0)->basis().comm())
ft_[ikp] = new FourierTransform(wb,np0v,np1v,np2v);
}
}
// initialize core density ptr to null ptr
rhocore_r = 0;
}
////////////////////////////////////////////////////////////////////////////////
......@@ -171,8 +173,15 @@ void ChargeDensity::update_density(void)
tmap["charge_vft"].start();
vft_->forward(&rhotmp[0],&rhog[ispin][0]);
tmap["charge_vft"].stop();
// add core correction charge
if ( rhocore_r )
for ( int i = 0; i < rhor_size; i++ )
rhor[ispin][i] += rhocore_r[i];
}
}
////////////////////////////////////////////////////////////////////////////////
void ChargeDensity::update_rhor(void)
{
......@@ -191,10 +200,17 @@ void ChargeDensity::update_rhor(void)
const int rhor_size = rhor[ispin].size();
double *const prhor = &rhor[ispin][0];
#pragma omp parallel for
for ( int i = 0; i < rhor_size; i++ )
if ( rhocore_r )
{
#pragma omp parallel for
for ( int i = 0; i < rhor_size; i++ )
prhor[i] = ( rhotmp[i].real() + rhocore_r[i] ) * omega_inv;
}
else
{
prhor[i] = rhotmp[i].real() * omega_inv;
#pragma omp parallel for
for ( int i = 0; i < rhor_size; i++ )
prhor[i] = rhotmp[i].real() * omega_inv;
}
}
}
......@@ -51,7 +51,8 @@ class ChargeDensity
std::vector<std::vector<double> > rhor; // rhor[ispin][i]
std::vector<std::vector<std::complex<double> > > rhog; // rhog[ispin][ig]
// core density ptr. If non-zero, contains the real-space core density
double* rhocore_r;
void update_density(void);
void update_rhor(void);
......
......@@ -34,11 +34,12 @@
#include <iostream>
#include <iomanip>
#include <algorithm> // fill()
#include <algorithm> // fill(), copy()
#include <functional>
using namespace std;
////////////////////////////////////////////////////////////////////////////////
EnergyFunctional::EnergyFunctional( Sample& s, const ChargeDensity& cd)
EnergyFunctional::EnergyFunctional( Sample& s, ChargeDensity& cd)
: s_(s), cd_(cd)
{
const AtomSet& atoms = s_.atoms;
......@@ -121,6 +122,8 @@ EnergyFunctional::EnergyFunctional( Sample& s, const ChargeDensity& cd)
eself_ = 0.0;
// core_charge_: true if any species has a core charge density
core_charge_ = false;
for ( int is = 0; is < nsp_; is++ )
{
Species *s = atoms.species_list[is];
......@@ -133,8 +136,22 @@ EnergyFunctional::EnergyFunctional( Sample& s, const ChargeDensity& cd)
eself_ += na * s->eself();
na_[is] = na;
zv_[is] = s->zval();
zv_[is] = s->ztot();
rcps_[is] = s->rcps();
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];
}
// FT for interpolation of wavefunctions on the fine grid
......@@ -229,11 +246,32 @@ void EnergyFunctional::update_vhxc(bool compute_stress)
// compute xc energy, update self-energy operator and potential
tmap["exc"].start();
for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
fill(v_r[ispin].begin(),v_r[ispin].end(),0.0);
memset((void*)&v_r[ispin][0], 0, vft->np012loc()*sizeof(double));
xco->update(v_r, compute_stress);
exc_ = xco->exc();
if ( compute_stress )
xco->compute_stress(sigma_exc);
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] );
}
}
}
tmap["exc"].stop();
// compute local potential energy:
......@@ -535,11 +573,22 @@ double EnergyFunctional::energy(bool compute_hpsi, Wavefunction& dwf,
{
const complex<double> sg = s[ig];
const complex<double> rg = rhogt[ig];
Species *s = s_.atoms.species_list[is];
const double g = vbasis_->g_ptr()[ig];
const double gi = vbasis_->gi_ptr()[ig];
// next line: keep only real part
const double temp = fac2 *
( rg.real() * sg.real() +
rg.imag() * sg.imag() )
* rhops[is][ig] * g2i[ig];
// 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;
}
const double tgx = g_x[ig];
const double tgy = g_y[ig];
......@@ -670,6 +719,8 @@ double EnergyFunctional::energy(bool compute_hpsi, Wavefunction& dwf,
{
double tmp = fpi * rhops[is][ig] * g2i[ig];
vtemp[ig] = tmp * conj(rhogt[ig]) + vps[is][ig] * conj(rhoelg[ig]);
if ( core_charge_ )
vtemp[ig] += conj(vxc_g[ig]) * rhocore_sp_g[is][ig];
}
fion_nl.resize(3*na_[is]);
fion_nl = 0.0;
......@@ -863,6 +914,7 @@ double EnergyFunctional::energy(bool compute_hpsi, Wavefunction& dwf,
void EnergyFunctional::atoms_moved(void)
{
const AtomSet& atoms = s_.atoms;
const Wavefunction& wf = s_.wf;
int ngloc = vbasis_->localsize();
// fill tau0 with values in atom_list
......@@ -874,6 +926,29 @@ void EnergyFunctional::atoms_moved(void)
memset( (void*)&vion_local_g[0], 0, 2*ngloc*sizeof(double) );
memset( (void*)&dvion_local_g[0], 0, 2*ngloc*sizeof(double) );
memset( (void*)&rhopst[0], 0, 2*ngloc*sizeof(double) );
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();
}
for ( int is = 0; is < atoms.nsp(); is++ )
{
complex<double> *s = &sf.sfac[is][0];
......@@ -1007,13 +1082,18 @@ void EnergyFunctional::cell_moved(void)
{
Species *s = atoms.species_list[is];
const double * const g = vbasis_->g_ptr();
double v,dv;
double v,dv,rhoc;
for ( int ig = 0; ig < ngloc; ig++ )
{
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;
if ( core_charge_ )
{
s->rhocoreg(g[ig],rhoc);
rhocore_sp_g[is][ig] = rhoc * omega_inv;
}
}
}
......
......@@ -45,7 +45,7 @@ class EnergyFunctional
private:
Sample& s_;
const ChargeDensity& cd_;
ChargeDensity& cd_;
Basis* vbasis_;
FourierTransform *vft;
std::vector<FourierTransform*> ft;
......@@ -54,9 +54,10 @@ class EnergyFunctional
std::vector<std::vector<NonLocalPotential*> > nlp; // nlp[ispin][ikp]
std::vector<ConfinementPotential*> cfp; // cfp[ikp]
std::vector<std::vector<double> > vps, dvps, rhops;
std::vector<std::vector<double> > vps, dvps, rhops, rhocore_sp_g;
std::vector<double> rhocore_r;
std::vector<std::complex<double> > tmp_r, vion_local_g,
dvion_local_g, vlocal_g, rhopst, rhogt, rhoelg, vtemp;
dvion_local_g, vxc_g, vlocal_g, rhopst, rhogt, rhoelg, vtemp, rhocore_g;
std::vector<std::vector<double> > tau0, fion_esr;
std::vector<std::vector<double> > fext;
......@@ -68,6 +69,8 @@ class EnergyFunctional
std::valarray<double> sigma_ekin,sigma_econf,sigma_eps,sigma_ehart,sigma_exc,
sigma_enl, sigma_esr, sigma;
bool core_charge_;
public:
std::vector<std::vector<double> > v_r;
......@@ -99,7 +102,7 @@ class EnergyFunctional
void print(std::ostream& os) const;
EnergyFunctional(Sample& s, const ChargeDensity& cd);
EnergyFunctional(Sample& s, ChargeDensity& cd);
~EnergyFunctional();
};
std::ostream& operator << ( std::ostream& os, const EnergyFunctional& e );
......
......@@ -57,7 +57,7 @@ void NonLocalPotential::init(void)
nsp = atoms_.nsp();
lmax.resize(nsp);
nop.resize(nsp);
lloc.resize(nsp);
lproj.resize(nsp);
na.resize(nsp);
......@@ -83,7 +83,7 @@ void NonLocalPotential::init(void)
{
nspnl++;
na[is] = atoms_.na(is);
lmax[is] = s->lmax();
nop[is] = s->nop();
lloc[is] = s->llocal();
nquad[is] = s->nquad();
......@@ -188,8 +188,9 @@ void NonLocalPotential::init(void)
// compute lproj[is][ipr]
int ipr_base = 0;
for ( int l = 0; l <= lmax[is]; l++ )
for ( int iop = 0; iop < nop[is]; iop++ )
{
int l = s->l(iop);
if ( l != lloc[is] )
{
if ( nquad[is] == 0 )
......@@ -200,7 +201,7 @@ void NonLocalPotential::init(void)
for ( int m = 0; m < 2*l+1; m++ )
{
const int ipr = ipr_base + m;
wt[is][ipr] = s->wsg(l);
wt[is][ipr] = s->wsg(iop);
lproj[is][ipr] = l;
}
ipr_base += 2*l+1;
......@@ -246,9 +247,13 @@ void NonLocalPotential::update_twnl(void)
const double s14pi = sqrt(1.0/fpi);
const double s34pi = sqrt(3.0/fpi);
const double s54pi = sqrt(5.0/fpi);
const double s74pi = sqrt(7.0/fpi);
const double s20pi = sqrt(20.0*pi);
const double s20pi3 = sqrt(20.0*pi/3.0);
const double s3 = sqrt(3.0);
const double s32 = sqrt(1.5);
const double s52 = sqrt(2.5);
const double s15 = 0.5 * s32 * s52;
const double *kpg = basis_.kpg_ptr();
const double *kpgi = basis_.kpgi_ptr();
......@@ -262,8 +267,9 @@ void NonLocalPotential::update_twnl(void)
Species *s = atoms_.species_list[is];
int ilm = 0;
for ( int l = 0; l <= lmax[is]; l++ )
for ( int iop = 0; iop < nop[is]; iop++ )
{
int l = s->l(iop);
if ( l != lloc[is] )
{
if ( l == 0 )
......@@ -273,24 +279,24 @@ void NonLocalPotential::update_twnl(void)
// Kleinman-Bylander
// twnl[is][ipr][ig]
// ipr = ilm = 0
const int ipr = ilm;
// index = ig + ngwl*ipr, i.e. index = ig
double *t0 = &twnl[is][0];
double *t0 = &twnl[is][ngwl*ipr];
// dtwnl[is][ipr][ij][ngwl]
// index = ig + ngwl * ( ij + 6 * ipr ), ipr = 0
// i.e. index = ig + ij * ngwl
double *dt0_xx = &dtwnl[is][0*ngwl];
double *dt0_yy = &dtwnl[is][1*ngwl];
double *dt0_zz = &dtwnl[is][2*ngwl];
double *dt0_xy = &dtwnl[is][3*ngwl];
double *dt0_yz = &dtwnl[is][4*ngwl];
double *dt0_xz = &dtwnl[is][5*ngwl];
double *dt0_xx = &dtwnl[is][ngwl*(0+6*ipr)];
double *dt0_yy = &dtwnl[is][ngwl*(1+6*ipr)];
double *dt0_zz = &dtwnl[is][ngwl*(2+6*ipr)];
double *dt0_xy = &dtwnl[is][ngwl*(3+6*ipr)];
double *dt0_yz = &dtwnl[is][ngwl*(4+6*ipr)];
double *dt0_xz = &dtwnl[is][ngwl*(5+6*ipr)];
// Special case k=G=0 is ok since kpgi[0] = 0.0 at k=G=0
for ( int ig = 0; ig < ngwl; ig++ )
{
double v,dv;
s->dvnlg(0,kpg[ig],v,dv);
s->dvnlg(iop,kpg[ig],v,dv);
t0[ig] = s14pi * v;
......@@ -460,7 +466,7 @@ void NonLocalPotential::update_twnl(void)
{
double v,dv;
const double tg = kpg[ig];
s->dvnlg(l,tg,v,dv);
s->dvnlg(iop,tg,v,dv);
const double tgx = kpg_x[ig];
const double tgy = kpg_y[ig];
......@@ -675,7 +681,7 @@ void NonLocalPotential::update_twnl(void)
double v,dv;
const double tg = kpg[ig];
s->dvnlg(l,tg,v,dv);
s->dvnlg(iop,tg,v,dv);
const double tgx = kpg_x[ig];
const double tgy = kpg_y[ig];
......@@ -952,6 +958,289 @@ void NonLocalPotential::update_twnl(void)
}
ilm += 2*l+1;
}
else if ( l == 3 )
{
// only implemented for Kleiman-Bylander type
assert(nquad[is] == 0);
// Kleinman-Bylander
const int ipr09 = ilm;
const int ipr10 = ilm + 1;
const int ipr11 = ilm + 2;
const int ipr12 = ilm + 3;
const int ipr13 = ilm + 4;
const int ipr14 = ilm + 5;
const int ipr15 = ilm + 6;
double *t09 = &twnl[is][ngwl * ipr09];
double *t10 = &twnl[is][ngwl * ipr10];
double *t11 = &twnl[is][ngwl * ipr11];
double *t12 = &twnl[is][ngwl * ipr12];
double *t13 = &twnl[is][ngwl * ipr13];
double *t14 = &twnl[is][ngwl * ipr14];
double *t15 = &twnl[is][ngwl * ipr15];
// dtwnl[is][ipr][ij][ngwl]
// index = ig + ngwl * ( ij + 6 * ipr )
double *dt09_xx = &dtwnl[is][ngwl * ( 0 + 6 * ipr09 )];
double *dt09_yy = &dtwnl[is][ngwl * ( 1 + 6 * ipr09 )];
double *dt09_zz = &dtwnl[is][ngwl * ( 2 + 6 * ipr09 )];
double *dt09_xy = &dtwnl[is][ngwl * ( 3 + 6 * ipr09 )];
double *dt09_yz = &dtwnl[is][ngwl * ( 4 + 6 * ipr09 )];
double *dt09_xz = &dtwnl[is][ngwl * ( 5 + 6 * ipr09 )];
double *dt10_xx = &dtwnl[is][ngwl * ( 0 + 6 * ipr10 )];
double *dt10_yy = &dtwnl[is][ngwl * ( 1 + 6 * ipr10 )];
double *dt10_zz = &dtwnl[is][ngwl * ( 2 + 6 * ipr10 )];
double *dt10_xy = &dtwnl[is][ngwl * ( 3 + 6 * ipr10 )];
double *dt10_yz = &dtwnl[is][ngwl * ( 4 + 6 * ipr10 )];
double *dt10_xz = &dtwnl[is][ngwl * ( 5 + 6 * ipr10 )];
double *dt11_xx = &dtwnl[is][ngwl * ( 0 + 6 * ipr11 )];
double *dt11_yy = &dtwnl[is][ngwl * ( 1 + 6 * ipr11 )];
double *dt11_zz = &dtwnl[is][ngwl * ( 2 + 6 * ipr11 )];
double *dt11_xy = &dtwnl[is][ngwl * ( 3 + 6 * ipr11 )];
double *dt11_yz = &dtwnl[is][ngwl * ( 4 + 6 * ipr11 )];
double *dt11_xz = &dtwnl[is][ngwl * ( 5 + 6 * ipr11 )];
double *dt12_xx = &dtwnl[is][ngwl * ( 0 + 6 * ipr12 )];
double *dt12_yy = &dtwnl[is][ngwl * ( 1 + 6 * ipr12 )];
double *dt12_zz = &dtwnl[is][ngwl * ( 2 + 6 * ipr12 )];
double *dt12_xy = &dtwnl[is][ngwl * ( 3 + 6 * ipr12 )];
double *dt12_yz = &dtwnl[is][ngwl * ( 4 + 6 * ipr12 )];
double *dt12_xz = &dtwnl[is][ngwl * ( 5 + 6 * ipr12 )];
double *dt13_xx = &dtwnl[is][ngwl * ( 0 + 6 * ipr13 )];
double *dt13_yy = &dtwnl[is][ngwl * ( 1 + 6 * ipr13 )];
double *dt13_zz = &dtwnl[is][ngwl * ( 2 + 6 * ipr13 )];
double *dt13_xy = &dtwnl[is][ngwl * ( 3 + 6 * ipr13 )];
double *dt13_yz = &dtwnl[is][ngwl * ( 4 + 6 * ipr13 )];
double *dt13_xz = &dtwnl[is][ngwl * ( 5 + 6 * ipr13 )];
double *dt14_xx = &dtwnl[is][ngwl * ( 0 + 6 * ipr14 )];
double *dt14_yy = &dtwnl[is][ngwl * ( 1 + 6 * ipr14 )];
double *dt14_zz = &dtwnl[is][ngwl * ( 2 + 6 * ipr14 )];
double *dt14_xy = &dtwnl[is][ngwl * ( 3 + 6 * ipr14 )];
double *dt14_yz = &dtwnl[is][ngwl * ( 4 + 6 * ipr14 )];
double *dt14_xz = &dtwnl[is][ngwl * ( 5 + 6 * ipr14 )];
double *dt15_xx = &dtwnl[is][ngwl * ( 0 + 6 * ipr15 )];
double *dt15_yy = &dtwnl[is][ngwl * ( 1 + 6 * ipr15 )];
double *dt15_zz = &dtwnl[is][ngwl * ( 2 + 6 * ipr15 )];
double *dt15_xy = &dtwnl[is][ngwl * ( 3 + 6 * ipr15 )];
double *dt15_yz = &dtwnl[is][ngwl * ( 4 + 6 * ipr15 )];
double *dt15_xz = &dtwnl[is][ngwl * ( 5 + 6 * ipr15 )];
for ( int ig = 0; ig < ngwl; ig++ )
{
double v, dv;
const double tg = kpg[ig];
s->dvnlg(iop,tg,v,dv);
const double tgx = kpg_x[ig];
const double tgy = kpg_y[ig];
const double tgz = kpg_z[ig];
const double tgx2 = tgx * tgx;
const double tgy2 = tgy * tgy;
const double tgz2 = tgz * tgz;
const double tgx3 = tgx * tgx2;
const double tgy3 = tgy * tgy2;
const double tgz3 = tgz * tgz2;
const double tgi = kpgi[ig];
const double tgi2 = tgi * tgi;
const double tgi3 = tgi * tgi2;
const double tgxx = tgx2 * tgi2;
const double tgyy = tgy2 * tgi2;
const double tgzz = tgz2 * tgi2;
const double tgxy = tgx * tgy * tgi2;
const double tgyz = tgy * tgz * tgi2;
const double tgxz = tgx * tgz * tgi2;
const double tgxxx = tgx3 * tgi3;
const double tgyyy = tgy3 * tgi3;
const double tgzzz = tgz3 * tgi3;
const double tgxyy = tgx * tgy2 * tgi3;
const double tgxzz = tgx * tgz2 * tgi3;
const double tgyxx = tgy * tgx2 * tgi3;
const double tgyzz = tgy * tgz2 * tgi3;
const double tgzxx = tgz * tgx2 * tgi3;
const double tgzyy = tgz * tgy2 * tgi3;
const double tgxyz = tgx * tgy * tgz * tgi3;
const double factor = 0.5 * s74pi;
const double y09 = factor * s52 * ( 3.0 * tgyxx - tgyyy );
const double y10 = factor * s15 * 2.0 * tgxyz;
const double y11 = factor * s32 * ( 4.0 * tgyzz - tgyxx - tgyyy );
const double y12 = factor
* ( 2.0 * tgzzz - 3.0 * ( tgzxx + tgzyy ) );
const double y13 = factor * s32 * ( 4.0 * tgxzz - tgxxx - tgxyy );
const double y14 = factor * s15 * ( tgzxx - tgzyy );
const double y15 = factor * s52 * ( tgxxx - 3.0 * tgxyy );
// derivative of x^3/r^3 w.r.t. x, y, and z
const double dx_xxx = 3.0 * tgxx * ( 1.0 - tgxx ) * tgi;
const double dy_xxx = -3.0 * tgxx * tgxy * tgi;
const double dz_xxx = -3.0 * tgxx * tgxz * tgi;
// derivative of y^3/r^3 w.r.t. x, y, and z
const double dx_yyy = -3.0 * tgyy * tgxy * tgi;
const double dy_yyy = 3.0 * tgyy * ( 1.0 - tgyy ) * tgi;
const double dz_yyy = -3.0 * tgyy * tgyz * tgi;
// derivative of z^3/r^3 w.r.t. x, y, and z
const double dx_zzz = -3.0 * tgzz * tgxz * tgi;
const double dy_zzz = -3.0 * tgzz * tgyz * tgi;
const double dz_zzz = 3.0 * tgzz * ( 1.0 - tgzz ) * tgi;
// derivative of xy^2/r^3 w.r.t. x, y, and z
const double dx_xyy = tgyy * ( 1.0 - 3.0 * tgxx ) * tgi;
const double dy_xyy = tgxy * ( 2.0 - 3.0 * tgyy ) * tgi;
const double dz_xyy = -3.0 * tgyy * tgxz * tgi;
// derivative of xz^2/r^3 w.r.t. x, y, and z
const double dx_xzz = tgzz * ( 1.0 - 3.0 * tgxx ) * tgi;
const double dy_xzz = -3.0 * tgzz * tgxy * tgi;
const double dz_xzz = tgxz * ( 2.0 - 3.0 * tgzz ) * tgi;
// derivative of yx^2/r^3 w.r.t. x, y, and z
const double dx_yxx = tgxy * ( 2.0 - 3.0 * tgxx ) * tgi;
const double dy_yxx = tgxx * ( 1.0 - 3.0 * tgyy ) * tgi;
const double dz_yxx = -3.0 * tgxx * tgyz * tgi;
// derivative of yz^2/r^3 w.r.t. x, y, and z
const double dx_yzz = -3.0 * tgzz * tgxy * tgi;
const double dy_yzz = tgzz * ( 1.0 - 3.0 * tgyy ) * tgi;
const double dz_yzz = tgyz * ( 2.0 - 3.0 * tgzz ) * tgi;
// derivative of zx^2/r^3 w.r.t. x, y, and z
const double dx_zxx = tgxz * ( 2.0 - 3.0 * tgxx ) * tgi;
const double dy_zxx = -3.0 * tgxx * tgyz * tgi;
const double dz_zxx = tgxx * ( 1.0 - 3.0 * tgzz ) * tgi;
// derivative of zy^2/r^3 w.r.t. x, y, and z
const double dx_zyy = -3.0 * tgyy * tgxz * tgi;
const double dy_zyy = tgyz * ( 2.0 - 3.0 * tgyy ) * tgi;
const double dz_zyy = tgyy * ( 1.0 - 3.0 * tgzz ) * tgi;
// derivative of xyz/r^3 w.r.t. x, y, and z
const double dx_xyz = tgyz * ( 1 - 3.0 * tgxx ) * tgi;
const double dy_xyz = tgxz * ( 1 - 3.0 * tgyy ) * tgi;
const double dz_xyz = tgxy * ( 1 - 3.0 * tgzz ) * tgi;
// derivatives of spherical harmonics
// y9 = factor * s52 * ( 3.0 * tgyxx - tgyyy );
const double dx_y09 = factor * s52 * ( 3.0 * dx_yxx - dx_yyy );
const double dy_y09 = factor * s52 * ( 3.0 * dy_yxx - dy_yyy );
const double dz_y09 = factor * s52 * ( 3.0 * dz_yxx - dz_yyy );
// y10 = factor * s15 * 2.0 * tgxyz;
const double dx_y10 = factor * s15 * 2.0 * dx_xyz;
const double dy_y10 = factor * s15 * 2.0 * dy_xyz;
const double dz_y10 = factor * s15 * 2.0 * dz_xyz;
// y11 = factor * s32 * ( 4.0 * tgyzz - tgyxx - tgyyy );
const double dx_y11 = factor * s32 * ( 4.0 * dx_yzz - dx_yxx - dx_yyy );
const double dy_y11 = factor * s32 * ( 4.0 * dy_yzz - dy_yxx - dy_yyy );
const double dz_y11 = factor * s32 * ( 4.0 * dz_yzz - dz_yxx - dz_yyy );
// y12 = factor * ( 2.0 * tgzzz - 3.0 * ( tgzxx + tgzyy ) );
const double dx_y12 = factor
* ( 2.0 * dx_zzz - 3.0 * ( dx_zxx + dx_zyy ) );
const double dy_y12 = factor
* ( 2.0 * dy_zzz - 3.0 * ( dy_zxx + dy_zyy ) );
const double dz_y12 = factor
* ( 2.0 * dz_zzz - 3.0 * ( dz_zxx + dz_zyy ) );
// y13 = factor * s32 * ( 4.0 * tgxzz - tgxxx - tgxyy );
const double dx_y13 = factor * s32 * ( 4.0 * dx_xzz - dx_xxx - dx_xyy );
const double dy_y13 = factor * s32 * ( 4.0 * dy_xzz - dy_xxx - dy_xyy );
const double dz_y13 = factor * s32 * ( 4.0 * dz_xzz - dz_xxx - dz_xyy );
// y14 = factor * s15 * ( tgzxx - tgzyy );
const double dx_y14 = factor * s15 * ( dx_zxx - dx_zyy );
const double dy_y14 = factor * s15 * ( dy_zxx - dy_zyy );
const double dz_y14 = factor * s15 * ( dz_zxx - dz_zyy );
// y15 = factor * s52 * ( tgxxx - 3.0 * tgxyy );
const double dx_y15 = factor * s52 * ( dx_xxx - 3.0 * dz_xyy );
const double dy_y15 = factor * s52 * ( dy_xxx - 3.0 * dy_xyy );
const double dz_y15 = factor * s52 * ( dz_xxx - 3.0 * dz_xyy );
t09[ig] = y09 * v;
t10[ig] = y10 * v;
t11[ig] = y11 * v;
t12[ig] = y12 * v;
t13[ig] = y13 * v;
t14[ig] = y14 * v;
t15[ig] = y15 * v;
// contribution to stress tensor
// sigma_ij = 0.5 * (xi * dj + xj * di)
dt09_xx[ig] = -( v * dx_y09 * tgx - y09 * dv * tg * tgxx );
dt09_yy[ig] = -( v * dy_y09 * tgy - y09 * dv * tg * tgyy );
dt09_zz[ig] = -( v * dz_y09 * tgz - y09 * dv * tg * tgzz );
dt09_xy[ig] = -( 0.5 * v * ( dx_y09 * tgy + dy_y09 * tgx ) //
- y09 * dv * tg * tgxy );
dt09_yz[ig] = -( 0.5 * v * ( dy_y09 * tgz + dz_y09 * tgy ) //
- y09 * dv * tg * tgyz );
dt09_xz[ig] = -( 0.5 * v * ( dx_y09 * tgz + dz_y09 * tgx ) //
- y09 * dv * tg * tgxz );
dt10_xx[ig] = -( v * dx_y10 * tgx - y10 * dv * tg * tgxx );
dt10_yy[ig] = -( v * dy_y10 * tgy - y10 * dv * tg * tgyy );
dt10_zz[ig] = -( v * dz_y10 * tgz - y10 * dv * tg * tgzz );
dt10_xy[ig] = -( 0.5 * v * ( dx_y10 * tgy + dy_y10 * tgx ) //
- y10 * dv * tg * tgxy );
dt10_yz[ig] = -( 0.5 * v * ( dy_y10 * tgz + dz_y10 * tgy ) //
- y10 * dv * tg * tgyz );
dt10_xz[ig] = -( 0.5 * v * ( dx_y10 * tgz + dz_y10 * tgx ) //
- y10 * dv * tg * tgxz );
dt11_xx[ig] = -( v * dx_y11 * tgx - y11 * dv * tg * tgxx );
dt11_yy[ig] = -( v * dy_y11 * tgy - y11 * dv * tg * tgyy );
dt11_zz[ig] = -( v * dz_y11 * tgz - y11 * dv * tg * tgzz );
dt11_xy[ig] = -( 0.5 * v * ( dx_y11 * tgy + dy_y11 * tgx ) //
- y11 * dv * tg * tgxy );
dt11_yz[ig] = -( 0.5 * v * ( dy_y11 * tgz + dz_y11 * tgy ) //
- y11 * dv * tg * tgyz );
dt11_xz[ig] = -( 0.5 * v * ( dx_y11 * tgz + dz_y11 * tgx ) //
- y11 * dv * tg * tgxz );
dt12_xx[ig] = -( v * dx_y12 * tgx - y12 * dv * tg * tgxx );
dt12_yy[ig] = -( v * dy_y12 * tgy - y12 * dv * tg * tgyy );
dt12_zz[ig] = -( v * dz_y12 * tgz - y12 * dv * tg * tgzz );
dt12_xy[ig] = -( 0.5 * v * ( dx_y12 * tgy + dy_y12 * tgx ) //
- y12 * dv * tg * tgxy );
dt12_yz[ig] = -( 0.5 * v * ( dy_y12 * tgz + dz_y12 * tgy ) //
- y12 * dv * tg * tgyz );
dt12_xz[ig] = -( 0.5 * v * ( dx_y12 * tgz + dz_y12 * tgx ) //
- y12 * dv * tg * tgxz );
dt13_xx[ig] = -( v * dx_y13 * tgx - y13 * dv * tg * tgxx );
dt13_yy[ig] = -( v * dy_y13 * tgy - y13 * dv * tg * tgyy );
dt13_zz[ig] = -( v * dz_y13 * tgz - y13 * dv * tg * tgzz );
dt13_xy[ig] = -( 0.5 * v * ( dx_y13 * tgy + dy_y13 * tgx ) //
- y13 * dv * tg * tgxy );
dt13_yz[ig] = -( 0.5 * v * ( dy_y13 * tgz + dz_y13 * tgy ) //
- y13 * dv * tg * tgyz );
dt13_xz[ig] = -( 0.5 * v * ( dx_y13 * tgz + dz_y13 * tgx ) //
- y13 * dv * tg * tgxz );
dt14_xx[ig] = -( v * dx_y14 * tgx - y14 * dv * tg * tgxx );
dt14_yy[ig] = -( v * dy_y14 * tgy - y14 * dv * tg * tgyy );
dt14_zz[ig] = -( v * dz_y14 * tgz - y14 * dv * tg * tgzz );
dt14_xy[ig] = -( 0.5 * v * ( dx_y14 * tgy + dy_y14 * tgx ) //
- y14 * dv * tg * tgxy );
dt14_yz[ig] = -( 0.5 * v * ( dy_y14 * tgz + dz_y14 * tgy ) //
- y14 * dv * tg * tgyz );
dt14_xz[ig] = -( 0.5 * v * ( dx_y14 * tgz + dz_y14 * tgx ) //
- y14 * dv * tg * tgxz );
dt15_xx[ig] = -( v * dx_y15 * tgx - y15 * dv * tg * tgxx );
dt15_yy[ig] = -( v * dy_y15 * tgy - y15 * dv * tg * tgyy );
dt15_zz[ig] = -( v * dz_y15 * tgz - y15 * dv * tg * tgzz );
dt15_xy[ig] = -( 0.5 * v * ( dx_y15 * tgy + dy_y15 * tgx ) //
- y15 * dv * tg * tgxy );
dt15_yz[ig] = -( 0.5 * v * ( dy_y15 * tgz + dz_y15 * tgy ) //
- y15 * dv * tg * tgyz );