Skip to content

Fix fragment parsing and placement bugs; neighbor lists, load balancing, ligand options - #28

Merged
slipchenko merged 19 commits into
masterfrom
lyuda_bugfix_Oct2026
Oct 20, 2026
Merged

Fix fragment parsing and placement bugs; neighbor lists, load balancing, ligand options#28
slipchenko merged 19 commits into
masterfrom
lyuda_bugfix_Oct2026

Conversation

@slipchenko

@slipchenko slipchenko commented Oct 20, 2026

Copy link
Copy Markdown

Bug fixes, new options and performance work found while running EFP on a 151,550-fragment photosystem I model (polarization did not converge on master). Full details are in CHANGELOG.md.

Bug fixes

  • Polarizable points from CT100 on were misplaced (fixed 4-character label skip; GAMESS writes CT123). Affected every fragment with >99 polarizable points and caused polarization catastrophes.
  • Labels touching negative coordinates (CT1-132.28…, BO123124-132.28…) lost the sign or the coordinate.
  • Dipole-quadrupole points stored at the wrong (out-of-bounds) index.
  • Dispersion crashed for fragments without LMOs (uninitialized overlap arrays).
  • Fragments not in their center-of-mass frame were misplaced with coord atoms/xyzabc, and cutoffs/PBC used the frame origin instead of the molecule.
  • Bond midpoints with atom numbers >99 decoded wrongly (special/flexible fragments).
  • Symmetry class multiplicities (n_symm_frag) were wrong unless the input started with one fragment per class.
  • Memory leaks in the symmetric-crystal polarization.

New options

ligand_cutoff, pairwise_ligand_only, print_geometry, print_pairwise_geometry (efpmd keywords / efp_opts fields), and TIMING lines with print 1. Website documentation: libefp2/libefp2.github.io#31.

Behavior changes reviewers should check

  • Fragments are moved to their center-of-mass frame when parameters are read (no change for MAKEFP potentials; coord points placement unchanged). MD/Hessian warn for non-principal-axes frames.
  • efp_compute(efp, do_gradient) now honors do_gradient (it was forced to 1). Gradient getters fail after an energy-only call. This may matter for interfaces (Q-Chem, Psi4) that read gradients after efp_compute(efp, 0).
  • etest computes the energy without gradients.
  • symm_2, symm_2full, symm_2pw reference energies changed (-2.4694268198 → -5.5786738865): the xxii_exp_n2_* frame origin is 14 bohr from the molecule and was used as the fragment center for the cutoff. Energies without cutoff are unchanged.

Performance

Linear-time parsing of large fragments (493k-point MM fragment: 25 min → 2 s); sparse skip list and cell-list neighbor lists instead of the dense n_frag² skip list; per-thread stress tensor; OpenMP load balancing (polarization over points, two-body over partner blocks, chunked heavy MM pairs); parallel symmetric-crystal path. Example, 151,550 fragments on 64 threads: 15 Å cutoff single point 53 s; EFP/MM with a 30 Å boundary 2.5 min (was 1 h 12 min).

Testing

  • ctest without Torch: 200/200 (master: 156/156); 44 new regression tests, each verified to fail without its fix where applicable.
  • ctest with Torch (LibTorch 2.5.0, macOS arm64): 216/216 incl. all torch_* tests (master: 172/172).
  • pylibefp pytest (CI commands, Psi4 1.10): 60/60 passed, 0 skipped (same as master).
  • Full PSI system on Negishi: energies identical across all performance commits.

🤖 Generated with Claude Code

slipchenko and others added 19 commits October 19, 2026 09:49
@slipchenko @claude
Fix fragment file parsing of GAMESS labels; fix disp crash without LMOs
7312a4c 
Parser (src/parse.c):
- Polarizable points: GAMESS writes the labels as 'CT',I1/I2/I3, so labels
  from CT100 on are five characters long; the fixed 4-character skip read the
  last label digit as x and dropped z. This misplaced the points of every
  fragment with more than 99 polarizable points (chlorophylls, carotenoids)
  and caused polarization catastrophes.
- New skip_ct_label() (skip "CT", blanks, index digits) for static, dynamic
  and dipole-quadrupole polarizable points and LMO centroids. It also handles
  negative coordinates that fill their F15.10 field and run into the label
  ("CT1-132.2805782867"), where the fixed skips ate the sign and the token
  skip of the LMO centroids swallowed the coordinate.
