Commit c4fb722b by Francois Gygi

Merge branch 'develop'

parents f10d1877 b2326e0c
......@@ -253,7 +253,7 @@ double ExchangeOperator::update_energy(bool compute_stress)
if ( gamma_only_ )
return eex_ = compute_exchange_at_gamma_(s_.wf, 0, compute_stress);
else
return eex_ = compute_exchange_for_general_case_(&s_, 0, compute_stress);
return eex_ = compute_exchange_for_general_case_(s_.wf, 0, compute_stress);
}
////////////////////////////////////////////////////////////////////////////////
......@@ -265,7 +265,7 @@ double ExchangeOperator::update_operator(bool compute_stress)
if ( gamma_only_ )
eex_ = compute_exchange_at_gamma_(s_.wf, &dwf0_, compute_stress);
else
eex_ = compute_exchange_for_general_case_(&s_, &dwf0_, compute_stress);
eex_ = compute_exchange_for_general_case_(s_.wf, &dwf0_, compute_stress);
// wf0_ is kept as a reference state
wf0_ = s_.wf;
......@@ -288,12 +288,6 @@ void ExchangeOperator::apply_VXC_(double mix, Wavefunction& wf_ref,
const Context &ctxt = s_.wf.sd(ispin,ikp)->c().context();
if ( s_.wf.sd(ispin,ikp)->basis().real() )
{
#if 0
update_sigma();
DoubleMatrix dc_proxy(dwf.sd(ispin,ikp)->c());
DoubleMatrix dcref_proxy(dwf_ref.sd(ispin,ikp)->c());
dc_proxy += dcref_proxy;
#else
DoubleMatrix c_proxy(s_.wf.sd(ispin,ikp)->c());
DoubleMatrix dc_proxy(dwf.sd(ispin,ikp)->c());
DoubleMatrix cref_proxy(wf_ref.sd(ispin,ikp)->c());
......@@ -324,7 +318,6 @@ void ExchangeOperator::apply_VXC_(double mix, Wavefunction& wf_ref,
// |dwf> += mix * |wf_ref> * matproj2
dc_proxy.gemm('n','n',mix,cref_proxy,matproj2,1.0);
#endif
}
else // complex wave functions
{
......@@ -361,12 +354,12 @@ void ExchangeOperator::apply_VXC_(double mix, Wavefunction& wf_ref,
}
////////////////////////////////////////////////////////////////////////////////
double ExchangeOperator::apply_operator(Wavefunction& dwf)
void ExchangeOperator::apply_operator(Wavefunction& dwf)
{
cout << "ExchangeOperator::apply_operator" << endl;
// apply sigmaHF to s_.wf and store result in dwf
// use the reference function wf0_ and reference sigma(wf) dwf0_
apply_VXC_(1.0, wf0_, dwf0_, dwf);
return eex_;
}
////////////////////////////////////////////////////////////////////////////////
......@@ -387,14 +380,12 @@ void ExchangeOperator::cell_moved(void)
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
double ExchangeOperator::compute_exchange_for_general_case_(const Sample* s,
Wavefunction* dwf, bool compute_stress)
double ExchangeOperator::compute_exchange_for_general_case_
(const Wavefunction& wf, Wavefunction* dwf, bool compute_stress)
{
Timer tm;
tm.start();
const Wavefunction& wf = s->wf;
const double omega = wf.cell().volume();
const int nkpoints = wf.nkp();
const int nspin = wf.nspin();
......
......@@ -43,8 +43,8 @@ class ExchangeOperator
// copy of wf
Wavefunction wfc_;
double compute_exchange_for_general_case_(const Sample* s, Wavefunction* dwf,
bool compute_stress);
double compute_exchange_for_general_case_(const Wavefunction& wf,
Wavefunction* dwf, bool compute_stress);
double compute_exchange_at_gamma_(const Wavefunction &wf, Wavefunction* dwf,
bool compute_stress);
void apply_VXC_(double mix, Wavefunction& wf_ref,
......@@ -190,7 +190,7 @@ class ExchangeOperator
double eex() { return eex_; };
double update_energy(bool compute_stress);
double update_operator(bool compute_stress);
double apply_operator(Wavefunction& dwf);
void apply_operator(Wavefunction& dwf);
void add_stress (std::valarray<double> & sigma_exc);
void cell_moved(void);
};
......
......@@ -433,9 +433,9 @@ void XCPotential::update(vector<vector<double> >& vr)
const double *const v3 = xcf_->vxc3;
const double *const tau = xcf_->tau;
for ( int ir = 0; ir < np012loc_; ir++ )
{
sum += tau[ir] * v3[ir];
}
{
sum += tau[ir] * v3[ir];
}
sum *= vbasis_.cell().volume() / vft_.np012();
double tsum = 0.0;
MPI_Allreduce(&sum,&tsum,1,MPI_DOUBLE,MPI_SUM,vbasis_.comm());
......@@ -450,10 +450,10 @@ void XCPotential::update(vector<vector<double> >& vr)
const double *const tau_up = xcf_->tau_up;
const double *const tau_dn = xcf_->tau_dn;
for ( int ir = 0; ir < np012loc_; ir++ )
{
sum_up += tau_up[ir] * v3_up[ir];
sum_dn += tau_dn[ir] * v3_dn[ir];
}
{
sum_up += tau_up[ir] * v3_up[ir];
sum_dn += tau_dn[ir] * v3_dn[ir];
}
sum_up *= vbasis_.cell().volume() / vft_.np012();
sum_dn *= vbasis_.cell().volume() / vft_.np012();
double tsum_up = 0.0;
......@@ -464,6 +464,7 @@ void XCPotential::update(vector<vector<double> >& vr)
}
}
}
////////////////////////////////////////////////////////////////////////////////
void XCPotential::compute_stress(valarray<double>& sigma_exc)
{
......@@ -1265,102 +1266,99 @@ void XCPotential::apply_meta_operator(Wavefunction& dwf)
{
const Wavefunction& wf0 = s_.wf;
if (wf0.nspin() == 1)
if ( wf0.nspin() == 1 )
{
for ( int ispin = 0; ispin < wf0.nspin(); ispin++ )
for ( int ikp = 0; ikp < wf0.nkp(); ikp++ )
{
for ( int ikp = 0; ikp < wf0.nkp(); ikp++ )
if ( wf0.sd(0,ikp)->basis().real() )
{
if ( wf0.sd(ispin,ikp)->basis().real() )
const int ngwloc = wf0.sd(0,ikp)->basis().localsize();
vector<complex<double> > tmp0(ngwloc);
const int mloc = wf0.sd(0,ikp)->c().mloc();
const complex<double>* p = wf0.sd(0,ikp)->c().cvalptr();
complex<double>* dp = dwf.sd(0,ikp)->c().valptr();
for ( int n = 0; n < wf0.sd(0,ikp)->nstloc()-1; n++, n++ )
{
const int ngwloc = wf0.sd(ispin,ikp)->basis().localsize();
vector<complex<double> > tmp0(ngwloc);
const int mloc = wf0.sd(ispin,ikp)->c().mloc();
const complex<double>* p = wf0.sd(ispin,ikp)->c().cvalptr();
complex<double>* dp = dwf.sd(ispin,ikp)->c().valptr();
for ( int n = 0; n < wf0.sd(ispin,ikp)->nstloc()-1; n++, n++ )
for ( int j = 0; j < 3; j++ )
{
for ( int j = 0; j < 3; j++ )
// Compute Grad_j psi_n(ikp)
const double *const gxj = wf0.sd(0,ikp)->basis().gx_ptr(j);
for ( int ig = 0; ig < ngwloc; ig++ )
{
// Compute Grad_j psi_n(ikp)
const double *const gxj = wf0.sd(ispin,ikp)->basis().gx_ptr(j);
for ( int ig = 0; ig < ngwloc; ig++ )
{
/* i*G_j*c(G) */
tmp0[ig] = complex<double>(0.0,gxj[ig]) * p[ig+n*mloc];
tmp1[ig] = complex<double>(0.0,gxj[ig]) * p[ig+(n+1)*mloc];
}
cd_.ft(ikp)->backward(&tmp0[0],&tmp1[0],&tmpr[0]);
// Compute V3 * Grad_j psi_n(ikp)
for ( int i = 0; i < np012loc_; i++ )
tmpr[i] *= xcf_->vxc3[i];
// Transform to k-space
cd_.ft(ikp)->forward(&tmpr[0],&tmp0[0],&tmp1[0]);
// Compute Div_j[V3 * Grad_j psi_n(ikp)]
// Note -1/2 comes from definition of V3
for ( int ig = 0; ig < ngwloc; ig++ )
{
/* i*G_j*c(G) */
dp[ig+n*mloc] += -0.5 * complex<double>(0.0,gxj[ig]) * tmp0[ig];
dp[ig+(n+1)*mloc] += -0.5 * complex<double>(0.0,gxj[ig]) *
tmp1[ig];
}
/* i*G_j*c(G) */
tmp0[ig] = complex<double>(0.0,gxj[ig]) * p[ig+n*mloc];
tmp1[ig] = complex<double>(0.0,gxj[ig]) * p[ig+(n+1)*mloc];
}
cd_.ft(ikp)->backward(&tmp0[0],&tmp1[0],&tmpr[0]);
// Compute V3 * Grad_j psi_n(ikp)
for ( int i = 0; i < np012loc_; i++ )
tmpr[i] *= xcf_->vxc3[i];
// Transform to k-space
cd_.ft(ikp)->forward(&tmpr[0],&tmp0[0],&tmp1[0]);
// Compute Div_j[V3 * Grad_j psi_n(ikp)]
// Note -1/2 comes from definition of V3
for ( int ig = 0; ig < ngwloc; ig++ )
{
/* i*G_j*c(G) */
dp[ig+n*mloc] += -0.5 * complex<double>(0.0,gxj[ig]) * tmp0[ig];
dp[ig+(n+1)*mloc] += -0.5 * complex<double>(0.0,gxj[ig]) *
tmp1[ig];
}
}
if ( wf0.sd(ispin,ikp)->nstloc() % 2 != 0 )
}
if ( wf0.sd(0,ikp)->nstloc() % 2 != 0 )
{
const int n = wf0.sd(0,ikp)->nstloc()-1;
for ( int j = 0; j < 3; j++ )
{
const int n = wf0.sd(ispin,ikp)->nstloc()-1;
for ( int j = 0; j < 3; j++ )
const double *const gxj = wf0.sd(0,ikp)->basis().gx_ptr(j);
for ( int ig = 0; ig < ngwloc; ig++ )
{
const double *const gxj = wf0.sd(ispin,ikp)->basis().gx_ptr(j);
for ( int ig = 0; ig < ngwloc; ig++ )
{
tmp0[ig] = complex<double>(0.0,gxj[ig]) * p[ig+n*mloc];
}
cd_.ft(ikp)->backward(&tmp0[0],&tmpr[0]);
for ( int i = 0; i < np012loc_; i++ )
tmpr[i] *= xcf_->vxc3[i];
cd_.ft(ikp)->forward(&tmpr[0],&tmp0[0]);
for ( int ig = 0; ig < ngwloc; ig++ )
{
dp[ig+n*mloc] += -0.5 * complex<double>(0.0,gxj[ig]) * tmp0[ig];
}
tmp0[ig] = complex<double>(0.0,gxj[ig]) * p[ig+n*mloc];
}
cd_.ft(ikp)->backward(&tmp0[0],&tmpr[0]);
for ( int i = 0; i < np012loc_; i++ )
tmpr[i] *= xcf_->vxc3[i];
cd_.ft(ikp)->forward(&tmpr[0],&tmp0[0]);
for ( int ig = 0; ig < ngwloc; ig++ )
{
dp[ig+n*mloc] += -0.5 * complex<double>(0.0,gxj[ig]) * tmp0[ig];
}
}
}
else
}
else
{
const int ngwloc = wf0.sd(0,ikp)->basis().localsize();
vector<complex<double> > tmp0(ngwloc);
const int mloc = wf0.sd(0,ikp)->c().mloc();
const complex<double>* p = wf0.sd(0,ikp)->c().cvalptr();
complex<double>* dp = dwf.sd(0,ikp)->c().valptr();
for ( int n = 0; n < wf0.sd(0,ikp)->nstloc(); n++ )
{
const int ngwloc = wf0.sd(ispin,ikp)->basis().localsize();
vector<complex<double> > tmp0(ngwloc);
const int mloc = wf0.sd(ispin,ikp)->c().mloc();
const complex<double>* p = wf0.sd(ispin,ikp)->c().cvalptr();
complex<double>* dp = dwf.sd(ispin,ikp)->c().valptr();
for ( int n = 0; n < wf0.sd(ispin,ikp)->nstloc(); n++ )
for ( int j = 0; j < 3; j++ )
{
for ( int j = 0; j < 3; j++ )
// Compute Grad_j psi_n(ikp)
const double *const kpgxj =
wf0.sd(0,ikp)->basis().kpgx_ptr(j);
for ( int ig = 0; ig < ngwloc; ig++ )
{
// Compute Grad_j psi_n(ikp)
const double *const kpgxj =
wf0.sd(ispin,ikp)->basis().kpgx_ptr(j);
for ( int ig = 0; ig < ngwloc; ig++ )
{
// i*(k+G)_j*c(G)
tmp0[ig] = complex<double>(0.0,kpgxj[ig]) * p[ig+n*mloc];
}
cd_.ft(ikp)->backward(&tmp0[0],&tmpr[0]);
// Compute V3 * Grad_j psi_n(ikp)
for ( int i = 0; i < np012loc_; i++ )
tmpr[i] *= xcf_->vxc3[i];
// Transform to k-space
cd_.ft(ikp)->forward(&tmpr[0],&tmp0[0]);
// Compute Div_j[V3 * Grad_j psi_n(ikp)]
// Note -1/2 comes from definition of V3
for ( int ig = 0; ig < ngwloc; ig++ )
{
// i*(k+G)_j*c(G)
dp[ig+n*mloc] += -0.5 * complex<double>(0.0,kpgxj[ig]) *
tmp0[ig];
}
// i*(k+G)_j*c(G)
tmp0[ig] = complex<double>(0.0,kpgxj[ig]) * p[ig+n*mloc];
}
cd_.ft(ikp)->backward(&tmp0[0],&tmpr[0]);
// Compute V3 * Grad_j psi_n(ikp)
for ( int i = 0; i < np012loc_; i++ )
tmpr[i] *= xcf_->vxc3[i];
// Transform to k-space
cd_.ft(ikp)->forward(&tmpr[0],&tmp0[0]);
// Compute Div_j[V3 * Grad_j psi_n(ikp)]
// Note -1/2 comes from definition of V3
for ( int ig = 0; ig < ngwloc; ig++ )
{
// i*(k+G)_j*c(G)
dp[ig+n*mloc] += -0.5 * complex<double>(0.0,kpgxj[ig]) *
tmp0[ig];
}
}
}
......
......@@ -356,7 +356,7 @@ int main(int argc, char **argv, char **envp)
// first argument is "-server"
string inputfilename(argv[2]);
string outputfilename(argv[3]);
bool echo = false;
bool echo = true;
ui.processCmdsServer(inputfilename, outputfilename, "[qbox]", echo);
}
else
......
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