(* Kimball model and energy derived from PDF ES 07 May 2013 : same ERROR as earlier:
C38 has only three bonds. There should be a double bond between C37-C38, my numeration *) 
Clear[tw,n2,n,nb,bc,bh,u,u1,u2,Vee,Vne,Vnn];

t1=1.889726*{{5.814,7.901,-30.790},{3.149,9.127,-32.553},{7.561,7.561,-33.767},{8.734,7.364,-29.079},{4.284,8.853,-27.804},{5.414,8.329,-32.810},{4.176,8.724,-33.336},{4.234,8.885,-34.812},{5.508,8.110,-35.193},{6.231,7.974,-33.857},{5.069,6.778,-35.737},{4.414,10.364,-34.998},{5.765,10.849,-34.445},{5.916,12.319,-34.970},{5.142,12.710,-35.755},{6.849,13.100,-34.365},{7.823,7.602,-31.299},{8.325,7.402,-32.599},{9.709,6.944,-32.504},{10.063,6.898,-31.198},{8.840,7.323,-30.421},{10.619,6.520,-33.613},{11.404,6.476,-30.694},{12.230,5.860,-31.424},{11.980,6.812,-29.317},{6.414,8.128,-28.791},{7.622,7.715,-28.293},{7.607,7.550,-26.801},{6.354,8.377,-26.431},{5.580,8.411,-27.700},{7.295,6.054,-26.600},{6.117,9.115,-25.348},{6.780,8.846,-23.951},{4.151,8.893,-30.252},{3.559,9.101,-29.003},{2.191,9.503,-29.146},{1.987,9.548,-30.540},{3.201,9.181,-31.159},{1.206,9.693,-28.061},{1.006,9.719,-31.567},{-0.232,9.845,-31.472},{1.731,9.571,-32.936},{0.954,8.467,-33.637},{0.490,8.533,-34.647},{0.870,7.311,-32.986},{0.000,6.330,-33.596},{7.224,14.416,-34.971},{8.030,14.064,-36.222},{9.229,14.456,-36.548},{10.000,15.424,-35.643},{9.845,13.847,-37.852},{10.210,14.811,-39.006},{10.445,13.908,-40.267},{11.728,14.401,-40.937},{11.656,15.933,-41.054},{11.893,13.769,-42.326},{13.295,14.105,-42.858},{13.478,13.468,-44.244},{14.922,13.678,-44.734},{15.274,15.178,-44.701},{15.069,13.129,-46.168},{16.549,13.150,-46.600},{16.688,12.603,-48.036},{18.176,12.549,-48.442},{18.306,12.130,-49.919},{18.932,11.550,-47.546},{8.063,7.337,-34.708},{9.632,7.088,-28.527},{3.754,9.033,-26.869},{3.346,8.469,-35.288},{6.100,8.660,-35.924},{5.945,6.194,-36.020},{4.506,6.241,-34.974},{4.437,6.933,-36.612},{3.637,10.934,-34.488},{4.379,10.647,-36.050},{6.503,10.162,-34.858},{5.668,10.782,-33.361},{10.549,7.234,-34.434},{11.646,6.486,-33.248},{10.326,5.531,-33.966},{12.343,7.840,-29.315},{11.203,6.700,-28.561},{12.805,6.136,-29.093},{8.518,7.936,-26.344},{7.262,5.829,-25.534},{6.330,5.820,-27.050},{8.072,5.454,-27.073},{5.046,9.035,-25.159},{6.425,10.129,-25.603},{6.426,9.584,-23.232},{6.512,7.847,-23.607},{7.864,8.920,-24.042},{1.096,10.757,-27.852},{0.244,9.282,-28.367},{1.551,9.181,-27.163},{1.753,10.514,-33.483},{-0.003,5.422,-32.993},{-1.013,6.729,-33.656},{0.359,6.099,-34.599},{7.811,14.925,-34.206},{6.277,14.909,-35.190},{7.425,14.425,-37.054},{8.131,12.979,-36.191},{10.832,15.857,-36.199},{9.333,16.219,-35.309},{10.384,14.885,-34.777},{9.116,13.146,-38.259},{10.770,13.344,-37.570},{11.113,15.333,-38.688},{9.357,15.480,-39.121},{9.577,14.032,-40.914},{10.542,12.882,-39.911},{12.575,14.104,-40.318},{12.563,16.305,-41.530},{10.791,16.212,-41.655},{11.563,16.370,-40.060},{11.131,14.178,-42.990},{11.774,12.689,-42.236},{14.045,13.711,-42.172},{13.401,15.187,-42.935},{12.793,13.925,-44.958},{13.274,12.398,-44.193},{15.615,13.136,-44.090},{16.297,15.319,-45.049},{14.591,15.726,-45.350},{15.183,15.551,-43.681},{14.494,13.737,-46.867},{14.706,12.102,-46.219},{17.143,12.531,-45.927},{16.932,14.170,-46.572},{16.156,13.247,-48.736},{16.271,11.598,-48.097},{18.635,13.529,-48.312},{19.360,12.095,-50.196},{17.861,11.145,-50.058},{17.789,12.854,-50.549},{19.981,11.519,-47.840},{18.856,11.865,-46.505},{18.494,10.558,-47.657}};

