diff --git a/.gitignore b/.gitignore index c8e3819..3e1f83b 100644 --- a/.gitignore +++ b/.gitignore @@ -35,3 +35,8 @@ /build*/ /install*/ /_*/ + +# project-specific +/scripts/old/* +/scripts/*.h +/scripts/*.cpp diff --git a/include/dftd_ncoord.h b/include/dftd_ncoord.h index 07d71d9..6e112d6 100644 --- a/include/dftd_ncoord.h +++ b/include/dftd_ncoord.h @@ -200,4 +200,4 @@ extern inline double log_cn_cut(const double cn_max, const double cn); extern inline double dlog_cn_cut(const double cn_max, const double cn); -}; // namespace dftd4 +} // namespace dftd4 diff --git a/include/dftd_parameters.h b/include/dftd_parameters.h index b43dbec..4e9d235 100644 --- a/include/dftd_parameters.h +++ b/include/dftd_parameters.h @@ -63,10 +63,10 @@ static const double zeff[119]{ * These values are actually never used in the code. Only r4r2 is used. */ static const double r2r4[119]{ - 0.0, // dummy - 8.0589, 3.4698, // H,He + 0.0, // dummy + 8.0589, 3.4698, // H,He 29.0974, 14.8517, 11.8799, 7.8715, 5.5588, 4.7566, 3.8025, - 3.1036, // Li-Ne + 3.1036, // Li-Ne 26.1552, 17.2304, 17.7210, 12.7442, 9.5361, 8.1652, 6.7463, 5.6004, // Na-Ar 29.2012, 22.3934, // K,Ca diff --git a/scripts/README.md b/scripts/README.md new file mode 100644 index 0000000..db51510 --- /dev/null +++ b/scripts/README.md @@ -0,0 +1,3 @@ +# Scripts + +Scripts for creating the ORCA source files. This mainly concatenates the corresponding files and adds headers and footers. diff --git a/scripts/make-qcdftd4.sh b/scripts/make-qcdftd4.sh new file mode 100755 index 0000000..23db57b --- /dev/null +++ b/scripts/make-qcdftd4.sh @@ -0,0 +1,317 @@ +#!/bin/bash + +INCLUDE="../include" +SRC="../src" + +cat > qcdftd4param.h << 'EOT' +/* + * This is the DFT-D4 equivalent of copyc6, since this is an *empirical* + * dispersion correction we can hardly avoid having this big, scary blob + * of data lying around somewhere. + * + * Files like this are usually computer generated so avoid modifying it, + * but have the appropriate tool generate it from dftd4. + * + * Responsible for this mess: + * Sebastian Ehlert (SAW190521) + * + * Update/Fix for parameters from periodic extension by: + * Marvin Friede (MF161222) + * + * Extension for Fr, Ra and Actinides by: + * Marvin Friede (MF121223) + */ + +#ifndef QCDFTD4PARAM_H +#define QCDFTD4PARAM_H + +EOT + +sed -n '/namespace/,$p' "${INCLUDE}/dftd_parameters.h" >> qcdftd4param.h + +# Append the Footer: +printf "\n#endif // QCDFTD4PARAM_H\n" >> qcdftd4param.h + + +################################################################## +################################################################## + +dosed(){ + sed -n '/namespace dftd4 {/,/} \/\/ namespace dftd4/p' $1 | \ + sed '/namespace dftd4 {/d; /} \/\/ namespace dftd4/d' | \ + awk '/./,EOF' > $2 +} + +dosed "${INCLUDE}/dftd_cutoff.h" cutoff.txt +dosed "${INCLUDE}/dftd_model.h" model.txt +dosed "${INCLUDE}/dftd_dispersion.h" dispersion.txt +dosed "${INCLUDE}/damping/dftd_rational.h" rational.txt +dosed "${INCLUDE}/damping/dftd_atm.h" atm.txt + +cat > qcdftd4.h << 'EOT' +/* + * D4(EEQ)-ATM implementation + * + * Sebastian Ehlert (SAW190521) + * Marvin Friede (MF161222) + */ + +#ifndef __QCDFTD4_H +#define __QCDFTD4_H + +#include "qcinpdat.h" // TGeomInput class +#include "qcmat1.h" // TRVector and TRMatrix class + +namespace dftd4 { + +EOT + +###### + +cat >> qcdftd4.h << 'EOT' +/* -------------------------------------------------------------------------- + + Cutoff (dftd_cutoff.h) + https://github.com/dftd4/cpp-d4/blob/main/include/dftd_cutoff.h + + -------------------------------------------------------------------------- */ + +EOT + +cat cutoff.txt >> qcdftd4.h +rm cutoff.txt + +###### + +cat >> qcdftd4.h << 'EOT' +/* -------------------------------------------------------------------------- + + D4 model (dftd_model.h) + https://github.com/dftd4/cpp-d4/blob/main/include/dftd_model.h + + -------------------------------------------------------------------------- */ + +EOT + +cat model.txt >> qcdftd4.h +rm model.txt + +###### + +cat >> qcdftd4.h << 'EOT' +/* -------------------------------------------------------------------------- + + Dispersion (dftd_dispersion.h) + https://github.com/dftd4/cpp-d4/blob/main/include/dftd_dispersion.h + + -------------------------------------------------------------------------- */ + +EOT + +cat dispersion.txt >> qcdftd4.h +rm dispersion.txt + +###### + +cat >> qcdftd4.h << 'EOT' +/* -------------------------------------------------------------------------- + + Dispersion (damping/dftd_rational.h) + https://github.com/dftd4/cpp-d4/blob/main/include/damping/dftd_rational.h + + -------------------------------------------------------------------------- */ + +EOT + +cat rational.txt >> qcdftd4.h +rm rational.txt + +###### + +cat >> qcdftd4.h << 'EOT' +/* -------------------------------------------------------------------------- + + Dispersion (damping/dftd_atm.h) + https://github.com/dftd4/cpp-d4/blob/main/include/damping/dftd_atm.h + + -------------------------------------------------------------------------- */ + +EOT + +cat atm.txt >> qcdftd4.h +rm atm.txt + +###### + +cat >> qcdftd4.h << 'EOT' +} // namespace dftd4 + +/* -------------------------------------------------------------------------------------- +// Calculates the EEQ charges according to the D4 paper +// https://doi.org/10.1063/1.5090222 +// Bernardo de Souza, 14/09/2023 + -------------------------------------------------------------------------------------- */ +int CalcEEQCharges(TRMatrix &XYZ, TIVector &ATNO, int NAtoms, int totalcharge, TRVector &q, bool printerror=true); + +#endif // __QCDFTD4_H + +EOT + + +################################################################## +################################################################## + + +dosed "${SRC}/dftd_cutoff.cpp" cutoff.txt +dosed "${SRC}/dftd_model.cpp" model.txt +dosed "${SRC}/dftd_dispersion.cpp" dispersion.txt +dosed "${SRC}/damping/dftd_rational.cpp" rational.txt +dosed "${SRC}/damping/dftd_atm.cpp" atm.txt + +cat > qcdftd4.cpp << 'EOT' +/* + * D4(EEQ)-ATM implementation + * + * Sebastian Ehlert (SAW190521) + * Marvin Friede (MF161222) + */ +#include +#include + +#include "qcinpdat.h" // TGeomInput class +#include "qcmat2.h" // BLAS routines +#include "qcmath.h" // TVector and TMatrix class + +// we cannot avoid one big, scary parameter file... +#include "qcdftd4param.h" + +#include "qceeq.h" +#include "qcncoord.h" + +// always include self +#include "qcdftd4.h" + +/* -------------------------------------------------------------------------------------- +// Calculates the EEQ charges according to the D4 paper +// https://doi.org/10.1063/1.5090222 +// Bernardo de Souza, 14/09/2023 + -------------------------------------------------------------------------------------- */ +int CalcEEQCharges(TRMatrix &XYZ, TIVector &ATNO, int NAtoms, int totalcharge, TRVector &q, bool printerror){ + + // initialize the charges to zero + q.Init(); + + // check atomic numbers to guarantee we have all parameters + // ghost atoms (ATNO=0) will have no charge + for (int i=0;i86){ + if (printerror) printMessage("Atomic number %d detected. EEQ charges can not be calculated.\n",ATNO(i)); + return 1; + } + + TGeomInput molecule; + molecule.NAtoms=NAtoms; + molecule.CC.CopyMat(XYZ); + molecule.ATNO.CopyVec(ATNO); + TIVector Index(NAtoms); + for (int i=0;i> qcdftd4.cpp << 'EOT' +/* -------------------------------------------------------------------------- + + Cutoff (dftd_cutoff.cpp) + https://github.com/dftd4/cpp-d4/blob/main/src/dftd_cutoff.cpp + + -------------------------------------------------------------------------- */ + +EOT + +cat cutoff.txt >> qcdftd4.cpp +rm cutoff.txt + +###### + +cat >> qcdftd4.cpp << 'EOT' +/* -------------------------------------------------------------------------- + + D4 model (dftd_model.cpp) + https://github.com/dftd4/cpp-d4/blob/main/src/dftd_model.cpp + + -------------------------------------------------------------------------- */ + +EOT + +cat model.txt >> qcdftd4.cpp +rm model.txt + +###### + +cat >> qcdftd4.cpp << 'EOT' +/* -------------------------------------------------------------------------- + + Dispersion (dftd_dispersion.cpp) + https://github.com/dftd4/cpp-d4/blob/main/src/dftd_dispersion.cpp + + -------------------------------------------------------------------------- */ + +EOT + +cat dispersion.txt >> qcdftd4.cpp +rm dispersion.txt + +###### + +cat >> qcdftd4.cpp << 'EOT' +/* -------------------------------------------------------------------------- + + Dispersion (damping/dftd_rational.cpp) + https://github.com/dftd4/cpp-d4/blob/main/src/damping/dftd_rational.cpp + + -------------------------------------------------------------------------- */ + +EOT + +cat rational.txt >> qcdftd4.cpp +rm rational.txt + +###### + +cat >> qcdftd4.cpp << 'EOT' +/* -------------------------------------------------------------------------- + + Dispersion (damping/dftd_atm.cpp) + https://github.com/dftd4/cpp-d4/blob/main/src/damping/dftd_atm.cpp + + -------------------------------------------------------------------------- */ + +EOT + +cat atm.txt >> qcdftd4.cpp +rm atm.txt + +###### + +cat >> qcdftd4.cpp << 'EOT' +} // namespace dftd4 +EOT + + +########################################################################## + + +sed -i "s/TMolecule/TGeomInput/" *.h +sed -i "s/TMolecule/TGeomInput/" *.cpp diff --git a/scripts/make-qcdftd4param.sh b/scripts/make-qcdftd4param.sh new file mode 100755 index 0000000..76467ce --- /dev/null +++ b/scripts/make-qcdftd4param.sh @@ -0,0 +1,34 @@ +#!/bin/bash + +INCLUDE="../include" + +cat > qcdftd4param.h << 'EOT' +/* + * This is the DFT-D4 equivalent of copyc6, since this is an *empirical* + * dispersion correction we can hardly avoid having this big, scary blob + * of data lying around somewhere. + * + * Files like this are usually computer generated so avoid modifying it, + * but have the appropriate tool generate it from dftd4. + * + * Responsible for this mess: + * Sebastian Ehlert (SAW190521) + * + * Update/Fix for parameters from periodic extension by: + * Marvin Friede (MF161222) + * + * Extension for Fr, Ra and Actinides by: + * Marvin Friede (MF121223) + */ + +#ifndef QCDFTD4PARAM_H +#define QCDFTD4PARAM_H + +EOT + +sed -n '/namespace/,$p' "${INCLUDE}/dftd_parameters.h" >> qcdftd4param.h + +printf "\n#endif // QCDFTD4PARAM_H\n" >> qcdftd4param.h + + +sed -i "s/TMolecule/TGeomInput/" qcdftd4param.h diff --git a/scripts/make-qceeq.sh b/scripts/make-qceeq.sh new file mode 100755 index 0000000..4a6cf63 --- /dev/null +++ b/scripts/make-qceeq.sh @@ -0,0 +1,91 @@ +#!/bin/bash + +INCLUDE="../include" +SRC="../src" + +dosed(){ + sed -n '/namespace dftd4 {/,/} \/\/ namespace dftd4/p' $1 | \ + sed '/namespace dftd4 {/d; /} \/\/ namespace dftd4/d' | \ + awk '/./,EOF' > $2 +} + +################################################################## + +dosed "${INCLUDE}/dftd_eeq.h" eeq.txt + +cat > qceeq.h << 'EOT' +/* + * Electronegativity equilibration (EEQ) model for DFT-D4. + * Contains only the essential parts for DFT-D4. + * + * Sebastian Ehlert (SAW190521) + * Marvin Friede (MF161222) + * + * Extension for Fr, Ra and Actinides by: + * Marvin Friede (MF121223) + */ + +#ifndef QCEEQ_H +#define QCEEQ_H + +#include "qcinpdat.h" // TGeomInput class +#include "qcmat1.h" // TRVector and TRMatrix class + +namespace dftd4 { + +EOT + +cat eeq.txt >> qceeq.h +rm eeq.txt + +cat >> qceeq.h << 'EOT' +} // namespace dftd4 + +#endif // QCEEQ_H +EOT + +###### + +dosed "${SRC}/dftd_eeq.cpp" eeq.txt + +cat > qceeq.cpp << 'EOT' +/* + * Electronegativity equilibration (EEQ) model for DFT-D4. + * This implementation contains only the essential parts for DFT-D4. + * + * Sebastian Ehlert (SAW190521) + * Marvin Friede (MF161222) + * + * Extension for Fr, Ra and Actinides by: + * Marvin Friede (MF121223) + */ + +#include + +#include "qcinpdat.h" // TGeomInput class +#include "qclineq.h" // Linear Algebra routines +#include "qcmat1.h" // TRVector and TRMatrix class +#include "qcmat2.h" // BLAS routines + +// CN-related logic +#include "qcncoord.h" + +// always include self +#include "qceeq.h" + +// wrap everything in the dftd namespace to keep it nicely confined +namespace dftd4 { + +EOT + +cat eeq.txt >> qceeq.cpp +rm eeq.txt + +cat >> qceeq.cpp << 'EOT' +} // namespace dftd4 +EOT + +########################################################################## + +sed -i "s/TMolecule/TGeomInput/" *.h +sed -i "s/TMolecule/TGeomInput/" *.cpp diff --git a/scripts/make-qcncoord.sh b/scripts/make-qcncoord.sh new file mode 100755 index 0000000..c0a18f3 --- /dev/null +++ b/scripts/make-qcncoord.sh @@ -0,0 +1,89 @@ +#!/bin/bash + +INCLUDE="../include" +SRC="../src" + +dosed(){ + sed -n '/namespace dftd4 {/,/} \/\/ namespace dftd4/p' $1 | \ + sed '/namespace dftd4 {/d; /} \/\/ namespace dftd4/d' | \ + awk '/./,EOF' > $2 +} + +################################################################## + +dosed "${INCLUDE}/dftd_ncoord.h" ncoord.txt + +cat > qcncoord.h << 'EOT' +/* + * Definition of the coordination numbers used in DFT-D4 and the + * electronegativity equilibration (EEQ) model. + * + * Sebastian Ehlert (SAW190521) + * Marvin Friede (MF161222) + * + * Extension for Fr, Ra and Actinides by: + * Marvin Friede (MF121223) + */ + +#ifndef QCNCOORD_H +#define QCNCOORD_H + +#include "qcinpdat.h" // TGeomInput class +#include "qcmat1.h" // TRVector and TRMatrix class + +namespace dftd4 { + +EOT + +cat ncoord.txt >> qcncoord.h +rm ncoord.txt + +cat >> qcncoord.h << 'EOT' +} // namespace dftd4 + +#endif // QCNCOORD_H +EOT + +###### + +dosed "${SRC}/dftd_ncoord.cpp" ncoord.txt + +cat > qcncoord.cpp << 'EOT' +/* + * Definition of the coordination numbers used in DFT-D4 and the + * electronegativity equilibration (EEQ) model + * + * This module works on a distance matrix to avoid recalculating + * the distances every time. + * + * Sebastian Ehlert (SAW190521) + * Marvin Friede (MF161222) + * + * Extension for Fr, Ra and Actinides by: + * Marvin Friede (MF121223) + */ + +#include // erf, exp, log, sqrt + +#include "qcinpdat.h" // TGeomInput class +#include "qcmat1.h" // TRVector and TRMatrix class + +// always include self +#include "qcncoord.h" + +// wrap everything in the dftd namespace to keep it nicely confined +namespace dftd4 { + +EOT + +cat ncoord.txt >> qcncoord.cpp +rm ncoord.txt + +cat >> qcncoord.cpp << 'EOT' +} // namespace dftd4 +EOT + +########################################################################## + +sed -i "s/TMolecule/TGeomInput/" *.h +sed -i "s/TMolecule/TGeomInput/" *.cpp diff --git a/src/damping/dftd_atm.cpp b/src/damping/dftd_atm.cpp index ff1ed0a..f16b5f0 100644 --- a/src/damping/dftd_atm.cpp +++ b/src/damping/dftd_atm.cpp @@ -301,15 +301,15 @@ int get_atm_dispersion_derivs( dgyjk = c9ijk * (-dang * fdmp + ang * dfdmp) / r2jk * yjk; dgzjk = c9ijk * (-dang * fdmp + ang * dfdmp) / r2jk * zjk; - gradient(3 * iat + 0) += -dgxij - dgxik; - gradient(3 * iat + 1) += -dgyij - dgyik; - gradient(3 * iat + 2) += -dgzij - dgzik; - gradient(3 * jat + 0) += dgxij - dgxjk; - gradient(3 * jat + 1) += dgyij - dgyjk; - gradient(3 * jat + 2) += dgzij - dgzjk; - gradient(3 * kat + 0) += dgxik + dgxjk; - gradient(3 * kat + 1) += dgyik + dgyjk; - gradient(3 * kat + 2) += dgzik + dgzjk; + gradient(3 * ii + 0) += -dgxij - dgxik; + gradient(3 * ii + 1) += -dgyij - dgyik; + gradient(3 * ii + 2) += -dgzij - dgzik; + gradient(3 * jj + 0) += dgxij - dgxjk; + gradient(3 * jj + 1) += dgyij - dgyjk; + gradient(3 * jj + 2) += dgzij - dgzjk; + gradient(3 * kk + 0) += dgxik + dgxjk; + gradient(3 * kk + 1) += dgyik + dgyjk; + gradient(3 * kk + 2) += dgzik + dgzjk; dEdcn(ii) -= e * 0.5 * (dc6dcn(ii, jj) / c6ij + dc6dcn(ii, kk) / c6ik); dEdcn(jj) -= e * 0.5 * (dc6dcn(jj, ii) / c6ij + dc6dcn(jj, kk) / c6jk); diff --git a/src/damping/dftd_rational.cpp b/src/damping/dftd_rational.cpp index 906c7d2..23f10b6 100644 --- a/src/damping/dftd_rational.cpp +++ b/src/damping/dftd_rational.cpp @@ -18,6 +18,8 @@ #include +#include "damping/dftd_atm.h" +#include "damping/dftd_rational.h" #include "dftd_cblas.h" #include "dftd_dispersion.h" #include "dftd_eeq.h" @@ -25,9 +27,6 @@ #include "dftd_matrix.h" #include "dftd_ncoord.h" #include "dftd_parameters.h" -#include "damping/dftd_atm.h" -#include "damping/dftd_rational.h" - namespace dftd4 { @@ -119,8 +118,8 @@ int get_dispersion2_energy( } e = -c6ij * edisp * 0.5; - energy(iat) += e; - energy(jat) += e; + energy(ii) += e; + energy(jj) += e; } } @@ -164,7 +163,7 @@ int get_dispersion2_derivs( r4r2ij = 3.0 * r4r2[izp] * r4r2[jzp]; r0ij = par.a1 * sqrt(r4r2ij) + par.a2; - c6ij = c6(iat, jat); + c6ij = c6(ii, jj); t6 = fdmpr_bj(6, r, r0ij); t8 = fdmpr_bj(8, r, r0ij); @@ -199,12 +198,12 @@ int get_dispersion2_derivs( dEdq(ii) -= dc6dq(ii, jj) * edisp; dEdq(jj) -= dc6dq(jj, ii) * edisp; - gradient(3 * iat) += dgx; - gradient(3 * iat + 1) += dgy; - gradient(3 * iat + 2) += dgz; - gradient(3 * jat) -= dgx; - gradient(3 * jat + 1) -= dgy; - gradient(3 * jat + 2) -= dgz; + gradient(3 * ii) += dgx; + gradient(3 * ii + 1) += dgy; + gradient(3 * ii + 2) += dgz; + gradient(3 * jj) -= dgx; + gradient(3 * jj + 1) -= dgy; + gradient(3 * jj + 2) -= dgz; } } return EXIT_SUCCESS; diff --git a/src/dftd_dispersion.cpp b/src/dftd_dispersion.cpp index 2f1edae..0bba0cd 100644 --- a/src/dftd_dispersion.cpp +++ b/src/dftd_dispersion.cpp @@ -65,9 +65,9 @@ int get_dispersion( cn.NewVector(nat); q.NewVector(nat); if (lgrad) { - dcndr.NewMatrix(3 * mol.NAtoms, nat); - dqdr.NewMatrix(3 * mol.NAtoms, nat); - gradient.NewVector(3 * mol.NAtoms); + dcndr.NewMatrix(3 * nat, nat); + dqdr.NewMatrix(3 * nat, nat); + gradient.NewVector(3 * nat); } // calculate partial charges from EEQ model @@ -138,8 +138,7 @@ int get_dispersion( ); if (info != EXIT_SUCCESS) return info; - if (lgrad) { BLAS_Add_Mat_x_Vec(gradient, dqdr, dEdq, false, 1.0); } - + if (lgrad) BLAS_Add_Mat_x_Vec(gradient, dqdr, dEdq, false, 1.0); dqdr.DelMat(); // ---------------------------- @@ -224,9 +223,15 @@ int get_dispersion( // write to input gradient if (lgrad) { - for (int i = 0; i != 3 * mol.NAtoms; i++) { - GRAD[i] = gradient(i); + for (int i = 0, ii = 0; i != mol.NAtoms; i++) { + ii = realIdx(i); + if (ii < 0) continue; + + GRAD[3 * i] = gradient(3 * ii); + GRAD[3 * i + 1] = gradient(3 * ii + 1); + GRAD[3 * i + 2] = gradient(3 * ii + 2); } + gradient.DelVec(); } diff --git a/src/dftd_eeq.cpp b/src/dftd_eeq.cpp index 24deec5..c0370ca 100644 --- a/src/dftd_eeq.cpp +++ b/src/dftd_eeq.cpp @@ -34,103 +34,90 @@ namespace dftd4 { constexpr int MAXELEMENT = 104; // 103 + dummy static const double xi[MAXELEMENT]{ - +0.00000000, - +1.23695041, +1.26590957, +0.54341808, +0.99666991, +1.26691604, - +1.40028282, +1.55819364, +1.56866440, +1.57540015, +1.15056627, - +0.55936220, +0.72373742, +1.12910844, +1.12306840, +1.52672442, - +1.40768172, +1.48154584, +1.31062963, +0.40374140, +0.75442607, - +0.76482096, +0.98457281, +0.96702598, +1.05266584, +0.93274875, - +1.04025281, +0.92738624, +1.07419210, +1.07900668, +1.04712861, - +1.15018618, +1.15388455, +1.36313743, +1.36485106, +1.39801837, - +1.18695346, +0.36273870, +0.58797255, +0.71961946, +0.96158233, - +0.89585296, +0.81360499, +1.00794665, +0.92613682, +1.09152285, - +1.14907070, +1.13508911, +1.08853785, +1.11005982, +1.12452195, - +1.21642129, +1.36507125, +1.40340000, +1.16653482, +0.34125098, - +0.58884173, +0.68441115, +0.56999999, +0.56999999, +0.56999999, - +0.56999999, +0.56999999, +0.56999999, +0.56999999, +0.56999999, - +0.56999999, +0.56999999, +0.56999999, +0.56999999, +0.56999999, - +0.56999999, +0.87936784, +1.02761808, +0.93297476, +1.10172128, - +0.97350071, +1.16695666, +1.23997927, +1.18464453, +1.14191734, - +1.12334192, +1.01485321, +1.12950808, +1.30804834, +1.33689961, - +1.27465977, +1.06598299, +0.68184178, +1.25994649, +0.40030316, - +0.29458157, +0.75460678, +0.42330605, +0.53613714, +0.58870880, - +0.50349105, +0.02647510, +0.45350118, +0.63133032, +0.45995218, - +1.01557723, +1.17163176, +0.23501195, + +0.00000000, +1.23695041, +1.26590957, +0.54341808, +0.99666991, +1.26691604, + +1.40028282, +1.55819364, +1.56866440, +1.57540015, +1.15056627, +0.55936220, + +0.72373742, +1.12910844, +1.12306840, +1.52672442, +1.40768172, +1.48154584, + +1.31062963, +0.40374140, +0.75442607, +0.76482096, +0.98457281, +0.96702598, + +1.05266584, +0.93274875, +1.04025281, +0.92738624, +1.07419210, +1.07900668, + +1.04712861, +1.15018618, +1.15388455, +1.36313743, +1.36485106, +1.39801837, + +1.18695346, +0.36273870, +0.58797255, +0.71961946, +0.96158233, +0.89585296, + +0.81360499, +1.00794665, +0.92613682, +1.09152285, +1.14907070, +1.13508911, + +1.08853785, +1.11005982, +1.12452195, +1.21642129, +1.36507125, +1.40340000, + +1.16653482, +0.34125098, +0.58884173, +0.68441115, +0.56999999, +0.56999999, + +0.56999999, +0.56999999, +0.56999999, +0.56999999, +0.56999999, +0.56999999, + +0.56999999, +0.56999999, +0.56999999, +0.56999999, +0.56999999, +0.56999999, + +0.87936784, +1.02761808, +0.93297476, +1.10172128, +0.97350071, +1.16695666, + +1.23997927, +1.18464453, +1.14191734, +1.12334192, +1.01485321, +1.12950808, + +1.30804834, +1.33689961, +1.27465977, +1.06598299, +0.68184178, +1.25994649, + +0.40030316, +0.29458157, +0.75460678, +0.42330605, +0.53613714, +0.58870880, + +0.50349105, +0.02647510, +0.45350118, +0.63133032, +0.45995218, +1.01557723, + +1.17163176, +0.23501195, }; -static const double gam[MAXELEMENT]{ // eta (hardness) - +0.00000000, - -0.35015861, +1.04121227, +0.09281243, +0.09412380, +0.26629137, - +0.19408787, +0.05317918, +0.03151644, +0.32275132, +1.30996037, - +0.24206510, +0.04147733, +0.11634126, +0.13155266, +0.15350650, - +0.15250997, +0.17523529, +0.28774450, +0.42937314, +0.01896455, - +0.07179178, -0.01121381, -0.03093370, +0.02716319, -0.01843812, - -0.15270393, -0.09192645, -0.13418723, -0.09861139, +0.18338109, - +0.08299615, +0.11370033, +0.19005278, +0.10980677, +0.12327841, - +0.25345554, +0.58615231, +0.16093861, +0.04548530, -0.02478645, - +0.01909943, +0.01402541, -0.03595279, +0.01137752, -0.03697213, - +0.08009416, +0.02274892, +0.12801822, -0.02078702, +0.05284319, - +0.07581190, +0.09663758, +0.09547417, +0.07803344, +0.64913257, - +0.15348654, +0.05054344, +0.11000000, +0.11000000, +0.11000000, - +0.11000000, +0.11000000, +0.11000000, +0.11000000, +0.11000000, - +0.11000000, +0.11000000, +0.11000000, +0.11000000, +0.11000000, - +0.11000000, -0.02786741, +0.01057858, -0.03892226, -0.04574364, - -0.03874080, -0.03782372, -0.07046855, +0.09546597, +0.21953269, - +0.02522348, +0.15263050, +0.08042611, +0.01878626, +0.08715453, - +0.10500484, +0.10034731, +0.15801991, -0.00995585, +0.13654805, - +0.15968516, +0.12442023, +0.13340491, +0.11523968, +0.12165907, - +0.08181387, +0.65060023, +0.18634958, +0.16685501, +0.15080023, - +0.10094112, +0.10503808, +0.68963544, +static const double gam[MAXELEMENT]{ + // eta (hardness) + +0.00000000, -0.35015861, +1.04121227, +0.09281243, +0.09412380, +0.26629137, + +0.19408787, +0.05317918, +0.03151644, +0.32275132, +1.30996037, +0.24206510, + +0.04147733, +0.11634126, +0.13155266, +0.15350650, +0.15250997, +0.17523529, + +0.28774450, +0.42937314, +0.01896455, +0.07179178, -0.01121381, -0.03093370, + +0.02716319, -0.01843812, -0.15270393, -0.09192645, -0.13418723, -0.09861139, + +0.18338109, +0.08299615, +0.11370033, +0.19005278, +0.10980677, +0.12327841, + +0.25345554, +0.58615231, +0.16093861, +0.04548530, -0.02478645, +0.01909943, + +0.01402541, -0.03595279, +0.01137752, -0.03697213, +0.08009416, +0.02274892, + +0.12801822, -0.02078702, +0.05284319, +0.07581190, +0.09663758, +0.09547417, + +0.07803344, +0.64913257, +0.15348654, +0.05054344, +0.11000000, +0.11000000, + +0.11000000, +0.11000000, +0.11000000, +0.11000000, +0.11000000, +0.11000000, + +0.11000000, +0.11000000, +0.11000000, +0.11000000, +0.11000000, +0.11000000, + -0.02786741, +0.01057858, -0.03892226, -0.04574364, -0.03874080, -0.03782372, + -0.07046855, +0.09546597, +0.21953269, +0.02522348, +0.15263050, +0.08042611, + +0.01878626, +0.08715453, +0.10500484, +0.10034731, +0.15801991, -0.00995585, + +0.13654805, +0.15968516, +0.12442023, +0.13340491, +0.11523968, +0.12165907, + +0.08181387, +0.65060023, +0.18634958, +0.16685501, +0.15080023, +0.10094112, + +0.10503808, +0.68963544, }; -static const double kappa[MAXELEMENT]{ // kcn - +0.00000000, - +0.04916110, +0.10937243, -0.12349591, -0.02665108, -0.02631658, - +0.06005196, +0.09279548, +0.11689703, +0.15704746, +0.07987901, - -0.10002962, -0.07712863, -0.02170561, -0.04964052, +0.14250599, - +0.07126660, +0.13682750, +0.14877121, -0.10219289, -0.08979338, - -0.08273597, -0.01754829, -0.02765460, -0.02558926, -0.08010286, - -0.04163215, -0.09369631, -0.03774117, -0.05759708, +0.02431998, - -0.01056270, -0.02692862, +0.07657769, +0.06561608, +0.08006749, - +0.14139200, -0.05351029, -0.06701705, -0.07377246, -0.02927768, - -0.03867291, -0.06929825, -0.04485293, -0.04800824, -0.01484022, - +0.07917502, +0.06619243, +0.02434095, -0.01505548, -0.03030768, - +0.01418235, +0.08953411, +0.08967527, +0.07277771, -0.02129476, - -0.06188828, -0.06568203, -0.11000000, -0.11000000, -0.11000000, - -0.11000000, -0.11000000, -0.11000000, -0.11000000, -0.11000000, - -0.11000000, -0.11000000, -0.11000000, -0.11000000, -0.11000000, - -0.11000000, -0.03585873, -0.03132400, -0.05902379, -0.02827592, - -0.07606260, -0.02123839, +0.03814822, +0.02146834, +0.01580538, - -0.00894298, -0.05864876, -0.01817842, +0.07721851, +0.07936083, - +0.05849285, +0.00013506, -0.00020631, +0.00328823, +0.00419390, - +0.00429264, +0.00193300, +0.00478177, +0.00264040, +0.00418168, - +0.00399258, -0.00293781, -0.00195990, +0.00148017, -0.00011254, - +0.00023249, -0.00006144, -0.02459107, +static const double kappa[MAXELEMENT]{ + // kcn + +0.00000000, +0.04916110, +0.10937243, -0.12349591, -0.02665108, -0.02631658, + +0.06005196, +0.09279548, +0.11689703, +0.15704746, +0.07987901, -0.10002962, + -0.07712863, -0.02170561, -0.04964052, +0.14250599, +0.07126660, +0.13682750, + +0.14877121, -0.10219289, -0.08979338, -0.08273597, -0.01754829, -0.02765460, + -0.02558926, -0.08010286, -0.04163215, -0.09369631, -0.03774117, -0.05759708, + +0.02431998, -0.01056270, -0.02692862, +0.07657769, +0.06561608, +0.08006749, + +0.14139200, -0.05351029, -0.06701705, -0.07377246, -0.02927768, -0.03867291, + -0.06929825, -0.04485293, -0.04800824, -0.01484022, +0.07917502, +0.06619243, + +0.02434095, -0.01505548, -0.03030768, +0.01418235, +0.08953411, +0.08967527, + +0.07277771, -0.02129476, -0.06188828, -0.06568203, -0.11000000, -0.11000000, + -0.11000000, -0.11000000, -0.11000000, -0.11000000, -0.11000000, -0.11000000, + -0.11000000, -0.11000000, -0.11000000, -0.11000000, -0.11000000, -0.11000000, + -0.03585873, -0.03132400, -0.05902379, -0.02827592, -0.07606260, -0.02123839, + +0.03814822, +0.02146834, +0.01580538, -0.00894298, -0.05864876, -0.01817842, + +0.07721851, +0.07936083, +0.05849285, +0.00013506, -0.00020631, +0.00328823, + +0.00419390, +0.00429264, +0.00193300, +0.00478177, +0.00264040, +0.00418168, + +0.00399258, -0.00293781, -0.00195990, +0.00148017, -0.00011254, +0.00023249, + -0.00006144, -0.02459107, }; -static const double alp[MAXELEMENT]{ // rad - +0.00000000, - +0.55159092, +0.66205886, +0.90529132, +1.51710827, +2.86070364, - +1.88862966, +1.32250290, +1.23166285, +1.77503721, +1.11955204, - +1.28263182, +1.22344336, +1.70936266, +1.54075036, +1.38200579, - +2.18849322, +1.36779065, +1.27039703, +1.64466502, +1.58859404, - +1.65357953, +1.50021521, +1.30104175, +1.46301827, +1.32928147, - +1.02766713, +1.02291377, +0.94343886, +1.14881311, +1.47080755, - +1.76901636, +1.98724061, +2.41244711, +2.26739524, +2.95378999, - +1.20807752, +1.65941046, +1.62733880, +1.61344972, +1.63220728, - +1.60899928, +1.43501286, +1.54559205, +1.32663678, +1.37644152, - +1.36051851, +1.23395526, +1.65734544, +1.53895240, +1.97542736, - +1.97636542, +2.05432381, +3.80138135, +1.43893803, +1.75505957, - +1.59815118, +1.76401732, +1.63999999, +1.63999999, +1.63999999, - +1.63999999, +1.63999999, +1.63999999, +1.63999999, +1.63999999, - +1.63999999, +1.63999999, +1.63999999, +1.63999999, +1.63999999, - +1.63999999, +1.47055223, +1.81127084, +1.40189963, +1.54015481, - +1.33721475, +1.57165422, +1.04815857, +1.78342098, +2.79106396, - +1.78160840, +2.47588882, +2.37670734, +1.76613217, +2.66172302, - +2.82773085, +1.04059593, +0.60550051, +1.93854984, +0.50189075, - +0.58180664, +0.73094166, +0.49548126, +0.67685715, +0.59573917, - +0.35345732, +0.29902232, +0.49626064, +0.55816329, +0.49019371, - +1.05120718, +0.95651052, +0.35885251, +static const double alp[MAXELEMENT]{ + // rad + +0.00000000, +0.55159092, +0.66205886, +0.90529132, +1.51710827, +2.86070364, + +1.88862966, +1.32250290, +1.23166285, +1.77503721, +1.11955204, +1.28263182, + +1.22344336, +1.70936266, +1.54075036, +1.38200579, +2.18849322, +1.36779065, + +1.27039703, +1.64466502, +1.58859404, +1.65357953, +1.50021521, +1.30104175, + +1.46301827, +1.32928147, +1.02766713, +1.02291377, +0.94343886, +1.14881311, + +1.47080755, +1.76901636, +1.98724061, +2.41244711, +2.26739524, +2.95378999, + +1.20807752, +1.65941046, +1.62733880, +1.61344972, +1.63220728, +1.60899928, + +1.43501286, +1.54559205, +1.32663678, +1.37644152, +1.36051851, +1.23395526, + +1.65734544, +1.53895240, +1.97542736, +1.97636542, +2.05432381, +3.80138135, + +1.43893803, +1.75505957, +1.59815118, +1.76401732, +1.63999999, +1.63999999, + +1.63999999, +1.63999999, +1.63999999, +1.63999999, +1.63999999, +1.63999999, + +1.63999999, +1.63999999, +1.63999999, +1.63999999, +1.63999999, +1.63999999, + +1.47055223, +1.81127084, +1.40189963, +1.54015481, +1.33721475, +1.57165422, + +1.04815857, +1.78342098, +2.79106396, +1.78160840, +2.47588882, +2.37670734, + +1.76613217, +2.66172302, +2.82773085, +1.04059593, +0.60550051, +1.93854984, + +0.50189075, +0.58180664, +0.73094166, +0.49548126, +0.67685715, +0.59573917, + +0.35345732, +0.29902232, +0.49626064, +0.55816329, +0.49019371, +1.05120718, + +0.95651052, +0.35885251, }; static const double small = 1e-14; @@ -156,7 +143,7 @@ int get_charges( TMatrix dcndr; // Derivative of EEQ-CN cn.NewVec(nat); - if (lgrad) dcndr.NewMat(nat, 3 * mol.NAtoms); + if (lgrad) dcndr.NewMat(nat, 3 * nat); // get the EEQ coordination number info = get_ncoord_erf(mol, realIdx, dist, cutoff, cn, dcndr, lgrad); @@ -288,7 +275,7 @@ int get_damat_0d( gam = 1.0 / std::sqrt((alphai + alphaj)); arg = gam * gam * r2; - dtmp = 2.0 * gam * std::exp(-arg) / (sqrtpi * r2) - Amat(j, i) / r2; + dtmp = 2.0 * gam * std::exp(-arg) / (sqrtpi * r2) - Amat(jj, ii) / r2; dgx = dtmp * rx; dgy = dtmp * ry; dgz = dtmp * rz; @@ -300,12 +287,12 @@ int get_damat_0d( atrace(jj, 1) -= dgy * q(ii); atrace(jj, 2) -= dgz * q(ii); - dAmat(3 * i, jj) = dgx * q(ii); - dAmat(3 * i + 1, jj) = dgy * q(ii); - dAmat(3 * i + 2, jj) = dgz * q(i); - dAmat(3 * j, ii) = -dgx * q(jj); - dAmat(3 * j + 1, ii) = -dgy * q(jj); - dAmat(3 * j + 2, ii) = -dgz * q(jj); + dAmat(3 * ii, jj) = dgx * q(ii); + dAmat(3 * ii + 1, jj) = dgy * q(ii); + dAmat(3 * ii + 2, jj) = dgz * q(i); + dAmat(3 * jj, ii) = -dgx * q(jj); + dAmat(3 * jj + 1, ii) = -dgy * q(jj); + dAmat(3 * jj + 2, ii) = -dgz * q(jj); } } @@ -328,7 +315,6 @@ int eeq_chrgeq( int info = 0; int n = realIdx.Max() + 1; int m = n + 1; - int mm = m - 1; TMatrix Amat; // Coulomb matrix TVector xvec; // x (chi) vector @@ -336,7 +322,7 @@ int eeq_chrgeq( xvec.NewVec(m); TVector dxdcn; // Derivative of chi vector w.r.t. CN - if (lgrad) { dxdcn.NewVec(m); } + if (lgrad) dxdcn.NewVec(m); info = get_vrhs(mol, realIdx, charge, cn, xvec, dxdcn, lgrad); if (info != EXIT_SUCCESS) return info; @@ -353,15 +339,16 @@ int eeq_chrgeq( // solve: A Q = X (Eq.4) -> Q = Ainv X info = BLAS_InvertMatrix(Ainv); - - if (info != EXIT_SUCCESS) return info; - info = BLAS_Add_Mat_x_Vec(vrhs, Ainv, xvec, false, 1.0); if (info != EXIT_SUCCESS) return info; + // no return in ORCA + BLAS_Add_Mat_x_Vec(vrhs, Ainv, xvec, false, 1.0); + // if (info != EXIT_SUCCESS) return info; + // remove charge constraint (make vector smaller by one) and total charge qtotal = 0.0; - for (int i = 0, ii = 0; i != mm; i++) { - ii = realIdx(i); + for (int iat = 0, ii = 0; iat != mol.NAtoms; iat++) { + ii = realIdx(iat); if (ii < 0) continue; qvec(ii) = vrhs(ii); @@ -393,10 +380,8 @@ int eeq_chrgeq( // Gradient (note that the corresponding gradient flag in Fortran is `cpq`) if (lgrad) { - int ThreeN = 3 * mol.NAtoms; - TMatrix dAmat; - dAmat.NewMat(ThreeN, m); + dAmat.NewMat(3 * n, m); TMatrix atrace; atrace.NewMat(m, 3); @@ -407,17 +392,17 @@ int eeq_chrgeq( ii = realIdx(i); if (ii < 0) continue; - dAmat(3 * i, ii) += atrace(ii, 0); - dAmat(3 * i + 1, ii) += atrace(ii, 1); - dAmat(3 * i + 2, ii) += atrace(ii, 2); + dAmat(3 * ii, ii) += atrace(ii, 0); + dAmat(3 * ii + 1, ii) += atrace(ii, 1); + dAmat(3 * ii + 2, ii) += atrace(ii, 2); for (int j = 0, jj = 0; j != mol.NAtoms; j++) { jj = realIdx(j); if (jj < 0) continue; - dAmat(3 * j, ii) -= dcndr(jj, 3 * i) * dxdcn(ii); - dAmat(3 * j + 1, ii) -= dcndr(jj, 3 * i + 1) * dxdcn(ii); - dAmat(3 * j + 2, i) -= dcndr(jj, 3 * i + 2) * dxdcn(ii); + dAmat(3 * jj, ii) -= dcndr(jj, 3 * ii) * dxdcn(ii); + dAmat(3 * jj + 1, ii) -= dcndr(jj, 3 * ii + 1) * dxdcn(ii); + dAmat(3 * jj + 2, ii) -= dcndr(jj, 3 * ii + 2) * dxdcn(ii); } } @@ -434,8 +419,9 @@ int eeq_chrgeq( } } - info = BLAS_Add_Mat_x_Mat(dqdr, dAmat, A, false, false, -1.0); - if (info != EXIT_SUCCESS) return info; + // no return in ORCA + BLAS_Add_Mat_x_Mat(dqdr, dAmat, A, false, false, -1.0); + // if (info != EXIT_SUCCESS) return info; dAmat.DelMat(); } diff --git a/src/dftd_model.cpp b/src/dftd_model.cpp index 8df443d..e404518 100644 --- a/src/dftd_model.cpp +++ b/src/dftd_model.cpp @@ -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*/, @@ -126,9 +127,8 @@ int TD4Model::weight_references( } } - gwvec(iref, iat) = gwk * zeta(ga, gi, refq[izp][iref] + zi, q(ii) + zi); - dgwdq(iref, iat) = - gwk * dzeta(ga, gi, refq[izp][iref] + zi, q(ii) + zi); + gwvec(iref, ii) = gwk * zeta(ga, gi, refq[izp][iref] + zi, q(ii) + zi); + dgwdq(iref, ii) = gwk * dzeta(ga, gi, refq[izp][iref] + zi, q(ii) + zi); dgwk = norm * (dexpw - expw * dnorm * norm); if (is_exceptional(dgwk)) { dgwk = 0.0; } diff --git a/src/dftd_ncoord.cpp b/src/dftd_ncoord.cpp index c41800b..9cffa82 100644 --- a/src/dftd_ncoord.cpp +++ b/src/dftd_ncoord.cpp @@ -27,8 +27,8 @@ #include #include "dftd_geometry.h" -#include "dftd_ncoord.h" #include "dftd_matrix.h" +#include "dftd_ncoord.h" namespace dftd4 { @@ -383,7 +383,7 @@ int get_ncoord_d4( TMatrix &dcndr, bool lgrad ) { - if (lgrad) { return dncoord_d4(mol, realIdx, dist, cutoff, cn, dcndr); } + if (lgrad) return dncoord_d4(mol, realIdx, dist, cutoff, cn, dcndr); return ncoord_d4(mol, realIdx, dist, cutoff, cn); }; @@ -462,18 +462,18 @@ int dncoord_d4( cn(jj) += countf; dcountf = den * derf_count(kn, rr) / rcovij; - dcndr(3 * j, jj) += dcountf * rx; - dcndr(3 * j + 1, jj) += dcountf * ry; - dcndr(3 * j + 2, jj) += dcountf * rz; - dcndr(3 * j, ii) = dcountf * rx; - dcndr(3 * j + 1, ii) = dcountf * ry; - dcndr(3 * j + 2, ii) = dcountf * rz; - dcndr(3 * i, jj) = -dcountf * rx; - dcndr(3 * i + 1, jj) = -dcountf * ry; - dcndr(3 * i + 2, jj) = -dcountf * rz; - dcndr(3 * i, ii) += -dcountf * rx; - dcndr(3 * i + 1, ii) += -dcountf * ry; - dcndr(3 * i + 2, ii) += -dcountf * rz; + dcndr(3 * jj, jj) += dcountf * rx; + dcndr(3 * jj + 1, jj) += dcountf * ry; + dcndr(3 * jj + 2, jj) += dcountf * rz; + dcndr(3 * jj, ii) = dcountf * rx; + dcndr(3 * jj + 1, ii) = dcountf * ry; + dcndr(3 * jj + 2, ii) = dcountf * rz; + dcndr(3 * ii, jj) = -dcountf * rx; + dcndr(3 * ii + 1, jj) = -dcountf * ry; + dcndr(3 * ii + 2, jj) = -dcountf * rz; + dcndr(3 * ii, ii) += -dcountf * rx; + dcndr(3 * ii + 1, ii) += -dcountf * ry; + dcndr(3 * ii + 2, ii) += -dcountf * rz; } } return EXIT_SUCCESS; @@ -546,21 +546,21 @@ int dncoord_erf( cn(jj) += countf; dcountf = derf_count(kn, rr) / rcovij; - dcndr(jj, 3 * j) += dcountf * rx; - dcndr(jj, 3 * j + 1) += dcountf * ry; - dcndr(jj, 3 * j + 2) += dcountf * rz; + dcndr(jj, 3 * jj) += dcountf * rx; + dcndr(jj, 3 * jj + 1) += dcountf * ry; + dcndr(jj, 3 * jj + 2) += dcountf * rz; - dcndr(jj, 3 * i) += dcountf * rx; - dcndr(jj, 3 * i + 1) += dcountf * ry; - dcndr(jj, 3 * i + 2) += dcountf * rz; + dcndr(jj, 3 * ii) += dcountf * rx; + dcndr(jj, 3 * ii + 1) += dcountf * ry; + dcndr(jj, 3 * ii + 2) += dcountf * rz; - dcndr(ii, 3 * j) -= dcountf * rx; - dcndr(ii, 3 * j + 1) -= dcountf * ry; - dcndr(ii, 3 * j + 2) -= dcountf * rz; + dcndr(ii, 3 * jj) -= dcountf * rx; + dcndr(ii, 3 * jj + 1) -= dcountf * ry; + dcndr(ii, 3 * jj + 2) -= dcountf * rz; - dcndr(ii, 3 * i) -= dcountf * rx; - dcndr(ii, 3 * i + 1) -= dcountf * ry; - dcndr(ii, 3 * i + 2) -= dcountf * rz; + dcndr(ii, 3 * ii) -= dcountf * rx; + dcndr(ii, 3 * ii + 1) -= dcountf * ry; + dcndr(ii, 3 * ii + 2) -= dcountf * rz; } } diff --git a/test/molecules.h b/test/molecules.h index 05de81d..5168c8e 100644 --- a/test/molecules.h +++ b/test/molecules.h @@ -103,6 +103,38 @@ static const double water_dimer_coord[water_dimer_n * 3]{ -1.42952433403007, }; +// Water ghost +static const int water_ghost_n{6}; +static const int water_ghost_charge{0}; +static const char water_ghost_atoms[water_ghost_n][4]{ + "X", + "O", + "X", + "H", + "H", + "X", +}; +static const double water_ghost_coord[water_ghost_n * 3]{ + +2.36992046172956, // ghost + +0.05080150740831, // ghost + +0.00000000000000, // ghost + -3.24461780406248, // O + -0.10530990158270, // O + +0.00000000000000, // O + +3.07751820457461, // ghost + -0.81358303254468, // ghost + +1.42952433403007, // ghost + -3.84932865216183, // H1 + +1.60128210841514, // H1 + +0.00000000000000, // H1 + -1.43101041465447, // H2 + +0.08039197290338, // H2 + +0.00000000000000, // H2 + +3.07751820457461, // ghost + -0.81358303254468, // ghost + -1.42952433403007, // ghost +}; + // MB16_43: 01 static const int mb16_43_01_n{16}; static const int mb16_43_01_charge{0}; diff --git a/test/test_ghost.cpp b/test/test_ghost.cpp index 35aaa98..de75b8b 100644 --- a/test/test_ghost.cpp +++ b/test/test_ghost.cpp @@ -31,7 +31,13 @@ using namespace dftd4; -int test_water(const int n, const char atoms[][4], const double coord[]) { +int test_water( + const int n, + const char atoms[][4], + const double coord[], + const double ref_grad[], + const TIVector &realIdx +) { int info; // assemble molecule @@ -40,16 +46,6 @@ int test_water(const int n, const char atoms[][4], const double coord[]) { info = get_molecule(n, atoms, coord, mol); if (info != EXIT_SUCCESS) return info; - // dummy for masking of ghost atoms - TIVector realIdx; - realIdx.NewVec(mol.NAtoms); - realIdx(0) = 0; - realIdx(1) = 1; - realIdx(2) = 2; - realIdx(3) = -1; - realIdx(4) = -1; - realIdx(5) = -1; - // number of real atoms int nat = realIdx.Max() + 1; @@ -116,9 +112,6 @@ int test_water(const int n, const char atoms[][4], const double coord[]) { info = get_dispersion(mol, realIdx, charge, d4, par, cutoff, energy, nullptr); if (info != EXIT_SUCCESS) return info; - - - if (check(energy, ref, 1e-8) == EXIT_FAILURE) { print_fail("GHOST: Two-body Energy", energy, ref); return EXIT_FAILURE; @@ -153,11 +146,11 @@ int test_water(const int n, const char atoms[][4], const double coord[]) { d4grad[i] = 0.0; } - dcndr.NewMatrix(3 * mol.NAtoms, nat); + dcndr.NewMatrix(3 * nat, nat); info = get_ncoord_d4(mol, realIdx, dist, cutoff.cn, cn, dcndr, lgrad); if (info != EXIT_SUCCESS) return info; - dqdr.NewMatrix(3 * mol.NAtoms, nat); + dqdr.NewMatrix(3 * nat, nat); info = get_charges(mol, realIdx, dist, charge, cutoff.cn_eeq, q, dqdr, lgrad); if (info != EXIT_SUCCESS) return info; @@ -165,8 +158,8 @@ int test_water(const int n, const char atoms[][4], const double coord[]) { if (info != EXIT_SUCCESS) return info; for (int i = 0; i < 3 * mol.NAtoms; i++) { - if (check(d4grad[i], water_dimer_ref_grad[i], 1e-8) == EXIT_FAILURE) { - print_fail("GHOST: Gradient", d4grad[i], water_dimer_ref_grad[i]); + if (check(d4grad[i], ref_grad[i], 1e-8) == EXIT_FAILURE) { + print_fail("GHOST: Gradient", d4grad[i], ref_grad[i]); return EXIT_FAILURE; } } @@ -177,7 +170,39 @@ int test_water(const int n, const char atoms[][4], const double coord[]) { int test_ghost() { int info; - info = test_water(water_dimer_n, water_dimer_atoms, water_dimer_coord); + // dummy for masking of ghost atoms + TIVector realidx; + realidx.NewVec(water_dimer_n); + realidx(0) = 0; + realidx(1) = 1; + realidx(2) = 2; + realidx(3) = -1; + realidx(4) = -1; + realidx(5) = -1; + + info = test_water( + water_dimer_n, + water_dimer_atoms, + water_dimer_coord, + water_dimer_ref_grad, + realidx + ); + if (info != EXIT_SUCCESS) return info; + + realidx(0) = -1; + realidx(1) = 0; + realidx(2) = -1; + realidx(3) = 1; + realidx(4) = 2; + realidx(5) = -1; + + info = test_water( + water_ghost_n, + water_ghost_atoms, + water_ghost_coord, + water_ghost_ref_grad, + realidx + ); if (info != EXIT_SUCCESS) return info; return EXIT_SUCCESS; diff --git a/test/test_ghost.h b/test/test_ghost.h index 3d5c12f..3d66532 100644 --- a/test/test_ghost.h +++ b/test/test_ghost.h @@ -18,6 +18,10 @@ #ifndef TEST_GHOST_H #define TEST_GHOST_H +#include + +using namespace dftd4; + static const double water_dimer_ref_cn[3]{ +1.6104087812936088, +0.80554651093269469, @@ -52,7 +56,34 @@ static const double water_dimer_ref_grad[3 * 6]{ +0.0000000000000000E+00, }; -extern int test_water(const int n, const char atoms[][4], const double coord[]); +static const double water_ghost_ref_grad[3 * 6]{ + +0.0000000000000000E+00, // ghost + +0.0000000000000000E+00, // ghost + +0.0000000000000000E+00, // ghost + +3.6769324466783607E-05, // O + +5.8683292759696172E-05, // O + +0.0000000000000000E+00, // O + +0.0000000000000000E+00, // ghost + +0.0000000000000000E+00, // ghost + +0.0000000000000000E+00, // ghost + +4.8209580157792990E-06, // H1 + -4.4217699200268973E-05, // H1 + +0.0000000000000000E+00, // H1 + -4.1590282482562915E-05, // H2 + -1.4465593559427221E-05, // H2 + +0.0000000000000000E+00, // H2 + +0.0000000000000000E+00, // ghost + +0.0000000000000000E+00, // ghost + +0.0000000000000000E+00, // ghost +}; + +extern int test_water( + const int n, + const char atoms[][4], + const double coord[], + const double ref_grad[], + const TIVector &realIdx +); extern int test_ghost(void);