...
 
Commits (54)
......@@ -28,4 +28,4 @@
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.
Michael LaCount implemented the SCAN meta-GGA density functional.
......@@ -74,7 +74,8 @@ int walsh(int l, int n, int i)
}
////////////////////////////////////////////////////////////////////////////////
Bisection::Bisection(const SlaterDet& sd, int nlevels[3]) : ctxt_(sd.context())
Bisection::Bisection(const SlaterDet& sd, const int nlevels[3])
: ctxt_(sd.context())
{
// localization vectors are long int
assert(sizeof(long int) >= 4);
......
......@@ -69,7 +69,7 @@ class Bisection
public:
Bisection(const SlaterDet& sd, int nlevels[3]);
Bisection(const SlaterDet& sd, const int nlevels[3]);
void compute_transform(const SlaterDet& sd, int maxsweep, double tol);
void compute_localization(double epsilon);
void forward(SlaterDet& sd);
......
......@@ -124,7 +124,6 @@ void ChargeDensity::update_density(void)
{
assert(rhor.size() == wf_.nspin());
const double omega = vbasis_->cell().volume();
for ( int ispin = 0; ispin < wf_.nspin(); ispin++ )
{
assert(rhor[ispin].size() == vft_->np012loc() );
......@@ -209,6 +208,52 @@ void ChargeDensity::update_rhor(void)
}
////////////////////////////////////////////////////////////////////////////////
void ChargeDensity::update_taur(double* taur) const
{
memset( (void*)taur, 0, vft_->np012loc()*sizeof(double) );
tmap["update_taur"].start();
for ( int ispin = 0; ispin < wf_.nspin(); ispin++ )
{
for ( int ikp = 0; ikp < wf_.nkp(); ikp++ )
{
wf_.sd(ispin,ikp)->compute_tau(*ft_[ikp], wf_.weight(ikp), taur);
}
}
// sum along columns of spincontext
wf_.kpcontext()->dsum('r',vft_->np012loc(),1,&taur[0],vft_->np012loc());
tmap["update_taur"].stop();
// stop if computing taur with NLCCs
if ( !rhocore_r.empty() )
assert(!"ChargeDensity: Cannot compute taur with NLCCs");
}
////////////////////////////////////////////////////////////////////////////////
void ChargeDensity::update_taur(double* taur_up, double* taur_dn) const
{
memset( (void*)taur_up, 0, vft_->np012loc()*sizeof(double) );
memset( (void*)taur_dn, 0, vft_->np012loc()*sizeof(double) );
tmap["update_taur"].start();
for ( int ikp = 0; ikp < wf_.nkp(); ikp++ )
{
wf_.sd(0,ikp)->compute_tau(*ft_[ikp], wf_.weight(ikp), taur_up);
}
for ( int ikp = 0; ikp < wf_.nkp(); ikp++ )
{
wf_.sd(1,ikp)->compute_tau(*ft_[ikp], wf_.weight(ikp), taur_dn);
}
// sum along columns of spincontext
wf_.kpcontext()->dsum('r',vft_->np012loc(),1,&taur_up[0],vft_->np012loc());
wf_.kpcontext()->dsum('r',vft_->np012loc(),1,&taur_dn[0],vft_->np012loc());
tmap["update_taur"].stop();
// stop if computing taur with NLCCs
if ( !rhocore_r.empty() )
assert(!"ChargeDensity: Cannot compute taur with NLCCs");
}
////////////////////////////////////////////////////////////////////////////////
double ChargeDensity::total_charge(void) const
{
assert((wf_.nspin()==1)||(wf_.nspin()==2));
......
......@@ -59,6 +59,8 @@ class ChargeDensity
std::vector<std::complex<double> > rhocore_g;
void update_density(void);
void update_rhor(void);
void update_taur(double* taur) const;
void update_taur(double* taur_up, double* taur_dn) const;
void update_rhog(void);
const Context& context(void) const { return ctxt_; }
......
......@@ -387,7 +387,7 @@ void ExchangeOperator::cell_moved(void)
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
double ExchangeOperator::compute_exchange_for_general_case_( Sample* s,
double ExchangeOperator::compute_exchange_for_general_case_(const Sample* s,
Wavefunction* dwf, bool compute_stress)
{
Timer tm;
......
......@@ -28,7 +28,7 @@ class ExchangeOperator
{
private:
Sample& s_;
const Sample& s_;
// exchange energy
double eex_;
// constant of support function for exchange integration
......@@ -43,7 +43,7 @@ class ExchangeOperator
// copy of wf
Wavefunction wfc_;
double compute_exchange_for_general_case_(Sample* s, Wavefunction* dwf,
double compute_exchange_for_general_case_(const Sample* s, Wavefunction* dwf,
bool compute_stress);
double compute_exchange_at_gamma_(const Wavefunction &wf, Wavefunction* dwf,
bool compute_stress);
......
......@@ -48,16 +48,10 @@ class HSEFunctional : public XCFunctional
HSEFunctional(const std::vector<std::vector<double> > &rhoe);
// HSE's local part is derived from PBE
bool isGGA() const
{
return true;
}
bool isGGA() const { return true; }
// return the name of the functional
std::string name() const
{
return "HSE";
}
std::string name() const { return "HSE"; }
void setxc(void);
};
......
......@@ -62,7 +62,6 @@ class LDAFunctional : public XCFunctional
}
};
bool isGGA() const { return false; };
std::string name() const { return "LDA"; };
void setxc(void);
};
......
......@@ -30,8 +30,8 @@ OBJECTS=qb.o AtomSet.o Atom.o Species.o \
SpeciesCmd.o SpeciesReader.o SpeciesHandler.o \
XCOperator.o ExchangeOperator.o Bisection.o \
XCPotential.o LDAFunctional.o VWNFunctional.o \
PBEFunctional.o BLYPFunctional.o B3LYPFunctional.o \
BHandHLYPFunctional.o \
PBEFunctional.o BLYPFunctional.o SCANFunctional.o \
B3LYPFunctional.o BHandHLYPFunctional.o \
ExponentialIntegral.o HSEFunctional.o \
RSHFunctional.o \
NonLocalPotential.o SampleReader.o StructuredDocumentHandler.o \
......@@ -491,6 +491,17 @@ JDWavefunctionStepper.o: D3vector.h UnitCell.h
KpointCmd.o: UserInterface.h D3vector.h Sample.h AtomSet.h Context.h blacs.h
KpointCmd.o: Atom.h UnitCell.h D3tensor.h blas.h ConstraintSet.h
KpointCmd.o: ExtForceSet.h Wavefunction.h Control.h
LBFGS_IonicStepper.o: LBFGS_IonicStepper.h IonicStepper.h Sample.h AtomSet.h
LBFGS_IonicStepper.o: Context.h blacs.h Atom.h D3vector.h UnitCell.h
LBFGS_IonicStepper.o: D3tensor.h blas.h ConstraintSet.h ExtForceSet.h
LBFGS_IonicStepper.o: Wavefunction.h Control.h Species.h
LBFGS_IonicStepper.o: IonicStepper.h Sample.h AtomSet.h Context.h blacs.h
LBFGS_IonicStepper.o: Atom.h D3vector.h UnitCell.h D3tensor.h blas.h
LBFGS_IonicStepper.o: ConstraintSet.h ExtForceSet.h Wavefunction.h Control.h
LBFGS_IonicStepper.o: Species.h
LBFGS_hessian.o: Sample.h AtomSet.h Context.h blacs.h Atom.h D3vector.h
LBFGS_hessian.o: UnitCell.h D3tensor.h blas.h ConstraintSet.h ExtForceSet.h
LBFGS_hessian.o: Wavefunction.h Control.h
LDAFunctional.o: LDAFunctional.h XCFunctional.h
LDAFunctional.o: XCFunctional.h
LineMinimizer.o: LineMinimizer.h
......@@ -638,6 +649,8 @@ RunCmd.o: Context.h blacs.h D3vector.h Wavefunction.h UnitCell.h SlaterDet.h
RunCmd.o: Basis.h Sample.h AtomSet.h Atom.h D3tensor.h blas.h ConstraintSet.h
RunCmd.o: ExtForceSet.h Control.h ChargeDensity.h CPSampleStepper.h
RunCmd.o: UserInterface.h
SCANFunctional.o: SCANFunctional.h XCFunctional.h
SCANFunctional.o: XCFunctional.h
SDAIonicStepper.o: SDAIonicStepper.h IonicStepper.h Sample.h AtomSet.h
SDAIonicStepper.o: Context.h blacs.h Atom.h D3vector.h UnitCell.h D3tensor.h
SDAIonicStepper.o: blas.h ConstraintSet.h ExtForceSet.h Wavefunction.h
......@@ -793,15 +806,15 @@ XCOperator.o: ExtForceSet.h Wavefunction.h Control.h ChargeDensity.h Timer.h
XCOperator.o: XCPotential.h ExchangeOperator.h SlaterDet.h Basis.h Matrix.h
XCOperator.o: FourierTransform.h HSEFunctional.h XCFunctional.h
XCOperator.o: RSHFunctional.h
XCOperator.o: Sample.h AtomSet.h Context.h blacs.h Atom.h D3vector.h
XCOperator.o: UnitCell.h D3tensor.h blas.h ConstraintSet.h ExtForceSet.h
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
XCPotential.o: BHandHLYPFunctional.h Basis.h UnitCell.h FourierTransform.h
XCPotential.o: blas.h
XCPotential.o: SCANFunctional.h HSEFunctional.h RSHFunctional.h
XCPotential.o: B3LYPFunctional.h BHandHLYPFunctional.h Sample.h AtomSet.h
XCPotential.o: Atom.h UnitCell.h D3tensor.h blas.h ConstraintSet.h
XCPotential.o: ExtForceSet.h Wavefunction.h SlaterDet.h Basis.h Matrix.h
XCPotential.o: StructureFactor.h XCOperator.h NonLocalPotential.h
XCPotential.o: ConfinementPotential.h ElectricEnthalpy.h FourierTransform.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
......
......@@ -74,16 +74,10 @@ class RSHFunctional : public XCFunctional
RSHFunctional(const std::vector<std::vector<double> > &rhoe,
double alpha_RSH, double beta_RSH, double mu_RSH);
bool isGGA() const
{
return true;
}
bool isGGA() const { return true; }
// return the name of the functional
std::string name() const
{
return "RSH";
}
std::string name() const { return "RSH"; }
void setxc(void);
};
......
////////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2008 The Regents of the University of California
//
// This file is part of Qbox
//
// 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/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// SCANFunctional.h
//
////////////////////////////////////////////////////////////////////////////////
#ifndef SCANFUNCTIONAL_H
#define SCANFUNCTIONAL_H
#include "XCFunctional.h"
#include <vector>
class SCANFunctional : public XCFunctional
{
double x_coeff_, c_coeff_;
std::vector<double> exc_, exc_up_, exc_dn_;
std::vector<double> vxc1_, vxc1_up_, vxc1_dn_,
vxc2_, vxc2_upup_, vxc2_updn_, vxc2_dnup_, vxc2_dndn_,
vxc3_,vxc3_up_,vxc3_dn_;
std::vector<double> grad_rho_[3], grad_rho_up_[3], grad_rho_dn_[3];
std::vector<double> tau_, tau_up_, tau_dn_;
void gPW92(double alpha, double beta0, double beta1, double beta2,
double beta3, double beta4, double rtrs, double *gg, double *dgdrs);
void excSCAN(double rho, double grad, double tau, double *exc,
double *vxc1, double *vxc2, double *vxc3);
void excSCAN_sp(double rho_up, double rho_dn, double grad_up, double grad_dn,
double grad, double tau_up, double tau_dn, double *exc_up, double *exc_dn,
double *vxc1_up, double *vxc1_dn, double *vxc2_upup, double *vxc2_dndn,
double *vxc2_updn, double *vxc2_dnup, double *vxc3_up, double *vxc3_dn);
public:
// constructor with variable coefficients for exchange and correlation
// with default values 1.0
SCANFunctional(const std::vector<std::vector<double> > &rhoe,
double x_coeff=1.0, double c_coeff=1.0);
bool isGGA() const { return true; };
bool isMeta() const { return true; };
std::string name() const { return "SCAN"; };
void setxc(void);
};
#endif
......@@ -358,6 +358,115 @@ void SlaterDet::compute_density(FourierTransform& ft,
}
////////////////////////////////////////////////////////////////////////////////
void SlaterDet::compute_tau(FourierTransform& ft,
double weight, double* taur) const
{
// compute tau of the states residing on my column of ctxt_
assert(occ_.size() == c_.n());
vector<complex<double> > tmp(ft.np012loc());
const int ngwloc = basis_->localsize();
const int mloc = c_.mloc();
assert(basis_->cell().volume() > 0.0);
const double omega_inv = 1.0 / basis_->cell().volume();
const int np012loc = ft.np012loc();
if ( basis_->real() )
{
vector<complex<double> > taug0(ngwloc), taug1(ngwloc);
// transform two states at a time
for ( int n = 0; n < nstloc()-1; n++, n++ )
{
const int nn = ctxt_.mycol() * c_.nb() + n;
// next line: factor 0.5 from definition of tau
const double fac1 = 0.5 * weight * omega_inv * occ_[nn];
const double fac2 = 0.5 * weight * omega_inv * occ_[nn+1];
const complex<double>* p = c_.cvalptr();
if ( fac1 + fac2 > 0.0 )
{
for ( int j = 0; j < 3; j++ )
{
const double *const gxj = basis_->gx_ptr(j);
for ( int ig = 0; ig < ngwloc; ig++ )
{
// i*G_j*c(G)
taug0[ig] = complex<double>(0.0,gxj[ig]) * p[ig+n*mloc];
taug1[ig] = complex<double>(0.0,gxj[ig]) * p[ig+(n+1)*mloc];
}
ft.backward(&taug0[0],&taug1[0],&tmp[0]);
const double* gpsi = (double*) &tmp[0];
int ii = 0;
for ( int i = 0; i < np012loc; i++ )
{
const double gpsi1 = gpsi[ii];
const double gpsi2 = gpsi[ii+1];
taur[i] += fac1 * gpsi1 * gpsi1 + fac2 * gpsi2 * gpsi2;
ii++; ii++;
}
}
}
}
if ( nstloc() % 2 != 0 )
{
const int n = nstloc()-1;
// global n index
const int nn = ctxt_.mycol() * c_.nb() + n;
const double fac1 = 0.5 * weight * omega_inv * occ_[nn];
const complex<double>* p = c_.cvalptr();
if ( fac1 > 0.0 )
{
for ( int j = 0; j < 3; j++ )
{
const double *const gxj = basis_->gx_ptr(j);
for ( int ig = 0; ig < ngwloc; ig++ )
{
// i*G_j*c(G)
taug1[ig] = complex<double>(0.0,gxj[ig]) * p[ig+n*mloc];
}
ft.backward(&taug1[0],&tmp[0]);
const double* gpsi = (double*) &tmp[0];
int ii = 0;
for ( int i = 0; i < np012loc; i++ )
{
const double gpsi1 = gpsi[ii];
taur[i] += fac1 * gpsi1 * gpsi1;
ii++; ii++;
}
}
}
}
}
else
{
vector<complex<double> > taug(ngwloc);
for ( int n = 0; n < nstloc(); n++ )
{
const int nn = ctxt_.mycol() * c_.nb() + n;
// next line: factor 0.5 from definition of tau
const double fac1 = 0.5 * weight * omega_inv * occ_[nn];
const complex<double>* p = c_.cvalptr();
if ( fac1 > 0.0 )
{
for ( int j = 0; j < 3; j++ )
{
const double *const kpgxj = basis_->kpgx_ptr(j);
for ( int ig = 0; ig < ngwloc; ig++ )
{
// i*(k+G)_j*c(G)
taug[ig] = complex<double>(0.0,kpgxj[ig]) * p[ig+n*mloc];
}
ft.backward(&taug[0],&tmp[0]);
for ( int i = 0; i < np012loc; i++ )
{
taur[i] += fac1 * norm(tmp[i]);
}
}
}
}
}
}
////////////////////////////////////////////////////////////////////////////////
void SlaterDet::rs_mul_add(FourierTransform& ft,
const double* v, SlaterDet& sdp) const
{
......
......@@ -67,6 +67,7 @@ class SlaterDet
void resize(const UnitCell& cell, const UnitCell& refcell,
double ecut, int nst);
void compute_density(FourierTransform& ft, double weight, double* rho) const;
void compute_tau(FourierTransform& ft, double weight, double* taur) const;
void rs_mul_add(FourierTransform& ft, const double* v, SlaterDet& sdp) const;
void randomize(double amplitude);
void cleanup(void);
......
......@@ -74,7 +74,6 @@ class VWNFunctional : public XCFunctional
static void ecvwn_sp(double roe_up, double roe_dn,
double &ec, double &vc_up, double &vc_dn);
bool isGGA() const { return false; };
std::string name() const { return "VWN"; };
void setxc(void);
};
......
......@@ -65,14 +65,20 @@ class XCFunctional
const double *rho, *rho_up, *rho_dn;
double *grad_rho[3], *grad_rho_up[3], *grad_rho_dn[3];
double *tau, *tau_up, *tau_dn;
double *exc, *exc_up, *exc_dn;
double *vxc1, *vxc1_up, *vxc1_dn;
double *vxc2, *vxc2_upup, *vxc2_dndn, *vxc2_updn, *vxc2_dnup;
double *vxc3,*vxc3_up,*vxc3_dn;
virtual bool isGGA(void) const = 0;
// default functional is not GGA, not Meta
// GGA functionals must override isGGA()
// meta functionals must override isMeta()
virtual bool isGGA(void) const { return false; }
virtual bool isMeta(void) const { return false; }
virtual std::string name(void) const = 0;
int np(void) const { return _np; };
int nspin(void) const { return _nspin; };
int np(void) const { return _np; }
int nspin(void) const { return _nspin; }
XCFunctional()
{
......@@ -80,9 +86,11 @@ class XCFunctional
grad_rho[0] = grad_rho[1] = grad_rho[2] = 0;
grad_rho_up[0] = grad_rho_up[1] = grad_rho_up[2] = 0;
grad_rho_dn[0] = grad_rho_dn[1] = grad_rho_dn[2] = 0;
tau = tau_up = tau_dn = 0;
exc = exc_up = exc_dn = 0;
vxc1 = vxc1_up = vxc1_dn = 0;
vxc2 = vxc2_upup = vxc2_dndn = vxc2_updn = vxc2_dnup = 0;
vxc3 = 0;
}
// virtual destructor needed to ensure proper deallocation
......
......@@ -16,6 +16,7 @@
//
////////////////////////////////////////////////////////////////////////////////
#include "XCOperator.h"
#include "Sample.h"
#include "ChargeDensity.h"
#include "XCPotential.h"
#include "ExchangeOperator.h"
......@@ -24,13 +25,15 @@
using namespace std;
////////////////////////////////////////////////////////////////////////////////
XCOperator::XCOperator(Sample& s, const ChargeDensity& cd) :cd_(cd)
XCOperator::XCOperator(Sample& s, const ChargeDensity& cd) :cd_(cd), s_(s)
{
// set initial values
xcp_ = 0;
xop_ = 0;
exc_ = 0.0 ;
dxc_ = 0.0 ;
hasHF_ = false;
hasMeta_ = false;
sigma_exc_.resize(6);
......@@ -43,7 +46,7 @@ XCOperator::XCOperator(Sample& s, const ChargeDensity& cd) :cd_(cd)
( functional_name == "BLYP" ) )
{
// create only an xc potential
xcp_ = new XCPotential(cd, functional_name, s.ctrl);
xcp_ = new XCPotential(cd, functional_name, s_);
hasPotential_ = true;
hasGGA_ = xcp_->isGGA();
hasHF_ = false;
......@@ -59,7 +62,7 @@ XCOperator::XCOperator(Sample& s, const ChargeDensity& cd) :cd_(cd)
else if ( functional_name == "PBE0" )
{
// create an exchange potential
xcp_ = new XCPotential(cd, functional_name, s.ctrl);
xcp_ = new XCPotential(cd, functional_name, s_);
// create the exchange operator with mixing coeff=0.25
xop_ = new ExchangeOperator(s, s.ctrl.alpha_PBE0, s.ctrl.alpha_PBE0, 0.0);
......@@ -70,7 +73,7 @@ XCOperator::XCOperator(Sample& s, const ChargeDensity& cd) :cd_(cd)
else if ( functional_name == "HSE" )
{
// create an exchange potential
xcp_ = new XCPotential(cd, functional_name, s.ctrl);
xcp_ = new XCPotential(cd, functional_name, s_);
// create the exchange operator with mixing coeff=0.25
xop_ = new ExchangeOperator(s, 0.0, 0.25, 0.11);
......@@ -81,7 +84,7 @@ XCOperator::XCOperator(Sample& s, const ChargeDensity& cd) :cd_(cd)
else if ( functional_name == "RSH" )
{
// create an exchange potential
xcp_ = new XCPotential(cd, functional_name, s.ctrl);
xcp_ = new XCPotential(cd, functional_name, s_);
// create the exchange operator with mixing coeff=beta_RSH
xop_ = new ExchangeOperator(s, s.ctrl.alpha_RSH, s.ctrl.beta_RSH,
......@@ -93,7 +96,7 @@ XCOperator::XCOperator(Sample& s, const ChargeDensity& cd) :cd_(cd)
else if ( functional_name == "B3LYP" )
{
// create an exchange potential
xcp_ = new XCPotential(cd, functional_name, s.ctrl);
xcp_ = new XCPotential(cd, functional_name, s_);
// create the exchange operator with mixing coeff=0.20
xop_ = new ExchangeOperator(s, 0.20, 0.20, 0.0);
......@@ -104,7 +107,7 @@ XCOperator::XCOperator(Sample& s, const ChargeDensity& cd) :cd_(cd)
else if ( functional_name == "BHandHLYP" )
{
// create an exchange potential
xcp_ = new XCPotential(cd, functional_name, s.ctrl);
xcp_ = new XCPotential(cd, functional_name, s_);
// create the exchange operator with mixing coeff=0.50
xop_ = new ExchangeOperator(s, 0.50, 0.50, 0.0);
......@@ -112,6 +115,14 @@ XCOperator::XCOperator(Sample& s, const ChargeDensity& cd) :cd_(cd)
hasGGA_ = xcp_->isGGA();
hasHF_ = true;
}
else if ( functional_name == "SCAN" )
{
// create an exchange potential
xcp_ = new XCPotential(cd, functional_name, s_);
hasPotential_ = true;
hasGGA_ = xcp_->isGGA();
hasMeta_ = true;
}
else
{
throw XCOperatorException(
......@@ -133,7 +144,7 @@ void XCOperator::update(std::vector<std::vector<double> >& vr, bool compute_stre
// compute vxc potential and energy
if ( hasPotential_ )
{
// update LDA/GGA xc potential
// update LDA/GGA/MetaGGA xc potential
xcp_->update( vr );
// LDA/GGA exchange energy
......@@ -165,6 +176,8 @@ void XCOperator::apply_self_energy(Wavefunction &dwf)
{
if ( hasHF() )
xop_->apply_operator(dwf);
if ( hasMeta() )
xcp_->apply_meta_operator(dwf);
}
////////////////////////////////////////////////////////////////////////////////
......
......@@ -19,9 +19,12 @@
#ifndef XCOPERATOR_H
#define XCOPERATOR_H
#include "Sample.h"
#include <valarray>
#include <vector>
#include <string>
class Sample;
class Wavefunction;
class ChargeDensity;
class XCPotential;
class ExchangeOperator;
......@@ -32,6 +35,7 @@ class XCOperator
XCPotential* xcp_;
ExchangeOperator* xop_;
const Sample& s_;
const ChargeDensity& cd_;
double exc_; // XC energy: includes local and HF terms
double dxc_;
......@@ -41,6 +45,7 @@ class XCOperator
bool hasPotential_;
bool hasGGA_;
bool hasHF_;
bool hasMeta_;
public:
......@@ -58,6 +63,7 @@ class XCOperator
bool hasGGA(void) { return hasGGA_; };
bool hasHF(void) { return hasHF_; };
bool hasMeta(void) { return hasMeta_; };
void update(std::vector<std::vector<double> >& vr, bool compute_stress);
void apply_self_energy(Wavefunction &dwf);
......
......@@ -19,7 +19,6 @@
#ifndef XCPOTENTIAL_H
#define XCPOTENTIAL_H
#include "Control.h"
#include "ChargeDensity.h"
#include <string>
#include <vector>
......@@ -29,11 +28,14 @@
class Basis;
class FourierTransform;
class XCFunctional;
class Sample;
class Wavefunction;
class XCPotential
{
private:
const Sample& s_;
const ChargeDensity& cd_;
XCFunctional* xcf_;
......@@ -58,11 +60,13 @@ class XCPotential
const XCFunctional* xcf() { return xcf_; }
bool isGGA(void);
bool isMeta(void);
XCPotential(const ChargeDensity& cd, const std::string functional_name,
const Control& ctrl);
const Sample& s);
~XCPotential();
void update(std::vector<std::vector<double> >& vr);
void compute_stress(std::valarray<double>& sigma_exc);
void apply_meta_operator(Wavefunction& dwf);
double exc(void) { return exc_; }
double dxc(void) { return dxc_; }
};
......
......@@ -48,6 +48,7 @@ class Xc : public Var
v == "VWN" ||
v == "PBE" ||
v == "BLYP" ||
v == "SCAN" ||
v == "HF" ||
v == "PBE0" ||
v == "HSE" ||
......@@ -56,8 +57,8 @@ class Xc : public Var
v == "BHandHLYP" ) )
{
if ( ui->onpe0() )
cout << " xc must be LDA, VWN, PBE, BLYP, "
<< "HF, PBE0, HSE, RSH, B3LYP or BHandHLYP" << endl;
cout << " xc must be LDA, VWN, PBE, BLYP, SCAN,\n"
<< " HF, PBE0, HSE, RSH, B3LYP or BHandHLYP" << endl;
return 1;
}
......
--------------------------------------------------------------------------------
rel1_68_0
--------------------------------------------------------------------------------
05cd227 Merge branch 'scan' into develop
5c8ff65 Make Sample const in ExchangeOperator
15649b6 Make nlevels const in Bisection
8afc2a0 Add qbox_species_temp.sh utility script
4f46485 Add scripts to extract position,velocity,force
--------------------------------------------------------------------------------
rel1_67_4
--------------------------------------------------------------------------------
6e23b88 Merge branch 'develop'
......
......@@ -19,5 +19,5 @@
#include "release.h"
std::string release(void)
{
return std::string("1.67.4");
return std::string("1.68.0");
}
#!/usr/bin/python
# qbox_force.py
# extract force of an atom from Qbox output
# use: qbox_force.py atom_name file.r
import xml.sax
import sys
import math
if len(sys.argv) != 3:
print "use: ",sys.argv[0]," atom_name file.r"
sys.exit()
# Qbox output handler to extract and process <atomset>
class QboxOutputHandler(xml.sax.handler.ContentHandler):
def __init__(self):
self.readForce = 0
self.buffer = ""
def startElement(self, name, attributes):
if name == "atomset":
self.buffer = ""
elif name == "atom":
if attributes["name"] == atom_name:
self.readForce = 1
elif name == "force":
self.buffer = ""
def characters(self, data):
if self.readForce:
self.buffer += data
def endElement(self, name):
if name == "force":
if self.readForce:
force = self.buffer.split()
fx = float(force[0])
fy = float(force[1])
fz = float(force[2])
print '%.8f'%fx,'%.8f'%fy,'%.8f'%fz
self.readForce = 0
atom_name = sys.argv[1]
filename = sys.argv[2]
parser = xml.sax.make_parser()
handler = QboxOutputHandler()
parser.setContentHandler(handler)
parser.parse(sys.argv[2])
#!/usr/bin/python
# qbox_position.py
# extract position of an atom from Qbox output
# use: qbox_position.py atom_name file.r
import xml.sax
import sys
import math
if len(sys.argv) != 3:
print "use: ",sys.argv[0]," atom_name file.r"
sys.exit()
# Qbox output handler to extract and process <atomset>
class QboxOutputHandler(xml.sax.handler.ContentHandler):
def __init__(self):
self.readPos = 0
self.buffer = ""
def startElement(self, name, attributes):
if name == "atomset":
self.buffer = ""
elif name == "atom":
if attributes["name"] == atom_name:
self.readPos = 1
elif name == "position":
self.buffer = ""
def characters(self, data):
if self.readPos:
self.buffer += data
def endElement(self, name):
if name == "position":
if self.readPos:
pos = self.buffer.split()
rx = float(pos[0])
ry = float(pos[1])
rz = float(pos[2])
print '%.8f'%rx,'%.8f'%ry,'%.8f'%rz
self.readPos = 0
atom_name = sys.argv[1]
filename = sys.argv[2]
parser = xml.sax.make_parser()
handler = QboxOutputHandler()
parser.setContentHandler(handler)
parser.parse(sys.argv[2])
#!/bin/bash
#
# qbox_species_temp.sh
#
# compute the temperature of each atom of a given species in an MD run
#
# use: qbox_species_temp.sh species_name file [file ...]
#
if (($#<2))
then echo " use: qbox_species_temp.sh species file [file ...]"
exit 1
fi
species=$1
shift 1
echo "# temperature of atoms of species: $species, file: ${*}"
# extract mass of species
xpath='//species[@name="'$species'"]/mass'
elem=$(basename $xpath)
mass=$(xmllint --xpath "$xpath" ${*} |sed "s/<$elem>//g" |sed "s/<\/$elem>/>/g" |tr '>' '\n' | sed "/^$/d")
# Boltzmann constant in hartree/K
kB=3.1667907e-06
# nucleon mass
nmass=1822.89
# conversion factor from velocity^2 (a.u.) to K
# fac = (2/3) * (1/2) * nmass * mass * (1/kB)
# extract velocity of each atom of the given species
xpath='//atomset/atom[@species="'$species'"]/velocity'
elem=$(basename $xpath)
xmllint --xpath "$xpath" ${*} |sed "s/<$elem>//g" |sed "s/<\/$elem>/>/g" |tr '>' '\n' | sed "/^$/d" | \
awk -v kB=$kB -v nmass=$nmass -v mass=$mass \
'{v2=$1*$1+$2*$2+$3*$3; print (2.0/3.0)*0.5*nmass*mass*v2/kB}' -
#!/usr/bin/python
# qbox_velocity.py
# extract velocity of an atom from Qbox output
# use: qbox_velocity.py atom_name file.r
import xml.sax
import sys
import math
if len(sys.argv) != 3:
print "use: ",sys.argv[0]," atom_name file.r"
sys.exit()
# Qbox output handler to extract and process <atomset>
class QboxOutputHandler(xml.sax.handler.ContentHandler):
def __init__(self):
self.readVelocity = 0
self.buffer = ""
def startElement(self, name, attributes):
if name == "atomset":
self.buffer = ""
elif name == "atom":
if attributes["name"] == atom_name:
self.readVelocity = 1
elif name == "velocity":
self.buffer = ""
def characters(self, data):
if self.readVelocity:
self.buffer += data
def endElement(self, name):
if name == "velocity":
if self.readVelocity:
velocity = self.buffer.split()
vx = float(velocity[0])
vy = float(velocity[1])
vz = float(velocity[2])
print '%.8f'%vx,'%.8f'%vy,'%.8f'%vz
self.readVelocity = 0
atom_name = sys.argv[1]
filename = sys.argv[2]
parser = xml.sax.make_parser()
handler = QboxOutputHandler()
parser.setContentHandler(handler)
parser.parse(sys.argv[2])