(* H4++ with free start of protons *) 
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];
d14  =  
        Sqrt[(x1 - x4)^2 + (y1 - y4)^2 + (z1 - 
                  z4)^2];
d23  =  
        Sqrt[(x3 - x2)^2 + (y3 - y2)^2 + (z3 - 
                  z2)^2];
d24  =  
        Sqrt[(x4 - x2)^2 + (y4 - y2)^2 + (z4 - 
                  z2)^2];
d34  =  
        Sqrt[(x3 - x4)^2 + (y3 - y4)^2 + (z3 - 
                  z4)^2];

Vnn  =  
        1/d12 + 1/d13 + 1/d14 + 1/d23 + 1/d24 + 
          1/d34;


r1sq  =  x1^2 + y1^2 + z1^2;
r2sq  =  x2^2 + y2^2 + z2^2;
r3sq  =  x3^2 + y3^2 + z3^2;
r4sq  =  
        x4^2 + y4^2 + z4^2;

Vne  =  -(12  -  (r1sq + r2sq + r3sq + r4sq)/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}}, {x4, {0.5, -2}}, {y4, {0.4, -2}}, {z4, {-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
Sqrt[x4^2 + y4^2 + z4^2]  /.  u
  
(x1 + x2 + x3 + x4)/4  /.  u
(y1 + y2 + y3 + y4)/4  /.  u
(z1 + z2 + z3 + z4)/4  /.  u
  
d12  /. u
d13  /.  u
d14  /.  u
d23  /.  u
d24  /.  u
d34 /.u