Merge branch 'main' into fix-alchmtk-compute-distances

Author: shyuep

.github/workflows/lammps-build.yml

@@ -169,9 +169,9 @@ jobs:
           import torch
           from pymatgen.core import Lattice, Structure
           from pymatgen.optimization.neighbors import find_points_in_spheres
-          from matgl.apps._pes_pyg import Potential
-          from matgl.ext._lammps import LAMMPSMatGLModel
-          from matgl.models._tensornet_pyg import TensorNet
+          from matgl.apps.pes import Potential
+          from matgl.ext.lammps import LAMMPSMatGLModel
+          from matgl.models._tensornet import TensorNet
 
           torch.manual_seed(0)
           m = TensorNet(
@@ -310,13 +310,28 @@ jobs:
               sys.exit(1)
           PY
 
+      - name: Validate pressure sign by finite differences
+        shell: bash
+        env:
+          LIBTORCH: ${{ steps.libtorch.outputs.path }}
+        run: |
+          export LD_LIBRARY_PATH="${LIBTORCH}/lib:${LD_LIBRARY_PATH:-}"
+          cp lammps_artifacts/model.pt matgl/lammps/tests/model.pt
+          cd matgl/lammps/tests
+          LMP="${GITHUB_WORKSPACE}/lammps_build/lmp"
+          "${LMP}" -var scale 0.99 -log log.strain-minus.lammps -in in.matgl_finite_strain
+          "${LMP}" -var scale 1.00 -log log.strain-zero.lammps  -in in.matgl_finite_strain
+          "${LMP}" -var scale 1.01 -log log.strain-plus.lammps  -in in.matgl_finite_strain
+          python3 check_finite_difference_pressure.py \
+            log.strain-minus.lammps log.strain-zero.lammps log.strain-plus.lammps
+
       - name: Upload artifacts on failure
         if: failure()
         uses: actions/upload-artifact@v4
         with:
           name: lammps-debug
           path: |
-            matgl/lammps/tests/log.lammps
+            matgl/lammps/tests/log*.lammps
             lammps_artifacts/
             lammps_build/CMakeFiles/CMakeOutput.log
             lammps_build/CMakeFiles/CMakeError.log
@@ -443,9 +458,9 @@ jobs:
           import torch
           from pymatgen.core import Lattice, Structure
           from pymatgen.optimization.neighbors import find_points_in_spheres
-          from matgl.apps._pes_pyg import Potential
-          from matgl.ext._lammps import LAMMPSMatGLModel
-          from matgl.models._tensornet_pyg import TensorNet
+          from matgl.apps.pes import Potential
+          from matgl.ext.lammps import LAMMPSMatGLModel
+          from matgl.models._tensornet import TensorNet
 
           torch.manual_seed(0)
           m = TensorNet(
@@ -557,7 +572,9 @@ jobs:
           export LD_LIBRARY_PATH="${LIBTORCH}/lib:${LD_LIBRARY_PATH:-}"
           cp lammps_artifacts/model.pt matgl/lammps/tests/model.pt
           cd matgl/lammps/tests
-          "${GITHUB_WORKSPACE}/lammps_build/lmp" -k on g 1 -sf kk -in in.matgl_si | tee log.lammps
+          "${GITHUB_WORKSPACE}/lammps_build/lmp" \
+            -k on g 1 -sf kk -pk kokkos neigh half newton on \
+            -in in.matgl_si | tee log.lammps
 
       - name: Diff LAMMPS energy against Python reference
         shell: bash
@@ -589,13 +606,31 @@ jobs:
               sys.exit(1)
           PY
 
+      - name: Validate Kokkos pressure sign by finite differences
+        shell: bash
+        env:
+          LIBTORCH: ${{ steps.libtorch.outputs.path }}
+        run: |
+          export LD_LIBRARY_PATH="${LIBTORCH}/lib:${LD_LIBRARY_PATH:-}"
+          cp lammps_artifacts/model.pt matgl/lammps/tests/model.pt
+          cd matgl/lammps/tests
+          LMP="${GITHUB_WORKSPACE}/lammps_build/lmp"
+          "${LMP}" -k on g 1 -sf kk -pk kokkos neigh half newton on \
+            -var scale 0.99 -log log.strain-minus.lammps -in in.matgl_finite_strain
+          "${LMP}" -k on g 1 -sf kk -pk kokkos neigh half newton on \
+            -var scale 1.00 -log log.strain-zero.lammps -in in.matgl_finite_strain
+          "${LMP}" -k on g 1 -sf kk -pk kokkos neigh half newton on \
+            -var scale 1.01 -log log.strain-plus.lammps -in in.matgl_finite_strain
+          python3 check_finite_difference_pressure.py \
+            log.strain-minus.lammps log.strain-zero.lammps log.strain-plus.lammps
+
       - name: Upload artifacts on failure
         if: failure()
         uses: actions/upload-artifact@v4
         with:
           name: lammps-kokkos-debug
           path: |
-            matgl/lammps/tests/log.lammps
+            matgl/lammps/tests/log*.lammps
             lammps_artifacts/
             lammps_build/CMakeFiles/CMakeOutput.log
             lammps_build/CMakeFiles/CMakeError.log

lammps/src/KOKKOS/pair_matgl_kokkos.cpp

@@ -130,8 +130,7 @@ void PairMATGLKokkos<DeviceType>::compute(int eflag, int vflag)
   atomKK->sync(execution_space, datamask_read);
   atomKK->modified(execution_space, datamask_modify);
 
-  using AT_ = typename AT::t_x_array;  // (nall, 3) double on the device
-  AT_ x = atomKK->k_x.template view<DeviceType>();
+  auto x = atomKK->k_x.template view<DeviceType>();
   auto f = atomKK->k_f.template view<DeviceType>();
   auto type = atomKK->k_type.template view<DeviceType>();
 
@@ -173,13 +172,15 @@ void PairMATGLKokkos<DeviceType>::compute(int eflag, int vflag)
 
   // Fill Z + mask from atom type.
   const auto type_to_z = d_type_to_z_;
+  const auto d_atomic_numbers = d_atomic_numbers_;
+  const auto d_local_or_ghost = d_local_or_ghost_;
   Kokkos::parallel_for(
       "matgl_kk:fill_atoms",
       Kokkos::RangePolicy<DeviceType>(0, nall),
       KOKKOS_LAMBDA(const int i) {
         const int t = type(i);
-        d_atomic_numbers_(i) = type_to_z(t);
-        d_local_or_ghost_(i) = (i < nlocal);
+        d_atomic_numbers(i) = type_to_z(t);
+        d_local_or_ghost(i) = (i < nlocal);
       });
 
   // local_row_of_(j) = the owned row representing the same physical atom as
@@ -205,6 +206,7 @@ void PairMATGLKokkos<DeviceType>::compute(int eflag, int vflag)
   //    so initialize numneigh_short_ for ghost atoms to zero).
   Kokkos::deep_copy(d_numneigh_short_, 0);
   const double r_max_sq = r_max_squared_;
+  const auto d_numneigh_short = d_numneigh_short_;
 
   Kokkos::parallel_for(
       "matgl_kk:count_neigh",
@@ -224,16 +226,17 @@ void PairMATGLKokkos<DeviceType>::compute(int eflag, int vflag)
           const double rsq = dx * dx + dy * dy + dz * dz;
           if (rsq <= r_max_sq) ++nshort;
         }
-        d_numneigh_short_(i) = nshort;
+        d_numneigh_short(i) = nshort;
       });
 
   // 3) Exclusive prefix-sum into d_first_edge_ (length nall+1).
