Skip to content

Commit

Permalink
some fixes and changes
Browse files Browse the repository at this point in the history 
Signed-off-by: Conrad Hübler <Conrad.Huebler@gmx.net>
  • Loading branch information
@conradhuebler
conradhuebler committed Mar 22, 2024
1 parent 827d2f5 commit e91226a
Show file tree
Hide file tree
Showing 16 changed files with 266 additions and 197 deletions.
2 changes: 1 addition & 1 deletion CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -160,7 +160,7 @@ set(curcuma_core_SRC
src/capabilities/qmdfffit.cpp
src/core/energycalculator.cpp
src/core/molecule.cpp
#src/core/pseudoff.cpp
src/core/fileiterator.cpp
src/core/eigen_uff.cpp
src/core/qmdff.cpp
src/core/eht.cpp
Expand Down
1 change: 0 additions & 1 deletion src/capabilities/confscan.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -268,7 +268,6 @@ void ConfScan::LoadControlJson()
bool ConfScan::openFile()
{
bool xyzfile = std::string(m_filename).find(".xyz") != std::string::npos || std::string(m_filename).find(".trj") != std::string::npos;

if (xyzfile == false)
throw 1;

Expand Down
17 changes: 9 additions & 8 deletions src/capabilities/hessian.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,6 @@ HessianThread::HessianThread(const json& controller, int i, int j, int xi, int x
{
setAutoDelete(true);
m_method = m_controller["method"];

if (m_fullnumerical)
m_schema = [this, i, j, xi, xj]() {
this->Numerical();
Expand Down Expand Up @@ -120,7 +119,7 @@ void HessianThread::Seminumerical()
energy.updateGeometry(m_geom_im_jp);
energy.CalculateEnergy(true, false);
Matrix gradientm = energy.Gradient();
m_gradient = (gradientp - gradientm) / (2 * m_d) / au / au;
m_gradient = (gradientp - gradientm) / (2 * m_d);
}

Hessian::Hessian(const std::string& method, const json& controller, bool silent)
Expand All @@ -132,13 +131,13 @@ Hessian::Hessian(const std::string& method, const json& controller, bool silent)
UpdateController(controller);
m_threads = m_defaults["threads"];
/* Yeah, thats not really correct, but it works a bit */
if (m_method.compare("gfnff") == 0) {
if (m_method.compare("uff") == 0) {
m_scale_functions = [](double val) -> double {
return val * 2720.57 - 0.0338928;
return val * 5150.4 + 47.349;
};
} else {
m_scale_functions = [](double val) -> double {
return val * 2720.57 - 0.0338928;
return val * 5150.4 + 47.349;
};
}
}
Expand All @@ -148,14 +147,16 @@ Hessian::Hessian(const json& controller, bool silent)
, m_controller(controller)
{
UpdateController(controller);
m_threads = m_defaults["threads"];

/* Yeah, thats not really correct, but it works a bit */
if (m_method.compare("gfnff") == 0) {
if (m_method.compare("uff") == 0) {
m_scale_functions = [](double val) -> double {
return val * 2720.57 - 0.0338928;
return val * 5150.4 + 47.349;
};
} else {
m_scale_functions = [](double val) -> double {
return val * 2720.57 - 0.0338928;
return val * 5150.4 + 47.349;
};
}
}
Expand Down
14 changes: 12 additions & 2 deletions src/capabilities/simplemd.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -887,6 +887,18 @@ void SimpleMD::Rattle(double* coord, double* grad)
m_T = T;
}

void SimpleMD::Rattle_Verlet_First(double* coord, double* grad)
{
}

void SimpleMD::Rattle_Constrain_First(double* coord, double* grad)
{
}

void SimpleMD::Rattle_Verlet_Second(double* coord, double* grad)
{
}

double SimpleMD::ApplySphericLogFermiWalls(double* grad)
{
double potential = 0;
Expand Down Expand Up @@ -1092,14 +1104,12 @@ void SimpleMD::RemoveRotations(std::vector<double>& velo)

Position rlm = { 0, 0, 0 }, ram = { 0, 0, 0 };
for (int i : fragments[f]) {

rlm(0) = rlm(0) + m_mass[i] * velo[3 * i + 0];
rlm(1) = rlm(1) + m_mass[i] * velo[3 * i + 1];
rlm(2) = rlm(2) + m_mass[i] * velo[3 * i + 2];
}

for (int i : fragments[f]) {

ram(0) = (omega(1) * geom(i, 2) - omega(2) * geom(i, 1));
ram(1) = (omega(2) * geom(i, 0) - omega(0) * geom(i, 2));
ram(2) = (omega(0) * geom(i, 1) - omega(1) * geom(i, 0));
Expand Down
5 changes: 5 additions & 0 deletions src/capabilities/simplemd.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -143,6 +143,11 @@ class SimpleMD : public CurcumaMethod {
void Verlet(double* coord, double* grad);
void Rattle(double* coord, double* grad);

void Rattle_Verlet_First(double* coord, double* grad);
void Rattle_Constrain_First(double* coord, double* grad);
void Rattle_Verlet_Second(double* coord, double* grad);
void Rattle_Constrain_Second(double* coord, double* grad);

void RemoveRotation(std::vector<double>& velo);
void RemoveRotations(std::vector<double>& velo);

Expand Down
6 changes: 3 additions & 3 deletions src/core/eigen_uff.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1829,9 +1829,9 @@ double eigenUFF::Calculate(bool grd, bool verbose)
<< "HHRepCorrection " << m_final_factor * m_hh_scaling * energy_hh << " Eh" << std::endl
<< std::endl;

for (int i = 0; i < m_atom_types.size(); ++i) {
std::cout << m_gradient(i, 0) << " " << m_gradient(i, 1) << " " << m_gradient(i, 2) << std::endl;
}
// for (int i = 0; i < m_atom_types.size(); ++i) {
// std::cout << m_gradient(i, 0) << " " << m_gradient(i, 1) << " " << m_gradient(i, 2) << std::endl;
// }
}
return energy;
}
74 changes: 36 additions & 38 deletions src/core/energycalculator.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,11 @@
#include "src/core/xtbinterface.h"
#endif

#ifndef _WIN32
#include <filesystem>
namespace fs = std::filesystem;
#endif

#include <functional>

#include "forcefieldgenerator.h"
Expand All @@ -33,8 +38,16 @@

EnergyCalculator::EnergyCalculator(const std::string& method, const json& controller)
: m_method(method)
, m_controller(controller)

{
if (controller.contains("param_file")) {
m_param_file = controller["param_file"];
}

if (controller.contains("write_param")) {
m_writeparam = controller["write_param"];
}

m_charges = []() {
return std::vector<double>{};
};
Expand Down Expand Up @@ -120,17 +133,6 @@ EnergyCalculator::EnergyCalculator(const std::string& method, const json& contro

} else if (std::find(m_ff_methods.begin(), m_ff_methods.end(), m_method) != m_ff_methods.end()) { // Just D4 energy calculator requested
m_forcefield = new ForceField(controller);
/*
std::string file = controller["param"];
nlohmann::json parameter;
std::ifstream parameterfile(file);
try {
parameterfile >> parameter;
} catch (nlohmann::json::type_error& e) {
} catch (nlohmann::json::parse_error& e) {
}
ForceFieldGenerator ff;
m_forcefield->setParameter(parameter);*/
m_ecengine = [this](bool gradient, bool verbose) {
this->CalculateFF(gradient, verbose);
};
Expand Down Expand Up @@ -172,25 +174,11 @@ void EnergyCalculator::setMolecule(const Molecule& molecule)
{
m_atoms = molecule.AtomCount();

// m_atom_type[m_atoms];
std::vector<int> atoms = molecule.Atoms();
m_coord = new double[3 * m_atoms];
m_grad = new double[3 * m_atoms];
// std::vector<std::array<double, 3>> geom(m_atoms);
m_gradient = Eigen::MatrixXd::Zero(m_atoms, 3);
m_geometry = Eigen::MatrixXd::Zero(m_atoms, 3);
for (int i = 0; i < m_atoms; ++i) {
std::pair<int, Position> atom = molecule.Atom(i);
m_geometry(i, 0) = atom.second(0);
m_geometry(i, 1) = atom.second(1);
m_geometry(i, 2) = atom.second(2);
/*
geom[i][0] = atom.second(0);
geom[i][1] = atom.second(1);
geom[i][2] = atom.second(2);*/
// m_gradient.push_back({ 0, 0, 0 });
}
// m_geometry = geom;
m_geometry = molecule.getGeometry();
if (std::find(m_uff_methods.begin(), m_uff_methods.end(), m_method) != m_uff_methods.end()) { // UFF energy calculator requested
m_uff->setMolecule(atoms, m_geometry);
m_uff->Initialise(molecule.Bonds());
Expand Down Expand Up @@ -228,18 +216,28 @@ void EnergyCalculator::setMolecule(const Molecule& molecule)
m_qmdff->Initialise();

} else if (std::find(m_ff_methods.begin(), m_ff_methods.end(), m_method) != m_ff_methods.end()) { //

ForceFieldGenerator ff(m_controller);
ff.setMolecule(molecule);
ff.Generate();
json parameters = ff.getParameter();
// std::ofstream parameterfile("uff.json");
// parameterfile << parameters;
if (m_parameter.size() == 0) {
if (!std::filesystem::exists(m_param_file)) {
ForceFieldGenerator ff(m_controller);
ff.setMolecule(molecule);
ff.Generate();
m_parameter = ff.getParameter();
if (m_writeparam) {
std::ofstream parameterfile("ff_param.json");
parameterfile << m_parameter;
}
} else {
std::ifstream parameterfile(m_param_file);
try {
parameterfile >> m_parameter;
} catch (nlohmann::json::type_error& e) {
} catch (nlohmann::json::parse_error& e) {
}
}
}
m_forcefield->setAtomTypes(molecule.Atoms());
m_forcefield->setParameter(parameters);
// m_forcefield->setAtomCount(molecule.AtomCount());
// m_forcefield->setMolecule(atoms, geom);
// m_forcefield->Initialise();

m_forcefield->setParameter(m_parameter);

} else { // Fall back to UFF?
}
Expand Down
9 changes: 6 additions & 3 deletions src/core/energycalculator.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,10 @@
#include "json.hpp"
using json = nlohmann::json;

static json EnergyCalculatorJson{
{ "param_file", "none" }
};

class EnergyCalculator {
public:
EnergyCalculator(const std::string& method, const json& controller);
Expand All @@ -56,7 +60,6 @@ class EnergyCalculator {
void updateGeometry(const double* coord);
void updateGeometry(const std::vector<double>& geometry);

// void updateGeometry(const std::vector<std::array<double, 3>>& geometry);
void updateGeometry(const Matrix& geometry);
void updateGeometry(const Eigen::VectorXd& geometry);

Expand Down Expand Up @@ -168,9 +171,8 @@ class EnergyCalculator {
std::function<std::vector<double>()> m_charges, m_dipole;
std::function<std::vector<std::vector<double>>()> m_bonds;
json m_parameter;
std::string m_method;
std::string m_method, m_param_file;
Matrix m_geometry, m_gradient;
// Matrix m_eigen_geometry, m_eigen_gradient;
double m_energy;
double *m_coord, *m_grad;

Expand All @@ -181,4 +183,5 @@ class EnergyCalculator {
bool m_initialised = false;
bool m_containsNaN = false;
bool m_error = false;
bool m_writeparam = false;
};
Loading

0 comments on commit e91226a

Please sign in to comment.