求助关系python语言解决向量问题

少年爱吃地三鲜

求助各位老师，我是用MDAnalysis计算二面角的， 图二是 根据的公式， 图一是找到的算法代码。 请教几个小问题，1：请问b0 为啥要乘以负一  2： v和w是什么含义？ 麻烦各位老师了

Daniel_Arndt

我以前写过一段算二面角的awk代码，当时是照着stackoverflow上的一个回答写出来的（记不清具体网址了），用过多次。我用的方法应该是你的图二中由三个向量b1、b2和b3定义的方法。

function fdihedral(d1,d2,d3,d4,    x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,norm1,x5,y5,z5,norm2,x6,y6,z6)
{
    x1 = Atom[d1,1];
    y1 = Atom[d1,2];
    z1 = Atom[d1,3];
    x2 = Atom[d2,1];
    y2 = Atom[d2,2];
    z2 = Atom[d2,3];
    x3 = Atom[d3,1];
    y3 = Atom[d3,2];
    z3 = Atom[d3,3];
    x4 = Atom[d4,1];
    y4 = Atom[d4,2];
    z4 = Atom[d4,3];
    norm1 = sqrt((y1*(z2-z3)+y2*(z3-z1)+y3*(z1-z2))*(y1*(z2-z3)+y2*(z3-z1)+y3*(z1-z2))+(z1*(x2-x3)+z2*(x3-x1)+z3*(x1-x2))*(z1*(x2-x3)+z2*(x3-x1)+z3*(x1-x2))+(x1*(y2-y3)+x2*(y3-y1)+x3*(y1-y2))*(x1*(y2-y3)+x2*(y3-y1)+x3*(y1-y2)));
    x5 = (y1*(z2-z3)+y2*(z3-z1)+y3*(z1-z2))/norm1;
    y5 = (z1*(x2-x3)+z2*(x3-x1)+z3*(x1-x2))/norm1;
    z5 = (x1*(y2-y3)+x2*(y3-y1)+x3*(y1-y2))/norm1;
    norm2 = sqrt((x3-x2)*(x3-x2)+(y3-y2)*(y3-y2)+(z3-z2)*(z3-z2));
    x6 = (y5*(z3-z2)-(y3-y2)*z5)/norm2;
    y6 = (z5*(x3-x2)-(z3-z2)*x5)/norm2;
    z6 = (x5*(y3-y2)-(x3-x2)*y5)/norm2;
    return -57.29577951*atan2(x6*(y2*(z3-z4)+y3*(z4-z2)+y4*(z2-z3))+y6*(z2*(x3-x4)+z3*(x4-x2)+z4*(x2-x3))+z6*(x2*(y3-y4)+x3*(y4-y2)+x4*(y2-y3)),x5*(y2*(z3-z4)+y3*(z4-z2)+y4*(z2-z3))+y5*(z2*(x3-x4)+z3*(x4-x2)+z4*(x2-x3))+z5*(x2*(y3-y4)+x3*(y4-y2)+x4*(y2-y3)));
}

里面的x1、y1、z1、x2、y2、z2、x3、y3、z3、x4、y4、z4就是原子的直角坐标。最后的57.29577951是为了把弧度制转成角度制。

少年爱吃地三鲜

Daniel_Arndt 发表于 2021-1-7 07:29
我以前写过一段算二面角的awk代码，当时是照着stackoverflow上的一个回答写出来的（记不清具体网址了），用 ...

感谢您的回复！ 请问您能看出那个v 和w是什么意思吗

highlight

这里 https://stackoverflow.com/questi ... ordinates-in-python
乘负一的原因是

"Flip" the first vector so that eclipsing vectors have dihedral=0

至于 v 和 w，他用的不是法向量夹角的求法，而是两平面内垂直于交线的向量夹角

