Commit 9c0b5074 authored by Francois Gygi's avatar Francois Gygi
Browse files

Fixed calculation of the quadratic exchange correction contribution to stress.

Added cell_moved() member to update quantities that depend on cell size.


git-svn-id: http://qboxcode.org/svn/qb/trunk@1358 cba15fb0-1239-40c8-b417-11db7ca47a34
parent 1143f241
......@@ -369,6 +369,8 @@ void ExchangeOperator::add_stress(valarray<double>& sigma_exc)
void ExchangeOperator::cell_moved(void)
{
vbasis_->resize( s_.wf.cell(),s_.wf.refcell(),4.0*s_.wf.ecut());
KPGridStat_.cell_moved();
KPGridPerm_.cell_moved();
}
////////////////////////////////////////////////////////////////////////////////
......@@ -903,10 +905,6 @@ double ExchangeOperator::compute_exchange_for_general_case_( Sample* s,
div_corr += div_corr_4;
const double e_div_corr_4 = -0.5 * div_corr_4 * occ_ki_[i];
exchange_sum += e_div_corr_4;
//const double fac4 = ( 4.0 * M_PI / omega );
//sigma_exhf_[0] += ( e_div_corr_4 + fac4 * 2.0 * beta_x ) / omega;
//sigma_exhf_[1] += ( e_div_corr_4 + fac4 * 2.0 * beta_y ) / omega;
//sigma_exhf_[2] += ( e_div_corr_4 + fac4 * 2.0 * beta_z ) / omega;
} // if quad_correction
......@@ -2157,19 +2155,28 @@ double ExchangeOperator::compute_exchange_at_gamma_(const Wavefunction &wf,
double beta_z=(s1_z+s2_z-2.0*s0)/(d1_z*d1_z+d2_z*d2_z)*
KPGridPerm_.integral_kz(0);
#if 0
if (gcontext_.onpe0())
{
cout << "i=" << i << " beta=" << beta_x << " " << beta_y << " " << beta_z
<< " integral=" << KPGridPerm_.integral_kx(0) << " "
<< KPGridPerm_.integral_ky(0) << " "
<< KPGridPerm_.integral_kz(0) << endl;
}
#endif
// note: factor occ_ki_[i] * spinFactor already in beta
const double beta_sum = beta_x + beta_y + beta_z ;
const double div_corr_4 = (4.0 * M_PI / omega ) * beta_sum;
div_corr += div_corr_4;
const double e_div_corr_4 = -0.5 * div_corr_4 * occ_ki_[i];
exchange_sum += e_div_corr_4;
const double fac4 = ( 4.0 * M_PI / omega );
sigma_exhf_[0] += ( e_div_corr_4 + fac4 * 2.0 * beta_x ) / omega;
sigma_exhf_[1] += ( e_div_corr_4 + fac4 * 2.0 * beta_y ) / omega;
sigma_exhf_[2] += ( e_div_corr_4 + fac4 * 2.0 * beta_z ) / omega;
const double fac4 = ( 4.0 * M_PI / omega ) * occ_ki_[i];
sigma_exhf_[0] += ( e_div_corr_4 + fac4 * beta_x ) / omega;
sigma_exhf_[1] += ( e_div_corr_4 + fac4 * beta_y ) / omega;
sigma_exhf_[2] += ( e_div_corr_4 + fac4 * beta_z ) / omega;
}
// contribution of divergence corrections to forces on wave functions
if (dwf)
{
......@@ -2204,8 +2211,9 @@ double ExchangeOperator::compute_exchange_at_gamma_(const Wavefunction &wf,
// scale stress tensor with HF coefficient
sigma_exhf_ *= HFCoeff_;
// print result
tm.stop();
#ifdef DEBUG
if ( gcontext_.onpe0() )
{
cout << setprecision(10);
......@@ -2235,6 +2243,7 @@ double ExchangeOperator::compute_exchange_at_gamma_(const Wavefunction &wf,
<< endl;
}
}
#endif
// return total exchange in Hartree, scaled by HF coefficient
return extot;
......
......@@ -28,22 +28,20 @@ using namespace std;
#define TAG_Overlaps_local 12
////////////////////////////////////////////////////////////////////////////////
KPConnectivity::KPConnectivity(const Sample& s)
KPConnectivity::KPConnectivity(const Sample& s) :s_(s), kdist_tol_(1.e-5)
{
const double kdist_tol = 1.e-5;
int onpe0 = ( ( s.ctxt_.myrow()==0 ) && ( s.ctxt_.mycol()==0 ) );
#ifdef DEBUG
int onpe0 = ( ( s_.ctxt_.myrow()==0 ) && ( s_.ctxt_.mycol()==0 ) );
if ( onpe0 )
cout << "Building kp grid connection for quad HF exchange convergence\n";
#endif
// get number of kpoints
nkpoints_ = s.wf.nkp();
comm_ = s.wf.sd(0,0)->context().comm();
nkpoints_ = s_.wf.nkp();
comm_ = s_.wf.sd(0,0)->context().comm();
// get maximum number of local states
int nStateLoc_= s.wf.sd(0,0)->nstloc();
int nStateLoc_= s_.wf.sd(0,0)->nstloc();
MPI_Allreduce(&nStateLoc_,&nStateMax_,1,MPI_INT,MPI_MAX,comm_);
// allocate memory for kpoint weights
......@@ -120,23 +118,31 @@ KPConnectivity::KPConnectivity(const Sample& s)
for ( int iKp=0; iKp<nkpoints_; iKp++ )
{
// test direct difference
D3vector dk = s.wf.kpoint(iKp) - s.wf.kpoint(0);
if ( abs(dk.x)>kdist_tol ) DimX_=1;
if ( abs(dk.y)>kdist_tol ) DimY_=1;
if ( abs(dk.z)>kdist_tol ) DimZ_=1;
D3vector dk = s_.wf.kpoint(iKp) - s_.wf.kpoint(0);
if ( abs(dk.x)>kdist_tol_ ) DimX_=1;
if ( abs(dk.y)>kdist_tol_ ) DimY_=1;
if ( abs(dk.z)>kdist_tol_ ) DimZ_=1;
// test diference with symmetric
dk = s.wf.kpoint(iKp) + s.wf.kpoint(0);
if ( abs(dk.x)>kdist_tol ) DimX_=1;
if ( abs(dk.y)>kdist_tol ) DimY_=1;
if ( abs(dk.z)>kdist_tol ) DimZ_=1;
dk = s_.wf.kpoint(iKp) + s_.wf.kpoint(0);
if ( abs(dk.x)>kdist_tol_ ) DimX_=1;
if ( abs(dk.y)>kdist_tol_ ) DimY_=1;
if ( abs(dk.z)>kdist_tol_ ) DimZ_=1;
}
// get the length of each vector of the
// reciprocal cell
double length_kx=length(s.wf.cell().b(0));
double length_ky=length(s.wf.cell().b(1));
double length_kz=length(s.wf.cell().b(2));
cell_moved();
}
////////////////////////////////////////////////////////////////////////////////
void KPConnectivity::cell_moved(void)
{
// recompute all quantities that depend on the cell dimensions
int onpe0 = ( ( s_.ctxt_.myrow()==0 ) && ( s_.ctxt_.mycol()==0 ) );
// get the length of each vector of the reciprocal cell
double length_kx=length(s_.wf.cell().b(0));
double length_ky=length(s_.wf.cell().b(1));
double length_kz=length(s_.wf.cell().b(2));
// get the connectivity of the grid
for ( int iKpi=0; iKpi<nkpoints_; iKpi++ )
......@@ -166,14 +172,14 @@ KPConnectivity::KPConnectivity(const Sample& s)
second_distance_kz_[iKpi]=1;
// if this kpoint is part of the integration grid
if ( s.wf.weight(iKpi)!=0.0 )
if ( s_.wf.weight(iKpi)!=0.0 )
{
// explore the remaining k points
ConnectivityComplete_=1;
for ( int iKpj=0; iKpj<nkpoints_; iKpj++ )
{
// if this kpoint is part of the integration grid
if ( s.wf.weight(iKpj)!=0.0 )
if ( s_.wf.weight(iKpj)!=0.0 )
{
for ( int itx = -1; itx <=1; itx++ )
for ( int ity = -1; ity <=1; ity++ )
......@@ -181,14 +187,14 @@ KPConnectivity::KPConnectivity(const Sample& s)
{
D3vector T(itx, ity, itz);
// first test direct difference
D3vector dk = s.wf.kpoint(iKpi) - s.wf.kpoint(iKpj) - T;
D3vector dk = s_.wf.kpoint(iKpi) - s_.wf.kpoint(iKpj) - T;
// make sure that we are not just considering 0 vector
if ( dk.x!=0 || dk.y!=0 || dk.z!=0 )
{
// kx direction
if ( abs(dk.x)<dminx1+kdist_tol && abs(dk.y)<kdist_tol &&
abs(dk.z)<kdist_tol )
if ( abs(dk.x)<dminx1+kdist_tol_ && abs(dk.y)<kdist_tol_ &&
abs(dk.z)<kdist_tol_ )
{
// copy first neighbour in second
dminx2=dminx1;
......@@ -203,8 +209,8 @@ KPConnectivity::KPConnectivity(const Sample& s)
first_T_overlap_kx_[iKpi]=T;
first_distance_kx_[iKpi] =dk.x*length_kx;
}
else if ( abs(dk.x)<dminx2+kdist_tol && abs(dk.y)<kdist_tol &&
abs(dk.z)<kdist_tol )
else if ( abs(dk.x)<dminx2+kdist_tol_ && abs(dk.y)<kdist_tol_ &&
abs(dk.z)<kdist_tol_ )
{
// set new second neighbour
dminx2=abs(dk.x);
......@@ -215,8 +221,8 @@ KPConnectivity::KPConnectivity(const Sample& s)
}
// ky direction
if ( abs(dk.x)<kdist_tol && abs(dk.y)<dminy1+kdist_tol &&
abs(dk.z)<kdist_tol )
if ( abs(dk.x)<kdist_tol_ && abs(dk.y)<dminy1+kdist_tol_ &&
abs(dk.z)<kdist_tol_ )
{
// copy first neighbour in second
dminy2=dminy1;
......@@ -231,8 +237,8 @@ KPConnectivity::KPConnectivity(const Sample& s)
first_T_overlap_ky_[iKpi]=T;
first_distance_ky_[iKpi] =dk.y*length_ky;
}
else if ( abs(dk.x)<kdist_tol && abs(dk.y)<dminy2+kdist_tol &&
abs(dk.z)<kdist_tol )
else if ( abs(dk.x)<kdist_tol_ && abs(dk.y)<dminy2+kdist_tol_ &&
abs(dk.z)<kdist_tol_ )
{
// set new second neighbour
dminy2=abs(dk.y);
......@@ -243,8 +249,8 @@ KPConnectivity::KPConnectivity(const Sample& s)
}
// kz direction
if ( abs(dk.x)<kdist_tol && abs(dk.y)<kdist_tol &&
abs(dk.z)<dminz1+kdist_tol )
if ( abs(dk.x)<kdist_tol_ && abs(dk.y)<kdist_tol_ &&
abs(dk.z)<dminz1+kdist_tol_ )
{
// copy first neighbour in second
dminz2=dminz1;
......@@ -259,8 +265,8 @@ KPConnectivity::KPConnectivity(const Sample& s)
first_T_overlap_kz_[iKpi]=T;
first_distance_kz_[iKpi] =dk.z*length_kz;
}
else if ( abs(dk.x)<kdist_tol && abs(dk.y)<kdist_tol &&
abs(dk.z)<dminz2+kdist_tol )
else if ( abs(dk.x)<kdist_tol_ && abs(dk.y)<kdist_tol_ &&
abs(dk.z)<dminz2+kdist_tol_ )
{
// set new second neighbour
dminz2=abs(dk.z);
......@@ -273,21 +279,21 @@ KPConnectivity::KPConnectivity(const Sample& s)
// then test difference with symmetric
// (except for gamma and 0.5 vectors)
if ( ( s.wf.kpoint(iKpj).x!=0.0 &&
fabs(s.wf.kpoint(iKpj).x)!=0.5 ) ||
( s.wf.kpoint(iKpj).y!=0.0 &&
fabs(s.wf.kpoint(iKpj).y)!=0.5 ) ||
( s.wf.kpoint(iKpj).z!=0.0 &&
fabs(s.wf.kpoint(iKpj).z)!=0.5 ) )
if ( ( s_.wf.kpoint(iKpj).x!=0.0 &&
fabs(s_.wf.kpoint(iKpj).x)!=0.5 ) ||
( s_.wf.kpoint(iKpj).y!=0.0 &&
fabs(s_.wf.kpoint(iKpj).y)!=0.5 ) ||
( s_.wf.kpoint(iKpj).z!=0.0 &&
fabs(s_.wf.kpoint(iKpj).z)!=0.5 ) )
{
dk = s.wf.kpoint(iKpi) + s.wf.kpoint(iKpj) + T;
dk = s_.wf.kpoint(iKpi) + s_.wf.kpoint(iKpj) + T;
// make sure that we are not just considering 0 vector
if ( dk.x!=0 || dk.y!=0 || dk.z!=0 )
{
// kx direction
if ( abs(dk.x)<dminx1+kdist_tol && abs(dk.y)<kdist_tol &&
abs(dk.z)<kdist_tol )
if ( abs(dk.x)<dminx1+kdist_tol_ && abs(dk.y)<kdist_tol_ &&
abs(dk.z)<kdist_tol_ )
{
// copy first neighbour in second
dminx2=dminx1;
......@@ -302,8 +308,8 @@ KPConnectivity::KPConnectivity(const Sample& s)
first_T_overlap_kx_[iKpi]=T;
first_distance_kx_[iKpi] =dk.x*length_kx;
}
else if ( abs(dk.x)<dminx2+kdist_tol && abs(dk.y)<kdist_tol &&
abs(dk.z)<kdist_tol )
else if ( abs(dk.x)<dminx2+kdist_tol_ && abs(dk.y)<kdist_tol_ &&
abs(dk.z)<kdist_tol_ )
{
// set new second neighbour
dminx2=abs(dk.x);
......@@ -314,8 +320,8 @@ KPConnectivity::KPConnectivity(const Sample& s)
}
// ky direction
if ( abs(dk.x)<kdist_tol && abs(dk.y)<dminy1+kdist_tol &&
abs(dk.z)<kdist_tol )
if ( abs(dk.x)<kdist_tol_ && abs(dk.y)<dminy1+kdist_tol_ &&
abs(dk.z)<kdist_tol_ )
{
// copy first neighbour in second
dminy2=dminy1;
......@@ -330,8 +336,8 @@ KPConnectivity::KPConnectivity(const Sample& s)
first_T_overlap_ky_[iKpi]=T;
first_distance_ky_[iKpi] =dk.y*length_ky;
}
else if ( abs(dk.x)<kdist_tol && abs(dk.y)<dminy2+kdist_tol &&
abs(dk.z)<kdist_tol )
else if ( abs(dk.x)<kdist_tol_ && abs(dk.y)<dminy2+kdist_tol_ &&
abs(dk.z)<kdist_tol_ )
{
// set new second neighbour
dminy2=abs(dk.y);
......@@ -342,8 +348,8 @@ KPConnectivity::KPConnectivity(const Sample& s)
}
// kz direction
if ( abs(dk.x)<kdist_tol && abs(dk.y)<kdist_tol &&
abs(dk.z)<dminz1+kdist_tol )
if ( abs(dk.x)<kdist_tol_ && abs(dk.y)<kdist_tol_ &&
abs(dk.z)<dminz1+kdist_tol_ )
{
// copy first neighbour in second
dminz2=dminz1;
......@@ -358,8 +364,8 @@ KPConnectivity::KPConnectivity(const Sample& s)
first_T_overlap_kz_[iKpi]=T;
first_distance_kz_[iKpi] =dk.z*length_kz;
}
else if ( abs(dk.x)<kdist_tol && abs(dk.y)<kdist_tol &&
abs(dk.z)<dminz2+kdist_tol )
else if ( abs(dk.x)<kdist_tol_ && abs(dk.y)<kdist_tol_ &&
abs(dk.z)<dminz2+kdist_tol_ )
{
// set new second neighbour
dminz2=abs(dk.z);
......@@ -402,9 +408,9 @@ KPConnectivity::KPConnectivity(const Sample& s)
// check if the grid is irregular (not implemented yet)
// take into account only the kpoint with non zero weight
if (fabs(second_distance_kx_[iKpi]+first_distance_kx_[iKpi])>kdist_tol ||
fabs(second_distance_ky_[iKpi]+first_distance_ky_[iKpi])>kdist_tol ||
fabs(second_distance_kz_[iKpi]+first_distance_kz_[iKpi])>kdist_tol)
if (fabs(second_distance_kx_[iKpi]+first_distance_kx_[iKpi])>kdist_tol_ ||
fabs(second_distance_ky_[iKpi]+first_distance_ky_[iKpi])>kdist_tol_ ||
fabs(second_distance_kz_[iKpi]+first_distance_kz_[iKpi])>kdist_tol_)
{
if ( onpe0 )
{
......@@ -431,7 +437,7 @@ KPConnectivity::KPConnectivity(const Sample& s)
// compute the weight of this kpoint
// this weight should be non zero only for kpoints
// taken into account during the computation
if ( s.wf.weight(iKpi)!=0.0 )
if ( s_.wf.weight(iKpi)!=0.0 )
{
weight_[iKpi]=fabs((first_distance_kx_[iKpi]-second_distance_kx_[iKpi]) *
(first_distance_ky_[iKpi]-second_distance_ky_[iKpi]) *
......@@ -480,7 +486,7 @@ KPConnectivity::KPConnectivity(const Sample& s)
z[i]=zmin+deltaZ/(n+1)*i;
}
// integrate analytically over tha first coordinate
// integrate analytically over the first coordinate
// then integrate numerically
double IntX=0.0;
double IntY=0.0;
......@@ -506,7 +512,7 @@ KPConnectivity::KPConnectivity(const Sample& s)
IntY/=(n*n);
IntZ/=(n*n);
// the sum of the three integral should be equal to 2
// the sum of the three integrals should be equal to 2
// renormalize
//!! note: this step should not be necessary
//!! replace with an assert to check if the sum is correct
......
......@@ -30,12 +30,14 @@ class KPConnectivity
{
private:
const Sample &s_;
int DimX_;
int DimY_;
int DimZ_;
int nkpoints_;
int nStateMax_;
int ConnectivityComplete_;
double kdist_tol_;
std::vector<int> first_neighbour_kx_;
std::vector<int> first_neighbour_ky_;
......@@ -128,6 +130,7 @@ class KPConnectivity
void AddOverlap(int iKpi, int iKpj, int iLocStatei,
complex<double> *valueDirect, complex<double> *valueSymmetric,
double occupation);
void cell_moved(void);
void StartPermutation(int iKp, int iSendTo, int iRecvFr);
void EndPermutation();
......
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