#!/usr/bin/env python
"""Contains useful functions for the determination of CV values and gradients."""
import numpy as np
from numpy.linalg import norm
import itertools as it
import matplotlib.pyplot as plt
from .. import constants as c
from .. import atoms as at
[docs]def bond(coord):
ij = coord[0]-coord[1]
return norm(ij)
[docs]def dbond(coord):
ij = coord[0]-coord[1]
d_di = ij/norm(ij)
d_dj = -d_di
return np.array([d_di, d_dj])
[docs]def angle(coord):
ij = coord[0]-coord[1]
kj = coord[2]-coord[1]
return np.arccos(np.dot(ij, kj)/(norm(ij)*norm(kj)))*180/np.pi
[docs]def dangle(coord):
a = angle(coord)*np.pi/180
ij = coord[0]-coord[1]
kj = coord[2]-coord[1]
ij0 = ij/norm(ij)
kj0 = kj/norm(kj)
d_di = 1/np.sqrt(1-np.cos(a)**2)*(1/norm(ij))*(ij0*np.cos(a)-kj0)*180/np.pi
d_dk = 1/np.sqrt(1-np.cos(a)**2)*(1/norm(kj))*(kj0*np.cos(a)-ij0)*180/np.pi
d_dj = -d_di-d_dk
return np.array([d_di, d_dj, d_dk])
[docs]def torsion(coord):
ji = (coord[1]-coord[0])/norm(coord[1]-coord[0])
jk = (coord[1]-coord[2])/norm(coord[1]-coord[2])
lk = (coord[3]-coord[2])/norm(coord[3]-coord[2])
cross1 = np.cross(ji, jk)
cross2 = np.cross(jk, lk)
cross3 = np.cross(cross1, jk)
dot1 = np.dot(cross1, cross2)
dot2 = np.dot(cross3, cross2)
return np.arctan2(dot2, dot1)*180/np.pi
[docs]def dtorsion(coord):
"""
Implementation according to van Schalk et al, J. Mol. Biol. 234, 751 /
https://salilab.org/modeller/9v6/manual/node436.html
"""
ij = coord[0]-coord[1]
kj = coord[2]-coord[1]
kl = coord[2]-coord[3]
mj = np.cross(ij, kj)
nk = np.cross(kj, kl)
dot1 = np.dot(ij, kj)
dot2 = np.dot(kl, kj)
d_di = (norm(kj)/norm(mj)**2)*mj*180/np.pi
d_dl = -(norm(kj)/norm(nk)**2)*nk*180/np.pi
d_dj = ((dot1/norm(kj)**2)-1)*d_di-(dot2/norm(kj)**2)*d_dl
d_dk = ((dot2/norm(kj)**2)-1)*d_dl-(dot1/norm(kj)**2)*d_di
return np.array([d_di, d_dj, d_dk, d_dl])
[docs]def cn_i(r, d, n, m):
return (1-(r/d)**n)/(1-(r/d)**m)
[docs]def dcn_i(r, d, n, m):
t1 = m*(1-(r/d)**n)*(r/d)**m/(r*(1-(r/d)**m)**2)
t2 = n*(r/d)**n/(r*(1-(r/d)**m))
return t1-t2
[docs]def get_d(symb, i, j):
tol = 1.0
return (at.rcov[symb[i]]+at.rcov[symb[j]])*tol/c.a0
[docs]def cn(symb, coord, ndx1, n, m, d, ndx2=[]):
cn = 0.0
ds_dr = np.zeros((len(coord), 3))
for i, coord_i in enumerate(coord):
# if no set of reference atoms is defined or if i is in the list of reference atoms:
if (i != ndx1 and not ndx2) or (i != ndx1 and i in ndx2):
r = bond([coord[ndx1], coord_i])
# determine d for the current couple of atom types
if ndx2: # if reference atoms are defined, d contains only values for atoms in ndx2
d_i = d[ndx2.index(i)]
else: # if no reference atoms are defined, d contains values for all atoms
d_i = d[i]
cn += cn_i(r, d_i, n, m)
dcn = dcn_i(r, d_i, n, m)
dr = dbond([coord[ndx1], coord_i])
ds_dr[ndx1] += dcn*dr[0]
ds_dr[i] = dcn*dr[1]
return cn, ds_dr
[docs]def check_cn_parameters(n=6, m=12, d=1.):
"""
Plot the coordination number against the distance between two atoms.
Parameters
----------
n : int
first exponent
m : int
second exponent
d : float
distance threshold in angstrom
"""
d = d/c.a0
plt.title(r"$n = %i, m = %i, d = %.2f \mathrm{\AA}$" % (n, m, d*c.a0))
x = np.linspace(0., 4., 100)
cn = cn_i(x/c.a0, d, n, m)
plt.xlabel(r"bond length / $\mathrm{\AA}$")
plt.ylabel(r"coordination number")
plt.axvline(x=d*c.a0, color="grey", ls="dashed")
plt.plot(x, cn)
plt.tight_layout(pad=0.1)
plt.savefig("cn_parameters.png", dpi=600)
plt.show()
[docs]def wiener(symbols, coords, exclude_h=True):
nat = len(coords)
wiener = 0.0
dwiener = np.zeros_like(coords)
for p in it.combinations(range(nat), 2):
if exclude_h == True and ("H" in [symbols[p_i] for p_i in p]):
continue
wiener += bond(coords[[p]])
dwiener[[p]] += dbond(coords[[p]])
return wiener, dwiener