Commit 16e264d0 by Francois Gygi

rel1_18_0


git-svn-id: http://qboxcode.org/svn/qb/trunk@251 cba15fb0-1239-40c8-b417-11db7ca47a34
parent 869fe59a
......@@ -3,7 +3,7 @@
// AtomCmd.h:
//
////////////////////////////////////////////////////////////////////////////////
// $Id: AtomCmd.h,v 1.4 2003-09-16 16:24:26 fgygi Exp $
// $Id: AtomCmd.h,v 1.5 2004-09-14 22:24:11 fgygi Exp $
#ifndef ATOMCMD_H
#define ATOMCMD_H
......@@ -79,7 +79,7 @@ cout << " dbg check "<< __FILE__ <<" "<< __LINE__ <<" mype="<< mype << endl;
#endif
s->wf.set_nel(s->atoms.nel());
s->wf.update_occ();
s->wf.update_occ(0.0);
if ( s->wfv != 0 )
{
s->wfv->set_nel(s->atoms.nel());
......
......@@ -3,7 +3,9 @@
// AtomSetHandler.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: AtomSetHandler.C,v 1.2 2003-05-23 21:51:04 fgygi Exp $
// $Id: AtomSetHandler.C,v 1.3 2004-09-14 22:24:11 fgygi Exp $
#if USE_XERCES
#include "AtomSetHandler.h"
#include "AtomSet.h"
......@@ -193,4 +195,4 @@ void AtomSetHandler::endSubHandler(const XMLCh* const uri,
delete last;
}
#endif
......@@ -3,7 +3,7 @@
// BLYPFunctional.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: BLYPFunctional.C,v 1.2 2003-08-22 18:01:13 fgygi Exp $
// $Id: BLYPFunctional.C,v 1.3 2004-09-14 22:24:11 fgygi Exp $
#include <cmath>
#include <cassert>
......@@ -11,7 +11,7 @@
using namespace std;
////////////////////////////////////////////////////////////////////////////////
BLYPFunctional::BLYPFunctional(vector<vector<double> > &rhoe)
BLYPFunctional::BLYPFunctional(const vector<vector<double> > &rhoe)
{
_nspin = rhoe.size();
if ( _nspin > 1 ) assert(rhoe[0].size() == rhoe[1].size());
......
......@@ -3,7 +3,7 @@
// BLYPFunctional.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: BLYPFunctional.h,v 1.1 2003-04-10 19:11:38 fgygi Exp $
// $Id: BLYPFunctional.h,v 1.2 2004-09-14 22:24:11 fgygi Exp $
#ifndef BLYPFUNCTIONAL_H
#define BLYPFUNCTIONAL_H
......@@ -32,7 +32,7 @@ class BLYPFunctional : public XCFunctional
public:
BLYPFunctional(vector<vector<double> > &rhoe);
BLYPFunctional(const vector<vector<double> > &rhoe);
bool isGGA() { return true; };
string name() { return "BLYP"; };
......
......@@ -3,15 +3,15 @@
// BOSampleStepper.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: BOSampleStepper.h,v 1.2 2004-03-11 21:52:32 fgygi Exp $
// $Id: BOSampleStepper.h,v 1.3 2004-09-14 22:24:11 fgygi Exp $
#ifndef BOSAMPLESTEPPER_H
#define BOSAMPLESTEPPER_H
#include "SampleStepper.h"
#include "EnergyFunctional.h"
#include "Sample.h"
#include "Wavefunction.h"
class EnergyFunctional;
class WavefunctionStepper;
class IonicStepper;
using namespace std;
......@@ -22,8 +22,10 @@ class BOSampleStepper : public SampleStepper
Wavefunction dwf;
Wavefunction* wfv;
int nitscf_;
int nite_;
EnergyFunctional& ef_;
ChargeDensity cd_;
EnergyFunctional ef_;
WavefunctionStepper* wf_stepper;
IonicStepper* ionic_stepper;
......@@ -37,7 +39,7 @@ class BOSampleStepper : public SampleStepper
void step(int niter);
BOSampleStepper(Sample& s, EnergyFunctional& ef, int nite);
BOSampleStepper(Sample& s, int nitscf, int nite);
//~BOSampleStepper();
};
#endif
......@@ -3,10 +3,9 @@
// CPSampleStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: CPSampleStepper.C,v 1.6 2004-05-04 21:26:24 fgygi Exp $
// $Id: CPSampleStepper.C,v 1.7 2004-09-14 22:24:11 fgygi Exp $
#include "CPSampleStepper.h"
#include "EnergyFunctional.h"
#include "SlaterDet.h"
#include "MDWavefunctionStepper.h"
#include "MDIonicStepper.h"
......@@ -18,8 +17,8 @@
using namespace std;
////////////////////////////////////////////////////////////////////////////////
CPSampleStepper::CPSampleStepper(Sample& s, EnergyFunctional& ef) :
SampleStepper(s), ef_(ef), dwf(s.wf), wfv(s.wfv)
CPSampleStepper::CPSampleStepper(Sample& s) :
SampleStepper(s), cd_(s.wf), ef_(s,cd_), dwf(s.wf), wfv(s.wfv)
{
mdwf_stepper = new MDWavefunctionStepper(s,tmap);
assert(mdwf_stepper!=0);
......@@ -66,6 +65,8 @@ void CPSampleStepper::step(int niter)
Timer tm_iter;
cd_.update_density();
ef_.update_vhxc();
double energy =
ef_.energy(compute_hpsi,dwf,compute_forces,fion,compute_stress,sigma_eks);
......@@ -196,6 +197,8 @@ void CPSampleStepper::step(int niter)
ef_.atoms_moved();
}
cd_.update_density();
ef_.update_vhxc();
energy =
ef_.energy(compute_hpsi,dwf,compute_forces,fion,compute_stress,sigma_eks);
......
......@@ -3,15 +3,16 @@
// CPSampleStepper.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: CPSampleStepper.h,v 1.2 2004-03-11 21:52:31 fgygi Exp $
// $Id: CPSampleStepper.h,v 1.3 2004-09-14 22:24:11 fgygi Exp $
#ifndef CPSAMPLESTEPPER_H
#define CPSAMPLESTEPPER_H
#include "SampleStepper.h"
#include "EnergyFunctional.h"
#include "ChargeDensity.h"
#include "Sample.h"
#include "Wavefunction.h"
class EnergyFunctional;
class MDWavefunctionStepper;
class MDIonicStepper;
using namespace std;
......@@ -20,7 +21,8 @@ class CPSampleStepper : public SampleStepper
{
private:
EnergyFunctional& ef_;
ChargeDensity cd_;
EnergyFunctional ef_;
Wavefunction dwf;
Wavefunction* wfv;
......@@ -36,7 +38,7 @@ class CPSampleStepper : public SampleStepper
void step(int niter);
CPSampleStepper(Sample& s, EnergyFunctional& ef);
CPSampleStepper(Sample& s);
~CPSampleStepper();
};
#endif
......@@ -3,7 +3,7 @@
// ChargeDensity.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: ChargeDensity.C,v 1.7 2004-08-11 17:56:24 fgygi Exp $
// $Id: ChargeDensity.C,v 1.8 2004-09-14 22:24:11 fgygi Exp $
#include "ChargeDensity.h"
#include "Basis.h"
......@@ -93,6 +93,7 @@ ChargeDensity::~ChargeDensity(void)
for ( int ikp = 0; ikp < ft_.size(); ikp++ )
delete ft_[ikp];
}
////////////////////////////////////////////////////////////////////////////////
void ChargeDensity::update_density(void)
{
......@@ -149,5 +150,33 @@ void ChargeDensity::update_density(void)
}
}
}
////////////////////////////////////////////////////////////////////////////////
void ChargeDensity::update_rhor(void)
{
// recalculate rhor from rhog
assert(rhor.size() == wf_.nspin());
const double omega = vbasis_->cell().volume();
assert(omega!=0.0);
const double omega_inv = 1.0 / omega;
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];
#pragma ivdep
for ( int i = 0; i < rhor_size; i++ )
{
prhor[i] = rhotmp[i].real() * omega_inv;
}
}
}
}
......@@ -3,7 +3,7 @@
// ChargeDensity.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: ChargeDensity.h,v 1.1 2003-01-25 01:23:31 fgygi Exp $
// $Id: ChargeDensity.h,v 1.2 2004-09-14 22:24:11 fgygi Exp $
#ifndef CHARGEDENSITY_H
#define CHARGEDENSITY_H
......@@ -35,6 +35,7 @@ class ChargeDensity
vector<vector<complex<double> > > rhog; // rhog[ispin][ig]
void update_density(void);
void update_rhor(void);
Basis* vbasis(void) const { return vbasis_; }
const Context* vcontext(void) const { return vcontext_; }
......
......@@ -3,7 +3,7 @@
// Control.h:
//
////////////////////////////////////////////////////////////////////////////////
// $Id: Control.h,v 1.7 2004-03-11 21:52:32 fgygi Exp $
// $Id: Control.h,v 1.8 2004-09-14 22:24:11 fgygi Exp $
#ifndef CONTROL_H
#define CONTROL_H
......@@ -19,7 +19,6 @@ struct Control
int nite;
double emass; // electron mass
string fermi; // use Fermi distribution to fill states
double fermi_temp; // temperature of Fermi distribution
double ecutprec;
......
......@@ -3,7 +3,7 @@
// EnergyFunctional.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: EnergyFunctional.h,v 1.11 2004-03-11 21:52:32 fgygi Exp $
// $Id: EnergyFunctional.h,v 1.12 2004-09-14 22:24:11 fgygi Exp $
#ifndef ENERGYFUNCTIONAL_H
#define ENERGYFUNCTIONAL_H
......@@ -35,7 +35,7 @@ class EnergyFunctional
private:
const Sample& s_;
ChargeDensity cd_;
const ChargeDensity& cd_;
Basis* vbasis_;
FourierTransform *vft;
vector<FourierTransform*> ft;
......@@ -48,7 +48,6 @@ class EnergyFunctional
vector<complex<double> > tmp_r, vion_local_g, dvion_local_g, vlocal_g,
rhopst, rhogt, rhoelg, vtemp;
vector<double> ftmp;
vector<vector<double> > v_r;
vector<vector<double> > tau0, taum, fion_esr;
vector<double> zv_, rcps_;
......@@ -62,6 +61,7 @@ class EnergyFunctional
public:
vector<vector<double> > v_r;
mutable TimerMap tmap;
double energy(bool compute_hpsi, Wavefunction& dwf,
......@@ -81,12 +81,14 @@ class EnergyFunctional
const ConfinementPotential *confpot(int ikp) const { return cfp[ikp]; }
void update_vhxc(void);
void atoms_moved(void);
void cell_moved(void);
void print(ostream& os) const;
EnergyFunctional(const Sample& s);
EnergyFunctional(const Sample& s, const ChargeDensity& cd);
~EnergyFunctional();
};
ostream& operator << ( ostream& os, const EnergyFunctional& e );
......
......@@ -3,7 +3,7 @@
// LDAFunctional.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: LDAFunctional.h,v 1.2 2003-03-27 22:05:59 fgygi Exp $
// $Id: LDAFunctional.h,v 1.3 2004-09-14 22:24:11 fgygi Exp $
#ifndef LDAFUNCTIONAL_H
#define LDAFUNCTIONAL_H
......@@ -24,7 +24,7 @@ class LDAFunctional : public XCFunctional
public:
LDAFunctional(vector<vector<double> > &rhoe)
LDAFunctional(const vector<vector<double> > &rhoe)
{
_nspin = rhoe.size();
if ( _nspin > 1 ) assert(rhoe[0].size() == rhoe[1].size());
......
......@@ -3,7 +3,7 @@
// LoadCmd.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: LoadCmd.C,v 1.5 2003-10-02 17:37:49 fgygi Exp $
// $Id: LoadCmd.C,v 1.6 2004-09-14 22:24:11 fgygi Exp $
#include "LoadCmd.h"
#include "SampleReader.h"
......@@ -58,7 +58,7 @@ int LoadCmd::action(int argc, char **argv)
if ( s->wf.nel() != s->atoms.nel() )
{
s->wf.set_nel(s->atoms.nel());
s->wf.update_occ();
s->wf.update_occ(0.0);
}
//cout << " LoadCmd: atoms.nel() = " << s->atoms.nel() << endl;
//cout << " LoadCmd: wf.nel() = " << s->wf.nel() << endl;
......
#-------------------------------------------------------------------------------
# $Id: Makefile,v 1.26 2004-08-11 17:56:24 fgygi Exp $
# $Id: Makefile,v 1.27 2004-09-14 22:24:11 fgygi Exp $
#------------------------------------------------------------------------------
#
include $(TARGET).mk
......@@ -9,8 +9,7 @@ EXECS=qb testMatrix testFTGrid testGridFunction testBasis \
testBlacsContext testSlaterDet testEnergyFunctional testSample \
testChargeDensity testFourierTransform testSpecies testContext
CXXFLAGS += -DTARGET='"$(TARGET)"'
qb: qb.o AtomSet.o Atom.o Species.o \
OBJECTS=qb.o AtomSet.o Atom.o Species.o \
Wavefunction.o SlaterDet.o \
EnergyFunctional.o SampleStepper.o \
Basis.o FourierTransform.o Matrix.o Context.o \
......@@ -29,7 +28,11 @@ CXXFLAGS += -DTARGET='"$(TARGET)"'
PSDWavefunctionStepper.o PSDAWavefunctionStepper.o \
SDCellStepper.o ConfinementPotential.o Preconditioner.o \
release.o qbox_xmlns.o isodate.o $(PLTOBJECTS)
CXXFLAGS += -DTARGET='"$(TARGET)"'
qb: $(OBJECTS)
$(LD) $(DFLAGS) -o $@ $^ $(LDFLAGS)
lib: $(OBJECTS)
ar cq libqb.a $^
SamplePrint: SamplePrint.o SamplePrintHandlers.o
$(LD) $(DFLAGS) -o $@ $^ $(LDFLAGS)
testSample: testSample.o AtomSet.o Atom.o Species.o \
......@@ -157,7 +160,8 @@ BOSampleStepper.o: SDIonicStepper.h IonicStepper.h MDIonicStepper.h
BOSampleStepper.o: SDCellStepper.h CellStepper.h Preconditioner.h
BOSampleStepper.o: SampleStepper.h Sample.h AtomSet.h Context.h Atom.h
BOSampleStepper.o: D3vector.h Species.h Wavefunction.h UnitCell.h Control.h
BOSampleStepper.o: Timer.h
BOSampleStepper.o: Timer.h EnergyFunctional.h ChargeDensity.h
BOSampleStepper.o: StructureFactor.h
CellDyn.o: Sample.h AtomSet.h Context.h Atom.h D3vector.h Species.h
CellDyn.o: Wavefunction.h UnitCell.h Control.h SlaterDet.h Basis.h Matrix.h
CellDyn.o: Timer.h
......@@ -184,7 +188,8 @@ CPSampleStepper.o: MDIonicStepper.h IonicStepper.h SDCellStepper.h
CPSampleStepper.o: CellStepper.h
CPSampleStepper.o: SampleStepper.h Sample.h AtomSet.h Context.h Atom.h
CPSampleStepper.o: D3vector.h Species.h Wavefunction.h UnitCell.h Control.h
CPSampleStepper.o: Timer.h
CPSampleStepper.o: Timer.h EnergyFunctional.h ChargeDensity.h
CPSampleStepper.o: StructureFactor.h
Debug.o: Sample.h AtomSet.h Context.h Atom.h D3vector.h Species.h
Debug.o: Wavefunction.h UnitCell.h Control.h
Dt.o: Sample.h AtomSet.h Context.h Atom.h D3vector.h Species.h Wavefunction.h
......@@ -207,6 +212,8 @@ EnergyFunctional.o: ConfinementPotential.h blas.h
EnergyFunctional.o: ChargeDensity.h Context.h StructureFactor.h Timer.h
ExtStress.o: Sample.h AtomSet.h Context.h Atom.h D3vector.h Species.h
ExtStress.o: Wavefunction.h UnitCell.h Control.h
FermiTemp.o: Sample.h AtomSet.h Context.h Atom.h D3vector.h Species.h
FermiTemp.o: Wavefunction.h UnitCell.h Control.h
FourierTransform.o: FourierTransform.h Basis.h D3vector.h UnitCell.h
FourierTransform.o: Context.h Timer.h
HelpCmd.o: UserInterface.h Sample.h AtomSet.h Context.h Atom.h D3vector.h
......@@ -272,14 +279,16 @@ PSDWavefunctionStepper.o: Preconditioner.h
PSDWavefunctionStepper.o: WavefunctionStepper.h Sample.h AtomSet.h Context.h
PSDWavefunctionStepper.o: Atom.h D3vector.h Species.h Wavefunction.h
PSDWavefunctionStepper.o: UnitCell.h Control.h Timer.h
qb.o: isodate.h release.h Context.h UserInterface.h Sample.h AtomSet.h Atom.h
qb.o: D3vector.h Species.h Wavefunction.h UnitCell.h Control.h Timer.h
qb.o: AtomCmd.h HelpCmd.h ListAtomsCmd.h ListSpeciesCmd.h LoadCmd.h
qb.o: PrintCmd.h QuitCmd.h RandomizeWfCmd.h RunCmd.h SaveCmd.h SetCmd.h
qb.o: SpeciesCmd.h StatusCmd.h AtomsDyn.h Cell.h CellDyn.h SlaterDet.h
qb.o: Basis.h Matrix.h CellLock.h CellMass.h Debug.h Ecut.h Ecutprec.h
qb.o: Ecuts.h Emass.h ExtStress.h Dt.h Nempty.h Nrowmax.h RefCell.h Stress.h
qb.o: Thermostat.h ThTemp.h ThTime.h ThWidth.h WfDiag.h WfDyn.h Xc.h
qb.o: isodate.h release.h qbox_xmlns.h Context.h UserInterface.h Sample.h
qb.o: AtomSet.h Atom.h D3vector.h Species.h Wavefunction.h UnitCell.h
qb.o: Control.h Timer.h AtomCmd.h HelpCmd.h ListAtomsCmd.h ListSpeciesCmd.h
qb.o: LoadCmd.h PrintCmd.h QuitCmd.h RandomizeWfCmd.h RunCmd.h SaveCmd.h
qb.o: SetCmd.h SpeciesCmd.h StatusCmd.h AtomsDyn.h Cell.h CellDyn.h
qb.o: SlaterDet.h Basis.h Matrix.h CellLock.h CellMass.h Debug.h Ecut.h
qb.o: Ecutprec.h Ecuts.h Emass.h ExtStress.h FermiTemp.h Dt.h Nempty.h
qb.o: Nrowmax.h RefCell.h Stress.h Thermostat.h ThTemp.h ThTime.h ThWidth.h
qb.o: WfDiag.h WfDyn.h Xc.h
qbox_xmlns.o: qbox_xmlns.h
QuitCmd.o: UserInterface.h Sample.h AtomSet.h Context.h Atom.h D3vector.h
QuitCmd.o: Species.h Wavefunction.h UnitCell.h Control.h
RandomizeWfCmd.o: UserInterface.h Sample.h AtomSet.h Context.h Atom.h
......@@ -289,8 +298,8 @@ RefCell.o: Wavefunction.h UnitCell.h Control.h
release.o: release.h
RunCmd.o: RunCmd.h UserInterface.h Sample.h AtomSet.h Context.h Atom.h
RunCmd.o: D3vector.h Species.h Wavefunction.h UnitCell.h Control.h
RunCmd.o: EnergyFunctional.h ChargeDensity.h StructureFactor.h Timer.h
RunCmd.o: BOSampleStepper.h SampleStepper.h CPSampleStepper.h
RunCmd.o: BOSampleStepper.h SampleStepper.h Timer.h EnergyFunctional.h
RunCmd.o: ChargeDensity.h StructureFactor.h CPSampleStepper.h
RunCmd.o: UserInterface.h Sample.h AtomSet.h Context.h Atom.h D3vector.h
RunCmd.o: Species.h Wavefunction.h UnitCell.h Control.h
Sample.o: AtomSet.h Context.h Atom.h D3vector.h Species.h Wavefunction.h
......@@ -313,7 +322,7 @@ SampleStepper.o: Sample.h AtomSet.h Context.h Atom.h D3vector.h Species.h
SampleStepper.o: Wavefunction.h UnitCell.h Control.h Timer.h
SaveCmd.o: SaveCmd.h UserInterface.h Sample.h AtomSet.h Context.h Atom.h
SaveCmd.o: D3vector.h Species.h Wavefunction.h UnitCell.h Control.h isodate.h
SaveCmd.o: release.h
SaveCmd.o: release.h qbox_xmlns.h
SaveCmd.o: UserInterface.h Sample.h AtomSet.h Context.h Atom.h D3vector.h
SaveCmd.o: Species.h Wavefunction.h UnitCell.h Control.h
SDCellStepper.o: SDCellStepper.h CellStepper.h Sample.h AtomSet.h Context.h
......@@ -392,11 +401,13 @@ testSlaterDet.o: Timer.h FourierTransform.h
testSpecies.o: Species.h Context.h SpeciesReader.h
testStepper.o: Sample.h AtomSet.h Context.h Atom.h D3vector.h Species.h
testStepper.o: Wavefunction.h UnitCell.h Control.h BOSampleStepper.h
testStepper.o: SampleStepper.h Timer.h
testStepper.o: SampleStepper.h Timer.h EnergyFunctional.h ChargeDensity.h
testStepper.o: StructureFactor.h
testUnitCell.o: UnitCell.h D3vector.h
testWavefunction.o: Context.h Wavefunction.h D3vector.h UnitCell.h
testWavefunction.o: SlaterDet.h Basis.h Matrix.h Timer.h
testwrite.o: Context.h
testXCFunctional.o: LDAFunctional.h XCFunctional.h PBEFunctional.h Timer.h
testXMLGFPreprocessor.o: Context.h Matrix.h XMLGFPreprocessor.h
Thermostat.o: Sample.h AtomSet.h Context.h Atom.h D3vector.h Species.h
Thermostat.o: Wavefunction.h UnitCell.h Control.h
......@@ -440,3 +451,4 @@ xmlget.o: StrX.h
XMLGFPreprocessor.o: Timer.h Context.h Base64Transcoder.h Matrix.h
XMLGFPreprocessor.o: XMLGFPreprocessor.h
XMLGFPreprocessor.o: Matrix.h
xmlSpecies.o: qbox_xmlns.h
......@@ -3,7 +3,7 @@
// PBEFunctional.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: PBEFunctional.C,v 1.4 2004-05-20 00:19:26 fgygi Exp $
// $Id: PBEFunctional.C,v 1.5 2004-09-14 22:24:11 fgygi Exp $
#include "PBEFunctional.h"
#include <cmath>
......@@ -12,7 +12,7 @@
#include <vector>
using namespace std;
PBEFunctional::PBEFunctional(vector<vector<double> > &rhoe)
PBEFunctional::PBEFunctional(const vector<vector<double> > &rhoe)
{
_nspin = rhoe.size();
if ( _nspin > 1 ) assert(rhoe[0].size() == rhoe[1].size());
......
......@@ -3,7 +3,7 @@
// PBEFunctional.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: PBEFunctional.h,v 1.1 2003-01-25 01:23:31 fgygi Exp $
// $Id: PBEFunctional.h,v 1.2 2004-09-14 22:24:11 fgygi Exp $
#ifndef PBEFUNCTIONAL_H
#define PBEFUNCTIONAL_H
......@@ -36,7 +36,7 @@ class PBEFunctional : public XCFunctional
public:
PBEFunctional(vector<vector<double> > &rhoe);
PBEFunctional(const vector<vector<double> > &rhoe);
bool isGGA() { return true; };
string name() { return "PBE"; };
......
......@@ -3,12 +3,11 @@
// RunCmd.C:
//
////////////////////////////////////////////////////////////////////////////////
// $Id: RunCmd.C,v 1.4 2004-03-11 21:52:31 fgygi Exp $
// $Id: RunCmd.C,v 1.5 2004-09-14 22:24:11 fgygi Exp $
#include "RunCmd.h"
#include<iostream>
using namespace std;
#include "EnergyFunctional.h"
#include "BOSampleStepper.h"
#include "CPSampleStepper.h"
......@@ -18,41 +17,48 @@ using namespace std;
int RunCmd::action(int argc, char **argv)
{
if ( argc < 2 || argc > 3)
if ( argc < 2 || argc > 4)
{
if ( ui->onpe0() )
cout << " use: run niter" << endl;
cout << " run niter nitscf" << endl;
cout << " run niter nitscf nite" << endl;
return 1;
}
if ( s->wf.nst() == 0 )
{
if ( ui->onpe0() )
cout << " <qb:error> RunCmd: no states, cannot run </qb:error>" << endl;
cout << " <!-- RunCmd: no states, cannot run -->" << endl;
return 1;
}
if ( s->wf.ecut() == 0.0 )
{
if ( ui->onpe0() )
cout << " <qb:error> RunCmd: ecut = 0.0, cannot run </qb:error>" << endl;
cout << " <!-- RunCmd: ecut = 0.0, cannot run -->" << endl;
return 1;
}
EnergyFunctional ef(*s);
SampleStepper* stepper;
int niter = atoi(argv[1]);
int nite = 1;
int nitscf = 1;
if ( argc == 3 )
{
// run niter nite
nite = atoi(argv[2]);
// run niter nitscf
nitscf = atoi(argv[2]);
}
else if ( argc == 4 )
{
// run niter nitscf nite
nitscf = atoi(argv[2]);
nite = atoi(argv[3]);
}
if ( s->ctrl.wf_dyn == "MD" )
stepper = new CPSampleStepper(*s,ef);
stepper = new CPSampleStepper(*s);
else
stepper = new BOSampleStepper(*s,ef,nite);
stepper = new BOSampleStepper(*s,nitscf,nite);
assert(stepper!=0);
......
......@@ -3,7 +3,9 @@
// SampleHandler.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: SampleHandler.C,v 1.4 2003-09-23 19:03:21 fgygi Exp $
// $Id: SampleHandler.C,v 1.5 2004-09-14 22:24:11 fgygi Exp $
#if USE_XERCES
#include "SampleHandler.h"
#include "Sample.h"
......@@ -78,3 +80,4 @@ void SampleHandler::endSubHandler(const XMLCh* const uri,
// cout << " SampleHandler::endSubHandler " << StrX(qname) << endl;
delete subHandler;
}
#endif
......@@ -3,14 +3,12 @@
// SampleReader.C:
//
////////////////////////////////////////////////////////////////////////////////
// $Id: SampleReader.C,v 1.14 2004-03-19 00:56:04 fgygi Exp $
// $Id: SampleReader.C,v 1.15 2004-09-14 22:24:11 fgygi Exp $
#include "Sample.h"
#include "SampleReader.h"
#include "SpeciesReader.h"
#include "StructuredDocumentHandler.h"
#include "SampleHandler.h"
#include "Basis.h"
#include "FourierTransform.h"
#include "SlaterDet.h"
......@@ -25,6 +23,9 @@
#include <sys/stat.h>
using namespace std;
#if USE_XERCES
#include "SampleHandler.h"
#include "StructuredDocumentHandler.h"
#include <xercesc/util/XMLUniDefs.hpp>
#include <xercesc/sax2/Attributes.hpp>
#include <xercesc/util/PlatformUtils.hpp>
......@@ -33,6 +34,7 @@ using namespace std;
#include <xercesc/sax2/XMLReaderFactory.hpp>
#include <xercesc/framework/MemBufInputSource.hpp>
using namespace xercesc;
#endif
////////////////////////////////////////////////////////////////////////////////
SampleReader::SampleReader(const Context& ctxt) : ctxt_(ctxt) {}
......@@ -40,6 +42,7 @@ SampleReader::SampleReader(const Context& ctxt) : ctxt_(ctxt) {}
////////////////////////////////////////////////////////////////////////////////
void SampleReader::readSample (Sample& s, const string uri, bool serial)
{
#if USE_XERCES
const char* encodingName = "UTF-8";
//SAX2XMLReader::ValSchemes valScheme = SAX2XMLReader::Val_Auto;
SAX2XMLReader::ValSchemes valScheme = SAX2XMLReader::Val_Always;
......@@ -524,4 +527,12 @@ void SampleReader::readSample (Sample& s, const string uri, bool serial)
s.wfv = new Wavefunction(s.wf);
*s.wfv = wfvtmp;
}
#else
// USE_XERCES was not defined
if ( ctxt_.onpe0() )
{
cout << " <!-- SampleReader: could not read (parser not defined) -->"
<< endl;
}
#endif
}
......@@ -3,12 +3,17 @@
// SlaterDet.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: SlaterDet.C,v 1.27 2004-08-11 17:56:24 fgygi Exp $
// $Id: SlaterDet.C,v 1.28 2004-09-14 22:24:11 fgygi Exp $
#include "SlaterDet.h"
#include "FourierTransform.h"
#include "Context.h"
#include "blas.h" // daxpy
#include "Base64Transcoder.h"
// // XML transcoding for print function
// #include <xercesc/util/Base64.hpp>
// #include <xercesc/util/XMLString.hpp>
// using namespace xercesc;
#include <cstdlib>
#include <iostream>
......@@ -19,11 +24,6 @@
#endif
using namespace std;
// XML transcoding for print function
#include <xercesc/util/Base64.hpp>
#include <xercesc/util/XMLString.hpp>
using namespace xercesc;
////////////////////////////////////////////////////////////////////////////////
SlaterDet::SlaterDet(const Context& ctxt, D3vector kpoint) : ctxt_(ctxt),
basis_(Context(ctxt,'c',ctxt.mycol()),kpoint), c_(ctxt)
......@@ -40,23 +40,6 @@ SlaterDet::~SlaterDet()
}
////////////////////////////////////////////////////////////////////////////////
void SlaterDet::byteswap_double(size_t n, double* x)
{
if (n==0) return;
unsigned char* c = (unsigned char*) x;
while ( n-- > 0 )
{
unsigned char tmp;
tmp = c[7]; c[7] = c[0]; c[0] = tmp;
tmp = c[6]; c[6] = c[1]; c[1] = tmp;
tmp = c[5]; c[5] = c[2]; c[2] = tmp;
tmp = c[4]; c[4] = c[3]; c[3] = tmp;
c+=8;
}
}
////////////////////////////////////////////////////////////////////////////////
void SlaterDet::resize(const UnitCell& cell, const UnitCell& refcell,
double ecut, int nst)
{
......@@ -237,7 +220,7 @@ void SlaterDet::compute_density(FourierTransform& ft,
////////////////////////////////////////////////////////////////////////////////
void SlaterDet::rs_mul_add(FourierTransform& ft,
double* v, SlaterDet& sdp) const
const double* v, SlaterDet& sdp) const
{
// transform states to real space, multiply states by v[r] in real space
// transform back to reciprocal space and add to sdp
......@@ -686,7 +669,8 @@ void SlaterDet::align(const SlaterDet& sd)
// Compute the distance | c - sdc | before alignment
for ( int i = 0; i < c_proxy.size(); i++ )
c_tmp_proxy[i] = c_proxy[i] - sdc_proxy[i];
cout << " SlaterDet::align: distance before: "<< c_tmp_proxy.nrm2() << endl;
//cout << " SlaterDet::align: distance before: "
// << c_tmp_proxy.nrm2() << endl;
// compute the polar decomposition of B
double diff = 1.0;
......@@ -708,10 +692,10 @@ void SlaterDet::align(const SlaterDet& sd)
// Next lines: use t as temporary to compute || x - xp ||_F
for ( int i = 0; i < t.size(); i++ )
t[i] = x[i] - xp[i];
//for ( int i = 0; i < t.size(); i++ )
// t[i] = x[i] - xp[i];
diff = t.nrm2();
//diff = t.nrm2();
//cout << " SlaterDet::align: diff=" << diff << endl;
......@@ -731,9 +715,10 @@ void SlaterDet::align(const SlaterDet& sd)
c_proxy.gemm('n','n',1.0,c_tmp_proxy,x,0.0);
// Compute the distance | c - sdc | after alignment
for ( int i = 0; i < c_proxy.size(); i++ )
c_tmp_proxy[i] = c_proxy[i] - sdc_proxy[i];