(* N,N He shell*)
T = 2*2.25*k1/R1^2;
Vee=2*3.0*sig1/R1;
Vne=-2*3.0*z/R1;

a=Table[{0,0,0},{nc}];
b=Table[{0,0,0},{nc}];
c=Table[{0,0,0},{nc}];

(* bonding and lone pairs *)
T = T + 2.25*(2*k2/R2^2+k3/R3^2)+4*1.125*k4/R4^2);
Vee=Vee+3.0*(2*sig2/R2+sig3/R3);

c[[1,1]]=-2;      c[[1,2]]=z;  c[[1,3]]=R1;
c[[2,1]]=-2;      c[[2,2]]=z;  c[[2,3]]=R1;
c[[3,1]]=-2;      c[[3,2]]=0;  c[[3,3]]=R2;
c[[4,1]]=-2;      c[[4,2]]=0;  c[[4,3]]=R2;
c[[5,1]]=-2;      c[[5,2]]=0;  c[[5,3]]=R3;
c[[6,1]]=-1;      c[[6,2]]=0;  c[[6,3]]=R4;
c[[7,1]]=-1;      c[[7,2]]=0;  c[[7,3]]=R4;
c[[8,1]]=-1;      c[[8,2]]=0;  c[[8,3]]=R4;
c[[9,1]]=-1;      c[[9,2]]=0;  c[[9,3]]=R4;


d1=R1+R3;              b[[1,1]]=d1; 
d2=R1+R3+R1+R2;        b[[2,1]]=-d1;

                       a[[1,1]]=d1;
                       a[[2,1]]=-d1;
                       a[[3,1]]=d2;
		       a[[4,1]]=-d2;
                       a[[5,1]]=0.0;             
                       a[[6,1]]=d1;     a[[6,2]]=R4;
                       a[[7,1]]=d1;                 a[[7,3]]=-R4;
                       a[[8,1]]=-d1;    a[[8,2]]=R4;
		       a[[9,1]]=-d1;                 a[[9,3]]=-R4


{ww}
For[ i=1,i<nc,i++,
  For[ j=i+1,j<=nc,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;

{nn}
vnn=z*z/(2*d1);

{ne}
For[ i=1,i<=nc,i++,
  For[ j=i+1,j<=nc,j++,
    If[ i<>j,
      d=Sqrt[(a[[i,1]]-b[[j,1]])^2+(a[[i,2]]-b[[j,2]]^2)+(a[[i,3]]-b[[j,3]])^2];
      If[ d<>0, vne=vne+c[[i,1]]*c[[j,2]]/d;

Ekin=T;
Epot=vne+vee+vnn;
E=Ekin+Epot;