Skip to content

Commit

Permalink
stop md if there are too many errors
Browse files Browse the repository at this point in the history 
Signed-off-by: Conrad Hübler <[email protected]>
  • Loading branch information
@conradhuebler
conradhuebler committed Jan 30, 2024
1 parent d109900 commit 3309893
Show file tree
Hide file tree
Showing 6 changed files with 37 additions and 10 deletions.
40 changes: 31 additions & 9 deletions src/capabilities/simplemd.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -639,7 +639,7 @@ void SimpleMD::start()
WallPotential(gradient);
Integrator(coord, gradient);

if (m_unstable) {
if (m_unstable || m_interface->Error() || m_interface->HasNan()) {
PrintStatus();
fmt::print(fg(fmt::color::salmon) | fmt::emphasis::bold, "Simulation got unstable, exiting!\n");

Expand Down Expand Up @@ -689,9 +689,10 @@ void SimpleMD::start()
*/
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::ofstream restart_file("curcuma_final.json");
restart_file << WriteRestartInformation() << std::endl;

std::remove("curcuma_restart.json");
delete[] coord;
delete[] gradient;
}
Expand Down Expand Up @@ -747,9 +748,11 @@ void SimpleMD::Rattle(double* coord, double* grad)
* like dT^3 -> dT^2 and
* updated velocities of the second atom (minus instread of plus)
*/

// TriggerWriteRestart();
double m_dT_inverse = 1 / m_dT;
// double T_begin = 0, T_rattle_1 = 0, T_verlet = 0;
std::vector<int> moved(m_natoms, 0);
// double kin = 0;
for (int i = 0; i < m_natoms; ++i) {
coord[3 * i + 0] = m_current_geometry[3 * i + 0] + m_dT * m_velocities[3 * i + 0] - 0.5 * grad[3 * i + 0] * m_rmass[3 * i + 0] * m_dt2;
coord[3 * i + 1] = m_current_geometry[3 * i + 1] + m_dT * m_velocities[3 * i + 1] - 0.5 * grad[3 * i + 1] * m_rmass[3 * i + 1] * m_dt2;
Expand All @@ -758,9 +761,12 @@ void SimpleMD::Rattle(double* coord, double* grad)
m_velocities[3 * i + 0] -= 0.5 * m_dT * grad[3 * i + 0] * m_rmass[3 * i + 0];
m_velocities[3 * i + 1] -= 0.5 * m_dT * grad[3 * i + 1] * m_rmass[3 * i + 1];
m_velocities[3 * i + 2] -= 0.5 * m_dT * grad[3 * i + 2] * m_rmass[3 * i + 2];
// kin += m_mass[i] * (m_velocities[3 * i] * m_velocities[3 * i] + m_velocities[3 * i + 1] * m_velocities[3 * i + 1] + m_velocities[3 * i + 2] * m_velocities[3 * i + 2]);
}

// T_begin = kin / (kb_Eh * m_dof);
double iter = 0;
double max_mu = 10;

while (iter < m_rattle_maxiter) {
iter++;
for (auto bond : m_bond_constrained) {
Expand All @@ -785,6 +791,8 @@ void SimpleMD::Rattle(double* coord, double* grad)
moved[j] = 1;

double lambda = r / (1 * (m_rmass[i] + m_rmass[j]) * scalarproduct);
while (std::abs(lambda) > max_mu)
lambda /= 2;
coord[3 * i + 0] += dx * lambda * 0.5 * m_rmass[i];
coord[3 * i + 1] += dy * lambda * 0.5 * m_rmass[i];
coord[3 * i + 2] += dz * lambda * 0.5 * m_rmass[i];
Expand All @@ -800,17 +808,23 @@ void SimpleMD::Rattle(double* coord, double* grad)
m_velocities[3 * j + 0] -= dx * lambda * 0.5 * m_rmass[j] * m_dT_inverse;
m_velocities[3 * j + 1] -= dy * lambda * 0.5 * m_rmass[j] * m_dT_inverse;
m_velocities[3 * j + 2] -= dz * lambda * 0.5 * m_rmass[j] * m_dT_inverse;
// std::cout << i << " " << j << " " <<lambda << " :: " ;
}
}
}
}

// std::cout << std::endl;
m_Epot = Energy(coord, grad);
double ekin = 0.0;
// kin = 0;

for (int i = 0; i < m_natoms; ++i) {
// kin += m_mass[i] * (m_velocities[3 * i] * m_velocities[3 * i] + m_velocities[3 * i + 1] * m_velocities[3 * i + 1] + m_velocities[3 * i + 2] * m_velocities[3 * i + 2]);

m_velocities[3 * i + 0] -= 0.5 * m_dT * grad[3 * i + 0] * m_rmass[3 * i + 0];
m_velocities[3 * i + 1] -= 0.5 * m_dT * grad[3 * i + 1] * m_rmass[3 * i + 1];
m_velocities[3 * i + 2] -= 0.5 * m_dT * grad[3 * i + 2] * m_rmass[3 * i + 2];
// ekin += m_mass[i] * (m_velocities[3 * i] * m_velocities[3 * i] + m_velocities[3 * i + 1] * m_velocities[3 * i + 1] + m_velocities[3 * i + 2] * m_velocities[3 * i + 2]);

m_current_geometry[3 * i + 0] = coord[3 * i + 0];
m_current_geometry[3 * i + 1] = coord[3 * i + 1];
Expand All @@ -822,7 +836,9 @@ void SimpleMD::Rattle(double* coord, double* grad)
}
m_virial_correction = 0;
iter = 0;

