Merge pull request #1 from KULL-Centre/main

Update from head

.github/workflows/python-app.yml

@@ -0,0 +1,40 @@
+# This workflow will install Python dependencies, run tests and lint with a single version of Python
+# For more information see: https://docs.github.com/en/actions/automating-builds-and-tests/building-and-testing-python
+
+name: Python application
+
+on:
+  push:
+    branches: [ "main" ]
+  pull_request:
+    branches: [ "main" ]
+
+permissions:
+  contents: read
+
+jobs:
+  build:
+
+    runs-on: ubuntu-latest
+
+    steps:
+    - uses: actions/checkout@v4
+    - name: Set up Python 3.10
+      uses: actions/setup-python@v3
+      with:
+        python-version: "3.10"
+    - name: Install pip dependencies
+      run: |
+        python -m pip install --upgrade pip
+        pip install flake8 pytest
+        pip install .
+        # if [ -f requirements.txt ]; then pip install -r requirements.txt; fi
+    - name: Lint with flake8
+      run: |
+        # stop the build if there are Python syntax errors or undefined names
+        flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics
+        # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide
+        flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics
+    - name: Test with pytest
+      run: |
+        pytest

.readthedocs.yaml

@@ -0,0 +1,35 @@
+# Read the Docs configuration file for Sphinx projects
+# See https://docs.readthedocs.io/en/stable/config-file/v2.html for details
+
+# Required
+version: 2
+
+# Set the OS, Python version and other tools you might need
+build:
+  os: ubuntu-22.04
+  tools:
+    python: "3.12"
+    # You can also specify other tool versions:
+    # nodejs: "20"
+    # rust: "1.70"
+    # golang: "1.20"
+
+# Build documentation in the "docs/" directory with Sphinx
+sphinx:
+  configuration: docs/conf.py
+  # You can configure Sphinx to use a different builder, for instance use the dirhtml builder for simpler URLs
+  # builder: "dirhtml"
+  # Fail on all warnings to avoid broken references
+  # fail_on_warning: true
+
+# Optionally build your docs in additional formats such as PDF and ePub
+# formats:
+#   - pdf
+#   - epub
+
+# Optional but recommended, declare the Python requirements required
+# to build your documentation
+# See https://docs.readthedocs.io/en/stable/guides/reproducible-builds.html
+python:
+  install:
+    - requirements: docs/requirements.txt

README.md

