From be10d29b032813d7e8c4c2461ea0bae5d97d657b Mon Sep 17 00:00:00 2001 From: nithishwer Date: Wed, 12 Nov 2025 19:46:51 +0000 Subject: [PATCH 01/40] feat: add gh action for test --- .github/workflows/test.yml | 34 ++++++++++++++++++++++++++++++++++ 1 file changed, 34 insertions(+) create mode 100644 .github/workflows/test.yml diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml new file mode 100644 index 0000000..df69044 --- /dev/null +++ b/.github/workflows/test.yml @@ -0,0 +1,34 @@ +name: Tests + +on: + push: + branches: [ main ] + pull_request: + branches: [ main ] + +jobs: + test: + runs-on: ubuntu-latest + + strategy: + matrix: + python-version: ["3.9", "3.10", "3.11"] + + steps: + - name: Checkout code + uses: actions/checkout@v4 + + - name: Set up Python + uses: actions/setup-python@v5 + with: + python-version: ${{ matrix.python-version }} + + - name: Install dependencies + run: | + python -m pip install --upgrade pip + if [ -f requirements.txt ]; then pip install -r requirements.txt; fi + pip install pytest + + - name: Run tests + run: | + pytest -v --maxfail=1 --disable-warnings From ca69f0be290445dd44a48889f6312978ab698e45 Mon Sep 17 00:00:00 2001 From: nithishwer Date: Wed, 12 Nov 2025 19:48:40 +0000 Subject: [PATCH 02/40] feat: add gh action to develop too --- .github/workflows/test.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index df69044..c11cc27 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -2,9 +2,9 @@ name: Tests on: push: - branches: [ main ] + branches: [ main, develop ] pull_request: - branches: [ main ] + branches: [ main, develop ] jobs: test: From 240883ced6c0461dfcb0356845bab192a596438d Mon Sep 17 00:00:00 2001 From: nithishwer Date: Wed, 12 Nov 2025 19:49:50 +0000 Subject: [PATCH 03/40] feat: add gh action to doc too --- .github/workflows/test.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index c11cc27..aa55b56 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -2,9 +2,9 @@ name: Tests on: push: - branches: [ main, develop ] + branches: [ main, develop, feat/documentation ] pull_request: - branches: [ main, develop ] + branches: [ main, develop, feat/documentation ] jobs: test: From e66c110c484ac412db6c75a8324ac4700b17f26f Mon Sep 17 00:00:00 2001 From: nithishwer Date: Wed, 12 Nov 2025 19:52:06 +0000 Subject: [PATCH 04/40] feat: add fepa install stel --- .github/workflows/test.yml | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index aa55b56..9cc8d09 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -28,7 +28,12 @@ jobs: python -m pip install --upgrade pip if [ -f requirements.txt ]; then pip install -r requirements.txt; fi pip install pytest - + - name: Install dependencies + run: | + python -m pip install --upgrade pip + if [ -f requirements.txt ]; then pip install -r requirements.txt; fi + pip install -e . + pip install pytest - name: Run tests run: | pytest -v --maxfail=1 --disable-warnings From bc44422f0aca4d2a1d23f3f13b1af4f271d74489 Mon Sep 17 00:00:00 2001 From: nithishwer Date: Wed, 12 Nov 2025 19:53:11 +0000 Subject: [PATCH 05/40] feat: match python version --- .github/workflows/test.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 9cc8d09..4d6d855 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -12,7 +12,7 @@ jobs: strategy: matrix: - python-version: ["3.9", "3.10", "3.11"] + python-version: ["3.11.14"] steps: - name: Checkout code From 4bb171f456f749f9c5194debf81b125dff6c43ae Mon Sep 17 00:00:00 2001 From: nithishwer Date: Wed, 12 Nov 2025 19:58:01 +0000 Subject: [PATCH 06/40] feat: match python version --- .github/workflows/test.yml | 8 ++------ pyproject.toml | 15 ++++++--------- 2 files changed, 8 insertions(+), 15 deletions(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 4d6d855..c91d188 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -12,7 +12,7 @@ jobs: strategy: matrix: - python-version: ["3.11.14"] + python-version: ["3.10", "3.11"] # 👈 matches >=3.10.0 steps: - name: Checkout code @@ -23,17 +23,13 @@ jobs: with: python-version: ${{ matrix.python-version }} - - name: Install dependencies - run: | - python -m pip install --upgrade pip - if [ -f requirements.txt ]; then pip install -r requirements.txt; fi - pip install pytest - name: Install dependencies run: | python -m pip install --upgrade pip if [ -f requirements.txt ]; then pip install -r requirements.txt; fi pip install -e . pip install pytest + - name: Run tests run: | pytest -v --maxfail=1 --disable-warnings diff --git a/pyproject.toml b/pyproject.toml index 365b18f..260a2ad 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,25 +1,22 @@ [build-system] -requires = ["setuptools >= 61.0","wheel"] +requires = ["setuptools >= 61.0", "wheel"] build-backend = "setuptools.build_meta" [project] name = "fepa" version = "0.1.0" description = "A package for analyzing the binding pocket for Free Energy simulations" -authors = [{ name = "Nithishwer Mouroug Anand", email = "nithishwer@gmail.com" }] +authors = [ + { name = "Nithishwer Mouroug Anand", email = "nithishwer@gmail.com" }, +] readme = "README.md" license = { file = "LICENSE" } requires-python = ">=3.10.0" -dependencies = [ - "pandas>=2.2.3", - "pensa>=0.6.0", - "seaborn>=0.13.2", -] +dependencies = ["pandas>=2.2.3", "pensa>=0.6.0", "seaborn>=0.13.2"] [tool.setuptools] -packages = ["fepa"] # You can replace with `find:` for auto-discovery +packages = ["fepa"] # You can replace with `find:` for auto-discovery [project.urls] Homepage = "https://github.com/Nithishwer/FEPA" Repository = "https://github.com/Nithishwer/FEPA" - From 016f3533837a3a5bc25cda5ee800514932389e37 Mon Sep 17 00:00:00 2001 From: nithishwer Date: Wed, 12 Nov 2025 19:59:11 +0000 Subject: [PATCH 07/40] feat: match python 3.10 --- .github/workflows/test.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index c91d188..5421433 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -12,7 +12,7 @@ jobs: strategy: matrix: - python-version: ["3.10", "3.11"] # 👈 matches >=3.10.0 + python-version: ["3.10"] # 👈 matches >=3.10.0 steps: - name: Checkout code From 5078670735b01359b209e94bef4bff9ef9c05e12 Mon Sep 17 00:00:00 2001 From: nithishwer Date: Wed, 12 Nov 2025 20:01:34 +0000 Subject: [PATCH 08/40] fix: clear caceh --- .github/workflows/test.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 5421433..6fe57b8 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -27,7 +27,7 @@ jobs: run: | python -m pip install --upgrade pip if [ -f requirements.txt ]; then pip install -r requirements.txt; fi - pip install -e . + pip install --no-cache-dir -e . pip install pytest - name: Run tests From 76c7fdd63462dad5169a78ca233d0ee29101b7f1 Mon Sep 17 00:00:00 2001 From: nithishwer Date: Wed, 12 Nov 2025 20:06:47 +0000 Subject: [PATCH 09/40] fix: remove setup.py --- setup.py | 24 ------------------------ 1 file changed, 24 deletions(-) delete mode 100644 setup.py diff --git a/setup.py b/setup.py deleted file mode 100644 index 8fb6f35..0000000 --- a/setup.py +++ /dev/null @@ -1,24 +0,0 @@ -from setuptools import setup, find_packages - -setup( - name="fepa", - version="0.1.0", - description="A package for analyzing the binding pocket for Free Energy simulations", - author="Nithishwer Mouroug Anand", - author_email="nithishwer@gmail.com", - long_description=open("README.md").read(), - long_description_content_type="text/markdown", - license=open("LICENSE").read(), - python_requires="==3.10.0", - install_requires=[ - "pandas>=2.2.3", - "pensa>=0.6.0", - "seaborn>=0.13.2", - ], - packages=find_packages(), - url="https://github.com/Nithishwer/FEPA", - project_urls={ - "Homepage": "https://github.com/Nithishwer/FEPA", - "Repository": "https://github.com/Nithishwer/FEPA", - }, -) From e55e176afe10ec9d3a649ac5d93368caacaee4a2 Mon Sep 17 00:00:00 2001 From: nithishwer Date: Wed, 12 Nov 2025 20:10:29 +0000 Subject: [PATCH 10/40] fix: use conda env --- .github/workflows/test.yml | 18 +++++++----------- 1 file changed, 7 insertions(+), 11 deletions(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 6fe57b8..d8ab1d8 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -10,26 +10,22 @@ jobs: test: runs-on: ubuntu-latest - strategy: - matrix: - python-version: ["3.10"] # 👈 matches >=3.10.0 - steps: - name: Checkout code uses: actions/checkout@v4 - - name: Set up Python - uses: actions/setup-python@v5 + - name: Set up Conda + uses: conda-incubator/setup-miniconda@v2 with: - python-version: ${{ matrix.python-version }} + environment-file: requirements.yml # 👈 your Conda env file + activate-environment: fepa_env # 👈 env name from the yml + python-version: 3.10 # 👈 matches your pyproject.toml - - name: Install dependencies + - name: Install editable package run: | - python -m pip install --upgrade pip - if [ -f requirements.txt ]; then pip install -r requirements.txt; fi pip install --no-cache-dir -e . - pip install pytest - name: Run tests run: | pytest -v --maxfail=1 --disable-warnings + From 9eb92ae3da2cd8d6afc287e7e0729cbd079417b9 Mon Sep 17 00:00:00 2001 From: nithishwer Date: Wed, 12 Nov 2025 20:19:51 +0000 Subject: [PATCH 11/40] fix: add salilab conda channel --- .github/workflows/test.yml | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index d8ab1d8..04ed04a 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -17,9 +17,16 @@ jobs: - name: Set up Conda uses: conda-incubator/setup-miniconda@v2 with: - environment-file: requirements.yml # 👈 your Conda env file - activate-environment: fepa_env # 👈 env name from the yml - python-version: 3.10 # 👈 matches your pyproject.toml + python-version: 3.10 + activate-environment: fepa_env + environment-file: requirements.yml + auto-update-conda: true + + - name: Configure Conda channels + run: | + conda config --add channels salilab + conda config --add channels conda-forge + conda config --add channels defaults - name: Install editable package run: | @@ -29,3 +36,4 @@ jobs: run: | pytest -v --maxfail=1 --disable-warnings + From 99c2e49cd21587abb77d1f1f2328d99f1bf61140 Mon Sep 17 00:00:00 2001 From: nithishwer Date: Wed, 12 Nov 2025 20:25:08 +0000 Subject: [PATCH 12/40] fix: loosen modeller version req --- requirements.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.yml b/requirements.yml index d1b813a..22b24f6 100644 --- a/requirements.yml +++ b/requirements.yml @@ -177,7 +177,7 @@ dependencies: - mistune=3.1.3=pyh29332c3_0 - mkl=2024.2.2=ha957f24_16 - mmtf-python=1.1.3=pyhd8ed1ab_0 - - modeller=10.6=py310h9bf148f_0 + - modeller>=10.6,<11 - mpfr=4.2.1=h90cbb55_3 - mrcfile=1.5.4=pyhd8ed1ab_0 - msgpack-python=1.1.0=py310h3788b33_0 From 9df1c2fee81daaf3822821d91ec7168379a46ae6 Mon Sep 17 00:00:00 2001 From: nithishwer Date: Wed, 12 Nov 2025 20:33:18 +0000 Subject: [PATCH 13/40] fix: loosen python requirements --- requirements.yml | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/requirements.yml b/requirements.yml index 22b24f6..7d70fd1 100644 --- a/requirements.yml +++ b/requirements.yml @@ -243,13 +243,12 @@ dependencies: - pygments=2.19.1=pyhd8ed1ab_0 - pysocks=1.7.1=pyha55dd90_7 - pytables=3.10.1=py310h1affd9f_4 - - python=3.10.0=h543edf9_3_cpython + - python=3.10 - python-constraint=1.4.0=pyhff2d567_1 - python-dateutil=2.9.0.post0=pyhff2d567_1 - python-fastjsonschema=2.21.1=pyhd8ed1ab_0 - python-json-logger=2.0.7=pyhd8ed1ab_0 - python-tzdata=2025.2=pyhd8ed1ab_0 - - python_abi=3.10=7_cp310 - pytng=0.3.3=py310he95d665_1 - pyyaml=6.0.2=py310h89163eb_2 - pyzmq=26.4.0=py310h71f11fc_0 From 8e59511c3704b53a0b71874ad9195399a7215d92 Mon Sep 17 00:00:00 2001 From: nithishwer Date: Wed, 12 Nov 2025 20:36:28 +0000 Subject: [PATCH 14/40] fix: no python req --- requirements.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/requirements.yml b/requirements.yml index 7d70fd1..50eca8f 100644 --- a/requirements.yml +++ b/requirements.yml @@ -243,7 +243,6 @@ dependencies: - pygments=2.19.1=pyhd8ed1ab_0 - pysocks=1.7.1=pyha55dd90_7 - pytables=3.10.1=py310h1affd9f_4 - - python=3.10 - python-constraint=1.4.0=pyhff2d567_1 - python-dateutil=2.9.0.post0=pyhff2d567_1 - python-fastjsonschema=2.21.1=pyhd8ed1ab_0 From 5adb2258788fdf8327b79617cc4ddd48b2bd070c Mon Sep 17 00:00:00 2001 From: nithishwer Date: Wed, 12 Nov 2025 20:38:49 +0000 Subject: [PATCH 15/40] fix: add requirements yml --- requirements.yml | 433 ++----------------------------------------- requirements_old.yml | 424 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 437 insertions(+), 420 deletions(-) create mode 100644 requirements_old.yml diff --git a/requirements.yml b/requirements.yml index 50eca8f..80eafc4 100644 --- a/requirements.yml +++ b/requirements.yml @@ -1,424 +1,17 @@ -name: fepa_env +name: fepa channels: - salilab - conda-forge - - defaults + - omnia dependencies: - - _libgcc_mutex=0.1=conda_forge - - _openmp_mutex=4.5=3_kmp_llvm - - ambertools=24.8=cuda_None_nompi_py310h834fefc_101 - - annotated-types=0.7.0=pyhd8ed1ab_1 - - anyio=4.9.0=pyh29332c3_0 - - argon2-cffi=23.1.0=pyhd8ed1ab_1 - - argon2-cffi-bindings=21.2.0=py310ha75aee5_5 - - arpack=3.9.1=nompi_hf03ea27_102 - - arrow=1.3.0=pyhd8ed1ab_1 - - asttokens=3.0.0=pyhd8ed1ab_1 - - async-lru=2.0.5=pyh29332c3_0 - - attrs=25.3.0=pyh71513ae_0 - - babel=2.17.0=pyhd8ed1ab_0 - - beautifulsoup4=4.13.4=pyha770c72_0 - - bleach=6.2.0=pyh29332c3_4 - - bleach-with-css=6.2.0=h82add2a_4 - - blosc=1.21.6=he440d0b_1 - - brotli=1.1.0=hb9d3cd8_2 - - brotli-bin=1.1.0=hb9d3cd8_2 - - brotli-python=1.1.0=py310hf71b8c6_2 - - bson=0.5.10=pyhd8ed1ab_0 - - bzip2=1.0.8=h4bc722e_7 - - c-ares=1.34.5=hb9d3cd8_0 - - c-blosc2=2.15.2=h3122c55_1 - - ca-certificates=2025.6.15=hbd8a1cb_0 - - cached-property=1.5.2=hd8ed1ab_1 - - cached_property=1.5.2=pyha770c72_1 - - cachetools=5.5.2=pyhd8ed1ab_0 - - cairo=1.18.4=h3394656_0 - - certifi=2025.6.15=pyhd8ed1ab_0 - - cffi=1.17.1=py310h8deb56e_0 - - cftime=1.6.4=py310hf462985_1 - - chardet=5.2.0=pyhd8ed1ab_3 - - colorama=0.4.6=pyhd8ed1ab_1 - - comm=0.2.2=pyhd8ed1ab_1 - - cudatoolkit=11.8.0=h4ba93d1_13 - - cycler=0.12.1=pyhd8ed1ab_1 - - cyrus-sasl=2.1.27=h54b06d7_7 - - cython=3.0.12=py310had8cdd9_0 - - debugpy=1.8.14=py310hf71b8c6_0 - - decorator=5.2.1=pyhd8ed1ab_0 - - defusedxml=0.7.1=pyhd8ed1ab_0 - - exceptiongroup=1.2.2=pyhd8ed1ab_1 - - executing=2.2.0=pyhd8ed1ab_0 - - fasteners=0.19=pyhd8ed1ab_1 - - fftw=3.3.10=nompi_hf1063bd_110 - - filelock=3.18.0=pyhd8ed1ab_0 - - flexcache=0.3=pyhd8ed1ab_1 - - flexparser=0.4=pyhd8ed1ab_1 - - font-ttf-dejavu-sans-mono=2.37=hab24e00_0 - - font-ttf-inconsolata=3.000=h77eed37_0 - - font-ttf-source-code-pro=2.038=h77eed37_0 - - font-ttf-ubuntu=0.83=h77eed37_3 - - fontconfig=2.15.0=h7e30c49_1 - - fonts-conda-ecosystem=1=0 - - fonts-conda-forge=1=0 - - fqdn=1.5.1=pyhd8ed1ab_1 - - freetype=2.13.3=ha770c72_1 - - freetype-py=2.3.0=pyhd8ed1ab_0 - - gawk=5.3.1=hcd3d067_0 - - gmp=6.3.0=hac33072_2 - - greenlet=3.2.1=py310hf71b8c6_0 - - griddataformats=1.0.2=pyhd8ed1ab_1 - - gsd=3.4.2=py310ha75aee5_0 - - gsl=2.7=he838d99_0 - - h11=0.16.0=pyhd8ed1ab_0 - - h2=4.2.0=pyhd8ed1ab_0 - - h5py=3.13.0=nompi_py310h60e0fe6_100 - - hdf4=4.2.15=h2a13503_7 - - hdf5=1.14.3=nompi_h2d575fe_109 - - hdf5-1107=1.10.7=1 - - hpack=4.1.0=pyhd8ed1ab_0 - - httpcore=1.0.9=pyh29332c3_0 - - httpx=0.28.1=pyhd8ed1ab_0 - - hyperframe=6.1.0=pyhd8ed1ab_0 - - icu=75.1=he02047a_0 - - idna=3.10=pyhd8ed1ab_1 - - importlib-metadata=8.6.1=pyha770c72_0 - - importlib_resources=6.5.2=pyhd8ed1ab_0 - - ipykernel=6.29.5=pyh3099207_0 - - ipython=8.36.0=pyh907856f_0 - - ipywidgets=8.1.7=pyhd8ed1ab_0 - - isoduration=20.11.0=pyhd8ed1ab_1 - - jedi=0.19.2=pyhd8ed1ab_1 - - jinja2=3.1.6=pyhd8ed1ab_0 - - joblib=1.4.2=pyhd8ed1ab_1 - - json5=0.12.0=pyhd8ed1ab_0 - - jsonpointer=3.0.0=py310hff52083_1 - - jsonschema=4.23.0=pyhd8ed1ab_1 - - jsonschema-specifications=2025.4.1=pyh29332c3_0 - - jsonschema-with-format-nongpl=4.23.0=hd8ed1ab_1 - - jupyter-lsp=2.2.5=pyhd8ed1ab_1 - - jupyter_client=8.6.3=pyhd8ed1ab_1 - - jupyter_core=5.7.2=pyh31011fe_1 - - jupyter_events=0.12.0=pyh29332c3_0 - - jupyter_server=2.15.0=pyhd8ed1ab_0 - - jupyter_server_terminals=0.5.3=pyhd8ed1ab_1 - - jupyterlab=4.4.1=pyhd8ed1ab_0 - - jupyterlab_pygments=0.3.0=pyhd8ed1ab_2 - - jupyterlab_server=2.27.3=pyhd8ed1ab_1 - - jupyterlab_widgets=3.0.15=pyhd8ed1ab_0 - - kernel-headers_linux-64=3.10.0=he073ed8_18 - - keyutils=1.6.1=h166bdaf_0 - - krb5=1.21.3=h659f571_0 - - lcms2=2.17=h717163a_0 - - ld_impl_linux-64=2.43=h712a8e2_4 - - lerc=4.0.0=h0aef613_1 - - libaec=1.1.3=h59595ed_0 - - libasprintf=0.23.1=h8e693c7_0 - - libblas=3.9.0=31_hfdb39a5_mkl - - libboost=1.86.0=h6c02f8c_3 - - libboost-python=1.86.0=py310ha2bacc8_3 - - libbrotlicommon=1.1.0=hb9d3cd8_2 - - libbrotlidec=1.1.0=hb9d3cd8_2 - - libbrotlienc=1.1.0=hb9d3cd8_2 - - libcblas=3.9.0=31_h372d94f_mkl - - libcurl=8.13.0=h332b0f4_0 - - libdeflate=1.23=h86f0d12_0 - - libedit=3.1.20250104=pl5321h7949ede_0 - - libev=4.33=hd590300_2 - - libexpat=2.7.0=h5888daf_0 - - libffi=3.4.6=h2dba641_0 - - libfreetype=2.13.3=ha770c72_1 - - libfreetype6=2.13.3=h48d6fc4_1 - - libgcc=14.2.0=h767d61c_2 - - libgcc-ng=14.2.0=h69a702a_2 - - libgettextpo=0.23.1=h5888daf_0 - - libgfortran=14.2.0=h69a702a_2 - - libgfortran-ng=14.2.0=h69a702a_2 - - libgfortran5=14.2.0=hf1ad2bd_2 - - libglib=2.84.1=h3618099_1 - - libgomp=14.2.0=h767d61c_2 - - libhwloc=2.11.2=default_h0d58e46_1001 - - libiconv=1.18=h4ce23a2_1 - - libjpeg-turbo=3.1.0=hb9d3cd8_0 - - liblapack=3.9.0=31_hc41d3b0_mkl - - liblzma=5.6.4=hb9d3cd8_0 - - liblzma-devel=5.6.4=hb9d3cd8_0 - - libnetcdf=4.9.2=nompi_h00e09a9_116 - - libnghttp2=1.64.0=h161d5f1_0 - - libnsl=2.0.1=hd590300_0 - - libntlm=1.8=hb9d3cd8_0 - - libopenblas=0.3.29=pthreads_h94d23a6_0 - - libpng=1.6.47=h943b412_0 - - libpq=17.4=h27ae623_1 - - librdkit=2025.03.2=h84b0b3c_0 - - libsodium=1.0.20=h4ab18f5_0 - - libsqlite=3.49.1=hee588c1_1 - - libssh2=1.11.1=hcf80075_0 - - libstdcxx=14.2.0=h8f9b012_2 - - libstdcxx-ng=14.2.0=h4852527_2 - - libtiff=4.7.0=hd9ff511_3 - - libuuid=2.38.1=h0b41bf4_0 - - libwebp-base=1.5.0=h851e524_0 - - libxcb=1.17.0=h8a09558_0 - - libxcrypt=4.4.36=hd590300_1 - - libxml2=2.13.7=h8d12d68_0 - - libxslt=1.1.39=h76b75d6_0 - - libzip=1.11.2=h6991a6a_0 - - libzlib=1.3.1=hb9d3cd8_2 - - llvm-openmp=20.1.4=h024ca30_0 - - lxml=5.4.0=py310h490dddc_0 - - lz4-c=1.10.0=h5888daf_1 - - markupsafe=3.0.2=py310h89163eb_1 - - matplotlib-base=3.10.1=py310h68603db_0 - - matplotlib-inline=0.1.7=pyhd8ed1ab_1 - - mda-xdrlib=0.2.0=pyhd8ed1ab_1 - - mdaencore=1.0.0=py310hf462985_2 - - mdahole2-base=0.5.0=pyhd8ed1ab_1 - - mdtraj=1.10.3=py310h4cdbd58_0 - - mistune=3.1.3=pyh29332c3_0 - - mkl=2024.2.2=ha957f24_16 - - mmtf-python=1.1.3=pyhd8ed1ab_0 - - modeller>=10.6,<11 - - mpfr=4.2.1=h90cbb55_3 - - mrcfile=1.5.4=pyhd8ed1ab_0 - - msgpack-python=1.1.0=py310h3788b33_0 - - munkres=1.1.4=pyh9f0ad1d_0 - - nbclient=0.10.2=pyhd8ed1ab_0 - - nbconvert-core=7.16.6=pyh29332c3_0 - - nbformat=5.10.4=pyhd8ed1ab_1 - - ncurses=6.5=h2d0b736_3 - - nest-asyncio=1.6.0=pyhd8ed1ab_1 - - netcdf-fortran=4.6.1=nompi_h22f9119_108 - - netcdf4=1.7.2=nompi_py310h9f0ad05_101 - - networkx=3.4.2=pyh267e887_2 - - nglview=3.1.4=pyh620948e_0 - - notebook=7.4.1=pyhd8ed1ab_0 - - notebook-shim=0.2.4=pyhd8ed1ab_1 - - numexpr=2.10.2=mkl_py310ha004913_0 - - ocl-icd=2.3.3=hb9d3cd8_0 - - ocl-icd-system=1.0.0=1 - - opencl-headers=2024.10.24=h5888daf_0 - - openff-amber-ff-ports=0.0.4=pyhca7485f_0 - - openff-forcefields=2024.09.0=pyhff2d567_0 - - openff-interchange=0.4.2=pyhd8ed1ab_2 - - openff-interchange-base=0.4.2=pyhd8ed1ab_2 - - openff-toolkit=0.16.9=pyhd8ed1ab_0 - - openff-toolkit-base=0.16.9=pyhd8ed1ab_0 - - openff-units=0.3.0=pyhd8ed1ab_1 - - openff-utilities=0.1.15=pyhd8ed1ab_0 - - openjpeg=2.5.3=h5fbd93e_0 - - openldap=2.6.9=he970967_0 - - openmm=8.2.0=py310h30bdd6a_2 - - openmmforcefields=0.14.2=pyhd8ed1ab_0 - - openssl=3.5.1=h7b32b05_0 - - overrides=7.7.0=pyhd8ed1ab_1 - - pandas=2.2.3=py310h5eaa309_3 - - pandocfilters=1.5.0=pyhd8ed1ab_0 - - panedr=0.8.0=pyhd8ed1ab_1 - - parmed=4.3.0=py310h78e4988_2 - - parso=0.8.4=pyhd8ed1ab_1 - - patsy=1.0.1=pyhd8ed1ab_1 - - pcre2=10.45=hc749103_0 - - perl=5.32.1=7_hd590300_perl5 - - pexpect=4.9.0=pyhd8ed1ab_1 - - pickleshare=0.7.5=pyhd8ed1ab_1004 - - pillow=11.1.0=py310h7e6dc6c_0 - - pint=0.24.4=pyhd8ed1ab_1 - - pip=25.1.1=pyh8b19718_0 - - pixman=0.46.0=h29eaf8c_0 - - pkgutil-resolve-name=1.3.10=pyhd8ed1ab_2 - - platformdirs=4.3.7=pyh29332c3_0 - - plumed=2.9.2=nompi_hda5c677_100 - - prometheus_client=0.21.1=pyhd8ed1ab_0 - - prompt-toolkit=3.0.51=pyha770c72_0 - - psutil=7.0.0=py310ha75aee5_0 - - pthread-stubs=0.4=hb9d3cd8_1002 - - ptyprocess=0.7.0=pyhd8ed1ab_1 - - pure_eval=0.2.3=pyhd8ed1ab_1 - - py-cpuinfo=9.0.0=pyhd8ed1ab_1 - - pycairo=1.28.0=py310h516fd05_0 - - pycparser=2.22=pyh29332c3_1 - - pydantic=2.11.3=pyh3cfb1c2_0 - - pydantic-core=2.33.1=py310hc1293b2_0 - - pyedr=0.8.0=pyhd8ed1ab_1 - - pygments=2.19.1=pyhd8ed1ab_0 - - pysocks=1.7.1=pyha55dd90_7 - - pytables=3.10.1=py310h1affd9f_4 - - python-constraint=1.4.0=pyhff2d567_1 - - python-dateutil=2.9.0.post0=pyhff2d567_1 - - python-fastjsonschema=2.21.1=pyhd8ed1ab_0 - - python-json-logger=2.0.7=pyhd8ed1ab_0 - - python-tzdata=2025.2=pyhd8ed1ab_0 - - pytng=0.3.3=py310he95d665_1 - - pyyaml=6.0.2=py310h89163eb_2 - - pyzmq=26.4.0=py310h71f11fc_0 - - qhull=2020.2=h434a139_5 - - rdkit=2025.03.2=py310hcd13295_0 - - readline=8.2=h8c095d6_2 - - referencing=0.36.2=pyh29332c3_0 - - reportlab=4.4.0=py310ha75aee5_0 - - requests=2.32.3=pyhd8ed1ab_1 - - rfc3339-validator=0.1.4=pyhd8ed1ab_1 - - rfc3986-validator=0.1.1=pyh9f0ad1d_0 - - rlpycairo=0.2.0=pyhd8ed1ab_0 - - rpds-py=0.24.0=py310hc1293b2_0 - - scikit-learn=1.6.1=py310h27f47ee_0 - - seaborn=0.13.2=hd8ed1ab_3 - - seaborn-base=0.13.2=pyhd8ed1ab_3 - - send2trash=1.8.3=pyh0d859eb_1 - - setuptools=75.8.2=pyhff2d567_0 - - six=1.17.0=pyhd8ed1ab_0 - - smirnoff99frosst=1.1.0=pyh44b312d_0 - - snappy=1.2.1=h8bd8927_1 - - sniffio=1.3.1=pyhd8ed1ab_1 - - soupsieve=2.5=pyhd8ed1ab_1 - - sqlalchemy=2.0.40=py310ha75aee5_0 - - sqlite=3.49.1=h9eae976_1 - - stack_data=0.6.3=pyhd8ed1ab_1 - - statsmodels=0.14.4=py310hf462985_0 - - sysroot_linux-64=2.17=h0157908_18 - - tbb=2021.13.0=hceb3a55_1 - - terminado=0.18.1=pyh0d859eb_0 - - threadpoolctl=3.6.0=pyhecae5ae_0 - - tidynamics=1.1.2=pyhd8ed1ab_0 - - tinycss2=1.4.0=pyhd8ed1ab_0 - - tinydb=4.8.2=pyhd8ed1ab_1 - - tk=8.6.13=noxft_h4845f30_101 - - tomli=2.2.1=pyhd8ed1ab_1 - - tornado=6.4.2=py310ha75aee5_0 - - tqdm=4.67.1=pyhd8ed1ab_1 - - traitlets=5.14.3=pyhd8ed1ab_1 - - types-python-dateutil=2.9.0.20241206=pyhd8ed1ab_0 - - typing-extensions=4.13.2=h0e9735f_0 - - typing-inspection=0.4.0=pyhd8ed1ab_0 - - typing_extensions=4.13.2=pyh29332c3_0 - - typing_utils=0.1.0=pyhd8ed1ab_1 - - unicodedata2=16.0.0=py310ha75aee5_0 - - uri-template=1.3.0=pyhd8ed1ab_1 - - validators=0.35.0=pyhd8ed1ab_0 - - wcwidth=0.2.13=pyhd8ed1ab_1 - - webcolors=24.11.1=pyhd8ed1ab_0 - - webencodings=0.5.1=pyhd8ed1ab_3 - - websocket-client=1.8.0=pyhd8ed1ab_1 - - wheel=0.45.1=pyhd8ed1ab_1 - - widgetsnbextension=4.0.14=pyhd8ed1ab_0 - - xmltodict=0.14.2=pyhd8ed1ab_1 - - xorg-libice=1.1.2=hb9d3cd8_0 - - xorg-libsm=1.2.6=he73a12e_0 - - xorg-libx11=1.8.12=h4f16b4b_0 - - xorg-libxau=1.0.12=hb9d3cd8_0 - - xorg-libxdmcp=1.1.5=hb9d3cd8_0 - - xorg-libxext=1.3.6=hb9d3cd8_0 - - xorg-libxrender=0.9.12=hb9d3cd8_0 - - xorg-libxt=1.3.1=hb9d3cd8_0 - - xz=5.6.4=hbcc6ac9_0 - - xz-gpl-tools=5.6.4=hbcc6ac9_0 - - xz-tools=5.6.4=hb9d3cd8_0 - - yaml=0.2.5=h7f98852_2 - - zeromq=4.3.5=h3b0a872_7 - - zipp=3.21.0=pyhd8ed1ab_1 - - zlib=1.3.1=hb9d3cd8_2 - - zlib-ng=2.2.4=h7955e40_0 - - zstandard=0.23.0=py310ha75aee5_2 - - zstd=1.5.7=hb8e6e7a_2 - # --- move core scientific stack to conda to avoid conflicts --- - - numpy=1.26.* - - scipy=1.13.* - - mdanalysis=2.9.* - - umap-learn=0.5.* - - numba=0.61.* - - llvmlite=0.44.* - - pynndescent=0.5.13 - - pip: - - aiohappyeyeballs==2.6.1 - - aiohttp==3.12.13 - - aiosignal==1.3.2 - - amberutils==21.0 - - antlr4-python3-runtime==4.9.3 - - async-timeout==5.0.1 - - biopython==1.84 - - biotite==1.1.0 - - biotraj==1.2.2 - - boltz==2.1.1 - - charset-normalizer==3.4.1 - - chembl-structure-pipeline==1.2.2 - - click==8.1.7 - - coverage==7.6.12 - - cuequivariance==0.5.0 - - cuequivariance-ops-cu12==0.5.0 - - cuequivariance-ops-torch-cu12==0.5.0 - - cuequivariance-torch==0.5.0 - - deeptime==0.4.4 - - dm-tree==0.1.8 - - docker-pycreds==0.4.0 - - edgembar==3.0 - - einops==0.8.0 - - einx==0.3.0 - - fairscale==0.4.13 - - fonttools==4.56.0 - - frozendict==2.4.6 - - frozenlist==1.7.0 - - fsspec==2025.5.1 - - gemmi==0.6.5 - - gitdb==4.0.12 - - gitpython==3.1.44 - - gpcrmining==0.2.0 - - gromacswrapper==0.9.1 - - hydra-core==1.3.2 - - ihm==2.6 - - iniconfig==2.1.0 - - loguru==0.7.3 - - mashumaro==3.14 - - mdahole2==0.5.0 - - mmpbsa-py==16.0 - - modelcif==1.2 - - mpmath==1.3.0 - - multidict==6.5.0 - - numkit==1.3.1 - - nvidia-cublas-cu12==12.6.4.1 - - nvidia-cuda-cupti-cu12==12.6.80 - - nvidia-cuda-nvrtc-cu12==12.6.77 - - nvidia-cuda-runtime-cu12==12.6.77 - - nvidia-cudnn-cu12==9.5.1.17 - - nvidia-cufft-cu12==11.3.0.4 - - nvidia-cufile-cu12==1.11.1.6 - - nvidia-curand-cu12==10.3.7.77 - - nvidia-cusolver-cu12==11.7.1.2 - - nvidia-cusparse-cu12==12.5.4.2 - - nvidia-cusparselt-cu12==0.6.3 - - nvidia-ml-py==12.575.51 - - nvidia-nccl-cu12==2.26.2 - - nvidia-nvjitlink-cu12==12.6.85 - - nvidia-nvtx-cu12==12.6.77 - - omegaconf==2.3.0 - - packaging==24.2 - - packmol-memgen==2025.1.29 - - pathsimanalysis==1.2.0 - - pdb4amber==22.0 - - pensa==0.6.0 - - pluggy==1.6.0 - - propcache==0.3.2 - - protobuf==5.29.5 - - py==1.11.0 - - pymemento==0.1.1 - - pymsmt==22.0 - - pynvml==12.0.0 - - pyparsing==3.2.1 - - pytest==8.4.1 - - pytest-cov==3.0.0 - - pytorch-lightning==2.5.0 - - pytraj==2.0.6 - - pytz==2025.1 - - sander==22.0 - - sentry-sdk==2.30.0 - - setproctitle==1.3.6 - - smmap==5.0.2 - - sympy==1.14.0 - - torch==2.7.1 - - torchmetrics==1.7.3 - - triton==3.3.1 - - types-requests==2.32.4.20250611 - - tzdata==2025.1 - - urllib3==2.3.0 - - wandb==0.18.7 - - waterdynamics==1.2.0 - - yarl==1.20.1 + - python=3.10.0 + - plumed + - modeller + - mdaencore + - ipykernel + - numpy + - cython + - openmm + - openmmforcefields + - numba +prefix: /biggin/b211/reub0138/miniforge3/envs/fepa diff --git a/requirements_old.yml b/requirements_old.yml new file mode 100644 index 0000000..50eca8f --- /dev/null +++ b/requirements_old.yml @@ -0,0 +1,424 @@ +name: fepa_env +channels: + - salilab + - conda-forge + - defaults +dependencies: + - _libgcc_mutex=0.1=conda_forge + - _openmp_mutex=4.5=3_kmp_llvm + - ambertools=24.8=cuda_None_nompi_py310h834fefc_101 + - annotated-types=0.7.0=pyhd8ed1ab_1 + - anyio=4.9.0=pyh29332c3_0 + - argon2-cffi=23.1.0=pyhd8ed1ab_1 + - argon2-cffi-bindings=21.2.0=py310ha75aee5_5 + - arpack=3.9.1=nompi_hf03ea27_102 + - arrow=1.3.0=pyhd8ed1ab_1 + - asttokens=3.0.0=pyhd8ed1ab_1 + - async-lru=2.0.5=pyh29332c3_0 + - attrs=25.3.0=pyh71513ae_0 + - babel=2.17.0=pyhd8ed1ab_0 + - beautifulsoup4=4.13.4=pyha770c72_0 + - bleach=6.2.0=pyh29332c3_4 + - bleach-with-css=6.2.0=h82add2a_4 + - blosc=1.21.6=he440d0b_1 + - brotli=1.1.0=hb9d3cd8_2 + - brotli-bin=1.1.0=hb9d3cd8_2 + - brotli-python=1.1.0=py310hf71b8c6_2 + - bson=0.5.10=pyhd8ed1ab_0 + - bzip2=1.0.8=h4bc722e_7 + - c-ares=1.34.5=hb9d3cd8_0 + - c-blosc2=2.15.2=h3122c55_1 + - ca-certificates=2025.6.15=hbd8a1cb_0 + - cached-property=1.5.2=hd8ed1ab_1 + - cached_property=1.5.2=pyha770c72_1 + - cachetools=5.5.2=pyhd8ed1ab_0 + - cairo=1.18.4=h3394656_0 + - certifi=2025.6.15=pyhd8ed1ab_0 + - cffi=1.17.1=py310h8deb56e_0 + - cftime=1.6.4=py310hf462985_1 + - chardet=5.2.0=pyhd8ed1ab_3 + - colorama=0.4.6=pyhd8ed1ab_1 + - comm=0.2.2=pyhd8ed1ab_1 + - cudatoolkit=11.8.0=h4ba93d1_13 + - cycler=0.12.1=pyhd8ed1ab_1 + - cyrus-sasl=2.1.27=h54b06d7_7 + - cython=3.0.12=py310had8cdd9_0 + - debugpy=1.8.14=py310hf71b8c6_0 + - decorator=5.2.1=pyhd8ed1ab_0 + - defusedxml=0.7.1=pyhd8ed1ab_0 + - exceptiongroup=1.2.2=pyhd8ed1ab_1 + - executing=2.2.0=pyhd8ed1ab_0 + - fasteners=0.19=pyhd8ed1ab_1 + - fftw=3.3.10=nompi_hf1063bd_110 + - filelock=3.18.0=pyhd8ed1ab_0 + - flexcache=0.3=pyhd8ed1ab_1 + - flexparser=0.4=pyhd8ed1ab_1 + - font-ttf-dejavu-sans-mono=2.37=hab24e00_0 + - font-ttf-inconsolata=3.000=h77eed37_0 + - font-ttf-source-code-pro=2.038=h77eed37_0 + - font-ttf-ubuntu=0.83=h77eed37_3 + - fontconfig=2.15.0=h7e30c49_1 + - fonts-conda-ecosystem=1=0 + - fonts-conda-forge=1=0 + - fqdn=1.5.1=pyhd8ed1ab_1 + - freetype=2.13.3=ha770c72_1 + - freetype-py=2.3.0=pyhd8ed1ab_0 + - gawk=5.3.1=hcd3d067_0 + - gmp=6.3.0=hac33072_2 + - greenlet=3.2.1=py310hf71b8c6_0 + - griddataformats=1.0.2=pyhd8ed1ab_1 + - gsd=3.4.2=py310ha75aee5_0 + - gsl=2.7=he838d99_0 + - h11=0.16.0=pyhd8ed1ab_0 + - h2=4.2.0=pyhd8ed1ab_0 + - h5py=3.13.0=nompi_py310h60e0fe6_100 + - hdf4=4.2.15=h2a13503_7 + - hdf5=1.14.3=nompi_h2d575fe_109 + - hdf5-1107=1.10.7=1 + - hpack=4.1.0=pyhd8ed1ab_0 + - httpcore=1.0.9=pyh29332c3_0 + - httpx=0.28.1=pyhd8ed1ab_0 + - hyperframe=6.1.0=pyhd8ed1ab_0 + - icu=75.1=he02047a_0 + - idna=3.10=pyhd8ed1ab_1 + - importlib-metadata=8.6.1=pyha770c72_0 + - importlib_resources=6.5.2=pyhd8ed1ab_0 + - ipykernel=6.29.5=pyh3099207_0 + - ipython=8.36.0=pyh907856f_0 + - ipywidgets=8.1.7=pyhd8ed1ab_0 + - isoduration=20.11.0=pyhd8ed1ab_1 + - jedi=0.19.2=pyhd8ed1ab_1 + - jinja2=3.1.6=pyhd8ed1ab_0 + - joblib=1.4.2=pyhd8ed1ab_1 + - json5=0.12.0=pyhd8ed1ab_0 + - jsonpointer=3.0.0=py310hff52083_1 + - jsonschema=4.23.0=pyhd8ed1ab_1 + - jsonschema-specifications=2025.4.1=pyh29332c3_0 + - jsonschema-with-format-nongpl=4.23.0=hd8ed1ab_1 + - jupyter-lsp=2.2.5=pyhd8ed1ab_1 + - jupyter_client=8.6.3=pyhd8ed1ab_1 + - jupyter_core=5.7.2=pyh31011fe_1 + - jupyter_events=0.12.0=pyh29332c3_0 + - jupyter_server=2.15.0=pyhd8ed1ab_0 + - jupyter_server_terminals=0.5.3=pyhd8ed1ab_1 + - jupyterlab=4.4.1=pyhd8ed1ab_0 + - jupyterlab_pygments=0.3.0=pyhd8ed1ab_2 + - jupyterlab_server=2.27.3=pyhd8ed1ab_1 + - jupyterlab_widgets=3.0.15=pyhd8ed1ab_0 + - kernel-headers_linux-64=3.10.0=he073ed8_18 + - keyutils=1.6.1=h166bdaf_0 + - krb5=1.21.3=h659f571_0 + - lcms2=2.17=h717163a_0 + - ld_impl_linux-64=2.43=h712a8e2_4 + - lerc=4.0.0=h0aef613_1 + - libaec=1.1.3=h59595ed_0 + - libasprintf=0.23.1=h8e693c7_0 + - libblas=3.9.0=31_hfdb39a5_mkl + - libboost=1.86.0=h6c02f8c_3 + - libboost-python=1.86.0=py310ha2bacc8_3 + - libbrotlicommon=1.1.0=hb9d3cd8_2 + - libbrotlidec=1.1.0=hb9d3cd8_2 + - libbrotlienc=1.1.0=hb9d3cd8_2 + - libcblas=3.9.0=31_h372d94f_mkl + - libcurl=8.13.0=h332b0f4_0 + - libdeflate=1.23=h86f0d12_0 + - libedit=3.1.20250104=pl5321h7949ede_0 + - libev=4.33=hd590300_2 + - libexpat=2.7.0=h5888daf_0 + - libffi=3.4.6=h2dba641_0 + - libfreetype=2.13.3=ha770c72_1 + - libfreetype6=2.13.3=h48d6fc4_1 + - libgcc=14.2.0=h767d61c_2 + - libgcc-ng=14.2.0=h69a702a_2 + - libgettextpo=0.23.1=h5888daf_0 + - libgfortran=14.2.0=h69a702a_2 + - libgfortran-ng=14.2.0=h69a702a_2 + - libgfortran5=14.2.0=hf1ad2bd_2 + - libglib=2.84.1=h3618099_1 + - libgomp=14.2.0=h767d61c_2 + - libhwloc=2.11.2=default_h0d58e46_1001 + - libiconv=1.18=h4ce23a2_1 + - libjpeg-turbo=3.1.0=hb9d3cd8_0 + - liblapack=3.9.0=31_hc41d3b0_mkl + - liblzma=5.6.4=hb9d3cd8_0 + - liblzma-devel=5.6.4=hb9d3cd8_0 + - libnetcdf=4.9.2=nompi_h00e09a9_116 + - libnghttp2=1.64.0=h161d5f1_0 + - libnsl=2.0.1=hd590300_0 + - libntlm=1.8=hb9d3cd8_0 + - libopenblas=0.3.29=pthreads_h94d23a6_0 + - libpng=1.6.47=h943b412_0 + - libpq=17.4=h27ae623_1 + - librdkit=2025.03.2=h84b0b3c_0 + - libsodium=1.0.20=h4ab18f5_0 + - libsqlite=3.49.1=hee588c1_1 + - libssh2=1.11.1=hcf80075_0 + - libstdcxx=14.2.0=h8f9b012_2 + - libstdcxx-ng=14.2.0=h4852527_2 + - libtiff=4.7.0=hd9ff511_3 + - libuuid=2.38.1=h0b41bf4_0 + - libwebp-base=1.5.0=h851e524_0 + - libxcb=1.17.0=h8a09558_0 + - libxcrypt=4.4.36=hd590300_1 + - libxml2=2.13.7=h8d12d68_0 + - libxslt=1.1.39=h76b75d6_0 + - libzip=1.11.2=h6991a6a_0 + - libzlib=1.3.1=hb9d3cd8_2 + - llvm-openmp=20.1.4=h024ca30_0 + - lxml=5.4.0=py310h490dddc_0 + - lz4-c=1.10.0=h5888daf_1 + - markupsafe=3.0.2=py310h89163eb_1 + - matplotlib-base=3.10.1=py310h68603db_0 + - matplotlib-inline=0.1.7=pyhd8ed1ab_1 + - mda-xdrlib=0.2.0=pyhd8ed1ab_1 + - mdaencore=1.0.0=py310hf462985_2 + - mdahole2-base=0.5.0=pyhd8ed1ab_1 + - mdtraj=1.10.3=py310h4cdbd58_0 + - mistune=3.1.3=pyh29332c3_0 + - mkl=2024.2.2=ha957f24_16 + - mmtf-python=1.1.3=pyhd8ed1ab_0 + - modeller>=10.6,<11 + - mpfr=4.2.1=h90cbb55_3 + - mrcfile=1.5.4=pyhd8ed1ab_0 + - msgpack-python=1.1.0=py310h3788b33_0 + - munkres=1.1.4=pyh9f0ad1d_0 + - nbclient=0.10.2=pyhd8ed1ab_0 + - nbconvert-core=7.16.6=pyh29332c3_0 + - nbformat=5.10.4=pyhd8ed1ab_1 + - ncurses=6.5=h2d0b736_3 + - nest-asyncio=1.6.0=pyhd8ed1ab_1 + - netcdf-fortran=4.6.1=nompi_h22f9119_108 + - netcdf4=1.7.2=nompi_py310h9f0ad05_101 + - networkx=3.4.2=pyh267e887_2 + - nglview=3.1.4=pyh620948e_0 + - notebook=7.4.1=pyhd8ed1ab_0 + - notebook-shim=0.2.4=pyhd8ed1ab_1 + - numexpr=2.10.2=mkl_py310ha004913_0 + - ocl-icd=2.3.3=hb9d3cd8_0 + - ocl-icd-system=1.0.0=1 + - opencl-headers=2024.10.24=h5888daf_0 + - openff-amber-ff-ports=0.0.4=pyhca7485f_0 + - openff-forcefields=2024.09.0=pyhff2d567_0 + - openff-interchange=0.4.2=pyhd8ed1ab_2 + - openff-interchange-base=0.4.2=pyhd8ed1ab_2 + - openff-toolkit=0.16.9=pyhd8ed1ab_0 + - openff-toolkit-base=0.16.9=pyhd8ed1ab_0 + - openff-units=0.3.0=pyhd8ed1ab_1 + - openff-utilities=0.1.15=pyhd8ed1ab_0 + - openjpeg=2.5.3=h5fbd93e_0 + - openldap=2.6.9=he970967_0 + - openmm=8.2.0=py310h30bdd6a_2 + - openmmforcefields=0.14.2=pyhd8ed1ab_0 + - openssl=3.5.1=h7b32b05_0 + - overrides=7.7.0=pyhd8ed1ab_1 + - pandas=2.2.3=py310h5eaa309_3 + - pandocfilters=1.5.0=pyhd8ed1ab_0 + - panedr=0.8.0=pyhd8ed1ab_1 + - parmed=4.3.0=py310h78e4988_2 + - parso=0.8.4=pyhd8ed1ab_1 + - patsy=1.0.1=pyhd8ed1ab_1 + - pcre2=10.45=hc749103_0 + - perl=5.32.1=7_hd590300_perl5 + - pexpect=4.9.0=pyhd8ed1ab_1 + - pickleshare=0.7.5=pyhd8ed1ab_1004 + - pillow=11.1.0=py310h7e6dc6c_0 + - pint=0.24.4=pyhd8ed1ab_1 + - pip=25.1.1=pyh8b19718_0 + - pixman=0.46.0=h29eaf8c_0 + - pkgutil-resolve-name=1.3.10=pyhd8ed1ab_2 + - platformdirs=4.3.7=pyh29332c3_0 + - plumed=2.9.2=nompi_hda5c677_100 + - prometheus_client=0.21.1=pyhd8ed1ab_0 + - prompt-toolkit=3.0.51=pyha770c72_0 + - psutil=7.0.0=py310ha75aee5_0 + - pthread-stubs=0.4=hb9d3cd8_1002 + - ptyprocess=0.7.0=pyhd8ed1ab_1 + - pure_eval=0.2.3=pyhd8ed1ab_1 + - py-cpuinfo=9.0.0=pyhd8ed1ab_1 + - pycairo=1.28.0=py310h516fd05_0 + - pycparser=2.22=pyh29332c3_1 + - pydantic=2.11.3=pyh3cfb1c2_0 + - pydantic-core=2.33.1=py310hc1293b2_0 + - pyedr=0.8.0=pyhd8ed1ab_1 + - pygments=2.19.1=pyhd8ed1ab_0 + - pysocks=1.7.1=pyha55dd90_7 + - pytables=3.10.1=py310h1affd9f_4 + - python-constraint=1.4.0=pyhff2d567_1 + - python-dateutil=2.9.0.post0=pyhff2d567_1 + - python-fastjsonschema=2.21.1=pyhd8ed1ab_0 + - python-json-logger=2.0.7=pyhd8ed1ab_0 + - python-tzdata=2025.2=pyhd8ed1ab_0 + - pytng=0.3.3=py310he95d665_1 + - pyyaml=6.0.2=py310h89163eb_2 + - pyzmq=26.4.0=py310h71f11fc_0 + - qhull=2020.2=h434a139_5 + - rdkit=2025.03.2=py310hcd13295_0 + - readline=8.2=h8c095d6_2 + - referencing=0.36.2=pyh29332c3_0 + - reportlab=4.4.0=py310ha75aee5_0 + - requests=2.32.3=pyhd8ed1ab_1 + - rfc3339-validator=0.1.4=pyhd8ed1ab_1 + - rfc3986-validator=0.1.1=pyh9f0ad1d_0 + - rlpycairo=0.2.0=pyhd8ed1ab_0 + - rpds-py=0.24.0=py310hc1293b2_0 + - scikit-learn=1.6.1=py310h27f47ee_0 + - seaborn=0.13.2=hd8ed1ab_3 + - seaborn-base=0.13.2=pyhd8ed1ab_3 + - send2trash=1.8.3=pyh0d859eb_1 + - setuptools=75.8.2=pyhff2d567_0 + - six=1.17.0=pyhd8ed1ab_0 + - smirnoff99frosst=1.1.0=pyh44b312d_0 + - snappy=1.2.1=h8bd8927_1 + - sniffio=1.3.1=pyhd8ed1ab_1 + - soupsieve=2.5=pyhd8ed1ab_1 + - sqlalchemy=2.0.40=py310ha75aee5_0 + - sqlite=3.49.1=h9eae976_1 + - stack_data=0.6.3=pyhd8ed1ab_1 + - statsmodels=0.14.4=py310hf462985_0 + - sysroot_linux-64=2.17=h0157908_18 + - tbb=2021.13.0=hceb3a55_1 + - terminado=0.18.1=pyh0d859eb_0 + - threadpoolctl=3.6.0=pyhecae5ae_0 + - tidynamics=1.1.2=pyhd8ed1ab_0 + - tinycss2=1.4.0=pyhd8ed1ab_0 + - tinydb=4.8.2=pyhd8ed1ab_1 + - tk=8.6.13=noxft_h4845f30_101 + - tomli=2.2.1=pyhd8ed1ab_1 + - tornado=6.4.2=py310ha75aee5_0 + - tqdm=4.67.1=pyhd8ed1ab_1 + - traitlets=5.14.3=pyhd8ed1ab_1 + - types-python-dateutil=2.9.0.20241206=pyhd8ed1ab_0 + - typing-extensions=4.13.2=h0e9735f_0 + - typing-inspection=0.4.0=pyhd8ed1ab_0 + - typing_extensions=4.13.2=pyh29332c3_0 + - typing_utils=0.1.0=pyhd8ed1ab_1 + - unicodedata2=16.0.0=py310ha75aee5_0 + - uri-template=1.3.0=pyhd8ed1ab_1 + - validators=0.35.0=pyhd8ed1ab_0 + - wcwidth=0.2.13=pyhd8ed1ab_1 + - webcolors=24.11.1=pyhd8ed1ab_0 + - webencodings=0.5.1=pyhd8ed1ab_3 + - websocket-client=1.8.0=pyhd8ed1ab_1 + - wheel=0.45.1=pyhd8ed1ab_1 + - widgetsnbextension=4.0.14=pyhd8ed1ab_0 + - xmltodict=0.14.2=pyhd8ed1ab_1 + - xorg-libice=1.1.2=hb9d3cd8_0 + - xorg-libsm=1.2.6=he73a12e_0 + - xorg-libx11=1.8.12=h4f16b4b_0 + - xorg-libxau=1.0.12=hb9d3cd8_0 + - xorg-libxdmcp=1.1.5=hb9d3cd8_0 + - xorg-libxext=1.3.6=hb9d3cd8_0 + - xorg-libxrender=0.9.12=hb9d3cd8_0 + - xorg-libxt=1.3.1=hb9d3cd8_0 + - xz=5.6.4=hbcc6ac9_0 + - xz-gpl-tools=5.6.4=hbcc6ac9_0 + - xz-tools=5.6.4=hb9d3cd8_0 + - yaml=0.2.5=h7f98852_2 + - zeromq=4.3.5=h3b0a872_7 + - zipp=3.21.0=pyhd8ed1ab_1 + - zlib=1.3.1=hb9d3cd8_2 + - zlib-ng=2.2.4=h7955e40_0 + - zstandard=0.23.0=py310ha75aee5_2 + - zstd=1.5.7=hb8e6e7a_2 + # --- move core scientific stack to conda to avoid conflicts --- + - numpy=1.26.* + - scipy=1.13.* + - mdanalysis=2.9.* + - umap-learn=0.5.* + - numba=0.61.* + - llvmlite=0.44.* + - pynndescent=0.5.13 + - pip: + - aiohappyeyeballs==2.6.1 + - aiohttp==3.12.13 + - aiosignal==1.3.2 + - amberutils==21.0 + - antlr4-python3-runtime==4.9.3 + - async-timeout==5.0.1 + - biopython==1.84 + - biotite==1.1.0 + - biotraj==1.2.2 + - boltz==2.1.1 + - charset-normalizer==3.4.1 + - chembl-structure-pipeline==1.2.2 + - click==8.1.7 + - coverage==7.6.12 + - cuequivariance==0.5.0 + - cuequivariance-ops-cu12==0.5.0 + - cuequivariance-ops-torch-cu12==0.5.0 + - cuequivariance-torch==0.5.0 + - deeptime==0.4.4 + - dm-tree==0.1.8 + - docker-pycreds==0.4.0 + - edgembar==3.0 + - einops==0.8.0 + - einx==0.3.0 + - fairscale==0.4.13 + - fonttools==4.56.0 + - frozendict==2.4.6 + - frozenlist==1.7.0 + - fsspec==2025.5.1 + - gemmi==0.6.5 + - gitdb==4.0.12 + - gitpython==3.1.44 + - gpcrmining==0.2.0 + - gromacswrapper==0.9.1 + - hydra-core==1.3.2 + - ihm==2.6 + - iniconfig==2.1.0 + - loguru==0.7.3 + - mashumaro==3.14 + - mdahole2==0.5.0 + - mmpbsa-py==16.0 + - modelcif==1.2 + - mpmath==1.3.0 + - multidict==6.5.0 + - numkit==1.3.1 + - nvidia-cublas-cu12==12.6.4.1 + - nvidia-cuda-cupti-cu12==12.6.80 + - nvidia-cuda-nvrtc-cu12==12.6.77 + - nvidia-cuda-runtime-cu12==12.6.77 + - nvidia-cudnn-cu12==9.5.1.17 + - nvidia-cufft-cu12==11.3.0.4 + - nvidia-cufile-cu12==1.11.1.6 + - nvidia-curand-cu12==10.3.7.77 + - nvidia-cusolver-cu12==11.7.1.2 + - nvidia-cusparse-cu12==12.5.4.2 + - nvidia-cusparselt-cu12==0.6.3 + - nvidia-ml-py==12.575.51 + - nvidia-nccl-cu12==2.26.2 + - nvidia-nvjitlink-cu12==12.6.85 + - nvidia-nvtx-cu12==12.6.77 + - omegaconf==2.3.0 + - packaging==24.2 + - packmol-memgen==2025.1.29 + - pathsimanalysis==1.2.0 + - pdb4amber==22.0 + - pensa==0.6.0 + - pluggy==1.6.0 + - propcache==0.3.2 + - protobuf==5.29.5 + - py==1.11.0 + - pymemento==0.1.1 + - pymsmt==22.0 + - pynvml==12.0.0 + - pyparsing==3.2.1 + - pytest==8.4.1 + - pytest-cov==3.0.0 + - pytorch-lightning==2.5.0 + - pytraj==2.0.6 + - pytz==2025.1 + - sander==22.0 + - sentry-sdk==2.30.0 + - setproctitle==1.3.6 + - smmap==5.0.2 + - sympy==1.14.0 + - torch==2.7.1 + - torchmetrics==1.7.3 + - triton==3.3.1 + - types-requests==2.32.4.20250611 + - tzdata==2025.1 + - urllib3==2.3.0 + - wandb==0.18.7 + - waterdynamics==1.2.0 + - yarl==1.20.1 From 95cf1ec029622bbabef8dedbde31e92f9e3bf04f Mon Sep 17 00:00:00 2001 From: nithishwer Date: Wed, 12 Nov 2025 21:07:29 +0000 Subject: [PATCH 16/40] fix: edit python req --- requirements.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.yml b/requirements.yml index 80eafc4..93d3c05 100644 --- a/requirements.yml +++ b/requirements.yml @@ -4,7 +4,7 @@ channels: - conda-forge - omnia dependencies: - - python=3.10.0 + - python>=3.10 - plumed - modeller - mdaencore From 5615b36c2c3d59e0eec69f6aa32495b8d9b3e6be Mon Sep 17 00:00:00 2001 From: nithishwer Date: Wed, 12 Nov 2025 21:35:18 +0000 Subject: [PATCH 17/40] fix: better gh workflow + new how-to-guide --- .github/workflows/test.yml | 21 ++- docs/how-to-guides/guide-1.md | 104 ++++++++++++- site/how-to-guides/guide-1/index.html | 213 +++++++++++++++++++++++++- site/search.json | 2 +- 4 files changed, 321 insertions(+), 19 deletions(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 04ed04a..9d70c02 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -14,26 +14,25 @@ jobs: - name: Checkout code uses: actions/checkout@v4 - - name: Set up Conda - uses: conda-incubator/setup-miniconda@v2 - with: - python-version: 3.10 - activate-environment: fepa_env - environment-file: requirements.yml - auto-update-conda: true - - name: Configure Conda channels run: | conda config --add channels salilab conda config --add channels conda-forge conda config --add channels defaults - - name: Install editable package + - name: Set up Conda environment + uses: conda-incubator/setup-miniconda@v3 + with: + python-version: 3.10 + environment-file: requirements.yml + activate-environment: fepa_env + auto-update-conda: true + use-mamba: true + + - name: Install package in editable mode run: | pip install --no-cache-dir -e . - name: Run tests run: | pytest -v --maxfail=1 --disable-warnings - - diff --git a/docs/how-to-guides/guide-1.md b/docs/how-to-guides/guide-1.md index f2a3cc1..3f1fdf2 100644 --- a/docs/how-to-guides/guide-1.md +++ b/docs/how-to-guides/guide-1.md @@ -1,6 +1,106 @@ --- icon: lucide/compass -title: How-to-Guide +title: Binding Pocket Resolvation Analysis --- -Put the example analysis in the how to guide section here \ No newline at end of file +--- +icon: lucide/compass +title: Binding Pocket Resolvation Analysis +--- + +Visualizing binding pocket resolvation across ABFE simulations is straightforward with FEPA. We follow a workflow similar to Tutorial 1, but instead of `SelfDistanceFeaturizer`, we use `BPWaterFeaturizer` to quantify water occupancy. + +### Load Config and Prepare Paths + +```python +import logging, os +from fepa.utils.file_utils import load_config +from fepa.utils.path_utils import load_abfe_paths_for_compound + +logging.basicConfig(level=logging.INFO, format="%(asctime)s - %(levelname)s - %(message)s") + +config = load_config("../../config/config.json") +analysis_output_dir = "wdir" +cmp = config["compounds"][0] + +logging.info("Analyzing compound %s ...", cmp) +cmp_output_dir = os.path.join(analysis_output_dir, cmp) +os.makedirs(cmp_output_dir, exist_ok=True) + +logging.info("Loading paths for compound %s...", cmp) +path_dict = load_abfe_paths_for_compound( + config, + cmp, + van_list=[1, 2, 3], + leg_window_list=[f"coul.{i:02d}" for i in range(0, 11)] + + [f"rest.{i:02d}" for i in range(0, 12)] + + [f"vdw.{i:02d}" for i in range(0, 21)], + bp_selection_string="name CA and resid 57 58 61 64 83 84 87 88 91 92 173 177 218 221 235 238 239 242 243 246", + apo=False, +) +``` + +### Load trajectories and featurize waters + +```python +from fepa.core.ensemble_handler import EnsembleHandler +from fepa.core.featurizers import BPWaterFeaturizer + +logging.info("Loading trajectories for compound %s ...", cmp) +ensemble_handler = EnsembleHandler(path_dict) +ensemble_handler.make_universes() + +logging.info("Featurizing binding pocket waters ...") +bp_water_featurizer = BPWaterFeaturizer(ensemble_handler=ensemble_handler) +bp_water_featurizer.featurize(radius=10) # Count waters within 10 Å of pocket COM + +logging.info("Saving features for compound %s ...", cmp) +bp_water_featurizer.save_features(cmp_output_dir, overwrite=True) +``` + +### Plot Water Occupancy Over Time or Windows + +```python +import pandas as pd +import matplotlib.pyplot as plt +import seaborn as sns + +features_df = pd.read_csv(os.path.join(cmp_output_dir, "WaterOccupancy_features.csv")) + +for van in [1, 2, 3]: + van_features_df = features_df[features_df["ensemble"].str.contains(f"van_{van}")] + plt.figure(figsize=(12, 6)) + sns.lineplot(data=van_features_df, x="Time (ps)", y="occupancy", hue="ensemble") + plt.title(f"Water Occupancy for {cmp}") + plt.xlabel("Time (ps)") + plt.ylabel("Number of Waters") + plt.xlim(0, 20000) + plt.legend(title="Ensemble", bbox_to_anchor=(1.05, 1), loc="upper left", ncol=2) + plt.tight_layout() + plt.savefig(os.path.join(cmp_output_dir, f"{cmp}_water_occupancy_van{van}.png")) + +``` + +```python +# Extract vanilla replicate and window ID +features_df["van"] = features_df["ensemble"].str.extract(r"van_(\d)") +features_df["id"] = features_df["ensemble"].str.replace(r"_van_\d+", "", regex=True) + +# Compute average occupancy +avg_df = features_df.groupby(["id", "van"], as_index=False)["occupancy"].mean() +avg_df.to_csv(os.path.join(cmp_output_dir, "avg_water_occupancy.csv"), index=False) + +# Plot average occupancy +plt.figure(figsize=(12, 8)) +sns.lineplot(data=avg_df, x="id", y="occupancy", hue="van", palette="tab10") +plt.title(f"Average Water Occupancy Across Windows for {cmp}", fontsize=16, fontweight="bold") +plt.xlabel("Window ID", fontsize=14) +plt.ylabel("Average Number of Waters", fontsize=14) +plt.legend(title="Vanilla Repeat", title_fontsize=12, fontsize=10, loc="upper right") +plt.xticks(rotation=45, fontsize=10) +plt.yticks(fontsize=10) +plt.tight_layout() +plt.savefig(os.path.join(cmp_output_dir, f"{cmp}_water_occupancy_across_windows.png"), dpi=300) +plt.close() + +``` \ No newline at end of file diff --git a/site/how-to-guides/guide-1/index.html b/site/how-to-guides/guide-1/index.html index 4fa7b6a..76b5afc 100644 --- a/site/how-to-guides/guide-1/index.html +++ b/site/how-to-guides/guide-1/index.html @@ -27,7 +27,7 @@ - How-to-Guide - FEPA Documentation + Binding Pocket Resolvation Analysis - FEPA Documentation @@ -111,6 +111,11 @@
@@ -143,7 +148,7 @@
- How-to-Guide + Binding Pocket Resolvation Analysis
@@ -519,6 +524,26 @@ + + @@ -538,6 +563,57 @@ + + + + + @@ -599,6 +675,47 @@ + + +
@@ -660,9 +777,95 @@ -

How-to-Guide

- -

Put the example analysis in the how to guide section here

+

Binding Pocket Resolvation Analysis

+ +
+

icon: lucide/compass +title: Binding Pocket Resolvation Analysis

+
+

Visualizing binding pocket resolvation across ABFE simulations is straightforward with FEPA. We follow a workflow similar to Tutorial 1, but instead of SelfDistanceFeaturizer, we use BPWaterFeaturizer to quantify water occupancy.

+

Load Config and Prepare Paths

+
import logging, os
+from fepa.utils.file_utils import load_config
+from fepa.utils.path_utils import load_abfe_paths_for_compound
+
+logging.basicConfig(level=logging.INFO, format="%(asctime)s - %(levelname)s - %(message)s")
+
+config = load_config("../../config/config.json")
+analysis_output_dir = "wdir"
+cmp = config["compounds"][0]
+
+logging.info("Analyzing compound %s ...", cmp)
+cmp_output_dir = os.path.join(analysis_output_dir, cmp)
+os.makedirs(cmp_output_dir, exist_ok=True)
+
+logging.info("Loading paths for compound %s...", cmp)
+path_dict = load_abfe_paths_for_compound(
+    config,
+    cmp,
+    van_list=[1, 2, 3],
+    leg_window_list=[f"coul.{i:02d}" for i in range(0, 11)]
+                     + [f"rest.{i:02d}" for i in range(0, 12)]
+                     + [f"vdw.{i:02d}" for i in range(0, 21)],
+    bp_selection_string="name CA and resid 57 58 61 64 83 84 87 88 91 92 173 177 218 221 235 238 239 242 243 246",
+    apo=False,
+)
+
+

Load trajectories and featurize waters

+
from fepa.core.ensemble_handler import EnsembleHandler
+from fepa.core.featurizers import BPWaterFeaturizer
+
+logging.info("Loading trajectories for compound %s ...", cmp)
+ensemble_handler = EnsembleHandler(path_dict)
+ensemble_handler.make_universes()
+
+logging.info("Featurizing binding pocket waters ...")
+bp_water_featurizer = BPWaterFeaturizer(ensemble_handler=ensemble_handler)
+bp_water_featurizer.featurize(radius=10)  # Count waters within 10 Ã… of pocket COM
+
+logging.info("Saving features for compound %s ...", cmp)
+bp_water_featurizer.save_features(cmp_output_dir, overwrite=True)
+
+

Plot Water Occupancy Over Time or Windows

+
import pandas as pd
+import matplotlib.pyplot as plt
+import seaborn as sns
+
+features_df = pd.read_csv(os.path.join(cmp_output_dir, "WaterOccupancy_features.csv"))
+
+for van in [1, 2, 3]:
+    van_features_df = features_df[features_df["ensemble"].str.contains(f"van_{van}")]
+    plt.figure(figsize=(12, 6))
+    sns.lineplot(data=van_features_df, x="Time (ps)", y="occupancy", hue="ensemble")
+    plt.title(f"Water Occupancy for {cmp}")
+    plt.xlabel("Time (ps)")
+    plt.ylabel("Number of Waters")
+    plt.xlim(0, 20000)
+    plt.legend(title="Ensemble", bbox_to_anchor=(1.05, 1), loc="upper left", ncol=2)
+    plt.tight_layout()
+    plt.savefig(os.path.join(cmp_output_dir, f"{cmp}_water_occupancy_van{van}.png"))
+
+
# Extract vanilla replicate and window ID
+features_df["van"] = features_df["ensemble"].str.extract(r"van_(\d)")
+features_df["id"] = features_df["ensemble"].str.replace(r"_van_\d+", "", regex=True)
+
+# Compute average occupancy
+avg_df = features_df.groupby(["id", "van"], as_index=False)["occupancy"].mean()
+avg_df.to_csv(os.path.join(cmp_output_dir, "avg_water_occupancy.csv"), index=False)
+
+# Plot average occupancy
+plt.figure(figsize=(12, 8))
+sns.lineplot(data=avg_df, x="id", y="occupancy", hue="van", palette="tab10")
+plt.title(f"Average Water Occupancy Across Windows for {cmp}", fontsize=16, fontweight="bold")
+plt.xlabel("Window ID", fontsize=14)
+plt.ylabel("Average Number of Waters", fontsize=14)
+plt.legend(title="Vanilla Repeat", title_fontsize=12, fontsize=10, loc="upper right")
+plt.xticks(rotation=45, fontsize=10)
+plt.yticks(fontsize=10)
+plt.tight_layout()
+plt.savefig(os.path.join(cmp_output_dir, f"{cmp}_water_occupancy_across_windows.png"), dpi=300)
+plt.close()
+
diff --git a/site/search.json b/site/search.json index 930943a..b20bd14 100644 --- a/site/search.json +++ b/site/search.json @@ -1 +1 @@ -{"config":{"separator":"[\\s\\-_,:!=\\[\\]()\\\\\"`/]+|\\.(?!\\d)"},"items":[{"location":"","level":1,"title":"Getting Started with FEPA","text":"

FEPA (Free Energy Perturbation Analysis) is a Python package for analyzing molecular dynamics (MD) trajectories from FEP simulations, particularly ABFEs. FEPA allows you to visualize conformational changes and set up simulations to correct free energy estimates.

This guide covers installation, basic usage, and key workflows.

","path":["Getting Started"],"tags":[]},{"location":"#installation","level":2,"title":"Installation","text":"

FEPA is must be installed from GitHub:

git clone https://github.com/Nithishwer/FEPA.git\ncd FEPA\npip install -e .\n
","path":["Getting Started"],"tags":[]},{"location":"explanation/","level":1,"title":"Explanation section","text":"

Lorem Ipsum

","path":["Explanation"],"tags":[]},{"location":"reference/","level":1,"title":"Reference","text":"

This part of the project documentation focuses on an information-oriented approach. Use it as a reference for the technical implementation of the calculator project code.

::: fepa.core.analyzers ::: fepa.core.dim_reducers ::: fepa.core.ensemble_handler ::: fepa.core.featurizers ::: fepa.core.visualizers

","path":["Reference"],"tags":[]},{"location":"how-to-guides/guide-1/","level":1,"title":"How-to-Guide","text":"

Put the example analysis in the how to guide section here

","path":["How-to-guides","How-to-Guide"],"tags":[]},{"location":"how-to-guides/guide-2/","level":1,"title":"How-to-Guide","text":"

Put the example analysis in the how to guide section here

","path":["How-to-guides","How-to-Guide"],"tags":[]},{"location":"tutorials/tutorial-1/","level":1,"title":"Featurizing ABFEs","text":"

This tutorial is designed to the user to introduce fepa by analyzing MD data from a set of simulations of the ligand 42922 (from the Deflorian et al set 1 datset) bound to the Orexin 2 Receptor. Specifically, we will compare the holo, apo and ABFE trajectories. Through this comparison, we will demonstrate that the protein’s binding pocket adopts distinct conformations between the final lambda windows (where the ligand is fully annihilated) of the ABFE simulation and the long apo simulation of the receptor.

This tutorial is followed by Tutorial 2 where we will use FEPA to set up REUS simulations to estimate the free energy of this conformational change from the holo like to apo like states. This tutorial uses MDAnalysis, MDtraj and the PENSA package for analysis, GROMACS for simulation and Plumed for enhanced sampling. It is assumed that the user is familiar with setting up and analyzing MD simulations run with GROMACS and Plumed.

","path":["Tutorials","Featurizing ABFEs"],"tags":[]},{"location":"tutorials/tutorial-1/#loading-the-config-file","level":2,"title":"Loading the config file","text":"

ABFE simulations produce multiple MD trajectories over different lamba windows. In our case, we have 44 lambda windows (11 Coulomb, 12 Restrains and 22 Van-der-Walls) plus the holo and the apo simulations that we have to analyse. To make things easier when parsing paths to the topology and coordinate files of these trajectories, we use a config file. The config file is json-formatted and contains all the information necessary to read the simulation trajectories. Here is a sample config file that I use for this tutorial:

{\n    \"base_path\": \"deflorian_set_1_j13_v1\",\n    \"abfe_window_path_template\": \"deflorian_set_1_j13_v1/OX2_{CMP_NAME}/abfe_van{REP_NO}_hrex_r{ABFE_REP_NO}/complex/{LEG_WINDOW}/{STAGE}\",\n    \"vanilla_path_template\": \"deflorian_set_1_j13_v1/OX2_{CMP_NAME}/vanilla_rep_{REP_NO}\",\n    \"apo_path_template\": \"deflorian_set_1_j13_v1/apo_OX2_r{REP_NO}\",\n    \"compounds\": [\n        \"42922\",\n    ],\n    \"pocket_residues_string\": \"12  54  57  58  59  60  61  62  63  64  65  70  71  78  81  82  83  85  86  89 138 142 160 161 162 163 175 178 179 182 183 232 235 236 239 240 242 243 261 265 268 269\"\n}\n

The pocket_residues_string variable is a string that stores the residue ids of all the proteins residue that have any atom within 6 A of the ligand. The JSON file is then loaded into a dictionary with the function load_config from fepa.utils.file_utils

","path":["Tutorials","Featurizing ABFEs"],"tags":[]},{"location":"tutorials/tutorial-1/#loading-md-trajectories","level":2,"title":"Loading MD trajectories","text":"

Now that we have the template paths to all the simulations: Apo EQ, Holo EQ and ABFE, we can use the function load_abfe_paths_for_compound from fepa.utils.path_utils to generate a path_dict dictionary that contains all the MD file paths for a single compound as follows:

cmp = config['compounds'][0]\npath_dict = load_abfe_paths_for_compound(\n            config,\n            cmp,\n            van_list=[3], # Loading only vanilla simulation 3\n            leg_window_list=[f\"coul.{i:02d}\" for i in range(0, 11,2)] + [f'vdw.{i:02d}' for i in range(0, 21,2)] + [f'rest.{i:02d}' for i in range(0, 12,2)], # Loading only every other lambda window\n            bp_selection_string=\"name CA and resid \" + config[\"pocket_residues_string\"],\n            apo=True,\n        )\n

The load_abfe_paths_for_compound function has various arguments that allow the user control over what simulation paths are loaded. For more information, please refer to the API.

Now that we have the path_dict, it is time for us to load the trajectories themselves. We will be using the EnsembleHandler class from fepa.core.ensemble_handler to do this. EnsembleHandler is a neat way of storing and manipulating the trajectories from multiple ensembles with some built in functions for sanity checks. Internally EnsembleHandler stores the trajectories as dictionary of MDA universes.

TODO: Need to remove ensemblehandler and just use universe dict

# Load trajectories\nensemble_handler = EnsembleHandler(path_dict)\n# Make universes\nensemble_handler.make_universes()\nlogging.info(\"Making universes for compound %s...\", cmp)\n# Check for BP residue consistency across all the trajectories\nlogging.info(\"Checking residue consistency for compound %s...\", cmp)\n
","path":["Tutorials","Featurizing ABFEs"],"tags":[]},{"location":"tutorials/tutorial-1/#featurizing-md-trajectories","level":2,"title":"Featurizing MD trajectories","text":"

We will now be featurizing these trajectories by their pairwise C-alpha distances in the binding pocket using the function SelfDistanceFeaturizer from fepa.core.featurizers. SelfDistanceFeaturizer computes and stored all possible pairs of distances between the C-alpha atoms of the binding pocket residues. The class also has save_features and load_features functions that help save the features as a csv file to make sure the time consuming featurization step need not be repeated every run.

TODO: Make sure the binding pocket selection is consistent

# Make a folder for the analysis output\ncmp_run_dir = f'analysis/{cmp}/'\nif !(os.path.exists(cmp_run_dir)):\n    os.mkdir(cmp_run_dir)\n\n# Featurize and save features\nfeaturizer = SelfDistanceFeaturizer(ensemble_handler)\nfeaturizer.featurize()\nfeaturizer.save_features(input_dir=cmp_existing_run_dir)\n
","path":["Tutorials","Featurizing ABFEs"],"tags":[]},{"location":"tutorials/tutorial-1/#visualizing-ensembles","level":2,"title":"Visualizing ensembles","text":"

Saving the features in a csv format gives us the flexibility to analyse it as required. For the purpose of this tutorial, we will be looking at how our features capture the difference between the apo, the holo and the abfe ensembles. To do this we reduce the dimensions of the features data using the PCADimReducer class from fepa.core.dim_reducers. FEPA also supports other dimensionality reduction techniqeus like UMAP and tSNE. In fact UMAP is better able to resolve the differences in binding pocket configurations between different ensembles. We will be doing PCA here as it also doubles up as a nice CV to bias when performing umbrella sampling later.

# Dimensionality Reduction\nlogging.info(\"Performing dimensionality reduction for compound %s...\", cmp)\ndimreducer = PCADimReducer(featurizer.get_feature_df(), n_components=8)\ndimreducer.reduce_dimensions()\ndimreducer.calculate_projections()\ndimreducer.save_projection_df(\n    save_path=os.path.join(cmp_output_dir, \"pca_projection_df.csv\")\n)\n

First we plot the eigen values of all the PCs to understand what percentage of variance is captured by the first few PCs:

logging.info(\"Plotting PCA eigenvalues for compound %s...\", cmp)\nplot_eigenvalues(\n    pca_object=dimreducer.get_pca(),\n    n_components=8,\n    save_path=os.path.join(cmp_output_dir, \"eigenvalues.png\"),\n)\n

Figure 1: PCA eigenvalues for the binding-pocket Cα self-distance features.

Here, most of the variance is captured by PC1. However, this may not always be the case. Therefore, for this tutorial, we will use PC1 and PC2 together, as they collectively capture the majority of the variance. To improve visualization, we define a new column in projection_df called simtype, which groups the ensembles into broader categories: abfe_coul, abfe_vdw, abfe_rest, holo_equil, apo, and nvt.

# Further labelling the ensembles\nprojection_df = dimreducer.get_pca_projection_df()\n# Add another column called simtype based on ensemble\ndef get_simtype(ensemble_name: str) -> str:\n    if \"van\" in ensemble_name:\n        if 'nvt' in ensemble_name:\n            return \"nvt\"\n        if \"coul\" in ensemble_name:\n            return \"abfe_coul\"\n        elif \"vdw\" in ensemble_name:\n            return \"abfe_vdw\"\n        elif \"rest\" in ensemble_name:\n            return \"abfe_rest\"\n        else:\n            return \"holo_equil\"\n    elif \"apo\" in ensemble_name:\n        return \"apo\"\n    else:\n        return \"other\"\nprojection_df[\"simtype\"] = projection_df[\"ensemble\"].apply(get_simtype)\n

We can visualize the dimensionality-reduced data using the DimRedVisualizer class. In the example below, we plot the first two principal components, colored by simulation and time. We also highlight NVT, as it represents the initial crystal structure after energy minimization.

logging.info(\"Visualizing compound %s...\", cmp)\ndimred_visualizer = DimRedVisualizer(projection_df=projection_df, data_name=\"PCA\")\ndimred_visualizer.plot_dimred_sims(\n    save_path=os.path.join(cmp_run_dir, \"pca_components_ensemble_noapo.png\"),\n    highlights=[f\"{cmp}_nvt\"],\n)\ndimred_visualizer.plot_dimred_time(\n    save_path=os.path.join(cmp_run_dir, \"pca_components_time_noapo.png\")\n)\ndimred_visualizer.plot_dimred_sims(\n    column=\"simtype\",\n    save_path=os.path.join(cmp_run_dir, \"pca_components_simtype.png\"),\n)\n

Figure 2: Simulation frames in PC-space (PC1 vs PC2), colored by simtype

","path":["Tutorials","Featurizing ABFEs"],"tags":[]},{"location":"tutorials/tutorial-1/#clustering-ensembles","level":2,"title":"Clustering ensembles","text":"

In Figure 2, frames from the ABFE simulations closely resemble those from the holo equilibrium simulation, forming two distinct clusters: apo-like and holo-like. We will be clustering them using cluster_pca. After clustering, we can estimtate the data point that is closest to the center with the function make_ensemble_center_df which returns a DataFrame containing the details of the frame closest to the centroid in PC space. 

# Cluster the projection df\npca_projection_df_clustered = cluster_pca(\n    projection_df, n_clusters=3, n_components=3\n)\n# Ensemble center df\nensemble_center_df = make_ensemble_center_df(\n    pca_projection_df_clustered, key=\"cluster\"\n)\n

We then visualize the clustered data with ensemble centers as follows:

# Visualization of clustered data\ndimred_visualizer_clustered = DimRedVisualizer(\n    projection_df=pca_projection_df_clustered, data_name=\"PCA\"\n)\ndimred_visualizer_clustered.plot_dimred_cluster(\n    save_path=os.path.join(\n        cmp_run_dir, \"subset_pca_components_clusters_w_center.png\"\n    ),\n    centroid_df=ensemble_center_df,\n    cluster_column=\"cluster\",\n)\n

Figure 3: Clustered PCA data with ensemble centers marked by an X.

","path":["Tutorials","Featurizing ABFEs"],"tags":[]},{"location":"tutorials/tutorial-1/#visualizing-differences-across-clustered-ensembles","level":2,"title":"Visualizing differences across clustered ensembles","text":"

With two distinct clusters identified, we can compare them across different features by computing the Jensen–Shannon (JS) entropy for each feature. The code below plots feature-wise histograms for cluster 0 and cluster 1, annotated with their JS entropy values:

# Copy the cluster labels from pca_projection_df_clustered to feature_df\nfeature_df_w_clusters = featurizer.get_feature_df().copy()\nfeature_df_w_clusters[\"cluster\"] = pca_projection_df_clustered[\"cluster\"].values\n\n\n# Compute feature-level histograms and JS entropy between clusters 0 and 1\nhistograms = compute_histograms(\n    feature_df_w_clusters,\n    \"cluster\",\n    0,\n    1,\n    num_bins=50,\n    feature_column_keyword=\"DIST\",\n)\n\nrel_ent_dict = compute_relative_entropy(\n    feature_df_w_clusters,\n    ensemble1=0,\n    ensemble2=1,\n    num_bins=50,\n    key=\"cluster\",\n    feature_column_keyword=\"DIST\",\n)\n\nplot_jsd_histograms(\n    histograms,\n    rel_ent_dict,\n    top_n=16,\n    save_path=os.path.join(cmp_run_dir, \"jsd_cluster0_vs_cluster1.png\"),\n)\n

Figure 4: Histograms of the top 16 features showing the highest Jensen–Shannon divergence between clusters 0 and 1.

We then save the ensemble center frames as GRO files to compare structures in molecular visualization tools such as PyMOL:

# Save ensemble center frames as GRO files\nfor _, row in ensemble_center_df.iterrows():\n    center_ensemble = row[\"ensemble\"]\n    center_timestep = row[\"timestep\"]\n    print(f\"Ensemble: {center_ensemble}, Timestep: {center_timestep}\")\n\n    # Load trajectories\n    ensemble_handler.make_universes()\n\n    # Export the frame corresponding to the ensemble center\n    ensemble_handler.dump_frames(\n        ensemble=center_ensemble,\n        timestep=center_timestep,\n        save_path=f\"subset_cluster_{int(row['cluster'])}_center.gro\",\n    )\n
","path":["Tutorials","Featurizing ABFEs"],"tags":[]},{"location":"tutorials/tutorial-2/","level":1,"title":"Setting up REUS","text":"

In Tutorial 2, we build on the concepts learned in Tutorial 1, where the apo and holo states were shown to form distinct conformational clusters. Using the data provided in the [template], we will set up Replica Exchange Umbrella Sampling (REUS) simulations to estimate the free energy of the conformational transition from holo to apo, using the principal component (PC) 1 as the collective variable (CV). We will achieve this by morphing the protein from one state to another, equilibrating the system, and then running REUS to sample the transition.

","path":["Tutorials","Setting up REUS"],"tags":[]},{"location":"tutorials/tutorial-2/#preparing-the-plumeddat-file","level":2,"title":"Preparing the plumed.dat file","text":"

The first step in performing umbrella sampling with GROMACS and PLUMED is to create a plumed.dat file that defines the collective variable (CV), and the positions of the umbrella sampling restraints. To do this, we load the feature CSV files, initialize a DimReducer object, and extract the corresponding PCA object to define the CV based on pairwise Cα distances.

# Load Features df and reduce dimensions\nfeature_df = pd.read_csv(os.path.join(\"md_data\", \"top_features_apo_vdw.20.csv\"))\ndimreducer = PCADimReducer(feature_df, n_components=8)\ndimreducer.reduce_dimensions()\ndimreducer.calculate_projections()\npca_projection_df = dimreducer.get_pca_projection_df()\n# Use PCA from dimreducer to write plumed file\ntop_features = feature_df.filter(regex=\"DIST\", axis=1).columns\nwrite_plumed_file(\n    sdf_names=top_features,\n    top_features_pca=dimreducer.get_pca(),\n    save_path=\"plumed.dat\",\n    molinfo_structure=\"../reference.pdb\",  # fix molinfo here\n)\n

The PLUMED input file should look like this:

MOLINFO STRUCTURE=../reference.pdb\nd1: DISTANCE ATOMS=@CA-236,@CA-243\nd2: DISTANCE ATOMS=@CA-175,@CA-243\nd3: DISTANCE ATOMS=@CA-138,@CA-243\n.\n.\n.\nd200: DISTANCE ATOMS=@CA-61,@CA-70\n# Create the dot product\ndot: COMBINE ARG=d1,d2,d3...d200 COEFFICIENTS=0.095,0.152,0.133...0.057,0.054,0.032 PERIODIC=NO\nCV: MATHEVAL ARG=dot FUNC=10*x PERIODIC=NO\nPRINT ARG=CV FILE=COLVAR STRIDE=1\n# Put position of restraints here for each window\nrestraint: RESTRAINT ARG=CV AT=@replicas:$RESTRAINT_ARRAY KAPPA=$KAPPA\nPRINT ARG=restraint.* FILE=restr\n
Note that $RESTRAINT_ARRAY is a placeholder for the harmonic restraint positions in CV space. When defining these positions, ensure they align with the direction of the structural morph. By convention, we transition from vdw.20 to apo structures. To assign restraints correctly, identify the ensemble with the lowest PC value and order the restraints from min to max (or vice versa) accordingly.

# Pair of ensembles to compare\npair = (\"vdw.20\", \"apo\")\n# Get mean PC1 for the two ensembles\nmean_pc1_ensemble1 = pca_projection_df[pca_projection_df[\"state\"] == pair[0]][\n    \"PC1\"\n].mean()\nmean_pc1_ensemble2 = pca_projection_df[pca_projection_df[\"state\"] == pair[1]][\n    \"PC1\"\n].mean()\nlogging.info(\n    f\"Mean PC1 for {pair[0]}: {mean_pc1_ensemble1}, Mean PC1 for {pair[1]}: {mean_pc1_ensemble2}\"\n)\n# Get the restraint array based on the two ensembles\nPC1_min = pca_projection_df[\"PC1\"].min()\nPC1_max = pca_projection_df[\"PC1\"].max()\n# If min PC1 is from ensemble 1, then the restraint array should be from PC1_min to PC1_max\nif mean_pc1_ensemble1 < mean_pc1_ensemble2:\n    restraint_array = np.linspace(PC1_min, PC1_max, 24)\nelif mean_pc1_ensemble1 > mean_pc1_ensemble2:\n    restraint_array = np.linspace(PC1_max, PC1_min, 24)\n

Once the restraint array is prepared, we write it to the PLUMED file using the write_plumed_restraints function:

# Write the restraint array to the plumed file\nwrite_plumed_restraints(\n    plumed_file=\"plumed.dat\",\n    restraint_centers=restraint_array,\n    kappa=5,\n)\n

The restrain line on plumed.dat file must now be an array of distance restraints that looks like this:

restraint: RESTRAINT ARG=CV AT=@replicas:116.856,119.047,121.239...167.262 KAPPA=5\n
","path":["Tutorials","Setting up REUS"],"tags":[]},{"location":"tutorials/tutorial-2/#morphing-with-memento","level":2,"title":"Morphing with memento","text":"

For this tutorial, the GRO files representing the vdw.20 and apo states are provided in the data folder, prepared in the same way as in Tutorial 1. The plumed.dat file has also been prepared. The next step for REUS is to generate intermediate protein conformations using Memento (JCTC, 2023) . FEPA provides a class for easy access to Memento.

Using FEPA’s memento_workflow class, we can perform the protein morphing. A Memento directory is needed to store morphs for each entry in initial_target_gro_dict, and an apo_template path must be provided, containing the topology and MDP files required for equilibrium simulations.

# Declaring variables:\nmemento_dir = \"/biggin/b211/reub0138/Projects/orexin/deflorian_set_1_j13_v1_memento\"\n# Defining pairs\npair = (\"all_vdw20\", \"apo_3\")\nall_vdw20_gro_file = \"/biggin/b211/reub0138/Projects/orexin/deflorian_set_1_j13_v1/analysis/a12_dimreduce_vdw.20/subset_cluster_1_center.gro\"\napo_3_gro_file = \"/biggin/b211/reub0138/Projects/orexin/deflorian_set_1_j13_v1/analysis/a12_dimreduce_vdw.20/subset_cluster_2_center.gro\"\ntemplate_path = \"/biggin/b211/reub0138/Projects/orexin/deflorian_set_1_j13_v1_memento/apo_template\"  # Make sure no water and no ions in the template topology file\n# Setup Memento folders\nmemento_flow = memento_workflow(\n    memento_dir=memento_dir,\n    initial_gro=all_vdw20_gro_file,\n    target_gro=apo_3_gro_file,\n    initial_name=pair[0],\n    target_name=pair[1],\n    template_path=template_path,\n    run_name=\"memento_run_v1\",\n    n_residues=296,\n)\n

We can use the workflow functions to perform each step. First, prepare_memento sets up the files required by Memento. Then, run_memento executes the morphing. run_memento requires the protonation states of all histidines in the format expected by GROMACS pdb2gmx (0 for HID, 1 for HIE, 2 for HIP, 3 for HIS1). Additionally, the indices of CYX residues must be provided, as Memento cannot process them automatically.

# Preparing memento input\nmemento_flow.prepare_memento()\n# Running memento\nmemento_flow.run_memento(\n    template_path=template_path,\n    last_run=\"memento_run_v1\",\n    protonation_states=[1, 1, 1, 2, 1, 2],\n    cyx_residue_indices = []\n)\n
","path":["Tutorials","Setting up REUS"],"tags":[]},{"location":"tutorials/tutorial-2/#running-equilibration","level":2,"title":"Running Equilibration","text":"

This step may take some time. Once complete, we can set up equilibration simulations with Memento by providing a job script path, which can be modified to suit the specific HPC system.

# Running analysis\nmemento_flow.prepare_equil_simulations(\n    job_script_template=\"/biggin/b211/reub0138/Projects/orexin/lenselink_a2a_memento_v1/job_vanilla_ranv_equil_arr_template.sh\"\n)\n

After running prepare_equil_simulations, individual boxes for each morph are created. We then simulate each box to relax the side chains by submitting the job to our HPC. Once complete, we can proceed to the REUS simulations.

","path":["Tutorials","Setting up REUS"],"tags":[]},{"location":"tutorials/tutorial-2/#preparing-reus","level":2,"title":"Preparing REUS","text":"

Once the equilibration simulations are complete, use FEPA’s reus_umbrella_sampling_workflow class to set up umbrella sampling in the same working directory. The number of windows can be adjusted via n_windows (24 is typically sufficient). If the residue IDs in the plumed.dat file does not match those in the MD GRO files, FEPA requires the correct offset to be set via the plumed_resid_offset parameter.

import logging\nfrom pathlib import Path\nfrom fepa.flows.reus_flows import (\n    reus_umbrella_sampling_workflow,\n)\n\nsim_path = \"wdir/all_vdw20_apo_3/memento_run_v1/wdir/boxes/sim0\"\nwdir_path = Path(sim_path).parents[1]\nplumed_path = \"plumed.dat\"\nsimname = \"all_vdw20_apo_3\"\ninitial_gro = \"md_data/cluster_1_center.gro\"\nfinal_gro = \"md_data/cluster_2_center.gro\"\nsubmission_script_template_arr = \"md_data/job_reus_template.sh\"\nprint(f\"wdir_path: {wdir_path}, plumed_path: {plumed_path}\")\numbrella_sampler = reus_umbrella_sampling_workflow(\n    wdir_path=wdir_path,\n    plumed_path=plumed_path,\n    submission_script_template_arr=submission_script_template_arr,\n    start=simname.split(\"_\")[0] + \"_\" + simname.split(\"_\")[1],\n    end=simname.split(\"_\")[2] + \"_\" + simname.split(\"_\")[3],\n    reus_folder_name=\"reus_v1\",\n    n_windows=24,\n    plumed_resid_offset=0,\n    initial_gro=initial_gro,\n    target_gro=final_gro,\n)\numbrella_sampler.setup_simulations(exist_ok=True)\n
","path":["Tutorials","Setting up REUS"],"tags":[]},{"location":"tutorials/tutorial-2/#analyzing-reus","level":2,"title":"Analyzing REUS","text":"

Once your REUS system is set up, you can analyze the results using the workflow:

umbrella_sampler.prepare_wham()\numbrella_sampler.run_wham()\numbrella_sampler.analyse_us_hist(range=(90, 180), colvar_prefix=\"COLVAR\")\numbrella_sampler.get_initial_final_CVs()\numbrella_sampler.plot_free_energies(\n    units=\"kcal\",\n)\n

This prepares and runs WHAM on the histograms and generates the free energy curves.

Figure 1: Free energy curves from the tutorial simulations. The CV values of apo and holo structures are marked with gray dotted lines. Different lines represent curves computed using varying proportions of data to assess convergence.

","path":["Tutorials","Setting up REUS"],"tags":[]}]} \ No newline at end of file +{"config":{"separator":"[\\s\\-_,:!=\\[\\]()\\\\\"`/]+|\\.(?!\\d)"},"items":[{"location":"","level":1,"title":"Getting Started with FEPA","text":"

FEPA (Free Energy Perturbation Analysis) is a Python package for analyzing molecular dynamics (MD) trajectories from FEP simulations, particularly ABFEs. FEPA allows you to visualize conformational changes and set up simulations to correct free energy estimates.

This guide covers installation, basic usage, and key workflows.

","path":["Getting Started"],"tags":[]},{"location":"#installation","level":2,"title":"Installation","text":"

FEPA is must be installed from GitHub:

git clone https://github.com/Nithishwer/FEPA.git\ncd FEPA\npip install -e .\n
","path":["Getting Started"],"tags":[]},{"location":"explanation/","level":1,"title":"Explanation section","text":"

Lorem Ipsum

","path":["Explanation"],"tags":[]},{"location":"reference/","level":1,"title":"Reference","text":"

This part of the project documentation focuses on an information-oriented approach. Use it as a reference for the technical implementation of the calculator project code.

::: fepa.core.analyzers ::: fepa.core.dim_reducers ::: fepa.core.ensemble_handler ::: fepa.core.featurizers ::: fepa.core.visualizers

","path":["Reference"],"tags":[]},{"location":"how-to-guides/guide-1/","level":1,"title":"Binding Pocket Resolvation Analysis","text":"

icon: lucide/compass title: Binding Pocket Resolvation Analysis

Visualizing binding pocket resolvation across ABFE simulations is straightforward with FEPA. We follow a workflow similar to Tutorial 1, but instead of SelfDistanceFeaturizer, we use BPWaterFeaturizer to quantify water occupancy.

","path":["How-to-guides","Binding Pocket Resolvation Analysis"],"tags":[]},{"location":"how-to-guides/guide-1/#load-config-and-prepare-paths","level":3,"title":"Load Config and Prepare Paths","text":"
import logging, os\nfrom fepa.utils.file_utils import load_config\nfrom fepa.utils.path_utils import load_abfe_paths_for_compound\n\nlogging.basicConfig(level=logging.INFO, format=\"%(asctime)s - %(levelname)s - %(message)s\")\n\nconfig = load_config(\"../../config/config.json\")\nanalysis_output_dir = \"wdir\"\ncmp = config[\"compounds\"][0]\n\nlogging.info(\"Analyzing compound %s ...\", cmp)\ncmp_output_dir = os.path.join(analysis_output_dir, cmp)\nos.makedirs(cmp_output_dir, exist_ok=True)\n\nlogging.info(\"Loading paths for compound %s...\", cmp)\npath_dict = load_abfe_paths_for_compound(\n    config,\n    cmp,\n    van_list=[1, 2, 3],\n    leg_window_list=[f\"coul.{i:02d}\" for i in range(0, 11)]\n                     + [f\"rest.{i:02d}\" for i in range(0, 12)]\n                     + [f\"vdw.{i:02d}\" for i in range(0, 21)],\n    bp_selection_string=\"name CA and resid 57 58 61 64 83 84 87 88 91 92 173 177 218 221 235 238 239 242 243 246\",\n    apo=False,\n)\n
","path":["How-to-guides","Binding Pocket Resolvation Analysis"],"tags":[]},{"location":"how-to-guides/guide-1/#load-trajectories-and-featurize-waters","level":3,"title":"Load trajectories and featurize waters","text":"
from fepa.core.ensemble_handler import EnsembleHandler\nfrom fepa.core.featurizers import BPWaterFeaturizer\n\nlogging.info(\"Loading trajectories for compound %s ...\", cmp)\nensemble_handler = EnsembleHandler(path_dict)\nensemble_handler.make_universes()\n\nlogging.info(\"Featurizing binding pocket waters ...\")\nbp_water_featurizer = BPWaterFeaturizer(ensemble_handler=ensemble_handler)\nbp_water_featurizer.featurize(radius=10)  # Count waters within 10 Ã… of pocket COM\n\nlogging.info(\"Saving features for compound %s ...\", cmp)\nbp_water_featurizer.save_features(cmp_output_dir, overwrite=True)\n
","path":["How-to-guides","Binding Pocket Resolvation Analysis"],"tags":[]},{"location":"how-to-guides/guide-1/#plot-water-occupancy-over-time-or-windows","level":3,"title":"Plot Water Occupancy Over Time or Windows","text":"
import pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nfeatures_df = pd.read_csv(os.path.join(cmp_output_dir, \"WaterOccupancy_features.csv\"))\n\nfor van in [1, 2, 3]:\n    van_features_df = features_df[features_df[\"ensemble\"].str.contains(f\"van_{van}\")]\n    plt.figure(figsize=(12, 6))\n    sns.lineplot(data=van_features_df, x=\"Time (ps)\", y=\"occupancy\", hue=\"ensemble\")\n    plt.title(f\"Water Occupancy for {cmp}\")\n    plt.xlabel(\"Time (ps)\")\n    plt.ylabel(\"Number of Waters\")\n    plt.xlim(0, 20000)\n    plt.legend(title=\"Ensemble\", bbox_to_anchor=(1.05, 1), loc=\"upper left\", ncol=2)\n    plt.tight_layout()\n    plt.savefig(os.path.join(cmp_output_dir, f\"{cmp}_water_occupancy_van{van}.png\"))\n
# Extract vanilla replicate and window ID\nfeatures_df[\"van\"] = features_df[\"ensemble\"].str.extract(r\"van_(\\d)\")\nfeatures_df[\"id\"] = features_df[\"ensemble\"].str.replace(r\"_van_\\d+\", \"\", regex=True)\n\n# Compute average occupancy\navg_df = features_df.groupby([\"id\", \"van\"], as_index=False)[\"occupancy\"].mean()\navg_df.to_csv(os.path.join(cmp_output_dir, \"avg_water_occupancy.csv\"), index=False)\n\n# Plot average occupancy\nplt.figure(figsize=(12, 8))\nsns.lineplot(data=avg_df, x=\"id\", y=\"occupancy\", hue=\"van\", palette=\"tab10\")\nplt.title(f\"Average Water Occupancy Across Windows for {cmp}\", fontsize=16, fontweight=\"bold\")\nplt.xlabel(\"Window ID\", fontsize=14)\nplt.ylabel(\"Average Number of Waters\", fontsize=14)\nplt.legend(title=\"Vanilla Repeat\", title_fontsize=12, fontsize=10, loc=\"upper right\")\nplt.xticks(rotation=45, fontsize=10)\nplt.yticks(fontsize=10)\nplt.tight_layout()\nplt.savefig(os.path.join(cmp_output_dir, f\"{cmp}_water_occupancy_across_windows.png\"), dpi=300)\nplt.close()\n
","path":["How-to-guides","Binding Pocket Resolvation Analysis"],"tags":[]},{"location":"how-to-guides/guide-2/","level":1,"title":"How-to-Guide","text":"

Put the example analysis in the how to guide section here

","path":["How-to-guides","How-to-Guide"],"tags":[]},{"location":"tutorials/tutorial-1/","level":1,"title":"Featurizing ABFEs","text":"

This tutorial is designed to the user to introduce fepa by analyzing MD data from a set of simulations of the ligand 42922 (from the Deflorian et al set 1 datset) bound to the Orexin 2 Receptor. Specifically, we will compare the holo, apo and ABFE trajectories. Through this comparison, we will demonstrate that the protein’s binding pocket adopts distinct conformations between the final lambda windows (where the ligand is fully annihilated) of the ABFE simulation and the long apo simulation of the receptor.

This tutorial is followed by Tutorial 2 where we will use FEPA to set up REUS simulations to estimate the free energy of this conformational change from the holo like to apo like states. This tutorial uses MDAnalysis, MDtraj and the PENSA package for analysis, GROMACS for simulation and Plumed for enhanced sampling. It is assumed that the user is familiar with setting up and analyzing MD simulations run with GROMACS and Plumed.

","path":["Tutorials","Featurizing ABFEs"],"tags":[]},{"location":"tutorials/tutorial-1/#loading-the-config-file","level":2,"title":"Loading the config file","text":"

ABFE simulations produce multiple MD trajectories over different lamba windows. In our case, we have 44 lambda windows (11 Coulomb, 12 Restrains and 22 Van-der-Walls) plus the holo and the apo simulations that we have to analyse. To make things easier when parsing paths to the topology and coordinate files of these trajectories, we use a config file. The config file is json-formatted and contains all the information necessary to read the simulation trajectories. Here is a sample config file that I use for this tutorial:

{\n    \"base_path\": \"deflorian_set_1_j13_v1\",\n    \"abfe_window_path_template\": \"deflorian_set_1_j13_v1/OX2_{CMP_NAME}/abfe_van{REP_NO}_hrex_r{ABFE_REP_NO}/complex/{LEG_WINDOW}/{STAGE}\",\n    \"vanilla_path_template\": \"deflorian_set_1_j13_v1/OX2_{CMP_NAME}/vanilla_rep_{REP_NO}\",\n    \"apo_path_template\": \"deflorian_set_1_j13_v1/apo_OX2_r{REP_NO}\",\n    \"compounds\": [\n        \"42922\",\n    ],\n    \"pocket_residues_string\": \"12  54  57  58  59  60  61  62  63  64  65  70  71  78  81  82  83  85  86  89 138 142 160 161 162 163 175 178 179 182 183 232 235 236 239 240 242 243 261 265 268 269\"\n}\n

The pocket_residues_string variable is a string that stores the residue ids of all the proteins residue that have any atom within 6 A of the ligand. The JSON file is then loaded into a dictionary with the function load_config from fepa.utils.file_utils

","path":["Tutorials","Featurizing ABFEs"],"tags":[]},{"location":"tutorials/tutorial-1/#loading-md-trajectories","level":2,"title":"Loading MD trajectories","text":"

Now that we have the template paths to all the simulations: Apo EQ, Holo EQ and ABFE, we can use the function load_abfe_paths_for_compound from fepa.utils.path_utils to generate a path_dict dictionary that contains all the MD file paths for a single compound as follows:

cmp = config['compounds'][0]\npath_dict = load_abfe_paths_for_compound(\n            config,\n            cmp,\n            van_list=[3], # Loading only vanilla simulation 3\n            leg_window_list=[f\"coul.{i:02d}\" for i in range(0, 11,2)] + [f'vdw.{i:02d}' for i in range(0, 21,2)] + [f'rest.{i:02d}' for i in range(0, 12,2)], # Loading only every other lambda window\n            bp_selection_string=\"name CA and resid \" + config[\"pocket_residues_string\"],\n            apo=True,\n        )\n

The load_abfe_paths_for_compound function has various arguments that allow the user control over what simulation paths are loaded. For more information, please refer to the API.

Now that we have the path_dict, it is time for us to load the trajectories themselves. We will be using the EnsembleHandler class from fepa.core.ensemble_handler to do this. EnsembleHandler is a neat way of storing and manipulating the trajectories from multiple ensembles with some built in functions for sanity checks. Internally EnsembleHandler stores the trajectories as dictionary of MDA universes.

TODO: Need to remove ensemblehandler and just use universe dict

# Load trajectories\nensemble_handler = EnsembleHandler(path_dict)\n# Make universes\nensemble_handler.make_universes()\nlogging.info(\"Making universes for compound %s...\", cmp)\n# Check for BP residue consistency across all the trajectories\nlogging.info(\"Checking residue consistency for compound %s...\", cmp)\n
","path":["Tutorials","Featurizing ABFEs"],"tags":[]},{"location":"tutorials/tutorial-1/#featurizing-md-trajectories","level":2,"title":"Featurizing MD trajectories","text":"

We will now be featurizing these trajectories by their pairwise C-alpha distances in the binding pocket using the function SelfDistanceFeaturizer from fepa.core.featurizers. SelfDistanceFeaturizer computes and stored all possible pairs of distances between the C-alpha atoms of the binding pocket residues. The class also has save_features and load_features functions that help save the features as a csv file to make sure the time consuming featurization step need not be repeated every run.

TODO: Make sure the binding pocket selection is consistent

# Make a folder for the analysis output\ncmp_run_dir = f'analysis/{cmp}/'\nif !(os.path.exists(cmp_run_dir)):\n    os.mkdir(cmp_run_dir)\n\n# Featurize and save features\nfeaturizer = SelfDistanceFeaturizer(ensemble_handler)\nfeaturizer.featurize()\nfeaturizer.save_features(input_dir=cmp_existing_run_dir)\n
","path":["Tutorials","Featurizing ABFEs"],"tags":[]},{"location":"tutorials/tutorial-1/#visualizing-ensembles","level":2,"title":"Visualizing ensembles","text":"

Saving the features in a csv format gives us the flexibility to analyse it as required. For the purpose of this tutorial, we will be looking at how our features capture the difference between the apo, the holo and the abfe ensembles. To do this we reduce the dimensions of the features data using the PCADimReducer class from fepa.core.dim_reducers. FEPA also supports other dimensionality reduction techniqeus like UMAP and tSNE. In fact UMAP is better able to resolve the differences in binding pocket configurations between different ensembles. We will be doing PCA here as it also doubles up as a nice CV to bias when performing umbrella sampling later.

# Dimensionality Reduction\nlogging.info(\"Performing dimensionality reduction for compound %s...\", cmp)\ndimreducer = PCADimReducer(featurizer.get_feature_df(), n_components=8)\ndimreducer.reduce_dimensions()\ndimreducer.calculate_projections()\ndimreducer.save_projection_df(\n    save_path=os.path.join(cmp_output_dir, \"pca_projection_df.csv\")\n)\n

First we plot the eigen values of all the PCs to understand what percentage of variance is captured by the first few PCs:

logging.info(\"Plotting PCA eigenvalues for compound %s...\", cmp)\nplot_eigenvalues(\n    pca_object=dimreducer.get_pca(),\n    n_components=8,\n    save_path=os.path.join(cmp_output_dir, \"eigenvalues.png\"),\n)\n

Figure 1: PCA eigenvalues for the binding-pocket Cα self-distance features.

Here, most of the variance is captured by PC1. However, this may not always be the case. Therefore, for this tutorial, we will use PC1 and PC2 together, as they collectively capture the majority of the variance. To improve visualization, we define a new column in projection_df called simtype, which groups the ensembles into broader categories: abfe_coul, abfe_vdw, abfe_rest, holo_equil, apo, and nvt.

# Further labelling the ensembles\nprojection_df = dimreducer.get_pca_projection_df()\n# Add another column called simtype based on ensemble\ndef get_simtype(ensemble_name: str) -> str:\n    if \"van\" in ensemble_name:\n        if 'nvt' in ensemble_name:\n            return \"nvt\"\n        if \"coul\" in ensemble_name:\n            return \"abfe_coul\"\n        elif \"vdw\" in ensemble_name:\n            return \"abfe_vdw\"\n        elif \"rest\" in ensemble_name:\n            return \"abfe_rest\"\n        else:\n            return \"holo_equil\"\n    elif \"apo\" in ensemble_name:\n        return \"apo\"\n    else:\n        return \"other\"\nprojection_df[\"simtype\"] = projection_df[\"ensemble\"].apply(get_simtype)\n

We can visualize the dimensionality-reduced data using the DimRedVisualizer class. In the example below, we plot the first two principal components, colored by simulation and time. We also highlight NVT, as it represents the initial crystal structure after energy minimization.

logging.info(\"Visualizing compound %s...\", cmp)\ndimred_visualizer = DimRedVisualizer(projection_df=projection_df, data_name=\"PCA\")\ndimred_visualizer.plot_dimred_sims(\n    save_path=os.path.join(cmp_run_dir, \"pca_components_ensemble_noapo.png\"),\n    highlights=[f\"{cmp}_nvt\"],\n)\ndimred_visualizer.plot_dimred_time(\n    save_path=os.path.join(cmp_run_dir, \"pca_components_time_noapo.png\")\n)\ndimred_visualizer.plot_dimred_sims(\n    column=\"simtype\",\n    save_path=os.path.join(cmp_run_dir, \"pca_components_simtype.png\"),\n)\n

Figure 2: Simulation frames in PC-space (PC1 vs PC2), colored by simtype

","path":["Tutorials","Featurizing ABFEs"],"tags":[]},{"location":"tutorials/tutorial-1/#clustering-ensembles","level":2,"title":"Clustering ensembles","text":"

In Figure 2, frames from the ABFE simulations closely resemble those from the holo equilibrium simulation, forming two distinct clusters: apo-like and holo-like. We will be clustering them using cluster_pca. After clustering, we can estimtate the data point that is closest to the center with the function make_ensemble_center_df which returns a DataFrame containing the details of the frame closest to the centroid in PC space. 

# Cluster the projection df\npca_projection_df_clustered = cluster_pca(\n    projection_df, n_clusters=3, n_components=3\n)\n# Ensemble center df\nensemble_center_df = make_ensemble_center_df(\n    pca_projection_df_clustered, key=\"cluster\"\n)\n

We then visualize the clustered data with ensemble centers as follows:

# Visualization of clustered data\ndimred_visualizer_clustered = DimRedVisualizer(\n    projection_df=pca_projection_df_clustered, data_name=\"PCA\"\n)\ndimred_visualizer_clustered.plot_dimred_cluster(\n    save_path=os.path.join(\n        cmp_run_dir, \"subset_pca_components_clusters_w_center.png\"\n    ),\n    centroid_df=ensemble_center_df,\n    cluster_column=\"cluster\",\n)\n

Figure 3: Clustered PCA data with ensemble centers marked by an X.

","path":["Tutorials","Featurizing ABFEs"],"tags":[]},{"location":"tutorials/tutorial-1/#visualizing-differences-across-clustered-ensembles","level":2,"title":"Visualizing differences across clustered ensembles","text":"

With two distinct clusters identified, we can compare them across different features by computing the Jensen–Shannon (JS) entropy for each feature. The code below plots feature-wise histograms for cluster 0 and cluster 1, annotated with their JS entropy values:

# Copy the cluster labels from pca_projection_df_clustered to feature_df\nfeature_df_w_clusters = featurizer.get_feature_df().copy()\nfeature_df_w_clusters[\"cluster\"] = pca_projection_df_clustered[\"cluster\"].values\n\n\n# Compute feature-level histograms and JS entropy between clusters 0 and 1\nhistograms = compute_histograms(\n    feature_df_w_clusters,\n    \"cluster\",\n    0,\n    1,\n    num_bins=50,\n    feature_column_keyword=\"DIST\",\n)\n\nrel_ent_dict = compute_relative_entropy(\n    feature_df_w_clusters,\n    ensemble1=0,\n    ensemble2=1,\n    num_bins=50,\n    key=\"cluster\",\n    feature_column_keyword=\"DIST\",\n)\n\nplot_jsd_histograms(\n    histograms,\n    rel_ent_dict,\n    top_n=16,\n    save_path=os.path.join(cmp_run_dir, \"jsd_cluster0_vs_cluster1.png\"),\n)\n

Figure 4: Histograms of the top 16 features showing the highest Jensen–Shannon divergence between clusters 0 and 1.

We then save the ensemble center frames as GRO files to compare structures in molecular visualization tools such as PyMOL:

# Save ensemble center frames as GRO files\nfor _, row in ensemble_center_df.iterrows():\n    center_ensemble = row[\"ensemble\"]\n    center_timestep = row[\"timestep\"]\n    print(f\"Ensemble: {center_ensemble}, Timestep: {center_timestep}\")\n\n    # Load trajectories\n    ensemble_handler.make_universes()\n\n    # Export the frame corresponding to the ensemble center\n    ensemble_handler.dump_frames(\n        ensemble=center_ensemble,\n        timestep=center_timestep,\n        save_path=f\"subset_cluster_{int(row['cluster'])}_center.gro\",\n    )\n
","path":["Tutorials","Featurizing ABFEs"],"tags":[]},{"location":"tutorials/tutorial-2/","level":1,"title":"Setting up REUS","text":"

In Tutorial 2, we build on the concepts learned in Tutorial 1, where the apo and holo states were shown to form distinct conformational clusters. Using the data provided in the [template], we will set up Replica Exchange Umbrella Sampling (REUS) simulations to estimate the free energy of the conformational transition from holo to apo, using the principal component (PC) 1 as the collective variable (CV). We will achieve this by morphing the protein from one state to another, equilibrating the system, and then running REUS to sample the transition.

","path":["Tutorials","Setting up REUS"],"tags":[]},{"location":"tutorials/tutorial-2/#preparing-the-plumeddat-file","level":2,"title":"Preparing the plumed.dat file","text":"

The first step in performing umbrella sampling with GROMACS and PLUMED is to create a plumed.dat file that defines the collective variable (CV), and the positions of the umbrella sampling restraints. To do this, we load the feature CSV files, initialize a DimReducer object, and extract the corresponding PCA object to define the CV based on pairwise Cα distances.

# Load Features df and reduce dimensions\nfeature_df = pd.read_csv(os.path.join(\"md_data\", \"top_features_apo_vdw.20.csv\"))\ndimreducer = PCADimReducer(feature_df, n_components=8)\ndimreducer.reduce_dimensions()\ndimreducer.calculate_projections()\npca_projection_df = dimreducer.get_pca_projection_df()\n# Use PCA from dimreducer to write plumed file\ntop_features = feature_df.filter(regex=\"DIST\", axis=1).columns\nwrite_plumed_file(\n    sdf_names=top_features,\n    top_features_pca=dimreducer.get_pca(),\n    save_path=\"plumed.dat\",\n    molinfo_structure=\"../reference.pdb\",  # fix molinfo here\n)\n

The PLUMED input file should look like this:

MOLINFO STRUCTURE=../reference.pdb\nd1: DISTANCE ATOMS=@CA-236,@CA-243\nd2: DISTANCE ATOMS=@CA-175,@CA-243\nd3: DISTANCE ATOMS=@CA-138,@CA-243\n.\n.\n.\nd200: DISTANCE ATOMS=@CA-61,@CA-70\n# Create the dot product\ndot: COMBINE ARG=d1,d2,d3...d200 COEFFICIENTS=0.095,0.152,0.133...0.057,0.054,0.032 PERIODIC=NO\nCV: MATHEVAL ARG=dot FUNC=10*x PERIODIC=NO\nPRINT ARG=CV FILE=COLVAR STRIDE=1\n# Put position of restraints here for each window\nrestraint: RESTRAINT ARG=CV AT=@replicas:$RESTRAINT_ARRAY KAPPA=$KAPPA\nPRINT ARG=restraint.* FILE=restr\n
Note that $RESTRAINT_ARRAY is a placeholder for the harmonic restraint positions in CV space. When defining these positions, ensure they align with the direction of the structural morph. By convention, we transition from vdw.20 to apo structures. To assign restraints correctly, identify the ensemble with the lowest PC value and order the restraints from min to max (or vice versa) accordingly.

# Pair of ensembles to compare\npair = (\"vdw.20\", \"apo\")\n# Get mean PC1 for the two ensembles\nmean_pc1_ensemble1 = pca_projection_df[pca_projection_df[\"state\"] == pair[0]][\n    \"PC1\"\n].mean()\nmean_pc1_ensemble2 = pca_projection_df[pca_projection_df[\"state\"] == pair[1]][\n    \"PC1\"\n].mean()\nlogging.info(\n    f\"Mean PC1 for {pair[0]}: {mean_pc1_ensemble1}, Mean PC1 for {pair[1]}: {mean_pc1_ensemble2}\"\n)\n# Get the restraint array based on the two ensembles\nPC1_min = pca_projection_df[\"PC1\"].min()\nPC1_max = pca_projection_df[\"PC1\"].max()\n# If min PC1 is from ensemble 1, then the restraint array should be from PC1_min to PC1_max\nif mean_pc1_ensemble1 < mean_pc1_ensemble2:\n    restraint_array = np.linspace(PC1_min, PC1_max, 24)\nelif mean_pc1_ensemble1 > mean_pc1_ensemble2:\n    restraint_array = np.linspace(PC1_max, PC1_min, 24)\n

Once the restraint array is prepared, we write it to the PLUMED file using the write_plumed_restraints function:

# Write the restraint array to the plumed file\nwrite_plumed_restraints(\n    plumed_file=\"plumed.dat\",\n    restraint_centers=restraint_array,\n    kappa=5,\n)\n

The restrain line on plumed.dat file must now be an array of distance restraints that looks like this:

restraint: RESTRAINT ARG=CV AT=@replicas:116.856,119.047,121.239...167.262 KAPPA=5\n
","path":["Tutorials","Setting up REUS"],"tags":[]},{"location":"tutorials/tutorial-2/#morphing-with-memento","level":2,"title":"Morphing with memento","text":"

For this tutorial, the GRO files representing the vdw.20 and apo states are provided in the data folder, prepared in the same way as in Tutorial 1. The plumed.dat file has also been prepared. The next step for REUS is to generate intermediate protein conformations using Memento (JCTC, 2023) . FEPA provides a class for easy access to Memento.

Using FEPA’s memento_workflow class, we can perform the protein morphing. A Memento directory is needed to store morphs for each entry in initial_target_gro_dict, and an apo_template path must be provided, containing the topology and MDP files required for equilibrium simulations.

# Declaring variables:\nmemento_dir = \"/biggin/b211/reub0138/Projects/orexin/deflorian_set_1_j13_v1_memento\"\n# Defining pairs\npair = (\"all_vdw20\", \"apo_3\")\nall_vdw20_gro_file = \"/biggin/b211/reub0138/Projects/orexin/deflorian_set_1_j13_v1/analysis/a12_dimreduce_vdw.20/subset_cluster_1_center.gro\"\napo_3_gro_file = \"/biggin/b211/reub0138/Projects/orexin/deflorian_set_1_j13_v1/analysis/a12_dimreduce_vdw.20/subset_cluster_2_center.gro\"\ntemplate_path = \"/biggin/b211/reub0138/Projects/orexin/deflorian_set_1_j13_v1_memento/apo_template\"  # Make sure no water and no ions in the template topology file\n# Setup Memento folders\nmemento_flow = memento_workflow(\n    memento_dir=memento_dir,\n    initial_gro=all_vdw20_gro_file,\n    target_gro=apo_3_gro_file,\n    initial_name=pair[0],\n    target_name=pair[1],\n    template_path=template_path,\n    run_name=\"memento_run_v1\",\n    n_residues=296,\n)\n

We can use the workflow functions to perform each step. First, prepare_memento sets up the files required by Memento. Then, run_memento executes the morphing. run_memento requires the protonation states of all histidines in the format expected by GROMACS pdb2gmx (0 for HID, 1 for HIE, 2 for HIP, 3 for HIS1). Additionally, the indices of CYX residues must be provided, as Memento cannot process them automatically.

# Preparing memento input\nmemento_flow.prepare_memento()\n# Running memento\nmemento_flow.run_memento(\n    template_path=template_path,\n    last_run=\"memento_run_v1\",\n    protonation_states=[1, 1, 1, 2, 1, 2],\n    cyx_residue_indices = []\n)\n
","path":["Tutorials","Setting up REUS"],"tags":[]},{"location":"tutorials/tutorial-2/#running-equilibration","level":2,"title":"Running Equilibration","text":"

This step may take some time. Once complete, we can set up equilibration simulations with Memento by providing a job script path, which can be modified to suit the specific HPC system.

# Running analysis\nmemento_flow.prepare_equil_simulations(\n    job_script_template=\"/biggin/b211/reub0138/Projects/orexin/lenselink_a2a_memento_v1/job_vanilla_ranv_equil_arr_template.sh\"\n)\n

After running prepare_equil_simulations, individual boxes for each morph are created. We then simulate each box to relax the side chains by submitting the job to our HPC. Once complete, we can proceed to the REUS simulations.

","path":["Tutorials","Setting up REUS"],"tags":[]},{"location":"tutorials/tutorial-2/#preparing-reus","level":2,"title":"Preparing REUS","text":"

Once the equilibration simulations are complete, use FEPA’s reus_umbrella_sampling_workflow class to set up umbrella sampling in the same working directory. The number of windows can be adjusted via n_windows (24 is typically sufficient). If the residue IDs in the plumed.dat file does not match those in the MD GRO files, FEPA requires the correct offset to be set via the plumed_resid_offset parameter.

import logging\nfrom pathlib import Path\nfrom fepa.flows.reus_flows import (\n    reus_umbrella_sampling_workflow,\n)\n\nsim_path = \"wdir/all_vdw20_apo_3/memento_run_v1/wdir/boxes/sim0\"\nwdir_path = Path(sim_path).parents[1]\nplumed_path = \"plumed.dat\"\nsimname = \"all_vdw20_apo_3\"\ninitial_gro = \"md_data/cluster_1_center.gro\"\nfinal_gro = \"md_data/cluster_2_center.gro\"\nsubmission_script_template_arr = \"md_data/job_reus_template.sh\"\nprint(f\"wdir_path: {wdir_path}, plumed_path: {plumed_path}\")\numbrella_sampler = reus_umbrella_sampling_workflow(\n    wdir_path=wdir_path,\n    plumed_path=plumed_path,\n    submission_script_template_arr=submission_script_template_arr,\n    start=simname.split(\"_\")[0] + \"_\" + simname.split(\"_\")[1],\n    end=simname.split(\"_\")[2] + \"_\" + simname.split(\"_\")[3],\n    reus_folder_name=\"reus_v1\",\n    n_windows=24,\n    plumed_resid_offset=0,\n    initial_gro=initial_gro,\n    target_gro=final_gro,\n)\numbrella_sampler.setup_simulations(exist_ok=True)\n
","path":["Tutorials","Setting up REUS"],"tags":[]},{"location":"tutorials/tutorial-2/#analyzing-reus","level":2,"title":"Analyzing REUS","text":"

Once your REUS system is set up, you can analyze the results using the workflow:

umbrella_sampler.prepare_wham()\numbrella_sampler.run_wham()\numbrella_sampler.analyse_us_hist(range=(90, 180), colvar_prefix=\"COLVAR\")\numbrella_sampler.get_initial_final_CVs()\numbrella_sampler.plot_free_energies(\n    units=\"kcal\",\n)\n

This prepares and runs WHAM on the histograms and generates the free energy curves.

Figure 1: Free energy curves from the tutorial simulations. The CV values of apo and holo structures are marked with gray dotted lines. Different lines represent curves computed using varying proportions of data to assess convergence.

","path":["Tutorials","Setting up REUS"],"tags":[]}]} \ No newline at end of file From 74e07dbaaf8825baad92c21396d65b7b5e599f19 Mon Sep 17 00:00:00 2001 From: nithishwer Date: Thu, 13 Nov 2025 10:44:05 +0000 Subject: [PATCH 18/40] fix: remove python from requirement.yml --- requirements.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/requirements.yml b/requirements.yml index 93d3c05..f464895 100644 --- a/requirements.yml +++ b/requirements.yml @@ -4,7 +4,6 @@ channels: - conda-forge - omnia dependencies: - - python>=3.10 - plumed - modeller - mdaencore From f8a9a64351ab9e4c5adf4d9be78623bf4faced0c Mon Sep 17 00:00:00 2001 From: nithishwer Date: Thu, 13 Nov 2025 10:45:13 +0000 Subject: [PATCH 19/40] fix: remove prefix from requirement.yml --- requirements.yml | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/requirements.yml b/requirements.yml index f464895..3ccee4e 100644 --- a/requirements.yml +++ b/requirements.yml @@ -12,5 +12,4 @@ dependencies: - cython - openmm - openmmforcefields - - numba -prefix: /biggin/b211/reub0138/miniforge3/envs/fepa + - numba \ No newline at end of file From 446f28509145cd444dd2324a699dda390dd5cf83 Mon Sep 17 00:00:00 2001 From: nithishwer Date: Thu, 13 Nov 2025 11:39:00 +0000 Subject: [PATCH 20/40] fix: remove python vers from test.yml --- .github/workflows/test.yml | 1 - docs/how-to-guides/guide-1.md | 5 ----- docs/how-to-guides/guide-2.md | 2 +- site/how-to-guides/guide-1/index.html | 4 ---- site/how-to-guides/guide-2/index.html | 6 +++--- site/search.json | 2 +- 6 files changed, 5 insertions(+), 15 deletions(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 9d70c02..d124fc5 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -23,7 +23,6 @@ jobs: - name: Set up Conda environment uses: conda-incubator/setup-miniconda@v3 with: - python-version: 3.10 environment-file: requirements.yml activate-environment: fepa_env auto-update-conda: true diff --git a/docs/how-to-guides/guide-1.md b/docs/how-to-guides/guide-1.md index 3f1fdf2..b7c2526 100644 --- a/docs/how-to-guides/guide-1.md +++ b/docs/how-to-guides/guide-1.md @@ -3,11 +3,6 @@ icon: lucide/compass title: Binding Pocket Resolvation Analysis --- ---- -icon: lucide/compass -title: Binding Pocket Resolvation Analysis ---- - Visualizing binding pocket resolvation across ABFE simulations is straightforward with FEPA. We follow a workflow similar to Tutorial 1, but instead of `SelfDistanceFeaturizer`, we use `BPWaterFeaturizer` to quantify water occupancy. ### Load Config and Prepare Paths diff --git a/docs/how-to-guides/guide-2.md b/docs/how-to-guides/guide-2.md index f2a3cc1..fccbd24 100644 --- a/docs/how-to-guides/guide-2.md +++ b/docs/how-to-guides/guide-2.md @@ -1,6 +1,6 @@ --- icon: lucide/compass -title: How-to-Guide +title: Side Chain Torsion Analysis --- Put the example analysis in the how to guide section here \ No newline at end of file diff --git a/site/how-to-guides/guide-1/index.html b/site/how-to-guides/guide-1/index.html index 76b5afc..a3a6864 100644 --- a/site/how-to-guides/guide-1/index.html +++ b/site/how-to-guides/guide-1/index.html @@ -779,10 +779,6 @@

Binding Pocket Resolvation Analysis

-
-

icon: lucide/compass -title: Binding Pocket Resolvation Analysis

-

Visualizing binding pocket resolvation across ABFE simulations is straightforward with FEPA. We follow a workflow similar to Tutorial 1, but instead of SelfDistanceFeaturizer, we use BPWaterFeaturizer to quantify water occupancy.

Load Config and Prepare Paths

import logging, os
diff --git a/site/how-to-guides/guide-2/index.html b/site/how-to-guides/guide-2/index.html
index 7e72c81..9b5ae4a 100644
--- a/site/how-to-guides/guide-2/index.html
+++ b/site/how-to-guides/guide-2/index.html
@@ -25,7 +25,7 @@
     
     
       
-        How-to-Guide - FEPA Documentation
+        Side Chain Torsion Analysis - FEPA Documentation
       
     
     
@@ -141,7 +141,7 @@
         

           
             
-              How-to-Guide
+              Side Chain Torsion Analysis
             
           
         

@@ -658,7 +658,7 @@
                 
                   
 
-  
How-to-Guide

+  
Side Chain Torsion Analysis

 
 
Put the example analysis in the how to guide section here

 
diff --git a/site/search.json b/site/search.json
index b20bd14..6711c71 100644
--- a/site/search.json
+++ b/site/search.json
@@ -1 +1 @@
-{"config":{"separator":"[\\s\\-_,:!=\\[\\]()\\\\\"`/]+|\\.(?!\\d)"},"items":[{"location":"","level":1,"title":"Getting Started with FEPA","text":"
FEPA (Free Energy Perturbation Analysis) is a Python package for analyzing molecular dynamics (MD) trajectories from FEP simulations, particularly ABFEs. FEPA allows you to visualize conformational changes and set up simulations to correct free energy estimates.
 
This guide covers installation, basic usage, and key workflows.
","path":["Getting Started"],"tags":[]},{"location":"#installation","level":2,"title":"Installation","text":"
FEPA is must be installed from GitHub:
 
git clone https://github.com/Nithishwer/FEPA.git\ncd FEPA\npip install -e .\n
","path":["Getting Started"],"tags":[]},{"location":"explanation/","level":1,"title":"Explanation section","text":"
Lorem Ipsum
","path":["Explanation"],"tags":[]},{"location":"reference/","level":1,"title":"Reference","text":"
This part of the project documentation focuses on an information-oriented approach. Use it as a reference for the technical implementation of the calculator project code.
 
::: fepa.core.analyzers ::: fepa.core.dim_reducers ::: fepa.core.ensemble_handler ::: fepa.core.featurizers ::: fepa.core.visualizers
","path":["Reference"],"tags":[]},{"location":"how-to-guides/guide-1/","level":1,"title":"Binding Pocket Resolvation Analysis","text":"
icon: lucide/compass title: Binding Pocket Resolvation Analysis
 
Visualizing binding pocket resolvation across ABFE simulations is straightforward with FEPA. We follow a workflow similar to Tutorial 1, but instead of SelfDistanceFeaturizer, we use BPWaterFeaturizer to quantify water occupancy.
","path":["How-to-guides","Binding Pocket Resolvation Analysis"],"tags":[]},{"location":"how-to-guides/guide-1/#load-config-and-prepare-paths","level":3,"title":"Load Config and Prepare Paths","text":"
import logging, os\nfrom fepa.utils.file_utils import load_config\nfrom fepa.utils.path_utils import load_abfe_paths_for_compound\n\nlogging.basicConfig(level=logging.INFO, format=\"%(asctime)s - %(levelname)s - %(message)s\")\n\nconfig = load_config(\"../../config/config.json\")\nanalysis_output_dir = \"wdir\"\ncmp = config[\"compounds\"][0]\n\nlogging.info(\"Analyzing compound %s ...\", cmp)\ncmp_output_dir = os.path.join(analysis_output_dir, cmp)\nos.makedirs(cmp_output_dir, exist_ok=True)\n\nlogging.info(\"Loading paths for compound %s...\", cmp)\npath_dict = load_abfe_paths_for_compound(\n    config,\n    cmp,\n    van_list=[1, 2, 3],\n    leg_window_list=[f\"coul.{i:02d}\" for i in range(0, 11)]\n                     + [f\"rest.{i:02d}\" for i in range(0, 12)]\n                     + [f\"vdw.{i:02d}\" for i in range(0, 21)],\n    bp_selection_string=\"name CA and resid 57 58 61 64 83 84 87 88 91 92 173 177 218 221 235 238 239 242 243 246\",\n    apo=False,\n)\n
","path":["How-to-guides","Binding Pocket Resolvation Analysis"],"tags":[]},{"location":"how-to-guides/guide-1/#load-trajectories-and-featurize-waters","level":3,"title":"Load trajectories and featurize waters","text":"
from fepa.core.ensemble_handler import EnsembleHandler\nfrom fepa.core.featurizers import BPWaterFeaturizer\n\nlogging.info(\"Loading trajectories for compound %s ...\", cmp)\nensemble_handler = EnsembleHandler(path_dict)\nensemble_handler.make_universes()\n\nlogging.info(\"Featurizing binding pocket waters ...\")\nbp_water_featurizer = BPWaterFeaturizer(ensemble_handler=ensemble_handler)\nbp_water_featurizer.featurize(radius=10)  # Count waters within 10 Ã… of pocket COM\n\nlogging.info(\"Saving features for compound %s ...\", cmp)\nbp_water_featurizer.save_features(cmp_output_dir, overwrite=True)\n
","path":["How-to-guides","Binding Pocket Resolvation Analysis"],"tags":[]},{"location":"how-to-guides/guide-1/#plot-water-occupancy-over-time-or-windows","level":3,"title":"Plot Water Occupancy Over Time or Windows","text":"
import pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nfeatures_df = pd.read_csv(os.path.join(cmp_output_dir, \"WaterOccupancy_features.csv\"))\n\nfor van in [1, 2, 3]:\n    van_features_df = features_df[features_df[\"ensemble\"].str.contains(f\"van_{van}\")]\n    plt.figure(figsize=(12, 6))\n    sns.lineplot(data=van_features_df, x=\"Time (ps)\", y=\"occupancy\", hue=\"ensemble\")\n    plt.title(f\"Water Occupancy for {cmp}\")\n    plt.xlabel(\"Time (ps)\")\n    plt.ylabel(\"Number of Waters\")\n    plt.xlim(0, 20000)\n    plt.legend(title=\"Ensemble\", bbox_to_anchor=(1.05, 1), loc=\"upper left\", ncol=2)\n    plt.tight_layout()\n    plt.savefig(os.path.join(cmp_output_dir, f\"{cmp}_water_occupancy_van{van}.png\"))\n
 
# Extract vanilla replicate and window ID\nfeatures_df[\"van\"] = features_df[\"ensemble\"].str.extract(r\"van_(\\d)\")\nfeatures_df[\"id\"] = features_df[\"ensemble\"].str.replace(r\"_van_\\d+\", \"\", regex=True)\n\n# Compute average occupancy\navg_df = features_df.groupby([\"id\", \"van\"], as_index=False)[\"occupancy\"].mean()\navg_df.to_csv(os.path.join(cmp_output_dir, \"avg_water_occupancy.csv\"), index=False)\n\n# Plot average occupancy\nplt.figure(figsize=(12, 8))\nsns.lineplot(data=avg_df, x=\"id\", y=\"occupancy\", hue=\"van\", palette=\"tab10\")\nplt.title(f\"Average Water Occupancy Across Windows for {cmp}\", fontsize=16, fontweight=\"bold\")\nplt.xlabel(\"Window ID\", fontsize=14)\nplt.ylabel(\"Average Number of Waters\", fontsize=14)\nplt.legend(title=\"Vanilla Repeat\", title_fontsize=12, fontsize=10, loc=\"upper right\")\nplt.xticks(rotation=45, fontsize=10)\nplt.yticks(fontsize=10)\nplt.tight_layout()\nplt.savefig(os.path.join(cmp_output_dir, f\"{cmp}_water_occupancy_across_windows.png\"), dpi=300)\nplt.close()\n
","path":["How-to-guides","Binding Pocket Resolvation Analysis"],"tags":[]},{"location":"how-to-guides/guide-2/","level":1,"title":"How-to-Guide","text":"
Put the example analysis in the how to guide section here
","path":["How-to-guides","How-to-Guide"],"tags":[]},{"location":"tutorials/tutorial-1/","level":1,"title":"Featurizing ABFEs","text":"
This tutorial is designed to the user to introduce fepa by analyzing MD data from a set of simulations of the ligand 42922 (from the Deflorian et al set 1 datset) bound to the Orexin 2 Receptor. Specifically, we will compare the holo, apo and ABFE trajectories. Through this comparison, we will demonstrate that the protein’s binding pocket adopts distinct conformations between the final lambda windows (where the ligand is fully annihilated) of the ABFE simulation and the long apo simulation of the receptor.
 
This tutorial is followed by Tutorial 2 where we will use FEPA to set up REUS simulations to estimate the free energy of this conformational change from the holo like  to apo like states. This tutorial uses MDAnalysis, MDtraj and the PENSA package for analysis, GROMACS for simulation and Plumed for enhanced sampling. It is assumed that the user is familiar with setting up and analyzing MD simulations run with GROMACS and Plumed. 
","path":["Tutorials","Featurizing ABFEs"],"tags":[]},{"location":"tutorials/tutorial-1/#loading-the-config-file","level":2,"title":"Loading the config file","text":"
ABFE simulations produce multiple MD trajectories over different lamba windows. In our case, we have 44 lambda windows (11 Coulomb, 12 Restrains and 22 Van-der-Walls) plus the holo and the apo simulations that we have to analyse. To make things easier when parsing paths to the topology and coordinate files of these trajectories, we use a config file. The config file is json-formatted and contains all the information necessary to read the simulation trajectories. Here is a sample config file that I use for this tutorial:
 
{\n    \"base_path\": \"deflorian_set_1_j13_v1\",\n    \"abfe_window_path_template\": \"deflorian_set_1_j13_v1/OX2_{CMP_NAME}/abfe_van{REP_NO}_hrex_r{ABFE_REP_NO}/complex/{LEG_WINDOW}/{STAGE}\",\n    \"vanilla_path_template\": \"deflorian_set_1_j13_v1/OX2_{CMP_NAME}/vanilla_rep_{REP_NO}\",\n    \"apo_path_template\": \"deflorian_set_1_j13_v1/apo_OX2_r{REP_NO}\",\n    \"compounds\": [\n        \"42922\",\n    ],\n    \"pocket_residues_string\": \"12  54  57  58  59  60  61  62  63  64  65  70  71  78  81  82  83  85  86  89 138 142 160 161 162 163 175 178 179 182 183 232 235 236 239 240 242 243 261 265 268 269\"\n}\n
 
The pocket_residues_string variable is a string that stores the residue ids of all the proteins residue that have any atom within 6 A of the ligand. The JSON file is then loaded into a dictionary with the function load_config from fepa.utils.file_utils
","path":["Tutorials","Featurizing ABFEs"],"tags":[]},{"location":"tutorials/tutorial-1/#loading-md-trajectories","level":2,"title":"Loading MD trajectories","text":"
Now that we have the template paths to all the simulations: Apo EQ, Holo EQ and ABFE, we can use the function load_abfe_paths_for_compound from fepa.utils.path_utils to generate a path_dict dictionary that contains all the MD file paths for a single compound as follows:
 
cmp = config['compounds'][0]\npath_dict = load_abfe_paths_for_compound(\n            config,\n            cmp,\n            van_list=[3], # Loading only vanilla simulation 3\n            leg_window_list=[f\"coul.{i:02d}\" for i in range(0, 11,2)] + [f'vdw.{i:02d}' for i in range(0, 21,2)] + [f'rest.{i:02d}' for i in range(0, 12,2)], # Loading only every other lambda window\n            bp_selection_string=\"name CA and resid \" + config[\"pocket_residues_string\"],\n            apo=True,\n        )\n
 
The load_abfe_paths_for_compound function has various arguments that allow the user control over what simulation paths are loaded. For more information, please refer to the API. 
 
Now that we have the path_dict, it is time for us to load the trajectories themselves. We will be using the EnsembleHandler class from fepa.core.ensemble_handler to do this. EnsembleHandler is a neat way of storing and manipulating the trajectories from multiple ensembles with some built in functions for sanity checks. Internally EnsembleHandler stores the trajectories as dictionary of MDA universes.  
 
TODO: Need to remove ensemblehandler and just use universe dict
 
# Load trajectories\nensemble_handler = EnsembleHandler(path_dict)\n# Make universes\nensemble_handler.make_universes()\nlogging.info(\"Making universes for compound %s...\", cmp)\n# Check for BP residue consistency across all the trajectories\nlogging.info(\"Checking residue consistency for compound %s...\", cmp)\n
","path":["Tutorials","Featurizing ABFEs"],"tags":[]},{"location":"tutorials/tutorial-1/#featurizing-md-trajectories","level":2,"title":"Featurizing MD trajectories","text":"
We will now be featurizing these trajectories by their pairwise C-alpha distances in the binding pocket using the function SelfDistanceFeaturizer from fepa.core.featurizers. SelfDistanceFeaturizer computes and stored all possible pairs of distances between the C-alpha atoms of the binding pocket residues. The class also has save_features and load_features functions that help save the features as a csv file to make sure the time consuming featurization step need not be repeated every run.
 
TODO: Make sure the binding pocket selection is consistent
 
# Make a folder for the analysis output\ncmp_run_dir = f'analysis/{cmp}/'\nif !(os.path.exists(cmp_run_dir)):\n    os.mkdir(cmp_run_dir)\n\n# Featurize and save features\nfeaturizer = SelfDistanceFeaturizer(ensemble_handler)\nfeaturizer.featurize()\nfeaturizer.save_features(input_dir=cmp_existing_run_dir)\n
","path":["Tutorials","Featurizing ABFEs"],"tags":[]},{"location":"tutorials/tutorial-1/#visualizing-ensembles","level":2,"title":"Visualizing ensembles","text":"
Saving the features in a csv format gives us the flexibility to analyse it as required. For the purpose of this tutorial, we will be looking at how our features capture the difference between the apo, the holo and the abfe ensembles. To do this we reduce the dimensions of the features data using the PCADimReducer class from fepa.core.dim_reducers. FEPA also supports other dimensionality reduction techniqeus like UMAP and tSNE. In fact UMAP is better able to resolve the differences in binding pocket configurations between different ensembles. We will be doing PCA here as it also doubles up as a nice CV to bias when performing umbrella sampling later.
 
# Dimensionality Reduction\nlogging.info(\"Performing dimensionality reduction for compound %s...\", cmp)\ndimreducer = PCADimReducer(featurizer.get_feature_df(), n_components=8)\ndimreducer.reduce_dimensions()\ndimreducer.calculate_projections()\ndimreducer.save_projection_df(\n    save_path=os.path.join(cmp_output_dir, \"pca_projection_df.csv\")\n)\n
 
First we plot the eigen values of all the PCs to understand what percentage of variance is captured by the first few PCs:
 
logging.info(\"Plotting PCA eigenvalues for compound %s...\", cmp)\nplot_eigenvalues(\n    pca_object=dimreducer.get_pca(),\n    n_components=8,\n    save_path=os.path.join(cmp_output_dir, \"eigenvalues.png\"),\n)\n
 
 
Figure 1: PCA eigenvalues for the binding-pocket Cα self-distance features.
 
Here, most of the variance is captured by PC1. However, this may not always be the case. Therefore, for this tutorial, we will use PC1 and PC2 together, as they collectively capture the majority of the variance. To improve visualization, we define a new column in projection_df called simtype, which groups the ensembles into broader categories: abfe_coul, abfe_vdw, abfe_rest, holo_equil, apo, and nvt.
 
# Further labelling the ensembles\nprojection_df = dimreducer.get_pca_projection_df()\n# Add another column called simtype based on ensemble\ndef get_simtype(ensemble_name: str) -> str:\n    if \"van\" in ensemble_name:\n        if 'nvt' in ensemble_name:\n            return \"nvt\"\n        if \"coul\" in ensemble_name:\n            return \"abfe_coul\"\n        elif \"vdw\" in ensemble_name:\n            return \"abfe_vdw\"\n        elif \"rest\" in ensemble_name:\n            return \"abfe_rest\"\n        else:\n            return \"holo_equil\"\n    elif \"apo\" in ensemble_name:\n        return \"apo\"\n    else:\n        return \"other\"\nprojection_df[\"simtype\"] = projection_df[\"ensemble\"].apply(get_simtype)\n
 
We can visualize the dimensionality-reduced data using the DimRedVisualizer class. In the example below, we plot the first two principal components, colored by simulation and time. We also highlight NVT, as it represents the initial crystal structure after energy minimization.
 
logging.info(\"Visualizing compound %s...\", cmp)\ndimred_visualizer = DimRedVisualizer(projection_df=projection_df, data_name=\"PCA\")\ndimred_visualizer.plot_dimred_sims(\n    save_path=os.path.join(cmp_run_dir, \"pca_components_ensemble_noapo.png\"),\n    highlights=[f\"{cmp}_nvt\"],\n)\ndimred_visualizer.plot_dimred_time(\n    save_path=os.path.join(cmp_run_dir, \"pca_components_time_noapo.png\")\n)\ndimred_visualizer.plot_dimred_sims(\n    column=\"simtype\",\n    save_path=os.path.join(cmp_run_dir, \"pca_components_simtype.png\"),\n)\n
 
 Figure 2: Simulation frames in PC-space (PC1 vs PC2), colored by simtype
","path":["Tutorials","Featurizing ABFEs"],"tags":[]},{"location":"tutorials/tutorial-1/#clustering-ensembles","level":2,"title":"Clustering ensembles","text":"
In Figure 2, frames from the ABFE simulations closely resemble those from the holo equilibrium simulation, forming two distinct clusters: apo-like and holo-like. We will be clustering them using cluster_pca. After clustering, we can estimtate the data point that is closest to the center with the function make_ensemble_center_df which returns a DataFrame containing the details of the frame closest to the centroid in PC space. 
 
# Cluster the projection df\npca_projection_df_clustered = cluster_pca(\n    projection_df, n_clusters=3, n_components=3\n)\n# Ensemble center df\nensemble_center_df = make_ensemble_center_df(\n    pca_projection_df_clustered, key=\"cluster\"\n)\n
 
We then visualize the clustered data with ensemble centers as follows:
 
# Visualization of clustered data\ndimred_visualizer_clustered = DimRedVisualizer(\n    projection_df=pca_projection_df_clustered, data_name=\"PCA\"\n)\ndimred_visualizer_clustered.plot_dimred_cluster(\n    save_path=os.path.join(\n        cmp_run_dir, \"subset_pca_components_clusters_w_center.png\"\n    ),\n    centroid_df=ensemble_center_df,\n    cluster_column=\"cluster\",\n)\n
 
 
Figure 3: Clustered PCA data with ensemble centers marked by an X.
","path":["Tutorials","Featurizing ABFEs"],"tags":[]},{"location":"tutorials/tutorial-1/#visualizing-differences-across-clustered-ensembles","level":2,"title":"Visualizing differences across clustered ensembles","text":"
With two distinct clusters identified, we can compare them across different features by computing the Jensen–Shannon (JS) entropy for each feature. The code below plots feature-wise histograms for cluster 0 and cluster 1, annotated with their JS entropy values:
 
# Copy the cluster labels from pca_projection_df_clustered to feature_df\nfeature_df_w_clusters = featurizer.get_feature_df().copy()\nfeature_df_w_clusters[\"cluster\"] = pca_projection_df_clustered[\"cluster\"].values\n\n\n# Compute feature-level histograms and JS entropy between clusters 0 and 1\nhistograms = compute_histograms(\n    feature_df_w_clusters,\n    \"cluster\",\n    0,\n    1,\n    num_bins=50,\n    feature_column_keyword=\"DIST\",\n)\n\nrel_ent_dict = compute_relative_entropy(\n    feature_df_w_clusters,\n    ensemble1=0,\n    ensemble2=1,\n    num_bins=50,\n    key=\"cluster\",\n    feature_column_keyword=\"DIST\",\n)\n\nplot_jsd_histograms(\n    histograms,\n    rel_ent_dict,\n    top_n=16,\n    save_path=os.path.join(cmp_run_dir, \"jsd_cluster0_vs_cluster1.png\"),\n)\n
 
 Figure 4: Histograms of the top 16 features showing the highest Jensen–Shannon divergence between clusters 0 and 1.
 
We then save the ensemble center frames as GRO files to compare structures in molecular visualization tools such as PyMOL:
 
# Save ensemble center frames as GRO files\nfor _, row in ensemble_center_df.iterrows():\n    center_ensemble = row[\"ensemble\"]\n    center_timestep = row[\"timestep\"]\n    print(f\"Ensemble: {center_ensemble}, Timestep: {center_timestep}\")\n\n    # Load trajectories\n    ensemble_handler.make_universes()\n\n    # Export the frame corresponding to the ensemble center\n    ensemble_handler.dump_frames(\n        ensemble=center_ensemble,\n        timestep=center_timestep,\n        save_path=f\"subset_cluster_{int(row['cluster'])}_center.gro\",\n    )\n
","path":["Tutorials","Featurizing ABFEs"],"tags":[]},{"location":"tutorials/tutorial-2/","level":1,"title":"Setting up REUS","text":"
In Tutorial 2, we build on the concepts learned in Tutorial 1, where the apo and holo states were shown to form distinct conformational clusters. Using the data provided in the [template], we will set up Replica Exchange Umbrella Sampling (REUS) simulations to estimate the free energy of the conformational transition from holo to apo, using the principal component (PC) 1 as the collective variable (CV). We will achieve this by morphing the protein from one state to another, equilibrating the system, and then running REUS to sample the transition.
","path":["Tutorials","Setting up REUS"],"tags":[]},{"location":"tutorials/tutorial-2/#preparing-the-plumeddat-file","level":2,"title":"Preparing the plumed.dat file","text":"
The first step in performing umbrella sampling with GROMACS and PLUMED is to create a plumed.dat file that defines the collective variable (CV), and the positions of the umbrella sampling restraints. To do this, we load the feature CSV files, initialize a DimReducer object, and extract the corresponding PCA object to define the CV based on pairwise Cα distances.
 
# Load Features df and reduce dimensions\nfeature_df = pd.read_csv(os.path.join(\"md_data\", \"top_features_apo_vdw.20.csv\"))\ndimreducer = PCADimReducer(feature_df, n_components=8)\ndimreducer.reduce_dimensions()\ndimreducer.calculate_projections()\npca_projection_df = dimreducer.get_pca_projection_df()\n# Use PCA from dimreducer to write plumed file\ntop_features = feature_df.filter(regex=\"DIST\", axis=1).columns\nwrite_plumed_file(\n    sdf_names=top_features,\n    top_features_pca=dimreducer.get_pca(),\n    save_path=\"plumed.dat\",\n    molinfo_structure=\"../reference.pdb\",  # fix molinfo here\n)\n
 
The PLUMED input file should look like this:
 
MOLINFO STRUCTURE=../reference.pdb\nd1: DISTANCE ATOMS=@CA-236,@CA-243\nd2: DISTANCE ATOMS=@CA-175,@CA-243\nd3: DISTANCE ATOMS=@CA-138,@CA-243\n.\n.\n.\nd200: DISTANCE ATOMS=@CA-61,@CA-70\n# Create the dot product\ndot: COMBINE ARG=d1,d2,d3...d200 COEFFICIENTS=0.095,0.152,0.133...0.057,0.054,0.032 PERIODIC=NO\nCV: MATHEVAL ARG=dot FUNC=10*x PERIODIC=NO\nPRINT ARG=CV FILE=COLVAR STRIDE=1\n# Put position of restraints here for each window\nrestraint: RESTRAINT ARG=CV AT=@replicas:$RESTRAINT_ARRAY KAPPA=$KAPPA\nPRINT ARG=restraint.* FILE=restr\n
 Note that $RESTRAINT_ARRAY is a placeholder for the harmonic restraint positions in CV space. When defining these positions, ensure they align with the direction of the structural morph. By convention, we transition from vdw.20 to apo structures. To assign restraints correctly, identify the ensemble with the lowest PC value and order the restraints from min to max (or vice versa) accordingly.
 
# Pair of ensembles to compare\npair = (\"vdw.20\", \"apo\")\n# Get mean PC1 for the two ensembles\nmean_pc1_ensemble1 = pca_projection_df[pca_projection_df[\"state\"] == pair[0]][\n    \"PC1\"\n].mean()\nmean_pc1_ensemble2 = pca_projection_df[pca_projection_df[\"state\"] == pair[1]][\n    \"PC1\"\n].mean()\nlogging.info(\n    f\"Mean PC1 for {pair[0]}: {mean_pc1_ensemble1}, Mean PC1 for {pair[1]}: {mean_pc1_ensemble2}\"\n)\n# Get the restraint array based on the two ensembles\nPC1_min = pca_projection_df[\"PC1\"].min()\nPC1_max = pca_projection_df[\"PC1\"].max()\n# If min PC1 is from ensemble 1, then the restraint array should be from PC1_min to PC1_max\nif mean_pc1_ensemble1 < mean_pc1_ensemble2:\n    restraint_array = np.linspace(PC1_min, PC1_max, 24)\nelif mean_pc1_ensemble1 > mean_pc1_ensemble2:\n    restraint_array = np.linspace(PC1_max, PC1_min, 24)\n
 
Once the restraint array is prepared, we write it to the PLUMED file using the write_plumed_restraints function:
 
# Write the restraint array to the plumed file\nwrite_plumed_restraints(\n    plumed_file=\"plumed.dat\",\n    restraint_centers=restraint_array,\n    kappa=5,\n)\n
 
The restrain line on plumed.dat file must now be an array of distance restraints that looks like this:
 
restraint: RESTRAINT ARG=CV AT=@replicas:116.856,119.047,121.239...167.262 KAPPA=5\n
","path":["Tutorials","Setting up REUS"],"tags":[]},{"location":"tutorials/tutorial-2/#morphing-with-memento","level":2,"title":"Morphing with memento","text":"
For this tutorial, the GRO files representing the vdw.20 and apo states are provided in the data folder, prepared in the same way as in Tutorial 1. The plumed.dat file has also been prepared. The next step for REUS is to generate intermediate protein conformations using Memento (JCTC, 2023) . FEPA provides a class for easy access to Memento.
 
Using FEPA’s memento_workflow class, we can perform the protein morphing. A Memento directory is needed to store morphs for each entry in initial_target_gro_dict, and an apo_template path must be provided, containing the topology and MDP files required for equilibrium simulations.
 
# Declaring variables:\nmemento_dir = \"/biggin/b211/reub0138/Projects/orexin/deflorian_set_1_j13_v1_memento\"\n# Defining pairs\npair = (\"all_vdw20\", \"apo_3\")\nall_vdw20_gro_file = \"/biggin/b211/reub0138/Projects/orexin/deflorian_set_1_j13_v1/analysis/a12_dimreduce_vdw.20/subset_cluster_1_center.gro\"\napo_3_gro_file = \"/biggin/b211/reub0138/Projects/orexin/deflorian_set_1_j13_v1/analysis/a12_dimreduce_vdw.20/subset_cluster_2_center.gro\"\ntemplate_path = \"/biggin/b211/reub0138/Projects/orexin/deflorian_set_1_j13_v1_memento/apo_template\"  # Make sure no water and no ions in the template topology file\n# Setup Memento folders\nmemento_flow = memento_workflow(\n    memento_dir=memento_dir,\n    initial_gro=all_vdw20_gro_file,\n    target_gro=apo_3_gro_file,\n    initial_name=pair[0],\n    target_name=pair[1],\n    template_path=template_path,\n    run_name=\"memento_run_v1\",\n    n_residues=296,\n)\n
 
We can use the workflow functions to perform each step. First, prepare_memento sets up the files required by Memento. Then, run_memento executes the morphing. run_memento requires the protonation states of all histidines in the format expected by GROMACS pdb2gmx (0 for HID, 1 for HIE, 2 for HIP, 3 for HIS1). Additionally, the indices of CYX residues must be provided, as Memento cannot process them automatically.
 
# Preparing memento input\nmemento_flow.prepare_memento()\n# Running memento\nmemento_flow.run_memento(\n    template_path=template_path,\n    last_run=\"memento_run_v1\",\n    protonation_states=[1, 1, 1, 2, 1, 2],\n    cyx_residue_indices = []\n)\n
","path":["Tutorials","Setting up REUS"],"tags":[]},{"location":"tutorials/tutorial-2/#running-equilibration","level":2,"title":"Running Equilibration","text":"
This step may take some time. Once complete, we can set up equilibration simulations with Memento by providing a job script path, which can be modified to suit the specific HPC system.
 
# Running analysis\nmemento_flow.prepare_equil_simulations(\n    job_script_template=\"/biggin/b211/reub0138/Projects/orexin/lenselink_a2a_memento_v1/job_vanilla_ranv_equil_arr_template.sh\"\n)\n
 
After running prepare_equil_simulations, individual boxes for each morph are created. We then simulate each box to relax the side chains by submitting the job to our HPC. Once complete, we can proceed to the REUS simulations.
","path":["Tutorials","Setting up REUS"],"tags":[]},{"location":"tutorials/tutorial-2/#preparing-reus","level":2,"title":"Preparing REUS","text":"
Once the equilibration simulations are complete, use FEPA’s reus_umbrella_sampling_workflow class to set up umbrella sampling in the same working directory. The number of windows can be adjusted via n_windows (24 is typically sufficient). If the residue IDs in the plumed.dat file does not match those in the MD GRO files, FEPA requires the correct offset to be set via the plumed_resid_offset parameter.
 
import logging\nfrom pathlib import Path\nfrom fepa.flows.reus_flows import (\n    reus_umbrella_sampling_workflow,\n)\n\nsim_path = \"wdir/all_vdw20_apo_3/memento_run_v1/wdir/boxes/sim0\"\nwdir_path = Path(sim_path).parents[1]\nplumed_path = \"plumed.dat\"\nsimname = \"all_vdw20_apo_3\"\ninitial_gro = \"md_data/cluster_1_center.gro\"\nfinal_gro = \"md_data/cluster_2_center.gro\"\nsubmission_script_template_arr = \"md_data/job_reus_template.sh\"\nprint(f\"wdir_path: {wdir_path}, plumed_path: {plumed_path}\")\numbrella_sampler = reus_umbrella_sampling_workflow(\n    wdir_path=wdir_path,\n    plumed_path=plumed_path,\n    submission_script_template_arr=submission_script_template_arr,\n    start=simname.split(\"_\")[0] + \"_\" + simname.split(\"_\")[1],\n    end=simname.split(\"_\")[2] + \"_\" + simname.split(\"_\")[3],\n    reus_folder_name=\"reus_v1\",\n    n_windows=24,\n    plumed_resid_offset=0,\n    initial_gro=initial_gro,\n    target_gro=final_gro,\n)\numbrella_sampler.setup_simulations(exist_ok=True)\n
","path":["Tutorials","Setting up REUS"],"tags":[]},{"location":"tutorials/tutorial-2/#analyzing-reus","level":2,"title":"Analyzing REUS","text":"
Once your REUS system is set up, you can analyze the results using the workflow:
 
umbrella_sampler.prepare_wham()\numbrella_sampler.run_wham()\numbrella_sampler.analyse_us_hist(range=(90, 180), colvar_prefix=\"COLVAR\")\numbrella_sampler.get_initial_final_CVs()\numbrella_sampler.plot_free_energies(\n    units=\"kcal\",\n)\n
 
This prepares and runs WHAM on the histograms and generates the free energy curves.
 
 
Figure 1: Free energy curves from the tutorial simulations. The CV values of apo and holo structures are marked with gray dotted lines. Different lines represent curves computed using varying proportions of data to assess convergence.
","path":["Tutorials","Setting up REUS"],"tags":[]}]}
\ No newline at end of file
+{"config":{"separator":"[\\s\\-_,:!=\\[\\]()\\\\\"`/]+|\\.(?!\\d)"},"items":[{"location":"","level":1,"title":"Getting Started with FEPA","text":"
FEPA (Free Energy Perturbation Analysis) is a Python package for analyzing molecular dynamics (MD) trajectories from FEP simulations, particularly ABFEs. FEPA allows you to visualize conformational changes and set up simulations to correct free energy estimates.
 
This guide covers installation, basic usage, and key workflows.
","path":["Getting Started"],"tags":[]},{"location":"#installation","level":2,"title":"Installation","text":"
FEPA is must be installed from GitHub:
 
git clone https://github.com/Nithishwer/FEPA.git\ncd FEPA\npip install -e .\n
","path":["Getting Started"],"tags":[]},{"location":"explanation/","level":1,"title":"Explanation section","text":"
Lorem Ipsum
","path":["Explanation"],"tags":[]},{"location":"reference/","level":1,"title":"Reference","text":"
This part of the project documentation focuses on an information-oriented approach. Use it as a reference for the technical implementation of the calculator project code.
 
::: fepa.core.analyzers ::: fepa.core.dim_reducers ::: fepa.core.ensemble_handler ::: fepa.core.featurizers ::: fepa.core.visualizers
","path":["Reference"],"tags":[]},{"location":"how-to-guides/guide-1/","level":1,"title":"Binding Pocket Resolvation Analysis","text":"
Visualizing binding pocket resolvation across ABFE simulations is straightforward with FEPA. We follow a workflow similar to Tutorial 1, but instead of SelfDistanceFeaturizer, we use BPWaterFeaturizer to quantify water occupancy.
","path":["How-to-guides","Binding Pocket Resolvation Analysis"],"tags":[]},{"location":"how-to-guides/guide-1/#load-config-and-prepare-paths","level":3,"title":"Load Config and Prepare Paths","text":"
import logging, os\nfrom fepa.utils.file_utils import load_config\nfrom fepa.utils.path_utils import load_abfe_paths_for_compound\n\nlogging.basicConfig(level=logging.INFO, format=\"%(asctime)s - %(levelname)s - %(message)s\")\n\nconfig = load_config(\"../../config/config.json\")\nanalysis_output_dir = \"wdir\"\ncmp = config[\"compounds\"][0]\n\nlogging.info(\"Analyzing compound %s ...\", cmp)\ncmp_output_dir = os.path.join(analysis_output_dir, cmp)\nos.makedirs(cmp_output_dir, exist_ok=True)\n\nlogging.info(\"Loading paths for compound %s...\", cmp)\npath_dict = load_abfe_paths_for_compound(\n    config,\n    cmp,\n    van_list=[1, 2, 3],\n    leg_window_list=[f\"coul.{i:02d}\" for i in range(0, 11)]\n                     + [f\"rest.{i:02d}\" for i in range(0, 12)]\n                     + [f\"vdw.{i:02d}\" for i in range(0, 21)],\n    bp_selection_string=\"name CA and resid 57 58 61 64 83 84 87 88 91 92 173 177 218 221 235 238 239 242 243 246\",\n    apo=False,\n)\n
","path":["How-to-guides","Binding Pocket Resolvation Analysis"],"tags":[]},{"location":"how-to-guides/guide-1/#load-trajectories-and-featurize-waters","level":3,"title":"Load trajectories and featurize waters","text":"
from fepa.core.ensemble_handler import EnsembleHandler\nfrom fepa.core.featurizers import BPWaterFeaturizer\n\nlogging.info(\"Loading trajectories for compound %s ...\", cmp)\nensemble_handler = EnsembleHandler(path_dict)\nensemble_handler.make_universes()\n\nlogging.info(\"Featurizing binding pocket waters ...\")\nbp_water_featurizer = BPWaterFeaturizer(ensemble_handler=ensemble_handler)\nbp_water_featurizer.featurize(radius=10)  # Count waters within 10 Ã… of pocket COM\n\nlogging.info(\"Saving features for compound %s ...\", cmp)\nbp_water_featurizer.save_features(cmp_output_dir, overwrite=True)\n
","path":["How-to-guides","Binding Pocket Resolvation Analysis"],"tags":[]},{"location":"how-to-guides/guide-1/#plot-water-occupancy-over-time-or-windows","level":3,"title":"Plot Water Occupancy Over Time or Windows","text":"
import pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nfeatures_df = pd.read_csv(os.path.join(cmp_output_dir, \"WaterOccupancy_features.csv\"))\n\nfor van in [1, 2, 3]:\n    van_features_df = features_df[features_df[\"ensemble\"].str.contains(f\"van_{van}\")]\n    plt.figure(figsize=(12, 6))\n    sns.lineplot(data=van_features_df, x=\"Time (ps)\", y=\"occupancy\", hue=\"ensemble\")\n    plt.title(f\"Water Occupancy for {cmp}\")\n    plt.xlabel(\"Time (ps)\")\n    plt.ylabel(\"Number of Waters\")\n    plt.xlim(0, 20000)\n    plt.legend(title=\"Ensemble\", bbox_to_anchor=(1.05, 1), loc=\"upper left\", ncol=2)\n    plt.tight_layout()\n    plt.savefig(os.path.join(cmp_output_dir, f\"{cmp}_water_occupancy_van{van}.png\"))\n
 
# Extract vanilla replicate and window ID\nfeatures_df[\"van\"] = features_df[\"ensemble\"].str.extract(r\"van_(\\d)\")\nfeatures_df[\"id\"] = features_df[\"ensemble\"].str.replace(r\"_van_\\d+\", \"\", regex=True)\n\n# Compute average occupancy\navg_df = features_df.groupby([\"id\", \"van\"], as_index=False)[\"occupancy\"].mean()\navg_df.to_csv(os.path.join(cmp_output_dir, \"avg_water_occupancy.csv\"), index=False)\n\n# Plot average occupancy\nplt.figure(figsize=(12, 8))\nsns.lineplot(data=avg_df, x=\"id\", y=\"occupancy\", hue=\"van\", palette=\"tab10\")\nplt.title(f\"Average Water Occupancy Across Windows for {cmp}\", fontsize=16, fontweight=\"bold\")\nplt.xlabel(\"Window ID\", fontsize=14)\nplt.ylabel(\"Average Number of Waters\", fontsize=14)\nplt.legend(title=\"Vanilla Repeat\", title_fontsize=12, fontsize=10, loc=\"upper right\")\nplt.xticks(rotation=45, fontsize=10)\nplt.yticks(fontsize=10)\nplt.tight_layout()\nplt.savefig(os.path.join(cmp_output_dir, f\"{cmp}_water_occupancy_across_windows.png\"), dpi=300)\nplt.close()\n
","path":["How-to-guides","Binding Pocket Resolvation Analysis"],"tags":[]},{"location":"how-to-guides/guide-2/","level":1,"title":"Side Chain Torsion Analysis","text":"
Put the example analysis in the how to guide section here
","path":["How-to-guides","Side Chain Torsion Analysis"],"tags":[]},{"location":"tutorials/tutorial-1/","level":1,"title":"Featurizing ABFEs","text":"
This tutorial is designed to the user to introduce fepa by analyzing MD data from a set of simulations of the ligand 42922 (from the Deflorian et al set 1 datset) bound to the Orexin 2 Receptor. Specifically, we will compare the holo, apo and ABFE trajectories. Through this comparison, we will demonstrate that the protein’s binding pocket adopts distinct conformations between the final lambda windows (where the ligand is fully annihilated) of the ABFE simulation and the long apo simulation of the receptor.
 
This tutorial is followed by Tutorial 2 where we will use FEPA to set up REUS simulations to estimate the free energy of this conformational change from the holo like  to apo like states. This tutorial uses MDAnalysis, MDtraj and the PENSA package for analysis, GROMACS for simulation and Plumed for enhanced sampling. It is assumed that the user is familiar with setting up and analyzing MD simulations run with GROMACS and Plumed. 
","path":["Tutorials","Featurizing ABFEs"],"tags":[]},{"location":"tutorials/tutorial-1/#loading-the-config-file","level":2,"title":"Loading the config file","text":"
ABFE simulations produce multiple MD trajectories over different lamba windows. In our case, we have 44 lambda windows (11 Coulomb, 12 Restrains and 22 Van-der-Walls) plus the holo and the apo simulations that we have to analyse. To make things easier when parsing paths to the topology and coordinate files of these trajectories, we use a config file. The config file is json-formatted and contains all the information necessary to read the simulation trajectories. Here is a sample config file that I use for this tutorial:
 
{\n    \"base_path\": \"deflorian_set_1_j13_v1\",\n    \"abfe_window_path_template\": \"deflorian_set_1_j13_v1/OX2_{CMP_NAME}/abfe_van{REP_NO}_hrex_r{ABFE_REP_NO}/complex/{LEG_WINDOW}/{STAGE}\",\n    \"vanilla_path_template\": \"deflorian_set_1_j13_v1/OX2_{CMP_NAME}/vanilla_rep_{REP_NO}\",\n    \"apo_path_template\": \"deflorian_set_1_j13_v1/apo_OX2_r{REP_NO}\",\n    \"compounds\": [\n        \"42922\",\n    ],\n    \"pocket_residues_string\": \"12  54  57  58  59  60  61  62  63  64  65  70  71  78  81  82  83  85  86  89 138 142 160 161 162 163 175 178 179 182 183 232 235 236 239 240 242 243 261 265 268 269\"\n}\n
 
The pocket_residues_string variable is a string that stores the residue ids of all the proteins residue that have any atom within 6 A of the ligand. The JSON file is then loaded into a dictionary with the function load_config from fepa.utils.file_utils
","path":["Tutorials","Featurizing ABFEs"],"tags":[]},{"location":"tutorials/tutorial-1/#loading-md-trajectories","level":2,"title":"Loading MD trajectories","text":"
Now that we have the template paths to all the simulations: Apo EQ, Holo EQ and ABFE, we can use the function load_abfe_paths_for_compound from fepa.utils.path_utils to generate a path_dict dictionary that contains all the MD file paths for a single compound as follows:
 
cmp = config['compounds'][0]\npath_dict = load_abfe_paths_for_compound(\n            config,\n            cmp,\n            van_list=[3], # Loading only vanilla simulation 3\n            leg_window_list=[f\"coul.{i:02d}\" for i in range(0, 11,2)] + [f'vdw.{i:02d}' for i in range(0, 21,2)] + [f'rest.{i:02d}' for i in range(0, 12,2)], # Loading only every other lambda window\n            bp_selection_string=\"name CA and resid \" + config[\"pocket_residues_string\"],\n            apo=True,\n        )\n
 
The load_abfe_paths_for_compound function has various arguments that allow the user control over what simulation paths are loaded. For more information, please refer to the API. 
 
Now that we have the path_dict, it is time for us to load the trajectories themselves. We will be using the EnsembleHandler class from fepa.core.ensemble_handler to do this. EnsembleHandler is a neat way of storing and manipulating the trajectories from multiple ensembles with some built in functions for sanity checks. Internally EnsembleHandler stores the trajectories as dictionary of MDA universes.  
 
TODO: Need to remove ensemblehandler and just use universe dict
 
# Load trajectories\nensemble_handler = EnsembleHandler(path_dict)\n# Make universes\nensemble_handler.make_universes()\nlogging.info(\"Making universes for compound %s...\", cmp)\n# Check for BP residue consistency across all the trajectories\nlogging.info(\"Checking residue consistency for compound %s...\", cmp)\n
","path":["Tutorials","Featurizing ABFEs"],"tags":[]},{"location":"tutorials/tutorial-1/#featurizing-md-trajectories","level":2,"title":"Featurizing MD trajectories","text":"
We will now be featurizing these trajectories by their pairwise C-alpha distances in the binding pocket using the function SelfDistanceFeaturizer from fepa.core.featurizers. SelfDistanceFeaturizer computes and stored all possible pairs of distances between the C-alpha atoms of the binding pocket residues. The class also has save_features and load_features functions that help save the features as a csv file to make sure the time consuming featurization step need not be repeated every run.
 
TODO: Make sure the binding pocket selection is consistent
 
# Make a folder for the analysis output\ncmp_run_dir = f'analysis/{cmp}/'\nif !(os.path.exists(cmp_run_dir)):\n    os.mkdir(cmp_run_dir)\n\n# Featurize and save features\nfeaturizer = SelfDistanceFeaturizer(ensemble_handler)\nfeaturizer.featurize()\nfeaturizer.save_features(input_dir=cmp_existing_run_dir)\n
","path":["Tutorials","Featurizing ABFEs"],"tags":[]},{"location":"tutorials/tutorial-1/#visualizing-ensembles","level":2,"title":"Visualizing ensembles","text":"
Saving the features in a csv format gives us the flexibility to analyse it as required. For the purpose of this tutorial, we will be looking at how our features capture the difference between the apo, the holo and the abfe ensembles. To do this we reduce the dimensions of the features data using the PCADimReducer class from fepa.core.dim_reducers. FEPA also supports other dimensionality reduction techniqeus like UMAP and tSNE. In fact UMAP is better able to resolve the differences in binding pocket configurations between different ensembles. We will be doing PCA here as it also doubles up as a nice CV to bias when performing umbrella sampling later.
 
# Dimensionality Reduction\nlogging.info(\"Performing dimensionality reduction for compound %s...\", cmp)\ndimreducer = PCADimReducer(featurizer.get_feature_df(), n_components=8)\ndimreducer.reduce_dimensions()\ndimreducer.calculate_projections()\ndimreducer.save_projection_df(\n    save_path=os.path.join(cmp_output_dir, \"pca_projection_df.csv\")\n)\n
 
First we plot the eigen values of all the PCs to understand what percentage of variance is captured by the first few PCs:
 
logging.info(\"Plotting PCA eigenvalues for compound %s...\", cmp)\nplot_eigenvalues(\n    pca_object=dimreducer.get_pca(),\n    n_components=8,\n    save_path=os.path.join(cmp_output_dir, \"eigenvalues.png\"),\n)\n
 
 
Figure 1: PCA eigenvalues for the binding-pocket Cα self-distance features.
 
Here, most of the variance is captured by PC1. However, this may not always be the case. Therefore, for this tutorial, we will use PC1 and PC2 together, as they collectively capture the majority of the variance. To improve visualization, we define a new column in projection_df called simtype, which groups the ensembles into broader categories: abfe_coul, abfe_vdw, abfe_rest, holo_equil, apo, and nvt.
 
# Further labelling the ensembles\nprojection_df = dimreducer.get_pca_projection_df()\n# Add another column called simtype based on ensemble\ndef get_simtype(ensemble_name: str) -> str:\n    if \"van\" in ensemble_name:\n        if 'nvt' in ensemble_name:\n            return \"nvt\"\n        if \"coul\" in ensemble_name:\n            return \"abfe_coul\"\n        elif \"vdw\" in ensemble_name:\n            return \"abfe_vdw\"\n        elif \"rest\" in ensemble_name:\n            return \"abfe_rest\"\n        else:\n            return \"holo_equil\"\n    elif \"apo\" in ensemble_name:\n        return \"apo\"\n    else:\n        return \"other\"\nprojection_df[\"simtype\"] = projection_df[\"ensemble\"].apply(get_simtype)\n
 
We can visualize the dimensionality-reduced data using the DimRedVisualizer class. In the example below, we plot the first two principal components, colored by simulation and time. We also highlight NVT, as it represents the initial crystal structure after energy minimization.
 
logging.info(\"Visualizing compound %s...\", cmp)\ndimred_visualizer = DimRedVisualizer(projection_df=projection_df, data_name=\"PCA\")\ndimred_visualizer.plot_dimred_sims(\n    save_path=os.path.join(cmp_run_dir, \"pca_components_ensemble_noapo.png\"),\n    highlights=[f\"{cmp}_nvt\"],\n)\ndimred_visualizer.plot_dimred_time(\n    save_path=os.path.join(cmp_run_dir, \"pca_components_time_noapo.png\")\n)\ndimred_visualizer.plot_dimred_sims(\n    column=\"simtype\",\n    save_path=os.path.join(cmp_run_dir, \"pca_components_simtype.png\"),\n)\n
 
 Figure 2: Simulation frames in PC-space (PC1 vs PC2), colored by simtype
","path":["Tutorials","Featurizing ABFEs"],"tags":[]},{"location":"tutorials/tutorial-1/#clustering-ensembles","level":2,"title":"Clustering ensembles","text":"
In Figure 2, frames from the ABFE simulations closely resemble those from the holo equilibrium simulation, forming two distinct clusters: apo-like and holo-like. We will be clustering them using cluster_pca. After clustering, we can estimtate the data point that is closest to the center with the function make_ensemble_center_df which returns a DataFrame containing the details of the frame closest to the centroid in PC space. 
 
# Cluster the projection df\npca_projection_df_clustered = cluster_pca(\n    projection_df, n_clusters=3, n_components=3\n)\n# Ensemble center df\nensemble_center_df = make_ensemble_center_df(\n    pca_projection_df_clustered, key=\"cluster\"\n)\n
 
We then visualize the clustered data with ensemble centers as follows:
 
# Visualization of clustered data\ndimred_visualizer_clustered = DimRedVisualizer(\n    projection_df=pca_projection_df_clustered, data_name=\"PCA\"\n)\ndimred_visualizer_clustered.plot_dimred_cluster(\n    save_path=os.path.join(\n        cmp_run_dir, \"subset_pca_components_clusters_w_center.png\"\n    ),\n    centroid_df=ensemble_center_df,\n    cluster_column=\"cluster\",\n)\n
 
 
Figure 3: Clustered PCA data with ensemble centers marked by an X.
","path":["Tutorials","Featurizing ABFEs"],"tags":[]},{"location":"tutorials/tutorial-1/#visualizing-differences-across-clustered-ensembles","level":2,"title":"Visualizing differences across clustered ensembles","text":"
With two distinct clusters identified, we can compare them across different features by computing the Jensen–Shannon (JS) entropy for each feature. The code below plots feature-wise histograms for cluster 0 and cluster 1, annotated with their JS entropy values:
 
# Copy the cluster labels from pca_projection_df_clustered to feature_df\nfeature_df_w_clusters = featurizer.get_feature_df().copy()\nfeature_df_w_clusters[\"cluster\"] = pca_projection_df_clustered[\"cluster\"].values\n\n\n# Compute feature-level histograms and JS entropy between clusters 0 and 1\nhistograms = compute_histograms(\n    feature_df_w_clusters,\n    \"cluster\",\n    0,\n    1,\n    num_bins=50,\n    feature_column_keyword=\"DIST\",\n)\n\nrel_ent_dict = compute_relative_entropy(\n    feature_df_w_clusters,\n    ensemble1=0,\n    ensemble2=1,\n    num_bins=50,\n    key=\"cluster\",\n    feature_column_keyword=\"DIST\",\n)\n\nplot_jsd_histograms(\n    histograms,\n    rel_ent_dict,\n    top_n=16,\n    save_path=os.path.join(cmp_run_dir, \"jsd_cluster0_vs_cluster1.png\"),\n)\n
 
 Figure 4: Histograms of the top 16 features showing the highest Jensen–Shannon divergence between clusters 0 and 1.
 
We then save the ensemble center frames as GRO files to compare structures in molecular visualization tools such as PyMOL:
 
# Save ensemble center frames as GRO files\nfor _, row in ensemble_center_df.iterrows():\n    center_ensemble = row[\"ensemble\"]\n    center_timestep = row[\"timestep\"]\n    print(f\"Ensemble: {center_ensemble}, Timestep: {center_timestep}\")\n\n    # Load trajectories\n    ensemble_handler.make_universes()\n\n    # Export the frame corresponding to the ensemble center\n    ensemble_handler.dump_frames(\n        ensemble=center_ensemble,\n        timestep=center_timestep,\n        save_path=f\"subset_cluster_{int(row['cluster'])}_center.gro\",\n    )\n
","path":["Tutorials","Featurizing ABFEs"],"tags":[]},{"location":"tutorials/tutorial-2/","level":1,"title":"Setting up REUS","text":"
In Tutorial 2, we build on the concepts learned in Tutorial 1, where the apo and holo states were shown to form distinct conformational clusters. Using the data provided in the [template], we will set up Replica Exchange Umbrella Sampling (REUS) simulations to estimate the free energy of the conformational transition from holo to apo, using the principal component (PC) 1 as the collective variable (CV). We will achieve this by morphing the protein from one state to another, equilibrating the system, and then running REUS to sample the transition.
","path":["Tutorials","Setting up REUS"],"tags":[]},{"location":"tutorials/tutorial-2/#preparing-the-plumeddat-file","level":2,"title":"Preparing the plumed.dat file","text":"
The first step in performing umbrella sampling with GROMACS and PLUMED is to create a plumed.dat file that defines the collective variable (CV), and the positions of the umbrella sampling restraints. To do this, we load the feature CSV files, initialize a DimReducer object, and extract the corresponding PCA object to define the CV based on pairwise Cα distances.
 
# Load Features df and reduce dimensions\nfeature_df = pd.read_csv(os.path.join(\"md_data\", \"top_features_apo_vdw.20.csv\"))\ndimreducer = PCADimReducer(feature_df, n_components=8)\ndimreducer.reduce_dimensions()\ndimreducer.calculate_projections()\npca_projection_df = dimreducer.get_pca_projection_df()\n# Use PCA from dimreducer to write plumed file\ntop_features = feature_df.filter(regex=\"DIST\", axis=1).columns\nwrite_plumed_file(\n    sdf_names=top_features,\n    top_features_pca=dimreducer.get_pca(),\n    save_path=\"plumed.dat\",\n    molinfo_structure=\"../reference.pdb\",  # fix molinfo here\n)\n
 
The PLUMED input file should look like this:
 
MOLINFO STRUCTURE=../reference.pdb\nd1: DISTANCE ATOMS=@CA-236,@CA-243\nd2: DISTANCE ATOMS=@CA-175,@CA-243\nd3: DISTANCE ATOMS=@CA-138,@CA-243\n.\n.\n.\nd200: DISTANCE ATOMS=@CA-61,@CA-70\n# Create the dot product\ndot: COMBINE ARG=d1,d2,d3...d200 COEFFICIENTS=0.095,0.152,0.133...0.057,0.054,0.032 PERIODIC=NO\nCV: MATHEVAL ARG=dot FUNC=10*x PERIODIC=NO\nPRINT ARG=CV FILE=COLVAR STRIDE=1\n# Put position of restraints here for each window\nrestraint: RESTRAINT ARG=CV AT=@replicas:$RESTRAINT_ARRAY KAPPA=$KAPPA\nPRINT ARG=restraint.* FILE=restr\n
 Note that $RESTRAINT_ARRAY is a placeholder for the harmonic restraint positions in CV space. When defining these positions, ensure they align with the direction of the structural morph. By convention, we transition from vdw.20 to apo structures. To assign restraints correctly, identify the ensemble with the lowest PC value and order the restraints from min to max (or vice versa) accordingly.
 
# Pair of ensembles to compare\npair = (\"vdw.20\", \"apo\")\n# Get mean PC1 for the two ensembles\nmean_pc1_ensemble1 = pca_projection_df[pca_projection_df[\"state\"] == pair[0]][\n    \"PC1\"\n].mean()\nmean_pc1_ensemble2 = pca_projection_df[pca_projection_df[\"state\"] == pair[1]][\n    \"PC1\"\n].mean()\nlogging.info(\n    f\"Mean PC1 for {pair[0]}: {mean_pc1_ensemble1}, Mean PC1 for {pair[1]}: {mean_pc1_ensemble2}\"\n)\n# Get the restraint array based on the two ensembles\nPC1_min = pca_projection_df[\"PC1\"].min()\nPC1_max = pca_projection_df[\"PC1\"].max()\n# If min PC1 is from ensemble 1, then the restraint array should be from PC1_min to PC1_max\nif mean_pc1_ensemble1 < mean_pc1_ensemble2:\n    restraint_array = np.linspace(PC1_min, PC1_max, 24)\nelif mean_pc1_ensemble1 > mean_pc1_ensemble2:\n    restraint_array = np.linspace(PC1_max, PC1_min, 24)\n
 
Once the restraint array is prepared, we write it to the PLUMED file using the write_plumed_restraints function:
 
# Write the restraint array to the plumed file\nwrite_plumed_restraints(\n    plumed_file=\"plumed.dat\",\n    restraint_centers=restraint_array,\n    kappa=5,\n)\n
 
The restrain line on plumed.dat file must now be an array of distance restraints that looks like this:
 
restraint: RESTRAINT ARG=CV AT=@replicas:116.856,119.047,121.239...167.262 KAPPA=5\n
","path":["Tutorials","Setting up REUS"],"tags":[]},{"location":"tutorials/tutorial-2/#morphing-with-memento","level":2,"title":"Morphing with memento","text":"
For this tutorial, the GRO files representing the vdw.20 and apo states are provided in the data folder, prepared in the same way as in Tutorial 1. The plumed.dat file has also been prepared. The next step for REUS is to generate intermediate protein conformations using Memento (JCTC, 2023) . FEPA provides a class for easy access to Memento.
 
Using FEPA’s memento_workflow class, we can perform the protein morphing. A Memento directory is needed to store morphs for each entry in initial_target_gro_dict, and an apo_template path must be provided, containing the topology and MDP files required for equilibrium simulations.
 
# Declaring variables:\nmemento_dir = \"/biggin/b211/reub0138/Projects/orexin/deflorian_set_1_j13_v1_memento\"\n# Defining pairs\npair = (\"all_vdw20\", \"apo_3\")\nall_vdw20_gro_file = \"/biggin/b211/reub0138/Projects/orexin/deflorian_set_1_j13_v1/analysis/a12_dimreduce_vdw.20/subset_cluster_1_center.gro\"\napo_3_gro_file = \"/biggin/b211/reub0138/Projects/orexin/deflorian_set_1_j13_v1/analysis/a12_dimreduce_vdw.20/subset_cluster_2_center.gro\"\ntemplate_path = \"/biggin/b211/reub0138/Projects/orexin/deflorian_set_1_j13_v1_memento/apo_template\"  # Make sure no water and no ions in the template topology file\n# Setup Memento folders\nmemento_flow = memento_workflow(\n    memento_dir=memento_dir,\n    initial_gro=all_vdw20_gro_file,\n    target_gro=apo_3_gro_file,\n    initial_name=pair[0],\n    target_name=pair[1],\n    template_path=template_path,\n    run_name=\"memento_run_v1\",\n    n_residues=296,\n)\n
 
We can use the workflow functions to perform each step. First, prepare_memento sets up the files required by Memento. Then, run_memento executes the morphing. run_memento requires the protonation states of all histidines in the format expected by GROMACS pdb2gmx (0 for HID, 1 for HIE, 2 for HIP, 3 for HIS1). Additionally, the indices of CYX residues must be provided, as Memento cannot process them automatically.
 
# Preparing memento input\nmemento_flow.prepare_memento()\n# Running memento\nmemento_flow.run_memento(\n    template_path=template_path,\n    last_run=\"memento_run_v1\",\n    protonation_states=[1, 1, 1, 2, 1, 2],\n    cyx_residue_indices = []\n)\n
","path":["Tutorials","Setting up REUS"],"tags":[]},{"location":"tutorials/tutorial-2/#running-equilibration","level":2,"title":"Running Equilibration","text":"
This step may take some time. Once complete, we can set up equilibration simulations with Memento by providing a job script path, which can be modified to suit the specific HPC system.
 
# Running analysis\nmemento_flow.prepare_equil_simulations(\n    job_script_template=\"/biggin/b211/reub0138/Projects/orexin/lenselink_a2a_memento_v1/job_vanilla_ranv_equil_arr_template.sh\"\n)\n
 
After running prepare_equil_simulations, individual boxes for each morph are created. We then simulate each box to relax the side chains by submitting the job to our HPC. Once complete, we can proceed to the REUS simulations.
","path":["Tutorials","Setting up REUS"],"tags":[]},{"location":"tutorials/tutorial-2/#preparing-reus","level":2,"title":"Preparing REUS","text":"
Once the equilibration simulations are complete, use FEPA’s reus_umbrella_sampling_workflow class to set up umbrella sampling in the same working directory. The number of windows can be adjusted via n_windows (24 is typically sufficient). If the residue IDs in the plumed.dat file does not match those in the MD GRO files, FEPA requires the correct offset to be set via the plumed_resid_offset parameter.
 
import logging\nfrom pathlib import Path\nfrom fepa.flows.reus_flows import (\n    reus_umbrella_sampling_workflow,\n)\n\nsim_path = \"wdir/all_vdw20_apo_3/memento_run_v1/wdir/boxes/sim0\"\nwdir_path = Path(sim_path).parents[1]\nplumed_path = \"plumed.dat\"\nsimname = \"all_vdw20_apo_3\"\ninitial_gro = \"md_data/cluster_1_center.gro\"\nfinal_gro = \"md_data/cluster_2_center.gro\"\nsubmission_script_template_arr = \"md_data/job_reus_template.sh\"\nprint(f\"wdir_path: {wdir_path}, plumed_path: {plumed_path}\")\numbrella_sampler = reus_umbrella_sampling_workflow(\n    wdir_path=wdir_path,\n    plumed_path=plumed_path,\n    submission_script_template_arr=submission_script_template_arr,\n    start=simname.split(\"_\")[0] + \"_\" + simname.split(\"_\")[1],\n    end=simname.split(\"_\")[2] + \"_\" + simname.split(\"_\")[3],\n    reus_folder_name=\"reus_v1\",\n    n_windows=24,\n    plumed_resid_offset=0,\n    initial_gro=initial_gro,\n    target_gro=final_gro,\n)\numbrella_sampler.setup_simulations(exist_ok=True)\n
","path":["Tutorials","Setting up REUS"],"tags":[]},{"location":"tutorials/tutorial-2/#analyzing-reus","level":2,"title":"Analyzing REUS","text":"
Once your REUS system is set up, you can analyze the results using the workflow:
 
umbrella_sampler.prepare_wham()\numbrella_sampler.run_wham()\numbrella_sampler.analyse_us_hist(range=(90, 180), colvar_prefix=\"COLVAR\")\numbrella_sampler.get_initial_final_CVs()\numbrella_sampler.plot_free_energies(\n    units=\"kcal\",\n)\n
 
This prepares and runs WHAM on the histograms and generates the free energy curves.
 
 
Figure 1: Free energy curves from the tutorial simulations. The CV values of apo and holo structures are marked with gray dotted lines. Different lines represent curves computed using varying proportions of data to assess convergence.
","path":["Tutorials","Setting up REUS"],"tags":[]}]}
\ No newline at end of file

From f3dd7a5b54bc607da78e5ab46a93330fc305ef0d Mon Sep 17 00:00:00 2001
From: nithishwer 
Date: Thu, 13 Nov 2025 11:54:50 +0000
Subject: [PATCH 21/40] fix: add python 3.10 to req

---
 requirements.yml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/requirements.yml b/requirements.yml
index 3ccee4e..1698643 100644
--- a/requirements.yml
+++ b/requirements.yml
@@ -4,6 +4,7 @@ channels:
   - conda-forge
   - omnia
 dependencies:
+  - python=3.10
   - plumed
   - modeller
   - mdaencore

From 45e55baa410d7f8357d1e572b16fb01fe3a661a0 Mon Sep 17 00:00:00 2001
From: nithishwer 
Date: Thu, 13 Nov 2025 12:07:51 +0000
Subject: [PATCH 22/40] fix: add fepa env activate with test

---
 .github/workflows/test.yml | 16 +++++++++++-----
 1 file changed, 11 insertions(+), 5 deletions(-)

diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index d124fc5..dddc0fb 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -14,11 +14,13 @@ jobs:
       - name: Checkout code
         uses: actions/checkout@v4
 
-      - name: Configure Conda channels
-        run: |
-          conda config --add channels salilab
-          conda config --add channels conda-forge
-          conda config --add channels defaults
+      - name: Cache Conda packages
+        uses: actions/cache@v4
+        with:
+          path: ~/.conda/pkgs
+          key: conda-${{ runner.os }}-${{ hashFiles('requirements.yml') }}
+          restore-keys: |
+            conda-${{ runner.os }}-
 
       - name: Set up Conda environment
         uses: conda-incubator/setup-miniconda@v3
@@ -29,9 +31,13 @@ jobs:
           use-mamba: true
 
       - name: Install package in editable mode
+        shell: bash -l {0}
         run: |
+          conda activate fepa_env
           pip install --no-cache-dir -e .
 
       - name: Run tests
+        shell: bash -l {0}
         run: |
+          conda activate fepa_env
           pytest -v --maxfail=1 --disable-warnings

From 21fa9b7621ee2e1bcf0132d88b3bd22cfd0e6b25 Mon Sep 17 00:00:00 2001
From: nithishwer 
Date: Thu, 13 Nov 2025 13:39:28 +0000
Subject: [PATCH 23/40] fix: using conda for test

---
 .github/workflows/test.yml | 10 ++--------
 tests/conftest.py          | 22 ++++++++++++----------
 2 files changed, 14 insertions(+), 18 deletions(-)

diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index dddc0fb..5bb5d3e 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -31,13 +31,7 @@ jobs:
           use-mamba: true
 
       - name: Install package in editable mode
-        shell: bash -l {0}
-        run: |
-          conda activate fepa_env
-          pip install --no-cache-dir -e .
+        run: conda run -n fepa_env pip install --no-cache-dir -e .
 
       - name: Run tests
-        shell: bash -l {0}
-        run: |
-          conda activate fepa_env
-          pytest -v --maxfail=1 --disable-warnings
+        run: conda run -n fepa_env pytest -v --maxfail=1 --disable-warnings
diff --git a/tests/conftest.py b/tests/conftest.py
index e359a73..295d4e6 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -1,11 +1,13 @@
 from pathlib import Path
 import pytest
+
 # Assuming this function loads the JSON file into a Python dictionary
-from fepa.utils.file_utils import load_config 
+from fepa.utils.file_utils import load_config
+
 
 def make_paths_absolute(config: dict, root: Path) -> dict:
     """
-    Converts relative path templates in the config to absolute paths 
+    Converts relative path templates in the config to absolute paths
     based on the project root.
     """
     # Keys in your config that contain relative paths
@@ -13,9 +15,9 @@ def make_paths_absolute(config: dict, root: Path) -> dict:
         "abfe_window_path_template",
         "vanilla_path_template",
         "vanilla_path_template_old",
-        "apo_path_template"
+        "apo_path_template",
     ]
-    
+
     # Iterate and update the paths
     for key in path_keys:
         if key in config:
@@ -23,24 +25,24 @@ def make_paths_absolute(config: dict, root: Path) -> dict:
             # We use Path() for safe, cross-platform path joining.
             relative_path = Path(config[key])
             config[key] = str(root / relative_path)
-            
+
     return config
 
 
 @pytest.fixture
 def test_env():
     """Provide the test environment with project root and loaded test configuration."""
-    
+
     # 1. Determine the project root (assuming tests/ is one level down)
-    # Path(__file__).resolve().parents[1] correctly points to the directory 
+    # Path(__file__).resolve().parents[1] correctly points to the directory
     # above 'tests', which is your project root.
     root = Path(__file__).resolve().parents[1]
-    
+
     # 2. Load the configuration file
     config_path = root / "tests" / "test_config" / "config.json"
     config = load_config(config_path)
 
     # 3. Convert relative paths inside the loaded config to absolute paths
     config = make_paths_absolute(config, root)
-    
-    return {"root": root, "config": config}
\ No newline at end of file
+
+    return {"root": root, "config": config}

From df8f9ed124380df479c76eff8a760e30d8a875d4 Mon Sep 17 00:00:00 2001
From: nithishwer 
Date: Thu, 13 Nov 2025 13:52:43 +0000
Subject: [PATCH 24/40] fix: add pytest to requirements.yml

---
 requirements.yml | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/requirements.yml b/requirements.yml
index 1698643..3fdb99a 100644
--- a/requirements.yml
+++ b/requirements.yml
@@ -13,4 +13,6 @@ dependencies:
   - cython
   - openmm
   - openmmforcefields
-  - numba
\ No newline at end of file
+  - numba
+  - pytest
+  
\ No newline at end of file

From 960ab005a232d403ebfb8a35357674c0ebe1ba30 Mon Sep 17 00:00:00 2001
From: nithishwer 
Date: Thu, 13 Nov 2025 13:58:13 +0000
Subject: [PATCH 25/40] fix: add umap to requirements,yml

---
 requirements.yml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/requirements.yml b/requirements.yml
index 3fdb99a..d6b7bbd 100644
--- a/requirements.yml
+++ b/requirements.yml
@@ -15,4 +15,5 @@ dependencies:
   - openmmforcefields
   - numba
   - pytest
+  - umap-learn
   
\ No newline at end of file

From 2142cca1aa4748bf1a447b083c571b74de5ac756 Mon Sep 17 00:00:00 2001
From: nithishwer 
Date: Thu, 13 Nov 2025 14:09:18 +0000
Subject: [PATCH 26/40] fix: add pymemento to requirement

---
 requirements.yml | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/requirements.yml b/requirements.yml
index d6b7bbd..f6c1624 100644
--- a/requirements.yml
+++ b/requirements.yml
@@ -16,4 +16,5 @@ dependencies:
   - numba
   - pytest
   - umap-learn
-  
\ No newline at end of file
+  - pip:
+    - PyMEMENTO
\ No newline at end of file

From 65da581cd38f356d35b443417f057aa257020d0d Mon Sep 17 00:00:00 2001
From: nithishwer 
Date: Thu, 13 Nov 2025 14:36:06 +0000
Subject: [PATCH 27/40] fix: fepa test import

---
 .../a1_test/1_binding_pocket_analysis_test.py | 13 +---------
 ...rs_across_repeats_v2_analysisclass_test.py | 19 ++++++++------
 tests/a5_test/1_sc_torsions_featurize_test.py | 25 ++++++++++---------
 3 files changed, 25 insertions(+), 32 deletions(-)

diff --git a/tests/a1_test/1_binding_pocket_analysis_test.py b/tests/a1_test/1_binding_pocket_analysis_test.py
index 3e891b7..8982ef5 100644
--- a/tests/a1_test/1_binding_pocket_analysis_test.py
+++ b/tests/a1_test/1_binding_pocket_analysis_test.py
@@ -24,15 +24,10 @@
 warnings.filterwarnings("ignore", category=DeprecationWarning, module=".*importlib.*")
 
 import os
-from pathlib import Path
 import builtins
 from typing import Literal
 
 import pytest
-import shutil
-import numpy as np
-import pandas as pd
-import pandas.testing as pdt
 
 # Default numeric tolerances for floating-point comparisons
 DECIMALS = 6
@@ -43,20 +38,14 @@
 builtins.Literal = Literal
 
 # --- FEPA imports ---
-from fepa.utils.file_utils import load_config
 from fepa.core.ensemble_handler import EnsembleHandler
 from fepa.utils.path_utils import load_abfe_paths_for_compound
 from fepa.flows import binding_pocket_analysis_workflow
 
 # --- Test utilities and fixtures ---
-from fepa.tests.utils import (
-    round_numeric,
-    sort_by,
-    align_pc_signs,
-    first_existing,
+from tests.utils import (
     check_csv_equality,
 )
-from fepa.tests.conftest import test_env as test_env_fixture
 
 
 @pytest.mark.integration
diff --git a/tests/a2_test/1_plot_bp_waters_across_repeats_v2_analysisclass_test.py b/tests/a2_test/1_plot_bp_waters_across_repeats_v2_analysisclass_test.py
index 99e9395..21976db 100644
--- a/tests/a2_test/1_plot_bp_waters_across_repeats_v2_analysisclass_test.py
+++ b/tests/a2_test/1_plot_bp_waters_across_repeats_v2_analysisclass_test.py
@@ -12,17 +12,12 @@
 warnings.filterwarnings("ignore", category=DeprecationWarning, module="MDAnalysis.*")
 warnings.filterwarnings("ignore", category=DeprecationWarning, module=".*importlib.*")
 
-from pathlib import Path
-import pandas as pd
-import pandas.testing as pdt
 import pytest
 
-from fepa.utils.file_utils import load_config
 from fepa.core.ensemble_handler import EnsembleHandler
 from fepa.utils.path_utils import load_abfe_paths_for_compound
 from fepa.core.featurizers import BPWaterFeaturizer
-from fepa.tests.utils import round_numeric, align_on_common_keys, check_csv_equality
-from fepa.tests.conftest import test_env
+from tests.utils import check_csv_equality
 
 DECIMALS = 6
 RTOL = 1e-6
@@ -34,6 +29,7 @@
     "173 176 177 180 217 218 221 225 235 238 239 240 241 242 243 244 245 246 247"
 )
 
+
 @pytest.mark.slow
 @pytest.mark.integration
 def test_bp_waters_minimal_against_expected(tmp_path, test_env):
@@ -65,8 +61,15 @@ def test_bp_waters_minimal_against_expected(tmp_path, test_env):
 
     # --- Comparison ---
     expected_csv = (
-        repo_root / "tests" / "test_data" / "2_expected" / cmp_name / "WaterOccupancy_features.csv"
+        repo_root
+        / "tests"
+        / "test_data"
+        / "2_expected"
+        / cmp_name
+        / "WaterOccupancy_features.csv"
     )
     assert expected_csv.exists(), f"Expected CSV not found: {expected_csv}"
 
-    check_csv_equality(act_csv, expected_csv, label="Water occupancy CSV", rtol=RTOL, atol=ATOL)
+    check_csv_equality(
+        act_csv, expected_csv, label="Water occupancy CSV", rtol=RTOL, atol=ATOL
+    )
diff --git a/tests/a5_test/1_sc_torsions_featurize_test.py b/tests/a5_test/1_sc_torsions_featurize_test.py
index 563a3ef..f537876 100644
--- a/tests/a5_test/1_sc_torsions_featurize_test.py
+++ b/tests/a5_test/1_sc_torsions_featurize_test.py
@@ -8,6 +8,7 @@
 import warnings
 from Bio import BiopythonDeprecationWarning
 import sys
+
 print("Python executable:", sys.executable)
 print("Python version:", sys.version.splitlines()[0])
 print("sys.path:")
@@ -18,26 +19,19 @@
 warnings.filterwarnings("ignore", category=DeprecationWarning, module="MDAnalysis.*")
 warnings.filterwarnings("ignore", category=DeprecationWarning, module=".*importlib.*")
 
-from pathlib import Path
 import pytest
 import shutil
-import pandas as pd
-import pandas.testing as pdt
 
-from fepa.utils.file_utils import load_config
 from fepa.core.ensemble_handler import EnsembleHandler
 from fepa.utils.path_utils import load_abfe_paths_for_compound
 from fepa.core.featurizers import SideChainTorsionsFeaturizer
-from fepa.tests.utils import check_csv_equality
-from fepa.tests.conftest import test_env
+from tests.utils import check_csv_equality
 import sys
 
 RTOL = 1e-6
 ATOL = 1e-8
 
-BP_SELECTION_STRING = (
-    "name CA and resid 54 55 56 57" 
-) # Doesnt matter what residues we pick here, as sidechain torsions dont use this selection
+BP_SELECTION_STRING = "name CA and resid 54 55 56 57"  # Doesnt matter what residues we pick here, as sidechain torsions dont use this selection
 
 
 @pytest.mark.slow
@@ -70,12 +64,19 @@ def test_sidechain_torsions_minimal_against_est(tmp_path, test_env):
     assert act_csv.exists(), f"Expected feature file not found: {act_csv}"
 
     # Copy actual CSV to a known location for easier debugging if needed
-    shutil.copy(act_csv, repo_root / "tests" / "test_data" )
+    shutil.copy(act_csv, repo_root / "tests" / "test_data")
 
     # --- Comparison ---
     expected_csv = (
-        repo_root / "tests" / "test_data" / "5_expected" / cmp_name / "SideChainTorsions_features.csv"
+        repo_root
+        / "tests"
+        / "test_data"
+        / "5_expected"
+        / cmp_name
+        / "SideChainTorsions_features.csv"
     )
     assert expected_csv.exists(), f"Golden CSV not found: {expected_csv}"
 
-    check_csv_equality(act_csv, expected_csv, label="Side-chain torsions CSV", rtol=RTOL, atol=ATOL)
+    check_csv_equality(
+        act_csv, expected_csv, label="Side-chain torsions CSV", rtol=RTOL, atol=ATOL
+    )

From a20b595af01b8f41e54a79a448e9ba39ca1f30b2 Mon Sep 17 00:00:00 2001
From: nithishwer 
Date: Thu, 13 Nov 2025 14:51:35 +0000
Subject: [PATCH 28/40] fix: change util import in tests

---
 tests/a1_test/1_binding_pocket_analysis_test.py                 | 2 +-
 .../1_plot_bp_waters_across_repeats_v2_analysisclass_test.py    | 2 +-
 tests/a5_test/1_sc_torsions_featurize_test.py                   | 2 +-
 3 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/tests/a1_test/1_binding_pocket_analysis_test.py b/tests/a1_test/1_binding_pocket_analysis_test.py
index 8982ef5..2bbd157 100644
--- a/tests/a1_test/1_binding_pocket_analysis_test.py
+++ b/tests/a1_test/1_binding_pocket_analysis_test.py
@@ -43,7 +43,7 @@
 from fepa.flows import binding_pocket_analysis_workflow
 
 # --- Test utilities and fixtures ---
-from tests.utils import (
+from ..tests.utils import (
     check_csv_equality,
 )
 
diff --git a/tests/a2_test/1_plot_bp_waters_across_repeats_v2_analysisclass_test.py b/tests/a2_test/1_plot_bp_waters_across_repeats_v2_analysisclass_test.py
index 21976db..bf547ff 100644
--- a/tests/a2_test/1_plot_bp_waters_across_repeats_v2_analysisclass_test.py
+++ b/tests/a2_test/1_plot_bp_waters_across_repeats_v2_analysisclass_test.py
@@ -17,7 +17,7 @@
 from fepa.core.ensemble_handler import EnsembleHandler
 from fepa.utils.path_utils import load_abfe_paths_for_compound
 from fepa.core.featurizers import BPWaterFeaturizer
-from tests.utils import check_csv_equality
+from ..tests.utils import check_csv_equality
 
 DECIMALS = 6
 RTOL = 1e-6
diff --git a/tests/a5_test/1_sc_torsions_featurize_test.py b/tests/a5_test/1_sc_torsions_featurize_test.py
index f537876..8c9196d 100644
--- a/tests/a5_test/1_sc_torsions_featurize_test.py
+++ b/tests/a5_test/1_sc_torsions_featurize_test.py
@@ -25,7 +25,7 @@
 from fepa.core.ensemble_handler import EnsembleHandler
 from fepa.utils.path_utils import load_abfe_paths_for_compound
 from fepa.core.featurizers import SideChainTorsionsFeaturizer
-from tests.utils import check_csv_equality
+from ..tests.utils import check_csv_equality
 import sys
 
 RTOL = 1e-6

From 0e91b4f903a83f7067ecdcfe531332ea642cefc6 Mon Sep 17 00:00:00 2001
From: nithishwer 
Date: Thu, 13 Nov 2025 15:28:50 +0000
Subject: [PATCH 29/40] feat: import  utils from tests

---
 .../1_plot_bp_waters_across_repeats_v2_analysisclass_test.py    | 2 +-
 tests/a5_test/1_sc_torsions_featurize_test.py                   | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/tests/a2_test/1_plot_bp_waters_across_repeats_v2_analysisclass_test.py b/tests/a2_test/1_plot_bp_waters_across_repeats_v2_analysisclass_test.py
index bf547ff..21976db 100644
--- a/tests/a2_test/1_plot_bp_waters_across_repeats_v2_analysisclass_test.py
+++ b/tests/a2_test/1_plot_bp_waters_across_repeats_v2_analysisclass_test.py
@@ -17,7 +17,7 @@
 from fepa.core.ensemble_handler import EnsembleHandler
 from fepa.utils.path_utils import load_abfe_paths_for_compound
 from fepa.core.featurizers import BPWaterFeaturizer
-from ..tests.utils import check_csv_equality
+from tests.utils import check_csv_equality
 
 DECIMALS = 6
 RTOL = 1e-6
diff --git a/tests/a5_test/1_sc_torsions_featurize_test.py b/tests/a5_test/1_sc_torsions_featurize_test.py
index 8c9196d..f537876 100644
--- a/tests/a5_test/1_sc_torsions_featurize_test.py
+++ b/tests/a5_test/1_sc_torsions_featurize_test.py
@@ -25,7 +25,7 @@
 from fepa.core.ensemble_handler import EnsembleHandler
 from fepa.utils.path_utils import load_abfe_paths_for_compound
 from fepa.core.featurizers import SideChainTorsionsFeaturizer
-from ..tests.utils import check_csv_equality
+from tests.utils import check_csv_equality
 import sys
 
 RTOL = 1e-6

From 9bd6fb256c5d9732f85f8f7287eb2c0d7ac0755d Mon Sep 17 00:00:00 2001
From: nithishwer 
Date: Thu, 13 Nov 2025 15:34:20 +0000
Subject: [PATCH 30/40] feat: import utils from tests a1

---
 tests/a1_test/1_binding_pocket_analysis_test.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/a1_test/1_binding_pocket_analysis_test.py b/tests/a1_test/1_binding_pocket_analysis_test.py
index 2bbd157..8982ef5 100644
--- a/tests/a1_test/1_binding_pocket_analysis_test.py
+++ b/tests/a1_test/1_binding_pocket_analysis_test.py
@@ -43,7 +43,7 @@
 from fepa.flows import binding_pocket_analysis_workflow
 
 # --- Test utilities and fixtures ---
-from ..tests.utils import (
+from tests.utils import (
     check_csv_equality,
 )
 

From ea04ee1fbb5473d57fe1fb8bc84ede95b3d041e9 Mon Sep 17 00:00:00 2001
From: nithishwer 
Date: Thu, 13 Nov 2025 15:59:09 +0000
Subject: [PATCH 31/40] fix: add init.py to test subfolders

---
 tests/a1_test/__init__.py | 0
 tests/a2_test/__init__.py | 0
 tests/a5_test/__init__.py | 0
 tests/a6_test/__init__.py | 0
 4 files changed, 0 insertions(+), 0 deletions(-)
 create mode 100644 tests/a1_test/__init__.py
 create mode 100644 tests/a2_test/__init__.py
 create mode 100644 tests/a5_test/__init__.py
 create mode 100644 tests/a6_test/__init__.py

diff --git a/tests/a1_test/__init__.py b/tests/a1_test/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/tests/a2_test/__init__.py b/tests/a2_test/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/tests/a5_test/__init__.py b/tests/a5_test/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/tests/a6_test/__init__.py b/tests/a6_test/__init__.py
new file mode 100644
index 0000000..e69de29

From 99308db8efec02f5acd3059e5ed5cde9039ebcb0 Mon Sep 17 00:00:00 2001
From: nithishwer 
Date: Thu, 13 Nov 2025 16:00:00 +0000
Subject: [PATCH 32/40] fix: remove steus utils

---
 fepa/utils/steus_utils.py | 512 --------------------------------------
 1 file changed, 512 deletions(-)
 delete mode 100644 fepa/utils/steus_utils.py

diff --git a/fepa/utils/steus_utils.py b/fepa/utils/steus_utils.py
deleted file mode 100644
index 57fbcbc..0000000
--- a/fepa/utils/steus_utils.py
+++ /dev/null
@@ -1,512 +0,0 @@
-import pandas as pd
-import logging
-import shutil
-import gromacs
-import pyedr
-import os
-
-
-def make_test_mdp():
-    mdp = """
-    integrator              = md
-    dt                      = 0.002
-    nsteps                  = 250000  ; 500 ps
-    nstxout-compressed      = 5000
-    nstxout                 = 0
-    nstvout                 = 0
-    nstfout                 = 0
-    nstcalcenergy           = 50
-    nstenergy               = 50
-    nstlog                  = 5000
-    ;
-    cutoff-scheme           = Verlet
-    nstlist                 = 20
-    rlist                   = 0.9
-    vdwtype                 = Cut-off
-    vdw-modifier            = None
-    DispCorr                = EnerPres
-    rvdw                    = 0.9
-    coulombtype             = PME
-    rcoulomb                = 0.9
-    ;
-    tcoupl                  = v-rescale
-    tc_grps                 = POPC_FIP SOLV
-    tau_t                   = 1.0 1.0
-    ref_t                   = {T} {T}
-    ;
-    pcoupl                  = C-rescale
-    pcoupltype              = semiisotropic 
-    tau_p                   = 5.0
-    compressibility         = 4.5e-5  4.5e-5
-    ref_p                   = 1.0     1.0
-    refcoord_scaling        = com
-    ;
-    constraints             = h-bonds
-    constraint_algorithm    = LINCS
-    continuation            = no
-    gen-vel                 = yes
-    gen-temp                = {T}
-    ;
-    nstcomm                 = 100
-    comm_mode               = linear
-    comm_grps               = POPC_FIP SOLV
-
-    ; Pull code
-    pull                     = yes
-    pull_ncoords             = 1
-    pull_ngroups             = 2
-    pull_group1_name         = FIP
-    pull_group2_name         = POPC
-    pull_group2_pbcatom      = 995
-    pull_pbc_ref_prev_step_com = yes
-    pull_coord1_type         = umbrella
-    pull_coord1_geometry     = direction
-    pull_coord1_vec          = 0 0 1
-    pull_coord1_dim          = N N Y  ; pull along z
-    pull_coord1_groups       = 1 2
-    pull_coord1_start        = yes
-    pull_coord1_rate         = 0.0   ; restrain in place
-    pull_coord1_k            = 1000  ; kJ mol^-1 nm^-2
-    pull_nstfout             = 50
-
-    """
-
-    # These sometimes very small steps are tested because I might want to run STeUS
-    # either with delta T of 4 or of 5, and I'll need these simulations to calculate
-    # the weights, so I might as well run them all now while testing stability.
-    temps = ["310", "316", "322", "328", "334", "340", "346", "352", "358"]
-
-    for temp in temps:
-        with open(f"T{temp}.mdp", "w") as f:
-            f.write(mdp.format(T=temp))
-        os.system(
-            f"gmx grompp -f T{temp}.mdp -c ../../window_prep/alchembed13.gro -r ../../window_prep/alchembed13.gro -n ../../window_prep/index.ndx -p ../../window_prep/topol.top -o runs/T{temp}"
-        )
-
-
-def filter_colvar_by_temp(colvar_file, temp_time_file, T=310.0):
-    """
-    Reads a colvar file and a temp_time.csv file using pandas,
-    filters the colvar file based on the integer part of the time column
-    and the temperature in the temp_time file, and saves the filtered
-    data to a new colvar file.
-
-    Args:
-        colvar_file (str): Path to the input colvar file.
-        temp_time_file (str): Path to the input temp_time.csv file.
-        output_colvar_file (str): Path to the output filtered colvar file.
-        T (float): Temperature to filter the temp_time file. Default is 310.0.
-    """
-
-    # Read temp_time.csv using pandas
-    temp_df = pd.read_csv(temp_time_file)
-    # Filter temp_time DataFrame for Temperature (K) == T
-    filtered_temp_df = temp_df[temp_df["Temperature (K)"] == T]
-    # Get the set of integer times from the filtered temp_time DataFrame
-    valid_times = set(filtered_temp_df["Time (ps)"].astype(int))
-    # Read the colvar file using pandas
-    colvar_df = pd.read_csv(colvar_file, comment="#", sep="\s+", names=["time", "CV"])
-    # Filter the colvar DataFrame
-    filtered_colvar_df = colvar_df[colvar_df["time"].astype(int).isin(valid_times)]
-
-    # outfile name
-    output_colvar_file = colvar_file + ("_groundstate")
-    # Save the filtered colvar DataFrame to a new file
-    with open(output_colvar_file, "w") as outfile:
-        # Write the header
-        with open(colvar_file, "r") as infile:
-            for line in infile:
-                if line.startswith("#"):
-                    outfile.write(line)
-                    break
-        # Write the data
-        filtered_colvar_df.to_csv(outfile, sep=" ", index=False, header=False)
-
-    logging.info(f"Filtered colvar data saved to {output_colvar_file}")
-
-
-def filter_colvar_pandas(colvar_file, temp_time_file):
-    """
-    Reads a colvar file and a temp_time.csv file using pandas,
-    filters the colvar file based on the integer part of the time column
-    and the temperature in the temp_time file, and saves the filtered
-    data to a new colvar file.
-
-    Args:
-        colvar_file (str): Path to the input colvar file.
-        temp_time_file (str): Path to the input temp_time.csv file.
-        output_colvar_file (str): Path to the output filtered colvar file.
-    """
-
-    # Read temp_time.csv using pandas
-    temp_df = pd.read_csv(temp_time_file)
-    # Filter temp_time DataFrame for Temperature (K) == 310.0
-    filtered_temp_df = temp_df[temp_df["Temperature (K)"] == 310.0]
-    # Get the set of integer times from the filtered temp_time DataFrame
-    valid_times = set(filtered_temp_df["Time (ps)"].astype(int))
-    # Read the colvar file using pandas
-    colvar_df = pd.read_csv(colvar_file, comment="#", sep="\s+", names=["time", "CV"])
-    # Filter the colvar DataFrame
-    filtered_colvar_df = colvar_df[colvar_df["time"].astype(int).isin(valid_times)]
-
-    # outfile name
-    output_colvar_file = colvar_file + ("_groundstate")
-    # Save the filtered colvar DataFrame to a new file
-    with open(output_colvar_file, "w") as outfile:
-        # Write the header
-        with open(colvar_file, "r") as infile:
-            for line in infile:
-                if line.startswith("#"):
-                    outfile.write(line)
-                    break
-        # Write the data
-        filtered_colvar_df.to_csv(outfile, sep=" ", index=False, header=False)
-
-    logging.info(f"Filtered colvar data saved to {output_colvar_file}")
-
-
-def extract_temp_at_marker(filename):
-    """
-    Extracts the temperature value from the line preceding the "<<" marker
-    for each time step in the given gmx log file.
-
-    Args:
-        filename (str): The path to the gmx log file.
-
-    Returns:
-        list: A list of temperature values (floats) corresponding to each
-              occurrence of the "<<" marker.
-    """
-
-    temperatures = []
-    with open(filename, "r") as file:
-        lines = file.readlines()
-
-    for i, line in enumerate(lines):
-        if "<<" in line:
-            # Extract the temperature from the line before the marker
-            temp = float(line.split()[1])  # Assuming Temp is the 2nd value
-            temperatures.append(temp)
-    return temperatures
-
-
-def extract_time_at_marker(filename):
-    """
-    Extracts the time value before the line with the "<<" marker
-    for each time step in the given gromacs log file.
-
-    Args:
-        filename (str): The path to the gmx log file.
-
-    Returns:
-        list: A list of time values (floats) corresponding to each
-              occurrence of the "<<" marker.
-    """
-
-    times = []
-    with open(filename, "r") as file:
-        lines = file.readlines()
-
-    for i, line in enumerate(lines):
-        if "<<" in line:
-            # Look in reverse for a line with 'Time'
-            for j in range(i - 1, -1, -1):
-                if "Time" in lines[j]:
-                    # Extract the time value
-                    time = lines[j + 1].split()[
-                        1
-                    ]  # Assuming Time is the 2nd value in the line following 'Time'
-                    times.append(time)
-                    break
-    return times
-
-
-def write_temp_MDP(temp):
-    script_dir = os.path.dirname(
-        os.path.abspath(__file__)
-    )  # Get the directory of the script
-    template_path = os.path.join(
-        script_dir, "../data/template_steus_weight_sampling.mdp"
-    )
-    with open(template_path, "r") as template_file:
-        sample_mdp = template_file.read()
-    with open(f"T{temp}.mdp", "w") as f:
-        f.write(sample_mdp.format(T=temp))
-
-
-def relative_weight(i, betas, e_pots):
-    """
-    Code from Bjarne Feddersen
-    Assumes that i is the index of window n, and j is the index
-    of window n + 1. Returns the relative weight g(n+1) - g(n) as in eq. 6
-    of Sousa et al. https://doi.org/10.1021/acs.jctc.2c01162
-    """
-    j = i + 1
-    delta_beta = betas[j] - betas[i]
-    avg_e_pot = (e_pots[i] + e_pots[j]) / 2
-    return delta_beta * avg_e_pot
-
-
-def get_window_weights(temps, weight_sampling_path):
-    """
-    Code from Bjarne Feddersen
-
-    Get the window weights for the steus simulations from weight sampling
-    Args:
-        temps (list): List of temperatures
-        weight_sampling_path (str): Path to the weight sampling folder
-    Returns:
-        str: String of weights that can be copied into the mdp file
-    """
-    e_pots = []
-    R = 0.0083144621
-    betas = [1 / (R * int(T)) for T in temps]
-
-    for temp in temps:
-        edr_path = os.path.join(weight_sampling_path, f"T{temp}", f"T{temp}.edr")
-        edr = pyedr.edr_to_dict(edr_path)
-        avg_e_pot = edr["Potential"].mean()
-        e_pots.append(avg_e_pot)
-
-    rel_weights = [relative_weight(i, betas, e_pots) for i in range(len(e_pots) - 1)]
-
-    # first weight initially assumed to be 0
-    window_weights = [0]
-
-    # for each delta weight, add the delta to the last absolute weight to obtain
-    # the next absolute weight
-    for rel_weight in rel_weights:
-        window_weights.append(window_weights[-1] + rel_weight)
-
-    # subtract the average weight from each of the absolute weights
-    # to make the sum of the weights be equal to 0
-    # as seen in http://dx.doi.org/10.1103/PhysRevE.76.016703 below table 1
-    avg_weight = sum(window_weights) / len(window_weights)
-    window_weights = [weight - avg_weight for weight in window_weights]
-
-    weight_list = [(round(weight)) for weight in window_weights]
-    # write out string of weights that can be copied into the mdp file
-    return weight_list
-
-
-def add_steus_code_to_mdp(mdp_path, temp_weights, temp_lambdas):
-    steus_code = """\n
-;----------------------------------------------------
-; EE/STEUS Stuff
-;----------------------------------------------------
-free-energy              = expanded
-nstexpanded              = 500  ; attempt temperature change every 500 simulation steps
-lmc-stats                = wang-landau  ; update expanded ensemble weights with wang-landau algorithm
-lmc-move                 = metropolis
-init-wl-delta            = 1 ; initial value of the Wang-Landay incrementor in kT
-init-lambda-state        = 0
-init-lambda-weights      = {temp_weights} ; vector of the initial weights (free energies in kT) used for the expanded ensemble states. Vector of floats, length must match lambda vector lengths
-lmc-weights-equil        = wl-delta
-weight-equil-wl-delta    = 0.0001
-wl-ratio                 = 0.7
-wl-scale                 = 0.8
-wl-oneovert              = yes
-simulated-tempering      = yes
-sim-temp-low             = 310
-sim-temp-high            = 358
-temperature-lambdas      = {temp_lambdas}
-simulated-tempering-scaling = linear  ; linearly interpolates the temperatures using the values of temperature-lambdas
-"""
-
-    # Read the existing file content
-    with open(mdp_path, "r") as f:
-        lines = f.readlines()
-
-    # Find the first line that starts with 'free-energy' and truncate everything after it
-    new_lines = []
-    found = False
-    for line in lines:
-        if not found and line.strip().startswith("free-energy"):
-            found = True
-        if found:
-            break
-        new_lines.append(line)
-
-    # Write the truncated content back to the file
-    with open(mdp_path, "w") as f:
-        f.writelines(new_lines)
-
-    # Add the steus code to the end of the mdp file
-    with open(mdp_path, "a") as f:
-        f.write(
-            steus_code.format(
-                temp_weights=" ".join([str(weight) for weight in temp_weights]),
-                temp_lambdas=" ".join(temp_lambdas),
-            )
-        )
-
-
-def make_steus_plumed_dat(reus_plumed_dat, steus_plumed_dat, i):
-    """
-    Creates a modified version of the input file with only the i-th restraint value.
-
-    Args:
-        reus_plumed_dat (str): Path to the input file
-        steus_plumed_dat (str): Path where to save the modified file
-        i (int): Index of the restraint value to keep (0-based)
-
-    Returns:
-        str: Path to the created output file
-    """
-    with open(reus_plumed_dat, "r") as f:
-        content = f.readlines()
-
-    # Find and modify the restraint line
-    for idx, line in enumerate(content):
-        if line.strip().startswith("MOLINFO"):
-            parts = line.split()
-            for j, part in enumerate(parts):
-                if part == "STRUCTURE=../reference.pdb":
-                    # Modify to keep only the i-th value
-                    parts[j] = "STRUCTURE=reference.pdb"
-                    content[idx] = " ".join(parts) + "\n"
-        if line.strip().startswith("restraint: RESTRAINT"):
-            parts = line.split()
-
-            for j, part in enumerate(parts):
-                if part.startswith("AT=@replicas:"):
-                    # Extract and validate values
-                    values_str = part[len("AT=@replicas:") :]
-                    values = values_str.split(",")
-
-                    if i >= len(values):
-                        raise ValueError(
-                            f"Index {i} out of range. File has {len(values)} restraint values."
-                        )
-
-                    # Modify to keep only the i-th value
-                    parts[j] = f"AT={values[i]}"
-                    content[idx] = " ".join(parts) + "\n"
-                    break
-            break
-
-    # Write the modified content
-    with open(steus_plumed_dat, "w", encoding="UTF-8") as f:
-        f.writelines(content)
-
-    return values
-
-
-def setup_steus(
-    wdir_path,
-    plumed_path,
-    submission_script_template_arr,
-    template_mdp,
-    temp_weights,
-    temp_lambdas,
-    steus_name="steus_v1",
-    exist_ok=False,
-):
-    # Get name of parent dir of parent dir or wdir
-    run_name = os.path.basename(os.path.dirname(os.path.dirname(wdir_path)))
-
-    logging.info(f"Setting up STEUS for : {run_name}")
-
-    # Define the path to the 'boxes' folder
-    boxes_path = os.path.join(wdir_path, "boxes")
-
-    # Create 'steus' folder inside the working directory if it doesn't exist
-    steus_path = os.path.join(wdir_path, steus_name)
-    os.makedirs(steus_path, exist_ok=exist_ok)
-
-    # Get the list of folder names inside 'boxes'
-    sim_folders = [item for item in os.listdir(boxes_path) if item.startswith("sim")]
-
-    # Copy the array job script template
-    shutil.copy(
-        submission_script_template_arr,
-        os.path.join(steus_path, "job_ranv_steus_arr.sh"),
-    )
-
-    # Change job name in submission scripts
-    with open(
-        os.path.join(steus_path, "job_ranv_steus_arr.sh"),
-        "r",
-        encoding="utf-8",
-    ) as f:
-        lines = f.read()
-    with open(
-        os.path.join(steus_path, "job_ranv_steus_arr.sh"),
-        "w",
-        encoding="utf-8",
-    ) as f:
-        new_lines = lines.replace("$JOBNAME", run_name)
-        f.write(new_lines)
-
-    # Iterate through each folder in 'boxes' and replicate the folder structure in 'steus'
-    for sim_folder in sim_folders:
-        logging.info(f"Setting up STEUS for : {sim_folder}")
-        # Define the new folder path inside 'steus'
-        steus_sim_folder_path = os.path.join(steus_path, sim_folder)
-        os.makedirs(steus_sim_folder_path, exist_ok=True)
-
-        # Paths to the folders to be copied
-        toppar_source = os.path.join(boxes_path, sim_folder, "toppar")
-        amber_source = os.path.join(
-            boxes_path, sim_folder, "amber99sb-star-ildn-mut.ff"
-        )
-
-        # Copy 'toppar' and 'amber99sb-star-ildn-mut.ff' if they exist
-        if os.path.exists(toppar_source):
-            shutil.copytree(
-                toppar_source,
-                os.path.join(steus_sim_folder_path, "toppar"),
-                dirs_exist_ok=True,
-            )
-        if os.path.exists(amber_source):
-            shutil.copytree(
-                amber_source,
-                os.path.join(steus_sim_folder_path, "amber99sb-star-ildn-mut.ff"),
-                dirs_exist_ok=True,
-            )
-
-        # copy other inout files
-        shutil.copy(
-            os.path.join(boxes_path, sim_folder, "index.ndx"), steus_sim_folder_path
-        )
-        shutil.copy(
-            os.path.join(boxes_path, sim_folder, "prod.gro"),
-            os.path.join(steus_sim_folder_path, "equilibrated.gro"),
-        )
-        shutil.copy(
-            os.path.join(boxes_path, sim_folder, "topol.top"),
-            os.path.join(steus_sim_folder_path, "topol.top"),
-        )
-        # Copy template mdps
-        shutil.copy(
-            template_mdp,
-            os.path.join(steus_sim_folder_path, "prod.mdp"),
-        )
-
-        # Make steus plumed.dat
-        make_steus_plumed_dat(
-            plumed_path,
-            os.path.join(steus_sim_folder_path, "steus_plumed.dat"),
-            int(sim_folder.split("sim")[-1]),
-        )
-
-        logging.info(f"temp_weights: {temp_weights}")
-        logging.info(f"temp_lambdas: {temp_lambdas}")
-
-        # Add steus code to mdp
-        add_steus_code_to_mdp(
-            os.path.join(steus_sim_folder_path, "prod.mdp"),
-            temp_weights=temp_weights,
-            temp_lambdas=temp_lambdas,
-        )
-
-        # Make reference from equilibrated.gro
-        gromacs.editconf(
-            f=os.path.join(steus_sim_folder_path, "equilibrated.gro"),
-            o=os.path.join(steus_sim_folder_path, "reference.pdb"),
-        )
-
-        # Make tpr file
-
-    # Success
-    print(f"'steus' folder structure created successfully at {steus_path}")

From df28023a2be7bd7913981a262fffe0f6ca45a58a Mon Sep 17 00:00:00 2001
From: nithishwer 
Date: Thu, 13 Nov 2025 16:00:54 +0000
Subject: [PATCH 33/40] fix: remove steus mdp

---
 fepa/data/template_steus_weight_sampling.mdp | 82 --------------------
 1 file changed, 82 deletions(-)
 delete mode 100644 fepa/data/template_steus_weight_sampling.mdp

diff --git a/fepa/data/template_steus_weight_sampling.mdp b/fepa/data/template_steus_weight_sampling.mdp
deleted file mode 100644
index 8e6518c..0000000
--- a/fepa/data/template_steus_weight_sampling.mdp
+++ /dev/null
@@ -1,82 +0,0 @@
-;====================================================
-; NPT equilibration for vanilla
-;====================================================
-
-; RUN CONTROL
-;----------------------------------------------------
-define       = -DPOSRES 
-integrator   = md            ; stochastic leap-frog integrator
-nsteps       = 250000        ; 250000 x 0.002 = 500ps
-dt           = 0.002         ; 2 fs
-comm-mode    = Linear        ; remove center of mass translation
-nstcomm      = 100           ; frequency for center of mass motion removal
-comm-grps   = Protein_PA_PC_OL Water_and_ions
-
-; OUTPUT CONTROL
-;----------------------------------------------------
-nstxout          = 0      ; Dont save coordinates to .trr every 100 ps
-nstvout          = 0          ; don't save velocities to .trr
-nstfout          = 0          ; don't save forces to .trr
-refcoord_scaling = all        ; The reference coordinates are scaled with the scaling matrix of the pressure coupling
-nstxout-compressed     = 100000        ; xtc compressed trajectory output every 2 ps
-compressed-x-precision = 1000       ; precision with which to write to the compressed trajectory file
-nstlog                 = 1000       ; update log file every 2 ps
-nstenergy              = 1000       ; save energies every 2 ps
-nstcalcenergy          = 100        ; calculate energies every 200 fs
-
-; BONDS
-;----------------------------------------------------
-constraint_algorithm   = lincs    ; holonomic constraints
-constraints            = h-bonds  ; hydrogens only are constrained
-lincs-iter             = 1        ; accuracy of LINCS (1 is default)
-lincs-order            = 4        ; also related to accuracy (4 is default) 6 is needed for large time-steps with virtual sites or BD.
-lincs-warnangle        = 30       ; maximum angle that a bond can rotate before LINCS will complain (30 is default)
-continuation           = yes      ; formerly known as 'unconstrained-start' - yes for exact continuations and reruns
-
-; NEIGHBOR SEARCHING
-;----------------------------------------------------
-cutoff-scheme       = Verlet ; Default value
-ns-type             = grid   ; search neighboring grid cells
-nstlist             = 20     ; Frequency to update the neighbor list 20 fs (default is 10; 20 gives best performance)
-rlist               = 1.2    ; short-range neighborlist cutoff (in nm)
-pbc                 = xyz    ; 3D PBC
-
-; ELECTROSTATICS & EWALD
-;----------------------------------------------------
-coulombtype      = PME      ; Particle Mesh Ewald for long-range electrostatics
-rcoulomb         = 1.2      ; The distance for the Coulomb cut-off (in nm default is 1; lemkul lab uses 1.2)
-ewald_geometry   = 3d       ; Ewald sum is performed in all three dimensions (Suitable for non cubic boxes)
-pme-order        = 4        ; interpolation order for PME (default is 4)
-fourierspacing   = 0.12     ; grid spacing for FFT
-ewald-rtol       = 1e-6     ; relative strength of the Ewald-shifted direct potential at rcoulomb (default is 1e-5 Decreasing this will give a more accurate direct sum, but expensive)
-
-; VAN DER WAALS
-;----------------------------------------------------
-vdwtype                 = Cut-off                ; Default 
-vdw-modifier            = Potential-shift-Verlet ; (Default is none) Potential-shift-Verlet often used in conjunction with the Verlet neighbor search algorithm
-verlet-buffer-tolerance = 0.005                  ; Default
-rvdw                    = 1.2                    ; short-range van der Waals cutoff (in nm) (Default is 1)
-DispCorr                = EnerPres               ; apply long range dispersion corrections for Energy and Pressure (lemkul labs defualt)
-
-;----------------------------------------------------
-; TEMPERATURE & PRESSURE COUPL
-;----------------------------------------------------
-tcoupl           = V-rescale
-tc-grps          = Protein_PA_PC_OL Water_and_ions
-tau-t            = 2.0 2.0 ;
-ref-t            = {T} {T} ; 
-pcoupl           = Parrinello-Rahman   ; Better than berendson
-pcoupltype       = semiisotropic
-tau-p            = 2.0                         ; time constant (ps)
-ref-p            = 1.01325 1.01325             ; reference pressure (bar)
-compressibility  = 4.5e-05 4.5e-05             ; semiisothermal compressibility of water (bar^-1)
-
-;----------------------------------------------------
-; VELOCITY GENERATION
-;----------------------------------------------------
-gen_vel      = no      ; Velocity generation is off
-gen-temp     = {T}   ; Temperature for Maxwell distribution
-;----------------------------------------------------
-; FREE ENERGY
-;----------------------------------------------------
-free-energy              = no
\ No newline at end of file

From f755d69f32e44880f3aca526eb23fe9a64985efa Mon Sep 17 00:00:00 2001
From: nithishwer 
Date: Thu, 13 Nov 2025 16:09:07 +0000
Subject: [PATCH 34/40] fix: add new runtest in tests.yml actionm

---
 .github/workflows/test.yml | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index 5bb5d3e..67bdb61 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -34,4 +34,7 @@ jobs:
         run: conda run -n fepa_env pip install --no-cache-dir -e .
 
       - name: Run tests
-        run: conda run -n fepa_env pytest -v --maxfail=1 --disable-warnings
+        run: |
+          export PYTHONPATH=$(pwd)
+          conda run -n fepa_env pytest -v --maxfail=1 --disable-warnings
+

From 1cbb22253419c2f1096babc3dc3315c39d6e3c1b Mon Sep 17 00:00:00 2001
From: nithishwer 
Date: Thu, 13 Nov 2025 16:15:45 +0000
Subject: [PATCH 35/40] fix: remove steus flow

---
 fepa/flows/steus_flows.py | 317 --------------------------------------
 1 file changed, 317 deletions(-)
 delete mode 100644 fepa/flows/steus_flows.py

diff --git a/fepa/flows/steus_flows.py b/fepa/flows/steus_flows.py
deleted file mode 100644
index 36bdbad..0000000
--- a/fepa/flows/steus_flows.py
+++ /dev/null
@@ -1,317 +0,0 @@
-import os
-import gromacs
-import pandas as pd
-import logging
-import shutil
-from pathlib import Path
-from fepa.utils.wham_utils import (
-    parse_us_plumed_file,
-    process_wham_path,
-    plot_free_combined,
-    analyse_us_hist,
-)
-from fepa.utils.steus_utils import (
-    write_temp_MDP,
-    get_window_weights,
-    setup_steus,
-    extract_temp_at_marker,
-    extract_time_at_marker,
-    filter_colvar_pandas,
-)
-from fepa.utils.wham_utils import (
-    create_colvar_chunks,
-    generate_metadata,
-    parse_us_plumed_file,
-    plot_colvars,
-    plot_histogram,
-    process_colvars,
-)
-from multiprocessing import Pool
-
-
-class steus_weight_sampling_workflow:
-    """
-    Class to handle sampling weight for steus simulations.
-    """
-
-    def __init__(
-        self,
-        sim_paths,
-        temps=["310", "316", "322", "328", "334", "340", "346", "352", "358"],
-        weights_sampling_folder="weight_sampling",
-    ):
-        self.sim_paths = sim_paths
-        self.temps = temps
-        self.weights_sampling_folder = weights_sampling_folder
-
-    def setup_simulations(self, exist_ok=False):
-        """
-        Sets up the simulations for each sim path in sim_paths
-        by creating a folder called weight_sampling with subfolders
-        for each temperature in temps. It writes the MDP file
-        for each temperature and runs grompp for each
-        temperature.
-        """
-        for sim_path in self.sim_paths:
-            cwd = os.getcwd()
-            os.chdir(sim_path)
-
-            for temp in self.temps:
-                temp_folder = f"{self.weights_sampling_folder}/T{temp}"
-                os.makedirs(temp_folder, exist_ok=exist_ok)
-                logging.info(f"Creating {temp_folder}")
-                write_temp_MDP(temp)
-                logging.info(f"Writing MDP file for {temp}")
-                gromacs.grompp(
-                    f=f"T{temp}.mdp",
-                    c="prod.gro",
-                    r="prod.gro",
-                    n="index.ndx",
-                    p="topol.top",
-                    o=f"{temp_folder}/T{temp}",
-                )
-            os.chdir(cwd)
-
-    def run_simulations(self):
-        """
-        Runs the weight sampling simulations for each sim path
-        in sim_paths for each temperature in temps.
-        """
-        for sim_path in self.sim_paths:
-            cwd = os.getcwd()
-            os.chdir(sim_path)
-
-            for temp in self.temps:
-                temp_folder = f"{self.weights_sampling_folder}/T{temp}"
-                gromacs.mdrun(deffnm=f"{temp_folder}/T{temp}", v=True)
-                logging.info(f"Finished MD simulation for {temp}")
-            os.chdir(cwd)
-
-    def get_weights_from_simulations(self):
-        """
-        Gets the weights from the simulations for each sim path
-        in sim_paths
-        """
-        weight_lists = []  # List to store weight lists for each simulation; This is a list of lists
-        path_strings = []
-        for sim_path in self.sim_paths:
-            weight_sampling_path = os.path.join(sim_path, self.weights_sampling_folder)
-            weight_lists.append(get_window_weights(self.temps, weight_sampling_path))
-            path_strings.append(weight_sampling_path)
-        logging.info(f"Weight strings: {weight_lists}")
-        logging.info(f"Weight sampling paths: {path_strings}")
-        return weight_lists
-
-
-class steus_umbrella_sampling_workflow:
-    """
-    Class to handle umbrella sampling for steus simulations.
-    """
-
-    def __init__(
-        self,
-        wdir_path,
-        plumed_path,
-        submission_script_template_arr,
-        template_mdp,
-        temp_weights,
-        temp_lambdas,
-        start,
-        end,
-        steus_folder_name="steus_v1",
-        n_windows=24,
-        wham_dirname="wham",
-    ):
-        logging.info("Initializing steus_umbrella_sampling_workflow...")
-        self.wdir_path = wdir_path
-        self.plumed_path = plumed_path
-        self.submission_script_template_arr = submission_script_template_arr
-        self.template_mdp = template_mdp
-        self.temp_weights = temp_weights
-        self.temp_lambdas = temp_lambdas
-        self.steus_folder_name = steus_folder_name
-        self.n_windows = n_windows
-        self.steus_path = os.path.join(self.wdir_path, self.steus_folder_name)
-        self.wham_path = os.path.join(self.wdir_path, self.steus_folder_name, wham_dirname)
-        self.start = start
-        self.end = end
-
-    def setup_simulations(self, exist_ok=False):
-        setup_steus(
-            self.wdir_path,
-            self.plumed_path,
-            self.submission_script_template_arr,
-            self.template_mdp,
-            self.temp_weights,
-            self.temp_lambdas,
-            steus_name=self.steus_folder_name,
-            exist_ok=exist_ok,
-        )
-
-    def extract_ground_state_colvar(self):
-        logging.info("Extracting ground state colvar data...")
-        # For each window
-        for i in range(self.n_windows):
-            logging.info(f"Processing simulation {i}")
-            filename = f"{self.wdir_path}/{self.steus_folder_name}/sim{i}/sim{i}.log"  # Replace with the actual filename
-            # Extract temperature and time values from the gmx log file
-            temperatures = extract_temp_at_marker(filename)
-            times = extract_time_at_marker(filename)
-
-            logging.info(f"Extracted {len(temperatures)} temperature values.")
-            logging.info(f"Extracted {len(times)} time values.")
-
-            # Create a DataFrame from the extracted data
-            df = pd.DataFrame({"Time (ps)": times, "Temperature (K)": temperatures})
-
-            # Save the DataFrame to a CSV file
-            temp_file_name = filename.replace(".log", "_temp.csv")
-            df.to_csv(temp_file_name, index=False)
-            logging.info(f"DataFrame saved to {temp_file_name}.")
-
-            colvar_file = f"{self.wdir_path}/{self.steus_folder_name}/sim{i}/COLVAR"
-
-            # Filter the colvar file for the ground state temperature
-            filter_colvar_pandas(colvar_file, temp_file_name)
-
-    def prepare_wham(self):
-        # Create wham folder in steus path
-        logging.info("Preparing WHAM folder...")
-        os.makedirs(self.wham_path, exist_ok=True)
-        logging.info(f"Created WHAM folder at: {self.wham_path}")
-        colvar_path = os.path.join(self.wham_path, "colvars_100_pct_forward")
-        os.makedirs(colvar_path, exist_ok=True)
-        logging.info(f"Created COLVAR folder at: {colvar_path}")
-        # Copy all the colvar_groundstate files from all sims folder in steus_path to colvar_path
-        for i in range(0, self.n_windows):
-            sim_path = os.path.join(self.steus_path, f"sim{i}")
-            colvar_file = os.path.join(sim_path, "COLVAR_groundstate")
-            destination_path = os.path.join(colvar_path, f"COLVAR_groundstate_{i}")
-            if os.path.exists(colvar_file):
-                shutil.copy(colvar_file, destination_path)
-                logging.info(f"Copied {colvar_file} to {destination_path}")
-        # Plot colvars
-        plot_colvars(colvar_path)
-        # Process colvars
-        process_colvars(colvar_path=colvar_path, relaxation_time=500)
-        # Plot processed colvars
-        plot_histogram(colvar_path)
-        kappa_values = []
-        at_values = []
-        # Get kappa and at values from plumed files
-        for i in range(24):
-            plumed_file = os.path.join(self.steus_path, f"sim{i}", "steus_plumed.dat")
-            kappa, at = parse_us_plumed_file(plumed_file)
-            kappa_values.append(kappa)
-            at_values.append(at)
-        logging.info(f"Kappa values: {kappa_values}")
-        logging.info(f"AT values: {at_values}")
-        # Ensure all kappa values are the same
-        assert all(k == kappa_values[0] for k in kappa_values), (
-            "Kappa values are not consistent."
-        )
-        generate_metadata(
-            colvar_path,
-            kappa_values[0],
-            at_values,
-            output_file="metadata.dat",
-            decorrelation_time=0.02,
-        )
-        ## Process each chunk size in both forward and reverse directions
-        for chunk_size in [20, 40, 60, 80]:
-            print(f"Processing chunk size: {chunk_size}% for {colvar_path}")
-            # Forward direction
-            forward_chunk_path = create_colvar_chunks(
-                wham_path=self.wham_path,
-                colvar_100_pct_path=colvar_path,
-                chunk_size_percentage=chunk_size,
-                direction="forward",
-            )
-            # Call the plotting function
-            plot_histogram(forward_chunk_path)
-            # Create metadata file
-            generate_metadata(
-                forward_chunk_path,
-                kappa_values[0],
-                at_values,
-                output_file="metadata.dat",
-                decorrelation_time=0.02,
-            )
-            # Reverse direction
-            reverse_chunk_path = create_colvar_chunks(
-                wham_path=self.wham_path,
-                colvar_100_pct_path=colvar_path,
-                chunk_size_percentage=chunk_size,
-                direction="reverse",
-            )
-            # Call the plotting function
-            plot_histogram(reverse_chunk_path)
-            # Create metadata file
-            generate_metadata(
-                reverse_chunk_path,
-                kappa_values[0],
-                at_values,
-                output_file="metadata.dat",
-                decorrelation_time=0.02,
-            )
-
-    def run_wham(self):
-        if not os.path.exists(self.wham_path):
-            # error
-            raise FileNotFoundError(f"WHAM path does not exist: {self.wham_path}")
-        # Create a list to store prepared WHAM paths
-        prepared_wham_paths = []
-        prepared_wham_paths.append(
-            os.path.join(self.wham_path, "colvars_100_pct_forward")
-        )
-        for pct in ["20", "40", "60", "80"]:
-            prepared_wham_paths.append(
-                os.path.join(self.wham_path, f"colvars_{pct}_pct_forward")
-            )
-            prepared_wham_paths.append(
-                os.path.join(self.wham_path, f"colvars_{pct}_pct_reverse")
-            )
-        # Print prepared WHAM paths
-        logging.info(f"Prepared WHAM paths: {prepared_wham_paths}")
-
-        # Get at values
-        kappa_values = []
-        at_values = []
-        # Get kappa and at values from plumed files
-        for i in range(self.n_windows):
-            plumed_file = os.path.join(self.steus_path, f"sim{i}", "steus_plumed.dat")
-            kappa, at = parse_us_plumed_file(plumed_file)
-            kappa_values.append(kappa)
-            at_values.append(at)
-
-        # Use Pool to parallelize WHAM runs
-        with Pool(processes=2) as pool:
-            # Map each wham_path to the process_wham_path function
-            pool.starmap(
-                process_wham_path,
-                [(wham_path, at_values) for wham_path in prepared_wham_paths],
-            )
-
-    def analyse_us_hist(self, range=(90, 180), colvar_filename="COLVAR"):
-        # Make us_path
-        us_path = os.path.join(self.wdir_path, self.steus_folder_name)
-
-        # Analyse
-        analyse_us_hist(
-            us_path,
-            colvar_filename=colvar_filename,
-            label=f"{self.start}_{self.end}",
-            range=range,
-        )
-
-    def plot_free_energies(
-        self,
-        units="kcal",
-    ):
-        plot_free_combined(
-            self.wham_path,
-            structure_1=self.start,
-            structure_2=self.end,
-            units=units,
-            box_CV_means_csv=os.path.join(self.steus_path, "mean_values.csv"),
-        )

From 9b53b99aad0897335da696e203d4ed4a33bd2bd1 Mon Sep 17 00:00:00 2001
From: nithishwer 
Date: Thu, 13 Nov 2025 16:20:36 +0000
Subject: [PATCH 36/40] fix: remove steus flow from init

---
 fepa/flows/__init__.py | 5 -----
 1 file changed, 5 deletions(-)

diff --git a/fepa/flows/__init__.py b/fepa/flows/__init__.py
index 3e85e6b..56020d4 100644
--- a/fepa/flows/__init__.py
+++ b/fepa/flows/__init__.py
@@ -5,11 +5,6 @@
     memento_workflow,
 )
 
-from .steus_flows import (
-    steus_weight_sampling_workflow,
-    steus_umbrella_sampling_workflow,
-)
-
 from .torsions_flow import (
     torsions_analysis_workflow,
 )

From 2c8bddbaa2470a8cc970541a672b05f06d88ec2f Mon Sep 17 00:00:00 2001
From: nithishwer 
Date: Sun, 23 Nov 2025 15:57:50 +0000
Subject: [PATCH 37/40] fix: reduce comp tolerance

---
 .../a1_test/1_binding_pocket_analysis_test.py |  5 ----
 tests/utils.py                                | 27 ++++++++++++-------
 2 files changed, 18 insertions(+), 14 deletions(-)

diff --git a/tests/a1_test/1_binding_pocket_analysis_test.py b/tests/a1_test/1_binding_pocket_analysis_test.py
index 8982ef5..efa7b01 100644
--- a/tests/a1_test/1_binding_pocket_analysis_test.py
+++ b/tests/a1_test/1_binding_pocket_analysis_test.py
@@ -29,11 +29,6 @@
 
 import pytest
 
-# Default numeric tolerances for floating-point comparisons
-DECIMALS = 6
-RTOL = 1e-6
-ATOL = 1e-8
-
 # Ensure Literal works under older Python versions
 builtins.Literal = Literal
 
diff --git a/tests/utils.py b/tests/utils.py
index 1b654eb..0fe27d2 100644
--- a/tests/utils.py
+++ b/tests/utils.py
@@ -4,9 +4,9 @@
 import pandas.testing as pdt
 
 # === numeric comparison tolerances (module-level defaults)
-DECIMALS = 6
-RTOL = 1e-6
-ATOL = 1e-8
+DECIMALS = 3
+RTOL = 1e-3
+ATOL = 1e-3
 
 
 # === small helpers for deterministic comparisons
@@ -25,6 +25,7 @@ def first_existing(*paths: Path) -> Path:
             return p
     raise FileNotFoundError(f"None of the candidate files exist: {paths!r}")
 
+
 def align_on_common_keys(
     left: pd.DataFrame,
     right: pd.DataFrame,
@@ -34,7 +35,9 @@ def align_on_common_keys(
     keys = [k for k in keys if k in left.columns and k in right.columns]
     if not keys:
         n = min(len(left), len(right))
-        return left.iloc[:n].reset_index(drop=True), right.iloc[:n].reset_index(drop=True)
+        return left.iloc[:n].reset_index(drop=True), right.iloc[:n].reset_index(
+            drop=True
+        )
     merged_keys = right[keys].drop_duplicates()
     left2 = left.merge(merged_keys, on=keys, how="inner")
     right2 = right.merge(left2[keys].drop_duplicates(), on=keys, how="inner")
@@ -42,16 +45,21 @@ def align_on_common_keys(
     right2 = sort_by(right2, keys)
     return left2.reset_index(drop=True), right2.reset_index(drop=True)
 
+
 def sort_by(df: pd.DataFrame, keys=("ensemble", "timestep", "frame")) -> pd.DataFrame:
     """Sort a DataFrame by available keys, resetting the index."""
     use = [k for k in keys if k in df.columns]
     return df.sort_values(use).reset_index(drop=True) if use else df
 
 
-def align_pc_signs(actual: pd.DataFrame, expected: pd.DataFrame, pc_prefix: str = "PC") -> pd.DataFrame:
+def align_pc_signs(
+    actual: pd.DataFrame, expected: pd.DataFrame, pc_prefix: str = "PC"
+) -> pd.DataFrame:
     """Flip PC signs in 'expected' to match 'actual' for consistent orientation."""
     exp = expected.copy()
-    pc_cols = [c for c in exp.columns if c.startswith(pc_prefix) and c in actual.columns]
+    pc_cols = [
+        c for c in exp.columns if c.startswith(pc_prefix) and c in actual.columns
+    ]
     for c in pc_cols:
         a = actual[c].to_numpy()
         b = exp[c].to_numpy()
@@ -87,8 +95,7 @@ def check_csv_equality(
         f"Expected columns: {list(exp.columns)}"
     )
     assert len(act) == len(exp), (
-        f"{label}: row count mismatch "
-        f"(actual={len(act)}, expected={len(exp)})."
+        f"{label}: row count mismatch (actual={len(act)}, expected={len(exp)})."
     )
 
     # --- Numeric comparison
@@ -102,4 +109,6 @@ def check_csv_equality(
             check_dtype=False,
         )
     except AssertionError as e:
-        raise AssertionError(f"{label}: numeric values differ beyond tolerance.\n{e}") from e
+        raise AssertionError(
+            f"{label}: numeric values differ beyond tolerance.\n{e}"
+        ) from e

From 39f29dc53939cf8c41c0c2b659b2504925b82134 Mon Sep 17 00:00:00 2001
From: nithishwer 
Date: Sun, 23 Nov 2025 16:07:51 +0000
Subject: [PATCH 38/40] fix: add badge to readme

---
 README.md | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/README.md b/README.md
index d4961b2..f26c5e4 100644
--- a/README.md
+++ b/README.md
@@ -1,6 +1,9 @@
 
 # Getting Started with FEPA
 
+![Tests](https://github.com/Nithishwer/FEPA/actions/workflows/tests.yml/badge.svg)
+
+
 **FEPA** (Free Energy Perturbation Analysis) is a Python package for analyzing molecular dynamics (MD) trajectories from FEP simulations, particularly **ABFEs**. FEPA allows you to visualize conformational changes and set up simulations to correct free energy estimates.
 
 This guide covers installation, basic usage, and key workflows.

From 43d0fb9a6aef7459c4740d3c36ba4d05a837b63c Mon Sep 17 00:00:00 2001
From: nithishwer 
Date: Sun, 23 Nov 2025 16:09:47 +0000
Subject: [PATCH 39/40] fix: add badge to readme

---
 README.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/README.md b/README.md
index f26c5e4..ef13282 100644
--- a/README.md
+++ b/README.md
@@ -1,7 +1,7 @@
 
 # Getting Started with FEPA
 
-![Tests](https://github.com/Nithishwer/FEPA/actions/workflows/tests.yml/badge.svg)
+[![Tests](https://github.com/Nithishwer/FEPA/actions/workflows/test.yml/badge.svg)](https://github.com/Nithishwer/FEPA/actions/workflows/test.yml)
 
 
 **FEPA** (Free Energy Perturbation Analysis) is a Python package for analyzing molecular dynamics (MD) trajectories from FEP simulations, particularly **ABFEs**. FEPA allows you to visualize conformational changes and set up simulations to correct free energy estimates.

From a4a414e3362c11792ce4b2f5553019265e32547d Mon Sep 17 00:00:00 2001
From: Nithishwer Mouroug Anand <31736866+Nithishwer@users.noreply.github.com>
Date: Sun, 23 Nov 2025 16:15:47 +0000
Subject: [PATCH 40/40] Add LICENSE file

---
 LICENSE | 695 ++++++++++++++++++++++++++++++++++++++++++++++++++++++--
 1 file changed, 674 insertions(+), 21 deletions(-)

diff --git a/LICENSE b/LICENSE
index ea32032..f288702 100644
--- a/LICENSE
+++ b/LICENSE
@@ -1,21 +1,674 @@
-MIT License
-
-Copyright (c) 2025 Nithishwer Mouroug Anand
-
-Permission is hereby granted, free of charge, to any person obtaining a copy
-of this software and associated documentation files (the "Software"), to deal
-in the Software without restriction, including without limitation the rights
-to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
-copies of the Software, and to permit persons to whom the Software is
-furnished to do so, subject to the following conditions:
-
-The above copyright notice and this permission notice shall be included in all
-copies or substantial portions of the Software.
-
-THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
-SOFTWARE.
+                    GNU GENERAL PUBLIC LICENSE
+                       Version 3, 29 June 2007
+
+ Copyright (C) 2007 Free Software Foundation, Inc. 
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+                            Preamble
+
+  The GNU General Public License is a free, copyleft license for
+software and other kinds of works.
+
+  The licenses for most software and other practical works are designed
+to take away your freedom to share and change the works.  By contrast,
+the GNU General Public License is intended to guarantee your freedom to
+share and change all versions of a program--to make sure it remains free
+software for all its users.  We, the Free Software Foundation, use the
+GNU General Public License for most of our software; it applies also to
+any other work released this way by its authors.  You can apply it to
+your programs, too.
+
+  When we speak of free software, we are referring to freedom, not
+price.  Our General Public Licenses are designed to make sure that you
+have the freedom to distribute copies of free software (and charge for
+them if you wish), that you receive source code or can get it if you
+want it, that you can change the software or use pieces of it in new
+free programs, and that you know you can do these things.
+
+  To protect your rights, we need to prevent others from denying you
+these rights or asking you to surrender the rights.  Therefore, you have
+certain responsibilities if you distribute copies of the software, or if
+you modify it: responsibilities to respect the freedom of others.
+
+  For example, if you distribute copies of such a program, whether
+gratis or for a fee, you must pass on to the recipients the same
+freedoms that you received.  You must make sure that they, too, receive
+or can get the source code.  And you must show them these terms so they
+know their rights.
+
+  Developers that use the GNU GPL protect your rights with two steps:
+(1) assert copyright on the software, and (2) offer you this License
+giving you legal permission to copy, distribute and/or modify it.
+
+  For the developers' and authors' protection, the GPL clearly explains
+that there is no warranty for this free software.  For both users' and
+authors' sake, the GPL requires that modified versions be marked as
+changed, so that their problems will not be attributed erroneously to
+authors of previous versions.
+
+  Some devices are designed to deny users access to install or run
+modified versions of the software inside them, although the manufacturer
+can do so.  This is fundamentally incompatible with the aim of
+protecting users' freedom to change the software.  The systematic
+pattern of such abuse occurs in the area of products for individuals to
+use, which is precisely where it is most unacceptable.  Therefore, we
+have designed this version of the GPL to prohibit the practice for those
+products.  If such problems arise substantially in other domains, we
+stand ready to extend this provision to those domains in future versions
+of the GPL, as needed to protect the freedom of users.
+
+  Finally, every program is threatened constantly by software patents.
+States should not allow patents to restrict development and use of
+software on general-purpose computers, but in those that do, we wish to
+avoid the special danger that patents applied to a free program could
+make it effectively proprietary.  To prevent this, the GPL assures that
+patents cannot be used to render the program non-free.
+
+  The precise terms and conditions for copying, distribution and
+modification follow.
+
+                       TERMS AND CONDITIONS
+
+  0. Definitions.
+
+  "This License" refers to version 3 of the GNU General Public License.
+
+  "Copyright" also means copyright-like laws that apply to other kinds of
+works, such as semiconductor masks.
+
+  "The Program" refers to any copyrightable work licensed under this
+License.  Each licensee is addressed as "you".  "Licensees" and
+"recipients" may be individuals or organizations.
+
+  To "modify" a work means to copy from or adapt all or part of the work
+in a fashion requiring copyright permission, other than the making of an
+exact copy.  The resulting work is called a "modified version" of the
+earlier work or a work "based on" the earlier work.
+
+  A "covered work" means either the unmodified Program or a work based
+on the Program.
+
+  To "propagate" a work means to do anything with it that, without
+permission, would make you directly or secondarily liable for
+infringement under applicable copyright law, except executing it on a
+computer or modifying a private copy.  Propagation includes copying,
+distribution (with or without modification), making available to the
+public, and in some countries other activities as well.
+
+  To "convey" a work means any kind of propagation that enables other
+parties to make or receive copies.  Mere interaction with a user through
+a computer network, with no transfer of a copy, is not conveying.
+
+  An interactive user interface displays "Appropriate Legal Notices"
+to the extent that it includes a convenient and prominently visible
+feature that (1) displays an appropriate copyright notice, and (2)
+tells the user that there is no warranty for the work (except to the
+extent that warranties are provided), that licensees may convey the
+work under this License, and how to view a copy of this License.  If
+the interface presents a list of user commands or options, such as a
+menu, a prominent item in the list meets this criterion.
+
+  1. Source Code.
+
+  The "source code" for a work means the preferred form of the work
+for making modifications to it.  "Object code" means any non-source
+form of a work.
+
+  A "Standard Interface" means an interface that either is an official
+standard defined by a recognized standards body, or, in the case of
+interfaces specified for a particular programming language, one that
+is widely used among developers working in that language.
+
+  The "System Libraries" of an executable work include anything, other
+than the work as a whole, that (a) is included in the normal form of
+packaging a Major Component, but which is not part of that Major
+Component, and (b) serves only to enable use of the work with that
+Major Component, or to implement a Standard Interface for which an
+implementation is available to the public in source code form.  A
+"Major Component", in this context, means a major essential component
+(kernel, window system, and so on) of the specific operating system
+(if any) on which the executable work runs, or a compiler used to
+produce the work, or an object code interpreter used to run it.
+
+  The "Corresponding Source" for a work in object code form means all
+the source code needed to generate, install, and (for an executable
+work) run the object code and to modify the work, including scripts to
+control those activities.  However, it does not include the work's
+System Libraries, or general-purpose tools or generally available free
+programs which are used unmodified in performing those activities but
+which are not part of the work.  For example, Corresponding Source
+includes interface definition files associated with source files for
+the work, and the source code for shared libraries and dynamically
+linked subprograms that the work is specifically designed to require,
+such as by intimate data communication or control flow between those
+subprograms and other parts of the work.
+
+  The Corresponding Source need not include anything that users
+can regenerate automatically from other parts of the Corresponding
+Source.
+
+  The Corresponding Source for a work in source code form is that
+same work.
+
+  2. Basic Permissions.
+
+  All rights granted under this License are granted for the term of
+copyright on the Program, and are irrevocable provided the stated
+conditions are met.  This License explicitly affirms your unlimited
+permission to run the unmodified Program.  The output from running a
+covered work is covered by this License only if the output, given its
+content, constitutes a covered work.  This License acknowledges your
+rights of fair use or other equivalent, as provided by copyright law.
+
+  You may make, run and propagate covered works that you do not
+convey, without conditions so long as your license otherwise remains
+in force.  You may convey covered works to others for the sole purpose
+of having them make modifications exclusively for you, or provide you
+with facilities for running those works, provided that you comply with
+the terms of this License in conveying all material for which you do
+not control copyright.  Those thus making or running the covered works
+for you must do so exclusively on your behalf, under your direction
+and control, on terms that prohibit them from making any copies of
+your copyrighted material outside their relationship with you.
+
+  Conveying under any other circumstances is permitted solely under
+the conditions stated below.  Sublicensing is not allowed; section 10
+makes it unnecessary.
+
+  3. Protecting Users' Legal Rights From Anti-Circumvention Law.
+
+  No covered work shall be deemed part of an effective technological
+measure under any applicable law fulfilling obligations under article
+11 of the WIPO copyright treaty adopted on 20 December 1996, or
+similar laws prohibiting or restricting circumvention of such
+measures.
+
+  When you convey a covered work, you waive any legal power to forbid
+circumvention of technological measures to the extent such circumvention
+is effected by exercising rights under this License with respect to
+the covered work, and you disclaim any intention to limit operation or
+modification of the work as a means of enforcing, against the work's
+users, your or third parties' legal rights to forbid circumvention of
+technological measures.
+
+  4. Conveying Verbatim Copies.
+
+  You may convey verbatim copies of the Program's source code as you
+receive it, in any medium, provided that you conspicuously and
+appropriately publish on each copy an appropriate copyright notice;
+keep intact all notices stating that this License and any
+non-permissive terms added in accord with section 7 apply to the code;
+keep intact all notices of the absence of any warranty; and give all
+recipients a copy of this License along with the Program.
+
+  You may charge any price or no price for each copy that you convey,
+and you may offer support or warranty protection for a fee.
+
+  5. Conveying Modified Source Versions.
+
+  You may convey a work based on the Program, or the modifications to
+produce it from the Program, in the form of source code under the
+terms of section 4, provided that you also meet all of these conditions:
+
+    a) The work must carry prominent notices stating that you modified
+    it, and giving a relevant date.
+
+    b) The work must carry prominent notices stating that it is
+    released under this License and any conditions added under section
+    7.  This requirement modifies the requirement in section 4 to
+    "keep intact all notices".
+
+    c) You must license the entire work, as a whole, under this
+    License to anyone who comes into possession of a copy.  This
+    License will therefore apply, along with any applicable section 7
+    additional terms, to the whole of the work, and all its parts,
+    regardless of how they are packaged.  This License gives no
+    permission to license the work in any other way, but it does not
+    invalidate such permission if you have separately received it.
+
+    d) If the work has interactive user interfaces, each must display
+    Appropriate Legal Notices; however, if the Program has interactive
+    interfaces that do not display Appropriate Legal Notices, your
+    work need not make them do so.
+
+  A compilation of a covered work with other separate and independent
+works, which are not by their nature extensions of the covered work,
+and which are not combined with it such as to form a larger program,
+in or on a volume of a storage or distribution medium, is called an
+"aggregate" if the compilation and its resulting copyright are not
+used to limit the access or legal rights of the compilation's users
+beyond what the individual works permit.  Inclusion of a covered work
+in an aggregate does not cause this License to apply to the other
+parts of the aggregate.
+
+  6. Conveying Non-Source Forms.
+
+  You may convey a covered work in object code form under the terms
+of sections 4 and 5, provided that you also convey the
+machine-readable Corresponding Source under the terms of this License,
+in one of these ways:
+
+    a) Convey the object code in, or embodied in, a physical product
+    (including a physical distribution medium), accompanied by the
+    Corresponding Source fixed on a durable physical medium
+    customarily used for software interchange.
+
+    b) Convey the object code in, or embodied in, a physical product
+    (including a physical distribution medium), accompanied by a
+    written offer, valid for at least three years and valid for as
+    long as you offer spare parts or customer support for that product
+    model, to give anyone who possesses the object code either (1) a
+    copy of the Corresponding Source for all the software in the
+    product that is covered by this License, on a durable physical
+    medium customarily used for software interchange, for a price no
+    more than your reasonable cost of physically performing this
+    conveying of source, or (2) access to copy the
+    Corresponding Source from a network server at no charge.
+
+    c) Convey individual copies of the object code with a copy of the
+    written offer to provide the Corresponding Source.  This
+    alternative is allowed only occasionally and noncommercially, and
+    only if you received the object code with such an offer, in accord
+    with subsection 6b.
+
+    d) Convey the object code by offering access from a designated
+    place (gratis or for a charge), and offer equivalent access to the
+    Corresponding Source in the same way through the same place at no
+    further charge.  You need not require recipients to copy the
+    Corresponding Source along with the object code.  If the place to
+    copy the object code is a network server, the Corresponding Source
+    may be on a different server (operated by you or a third party)
+    that supports equivalent copying facilities, provided you maintain
+    clear directions next to the object code saying where to find the
+    Corresponding Source.  Regardless of what server hosts the
+    Corresponding Source, you remain obligated to ensure that it is
+    available for as long as needed to satisfy these requirements.
+
+    e) Convey the object code using peer-to-peer transmission, provided
+    you inform other peers where the object code and Corresponding
+    Source of the work are being offered to the general public at no
+    charge under subsection 6d.
+
+  A separable portion of the object code, whose source code is excluded
+from the Corresponding Source as a System Library, need not be
+included in conveying the object code work.
+
+  A "User Product" is either (1) a "consumer product", which means any
+tangible personal property which is normally used for personal, family,
+or household purposes, or (2) anything designed or sold for incorporation
+into a dwelling.  In determining whether a product is a consumer product,
+doubtful cases shall be resolved in favor of coverage.  For a particular
+product received by a particular user, "normally used" refers to a
+typical or common use of that class of product, regardless of the status
+of the particular user or of the way in which the particular user
+actually uses, or expects or is expected to use, the product.  A product
+is a consumer product regardless of whether the product has substantial
+commercial, industrial or non-consumer uses, unless such uses represent
+the only significant mode of use of the product.
+
+  "Installation Information" for a User Product means any methods,
+procedures, authorization keys, or other information required to install
+and execute modified versions of a covered work in that User Product from
+a modified version of its Corresponding Source.  The information must
+suffice to ensure that the continued functioning of the modified object
+code is in no case prevented or interfered with solely because
+modification has been made.
+
+  If you convey an object code work under this section in, or with, or
+specifically for use in, a User Product, and the conveying occurs as
+part of a transaction in which the right of possession and use of the
+User Product is transferred to the recipient in perpetuity or for a
+fixed term (regardless of how the transaction is characterized), the
+Corresponding Source conveyed under this section must be accompanied
+by the Installation Information.  But this requirement does not apply
+if neither you nor any third party retains the ability to install
+modified object code on the User Product (for example, the work has
+been installed in ROM).
+
+  The requirement to provide Installation Information does not include a
+requirement to continue to provide support service, warranty, or updates
+for a work that has been modified or installed by the recipient, or for
+the User Product in which it has been modified or installed.  Access to a
+network may be denied when the modification itself materially and
+adversely affects the operation of the network or violates the rules and
+protocols for communication across the network.
+
+  Corresponding Source conveyed, and Installation Information provided,
+in accord with this section must be in a format that is publicly
+documented (and with an implementation available to the public in
+source code form), and must require no special password or key for
+unpacking, reading or copying.
+
+  7. Additional Terms.
+
+  "Additional permissions" are terms that supplement the terms of this
+License by making exceptions from one or more of its conditions.
+Additional permissions that are applicable to the entire Program shall
+be treated as though they were included in this License, to the extent
+that they are valid under applicable law.  If additional permissions
+apply only to part of the Program, that part may be used separately
+under those permissions, but the entire Program remains governed by
+this License without regard to the additional permissions.
+
+  When you convey a copy of a covered work, you may at your option
+remove any additional permissions from that copy, or from any part of
+it.  (Additional permissions may be written to require their own
+removal in certain cases when you modify the work.)  You may place
+additional permissions on material, added by you to a covered work,
+for which you have or can give appropriate copyright permission.
+
+  Notwithstanding any other provision of this License, for material you
+add to a covered work, you may (if authorized by the copyright holders of
+that material) supplement the terms of this License with terms:
+
+    a) Disclaiming warranty or limiting liability differently from the
+    terms of sections 15 and 16 of this License; or
+
+    b) Requiring preservation of specified reasonable legal notices or
+    author attributions in that material or in the Appropriate Legal
+    Notices displayed by works containing it; or
+
+    c) Prohibiting misrepresentation of the origin of that material, or
+    requiring that modified versions of such material be marked in
+    reasonable ways as different from the original version; or
+
+    d) Limiting the use for publicity purposes of names of licensors or
+    authors of the material; or
+
+    e) Declining to grant rights under trademark law for use of some
+    trade names, trademarks, or service marks; or
+
+    f) Requiring indemnification of licensors and authors of that
+    material by anyone who conveys the material (or modified versions of
+    it) with contractual assumptions of liability to the recipient, for
+    any liability that these contractual assumptions directly impose on
+    those licensors and authors.
+
+  All other non-permissive additional terms are considered "further
+restrictions" within the meaning of section 10.  If the Program as you
+received it, or any part of it, contains a notice stating that it is
+governed by this License along with a term that is a further
+restriction, you may remove that term.  If a license document contains
+a further restriction but permits relicensing or conveying under this
+License, you may add to a covered work material governed by the terms
+of that license document, provided that the further restriction does
+not survive such relicensing or conveying.
+
+  If you add terms to a covered work in accord with this section, you
+must place, in the relevant source files, a statement of the
+additional terms that apply to those files, or a notice indicating
+where to find the applicable terms.
+
+  Additional terms, permissive or non-permissive, may be stated in the
+form of a separately written license, or stated as exceptions;
+the above requirements apply either way.
+
+  8. Termination.
+
+  You may not propagate or modify a covered work except as expressly
+provided under this License.  Any attempt otherwise to propagate or
+modify it is void, and will automatically terminate your rights under
+this License (including any patent licenses granted under the third
+paragraph of section 11).
+
+  However, if you cease all violation of this License, then your
+license from a particular copyright holder is reinstated (a)
+provisionally, unless and until the copyright holder explicitly and
+finally terminates your license, and (b) permanently, if the copyright
+holder fails to notify you of the violation by some reasonable means
+prior to 60 days after the cessation.
+
+  Moreover, your license from a particular copyright holder is
+reinstated permanently if the copyright holder notifies you of the
+violation by some reasonable means, this is the first time you have
+received notice of violation of this License (for any work) from that
+copyright holder, and you cure the violation prior to 30 days after
+your receipt of the notice.
+
+  Termination of your rights under this section does not terminate the
+licenses of parties who have received copies or rights from you under
+this License.  If your rights have been terminated and not permanently
+reinstated, you do not qualify to receive new licenses for the same
+material under section 10.
+
+  9. Acceptance Not Required for Having Copies.
+
+  You are not required to accept this License in order to receive or
+run a copy of the Program.  Ancillary propagation of a covered work
+occurring solely as a consequence of using peer-to-peer transmission
+to receive a copy likewise does not require acceptance.  However,
+nothing other than this License grants you permission to propagate or
+modify any covered work.  These actions infringe copyright if you do
+not accept this License.  Therefore, by modifying or propagating a
+covered work, you indicate your acceptance of this License to do so.
+
+  10. Automatic Licensing of Downstream Recipients.
+
+  Each time you convey a covered work, the recipient automatically
+receives a license from the original licensors, to run, modify and
+propagate that work, subject to this License.  You are not responsible
+for enforcing compliance by third parties with this License.
+
+  An "entity transaction" is a transaction transferring control of an
+organization, or substantially all assets of one, or subdividing an
+organization, or merging organizations.  If propagation of a covered
+work results from an entity transaction, each party to that
+transaction who receives a copy of the work also receives whatever
+licenses to the work the party's predecessor in interest had or could
+give under the previous paragraph, plus a right to possession of the
+Corresponding Source of the work from the predecessor in interest, if
+the predecessor has it or can get it with reasonable efforts.
+
+  You may not impose any further restrictions on the exercise of the
+rights granted or affirmed under this License.  For example, you may
+not impose a license fee, royalty, or other charge for exercise of
+rights granted under this License, and you may not initiate litigation
+(including a cross-claim or counterclaim in a lawsuit) alleging that
+any patent claim is infringed by making, using, selling, offering for
+sale, or importing the Program or any portion of it.
+
+  11. Patents.
+
+  A "contributor" is a copyright holder who authorizes use under this
+License of the Program or a work on which the Program is based.  The
+work thus licensed is called the contributor's "contributor version".
+
+  A contributor's "essential patent claims" are all patent claims
+owned or controlled by the contributor, whether already acquired or
+hereafter acquired, that would be infringed by some manner, permitted
+by this License, of making, using, or selling its contributor version,
+but do not include claims that would be infringed only as a
+consequence of further modification of the contributor version.  For
+purposes of this definition, "control" includes the right to grant
+patent sublicenses in a manner consistent with the requirements of
+this License.
+
+  Each contributor grants you a non-exclusive, worldwide, royalty-free
+patent license under the contributor's essential patent claims, to
+make, use, sell, offer for sale, import and otherwise run, modify and
+propagate the contents of its contributor version.
+
+  In the following three paragraphs, a "patent license" is any express
+agreement or commitment, however denominated, not to enforce a patent
+(such as an express permission to practice a patent or covenant not to
+sue for patent infringement).  To "grant" such a patent license to a
+party means to make such an agreement or commitment not to enforce a
+patent against the party.
+
+  If you convey a covered work, knowingly relying on a patent license,
+and the Corresponding Source of the work is not available for anyone
+to copy, free of charge and under the terms of this License, through a
+publicly available network server or other readily accessible means,
+then you must either (1) cause the Corresponding Source to be so
+available, or (2) arrange to deprive yourself of the benefit of the
+patent license for this particular work, or (3) arrange, in a manner
+consistent with the requirements of this License, to extend the patent
+license to downstream recipients.  "Knowingly relying" means you have
+actual knowledge that, but for the patent license, your conveying the
+covered work in a country, or your recipient's use of the covered work
+in a country, would infringe one or more identifiable patents in that
+country that you have reason to believe are valid.
+
+  If, pursuant to or in connection with a single transaction or
+arrangement, you convey, or propagate by procuring conveyance of, a
+covered work, and grant a patent license to some of the parties
+receiving the covered work authorizing them to use, propagate, modify
+or convey a specific copy of the covered work, then the patent license
+you grant is automatically extended to all recipients of the covered
+work and works based on it.
+
+  A patent license is "discriminatory" if it does not include within
+the scope of its coverage, prohibits the exercise of, or is
+conditioned on the non-exercise of one or more of the rights that are
+specifically granted under this License.  You may not convey a covered
+work if you are a party to an arrangement with a third party that is
+in the business of distributing software, under which you make payment
+to the third party based on the extent of your activity of conveying
+the work, and under which the third party grants, to any of the
+parties who would receive the covered work from you, a discriminatory
+patent license (a) in connection with copies of the covered work
+conveyed by you (or copies made from those copies), or (b) primarily
+for and in connection with specific products or compilations that
+contain the covered work, unless you entered into that arrangement,
+or that patent license was granted, prior to 28 March 2007.
+
+  Nothing in this License shall be construed as excluding or limiting
+any implied license or other defenses to infringement that may
+otherwise be available to you under applicable patent law.
+
+  12. No Surrender of Others' Freedom.
+
+  If conditions are imposed on you (whether by court order, agreement or
+otherwise) that contradict the conditions of this License, they do not
+excuse you from the conditions of this License.  If you cannot convey a
+covered work so as to satisfy simultaneously your obligations under this
+License and any other pertinent obligations, then as a consequence you may
+not convey it at all.  For example, if you agree to terms that obligate you
+to collect a royalty for further conveying from those to whom you convey
+the Program, the only way you could satisfy both those terms and this
+License would be to refrain entirely from conveying the Program.
+
+  13. Use with the GNU Affero General Public License.
+
+  Notwithstanding any other provision of this License, you have
+permission to link or combine any covered work with a work licensed
+under version 3 of the GNU Affero General Public License into a single
+combined work, and to convey the resulting work.  The terms of this
+License will continue to apply to the part which is the covered work,
+but the special requirements of the GNU Affero General Public License,
+section 13, concerning interaction through a network will apply to the
+combination as such.
+
+  14. Revised Versions of this License.
+
+  The Free Software Foundation may publish revised and/or new versions of
+the GNU General Public License from time to time.  Such new versions will
+be similar in spirit to the present version, but may differ in detail to
+address new problems or concerns.
+
+  Each version is given a distinguishing version number.  If the
+Program specifies that a certain numbered version of the GNU General
+Public License "or any later version" applies to it, you have the
+option of following the terms and conditions either of that numbered
+version or of any later version published by the Free Software
+Foundation.  If the Program does not specify a version number of the
+GNU General Public License, you may choose any version ever published
+by the Free Software Foundation.
+
+  If the Program specifies that a proxy can decide which future
+versions of the GNU General Public License can be used, that proxy's
+public statement of acceptance of a version permanently authorizes you
+to choose that version for the Program.
+
+  Later license versions may give you additional or different
+permissions.  However, no additional obligations are imposed on any
+author or copyright holder as a result of your choosing to follow a
+later version.
+
+  15. Disclaimer of Warranty.
+
+  THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY
+APPLICABLE LAW.  EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT
+HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY
+OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO,
+THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+PURPOSE.  THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM
+IS WITH YOU.  SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF
+ALL NECESSARY SERVICING, REPAIR OR CORRECTION.
+
+  16. Limitation of Liability.
+
+  IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS
+THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY
+GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE
+USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF
+DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD
+PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS),
+EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF
+SUCH DAMAGES.
+
+  17. Interpretation of Sections 15 and 16.
+
+  If the disclaimer of warranty and limitation of liability provided
+above cannot be given local legal effect according to their terms,
+reviewing courts shall apply local law that most closely approximates
+an absolute waiver of all civil liability in connection with the
+Program, unless a warranty or assumption of liability accompanies a
+copy of the Program in return for a fee.
+
+                     END OF TERMS AND CONDITIONS
+
+            How to Apply These Terms to Your New Programs
+
+  If you develop a new program, and you want it to be of the greatest
+possible use to the public, the best way to achieve this is to make it
+free software which everyone can redistribute and change under these terms.
+
+  To do so, attach the following notices to the program.  It is safest
+to attach them to the start of each source file to most effectively
+state the exclusion of warranty; and each file should have at least
+the "copyright" line and a pointer to where the full notice is found.
+
+    
+    Copyright (C)   
+
+    This program 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.
+
+    This program 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 this program.  If not, see .
+
+Also add information on how to contact you by electronic and paper mail.
+
+  If the program does terminal interaction, make it output a short
+notice like this when it starts in an interactive mode:
+
+      Copyright (C)   
+    This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
+    This is free software, and you are welcome to redistribute it
+    under certain conditions; type `show c' for details.
+
+The hypothetical commands `show w' and `show c' should show the appropriate
+parts of the General Public License.  Of course, your program's commands
+might be different; for a GUI interface, you would use an "about box".
+
+  You should also get your employer (if you work as a programmer) or school,
+if any, to sign a "copyright disclaimer" for the program, if necessary.
+For more information on this, and how to apply and follow the GNU GPL, see
+.
+
+  The GNU General Public License does not permit incorporating your program
+into proprietary programs.  If your program is a subroutine library, you
+may consider it more useful to permit linking proprietary applications with
+the library.  If this is what you want to do, use the GNU Lesser General
+Public License instead of this License.  But first, please read
+.