From a8af13f2edd41a7e8af2743b6ad32612bb4b1f68 Mon Sep 17 00:00:00 2001 From: Gerd Gehrisch Date: Tue, 29 Oct 2024 18:00:12 +0100 Subject: [PATCH] dipole Calculation with linear and nonlinear scaling vector --- external/CxxThreadPool | 2 +- external/tblite | 2 +- external/xtb | 2 +- .../optimiser/OptimiseDipoleScaling.h | 36 +++--- src/capabilities/simplemd.cpp | 103 ++++++++--------- src/capabilities/simplemd.h | 4 +- src/main.cpp | 106 +++++++++--------- 7 files changed, 122 insertions(+), 133 deletions(-) diff --git a/external/CxxThreadPool b/external/CxxThreadPool index e0f2fec..ca5058b 160000 --- a/external/CxxThreadPool +++ b/external/CxxThreadPool @@ -1 +1 @@ -Subproject commit e0f2fec388d0fb1f58f6da1894729a04aa407cbb +Subproject commit ca5058bd27ec88cb7291781068a8635eab74bd32 diff --git a/external/tblite b/external/tblite index 9075e2e..c5844a4 160000 --- a/external/tblite +++ b/external/tblite @@ -1 +1 @@ -Subproject commit 9075e2e34d5e5b8a2db2d740c3310a0b35c46106 +Subproject commit c5844a4d34d30a03118e54a87e3f9fd46b915ee9 diff --git a/external/xtb b/external/xtb index d7a364b..26b2801 160000 --- a/external/xtb +++ b/external/xtb @@ -1 +1 @@ -Subproject commit d7a364b353c1f845b22192e0ab505f1c74cea209 +Subproject commit 26b28010e805f7d1aeeef39813feb473e69cc4be diff --git a/src/capabilities/optimiser/OptimiseDipoleScaling.h b/src/capabilities/optimiser/OptimiseDipoleScaling.h index d5e9954..f1cfc14 100644 --- a/src/capabilities/optimiser/OptimiseDipoleScaling.h +++ b/src/capabilities/optimiser/OptimiseDipoleScaling.h @@ -68,6 +68,7 @@ struct OptDipoleFunctor : TFunctor { 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; } @@ -79,6 +80,7 @@ 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); @@ -108,33 +110,29 @@ inline Vector OptimiseDipoleScaling(const std::vector& conformers, Vec return scaling; } -inline Matrix DipoleScalingCalculation(const std::vector& conformers) +inline Vector DipoleScalingCalculation(const std::vector& conformers) { - std::vector y; //Dipole moments - std::vector F; // Geometry multiplied with partial Charge const auto para_size = conformers[0].AtomCount(); const auto conformer_size = conformers.size(); + Matrix F(3*conformer_size,para_size); // Geometry multiplied with partial Charge + Matrix y(3*conformer_size,1); //Dipoles Matrix FTF = Matrix::Zero(para_size, para_size); Matrix FTy = Matrix::Zero(para_size, 1); for (int i = 0; i < conformer_size; ++i) { - y.push_back(conformers[i].getDipole()); // in eA - F.push_back(conformers[i].ChargeDistribution()); - } - for (auto k : F) { - for (int i = 0; i < para_size; ++i) { - for (int j = 0; j < para_size; ++j) { - FTF(i, j) += k(i, 0) * k(j, 0) + k(i, 1) * k(j, 1) + k(i, 2) * k(j, 2); - } - } - } - - for (int i = 0; i < conformer_size; ++i) { + y(3*i,0) = conformers[i].getDipole()[0]; + y(3*i+1,0) = conformers[i].getDipole()[1]; + y(3*i+2,0) = conformers[i].getDipole()[2]; + const auto& f = conformers[i].ChargeDistribution(); for (int j = 0; j < para_size; ++j) { - FTy(j) += y[i](0) * F[i](j, 0) + y[i](1) * F[i](j, 1) + y[i](2) * F[i](j, 2); + F(3*i,j) = f(j,0); + F(3*i+1,j) = f(j,1); + F(3*i+2,j) = f(j,2); } } - - Matrix Theta = FTF.colPivHouseholderQr().solve(FTy); - + const Vector Theta = (F.transpose()*F).colPivHouseholderQr().solve(F.transpose()*y); + //const Matrix H = (F*(F.transpose()*F).inverse()*F.transpose()).diagonal(); + //std::cout << "diag(H) x y z:" << std::endl; + //for (int i = 0; i < H.rows()/3; ++i) + // std::cout << H(3*i,0) << " " << H(3*i+1,0) << " " << H(3*i+2,0) << " " << std::endl; return Theta; } \ No newline at end of file diff --git a/src/capabilities/simplemd.cpp b/src/capabilities/simplemd.cpp index db40c4a..f6052a3 100644 --- a/src/capabilities/simplemd.cpp +++ b/src/capabilities/simplemd.cpp @@ -339,7 +339,8 @@ bool SimpleMD::Initialise() m_natoms = m_molecule.AtomCount(); m_start_fragments = m_molecule.GetFragments(); - m_scaling_vector = std::vector(m_natoms, 1); + m_scaling_vector_linear = std::vector(m_natoms, 1); + m_scaling_vector_nonlinear = std::vector(m_natoms, 1); if (m_scaling_json != "none") { json scaling; std::ifstream file(m_scaling_json); @@ -350,13 +351,17 @@ bool SimpleMD::Initialise() } catch ([[maybe_unused]] nlohmann::json::parse_error& e) { throw 404; } - std::string scaling_vector; + std::string scaling_vector_linear, scaling_vector_nonlinear; try { - scaling_vector = scaling["scaling_vector"]; + scaling_vector_linear = scaling["scaling_vector_linear"]; + scaling_vector_nonlinear = scaling["scaling_vector_nonlinear"]; } catch ([[maybe_unused]] json::type_error& e) { } - if (!scaling_vector.empty()) { - m_scaling_vector = Tools::String2DoubleVec(scaling_vector, "|"); + if (!scaling_vector_linear.empty()) { + m_scaling_vector_linear = Tools::String2DoubleVec(scaling_vector_linear, "|"); + } + if (!scaling_vector_nonlinear.empty()) { + m_scaling_vector_nonlinear = Tools::String2DoubleVec(scaling_vector_nonlinear, "|"); } } @@ -746,8 +751,6 @@ nlohmann::json SimpleMD::WriteRestartInformation() restart["average_Virial"] = m_average_virial_correction; restart["average_Wall"] = m_average_wall_potential; - restart["scaling_vector"] = Tools::DoubleVector2String(m_scaling_vector); - restart["rattle"] = m_rattle; restart["rattle_maxiter"] = m_rattle_maxiter; restart["rattle_dynamic_tol"] = m_rattle_tolerance; @@ -824,7 +827,7 @@ bool SimpleMD::LoadRestartInformation() bool SimpleMD::LoadRestartInformation(const json& state) { - std::string geometry, scaling_vector, velocities, constrains, xi, Q; + std::string geometry, velocities, constrains, xi, Q; try { m_method = state["method"]; @@ -913,11 +916,6 @@ bool SimpleMD::LoadRestartInformation(const json& state) } catch ([[maybe_unused]] json::type_error& e) { } - try { - scaling_vector = state["scaling_vector"]; - } catch ([[maybe_unused]] json::type_error& e) { - } - try { velocities = state["velocities"]; } catch ([[maybe_unused]] json::type_error& e) { @@ -980,9 +978,6 @@ bool SimpleMD::LoadRestartInformation(const json& state) if (!geometry.empty()) { m_current_geometry = Tools::String2DoubleVec(geometry, "|"); } - if (!scaling_vector.empty()) { - m_scaling_vector = Tools::String2DoubleVec(scaling_vector, "|"); - } if (!velocities.empty()) { m_velocities = Tools::String2DoubleVec(velocities, "|"); } @@ -1109,7 +1104,7 @@ void SimpleMD::start() PrintStatus(); /* Start MD Lopp here */ - for (; m_currentStep < m_maxtime;) { + while (m_currentStep < m_maxtime) { auto step0 = std::chrono::system_clock::now(); if (CheckStop() == true) { @@ -1150,28 +1145,30 @@ void SimpleMD::start() /////////// Dipole if (m_dipole && m_method == "gfn2") { - - //json blob; - // calc partialCharge with gnf2-xtb - //EnergyCalculator interface("gfn2", blob); // set method to gfn2-xtb - //interface.setMolecule(m_molecule); // set molecule for calc - //interface.CalculateEnergy(false, true); // calc energy and charges and dipole moment - //m_molecule.setPartialCharges(interface.Charges()); // calc Partial Charges and give it to mol - //m_molecule.setDipole(interface.Dipole()*au); - - m_curr_dipoles = m_molecule.CalculateDipoleMoments(m_scaling_vector, m_start_fragments); + //linear Dipoles + auto curr_dipoles_lin = m_molecule.CalculateDipoleMoments(m_scaling_vector_linear, m_start_fragments); + std::cout << m_scaling_vector_linear[0] << std::endl; std::ofstream file; - file.open(Basename() + "_dipole.out", std::ios_base::app); - + file.open(Basename() + "_dipole_linear.out", std::ios_base::app); Position d = {0,0,0}; - for (const auto& dipole : m_curr_dipoles) { - d += dipole; - file << dipole.norm() << ", "; + for (const auto& dipole_lin : curr_dipoles_lin) { + d += dipole_lin; + file << dipole_lin[0] << " " << dipole_lin[1] << " " << dipole_lin[2] << " " << dipole_lin.norm() << ", "; } - - file << d[0] << " " << d[1] << " " << d[2] << " " << m_molecule.getDipole()[0] << " " << m_molecule.getDipole()[1] << " " << m_molecule.getDipole()[2] << std::endl; + file << d[0] << " " << d[1] << " " << d[2] << ", " << m_molecule.getDipole()[0] << " " << m_molecule.getDipole()[1] << " " << m_molecule.getDipole()[2] << std::endl; file.close(); - + //nonlinear Dipoles + auto curr_dipoles_nlin = m_molecule.CalculateDipoleMoments(m_scaling_vector_nonlinear, m_start_fragments); + std::cout << m_scaling_vector_nonlinear[0] << std::endl; + std::ofstream file2; + file2.open(Basename() + "_dipole_nonlinear.out", std::ios_base::app); + Position sum = {0,0,0}; + for (const auto& dipole_nlin : curr_dipoles_nlin) { + sum += dipole_nlin; + file2 << dipole_nlin[0] << " " << dipole_nlin[1] << " " << dipole_nlin[2] << " " << dipole_nlin.norm() <<", "; + } + file2 << sum[0] << " " << sum[1] << " " << sum[2] << ", " << m_molecule.getDipole()[0] << " " << m_molecule.getDipole()[1] << " " << m_molecule.getDipole()[2] << std::endl; + file2.close(); } //////////// Dipole @@ -1222,12 +1219,12 @@ void SimpleMD::start() } if (m_writerestart > -1 && m_step % m_writerestart == 0) { - std::ofstream restart_file("curcuma_step_" + std::to_string(int(m_step * m_dT)) + ".json"); + std::ofstream restart_file("curcuma_step_" + std::to_string(static_cast(m_step * m_dT)) + ".json"); json restart; restart[MethodName()[0]] = WriteRestartInformation(); restart_file << restart << std::endl; } - if ((m_step && int(m_step * m_dT) % m_print == 0)) { + if ((m_step && static_cast(m_step * m_dT) % m_print == 0)) { m_Etot = m_Epot + m_Ekin; PrintStatus(); m_time_step = 0; @@ -1257,14 +1254,11 @@ void SimpleMD::start() if (m_thermostat == "csvr") std::cout << "Exchange with heat bath " << m_Ekin_exchange << "Eh" << std::endl; if (m_dipole) { - /* + double dipole = 0.0; - for( auto d : m_collected_dipole) - dipole += d; - dipole /= m_collected_dipole.size(); - std::cout << dipole*2.5418 << " average dipole in Debye and " << dipole*2.5418*3.3356e-30 << " Cm" << std::endl; - */ - std::cout << "Calculated averaged dipole moment " << m_aver_dipol * 2.5418 << " Debye and " << m_aver_dipol * 2.5418 * 3.3356 << " Cm [e-30]" << std::endl; + //std::cout << dipole*2.5418 << " average dipole in Debye and " << dipole*2.5418*3.3356e-30 << " Cm" << std::endl; + + std::cout << "Calculated averaged dipole moment " << m_aver_dipol_linear * 2.5418 << " Debye and " << m_aver_dipol_linear * 2.5418 * 3.3356 << " Cm [e-30]" << std::endl; } #ifdef USE_Plumed @@ -1277,7 +1271,7 @@ void SimpleMD::start() for (int i = 0; i < m_bias_threads.size(); ++i) { auto structures = m_bias_threads[i]->getBiasStructure(); for (int j = 0; j < structures.size(); ++j) { - std::cout << structures[j].rmsd_reference << "\t" << structures[j].energy << "\t" << structures[j].counter / double(m_colvar_incr) * 100 << std::endl; + std::cout << structures[j].rmsd_reference << "\t" << structures[j].energy << "\t" << structures[j].counter / static_cast(m_colvar_incr) * 100 << std::endl; m_rmsd_mtd_molecule.setGeometry(structures[j].geometry); m_rmsd_mtd_molecule.setEnergy(structures[j].energy); @@ -1300,7 +1294,7 @@ void SimpleMD::start() void SimpleMD::AdjustRattleTolerance() { - m_aver_rattle_Temp /= double(m_rattle_counter); + m_aver_rattle_Temp /= static_cast(m_rattle_counter); // std::pair pair(m_rattle_tolerance, m_aver_Temp); @@ -1993,11 +1987,10 @@ void SimpleMD::PrintStatus() const { const auto unix_timestamp = std::chrono::seconds(std::time(NULL)); - int current = std::chrono::milliseconds(unix_timestamp).count(); - double duration = (current - m_unix_started) / (1000.0 * double(m_currentStep)); + const int current = std::chrono::milliseconds(unix_timestamp).count(); + const double duration = (current - m_unix_started) / (1000.0 * static_cast(m_currentStep)); double remaining; - double tmp = (m_maxtime - m_currentStep) * duration / 60; - if (tmp >= 1) + if (const double tmp = (m_maxtime - m_currentStep) * duration / 60; tmp >= 1) remaining = tmp; else remaining = (m_maxtime - m_currentStep) * duration; @@ -2014,7 +2007,7 @@ void SimpleMD::PrintStatus() const #ifdef GCC if (m_dipole) std::cout << fmt::format("{1: ^{0}f} {2: ^{0}f} {3: ^{0}f} {4: ^{0}f} {5: ^{0}f} {6: ^{0}f} {7: ^{0}f} {8: ^{0}f} {9: ^{0}f} {10: ^{0}f} {11: ^{0}f} {12: ^{0}f} {13: ^{0}f} {14: ^{0}f} {15: ^{0}f} {16: ^{0}f}\n", 15, - m_currentStep / 1000, m_Epot, m_aver_Epot, m_Ekin, m_aver_Ekin, m_Etot, m_aver_Etot, m_T, m_aver_Temp, m_wall_potential, m_average_wall_potential, m_aver_dipol * 2.5418 * 3.3356, m_virial_correction, m_average_virial_correction, remaining, m_time_step / 1000.0); + m_currentStep / 1000, m_Epot, m_aver_Epot, m_Ekin, m_aver_Ekin, m_Etot, m_aver_Etot, m_T, m_aver_Temp, m_wall_potential, m_average_wall_potential, m_aver_dipol_linear * 2.5418 * 3.3356, m_virial_correction, m_average_virial_correction, remaining, m_time_step / 1000.0); else std::cout << fmt::format("{1: ^{0}f} {2: ^{0}f} {3: ^{0}f} {4: ^{0}f} {5: ^{0}f} {6: ^{0}f} {7: ^{0}f} {8: ^{0}f} {9: ^{0}f} {10: ^{0}f} {11: ^{0}f} {12: ^{0}f} {13: ^{0}f} {14: ^{0}f} {15: ^{0}f}\n", 15, m_currentStep / 1000, m_Epot, m_aver_Epot, m_Ekin, m_aver_Ekin, m_Etot, m_aver_Etot, m_T, m_aver_Temp, m_wall_potential, m_average_wall_potential, m_virial_correction, m_average_virial_correction, remaining, m_time_step / 1000.0); @@ -2041,10 +2034,10 @@ double SimpleMD::CleanEnergy(double* grad) interface.setMolecule(m_molecule); interface.updateGeometry(m_current_geometry); - double Energy = interface.CalculateEnergy(true); + const double Energy = interface.CalculateEnergy(true); interface.getGradient(grad); if (m_dipole && m_method == "gfn2") { - m_molecule.setDipole(interface.Dipole()*au); + m_molecule.setDipole(interface.Dipole()*au);//in eA m_molecule.setPartialCharges(interface.Charges()); } return Energy; @@ -2054,10 +2047,10 @@ double SimpleMD::FastEnergy(double* grad) { m_interface->updateGeometry(m_current_geometry); - double Energy = m_interface->CalculateEnergy(true); + const double Energy = m_interface->CalculateEnergy(true); m_interface->getGradient(grad); if (m_dipole && m_method == "gfn2") { - m_molecule.setDipole(m_interface->Dipole()*au); + m_molecule.setDipole(m_interface->Dipole()*au);// in eA m_molecule.setPartialCharges(m_interface->Charges()); } return Energy; diff --git a/src/capabilities/simplemd.h b/src/capabilities/simplemd.h index ee28db8..9af4c5b 100644 --- a/src/capabilities/simplemd.h +++ b/src/capabilities/simplemd.h @@ -292,7 +292,7 @@ class SimpleMD : public CurcumaMethod { int m_natoms = 0; int m_dump = 1; - double m_T = 0, m_Epot = 0, m_aver_Epot = 0, m_Ekin = 0, m_aver_Ekin = 0, m_Etot = 0, m_aver_Etot = 0, m_aver_dipol = 0; + double m_T = 0, m_Epot = 0, m_aver_Epot = 0, m_Ekin = 0, m_aver_Ekin = 0, m_Etot = 0, m_aver_Etot = 0, m_aver_dipol_linear = 0; double m_rm_COM = 100; int m_rm_COM_step = -1; int m_hmass = 1; @@ -302,7 +302,7 @@ class SimpleMD : public CurcumaMethod { double m_T0 = 298.15, m_aver_Temp = 0, m_aver_rattle_Temp = 0, m_rmsd = 1.5; double m_x0 = 0, m_y0 = 0, m_z0 = 0; double m_Ekin_exchange = 0.0; - std::vector m_current_geometry, m_mass, m_velocities, m_gradient, m_rmass, m_virial, m_gradient_bias, m_scaling_vector; + std::vector m_current_geometry, m_mass, m_velocities, m_gradient, m_rmass, m_virial, m_gradient_bias, m_scaling_vector_linear, m_scaling_vector_nonlinear; std::vector m_curr_dipoles; std::vector m_atomtype; Molecule m_molecule, m_reference, m_target, m_rmsd_mtd_molecule; diff --git a/src/main.cpp b/src/main.cpp index 32e2ffd..49b194c 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -931,79 +931,77 @@ int main(int argc, char **argv) { std::vector conformers; while (!file.AtEnd()) { // calculation and output dipole moment Molecule mol = file.Next(); // load Molecule - mol.Center(false); - - EnergyCalculator interface("gfn2", blob); // set method to gfn2-xtb and give - interface.setMolecule(mol); // set molecule for calc - interface.CalculateEnergy(false, true); // calc energy and charges and dipole moment - - mol.setPartialCharges(interface.Charges()); // calc Partial Charges and give it to mol - mol.setDipole(interface.Dipole() * au); // in eA - + mol.Center(false); //sets the Centroid to the origin + EnergyCalculator interface("gfn2", blob); // set method to gfn2-xtb + interface.setMolecule(mol); // set molecule + interface.CalculateEnergy(false, true); // calc energy and Wave function + mol.setPartialCharges(interface.Charges()); // calc partial Charges and set it to mol + mol.setDipole(interface.Dipole() * au); //calc dipole moments and set it to mol in eA conformers.push_back(mol); } Molecule mol = conformers.at(0); // maybe molecule in ground state + //Calculation of the scaling vector linear and nonlinear + const auto linear_vector = DipoleScalingCalculation(conformers); //linear + const auto nonlinear_vector = OptimiseDipoleScaling(conformers, linear_vector); //nonlinear + // output scaling vector as JSON + std::vector vec_linear_scaling(linear_vector.data(), linear_vector.data() + linear_vector.rows() * linear_vector.cols()); + std::vector vec_nonlinear_scaling(nonlinear_vector.data(), nonlinear_vector.data() + nonlinear_vector.rows() * nonlinear_vector.cols()); + json scaling_vector; + scaling_vector["scaling_vector_linear"] = Tools::DoubleVector2String(vec_linear_scaling); + scaling_vector["scaling_vector_nonlinear"] = Tools::DoubleVector2String(vec_nonlinear_scaling); + std::ofstream out(lm_basename + "_scaling_vector.json"); + out << scaling_vector << std::endl; - std::cout << "\n xtb2-Dipole: " << mol.getDipole().norm() << std::endl << std::endl; - - const auto theta = DipoleScalingCalculation(conformers); - std::cout << "Analytical Scaling vector:\n" - << theta << "\n" - << "Dipole: " << mol.CalculateDipoleMoment(theta).norm() << std::endl - << std::endl; - - const Vector one = VectorXd::Ones(mol.AtomCount()); - - const auto result = OptimiseDipoleScaling(conformers, one); - std::cout << "nonlinear Scaling vector:\n" - << result << "\n" - << "Dipole: " << mol.CalculateDipoleMoment(result).norm() << std::endl - << std::endl; - - //TODO make it faster... - + double mean_dipole_gfn2 = 0; + double mean_dipole_nonlinear = 0; + double mean_dipole_linear = 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) - 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 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; - + //output Dipole moments + Calculation of Mean Dipole 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; + file_dipole << "linear Dipole (x y z magn.); nonlinear Dipole (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_lin = conf.CalculateDipoleMoment(vec_linear_scaling); const auto dipole_nlin = conf.CalculateDipoleMoment(vec_nonlinear_scaling); const auto dipole_gfn2 = conf.getDipole(); + mean_dipole_linear += dipole_lin.norm()/ conformers.size(); + mean_dipole_nonlinear += dipole_nlin.norm()/ conformers.size(); + mean_dipole_gfn2 += dipole_gfn2.norm()/ conformers.size(); 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; + const double residual = (conf.CalculateDipoleMoment(linear_vector) - conf.getDipole()).norm(); + r2_lin += residual * residual; + const double residual_1 = (conf.CalculateDipoleMoment(nonlinear_vector) - conf.getDipole()).norm(); + r2_nlin += residual_1 * residual_1; + const double residual_2 = conf.CalculateDipoleMoment(linear_vector).norm() - conf.getDipole().norm(); + r2_lin_diffofnorm += residual_2 * residual_2; + const double residual_3 = conf.CalculateDipoleMoment(nonlinear_vector).norm() - conf.getDipole().norm(); + r2_nlin_diffofnorm += residual_3 * residual_3; }; file_dipole.close(); + std::cout << "\nMean xtb2-Dipole: " << mean_dipole_gfn2 << " [eA], " << mean_dipole_gfn2/0.2082 << " [D]" << std::endl; + std::cout << "Mean linear Dipole: " << mean_dipole_linear << " [eA], " << mean_dipole_linear/0.2082 << " [D]" << std::endl + << "Mean nonlinear Dipole: " << mean_dipole_nonlinear << " [eA], " << mean_dipole_nonlinear/0.2082 << " [D]" << std::endl + << std::endl; + + std::cout << "linear Scaling vector:\n" + << linear_vector << "\n" + << "nonlinear Scaling vector:\n" + << nonlinear_vector << "\n" + << std::endl; + + std::cout << "Square Sum of Residuals of Components:" << std::endl + << "linear: " << r2_lin << std::endl + << "nonlinear " << r2_nlin << std::endl; + std::cout << "Square Sum of Residuals of Magnitudes" << std::endl + << "linear: " << r2_lin_diffofnorm << std::endl + << "nonlinear: " << r2_nlin_diffofnorm << std::endl; + } else if (strcmp(argv[1], "-stride") == 0) {