@@ -2,10 +2,11 @@
 [![Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/KULL-Centre/_2023_Tesei_IDRome/blob/main/IDRLab.ipynb)
 [![CALVADOS Video](http://img.shields.io/badge/►-Video-FF0000.svg)](https://youtu.be/r-eFzoBiQZ4)
 [![IDRome Video](http://img.shields.io/badge/►-Video-FF0000.svg)](https://youtu.be/kL3-cusHgzM)
+[![Python application](https://github.com/KULL-Centre/CALVADOS/actions/workflows/python-app.yml/badge.svg)](https://github.com/KULL-Centre/CALVADOS/actions/workflows/python-app.yml)
 
 # CALVADOS
 
-Coarse-grained implicit-solvent simulations of biomolecules in the openMM framework.
+Coarse-grained implicit-solvent simulations of biomolecules in the OpenMM framework.
 Earlier implementations of the code are available on [Zenodo](https://zenodo.org/search?q=metadata.subjects.subject%3A%22CALVADOS%22&l=list&p=1&s=10&sort=bestmatch) ([DOI: 10.5281/zenodo.13754000](https://doi.org/10.5281/zenodo.13754000)).
 
 Please cite the following references when using the software:
@@ -21,11 +22,9 @@ Please cite the following references when using the software:
 conda create -n calvados python=3.10
 conda activate calvados
 ```
-2. Install numba, mdtraj with conda and openmm (they have caused issues with pip install)
+(2. Only needed when planning to use GPUs: Install openmm via conda-force with cudatoolkit. This step can be skipped if running on CPU only.)
 ```
-conda install numba
-conda install -c conda-forge mdtraj
-conda install -c conda-forge openmm cudatoolkit=11.2
+conda install -c conda-forge openmm=8.2.0 cudatoolkit=11.8
 ```
 3. Clone package and install CALVADOS and its dependencies using pip
 ``` 
@@ -34,10 +33,14 @@ cd CALVADOS
 pip install .
 (or pip install -e .)
 ```
-4. Clean up faulty pip install of scipy:
-```
-conda install scipy
+
+## Testing
+
+```bash
+
+  python -m pytest
 ```
+The test `test_ah_dh_interactions` simulates two free amino acids, calculates the potential energies based on the saved trajectory and compares these values with those in the OpenMM log file.
 
 ## Contact
 
@@ -52,5 +55,7 @@ For further questions and inquiries, please contact us.
 
 [Fan Cao (@fancaoErik)](https://github.com/fancaoErik)
 
+[Ikki Yasuda (@iyasuda)](https://github.com/iyasuda)
+
 [Kresten Lindorff-Larsen (@lindorff-larsen)](https://github.com/lindorff-larsen)
 

calvados/__init__.py

@@ -5,8 +5,9 @@
 import calvados.postprocess
 import calvados.sequence
 import calvados.sim
-import calvados.slab_conc
+# import calvados.slab_conc
 import calvados.utilities
 import calvados.components
 
-
+# from pathlib import Path
+# PACKAGEDIR = Path(__file__).parent.absolute()

calvados/analysis.py

@@ -94,15 +94,21 @@ def calc_raw_dmap(pos0,pos1):
     dmap = distances.distance_array(pos0,pos1)
     return dmap
 
-def self_distances(pos):
+def self_distances(pos,box=None):
     """ Self distance map for matrix of positions
 
-    Input: Matrix of positions
+    If box dimensions are provided, distances are
+    calculated using minimum image convention
+
+    Input: Matrix of positions and (optional) box dimensions
     Output: Self distance map
     """
     N = len(pos)
     dmap = np.zeros((N,N))
-    d = distances.self_distance_array(pos)
+    if box is not None:
+        d = distances.self_distance_array(pos,box)
+    else:
+        d = distances.self_distance_array(pos)
     k = 0
     for i in range(N):
         for j in range(i + 1, N):
@@ -468,7 +474,7 @@ def center_slab(path,name,ref_atoms='all',start=None,end=None,step=1,input_pdb='
             W.write(ag)
     return hs, z/10.
 
-def calc_slab_profiles(path,name,ref_atoms,sel_atoms_list,output_folder,start=None,end=None,step=1,input_pdb='top.pdb'):
+def calc_slab_profiles(path,name,output_folder,ref_atoms,sel_atoms_list=[],start=None,end=None,step=1,input_pdb='top.pdb'):
     """
     path: path where trajectory and pdb are saved
     ref_atoms: reference atoms to shift to the middle of the box
@@ -487,15 +493,16 @@ def calc_slab_profiles(path,name,ref_atoms,sel_atoms_list,output_folder,start=No
     volume = u.dimensions[0]*u.dimensions[1]*binwidth/1e3
     # density profile for ref_atoms
     edges = np.arange(0,z.size+binwidth,binwidth)
+    n_bins = edges.size - 1
     np.save(output_folder+f'/{name:s}_ref_profile.npy',np.c_[h_ref/volume.mean()])
     h_ref_mean = h_ref.mean(axis=0)/volume.mean() # number of beads per nm3
-    n_bins = edges.size - 1
     all_profiles = np.c_[z,h_ref_mean]
     # density profile: selected atoms
     for i,sel_atoms in enumerate(sel_atoms_list):
         ag_sel = u.select_atoms(sel_atoms)
         h_sel = np.zeros((n_frames,n_bins))
         for t,ts in enumerate(u.trajectory[start:end:step]):
+            ts = transformations.wrap(ag_sel)(ts)
             zpos = ag_sel.positions.T[2]
             h, e = np.histogram(zpos,bins=edges)
             h_sel[t] = h

calvados/build.py

@@ -81,17 +81,21 @@ def draw_starting_vec(box):
     vec = np.random.random(size=3) * box
     return vec
 
-def build_linear(bondlengths):
+def build_linear(z_bondlengths, n_per_res=1, ys=None):
     """
-    linear chain growing in z direction
+    linear chain growing in z direction with possible extra beads at y offset n_per_res*ys
     centered at [0,0,0]
     """
-    z = 0
-    N = len(bondlengths)
-    coords = np.zeros((N,3))
-    for idx in range(N-1):
-        z += 0.5 * (bondlengths[idx] + bondlengths[idx+1])
-        coords[idx+1,2] = z
+    N = len(z_bondlengths)
+    zs = np.zeros(N)
+    zs[1:] = [(z_bondlengths[idx] + z_bondlengths[idx+1]) / 2 for idx in range(N-1)]
+    zs = np.cumsum(zs)
+    coords = np.zeros((N*n_per_res,3))
+    if ys == None:
+        ys = np.zeros(N)
+    for idx, (z, y) in enumerate(zip(zs,ys)):
+        for jdx in range(n_per_res):
+            coords[idx*n_per_res + jdx] = [0,jdx*y,z]
     coords[:,2] -= coords[:,2].mean()
     return coords
 
@@ -101,7 +105,7 @@ def p2c(r, phi):
     """
     return (r * np.cos(phi), r * np.sin(phi))
 
-def build_spiral(n, delta=0, arc=.38, separation=.7):
+def build_spiral(bondlengths, delta=0, arc=.38, separation=.7, n_per_res=1):
     """
     create points on an Archimedes' spiral
     with `arc` giving the length of arc between two points
@@ -112,12 +116,41 @@ def build_spiral(n, delta=0, arc=.38, separation=.7):
     b = separation / (2 * np.pi)
     phi = float(r) / b
     coords = []
-    for i in range(n):
-        coords.append(list(p2c(r, phi))+[0])
+    for i,z in enumerate(bondlengths):
+        for j in range(n_per_res): # number of beads per residue (placed along z with bondlength)
+            coords.append(list(p2c(r, phi))+[j*z]) # j*z = 0 for n_per_res=1
         phi += float(arc) / r
         r = b * phi
     return np.array(coords)+delta
 
+def build_compact(nbeads, d=0.38, verbose=False):
+    N = int(np.ceil(np.cbrt(nbeads)) - 1)
+    if verbose:
+        print(f'Building {N+1} * {N+1} * {N+1} grid.')
+    xs = []
+    i, j, k = 0, 0, 0
+    di, dj, dk = 1, 1, 1 # direction
+    cti, ctj, ctk = 0, 0, 0
+
+    for idx in range(nbeads):
+        xs.append([i,j,k])
+        if ctk == N:
+            if ctj == N:
+                i += di
+                cti += 1
+                ctj = 0
+                dj *= -1
+            else:
+                j += dj
+                ctj += 1
+            ctk = 0
+            dk *= -1
+        else:
+            k += dk
+            ctk += 1
+    xs = (np.array(xs) - 0.5*N) * d
+    return xs
+
 def random_placement(box,xs_others,xinit,ntries=10000):
     ntry = 0
     while True: # random placement
@@ -277,7 +310,7 @@ def geometry_from_pdb(pdb,use_com=False):
     else:
         cas = u.select_atoms('name CA')
         pos = cas.positions / 10.
-    return pos
+    return pos, u.dimensions
 
 def bfac_from_pdb(pdb,confidence=70.):
     """ get pLDDT encoded in pdb b-factor column """

calvados/cfg.py

@@ -34,7 +34,7 @@ def load_default_config(pkg_base):
         """ Load default config. """
         with open(f'{pkg_base}/data/default_config.yaml','r') as stream:
             default_config = yaml.safe_load(stream)
-        default_config['fresidues'] = f'{pkg_base}/data/residues.csv'
+        # default_config['fresidues'] = f'{pkg_base}/data/residues.csv'
         return default_config
 
     @staticmethod