sq3=Sqrt[3];
nc=9;
a=Table[{0,0,0},{nc}];
b=a; c=a; g=a; sig=a;
k1=1.0;k2=1.0;k3=1.0;k4=1.0;k5=1.0;
k6=1.0;k7=1.0;
sig1=0.31403448;sig2=0.30821367;sig3=0.30820426;sig4=0.3;sig5=0.3;
sig6=0.3;sig7=0.3;
For[ i=1,i<10,i++, g[[i]]=1.0];
g[[1]]=2.35018989;g[[2]]=0.95875665;g[[6]]=g[[2]];g[[3]]=1.21283096;
g[[4]]=g[[3]];g[[5]]=g[[3]];g[[8]]=g[[3]];g[[9]]=g[[3]];g[[7]]=g[[3]];
sig[[1]]=0.31403448;sig[[2]]=0.30821367;sig[[6]]=sig[[2]];sig[[3]]=0.30820426;
sig[[4]]=sig[[3]];sig[[5]]=sig[[3]];sig[[8]]=sig[[3]];sig[[9]]=sig[[3]];sig[[7]]=sig[[3]];

ss=Sin[x4];cs=Cos[x4];
c[[1,1]]=-2;c[[1,2]]=0;c[[1,3]]=x1;
(*F right*)
c[[2,1]]=-2;c[[2,2]]=9;c[[2,3]]=x2;

c[[3,1]]=-2;c[[3,2]]=0;c[[3,3]]=x3;
c[[4,1]]=-2;c[[4,2]]=0;c[[4,3]]=x3;
c[[5,1]]=-2;c[[5,2]]=0;c[[5,3]]=x3;
(*F left*)
c[[6,1]]=-2;c[[6,2]]=9;c[[6,3]]=x2;

c[[7,1]]=-2;c[[7,2]]=0;c[[7,3]]=x3;
c[[8,1]]=-2;c[[8,2]]=0;c[[8,3]]=x3;
c[[9,1]]=-2;c[[9,2]]=0;c[[9,3]]=x3;

(*coord with zero value not explicitely written!*)
(*nuclei*)
b[[2,1]]=x1+x2;
b[[6,1]]=-b[[2,1]];
(*clouds*)
a[[2,1]]=b[[2,1]];

a[[3,1]]=b[[2,1]]+(x2+x3)*cs;a[[3,2]]=-(x2+x3)*ss;
a[[4,1]]=a[[3,1]];a[[4,2]]=-0.5*a[[3,2]];a[[4,3]]=0.5*sq3*a[[3,2]];
a[[5,1]]=a[[3,1]];a[[5,2]]=a[[4,2]];a[[5,3]]=-a[[4,3]];

a[[6,1]]=-b[[2,1]];
a[[7,1]]=-a[[3,1]];a[[7,2]]=-a[[3,2]];
a[[8,1]]=-a[[3,1]];a[[8,2]]=-a[[4,2]];a[[8,3]]=a[[4,3]];
a[[9,1]]=-a[[3,1]];a[[9,2]]=-a[[4,2]];a[[9,3]]=a[[5,3]];

(*CALC.OF ENERGY-TERMS*)
(*KI/KJ*)
vnn=9*9/(2*(x1+x2));   vee=0.0;vne=0.0;ekin=0.0;
(*WI/WJ EXCL.I=J*)
For[ i=1,i< nc ,i++,
For[ j=i+1,j<nc+1,j++,
d=Sqrt[(a[[i,1]]-a[[j,1]])^2+(a[[i,2]]-a[[j,2]])^2+(a[[i,3]]-a[[j,3]])^2];
If[ d>0 ,vee=vee+c[[i,1]]*c[[j,1]]/d]]]
(*KJ/WI I=J*)
vne=vne+1.5*c[[2,1]]*c[[2,2]]/c[[2,3]];
vne=vne+1.5*c[[6,1]]*c[[6,2]]/c[[6,3]];

(*KJ/WI I<>J inc;*)
For[ i=1,i<nc+1,i++,
If[ i≠2 ,
d=Sqrt[(a[[i,1]]-b[[2,1]])^2+(a[[i,2]])^2+(a[[i,3]])^2];
If[ d>0,vne=vne+c[[i,1]]*c[[2,2]]/d]];
If[ i≠6,
d=Sqrt[(a[[i,1]]-b[[6,1]])^2+(a[[i,2]])^2+(a[[i,3]])^2];
If[ d>0 ,vne=vne+c[[i,1]]*c[[6,2]]/d]]]

For[ i=1,i< nc+1,i++,
(*T*)
ekin=ekin-1.125*g[[i]]*c[[i,1]]/c[[i,3]]^2;
(*WI/WI*)
vee=vee+3.0*sig[[i]]/c[[i,3]]]


Epot=vne+vee+vnn;
func=Ekin+Epot;

t=FindMinimum[func,{x1,1.17843},{x2,0.1595},{x3,0.8044},{x4,1.34},Method->"Newton",MaxIterations->500]