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

.github/workflows/ci.yml

@@ -296,7 +296,15 @@ jobs:
       # psi4 with MKL avail for qmefp test
       - name: Install Psi4 for QM/EFP testing
         if: ${{ matrix.cfg.blas == 'MKL' }}
-        run: conda install psi4 -c conda-forge
+        # einsums 1.1.6 lacks a symbol of the conda-forge psi4 1.11 build (import fails with
+        # undefined symbol einsums::random_engine); pin einsums until psi4 is rebuilt.
+        # psi4 on Windows does not use einsums (no win-64 package).
+        run: |
+          if [[ "${{ runner.os }}" == "Windows" ]]; then
+            conda install psi4 -c conda-forge
+          else
+            conda install psi4 einsums=1.1.5 -c conda-forge
+          fi
 
       - name: Test (pytest) -- unit tests Python bindings
         run: |

CHANGELOG.md

@@ -0,0 +1,78 @@
+# Changelog
+
+## Unreleased (branch `lyuda_bugfix_Oct2026`)
+
+### Bug fixes
+
+- **Polarizable points from CT100 on were misplaced.** GAMESS writes polarizable point labels as `CT1`,
+  `CT12`, `CT123`; LibEFP skipped a fixed 4 characters, so from the 100th point on the last label digit was
+  read as x and z was lost. Every fragment with more than 99 polarizable points (chlorophylls, carotenoids,
+  large ligands) was affected, typically causing a polarization catastrophe.
+- **Labels touching negative coordinates.** A coordinate that fills its F15.10 field runs into the label
+  (`CT1-132.2805782867`, `BO123124-132.2805782867`); the sign was lost or the coordinate was read as part of
+  the label. Fixed for polarizable points, dynamic and dipole-quadrupole polarizable points, LMO centroids,
+  coordinates and multipole sections.
+- **Dipole-quadrupole polarizable points** were stored at the wrong index (out of bounds when no dynamic
+  polarizable points had been read).
+- **Dispersion crashed for fragments without LMOs** (no exchange-repulsion data): uninitialized overlap
+  arrays were read also with Tang-Toennies damping.
+- **Fragments not written in their center-of-mass frame were misplaced** with atom coordinates (`coord atoms`)
+  and `xyzabc`, and cutoffs and periodic images were computed from the origin of the fragment frame instead
+  of the molecule. See *Behavior changes*.
+- **Bond midpoints with atom numbers above 99** were decoded from the label by its length (ambiguous, wrong
+  for `BO10074`, uninitialized indices for long labels) when atom positions were given (special/flexible
+  fragments). Points that are neither atoms nor bond midpoints were not moved at all.
+- **Symmetry class multiplicities** (`n_symm_frag`) were taken from the classes of the first fragments of the
+  input, so symmetric-crystal energies were wrong unless the input started with one fragment per class.
+- Memory leaks in the symmetric-crystal polarization (every iteration).
+
+### New options
+
+| efpmd keyword / `efp_opts` field | Default | Description |
+|---|---|---|
+| `ligand_cutoff` | 0 (= `swf_cutoff`) | Separate cutoff for the pairs with the ligand (two-body terms and polarization fields); all other pairs use `swf_cutoff`. Requires `enable_cutoff` and `ligand` >= 0; not with symmetry. |
+| `pairwise_ligand_only` | false | With `enable_pairwise`: two-body terms only for the pairs with the ligand; polarization of the whole system as usual. The total energy is then not the energy of the system; no gradients. |
+| `print_geometry` (efpmd) | true | `false` skips the geometry printout (large systems). |
+| `print_pairwise_geometry` (efpmd) | true | `false` prints one line and the energies per fragment in the pairwise analysis, the ligand geometry once. |
+
+With `print 1` efpmd now prints `TIMING` lines for the main phases (reading parameters, creating fragments,
+neighbor lists, two-body terms, polarization: static field, induced dipoles, energy, gradient, pairwise
+printout).
+
+### Behavior changes
+
+- **Fragments are moved to their center-of-mass frame when the parameters are read.** The `.efp` file may use
+  any frame (e.g., the frame of an MD snapshot); placement with `coord points` is unchanged. `xyzabc`
+  positions refer to the center of mass, and cutoffs and periodic images are measured between centers of
+  mass. For fragments written by MAKEFP (center of mass at the origin) nothing changes.
+- MD and the Hessian assume that the fragment frame axes are principal axes of inertia; efpmd now warns for
+  fragments in other frames (fragments are not rotated, which would change the meaning of `xyzabc` Euler
+  angles).
+- **`efp_compute(efp, do_gradient)` honors `do_gradient`** (it used to compute gradients always).
+  `efp_get_gradient()`, `efp_get_point_charge_gradient()`, `efp_get_stress_tensor()` and the atomic gradients
+  return an error after `efp_compute(efp, 0)`; programs that read gradients must request them.
+- `run_type etest` computes the energy without gradients.
+- Reference energies of `symm_2`, `symm_2full`, `symm_2pw` changed (-2.4694268198 -> -5.5786738865): the
+  internal origin of `xxii_exp_n2_*.efp` is 14.1 bohr away from the molecule and was used as the fragment
+  center for the cutoff and periodic images. Energies without cutoff are unchanged.
+
+### Performance
+
+- Fragment files with many points are parsed in linear time (hash index of labels instead of a search per
+  line): a 493,567-point MM fragment loads in 2 s instead of 25 min.
+- The dense n_frag x n_frag skip list (23 GB for 1.5e5 fragments) is replaced by a default rule plus a hash
+  map of explicitly skipped pairs; with a cutoff, `efp_compute()` builds neighbor lists (cell lists) and all
+  pair loops run over them.
+- Stress tensor accumulated per thread instead of with atomic updates per pair.
+- OpenMP load balancing: the polarization loops run over polarizable points; the two-body loop runs over
+  blocks of partners; electrostatics of pairs with more than 10^6 point pairs (MM superfragments) is computed
+  in chunks of points. The symmetric-crystal path is parallelized as well.
+- Example (151,550 fragments, 64 threads, single point, elec + pol): 15 A cutoff 53 s (21.5 min after the
+  parser fixes alone; on the original master the polarization did not converge); EFP/MM with a 30 A boundary
+  2.5 min (1 h 12 min before the load balancing).
+
+### Tests
+
+New: `parse_labels_1a/b ... 4a/b`, `disp_no_lmo_1a/b`, `coord_com_1a/b`, `bo_labels_1a/b`, `nbr_1`, `nbr_2`,
+`pairwise_ligand_1`, `ligand_cutoff_1`, `ligand_cutoff_2`, `heavy_pair_1`, `symm_classes_a/b`. Test fragments
+are generated from the fraglib potentials by `tests/tools/make_parse_test_fragments.py`.

