Commit 0e1b782c by Francois Gygi

kpoints modifs


git-svn-id: http://qboxcode.org/svn/qb/trunk@521 cba15fb0-1239-40c8-b417-11db7ca47a34
parent d75be547
......@@ -3,7 +3,7 @@
// Basis.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: Basis.C,v 1.14 2007-10-19 16:24:04 fgygi Exp $
// $Id: Basis.C,v 1.15 2007-10-19 17:10:58 fgygi Exp $
#include "Basis.h"
#include "Context.h"
......@@ -45,14 +45,16 @@ struct BasisImpl
vector<double> g_; // norm of g vectors g[localsize]
vector<double> kpg_; // norm of g vectors g[localsize]
vector<double> gi_; // inverse norm of g vectors gi[localsize]
vector<double> kpgi_; // inverse norm of k+g vectors kpgi[localsize]
vector<double> g2_; // 2-norm of g vectors g2[localsize]
vector<double> kpg2_; // 2-norm of g vectors g2[localsize]
vector<double> g2i_; // inverse square norm of g vectors g2i[localsize]
vector<double> kpg2i_; // inverse square norm of k+g vectors kpg2i[localsize]
int np_[3]; // cache for the function np
vector<double> gx_; // g vectors components gx[j*localsize+i], j=0,1,2
vector<double> gx2_; // g vectors components^2 gx2[j*localsize+i], j=0,1,2
vector<double> kpgx_; // k+g vectors components kpgx[j*localsize+i], j=0,1,2
vector<int> isort_loc; // index array to access locally sorted vectors
// g2_[isort_loc[i]] < g2_[isort_loc[j]] if i < j
// kpg2_[isort_loc[i]] < kpg2_[isort_loc[j]] if i < j
bool real_; // true if k=0
bool resize(const UnitCell& cell, const UnitCell& refcell, double ecut);
......@@ -119,27 +121,31 @@ int Basis::rod_first(int ipe, int irod) const
{ return pimpl_->rod_first_[ipe][irod]; }
int Basis::idx(int i) const { return pimpl_->idx_[i]; }
double Basis::g(int i) const { return pimpl_->g_[i]; }
double Basis::kpg(int i) const { return pimpl_->kpg_[i]; }
double Basis::gi(int i) const { return pimpl_->gi_[i]; }
double Basis::g2(int i) const { return pimpl_->g2_[i]; }
double Basis::kpg2(int i) const { return pimpl_->kpg2_[i]; }
double Basis::g2i(int i) const { return pimpl_->g2i_[i]; }
double Basis::gx(int i) const { return pimpl_->gx_[i]; }
double Basis::gx2(int i) const { return pimpl_->gx2_[i]; }
double Basis::g(int i) const { return pimpl_->g_[i]; }
double Basis::kpg(int i) const { return pimpl_->kpg_[i]; }
double Basis::gi(int i) const { return pimpl_->gi_[i]; }
double Basis::kpgi(int i) const { return pimpl_->kpgi_[i]; }
double Basis::g2(int i) const { return pimpl_->g2_[i]; }
double Basis::kpg2(int i) const { return pimpl_->kpg2_[i]; }
double Basis::g2i(int i) const { return pimpl_->g2i_[i]; }
double Basis::kpg2i(int i) const { return pimpl_->kpg2i_[i]; }
double Basis::gx(int i) const { return pimpl_->gx_[i]; }
double Basis::kpgx(int i) const { return pimpl_->kpgx_[i]; }
int Basis::isort(int i) const { return pimpl_->isort_loc[i]; }
const int* Basis::idx_ptr(void) const { return &(pimpl_->idx_[0]); }
const double* Basis::g_ptr(void) const { return &(pimpl_->g_[0]); }
const double* Basis::kpg_ptr(void) const { return &(pimpl_->kpg_[0]); }
const double* Basis::gi_ptr(void) const { return &(pimpl_->gi_[0]); }
const double* Basis::g2_ptr(void) const { return &(pimpl_->g2_[0]); }
const double* Basis::kpg2_ptr(void) const { return &(pimpl_->kpg2_[0]); }
const double* Basis::g2i_ptr(void) const { return &(pimpl_->g2i_[0]); }
const double* Basis::g_ptr(void) const { return &(pimpl_->g_[0]); }
const double* Basis::kpg_ptr(void) const { return &(pimpl_->kpg_[0]); }
const double* Basis::gi_ptr(void) const { return &(pimpl_->gi_[0]); }
const double* Basis::kpgi_ptr(void) const { return &(pimpl_->kpgi_[0]); }
const double* Basis::g2_ptr(void) const { return &(pimpl_->g2_[0]); }
const double* Basis::kpg2_ptr(void) const { return &(pimpl_->kpg2_[0]); }
const double* Basis::g2i_ptr(void) const { return &(pimpl_->g2i_[0]); }
const double* Basis::kpg2i_ptr(void) const { return &(pimpl_->kpg2i_[0]); }
const double* Basis::gx_ptr(int j) const
{ return &(pimpl_->gx_[j*pimpl_->localsize_[pimpl_->myrow_]]); }
const double* Basis::gx2_ptr(int j) const
{ return &(pimpl_->gx2_[j*pimpl_->localsize_[pimpl_->myrow_]]); }
const double* Basis::kpgx_ptr(int j) const
{ return &(pimpl_->kpgx_[j*pimpl_->localsize_[pimpl_->myrow_]]); }
////////////////////////////////////////////////////////////////////////////////
inline bool factorizable(int n)
......@@ -345,11 +351,13 @@ bool BasisImpl::resize(const UnitCell& cell, const UnitCell& refcell,
g_.resize(localsize_[myrow_]);
kpg_.resize(localsize_[myrow_]);
gi_.resize(localsize_[myrow_]);
kpgi_.resize(localsize_[myrow_]);
g2_.resize(localsize_[myrow_]);
kpg2_.resize(localsize_[myrow_]);
g2i_.resize(localsize_[myrow_]);
kpg2i_.resize(localsize_[myrow_]);
gx_.resize(3*localsize_[myrow_]);
gx2_.resize(3*localsize_[myrow_]);
kpgx_.resize(3*localsize_[myrow_]);
isort_loc.resize(localsize_[myrow_]);
return true;
}
......@@ -654,7 +662,7 @@ bool BasisImpl::resize(const UnitCell& cell, const UnitCell& refcell,
}
}
// local arrays idx, g, gi, g2i, g2, gx, gx2
// local arrays idx, g, gi, g2i, g2, gx
idx_.resize(3*localsize_[myrow_]);
int i = 0;
for ( int irod = 0; irod < nrod_loc_[myrow_]; irod++ )
......@@ -672,11 +680,13 @@ bool BasisImpl::resize(const UnitCell& cell, const UnitCell& refcell,
g_.resize(localsize_[myrow_]);
kpg_.resize(localsize_[myrow_]);
gi_.resize(localsize_[myrow_]);
kpgi_.resize(localsize_[myrow_]);
g2_.resize(localsize_[myrow_]);
kpg2_.resize(localsize_[myrow_]);
g2i_.resize(localsize_[myrow_]);
kpg2i_.resize(localsize_[myrow_]);
gx_.resize(3*localsize_[myrow_]);
gx2_.resize(3*localsize_[myrow_]);
kpgx_.resize(3*localsize_[myrow_]);
isort_loc.resize(localsize_[myrow_]);
update_g();
......@@ -689,7 +699,7 @@ bool BasisImpl::resize(const UnitCell& cell, const UnitCell& refcell,
////////////////////////////////////////////////////////////////////////////////
void BasisImpl::update_g(void)
{
// compute the values of g, kpg, gi, g2i, g2, kpg2, gx, gx2
// compute the values of g, kpg, gi, g2i, g2, kpg2, gx
// N.B. use the values of cell (not defcell)
const int locsize = localsize_[myrow_];
......@@ -706,6 +716,9 @@ void BasisImpl::update_g(void)
gx_[i] = gt.x;
gx_[locsize+i] = gt.y;
gx_[locsize+locsize+i] = gt.z;
kpgx_[i] = kpgt.x;
kpgx_[locsize+i] = kpgt.y;
kpgx_[locsize+locsize+i] = kpgt.z;
g2_[i] = norm(gt);
g_[i] = sqrt( g2_[i] );
......@@ -714,7 +727,9 @@ void BasisImpl::update_g(void)
kpg_[i] = sqrt( kpg2_[i] );
gi_[i] = g_[i] > 0.0 ? 1.0 / g_[i] : 0.0;
kpgi_[i] = kpg_[i] > 0.0 ? 1.0 / kpg_[i] : 0.0;
g2i_[i] = gi_[i] * gi_[i];
kpg2i_[i] = kpgi_[i] * kpgi_[i];
isort_loc[i] = i;
}
......
......@@ -3,7 +3,7 @@
// Basis.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: Basis.h,v 1.7 2007-10-19 16:24:04 fgygi Exp $
// $Id: Basis.h,v 1.8 2007-10-19 17:10:58 fgygi Exp $
#ifndef BASIS_H
#define BASIS_H
......@@ -63,11 +63,13 @@ class Basis
double g(int i) const; // norm of g vectors g[i]
double kpg(int i) const; // norm of k+g vectors kpg[i]
double gi(int i) const; // inverse norm of g vectors gi[i]
double kpgi(int i) const; // inverse norm of k+g vectors kpgi[i]
double g2(int i) const; // 2-norm of g vectors g2[i]
double kpg2(int i) const; // 2-norm of k+g vectors kpg[i]
double g2i(int i) const; // inverse square norm of g g2i[i]
double kpg2i(int i) const; // inverse square norm of k+g kpg2i[i]
double gx(int i) const; // g vectors gx[i+localsize*j],j=0,1,2
double gx2(int i) const; // g vectors components^2 gx2[i+localsize*j]
double kpgx(int i) const; // k+g vectors kpgx[i+localsize*j],j=0,1,2
int isort(int i) const; // index of vectors locally sorted by norm
......@@ -75,11 +77,13 @@ class Basis
const double* g_ptr(void) const;
const double* kpg_ptr(void) const;
const double* gi_ptr(void) const;
const double* kpgi_ptr(void) const;
const double* g2_ptr(void) const;
const double* kpg2_ptr(void) const;
const double* g2i_ptr(void) const;
const double* kpg2i_ptr(void) const;
const double* gx_ptr(int j) const;
const double* gx2_ptr(int j) const;
const double* kpgx_ptr(int j) const;
double memsize(void) const;
double localmemsize(void) const;
......
......@@ -3,7 +3,7 @@
// ChargeDensity.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: ChargeDensity.C,v 1.12 2007-10-19 16:24:04 fgygi Exp $
// $Id: ChargeDensity.C,v 1.13 2007-10-19 17:10:58 fgygi Exp $
#include "ChargeDensity.h"
#include "Basis.h"
......@@ -44,11 +44,18 @@ wf_(wf), vcontext_(wf.sd(0,0)->basis().context())
{
for ( int ikp = 0; ikp < wf.nkp(); ikp++ )
{
if ( wf.sd(ispin,ikp) != 0 )
ft_[ikp] =
new FourierTransform(wf.sd(ispin,ikp)->basis(),np0v,np1v,np2v);
else
ft_[ikp] = 0;
// check that twice the wf basis fits into np0v,np1v,np2v at each kpoint
const Basis& b = wf.sd(ispin,ikp)->basis();
#if DEBUG
cout << " ChargeDensity: ikp=" << ikp << " basis: " << endl;
cout << " idxmin: " << b.idxmin(0) << "/" << b.idxmin(1)
<< "/" << b.idxmin(2) << endl;
cout << " idxmax: " << b.idxmax(0) << "/" << b.idxmax(1)
<< "/" << b.idxmax(2) << endl;
cout << " ftgrid: " << np0v << "/" << np1v
<< "/" << np2v << endl;
#endif
ft_[ikp] = new FourierTransform(b,np0v,np1v,np2v);
}
}
}
......@@ -85,58 +92,53 @@ void ChargeDensity::update_density(void)
for ( int ispin = 0; ispin < wf_.nspin(); ispin++ )
{
if ( wf_.spincontext(ispin)->active())
assert(rhor[ispin].size() == vft_->np012loc() );
assert(rhotmp.size() == vft_->np012loc() );
const int rhor_size = rhor[ispin].size();
tmap["charge_compute"].start();
double *p = &rhor[ispin][0];
for ( int i = 0; i < rhor_size; i++ )
p[i] = 0.0;
for ( int ikp = 0; ikp < wf_.nkp(); ikp++ )
{
assert(rhor[ispin].size() == vft_->np012loc() );
assert(rhotmp.size() == vft_->np012loc() );
tmap["charge_compute"].start();
for ( int ikp = 0; ikp < wf_.nkp(); ikp++ )
{
if ( wf_.sd(ispin,ikp) != 0 && wf_.sdcontext(ispin,ikp)->active() )
{
wf_.sd(ispin,ikp)->compute_density(*ft_[ikp],
wf_.weight(ikp), &rhor[ispin][0]);
}
}
tmap["charge_compute"].stop();
// sum across rows of spincontext[ispin] context
tmap["charge_rowsum"].start();
wf_.spincontext(ispin)->dsum('r',vft_->np012loc(),1,
&rhor[ispin][0],vft_->np012loc());
tmap["charge_rowsum"].stop();
// check integral of charge density
// compute Fourier coefficients of the charge density
double sum = 0.0;
const int rhor_size = rhor[ispin].size();
const double *const prhor = &rhor[ispin][0];
tmap["charge_integral"].start();
#pragma ivdep
for ( int i = 0; i < rhor_size; i++ )
{
const double prh = prhor[i];
sum += prh;
rhotmp[i] = complex<double>(omega * prh, 0.0);
}
sum *= omega / vft_->np012();
wf_.spincontext(ispin)->dsum('c',1,1,&sum,1);
tmap["charge_integral"].stop();
if ( ctxt_.onpe0() )
{
cout.setf(ios::fixed,ios::floatfield);
cout.setf(ios::right,ios::adjustfield);
cout << " <!-- total_electronic_charge: " << setprecision(8) << sum
<< " -->" << endl;
}
tmap["charge_vft"].start();
vft_->forward(&rhotmp[0],&rhog[ispin][0]);
tmap["charge_vft"].stop();
wf_.sd(ispin,ikp)->compute_density(*ft_[ikp],
wf_.weight(ikp), &rhor[ispin][0]);
}
tmap["charge_compute"].stop();
// sum over kpoints: sum along rows of kpcontext
tmap["charge_rowsum"].start();
wf_.kpcontext()->dsum('r',vft_->np012loc(),1,
&rhor[ispin][0],vft_->np012loc());
tmap["charge_rowsum"].stop();
// check integral of charge density
// compute Fourier coefficients of the charge density
double sum = 0.0;
const double *const prhor = &rhor[ispin][0];
tmap["charge_integral"].start();
for ( int i = 0; i < rhor_size; i++ )
{
const double prh = prhor[i];
sum += prh;
rhotmp[i] = complex<double>(omega * prh, 0.0);
}
sum *= omega / vft_->np012();
// sum on all indices except spin: sum along columns of spincontext
wf_.spincontext()->dsum('c',1,1,&sum,1);
tmap["charge_integral"].stop();
if ( ctxt_.onpe0() )
{
cout.setf(ios::fixed,ios::floatfield);
cout.setf(ios::right,ios::adjustfield);
cout << " <!-- total_electronic_charge: " << setprecision(8) << sum
<< " -->" << endl;
}
tmap["charge_vft"].start();
vft_->forward(&rhotmp[0],&rhog[ispin][0]);
tmap["charge_vft"].stop();
}
}
////////////////////////////////////////////////////////////////////////////////
......@@ -150,20 +152,16 @@ void ChargeDensity::update_rhor(void)
for ( int ispin = 0; ispin < wf_.nspin(); ispin++ )
{
if ( wf_.spincontext(ispin)->active())
assert(rhor[ispin].size() == vft_->np012loc() );
assert(rhotmp.size() == vft_->np012loc() );
vft_->backward(&rhog[ispin][0],&rhotmp[0]);
const int rhor_size = rhor[ispin].size();
double *const prhor = &rhor[ispin][0];
for ( int i = 0; i < rhor_size; i++ )
{
assert(rhor[ispin].size() == vft_->np012loc() );
assert(rhotmp.size() == vft_->np012loc() );
vft_->backward(&rhog[ispin][0],&rhotmp[0]);
const int rhor_size = rhor[ispin].size();
double *const prhor = &rhor[ispin][0];
#pragma ivdep
for ( int i = 0; i < rhor_size; i++ )
{
prhor[i] = rhotmp[i].real() * omega_inv;
}
prhor[i] = rhotmp[i].real() * omega_inv;
}
}
}
......
......@@ -3,7 +3,7 @@
// EnergyFunctional.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: EnergyFunctional.h,v 1.16 2007-10-19 16:24:04 fgygi Exp $
// $Id: EnergyFunctional.h,v 1.17 2007-10-19 17:10:58 fgygi Exp $
#ifndef ENERGYFUNCTIONAL_H
#define ENERGYFUNCTIONAL_H
......@@ -40,7 +40,7 @@ class EnergyFunctional
std::vector<FourierTransform*> ft;
StructureFactor sf;
XCPotential* xcp;
NonLocalPotential* nlp;
std::vector<NonLocalPotential*> nlp; // nlp[ikp]
std::vector<ConfinementPotential*> cfp; // cfp[ikp]
std::vector<std::vector<double> > vps, dvps, rhops;
......
......@@ -3,7 +3,7 @@
// MDWavefunctionStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: MDWavefunctionStepper.C,v 1.8 2007-10-19 16:24:04 fgygi Exp $
// $Id: MDWavefunctionStepper.C,v 1.9 2007-10-19 17:12:02 fgygi Exp $
#include "MDWavefunctionStepper.h"
#include "Wavefunction.h"
......@@ -28,10 +28,6 @@ void MDWavefunctionStepper::update(Wavefunction& dwf)
{
for ( int ikp = 0; ikp < wf_.nkp(); ikp++ )
{
if ( wf_.sd(ispin,ikp) != 0 )
{
if ( wf_.sdcontext(ispin,ikp)->active() )
{
tmap_["md_update_wf"].start();
// Verlet update of wf
// cp = c + (c - cm) - dt2/m * hpsi
......@@ -74,8 +70,6 @@ void MDWavefunctionStepper::update(Wavefunction& dwf)
tmap_["riccati"].start();
wf_.sd(ispin,ikp)->riccati(*(wfv_->sd(ispin,ikp)));
tmap_["riccati"].stop();
}
}
}
}
ekin_em_ = ekin_ep_;
......@@ -96,10 +90,6 @@ void MDWavefunctionStepper::compute_wfm(Wavefunction& dwf)
{
for ( int ikp = 0; ikp < wf_.nkp(); ikp++ )
{
if ( wf_.sd(ispin,ikp) != 0 )
{
if ( wf_.sdcontext(ispin,ikp)->active() )
{
SlaterDet* sd = wf_.sd(ispin,ikp);
double* cptr = (double*) sd->c().valptr();
......@@ -133,8 +123,6 @@ void MDWavefunctionStepper::compute_wfm(Wavefunction& dwf)
tmap_["riccati"].start();
wfv_->sd(ispin,ikp)->riccati(*wf_.sd(ispin,ikp));
tmap_["riccati"].stop();
}
}
}
}
ekin_em_ = 0.0;
......@@ -154,10 +142,6 @@ void MDWavefunctionStepper::compute_wfv(Wavefunction& dwf)
{
for ( int ikp = 0; ikp < wf_.nkp(); ikp++ )
{
if ( wf_.sd(ispin,ikp) != 0 )
{
if ( wf_.sdcontext(ispin,ikp)->active() )
{
// compute final velocity wfv
// v = ( c - cm ) / dt - 0.5 * dt/m * hpsi
......@@ -219,8 +203,6 @@ void MDWavefunctionStepper::compute_wfv(Wavefunction& dwf)
}
}
// Note: *wfv_ now contains the wavefunction velocity
}
}
}
}
}
......@@ -236,10 +218,6 @@ double MDWavefunctionStepper::ekin_eh(void)
{
for ( int ikp = 0; ikp < wf_.nkp(); ikp++ )
{
if ( wf_.sd(ispin,ikp) != 0 )
{
if ( wf_.sdcontext(ispin,ikp)->active() )
{
SlaterDet* sd = wf_.sd(ispin,ikp);
double* cptr = (double*) sd->c().valptr();
double* cptrm = (double*) wfv_->sd(ispin,ikp)->c().valptr();
......@@ -277,8 +255,6 @@ double MDWavefunctionStepper::ekin_eh(void)
// Note: 2 in next line: from (G,-G)
ekin_e += ( 2.0 * occn / dt2bye_ ) * tmpsum;
}
}
}
}
}
wf_.context().dsum(1,1,&ekin_e,1);
......
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