Skip to content

feat: add CPU SNAP potential support - #8235

Open
19hello wants to merge 9 commits into
deepmodeling:developfrom
19hello:feature/snap-cpu-develop
Open

19hello wants to merge 9 commits into
deepmodeling:developfrom
19hello:feature/snap-cpu-develop

Conversation

@19hello

@19hello 19hello commented Oct 29, 2026

Copy link
Copy Markdown
Collaborator

Reminder

  • I have read AGENTS.md and docs/developers_guide/agent_governance.md.
  • I have explained why this PR does not need a linked issue.
  • I have added unit tests and a CPU case test.
  • I have listed the exact verification commands and results.
  • I have described user-visible behavior changes, including INPUT usage.
  • I have explained the core-module impact.
  • I have requested the governance exceptions below.

Linked Issue

This PR does not require a linked issue because it adds a standalone CPU SNAP potential implementation and corresponding tests.

Unit Tests and/or Case Tests for my changes

Commands run

cmake --build /tmp/abacus-snap-check2 -j2

ctest --test-dir /tmp/abacus-snap-check2 \
  --output-on-failure \
  -R 'MODULE_MD_SNAP_(model|clebsch_gordan|descriptor|calculator)|MODULE_IO_(read_input_serial|read_item_serial)'

python3 tools/03_code_analysis/agent_governance_check.py --staged

git diff --cached --check

LAMMPS cross-checks were also performed using:

/tmp/lammps-snap-build/lmp

The following small configurations were compared between LAMMPS and ABACUS:

- Ta linear SNAP
- C quadratic SNAP
- In/P chemical SNAP

#### Result summary

- ABACUS build passed.
- 4 SNAP unit tests passed.
- 2 INPUT parsing tests passed.
- Minimal ABACUS SNAP case passed.
- Linear, quadratic, and chemical SNAP small-cell results agreed with LAMMPS within the observed numerical precision.
- git diff --cached --check passed.
- Governance check reported no blocking error. Remaining findings are documented below.

#### Checks not run, with reason

- Full MPI scalability validation was not run because this PR focuses on the CPU SNAP implementation and small regression cases.
- GPU SNAP validation was not run because this PR does not add a GPU SNAP implementation.
- pre-commit run abacus-agent-governance --all-files was not run because pre-commit is not installed in the local environment.

### What's changed?

- Add CPU SNAP energy, force, and virial evaluation to ABACUS.
- Support LAMMPS-style SNAP coefficient and parameter files.
- Support:
  - linear SNAP;
  - quadratic SNAP;
  - multi-element chemical SNAP;
  - bzeroflag;
  - bnormflag;
  - wselfallflag;
  - chemflag;
  - standard cutoff switching through switchflag.
- Integrate SNAP with the ABACUS energy-solver factory and MD workflow.
- Add SNAP model, descriptor, Clebsch–Gordan, calculator, and force tests.
- Add a minimal CPU regression case under tests/04_FF/102_SNAP_minimal.

Users can enable SNAP with:

esolver_type    snap
pot_file        snapcoeff
snap_param_file snapparam

The coefficient file is specified through pot_file, and the SNAP parameter file is specified through snap_param_file.

### Governance Notes

#### INPUT/docs changes

This PR enables the existing esolver_type=snap interface and uses the existing pot_file and snap_param_file INPUT parameters. No new INPUT keyword is introduced.

The minimal regression case documents the required files and usage in:

tests/04_FF/102_SNAP_minimal/README

#### Core module impact

The core-module impact is limited to ESolver:

- add ESolver_SNAP;
- connect SNAP to esolver_factory.cpp;
- evaluate SNAP energy, force, and virial during MD and unit-cell calculations.

There are no changes to HSolver, ElecState, Hamilt, Operator, or Psi.

#### Exceptions requested

1. source/source_esolver/esolver_factory.cpp

   The SNAP solver needs access to the existing running log stream. The implementation passes inp.mdp values explicitly and injects the log stream into ESolver_SNAP. The overall global dependency delta is negative (net_delta=-2), but one GlobalV::ofs_running reference remains for logger injection.
2. source/source_esolver/esolver_snap.h

   <iosfwd> is included because the public constructor stores an std::ofstream&. This is the minimal declaration dependency required by the interface.
3. Documentation synchronization

   The user-facing behavior is documented in the regression-case README and in this PR description. A larger user documentation update is not included because the INPUT interface already exists.

@mohanchen mohanchen added MD & LAM MD and Larege Atomic Models Features Needed The features are indeed needed, and developers should have sophisticated knowledge Machine Learning & AI Issues related to AI labels Oct 29, 2026

@ABACUS-PR-BOT ABACUS-PR-BOT left a comment

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The SNAP MD path mixes ghost and owned force units, and accepted model settings can produce inconsistent or non-finite results. Please correct the force handoff and validate or handle those model settings.

Static source and analytic review at 47543b3 against e5ad530, with separate compatibility inspection of develop 5a7e6a6. No build, tests, numerical/MPI runs or CI diagnostics were executed; the author’s reported checks were not independently verified.

