Skip to content

qbox-public
qbox-public
Commit e39dbb70 by mahe
Options

merged origin/develop, to be tested

parents 0d62c878 e9d66bfe
Showing with 1107 additions and 392 deletions
build/centos6.mk
......@@ -60,11 +60,11 @@ ifeq ($(FFT),FFTW3)
PLTFLAGS += -DFFTWMEASURE
#PLTFLAGS += -DFFTW_TRANSPOSE
PLTFLAGS += -DFFTW3_2D
FFTWDIR=$(HOME)/software/fftw/fftw-3.3.4
FFTWINCLUDEDIR=$(FFTWDIR)/api
FFTWLIBDIR=$(FFTWDIR)/.libs
INCLUDE += -I$(FFTWINCLUDEDIR)
LIBPATH += -L$(FFTWLIBDIR)
#FFTWDIR=$(HOME)/software/fftw/fftw-3.3.4
#FFTWINCLUDEDIR=$(FFTWDIR)/api
#FFTWLIBDIR=$(FFTWDIR)/.libs
#INCLUDE += -I$(FFTWINCLUDEDIR)
#LIBPATH += -L$(FFTWLIBDIR)
LIBS += -lfftw3
endif
......
src/AngleConstraint.C
......@@ -407,9 +407,10 @@ ostream& AngleConstraint::print( ostream &os )
os << name2_ << " " << name3_ << "\"\n";
os.setf(ios::fixed,ios::floatfield);
os.setf(ios::right,ios::adjustfield);
os << " value=\"" << setprecision(6) << angle_;
os << "\" velocity=\"" << setprecision(6) << velocity_ << "\"\n";
os << " force=\"" << setprecision(6) << force_;
os << "\" weight=\"" << setprecision(6) << weight_ << "\"/>";
os << " velocity=\"" << setprecision(6) << velocity_ << "\"";
os << " weight=\"" << setprecision(6) << weight_ << "\">\n";
os << " <value> " << setprecision(6) << angle_ << " </value>";
os << " <force> " << setprecision(6) << force_ << " </force>\n";
os << " </constraint>";
return os;
}
src/BOSampleStepper.C
......@@ -347,6 +347,9 @@ void BOSampleStepper::step(int niter)
wkerker[i] = 1.0;
}
if ( onpe0 )
cout << "<net_charge> " << atoms.nel()-wf.nel() << " </net_charge>\n";
// Next line: special case of niter=0: compute GS only
for ( int iter = 0; iter < max(niter,1); iter++ )
{
......@@ -377,6 +380,7 @@ void BOSampleStepper::step(int niter)
if ( onpe0 )
{
cout << cd_;
cout << ef_;
if ( ef_.el_enth() )
cout << *ef_.el_enth();
......@@ -758,6 +762,9 @@ void BOSampleStepper::step(int niter)
cd_.update_density();
tmap["charge"].stop();
if ( onpe0 )
cout << cd_;
// charge mixing
if ( nite_ > 0 )
{
......@@ -1191,10 +1198,10 @@ void BOSampleStepper::step(int niter)
tmap["energy"].stop();
if ( onpe0 )
{
cout << cd_;
cout << ef_;
if ( ef_.el_enth() )
cout << *ef_.el_enth();
}
}
......
src/ChargeDensity.C
......@@ -55,6 +55,7 @@ wf_(wf), vcomm_(wf.sd(0,0)->basis().comm())
#endif
vft_ = new FourierTransform(*vbasis_,np0v,np1v,np2v);
total_charge_.resize(wf.nspin());
rhor.resize(wf.nspin());
rhog.resize(wf.nspin());
for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
......@@ -162,13 +163,7 @@ void ChargeDensity::update_density(void)
// 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;
}
total_charge_[ispin] = sum;
tmap["charge_vft"].start();
vft_->forward(&rhotmp[0],&rhog[ispin][0]);
......@@ -272,4 +267,4 @@ void ChargeDensity::update_rhog(void)
vft_->forward(&rhotmp[0],&rhog[ispin][0]);
}
}
}
\ No newline at end of file
src/ChargeDensity.h
......@@ -44,6 +44,7 @@ class ChargeDensity
FourierTransform* vft_;
std::vector<FourierTransform*> ft_; // ft_[ikp];
std::valarray<std::complex<double> > rhotmp;
std::vector<double> total_charge_;
public:
......@@ -62,8 +63,13 @@ class ChargeDensity
Basis* vbasis(void) const { return vbasis_; }
FourierTransform* vft(void) const { return vft_; }
FourierTransform* ft(int ikp) const { return ft_[ikp]; }
double total_charge(int ispin) const { return total_charge_[ispin]; }
double total_charge(void) const;
void print(std::ostream& os) const;
ChargeDensity(const Wavefunction& wf);
~ChargeDensity();
};
std::ostream& operator << ( std::ostream& os, const ChargeDensity& cd );
#endif
src/ConstraintCmd.h
......@@ -45,6 +45,7 @@ class ConstraintCmd : public Cmd
" constraint list\n"
" constraint enforce\n\n"
" Constraints are enforced at each MD step if ions are allowed to move.\n"
" If a distance or angle is replaced by '*', the current value is used.\n"
" Velocity parameters are optional.\n\n";
}
......
src/ConstraintSet.C
......@@ -23,6 +23,7 @@
#include "TorsionConstraint.h"
#include "Atom.h"
#include "AtomSet.h"
#include <cstring>
#include "Context.h"
#include <iostream>
#include <iomanip>
......@@ -192,7 +193,18 @@ bool ConstraintSet::define_constraint(AtomSet &atoms, int argc, char **argv)
return false;
}
distance = atof(argv[6]);
// define distance value
if ( !strcmp(argv[6],"*") )
{
// use current distance
distance = length(a1->position()-a2->position());
if ( onpe0 )
cout << "ConstraintSet::define_constraint: using current distance "
<< distance << endl;
}
else
distance = atof(argv[6]);
if ( argc == 8 )
{
velocity = atof(argv[7]);
......@@ -223,8 +235,10 @@ bool ConstraintSet::define_constraint(AtomSet &atoms, int argc, char **argv)
if ( found )
{
if ( onpe0 )
cout << " ConstraintSet: constraint is already defined:\n"
<< " cannot define constraint" << endl;
cout << " ConstraintSet: a distance constraint named " << name << endl
<< " or involving atoms " << name1 << " " << name2 << endl
<< " is already defined" << endl
<< " ConstraintSet: cannot define constraint" << endl;
return false;
}
else
......@@ -280,7 +294,33 @@ bool ConstraintSet::define_constraint(AtomSet &atoms, int argc, char **argv)
return false;
}
const double angle = atof(argv[7]);
// define angle value
double angle = 0.0;
if ( !strcmp(argv[7],"*") )
{
// use current angle
D3vector r12(a1->position()-a2->position());
D3vector r32(a3->position()-a2->position());
if ( norm2(r12) == 0.0 || norm2(r32) == 0.0 )
{
if ( onpe0 )
{
cout << " ConstraintSet: cannot define angle constraint." << endl;
cout << " ConstraintSet: atoms are too close" << endl;
return false;
}
}
const double sp = normalized(r12) * normalized(r32);
const double c = max(-1.0,min(1.0,sp));
angle = (180.0/M_PI)*acos(c);
if ( onpe0 )
cout << "ConstraintSet::define_constraint: using current angle "
<< angle << endl;
}
else
angle = atof(argv[7]);
double velocity = 0.0;
if ( argc == 9 )
{
......@@ -317,9 +357,10 @@ bool ConstraintSet::define_constraint(AtomSet &atoms, int argc, char **argv)
if ( found )
{
if ( onpe0 )
cout << " ConstraintSet:set_constraint: an angle constraint "
<< name1 << " " << name2 << " " << name3
<< " was found" << endl
cout << " ConstraintSet: an angle constraint named " << name << endl
<< " or involving atoms "
<< name1 << " " << name2 << " " << name3 << endl
<< " is already defined" << endl
<< " ConstraintSet: cannot define constraint" << endl;
return false;
}
......@@ -380,7 +421,48 @@ bool ConstraintSet::define_constraint(AtomSet &atoms, int argc, char **argv)
return false;
}
double angle = atof(argv[8]);
// define angle value
double angle = 0.0;
if ( !strcmp(argv[8],"*") )
{
// use current angle
D3vector r12(a1->position()-a2->position());
D3vector r32(a3->position()-a2->position());
D3vector r43(a4->position()-a3->position());
if ( norm2(r12) == 0.0 || norm2(r32) == 0.0 || norm2(r43) == 0.0 )
{
if ( onpe0 )
{
cout << " ConstraintSet: cannot define torsion constraint." << endl;
cout << " ConstraintSet: atoms are too close" << endl;
return false;
}
}
D3vector e12(normalized(r12));
D3vector e32(normalized(r32));
D3vector e23(-e32);
D3vector e43(normalized(r43));
const double sin123 = length(e12^e32);
const double sin234 = length(e23^e43);
if ( sin123 != 0.0 && sin234 != 0.0 )
{
D3vector e123 = normalized(e12^e32);
D3vector e234 = normalized(e23^e43);
double cc = max(min(e123*e234,1.0),-1.0);
double ss = max(min((e123^e234)*e32,1.0),-1.0);
angle = (180.0/M_PI) * atan2(ss,cc);
}
if ( onpe0 )
cout << "ConstraintSet::define_constraint: using current angle "
<< angle << endl;
}
else
angle = atof(argv[8]);
if ( angle > 180.0 )
while ( angle > 180.0 ) angle -= 360.0;
else if ( angle < -180.0 )
......@@ -401,7 +483,7 @@ bool ConstraintSet::define_constraint(AtomSet &atoms, int argc, char **argv)
assert(pc != 0);
// check if an equivalent constraint (name1,name2,name3,name4) or
// (name4,name3,name2,name1) is defined
if ( pc->type() == "angle" )
if ( pc->type() == "torsion" )
found = ( pc->names(0) == name1 &&
pc->names(1) == name2 &&
pc->names(2) == name3 &&
......@@ -416,8 +498,9 @@ bool ConstraintSet::define_constraint(AtomSet &atoms, int argc, char **argv)
if ( found )
{
if ( onpe0 )
cout << " ConstraintSet: a torsion constraint "
<< name1 << " " << name2 << " " << name3 << " " << name4
cout << " ConstraintSet: a torsion constraint named " << name << endl
<< " or involving atoms "
<< name1 << " " << name2 << " " << name3 << " " << name4 << endl
<< " is already defined" << endl
<< " ConstraintSet: cannot define constraint" << endl;
return false;
......
src/DistanceConstraint.C
......@@ -228,11 +228,10 @@ ostream& DistanceConstraint::print( ostream &os )
os << "\" atoms=\"" << name1_ << " " << name2_ << "\"\n";
os.setf(ios::fixed,ios::floatfield);
os.setf(ios::right,ios::adjustfield);
os << " value=\"" << setprecision(6) << distance_;
os << "\" velocity=\"" << setprecision(6) << velocity_ << "\"\n";
os << " force=\"" << setprecision(6) << force_;
os << "\" weight=\"" << setprecision(6) << weight_ << "\"/>";
os << " velocity=\"" << setprecision(6) << velocity_ << "\"";
os << " weight=\"" << setprecision(6) << weight_ << "\">\n";
os << " <value> " << setprecision(6) << distance_ << " </value>";
os << " <force> " << setprecision(6) << force_ << " </force>\n";
os << " </constraint>";
return os;
}
src/ExchangeOperator.C
......@@ -901,7 +901,7 @@ double ExchangeOperator::compute_exchange_at_gamma_(const Wavefunction &wf,
<< maxsweep << endl;
assert(maxsweep >= 0);
}
double tol = 1.0;
double tol = 0.001;
if ( s_.ctrl.debug.find("BISECTION_TOL") != string::npos )
{
// override tolerance for bisection
......
src/KpointCmd.h
......@@ -39,6 +39,7 @@ class KpointCmd : public Cmd
" syntax:\n\n"
" kpoint add kx ky kz weight\n"
" kpoint delete kx ky kz \n"
" kpoint move kx ky kz new_kx new_ky new_kz\n"
" kpoint list\n\n";
}
......@@ -79,6 +80,22 @@ class KpointCmd : public Cmd
double kz = atof(argv[4]);
s->wf.del_kpoint(D3vector(kx,ky,kz));
}
else if ( subcmd == "move" )
{
if ( argc != 8 )
{
if ( onpe0 )
cout << help_msg();
return 1;
}
double kx = atof(argv[2]);
double ky = atof(argv[3]);
double kz = atof(argv[4]);
double newkx = atof(argv[5]);
double newky = atof(argv[6]);
double newkz = atof(argv[7]);
s->wf.move_kpoint(D3vector(kx,ky,kz),D3vector(newkx,newky,newkz));
}
else if ( subcmd == "list" )
{
if ( argc != 2 )
......@@ -89,7 +106,7 @@ class KpointCmd : public Cmd
}
if ( onpe0 )
{
cout << " <-- kpoint list: reciprocal lattice coordinates" << endl;
cout << " <!-- kpoint list: reciprocal lattice coordinates" << endl;
for ( int ikp = 0; ikp < s->wf.nkp(); ikp++ )
{
cout << " "
......
src/Makefile
......@@ -72,6 +72,19 @@ CXXFLAGS += -DTARGET='"$(TARGET)"'
ExtForceSet.o ExtForce.o PairExtForce.o AtomicExtForce.o \
GlobalExtForce.o sampling.o
$(LD) $(DFLAGS) -o $@ $^ $(LDFLAGS)
testSampleReader: testSampleReader.o AtomSet.o Atom.o Species.o \
Wavefunction.o SlaterDet.o \
Basis.o FourierTransform.o Matrix.o Context.o \
sinft.o spline.o UnitCell.o \
Base64Transcoder.o Constraint.o ConstraintSet.o DistanceConstraint.o \
AngleConstraint.o TorsionConstraint.o PositionConstraint.o \
ExtForceSet.o ExtForce.o PairExtForce.o AtomicExtForce.o \
GlobalExtForce.o sampling.o \
SampleReader.o StructuredDocumentHandler.o \
SampleHandler.o AtomSetHandler.o WavefunctionHandler.o \
SpeciesReader.o SpeciesHandler.o \
XMLGFPreprocessor.o Base64Transcoder.o
$(LD) $(DFLAGS) -o $@ $^ $(LDFLAGS)
testChargeDensity: testChargeDensity.o ChargeDensity.o \
Wavefunction.o SlaterDet.o \
Basis.o FourierTransform.o Matrix.o UnitCell.o Context.o \
......@@ -693,7 +706,8 @@ Vext.o: Sample.h AtomSet.h Context.h blacs.h Atom.h D3vector.h UnitCell.h
Vext.o: D3tensor.h blas.h ConstraintSet.h ExtForceSet.h Wavefunction.h
Vext.o: Control.h ExternalPotential.h ChargeDensity.h Timer.h
Wavefunction.o: Wavefunction.h D3vector.h UnitCell.h SlaterDet.h Context.h
Wavefunction.o: blacs.h Basis.h Matrix.h Timer.h jacobi.h SharedFilePtr.h
Wavefunction.o: blacs.h Basis.h Matrix.h Timer.h FourierTransform.h jacobi.h
Wavefunction.o: SharedFilePtr.h
Wavefunction.o: D3vector.h UnitCell.h
WavefunctionHandler.o: WavefunctionHandler.h StructureHandler.h UnitCell.h
WavefunctionHandler.o: D3vector.h Wavefunction.h SlaterDet.h Context.h
......@@ -778,6 +792,10 @@ testMatrix.o: Timer.h Context.h blacs.h Matrix.h
testSample.o: Context.h blacs.h SlaterDet.h Basis.h D3vector.h UnitCell.h
testSample.o: Matrix.h Timer.h Sample.h AtomSet.h Atom.h D3tensor.h blas.h
testSample.o: ConstraintSet.h ExtForceSet.h Wavefunction.h Control.h
testSampleReader.o: Context.h blacs.h SlaterDet.h Basis.h D3vector.h
testSampleReader.o: UnitCell.h Matrix.h Timer.h Sample.h AtomSet.h Atom.h
testSampleReader.o: D3tensor.h blas.h ConstraintSet.h ExtForceSet.h
testSampleReader.o: Wavefunction.h Control.h SampleReader.h
testSlaterDet.o: Context.h blacs.h SlaterDet.h Basis.h D3vector.h UnitCell.h
testSlaterDet.o: Matrix.h Timer.h FourierTransform.h
testSpecies.o: Species.h SpeciesReader.h
......
src/Matrix.C
......@@ -56,6 +56,7 @@ using namespace std;
#define pdtrsm pdtrsm_
#define pztrsm pztrsm_
#define pdtrtrs pdtrtrs_
#define pztrtrs pztrtrs_
#define pdpotrf pdpotrf_
#define pzpotrf pzpotrf_
#define pdpotri pdpotri_
......@@ -67,6 +68,7 @@ using namespace std;
#define pzheev pzheev_
#define pzheevd pzheevd_
#define pdtrtri pdtrtri_
#define pztrtri pztrtri_
#define pdlatra pdlatra_
#define pdlacp2 pdlacp2_
#define pdlacp3 pdlacp3_
......@@ -102,6 +104,7 @@ using namespace std;
#define dtrmm dtrmm_
#define dtrsm dtrsm_
#define dtrtri dtrtri_
#define ztrtri ztrtri_
#define ztrsm ztrsm_
#define dtrtrs dtrtrs_
#define dpotrf dpotrf_
......@@ -197,6 +200,9 @@ extern "C"
void pdtrtrs(const char*, const char*, const char*, const int*, const int*,
const double*, const int*, const int*, const int*,
double*, const int*, const int*, const int*, int*);
void pztrtrs(const char*, const char*, const char*, const int*, const int*,
const complex<double>*, const int*, const int*, const int*,
complex<double>*, const int*, const int*, const int*, int*);
void pigemr2d(const int*,const int*,
const int*,const int*,const int*, const int*,
int*,const int*,const int*,const int*,const int*);
......@@ -249,6 +255,8 @@ extern "C"
int* iwork, int* liwork, int* info);
void pdtrtri(const char*, const char*, const int*, double*,
const int*, const int*, const int*, int*);
void pztrtri(const char*, const char*, const int*, complex<double>*,
const int*, const int*, const int*, int*);
void pdgetrf(const int* m, const int* n, double* val,
int* ia, const int* ja, const int* desca, int* ipiv, int* info);
void pzgetrf(const int* m, const int* n, complex<double>* val,
......@@ -1882,6 +1890,41 @@ void DoubleMatrix::trtrs(char uplo, char trans, char diag,
}
////////////////////////////////////////////////////////////////////////////////
// Solves a triangular system of the form A * X = B or
// A**T * X = B, where A is a triangular matrix of order N,
// and B is an N-by-NRHS matrix.
// Output in B.
////////////////////////////////////////////////////////////////////////////////
void ComplexMatrix::trtrs(char uplo, char trans, char diag,
ComplexMatrix& b) const
{
int info;
if ( active() )
{
assert(m_==n_);
#ifdef SCALAPACK
int ione=1;
pztrtrs(&uplo, &trans, &diag, &m_, &b.n_,
val, &ione, &ione, desc_,
b.val, &ione, &ione, b.desc_, &info);
#else
ztrtrs(&uplo, &trans, &diag, &m_, &b.n_, val, &m_,
b.val, &b.m_, &info);
#endif
if(info!=0)
{
cout <<" ComplexMatrix::trtrs, info=" << info << endl;
#ifdef USE_MPI
MPI_Abort(MPI_COMM_WORLD, 2);
#else
exit(2);
#endif
}
}
}
////////////////////////////////////////////////////////////////////////////////
// LU decomposition of a double matrix
////////////////////////////////////////////////////////////////////////////////
void DoubleMatrix::lu(valarray<int>& ipiv)
......@@ -2269,6 +2312,31 @@ void DoubleMatrix::trtri(char uplo, char diag)
}
}
void ComplexMatrix::trtri(char uplo, char diag)
{
int info;
if ( active() )
{
assert(m_==n_);
#ifdef SCALAPACK
int ione=1;
pztrtri(&uplo, &diag, &m_, val, &ione, &ione, desc_, &info);
#else
ztrtri(&uplo, &diag, &m_, val, &m_, &info);
#endif
if(info!=0)
{
cout << " Matrix::trtri, info=" << info << endl;
#ifdef USE_MPI
MPI_Abort(MPI_COMM_WORLD, 2);
#else
exit(2);
#endif
}
}
}
////////////////////////////////////////////////////////////////////////////////
// Polar decomposition A = UH
// Replace *this with its orthogonal polar factor U
......@@ -2352,6 +2420,88 @@ void DoubleMatrix::polar(double tol, int maxiter)
}
////////////////////////////////////////////////////////////////////////////////
// Polar decomposition A = UH (complex case)
// Replace *this with its unitary polar factor U
// return when iter > maxiter or ||I - X^H*X|| < tol
////////////////////////////////////////////////////////////////////////////////
void ComplexMatrix::polar(double tol, int maxiter)
{
ComplexMatrix x(ctxt_,m_,n_,mb_,nb_);
ComplexMatrix xp(ctxt_,m_,n_,mb_,nb_);
ComplexMatrix q(ctxt_,n_,n_,nb_,nb_);
ComplexMatrix qt(ctxt_,n_,n_,nb_,nb_);
ComplexMatrix t(ctxt_,n_,n_,nb_,nb_);
#ifdef SCALAPACK
double qnrm2 = numeric_limits<double>::max();
int iter = 0;
x = *this;
while ( iter < maxiter && qnrm2 > tol )
{
// q = I - x^T * x
q.identity();
q.herk('l','c',-1.0,x,1.0);
q.symmetrize('l');
double qnrm2 = q.nrm2();
#ifdef DEBUG
if ( ctxt_.onpe0() )
cout << " ComplexMatrix::polar: qnrm2 = " << qnrm2 << endl;
#endif
// choose Bjork-Bowie or Higham iteration depending on q.nrm2
// threshold value
// see A. Bjork and C. Bowie, SIAM J. Num. Anal. 8, 358 (1971) p.363
if ( qnrm2 < 1.0 )
{
// Bjork-Bowie iteration
// compute xp = x * ( I + 0.5*q * ( I + 0.75 * q ) )
// t = ( I + 0.75 * q )
t.identity();
t.axpy(0.75,q);
// compute q*t
qt.gemm('n','n',1.0,q,t,0.0);
// xp = x * ( I + 0.5*q * ( I + 0.75 * q ) )
// = x * ( I + 0.5 * qt )
// Use t to store (I + 0.5 * qt)
t.identity();
t.axpy(0.5,qt);
// t now contains (I + 0.5 * qt)
// xp = x * t
xp.gemm('n','n',1.0,x,t,0.0);
// update x
x = xp;
}
else
{
// Higham iteration
assert(m_==n_);
//if ( ctxt_.onpe0() )
// cout << " ComplexMatrix::polar: using Higham algorithm" << endl;
// t = X^H
t.transpose(1.0,x,0.0);
t.inverse();
// t now contains X^-H
// xp = 0.5 * ( x + x^-H );
for ( int i = 0; i < x.size(); i++ )
x[i] = 0.5 * ( x[i] + t[i] );
}
iter++;
}
*this = x;
#else
#error "ComplexMatrix::polar only implemented with SCALAPACK"
#endif
}
////////////////////////////////////////////////////////////////////////////////
// estimate the reciprocal of the condition number (in the 1-norm) of a
// real symmetric positive definite matrix using the Cholesky factorization
// A = U**T*U or A = L*L**T computed by DoubleMatrix::potrf
......
src/Matrix.h
......@@ -537,6 +537,9 @@ class ComplexMatrix
// Inverse of a symmetric matrix from Cholesky factor
void potri(char uplo);
// Polar decomposition, tolerance ||I-X^H*X||<tol or iter<maxiter
void polar(double tol, int maxiter);
// LU decomposition
void lu(std::valarray<int>& ipiv);
......
src/PSDAWavefunctionStepper.C
......@@ -62,9 +62,12 @@ void PSDAWavefunctionStepper::update(Wavefunction& dwf)
}
else
{
// not implemented in the complex case
cout << "PSDA is not implemented for complex wave functions" << endl;
assert(false);
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);
}
}
}
......@@ -123,21 +126,25 @@ void PSDAWavefunctionStepper::update(Wavefunction& dwf)
a += f * delta_f;
b += delta_f * delta_f;
}
// correct for double counting of asum and bsum on first row
// factor 2.0: G and -G
a *= 2.0;
b *= 2.0;
if ( wf_.sdcontext()->myrow() == 0 )
if ( wf_.sd(ispin,ikp)->basis().real() )
{
for ( int n = 0; n < nloc; n++ )
// correct for double counting of asum and bsum on first row
// factor 2.0: G and -G
a *= 2.0;
b *= 2.0;
if ( wf_.sdcontext()->myrow() == 0 )
{
const int i = 2*mloc*n;
const double f0 = dc[i];
const double f1 = dc[i+1];
const double delta_f0 = f0 - dc_last[i];
const double delta_f1 = f1 - dc_last[i+1];
a -= f0 * delta_f0 + f1 * delta_f1;
b -= delta_f0 * delta_f0 + delta_f1 * delta_f1;
for ( int n = 0; n < nloc; n++ )
{
const int i = 2*mloc*n;
const double f0 = dc[i];
const double f1 = dc[i+1];
const double delta_f0 = f0 - dc_last[i];
const double delta_f1 = f1 - dc_last[i+1];
a -= f0 * delta_f0 + f1 * delta_f1;
b -= delta_f0 * delta_f0 + delta_f1 * delta_f1;
}
}
}
......@@ -151,9 +158,6 @@ void PSDAWavefunctionStepper::update(Wavefunction& dwf)
if ( b != 0.0 )
theta = - a / b;
if ( wf_.sdcontext()->onpe0() )
cout << " Anderson extrapolation: theta=" << theta;
if ( theta < -1.0 )
{
theta = 0.0;
......@@ -161,9 +165,6 @@ void PSDAWavefunctionStepper::update(Wavefunction& dwf)
theta = min(2.0,theta);
if ( wf_.sdcontext()->onpe0() )
cout <<" (" << theta << ")" << endl;
// extrapolation
for ( int i = 0; i < 2*mloc*nloc; i++ )
{
......
src/PositionConstraint.C
......@@ -102,9 +102,10 @@ ostream& PositionConstraint::print( ostream &os )
os << "\" atoms=\"" << name1_ << "\"\n";
os.setf(ios::fixed,ios::floatfield);
os.setf(ios::right,ios::adjustfield);
os << " value=\"" << setprecision(6) << 0;
os << "\" velocity=\"" << setprecision(6) << 0 << "\"\n";
os << " force=\"" << setprecision(6) << force_;
os << "\" weight=\"" << setprecision(6) << weight_ << "\"/>";
os << " velocity=\"" << setprecision(6) << 0 << "\"";
os << " weight=\"" << setprecision(6) << weight_ << "\">\n";
os << " <value> " << setprecision(6) << 0 << " </value>";
os << " <force> " << setprecision(6) << force_ << " </force>\n";
os << " </constraint>";
return os;
}
src/RseedCmd.h
......@@ -41,7 +41,7 @@ class RseedCmd : public Cmd
return
"\n rseed\n\n"
" syntax: rseed [seed_value]\n\n"
" The rseed command initializes the random number generator."
" The rseed command initializes the random number generator.\n"
" If no argument is given, the time() function is used as a seed value"
"\n\n";
}
......
src/RunCmd.C
......@@ -50,6 +50,12 @@ int RunCmd::action(int argc, char **argv)
cout << " RunCmd: ecut = 0.0, cannot run" << endl;
return 1;
}
if ( s->wf.cell().volume() == 0.0 )
{