- tok_label(): end an A8 label that runs into a negative number
  ("BO123124-132.2805782867").
- parse_dipquad_polarizable_pts(): store points at
  n_dipquad_polarizable_pts - 1 (was n_dynamic_polarizable_pts - 1, out of
  bounds when no dynamic points were read before).
- Projection basis/wavefunction keep their fixed widths (correct for the
  GAMESS A8,A2 and I2,I3 formats); documented.

Dispersion (src/efp.c, src/disp.c): the LMO overlap arrays s/ds were left
uninitialized for pairs with a fragment without LMOs but were read for every
pair of dynamic polarizable points, also with Tang-Toennies damping
(segmentation fault). Initialize them to NULL and read them only when present.

Tests: parse_labels_1a/b (labels CT100+), parse_labels_2a/b (labels glued to
negative coordinates, all sections), parse_labels_3a/b (dipole-quadrupole
section), disp_no_lmo_1a/b (dispersion without LMOs). The b tests use
fragments generated from pentacene/water by tests/tools/
make_parse_test_fragments.py and must reproduce the energies of the a tests
with the library fragments; all b tests fail without these fixes.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@slipchenko @claude
Center fragment frames at the center of mass; robust bond mid-point a…
01e8c84 
…toms

Center of mass (src/parse.c, src/efp.c):
- Fragments are placed as frag->x + R * p_lib, and frag->x is the center
  of mass for the xyzabc and atoms coordinate types and for cutoffs and
  periodic images. Library fragments written in another frame than the
  center-of-mass frame of GAMESS (e.g. in MD coordinates) were therefore
  misplaced with atom coordinates (all polarizable, dispersion and exchange-
  repulsion points shifted by R * COM_lib), and cutoffs/periodic images were
  computed from the internal origin instead of the molecule. center_fragment()
  now moves every library fragment to its center-of-mass frame when it is
  loaded. Placement with coord points is frame independent and unchanged.
- center_fragment() also records whether the axes are principal axes of
  inertia; efp_get_frag_inertia() (used by MD and the Hessian, which assume
  principal axes) warns otherwise. Fragments are not rotated, because that
  would change the meaning of xyzabc Euler angles.

Bond mid-points (src/parse.c, src/elec.c): efp_update_elec_special() decoded
the atoms of "BO<i><j>" labels by string length, which is ambiguous ("BO1211"
is atoms 12 and 11 or 121 and 1), wrong for atom numbers above 99, and left
the indices uninitialized for longer labels. assign_multipole_atoms() now
finds the atom, or the atom pair of a bond mid-point, of every multipole point
when the fragment is loaded: label splits are checked against the coordinates
(linear in the number of points). Points that are neither move rigidly with
the fragment; before, they were not moved at all.

Tests:
- coord_com_1a/b: atom coordinates with a fragment whose frame is 500 bohr
  away from the molecule (reproduces parse_labels_2a, which uses coord points).
- bo_labels_1a/b: 108-atom fragment with labels like BO10075, atom
  coordinates vs coord points.
Both b tests fail without these changes.
- symm_2, symm_2full, symm_2pw: references updated. The internal origin of
  xxii_exp_n2_*.efp is 14.1 bohr away from the molecule and was the fragment
  center for the cutoff and periodic images; energies without cutoff are
  unchanged, and symmetric and full calculations still agree.
- parse_labels_1a/b, disp_no_lmo_1a/b: smaller num_step_angle, as the center
  of mass of pentacene_2x is 200 bohr away from the interacting molecule.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@slipchenko @claude
Speed up parsing of large fragments (hash index of point labels)
4c49d8a 
The MONOPOLES, DIPOLES, QUADRUPOLES, OCTUPOLES and SCREEN sections give one
line per multipole point, and every line searched all points for its label:
quadratic in the number of points. SCREEN3 lines did not even stop at the
match. COORDINATES grew its arrays by one element per point.

- find_multipole_pt(): open-addressing hash index label -> first point with
  that label (same result as the linear search, including duplicated
  labels), built once per fragment and freed after parsing.
- parse_coordinates(): geometric growth of the atom and point arrays.

A 493,567-point MM superfragment (point charges of a protein environment)
plus a water molecule: 1519 s -> 2.0 s for the whole efpmd run, identical
energy.

