(*H2O with "orthogonal" σ,π LPs, 12.01.2013, not final optimization: slight overlaps *)
Clear[z,sig1,sig2,sig3,k1,k2,k3,k4,nc,R1,R2,R3,R4,w,w1,p,vee,vne,vnn,
xn,yn,zn,xc,yc,zc,c1,c2,c3,cs,ss,d1,d2,d3,d4,d5,d6,pi,i,j,d,t];
z=8.0; sig1=0.43; sig2=0.318; sig3=0.36; nc=6; pi=0.0;
k1=0.96; k2=1.355; k3=1.1; k4=1.9; bohr=0.529177; rad=57.29578;

(* C He shell, i=j each *)
Ekin = 2.25*k1/R1^2;
vee=3.0*sig1/R1;
(* vne=-3.0*z/R1; *) vne=0.0;

(* bonding and lone pairs *)
Ekin = Ekin + 2.25*(2*k2/R2^2+k3/R3^2+k4/R4^2);
vee=vee+3.0*(2*sig2/R2+sig3/R3);
oc={-2,-2,-2,-2,-1,-1};
ch={z,1,1,0,0,0};
rr={R1,R2,R2,R3,R4,R4}; w1=w;
cs=Cos[w/2];  ss=Sin[w/2]; cs1=Cos[w1/2]; ss1=Sin[w1/2];
d1=(R1+R2+p)*cs1; d2=(R1+R2+p)*ss1; d3=(R1+R2)*cs; d4=(R1+R2)*ss; d5=R1+R3; d6=R4;
xn={d0,-d1,-d1,0,0,0};
yn={0,d2,-d2,0,0,0};
zn={0,0,0,0,0,0};        
xc={0,-d3,-d3,d5,d0,d0};
yc={0,d4,-d4,0,0,0};
zc={0,0,0,0,d6,-d6};
vne=vne-2*(3-(p/R2)^2)/R2;
vne=vne-z*(3-(d0/R1)^2)/R1;
vnn = 0.0;

(*ww*)
For[i = 1, i < nc, i++,
  For[j = i+1, j < nc+1, j++,
    vee = vee + oc[[i]]*oc[[j]]/Sqrt[(xc[[i]]-xc[[j]])^2+(yc[[i]]-yc[[j]])^2+(zc[[i]]-zc[[j]])^2]]]

(*nn*)
For[i = 1, i < 3, i++,
  For[j = i+1, j < 4, j++,
    vnn = vnn + ch[[i]]*ch[[j]]/Sqrt[(xn[[i]]-xn[[j]])^2+(yn[[i]]-yn[[j]])^2+(zn[[i]]-zn[[j]])^2]]]

(*nw*)
For[i = 1, i < nc+1, i++,
  For[j = 1, j < 4, j++,
    If[i ≠ j,
      vne = vne + oc[[i]]*ch[[j]]/Sqrt[(xc[[i]]-xn[[j]])^2+(yc[[i]]-yn[[j]])^2+(zc[[i]]-zn[[j]])^2]]]]

Epot=vne+vee+vnn;
func=Ekin+Epot;
t = FindMinimum[func,{R1,0.186},{d0,0.001},{R2,1.04},{R3,0.98},{R4,0.93},{w,1.9},{p,0.58},{Method -> "Newton"}, {MaxIterations -> 500}]

(* N[func, 14] *)
vne /. t[[2]]
vee /. t[[2]]
vnn /. t[[2]]
-Epot/Ekin /. t[[2]]
(R1+R2+p)*bohr /. t[[2]]
w*rad /. t[[2]]
w1*rad /. t[[2]]

(* Ovnfx: sum of forces between O1 - nuclei, x-direction *)
Ovnfx=0.0;
For[j=2,j<4,j++,
Ovnfx=Ovnfx+(xn[[j]]-xn[[1]])*ch[[1]]*ch[[j]]/((xn[[1]]-xn[[j]])^2+(yn[[1]]-yn[[j]])^2+
(zn[[1]]-zn[[j]])^2)^(3/2)];

(* Ovefx: sum of forces between O1 - clouds, x-direction *)
Ovefx=0.0;
For[j=2,j<nc+1,j++,
Ovefx=Ovefx+(xc[[j]]-xn[[1]])*oc[[j]]*ch[[1]]/((xc[[j]]-xn[[1]])^2+(yc[[j]]-yn[[1]])^2+
(zc[[j]]-zn[[1]])^2)^(3/2)];
Ovefx=Ovefx+2*z*d0/R1^3*cs/.t[[2]];
Oforcetotx=Ovnfx+Ovefx/.t[[2]]
Oforcennx=Ovnfx/.t[[2]]
OForcenex=Ovefx/.t[[2]]

