(* Kimball model and energy derived from Facio-PDB source, ES 10 May 2013 *) 
(* Graphite stack of 4 layers 50 atoms each, dangling bonds on the outside *)
Clear[tw,n1,n2,n3,n,nb,bc,bh,u,u1,u2,Vee,Vne,Vnn];

t1=1.889726*{{1.339,5.378,-4.149},{-1.079,1.746,-3.946},{1.441,1.816,-3.988},{0.149,3.619,-2.788},{-1.023,-0.984,-5.086},{0.206,0.890,-3.927},{-0.977,-1.816,-3.786},{1.543,-1.746,-3.828},{0.251,0.057,-2.628},{-0.921,-4.546,-4.925},{0.308,-2.672,-3.766},{-0.876,-5.378,-3.626},{0.353,-3.505,-2.468},{-1.143,5.422,-1.587},{-3.607,2.623,-2.684},{-1.085,2.693,-2.726},{-2.378,4.496,-1.526},{-3.562,1.790,-1.385},{-1.041,1.860,-1.427},{-2.333,3.663,-0.227},{-2.372,3.548,-2.746},{-2.315,0.820,-3.885},{-3.505,-0.939,-2.524},{-0.984,-0.869,-2.566},{-2.276,0.934,-1.365},{-3.460,-1.772,-1.225},{-0.939,-1.702,-1.267},{-2.231,0.101,-0.067},{-2.270,-0.014,-2.585},{-2.213,-2.742,-3.725},{-3.403,-4.501,-2.364},{-0.882,-4.431,-2.406},{-2.174,-2.628,-1.205},{-3.358,-5.334,-1.065},{-0.838,-5.264,-1.107},{-2.129,-3.461,0.094},{-2.168,-3.575,-2.425},{-3.625,5.466,0.974},{-6.089,2.667,-0.123},{-3.568,2.737,-0.165},{-4.860,4.541,1.036},{-6.044,1.834,1.176},{-3.524,1.904,1.134},{-4.815,3.707,2.334},{-4.854,3.593,-0.184},{-4.797,0.864,-1.324},{-5.987,-0.895,0.037},{-3.466,-0.825,-0.005},{-4.758,0.979,1.196},{-5.942,-1.728,1.336},{-3.422,-1.657,1.294},{-4.713,0.145,2.495},{-4.752,0.031,-0.024},{-4.695,-2.698,-1.164},{-3.364,-4.387,0.155},{-4.656,-2.583,1.356},{-3.320,-5.219,1.455},{-4.612,-3.417,2.655},{-4.650,-3.531,0.136},{-6.050,2.781,2.396},{-6.006,1.949,3.695},{-5.948,-0.781,2.556},{-5.904,-1.613,3.856},{1.379,5.492,-1.629},{3.899,5.562,-1.671},{2.625,4.522,-4.129},{1.435,2.763,-2.768},{3.956,2.833,-2.810},{2.664,4.637,-1.609},{1.480,1.930,-1.469},{4.001,2.001,-1.511},{2.709,3.803,-0.311},{2.670,3.689,-2.829},{2.727,0.960,-3.969},{1.537,-0.799,-2.608},{4.058,-0.729,-2.650},{2.766,1.075,-1.449},{1.582,-1.632,-1.309},{4.102,-1.561,-1.351},{2.811,0.241,-0.150},{2.772,0.127,-2.669},{2.829,-2.602,-3.809},{1.638,-4.361,-2.448},{2.868,-2.487,-1.289},{1.684,-5.193,-1.149},{2.913,-3.321,0.010},{2.874,-3.435,-2.509},{-1.104,5.537,0.932},{1.417,5.607,0.890},{0.143,4.566,-1.568},{-1.047,2.807,-0.207},{1.474,2.877,-0.249},{0.182,4.681,0.952},{-1.002,1.975,1.092},{1.518,2.045,1.050},{0.227,3.847,2.251},{0.188,3.733,-0.268},{0.245,1.005,-1.408},{-0.946,-0.755,-0.047},{1.576,-0.685,-0.089},{0.284,1.119,1.112},{-0.900,-1.587,1.253},{1.620,-1.517,1.211},{0.328,0.285,2.411},{0.290,0.171,-0.108},{0.347,-2.557,-1.248},{-0.844,-4.317,0.113},{1.678,-4.247,0.071},{0.386,-2.443,1.272},{-0.798,-5.149,1.413},{1.722,-5.079,1.371},{0.430,-3.276,2.571},{0.392,-3.391,0.052},{-3.586,5.581,3.493},{-1.066,5.651,3.451},{-2.339,4.611,0.993},{-3.530,2.851,2.354},{-1.008,2.922,2.312},{-2.300,4.725,3.513},{-3.484,2.019,3.654},{-0.964,2.089,3.612},{-2.256,3.892,4.812},{-2.294,3.777,2.293},{-2.237,1.049,1.153},{-3.428,-0.710,2.514},{-0.906,-0.640,2.473},{-2.199,1.163,3.673},{-3.382,-1.543,3.814},{-0.862,-1.473,3.772},{-2.154,0.330,4.972},{-2.192,0.215,2.453},{-2.135,-2.513,1.314},{-3.326,-4.272,2.675},{-0.804,-4.202,2.633},{-2.097,-2.399,3.834},{-3.280,-5.105,3.974},{-0.760,-5.035,3.932},{-2.052,-3.232,5.132},{-2.091,-3.347,2.614},{-3.548,5.695,6.013},{-4.821,4.655,3.554},{-6.012,2.896,4.915},{-3.490,2.966,4.874},{-4.783,4.770,6.074},{-5.966,2.063,6.215},{-3.446,2.133,6.173},{-4.738,3.936,7.373},{-4.776,3.822,4.854},{-4.719,1.093,3.715},{-5.910,-0.666,5.076},{-3.388,-0.596,5.034},{-4.681,1.208,6.235},{-5.864,-1.499,6.375},{-3.344,-1.429,6.333},{-4.636,0.374,7.533},{-4.675,0.260,5.015},{-4.618,-2.469,3.875},{-3.287,-4.158,5.194},{-4.579,-2.354,6.395},{-3.242,-4.991,6.493},{-4.534,-3.188,7.694},{-4.573,-3.302,5.175},{3.938,5.677,0.848},{5.185,4.707,-1.652},{3.994,2.948,-0.291},{6.516,3.018,-0.333},{5.224,4.821,0.868},{4.040,2.115,1.008},{6.560,2.185,0.967},{5.269,3.988,2.167},{5.230,3.873,-0.352},{5.287,1.145,-1.492},{4.096,-0.614,-0.131},{6.618,-0.544,-0.173},{5.326,1.259,1.028},{4.142,-1.447,1.169},{6.662,-1.377,1.127},{5.370,0.426,2.327},{5.332,0.311,-0.192},{5.389,-2.417,-1.332},{4.198,-4.176,0.030},{5.427,-2.303,1.188},{4.244,-5.009,1.329},{5.472,-3.136,2.487},{5.433,-3.251,-0.031},{1.456,5.721,3.409},{3.976,5.791,3.368},{2.703,4.751,0.909},{1.512,2.992,2.270},{4.034,3.062,2.228},{2.742,4.866,3.429},{1.558,2.159,3.570},{4.078,2.229,3.528},{2.786,4.032,4.728},{2.748,3.918,2.209},{2.805,1.189,1.069},{1.614,-0.570,2.431},{4.136,-0.500,2.389},{2.843,1.304,3.590},{1.660,-1.403,3.730},{4.180,-1.333,3.688},{2.888,0.470,4.888},{2.849,0.356,2.370},{2.906,-2.373,1.230},{1.716,-4.132,2.591},{4.238,-4.062,2.549},{2.945,-2.258,3.750},{1.761,-4.965,3.890},{4.282,-4.894,3.848},{2.990,-3.092,5.048},{2.951,-3.206,2.530},{-1.026,5.765,5.971},{1.494,5.836,5.929},{0.221,4.795,3.470},{-0.970,3.036,4.832},{1.552,3.106,4.790},{0.259,4.910,5.991},{-0.924,2.204,6.131},{1.596,2.274,6.089},{0.304,4.076,7.289},{0.265,3.962,4.771},{0.322,1.233,3.631},{-0.868,-0.526,4.992},{1.653,-0.456,4.950},{0.361,1.348,6.151},{-0.823,-1.358,6.291},{1.698,-1.288,6.249},{0.406,0.514,7.449},{0.367,0.400,4.931},{0.424,-2.328,3.791},{-0.766,-4.088,5.152},{1.755,-4.018,5.110},{0.463,-2.214,6.311},{-0.721,-4.920,6.451},{1.799,-4.850,6.410},{0.508,-3.048,7.610},{0.469,-3.162,5.091},{-0.988,5.880,8.490},{-2.262,4.840,6.032},{-3.452,3.080,7.393},{-0.931,3.151,7.351},{-2.223,4.954,8.552},{-3.407,2.248,8.692},{-0.886,2.318,8.650},{-2.178,4.121,9.851},{-2.217,4.006,7.332},{-2.160,1.278,6.192},{-3.350,-0.482,7.553},{-0.829,-0.411,7.511},{-2.121,1.392,8.712},{-3.305,-1.314,8.852},{-0.785,-1.244,8.811},{-2.076,0.559,10.011},{-2.115,0.444,7.492},{-2.058,-2.284,6.352},{-3.249,-4.044,7.713},{-0.727,-3.973,7.671},{-2.019,-2.170,8.872},{-3.203,-4.876,9.013},{-0.683,-4.806,8.971},{-1.974,-3.003,10.171},{-2.013,-3.118,7.652},{6.554,3.132,2.187},{6.600,2.300,3.486},{6.656,-0.430,2.347},{6.701,-1.262,3.646},{4.016,5.906,5.887},{5.262,4.936,3.387},{4.072,3.176,4.748},{6.594,3.247,4.706},{5.301,5.050,5.907},{4.117,2.344,6.047},{6.638,2.414,6.005},{5.346,4.217,7.205},{5.307,4.102,4.687},{5.364,1.374,3.547},{4.174,-0.386,4.908},{6.695,-0.315,4.866},{5.403,1.488,6.067},{4.219,-1.218,6.207},{6.739,-1.148,6.165},{5.448,0.655,7.366},{5.409,0.540,4.847},{5.466,-2.188,3.707},{4.276,-3.947,5.068},{5.505,-2.074,6.227},{4.321,-4.780,6.368},{5.550,-2.907,7.526},{5.511,-3.022,5.007},{1.533,5.950,8.448},{2.780,4.980,5.948},{1.590,3.221,7.309},{4.111,3.291,7.267},{2.819,5.094,8.468},{1.635,2.388,8.608},{4.155,2.458,8.566},{2.864,4.261,9.767},{2.825,4.146,7.248},{2.882,1.418,6.108},{1.691,-0.341,7.469},{4.213,-0.271,7.427},{2.921,1.533,8.628},{1.737,-1.174,8.769},{4.257,-1.104,8.727},{2.966,0.699,9.927},{2.927,0.585,7.408},{2.984,-2.144,6.268},{1.793,-3.903,7.630},{4.315,-3.833,7.588},{3.023,-2.029,8.788},{1.839,-4.736,8.929},{3.067,-2.863,10.087},{3.029,-2.977,7.568},{0.298,5.024,8.509},{-0.893,3.265,9.870},{1.629,3.335,9.828},{-0.847,2.432,11.170},{1.673,2.503,11.128},{0.343,4.191,9.809},{0.400,1.462,8.669},{-0.791,-0.297,10.031},{1.731,-0.227,9.989},{0.439,1.577,11.189},{-0.745,-1.130,11.330},{1.775,-1.059,11.288},{0.445,0.629,9.970},{0.502,-2.100,8.830},{-0.689,-3.859,10.191},{1.833,-3.789,10.149},{0.540,-1.985,11.350},{0.546,-2.933,10.130}};


