Merge branch 'hse-dev' into develop

src/ExchangeOperator.h

@@ -19,6 +19,7 @@
 #include "Sample.h"
 #include "SlaterDet.h"
 #include "FourierTransform.h"
+#include "InteractionPotential.h"
 
 #ifndef EXCHANGEOPERATOR_H
 #define EXCHANGEOPERATOR_H
@@ -74,10 +75,8 @@ class ExchangeOperator
   // G vectors
   valarray<double> qpG21_;
   valarray<double> qpG22_;
-  valarray<double> qpG2i1_;
-  valarray<double> qpG2i2_;
-  valarray<double> G2_;
-  valarray<double> G2i_;
+  valarray<double> int_pot1_;
+  valarray<double> int_pot2_;
 
   // numbers of states
   int nLocalStates_;
@@ -173,10 +172,17 @@ class ExchangeOperator
   vector<DoubleMatrix*> uc_;
   vector<long int> localization_;
 
+  // Fourier transform of interaction potential
+  const InteractionPotential interaction_potential_;
+
+  // coulomb potential
+  bool coulomb_;
+
   public:
 
   // constructor
-  ExchangeOperator(Sample& s_, double HFCoeff);
+  ExchangeOperator(Sample& s_, double HFCoeff,
+    const InteractionPotential& interaction_potential = InteractionPotential());
 
   // destructor
   ~ExchangeOperator();

src/ExponentialIntegral.C

+////////////////////////////////////////////////////////////////////////////////
+//
+// 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/>.
+//
+////////////////////////////////////////////////////////////////////////////////
+//
+// ExponentialIntegral.C
+//
+////////////////////////////////////////////////////////////////////////////////
+//
+// Calculate exponential integral using the algorithm of
+// Tseng, Lee, Journal of Hydrology, 205 (1998) 38-51
+//
+////////////////////////////////////////////////////////////////////////////////
+
+#include "ExponentialIntegral.h"
+#include <cmath>
+#include <algorithm>
+#include <vector>
+
+using namespace std;
+
+namespace util
+{
+
+// Calculate the exponential integral E_1(x):
+//
+//        inf
+//          /     -t
+//         |     e
+// E (x) = | dt -----
+//  1      |      t
+//        /
+//         x
+//
+// Input:  x - position at which exponential integral is evaluated (x > 0)
+// Return: E1(x)
+double E1(const double x)
+{
+
+  if (x < series_cutoff)
+  {
+    // use series expansion
+    return E1_series(x);
+  }
+  else
+  {
+    // use gauss_laguerre expression
+    return exp(-x) * gauss_laguerre(x);
+  }
+}
+
+// Series expansion of the exponential integral
+//
+//                          n_cut
+//                          -----     n  n
+//                           \    (-1)  x
+// E (x) = -gamma - ln(x) -   )   --------
+//  1                        /     n * n!
+//                          -----
+//                          n = 1
+//
+// where gamma is the Euler constant.
+// n_cut is set to 25
+// Input:  x - position at which exponential integral is evaluated (x > 0)
+// Return: approximation by series expansion for E_1(x)
+double E1_series(const double x)
+{
+  // Euler constant
+  const double EULER_GAMMA = 0.57721566490153286060651209008241;
+
+  // Cutoff for series expansion
+  const int itmax = 25;
+
+  // initialize summation result
+  double res = 0.0;
+
+  // perform the summation
+  for (int it = itmax; it > 1; it--)
+  {
+    // calculate 1/n
+    const double fact = 1.0 / it;
+    // add next term of summation
+    res = x * fact * (fact - res);
+  }
+
+  // add everything up
+  return -EULER_GAMMA - log(x) + x * (1.0 - res);
+}
+
+// The Gauss Laguerre expansion of the exponential integral can be written as
+//
+//             N
+// E (x0)    -----     a
+//  1         \         n
+// ------ =    )   ---------
+//   -x0      /     x  + x0
+//  e        -----   n
+//            n=1
+//
+// where the a_n and x_n are determined by least quadrature (see reference)
+// Input: x0 - point at which Gaussian Laguerre quadrature is calculated
+// Return: E_1(x0) / exp(-x0) in this approximation
+double gauss_laguerre(const double x0)
+{
+  // initialize constants a_n and x_n
+  const double size = 15;
+  const double a[] = { 0.2182348859400869e+00, 0.3422101779228833e+00,
+      0.2630275779416801e+00, 0.1264258181059305e+00, 0.4020686492100091e-01,
+      0.8563877803611838e-02, 0.1212436147214252e-02, 0.1116743923442519e-03,
+      0.6459926762022901e-05, 0.2226316907096273e-06, 0.4227430384979365e-08,
+      0.3921897267041089e-10, 0.1456515264073126e-12, 0.1483027051113301e-15,
+      0.1600594906211133e-19 };
+  const double x[] = { 0.9330781201728180e-01, 0.4926917403018839e+00,
+      0.1215595412070949e+01, 0.2269949526203743e+01, 0.3667622721751437e+01,
+      0.5425336627413553e+01, 0.7565916226613068e+01, 0.1012022856801911e+02,
+      0.1313028248217572e+02, 0.1665440770832996e+02, 0.2077647889944877e+02,
+      0.2562389422672878e+02, 0.3140751916975394e+02, 0.3853068330648601e+02,
+      0.4802608557268579e+02 };
+
+  // initialize
+  double res = 0.0;
+
+  // evaluate a_n / ( x_n + x0 )
+  for ( int i = 0; i < size; i++ )
+  {
+    res += a[i] / (x[i] + x0);
+  }
+
+  return res;
+}
+
+}