MatrixForm[t1]

u={"Mg","C","C","C","C","N","C","C","C","C","C","C","C","C","O","O","N","C","C","C","C","C","C","O","C","N","C","C","C","C","C","C","C","N","C","C","C","C","C","C","O","C","C","O","O","C","C","C","C","C","C","C","C","C","C","C","C","C","C","C","C","C","C","C","C","C","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H","H"};

n=Length[t1]

un=Table[If[u[[i]]=="H","Z",u[[i]]],{i,1,n}]

num1=Table[i,{i,1,n}];

um=Table[{un[[i]],num1[[i]]},{i,1,n}];

us=Sort[um]

num=Table[us[[i,2]],{i,1,n}]

For[i=1,i<n+1,i++,If[us[[i,1]]=="Z",us[[i,1]]="H"]];

u1=Table[If[us[[i,1]]=="C",0.255,If[us[[i,1]]=="H",0.12,If[us[[i,1]]=="N",0.24,If[us[[i,1]]=="O",0.178,If[us[[i,1]]=="Mg",0.5]]]]],{i,1,n}]
u2=Table[If[us[[i,1]]=="C",0.255,If[us[[i,1]]=="H",0.6,If[us[[i,1]]=="N",0.22522,If[us[[i,1]]=="O",0.178,If[us[[i,1]]=="Mg",0.4507]]]]],{i,1,n}]

u=Table[us[[i,1]],{i,1,n}]

cn=Table[If[u[[i]]=="C",6.0,If[u[[i]]=="H",1.0,If[u[[i]]=="N",7.0,If[u[[i]]=="O",8.0,If[u[[i]]=="Mg",2.0]]]]],{i,1,n}]

bc=0;
For[i=1, i<n+1,i++,If[u[[i]]≠"H", bc=bc+1]];
bc
bh=n-bc

DeleteFile[ "C:\\Kimb_Mathem\\tmp3"];
For[i=1,i<n+1,i++,
{u[[i]],Extract[t1,{i}]}>>> "C:\\Kimb_Mathem\\tmp3"]