u=Table["C",{i,1,331}];

n=Length[t1]

u1=Table[0.255,{i,1,n}]
u2=Table[0.25384,{i,1,n}]

cn=Table[6.0,{i,1,n}]

bc=n;
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=t1;

x=Table[t[[i,1]],{i,1,n}];
y=Table[t[[i,2]],{i,1,n}];
z=Table[t[[i,3]],{i,1,n}];
(*m=Table[If[u[[i]]=="C",12.0107,If[u[[i]]=="H",1.00794]],{i,1,n}];
M=Sum[m[[i]],{i,n}]*)

(*coords=(m.t)/M*)
(*no=Table[1,{i,1,200}];
coords=(no.t)/n
x=Table[t[[i,1]]-coords[[1]],{i,1,n}]; 
y=Table[t[[i,2]]-coords[[2]],{i,1,n}];
z=Table[t[[i,3]]-coords[[3]],{i,1,n}];
t=Table[{x[[i]],y[[i]],z[[i]]},{i,1,n}];*)

(*w3=Sqrt[3.0]/2; ws=0.5;
x=Table[t[[i,1]]*w3-ws*t[[i,2]],{i,1,n}]; 
y=Table[t[[i,1]]*ws+w3*t[[i,2]]-2.8727498,{i,1,n}];
t=Table[{x[[i]],y[[i]],z[[i]]},{i,1,n}];
MatrixForm[t]*)

