Skip to content

Add simple molecular dynamics #2527

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
ghutchis merged 19 commits into from
Feb 6, 2026
+764 −278
Merged

Add simple molecular dynamics #2527

Show file tree
Hide file tree
Changes from 17 commits
Commits
Show all changes
19 commits
Select commit Hold shift + click to select a range
340a7e7
Initial implementation of dynamics - no thermostat
ghutchis Oct 8, 2025
faf7887
Add a thermostat - needs proper unit conversions
ghutchis Oct 8, 2025
9c343b0
Some work to do unit conversions
ghutchis Oct 8, 2025
84172c5
Switch to using kJ/mol for internal energies
ghutchis Jan 17, 2026
f9c4502
Try a more complete velocity Verlet implementation
ghutchis Jan 18, 2026
c2ddc07
Make sure Open Babel force fields all output kJ/mol internally
ghutchis Jan 18, 2026
9daf9a0
Finish a velocity Verlet integration method
ghutchis Jan 18, 2026
0273a4b
Add the thermostat
ghutchis Jan 18, 2026
e6063fe
Clamp large gradients - helps with interactive manipulation
ghutchis Jan 18, 2026
4d279e5
Ensure all masses are at least 1.0 amu
ghutchis Jan 18, 2026
3bdf0c6
Clean up mass assignment and ensure velocities and accel obey mask
ghutchis Jan 18, 2026
077392e
Revert "Migrate oopGradient to gradients.h and use them in uff.cpp"
ghutchis Jan 18, 2026
f3d1be1
More work on the thermostat
ghutchis Jan 18, 2026
7bf3629
Merge branch 'master' of github.com:OpenChemistry/avogadrolibs into e…
ghutchis Feb 5, 2026
996d41f
Clamp displacements instead of gradients
ghutchis Feb 6, 2026
0729fc1
Change atom displacement clamping to vary based on mass and temp
ghutchis Feb 6, 2026
7b89a26
Clean up some of the forcefield use of Eigen
ghutchis Feb 6, 2026
193b084
Fix review comments
ghutchis Feb 6, 2026
c60bd6a
Use chi_squared_distribution
ghutchis Feb 6, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
240 changes: 66 additions & 174 deletions avogadro/calc/uff.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -88,6 +88,8 @@ class UFFVdW
Index _atom2;
Real _depth;
Real _x;
Real _x6; // precomputed x^6
Real _x12; // precomputed x^12
};

class UFFPrivate // track the particular calculations for a molecule
Expand Down Expand Up @@ -627,6 +629,9 @@ class UFFPrivate // track the particular calculations for a molecule
sqrt(uffparams[m_atomTypes[i]].D1 * uffparams[m_atomTypes[j]].D1);
v._x =
sqrt(uffparams[m_atomTypes[i]].x1 * uffparams[m_atomTypes[j]].x1);
Real x3 = v._x * v._x * v._x;
v._x6 = x3 * x3;
v._x12 = v._x6 * v._x6;
m_vdws.push_back(v);
}
}
Expand All @@ -638,24 +643,12 @@ class UFFPrivate // track the particular calculations for a molecule
Real energy = 0.0;

