Commit 8584567d by Francois Gygi

charge mixing of both spins together

git-svn-id: http://qboxcode.org/svn/qb/trunk@1290 cba15fb0-1239-40c8-b417-11db7ca47a34
parent 102f1a88
......@@ -285,26 +285,17 @@ void BOSampleStepper::step(int niter)
mlwft[ispin] = new MLWFTransform(*wf.sd(ispin,0));
}
// Charge mixing variables
vector<vector<complex<double> > > rhog_in(nspin);
vector<vector<complex<double> > > drhog(nspin);
vector<vector<complex<double> > > rhobar(nspin);
vector<vector<complex<double> > > drhobar(nspin);
for ( int ispin = 0; ispin < nspin; ispin++ )
{
rhog_in[ispin].resize(cd_.rhog[ispin].size());
drhog[ispin].resize(rhog_in[ispin].size());
rhobar[ispin].resize(rhog_in[ispin].size());
drhobar[ispin].resize(rhog_in[ispin].size());
}
// Charge mixing variables: include both spins in the same vector
if ( nspin > 1 ) assert(cd_.rhog[0].size()==cd_.rhog[1].size());
const int ng = cd_.rhog[0].size();
vector<complex<double> > rhog_in(nspin*ng), drhog(nspin*ng),
rhobar(nspin*ng), drhobar(nspin*ng);
const int anderson_ndim = s_.ctrl.charge_mix_ndim;
vector<double> wkerker(cd_.rhog[0].size());
vector<double> wls(cd_.rhog[0].size());
vector<double> wkerker(ng), wls(ng);
vector<AndersonMixer*> mixer(nspin);
for ( int ispin = 0; ispin < nspin; ispin++ )
mixer[ispin] =
new AndersonMixer(2*rhog_in[ispin].size(),anderson_ndim,&cd_.vcontext());
// Anderson charge mixer: include both spins in the same vector
// Factor of 2: complex coeffs stored as double
AndersonMixer mixer(2*nspin*ng,anderson_ndim,&cd_.vcontext());
// compute Kerker preconditioning
// real space Kerker cutoff in a.u.
......@@ -741,10 +732,7 @@ void BOSampleStepper::step(int niter)
{
wf_stepper->preprocess();
if ( anderson_charge_mixing )
{
for ( int ispin = 0; ispin < nspin; ispin++ )
mixer[ispin]->restart();
}
mixer.restart();
double ehart, ehart_m;
......@@ -770,7 +758,8 @@ void BOSampleStepper::step(int niter)
// memorize rhog in rhog_in for next iteration
for ( int ispin = 0; ispin < nspin; ispin++ )
rhog_in[ispin] = cd_.rhog[ispin];
for ( int i = 0; i < ng; i++ )
rhog_in[i+ng*ispin] = cd_.rhog[ispin][i];
}
else
{
......@@ -778,9 +767,9 @@ void BOSampleStepper::step(int niter)
// compute unscreened correction drhog
for ( int ispin = 0; ispin < nspin; ispin++ )
{
for ( int i = 0; i < rhog_in[ispin].size(); i++ )
for ( int i = 0; i < ng; i++ )
{
drhog[ispin][i] = (cd_.rhog[ispin][i] - rhog_in[ispin][i]);
drhog[i+ng*ispin] = (cd_.rhog[ispin][i]-rhog_in[i+ng*ispin]);
}
}
......@@ -791,58 +780,49 @@ void BOSampleStepper::step(int niter)
// row weighting of LS calculation
for ( int ispin = 0; ispin < nspin; ispin++ )
{
for ( int i = 0; i < drhog[ispin].size(); i++ )
drhog[ispin][i] /= wls[i];
for ( int i = 0; i < ng; i++ )
drhog[i+ng*ispin] /= wls[i];
}
const Context * const kpctxt = s_.wf.kpcontext();
if ( kpctxt->mycol() == 0 )
{
// use AndersonMixer on first column only and bcast results
for ( int ispin = 0; ispin < nspin; ispin++ )
{
mixer[ispin]->update((double*)&rhog_in[ispin][0],
(double*)&drhog[ispin][0],
(double*)&rhobar[ispin][0],
(double*)&drhobar[ispin][0]);
const int n = 2*rhobar[ispin].size();
kpctxt->dbcast_send('r',n,1,(double*)&rhobar[ispin][0],n);
kpctxt->dbcast_send('r',n,1,(double*)&drhobar[ispin][0],n);
}
mixer.update((double*)&rhog_in[0], (double*)&drhog[0],
(double*)&rhobar[0], (double*)&drhobar[0]);
const int n = 2*nspin*ng;
kpctxt->dbcast_send('r',n,1,(double*)&rhobar[0],n);
kpctxt->dbcast_send('r',n,1,(double*)&drhobar[0],n);
}
else
{
for ( int ispin = 0; ispin < nspin; ispin++ )
{
const int n = 2*rhobar[ispin].size();
kpctxt->dbcast_recv('r',n,1,
(double*)&rhobar[ispin][0],n,-1,0);
kpctxt->dbcast_recv('r',n,1,
(double*)&drhobar[ispin][0],n,-1,0);
}
const int n = 2*nspin*ng;
kpctxt->dbcast_recv('r',n,1,(double*)&rhobar[0],n,-1,0);
kpctxt->dbcast_recv('r',n,1,(double*)&drhobar[0],n,-1,0);
}
for ( int ispin = 0; ispin < nspin; ispin++ )
{
for ( int i = 0; i < drhog[ispin].size(); i++ )
drhobar[ispin][i] *= wls[i];
for ( int i = 0; i < rhog_in[ispin].size(); i++ )
rhog_in[ispin][i] = rhobar[ispin][i] +
alpha * drhobar[ispin][i] * wkerker[i];
for ( int i = 0; i < ng; i++ )
drhobar[i+ng*ispin] *= wls[i];
for ( int i = 0; i < ng; i++ )
rhog_in[i+ng*ispin] = rhobar[i+ng*ispin] +
alpha * drhobar[i+ng*ispin] * wkerker[i];
}
}
else
{
for ( int ispin = 0; ispin < nspin; ispin++ )
{
for ( int i = 0; i < rhog_in[ispin].size(); i++ )
rhog_in[ispin][i] += alpha * drhog[ispin][i] * wkerker[i];
for ( int i = 0; i < ng; i++ )
rhog_in[i+ng*ispin] += alpha * drhog[i+ng*ispin] * wkerker[i];
}
}
for ( int ispin = 0; ispin < nspin; ispin++ )
{
cd_.rhog[ispin] = rhog_in[ispin];
for ( int i = 0; i < ng; i++ )
cd_.rhog[ispin][i] = rhog_in[i+ng*ispin];
}
cd_.update_rhor();
}
......@@ -1234,11 +1214,9 @@ void BOSampleStepper::step(int niter)
s_.wfv = 0;
}
for ( int ispin = 0; ispin < nspin; ispin++ )
{
delete mlwft[ispin];
delete mixer[ispin];
}
// delete steppers
......
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