Test parse_labels_4a/b: multipole and screening sections listed in reverse
order relative to COORDINATES (h2o_shuffled.efp) must reproduce h2o_l.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@slipchenko @claude
Accumulate the stress tensor only when gradients are computed
22f95d4 
efp_add_stress() adds to one shared tensor with nine omp atomic updates per
fragment pair. The stress tensor can only be read with gradients
(efp_get_stress_tensor), but it was accumulated in every calculation; in a
151,550-fragment single point most of the electrostatics time was spent in
these contended atomics. Skip it without gradients. Runs with gradients
still accumulate it with atomics (a per-thread reduction would remove that
contention too).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@slipchenko @claude
Sparse skipped pairs, neighbor lists for cutoffs, per-thread stress
252ad53 
For large systems the fragment pair bookkeeping dominated: the skip-list
was a dense n_frag x n_frag byte array (23 GB for 151,550 fragments), and
every pair loop (two-body terms, static and induced fields of every
polarization iteration) tested all n_frag^2 pairs against it and against
the cutoff.

Skipped pairs (src/pairs.c): a default rule plus a hash map of the pairs set
with efp_skip_fragments(). For symmetric crystals efp_set_symmlist() sets the
rule "skip pairs without a symmetry-unique fragment" instead of n_frag^2
entries (same pairs as before, also replacing earlier explicit settings).
efp->prepared replaces the use of the skip-list as "efp_prepare was called".

Neighbor lists (src/pairs.c): with a cutoff, efp_compute() builds the list
of fragments within the cutoff of every fragment with cell lists (open
boundaries and orthogonal periodic boxes with at least 3 cells per
direction; brute force otherwise), using the same minimum-image test as
efp_skip_frag_pair(). The lists are sorted and valid only inside
efp_compute(); the pair loops in compute_two_body_range() (same pair
ownership as before), compute_two_body_crystal() and pol.c iterate over
them and keep their own checks. Outside efp_compute() (e.g.
efp_get_electric_field from QM codes) they loop over all fragments as before.

Stress tensor: every thread accumulates its own padded copy during
efp_compute(), summed before the MPI reduction, instead of nine omp atomic
updates of one shared tensor per pair. This matters because efp_compute()
currently always computes gradients (efp->do_gradient = 1).

Checks: 3375 water box with 12 A cutoff (open) and 10 A cutoff (periodic
45 A cube): energies agree with the previous code to 6e-10, all 20250
gradient components to 1e-10 (printed precision); NPT pressures agree to
3e-11 relative. New tests nbr_1 (open boundaries) and nbr_2 (periodic,
cell lists) with 512 water molecules and references from the previous code.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@slipchenko @claude
Honor the do_gradient argument of efp_compute()
b47c7a2 
efp_compute() always set efp->do_gradient = 1, so every energy-only
calculation also computed all gradients (e.g. the polarization gradient of
every polarizable point against all point charges of an MM superfragment).
Use the argument again. All tests pass; note that efp_get_gradient(),
efp_get_point_charge_gradient(), efp_get_stress_tensor() and the atomic
gradients now fail after efp_compute(efp, 0), as they are documented to;
programs that ask for gradients after an energy-only call (to be checked
for the Q-Chem interface) must request them in efp_compute().

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@slipchenko @claude
Timing printout of the main phases; efpmd option print_geometry
5d8cac4 
With print >= 1, efpmd prints the wall time of reading the fragment
parameters, creating/placing the fragments and printing the geometry, and
efp_compute() prints the time of the neighbor lists, the two-body terms and
the polarization phases (static field, induced dipoles, energy, gradient).
The new efpmd option print_geometry (default true) turns off the geometry
printout of large systems.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@slipchenko @claude
efpmd: time the pairwise energy printout (print >= 1)
199c2b9 
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@slipchenko @claude
Option pairwise_ligand_only: two-body terms only for pairs with the l…
524b920 
…igand

In a pairwise analysis (enable_pairwise) only the interactions with the
ligand are needed. With the new option pairwise_ligand_only (efp_opts and
efpmd keyword, default off) the two-body terms (electrostatics, dispersion,
exchange repulsion, charge penetration) are computed only for the fragment
pairs with the ligand; polarization is computed for the whole system as
before (static field from all fragments, all induced dipoles), so the
pairwise polarization energies are unchanged. The total energy is then not
the energy of the system (a note is printed), and gradients are refused.

