Commit 85b9d4de by Francois Gygi

rel1_48_1


git-svn-id: http://qboxcode.org/svn/qb/trunk@710 cba15fb0-1239-40c8-b417-11db7ca47a34
parent eb5f565d
......@@ -15,7 +15,7 @@
// AndersonMixer.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: AndersonMixer.C,v 1.8 2009-03-08 01:10:30 fgygi Exp $
// $Id: AndersonMixer.C,v 1.9 2009-08-14 17:06:43 fgygi Exp $
#include "AndersonMixer.h"
#include "blas.h"
......@@ -38,13 +38,17 @@ void AndersonMixer::update(double* x, double* f, double* xbar, double* fbar)
//
// Computes the pair (xbar,fbar) using pairs (x,f) used
// in previous updates, according to the Anderson algorithm.
//
// Tikhonov regularization parameter
const double tikhonov_parameter = 0.2;
// increment index of current vector
k_ = ( k_ + 1 ) % ( nmax_ + 1 );
// increment current number of vectors
if ( n_ < nmax_ ) n_++;
// save vectors
for ( int i = 0; i < m_; i++ )
{
x_[k_][i] = x[i];
......@@ -78,95 +82,146 @@ void AndersonMixer::update(double* x, double* f, double* xbar, double* fbar)
sum += (f_[k_][l] - f_[kmi][l]) * (f_[k_][l] - f_[kmj][l]);
a[i+j*n_] = sum;
}
double bsum = 0.0;
for ( int l = 0; l < m_; l++ )
bsum += ( f_[k_][l] - f_[kmi][l] ) * f_[k_][l];
b[i] = bsum;
}
#if 0
// print matrix a and rhs b
for ( int i = 0; i < n_; i++ )
for ( int j = 0; j <=i; j++ )
cout << "a("<<i<<","<<j<<")=" << a[i+j*n_] << endl;
for ( int i = 0; i < n_; i++ )
cout << "b("<<i<<")=" << b[i] << endl;
#endif
if ( pctxt_ != 0 )
{
pctxt_->dsum(n_*n_,1,&a[0],n_*n_);
pctxt_->dsum(n_,1,&b[0],n_);
}
// solve the linear system a * theta = b
const bool diag = false;
if ( diag )
{
// solve the linear system using eigenvalues and eigenvectors
// compute eigenvalues of a
char jobz = 'V';
char uplo = 'L';
valarray<double> w(n_);
int lwork = 3*n_;
valarray<double> work(lwork);
int info;
dsyev(&jobz,&uplo,&n_,&a[0],&n_,&w[0],&work[0],&lwork,&info);
for ( int i = 0; i < n_; i++ )
cout << w[i] << " ";
cout << endl;
if ( info != 0 )
{
cerr << " AndersonMixer: Error in dsyev" << endl;
cerr << " info = " << info << endl;
exit(1);
}
for ( int i = 0; i < n_; i++ )
a[i+i*n_] += tikhonov_parameter;
// solve for theta
// theta_i = sum_j
for ( int k = 0; k < n_; k++ )
#if 0
// print matrix a and rhs b
if ( pctxt_ != 0 )
{
for ( int ip =0; ip < pctxt_->size(); ip++ )
{
// correct only if eigenvalue w[k] is large enough compared to the
// largest eigenvalue
if ( w[n_-1] < 5.0 * w[k] )
pctxt_->barrier();
if ( pctxt_->mype() == ip )
{
const double fac = 1.0 / w[k];
// print matrix a and rhs b
double anrm = 0.0;
for ( int i = 0; i < n_; i++ )
for ( int j = 0; j < n_; j++ )
theta[i] += fac * a[i+k*n_] * a[j+k*n_] * b[j];
for ( int j = 0; j <=i; j++ )
{
cout << pctxt_->mype() << ": "
<< "a("<<i<<","<<j<<")=" << a[i+j*n_] << endl;
anrm += a[i+j*n_]*a[i+j*n_];
}
for ( int i = 0; i < n_; i++ )
cout << pctxt_->mype() << ": "
<< "b("<<i<<")=" << b[i] << endl;
//cout << " AndersonMixer: n=" << n_ << " anorm = " << anrm << endl;
}
}
}
else
#endif
// solve the linear system a * theta = b
// solve on task 0 and bcast result
if ( pctxt_ == 0 || pctxt_->onpe0() )
{
// solve the linear system directly
char uplo = 'L';
int nrhs = 1;
valarray<int> ipiv(n_);
valarray<double> work(n_);
int info;
dsysv(&uplo,&n_,&nrhs,&a[0],&n_,&ipiv[0],&b[0],&n_,&work[0],&n_,&info);
if ( info != 0 )
const bool diag = false;
if ( diag )
{
cerr << " AndersonMixer: Error in dsysv" << endl;
cerr << " info = " << info << endl;
exit(1);
// solve the linear system using eigenvalues and eigenvectors
// compute eigenvalues of a
char jobz = 'V';
char uplo = 'L';
valarray<double> w(n_);
int lwork = 3*n_;
valarray<double> work(lwork);
int info;
dsyev(&jobz,&uplo,&n_,&a[0],&n_,&w[0],&work[0],&lwork,&info);
if ( pctxt_ != 0 )
{
if ( pctxt_->onpe0() )
{
cout << "AndersonMixer: eigenvalues: ";
for ( int i = 0; i < n_; i++ )
cout << w[i] << " ";
cout << endl;
}
}
if ( info != 0 )
{
cerr << " AndersonMixer: Error in dsyev" << endl;
cerr << " info = " << info << endl;
exit(1);
}
// solve for theta
// theta_i = sum_j
for ( int k = 0; k < n_; k++ )
{
// correct only if eigenvalue w[k] is large enough compared to the
// largest eigenvalue
if ( w[n_-1] < 5.0 * w[k] )
{
const double fac = 1.0 / w[k];
for ( int i = 0; i < n_; i++ )
for ( int j = 0; j < n_; j++ )
theta[i] += fac * a[i+k*n_] * a[j+k*n_] * b[j];
}
}
}
// the vector b now contains the solution
theta = b;
if ( pctxt_ != 0 )
else
{
if ( pctxt_->onpe0() )
// solve the linear system directly
char uplo = 'L';
int nrhs = 1;
valarray<int> ipiv(n_);
valarray<double> work(n_);
int info;
dsysv(&uplo,&n_,&nrhs,&a[0],&n_,&ipiv[0],&b[0],&n_,&work[0],&n_,&info);
if ( info != 0 )
{
cout << " AndersonMixer: theta = ";
for ( int i = 0; i < theta.size(); i++ )
cout << theta[i] << " ";
cout << endl;
cerr << " AndersonMixer: Error in dsysv" << endl;
cerr << " info = " << info << endl;
exit(1);
}
// the vector b now contains the solution
theta = b;
}
cout << " AndersonMixer: theta = ";
for ( int i = 0; i < theta.size(); i++ )
cout << theta[i] << " ";
cout << endl;
}
// broadcast theta from task 0
if ( pctxt_ != 0 )
{
if ( pctxt_->onpe0() )
pctxt_->dbcast_send(n_,1,&theta[0],n_);
else
pctxt_->dbcast_recv(n_,1,&theta[0],n_,0,0);
}
#if 0
for ( int ip = 0; ip < pctxt_->size(); ip++ )
{
pctxt_->barrier();
if ( pctxt_->mype() == ip )
{
cout << pctxt_->mype() << ": ";
cout << " AndersonMixer: theta = ";
for ( int i = 0; i < theta.size(); i++ )
cout << theta[i] << " ";
cout << endl;
}
}
#endif
} // if n_ > 0
// fbar = f[k] + sum_i theta_i * ( f[k] - f[kmi] )
......
......@@ -15,7 +15,7 @@
// AtomCmd.h:
//
////////////////////////////////////////////////////////////////////////////////
// $Id: AtomCmd.h,v 1.10 2008-09-08 15:56:17 fgygi Exp $
// $Id: AtomCmd.h,v 1.11 2009-08-14 17:06:43 fgygi Exp $
#ifndef ATOMCMD_H
#define ATOMCMD_H
......@@ -42,8 +42,8 @@ class AtomCmd : public Cmd
" syntax: atom name species x y z [vx vy vz]\n\n"
" The atom command defines a new atom and adds it to the atom list.\n"
" The name can be any character string, the species must be the name\n"
" of a file containing the definition of the pseudopotential for this\n"
" atomic species. The position of the atom is specified by x y and z.\n"
" of a species previously defined using the species command.\n"
" The position of the atom is specified by x y and z in atomic units.\n"
" Optionally, the atom velocity can be specified by vx vy and vz.\n\n";
}
......
......@@ -15,7 +15,7 @@
// BOSampleStepper.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: BOSampleStepper.C,v 1.50 2009-05-15 04:37:42 fgygi Exp $
// $Id: BOSampleStepper.C,v 1.51 2009-08-14 17:06:43 fgygi Exp $
#include "BOSampleStepper.h"
#include "EnergyFunctional.h"
......@@ -219,8 +219,8 @@ void BOSampleStepper::step(int niter)
const double q0_kerker2 = q0_kerker * q0_kerker;
for ( int i=0; i < wkerker.size(); i++ )
{
wkerker[i] = g2[i] / ( g2[i] + q0_kerker2 );
wkerker_inv[i] = g2i[i] * ( g2[i] + q0_kerker2 );
wkerker[i] = sqrt(g2[i] / ( g2[i] + q0_kerker2 ));
wkerker_inv[i] = sqrt(g2i[i] * ( g2[i] + q0_kerker2 ));
}
}
else
......@@ -579,6 +579,9 @@ void BOSampleStepper::step(int niter)
assert(cd_.rhog.size()==1); // works for nspin=1 only
wf_stepper->preprocess();
if ( anderson_charge_mixing )
mixer.restart();
for ( int itscf = 0; itscf < nitscf_; itscf++ )
{
if ( nite_ > 1 && onpe0 )
......@@ -609,14 +612,18 @@ void BOSampleStepper::step(int niter)
// Anderson acceleration
if ( anderson_charge_mixing )
{
#if 1
// row weighting of LS calculation (preconditioning)
for ( int i=0; i < drhog.size(); i++ )
drhog[i] *= wkerker_inv[i];
#endif
mixer.update((double*)&rhog_in[0],(double*)&drhog[0],
(double*)&rhobar[0],(double*)&drhobar[0]);
#if 1
for ( int i=0; i < drhog.size(); i++ )
drhobar[i] *= wkerker[i];
drhobar[i] *= wkerker[i]*wkerker[i];
#endif
for ( int i=0; i < rhog_in.size(); i++ )
rhog_in[i] = rhobar[i] + alpha * drhobar[i];
......@@ -624,7 +631,7 @@ void BOSampleStepper::step(int niter)
else
{
for ( int i=0; i < rhog_in.size(); i++ )
rhog_in[i] += alpha * drhog[i] * wkerker[i];
rhog_in[i] += alpha * drhog[i] * wkerker[i] * wkerker[i];
}
cd_.rhog[0] = rhog_in;
......
......@@ -11,7 +11,7 @@
# or <http://www.gnu.org/licenses/>.
#
#-------------------------------------------------------------------------------
# $Id: Makefile,v 1.63 2009-06-29 10:08:59 fgygi Exp $
# $Id: Makefile,v 1.64 2009-08-14 17:06:43 fgygi Exp $
#------------------------------------------------------------------------------
#
include $(TARGET).mk
......@@ -92,7 +92,8 @@ CXXFLAGS += -DTARGET='"$(TARGET)"'
Species.o SpeciesReader.o sinft.o spline.o SpeciesHandler.o \
SampleHandler.o WavefunctionHandler.o AtomSetHandler.o \
Base64Transcoder.o ConstraintSet.o DistanceConstraint.o \
AngleConstraint.o TorsionConstraint.o XMLGFPreprocessor.o \
AngleConstraint.o TorsionConstraint.o PositionConstraint.o \
XMLGFPreprocessor.o \
SampleWriter.o qbox_xmlns.o isodate.o
$(LD) $(DFLAGS) -o $@ $^ $(LDFLAGS)
plotSample: plotSample.o \
......@@ -159,6 +160,8 @@ CXXFLAGS += -DTARGET='"$(TARGET)"'
$(LD) $(DFLAGS) -o $@ $^ $(LDFLAGS)
testMemParse: testMemParse.o
$(LD) $(DFLAGS) -o $@ $^ $(LDFLAGS)
upf2qso: upf2qso.o PeriodicTable.o spline.o
$(LD) $(DFLAGS) -o $@ $^ $(LDFLAGS)
xmlSpecies: xmlSpecies.o qbox_xmlns.o
$(LD) $(DFLAGS) -o $@ $^ $(LDFLAGS)
xmlget: xmlget.o
......
......@@ -15,7 +15,7 @@
// PlotCmd.C:
//
////////////////////////////////////////////////////////////////////////////////
// $Id: PlotCmd.C,v 1.2 2009-06-29 10:08:20 fgygi Exp $
// $Id: PlotCmd.C,v 1.3 2009-08-14 17:06:43 fgygi Exp $
#include "PlotCmd.h"
#include "isodate.h"
......@@ -37,8 +37,6 @@
#include "Atom.h"
#include "ChargeDensity.h"
//#include "atomic_properties.h"
using namespace std;
////////////////////////////////////////////////////////////////////////////////
......
......@@ -15,7 +15,7 @@
// SlaterDet.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: SlaterDet.C,v 1.54 2009-06-29 09:59:06 fgygi Exp $
// $Id: SlaterDet.C,v 1.55 2009-08-14 17:06:43 fgygi Exp $
#include "SlaterDet.h"
#include "FourierTransform.h"
......@@ -1629,7 +1629,7 @@ void SlaterDet::write(SharedFilePtr& sfp, string encoding, double weight, int is
<< " encoding=\"text\">" << endl;
}
int count = 0;
for ( int k = 0; k < ft.np2(); k++ )
for ( int k = 0; k < ft.np2_loc(); k++ )
for ( int j = 0; j < ft.np1(); j++ )
for ( int i = 0; i < ft.np0(); i++ )
{
......
To do:
Fix the ndim>0 Anderson charge mixing to avoid singular systems
1.48.1_beta_x try Tikhonov regularization of the Anderson LS problem
--------------------------------------------------------------------------------
rel1_48_1 candidate
SlaterDet.C: fix bug in SlaterDet::write for -text mode
BOSampleStepper.C: modified Kerker weights. Added mixer.reset() call.
AtomCmd.h: fixed text of help message.
AndersonMixer.C: added Tikhonov regularization. Bcast results from task 0.
PlotCmd.C: cleanup include.
--------------------------------------------------------------------------------
rel1_48_0
PlotCmd.[Ch]: added plot command. Generates cube format or xyz
......
......@@ -15,10 +15,10 @@
// release.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: release.C,v 1.72 2009-06-29 10:00:19 fgygi Exp $
// $Id: release.C,v 1.73 2009-08-14 17:06:43 fgygi Exp $
#include "release.h"
std::string release(void)
{
return std::string("1.48.0");
return std::string("1.48.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