src/ExponentialIntegral.h

+////////////////////////////////////////////////////////////////////////////////
+//
+// 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/>.
+//
+////////////////////////////////////////////////////////////////////////////////
+//
+// ExponentialIntegral.h
+//
+////////////////////////////////////////////////////////////////////////////////
+//
+// Calculate exponential integral using the algorithm of
+// Tseng, Lee, Journal of Hydrology, 205 (1998) 38-51
+//
+////////////////////////////////////////////////////////////////////////////////
+#ifndef EXPONENTIAL_INTEGRAL_H
+#define EXPONENTIAL_INTEGRAL_H
+
+namespace util {
+
+// constant at which we switch from series to gauss laguerre expansion
+static const double series_cutoff = 4.0;
+
+// Calculate the exponential integral E_1(x):
+//
+//        inf
+//          /     -t
+//         |     e
+// E (x) = | dt -----
+//  1      |      t
+//        /
+//         x
+//
+// Input:  x - position at which exponential integral is evaluated (x > 0)
+// Return: E1(x)
+double E1(const double x);
+
+// Series expansion of the exponential integral
+//
+//                          n_cut
+//                          -----     n  n
+//                           \    (-1)  x
+// E (x) = -gamma - ln(x) -   )   --------
+//  1                        /     n * n!
+//                          -----
+//                          n = 1
+//
+// where gamma is the Euler constant.
+// n_cut is set to 25
+// Input:  x - position at which exponential integral is evaluated (x > 0)
+// Return: approximation by series expansion for E_1(x)
+double E1_series(const double x);
+
+// The Gauss Laguerre expansion of the exponential integral can be written as
+//
+//             N
+// E (x0)    -----     a
+//  1         \         n
+// ------ =    )   ---------
+//   -x0      /     x  + x0
+//  e        -----   n
+//            n=1
+//
+// where the a_n and x_n are determined by least quadrature (see reference)
+// Input: x0 - point at which Gaussian Laguerre quadrature is calculated
+// Return: E_1(x0) / exp(-x0) in this approximation
+double gauss_laguerre(const double x0);
+
+}
+
+#endif

src/HSEFunctional.h

+////////////////////////////////////////////////////////////////////////////////
+//
+// 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/>.
+//
+////////////////////////////////////////////////////////////////////////////////
+//
+// HSEFunctional.h
+//
+////////////////////////////////////////////////////////////////////////////////
+//
+// Implements the HSE hybrid functional
+// Heyd et al.,   J. Chem. Phys. 118, 8207 (2003)
+// Heyd et al.,   J. Chem. Phys. 120, 7274 (2004)
+// Krukau et al., J. Chem. Phys. 125, 224106 (2006)
+//
+////////////////////////////////////////////////////////////////////////////////
+
+#ifndef HSEFUNCTIONAL_H
+#define HSEFUNCTIONAL_H
+
+#include "XCFunctional.h"
+#include "InteractionPotential.h"
+#include <vector>
+
+class HSEFunctional : public XCFunctional
+{
+  const double x_coeff_, c_coeff_;
+
+  // screening parameter of the HSE functional
+  static const double omega = 0.11;
+
+  // vectors common to all GGA exchange functionals
+  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;
+  std::vector<double> _grad_rho[3], _grad_rho_up[3], _grad_rho_dn[3];
+
+  public:
+
+  // constructor
+  HSEFunctional(const std::vector<std::vector<double> > &rhoe);
+
+  // HSE's local part is derived from PBE
+  bool isGGA() const
+  {
+    return true;
+  }
+
+  // return the name of the functional
+  std::string name() const
+  {
+    return "HSE";
+  }
+
+  void setxc(void);
+
+  // evaluate fourier transform of nonlocal potential for given G vector
+  // input g2 = G^2
+  static double interaction_potential(const double& g2);
+  // derivative of interaction potential
+  static double derivative_interaction_potential(const double& g2);
+  // scaling of the divergence correction relative to the Coulomb potential
+  static double divergence_scaling(const double& rcut);
+
+  // construct interaction potential class
+  static const InteractionPotential make_interaction_potential()
+  {
+    return InteractionPotential(&interaction_potential,
+      &derivative_interaction_potential,&divergence_scaling);
+  }
+
+};
+
+#endif

