Commit 612ee413 by Francois Gygi

Merge branch 'develop'

parents e3136eb6 a23e6450
This diff is collapsed. Click to expand it.
Qbox First-Principles Molecular Dynamics
========================================
Qbox is a C++/MPI/OpenMP implementation of First-Principles Molecular Dynamics.
It implements electronic structure calculations within the framework
of Density Functional Theory, using a plane-wave basis set and
norm-conserving pseudopotentials.
For documentation, examples and installation instructions
see http://qboxcode.org
Qbox is developed in the Gygi group at the University of California Davis.
Qbox is distributed under the terms of the GNU General Public License
as published by the Free Software Foundation, either version 2 of
the License, or (at your option) any later version.
See the file COPYING in the root directory of this distribution
or <http://www.gnu.org/licenses/>.
Qbox was written by Francois Gygi.
Ivan Duchemin developed the Hartree-Fock exchange implementation and the
client-server interface.
Jun Wu developed the implementation of van der Waals density functionals.
Quan (Andy) Wan implemented the FFTW3 interface.
William Dawson contributed to the implementation of Hartree-Fock exchange.
Martin Schlipf implemented Optimized Norm-Conserving Vanderbilt potentials.
Martin Schlipf implemented the HSE hybrid density functional
Quan (Andy) Wan implemented the electric field.
Ma He contributed to the implementation of the response functionality.
......@@ -136,6 +136,14 @@ bool Basis::factorizable(int n) const
int Basis::np(int i) const { return np_[i]; }
////////////////////////////////////////////////////////////////////////////////
bool Basis::fits_in_grid(int np0, int np1, int np2) const
{ return
( idxmax_[0] < np0/2 ) && ( idxmin_[0] >= -np0/2 ) &&
( idxmax_[1] < np1/2 ) && ( idxmin_[1] >= -np1/2 ) &&
( idxmax_[2] < np2/2 ) && ( idxmin_[2] >= -np2/2 );
}
////////////////////////////////////////////////////////////////////////////////
const D3vector Basis::kpoint(void) const { return kpoint_; }
////////////////////////////////////////////////////////////////////////////////
......
......@@ -80,6 +80,7 @@ class Basis
const UnitCell& refcell() const;// reference cell dimensions
const D3vector kpoint() const; // k-point in units of b0,b1,b2
int np(int i) const; // good size of FFT grid in direction i
bool fits_in_grid(int np0, int np1, int np2) const;
bool factorizable(int n) const;// check if n is factorizable with low factors
int idxmin(int i) const; // smallest index in direction i
int idxmax(int i) const; // largest index in direction i
......
......@@ -74,6 +74,15 @@ class BisectionCmd : public Cmd
return 1;
}
if ( wf.nkp() > 1 )
{
if ( ui->onpe0() )
{
cout << " BisectionCmd: only implemented for k=0" << endl;
}
return 1;
}
if ( nLevels[0] < 0 || nLevels[0] > 5 ||
nLevels[1] < 0 || nLevels[1] > 5 ||
nLevels[2] < 0 || nLevels[2] > 5 )
......
......@@ -73,9 +73,10 @@ EnergyFunctional::EnergyFunctional( Sample& s, ChargeDensity& cd)
if ( s_.ctxt_.onpe0() )
{
cout << " EnergyFunctional: np0v,np1v,np2v: " << np0v << " "
<< np1v << " " << np2v << endl;
cout << " EnergyFunctional: vft->np012(): "
cout << " EnergyFunctional: <np0v> " << np0v << " </np0v> "
<< "<np1v> " << np1v << " </np1v> "
<< "<np2v> " << np2v << " </np2v>" << endl;
cout << " EnergyFunctional: vft->np012(): "
<< vft->np012() << endl;
}
......
......@@ -248,6 +248,10 @@ FourierTransform::FourierTransform (const Basis &basis,
rdispl[iproc] = rdispl[iproc-1] + rcounts[iproc-1];
}
// check if the basis_ fits in the grid np0, np1, np2
basis_fits_in_grid_ = basis_.fits_in_grid(np0,np1,np2);
assert(basis_fits_in_grid_);
if ( basis_.real() )
{
// compute index arrays ifftp_ and ifftm_ for mapping vector->zvec
......@@ -700,15 +704,7 @@ void FourierTransform::bwd(complex<double>* val)
// copy from rbuf to val
// scatter index array iunpack
{
const int len = np012loc() * 2;
double* const pv = (double*) &val[0];
#pragma omp parallel for
for ( int i = 0; i < len; i++ )
{
pv[i] = 0.0;
}
}
memset((void*)&val[0],0,2*np012loc()*sizeof(double));
#if TIMING
tm_b_zero.stop();
......@@ -1594,15 +1590,10 @@ void FourierTransform::init_lib(void)
////////////////////////////////////////////////////////////////////////////////
void FourierTransform::vector_to_zvec(const complex<double> *c)
{
// map one real or complex function to zvec
memset((void*)&zvec_[0],0,zvec_.size()*sizeof(complex<double>));
const int ng = basis_.localsize();
const int zvec_size = zvec_.size();
double* const pz = (double*) &zvec_[0];
#pragma omp parallel for
for ( int i = 0; i < zvec_size; i++ )
{
pz[2*i] = 0.0;
pz[2*i+1] = 0.0;
}
const double* const pc = (double*) &c[0];
if ( basis_.real() )
{
......@@ -1655,16 +1646,10 @@ void FourierTransform::zvec_to_vector(complex<double> *c)
void FourierTransform::doublevector_to_zvec(const complex<double> *c1,
const complex<double> *c2)
{
// Mapping of two real functions onto zvec
// map two real functions to zvec
assert(basis_.real());
const int zvec_size = zvec_.size();
memset((void*)&zvec_[0],0,zvec_.size()*sizeof(complex<double>));
double* const pz = (double*) &zvec_[0];
#pragma omp parallel for
for ( int i = 0; i < zvec_size; i++ )
{
pz[2*i] = 0.0;
pz[2*i+1] = 0.0;
}
const int ng = basis_.localsize();
const double* const pc1 = (double*) &c1[0];
const double* const pc2 = (double*) &c2[0];
......
......@@ -65,6 +65,7 @@ class FourierTransform
int ntrans0_,ntrans1_,ntrans2_;
int nvec_;
bool basis_fits_in_grid_;
std::vector<int> np2_loc_; // np2_loc_[iproc], iproc=0, nprocs_-1
std::vector<int> np2_first_; // np2_first_[iproc], iproc=0, nprocs_-1
......
......@@ -203,13 +203,13 @@ void approximateIntegral(const double omega_kF, const double Hs2,
for ( int i = 1; i < no_integral / 2; i++ )
{
term1 = term1 / bw2_D_term * static_cast<double> (i);
term1 = i * term1 / bw2_D_term;
factor2 = -factor2 * r9_4A;
term2 = term2 + arg_n;
integral[2 * ( i + 1 )] = term1 + factor2 * term2;
arg_n = -arg_n * static_cast<double> (i) / arg;
arg_n = -arg_n * i / arg;
}
// The n = 1 integral is
......@@ -537,9 +537,8 @@ void HSE_enhance(const double s_inp, const double kF, const double w,
{
// update values
prefactor *= static_cast<double> (i) * r1_D_term;
summand *= static_cast<double> (2 * i - 1) / static_cast<double> (2 * i)
* q_q_1;
prefactor *= i * r1_D_term;
summand *= (2.0 * i - 1.0) / (2.0 * i) * q_q_1;
sum -= summand;
intYExpErfc.push_back(prefactor * sum);
......@@ -682,7 +681,7 @@ void gcor2(double a, double a1, double b1, double b2, double b3, double b4,
////////////////////////////////////////////////////////////////////////////////
//
// calculate correletation energy of PBE functional
// calculate correlation energy of the PBE functional
// K.Burke's modification of PW91 codes, May 14, 1996.
// Modified again by K.Burke, June 29, 1996, with simpler Fx(s)
// Translated into C and modified by F.Gygi, Dec 9, 1996.
......@@ -900,7 +899,6 @@ void HSEFunctional::setxc(void)
if ( _np == 0 ) return;
if ( _nspin == 1 )
{
// test for void pointer
assert(rho != 0);
assert(grad_rho[0] != 0 && grad_rho[1] != 0 && grad_rho[2] != 0);
assert(exc != 0);
......@@ -928,7 +926,6 @@ void HSEFunctional::setxc(void)
}
else
{
// test for void pointer
assert(rho_up != 0);
assert(rho_dn != 0);
assert(grad_rho_up[0] != 0 && grad_rho_up[1] != 0 && grad_rho_up[2] != 0);
......
......@@ -99,8 +99,19 @@ CXXFLAGS += -DTARGET='"$(TARGET)"'
testEnergyFunctional: testEnergyFunctional.o EnergyFunctional.o Basis.o \
SlaterDet.o Matrix.o UnitCell.o Context.o FourierTransform.o \
Wavefunction.o Species.o Atom.o AtomSet.o StructureFactor.o \
ChargeDensity.o \
sinft.o spline.o
ExtForceSet.o ExtForce.o PairExtForce.o AtomicExtForce.o \
GlobalExtForce.o \
ConstraintSet.o Constraint.o DistanceConstraint.o \
AngleConstraint.o TorsionConstraint.o PositionConstraint.o \
NonLocalPotential.o sampling.o Base64Transcoder.o \
ChargeDensity.o sinft.o spline.o \
XCOperator.o ExchangeOperator.o Bisection.o \
XCPotential.o LDAFunctional.o VWNFunctional.o \
PBEFunctional.o BLYPFunctional.o B3LYPFunctional.o \
BHandHLYPFunctional.o \
ExponentialIntegral.o HSEFunctional.o RSHFunctional.o \
ConfinementPotential.o ElectricEnthalpy.o MLWFTransform.o \
jade.o BasisMapping.o
$(LD) $(DFLAGS) -o $@ $^ $(LDFLAGS)
testSlaterDet: testSlaterDet.o SlaterDet.o FourierTransform.o \
Basis.o UnitCell.o Matrix.o Context.o Base64Transcoder.o
......@@ -113,6 +124,10 @@ CXXFLAGS += -DTARGET='"$(TARGET)"'
$(LD) $(DFLAGS) -o $@ $^ $(LDFLAGS)
testXCFunctional: testXCFunctional.o LDAFunctional.o PBEFunctional.o
$(LD) $(DFLAGS) -o $@ $^ $(LDFLAGS)
testLDAFunctional: testLDAFunctional.o LDAFunctional.o
$(LD) $(DFLAGS) -o $@ $^ $(LDFLAGS)
testVWN: testVWN.o VWNFunctional.o LDAFunctional.o
$(LD) $(DFLAGS) -o $@ $^ $(LDFLAGS)
testMatrix: testMatrix.o Matrix.o Context.o
$(LD) $(DFLAGS) -o $@ $^ $(LDFLAGS)
testjacobi: testjacobi.o jacobi.o Matrix.o Context.o
......@@ -131,6 +146,8 @@ CXXFLAGS += -DTARGET='"$(TARGET)"'
testXMLGFPreprocessor: testXMLGFPreprocessor.o XMLGFPreprocessor.o Context.o \
Base64Transcoder.o Matrix.o
$(LD) $(DFLAGS) -o $@ $^ $(LDFLAGS)
testBase64Transcoder: testBase64Transcoder.o Base64Transcoder.o
$(LD) $(DFLAGS) -o $@ $^ $(LDFLAGS)
test_fftw: test_fftw.o $(PLTOBJECTS)
$(LD) $(DFLAGS) -o $@ $^ $(LDFLAGS)
#------------------------------------------------------------------------------
......@@ -197,6 +214,8 @@ AtomsDyn.o: Control.h
B3LYPFunctional.o: B3LYPFunctional.h XCFunctional.h BLYPFunctional.h
B3LYPFunctional.o: VWNFunctional.h
B3LYPFunctional.o: XCFunctional.h
BHandHLYPFunctional.o: BHandHLYPFunctional.h XCFunctional.h BLYPFunctional.h
BHandHLYPFunctional.o: XCFunctional.h
BLYPFunctional.o: BLYPFunctional.h XCFunctional.h
BLYPFunctional.o: XCFunctional.h
BMDIonicStepper.o: BMDIonicStepper.h IonicStepper.h Sample.h AtomSet.h
......@@ -664,7 +683,8 @@ SpeciesReader.o: Species.h SpeciesReader.h StructuredDocumentHandler.h StrX.h
SpeciesReader.o: StructureHandler.h SpeciesHandler.h
StatusCmd.o: UserInterface.h Sample.h AtomSet.h Context.h blacs.h Atom.h
StatusCmd.o: D3vector.h UnitCell.h D3tensor.h blas.h ConstraintSet.h
StatusCmd.o: ExtForceSet.h Wavefunction.h Control.h
StatusCmd.o: ExtForceSet.h Wavefunction.h Control.h ChargeDensity.h Timer.h
StatusCmd.o: FourierTransform.h
StrainCmd.o: UserInterface.h Sample.h AtomSet.h Context.h blacs.h Atom.h
StrainCmd.o: D3vector.h UnitCell.h D3tensor.h blas.h ConstraintSet.h
StrainCmd.o: ExtForceSet.h Wavefunction.h Control.h
......@@ -729,8 +749,9 @@ XCOperator.o: Wavefunction.h Control.h
XCPotential.o: XCPotential.h Control.h D3vector.h ChargeDensity.h Timer.h
XCPotential.o: Context.h blacs.h LDAFunctional.h XCFunctional.h
XCPotential.o: VWNFunctional.h PBEFunctional.h BLYPFunctional.h
XCPotential.o: HSEFunctional.h RSHFunctional.h B3LYPFunctional.h Basis.h
XCPotential.o: UnitCell.h FourierTransform.h blas.h
XCPotential.o: HSEFunctional.h RSHFunctional.h B3LYPFunctional.h
XCPotential.o: BHandHLYPFunctional.h Basis.h UnitCell.h FourierTransform.h
XCPotential.o: blas.h
XCPotential.o: Control.h D3vector.h ChargeDensity.h Timer.h Context.h blacs.h
XMLGFPreprocessor.o: Timer.h Context.h blacs.h Base64Transcoder.h Matrix.h
XMLGFPreprocessor.o: XMLGFPreprocessor.h
......@@ -749,16 +770,17 @@ qb.o: Control.h Timer.h AngleCmd.h AtomCmd.h ComputeMLWFCmd.h MLWFTransform.h
qb.o: BasisMapping.h ConstraintCmd.h DistanceCmd.h ExtForceCmd.h
qb.o: FoldInWsCmd.h HelpCmd.h KpointCmd.h ListAtomsCmd.h ListSpeciesCmd.h
qb.o: LoadCmd.h MoveCmd.h PartialChargeCmd.h PlotCmd.h PrintCmd.h QuitCmd.h
qb.o: RandomizeRCmd.h RandomizeVCmd.h RandomizeWfCmd.h ResetVcmCmd.h
qb.o: RescaleVCmd.h RseedCmd.h RunCmd.h SaveCmd.h SetCmd.h SetVelocityCmd.h
qb.o: SpeciesCmd.h StatusCmd.h StrainCmd.h TorsionCmd.h BisectionCmd.h
qb.o: Bisection.h SlaterDet.h Basis.h Matrix.h AlphaPBE0.h AlphaRSH.h
qb.o: AtomsDyn.h BetaRSH.h BlHF.h BtHF.h Cell.h CellDyn.h CellLock.h
qb.o: CellMass.h ChargeMixCoeff.h ChargeMixNdim.h ChargeMixRcut.h Debug.h
qb.o: Dspin.h Ecut.h Ecutprec.h Ecuts.h Efield.h Polarization.h Emass.h
qb.o: ExtStress.h FermiTemp.h IterCmd.h IterCmdPeriod.h Dt.h MuRSH.h Nempty.h
qb.o: NetCharge.h Nrowmax.h Nspin.h RefCell.h ScfTol.h Stress.h Thermostat.h
qb.o: ThTemp.h ThTime.h ThWidth.h WfDiag.h WfDyn.h Xc.h
qb.o: RandomizeRCmd.h RandomizeVCmd.h RandomizeWfCmd.h ResetRotationCmd.h
qb.o: ResetVcmCmd.h RescaleVCmd.h RseedCmd.h RunCmd.h SaveCmd.h SetCmd.h
qb.o: SetVelocityCmd.h SpeciesCmd.h StatusCmd.h ChargeDensity.h
qb.o: FourierTransform.h StrainCmd.h TorsionCmd.h BisectionCmd.h Bisection.h
qb.o: SlaterDet.h Basis.h Matrix.h AlphaPBE0.h AlphaRSH.h AtomsDyn.h
qb.o: BetaRSH.h BlHF.h BtHF.h Cell.h CellDyn.h CellLock.h CellMass.h
qb.o: ChargeMixCoeff.h ChargeMixNdim.h ChargeMixRcut.h Debug.h Dspin.h Ecut.h
qb.o: Ecutprec.h Ecuts.h Efield.h Polarization.h Emass.h ExtStress.h
qb.o: FermiTemp.h IterCmd.h IterCmdPeriod.h Dt.h MuRSH.h Nempty.h NetCharge.h
qb.o: Nrowmax.h Nspin.h RefCell.h ScfTol.h Stress.h Thermostat.h ThTemp.h
qb.o: ThTime.h ThWidth.h WfDiag.h WfDyn.h Xc.h
qbox_xmlns.o: qbox_xmlns.h
release.o: release.h
sinft.o: sinft.h
......@@ -776,9 +798,9 @@ testContext.o: Context.h blacs.h
testEnergyFunctional.o: Context.h blacs.h Sample.h AtomSet.h Atom.h
testEnergyFunctional.o: D3vector.h UnitCell.h D3tensor.h blas.h
testEnergyFunctional.o: ConstraintSet.h ExtForceSet.h Wavefunction.h
testEnergyFunctional.o: Control.h EnergyFunctional.h StructureFactor.h
testEnergyFunctional.o: ElectricEnthalpy.h Matrix.h SlaterDet.h Basis.h
testEnergyFunctional.o: Timer.h
testEnergyFunctional.o: Control.h ChargeDensity.h Timer.h EnergyFunctional.h
testEnergyFunctional.o: StructureFactor.h ElectricEnthalpy.h Matrix.h
testEnergyFunctional.o: SlaterDet.h Basis.h
testFourierTransform.o: Basis.h D3vector.h UnitCell.h FourierTransform.h
testFourierTransform.o: Timer.h
testLDAFunctional.o: LDAFunctional.h XCFunctional.h
......
......@@ -212,13 +212,13 @@ void RSHFunctional::approximateIntegral(const double omega_kF, const double Hs2,
for ( int i = 1; i < no_integral / 2; i++ )
{
term1 = term1 / bw2_D_term * static_cast<double> (i);
term1 = i * term1 / bw2_D_term;
factor2 = -factor2 * r9_4A;
term2 = term2 + arg_n;
integral[2 * ( i + 1 )] = term1 + factor2 * term2;
arg_n = -arg_n * static_cast<double> (i) / arg;
arg_n = -arg_n * i / arg;
}
// The n = 1 integral is
......@@ -546,9 +546,8 @@ void RSHFunctional::RSH_enhance(const double s_inp, const double kF,
{
// update values
prefactor *= static_cast<double> (i) * r1_D_term;
summand *= static_cast<double> (2 * i - 1) / static_cast<double> (2 * i)
* q_q_1;
prefactor *= i * r1_D_term;
summand *= (2.0 * i - 1.0) / (2.0 * i) * q_q_1;
sum -= summand;
intYExpErfc.push_back(prefactor * sum);
......@@ -691,7 +690,7 @@ void RSHFunctional::gcor2(double a, double a1, double b1, double b2,
////////////////////////////////////////////////////////////////////////////////
//
// calculate correletation energy of PBE functional
// calculate correlation energy of the PBE functional
// K.Burke's modification of PW91 codes, May 14, 1996.
// Modified again by K.Burke, June 29, 1996, with simpler Fx(s)
// Translated into C and modified by F.Gygi, Dec 9, 1996.
......@@ -909,7 +908,6 @@ void RSHFunctional::setxc(void)
if ( _np == 0 ) return;
if ( _nspin == 1 )
{
// test for void pointer
assert(rho != 0);
assert(grad_rho[0] != 0 && grad_rho[1] != 0 && grad_rho[2] != 0);
assert(exc != 0);
......@@ -937,7 +935,6 @@ void RSHFunctional::setxc(void)
}
else
{
// test for void pointer
assert(rho_up != 0);
assert(rho_dn != 0);
assert(grad_rho_up[0] != 0 && grad_rho_up[1] != 0 && grad_rho_up[2] != 0);
......
......@@ -22,6 +22,8 @@
#include <iostream>
#include "UserInterface.h"
#include "Sample.h"
#include "ChargeDensity.h"
#include "FourierTransform.h"
class StatusCmd : public Cmd
{
......@@ -47,8 +49,22 @@ class StatusCmd : public Cmd
int action(int argc, char **argv)
{
// compute the size of the potential grid
int np0v = 0;
int np1v = 0;
int np2v = 0;
if ( s->wf.ecut() > 0 && s->wf.cell().volume() > 0 )
{
ChargeDensity cd(s->wf);
np0v = cd.vft()->np0();
np1v = cd.vft()->np1();
np2v = cd.vft()->np2();
}
if ( ui->onpe0() )
{
cout << "<np0v> " << np0v << " </np0v> "
<< "<np1v> " << np1v << " </np1v> "
<< "<np2v> " << np2v << " </np2v>" << endl;
s->wf.info(cout,"wf");
if ( s->wfv != 0 )
s->wfv->info(cout,"wfv");
......
......@@ -15,17 +15,15 @@
// UserInterface.h:
//
////////////////////////////////////////////////////////////////////////////////
// $ Id: $
#ifndef USER_INTERFACE_H
#define USER_INTERFACE_H
#include <iostream>
#include <string>
#include <cstring>
#include <iomanip>
#include <cstring> // strncpy(), strtok()
#include <cstdlib> // free(), system()
#include <list>
#include <algorithm>
class UserInterface;
......@@ -53,7 +51,8 @@ class UserInterface
{
private:
int readCmd(char *s, int max, std::istream &fp, bool echo);
int readCmd(std::string& s, std::istream &is, bool echo);
void execCmd(std::string s, std::string prompt);
bool terminate_;
bool onpe0_;
......@@ -108,9 +107,10 @@ class UserInterface
}
};
void processCmds(std::istream &cmdstream, const char *prompt, bool echo);
void processCmds(std::istream &cmdstream, const std::string prompt,
bool echo);
void processCmdsServer(std::string inputfilename, std::string outputfilename,
const char *prompt, bool echo);
const std::string prompt, bool echo);
void terminate(void) { terminate_ = true; }
......
......@@ -24,9 +24,6 @@ using namespace std;
#include <unistd.h>
#include <cstdlib>
#include <fstream>
#if AIX
#include<filehdr.h>
#endif
#include "isodate.h"
#include "release.h"
......@@ -119,10 +116,6 @@ using namespace std;
#include "WfDyn.h"
#include "Xc.h"
#if BGLDEBUG
#include <rts.h>
#endif
int main(int argc, char **argv, char **envp)
{
Timer tm;
......@@ -132,21 +125,6 @@ int main(int argc, char **argv, char **envp)
MPI_Init(&argc,&argv);
#endif
#if BGLDEBUG
{
int myrank,mysize;
BGLPersonality personality;
rts_get_personality (&personality, sizeof(personality));
MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
MPI_Comm_size(MPI_COMM_WORLD, &mysize);
cout << myrank << ": at "
<< personality.xCoord << " "
<< personality.yCoord << " "
<< personality.zCoord << endl;
}
#endif
{
Context ctxt(MPI_COMM_WORLD);
......@@ -176,25 +154,10 @@ int main(int argc, char **argv, char **envp)
cout << " ============================\n\n";
cout << "\n";
cout << "<release> " << release() << " " << TARGET << " </release>" << endl;
#ifdef SVN_VERSION
cout << "<svn_version> " << SVN_VERSION << " </svn_version>" << endl;
#endif
// Identify executable name, checksum, size and link date
if ( getlogin() != 0 )
cout << "<user> " << getlogin() << " </user>" << endl;
#if AIX || OSF1
// read filehdr for link time
filehdr hdr;
FILE *fx = fopen(argv[0],"r");
if ( fx != 0 )
{
size_t sz = fread((void*)&hdr,sizeof(filehdr),1,fx);
fclose(fx);
string s = ctime((time_t*)&hdr.f_timdat);
cout << "<linktime> " << s << " </linktime>" << endl;
}
#endif
// Identify platform
{
......@@ -346,15 +309,26 @@ int main(int argc, char **argv, char **envp)
string inputfilename(argv[1]);
string outputfilename("stdout");
ifstream in;
int file_ok = 0;
if ( ctxt.onpe0() )
{
in.open(argv[1],ios::in);
if ( in )
if ( in )
{
// file was opened on process 0
file_ok = 1;
}
}
MPI_Bcast(&file_ok,1,MPI_INT,0,MPI_COMM_WORLD);
if ( file_ok )
{
ui.processCmds(in, "[qbox]", echo);
}
else
{
cout << " qbox: could not open input file "
<< argv[1] << endl;
ctxt.abort(1);
if ( ctxt.onpe0() )
cout << " Could not open input file " << argv[1] << endl;
}
}
else if ( argc == 4 )
......
......@@ -18,6 +18,7 @@
#include "Base64Transcoder.h"
#include <iostream>
#include <cassert>
using namespace std;
int main()
......
......@@ -26,6 +26,7 @@
#include <fstream>
#include <iomanip>
#include <cassert>
#include <cstdlib>
using namespace std;
#ifdef USE_MPI
......@@ -39,7 +40,6 @@ int main(int argc, char **argv)
#endif
{
// use: testBasisMapping a0 a1 a2 b0 b1 b2 c0 c1 c2 ecut npr npc
double err;
assert(argc==13);
D3vector a(atof(argv[1]),atof(argv[2]),atof(argv[3]));
D3vector b(atof(argv[4]),atof(argv[5]),atof(argv[6]));
......
......@@ -40,7 +40,13 @@ int main(int argc, char **argv)
{
// use:
// testChargeDensity a0 a1 a2 b0 b1 b2 c0 c1 c2 ecut nel nspin nkp
assert(argc==14);
if ( argc !=14 )
{
cout <<
" use: testChargeDensity a0 a1 a2 b0 b1 b2 c0 c1 c2 ecut nel nspin nkp"
<< endl;
return 1;
}
D3vector a(atof(argv[1]),atof(argv[2]),atof(argv[3]));
D3vector b(atof(argv[4]),atof(argv[5]),atof(argv[6]));
D3vector c(atof(argv[7]),atof(argv[8]),atof(argv[9]));
......
......@@ -32,6 +32,11 @@ int main(int argc, char **argv)
MPI_Comm_size(MPI_COMM_WORLD, &npes);
MPI_Comm_rank(MPI_COMM_WORLD, &mype);
if ( argc != 3 )
{
cout << "use: testContext nrow ncol" << endl;
return 1;
}
int nr = atoi(argv[1]);
int nc = atoi(argv[2]);
......
......@@ -19,6 +19,7 @@
#include "Context.h"
#include "Sample.h"
#include "Wavefunction.h"
#include "ChargeDensity.h"
#include "EnergyFunctional.h"
#include "Timer.h"
......@@ -38,7 +39,12 @@ int main(int argc, char **argv)
#endif
{
// use: testEnergyFunctional a0 a1 a2 b0 b1 b2 c0 c1 c2 ecut nel
assert(argc==12);
if ( argc != 12 )
{
cout << "use: testEnergyFunctional a0 a1 a2 b0 b1 b2 c0 c1 c2 ecut nel"
<< endl;
return 1;
}
D3vector a(atof(argv[1]),atof(argv[2]),atof(argv[3]));
D3vector b(atof(argv[4]),atof(argv[5]),atof(argv[6]));
D3vector c(atof(argv[7]),atof(argv[8]),atof(argv[9]));
......@@ -48,21 +54,17 @@ int main(int argc, char **argv)
Timer tm;
Context ctxt;
Context ctxt(MPI_COMM_WORLD);
cout << " initial context: " << ctxt;
Sample s(ctxt);
s.atoms.addAtom(
new Atom("G","gaussium",D3vector(0.2,0.3,0),D3vector(0,0,0)));
s.atoms.listAtoms();
s.atoms.listSpecies();
s.wf.resize(cell,cell,ecut);
s.wf.set_nel(nel);
if ( ctxt.onpe0() ) cout << " nel: " << s.wf.nel() << endl;
s.wf.update_occ();
//s.wf.randomize(0.05);
s.wf.update_occ(0.0);
s.wf.randomize(0.05);
tm.reset();
tm.start();
......@@ -70,16 +72,31 @@ int main(int argc, char **argv)
tm.stop();
cout << " Gram: CPU/Real: " << tm.cpu() << " / " << tm.real() << endl;
ChargeDensity cd(s.wf);
tm.reset();
tm.start();
EnergyFunctional ef(s);
cout << " ChargeDensity::update_density..." << endl;
cd.update_density();
tm.stop();
cout << " ChargeDensity::update_density: CPU/Real: "
<< tm.cpu() << " / " << tm.real() << endl;
s.ctrl.xc = "LDA";
s.ctrl.polarization = "OFF";
tm.reset();
tm.start();
EnergyFunctional ef(s,cd);
tm.stop();
cout << " EnergyFunctional:ctor: CPU/Real: "
<< tm.cpu() << " / " << tm.real() << endl;
tm.reset();
tm.start();
cout << " ef.energy(): " << ef.energy() << endl;
Wavefunction dwf(s.wf);
vector<vector<double> > fion;
valarray<double> sigma(6);
double e = ef.energy(true,dwf,false,fion,false,sigma);
cout << " ef.energy(): " << e << endl;
tm.stop();
cout << " EnergyFunctional:energy: CPU/Real: "
<< tm.cpu() << " / " << tm.real() << endl;
......
......@@ -23,27 +23,26 @@
// dExc/da must be 0.911682
#include<iostream>
#include<vector>
#include "LDAFunctional.h"
#include <cassert>
#include <cmath>
#include <cstdlib>
using namespace std;
int main(int argc, char **argv)
{
// use: testxcf alat np
if ( argc != 3 )
{
cout << " use: testLDAFunctional alat np" << endl;
return 0;
cout << " use: testLDAFunctional alat n" << endl;
return 1;
}
assert(argc==3);
double a = atof(argv[1]);
double omega = a*a*a;
int n = atoi(argv[2]);
int n3 = n*n*n;
double *rh = new double[n3];
double *vxc = new double[n3];
double *exc = new double[n3];
vector<vector<double> > rh(1);
rh[0].resize(n3);
double excsum = 0.0, dxcsum = 0.0;
double rc = 0.1 * a;
......@@ -62,8 +61,8 @@ int main(int argc, char **argv)
double z = ( k * a ) / n - a/2;
double r2 = x*x + y*y + z*z;
int ii = i + n * ( j + n * k );
rh[ii] = fac * exp( -r2 / (rc*rc) );
sum += rh[ii];
rh[0][ii] = fac * exp( -r2 / (rc*rc) );
sum += rh[0][ii];
}