+  const auto d_first_edge = d_first_edge_;
   Kokkos::parallel_scan(
       "matgl_kk:scan_edges",
       Kokkos::RangePolicy<DeviceType>(0, nall + 1),
       KOKKOS_LAMBDA(const int i, int &update, const bool final) {
-        const int v = (i < nall) ? d_numneigh_short_(i) : 0;
-        if (final) d_first_edge_(i) = update;
+        const int v = (i < nall) ? d_numneigh_short(i) : 0;
+        if (final) d_first_edge(i) = update;
         update += v;
       });
 
@@ -256,6 +259,8 @@ void PairMATGLKokkos<DeviceType>::compute(int eflag, int vflag)
   //    need one consistent row per physical atom (periodicity goes through
   //    unit_shifts, not through ghost-row duplication; see pair_matgl.cpp).
   const auto d_local_row_of = d_local_row_of_;
+  const auto d_edge_index = d_edge_index_;
+  const auto d_unit_shifts = d_unit_shifts_;
   const double *const h_inv_host = domain->h_inv;
   const double hinv0 = h_inv_host[0], hinv1 = h_inv_host[1], hinv2 = h_inv_host[2];
   const double hinv3 = h_inv_host[3], hinv4 = h_inv_host[4], hinv5 = h_inv_host[5];
