Commit e4447360 by Martin Schlipf

implement InteractionPotential class

for stress calculations not only the potential but the derivative is needed as well
create class that encapsulates potential and derivative
generator function in HSEFunctional to create the new type

git-svn-id: http://qboxcode.org/svn/qb/branches/hse-dev@1400 cba15fb0-1239-40c8-b417-11db7ca47a34
parent 1d033ec8
......@@ -964,3 +964,28 @@ double HSEFunctional::interaction_potential(const double& g2) {
}
}
// derivative of interaction potential
// input g2 = G^2
// exp( -g2 / 4 w^2 ) V(g2)
// ------------------ - -----
// 4 g2 w^2 g2
double HSEFunctional::derivative_interaction_potential(const double& g2) {
// helper variable
const double r1_4w2 = 0.25 / (omega * omega);
const double x = g2 * r1_4w2;
const double third = 1.0 / 3.0;
if (g2 == 0) {
// trivial limit for g2 = 0
return -0.5 * r1_4w2 * r1_4w2;
} else if (g2 < 1e-6) {
// taylor expansion near origin
return (-0.5 + x * (third - 0.125 * x)) * r1_4w2 * r1_4w2;
} else {
// exact derivative
return (exp(-x) * r1_4w2 - interaction_potential(g2)) / g2;
}
}
......@@ -30,6 +30,7 @@
#define HSEFUNCTIONAL_H
#include "XCFunctional.h"
#include "InteractionPotential.h"
#include <vector>
class HSEFunctional : public XCFunctional
......@@ -66,7 +67,16 @@ class HSEFunctional : public XCFunctional
// evaluate fourier transform of nonlocal potential for given G vector
// input g2 = G^2
static double interaction_potential( const double& g2 );
static double interaction_potential(const double& g2);
// derivative of interaction potential
static double derivative_interaction_potential(const double& g2);
// construct interaction potential class
static const InteractionPotential make_interaction_potential()
{
return InteractionPotential(&interaction_potential,
&derivative_interaction_potential);
}
};
......
////////////////////////////////////////////////////////////////////////////////
//
// 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
//
// Author: Martin Schlipf (2013)
// Contact: martin.schlipf@gmail.com
//
////////////////////////////////////////////////////////////////////////////////
#ifndef INTERACTIONPOTENTIAL_H
#define INTERACTIONPOTENTIAL_H
class InteractionPotential
{
public:
// constructor - define function and derivative
InteractionPotential(double(*V)(const double&), double(*dV)(const double&)) :
V_(V), dV_(dV)
{
}
// evaluate the interaction potential for given norm of G vector
inline double operator()(const double G2) const
{
return V_(G2);
}
// evaluate the derivative of the interaction potential w.r.t. G^2
inline double derivative(const double G2) const
{
return dV_(G2);
}
private:
double (*V_)(const double&);
double (*dV_)(const double&);
};
#endif
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