Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Minor fixes #32

Merged
merged 2 commits into from
Aug 1, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 4 additions & 1 deletion include/dftd_model.h
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,8 @@ class TD4Model {
double wf_scale = wf_default
);

virtual ~TD4Model() {};

int weight_references(
const TMolecule &mol,
const TVector<double> &cn,
Expand All @@ -68,7 +70,8 @@ class TD4Model {

virtual int set_refq_eeq(const TMolecule &mol, TMatrix<double> &refq) const;

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

extern inline double trapzd(const double a[23], const double b[23]);
Expand Down
3 changes: 2 additions & 1 deletion src/dftd_dispersion.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -174,7 +174,8 @@ int get_dispersion(
dc6dcn.NewMat(mol.NAtoms, mol.NAtoms);
dc6dq.NewMat(mol.NAtoms, mol.NAtoms);
}
info = d4.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
14 changes: 7 additions & 7 deletions src/dftd_model.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -57,7 +57,8 @@ static const double weights[23]{
(freq[19] - freq[18]) + (freq[20] - freq[19]),
(freq[20] - freq[19]) + (freq[21] - freq[20]),
(freq[21] - freq[20]) + (freq[22] - freq[21]),
(freq[22] - freq[21])};
(freq[22] - freq[21])
};

TD4Model::TD4Model(
double ga_scale /*= ga_default*/,
Expand Down Expand Up @@ -123,14 +124,14 @@ int TD4Model::weight_references(
}

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

dgwk = norm * (dexpw - expw * dnorm * norm);
if (is_exceptional(dgwk)) { dgwk = 0.0; }
dgwdcn(iref, iat) =
dgwk * zeta(ga, gi, refq(iref,iat) + zi, q(iat) + zi);
dgwk * zeta(ga, gi, refq(iref, iat) + zi, q(iat) + zi);
}
}
} else {
Expand Down Expand Up @@ -165,7 +166,7 @@ int TD4Model::weight_references(
}

gwvec(iref, iat) =
gwk * zeta(ga, gi, refq(iref,iat) + zi, q(iat) + zi);
gwk * zeta(ga, gi, refq(iref, iat) + zi, q(iat) + zi);
}
}
}
Expand Down Expand Up @@ -257,8 +258,7 @@ int TD4Model::get_atomic_c6(
return EXIT_SUCCESS;
}

int TD4Model::set_refq_eeq(const TMolecule &mol, TMatrix<double> &refq)
const {
int TD4Model::set_refq_eeq(const TMolecule &mol, TMatrix<double> &refq) const {
int izp;
for (int iat = 0; iat != mol.NAtoms; iat++) {
izp = mol.at(iat);
Expand Down
26 changes: 14 additions & 12 deletions test/test_grad.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -48,17 +48,19 @@ int test_numgrad(TMolecule &mol, const int charge, const dparam &par) {
mol.xyz(i, c) += step;
get_dispersion(mol, charge, d4, par, cutoff, er, nullptr);

mol.xyz(i, c) = mol.xyz(i, c) - 2*step;
mol.xyz(i, c) = mol.xyz(i, c) - 2 * step;
get_dispersion(mol, charge, d4, par, cutoff, el, nullptr);

mol.xyz(i, c) = mol.xyz(i, c) + step;
numgrad(i, c) = 0.5 * (er - el) / step;
}
}
}

// analytical gradient
double* d4grad = new double[3*mol.NAtoms];
for (int i = 0; i < 3*mol.NAtoms; i++) d4grad[i] = 0.0;
double *d4grad = new double[3 * mol.NAtoms];
for (int i = 0; i < 3 * mol.NAtoms; i++) {
d4grad[i] = 0.0;
}
info = get_dispersion(mol, charge, d4, par, cutoff, energy, d4grad);
if (!info == EXIT_SUCCESS) return info;

Expand All @@ -68,8 +70,8 @@ int test_numgrad(TMolecule &mol, const int charge, const dparam &par) {
// compare against numerical gradient
for (int i = 0; i < mol.NAtoms; i++) {
for (int c = 0; c < 3; c++) {
if (check(d4grad[3*i + c], numgrad(i, c), thr) != EXIT_SUCCESS) {
print_fail("Gradient mismatch", d4grad[3*i + c], numgrad(i, c));
if (check(d4grad[3 * i + c], numgrad(i, c), thr) != EXIT_SUCCESS) {
print_fail("Gradient mismatch", d4grad[3 * i + c], numgrad(i, c));
return EXIT_FAILURE;
}
}
Expand All @@ -80,14 +82,14 @@ int test_numgrad(TMolecule &mol, const int charge, const dparam &par) {
return EXIT_SUCCESS;
}

int is_trans_invar(const TMolecule& mol, double gradient[]) {
int is_trans_invar(const TMolecule &mol, double gradient[]) {
double xsum{0.0};
double ysum{0.0};
double zsum{0.0};
for (int i = 0; i < mol.NAtoms; i++) {
xsum += gradient[3*i];
ysum += gradient[3*i + 1];
zsum += gradient[3*i + 2];
xsum += gradient[3 * i];
ysum += gradient[3 * i + 1];
zsum += gradient[3 * i + 2];
}

if (check(xsum, 0.0) != EXIT_SUCCESS) {
Expand All @@ -106,7 +108,6 @@ int is_trans_invar(const TMolecule& mol, double gradient[]) {
return EXIT_SUCCESS;
}


int test_pbed4_mb01() {
// PBE-D4(EEQ) parameters
dparam par;
Expand All @@ -120,7 +121,8 @@ int test_pbed4_mb01() {
// assemble molecule
int charge = mb16_43_01_charge;
TMolecule mol;
int info = get_molecule(mb16_43_01_n, mb16_43_01_atoms, mb16_43_01_coord, mol);
int info =
get_molecule(mb16_43_01_n, mb16_43_01_atoms, mb16_43_01_coord, mol);
if (!info == EXIT_SUCCESS) return info;

return test_numgrad(mol, charge, par);
Expand Down
Loading