(* Hvnfx: sum of forces between H2 - nuclei, x-direction *)
Hvnfx = 0.0;
  For[j = 1, j < 2, j++,
    If[j ≠ 2,
    Hvnfx = Hvnfx + (xn[[j]]-xn[[2]])*ch[[2]]*ch[[j]]/((xn[[2]]-xn[[j]])^2+(yn[[2]]-yn[[j]])^2+(zn[[2]]-zn[[j]])^2)^(3/2)]]

(* Hvefx: sum of forces between H2 - clouds, x-direction *)
Hvefx = 0.0;
  For[j = 1, j < nc+1, j++,
    If[j ≠ 2,
    Hvefx = Hvefx + (xc[[j]]-xn[[2]])*oc[[j]]*ch[[2]]/((xc[[j]]-xn[[2]])^2+(yc[[j]]-yn[[2]])^2+(zc[[j]]-zn[[2]])^2)^(3/2)]]
Hvefx = Hvefx - 2*p/R2^3*cs /. t[[2]];
Hforcetotx = Hvnfx+Hvefx/. t[[2]]
Hforcennx=Hvnfx /. t[[2]]
HForcenex=Hvefx /. t[[2]]

(* Hvnfy: sum of forces between H2-nuclei, y-direction *)
Hvnfy = 0.0;
  For[j = 1, j < 4, j++,
    If[j ≠ 2,
    Hvnfy = Hvnfy - (yn[[j]]-yn[[2]])*ch[[j]]*ch[[2]]/((xn[[2]]-xn[[j]])^2+(yn[[2]]-yn[[j]])^2+(zn[[2]]-zn[[j]])^2)^(3/2)]]

(* Hvefy: sum of forces between H2 - clouds, y-direction *)
Hvefy = 0.0;
  For[j = 1, j < nc+1, j++,
    If[j ≠ 2,
    Hvefy = Hvefy - (yc[[j]]-yn[[2]])*oc[[j]]*ch[[2]]/((xc[[j]]-xn[[2]])^2+(yc[[j]]-yn[[2]])^2+(zc[[j]]-zn[[2]])^2)^(3/2)]]
Hvefy = Hvefy - 2*p/R2^3*ss /. t[[2]];
Hforcetoty = Hvnfy+Hvefy /. t[[2]]
Hforcenny=Hvnfy /. t[[2]]
HForceney=Hvefy /. t[[2]]

plot1=Graphics[{Circle[{xc[[1]],yc[[1]]},R1],      Circle[{xc[[2]],yc[[2]]},R2],Circle[{xc[[3]],yc[[3]]},R2],Circle[{xc[[4]],yc[[4]]},R3],Circle[{xc[[5]],yc[[5]]},R4],Disk[{xn[[1]],yn[[1]]},0.08], Disk[{xn[[2]],yn[[2]]},0.08], Disk[{xn[[3]],yn[[3]]},0.08]} ] /. t[[2]];

Show[plot1,{AspectRatio -> Automatic,Axes -> True,GridLines -> Automatic, PlotRange -> {{-3,3},{-3,3}}, Frame -> True}]

plot2=Graphics[{Circle[{xc[[1]],zc[[1]]},R1],      Circle[{xc[[2]],zc[[2]]},R2],Circle[{xc[[3]],zc[[3]]},R2],Circle[{xc[[4]],zc[[4]]},R3],Circle[{xc[[5]],zc[[5]]},R4],Circle[{xc[[6]],zc[[6]]},R4],Disk[{xn[[1]],zn[[1]]},0.08], Disk[{xn[[2]],zn[[2]]},0.08], Disk[{xn[[3]],zn[[3]]},0.08]} ] /. t[[2]];

Show[plot2,{AspectRatio -> Automatic,Axes -> True,GridLines -> Automatic, PlotRange -> {{-3,3},{-3,3}}, Frame -> True}]

plot3=Graphics3D[{{Opacity[0.6],Sphere[{xc[[1]],yc[[1]],zc[[1]]},R1],      Sphere[{xc[[2]],yc[[2]],zc[[2]]},R2],Sphere[{xc[[3]],yc[[3]],zc[[3]]},R2],{Lighter[Green],Sphere[{xc[[4]],yc[[4]],zc[[4]]},R3]},{Lighter[Yellow],Sphere[{xc[[5]],yc[[5]],zc[[5]]},R4],Sphere[{xc[[6]],yc[[6]],zc[[6]]},R4]}},{Darker[Red,1],Sphere[{xn[[1]],yc[[1]],zn[[1]]},0.08],Sphere[{xn[[2]],yc[[2]],zn[[2]]},0.08],Sphere[{xn[[3]],yc[[3]],zn[[3]]},0.08]}} ] /. t[[2]];

Show[plot3,{AspectRatio -> Automatic,Boxed->False}]