CMakeLists.txt

@@ -184,11 +184,13 @@ set(raw_sources_list
   electerms.c
   int.c
   log.c
+  pairs.c
   parse.c
   pol.c
   poldirect.c
   stream.c
   swf.c
+  timer.c
   util.c
   xr.c
   )

efpmd/CMakeLists.txt

@@ -44,6 +44,7 @@ if(TORCH_SWITCH)
     src/parse.c
     src/rand.c
     src/sp.c
+    src/timer.c
     src/torch.c
   )
 else()
@@ -64,6 +65,7 @@ else()
     src/parse.c
     src/rand.c
     src/sp.c
+    src/timer.c
   )
 endif()
 

efpmd/src/common.c

@@ -370,6 +370,18 @@ void print_pair_energy(struct state *state) {
     char frag_name[32];
     size_t frag_atoms;
     double lattice_energy[6] = {0.0};
+    /* without geometries every fragment block is one line plus the energies;
+     * the ligand geometry is printed once */
+    bool print_geometries = cfg_get_bool(state->cfg, "print_pairwise_geometry");
+
+    if (!print_geometries) {
+        msg("ligand %zu fragment %s\n", ligand_index, ligand);
+        for (size_t a = 0; a < lig_atoms; a++)
+            msg("   %-16s %12.6lf %12.6lf %12.6lf\n", latoms[a].label,
+                latoms[a].x * BOHR_RADIUS, latoms[a].y * BOHR_RADIUS,
+                latoms[a].z * BOHR_RADIUS);
+        msg("\n");
+    }
 
 	for (size_t i=0; i <n_frags; i++){
         check_fail(efp_get_frag_name(state->efp, i, sizeof(frag_name),frag_name));
@@ -379,6 +391,7 @@ void print_pair_energy(struct state *state) {
         check_fail(efp_get_frag_atoms(state->efp, i, frag_atoms, atoms));
 
         msg("   PAIRWISE ENERGY BETWEEN FRAGMENT %zu (%s) AND LIGAND %zu (%s) \n", i, frag_name, ligand_index, ligand);
+        if (print_geometries) {
         msg("fragment %s\n", frag_name);
         for (size_t a = 0; a < frag_atoms; a++) {
             double x = atoms[a].x * BOHR_RADIUS;
@@ -394,10 +407,11 @@ void print_pair_energy(struct state *state) {
             double y = latoms[a].y * BOHR_RADIUS;
             double z = latoms[a].z * BOHR_RADIUS;
             msg("   %-16s %12.6lf %12.6lf %12.6lf\n", latoms[a].label, x, y, z);
+        }
         }
 
 		free(atoms);
-        
+
 		msg("\n");
         msg("    PAIRWISE ENERGY COMPONENTS (ATOMIC UNITS)\n");
         msg("%40s %16.10lf\n", "PAIRWISE ELECTROSTATIC ENERGY", energies[i].electrostatic);

efpmd/src/common.h

@@ -123,6 +123,7 @@ void *xrealloc(void *, size_t);
 
 void print_vec(const double *);
 void print_geometry(struct efp *);
+double efpmd_wall_time(void); /* timer.c */
 void print_geometry_pbc(struct efp *, int);
 void print_energy(struct state *);
 void print_gradient(struct state *);

efpmd/src/energy.c

@@ -42,7 +42,7 @@ void compute_energy(struct state *state, bool do_grad)
 	size_t ifrag, nfrag, iatom, natom, spec_frag, n_special_atoms;
 	int itotal;
 
-    if (cfg_get_int(state->cfg, "print")>0) {
+    if (cfg_get_int(state->cfg, "print") > 0 && cfg_get_bool(state->cfg, "print_geometry")) {
         print_geometry(state->efp);
     }
 

efpmd/src/gtest.c

@@ -397,7 +397,8 @@ void sim_etest(struct state *state)
     msg("ENERGY TEST JOB\n\n\n");
 
     print_geometry(state->efp);
-    compute_energy(state, 1);
+    /* only the energy is tested */
+    compute_energy(state, 0);
     print_energy(state);
     test_energy(state);
 

efpmd/src/main.c

@@ -136,6 +136,7 @@ static struct cfg *make_cfg(void)
 	cfg_add_string(cfg, "efp_params_file", "params.efp");
 	cfg_add_bool(cfg, "enable_cutoff", false);
 	cfg_add_double(cfg, "swf_cutoff", 10.0);
+	cfg_add_double(cfg, "ligand_cutoff", 0.0);
     cfg_add_double(cfg, "xr_cutoff", 0.0);
 	cfg_add_int(cfg, "max_steps", 100);
 	cfg_add_int(cfg, "multistep_steps", 1);
@@ -170,7 +171,10 @@ static struct cfg *make_cfg(void)
 
 	cfg_add_int(cfg, "ligand", -100);
     cfg_add_bool(cfg, "enable_pairwise", false);
+    cfg_add_bool(cfg, "pairwise_ligand_only", false);
     cfg_add_bool(cfg, "print_pbc", false);
+    cfg_add_bool(cfg, "print_geometry", true);
+    cfg_add_bool(cfg, "print_pairwise_geometry", true);
     cfg_add_bool(cfg, "symmetry", false);
 
     cfg_add_int(cfg, "special_fragment", -100);
@@ -363,9 +367,11 @@ static struct efp *create_efp(const struct cfg *cfg, const struct sys *sys)
 		.enable_elpot = cfg_get_bool(cfg, "enable_elpot"),
 		.enable_cutoff = cfg_get_bool(cfg, "enable_cutoff"),
 		.swf_cutoff = cfg_get_double(cfg, "swf_cutoff"),
+		.ligand_cutoff = cfg_get_double(cfg, "ligand_cutoff"),
 		.xr_cutoff = cfg_get_double(cfg, "xr_cutoff"),
         .enable_pairwise = cfg_get_bool(cfg, "enable_pairwise"), 
         .ligand = cfg_get_int(cfg, "ligand"),
+        .pairwise_ligand_only = cfg_get_bool(cfg, "pairwise_ligand_only"),
         .special_fragment = cfg_get_int(cfg, "special_fragment"),
         .symmetry = cfg_get_bool(cfg, "symmetry"),
         .symm_frag = cfg_get_enum(cfg, "symm_frag"),
@@ -408,11 +414,15 @@ static struct efp *create_efp(const struct cfg *cfg, const struct sys *sys)
 
     check_fail(efp_set_opts(efp, &opts));
 
+	double t0 = efpmd_wall_time();
+
 	if (cfg_get_bool(cfg, "single_params_file"))
 		check_fail(efp_add_potential(efp, cfg_get_string(cfg, "efp_params_file")));
 	else
 		add_potentials(efp, cfg, sys);
 
+	double t1 = efpmd_wall_time();
+
 	for (size_t i = 0; i < sys->n_frags; i++)
 	    check_fail(efp_add_fragment(efp, sys->frags[i].name));
 
@@ -430,6 +440,12 @@ static struct efp *create_efp(const struct cfg *cfg, const struct sys *sys)
 	for (size_t i = 0; i < sys->n_frags; i++)
         check_fail(efp_set_frag_coordinates(efp, i, coord_type, sys->frags[i].coord));
 
+	if (cfg_get_int(cfg, "print") > 0) {
+		msg(" TIMING reading fragment parameters %10.2f s\n", t1 - t0);
+		msg(" TIMING creating and placing %zu fragments %10.2f s\n",
+		    sys->n_frags, efpmd_wall_time() - t1);
+	}
+
 	/*
 	// LVS: need to coordinate this function with update_fragment() in efp.c
 	// copying atomic coordinates of fragments
@@ -592,6 +608,8 @@ static void convert_units(struct cfg *cfg, struct sys *sys)
 		cfg_get_double(cfg, "pressure") * BAR_TO_AU);
 	cfg_set_double(cfg, "swf_cutoff",
 		cfg_get_double(cfg, "swf_cutoff") / BOHR_RADIUS);
+	cfg_set_double(cfg, "ligand_cutoff",
+		cfg_get_double(cfg, "ligand_cutoff") / BOHR_RADIUS);
     cfg_set_double(cfg, "xr_cutoff",
             cfg_get_double(cfg, "xr_cutoff") / BOHR_RADIUS);
     cfg_set_double(cfg, "num_step_dist",

efpmd/src/sp.c

@@ -32,7 +32,13 @@ void sim_sp(struct state *state)
 {
 	msg("SINGLE POINT ENERGY JOB\n\n\n");
 
-	print_geometry(state->efp);
+	if (cfg_get_bool(state->cfg, "print_geometry")) {
+		double t0 = efpmd_wall_time();
+
+		print_geometry(state->efp);
+		if (cfg_get_int(state->cfg, "print") > 0)
+			msg(" TIMING printing geometry %10.2f s\n", efpmd_wall_time() - t0);
+	}
 	if (cfg_get_bool(state->cfg, "print_pbc") && cfg_get_bool(state->cfg, "enable_pbc") ) {
 		if (cfg_get_int(state->cfg, "ligand") > -1)
 	    	print_geometry_pbc(state->efp,cfg_get_int(state->cfg, "ligand"));
@@ -42,7 +48,11 @@ void sim_sp(struct state *state)
 	compute_energy(state, false);
 
     if (cfg_get_bool(state->cfg, "enable_pairwise") && cfg_get_int(state->cfg, "ligand") > -1){
+        double t0 = efpmd_wall_time();
+
         print_pair_energy(state);
+        if (cfg_get_int(state->cfg, "print") > 0)
+            msg(" TIMING printing pairwise energies %10.2f s\n", efpmd_wall_time() - t0);
     }
 	print_energy(state);
 

efpmd/src/timer.c

@@ -0,0 +1,37 @@
+/*
+ * Wall clock time for the timing printouts. The library is compiled as strict
+ * C99, where the POSIX clock_gettime() is not declared unless requested, so
+ * this file requests it before any header (and uses omp_get_wtime() with
+ * OpenMP, QueryPerformanceCounter() on Windows).
+ */
+
+#if !defined(_OPENMP) && !defined(_WIN32) && !defined(_POSIX_C_SOURCE)
+#define _POSIX_C_SOURCE 200809L
+#endif
+
+#if defined(_OPENMP)
+#include <omp.h>
+#elif defined(_WIN32)
+#include <windows.h>
+#else
+#include <time.h>
+#endif
+
+double
+efpmd_wall_time(void)
+{
+#if defined(_OPENMP)
+	return omp_get_wtime();
+#elif defined(_WIN32)
+	LARGE_INTEGER freq, count;
+
+	QueryPerformanceFrequency(&freq);
+	QueryPerformanceCounter(&count);
+	return (double)count.QuadPart / (double)freq.QuadPart;
+#else
+	struct timespec ts;
+
+	clock_gettime(CLOCK_MONOTONIC, &ts);
+	return (double)ts.tv_sec + 1.0e-9 * (double)ts.tv_nsec;
+#endif
+}

src/disp.c

@@ -102,7 +102,8 @@ disp_tt(struct efp *efp, size_t fr_i_idx, size_t fr_j_idx, size_t pt_i_idx,
 		    CVEC(pt_i->x), &force, NULL);
 		efp_sub_force(efp->grad + fr_j_idx, CVEC(fr_j->x),
 		    CVEC(pt_j->x), &force, NULL);
-		efp_add_stress(&swf->dr, &force, &efp->stress);
+		if (efp->do_gradient) /* the stress tensor is available only with gradients */
+			efp_add_thread_stress(efp, &swf->dr, &force);
 	}
 	return energy;
 }
@@ -175,7 +176,8 @@ disp_overlap(struct efp *efp, size_t fr_i_idx, size_t fr_j_idx,
 		six_atomic_add_abc(efp->grad + fr_i_idx, &torque_i);
 		six_atomic_sub_xyz(efp->grad + fr_j_idx, &force);
 		six_atomic_sub_abc(efp->grad + fr_j_idx, &torque_j);
-		efp_add_stress(&swf->dr, &force, &efp->stress);
+		if (efp->do_gradient) /* the stress tensor is available only with gradients */
+			efp_add_thread_stress(efp, &swf->dr, &force);
 	}
 	return energy;
 }
@@ -218,7 +220,8 @@ disp_off(struct efp *efp, size_t fr_i_idx, size_t fr_j_idx, size_t pt_i_idx,
 		    CVEC(pt_i->x), &force, NULL);
 		efp_sub_force(efp->grad + fr_j_idx, CVEC(fr_j->x),
 		    CVEC(pt_j->x), &force, NULL);
-		efp_add_stress(&swf->dr, &force, &efp->stress);
+		if (efp->do_gradient) /* the stress tensor is available only with gradients */
+			efp_add_thread_stress(efp, &swf->dr, &force);
 	}
 	return energy;
 }
