diff --git a/external/CxxThreadPool b/external/CxxThreadPool index ca5058b..e0f2fec 160000 --- a/external/CxxThreadPool +++ b/external/CxxThreadPool @@ -1 +1 @@ -Subproject commit ca5058bd27ec88cb7291781068a8635eab74bd32 +Subproject commit e0f2fec388d0fb1f58f6da1894729a04aa407cbb diff --git a/external/tblite b/external/tblite index c5844a4..9075e2e 160000 --- a/external/tblite +++ b/external/tblite @@ -1 +1 @@ -Subproject commit c5844a4d34d30a03118e54a87e3f9fd46b915ee9 +Subproject commit 9075e2e34d5e5b8a2db2d740c3310a0b35c46106 diff --git a/external/xtb b/external/xtb index 26b2801..d7a364b 160000 --- a/external/xtb +++ b/external/xtb @@ -1 +1 @@ -Subproject commit 26b28010e805f7d1aeeef39813feb473e69cc4be +Subproject commit d7a364b353c1f845b22192e0ab505f1c74cea209 diff --git a/src/capabilities/optimiser/OptimiseDipoleScaling.h b/src/capabilities/optimiser/OptimiseDipoleScaling.h index 44e36a2..d5e9954 100644 --- a/src/capabilities/optimiser/OptimiseDipoleScaling.h +++ b/src/capabilities/optimiser/OptimiseDipoleScaling.h @@ -60,15 +60,14 @@ struct OptDipoleFunctor : TFunctor { int operator()(const Vector& scaling, Eigen::VectorXd& fvec) const { // calculation of residuals for (int i = 0; i < m_conformers.size(); ++i) { - auto conf = m_conformers.at(i); - fvec(i) = (conf.getDipole() - conf.CalculateDipoleMoment(scaling, m_bond)).norm(); + const auto& conf = m_conformers.at(i); + fvec(i) = conf.getDipole().norm() - conf.CalculateDipoleMoment(scaling).norm(); } return 0; } int no_parameter; int no_points; std::vector m_conformers; - bool m_bond; int inputs() const { return no_parameter; } int values() const { return no_points; } @@ -80,7 +79,6 @@ inline Vector OptimiseDipoleScaling(const std::vector& conformers, Vec OptDipoleFunctor functor(2, conformers.size()); functor.m_conformers = conformers; - functor.m_bond = bond; Eigen::NumericalDiff numDiff(functor); Eigen::LevenbergMarquardt lm(numDiff); diff --git a/src/main.cpp b/src/main.cpp index 7023bdb..32e2ffd 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -955,50 +955,55 @@ int main(int argc, char **argv) { const Vector one = VectorXd::Ones(mol.AtomCount()); const auto result = OptimiseDipoleScaling(conformers, one); - std::cout << "nonlinear Scaling vector, 1 as initial guess:\n" + std::cout << "nonlinear Scaling vector:\n" << result << "\n" << "Dipole: " << mol.CalculateDipoleMoment(result).norm() << std::endl << std::endl; - const auto result2 = OptimiseDipoleScaling(conformers, theta, true); - std::cout << "nonlinear Scaling vector, with bond true:\n" - << result2 << "\n" - << "Dipole: " << mol.CalculateDipoleMoment(result2, true).norm() << std::endl - << std::endl; - //TODO make it faster... - double r2_anal = 0; - double r2_LM = 0; - double r2_LM2 = 0; - //double r2_anal_dn = 0; - //double r2_LM_dn = 0; - //double r2_LM1_dn = 0; - //double r2_LM2_dn = 0; + double r2_lin = 0; + double r2_nlin = 0; + double r2_lin_diffofnorm = 0; + double r2_nlin_diffofnorm = 0; for (const auto& conf : conformers) { - const double residual = conf.CalculateDipoleMoment(theta).norm() - conf.getDipole().norm(); - r2_anal += residual * residual; - const double residual_1 = conf.CalculateDipoleMoment(result).norm() - conf.getDipole().norm(); - r2_LM += residual_1 * residual_1; - const double residual_2 = conf.CalculateDipoleMoment(result2).norm() - conf.getDipole().norm(); - r2_LM2 += residual_2 * residual_2; - //r2_anal_dn += std::pow((conf.CalculateDipoleMoment(theta) - conf.getDipole()).norm(), 2); - //r2_LM_dn += std::pow((conf.CalculateDipoleMoment(result) - conf.getDipole()).norm(), 2); - //r2_LM1_dn += std::pow((conf.CalculateDipoleMoment(result1) - conf.getDipole()).norm(), 2); - //r2_LM2_dn += std::pow((conf.CalculateDipoleMoment(result2) - conf.getDipole()).norm(), 2); + const double residual = (conf.CalculateDipoleMoment(theta) - conf.getDipole()).norm(); + r2_lin += residual * residual; + const double residual_1 = (conf.CalculateDipoleMoment(result) - conf.getDipole()).norm(); + r2_nlin += residual_1 * residual_1; + const double residual_2 = conf.CalculateDipoleMoment(theta).norm() - conf.getDipole().norm(); + r2_lin_diffofnorm += residual_2 * residual_2; + const double residual_3 = conf.CalculateDipoleMoment(result).norm() - conf.getDipole().norm(); + r2_nlin_diffofnorm += residual_3 * residual_3; } - std::cout << "Residuals:" << std::endl - << "anal: " << r2_anal << std::endl - << "LM 1 as initial guess: " << r2_LM << std::endl - << "LM bond true: " << r2_LM2 << std::endl; + std::cout << "Residuals magnitute of diffrence of Vector:" << std::endl + << "linear: " << r2_lin << std::endl + << "nonlinear " << r2_nlin << std::endl; + std::cout << "Residuals diffrence of magnitute of Vector:" << std::endl + << "linear: " << r2_lin_diffofnorm << std::endl + << "nonlinear: " << r2_lin_diffofnorm << std::endl; std::vector vec_theta(theta.data(), theta.data() + theta.rows() * theta.cols()); + std::vector vec_nonlinear_scaling(result.data(), result.data() + result.rows() * result.cols()); json scaling_vector; scaling_vector["scaling_vector"] = Tools::DoubleVector2String(vec_theta); std::ofstream out(lm_basename + "_scaling_vector.json"); out << scaling_vector << std::endl; + std::ofstream file_dipole; + file_dipole.open(lm_basename + "_dipole.out", std::ios_base::app); + file_dipole << "linear Dipole of Fragments (x y z magn.); nonlinear Dipole of Fragments (x y z magn.); gfn2 Dipoles (x y z magn.)" << std::endl; + for (const auto& conf : conformers){ + const auto dipole_lin = conf.CalculateDipoleMoment(vec_theta); + const auto dipole_nlin = conf.CalculateDipoleMoment(vec_nonlinear_scaling); + const auto dipole_gfn2 = conf.getDipole(); + file_dipole << dipole_lin[0] << " " << dipole_lin[1] << " " << dipole_lin[2] << " " << dipole_lin.norm() << "; "; + file_dipole << dipole_nlin[0] << " " << dipole_nlin[1] << " " << dipole_nlin[2] << " " << dipole_nlin.norm() << "; "; + file_dipole << dipole_gfn2[0] << " " << dipole_gfn2[1] << " " << dipole_gfn2[2] << " " << dipole_gfn2.norm() << std::endl; + }; + file_dipole.close(); + } else if (strcmp(argv[1], "-stride") == 0) {