for (const UFFBond& bond : m_bonds) {
Index i = bond._atom1;
Index j = bond._atom2;
Real r0 = bond._r0;
Real kb = bond._kb;

Real dx = x[3 * i] - x[3 * j];
Real dy = x[3 * i + 1] - x[3 * j + 1];
Real dz = x[3 * i + 2] - x[3 * j + 2];
Real r = std::hypot(dx, dy, dz);
Real dr = r - r0;

/*
std::cout << " Bond " << i << " " << j << " " << r0 << " " << r << " "
<< dr << std::endl;
*/
Vector3d diff =
x.segment<3>(3 * bond._atom1) - x.segment<3>(3 * bond._atom2);
Real dr = diff.norm() - bond._r0;

// the 0.5 * kb is already in the kb to save a multiplication
energy += kb * dr * dr;
energy += bond._kb * dr * dr;
}
return energy;
}
Expand All @@ -664,22 +657,15 @@ class UFFPrivate // track the particular calculations for a molecule
{
Real energy = 0.0;
for (const UFFAngle& angle : m_angles) {
Index i = angle._atom1;
Index j = angle._atom2;
Index k = angle._atom3;
Real theta0 = angle._theta0 * DEG_TO_RAD;
Real kijk = angle._kijk;

Real dx1 = x[3 * i] - x[3 * j];
Real dy1 = x[3 * i + 1] - x[3 * j + 1];
Real dz1 = x[3 * i + 2] - x[3 * j + 2];
Real dx2 = x[3 * k] - x[3 * j];
Real dy2 = x[3 * k + 1] - x[3 * j + 1];
Real dz2 = x[3 * k + 2] - x[3 * j + 2];
Real r1 = std::hypot(dx1, dy1, dz1);
Real r2 = std::hypot(dx2, dy2, dz2);
Real dot = dx1 * dx2 + dy1 * dy2 + dz1 * dz2;
Real theta = acos(std::clamp(dot / (r1 * r2), -1.0, 1.0));
Vector3d ij = x.segment<3>(3 * angle._atom1) - x.segment<3>(3 * j);
Vector3d kj = x.segment<3>(3 * angle._atom3) - x.segment<3>(3 * j);
Real r1 = ij.norm();
Real r2 = kj.norm();
Real theta = acos(std::clamp(ij.dot(kj) / (r1 * r2), -1.0, 1.0));

/*
std::cout << " Angle " << angle.coordination << " " << i << " " << j
Expand Down Expand Up @@ -747,10 +733,10 @@ class UFFPrivate // track the particular calculations for a molecule
Real c1 = oop._c1;
Real c2 = oop._c2;

Vector3d vi(x[3 * i], x[3 * i + 1], x[3 * i + 2]);
Vector3d vj(x[3 * j], x[3 * j + 1], x[3 * j + 2]);
Vector3d vk(x[3 * k], x[3 * k + 1], x[3 * k + 2]);
Vector3d vl(x[3 * l], x[3 * l + 1], x[3 * l + 2]);
Vector3d vi(x.segment<3>(3 * i));
Vector3d vj(x.segment<3>(3 * j));
Vector3d vk(x.segment<3>(3 * k));
Vector3d vl(x.segment<3>(3 * l));

// use outOfPlaneAngle() from angletools.h
Real angle = outOfPlaneAngle(vi, vj, vk, vl) * DEG_TO_RAD;
Expand All @@ -764,15 +750,10 @@ class UFFPrivate // track the particular calculations for a molecule
{
Real energy = 0.0;
for (const UFFTorsion& torsion : m_torsions) {
Index i = torsion._atom1;
Index j = torsion._atom2;
Index k = torsion._atom3;
Index l = torsion._atom4;

Vector3d vi(x[3 * i], x[3 * i + 1], x[3 * i + 2]);
Vector3d vj(x[3 * j], x[3 * j + 1], x[3 * j + 2]);
Vector3d vk(x[3 * k], x[3 * k + 1], x[3 * k + 2]);
Vector3d vl(x[3 * l], x[3 * l + 1], x[3 * l + 2]);
Vector3d vi(x.segment<3>(3 * torsion._atom1));
Vector3d vj(x.segment<3>(3 * torsion._atom2));
Vector3d vk(x.segment<3>(3 * torsion._atom3));
Vector3d vl(x.segment<3>(3 * torsion._atom4));

Real phi = calculateDihedral(vi, vj, vk, vl) * DEG_TO_RAD;

Expand All @@ -791,27 +772,20 @@ class UFFPrivate // track the particular calculations for a molecule
{
Real energy = 0.0;
for (const UFFVdW& vdw : m_vdws) {
Index i = vdw._atom1;
Index j = vdw._atom2;
Real depth = vdw._depth;
Real xij = vdw._x;
Real x6 = xij * xij * xij * xij * xij * xij;
Real x12 = x6 * x6;

Vector3 atom_i(x[3 * i], x[3 * i + 1], x[3 * i + 2]);
Vector3 atom_j(x[3 * j], x[3 * j + 1], x[3 * j + 2]);
// if the cell is nullptr, we can't do periodic boundary conditions
Vector3 diff = Vector3(x.segment<3>(3 * vdw._atom1)) -
Vector3(x.segment<3>(3 * vdw._atom2));
Real r2;
if (m_cell == nullptr) {
r2 = (atom_i - atom_j).squaredNorm();
r2 = diff.squaredNorm();
} else {
r2 = m_cell->distanceSquared(atom_i, atom_j);
r2 = m_cell->distanceSquared(Vector3(x.segment<3>(3 * vdw._atom1)),
Vector3(x.segment<3>(3 * vdw._atom2)));
}

// we don't need a square root since 6 and 12 are even powers
Real r6 = r2 * r2 * r2;
Real r12 = r6 * r6;
energy += depth * (x12 / r12 - 2 * x6 / r6);
energy += vdw._depth * (vdw._x12 / r12 - 2 * vdw._x6 / r6);
}
return energy;
}
Expand All @@ -821,23 +795,12 @@ class UFFPrivate // track the particular calculations for a molecule
for (const UFFBond& bond : m_bonds) {
Index i = bond._atom1;
Index j = bond._atom2;
Real r0 = bond._r0;
Real kb = bond._kb;

Real dx = x[3 * i] - x[3 * j];
Real dy = x[3 * i + 1] - x[3 * j + 1];
Real dz = x[3 * i + 2] - x[3 * j + 2];

Real r = std::hypot(dx, dy, dz);
Real dr = r - r0;
Real f = 2.0 * kb * dr / r;
grad[3 * i] += f * dx;
grad[3 * i + 1] += f * dy;
grad[3 * i + 2] += f * dz;

grad[3 * j] -= f * dx;
grad[3 * j + 1] -= f * dy;
grad[3 * j + 2] -= f * dz;

Vector3d diff = x.segment<3>(3 * i) - x.segment<3>(3 * j);
Real r = diff.norm();
Vector3d force = 2.0 * bond._kb * (r - bond._r0) / r * diff;
grad.segment<3>(3 * i) += force;
grad.segment<3>(3 * j) -= force;
}
}

Expand All @@ -851,12 +814,8 @@ class UFFPrivate // track the particular calculations for a molecule
Real theta0 = angle._theta0 * DEG_TO_RAD;
Real kijk = angle._kijk;

const Vector3d vi(x[3 * i], x[3 * i + 1], x[3 * i + 2]);
const Vector3d vj(x[3 * j], x[3 * j + 1], x[3 * j + 2]);
const Vector3d vk(x[3 * k], x[3 * k + 1], x[3 * k + 2]);

const Vector3d ij = vi - vj;
const Vector3d kj = vk - vj;
const Vector3d ij = x.segment<3>(3 * i) - x.segment<3>(3 * j);
const Vector3d kj = x.segment<3>(3 * k) - x.segment<3>(3 * j);

Real rij = ij.norm();
Real rkj = kj.norm();
Expand Down Expand Up @@ -962,17 +921,9 @@ class UFFPrivate // track the particular calculations for a molecule
f * (grad_cross_k * dot - crossNorm * grad_dot_k) / denom;

// Add the gradients to the total gradients for each atom
grad[3 * i] += grad_i[0];
grad[3 * i + 1] += grad_i[1];
grad[3 * i + 2] += grad_i[2];

grad[3 * j] += grad_j[0];
grad[3 * j + 1] += grad_j[1];
grad[3 * j + 2] += grad_j[2];

grad[3 * k] += grad_k[0];
grad[3 * k + 1] += grad_k[1];
grad[3 * k + 2] += grad_k[2];
grad.segment<3>(3 * i) += grad_i;
grad.segment<3>(3 * j) += grad_j;
grad.segment<3>(3 * k) += grad_k;
}
}

Expand All @@ -990,10 +941,10 @@ class UFFPrivate // track the particular calculations for a molecule
Real c1 = oop._c1;
Real c2 = oop._c2;

Vector3d vi(x[3 * i], x[3 * i + 1], x[3 * i + 2]);
Vector3d vj(x[3 * j], x[3 * j + 1], x[3 * j + 2]);
Vector3d vk(x[3 * k], x[3 * k + 1], x[3 * k + 2]);
Vector3d vl(x[3 * l], x[3 * l + 1], x[3 * l + 2]);
Vector3d vi(x.segment<3>(3 * i));
Vector3d vj(x.segment<3>(3 * j));
Vector3d vk(x.segment<3>(3 * k));
Vector3d vl(x.segment<3>(3 * l));

// use outOfPlaneAngle() from angletools.h
Real angle = outOfPlaneAngle(vi, vj, vk, vl) * DEG_TO_RAD;
Expand Down Expand Up @@ -1038,49 +989,16 @@ class UFFPrivate // track the particular calculations for a molecule
[[maybe_unused]] Real numerator = cosTheta * sinAngle / sinTheta;

// get the forces on the atoms
Real dj0 =
-dE *
(ik_cross_il[0] - ij[0] + (ik[0] * cosTheta * sinAngle / sinTheta)) /
(rij * sinTheta);
Real dj1 =
-dE *
(ik_cross_il[1] - ij[1] + (ik[1] * cosTheta * sinAngle / sinTheta)) /
(rij * sinTheta);
Real dj2 =
-dE *
(ik_cross_il[2] - ij[2] + (ik[2] * cosTheta * sinAngle / sinTheta)) /
(rij * sinTheta);
grad[3 * j] += dj0;
grad[3 * j + 1] += dj1;
grad[3 * j + 2] += dj2;

Real dk0 =
-dE *
(ij_cross_il[0] - ik[0] + (ij[0] * cosTheta * sinAngle / sinTheta)) /
(rik * sinTheta);
Real dk1 =
-dE *
(ij_cross_il[1] - ik[1] + (ij[1] * cosTheta * sinAngle / sinTheta)) /
(rik * sinTheta);
Real dk2 =
-dE *
(ij_cross_il[2] - ik[2] + (ij[2] * cosTheta * sinAngle / sinTheta)) /
(rik * sinTheta);
grad[3 * k] += dk0;
grad[3 * k + 1] += dk1;
grad[3 * k + 2] += dk2;

Real dl0 = -dE * (-ij_cross_il[0] / sinTheta - il[0] * sinAngle) / ril;
Real dl1 = -dE * (-ij_cross_il[1] / sinTheta - il[1] * sinAngle) / ril;
Real dl2 = -dE * (-ij_cross_il[2] / sinTheta - il[2] * sinAngle) / ril;
grad[3 * l] += dl0;
grad[3 * l + 1] += dl1;
grad[3 * l + 2] += dl2;

Real ratio = cosTheta * sinAngle / sinTheta;
Vector3d dj = -dE / (rij * sinTheta) * (ik_cross_il - ij + ik * ratio);
Vector3d dk = -dE / (rik * sinTheta) * (ij_cross_il - ik + ij * ratio);
Vector3d dl = -dE / ril * (-ij_cross_il / sinTheta - il * sinAngle);

grad.segment<3>(3 * j) += dj;
grad.segment<3>(3 * k) += dk;
grad.segment<3>(3 * l) += dl;
// i is the central atom, so add the other forces
grad[3 * i] -= dj0 + dk0 + dl0;
grad[3 * i + 1] -= dj1 + dk1 + dl1;
grad[3 * i + 2] -= dj2 + dk2 + dl2;
grad.segment<3>(3 * i) -= dj + dk + dl;
}
}

Expand All @@ -1092,10 +1010,10 @@ class UFFPrivate // track the particular calculations for a molecule
Index k = torsion._atom3;
Index l = torsion._atom4;

Vector3d vi(x[3 * i], x[3 * i + 1], x[3 * i + 2]);
Vector3d vj(x[3 * j], x[3 * j + 1], x[3 * j + 2]);
Vector3d vk(x[3 * k], x[3 * k + 1], x[3 * k + 2]);
Vector3d vl(x[3 * l], x[3 * l + 1], x[3 * l + 2]);
Vector3d vi(x.segment<3>(3 * i));
Vector3d vj(x.segment<3>(3 * j));
Vector3d vk(x.segment<3>(3 * k));
Vector3d vl(x.segment<3>(3 * l));

// get the bond vectors
Vector3d ij = vj - vi;
Expand Down Expand Up @@ -1163,21 +1081,10 @@ class UFFPrivate // track the particular calculations for a molecule
Vector3d grad_k = -(grad_i + grad_l + grad_j);

// add the gradients to the total gradients for each atom
grad[3 * i] += dE * grad_i[0];
grad[3 * i + 1] += dE * grad_i[1];
grad[3 * i + 2] += dE * grad_i[2];

grad[3 * j] += dE * grad_j[0];
grad[3 * j + 1] += dE * grad_j[1];
grad[3 * j + 2] += dE * grad_j[2];

grad[3 * k] += dE * grad_k[0];
grad[3 * k + 1] += dE * grad_k[1];
grad[3 * k + 2] += dE * grad_k[2];

grad[3 * l] += dE * grad_l[0];
grad[3 * l + 1] += dE * grad_l[1];
grad[3 * l + 2] += dE * grad_l[2];
grad.segment<3>(3 * i) += dE * grad_i;
grad.segment<3>(3 * j) += dE * grad_j;
grad.segment<3>(3 * k) += dE * grad_k;
grad.segment<3>(3 * l) += dE * grad_l;
}
}

Expand All @@ -1186,39 +1093,25 @@ class UFFPrivate // track the particular calculations for a molecule
for (const UFFVdW& vdw : m_vdws) {
Index i = vdw._atom1;
Index j = vdw._atom2;
Real depth = vdw._depth;
Real xij = vdw._x;

// dE / dr for a Lennard-Jones potential
// E = depth * (x^12 / r^12 - 2 * x^6 / r^6)
// dE / dr = -12 * depth * x^12 / r^13 + 12 * depth * x^6 / r^7
// = -12 * depth * x^6 / r^7 * (x^6 / r^6 - 1)
// dE / dr = -12 * depth * x^6 / r^7 * (x^6 / r^6 - 1)

// TODO: handle unit cells and periodic boundary conditions
Vector3 atom_i(x[3 * i], x[3 * i + 1], x[3 * i + 2]);
Vector3 atom_j(x[3 * j], x[3 * j + 1], x[3 * j + 2]);
Vector3 r = atom_i - atom_j;
Vector3 r = Vector3(x.segment<3>(3 * i)) - Vector3(x.segment<3>(3 * j));
if (m_cell != nullptr) {
r = m_cell->minimumImage(r);
}
Real r2 = r.squaredNorm();

Real dx = r[0];
Real dy = r[1];
Real dz = r[2];

Real r6 = r2 * r2 * r2;
Real r7 = r6 * sqrt(r2);
Real x6 = xij * xij * xij * xij * xij * xij;
Real dE = 12 * depth * x6 / r7 * (1 - x6 / r6);

grad[3 * i] += dE * dx;
grad[3 * i + 1] += dE * dy;
grad[3 * i + 2] += dE * dz;
Real dE = 12 * vdw._depth * vdw._x6 / r7 * (1 - vdw._x6 / r6);

grad[3 * j] -= dE * dx;
grad[3 * j + 1] -= dE * dy;
grad[3 * j + 2] -= dE * dz;
Vector3 force = dE * r;
grad.segment<3>(3 * i) += force;
grad.segment<3>(3 * j) -= force;
}
}
};
Expand Down Expand Up @@ -1278,7 +1171,6 @@ Real UFF::value(const Eigen::VectorXd& x)
energy += constraintEnergies(x);

return energy * KCAL_TO_KJ;
;
}

Real UFF::bondEnergy(const Eigen::VectorXd& x)
Expand Down
Loading
Loading