Non-blocking: Correct the PR INPUT notes: snap_param_file is newly added, and snap becomes an accepted esolver_type value. The C++ metadata registrations are already present.

Comment on lines +155 to +157
for (std::size_t i = 0; i < mdcell.ghost_atoms().size(); ++i)
{
mdcell.ghost_atoms()[i].force = all_force[nowned + i];

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Convert ghost forces with the same eV/Angstrom-to-Ry/Bohr factor used for owned atoms below. MD_func immediately returns these values through accumulate_ghost_forces(), which adds them directly to owned forces, so a nonzero periodic-image or cross-rank contribution produces incorrect MD forces. Please cover the ghost-return path with a boundary-crossing regression.

Comment on lines +100 to +102
const double z0 = distance / std::tan(theta0);
const double r0inv = 1.0 / std::sqrt(rsq + z0 * z0);
const double ar = r0inv * z0;

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Handle the finite radial limit at distance == rmin0 in both descriptor paths. A positive rmin0 is accepted, but this expression then forms r/tan(0) followed by 0*infinity, producing non-finite higher-order components; the zero-distance guard does not cover it. Please add value-and-force checks at and near a positive rmin0.

Comment on lines +307 to +310
maximum_cutoff = std::max(maximum_cutoff, 2.0 * parameters.element_radii[i] * parameters.rcutfac);
}
parameters.rcut = maximum_cutoff;
return parameters.valid(error);

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Validate rmin0 against every active pair cutoff. Radii 0.5 and 2 with rcutfac=1 and rmin0=1.5 pass this maximum-cutoff check, but the small-small cutoff is 1: the value path retains a neighbor at distance 1.2 while the derivative path excludes it. Please reject invalid pair intervals or apply one consistent radial-domain policy.

Comment on lines +296 to +298
parameters.chem_flag = value != 0;
}
}

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Reject unsupported semantic SNAP parameters explicitly. A file containing switchinnerflag 1 plus sinner/dinner is accepted here, but those settings never reach the descriptor, so the requested short-range switching is omitted from its energy and forces. Please diagnose unsupported options, or implement their value and derivative behavior.

found[element_index] = true;
parameters.element_radii[element_index] = radius;
parameters.element_weights[element_index] = weight;
coefficients[element_index] = values;

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Non-blocking: Document the chemical-channel order expected in STRU when importing a model, and add a test of that contract. Coefficient rows are matched by name, but channel axes follow STRU order, so changing species order can reinterpret the same flat coefficient vector.

Comment on lines +160 to +161
const ModuleMD::SNAP::ForceResult analytic = calculator.compute_with_analytic_force(neighbors, 0);
const ModuleMD::SNAP::ForceResult finite = calculator.compute_with_finite_difference_force(neighbors, 0, 1.0e-6);

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Non-blocking: Add checked-in force comparisons for a non-collinear, multi-neighbor model with twojmax >= 2 and nonzero angular coefficients, plus a quadratic model and a strain check for the virial. The current calculator finite-difference tests use twojmax=0, while the fixed descriptor reference checks values only.

{
const UWithDerivative u_derivative =
descriptor_evaluator_.compute_u_with_derivative(neighbors[neighbor_index], direction);
std::vector<std::complex<double> > u = descriptor_evaluator_.accumulate_u(neighbors, central_element);

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Non-blocking: Accumulate u once before these loops and reuse it for each derivative. The environment is unchanged, but accumulate_u traverses all neighbors on each of the 3*N iterations, adding redundant recurrences and allocations to the force path.

return false;
}

std::ifstream coeff_stream(coefficient_file.c_str());

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Non-blocking: Record the full coefficient and parameter paths after successful reads in stdout and the running log, as required by the file-I/O rule. Only failures identify the files currently; the existing injected running_log can support successful-read confirmation in the solver wrapper.

Comment on lines +149 to +150
snap_virial_(row, column) += relative[row] * pair_force[column]
/ (ModuleBase::Ry_to_eV * mdcell.omega());

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Non-blocking: Move the common cell-level force and virial assembly into module_snap and call it as a named solver step, following the ESolver module rule. The MDCell and UnitCell loops currently duplicate physical assembly and unit conversion outside the calculator tests.

Comment on lines +120 to +122
const Dual x(neighbor.x, direction == 0 ? 1.0 : 0.0);
const Dual y(neighbor.y, direction == 1 ? 1.0 : 0.0);
const Dual z(neighbor.z, direction == 2 ? 1.0 : 0.0);

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Non-blocking: Assign these conditional derivative seeds to named locals before passing them to Dual, following the inspectable-argument rule. Apply the same cleanup to the new numerical-call conversions, and move the computed twojmax member clamp into the constructor body.

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

Features Needed The features are indeed needed, and developers should have sophisticated knowledge Machine Learning & AI Issues related to AI MD & LAM MD and Larege Atomic Models

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants