#
#
#
                                                                                                                                
__author__ = "Randall J. Radmer"
__version__ = "1.0"
__doc__ = """Chemical Structure Cals, Including Coordinate Modification."""
                                                                                                                                
                                                                                                                                
import math

import simtk.utils.vecOps as vecOps

radPerDeg = math.pi/180
degPerRad = 180/math.pi


def resetBondLength(atom1, atom2, newLength):
    bondVec = vecOps.delta3(atom1, atom2)
    scaleLength = newLength/vecOps.mag3(bondVec)
    return vecOps.scale3(bondVec, scaleLength, atom2)

def addAtom(atom1, atom2, atom3,
            dihedralAngle=180.0, angle=109.5, bondLength=1.09):
    normal123 = vecOps.normalVec3(atom1, atom2, atom3)
    normal012 = vecOps.rotVec3(atom1,
                                       atom2,
                                       radPerDeg*dihedralAngle,
                                       normal123)
    bondVec12 = vecOps.unitVec3(vecOps.delta3(atom1, atom2))
    atom0 = vecOps.trans3(atom1, bondVec12)
    atom0 = resetBondLength(atom0, atom1, bondLength)
    atom0 = vecOps.trans3(atom0, atom1, -1)
    atom0 = vecOps.rotVec3(atom1,
                                   vecOps.trans3(atom1, normal012),
                                   -radPerDeg*(180-angle),
                                   atom0)
    atom0 = vecOps.trans3(atom0, atom1)
    return atom0
 
def calcCenterGeom(coordX, coordY, coordZ):
    xTot = 0.0
    yTot = 0.0
    zTot = 0.0
    n=len(coordX)
    for i in range(n):
        xTot += coordX[i]
        yTot += coordY[i]
        zTot += coordZ[i]
    return (xTot/n, yTot/n, zTot/n)

def getBondLength(atomA, atomB):
    return vecOps.mag3(vecOps.delta3(atomA, atomB))

def getBondAngle(atomA, atomB, atomC):
    dAB = vecOps.unitVec3(vecOps.delta3(atomA, atomB))
    dCB = vecOps.unitVec3(vecOps.delta3(atomC, atomB))
    cosAC = vecOps.dotProd(dAB, dCB)
    return degPerRad*math.acos(cosAC)

def getDihedralAngle(atomA, atomB, atomC, atomD):
    crossABC=vecOps.normalVec3(atomA, atomB, atomC)
    crossBCD=vecOps.normalVec3(atomB, atomC, atomD)
    cosDihedralAngle = min(max(vecOps.dotProd(crossABC, crossBCD),-1.0), 1.0)
    dihedralAngle = degPerRad*math.acos(cosDihedralAngle)

    dCB = vecOps.unitVec3(vecOps.delta3(atomC, atomB))
    sinDihedralAngleVec = vecOps.crossProd3(crossABC, crossBCD)
    dotSinDihedralAngle = vecOps.dotProd(sinDihedralAngleVec, dCB)
    try:
        signDihedralAngle =  dotSinDihedralAngle/abs(dotSinDihedralAngle)
    except ZeroDivisionError:
        signDihedralAngle =  1

    return signDihedralAngle*dihedralAngle

def getChirality(atom0, atomA, atomB, atomX):
    crossA0B=vecOps.normalVec3(atomA, atom0, atomB)
    dX0 = vecOps.delta3(atomX, atom0)
    return vecOps.dotProd(crossA0B, dX0)

def getHBondEnergy(vAtomN, vAtomH, vAtomO, vAtomC):
    rON = vecOps.mag3(vecOps.delta3(vAtomO, vAtomN))
    rCN = vecOps.mag3(vecOps.delta3(vAtomC, vAtomN))
    rOH = vecOps.mag3(vecOps.delta3(vAtomO, vAtomH))
    rCH = vecOps.mag3(vecOps.delta3(vAtomC, vAtomH))
    scale = 332*0.42*0.2
    return scale*(+1/rON +1/rCH -1/rOH -1/rCN)

if __name__=='__main__':
    atomRef = vecOps.unitZ
    atomC0 = [0.0, 0.0, 0.0]
    atomH0 = vecOps.trans3(atomC0, vecOps.unitX)
    print 'pos C0 = %s' % str(atomC0)
    print 'pos H0 = %s' % str(atomH0)
    print 'extended bond = %s' % str(resetBondLength(atomH0, atomC0, 1.5))

    coordX=[atomRef[0], atomC0[0], atomH0[0]]
    coordY=[atomRef[1], atomC0[1], atomH0[1]]
    coordZ=[atomRef[2], atomC0[2], atomH0[2]]
    (coordX, coordY, coordZ) = \
        changeHtoCH3(coordX, coordY, coordZ)

    atomC1 = [coordX[1], coordY[1], coordZ[1]]
    atomH1 = [coordX[2], coordY[2], coordZ[2]]
    atomH2 = [coordX[3], coordY[3], coordZ[3]]
    atomH3 = [coordX[4], coordY[4], coordZ[4]]

    print 'pos C0  = %s' % str(atomC0)
    print 'pos C1  = %s' % str(atomC1)
    print 'pos HC1 = %s' % str(atomH1)
    print 'pos HC2 = %s' % str(atomH2)
    print 'pos HC3 = %s' % str(atomH3)

    print 'delta3 C1, C0 = %s' \
        % vecOps.mag3(vecOps.delta3(atomC1, atomC0))
    print 'delta3 HC1, C1 = %s' \
        % vecOps.mag3(vecOps.delta3(atomH1, atomC1))
    print 'delta3 HC2, C1 = %s' \
        % vecOps.mag3(vecOps.delta3(atomH2, atomC1))
    print 'delta3 HC3, C1 = %s' \
        % vecOps.mag3(vecOps.delta3(atomH3, atomC1))
    print 'angle C0, C1, HC1 = %s' \
        % getBondAngle(atomC0, atomC1, atomH1)
    print 'angle C0, C1, HC2 = %s' \
        % getBondAngle(atomC0, atomC1, atomH2)
    print 'angle C0, C1, HC3 = %s' \
        % getBondAngle(atomC0, atomC1, atomH3)
    print 'angle HC2, C1, HC1 = %s' \
        % getBondAngle(atomH2, atomC1, atomH1)
    print 'angle HC3, C1, HC1 = %s' \
        % getBondAngle(atomH3, atomC1, atomH1)
    print 'angle HC3, C1, HC2 = %s' \
        % getBondAngle(atomH3, atomC1, atomH2)
    print 'atomH1, atomC1, atomC0, atomRef'
    print atomH1, atomC1, atomC0, atomRef
    print 'dihedralAngle HC1, C1, C0, Ref = %0.2f' \
        % getDihedralAngle(atomH1, atomC1, atomC0, atomRef)
    print 'dihedralAngle HC2, C1, C0, Ref = %0.2f' \
        % getDihedralAngle(atomH2, atomC1, atomC0, atomRef)
    print 'dihedralAngle HC3, C1, C0, Ref = %0.2f' \
        % getDihedralAngle(atomH3, atomC1, atomC0, atomRef)