etest jobs compute the energy without gradients (only the energy is
compared).

Test pairwise_ligand_1: pairwise_1 with the option; the reference is the
sum of the pairwise elec/disp/xr energies of pairwise_1 plus its
polarization energy.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@slipchenko @claude
Option ligand_cutoff: separate cutoff for the pairs with the ligand
7e29dec 
With enable_cutoff, the pairs with the ligand (fragment opts.ligand >= 0)
now use ligand_cutoff (efp_opts and efpmd keyword; 0 = swf_cutoff) and all
other pairs swf_cutoff. This applies to the two-body terms and to the
polarization fields between the ligand and the other fragments, so the
interaction with the ligand can be converged with distance while the
static and induced fields in the environment stay at a smaller cutoff.
Every pair keeps one cutoff for all its interactions, so the polarization
problem stays symmetric. Exchange repulsion keeps xr_cutoff.

efp_pair_cutoff() gives the cutoff of a pair; efp_skip_frag_pair(),
efp_make_swf() and the neighbor lists use it (the ligand is tested
explicitly). ligand_cutoff requires enable_cutoff and a ligand fragment,
not with symmetry, and must fit in the periodic box.

Tests:
- ligand_cutoff_1 (gradient test): 4 water molecules on a line with
  swf_cutoff 10 and ligand_cutoff 16 have the same pairs and switching
  factors as a standard run with swf_cutoff 16, which gives the reference.
- ligand_cutoff_2: 512 water molecules (cell lists), swf_cutoff 7 and
  ligand_cutoff 12; reference from standard runs: E(7 A) minus the ligand
  pairwise energies at 7 A plus those at 12 A.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@slipchenko @claude
efpmd: option print_pairwise_geometry
6a4d635 
With print_pairwise_geometry false, the pairwise energy printout gives one
header line and the energies per fragment and the ligand geometry once,
instead of the coordinates of the fragment and of the ligand for every
fragment (2.6 GB for 151,550 fragments). Default true (unchanged output).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@slipchenko @claude
Parallelize the polarization loops over polarizable points
04595d5 
The static field, the ligand field, the induced dipole iterations and the
polarization gradient were parallelized over fragments, so all polarizable
points of one fragment were computed by one thread. A fragment interacting
with very many others (e.g. the ligand with a large ligand_cutoff) then
serialized the calculation: with ligand_cutoff 100 A in a 151,550-fragment
system the induced dipoles took 536 s. The loops now run over the
polarizable points of the fragment range (efp->pol_pt_frag maps a point to
its fragment); gradient contributions were already accumulated atomically.
The polarization energy loop, which adds to per-fragment pairwise energies,
stays per fragment.

Energies and all 20250 gradient components of a 3375 water box (12 A
cutoff, and 10 A cutoff in a periodic box) agree with the previous code to
6e-10 and 1e-10 (printed precision).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@slipchenko @claude
Load balancing of the two-body loop: partner blocks and heavy-pair ch…
f85c6a0 
…unks

The two-body loop gave every fragment's pairs to one thread. Two cases
serialized the calculation:
- a fragment with very many partners (the ligand with a large
  ligand_cutoff owns ~75,000 pairs);
- pairs with a very large fragment (an MM superfragment of 10^5 - 10^6
  point charges with every inner fragment: 71 min in a 30 A EFP/MM run).

With neighbor lists the loop now runs over work items of up to PAIR_BLOCK
(256) partners of one fragment, with the same pair ownership. The body of
the loop is compute_pair(), unchanged. The electrostatics of pairs with more
than HEAVY_PAIR_POINTS (10^6) multipole point pairs is computed after the
pair loop in chunks of HEAVY_PAIR_CHUNK (2048) points of the larger fragment
(efp_frag_frag_elec_range(), the point range version of
efp_frag_frag_elec()), shared among threads; their pairwise energies are
accumulated atomically. Other terms of such pairs (normally absent for MM
fragments) and pairs with a special fragment stay in the pair loop.

