From b772a35e33cee205a3dad1657bb181edd34abf28 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=B6ren=20von=20B=C3=BClow?= Date: Mon, 23 Feb 2026 17:46:06 +0000 Subject: [PATCH] Fix memory issues in slab analysis --- calvados/analysis.py | 16 +++++++++------- examples/slab_IDR_MDP/prepare.py | 4 ++++ 2 files changed, 13 insertions(+), 7 deletions(-) diff --git a/calvados/analysis.py b/calvados/analysis.py index b6e86c3..367062d 100644 --- a/calvados/analysis.py +++ b/calvados/analysis.py @@ -963,13 +963,15 @@ 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)#, 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) + np.save(output_path+f'/{sysname:s}_{name_1:s}_{name_2:s}_cmap.npy',cmap) \ No newline at end of file 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)