src/InteractionPotential.h

+////////////////////////////////////////////////////////////////////////////////
+//
+// 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/>.
+//
+////////////////////////////////////////////////////////////////////////////////
+//
+// InteractionPotential.h
+//
+////////////////////////////////////////////////////////////////////////////////
+//
+// Implements InteractionPotential class that evaluates the potential for given
+// norm of G vectors and the derivative w.r.t. this argument
+//
+////////////////////////////////////////////////////////////////////////////////
+
+#ifndef INTERACTIONPOTENTIAL_H
+#define INTERACTIONPOTENTIAL_H
+
+#include <cassert>
+
+class InteractionPotential
+{
+  public:
+
+  // default constructor = Coulomb potential
+  InteractionPotential() : coulomb_(true) {}
+
+  // constructor - define function and derivative
+  InteractionPotential(double(*V)(const double&), double(*dV)(const double&),
+    double(*div_scale)(const double&)) :
+    V_(V), dV_(dV), div_scale_(div_scale), coulomb_(false) {}
+
+  // is the interaction potential a coulomb potential?
+  inline bool coulomb() const
+  {
+    return coulomb_;
+  }
+
+  // evaluate the interaction potential for given norm of G vector
+  inline double operator()(const double G2) const
+  {
+    // the current implementation expects that the coulomb potential
+    // is treated externaly
+    assert(not coulomb_);
+    return V_(G2);
+  }
+
+  // evaluate the derivative of the interaction potential w.r.t. G^2
+  inline double derivative(const double G2) const
+  {
+    // the current implementation expects that the coulomb potential
+    // is treated externaly
+    assert(not coulomb_);
+    return dV_(G2);
+  }
+
+  inline double divergence_scaling(const double rcut) const
+  {
+    // the current implementation expects that the coulomb potential
+    // is treated externaly
+    assert(not coulomb_);
+    return div_scale_(rcut);
+  }
+
+  private:
+
+  const bool coulomb_;
+  double (*V_)(const double&);
+  double (*dV_)(const double&);
+  double (*div_scale_)(const double&);
+
+};
+
+#endif

src/Makefile

@@ -30,6 +30,7 @@ OBJECTS=qb.o AtomSet.o Atom.o Species.o \
         XCOperator.o ExchangeOperator.o Bisection.o \
 	XCPotential.o LDAFunctional.o VWNFunctional.o \
         PBEFunctional.o BLYPFunctional.o B3LYPFunctional.o  \
+	ExponentialIntegral.o HSEFunctional.o \
 	NonLocalPotential.o SampleReader.o StructuredDocumentHandler.o \
 	SampleHandler.o AtomSetHandler.o WavefunctionHandler.o \
         XMLGFPreprocessor.o Base64Transcoder.o \

src/XCOperator.C

@@ -19,6 +19,7 @@
 #include "ChargeDensity.h"
 #include "XCPotential.h"
 #include "ExchangeOperator.h"
+#include "HSEFunctional.h"
 using namespace std;
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -68,6 +69,19 @@ XCOperator::XCOperator(Sample& s, const ChargeDensity& cd) :cd_(cd)
     hasHF_ = true;
     HFmixCoeff_ = s.ctrl.alpha_PBE0;;
   }
