Willard Chandler Calculation

[octupole]

Hi,

I am working on some interdigitated membranes formed by C18:0 and C18:1 glycolipids. Since these systems often exhibit very corrugated surfaces, I started using the Willard–Chandler approach to compute the surface of each leaflet in order to analyze properties such as surface corrugation, skewness, and bending rigidity (for which experimental values exist). This is why I began using pytim.

However, my systems are relatively large (100–120k atoms), and pytim.WillardChandler proved to be quite slow. To address this, I replaced the neighbour-list construction and gaussian_kde_pbc.py with a C++17/pybind11 implementation, adding simple OpenMP parallelism where possible. The new code is about 9× faster than the original. Results from the accelerated version are nearly identical to SciPy’s, though I do observe very small differences in the computed Willard–Chandler surfaces. I suspect these stem from algorithmic details in the neighbour search, but I cannot be completely certain.

I also tested a few branches from a-ws-m, but did not find measurable performance improvements compared to the original pytim. In a previous issue he was talking about JAX, but I couldn't find this extension in his code. Keep in mind that part of the C++ implementation was generated with Codex, but I carefully verified the code and the results for consistency. If you could review the changes and perhaps test them against your own benchmarks, that would help confirm nothing fundamental is wrong. I find that AI is very usefull to understand a code you haven't written, although its hallucinations are sometime quite unsettling.

Here is a summary of what has been changed:

• Replaced the Python/SciPy neighbour search with a C++17/pybind11 module (pytim/_wc_kde.cpp) that bins grid points into periodic cells and accumulates the Gaussian kernel in an OpenMP-parallel loop.
• Updated pytim/gaussian_kde_pbc.py so that it uses the new neighbour-search path automatically, falling back to the SciPy/KD-tree routine if the extension is not available.
• Added a regression benchmark (benchmarks/willard_chandler_kde.py), integrated it into CI to enforce a minimum speed-up, and updated documentation (docs/source/quick.rst) and the changelog.

And how it was done

• The extension is built with C++17 + OpenMP flags (via setup.py), releasing the GIL inside the loops; atomic updates ensure thread safety.
• The neighbour cutoff remains at 2.5·σ, so behavior "should" be unchanged, but faster.
• Documentation now highlights the new build requirements (pybind11, OpenMP) and describes how to run the benchmark.

Best regards,

Massimo MARCHI

[Marcello-Sega]

Dear Massimo,

Thanks a lot for your message. I had a look at your repo, the patch looks very neat, is well documented and also comes with a benchmark script, amazing!

Is the 9x improvement obtained with one thread or, if not, how many? — I'd be curious to see the hard scaling with the number of threads.

I'm happy to put your patch in a new branch here, but I'd like to run a couple of test myself before merging into master. I'm not concerned about the small differences — the order in which the gaussians are added is surely different. Mainly, I'd like to see how performances are impacted by various ratio of node vs atom densities and box occupancy. It might take some time, as I'm quite busy at the moment, but I promise to do my best.

BTW, I'm not familiar with pybind, and in general I'd tend to keep as few requirements as possible — do you think cython would be an option? On the other hand, it looks extremely effective and easy to use, I might actually consider dropping the cython layer for pybind, if it does not degrade perfomances. Do you have any idea about this? Is it pybind widespread?

In general, before importing your branch here, could you please:

• remove .gitignore from the tracked files, you don't want to overwrite othere's settings
• move the pbc search: since it might prove useful for other parts of the code, we could expose it somewhere else and import it within the WC kernel.

On a separate note: I'm quite intrigued by the glycolipid topic, as I did my PhD on GM3 membranes, it would be nice to follow what you're doing!

[octupole]

Hi,

I will answer your questions:

i) Multithreading gives only about a 1.5× speedup, and I haven’t seen the code using more than 300% of the machine’s capacity. The biggest improvement comes from calculating the neighbor list in C++.

ii) Please go ahead with your tests; they will tell us whether the code is fine or not.

iii) Concerning pybind11: I first used it a few years ago to interface some C++ code computing Voronoi polyhedra with MDAnalysis. Before that, I had everything inside C++ code, but it was rather cumbersome to access the molecular topology (proteins or lipids) andI find Python much easier for students. I remember trying Cython, but I couldn’t get it to work conveniently. Pybind11, being modern and header-only, seemed like a natural way to interface Python with C++. Cython is still the most popular option, mainly because a lot of libraries and legacy code use it, so it is downloaded more often than pybind11.

iv) Ok for removing .gitignore. I usually tracked it because people typically do a git clone into a new directory without their own .gitignore. As for moving the PBC search, I don’t fully understand what you mean. As far as I can tell, it is used only for the Willard–Chandler calculation. I could create a src directory and place it there— is that what you mean?

I see you worked with Vallauri (whom I know by name) and Melchionna. I interacted with the latter when he was a student of the great Ciccotti back in the ’90s. Working on glycolipids wasn’t originally my project: a colleague of mine had obtained funding from the French Ministry of Research for a collaboration with an experimental group doing SAXS and other structural studies on C18:0 and C18:1. He passed away two years ago, and while he was ill he asked me to take over the project until he recovered—sadly, he never did. Still, these systems are quite interesting. The saturated system is extremely corrugated and rigid: it self-aggregates after almost 2 microseconds of simulation, whereas C18:1 forms membranes in only a few hundred nanoseconds.Also, I computed the SAXS of these two systems from simulation, which matches extremely well with the experimental SAXS. Interestingly, the thickness derived from the Willard–Chandler surfaces on both sides also agrees very closely with the values estimated from SAXS, while the thickness computed from Voronoi was 2-3 Angstroms off.

I might also run some simulations on sphingolipids. The person I am currently collaborating with (Nikki Baccile) has done quite a few experiments on these molecules. Their main interest is in the interaction between DOPC and C18:1, though that is too complicated to study fully at the atomic level.
Best regards
Massimo