From 021f460dd1058c86182fb770009eee921d681697 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 17 Feb 2026 14:40:32 +0000 Subject: [PATCH 01/36] Add 2026.1.0 placeholder to CHANGELOG. --- CHANGELOG.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 6a4b6be..15c2eaa 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,6 +1,11 @@ Changelog ========= +[2026.1.0](https://github.com/openbiosim/loch/compare/2025.2.0...2026.1.0) - ******** +------------------------------------------------------------------------------------- + +* Please add an item to this CHANGELOG for any new features or bug fixes when creating a PR. + [2025.2.0](https://github.com/openbiosim/loch/compare/2025.1.0...2025.2.0) - Feb 2026 ------------------------------------------------------------------------------------- From d54fa78e9cef74c875d3e6f171e979d654472a87 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 4 Mar 2026 19:06:52 +0000 Subject: [PATCH 02/36] Update OpenCL instructions. [ci skip] --- README.md | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 226d97b..ac70f87 100644 --- a/README.md +++ b/README.md @@ -51,11 +51,12 @@ pip install -e . > [!Note] > Pixi does not run conda post-link scripts, so the `ocl-icd-system` > symlink needed for OpenCL won't be created automatically. After -> creating the environment, run the following once to fix this: +> creating the environment (or after a pixi update), run the following +> to fix this: > > ```bash > pixi shell -> ln -s /etc/OpenCL/vendors "${CONDA_PREFIX}/etc/OpenCL/vendors/ocl-icd-system" +> ln -sfn /etc/OpenCL/vendors "${CONDA_PREFIX}/etc/OpenCL/vendors/ocl-icd-system" > ``` ### Installing from source (full OpenBioSim development) From 99a261d75403361c717cb43724903b57b9681c64 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 5 Mar 2026 09:29:59 +0000 Subject: [PATCH 03/36] Update pre-commit. [ci skip] --- .pre-commit-config.yaml | 2 +- src/loch/_sampler.py | 2 +- tests/test_energy.py | 30 +++++++++++++++--------------- 3 files changed, 17 insertions(+), 17 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 0043ca4..61e8edf 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -14,7 +14,7 @@ repos: # Python formatting and linting - repo: https://github.com/astral-sh/ruff-pre-commit - rev: v0.8.4 + rev: v0.15.4 hooks: # Run the formatter - id: ruff-format diff --git a/src/loch/_sampler.py b/src/loch/_sampler.py index ba36f95..013bb4e 100644 --- a/src/loch/_sampler.py +++ b/src/loch/_sampler.py @@ -1946,7 +1946,7 @@ def _initialise_gpu_memory(self): charges[i] = q # Rescale and convert units. sigmas[i] = _sr.u(f"{2.0 * half_sigma} nm").to("angstrom") - epsilons[i] = _sr.u(f"{(0.5 * two_sqrt_epsilon)**2} kJ/mol").to( + epsilons[i] = _sr.u(f"{(0.5 * two_sqrt_epsilon) ** 2} kJ/mol").to( "kcal/mol" ) # Store the softening parameter. diff --git a/tests/test_energy.py b/tests/test_energy.py index b2721dc..50e115b 100644 --- a/tests/test_energy.py +++ b/tests/test_energy.py @@ -264,9 +264,9 @@ def test_platform_consistency(fixture, request): relative_diff = abs(cuda_energy - opencl_energy) / max( abs(cuda_energy), abs(opencl_energy), 1.0 ) - assert ( - relative_diff < 0.001 - ), f"Platform energies differ: CUDA={cuda_energy:.6f}, OpenCL={opencl_energy:.6f}, relative_diff={relative_diff:.6f}" + assert relative_diff < 0.001, ( + f"Platform energies differ: CUDA={cuda_energy:.6f}, OpenCL={opencl_energy:.6f}, relative_diff={relative_diff:.6f}" + ) # Reference energy values captured with seed=42 on the original kernel implementation. @@ -351,12 +351,12 @@ def test_energy_regression(fixture, platform, request): # Check against reference values. ref = _REFERENCE_ENERGIES[fixture] - assert math.isclose( - energy_coul, ref["energy_coul"], abs_tol=1e-4 - ), f"Coulomb energy changed: {energy_coul!r} != {ref['energy_coul']!r}" - assert math.isclose( - energy_lj, ref["energy_lj"], abs_tol=1e-4 - ), f"LJ energy changed: {energy_lj!r} != {ref['energy_lj']!r}" + assert math.isclose(energy_coul, ref["energy_coul"], abs_tol=1e-4), ( + f"Coulomb energy changed: {energy_coul!r} != {ref['energy_coul']!r}" + ) + assert math.isclose(energy_lj, ref["energy_lj"], abs_tol=1e-4), ( + f"LJ energy changed: {energy_lj!r} != {ref['energy_lj']!r}" + ) @pytest.mark.skipif( @@ -438,9 +438,9 @@ def _create_and_run(seed): energy2_coul = sampler2._debug["energy_coul"] energy2_lj = sampler2._debug["energy_lj"] - assert math.isclose( - energy1_coul, energy2_coul, abs_tol=1e-4 - ), f"Coulomb energy mismatch: {energy1_coul!r} vs {energy2_coul!r}" - assert math.isclose( - energy1_lj, energy2_lj, abs_tol=1e-4 - ), f"LJ energy mismatch: {energy1_lj!r} vs {energy2_lj!r}" + assert math.isclose(energy1_coul, energy2_coul, abs_tol=1e-4), ( + f"Coulomb energy mismatch: {energy1_coul!r} vs {energy2_coul!r}" + ) + assert math.isclose(energy1_lj, energy2_lj, abs_tol=1e-4), ( + f"LJ energy mismatch: {energy1_lj!r} vs {energy2_lj!r}" + ) From e5c178618aa3748c738d07a6ad17d9136d59fc27 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 9 Mar 2026 09:36:46 +0000 Subject: [PATCH 04/36] Allow GeneralUnit inputs. --- src/loch/_sampler.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/loch/_sampler.py b/src/loch/_sampler.py index 013bb4e..2bd98a0 100644 --- a/src/loch/_sampler.py +++ b/src/loch/_sampler.py @@ -1626,8 +1626,8 @@ def _validate_sire_unit(parameter: str, value: str, unit: _Any) -> _Any: parameter: str The name of the parameter. - value: str - The value to validate. + value: str, sire.units.GeneralUnit + The value or GeneralUnit to validate. unit: str The unit to validate. @@ -1639,8 +1639,10 @@ def _validate_sire_unit(parameter: str, value: str, unit: _Any) -> _Any: The validated unit. """ - if not isinstance(value, str): - raise ValueError(f"'{parameter}' must be of type 'str'") + if not isinstance(value, (str, _sr.units.GeneralUnit)): + raise ValueError( + f"'{parameter}' must be of type 'str' or 'sire.units.GeneralUnit'" + ) try: u = _sr.u(value) From 4156f5b60c168b1b8af14d1528f6666ffe271e1c Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Fri, 13 Mar 2026 16:49:16 +0000 Subject: [PATCH 05/36] Update logo. [ci skip] --- .img/loch.png | Bin 22410 -> 25531 bytes 1 file changed, 0 insertions(+), 0 deletions(-) diff --git a/.img/loch.png b/.img/loch.png index df1fbc230b34156437da1f60296d0cbe93804e16..09622d0b08eafa841dd4d44096d024f88bd3baa6 100644 GIT binary patch literal 25531 zcmc$Eg>xK1ujiWAyyhgPH8VrZ%*+%sGqW8tGh-YxGc(&UbIcS|%*?o#@4b6hb@vaP zrmIJz>Yt=hOEXp79jT-siHh_E2><|~N=u2U001yx000^X2=h-uWN&lz??&w^uH|ZK zAIze`fM%a*#;N``1BKuX2#qKzThQg>16j};qZJc&Ca&o|}Jn{bW$PzwmOpa+7OZ@luIf?=nNJ47SoLo4vZ z-wu+j@M!Ax8e|Po4To4czG!4X>}{XF`Mv*?xDC^KP>O;4bKceb$<+m-5_V$|0P%Hz zKyz4lJ3`!DAr{WBS|+15_Ro%1Pyb#)$f#Y@sa!7#V&e<3_ky5uny5rVl%ns|owjsh zAu+1=Rl;`>W)K7KH@)xg(3JA(X|GCQ58n+Q;DtSf4BOSCp0y%hOhO)mDka5u2$85OkU>K( z38J3*rjhZX;C}&6t>hxVrWEm@^8HpRQY4iR zqEm@fiFi_pdxoKq{pXUKo0HNchNu@DrJQQwi;vb{8}lP%VRxmdM^B|iSaKO? zW_=VPca_K|4WnP0DIfNB*DeOfa1=6*di#h1&MNUQs*%q)qMpu1C(uM<2sBCUlqW68F*yc?nTF z3sXAsws{UVdqEfTMdh>i)V;8=ya_XWbTGehHM|Vax&BWk3%LHvz&Pc{T#0-2Oh~lq zy~YoSdKSdb;IUBr0h-4Vn#*Q(1Q7nOcy3))6jT8p?;ntl|4r}z!+iWl_WobjKklFS zUnBd^Fv9;D0rRha{~!5(_>T|L{~m%RSibb+jt;XxCPkln$@?w0B z>}|%X)xq(1kvFxnlqr{!rIDlCsqxL{!{fW_-f{MSA1E+DT1;5YbM0(nMi$=(Z{(V{ zkWDTOTMm~o`5-SlX!gDeI+Q-d3eQ+?Q)V5R0p9kqnMTH8SO-zLEH!|!2pyApB{Z6` z2pMl?9#Q}q5CdW-9sal++G^3Xt3Ztx++y+{icoNU@^hZf{QQ3Y^MA^Q3k28j%~0dVo3YacT9ZA#>q+F!E_Ggk9!^~Y<&)`u1~-PX+9p|F&PQ# z%Hs7iLR~jSq7C`*pG;8FRL`@HmbwMDZh^|Mn(KSVGCE=Cg5}_(sA)S&eDPJ=Zq_O4 zph`Yp+vSCM{bGJ526m#=!2r|7c3u98UjOHy_Yk^ccPlq;yAUg;_G*EZ3bTQtx$#Kb zl^cz7-}{Ha9JyBSv`TE;f?NowgZY#5jocM26mw z)6D&9J1)_ux0eh@oLi?(uC(1Ujo(3P=XmyBMHij9E;sEpd~;&?Ss7rTc7v9G-sl-rI^N z*~`ST?iL%t{*9;h6dBKQhxU`+M$PHj9AbErmXX767A}p5$YUm(wJIFZ6n07 ze3PkqNlr6$$&tNyWLzyiyRVi>?S8GwHe4S#Lmeb5ztlA(8n5G}BO=12w@yn4IVi3^ z7-S=D*QHIDIrt>*=(YKIJ3%YN`;PYaGvD3GI#sR5`N({xF5{S#3}#XVqvI%E`dsYT zqM21?*X$z?%;?{7 zX({Z*7QD~mj_qObQ-^>1)Kz$kCS;a7b17xMaeSV0NwjKEDfd13b8PcrOPFr?q_q)o zH{_+dFo$w5Sk~}>rFk63*fW#ONQWpTP+JK67ZB&!N_Pq@?~B2&7}xV2{_gXBT@RA9 z_-3l*Fv*~D%I9%P3voxp4c)d%gc{S2qdU#1gT?uGRofEU`lA%kJGJJvq0Wf-qh|-^ zvF^pIUzvkdZI4rj{|K{=&i5=TrEY^{-XGc`&j5oVtk&L z7eiA)KtO0#S0LCe=upMDp^81t{3pv7g)dvfBYAERqk@a7i;h87a_N3D?+?R zU)wGDaKBp}Oi(%5)Ybc==J`j{{L(^j(Fp6FclS@hvXqtOwupsi7pA)`Ob{fC_@utl zfK@GCp{J&}DQv%Z@XO>O!Yvll>T+AaPrgSE@VuVh!h}a9a_>>y8S-cG$ufa5cz8r0 zqb1zjp2nK~7=d$zXy=Gmm-b=+0ScP6)D0{wUULxn?Z~0y-ZD{qx1KjY3E@(pMj{_n zKOAbrq{H5i-lQ(AR}&%g1YzhPWL;sNbV|lMjP8)PKn)I3$#+!SDADR}vzTUbFcA6AcEkZP-ak7x+vjgO; zg{Bu!SX%`lsX3AsTYru&-a5k_N=nArB5!_=kF=h;G+K1H8`|2xpO>aG;WMO{>VNla z74#Lk`SrK=elWrB5tSGAst9TRZm_r1M@txfS|0AXj{P>?f^N`aeeHeoU9dDQEUG%l z4)3fRTqk*ESUby@rYUfC`ST{?Phn1L?Kp&{1vylt&H};Mp%%RUgQnN}tA}nuxIg6N zk#^?oV>zDaAlo5B`-MvuH!fHB4%>SMnw6QEgz3->-iynW+v#jH+fqJ<23*LCt{JxG zRA{uCb};Vkj@)uEnr$~y{V+d1_^gUuws{}VCVH-+TSvbd*G1P>V`q%1RlPeoa=&V_ zIp-4bmW8;gri>X-H1nM&e!E-acbnm9T(2LoEg78n;4R|PS67UCyUH#jZ8Yo8? zo}KKF0_+e@KNenPPD8E$Ifh~7fmoH;-N(bW?2XQ|+e|fM#!<*>=iW%_$ksmi7evKW z*P;wRLNhqE(14R-_Lo>WhOD*vyCtp{-(tyVBec&^OTiRHf=86#Bs%Xus zp9?kKI*oNxxx~=fXq+a|$j$d1+s&J_CnfyRKV2Q%@kbz!WB5k4b`;-P1FpB?c%_GS z7J)l{^T^d#Pk^BTV_{((Qh}wz53GCr)&lr+nDAbM71|wn>9i%_j-9@kJ5poYzITx= z1b_SCy{wMmbyu^2CZ~FE(%4aTX9}N4q~ODFTz@7I8|T;1#_IoS0Od5vL@6k`QRjbY z4xhMtcz1t4JfDG)utyh>{+6j$OBElOb*$4Z)ad z7rjY4R<_?L?GaTHQ(dsQ`)pq?QoaiOmTKj_rl&M^TY7d`>6!9op30YZHHCqS$0ZqAouKk1_F}A7g>jq z9@L0RGH+2LHCmHdc_DV&SCV% z?I~~^k#|CPOlK814pEihE7zpYt};D;#Nyd1it%)C#4TG<_DjC#OF7!@MO72SGW|}N z>g=G<{PIfR&&pX7`J80DwRWx=1vqaygAMN%zt^EyTr<%ZRh%Ye7n9b!Z|li1*rq?{ zpfes@+`~7pn@)T4fK{SWo!bIZ=Zp5BnxX?@)b-6dPd*|yghBRe zlrq3gO-AFV-AmTbKlnN|B+E39YoTZj-)N)-%M#HFwBU{4+Sw7fdE#Rb%!Mm| zTNHq3GhyFvr-!d^fFliGx$*C@Dl@ZyWs58}yF3Kmg03& zEQ6eC4 zkRGS(EKDVHS!<)+`0>CHOB!kRYl=jsJ_KK^?haP+68J)n_8mc?Rcpe-d9IUqKWi5j z4pgZ&d9y#CyY|mU9{!Xv-}AcECw#PX5Fjei{W9o9GpK_l?9t%yb$FaK0Z#g|8wGS6 zSSR$&N~Ov-DUL@Ym7MUY3|Tr?LV00hbq5TYxyRAAyZ!z|1FUcicmaC>^b* zA&m%b&$Ty&qw$mlGbWEYGFqXX;E1?0`z2M|YH)CXzyEth@^xVhNN?{ zR?RS#U34$p-CDd~%sFtQiS*&Aj2dHzHOwOVvYQY42Iw~H)Y7TV+YAe9Nr#)2vm-q3 zFkj`dI0Ncauqn=?@$~DDO1=zq!$H=8RcM|z`A0cEAk;er# zcFUDmb)nP1p6E(@S#u<-sqo@avz|W~J5tMQ?8y=?fZBv ztbPi1+XMF?L${Xq))q4@(&G>SD2=Ne38W zMZu9$VacBrIz%vu%%;^~ls?Qy^{-K*u2NOx!SrkUSw)W}{5P;^Ltv-A($`V!XW2UG*PiP^ssDQOr(w_0Zm_ZHTAE}`x zRDVfyt!bw%`#cN<9ft;oCQg4Ey;Cr`cCJHBLLWT)(p}6Hy~BI=wH+!-#ZRXS7L+OP z1}kJZzI!c!NQj#&;_b*(+-p7;OZMn~UtWk=IyO(Y~P%M*yIgi?~zck3%3Db@kNv zw^1~T4c5DOB@|@CwBi8?>Pf|VkT7S(83K11j!=!{z1BCrAUX5-n*^914k`ng`cMVK zI%f`k>SUyu?Ja>^d0IHx*22P6z`$q2PX)xE`o#vGYUtsqbe_&dYaIJSOE%<4*|lht z62gM_VN6?f@*U~$r@!F^6!FPJwKx{L8U;)uf(76&jc-zNiPAbKHC=E*5vBcyma`>p z2GtgfIe-2oplg4t!xR#t6d}Sk@+xAoeMW(?j|qQ80qo#PT2i19+gC`d{+)^iS68T3 zRcEQ!V|RM(v8vOnd6@uUH4FuL=-j5pW`{?+G-2TsT0BprgmyJ9@}eg?GgZt8r8Cp zq-j9Cf^59F8fQ>jHKy^W$;%S0#1T#pG!i-(v3__m-RSUdFVmCL7{f59m?<{g3q^mv zl4Yc_lfRpxIjegKj+%c56xeDEsZexj-3kpmGem>^)Z&UPSRq6TO3ie{3sBjvXL- zqtm8g=PaqnB(&WdG%i!mlN%R!6b}5I7J(RD`6OpXbei>C)kvx9{7?(0?SfiWX4%Z1 zDW+cfEpt+qjiZ6sRXbp&BCSssJzEuy4WwWH^XDgu*ZpIr$~WBc=_UvFGp2cSZn};_ zj?l{)TOgHe$yZVt3v=dTLjKSb-Xm}nl?lt0x>S*<%}#abPSzo3*$p$Z{`GIquWFy1 zNiB?HwOIs%u57nrkxDxq7O97`B(m6IlitB2RfFxep{!3T1&>?lYEQNrVnCKwSaX3r z6rZpbje0C$j)8tx=aic&tHbyn&L(Yt#?Cx2QDWqmj z1NDyfYO!Ed>hj({!v2cye=dg}9@#Ng|0Q4_rids1s+4$>=qG}MH7w{8@kD^|F z`J$<_t|dI-%&3LkWaZ_~24&dKTTZu8lT|=8V~yiflhTaIKK-f;u8uCkt{wsk z$8GcVFInoZ!p1}dBPVS%U4dxs?WzAli4Lyy_-5--OFsN^JYxG8{M8@Q_P9CS0XeO` z?9B`|Y`{^3P&ZGUn})%85+V#;6#vv|?4+EA6DWti5wS*r?}PJC$CO_|m{r2~`?7N4EclM5LLo#Z|A zWn+n?%#76gYm@9rmw9yXI(7Nh<1&4zC2*vNEN#t#HnN7*@q_CK5 zpVlNv|6;s2m&rVyAf*sovD0!GkNnJ z{G2+j0uR{KdtT-c?zPz%VUoCNMe`I*SIG6Y6e0VrRrm8dEeYfZZ#mK3y5FY;?y3UF zgNpULW{18A!45{Q?~JlY>A*v|VsmxtUus)yPb_VF93Zq3{q3B5Rsn#XTrMBV1^XcG zi-VfayBn2fVmxfnM7$&#_Zi=e_W4N-_fiVJ$)Co?n0SL!(NYu?>uN86lE+EEnJFS> zHrPr5!28NFHA&;K#UixQJ0vf)t_-Z>V}J_a8fV%FCPr>bzxLT~ZqdlECpw*WsY9sP zpw035rn@BfzBGx{VEZ!r>B=*#=X>4zLIFI6+D1Iuw=esXzO()EdkLtdaq!WU4#O~g ze7{knQA2}v?V};@P;SufMSG8=SHI^o`Gw?CgY)w^;QMdw!>q)GN1@0$gg&kduHyqi zhqrjWqBO4ez69<&1E5nrckCrhA;WXK?^8ObzWUoWcWU0Wav`1%n*L9XfQY`>J=;vp zZM}IMze(B3K-8vho|@|T-nU4STV0VSenh-uUV=b0SDU(DvL7ed01{%=vm8x&C|HxK zh+r9f7l3qos)Kx>j39EnAPB5q`uH~)p%jl8q}M(TBg8DOb?Ur04-MZh$a_kmt>=z8 zgs@BImwM@M&_ZZz^@Z5qaH6Zjo1q9X#)BM#c9fv_sm>HN)cq&Zd<<7Mg7anA$Ib&OG# z!VFe&#&T*Kotd2K_fnv4H^=cz+WhrdmdDNJ1mm_&RYbpnFva>tpUwRGI4nSrIYUKA z=^R@9bcL6j{aaut^{I640x|LzV!pHkz+)b2|8WCU>${B52^x_O6cx zd=Pa*C17ffXJaa<`hD=A_%jfs0tO@wt`|~B(c&4S%!kBxR}*sz{2`1$7eWA&BNE&# zyIh@%rkVxJhR}Q&C|}KD4n=4wSPcsQYEa(antk#9k!hHs>oI14fL>3~f3a6)lW?7} z?0A@n-_$3?qrLb!P9Qizj#%rK$_8Xg^7g4RQiS$z>^{ks3t$Y};(k&5Xe`wIRZ^THzFgW_|<6Dtsg*W`A9v8F}R#ncULG^%%eLdzXOnP9PGv@KZ8(@14~X-Q-N;D?n`#p>VN9gddgJ=JJZSj-lOAz<-6^UwXpUaw4) zOggay5V_Fs?4q&6@B#F5&mQrF!xUwXlLnQEi6T=Gn?H9B?yam?{Ptr$8;$sx0z9+r zH=&m<+;R<49wOA|eP%wDb)g=YtQ{2J?7o%$LLSunseg`C2o?scQe!lph0Q(jx`09a zbnx4@J9fIgrj%GuoB$jSoA7qqiUn?YGo1+l8U!T9PHHC}EapeaV3D|nsDGRXL*!2H zN~oYnix09jihwK;eNsgx91a`?v3F`qE5P`ZC%2Lg9)Zytnm!=(U}O*qRXt5qty$Ik zO1;gZ?>8Y$({s%#)6N^sidnsTH-b z1b9i9KBlLjeAhQGXR(2x=`k4}Dyke_f)rqnQBOlX_h)5O3^k~vLc^v$=$GrmkIv2m zbq$yIz2CLa0!~O1jY2VRGi!;L!Bp?!hr57X9$zje$ftMqo|ETQ*C8-$QqE(!Xdz=n z4ej0)6mV4#6aS;N65-&wvgz(JMkU01o`Y|Az!u24)$Y}rb6IxEXcUV+pfFWi!Ch)oWuOb79o(gv1YXXw|)na$b?&PU#Z5-ti9|Dc4iddRv z2Gi_1XX7s>o<@AMr{Af3@?^YvzuzmsG<#F(PADSBQuMj4TtXv2$W7Bm2o0~w{{5c6 z6UDA8W~$uOgJ~CoNB0?-qhmnZ9~Qn0Qcr8W?MO0^mS)YHuD=JXzJC?kW_b(-uweuw zq3L9S3>J?Nl$lDOm(xx(|GM=S38hpqh-^YimW&?L0pvGQK6`+F_S*G&LyZVMWgcM=x3}$ zv<9WEf1br&M4l>Mv^Dd`i0*s^a(evhn(+$Jp+{3I@m)Bm%mrwt%w|V+H}>m6whJj> zAiCC9+<&W9W}pp2ISSMLx#bYJAh;e(!A}sNqF1Dem2g@3XzT1On0-7}#&>K>1_vj! zg{!>;uZ~LZf6fB*ef@MG3}EVU%oyPDUxEmk~ljf$xT2r0L$Ow zj2w+Kwi1WQr-5I>p1S_AKfZY`op{JMM4Djez5mxzI*WXK$Z@B=@T%fS*U$gOIQ7>H zY7j#BC{D=etc0`t;{KgGCLisC;t5pw$IJI^@&xC&gy7+?99ExTTQka3$`pi{?Z?_K zk5WJM$o2u|H%XuNHGBX(dJoXGG+QCPfIGDDVsq80*{8{3dnQ^Dz$q79^r2IUig+G?Hd$P$BjI zkm}!CE~m~nZSp^vFBUrVyiLv%)R7$zR^r9$ z!BDKH2ZX>!T5{i>&^MBX!i8Y)-Ma-5WM=V8tGV)jkIkUoeN$Fr=Qw!r+?07JGW|Yc z3K~GDI@i)DdO1Cr^}pPR6a_}b|CT8;w`coB?7R8ItBz{Yn4M{aq;bOXvwtZLm&}^Q zgy~n3C_Vi*+c})3n~-nZCE~O52uvcHi;fT^zWi7Z?63KDVHy7 zh+iIN6P!OxvoimQ3SEAeG}H!k&m{yV)%k*_@uW?M2||hm9)Ag26$O41bA1{jMZ?Lq zJvFDtt~ba=ENd7&ZMb3u3cfHUQjfs`N%6~uH?y%Dj3VHO<4KZDtBs;_%Xql&v z)3meGPbEhZHSry!#tj0;5}T5M5>lLk6GmH_CV(cZ9Yrc8l{q>RI+c4U#<0z_)^jpc zOU1>!TCFHLv}zOc=S#;a3Z*Tg`6yMq1~2{7UqEg%+%vHyoNF1ng|RbZ8LVm!yXWr* z{vU_jkK7ArbkGfQqoq(yP1~=B3t#owS6rkIg0~nUNFqdHD_(4A);f)oKO3>5YrKRl zj0@kO8crwjvx;?Jszw~hnc_5A-}cun)o@=nM!M+eY~25x_HlYL>%)>IiN7ZAlC1h3 z(F;pBhUP#8C1V27I+b6?f%Hd#A7(vsOv)qY@Tcp$O zfNT>6y+v~LO2;9q+y}hs=WHLxd{8nvC7TR3J-A@BJ@&igKZpFNaTmj`gbqBnhNk2< zPpd5oHxyUk$b`izP%2MInr=brR9bVZEOUi31yES4J2QW666HLABS*jLs{;@Yg#$*< z0~{*PJG{OnBw)VP7LhZYPc1BZEky<%L1(c6iSOuFqiC+uYJriKJBlx|Kd{y>1O;h3 zxL)K*ybh`R{p;c#?#9%vhleV!_9NNwxrv*}I6}4k1A~9DF6N1n1^FOdY;_|C44LAX z=63BKp5S4_TC9^^1gh5;ARh$-{Gzqdq#Tvi{|;EE4GbmGr;aa0Zc=Hjsn5-GUl}zD zmn`byZhVhW7wpEyCcgxc7zb)$Bpln|w>MWn*05^5>!d2H9+83PZu+2aeb5<(E@=hE ze>!5IUI|5kbKBK$qcwU5c$GE{cviAWmIbkl3IxouVLRu>?h~nO_saU7(sf2a*3M8U~XikdM{{tI<OYlmf$Ek2cl-ZjSvCzPbzAv1lc3YdlHGR#v^OS2%qB z40@4iSdJCBl$Z&HIX^hnofa+RfRSyi5)AplU8A=zfk%&v_B+NMp*ecy#_o@aHqU+mbbohvOqx;?71R^zjNw_#kOiZ-i_I!kBj7Z(5y2Yw{ z**%`7%^)-&r{m(`S>YcI7oP=YJ)B&u<+VD{44C}H4-hr(9JhC$9bnRt2xqK$=7w}X zel_Zy9z{4Ezk(t2q&=$a@Gkpda4amTV9;nRE@sO)US1a35oS@%xmkCiYV`u4|V)aY{r6b zdXKZGwW}=O)+8mvI(lE7yKSC4xx1GSo(%O0pEkE|xq2r>BYf1&rY5*?i_|P9yF*;l z)FZEIvnfAXE1;q}(^$l!13HFT-55kP%~c~d7m*h>S61t7f_aJ7O-=RAn|Kz3ELi3;EbDe9{DU8AZVi1cw_4h2SHJLz5@N z43G@~garuxX|lat$$FiD?+@JnkJU%6H7(q?ly#6ik7dMSKh#crY+=DvF`E0hT?CSG zK7L7CknLcd$7bzijxQ9;^V!^cO#z_~)@y&53V~l9(0bgiEWj{TvrE59ZY!EX0sB!f z+nbA1;7m?RaoN(fW9PcQek3`*!~IWrS*Pn((=vG>V(9LsVkZ_E?g}=a_ zD-RZ?D^IKhTPaHhD*RxOFdHA|eVL_24$@xz60%IDbeK!A<+z=RmSP|HdIwaAN0$sk zDstERv~#E^;KZ|s<(lORJJctFoM1UZtnhG|=c25xQK2?Y9u)8gMX`=_*a$X6h2H5u zvl0mUaJIsS*eZE~c9+^7!4E-@gP_=L+g4kW_FG>%Q0=GM<`1Rmm#PowfF`rmkHU9LI z0R-S9c%Aw}5W+2p0f-#AlL5d!2M7t>fplVaMjpxD)&4XQWJypbN}IA~3IXko?N!=z z;#UXZfTX1QoxBL-A=D>1(#||sfeD^>hkfDNCY;GZ^q+hI(xJ82Ks*fkbJ~7CLvxD4 z=N7Fa4n$RS@Cp-Fx8hB`11sGT^1aVB=c zn4=kBPn}Ash7KGKHUC1lUGnx{p$~4NML|WC+HT^r1kHge0)7_D z)D1NB+-l3J2JBZm(6t{3P3r2YDdAhrR*s`dMmqMv|l-E1Be zkU}3D#D-Ik+q5s-1WhbriqadwEaIrdP^>Y8z?h%yYttdlb&w_@yHf{~y>2eM6@!ZN z9cM_SAIJ3cIf4u>A4x;R8(Q6yBdeWTV?Sq*L!>t7D7~=?NLKV|a-<9PkWbc>C{x;s zbu|m$Sp!iHD;8GSCtO-u)!x>a>akCFOiv6r~}o1~NVe7GqJ zN$aZVQ6zqQdZ^&nX`6O)6c9>CrRd^sn0y22uKilyKfU!C&;q+RSjo)2#-Gz93b9Zk zJ_8fk#84)P#6X}BaxxQoQryMpIAQl$kyKp1o%GG zLG}a@a5VMu@NcgU@_R5Nhng>{%7RfX64AawC(+AGOD9HgWLH~ z;=Z#%qIp5-DG~Bo5eQ#gb;z;P=1AI;@P)*ioS%`1z0-`Kz<}|+`0sgl$2}^Mn0*-t zB4r(?P6mQ49|!h)OKau7#*WIsG*S_FXXNZU5v8G-2?@V8)=sImik#pOS^J|z5t0-uL*#_(e$lbhW{m}p zmTzu?@pRI(#r3TZqH!C2oXZ9PnhpH_bUO+^acY2Y=Z279Rq z2?=OLh2rDa#fTsEb)(7&=cP_V%~=A$j1c*%d!1|pGy-9>^dSuB0&AHS@iv5Y<{2KY zS;7PZhf_4?brr_vHA>BkXi_98`KhHLEzDJQ>;#|fw0#p|{}^8@p=h)K!(^tmuQ!}EL5NX$u&_`RPpkl;7i{5ulF-Oe^df-C4M!fZFNHNj~q4Gu9MZ< zb$YbfI^r&w`u)x;_s*J{3xW5706xQJ zO>_P;6NM&Pc|_)9EZ|1`!9n~FQ4})5u>xoji~2GIVL~Xx>#e4+j#sBg2Dl0d@*wRr z+GCID3F2;FI72XLJGtvvy6+g|;G;=ha|Fr=cb^0S*4WaXoiXlI(ia#S5&u3#(dm1( z7mvyamJ@Dpj)0jP;Bd7AL`Mmq;RMC?rLHgd{xnv`g5HZNBlb?GWA;n-bF>l6iU(pA zKnpQi{TSIZE1@P1!2*PMCT@~@)doS!5)OdJt# z5G(WX2tKaui{_3J|2o`*r#9Mrj=f40N5C<{!G#KG9>IIlI^*XX5~=+oj3UR&ek82R zn7&oqD=e@%E!&w0Oa-Z%f`>h(Hyu8H8%4S0O$^sWFt*w;JMEqop01-Iz`=eZslR`C zA6VxvyLE=E*P<_3tgITHME|?@bB4@9Gje#o)loEq2B?Rl8B!?`D`DSZ*H6l}kvZU_ zsw`9nmmhFs)n-0oKOV1lUxj$BL`^D*A=QXIkW;2cFMsL+8$t_YU(P{9A4u=obpOrG zMFIws=QI?4L5T#G*Jsd|$i~0cSv1!t-V1ypf5kGniBx1x!v_-oE<9v`qfsHooC1Zd zOr&=yjqtDsRS;kbp`4vCc!My&Vq$3yztB-T-Pz_`0bh_pWt>iD5RO+|%xwIwWu)p| zrKF65uPzjMVNu@P9iOj$-%i@tmq}pbCet)=c?IUy#R*xp*4MGnlDuZ$erEV*mOczZ z361$n;N>#G_}|7lyYO$Q9R~p$q1aS5mPq(@kBb}bxsFICRS%fY~!P-iL^A+gB(9oF1?8L07k(@}8A=nX6 z!Xo6CA{+6Yk){buXi{RS4~zWqAKn9Eur3Uxs*fk)x2bvyJ2+Zf?zaVA%NcXvvz5pX%(#d}FMcv~4)`OjV&@J`=|zC~Ql0`_3mLziNU& zlXVxYP;;2*1)TU6CTvSyg8qQOb0Q8$AwyhbilI;>(gSzz#YV~kf=mqc5N%iSnag+{$QZxT-{Fh|7oR0D!?Zo=XMHNf} zbnw2vv(?woiv7PVF#FSbM_ABcaqC^-lzukIfuuuEqISXpH5gu?iaMqKt|m`x@$EQx zxycyLv(GY8Wbu1-ZFzGZ&e(xxKm;GR2Z(420a;>nN!>$6r-unmIL>UBRK=*`O5W2K z2d5aCq)*b)MOfhP*vQH@YyvA3sK1ELzf3|q_}3?N5qaQE{v4i|e1}?gC@1xK(0y6R zm|l`exxF65=f~1tJOx4y68b2Q%Walmc!Qn=+lhT#C85Vr+o|51;^pheR2Qgk+eYvcgTs3Hld7+Kb^^ z3SN$#BEtF)9vLYoCp#xQJ3GlLocf1awboWiIM^#K;#h|a$Bjm7%%EK;<$gXXqx>1N z?D&*j7~NE8Awv^_o5bW$#hTCEf}C&O?*s{|`hOWfkBFH;2?VmvjRhAYo0($l3VAkY z>tZ-mt$qUUsrM)p$VB*lDRkVRc?s{7eJT{GBd zOfZ8!YvA@|n%w1@klPAZYuNB{{QQf5t6iM#1wvg0QCqaXCA~x=PaQ~|%=T;O47;=%!==i7XT&(YSI*PO`W3~KVN*rdR!Q}$2`W8|AastuGRd%5y5Uku6 z@-17TBQ!tzS_L!5se@Tb)P%{^kIt&OoYHvs4&NGsP;-$J<$-@BqGcI1w>p%nxf10dI|+%{Am7BU-Gyu}fC|BCB0#Xrf%sI+b-;59 zZ^F?HSL{l0r|yIbH0%N`wB$ZH=4a34H{{NYyJZ%Gjw_s3K{IlL6CLpGL-0Of2D!hkQ^3?m$8F7_SJgfOLiW3hjp~>0*Pn{1=HDy!kX)d0rX_JzIm*Il@1F< zlaB;F_g-?6?*AjLW=T&-Y_wKMqy=Gt_1ZW$^D1U6oF~?39c9E5)D{;Q$hUXWgbK(Z~#)Iok77=hzOngkLVRd=~U{N1=O<)ZTToIx9*bP!s7ft@f> zGi^Wl^BN4tLa8TXrhf&_FnTvpq9rcx>6OmBh;!2=q_7n%?c`5fDi zKlia`YltCWvQ?FNJW~|ik`F`ksYQ(=Bt6S9RC6W6MHDHmrS83mS_aoi@YuVReovrN zGkn^kr`KS&e1cYqEsLNi)~-sT^mb8_2`~kKM@KH2^@p}3H1*pn zz6%_MFBUc-l6JWBGE@^T<$;gpis6Cg+IAE-Bq`2v4cc3*M%dAA^HVZ|si|c*#Q~4l z2PQl`^WTIBB8pn*NWlqa(nz2!WrxZn^!$`P8*R~GuHm58Uod41q>gDH_%mhn52Wbonr%;o_}!ilpuOsGa{p z3WpwKMvqvcrBI3z(PYAgH0l7GK<~MRt3uzXViobYd$uCsxP<&buom)>_nXLE7%q_r zJZ=1gD-kM@3v7e|Zi)*14g!$dSK2WV8@zl$nVg0XGB~8roVN&|{+QW%4i(F+b?V}M z>B_Q27j36Pcz<`7f-60?)9YYT-An*YJQ91-^MM_dZ zx}>{X>hXQ=e|Y!DyVhO%+;i65>z;M*-eC{p5~jCL1X@WU_P)XeOuthS)Ql2X z@1OfG7`=f;N#yvhHMDbARTS$j49g9lUp^(r-3)~6Oa$ZVH%z5=9u_`k=GgGT`LiyP zG^ja-jQ*h_U3=9^bm>BP`9KAV#e!V~T|Qg>1bg|H8~J${!I7uxW`el-@6V zm)(EccL0UC^FD#pSX`FI7?fvvrau%-X4Q%%WGLx|iPggT*ZnI<5!;u)wzY*iM;OBE z<;|wE)>NRO0pYizvIPVP`vFYRX0&IyQprpNzH3$koz6Nke9Y8 z#zRd%kpOm}%Jeh49L^#FVhXFjs&Ug{sud0M3Bqf%%6-x?X?O^XJsvd*RqVSTiV7}N zA1tKxfa;fsL#)F>UT_b-*o*P}q{Ln#|o#jb!bGa^WjCsP%%WMD~g(!yD62e!)a6Dzl`^Zj*4iPO2t&?*O_u1mjYqT}Q zM;;xGPy=5T1KFGt)2&JinroFL`2JOq*VxuS10qMhN+*nap94e*uth*NCmq^Zb8}5b0+R&u_hsu{vXvZT@M3Hu5@;IjF=o` zBz{;WrwvF}h2etl>#wq|8|p@;&rh1EDrJbeKa|G+U^ngtBJRLHstE>6sp;LDtLyc>r(k0iN-}LaO z=S88}yG35N!#ez{OX6oi%<0<4`t=*Kgt+onG!u^38nDP7PK}8CKCAjoq>hlyG11WJ zh4ySDGtx>hJ#jdUM@Xa#m|vvCa>fXMaIdM@RqCCNz$fp`vR#5>_@(r)Z$|R9?XrfS z1Hcq)W3kp~^)+N20h^PgC3`2a{~3k``=FtQ921#4v7&$g`&Yla69U07_E40;hjvQx zPBU{e{*e{|R#M!1)=YL{Y>0?itY?j(sIBWDrfV73b<|%(_AAyp(xdHvTfOWQ!zOi# z?OT&U^tL-W;14+XRbBtqVi@+BxMLl|MM-< z&H=`ZNlHiHIy-z}dKzyCid%w9+>{AzM}7ipCPYZdAlBK{&H|{}!4EZK%62eT2i$wt zZV!O84l(Qp0ObjfEqerF&DN%d`wtZ=^2ijkg!?>REF$p=XRwZiYW>!Agw?2^m)+l^0T( zVkuX4lI5J4>oQDP8qpt$S`d`O>d-(@dwq}K#9mmHaU1xKgHhpAGX2luYQ=!>|D(0C z_!(5#ODxUVS(cB+rY_H);2qL?o0AP0h#8RxZM3pgzT5d)Ia!N7j~JGLg(D9mDfkuE zDwR+&XA&*N!T~3b&Ga-T)FX$*?Czzzh_gApQTWQj@g?#Yf%$?F+A!Eu@cJJsV>Wzo z<>O=Y@lRwOwQcTl=T~M*D+OT+uOv%+OMQQ z4z`IEfU26Q#%#42=-pJ(SeBg#|&()@^ z{FI-nq#_4pg$7~L5abml$DUO+7%>o{x#fM%Ts2cg`!ou@As@MLA~EI*0HWv+%;ozZ zS#6(&a?l}iVI&tk>){6>u#Z#9X)flZH6@P{MuUa?ZwfwDR76i1LPvf#wOkJ$mf0m$ za|0kY)o#VQuntU_zO?YFXG^KDPz1vYv@YGqGXLOQ!Gn?c*67gz@D}gV4v%Pr+^5da zP#_~sOR=}3*{=>9AVjOb`$wsYlCnoMR~w7pzpQC7w4mrl9FeD28paWH9vm0fR&E`P zV1uY-)8cAYYuehzTt)8c%?&sxX=pTyjEs7eu?9oHB+1G5&wd3yOUm^0E9d%J%SH=` zJ|m5fk3Ze{Mm!yI*zLE{)q;q#}ZC|l3hijILAlWyRV7PZ5 zl<3$G3kpE6P4;+$)o1s22?lYp4#HwJJdcZkY?<5Z8bHbI`sKt-qs?olwtV1WN$E0$ z;oBP=0zfDNq$|=;z3oiIlT$Yp{KSAqjj#vPI`#s#FlJTQ+yy)sz)MF_CBs{SkR8+*6_wo#I%)sLB1O5J^nJj#^ ziQq(C!4%0*L1Hox}dv1 zZj3;vw=YFqgNq5H&39Osz(@x0J^e_r?lCi;D}HHV&T-wiU0~@17U}~)@p`eD8u0jY z>EMIN`BbAY00O^87t6viA*;ngZq~m(Q!0(v^t_GVemvzXh&Q1VU5nuV7HUZKqe%jJ z80CGltt5E=-i|`b=(5PsU)4=0cASvCVS}~bm|(})~6_O{j&~iJ3wsHISW5bj?muHhMUnp38sfpLWGnB5|iVq zv@Pht$R=y-7kU^MG`&aBerx})a%?uh?Ml;!nj2{E1^g5|^qmyn_>BSpAs5f-r=<^3 zYt*#tA_-+d2KJt6Q(Z@#v#O#lr0v|ejM`@^;OMnppakWgN!dg|YhgA5dgr&0~UyD37 zPr6Oh1YE<#Bp-{nI!H7Pc2yF{zM0izP4I&KWXr{oZc&vsiek z(o;b=4h+ishED|kiuk?C4X3r8 z$_1&h^81&Ufz!~swBg-h7V(YvN2~rr$8!5uE$iT)8r%Iw`b%S#^QwKY8Ujwvriu5w;@Hgg^WHQ__8h zWBOZsB%>CKQZIJRyV@^by!PurOhy#Ggn}RxbYM-Hpp)UtV+(l${K2e>9oL9#)E@zp zKO>(VWi`_F!?>hP<{i-Xv`wrQ&rdU^Z_X_HY7C8@%E`FA&MT&W=(f=*H+=Qg%liwS z4B%~MSxi?b*OI=+=X(bRaH37;5530Pm<+$3dH`awu&1g@o%}mtkQm$tXGeewt~}C! zCF=>9Zjb1LzCT~;d&4t}MQ-sXm~lryfBah-Q43?;9&dUT3I)Nk>Nr_)s59u00YxxP zgyA+lTQ|Ep=64|x>1@ptKMf5Ku(x1@dPMxg#5O$=b0JdsCoGjMo`%`~L~6+L?VU1> zrejk@OZRHc%Wwo^X#Pl{l~G~Ue{=R?`oVSwuh3enG`H6j*3$IIvK?SL2<#sw@8aq@ z&ka`wV>oFqAI)$Ps396z(NN=08w9bgu7y7|~ke7+MDUYuZZfVrlvUSBkEx4z;SvioyDQqX)?8rR=b=@A2EK`i`5q|+h%h;i>E z6@1nTF(8=x^NaV?2wFKn2lRZGXIeY&{U`|_v!hE?{Xyq|D!DwwV#iRC2ps61$NL@ z+co5BNRO!ZF!~Z}-wi*!hzjNIdSXAbZZ*UCYj?Pr9!SB@53JOy1S=)t3O(J?}9 zoBB>(%=Jb*eWRKw6It)A_2)YYo{m9^0}f()YXMtH=$94!0_Of zn4a}AzK_-2^@f^-9+4lo*n?Y!lYnYs9PswLkGv&DSr77D889U!P&P>3id{75JRgWh zNytM)>J@Pw=Pel=tg53~2x!D!g$F0fgK)dD^LxxBh6VHBh%c;RWe}={ZwP9HY?8)U zHP^bnk^UzlHo&XZ+6$%zLpC6UY4#xF@Gk-6d2Vj%YCWa+@esm6M=NKT__>^B}hY;``jL}%-gil+v<{a^r*O*iaMR+Ag&b{U&?VyI|iDhQ7H%ld< zKc=UBVvKjrH{zfqp$n%{=hyG1p3ykiS1FqNZeVaUc69;s&8wf1=C^ zMlhUAD77bW_8o2tL1B)?b}VLKSPBGTp6yBa0H%b==}+;4Q5v-H0K%P?0We*35&@^_ zzY9&?KkJ6*2jNs>zMPb?62#xb5Cq1&q-4P_bBembp#adUCtqaR+y6YJm=W3$n=bGP zMqn>|`So19WHRxSK+07JE~*}m^z$f#N9giURmF!XVv8T)-*?x8zgg+$32-iMjndqA6%k%fu;`qb%#IaunruVM}Kd5*qn(3?N6=@fmINx^VFT zC?eFJL|obH4@<~E6a7qnib9Wri;D{q?^4T&^9_vuG(8=ZkziM|f7?Q%?omfaB=qhM z@_wX(G4sqInusIu91*9ae@Rv3-%3lD@L;J#c9&{Rqe208`U0IUnR$)?;2vk7&wBo= z6p|icdlaI?Va*>FV*VU!$J=|?7 zkBdArsPDl;`VdI>i$c`$PJLO)6Kf@d!^?6txqo~Q`^FWxQSRx#Ukrc| zg>m?ZU>E70sEHM7rH>e})gx9*$0OSxUm%0{O=4BD?-tn~=^K?xWRoq2C95GiDWMm1 z<^pB8`G7;Gu)?ep{}8k~%w`P1E;9{-{Z$3eHs#8hL4ifdoFvRI2c;F`90AO{b# zI{N+K}XP!jt?=krG{iSXO$sy<$H03&!fn zT}ku7bxB&{i6maPZYBotDbDY>%(yaEB|_)sMz$NIq>@5WDZ;(}M5sPe$JQg%=klwg z);Qcne#dumB(v@Tez3A~cw^D$8b&`)3Yi)w2qH|v@|o|*-6DbzN7ddapYDZEQh7WS z%>?oB3|v#DW>~8?rg8MGH(Kk9I5D#D>q0bSTQFDsuxJ%ki~^v7>3>8sLU0Vm&50eS9`J&p0zEKvZX_we*OBze0nkNgANAMlKvJ)Augv2-xE4*dz861UjxA?I?Z4%UpV+n@&cVFaOHosu9h33NqOravH4*zqq#D22rBWZl1omGpLDa#ZGcV9?$nBcFL}l8w+fma+txw+ z^hqfGe>1w>it z_0H#1IU#V$2Z0h=o}=>4EY`19eu_2HuZcYncaQ)|{!ele!dVhrGCRQyxj;wzy6AtGs+)S+H zNhmuyzP{NT*cqA3QC;5r{MLIX#KKi#;bbrRdU?0IJKEowl~66L^Y;O+qs`K)WEM6+6+)=70Ga@B=Ny0T}TQzL)r0X!TLz0_?x!p&Esin4nQ%UOZ zI(BSV&smd;N>C4 zqorydmVV-+Q;Wc60Sm=X;uVb6#kmz{_A1A`l$gmedjEjA$fNHbI+vt=hIH(#v^mx7 z&NNA|9>qH~g^AmG1DQT$uVbwp;h^J&tJrX-*gz}W^NKRYf}4#;TdMn;wSwbzj*^2j znTvILmR$dc{dwg#-Ykpse&q9Qljp1b5xpfwtN5f=w<}3AM#OlIUnu^CFkXPZla+!M z-_dzb>#hJ>>+byA+MnQW<*9|dd5~812 zeCNKB?s7GHH0ERlvca8mb1arV>DS(%nXn?S9bbz-`&LK!istfN_KR$Xi>}&_{St38 za%og6ChpJ%G!HyYC!C#U$U)8TRE_BLOxC=K^te=*Uu0A^KA7Wo zRvnmFlxto|aG><5dH?V;)$?(Kel=dl%z&m^lRNSAsgD+cZE}`869%z!r9G@0@_M#9 z5JOXfie>ta;f~VG^=a$2pdxO0Qn@KaY!ZEwlBGF2hvnSB*0}KVF}IzaBW|YRRgp3c z^q(}ZhFx7kgtI$@8^`2SA#y4Lh* zmaYHu^b=!&M*TQ^9`qm{T57i^lRqkt&=z1CcfXx}TsKK>tqhN1ADLbm|LtpG_j%LC zblPRoTC*hq!I})cl;}5HacTJF*Ai*8J0D)>=w|cTv?V&}6E^g}KJAZoGb*K2D&De{ z7I}%PSB99iM##wO*?3)v-%a~j>9@>Bl=ls<>2u2Fq0h6^6Lk#)oXie>`qph`QUzz8 zP(w1w5;29zl+2ycy~rPP?H*%$Gz{-&p{B#xZ`m;-()n~4lu=ONDXS~Xy_ZUgj^Cwg zKCjMHD05aDv%n1d6t1MEn5LZh!k8&mC`O97$%RuLljX~A_+IHnZ&NjCl#iNS2GQ7Z zMtN`E{`^AR*w|?##~Ge5g@3NCr$2CsP)}{GhE=U{k(R6CS~49mG4##x~{###~%GQWYy8hb;4n z{XV&ma$XBhCA0s7SRYsNnyIA3cZuLHwWnup`lYGlClwV7c2HPxZGI!GXkU zt&3tmYp3DT?#jCVY6?!S-X5MCp8s%K^K}t!_)siT{;H2L=F7H!4ay*ABggEq%*aXM zc2B>|=WQOy3hRPkwq_ZS%8ZfB=YhMKUPwXTXSQgoFcFeT%*&dHz6q-_*-Sg3HMPxLolWAFGyQGHovy%~=#}%V z&Bwyf*yiut(R+5@T2u1j6c#aj0A0;#qNRmoKyzN>PdQ%KK_>9hX(P36x$C=GYzE* zQ+U=McR4WsBzaG+b*y`B`SusB=w>|Q?r{j_bo@Wb`P!`U2dv*d-K1{PHgqJ%&V@)i z^dwD9YA=xRR(SVB#^2O3TqL>JB<*Ta4<|fio%|Y%9#e##jtX2eYn5{h7yP&Dk|i>& zxjxrHa0d4k%Omv4H1Ks1Sn*A}Ols-k6N5K|OL^NI8om;ZR}N!yg&xX%q>9X6r;nW8 z?{^BIhU?kXeC>FWK~c&Z)-}jhJEeLQ2}`2AJf*D9Tb0R?t}61%C=|2ay!v91K1WNv zY0}g+7XMAKaP~2zB}5+x`2R5W^X2D3_c(CJ{|h#)_5UzibZ$v}kKdD$0ki Ltz3b68Tx+!x9eRq literal 22410 zcmZ5{XIN8Bv^7ntfDnpQMIoW9fFbmkBGM8R6zQFW5+L;6n{+AC5*-1|J=pYQxR&nYu|&zhM%v(DNl_66)YJuN3K2?+_kriPjx2?-g7goN}N zHTlgQZyoIJ%|)mGD??8`E1&x)cUQOr((b;eAIk2&gN+Rd2_&&Z-z15f@Akh<7MbwF zdKm`tdA+@mu)DVjvojiKL&pMfn{P;VypI$9*W;%;k8g|E>yNom=m5s5%4!Slb3F{_ z)+>Gxn%gsfz1}sv)iE0lrsu_SuGo&LH*yju3|I=yFSFyd5ITd8&N{p$oNT07In zsz<%;BM|Nu?qawi-?Xm`+3-#vn zH#~p!rc+5h^q#AbRN=VsHy0ELWvDU8orvt-+oxAw$% z1kvAx80kn%vLk*pCpdT$i!2E5?1)*`$H{g?9~AMu&B=S~V~o{t5}X)^Ao{ryW9?26 zKEwnBG02%15;2y<6kFmu+Y?)t%RH-NH^*~N z$MY~JqLs^Kj6E?JNklkbMBTI$ern-*nPGF13@0Ypo|IS+g6+@BEeJl2=eZWdO0%O_ z1hLSJ=!^ddStoxgGW8|_HUHz%gq5F-(17z?5=isVxi@+3-Th}nwS75dbttN*2ku< zmyX`VSO=oF!$L|0FupA#_xewuGbjI$=% z_z{Cp*Tt5{>E^^l4`PrXG1r)w>_p6UBnEg8KNt{8Oo^U8L=^HO!S=+5Z5K^Ua5{YtC;HnGjRJ_NE~hC*#P<#-B{s(x#BrrHA=mNbIi&7p zE|8FV>cR9${@a*_laSme(Nt4@RZKy^Tkr-dN^y!u#11VoIzc@cBR}rk6MPHbCbgT4MkbE)3nIrO7)JUyL(7^_) z+aWkvJbGqGW^-X*yYNYW=6++w&;>|6zTY|E{7`9gZa?_USX7iPzQ?-O@c%nav5(#B zI)pqAJ7L{w@|EX}>$zUisA&U-IrXmA8Va%*uUpo%0edt;H2 zTioySS$}^lRn{0Z((hZ>=6i4ASLT0 zCr12HU*gHd-}?fssOPi6eP>W-UGLn+&44+jvufa z-wx_LvF}=K03-??ahWQM-=tz z%}|}BJgykmV|2dSfX5E=K({P0{~<0FpQo2d6veyOKVteduV|=& zJ>MdydS6n`SeU*5dQ|2_K(7bQdQE`!D;mhV-Ry*vCHV%i-3K-1E*VPU8<@8OfhF83 z!AWA;QLd-4&1YZR>U*`O+>@fM?{4HkNhlr-U-U$5^AIpLbwj>+pdTpahR* zbXrOIiq7rr?;|LIuGXT2s=vUc?6v7S9brBv*si~SPg8MOlpAxl=8vx7ry%WT#_ZCz zcQEkiu>d}et$zNleu-TSO^dbAlW_f;YT1nxC&IlH-bt;{XO6v4Tw;7b6 zD|11|Em4T4E#+SOCy5OPVudpKdRSyYoS=h#RNau=W=*S^M$c6A2(8t38v7&(j3|LV z#al+9m;}e)9?Oad3c7;cloQ3wyyvC$o~QWTt*&51tVZ7k4OXsE#yn&H&V|D9LqH#l zLeiFPvopSmvXntvhLR~ca@}R+K-!5Oy=%QRv`gLL^@_1=Mka%`c0wrJySW!8Q65u& zYdOdut{%u6m=YXJzJEf&!hkjXxPA%xGEZz)Y6X7{6dCREYs%dik)0uL^56iiaMnA6nFp|CuyB5B&CjU32pcY^a$6apQOHDA+t?vOV?$_<;j;<8CR}gIy$FvM_`^f}{+6)9h$TCk!ov4f*fe{}t^~ zt|qu6jewJ^19^58?J{A|18nGT9_Q8Rx}nZ&k3a8eAqBG(s|k%?qGTFH>*-&froX<0 z({z?ShnqAcj>~FCc`#@5vOWdjTw99%1Kwcx+L*n}Jo_HKhip2U1d)}SX2%OlhUW?| zaeGp=OLy7l-n7dIqQHxVu(snFEcKqR9$!EjMbnP{Xqv@(Si@ZDRwJBX8F3BbOP&=d z9Jm4Mdmae*-$I+mw!>v$C1Vzf?=Pm=^DH4ov@Kp&zK(>V6u0Q#q0k#hg*hK2tCx{& zMp{d5(ZvJ-08IZkj5KXOhNZz53l<_lUHy71m?K}ur#PA_8vDz;u#Ai7pGEnYs&<*6 zml8M){xNeix67c_-QSB?#`+%@_H>?Kacl3X%AQAz+ktH@(e@ZX%G>JiS2pioG=Hs^ zWW9~U;V2_IDk^e&dJ+dHjg{o(8-DL5^UmzTzFigYTz?(-_dayIWN`Z=M-}k1IAJfo zp{Q4vphOAFEQ*Wq<1#fdF*Y{NV-3IeEBt*z;}86te`H-6ru=hh+MK1=!%v05S2hMp zZzM5C$BSB2w9ku1$|IchV8pj5^Eb+yWwc5@&dz@HUm6-S`VB=qp%iJNk{Z=3efMuW zyc@=fDVO#Fv&-r$@Ox7oR(#-?)6=|n8JWKpSh>#cEfPsCX-@XHBJ>|YzF+p^-`}pM z4{BjF;d94%TVYmzL04x+o4=XNfcA@)!;oO~({YP_>S3SK43eKqW6XRsGG>2W>Y*oG z?F8v>6#6xIPg)52KSt&4TT45;yFn+kORe6{&4QQ^nSj%sMD5Hf^Y6w4C)Q5CZph03 zb*hUyC3IWUbkAPcqIW}qOAgG?z#W4P3^U9Ln8r~aiLEFR$aO{2w~=QQyDB$cTMUSU zcgG~WeycrYLi@aoxwWje+RBI$XXE-C!xg$w^6BZ*r;D$jU}TW~<9Np2IB2V=+Q|k* z^Tu;)4-AIjzwR5nHGNB#z2zIbxF;9!I9(IM{tL;hg%;1n2EQo%(i+xyH0Gy>W{kys zkC-e%DUMEy(tzqNzLstbWp71EQeo5Chlx*LAmdKZilg*+?8y0vn+l-K+~D}p@7@2# z=gg*e_x9RQ%PP-M)-WmvS7KK)`3Civ1zu?4k<{ZFv)&No)3OQAIM~~dP_Z(lRynDO z03X*AK@7jG$7g=bn|n%>;Ne$RSWRDlGe6`5%=*eqY4u{U`A^UQ)_dR&^Kf2X9<78w z)g^qBtj3SjYkMk0By2nm3FmwxY)LTr^pE+`aRlpV!|7QQ!|h58`FoeK-3$7b(n?o` z0tVX7SY$@He$B;vwj}!}W>~zyAY6D5KCmoU#P~pJ-Hz`uRznWF^};2|$(9+dprD&r zG&4ifJM*1RGj|%(>NR-R+z7%&-_@MPBxua^+@Ei-p|Q}ev5^c-{&S?#O=~lIHFMDJkN?UT!$Dg196QzLc-A7L!e<@xV zKWHF*<#Mkf-d+#QfW0$`w#B4Lrb!ag)R!h^Fp`v;)*|6q)z_5>EvdfV-YV0IfR~>B zgqM7)$U1EauSI!Tg~E>JR(ZK^uUK%}WH3uY^&}l=&HFwe5Dv;PWx}!@EcdwLVK$vHEI&J&Vw*dgLx>B8jINwgUq!VIZ`M5_ zT%RO;I6e9MEAi(>Trb>geo$iLT^-7WhnpK&wO3u$b@q7-brWc`JB%W7FMC`+a>3g5 zQQ2Rh^5hWw8N2_!W%?c)+3M|he|>YdvEWDgEyel}ger%x`QG$zd|_HP2VrG?5i zhqMtQI7$f|WkJ672O(3BrZifZ8X<$_gN%%l5}7}E%#;nofj+=C%h;gAnbQ&t^;fnB zUjCZTkl|^Uw|H~5@^3ILtp(%j6L;}&v#Y3nH>MlDKk}^;Wv$o58l#V>zg%tz9O3j@ zMr;0%AC$;CD^>Pkk>f)~bxRaei>c1y?uqn0?E#!@EcAz5Eb^mYxsJ`r&6ZY5k3zzR zE>aNKVdET=0VhW_7*!(2+m1LKaO>nz-@-xQ4rzO=(AB>S_76Yu-%T{R>i{G#;bGFw zI|}$9_Z>2UcA*4-e~XR-Bo@0BGj44JlNo_r9W5Qd{9%D1#q1`a9`Q@h_eD3^zH`$1 z?uRVTWX6?oW}hw0ee>DJ(ul|?nzhX-dwRc>DR#*oduiyRqJB?#GjM8dE^rg3tu2Fs zQ=l1n?}ie-WIfD@sHS4X z9dc;>r^G+yu@+wH?p8yI)B?;*NJ^>_?}Gc()_3iwj;*DXB_2XQiu>0}7dQZlm$WpH zXo%)h&DuoxW0ravDAML>usf9(%}S;nyrEnv|a2n#>r`L7zDOp zgwPD@{wJ&aREepYBsgY4HI&Nh&`JY2@r(%UKK`ojriS=0d4%MId;jgTv}c6-o9Ozw zaQ@({C-{vKgarfr?#wr;uX6!5`vC&d<^o}95wcdL|6D^n4I;+ZjT19<6^}Tl<*x#Z zwDWRiTC>_GEavXQP^zxL7|DA#T*!hQ!lbJqZB>ID`zE>=?<8U zBAJ<0^#qw#D66%3?#q%pDTtoTa^YP&0*cdpcQu6(_xn`y3tlb~)_K#RBRJ-}5dt*J zm8mHI@8Vuyv#KYV;1%;GQ!3k8zlKN7uVR19*p}b_e)YNJx2_eB3HtqweJ988YI@9X zUhB44gFlPsZU6Z)Bmuhf?>1n4QTv*=29nhIWU~2OCiXMTn1Sfyi)q$ z-?(ie)Bwf?vtWj#V0U(Y(ZNva_LKc|hl5cXkmLfFu?*$b1T5I6gUl$ z@DrI43*y3BEEB6;9M+IufG!#r!Ecs;rGJhm0e3%HtL6RhWxDXETfK;GP?`N}vbKya z3?V$6#Ns>4ALda2>t>$F>#RSt^vDv{P=QoLAG>!s+%K0zOTWXi7JGQ28d?Kwh73x$ z>c|jzgX=xpgsrVDHfHZ`GY5Ni`<4LUruDXl8juwdA6MzZtp-slN%~Qj8@r6O0h8h? zrEa_kn*Rl=yFV~`nCEwu6cr^hXu>w^v+&`K(F1hfPdX@+XCe6tt^Z0@Pl*A)&a^zy z_v{6g)F&Dpg>sK^xt!zk@8|%b@_PGJBtvKZ(a{^XWtKA|R@elzy^cXdEv0+dPA1gy zbSVqd&@++ZV(mfI#KP|7T7P?U^Q#=ij1DD&J-Z+$=Q?g;IdDxbB+vg&aR_(2t9|a` zezMwwH;q;PE}5PEG@v!KzGpd=KcDvQks!Lo#rqYbWL8(Ij^k-o(7ho4LWgTA{pd?e4^P#@9g}sr-8!;OkLupihpB7Ee_h}n?(XD}5JzVE9)}L23>I8qIAB_b z2S-@U&s3l`b|arrgIrzzcmGJf#x6I-5}+zeePh?!$R4|;^Sn3tXJCCqxJ-QM?X5dw zh}zRUBxi$4IXZ$GnBe{&;CZFH-dlI0aLH8APfv*aaGGD4fb=I|Qh71q!B@UgfXd&P zUt^*5!Zeuowv$XaS3SShALr#y@BB7dV7?w*m#EQ1%aQMz7$^?jfd*<)I+oh2UIqvw ze$J11dieSp7W%zlh{J0IL}`GvDc{^LPIEk>Zpx{?w0=L+Ew^j&A__2b+@9QLqf5Ej zD8K?m?{>!mW=65;uN0Y4-XTsGFtg`!k?8kb(@X$;&MDS#i~e+4r-&1KRT@Xyt)7>xw-j9~}Gn$Fr%=EtDV{goOlzjyBxL$cU9S zSEyZc#k91ee~_`el)bI^C$L@fuUuIsudA6n~XaKBpX`Y3w(#Km6p1X8Gr2GSoM_IHV-3v@%Xb} zd33QIdXnnPH8lQPhy@LZf97^mG#5KN7R6wX$}>Q-^dBG8>8h%ZeF|wLeSdp%)KUu5 z${w-wkB=3SjTQb$2I%78htO$ji*@;B4v_L=Twk9&L5pWnZx!cwGj=q0qj1<0q>=1B zvTV?q?&lJ@VQn^?qu6RO7GFOF8nf|tl4_AIap2@$ zqQDL)Z&z(n93A%r6D|^etA%d_%d*t%EJ|IbtLD{TTKV8YA@Ck5b4;pB=InecqsPP9 z(ZhVu0DHq*q*Y?{&>F+iU7rbwg!sgb8XVt8hI;>}kZKalzyDr{q*^4ScoqB`NYoF) z4%wCFkYie!^Ga9CD1-YLm=h{*#q7#%ZR#obJ3ODE$vU=wsGk(PIXeq%vb!8fkc~N3 z{ZS$hphJkzB5Zt1-Fk1A3C8+PEH@BFRo6&TEkde8+bhGES`i=lHLH_;Pgj;R^NRSZ$jFFUol}G=qhH2OP$g^snj@2V)r! zV>>aPBQnIqvoo@yQzWz59RT%_{cwZTe>UUai^&j);31UxxF=vqr=l|=b>mZ^SzcAx zBzyE>y<(nxh3Gq^tXPNNp5EJt7E$*N+Vx@qYNe8ep$!)R`r^Xu%$#aH{R)kToqqiC z#jUOON$WI3QW8K1>I1!oTo{TraY>5af<#>XvqG^0lo#H6S23RqQYA2IK}9+Z`$-v=6e+t?y`K!=-%iSe+>y4Wrm?lDIcvBZiJlT3+e2CMVs5}(Q#1(QmH#oS5xQ#r3Qtq{TSx=< z43Tun2z<9E;T=R*?rLV5R$1qfxS~yT+Gp7fm?O&k{Y@PPLYqvO6YH!cQ7v{D?H2G6 zOD2}>{J`uZs*XdVso%jwKeeVJO!H%Ivkd6Ob*z9FO%AH{D?M+dfqmw2Gv15%m7?n zGr12i7LPF!i)!;jDadeS9pC|EJ0KthdJ~AERAe#!AiZ3}Fn5i4Pi-6fID7(}kE?c{ z*A|Fr#ilYw&r?8Nf?rCSIB6OdBXVPLHeS82IcKyh-cn#ja zc#3DA18^C!G`58H!&O?;!qbe!b#L%ifQ6yC`N0dQq4;9F7hfA+4q#wxS#)Pur%7?} z_xNFG*D>PZq`U`L$hCEJ8W{2P!tRj#mvccz<|$kV^v+QvI}sc_zna0#LUR~Yh5iE>TTy_^Y8ep4w0 zL=!jnOl=Bf$03To9fglSyt=hYudg#@(8{cmsHH-g^R~c2h5e8$_2VnmNgaa+RYm<% zPeXg1t$`U|MyIkKe{s}&Q6!S(QYqFOcMFhM@B^fxK z$)V`ki|^0AN8CJ8+$Y$hHm)r?F9Prn&KA%M?XfMfpJppDPEpXafeT~;uC!iU-p@<} zRoW%Y`?yU-zyoauR5nvV1-ET}RKUbR*-_qkt&9P}!&nFW zJQ;V86I()~yWWxFa%lhhrqGE+!<_Xleu2yFnzFaYkSv#D1ROgX*G=OTzor{TK(uWS?uUpLcb%c}7g1;oA zdhw3@^S%JR2C`;p^;mYK?D!nSaiHQv^ppgb7J|6|{M^7h8Ui5|{$bmTrg`F!I)^lu zPRz2-@JL2f0=h0BmLQ>yx;#3{y%i(VC*xw8EQYb1pfgeR<-?&Cil!{n92VUb@9BYx zT#(LV5gm6rj5wUfJ#Jx1ZP-W%cYKh8b;?1bML0#&$#K$8m$|2O8Xj&L^jm zR)E1lo^S;mf^%iI5~y;21HNm~@mfqm@`?%j{xU6zi3I! zBAN@PaT)ZHq3p>w&9+iwv#=_e{!Az{j62nh1#jA@ID~rTHUTZbA3zZZDV>A%QS$w}xNgkviJkUJu#>&Ol+Kk{uEGs4hmMo+2 z3Ml>hH)pCdUm}zoV`dg+Yh7o|*#!oJ-N&uDm+N-ypVwh!(Z^QUx`r^o6;!Rhn-JDveV7+Q6W9T>E-%7JxmFVNZ-F}0L9*(*N z!+Q$iZ-!#fy7aUW7~4L=q(&t0Fo$cu0!6&XCgJg{>)5V5ilqRDc#dg|Wzf z16V?oePlcpDpeB%cG5xZ`MsHB_!eCsLdE{7o&ANm%l->`wXv3~P&?ztV%Ytrwo6CJ z08DIW_9_Kv_BbIy3IX}DAvxi{LH8Mkn|;`{kN4gE`7Yzx4G&^-fcw3ntDPBpFABLg z&P}G~k+|jh_Vc^1h{H=I5EE{72es1PUfQjc=&tb=D32YODZ&Jt^%?oB4G0~OEChY? zIi>ve?JC%pG9ox6;$+e(b)6rw?M=@Gvzzo-xb9`cSvTN!EFkUucpu?f>79y4iF*@- zbOaX@#x$;2Yp;Sf?%L%>?RbaH$lx!L^=r6ufp+|QC;4al9vXu9B37Ze8DKw-&XEtT zH=a38Wv+K0+qfg}Y_4l2fADpx?~&C73tGC=>mmBOu1@0LJkPC7EY^C;t;O#6t|$nS z30pz^WWrh3;X5oxNO6N{T~<2}9^id)5D{`O^3}Y)H48iNM+2Pwwc_H&Ik4{CdQ)k( zADVN#Jm~xISc(n48od15&`o_M7#pU~!Q65@6V!H00kf*_B?I+rExCu&?sze7#<*MF z)_!U+>x;3wbo`SPoICpJ;ini5X0snh(REue6ubuN3e^hzf4KluCyg-xdEDcsuPGY2 zPK&?xVtKO<6F1yot79*pY}NnFTl&r?+7zu)@L5bO546Gvcju4niZTR`&dYb#M@_Aw z3Rs8pSC5Zvg+TBW`_vLDRFy5-Twh<$*fFv`Mmn}X#GSs3aVGHjG=_Q6{^T+d#h~_M zIVTMEyor!0xE>e*N3DMiOGCCS(%^89y|i#+%cfe1?l$c+XF~b*gKzAxt~BU+9i+=M zSykzbfJf+0J;nMMaUjo#yyu_#VKuOof~$zIR@`vx01YzFz6CN&e#1{wbH>J=#p|Hl zeG&$&gvy55Ey0Hyx=+|UIcAV|#HE+wVtb&8Fn77#2EuibAu1IWe*C+u{~CDNrV8Wm zrHN_ek{#}j9REdS-s1fsaTaK^Abr0@$dVRQAe2{_(htl<`b%QJfBfJOHXV`%Dy!}z zovewAn~j^@^JN6#jCyA9=Qlkwf@Llquf$;1i}l&NLJ5=(TdQKB^P2ZHv`#xTeTrqcJz&?>(7KGQw8o>vCYuK^8YF|4XF`r-rf2M zvH@tbAz-M{f5>RsRy9feG-zH*SXWegUXhzj9D&SjMEK)I9ILalb7a}sfZvMmI1MNJ zpG5b%1(WuMTnb z?kREm9F_FhlN59_luid`+{ATa!e%~7qJ7Jiz?il&$JH@vM02Ey2@Q@cOmdBbjfc?= z$OUU#aprTEUbjmcHv(%)ACl@0FC$l$@r0wJ8`_bhcOb2;U9`llWG=9CdP*e!1qh-q zpLDut{q@R(!RPxs=DAVN({RTOF>vMdFE6a+4HGtLvPRbV%=L-8x-Au}Iy0bF*hql& z|LnX`3#-xj%nUq=Lnyt)>@3u2;^?@B@GI%*<8U%F8CCP1?Z1+4z~7v3z-Z*XgQxUg zGzJRE(I2v%&*=_Kn{XyKoIV~hFh&mPD9jW!rs$NhiC$K^rwr*c(1TIXzjE2z{VhGw zA|fIoC3XCX6*poKoi=UC8f!Q+*OdHj&UzR~x#JZFjRZ_ZgqQ>lCjS7c%3Wm|qK>RA z|B@f~aNj0Ha|JgJw-3LUAM>1fmHGMaH0ts)&y7<^K>C>vz`cU z&!TaMECFNBgmICa&-R@y?@{O`1T_GEc22*MzDs5sn_ zEA6@ACH1gF5E4{CrpWP1bXncl$7jDibXa=iw^(Yk9iRC5ci-jC^EW_;V%gBt<+)}6 zRfW=@KYza7W6)RiwgoH&dfqrO*)VtuV`oQyQ1cqLz&%@~;>D{Uy=lKGWssj%a&`g|3jOPG8eV?|o<8m02n$IvWQ=@61!j z8rZFE#N8WR?cW7qvv{8peEJ^*^v|+!``4oi%6m4L<@AJYN)UwpCWfbNX}PbQqW_eg z<|D`b-h_|7D$F)dh*6b)fx8^?=6s>>5)Vu!$-B3I^~_mIh6A_bB_%cgP`+cH0-Z-? z?_Y>1MpJJlyee{6YaepoP9d834vUm?XOrvty0JyOiolYfiS z8|DI_@fwH~)ShegugJ7{XpWtO$%)Icv?ipUMmISfpmW*tQ8+nWe_AKXE@3KPZ&WmL zi^V(6t5aic$&J#c3m;%q*`>MIrsa?`Q~BAe<*S)GIP)>)C*1I|l?ej?23iT}wI1Jm zl5ifZq#l&r_<~&8`I$=0+)C)ia+Wv6KqvKy3>Cb$m#Xe`R6KZA`M^#B5J9jC2vWE} zzltc{k`}=**qWw5%n)gRj>(2g{VXRwZ7dHXA1#RrOkz@64zoC5A0-#kAF;;zW-bFo z9<>=xC(xAbi<@`_GOBlSRTdhVm~2Nw7n^iIF{;`b_MIFCzn^x)C$bu~W!g8G9039C z-UL~dH!;7OhADBdJZx=l-#z=ze$*@`J|y*T0W~N}<(B7I+BTXs%U9j5+hTt8^3wYq zG{U70)gricM{2ZKM55je>E8(<&9NUmrP+#QtYI={0^3sh3{D3h!prtu^T~a?@Ns2Pw^OXO!9HZ zlaHJZyAZ^Bq!&6|gB;^7&fnF(@phR-O6HqpKn9?KwtxNy5HegfKJ)L%-#+KxjaMM$ zy!Lsuc*Spy;6xbF)!M%BCDlBIe6q9D1?E-1(C7=TIO?*2hYrZzuyN^;98we++EywN zIW=&Z`As{pX3%l!7R>N2y;)W;-`2>sfW%gQ5{IPoei&$P$rMS2$WWYoJN-{>>+t?o z-6HSFiHC=A|95tccQ6!lRti(s%)$c0HC(R0biy?i*>WAa&0LN^h%8g!Bl9Vu>S?)Q z4fuLQ$QjR(%#ope$_JVjznzes+z_3KGmbveT9H?_Pk6l4s#j6W#u81`8_56oMy5#M z{Y)cwp=op#5cNz;S%?8rqe_`Zo%lF$BaoZy)VWfXne!edHmgY{&>yI z6oIF;&Ql6M6c_*Ib-FFC2Ap545;SI&wA#+P-YgC$4-@?I%$~ULEd3K;tZ&%TzgWF3)4FOW4_HfhVWew*e!4m84BzxmvEg zEY0Ti{9i6k&eU`)S&}pt*GtJh1px1VfM5D@K_kR-LwlYLmAfvyWa|yzPObx_B;NgP zJHV{kwHV&9(;t;#-G41!juaG1cIhIttWPZ#MgB2!-wtT|D)*`D=+Uz>-l0J>J8a}6 zSFVKN|XIPUNGTzBOfiW?7uo^E-sx1LQ( z%_5M**QfL|6Vmz_SD?LZj6m772N#+Pj68qCY_Ss46mlqtQR1GL{JVcm7(Q!X``#xe z0EI%$o88wruT9Z^Dd3*YqG<9stbleVzk3;xpqHa+%x33d^9u4D(>b$C_fB^cCHZ@6 zQD=-b+#$`52 z*i>pLnlVtk!Rllr4Fv5cVF3=xC1i6Z3Kp`f6rigpxk_jg~{vn;3_CM{`U$nF*N zpK_-7A#g*5?mRG^mL~53hCVkpH|b#Kk07;2D~+kjR~lX(Q`vvY%lxu*fY(?pHNfw8 zM?wqo-JYc!=*ds}R*-;z06Si1AoLL%Bt$BA+J?8)yw|hNrO>v}&g|=!;57F!eYm}z zaPa{8K1RuC`Z@pT9K9`V4*j!s;6b+UdakU3#?$8L$N|kxAzL;!bw{VokQY5kv<%Xd zwRd9B-)={FWTh|GCK7I$v}+>HAW?)*6!{*O`9rt6LiL~RK1>?p%V+H91`Pfu-o)lX z+~`!|AEmJuzVSDnxrp$a!9~TfSS^S<=#oMm&5wy1aD4bBrS_ch@kl^rK_nZF@*N0uPm)Y;X>i$?I&_Bp~N54*4Oh;ok z{)O&OIQ*q5$VwAs!1@|n|F>z?@au=JZF+tB*dW;rd|rdpfBF~4)c#&`xxxM}@)s>H zzR0o_$_M@8TOC>b?GYI4%JcC@IzE{eq28-o#-q4@UBaX-&Ew~(2Ji$#dDL$$?zj`u zjKy|;_$tPJ>@|32QBA5(2Krp~Vq$C{idGps_|WF<2I}Hq$9@G;tR?iFYer`x z_=H?@((dX&JD8}o=a(8~8GpbVX#Z#U$wtHbh529}PyDk=%z^F*)HK)pwdAkrdtGfl zuB7+Mq{lk$tmC(SgxsT-erU5a@o+u)Gi;NECmh<6VpW^%I-+#RG3}s<}prQvY&vBd=@Y7E`hDN^p@G7GX~2VrI8Bwbt@}XMzT2`FMZ#SqfBRg8bFI?R@1RZ1e;6LI2GI zbT=Fuiw5zUkPa0AWIGqg_KrUQw5h@Oc2(m(=q6jX%g7IDhb|U6XLu*$#@~H>k5aS? z$Erc68#-cRw#X2V?~3MK%dgvqyLj6GnAI@N=!+o~jH`?0mi z$MYl`ruMy(Zvf1`o~S=%%Z!5?ADd8ozkBuzJAl3?ThsFMmYIT>cDU=ls2XvoL zY_%urHx^O^(Wi9z~(d*(YzmI1G9>(wd?-$?{3K zZUyRzZ+QH;#psV9PMV^~w?1WA?QifNm3p)QdnarA{I{ z1$@Qq7EO(7ZCyETiV*T&=X1j^4N+D6ib#%HFiX&d4rqQG|E%!o=8veE@S=v~``LG# zqraP5jjL{|xaETmfpnX~U(RXm073o}6-svDP^wcmsn?sb%oLcbJ1D&$g67vwP+e@? zgaUO~_0((BhZNswRd-k<8bz8sADPs;BaPD%~xLFe^mA5zFJXoBB- zXGEH}YWy{6)NVEEfaTgsE|P3^@8#X(Ykonfrr)-!(Jfow#XF7Zhi)A5i7XV`@HJ!F z-sbO)oQwo;$FTZm3uZ z@mUE=<()s{V}ZT9G3F`Evb!kNQK~Oj{QwshrgMO&K&Zchn7CE;AzY}?eaU|-ebfF* zzm#?ry)Ij96l)QG;ccrS&R^oQ|FvEY*)XJc+Jq{sVi-?E$~Ku~pd>I@ZRE`=f0+aP zwcKDk{(KUehdHoSB-T4(x7g;?LMO2{R$GhdtH*9J!}={o{hdainY~nv-BHE@F!`@L zO_`8*hmYQAztQ6{NhZss$F7P(Qlm{Ct+8@5|3NfmK3`*ttFGfvV-+t-Joz|qsgdJb z>4d+B#-KM4Lb&o8KG9z^KC+^|uL8uOywX;& zUbU+6i2|e`J5`Z&$LJ&c+QcvUAO2ox+bGKN^9Ae`x>@&98=~@;ObDtJ{XAo##>cAb zovkM2nJeed)Gs0A=rYJs`jMOv=+cu7X1e0Or`@z{JJTJHjA-Tzq+@)-n(DzP`dz=h z`2)xRSNFC(%LAfIgLwp&E;3=PHQC&q|5{sB+#LiVHz1Y*7s(A4D`;Z7uu+VVxx+u4 zWHw2=`fRW;eVmxk5`R9-q(kf4f6@6y3}wDJ9&8C@;u-s)Y3o1lTd3?&;GH+9K`#6n zag!-y00%6KDlAE0NLxhJuFyW77Z2>Gl4-jx1yAV(-G7wJGS*-vilf}rmMC}2JCxTW-KXbWj&vb#Hd?4>k4TTZ)7I|^(f zcZV9L6Ro8s$W_=;!!;iW0J8PQu2mJv23HE0)%f<+qnww^rhhEDd?-KG9<~^V5EV zsvjPS`nNE*L9X-^^P(c#CK;fc*e8%<~Bh^A9m8nT7_?fJA>aU|z^q6eFBgHpERR}NcpPLk+MrNZ^u zm~?z+>Di5^pEJ7~Ga=pu$@+V-qB54HcKhjEgv)-y%~aw=+52|lsc*9zQ*VZftu>Tw zWic$%deK0nojKhAxE-ulCED3}Hjq`aE{5b+)7{6OHPRAP5}$+lU$62*nwKN$3%KV9 zcwtPP{P)kFwb9iorEKGzWaQBdr$UAJO|K|Ew=K_)H(G1E=5Ld^sI>4@Gv;bel-*?2Gz^uoTeK@?WE3#mLJfii^a->(pKuCoQdIr?Z+S;Z z%JR?h)D^iPMjMIJ>Z7_za$<(~g~%2fg}gvl7n+@#FT=n$rW@Zsz`->q;~g#&=wYTD z8Yb0sr!91#4ZK=V4v!oT>@%g`CSjNg9a9KG4rO>G~o4|?)(i`Gy_2b{DHn86RQHvw>NS` zJVstPV}l2*@FPdxYv;@`S^y^+>ERb2k5VaI@H()(QqG<8|E4l7SZ0872X}WtWNx_p z^3JwNv`AmH!=BVtN|5?Cfq+mF`e1N^6kH?*tg0#xNO%Rh-quBKKNX zW9m%69Pc?z{G^l~ZLEzEWoKQ?>c30~klZ;KL&N~s+G{?g;A(t2-kA15#0Njgoi@|c3lH4>wJkf_(cA^kAD@Ay`3Bv^ zeL@h<12<|FM17jlBJ*+lyV=Fuiqul1)T?S{n4t6%I?`OPqD#R|*R`A(L|Oq8hH2tC zt*gsMz_;*SpE<6$k~n<<%qy4c)197ruOhFOj#yVC+_bHI&B9~Mu_xucaAw=&9asgr z)`Mn!{p7?YUj&XuoigFRFVZk#Q3k?NKd;V2t z>D0X0B_(aawl@ne(cS>_(0fxUoAodnc_Ij5@k$r|KX z_$X6-lOz|j&Z%iZoii18aPc3|HR@Cb2A<&_lwwy4LzQM8#y@zEI>3hzl&>#mzGK=; z9<;{0=zqb2kjYM(TcB6yBFn9w8=!7apwC7o5Y8sQEMr%<3X@rq6mW60Gpj6J8kA`9 z_cuR5+$;aYbyrY)>bcC>%kMJVC;qO}NW*(EWN=8u3X9LcuHwB&Ky?Et zNc?@w%}V1l70@N}r(wCqKj?#WxC{e_E&~TaV`mF;po{_wMgbxLBcgYwv}6$U%bwFH z!|x0CG3wk}nQlM-w5xXD%fQTR7_DVfE8TnIcX_9%V2_^h{a4>-B~rG&{qhcM^!Y$n zA5Af;cL2;xxB2#Yfjh0#VJ289$tb*@>JN9*`|+(4sdfR7!W~gGrti;&I1Q*q^U5dG z!?t_8aHlAKW+H};R!R?Mg%`wtF$UFtVR>QFOG7r09-ezrHakE_ZQ|PDjgup$eCrcb zIHNJ>eUCW<*r#7v@*g@-*1QrOLZ-Rk$tNF0qCKLvq3N*Pbbyx!_d&yQgRfWrZh?zp z6)gG7XnO5j^D(-x#TsHHh1V8=@U(Ay?Rw&}YPBouesnv~#|L%}ZdUM2@QGPTOkvEg z9$)Va2TVbu8W}9D`zu&6FPvqUoaPRTK4aJ(x&2r4`*zvTR_kY6=Zuh!nE-||yUVx( z|A%jxyUBTe<)k}Q5+qLzVi@d+EjQ<#`^pbiP!fu)|LME~(_?)9Q27fp@!IF4qV}b+vF72<5ER z2~H?CQ8*7?-($mR=`zk;;=@B2hnP)-lNw2Lb4JNY!l@cgq7%_fO43L<#~})#sqeAj zv~tM~JiUs=m%!)K%PvNfQx41FI1SQvUzZGQo(cb3$Ugo4CTb+n^Db0!9?0VZ8t26+lY}#o(olOYMRiedGMf~k{}6(2@X1)UQ$A9V6Wks$ zm7EEtjdi1w^iCz5hN&+$>7$&bRBs}kn>A(gB+41E*-k!7v>eGila{{sTSXu6_KPis zaNgONNNuEC*+$}zgER!fSpiq8zO~|CYI1t^7yWH*3z&weCB)da000s-Nklt{x^WeKDy?C*>5f zxmiy*e_p`OV@Bc$PD5Q>ai*l4*_K|B-@9AEC%jrB;hYe{nF0R|z|*Vb-VnPwvWWqj zcv0@^UDZ-KjX1fI5(gLE)qt za6k1N)Ar^c+uj?u_K*1Ga&~{CtT@Tl3C6@dXG=#Gan7K+;+g^n;7)4Uau;fHI zzW`N4J|EL9DfzR~@A3Cz5>E5YhkQUojwE4L z_XUM>xkWvkES$LXffI8&%BGfZzJD=+ZE=7Tg|kvdb;!2C{0^+XcR4ctarzfPDsvB?8$0Lx3I-(+iKHFlB!&P?-_u`;Wdn74<5Y%VXqQwrxK z%O(itqE2w`5W;y_J2v~#ZLjua9M8h2-pXkxITM95jfAtX@b>M5wRMp$o8Y`deWYvc z+<)()y`kjHl$4-w#>F8kFO^DX48xfmAv^uzGKK%OcWDw%rEG$7Hj!>Z$*GsE6Pi(n zFKJhB5;FiN#34v1tsrGi>TwdA=X!hNfXxp7n%D&4%KpRl96SO>PMgi#EKK(1U7H`{pD*)Lqa$wICzOTv@3=-)xRWD=crx*KHH0?-r@Jx z*Q3* zIAsUkys>T`NdiyzD0cNQ5>CM;EIletp_?e2n!)+LvA+$NAy;pwtPc7j`#KPZfXzY> zPQ2*%W9fYROR(gemfPy)KGLc-62ckVAEXMM-9gtm4nllDC#k*L+ zR#HpLvZqDyPg^AQp&DYdrskj4Te^tp046NW0@3^ae(ZAgZN#izwTp1s9~91BDV)B} zk&it~^LQ!PNWwwR#H=jZRp*kEIs+*J&T|AOv00C=S*U^Y$+f}%2Pa>Wgd<53&IYyj z%SkvvH;067*3_Ijb;{x>1iBf8G5s(Jr)IYPIGf;{?!zW$I|^sDbN&8gZX;CyoS&Ab+h$9lHF3tnHnwTk^a-(xyM9x<#Bvk zmD=cgLd175V&`X-CvaB672xf6tQ@y^BY3Lcv5)AtG4%rGaM?A5Wh zvr2e%NYxWxm#I2LoVcx(Z#r^c!1=Zw&g}^2)K#kvP6M6J2wB^8Id?ewjkw`Ix_tUh^yviX z;?=9i`MKlxCdy_=h$fsFrW4L?f^#{-iG?J5f(MT$)5u3+qNs6Zv2gz0$326b9w8^L z&miH%a{obEoq=qEaOzTZ@on&XSK*`8SRr$R^H~o#dDB>`a>~+htoMILrz{13~*}8Ggk5VlZ5j_eX1_tJmMdm^$2GrFF)2Q1KJtubB4s6`SRd4 z%*L(-2u{FLswTde&5n+Qj`2}VT9FaPbcUGJqqD3*EC`_>dG(?hD`@RdGNwZdUP{D~ z4|n181(a?IIR%^?a*Fnl$kknJZr(F<<{zFH!kH`B^fTscdqF7YdfB}caEfKYb98Y2 zQOni|PLEW5px40r{v&(!#X>j15)u+RI#L7vz#k)<(xr!m$gZUc(F9$)DDokCgbXKn zb%cfLNVBwv(+D_4m2>-&6)Vr@$_ zQ#8b+@i1`mN!Gc!mVgN-n#j|IZ~~a{?{I4B;lMucI7X94J-G0WUfp6&ySVpl2`A)3 z04MirxN`|`?r_0*iQfhugCFA7!LR=rh0~*1+uv{Wt>3AW&5i^@Gc`0c^ikj*c1&D)`ZCw_CXONP z+m~TZ*(ATd4d6Thkr~xnxYNB*6diGpbBCkbW0p0-~kwiJN`<;8=?k(Po5 z3%KA9UOjdPSFUoUaQeUavG($Y6>@PrZ``m|Q_iI;b=$^LB`F`Z6W-smacVrAmIs)p zW10B2m)r?CQi-lx))0}ybJ|zN0JMq2?&M91hXRjWsgA-8ksa@m5 zu@4uVxLfQJxgvPlaq8j&DpXLXvTCS}Yef@IcGU63(qUIJL8aMQ@$pEPmDE z>D31UaPwI=Q8rOG5zTdAr2)ZHtWM{4x>?7NJgxl;4>&UmC)i?Dpy(@VGx{_u2Z}+O1I^WT-_#1Nwat=DW8?o)2RnSI2q3l zL^DCM8H#AatJHXpF2l|EFsYmH`o&^$t2=>V zrtN2QNnB%d3U%?H;Z&XKNP0L8Q%!f8IvyJO%ZN}OaPCrS>MsNmLjYL@|-bEuPvpz4p%uh zZSrqAsD}$yyj*=d@5qS}Cwp~DhiH)}gb^Na8hMpxKL3cfF(J3UmKJHa&e5P{>ooiE z+AwVP2lDodJbj(e&9bu63V^dB&oJ3Ur>8v ziC+Eo?c3L`-SBPX*Xb77-jc#e%he^EJdY&S0}D2Fa2h>MX*vGQH`x9f8BTcjVFvyO z2LpWbLxod^fU^vQvm&ny*sKb#G6GH7L{8Th%E&;UZkdSob!)+w=fJ)KJ9gHu-L#zU zvBUe-xHIkc?cC^|Dof%8Zf$i9&SGDy1EXHxuKDxlCqlMPtJ{|Hp_j4Zn3t?!*UEPf zb_e{vi+B^zgov{&4=g+{JUqgGG(k7h(|LzOPK7W}w@g^;Z2gYR?|t6f>8x#Q*&;?h zcxO2pGp1oVs=!OHjt!#XTz?~pCthFc*0)}H<)=Bk_;MzeS1If~d$wt~Yj!qwkHcMm zguy1zSpo4Tcy$!ch={6)6hHpt<1}#$fvznSV!rjiL9%+!lFsFump5#HY-p+&OBYkd zuj|@FSy?^@5@WH%d*B}Spy5~`*zMc@8DRyV`>8{N;!$K5DPnGrT?HTX41eXrNiQ#5M1nH)&5&ml| z^!_g#O4A9aXbVBO`HOm-!87NO-}v2gtlQQTgbRUivx=1SyzK`>voWX9mXlNXtG8&< zQI>JKjwK1F0wkDDy>95Q{2V1}ztt0z3x0(0CWo8n*{v5A<`hOodV?k%o=(-xI6jn( zOGt{Mz?io=lJn%y-}U;jx}xusPwIko%_*3u|G?p9V@@IcM0$ZIp&8>k%rf4*hh!)N z@$^+a>)Zd@9>mWrWU~=vvoQyClPqa8*^J@ysnexX2Q*c%86bPr|3A_1%e-4(b#--L za4+2^JUS_7WTe2;KC)_B(^xp;bh>HYQ`sgPnpa2oeeYPuirTVg@2?N~)&HEs!l`F7 zGAU_z@qBExf^f={)6=z8mI@++ljJdGLf$y{{sHf!&)TEe+57v8t^>KlVhkoa0ThbS!Vu{N_`jbXVsU{oB~T==6`vVxRZr>>A~BbvM&e zHQ)XDXy3AKhr`hvwPwxUp1QjFK9@k&cD1z}IPhF$OIw@So=&7x&ZHzepji`r`gHWC zU>cj{<%DLMEFf{u@@Weex88gB@cZu{9`k?t!a0(5HvErm;Thd8um9#NcwaVu<4o_q zWsYy(c;lb9s?lz`>~fnqpNn+r_N3?<`)T{9Y`b?8npdnGZ%W}LGM8DDJ%LWgXb^R_ zw!Sq#@!Yu&H~;jVk+Vo=##%J09bM3SO@L<-g`4(h`)P3M$R^OayWEODhcq-}L|aG! z_(ZGorn}O4=gyr@qu;sH_=^cToR(iP-UP1>K3&MB6<*-auasM5brZdMJg#I@^o>gA zx4}wzJI_6JY=T9DZ0dbFm=)<1B_uJ16tB?q&f1p;t&rVhD{v zSMZglbed%CoyVC+z3%9WdDkqHxn%zMD!&hq#2qD>9U7j@>p+ z$t>7n_=hQ67Dg znX}uhpNw^r@^mh5s)senY+%5|9}dDzZAnIqdRCK7HV1pn{Q<37xTzl2Wb^1w^M6PV zIMu^a(WV(qp*i_1s!rVlWoERxlj Date: Mon, 23 Mar 2026 16:03:29 +0000 Subject: [PATCH 06/36] Standardise docstring formatting. [ci skip] --- src/loch/_platforms/_base.py | 11 +++++++++++ src/loch/_platforms/_cuda.py | 9 +++++++++ src/loch/_platforms/_opencl.py | 9 +++++++++ src/loch/_platforms/_rng.py | 6 ++++++ src/loch/_utils.py | 10 +++++----- 5 files changed, 40 insertions(+), 5 deletions(-) diff --git a/src/loch/_platforms/_base.py b/src/loch/_platforms/_base.py index 1584a4f..27b5380 100644 --- a/src/loch/_platforms/_base.py +++ b/src/loch/_platforms/_base.py @@ -55,6 +55,7 @@ def __init__( Parameters ---------- + device : int The GPU device index to use. @@ -89,6 +90,7 @@ def compile_kernels(self) -> _Dict[str, _Callable]: Returns ------- + dict Dictionary mapping kernel names to callable kernel functions. Expected keys: 'update_water', 'deletion', 'water', 'energy', @@ -103,11 +105,13 @@ def to_gpu(self, array: _np.ndarray) -> _Any: Parameters ---------- + array : numpy.ndarray Array to transfer to GPU. Returns ------- + GPU buffer Platform-specific GPU buffer object. """ @@ -120,6 +124,7 @@ def empty(self, shape, dtype) -> _Any: Parameters ---------- + shape : tuple Shape of the array to allocate. @@ -128,6 +133,7 @@ def empty(self, shape, dtype) -> _Any: Returns ------- + GPU buffer Platform-specific GPU buffer object. """ @@ -140,11 +146,13 @@ def from_gpu(self, buffer: _Any) -> _np.ndarray: Parameters ---------- + buffer : GPU buffer Platform-specific GPU buffer to transfer from. Returns ------- + numpy.ndarray Array containing the data from GPU. """ @@ -185,6 +193,7 @@ def platform_name(self) -> str: Returns ------- + str Platform name ('cuda' or 'opencl'). """ @@ -197,6 +206,7 @@ def cache_hit(self) -> bool: Returns ------- + bool True if kernels were loaded from cache, False if freshly compiled. """ @@ -209,6 +219,7 @@ def compiler_log(self) -> str: Returns ------- + str Compiler log output, or empty string if no warnings/messages. """ diff --git a/src/loch/_platforms/_cuda.py b/src/loch/_platforms/_cuda.py index 4047115..fbe298c 100644 --- a/src/loch/_platforms/_cuda.py +++ b/src/loch/_platforms/_cuda.py @@ -68,6 +68,7 @@ def __init__( Parameters ---------- + device : int The CUDA device index to use. @@ -134,6 +135,7 @@ def compile_kernels(self) -> _Dict[str, _Callable]: Returns ------- + dict Dictionary mapping kernel names to callable kernel functions. """ @@ -198,11 +200,13 @@ def to_gpu(self, array: _np.ndarray) -> _Any: Parameters ---------- + array : numpy.ndarray Array to transfer to GPU. Returns ------- + pycuda.gpuarray.GPUArray GPU array containing the data. """ @@ -214,6 +218,7 @@ def empty(self, shape, dtype) -> _Any: Parameters ---------- + shape : tuple Shape of the array to allocate. @@ -222,6 +227,7 @@ def empty(self, shape, dtype) -> _Any: Returns ------- + pycuda.gpuarray.GPUArray Allocated GPU array. """ @@ -233,11 +239,13 @@ def from_gpu(self, buffer: _Any) -> _np.ndarray: Parameters ---------- + buffer : pycuda.gpuarray.GPUArray GPU array to transfer from. Returns ------- + numpy.ndarray Array containing the data from GPU. """ @@ -273,6 +281,7 @@ def platform_name(self) -> str: Returns ------- + str Platform name ('cuda'). """ diff --git a/src/loch/_platforms/_opencl.py b/src/loch/_platforms/_opencl.py index 6395c43..29f4e31 100644 --- a/src/loch/_platforms/_opencl.py +++ b/src/loch/_platforms/_opencl.py @@ -64,6 +64,7 @@ def __init__( Parameters ---------- + device : int The OpenCL device index to use. @@ -131,6 +132,7 @@ def compile_kernels(self) -> _Dict[str, _Callable]: Returns ------- + dict Dictionary mapping kernel names to callable kernel functions. """ @@ -240,11 +242,13 @@ def to_gpu(self, array: _np.ndarray) -> _Any: Parameters ---------- + array : numpy.ndarray Array to transfer to GPU. Returns ------- + pyopencl.array.Array GPU array containing the data. """ @@ -256,6 +260,7 @@ def empty(self, shape, dtype) -> _Any: Parameters ---------- + shape : tuple Shape of the array to allocate. @@ -264,6 +269,7 @@ def empty(self, shape, dtype) -> _Any: Returns ------- + pyopencl.array.Array Allocated GPU array. """ @@ -275,11 +281,13 @@ def from_gpu(self, buffer: _Any) -> _np.ndarray: Parameters ---------- + buffer : pyopencl.array.Array GPU array to transfer from. Returns ------- + numpy.ndarray Array containing the data from GPU. """ @@ -299,6 +307,7 @@ def platform_name(self) -> str: Returns ------- + str Platform name ('opencl'). """ diff --git a/src/loch/_platforms/_rng.py b/src/loch/_platforms/_rng.py index 42699f4..5b80e66 100644 --- a/src/loch/_platforms/_rng.py +++ b/src/loch/_platforms/_rng.py @@ -38,12 +38,16 @@ class BatchRandoms: Attributes ---------- + rotation : numpy.ndarray Shape (batch_size, 3) array of uniform [0,1) randoms for water rotation. + position : numpy.ndarray Shape (batch_size, 3) array of normal randoms for position direction. + radius : numpy.ndarray Shape (batch_size,) array of uniform [0,1) randoms for radial distance. + acceptance : numpy.ndarray Shape (batch_size,) array of uniform [0,1) randoms for Metropolis acceptance. """ @@ -75,6 +79,7 @@ def __init__(self, batch_size: int, seed: _Optional[int] = None) -> None: Parameters ---------- + batch_size : int Number of parallel trials in a batch. @@ -117,6 +122,7 @@ def get_batch_randoms(self) -> BatchRandoms: Returns ------- + BatchRandoms Container with rotation, position, radius, and acceptance arrays. """ diff --git a/src/loch/_utils.py b/src/loch/_utils.py index 3defbac..c84b1c7 100644 --- a/src/loch/_utils.py +++ b/src/loch/_utils.py @@ -239,19 +239,19 @@ def standard_volume( system: sire.system.System The bulk water system. - temperature : str, optional + temperature: str, optional Temperature of the system (default is "298 K"). - pressure : str, optional + pressure: str, optional Pressure of the system (default is "1 bar"). - cutoff : str, optional + cutoff: str, optional Non-bonded interaction cutoff distance (default is "10 A"). - num_samples : int, optional + num_samples: int, optional Number of volume samples to collect (default is 5000). - sample_interval : str, optional + sample_interval: str, optional Interval at which to sample the volume (default is "1 ps"). Returns From ff787c17314762cf467b2b11e147bd65996ea08a Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 24 Mar 2026 12:10:29 +0000 Subject: [PATCH 07/36] Fix CUDA context stack leak on crash by tracking outstanding push count. --- src/loch/_platforms/_cuda.py | 15 ++++-- tests/test_platform.py | 93 ++++++++++++++++++++++++++++++++++++ 2 files changed, 103 insertions(+), 5 deletions(-) create mode 100644 tests/test_platform.py diff --git a/src/loch/_platforms/_cuda.py b/src/loch/_platforms/_cuda.py index fbe298c..ee220f2 100644 --- a/src/loch/_platforms/_cuda.py +++ b/src/loch/_platforms/_cuda.py @@ -114,6 +114,7 @@ def __init__( # Use the primary context (shared with OpenMM and other CUDA users). self._pycuda_context = self._cuda_device.retain_primary_context() self._pycuda_context.push() + self._push_count = 1 self._device = self._pycuda_context.get_device() @@ -256,22 +257,26 @@ def push_context(self): Push the primary context onto the calling thread's context stack. """ self._pycuda_context.push() + self._push_count += 1 def pop_context(self): """ Pop the primary context from the calling thread's context stack. """ self._pycuda_context.pop() + self._push_count -= 1 def cleanup(self): """ - Clean up CUDA resources and pop the context pushed during __init__. + Clean up CUDA resources and pop all outstanding context pushes. """ if self._pycuda_context is not None: - try: - self._pycuda_context.pop() - except Exception: - pass + for _ in range(self._push_count): + try: + self._pycuda_context.pop() + except Exception: + pass + self._push_count = 0 self._pycuda_context = None @property diff --git a/tests/test_platform.py b/tests/test_platform.py new file mode 100644 index 0000000..81649b4 --- /dev/null +++ b/tests/test_platform.py @@ -0,0 +1,93 @@ +from unittest.mock import MagicMock, patch + +import pytest + +# Skip the entire module if PyCUDA is not installed. +pytest.importorskip("pycuda") + + +def _make_backend(mock_driver): + """Instantiate CUDAPlatform with a mocked PyCUDA driver.""" + from loch._platforms._cuda import CUDAPlatform + + mock_driver.Device.count.return_value = 1 + mock_context = MagicMock() + mock_driver.Device.return_value.retain_primary_context.return_value = mock_context + + backend = CUDAPlatform( + device=0, + num_points=3, + num_batch=10, + num_waters=5, + num_atoms=100, + num_threads=32, + ) + return backend, mock_context + + +class TestCUDAPushCount: + """Tests for CUDAPlatform context push-count tracking.""" + + def test_initial_push_count(self): + """Push count starts at 1 after __init__ (one push for the lifetime context).""" + with patch("loch._platforms._cuda._cuda") as mock_driver: + backend, mock_context = _make_backend(mock_driver) + assert backend._push_count == 1 + mock_context.push.assert_called_once() + + def test_push_increments_count(self): + """push_context() increments _push_count.""" + with patch("loch._platforms._cuda._cuda") as mock_driver: + backend, _ = _make_backend(mock_driver) + backend.push_context() + assert backend._push_count == 2 + backend.push_context() + assert backend._push_count == 3 + + def test_pop_decrements_count(self): + """pop_context() decrements _push_count.""" + with patch("loch._platforms._cuda._cuda") as mock_driver: + backend, _ = _make_backend(mock_driver) + backend.push_context() + backend.pop_context() + assert backend._push_count == 1 + + def test_cleanup_pops_once_normally(self): + """cleanup() pops exactly once when no extra pushes are outstanding.""" + with patch("loch._platforms._cuda._cuda") as mock_driver: + backend, mock_context = _make_backend(mock_driver) + backend.cleanup() + assert mock_context.pop.call_count == 1 + assert backend._push_count == 0 + assert backend._pycuda_context is None + + def test_cleanup_pops_all_outstanding(self): + """cleanup() pops all outstanding pushes, simulating a crash mid-move.""" + with patch("loch._platforms._cuda._cuda") as mock_driver: + backend, mock_context = _make_backend(mock_driver) + # Simulate two push_context() calls that were never popped (e.g. + # two GCMC moves crashed before their paired pop_context()). + backend.push_context() + backend.push_context() + assert backend._push_count == 3 + backend.cleanup() + assert mock_context.pop.call_count == 3 + assert backend._push_count == 0 + assert backend._pycuda_context is None + + def test_cleanup_handles_pop_exception(self): + """cleanup() continues safely if pop() raises (e.g. stack already empty).""" + with patch("loch._platforms._cuda._cuda") as mock_driver: + backend, mock_context = _make_backend(mock_driver) + mock_context.pop.side_effect = Exception("context stack is empty") + backend.cleanup() + assert backend._pycuda_context is None + + def test_cleanup_idempotent(self): + """Calling cleanup() a second time is a no-op.""" + with patch("loch._platforms._cuda._cuda") as mock_driver: + backend, mock_context = _make_backend(mock_driver) + backend.cleanup() + pop_count_after_first = mock_context.pop.call_count + backend.cleanup() + assert mock_context.pop.call_count == pop_count_after_first From 065ef75579fc04735745e1190d28ab08b4a5fa13 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Sun, 29 Mar 2026 10:03:34 +0100 Subject: [PATCH 08/36] Add support for different alchemical softcore forms. --- src/loch/__init__.py | 1 + src/loch/_kernels.py | 31 +++++++++++++++++++++++------ src/loch/_sampler.py | 45 +++++++++++++++++++++++++++++++++++++++++++ src/loch/_softcore.py | 31 +++++++++++++++++++++++++++++ tests/test_energy.py | 21 +++++++++++++++++--- 5 files changed, 120 insertions(+), 9 deletions(-) create mode 100644 src/loch/_softcore.py diff --git a/src/loch/__init__.py b/src/loch/__init__.py index 45a2c0c..2c8d0e7 100644 --- a/src/loch/__init__.py +++ b/src/loch/__init__.py @@ -20,5 +20,6 @@ ##################################################################### from ._sampler import * +from ._softcore import * from ._utils import * from ._version import __version__ diff --git a/src/loch/_kernels.py b/src/loch/_kernels.py index ba20ebb..83c7114 100644 --- a/src/loch/_kernels.py +++ b/src/loch/_kernels.py @@ -394,9 +394,11 @@ float rf_cutoff, float rf_kappa, float rf_correction, + int softcore_form, float sc_coulomb_power, float sc_shift_coulomb, - float sc_shift_delta) + float sc_shift_delta, + int sc_taylor_power) { // Work out the atom index. const int idx_atom = GET_GLOBAL_ID(0); @@ -538,7 +540,7 @@ energy_coul[idx] += (q0 * q1) * (r_inv + (rf_kappa * r2) - rf_correction); } - // Zacharias soft-core potential. + // Soft-core potential for ghost atoms. else { // Store required parameters. @@ -558,10 +560,27 @@ r = 0.001f; } - // Compute the Lennard-Jones interaction. - const float delta_lj = sc_shift_delta * a; - const float s6 = powf(s, 6.0f) / powf((s * delta_lj) + (r * r), 3.0f); - energy_lj[idx] += 4.0f * e * s6 * (s6 - 1.0f); + // Compute the LJ interaction using the chosen soft-core form. + float sig6; + if (softcore_form == 1) + { + // Taylor soft-core LJ: + // sig6 = sigma^6 / (alpha^m * sigma^6 + r^6) + const float s6_val = powf(s, 6.0f); + const float r6 = r * r * r * r * r * r; + const float alpha_m = (sc_taylor_power == 1) + ? a : powf(a, (float)sc_taylor_power); + sig6 = s6_val / (alpha_m * s6_val + r6); + } + else + { + // Zacharias soft-core LJ: + // sig6 = sigma^6 / (sigma*delta + r^2)^3 + // delta = shift_delta * alpha + const float delta_lj = sc_shift_delta * a; + sig6 = powf(s, 6.0f) / powf((s * delta_lj) + (r * r), 3.0f); + } + energy_lj[idx] += 4.0f * e * sig6 * (sig6 - 1.0f); // Compute the Coulomb power expression. float cpe; diff --git a/src/loch/_sampler.py b/src/loch/_sampler.py index 2bd98a0..2fe278f 100644 --- a/src/loch/_sampler.py +++ b/src/loch/_sampler.py @@ -37,6 +37,7 @@ from ._platforms import create_backend as _create_backend from ._platforms._rng import RNGManager as _RNGManager +from ._softcore import SoftcoreForm as _SoftcoreForm def _as_float32(arr: _np.ndarray) -> _np.ndarray: @@ -81,6 +82,8 @@ def __init__( coulomb_power: float = 0.0, shift_coulomb: str = "1 A", shift_delta: str = "1.5 A", + softcore_form: str = "zacharias", + taylor_power: int = 1, swap_end_states: bool = False, restart: bool = False, overwrite: bool = False, @@ -207,6 +210,17 @@ def __init__( The soft-core shift-delta parameter. This is used to soften the Lennard-Jones interaction. + softcore_form : str + The soft-core potential form to use for alchemical interactions. + This can be either 'zacharias' or 'taylor'. The default is + 'zacharias'. + + taylor_power : int + The power to use for the alpha term in the Taylor soft-core LJ + expression, i.e. sig6 = sigma^6 / (alpha^m * sigma^6 + r^6). + Must be between 0 and 4. The default is 1. Only used when + softcore_form is 'taylor'. + swap_end_states: bool Whether to swap the end states of the alchemical systems. @@ -495,6 +509,26 @@ def __init__( except Exception as e: raise ValueError(f"Could not validate the 'shift_delta': {e}") + if not isinstance(softcore_form, str): + raise TypeError("'softcore_form' must be of type 'str'") + softcore_form = softcore_form.lower().replace(" ", "") + _valid_softcore_forms = {m.name.lower(): m for m in _SoftcoreForm} + if softcore_form not in _valid_softcore_forms: + raise ValueError( + f"'softcore_form' not recognised. Valid forms are: " + f"{', '.join(_valid_softcore_forms)}" + ) + self._softcore_form = _valid_softcore_forms[softcore_form] + + if not isinstance(taylor_power, int): + try: + taylor_power = int(taylor_power) + except Exception: + raise ValueError("'taylor_power' must be of type 'int'") + if not 0 <= taylor_power <= 4: + raise ValueError("'taylor_power' must be between 0 and 4") + self._taylor_power = taylor_power + if not isinstance(swap_end_states, bool): raise ValueError("'swap_end_states' must be of type 'bool'") self._swap_end_states = swap_end_states @@ -1323,9 +1357,11 @@ def move(self, context: _openmm.Context) -> list[int]: self._rf_cutoff, self._rf_kappa, self._rf_correction, + self._sc_softcore_form, self._sc_coulomb_power, self._sc_shift_coulomb, self._sc_shift_delta, + self._sc_taylor_power, block=(self._num_threads, 1, 1), grid=(self._atom_blocks, self._batch_size, 1), ) @@ -1902,6 +1938,12 @@ def _initialise_gpu_memory(self): # Link to the reference state. mols = _sr.morph.link_to_reference(self._system) + # Build map of extra options for the dynamics object. + _map = {} + if self._softcore_form == _SoftcoreForm.TAYLOR: + _map["use_taylor_softening"] = True + _map["taylor_power"] = self._taylor_power + # Create a dynamics object. d = mols.dynamics( cutoff_type=self._cutoff, @@ -1916,6 +1958,7 @@ def _initialise_gpu_memory(self): rest2_selection=self._rest2_selection, swap_end_states=self._swap_end_states, platform="cpu", + map=_map, ) # Flags for the required force. @@ -2059,9 +2102,11 @@ def _initialise_gpu_memory(self): ) # Store soft-core parameters as scalars. + self._sc_softcore_form = _np.int32(int(self._softcore_form)) self._sc_coulomb_power = _np.float32(self._coulomb_power) self._sc_shift_coulomb = _np.float32(self._shift_coulomb.value()) self._sc_shift_delta = _np.float32(self._shift_delta.value()) + self._sc_taylor_power = _np.int32(self._taylor_power) # Store immutable per-atom buffers on GPU. self._gpu_sigma = sigmas diff --git a/src/loch/_softcore.py b/src/loch/_softcore.py new file mode 100644 index 0000000..6d165af --- /dev/null +++ b/src/loch/_softcore.py @@ -0,0 +1,31 @@ +###################################################################### +# Loch: GPU accelerated GCMC water sampling engine. +# +# Copyright: 2025-2026 +# +# Authors: The OpenBioSim Team +# +# Loch is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# Loch is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with Loch. If not, see . +##################################################################### + +__all__ = ["SoftcoreForm"] + +from enum import IntEnum + + +class SoftcoreForm(IntEnum): + """Enum for the soft-core potential form used for alchemical interactions.""" + + ZACHARIAS = 0 + TAYLOR = 1 diff --git a/tests/test_energy.py b/tests/test_energy.py index 50e115b..9235cca 100644 --- a/tests/test_energy.py +++ b/tests/test_energy.py @@ -5,7 +5,7 @@ import sire as sr -from loch import GCMCSampler +from loch import GCMCSampler, SoftcoreForm @pytest.mark.skipif( @@ -13,8 +13,16 @@ reason="Requires CUDA enabled GPU.", ) @pytest.mark.parametrize("platform", ["cuda", "opencl"]) -@pytest.mark.parametrize("fixture", ["water_box", "bpti", "sd12"]) -def test_energy(fixture, platform, request): +@pytest.mark.parametrize( + "fixture,softcore_form", + [ + ("water_box", "zacharias"), + ("bpti", "zacharias"), + ("sd12", "zacharias"), + ("sd12", "taylor"), + ], +) +def test_energy(fixture, softcore_form, platform, request): """ Test that the RF energy difference agrees with OpenMM. """ @@ -36,6 +44,7 @@ def test_energy(fixture, platform, request): reference=reference, lambda_schedule=schedule, lambda_value=lambda_value, + softcore_form=softcore_form, log_level="debug", ghost_file=None, log_file=None, @@ -43,6 +52,11 @@ def test_energy(fixture, platform, request): platform=platform, ) + # Build map of extra options for the dynamics object. + dyn_map = {} + if sampler._softcore_form == SoftcoreForm.TAYLOR: + dyn_map["use_taylor_softening"] = True + # Create a dynamics object using the modified GCMC system. d = sampler.system().dynamics( cutoff_type="rf", @@ -57,6 +71,7 @@ def test_energy(fixture, platform, request): shift_coulomb=str(sampler._shift_coulomb), shift_delta=str(sampler._shift_delta), platform=platform, + map=dyn_map, ) # Loop until we accept an insertion move. From a6992fccc9093638bb09702406b11ddc1eecac2e Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Sun, 29 Mar 2026 10:21:30 +0100 Subject: [PATCH 09/36] Optimise taylor_power branch. --- src/loch/_kernels.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/loch/_kernels.py b/src/loch/_kernels.py index 83c7114..8165486 100644 --- a/src/loch/_kernels.py +++ b/src/loch/_kernels.py @@ -568,8 +568,9 @@ // sig6 = sigma^6 / (alpha^m * sigma^6 + r^6) const float s6_val = powf(s, 6.0f); const float r6 = r * r * r * r * r * r; - const float alpha_m = (sc_taylor_power == 1) - ? a : powf(a, (float)sc_taylor_power); + const float alpha_m = (sc_taylor_power == 1) ? a + : (sc_taylor_power == 0) ? 1.0f + : powf(a, (float)sc_taylor_power); sig6 = s6_val / (alpha_m * s6_val + r6); } else From aa87449fe9edc6ca66022eb2340a12f67e2a0b95 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 31 Mar 2026 14:30:40 +0100 Subject: [PATCH 10/36] Clarify move acceptance probability. [ci skip] --- src/loch/_sampler.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/loch/_sampler.py b/src/loch/_sampler.py index 2fe278f..ead70d3 100644 --- a/src/loch/_sampler.py +++ b/src/loch/_sampler.py @@ -1041,7 +1041,8 @@ def water_state(self) -> _np.ndarray: def move_acceptance_probability(self) -> float: """ - Return the acceptance probability. + Return the acceptance probability. Note that this can be greater than + 1, since multiple insertions/deletions can be accepter per move. Returns ------- From 537571d2f9f05d35415357822b06d0acef57241a Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 1 Apr 2026 11:56:02 +0100 Subject: [PATCH 11/36] Add methods to get and restore sampling statistics. --- src/loch/_sampler.py | 38 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 38 insertions(+) diff --git a/src/loch/_sampler.py b/src/loch/_sampler.py index ead70d3..23698bd 100644 --- a/src/loch/_sampler.py +++ b/src/loch/_sampler.py @@ -1102,6 +1102,44 @@ def reset(self) -> None: self._num_accepted_insertions = 0 self._num_accepted_deletions = 0 + def restore_stats(self, stats: dict) -> None: + """ + Restore sampler statistics from a dictionary. + + Parameters + ---------- + + stats : dict + Dictionary of sampler statistics as returned by ``get_stats()``. + """ + self._num_moves = stats["num_moves"] + self._num_accepted = stats["num_accepted"] + self._num_insertions = stats["num_insertions"] + self._num_deletions = stats["num_deletions"] + self._num_accepted_attempts = stats["num_accepted_attempts"] + self._num_accepted_insertions = stats["num_accepted_insertions"] + self._num_accepted_deletions = stats["num_accepted_deletions"] + + def get_stats(self) -> dict: + """ + Return the current sampler statistics as a dictionary. + + Returns + ------- + + dict + Dictionary of sampler statistics. + """ + return { + "num_moves": self._num_moves, + "num_accepted": self._num_accepted, + "num_insertions": self._num_insertions, + "num_deletions": self._num_deletions, + "num_accepted_attempts": self._num_accepted_attempts, + "num_accepted_insertions": self._num_accepted_insertions, + "num_accepted_deletions": self._num_accepted_deletions, + } + # Clear the forces. self._nonbonded_force = None self._custom_nonbonded_force = None From 40098c39e4bd44bdc61075fc097997457520341f Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 2 Apr 2026 15:24:39 +0100 Subject: [PATCH 12/36] Update CHANGELOG. [ci skip] --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 15c2eaa..7b6c5ae 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,7 +4,7 @@ Changelog [2026.1.0](https://github.com/openbiosim/loch/compare/2025.2.0...2026.1.0) - ******** ------------------------------------------------------------------------------------- -* Please add an item to this CHANGELOG for any new features or bug fixes when creating a PR. +* Add support for getting and restoring sampling statistics. [2025.2.0](https://github.com/openbiosim/loch/compare/2025.1.0...2025.2.0) - Feb 2026 ------------------------------------------------------------------------------------- From 815186685479697ccc0c27796a8b28cc89c41da3 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 16 Apr 2026 11:26:32 +0100 Subject: [PATCH 13/36] Handle XED force field virtual sites. --- CHANGELOG.md | 1 + src/loch/_sampler.py | 135 ++++++++++++++++++++++++++++++++---- tests/test_vsites.py | 160 +++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 283 insertions(+), 13 deletions(-) create mode 100644 tests/test_vsites.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 7b6c5ae..d5cee1e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,7 @@ Changelog ------------------------------------------------------------------------------------- * Add support for getting and restoring sampling statistics. +* Handle XED force field virtual sites. [2025.2.0](https://github.com/openbiosim/loch/compare/2025.1.0...2025.2.0) - Feb 2026 ------------------------------------------------------------------------------------- diff --git a/src/loch/_sampler.py b/src/loch/_sampler.py index 23698bd..b268716 100644 --- a/src/loch/_sampler.py +++ b/src/loch/_sampler.py @@ -566,6 +566,33 @@ def __init__( except Exception as e: raise ValueError(f"Could not prepare the system for GCMC sampling: {e}") + # Compute per-molecule virtual site information. Virtual sites are + # appended after the real atoms of each molecule in the OpenMM system, + # so all subsequent molecules have their OpenMM particle indices shifted + # by the cumulative number of virtual sites in preceding molecules. + ( + self._total_vsites, + self._vsite_atom_offsets, + self._mol_vsite_charges, + ) = self._get_vsite_offsets(self._system) + + if self._total_vsites > 0: + # Offset water oxygen indices from Sire atom indices to OpenMM + # particle indices. + self._water_indices = ( + self._water_indices + self._vsite_atom_offsets[self._water_indices] + ) + + # Apply the same correction to the reference atom indices. + if self._reference is not None: + self._reference_indices = ( + self._reference_indices + + self._vsite_atom_offsets[self._reference_indices] + ) + + # Update the total atom count to include virtual site particles. + self._num_atoms = self._system.num_atoms() + self._total_vsites + # Validate the platform parameter. valid_platforms = {"auto", "cuda", "opencl"} @@ -1102,6 +1129,13 @@ def reset(self) -> None: self._num_accepted_insertions = 0 self._num_accepted_deletions = 0 + # Clear the forces. + self._nonbonded_force = None + self._custom_nonbonded_force = None + + # Clear the OpenMM context. + self._openmm_context = None + def restore_stats(self, stats: dict) -> None: """ Restore sampler statistics from a dictionary. @@ -1140,17 +1174,11 @@ def get_stats(self) -> dict: "num_accepted_deletions": self._num_accepted_deletions, } - # Clear the forces. - self._nonbonded_force = None - self._custom_nonbonded_force = None - - # Clear the OpenMM context. - self._openmm_context = None - def ghost_residues(self) -> _np.ndarray: """ - Return the current indices of the ghost water residues in the OpenMM - context. + Return the residue indices of the current ghost waters in the input + topology. These are Sire/BioSimSpace residue indices and do not + include any virtual site particles that were added on context creation. Returns ------- @@ -1785,6 +1813,75 @@ def _get_box_information(self, space): return cell_matrix, cell_matrix_inverse, M + @staticmethod + def _get_vsite_offsets(system): + """ + Compute per-atom OpenMM index offsets due to virtual sites. + + In OpenMM, virtual site particles are appended after the real atoms + of each molecule. Molecules that appear after a molecule with virtual + sites therefore have their OpenMM particle indices shifted relative + to their Sire atom indices. + + Parameters + ---------- + + system: sire.system.System + The molecular system. + + Returns + ------- + + total_vsites: int + Total number of virtual site particles in the system. + + atom_offsets: numpy.ndarray + Array of shape (num_sire_atoms,) where entry i is the + cumulative number of virtual sites in all molecules that + precede the molecule containing Sire atom i. Adding this + offset to a Sire atom index yields the corresponding OpenMM + particle index. + + mol_vsite_charges: dict + Mapping from molecule number to a list of virtual site charges in + units of elementary charge. Only molecules that carry virtual sites + appear as keys; all other molecules are absent from the dict. + """ + n_sire_atoms = system.num_atoms() + atom_offsets = _np.zeros(n_sire_atoms, dtype=_np.int32) + mol_vsite_charges = {} + total_vsites = 0 + + try: + vsite_mols = system["property n_virtual_sites"].molecules() + except Exception: + # No molecules carry virtual sites. + return 0, atom_offsets, mol_vsite_charges + + all_atoms = system.atoms() + + for mol in vsite_mols: + n_vs = int(mol.property("n_virtual_sites")) + if n_vs <= 0: + continue + + # Locate where this molecule's atoms sit in the global index space, + # then shift every subsequent atom's offset by n_vs in one operation. + mol_start = int(_np.array(all_atoms.find(mol.atoms()))[0]) + mol_end = mol_start + mol.num_atoms() + atom_offsets[mol_end:] += n_vs + total_vsites += n_vs + + try: + raw_charges = mol.property("vs_charges") + vs_charges = [float(raw_charges[k]) for k in range(n_vs)] + except Exception: + vs_charges = [0.0] * n_vs + + mol_vsite_charges[mol.number()] = vs_charges + + return total_vsites, atom_offsets, mol_vsite_charges + @staticmethod def _get_reference_indices(system, reference): """ @@ -1937,6 +2034,10 @@ def _initialise_gpu_memory(self): for q in mol.property("charge"): charges[i] = q.value() i += 1 + # Append virtual site charges (zero LJ, non-zero charge). + for vc in self._mol_vsite_charges.get(mol.number(), []): + charges[i] = vc + i += 1 # Convert to a GPU array. charges = self._backend.to_gpu(charges.astype(_np.float32)) @@ -1954,6 +2055,12 @@ def _initialise_gpu_memory(self): sigmas[i] = lj.sigma().value() epsilons[i] = lj.epsilon().value() i += 1 + # Virtual sites have zero LJ. Use sigma=1.0 Å as a + # nominal placeholder (epsilon=0 so it has no effect). + for _ in self._mol_vsite_charges.get(mol.number(), []): + sigmas[i] = 1.0 + epsilons[i] = 0.0 + i += 1 # Convert to GPU arrays. sigmas = self._backend.to_gpu(sigmas.astype(_np.float32)) @@ -2063,8 +2170,9 @@ def _initialise_gpu_memory(self): # This is a null LJ parameter. if _np.isclose(lj.epsilon().value(), 0.0): - idx = atoms.find(atom) - is_ghost_fep[idx] = 1 + sire_idx = atoms.find(atom) + omm_idx = sire_idx + int(self._vsite_atom_offsets[sire_idx]) + is_ghost_fep[omm_idx] = 1 # The charge at the perturbed state is zero. elif _np.isclose(charge1, 0.0): @@ -2073,8 +2181,9 @@ def _initialise_gpu_memory(self): # This is a null LJ parameter. if _np.isclose(lj.epsilon().value(), 0.0): - idx = atoms.find(atom) - is_ghost_fep[idx] = 1 + sire_idx = atoms.find(atom) + omm_idx = sire_idx + int(self._vsite_atom_offsets[sire_idx]) + is_ghost_fep[omm_idx] = 1 # Convert to GPU array. is_ghost_fep = self._backend.to_gpu(is_ghost_fep.astype(_np.int32)) diff --git a/tests/test_vsites.py b/tests/test_vsites.py new file mode 100644 index 0000000..b28519d --- /dev/null +++ b/tests/test_vsites.py @@ -0,0 +1,160 @@ +import os + +import pytest +import sire as sr + +from loch import GCMCSampler + + +# --------------------------------------------------------------------------- +# Helpers +# --------------------------------------------------------------------------- + + +def _add_vsites_to_mol(mol, n_vs): + """ + Return a copy of mol with n_vs zero-charge virtual sites all parented to + atom 0. + + Using atom 0 for all parent references works for any molecule regardless + of its size. Virtual sites have zero charge so they do not affect energies, + making it possible to compare sampler results with and without virtual sites. + """ + n_atoms = mol.num_atoms() + + # LocalCoordinatesSite requires three parent atoms; clamp to valid range so + # this helper works even for single-atom molecules (e.g. ions). + p1 = min(1, n_atoms - 1) + p2 = min(2, n_atoms - 1) + + vsite_dict = { + str(k): { + "vs_indices": [0, p1, p2], + "vs_ows": [1, 0, 0], + "vs_xs": [1, -1, 0], + "vs_ys": [0, 1, -1], + # Offset each vsite slightly so their positions are distinct. + "vs_local": [(k + 1) * 0.03, 0, 0], + } + for k in range(n_vs) + } + + # All vsites are children of atom 0. + parents = {str(i): [] for i in range(n_atoms)} + parents["0"] = list(range(n_vs)) + + cursor = mol.cursor() + cursor.set("n_virtual_sites", n_vs) + cursor.set("vs_charges", [0.0] * n_vs) + cursor.set("virtual_sites", vsite_dict) + cursor.set("parents", parents) + return cursor.commit() + + +# --------------------------------------------------------------------------- +# Unit tests for _get_vsite_offsets (no GPU required) +# --------------------------------------------------------------------------- + + +def test_get_vsite_offsets_no_vsites(): + """ + _get_vsite_offsets on a system with no virtual sites returns all-zero + offsets and empty per-molecule charge lists. + """ + mols = sr.load_test_files("bpti.prm7", "bpti.rst7") + + total_vsites, offsets, mol_vsite_charges = GCMCSampler._get_vsite_offsets(mols) + + assert total_vsites == 0 + assert offsets.shape == (mols.num_atoms(),) + assert (offsets == 0).all() + assert mol_vsite_charges == {} + + +def test_get_vsite_offsets_with_vsites(): + """ + Adding N virtual sites to molecule 0 gives a zero offset for every atom + inside that molecule and an offset of N for all atoms in subsequent + molecules. + """ + mols = sr.load_test_files("bpti.prm7", "bpti.rst7") + + n_vs = 2 + n_atoms_mol0 = mols[0].num_atoms() + + mols_with_vs = mols.clone() + mols_with_vs.update(_add_vsites_to_mol(mols_with_vs[0], n_vs)) + + total_vsites, offsets, mol_vsite_charges = GCMCSampler._get_vsite_offsets( + mols_with_vs + ) + + assert total_vsites == n_vs + + # Atoms in molecule 0 precede no vsite-bearing molecule, so offset is 0. + assert (offsets[:n_atoms_mol0] == 0).all() + + # All atoms in molecules 1..N follow molecule 0 and are shifted by n_vs. + assert (offsets[n_atoms_mol0:] == n_vs).all() + + # Only molecule 0 appears in the dict, with n_vs zero charges. + assert len(mol_vsite_charges) == 1 + assert list(mol_vsite_charges.values())[0] == [0.0] * n_vs + + +# --------------------------------------------------------------------------- +# Integration tests (GPU required) +# --------------------------------------------------------------------------- + + +@pytest.mark.skipif( + "CUDA_VISIBLE_DEVICES" not in os.environ, + reason="Requires CUDA enabled GPU.", +) +@pytest.mark.parametrize("platform", ["cuda", "opencl"]) +def test_vsite_index_offsets(bpti, platform): + """ + GCMCSampler correctly offsets _water_indices and _reference_indices when + virtual sites are present on a preceding molecule. + + In BPTI the protein is molecule 0 and all water molecules follow it. + Adding N vsites to the protein shifts every water's OpenMM particle index + by N. The reference atoms also live in the protein, so their OpenMM + indices are not shifted (no vsites precede molecule 0). + """ + mols, reference = bpti + + common_kwargs = dict( + reference=reference, + cutoff_type="rf", + cutoff="10 A", + ghost_file=None, + log_file=None, + test=True, + platform=platform, + seed=42, + ) + + # Baseline: no virtual sites. + baseline = GCMCSampler(mols, **common_kwargs) + + # Add 2 zero-charge virtual sites to molecule 0 (the protein). All + # water molecules follow the protein, so their OpenMM indices shift by 2. + n_vs = 2 + mols_with_vs = mols.clone() + mols_with_vs.update(_add_vsites_to_mol(mols_with_vs[0], n_vs)) + + sampler = GCMCSampler(mols_with_vs, **common_kwargs) + + # Total vsite count must match what we added. + assert sampler._total_vsites == n_vs + + # _num_atoms must include the virtual site particles. + assert sampler._num_atoms == baseline._num_atoms + n_vs + + # Every water oxygen index is shifted by n_vs (waters follow the protein). + assert (sampler._water_indices == baseline._water_indices + n_vs).all() + + # Reference atoms are inside molecule 0 (the protein), which has no + # preceding vsites, so their OpenMM indices are unchanged. + assert (sampler._reference_indices == baseline._reference_indices).all() From d90aae6ee0b625a945769aefc3465dc146e68f31 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 21 Apr 2026 18:20:00 +0100 Subject: [PATCH 14/36] Fix vsite index offsets for OpenMM context vs Sire topology operations. --- src/loch/_sampler.py | 11 ++++++++-- tests/test_vsites.py | 50 ++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 59 insertions(+), 2 deletions(-) diff --git a/src/loch/_sampler.py b/src/loch/_sampler.py index b268716..ea2b146 100644 --- a/src/loch/_sampler.py +++ b/src/loch/_sampler.py @@ -576,6 +576,12 @@ def __init__( self._mol_vsite_charges, ) = self._get_vsite_offsets(self._system) + # Keep a copy of the Sire atom indices before applying the OpenMM + # offset. These are needed by any operation that works on the Sire + # topology (e.g. _flag_ghost_waters, ghost_residues), which has no + # knowledge of virtual site particles. + self._water_indices_sire = self._water_indices.copy() + if self._total_vsites > 0: # Offset water oxygen indices from Sire atom indices to OpenMM # particle indices. @@ -2925,8 +2931,9 @@ def _flag_ghost_waters(self, system): if not isinstance(system, _sr.system.System): raise ValueError("'system' must be a Sire system") - # Now extract the oxygen indices using cached ghost water indices. - ghost_oxygens = self._water_indices[self._get_ghost_waters()] + # Use the Sire atom indices (no vsite offset) so that lookups into the + # input topology are correct regardless of virtual sites in the context. + ghost_oxygens = self._water_indices_sire[self._get_ghost_waters()] # Loop over the ghost waters and set the is_ghost property. for i in ghost_oxygens: diff --git a/tests/test_vsites.py b/tests/test_vsites.py index b28519d..225b491 100644 --- a/tests/test_vsites.py +++ b/tests/test_vsites.py @@ -158,3 +158,53 @@ def test_vsite_index_offsets(bpti, platform): # Reference atoms are inside molecule 0 (the protein), which has no # preceding vsites, so their OpenMM indices are unchanged. assert (sampler._reference_indices == baseline._reference_indices).all() + + +@pytest.mark.skipif( + "CUDA_VISIBLE_DEVICES" not in os.environ, + reason="Requires CUDA enabled GPU.", +) +@pytest.mark.parametrize("platform", ["cuda", "opencl"]) +def test_flag_ghost_waters_sire_indices(bpti, platform): + """ + _flag_ghost_waters must use Sire atom indices, not OpenMM particle indices. + + With virtual sites present, OpenMM indices are larger than the corresponding + Sire atom indices. Using the OpenMM indices to look up atoms in the Sire + topology raises SireError::invalid_index. This test confirms that the + Sire-side indices are used and that _flag_ghost_waters completes without + error. + """ + mols, reference = bpti + + n_vs = 2 + mols_with_vs = mols.clone() + mols_with_vs.update(_add_vsites_to_mol(mols_with_vs[0], n_vs)) + + sampler = GCMCSampler( + mols_with_vs, + reference=reference, + cutoff_type="rf", + cutoff="10 A", + ghost_file=None, + log_file=None, + test=True, + platform=platform, + seed=42, + ) + + n_sire_atoms = sampler.system().num_atoms() + + # Sire indices must be within the topology's atom range. + assert (sampler._water_indices_sire < n_sire_atoms).all() + + # OpenMM indices exceed the Sire range by n_vs, which is exactly the + # bug that would have caused SireError::invalid_index. + assert (sampler._water_indices == sampler._water_indices_sire + n_vs).all() + + # _flag_ghost_waters must complete without raising an index error. + flagged = sampler._flag_ghost_waters(sampler.system()) + + # Every initially-ghost water molecule should be flagged. + n_flagged = len(flagged["property is_ghost_water"].molecules()) + assert n_flagged == sampler._num_ghost_waters From 6cc74f04ea2ea538091e1d74410fae4aa267b20a Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 21 Apr 2026 19:30:21 +0100 Subject: [PATCH 15/36] Fix link. [ci skip] --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index ac70f87..e6691e0 100644 --- a/README.md +++ b/README.md @@ -384,5 +384,5 @@ inpcrd = AmberInpcrdFile("gcmc_system.inpcrd") * We thank the [Essex Lab](https://essexgroup.soton.ac.uk/) and [grand](https://github.com/essex-lab/grand) for the inspiration. -* Many thanks to [Gregory Ross](https://github.com/GregRRoss) for clarifying +* Many thanks to [Gregory Ross](https://github.com/gregoryross) for clarifying the parallelisation scheme described [here](https://doi.org/10.1021/acs.jctc.0c00660). From 0940cf312aadc07cf89e4c94907c6c7b5b6e5122 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Fri, 24 Apr 2026 18:34:32 +0100 Subject: [PATCH 16/36] Allow optional context to be passed to num_waters. --- src/loch/_sampler.py | 29 ++++++++++++++++++++--------- 1 file changed, 20 insertions(+), 9 deletions(-) diff --git a/src/loch/_sampler.py b/src/loch/_sampler.py index ea2b146..0b9e04e 100644 --- a/src/loch/_sampler.py +++ b/src/loch/_sampler.py @@ -979,7 +979,7 @@ def delete_waters(self, context: _openmm.Context) -> None: # Set the number of waters in the GCMC sphere to zero. self._N = 0 - def num_waters(self) -> int: + def num_waters(self, context=None) -> int: """ Return the number of waters in the GCMC region. @@ -988,18 +988,29 @@ def num_waters(self) -> int: num_waters: int The number of waters. + + context: openmm.Context, optional + The OpenMM context to use for counting the waters. If None, then the + internal context will be used if available. """ - # The last move was a bulk sampling move, so we need to recalculate - # the number of waters in the GCMC sphere. - if self._reference is not None and self._is_bulk: - if not self._openmm_context: - msg = "OpenMM context is not set!" - _logger.error(msg) - raise RuntimeError(msg) + # Whether we need to recalculate the number of waters in the GCMC sphere. + recalculate = context is not None or ( + self._reference is not None and self._is_bulk + ) + + # We need to recalculate the number of waters. + if recalculate: + if context is None: + if not self._openmm_context: + msg = "OpenMM context is not set!" + _logger.error(msg) + raise RuntimeError(msg) + else: + context = self._openmm_context # Get the OpenMM state. - state = self._openmm_context.getState(getPositions=True) + state = context.getState(getPositions=True) # Get the current positions in Angstrom. positions = state.getPositions(asNumpy=True) / _openmm.unit.angstrom From 40f79b3becd613f21be4205424bd536935e17507 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 29 Apr 2026 14:48:43 +0100 Subject: [PATCH 17/36] Add support for long-range LJ dispersion correction. --- src/loch/_sampler.py | 183 +++++++++++++++++++++++++++++++++++-------- 1 file changed, 150 insertions(+), 33 deletions(-) diff --git a/src/loch/_sampler.py b/src/loch/_sampler.py index 0b9e04e..992b3f9 100644 --- a/src/loch/_sampler.py +++ b/src/loch/_sampler.py @@ -287,6 +287,11 @@ def __init__( else: self._is_pme = False + # LRC state: initialised lazily on first move() call. + self._has_gcmc_lrc = False + self._lrc_w_solute = 0.0 + self._lrc_ww_half = 0.0 + try: self._radius = self._validate_sire_unit("radius", radius, _sr.u("A")) except Exception as e: @@ -707,7 +712,8 @@ def __init__( # Null the nonbonded forces. self._nonbonded_force = None - self._custom_nonbonded_force = None + self._custom_coulomb_force = None + self._custom_lj_force = None # Flag for whether the last move was a bulk sampling move. self._is_bulk = False @@ -1148,7 +1154,8 @@ def reset(self) -> None: # Clear the forces. self._nonbonded_force = None - self._custom_nonbonded_force = None + self._custom_coulomb_force = None + self._custom_lj_force = None # Clear the OpenMM context. self._openmm_context = None @@ -1273,6 +1280,10 @@ def move(self, context: _openmm.Context) -> list[int]: # We only need to get the positions and initial energy for the first # batch. These will be updated dynamically as moves are accepted. if num_batches == 1: + # Detect GCMC LRC parameters from context on first call. + if not self._has_gcmc_lrc: + self._init_gcmc_lrc(context) + # Get the OpenMM state. state = context.getState(getPositions=True, getEnergy=self._is_pme) @@ -1286,6 +1297,11 @@ def move(self, context: _openmm.Context) -> list[int]: else: initial_energy = None + # Cache the box volume (NVT, so constant throughout). + if self._has_gcmc_lrc: + box = state.getPeriodicBoxVectors(asNumpy=True) + v_nm3 = _np.linalg.det(box / _openmm.unit.nanometer) + # Sample within the GCMC sphere. if self._reference is not None and not self._is_bulk: target = self._get_target_position(positions_angstrom).astype( @@ -1521,6 +1537,10 @@ def move(self, context: _openmm.Context) -> list[int]: # Insertion move. if is_deletion[idx] == 0: + # Capture n_w before the insertion for LRC delta calculation. + if self._has_gcmc_lrc: + n_w_before_insert = context.getParameter("n_w") + # Accept the move. self._accept_insertion( idx, idx_water, positions_openmm, positions_angstrom, context @@ -1548,6 +1568,19 @@ def move(self, context: _openmm.Context) -> list[int]: getEnergy=True ).getPotentialEnergy() + # Add the analytic LRC delta so the PME correction sees only + # the RF→PME electrostatic difference. + if self._has_gcmc_lrc: + dLRC = ( + ( + self._lrc_w_solute + + 2.0 * n_w_before_insert * self._lrc_ww_half + ) + / v_nm3 + * _openmm.unit.kilojoules_per_mole + ) + dE_RF += dLRC + # Compute the PME acceptance correction. acc_prob = _np.exp( min( @@ -1608,6 +1641,10 @@ def move(self, context: _openmm.Context) -> list[int]: # Deletion move. else: + # Capture n_w before the deletion for LRC delta calculation. + if self._has_gcmc_lrc: + n_w_before_delete = context.getParameter("n_w") + # Accept the move. self._accept_deletion(candidates[idx], context) @@ -1633,6 +1670,20 @@ def move(self, context: _openmm.Context) -> list[int]: getEnergy=True ).getPotentialEnergy() + # Add the analytic LRC delta. + if self._has_gcmc_lrc: + dLRC = ( + -( + self._lrc_w_solute + + 2.0 + * (n_w_before_delete - 1.0) + * self._lrc_ww_half + ) + / v_nm3 + * _openmm.unit.kilojoules_per_mole + ) + dE_RF += dLRC + # Compute the PME acceptance correction. acc_prob = _np.exp( min( @@ -2125,31 +2176,37 @@ def _initialise_gpu_memory(self): ) # Flags for the required force. - has_gng = False + has_gng_coul = False + has_gng_lj = False # Find the required forces. for force in d.context().getSystem().getForces(): - if force.getName() == "GhostNonGhostNonbondedForce": - gng_force = force - has_gng = True - break + if force.getName() == "GhostNonGhostCoulombForce": + gng_coul_force = force + has_gng_coul = True + elif force.getName() == "GhostNonGhostLJForce": + gng_lj_force = force + has_gng_lj = True # Make sure the force was found. - if not has_gng: + if not has_gng_coul: raise ValueError( - "Could not find the GhostNonGhostNonbondedForce in the system" + "Could not find the GhostNonGhostCoulombForce in the system" + ) + if not has_gng_lj: + raise ValueError( + "Could not find the GhostNonGhostLJForce in the system" ) - # Get the parameters for the GhostNonGhostNonbondedForce. + # Get the parameters for the GhostNonGhost nonbonded forces. charges = _np.zeros(self._num_atoms, dtype=_np.float32) sigmas = _np.zeros(self._num_atoms, dtype=_np.float32) epsilons = _np.zeros(self._num_atoms, dtype=_np.float32) alphas = _np.zeros(self._num_atoms, dtype=_np.float32) - for i in range(gng_force.getNumParticles()): + for i in range(gng_coul_force.getNumParticles()): # Custom force parameters are returned as floats. - q, half_sigma, two_sqrt_epsilon, alpha, _ = ( - gng_force.getParticleParameters(i) - ) + q, alpha, _ = gng_coul_force.getParticleParameters(i) + half_sigma, two_sqrt_epsilon, _ = gng_lj_force.getParticleParameters(i) # Charge in |e|, sigma in nm, epsilon in kJ/mol. charges[i] = q # Rescale and convert units. @@ -2321,6 +2378,15 @@ def _initialise_gpu_memory(self): (1, self._num_waters), _np.int32 ) + def _init_gcmc_lrc(self, context): + """Detect and cache GCMC LRC parameters from the OpenMM context.""" + try: + self._lrc_w_solute = context.getParameter("lrc_w_solute") + self._lrc_ww_half = context.getParameter("lrc_ww_half") + self._has_gcmc_lrc = True + except Exception: + self._has_gcmc_lrc = False + def _accept_insertion( self, idx, idx_water, positions_openmm, positions_angstrom, context ): @@ -2374,14 +2440,20 @@ def _accept_insertion( ) # Update the custom NonBondedForce parameters. if self._is_fep: - self._custom_nonbonded_force.setParticleParameters( + self._custom_coulomb_force.setParticleParameters( start_idx + i, ( self._water_charge[i], + 0.0, + 0.0, + ), + ) + self._custom_lj_force.setParticleParameters( + start_idx + i, + ( self._water_sigma_custom[i], self._water_epsilon_custom[i], 0.0, - 0.0, ), ) @@ -2393,7 +2465,8 @@ def _accept_insertion( # Update the CustomNonbondedForce parameters in the context. if self._is_fep: - self._custom_nonbonded_force.updateParametersInContext(context) + self._custom_coulomb_force.updateParametersInContext(context) + self._custom_lj_force.updateParametersInContext(context) # Update the state of the water on the GPU. self._kernels["update_water"]( @@ -2417,6 +2490,10 @@ def _accept_insertion( # Update the number of waters in the sampling volume. self._N += 1 + # Update the GCMC LRC water count in the context. + if self._has_gcmc_lrc: + context.setParameter("n_w", context.getParameter("n_w") + 1.0) + def _accept_deletion(self, idx, context): """ Accept a deletion move. @@ -2445,13 +2522,19 @@ def _accept_deletion(self, idx, context): ) # Update the CustomNonBondedForce parameters. if self._is_fep: - self._custom_nonbonded_force.setParticleParameters( + self._custom_coulomb_force.setParticleParameters( start_idx + i, ( 0.0, - self._water_sigma_custom[i], 0.0, 0.0, + ), + ) + self._custom_lj_force.setParticleParameters( + start_idx + i, + ( + self._water_sigma_custom[i], + 0.0, 0.0, ), ) @@ -2461,7 +2544,8 @@ def _accept_deletion(self, idx, context): # Update the CustomNonbondedForce parameters in the context. if self._is_fep: - self._custom_nonbonded_force.updateParametersInContext(context) + self._custom_coulomb_force.updateParametersInContext(context) + self._custom_lj_force.updateParametersInContext(context) # Update the state of the water on the GPU. self._kernels["update_water"]( @@ -2487,6 +2571,10 @@ def _accept_deletion(self, idx, context): # Update the number of waters in the sampling volume. self._N -= 1 + # Update the GCMC LRC water count in the context. + if self._has_gcmc_lrc: + context.setParameter("n_w", context.getParameter("n_w") - 1.0) + def _reject_deletion(self, idx, context): """ Reject a deletion move. @@ -2518,14 +2606,20 @@ def _reject_deletion(self, idx, context): ) # Update the CustomNonBondedForce parameters. if self._is_fep: - self._custom_nonbonded_force.setParticleParameters( + self._custom_coulomb_force.setParticleParameters( start_idx + i, ( self._water_charge[i], + 0.0, + 0.0, + ), + ) + self._custom_lj_force.setParticleParameters( + start_idx + i, + ( self._water_sigma_custom[i], self._water_epsilon_custom[i], 0.0, - 0.0, ), ) @@ -2534,7 +2628,8 @@ def _reject_deletion(self, idx, context): # Update the CustomNonbondedForce parameters in the context. if self._is_fep: - self._custom_nonbonded_force.updateParametersInContext(context) + self._custom_coulomb_force.updateParametersInContext(context) + self._custom_lj_force.updateParametersInContext(context) # Update the state of the water on the GPU. self._kernels["update_water"]( @@ -2560,6 +2655,10 @@ def _reject_deletion(self, idx, context): # Update the number of waters in the sampling volume. self._N += 1 + # Update the GCMC LRC water count in the context. + if self._has_gcmc_lrc: + context.setParameter("n_w", context.getParameter("n_w") + 1.0) + def _set_water_state(self, context, indices=None, states=None, force=False): """ Update the state for a list of waters. This can be used by external @@ -2597,7 +2696,8 @@ def _set_water_state(self, context, indices=None, states=None, force=False): # Assume the context has been recreated, so we need to get the # new forces. self._nonbonded_force = None - self._custom_nonbonded_force = None + self._custom_coulomb_force = None + self._custom_lj_force = None # Update even if the state is unchanged. force = True @@ -2627,13 +2727,19 @@ def _set_water_state(self, context, indices=None, states=None, force=False): ) # Update the CustomNonbondedForce parameters. if self._is_fep: - self._custom_nonbonded_force.setParticleParameters( + self._custom_coulomb_force.setParticleParameters( start_idx + i, ( 0.0, - self._water_sigma_custom[i], 0.0, 0.0, + ), + ) + self._custom_lj_force.setParticleParameters( + start_idx + i, + ( + self._water_sigma_custom[i], + 0.0, 0.0, ), ) @@ -2674,14 +2780,20 @@ def _set_water_state(self, context, indices=None, states=None, force=False): ) # Update the CustomNonBondedForce parameters. if self._is_fep: - self._custom_nonbonded_force.setParticleParameters( + self._custom_coulomb_force.setParticleParameters( start_idx + i, ( self._water_charge[i], + 0.0, + 0.0, + ), + ) + self._custom_lj_force.setParticleParameters( + start_idx + i, + ( self._water_sigma_custom[i], self._water_epsilon_custom[i], 0.0, - 0.0, ), ) @@ -2717,7 +2829,8 @@ def _set_water_state(self, context, indices=None, states=None, force=False): # Update the CustomNonbondedForce parameters in the context. if self._is_fep: - self._custom_nonbonded_force.updateParametersInContext(context) + self._custom_coulomb_force.updateParametersInContext(context) + self._custom_lj_force.updateParametersInContext(context) def _set_nonbonded_forces(self, context): """ @@ -2730,13 +2843,15 @@ def _set_nonbonded_forces(self, context): The OpenMM context to use. """ if self._nonbonded_force is None or ( - self._is_fep and self._custom_nonbonded_force is None + self._is_fep and self._custom_coulomb_force is None ): for force in context.getSystem().getForces(): if isinstance(force, _openmm.NonbondedForce): self._nonbonded_force = force - elif self._is_fep and force.getName() == "GhostNonGhostNonbondedForce": - self._custom_nonbonded_force = force + elif self._is_fep and force.getName() == "GhostNonGhostCoulombForce": + self._custom_coulomb_force = force + elif self._is_fep and force.getName() == "GhostNonGhostLJForce": + self._custom_lj_force = force elif "Barostat" in force.getName(): msg = ( f"GCMC must be used at constant volume: " @@ -2750,7 +2865,9 @@ def _set_nonbonded_forces(self, context): _logger.error(msg) raise ValueError(msg) - if self._is_fep and self._custom_nonbonded_force is None: + if self._is_fep and ( + self._custom_coulomb_force is None or self._custom_lj_force is None + ): msg = "Could not find a CustomNonbondedForce in the system" _logger.error(msg) raise ValueError(msg) From 23bb1cd70dcd68402a17a2d563729cc5f8bfe532 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 30 Apr 2026 10:15:51 +0100 Subject: [PATCH 18/36] Add Beutler soft-core form for ABFE. --- src/loch/_kernels.py | 30 +++++++++++++++--------------- src/loch/_sampler.py | 33 +++++++++++++++++++++------------ src/loch/_softcore.py | 1 + tests/test_energy.py | 9 ++++----- 4 files changed, 41 insertions(+), 32 deletions(-) diff --git a/src/loch/_kernels.py b/src/loch/_kernels.py index 8165486..d402b89 100644 --- a/src/loch/_kernels.py +++ b/src/loch/_kernels.py @@ -395,10 +395,10 @@ float rf_kappa, float rf_correction, int softcore_form, - float sc_coulomb_power, float sc_shift_coulomb, float sc_shift_delta, - int sc_taylor_power) + int sc_taylor_power, + float sc_beutler_alpha) { // Work out the atom index. const int idx_atom = GET_GLOBAL_ID(0); @@ -562,6 +562,7 @@ // Compute the LJ interaction using the chosen soft-core form. float sig6; + float lj_prefactor = 1.0f; if (softcore_form == 1) { // Taylor soft-core LJ: @@ -573,6 +574,16 @@ : powf(a, (float)sc_taylor_power); sig6 = s6_val / (alpha_m * s6_val + r6); } + else if (softcore_form == 2) + { + // Beutler soft-core LJ: + // sig6 = sigma^6 / (sc_beutler_alpha * sigma^6 * alpha + r^6) + // V_LJ = (1 - alpha) * 4 * epsilon * sig6 * (sig6 - 1) + const float s6_val = powf(s, 6.0f); + const float r6 = r * r * r * r * r * r; + sig6 = s6_val / (sc_beutler_alpha * s6_val * a + r6); + lj_prefactor = 1.0f - a; + } else { // Zacharias soft-core LJ: @@ -581,22 +592,11 @@ const float delta_lj = sc_shift_delta * a; sig6 = powf(s, 6.0f) / powf((s * delta_lj) + (r * r), 3.0f); } - energy_lj[idx] += 4.0f * e * sig6 * (sig6 - 1.0f); - - // Compute the Coulomb power expression. - float cpe; - if (sc_coulomb_power == 0.0f) - { - cpe = 1.0f; - } - else - { - cpe = powf((1.0f - a), sc_coulomb_power); - } + energy_lj[idx] += lj_prefactor * 4.0f * e * sig6 * (sig6 - 1.0f); // Compute the Coulomb interaction. energy_coul[idx] += (q0 * q1) * - ((cpe / sqrtf((sc_shift_coulomb * sc_shift_coulomb * a) + ((1.0f / sqrtf((sc_shift_coulomb * sc_shift_coulomb * a) + (r * r))) + (rf_kappa * r2) - rf_correction); } diff --git a/src/loch/_sampler.py b/src/loch/_sampler.py index 992b3f9..d2e7ac0 100644 --- a/src/loch/_sampler.py +++ b/src/loch/_sampler.py @@ -79,11 +79,11 @@ def __init__( lambda_value: float = 0.0, rest2_scale: float = 1.0, rest2_selection: _Optional[str] = None, - coulomb_power: float = 0.0, shift_coulomb: str = "1 A", shift_delta: str = "1.5 A", softcore_form: str = "zacharias", taylor_power: int = 1, + beutler_alpha: float = 0.5, swap_end_states: bool = False, restart: bool = False, overwrite: bool = False, @@ -212,8 +212,8 @@ def __init__( softcore_form : str The soft-core potential form to use for alchemical interactions. - This can be either 'zacharias' or 'taylor'. The default is - 'zacharias'. + Valid options are 'zacharias' (default), 'taylor', and 'beutler'. + The Beutler form is recommended for ABFE calculations. taylor_power : int The power to use for the alpha term in the Taylor soft-core LJ @@ -221,6 +221,11 @@ def __init__( Must be between 0 and 4. The default is 1. Only used when softcore_form is 'taylor'. + beutler_alpha : float + The dimensionless scale factor for the r^6 shift in the Beutler + soft-core form. Must be >= 0. The default is 0.5. Only used when + softcore_form is 'beutler'. + swap_end_states: bool Whether to swap the end states of the alchemical systems. @@ -494,12 +499,6 @@ def __init__( ) self._rest2_selection = rest2_selection - try: - coulomb_power = float(coulomb_power) - except Exception: - raise ValueError("'coulomb_power' must be of type 'float'") - self._coulomb_power = float(coulomb_power) - try: self._shift_coulomb = self._validate_sire_unit( "shift_coulomb", shift_coulomb, _sr.u("A") @@ -534,6 +533,14 @@ def __init__( raise ValueError("'taylor_power' must be between 0 and 4") self._taylor_power = taylor_power + try: + beutler_alpha = float(beutler_alpha) + except Exception: + raise ValueError("'beutler_alpha' must be of type 'float'") + if beutler_alpha < 0.0: + raise ValueError("'beutler_alpha' must be >= 0") + self._beutler_alpha = beutler_alpha + if not isinstance(swap_end_states, bool): raise ValueError("'swap_end_states' must be of type 'bool'") self._swap_end_states = swap_end_states @@ -781,7 +788,6 @@ def __str__(self) -> str: f"restart={self._restart}, " f"rest2_scale={self._rest2_scale}, " f"rest2_selection={self._rest2_selection}, " - f"coulomb_power={self._coulomb_power}, " f"shift_coulomb={self._shift_coulomb}, " f"shift_delta={self._shift_delta}, " f"overwrite={self._overwrite}, " @@ -1458,10 +1464,10 @@ def move(self, context: _openmm.Context) -> list[int]: self._rf_kappa, self._rf_correction, self._sc_softcore_form, - self._sc_coulomb_power, self._sc_shift_coulomb, self._sc_shift_delta, self._sc_taylor_power, + self._sc_beutler_alpha, block=(self._num_threads, 1, 1), grid=(self._atom_blocks, self._batch_size, 1), ) @@ -2157,6 +2163,9 @@ def _initialise_gpu_memory(self): if self._softcore_form == _SoftcoreForm.TAYLOR: _map["use_taylor_softening"] = True _map["taylor_power"] = self._taylor_power + elif self._softcore_form == _SoftcoreForm.BEUTLER: + _map["use_beutler_softening"] = True + _map["beutler_alpha"] = self._beutler_alpha # Create a dynamics object. d = mols.dynamics( @@ -2325,10 +2334,10 @@ def _initialise_gpu_memory(self): # Store soft-core parameters as scalars. self._sc_softcore_form = _np.int32(int(self._softcore_form)) - self._sc_coulomb_power = _np.float32(self._coulomb_power) self._sc_shift_coulomb = _np.float32(self._shift_coulomb.value()) self._sc_shift_delta = _np.float32(self._shift_delta.value()) self._sc_taylor_power = _np.int32(self._taylor_power) + self._sc_beutler_alpha = _np.float32(self._beutler_alpha) # Store immutable per-atom buffers on GPU. self._gpu_sigma = sigmas diff --git a/src/loch/_softcore.py b/src/loch/_softcore.py index 6d165af..40af129 100644 --- a/src/loch/_softcore.py +++ b/src/loch/_softcore.py @@ -29,3 +29,4 @@ class SoftcoreForm(IntEnum): ZACHARIAS = 0 TAYLOR = 1 + BEUTLER = 2 diff --git a/tests/test_energy.py b/tests/test_energy.py index 9235cca..5493069 100644 --- a/tests/test_energy.py +++ b/tests/test_energy.py @@ -20,6 +20,7 @@ ("bpti", "zacharias"), ("sd12", "zacharias"), ("sd12", "taylor"), + ("sd12", "beutler"), ], ) def test_energy(fixture, softcore_form, platform, request): @@ -56,6 +57,9 @@ def test_energy(fixture, softcore_form, platform, request): dyn_map = {} if sampler._softcore_form == SoftcoreForm.TAYLOR: dyn_map["use_taylor_softening"] = True + elif sampler._softcore_form == SoftcoreForm.BEUTLER: + dyn_map["use_beutler_softening"] = True + dyn_map["beutler_alpha"] = sampler._beutler_alpha # Create a dynamics object using the modified GCMC system. d = sampler.system().dynamics( @@ -67,7 +71,6 @@ def test_energy(fixture, softcore_form, platform, request): timestep="2 fs", schedule=schedule, lambda_value=lambda_value, - coulomb_power=sampler._coulomb_power, shift_coulomb=str(sampler._shift_coulomb), shift_delta=str(sampler._shift_delta), platform=platform, @@ -227,7 +230,6 @@ def test_platform_consistency(fixture, request): timestep="2 fs", schedule=schedule, lambda_value=lambda_value, - coulomb_power=cuda_sampler._coulomb_power, shift_coulomb=str(cuda_sampler._shift_coulomb), shift_delta=str(cuda_sampler._shift_delta), platform="cuda", @@ -242,7 +244,6 @@ def test_platform_consistency(fixture, request): timestep="2 fs", schedule=schedule, lambda_value=lambda_value, - coulomb_power=opencl_sampler._coulomb_power, shift_coulomb=str(opencl_sampler._shift_coulomb), shift_delta=str(opencl_sampler._shift_delta), platform="opencl", @@ -347,7 +348,6 @@ def test_energy_regression(fixture, platform, request): timestep="2 fs", schedule=schedule, lambda_value=lambda_value, - coulomb_power=sampler._coulomb_power, shift_coulomb=str(sampler._shift_coulomb), shift_delta=str(sampler._shift_delta), platform=platform, @@ -414,7 +414,6 @@ def _create_and_run(seed): timestep="2 fs", schedule=schedule, lambda_value=0.5, - coulomb_power=sampler._coulomb_power, shift_coulomb=str(sampler._shift_coulomb), shift_delta=str(sampler._shift_delta), platform=platform, From c16f401c87af1161db538f66ff6819629cf4deb4 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Fri, 1 May 2026 09:56:06 +0100 Subject: [PATCH 19/36] Optimise softcore. --- src/loch/_kernels.py | 22 +++++++++------------- 1 file changed, 9 insertions(+), 13 deletions(-) diff --git a/src/loch/_kernels.py b/src/loch/_kernels.py index d402b89..146accc 100644 --- a/src/loch/_kernels.py +++ b/src/loch/_kernels.py @@ -551,14 +551,13 @@ const float e = sqrtf(e0 * e1); const float a = alpha[idx_atom]; - // Compute the distance between the atoms. - float r = sqrtf(r2); + // Clamp r2 to avoid singularities. + const float r2_sc = (r2 < 1e-6f) ? 1e-6f : r2; - // Truncate the distance. - if (r < 0.001f) - { - r = 0.001f; - } + // Precompute r^6 and sigma^6 using r2 directly (avoids sqrtf and powf). + const float r6 = r2_sc * r2_sc * r2_sc; + const float s2 = s * s; + const float s6_val = s2 * s2 * s2; // Compute the LJ interaction using the chosen soft-core form. float sig6; @@ -567,8 +566,6 @@ { // Taylor soft-core LJ: // sig6 = sigma^6 / (alpha^m * sigma^6 + r^6) - const float s6_val = powf(s, 6.0f); - const float r6 = r * r * r * r * r * r; const float alpha_m = (sc_taylor_power == 1) ? a : (sc_taylor_power == 0) ? 1.0f : powf(a, (float)sc_taylor_power); @@ -579,8 +576,6 @@ // Beutler soft-core LJ: // sig6 = sigma^6 / (sc_beutler_alpha * sigma^6 * alpha + r^6) // V_LJ = (1 - alpha) * 4 * epsilon * sig6 * (sig6 - 1) - const float s6_val = powf(s, 6.0f); - const float r6 = r * r * r * r * r * r; sig6 = s6_val / (sc_beutler_alpha * s6_val * a + r6); lj_prefactor = 1.0f - a; } @@ -590,14 +585,15 @@ // sig6 = sigma^6 / (sigma*delta + r^2)^3 // delta = shift_delta * alpha const float delta_lj = sc_shift_delta * a; - sig6 = powf(s, 6.0f) / powf((s * delta_lj) + (r * r), 3.0f); + const float denom = (s * delta_lj) + r2_sc; + sig6 = s6_val / (denom * denom * denom); } energy_lj[idx] += lj_prefactor * 4.0f * e * sig6 * (sig6 - 1.0f); // Compute the Coulomb interaction. energy_coul[idx] += (q0 * q1) * ((1.0f / sqrtf((sc_shift_coulomb * sc_shift_coulomb * a) - + (r * r))) + (rf_kappa * r2) - rf_correction); + + r2_sc)) + (rf_kappa * r2) - rf_correction); } } From 474a24c0301b8739f5f8119cf6ac8e4fb560bb53 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Fri, 1 May 2026 09:58:42 +0100 Subject: [PATCH 20/36] Optimise softcore. --- src/loch/_kernels.py | 20 +++++++++----------- 1 file changed, 9 insertions(+), 11 deletions(-) diff --git a/src/loch/_kernels.py b/src/loch/_kernels.py index 8165486..60abc96 100644 --- a/src/loch/_kernels.py +++ b/src/loch/_kernels.py @@ -551,14 +551,13 @@ const float e = sqrtf(e0 * e1); const float a = alpha[idx_atom]; - // Compute the distance between the atoms. - float r = sqrtf(r2); + // Clamp r2 to avoid singularities. + const float r2_sc = (r2 < 1e-6f) ? 1e-6f : r2; - // Truncate the distance. - if (r < 0.001f) - { - r = 0.001f; - } + // Precompute r^6 and sigma^6 using r2 directly (avoids sqrtf and powf). + const float r6 = r2_sc * r2_sc * r2_sc; + const float s2 = s * s; + const float s6_val = s2 * s2 * s2; // Compute the LJ interaction using the chosen soft-core form. float sig6; @@ -566,8 +565,6 @@ { // Taylor soft-core LJ: // sig6 = sigma^6 / (alpha^m * sigma^6 + r^6) - const float s6_val = powf(s, 6.0f); - const float r6 = r * r * r * r * r * r; const float alpha_m = (sc_taylor_power == 1) ? a : (sc_taylor_power == 0) ? 1.0f : powf(a, (float)sc_taylor_power); @@ -579,7 +576,8 @@ // sig6 = sigma^6 / (sigma*delta + r^2)^3 // delta = shift_delta * alpha const float delta_lj = sc_shift_delta * a; - sig6 = powf(s, 6.0f) / powf((s * delta_lj) + (r * r), 3.0f); + const float denom = (s * delta_lj) + r2_sc; + sig6 = s6_val / (denom * denom * denom); } energy_lj[idx] += 4.0f * e * sig6 * (sig6 - 1.0f); @@ -597,7 +595,7 @@ // Compute the Coulomb interaction. energy_coul[idx] += (q0 * q1) * ((cpe / sqrtf((sc_shift_coulomb * sc_shift_coulomb * a) - + (r * r))) + (rf_kappa * r2) - rf_correction); + + r2_sc)) + (rf_kappa * r2) - rf_correction); } } From 0845bff4e29b4ca663e2cf8a9941b4f438fb7cb7 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 6 May 2026 15:02:51 +0100 Subject: [PATCH 21/36] Recombine ghost nonbonded forces for performance. --- src/loch/_sampler.py | 115 +++++++++++++------------------------------ 1 file changed, 34 insertions(+), 81 deletions(-) diff --git a/src/loch/_sampler.py b/src/loch/_sampler.py index d2e7ac0..d367194 100644 --- a/src/loch/_sampler.py +++ b/src/loch/_sampler.py @@ -719,8 +719,7 @@ def __init__( # Null the nonbonded forces. self._nonbonded_force = None - self._custom_coulomb_force = None - self._custom_lj_force = None + self._custom_nonbonded_force = None # Flag for whether the last move was a bulk sampling move. self._is_bulk = False @@ -1160,8 +1159,7 @@ def reset(self) -> None: # Clear the forces. self._nonbonded_force = None - self._custom_coulomb_force = None - self._custom_lj_force = None + self._custom_nonbonded_force = None # Clear the OpenMM context. self._openmm_context = None @@ -2184,38 +2182,32 @@ def _initialise_gpu_memory(self): map=_map, ) - # Flags for the required force. - has_gng_coul = False - has_gng_lj = False + # Flag for the required force. + has_gng = False # Find the required forces. for force in d.context().getSystem().getForces(): - if force.getName() == "GhostNonGhostCoulombForce": - gng_coul_force = force - has_gng_coul = True - elif force.getName() == "GhostNonGhostLJForce": - gng_lj_force = force - has_gng_lj = True + if force.getName() == "GhostNonGhostNonbondedForce": + gng_force = force + has_gng = True + break # Make sure the force was found. - if not has_gng_coul: - raise ValueError( - "Could not find the GhostNonGhostCoulombForce in the system" - ) - if not has_gng_lj: + if not has_gng: raise ValueError( - "Could not find the GhostNonGhostLJForce in the system" + "Could not find the GhostNonGhostNonbondedForce in the system" ) - # Get the parameters for the GhostNonGhost nonbonded forces. + # Get the parameters for the GhostNonGhostNonbondedForce. charges = _np.zeros(self._num_atoms, dtype=_np.float32) sigmas = _np.zeros(self._num_atoms, dtype=_np.float32) epsilons = _np.zeros(self._num_atoms, dtype=_np.float32) alphas = _np.zeros(self._num_atoms, dtype=_np.float32) - for i in range(gng_coul_force.getNumParticles()): + for i in range(gng_force.getNumParticles()): # Custom force parameters are returned as floats. - q, alpha, _ = gng_coul_force.getParticleParameters(i) - half_sigma, two_sqrt_epsilon, _ = gng_lj_force.getParticleParameters(i) + q, half_sigma, two_sqrt_epsilon, alpha, _ = ( + gng_force.getParticleParameters(i) + ) # Charge in |e|, sigma in nm, epsilon in kJ/mol. charges[i] = q # Rescale and convert units. @@ -2449,20 +2441,14 @@ def _accept_insertion( ) # Update the custom NonBondedForce parameters. if self._is_fep: - self._custom_coulomb_force.setParticleParameters( + self._custom_nonbonded_force.setParticleParameters( start_idx + i, ( self._water_charge[i], - 0.0, - 0.0, - ), - ) - self._custom_lj_force.setParticleParameters( - start_idx + i, - ( self._water_sigma_custom[i], self._water_epsilon_custom[i], 0.0, + 0.0, ), ) @@ -2474,8 +2460,7 @@ def _accept_insertion( # Update the CustomNonbondedForce parameters in the context. if self._is_fep: - self._custom_coulomb_force.updateParametersInContext(context) - self._custom_lj_force.updateParametersInContext(context) + self._custom_nonbonded_force.updateParametersInContext(context) # Update the state of the water on the GPU. self._kernels["update_water"]( @@ -2531,20 +2516,14 @@ def _accept_deletion(self, idx, context): ) # Update the CustomNonBondedForce parameters. if self._is_fep: - self._custom_coulomb_force.setParticleParameters( + self._custom_nonbonded_force.setParticleParameters( start_idx + i, ( 0.0, - 0.0, - 0.0, - ), - ) - self._custom_lj_force.setParticleParameters( - start_idx + i, - ( self._water_sigma_custom[i], 0.0, 0.0, + 0.0, ), ) @@ -2553,8 +2532,7 @@ def _accept_deletion(self, idx, context): # Update the CustomNonbondedForce parameters in the context. if self._is_fep: - self._custom_coulomb_force.updateParametersInContext(context) - self._custom_lj_force.updateParametersInContext(context) + self._custom_nonbonded_force.updateParametersInContext(context) # Update the state of the water on the GPU. self._kernels["update_water"]( @@ -2615,20 +2593,14 @@ def _reject_deletion(self, idx, context): ) # Update the CustomNonBondedForce parameters. if self._is_fep: - self._custom_coulomb_force.setParticleParameters( + self._custom_nonbonded_force.setParticleParameters( start_idx + i, ( self._water_charge[i], - 0.0, - 0.0, - ), - ) - self._custom_lj_force.setParticleParameters( - start_idx + i, - ( self._water_sigma_custom[i], self._water_epsilon_custom[i], 0.0, + 0.0, ), ) @@ -2637,8 +2609,7 @@ def _reject_deletion(self, idx, context): # Update the CustomNonbondedForce parameters in the context. if self._is_fep: - self._custom_coulomb_force.updateParametersInContext(context) - self._custom_lj_force.updateParametersInContext(context) + self._custom_nonbonded_force.updateParametersInContext(context) # Update the state of the water on the GPU. self._kernels["update_water"]( @@ -2705,8 +2676,7 @@ def _set_water_state(self, context, indices=None, states=None, force=False): # Assume the context has been recreated, so we need to get the # new forces. self._nonbonded_force = None - self._custom_coulomb_force = None - self._custom_lj_force = None + self._custom_nonbonded_force = None # Update even if the state is unchanged. force = True @@ -2736,20 +2706,14 @@ def _set_water_state(self, context, indices=None, states=None, force=False): ) # Update the CustomNonbondedForce parameters. if self._is_fep: - self._custom_coulomb_force.setParticleParameters( + self._custom_nonbonded_force.setParticleParameters( start_idx + i, ( 0.0, - 0.0, - 0.0, - ), - ) - self._custom_lj_force.setParticleParameters( - start_idx + i, - ( self._water_sigma_custom[i], 0.0, 0.0, + 0.0, ), ) @@ -2789,20 +2753,14 @@ def _set_water_state(self, context, indices=None, states=None, force=False): ) # Update the CustomNonBondedForce parameters. if self._is_fep: - self._custom_coulomb_force.setParticleParameters( + self._custom_nonbonded_force.setParticleParameters( start_idx + i, ( self._water_charge[i], - 0.0, - 0.0, - ), - ) - self._custom_lj_force.setParticleParameters( - start_idx + i, - ( self._water_sigma_custom[i], self._water_epsilon_custom[i], 0.0, + 0.0, ), ) @@ -2838,8 +2796,7 @@ def _set_water_state(self, context, indices=None, states=None, force=False): # Update the CustomNonbondedForce parameters in the context. if self._is_fep: - self._custom_coulomb_force.updateParametersInContext(context) - self._custom_lj_force.updateParametersInContext(context) + self._custom_nonbonded_force.updateParametersInContext(context) def _set_nonbonded_forces(self, context): """ @@ -2852,15 +2809,13 @@ def _set_nonbonded_forces(self, context): The OpenMM context to use. """ if self._nonbonded_force is None or ( - self._is_fep and self._custom_coulomb_force is None + self._is_fep and self._custom_nonbonded_force is None ): for force in context.getSystem().getForces(): if isinstance(force, _openmm.NonbondedForce): self._nonbonded_force = force - elif self._is_fep and force.getName() == "GhostNonGhostCoulombForce": - self._custom_coulomb_force = force - elif self._is_fep and force.getName() == "GhostNonGhostLJForce": - self._custom_lj_force = force + elif self._is_fep and force.getName() == "GhostNonGhostNonbondedForce": + self._custom_nonbonded_force = force elif "Barostat" in force.getName(): msg = ( f"GCMC must be used at constant volume: " @@ -2874,9 +2829,7 @@ def _set_nonbonded_forces(self, context): _logger.error(msg) raise ValueError(msg) - if self._is_fep and ( - self._custom_coulomb_force is None or self._custom_lj_force is None - ): + if self._is_fep and self._custom_nonbonded_force is None: msg = "Could not find a CustomNonbondedForce in the system" _logger.error(msg) raise ValueError(msg) From 12a8c94a955cb909a59eb083cc872d7203f803db Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 6 May 2026 16:04:06 +0100 Subject: [PATCH 22/36] Use Beutler LJ soft-core for decoupling and add support for LRC. --- src/loch/_utils.py | 24 ++++++++++++++++++++++-- 1 file changed, 22 insertions(+), 2 deletions(-) diff --git a/src/loch/_utils.py b/src/loch/_utils.py index c84b1c7..2da9cbd 100644 --- a/src/loch/_utils.py +++ b/src/loch/_utils.py @@ -36,6 +36,7 @@ def excess_chemical_potential( runtime: str = "5 ns", num_lambda: int = 24, replica_exchange: bool = False, + use_dispersion_correction: bool = False, work_dir: _Optional[str] = None, ) -> float: """ @@ -66,6 +67,9 @@ def excess_chemical_potential( replica_exchange: bool, optional Whether to use replica exchange during the calculation (default is False). + use_dispersion_correction: bool, optional + Whether to include the long-range dispersion correction (default is False). + work_dir: str, optional Working directory for the decoupling simulation (default is None, which uses a temporary directory). @@ -120,6 +124,9 @@ def excess_chemical_potential( if not isinstance(replica_exchange, bool): raise TypeError("'replica_exchange' must be a of type 'bool'.") + if not isinstance(use_dispersion_correction, bool): + raise TypeError("'use_dispersion_correction' must be a of type 'bool'.") + if work_dir is not None: if not isinstance(work_dir, str): raise TypeError("'work_dir' must be a of type 'str'.") @@ -154,7 +161,8 @@ def excess_chemical_potential( runtime=runtime, timestep="2 fs", h_mass_factor=1, - shift_delta="2.25 A", + soft_core_form="beutler", + use_dispersion_correction=use_dispersion_correction, output_directory=work_dir, ) except Exception as e: @@ -229,6 +237,7 @@ def standard_volume( cutoff: str = "10 A", num_samples: int = 5000, sample_interval: str = "1 ps", + use_dispersion_correction: bool = False, ) -> float: """ Calculate the standard volume of water at the given temperature and pressure. @@ -254,6 +263,9 @@ def standard_volume( sample_interval: str, optional Interval at which to sample the volume (default is "1 ps"). + use_dispersion_correction: bool, optional + Whether to include the long-range dispersion correction (default is False). + Returns ------- @@ -304,6 +316,9 @@ def standard_volume( if not u.has_same_units(sr.units.picosecond): raise ValueError("'sample_interval' has incorrect units.") + if not isinstance(use_dispersion_correction, bool): + raise TypeError("'use_dispersion_correction' must be a of type 'bool'.") + # Disable the dynamics progress bar. sr.base.ProgressBar.set_silent() @@ -319,7 +334,12 @@ def standard_volume( # Set up the NPT simulation. try: - d = system.dynamics(temperature=temperature, pressure=pressure, timestep="2 fs") + d = system.dynamics( + temperature=temperature, + pressure=pressure, + timestep="2 fs", + map={"use_dispersion_correction": use_dispersion_correction}, + ) except Exception as e: raise ValueError(f"Unable to set up NPT dynamics: {e}") From 80636a3d96a7b1cf8054c150c213b1d40afadc66 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 7 May 2026 15:08:34 +0100 Subject: [PATCH 23/36] Reverse lambda schedule when swap_end_states=True. --- src/loch/_sampler.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/loch/_sampler.py b/src/loch/_sampler.py index d367194..297bf9a 100644 --- a/src/loch/_sampler.py +++ b/src/loch/_sampler.py @@ -545,6 +545,9 @@ def __init__( raise ValueError("'swap_end_states' must be of type 'bool'") self._swap_end_states = swap_end_states + if swap_end_states and self._lambda_schedule is not None: + self._lambda_schedule = self._lambda_schedule.reverse() + # Check for waters and validate the template. try: self._water_template = system["water and not property is_perturbable"][0] From fcd8e96ad5cd545024d169fcde4357bc72f19786 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 7 May 2026 15:14:08 +0100 Subject: [PATCH 24/36] Update CHANGELOG. --- CHANGELOG.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index d5cee1e..319433c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,8 @@ Changelog * Add support for getting and restoring sampling statistics. * Handle XED force field virtual sites. +* Add support for long-range Lennard-Jones dispersion correction. +* Add support for Beutler soft-core Lennard-Jones form. [2025.2.0](https://github.com/openbiosim/loch/compare/2025.1.0...2025.2.0) - Feb 2026 ------------------------------------------------------------------------------------- From 87b788261a34353873da24210a725abb641c5a10 Mon Sep 17 00:00:00 2001 From: Zheng Gong Date: Tue, 26 May 2026 00:04:25 +0800 Subject: [PATCH 25/36] fix water template for dry system --- src/loch/_sampler.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/loch/_sampler.py b/src/loch/_sampler.py index 297bf9a..2724506 100644 --- a/src/loch/_sampler.py +++ b/src/loch/_sampler.py @@ -558,7 +558,7 @@ def __init__( "Please provide a water template." ) else: - if not isinstance(water_template, _sr.molecule.Molecule): + if not isinstance(water_template, _sr.mol.Molecule): raise ValueError( "'water_template' must be of type 'sire.mol.Molecule'" ) From 3f0160c1e35bdfde521a94e4bbd803a6d8441c3d Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 25 May 2026 17:54:14 +0100 Subject: [PATCH 26/36] Update CHANGELOG. [ci skip] --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 319433c..4c77fb9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,7 @@ Changelog * Handle XED force field virtual sites. * Add support for long-range Lennard-Jones dispersion correction. * Add support for Beutler soft-core Lennard-Jones form. +* Fixed type check for ``water_template``. [2025.2.0](https://github.com/openbiosim/loch/compare/2025.1.0...2025.2.0) - Feb 2026 ------------------------------------------------------------------------------------- From 59e3195748313b2c60faa491bf55b20669520eef Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 26 May 2026 08:28:18 +0100 Subject: [PATCH 27/36] Fix softcore_form kwarg and use native SOMD2 decoupling schedule. --- src/loch/_utils.py | 30 ++---------------------------- 1 file changed, 2 insertions(+), 28 deletions(-) diff --git a/src/loch/_utils.py b/src/loch/_utils.py index 2da9cbd..5f6fe85 100644 --- a/src/loch/_utils.py +++ b/src/loch/_utils.py @@ -161,7 +161,8 @@ def excess_chemical_potential( runtime=runtime, timestep="2 fs", h_mass_factor=1, - soft_core_form="beutler", + lambda_schedule="decouple", + softcore_form="beutler", use_dispersion_correction=use_dispersion_correction, output_directory=work_dir, ) @@ -181,33 +182,6 @@ def excess_chemical_potential( system.update(mol) system = sr.morph.link_to_reference(system) - # Get the lambda schedule from the molecule. - s = mol.property("schedule") - - # Avoid scaling kappa during decouple stage. - s.set_equation(stage="decouple", lever="kappa", force="ghost/ghost", equation=0) - s.set_equation(stage="decouple", lever="kappa", force="ghost-14", equation=0) - - # Add new discharging stage. - s.set_equation(stage="decouple", lever="charge", equation=s.final()) - s.prepend_stage("decharge", s.initial()) - s.set_equation( - stage="decharge", - lever="charge", - equation=s.lam() * s.final() + s.initial() * (1 - s.lam()), - ) - s.set_equation(stage="decharge", force="ghost/ghost", equation=s.initial()) - s.set_equation(stage="decharge", force="ghost-14", equation=s.initial()) - s.set_equation( - stage="decharge", lever="kappa", force="ghost/ghost", equation=-s.lam() + 1 - ) - s.set_equation( - stage="decharge", lever="kappa", force="ghost-14", equation=-s.lam() + 1 - ) - - # Update the schedule in the configuration. - config.lambda_schedule = s - # Set up the runner. try: runner = Runner(system, config) From 43349d2cf0f39e163315c4159215f6e6534f9a7e Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 26 May 2026 21:49:42 +0100 Subject: [PATCH 28/36] Add support for the osmotic ensemble. --- CHANGELOG.md | 1 + src/loch/_sampler.py | 139 ++++++++++++++++++++++++++++++------------ tests/test_osmotic.py | 70 +++++++++++++++++++++ 3 files changed, 172 insertions(+), 38 deletions(-) create mode 100644 tests/test_osmotic.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 4c77fb9..b3e4c90 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,7 @@ Changelog * Add support for long-range Lennard-Jones dispersion correction. * Add support for Beutler soft-core Lennard-Jones form. * Fixed type check for ``water_template``. +* Add support for simulations in the osmotic ensemble. [2025.2.0](https://github.com/openbiosim/loch/compare/2025.1.0...2025.2.0) - Feb 2026 ------------------------------------------------------------------------------------- diff --git a/src/loch/_sampler.py b/src/loch/_sampler.py index 2724506..79166f0 100644 --- a/src/loch/_sampler.py +++ b/src/loch/_sampler.py @@ -39,6 +39,16 @@ from ._platforms._rng import RNGManager as _RNGManager from ._softcore import SoftcoreForm as _SoftcoreForm +import sys as _sys + +try: + "μ".encode(_sys.stdout.encoding or "utf-8") + _mu_sym = "μ" +except (UnicodeEncodeError, LookupError): + _mu_sym = "mu" + +del _sys + def _as_float32(arr: _np.ndarray) -> _np.ndarray: """Convert array to float32 only if not already float32.""" @@ -65,6 +75,7 @@ def __init__( excess_chemical_potential: str = "-6.09 kcal/mol", standard_volume: str = "30.543 A^3", temperature: str = "298 K", + pressure: _Optional[str] = None, adams_shift: _Union[int, float] = 0.0, num_ghost_waters: int = 20, batch_size: int = 1000, @@ -126,6 +137,12 @@ def __init__( temperature: str The temperature of the system. + pressure: str + The pressure for NPT simulation, e.g. "1 atm". If None (default), + the simulation runs in the μVT ensemble with a fixed box. When set, + the box is updated from the context at each move and the Adams value + is recomputed accordingly, sampling the μpT (osmotic) ensemble. + adams_shift: float The Adams shift. @@ -330,6 +347,16 @@ def __init__( except Exception as e: raise ValueError(f"Could not validate the 'temperature': {e}") + if pressure is not None: + try: + self._pressure = self._validate_sire_unit( + "pressure", pressure, _sr.u("atm") + ) + except Exception as e: + raise ValueError(f"Could not validate the 'pressure': {e}") + else: + self._pressure = None + if not isinstance(num_ghost_waters, int): raise ValueError("'num_ghost_waters' must be of type 'int'") self._num_ghost_waters = num_ghost_waters @@ -685,27 +712,8 @@ def __init__( * _openmm.unit.kelvin ) - # Work out the volume of the system and GCMC sphere. - volume = self._space.volume().value() - gcmc_volume = (4.0 * _np.pi * self._radius.value() ** 3) / 3.0 - - # Work out the Adams value. - B = ( - self._beta * self._excess_chemical_potential.value() - + _np.log(gcmc_volume / self._standard_volume.value()) - ) + self._adams_shift - - # Work out the bulk Adams value. - B_bulk = ( - self._beta * self._excess_chemical_potential.value() - + _np.log(volume / self._standard_volume.value()) - ) + self._adams_shift - - # Store the exponentials for the Adams values. - self._exp_B = _np.exp(B) - self._exp_minus_B = _np.exp(-B) - self._exp_B_bulk = _np.exp(B_bulk) - self._exp_minus_B_bulk = _np.exp(-B_bulk) + # Compute the Adams values from the current box and sphere volume. + self._compute_adams_values(self._system) # Coulomb energy prefactor. self._prefactor = 1.0 / (4.0 * _np.pi * _sr.units.epsilon0.value()) @@ -740,9 +748,6 @@ def __init__( self._log_file, level=self._log_level.upper(), filter="loch" ) - # Log the Adams value. - _logger.debug(f"Adams value: {B:.6f}") - import atexit # Register the cleanup function. @@ -772,6 +777,7 @@ def __str__(self) -> str: f"excess_chemical_potential={self._excess_chemical_potential}, " f"standard_volume={self._standard_volume}, " f"temperature={self._temperature}, " + f"pressure={self._pressure}, " f"num_ghost_waters={self._num_ghost_waters}, " f"adams_shift={self._adams_shift}, " f"batch_size={self._batch_size}, " @@ -872,6 +878,55 @@ def compiler_log(self) -> str: """ return self._backend.compiler_log + def _compute_adams_values(self, system: _Any) -> None: + """ + Compute the Adams values from the current box and sphere volume. + + Updates self._exp_B, self._exp_minus_B, self._exp_B_bulk, and + self._exp_minus_B_bulk. + + Parameters + ---------- + + system: sire.system.System, openmm.Context, openmm.State + The molecular system, OpenMM context, or OpenMM state. + """ + if isinstance(system, _sr.system.System): + volume = system.property("space").volume().value() + elif isinstance(system, _openmm.Context): + box = system.getState().getPeriodicBoxVectors(asNumpy=True) + volume = _np.linalg.det(box / _openmm.unit.angstrom) + elif isinstance(system, _openmm.State): + box = system.getPeriodicBoxVectors(asNumpy=True) + volume = _np.linalg.det(box / _openmm.unit.angstrom) + else: + raise ValueError( + "'system' must be of type 'sire.system.System', 'openmm.Context', " + "or 'openmm.State'" + ) + + gcmc_volume = (4.0 * _np.pi * self._radius.value() ** 3) / 3.0 + + B = ( + self._beta * self._excess_chemical_potential.value() + + _np.log(gcmc_volume / self._standard_volume.value()) + ) + self._adams_shift + + B_bulk = ( + self._beta * self._excess_chemical_potential.value() + + _np.log(volume / self._standard_volume.value()) + ) + self._adams_shift + + self._exp_B = _np.exp(B) + self._exp_minus_B = _np.exp(-B) + self._exp_B_bulk = _np.exp(B_bulk) + self._exp_minus_B_bulk = _np.exp(-B_bulk) + + _logger.debug(f"Adams value: {B:.6f}") + + # Store the box volume in nm^3 for reuse in LRC corrections. + self._v_nm3 = volume / 1000.0 + def set_box(self, system: _Any) -> None: """ Set the box information. @@ -879,8 +934,8 @@ def set_box(self, system: _Any) -> None: Parameters ---------- - system: sire.system.System, openmm.Context - The molecular system, or OpenMM context. + system: sire.system.System, openmm.Context, openmm.State + The molecular system, OpenMM context, or OpenMM state. """ # Get the space property from the system. @@ -890,8 +945,12 @@ def set_box(self, system: _Any) -> None: except Exception: raise ValueError("'system' must contain a 'space' property") # Create a Sire TriclinicBox from the OpenMM box vectors. - elif isinstance(system, _openmm.Context): - box = system.getState().getPeriodicBoxVectors() + elif isinstance(system, (_openmm.Context, _openmm.State)): + box = ( + system.getState().getPeriodicBoxVectors() + if isinstance(system, _openmm.Context) + else system.getPeriodicBoxVectors() + ) v0 = [10 * box[0].x, 10 * box[0].y, 10 * box[0].z] v1 = [10 * box[1].x, 10 * box[1].y, 10 * box[1].z] v2 = [10 * box[2].x, 10 * box[2].y, 10 * box[2].z] @@ -900,7 +959,8 @@ def set_box(self, system: _Any) -> None: ) else: raise ValueError( - "'system' must be of type 'sire.system.System' or 'openmm.Context'" + "'system' must be of type 'sire.system.System', 'openmm.Context', " + "or 'openmm.State'" ) # Get the box information. @@ -1304,10 +1364,12 @@ def move(self, context: _openmm.Context) -> list[int]: else: initial_energy = None - # Cache the box volume (NVT, so constant throughout). - if self._has_gcmc_lrc: - box = state.getPeriodicBoxVectors(asNumpy=True) - v_nm3 = _np.linalg.det(box / _openmm.unit.nanometer) + # For NPT (μpT), update the box and recompute Adams values + # from the current context state so that volume fluctuations + # driven by the barostat are reflected in the acceptance criterion. + if self._pressure is not None: + self.set_box(state) + self._compute_adams_values(state) # Sample within the GCMC sphere. if self._reference is not None and not self._is_bulk: @@ -1583,7 +1645,7 @@ def move(self, context: _openmm.Context) -> list[int]: self._lrc_w_solute + 2.0 * n_w_before_insert * self._lrc_ww_half ) - / v_nm3 + / self._v_nm3 * _openmm.unit.kilojoules_per_mole ) dE_RF += dLRC @@ -1686,7 +1748,7 @@ def move(self, context: _openmm.Context) -> list[int]: * (n_w_before_delete - 1.0) * self._lrc_ww_half ) - / v_nm3 + / self._v_nm3 * _openmm.unit.kilojoules_per_mole ) dE_RF += dLRC @@ -2819,10 +2881,11 @@ def _set_nonbonded_forces(self, context): self._nonbonded_force = force elif self._is_fep and force.getName() == "GhostNonGhostNonbondedForce": self._custom_nonbonded_force = force - elif "Barostat" in force.getName(): + elif self._pressure is None and "Barostat" in force.getName(): msg = ( - f"GCMC must be used at constant volume: " - f"'{force.getName()}' is not supported." + f"A barostat was detected but no pressure was set: " + f"'{force.getName()}' is not supported in the {_mu_sym}VT ensemble. " + f"Pass pressure= to enable {_mu_sym}pT (osmotic) sampling." ) _logger.error(msg) raise TypeError(msg) diff --git a/tests/test_osmotic.py b/tests/test_osmotic.py new file mode 100644 index 0000000..840801a --- /dev/null +++ b/tests/test_osmotic.py @@ -0,0 +1,70 @@ +import os + +import openmm +import pytest + +from loch import GCMCSampler + + +@pytest.mark.skipif( + "CUDA_VISIBLE_DEVICES" not in os.environ, + reason="Requires CUDA enabled GPU.", +) +@pytest.mark.parametrize("platform", ["cuda", "opencl"]) +def test_osmotic_ensemble(water_box, platform): + """ + When pressure is set, move() updates box parameters and Adams values + from the current OpenMM context state after each call. + + Two moves are performed with the box scaled between them. After the + second move, _v_nm3, _exp_B_bulk, and _exp_minus_B_bulk must all + reflect the new volume. + """ + mols, _ = water_box + + sampler = GCMCSampler( + mols, + cutoff_type="rf", + cutoff="10 A", + pressure="1 atm", + ghost_file=None, + log_file=None, + test=True, + platform=platform, + seed=42, + ) + + # NPT dynamics so the context carries a barostat (bypasses the muVT guard). + d = sampler.system().dynamics( + cutoff_type="rf", + cutoff="10 A", + temperature="298 K", + pressure="1 atm", + constraint="h_bonds", + timestep="2 fs", + platform=platform, + ) + context = d.context() + + # First move: record initial box-dependent values. + sampler.move(context) + v_nm3_before = sampler._v_nm3 + exp_B_bulk_before = sampler._exp_B_bulk + exp_minus_B_bulk_before = sampler._exp_minus_B_bulk + + # Scale the box uniformly to simulate a barostat volume move. + scale = 1.1 + box = context.getState().getPeriodicBoxVectors(asNumpy=True) + box_nm = box.value_in_unit(openmm.unit.nanometer) + context.setPeriodicBoxVectors( + openmm.Vec3(*box_nm[0] * scale) * openmm.unit.nanometer, + openmm.Vec3(*box_nm[1] * scale) * openmm.unit.nanometer, + openmm.Vec3(*box_nm[2] * scale) * openmm.unit.nanometer, + ) + + # Second move: box parameters must now reflect the scaled volume. + sampler.move(context) + + assert sampler._v_nm3 == pytest.approx(v_nm3_before * scale**3, rel=1e-5) + assert sampler._exp_B_bulk != exp_B_bulk_before + assert sampler._exp_minus_B_bulk != exp_minus_B_bulk_before From 2077da60209ff10ee106561e0379db0059233577 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 28 May 2026 08:23:21 +0100 Subject: [PATCH 29/36] Fix non-uniform bulk GCMC insertion positions caused by wrong RNG distribution --- CHANGELOG.md | 1 + src/loch/_kernels.py | 15 ++++++++------- src/loch/_platforms/_rng.py | 15 +++++++++++---- src/loch/_sampler.py | 8 ++++++-- tests/test_energy.py | 10 +++++----- tests/test_rng.py | 18 ++++++++++++------ 6 files changed, 43 insertions(+), 24 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index b3e4c90..2b17936 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,6 +10,7 @@ Changelog * Add support for Beutler soft-core Lennard-Jones form. * Fixed type check for ``water_template``. * Add support for simulations in the osmotic ensemble. +* Fixed non-uniform bulk insertion positions caused by use of normal rather than uniform random numbers. [2025.2.0](https://github.com/openbiosim/loch/compare/2025.1.0...2025.2.0) - Feb 2026 ------------------------------------------------------------------------------------- diff --git a/src/loch/_kernels.py b/src/loch/_kernels.py index 146accc..8959791 100644 --- a/src/loch/_kernels.py +++ b/src/loch/_kernels.py @@ -272,7 +272,8 @@ GLOBAL float* water_position, int is_target, GLOBAL float* randoms_rotation, - GLOBAL float* randoms_position, + GLOBAL float* randoms_position_sphere, + GLOBAL float* randoms_position_bulk, GLOBAL float* randoms_radius, GLOBAL const float* cell_matrix) { @@ -319,9 +320,9 @@ if (is_target == 1) { // Generate a random position around the target using pre-generated normals. - xyz[0] = randoms_position[tidx * 3]; - xyz[1] = randoms_position[tidx * 3 + 1]; - xyz[2] = randoms_position[tidx * 3 + 2]; + xyz[0] = randoms_position_sphere[tidx * 3]; + xyz[1] = randoms_position_sphere[tidx * 3 + 1]; + xyz[2] = randoms_position_sphere[tidx * 3 + 2]; float norm = sqrtf(xyz[0] * xyz[0] + xyz[1] * xyz[1] + xyz[2] * xyz[2]); xyz[0] /= norm; @@ -337,9 +338,9 @@ { // Use pre-generated uniform randoms for bulk sampling. float r[3]; - r[0] = randoms_position[tidx * 3]; - r[1] = randoms_position[tidx * 3 + 1]; - r[2] = randoms_position[tidx * 3 + 2]; + r[0] = randoms_position_bulk[tidx * 3]; + r[1] = randoms_position_bulk[tidx * 3 + 1]; + r[2] = randoms_position_bulk[tidx * 3 + 2]; for (int i = 0; i < 3; i++) { diff --git a/src/loch/_platforms/_rng.py b/src/loch/_platforms/_rng.py index 5b80e66..26ae59c 100644 --- a/src/loch/_platforms/_rng.py +++ b/src/loch/_platforms/_rng.py @@ -42,8 +42,11 @@ class BatchRandoms: rotation : numpy.ndarray Shape (batch_size, 3) array of uniform [0,1) randoms for water rotation. - position : numpy.ndarray - Shape (batch_size, 3) array of normal randoms for position direction. + position_sphere : numpy.ndarray + Shape (batch_size, 3) array of normal randoms for sphere position direction. + + position_bulk : numpy.ndarray + Shape (batch_size, 3) array of uniform [0,1) randoms for bulk box sampling. radius : numpy.ndarray Shape (batch_size,) array of uniform [0,1) randoms for radial distance. @@ -53,7 +56,8 @@ class BatchRandoms: """ rotation: _np.ndarray - position: _np.ndarray + position_sphere: _np.ndarray + position_bulk: _np.ndarray radius: _np.ndarray acceptance: _np.ndarray @@ -101,7 +105,10 @@ def _generate_batch(self) -> BatchRandoms: rotation=self._rng.uniform(0, 1, size=(self._batch_size, 3)).astype( _np.float32 ), - position=self._rng.normal(0, 1, size=(self._batch_size, 3)).astype( + position_sphere=self._rng.normal(0, 1, size=(self._batch_size, 3)).astype( + _np.float32 + ), + position_bulk=self._rng.uniform(0, 1, size=(self._batch_size, 3)).astype( _np.float32 ), radius=self._rng.uniform(0, 1, size=self._batch_size).astype(_np.float32), diff --git a/src/loch/_sampler.py b/src/loch/_sampler.py index 79166f0..0fec54b 100644 --- a/src/loch/_sampler.py +++ b/src/loch/_sampler.py @@ -1479,7 +1479,10 @@ def move(self, context: _openmm.Context) -> list[int]: # Get pre-computed random numbers for this batch. batch_randoms = self._rng_manager.get_batch_randoms() randoms_rotation = self._backend.to_gpu(batch_randoms.rotation) - randoms_position = self._backend.to_gpu(batch_randoms.position) + randoms_position_sphere = self._backend.to_gpu( + batch_randoms.position_sphere + ) + randoms_position_bulk = self._backend.to_gpu(batch_randoms.position_bulk) randoms_radius = self._backend.to_gpu(batch_randoms.radius) # Generate the random water positions and orientations. @@ -1492,7 +1495,8 @@ def move(self, context: _openmm.Context) -> list[int]: self._water_positions, is_target, randoms_rotation, - randoms_position, + randoms_position_sphere, + randoms_position_bulk, randoms_radius, self._gpu_cell_matrix, block=(self._num_threads, 1, 1), diff --git a/tests/test_energy.py b/tests/test_energy.py index 5493069..1dc878b 100644 --- a/tests/test_energy.py +++ b/tests/test_energy.py @@ -285,17 +285,17 @@ def test_platform_consistency(fixture, request): ) -# Reference energy values captured with seed=42 on the original kernel implementation. +# Reference energy values captured with seed=42. # These anchor the kernel output to exact values so that refactors (e.g. moving from # __device__ static arrays to buffer arguments) can be validated. _REFERENCE_ENERGIES = { "water_box": { - "energy_coul": -9.45853172201302, - "energy_lj": 3.2191088, + "energy_coul": -4.674884, + "energy_lj": 0.82380486, }, "bpti": { - "energy_coul": -15.377882774621897, - "energy_lj": -0.58867246, + "energy_coul": -13.205343, + "energy_lj": 5.061536, }, } diff --git a/tests/test_rng.py b/tests/test_rng.py index cbb514f..6dda0a8 100644 --- a/tests/test_rng.py +++ b/tests/test_rng.py @@ -11,12 +11,14 @@ def test_dataclass_fields(self): """Test that BatchRandoms has the expected fields.""" batch = BatchRandoms( rotation=np.zeros((10, 3)), - position=np.zeros((10, 3)), + position_sphere=np.zeros((10, 3)), + position_bulk=np.zeros((10, 3)), radius=np.zeros(10), acceptance=np.zeros(10), ) assert hasattr(batch, "rotation") - assert hasattr(batch, "position") + assert hasattr(batch, "position_sphere") + assert hasattr(batch, "position_bulk") assert hasattr(batch, "radius") assert hasattr(batch, "acceptance") @@ -32,7 +34,8 @@ def test_batch_shapes(self): batch = rng.get_batch_randoms() assert batch.rotation.shape == (batch_size, 3) - assert batch.position.shape == (batch_size, 3) + assert batch.position_sphere.shape == (batch_size, 3) + assert batch.position_bulk.shape == (batch_size, 3) assert batch.radius.shape == (batch_size,) assert batch.acceptance.shape == (batch_size,) @@ -45,7 +48,8 @@ def test_batch_dtypes(self): batch = rng.get_batch_randoms() assert batch.rotation.dtype == np.float32 - assert batch.position.dtype == np.float32 + assert batch.position_sphere.dtype == np.float32 + assert batch.position_bulk.dtype == np.float32 assert batch.radius.dtype == np.float32 assert batch.acceptance.dtype == np.float32 @@ -60,6 +64,7 @@ def test_uniform_range(self): batch = rng.get_batch_randoms() assert np.all(batch.rotation >= 0) and np.all(batch.rotation < 1) + assert np.all(batch.position_bulk >= 0) and np.all(batch.position_bulk < 1) assert np.all(batch.radius >= 0) and np.all(batch.radius < 1) assert np.all(batch.acceptance >= 0) and np.all(batch.acceptance < 1) @@ -73,7 +78,7 @@ def test_normal_distribution(self): samples = [] for _ in range(10): batch = rng.get_batch_randoms() - samples.append(batch.position.flatten()) + samples.append(batch.position_sphere.flatten()) all_samples = np.concatenate(samples) @@ -92,7 +97,8 @@ def test_reproducibility_with_seed(self): batch2 = rng2.get_batch_randoms() np.testing.assert_array_equal(batch1.rotation, batch2.rotation) - np.testing.assert_array_equal(batch1.position, batch2.position) + np.testing.assert_array_equal(batch1.position_sphere, batch2.position_sphere) + np.testing.assert_array_equal(batch1.position_bulk, batch2.position_bulk) np.testing.assert_array_equal(batch1.radius, batch2.radius) np.testing.assert_array_equal(batch1.acceptance, batch2.acceptance) From be2218ac2d907fe334b148e3d862d68df5cbbd29 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 15 Jun 2026 14:44:56 +0100 Subject: [PATCH 30/36] Add method to update system with current water state and positions. --- src/loch/_sampler.py | 83 +++++++++++++++++++++++++++++++++++++ tests/test_energy.py | 97 ++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 180 insertions(+) diff --git a/src/loch/_sampler.py b/src/loch/_sampler.py index 0fec54b..d0cd3b0 100644 --- a/src/loch/_sampler.py +++ b/src/loch/_sampler.py @@ -863,6 +863,86 @@ def system(self) -> _Any: """ return self._system.clone() + def update_system(self, context: _Optional[_openmm.Context] = None) -> _Any: + """ + Update the system with the current water state and (optionally) positions. + + Parameters + ---------- + + context: openmm.Context + The OpenMM context containing the current positions. + + Returns + ------- + + system: sire.system.System + The updated GCMC system. + """ + + # Clone the system. + system = self._system.clone() + + if context is not None: + if not isinstance(context, _openmm.Context): + raise ValueError("'context' must be of type 'openmm.Context'") + + from sire.legacy.IO import setCoordinates + + # Get the current OpenMM positions in Angstrom. + positions = ( + context.getState(getPositions=True).getPositions(asNumpy=True) + / _openmm.unit.angstrom + ).tolist() + + try: + # Update the system with the current positions. + system._system = setCoordinates( + self._system._system, + positions, + ) + except Exception as e: + raise ValueError( + f"Could not update the system with the current positions: {e}" + ) + + # Flag the ghost waters in the system according to the current water state. + system = self._flag_ghost_waters(system) + + # Turn OFF interactions for ghost waters in the system. This is done + # by setting the charges and Lennard-Jones parameters of ghost waters + # to zero. + for mol in system["property is_ghost_water"].molecules(): + cursor = mol.cursor() + for atom in cursor.atoms(): + LJ = atom["LJ"] + atom["charge"] = 0.0 * _sr.units.mod_electron + atom["LJ"] = _sr.legacy.MM.LJParameter( + LJ.sigma(), 0.0 * _sr.units.kcal_per_mol + ) + mol = cursor.commit() + system.update(mol) + + # Turn ON any ghost buffer waters that were inserted during sampling. + # These are waters in the last _num_ghost_waters slots that now have state=1. + first_ghost_buffer = self._num_waters - self._num_ghost_waters + for i in range(first_ghost_buffer, self._num_waters): + if self._water_state[i] == 1: + mol = system[ + system.atoms()[int(self._water_indices_sire[i])].molecule() + ] + cursor = mol.cursor() + for j, atom in enumerate(cursor.atoms()): + atom["charge"] = self._water_charge[j] * _sr.units.mod_electron + atom["LJ"] = _sr.legacy.MM.LJParameter( + self._water_sigma[j] * _sr.units.angstrom, + self._water_epsilon[j] * _sr.units.kcal_per_mol, + ) + mol = cursor.commit() + system.update(mol) + + return system + def compiler_log(self) -> str: """ Return the GPU kernel compiler log. @@ -3091,6 +3171,9 @@ def _flag_ghost_waters(self, system): if not isinstance(system, _sr.system.System): raise ValueError("'system' must be a Sire system") + # Clone the system so that we don't modify the original. + system = system.clone() + # Use the Sire atom indices (no vsite offset) so that lookups into the # input topology are correct regardless of virtual sites in the context. ghost_oxygens = self._water_indices_sire[self._get_ghost_waters()] diff --git a/tests/test_energy.py b/tests/test_energy.py index 1dc878b..c0e2f93 100644 --- a/tests/test_energy.py +++ b/tests/test_energy.py @@ -1,8 +1,10 @@ import math import openmm +import openmm.unit as omm_unit import os import pytest +import numpy as np import sire as sr from loch import GCMCSampler, SoftcoreForm @@ -458,3 +460,98 @@ def _create_and_run(seed): assert math.isclose(energy1_lj, energy2_lj, abs_tol=1e-4), ( f"LJ energy mismatch: {energy1_lj!r} vs {energy2_lj!r}" ) + + +@pytest.mark.skipif( + "CUDA_VISIBLE_DEVICES" not in os.environ, + reason="Requires CUDA enabled GPU.", +) +@pytest.mark.parametrize("platform", ["cuda", "opencl"]) +def test_update_system_insertion(platform, water_box): + """ + After an accepted insertion of a ghost buffer water, update_system() must: + - restore real LJ/charge parameters on the inserted water, and + - reflect the current OpenMM coordinates for that water. + """ + mols, reference = water_box + + sampler = GCMCSampler( + mols, + cutoff_type="rf", + cutoff="10 A", + reference=reference, + log_level="debug", + ghost_file=None, + log_file=None, + test=True, + platform=platform, + ) + + d = sampler.system().dynamics( + cutoff_type="rf", + cutoff="10 A", + temperature="298 K", + pressure=None, + constraint="h_bonds", + timestep="2 fs", + platform=platform, + ) + + ctx = d.context() + first_ghost_buffer = sampler._num_waters - sampler._num_ghost_waters + + # Loop until a ghost buffer water is inserted. + inserted_idx = None + while inserted_idx is None: + state_before = sampler.water_state() + moves = sampler.move(ctx) + if len(moves) == 0 or moves[0] != 0: + continue + state_after = sampler.water_state() + # Find a buffer-slot water that changed from 0 → 1. + for i in range(first_ghost_buffer, sampler._num_waters): + if state_before[i] == 0 and state_after[i] == 1: + inserted_idx = i + break + + assert inserted_idx is not None, "No buffer water was inserted" + + # Call update_system and get the returned system. + system = sampler.update_system(ctx) + + # Get the sire atom index (no vsite offset) of the inserted water's oxygen. + o_idx = int(sampler._water_indices_sire[inserted_idx]) + inserted_mol = system[system.atoms()[o_idx].molecule()] + + # Check parameters are real (non-zero). + for j, atom in enumerate(inserted_mol.atoms()): + charge = atom.property("charge").value() + lj = atom.property("LJ") + epsilon = lj.epsilon().value() + + assert charge == pytest.approx(sampler._water_charge[j], abs=1e-6), ( + f"Atom {j} charge not restored: got {charge}, expected {sampler._water_charge[j]}" + ) + assert epsilon == pytest.approx(sampler._water_epsilon[j], abs=1e-6), ( + f"Atom {j} epsilon not restored: got {epsilon}, expected {sampler._water_epsilon[j]}" + ) + + # Check coordinates match the OpenMM context. + omm_positions = ( + ctx.getState(getPositions=True, enforcePeriodicBox=False) + .getPositions(asNumpy=True) + .value_in_unit(omm_unit.angstrom) + ) + + # _water_indices uses the vsite-offset index into the OpenMM atom list. + omm_o_idx = int(sampler._water_indices[inserted_idx]) + + for j, atom in enumerate(inserted_mol.atoms()): + sire_coord = atom.property("coordinates") + sire_xyz = np.array( + [sire_coord.x().value(), sire_coord.y().value(), sire_coord.z().value()] + ) + omm_xyz = omm_positions[omm_o_idx + j] + assert np.allclose(sire_xyz, omm_xyz, atol=1e-3), ( + f"Atom {j} coordinates mismatch: sire={sire_xyz}, omm={omm_xyz}" + ) From 482a84d2582b87149ad38d75229e3890581d2581 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 15 Jun 2026 16:08:29 +0100 Subject: [PATCH 31/36] Add method to return system without ghost waters. --- src/loch/_sampler.py | 56 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 56 insertions(+) diff --git a/src/loch/_sampler.py b/src/loch/_sampler.py index d0cd3b0..302577a 100644 --- a/src/loch/_sampler.py +++ b/src/loch/_sampler.py @@ -943,6 +943,62 @@ def update_system(self, context: _Optional[_openmm.Context] = None) -> _Any: return system + def finalise_system(self, context: _Optional[_openmm.Context] = None) -> _Any: + """ + Return the system with ghost waters removed and (optionally) positions + updated from an OpenMM context. + + Unlike :meth:`update_system`, which retains ghost waters with zeroed + parameters, this method deletes them entirely so the returned system is + suitable for use as input to non-GCMC simulations or analyses. + + Parameters + ---------- + + context: openmm.Context + The OpenMM context containing the current positions. + + Returns + ------- + + system: sire.system.System + The system with ghost waters removed. + """ + + # Clone the system. + system = self._system.clone() + + if context is not None: + if not isinstance(context, _openmm.Context): + raise ValueError("'context' must be of type 'openmm.Context'") + + from sire.legacy.IO import setCoordinates + + # Get the current OpenMM positions in Angstrom. + positions = ( + context.getState(getPositions=True).getPositions(asNumpy=True) + / _openmm.unit.angstrom + ).tolist() + + try: + system._system = setCoordinates( + self._system._system, + positions, + ) + except Exception as e: + raise ValueError( + f"Could not update the system with the current positions: {e}" + ) + + # Collect all ghost water molecules before removing any, since each + # removal shifts atom indices. + ghost_oxygens = self._water_indices_sire[self._get_ghost_waters()] + ghost_mols = [system[system.atoms()[int(i)].molecule()] for i in ghost_oxygens] + for mol in ghost_mols: + system.remove(mol) + + return system + def compiler_log(self) -> str: """ Return the GPU kernel compiler log. From ff7efee41795e74ab6f24a1c65cb83a8aee935ca Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 15 Jun 2026 16:31:11 +0100 Subject: [PATCH 32/36] Update space property. --- src/loch/_sampler.py | 90 +++++++++++++++++++++----------------------- 1 file changed, 42 insertions(+), 48 deletions(-) diff --git a/src/loch/_sampler.py b/src/loch/_sampler.py index 302577a..0c17461 100644 --- a/src/loch/_sampler.py +++ b/src/loch/_sampler.py @@ -863,24 +863,11 @@ def system(self) -> _Any: """ return self._system.clone() - def update_system(self, context: _Optional[_openmm.Context] = None) -> _Any: + def _system_from_context(self, context: _Optional[_openmm.Context]) -> _Any: """ - Update the system with the current water state and (optionally) positions. - - Parameters - ---------- - - context: openmm.Context - The OpenMM context containing the current positions. - - Returns - ------- - - system: sire.system.System - The updated GCMC system. + Clone the internal system and, if a context is provided, update + atomic coordinates and the periodic box from it. """ - - # Clone the system. system = self._system.clone() if context is not None: @@ -889,23 +876,52 @@ def update_system(self, context: _Optional[_openmm.Context] = None) -> _Any: from sire.legacy.IO import setCoordinates - # Get the current OpenMM positions in Angstrom. + state = context.getState(getPositions=True) positions = ( - context.getState(getPositions=True).getPositions(asNumpy=True) - / _openmm.unit.angstrom + state.getPositions(asNumpy=True) / _openmm.unit.angstrom ).tolist() try: - # Update the system with the current positions. - system._system = setCoordinates( - self._system._system, - positions, - ) + system._system = setCoordinates(self._system._system, positions) except Exception as e: raise ValueError( f"Could not update the system with the current positions: {e}" ) + # Update the periodic box from the context. + box = state.getPeriodicBoxVectors() + v0 = [10 * box[0].x, 10 * box[0].y, 10 * box[0].z] + v1 = [10 * box[1].x, 10 * box[1].y, 10 * box[1].z] + v2 = [10 * box[2].x, 10 * box[2].y, 10 * box[2].z] + space = _sr.vol.TriclinicBox( + _sr.maths.Vector(*v0), + _sr.maths.Vector(*v1), + _sr.maths.Vector(*v2), + ) + system.set_property("space", space) + + return system + + def update_system(self, context: _Optional[_openmm.Context] = None) -> _Any: + """ + Update the system with the current water state and (optionally) positions. + + Parameters + ---------- + + context: openmm.Context + The OpenMM context containing the current positions. + + Returns + ------- + + system: sire.system.System + The updated GCMC system. + """ + + # Clone the system and update coordinates and box from the context. + system = self._system_from_context(context) + # Flag the ghost waters in the system according to the current water state. system = self._flag_ghost_waters(system) @@ -965,30 +981,8 @@ def finalise_system(self, context: _Optional[_openmm.Context] = None) -> _Any: The system with ghost waters removed. """ - # Clone the system. - system = self._system.clone() - - if context is not None: - if not isinstance(context, _openmm.Context): - raise ValueError("'context' must be of type 'openmm.Context'") - - from sire.legacy.IO import setCoordinates - - # Get the current OpenMM positions in Angstrom. - positions = ( - context.getState(getPositions=True).getPositions(asNumpy=True) - / _openmm.unit.angstrom - ).tolist() - - try: - system._system = setCoordinates( - self._system._system, - positions, - ) - except Exception as e: - raise ValueError( - f"Could not update the system with the current positions: {e}" - ) + # Clone the system and update coordinates and box from the context. + system = self._system_from_context(context) # Collect all ghost water molecules before removing any, since each # removal shifts atom indices. From b2d10223f6e88bacff34004eb7afaeb0a5017f36 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 15 Jun 2026 16:32:28 +0100 Subject: [PATCH 33/36] Add unit test for finalise_system. --- tests/test_energy.py | 61 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 61 insertions(+) diff --git a/tests/test_energy.py b/tests/test_energy.py index c0e2f93..703a8a2 100644 --- a/tests/test_energy.py +++ b/tests/test_energy.py @@ -555,3 +555,64 @@ def test_update_system_insertion(platform, water_box): assert np.allclose(sire_xyz, omm_xyz, atol=1e-3), ( f"Atom {j} coordinates mismatch: sire={sire_xyz}, omm={omm_xyz}" ) + + +@pytest.mark.skipif( + "CUDA_VISIBLE_DEVICES" not in os.environ, + reason="Requires CUDA enabled GPU.", +) +@pytest.mark.parametrize("platform", ["cuda", "opencl"]) +def test_finalise_system(platform, water_box): + """ + After an accepted insertion of a ghost buffer water, finalise_system() must + return a system whose water count equals the number of real (non-ghost) + waters in the sampler. + """ + mols, reference = water_box + + sampler = GCMCSampler( + mols, + cutoff_type="rf", + cutoff="10 A", + reference=reference, + log_level="debug", + ghost_file=None, + log_file=None, + test=True, + platform=platform, + ) + + d = sampler.system().dynamics( + cutoff_type="rf", + cutoff="10 A", + temperature="298 K", + pressure=None, + constraint="h_bonds", + timestep="2 fs", + platform=platform, + ) + + ctx = d.context() + first_ghost_buffer = sampler._num_waters - sampler._num_ghost_waters + + # Loop until a ghost buffer water is inserted. + inserted = False + while not inserted: + state_before = sampler.water_state() + moves = sampler.move(ctx) + if len(moves) == 0 or moves[0] != 0: + continue + state_after = sampler.water_state() + for i in range(first_ghost_buffer, sampler._num_waters): + if state_before[i] == 0 and state_after[i] == 1: + inserted = True + break + + system = sampler.finalise_system(ctx) + + expected_waters = int(sampler._get_non_ghost_waters().size) + actual_waters = system["water"].num_molecules() + + assert actual_waters == expected_waters, ( + f"Water count mismatch: got {actual_waters}, expected {expected_waters}" + ) From b3bbae22259547c85a762a288878722a8dbd32e8 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 15 Jun 2026 16:35:39 +0100 Subject: [PATCH 34/36] Update CHANGELOG. --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2b17936..392c880 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,7 @@ Changelog * Fixed type check for ``water_template``. * Add support for simulations in the osmotic ensemble. * Fixed non-uniform bulk insertion positions caused by use of normal rather than uniform random numbers. +* Add methods to update the system with the current water state and return system without ghost waters. [2025.2.0](https://github.com/openbiosim/loch/compare/2025.1.0...2025.2.0) - Feb 2026 ------------------------------------------------------------------------------------- From 3ebc5f191d15aaf72db94af4a5f5d833cd507603 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 29 Jun 2026 13:12:36 +0100 Subject: [PATCH 35/36] Add BioSimSpace compatibility pin. --- pixi.toml | 5 ++++- recipes/loch/recipe.yaml | 5 ++++- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/pixi.toml b/pixi.toml index 446a378..16c6cff 100644 --- a/pixi.toml +++ b/pixi.toml @@ -5,7 +5,10 @@ platforms = ["linux-64", "osx-arm64", "win-64"] [dependencies] python = ">=3.10" -biosimspace = "*" +# main +biosimspace = ">=2026.1.0,<2026.2.0" +# devel +#biosimspace = "==2026.2.0.dev" loguru = "*" pyopencl = "*" diff --git a/recipes/loch/recipe.yaml b/recipes/loch/recipe.yaml index d840361..1e854de 100644 --- a/recipes/loch/recipe.yaml +++ b/recipes/loch/recipe.yaml @@ -19,7 +19,10 @@ requirements: - setuptools - versioningit run: - - biosimspace + # main + - biosimspace >=2026.1.0,<2026.2.0 + # devel + #- biosimspace ==2026.2.0.dev - loguru - pyopencl - python From 9fc660cb0d4bbfd983f39fa793d96d5f0371e907 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 29 Jun 2026 13:12:50 +0100 Subject: [PATCH 36/36] Update CHANGELOG for the 2026.1.0 release. --- CHANGELOG.md | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 392c880..50de0be 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,17 +1,17 @@ Changelog ========= -[2026.1.0](https://github.com/openbiosim/loch/compare/2025.2.0...2026.1.0) - ******** +[2026.1.0](https://github.com/openbiosim/loch/compare/2025.2.0...2026.1.0) - Jun 2026 ------------------------------------------------------------------------------------- * Add support for getting and restoring sampling statistics. -* Handle XED force field virtual sites. -* Add support for long-range Lennard-Jones dispersion correction. -* Add support for Beutler soft-core Lennard-Jones form. -* Fixed type check for ``water_template``. -* Add support for simulations in the osmotic ensemble. -* Fixed non-uniform bulk insertion positions caused by use of normal rather than uniform random numbers. -* Add methods to update the system with the current water state and return system without ghost waters. +* Handle XED force field virtual sites [#17](https://github.com/OpenBioSim/loch/pull/17). +* Add support for long-range Lennard-Jones dispersion correction [#18](https://github.com/OpenBioSim/loch/pull/18). +* Add support for Beutler soft-core Lennard-Jones form [#18](https://github.com/OpenBioSim/loch/pull/18). +* Fixed type check for ``water_template`` [#19](https://github.com/OpenBioSim/loch/pull/19). +* Add support for simulations in the osmotic ensemble [#22](https://github.com/OpenBioSim/loch/pull/22). +* Fixed non-uniform bulk insertion positions caused by use of normal rather than uniform random numbers [#24](https://github.com/OpenBioSim/loch/pull/24). +* Add methods to update the system with the current water state and return system without ghost waters [#26](https://github.com/OpenBioSim/loch/pull/26). [2025.2.0](https://github.com/openbiosim/loch/compare/2025.1.0...2025.2.0) - Feb 2026 -------------------------------------------------------------------------------------