+  else if ( functional_name == "HSE" )
+  {
+    // create an exchange potential
+    xcp_ = new XCPotential(cd, functional_name, s.ctrl);
+
+    // create the exchange operator with mixing coeff=0.25
+    xop_ = new ExchangeOperator(s, 0.25,
+      HSEFunctional::make_interaction_potential() );
+    hasPotential_ = true;
+    hasGGA_ = xcp_->isGGA();
+    hasHF_ = true;
+    HFmixCoeff_ = 0.25;
+  }
   else if ( functional_name == "B3LYP" )
   {
     // create an exchange potential

src/XCPotential.C

@@ -21,6 +21,7 @@
 #include "VWNFunctional.h"
 #include "PBEFunctional.h"
 #include "BLYPFunctional.h"
+#include "HSEFunctional.h"
 #include "B3LYPFunctional.h"
 #include "Basis.h"
 #include "FourierTransform.h"
@@ -54,6 +55,10 @@ XCPotential::XCPotential(const ChargeDensity& cd, const string functional_name,
     const double c_coeff = 1.0;
     xcf_ = new PBEFunctional(cd_.rhor,x_coeff,c_coeff);
   }
+  else if ( functional_name == "HSE" )
+  {
+    xcf_ = new HSEFunctional(cd_.rhor);
+  }
   else if ( functional_name == "B3LYP" )
   {
     xcf_ = new B3LYPFunctional(cd_.rhor);

src/XCPotential.h

 ////////////////////////////////////////////////////////////////////////////////
+//
 // Copyright (c) 2008 The Regents of the University of California
 //
 // This file is part of Qbox

src/Xc.h

@@ -50,10 +50,12 @@ class Xc : public Var
             v == "BLYP" ||
             v == "HF" ||
             v == "PBE0" ||
+            v == "HSE" ||
             v == "B3LYP" ) )
     {
       if ( ui->onpe0() )
-        cout << " xc must be LDA, VWN, PBE, BLYP, HF, PBE0 or B3LYP" << endl;
+        cout << " xc must be LDA, VWN, PBE, BLYP, "
+             << "HF, PBE0, HSE or B3LYP" << endl;
       return 1;
     }
 

util/qbox_dos/qbox_dos.py

+#!/usr/bin/python
+# qbox_dos.py: extract dos from Qbox output
+# generate dos plot in gnuplot format
+# use: qbox_dos.py emin emax width file.r
+# the DOS is accumulated separately for each spin
+
+import xml.sax
+import sys
+import math
+
+if len(sys.argv) != 5:
+  print "use: ",sys.argv[0]," emin emax width file.r"
+  sys.exit()
+
+emin = float(sys.argv[1])
+emax = float(sys.argv[2])
+width = float(sys.argv[3])
+infile = sys.argv[4]
+
+ndos = 501
+de = (emax - emin)/(ndos-1)
+
+# normalized gaussian distribution in one dimension
+# f(x) = 1/(sqrt(pi)*width) * exp(-(x/width)^2 )
+def gauss(x, width):
+  return (1.0/(math.sqrt(math.pi)*width)) * math.exp(-(x/width)**2)
+  
+# Qbox output handler to extract and process data
+class QboxOutputHandler(xml.sax.handler.ContentHandler):
+  def __init__(self):
+    self.nspin = 1
+    self.readData = 0
+    self.dos_up = [0] * ndos
+    self.dos_dn = [0] * ndos
+
+  def startElement(self, name, attributes):
+    if name == "eigenvalues":
+      self.n = attributes["n"]
+      self.spin = int(attributes["spin"])
+      self.kpoint = attributes["kpoint"]
+      self.weight = float(attributes["weight"])
+      self.readData = 1
+      self.buffer = ""
+      if self.spin == 1:
+        self.nspin = 2
+
+  def characters(self, data):
+    if self.readData:
+      self.buffer += data
+
+  def endElement(self, name):
+    if name == "eigenvalues":
+      self.readData = 0
+      self.accumulate_dos()
+
+  def accumulate_dos(self):
+    self.e = self.buffer.split()
+    if self.spin == 0:
+      for i in range(len(self.e)):
+        for j in range(ndos):
+          ej = emin + j * de
+          self.dos_up[j] += gauss(float(self.e[i])-ej, width ) * self.weight
+    if self.spin == 1:
+      for i in range(len(self.e)):
+        for j in range(ndos):
+          ej = emin + j * de
+          self.dos_dn[j] += gauss(float(self.e[i])-ej, width ) * self.weight
+
+  def print_dos(self):
+    print "# ",infile," spin=0 width=",width
+    for j in range(ndos):
+      ej = emin + j * de
+      print ej, self.dos_up[j]
+    if self.nspin == 2:
+      print
+      print
+      print "# ",infile," spin=1 width=",width
+      for j in range(ndos):
+        ej = emin + j * de
+        print ej, self.dos_dn[j]
+
+parser = xml.sax.make_parser()
+handler = QboxOutputHandler()
+parser.setContentHandler(handler)
+parser.parse(infile)
+handler.print_dos()

util/qbox_maxforce.sh

+#!/bin/bash
+# get the largest force component in a file of "<force> fx fy fz </force>"
+# use: qbox_maxforce nat file.r
+grep '<force>' $2 | tail -$1 | \
+awk '{
+  if ( mx*mx < $2*$2 ) mx = $2;
+  if ( my*my < $3*$3 ) my = $3;
+  if ( mz*mz < $4*$4 ) mz = $4;
+} END {printf("%9.2e %9.2e %9.2e\n", mx, my, mz)}' -