def new_dihedral(p):
    """Praxeolitic formula
    1 sqrt, 1 cross product"""
    p0 = p[0]
    p1 = p[1]
    p2 = p[2]
    p3 = p[3]

    b0 = -1.0*(p1 - p0)
    b1 = p2 - p1
    b2 = p3 - p2

    # normalize b1 so that it does not influence magnitude of vector
    # rejections that come next
    b1 /= np.linalg.norm(b1)

    # vector rejections
    # v = projection of b0 onto plane perpendicular to b1
    #   = b0 minus component that aligns with b1
    # w = projection of b2 onto plane perpendicular to b1
    #   = b2 minus component that aligns with b1
    v = b0 - np.dot(b0, b1)*b1
    w = b2 - np.dot(b2, b1)*b1

    # angle between v and w in a plane is the torsion angle
    # v and w may not be normalized but that's fine since tan is y/x
    x = np.dot(v, w)
    y = np.dot(np.cross(b1, v), w)
    return np.degrees(np.arctan2(y, x))

PS：有注释的代码才是真生产力啊~

liyuanhe211

代码截图的分辨率起码得能看得清是等号还是减号吧。。。

令二面角是A1-A2-A3-A4组成的，则算二面角的俩向量应该是向量A2->A1和向量A3->A4之间的二面角，而不是向量A1->A2和向量A3->A4的二面角，所以它乘个负一。。。为啥要减个错的再乘负一我也不知道

设有一单位向量a，(b·a)a 的几何意义是求b在平行于a方向的分向量。b - (b·a)a 就是求b在垂直于a方向的分向量.

少年爱吃地三鲜

liyuanhe211 发表于 2021-1-7 09:45
代码截图的分辨率起码得能看得清是等号还是减号吧。。。

令二面角是A1-A2-A3-A4组成的，则算二面角的俩 ...

感谢您的教导!!!!

liuyuje714

按照图二的公式写就好了。对于做MD轨迹分析的话你还是应该看看一些源码的，比如mdanalysis的源中就有相应的二面角计算函数，你主要还应该考虑跨周期性分子的问题。物理知识只是告诉你如何做和原理是什么，真正用到实际问题考虑要更加全面。

少年爱吃地三鲜

def dihedral(p1,p2,p3,p4):
    '''p1,p2,p3,p4 are atoms'''
    c1=p1.position
    c2=p2.position
    c3=p3.position
    c4=p4.position
    b0 = -1.0*(c2-c1)
    b1 = c3 - c2
    b2 = c4 - c3
    b1 /= np.linalg.norm(b1)
    v = b0 - np.dot(b0, b1)*b1
    w = b2 - np.dot(b2, b1)*b1
    x = np.dot(v, w)
    y = np.dot(np.cross(b1, v), w)   
    return np.degrees(np.arctan2(y, x))

少年爱吃地三鲜

highlight 发表于 2021-1-7 09:41
这里 https://stackoverflow.com/questi ... ordinates-in-python
乘负一的原因是

非常感谢您!

我本是个娃娃

2019年7月写的，当时是为了解决另一个帖子的关于乙烷的二面角问题。

diangle.py (1.28 KB, 下载次数 Times of downloads: 14)

2021-1-7 12:09 上传 Uploaded

点击下载
Click to download

liuyuje714

我本是个娃娃 发表于 2021-1-7 12:10
2019年7月写的，当时是为了解决另一个帖子的关于乙烷的二面角问题。

你这个实际用起来应该不一定对，你没有考虑正负，正常情况应该是[-PI, PI].

少年爱吃地三鲜

我本是个娃娃 发表于 2021-1-7 12:10
2019年7月写的，当时是为了解决另一个帖子的关于乙烷的二面角问题。

liyuanhe211

我本是个娃娃 发表于 2021-1-7 12:10
2019年7月写的，当时是为了解决另一个帖子的关于乙烷的二面角问题。

学一下numpy吧，很好用的。

我本是个娃娃

liyuanhe211 发表于 2021-1-7 15:22
学一下numpy吧，很好用的。

学过。但是脑子一抽，就想不起来用……

少年爱吃地三鲜

liyuanhe211 发表于 2021-1-7 15:22
学一下numpy吧，很好用的。

老师，我想问一下您，我们自己构建函数时候用到的原子坐标可能随着轨迹突然跳到另一端，这个时候我们编写程序应该如何定义呢？ 还是说只能在gomacs 中 使用nojump？ 这个周期性边界自己写该如何考虑，麻烦您了