(*coords=(m.t)/M*)

dc=Table[Norm[t[[i]]-t[[j]]],{i,1,n},{j,1,n}];

(* Distance Matrix Bohr *)
Dij=Table[If[(dc[[i,j]]<3.0 ||(u[[i]]=="Cl" && j==10 && i== 17)),dc[[i,j]],0],{i,1,n},{j,1,n}];
(*MatrixForm[Dij]*)

For[i=1 ,i<n+1,i++,
For[j=1, j<i,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

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];*)

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];*)

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,{bc}];
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 *)
con1=con2=con3=con4=Table[0,{n}];

j=0;
For[i=1,i≤n,i++,If[Length[conn[[i]]]==2(* && u[[conn[[i,4]]]]≠"O")*),j=j+1;con1[[j]]=i]];
con1=Take[con1,{1,j}]
n1=Length[con1]

j=0;
For[i=1,i≤n,i++,If[Length[conn[[i]]]==3(* && u[[conn[[i,4]]]]≠"O")*),j=j+1;con2[[j]]=i]];
con2=Take[con2,{1,j}]
n2=Length[con2]

j=0;
For[i=1,i≤n,i++,If[Length[conn[[i]]]==4(* && u[[conn[[i,4]]]]≠"O")*),j=j+1;con3[[j]]=i]];
con3=Take[con3,{1,j}]
n3=Length[con3]

