Commit bd8bc088 by Francois Gygi

Merge branch 'develop'

parents 13d04673 23da6945
......@@ -26,19 +26,19 @@
CXX=mpiCC
LD=$(CXX)
PLTFLAGS += -DIA32 -DUSE_MPI -DUSE_FFTW2 -D_LARGEFILE_SOURCE \
PLTFLAGS += -DIA32 -DUSE_MPI -DUSE_FFTW3 -D_LARGEFILE_SOURCE \
-D_FILE_OFFSET_BITS=64 -DADD_ \
-DAPP_NO_THREADS -DXML_USE_NO_THREADS -DUSE_XERCES -DXERCESC_3 \
-DSCALAPACK -DUSE_UUID
INCLUDE = -I$(FFTWDIR)
CXXFLAGS= -g -O3 -Wunused -D$(PLT) $(INCLUDE) $(PLTFLAGS) $(DFLAGS)
CXXFLAGS= -g -O3 -D$(PLT) $(INCLUDE) $(PLTFLAGS) $(DFLAGS)
LIBPATH = -L$(SCALAPACKDIR) -L$(LAPACKDIR) -L$(BLASDIR)
LIBS = -lfftw -lscalapack -llapack -lblas -lm -lgfortran \
-lxerces-c -lpthread -lmpi_cxx -lstdc++ -lcurl -liconv
LIBS = -lfftw3 -lscalapack -llapack -lblas -lm -lgfortran \
-lxerces-c -lpthread -lstdc++
LDFLAGS = $(LIBPATH) $(LIBS)
#-------------------------------------------------------------------------------
......@@ -76,27 +76,27 @@ void BOSampleStepper::initialize_density(void)
double atom_radius[] =
{
0.00, 0.80, 0.59, 3.16, 2.12, // null H He Li Be
1.64, 1.27, 1.06, 0.91, 0.79, // B C N O F
0.72, 3.59, 2.74, 2.23, 2.10, // Ne Na Mg Al Si
1.85, 1.66, 1.49, 1.34, 4.59, // P S Cl Ar K
3.67, 3.48, 3.33, 3.23, 3.14, // Ca Sc Ti V Cr
3.04, 2.95, 2.87, 2.82, 2.74, // Mn Fe Co Ni Cu
2.68, 2.57, 2.36, 2.15, 1.95, // Zn Ga Ge As Se
1.78, 1.66, 5.01, 4.14, 4.01, // Br Kr Rb Sr Y
3.89, 3.74, 3.59, 3.46, 3.36, // Zr Nb Mo Tc Ru
3.27, 3.19, 3.12, 3.04, 2.95, // Rh Pd Ag Cd In
2.74, 2.51, 2.32, 2.17, 2.04, // Sn Sb Te I Xe
5.63, 4.78, 0.00, 0.00, 4.67, // Cs Ba La Ce Pr
3.89, 3.87, 4.50, 4.37, 4.40, // Nd Pm Sm Eu Gd
4.25, 4.31, 0.00, 4.27, 4.20, // Tb Dy Ho Er Tm
4.20, 4.10, 3.93, 3.78, 3.65, // Yb Lu Hf Ta W
3.55, 3.50, 3.40, 3.34, 3.29, // Re Os Ir Pt Au
3.23, 2.95, 2.91, 2.70, 2.55, // Hg Tl Pb Bi Po
4.00, 2.27, 4.00, 4.00, 4.00, // At Rn Fr Ra Ac
4.00, 4.00, 4.00, 4.00, 4.00, // Th Pa U Np Pu
4.00, 4.00, 4.00, 4.00, 4.00, // Am Cm Bk Cf Es
4.00, 4.00, 4.00, 4.00 // Fm Md No Lr
0.000, 1.542, 0.931, 1.727, 1.636, // null H He Li Be
1.845, 1.538, 1.311, 1.141, 1.014, // B C N O F
0.983, 1.238, 1.347, 1.376, 2.151, // Ne Na Mg Al Si
1.927, 1.733, 1.563, 1.429, 1.546, // P S Cl Ar K
1.663, 1.604, 1.516, 1.434, 1.377, // Ca Sc Ti V Cr
1.310, 1.245, 1.201, 1.162, 1.098, // Mn Fe Co Ni Cu
1.077, 1.331, 1.415, 2.015, 1.880, // Zn Ga Ge As Se
1.749, 1.630, 1.705, 1.819, 1.794, // Br Kr Rb Sr Y
1.728, 1.664, 1.589, 1.523, 1.461, // Zr Nb Mo Tc Ru
1.410, 1.348, 1.306, 1.303, 1.554, // Rh Pd Ag Cd In
1.609, 1.611, 1.530, 1.514, 1.464, // Sn Sb Te I Xe
1.946, 1.967, 1.943, 1.930, 1.920, // Cs Ba La Ce Pr
1.910, 1.900, 1.890, 1.880, 1.870, // Nd Pm Sm Eu Gd
1.860, 1.850, 1.840, 1.830, 1.820, // Tb Dy Ho Er Tm
1.810, 1.800, 1.701, 1.658, 1.606, // Yb Lu Hf Ta W
1.550, 1.500, 1.446, 1.398, 1.355, // Re Os Ir Pt Au
1.314, 1.624, 1.659, 1.634, 1.620, // Hg Tl Pb Bi Po
1.600, 1.500, 1.600, 1.600, 1.600, // At Rn Fr Ra Ac
1.600, 1.600, 1.600, 1.600, 1.600, // Th Pa U Np Pu
1.600, 1.600, 1.600, 1.600, 1.600, // Am Cm Bk Cf Es
1.600, 1.600, 1.600, 1.600 // Fm Md No Lr
};
const AtomSet& atoms = s_.atoms;
......@@ -115,8 +115,7 @@ void BOSampleStepper::initialize_density(void)
const int zval = s->zval();
const int atomic_number = s->atomic_number();
assert(atomic_number < sizeof(atom_radius)/sizeof(double));
// scaling factor 2.0 in next line is empirically adjusted
double rc = 2.0 * atom_radius[atomic_number];
double rc = atom_radius[atomic_number];
for ( int ig = 0; ig < ngloc; ig++ )
{
double arg = 0.25 * rc * rc * g2[ig];
......@@ -745,6 +744,7 @@ void BOSampleStepper::step(int niter)
mixer.restart();
double ehart, ehart_m;
double delta_ehart;
bool scf_converged = false;
int itscf = 0;
double etotal = 0.0, etotal_m = 0.0, etotal_mm = 0.0;
......@@ -862,16 +862,36 @@ void BOSampleStepper::step(int niter)
// non-self-consistent loop
// repeat until the change in eigenvalue_sum is smaller than a
// fraction of the change in Hartree energy in the last scf iteration
// fraction delta_ratio of the change in Hartree energy delta_ehart
// in the last scf iteration
bool nonscf_converged = false;
if ( itscf > 0 )
if ( itscf == 0 )
{
ehart = ef_.ehart();
delta_ehart = 0.0;
}
else if ( itscf == 1 )
{
ehart_m = ehart;
ehart = ef_.ehart();
double delta_ehart = 0.0;
if ( itscf > 0 )
ehart = ef_.ehart();
delta_ehart = fabs(ehart - ehart_m);
// if ( onpe0 && nite_ > 0 )
// cout << " delta_ehart = " << delta_ehart << endl;
}
else
{
// itscf > 1
// only allow decrease in delta_ehart
ehart_m = ehart;
ehart = ef_.ehart();
delta_ehart = min(delta_ehart,fabs(ehart - ehart_m));
}
#if DEBUG
if ( onpe0 )
{
cout << " BOSampleStepper::step: delta_ehart: "
<< delta_ehart << endl;
}
#endif
int ite = 0;
double etotal_int;
......@@ -897,7 +917,7 @@ void BOSampleStepper::step(int niter)
eigenvalue_sum = real(s_.wf.dot(dwf));
if ( onpe0 )
cout << " <eigenvalue_sum> "
cout << " <eigenvalue_sum> " << setprecision(8)
<< eigenvalue_sum << " </eigenvalue_sum>" << endl;
tmap["wf_update"].start();
......@@ -907,15 +927,14 @@ void BOSampleStepper::step(int niter)
// compare delta_eig_sum only after first iteration
if ( ite > 0 )
{
const double delta_ratio = 0.1;
double delta_eig_sum = fabs(eigenvalue_sum - eigenvalue_sum_m);
nonscf_converged |= (delta_eig_sum < 0.01 * delta_ehart);
#ifdef DEBUG
nonscf_converged |= (delta_eig_sum < delta_ratio * delta_ehart);
#if DEBUG
if ( onpe0 )
{
cout << " BOSampleStepper::step delta_eig_sum: "
<< delta_eig_sum << endl;
cout << " BOSampleStepper::step: delta_ehart: "
<< delta_ehart << endl;
}
#endif
}
......
......@@ -47,6 +47,8 @@ class CPSampleStepper : public SampleStepper
mutable TimerMap tmap;
void step(int niter);
ChargeDensity& cd(void) { return cd_; }
EnergyFunctional& ef(void) { return ef_; }
CPSampleStepper(Sample& s);
~CPSampleStepper();
......
......@@ -724,7 +724,8 @@ double EnergyFunctional::energy(bool compute_hpsi, Wavefunction& dwf,
if ( compute_stress )
{
valarray<double> sigma_ext(s_.ctrl.ext_stress,6);
const double pext = (sigma_ext[0]+sigma_ext[1]+sigma_ext[2])/3.0;
const double gpa = 29421.5;
const double pext = (sigma_ext[0]+sigma_ext[1]+sigma_ext[2])/(3.0*gpa);
epv_ = pext * omega;
enthalpy_ += epv_;
}
......@@ -880,7 +881,6 @@ double EnergyFunctional::energy(bool compute_hpsi, Wavefunction& dwf,
if ( debug_stress && s_.ctxt_.onpe0() )
{
//const double gpa = 29421.5;
cout.setf(ios::fixed,ios::floatfield);
cout.setf(ios::right,ios::adjustfield);
cout << setprecision(8);
......
......@@ -56,8 +56,8 @@ class SampleStepper
void set_iter_cmd(std::string s) { iter_cmd_ = s; }
void set_iter_cmd_period(int i) { iter_cmd_period_ = i; }
virtual EnergyFunctional& ef(void) {}
virtual ChargeDensity& cd(void) {}
virtual EnergyFunctional& ef(void) = 0;
virtual ChargeDensity& cd(void) = 0;
SampleStepper(Sample& s);
virtual ~SampleStepper(void);
......
......@@ -77,8 +77,10 @@ SlaterDet::SlaterDet(const Context& ctxt, D3vector kpoint) : ctxt_(ctxt),
////////////////////////////////////////////////////////////////////////////////
SlaterDet::SlaterDet(const SlaterDet& rhs) : ctxt_(rhs.context()),
basis_(new Basis(*(rhs.basis_))),
my_col_comm_(rhs.my_col_comm_), c_(rhs.c_){}
basis_(new Basis(*(rhs.basis_))), c_(rhs.c_)
{
MPI_Comm_dup(rhs.my_col_comm_,&my_col_comm_);
}
////////////////////////////////////////////////////////////////////////////////
SlaterDet::~SlaterDet()
......@@ -103,6 +105,7 @@ SlaterDet::~SlaterDet()
}
}
#endif
MPI_Comm_free(&my_col_comm_);
}
////////////////////////////////////////////////////////////////////////////////
......
//
// kpgen.C: generate a kpoint list for an arbitrary cell
// use: kpgen nx ny nz sx sy sz a0x a0y a0z a1x a1y a1z a2x a2y a2z
//
// where a0x a0y a0z = first basis vector of the unit cell
// a1x a1y a1z = second basis vector of the unit cell
// a2x a2y a2z = third basis vector of the unit cell
//
// nx,ny,nz: number of kpoints in each direction
// sx,sy,sz: shift in each direction (floating point)
// shift: 0: symmetric set: boundary points not included
// even-numbered sets do not include the gamma point
// odd-numbered sets include the gamma point
//
#include<iostream>
#include<fstream>
#include<iomanip>
#include<vector>
#include<cassert>
#include<cstdlib>
#include<list>
#include "D3vector.h"
using namespace std;
// ib_BZ: test if the vector k is in the BZ defined by b0,b1,b2
bool in_BZ(D3vector k, D3vector b0, D3vector b1, D3vector b2)
{
const double epsilon = 1.e-8;
D3vector g;
// check projection of kpoint k along all 26 reciprocal lattice vectors
// that are nearest g=0
bool in_bz = true;
for ( int i0 = -1; i0 <= 1; i0++ )
for ( int i1 = -1; i1 <= 1; i1++ )
for ( int i2 = -1; i2 <= 1; i2++ )
if ( !(i0 == 0 && i1 == 0 && i2 == 0) )
{
D3vector g = i0 * b0 + i1 * b1 + i2 * b2;
if ( k*g > 0.5 * g*g + epsilon )
in_bz = false;
}
return in_bz;
}
int main(int argc, char** argv)
{
cout << "# kpgen-1.0" << endl;
if ( argc != 16 )
{
cerr << " use: " << argv[0] << " nx ny nz shiftx shifty shiftz {cell}"
<< endl;
cerr << " shift==0: symmetric set, zone boundary not included" << endl;
return 1;
}
int nx = atoi(argv[1]);
int ny = atoi(argv[2]);
int nz = atoi(argv[3]);
double sx = atof(argv[4]);
double sy = atof(argv[5]);
double sz = atof(argv[6]);
if ( nx <= 0 || ny <= 0 || nz <=0 )
{
cerr << " use: " << argv[0] << " nx ny nz shiftx shifty shiftz {cell}"
<< endl;
cerr << " nx, ny, nz must be positive" << endl;
return 1;
}
if ( sx < 0 || sx > 1 ||
sy < 0 || sy > 1 ||
sz < 0 || sz > 1 )
{
cerr << " use: " << argv[0] << " nx ny nz shiftx shifty shiftz {cell}"
<< endl;
cerr << " shifts must be in [0,1]" << endl;
return 1;
}
D3vector a0(atof(argv[7]),atof(argv[8]),atof(argv[9]));
D3vector a1(atof(argv[10]),atof(argv[11]),atof(argv[12]));
D3vector a2(atof(argv[13]),atof(argv[14]),atof(argv[15]));
const double volume = a0 * ( a1 ^ a2 );
if ( volume == 0.0 )
{
cout << " cell volume is zero" << endl;
return 1;
}
double fac = 2.0 * M_PI / volume;
D3vector b0 = fac * a1 ^ a2;
D3vector b1 = fac * a2 ^ a0;
D3vector b2 = fac * a0 ^ a1;
list<vector<int> > kplist;
vector<int> kpint(4);
// scan volume enclosing the BZ
for ( int ii = -2; ii <= 2; ii++ )
for ( int jj = -2; jj <= 2; jj++ )
for ( int kk = -1; kk <= 2; kk++ )
for ( int i = 0; i < nx; i++ )
{
for ( int j = 0; j < ny; j++ )
{
for ( int k = 0; k < nz; k++ )
{
kpint[0] = ii*2*nx + 2*i-nx+1;
kpint[1] = jj*2*ny + 2*j-ny+1;
kpint[2] = kk*2*nz + 2*k-nz+1;
kpint[3] = 1;
D3vector k = kpint[0]/(2.0*nx) * b0 +
kpint[1]/(2.0*ny) * b1 +
kpint[2]/(2.0*nz) * b2;
if ( in_BZ(k,b0,b1,b2) )
kplist.push_back(kpint);
}
}
}
int total_weight = kplist.size();
#if 1
// remove -k
for (list<vector<int> >::iterator i = kplist.begin(); i != kplist.end(); i++)
{
// test if -k is in the set (and is not 0 0 0)
kpint[0] = (*i)[0];
kpint[1] = (*i)[1];
kpint[2] = (*i)[2];
kpint[3] = (*i)[3];
if ( kpint[0]*kpint[0]+kpint[1]*kpint[1]+kpint[2]*kpint[2] != 0 )
{
// look for -k in the rest of the list
for ( list<vector<int> >::iterator j = i; j != kplist.end(); j++ )
{
if ( (*j)[0]==-kpint[0] && (*j)[1]==-kpint[1] && (*j)[2]==-kpint[2] )
{
// transfer weight to (*i)
(*i)[3] += (*j)[3];
(*j)[3] = 0;
//cout << " erasing " << "(" << kpint[0] << ","
// << kpint[1] << "," << kpint[2] << ") == -("
// << (*j)[0] << "," << (*j)[1] << "," << (*j)[2] << ")" << endl;
}
}
}
}
#endif
#if 1
// remove equivalent points
for (list<vector<int> >::iterator i = kplist.begin(); i != kplist.end(); i++)
{
D3vector ki = (*i)[0]/(2.0*nx) * b0 +
(*i)[1]/(2.0*ny) * b1 +
(*i)[2]/(2.0*nz) * b2;
// look for equivalent points in the rest of the list
for ( list<vector<int> >::iterator j = i; j != kplist.end(); j++ )
{
if ( j != i )
{
D3vector kj = (*j)[0]/(2.0*nx) * b0 +
(*j)[1]/(2.0*ny) * b1 +
(*j)[2]/(2.0*nz) * b2;
if ( length(ki-kj) < 1.e-5 )
{
// transfer the weight of kj to ki
(*i)[3] += (*j)[3];
(*j)[3] = 0;
//cout << " erasing equivalent point " << "(" << (*j)[0] << ","
// << (*j)[1] << "," << (*j)[2] << ") == ("
// << (*i)[0] << "," << (*i)[1] << "," << (*i)[2] << ")" << endl;
}
}
}
}
#endif
// remove elements with zero weight
for (list<vector<int> >::iterator i = kplist.begin();
i != kplist.end(); /* nothing */ )
{
if ( (*i)[3] == 0 )
{
kplist.erase(i++);
}
else
{
i++;
}
}
#if 1
// make first index positive
for (list<vector<int> >::iterator i = kplist.begin(); i != kplist.end(); i++)
{
D3vector ki = (*i)[0]/(2.0*nx) * b0 +
(*i)[1]/(2.0*ny) * b1 +
(*i)[2]/(2.0*nz) * b2;
if ( ki.x < 0 )
{
(*i)[0] *= -1;
(*i)[1] *= -1;
(*i)[2] *= -1;
}
}
#endif
#if 1
// check that the sum of weights is one
int sum = 0;
for (list<vector<int> >::iterator i = kplist.begin(); i != kplist.end(); i++)
{
sum += (*i)[3];
}
assert(sum==total_weight);
#endif
#if 1
// output list
// traverse list backwards to have increasing indices
// kpoints are output in reciprocal lattice coordinates
cout.setf(ios::right,ios::adjustfield);
cout << "# nx,ny,nz: " << nx << " " << ny << " " << nz << endl;
cout << "# sx,sy,sz: " << sx << " " << sy << " " << sz << endl;
cout << "# a0: " << a0 << endl;
cout << "# a1: " << a1 << endl;
cout << "# a2: " << a2 << endl;
cout << "# b0/(2pi): " << b0/(2*M_PI) << endl;
cout << "# b1/(2pi): " << b1/(2*M_PI) << endl;
cout << "# b2/(2pi): " << b2/(2*M_PI) << endl;
cout << "# " << kplist.size() << " k-points" << endl;
cout << " kpoint delete 0 0 0" << endl;
for (list<vector<int> >::reverse_iterator i = kplist.rbegin();
i != kplist.rend(); i++)
{
cout.setf(ios::fixed,ios::floatfield);
cout << " kpoint add "
<< setprecision(10)
<< setw(13) << ((*i)[0]+sx)/(2.0*nx) << " "
<< setw(13) << ((*i)[1]+sy)/(2.0*ny) << " "
<< setw(13) << ((*i)[2]+sz)/(2.0*nz) << " ";
cout.setf(ios::scientific,ios::floatfield);
cout << setprecision(14)
<< setw(16) << (*i)[3]/((double) total_weight)
<< endl;
}
// output list in absolute coordinates for plot
ofstream plotfile("kpoint.plt");
for (list<vector<int> >::reverse_iterator i = kplist.rbegin();
i != kplist.rend(); i++)
{
D3vector k = (((*i)[0]+sx)/(2.0*nx)) * b0 / (2*M_PI) +
(((*i)[1]+sy)/(2.0*ny)) * b1 / (2*M_PI) +
(((*i)[2]+sz)/(2.0*nz)) * b2 / (2*M_PI);
plotfile.setf(ios::fixed,ios::floatfield);
plotfile << setprecision(8)
<< setw(13) << k.x << " "
<< setw(13) << k.y << " "
<< setw(13) << k.z << endl;
}
#endif
}
......@@ -17,7 +17,6 @@
// double 3-vectors
//
////////////////////////////////////////////////////////////////////////////////
// $Id: D3vector.h,v 1.8 2008-09-08 15:56:18 fgygi Exp $
#ifndef D3VECTOR_H
#define D3VECTOR_H
......@@ -46,6 +45,14 @@ class D3vector
else return z;
}
double operator[] (const int &i) const
{
assert(i>=0 && i <3);
if ( i == 0 ) return x;
else if ( i == 1 ) return y;
else return z;
}
bool operator==(const D3vector &rhs) const
{
return x == rhs.x && y == rhs.y && z == rhs.z;
......@@ -142,7 +149,7 @@ class D3vector
return sqrt( a.x * a.x + a.y * a.y + a.z * a.z );
}
friend double norm( const D3vector& a )
friend double norm2( const D3vector& a )
{
return a.x * a.x + a.y * a.y + a.z * a.z;
}
......
kpgen: kpgen.cpp D3vector.h
g++ -o $@ $^
g++ -o $@ kpgen.cpp
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