Skip to content
Snippets Groups Projects
Commit 1ec1c794 authored by Javier González-Delgado's avatar Javier González-Delgado
Browse files

Update get_contacts.ipynb

parent 53568a6b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
concha/get_contacts.ipynb
+ 2
− 2
View file @ 1ec1c794
%% Cell type:code id:5417c4b3 tags:
``` python
import numpy as np
import os
import mdtraj as md
import itertools
import pandas as pd
import warnings #Optional
warnings.filterwarnings("ignore") #Optional
```
%% Cell type:code id:57a23fbe tags:
``` python
def get_contacts(coor_conf, threshold_file, assort = False):
L = int(0.5*(1 + np.sqrt(1 + 8*np.shape(coor_conf)[0])))
pos_pairs = np.array(list(itertools.combinations(range(L), 2)))
contact_thresholds = pd.read_csv(threshold_file, sep=" ", header=0)
contact_thresholds['th11'] = np.deg2rad(contact_thresholds['th11'])
contact_thresholds['th12'] = np.deg2rad(contact_thresholds['th12'])
contact_thresholds['th13'] = np.deg2rad(contact_thresholds['th13'])
contact_thresholds['th21'] = np.deg2rad(contact_thresholds['th21'])
contact_thresholds['th22'] = np.deg2rad(contact_thresholds['th22'])
contact_thresholds['th23'] = np.deg2rad(contact_thresholds['th23'])
add = pd.DataFrame(contact_thresholds)
add.columns = contact_thresholds.columns
add = add.loc[(add.AA1 - add.AA2 != 0)]
add[['AA1','AA2']] = add[['AA2','AA1']].values
contact_thresholds = pd.concat([contact_thresholds,add], ignore_index = True)
contact_thresholds['AA1'] = contact_thresholds.AA1.astype('int')
contact_thresholds['AA2'] = contact_thresholds.AA2.astype('int')
contact_thresholds['range'] = contact_thresholds.range.astype('int')
contact_thresholds['AA1-AA2-range'] = contact_thresholds.AA1.astype(str) + '-' + contact_thresholds.AA2.astype(str) + '-' + contact_thresholds.range.astype(str)
contact_thresholds = contact_thresholds[['AA1-AA2-range','delta_min','delta_max','delta', 'th11', 'th12', 'th13', 'th21', 'th22', 'th23', 'delta_se3_min', 'delta_se3_max']]
mins = np.minimum(contact_thresholds.delta_min.values,contact_thresholds.delta_max.values)
maxs = np.maximum(contact_thresholds.delta_min.values,contact_thresholds.delta_max.values)
contact_thresholds.delta_min = mins
contact_thresholds.delta_max = maxs
coor_conf = pd.DataFrame(np.concatenate([coor_conf, pos_pairs], axis = 1),
columns = ['coor_x','coor_y','coor_z','or1_x','or1_y','or1_z','or2_x','or2_y','or2_z','AA1','AA2','pos1','pos2'])
coor_conf.range = np.abs(coor_conf['pos1'] - coor_conf['pos2'])
coor_conf.range = (coor_conf.range*(coor_conf.range<5) + 5*(coor_conf.range>=5)).astype('int')
coor_conf['AA1'] = coor_conf.AA1.astype('int')
coor_conf['AA2'] = coor_conf.AA2.astype('int')
coor_conf['AA1-AA2-range'] = coor_conf.AA1.astype(str) + '-' + coor_conf.AA2.astype(str) + '-' + coor_conf.range.astype(str)
#coor_conf = coor_conf.join(vaex.from_pandas(contact_thresholds), left_on = 'AA1-AA2-range', right_on = 'AA1-AA2-range', how = 'left')
coor_conf = coor_conf.merge(contact_thresholds, on = 'AA1-AA2-range', how = 'left')
# For range <= 4, we correct distance by admissible orientations
coor_conf['min_th1'] = np.nanmin([np.abs(np.arccos(coor_conf['or1_x']) - coor_conf['th11']),
np.abs(np.arccos(coor_conf['or1_x']) - coor_conf['th12']),
np.abs(np.arccos(coor_conf['or1_x']) - coor_conf['th13'])], axis = 0)
coor_conf['min_th2'] = np.nanmin([np.abs(np.arccos(coor_conf['or2_z']) - coor_conf['th21']),
np.abs(np.arccos(coor_conf['or2_z']) - coor_conf['th22']),
np.abs(np.arccos(coor_conf['or2_z']) - coor_conf['th23'])], axis = 0)
coor_conf['dis_th1'] = 0.5*(np.sin(coor_conf.min_th1)**2*(coor_conf.min_th1 < np.pi/2) + (1 + np.cos(coor_conf.min_th1)**2)*(coor_conf.min_th1 >= np.pi/2))
coor_conf['dis_th2'] = 0.5*(np.sin(coor_conf.min_th2)**2*(coor_conf.min_th2 < np.pi/2) + (1 + np.cos(coor_conf.min_th2)**2)*(coor_conf.min_th2 >= np.pi/2))
coor_conf['dis_th1'] = 0.5*(np.sin(coor_conf.min_th1)**2*(coor_conf.min_th1 < np.pi/2) + (1 - np.cos(coor_conf.min_th1)**2)*(coor_conf.min_th1 >= np.pi/2))
coor_conf['dis_th2'] = 0.5*(np.sin(coor_conf.min_th2)**2*(coor_conf.min_th2 < np.pi/2) + (1 - np.cos(coor_conf.min_th2)**2)*(coor_conf.min_th2 >= np.pi/2))
alpha = beta = 0.5
coor_conf['dis_r3'] = np.sqrt(coor_conf.coor_x**2 + coor_conf.coor_y**2 + coor_conf.coor_z**2)
coor_conf['dis_or'] = np.sqrt(alpha*coor_conf['dis_th1']**2 + beta*coor_conf['dis_th1']**2)
argtanh = lambda x: 0.5*np.log((1+x)/(1-x))
coor_conf[coor_conf.delta_min < 2].delta_min = 2
coor_conf[coor_conf.delta_max <= 2].delta_max = 3
coor_conf['d'] = np.log(argtanh(1/coor_conf.delta_min))/np.log(coor_conf.delta_min/coor_conf.delta_max)
coor_conf['w_or_pos'] = 1-np.tanh((coor_conf.dis_r3/coor_conf.delta_max)**coor_conf.d)
coor_conf['a'] = 0.5*np.sqrt(argtanh(1-coor_conf.w_or_pos))
coor_conf['w_or_or'] = 1-np.tanh((2*(coor_conf.dis_or+coor_conf.a))**2)
coor_conf[coor_conf.w_or_pos == 0].w_or_or = 0
coor_conf['dis_or'] = coor_conf['dis_or'].fillna(0)
coor_conf['w_or_or'] = coor_conf['w_or_or'].fillna(0)
coor_conf['dis_se3'] = np.sqrt((1-coor_conf.w_or_or)**2*coor_conf.dis_r3**2 + coor_conf.w_or_or**2*coor_conf.dis_or**2)
coor_conf.delta_se3_min[coor_conf['delta_se3_min'].isnull()] = coor_conf[coor_conf['delta_se3_min'].isnull()]['delta_min']
coor_conf.delta_se3_min[coor_conf['delta_se3_min'] < 2] = 2
coor_conf.delta_se3_max[coor_conf['delta_se3_max'].isnull()] = coor_conf[coor_conf['delta_se3_max'].isnull()]['delta_max']
coor_conf.delta_se3_max[coor_conf['delta_se3_max'] <= 2] = 3
coor_conf['d_se3'] = np.log(argtanh(1/coor_conf.delta_se3_min))/np.log(coor_conf.delta_se3_min/coor_conf.delta_se3_max)
coor_conf['w_dis_se3'] = 1-np.tanh((coor_conf.dis_se3/coor_conf.delta_se3_max)**coor_conf.d_se3)
coor_conf = coor_conf[['pos1','pos2','w_dis_se3','AA1', 'AA2']]
coor_conf['pos1'] = coor_conf.pos1.astype('int') + 1
coor_conf['pos2'] = coor_conf.pos2.astype('int') + 1
if assort:
return coor_conf
else:
return coor_conf.w_dis_se3
```
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment