Commit 7e5c5e05 by Francois Gygi

Implemented calculation of tau in SlaterDet

parent fbd76df5
......@@ -211,6 +211,31 @@ void ChargeDensity::update_rhor(void)
}
////////////////////////////////////////////////////////////////////////////////
void ChargeDensity::update_taur(void) const
{
if (taur.size() != vft_->np012loc())
{
taur.resize(vft_->np012loc());
}
for ( int ispin = 0; ispin < wf_.nspin(); ispin++ )
{
tmap["update_taur"].start();
fill(taur.begin(),taur.end(),0.0);
for ( int ikp = 0; ikp < wf_.nkp(); ikp++ )
{
assert(taur.size()==ft_[ikp]->np012loc());
wf_.sd(ispin,ikp)->compute_tau(*ft_[ikp],
wf_.weight(ikp), &taur[0]);
}
tmap["update_taur"].stop();
if ( rhocore_r )
assert(false);
}
}
////////////////////////////////////////////////////////////////////////////////
double ChargeDensity::total_charge(void) const
{
assert((wf_.nspin()==1)||(wf_.nspin()==2));
......
......@@ -51,11 +51,13 @@ class ChargeDensity
mutable TimerMap tmap;
std::vector<std::vector<double> > rhor; // rhor[ispin][i]
mutable std::vector<double> taur; // taur[i]
std::vector<std::vector<std::complex<double> > > rhog; // rhog[ispin][ig]
// core density ptr. If non-zero, contains the real-space core density
double* rhocore_r;
void update_density(void);
void update_rhor(void);
void update_taur(void) const;
const Context& context(void) const { return ctxt_; }
MPI_Comm vcomm(void) const { return vcomm_; }
......
......@@ -245,7 +245,7 @@ void SCANFunctional::excSCAN(double rho, double grad, double tau, double *exc,
// SCAN exchange enhancement factor
tmp0 = (b1 * s2 + b2 * oneMalpha * exp(-b3 * oneMalpha * oneMalpha));
x = mu_AK * s2 * (1.0 + (b4 * s2 / mu_AK) * exp(-b4 * s2 / mu_AK)) +
x = mu_AK * s2 * (1.0 + (b4 * s2 / mu_AK) * exp(-b4 * s2 / mu_AK)) +
tmp0 * tmp0;
gx = 1.0 - exp(-a1 / sqrt(s));
hx1 = 1.0 + k1 - k1/(1.0 + x / k1);
......@@ -343,10 +343,10 @@ void SCANFunctional::excSCAN(double rho, double grad, double tau, double *exc,
dbeta1drs = 0.066725 * (0.1 - 0.1778) / tmp2 / tmp2;
dw1drs = - 1.0 * (1.0 + w1) * decLSDAdrs / gamma;
dA1drs = dbeta1drs / (gamma * w1) - beta1 * dw1drs / (gamma * w1 * w1);
dt1drs = -1.0 * pow(3.0 * pi * pi / 16.0, 1.0/3.0) * s / (2.0 * rtrs * rtrs *
dt1drs = -1.0 * pow(3.0 * pi * pi / 16.0, 1.0/3.0) * s / (2.0 * rtrs * rtrs *
rtrs);
dH1drs = (1.0 - g1) * gamma / (1.0 + w1 * (1.0 - g1)) * dw1drs + w1 * gamma /
(1.0 + w1 * (1.0 - g1)) * (t1 * t1 * g1 * g1 * g1 * g1 * g1 *
dH1drs = (1.0 - g1) * gamma / (1.0 + w1 * (1.0 - g1)) * dw1drs + w1 * gamma /
(1.0 + w1 * (1.0 - g1)) * (t1 * t1 * g1 * g1 * g1 * g1 * g1 *
dA1drs + 2.0 * A1 * t1 * g1 * g1 * g1 * g1 * g1 * dt1drs);
dec1drs = decLSDAdrs + dH1drs;
......@@ -378,8 +378,8 @@ void SCANFunctional::excSCAN(double rho, double grad, double tau, double *exc,
dsdn = -4.0 * s / (3.0 * rho);
dalphadn = (tau_W / tau_unif - 5.0 * XCalpha / 3.0) / rho;
decdn = (dec1drs * drsdn + dec1ds * dsdn) +
dfcdalpha * dalphadn * (ec0 - ec1) + fc *
((dec0drs * drsdn + dec0ds * dsdn) -
dfcdalpha * dalphadn * (ec0 - ec1) + fc *
((dec0drs * drsdn + dec0ds * dsdn) -
(dec1drs * drsdn + dec1ds * dsdn));
vc1 = ec + rho * decdn;
......
......@@ -355,6 +355,55 @@ void SlaterDet::compute_density(FourierTransform& ft,
}
////////////////////////////////////////////////////////////////////////////////
void SlaterDet::compute_tau(FourierTransform& ft,
double weight, double* taur) const
{
// compute tau of the states residing on my column of ctxt_
assert(occ_.size() == c_.n());
vector<complex<double> > tmp(ft.np012loc());
const int ngwloc = basis_->localsize();
vector<complex<double> > taug(ngwloc);
const int mloc = c_.mloc();
assert(basis_->cell().volume() > 0.0);
const double omega_inv = 1.0 / basis_->cell().volume();
const int np012loc = ft.np012loc();
if ( basis_->real() )
{
for ( int n = 0; n < nstloc(); n++ )
{
const int nn = ctxt_.mycol() * c_.nb() + n;
// next line: factor 0.5 from definition of tau
const double fac1 = 0.5 * weight * omega_inv * occ_[nn];
const complex<double>* p = c_.cvalptr();
if ( fac1 > 0.0 )
{
for ( int j = 0; j < 3; j++ )
{
const double *const gxj = basis_->gx_ptr(j);
for ( int ig = 0; ig < ngwloc; ig++ )
{
/* i*G_j*c(G) */
taug[ig] = complex<double>(0.0,gxj[ig]) * p[ig+n*mloc];
}
ft.backward(&taug[0],&tmp[0]);
for ( int i = 0; i < np012loc; i++ )
{
taur[i] += fac1 * norm(tmp[i]);
}
}
}
}
}
else
{
assert(0);
}
}
////////////////////////////////////////////////////////////////////////////////
void SlaterDet::rs_mul_add(FourierTransform& ft,
const double* v, SlaterDet& sdp) const
{
......
......@@ -67,6 +67,7 @@ class SlaterDet
void resize(const UnitCell& cell, const UnitCell& refcell,
double ecut, int nst);
void compute_density(FourierTransform& ft, double weight, double* rho) const;
void compute_tau(FourierTransform& ft, double weight, double* taur) const;
void rs_mul_add(FourierTransform& ft, const double* v, SlaterDet& sdp) const;
void randomize(double amplitude);
void cleanup(void);
......
......@@ -237,69 +237,10 @@ void XCPotential::update(vector<vector<double> >& vr)
} // j
}
//!! requires cleanning up
//!! I think the basic structure is right but needs to be checked
//if ( xcf_->isMeta() )
if (0)
if ( xcf_->isMeta() )
{
// compute tau
std::vector<std::complex<double> > tmptau1;
std::vector<double> tmptau2, sum;
tmptau1.resize(np012loc_);
tmptau2.resize(np012loc_);
sum.resize(np012loc_);
double *taur = xcf_->tau;
const Wavefunction& wf = s_.wf;
#pragma omp parallel for
for ( int i = 0; i < np012loc_; i++ )
{
tmptau2[i] = 0.0;
}
for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
{
for ( int ikp = 0; ikp < wf.nkp(); ikp++ )
{
const double weight = wf.weight(ikp);
const SlaterDet& sd = *(wf.sd(ispin,ikp));
const Basis& wfbasis = sd.basis();
// factor fac in next lines: 2.0 for G and -G (if basis is real) and
// 0.5 from 1/(2m)
const double fac = wfbasis.real() ? 1.0 : 0.5;
const ComplexMatrix& c = sd.c();
const Context& sdctxt = sd.context();
const int ngwloc = wfbasis.localsize();
const complex<double>* p = c.cvalptr();
const int pmloc = c.mloc();
const int pnloc = c.nloc();
const int nnbase = sdctxt.mycol() * c.nb();
const double * const occ = sd.occ_ptr(nnbase);
const double omega_inv = 1.0 / vbasis_.cell().volume();
for ( int j = 0; j < 3; j++ )
{
const double *const gxj = vbasis_.gx_ptr(j);
for ( int n = 0; n < pnloc; n++ )
{
for ( int ig = 0; ig < ngwloc; ig++ )
{
/* i*G_j*c(G) */
tmp1[ig] = fac * occ[n] * complex<double>(0.0,omega_inv*gxj[ig])
* p[ig+n*pmloc];
}
vft_.backward(&tmp1[0],&tmptau1[0]);
for ( int i = 0; i < np012loc_; i++ )
{
tmptau2[i] += abs(tmptau1[i]) * abs(tmptau1[i]);
}
}
}
for ( int i = 0; i < np012loc_; i++ )
{
sum[i] += weight * tmptau2[i];
}
}
}
*taur = sum[0];
cd_.update_taur();
}
xcf_->setxc();
......
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