@@ -371,12 +371,51 @@ void DeepSpinPTExpt::compute(ENERGYVTYPE& ener,
371371 int nloc = nall_real - nghost_real;
372372 int nframes = 1 ;
373373
374- // Build spin tensor for real atoms using bkw_map
375- std::vector<VALUETYPE > dspin (static_cast <size_t >(nall_real) * 3 );
376- for (int ii = 0 ; ii < nall_real; ++ii) {
374+ // Phantom-atom padding for the empty-subdomain corner case
375+ // (``nloc_real == 0``). Multi-rank spin MD can land a rank with zero
376+ // real local atoms when atoms migrate to other subdomains. The
377+ // with-comm AOTI artifact, traced with ``nloc_min=1`` and lowered by
378+ // inductor with an even stricter ``nloc >= 2`` runtime-check
379+ // (silently bypassed because ``AOTI_RUNTIME_CHECK_INPUTS`` is unset by
380+ // default), then SIGFPEs at runtime with an "integer divide by zero"
381+ // inside inductor-generated shape arithmetic that uses ``nloc`` as a
382+ // divisor. The failure is intermittent because inductor re-codegens
383+ // across runs and only some compiles emit the offending divide.
384+ //
385+ // Fix: prepend two phantom atoms with no neighbours so the AOTI graph
386+ // runs with ``nloc == 2``. The phantoms have an empty nlist row and
387+ // therefore contribute zero atomic energy / force / virial, preserving
388+ // the physically-correct "this rank has no real atoms" semantics.
389+ // ``nlocal`` in the comm tensors is set to ``2`` so border_op writes
390+ // received ghost features past the phantom slots; outputs are stripped
391+ // of the phantom prefix before being scattered back to LAMMPS atoms
392+ // via ``select_map``.
393+ const int phantom_n = (nloc_real == 0 && nall_real > 0 ) ? 2 : 0 ;
394+ if (phantom_n > 0 ) {
395+ dcoord.insert (dcoord.begin (), static_cast <size_t >(phantom_n) * 3 ,
396+ static_cast <VALUETYPE >(0 ));
397+ datype.insert (datype.begin (), static_cast <size_t >(phantom_n), 0 );
398+ // Keep aparam_ aligned with the padded local atoms: the phantom atoms
399+ // get zero-valued atomic-parameter rows so the aparam tensor built below
400+ // (shape {1, nloc, daparam}) stays consistent with the padded ``nloc``.
401+ // (aparam_nall is false here, so aparam_ is a per-local-atom buffer.)
402+ if (daparam > 0 ) {
403+ aparam_.insert (aparam_.begin (), static_cast <size_t >(phantom_n) * daparam,
404+ static_cast <VALUETYPE >(0 ));
405+ }
406+ nall_real += phantom_n;
407+ nloc_real = phantom_n;
408+ nloc = nall_real - nghost_real;
409+ }
410+
411+ // Build spin tensor for real atoms using bkw_map (skip phantom prefix
412+ // which keeps zero spin).
413+ std::vector<VALUETYPE > dspin (static_cast <size_t >(nall_real) * 3 ,
414+ static_cast <VALUETYPE >(0 ));
415+ for (int ii = phantom_n; ii < nall_real; ++ii) {
377416 for (int dd = 0 ; dd < 3 ; ++dd) {
378417 dspin[static_cast <size_t >(ii) * 3 + dd] =
379- spin[static_cast <size_t >(bkw_map[ii]) * 3 + dd];
418+ spin[static_cast <size_t >(bkw_map[ii - phantom_n ]) * 3 + dd];
380419 }
381420 }
382421
@@ -445,11 +484,25 @@ void DeepSpinPTExpt::compute(ENERGYVTYPE& ener,
445484 nlist_data.shuffle_exclude_empty (fwd_map);
446485 nlist_data.padding ();
447486
448- // Rebuild mapping tensor
487+ // Rebuild mapping tensor. Phantom slots (when phantom_n > 0) get
488+ // identity entries — they index into their own row and never appear
489+ // in any other atom's nlist (their nlist rows are all -1 below).
449490 if (lmp_list.mapping ) {
450491 std::vector<std::int64_t > mapping (nall_real);
451- for (int ii = 0 ; ii < nall_real; ii++) {
452- mapping[ii] = fwd_map[lmp_list.mapping [bkw_map[ii]]];
492+ for (int ii = 0 ; ii < phantom_n; ii++) {
493+ mapping[ii] = ii;
494+ }
495+ for (int ii = phantom_n; ii < nall_real; ii++) {
496+ // Defensive: this branch (lmp_list.mapping != nullptr) is single-rank
497+ // only (set_mapping is gated on comm->nprocs==1 in pair_deepspin /
498+ // pair_deepmd), while phantom_n>0 only occurs on a multi-rank empty
499+ // subdomain, so the two cannot currently co-occur and the +phantom_n
500+ // term is a no-op (phantom_n==0) on every reachable path. It is kept
501+ // so the mapping stays correct -- resolving fwd_map's pre-padding local
502+ // index into the post-padding local index space -- if that invariant
503+ // ever changes.
504+ mapping[ii] =
505+ fwd_map[lmp_list.mapping [bkw_map[ii - phantom_n]]] + phantom_n;
453506 }
454507 mapping_tensor =
455508 torch::from_blob (mapping.data (), {1 , nall_real}, int_option)
@@ -472,8 +525,16 @@ void DeepSpinPTExpt::compute(ENERGYVTYPE& ener,
472525 }
473526
474527 // Flatten raw nlist — the .pt2 model sorts by distance on-device.
528+ // Phantom rows (all -1) are prepended below so the AOTI graph sees
529+ // nloc == phantom_n + nloc_real_orig instead of 0.
475530 firstneigh_tensor =
476531 createNlistTensor (nlist_data.jlist , nnei).to (torch::kInt64 ).to (device);
532+ if (phantom_n > 0 ) {
533+ auto phantom_rows = torch::full (
534+ {1 , phantom_n, nnei}, static_cast <std::int64_t >(-1 ),
535+ torch::TensorOptions ().dtype (torch::kInt64 ).device (device));
536+ firstneigh_tensor = torch::cat ({phantom_rows, firstneigh_tensor}, 1 );
537+ }
477538 }
478539
479540 // Build fparam/aparam tensors
@@ -566,6 +627,23 @@ void DeepSpinPTExpt::compute(ENERGYVTYPE& ener,
566627 ener.assign (flat_energy_.data_ptr <ENERGYTYPE >(),
567628 flat_energy_.data_ptr <ENERGYTYPE >() + flat_energy_.numel ());
568629
630+ // Zero the reduced energy on an empty rank. Phantoms have constant
631+ // atomic outputs (per-type bias + zero-neighbour MLP) that flow into
632+ // ``energy_redu`` -- and on the spin path the SpinModel doubles atoms
633+ // so the bias contribution appears for both real and spin phantom
634+ // halves; subtracting only the real-half exposed by
635+ // ``output_map["energy"]`` after the ``[:, :nloc]`` slice leaves the
636+ // spin-half leaking into the MPI-reduced LAMMPS total. The physical
637+ // contribution of a rank with no real local atoms is zero by
638+ // definition, so just clear ``ener`` directly.
639+ //
640+ // Forces, force_mag, and virial are unaffected because phantom atomic
641+ // outputs are coord-independent (no neighbours) so their derivatives
642+ // are zero -- no analogous correction is needed.
643+ if (phantom_n > 0 ) {
644+ std::fill (ener.begin (), ener.end (), static_cast <ENERGYTYPE >(0 ));
645+ }
646+
569647 // Extract force: energy_derv_r (nf, nall, 1, 3) -> (nf, nall, 3)
570648 torch::Tensor force_tensor =
571649 output_map[" energy_derv_r" ].squeeze (-2 ).view ({-1 }).to (floatType);
@@ -588,6 +666,17 @@ void DeepSpinPTExpt::compute(ENERGYVTYPE& ener,
588666 virial.assign (cpu_virial_.data_ptr <VALUETYPE >(),
589667 cpu_virial_.data_ptr <VALUETYPE >() + cpu_virial_.numel ());
590668
669+ // Strip the phantom prefix (see phantom-atom padding comment near
670+ // ``select_real_atoms_coord``) so the ``bkw_map`` lookup below sees
671+ // only the real / ghost atoms it was built for. The phantom slots
672+ // carry zero forces because their nlist rows were all -1 — they
673+ // produce no neighbour contributions, so dropping them is exact.
674+ if (phantom_n > 0 ) {
675+ dforce.erase (dforce.begin (), dforce.begin () + phantom_n * 3 );
676+ dforce_mag.erase (dforce_mag.begin (), dforce_mag.begin () + phantom_n * 3 );
677+ nall_real -= phantom_n;
678+ }
679+
591680 // bkw map: map force from real atoms back to full atom list
592681 force.resize (static_cast <size_t >(nframes) * fwd_map.size () * 3 );
593682 force_mag.resize (static_cast <size_t >(nframes) * fwd_map.size () * 3 );
@@ -612,6 +701,16 @@ void DeepSpinPTExpt::compute(ENERGYVTYPE& ener,
612701 cpu_atom_virial_.data_ptr <VALUETYPE >(),
613702 cpu_atom_virial_.data_ptr <VALUETYPE >() + cpu_atom_virial_.numel ());
614703
704+ // Strip the phantom prefix from atomic outputs as well (see force
705+ // block above). Phantom slots carry zero atomic energy / virial
706+ // because their nlist rows were all -1.
707+ if (phantom_n > 0 ) {
708+ datom_energy.erase (datom_energy.begin (),
709+ datom_energy.begin () + phantom_n);
710+ datom_virial.erase (datom_virial.begin (),
711+ datom_virial.begin () + phantom_n * 9 );
712+ }
713+
615714 atom_energy.resize (static_cast <size_t >(nframes) * fwd_map.size ());
616715 atom_virial.resize (static_cast <size_t >(nframes) * fwd_map.size () * 9 );
617716 select_map<VALUETYPE >(atom_energy, datom_energy, bkw_map, 1 , nframes,
0 commit comments