Skip to content

qbox-public
qbox-public
Commit 16b5cf07 by Francois Gygi
Options

kpoints modifs 


git-svn-id: http://qboxcode.org/svn/qb/trunk@522 cba15fb0-1239-40c8-b417-11db7ca47a34
parent 0e1b782c
Showing with 166 additions and 123 deletions
src/PSDWavefunctionStepper.C
......@@ -3,7 +3,7 @@
// PSDWavefunctionStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: PSDWavefunctionStepper.C,v 1.8 2007-10-19 16:24:04 fgygi Exp $
// $Id: PSDWavefunctionStepper.C,v 1.9 2007-10-19 17:37:06 fgygi Exp $
#include "PSDWavefunctionStepper.h"
#include "Wavefunction.h"
......@@ -25,69 +25,67 @@ void PSDWavefunctionStepper::update(Wavefunction& dwf)
{
for ( int ikp = 0; ikp < wf_.nkp(); ikp++ )
{
if ( wf_.sd(ispin,ikp) != 0 )
{
if ( wf_.sdcontext(ispin,ikp)->active() )
{
// compute A = V^T H V and descent direction HV - VA
tmap_["psd_residual"].start();
if ( wf_.sd(ispin,ikp)->basis().real() )
{
// proxy real matrices c, cp
DoubleMatrix c(wf_.sd(ispin,ikp)->c());
DoubleMatrix cp(dwf.sd(ispin,ikp)->c());
// compute A = V^T H V and descent direction HV - VA
DoubleMatrix a(c.context(),c.n(),c.n(),c.nb(),c.nb());
tmap_["psd_residual"].start();
if ( wf_.sd(ispin,ikp)->basis().real() )
{
// proxy real matrices c, cp
DoubleMatrix c(wf_.sd(ispin,ikp)->c());
DoubleMatrix cp(dwf.sd(ispin,ikp)->c());
// factor 2.0 in next line: G and -G
a.gemm('t','n',2.0,c,cp,0.0);
// rank-1 update correction
a.ger(-1.0,c,0,cp,0);
DoubleMatrix a(c.context(),c.n(),c.n(),c.nb(),c.nb());
// cp = cp - c * a
cp.gemm('n','n',-1.0,c,a,1.0);
}
else
{
// not implemented in the complex case
assert(false);
}
tmap_["psd_residual"].stop();
// factor 2.0 in next line: G and -G
a.gemm('t','n',2.0,c,cp,0.0);
// rank-1 update correction
a.ger(-1.0,c,0,cp,0);
// dwf.sd->c() now contains the descent direction (HV-VA)
// cp = cp - c * a
cp.gemm('n','n',-1.0,c,a,1.0);
}
else
{
ComplexMatrix c(wf_.sd(ispin,ikp)->c());
ComplexMatrix cp(dwf.sd(ispin,ikp)->c());
ComplexMatrix a(c.context(),c.n(),c.n(),c.nb(),c.nb());
a.gemm('c','n',1.0,c,cp,0.0);
// cp = cp - c * a
cp.gemm('n','n',-1.0,c,a,1.0);
}
tmap_["psd_residual"].stop();
tmap_["psd_update_wf"].start();
const valarray<double>& diag = prec_.diag(ispin,ikp);
// dwf.sd->c() now contains the descent direction (HV-VA)
double* coeff = (double*) wf_.sd(ispin,ikp)->c().valptr();
const double* dcoeff =
(const double*) dwf.sd(ispin,ikp)->c().cvalptr();
const int mloc = wf_.sd(ispin,ikp)->c().mloc();
const int ngwl = wf_.sd(ispin,ikp)->basis().localsize();
const int nloc = wf_.sd(ispin,ikp)->c().nloc();
for ( int n = 0; n < nloc; n++ )
{
// note: double mloc length for complex<double> indices
double* c = &coeff[2*mloc*n];
const double* dc = &dcoeff[2*mloc*n];
// loop to ngwl only since diag[i] is not defined on [0:mloc-1]
for ( int i = 0; i < ngwl; i++ )
{
const double fac = diag[i];
const double delta_re = fac * dc[2*i];
const double delta_im = fac * dc[2*i+1];
c[2*i] -= delta_re;
c[2*i+1] -= delta_im;
}
}
tmap_["psd_update_wf"].stop();
tmap_["psd_update_wf"].start();
const valarray<double>& diag = prec_.diag(ispin,ikp);
tmap_["gram"].start();
wf_.sd(ispin,ikp)->gram();
tmap_["gram"].stop();
double* coeff = (double*) wf_.sd(ispin,ikp)->c().valptr();
const double* dcoeff =
(const double*) dwf.sd(ispin,ikp)->c().cvalptr();
const int mloc = wf_.sd(ispin,ikp)->c().mloc();
const int ngwl = wf_.sd(ispin,ikp)->basis().localsize();
const int nloc = wf_.sd(ispin,ikp)->c().nloc();
for ( int n = 0; n < nloc; n++ )
{
// note: double mloc length for complex<double> indices
double* c = &coeff[2*mloc*n];
const double* dc = &dcoeff[2*mloc*n];
// loop to ngwl only since diag[i] is not defined on [0:mloc-1]
for ( int i = 0; i < ngwl; i++ )
{
const double fac = diag[i];
const double delta_re = fac * dc[2*i];
const double delta_im = fac * dc[2*i+1];
c[2*i] -= delta_re;
c[2*i+1] -= delta_im;
}
}
tmap_["psd_update_wf"].stop();
tmap_["gram"].start();
wf_.sd(ispin,ikp)->gram();
tmap_["gram"].stop();
}
}
}
src/Preconditioner.C
......@@ -3,7 +3,7 @@
// Preconditioner.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: Preconditioner.C,v 1.5 2007-10-19 16:24:04 fgygi Exp $
// $Id: Preconditioner.C,v 1.6 2007-10-19 17:37:06 fgygi Exp $
#include "Preconditioner.h"
#include "EnergyFunctional.h"
......@@ -35,30 +35,27 @@ void Preconditioner::update(void)
diag_[ispin].resize(wf.nkp());
for ( int ikp = 0; ikp < wf.nkp(); ikp++ )
{
if ( wf.sd(ispin,ikp) != 0 && wf.sdcontext(ispin,ikp)->active() )
{
// Only resize and initialize diag_ if ikp is active on this task
const Basis& basis = wf.sd(ispin,ikp)->basis();
const int ngwloc = basis.localsize();
diag_[ispin][ikp].resize(ngwloc);
const double *kpg2_ptr = basis.kpg2_ptr();
// Only resize and initialize diag_ if ikp is active on this task
const Basis& basis = wf.sd(ispin,ikp)->basis();
const int ngwloc = basis.localsize();
diag_[ispin][ikp].resize(ngwloc);
const double *kpg2_ptr = basis.kpg2_ptr();
if ( use_confinement )
if ( use_confinement )
{
const valarray<double>& fstress = ef_.confpot(ikp)->fstress();
for ( int ig = 0; ig < ngwloc; ig++ )
{
const valarray<double>& fstress = ef_.confpot(ikp)->fstress();
for ( int ig = 0; ig < ngwloc; ig++ )
{
double e = 0.5 * ( kpg2_ptr[ig] + fstress[ig] );
diag_[ispin][ikp][ig] = ( e < ecutpr ) ? 0.5 / ecutpr : 0.5 / e;
}
double e = 0.5 * ( kpg2_ptr[ig] + fstress[ig] );
diag_[ispin][ikp][ig] = ( e < ecutpr ) ? 0.5 / ecutpr : 0.5 / e;
}
else
}
else
{
for ( int ig = 0; ig < ngwloc; ig++ )
{
for ( int ig = 0; ig < ngwloc; ig++ )
{
double e = 0.5 * kpg2_ptr[ig];
diag_[ispin][ikp][ig] = ( e < ecutpr ) ? 0.5 / ecutpr : 0.5 / e;
}
double e = 0.5 * kpg2_ptr[ig];
diag_[ispin][ikp][ig] = ( e < ecutpr ) ? 0.5 / ecutpr : 0.5 / e;
}
}
}
......
src/SDWavefunctionStepper.C
......@@ -3,7 +3,7 @@
// SDWavefunctionStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: SDWavefunctionStepper.C,v 1.5 2007-10-19 16:24:04 fgygi Exp $
// $Id: SDWavefunctionStepper.C,v 1.6 2007-10-19 17:37:06 fgygi Exp $
#include "SDWavefunctionStepper.h"
#include "Wavefunction.h"
......@@ -25,19 +25,13 @@ void SDWavefunctionStepper::update(Wavefunction& dwf)
{
for ( int ikp = 0; ikp < wf_.nkp(); ikp++ )
{
if ( wf_.sd(ispin,ikp) != 0 )
{
if ( wf_.sdcontext(ispin,ikp)->active() )
{
// c = c - dt2bye * hpsi
tmap_["sd_update_wf"].start();
wf_.sd(ispin,ikp)->c().axpy(-alpha_,dwf.sd(ispin,ikp)->c());
tmap_["sd_update_wf"].stop();
tmap_["gram"].start();
wf_.sd(ispin,ikp)->gram();
tmap_["gram"].stop();
}
}
// c = c - dt2bye * hpsi
tmap_["sd_update_wf"].start();
wf_.sd(ispin,ikp)->c().axpy(-alpha_,dwf.sd(ispin,ikp)->c());
tmap_["sd_update_wf"].stop();
tmap_["gram"].start();
wf_.sd(ispin,ikp)->gram();
tmap_["gram"].stop();
}
}
}
src/SampleHandler.C
......@@ -3,7 +3,7 @@
// SampleHandler.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: SampleHandler.C,v 1.7 2007-10-19 16:24:04 fgygi Exp $
// $Id: SampleHandler.C,v 1.8 2007-10-19 17:37:06 fgygi Exp $
#if USE_XERCES
......@@ -19,8 +19,10 @@ using namespace std;
////////////////////////////////////////////////////////////////////////////////
SampleHandler::SampleHandler(Sample& s, DoubleMatrix& gfdata,
vector<vector<vector<double> > > &dmat,
Wavefunction& wfvtmp) :
s_(s), gfdata_(gfdata), read_wf(false), read_wfv(false), wfvtmp_(wfvtmp) {}
s_(s), gfdata_(gfdata), dmat_(dmat), read_wf(false), read_wfv(false),
wfvtmp_(wfvtmp) {}
////////////////////////////////////////////////////////////////////////////////
SampleHandler::~SampleHandler() {}
......@@ -59,12 +61,12 @@ StructureHandler* SampleHandler::startSubHandler(const XMLCh* const uri,
else if ( qnm == "wavefunction" )
{
read_wf = true;
return new WavefunctionHandler(s_.wf,gfdata_);
return new WavefunctionHandler(s_.wf,gfdata_,dmat_);
}
else if ( qnm == "wavefunction_velocity" )
{
read_wfv = true;
return new WavefunctionHandler(wfvtmp_,gfdata_);
return new WavefunctionHandler(wfvtmp_,gfdata_,dmat_);
}
else
{
......
src/SampleHandler.h
......@@ -3,13 +3,14 @@
// SampleHandler.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: SampleHandler.h,v 1.6 2007-10-19 16:24:04 fgygi Exp $
// $Id: SampleHandler.h,v 1.7 2007-10-19 17:37:06 fgygi Exp $
#ifndef SAMPLEHANDLER_H
#define SAMPLEHANDLER_H
#include "StructureHandler.h"
#include <string>
#include <vector>
class DoubleMatrix;
class Sample;
class Wavefunction;
......@@ -19,6 +20,7 @@ class SampleHandler : public StructureHandler
private:
Sample& s_;
std::vector<std::vector<std::vector<double> > > &dmat_;
DoubleMatrix& gfdata_;
Wavefunction& wfvtmp_;
......@@ -44,7 +46,9 @@ class SampleHandler : public StructureHandler
const XMLCh* const localname, const XMLCh* const qname,
const StructureHandler* const subHandler);
SampleHandler(Sample& s, DoubleMatrix& gfdata, Wavefunction& wfvtmp);
SampleHandler(Sample& s, DoubleMatrix& gfdata,
std::vector<std::vector<std::vector<double> > > &dmat,
Wavefunction& wfvtmp);
~SampleHandler();
};
#endif
src/SampleReader.h
......@@ -3,12 +3,16 @@
// SampleReader.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: SampleReader.h,v 1.4 2007-10-19 16:24:05 fgygi Exp $
// $Id: SampleReader.h,v 1.5 2007-10-19 17:37:06 fgygi Exp $
#ifndef SAMPLEREADER_H
#define SAMPLEREADER_H
enum event_type { species, atom, wavefunction, wavefunction_velocity,
slater_determinant, end, invalid };
class Context;
class Sample;
class SampleReader
{
......
src/SlaterDet.h
......@@ -3,7 +3,7 @@
// SlaterDet.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: SlaterDet.h,v 1.20 2007-10-19 16:24:05 fgygi Exp $
// $Id: SlaterDet.h,v 1.21 2007-10-19 17:37:06 fgygi Exp $
#ifndef SLATERDET_H
#define SLATERDET_H
......@@ -63,11 +63,11 @@ class SlaterDet
void cleanup(void);
void reset(void);
void gram(void);
void riccati(SlaterDet& sd);
void riccati(const SlaterDet& sd);
void lowdin(void);
void align(const SlaterDet& sd);
void ortho_align(const SlaterDet& sd);
double dot(const SlaterDet& sd) const;
std::complex<double> dot(const SlaterDet& sd) const;
double total_charge(void);
void update_occ(int nel, int nspin);
void update_occ(int nspin, double mu, double temp);
......@@ -91,9 +91,11 @@ class SlaterDet
double memsize(void) const;
double localmemsize(void) const;
SlaterDet& operator=(SlaterDet& rhs);
void print(std::ostream& os, std::string encoding);
void print(std::ostream& os, std::string encoding, double weight, int ispin,
int nspin);
#if USE_CSTDIO_LFS
void write(FILE* outfile, std::string encoding);
void write(FILE* outfile, std::string encoding, double weight, int ispin,
int nspin);
#endif
void info(std::ostream& os);
};
......
src/Wavefunction.h
......@@ -3,7 +3,7 @@
// Wavefunction.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: Wavefunction.h,v 1.16 2007-10-19 16:24:05 fgygi Exp $
// $Id: Wavefunction.h,v 1.17 2007-10-19 17:37:06 fgygi Exp $
#ifndef WAVEFUNCTION_H
#define WAVEFUNCTION_H
......@@ -11,6 +11,7 @@
#include "D3vector.h"
#include "UnitCell.h"
#include <vector>
#include <complex>
#if USE_CSTDIO_LFS
#include <cstdio>
#endif
......@@ -38,10 +39,11 @@ class Wavefunction
std::vector<double> weight_; // weight[ikp]
std::vector<D3vector> kpoint_; // kpoint[ikp]
std::vector<int> nst_; // nst_[ispin]
std::vector<Context*> spincontext_; // spincontext[ispin]
std::vector<std::vector<Context*> > sdcontext_; // sdcontext_[ispin][ikp]
std::vector<std::vector<SlaterDet*> > sd_; // sd[ispin][ikp]
std::vector<int> nst_; // nst_[ispin]
const Context* spincontext_; // context used for spin reductions
const Context* kpcontext_; // context used for kp reductions
const Context* sdcontext_; // context of local SlaterDet instances
std::vector<std::vector<SlaterDet*> > sd_; // local SlaterDets sd_[ispin][ikp]
void allocate(); // create contexts and allocate SlaterDet's
void deallocate();
......@@ -62,10 +64,10 @@ class Wavefunction
double weight(int ikp) const { return weight_[ikp]; }
double ecut(void) const { return ecut_; }
SlaterDet* sd(int ispin, int ikp) const { return sd_[ispin][ikp]; }
const Context* sdcontext(int ispin, int ikp) const
{ return sdcontext_[ispin][ikp]; }
const Context* spincontext(int ispin) const
{ return spincontext_[ispin]; }
const Context* spincontext(void) const { return spincontext_; }
const Context* kpcontext(void) const { return kpcontext_; }
const Context* sdcontext(void) const { return sdcontext_; }
int nkp(void) const; // number of k points
int nel(void) const; // total number of electrons
int nst(int ispin) const; // number of states of spin ispin
......@@ -98,7 +100,7 @@ class Wavefunction
void align(Wavefunction& wf);
void diag(Wavefunction& dwf, bool eigvec);
double dot(const Wavefunction& wf) const;
std::complex<double> dot(const Wavefunction& wf) const;
void print(std::ostream& os, std::string encoding, std::string tag) const;
#if USE_CSTDIO_LFS
......
src/WavefunctionHandler.h
......@@ -3,7 +3,7 @@
// WavefunctionHandler.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: WavefunctionHandler.h,v 1.9 2007-10-19 16:24:05 fgygi Exp $
// $Id: WavefunctionHandler.h,v 1.10 2007-10-19 17:37:06 fgygi Exp $
#ifndef WavefunctionHANDLER_H
#define WavefunctionHANDLER_H
......@@ -24,11 +24,15 @@ class WavefunctionHandler : public StructureHandler
UnitCell uc;
UnitCell ruc;
double ecut;
std::vector<double> dmat;
// dmat[ispin][ikp][i]
std::vector<std::vector<std::vector<double> > > &dmat_;
int nx,ny,nz;
int current_gf_nx,current_gf_ny,current_gf_nz;
std::string current_gf_encoding;
int current_ispin,current_ikp,current_n,current_igf;
std::vector<double> dmat_tmp;
double current_kx, current_ky, current_kz, current_weight;
int current_size;
int read_from_gfdata;
FourierTransform* ft;
std::vector<std::complex<double> > wftmp;
......@@ -55,7 +59,8 @@ class WavefunctionHandler : public StructureHandler
const XMLCh* const localname, const XMLCh* const qname,
const StructureHandler* const subHandler);
WavefunctionHandler(Wavefunction& wf, DoubleMatrix& gfdata);
WavefunctionHandler(Wavefunction& wf, DoubleMatrix& gfdata,
std::vector<std::vector<std::vector<double> > > &dmat);
~WavefunctionHandler();
};
#endif
src/XMLGFPreprocessor.C
......@@ -3,7 +3,7 @@
// XMLGFPreprocessor.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: XMLGFPreprocessor.C,v 1.10 2007-10-19 16:24:05 fgygi Exp $
// $Id: XMLGFPreprocessor.C,v 1.11 2007-10-19 17:37:06 fgygi Exp $
#include <cassert>
#include <iostream>
......@@ -711,7 +711,7 @@ void XMLGFPreprocessor::process(const char* const filename,
#if DEBUG
for ( int iseg = 0; iseg < seg_start.size(); iseg++ )
{
cout << rctxt.mype() << ": seg_data[iseg=" << iseg << "]: ";
cout << rctxt.mype() << ": seg_data[iseg=" << iseg << "]: size=" << dbuf[iseg].size() << " ";
if ( dbuf[iseg].size() >=3 )
cout << dbuf[iseg][0] << " " << dbuf[iseg][1]
<< " " << dbuf[iseg][2];
......@@ -752,9 +752,36 @@ void XMLGFPreprocessor::process(const char* const filename,
int gfnb = ngf / ctxt.npcol() +
( ngf % ctxt.npcol() != 0 ? 1 : 0 );
gfdata.resize(maxgfsize,ngf,gfmb,gfnb);
//cout << ctxt.mype() << ": gfdata resized: (" << maxgfsize << "x" << ngf
//<< ") (" << gfmb << "x" << gfnb << ") blocks" << endl;
//cout << ctxt.mype() << ": gfdata.context(): " << gfdata.context();
#if DEBUG
cout << ctxt.mype() << ": gfdata resized: (" << maxgfsize << "x" << ngf
<< ") (" << gfmb << "x" << gfnb << ") blocks" << endl;
cout << ctxt.mype() << ": gfdata.context(): " << gfdata.context();
#endif
// If the data in dbuf segments corresponds to real functions (kpoint==0)
// and if there are also complex functions (kpoint!=0), then
// for real functions: gfsize[igf[iseg]] == maxgfsize/2
// for complex functions: gfsize[igf[iseg] == maxgfsize
// If both real and complex functions are present, the real functions
// are converted to complex functions before redistribution of data
// This is done by resizing dbuf[iseg] and copying the data with stride 2
for ( int iseg = 0; iseg < seg_start.size(); iseg++ )
{
if ( maxgfsize == 2*gfsize[igf[iseg]] )
{
// the function is real and there are also complex functions
// resize dbuf[iseg] to double its size
const int ndoubles = dbuf[iseg].size();
dbuf[iseg].resize(2*ndoubles);
// redistribute the data in the array
for ( int i = ndoubles-1; i >= 0; i-- )
{
dbuf[iseg][2*i] = dbuf[iseg][i];
dbuf[iseg][2*i+1] = 0.0;
}
}
}
// prepare buffer sbuf for all_to_all operation
int sbufsize = 0;
......
src/notes
......@@ -33,6 +33,12 @@ with srand48(ctxt.myproc()).
Could reinitialize the seed in MDIonicsStepper.C using a common value identical
on all processors.
--------------------------------------------------------------------------------
rel1_34_0
Added kpoints
--------------------------------------------------------------------------------
rel1_33_5
removed trailing blanks in *.[Ch] *.mk Makefile
--------------------------------------------------------------------------------
rel1_33_4
ComputeMLWFCmd.C: fixed memory leak.
blas.h: added zgemm, dnrm2.
......
src/qb.C
......@@ -3,7 +3,7 @@
// qb.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: qb.C,v 1.55 2007-10-19 16:24:06 fgygi Exp $
// $Id: qb.C,v 1.56 2007-10-19 17:37:06 fgygi Exp $
#include <iostream>
#include <string>
......@@ -35,6 +35,7 @@ using namespace std;
#include "ConstraintCmd.h"
#include "DistanceCmd.h"
#include "HelpCmd.h"
#include "KpointCmd.h"
#include "ListAtomsCmd.h"
#include "ListSpeciesCmd.h"
#include "LoadCmd.h"
......@@ -217,6 +218,7 @@ int main(int argc, char **argv, char **envp)
ui.addCmd(new ConstraintCmd(s));
ui.addCmd(new DistanceCmd(s));
ui.addCmd(new HelpCmd(s));
ui.addCmd(new KpointCmd(s));
ui.addCmd(new ListAtomsCmd(s));
ui.addCmd(new ListSpeciesCmd(s));
ui.addCmd(new LoadCmd(s));
......
src/release.C
......@@ -3,10 +3,10 @@
// release.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: release.C,v 1.37 2007-10-16 18:23:21 fgygi Exp $
// $Id: release.C,v 1.38 2007-10-19 17:37:06 fgygi Exp $
#include "release.h"
std::string release(void)
{
return std::string("1.33.4");
return std::string("1.34.0");
}
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