Enable PDBx/mmCIF support

.github/workflows/python-app.yml

@@ -19,10 +19,10 @@ jobs:
 
     steps:
     - uses: actions/checkout@v4
-    - name: Set up Python 3.10
+    - name: Set up Python 3.13
       uses: actions/setup-python@v3
       with:
-        python-version: "3.10"
+        python-version: "3.13"
     - name: Install pip dependencies
       run: |
         python -m pip install --upgrade pip

README.md

@@ -26,14 +26,18 @@ The examples described in the paper can be found in the `examples` folder.
 
 ## Installation Instructions
 
-1. Make new conda environment for calvados
+1. Create a new conda environment for calvados
 ``` 
-conda create -n calvados python=3.10
+conda create -n calvados python=3.13
 conda activate calvados
 ```
-(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.)
+2. Install openmm via conda-force with cudatoolkit, if you are planning to run on GPUs with CUDA
 ```
-conda install -c conda-forge openmm=8.2.0 cudatoolkit=11.8
+conda install -c conda-forge openmm=8.2.0 cudatoolkit=11.8 mdanalysis=2.9 mdtraj=1.11
+```
+or without cudatoolkit if you are running on CPU or OpenCL
+```
+conda install -c conda-forge openmm=8.2.0 mdanalysis=2.9 mdtraj=1.11
 ```
 3. Clone package and install CALVADOS and its dependencies using pip
 ``` 

calvados/analysis.py

@@ -28,9 +28,7 @@
 
 def center_traj(pdb,traj,start=None,stop=None,step=1):
     """ Center trajectory """
-
     u = mda.Universe(pdb,traj)
-
     with mda.Writer(f'{traj[:-4]}_c.dcd', len(u.atoms)) as W:
         for ts in u.trajectory[start:stop:step]:
             u.atoms.translate(-u.atoms.center_of_geometry() + 0.5 * u.dimensions[:3])
@@ -965,5 +963,6 @@ def calc_contact_map(path,sysname,output_path,chainid_dict={},is_slab=False,inpu
             cmap += (.5-.5*np.tanh((d-1.)/.3)).reshape(mask_frames.sum(),
                         N_res_1,-1,N_res_2).sum(axis=0).sum(axis=1)
     cmap /= traj.n_frames
+
     # save energy and contact maps
     np.save(output_path+f'/{sysname:s}_{name_1:s}_{name_2:s}_cmap.npy',cmap)

calvados/build.py

@@ -1,7 +1,7 @@
 from scipy import constants
 import numpy as np
 import pandas as pd
-from openmm import unit
+from openmm import app, unit
 
 from MDAnalysis import Universe
 from MDAnalysis.analysis import distances
@@ -223,17 +223,17 @@ def build_xyzgrid(N,box):
     """ 3D grid """
 
     r = box / np.sum(box)
-    a = np.cbrt(N / np.product(r))
+    a = np.cbrt(N / np.prod(r))
     n = a * r
     nxyz = np.floor(n)
-    while np.product(nxyz) < N:
+    while np.prod(nxyz) < N:
         ndeviation = n / nxyz
         devmax = np.argmax(ndeviation)
         nxyz[devmax] += 1
-    while np.product(nxyz) > N:
+    while np.prod(nxyz) > N:
         nmax = np.argmax(nxyz)
         nxyz[nmax] -= 1
-        if np.product(nxyz) < N:
+        if np.prod(nxyz) < N:
             nxyz[nmax] += 1
             break
 
@@ -298,7 +298,11 @@ def geometry_from_pdb(pdb,use_com=False):
     """ positions in nm"""
     with catch_warnings():
         simplefilter("ignore")
-        u = Universe(pdb)
+        if pdb.lower().endswith('.cif'):
+            pdbx = app.pdbxfile.PDBxFile(pdb)
+            u = Universe(pdbx)
+        else:
+            u = Universe(pdb)
     ag = u.atoms
     ag.translate(-ag.center_of_mass())
     if use_com:

calvados/components.py

@@ -8,6 +8,8 @@
 
 import numpy as np
 
+import os
+
 class Component:
     """ Generic component of the system. """
 
@@ -33,7 +35,12 @@ def calc_comp_seq(self):
         """ Calculate sequence of component. """
 
         if self.restraint:
-            self.seq, self.n_termini, self.c_termini = seq_from_pdb(f'{self.pdb_folder}/{self.name}.pdb')
+            input_pdb = f'{self.pdb_folder}/{self.name}.pdb'
+            input_pdbx = f'{self.pdb_folder}/{self.name}.cif'
+            if os.path.isfile(input_pdbx):
+                self.seq, self.n_termini, self.c_termini = seq_from_pdb(input_pdbx)
+            else:
+                self.seq, self.n_termini, self.c_termini = seq_from_pdb(input_pdb)
         else:
             records = read_fasta(self.ffasta)
             self.seq = str(records[self.name].seq)
@@ -116,7 +123,12 @@ def calc_x_from_pdb(self):
         """ Calculate protein positions from pdb. """
 
         input_pdb = f'{self.pdb_folder}/{self.name}.pdb'
-        self.xinit, self.dimensions = build.geometry_from_pdb(input_pdb,use_com=self.use_com) # read from pdb
+        input_pdbx = f'{self.pdb_folder}/{self.name}.cif'
+
+        if os.path.isfile(input_pdbx):
+            self.xinit, self.dimensions = build.geometry_from_pdb(input_pdbx,use_com=self.use_com) # read from pdb
+        else:
+            self.xinit, self.dimensions = build.geometry_from_pdb(input_pdb,use_com=self.use_com) # read from pdbx
 
     def calc_ssdomains(self):
         """ Get bounds for restraints (harmonic). """

calvados/sequence.py

@@ -26,6 +26,8 @@
 
 from calvados import analysis, interactions
 
+from openmm import app
+
 ### SEQUENCE INPUT / OUTPUT
 def read_fasta(ffasta):
     records = SeqIO.to_dict(SeqIO.parse(ffasta, "fasta"))
@@ -36,7 +38,11 @@ def seq_from_pdb(pdb,selection='all',fmt='string'):
 
     with warnings.catch_warnings():
         warnings.simplefilter("ignore")
-        u = Universe(pdb)
+        if pdb.lower().endswith('.cif'):
+            pdbx = app.pdbxfile.PDBxFile(pdb)
+            u = Universe(pdbx)
+        else:
+            u = Universe(pdb)
 
     # we do not assume residues in the PDB file are numbered from 1
     n_termini = [0]

calvados/sim.py

@@ -215,9 +215,11 @@ def build_system(self):
             self.add_custom_restraints()
 
         self.pdb_cg = f'{self.path}/top.pdb'
+        self.cif_cg = f'{self.path}/top.cif'
         a = md.Trajectory(self.pos, self.top, 0, self.box, [90,90,90])
         if self.restart != 'pdb': # only save new topology if no system pdb is given
             a.save_pdb(self.pdb_cg)
+            a.save_cif(self.cif_cg)
 
         self.add_forces_to_system()
         self.print_system_summary()
@@ -507,10 +509,23 @@ def simulate(self):
         fcheck_out = f'{self.path}/restart.chk'
         append = False
 
-        if self.restart == 'pdb' and os.path.isfile(fcheck_in):
-            pdb = app.pdbfile.PDBFile(fcheck_in)
+        if self.restart == 'pdb':
+            if os.path.isfile(fcheck_in):
+                pdb = app.pdbfile.PDBFile(fcheck_in)
+            else:
+                pdb = app.pdbfile.PDBFile(self.pdb_cg)
+        elif self.restart == 'cif':
+            if os.path.isfile(fcheck_in):
+                pdb = app.pdbxfile.PDBxFile(fcheck_in)
+            else:
+                pdb = app.pdbxfile.PDBxFile(self.cif_cg)
+        elif self.restart == 'checkpoint':
+            if os.path.isfile(fcheck_in):
+                pdb = app.pdbxfile.PDBxFile(fcheck_in)
+            else:
+                pdb = app.pdbxfile.PDBxFile(self.cif_cg)
         else:
-            pdb = app.pdbfile.PDBFile(self.pdb_cg)
+            pdb = app.pdbxfile.PDBxFile(self.cif_cg)
 
         # use langevin integrator
         integrator = openmm.openmm.LangevinMiddleIntegrator(self.temp*unit.kelvin,self.friction_coeff/unit.picosecond,0.01*unit.picosecond)
@@ -537,7 +552,7 @@ def simulate(self):
             print(f'Appending log file to {self.path}/{self.sysname:s}.log')
             simulation.loadCheckpoint(fcheck_in)
         else:
-            if self.restart == 'pdb':
+            if self.restart in ['pdb','cif']:
                 print(f'Reading in system configuration {self.frestart}')
             elif self.restart == 'checkpoint':
                 print(f'No checkpoint file {self.frestart} found: Starting from new system configuration')
@@ -561,9 +576,11 @@ def simulate(self):
             simulation.reporters.append(app.dcdreporter.DCDReporter(f'{self.path}/equilibration_{self.sysname:s}.dcd',self.wfreq,append=append))
             simulation.step(self.steps_eq)
             state_final = simulation.context.getState(getPositions=True)
-            rep = app.pdbreporter.PDBReporter(f'{self.path}/equilibration_final.pdb',0)
-            rep.report(simulation,state_final)
-            pdb = app.pdbfile.PDBFile(f'{self.path}/equilibration_final.pdb')
+            with open(f'{self.path}/equilibration_final.pdb', 'w') as f:
+                app.PDBFile.writeFile(simulation.topology, state_final.getPositions(), f)
+            with open(f'{self.path}/equilibration_final.cif', 'w') as f:
+                app.PDBxFile.writeFile(simulation.topology, state_final.getPositions(), f)
+            pdb = app.pdbxfile.PDBxFile(f'{self.path}/equilibration_final.cif')
 
             for index, force in enumerate(self.system.getForces()):
                 if isinstance(force, openmm.CustomExternalForce):
@@ -584,9 +601,11 @@ def simulate(self):
             simulation.reporters.append(app.dcdreporter.DCDReporter(f'{self.path}/equilibration_{self.sysname:s}.dcd',self.wfreq,append=append))
             simulation.step(self.steps_eq)
             state_final = simulation.context.getState(getPositions=True,enforcePeriodicBox=True)
-            rep = app.pdbreporter.PDBReporter(f'{self.path}/equilibration_final.pdb',0)
-            rep.report(simulation,state_final)
-            pdb = app.pdbfile.PDBFile(f'{self.path}/equilibration_final.pdb')
+            with open(f'{self.path}/equilibration_final.pdb', 'w') as f:
+                app.PDBFile.writeFile(simulation.topology, state_final.getPositions(), f)
+            with open(f'{self.path}/equilibration_final.cif', 'w') as f:
+                app.PDBxFile.writeFile(simulation.topology, state_final.getPositions(), f)
+            pdb = app.pdbxfile.PDBxFile(f'{self.path}/equilibration_final.cif')
             topology = pdb.getTopology()
             a, b, c = state_final.getPeriodicBoxVectors()
             topology.setPeriodicBoxVectors(state_final.getPeriodicBoxVectors())
@@ -632,10 +651,10 @@ def simulate(self):
         dt_string = now.strftime("%Y%d%m_%Hh%Mm%Ss")
 
         state_final = simulation.context.getState(getPositions=True,enforcePeriodicBox=True)
-        rep = app.pdbreporter.PDBReporter(f'{self.path}/{self.sysname}_{dt_string}.pdb',0)
-        rep.report(simulation,state_final)
-        rep = app.pdbreporter.PDBReporter(f'{self.path}/checkpoint.pdb',0)
-        rep.report(simulation,state_final)
+        with open(f'{self.path}/checkpoint.pdb', 'w') as f:
+            app.PDBFile.writeFile(simulation.topology, state_final.getPositions(), f)
+        with open(f'{self.path}/checkpoint.cif', 'w') as f:
+            app.PDBxFile.writeFile(simulation.topology, state_final.getPositions(), f)
 
 def run(path='.',fconfig='config.yaml',fcomponents='components.yaml'):
     with open(f'{path}/{fconfig}','r') as stream:

examples/slab_IDR/README.md

@@ -1,7 +1,7 @@
 The lines below run direct coexistence simulations of 100 copies of a single IDR:
 
 ```bash