t=Table[t1[[num[[i]]]],{i,1,n}];
Plot1o=Table[Graphics3D[{Opacity[0.8],If[u[[i]]=="N",Lighter[Blue],If[u[[i]]=="H",Lighter[White],If[u[[i]]=="O",Lighter[Red],If[u[[i]]=="Cl",Lighter[Green],If[u[[i]]=="Mg",Lighter[Orange],Darker[Brown]]]]]],Sphere[{t[[i,1]],t[[i,2]],t[[i,3]]},u1[[i]]]}],{i,1,n}];
Plot2=Table[Graphics3D[Text[i,{t[[i,1]]-0.6,t[[i,2]],t[[i,3]]} ,{1,0}]],{i,1,n}];
Show[Plot1o,Plot2,{AspectRatio->Automatic,Boxed->True, Axes->Automatic,AxesLabel->{xl,yl,zl},SphericalRegion->True}]

(* Distance Matrix Bohr *)
dc=Table[Norm[t[[i]]-t[[j]]],{i,1,n},{j,1,n}];
MatrixForm[Prepend[dc,num1]]

x=Table[t[[i,1]],{i,1,n}];
y=Table[t[[i,2]],{i,1,n}];
z=Table[t[[i,3]],{i,1,n}];

(* connected pairs: distance criteria < 3.0 bohr *)
Dij=Table[If[(dc[[i,j]]<3.0||(u[[i]]=="Cl" && dc[[i,j]]<3.24)),dc[[i,j]],0],{i,1,n},{j,1,n}];
MatrixForm[Prepend[Dij,num1]]

(* count C atoms *)
nc=0;  
For[i=1,i<bc+1,i++,If[u[[i]]=="C",nc=nc+1,nc]]
cnt=Table[0,{nc}];
pri=cnt; sec=cnt;ter=cnt; qua=cnt;
nc

(* count H atoms bonded to each C atom *)
For[j=1,j≤nc,j++,
For[i=bc+1,i≤n,i++,If[ Dij[[i,j]]≠0 ,cnt[[j]]=cnt[[j]]+1,0]]];
cnt

(* parse into prim, sec, tert, and quart C cores with different radii *)
n1=0; n2=0; n3=0; n4=0;
For[i=1,i≤nc,i++,
If[cnt[[i]]==3,n1=n1+1;pri[[n1]]=i,
If[cnt[[i]]==2,n2=n2+1;sec[[n2]]=i,
If[cnt[[i]]==1,n3=n3+1;ter[[n3]]=i,
If[cnt[[i]]=0,n4=n4+1;qua[[n4]]=i]]]]];
For [i=1, i≤nc,i++,If[pri[[i]]≠0,u2[[pri[[i]]]]=0.2572]]
For [i=1, i≤nc,i++,If[sec[[i]]≠0,u2[[sec[[i]]]]=0.2556]]
For [i=1, i≤nc,i++,If[ter[[i]]≠0,u2[[ter[[i]]]]=0.257]]
For [i=1, i≤nc,i++,If[qua[[i]]≠0,u2[[qua[[i]]]]=0.257]]

