Skip to content

Commit

Permalink
Add TD4Model class (#17)
Browse files Browse the repository at this point in the history 
- add class for D4 model consistent with Fortran
- allows passing scaling parameters for zeta function
  • Loading branch information
@marvinfriede
marvinfriede authored Jan 18, 2023
1 parent cb34682 commit 963fdc5
Show file tree
Hide file tree
Showing 13 changed files with 284 additions and 126 deletions.
3 changes: 3 additions & 0 deletions include/dftd_dispersion.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,7 @@

#include "dftd_cutoff.h"
#include "dftd_geometry.h"
#include "dftd_model.h"

namespace dftd4 {

Expand All @@ -43,6 +44,7 @@ class dparam {
* @param mol Molecular geometry.
* @param charge Molecular charge.
* @param par DFT-D4 parameters.
* @param d4 Base D4 dispersion model.
* @param cutoff Real-space cutoffs for CN and dispersion.
* @param energy Dispersion energy.
* @param GRAD Dispersion gradient.
Expand All @@ -51,6 +53,7 @@ class dparam {
extern int get_dispersion(
const TMolecule &mol,
const int charge,
TD4Model &d4,
const dparam &par,
TCutoff cutoff,
double &energy,
Expand Down
68 changes: 45 additions & 23 deletions include/dftd_model.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,51 @@

namespace dftd4 {

// Default weighting factor for coordination number interpolation
static const double wf_default = 6.0;

// Default maximum charge scaling height for partial charge extrapolation
static const double ga_default = 3.0;

// Default charge scaling steepness for partial charge extrapolation
static const double gc_default = 2.0;

class TD4Model {
public:
double wf;
double ga;
double gc;

explicit TD4Model(
double wf_scale = wf_default,
double ga_scale = ga_default,
double gc_scale = gc_default
);

int weight_references(
const TMolecule &mol,
const TVector<double> &cn,
const TVector<double> &q,
TMatrix<double> &gwvec,
TMatrix<double> &dgwdcn,
TMatrix<double> &dgwdq,
bool lgrad = false
);

int get_atomic_c6(
const TMolecule &mol,
const TMatrix<double> &gwvec,
const TMatrix<double> &dgwdcn,
const TMatrix<double> &dgwdq,
TMatrix<double> &c6,
TMatrix<double> &dc6dcn,
TMatrix<double> &dc6dq,
bool lgrad = false
);

int set_refalpha_eeq(const TMolecule &mol, TMatrix<double> &alpha);
};

extern inline double trapzd(const double a[23], const double b[23]);

extern inline double
Expand All @@ -35,29 +80,6 @@ extern inline double

extern int get_max_ref(const TMolecule &mol, int &mref);

extern int weight_references(
const TMolecule &mol,
const TVector<double> &cn,
const TVector<double> &q,
TMatrix<double> &gwvec,
TMatrix<double> &dgwdcn,
TMatrix<double> &dgwdq,
bool lgrad = false
);

extern int get_atomic_c6(
const TMolecule &mol,
const TMatrix<double> &gwvec,
const TMatrix<double> &dgwdcn,
const TMatrix<double> &dgwdq,
TMatrix<double> &c6,
TMatrix<double> &dc6dcn,
TMatrix<double> &dc6dq,
bool lgrad = false
);

extern int set_refalpha_eeq(const TMolecule &mol, TMatrix<double> &alpha);

extern bool is_exceptional(double val);

} // namespace dftd4
9 changes: 5 additions & 4 deletions src/dftd_dispersion.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,7 @@ namespace dftd4 {
int get_dispersion(
const TMolecule &mol,
const int charge,
TD4Model &d4,
const dparam &par,
TCutoff cutoff,
double &energy,
Expand Down Expand Up @@ -87,7 +88,7 @@ int get_dispersion(
dgwdcn.NewMat(mref, mol.NAtoms);
dgwdq.NewMat(mref, mol.NAtoms);
}
info = weight_references(mol, cn, q, gwvec, dgwdcn, dgwdq, lgrad);
info = d4.weight_references(mol, cn, q, gwvec, dgwdcn, dgwdq, lgrad);
if (!info == EXIT_SUCCESS) return info;

TMatrix<double> c6;
Expand All @@ -99,7 +100,7 @@ int get_dispersion(
dc6dq.NewMat(mol.NAtoms, mol.NAtoms);
}

info = get_atomic_c6(mol, gwvec, dgwdcn, dgwdq, c6, dc6dcn, dc6dq, lgrad);
info = d4.get_atomic_c6(mol, gwvec, dgwdcn, dgwdq, c6, dc6dcn, dc6dq, lgrad);
if (!info == EXIT_SUCCESS) return info;

// --------------------------
Expand Down Expand Up @@ -154,7 +155,7 @@ int get_dispersion(
dgwdcn.NewMat(mref, mol.NAtoms);
dgwdq.NewMat(mref, mol.NAtoms);
}
info = weight_references(mol, cn, q, gwvec, dgwdcn, dgwdq, lgrad);
info = d4.weight_references(mol, cn, q, gwvec, dgwdcn, dgwdq, lgrad);
if (!info == EXIT_SUCCESS) return info;

cn.Delete();
Expand All @@ -166,7 +167,7 @@ int get_dispersion(
dc6dcn.NewMat(mol.NAtoms, mol.NAtoms);
dc6dq.NewMat(mol.NAtoms, mol.NAtoms);
}
info = get_atomic_c6(mol, gwvec, dgwdcn, dgwdq, c6, dc6dcn, dc6dq, lgrad);
info = d4.get_atomic_c6(mol, gwvec, dgwdcn, dgwdq, c6, dc6dcn, dc6dq, lgrad);
if (!info == EXIT_SUCCESS) return info;

gwvec.Delete();
Expand Down
119 changes: 62 additions & 57 deletions src/dftd_model.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,11 +26,6 @@

namespace dftd4 {

// Parameters
const double wf = 6.0;
const double g_a = 3.0;
const double g_c = 2.0;

static const double pi = 3.141592653589793238462643383279502884197;
static const double thopi = 3.0 / pi;

Expand Down Expand Up @@ -64,50 +59,17 @@ static const double weights[23]{
(freq[21] - freq[20]) + (freq[22] - freq[21]),
(freq[22] - freq[21])};

inline double trapzd(const double a[23], const double b[23]) {
double c6 = 0.0;
for (int w = 0; w != 23; w++)
c6 += weights[w] * a[w] * b[w];
return 0.5 * c6;
}

inline double weight_cn(const double wf, const double cn, const double cnref) {
double dcn = cn - cnref;
double arg = wf * dcn * dcn;
return exp(-arg);
}

inline double
zeta(const double a, const double c, const double qref, const double qmod) {
if (qmod <= 0.0) {
return exp(a);
} else {
return exp(a * (1.0 - exp(c * (1.0 - qref / qmod))));
}
}

inline double
dzeta(const double a, const double c, const double qref, const double qmod) {
if (qmod <= 0.0) {
return 0.0;
} else {
return -a * c * exp(c * (1.0 - qref / qmod)) * zeta(a, c, qref, qmod) *
qref / pow(qmod, 2);
}
}

int get_max_ref(const TMolecule &mol, int &mref) {
mref = refn[mol.at(0)];
for (int i = 1; i != mol.NAtoms; i++) {
int val = refn[mol.at(i)];
if (val > mref) mref = val;
}

if (mref == 0.0) return EXIT_FAILURE;
return EXIT_SUCCESS;
}
TD4Model::TD4Model(
double wf_scale /*= wf_default*/,
double ga_scale /*= ga_default*/,
double gc_scale /*= gc_default*/
) {
wf = wf_scale;
ga = ga_scale;
gc = gc_scale;
};

int weight_references(
int TD4Model::weight_references(
const TMolecule &mol,
const TVector<double> &cn,
const TVector<double> &q,
Expand All @@ -126,7 +88,7 @@ int weight_references(
for (int iat = 0; iat != mol.NAtoms; iat++) {
izp = mol.at(iat);
zi = zeff[izp];
gi = gam[izp] * g_c;
gi = gam[izp] * gc;

norm = 0.0;
dnorm = 0.0;
Expand Down Expand Up @@ -160,21 +122,21 @@ int weight_references(
}

gwvec(iref, iat) =
gwk * zeta(g_a, gi, refq[izp][iref] + zi, q(iat) + zi);
gwk * zeta(ga, gi, refq[izp][iref] + zi, q(iat) + zi);
dgwdq(iref, iat) =
gwk * dzeta(g_a, gi, refq[izp][iref] + zi, q(iat) + zi);
gwk * dzeta(ga, gi, refq[izp][iref] + zi, q(iat) + zi);

dgwk = norm * (dexpw - expw * dnorm * norm);
if (is_exceptional(dgwk)) { dgwk = 0.0; }
dgwdcn(iref, iat) =
dgwk * zeta(g_a, gi, refq[izp][iref] + zi, q(iat) + zi);
dgwk * zeta(ga, gi, refq[izp][iref] + zi, q(iat) + zi);
}
}
} else {
for (int iat = 0; iat != mol.NAtoms; iat++) {
izp = mol.at(iat);
zi = zeff[izp];
gi = gam[izp] * g_c;
gi = gam[izp] * gc;

norm = 0.0;
maxcn = 0.0;
Expand Down Expand Up @@ -202,15 +164,15 @@ int weight_references(
}

gwvec(iref, iat) =
gwk * zeta(g_a, gi, refq[izp][iref] + zi, q(iat) + zi);
gwk * zeta(ga, gi, refq[izp][iref] + zi, q(iat) + zi);
}
}
}

return EXIT_SUCCESS;
}

int get_atomic_c6(
int TD4Model::get_atomic_c6(
const TMolecule &mol,
const TMatrix<double> &gwvec,
const TMatrix<double> &dgwdcn,
Expand Down Expand Up @@ -294,7 +256,7 @@ int get_atomic_c6(
return EXIT_SUCCESS;
}

int set_refalpha_eeq(const TMolecule &mol, TMatrix<double> &alpha) {
int TD4Model::set_refalpha_eeq(const TMolecule &mol, TMatrix<double> &alpha) {
int iat{0}, is{0};
double iz{0.0}, aiw{0.0};

Expand All @@ -305,7 +267,7 @@ int set_refalpha_eeq(const TMolecule &mol, TMatrix<double> &alpha) {
iz = zeff[is];
for (int k = 0; k != 23; k++) {
aiw = secscale[is] * secalpha[is][k] *
zeta(g_a, gam[is] * g_c, iz, refsq[iat][ir] + iz);
zeta(ga, gam[is] * gc, iz, refsq[iat][ir] + iz);
alpha(i, 23 * ir + k) = std::max(
0.0,
refascale[iat][ir] *
Expand All @@ -318,6 +280,49 @@ int set_refalpha_eeq(const TMolecule &mol, TMatrix<double> &alpha) {
return EXIT_SUCCESS;
}

inline double trapzd(const double a[23], const double b[23]) {
double c6 = 0.0;
for (int w = 0; w != 23; w++)
c6 += weights[w] * a[w] * b[w];
return 0.5 * c6;
}

inline double weight_cn(const double wf, const double cn, const double cnref) {
double dcn = cn - cnref;
double arg = wf * dcn * dcn;
return exp(-arg);
}

inline double
zeta(const double a, const double c, const double qref, const double qmod) {
if (qmod <= 0.0) {
return exp(a);
} else {
return exp(a * (1.0 - exp(c * (1.0 - qref / qmod))));
}
}

inline double
dzeta(const double a, const double c, const double qref, const double qmod) {
if (qmod <= 0.0) {
return 0.0;
} else {
return -a * c * exp(c * (1.0 - qref / qmod)) * zeta(a, c, qref, qmod) *
qref / pow(qmod, 2);
}
}

int get_max_ref(const TMolecule &mol, int &mref) {
mref = refn[mol.at(0)];
for (int i = 1; i != mol.NAtoms; i++) {
int val = refn[mol.at(i)];
if (val > mref) mref = val;
}

if (mref == 0.0) return EXIT_FAILURE;
return EXIT_SUCCESS;
}

bool is_exceptional(double val) {
return std::isnan(val) || (fabs(val) > std::numeric_limits<double>::max());
}
Expand Down
6 changes: 4 additions & 2 deletions src/program_dftd.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@
#include "dftd_damping.h"
#include "dftd_dispersion.h"
#include "dftd_geometry.h"
#include "dftd_model.h"
#include "dftd_matrix.h"
#include "dftd_readxyz.h"

Expand Down Expand Up @@ -144,10 +145,11 @@ int main(int argc, char **argv) {
// readin the geometry file
read_xyzfile(fname, mol);

// initialize default cutoffs
// initialize default cutoffs and default D4 model
dftd4::TCutoff cutoff;
dftd4::TD4Model d4;

info = dftd4::get_dispersion(mol, charge, par, cutoff, energy, nullptr);
info = dftd4::get_dispersion(mol, charge, d4, par, cutoff, energy, nullptr);
if (info != 0) return EXIT_FAILURE;

std::cout << "Dispersion energy: " << energy << " Eh\n";
Expand Down
Loading

0 comments on commit 963fdc5

Please sign in to comment.