IPKM/mdevaluate
6
1
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Added docu to hbonds

Browse Source
This commit is contained in:
sebastiankloth 2022-11-04 08:46:31 +01:00
parent a38e601fcb
commit 88e4fb964b
1 changed files with 26 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

33
mdevaluate/distribution.py
View File

@ -360,10 +360,29 @@ def next_neighbor_distribution(atoms, reference=None, number_of_neighbors=4, bin
def hbonds(D, H, A, box, DA_lim=0.35, HA_lim=0.35, min_cos=np.cos(30*np.pi/180), full_output=False):
"""
Compute h-bond pairs
Args:
D: Set of coordinates for donators.
H: Set of coordinates for hydrogen atoms. Should have the same
length as D.
A: Set of coordinates for acceptors.
DA_lim (opt.): Minimum distance beteen donator and acceptor.
HA_lim (opt.): Minimum distance beteen hydrogen and acceptor.
min_cos (opt.): Minimum cosine for the HDA angle. Default is
equivalent to a maximum angle of 30 degree.
full_output (opt.): Returns additionally the cosine of the
angles and the DA distances
Return:
List of (D,A)-pairs in hbonds.
"""
def dist_DltA(D, H, A, box, max_dist=0.35):
ppoints, pind = pbc_points(D, box, thickness=max_dist+0.1, index=True)
Dtree = cKDTree(ppoints)
Atree = cKDTree(A)
Dtree = spatial.cKDTree(ppoints)
Atree = spatial.cKDTree(A)
pairs = Dtree.sparse_distance_matrix(Atree, max_dist, output_type='ndarray')
pairs = np.asarray(pairs.tolist())
pairs = np.int_(pairs[pairs[:,2] > 0][:,:2])
@ -372,8 +391,8 @@ def hbonds(D, H, A, box, DA_lim=0.35, HA_lim=0.35, min_cos=np.cos(30*np.pi/180),
def dist_AltD(D, H, A, box, max_dist=0.35):
ppoints, pind = pbc_points(A, box, thickness=max_dist+0.1, index=True)
Atree = cKDTree(ppoints)
Dtree = cKDTree(D)
Atree = spatial.cKDTree(ppoints)
Dtree = spatial.cKDTree(D)
pairs = Atree.sparse_distance_matrix(Dtree, max_dist, output_type='ndarray')
pairs = np.asarray(pairs.tolist())
pairs = np.int_(pairs[pairs[:,2] > 0][:,:2])
@ -386,9 +405,9 @@ def hbonds(D, H, A, box, DA_lim=0.35, HA_lim=0.35, min_cos=np.cos(30*np.pi/180),
else:
pairs = dist_AltD(D,H,A,box,DA_lim)
vDH = md.pbc.pbc_diff(D[pairs[:,0]], H[pairs[:,0]], box)
vDA = md.pbc.pbc_diff(D[pairs[:,0]], A[pairs[:,1]], box)
vHA = md.pbc.pbc_diff(H[pairs[:,0]], A[pairs[:,1]], box)
vDH = pbc_diff(D[pairs[:,0]], H[pairs[:,0]], box)
vDA = pbc_diff(D[pairs[:,0]], A[pairs[:,1]], box)
vHA = pbc_diff(H[pairs[:,0]], A[pairs[:,1]], box)
angles_cos = np.clip(np.einsum('ij,ij->i', vDH, vDA)/
np.linalg.norm(vDH,axis=1)/
np.linalg.norm(vDA,axis=1), -1, 1)