Eij=Table[If[i≥j,Dij[[i,j]],0],{i,1,nc},{j,1,nc}];
CD=Select[Flatten[Eij],(#≥2.4 &&#≤2.73) &]
Mean[CD] ± StandardDeviation[CD]

CC=Select[Flatten[Eij],(#≥2.73) &]
Mean[CC] ± StandardDeviation[CC]

Eij=Table[If[i≥j,Dij[[i,j]],0],{i,bc+1,n},{j,1,nc}];
CH=Select[Flatten[Eij],(#>0 &&#≤2.1) &]
Mean[CH] ± StandardDeviation[CH]

For[i=1 ,i<n+1,i++,
For[j=1, j<bc+1,j++,
If[(Dij[[i,j]]≠0  && u[[i]]== "H"),If[u[[j]]=="C",u2[[i]]=0.5693,If[u[[j]]=="F",u2[[i]]=0.5,If[u[[j]]=="N",u2[[i]]=0.5541,If[u[[j]]=="O",u2[[i]]=0.5241]]]]]]]
u2

(* subtract core radii from pair distances *)
Fij=Dij-u2;
(*MatrixForm[Fij];*)
Gij=Table[If[i≥j,If[Fij[[i,j]]>0,Fij[[i,j]],0]," "],{i,1,n},{j,1,n}];
(*MatrixForm[Gij];*)

Hij=Transpose[Gij];
(*MatrixForm[Hij];*)

Kij=Hij-u2;
(*MatrixForm[Kij];*)
Lij=Table[If[i≤j,If[Kij[[i,j]]>0,Kij[[i,j]],0]," "],{i,1,n},{j,1,n}];
(*MatrixForm[Lij];*)

(* determine bond cloud radii *)
Lij=Table[If[i≤j,If[Lij[[i,j]]>0,If[u[[j]]≠"H",0.5*Lij[[i,j]],Lij[[i,j]]],0]," "],{i,1,n},{j,1,n}];  
Qij=Lij;
(*MatrixForm[Lij];*)

(* determine kinetic energy *)
Nij=Table[If[i≤j,If[Lij[[i,j]]>0,2.25/Lij[[i,j]]^2,0]," "],{i,1,n},{j,1,n}];
(*MatrixForm[Nij]*)

Tval=Sum[If[i≤j,Nij[[i,j]],0],{i,1,n},{j,1,n}]

Tcor=Sum[If[u[[i]]≠"H",2.25/u2[[i]]^2,0],{i,1,n}]

Ttot=Tval+Tcor  (* plus π and lone pairs: + 8.7616 *)

TvalH=Sum[If[(i≤j && u[[j]]== "H"),Nij[[i,j]],0],{i,1,n},{j,1,n}]

TvalC=Tval-TvalH

Mij=Prepend[Qij,u];
Oij=Transpose[Mij];
Pij=Prepend[Oij,Prepend[u," "]];

con=Table[{i},{i,1,n}]; conj=Table[0,{nc}];
For[i=1,i≤n,i++,
For[j=1,j≤n,j++,If[ Dij[[i,j]]≠0 ,con[[i]]={con[[i]],j}]]];
conn=Table[Flatten[con[[i]]],{i,1,n}];
(* identify C atoms at double bonds *)
j=0;
For[i=1,i≤nc,i++,If[Length[conn[[i]]]==4,j=j+1;conj[[j]]=i]];
conj=Take[conj,{1,j}]
nj=Length[conj];
MatrixForm[conn, TableAlignments->Left]

Clear[c,v,p,Q,P,v1,v2,nk,nl]; (*c1=20; c2=21; c3=51;*) 
rpo=1.4;
rpi=1.4;
tpo=Table[{0,0,0},{nj+8}];
tpu=Table[{0,0,0},{nj+8}];
P=3; c={0,0,0,0,0}; 
(*--------------------------------------------------------------------------*)
(* construct the π clouds for C atoms *)

Do[(c=conn[[conj[[k]]]];
c=Take[c,P];
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];
mk=p;
ml=p;
mk[[1,3]]=rpo;
ml[[1,3]]=-rpi;
U=p;
U1=mk;
U2=ml;
T=RotationMatrix[-T3,{1,0,0}];
p=U.T;
mk=U1.T;
ml=U2.T;
U=p;
U1=mk;
U2=ml;
T=RotationMatrix[-T2,{0,0,1}];
p=U.T;
mk=U1.T;
ml=U2.T;
U=p;
U1=mk;
U2=ml;
T=RotationMatrix[-T1,{1,0,0}];
p=U.T;
mk=U1.T;
ml=U2.T;
Clear[v,v1,v2];
v  =Table[Join[p[[i]],{1}],{i,1,P}];
v1=Table[Join[mk[[i]],{1}],{i,1,P}];
v2=Table[Join[ml[[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;
mk=v1.Q;
ml=v2.Q;
p=Table[Take[p[[i]],3],{i,1,P}];
mk=Table[Take[mk[[i]],3],{i,1,P}];
ml=Table[Take[ml[[i]],3],{i,1,P}];
tpo[[k]]=mk[[1]];
tpu[[k]]=ml[[1]]),{k,nj}];

(*-------------------------------------------------------------------------*)
(* construct the π clouds for the 4 N atoms *)

k=1;
Clear[conn];
conn={{57,5,8},{58,13,16},{59,20,23},{60,27,30}};

Do[(c=conn[[k]];
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=p;
nl=p;
nk[[1,3]]=rpo;
nl[[1,3]]=-rpi;
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}];
tpo[[nj+k]]=nk[[1]];
tpu[[nj+k]]=nl[[1]]),{k,4}];

(*-------------------------------------------------------------------------*)
(* construct the π clouds for the 4 O atoms in C=O *)
k=1;
Clear[conn];
conn={{61,12,11},{63,18,15},{64,32,33},{66,34,35}};

Do[(c=conn[[k]];
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=p;
nl=p;
nk[[1,3]]=rpo;
nl[[1,3]]=-rpi;
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}];
tpo[[nj+4+k]]=nk[[1]];
tpu[[nj+4+k]]=nl[[1]]),{k,4}];

(* ---------------------------------------------------------------------------------------------
O atoms 62,66 LP *)

k=1;
Clear[conn];
conn={{36,62,12},{35,66,34}};
rlo=1.125;
tlo={{0.0,0.0,0.0},{0.0,0.0,0.0}};
tlu={{0.0,0.0,0.0},{0.0,0.0,0.0}};
Do[(c=conn[[k]];
a1=dc[[c[[2]],c[[1]]]]; a2=dc[[c[[2]],c[[3]]]]; d1=dc[[c[[1]],c[[3]]]];
α=ArcCos[(a1^2+a2^2-d1^2)/(2*a1*a2)];
L=(u2[[c[[2]]]]+rlo);
β=ArcSin[rlo/L];
If[k==1,γ=π/2-α/2,γ=3*π/2+α/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=p;
nl=p;
nk[[2,2]]=L*Cos[γ]*Cos[β]; nk[[2,1]]=-L*Sin[γ]*Cos[β];  nk[[2,3]]=rlo; 
nl[[2,2]]=nk[[2,2]]; nl[[2,1]]=nk[[2,1]];  nl[[2,3]]=-rlo;
Plota = Table[Graphics3D[{Opacity[0.5],Sphere[p[[i]],u2[[c[[i]]]]]}],{i,1,P}];
Plotc = Table[Graphics3D[{Tube[{p[[i]],p[[j]]},0.02]}],{i,1,P-1},{j,i+1,P}];
Plotd=Graphics3D[{Opacity[0.5],Lighter[Red],Sphere[nk[[2]],rlo]}];
Plote=Graphics3D[{Opacity[0.5],Lighter[Red],Sphere[nl[[2]],rlo]}];
Plotb0=Table[Graphics3D[Text[c[[i]],p[[i]],{1,0}]],{i,1,P}];
Show[Plotd,Plote,Plota,Plotb0,(*Plotb1,Plotb2,*)Plotc,{Axes->Automatic,AxesLabel->{xl,yl,zl},AspectRatio->Automatic,Boxed->False,PlotRange->Automatic, SphericalRegion->True}]
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}];
tlo[[k]]=nk[[2]];
tlu[[k]]=nl[[2]]),{k,2}];

(*-------------------------------------------------------------------------*)
(* next: N57, N58, N59, N60 one doubly occ LP as sp2 ligand *)
k=1;
Clear[conn];
tpln={{0.0,0.0,0.0},{0.0,0.0,0.0},{0.0,0.0,0.0},{0.0,0.0,0.0}};
conn={{5,57,8},{13,58,16},{20,59,23},{27,60,30}};
rln=1.246;
Do[(c=conn[[k]];
a1=dc[[c[[2]],c[[1]]]]; a2=dc[[c[[2]],c[[3]]]]; d1=dc[[c[[1]],c[[3]]]];
α=ArcCos[(a1^2+a2^2-d1^2)/(2*a1*a2)];
β=ArcSin[rln/(u2[[c[[2]]]]+rln)];
If[(k==1 || k==3),γ=π-α/4,γ=3*π/2+α/2];
L=(u2[[c[[2]]]]+rln);
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=p;
nk[[2,2]]=L*Cos[γ]; nk[[2,1]]=-L*Sin[γ];  nk[[2,3]]=0; 
(*Plota = Table[Graphics3D[{Opacity[0.5],Sphere[p[[i]],u2[[c[[i]]]]]}],{i,1,P}];
Plotc = Table[Graphics3D[{Tube[{p[[i]],p[[j]]},0.02]}],{i,1,P-1},{j,i+1,P}];
Plotd=Graphics3D[{Opacity[0.5],Lighter[Blue],Sphere[nk[[2]],rln]}];
Plotb0=Table[Graphics3D[Text[c[[i]],p[[i]],{1,0}]],{i,1,P}];
Show[Plotd,(*Plote,*)Plota,Plotb0,(*Plotb1,Plotb2,*)Plotc,{Axes->Automatic,AxesLabel->{xl,yl,zl},AspectRatio->Automatic,Boxed->False,PlotRange->All, SphericalRegion->True}]*)
U=p;
U1=nk;
T=RotationMatrix[-T3,{1,0,0}];
p=U.T;
nk=U1.T;
U=p;
U1=nk;
T=RotationMatrix[-T2,{0,0,1}];
p=U.T;
nk=U1.T;
U=p;
U1=nk;
T=RotationMatrix[-T1,{1,0,0}];
p=U.T;
nk=U1.T;
Clear[v,v1,v2];
v  =Table[Join[p[[i]],{1}],{i,1,P}];
v1=Table[Join[nk[[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;
p=Table[Take[p[[i]],3],{i,1,P}];
nk=Table[Take[nk[[i]],3],{i,1,P}];
tpln[[k]]=nk[[2]]),{k,4}];

DeleteFile[ "C:\\Kimb_Mathem\\tmp4"];
For[i=2,i<n+2,i++,
For[j=2,j≤i,j++,
If[Pij[[i,j]]>0,{i-1,j-1,Extract[Pij,{i,j}]}>>> "C:\\Kimb_Mathem\\tmp4"]]]

a=ReadList["C:\\Kimb_Mathem\\tmp4"]
nb=Length[a]
cw=nb-bh

wxc=Table[(a[[i,3]]+u2[[a[[i,1]]]])/(u2[[a[[i,1]]]]+2*a[[i,3]]+u2[[a[[i,2]]]])*(x[[a[[i,2]]]]-x[[a[[i,1]]]])+x[[a[[i,1]]]],{i,1,cw}];
wyc=Table[(a[[i,3]]+u2[[a[[i,1]]]])/(u2[[a[[i,1]]]]+2*a[[i,3]]+u2[[a[[i,2]]]])*(y[[a[[i,2]]]]-y[[a[[i,1]]]])+y[[a[[i,1]]]],{i,1,cw}];   
wzc=Table[(a[[i,3]]+u2[[a[[i,1]]]])/(u2[[a[[i,1]]]]+2*a[[i,3]]+u2[[a[[i,2]]]])*(z[[a[[i,2]]]]-z[[a[[i,1]]]])+z[[a[[i,1]]]],{i,1,cw}];

Plot7=Table[Graphics3D[{Darker[Brown,0.4],Opacity[0.5],If[wxc[[i]]≠0,Sphere[{wxc[[i]],wyc[[i]],wzc[[i]]},a[[i,3]]]]}],{i,1,cw}];

wxh=Table[(a[[i,3]]+u2[[a[[i,2]]]])/(u2[[a[[i,1]]]]+a[[i,3]]+u2[[a[[i,2]]]])*(x[[a[[i,1]]]]-x[[a[[i,2]]]])+x[[a[[i,2]]]],{i,cw+1,nb}];
wyh=Table[(a[[i,3]]+u2[[a[[i,2]]]])/(u2[[a[[i,1]]]]+a[[i,3]]+u2[[a[[i,2]]]])*(y[[a[[i,1]]]]-y[[a[[i,2]]]])+y[[a[[i,2]]]],{i,cw+1,nb}];
wzh=Table[(a[[i,3]]+u2[[a[[i,2]]]])/(u2[[a[[i,1]]]]+a[[i,3]]+u2[[a[[i,2]]]])*(z[[a[[i,1]]]]-z[[a[[i,2]]]])+z[[a[[i,2]]]],{i,cw+1,nb}];

Plot8=Table[Graphics3D[{Opacity[0.5],Sphere[{wxh[[i]],wyh[[i]],wzh[[i]]},a[[i+cw,3]]]}],{i,1,bh}];

Plot8a=Table[Graphics3D[{Opacity[0.3],Lighter[Cyan],Sphere[tpo[[i]],rpo]}],{i,1,nj+8}];
Plot8b=Table[Graphics3D[{Opacity[0.3],Lighter[Yellow],Sphere[tpu[[i]],rpi]}],{i,1,nj+8}];
Plot8c=Table[Graphics3D[{Opacity[0.4],Lighter[Red],Sphere[tlo[[i]],rlo]}],{i,1,2}];
Plot8d=Table[Graphics3D[{Opacity[0.4],Lighter[Red],Sphere[tlu[[i]],rlo]}],{i,1,2}];
Plot8e=Table[Graphics3D[{Opacity[0.4],Lighter[Blue],Sphere[tpln[[i]],rln]}],{i,1,4}];

Plot9= Table[Graphics3D[{Lighter[Orange,0.5],Tube[{{x[[a[[i,1]]]],y[[a[[i,1]]]],z[[a[[i,1]]]]},{x[[a[[i,2]]]],y[[a[[i,2]]]],z[[a[[i,2]]]]}},0.1]}],{i,1,nb}];
(*Plot9a= Graphics3D[Tube[{{x[[21]],y[[21]],z[[21]]},{nk[[2,1]],nk[[2,2]],nk[[2,3]]}},0.06]];
Plot9b= Graphics3D[Tube[{{x[[21]],y[[21]],z[[21]]},{nl[[2,1]],nl[[2,2]],nl[[2,3]]}},0.06]];*)

Show[Plot1o,(*Plot2,*)Plot7,Plot8,Plot8a,Plot8b,Plot8c,Plot8d,Plot8e,Plot9,{AspectRatio->Automatic,Boxed->False,(* Axes->Automatic,AxesLabel->{xl,yl,zl},*)SphericalRegion->True}]

Show[Plot1o,Plot2,Plot9,(*Plot9a,Plot9b,*){Axes->Automatic,AxesLabel->{xl,yl,zl},AspectRatio->Automatic,Boxed->True, SphericalRegion->True}]

tc=Table[{x[[i]],y[[i]],z[[i]]},{i,1,bc}];   (* core clouds *)
twc=Table[{wxc[[i]],wyc[[i]],wzc[[i]]},{i,1,cw}];  (* Bond cöouds between cores *)
twh=Table[{wxh[[i]],wyh[[i]],wzh[[i]]},{i,1,bh}]; (* H-bonding clouds *) 
tw=Join[tc,twc,twh,tpo,tpu,tlo,tlu,tpln];
bt=Length[tw]
(*MatrixForm[tw]*)
(* coordinates of clouds *)

(* cloud-cloud distances *)
dw=Table[Norm[tw[[i]]-tw[[j]]],{i,1,bt},{j,1,bt}];
(*MatrixForm[dw]*)
Length[dw]

(* nucleus-nucleus distances *)
(*MatrixForm[dc] *)

tv1=Table[t[[i]],{i,1,bc}];
tv2=Table[{0,0,0},{i,1,cw}];
tv3= Table[t[[i]],{i,bc+1,n}];
tv4 =Table[{0,0,0},{i,1,nj+4}]; 
tv5 =Table[{0,0,0},{i,1,nj+4}]; 
tv6 =Table[{0,0,0},{i,1,16}]; 
tv  = Join[tv1,tv2,tv3,tv4,tv5,tv6];
MatrixForm[tv]
Length[tv]

ow=Join[Table[-2.0,{i,1,bc+nb}],Table[-0.5,{i,1,2*nj}],Table[-0.5,{i,1,6}],Table[-1.0,{i,1,4}],Table[-2.0,{i,1,16}]];
lw=Length[ow]
cnv=Join[Table[cn[[i]],{i,1,bc}],Table[0.0,{i,1,cw}],Table[cn[[i]],{i,bc+1,n}],Table[0.0,{i,1,2*(nj+5)}],Table[0.0,{i,1,16}]];
lc=Length[cnv]

echarg=Sum[ow[[i]],{i,1,lw}]
pcharg=Sum[cnv[[i]],{i,1,lc}]

(* nucleus-cloud distances *)
dcw=Table[sq=Norm[tw[[i]]-tv[[j]]];If[(sq> 0.1 && i≠j),1/sq,0],{i,1,bt},{j,1,bt}];
Length[dcw]
(*MatrixForm[dcw]*)

Vee=0.0; Vnn=0.0;

(* all clouds repulsion *)
For[i=1,i<bt,i++,
For[j=i+1,j<bt+1,j++,
Vee=Vee+ow[[i]]*ow[[j]]/dw[[i,j]]]]
Vee

(* half cloud with other common halfcloud: subtract WW *)
Vee=Vee-2/(rpo+rpi)]; (* 5*1/4 C + 3*1/4 N *)
Vee

sc=0.372;   (* core clouds *)
For[i=1 ,i<bc+1,i++,Vee=Vee+3*sc/u2[[i]]]
sw=0.378;  (* bonding clouds *)
For[i=1 ,i<nb+1,i++,Vee=Vee+3*sw/a[[i,3]]]
For[i=1 ,i<9,i++,Vee=Vee+3*sw/rlo]
For[i=1 ,i<4,i++,Vee=Vee+3*sw/rln]
Vee

(* all nuclei with all clouds, except i = j *)
Vne=cnv.(ow.dcw)

(* core nuclei in their cloud centers *)

For[i=1,i<bc+1,i++,Vne=Vne-3*cn[[i]]/u2[[i]]]
Vne

(* protons in their bonding cloud, distance u2[[a[[i,1]]]] to cloud center *)
For[i=cw+1,i<nb+1,i++,Vne=Vne-(3-(u2[[a[[i,1]]]]/a[[i,3]])^2)/a[[i,3]]]
Vne

(* nucleus-nucleus *)
For[i=1,i<n,i++,
For[j=i+1,j<n+1,j++,Vnn=Vnn+cn[[i]]*cn[[j]]/dc[[i,j]]]]
Vnn

Vee
Vnn
Vne
Tpi   = 2.296*(nj+3)*0.5*1.125*(1/rpo^2+1/rpi^2) (* pi kin e added *)
Tpi   = Tpi+ 2.296*4*1.0*1.125*(1/rpo^2+1/rpi^2) (* pi kin e added *)
Tlo   = 8*2*1.125*(1/rlo^2) (* lp kin e of 4 O added *)
Tln   = 3*2*1.125*(1/rln^2) (* lp kin e of 3 N added *)
Vtot=Vne+Vee+Vnn
Ttot=Ttot+Tpi+Tlo+Tln
Vir   =-Vtot/Ttot
Etot =Ttot+Vtot