(* H3 +  
        with protons at equilateral triangle trough origin *) 
func  =  
        2*9/(8*R^2) + 6/(5*R) + 
          3*1/(r*Sqrt[3]) - 
3*2*(3/2 - 
                  1/2*(r^2/R^2))/R;
N[t = FindMinimum[func,  {R, 1}, {r, 1}], 
      10]


  
T  =  2.25/R^2;
Vee  =  1.2/R;
d12  =  
        Sqrt[(x1 - x2)^2 + (y1 - y2)^2 + (z1 - 
                  z2)^2];
d13  =  
        Sqrt[(x1 - x3)^2 + (y1 - y3)^2 + (z1 - 
                  z3)^2];
d32  =  
        Sqrt[(x3 - x2)^2 + (y3 - y2)^2 + (z3 - 
                  z2)^2];
Vnn  =  1/d12 + 1/d13 + 1/d32;
r1sq  =  x1^2 + y1^2 + z1^2;
r2sq  =  x2^2 + y2^2 + z2^2;
r3sq  =  
        x3^2 + y3^2 + z3^2;

Vne  =  -(9  -  (r1sq + r2sq + r3sq)/R^2)/
          R;

func  =  T  +  Vee  +  Vnn  +  Vne;
   
N[t  =  
        FindMinimum[
          func,  {R, {1.0, 1.5}}, 
     {x1, {0.4, 
              2}}, {y1, {0.1, 2}}, {z1, {0.2, 
              2}}, {x2, {-0.3, -2}}, {y2, {0.8, 
              2}}, {z2, {-0.3, -2}}, {x3, {-0.5, -2}}, {y3, {-0.4, -2}}, {z3, {0.5, 2}}, {MaxIterations  ->  100}], 
      20]  
u  =  t[[2]];
  
Sqrt[ x1^2 + y1^2 + z1^2]  /.  u
Sqrt[x2^2 + y2^2 + z2^2]  /.  u
Sqrt[x3^2 + y3^2 + z3^2]  /.  u
  
(x1 + x2 + x3)/3  /.  u
(y1 + y2 + y3)/3  /.  u
(z1 + z2 + z3)/
        3 /. u