j=0;
For[i=1,i≤n,i++,If[Length[conn[[i]]]==5(* && u[[conn[[i,4]]]]≠"O")*),j=j+1;con4[[j]]=i]];
con4=Take[con4,{1,j}]
n4=Length[con4]

MatrixForm[conn, TableAlignments->Left]

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

data1=Table[a[[i,3]],{i,1,nb}];
Histogram[data1,{1.20,1.21,0.0005}]

rcc=1.093;  L=u2[[1]]+rcc;
db1=Table[{0,0,0},{i,1,n1}];
db2=Table[{0,0,0},{i,1,n2}];
db3=Table[{0,0,0},{i,1,n3}];
For[i=1,i≤n1,i++,db1[[i]]=t[[con1[[i]]]]];
For[i=1,i≤n2,i++,db2[[i]]=t[[con2[[i]]]]];
For[i=1,i≤n3,i++,db3[[i]]=t[[con3[[i]]]]];

(*onebond, two dangling bonds with one electron, each *) (*od={1,41,81,121,161,5,45,85,125,165,34,74,114,154,194,38,78,118,158,198,3,11,19,27,35,7,15,23,31,39,164,172,180,188,196,168,176,184,192,200};
(*twobond, one dangling bond with one elctron *) 
td={9,17,25,33,13,21,29,37,162,170,178,186,166,174,182,190,43,83,123,163,47,87,127,167,36,76,116,156,40,80,120,160};
rcc=1.093; L=0.25514+rcc; Lc=L*Cos[0.5236]; Ls=0.5*L; 
db11=Table[{0,0,0},{i,1,40}];
db12=Table[{0,0,0},{i,1,40}];
db2=Table[{0,0,0},{i,1,32}];
For[i=1,  i≤10,i++, db11[[i,1]]=t[[od[[i]],1]]-Lc;  db11[[i,2]]=t[[od[[i]],2]]+Ls;    db11[[i,3]]=t[[od[[i]],3]]];
For[i=11,i≤20,i++,db11[[i,1]]=t[[od[[i]],1]]-Lc;   db11[[i,2]]=t[[od[[i]],2]]-Ls;    db11[[i,3]]=t[[od[[i]],3]]];
For[i=21,i≤30,i++,db11[[i,1]]=t[[od[[i]],1]];            db11[[i,2]]=t[[od[[i]],2]]+L;      db11[[i,3]]=t[[od[[i]],3]]];
For[i=31,i≤40,i++,db11[[i,1]]=t[[od[[i]],1]]+Lc;   db11[[i,2]]=t[[od[[i]],2]]+Ls;    db11[[i,3]]=t[[od[[i]],3]]];

For[i=1,  i≤10,i++, db12[[i,1]]=t[[od[[i]],1]]+Lc;  db12[[i,2]]=t[[od[[i]],2]]+Ls;  db12[[i,3]]=t[[od[[i]],3]]];
For[i=11,i≤20,i++,db12[[i,1]]=t[[od[[i]],1]]+Lc;  db12[[i,2]]=t[[od[[i]],2]]-Ls;  db12[[i,3]]=t[[od[[i]],3]]];
For[i=21,i≤30,i++,db12[[i,1]]=t[[od[[i]],1]]-Lc;  db12[[i,2]]=t[[od[[i]],2]]-Ls;  db12[[i,3]]=t[[od[[i]],3]]];
For[i=31,i≤40,i++,db12[[i,1]]=t[[od[[i]],1]];           db12[[i,2]]=t[[od[[i]],2]]-L;     db12[[i,3]]=t[[od[[i]],3]]];

For[i=1,i≤8,    i++,db2[[i,1]]=t[[td[[i]],1]]-Lc;db2[[i,2]]=t[[td[[i]],2]]+Ls;db2[[i,3]]=t[[td[[i]],3]]];
For[i=9,i≤16,  i++,db2[[i,1]]=t[[td[[i]],1]]+Lc;db2[[i,2]]=t[[td[[i]],2]]-Ls;db2[[i,3]]=t[[td[[i]],3]]];
For[i=17,i≤24,i++,db2[[i,1]]=t[[td[[i]],1]];         db2[[i,2]]=t[[td[[i]],2]]+L;  db2[[i,3]]=t[[td[[i]],3]]];
For[i=25,i≤32,i++,db2[[i,1]]=t[[td[[i]],1]];         db2[[i,2]]=t[[td[[i]],2]]-L;  db2[[i,3]]=t[[td[[i]],3]]];
*)