@@ -268,7 +273,7 @@ void PairMATGLKokkos<DeviceType>::compute(int eflag, int vflag)
         const double yi = x(i, 1);
         const double zi = x(i, 2);
         const int jnum = d_numneigh(i);
-        int e = d_first_edge_(i);
+        int e = d_first_edge(i);
         for (int jj = 0; jj < jnum; ++jj) {
           const int j = d_neighbors(i, jj) & NEIGHMASK;
           const double dx = x(j, 0) - xi;
@@ -289,11 +294,11 @@ void PairMATGLKokkos<DeviceType>::compute(int eflag, int vflag)
           const double ly = hinv1 * ddy + hinv3 * ddz;
           const double lx = hinv0 * ddx + hinv5 * ddy + hinv4 * ddz;
 
-          d_edge_index_(0, e) = i;
-          d_edge_index_(1, e) = j_local;
-          d_unit_shifts_(e, 0) = static_cast<int64_t>(Kokkos::round(lx));
-          d_unit_shifts_(e, 1) = static_cast<int64_t>(Kokkos::round(ly));
-          d_unit_shifts_(e, 2) = static_cast<int64_t>(Kokkos::round(lz));
+          d_edge_index(0, e) = i;
+          d_edge_index(1, e) = j_local;
+          d_unit_shifts(e, 0) = static_cast<int64_t>(Kokkos::round(lx));
+          d_unit_shifts(e, 1) = static_cast<int64_t>(Kokkos::round(ly));
+          d_unit_shifts(e, 2) = static_cast<int64_t>(Kokkos::round(lz));
           ++e;
         }
       });

lammps/src/ML-MATGL/pair_matgl.cpp

@@ -496,17 +496,17 @@ void PairMATGL::compute(int eflag, int vflag)
     f[i][2] += fa[i][2];
   }
 
-  // 6) Virial — the model returns a 3x3 tensor with the LAMMPS sign
-  //    convention (V_ij = sum r_i F_j). LAMMPS stores 6 Voigt components
-  //    in `virial`: xx, yy, zz, xy, xz, yz.
+  // 6) Virial — the model returns dE/dstrain. LAMMPS uses
+  //    W = sum_i r_i (x) f_i = -dE/dstrain. Store the six Voigt components
+  //    in the order xx, yy, zz, xy, xz, yz.
   if (vflag_global) {
     auto vir_t = out.at("virials").toTensor().to(torch::kFloat64);
     auto va = vir_t.accessor<double, 2>();
-    virial[0] += va[0][0];
-    virial[1] += va[1][1];
-    virial[2] += va[2][2];
-    virial[3] += 0.5 * (va[0][1] + va[1][0]);
-    virial[4] += 0.5 * (va[0][2] + va[2][0]);
-    virial[5] += 0.5 * (va[1][2] + va[2][1]);
+    virial[0] -= va[0][0];
+    virial[1] -= va[1][1];
+    virial[2] -= va[2][2];
+    virial[3] -= 0.5 * (va[0][1] + va[1][0]);
+    virial[4] -= 0.5 * (va[0][2] + va[2][0]);
+    virial[5] -= 0.5 * (va[1][2] + va[2][1]);
   }
 }

lammps/tests/check_finite_difference_pressure.py

