q_ffld2q.py 文件源码

def calc_angle(ac1, ac2, ac3, fk, th0):
    th0 = math.pi/180.0 * th0   # degrees to radians
    rji1 = ac1[0] - ac2[0]
    rji2 = ac1[1] - ac2[1]
    rji3 = ac1[2] - ac2[2]
    rjk1 = ac3[0] - ac2[0]
    rjk2 = ac3[1] - ac2[1]
    rjk3 = ac3[2] - ac2[2]

    # get the angle from the dot product equation ( A*B = |A|*|B|*cos(theta) )
    # where A and B are vectors bewteen atoms (1->2 and 2->3)
    bji2inv = 1./(rji1**2 + rji2**2 + rji3**2)
    bjk2inv = 1./(rjk1**2 + rjk2**2 + rjk3**2)
    bjiinv = math.sqrt(bji2inv)
    bjkinv = math.sqrt(bjk2inv)
    scp = (rji1*rjk1 + rji2*rjk2 + rji3*rjk3)
    scp = scp * bjiinv* bjkinv
    if scp > 1.0:
        scp =  1.0
    elif scp < -1.0:
        scp = -1.0
    theta = math.acos(scp)   # in radians
    dtheta = (theta-th0)  
    en = 0.5*fk*dtheta**2 
    return en,theta*180.0/math.pi    # kcal/mol and degrees