// T_rattle_1 = kin / (kb_Eh * m_dof);
// T_verlet = ekin / (kb_Eh * m_dof);
ekin = 0.0;
while (iter < m_rattle_maxiter) {
iter++;
for (auto bond : m_bond_constrained) {
Expand All @@ -840,7 +856,9 @@ void SimpleMD::Rattle(double* coord, double* grad)
double r = (dx) * (dvx) + (dy) * (dvy) + (dz) * (dvz);

double mu = -1 * r / ((m_rmass[i] + m_rmass[j]) * distance);
if (std::abs(mu) > m_rattle_tolerance) {
while (std::abs(mu) > max_mu)
mu /= 2;
if (std::abs(mu) > m_rattle_tolerance && std::abs(mu) < max_mu) {
m_virial_correction += mu * distance;
m_velocities[3 * i + 0] += dx * mu * m_rmass[i];
m_velocities[3 * i + 1] += dy * mu * m_rmass[i];
Expand All @@ -849,6 +867,7 @@ void SimpleMD::Rattle(double* coord, double* grad)
m_velocities[3 * j + 0] -= dx * mu * m_rmass[j];
m_velocities[3 * j + 1] -= dy * mu * m_rmass[j];
m_velocities[3 * j + 2] -= dz * mu * m_rmass[j];
// std::cout << i << " " << j << " " << mu << " :: " ;
}
}
}
Expand All @@ -859,7 +878,10 @@ void SimpleMD::Rattle(double* coord, double* grad)
}
ekin *= 0.5;
double T = 2.0 * ekin / (kb_Eh * m_dof);
m_unstable = T > 1000 * m_T;
m_unstable = T > 10000 * m_T;
// std::cout << std::endl;

// std::cout << T_begin << " " << T_rattle_1 << " " << T_verlet << " " << T << std::endl;
m_T = T;
}

Expand Down Expand Up @@ -1267,6 +1289,7 @@ bool SimpleMD::WriteGeometry()
geometry(i, 1) = m_current_geometry[3 * i + 1];
geometry(i, 2) = m_current_geometry[3 * i + 2];
}
TriggerWriteRestart();
// int f1 = m_molecule.GetFragments().size();
m_molecule.setGeometry(geometry);
auto m = m_molecule.DistanceMatrix();
Expand Down Expand Up @@ -1320,7 +1343,6 @@ void SimpleMD::CSVR()
double alpha2 = c + (1 - c) * (SNf + R * R) * Ekin_target / (m_dof * m_Ekin) + 2 * R * sqrt(c * (1 - c) * Ekin_target / (m_dof * m_Ekin));
m_Ekin_exchange += m_Ekin * (alpha2 - 1);
double alpha = sqrt(alpha2);

for (int i = 0; i < 3 * m_natoms; ++i) {
m_velocities[i] *= alpha;
}
Expand Down
1 change: 0 additions & 1 deletion src/capabilities/simplemd.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,7 +60,6 @@ static json CurcumaMDJson{
{ "impuls_scaling", 0.75 },
{ "writeinit", false },
{ "initfile", "none" },
{ "constrain_file", "none" },
{ "norestart", false },
{ "writerestart", 1000 },
{ "rattle", false },
Expand Down
1 change: 1 addition & 0 deletions src/core/energycalculator.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,7 @@ EnergyCalculator::EnergyCalculator(const std::string& method, const json& contro
m_tblite = new TBLiteInterface(controller);
m_ecengine = [this](bool gradient, bool verbose) {
this->CalculateTBlite(gradient, verbose);
m_error = this->m_tblite->Error();
};
m_charges = [this]() {
return this->m_tblite->Charges();
Expand Down
2 changes: 2 additions & 0 deletions src/core/energycalculator.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,6 +64,7 @@ class EnergyCalculator {

bool HasNan() const { return m_containsNaN; }

bool Error() const { return m_error; }
#ifdef USE_TBLITE
TBLiteInterface* getTBLiterInterface() const
{
Expand Down Expand Up @@ -151,4 +152,5 @@ class EnergyCalculator {

bool m_initialised = false;
bool m_containsNaN = false;
bool m_error = false;
};
1 change: 1 addition & 0 deletions src/core/tbliteinterface.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -260,6 +260,7 @@ double TBLiteInterface::GFNCalculation(int parameter, double* grad)
tblite_get_singlepoint(m_ctx, m_tblite_mol, m_tblite_calc, m_tblite_res);
tblite_get_result_energy(m_error, m_tblite_res, &energy);
if (tblite_check_context(m_ctx)) {
m_error_count++;
std::cerr << "Error during SCF Calculation, reset wavefunction ..." << std::endl;
tbliteContextError();
tbliteError();
Expand Down
2 changes: 2 additions & 0 deletions src/core/tbliteinterface.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,6 +54,7 @@ class TBLiteInterface {
bool UpdateMolecule(const Molecule& molecule);
bool UpdateMolecule(const double* coord);

bool Error() { return m_error_count >= 10; }
double GFNCalculation(int parameter = 2, double* grad = 0);

void clear();
Expand All @@ -78,6 +79,7 @@ class TBLiteInterface {
int m_maxiter = 100;
int m_verbose = 0;
int m_guess = 0;
int m_error_count = 0;
double m_damping = 0.5;
double m_temp = 1000;
double m_cpcm_eps = -1, m_alpb_eps = -1;
Expand Down

0 comments on commit 3309893

Please sign in to comment.