@@ -0,0 +1,89 @@
+"""Validate LAMMPS hydrostatic pressure against a finite energy derivative."""
+
+from __future__ import annotations
+
+import argparse
+import math
+from pathlib import Path
+
+
+EV_PER_ANGSTROM3_TO_BAR = 1.602176634e6
+
+
+def read_thermo(path: Path) -> dict[str, float]:
+    """Return the last thermo row containing PotEng and pressure columns."""
+    header: list[str] | None = None
+    result: dict[str, float] | None = None
+
+    for line in path.read_text().splitlines():
+        fields = line.split()
+        if fields and fields[0] == "Step" and "PotEng" in fields and "Press" in fields:
+            header = fields
+            continue
+        if header is None or len(fields) != len(header):
+            continue
+        try:
+            values = [float(value) for value in fields]
+        except ValueError:
+            continue
+        result = dict(zip(header, values))
+
+    if result is None:
+        raise ValueError(f"Could not find a pressure thermo row in {path}")
+    return result
+
+
+def main() -> None:
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("minus", type=Path, help="LAMMPS log at scale 1-epsilon")
+    parser.add_argument("zero", type=Path, help="LAMMPS log at scale 1")
+    parser.add_argument("plus", type=Path, help="LAMMPS log at scale 1+epsilon")
+    parser.add_argument("--epsilon", type=float, default=0.01)
+    parser.add_argument("--relative-tolerance", type=float, default=0.05)
+    parser.add_argument("--absolute-tolerance-bar", type=float, default=10.0)
+    args = parser.parse_args()
+
+    minus = read_thermo(args.minus)
+    zero = read_thermo(args.zero)
+    plus = read_thermo(args.plus)
+
+    # lambda uniformly scales all three cell vectors and fractional atomic
+    # coordinates.  At lambda=1, dE/dlambda is the trace of dE/dstrain.
+    # LAMMPS pressure is -trace(dE/dstrain)/(3V).
+    energy_derivative = (plus["PotEng"] - minus["PotEng"]) / (2 * args.epsilon)
+    finite_difference_pressure = (
+        -energy_derivative / (3 * zero["Volume"]) * EV_PER_ANGSTROM3_TO_BAR
+    )
+    diagonal_pressure = (zero["Pxx"] + zero["Pyy"] + zero["Pzz"]) / 3
+
+    print(f"E(1-epsilon) = {minus['PotEng']:.12g} eV")
+    print(f"E(1+epsilon) = {plus['PotEng']:.12g} eV")
+    print(f"dE/dlambda   = {energy_derivative:.12g} eV")
+    print(f"FD pressure  = {finite_difference_pressure:.12g} bar")
+    print(f"LAMMPS press = {zero['Press']:.12g} bar")
+
+    if not math.isclose(
+        zero["Press"],
+        diagonal_pressure,
+        rel_tol=1e-10,
+        abs_tol=1e-6,
+    ):
+        raise AssertionError(
+            f"Press ({zero['Press']}) does not equal mean diagonal pressure "
+            f"({diagonal_pressure})"
+        )
+
+    if not math.isclose(
+        zero["Press"],
+        finite_difference_pressure,
+        rel_tol=args.relative_tolerance,
+        abs_tol=args.absolute_tolerance_bar,
+    ):
+        raise AssertionError(
+            "LAMMPS pressure does not match the finite-difference energy "
+            f"derivative: {zero['Press']} vs {finite_difference_pressure} bar"
+        )
+
+
+if __name__ == "__main__":
+    main()

lammps/tests/in.matgl_finite_strain

@@ -0,0 +1,40 @@
+# in.matgl_finite_strain
+#
+# Single-point energy and pressure for a uniformly scaled version of the
+# four-atom Mo-S fixture in in.matgl_si.  The scale must be supplied on the
+# command line, for example:
+#
+#   lmp -var scale 0.99 -in in.matgl_finite_strain
+#
+# Running at 1-epsilon, 1, and 1+epsilon lets
+# check_finite_difference_pressure.py validate the pressure/virial sign from
+# the energy derivative without using the model-reported stress as reference.
+
+units           metal
+atom_style      atomic
+boundary        p p p
+
+atom_modify     map yes
+newton          on
+
+lattice         custom ${scale} a1 4.5 0.0 0.0 a2 0.0 4.5 0.0 a3 0.0 0.0 4.5 &
+                origin 0 0 0 &
+                basis 0.00 0.00 0.00 &
+                basis 0.50 0.50 0.50 &
+                basis 0.50 0.00 0.25 &
+                basis 0.00 0.50 0.75
+region          box prism 0 1 0 1 0 1 0 0 0 units lattice
+create_box      2 box
+create_atoms    1 box basis 1 1 basis 2 2 basis 3 1 basis 4 2
+
+mass            1 95.95
+mass            2 32.06
+
+pair_style      matgl
+pair_coeff      * * model.pt Mo S
+
+thermo_style    custom step pe press pxx pyy pzz vol
+thermo_modify   format float %.15g
+thermo          1
+
+run             0