Commit 2ff91c49 by Francois Gygi

implementation of stress (s-only n-l projectors)


git-svn-id: http://qboxcode.org/svn/qb/trunk@173 cba15fb0-1239-40c8-b417-11db7ca47a34
parent fe82c329
......@@ -3,7 +3,7 @@
// BOSampleStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: BOSampleStepper.C,v 1.3 2003-12-01 17:57:44 fgygi Exp $
// $Id: BOSampleStepper.C,v 1.4 2004-02-04 19:55:16 fgygi Exp $
#include "BOSampleStepper.h"
#include "EnergyFunctional.h"
......@@ -34,7 +34,7 @@ void BOSampleStepper::step(EnergyFunctional& e, int niter)
{
AtomSet& atoms = s_.atoms;
Wavefunction& wf = s_.wf;
UnitCell dcell;
valarray<double> sigma(6);
atoms.get_positions(tau0);
atoms.get_velocities(vel);
......@@ -43,11 +43,13 @@ void BOSampleStepper::step(EnergyFunctional& e, int niter)
const string wf_dyn = s_.ctrl.wf_dyn;
const string atoms_dyn = s_.ctrl.atoms_dyn;
const string cell_dyn = s_.ctrl.cell_dyn;
const bool compute_hpsi = ( wf_dyn != "LOCKED" );
const bool compute_forces = ( atoms_dyn != "LOCKED" );
const bool compute_stress = ( s_.ctrl.stress == "ON" );
const bool move_cell = ( cell_dyn != "LOCKED" );
Timer tm_iter;
WavefunctionStepper* wf_stepper = 0;
......@@ -84,19 +86,22 @@ void BOSampleStepper::step(EnergyFunctional& e, int niter)
if ( s_.ctxt_.onpe0() )
cout << " <iteration count=\"" << iter+1 << "\">\n";
if ( compute_forces )
double energy = 0.0;
if ( compute_forces || compute_stress )
{
// compute energy and ionic forces using existing wavefunction
double energy =
e.energy(false,dwf,compute_forces,fion,compute_stress,dcell);
energy =
e.energy(false,dwf,compute_forces,fion,compute_stress,sigma);
if ( s_.ctxt_.onpe0() )
{
cout.setf(ios::fixed,ios::floatfield);
cout.setf(ios::right,ios::adjustfield);
cout << " <ekin> " << setprecision(8)
<< setw(15) << e.ekin() << " </ekin>\n"
<< " <eps> " << setw(15) << e.eps() << " </eps>\n"
<< setw(15) << e.ekin() << " </ekin>\n";
if ( compute_stress )
cout << " <econf> " << setw(15) << e.econf() << " </econf>\n";
cout << " <eps> " << setw(15) << e.eps() << " </eps>\n"
<< " <enl> " << setw(15) << e.enl() << " </enl>\n"
<< " <ecoul> " << setw(15) << e.ecoul() << " </ecoul>\n"
<< " <exc> " << setw(15) << e.exc() << " </exc>\n"
......@@ -105,7 +110,10 @@ void BOSampleStepper::step(EnergyFunctional& e, int niter)
<< " <etotal> " << setw(15) << e.etotal() << " </etotal>\n"
<< flush;
}
}
if ( compute_forces )
{
if ( iter == 0 )
ionic_stepper->preprocess(fion);
......@@ -221,7 +229,7 @@ void BOSampleStepper::step(EnergyFunctional& e, int niter)
for ( int ite = 0; ite < nite_; ite++ )
{
// at the last nite iteration, compute ionic forces for the last
double energy = e.energy(true,dwf,false,fion,false,dcell);
double energy = e.energy(true,dwf,false,fion,false,sigma);
wf_stepper->update(dwf);
......@@ -230,8 +238,13 @@ void BOSampleStepper::step(EnergyFunctional& e, int niter)
cout.setf(ios::fixed,ios::floatfield);
cout.setf(ios::right,ios::adjustfield);
cout << " <ekin_int> " << setprecision(8)
<< setw(15) << e.ekin() << " </ekin_int>\n"
<< " <eps_int> " << setw(15) << e.eps() << " </eps_int>\n"
<< setw(15) << e.ekin() << " </ekin_int>\n";
if ( compute_stress )
{
cout << " <econf_int> " << setw(15) << e.econf()
<< " </econf_int>\n";
}
cout << " <eps_int> " << setw(15) << e.eps() << " </eps_int>\n"
<< " <enl_int> " << setw(15) << e.enl() << " </enl_int>\n"
<< " <ecoul_int> " << setw(15) << e.ecoul() << " </ecoul_int>\n"
<< " <exc_int> " << setw(15) << e.exc() << " </exc_int>\n"
......@@ -293,7 +306,7 @@ void BOSampleStepper::step(EnergyFunctional& e, int niter)
// compute ionic forces at last position for postprocessing of ionic
// positions (Stoermer end step)
double energy =
e.energy(false,dwf,compute_forces,fion,compute_stress,dcell);
e.energy(false,dwf,compute_forces,fion,compute_stress,sigma);
ionic_stepper->postprocess(fion);
}
else
......
......@@ -3,7 +3,7 @@
// SampleStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: CPSampleStepper.C,v 1.1 2003-11-21 20:01:06 fgygi Exp $
// $Id: CPSampleStepper.C,v 1.2 2004-02-04 19:55:16 fgygi Exp $
#include "CPSampleStepper.h"
#include "EnergyFunctional.h"
......@@ -35,7 +35,7 @@ void CPSampleStepper::step(EnergyFunctional& e, int niter)
{
AtomSet& atoms = s_.atoms;
Wavefunction& wf = s_.wf;
UnitCell dcell;
valarray<double> sigma(6);
const double dt = s_.ctrl.dt;
double ekin_ion=0.0,ekin_e, temp_ion, eta;
......@@ -51,7 +51,7 @@ void CPSampleStepper::step(EnergyFunctional& e, int niter)
Timer tm_iter;
double energy =
e.energy(compute_hpsi,dwf,compute_forces,fion,compute_stress,dcell);
e.energy(compute_hpsi,dwf,compute_forces,fion,compute_stress,sigma);
for ( int iter = 0; iter < niter; iter++ )
{
......@@ -133,7 +133,7 @@ void CPSampleStepper::step(EnergyFunctional& e, int niter)
}
energy =
e.energy(compute_hpsi,dwf,compute_forces,fion,compute_stress,dcell);
e.energy(compute_hpsi,dwf,compute_forces,fion,compute_stress,sigma);
if ( s_.ctxt_.mype() == 0 )
cout << " </iteration>" << endl;
......
......@@ -3,25 +3,26 @@
// Control.h:
//
////////////////////////////////////////////////////////////////////////////////
// $Id: Control.h,v 1.5 2003-12-02 20:24:27 fgygi Exp $
// $Id: Control.h,v 1.6 2004-02-04 19:55:17 fgygi Exp $
#ifndef CONTROL_H
#define CONTROL_H
#include <string>
#include <vector>
struct Control
{
// control variables
string wf_dyn, atoms_dyn, cell_dyn; // dynamics string flags
string wf_dyn, atoms_dyn; // dynamics string flags
int nite;
double emass; // electron mass
string fermi; // use Fermi distribution to fill states
double fermi_temp; // temperature of Fermi distribution
double ecutprec;
double ecuts,sigmas,facs; // confinement energy parameters
double prefmbar; // externally applied pressure (Mbar)
string wf_diag;
string lock_abc;
string tcp;
double tcp_rcut;
......@@ -29,10 +30,16 @@ struct Control
double gms_mix; // mixing factor for generalized minimum spread functions
string solvation; // continuum dielectric model for solvent
string thermostat;
double th_temp,th_time; // thermostat control
string stress;
string cell_dyn;
string cell_lock;
double cell_mass;
double ecuts,sigmas,facs; // confinement energy parameters
double ext_stress[6]; // external stress tensor: xx,yy,zz,xy,yz,xz
string xc;
string spin;
int delta_spin;
......@@ -40,11 +47,5 @@ struct Control
double dt;
int iprint;
int timeout;
double th_temp,th_time; // thermostat control
double fermi_temp; // temperature of Fermi distribution
double emass; // electron mass
double vmass; // cell mass
};
#endif
......@@ -3,7 +3,7 @@
// EnergyFunctional.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: EnergyFunctional.C,v 1.15 2003-10-02 17:40:34 fgygi Exp $
// $Id: EnergyFunctional.C,v 1.16 2004-02-04 19:55:16 fgygi Exp $
#include "EnergyFunctional.h"
#include "Sample.h"
......@@ -28,9 +28,15 @@ EnergyFunctional::EnergyFunctional(const Sample& s) : s_(s), cd_(s.wf)
{
const AtomSet& atoms = s_.atoms;
const Wavefunction& wf = s_.wf;
const UnitCell& cell(wf.cell());
double omega = cell.volume();
double ecutv = 4 * wf.ecut();
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);
vbasis_ = cd_.vbasis();
//cout << vbasis_->context().mype() << ": vbasis_->context() = "
......@@ -58,7 +64,7 @@ EnergyFunctional::EnergyFunctional(const Sample& s) : s_(s), cd_(s.wf)
<< vft->np012() << " -->" << endl;
}
int ngloc = vbasis_->localsize();
const int ngloc = vbasis_->localsize();
//cout << " EnergyFunctional: ngloc: " << ngloc << endl;
nsp_ = atoms.nsp();
......@@ -92,7 +98,6 @@ EnergyFunctional::EnergyFunctional(const Sample& s) : s_(s), cd_(s.wf)
{
ft[ikp] = cd_.ft(ikp);
}
}
////////////////////////////////////////////////////////////////////////////////
......@@ -121,10 +126,9 @@ EnergyFunctional::~EnergyFunctional(void)
////////////////////////////////////////////////////////////////////////////////
double EnergyFunctional::energy(bool compute_hpsi, Wavefunction& dwf,
bool compute_forces, vector<vector<double> >& fion,
bool compute_stress, UnitCell& dcell)
bool compute_stress, valarray<double>& sigma)
{
const double fpi = 4.0 * M_PI;
const double *g2i = vbasis_->g2i_ptr();
const Wavefunction& wf = s_.wf;
const UnitCell& cell(wf.cell());
......@@ -134,39 +138,147 @@ double EnergyFunctional::energy(bool compute_hpsi, Wavefunction& dwf,
vector<double> temp(ngloc);
assert(wf.nspin()==1);
if ( compute_forces )
{
for ( int is = 0; is < fion.size(); is++ )
{
for ( int i = 0; i < fion[is].size(); i++ )
{
fion[is][i] = 0.0;
}
}
}
const bool use_confinement = s_.ctrl.ecuts > 0.0;
for ( int is = 0; is < fion.size(); is++ )
for ( int ia = 0; ia < fion[is].size(); ia++ )
fion[is][ia] = 0.0;
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();
}
}
}
// kinetic energy
tmap["ekin"].start();
// 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
ekin_ = 0.0;
econf_ = 0.0;
sigma_ekin = 0.0;
sigma_econf = 0.0;
valarray<double> sum(0.0,14), tsum(0.0,14);
for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
{
for ( int ikp = 0; ikp < wf.nkp(); ikp++ )
{
ekin_ += wf.sd(ispin,ikp)->ekin() * wf.weight(ikp);
}
}
const double weight = wf.weight(ikp);
if ( wf.sd(ispin,ikp) != 0 && wf.sdcontext(ispin,ikp)->active() )
{
const SlaterDet& sd = *(wf.sd(ispin,ikp));
const Basis& wfbasis = wf.sd(ispin,ikp)->basis();
// factor fac in next lines: 2.0 for G and -G (if basis is real) and
// 0.5 from 1/(2m)
// note: if basis is real, the factor of 2.0 for G=0 need not be
// corrected since G^2 = 0
// Note: the calculation of fac in next line is valid only for nkp=1
// If k!=0, kpg2(0) !=0 and the ig=0 coefficient must be dealt with
// separately
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++ )
{
double tmpsum = 0.0;
for ( int n = 0; n < nloc; n++ )
{
const double psi2 = norm(p[ig+n*mloc]);
tmpsum += fac * occ[n] * psi2;
}
psi2sum[ig] = tmpsum;
}
// accumulate contributions to ekin,econf,sigma_ekin,sigma_econf in tsum
tsum = 0.0;
// Note: next lines to be changed to kpg_ptr for nkp>1
const double *const g2 = wfbasis.g2_ptr();
const double *const g_x = wfbasis.gx_ptr(0);
const double *const g_y = wfbasis.gx_ptr(1);
const double *const g_z = wfbasis.gx_ptr(2);
for ( int ig = 0; ig < ngwloc; ig++ )
{
const double psi2s = psi2sum[ig];
// tsum[0]: ekin partial sum
tsum[0] += psi2s * g2[ig];
if ( use_confinement )
{
// tsum[1]: econf partial sum
tsum[1] += psi2s * fstress_[ig];
}
if ( compute_stress )
{
const double tgx = g_x[ig];
const double tgy = g_y[ig];
const double tgz = g_z[ig];
const double fac_ekin = 2.0 * psi2s;
const double fac_econf = psi2s * dfstress_[ig];
tsum[2] += fac_ekin * tgx * tgx;
tsum[3] += fac_ekin * tgy * tgy;
tsum[4] += fac_ekin * tgz * tgz;
tsum[5] += fac_ekin * tgx * tgy;
tsum[6] += fac_ekin * tgy * tgz;
tsum[7] += fac_ekin * tgx * tgz;
tsum[8] += fac_econf * tgx * tgx;
tsum[9] += fac_econf * tgy * tgy;
tsum[10] += fac_econf * tgz * tgz;
tsum[11] += fac_econf * tgx * tgy;
tsum[12] += fac_econf * tgy * tgz;
tsum[13] += fac_econf * tgx * tgz;
}
// tsum contains the contributions to
// ekin,econf,sigma_ekin,sigma_econf from vector ig
} // ig
sum += weight * tsum;
}
} // ikp
} // ispin
// sum contains the contributions to ekin, etc.. from this task
wf.context().dsum(14,1,&sum[0],1);
ekin_ = sum[0];
econf_ = sum[1];
sigma_ekin[0] = sum[2];
sigma_ekin[1] = sum[3];
sigma_ekin[2] = sum[4];
sigma_ekin[3] = sum[5];
sigma_ekin[4] = sum[6];
sigma_ekin[5] = sum[7];
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];
sigma_ekin *= omega_inv;
sigma_econf *= omega_inv;
tmap["ekin"].stop();
tmap["density"].start();
......@@ -190,35 +302,142 @@ double EnergyFunctional::energy(bool compute_hpsi, Wavefunction& dwf,
}
// potential energy: integral of electronic charge times ionic local pot.
tsum = 0.0;
int len=2*ngloc,inc1=1;
eps_ = 2.0 * ddot_(&len,(double*)&rhoelg[0],&inc1,
tsum[0] = 2.0 * ddot_(&len,(double*)&rhoelg[0],&inc1,
(double*)&vion_local_g[0],&inc1);
// remove double counting for G=0
if ( vbasis_->context().myrow() == 0 )
{
eps_ -= real(conj(rhoelg[0])*vion_local_g[0]);
tsum[0] -= real(conj(rhoelg[0])*vion_local_g[0]);
}
tsum[0] *= omega; // tsum[0] contains eps
// Stress from Eps
sigma_eps = 0.0;
if ( compute_stress )
{
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)
const double fac = 2.0 * gi[ig] *
real( conj(rhoelg[ig]) * dvion_local_g[ig] );
const double tgx = g_x[ig];
const double tgy = g_y[ig];
const double tgz = g_z[ig];
tsum[1] += fac * tgx * tgx;
tsum[2] += fac * tgy * tgy;
tsum[3] += fac * tgz * tgz;
tsum[4] += fac * tgx * tgy;
tsum[5] += fac * tgy * tgz;
tsum[6] += fac * tgx * tgz;
}
}
eps_ *= omega;
vbasis_->context().dsum(1,1,&eps_,1);
vbasis_->context().dsum(7,1,&tsum[0],7);
eps_ = tsum[0];
sigma_eps[0] = eps_ * omega_inv + tsum[1];
sigma_eps[1] = eps_ * omega_inv + tsum[2];
sigma_eps[2] = eps_ * omega_inv + tsum[3];
sigma_eps[3] = tsum[4];
sigma_eps[4] = tsum[5];
sigma_eps[5] = tsum[6];
// Hartree energy
tsum = 0.0;
ehart_ = 0.0;
const double *const g2i = vbasis_->g2i_ptr();
double ehsum = 0.0;
for ( int ig = 0; ig < ngloc; ig++ )
{
const complex<double> tmp = rhoelg[ig] + rhopst[ig];
ehart_ += norm(tmp) * g2i[ig];
ehsum += norm(tmp) * g2i[ig];
rhogt[ig] = tmp;
}
// factor 1/2 from definition of Ehart cancels with half sum over G
ehart_ *= omega * 4.0 * M_PI;
vbasis_->context().dsum(1,1,&ehart_,1);
// Note: rhogt[ig] includes a factor 1/Omega
// Factor omega in next line yields prefactor 4 pi / omega in
tsum[0] = omega * fpi * ehsum;
// tsum[0] contains ehart
// Stress from Hartree energy
if ( compute_stress )
{
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++ )
{
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];
tsum[1] += temp * tgx * tgx;
tsum[2] += temp * tgy * tgy;
tsum[3] += temp * tgz * tgz;
tsum[4] += temp * tgx * tgy;
tsum[5] += temp * tgy * tgz;
tsum[6] += temp * tgx * tgz;
}
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];
// next line: keep only real part
const double temp = fac2 *
( rg.real() * sg.real() +
rg.imag() * sg.imag() )
* rhops[is][ig] * g2i[ig];
const double tgx = g_x[ig];
const double tgy = g_y[ig];
const double tgz = g_z[ig];
tsum[1] += temp * tgx * tgx;
tsum[2] += temp * tgy * tgy;
tsum[3] += temp * tgz * tgz;
tsum[4] += temp * tgx * tgy;
tsum[5] += temp * tgy * tgz;
tsum[6] += temp * tgx * tgz;
}
}
} // compute_stress
vbasis_->context().dsum(7,1,&tsum[0],7);
ehart_ = tsum[0];
// 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[1];
sigma_ehart[1] = ehart_ * omega_inv - 2.0 * fpi * tsum[2];
sigma_ehart[2] = ehart_ * omega_inv - 2.0 * fpi * tsum[3];
sigma_ehart[3] = - 2.0 * fpi * tsum[4];
sigma_ehart[4] = - 2.0 * fpi * tsum[5];
sigma_ehart[5] = - 2.0 * fpi * tsum[6];
tmap["exc"].start();
// XC energy and potential
for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
fill(v_r[ispin].begin(),v_r[ispin].end(),0.0);
// boolean in next line: compute_stress
xcp->update(v_r,false);
xcp->update(v_r,compute_stress,sigma_exc);
exc_ = xcp->exc();
tmap["exc"].stop();
......@@ -226,11 +445,11 @@ double EnergyFunctional::energy(bool compute_hpsi, Wavefunction& dwf,
// Note: next line for nspin==0, nkp==0 only
tmap["nonlocal"].start();
enl_ = nlp->energy(compute_hpsi,*dwf.sd(0,0),
compute_forces, fion, compute_stress, dcell);
compute_forces, fion, compute_stress, sigma_enl);
tmap["nonlocal"].stop();
ecoul_ = ehart_ + esr_ - eself_;
etotal_ = ekin_ + eps_ + enl_ + ecoul_ + exc_;
etotal_ = ekin_ + econf_ + eps_ + enl_ + ecoul_ + exc_;
if ( compute_hpsi )
{
......@@ -279,13 +498,26 @@ double EnergyFunctional::energy(bool compute_hpsi, Wavefunction& dwf,
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();
for ( int n = 0; n < sd.nstloc(); n++ )
{
// Laplacian
for ( int ig = 0; ig < wfbasis.localsize(); ig++ )
if ( use_confinement )
{
for ( int ig = 0; ig < ngwloc; ig++ )
{
const double fac = 0.5 * ( 1.0 + fstress_[ig] );
cp[ig+mloc*n] += fac * kpg2[ig] * c[ig+mloc*n];
}
}
else
{
cp[ig+mloc*n] += 0.5 * wfbasis.kpg2(ig) * c[ig+mloc*n];
for ( int ig = 0; ig < ngwloc; ig++ )
{
cp[ig+mloc*n] += 0.5 * kpg2[ig] * c[ig+mloc*n];
}
}
}
sd.rs_mul_add(*ft[ikp], &v_r[ispin][0], sdp);
......@@ -373,6 +605,74 @@ double EnergyFunctional::energy(bool compute_hpsi, Wavefunction& dwf,
}
}
sigma = sigma_ekin + sigma_econf + sigma_eps + sigma_enl +
sigma_ehart + sigma_exc + sigma_esr;
if ( compute_stress && s_.ctxt_.onpe0() )
{
const double gpa = 29421.5;
cout.setf(ios::fixed,ios::floatfield);
cout.setf(ios::right,ios::adjustfield);
cout << setprecision(8);
cout << " <stress_tensor>\n"
<< " <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"
<< endl
<< " <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"
<< endl
<< " <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"
<< endl
<< " <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"
<< endl
<< " <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"
<< endl
<< " <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"
<< endl
<< " <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"
<< endl
<< " <sigma_xx> " << setw(12) << sigma[0] << " </sigma_xx>\n"
<< " <sigma_yy> " << setw(12) << sigma[1] << " </sigma_yy>\n"
<< " <sigma_zz> " << setw(12) << sigma[2] << " </sigma_zz>\n"
<< " <sigma_xy> " << setw(12) << sigma[3] << " </sigma_xy>\n"
<< " <sigma_yz> " << setw(12) << sigma[4] << " </sigma_yz>\n"
<< " <sigma_xz> " << setw(12) << sigma[5] << " </sigma_xz>\n"
<< " </stress_tensor>" << endl;
}
return etotal_;
}
......@@ -381,6 +681,7 @@ void EnergyFunctional::init(void)
{
int ngloc = vbasis_->localsize();
vion_local_g.resize(ngloc);
dvion_local_g.resize(ngloc);
vlocal_g.resize(ngloc);
vtemp.resize(ngloc);
rhoelg.resize(ngloc);
......@@ -400,9 +701,14 @@ void EnergyFunctional::init(void)
tau0.resize(nsp_);
taum.resize(nsp_);