(* ::Package:: *)

OBtrans::usage = "OBtrans takes as input C-O-C or C-O-H
connectivity triples from arbitrary molecular structures. It translates the
triangles to x=y=z=0 for the O atom; then rotaton around the x-axis turns 
the second atom into the xy plane, a rotation around the z-axis puts it onto 
the x-axis, a third rotation around the x-axis moves the third atom into
the xy plane. Now the lone pairs of the O atom are attached perpendicular 
to the xy plane and plotted. Then the four transformations are reversed, 
including the coordinates of the lone pairs. This puts the three atoms 
exactly back into their position in the structure with lone pairs added. 
The input requires:
1) list Ok composed of {{1,1,1},{1,-1,1}...}}, one sublist {1,1,1} for every
O atom. It is used to correct the sign of the x,y coordinates of the lone
pairs which sometimes is not chosen correctly.
2) list rlos of O,S,Se .. radii." 


k=1;
Plota= Table[{0}, {i, 1, hj}]; 
Plotc= Plota;
Plotd =Plota; 
Plote= Plota;
Plotb0 =Plota; 


Do[(c = coch[[k]];
e1=c[[1]]; e2=c[[2]]; c[[2]]=e1; c[[1]]=e2;  
a1 = dc[[c[[2]],c[[1]]]]; 
a2 = dc[[c[[2]],c[[3]]]]; 
d1 = dc[[c[[1]],c[[3]]]]; 
\[Alpha] = ArcCos[(a1^2 + a2^2 - d1^2)/(2*a1*a2)];
\[Delta]=0.025*rlos[[k]];
\[Beta] = ArcSin[(rlos[[k]]+\[Delta])/(u2[[2]] + rlos[[k]])]; 
L = (u2[[2]] + rlos[[k]])*Cos[\[Beta]];
\[Gamma] = \[Pi]/2 - \[Alpha]/2;  
v = Table[Join[t[[c[[i]]]], {1}], {i, 1, P}]; 
x2 = t[[c[[2]],1]]; 
y2 = t[[c[[2]],2]]; 
z2 = t[[c[[2]],3]]; 
Q = {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {-x2, -y2, -z2, 1}}; 
p = v . Q; 
p = Table[Take[p[[i]], 3], {i, 1, P}]; 
T1 = ArcTan[p[[1,3]]/p[[1,2]]]; 
U = p; 
T = RotationMatrix[T1, {1, 0, 0}]; 
p = U . T; 
T2 = ArcTan[p[[1,2]]/p[[1,1]]]; 
U = p; 
T = RotationMatrix[T2, {0, 0, 1}]; 
p = U . T; 
T3 = ArcTan[p[[3,3]]/p[[3,2]]]; 
U = p; 
T = RotationMatrix[T3, {1, 0, 0}]; 
p = U . T; 
nk = p; nl = p; 
nk[[2,1]] = L*Sin[\[Gamma]];  nk[[2,2]] = L*Cos[\[Gamma]];  nk[[2,3]] =  rlos[[k]]+\[Delta]; 
nl[[2,1]] = nk[[2,1]]; nl[[2,2]] = nk[[2,2]]; nl[[2,3]] = -rlos[[k]]-\[Delta];
nk[[2]]=Ok[[k]]*nk[[2]]; nl[[2]]=Ok[[k]]*nl[[2]];
Plota[[k]] = Table[Graphics3D[{Opacity[0.5], Sphere[p[[i]], u2[[c[[i]]]]]}], {i, 1, P}]; 
Plotc[[k]] = Table[Graphics3D[{Tube[{p[[i]], p[[j]]}, 0.02]}], {i, 1, P - 1}, {j, i + 1, P}]; 
Plotd[[k]]= Graphics3D[{Opacity[0.5], Lighter[Red], Sphere[nk[[2]], rlos[[k]]]}]; 
Plote[[k]]= Graphics3D[{Opacity[0.5], Lighter[Red], Sphere[nl[[2]], rlos[[k]]]}]; 
Plotb0[[k]] = Table[Graphics3D[Text[c[[i]], p[[i]], {1, 0}]], {i, 1, P}];
U = p; 
U1 = nk; 
U2 = nl; 
T = RotationMatrix[-T3, {1, 0, 0}]; 
p = U . T; 
nk = U1 . T; 
nl = U2 . T; 
U = p; 
U1 = nk; 
U2 = nl; 
T = RotationMatrix[-T2, {0, 0, 1}]; 
p = U . T; 
nk = U1 . T; 
nl = U2 . T; 
U = p; 
U1 = nk; 
U2 = nl; 
T = RotationMatrix[-T1, {1, 0, 0}]; 
p = U . T;
nk = U1 . T; 
nl = U2 . T; 
Clear[v, v1, v2]; 
v = Table[Join[p[[i]], {1}], {i, 1, P}]; 
v1 = Table[Join[nk[[i]], {1}], {i, 1, P}]; 
v2 = Table[Join[nl[[i]], {1}], {i, 1, P}]; 
Q = {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {x2, y2, z2, 1}}; 
p = v . Q; 
nk = v1 . Q; 
nl = v2 . Q; 
p = Table[Take[p[[i]], 3], {i, 1, P}]; 
nk = Table[Take[nk[[i]], 3], {i, 1, P}]; 
nl = Table[Take[nl[[i]], 3], {i, 1, P}]; 
tpso[[k]] = nk[[2]]; 
tpsu[[k]] = nl[[2]]), {k, hj}]; 


Table[Show[Plotd[[i]],Plote[[i]], Plota[[i]], Plotb0[[i]], Plotc[[i]], 
{Axes -> Automatic, AxesLabel -> {xl, yl, zl}, 
AspectRatio -> Automatic, Boxed -> False, 
PlotRange -> Automatic, SphericalRegion -> True}],{i,1,hj}]