Energies and all 20250 gradient components of a 3375 water box (12 A
cutoff, and 10 A cutoff in a periodic box) agree with the previous code to
6e-10 and 1e-10. New test heavy_pair_1: an MM fragment of 15000 charges
(mm_box.efp) with pentacene and water, pairwise; reference from the code
before chunking; the analytic gradient is unchanged (checked separately).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@slipchenko @claude
Fix the symmetry class multiplicities (n_symm_frag)
0ab1d40 
n_symm_frag() counted the fragments of the class of fragment i for
i = 0 ... nsymm_frag - 1 instead of class i + 1, so the multiplicities that
scale the two-body and polarization energies of symmetric crystals were
wrong unless the first fragments of the input were one per class in class
order. symm_2 starts with two fragments of the same class but has 32
fragments in both classes, which hid the error.

Tests symm_classes_a/b: symm_2 without 8 fragments of the second class
(32 + 24) in two orders with the same symmetry-unique fragments; before
the fix they gave -5.0784182113 and -4.4584671219.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@slipchenko @claude
Parallelize the symmetric-crystal path
6eedcfe 
compute_two_body_crystal(), compute_elec_field_crystal(),
compute_id_crystal() and compute_energy_crystal() were serial. They now
use OpenMP like the regular path: the two-body terms over blocks of
partners of the symmetry-unique fragments, the static field and the
induced dipoles over the polarizable points of the symmetry-unique
fragments, the broadcast of the induced dipoles over the fragments (each
written by one thread, its unique fragment from its symmetry class), and
the pairwise parts over the other fragment of each pair. The class arrays
were not freed in compute_elec_field_crystal(), compute_id_crystal() (every
polarization iteration) and compute_energy_crystal(); they are now.

symm_1, symm_2, symm_2full, symm_2pw, symm_2pw_printout and symm_classes_a
give the same energies as before to all printed digits, and all 384
pairwise/lattice energy lines of the pairwise tests are identical.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@slipchenko @claude
Add CHANGELOG.md for the bugfix branch
125c6b4 
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@slipchenko @claude
Portable wall-clock timer for the timing printouts
55149ae 
The library and efpmd are compiled as strict C99 (C_EXTENSIONS OFF), where
the POSIX clock_gettime() and struct timespec are not declared on Linux;
the Intel CI build failed (GCC only warned). The timers are now in
src/timer.c (efp_wall_time) and efpmd/src/timer.c (efpmd_wall_time), which
use omp_get_wtime() with OpenMP, QueryPerformanceCounter() on Windows, and
otherwise clock_gettime() with _POSIX_C_SOURCE requested in that file only.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@slipchenko @claude
CI: pin einsums 1.1.5 with psi4
1599b67 
The conda-forge psi4 1.11 build (py310hb493b4b_1) allows einsums
>=1.1.5,<1.2, but einsums 1.1.6 lacks einsums::random_engine, so
"import psi4" failed in the MKL jobs (also on master). Checked on Linux:
psi4 1.11 imports with einsums 1.1.5 and fails with 1.1.6.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@slipchenko @claude
CI: pin einsums only where psi4 uses it (not on Windows)
b2892a3 
There is no win-64 einsums package; psi4 on Windows does not depend on it.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@slipchenko

slipchenko commented Oct 20, 2026

Copy link
Copy Markdown
Author

Three commits added after the CI run on this PR:

  • 55149ae Portable wall-clock timer: the timing printout used clock_gettime(), which is not declared on Linux in the strict C99 mode the library is built with (the Intel build failed). Now in src/timer.c / efpmd/src/timer.c (omp_get_wtime() with OpenMP, QueryPerformanceCounter() on Windows, clock_gettime() otherwise).
  • 1599b67, b2892a3 CI: pin einsums=1.1.5 when installing psi4 (not on Windows, where psi4 does not use einsums). The current conda-forge psi4 1.11 build allows einsums 1.1.6, which lacks einsums::random_engine, so import psi4 failed in the MKL jobs; this affects master too. Checked on Linux: psi4 1.11 imports with einsums 1.1.5 and fails with 1.1.6.

All 7 jobs pass now (ctest and pytest on Linux GCC/Intel/Torch, macOS Clang, Windows MinGW).

🤖 Generated with Claude Code

@slipchenko
slipchenko merged commit 2c1d8fe into master Oct 20, 2026
7 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Reviewers

No reviews

Assignees

No one assigned

Labels

None yet

Projects

None yet

Milestone

No milestone

Development

Successfully merging this pull request may close these issues.

1 participant

@slipchenko