# -*- coding: utf-8 -*-
#!/usr/bin/env python2
"""
Functions for the calculation of Cremer-Pople parameters and their gradients.
@author: Karina Sultangaleeva
"""
import numpy as np
N=6
[docs]def schwerpunkt(coord_old):
sp=0
for i in coord_old:
sp+=i/N
coord_new=coord_old-sp
return coord_new
[docs]def drehung(coord_old):
n=find_n(coord_old)
coord_new=[np.dot(i,n) for i in coord_old]
return coord_new
[docs]def find_n(coord_old):
Rs1=0
Rs2=0
coord_new=np.copy(coord_old)
for i in range(N):
Rs1+=coord_new[i]*np.sin(i*2*np.pi/N)
Rs2+=coord_new[i]*np.cos(i*2*np.pi/N)
l=np.cross(Rs1,Rs2)
k=np.sqrt(sum(p*p for p in l))
n=l/k
return n
[docs]def find_phi_cos(coord_old):
x=sum([b*np.cos(a*np.pi*4./N) for a,b in enumerate(drehung(coord_old))])/find_q2(coord_old)
phi=np.arccos(x)
return phi
[docs]def find_phi_sin(coord_old):
x=sum([-b*np.sin(a*np.pi*4./N) for a,b in enumerate(drehung(coord_old))])/find_q2(coord_old)
phi=np.arcsin(x)
return phi
[docs]def find_phi(coord_old):
x=sum([b*np.cos(a*np.pi*4./N) for a,b in enumerate(drehung(coord_old))])/find_q2(coord_old)
y=-sum([b*np.sin(4*a*np.pi/N) for a,b in enumerate(drehung(coord_old))])/find_q2(coord_old)
phi=np.arctan2(y,x)
if phi<0:
phi+=2*np.pi
return phi*180./np.pi
[docs]def find_Q(coord_old):
coord_old=schwerpunkt(coord_old)
quadr_list=[i*i for i in drehung(coord_old)]
Q=(np.sum(quadr_list))**0.5
return Q
[docs]def find_q3(coord_old):
q3=np.sum([(N**(-0.5))*b*(-1)**a for a,b in enumerate(drehung(coord_old))])
return q3
[docs]def find_q2(coord_old):
q2=np.sqrt(find_Q(coord_old)**2-find_q3(coord_old)**2)
return q2
[docs]def find_theta(coord_old):
coord_old=schwerpunkt(coord_old)
y=find_q2(coord_old)/find_Q(coord_old)
x=find_q3(coord_old)/find_Q(coord_old)
atan=np.arctan2(y,x)
return atan*180./np.pi
[docs]def find_A(coord_old):
Rs1=0
Rs2=0
coord_new=np.copy(coord_old)
for i in range(N):
Rs1=sum([b*np.sin(2*a*np.pi/N) for a,b in enumerate(coord_new)])
Rs2=sum([b*np.cos(2*a*np.pi/N) for a,b in enumerate(coord_new)])
A=np.cross(Rs1,Rs2)
return A
[docs]def find_B(coord_old):
kp=find_A(coord_old)
betrag=np.sqrt(np.sum(i*i for i in kp))
return betrag
[docs]def reshape(func):
end=np.transpose(func,axes=[1,0,2])
return end
##################numeric derivative##################
[docs]def numeric_derive(func,coord_old,delta):
old=np.array(func(coord_old))
shape=np.array(np.zeros((old.size,N,3)))
for i in range(N):
for j in range(3):
coord_with_d=np.copy(coord_old)
coord_with_d[i,j]=coord_with_d[i,j]+delta
shape[:,i,j]=(func(coord_with_d)-old)/delta
return shape
[docs]def n_d(func,coord_old,delta):
old=np.array(func(coord_old))
shape=np.array(np.zeros((N)))
for i in range(N):
coord_with_d=np.copy(coord_old)
coord_with_d[i]=coord_with_d[i]+delta
shape[i]=(func(coord_with_d)-old)/delta
return shape
##################analytical derivative####################
[docs]def dn(coord):
shape=np.zeros((3,N,3))
A=find_A(coord)
B=find_B(coord)
dA=-d_A(coord)
dB=d_B(coord)
AdB=np.zeros((3,N,3))
for i in range(3):
AdB[i]=A[i]*dB
shape=(dA*B-AdB)/(B**2)
return shape
[docs]def d_B(coord):
dB=np.zeros((N,3))
Rs1=0
Rs2=0
coord_old=np.copy(coord)
for e in range(N):
Rs1+=coord_old[e]*np.sin(e*2*np.pi/N)
Rs2+=coord_old[e]*np.cos(e*2*np.pi/N)
A=np.cross(Rs1,Rs2)
B=np.sqrt(np.sum(p*p for p in A))
for t in range(N):
beta=t*np.pi*2./N
sin=np.sin(beta)
cos=np.cos(beta)
dRs1=np.zeros((3,3))
dRs2=np.zeros((3,3))
for w in range(3):
dRs1[w,w]=sin
dRs2[w,w]=cos
dB[t,0]=1./B*(A[1]*(Rs1[2]*cos-Rs2[2]*sin)+A[2]*(Rs2[1]*sin-Rs1[1]*cos))
dB[t,1]=1./B*(A[0]*(Rs2[2]*sin-Rs1[2]*cos)+A[2]*(Rs1[0]*cos-Rs2[0]*sin))
dB[t,2]=1./B*(A[0]*(Rs1[1]*cos-Rs2[1]*sin)+A[1]*(Rs2[0]*sin-Rs1[0]*cos))
return dB
[docs]def d_A(coord):
dA=np.zeros((3,N,3))
Rs1=0
Rs2=0
coord_old=np.copy(coord)
for e in range(N):
Rs1+=coord_old[e]*np.sin(e*2*np.pi/N)
Rs2+=coord_old[e]*np.cos(e*2*np.pi/N)
for t in range(N):
beta=t*np.pi*2./N
sin=np.sin(beta)
cos=np.cos(beta)
dRs1=np.zeros((3,3))
dRs2=np.zeros((3,3))
for w in range(3):
dRs1[w,w]=sin
dRs2[w,w]=cos
dA[:,t]=np.cross(dRs1,Rs2)+np.cross(Rs1,dRs2)
return dA
[docs]def Matrix():
shape=np.zeros((N,3,N,3))
for i in range(N):
for j in range(N):
if i==j:
shape[i,:,j]=np.diag(np.ones((3))*1-1./N)
else:
shape[i,:,j]=np.diag(np.ones((3))*(-1./N))
return shape
[docs]def dZ_dR(coord_old):
dR=np.zeros((N,3,N,3))
for i in range(N):
for k in range(N):
for l in range(3):
for m in range(3):
if i==k and l==m:
dR[i,l,k,m]=1
dern=dn(coord_old)
n=find_n(coord_old)
r=np.copy(coord_old)
result=np.tensordot(dR,n,axes=(1,0))+np.tensordot(r,dern,axes=(1,0))
return result
[docs]def dQ_dZ(coord_old):
shape=np.zeros((N))
for v in range(N):
shape[v]=drehung(coord_old)[v]/find_Q(coord_old)
return shape
[docs]def dQ(coord_old):
coord_old=schwerpunkt(coord_old)
shape=np.tensordot(dQ_dZ(coord_old),dZ_dR(coord_old),axes=(0,0))
return shape
[docs]def dPhi_dZ(coord_old):
up=sum([b*np.cos(a*np.pi*4./N) for a,b in enumerate(drehung(coord_old))])
down=-sum([b*np.sin(4*a*np.pi/N) for a,b in enumerate(drehung(coord_old))])
shape=np.zeros((N))
for u in range(N):
beta=np.pi*u*4./N
shape[-u]=(np.sin(beta)*up-down*np.cos(beta))/(up**2)/(1+(up/down)**(-2))
return shape
[docs]def dPhi(coord_old):
coord_old=schwerpunkt(coord_old)
shape=dPhi_dZ(coord_old)
end=np.tensordot(shape,dZ_dR(coord_old),axes=(0,0))
return end*180./np.pi
[docs]def dq3(coord_old):
shape=np.zeros((6))
for t in range(N):
shape[t]=N**(-0.5)*(-1)**t
end=np.tensordot(shape,dZ_dR(coord_old),axes=(0,0))
return end
[docs]def dq2(coord_old):
dq2=((find_Q(coord_old))**2-(find_q3(coord_old))**2)**(-0.5)*(find_Q(coord_old)*dQ(coord_old)-find_q3(coord_old)*dq3(coord_old))
return dq2
[docs]def dTheta(coord_old):
coord_old=schwerpunkt(coord_old)
A=find_q3(coord_old)
dA=np.zeros((6))
for t in range(N):
dA[t]=N**(-0.5)*(-1)**t
B=find_Q(coord_old)
dB=dQ_dZ(coord_old)
dtheta_z=(-(1-(A/B)**2)**(-0.5))*(dA*B-dB*A)/B/B
end=np.tensordot(dtheta_z,dZ_dR(coord_old), axes=(0,0))
return end*180./np.pi