-python prepare.py --name <protein name>  --replica <replica>
+python prepare.py --name <protein name>  --replica <replica> --gpu_id <gpu_id>
 python <protein name>_<replica>/run.py --path <protein name>_<replica>
 ```
 

examples/slab_IDR/prepare.py

@@ -34,7 +34,7 @@
   wfreq = N_save, # dcd writing frequency, 1 = 10fs
   steps = N_frames*N_save, # number of simulation steps
   runtime = 0, # overwrites 'steps' keyword if > 0
-  platform = 'CUDA', # 'CUDA'
+  platform = 'CPU', # 'CUDA'
   restart = 'checkpoint',
   frestart = 'restart.chk',
   verbose = True,

examples/slab_IDR_MDP/prepare.py

@@ -68,6 +68,10 @@
 # calculate contact map between name_1 and name_1
 chainid_dict = dict({args.name_1:s} = (0,99))
 calc_contact_map(path="{path:s}",sysname="{sysname:s}",output_path="data",chainid_dict=chainid_dict,is_slab=True)
+
+# calculate contact map between name_2 and name_2 (note that this is not always meaningful)
+# chainid_dict = dict({args.name_2:s} = (100,199))
+# calc_contact_map(path="{path:s}",sysname="{sysname:s}",output_path="data",chainid_dict=chainid_dict,is_slab=True)
 """
 
 config.write(path,name='config.yaml',analyses=analyses)

setup.py

@@ -11,7 +11,7 @@
 
 setup(
     name='calvados',
-    version='0.8.0',
+    version='0.8.1',
     description='Coarse-grained implicit-solvent simulations of biomolecules',
     url='https://github.com/KULL-Centre/CALVADOS',
     authors=[
@@ -23,29 +23,30 @@
     license='GNU GPL3',
     packages=find_packages(),
     install_requires=[
-        'OpenMM==8.2.0',
-        'numpy==1.24',
-        'pandas==2.1.1',
-        'MDAnalysis==2.6.1',
-        'biopython==1.81',
-        'Jinja2==3.1.2',
-        'tqdm==4.66',
-        'matplotlib==3.8',
-        'PyYAML==6.0',
-        'statsmodels==0.14',
-        'localcider==0.1.21',
-        'pytest==8.3.3',
-        'numba==0.60',
-        'scipy==1.13',
-        'mdtraj==1.10',
+        "OpenMM>=8.4,<8.5",
+        "MDAnalysis>=2.10,<2.11",
+        "mdtraj>=1.11,<1.12",
+        'numpy',
+        'pandas',
+        'biopython',
+        'Jinja2',
+        'tqdm',
+        'matplotlib',
+        'PyYAML',
+        'statsmodels',
+        'localcider',
+        'pytest',
+        'numba',
+        'scipy'
     ],
+    python_requires=">=3.13,<3.14",
 
     # include_package_data=True,
     package_data={'' : ['data/*.csv', 'data/*.yaml', 'data/templates/*']},
 
     classifiers=[
         'Intended Audience :: Science/Research',
         'Operating System :: POSIX :: Linux',
-        'Programming Language :: Python :: 3>=3.7,<3.11',
+        'Programming Language :: Python :: 3.13',
     ],
 )