Plot1=Table[Graphics3D[{If[u[[i]]=="N",Lighter[Blue],If[u[[i]]=="O",Lighter[Red],If[u[[i]]=="Cl",Lighter[Green],If[u[[i]]=="B",Lighter[Orange],Lighter[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}];

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[Orange,0.5],{Opacity[0.4],If[wxc[[i]]≠0,Sphere[{wxc[[i]],wyc[[i]],wzc[[i]]},a[[i,3]]]]}}],{i,1,cw}];

Plot8c=Table[Graphics3D[{Darker[Orange,0.3],{Opacity[0.2],Sphere[db1[[i]],rcc]}}],{i,1,n1}];
Plot8d=Table[Graphics3D[{Darker[Orange,0.3],{Opacity[0.2],Sphere[db2[[i]],rcc]}}],{i,1,n2}];
Plot8e=Table[Graphics3D[{Darker[Orange,0.3],{Opacity[0.2],Sphere[db3[[i]],rcc]}}],{i,1,n3}];

Plot9= Table[Graphics3D[{(*Darker[Red,0.8],*)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.15]}],{i,1,nb}];

(*Plot10=Graphics3D[{Darker[Red,0.5],Sphere[{coords[[1]],coords[[2]],coords[[3]]},0.2]}];*)

Show[Plot1,Plot7,(*Plot8,*)(*Plot8a,Plot8b,*)Plot8c,Plot8d,Plot8e,Plot9,{AspectRatio->Automatic,Boxed->False, SphericalRegion->True}]
Show[Plot1,Plot2,Plot8c,Plot8d,Plot8e,Plot9,(*Plot10,*){AspectRatio->Automatic,Axes->Automatic,Boxed->True, AxesLabel -> {xl, yl, zl}, 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 *)
(*tpo=Table[{px[[i]],py[[i]],pz[[i]]},{i,1,bc}]; (* pi clouds above *)
tpu=Table[{qx[[i]],qy[[i]],qz[[i]]},{i,1,bc}]; (* pi clouds below *) *)
tw=Join[tc,twc,db1,db2,db3];
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 *)
Length[dc]

(*MatrixForm[dc] *)

tv1=Table[{x[[i]],y[[i]],z[[i]]},{i,1,bc}];
tv2=Table[{0,0,0},{i,1,cw}];
(*tv4 =Table[{0,0,0},{i,1,bc}];*) 
tv3=Table[{0,0,0},{i,1,n1}]; 
tv4=Table[{0,0,0},{i,1,n2}]; 
tv5=Table[{0,0,0},{i,1,n3}]; 
tv  = Join[tv1,tv2,tv3,tv4,tv5];
Length[tv]
(*MatrixForm[tv]*)
(*For[i=1,i≤40,i++,cn[[onebond[[i]]]]=4.0];
For[i=1,i≤32,i++,cn[[twobond[[i]]]]=5.0];*)

ow=Join[Table[-2.0,{i,1,bc+nb}],Table[-3.0,{i,1,n1}],Table[-2.0,{i,1,n2}],Table[-1.0,{i,1,n3}]]
cnv=Join[Table[cn[[i]],{i,1,bc}],Table[0.0,{i,1,cw}],Table[0.0,{i,1,n1+n2+n3}]];

(* 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]

(*cnv=Join[Table[cn[[i]],{i,1,bc}],Table[0.0,{i,1,cw}],Table[cn[[i]],{i,bc+1,n}]]*)
Elcharge=Sum[ow[[i]],{i,1,Length[ow]}]
Length[ow]
Nuccharge=Sum[cnv[[i]],{i,1,Length[cnv]}]
Length[cnv]

Vee=0.0; Vnn=0.0;

(* all clouds repulsion *)
For[i=1,i<bt,i++,
For[j=i+1,j≤bt,j++,
If[dw[[i,j]]≠ 0,Vee=Vee+ow[[i]]*ow[[j]]/dw[[i,j]],Vee]]]
Vee=Vee/1.0027
Vee
(* half cloud with other common halfcloud: subtract WW *)
(*For [i=nb+bc+1,i<nb+2*bc+1,i++,Vee=Vee-ow[[i]]*ow[[i+bc]]/dw[[i,i+bc]]];*)
(*Vee=Vee-bc*0.5*0.25/rpo;
Vee*0

sc=0.385;   (* core clouds *)
For[i=1 ,i<bc+1,i++,Vee=Vee+3*sc/u2[[i]]]
sw=0.385;  (* bonding clouds *)
For[i=1 ,i<nb+1,i++,Vee=Vee+3*sw/a[[i,3]]]
Vee

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

(*Vne=Vne-(2*n1+n2)*4.0/L *)(* "ficticious" interaction with own core *)
Vne=Vne/0.999146

(* 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,j++,Vnn=Vnn+cn[[i]]*cn[[j]]/dc[[i,j]]]]
Vnn

Vee
Vnn
Vne
(*Tpi =2.15*bc*0.5*1.125*(1/rpo^2+1/rpi^2)*) (* pi kin e added *)
Tdb = (3*n1+2*n2+n3)*1.125/rcc^2
Vtot=Vne+Vee+Vnn
Ttot =Ttot+Tdb
Vir    = -Vtot/Ttot
Etot=Ttot+Vtot