Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F81604385
KernelFunction.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sat, Sep 7, 00:05
Size
23 KB
Mime Type
text/x-c
Expires
Mon, Sep 9, 00:05 (2 d)
Engine
blob
Format
Raw Data
Handle
20545155
Attached To
rLAMMPS lammps
KernelFunction.cpp
View Options
#include "KernelFunction.h"
#include "math.h"
#include <vector>
#include "ATC_Error.h"
#include "Quadrature.h"
#include "Utility.h"
using namespace std;
using namespace ATC_Utility;
static const double tol = 1.0e-8;
static const int line_ngauss = 10;
static double line_xg[line_ngauss], line_wg[line_ngauss];
namespace ATC {
//========================================================================
// KernelFunctionMgr
//========================================================================
KernelFunctionMgr * KernelFunctionMgr::myInstance_ = NULL;
//------------------------------------------------------------------------
// instance
//------------------------------------------------------------------------
KernelFunctionMgr * KernelFunctionMgr::instance()
{
if (myInstance_ == NULL) {
myInstance_ = new KernelFunctionMgr();
}
return myInstance_;
}
//------------------------------------------------------------------------
// get function from args
//------------------------------------------------------------------------
KernelFunction* KernelFunctionMgr::function(char ** arg, int narg)
{
/*! \page man_kernel_function fix_modify AtC kernel
\section syntax
fix_modify AtC kernel <type> <parameters>
- type (keyword) = step, cell, cubic_bar, cubic_cylinder, cubic_sphere,
quartic_bar, quartic_cylinder, quartic_sphere \n
- parameters :\n
step = radius (double) \n
cell = hx, hy, hz (double) or h (double) \n
cubic_bar = half-width (double) \n
cubic_cylinder = radius (double) \n
cubic_sphere = radius (double) \n
quartic_bar = half-width (double) \n
quartic_cylinder = radius (double) \n
quartic_sphere = radius (double) \n
\section examples
fix_modify AtC kernel cell 1.0 1.0 1.0
fix_modify AtC kernel quartic_sphere 10.0
\section description
\section restrictions
Must be used with the hardy AtC fix \n
For bar kernel types, half-width oriented along x-direction \n
For cylinder kernel types, cylindrical axis is assumed to be in z-direction \n
( see \ref man_fix_atc )
\section related
\section default
No default
*/
int argIdx = 0;
KernelFunction * ptr = NULL;
char* type = arg[argIdx++];
if (strcmp(type,"step")==0) {
double parameters[1] = {atof(arg[argIdx])}; // cutoff radius
ptr = new KernelFunctionStep(1,parameters);
}
else if (strcmp(type,"cell")==0) {
double parameters[3];
parameters[0] = parameters[1] = parameters[2] = atof(arg[argIdx++]);
if (narg > argIdx) { // L_x, L_y, L_z
for (int i = 1; i < 3; i++) { parameters[i] = atof(arg[argIdx++]); }
}
ptr = new KernelFunctionCell(2,parameters);
}
else if (strcmp(type,"cubic_bar")==0) {
double parameters[1] = {atof(arg[argIdx])}; // cutoff half-length
ptr = new KernelFunctionCubicBar(1,parameters);
}
else if (strcmp(type,"linear_bar")==0) {
double parameters[1] = {atof(arg[argIdx])}; // cutoff half-length
ptr = new KernelFunctionLinearBar(1,parameters);
}
else if (strcmp(type,"cubic_cylinder")==0) {
double parameters[1] = {atof(arg[argIdx])}; // cutoff radius
ptr = new KernelFunctionCubicCyl(1,parameters);
}
else if (strcmp(type,"cubic_sphere")==0) {
double parameters[1] = {atof(arg[argIdx])}; // cutoff radius
ptr = new KernelFunctionCubicSphere(1,parameters);
}
else if (strcmp(type,"quartic_bar")==0) {
double parameters[1] = {atof(arg[argIdx])}; // cutoff half-length
ptr = new KernelFunctionQuarticBar(1,parameters);
}
else if (strcmp(type,"quartic_cylinder")==0) {
double parameters[1] = {atof(arg[argIdx])}; // cutoff radius
ptr = new KernelFunctionQuarticCyl(1,parameters);
}
else if (strcmp(type,"quartic_sphere")==0) {
double parameters[1] = {atof(arg[argIdx])}; // cutoff radius
ptr = new KernelFunctionQuarticSphere(1,parameters);
}
pointerSet_.insert(ptr);
return ptr;
}
// Destructor
KernelFunctionMgr::~KernelFunctionMgr()
{
set<KernelFunction * >::iterator it;
for (it = pointerSet_.begin(); it != pointerSet_.end(); it++)
if (*it) delete *it;
}
//------------------------------------------------------------------------
// KernelFunction
//------------------------------------------------------------------------
// constructor
KernelFunction::KernelFunction(int nparameters, double* parameters):
Rc_(0),invRc_(0),nsd_(3),
lammpsInterface_(LammpsInterface::instance())
{
Rc_ = parameters[0];
invRc_ = 1.0/Rc_;
Rc_ = parameters[0];
invRc_ = 1.0/Rc_;
invVol_ = 1.0/(4.0/3.0*Pi_*pow(Rc_,3));
ATC::Quadrature::instance()->set_line_quadrature(line_ngauss,line_xg,line_wg);
// get periodicity and box bounds/lengths
lammpsInterface_->box_periodicity(periodicity[0],
periodicity[1],periodicity[2]);
lammpsInterface_->box_bounds(box_bounds[0][0],box_bounds[1][0],
box_bounds[0][1],box_bounds[1][1],
box_bounds[0][2],box_bounds[1][2]);
for (int k = 0; k < 3; k++) {
box_length[k] = box_bounds[1][k] - box_bounds[0][k];
}
}
// does an input node's kernel intersect bonds on this processor
bool KernelFunction::node_contributes(DENS_VEC node) const
{
DENS_VEC ghostNode = node;
bool contributes = true;
bool ghostContributes = lammpsInterface_->nperiodic();
double kernel_bounds[2][3];
lammpsInterface_->sub_bounds(kernel_bounds[0][0],kernel_bounds[1][0],
kernel_bounds[0][1],kernel_bounds[1][1],
kernel_bounds[0][2],kernel_bounds[1][2]);
for (int i=0; i<3; ++i) {
if (i < nsd_) {
kernel_bounds[0][i] -= (Rc_+lammpsInterface_->pair_cutoff());
kernel_bounds[1][i] += (Rc_+lammpsInterface_->pair_cutoff());
contributes = contributes && (node(i) >= kernel_bounds[0][i])
&& (node(i) < kernel_bounds[1][i]);
if (periodicity[i]) {
if (node[i] <= box_bounds[0][i] + box_length[i]/2) {
ghostNode[i] += box_length[i];
} else {
ghostNode[i] -= box_length[i];
}
ghostContributes = ghostContributes
&& ((ghostNode(i) >= kernel_bounds[0][i]) ||
(node(i) >= kernel_bounds[0][i]))
&& ((ghostNode(i) < kernel_bounds[1][i]) ||
(node(i) < kernel_bounds[1][i]));
}
}
if (!(contributes || ghostContributes)) break;
}
return true;
}
bool KernelFunction::in_support(DENS_VEC dx) const
{
if (dx.norm() > Rc_) return false;
return true;
}
// bond function value via quadrature
double KernelFunction::bond(DENS_VEC& xa, DENS_VEC&xb, double lam1, double lam2) const
{
DENS_VEC xab(nsd_), q(nsd_);
double lamg;
double bhsum=0.0;
xab = xa - xb;
for (int i = 0; i < line_ngauss; i++) {
lamg=0.5*((lam2-lam1)*line_xg[i]+(lam2+lam1));
q = lamg*xab + xb;
double locg_value=this->value(q);
bhsum+=locg_value*line_wg[i];
}
return 0.5*(lam2-lam1)*bhsum;
}
// localization-volume intercepts for bond calculation
// bond intercept values assuming spherical support
void KernelFunction::bond_intercepts(DENS_VEC& xa,
DENS_VEC& xb, double &lam1, double &lam2) const
{
if (nsd_ == 2) {// for cylinders, axis is always z!
const int iaxis = 2;
xa[iaxis] = 0.0;
xb[iaxis] = 0.0;
} else if (nsd_ == 1) {// for bars, analysis is 1D in x
xa[1] = 0.0;
xa[2] = 0.0;
xb[1] = 0.0;
xb[2] = 0.0;
}
lam1 = lam2 = -1;
double ra_n = invRc_*xa.norm(); // lambda = 1
double rb_n = invRc_*xb.norm(); // lambda = 0
bool a_in = (ra_n <= 1.0);
bool b_in = (rb_n <= 1.0);
if (a_in && b_in) {
lam1 = 0.0;
lam2 = 1.0;
return;
}
DENS_VEC xab = xa - xb;
double rab_n = invRc_*xab.norm();
double a = rab_n*rab_n; // always at least an interatomic distance
double b = 2.0*invRc_*invRc_*xab.dot(xb);
double c = rb_n*rb_n - 1.0;
double discrim = b*b - 4.0*a*c; // discriminant
if (discrim < 0) return; // no intersection
// num recipes:
double s1, s2;
if (b < 0.0) {
double aux = -0.5*(b-sqrt(discrim));
s1 = c/aux;
s2 = aux/a;
}
else {
double aux = -0.5*(b+sqrt(discrim));
s1 = aux/a;
s2 = c/aux;
}
if (a_in && !b_in) {
lam1 = s1;
lam2 = 1.0;
}
else if (!a_in && b_in) {
lam1 = 0.0;
lam2 = s2;
}
else {
if (s1 >= 0.0 && s2 <= 1.0) {
lam1 = s1;
lam2 = s2;
}
}
}
//------------------------------------------------------------------------
// constructor
KernelFunctionStep::KernelFunctionStep
(int nparameters, double* parameters):
KernelFunction(nparameters, parameters)
{
for (int k = 0; k < nsd_; k++ ) {
if ((bool) periodicity[k]) {
if (Rc_ > 0.5*box_length[k]) {
throw ATC_Error("Size of localization volume is too large for periodic boundary condition");
}
}
}
}
// function value
double KernelFunctionStep::value(DENS_VEC& x_atom) const
{
double rn=invRc_*x_atom.norm();
if (rn <= 1.0) { return 1.0; }
else { return 0.0; }
}
// function derivative value
void KernelFunctionStep::derivative(const DENS_VEC& x_atom, DENS_VEC& deriv) const
{
deriv.reset(nsd_);
}
//------------------------------------------------------------------------
/** a step with rectangular support suitable for a rectangular grid */
// constructor
KernelFunctionCell::KernelFunctionCell
(int nparameters, double* parameters):
KernelFunction(nparameters, parameters)
{
hx = parameters[0];
hy = parameters[1];
hz = parameters[2];
invVol_ = 1.0/8.0/(hx*hy*hz);
cellBounds_.reset(6);
cellBounds_(0) = -hx;
cellBounds_(1) = hx;
cellBounds_(2) = -hy;
cellBounds_(3) = hy;
cellBounds_(4) = -hz;
cellBounds_(5) = hz;
for (int k = 0; k < nsd_; k++ ) {
if ((bool) periodicity[k]) {
if (parameters[k] > 0.5*box_length[k]) {
throw ATC_Error("Size of localization volume is too large for periodic boundary condition");
}
}
}
}
// does an input node's kernel intersect bonds on this processor
bool KernelFunctionCell::node_contributes(DENS_VEC node) const
{
DENS_VEC ghostNode = node;
bool contributes = true;
bool ghostContributes = lammpsInterface_->nperiodic();
double kernel_bounds[2][3];
lammpsInterface_->sub_bounds(kernel_bounds[0][0],kernel_bounds[1][0],
kernel_bounds[0][1],kernel_bounds[1][1],
kernel_bounds[0][2],kernel_bounds[1][2]);
for (int i=0; i<3; ++i) {
kernel_bounds[0][i] -= (cellBounds_(i*2+1) +
lammpsInterface_->pair_cutoff());
kernel_bounds[1][i] += (cellBounds_(i*2+1) +
lammpsInterface_->pair_cutoff());
contributes = contributes && (node(i) >= kernel_bounds[0][i])
&& (node(i) < kernel_bounds[1][i]);
if (periodicity[i]) {
if (node[i] <= box_bounds[0][i] + box_length[i]/2) {
ghostNode[i] += box_length[i];
} else {
ghostNode[i] -= box_length[i];
}
ghostContributes = ghostContributes
&& ((ghostNode(i) >= kernel_bounds[0][i]) ||
(node(i) >= kernel_bounds[0][i]))
&& ((ghostNode(i) < kernel_bounds[1][i]) ||
(node(i) < kernel_bounds[1][i]));
}
if (!(contributes || ghostContributes)) break;
}
return true;
}
bool KernelFunctionCell::in_support(DENS_VEC dx) const
{
if (dx(0) < cellBounds_(0)
|| dx(0) > cellBounds_(1)
|| dx(1) < cellBounds_(2)
|| dx(1) > cellBounds_(3)
|| dx(2) < cellBounds_(4)
|| dx(2) > cellBounds_(5) ) return false;
return true;
}
// function value
double KernelFunctionCell::value(DENS_VEC& x_atom) const
{
if ((cellBounds_(0) <= x_atom(0)) && (x_atom(0) < cellBounds_(1))
&& (cellBounds_(2) <= x_atom(1)) && (x_atom(1) < cellBounds_(3))
&& (cellBounds_(4) <= x_atom(2)) && (x_atom(2) < cellBounds_(5))) {
return 1.0;
}
else {
return 0.0;
}
}
// function derivative value
void KernelFunctionCell::derivative(const DENS_VEC& x_atom, DENS_VEC& deriv) const
{
deriv.reset(nsd_);
}
// bond intercept values for rectangular region : origin is the node position
void KernelFunctionCell::bond_intercepts(DENS_VEC& xa,
DENS_VEC& xb, double &lam1, double &lam2) const
{
lam1 = 0.0; // start
lam2 = 1.0; // end
bool a_in = (value(xa) > 0.0);
bool b_in = (value(xb) > 0.0);
// (1) both in, no intersection needed
if (a_in && b_in) {
return;
}
// (2) for one in & one out -> one plane intersection
// determine the points of intersection between the line joining
// atoms a and b and the bounding planes of the localization volume
else if (a_in || b_in) {
DENS_VEC xab = xa - xb;
for (int i = 0; i < nsd_; i++) {
// check if segment is parallel to face
if (fabs(xab(i)) > tol) {
for (int j = 0; j < 2; j++) {
double s = (cellBounds_(2*i+j) - xb(i))/xab(i);
// check if between a & b
if (s >= 0 && s <= 1) {
bool in_bounds = false;
DENS_VEC x = xb + s*xab;
if (i == 0) {
if ((cellBounds_(2) <= x(1)) && (x(1) <= cellBounds_(3))
&& (cellBounds_(4) <= x(2)) && (x(2) <= cellBounds_(5))) {
in_bounds = true;
}
}
else if (i == 1) {
if ((cellBounds_(0) <= x(0)) && (x(0) <= cellBounds_(1))
&& (cellBounds_(4) <= x(2)) && (x(2) <= cellBounds_(5))) {
in_bounds = true;
}
}
else if (i == 2) {
if ((cellBounds_(0) <= x(0)) && (x(0) <= cellBounds_(1))
&& (cellBounds_(2) <= x(1)) && (x(1) <= cellBounds_(3))) {
in_bounds = true;
}
}
if (in_bounds) {
if (a_in) { lam1 = s;}
else { lam2 = s;}
return;
}
}
}
}
}
throw ATC_Error("logic failure in HardyKernel Cell for single intersection\n");
}
// (3) both out -> corner intersection
else {
lam2 = lam1; // default to no intersection
DENS_VEC xab = xa - xb;
double ss[6] = {-1,-1,-1,-1,-1,-1};
int is = 0;
for (int i = 0; i < nsd_; i++) {
// check if segment is parallel to face
if (fabs(xab(i)) > tol) {
for (int j = 0; j < 2; j++) {
double s = (cellBounds_(2*i+j) - xb(i))/xab(i);
// check if between a & b
if (s >= 0 && s <= 1) {
// check if in face
DENS_VEC x = xb + s*xab;
if (i == 0) {
if ((cellBounds_(2) <= x(1)) && (x(1) <= cellBounds_(3))
&& (cellBounds_(4) <= x(2)) && (x(2) <= cellBounds_(5))) {
ss[is++] = s;
}
}
else if (i == 1) {
if ((cellBounds_(0) <= x(0)) && (x(0) <= cellBounds_(1))
&& (cellBounds_(4) <= x(2)) && (x(2) <= cellBounds_(5))) {
ss[is++] = s;
}
}
else if (i == 2) {
if ((cellBounds_(0) <= x(0)) && (x(0) <= cellBounds_(1))
&& (cellBounds_(2) <= x(1)) && (x(1) <= cellBounds_(3))) {
ss[is++] = s;
}
}
}
}
}
}
if (is == 1) {
// intersection occurs at a box edge - leave lam1 = lam2
}
else if (is == 2) {
lam1 = min(ss[0],ss[1]);
lam2 = max(ss[0],ss[1]);
}
else if (is == 3) {
// intersection occurs at a box vertex - leave lam1 = lam2
}
else {
if (is != 0) throw ATC_Error("logic failure in HardyKernel Cell for corner intersection\n");
}
}
}
//------------------------------------------------------------------------
// constructor
KernelFunctionCubicSphere::KernelFunctionCubicSphere
(int nparameters, double* parameters):
KernelFunction(nparameters, parameters)
{
for (int k = 0; k < nsd_; k++ ) {
if ((bool) periodicity[k]) {
if (Rc_ > 0.5*box_length[k]) {
throw ATC_Error("Size of localization volume is too large for periodic boundary condition");
};
};
};
}
// function value
double KernelFunctionCubicSphere::value(DENS_VEC& x_atom) const
{
double r=x_atom.norm();
double rn=r/Rc_;
if (rn < 1.0) { return 5.0*(1.0-3.0*rn*rn+2.0*rn*rn*rn); }
else { return 0.0; }
}
// function derivative value
void KernelFunctionCubicSphere::derivative(const DENS_VEC& x_atom, DENS_VEC& deriv) const
{
deriv.reset(nsd_);
}
//------------------------------------------------------------------------
// constructor
KernelFunctionQuarticSphere::KernelFunctionQuarticSphere
(int nparameters, double* parameters):
KernelFunction(nparameters, parameters)
{
for (int k = 0; k < nsd_; k++ ) {
if ((bool) periodicity[k]) {
if (Rc_ > 0.5*box_length[k]) {
throw ATC_Error("Size of localization volume is too large for periodic boundary condition");
};
};
};
}
// function value
double KernelFunctionQuarticSphere::value(DENS_VEC& x_atom) const
{
double r=x_atom.norm();
double rn=r/Rc_;
if (rn < 1.0) { return 35.0/8.0*pow((1.0-rn*rn),2); }
else { return 0.0; }
}
// function derivative value
void KernelFunctionQuarticSphere::derivative(const DENS_VEC& x_atom, DENS_VEC& deriv) const
{
deriv.reset(nsd_);
}
//------------------------------------------------------------------------
// constructor
KernelFunctionCubicCyl::KernelFunctionCubicCyl
(int nparameters, double* parameters):
KernelFunction(nparameters, parameters)
{
nsd_ = 2;
double Lz = box_length[2];
invVol_ = 1.0/(Pi_*pow(Rc_,2)*Lz);
for (int k = 0; k < nsd_; k++ ) {
if ((bool) periodicity[k]) {
if (Rc_ > 0.5*box_length[k]) {
throw ATC_Error("Size of localization volume is too large for periodic boundary condition");
};
};
};
}
// function value
double KernelFunctionCubicCyl::value(DENS_VEC& x_atom) const
{
double r=sqrt(pow(x_atom(0),2)+pow(x_atom(1),2));
double rn=r/Rc_;
if (rn < 1.0) { return 10.0/3.0*(1.0-3.0*rn*rn+2.0*rn*rn*rn); }
else { return 0.0; }
}
// function derivative value
void KernelFunctionCubicCyl::derivative(const DENS_VEC& x_atom, DENS_VEC& deriv) const
{
deriv.reset(nsd_);
}
//------------------------------------------------------------------------
// constructor
KernelFunctionQuarticCyl::KernelFunctionQuarticCyl
(int nparameters, double* parameters):
KernelFunction(nparameters, parameters)
{
nsd_ = 2;
double Lz = box_length[2];
invVol_ = 1.0/(Pi_*pow(Rc_,2)*Lz);
for (int k = 0; k < nsd_; k++ ) {
if ((bool) periodicity[k]) {
if (Rc_ > 0.5*box_length[k]) {
throw ATC_Error("Size of localization volume is too large for periodic boundary condition");
};
};
};
}
// function value
double KernelFunctionQuarticCyl::value(DENS_VEC& x_atom) const
{
double r=sqrt(pow(x_atom(0),2)+pow(x_atom(1),2));
double rn=r/Rc_;
if (rn < 1.0) { return 3.0*pow((1.0-rn*rn),2); }
else { return 0.0; }
}
// function derivative value
void KernelFunctionQuarticCyl::derivative(const DENS_VEC& x_atom, DENS_VEC& deriv) const
{
deriv.reset(nsd_);
}
//------------------------------------------------------------------------
// constructor
KernelFunctionCubicBar::KernelFunctionCubicBar
(int nparameters, double* parameters):
KernelFunction(nparameters, parameters)
{
// Note: Bar is assumed to be oriented in the x(0) direction
nsd_ = 1;
double Ly = box_length[1];
double Lz = box_length[2];
invVol_ = 1.0/(2*Rc_*Ly*Lz);
if ((bool) periodicity[0]) {
if (Rc_ > 0.5*box_length[0]) {
throw ATC_Error("Size of localization volume is too large for periodic boundary condition");
};
};
}
// function value
double KernelFunctionCubicBar::value(DENS_VEC& x_atom) const
{
double r=abs(x_atom(0));
double rn=r/Rc_;
if (rn < 1.0) { return 2.0*(1.0-3.0*rn*rn+2.0*rn*rn*rn); }
else { return 0.0; }
}
// function derivative value
void KernelFunctionCubicBar::derivative(const DENS_VEC& x_atom, DENS_VEC& deriv) const
{
deriv.reset(nsd_);
}
//------------------------------------------------------------------------
// constructor
KernelFunctionLinearBar::KernelFunctionLinearBar
(int nparameters, double* parameters):
KernelFunction(nparameters, parameters)
{
// Note: Bar is assumed to be oriented in the z(0) direction
double Lx = box_length[0];
double Ly = box_length[1];
invVol_ = 1.0/(Lx*Ly*Rc_);
if ((bool) periodicity[2]) {
if (Rc_ > 0.5*box_length[2]) {
throw ATC_Error("Size of localization volume is too large for periodic boundary condition");
};
};
}
// function value
double KernelFunctionLinearBar::value(DENS_VEC& x_atom) const
{
double r=abs(x_atom(2));
double rn=r/Rc_;
if (rn < 1.0) { return 1.0-rn; }
else { return 0.0; }
}
// function derivative value
void KernelFunctionLinearBar::derivative(const DENS_VEC& x_atom, DENS_VEC& deriv) const
{
deriv.reset(nsd_);
deriv(0) = 0.0;
deriv(1) = 0.0;
double r=abs(x_atom(2));
double rn=r/Rc_;
if (rn < 1.0 && x_atom(2) <= 0.0) { deriv(2) = 1.0/Rc_; }
else if (rn < 1.0 && x_atom(2) > 0.0) { deriv(2) = -1.0/Rc_; }
else { deriv(2) = 0.0; }
}
//------------------------------------------------------------------------
// constructor
KernelFunctionQuarticBar::KernelFunctionQuarticBar
(int nparameters, double* parameters):
KernelFunction(nparameters, parameters)
{
// Note: Bar is assumed to be oriented in the x(0) direction
nsd_ = 1;
double Ly = box_length[1];
double Lz = box_length[2];
invVol_ = 1.0/(2*Rc_*Ly*Lz);
if ((bool) periodicity[0]) {
if (Rc_ > 0.5*box_length[0]) {
throw ATC_Error("Size of localization volume is too large for periodic boundary condition");
};
};
}
// function value
double KernelFunctionQuarticBar::value(DENS_VEC& x_atom) const
{
double r=abs(x_atom(0));
double rn=r/Rc_;
// if (rn < 1.0) { return 5.0/2.0*(1.0-6*rn*rn+8*rn*rn*rn-3*rn*rn*rn*rn); } - alternative quartic
if (rn < 1.0) { return 15.0/8.0*pow((1.0-rn*rn),2); }
else { return 0.0; }
}
// function derivative value
void KernelFunctionQuarticBar::derivative(const DENS_VEC& x_atom, DENS_VEC& deriv) const
{
deriv.reset(nsd_);
}
};
Event Timeline
Log In to Comment