Add simple molecular dynamics#2527
Conversation
Signed-off-by: Geoff Hutchison <geoff.hutchison@gmail.com>
Signed-off-by: Geoff Hutchison <geoff.hutchison@gmail.com>
Signed-off-by: Geoff Hutchison <geoff.hutchison@gmail.com>
Signed-off-by: Geoff Hutchison <geoff.hutchison@gmail.com>
Signed-off-by: Geoff Hutchison <geoff.hutchison@gmail.com>
MMFF94 and UFF use kcal, GAFF uses kJ/mol Signed-off-by: Geoff Hutchison <geoff.hutchison@gmail.com>
Signed-off-by: Geoff Hutchison <geoff.hutchison@gmail.com>
|
This seems to mostly work now, although dragging atoms or otherwise generating large gradients leads to instability. |
|
@coderabbitai review |
✅ Actions performedReview triggered.
|
WalkthroughThis PR introduces molecular dynamics simulation support to the AutoOpt plugin via a new CSVR thermostat and velocity-Verlet integrator, while refactoring multiple force field components to use vectorized Eigen operations for coordinate handling and precomputed van der Waals terms. Changes
Sequence Diagram(s)sequenceDiagram
participant UI as User Interface
participant AutoOpt as AutoOpt Plugin
participant ForceField as ForceField Engine
participant Thermostat as CSVR Thermostat
participant Molecule as Molecule State
rect rgba(100, 150, 255, 0.5)
Note over UI,Molecule: Dynamics Task Initialization
UI->>AutoOpt: select task=Dynamics
AutoOpt->>AutoOpt: taskChanged(): enable dynamics UI
UI->>AutoOpt: set temperature & timestep
AutoOpt->>AutoOpt: setMasses from atoms
AutoOpt->>Thermostat: initializeVelocities(mass)
AutoOpt->>Molecule: begin merge mode
end
rect rgba(150, 200, 150, 0.5)
Note over AutoOpt,Thermostat: Velocity-Verlet Dynamics Loop
loop Each dynamics step
AutoOpt->>Molecule: update positions (v·Δt + ½a·Δt²)
AutoOpt->>ForceField: compute forces from gradient
AutoOpt->>ForceField: gradient(positions)
ForceField-->>AutoOpt: return gradient
AutoOpt->>AutoOpt: compute acceleration (F/m)
AutoOpt->>AutoOpt: update velocities (v + ½(a_old + a_new)·Δt)
AutoOpt->>Thermostat: apply(velocities, masses)
Thermostat->>Thermostat: rescale velocities (CSVR)
Thermostat-->>AutoOpt: rescaled velocities
AutoOpt->>Molecule: update atom positions
UI->>AutoOpt: draw overlay (show T instead of ΔE)
end
end
rect rgba(200, 150, 100, 0.5)
Note over AutoOpt,Molecule: Task Completion
UI->>AutoOpt: stop simulation
AutoOpt->>Molecule: end merge mode
AutoOpt->>AutoOpt: cleanup thermostat
end
Estimated code review effort🎯 4 (Complex) | ⏱️ ~60 minutes Possibly related PRs
Suggested reviewers
🚥 Pre-merge checks | ✅ 2 | ❌ 1❌ Failed checks (1 warning)
✅ Passed checks (2 passed)
✏️ Tip: You can configure your own custom pre-merge checks in the settings. ✨ Finishing touches
🧪 Generate unit tests (beta)
Thanks for using CodeRabbit! It's free for OSS, and your support helps us grow. If you like it, consider giving us a shout-out. Comment |
![coderabbitai[bot]](https://avatars.githubusercontent.com/in/347564?s=60&v=4){kind=small}
There was a problem hiding this comment.
Actionable comments posted: 2
🤖 Fix all issues with AI agents
In `@avogadro/qtplugins/autoopt/autoopt.cpp`:
- Around line 462-470: The code populates m_masses from m_molecule but doesn't
guard against zero masses, leading to infinities when later dividing via
gradient.array() / m_masses.array(); change the population of m_masses in the
block that sets it (the loop using m_molecule->atom(i).mass()) to replace any
returned mass <= 0 with a small positive epsilon (e.g., 1e-8) or a configurable
minimum mass, so subsequent operations that use m_masses (including the
gradient.array() / m_masses.array() division) never divide by zero; ensure you
update m_masses only in that section and document the chosen epsilon where
m_masses is set.
- Around line 472-478: The code fails to compile because units::FORCE_CONVERSION
is undefined; add a units namespace (or extend the existing units header) and
define a constexpr double FORCE_CONVERSION that converts gradient in kJ/mol/Å to
acceleration in Å/fs² (apply kJ→J, per-mol→per-molecule via Avogadro, Å→m,
amu→kg, and fs→s² conversions) and use the computed numeric constant (~0.01);
place the constant in a header included by autoopt.cpp (so m_method->gradient,
m_acceleration and m_masses usages can reference units::FORCE_CONVERSION).
Ensure the constant is documented briefly and has internal linkage/constexpr
storage to avoid ODR issues.
🧹 Nitpick comments (4)
avogadro/core/atom.h (1)
293-298: Consider implementing isotope handling for consistency withMolecule::mass().The
Molecule::mass()method inmolecule.cpp(lines 969-979) handles isotopes by checking if an isotope is set and usingElements::isotopeMass(). For molecular dynamics accuracy, the atom-level mass should follow the same pattern sinceisotope()is already available on this class.♻️ Suggested implementation matching Molecule::mass() pattern
template <class Molecule_T> Real AtomTemplate<Molecule_T>::mass() const { - // TODO: handle isotopes - return Elements::mass(atomicNumber()); + unsigned short iso = isotope(); + if (iso > 0) + return Elements::isotopeMass(atomicNumber(), iso); + return Elements::mass(atomicNumber()); }Would you like me to open an issue to track this if you prefer to defer it?
avogadro/qtplugins/autoopt/autoopt.cpp (2)
214-236: Missing null checks on findChild results.If
taskChangedis ever called before the tool widget is fully initialized, thefindChildcalls could return nullptr, causing a crash onsetEnabled(). Consider adding defensive null checks.♻️ Suggested defensive null checks
void AutoOpt::taskChanged(int index) { m_task = index; auto temperatureSpinBox = m_toolWidget->findChild<QDoubleSpinBox*>("temperatureSpinBox"); auto timeStepSpinBox = m_toolWidget->findChild<QDoubleSpinBox*>("timeStepSpinBox"); + if (temperatureSpinBox == nullptr || timeStepSpinBox == nullptr) + return; + bool enabled = (index == 1); temperatureSpinBox->setEnabled(enabled); timeStepSpinBox->setEnabled(enabled);
297-304: Degrees of freedom calculation may need adjustment for nonlinear molecules.The DOF is set to
3n - 3, which correctly removes translational degrees of freedom. However, for nonlinear molecules, rotational DOF should also be subtracted (3n - 6). For linear molecules it would be3n - 5. This affects temperature calculation accuracy.This is a known simplification for interactive dynamics and may be acceptable for the current use case, but could affect temperature stability for small molecules.
avogadro/qtplugins/forcefield/obenergy.cpp (1)
246-253: Inconsistent ordering compared to obmmenergy.cpp.In
obenergy.cpp, the unit conversion happens beforecleanGradients(), while inobmmenergy.cppit happens after. This may cause subtle differences in how frozen atom gradients are handled. Consider aligning the ordering for consistency.In
obmmenergy.cpp:cleanGradients(grad); if (m_identifier != "GAFF") grad *= Calc::KCAL_TO_KJ; constraintGradients(x, grad);In
obenergy.cpp:grad *= -1; if (m_identifier != "GAFF") grad *= Calc::KCAL_TO_KJ; cleanGradients(grad); // <-- after conversion constraintGradients(x, grad);Consider making these consistent by moving
cleanGradientsbefore the conversion inobenergy.cpp:♻️ Suggested fix for consistent ordering
grad *= -1; // OpenBabel outputs forces, not grads + cleanGradients(grad); // if method is not GAFF, convert to kJ/mol if (m_identifier != "GAFF") grad *= Calc::KCAL_TO_KJ; - cleanGradients(grad); // add in any constraints constraintGradients(x, grad);
Signed-off-by: Geoff Hutchison <geoff.hutchison@gmail.com>
Signed-off-by: Geoff Hutchison <geoff.hutchison@gmail.com>
Signed-off-by: Geoff Hutchison <geoff.hutchison@gmail.com>
|
@coderabbitai review |
✅ Actions performedReview triggered.
|
There was a problem hiding this comment.
Actionable comments posted: 4
🤖 Fix all issues with AI agents
In `@avogadro/qtplugins/autoopt/autoopt.cpp`:
- Around line 111-121: The mass minimum check (setting masses < 0.5 to 1.0) is
applied in one place but not in start() or dynamicsStep(), causing inconsistent
handling of dummy atoms; extract the logic into a helper (e.g.
AutoOpt::initializeMasses) that reads m_molecule->atomCount() and
m_molecule->atom(i).mass() and fills m_masses (size n*3) using a min mass (use
std::max(mass, 1.0) or equivalent) for all three components per atom, then call
this helper from the constructor/initializer, start(), and dynamicsStep()
locations so m_masses is initialized consistently across the class.
- Around line 288-294: The mass initialization for m_masses in autoopt.cpp
currently writes raw atom masses (m_molecule->atom(i).mass()) into the 3N array
without the minimum-mass guard used in setMolecule(); modify the loop that fills
m_masses so each atom mass is checked and if mass < 0.5 it is replaced with 1.0
before assigning to m_masses[i*3], m_masses[i*3+1], and m_masses[i*3+2],
ensuring the same minimum-mass rule as setMolecule() to avoid division-by-zero
in dynamics.
- Around line 456-464: The frozen-atom mask (from
m_molecule->molecule().frozenAtomMask(), stored in mask and passed to
m_method->setMask(mask)) is only used for gradients but not enforced in the
integrator; update the integration code so that after the thermostat/velocity
update you zero velocities for indices where mask[i]==0.0, and after the
position update loop restore original positions for those same mask==0.0 degrees
of freedom (or skip updating them entirely during the position step). Use the
existing mask vector (and m_molecule/m_method context) to identify frozen DOFs
(0.0 = frozen, 1.0 = free) and apply the velocity-zeroing and
position-restoration logic in the integrator where velocities and positions are
updated.
In `@avogadro/qtplugins/autoopt/csvrthermostat.h`:
- Around line 172-175: The CSVR kinetic energy update in csvrthermostat.h is
implemented incorrectly; update the ke_new computation (the line assigning
ke_new) to match Eq. A7 from Bussi et al.: ke_new = c*ke_current + (1.0 -
c)*ke_target + 2.0*sqrt(c*(1.0 - c)*ke_current*ke_target/ n_dof)*R1 + (1.0 -
c)*ke_target/ n_dof * sum_Rsq, using the existing variables ke_new, c,
ke_current, ke_target, n_dof, R1 and sum_Rsq and ensure correct parenthesization
and floating-point promotion to avoid integer division.
♻️ Duplicate comments (1)
avogadro/qtplugins/autoopt/autoopt.cpp (1)
466-474: Apply minimum mass guard to prevent division by zero.This mass initialization doesn't guard against zero masses, which would cause infinity in the acceleration calculation at line 518 (
gradient.array() / m_masses.array()).
🧹 Nitpick comments (2)
avogadro/qtplugins/autoopt/csvrthermostat.h (1)
87-93: Consider documenting the fixed random seed.The default
seed = 12345provides reproducibility but will produce identical trajectories across runs. This is likely intentional for debugging, but worth noting in a comment or considering a time-based seed option for production use.avogadro/qtplugins/autoopt/autoopt.cpp (1)
178-204: Consider tightening timestep range for stability.The timestep range of 0-1000 fs is very permissive. Typical MD simulations use 0.5-2 fs timesteps for stability with hydrogen-containing molecules. Large timesteps (>5 fs) will cause integration instability.
💡 Suggested range adjustment
QDoubleSpinBox* timeStepSpinBox = new QDoubleSpinBox(); timeStepSpinBox->setObjectName("timeStepSpinBox"); - timeStepSpinBox->setRange(0.0, 1000.0); - timeStepSpinBox->setSingleStep(1.0); - timeStepSpinBox->setDecimals(1); + timeStepSpinBox->setRange(0.1, 10.0); + timeStepSpinBox->setSingleStep(0.5); + timeStepSpinBox->setDecimals(2);
Signed-off-by: Geoff Hutchison <geoff.hutchison@gmail.com>
This reverts commit c85ddb8. Signed-off-by: Geoff Hutchison <geoff.hutchison@gmail.com>
|
This pull request has been mentioned on Avogadro Discussion. There might be relevant details there: |
Signed-off-by: Geoff Hutchison <geoff.hutchison@gmail.com>
|
This pull request has been mentioned on Avogadro Discussion. There might be relevant details there: https://discuss.avogadro.cc/t/january-2026-development-updates/7272/1 |
…nh-add-dynamics Signed-off-by: Geoff Hutchison <geoff.hutchison@gmail.com>
Also fix (and double-check) CSVR thermostat equations Now has stable temperatures Signed-off-by: Geoff Hutchison <geoff.hutchison@gmail.com>
ghutchis
changed the title
Also fix some Eigen use based on suggestions from Claude Code Signed-off-by: Geoff Hutchison <geoff.hutchison@gmail.com>
Improved through some suggestions from Claude Opus 4.6 Signed-off-by: Geoff Hutchison <geoff.hutchison@gmail.com>
|
@coderabbitai review |
✅ Actions performedReview triggered.
|
There was a problem hiding this comment.
Actionable comments posted: 6
Caution
Some comments are outside the diff and can’t be posted inline due to platform limitations.
⚠️ Outside diff range comments (2)
avogadro/qtplugins/forcefield/obenergy.cpp (1)
194-216:⚠️ Potential issue | 🟡 MinorMissing null guard on
d->m_obmolbefore use.
value()checksm_molecule == nullptrat Line 196 but does not checkd->m_obmol. IfsetMoleculewas never called (or was called with a valid molecule butd->m_obmolallocation failed), Line 200 will dereference a null pointer. The same applies togradient()at Line 224.The
d->setupflag exists but is also not checked — if the force field setup failed,Energy()may produce meaningless results. Consider guarding ond->m_obmol != nullptr && d->setup:Proposed fix
Real OBEnergy::value(const Eigen::VectorXd& x) { - if (m_molecule == nullptr || m_molecule->atomCount() == 0) + if (m_molecule == nullptr || m_molecule->atomCount() == 0 || d->m_obmol == nullptr || !d->setup) return 0.0; // nothing to do(Same for
gradient().)avogadro/calc/uff.cpp (1)
970-980:⚠️ Potential issue | 🟠 MajorBug:
cosThetadivided by norms twice (double-normalization).At lines 971–973,
ij,ik, andilare normalized to unit vectors. Then line 976 computes:Real cosTheta = ij.dot(ik) / (rij * rik);Since
ijandikare already unit vectors,ij.dot(ik)already yields the cosine. Dividing by(rij * rik)again is incorrect—it effectively divides by the original norms twice, producing a wrong angle and cascading into wrongsinTheta,ratio, and all downstream gradient terms.Currently mitigated because OOP energy/gradient calls are commented out (lines 1165 and 1263), but this will silently produce wrong results when OOP is re-enabled.
🐛 Proposed fix
- Real cosTheta = ij.dot(ik) / (rij * rik); + Real cosTheta = ij.dot(ik);
🤖 Fix all issues with AI agents
In `@avogadro/calc/uff.cpp`:
- Around line 798-803: The bond gradient calculation lacks a near-zero guard for
r and can divide by zero; update the block that computes diff, r, and force (the
lines computing 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;) to skip the
division or use a safe branch when r is below a small epsilon (e.g., if r < eps
set force to zero or use diff scaled by eps) so you avoid r in the denominator
and prevent Inf/NaN propagation; mirror the same near-zero pattern used in the
angle/torsion/OOP gradient helpers.
In `@avogadro/qtplugins/autoopt/autoopt.cpp`:
- Around line 418-420: The DOF calculation uses all atoms
(m_thermostat->setDegreesOfFreedom(3 * n - 3)) but ignores frozen atoms masked
out later; compute the number of frozen atoms from the mask (e.g., int
frozenCount = (mask.array() < 0.5).count() / 3) and pass 3*(n - frozenCount) - 3
to m_thermostat->setDegreesOfFreedom instead, guarding the result to be
non-negative; update the code around m_thermostat, setDegreesOfFreedom, mask and
n to use this adjusted DOF.
- Around line 495-537: The dynamicsStep code must guard against NaN/Inf in
newPositions before committing them to the molecule: check
newPositions.allFinite() (Eigen::VectorXd::allFinite()) right before the
Eigen::Map assignment / m_molecule->setAtomPositions3d call and only perform the
write when true; if not finite, log a warning (e.g. qWarning()) and skip the
position update (optionally revert m_velocities/m_acceleration to their previous
safe state or bail out of the step) so the molecule state is not corrupted by
NaN/Inf values coming from gradient→acceleration→velocity→position propagation.
- Line 283: The moleculeChanged() and methodChanged() paths can call start()
while m_running is true, which triggers another
m_molecule->beginMergeMode("AutoOpt") without a matching endMergeMode(); to fix,
ensure you call stop() before any start() invocation when restarting the
optimizer: update moleculeChanged() and methodChanged() to check m_running (or
unconditionally) call stop() before calling start(), so
beginMergeMode()/endMergeMode() remain paired and nested macros are avoided.
In `@avogadro/qtplugins/autoopt/csvrthermostat.h`:
- Around line 168-176: compute_temperature currently divides by n_dof without
checking for zero; update compute_temperature to guard against n_dof <= 0 (or
nearly zero) before performing the division and return a safe value (e.g., 0.0)
when degrees of freedom is zero to avoid Inf/CRASH in the UI; locate
compute_temperature and the call to compute_kinetic_energy and add a simple
check like if (n_dof <= 0) return 0.0, otherwise compute and return 2.0 * ke /
(n_dof * kB_SI).
- Around line 56-58: The KINETIC_CONVERSION constant has a factor-of-1000 error:
replace the current constexpr double KINETIC_CONVERSION value with 10000.0 (or
1.0e4) so that amu·(Å/fs)² converts correctly to kJ/mol, and update the adjacent
comment to reflect the correct molar conversion (1.66054e-17 J × 6.022e23 =
1.0e7 J/mol = 1.0e4 kJ/mol); reference the KINETIC_CONVERSION symbol in
csvrthermostat.h when making this change.
🧹 Nitpick comments (2)
avogadro/qtplugins/forcefield/obenergy.cpp (1)
232-236: Sign convention is correct. The old per-atom loop implementation usedGetGradient(atom)with the same negation (grad *= -1; // OpenBabel outputs forces, not grads), confirming thatGetGradientPtrcorrectly returns forces and the negation is appropriate for gradient-based optimization.Consider using
d->m_obmol->NumAtoms()instead ofm_molecule->atomCount()at line 233 for defensive alignment with the gradient's source object, though both should work correctly in practice.avogadro/qtplugins/autoopt/csvrthermostat.h (1)
92-100: Considerstd::chi_squared_distributionfor large molecule performance.
sum_noises_squared(n_dof - 1)drawsn_dof - 1individual Gaussians in a loop. For a 1000-atom system this means ~3000 RNG calls per thermostat application.std::chi_squared_distribution<double>(n_dof - 1)uses optimized algorithms (typically Marsaglia/Tsang for large DOF) and would be significantly faster.Suggested replacement
- double sum_noises_squared(int n) - { - double sum = 0.0; - for (int i = 0; i < n; ++i) { - double r = gaussian_random(); - sum += r * r; - } - return sum; - } + double sum_noises_squared(int n) + { + if (n <= 0) + return 0.0; + std::chi_squared_distribution<double> chi2(n); + return chi2(rng); + }
Signed-off-by: Geoff Hutchison <geoff.hutchison@gmail.com>
Signed-off-by: Geoff Hutchison <geoff.hutchison@gmail.com>
|
This pull request has been mentioned on Avogadro Discussion. There might be relevant details there: |
github-project-automation
bot
moved this from Backlog
to Done
in Avogadro 2.1 - 02 June 2026
2 participants
Developer Certificate of Origin
Version 1.1
Copyright (C) 2004, 2006 The Linux Foundation and its contributors.
1 Letterman Drive
Suite D4700
San Francisco, CA, 94129
Everyone is permitted to copy and distribute verbatim copies of this
license document, but changing it is not allowed.
Developer's Certificate of Origin 1.1
By making a contribution to this project, I certify that:
(a) The contribution was created in whole or in part by me and I
have the right to submit it under the open source license
indicated in the file; or
(b) The contribution is based upon previous work that, to the best
of my knowledge, is covered under an appropriate open source
license and I have the right under that license to submit that
work with modifications, whether created in whole or in part
by me, under the same open source license (unless I am
permitted to submit under a different license), as indicated
in the file; or
(c) The contribution was provided directly to me by some other
person who certified (a), (b) or (c) and I have not modified
it.
(d) I understand and agree that this project and the contribution
are public and that a record of the contribution (including all
personal information I submit with it, including my sign-off) is
maintained indefinitely and may be redistributed consistent with
this project or the open source license(s) involved.
Summary by CodeRabbit
New Features
Performance