From 49524a089c2786720ec46782993601152c5fe8ba Mon Sep 17 00:00:00 2001 From: gitesei Date: Thu, 12 Feb 2026 11:20:15 +0100 Subject: [PATCH 1/7] Enable PDBx/mmCIF support --- README.md | 6 +++--- calvados/build.py | 8 ++++++-- calvados/components.py | 9 ++++++++- calvados/sim.py | 42 +++++++++++++++++++++++++++--------------- setup.py | 27 ++++++++++++--------------- 5 files changed, 56 insertions(+), 36 deletions(-) diff --git a/README.md b/README.md index 675f4a1..d37579b 100644 --- a/README.md +++ b/README.md @@ -28,12 +28,12 @@ The examples described in the paper can be found in the `examples` folder. 1. Make 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) ``` -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 ``` 3. Clone package and install CALVADOS and its dependencies using pip ``` diff --git a/calvados/build.py b/calvados/build.py index 2f9311b..cd6a96c 100644 --- a/calvados/build.py +++ b/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 @@ -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: diff --git a/calvados/components.py b/calvados/components.py index 94f4362..561dc2f 100644 --- a/calvados/components.py +++ b/calvados/components.py @@ -8,6 +8,8 @@ import numpy as np +import os + class Component: """ Generic component of the system. """ @@ -116,7 +118,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_pdb): + self.xinit, self.dimensions = build.geometry_from_pdb(input_pdb,use_com=self.use_com) # read from pdb + else: + self.vxinit, self.dimensions = build.geometry_from_pdb(input_pdbx,use_com=self.use_com) # read from pdbx def calc_ssdomains(self): """ Get bounds for restraints (harmonic). """ diff --git a/calvados/sim.py b/calvados/sim.py index bebc9ef..f974800 100644 --- a/calvados/sim.py +++ b/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,16 @@ 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) - else: - pdb = app.pdbfile.PDBFile(self.pdb_cg) + 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) # use langevin integrator integrator = openmm.openmm.LangevinMiddleIntegrator(self.temp*unit.kelvin,self.friction_coeff/unit.picosecond,0.01*unit.picosecond) @@ -537,7 +545,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 +569,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 +594,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 +644,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: diff --git a/setup.py b/setup.py index 1adde58..4979576 100644 --- a/setup.py +++ b/setup.py @@ -23,21 +23,18 @@ 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', + 'numpy', + 'pandas', + 'biopython', + 'Jinja2', + 'tqdm', + 'matplotlib', + 'PyYAML', + 'statsmodels', + 'localcider', + 'pytest', + 'numba', + 'scipy' ], # include_package_data=True, From 6426099ed55c7e71d01bff702518196693c51e7f Mon Sep 17 00:00:00 2001 From: gitesei Date: Wed, 18 Feb 2026 12:13:01 +0100 Subject: [PATCH 2/7] Fix deprecated numpy product and restart options in simulate --- README.md | 9 +++++++-- calvados/analysis.py | 12 +++++++++--- calvados/build.py | 8 ++++---- calvados/sim.py | 5 +++++ examples/slab_IDR/README.md | 2 +- examples/slab_IDR/prepare.py | 2 +- 6 files changed, 27 insertions(+), 11 deletions(-) diff --git a/README.md b/README.md index d37579b..15b6170 100644 --- a/README.md +++ b/README.md @@ -26,19 +26,24 @@ 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.13 conda activate calvados ``` -2. Install openmm via conda-force (with cudatoolkit if you are planning to run on GPUs) +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 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 ``` git clone https://github.com/KULL-Centre/CALVADOS.git cd CALVADOS +git checkout pdbx pip install . (or pip install -e .) ``` diff --git a/calvados/analysis.py b/calvados/analysis.py index 8a7df19..ce807f5 100644 --- a/calvados/analysis.py +++ b/calvados/analysis.py @@ -28,9 +28,9 @@ 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]) @@ -922,6 +922,11 @@ def calc_contact_map(path,sysname,output_path,chainid_dict={},is_slab=False,inpu name_2 = name_1 chainid_dict[name_1] = np.arange(traj.top.n_chains) + print(name_1) + print(chainid_dict[name_1]) + print(name_2) + print(chainid_dict[name_2]) + N_res_1 = traj.top.chain(chainid_dict[name_1][0]).n_residues N_res_2 = traj.top.chain(chainid_dict[name_2][0]).n_residues @@ -950,7 +955,8 @@ def calc_contact_map(path,sysname,output_path,chainid_dict={},is_slab=False,inpu chainid_dict[name_2] = chainid_dict[name_1] cm_z = cmtraj.xyz[:,chainid_dict[name_1],2] # per-frame central-chain indices - chainid_dict[name_1] = np.argmin(np.abs(cm_z),axis=1) + ids_central = np.argmin(np.abs(cm_z),axis=1) + chainid_dict[name_1] = np.array([chainid_dict[name_1][idx] for idx in ids_central]) cmap = np.zeros((N_res_1,N_res_2)) for chain_1 in np.unique(chainid_dict[name_1]): diff --git a/calvados/build.py b/calvados/build.py index cd6a96c..10e21e0 100644 --- a/calvados/build.py +++ b/calvados/build.py @@ -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 diff --git a/calvados/sim.py b/calvados/sim.py index f974800..0b84573 100644 --- a/calvados/sim.py +++ b/calvados/sim.py @@ -519,6 +519,11 @@ def simulate(self): 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) # use langevin integrator integrator = openmm.openmm.LangevinMiddleIntegrator(self.temp*unit.kelvin,self.friction_coeff/unit.picosecond,0.01*unit.picosecond) diff --git a/examples/slab_IDR/README.md b/examples/slab_IDR/README.md index a867a5b..e5b9fd2 100644 --- a/examples/slab_IDR/README.md +++ b/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 --replica +python prepare.py --name --replica --gpu_id python _/run.py --path _ ``` diff --git a/examples/slab_IDR/prepare.py b/examples/slab_IDR/prepare.py index 68d395a..5ee71cf 100644 --- a/examples/slab_IDR/prepare.py +++ b/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, From 646dfc09875fca8e728d52cf4359a6160634d0e9 Mon Sep 17 00:00:00 2001 From: gitesei Date: Wed, 18 Feb 2026 22:09:04 +0100 Subject: [PATCH 3/7] Fix PDBx support --- calvados/components.py | 13 +++++++++---- calvados/sequence.py | 8 +++++++- 2 files changed, 16 insertions(+), 5 deletions(-) diff --git a/calvados/components.py b/calvados/components.py index 561dc2f..af6ad58 100644 --- a/calvados/components.py +++ b/calvados/components.py @@ -35,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) @@ -120,10 +125,10 @@ def calc_x_from_pdb(self): input_pdb = f'{self.pdb_folder}/{self.name}.pdb' input_pdbx = f'{self.pdb_folder}/{self.name}.cif' - if os.path.isfile(input_pdb): - self.xinit, self.dimensions = build.geometry_from_pdb(input_pdb,use_com=self.use_com) # read from pdb + 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.vxinit, self.dimensions = build.geometry_from_pdb(input_pdbx,use_com=self.use_com) # read from pdbx + 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). """ diff --git a/calvados/sequence.py b/calvados/sequence.py index dc1c389..80aeffa 100644 --- a/calvados/sequence.py +++ b/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] From c9558fade1e4b9adbfb63ecfc85eb6fc23397a8e Mon Sep 17 00:00:00 2001 From: gitesei Date: Tue, 24 Feb 2026 23:31:05 -0800 Subject: [PATCH 4/7] Fix topology input selection --- calvados/analysis.py | 15 +++++++++------ calvados/sim.py | 2 ++ 2 files changed, 11 insertions(+), 6 deletions(-) diff --git a/calvados/analysis.py b/calvados/analysis.py index ce807f5..1dd8a58 100644 --- a/calvados/analysis.py +++ b/calvados/analysis.py @@ -963,13 +963,16 @@ def calc_contact_map(path,sysname,output_path,chainid_dict={},is_slab=False,inpu surrounding_chains = traj.top.select(' or '.join([f'chainid {i:d}' for i in chainid_dict[name_2] if i != chain_1])) pair_indices = traj.top.select_pairs(f'chainid {chain_1:d}',surrounding_chains) if is_slab: - mask_frames = chainid_dict[name_1] == chain_1 + mask_frames = np.where(chainid_dict[name_1] == chain_1)[0] else: - mask_frames = np.full(traj.n_frames, True, dtype=bool) - if np.any(mask_frames): - d = md.compute_distances(traj[mask_frames],pair_indices) - 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) + mask_frames = np.arange(traj.n_frames) + if len(mask_frames) > 0: + for mf in mask_frames: + d = md.compute_distances(traj[mf],pair_indices)[0] + cm = (.5-.5*np.tanh((d-1.)/.3)).reshape(N_res_1,-1,N_res_2) + cm = np.sum(cm,axis=1) + cmap += cm + 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) diff --git a/calvados/sim.py b/calvados/sim.py index 0b84573..110d773 100644 --- a/calvados/sim.py +++ b/calvados/sim.py @@ -524,6 +524,8 @@ def simulate(self): pdb = app.pdbxfile.PDBxFile(fcheck_in) else: pdb = app.pdbxfile.PDBxFile(self.cif_cg) + else: + 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) From 3ad7f9fc1f8d1c4b5a065b4b7e472321ee399483 Mon Sep 17 00:00:00 2001 From: gitesei Date: Sat, 7 Mar 2026 13:02:36 +0100 Subject: [PATCH 5/7] fix memory cmap calc --- calvados/analysis.py | 4 ++-- examples/slab_IDR_MDP/prepare.py | 4 ++++ 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/calvados/analysis.py b/calvados/analysis.py index 1dd8a58..dc0575a 100644 --- a/calvados/analysis.py +++ b/calvados/analysis.py @@ -965,14 +965,14 @@ def calc_contact_map(path,sysname,output_path,chainid_dict={},is_slab=False,inpu if is_slab: mask_frames = np.where(chainid_dict[name_1] == chain_1)[0] else: - mask_frames = np.arange(traj.n_frames) + mask_frames = np.arange(traj.n_frames)#, True, dtype=bool) if len(mask_frames) > 0: for mf in mask_frames: d = md.compute_distances(traj[mf],pair_indices)[0] cm = (.5-.5*np.tanh((d-1.)/.3)).reshape(N_res_1,-1,N_res_2) cm = np.sum(cm,axis=1) cmap += cm - 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) diff --git a/examples/slab_IDR_MDP/prepare.py b/examples/slab_IDR_MDP/prepare.py index 658dd6e..bd9fac1 100644 --- a/examples/slab_IDR_MDP/prepare.py +++ b/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) From d075cbe86c0939e1bea9b5f67e9e21e67e380abd Mon Sep 17 00:00:00 2001 From: gitesei Date: Sat, 7 Mar 2026 13:04:19 +0100 Subject: [PATCH 6/7] Update README --- README.md | 1 - 1 file changed, 1 deletion(-) diff --git a/README.md b/README.md index 15b6170..77d5d6e 100644 --- a/README.md +++ b/README.md @@ -43,7 +43,6 @@ conda install -c conda-forge openmm=8.2.0 mdanalysis=2.9 mdtraj=1.11 ``` git clone https://github.com/KULL-Centre/CALVADOS.git cd CALVADOS -git checkout pdbx pip install . (or pip install -e .) ``` From dc9b38c63aa805b11a7bbb81289899373fe46e83 Mon Sep 17 00:00:00 2001 From: gitesei Date: Sat, 7 Mar 2026 13:25:45 +0100 Subject: [PATCH 7/7] Fix dependencies --- .github/workflows/python-app.yml | 4 ++-- calvados/analysis.py | 24 +++++++----------------- setup.py | 8 ++++++-- 3 files changed, 15 insertions(+), 21 deletions(-) diff --git a/.github/workflows/python-app.yml b/.github/workflows/python-app.yml index 0603e70..fec5f40 100644 --- a/.github/workflows/python-app.yml +++ b/.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 diff --git a/calvados/analysis.py b/calvados/analysis.py index dc0575a..19af9ab 100644 --- a/calvados/analysis.py +++ b/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]) @@ -922,11 +920,6 @@ def calc_contact_map(path,sysname,output_path,chainid_dict={},is_slab=False,inpu name_2 = name_1 chainid_dict[name_1] = np.arange(traj.top.n_chains) - print(name_1) - print(chainid_dict[name_1]) - print(name_2) - print(chainid_dict[name_2]) - N_res_1 = traj.top.chain(chainid_dict[name_1][0]).n_residues N_res_2 = traj.top.chain(chainid_dict[name_2][0]).n_residues @@ -955,23 +948,20 @@ def calc_contact_map(path,sysname,output_path,chainid_dict={},is_slab=False,inpu chainid_dict[name_2] = chainid_dict[name_1] cm_z = cmtraj.xyz[:,chainid_dict[name_1],2] # per-frame central-chain indices - ids_central = np.argmin(np.abs(cm_z),axis=1) - chainid_dict[name_1] = np.array([chainid_dict[name_1][idx] for idx in ids_central]) + chainid_dict[name_1] = np.argmin(np.abs(cm_z),axis=1) cmap = np.zeros((N_res_1,N_res_2)) for chain_1 in np.unique(chainid_dict[name_1]): surrounding_chains = traj.top.select(' or '.join([f'chainid {i:d}' for i in chainid_dict[name_2] if i != chain_1])) pair_indices = traj.top.select_pairs(f'chainid {chain_1:d}',surrounding_chains) if is_slab: - mask_frames = np.where(chainid_dict[name_1] == chain_1)[0] + mask_frames = chainid_dict[name_1] == chain_1 else: - mask_frames = np.arange(traj.n_frames)#, True, dtype=bool) - if len(mask_frames) > 0: - for mf in mask_frames: - d = md.compute_distances(traj[mf],pair_indices)[0] - cm = (.5-.5*np.tanh((d-1.)/.3)).reshape(N_res_1,-1,N_res_2) - cm = np.sum(cm,axis=1) - cmap += cm + mask_frames = np.full(traj.n_frames, True, dtype=bool) + if np.any(mask_frames): + d = md.compute_distances(traj[mask_frames],pair_indices) + 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 diff --git a/setup.py b/setup.py index 4979576..1cc033c 100644 --- a/setup.py +++ b/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,6 +23,9 @@ license='GNU GPL3', packages=find_packages(), install_requires=[ + "OpenMM>=8.4,<8.5", + "MDAnalysis>=2.10,<2.11", + "mdtraj>=1.11,<1.12", 'numpy', 'pandas', 'biopython', @@ -36,6 +39,7 @@ 'numba', 'scipy' ], + python_requires=">=3.13,<3.14", # include_package_data=True, package_data={'' : ['data/*.csv', 'data/*.yaml', 'data/templates/*']}, @@ -43,6 +47,6 @@ classifiers=[ 'Intended Audience :: Science/Research', 'Operating System :: POSIX :: Linux', - 'Programming Language :: Python :: 3>=3.7,<3.11', + 'Programming Language :: Python :: 3.13', ], )