(* ::Package:: *)

OTtrans::usage = "OTtrans attaches three lone pairs to O- in -X-O- of arbitrary molecules, given tj sets of three
atoms with the target as second. The triangle defines the plane in which one of the LP's sit. The
two others are equally far from O but -60 and +60 deg below or above the plane in the z-direction. Since
the routine sometimes does not choose the correct sign of the x or y coordinates, lists Okk, Okl, and Okm
must be given as input with {{1,1,1},{-1,1,1}...} for every lobe accepting or changing a sign."


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


Do[(c = conn[[k]]; 
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)]; 
L = u2[[c[[2]]]] + rlo; 
(*\[Beta] = ArcSin[rlo/L]; *)
\[Gamma] = 3*\[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 = Chop[U.T]; 
nk = nl = nm = p; 
nk[[2,1]] = p[[1]]+L/3; nk[[2,2]] = p[[2]]+0.942809*L; nk[[2,3]] = 0; 
nl[[2,1]] = nk[[2,1]]; nl[[2,2]] = p[[2]]-0.471405*L; nl[[2,3]] = 0.816496*L;
nm[[2,1]] = nk[[2,1]]; nm[[2,2]] = nl[[2,2]]; nm[[2,3]] = -nl[[2,3]];
nk[[2]]=Okk[[k]]*nk[[2]]; nl[[2]]=Okl[[k]]*nl[[2]]; nm[[2]]=Okm[[k]]*nm[[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]], rlo],Sphere[nl[[2]], rlo],Sphere[nm[[2]], rlo]}]; 
(*Plote[[k]]  = Graphics3D[{Opacity[0.5], Lighter[Red], Sphere[nl[[2]], rlo]}]; 
Plotf[[k]]  = Graphics3D[{Opacity[0.5], Lighter[Red], Sphere[nm[[2]], rlo]}]; *)
Plotb0[[k]]  = Table[Graphics3D[Text[c[[i]], p[[i]], {1, 0}]], {i, 1, P}];
U1 = nk; 
U2 = nl; 
U3 = nm;
T = RotationMatrix[-T3, {1, 0, 0}];  
nk = U1.T; 
nl = U2.T; 
nm = U3.T; 
U1 = nk; 
U2 = nl;
U3 = nm; 
T = RotationMatrix[-T2, {0, 0, 1}]; 
nk = U1.T; 
nl = U2.T; 
nm = U3.T; 
U1 = nk; 
U2 = nl;
U3 = nm; 
T = RotationMatrix[-T1, {1, 0, 0}]; 
nk = U1.T; 
nl = U2.T; 
nm = U3.T; 
v1 = Table[Join[nk[[i]], {1}], {i, 1, P}]; 
v2 = Table[Join[nl[[i]], {1}], {i, 1, P}]; 
v3 = Table[Join[nm[[i]], {1}], {i, 1, P}]; 
Q = {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {x2, y2, z2, 1}}; 
U1 = v1.Q; 
U2 = v2.Q;
U3 = v3.Q; 
nk = Table[Take[U1, 3], {i, 1, P}]; 
nl = Table[Take[U2, 3], {i, 1, P}]; 
nm = Table[Take[U3, 3], {i, 1, P}]; 
tplk[[k]]= nk[[2]]; 
tpll[[k]]= nl[[2]];
tplm[[k]]= nm[[2]]),{k,tj}]; 


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