@@ -290,10 +293,16 @@ efp_frag_frag_disp(struct efp *efp, size_t frag_i, size_t frag_j,
         size_t n_disp_j = fr_j->n_dynamic_polarizable_pts;
 
 
+        /* s and ds (LMO overlaps, one LMO per dynamic polarizable point) exist
+         * only when both fragments have LMOs; they are used only for overlap
+         * damping, which requires n_lmo == n_dynamic_polarizable_pts */
+        six_t ds_zero = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
+
         for (size_t ii = 0, idx = 0; ii < n_disp_i; ii++)
             for (size_t jj = 0; jj < n_disp_j; jj++, idx++)
                 energy += point_point_disp(efp, frag_i, frag_j, ii, jj,
-                                           s[idx], ds[idx], &swf);
+                                           s ? s[idx] : 0.0,
+                                           ds ? ds[idx] : ds_zero, &swf);
 
         vec_t force = {
                 swf.dswf.x * energy,
@@ -303,7 +312,8 @@ efp_frag_frag_disp(struct efp *efp, size_t frag_i, size_t frag_j,
 
         six_atomic_add_xyz(efp->grad + frag_i, &force);
         six_atomic_sub_xyz(efp->grad + frag_j, &force);
-        efp_add_stress(&swf.dr, &force, &efp->stress);
+        if (efp->do_gradient) /* the stress tensor is available only with gradients */
+        	efp_add_thread_stress(efp, &swf.dr, &force);
         return energy * swf.swf;
     }
 }