Commit 204d57f4 authored by Francois Gygi's avatar Francois Gygi
Browse files

update to use Basis without Context

git-svn-id: http://qboxcode.org/svn/qb/trunk@1340 cba15fb0-1239-40c8-b417-11db7ca47a34
parent c653e15b
......@@ -87,14 +87,12 @@ EnergyFunctional::EnergyFunctional( Sample& s, const ChargeDensity& cd)
zv_.resize(nsp_);
rcps_.resize(nsp_);
na_.resize(nsp_);
namax_ = 0;
for ( int is = 0; is < nsp_; is++ )
{
vps[is].resize(ngloc);
dvps[is].resize(ngloc);
rhops[is].resize(ngloc);
if ( atoms.na(is) > namax_ ) namax_ = atoms.na(is);
}
xco = new XCOperator(s_, cd_);
......@@ -120,7 +118,6 @@ EnergyFunctional::EnergyFunctional( Sample& s, const ChargeDensity& cd)
tau0.resize(nsp_);
fion_esr.resize(nsp_);
fext.resize(nsp_);
ftmp.resize(3*namax_);
eself_ = 0.0;
......@@ -210,7 +207,7 @@ void EnergyFunctional::update_vhxc(bool compute_stress)
const double *const g2i = vbasis_->g2i_ptr();
const double fpi = 4.0 * M_PI;
const int ngloc = vbasis_->localsize();
double tsum[2];
double sum[2], tsum[2];
// compute total electronic density: rhoelg = rho_up + rho_dn
if ( wf.nspin() == 1 )
......@@ -242,14 +239,14 @@ void EnergyFunctional::update_vhxc(bool compute_stress)
// compute local potential energy:
// integral of el. charge times ionic local pot.
int len=2*ngloc,inc1=1;
tsum[0] = 2.0 * ddot(&len,(double*)&rhoelg[0],&inc1,
sum[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 )
if ( vbasis_->mype() == 0 )
{
tsum[0] -= real(conj(rhoelg[0])*vion_local_g[0]);
sum[0] -= real(conj(rhoelg[0])*vion_local_g[0]);
}
tsum[0] *= omega; // tsum[0] contains eps
sum[0] *= omega; // sum[0] contains eps
// Hartree energy
ehart_ = 0.0;
......@@ -263,10 +260,10 @@ void EnergyFunctional::update_vhxc(bool compute_stress)
// 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
tsum[1] = omega * fpi * ehsum;
sum[1] = omega * fpi * ehsum;
// tsum[1] contains ehart
vbasis_->context().dsum(2,1,&tsum[0],2);
MPI_Allreduce(sum,tsum,2,MPI_DOUBLE,MPI_SUM,vbasis_->comm());
eps_ = tsum[0];
ehart_ = tsum[1];
......@@ -472,7 +469,7 @@ double EnergyFunctional::energy(bool compute_hpsi, Wavefunction& dwf,
sigma_eps = 0.0;
if ( compute_stress )
{
tsum = 0.0;
sum = 0.0;
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);
......@@ -488,14 +485,14 @@ double EnergyFunctional::energy(bool compute_hpsi, Wavefunction& dwf,
const double tgy = g_y[ig];
const double tgz = g_z[ig];
tsum[0] += fac * tgx * tgx;
tsum[1] += fac * tgy * tgy;
tsum[2] += fac * tgz * tgz;
tsum[3] += fac * tgx * tgy;
tsum[4] += fac * tgy * tgz;
tsum[5] += fac * tgx * tgz;
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;
}
vbasis_->context().dsum(6,1,&tsum[0],6);
MPI_Allreduce(&sum[0],&tsum[0],6,MPI_DOUBLE,MPI_SUM,vbasis_->comm());
sigma_eps[0] = eps_ * omega_inv + tsum[0];
sigma_eps[1] = eps_ * omega_inv + tsum[1];
......@@ -508,7 +505,7 @@ double EnergyFunctional::energy(bool compute_hpsi, Wavefunction& dwf,
// Stress from Hartree energy
if ( compute_stress )
{
tsum = 0.0;
sum = 0.0;
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);
......@@ -520,12 +517,12 @@ double EnergyFunctional::energy(bool compute_hpsi, Wavefunction& dwf,
const double tgy = g_y[ig];
const double tgz = g_z[ig];
tsum[0] += temp * tgx * tgx;
tsum[1] += temp * tgy * tgy;
tsum[2] += temp * tgz * tgz;
tsum[3] += temp * tgx * tgy;
tsum[4] += temp * tgy * tgz;
tsum[5] += temp * tgx * tgz;
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;
}
for ( int is = 0; is < nsp_; is++ )
......@@ -548,15 +545,15 @@ double EnergyFunctional::energy(bool compute_hpsi, Wavefunction& dwf,
const double tgy = g_y[ig];
const double tgz = g_z[ig];
tsum[0] += temp * tgx * tgx;
tsum[1] += temp * tgy * tgy;
tsum[2] += temp * tgz * tgz;
tsum[3] += temp * tgx * tgy;
tsum[4] += temp * tgy * tgz;
tsum[5] += temp * tgx * tgz;
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;
}
}
vbasis_->context().dsum(6,1,&tsum[0],6);
MPI_Allreduce(&sum[0],&tsum[0],6,MPI_DOUBLE,MPI_SUM,vbasis_->comm());
// Factor in next line:
// factor 2 from G and -G
// factor fpi from definition of sigma
......@@ -665,6 +662,7 @@ double EnergyFunctional::energy(bool compute_hpsi, Wavefunction& dwf,
const double* gx0 = vbasis_->gx_ptr(0);
const double* gx1 = vbasis_->gx_ptr(1);
const double* gx2 = vbasis_->gx_ptr(2);
valarray<double> fion_nl, fion_nl_tmp;
for ( int is = 0; is < nsp_; is++ )
{
......@@ -673,7 +671,9 @@ 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]);
}
memset((void*)&ftmp[0],0,3*namax_*sizeof(double));
fion_nl.resize(3*na_[is]);
fion_nl = 0.0;
fion_nl_tmp.resize(3*na_[is]);
// loop over atoms of species is
for ( int ia = 0; ia < na_[is]; ia++ )
{
......@@ -720,19 +720,19 @@ double EnergyFunctional::energy(bool compute_hpsi, Wavefunction& dwf,
sum1 += tmp * gx1[ig];
sum2 += tmp * gx2[ig];
}
ftmp[3*ia] = sum0;
ftmp[3*ia+1] = sum1;
ftmp[3*ia+2] = sum2;
fion_nl[3*ia] = sum0;
fion_nl[3*ia+1] = sum1;
fion_nl[3*ia+2] = sum2;
}
int len = 3*na_[is];
vbasis_->context().dsum(len,1,&ftmp[0],len);
MPI_Allreduce(&fion_nl[0],&fion_nl_tmp[0],3*na_[is],
MPI_DOUBLE,MPI_SUM,vbasis_->comm());
double fac = -2.0 * omega;
for ( int i = 0; i < 3*na_[is]; i++ )
{
fion[is][i] += fion_esr[is][i] + fac * ftmp[i];
fion[is][i] += fion_esr[is][i] + fac * fion_nl_tmp[i];
}
// add external forces
......
......@@ -57,13 +57,11 @@ class EnergyFunctional
std::vector<std::vector<double> > vps, dvps, rhops;
std::vector<std::complex<double> > tmp_r, vion_local_g,
dvion_local_g, vlocal_g, rhopst, rhogt, rhoelg, vtemp;
std::vector<double> ftmp;
std::vector<std::vector<double> > tau0, fion_esr;
std::vector<std::vector<double> > fext;
std::vector<double> zv_, rcps_;
std::vector<int> na_;
int namax_;
int nsp_;
double ekin_, econf_, eps_, enl_, ehart_,
ecoul_, exc_, esr_, eself_, ets_, eexf_, etotal_;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment